FEATool Multiphysics
v1.17.0
Finite Element Analysis Toolbox
|
EX_POISSON1 1D Poisson equation example.
[ FEA, OUT ] = EX_POISSON1( VARARGIN ) Poisson equation on a line with a constant source term equal to 1, homogenous boundary conditions, and exact solution (-x^2+x)/2. Accepts the following property/value pairs.
Input Value/{Default} Description ----------------------------------------------------------------------------------- hmax scalar {1/10} Grid cell size sfun string {sflag1} Finite element shape function iphys scalar 0/{1} Use physics mode to define problem (=1) or directly define fea.eqn/bdr fields (=0) or use core assembly functions (<0) iplot scalar 0/{1} Plot solution (=1) . Output Value/(Size) Description ----------------------------------------------------------------------------------- fea struct Problem definition struct out struct Output struct
cOptDef = { 'hmax', 1/10; 'sfun', 'sflag1'; 'refsol', '(-x^2+x)/2'; 'fsrc', '1'; 'iphys', 1; 'icub', 2; 'iplot', 1; 'tol', 2e-2; 'fid', 1 }; [got,opt] = parseopt(cOptDef,varargin{:}); fid = opt.fid; % Grid generation. if( opt.hmax>0 ) nx = round( 1/opt.hmax ); fea.grid = linegrid( nx, 0, 1 ); else % Scrambled testing grid. fea.grid.p = [ 0 1/10 4/10 1/3 1 1-1/3 ]; fea.grid.c = [ 1 4 2 6 3 ; 2 3 4 5 6 ]; fea.grid.a = [ 0 2 1 4 3 ; 2 4 3 0 5 ]; fea.grid.b = [ 1 1 1 -1 ; 4 2 2 1 ]'; fea.grid.s = ones(1,5); end n_bdr = 2; % Number of boundaries. % Problem definition. fea.sdim = { 'x' }; % Coordinate name. switch opt.iphys case 0 % Directly define fea.eqn/bdr fields. fea.dvar = { 'u' }; % Dependent variable name. fea.sfun = { opt.sfun }; % Shape function. % Define equation system. fea.eqn.a.form = { [2;2] }; % First row indicates test function space (2=x-derivative), % second row indicates trial function space (2=x-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 = { opt.fsrc }; % Coefficient used in right hand side. % Define boundary conditions. if( strcmp(opt.sfun(end-1:end),'H3') ) % Prescribed derivatives at end points for Hermite elements. fea.bdr.d = {{ 0 0 ; 1/2 -1/2 }}; else fea.bdr.d = cell(1,n_bdr); [fea.bdr.d{:}] = deal(0); % Assign zero to all boundaries (homogenous Dirichlet conditions). end fea.bdr.n = cell(1,n_bdr); % No Neumann boundaries ('fea.bdr.n' empty). % 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. case 1 % Use physics mode. 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} = { opt.fsrc }; % 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( strcmp(opt.sfun(end-1:end),'H3') ) % Prescribed derivatives at end points for Hermite elements. fea.bdr.d = {{ 0 0 ; 1/2 -1/2 }}; 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. otherwise % Use core assembly functions. fea.dvar = { 'u' }; % Dependent variable name. fea.sfun = { opt.sfun }; % Shape function. fea = parseprob(fea); % Check and parse problem struct. % Assemble stiffness matrix. form = [2;2]; sfun = {opt.sfun;opt.sfun}; coefa = 1; sind = 1; i_cub = opt.icub; [vRowInds,vColInds,vAvals,n_rows,n_cols] = ... assemblea(form,sfun,coefa,i_cub,fea.grid.p,fea.grid.c,fea.grid.a,fea.grid.s,[]); A = sparse(vRowInds,vColInds,vAvals,n_rows,n_cols); % Check and compare with finite difference stencil. if (strcmp(opt.sfun,'sflag1')) h = 1/nx; n = nx+1; e = ones(n,1); A_ref = 1/h*spdiags([-e 2*e -e], -1:1, n, n); A_ref(1) = A_ref(1)/2; A_ref(end) = A_ref(end)/2; err = norm(A(:)-A_ref(:)); if err>opt.tol out.err = err; out.pass = -1; return end end form = 1; sfun = sfun{1}; coeff = 1; f = assemblef(form,sfun,coeff,i_cub,fea.grid.p,fea.grid.c,fea.grid.a,fea.grid.s,[]); % Check and compare with finite difference stencil. if (strcmp(opt.sfun,'sflag1')) f_ref = coeff*h*ones(n,1); f_ref([1 end]) = coeff*1/2*h; err = norm(f-f_ref); if err>1e-6 out.err = err; out.pass = -2; return end end % Set homogenous Dirichlet boundary conditions on first and last dof/node. bind = [1 nx+1]; A = A'; %' A(:,bind) = 0; % Zero out Dirichlet BC rows. for i=1:length(bind) % Loop to set diagonal entry to 1. i_a = bind(i); A(i_a,i_a) = 1; end A = A'; %' f(bind) = 0; % Set corresponding source term entries to Dirichlet BC values. % Solve problem. fea.sol.u = A\f; end % Postprocessing. if ( opt.iplot>0 ) x = linspace( 0, 1, 41 ); u = evalexpr( 'u', x, fea )'; figure subplot(3,1,1) plot( x, u ) axis( [0 1 0 0.2]) grid on title('Solution u') subplot(3,1,2) ux = (-x.^2+x)/2; plot( x, ux ) axis( [0 1 0 0.2]) grid on title('Exact solution') subplot(3,1,3) plot( x, abs(ux-u) ) title('Error') end % Error checking. xi = [1/2; 1/2]; s_err = ['abs(',opt.refsol,'-u)']; err = evalexpr0(s_err,xi,1,1:size(fea.grid.c,2),[],fea); ref = evalexpr0('u',xi,1,1:size(fea.grid.c,2),[],fea); err = sqrt(sum(err.^2)/sum(ref.^2)); if( ~isempty(fid) ) fprintf(fid,'\nL2 Error: %e\n',err) fprintf(fid,'\n\n') end out.err = err; out.tol = opt.tol; out.pass = out.err<out.tol; if ( nargout==0 ) clear fea out end