FaceRestriction numbering inconsistencies

[samuelpmishLLNL]

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:

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):

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:

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)?

[rw-anderson]

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

[samuelpmishLLNL]

@rw-anderson it's just a mathematica notebook: node_numbering.pdf

[YohannDudouit]

I'm not sure I'm following why you're saying that the FaceRestriction ordering is not lexicographic, as your animation is showing, it is the only one in lexicographic ordering (x, then y). GetBdrElementDofs + GetDofMap is not resulting in a lexicographic ordering (x, then -y).

[YohannDudouit]

Regarding the ordering of the faces, the only one that seems consistent with what is described in geom.cpp:990 is the one from FaceRestriction.

[samuelpmishLLNL]

I didn't say that the face restriction isn't lexicographic, I said that it doesn't appear to be consistent with mfem::FiniteElementSpace's in lexicographical ordering of the dofs. The point is that mfem is providing several ways to get (what should be) the same information and they don't agree.

From your comment it sounds like maybe you are thinking that the lexicographic ordering for a face means "the nodes should be ordered lexicographically w.r.t. (x, y, z) on their associated volume parent element". To me, it means that "the nodes of a face should be ordered lexicographically w.r.t. (xi, eta) in their parent element" (and this seems to be consistent with the FiniteElementSpace approach).

Perhaps the 6 faces are also reordered in the FaceRestriction in some lexicographical way (although I don't know what that would mean right now), but I don't believe this is desirable.

[YohannDudouit]

If you call "lexicographic" (xi,eta), then you need to apply permutations either on the interpolation operators or on the dofs, the purpose of the lexicographic ordering is to avoid that.

[samuelpmishLLNL]

I don't understand what you mean

[samuelpmishLLNL]

The mesh file appears to manually define the boundary faces in the following way: which means GetBdrElementDofs is numbering the faces (what I would guess is) correctly.

But it's still not clear why the FaceRestriction face and node orderings are not in agreement with mfem::FiniteElmentSpace::GetBdrElementDofs

[YohannDudouit]

The FaceRestriction is used to compute dofs on faces in order to compute values on faces at quadrature points, in order to do that one as to apply an interpolation operator, the interpolation operator expects dofs to be ordered lexicographically, if you use another ordering then you need to apply a permutation either on the interpolation operator or on the dofs. Since most basis functions are symmetric, we could exploit that these permutations are identity, but this would be a problem for non-symmetric basis functions. If you consider basis that are non-symmetric, e.g. RT or ND, then these permutations are not any more identity for instance. Finally, defining lexicographic as (x,y,z) is non-ambiguous, where (xi,eta) is ambiguous because the transformation is arbitrary.

[samuelpmishLLNL]

Finally, defining lexicographic as (x,y,z) is non-ambiguous, where (xi,eta) is ambiguous because the transformation is arbitrary.

There is no ambiguity in applying the lexicographical order to (xi, eta), provided that the faces have a clear ordering/winding (and they do, as you previously mentioned, here: https://github.com/mfem/mfem/blob/master/fem/geom.cpp#L990-L994)

[YohannDudouit]

There is ambiguity because it's arbitrary, you need to go read a hard-coded ordering that you need to reflect everywhere in your code, this is not the case with (x,y,z) ordering, you never need to go dive into this kind of considerations.

[samuelpmishLLNL]

There is ambiguity because it's arbitrary, you need to go read a hard-coded ordering that you need to reflect everywhere in your code, this is not the case with (x,y,z) ordering, you never need to go dive into this kind of considerations.

I agree that you need some hard coded ordering somewhere, but https://github.com/mfem/mfem/blob/master/fem/geom.cpp#L990-L994 is sufficient to define that ordering. But this discussion is veering away from the main point:

GetBdrElementDofs() together with GetDofMap() seem to work as I expected, and FaceRestriction is doing something else. Is this intentional or a bug?

[YohannDudouit]

I think I did my best to explain why it is intentional.

[samuelpmishLLNL]

To make matters even worse, after looking at the quadrature point data from GetFaceGeometricFactors:

It seems the quadrature point data is ordered in the same way as the FaceRestriction. So, my workaround plan of abandoning the use of mfem::FaceRestriction and using GetBdrElementDofs exclusively also doesn't work because the quadrature data is inconsistent.

[v-dobrev]

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:

1. The standard interface that supports all types of meshes. It works primarily at the level of individual face and boundary elements.
2. 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).

[samuelpmishLLNL]

(edited)

Thank you for the clarification.

I hope that, in the future, information like this will be made available directly in the relevant documentation to avoid unnecessary confusion.