dom_dec_gen_reduced_data.m

function reduced_data = dom_dec_gen_reduced_data(model, ...
                                                 detailed_data)
% function reduced_data = dom_dec_reduced_data(model,detailed_data)
%
% all required offline data for the rb-simulation is computed.
% there are two ingredients for the error estimator - the jump
% across `\Gamma` and the residuum.
%
% the computation of residuum-components needs
% Q_f+N1*Q_A1+N2*Q_A2 detailed simulations. the number of
% iterations is very high there, because the accuracy of the error
% estimator is good then. therefore very expensive computations.
%
% Generated fields of reduced_data:
%  AN_comp: nested cell containing the components of the system matrizes
%  for the subproblems
%  fN_comp: nested cell containing the components of the right hand
%  side vectors for the subproblems
%  R_11, R_12, R_22: matrizes for the computation of the jump
%  across `\Gamma` component of the error estimator
%  G: cell containing the matrizes for the residuum component
%  N: cell with the dimensions of the RB spaces
%  I_0, I_G: cells with index vectors of `0` and `\Gamma` spaces

% I.Maier, 19.07.2011

dir = model.dirichlet_side;
no_dir = mod(dir,2)+1;

reduced_data = [];

model.decomp_mode = 1;

% A_comp and f_comp are nested cells!
[A_comp, f_comp, dummy1, dummy2] = ...
    model.operators(model,detailed_data);

% number of matrix/vector components
Q_A = cell(1,2);
Q_f = cell(1,2);
N = cell(1,2);
for i = 1:2
  Q_A{i} = length(A_comp{i});
  Q_f{i} = length(f_comp{i});
  
  N{i} = size(detailed_data.RB{i},2);
end;

reduced_data.N = N;
reduced_data.I_0 = detailed_data.I_0;
reduced_data.I_G = detailed_data.I_G;

reduced_data.AN_comp = cell(1,2);
reduced_data.fN_comp = cell(1,2);

for i = 1:2
  reduced_data.AN_comp{i} = cell(1,Q_A{i});
  for q = 1:Q_A{i}
    reduced_data.AN_comp{i}{q} = detailed_data.RB{i}' * ...
      A_comp{i}{q} * detailed_data.RB{i};
  end;

  reduced_data.fN_comp{i} = cell(1,Q_f{i});
  for q = 1:Q_f{i}
    reduced_data.fN_comp{i}{q} = detailed_data.RB{i}' * ...
        f_comp{i}{q};
  end;
end;

% error estimation quantities
K = model.get_inner_product_matrices(detailed_data);

%%%%%%%%%%% computation of extension matrizes

RB_ext = cell(1,2);
% trivial finite element extension
for i = 1:2
  RB_ext{i} = zeros(size(detailed_data.RB{dir},1), ...
                    size(detailed_data.RB{i},2));
  RB_ext{i}(detailed_data.gamma_dofs{dir},:) = ...
      detailed_data.RB{i}(detailed_data.gamma_dofs{i},:);
end;
reduced_data.R_11 = RB_ext{dir}' * K{dir} * RB_ext{dir};
reduced_data.R_12 = RB_ext{dir}' * K{dir} * RB_ext{no_dir};
reduced_data.R_22 = RB_ext{no_dir}' * K{dir} * RB_ext{no_dir};

%%%%%%%%%%%%% computation of residuum-components
% perform a detailed simulation for each component

model.verbose = 0;

res_mat = cell(1,2);
res_rhs = cell(1,2);
for i = 1:2
  res_mat{i} = K{i};
  dir_gids = detailed_data.df_infos{i}.dirichlet_gids;
  mat_dir = sparse(dir_gids,dir_gids,ones(size(dir_gids)), ...
                   detailed_data.df_infos{i}.ndofs, ...
                   detailed_data.df_infos{i}.ndofs);
  res_mat{i}(dir_gids,:) = 0;
  res_mat{i} = res_mat{i} + mat_dir;
end;
clear('mat_dir');

% assume Q_f{1}=Q_f{2}
% v_f has one part on each domain
v_f = cell(1,2);
for i = 1:2
  v_f{i} = zeros(size(res_mat{i},1),Q_f{i});
end;

for q = 1:Q_f{1}
  if any(f_comp{1}{q}) || any(f_comp{2}{q})
      
    % prepare detailed simulation
    detailed_data.operators = @(j) deal(res_mat{j}, f_comp{j}{q});

    % detailed simulation
    sim_data = detailed_simulation(model,detailed_data);
    if sim_data.reached_maxiter
      warning('detailed simulation may have diverged!');
    end;
    for i = 1:2
      v_f{i}(:,q) = sim_data.uh{i}.dofs;
    end;
    
  end;
end;

% two different right hand sides for v_a
% (v_a{1}{m},phi) = a_1(phi_N^m,phi)
% (v_a{2}{m},phi) = a_2(phi_N^m,phi)
v_a = cell(1,2);
for i = 1:2
  
  v_a{i} = cell(1,N{i});
  for m = 1:N{i}
  
    % v_a{i}{m} has one part on each domain
    v_a{i}{m} = cell(1,2);
    for j = 1:2
      v_a{i}{m}{j} = zeros(size(res_mat{j},1),Q_A{i});
    end;
    
    for q = 1:Q_A{i}
      if any(any(A_comp{i}{q}))

        % prepare detailed simulation
        res_rhs{i} = A_comp{i}{q} * detailed_data.RB{i}(:,m);
        res_rhs{mod(i,2)+1} = zeros(size(res_mat{mod(i,2)+1},1),1);
        detailed_data.operators = @(j) deal(res_mat{j}, ...
                                            res_rhs{j});

        % detailed simulation
        sim_data = detailed_simulation(model,detailed_data);
        if sim_data.reached_maxiter
          warning('detailed simulation may have diverged!');
        end;
        for j = 1:2
          v_a{i}{m}{j}(:,q) = sim_data.uh{j}.dofs;
        end;
      
      end;
    end;

  end;
end;

reduced_data.G = cell(1,2);

% two matrizes G (same procedure as in lin_stat_gen_reduced_data)
for i = 1:2

  Gff = v_f{i}' * K{i} * v_f{i};
  
  Gfa1 = cell(1,N{1});
  for n = 1:N{1}
    Gfa1{n} = v_f{i}' * K{i} * v_a{1}{n}{i};
  end;
  Gfa2 = cell(1,N{2});
  for n = 1:N{2}
    Gfa2{n} = v_f{i}' * K{i} * v_a{2}{n}{i};
  end;
  
  Ga1a1 = cell(N{1},N{1});
  for n1 = 1:N{1}
    for n2 = 1:N{1}
      Ga1a1{n1,n2} = v_a{1}{n1}{i}' * K{i} * v_a{1}{n2}{i};
    end;
  end;
  Ga1a2 = cell(N{1},N{2});
  for n1 = 1:N{1}
    for n2 = 1:N{2}
      Ga1a2{n1,n2} = v_a{1}{n1}{i}' * K{i} * v_a{2}{n2}{i};
    end;
  end;
  Ga2a2 = cell(N{2},N{2});
  for n1 = 1:N{2}
    for n2 = 1:N{2}
      Ga2a2{n1,n2} = v_a{2}{n1}{i}' * K{i} * v_a{2}{n2}{i};
    end;
  end;
  
  Q_r = Q_f{1}+N{1}*Q_A{1}+N{2}*Q_A{2};
  G = zeros(Q_r,Q_r);
  G(1:Q_f{1},1:Q_f{1}) = Gff;
  for n = 1:N{1}
    G(1:Q_f{1},Q_f{1}+(n-1)*Q_A{1}+(1:Q_A{1})) = Gfa1{n};
    G(Q_f{1}+(n-1)*Q_A{1}+(1:Q_A{1}),1:Q_f{1}) = Gfa1{n}';
  end;
  for n = 1:N{2}
    G(1:Q_f{1},Q_f{1}+N{1}*Q_A{1}+(n-1)*Q_A{2}+(1:Q_A{2})) = ...
        Gfa2{n};
    G(Q_f{1}+N{1}*Q_A{1}+(n-1)*Q_A{2}+(1:Q_A{2}),1:Q_f{1}) = ...
        Gfa2{n}';
  end;
  for n1 = 1:N{1}
    for n2 = 1:N{1}
      G(Q_f{1}+(n1-1)*Q_A{1}+(1:Q_A{1}),Q_f{1}+(n2-1)*Q_A{1}+(1:Q_A{1})) ...
          = Ga1a1{n1,n2};
    end;
  end;

  for n1 = 1:N{1}
    for n2 = 1:N{2}
      G(Q_f{1}+(n1-1)*Q_A{1}+(1:Q_A{1}),Q_f{1}+N{1}*Q_A{1}+ ...
        (n2-1)*Q_A{2}+(1:Q_A{2})) = Ga1a2{n1,n2};
      G(Q_f{1}+N{1}*Q_A{1}+(n2-1)*Q_A{2}+(1:Q_A{2}),Q_f{1}+ ...
        (n1-1)*Q_A{1}+(1:Q_A{1})) = Ga1a2{n1,n2}';
    end;
  end;

  for n1 = 1:N{2}
    for n2 = 1:N{2}
      G(Q_f{1}+N{1}*Q_A{1}+(n1-1)*Q_A{2}+(1:Q_A{2}),Q_f{1}+ ...
        N{1}*Q_A{1}+(n2-1)*Q_A{2}+(1:Q_A{2})) = Ga2a2{n1,n2};
    end;
  end;
    
  reduced_data.G{i} = G;
  
end;