FEATool Multiphysics  v1.16.4 Finite Element Analysis Toolbox
ex_linearelasticity1.m File Reference

## Description

EX_LINEARELASTICITY1 Example for solid stress-strain on a cube.

[ FEA, OUT ] = EX_LINEARELASTICITY1( VARARGIN ) Example to calculate displacements and stresses for a cube stretched in axis-aligned directions.

Accepts the following property/value pairs.

Input       Value/{Default}        Description
-----------------------------------------------------------------------------------
E           scalar {1}             Modulus of elasticity
nu          scalar {0.3}           Poissons ratio
idir        scalar 1/{2,3}         Stress direction
igrid       scalar 1/{0}           Cell type (0=quadrilaterals, 1=triangles)
hmax        scalar {0.1}           Max grid cell size
sfun        string {sflag1}        Shape function for displacements
iplot       scalar 0/{1}           Plot solution (=1)
.
Output      Value/(Size)           Description
-----------------------------------------------------------------------------------
fea         struct                 Problem definition struct
out         struct                 Output struct


# Code listing

 cOptDef = { 'E',        1;
'nu',       0.3;
'force',    1;
'idir',     1;
'igrid',    0;
'hmax',     0.1;
'sfun',     'sflag1';
'iplot',    1;
'psecheck', 0;
'tol',      sqrt(eps);
'fid',      1 };
[got,opt] = parseopt(cOptDef,varargin{:});
fid       = opt.fid;

% Geometry definition.
fea.geom.objects = { gobj_block() };
fea.sdim = { 'x' 'y' 'z' };   % Coordinate names.

% Grid generation.
switch opt.igrid
case -1
fea.grid = blockgrid( round(1/opt.hmax) );
fea.grid = hex2tet( fea.grid );
case 0
fea.grid = blockgrid( round(1/opt.hmax) );
case 1
fea.grid = gridgen( fea, 'hmax', opt.hmax, 'fid', fid );
end
n_bdr = max(fea.grid.b(3,:));   % Number of boundaries.

% Boundary conditions.
s = sign(opt.idir); opt.idir = abs(opt.idir);
dtol     = opt.hmax;
fixbdr   = findbdr( fea, [fea.sdim{opt.idir},'<',num2str(dtol)] );     % "Right" boundary number.
forcebdr = findbdr( fea, [fea.sdim{opt.idir},'>',num2str(1-dtol)] );   % "Left" boundary number.

% Problem definition.
fea.phys.el.eqn.coef{1,end} = { opt.nu };
fea.phys.el.eqn.coef{2,end} = { opt.E  };
fea.phys.el.sfun            = { opt.sfun opt.sfun opt.sfun };

% Fix first boundary (set idir displacement to zero, Dirichlet BC).
bctype = mat2cell( zeros(3,n_bdr), [1 1 1], ones(1,n_bdr) );
bctype{opt.idir,fixbdr} = 1;
fea.phys.el.bdr.coef{1,5} = bctype;

% Apply load to second opposite boundary.
bccoef = mat2cell( zeros(3,n_bdr), [1 1 1], ones(1,n_bdr) );
bccoef{opt.idir,forcebdr} = s*opt.force;
fea.phys.el.bdr.coef{1,end} = bccoef;

% Parse and solve problem.
fea = parsephys( fea );
fea = parseprob( fea );
fea.sol.u = solvestat( fea, 'fid', fid, 'icub', 1+str2num(strrep(opt.sfun,'sflag','')) );

% Postprocessing.
if( opt.iplot>0 )
figure
postplot( fea, 'surfexpr', fea.dvar{opt.idir} )
title( [fea.sdim{opt.idir},'-displacement'] )
end

% Error checking.
x = linspace(0,1,7);
x = x(2:end-1)';
[x,y,z] = ndgrid(x,x,x);
p = [x(:) y(:) z(:)]';
ui = evalexpr( fea.dvar{opt.idir}, p, fea );
ui_ref = s*p(opt.idir,:)';
out.ui_err = norm( ui - ui_ref );
out.pserr  = [];
out.peerr  = [];
for i=1:6
out.serr(i) = norm( abs( evalexpr( fea.phys.el.eqn.vars{5+i, end}, p, fea ) - s*(i==opt.idir) ) );
out.eerr(i) = norm( abs( evalexpr( fea.phys.el.eqn.vars{14+i,end}, p, fea ) - s*(i==opt.idir) + s*0.3*(i<=3 & i~=opt.idir)) ) ;
if( opt.psecheck && i<=3 )
out.pserr(i) = norm( abs( evalexpr( fea.phys.el.eqn.vars{11+i,end}, p, fea ) - (i==1) ) );
out.peerr(i) = norm( abs( evalexpr( fea.phys.el.eqn.vars{20+i,end}, p, fea ) - (i==1) + 0.3*(i~=1) ) );
end
end
out.serr = [ out.serr norm( abs( evalexpr( fea.phys.el.eqn.vars{1, end}, p, fea ) - 1 ) ) ];
out.pass = all( [out.ui_err out.serr out.eerr out.pserr out.peerr] < opt.tol );

if ( nargout==0 )
clear fea out
end