Description

EX_POISSON3 2D Poisson equation example on a unit square.

[ FEA, OUT ] = EX_POISSON3( VARARGIN ) Poisson equation on a [0..1]^2 unit square.

Input       Value/{Default}        Description
-----------------------------------------------------------------------------------
igrid       scalar 1/{0}           Cell type (0=quadrilaterals, 1=triangles)
hmax        scalar {0.1}           Max grid cell size
sfun        string {sflag1}        Shape function
iphys       scalar 0/{1}           Use physics mode to define problem    (=1)
                                   or directly define fea.eqn/bdr fields (=0)
iplot       scalar 0/{1}           Plot solution (=1)
                                                                                  .
Output      Value/(Size)           Description
-----------------------------------------------------------------------------------
fea         struct                 Problem definition struct
out         struct                 Output struct

Code listing

 cOptDef = { ...
   'igrid',    0; ...
   'hmax',     0.1; ...
   'sfun',     'sflag1'; ...
   'iphys',    1; ...
   'icub',     2; ...
   'iplot',    1; ...
   'tol',      0.01;
   'fid',      1 };
 [got,opt] = parseopt(cOptDef,varargin{:});
 fid       = opt.fid;


% Geometry definition.
 fea.geom.objects = { gobj_rectangle() };


% Grid generation.
 switch opt.igrid
   case -2
     n  = round(1/opt.hmax);
     fea.grid = rectgrid(n,n,[0 1;0 1]);
     ix = setdiff( 1:(n+1)^2, [1:(n+1) (n+1):(n+1):(n+1)^2 (n+1)^2:-1:(n+1)^2-(n+1)+1 (n+1)^2-(n+1)+1:-(n+1):1 ] );
     h  = 0.3*1/n;
     fea.grid.p(1,ix) = fea.grid.p(1,ix) + h*(2*rand(1,numel(ix)) - 1);
     fea.grid.p(2,ix) = fea.grid.p(2,ix) + h*(2*rand(1,numel(ix)) - 1);
   case -1
     fea.grid = rectgrid(round(1/opt.hmax),round(1/opt.hmax),[0 1;0 1]);
     fea.grid = quad2tri(fea.grid);
   case 0
     fea.grid = rectgrid(round(1/opt.hmax),round(1/opt.hmax),[0 1;0 1]);
   case 1
     fea.grid = gridgen(fea,'hmax',opt.hmax,'fid',fid);
 end
 n_bdr = max(fea.grid.b(3,:));           % Number of boundaries.


% Problem definition.
 fea.sdim  = { 'x' 'y' };                % Coordinate names.
 if ( opt.iphys==1 )

   fea = addphys(fea,@poisson);          % Add Poisson equation physics mode.
   fea.phys.poi.sfun = { opt.sfun };     % Set shape function.
   fea.phys.poi.eqn.coef{3,4} = { 1 };   % Set source term coefficient.
   fea.phys.poi.bdr.coef{1,end} = repmat({0},1,n_bdr);   % Set Dirichlet boundary coefficient to zero.
   fea = parsephys(fea);                 % Check and parse physics modes.

   if( any(strcmp(opt.sfun,{'sf_tri_H3','sf_quad_H3'})) )   % Prescribed derivatives at end points for Hermite elements.
     fea.bdr.d = {{ 0 0 0 0 ;
                    0 '-ex_poisson3_derivative(y)' 0 'ex_poisson3_derivative(y)' ;
                    'ex_poisson3_derivative(x)' 0 '-ex_poisson3_derivative(x)' 0 }};
   end

 else

   fea.dvar  = { 'u' };                  % Dependent variable name.
   fea.sfun  = { opt.sfun };             % Shape function.

% Define equation system.
   fea.eqn.a.form = { [2 3;2 3] };       % First row indicates test function space   (2=x-derivative + 3=y-derivative),
% second row indicates trial function space (2=x-derivative + 3=y-derivative).
   fea.eqn.a.coef = { 1 };               % Coefficient used in assembling stiffness matrix.

   fea.eqn.f.form = { 1 };               % Test function space to evaluate in right hand side (1=function values).
   fea.eqn.f.coef = { 1 };               % Coefficient used in right hand side.

% Define boundary conditions.
   if( any(strcmp(opt.sfun,{'sf_tri_H3','sf_quad_H3'})) )   % Prescribed derivatives at end points for Hermite elements.
     fea.bdr.d = {{ 0 0 0 0 ;
                    0 '-ex_poisson3_derivative(y)' 0 'ex_poisson3_derivative(y)' ;
                    'ex_poisson3_derivative(x)' 0 '-ex_poisson3_derivative(x)' 0 }};
   else
     fea.bdr.d     = cell(1,n_bdr);
    [fea.bdr.d{:}] = deal(0);          % Assign zero to all boundaries (Dirichlet).
   end
   fea.bdr.n     = cell(1,n_bdr);    % No Neumann boundaries ('fea.bdr.n' empty).

 end


% Parse and solve problem.
 fea       = parseprob(fea);             % Check and parse problem struct.
 fea.sol.u = solvestat(fea,'fid',fid,'icub',opt.icub);   % Call to stationary solver.


% Postprocessing.
 if( opt.iplot>0 )
   figure
   g = rectgrid( 20 );
   u = evalexpr( 'u', g.p, fea );
   fv.faces    = g.c';
   fv.vertices = [g.p' u];
   fv.facevertexcdata = u;
   fv.facecolor = 'interp';
   patch( fv )
   grid on, axis normal, view(3)
   xlabel( 'x' )
   ylabel( 'y' )
   title('Solution u')
 end


% Error checking.
 x = linspace( 0, 1, 11 );
 [x,y] = meshgrid(x,x);
 u = evalexpr( 'u', [x(:) y(:)]', fea );
 u_ref = l_poisol2( x(:), y(:), 6 );
 u_diff = u - u_ref;
 err = norm(u_diff);
 if( ~isempty(fid) )
   fprintf(fid,'\nL2  Error: %f\n',err)
   fprintf(fid,'\nL00 Error: %f\n',max(abs(u_diff)))
   fprintf(fid,'\n\n')
 end


 out.err  = err;
 out.tol  = opt.tol;
 out.pass = out.err<opt.tol;
 if ( nargout==0 )
   clear fea out
 end


%   -----------------------------------------------------------------------------