demo_fem.m

Go to the documentation of this file.

% demo of a finite element simulation 
%
% Script solving a laplace problem (-laplace u = 0) with suitable
% general dirichlet and neuman boundary conditions
% Currently only cartesian grids are supported in the underlying
% fem routine

% Bernard Haasdonk 5.7.2006

%| \todo make work again for release

disp('Demonstration that solves a laplace-problem -laplace u = 0:')

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);

% 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);

  




% TO BE ADJUSTED TO NEW SYNTAX
%| \docupdate