mfem: FaceRestriction numbering inconsistencies

I’ve been trying to debug this boundary dof issue for days now, and I need some help making sense of what I’m seeing.

I’m using mfem::FiniteElementSpace::GetFaceRestriction to get an operator that would take me from an “L-vector” to a boundary “E-vector”. Our shape functions are ordered lexicographically, so I call it with the following arguments:

test_space_->GetFaceRestriction(mfem::ElementDofOrdering::LEXICOGRAPHIC, mfem::FaceType::Boundary, mfem::L2FaceValues::SingleValued)

I need the dof ids as well if I’m going to add the boundary-element jacobian (e.g. stiffness) matrix contributions back into their appropriate locations in the global jacobian matrix. I don’t see an obvious way to get that information out of the restriction operator itself, so I looked to mfem::FiniteElementSpace::GetBdrElementDofs(). It appears that the output of GetBdrElementDofs is in NATIVE ordering, so I apply the permutation from calling GetDofMap() to get the dofs in LEXICOGRAPHIC ordering. When I compare this with the dofs used by FaceRestriction I see the following:

boundary_element_node_numbering

the FaceRestriction’s node numbering within a face appears to be inconsistent with both the NATIVE and LEXICOGRAPHIC options. In addition, the numbering of the faces in the FaceRestriction is also inconsistent with GetBdrElementDofs (face number indicated in parentheses):

boundary_element_face_numbering

Also, the windings of the FaceRestriction node numbering are also inconsistent: for faces 1, 4, and 5, the normal directions (n := cross(dx_dxi, dx_deta)) are inward-facing while faces 2, 3, and 6 have outward-facing normals. As a sanity check, I compared GetElementDofs against GetElementRestriction’s numbering as well:

element_node_numbering

So, there seems to be agreement between GetElementDofs and GetElementRestriction, but not GetFaceRestriction and GetBdrElementDofs.

Is this behavior a bug, or is it intentional (and if so, is there any documentation or explanation available)?

About this issue

Original URL
State: closed
Created 3 years ago
Comments: 16 (16 by maintainers)

Most upvoted comments

Let me try to summarize the functionality in MFEM that has to do with faces and boundaries. There are basically two separate interfaces for this:

The standard interface that supports all types of meshes. It works primarily at the level of individual face and boundary elements.
The new interface that currently supports only meshes with tensor product elements. This interface works on subsets of faces: either the set of all interior faces (faces that have 2 adjacent volume elements) or the set of all boundary faces (faces that have only 1 adjacent volume element).

In the context of (1) there are two sets of elements of dimension (dim-1), where dim = Mesh::Dimension(): (A) the set of all face elements which is auto-generated from the volume elements, and (B) the set of all “boundary” elements which is typically read from the mesh file. Note that the boundary elements (B) contain the boundary attributes. Also, note that the boundary elements can describe actual interior faces (interfaces). Every element in (B) can be identified uniquely as a face in set (A), however, the boundary element may have an orientation (vertex ordering) different from that of its corresponding face in (A). The number of elements in (A) and (B) are given by Mesh::GetNumFaces() and Mesh::GetNBE(), respectively.

In the context of (2), the set (A) is decomposed in two parts: (Ai) the set of all real interior faces (faces that have 2 adjacent volume elements), enumerated contiguously following the ordering from (A); and (Ab) the set of all real boundary faces (faces that have 1 adjacent volume element), also enumerated contiguously following the ordering from (A). The number of elements in (Ai) and (Ab) are given by Mesh::GetNFbyType(FaceType::Interior) and Mesh::GetNFbyType(FaceType::Boundary), respectively. An important point is that in the context (2), both the elements in (Ai) and those in (Ab) have orientations (vertex ordering) different from the one in (A). The main reason for using different orientations for (Ai) and (Ab) is to simplify the implementation of the tensor product evaluations.

Now, given the above description, we have the following:

The methods FiniteElementSpace::GetFaceRestriction, Mesh::GetFaceGeometricFactors, and FiniteElementSpace::GetFaceQuadratureInterpolator are part of (2) and they work on the sets (Ai) or (Ab) depending on their parameter of type FaceType.
Methods like Mesh::GetBdrElementVertices, Mesh::GetBdrElementTransformation, and FiniteElementSpace::GetBdrElementDofs are part of (1) and they work with the set (B).

Hopefully this helps explain why you see the inconsistencies you observed.

I think there are different paths forward: either using (1), or using (2) with some functionality that we can add in MFEM, e.g. a method similar to FiniteElementSpace::GetBdrElementDofs that works with the set (Ab) instead of (B).

v-dobrev on Sep 2, 2021

I have no idea about these details, I just want to know how you made these animations. 😃

rw-anderson on Aug 30, 2021