FEATool
v1.6
Finite Element Analysis Toolbox

A selection of tutorial models and examples are presented in this section.
FEATool supports modeling heat transfer through both conduction, that is heat transported by a diffusion process, and also convection, which is heat transported through a fluid through convection by a velocity field. The heat transfer physics mode supports both these processes, and defines the following equation
where is the density, the heat capacity, is the thermal conductivity, heat source term, and a vector valued convective velocity field.
This example models heat conduction in the form of transient cooling for shrink fitting of a two part assembly. A tungsten rod heated to 84 C is inserted into a 10 C chilled steel frame part. The time when the maximum temperature has cooled to 70 C should be determined. The assembly is cooled due to convection through a surrounding medium kept at T_{inf} = 17 C and a heat transfer coefficient of h = 750 W/m^{2}K. The surrounding cooling medium is not modeled directly, thus the convective term is omitted, but the effects are incorporated into the model by the use of natural convection boundary conditions [1].
This section describes how to set up and solve the thermal shrink fitting example with the FEATool graphical user interface (GUI).
cd C:\featool
In the command window type
featool
to start the graphical user interface (GUI).
To see this more clearly we can process the simulation data on the command line by using the Export > FEA Problem Struct To Main Workspace from the File menu. Then switch to the main command window and use the following script commands to see that T = 70 is reached at about t = 11.
tlist = fea.sol.t; for i=1:numel(tlist) T_min(i) = min(fea.sol.u(:,i)); T_max(i) = max(fea.sol.u(:,i)); end ix = find(T_max<70); i1 = ix(1); i2 = i1  1; s = ( T_max(i2)  70 )/( T_max(i2)  T_max(i1) ); t_70 = tlist(i2) + s*( tlist(i1)  tlist(i2) ) u_70 = fea.sol.u(:,i2) + ... s*( fea.sol.u(:,i1)  fea.sol.u(:,i2) ) figure plot( tlist, T_min, 'b' ) hold on plot( tlist, T_max, 'r' ) grid on title('Maximum and minimum temperatures') ylabel('Temperature [C]') xlabel('Time [s]')
How to set up and solve the thermal shrink fitting problem on the command line interface is illustrated in the ex_heattransfer5 script file which can be found in the FEATool examples directory.
[1] Krysl P. A Pragmatic Introduction to the Finite Element Method for Thermal and Stress Analysis. Pressure Cooker Press, USA, 2005.
A well known benchmark, test, and validation problem suite for incompressible fluid flows are the DFG cylinder benchmark problems [2]. Although it is not possible to derive analytical solutions to these test cases, accurate numerical solutions to benchmark reference quantities have been established for a number of configurations [3], [4].
The test configuration used in the following places a solid cylinder centered at (0.2, 0.2) with diameter in a by rectangular channel. The fluid density is taken as 1 and the viscosity is 0.001. A fully developed parabolic velocity profile is prescribed at the inlet, , with the maximum velocity . This results in a mean velocity and a Reynolds number and thus the flow field will be laminar.
As the fluid is considered incompressible the problem is governed by the NavierStokes equations. That is,
where in this case the time dependent term can be neglected. The benchmark quantities that should be computed include the pressure difference between the front and rear of the cylinder , and the coefficients of drag and lift , defined as
The drag and lift forces, and , can be computed as
where is the velocity in the tangential direction .
This section describes how to set up and solve the Flow Around a Cylinder example with the FEATool graphical user interface (GUI).
cd C:\featool
In the command window type
featool
to start the graphical user interface (GUI).
The process to set up and solve the flow around cylinder test problem on the command line interface is illustrated in the ex_navierstokes3 script file which can be found in the examples directory.
In addition to the stationary test case described above. An instationary benchmark test case is also available (ex_navierstokes6). This test case uses the same geometry but instead applies an inflow condition that varies with time (u_inflow = 6*sin(pi*t/8)*(y*(0.41y))/0.41^2) so that 0<Re(t)<100 [5]. Computations with FEATool show that the drag and lift coefficients, and pressure difference between front and rear of the cylinder agrees very well with the reference values [6].
[2] Schäfer M, Turek S. Benchmark computations of laminar flow around cylinder. Flow Simulation with High–Performance Computers II (Notes on Numerical Fluid Mechanics), 52, 547–566, Vieweg, 1996.
[3] Nabh G. On higher order methods for the stationary incompressible NavierStokes equations. PhD Thesis 1998; Universität Heidelberg. Preprint 42/98, 1998.
[4] John V, Matthies G. Higherorder finite element discretizations in a benchmark problem for incompressible flows. International Journal for Numerical Methods in Fluids 2001; 37(8):885–903, DOI: 10.1002/fld.195.
[5] John V. Reference values for drag and lift of a twodimensional timedependent flow around a cylinder. International Journal for Numerical Methods in Fluids 2004; 44:777788.
[6] John V, Rang J. Adaptive time step control for the incompressible Navier–Stokes equations. Comput. Methods Appl. Mech. Engrg. 199 (2010) 514–524.
The Equation Editing Example  Axisymmetric Fluid Flow example in the quickstart guide shows how to edit predefined equations in the physics modes to model rotationally symmetric or axisymmetric flow in a narrowing pipe.
Flow over a backwards facing step is another classic computational fluid dynamics test problem which has in various forms been used quite extensively for validation of simulation codes. The test problem essentially consists of studying how a fully developed flow profile reacts to a sudden expansion or step in a channel. The expansion will cause a break in the flow and a recirculation or separation zone will form. To measure and compare results the resulting length of the recirculation or separation zone is most often used.
The stationary incompressible NavierStokes equations are to be applied corresponding to a Reynold number of . The size domain is chosen according to the data below. The inlet velocity function is given as where is the channel height, and maximum velocity Noslip zero velocity conditions are applied to all walls and a suitable outflow condition must also be applied. The reference recirculation zone length found in reference [7] and [8] is estimated to be 7.93 length units.
This section describes how to set up and solve the Flow over a Backwards Facing Step example with the FEATool graphical user interface (GUI) and using the FeatFlow external CFD Solver.
cd C:\featool
In the command window type
featool
to start the graphical user interface (GUI).
The process to set up and solve the flow over a backwards facing step test problem on the command line interface is illustrated in the ex_navierstokes4 script file which can be found in the examples directory.
[7] P.M. Gresho and R.L. Sani, Incompressible Flow and the Finite Element Method, Volume 1 & 2, John Wiley & Sons, New York, 2000.
[8] A. Rose and B. Simpson: Laminar, ConstantTemperature Flow Over a Backward Facing Step, 1st NAFEMS Workbook of CFD Examples, Glasgow, UK, 2000.
An example of setting up and solving the Thin Plate with Hole example with the FEATool graphical user interface (GUI) is described in the FEATool quickstart guide Structural Mechanics Example  Thin Plate with Hole.
In this section a three dimensional model of a bracket with a hole will be described. The vertical side of the bracket is fixed to a wall (thus zero displacement boundary conditions are suitable here), while the outer end of the bracket is loaded with a negative force of 10000 in the zdirection. The resulting displacements will be studied together with creating a displacement plot visualizing the deformation.
This section describes how to set up and solve the bracket stress problem with the FEATool graphical user interface (GUI).
cd C:\featool
In the command window type
featool
to start the graphical user interface (GUI).
scaling = 5e3; dp = zeros(size(fea.grid.p)); for i=1:3 dp(i,:) = scaling*evalexpr( fea.dvar{i}, ... fea.grid.p, fea ); end fea_disp = fea; fea_disp.grid.p = fea_disp.grid.p + dp; f = figure; left_cells = selcells(fea,['x>',num2str(1.5*0.02)]); plotgrid( fea, 'facecolor', [.95 .95 .95], ... 'edgecolor', [.8 .8 1], ... 'selcells', left_cells, ... 'parent', gca ) hold on plotgrid( fea_disp, 'parent', gca ) title( 'Displacement plot' ) view( [30 20] )
The process to set up and solve the bracket stress test problem on the command line interface is illustrated in the ex_linearelasticity2 script file which can be found in the examples directory.
Parametric studies can be very useful tools to clearly see how a model will behave under varying conditions. The following section presents a parametric study of the bracket model varying both the geometry, in this case plate thickness, and also the applied load force.
thickness = [ 0.02 0.015 0.01 0.008 ]; loads = [ 1e4 2e4 1e5 1e6 ];
for i=1:length(thickness) % Geometry operations. fea.geom = struct; gobj = gobj_block( 0, 0.02, 0, 0.2, 0, 0.2, 'B1' ); fea = geom_add_gobj( fea, gobj ); gobj = gobj_block( 0, 0.2, 0, 0.2, ... 0.1thickness(i)/2, 0.1+thickness(i)/2, 'B1' ); fea = geom_add_gobj( fea, gobj ); gobj = gobj_cylinder( [ 0.1, 0.1, 0.08 ], 0.08, 0.04,3,'C1'); fea = geom_add_gobj( fea, gobj ); fea = geom_apply_formula( fea, 'B1+B2C1' ); % Grid generation. fea.grid = gridgen( fea, 'hmax', thickness(i)/2 );
As second loop can now be introduced where the varying load is extracted and prescribed by using a constant variable. The equation settings can be kept from the saved model.
for j=1:length(loads) % Constants and expressions. fea.expr = { 'load', loads(j) }; % Equation settings. fea.phys.el.dvar = { 'u', 'v', 'w' }; fea.phys.el.sfun = { 'sflag1', 'sflag1', 'sflag1' }; fea.phys.el.eqn.coef = { 'nu_el', ', '', { '0.3' }; ...
Two loops have now been set up around the model, the first creating a geometry with varying thickness, and the second will prescribe the varying load. Note that setting up the loops in this way allows reusing the grid in the inner loop and saving computational time.
% Boundary conditions. n_bdr = max(fea.grid.b(3,:)); ind_fixed = findbdr( fea, 'x<0.005' ); ind_load = findbdr( fea, 'x>0.19' ); fea.phys.el.bdr.sel = ones(1,n_bdr); % Define Homogeneous Neumann BC type. n_bdr = max(fea.grid.b(3,:)); bctype = num2cell( zeros(3,n_bdr) ); % Set Dirichlet BCs for the fixed boundary. [bctype{:,ind_fixed}] = deal( 1 ); fea.phys.el.bdr.coef{1,end2} = bctype; % Define all zero BC coefficients. bccoef = num2cell( zeros(3,n_bdr) ); % Apply negative zload to outer edge. [bccoef{3,ind_load}] = deal( 'load' ); % Store BC coefficient array. fea.phys.el.bdr.coef{1,end} = bccoef;
% Solver call. fea = parsephys( fea ); fea = parseprob( fea ); fea.sol.u = solvestat( fea ); % Postprocessing. s_vm = fea.phys.el.eqn.vars{1,2}; vm_stress = evalexpr( s_vm, fea.grid.p, fea ); max_vm_stress(i,j) = max(vm_stress); end,end
% Visualization. [t,l] = meshgrid(thickness,loads); surf( t, l, max_vm_stress, log(max_vm_stress) ) xlabel( 'Thickness' ) ylabel( 'Load' ) zlabel( 'Maximum stress' ) view( 45, 30 )
The complete parameter study model can be found as the ex_linearelasticity3 script file in the examples directory.
This section contains an example how to easily define and couple multiphysics models in FEATool. The chosen example is a benchmark model by De Vahl [9], [10] which simulates natural convection in a unit square through the Boussinesq approximation. The model consists of a NavierStokes equations physics mode representing the fluid flow with solid wall or noslip boundary conditions everywhere. In addition a heat transfer physics mode is added for the temperature. The top and bottom boundaries are perfectly insulated while the left boundary has a temperature of 1 and the right zero.
The physics modes are two way coupled through the vertical source term in the NavierStokes equations, Pr*Ra*T (here the Prandtl and Rayleigh numbers are, Pr=0.71 and Ra=1e3 respectively), and the velocities transporting the temperature coming directly from the fluid flow. The references contain benchmark reference and comparison results for a number of quantities such as maximum velocities and the Nusselt number.
In the following it is described how to set up and solve a multiphysics natural convection benchmark problem the FEATool graphical user interface (GUI).
cd C:\featool
In the command window type
featool
to start the graphical user interface (GUI).
The lower command prompt can be used to input Octave and Matlab commands for advanced processing (or alternatively the whole FEATool fea struct can be exported to the main workspace for further processing). The following code shows how to calculate the maximum velocity and position along the vertical center line
y = [0:0.01:1]; x = 0.5*ones(size(y)); u_eval = evalexpr( 'u', [x;y], fea ); [u_max,ix] = max( u_eval ) y_max = y( ix )
The resulting values of u_max and y_max are here 2.6819 and 0.81, respectively, which is a little low compared to the reference solution (3.649 at a height of 0.813) which can indicate the the grid resolution is somewhat too coarse.
The process to set up and solve the natural convection problem on the command line interface is illustrated in the ex_natural_convection script file which can be found in the examples directory.
[9] D. de Vahl Davis, Natural Convection of Air in a Square Cavity: A Bench Mark Solution, Int. J. Numer. Meth. Fluids, vol. 3, pp. 249264, 1983.
[10] D. de Vahl Davis and I. P. Jones, Natural Convection of Air in a Square Cavity: A Comparison Exercise, Int. J. Numer. Meth. Fluids, vol. 3, pp. 227248, 1983.
An example of setting up and modeling a heat exchanger with coupled fluid flow and temperature with both free and forced convection is described in the FEATool quickstart guide Multiphysics Example  Heat Exchanger.
This example describes how to model piezoelectric bending of a beam on the command line as an mscript model in FEATool. The test case comes from a research paper by Hwang and Park [11]. For this problem, the following coupled system of equations should be solved
where is a stress tensor, represents body forces, electric displacement, and distributed charges.
By using the constitutive relations for a piezoelectric material [11], the stresses can in two dimensions be converted to the following form
where is an array of constitutive relations and material parameters. This leads to a system of three equations for the unknown displacements and , and the potential .
There are no readily available physics modes in FEATool yet so to model this we are defining the problem from scratch. In this case we will be looking at composition of two 1.2 cm long and 2 mm thick beams with opposite polarization.
First we assign names to the space coordinates (x and y), and create a grid. Note, that the shape is so simple here that we don't need to explicitly define a geometry, we can use the 'rectgrid' command directly
fea.sdim = { 'x' 'y' }; fea.grid = rectgrid( 20, 20, [ 0 12e3; 1e3 1e3] );
Next is to define the equation coefficients and constitutive relations. See the attached file for their definitions, one thing to note that to define the polarization the coefficient is multiplied with a switch expression 2*(y<0)1 which evaluates to 1 in the top half and 1 in the bottom (here we have also used the name y for the space coordinate in the 2nd direction as we defined earlier).
Now we define the dependent variables and assign 2nd order Lagrange finite element shape functions for them all
fea.dvar = { 'u' 'v' 'V' }; fea.sfun = { 'sflag2' 'sflag2' 'sflag2' };
The material parameters and constitutive relations are defined as
% Equation coefficients. Emod = 2e9; % Modulus of elasticity nu = 0.29; % Poissons ratio Gmod = 0.775e9; % Shear modulus d31 = 0.22e10; % Piezoelectric sn coefficient d33 = 0.3e10; % Piezoelectric sn coefficient prel = 12; % Relative electrical permittivity pvac = 0.88541878176211; % Elec. perm. of vacuum % Constitutive relations. constrel = [ Emod/(1nu^2) nu*Emod/(1nu^2) 0 ; nu*Emod/(1nu^2) Emod/(1nu^2) 0 ; 0 0 Gmod ]; piezoel_st = [ 0 d31 ; 0 d33 ; 0 0 ]; piezoel_sn = constrel*piezoel_st; dielmat_st = [ prel 0 ; 0 prel ]*pvac ; dielmat_sn = [ dielmat_st  piezoel_st'*piezoel_sn ]; % Populate coefficient matrices (negative sign % due to fem partial integration). c{1,1} = { constrel(1,1) constrel(1,3) ; constrel(1,3) constrel(3,3) }; c{1,2} = { constrel(1,3) constrel(1,2) ; constrel(3,3) constrel(2,3) }; c{2,1} = c{1,2}'; c{2,2} = { constrel(3,3) constrel(1,3) ; constrel(1,3) constrel(2,2) }; c{1,3} = { piezoel_sn(1,1) piezoel_sn(1,2) ; piezoel_sn(3,1) piezoel_sn(3,2) }; for i=1:4 c{1,3}{i} = [num2str(c{1,3}{i}),'*(2*(y<0)1)']; end c{3,1} = c{1,3}'; c{2,3} = { piezoel_sn(3,1) piezoel_sn(3,2); piezoel_sn(2,1) piezoel_sn(2,2) }; for i=1:4 c{2,3}{i} = [num2str(c{2,3}{i}),'*(2*(y<0)1)']; end c{3,2} = c{2,3}'; c{3,3} = { dielmat_sn(1,1) dielmat_sn(2,1) ; dielmat_sn(2,1) dielmat_sn(2,2) };
To define the finite element bilinear forms we use the following code
bilinear_form = [ 2 2 3 3 ; 2 3 2 3 ]; for i=1:length(fea.dvar) for j=1:length(fea.dvar) fea.eqn.a.form{i,j} = bilinear_form; fea.eqn.a.coef{i,j} = c{i,j}(:)'; end end
where each bilinear form has four terms, the top line in the definition defines the finite element shape for the dependent variable, and second line for the test function space (a 2 indicates xderivative, and 3 yderivative).
In this case all three source terms are zero, thus
fea.eqn.f.form = { 1 1 1 }; fea.eqn.f.coef = { 0 0 0 };
The boundary conditions are defined as follows
n_bdr = max(fea.grid.b(3,:)); fea.bdr.d = cell(3,n_bdr); fea.bdr.n = cell(3,n_bdr); i_top = 3; i_bottom = 1; i_left = 4; [fea.bdr.d{3,i_top}] = deal( 200 ); [fea.bdr.d{3,i_bottom}] = deal( 0 ); [fea.bdr.d{1:2,i_left}] = deal( 0 );
In this way we have set a potential difference of 200, and fixed the left hand side boundary. Now we parse the problem struct and solve the system with
fea = parseprob( fea ); fea.sol.u = solvestat( fea );
After the problem has been solved we can for example visualize the displacement in the ydirection with the following commands
plotsubd( fea, 'labels', 'off', 'setaxes', 'off' ) u = evalexpr( 'u', fea.grid.p, fea ); v = evalexpr( 'v', fea.grid.p, fea ); DSCALE = 5000; fea_disp = fea; fea_disp.grid.p = fea.grid.p + [ u v ]'*DSCALE; postplot( fea_disp, 'surfexpr', 'v', 'isoexpr', 'v' )
An example mscript file for this model can be found in the ex_piezoelectric1 mscript file.
[11] Hwang, W. S., Park, H. C. Finite Element Modeling of Piezoelectric Sensors and Actuators. AIAA Journal, 31(5), pp 930937, 1993.
This section contains an example how to easily define a custom equation or model in FEATool. In this case a BlackScholes model equation is studied which is often used in financial analytics. The equation to be modeled reads
with boundary conditions and on the left and right sides of the domain, respectively. Initial condition is taken as For this problem an exact analytical solution exists, namely . Consult the Custom Equation section for details on the equation syntax.
This section describes how to set up and solve the custom equation example with the FEATool graphical user interface (GUI).
cd C:\featool
In the command window type
featool
to start the graphical user interface (GUI).
The process to set up and solve the custom equation problem on the command line interface is illustrated in the ex_custom_equation1 script file which can be found in the examples directory.
In the quickstart guide there is described how to set up a model of the wave equation with the custom equations physics mode (Custom Equation Example  Wave Equation on a Circle). In this case the custom equation uses two equations/dependent variables do describe the system.
The classic Poisson equation is one of the most fundamental partial differential equations (PDEs). Although one of the simplest equations, it is a very good model for the process of diffusion and comes up again and again in many applications (for example fluid flow, heat transfer, and chemical transport). It is thus very fundamental to many simulation codes to be able to solve it correctly and efficiently.
This example shows how to up and solve the Poisson equation
for a scalar field on a circle with radius in two dimensions. Both the diffusion coefficient and right hand side source term are assumed constant and equal to 1. The Poisson problem is also considered stationary meaning the time dependent term can be neglected. With these assumptions equation (1) simplifies to
Moreover, homogeneous Dirichlet boundary conditions are prescribed on all boundaries of the domain, that is . The exact solution for this problem is which can be used to measure the accuracy of the computed solution.
This section describes how to set up and solve the Poisson equation (1) with the FEATool graphical user interface (GUI) which is available when using FEATool together with Octave version 4.0 or later and Matlab.
cd C:\featool
In the command window type
featool
to start the graphical user interface (GUI).
This section describes how to set up and solve the Poisson equation (1) with the Octave or Matlab command line interface (CLI). To make it easier and faster for users to set up model problems several physics modes have been predefined with equations and boundary conditions for physics such as convection, diffusion, and reaction of chemical species, heat transport through convection and conduction, and incompressible fluid flow with the NavierStokes equations. In this example the Poisson physics mode will be used.
Start Octave/Matlab and make sure you have either run the installation script to set up the FEATool directory paths or alternatively run
addpaths
on the command line from the FEATool installation directory before you start.
r = 1; % Circle radius. xc = 0; % xcenter coordinate. yc = 0; % ycenter coordinate. gobj = gobj_circle( [xc yc], r ); fea.geom.objects = { gobj };
hmax = 0.1; % Maximum grid cell size. fea.grid = gridgen( fea, 'hmax', hmax );
fea.sdim = { 'x' 'y' };
fea = addphys( fea, @poisson );
fea.phys.poi.eqn
fea.phys.poi.eqn.coef{2,4} = {1}; % Diffusion coefficient. fea.phys.poi.eqn.coef{3,4} = {1}; % Source term.
fea.phys.poi.bdr.coef
% Use first boundary condition specification. fea.phys.poi.bdr.sel = 1; % Set Dirichlet coefficient to zero. fea.phys.poi.bdr.coef{1,7} = {0}; % Set flag to use Dirichlet BC for first selection. fea.phys.poi.bdr.coef{1,5} = {1}; % Set flag not to use Neumann BC for first selection. fea.phys.poi.bdr.coef{2,5} = {0};
fea = parsephys( fea ); fea = parseprob( fea );
fea.sol.u = solvestat(fea);
postplot( fea, 'surfexpr', 'u', ... 'isoexpr', 'u', ... 'arrowexpr', {'ux','uy'} )
u_diff = 'u(1x^2y^2)/4'; clf postplot( fea, 'surfexpr', u_diff )
This section describes how to set up and solve the Poisson equation (1) with the Octave or Matlab command line interface (CLI) by directly setting up the necessary equation and boundary definition fields (that is without using the predefined Poisson physics mode).
Start Octave/Matlab and make sure you have either run the installation script to set up the FEATool directory paths or alternatively run
addpaths
on the command line from the FEATool installation directory before you start.
r = 1; % Circle radius. xc = 0; % xcenter coordinate. yc = 0; % ycenter coordinate. gobj = gobj_circle( [xc yc], r ); fea.geom.objects = { gobj };
hmax = 0.1; % Maximum grid cell size. fea.grid = gridgen( fea, 'hmax', hmax );
fea.sdim = { 'x' 'y' };
fea.dvar = { 'u' };
fea.sfun = { 'sflag1' };
% Define bilinear form. The first row indicates test function % space, and second row trial function space. Each column % defines a term in the bilinear form and 1 corresponds to % a basis function evaluation, 2 = xderivative, % 3 = yderivative, (d 4 = zderivative in 3D). fea.eqn.a.form = { [2 3; 2 3] }; % Define coefficients used in assembling the bilinear forms. fea.eqn.a.coef = { [1 1] }; % Test function space to evaluate in linear form. fea.eqn.f.form = { 1 }; % Coefficient value used in assembling the linear form. fea.eqn.f.coef = { 1 };
% Number of boundary segments. n_bdr = max(fea.grid.b(3,:)); % Allocate space for n_bdr boundaries. fea.bdr.d = cell(1,n_bdr); % Assign u=0 to all Dirichlet boundaries. [fea.bdr.d{:}] = deal(0); % No Neumann boundaries (fea.bdr.n{i} empty). fea.bdr.n = cell(1,n_bdr);
fea = parseprob( fea );
fea.sol.u = solvestat( fea );
postplot( fea, 'surfexpr', 'u', ... 'isoexpr', 'u', ... 'arrowexpr', {'ux','uy'} )
u_diff = 'u(1x^2y^2)/4'; clf postplot( fea, 'surfexpr', u_diff )
This final section shows how to define and solve the Poisson equation (1) by directly calling the core finite element library functions. The low level functions are used to assemble the system matrix and right hand side vector to which boundary conditions are applied, after which the system can be solved.
Start Octave/Matlab and make sure you have either run the installation script to set up the FEATool directory paths or alternatively run
addpaths
on the command line from the FEATool installation directory before you start.
r = 1; % Circle radius. xc = 0; % xcenter coordinate. yc = 0; % ycenter coordinate. gobj = gobj_circle( [xc yc], r ); fea.geom.objects = { gobj };
hmax = 0.1; % Maximum grid cell size. fea.grid = gridgen( fea, 'hmax', hmax );
% Bilinear form to assemble (u_x*v_x+u_y*v_y). form = [2 3;2 3]; % FEM shape functions used in test and trial function spaces % (here first order conforming P1 Lagrange functions). sfun = {'sflag1';'sflag1'}; % Coefficients to use for each term in the bilinear form. coef = [1 1]; i_cub = 3; % Numerical quadrature rule to use. [vRowInds,vColInds,vAvals,n_rows,n_cols] = ... assemblea( form, sfun, coef, i_cub, ... fea.grid.p, fea.grid.c, fea.grid.a ); % Construct sparse matrix. A = sparse( vRowInds, vColInds, vAvals, n_rows, n_cols );
% Linear form to assemble (1 = evaluate function values). form = [1]; % Finite element shape function. sfun = {'sflag1'}; % Coefficient to use in the linear form. coef = [1]; % Numerical quadrature rule to use. i_cub = 3; f = assemblef( form, sfun, coef, i_cub, ... fea.grid.p, fea.grid.c, fea.grid.a );
bind = []; for i_b=1:size(fea.grid.b,2) % Loop over boundary edges. i_c = fea.grid.b(1,i_b); % Cell number. i_edge = fea.grid.b(2,i_b); j_edge = mod( i_edge, size(fea.grid.c,1) ) + 1; bind = [bind fea.grid.c([i_edge j_edge],i_c)']; end bind = unique(bind);
A = A'; % Transpose matrix for faster row elimination. A(:,bind) = 0; % Zero/remove entire rows. for i=1:length(bind) % Set diagonal entries of Dirichlet % boundary condition rows to 1. i_a = bind(i); A(i_a,i_a) = 1; end A = A'; %' f(bind) = 0; % Set corresponding RHS entries % to the prescribed boundary values.
u = A\f;
fea.sdim = { 'x' 'y' }; fea.dvar = { 'u' }; fea.sfun = { 'sflag1' }; fea = parseprob( fea );
fea.sol.u = u; postplot( fea, 'surfexpr', 'u' )
u_diff = u  ( 1  fea.grid.p(1,:)'.^2  ... fea.grid.p(2,:)'.^2 )/4; clf fea.sol.u = u_diff; postplot( fea, 'surfexpr', 'u' )
The Classic Equation Example  Poisson Equation with a Point Source example in the quickstart guide shows how to model a circle with a point source at its center.
The following example illustrates how to set up and solve the Poisson equation on a unit sphere ( ) with homogeneous boundary conditions, that is with the FEATool FEM assembly primitives.
To derive the finite element discretization of the general Poisson equation first multiply the equation with an arbitrary function (called test function), and integrate over the whole domain , that is . By applying the Gauss theorem or partial integration to the left side we get
The boundary term can be neglected assuming that we only have prescribed or fixed value (Dirichlet) boundary conditions. In 3D the equation will look like
By discretizing the solution (trial function space) variable as and similarly test function , where are the finite element shape or basis functions (usually taken the same for and in the Galerkin approximation), and is the value of the solution at node . Inserting this we get
which after discretization and integration with numerical quadrature gives us a matrix system to solve
To start implementing the model pragmatically we first have to create a grid. Instead of creating a grid by using a geometry object we use a geometric sphere grid primitive directly
grid = spheregrid(); p = grid.p; % Extract grid points. c = grid.c; % Extract grid cells. a = grid.a; % Extract adjacency information.
To assemble the matrix A we can identify three bilinear terms (with a coefficient of one) . To specify this in FEATool first create a cell array with coefficient values (or expressions) for each term, in this case coef = {1 1 1};
. Function values and/or derivatives for the trial and test functions are specified in a 2 by nterms vector, here form = [2 3 4;2 3 4];
where the first row specifies trial functions to use for the three terms (1=function value, 2=xderivative, 3=yderivative, 4=zderivative), and the bottom row the test functions. Lastly, the actual shape or basis functions are specified in a cell array as sfun = {'sflag1';'sflag1'};
, which prescribes first order linear Lagrange shape functions for both the trial and test function spaces. Given a grid with grid points p and cell connectivity c the call to the fem assembly routine looks like
coef = { 1 1 1 }; form = [ 2 3 4 ; ... 2 3 4 ]; sfun = { 'sflag1' ; ... 'sflag1' }; i_cub = 3; % Numerical quadrature rule. cind = 1:size(c,2); % Index to cells to assemble. [i,j,s,m,n] = assemblea( form, sfun, coef, i_cub, p, c, a ); A = sparse( i, j, s, m, n );
Similarly, the right hand side linear form or load vector b can be assembled as
coef = { 1 }; form = [ 1 ]; sfun = { 'sflag1' }; i_cub = 3; % Numerical quadrature rule. cind = 1:size(c,2); % Index to cells to assemble. b = assemblef( form, sfun, coef, i_cub, p, c, a );
Prescribing homogeneous (zero value) Dirichlet boundary conditions on the boundary of the domain can be done by first finding the indexes to the nodes on the boundary (in this case the unit sphere), and then setting the corresponding values in the load vector to zero. Moreover, the corresponding rows in the matrix A must also be modified so that they are all zero except for the diagonal entries which are set to one (this preserves the boundary values specified in the b vector).
radi = sqrt(sum(p.^2,1)); % Radius of grid points. bind = find(radi>1sqrt(eps)); % Find grid points % on the boundary. b(bind) = 0; % Set 'bind' rhs entries to zero. A(bind,:) = 0; % Set 'bind' rows to zero. for i=1:length(bind) % Set 'bind' diagonal entries to 1 i_a = bind(i); A(i_a,i_a) = 1; end
The problem is now fully specified and can be solved as
u = A\b;
where u now contains the values of the solution in the grid points/nodes. To visualize the solution we put everything in a FEATool fea struct
fea.grid = grid; fea.dvar = { 'u' }; fea.sfun = sfun; fea.sol.u = u; fea = parseprob( fea ); postplot( fea, 'surfexpr', 'u', 'selcells', 'y>0.1' )
The following selection of command line examples and test cases can be found in the examples directory.
Field  Description 

ex_convdiff1  Convection and diffusion on a unit square. 
ex_convdiff2  Periodic convection and diffusion on a line with exact solution. 
ex_convdiff3  Infinite time dependent convection and diffusion on a line. 
ex_convdiff4  One dimensional Burgers equation with steady solution. 
ex_convdiff5  Dominating convection example requiring artificial stabilization. 
ex_convreact1  Time dependent one dimensional convection and reaction model. 
ex_custom_equation1  BlackScholes model equation implemented as a custom equation. 
ex_nonlinear_pde1  1D nonlinear PDE with analytical solution. 
ex_nonlinear_pde2  1D nonlinear PDE with analytical solution extended to 2D. 
ex_diffusion1  Diffusion equation on a unit square with different solutions. 
ex_diffusion2  Diffusion equation on a line with exact solutions. 
ex_heattransfer1  2D heat conduction with natural convection and radiation. 
ex_heattransfer2  One dimensional stationary heat transfer with radiation. 
ex_heattransfer3  One dimensional transient heat conduction. 
ex_heattransfer4  Two dimensional heat transfer with convective cooling. 
ex_heattransfer5  Two dimensional transient cooling shrink fitting example. 
ex_heattransfer6  Axisymmetric steady state heat conduction of a cylinder. 
ex_laplace1  Laplace equation on a unit square. 
ex_laplace2  Laplace equation on a unit circle. 
ex_linearelasticity1  Linear elasticity test case. 
ex_linearelasticity2  Three dimensional example of stress on a bracket. 
ex_linearelasticity3  Parametric study of the bracket deformation model. 
ex_linearelasticity4  Stress of loaded Ibeam supported by two brackets. 
ex_multiphase1  Multiphase flow example of a static bubble. 
ex_multiphase2  Multiphase flow example of a rising bubble. 
ex_multiphase3  Multiphase flow example of a breaking dam. 
ex_heat_exchanger1  Free and forced convection in a heat exchanger. 
ex_natural_convection  Natural convection in a square cavity. 
ex_navierstokes1  Incompressible fluid flow in a channel. 
ex_navierstokes2  Incompressible driven cavity flow. 
ex_navierstokes3  Incompressible fluid flow around a cylinder in a channel (Re=20). 
ex_navierstokes4  Incompressible fluid flow over a backwards facing step. 
ex_navierstokes5  Two dimensional decay of a standing vortex. 
ex_navierstokes6  Time dependent flow around a cylinder in a channel (0<=Re<=100). 
ex_navierstokes7  Laminar flow in a curved three dimensional pipe. 
ex_navierstokes8  Axisymmetric flow in a circular pipe. 
ex_navierstokes9  Axisymmetric impinging jet. 
ex_navierstokes10  3D flow in a pipe using the FeatFlow solver. 
ex_navierstokes11  3D cavity flow in a cube using the FeatFlow solver. 
ex_navierstokes12  3D flow over a backwards facing step using the FeatFlow solver. 
ex_navierstokes13  3D flow past a cylinder using the FeatFlow solver. 
ex_piezoelectric1  Bending of a beam due to piezoelectric effects. 
ex_planestress1  Stress computation for a hole in a thin plate. 
ex_planestress2  NAFEMS plane stress benchmark example. 
ex_planestress3  NAFEMS plane stress benchmark of a tapered plate. 
ex_planestress4  NAFEMS benchmark of plate deformation due to thermal stress. 
ex_planestress5  NAFEMS benchmark for plane stress in an elliptic membrane. 
ex_axistressstrain1  Axisymmetric stressstrain of a hollow cylinder. 
ex_axistressstrain2  Axisymmetric stressstrain of a hollow sphere. 
ex_axistressstrain3  Axisymmetric disc with fixed edge and point load. 
ex_poisson1  Poisson equation on a line. 
ex_poisson2  Poisson equation on a circle. 
ex_poisson3  Poisson equation on a unit square. 
ex_poisson4  Poisson equation on a rectangle with complex solution. 
ex_poisson5  Poisson equation on a sphere. 
ex_poisson6  Poisson equation on a unit cube. 
ex_poisson7  Poisson equation on a unit circle with a point constraint. 
ex_waveequation1  Wave equation on a unit circle. 