Revert changes in Mesh::GetLocalTriToTetTransformation() and

Mesh::GetLocalQuadToHexTransformation(). To address the issue seen when calling Mesh::GetBdrElementAdjacentElement() followed by calling Mesh::GetLocalFaceTransformation(), introduce a new method, Mesh::GetBdrElementAdjacentElement2() that returns the orientation of the face element w.r.t. the boundary element. This orientation then produces the desired result when used with Mesh::GetLocalFaceTransformation().
NorbertHofbauer · Sep 3, 2022 · 5180b3f · 5180b3f
1 parent f7fcfe1
commit 5180b3f
Show file tree

Hide file tree

Showing 3 changed files with 51 additions and 14 deletions.
diff --git a/mesh/mesh.cpp b/mesh/mesh.cpp
@@ -470,7 +470,7 @@ void Mesh::GetBdrElementTransformation(int i, IsoparametricTransformation* ElTr)
       else // L2 Nodes (e.g., periodic mesh)
       {
          int elem_id, face_info;
-         GetBdrElementAdjacentElement(i, elem_id, face_info);
+         GetBdrElementAdjacentElement2(i, elem_id, face_info);
 
          GetLocalFaceTransformation(GetBdrElementType(i),
                                     GetElementType(elem_id),
@@ -710,11 +710,10 @@ void Mesh::GetLocalTriToTetTransformation(
    locpm.SetSize(3, 3);
    for (int j = 0; j < 3; j++)
    {
-      int k = to[j];
-      const IntegrationPoint &vert = TetVert->IntPoint(tv[j]);
-      locpm(0, k) = vert.x;
-      locpm(1, k) = vert.y;
-      locpm(2, k) = vert.z;
+      const IntegrationPoint &vert = TetVert->IntPoint(tv[to[j]]);
+      locpm(0, j) = vert.x;
+      locpm(1, j) = vert.y;
+      locpm(2, j) = vert.z;
    }
 }
 
@@ -784,11 +783,10 @@ void Mesh::GetLocalQuadToHexTransformation(
    locpm.SetSize(3, 4);
    for (int j = 0; j < 4; j++)
    {
-      int k = qo[j];
-      const IntegrationPoint &vert = HexVert->IntPoint(hv[j]);
-      locpm(0, k) = vert.x;
-      locpm(1, k) = vert.y;
-      locpm(2, k) = vert.z;
+      const IntegrationPoint &vert = HexVert->IntPoint(hv[qo[j]]);
+      locpm(0, j) = vert.x;
+      locpm(1, j) = vert.y;
+      locpm(2, j) = vert.z;
    }
 }
 
@@ -6227,6 +6225,28 @@ void Mesh::GetBdrElementAdjacentElement(int bdr_el, int &el, int &info) const
    info = fi.Elem1Inf + ori;
 }
 
+void Mesh::GetBdrElementAdjacentElement2(int bdr_el, int &el, int &info) const
+{
+   int fid = GetBdrElementEdgeIndex(bdr_el);
+
+   const FaceInfo &fi = faces_info[fid];
+   MFEM_ASSERT(fi.Elem1Inf % 64 == 0, "internal error"); // orientation == 0
+
+   const int *fv = (Dim > 1) ? faces[fid]->GetVertices() : NULL;
+   const int *bv = boundary[bdr_el]->GetVertices();
+   int ori;
+   switch (GetBdrElementGeometry(bdr_el))
+   {
+      case Geometry::POINT:    ori = 0; break;
+      case Geometry::SEGMENT:  ori = (fv[0] == bv[0]) ? 0 : 1; break;
+      case Geometry::TRIANGLE: ori = GetTriOrientation(bv, fv); break;
+      case Geometry::SQUARE:   ori = GetQuadOrientation(bv, fv); break;
+      default: MFEM_ABORT("boundary element type not implemented"); ori = 0;
+   }
+   el   = fi.Elem1No;
+   info = fi.Elem1Inf + ori;
+}
+
 Element::Type Mesh::GetElementType(int i) const
 {
    return elements[i]->GetType();

diff --git a/mesh/mesh.hpp b/mesh/mesh.hpp
@@ -1151,9 +1151,25 @@ class Mesh
    int GetBdrElementEdgeIndex(int i) const;
 
    /** @brief For the given boundary element, bdr_el, return its adjacent
-       element and its info, i.e. 64*local_bdr_index+bdr_orientation. */
+       element and its info, i.e. 64*local_bdr_index+bdr_orientation.
+
+       The returned bdr_orientation is that of the boundary element relative to
+       the respective face element.
+
+       @sa GetBdrElementAdjacentElement2() */
    void GetBdrElementAdjacentElement(int bdr_el, int &el, int &info) const;
 
+   /** @brief For the given boundary element, bdr_el, return its adjacent
+       element and its info, i.e. 64*local_bdr_index+inverse_bdr_orientation.
+
+       The returned inverse_bdr_orientation is the inverse of the orientation of
+       the boundary element relative to the respective face element. In other
+       words this is the orientation of the face element relative to the
+       boundary element.
+
+       @sa GetBdrElementAdjacentElement() */
+   void GetBdrElementAdjacentElement2(int bdr_el, int &el, int &info) const;
+
    /// Returns the type of element i.
    Element::Type GetElementType(int i) const;
 

diff --git a/mesh/submesh/submesh_utils.cpp b/mesh/submesh/submesh_utils.cpp
@@ -116,8 +116,9 @@ void BuildVdofToVdofMap(const FiniteElementSpace& subfes,
             auto pm = parentfes.GetMesh();
 
             int face_info, parent_volel_id;
-            pm->GetBdrElementAdjacentElement(parent_element_ids[i], parent_volel_id,
-                                             face_info);
+            pm->GetBdrElementAdjacentElement2(parent_element_ids[i],
+                                              parent_volel_id,
+                                              face_info);
             pm->GetLocalFaceTransformation(
                pm->GetBdrElementType(parent_element_ids[i]),
                pm->GetElementType(parent_volel_id),