demo_fem.m

Go to the documentation of this file.

function demo_fem;
% demo of a finite element simulation 
%
% Script solving a laplace problem (-laplace u = 0) with suitable
% general dirichlet and neuman boundary conditions.

% Part 1:
%  direct assembly of Q1-FEM on cartesian grids 
% Part 2: 
%  arbitrary degree P_k triangular elements on triangular grid 
%  via RBmatlab fem routines
% Part 3: Convergence order study of RBmatlab P_k finite elements.
%  via RBmatlab fem routines
%
% See also demo_fem2 for examples using
% the more general RBmatlab +Fem Routines
% (including Stokes/Navier-Stokes examples)

% Bernard Haasdonk 5.7.2006

disp('===========================================')
disp('Part A: bilinear Q1-FEM on cartesian grid')

disp('Demonstration that solves a laplace-problem -laplace u = 0:')
disp('with nontrivial Dirichlet and Neumann values.');

params = [];

disp('setting parameters');

params.xrange = [0,1]; % grid x-coordinate range
params.yrange = [0,1]; % grid y-coordinate range
params.xnumintervals = 100; % number of grid cells in x-direction
params.ynumintervals = 50; % number of grid cells in y-direction

% the following exceeds the size of sparse-matrices, which can be
% linearly indexed (180000 elements, ... sec comp time)
%params.nx = 600;       % number of grid cells in x-direction
%params.ny = 300;       % number of grid cells in y-direction

% the following is trivially solvable by linear indices <= 1.6e9 
%params.nx = 200;       % number of grid cells in x-direction
%params.ny = 100;       % number of grid cells in y-direction

% set upper boundary to dirichlet, middle of upper and the lower boundary
% to neuman, set left and right to neuman
params.bnd_rect_corner1 = [params.xrange(1) , params.yrange(2) ; ...
                    params.xrange(1)+ 250e-6 , params.yrange(2); ...
                    params.xrange(1) , params.yrange(1); ...
                    params.xrange(1) , params.yrange(1); ...
                    params.xrange(2) , params.yrange(1);
                   ];
params.bnd_rect_corner1 = params.bnd_rect_corner1' - eps;
params.bnd_rect_corner2 = [params.xrange(2), params.yrange(2); ...
                    params.xrange(2)-250e-6, params.yrange(2); ...
                    params.xrange(2), params.yrange(1); ...
                    params.xrange(1), params.yrange(2); ...
                    params.xrange(2), params.yrange(2)
                   ];
params.bnd_rect_corner2 = params.bnd_rect_corner2' + eps;

params.bnd_rect_index = [-1, -2, -2, -2, -2];

params.gridtype = 'rectgrid';
grid = construct_grid(params);
plot(grid);
title('FEM grid');

% set data functions
params.name_dirichlet_values = 'leftright';
params.c_dir_left = 1.0; % pressure left
params.c_dir_right = 0.5; % pressure right
params.dir_middle = params.xrange(2)/2;

%params.name_dirichlet_values = 'homogeneous';
%params.c_dir = 1.0;
params.name_neuman_values = 'pressure_gdl';
params.c_neu_max = 1000; % maximum velocity in middle of left and right side

%params.name_neuman_values = 'homogeneous';
%params.c_neu = 0;

params.show_colorbar = 1;
params.nolines = 1;
params.verbose = 10;

disp('solving system');
p = fem_laplace(grid,params);

disp('plotting result');
figure;

plot_vertex_data(grid,p,params);
%plot_p1_data(p,params);
title('FEM solution of -Laplace u = 0 with nontrivial boundary values');

disp('press key to continue')
pause;




disp('===============================================================')
disp('Part B: P_k-FEM on triangular grid and convergence order study')

  % simple FEM simulation and plot
  % initialize model:
  params = [];
  params.pdeg = 2;
  params.numintervals = 2*2^4;
  model = poisson_model(params);
  % generate grid and fem matrices:
  model_data = gen_model_data(model);
  figure, plot(model_data.grid);
  axis equal; axis tight; title('FEM grid')
  % perform full simulation 
  sim_data = detailed_simulation(model, model_data); 
  % plot results and provide slider
  figure, plot_sim_data(model,model_data,sim_data,[]);
  title('FEM solution of -Laplace u = f for zero-boundary values')
  disp(['ndofs = ',num2str(sim_data.uh.df_info.ndofs)]);
 
  disp('===============================')
  disp('error convergence table:')
  
  i_s = [0,1,2,3,4];
  nis = length(i_s);
  l2errs = zeros(nis,1);
  h1errs = zeros(nis,1);
  ndofs = zeros(nis,1);
  
  for pdeg = 1:2
    disp('===================================================================')
    disp(['Convergence table for pdeg = ',num2str(pdeg)]);
    disp(' ');
    disp('  1/numintervals  |   ndofs   |    L2-error      |     H1-error   ')
    disp('-------------------------------------------------------------------')
    for i = 1:length(i_s)
      params.numintervals = 2*2^i_s(i);
      params.pdeg = pdeg; 
      model = poisson_model(params);
      model_data = gen_model_data(model);
      sim_data = detailed_simulation(model,model_data);
      
      % project exact solution onto higher degree polynomial fem func
      par.pdeg = 4;
      par.qdeg = 8;
      par.dimrange = 1;    
      p4_df_info = feminfo(par,model_data.grid);
      
      uexact_h = femdiscfunc([],p4_df_info);
      uexact_h = fem_interpol_global(model.solution,uexact_h);
      uh = femdiscfunc([],p4_df_info);
      u_local_eval = @(grid,elids,lcoord,params) ...
          my_uh_local_eval(grid,elids,lcoord,params,sim_data.uh); 
      uh = fem_interpol_local(u_local_eval,uh);
      
      err = uh - uexact_h;
      l2errs(i) = fem_l2_norm(err);
      h1errs(i) = fem_h1_norm(err);
      ndofs(i) = sim_data.uh.df_info.ndofs;
      disp(['   ',num2str(1/params.numintervals,'%10.5e'),'   |    ',...
            num2str(ndofs(i),'%10.4d'),'   |   ',...
            num2str(l2errs(i),'%10.5e'),'   |   ',...
            num2str(h1errs(i),'%10.5e')]);
    end;
  end;
  
function res = my_uh_local_eval(grid,elids,lcoord,params,df)
% dummy function used for evaluating a discrete function at finer
% lagrange-grid nodes
res = fem_evaluate(df,elids,lcoord,[],[]);