stokes_enrich_reduced_data.m

function reduced_data = stokes_enrich_reduced_data(model, ...
                                                   detailed_data)
%function reduced_data = stokes_enrich_reduced_data(model, detailed_data);
%
% only new residual components are computed / orthonormalized

model.decomp_mode = 1;

RB = detailed_data.RB;
N = size(RB, 2);
N_vel = size(get(RB, 1), 2);

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

% former residual basis available?
if ~isfield(detailed_data, 'SRDW')
    % compute only system components
    % (no enrichment)    
    
    Q_A = length(op.A_comp);
    reduced_data.Q_A = Q_A;
    Q_f = length(op.f_comp);
    reduced_data.Q_f = Q_f;
    
    % gram matrix
    K = model.get_inner_product_matrix(detailed_data);
    ndofs = detailed_data.df_info.ndofs;
    K = getMatrix(K) + eps*speye(ndofs);
    reduced_data.KN = RB' * K * 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

    reduced_data.N = N;
    
    % right-hand side
    reduced_data.fN_comp = cell(1, Q_f);
    for q = 1:Q_f
        reduced_data.fN_comp{q} = RB' * op.f_comp{q};
    end
    
    % system matrix
    reduced_data.AN_comp = cell(1, Q_A);
    for q = 1:Q_A
        reduced_data.AN_comp{q} = RB' * (op.A_comp{q} * RB);
    end
    
    % nonlinear parts
    N_nl = 0;
    if N > 0 && model.has_nonlinearity
        
        f_dir = ...
            cell2mat(detailed_data.bc_info.dirichlet_dof_vector_components);
        Q_f_dir = size(f_dir, 2);
        
        model.uh = Fem.DiscFunc([], detailed_data.df_info);
        
        N_nl = size(get(RB, 1), 2);
        Q_A_nl = cell(1, N_nl);
        AN_nl_comp = cell(1, N_nl);
        AN_comp_add = cell(1, N_nl);
        
        pool = gcp('nocreate');
        nWorkers = 0;
        if ~isempty(pool);
            nWorkers = pool.NumWorkers;
        end
        
        % evaluate nonlinearity for each basis function
        parfor (n = 1:N_nl, nWorkers)
            
            model_tmp = model;
            model_tmp.uh.dofs = RB(:, n);
            A_nl_comp = model.operators(model_tmp, detailed_data);
            Q_A_nl{n} = length(A_nl_comp);
            
            AN_comp_add{n} = zeros(N, Q_f_dir*Q_A_nl{n});
            for q = 1:Q_A_nl{n}
                AN_comp_add{n}(:, q:Q_A_nl{n}:Q_f_dir*Q_A_nl{n}) = ...
                    RB' * (A_nl_comp{q} * f_dir);
            end
            
            AN_nl_comp{n} = cell(1, Q_A_nl{n});
            for q = 1:Q_A_nl{n}
                AN_nl_comp{n}{q} = ...
                    RB(:, 1:N_nl)' * (A_nl_comp{q} * RB(:, 1:N_nl));
            end
        end
        
        reduced_data.Q_A_nl = Q_A_nl{1};
        reduced_data.AN_nl_comp = AN_nl_comp;
        clear('AN_nl_comp');
        
        for q = 1:Q_A_nl{1}
            for n = 1:N_nl
                reduced_data.AN_comp{q}(:, n) = ...
                    reduced_data.AN_comp{q}(:, n) + AN_comp_add{n}(:, ...
                                                                  q);
            end
        end
        clear('AN_comp_add');
    end
    reduced_data.N_nl = N_nl;
    
else
    % compute also residual components
    % advanced enrichment process

    reduced_data = detailed_data.SRDW;
    init = isempty(reduced_data);
    reduced_data.G = [];
    
    Q_A = length(op.A_comp);
    reduced_data.Q_A = Q_A;
    Q_f = length(op.f_comp);
    reduced_data.Q_f = Q_f;
    
    % gram matrix
    K = model.get_inner_product_matrix(detailed_data);
    ndofs = detailed_data.df_info.ndofs;
    K = getMatrix(K) + eps*speye(ndofs);
    reduced_data.KN = RB' * K * 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
    
    % determine basis offsets of new basis vectors
    N_end = [N_vel, N];
    if init
        N_start = [1, N_vel + 1];
    else
        N_start = [1 + reduced_data.N_vel, ...
                   1 + reduced_data.N + N_vel - reduced_data.N_vel];
    end
    reduced_data.N = N;
    reduced_data.N_vel = N_vel;
    
    % system matrix
    % residual components
    K_v_a = cell(1, 4*Q_A);
    for q = 1:Q_A
        for j = 1:2
            K_v_a{(q-1)*4+2*j-1} = op.A_comp{q} ...
                * RB(:, (j-1)*N_end(1)+1:N_start(j)-1);
            K_v_a{(q-1)*4+2*j} = op.A_comp{q} ...
                * RB(:, N_start(j):N_end(j));
        end    
    end
    
    % nonlinear parts
    K_v_nl = {};
    Q_A_nl = 0;
    reduced_data.N_nl = N_vel;
    
    if N > 0 && model.has_nonlinearity
        
        f_dir = detailed_data.bc_info.dirichlet_dof_vector_components;
        Q_f_dir = numel(f_dir);
        
        RB_nl = getMatrix(RB);
        RB_nl = RB_nl(:, 1:N_vel);
        model.uh = Fem.DiscFunc([], detailed_data.df_info);
        
        for n = 1:N_vel
            
            model.uh.dofs = RB_nl(:, n);
            A_nl_comp = model.operators(model, detailed_data);
            Q_A_nl = numel(A_nl_comp);
            
            % preallocate memory
            if numel(K_v_nl) == 0
                K_v_nl = cell(1, 3*Q_A_nl);
                for q = 1:3:3*Q_A_nl-2
                    K_v_nl{q} = zeros(ndofs, N_start(1)-1, N_start(1)- ...
                                      1);
                    K_v_nl{q+1} = zeros(ndofs, N_vel-N_start(1)+1, ...
                                        N_start(1)-1);
                    K_v_nl{q+2} = zeros(ndofs, N_vel, N_vel-N_start(1)+ ...
                                        1);
                end
            end
            
            % quadratic components
            for q = 1:3:3*Q_A_nl-2
                if n<N_start(1)
                    K_v_nl{q}(:, :, n) = A_nl_comp{(q+2)/3} * RB_nl(:, ...
                                                                    1:N_start(1)-1);
                    K_v_nl{q+1}(:, :, n) = A_nl_comp{(q+2)/3} * RB_nl(:, ...
                                                                      N_start(1):N_vel);
                else
                    K_v_nl{q+2}(:, :, n-N_start(1)+1) = A_nl_comp{(q+2)/3} ...
                        * RB_nl;
                end
            end
            
            % add to linear terms
            for qA = 1:Q_A_nl
                for qf = 1:Q_f_dir
                    q = (qf-1)*Q_A_nl+qA;
                    j = 1+(n>=N_start(1));
                    nj = n - (n>=N_start(1))*(N_start(1)-1);
                    K_v_a{(q-1)*4+j}(:, nj) = K_v_a{(q-1)*4+j}(:, nj) ...
                        + A_nl_comp{qA} * f_dir{qf};
                end            
            end
        end
        
        reduced_data.Q_A_nl = Q_A_nl;
    end
    
    % rhs
    % system components
    if init
        reduced_data.fN_comp = cell(1, Q_f);
    end
    for q = 1:Q_f
        reduced_data.fN_comp{q} = RB' * op.f_comp{q};
    end
    
    % system matrix
    % system components
    if init
        reduced_data.AN_comp = cell(1, Q_A);
    end
    for q = 1:Q_A
        reduced_data.AN_comp{q} = RB' * horzcat(K_v_a{4*(q-1)+1:4*q});
    end
    
    % nonlinear part
    if N > 0 && model.has_nonlinearity
        if init
            reduced_data.AN_nl_comp = cell(1, N_vel);
            for n = 1:N_vel
                reduced_data.AN_nl_comp{n} = cell(1, Q_A_nl);
            end
        end
        for q = 1:Q_A_nl
            for n = 1:N_vel
                if n<N_start(1)
                    reduced_data.AN_nl_comp{n}{q} = RB_nl' ...
                        * [K_v_nl{3*(q-1)+1}(:, :, n), ...
                           K_v_nl{3*(q-1)+2}(:, :, n)];
                else
                    reduced_data.AN_nl_comp{n}{q} = RB_nl' ...
                        * K_v_nl{3*q}(:, :, n-N_start(1)+1);
                end
            end
        end
    end
    
    K_v_lin = [op.f_comp, K_v_a];
    
    % extract new residual components
    constant_component_indices = [];
    if init
        reduced_data.v = [];
        constant_component_indices = 1:Q_f;
    end
    K_v_new = [K_v_lin{[constant_component_indices, Q_f+2:2:end]}, ...
               reshape(cat(2, K_v_nl{2:3:end}), ndofs, []), ...
               reshape(cat(2, K_v_nl{3:3:end}), ndofs, [])];
    
    if model.has_nonlinearity
        pool = gcp;
        reduced_data.v = VecMat.gram_schmidt_parallel(...
            [reduced_data.v, K \ K_v_new], K, 1000*eps, ...
            size(reduced_data.v, 2));
        delete(pool);
    else
        reduced_data.v = VecMat.gram_schmidt_reiterate(...
            [reduced_data.v, K \ K_v_new], K, 1000*eps, ...
            size(reduced_data.v, 2));
    end
    
    clear('K_v_new');
    
    % compute inner product matrix
    G = cell(1, 1+Q_A_nl);
    G{1} = reduced_data.v' * horzcat(K_v_lin{:});
    
    clear('K_v_lin', 'K_v_a');
    
    for q = 1:Q_A_nl
        G{q+1} = reduced_data.v' * ...
                 reshape(cat(3, cat(2, K_v_nl{3*(q-1)+[1, 2]}), K_v_nl{3*q}), ...
                         ndofs, []);
    end
    
    clear('K_v_nl');
    
    reduced_data.G = horzcat(G{:});
    reduced_data.Q_v = size(reduced_data.G, 2);
    
    clear('G');
    
end
    
% constants needed in error estimator
if isfield(detailed_data, 'infsup_constant')
    reduced_data.infsup_constant = detailed_data.infsup_constant;
else
    warning('rbmatlab:stokes:enrich_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:enrich_reduced_data', ...
                ['no sobolev embedding constant available, default ' ...
                 'value: 1/(2*sqrt(2))']);
        reduced_data.rho_sqr = 0.5/sqrt(2);
    end
end

end