Description

EX_PERIODIC2 2D Periodic Poisson equation example.

[ FEA, OUT ] = EX_PERIODIC2( VARARGIN ) 2D Periodic Poisson equation example. Accepts the following property/value pairs.

Input       Value/{Default}        Description
-----------------------------------------------------------------------------------
isol        scalar 0/{1}           Stationary/Time dependent solution
iplot       scalar 0/{1}           Plot solution and error (=1)
                                                                                  .
Output      Value/(Size)           Description
-----------------------------------------------------------------------------------
fea         struct                 Problem definition struct
out         struct                 Output struct

Code listing

 cOptDef = { 'isol',     0;
             'iplot',    1;
             'fid',      1 };
 [got,opt] = parseopt(cOptDef,varargin{:});
 fid       = opt.fid;


% Unit square grid.
 fea.sdim = {'x' 'y'};
 fea.grid = rectgrid( 40 );


% Define Poisson equation with custom physics mode.
 fea = addphys( fea, @customeqn, { 'u' } );
 fea.phys.ce.eqn.seqn = '- ux_x - uy_y = f + u*eps';


% Boundary conditions.
 DirBC = true;
 NeuBC = ~DirBC;
 NeuBC_coef = @periodic_bc;   % Assign periodic_bc solve hook function as boundary coefficient.
 fea.phys.ce.bdr.coef = { 'bcnd_ce', '', '', {}, ...
                          { DirBC, NeuBC,      DirBC, NeuBC }, [], ...
                          { 0    , NeuBC_coef, 0,     0     } };

% Source term.
 fea.expr = { 'f' 'x*sin(5*pi*y) + exp(-((x-0.5)^2 + (y-0.5)^2)/0.02)' };


% Parse and solve problem.
 fea = parsephys( fea );
 fea = parseprob( fea );
 if( opt.isol==0 )
   fea.sol.u = solvestat( fea, 'fid', opt.fid );
 else
   fea.sol.u = solvetime( fea, 'tstop', -eps, 'fid', opt.fid );
 end


% Postprocessing.
 y  = linspace( 0, 1, 100 );
 ul = evalexpr( 'u', [zeros(1,100);y], fea );
 ur = evalexpr( 'u', [ones(1,100);y],  fea );
 if( opt.iplot>0 )
   subplot(1,2,1)
   postplot(fea, 'surfexpr', 'u', 'surfhexpr', 'u' )

% Plot solution on left and right sides.
   hold on
   fea.grid.p(1,:) = fea.grid.p(1,:) - 1;
   postplot(fea, 'surfexpr', 'u', 'surfhexpr', 'u' )
   fea.grid.p(1,:) = fea.grid.p(1,:) + 2;
   postplot(fea, 'surfexpr', 'u', 'surfhexpr', 'u' )

   subplot(1,2,2)
   plot( y, ul, 'k-' )
   hold on
   plot( y, ur, 'r.' )
   grid on
   xlabel( 'y' )
   ylabel( 'u' )
   legend( 'u(x=0)', 'u(x=1)', 'Location', 'South' )
 end


% Error checking.
 err = norm( ur - ul );

 out.err  = err;
 out.pass = err<1e-6;
 if ( nargout==0 )
   clear fea out
 end


%