HowTo: Nonconforming and AMR meshes
The Mesh
class provides basic element refinement capabilities:
 All elements may be refined uniformly with
Mesh::UniformRefinement
.  Local refinement is supported, but only for simplex elements. The method
Mesh::GeneralRefinement
uses recursive bisection in this case.
These basic refinement methods preserve mesh conformity, i.e., no hanging nodes
are created. This also means that quadrilaterals and hexahedra cannot be refined
locally by the Mesh
class.
For more advanced AMR, MFEM has the class NCMesh
:
 Tensor product element refinement (quad, hex, prism) is supported, including anisotropic refinement. Hanging nodes are created and handled transparently.
 Triangles and tetrahedra use "red" (isotropic) refinement, also producing hanging nodes in this mode.
 Derefinement (coarsening) of previously refined elements is possible.
 In parallel, the mesh can be load balanced.
The user does not interact directly with the NCMesh
class — it is
created behind the scenes, and the Mesh
class in nonconforming mode, continually updated
to contain the finest elements of the refinement hierarchy, still serves as an interface
for the user and other MFEM classes.
To switch to the nonconforming mode (or convert and existing conforming Mesh
),
you need to call EnsureNCMesh
, typically at the beginning after loading the
mesh:
Mesh *mesh = new Mesh(mesh_file, 1, 1);
mesh>EnsureNCMesh(true);
The boolean parameter, if true
, forces simplex meshes to use nonconforming
refinement (the default is false
).
Nonconforming refinement
Once the Mesh
is in nonconforming mode, you can simply call Mesh::GeneralRefinement
to locally refine a subset of elements:
Array<int> refinement_list;
for (int i = 0; i < mesh>GetNE(); i++)
{
if (/*element i refinement condition*/)
{
refinement_list.Append(i);
}
}
mesh>GeneralRefinement(refinement_list);
The resulting hanging nodes will be treated transparently by the FiniteElementSpace
and
BilinearForm
classes:

FiniteElementSpace
will internally construct a conforming interpolation matrix $P$, that when applied to a vector of unconstrained ("true") DOFs, will augment the vector with interpolated constrained DOFs. 
Once the linear system $Ax = b$ is assembled,
BilinearForm::FormLinearSystem
will eliminate constrained nodes by transforming the linear system to $P^TAPx = P^Tb$ (seeex1.cpp
). 
After the reduced system is solved, the conforming solution on all nodes is recovered as $y = P x$ with
BilinearForm::RecoverFEMSolution
.
Limiting the level of hanging nodes
By default, MFEM does not limit the sizes of adjacent elements in nonconforming meshes. For some applications, it may be necessary to ensure that the refinement level of neighboring elements differs by at most one, for example.
The optional parameter nc_limit
of Mesh::GeneralRefinement
can be used to
control the maximum level of nonconformity. If nc_limit
is greater than zero, the
method will automatically perform additional refinements to make sure the difference of
refinement levels of adjacent elements is at most nc_limit
.
Anisotropic refinement
Uniquely, MFEM offers the capability to perform anisotropic refinement of tensor product
elements in both 2D and 3D. The method Mesh::GeneralRefinement
has two overloads, one
taking a simple list of elements to refine (as seen above), and the other taking a list
of struct Refinement { int index; char ref_type; }
, where one can specify a refinement
type for each element in the list:
Array<Refinement> refinement_list;
refinement_list.Append(Refinement(0, 2));
refinement_list.Append(Refinement(1, 4));
mesh>GeneralRefinement(refinement_list);
This code will refine the first element (index 0) of the mesh in the Y direction only (provided it is a quad or hex element) and the second element (index 1) in the Z direction only. The directions are assumed in the element reference coordinates and are encoded as follows:
Note that the refinement type is encoded as a 3bit number, where bits 0, 1, 2 correspond to the X, Y, Z directions, respectively. Other element geometries allow fewer but similar refinement types: triangle (3), quadrilateral (1, 2, 3), tetrahedron (7), prism (3, 4, 7).
In 3D meshes with anisotropic refinements it is easy to arrive at conflicting situations,
where the refined faces of adjacent elements are not subsets of each other. For example,
running the above code on a mesh with two hexahedra adjacent in the X direction will
create an interface that cannot be constrained correctly. In such cases, MFEM will
automatically adjust one side of the interface with additional refinements (called forced
refinements) to ensure that the mesh remains a valid FEM mesh. In pathological cases the
forced refinements may propagate. Using a reasonable nc_limit
may reduce this effect.
Nevertheless, a valid mesh is produced in all cases.
Derefinement
To coarsen elements, use the method Mesh::DerefineByError
.
The interface is different from refinement, because it is not possible to coarsen
arbitrary groups of fine elements: it is only possible to reintroduce previously existing
coarse elements by undoing their refinement (hence the term "derefinement").
Since one cannot supply the indices of elements that no longer exist in the Mesh
class
(the refinement trees are kept internal to NCMesh
), the method
DerefineByError
works indirectly by taking an array of "error" values corresponding to
each element of the current Mesh
. If the sum of error values of the children of some
coarse element is below a supplied threshold, the children are removed and the coarse
element is restored in Mesh
.
If the user specifies a nonzero nc_limit
, care is taken not to derefine elements that
are needed to keep the required level of nonconformity.
Note: derefinement is not yet supported for meshes containing 3D anisotropic refinements.
Parallel nonconforming meshes
Just as the Mesh
class has a parallel counterpart (ParMesh
), so does the NCMesh
class have a parallel descendant: ParNCMesh
. The parallel class is again kept
internal and the user can continue to interact with the standard ParMesh
class
(see examples ex1p
, ex6p
and ex15p
).
The refinement hierarchy in parallel NC mode is fully distributed and scales to billions
of elements and hundreds of thousands of MPI tasks. Ghost elements are automatically tracked
by the ParNCMesh
class, so that a parallel conforming interpolation matrix can be constructed
by ParFiniteElementSpace
. Depending on the assembly level, ParBilinearForm
will either explicitly assemble the parallel $P^TAP$ system using the Hypre library, or
the action of the $P$ matrix will be applied during solver iterations.
Parallel refinement is still done through Mesh::GeneralRefinement
inherited by the ParMesh
class. The method takes local element indices and works the same as in serial. All parallel
concerns such as keeping the ghost layers synchronized are handled internally in ParNCMesh
.
Note: parallel anisotropic refinement of 3D meshes is not supported yet.
After each mesh operation (refinement, derefinement, load balancing) the ParMesh
is
updated to reflect the current parallel mesh state (minus the ghost elements, which are not exported to ParMesh
). Communication groups, used in conforming mode for reductions/broadcasts over parallel solution vectors, are approximated in the NC mode as if the mesh was cut along the nonconforming interfaces.
Load balancing
In conforming mode, a serial Mesh
can only be partitioned statically (with METIS) when
constructing a ParMesh
. In nonconforming mode, the internal ParNCMesh
class is capable
of load balancing the distributed mesh at any time. This functionality is available to the
user through ParMesh::Rebalance
(see ex6p
and ex15p
).
The dynamic load balancing algorithm is based on partitioning a spacefilling curve (SFC) that naturally arises when traversing the distributed refinement trees. Compared to spectral partitioners like METIS the partitions are not as high quality but the process is extremely fast and scales to hundreds of thousands of processors.
For best results with SFCbased partitioning, one condition has to be met: the elements
of the coarse Mesh
from which the ParMesh
is constructed need to be
ordered, ideally as a sequence of faceneighbors. This makes it possible for ParNCMesh
to order the leaves of all refinement trees into a global linear sequence,
which when equipartitioned should produce compact (albeit not minimal surface) mesh partitions.
Take for example a coarse mesh produced by the polarnc
miniapp. Except for two
discontinuities, the elements are mostly ordered as a sequence of faceneighbors:
When we start refining elements (in both serial and parallel), MFEM will try to keep the
spacefilling curve continuous by inserting local Hilbert curves in the refined areas (press
Ctrl+O
in GLVis to visualize the ordering curve):
In a parallel computation, the global curve is then used for fast assignment of elements
to MPI ranks. In the following run of ex15p
, each processor is assigned the same number
of elements (+/ one element). Note that the last partition is discontinuous due to a jump
in ordering in the coarse mesh. This only affects the efficiency of MPI communication —
the numerical results will be the same regardless of the partitioning.
MFEM provides several methods to help with mesh ordering:

Procedurally generated rectangular grids (
Mesh::MakeCartesian2D
,Mesh::MakeCartesian3D
and alsoMFEM INLINE mesh v1.0
files) are by default ordered along a pseudoHilbert curve. Note that even grid dimensions are recommended, as explained here. 
General unstructured meshes may be ordered by a spatial sort algorithm (
Mesh::GetHilbertElementOrdering
). This is a fast method that will leave a number of jumps in complex meshes, but it is still highly recommended over not ordering the mesh at all. 
Highquality orderings of general meshes can be obtained with the Gecko library, now included directly in MFEM and available as
Mesh::GetGeckoElementOrdering
. The optimization algorithm used is more costly than a simple spatial sort, but it should produce better orderings for meshes with complex geometries. Beware the exponential cost of increasing thewindow
parameter. Large meshes should probably be ordered in a preprocessing step (you may use themeshexplorer
miniapp for that).
Nonconforming mesh I/O
Nonconforming meshes have their own file format MFEM NC mesh v1.0
, which supports all the additional internal structures (refinement trees, hanging nodes, etc.) and works for both serial and parallel NC meshes.
The method ParMesh::ParPrint
will automatically choose the right format and can be used to
save and restart an AMR computation, as demonstrated in example ex6p
.
ParMesh::ParPrint
should not be confused with the method ParMesh::Print
, an analog of
Mesh::Print
, which is only suitable for visualization, as it uses the serial MFEM mesh v1.0
format and only adds the parallel shared faces to the output.