dom_dec_RB_extension_eigenbasis.m

function RBext_G = dom_dec_RB_extension_eigenbasis(model, ...
                                                  detailed_data)
%function RBext_G = dom_dec_RB_extension_eigenbasis(model, ...
%                                                 detailed_data)
%

% IM, 17.04.2012

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

[RB_0,RB_G] = model.get_rb_from_detailed_data(detailed_data);
clear('RB_0');

RBext_G = cell(1,2);

A = detailed_data.harm_ext_matrices;
G = detailed_data.gamma_inner_product_matrices;
g_dofs = detailed_data.gamma_dofs;

for i = 1:2
  G{i} = G{i}(g_dofs{no_dir},g_dofs{no_dir});
end;

if (size(RB_G{no_dir},2) > 0)
  % compute orthonormalized interface bases
  warning off;
  RB_GG = orthonormalize(RB_G{no_dir}(g_dofs{no_dir},:),G{1},[],'qr');
  warning on;
  clear('RB_G');
  
  % compute orthogonal complement
  RB_GG_compl = eye(length(g_dofs{no_dir})) - ...
      RB_GG * (RB_GG' * G{1});
  clear('RB_GG');
  
  % eliminate linear dependent modes
  RB_GG_compl = orthonormalize(RB_GG_compl,G{1},[],'qr');
  
  % solve generalized eigenvalue problem
  [ep,vp] = eig(RB_GG_compl' * G{2} * RB_GG_compl,...
                eye(size(RB_GG_compl,2)));
  
  evalues = diag(vp);
  
  [evalues,indices] = sort(evalues,1,'ascend');
  RBext_GG = RB_GG_compl * ep(:,indices);
    
else
  % ............. see below
  [ep,vp] = eig(G{2},G{1});
  
  evalues = diag(vp);
  
  % important modes correspond to small eigenvalues ???????
  [evalues,indices] = sort(evalues,1,'ascend');
  RBext_GG = ep(:,indices);
  
end;

% harmonic extensions (XNG1, XNG2)
for i = 1:2
  zeroExt = zeros(detailed_data.df_infos{i}.ndofs, ...
                  size(RBext_GG,2));
  zeroExt(g_dofs{i},:) = RBext_GG;
  RBext_G{i} = A{i}\zeroExt;
end;