# Python FEM and Multiphysics Simulations with FEniCS and FEATool

FEniCS is a flexible and comprehensive finite element FEM and partial differential equation PDE modeling and simulation toolkit with Python and C++ interfaces along with many integrated solvers. As both FEATool and FEniCS discretize equations employing a weak finite element formulation it is quite straightforward to translate FEATool syntax and convert it to FEniCS python scripts. Similar to what has been done with the external CFD solvers, this post introduces the FEATool-FEniCS external solver integration allowing for easy conversion, exporting, solving, and importing FEATool Multiphysics models to FEniCS directly from the GUI, as well as the MATLAB command line interfaces.

In contrast to specialized solvers, FEniCS is aimed at
supporting and solving general systems of PDEs, such as found in
engineering and coupled multiphysics problems. As both FEATool and
FEniCS discretize the equations in the same way the solutions they
produce should be virtually identical. One of the many advantages of
using FEniCS then, in addition to now supporting the Python scripting
language as well as MATLAB, is that FEniCS solvers have
support for both distributed (MPI) and shared memory (OpenMP) parallel HPC
execution allowing for larger models to be solved significantly
faster (FEniCS has been tested on problem sizes up to 10^{8}
degrees of freedom run on 512 CPUs in parallel, where MATLAB and
Octave in contrast is limited to the serial sparse linear solvers, the
Umfpack direct solver and selection of iterative ones).

## GUI Usage

A button labeled **FEniCS** can be found in the *Solve Mode* toolbar
of the FEATool GUI. Pressing this button opens a dialog box where the
current model equations and parameters have been translated and
converted to a FEniCS python script. The script can be inspected and
edited, as well as changing the output file name and system shell
command to run it (which is *bash* by default).

The FEniCS solver dialog box also features **Solve**, **Auto**,
**Export**, and **Import** buttons. If the **Auto** toggle button is
engaged the **Solve** button will automatically try to **Export**,
**Solve**, and **Import** the computed solution. With the **Auto**
toggle disengaged the **Solve** button only execute the FEniCS python
script in a bash shell if present (not export or import). The
**Export** option saves the current grid in Dolfin XML format and
exports the FEniCS simulation script and **Import** will import an
existing solution if it matches with the FEATool problem definition.
Additional options include the base *File Name* and overriding the
default bash *Solve Command*.

## Command Line Usage

The installed `fenics`

MATLAB
function can also be used on the command line (CLI) to manually
perform the export, solve, and import actions. In the following is an
example to set up a simple heat transfer model on a unit circle with a
unit heat source term, and *T=0* fixed temperature on all the
boundaries. In the FEATool MATLAB m-script language the model will
look like the following

```
fea.sdim = {'x' 'y'};
fea.grid = quad2tri(circgrid());
fea = addphys(fea, @heattransfer);
fea.phys.ht.eqn.coef{6,end} = {1}; % Heat source term.
fea.phys.ht.bdr.sel(:) = 1; % Zero temperature boundary conditions.
fea = parsephys(fea); % Parse physics mode and fea problem struct.
fea = parseprob(fea);
fea.sol.u = solvestat(fea);
postplot(fea, 'surfexpr', 'T')
```

By running the command `fenics(fea, ‘mode’, ‘export’)`

a
grid file *featool_fenics_mesh.xml* and FEniCS simulation script file
*featool_fenics.py* will automatically be generated in the current
directory

```
from fenics import *
# Mesh and subdomains.
mesh = Mesh("featool_fenics_mesh.xml")
# Finite element and function spaces.
E0 = FiniteElement("P", mesh.ufl_cell(), 1)
V = FunctionSpace(mesh, E0)
T = Function(V)
T_t = TestFunction(V)
# Model constants and expressions.
rho = Constant(1)
cp = Constant(1)
k = Constant(1)
u = Constant(0)
v = Constant(0)
q = Constant(1)
# Bilinear forms.
a = ( (rho*cp*u)*T*T_t.dx(0) + (k)*T.dx(0)*T_t.dx(0) + \
(rho*cp*v)*T*T_t.dx(1) + (k)*T.dx(1)*T_t.dx(1) )*dx
# Linear forms.
f = q*T_t*dx
# Boundary conditions.
dbc0 = DirichletBC(V, Constant(0), 0)
dbc1 = DirichletBC(V, Constant(0), 1)
dbc2 = DirichletBC(V, Constant(0), 2)
dbc3 = DirichletBC(V, Constant(0), 3)
dbc = [dbc0, dbc1, dbc2, dbc3]
# Initial conditions.
assign(T, interpolate( Constant(0), V))
# Solve.
solve( a - f == 0, T, dbc )
# Output.
import numpy as np
T_h = T.compute_vertex_values(mesh)
np.savetxt("featool_fenics_sol.txt", np.column_stack((T_h)))
# Postprocessing.
plot(T, title = "T")
interactive()
```

The generated FEniCS python FEM script is longer and somewhat more verbose since all the FEATool physics mode defaults must be explicitly expressed.

Calling `fenics(fea, ‘mode’, ‘solve’)`

will attempt to
solve the problem if python and FEniCS is installed and set up
correctly. Alternatively, the script can be run by itself by a valid
FEniCS installation.

The solution process will generate the file *featool_fenics_sol.txt*
containing the nodal solution(s) which can be imported back into
FEATool with `fea = fenics(fea, ‘mode’, ‘import’)`

after
which it can be postprocessed and visualized with the usual FEATool
and MATLAB functions.

## Notes

The FEniCS-FEATool integration and problem file export should work for general multiphysics problems using both the GUI and command line. However, a subset of FEATool functionality is currently not supported (the known ones are listed here)

Non-linear Dirichlet fixed value boundary conditions are not supported. The exported FEniCS C++ boundary condition expressions do not currently allow for dependent variables / solution unknowns.

As FEniCS does not feature a built-in time dependent solver. Solvers for instationary problems must currently be implemented manually.

Switch and logical expressions (for example

*x>1 & y<=0*) must be manually rewritten to conform with the C++ expression syntax.The ordering of expressions and constants is arbitrary in FEATool while in FEniCS constants and expressions must be defined sequentially in the order they are used, and may therefore need to be rearranged in the auto generated FEniCS problem file.

Point source terms are not supported by the FEniCS non-linear solution form (which is used by default since it also handles linear problems).

Although FEATool and FEniCS uses the same finite element FEM basis functions the internal ordering is different, thus the FEniCS solution will be exported in the grid points (

*P1*space). Thus even if higher-order FEM spaces are used solution accuracy will be lost during FEATool import. This issue might not be that important for visualization but when calculating quantities, boundary and subdomain integrals, and expression evaluations it would currently be more accurate to perform that in FEniCS and python.FEniCS does not currently support quadrilateral or hexahedral grid cells and must be converted to triangles and tetrahedra, respectively. Grid conversion of these types of grids can be performed directly in the GUI or with the

`quad2tri`

and`hex2tet`

commands. Typically this will not be an issue since FEATool and the automatic mesh generator`gridgen`

by default creates simplex grids (line segments, triangles, and tetrahedra).In addition to supporting FEniCS, the exported python simulation scripts should also be compatible with the Firedrake project solver which also uses the FEniCS Unified Form Language (UFL) for problem definitions.

## Installation

FEATool-FEniCS functionality is integrated with the default FEATool Multiphysics distribution. However, FEATool does not include a Python interpreter or FEniCS itself which must be installed separately. The FEniCS homepage provides instructions how to install FEniCS on Linux systems and pre-configured Docker Linux images.

For systems running *Windows 10* FEniCS can be installed with the Ubuntu Bash Windows Subsystem for Linux by simply
opening a Windows Bash shell and running the FEniCS on Ubuntu commands
(which automatically also installs Python if required)

```
sudo add-apt-repository ppa:fenics-packages/fenics
sudo apt-get update
sudo apt-get install --no-install-recommends fenics
sudo apt-get dist-upgrade
```

To allow FEniCS and python plotting and visualization on Windows one must also install an X window server such as Xming or VcXsrv (note that plotting of the solution variables is commented and disabled in the FEATool-FEniCS scripts by default).