porsche_coercivity_alpha.m

function coercivity_alpha = porsche_coercivity_alpha(model)
%function coercivity_alpha = porsche_coercivity_alpha(model)
%
%
% function calculating a lower bound \alpha_{LB} for the coercivity
% constant \alpha
% The Lower bound is
% - according to "Evaluation of Flux Integral Outputs 
%   for the Reduced Basis Method" by Jens L. Eftang and Einar M. Rønquist,
%   2009 (equation (73)):
%   \alpha_{LB}(\mu) = \frac{\lambda^- (\mu)}{\lambda^+(\mu)}
% - according to "A Reduced Basis Element Method for the Steady Stokes 
%   Problem" by Lovgren, Maday, Ronquist, 2006 (equation (3.48)):
%   \alpha_{LB}(\mu) = \lambda_{min}
%
% Implemented here is: \alpha_{LB}(\mu) = \frac{\lambda^- (\mu)}{\lambda^+(\mu)}
%
% Oliver Zeeb, 12.05.11

%Calculate the diffusivity matrices (K) for every macro-transformation,
%then calculate the eigenvalues of theses matrices. See also
%porsche_diffusivity_tensor_coefficients, whic does similar things, but
%without calculation the eigenvalues


   % get the affine transformatins for every makrotriangle:
   % e.g. the transformation from reference to original domain.
   % macro points of the original domain:
   pmacro_ref=model.pmacro;     % macro points of the reference domain
   pmacro_glob = model.pmacro;  % macro points of the original/global domain
   % first 7 points in pmacro_glob are the car points, all the others are
   % points on the rectangular, so the first points have to be
   % changed according to the values in model.mu_names
   
   % set the macro_points according to the values in model.mu_names
   for k=1:length(model.mu_names)
       position_in_pmacro = str2double(model.mu_names{k}(2)); 
       x_or_y = model.mu_names{k}(1);
       if x_or_y == 'x'
           x_or_y = 1;
       elseif x_or_y == 'y'
           x_or_y = 2;
       else
           error('unknown entry in model.mu_names');
       end
       
       mu_k = getfield(model, model.mu_names{k});
       
       pmacro_glob(position_in_pmacro,x_or_y) = mu_k;
   end
   

   % get all the transformations of the macrotriangles to the
   % reference-triangle, calculate the K-Matrices and add them to the
   % coefficient-vector
   K_o = eye(2);
   for k=1:length(model.tmacro)       
       tria_pts_ref  = pmacro_ref(model.tmacro(k,:),:);
       tria_pts_glob = pmacro_glob(model.tmacro(k,:),:);
       [C, G] = aff_trafo_coef(tria_pts_ref, tria_pts_glob);
       det_G=det(G);
       G_inv = inv(G);
       
       %define the K-matrices (see equation (6.22) )
       % K_o,k = eyes for potential flow!
       K=det_G*G_inv*K_o*G_inv';
       eigenvalue_min_K(k)=min(eig(K));
       eigenvalue_max_K(k)=max(eig(K));
   end
   coercivity_alpha=min(eigenvalue_min_K)/max(eigenvalue_max_K);