FEATool Multiphysics  v1.16.6 Finite Element Analysis Toolbox
ex_linearelasticity5.m File Reference

## Description

EX_LINEARELASTICITY5 Vibration of a square plate.

[ FEA, OUT ] = EX_LINEARELASTICITY5( VARARGIN ) Vibration of a square plate (NAFEMS FV52 Benchmark).

Reference

[1] National Agency for Finite Element Methods and Standards. The Standard NAFEMS Benchmarks. Rev. 3. United Kingdom: NAFEMS, October 1990.

Accepts the following property/value pairs.

Input       Value/{Default}        Description
-----------------------------------------------------------------------------------
igrid       scalar 0/{1}           Cell type (>0=hexahedral, <0=tetrahedral)
sfun        string {sflag1}        Shape function
iplot       scalar 0/{1}           Plot solution (=1)
.
Output      Value/(Size)           Description
-----------------------------------------------------------------------------------
fea         struct                 Problem definition struct
out         struct                 Output struct


# Code listing

 cOptDef = { 'igrid',    2;
'sfun',     'sflag1';
'iplot',    1;
'tol',      0.02;
'fid',      1 };
[got,opt] = parseopt(cOptDef,varargin{:});
fid       = opt.fid;

E   = 200e9;
nu  = 0.3;
rho = 8000;

% Geometry definition.
fea.sdim = {'x' 'y' 'z'};
gobj = gobj_block( 0, 10, 0, 10, -0.5, 0.5, 'B1' );
fea.geom.objects = { gobj };

fea.grid = blockgrid(abs(opt.igrid)*10,abs(opt.igrid)*10,abs(opt.igrid),[0,10;0,10;-0.5,0.5]);
if( opt.igrid<0 )
fea.grid = hex2tet( fea.grid );
end

% Equations and problem definition.
fea = addphys( fea, @linearelasticity );
fea.phys.el.eqn.coef{1,end} = { nu  };
fea.phys.el.eqn.coef{2,end} = { E   };
fea.phys.el.eqn.coef{3,end} = { rho };
fea.phys.el.sfun            = { opt.sfun, opt.sfun, opt.sfun };

% Set/constrain w = 0 on lower edges (z = -0.5).
edg = [];
[be,e,ev] = gridbdre( fea.grid.b, fea.grid.c );
n_bdre = max(be(end,:));
for i_bdre=1:n_bdre
ix = be(end,:) == i_bdre;
ie = be(4,ix);
iv = unique([ev(ie,1);ev(ie,2)]);
pz = fea.grid.p(3,iv);
if( all( pz <= -0.5+sqrt(eps)) )
edg_i.type  = 'constraint';
edg_i.index = i_bdre;
edg_i.dvar  = 3;
edg_i.expr  = 0;

edg = [ edg, edg_i ];
end
end
fea.edg = edg;

% Solve problem.
fea = parsephys( fea );
fea = parseprob( fea );

[fea.sol.u,fea.sol.l] = solveeig( fea, 'neigs', 10, 'fid', fid );

% Postprocessing.
if( opt.iplot>0 )
postplot( fea, 'surfexpr', 'sqrt(u^2+v^2+w^2)', 'solnum', 4 )
end

out = [];
f = sqrt(max(0,fea.sol.l))/(2*pi);
f_ref = [0;0;0;45.897;109.44;109.44;167.89;193.59;206.19;206.19];
out.err  = norm(f_ref-f)/norm(f_ref);
out.pass = out.err < opt.tol;

if( nargout==0 )
clear fea out
end