fem2_poisson.m

function res = fem2_poisson(step,params)
%function res = fem2_poisson(step,params)
%
% Script demonstrating lagrange finite element functions, using
% them for function interpolation and finite element methods.
%
% Script realizing the simple poisson equation or a more complex
% elliptic problem with arbitrary diffusivity_tensor, advection,
% source, reaction, neumann , dirichlet and robin boundary
% conditions on the unit square
% with finite elements on triagrid discretization.
%
% step 1: some basic femdiscfunc demonstration and functionatity
% step 2: use femdiscfunc for interpolating functions by local and
%         global evaluations
% step 3: solve fem-problem, return errors
% step 4  error convergence investigation
% step 5: same as 3 by call of interface methods + Reduced basis steps
% step 6: demo_rb_gui

% B. Haasdonk 11.1.2011

% adapted to new assembly

% Immanuel Maier 12.04.2011

% poisson_model

if nargin == 0
  step = 1;
end;

if nargin < 2
  params = [];
end;

res = [];

switch step
  
 case 1 % elementary fem discrete function operations

  params = [];
  pdeg = 4;
  params.pdeg = pdeg;
  params.dimrange = 3;
  % params.dimrange = 1;
  params.debug = 1;
  params.numintervals = 5;
  model = poisson_model(params);
  % convert to local_model:
  model = elliptic_discrete_model(model);
  % convert to generic-fem model:
  model = generic_fem_model_adapter(model);
  grid = construct_grid(model);
  % generate function info (composite)
  df_info = Fem.Lagrange.DefaultInfo(model.pdeg, grid);
  ids = cell(1, params.dimrange);
  for i = 1:params.dimrange
      ids{i} = ['dim', num2str(i)];
  end
  df_info = Fem.CompositeFunctionSpace(ids, repmat({df_info}, 1, params.dimrange));
  %tmp = load('circle_grid');
  %grid = triagrid(tmp.p,tmp.t,[]);
  disp('model initialized');

  % without model_data, detailed_simulation, etc. but explicit
  % calls of fem operations
    
  % initialize vectorial discrete function, extract scalar
  % component and plot basis function
  df = Fem.DiscFunc([],df_info); % initialize zero df
  
  fprintf('\n');
  disp('initialization of femdiscfunc successful, display:');
  display(df);
  
  disp('press enter to continue');
  pause;
  
  % for check of dof consistency we have a plot routine:
  fprintf('\n');
  disp('plot of global_dof_index map for consistency check (scalar):');
  p = plot_dofmap(get(df_info, 1));
  
  disp('press enter to continue');
  pause;
  
  % global dof access:
  fprintf('\n');
  disp('example of dof access, scalar component extraction and plot:');    
  % set second vector component in 4th basis function nonzero;
  ncomp = 2;
  locbasisfunc_index = 4;
  df.dofs((ncomp-1)*df.ndofs/df.dimrange+locbasisfunc_index) = 1; 
  dfscalar = get_component(df,2); % should be nonzero function
  disp(['entry should be 1 : ',...
        num2str(dfscalar.dofs(locbasisfunc_index))]); 
  
  params = [];
  params.subsampling_level = 6;
  figure,plot(dfscalar,params);
  dfscalar.dofs(:) = 0; % now zero again.
%  keyboard;
  
  disp('press enter to continue');
  pause;
  
  fprintf('\n');
  disp('example of local dof access:');
  % local dof access via local to global map
  eind = 4; % set dof on element number 4
  lind = (pdeg+2); % local basis function index with lcoord (0,1/pdeg) 
  dfscalar.dofs(dfscalar.df_info.global_dof_index(eind,lind)) = 1;
  
  % example of local evaluation of femdiscfunc simultaneous on
  % several elements in the same local coordinate point
  elids = [4,6,7,10]; % evaluation on some elements
  lcoord = [0,1/pdeg]; % local coordinate vector == last point in all triangles
  f = evaluate(dfscalar,elids,lcoord); 
  disp(['first entry should be 1 : ',num2str(f(1))]); 
  % equivalent call (!) by () operator as abbreviation for local evaluation:
  f = dfscalar(elids,lcoord); 
  disp(['first entry should be 1 : ',num2str(f(1))]); 
  
  disp('press enter to continue');
  pause;

  disp('examples of norm computation:')
  params.dimrange = 1;
  params.pdeg = 1;
  dfinfo1 = Fem.Lagrange.DefaultInfo(params.pdeg,grid);
  df1 = Fem.DiscFunc([],dfinfo1);
  df1.dofs(:) = 1;
  disp(['L2-norm(f(x,y)=1) = ',num2str(l2_norm(df1))]);
  disp(['H10-norm(f(x,y)=1) = ',num2str(h10_norm(df1))]);
  df1.dofs(:) = df1.grid.X(:);
  disp(['L2-norm(f(x,y)=x) = ',num2str(l2_norm(df1))]);
  disp(['H10-norm(f(x,y)=x) = ',num2str(h10_norm(df1))]);
  df1.dofs(:) = df1.grid.Y(:);
  disp(['L2-norm(f(x,y)=y) = ',num2str(l2_norm(df1))]);
  disp(['H10-norm(f(x,y)=y) = ',num2str(h10_norm(df1))]);
  
  disp('press enter to continue');
  pause;
  
  % evaluate df in all lagrange nodes of element 4 by loop
  fprintf('\n');
  disp(['dfscalar on element 4 in all lagrange nodes,' ...
        'only (pdeg+2) entry should be 1:']);
  lagrange_nodes = Fem.Lagrange.nodes_lcoord(pdeg);
  elid = 4;
  for i = 1:size(lagrange_nodes,1);
    f = evaluate(dfscalar,elid,lagrange_nodes(i,:)); 
    disp(['f(l(',num2str(i),')) = ',num2str(f)]);
  end;

  disp('press enter to continue');
  pause;
  
  fprintf('\n');
  disp('example of requirement of subsampling in plot of discfuncs:');
  figure;
  subsamp_levels = [0,2,4,16];
  for i=1:length(subsamp_levels)
    subplot(2,2,i),
    params.axis_equal = 1;
    params.subsampling_level = subsamp_levels(i);
    params.clim = [-0.15    1.15]; % rough bounds
    plot(dfscalar,params);
    title(['subsampling level = ',num2str(subsamp_levels(i))]);
  end;

  disp('end of step 1, please inspect variables, if required.')
  keyboard;    

 case 2 

  params = [];
  pdeg = 4;
  params.pdeg = pdeg;
  params.dimrange = 1;
  params.numintervals = 5;
  model = poisson_model(params);
  % convert to local_model:
  model = elliptic_discrete_model(model);
  % convert to generic-fem model:
  model = generic_fem_model_adapter(model);
  grid = construct_grid(model);
  %tmp = load('circle_grid');
  %grid = triagrid(tmp.p,tmp.t,[]);
  disp('model initialized');
  
  % interpolate exact solution and other coefficient functions 
  % as fem-function and plot

  disp('examples of interpolation of analytical functions:');
  df_info=Fem.Lagrange.DefaultInfo(model.pdeg,grid);

  df = interpol_global(df_info, model.solution);  
  plot(df);
  title('analytical solution');

  % discretize source function and plot
  df = interpol_local(df_info, model.source);
  figure,plot(df);
  title('source function');

  disp('press enter to continue');
  pause;
  
  % discretize diffusivity and plot as 4-sequence
  % to be implemented for vectorial functions... 
  %params.dimrange = 4;
  %df4 = femdiscfunc([],grid,params);
  %df4 = fem_interpol_global(model.diffusivity_tensor, df)  
  % arrange as sequence of 4 scalar functions and plot
  %plot(df4);
  %title = 'diffusivity function';
  
  % example of local evaluation of vectorial basis function 
  % on reference triangle
  disp('local evaluation of vectorial basis functions on reference triangle:');
  lcoord = [0,1];
  params = [];
  params.dimrange = 3;
  params.pdeg = 2;
  df2 = Fem.DiscFunc([],df_info);
  res = evaluate_basis(df_info,lcoord)

  disp('local evaluation of scalar basis functions derivative on reference triangle:');
  res = evaluate_scalar_basis_derivative(df_info,lcoord)
  
  % later:
  % disp('local evaluation of vectorial basis function derivative on reference triangle:');

  % later: 
   % interpolation error convergence with grid refinement
 
 case 3 % finite element assembly, solution and plot.
  % different components assembled separately, as later
  % affine-decomposition will require such modularization
  
  % assemble discrete system, solve and plot results

  if isempty(params)
    params = [];
    params.dimrange = 1;
    pdeg = 1;
    qdeg = 3 * pdeg; % quadrature_degree
    params.pdeg = pdeg;
    params.qdeg = qdeg;
    params.numintervals = 20;
    params.no_plot = 0; %plotting active
  end;

  no_plot = params.no_plot;
  
  %model = poisson_model(params);
  model = elliptic_debug_model(params);
  model = elliptic_discrete_model(model); % convert to local_model
  % convert to generic-fem model:
  model = generic_fem_model_adapter(model);
  grid = construct_grid(model);
  disp(['grid.nelements=',num2str(grid.nelements)]);
  % for debugging: set all to neumann
  %  i = find(grid.NBI==-1);
  %  grid.NBI(i)=-2;
  
  df_info = Fem.Lagrange.DefaultInfo(model.pdeg,grid);
  % generate information for boundary conditions:
  bc_info = Fem.BcInfo(model, grid, df_info);
  disp(['df_info.ndofs=',num2str(df_info.ndofs)]);
  %tmp = load('circle_grid');
  %grid = triagrid(tmp.p,tmp.t,[]);
  disp('model initialized');
  
  %%%%%%%%%%%%%%% assembly of system %%%%%%%%%%%%%%%%%%%

  % assemble right hand side
  %[r_source, r_dirichlet, r_neumann, r_robin] = ...
  %    fem_rhs_parts_assembly(model,df_info);
  %r = r_source + r_neumann + r_robin + r_dirichlet;
  r_vol = Fem.Assembly.rhs_part(...
      @(varargin) Fem.Assembly.add_integral_kernels(...
          model.rhs_volume_int_kernel, varargin{:}), model, ...
      df_info);
  r_bnd = Fem.Assembly.rhs_part(model.rhs_boundary_int_kernel, ...
                                model, df_info, bc_info);
  r_dir = Fem.Assembly.dirichlet_dof_vector(model, df_info, ...
                                            bc_info);
  r_vol(bc_info.dirichlet_gids) = 0;
  r_bnd(bc_info.dirichlet_gids) = 0;
  r = r_vol + r_bnd + r_dir;
  disp('rhs assembled');
  
  % sparse system matrix:
  A = spalloc(df_info.ndofs,df_info.ndofs,10); % rough upper bound for number of nonzeros
  %[A_diff , A_adv, A_reac, A_dirichlet, A_neumann, A_robin] = ...
  %    fem_matrix_parts_assembly(model,df_info);  
  %A = A_diff + A_adv + A_reac + A_neumann + A_robin + A_dirichlet;
  A_vol = Fem.Assembly.matrix_part(...
      @(varargin) Fem.Assembly.add_integral_kernels(...
          model.matrix_volume_int_kernel, varargin{:}), model, ...
      df_info);
  A_bnd = Fem.Assembly.matrix_part(model.matrix_boundary_int_kernel, ...
                                   model, df_info, bc_info);
  dir_gids = bc_info.dirichlet_gids;
  A_dir= sparse(dir_gids, dir_gids, ones(size(dir_gids)), ...
                df_info.ndofs, df_info.ndofs);
  A_vol(dir_gids, :) = 0;
  A_bnd(dir_gids, :) = 0;
  A = A_vol + A_bnd + A_dir;
  disp('matrix assembled');
    
  %%%%%%%%%%%%%%% solve system %%%%%%%%%%%%%%%%%%%
  % solution variable
  u_h = Fem.DiscFunc([],df_info);
  u_h.dofs = A\r;    
  disp('system solved');
    
  %%%%%%%%%%%%%%% postprocessing %%%%%%%%%%%%%%%%%%%

  if ~no_plot
    disp('plotting...');
    % plot result
    params.subsampling_level = 10;
    figure, plot(u_h,params);
    title('discrete solution')
  end;   
  
  df = interpol_global(df_info, model.solution, model);  
  if ~no_plot
    figure, plot(df);
    title('analytical solution');
  end;
  
  % plot error of (interpolated) exact solution and discrete sol.
  % interpolation using same degree as discrete solution
  e_h = interpol_global(df_info, model.solution, model);  
  e_h.dofs = e_h.dofs -u_h.dofs;
  if ~no_plot
    figure, plot(e_h);
    title('rough error I\_rough(u)-u_h');
  end;
  
  if ~no_plot
    % better: choose high degree for interpolation of u
    % then require local_interpolation of u_h into high-resolved space
    params.pdeg = 4; % choose high degree
                     
    % map u_h to higher degree discretefunction
    %params.has_dirichlet_values = 10;
    df_info_fine = Fem.Lagrange.DefaultInfo(params.pdeg,grid);
    u_h_loc = @(grid,elids,lcoord,params) ...
              u_h(elids,lcoord,params); % no grid-input for femdiscfunc
    u_h_fine = interpol_local(df_info_fine, u_h_loc, model);  
    e_h_fine = interpol_global(df_info_fine, model.solution, model);  
    e_h_fine.dofs = e_h_fine.dofs -u_h_fine.dofs;
    figure, plot(e_h_fine);
    title('fine error I\_fine(u)-u_h');
  end;
  
  % determine error in L2 and H10 norm
  l2_err = l2_norm(e_h);
  h10_err = h10_norm(e_h);
  disp(['L2-error norm |u_h - I(u_exact)| = ',num2str(l2_err)]);
  disp(['H10-error norm |u_h - I(u_exact)| = ',num2str(h10_err)]);
  
  if ~no_plot
    disp('solution and errors plotted. Inspect workspace.')
    keyboard;
  end;
  
  res.l2_error = l2_err;
  res.h10_error = h10_err;
  
 case 4  % error convergence investigation
  
  params = [];
  params.dimrange = 1;
  pdeg = 2;
  qdeg = 3 * pdeg; % quadrature_degree
  params.pdeg = pdeg;
  params.qdeg = qdeg;
  params.no_plot = 1;
  
  numintervals = [2,4,8,16,32,64,128,256];
  l2_errs = zeros(length(numintervals),1);
  h10_errs = zeros(length(numintervals),1);
  for i = 1:length(numintervals)
    disp(['----------------------------------------------']);
    disp(['numintervals = ',num2str(numintervals(i))]);
    params.numintervals = numintervals(i);
    res = fem2_poisson(3,params);
    l2_errs(i) = res.l2_error;
    h10_errs(i) = res.h10_error;
  end
     
 case 5 % detailed and rb simulation by high-level methods
  
  disp('detailed simulation:');
  model = elliptic_debug_model;
  model = elliptic_discrete_model(model);
  model = generic_fem_model_adapter(model);
  model_data = gen_model_data(model);
  model = model.set_mu(model,[0,0,0]);
  sim_data = detailed_simulation(model,model_data);
  plot_params.title = 'detailed solution';
  figure, plot_sim_data(model,model_data,sim_data,plot_params);
 
  disp('offline computations (gen_detailed_data: basis):');
  detailed_data = gen_detailed_data(model,model_data);
  disp('offline computations II (gen_reduced_data: operator components):');
  reduced_data = gen_reduced_data(model,detailed_data);

  disp('online reduced simulation:');
  model = model.set_mu(model,[0,0,0]);
  rb_sim_data = rb_simulation(model,reduced_data);
  rb_sim_data = rb_reconstruction(model, detailed_data, rb_sim_data);
  plot_params.title = 'reduced solution';
  figure, plot_sim_data(model,model_data,rb_sim_data,plot_params);
  
  eh = sim_data.uh - rb_sim_data.uh;
  l2_err = l2_norm(eh);  
  h10_err = h10_norm(eh);  
  disp(['L2-error: ',num2str(l2_err)]);
  disp(['H10-error: ',num2str(h10_err)]);

 case 6 
  
  disp('demo_rb_gui:')  
  model = elliptic_debug_model;
  model = elliptic_discrete_model(model);
  model = generic_fem_model_adapter(model);
  model_data = gen_model_data(model);
  detailed_data = gen_detailed_data(model,model_data);
  plot_params.axis_tight = 1;
  demo_rb_gui(model,detailed_data,[],plot_params);
    
end

end