ChasteGuides/FiniteElementAssemblersAndSolvers

FEM assemblers and solvers

Note 1: New users should read the tutorial on solving linear PDEs before reading this; this page is more on how the assemblers and solvers work internally.

Note 2: this page only refers to the main assemblers and solvers, which use linear basis functions. The elasticity assemblers use quadratic basis functions and are outside this hierarchy.

For an example of how to convert an assembler from old to new design in the simplest way, see ConvertingAssemblersToUseNewDesign.

The main merge of the new assemblers and solvers occurred in r9856 (with final tidying, renames in r9858, r9863).

  1. FEM assemblers and solvers
  2. Quick Summary
  3. Full Description
    1. AbstractFeObjectAssembler
      1. Introduction
      2. ComputeMatrixTerm(), ComputeVectorTerm(), ComputeVectorSurfaceTerm()
      3. InterpolationLevel
    2. The solvers of linear PDEs
      1. Introduction
      2. AbstractAssemblerSolverHybrid
    3. Concrete solvers
      1. Non-heart
      2. Cardiac assemblers and solvers

Quick Summary

The old design involved classes that were both assemblers (of finite element matrices/vectors), and solvers. Assemblers and solvers are now completely seperate. See also the diagram at the bottom of this page.

Assemblers

The class AbstractFeObjectAssembler is a generic assembler class (linear bases only) and provides the basic functionality for assembling finite element matrices (for example, a stiffness matrix, or mass matrix, or the matrix A in Ax=b), or vectors (eg b in Ax=b). Concrete classes which inherit from AbstractFeObjectAssembler say whether they want to assemble a matrix or vector (or both); and if they are assembling a matrix, implement the pure-like method ComputeMatrixTerm() which says exactly what the matrix to be assembled should be. Similarly concrete assemblers that assemble vectors implement ComputeVectorTerm and maybe also ComputeVectorSurfaceTerm if boundary integral contributions need to be added to the vector.

Solvers

The main solvers now only solve linear PDEs. There is an abstract static an abstract dynamic solver, defining the following interfaces:

AbstractStaticLinearPdeSolver has

AbstractDynamicLinearPdeSolver has

Concrete classes need to implement the pure method SetupLinearSystem(), which should completely setup the linear system to be solved (each timestep). How they set it up will depend on both the PDE solved and the precise numerical scheme used. They can now make use of all the assemblers they want.

Hybrid

The class AbstractAssemblerSolverHybrid is a useful class for allowing concrete classes of the old design: classes that are both assemblers and solvers. AbstractAssemblerSolverHybrid inherits from AbstractFeObjectAssembler, and has a method for constructing a linear system using its parent assembler. Concrete classes should inherit from it and the static/dynamic solver (as appropriate). See below for more details.

Cardiac assemblers and solvers

See section in full description (more or less self-contained).

Full Description

Finite element solution of PDEs involves solving a linear system, whether it is (i) the single linear system solved in the FEM discretisation of a static linear problem; or (ii) the linear system solved each timestep of a FEM discretisation of a time-dependent linear problem; or (iii) the linear system solved each newton iteration of a nonlinear problem.

The previous assembler hierachy was actually an assembler-solver hierachy, and based on the notion that every problem involved defining how to do the finite element assembly of the left-hand-side (LHS) matrix and how to do the finite element assembly of the right-hand side (RHS) vector, ie the matrix A and vector b in the linear system Ax=b. Here, finite element assembly means looping over elements, doing numerical integration on each element, and adding small matrices/vectors to A or b. This led to the solver-assembler hierachy, where the concrete solver-assemblers were partly a PDE solver interface and partly an assembler of LHS matrix A and RHS vector b.

This worked fine for simple FEM implementations. However, for more advanced implementations, such as monodomain problems with matrix-based-RHS, where the linear system is of the form Ax=Mz, where A and M (mass matrix) need to be assembled in a finite element manner, but not z, this doesn't fit in very well. An even worse example is monodomain with a 'correction' term, where the linear system is Ax=Mz+c, where A,M and c need to be assembled in a finite element manner.

This led to the new assemblers and solvers, which is no longer a complex assembler-solver hierachy, but a simple assembler hierarchy, and a separate solver hierarchy.

AbstractFeObjectAssembler

Introduction

This class is used for creating any matrix or vector which needs to be assembled in a finite element manner, ie, the matrix A or vector b in Ax=b; or A or M or c in Ax=Mz+c. It currently uses linear basis functions only.

The class is templated over element, spatial and problem dimensions as usual, but also on two booleans, CAN_ASSEMBLE_VECTOR and CAN_ASSEMBLE_MATRIX: 

template<unsigned ELEM_DIM,unsigned SPACE_DIM,unsigned PROBLEM_DIM,bool CAN_ASSEMBLE_VECTOR, bool CAN_ASSEMBLE_MATRIX, InterpolationLevel INTERPOLATION_LEVEL>
class AbstractFeObjectAssembler

The methods which need to be implemented by the concrete classes are: 

ComputeMatrixTerm();
ComputeVectorTerm();
ComputeVectorSurfaceTerm();

These methods are essentially the pure methods which need to be implemented by the concrete class. Technically, though, they are not pure: they have a default implementation (which is to throw an error), so that the concrete classes just need to implement the ones that will be needed.

For example: To create a concrete class which is to be used for assembling the mass matrix M in 1d, the user has to write a concrete class inheriting from AbstractFeObjectAssembler<1,1,1,false,true,NORMAL>, and then implement the ComputeMatrixTerm() methods (the other two methods don't need overloading as no vectors are assembled). Similarly, to create a concrete class for assembling a vector, we inherit from AbstractFeObjectAssembler<1,1,1,true,false,NORMAL> and implement ComputeVectorTerm(), and maybe ComputeVectorSurfaceTerm() if boundary integrals are also added to the vector. The calling code then needs to provide a constructed Mat or Vec, and the assembler will assemble them. See pde/test/TestAbstractFeObjectAssembler.hpp for examples of concrete classes and usage.

ComputeMatrixTerm(), ComputeVectorTerm(), ComputeVectorSurfaceTerm()

These are the methods the concrete assembler has to implement. They need to provide the integrand of the intregral being evaluated at the current quadrature point, which will be a vector/matrix of size num_nodes_in_element

For example: Suppose we are in 2d, and want to set up the vector b, where 

b_i = \integral_{\Omega} f(x) \phi_i(x) dV

which is the RHS vector needed to solve linear elliptic equations with source term f.

We inherit from AbstractFeObjectAssembler<1,1,1,true,false,NORMAL> and implement ComputeVectorTerm(). This has parameters which include the quadrature point rX, and the vector of basis functions evaluated at that quad point rPhi ( = [\phi_1(x), \phi_2(x), \phi_3(x)]): 

c_vector<double,PROBLEM_DIM*(ELEMENT_DIM+1)> ComputeVectorTerm(
                c_vector<double, ELEMENT_DIM+1>& rPhi,
                c_matrix<double, SPACE_DIM, ELEMENT_DIM+1>& rGradPhi,
                ChastePoint<SPACE_DIM>& rX,
                c_vector<double,PROBLEM_DIM>& rU,
                c_matrix<double, PROBLEM_DIM, SPACE_DIM>& rGradU,
                Element<ELEMENT_DIM,SPACE_DIM>* pElement)

This implementation needs to return the vector [f(x)phi_1(x), f(x)phi_2(x), f(x) phi_3(x)]. It will therefore use rX and rPhi from the parameters, but not any of the others.

ComputeMatrixTerm() and ComputeVectorSurfaceTerm() are similar.

InterpolationLevel

In the above example the parameters rU, rGradPhi, rGradU are not needed by the concrete class, and therefore didn't need interpolating onto the current quadrature points. Some control is available over which quantities are interpolated through the final template parameter of the abstract class 

template<unsigned ELEM_DIM,unsigned SPACE_DIM,unsigned PROBLEM_DIM,bool CAN_ASSEMBLE_VECTOR, bool CAN_ASSEMBLE_MATRIX, InterpolationLevel INTERPOLATION_LEVEL>
class AbstractFeObjectAssembler

For example if InterpolationLevel is not equal to NONLINEAR, rGradPhi is not computed for vectors. See class documentation for available options. Non-interpolated parameters are just passed in which zeroed values, so choosing an incorrect InterpolationLevel will lead to the solver silently giving completely incorrect answers if a quantity is used when it hasn't been interpolated.

The solvers of linear PDEs

Introduction

There are three relatively simple solver classes for linear PDEs: AbstractStaticLinearPdeSolver for solving linear time-independent problems, AbstractDynamicLinearPdeSolver for solving linear time-dependent problems, and a parent which includes some common functionality 

                    AbstractLinearPdeSolver 
                        ^            ^
                        |            |
                        |            |
AbstractStaticLinearPdeSolver     AbstractDynamicLinearPdeSolver

AbstractStaticLinearPdeSolver has

AbstractDynamicLinearPdeSolver has

The three classes are abstract because of the following pure method in AbstractLinearPdeSolver

This is the method that needs to be implemented by a concrete solver. It is called by the two Solve methods, and should completely set up the linear system to be solved (that timestep), using whichever assemblers it needs and applying boundary conditions if necessary.

AbstractAssemblerSolverHybrid

Now, in principle, to solve a particular PDE, we need to create a concrete assembler, MyAssembler, for that PDE, and then create a concrete solver, MySolver, which implements SetupLinearSystem using MyAssembler. This would be a "has-a" solver-assembler relationship ("solver has an assembler"), and is the new design.

For convenience, the old "is-a" relationship ("solver is also an assembler") is still allowed, through the class AbstractAssemblerSolverHybrid. This can be used when an assembler will be written which will define the whole of the LHS matrix A, and the whole of the RHS vector b.

AbstractAssemblerSolverHybrid inherits from 

AbstractFeObjectAssembler<..,true,true,..>

and has a method SetupGivenLinearSystem(..,..,LinearSystem*) which sets up the given linear system using the parent-assembler part of the class.

Concrete classes should inherit from this and one of AbstractStaticLinearPdeSolver or AbstractDynamicLinearPdeSolver. They would then have the following pure methods

The concrete class then has to implement the Compute*Term methods as appropriate for its particular PDE. For SetupLinearSystem it simply has to have a one-line function which just calls SetupGivenLinearSystem. See SimpleLinearEllipticSolver as an example.

Concrete solvers

Non-heart

Cardiac assemblers and solvers

There are two types of cardiac solver, one which assembles the RHS vector (Ax=b, b assembled each timestep; slow), and one which uses matrix-based setup of the RHS vector (Ax=Mz, A,M assembled once; fast). The latter obviously must use a "has-a" solver-assembler approach (solver has two assemblers, in fact). The former could use an "is-a" approach, ie inherit from AbstractAssemblerSolverHybrid; however there are reasons to also use a "has-a" approach.

See also diagram below.

Monodomain

There is one main assembler class

There are two main solvers, inheriting from an abstract solver

Bidomain

Bidomain is a bit more complicated, due to the possibility of a bath and extra bidomain specific functions. There are again two main solvers, BasicBidomainSolver and MatrixBasedBidomainSolver, both of which can deal with both bath/non-bath problems. See class documentation/code for more details.

Assemblers:

Solvers: 

                                                AbstractLinearPdeSolver /* pure SetUpLinearSystem() method */
                                                    ^            ^
                                                    |            |
                                                    |            |
                            AbstractStaticLinearPdeSolver     AbstractDynamicLinearPdeSolver /* SetTimes(), SetInitCondition(), Solve() methods */
                                    ^  /* Solve() method */      ^                 ^     ^
                                    |                            |                 |     |
 SimpleLinearEllipticSolver ________|                            |                 |     |__________________ SimpleLinearParabolicSolver
   /* is a assembler, ie also */                                 |                 |                           /* is a assembler, ie also */ 
   /* inherits from hybrid    */                                 |                 |                           /* inherits from hybrid    */ 
                                                                 |                 |
                                           AbstractMonodomainSolver            AbstractBidomainSolver
                                              ^              ^                     ^             ^
                                              |              |                     |             |
                                              |              |                     |             |
                                BasicMonodomain     MatrixBasedMonodomain   BasicBidomain     MatrixBasedBidomain /* all implement SetUpLinearSystem() */
                                //uses 1 assembler     uses 2 asmblrs       uses 1 asmblr       uses 2 asmblrs