stokes_gen_reduced_data.m

function reduced_data = stokes_gen_reduced_data(model, detailed_data)
%function reduced_data = stokes_gen_reduced_data(model, detailed_data)
%
% Computes reduced data for stokes problems


reduced_data = [];

model.decomp_mode = 1;

pool = gcp('nocreate');
nWorkers = 0;
if ~isempty(pool)
    nWorkers = pool.NumWorkers;
end

N = model.get_rb_size(model, detailed_data);
reduced_data.N = N;

% get matrix and rhs components
op = detailed_data.operators;
if ~op.componentsAssembled
    op.assemble_components(model, detailed_data);
end

Q_A = length(op.A_comp);
reduced_data.Q_A = Q_A;
Q_f = length(op.f_comp);
reduced_data.Q_f = Q_f;

% compute components for residual error estimator
K = model.get_inner_product_matrix(detailed_data);
ndofs = detailed_data.df_info.ndofs;
K = K + eps*speye(ndofs);
reduced_data.KN = detailed_data.RB' * K * detailed_data.RB;

dirichlet_gids = detailed_data.bc_info.dirichlet_gids;

if ~isempty(dirichlet_gids)
    K(dirichlet_gids, :) = 0;
    K(:, dirichlet_gids) = 0;
    K_dirichlet = sparse(dirichlet_gids, dirichlet_gids, ...
        ones(size(dirichlet_gids)), ndofs, ndofs);
    K = K + K_dirichlet;
end
% system matrix complete.

K_v_f = cell2mat(op.f_comp);

% compute reduced rhs components
reduced_data.fN_comp = cell(1, Q_f);
for q = 1:Q_f
    reduced_data.fN_comp{q} = detailed_data.RB' * op.f_comp{q};
end

K_v_a = {zeros(size(K, 1), 0)};
for i = 1:N
    K_v_a{i} = zeros(size(K, 1), Q_A);
    for q = 1:Q_A
        K_v_a{i}(:, q) = op.A_comp{q} * detailed_data.RB(:, i);
    end
end

K_v_a_nl = {zeros(size(K, 1), 0)};
reduced_data.N_nl = 0;
Q_A_nl = 0;
if N > 0 && model.has_nonlinearity
    
    N_nl = size(get(detailed_data.RB, 1), 2);
    f_dir = detailed_data.bc_info.dirichlet_dof_vector_components;
    
    % compute components for each basis function
    model.uh = Fem.DiscFunc([], detailed_data.df_info);
    RB = detailed_data.RB;
    Q_A_nl = cell(1, N_nl);
    K_v_a_nl = cell(1, N_nl);
    
    parfor (i = 1:N_nl, nWorkers)
        
        tmodel = model;
        tmodel.uh.dofs = RB(:, i);
        A_nl_comp = model.operators(tmodel, detailed_data);
        Q_A_nl{i} = length(A_nl_comp);
        
        K_v_a_nl{i} = zeros(size(K, 1), N_nl*Q_A_nl{i});
        for q = 1:Q_A_nl{i}
            K_v_a_nl{i}(:, q:Q_A_nl{i}:N_nl*Q_A_nl{i}) = A_nl_comp{q} * RB(:, 1:N_nl);
        end
        
        for qA = 1:Q_A_nl{i}
            for qf = 1:length(f_dir)
                K_v_a{i}(:, (qf-1)*Q_A_nl{i}+qA) = K_v_a{i}(:, (qf-1)*Q_A_nl{i}+qA) + ...
                    A_nl_comp{qA} * f_dir{qf};
            end
        end
    end
    
    Q_A_nl = Q_A_nl{1};
    reduced_data.AN_nl_comp = cell(1, N_nl);
    for i = 1:N_nl
        reduced_data.AN_nl_comp{i} = cell(1, Q_A_nl);
        for q = 1:Q_A_nl
            reduced_data.AN_nl_comp{i}{q} = ...
                RB(:, 1:N_nl)' * K_v_a_nl{i}(:, q:Q_A_nl:N_nl*Q_A_nl);
        end
    end
    
    reduced_data.Q_A_nl = Q_A_nl;
    reduced_data.N_nl = N_nl;
    
end

K_v_a = cell2mat(K_v_a);
K_v_a_nl = cell2mat(K_v_a_nl);

reduced_data.AN_comp = cell(1, Q_A);
for q = 1:Q_A
    reduced_data.AN_comp{q} = detailed_data.RB' * K_v_a(:, q:Q_A:N*Q_A);
end
    
if isfield(model, 'RB_error_indicator') && ...
        strcmp(model.RB_error_indicator, 'estimator')
    
    K_v_tot = [K_v_f, K_v_a];
    v_tot = K \ K_v_tot;
    if model.has_nonlinearity
        n_nl = detailed_data.df_info.values{1}.ndofs;
        K1 = K(1:n_nl, 1:n_nl);
        v_tot_nl = K1 \ K_v_a_nl(1:n_nl, :);
    end
    
    if model.has_nonlinearity && nWorkers > 0
        v = VecMat.gram_schmidt_parallel(v_tot_nl, K1, 1000*eps);
        v = VecMat.gram_schmidt_parallel(...
            [[v; zeros(size(v_tot, 1)-n_nl, size(v, 2))], v_tot], ...
             K, 1000*eps, size(v, 2));
    else
        [v, ind] = VecMat.gram_schmidt_reiterate(v_tot, K, 1000*eps);
    end
    %v = v_tot;
    
    reduced_data.G = v' * [K_v_tot, K_v_a_nl];
    reduced_data.Q_v = size(reduced_data.G, 2);
    reduced_data.res_onb_ind = ind;
    
    % constants needed in error estimator
    if isfield(detailed_data, 'infsup_constant')
        reduced_data.infsup_constant = detailed_data.infsup_constant;
    else
        warning('rbmatlab:stokes:gen_reduced_data', ...
            'no infsup approximation available, default value: 1');
        reduced_data.infsup_constant = @(mu) 1;
    end
    
    if model.has_nonlinearity
        if isfield(detailed_data, 'rho_sqr')
            reduced_data.rho_sqr = detailed_data.rho_sqr;
        else
            warning('rbmatlab:stokes:gen_reduced_data', ...
                'no sobolev embedding constant available, default value: 1/(2*sqrt(2))');
            reduced_data.rho_sqr = 0.5/sqrt(2);
        end
    end
    
end

end