Modeling with Non-constant Coefficients in FEATool Multiphysics

Modeling with Non-constant Coefficients in FEATool Multiphysics

FEATool Multiphysics has been designed to be as easy to use as possible, and conveniently features built-in parsing and evaluation of non-constant variables and coefficients. This allows users to quickly and easily enter modeling expressions formulated much as one would write on paper, all without having to develop, write, and compile external user defined functions as other more complex simulation software often requires.

The following applies to all physics coefficients such as point constraints, boundary conditions, equation coefficients, initial conditions and postprocessing expressions, as well as command line functions such as evalexpr, intsubd, intbdr, minmaxsubd, minmaxbdr, and postplot. A coefficient may be composed as a string expression including any number of the various coefficient types listed below, combined with the operators +-*/^, parentheses pairs (), as well as numeric constants and scalars.

Space coordinates

The space coordinates x, y, and z for Cartesian coordinate systems, and r and z for axi-symmetric and cylindrical coordinate systems are valid coefficient expressions. (Note that the space coordinate names can be re-assigned from their default definitions when starting a New Model).

The y-coordinate is for example used to define a parabolic velocity inlet profile as uinlet = 4*0.3/0.41^2*y*(0.41-y) in the flow around a cylinder CFD tutorial example.

Dependent variables

Solution unknowns or dependent variables can also be used in coefficients, with their corresponding names given by the included physics modes, such as for example T for temperature. (As with the space coordinate names, dependent variable names may also be re-defined from the default names when a new physics mode is added. Also note that adding dependent variables in coefficients often make a problem more non-linear and harder to solve.)

The introductory heat transfer multiphysics model for example shows how the temperature T is coupled to the fluid flow via the source term alpha*g*rho*(T-Tc), where alpha, g, rho, and Tc, are model constants, while the flow velocities u and v in turn are coupled back and driving the temperature field through the convective terms.

Space derivatives of dependent variables

By appending a space coordinate postfix to a dependent variable, for example Tx, the derivatives in the corresponding direction can be evaluated (in the example here Tx is equivalent to $\frac{\partial T}{\partial x}$). The capacitance in the spherical capacitor model tutorial is for example calculated by integrating the expression eps0*epsr*(Vr^2+Vz^2)*pi*r, where V is the computed potential.

Second order derivatives, as for example Txy equivalent to $\frac{\partial T}{\partial x\partial y}$, can also be evaluated for dependent variables discretized by the second order Lagrange finite element basis functions.

Note that derivatives of dependent variables are evaluated by direct differentiation of the underlying FEM shape function polynomials and some overall accuracy may be lost in this process. For example, a variable discretized with a first order linear basis function, will feature a piecewise constant first derivative on each element, and zero second derivative (even if a global second derivative would exist analytically). For higher order accuracy with derivatives a projection or gradient reconstruction technique is necessary.


The variable name t is reserved for the current time in instationary simulations (evaluates to 0 for stationary problems). For example, the instationary flow around a cylinder example defines a parabolic and time dependent inlet velocity as uinlet = 6*sin(pi*t/8)*(y*(0.41-y))/0.41^2.


Boundary coefficients may use the external normals by prefixing a space dimension name with the character n, as for example nx for the normal in the x-direction. Boundary normals are evaluated and computed as the outward pointing unit vector from the center of each external cell edge or face.

Discontinuous functions

Logical switch expressions such as T>0 and T<=T0 are also valid syntax and will evaluate to 1 when true and 0 everywhere else. Valid expressions may include the characters <, >, <=, >=, &, and | (corresponding to the less than, greater than, less than or equal, greater than or equal, and, and or operators).

These types of switch expressions can be used to build quite complex expressions and coefficients, for example in the flow over a backwards facing step model the recirculation region is highlighted and visualized with the postprocessing expression x/hstep*(u<0)*(y<0). This expression effectively isolates the region in y < 0 where the x-velocity u is negative, and overlays it with the scaled coordinate x/hstep, resulting in the following plot

FEATool Multiphysics Discontinuous Visualization and Postprocessing

Abruptly varying and discontinuous functions like these can sometimes introduce numerical instabilities if used in equation and boundary conditions. From the perspective of the solver it might then be preferable to use a regularized or smoothed Heaviside function, such as

(abs(w)<1)*(1-(0.5*(1+w+1/pi*sin(pi*w)))) + (w<-1)

or a smoothed Dirac delta function

(w>-1 & w<1)*15/16*(1-2*w^2+w^4)/h_grid

where w is a weighting function defining the transition and width of the smoothed region. For example, with w = (sqrt(x^2+y^2)-0.5)/(2h_grid) a circle with radius 0.5 and transition region width 2h_grid is defined as shown below.

Visualization of Discontinuous Step and Delta Functions

Smoothed step and delta functions are commonly used in models with moving and immersed boundaries which cannot be resolved by the mesh. Step functions are used to define the discontinuous materials, and delta functions to define surface tension forces and reactions. For example see the multiphase flow models ex_multiphase1.

Built-in functions

All of the common built-in mathematical functions and constants in MATLAB, such as pi, eps, sqrt, sin, cos, log, exp, abs etc., may also be used to compose expressions.

Custom and user-defined functions

If a combination of the above methods is insufficient to implement an expression, or the expression can not be described analytically one can use a custom MATLAB m-file function. This could be the case if one wants to use tabulated or experimental data as input which must be interpolated, or complex coefficients with memory effects such as hysteresis. For example, if a coefficient is described as


FEATool will first evaluate the inner function arguments (in this case x, y, t, and T) and then call the function myfun with these arguments. The input arguments can be any combination of the valid coefficients described above. This assumes a function file myfun.m can be found on any of the MATLAB search paths and directories. The function must return a numeric array with the same size and dimensions as the input arguments. A simplified function m-file could for example look like the following

function [ result ] = myfun( x, y, t, T )

persistent data_x data_y data_values T0   % Persistent cache of data
if( isempty(data_x) )   % Load and define on first run.
   load( 'path_to_mydata_file', 'data_x', 'data_y', 'data_values' );
   T0 = 300;

interpolated_data = interp2( data_x, data_y, data_values, x, y );

result = T0 + (1-cos(2*pi*t))*T.^2.*interpolated_data;

index_limit = find( result < 0 );
result(index_limit) = 0;

In a custom function one can then use the usual MATLAB script functions and techniques such as load (to load external data), save, interp, interp2, persistent, varargin, system etc. (To see information and syntax about a MATLAB command one can enter help command in the command line interface.)

Note that depending on the grid size the function may be called several times per evaluation step for different grid cell blocks (for efficiency reasons the default setting is to assembly block sizes of 50000 grid cells or less per function call).

Reserved coefficient names

The coefficient name h_grid is reserved for the mean grid cell size or diameter, and is primarily used in the artificial and convective stabilization techniques.

Internal and local variables may be prefixed with double underscores __varname, and it is therefore not recommended to name variables with underscores to prevent potential name collisions.


  • An easy way to test and check an expression for correctness is to simply plot and visualize it in postprocessing mode.

  • For two-dimensional surface plots in the FEATool GUI one can simply click on a point in the plot to evaluate the visualized expression in the point.

  • When working in the FEATool GUI it is often convenient to use the Model Constants and Expressions dialog box to define and store modeling expressions. The named expressions are then available everywhere in the GUI from equation coefficients to postprocessing expressions.

    FEATool Multiphysics Model
Expressions and Constants Dialog Box

  • The more non-linearities that are introduced the harder it will be for the solver to converge to a solution. Highly non-linear problems typically require more relaxation (which can be achieved by decreasing the Non-linear relaxation parameter in the Solver Settings) and a very good initial guess. One can try to introduce a scalar weighting coefficient to the non-linear coefficients, and step by step increasing it from 0 (no contributions) to 1 (fully non-linear contributions) while using the previous solution as initial guess. Lastly, for stationary problems one can also employ pseudo time stepping with the Backward-Euler scheme to try to use the time dependent solver to reach a steady state.