Added unit test to check the presence of double faces in the space.

examples/ex3.cpp

@@ -106,19 +106,26 @@ int main(int argc, char *argv[])
    //    largest number that gives a final mesh with no more than 50,000
    //    elements.
    {
-      int ref_levels =
-         (int)floor(log(50000./mesh->GetNE())/log(2.)/dim);
+      int ref_levels = 1;
+         //(int)floor(log(50000./mesh->GetNE())/log(2.)/dim);
       for (int l = 0; l < ref_levels; l++)
       {
          mesh->UniformRefinement();
       }
    }
+#if 1
+   // boundary elements still need to be reoriented :-|
+   mesh->ReorientBdrElements();
+   mesh->EnsureNCMesh(true);
+#else
    mesh->ReorientTetMesh();
+#endif
 
    // 5. Define a finite element space on the mesh. Here we use the Nedelec
    //    finite elements of the specified order.
    FiniteElementCollection *fec = new ND_FECollection(order, dim);
    FiniteElementSpace *fespace = new FiniteElementSpace(mesh, fec);
+   cout << "Number of DOFs: " << fespace->GetVSize() << endl;
    cout << "Number of finite element unknowns: "
         << fespace->GetTrueVSize() << endl;
 
@@ -218,11 +225,42 @@ int main(int argc, char *argv[])
    // 15. Send the solution by socket to a GLVis server.
    if (visualization)
    {
+      // project to an L2^dim vector space
+      L2_FECollection d_fec(order, dim, 1);
+      FiniteElementSpace d_fes(mesh, &d_fec, 3);
+      GridFunction d_x(&d_fes);
+      x.ProjectVectorFieldOn(d_x);
+
       char vishost[] = "localhost";
       int  visport   = 19916;
       socketstream sol_sock(vishost, visport);
       sol_sock.precision(8);
-      sol_sock << "solution\n" << *mesh << x << flush;
+      sol_sock << "solution\n" << *mesh << d_x << flush;
+      sol_sock << "window_title 'solution'\n";
+
+      /*if (dim == 2)
+      {
+         sol_sock << "keys ARjlmevv\n";
+      }
+      else
+      {
+         sol_sock << "keys Amevv\n";
+      }
+      sol_sock << "pause\n" << flush;*/
+
+      // visualize basis functions
+      /*for (int i = 0; i < x.Size(); i++)
+      {
+         x = 0.0;
+         x(i) = 1.0;
+
+         stringstream title;
+         title << "basis " << i;
+
+         sol_sock << "solution\n" << *mesh << x;
+         sol_sock << "window_title '" << title.str() << "'\n";
+         sol_sock << "pause\n" << flush;
+      }*/
    }
 
    // 16. Free the used memory.

fem/fespace.cpp

@@ -768,6 +768,92 @@ int FiniteElementSpace::GetNumBorderDofs(Geometry::Type geom, int order) const
    return Geometry::NumVerts[geom] * (nv + ne);
 }
 
+void FiniteElementSpace::AddVarOrderDependencies(SparseMatrix &deps) const
+{
+   Array<int> master_dofs, slave_dofs;
+   IsoparametricTransformation T;
+   DenseMatrix I;
+
+   for (int entity = 1; entity < mesh->Dimension(); entity++)
+   {
+      const Table &ent_dofs = (entity == 1) ? var_edge_dofs : var_face_dofs;
+      int num_ent = (entity == 1) ? mesh->GetNEdges() : mesh->GetNFaces();
+      MFEM_ASSERT(ent_dofs.Size() == num_ent+1, "");
+
+      // add constraints within edges/faces holding multiple DOF sets
+      Geometry::Type last_geom = Geometry::INVALID;
+      for (int i = 0; i < num_ent; i++)
+      {
+         if (ent_dofs.RowSize(i) <= 1) { continue; }
+
+         Geometry::Type geom =
+            (entity == 1) ? Geometry::SEGMENT : mesh->GetFaceGeometry(i);
+
+         if (geom != last_geom)
+         {
+            T.SetIdentityTransformation(geom);
+            last_geom = geom;
+         }
+
+         // get lowest order variant DOFs and FE
+         int p = GetEntityDofs(entity, i, master_dofs, geom, 0);
+         const auto *master_fe = fec->GetFE(geom, p);
+
+         // constrain all higher order DOFs: interpolate lowest order function
+         for (int variant = 1; ; variant++)
+         {
+            int q = GetEntityDofs(entity, i, slave_dofs, geom, variant);
+            if (q < 0) { break; }
+
+            const auto *slave_fe = fec->GetFE(geom, q);
+            slave_fe->GetTransferMatrix(*master_fe, T, I);
+
+            AddDependencies(deps, master_dofs, slave_dofs, I);
+         }
+      }
+   }
+}
+
+static double tri_face_orient_pm[6][6] =
+{
+   {0,0, 1,0, 0,1}, {1,0, 0,0, 0,1}, {0,1, 0,0, 1,0},
+   {0,1, 1,0, 0,0}, {1,0, 0,1, 0,0}, {0,0, 0,1, 1,0}
+};
+
+void FiniteElementSpace::AddDoubleFaceDependencies(SparseMatrix &deps) const
+{
+   Array<int> master_dofs, slave_dofs;
+   IsoparametricTransformation T;
+   DenseMatrix I;
+
+   for (int i = 0; i < nd_double_faces.Size(); i++)
+   {
+      int face = nd_double_faces[i];
+      int inf1, inf2;
+      mesh->GetFaceInfos(face, &inf1, &inf2);
+
+      int order1, order2;
+      order1 = GetFaceDofs(face, master_dofs, 0); // elem1 side
+      order2 = GetFaceDofs(face, slave_dofs,  1); // elem2 side
+      MFEM_ASSERT(order1 == order2, "");
+
+      MFEM_ASSERT(mesh->GetFaceGeometry(face) == Geometry::TRIANGLE, "");
+      int nfdof = fec->GetNumDof(Geometry::TRIANGLE, order1);
+
+      int ori = inf2 % 64;
+      MFEM_ASSERT(ori >= 0 && ori < 6, "");
+      if (!i) { T.SetFE(&TriangleFE); }
+      T.SetPointMat(DenseMatrix(tri_face_orient_pm[ori], 2, 3));
+
+      auto *fe = fec->GetFE(Geometry::TRIANGLE, order1);
+      fe->GetLocalInterpolation(T, I);
+      // TODO: cache I matrices for (order,ori) pairs?
+
+      int nskip = master_dofs.Size() - nfdof;
+      AddDependencies(deps, master_dofs, slave_dofs, I, nskip);
+   }
+}
+
 int FiniteElementSpace::GetEntityDofs(int entity, int index, Array<int> &dofs,
                                       Geometry::Type master_geom,
                                       int variant) const
@@ -896,47 +982,15 @@ void FiniteElementSpace::BuildConformingInterpolation() const
       }
    }
 
-   // variable order spaces: enforce minimum rule on conforming edges/faces
    if (IsVariableOrder())
    {
-      for (int entity = 1; entity < mesh->Dimension(); entity++)
-      {
-         const Table &ent_dofs = (entity == 1) ? var_edge_dofs : var_face_dofs;
-         int num_ent = (entity == 1) ? mesh->GetNEdges() : mesh->GetNFaces();
-         MFEM_ASSERT(ent_dofs.Size() == num_ent+1, "");
-
-         // add constraints within edges/faces holding multiple DOF sets
-         Geometry::Type last_geom = Geometry::INVALID;
-         for (int i = 0; i < num_ent; i++)
-         {
-            if (ent_dofs.RowSize(i) <= 1) { continue; }
-
-            Geometry::Type geom =
-               (entity == 1) ? Geometry::SEGMENT : mesh->GetFaceGeometry(i);
-
-            if (geom != last_geom)
-            {
-               T.SetIdentityTransformation(geom);
-               last_geom = geom;
-            }
-
-            // get lowest order variant DOFs and FE
-            int p = GetEntityDofs(entity, i, master_dofs, geom, 0);
-            const auto *master_fe = fec->GetFE(geom, p);
-
-            // constrain all higher order DOFs: interpolate lowest order function
-            for (int variant = 1; ; variant++)
-            {
-               int q = GetEntityDofs(entity, i, slave_dofs, geom, variant);
-               if (q < 0) { break; }
-
-               const auto *slave_fe = fec->GetFE(geom, q);
-               slave_fe->GetTransferMatrix(*master_fe, T, I);
-
-               AddDependencies(deps, master_dofs, slave_dofs, I);
-            }
-         }
-      }
+      // variable order spaces: enforce minimum rule on conforming edges/faces
+      AddVarOrderDependencies(deps);
+   }
+   if (nd_double_faces.Size())
+   {
+      // Nédélec spaces with complex face orientations: constrain double faces
+      AddDoubleFaceDependencies(deps);
    }
 
    deps.Finalize();
@@ -1857,6 +1911,39 @@ void FiniteElementSpace::BuildNURBSFaceToDofTable() const
    face_dof = new Table(GetNF(), face_dof_list);
 }
 
+void FiniteElementSpace::GetDoubleFaces(Array<int> &double_faces) const
+{
+   double_faces.DeleteAll();
+
+   // TODO: explain
+
+   if (mesh->Dimension() == 3 &&
+       fec->GetContType() == FiniteElementCollection::TANGENTIAL &&
+       mesh->HasGeometry(Geometry::TRIANGLE) &&
+       GetMaxElementOrder() > 1)
+   {
+      for (int i = 0; i < mesh->GetNFaces(); i++)
+      {
+         if (mesh->GetFaceGeometry(i) == Geometry::TRIANGLE)
+         {
+            int elem1, elem2, inf1, inf2;
+            mesh->GetFaceElements(i, &elem1, &elem2);
+            mesh->GetFaceInfos(i, &inf1, &inf2);
+
+            int ori = inf2 % 64;
+            if (elem2 >= 0 && ori >= 1 && ori <= 4)
+            {
+               // TODO: check order
+               double_faces.Append(i);
+            }
+         }
+      }
+
+      // DEBUG
+      mfem::out << "### Number of double faces: " << double_faces.Size() << endl;
+   }
+}
+
 void FiniteElementSpace::Construct()
 {
    // This method should be used only for non-NURBS spaces.
@@ -1890,13 +1977,20 @@ void FiniteElementSpace::Construct()
    bool mixed_elements = (mesh->GetNumGeometries(dim) > 1);
    bool mixed_faces = (dim > 2 && mesh->GetNumGeometries(2) > 1);
 
+   // get a list of faces that need two sets of DOFs (Nedelec spaces only)
+   GetDoubleFaces(nd_double_faces);
+   MFEM_VERIFY(Nonconforming() || nd_double_faces.Size() == 0,
+               "H(curl) triangular faces of order >= 2 with orientations 1-4 are"
+               " only supported in NC meshes. Please use Mesh::ReorientTetMesh()"
+               " or Mesh::EnsureNCMesh().");
+
    Array<VarOrderBits> edge_orders, face_orders;
    if (IsVariableOrder())
    {
       // for variable order spaces, calculate orders of edges and faces
       CalcEdgeFaceVarOrders(edge_orders, face_orders);
    }
-   else if (mixed_faces)
+   else if (mixed_faces || nd_double_faces.Size())
    {
       // for mixed faces we also create the var_face_dofs table, see below
       face_orders.SetSize(mesh->GetNFaces());
@@ -1926,7 +2020,7 @@ void FiniteElementSpace::Construct()
    // assign face DOFs
    if (mesh->GetNFaces())
    {
-      if (IsVariableOrder() || mixed_faces)
+      if (IsVariableOrder() || mixed_faces || nd_double_faces.Size())
       {
          // NOTE: for simplicity, we also use Table var_face_dofs for mixed faces
          nfdofs = MakeDofTable(2, face_orders, var_face_dofs,
@@ -2106,6 +2200,10 @@ void FiniteElementSpace
    while (!done);
 }
 
+/// Return true iff 'n' is a power of 2.
+inline bool is_pow2(int n) { return !(n & (n-1)); }
+
+
 int FiniteElementSpace::MakeDofTable(int ent_dim,
                                      const Array<int> &entity_orders,
                                      Table &entity_dofs,
@@ -2134,11 +2232,11 @@ int FiniteElementSpace::MakeDofTable(int ent_dim,
    }
 
    // assign DOFs according to order bit masks
-   for (int i = 0; i < num_ent; i++)
+   for (int i = 0, j = 0; i < num_ent; i++)
    {
       auto geom = (ent_dim == 1) ? Geometry::SEGMENT : mesh->GetFaceGeometry(i);
 
-      VarOrderBits bits = entity_orders[i];
+      VarOrderBits bits = entity_orders[i], orig_bits = bits;
       for (int order = 0; bits != 0; order++, bits >>= 1)
       {
          if (bits & 1)
@@ -2147,6 +2245,21 @@ int FiniteElementSpace::MakeDofTable(int ent_dim,
             list.Append(Connection(i, total_dofs));
             total_dofs += dofs;
 
+            // add one more DOF set if 'i' is a double face with a single order
+            if ((ent_dim == 2) &&
+                (j < nd_double_faces.Size()) &&
+                (i == nd_double_faces[j]))
+            {
+               if (is_pow2(orig_bits))
+               {
+                  list.Append(Connection(i, total_dofs));
+                  total_dofs += dofs;
+                  if (var_ent_order) { var_ent_order->Append(order); }
+               }
+               j++;
+            }
+
+            // record the order of the DOF set variant
             if (var_ent_order) { var_ent_order->Append(order); }
          }
       }
@@ -2214,6 +2327,17 @@ int FiniteElementSpace::GetNVariants(int entity, int index) const
    return dof_table.GetRow(index + 1) - dof_table.GetRow(index);
 }
 
+bool FiniteElementSpace::IsDoubleFace(int face) const
+{
+   if (!nd_double_faces.Size()) { return false; }
+
+   const int size = var_face_dofs.RowSize(face);
+   const int *row = var_face_dofs.GetRow(face);
+
+   // a double face has exactly two DOF sets of the same size
+   return (size == 2) && ((row[1] - row[0]) == (row[2] - row[1]));
+}
+
 static const char* msg_orders_changed =
    "Element orders changed, you need to Update() the space first.";
 
@@ -2283,10 +2407,26 @@ void FiniteElementSpace::GetElementDofs(int elem, Array<int> &dofs) const
       for (int i = 0; i < F.Size(); i++)
       {
          auto fgeom = mesh->GetFaceGeometry(F[i]);
-         int nf = fec->GetNumDof(fgeom, order);
+         int fbase, nf = fec->GetNumDof(fgeom, order);
+         const int *ind;
 
-         int fbase = (var_face_dofs.Size() > 0) ? FindFaceDof(F[i], nf) : F[i]*nf;
-         const int *ind = fec->GetDofOrdering(fgeom, order, Fo[i]);
+         if (var_face_dofs.Size() <= 0) // simple constant-order space
+         {
+            fbase = F[i]*nf;
+            ind = fec->GetDofOrdering(fgeom, order, Fo[i]);
+         }
+         else if (IsDoubleFace(F[i])) // special disconnected face
+         {
+            int e1, e2;
+            mesh->GetFaceElements(F[i], &e1, &e2);
+            fbase = var_face_dofs.GetRow(F[i])[(elem == e1) ? 0 : 1];
+            ind = fec->GetDofOrdering(fgeom, order, 0);
+         }
+         else // mixed geometry or variable-order faces
+         {
+            fbase = FindFaceDof(F[i], nf);
+            ind = fec->GetDofOrdering(fgeom, order, Fo[i]);
+         }
 
          for (int j = 0; j < nf; j++)
          {
@@ -2313,8 +2453,12 @@ const FiniteElement *FiniteElementSpace::GetFE(int i) const
    MFEM_VERIFY(i < mesh->GetNE(),
                "Invalid element id " << i << ", maximum allowed " << mesh->GetNE()-1);
 
-   const FiniteElement *FE =
-      fec->GetFE(mesh->GetElementGeometry(i), GetElementOrderImpl(i));
+   Geometry::Type geom = mesh->GetElementGeometry(i);
+   const FiniteElement *FE = fec->GetFE(geom, GetElementOrderImpl(i));
+
+   MFEM_VERIFY(FE != NULL, "FiniteElementCollection " << fec->Name()
+               << " does not support element geometry " << geom
+               << " (" << Geometry::Name[geom] << ").");
 
    if (NURBSext)
    {
@@ -2427,7 +2571,21 @@ int FiniteElementSpace::GetFaceDofs(int face, Array<int> &dofs,
    int dim = mesh->Dimension();
    auto fgeom = (dim > 2) ? mesh->GetFaceGeometry(face) : Geometry::INVALID;
 
-   if (var_face_dofs.Size() > 0) // variable orders or *mixed* faces
+   if (var_face_dofs.Size() <= 0) // simple constant-order space
+   {
+      if (variant > 0) { return -1; }
+      order = fec->GetOrder();
+      nf = (dim > 2) ? fec->GetNumDof(fgeom, order) : 0;
+      fbase = face*nf;
+   }
+   else if (IsDoubleFace(face)) // special disconnected face
+   {
+      MFEM_ASSERT(variant >= 0 && variant < 2, "");
+      fbase = var_face_dofs.GetRow(face)[variant];
+      order = GetFaceOrder(face);
+      nf = fec->GetNumDof(fgeom, order);
+   }
+   else // variable-order or *mixed* faces
    {
       const int* beg = var_face_dofs.GetRow(face);
       const int* end = var_face_dofs.GetRow(face + 1);
@@ -2435,18 +2593,9 @@ int FiniteElementSpace::GetFaceDofs(int face, Array<int> &dofs,
 
       fbase = beg[variant];
       nf = beg[variant+1] - fbase;
-
-      order = !IsVariableOrder() ? fec->GetOrder() :
-              var_face_orders[var_face_dofs.GetI()[face] + variant];
+      order = var_face_orders[var_face_dofs.GetI()[face] + variant];
       MFEM_ASSERT(fec->GetNumDof(fgeom, order) == nf, "");
    }
-   else
-   {
-      if (variant > 0) { return -1; }
-      order = fec->GetOrder();
-      nf = (dim > 2) ? fec->GetNumDof(fgeom, order) : 0;
-      fbase = face*nf;
-   }
 
    // for 1D, 2D and 3D faces
    int nv = fec->GetNumDof(Geometry::POINT, order);

fem/fespace.hpp

@@ -118,12 +118,17 @@ protected:
    /** Variable order spaces only: DOF assignments for edges and faces, see
        docs in MakeDofTable. For constant order spaces the tables are empty. */
    Table var_edge_dofs;
-   Table var_face_dofs; ///< NOTE: also used for spaces with mixed faces
+   Table var_face_dofs; ///< NOTE: also used for spaces with mixed/double faces
 
    /** Additional data for the var_*_dofs tables: individual variant orders
        (these are basically alternate J arrays for var_edge/face_dofs). */
    Array<char> var_edge_orders, var_face_orders;
 
+   /** List of faces with two sets of DOFs. This is a special feature for
+       triangular Nedelec faces with orientations 1-4 and order >= 2, that we
+       disconnect and constrain with the P matrix. Normally the list is empty. */
+   Array<int> nd_double_faces;
+
    // precalculated DOFs for each element, boundary element, and face
    mutable Table *elem_dof; // owned (except in NURBS FE space)
    mutable Table *bdr_elem_dof; // owned (except in NURBS FE space)
@@ -272,9 +277,16 @@ protected:
                                 Array<int> &slave_dofs, int slave_face,
                                 const DenseMatrix *pm) const;
 
+   void AddVarOrderDependencies(SparseMatrix &deps) const;
+
+   void AddDoubleFaceDependencies(SparseMatrix &deps) const;
+
    /// Replicate 'mat' in the vector dimension, according to vdim ordering mode.
    void MakeVDimMatrix(SparseMatrix &mat) const;
 
+   void GetDoubleFaces(Array<int> &double_faces) const;
+   bool IsDoubleFace(int face) const;
+
    /// GridFunction interpolation operator applicable after mesh refinement.
    class RefinementOperator : public Operator
    {

mesh/mesh.cpp

@@ -5860,6 +5860,23 @@ void Mesh::ReorientTetMesh()
    }
 }
 
+void Mesh::ReorientBdrElements()
+{
+   // TODO: maybe there is a safer way to do this
+   for (int i = 0; i < boundary.Size(); i++)
+   {
+      int *bv = boundary[i]->GetVertices();
+      int *fv = faces[be_to_face[i]]->GetVertices();
+
+      // make sure boundary element ordering mirrors that of faces
+      int nv = boundary[i]->GetNVertices();
+      for (int j = 0; j < nv; j++)
+      {
+         bv[j] = fv[j];
+      }
+   }
+}
+
 int *Mesh::CartesianPartitioning(int nxyz[])
 {
    int *partitioning;

mesh/mesh.hpp

@@ -1128,6 +1128,9 @@ public:
        @note Refinement does not work after a call to this method! */
    virtual void ReorientTetMesh();
 
+   /** TODO */
+   virtual void ReorientBdrElements();
+
    int *CartesianPartitioning(int nxyz[]);
    int *GeneratePartitioning(int nparts, int part_method = 1);
    void CheckPartitioning(int *partitioning);

mesh/ncmesh.cpp

@@ -1328,7 +1328,7 @@ void NCMesh::RefineElement(int elem, char ref_type)
       //     +                         Faces: 0 back (1, 2, 3)
       //     |\\_                             1 left (0, 3, 2)
       //     ||  \_                           2 front (0, 1, 3)
-      //     | \   \_                         3 bottom (0, 1, 2)
+      //     | \   \_                         3 bottom (0, 2, 1)
       //     |  +__  \_
       //     | /2  \__ \_       Z  Y
       //     |/       \__\      | /

tests/unit/fem/test_double_faces.cpp

@@ -0,0 +1,67 @@
+// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "mfem.hpp"
+#include "unit_tests.hpp"
+
+namespace mfem
+{
+
+// Check complex orientation constraints for Nedelec elements
+// of order >= 2 on triangular faces with orientations 1-4.
+
+TEST_CASE("Nedelec arbitrary orientation tri faces",
+          "[FiniteElementSpace]"
+          "[NCMesh]")
+{
+   const int dim = 3;
+   const int order = GENERATE(1, 2, 3);
+
+   Mesh mesh(2, 2, 2, Element::TETRAHEDRON);
+   mesh.EnsureNCMesh(true);
+   // NOTE: no mesh.ReorientTetMesh()
+
+   int ncomplex = 0;
+   for (int i = 0; i < mesh.GetNumFaces(); i++)
+   {
+      int elem1, elem2, inf1, inf2;
+      mesh.GetFaceElements(i, &elem1, &elem2);
+      mesh.GetFaceInfos(i, &inf1, &inf2);
+
+      int ori = inf2 % 64;
+      if (elem2 >= 0 && ori >= 1 && ori <= 4)
+      {
+         //mfem::out << "Face orientation " << ori << " found\n";
+         ncomplex++;
+      }
+   }
+
+   // check that there are triangular faces with orientations 1-4 in the mesh
+   // -- if there aren't any it's not an error per se, but the following
+   // test won't work...
+   REQUIRE(ncomplex > 0);
+
+   ND_FECollection fec(order, dim);
+   FiniteElementSpace fespace(&mesh, &fec);
+
+   if (order < 2)
+   {
+      // no need for double faces at order 1
+      REQUIRE(fespace.GetVSize() == fespace.GetTrueVSize());
+   }
+   else
+   {
+      // there are extra DOFs in the space due to the double faces
+      REQUIRE(fespace.GetVSize() > fespace.GetTrueVSize());
+   }
+}
+
+} // namespace mfem