The FEATool Multiphysics model and simulation problem is defined in a MATLAB struct with fields representing data such as geometry, mesh, equation, boundary conditions, and solution. The toolbox GUI automatically manages this data internally, but it can also be exported from the FEATool toolbox with MATLAB, accessed, and processed with the toolbox API functions at the command line interface CLI (using the Export FEA Struct To MATLAB option from the File menu, option not available as stand-alone toolbox). Alternatively, models can also be saved and exported as MATLAB m-script FEATool models which define the model struct programmatically (select Save As MATLAB Script... from the File menu). The Poisson equation tutorial example shows how to manually define, set up the model definition struct, and solve a typical diffusion problem using the programming API.

Data struct fields

The model definition struct, usually named fea, typically consists of the following data fields

Field	Description	Data type	Required
geom	Geometry definition	struct
grid	Grid/mesh definition	struct	x
sdim	Space dimensions	cell array	x
phys	Physics mode definitions	struct
dvar	Dependent variables	cell array	x
sfun	Shape/FEM basis functions	cell array	x
coef	Coefficients and expressions	struct
eqn	Equation definitions	struct	x
bdr	Boundary condition definitions	struct	x
edg	Edge constraints	struct
pnt	Point constraints	struct
sol	Solution data	struct
vars	Variable data	struct

The fields marked as required are the fields necessary for calling a solver and computing a solution.

Geometry definition

A cell array of geometry objects is contained in the geom.objects field. The geometry is only used by the automatic grid generation process, and can be omitted if a grid is constructed by other means (for example with the grid function primitives or mesh import). The following code for example defines a geometry with a unit circle and a rectangle

fea.geom.objects = { gobj_circle(), geom_rectangle() };

or alternatively a geometry can be imported from a CAD file as

fea.geom = impexp_cad( 'path_to_my_cad_file.stp', 'import' );

Other geometry operations can be performed on the geometry object using the geometry API functions. (A field geom.tags may also be present, but is automatically generated from the top level geometry object tags.)

Grid/mesh definition

The grid field defines the grid and computational mesh on which to perform simulations (corresponding to the geometry if present). The grid data struct itself consists of five sub-fields

grid.p - grid points/vertices
grid.c - grid cell connectivity
grid.a - cell adjacencies
grid.b - boundary data
grid.s - subdomain indices

See the the grid format section for more detailed information on the mesh format.

To create a grid one can use the function grid primitives for creating structured grids, for example circgrid a unit circle

fea.grid = circgrid();

or use the gridgen function to generate a grid from the geometry defined in fea.geom

fea.grid = gridgen( fea, 'hmax', 0.1 );

Furthermore, the grid utility functions can be utilized to manipulate grids and meshes.

Space dimensions

The sdim field contains a cell array of strings with names of the space dimensions. These can then be used in functions such as coefficient and postprocessing expressions (for example x_coord = evalexpr( 'sin(x)', xp, fea )). This field will typically have the following default forms

fea.sdim = { 'x' };             % 1D.
fea.sdim = { 'x', 'y' };        % 2D.
fea.sdim = { 'r', 'z' };        % 2D - Axisymmetry.
fea.sdim = { 'x', 'y', 'z' };   % 3D.

but the strings can be substituted and replaced with other names if desired.

Physics mode definitions

The optional phys struct contains predefined physics modes that can be added with the addphys function, for example

fea = addphys( fea, 'heattransfer' );

fea = parsephys( fea );

A physics mode contains the sub-fields

fea.phys.(tag).name     % Physics mode name (string)
              .tag      % Physics mode tag (string)
              .descr    % Physics mode description (string)
              .dvar     % Physics mode dependent variables (cell array of strings)
              .defsfun  % Default shape functions (cell array of strings)
              .sfun     % Selected shape functions (cell array of strings)
              .eqn      % Physics mode equation definitions (struct)
              .bdr      % Boundary definitions (struct)
              .prop.    % Physics mode properties (struct)
                   .active   % Subdomain active flags (n_dvar x n_subd logical array)
                   .artstab  % Artificial stabilization options (struct)
                   .isaxi    % Axisymmetry flag (logical)
                   .turb     % Turbulence model options (struct)

The parsephys command is applied, after adding physics modes and setting the appropriate physics mode coefficients, to collect and expand all the phys fields to the parent global fea fields dvar, sfun, coef, eqn, and bdr.

Dependent variables

The dvar field contains a cell array of strings with names of the dependent (equation) variables which should be solved for. For example

fea.dvar = { 'u', 'v' };

defines two dependent variables, labeled u and v. This field is automatically generated from the physics modes if present when calling the parsephys function (fea = parsephys(fea)).

Shape/FEM basis functions

Similarly, the sfun field contains a cell array of strings with function names of the finite element shape/basis functions to use for each dependent variable (fea.dvar). Shape functions are defined by the functions in the ellib directory and determine the discretization order and potential accuracy. For example

fea.dvar = { 'sflag1', 'sflag2' };

defines first order conforming basis functions for the first dependent variable, and second order for the second one.

Coefficients and expressions

The coef field is a cell array of constants and coefficient expressions that can be used in boundary, equation definitions, and postprocessing.

It is a cell array of size n_coef x 4, where the first column specifies the coefficient names, second a short description, third a long description, and the fourth column is a sub cell array for the coefficient expression definitions in each subdomain. (If any of the field aliases const, expr, and/or vars exist and are defined they will be merged internally with the coef field.) The assignment below for example

fea.coef = { 'rho', '', '', { 1,  2            } ;
             'miu', '', '', { 3, '4*pi*sin(x)' } };

defines two coefficients rho with values 1 and 2 in two subdomains, and miu with values 3 and the expression 4*pi*sin(x). The descriptions are left blank here as they are only used by the GUI.

Coefficients must be named uniquely from dependent variable and space dimension names. Also internally reserved coefficients are t for simulation time, and h_grid for the mean grid cell diameter (h_grid = cell_volume^(1/n_sdim)).

Problem definition fields

The eqn, bdr, edg, and pnt fields are used to specify equations, boundary conditions, edge and point constraints. The contents and structure of these fields is explained in the following subsections.

Equation definitions

The equation struct eqn contains the following fields to define the global equations for subdomains

Field m	Description
eqn.m.form	Form specifications for temporal terms
eqn.m.coef	Coefficient specifications for temporal terms
	Shape functions are inherited from eqn.a.sfun

Field a	Description
eqn.a.form	Bilinear form specifications
eqn.a.coef	Coefficient specifications for bilinear forms
eqn.a.sfun	Shape function specification for bilinear forms

Field f	Description
eqn.f.form	Linear form specifications for right hand side/load vector
eqn.f.coef	Coefficient specifications for right hand side/load vector
eqn.f.sfun	Shape function specifications for right hand side/load vector

Field dofm/ndof	Description
eqn.dofm	Degree of freedom n numbering map for each cell
eqn.ndof	Numbers of degrees of freedom for each dependent variable

The eqn.m field contains specifications for the time dependent term (with time derivative). Similarly the eqn.a field contains specifications for the bilinear forms used in the iteration (stiffness) matrix, and the eqn.f field specifies the linear forms in the right hand side/load vector.

The above struct fields which contain the form field which specifies the (bi-)linear forms to build and assemble. For bilinear forms the first row corresponds to the trial function space, and the second row the test function space. Linear forms only need to contain one row. In the form specification a 1 indicates a function value, 2 x-derivative, 3 y-derivative, and 4 z-derivative. For example a form specification [2 3;2 3] indicates a bilinear form with two terms, one term with both x-derivatives for the test and trial function spaces, and one with y-derivatives for both spaces (which in this case is a typical two-dimensional diffusion operator). Second order derivatives can also be specified as 22 for the xx-derivative, 23 xy-derivative and so on.

The coef field, is a cell array with coefficient values or expressions used for each term in the form field.

The sfun field is a cell array of shape function names used in the form assembly. This field is usually automatically constructed when calling the parseprob function.

dofm is an array specifying the local to global degree of freedom numbering for each dependent variable (size n_ldof x n_c). The rows correspond to local degrees of freedoms on each cell and the columns give the cell numbers. (For linear conforming shape functions the dof mapping will be identical to the cell connectivity in grid.c.) This field is created when parseprob calls the mapdofbdr function.

ndof is a help array for the numbers of degrees of freedom for each dependent variable (ndof = cellfun(@(dofm) max(dofm(:)), eqn.dofm) and is also automatically generated by mapdofbdr).

Boundary condition definitions

The boundary struct bdr contains the following fields

Field	Description
bdr.d	Dirichlet boundary coefficients
bdr.n	Neumann (flux) boundary coefficients
bdr.bdrm	Boundary degree of freedom numbering maps

Dirichlet boundary conditions are used to prescribe and fix a specific value for the dependent variables, and Neumann conditions are used to represent inward or outward directed fluxes (which are functions of the gradients of the dependent variables). If a Dirichlet condition is prescribed for a boundary, the corresponding Neumann flux entry will be ignored. Alternatively, if a Dirichlet boundary coefficient entry is empty the Neumann flux contribution will be computed and used (however, all default Neumann contributions are zero corresponding to homogeneous do-nothing Neumann boundary conditions).

Boundary coefficients are specified as 1 x n_dvar cell arrays in the bdr.d/n fields, where the entry for each boundary is a n_bc_groups x n_bdr nested cell array containing the coefficients (n_bc_groups is 1, except for special element types such as Hermite basis functions where n_bc_groups > 1 are used to prescribe conditions on degrees of freedom corresponding to derivatives). The coefficient entries can be specified either as constant numeric values, or string expressions which will be evaluated during the simulation. A simplified syntax is also supported where d/n are given as cell arrays of size n_dvar x n_bdr (but does not support boundary groups). The Dirichlet boundary conditions are prescribed to a matrix and right hand side/load vector with the function bdrsetd, while a vector of Neumann flux contributions are computed by the function bdrneu.

In bdrm boundary condition mappings for each dependent variable is specified as a 1 x n_dvar cell array. In a bdrm array entry (size 5 + n_sdim x n_bdof) the first row gives the cell number, followed by edge/face, boundary, global and local degree of freedom numbers, and local coordinates on edges/faces.

The d and n fields must be prescribed by the user or derived from the physics modes with a call to parsephys, the bdrm field is automatically created when parseprob calls the mapdofbdr function. Moreover, if the default solvers solvestat and solvetime detects nonlinear Neumann boundary conditions they will attempt to linearize them by moving terms involving linear forms from the explicit right hand side to the implicit matrix (a boundary contribution to eqn.a instead of eqn.f).

As an example, given a two dimensional time dependent problem with one dependent variable/unknown u, then the standard boundary syntax will look like following

fea.bdr.d = { {  1,  'sin(pi*t*x)+u',  [],        []       } };
fea.bdr.n = { { [],         0,          2,   'nx*ux+ny*uy' } };

this will prescribe u constant u = 1 on the first boundary, the time dependent expression u(x,t) = sin(pi*t*x)+u on the second boundary, a constant Neumann/flux condition 2 on the third boundary, and lastly the flux expression nx*ux+ny*uy on the fourth and last boundary. Note that in string expressions it is perfectly valid to use constants, MATLAB functions, time t, function values and derivatives of dependent variables, as well as boundary normals (nx, ny, and nz).

Edge constraints

The edge struct is only applicable to 3D problems and contains the following fields to specify optional Dirichlet edge constraints

Field	Description
edg.index	Index to edge
edg.type	Must be constraint type
edg.dvar	Integer or string corresponding a dependent variable
edg.expr	scalar or string expression for constraint

Dirichlet edge constraints are specified in a edg struct. The type field must be set to the string constraint. Edge constraints are applied to both the load vector f and the global matrix A with the information in the finite element problem struct (edgset). The index field specifies the edge to apply the constraint to (as reconstructed by gridbdre and can be visualized with plotedg). The dvar field contains an integer or string expression pointing to a dependent variable in dvar. The constraints are specified in the expr field either as a string expression or numeric scalar.

Point constraints

The optional point struct contains the following fields

Field	Description
pnt.index	Index to grid point (in grid.p)
pnt.type	Specifies either sources or constraint
pnt.dvar	Integer or string corresponding a dependent variable
pnt.expr	scalar or string expression for constraint

Point sources and Dirichlet point constraints are specified in a pnt struct. The type field specifies either a point source or constraint for each point with a corresponding string value. Point sources are applied to the right hand side load vector f while constraints are applied to both f and the global matrix A with the information in the finite element problem struct (pntsetf and pntset). The sources or constraints will be applied to the degree of freedom closest to the grid point specified in the index field. The dvar field contains an integer or string expression pointing to a dependent variable in dvar. The sources or constraints are specified in the expr field either as a string expression or scalar so that pnt.dvar(grid.p(:,pnt.index)) = pnt.expr. (Note that a point constraint will override the usual Dirichlet boundary conditions at the given point if present.)

Solution data

After calling a solver (for example using any of the following solver calls below) and computing a solution to a problem (see the help/docstring for the corresponding solvers for description of the input arguments and output)

fea.sol.u = solvestat( fea );
[fea.sol.u,fea.sol.t] = solvetime( fea );
[fea.sol.u,fea.sol.l] = solveeig( fea );
[fea.sol.u,fea.sol.t] = fsisolve( fea );

fea = fenics( fea );
[fea.sol.u,fea.sol.t,fea.vars] = openfoam( fea );
[fea.sol.u,fea.sol.t,fea.vars] = su2( fea );

the fea.sol field will contain the solution vector(s) in sol.u, and additionally the output times/eigenvalues will be stored as a vector in either the sol.t or sol.l field, respectively.

The solution vector u is a numeric array with rows corresponding to the solution unknowns or degrees of freedom (dof) ordered sequentially according to the dependent variables listed in fea.dvar as shown below

sol.u(:,i_time) = [ dof   1 for fea.dvar{1}, ...
                    ...
                    dof m_1 for fea.dvar{1}, ...
                    ...

                    dof   1 for fea.dvar{n}, ...
                    ...
                    dof m_n for fea.dvar{n} ];

And for time dependent and eigenvalue problems the columns in u correspond to solutions at different times/eigenvalues. In cases were a physics mode is inactive in some subdomains, the corresponding solution values will be set to nan (not a number, which are not plotted visible during postprocessing).

After physics and equation parsing (parsephys/parseprob) the degrees of freedom for a given grid/mesh will be specified by the dof mapping in fea.eqn.dofm, and the number of dofs per dependent variable in fea.eqn.ndof (as computed with the function mapdofbdr). The function evalinit can be used to compute the actual coordinates of the degrees of freedom with the following m-code script

n_dvar = length(fea.dvar);
x0 = [];
for i=1:length(fea.sdim)
  init_expressions = repmat(fea.sdim(i), 1, n_dvar);
  x0 = [x0, evalinit(fea, init_expressions)];
end

Variable data

An optional field of variable data vars can also be present and have the sub-fields

name   % Variable name
type   % Variable type ('expr' or 'var')
descr  % Variable description (optional)
sfun   % Variable shape definition (default 'sflag1')
data   % Variable data

If the variable is of type expr the variable will be merged with all other coefficients and expressions in fea.coef (as fea.coef = [ fea.coef; { vars(i).name, vars(i).data } ], see also Coefficients and expressions).

However, if the type is var, then the variable will be treated similar to a dependent variable with shape function sfun (default sflag1 if unspecified) in equation and postprocessing expressions. In this case the data field should correspond to the same format and shape as an equivalent dependent variable (for example if vars(i).sfun = 'sflag1' then vars(i).data should have the shape number of grid points x 1 with values in fea.grid.p).