ExtendableFEM.jl

This package offers a toolkit to easily setup (mostly low-order, standard and non-standard) finite element methods for multiphysics problems in Julia and to run fixed-point iterations to solve them.

The implementation is based on ExtendableGrids.jl (for meshing and administration) and ExtendableFEMBase.jl (for quadrature and FEM basis functions).

Also note, that this package is part of the meta-package PDELIB.jl

Dependencies on other Julia packages

ExtendableGrids.jl
ExtendableSparse.jl
ExtendableFEMBase.jl
GridVisualize.jl
DocStringExtensions.jl
ForwardDiff.jl
DiffResults.jl

Getting started

The general work-flow is as follows:

1. Geometry description / meshing

The geometry description and meshing is not really separated. For meshes of rectangular domains, there are simple constructors available in ExtendableGrids.jl, e.g.

using ExtendableGrids
## unit square criss-cross into Triangles
xgrid1 = uniform_refine(grid_unitsquare(Triangle2D), 4)

## uniform rectangular grid
xgrid2 = simplexgrid(0:0.1:1, 0:0.2:2)

Note that these grids generate boundary regions from 1 to 4 (bottom, left, top, right) that can be used to assign boundary conditions.

More complex grids can be created via the mesh generators in SimplexGridFactory.jl, see e.g. Examples 245 or 265, or by loading a Gmsh grid file via the corresponding ExtendableGrids.jl extension.

2. Problem description

Before discretizing the user has the option to pose his problems in form of a Problem Description. Note, that usually no grid has to be defined at this point, but region numbers correspond to regions defined in the grid. Here is a short example:

# a simple Poisson problem with right-hand side f(x,y) = x*y and u = 0 along boundary
PD = ProblemDescription()
u = Unknown("u"; name = "potential")
assign_unknown!(PD, u)
assign_operator!(PD, BilinearOperator([grad(u)]; factor = 1e-3))
f! = (result, qpinfo) -> (result[1] = qpinfo.x[1] * qpinfo.x[2])
assign_operator!(PD, LinearOperator(f!, [id(u)]))
assign_operator!(PD, HomogeneousBoundaryData(u; regions = 1:4))

3. Discretization

In this step, the user chooses suitable finite element types for the unknowns of the problem, and generates finite element spaces on the grid (and probably already a solution vector to preoccupy it with an initial solution).

# cubic H1 element in 2D with one component
FES = FESpace{H1Pk{1,2,3}}(xgrid) 
# alternatively: create solution vector and tag blocks with problem unknowns
sol = FEVector(FES; tags = PD.unknowns) 
# fill block corresponding to unknown u with initial values
fill(sol[u], 1)
# interpolate some given function u!(result, qpinfo)
interpolate!(sol[u], u!)

4. Solve

Here, we solve the problem. If the problem is nonlinear, several additional arguments allow to steer the fixed-point iteration, see Stationary Solvers. In the simplest case, the user only needs to call:

# solve problem with finite element space(s)
# (in case of more than one unknown, provide a vector of FESpaces)
sol = solve(PD, FES; init = sol)

For time-dependent problem, the user can add the necessary operators for the time derivative manually. Alternatively, the problem description in space can be turned into an ODE and solve via DifferentialEquations.jl, see Time-dependent Solvers for details.

Also note, that the use can bypass the problem description and control the algebraic level manually via assembling the operators directly into an FEMatrix, see e.g. Example310. It is also possible to take control over the low-level assembly of the operators, see ExtendableFEMBase.jl for details.

5. Plot and postprocess

After solving, the user can postprocess the solution, calculate quantities of interest or plot components.

Gradient-robustness

This package offers some ingredients to build gradient-robust schemes via reconstruction operators or divergence-free elements. Gradient-robustness is a feature of discretisations that exactly balance gradient forces in the momentum balance. In the case of the incompressible Navier–Stokes equations this means that the discrete velocity does not depend on the exact pressure. Divergence-free finite element methods have this property but are usually expensive and difficult to contruct. However, also non-divergence-free classical finite element methods can be made pressure-robust with the help of reconstruction operators applied to testfunctions in certain terms of the momentum balance, see e.g. references [1,2] below.

Recently gradient-robustness was also connected to the design of well-balanced schemes e.g. in the context of (nearly) compressible flows, see e.g. reference [3] below.

References

• [1] "On the divergence constraint in mixed finite element methods for incompressible flows",
  V. John, A. Linke, C. Merdon, M. Neilan and L. Rebholz,
  SIAM Review 59(3) (2017), 492–544,
  >Journal-Link<, >Preprint-Link<
• [2] "Pressure-robustness and discrete Helmholtz projectors in mixed finite element methods for the incompressible Navier–Stokes equations",
  A. Linke and C. Merdon, Computer Methods in Applied Mechanics and Engineering 311 (2016), 304–326,
  >Journal-Link< >Preprint-Link<
• [3] "A gradient-robust well-balanced scheme for the compressible isothermal Stokes problem",
  M. Akbas, T. Gallouet, A. Gassmann, A. Linke and C. Merdon,
  Computer Methods in Applied Mechanics and Engineering 367 (2020),
  >Journal-Link< >Preprint-Link<