Mix formulation and integration on skeleton (h-adaptivity)

Table of Contents

Set by step example how to implement mix finite element and do hp-adaptivity

Here we present general procedure applied to h-adaptivity of mix problem. MoFEM is generic finite element code capable of solving problems from magnetics, mechanics of fluids and solids, and as a result of this operates on necessary complex data structures able to carry abstract problems. If you master this example, applying the same approach, you will be able to implement in similar way other problems, for instance, Discontinuous Galerkin Method, or thermo-mechanical coupled problem.

In this tutorial, we show how to implement h-adaptivity for mix transport formulation for stationary transport/heat-conduction problem (The Poisson equation). We focus on h-adaptivity, however, this implementation allows running problem with the arbitrary order of approximation and extension to p-/hp-adativity is possible.

Code relevant to this tutorial could be found in Mix transport element module. Plane source code for main program is available under link h_adaptive_transport.cpp, whereas main class implementation of finite element and finite element operators is available here MixTransportElement.hpp.

MoFEM is built on three pillars, PETSc library which manages abstraction related to the algebraic system of equations and delivers interface to several solvers, mesh-oriented database (MoAB) which manages complexities related to mesh topology and delivers various format readers and Boost library which through multi-indices manages access to containers with MoFEM data structures.

MoFEM handles complexities related to finite element method, degrees of freedom, approximation spaces and base functions. It provides tools and structures for implementation classical and nonclassical finite element codes and can handle linear, nonlinear and time dependent problems.

Here we are considering a problem

\[ A_{ij} \sigma_j + u_{,i} = 0 \; \textrm{on} \; \Omega \\ \sigma_{i,i} = f \; \textrm{on} \; \Omega \]

where first equation is a physical (constitutive) equation and tensor \(A_{ij}\) express material properties, for example, if \(A_{ij}\) is a diagonal matrix with single constant on diagonal we have an isotropic material. The second equation is conservation equation, which expresses law that mass or energy can not be created or destroyed (close system).

Multiplying the constitutive equation by a test field \(\tau\) and integrating by parts over \(\Omega\) (taking into account the inhomogeneous Dirichlet boundary condition), we obtain

\[ \int_\Omega A_{ij}\sigma_j \tau_i \textrm{d}x - \int_\Omega u \tau_{i,i} \textrm{d}x = \int_{\Gamma_\textrm{u}} \overline{u} \tau_{i} n_i \textrm{d}\Gamma ,\quad \forall \boldsymbol\tau \in H(\textrm{div},\Omega) \]

while from the equilibrium equation by multiplying by test field \(v\) we obtain

\[ \int_\Omega \sigma_{i,i}v \textrm{d}x = \int_\Omega f v \textrm{d}x, \quad \forall v \in L^2(\Omega). \]

Before we look into implementation details, you can run simple example

cd basic_finite_elements/mix_transport
mpirun -np 2 ./transport \
-my_file l-shape.msh -meshsets_config bc.cfg \
-ksp_type fgmres -pc_type lu -pc_factor_mat_solver_package mumps -ksp_monitor \
-my_order 1 -nb_levels 5 2>&1 | tee log_order1


Above code produces output files, out_0.h5m, out_1.h5m, ... where number of files is dependent on number of refinement levels. You can convert those files to VTK format running script

../nonlinear_elasticity/do_vtk.sh out_*.h5m

and post-process results in ParaView http://www.paraview.org.

As result of post-processing in ParaView we have created following output

Flux singularity
Convergence plot

You can note that subsequent mesh refinement reduces error and make it more uniformly distributed through the mesh. You can also observe that mesh is getting denser at the corner where singularity is located.

Should be implemented and tested problem from this article [18]. In this paper analytical solution is available, allowing to tested efficiency of evaluated error. You are very welcome to improve that documentation.

Explaining mesh BitRefLevel

MoFEM::BitRefLevel is a generic tool which allows working on several overlapping meshes, partially sharing entities and topology. In the context of this problem, we would use it to create and work with a hierarchy of mesh refinement.

For each entity in a database (vertex, edge, triangle, quad, ..., tetrahedra, brick, prism, ...) we attach bit tag, i.e. MoFEM::BitRefLevel. Using that tag we can make selections of entities and other operations on meshes.

It will be easier to understand if you look at the figure below

Mesh bit refinement levels (Example)

Imagine that we load the mesh with the topology like in the top-left figure. We seed the mesh (MoFEM::Interface::seed_ref_level) by tagging each entity with the first bit on, i.e. 001.

Next, we refine by splitting one of the edges, see middle-left figure. You can understand that as edge-splitting-mesh refinement, see [36]. As the result of this mesh refinement new vertex in the middle of the edge is created, four new edges and four triangles. The remaining of mesh entities is part of old and new mesh level. As a consequence entities which are shared between the second mesh and first mesh have two bites on, i.e. 011, since those are coexisting on two meshes. New entities have the only second bit on, i.e. 010.

In the next step, we made third mesh refinement. The procedure repeats, we have some new entities which have the only third bit on, i.e. 100, some entities are shared between refinement second and third, those have third and a second bit on, i.e. 110, and some entities are shared between all refinement levels, i.e. 111. In current implementation, we assume that user would not have more than 32 refinement levels. However bit squashing and bit level deletion allows form deep (more than 32) mesh refinements.

Physically all three meshes overlap, we do not delete overlapped entities or create new whole mesh each time. It is the reason behind that since we like to transfer data between different refinement levels without the redundancy and unnecessary copying and communication between processors if the mesh is distributed.

Now we present examples of mesh selections; you can practically do everything that you do with bits, see https://en.wikipedia.org/wiki/Bitwise_operation#Bit_shifts. For example;

That allows to project data between meshes, do calculations on the part of the mesh, to introduce cracks, make topological changes like face flipping. Note that in problem below is defined for particular bit levels. See link Get entities and adjacencies for details how to get entities and adjacent by bit refinement levels.

Setting up problem

We consider the problem which potentially can have complex geometry and a non-trivial set of boundary conditions. Here we use Gmsh (see http://gmsh.info) to generate the mesh. In this section, we focus attention on how to read mesh file and how to apply some arbitrary boundary conditions on that mesh.

The first step is to create mesh database, where information about the problem is stored. We create MoAB (see http://ftp.mcs.anl.gov/pub/fathom/moab-docs/) instance and interface to it

moab::Core mb_instance; ///< Create database
moab::Interface& moab = mb_instance; ///< Create interface to database

Once we have created mesh database, we get from command line name of file which we like to load

PetscBool flg = PETSC_TRUE; ///< Is true of option is set in line command
char mesh_file_name[255]; ///< File name
ierr = PetscOptionsGetString(PETSC_NULL,PETSC_NULL,"-my_file",mesh_file_name,255,&flg); CHKERRQ(ierr);

where here we are using PETSc function to do that, see http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Sys/PetscOptionsGetString.html for more details. Having mesh file name, we can load mesh, created in Gmsh, to database

const char *option;
option = "";
rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);

MoFEM can read various file format in this example we are using mesh from Gmsh, where example file geometry is stored in file l-shape.geo and mesh is in l-shape.msh. For more details about file formats please look into MoAB documentation.

Creating MoFEM database to store information about fields, finite elements and problem, we link it to MOAB database and create interface to manage data

MoFEM::Core core(moab); ///< Create database
MoFEM::Interface& m_field = core; ///< Create interface to database

In file l-shape.msh following geometry and mesh is described

Blocksets and GMesh mesh

where as result of reading of that mesh into MoFEM database, following information is printed on the screen

read cubit meshset 12682136550675316737 type BLOCKSET msId 2
read cubit meshset 12682136550675316738 type BLOCKSET msId 3
read cubit meshset 12682136550675316741 type BLOCKSET msId 1

The output to the screen indicates that three BLOCKSETs are available. Those BOLCKSETs are interpreted as Physical Volume and Physical Surfaces which were created in Gmsh. User can create more sophisticated geometry or use more complex set of boundary conditions. Interpretation of BLOCKSETs is a matter of separate config file, which describes problem specific types of boundary conditions.

The interpretation of BLOCKSETs is done by the meshset config file bc.cfg,

temperature_flag1=1 # 0: N/A, 1: temperature value applied
heatflux_flag1=1 # 0: N/A, 1: temperature value applied
name=MAT_THERMAL # Hast to be set for Thermal Mat

According to that file, BLOCKSET 3 is set to be NODESET 1001 where temperature is applied \(u=0\, x \in \partial \Omega_u\) and BLOCKSET 2 is set to be SIDESET 1002 where flux is applied \(n_i\sigma_i = 1\, x \in \partial \Omega_\sigma\). On BLOCKSET 1 is the domain \(\Omega\) where problem is defined.

In the above example, temperature type, heat flux type and temperature volume is used. This is a matter of interpretation as well, we are using this convention to be consistent with other meshing programs like CUBIT where boundary conditions like that can be set directly. However, MoFEM user can set own types and attach own interpretation if needed.

Boundary condition meshsets are managed in MoFEM by MoFEM::MeshsetsManager interface. Using that interface user can access information about boundary conditions and add a new one. All modules and software components using that interface can easily share information through it. To get access to interface, user needs to query MoFEM for it, as follows

MoFEM::MeshsetsManager *meshsets_manager_ptr; //< Click on class to see more details
ierr = m_field.getInterface(meshsets_manager_ptr); CHKERRQ(ierr);

using interface we can read file bc.cfg, parse it and add boundary conditions. All this is automatized for user convenience and done with single interface method

ierr = meshsets_manager_ptr->setMeshsetFromFile(); CHKERRQ(ierr); //<- Click on function to see more details

In general, user can query about range of interfaces giving access to specific tools. In this example we are using two types of interfaces, one to manage BLOCKSETs, i.e. MoFEM::MeshsetsManager and another for refining meshes, i.e. MoFEM::MeshRefinement.

MoFEM uses interfaces to do particular tasks, could have one or more interface for the same task, which differs in complexity, one is simple to do basic operation another can be for advanced users. It can have different versions of the interface. We design MoFEM like that to have a structure which can adapt in time to new technologies and programmers need.
Note that all commands return an error code and after each function, an error code is checked, in a case of the error, appropriate information is printed during execution time. It is important to follow that rule, of checking the error code, it simplifies code debugging and bug reporting.

Solving problem and calculating matrices and vectors

We start with calculating initial problem on mesh created with Gmsh. The problem is broken down into several stages, implementation of a mix transport problem is in ExampleMix which is derived from MixTransport::MixTransportElement.

We have the hierarchy of classes (structures), some of them are generic thus are part of MoFEM namespace and are implemented in the core library. Others are generic for a particular problem, for example, mix transport formulation (MixTransport::MixTransportElement) is implemented in users modules. And finally, some classes (ExampleMix) are for very specific example, where some close solution is needed.

In class MixTransport::MixTransportElement all operators for evaluating matrices, the right-hand side vector, error evaluation are implemented. Moreover methods for declaring and defining fields, finite elements and problems are added. This class could be used by other users, solving similar problems.

In derived class ExampleMix, particular functions for given transport problem with h-adaptivity are implemented. The method handles specific boundary conditions for this example.

The general procedure to solve problem looks like follows

BcFluxMap bc_flux_map;
ExampleMix ufe(m_field,bc_flux_map);
// Initially calculate problem on coarse mesh
ierr = ufe.addFields("VALUES","FLUXES",order); CHKERRQ(ierr);
ierr = ufe.addFiniteElements("FLUXES","VALUES"); CHKERRQ(ierr);
// Set boundary conditions
ierr = ufe.addBoundaryElements(ref_level);
ierr = ufe.buildProblem(ref_level); CHKERRQ(ierr);
ierr = ufe.createMatrices(); CHKERRQ(ierr);
ierr = ufe.solveLinearProblem(); CHKERRQ(ierr);
ierr = ufe.calculateResidual(); CHKERRQ(ierr);
ierr = ufe.evaluateError(); CHKERRQ(ierr);
ierr = ufe.destroyMatrices(); CHKERRQ(ierr);
ierr = ufe.postProc("out_0.h5m"); CHKERRQ(ierr);

such procedure will look very similar to other problems, where we set-up problem, apply boundary conditions, solve problem and postprocess results.

Some MoFEM functions and philosophy of building problem is explained in minimal_surface_equation, please look into that document for further explanation.

Adding fields to the database (ufe.addFields)

In MixTransport::MixTransportElement::addFields we declare and define fields and add entities to those on which given field is spanning.

PetscErrorCode MixTransport::MixTransportElement::addFields(const std::string &values,const std::string &fluxes,const int order) {
//meshset consisting all entities in mesh
EntityHandle root_set = mField.get_moab().get_root_set();
//add entities to field
ierr = mField.add_ents_to_field_by_type(root_set,MBTET,fluxes); CHKERRQ(ierr);
ierr = mField.add_ents_to_field_by_type(root_set,MBTET,values); CHKERRQ(ierr);
ierr = mField.set_field_order(root_set,MBTET,fluxes,order+1); CHKERRQ(ierr);
ierr = mField.set_field_order(root_set,MBTRI,fluxes,order+1); CHKERRQ(ierr);
ierr = mField.set_field_order(root_set,MBTET,values,order); CHKERRQ(ierr);

Note that in first two lines we declare two fields FLUXES and VALUES. We define fields, by setting space to HDIV and L2, respectively. Third argument in functions MoFEM::Interface::add_field sets number of coefficients. Since we have vector and scalar field, number of coefficients is 1. If for example we would like to have second rank field for HDIV space and vector field for L2 space, number of coefficients would be three. Note that in general case you can set the field on part of the mesh and on some part it could be 3d field, 2d or 1d.

Function MoFEM::Interface::add_ents_to_field_by_TETs add entities to a field, in this case, we like to solve the problem on whole mesh, so all entities from root meshset are added. Since we span that field on tetrahedra and all lower entities field is in 3d. The dimension of field is dictated by a dimension of entities.

Adding elements to the database (ufe.addFiniteElements)

In the following part of the code we declare finite elements, implementation of elements is done separately. By definition of finite elements, we are simply telling code on what fields given finite element operates. Domain of elements, i.e. entities on which elements integrate operators could be declared on the part of the mesh where the field is given. If you declare element on entities on which fields is no spanning, that element will have zero number of degrees of freedom.

In this example, we have four kinds of finite elements MIX, MIX_SKELETON, MIX_BCVALUE and MIX_BCFLUX.

const std::string &fluxes_name,
const std::string &values_name,
const std::string mesh_nodals_positions = "MESH_NODE_POSITIONS"
) {
// Set up volume element operators. Operators are used to calculate components
// of stiffness matrix & right hand side, in essence are used to do volume integrals over
// tetrahedral in this case.
// Declare element "MIX". Note that this element will work with fluxes_name and
// values_name. This reflect bilinear form for the problem
ierr = mField.add_finite_element("MIX",MF_ZERO); CHKERRQ(ierr);
// Define element
ierr = mField.modify_finite_element_add_field_row("MIX",fluxes_name); CHKERRQ(ierr);
ierr = mField.modify_finite_element_add_field_col("MIX",fluxes_name); CHKERRQ(ierr);
ierr = mField.modify_finite_element_add_field_row("MIX",values_name); CHKERRQ(ierr);
ierr = mField.modify_finite_element_add_field_col("MIX",values_name); CHKERRQ(ierr);
ierr = mField.modify_finite_element_add_field_data("MIX",fluxes_name); CHKERRQ(ierr);
ierr = mField.modify_finite_element_add_field_data("MIX",values_name); CHKERRQ(ierr);
// Declare finite element to integrate over skeleton, we need that to evaluate error
ierr = mField.add_finite_element("MIX_SKELETON",MF_ZERO); CHKERRQ(ierr);
// Define element
ierr = mField.modify_finite_element_add_field_row("MIX_SKELETON",fluxes_name); CHKERRQ(ierr);
ierr = mField.modify_finite_element_add_field_col("MIX_SKELETON",fluxes_name); CHKERRQ(ierr);
ierr = mField.modify_finite_element_add_field_data("MIX_SKELETON",fluxes_name); CHKERRQ(ierr);
// In some cases you like to use HO geometry to describe shape of the body, curved edges and faces, for
// example body is a sphere. HO geometry is approximated by a field, which can be hierarchical, so shape of
// the edges could be given by polynomial of arbitrary order.
// Check if field "mesh_nodals_positions" is declared, and if it is add that field to data of finite
// element. MoFEM will use that that to calculate Jacobian as result that geometry in nonlinear.
if(mField.check_field(mesh_nodals_positions)) {
ierr = mField.modify_finite_element_add_field_data("MIX",mesh_nodals_positions); CHKERRQ(ierr);
ierr = mField.modify_finite_element_add_field_data("MIX_SKELETON",mesh_nodals_positions); CHKERRQ(ierr);
// Look for all BLOCKSET which are MAT_THERMALSET, takes entities from those BLOCKSETS
// and add them to "MIX" finite element. In addition get data from that meshset
// and set cOnductivity which is used to calculate fluxes from gradients of concentration
// or gradient of temperature, depending how you interpret variables.
// cerr << *it << endl;
Mat_Thermal temp_data;
ierr = it->getAttributeDataStructure(temp_data); CHKERRQ(ierr);
setOfBlocks[it->getMeshsetId()].cOnductivity = temp_data.data.Conductivity;
setOfBlocks[it->getMeshsetId()].cApacity = temp_data.data.HeatCapacity;
rval = mField.get_moab().get_entities_by_type(
); CHKERRQ(ierr);
Range skeleton;
rval = mField.get_moab().get_adjacencies(
); CHKERRQ(ierr);
// Declare element to integrate natural boundary conditions, i.e. set values.
ierr = mField.add_finite_element("MIX_BCVALUE",MF_ZERO); CHKERRQ(ierr);
// Define element
ierr = mField.modify_finite_element_add_field_row("MIX_BCVALUE",fluxes_name); CHKERRQ(ierr);
ierr = mField.modify_finite_element_add_field_col("MIX_BCVALUE",fluxes_name); CHKERRQ(ierr);
ierr = mField.modify_finite_element_add_field_data("MIX_BCVALUE",fluxes_name); CHKERRQ(ierr);
ierr = mField.modify_finite_element_add_field_data("MIX_BCVALUE",values_name); CHKERRQ(ierr);
if(mField.check_field(mesh_nodals_positions)) {
ierr = mField.modify_finite_element_add_field_data("MIX_BCVALUE",mesh_nodals_positions); CHKERRQ(ierr);
// Declare element to apply essential boundary conditions.
ierr = mField.add_finite_element("MIX_BCFLUX",MF_ZERO); CHKERRQ(ierr);
// Define element
ierr = mField.modify_finite_element_add_field_row("MIX_BCFLUX",fluxes_name); CHKERRQ(ierr);
ierr = mField.modify_finite_element_add_field_col("MIX_BCFLUX",fluxes_name); CHKERRQ(ierr);
ierr = mField.modify_finite_element_add_field_data("MIX_BCFLUX",fluxes_name); CHKERRQ(ierr);
ierr = mField.modify_finite_element_add_field_data("MIX_BCFLUX",values_name); CHKERRQ(ierr);
if(mField.check_field(mesh_nodals_positions)) {
ierr = mField.modify_finite_element_add_field_data("MIX_BCFLUX",mesh_nodals_positions); CHKERRQ(ierr);

Function MoFEM::Interface::add_finite_element add the element of given name to the database, the second argument indicates that if the element of such name already exists do not throw the error. When the second argument is not given default behaviour is that error is thrown if the element already exists.

Function MoFEM::Interface::modify_finite_element_add_field_row, MoFEM::Interface::modify_finite_element_add_field_col and MoFEM::Interface::modify_finite_element_add_field_data define on what fields given element operates. You can think that this corresponds how you define the bilinear form for your operator, what is the first and second argument. Also, you specify data, for example, permeability could be given as a separate field which needs to be evaluated at integration points to assemble matrix. In that case you add that information using MoFEM::Interface::modify_finite_element_add_field_data.

Functiond MoFEM::Interface::add_ents_to_finite_element_by_type add tetrahedral and triangles to given element respectively. You can think about tetrahedral and triangles as domains on which you perform the integration. Note that field can be spanned on tetrahedral (3d entities) but you can integrate on triangles, then on such element trace of field on the triangle is available from that element.

Note part of the code

which iterates over all BLOCKSETs on which MAT_THERMALSET is defined, then it gets information about material parameters and range of tetrahedral entities in that block. Next, it adds those tetrahedral to MIX finite element. In addition takes all adjacent triangles to tetrahedra and add those entities to MIX_SKELETON finite element.

Adding elements to MIX_BCVALUE and MIX_BCFLUX is implemented in next section. In above code, those finite elements are only declared and defined, but not yet implemented. Declaration, definition and implementation are independent of each other. You can use two different implementations of the same element.

Adding boundary elements to the database (ufe.addBoundaryElements)

In ExampleMix::addBoundaryElements loop over boundary BLOCKSETs is made to get triangles in that BLOCKSET. Then appropriate triangles are added to MIX_BCVALUE and MIX_BCFLUX

PetscErrorCode ExampleMix::addBoundaryElements(BitRefLevel &ref_level) {
Range tets;
ierr = mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(ref_level,BitRefLevel().set(),MBTET,tets);
Skinner skin(&mField.get_moab());
Range skin_faces; // skin faces from 3d ents
rval = skin.find_skin(0,tets,false,skin_faces); CHKERRQ_MOAB(rval);
// note: what is essential (dirichlet) is natural (neumann) for mix-FE compared to classical FE
Range natural_bc;
Range tris;
ierr = it->getMeshsetIdEntitiesByDimension(mField.get_moab(),2,tris,true); CHKERRQ(ierr);
HeatFluxCubitBcData mydata;
ierr = it->getBcDataStructure(mydata); CHKERRQ(ierr);
if(mydata.data.flag1==1) {
Range tris;
ierr = it->getMeshsetIdEntitiesByDimension(mField.get_moab(),2,tris,true); CHKERRQ(ierr);
bcFluxMap[it->getMeshsetId()].eNts = tris;
bcFluxMap[it->getMeshsetId()].fLux = mydata.data.value1;
// cerr << bcFluxMap[it->getMeshsetId()].eNts << endl;
// cerr << bcFluxMap[it->getMeshsetId()].fLux << endl;
Range essential_bc = subtract(skin_faces,natural_bc);
Range bit_tris;
ierr = mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(ref_level,BitRefLevel().set(),MBTRI,bit_tris);
essential_bc = intersect(bit_tris,essential_bc);
natural_bc = intersect(bit_tris,natural_bc);
ierr = mField.add_ents_to_finite_element_by_type(essential_bc,MBTRI,"MIX_BCFLUX"); CHKERRQ(ierr);
ierr = mField.add_ents_to_finite_element_by_type(natural_bc,MBTRI,"MIX_BCVALUE"); CHKERRQ(ierr);
// ierr = mField.add_ents_to_finite_element_by_type(skin_faces,MBTRI,"MIX_BCVALUE"); CHKERRQ(ierr);

Note that at the at the end of the above code; the skin from full mesh is taken, i.e. triangle elements on boundary of meshet. Then, triangles on the boundary (skin) on which values or fluxes are not set are used to set zero-flux boundary condition.

For standard finite element formulation, boundary condition on fluxes is natural, so if you have zero flux on the boundary, you have zero value for DOFs on the right-hand side vector. In essence, you do nothing to enforce zero flux. This is not the case here.

In addition, function

PetscErrorCode ExampleMix::getBcOnValues(
const EntityHandle ent,
const double x,const double y,const double z,
double &value) {
value = 0;


PetscErrorCode ExampleMix::getBcOnFluxes(
const EntityHandle ent,
const double x,const double y,const double z,
double &flux) {
if(lastEnt==ent) {
flux = lastFlux;
} else {
flux = 0;
for(BcFluxMap::iterator mit = bcFluxMap.begin();mit!=bcFluxMap.end();mit++) {
Range &tris = mit->second.eNts;
if(tris.find(ent)!=tris.end()) {
flux = mit->second.fLux;
lastEnt = ent;
lastFlux = flux;

are implemented. Those functions are used to evaluate values and boundary conditions at integration points. For more details how those functions are used look into MixTransport::MixTransportElement::OpRhsBcOnValues and MixTransport::MixTransportElement::OpEvaluateBcOnFluxes, those classes are implementations of operations on particular elements which are later explained.

Declaring and building problem (ufe.buildProblem)

Once fields and finite elements are declared and defined, we build/rebuild data structures carrying information needed to do the analysis. Note that function below is general, takes into account that problem could be build on several refinement levels, and if only the necessary part of structures if rebuild in subsequent mesh refinements.

The mesh MoFEM::BitRefLevel and related details with mesh refinement are described in the following sections. Here we focus only on declaring, building and partitioning problem for parallelized calculations.

MoFEM can work on partitioned meshes and parallel algebra and calculations. However, here we assume that mesh is the same on all processors only calculations, matrices and vectors are distributed.

//build field
ierr = mField.build_fields(); CHKERRQ(ierr);
// get tetrahedrons which has been build previously and now in so called garbage bit level
Range done_tets;
ierr = mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
); CHKERRQ(ierr);
ierr = mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
); CHKERRQ(ierr);
// get tetrahedrons which belong to problem bit level
Range ref_tets;
ierr = mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
); CHKERRQ(ierr);
ref_tets = subtract(ref_tets,done_tets);
ierr = mField.build_finite_elements("MIX",&ref_tets,2); CHKERRQ(ierr);
// get triangles which has been build previously and now in so called garbage bit level
Range done_faces;
ierr = mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
); CHKERRQ(ierr);
ierr = mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
); CHKERRQ(ierr);
// get triangles which belong to problem bit level
Range ref_faces;
ierr = mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
); CHKERRQ(ierr);
ref_faces = subtract(ref_faces,done_faces);
//build finite elements structures
ierr = mField.build_finite_elements("MIX_BCFLUX",&ref_faces,2); CHKERRQ(ierr);
ierr = mField.build_finite_elements("MIX_BCVALUE",&ref_faces,2); CHKERRQ(ierr);
ierr = mField.build_finite_elements("MIX_SKELETON",&ref_faces,2); CHKERRQ(ierr);
//Build adjacencies of degrees of freedom and elements
ierr = mField.build_adjacencies(ref_level); CHKERRQ(ierr);
//Declare problem
ierr = mField.add_problem("MIX",MF_ZERO); CHKERRQ(ierr);
//set refinement level for problem
ierr = mField.modify_problem_ref_level_set_bit("MIX",ref_level); CHKERRQ(ierr);
// Add element to problem
ierr = mField.modify_problem_add_finite_element("MIX","MIX"); CHKERRQ(ierr);
ierr = mField.modify_problem_add_finite_element("MIX","MIX_SKELETON"); CHKERRQ(ierr);
// Boundary conditions
ierr = mField.modify_problem_add_finite_element("MIX","MIX_BCFLUX"); CHKERRQ(ierr);
ierr = mField.modify_problem_add_finite_element("MIX","MIX_BCVALUE"); CHKERRQ(ierr);
//build problem
ProblemsManager *prb_mng_ptr;
ierr = mField.getInterface(prb_mng_ptr); CHKERRQ(ierr);
ierr = prb_mng_ptr->buildProblems(); CHKERRQ(ierr);
//mesh partitioning
ierr = prb_mng_ptr->partitionProblem("MIX"); CHKERRQ(ierr);
ierr = prb_mng_ptr->partitionFiniteElements("MIX"); CHKERRQ(ierr);
//what are ghost nodes, see Petsc Manual
ierr = prb_mng_ptr->partitionGhostDofs("MIX"); CHKERRQ(ierr);

In above code the several key functions are used. Function MoFEM::Interface::build_fields build field data structure by construction degrees of freedom on each entity on which given field is spanning. For example for L2 space DOFs are only on tetrahedrons which were added to the field. For H-div space is the more complex case, degrees of freedom are associated with triangles and tetrahedrons, but a number of degrees of each type of entity depend on approximation order, for details look to [1].

At that point, field structure is created including DOFs on entities. Such structure can carry information about approximation order and values of a field. For example, you can do algebraic operations on field values directly Field Basic Algebra.

No indices are attached to DOFs since indices depend on a number of processors on which solution is carried on, partitioning method or directly type of problem which you like to solve.

Function MoFEM::Interface::build_finite_elements is used to build data structure. You can note that in above code in subsequent mesh refinements elements are added to the database, no whole base is rebuilt from scratch. At that point, information about degrees of freedom of fields acting on the element is aggregated. Finite element is problem independent.

Once fields and finite element structures are build, we can build problem. Domain of the problem, its size and structure of the matrix is uniquely defined by finite elements and fields on which element operates. First, a problem is declared using MoFEM::Interface::add_problem, where the argument of the function is the name of the problem which we like to add. Next we define problem, by setting finite elements, this is done using function MoFEM::Interface::modify_problem_add_finite_element where elements are added to the problem.

To build problem, we need to set MoFEM::BitRefLevel to the problem. In short MoFEM::BitRefLevel allows selecting a refined mesh on which we like to solve our problem. The details and unique flexibility of MoFEM::BitRefLevel are explained in Explaining mesh BitRefLevel.

Declaration, definition and building of data structures, are independent of each other. We can first declare fields, finite elements and problems, then define at the end build data structures. Implementation of a finite element is independent of its definition.

Having all those stages done, at that point structure of algebraic equations emerging from problem definition, finite elements, fields and domains on which finite elements and fields operate. That defines the structure of sparse matrix and structure of parallel vectors which ghost values.

MoFEM utilizes properties of approximation spaces, knowing space of approximation field, we know how DOFs are distributed on the mesh. With information about the bilinear form, which is set in the defined of a finite element, we can build adjacencies between DOFs. At last, with the definition of a problem, i.e. set of finite elements in the problem, we can build unique enumeration of DOFs and build the system of algebraic equations. Using spaces (L2,H1,H-DIV,H-CURL) and their properties we build a system of a data structure which abstracts and general.

Create vectors and matrices (ufe.createMatrices)

At that point, we have defined and the partitioned problem and that is enough to construct distributed matrices and vectors. The sparsity of the matrix, adjacencies and other bookkeeping operations are handled by MoFEM and PETSc.

Functions and use of below code was explained earlier in minimal_surface_equation;

ierr = mField.getInterface<VecManager>()->vecCreateGhost("MIX",COL,&D); CHKERRQ(ierr);
ierr = mField.getInterface<VecManager>()->vecCreateGhost("MIX",COL,&D0); CHKERRQ(ierr);
ierr = mField.getInterface<VecManager>()->vecCreateGhost("MIX",ROW,&F); CHKERRQ(ierr);

Assemble and solve

At that point, we have all data structures necessary to solve the problem. Now we focus attention on the implementation of finite element instances and solution procedure. In essence, we need to assemble matrices and vectors, enforce constraints, solve a linear system of equations and post-process results.

Finite element instances

Definition of the finite element allows building the structure (adjacency) of the matrix, however, assemble (calculate) matrix we need finite element class implementation and its instance which will do the work. MoFEM allows implementing finite element at a low level, however, except special cases, a generic finite element class implementations can be used, where some operations like calculation of Jacobian, face areas and normals are implemented. Generic finite element classes enable us to focusses attention on the implementation of the bilinear form and its algebraic representation; that is numerical integration and evaluation of physical equations, evaluation of base functions and fields. The implementation of generic finite element is by users data operators, see tutorials Hello world and Solving the Poisson equation. Class of users data operator is inherited, and problem specific application is by function overloading. From users data operator class, we have access to base functions, integration points, local and global indices and values of DOFs. etc. UserDataOperator is run for entities on the element, on which approximation base is constructed.

In MixTransport::MixTransportElement class we create two generic finite element instances, i.e. MixTransport::MixTransportElement::MyVolumeFE and MixTransport::MixTransportElement::MyTriFE.

int getRule(int order) { return 2*order+1; };
MyVolumeFE feVol; ///> Instance of volume element
MyTriFE(MoFEM::Interface &m_field): MoFEM::FaceElementForcesAndSourcesCore(m_field) {}
int getRule(int order) { return 2*order+1; };
MyTriFE feTri; ///< Instance of surface element

Those two finite element classes have method MoFEM::ForcesAndSourcesCore::getRule overladed, which sets rank of integration rule on the element. The programmer needs to set integration rule, since he knows what will be integrated.

MoFEM base on that function and higher order of approximation base functions sets precalculated in integration points and weights. In some classes, the programmer would like to use own integration rule, that could be easily done as well, but we postpone explanation how to do this to another example. If you like to see how it works in practice, pleas have a look here Post Process.

Now we can add tasks (UsersOperators) to finite element, which will do a sequence of calculations. However, before we show how to implement UsersOperators, we show first how finite elements objects are used.

Solve problem (ufe.solveLinearProblem)

In the previous step, we have constructed matrices and vectors, next step is to calculate matrices and vectors and to solve the linear system of equations. At that point, we not yet focus how to integrate local finite element matrices, showing only how a solution of this problem is constructed. Numeration of DOFs, looping over finite elements, parallelization of the problem or assembly itself is handled by the MoFEM.

The function

can be divided into several stages. First, we will focus on the calculation of essential boundary conditions. We already defined boundary where flux is applied, however unlike in standard FE formulation where nodal shape functions are used, for H-div space we need to apriori calculate values of degrees of freedom for essential boundary conditions. Exploiting properties of H-div space, trace of base functions on the boundary can be evaluated element by element, where for each boundary triangle small system of linear equations is solved.

feTri.getOpPtrVector().push_back(new OpEvaluateBcOnFluxes(*this,"FLUXES",D0));
ierr = mField.loop_finite_elements("MIX","MIX_BCFLUX",feTri); CHKERRQ(ierr);

Note that first line in the code above clean vector of operators on the element instance (element object could be used before to do other tasks), second line add operator MixTransport::MixTransportElement::OpEvaluateBcOnFluxes to the finite element instance. The operator has an argument vector D0 in which DOFs values for essential boundary conditions are stored. The third line makes the loop over all elements in the problem (name of the problem is the first argument), for the finite element (second argument) and evaluate operators of element instance for object MixTransport::MixTransportElement::feTri. It is clear here that definition of element is separated from its implementation.

The similar procedure is repeated to calculate matrix and the right-hand side vector

feVol.getOpPtrVector().push_back(new MixTransport::MixTransportElement::OpL2Source(*this,"VALUES",F));
feVol.getOpPtrVector().push_back(new MixTransport::MixTransportElement::OpFluxDivergenceAtGaussPts(*this,"FLUXES"));
feVol.getOpPtrVector().push_back(new MixTransport::MixTransportElement::OpValuesAtGaussPts(*this,"VALUES"));
feVol.getOpPtrVector().push_back(new MixTransport::MixTransportElement::OpDivTauU_HdivL2(*this,"FLUXES","VALUES",Aij,F));
feVol.getOpPtrVector().push_back(new MixTransport::MixTransportElement::OpTauDotSigma_HdivHdiv(*this,"FLUXES",Aij,F));
feVol.getOpPtrVector().push_back(new MixTransport::MixTransportElement::OpVDivSigma_L2Hdiv(*this,"VALUES","FLUXES",Aij,F));
ierr = mField.loop_finite_elements("MIX","MIX",feVol); CHKERRQ(ierr);

where the calculation of matrix and right-hand side vector is broken down into several stages, where each sub-matrix is calculated by an user operator acting on finite element entities. Note that this time integration is over volume element, i.e. MixTransport::MixTransportElement::feVol. Finally, this part of code ends with the loop over all MIX finite elements in the MIX problem. If calculations are done in parallel, MoFEM manages related complexities transparently.

To complete calculations we have to integrate natural boundary conditions, i.e. enforce Dirichlet boundary condition

feTri.getOpPtrVector().push_back(new OpRhsBcOnValues(*this,"FLUXES",F));
ierr = mField.loop_finite_elements("MIX","MIX_BCVALUE",feTri); CHKERRQ(ierr);

This procedure looks very similar what was before; we clear existing operators on MixTransport::MixTransportElement::feTri, add operator MixTransport::MixTransportElement::OpRhsBcOnValues and loop over are surface elements were Dirichlet boundary is set.

At this point, we have created and assembled linear system of equations, we have not given an explanation how to do integration at user operator level, this will be explained later. At this point, we will show how to enforce essential boundary conditions. We have calculated values of DOFs already; now we have to modify matrix and the right-hand side vector

ierr = getDirichletBCIndices(&essential_bc_ids); CHKERRQ(ierr);
ierr = MatZeroRowsIS(Aij,essential_bc_ids,1,D0,F); CHKERRQ(ierr);

while operator MixTransport::MixTransportElement::OpEvaluateBcOnFluxes is executed it records indices of DOFs for which essential boundary conditions are applied. This DOFs are accessed by MixTransport::MixTransportElement::getDirichletBCIndices. We use those indices to zero rows of matrix Aij and insert 1 on diagonal, in addition we fill the right-hand side vector with previously calculated vector D0. This operation is performed with PETSc function MatZeroRowsIS (http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatZeroRowsIS.html)

At that point we can solve system of equations

KSP solver;
ierr = KSPCreate(PETSC_COMM_WORLD,&solver); CHKERRQ(ierr);
ierr = KSPSetOperators(solver,Aij,Aij); CHKERRQ(ierr);
ierr = KSPSetFromOptions(solver); CHKERRQ(ierr);
ierr = KSPSetUp(solver); CHKERRQ(ierr);
ierr = KSPSolve(solver,F,D); CHKERRQ(ierr);
ierr = KSPDestroy(&solver); CHKERRQ(ierr);

For details abut above code we refer to PETSc user manual https://www.mcs.anl.gov/petsc/documentation/index.html.

The last operation at that stage is to store results on the mesh

ierr = mField.getInterface<VecManager>()->setGlobalGhostVector("MIX",COL,D,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);

You can think about MoFEM as a link between algebraic structures (e.g. vector) and mesh. Function MoFEM::Interface::set_global_ghost_vector is an example of this, using it you can transfer field values from mesh to vector and reverse of it.

Note that discretization field is abstract has no indices, where vector has indexing by its nature. You can imagine that you solve the thermal problem. As a result, you have calculated temperature field, in other problem, solved with different partitioning, on some subdomain of thermal problem, you like to calculate mechanical field using a staggered method where previously calculated temperature is used. In that case, an algebraic structure (enumeration of DOFs) is different. However, you use MoFEM::Interface::set_global_ghost_vector transfer of field from one problem to another is straightforward.

Implementation of user operator

In MoFEM we apply general procedure, recognising in the typical finite analysis to calculate matrices and vectors we make loops. Loop over finite elements on given partition, loop over entities in the elements (e.g. triangle is composed of entities like vertices, edges and triangle itself) and build on those entities base functions. Finally loop over integration points, where base functions are evaluated and local element matrix assembled.

Performing those loops, MoFEM manages associated complexities, like extraction of global and local indices or DOFs values, evaluation of base functions and derivatives and other typical operations.

All this work is done by generic finite element, the basic finite element instance is built on MoFEM::FEMethod, then more advanced structure recognising that finite element entity is composed of lower dimension entities is MoFEM::ForcesAndSourcesCore. MoFEM::ForcesAndSourcesCore allows to add so-called users operators, which are classes operating on each lower dimension entities and performing loops over integration points.

From class MoFEM::ForcesAndSourcesCore are derived more specialised classes like MoFEM::VolumeElementForcesAndSourcesCore, MoFEM::FaceElementForcesAndSourcesCore and others, see Forces and sources . Those elements recognise the fact that for example working with a volume (3d) finite element you would like to calculate volume or use appropriate integration rule. From other hand working with faces like triangles, MoFEM::FaceElementForcesAndSourcesCore, you would like to know what is the area of the face or its normal and tangent vectors. Derived classes like MoFEM::VolumeElementForcesAndSourcesCore delivers tools for finite element programmer to do the specific job.

In this example, we use functionality of MoFEM::VolumeElementForcesAndSourcesCore by creating new class MixTransport::MixTransportElement::MyVolumeFE that derives from the existing class. The derived class inherits the properties of the base class, and you can add or override methods and properties as required. Each class has method MoFEM::ForcesAndSourcesCore::getOpPtrVector, returning reference to vector of shared pointers to MoFEM::ForcesAndSourcesCore::UserDataOperator, i.e.

boost::ptr_vector<MoFEM::ForcesAndSourcesCore::UserDataOperator>& MoFEM::ForcesAndSourcesCore::getOpPtrVector ();

Instances of operator classes derived from MoFEM::DataOperator. MoFEM::ForcesAndSourcesCore::UserDataOperator do particular jobs, for example, calculate field values at integration points, integrate mass or stiffness matrix. Finite element instance can run the sequence of user operator in order you add them to shared common data structure, here MixTransport::MixTransportElement itself. So for example in one operator you can calculate gradients of displacement field at each integration point, in another, you can symmetrize gradient to calculate strain and in following calculate stress and another calculate the vector of internal forces. You break down the big job into smaller components, which you can use in several different contexts and easily test them one by one, making debugging easy. That is general philosophy of MoFEM, that complex and difficult problems are composition of small and simple tasks.

Class MoFEM::ForcesAndSourcesCore::UserDataOperator is inherited by MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator or MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator and other users operator classes associated with vertices and edges. User operator class has two types of functions which need to be overloaded to do calculations, and that is

virtual PetscErrorCode MoFEM::DataOperator::doWork( int row_side,int
col_side, EntityType row_type,EntityType col_type,
DataForcesAndSourcesCore::EntData &row_data, DataForcesAndSourcesCore::EntData
&col_data );


virtual PetscErrorCode
MoFEM::DataOperator::doWork( int side, EntityType type,
DataForcesAndSourcesCore::EntData &data )

You have seen how we use those operators in the previous section here we show how to implement particular cases. The first doWork(nt row_side,int col_side,...) function is usually overloaded when you assemble matrix, the second is often used to assemble vectors. However, use of users operators is beyond assembly of matrices and vectors. Note that fist doWork(int side,EntityType type,...) functions has as argument information about data in rows and columns, were second takes as an argument only rows or columns, depending on how user operator is set-up.

Assembling matrix

We will show how user operators work on two examples. First, we focus attention on term

\[ b_\mathbf{C}(\sigma,v) = \int_\Omega \textrm{div}[\boldsymbol\sigma] v \textrm{d}\Omega = \sum_\mathcal{T}^E \int_{\mathcal{T}^e} \textrm{div}[\boldsymbol\sigma] v \textrm{d}\mathcal{T} = \mathbf{q}_\textrm{h-div}^\textrm{T} \left( \sum_\mathcal{T}^E \sum_\mathcal{E_\textrm{row}}^e \sum_\mathcal{E_\textrm{col}}^e \sum_{g=0}^\textrm{N} w_g B^\textrm{h-div}_{i,i} B^\textrm{v} \right) \mathbf{q}_\textrm{v} = \left(\mathbf{q}_\textrm{h-div}\right)^\textrm{T} \mathbf{C} \mathbf{q}_\textrm{v} \]

were we implement matrix \(\mathbf{C}\)

\[ \mathbf{C} = \sum_\mathcal{T}^E \sum_\mathcal{E_\textrm{row}}^e \sum_\mathcal{E_\textrm{col}}^e \sum_{g=0}^\textrm{N} w_g B^\textrm{h-div}_{i,i} B^\textrm{v} \]

Note that in this case, where we have H-div space on rows and L2 space on columns, we have

\[ \mathcal{E}_\textrm{row}=\{ TRIANGLES\,\times\,4, TETRAHERAL \},\quad \mathcal{E}_\textrm{col} = \{ TETRAHEDRAL \} \]

You can notice three loops, over elements, entities on element and integration point. The first and second loop are managed by MoFEM the last one over integration points where a particular operator is evaluated is implemented by programmer, as shown below

Mat Aij;
Aij(aij) {
//this operator is not symmetric setting this variable makes element
//operator to loop over element entities without
//assumption that off-diagonal matrices are symmetric.
sYmm = false;
PetscErrorCode doWork(
int row_side,int col_side,
EntityType row_type,EntityType col_type,
DataForcesAndSourcesCore::EntData &row_data,
DataForcesAndSourcesCore::EntData &col_data
) {
int nb_row = row_data.getFieldData().size();
int nb_col = col_data.getFieldData().size();
if(nb_row==0) MoFEMFunctionReturnHot(0);
if(nb_col==0) MoFEMFunctionReturnHot(0);
int nb_gauss_pts = row_data.getHdivN().size1();
for(int gg = 0;gg<nb_gauss_pts;gg++) {
double w = getGaussPts()(3,gg)*getVolume();
); CHKERRQ(ierr);
noalias(C) += w*outer_prod(divVec,row_data.getN(gg));
ierr = MatSetValues(
); CHKERRQ(ierr);

Note that we integrate over tetrahedra. Thus operator MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator is inherited. The class constructor MixTransport::MixTransportElement::OpVDivSigma_L2Hdiv::OpVDivSigma_L2Hdiv takes an argument PETSc matrix and set variable sYmm to false, since matrix \(\mathbf{C}\) is non-symmetric. Moreover, constructing operator we indicate that in the first argument of bilinear form we take field FLUXES and in the second argument of bilinear form we take field VALUES. Finally, in the constructor of operator we indicate that we like to loop over all possible combinations of entities on rows and columns, by setting OPROWCOL from enum type MoFEM::ForcesAndSourcesCore::UserDataOperator::OpType.

The real work is done in

int row_side,int col_side,
EntityType row_type,EntityType col_type,
DataForcesAndSourcesCore::EntData &row_data,
DataForcesAndSourcesCore::EntData &col_data

MoFEM executes this method, in this particular case, for each combination of pairs

\[ \mathcal{E}_\textrm{row} \times \mathcal{E}_\textrm{col} = \{ \left(TRI_0, TET\right), \left(TRI_1, TET\right), \left(TRI_2, TET\right), \left(TRI_4, TET\right), \left(TET, TET\right) \} \]

Note that function MoFEM::DataOperator::doWork takes as argument row_side and col_side, which represent a local index of processed entity on rows and columns, respectively. Also, row_type and col_type shows the type of entity, for example, tetrahedra or triangle in this case. To see more details about a type of entities, local indices and canonical numeration see [42].

The row_data and col_data are data structures MoFEM::DataForcesAndSourcesCore::EntData carrying information about values of base functions, global and local DOFs indices, DOFs values on evaluated entities. For example, you can get indices on the row using row_data.getIndices(), or matrix of base functions at integration point as follows

MatrixDouble &base_row_hdiv = row_data.getHdivN(gg);

In this case, each row of the matrix base_row_hdiv represents base function, where columns represent values of the element of base function vector. Look to MoFEM::DataForcesAndSourcesCore::EntData to see a range of data which is accessible by this data structure.

From inside MoFEM::VolumeElementForcesAndSourcesCore::UserDataOperator we can access more finite element specific information, in that case volume element, for example,

double vol = getVolume(); // get volume of element
MarixDouble &coords_at_gauss_pts = getCoordsAtGaussPts(); // physical coordinates of integration points

and much more.

After above explanation, we focus on implementation itself, so finally we loop over integration points and calculate local matrix C, as follows

int nb_gauss_pts = row_data.getHdivN().size1();
for(int gg = 0;gg<nb_gauss_pts;gg++) {
double w = getGaussPts()(3,gg)*getVolume();
ierr = getDivergenceOfHDivBaseFunctions(
); CHKERRQ(ierr);
noalias(C) += w*outer_prod(divVec,row_data.getN(gg));

Here we are using Boost uBlas functions to do operations on local matrices. You can alternatively use tensor library implemented in MoFEM (FTensor) to do the job. Once local element matrix C is calculated, it is assembled into global matrix

ierr = MatSetValues(
); CHKERRQ(ierr);

where global indices on entities on rows and columns are easily retrieved by

VectorInt &global_row_indices = row_data.getIndices();
VectorInt &global_col_indices = col_data.getIndices();

where MatSetValues is a PETSc function documented here http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatSetValues.html

To solve problem other blocks of matrix need to be assembled, for details how is done please see MixTransport::MixTransportElement.

Assembling vector

At this stage you could see how to assemble vector, we focus attention on the assembly of Dirichlet boundary condition. Note the Dirichlet boundary conditions in this formulation, is natural, i.e. enforced in a weak sense when a system of linear equations is solved.

This time, we will calculate integral

\[ \mathbf{f} = \int_{\Gamma_\textrm{u}} \overline{u} \tau_{i} n_i \textrm{d}\Gamma = \sum_\mathcal{T}^E \int_{\partial\mathcal{T}_\textrm{u}} \overline{u} \tau_{i} n_i \textrm{d}\partial\mathcal{T} = \sum_\mathcal{T}^E \sum_\mathcal{E}^e \sum_{g=0}^N w_g \overline{u} B^\textrm{h-div}_{i} n_i \]

to do that we use following user operator class

Vec F;
F(f) {}
PetscErrorCode doWork(int side,EntityType type,DataForcesAndSourcesCore::EntData &data) {
if(data.getFieldData().size()==0) MoFEMFunctionReturnHot(0);
int nb_gauss_pts = data.getHdivN().size1();
for(int gg = 0;gg<nb_gauss_pts;gg++) {
const double x = getCoordsAtGaussPts()(gg,0);
const double y = getCoordsAtGaussPts()(gg,1);
const double z = getCoordsAtGaussPts()(gg,2);
double value;
ierr = cTx.getBcOnValues(x,y,z,value); CHKERRQ(ierr);
double w = getGaussPts()(2,gg)*0.5;
noalias(Nf) += w*prod(data.getHdivN(gg),getNormal())*value;
ierr = VecSetValues(F,data.getIndices().size(),&data.getIndices()[0],&Nf[0],ADD_VALUES); CHKERRQ(ierr);

Note that above class is derived from MoFEM::FaceElementForcesAndSourcesCore::UserDataOperator, since we integrate over triangles where Dirichlet boundary conditions are set.

In this case, the constructor has two arguments, reference to MixTransport::MixTransportElement object itself, and the right-hand side vector. We use cTx reference to access function evaluating boundary conditions. Constructing derived class we are indicating that use field FLUXES and loop only over entities on the row, by setting OPROW from enum type MoFEM::ForcesAndSourcesCore::UserDataOperator::OpType.

The real work is done in

int side,EntityType type,DataForcesAndSourcesCore::EntData &data

Since we loop only over row entities, the function doWork, looks different that that one to integrate matrix, it is simpler. It passes argument about side which is local number of entity for which operator is executed, type of entity and finally data structure

MoFEM executes this method, in this particular case, for only one entity, i.e.

\[ \mathcal{E}=\{ TRI \} \]

This is a consequence that testing functions on rows are in H-div space. Note that if the testing function would be in H-curl space, so for triangle loop on entities will involve three edges and triangle itself, from a testing function in H1 space loop will involve vertices, edges and triangle itself.

The integration takes place in this lines of the code

for(int gg = 0;gg<nb_gauss_pts;gg++) {
const double x = getCoordsAtGaussPts()(gg,0);
const double y = getCoordsAtGaussPts()(gg,1);
const double z = getCoordsAtGaussPts()(gg,2);
double value;
ierr = cTx.getBcOnValues(x,y,z,value); CHKERRQ(ierr);
double w = getGaussPts()(2,gg)*0.5;
noalias(Nf) += w*prod(data.getHdivN(gg),getNormal())*value;

where getNormal() get normal of the face and cTx.getBcOnValues(x,y,z,value) calculates value of boundary condition at coordinates of integration point.

Finally, local entity local vector is assembled,

ierr = VecSetValues(F,data.getIndices().size(),&data.getIndices()[0],&Nf[0],ADD_VALUES); CHKERRQ(ierr);

where VecSetValues is PETSc function, see http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Vec/VecSetValues.html for details.

Evaluate error

In this example to estimate posterior error we follow paper [11]. This error estimator works well for mesh-dependent norms, for which is reliable and efficient. However, it is not efficient in natural norm of \(H(\textrm{div},\Omega)\times L^2(\Omega)\) [13].

You can find two other papers, [13] and [4] which have efficient error estimator working in the natural norm. Those errors estimators are not very different from an implementation point of view and could be easily added here.
To do numerical experiment testing error bound we can solve the problem from [18] where analytical solutions is given. That will involve modification how boundary conditions are applied. Again, this is easy to implement, so if you are willing to add this to this documentation, it is very welcome.
Based on this work, you can write a good journal article about the efficiency of error estimator, taking into account two above remarks. If you are interested, we are happy to help. Extension to hp-adaptivity is straightforward.
A useful add-on is to construct preconditioner, which will take into account hierarchy created during hp-adaptivity, this involves some work, again we are happy to help.

Following paper [11], we consider bounds of the error estimator in respect to the mesh-dependent norms

\[ \| \boldsymbol\tau \|_{0,h} := \left( \| \boldsymbol\tau \|_0^2 + h \sum_{e\in\Gamma_h} \| \boldsymbol\tau \mathbf{n} \|^2_{0,e} \right)^{1/2} \]


\[ |v|_{1,h} := \left( \sum_{T\in\mathcal{T}_h} |v|^2_{1,T} + h^{-1}\sum_{e\in\Gamma_h} | J(v) |_{0,e}^2 \right)^{1/2} \]

where \(J(v)\) is a jump on the face \(e\) and h is is consistently replaced by the diameters of the actual elements and the lengths of the edges. \(\|\cdot\|_0\) and \(\|\cdot\|_1\) are Sobolev norms and semi-norm, respectively.

Our a posteriori estimates will be based on residuals. There will be contributions from the tetrahedrons as

\[ \begin{split} \eta_{T,R,1} &:= \| \mathbf{A} \boldsymbol\sigma + \textrm{grad}[\mathbf{u}] \|_{0,T} \\ \eta_{T,R,2} &:= h\| \textrm{div}[\boldsymbol\sigma] -\mathbf{f} \|_{0,T} \end{split} \]

and from the jumps on the inter-element boundaries

\[ \eta_{e,R} := h^{-1/2} \| J(u) \|_{0,e} \]

The local error estimator is a weighted combination

\[ \eta_{T,R} := \left( \eta^2_{T,R,1} + \eta^2_{T,R,2} + \sum_{e \in \partial T} \eta_{e,R}^2 \right)^{1/2} \]

For convenience, we will also refer to the global sum

\[ \eta_{T,R} := \left( \sum_{T \in \mathcal{T}_h} \left[ \eta^2_{T,R,1} + \eta^2_{T,R,2} \right] + \sum_{e \in \partial T} \eta_{e,R}^2 \right)^{1/2} \]

This error estimator is reliable and efficient, is upper bound is

\[ \|\boldsymbol\sigma^*-\boldsymbol\sigma\|_{0,h} + | u^* - u |_{1,h} \leq \frac{c}{1-\beta} \eta_R \]

where \(c\) depends on the shape parameter of \(\mathcal{T}_h\) and \(\beta\) is saturation parameter, see [11] for details. Lower bound estimates is

\[ \eta_{T,R} \leq \| \boldsymbol\sigma^* - \boldsymbol\sigma \|_{0,h} + | u^* - u |_{1,h} + h \| f - \textrm{div}[\boldsymbol\sigma] \|_0 \]

Calculate error and refine mesh

Once we have residual error estimator we could do h-adaptivity, it is done by executing following code

for(int ll = 1;ll!=nb_levels;ll++) {
const int nb_levels = ll;
ierr = ufe.squashBits(); CHKERRQ(ierr);
ierr = ufe.refineMesh(ufe,nb_levels,order); CHKERRQ(ierr);
ref_level = BitRefLevel().set(nb_levels);
ierr = ufe.addBoundaryElements(ref_level);
ierr = ufe.buildProblem(ref_level); CHKERRQ(ierr);
ierr = ufe.createMatrices(); CHKERRQ(ierr);
ierr = ufe.solveLinearProblem(); CHKERRQ(ierr);
ierr = ufe.calculateResidual(); CHKERRQ(ierr);
ierr = ufe.evaluateError(); CHKERRQ(ierr);
ierr = ufe.destroyMatrices(); CHKERRQ(ierr);
ierr = ufe.postProc(
(std::ostringstream().seekp(0) << "out_" << nb_levels << ".h5m").str()
); CHKERRQ(ierr);

where nb_levels is number of refinement levels.

We focus attention on mesh refinement procedures in next section, but first we will show how to calculate residual error estimator using users operators. Function evaluating error is implemented as follows

feTri.getOpPtrVector().push_back(new MixTransport::MixTransportElement::OpSkeleton(*this,"FLUXES"));
ierr = mField.loop_finite_elements("MIX","MIX_SKELETON",feTri,0,mField.get_comm_size()); CHKERRQ(ierr);
feVol.getOpPtrVector().push_back(new MixTransport::MixTransportElement::OpFluxDivergenceAtGaussPts(*this,"FLUXES"));
feVol.getOpPtrVector().push_back(new MixTransport::MixTransportElement::OpValuesGradientAtGaussPts(*this,"VALUES"));
feVol.getOpPtrVector().push_back(new MixTransport::MixTransportElement::OpError(*this,"VALUES"));
ierr = mField.loop_finite_elements("MIX","MIX",feVol,0,mField.get_comm_size()); CHKERRQ(ierr);
const Problem *problem_ptr;
ierr = mField.get_problem("MIX",&problem_ptr); CHKERRQ(ierr);
"Nb dofs %d error flux^2 = %6.4e error div^2 = %6.4e error jump^2 = %6.4e error tot^2 = %6.4e\n",

Once the problem is solved, first we need to calculate jump on the skeleton. We simply doing this applying the same procedure like before, adding appropriate entity operator to the element object and looping over all triangular elements in the problem, as follows

feTri.getOpPtrVector().push_back(new MixTransport::MixTransportElement::OpSkeleton(*this,"FLUXES"));
ierr = mField.loop_finite_elements("MIX","MIX_SKELETON",feTri,0,mField.get_comm_size()); CHKERRQ(ierr);

Note that entity operator MixTransport::MixTransportElement::OpSkeleton is executed for every element MIX_SKELETON. How this entity operator is implemented is described in following section.

Once the jump is calculated, we can evaluate error estimator, by running sequence of entity operators on each volume element, as follows,

feVol.getOpPtrVector().push_back(new MixTransport::MixTransportElement::OpFluxDivergenceAtGaussPts(*this,"FLUXES"));
feVol.getOpPtrVector().push_back(new MixTransport::MixTransportElement::OpValuesGradientAtGaussPts(*this,"VALUES"));
feVol.getOpPtrVector().push_back(new MixTransport::MixTransportElement::OpError(*this,"VALUES"));
ierr = mField.loop_finite_elements("MIX","MIX",feVol,0,mField.get_comm_size()); CHKERRQ(ierr);

Once we loop over all elements and calculate error, we aggregate all element contributions in sumErrorFlux and sumErrorDiv and sumErrorJump which are \(\eta^2_{T,R,1}\), \(\eta^2_{T,R,2}\) and \(\eta^2_{e,R}\), respectively.

Calculate jump on the skeleton

At this point, we will describe how to integrate over skeleton to get jumps on the faces. That involves loop over are triangles and then for each triangle visiting adjacent tetrahedral, one by one.

On each triangle, we set some integration rule, which integration points are projected on appropriate sides of two adjacent tetrahedra. Next, for each tetrahedron, base functions (and derivatives) are evaluated at integration points and values of function \(u\) calculated. Once values on each tetrahedron are computed we do loop over all integration points on face and calculate norm of the jump, i.e. \(\eta_{e,R}\).

MoFEM has some dedicated data structures and methods to do this type of integration. Those procedures were implemented to give user ability to implement Discontinuous (Petrov) Galerkin, contact mechanics or other approaches involving integration on skeleton where information from adjacent entities is needed.

* \brief volume element to get values from adjacent tets to face
/** store values at integration point, key of the map is sense of face in
* respect to adjacent tetrahedra
map<int,VectorDouble> valMap;
* \brief calculate values on adjacent tetrahedra to face
struct OpVolSide;
MixTransport::MixTransportElement &ctx,const std::string flux_name
cTx(ctx) {
volSideFe.getOpPtrVector().push_back(new OpSkeleton::OpVolSide(valMap));
PetscErrorCode doWork(int side,EntityType type,DataForcesAndSourcesCore::EntData &data);

User entity operator MixTransport::MixTransportElement::OpSkeleton has nested volume element and nested entity operator to integrate over it. Note that nested volume element is a class MoFEM::VolumeElementForcesAndSourcesCoreOnSide, that is inherited from MoFEM::VolumeElementForcesAndSourcesCore. This derived class has additional functionality, which if an instance of it is called from face operator, it project face integration points on the appropriate face of adjacent tetrahedra. As results volume, base functions are evaluated at integration points which come from currently evaluated skeleton face.

That nested volume element, here having instance volSideFe, run nested operator class MixTransport::MixTransportElement::OpSkeleton::OpVolSide, this class is derived form MoFEM::VolumeElementForcesAndSourcesCoreOnSide::UserDataOperator. Again this class give access to additional information useful when integration is over skeleton, like face normal vector, face area, etc. Note that nested operator MixTransport::MixTransportElement::OpSkeleton::OpVolSide is added to nested volume instance volSideFe, in constructor of face operator with the command line

volSideFe.getOpPtrVector().push_back(new MixTransport::MixTransportElement::OpSkeleton::OpSkeleton::OpVolSide(valMap));

The example of implementation of nested volume element operator looks as follows

map<int,VectorDouble> &valMap;
OpVolSide(map<int,VectorDouble> &val_map):
valMap(val_map) {
PetscErrorCode doWork(int side, EntityType type,DataForcesAndSourcesCore::EntData &data) {
if(data.getFieldData().size() == 0) MoFEMFunctionReturnHot(0);
int nb_gauss_pts = data.getN().size1();
for(int gg = 0;gg<nb_gauss_pts;gg++) {
valMap[getFaceSense()][gg] = inner_prod(trans(data.getN(gg)),data.getFieldData());

Implementation of doWork function is straightforward. It loops over integration points and computes function values using base functions and DOFs values. Results are stored in the map valMap. Note that internal face has two adjacent volume elements, whereas face on boundary only one. Each face has a unique orientation (sense) in respect to face, information about that sense is accessed by getFaceSense().

Once we evaluated approximation functions values at each face, we can calculate error \(\eta_{e,R}\), using face operator

PetscErrorCode MixTransport::MixTransportElement::OpSkeleton::doWork(int side,EntityType type,DataForcesAndSourcesCore::EntData &data) {
if(type == MBTRI) {
double def_val = 0;
Tag th_error_jump;
rval = cTx.mField.get_moab().tag_get_handle(
double* error_jump_ptr;
rval = cTx.mField.get_moab().tag_get_by_ptr(
th_error_jump,&fe_ent,1,(const void**)&error_jump_ptr
*error_jump_ptr = 0;
// check if this is essential boundary condition
EntityHandle essential_bc_meshset = cTx.mField.get_finite_element_meshset("MIX_BCFLUX");
if(cTx.mField.get_moab().contains_entities(essential_bc_meshset,&fe_ent,1)) {
// essential bc, np jump then, exit and go to next face
// calculate values form adjacent tets
ierr = loopSideVolumes("MIX",volSideFe); CHKERRQ(ierr);
int nb_gauss_pts = data.getHdivN().size1();
// it is only one face, so it has to be bc natural boundary condition
if(valMap.size()==1) {
if(valMap.begin()->second.size()!=nb_gauss_pts) {
SETERRQ(PETSC_COMM_WORLD,MOFEM_DATA_INCONSISTENCY,"wrong number of integration points");
for(int gg = 0;gg!=nb_gauss_pts;gg++) {
const double x = getCoordsAtGaussPts()(gg,0);
const double y = getCoordsAtGaussPts()(gg,1);
const double z = getCoordsAtGaussPts()(gg,2);
double value;
ierr = cTx.getBcOnValues(fe_ent,x,y,z,value); CHKERRQ(ierr);
double w = getGaussPts()(2,gg);
w *= getArea();
*error_jump_ptr += w*pow(value-valMap.begin()->second[gg],2);
} else if(valMap.size()==2) {
for(int gg = 0;gg!=nb_gauss_pts;gg++) {
double w = getGaussPts()(2,gg);
w *= getArea();
double delta = valMap.at(1)[gg]-valMap.at(-1)[gg];
*error_jump_ptr += w*pow(delta,2);

Note that we have to consider three scenarios, first, on the face is essential boundary condition. In that case, jump on face values is zero. If that is not a case, we run loop over adjacent volume elements

ierr = loopSideVolumes("MIX",volSideFe); CHKERRQ(ierr);

At this point nested volume element instance is called and then operator on this element executed. Results of this operation are stored in valMap.

If it is not case one, and if a face has only one volume in the side, thus on that face natural boundary is applied, then jump is a difference between boundary value and values evaluated by adjacent volume operator. If a face is in volume, i.e. third case, we calculate the difference between values from adjacent volumes. Results are stored on face tag values, here we are using MOAB to do that for us, using

rval = cTx.mField.get_moab().tag_get_by_ptr(
th_error_jump,&fe_ent,1,(const void**)&error_jump_ptr

which returns pointer to value stored (tagged) on face. Look here http://ftp.mcs.anl.gov/pub/fathom/moab-docs/contents.html#twofour to learn about tags in MOAB and how to use them.

Calculate residual error estimator

Once we have implemented how to calculate jumps between adjacent volume elements we can focus on computing residual error estimator. We do that in standard way using operator which acts on volume element

const std::string field_name
cTx(ctx) {}
PetscErrorCode doWork(
int side,EntityType type,DataForcesAndSourcesCore::EntData &data
) {
if(type != MBTET) MoFEMFunctionReturnHot(0);
int nb_gauss_pts = data.getN().size1();
double def_val = 0;
Tag th_error_flux,th_error_div;
rval = cTx.mField.get_moab().tag_get_handle(
double* error_flux_ptr;
rval = cTx.mField.get_moab().tag_get_by_ptr(
th_error_flux,&fe_ent,1,(const void**)&error_flux_ptr
rval = cTx.mField.get_moab().tag_get_handle(
double* error_div_ptr;
rval = cTx.mField.get_moab().tag_get_by_ptr(
th_error_div,&fe_ent,1,(const void**)&error_div_ptr
Tag th_error_jump;
rval = cTx.mField.get_moab().tag_get_handle(
double* error_jump_ptr;
rval = cTx.mField.get_moab().tag_get_by_ptr(
th_error_jump,&fe_ent,1,(const void**)&error_jump_ptr
*error_jump_ptr = 0;
/// characteristic size of the element
const double h = pow(getVolume()*12/sqrt(2),(double)1/3);
for(int ff = 0;ff!=4;ff++) {
rval = cTx.mField.get_moab().side_element(fe_ent,2,ff,face); CHKERRQ_MOAB(rval);
double* error_face_jump_ptr;
rval = cTx.mField.get_moab().tag_get_by_ptr(
th_error_jump,&face,1,(const void**)&error_face_jump_ptr
*error_face_jump_ptr = (1/sqrt(h))*sqrt(*error_face_jump_ptr);
*error_face_jump_ptr = pow(*error_face_jump_ptr,2);
*error_jump_ptr += *error_face_jump_ptr;
*error_flux_ptr = 0;
*error_div_ptr = 0;
for(int gg = 0;gg<nb_gauss_pts;gg++) {
double w = getGaussPts()(3,gg)*getVolume();
const double x = getCoordsAtGaussPts()(gg,0);
const double y = getCoordsAtGaussPts()(gg,1);
const double z = getCoordsAtGaussPts()(gg,2);
double flux;
ierr = cTx.getSource(fe_ent,x,y,z,flux); CHKERRQ(ierr);
*error_div_ptr += w*pow(cTx.divergenceAtGaussPts[gg]-flux,2);
ierr = cTx.getResistivity(fe_ent,x,y,z,invK); CHKERRQ(ierr);
noalias(deltaFlux) = prod(invK,cTx.fluxesAtGaussPts[gg])+cTx.valuesGradientAtGaussPts[gg];
*error_flux_ptr += w*inner_prod(deltaFlux,deltaFlux);
*error_div_ptr = h*sqrt(*error_div_ptr);
*error_div_ptr = pow(*error_div_ptr,2);
cTx.sumErrorFlux += *error_flux_ptr;
cTx.sumErrorDiv += *error_div_ptr;
// FIXME: Summation should be while skeleton is calculated
cTx.sumErrorJump += *error_jump_ptr; // FIXME: this need to be fixed
cTx.errorMap[sqrt(*error_flux_ptr+*error_div_ptr+*error_jump_ptr)] = fe_ent;

The scheme how we implement is similar to those above, we integrate values, save them on volume tags, and aggregate errors from individual elements. Note that we reading tags values of previously calculated error from adjacent to tetrahedra faces, summing them up to calculate \( \sum_{e \in \partial T} \eta_{e,R}^2 \).

Mesh refinement

We will build a refined mesh from scratch, each time we add the new level. Using bit refinement levels (see Explaining mesh BitRefLevel) we can do that efficiently.

Since we build mesh without hanging nodes, not all edges in coarse tetrahedra are always refined. Refining such tetrahedra in subsequent refinements would result to poorly shaped elements. We will show how to overcome this difficulty.

Implementation of mesh refinement is in ExampleMix::refineMesh. First, we have to create a set of edges which are going to be split. This set is built from edges which were refined in previous mesh refinement levels and new edges which are adjacent to tetrahedra which have error above threshold. We used the standard "greedy" strategy for h-adaptive refinements, i.e. all elements which contribute 25% of the maximum element contribution to the square of total error in energy norm, are marked for the refinement.

To get edges which have been refined in previous steps we run

const RefEntity_multiIndex *ref_ents_ptr;
// get pointer to multi-index in databse having all entities
ierr = mField.get_ref_ents(&ref_ents_ptr); CHKERRQ(ierr);
typedef RefEntity_multiIndex::index<Composite_EntType_and_ParentEntType_mi_tag>::type RefEntsByComposite;
const RefEntsByComposite &ref_ents = ref_ents_ptr->get<Composite_EntType_and_ParentEntType_mi_tag>();
RefEntsByComposite::iterator rit,hi_rit;
// get all entities which parenrt are vertices and parent is edges.
rit = ref_ents.lower_bound(boost::make_tuple(MBVERTEX,MBEDGE));
hi_rit = ref_ents.upper_bound(boost::make_tuple(MBVERTEX,MBEDGE));
// get parent edges, those edges were refined in previous steps
Range refined_edges;
// thist loop is over vertices which parent is edge
for(;rit!=hi_rit;rit++) {
refined_edges.insert((*rit)->getParentEnt()); // get parent edge

To get entities which are deemed to refine at this level, we first get all elements which have error of 25% of the maximum element contribution

// get tets which has large error
Range tets_to_refine;
const double max_error = ufe.errorMap.rbegin()->first;
map<double,EntityHandle>::iterator mit = ufe.errorMap.begin();
) {
if(mit->first<0.25*max_error) continue;

where ufe.errorMap is filled in in MixTransport::MixTransportElement::OpError .

Once we have elements we can find adjacent edges

Range tets_to_refine_edges;
rval = mField.get_moab().get_adjacencies(

and merge edges into the set of all edges which are refined or going to be refined


At that point we have edges to refine and elements, we can do mesh refinement. Mesh refinement is managed by interface MoFEM::MeshRefinement, we can query for that interface as follows,

ierr = mField.getInterface(refine_ptr); CHKERRQ(ierr);

with this interface at hand, we do

for(int ll = 0;ll!=nb_levels;ll++) {
// get edges from previous mesh level
Range edges;
ierr = mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
); CHKERRQ(ierr);
// get intersection of all edges to refine and edges of previous level
edges = intersect(edges,refined_edges);
// refine edges by creation of node in middle of edge, if node was previously created use that node, and set BitRefLevel to that node
ierr = refine_ptr->add_verices_in_the_middel_of_edges(edges,BitRefLevel().set(ll+1)); CHKERRQ(ierr);
// get all tetrahderals for current refined level
Range tets;
ierr = mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
); CHKERRQ(ierr);
// refine those tetrahedrals using nodes in middle of refined edges
ierr = refine_ptr->refine_TET(tets,BitRefLevel().set(ll+1)); CHKERRQ(ierr);
// update meshsets with boundary conditons, etc.
ierr = updateMeshsetsFieldsAndElements(ll+1); CHKERRQ(ierr);

Future work and contribution

If you like to contribute to this module (or MoFEM library), you are very welcome. Please feel free to add changes, including improving this documentation. If you like to work and modify this user module, you can fork this repository, make your changes, and do pull request. Contact us mofem-group@googlegroups.com or https://mofem.slack.com/ if you need any help or comments.

This is work in progress, and more stuff can be done. You could consider own paper or student project.

Make this documentation better
Problem tailored preconditioner for Krylov solver
Convergence study, estimation of constants