Merge remote-tracking branch 'origin/hypre-runtime-compute-policy' into apps-dev

.gitignore

@@ -266,7 +266,6 @@ miniapps/navier/*_output
 miniapps/nurbs/nurbs_ex1
 miniapps/nurbs/nurbs_ex1p
 miniapps/nurbs/nurbs_ex11p
-miniapps/nurbs/nurbs_printfunc
 miniapps/nurbs/nurbs_patch_ex1
 miniapps/nurbs/nurbs_curveint
 miniapps/nurbs/refined.mesh

CHANGELOG

@@ -11,15 +11,6 @@
 Version 4.6.1 (development)
 ===========================
 
-Discretization improvements
----------------------------
-- Introduced support for higher order non conformal Nedelec elements on
-  simplices in ParMesh.
-
-Miscellaneous
--------------
-- The ReadCubit Genesis mesh importer has been rewritten to improve readability.
-
 
 Version 4.6, released on September 27, 2023
 ===========================================
@@ -40,7 +31,6 @@ Meshing improvements
   * The edge to knot map for NURBS meshes can be determined automatically. It is
     no longer needed to specify this in the NURBS mesh.
   * Added curve interpolation method for NURBS.
-  * Added new small miniapp for printing of shape functions of a KnotVector
   * See miniapps/nurbs for example meshes and miniapps.
 
 Discretization improvements

config/defaults.mk

@@ -331,30 +331,16 @@ STRUMPACK_OPT = -I$(STRUMPACK_DIR)/include $(SCOTCH_OPT)
 STRUMPACK_LIB = -L$(STRUMPACK_DIR)/lib -lstrumpack $(MPI_FORTRAN_LIB)\
  $(SCOTCH_LIB) $(SCALAPACK_LIB)
 
-# Ginkgo library configuration
+# Ginkgo library configuration (currently not needed)
 GINKGO_DIR = @MFEM_DIR@/../ginkgo/install
-GINKGO_SEARCH_DIR = $(subst @MFEM_DIR@,$(MFEM_DIR),$(GINKGO_DIR))
 GINKGO_BUILD_TYPE=Release
 ifeq ($(MFEM_USE_GINKGO),YES)
    BASE_FLAGS = -std=c++14
 endif
 GINKGO_OPT = -isystem $(GINKGO_DIR)/include
-GINKGO_LIB_DIR = $(sort $(dir $(wildcard\
- $(GINKGO_SEARCH_DIR)/lib*/libginkgo*.a\
- $(GINKGO_SEARCH_DIR)/lib*/libginkgo*.so\
- $(GINKGO_SEARCH_DIR)/lib*/libginkgo*.dylib\
- $(GINKGO_SEARCH_DIR)/lib*/libginkgo*.dll)))
-GINKGO_LINK_LIB_DIR = $(GINKGO_DIR)$(subst $(GINKGO_SEARCH_DIR),,$(GINKGO_LIB_DIR))
-ALL_GINKGO_LIBS_DEBUG = $(notdir $(basename $(wildcard\
- $(GINKGO_SEARCH_DIR)/lib*/libginkgo*d.a\
- $(GINKGO_SEARCH_DIR)/lib*/libginkgo*d.so\
- $(GINKGO_SEARCH_DIR)/lib*/libginkgo*d.dylib\
- $(GINKGO_SEARCH_DIR)/lib*/libginkgo*d.dll)))
-ALL_GINKGO_LIBS = $(notdir $(basename $(wildcard\
- $(GINKGO_SEARCH_DIR)/lib*/libginkgo*.a\
- $(GINKGO_SEARCH_DIR)/lib*/libginkgo*.so\
- $(GINKGO_SEARCH_DIR)/lib*/libginkgo*.dylib\
- $(GINKGO_SEARCH_DIR)/lib*/libginkgo*.dll)))
+GINKGO_LIB_DIR = $(sort $(dir $(wildcard $(GINKGO_DIR)/lib*/libginkgo*.a  $(GINKGO_DIR)/lib*/libginkgo*.so $(GINKGO_DIR)/lib*/libginkgo*.dylib $(GINKGO_DIR)/lib*/libginkgo*.dll)))
+ALL_GINKGO_LIBS_DEBUG = $(notdir $(basename $(wildcard $(GINKGO_DIR)/lib*/libginkgo*d.a  $(GINKGO_DIR)/lib*/libginkgo*d.so $(GINKGO_DIR)/lib*/libginkgo*d.dylib $(GINKGO_DIR)/lib*/libginkgo*d.dll)))
+ALL_GINKGO_LIBS = $(notdir $(basename $(wildcard $(GINKGO_DIR)/lib*/libginkgo*.a  $(GINKGO_DIR)/lib*/libginkgo*.so $(GINKGO_DIR)/lib*/libginkgo*.dylib $(GINKGO_DIR)/lib*/libginkgo*.dll)))
 ALL_GINKGO_LIBS_RELEASE = $(filter-out $(ALL_GINKGO_LIBS_DEBUG),$(ALL_GINKGO_LIBS))
 GINKGO_LINK = $(subst libginkgo,-lginkgo,$(ALL_GINKGO_LIBS_RELEASE))
 ifeq ($(GINKGO_BUILD_TYPE),Debug)
@@ -363,8 +349,7 @@ ifeq ($(GINKGO_BUILD_TYPE),Debug)
   endif
 else
 endif
-GINKGO_LIB = $(XLINKER)-rpath,$(GINKGO_LINK_LIB_DIR) -L$(GINKGO_LINK_LIB_DIR)\
- $(GINKGO_LINK)
+GINKGO_LIB = $(XLINKER)-rpath,$(GINKGO_LIB_DIR) -L$(GINKGO_LIB_DIR) $(GINKGO_LINK)
 
 # AmgX library configuration
 AMGX_DIR = @MFEM_DIR@/../amgx

examples/ex13p.cpp

@@ -5,7 +5,6 @@
 // Sample runs:  mpirun -np 4 ex13p -m ../data/star.mesh
 //               mpirun -np 4 ex13p -m ../data/square-disc.mesh -o 2 -n 4
 //               mpirun -np 4 ex13p -m ../data/beam-tet.mesh
-//               mpirun -np 4 ex13p -m ../data/beam-tet.mesh -nc -o 2 -rs 1
 //               mpirun -np 4 ex13p -m ../data/beam-hex.mesh
 //               mpirun -np 4 ex13p -m ../data/escher.mesh
 //               mpirun -np 4 ex13p -m ../data/fichera.mesh
@@ -55,7 +54,6 @@ int main(int argc, char *argv[])
    int par_ref_levels = 1;
    int order = 1;
    int nev = 5;
-   bool nc = false;
    bool visualization = 1;
    const char *device_config = "cpu";
 
@@ -71,9 +69,6 @@ int main(int argc, char *argv[])
                   " isoparametric space.");
    args.AddOption(&nev, "-n", "--num-eigs",
                   "Number of desired eigenmodes.");
-   args.AddOption(&nc, "-nc", "--non-conforming", "-c",
-                  "--conforming",
-                  "Mark the mesh as nonconforming before partitioning.");
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
@@ -103,10 +98,6 @@ int main(int argc, char *argv[])
    //    and volume meshes with the same code.
    Mesh *mesh = new Mesh(mesh_file, 1, 1);
    int dim = mesh->Dimension();
-   if (nc)
-   {
-      mesh->EnsureNCMesh(true);
-   }
 
    // 5. Refine the serial mesh on all processors to increase the resolution. In
    //    this example we do 'ref_levels' of uniform refinement (2 by default, or

examples/ex15p.cpp

@@ -13,7 +13,6 @@
 //               mpirun -np 4 ex15p -m ../data/square-disc-nurbs.mesh
 //               mpirun -np 4 ex15p -m ../data/disc-nurbs.mesh
 //               mpirun -np 4 ex15p -m ../data/fichera.mesh -tf 0.5
-//               mpirun -np 4 ex15p -m ../data/fichera-mixed.mesh -tf 0.5
 //               mpirun -np 4 ex15p -m ../data/ball-nurbs.mesh -tf 0.5
 //               mpirun -np 4 ex15p -m ../data/mobius-strip.mesh
 //               mpirun -np 4 ex15p -m ../data/amr-quad.mesh

examples/ex2p.cpp

@@ -100,6 +100,21 @@ int main(int argc, char *argv[])
    Device device(device_config);
    if (myid == 0) { device.Print(); }
 
+   if (mfem::Device::Allows(mfem::Backend::DEVICE_MASK))
+   {
+      HYPRE_SetMemoryLocation(HYPRE_MEMORY_DEVICE);
+      HYPRE_SetExecutionPolicy(HYPRE_EXEC_DEVICE);
+      HYPRE_DeviceInitialize();
+   }
+   else
+   {
+      HYPRE_SetMemoryLocation(HYPRE_MEMORY_HOST);
+      HYPRE_SetExecutionPolicy(HYPRE_EXEC_HOST);
+   }
+
+   auto loc = mfem::GetHypreMemoryLocation();
+   auto exec = mfem::GetHypreExecutionPolicy();
+
    // 4. Read the (serial) mesh from the given mesh file on all processors.  We
    //    can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
    //    and volume meshes with the same code.

examples/ex3.cpp

@@ -63,7 +63,6 @@ int main(int argc, char *argv[])
    int order = 1;
    bool static_cond = false;
    bool pa = false;
-   bool nc = false;
    const char *device_config = "cpu";
    bool visualization = 1;
 
@@ -78,9 +77,6 @@ int main(int argc, char *argv[])
                   "--no-static-condensation", "Enable static condensation.");
    args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
                   "--no-partial-assembly", "Enable Partial Assembly.");
-   args.AddOption(&nc, "-nc", "--non-conforming", "-c",
-                  "--conforming",
-                  "Mark the mesh as nonconforming before partitioning.");
    args.AddOption(&device_config, "-d", "--device",
                   "Device configuration string, see Device::Configure().");
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
@@ -106,11 +102,6 @@ int main(int argc, char *argv[])
    Mesh *mesh = new Mesh(mesh_file, 1, 1);
    dim = mesh->Dimension();
    int sdim = mesh->SpaceDimension();
-   if (nc)
-   {
-      // Can set to false to use conformal refinement for simplices.
-      mesh->EnsureNCMesh(true);
-   }
 
    // 4. Refine the mesh to increase the resolution. In this example we do
    //    'ref_levels' of uniform refinement. We choose 'ref_levels' to be the

examples/ex3p.cpp

@@ -5,7 +5,6 @@
 // Sample runs:  mpirun -np 4 ex3p -m ../data/star.mesh
 //               mpirun -np 4 ex3p -m ../data/square-disc.mesh -o 2
 //               mpirun -np 4 ex3p -m ../data/beam-tet.mesh
-//               mpirun -np 4 ex3p -m ../data/beam-tet.mesh -nc -o 2
 //               mpirun -np 4 ex3p -m ../data/beam-hex.mesh
 //               mpirun -np 4 ex3p -m ../data/beam-hex.mesh -o 2 -pa
 //               mpirun -np 4 ex3p -m ../data/escher.mesh
@@ -71,7 +70,6 @@ int main(int argc, char *argv[])
    int order = 1;
    bool static_cond = false;
    bool pa = false;
-   bool nc = false;
    const char *device_config = "cpu";
    bool visualization = true;
 #ifdef MFEM_USE_AMGX
@@ -89,9 +87,6 @@ int main(int argc, char *argv[])
                   "--no-static-condensation", "Enable static condensation.");
    args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
                   "--no-partial-assembly", "Enable Partial Assembly.");
-   args.AddOption(&nc, "-nc", "--non-conforming", "-c",
-                  "--conforming",
-                  "Mark the mesh as nonconforming before partitioning.");
    args.AddOption(&device_config, "-d", "--device",
                   "Device configuration string, see Device::Configure().");
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
@@ -129,11 +124,6 @@ int main(int argc, char *argv[])
    Mesh *mesh = new Mesh(mesh_file, 1, 1);
    dim = mesh->Dimension();
    int sdim = mesh->SpaceDimension();
-   if (nc)
-   {
-      // Can set to false to use conformal refinement for simplices.
-      mesh->EnsureNCMesh(true);
-   }
 
    // 5. Refine the serial mesh on all processors to increase the resolution. In
    //    this example we do 'ref_levels' of uniform refinement. We choose

examples/ex9.cpp

@@ -450,8 +450,7 @@ int main(int argc, char *argv[])
 
 // Implementation of class FE_Evolution
 FE_Evolution::FE_Evolution(BilinearForm &M_, BilinearForm &K_, const Vector &b_)
-   : TimeDependentOperator(M_.FESpace()->GetTrueVSize()),
-     M(M_), K(K_), b(b_), z(height)
+   : TimeDependentOperator(M_.Height()), M(M_), K(K_), b(b_), z(M_.Height())
 {
    Array<int> ess_tdof_list;
    if (M.GetAssemblyLevel() == AssemblyLevel::LEGACY)

examples/ex9p.cpp

@@ -659,9 +659,9 @@ int main(int argc, char *argv[])
 // Implementation of class FE_Evolution
 FE_Evolution::FE_Evolution(ParBilinearForm &M_, ParBilinearForm &K_,
                            const Vector &b_, PrecType prec_type)
-   : TimeDependentOperator(M_.ParFESpace()->GetTrueVSize()), b(b_),
+   : TimeDependentOperator(M_.Height()), b(b_),
      M_solver(M_.ParFESpace()->GetComm()),
-     z(height)
+     z(M_.Height())
 {
    if (M_.GetAssemblyLevel()==AssemblyLevel::LEGACY)
    {

examples/petsc/ex5p.cpp

@@ -17,7 +17,7 @@
 //               finite elements (velocity u) and piecewise discontinuous
 //               polynomials (pressure p).
 //
-//               The example demonstrates the use of the BlockOperator class, as
+//               The example demonstrates the use of the BlockMatrix class, as
 //               well as the collective saving of several grid functions in a
 //               VisIt (visit.llnl.gov) visualization format.
 //

examples/petsc/ex9p.cpp

@@ -520,10 +520,10 @@ int main(int argc, char *argv[])
 // Implementation of class FE_Evolution
 FE_Evolution::FE_Evolution(ParBilinearForm &M_, ParBilinearForm &K_,
                            const Vector &b_,bool M_in_lhs)
-   : TimeDependentOperator(M_.ParFESpace()->GetTrueVSize(), 0.0,
+   : TimeDependentOperator(M_.Height(), 0.0,
                            M_in_lhs ? TimeDependentOperator::IMPLICIT
                            : TimeDependentOperator::EXPLICIT),
-     b(b_), comm(M_.ParFESpace()->GetComm()), M_solver(comm), z(height),
+     b(b_), comm(M_.ParFESpace()->GetComm()), M_solver(comm), z(M_.Height()),
      iJacobian(NULL), rJacobian(NULL)
 {
    MAlev = M_.GetAssemblyLevel();

examples/sundials/ex9.cpp

@@ -476,8 +476,7 @@ int main(int argc, char *argv[])
 
 // Implementation of class FE_Evolution
 FE_Evolution::FE_Evolution(BilinearForm &M_, BilinearForm &K_, const Vector &b_)
-   : TimeDependentOperator(M_.FESpace()->GetTrueVSize()),
-     M(M_), K(K_), b(b_), z(height)
+   : TimeDependentOperator(M_.Height()), M(M_), K(K_), b(b_), z(M_.Height())
 {
    Array<int> ess_tdof_list;
    if (M.GetAssemblyLevel() == AssemblyLevel::LEGACY)

examples/sundials/ex9p.cpp

@@ -679,10 +679,10 @@ int main(int argc, char *argv[])
 // Implementation of class FE_Evolution
 FE_Evolution::FE_Evolution(ParBilinearForm &M_, ParBilinearForm &K_,
                            const Vector &b_, PrecType prec_type)
-   : TimeDependentOperator(M_.ParFESpace()->GetTrueVSize()),
+   : TimeDependentOperator(M_.Height()),
      b(b_),
      M_solver(M_.ParFESpace()->GetComm()),
-     z(height)
+     z(M_.Height())
 {
    if (M_.GetAssemblyLevel()==AssemblyLevel::LEGACY)
    {

fem/coefficient.cpp

@@ -220,12 +220,12 @@ double TransformedCoefficient::Eval(ElementTransformation &T,
 {
    if (Q2)
    {
-      return Transform2(Q1->Eval(T, ip, GetTime()),
-                        Q2->Eval(T, ip, GetTime()));
+      return (*Transform2)(Q1->Eval(T, ip, GetTime()),
+                           Q2->Eval(T, ip, GetTime()));
    }
    else
    {
-      return Transform1(Q1->Eval(T, ip, GetTime()));
+      return (*Transform1)(Q1->Eval(T, ip, GetTime()));
    }
 }
 

fem/coefficient.hpp

@@ -422,15 +422,15 @@ class TransformedCoefficient : public Coefficient
 private:
    Coefficient * Q1;
    Coefficient * Q2;
-   std::function<double(double)> Transform1;
-   std::function<double(double, double)> Transform2;
+   double (*Transform1)(double);
+   double (*Transform2)(double,double);
 
 public:
-   TransformedCoefficient (Coefficient * q, std::function<double(double)> F)
-      : Q1(q), Transform1(std::move(F)) { Q2 = 0; Transform2 = 0; }
+   TransformedCoefficient (Coefficient * q,double (*F)(double))
+      : Q1(q), Transform1(F) { Q2 = 0; Transform2 = 0; }
    TransformedCoefficient (Coefficient * q1,Coefficient * q2,
-                           std::function<double(double, double)> F)
-      : Q1(q1), Q2(q2), Transform2(std::move(F)) { Transform1 = 0; }
+                           double (*F)(double,double))
+      : Q1(q1), Q2(q2), Transform2(F) { Transform1 = 0; }
 
    /// Set the time for internally stored coefficients
    void SetTime(double t);

fem/complex_fem.cpp

@@ -1243,25 +1243,28 @@ ParSesquilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list,
          HypreParMatrix * Ah;
          A_i.Get(Ah);
          hypre_ParCSRMatrix *Aih = *Ah;
-#if !defined(HYPRE_USING_GPU)
-         ess_tdof_list.HostRead();
-         for (int k = 0; k < n; k++)
+         if (!HypreUsingGPU())
          {
-            const int j = ess_tdof_list[k];
-            Aih->diag->data[Aih->diag->i[j]] = 0.0;
+            ess_tdof_list.HostRead();
+            for (int k = 0; k < n; k++)
+            {
+               const int j = ess_tdof_list[k];
+               Aih->diag->data[Aih->diag->i[j]] = 0.0;
+            }
          }
-#else
-         Ah->HypreReadWrite();
-         const int *d_ess_tdof_list =
-            ess_tdof_list.GetMemory().Read(MemoryClass::DEVICE, n);
-         const int *d_diag_i = Aih->diag->i;
-         double *d_diag_data = Aih->diag->data;
-         MFEM_GPU_FORALL(k, n,
+         else
          {
-            const int j = d_ess_tdof_list[k];
-            d_diag_data[d_diag_i[j]] = 0.0;
-         });
-#endif
+            Ah->HypreReadWrite();
+            const int *d_ess_tdof_list =
+               ess_tdof_list.GetMemory().Read(MemoryClass::DEVICE, n);
+            const int *d_diag_i = Aih->diag->i;
+            double *d_diag_data = Aih->diag->data;
+            MFEM_GPU_FORALL(k, n,
+            {
+               const int j = d_ess_tdof_list[k];
+               d_diag_data[d_diag_i[j]] = 0.0;
+            });
+         }
       }
       else
       {

fem/eltrans.cpp

@@ -492,7 +492,6 @@ int IsoparametricTransformation::OrderGrad(const FiniteElement *fe) const
 void IsoparametricTransformation::Transform (const IntegrationPoint &ip,
                                              Vector &trans)
 {
-   MFEM_ASSERT(FElem != nullptr, "Must provide a valid FiniteElement object!");
    shape.SetSize(FElem->GetDof());
    trans.SetSize(PointMat.Height());
 

fem/fespace.cpp

@@ -2428,7 +2428,6 @@ void FiniteElementSpace::Construct()
       {
          // the simple case: all edges are of the same order
          nedofs = mesh->GetNEdges() * fec->GetNumDof(Geometry::SEGMENT, order);
-         var_edge_dofs.Clear(); // ensure any old var_edge_dof table is dumped.
       }
    }
 
@@ -2447,7 +2446,6 @@ void FiniteElementSpace::Construct()
          // the simple case: all faces are of the same geometry and order
          uni_fdof = fec->GetNumDof(mesh->GetFaceGeometry(0), order);
          nfdofs = mesh->GetNFaces() * uni_fdof;
-         var_face_dofs.Clear(); // ensure any old var_face_dof table is dumped.
       }
    }
 
@@ -2658,6 +2656,7 @@ int FiniteElementSpace::MakeDofTable(int ent_dim,
             int dofs = fec->GetNumDof(geom, order);
             list.Append(Connection(i, total_dofs));
             total_dofs += dofs;
+
             if (var_ent_order) { var_ent_order->Append(order); }
          }
       }
@@ -2668,6 +2667,7 @@ int FiniteElementSpace::MakeDofTable(int ent_dim,
 
    // build the table
    entity_dofs.MakeFromList(num_ent+1, list);
+
    return total_dofs;
 }
 
@@ -2996,14 +2996,7 @@ int FiniteElementSpace::GetFaceDofs(int face, Array<int> &dofs,
 
       order = !IsVariableOrder() ? fec->GetOrder() :
               var_face_orders[var_face_dofs.GetI()[face] + variant];
-      MFEM_ASSERT(fec->GetNumDof(fgeom, order) == nf, [&]()
-      {
-         std::stringstream msg;
-         msg << "fec->GetNumDof(" << (fgeom == Geometry::SQUARE ? "square" : "triangle")
-             << ", " << order << ") = " << fec->GetNumDof(fgeom, order) << " nf " << nf;
-         msg << " face " << face << " variant " << variant << std::endl;
-         return msg.str();
-      }());
+      MFEM_ASSERT(fec->GetNumDof(fgeom, order) == nf, "");
    }
    else
    {

fem/gridfunc.cpp

@@ -38,9 +38,8 @@ GridFunction::GridFunction(Mesh *m, std::istream &input)
    // Grid functions are stored on the device
    UseDevice(true);
 
-   owned_fes.reset(new FiniteElementSpace);
-   fes = owned_fes.get();
-   fec.reset(fes->Load(m, input));
+   fes = new FiniteElementSpace;
+   fec = fes->Load(m, input);
 
    skip_comment_lines(input, '#');
    istream::int_type next_char = input.peek();
@@ -82,11 +81,10 @@ GridFunction::GridFunction(Mesh *m, GridFunction *gf_array[], int num_pieces)
    int vdim, ordering;
 
    fes = gf_array[0]->FESpace();
-   fec.reset(FiniteElementCollection::New(fes->FEColl()->Name()));
+   fec = FiniteElementCollection::New(fes->FEColl()->Name());
    vdim = fes->GetVDim();
    ordering = fes->GetOrdering();
-   owned_fes.reset(new FiniteElementSpace(m, fec.get(), vdim, ordering));
-   fes = owned_fes.get();
+   fes = new FiniteElementSpace(m, fec, vdim, ordering);
    SetSize(fes->GetVSize());
 
    if (m->NURBSext)
@@ -155,9 +153,12 @@ GridFunction::GridFunction(Mesh *m, GridFunction *gf_array[], int num_pieces)
 
 void GridFunction::Destroy()
 {
-   owned_fes.reset();
-   fec.reset();
-   fes = nullptr;
+   if (fec)
+   {
+      delete fes;
+      delete fec;
+      fec = NULL;
+   }
 }
 
 void GridFunction::Update()
@@ -719,6 +720,56 @@ void GridFunction::GetVectorValues(int i, const IntegrationRule &ir,
    GetVectorValues(*Tr, ir, vals);
 }
 
+void be_to_bfe(Geometry::Type geom, int o, const IntegrationPoint &ip,
+               IntegrationPoint &fip)
+{
+   if (geom == Geometry::TRIANGLE)
+   {
+      if (o == 2)
+      {
+         fip.x = 1.0 - ip.x - ip.y;
+         fip.y = ip.x;
+      }
+      else if (o == 4)
+      {
+         fip.x = ip.y;
+         fip.y = 1.0 - ip.x - ip.y;
+      }
+      else
+      {
+         fip.x = ip.x;
+         fip.y = ip.y;
+      }
+      fip.z = ip.z;
+   }
+   else
+   {
+      if (o == 2)
+      {
+         fip.x = ip.y;
+         fip.y = 1.0 - ip.x;
+      }
+      else if (o == 4)
+      {
+         fip.x = 1.0 - ip.x;
+         fip.y = 1.0 - ip.y;
+      }
+      else if (o == 6)
+      {
+         fip.x = 1.0 - ip.y;
+         fip.y = ip.x;
+      }
+      else
+      {
+         fip.x = ip.x;
+         fip.y = ip.y;
+      }
+      fip.z = ip.z;
+   }
+   fip.weight = ip.weight;
+   fip.index  = ip.index;
+}
+
 double GridFunction::GetValue(ElementTransformation &T,
                               const IntegrationPoint &ip,
                               int comp, Vector *tr) const
@@ -783,15 +834,18 @@ double GridFunction::GetValue(ElementTransformation &T,
             // boundary so we'll evaluate it in the neighboring element.
             FaceElementTransformations * FET =
                fes->GetMesh()->GetBdrFaceTransformations(T.ElementNo);
-            MFEM_ASSERT(FET != nullptr,
-                        "FaceElementTransformation must be valid for a boundary element");
 
-            // Boundary elements and boundary faces may have different
+            // Boundary elements and Boundary Faces may have different
             // orientations so adjust the integration point if necessary.
-            int f, o;
-            fes->GetMesh()->GetBdrElementFace(T.ElementNo, &f, &o);
-            IntegrationPoint fip =
-               Mesh::TransformBdrElementToFace(FET->GetGeometryType(), o, ip);
+            int o = 0;
+            if (fes->GetMesh()->Dimension() == 3)
+            {
+               int f;
+               fes->GetMesh()->GetBdrElementFace(T.ElementNo, &f, &o);
+            }
+
+            IntegrationPoint fip;
+            be_to_bfe(FET->GetGeometryType(), o, ip, fip);
 
             // Compute and set the point in element 1 from fip
             FET->SetAllIntPoints(&fip);
@@ -919,15 +973,18 @@ void GridFunction::GetVectorValue(ElementTransformation &T,
             // the boundary so we'll evaluate it in the neighboring element.
             FaceElementTransformations * FET =
                fes->GetMesh()->GetBdrFaceTransformations(T.ElementNo);
-            MFEM_ASSERT(FET != nullptr,
-                        "FaceElementTransformation must be valid for a boundary element");
 
-            // Boundary elements and boundary faces may have different
+            // Boundary elements and Boundary Faces may have different
             // orientations so adjust the integration point if necessary.
-            int f, o;
-            fes->GetMesh()->GetBdrElementFace(T.ElementNo, &f, &o);
-            IntegrationPoint fip =
-               Mesh::TransformBdrElementToFace(FET->GetGeometryType(), o, ip);
+            int o = 0;
+            if (fes->GetMesh()->Dimension() == 3)
+            {
+               int f;
+               fes->GetMesh()->GetBdrElementFace(T.ElementNo, &f, &o);
+            }
+
+            IntegrationPoint fip;
+            be_to_bfe(FET->GetGeometryType(), o, ip, fip);
 
             // Compute and set the point in element 1 from fip
             FET->SetAllIntPoints(&fip);
@@ -940,8 +997,6 @@ void GridFunction::GetVectorValue(ElementTransformation &T,
       {
          FaceElementTransformations * FET =
             dynamic_cast<FaceElementTransformations *>(&T);
-         MFEM_ASSERT(FET != nullptr,
-                     "FaceElementTransformation must be valid for a boundary element");
 
          // Evaluate in neighboring element for both continuous and
          // discontinuous fields (the integration point in T1 should have
@@ -1060,10 +1115,11 @@ int GridFunction::GetFaceVectorValues(
    int i, int side, const IntegrationRule &ir,
    DenseMatrix &vals, DenseMatrix &tr) const
 {
-   int di;
+   int n, di;
    FaceElementTransformations *Transf;
 
-   IntegrationRule eir(ir.GetNPoints());  // ---
+   n = ir.GetNPoints();
+   IntegrationRule eir(n);  // ---
    Transf = fes->GetMesh()->GetFaceElementTransformations(i, 0);
    if (side == 2)
    {
@@ -1085,14 +1141,12 @@ int GridFunction::GetFaceVectorValues(
    if (di == 0)
    {
       Transf = fes->GetMesh()->GetFaceElementTransformations(i, 5);
-      MFEM_ASSERT(Transf != nullptr, "FaceElementTransformation cannot be null!");
       Transf->Loc1.Transform(ir, eir);
       GetVectorValues(*Transf->Elem1, eir, vals, &tr);
    }
    else
    {
       Transf = fes->GetMesh()->GetFaceElementTransformations(i, 10);
-      MFEM_ASSERT(Transf != nullptr, "FaceElementTransformation cannot be null!");
       Transf->Loc2.Transform(ir, eir);
       GetVectorValues(*Transf->Elem2, eir, vals, &tr);
    }
@@ -1450,13 +1504,17 @@ double GridFunction::GetDivergence(ElementTransformation &T) const
          FaceElementTransformations * FET =
             fes->GetMesh()->GetBdrFaceTransformations(T.ElementNo);
 
-         // Boundary elements and boundary faces may have different
+         // Boundary elements and Boundary Faces may have different
          // orientations so adjust the integration point if necessary.
-         int f, o;
-         fes->GetMesh()->GetBdrElementFace(T.ElementNo, &f, &o);
-         IntegrationPoint fip =
-            Mesh::TransformBdrElementToFace(FET->GetGeometryType(), o,
-                                            T.GetIntPoint());
+         int o = 0;
+         if (fes->GetMesh()->Dimension() == 3)
+         {
+            int f;
+            fes->GetMesh()->GetBdrElementFace(T.ElementNo, &f, &o);
+         }
+
+         IntegrationPoint fip;
+         be_to_bfe(FET->GetGeometryType(), o, T.GetIntPoint(), fip);
 
          // Compute and set the point in element 1 from fip
          FET->SetAllIntPoints(&fip);
@@ -1543,13 +1601,17 @@ void GridFunction::GetCurl(ElementTransformation &T, Vector &curl) const
          FaceElementTransformations * FET =
             fes->GetMesh()->GetBdrFaceTransformations(T.ElementNo);
 
-         // Boundary elements and boundary faces may have different
+         // Boundary elements and Boundary Faces may have different
          // orientations so adjust the integration point if necessary.
-         int f, o;
-         fes->GetMesh()->GetBdrElementFace(T.ElementNo, &f, &o);
-         IntegrationPoint fip =
-            Mesh::TransformBdrElementToFace(FET->GetGeometryType(), o,
-                                            T.GetIntPoint());
+         int o = 0;
+         if (fes->GetMesh()->Dimension() == 3)
+         {
+            int f;
+            fes->GetMesh()->GetBdrElementFace(T.ElementNo, &f, &o);
+         }
+
+         IntegrationPoint fip;
+         be_to_bfe(FET->GetGeometryType(), o, T.GetIntPoint(), fip);
 
          // Compute and set the point in element 1 from fip
          FET->SetAllIntPoints(&fip);
@@ -1608,13 +1670,17 @@ void GridFunction::GetGradient(ElementTransformation &T, Vector &grad) const
          FaceElementTransformations * FET =
             fes->GetMesh()->GetBdrFaceTransformations(T.ElementNo);
 
-         // Boundary elements and boundary faces may have different
+         // Boundary elements and Boundary Faces may have different
          // orientations so adjust the integration point if necessary.
-         int f, o;
-         fes->GetMesh()->GetBdrElementFace(T.ElementNo, &f, &o);
-         IntegrationPoint fip =
-            Mesh::TransformBdrElementToFace(FET->GetGeometryType(), o,
-                                            T.GetIntPoint());
+         int o = 0;
+         if (fes->GetMesh()->Dimension() == 3)
+         {
+            int f;
+            fes->GetMesh()->GetBdrElementFace(T.ElementNo, &f, &o);
+         }
+
+         IntegrationPoint fip;
+         be_to_bfe(FET->GetGeometryType(), o, T.GetIntPoint(), fip);
 
          // Compute and set the point in element 1 from fip
          FET->SetAllIntPoints(&fip);
@@ -1690,13 +1756,17 @@ void GridFunction::GetVectorGradient(
          FaceElementTransformations * FET =
             fes->GetMesh()->GetBdrFaceTransformations(T.ElementNo);
 
-         // Boundary elements and boundary faces may have different
+         // Boundary elements and Boundary Faces may have different
          // orientations so adjust the integration point if necessary.
-         int f, o;
-         fes->GetMesh()->GetBdrElementFace(T.ElementNo, &f, &o);
-         IntegrationPoint fip =
-            Mesh::TransformBdrElementToFace(FET->GetGeometryType(), o,
-                                            T.GetIntPoint());
+         int o = 0;
+         if (fes->GetMesh()->Dimension() == 3)
+         {
+            int f;
+            fes->GetMesh()->GetBdrElementFace(T.ElementNo, &f, &o);
+         }
+
+         IntegrationPoint fip;
+         be_to_bfe(FET->GetGeometryType(), o, T.GetIntPoint(), fip);
 
          // Compute and set the point in element 1 from fip
          FET->SetAllIntPoints(&fip);

fem/gridfunc.hpp

@@ -20,7 +20,6 @@
 #include "../general/adios2stream.hpp"
 #endif
 #include <limits>
-#include <memory>
 #include <ostream>
 #include <string>
 
@@ -31,13 +30,14 @@ namespace mfem
 class GridFunction : public Vector
 {
 protected:
-   /// FE space on which the grid function lives.
+   /// FE space on which the grid function lives. Owned if #fec is not NULL.
    FiniteElementSpace *fes;
 
    /** @brief Used when the grid function is read from a file. It can also be
-       set explicitly, see MakeOwner(). */
-   std::shared_ptr<FiniteElementCollection> fec;
-   std::shared_ptr<FiniteElementSpace> owned_fes;
+       set explicitly, see MakeOwner().
+
+       If not NULL, this pointer is owned by the GridFunction. */
+   FiniteElementCollection *fec;
 
    long fes_sequence; // see FiniteElementSpace::sequence, Mesh::sequence
 
@@ -74,6 +74,11 @@ public:
 
    GridFunction() { fes = NULL; fec = NULL; fes_sequence = 0; UseDevice(true); }
 
+   /// Copy constructor. The internal true-dof vector #t_vec is not copied.
+   GridFunction(const GridFunction &orig)
+      : Vector(orig), fes(orig.fes), fec(NULL), fes_sequence(orig.fes_sequence)
+   { UseDevice(true); }
+
    /// Construct a GridFunction associated with the FiniteElementSpace @a *f.
    GridFunction(FiniteElementSpace *f) : Vector(f->GetVSize())
    { fes = f; fec = NULL; fes_sequence = f->GetSequence(); UseDevice(true); }
@@ -102,12 +107,21 @@ public:
 
    GridFunction(Mesh *m, GridFunction *gf_array[], int num_pieces);
 
+   /// Copy assignment. Only the data of the base class Vector is copied.
+   /** It is assumed that this object and @a rhs use FiniteElementSpace%s that
+       have the same size.
+
+       @note Defining this method overwrites the implicitly defined copy
+       assignment operator. */
+   GridFunction &operator=(const GridFunction &rhs)
+   { return operator=((const Vector &)rhs); }
+
    /// Make the GridFunction the owner of #fec and #fes.
    /** If the new FiniteElementCollection, @a fec_, is NULL, ownership of #fec
        and #fes is taken away. */
-   void MakeOwner(FiniteElementCollection *fec_) { fec.reset(fec_); }
+   void MakeOwner(FiniteElementCollection *fec_) { fec = fec_; }
 
-   FiniteElementCollection *OwnFEC() { return fec.get(); }
+   FiniteElementCollection *OwnFEC() { return fec; }
 
    int VectorDim() const;
    int CurlDim() const;
@@ -740,6 +754,9 @@ public:
    /** @brief Write the GridFunction in STL format. Note that the mesh dimension
        must be 2 and that quad elements will be broken into two triangles.*/
    void SaveSTL(std::ostream &out, int TimesToRefine = 1);
+
+   /// Destroys grid function.
+   virtual ~GridFunction() { Destroy(); }
 };
 
 

fem/gslib.cpp

@@ -10,7 +10,6 @@
 // CONTRIBUTING.md for details.
 
 #include "gslib.hpp"
-#include "geom.hpp"
 
 #ifdef MFEM_USE_GSLIB
 
@@ -239,8 +238,7 @@ void FindPointsGSLIB::FindPoints(const Vector &point_pos,
    }
 
    // Map element number for simplices, and ref_pos from [-1,1] to [0,1] for
-   // both simplices and quads. Also sets code to 1 for points found on element
-   // faces/edges.
+   // both simplices and quads.
    MapRefPosAndElemIndices();
 }
 
@@ -683,9 +681,6 @@ void FindPointsGSLIB::MapRefPosAndElemIndices()
    int nptorig = points_cnt,
        npt = points_cnt;
 
-   // tolerance for point to be marked as on element edge/face
-   double btol = 1e-12;
-
    GridFunction *gf_rst_map_temp = NULL;
    int nptsend = 0;
 
@@ -699,7 +694,7 @@ void FindPointsGSLIB::MapRefPosAndElemIndices()
 
    // Pack data to send via crystal router
    struct gslib::array *outpt = new gslib::array;
-   struct out_pt { double r[3]; uint index, el, proc, code; };
+   struct out_pt { double r[3]; uint index, el, proc; };
    struct out_pt *pt;
    array_init(struct out_pt, outpt, nptsend);
    outpt->n=nptsend;
@@ -717,12 +712,12 @@ void FindPointsGSLIB::MapRefPosAndElemIndices()
       pt->index = index;
       pt->proc  = gsl_proc[index];
       pt->el    = gsl_elem[index];
-      pt->code  = gsl_code[index];
       ++pt;
    }
 
    // Transfer data to target MPI ranks
    sarray_transfer(struct out_pt, outpt, proc, 1, cr);
+
    // Map received points
    npt = outpt->n;
    pt = (struct out_pt *)outpt->ptr;
@@ -736,13 +731,7 @@ void FindPointsGSLIB::MapRefPosAndElemIndices()
       const Geometry::Type gt = fe->GetGeomType();
       pt->el = mesh_elem;
 
-      if (gt == Geometry::SQUARE || gt == Geometry::CUBE)
-      {
-         // check if it is on element boundary
-         pt->code = Geometry::CheckPoint(gt, ip, -btol) ? 0 : 1;
-         ++pt;
-         continue;
-      }
+      if (gt == Geometry::SQUARE || gt == Geometry::CUBE) { ++pt; continue; }
       else if (gt == Geometry::TRIANGLE)
       {
          gf_rst_map_temp = gf_rst_map[0];
@@ -769,10 +758,6 @@ void FindPointsGSLIB::MapRefPosAndElemIndices()
       {
          pt->r[d] = mfem_ref(d);
       }
-
-      // check if point is on element boundary
-      ip.Set3(&pt->r[0]);
-      pt->code = Geometry::CheckPoint(gt, ip, -btol) ? 0 : 1;
       ++pt;
    }
 
@@ -789,7 +774,6 @@ void FindPointsGSLIB::MapRefPosAndElemIndices()
       {
          gsl_mfem_ref(d + pt->index*dim) = pt->r[d];
       }
-      gsl_code[pt->index] = pt->code;
       ++pt;
    }
    array_free(outpt);
@@ -800,22 +784,12 @@ void FindPointsGSLIB::MapRefPosAndElemIndices()
    {
       if (gsl_code[index] != 2 && gsl_proc[index] == gsl_comm->id)
       {
-
-         IntegrationPoint ip;
-         Vector mfem_ref(gsl_mfem_ref.GetData()+index*dim, dim);
-         ip.Set2(mfem_ref.GetData());
-         if (dim == 3) { ip.z = mfem_ref(2); }
-
          const int elem = gsl_elem[index];
          const int mesh_elem = split_element_map[elem];
          const FiniteElement *fe = mesh->GetNodalFESpace()->GetFE(mesh_elem);
          const Geometry::Type gt = fe->GetGeomType();
          gsl_mfem_elem[index] = mesh_elem;
-         if (gt == Geometry::SQUARE || gt == Geometry::CUBE)
-         {
-            gsl_code[index] = Geometry::CheckPoint(gt, ip, -btol) ? 0 : 1;
-            continue;
-         }
+         if (gt == Geometry::SQUARE || gt == Geometry::CUBE) { continue; }
          else if (gt == Geometry::TRIANGLE)
          {
             gf_rst_map_temp = gf_rst_map[0];
@@ -834,12 +808,11 @@ void FindPointsGSLIB::MapRefPosAndElemIndices()
          }
 
          int local_elem = split_element_index[elem];
-         gf_rst_map_temp->GetVectorValue(local_elem, ip, mfem_ref);
-
-         // Check if the point is on element boundary
+         IntegrationPoint ip;
+         Vector mfem_ref(gsl_mfem_ref.GetData()+index*dim, dim);
          ip.Set2(mfem_ref.GetData());
          if (dim == 3) { ip.z = mfem_ref(2); }
-         gsl_code[index]  = Geometry::CheckPoint(gt, ip, -btol) ? 0 : 1;
+         gf_rst_map_temp->GetVectorValue(local_elem, ip, mfem_ref);
       }
    }
 }

fem/lininteg.cpp

@@ -395,8 +395,8 @@ void VectorBoundaryLFIntegrator::AssembleRHSElementVect(
    {
       const IntegrationPoint &ip = ir->IntPoint(i);
 
-      Tr.SetIntPoint (&ip);
       Q.Eval(vec, Tr, ip);
+      Tr.SetIntPoint (&ip);
       vec *= Tr.Weight() * ip.weight;
       el.CalcShape(ip, shape);
       for (int k = 0; k < vdim; k++)

fem/lor/lor_ams.cpp

@@ -291,7 +291,8 @@ void BatchedLOR_AMS::FormCoordinateVectors(const Vector &X_vert)
    const auto ltdof_ldof = HypreRead(R->GetMemoryJ());
 
    // Go from E-vector format directly to T-vector format
-   MFEM_HYPRE_FORALL(i, ntdofs,
+   //MFEM_HYPRE_FORALL(i, ntdofs,
+   mfem::forall_switch(HypreUsingGPU(), ntdofs, [=] MFEM_HOST_DEVICE (int i)
    {
       const int j = d_offsets[ltdof_ldof[i]];
       for (int c = 0; c < sdim; ++c)

fem/pfespace.cpp

@@ -1946,8 +1946,9 @@ struct PMatrixRow
    void AddRow(const PMatrixRow &other, double coef)
    {
       elems.reserve(elems.size() + other.elems.size());
-      for (const PMatrixElement &oei : other.elems)
+      for (unsigned i = 0; i < other.elems.size(); i++)
       {
+         const PMatrixElement &oei = other.elems[i];
          elems.push_back(
             PMatrixElement(oei.column, oei.stride, coef * oei.value));
       }
@@ -2057,7 +2058,7 @@ void NeighborRowMessage::Encode(int rank)
    for (unsigned i = 0; i < rows.size(); i++)
    {
       const RowInfo &ri = rows[i];
-      const MeshId &id = *pncmesh->GetNCList(ri.entity).GetMeshIdAndType(ri.index).id;
+      const MeshId &id = pncmesh->GetNCList(ri.entity).LookUp(ri.index);
       ent_ids[ri.entity].Append(id);
       row_idx[ri.entity].Append(i);
       group_ids[ri.entity].Append(ri.group);
@@ -2131,33 +2132,25 @@ void NeighborRowMessage::Decode(int rank)
    rows.clear();
    rows.reserve(nrows);
 
-   // read rows ent = {0,1,2} means vertex, edge and face entity
+   // read rows
    for (int ent = 0, gi = 0; ent < 3; ent++)
    {
-      // extract the vertex list, edge list or face list.
       const Array<MeshId> &ids = ent_ids[ent];
       for (int i = 0; i < ids.Size(); i++)
       {
          const MeshId &id = ids[i];
-         // read the particular element dof value off the stream.
          int edof = bin_io::read<int>(stream);
 
-         // Handle orientation and sign change. This flips the sign on dofs
-         // where necessary, and for edges and faces also reorders if flipped,
-         // i.e. an edge with 1 -> 2 -> 3 -> 4 might become -4 -> -3 -> -2 -> -1
-         // This cannot treat all face dofs, as they can have rotations and
-         // reflections.
-         const int *ind = nullptr;
-         Geometry::Type geom = Geometry::Type::INVALID;
+         // handle orientation and sign change
+         const int *ind = NULL;
          if (ent == 1)
          {
-            // edge NC orientation is element defined.
             int eo = pncmesh->GetEdgeNCOrientation(id);
             ind = fec->DofOrderForOrientation(Geometry::SEGMENT, eo);
          }
          else if (ent == 2)
          {
-            geom = pncmesh->GetFaceGeometry(id.index);
+            Geometry::Type geom = pncmesh->GetFaceGeometry(id.index);
             int fo = pncmesh->GetFaceOrientation(id.index);
             ind = fec->DofOrderForOrientation(geom, fo);
          }
@@ -2172,14 +2165,13 @@ void NeighborRowMessage::Decode(int rank)
          // If edof arrived with a negative index, flip it, and the scaling.
          double s = (edof < 0) ? -1.0 : 1.0;
          edof = (edof < 0) ? -1 - edof : edof;
+
          if (ind && (edof = ind[edof]) < 0)
          {
             edof = -1 - edof;
             s *= -1.0;
          }
 
-         // Create a row for this entity, recording the index of the mesh
-         // element
          rows.push_back(RowInfo(ent, id.index, edof, group_ids[gi++]));
          rows.back().row.read(stream, s);
 
@@ -2189,69 +2181,6 @@ void NeighborRowMessage::Decode(int rank)
                    << rows.back().index << ", edof " << rows.back().edof
                    << std::endl;
 #endif
-
-         if (ent == 2 && fec->GetContType() == FiniteElementCollection::TANGENTIAL
-             && !Geometry::IsTensorProduct(geom))
-         {
-            // ND face dofs need to be processed together, as the transformation
-            // is given by a 2x2 matrix, so we manually apply an extra increment
-            // to the loop counter and add in a new row. Once these rows are
-            // placed, they represent the Identity transformation. To map across
-            // the processor boundary, we also need to apply a Primal
-            // Transformation (see doftrans.hpp) to a notional "global dof"
-            // orientation. For simplicity we perform the action of these 2x2
-            // matrices manually using the AddRow capability, followed by a
-            // Collapse.
-
-            // To perform the operations, we add and subtract initial versions
-            // of the rows, that represent [1 0; 0 1] in row major notation. The
-            // first row represents the 1 at (0,0) in [1 0; 0 1] The second row
-            // represents the 1 at (1,1) in [1 0; 0 1]
-
-            // We can safely bind this reference as rows was reserved above so
-            // there is no hidden copying that could result in a dangling
-            // reference.
-            auto &first_row = rows.back().row;
-            // This is the first "fundamental unit" used in the transformation.
-            const auto initial_first_row = first_row;
-            // Extract the next dof too, and apply any dof order transformation
-            // expected.
-            const MeshId &next_id = ids[++i];
-            const int fo = pncmesh->GetFaceOrientation(next_id.index);
-            ind = fec->DofOrderForOrientation(geom, fo);
-            edof = bin_io::read<int>(stream);
-
-            // If edof arrived with a negative index, flip it, and the scaling.
-            s = (edof < 0) ? -1.0 : 1.0;
-            edof = (edof < 0) ? -1 - edof : edof;
-            if (ind && (edof = ind[edof]) < 0)
-            {
-               edof = -1 - edof;
-               s *= -1.0;
-            }
-            rows.push_back(RowInfo(ent, next_id.index, edof, group_ids[gi++]));
-            rows.back().row.read(stream, s);
-            auto &second_row = rows.back().row;
-
-            // This is the second "fundamental unit" used in the transformation.
-            const auto initial_second_row = second_row;
-
-            const auto T = [&fo]()
-            {
-               auto T = ND_StatelessDofTransformation::GetFaceTransform(fo);
-               T(0,0) -= 1;
-               T(1,1) -= 1;
-               return T;
-            }();
-
-            first_row.AddRow(initial_first_row, T(0,0));
-            first_row.AddRow(initial_second_row, T(0,1));
-            second_row.AddRow(initial_first_row, T(1,0));
-            second_row.AddRow(initial_second_row, T(1,1));
-
-            first_row.Collapse();
-            second_row.Collapse();
-         }
       }
    }
 }
@@ -2379,6 +2308,12 @@ int ParFiniteElementSpace
                                        Array<int> *dof_tdof,
                                        bool partial) const
 {
+   // TODO: general face DOF transformations in NeighborRowMessage::Decode()
+   MFEM_VERIFY(!(fec->GetOrder() >= 2
+                 && pmesh->HasGeometry(Geometry::TETRAHEDRON)
+                 && fec->GetContType() == FiniteElementCollection::TANGENTIAL),
+               "Nedelec NC tets of order >= 2 are not supported yet.");
+
    const bool dg = (nvdofs == 0 && nedofs == 0 && nfdofs == 0);
 
 #ifdef MFEM_PMATRIX_STATS

fem/pfespace.hpp

@@ -190,15 +190,15 @@ private:
    /// Updates the internal mesh pointer. @warning @a new_mesh must be
    /// <b>topologically identical</b> to the existing mesh. Used if the address
    /// of the Mesh object has changed, e.g. in @a Mesh::Swap.
-   void UpdateMeshPointer(Mesh *new_mesh) override;
+   virtual void UpdateMeshPointer(Mesh *new_mesh);
 
    /// Copies the prolongation and restriction matrices from @a fes.
    ///
    /// Used for low order preconditioning on non-conforming meshes. If the DOFs
    /// require a permutation, it will be supplied by non-NULL @a perm. NULL @a
    /// perm indicates that no permutation is required.
-   void CopyProlongationAndRestriction(const FiniteElementSpace &fes,
-                                       const Array<int> *perm) override;
+   virtual void CopyProlongationAndRestriction(const FiniteElementSpace &fes,
+                                               const Array<int> *perm);
 
 public:
    // Face-neighbor data
@@ -282,32 +282,32 @@ public:
    { return Dof_TrueDof_Matrix()->GetGlobalNumCols(); }
 
    /// Return the number of local vector true dofs.
-   int GetTrueVSize() const override { return ltdof_size; }
+   virtual int GetTrueVSize() const { return ltdof_size; }
 
    /// Returns indexes of degrees of freedom in array dofs for i'th element.
-   DofTransformation *GetElementDofs(int i, Array<int> &dofs) const override;
+   virtual DofTransformation *GetElementDofs(int i, Array<int> &dofs) const;
 
    /// Returns indexes of degrees of freedom for i'th boundary element.
-   DofTransformation *GetBdrElementDofs(int i, Array<int> &dofs) const override;
+   virtual DofTransformation *GetBdrElementDofs(int i, Array<int> &dofs) const;
 
    /** Returns the indexes of the degrees of freedom for i'th face
        including the dofs for the edges and the vertices of the face. */
-   int GetFaceDofs(int i, Array<int> &dofs, int variant = 0) const override;
+   virtual int GetFaceDofs(int i, Array<int> &dofs, int variant = 0) const;
 
    /** Returns pointer to the FiniteElement in the FiniteElementCollection
        associated with i'th element in the mesh object. If @a i is greater than
        or equal to the number of local mesh elements, @a i will be interpreted
        as a shifted index of a face neighbor element. */
-   const FiniteElement *GetFE(int i) const override;
+   virtual const FiniteElement *GetFE(int i) const;
 
    /** Returns an Operator that converts L-vectors to E-vectors on each face.
        The parallel version is different from the serial one because of the
        presence of shared faces. Shared faces are treated as interior faces,
        the returned operator handles the communication needed to get the
        shared face values from other MPI ranks */
-   const FaceRestriction *GetFaceRestriction(
+   virtual const FaceRestriction *GetFaceRestriction(
       ElementDofOrdering f_ordering, FaceType type,
-      L2FaceValues mul = L2FaceValues::DoubleValued) const override;
+      L2FaceValues mul = L2FaceValues::DoubleValued) const;
 
    void GetSharedEdgeDofs(int group, int ei, Array<int> &dofs) const;
    void GetSharedTriangleDofs(int group, int fi, Array<int> &dofs) const;
@@ -347,15 +347,15 @@ public:
    void Synchronize(Array<int> &ldof_marker) const;
 
    /// Determine the boundary degrees of freedom
-   void GetEssentialVDofs(const Array<int> &bdr_attr_is_ess,
-                          Array<int> &ess_dofs,
-                          int component = -1) const override;
+   virtual void GetEssentialVDofs(const Array<int> &bdr_attr_is_ess,
+                                  Array<int> &ess_dofs,
+                                  int component = -1) const;
 
    /** Get a list of essential true dofs, ess_tdof_list, corresponding to the
        boundary attributes marked in the array bdr_attr_is_ess. */
-   void GetEssentialTrueDofs(const Array<int> &bdr_attr_is_ess,
-                             Array<int> &ess_tdof_list,
-                             int component = -1) override;
+   virtual void GetEssentialTrueDofs(const Array<int> &bdr_attr_is_ess,
+                                     Array<int> &ess_tdof_list,
+                                     int component = -1);
 
    /** If the given ldof is owned by the current processor, return its local
        tdof number, otherwise return -1 */
@@ -370,13 +370,13 @@ public:
    HYPRE_BigInt GetMyDofOffset() const;
    HYPRE_BigInt GetMyTDofOffset() const;
 
-   const Operator *GetProlongationMatrix() const override;
+   virtual const Operator *GetProlongationMatrix() const;
    /** Get an Operator that performs the action of GetRestrictionMatrix(),
        but potentially with a non-assembled optimized matrix-free
        implementation. */
-   const Operator *GetRestrictionOperator() const override;
+   virtual const Operator *GetRestrictionOperator() const;
    /// Get the R matrix which restricts a local dof vector to true dof vector.
-   const SparseMatrix *GetRestrictionMatrix() const override
+   virtual const SparseMatrix *GetRestrictionMatrix() const
    { Dof_TrueDof_Matrix(); return R; }
 
    // Face-neighbor functions
@@ -402,15 +402,15 @@ public:
    // Transfer parallel true-dof data from coarse_fes, defined on a coarse mesh,
    // to this FE space, defined on a refined mesh. See full documentation in the
    // base class, FiniteElementSpace::GetTrueTransferOperator.
-   void GetTrueTransferOperator(const FiniteElementSpace &coarse_fes,
-                                OperatorHandle &T) const override;
+   virtual void GetTrueTransferOperator(const FiniteElementSpace &coarse_fes,
+                                        OperatorHandle &T) const;
 
    /** Reflect changes in the mesh. Calculate one of the refinement/derefinement
        /rebalance matrices, unless want_transform is false. */
-   void Update(bool want_transform = true) override;
+   virtual void Update(bool want_transform = true);
 
    /// Free ParGridFunction transformation matrix (if any), to save memory.
-   void UpdatesFinished() override
+   virtual void UpdatesFinished()
    {
       FiniteElementSpace::UpdatesFinished();
       old_dof_offsets.DeleteAll();
@@ -442,9 +442,9 @@ public:
 
    const GroupCommunicator &GetGroupCommunicator() const;
 
-   void Mult(const Vector &x, Vector &y) const override;
+   virtual void Mult(const Vector &x, Vector &y) const;
 
-   void MultTranspose(const Vector &x, Vector &y) const override;
+   virtual void MultTranspose(const Vector &x, Vector &y) const;
 };
 
 /// Auxiliary device class used by ParFiniteElementSpace.
@@ -493,9 +493,9 @@ public:
 
    virtual ~DeviceConformingProlongationOperator();
 
-   void Mult(const Vector &x, Vector &y) const override;
+   virtual void Mult(const Vector &x, Vector &y) const;
 
-   void MultTranspose(const Vector &x, Vector &y) const override;
+   virtual void MultTranspose(const Vector &x, Vector &y) const;
 };
 
 }

fem/pgridfunc.cpp

@@ -39,11 +39,9 @@ ParGridFunction::ParGridFunction(ParMesh *pmesh, const GridFunction *gf,
 {
    const FiniteElementSpace *glob_fes = gf->FESpace();
    // duplicate the FiniteElementCollection from 'gf'
-   fec.reset(FiniteElementCollection::New(glob_fes->FEColl()->Name()));
+   fec = FiniteElementCollection::New(glob_fes->FEColl()->Name());
    // create a local ParFiniteElementSpace from the global one:
-   fes = pfes = new ParFiniteElementSpace(pmesh, glob_fes, partitioning,
-                                          fec.get());
-   owned_fes.reset(pfes);
+   fes = pfes = new ParFiniteElementSpace(pmesh, glob_fes, partitioning, fec);
    SetSize(pfes->GetVSize());
 
    if (partitioning)
@@ -83,10 +81,10 @@ ParGridFunction::ParGridFunction(ParMesh *pmesh, std::istream &input)
    : GridFunction(pmesh, input)
 {
    // Convert the FiniteElementSpace, fes, to a ParFiniteElementSpace:
-   pfes = new ParFiniteElementSpace(pmesh, fec.get(), fes->GetVDim(),
+   pfes = new ParFiniteElementSpace(pmesh, fec, fes->GetVDim(),
                                     fes->GetOrdering());
+   delete fes;
    fes = pfes;
-   owned_fes.reset(pfes);
 }
 
 void ParGridFunction::Update()
@@ -471,16 +469,17 @@ void ParGridFunction::GetVectorValue(ElementTransformation &T,
    }
 
    Array<int> vdofs;
-   DofTransformation * doftrans = pfes->GetFaceNbrElementVDofs(nbr_el_no, vdofs);
+   DofTransformation * doftrans = pfes->GetFaceNbrElementVDofs(nbr_el_no,
+                                                               vdofs);
+   const FiniteElement *fe = pfes->GetFaceNbrFE(nbr_el_no);
+
+   int dof = fe->GetDof();
    Vector loc_data;
    face_nbr_data.GetSubVector(vdofs, loc_data);
    if (doftrans)
    {
       doftrans->InvTransformPrimal(loc_data);
    }
-
-   const FiniteElement *fe = pfes->GetFaceNbrFE(nbr_el_no);
-   const int dof = fe->GetDof();
    if (fe->GetRangeType() == FiniteElement::SCALAR)
    {
       Vector shape(dof);

fem/pgridfunc.hpp

@@ -49,6 +49,10 @@ protected:
 public:
    ParGridFunction() { pfes = NULL; }
 
+   /// Copy constructor. The internal vector #face_nbr_data is not copied.
+   ParGridFunction(const ParGridFunction &orig)
+      : GridFunction(orig), pfes(orig.pfes) { }
+
    ParGridFunction(ParFiniteElementSpace *pf) : GridFunction(pf), pfes(pf) { }
 
    /// Construct a ParGridFunction using previously allocated array @a data.
@@ -89,6 +93,15 @@ public:
        constructed. The new ParGridFunction assumes ownership of both. */
    ParGridFunction(ParMesh *pmesh, std::istream &input);
 
+   /// Copy assignment. Only the data of the base class Vector is copied.
+   /** It is assumed that this object and @a rhs use ParFiniteElementSpace%s
+       that have the same size.
+
+       @note Defining this method overwrites the implicitly defined copy
+       assignment operator. */
+   ParGridFunction &operator=(const ParGridFunction &rhs)
+   { return operator=((const Vector &)rhs); }
+
    /// Assign constant values to the ParGridFunction data.
    ParGridFunction &operator=(double value)
    { GridFunction::operator=(value); return *this; }
@@ -199,7 +212,7 @@ public:
    virtual double GetValue(int i, const IntegrationPoint &ip,
                            int vdim = 1) const;
    double GetValue(ElementTransformation &T)
-   { return GetValue(T, T.GetIntPoint()); }
+   { return GetValue(T.ElementNo, T.GetIntPoint()); }
 
    // Redefine to handle the case when T describes a face-neighbor element
    virtual double GetValue(ElementTransformation &T, const IntegrationPoint &ip,
@@ -446,6 +459,8 @@ public:
 
    /// Merge the local grid functions
    void SaveAsOne(std::ostream &out = mfem::out) const;
+
+   virtual ~ParGridFunction() { }
 };
 
 

general/mem_manager.hpp

@@ -17,9 +17,6 @@
 #include <cstring> // std::memcpy
 #include <type_traits> // std::is_const
 #include <cstddef> // std::max_align_t
-#ifdef MFEM_USE_MPI
-#include <HYPRE_config.h> // HYPRE_USING_GPU
-#endif
 
 namespace mfem
 {
@@ -994,17 +991,7 @@ inline void Memory<T>::MakeAlias(const Memory &base, int offset, int size)
    h_ptr = base.h_ptr + offset;
    if (!(base.flags & Registered))
    {
-      if (
-#if !defined(HYPRE_USING_GPU)
-         // If the following condition is true then MemoryManager::Exists()
-         // should also be true:
-         IsDeviceMemory(MemoryManager::GetDeviceMemoryType())
-#else
-         // When HYPRE_USING_GPU is defined we always register the 'base' if
-         // the MemoryManager::Exists():
-         MemoryManager::Exists()
-#endif
-      )
+      if (IsDeviceMemory(MemoryManager::GetDeviceMemoryType()))
       {
          // Register 'base':
          MemoryManager::Register_(base.h_ptr, nullptr, base.capacity*sizeof(T),

general/table.cpp

@@ -218,8 +218,7 @@ void Table::SetIJ(int *newI, int *newJ, int newsize)
 
 int Table::Push(int i, int j)
 {
-   MFEM_ASSERT(i >=0 &&
-               i<size, "Index out of bounds.  i = " << i << " size " << size);
+   MFEM_ASSERT( i >=0 && i<size, "Index out of bounds.  i = "<<i);
 
    for (int k = I[i], end = I[i+1]; k < end; k++)
    {

linalg/ginkgo.cpp

@@ -21,7 +21,6 @@
 #include <iomanip>
 #include <algorithm>
 #include <cmath>
-#include <cstring>
 
 namespace mfem
 {
@@ -29,14 +28,8 @@ namespace mfem
 namespace Ginkgo
 {
 
-// Create a GinkgoExecutor of type exec_type.
 GinkgoExecutor::GinkgoExecutor(ExecType exec_type)
 {
-#if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
-   gko::version_info gko_version = gko::version_info::get();
-   bool gko_with_omp_support = (strcmp(gko_version.omp_version.tag,
-                                       "not compiled") != 0);
-#endif
    switch (exec_type)
    {
       case GinkgoExecutor::REFERENCE:
@@ -56,23 +49,13 @@ GinkgoExecutor::GinkgoExecutor(ExecType exec_type)
 #ifdef MFEM_USE_CUDA
             int current_device = 0;
             MFEM_GPU_CHECK(cudaGetDevice(&current_device));
-            if (gko_with_omp_support)
-            {
-               executor = gko::CudaExecutor::create(current_device,
-                                                    gko::OmpExecutor::create());
-            }
-            else
-            {
-               executor = gko::CudaExecutor::create(current_device,
-                                                    gko::ReferenceExecutor::create());
-            }
+            executor = gko::CudaExecutor::create(current_device,
+                                                 gko::OmpExecutor::create());
 #endif
          }
          else
-         {
             MFEM_ABORT("gko::CudaExecutor::get_num_devices() did not report "
                        "any valid devices.");
-         }
          break;
       }
       case GinkgoExecutor::HIP:
@@ -82,111 +65,24 @@ GinkgoExecutor::GinkgoExecutor(ExecType exec_type)
 #ifdef MFEM_USE_HIP
             int current_device = 0;
             MFEM_GPU_CHECK(hipGetDevice(&current_device));
-            if (gko_with_omp_support)
-            {
-               executor = gko::HipExecutor::create(current_device,
-                                                   gko::OmpExecutor::create());
-            }
-            else
-            {
-               executor = gko::HipExecutor::create(current_device,
-                                                   gko::ReferenceExecutor::create());
-            }
+            executor = gko::HipExecutor::create(current_device,
+                                                gko::OmpExecutor::create());
 #endif
          }
          else
-         {
-            MFEM_ABORT("gko::HipExecutor::get_num_devices() did not report "
-                       "any valid devices.");
-         }
+            mfem::err << "gko::HipExecutor::get_num_devices() did not report "
+                      << "any valid devices" << std::endl;
          break;
       }
       default:
-         MFEM_ABORT("Invalid ExecType specified");
+         mfem::err << "Invalid ExecType specified" << std::endl;
    }
 }
 
-// Create a GinkgoExecutor of type exec_type, with host_exec_type for the
-// related CPU Executor (only applicable to GPU backends).
-GinkgoExecutor::GinkgoExecutor(ExecType exec_type, ExecType host_exec_type)
-{
-   switch (exec_type)
-   {
-      case GinkgoExecutor::REFERENCE:
-      {
-         MFEM_WARNING("Parameter host_exec_type ignored for CPU GinkgoExecutor.");
-         executor = gko::ReferenceExecutor::create();
-         break;
-      }
-      case GinkgoExecutor::OMP:
-      {
-         MFEM_WARNING("Parameter host_exec_type ignored for CPU GinkgoExecutor.");
-         executor = gko::OmpExecutor::create();
-         break;
-      }
-      case GinkgoExecutor::CUDA:
-      {
-         if (gko::CudaExecutor::get_num_devices() > 0)
-         {
-#ifdef MFEM_USE_CUDA
-            int current_device = 0;
-            MFEM_GPU_CHECK(cudaGetDevice(&current_device));
-            if (host_exec_type == GinkgoExecutor::OMP)
-            {
-               executor = gko::CudaExecutor::create(current_device,
-                                                    gko::OmpExecutor::create());
-            }
-            else
-            {
-               executor = gko::CudaExecutor::create(current_device,
-                                                    gko::ReferenceExecutor::create());
-            }
-#endif
-         }
-         else
-         {
-            MFEM_ABORT("gko::CudaExecutor::get_num_devices() did not report "
-                       "any valid devices.");
-         }
-         break;
-      }
-      case GinkgoExecutor::HIP:
-      {
-         if (gko::HipExecutor::get_num_devices() > 0)
-         {
-#ifdef MFEM_USE_HIP
-            int current_device = 0;
-            MFEM_GPU_CHECK(hipGetDevice(&current_device));
-            if (host_exec_type == GinkgoExecutor::OMP)
-            {
-               executor = gko::HipExecutor::create(current_device,
-                                                   gko::OmpExecutor::create());
-            }
-            else
-            {
-               executor = gko::HipExecutor::create(current_device,
-                                                   gko::ReferenceExecutor::create());
-            }
-#endif
-         }
-         else
-         {
-            MFEM_ABORT("gko::HipExecutor::get_num_devices() did not report "
-                       "any valid devices.");
-         }
-         break;
-      }
-      default:
-         MFEM_ABORT("Invalid ExecType specified");
-   }
-}
-
-// Create a GinkgoExecutor to match MFEM's device configuration.
 GinkgoExecutor::GinkgoExecutor(Device &mfem_device)
 {
-   gko::version_info gko_version = gko::version_info::get();
-   bool gko_with_omp_support = (strcmp(gko_version.omp_version.tag,
-                                       "not compiled") != 0);
+
+   // Pick "best match" Executor based on MFEM device configuration.
    if (mfem_device.Allows(Backend::CUDA_MASK))
    {
       if (gko::CudaExecutor::get_num_devices() > 0)
@@ -194,23 +90,13 @@ GinkgoExecutor::GinkgoExecutor(Device &mfem_device)
 #ifdef MFEM_USE_CUDA
          int current_device = 0;
          MFEM_GPU_CHECK(cudaGetDevice(&current_device));
-         if (gko_with_omp_support)
-         {
-            executor = gko::CudaExecutor::create(current_device,
-                                                 gko::OmpExecutor::create());
-         }
-         else
-         {
-            executor = gko::CudaExecutor::create(current_device,
-                                                 gko::ReferenceExecutor::create());
-         }
+         executor = gko::CudaExecutor::create(current_device,
+                                              gko::OmpExecutor::create());
 #endif
       }
       else
-      {
          MFEM_ABORT("gko::CudaExecutor::get_num_devices() did not report "
                     "any valid devices.");
-      }
    }
    else if (mfem_device.Allows(Backend::HIP_MASK))
    {
@@ -219,123 +105,16 @@ GinkgoExecutor::GinkgoExecutor(Device &mfem_device)
 #ifdef MFEM_USE_HIP
          int current_device = 0;
          MFEM_GPU_CHECK(hipGetDevice(&current_device));
-         if (gko_with_omp_support)
-         {
-            executor = gko::HipExecutor::create(current_device,
-                                                gko::OmpExecutor::create());
-         }
-         else
-         {
-            executor = gko::HipExecutor::create(current_device,
-                                                gko::ReferenceExecutor::create());
-         }
+         executor = gko::HipExecutor::create(current_device, gko::OmpExecutor::create());
 #endif
       }
       else
-      {
          MFEM_ABORT("gko::HipExecutor::get_num_devices() did not report "
                     "any valid devices.");
-      }
    }
    else
    {
-      if (mfem_device.Allows(Backend::OMP_MASK))
-      {
-         // Also use OpenMP for Ginkgo, if Ginkgo supports it
-         if (gko_with_omp_support)
-         {
-            executor = gko::OmpExecutor::create();
-         }
-         else
-         {
-            executor = gko::ReferenceExecutor::create();
-         }
-      }
-      else
-      {
-         executor = gko::ReferenceExecutor::create();
-      }
-   }
-}
-
-// Create a GinkgoExecutor to match MFEM's device configuration, with
-// a specific host_exec_type for the associated CPU Executor (only
-// applicable to GPU backends).
-GinkgoExecutor::GinkgoExecutor(Device &mfem_device, ExecType host_exec_type)
-{
-
-   if (mfem_device.Allows(Backend::CUDA_MASK))
-   {
-      if (gko::CudaExecutor::get_num_devices() > 0)
-      {
-#ifdef MFEM_USE_CUDA
-         int current_device = 0;
-         MFEM_GPU_CHECK(cudaGetDevice(&current_device));
-         if (host_exec_type == GinkgoExecutor::OMP)
-         {
-            executor = gko::CudaExecutor::create(current_device,
-                                                 gko::OmpExecutor::create());
-         }
-         else
-         {
-            executor = gko::CudaExecutor::create(current_device,
-                                                 gko::ReferenceExecutor::create());
-         }
-#endif
-      }
-      else
-      {
-         MFEM_ABORT("gko::CudaExecutor::get_num_devices() did not report "
-                    "any valid devices.");
-      }
-   }
-   else if (mfem_device.Allows(Backend::HIP_MASK))
-   {
-      if (gko::HipExecutor::get_num_devices() > 0)
-      {
-#ifdef MFEM_USE_HIP
-         int current_device = 0;
-         MFEM_GPU_CHECK(hipGetDevice(&current_device));
-         if (host_exec_type == GinkgoExecutor::OMP)
-         {
-            executor = gko::HipExecutor::create(current_device,
-                                                gko::OmpExecutor::create());
-         }
-         else
-         {
-            executor = gko::HipExecutor::create(current_device,
-                                                gko::ReferenceExecutor::create());
-         }
-#endif
-      }
-      else
-      {
-         MFEM_ABORT("gko::HipExecutor::get_num_devices() did not report "
-                    "any valid devices.");
-      }
-   }
-   else
-   {
-      MFEM_WARNING("Parameter host_exec_type ignored for CPU GinkgoExecutor.");
-      if (mfem_device.Allows(Backend::OMP_MASK))
-      {
-         // Also use OpenMP for Ginkgo, if Ginkgo supports it
-         gko::version_info gko_version = gko::version_info::get();
-         bool gko_with_omp_support = (strcmp(gko_version.omp_version.tag,
-                                             "not compiled") != 0);
-         if (gko_with_omp_support)
-         {
-            executor = gko::OmpExecutor::create();
-         }
-         else
-         {
-            executor = gko::ReferenceExecutor::create();
-         }
-      }
-      else
-      {
-         executor = gko::ReferenceExecutor::create();
-      }
+      executor = gko::OmpExecutor::create();
    }
 }
 
@@ -416,7 +195,7 @@ const
                            gko::log::Logger::criterion_check_completed_mask);
 #endif
    residual_logger = std::make_shared<ResidualLogger<>>(executor,
-                                                        system_oper.get(),b);
+                                                        gko::lend(system_oper),b);
 
 }
 
@@ -555,7 +334,7 @@ GinkgoIterativeSolver::Mult(const Vector &x, Vector &y) const
 
    // Create the logger object to log some data from the solvers to confirm
    // convergence.
-   initialize_ginkgo_log(gko_x.get());
+   initialize_ginkgo_log(gko::lend(gko_x));
 
    MFEM_VERIFY(convergence_logger, "convergence logger not initialized" );
    if (print_level==1)
@@ -571,11 +350,7 @@ GinkgoIterativeSolver::Mult(const Vector &x, Vector &y) const
    combined_factory->add_logger(convergence_logger);
 
    // Finally, apply the solver to x and get the solution in y.
-#if MFEM_GINKGO_VERSION < 10600
    solver->apply(gko::lend(gko_x), gko::lend(gko_y));
-#else
-   solver->apply(gko_x, gko_y);
-#endif
 
    // Get the number of iterations taken to converge to the solution.
    final_iter = convergence_logger->get_num_iterations();
@@ -692,10 +467,6 @@ void GinkgoIterativeSolver::SetOperator(const Operator &op)
                        new OperatorWrapper(executor, op.Height(), &op));
    }
 
-   // Set MFEM Solver size values
-   height = op.Height();
-   width = op.Width();
-
    // Generate the solver from the solver using the system matrix or operator.
    solver = solver_gen->generate(system_oper);
 }
@@ -1107,11 +878,7 @@ GinkgoPreconditioner::Mult(const Vector &x, Vector &y) const
                             gko_array<double>::view(executor,
                                                     y.Size(),
                                                     y.ReadWrite(on_device)), 1);
-#if MFEM_GINKGO_VERSION < 10600
    generated_precond.get()->apply(gko::lend(gko_x), gko::lend(gko_y));
-#else
-   generated_precond.get()->apply(gko_x, gko_y);
-#endif
 }
 
 void GinkgoPreconditioner::SetOperator(const Operator &op)
@@ -1150,10 +917,6 @@ void GinkgoPreconditioner::SetOperator(const Operator &op)
 
    generated_precond = precond_gen->generate(gko::give(gko_matrix));
    has_generated_precond = true;
-
-   // Set MFEM Solver size values
-   height = op.Height();
-   width = op.Width();
 }
 
 

linalg/ginkgo.hpp

@@ -265,13 +265,9 @@ double compute_norm(const gko::matrix::Dense<ValueType> *b)
    // Initialize a result scalar containing the value 0.0.
    auto b_norm = gko::initialize<gko::matrix::Dense<ValueType>>({0.0}, exec);
    // Use the dense `compute_norm2` function to compute the norm.
-#if MFEM_GINKGO_VERSION < 10600
-   b->compute_norm2(gko::lend(b_norm));
-#else
-   b->compute_norm2(b_norm);
-#endif
+   b->compute_norm2(lend(b_norm));
    // Use the other utility function to return the norm contained in `b_norm``
-   return std::pow(get_norm(b_norm.get()),2);
+   return std::pow(get_norm(lend(b_norm)),2);
 }
 
 /**
@@ -323,43 +319,84 @@ struct ResidualLogger : gko::log::Logger
 
    using gko_dense = gko::matrix::Dense<ValueType>;
 
-   // Ginkgo 1.5 and older: version for solver that doesn't log implicit res norm
-   void on_iteration_complete(const gko::LinOp *op,
-                              const gko::size_type &iteration,
-                              const gko::LinOp *residual,
-                              const gko::LinOp *solution,
-                              const gko::LinOp *residual_norm) const override
-   {
-      iteration_complete_core(iteration, residual, solution, residual_norm,
-                              nullptr);
-   }
-   // Ginkgo 1.5 and older: version with implicit residual norm
-   void on_iteration_complete(const gko::LinOp *op,
+   // Customize the logging hook which is called every time an iteration is
+   // completed
+   void on_iteration_complete(const gko::LinOp *,
                               const gko::size_type &iteration,
                               const gko::LinOp *residual,
                               const gko::LinOp *solution,
                               const gko::LinOp *residual_norm,
                               const gko::LinOp *implicit_sq_residual_norm) const override
    {
-      iteration_complete_core(iteration, residual, solution, residual_norm,
-                              implicit_sq_residual_norm);
+      // If the solver shares the current solution vector and we want to
+      // compute the residual from that
+      if (solution && compute_real_residual)
+      {
+         // Store the matrix's executor
+         auto exec = matrix->get_executor();
+         // Compute the real residual vector by calling apply on the system
+         // First, compute res = A * x
+         matrix->apply(gko::lend(solution), gko::lend(res));
+         // Now do res = res - b, depending on which vector/oper type
+         // Check if b is a Ginkgo vector or wrapped MFEM Vector
+         if (dynamic_cast<const VectorWrapper*>(b))
+         {
+            const VectorWrapper *b_cast = gko::as<const VectorWrapper>(b);
+            // Copy the MFEM Vector stored in b
+            VectorWrapper *res_cast = gko::as<VectorWrapper>(res);
+            res_cast->get_mfem_vec_ref() -= b_cast->get_mfem_vec_const_ref();
+         }
+         else
+         {
+            // Create a scalar containing the value -1.0
+            auto neg_one = gko::initialize<gko_dense>({-1.0}, exec);
+            res->add_scaled(gko::lend(neg_one), gko::lend(b));
+         }
+
+         // Compute the norm of the residual vector and add it to the
+         // `residual_norms` vector
+         residual_norms.push_back(compute_norm(gko::lend(res)));
+      }
+      else
+      {
+         // If the solver shares an implicit or recurrent residual norm, log its value
+         if (implicit_sq_residual_norm)
+         {
+            auto dense_norm = gko::as<gko_dense>(implicit_sq_residual_norm);
+            // Add the norm to the `residual_norms` vector
+            residual_norms.push_back(get_norm(dense_norm));
+            // Otherwise, use the recurrent residual vector
+         }
+         else if (residual_norm)
+         {
+            auto dense_norm = gko::as<gko_dense>(residual_norm);
+            // Add the norm to the `residual_norms` vector
+            residual_norms.push_back(get_norm(dense_norm));
+            // Otherwise, use the recurrent residual vector
+         }
+         else
+         {
+            auto dense_residual = gko::as<gko_dense>(residual);
+            // Compute the residual vector's norm
+            auto norm = compute_norm(gko::lend(dense_residual));
+            // Add the computed norm to the `residual_norms` vector
+            residual_norms.push_back(norm);
+         }
+      }
+      // Add the current iteration number to the `iterations` vector
+      iterations.push_back(iteration);
    }
-#if MFEM_GINKGO_VERSION > 10500
-   // Ginkgo 1.6 and newer
+
+   // Version for solver that doesn't log implicit res norm
    void on_iteration_complete(const gko::LinOp *op,
-                              const gko::LinOp *rhs,
-                              const gko::LinOp *solution,
                               const gko::size_type &iteration,
                               const gko::LinOp *residual,
-                              const gko::LinOp *residual_norm,
-                              const gko::LinOp *implicit_sq_residual_norm,
-                              const gko::array<gko::stopping_status>* status,
-                              bool stopped) const override
+                              const gko::LinOp *solution,
+                              const gko::LinOp *residual_norm) const override
    {
-      iteration_complete_core(iteration, residual, solution, residual_norm,
-                              implicit_sq_residual_norm);
+      on_iteration_complete(op, iteration, residual, solution, residual_norm,
+                            nullptr);
    }
-#endif
 
    // Construct the logger and store the system matrix and b vectors
    ResidualLogger(std::shared_ptr<const gko::Executor> exec,
@@ -391,89 +428,6 @@ struct ResidualLogger : gko::log::Logger
    }
 
 private:
-   // Customize the logging hook which is called every time an iteration is
-   // completed.
-   void iteration_complete_core(const gko::size_type &iteration,
-                                const gko::LinOp *residual,
-                                const gko::LinOp *solution,
-                                const gko::LinOp *residual_norm,
-                                const gko::LinOp *implicit_sq_residual_norm) const
-   {
-      // If the solver shares the current solution vector and we want to
-      // compute the residual from that
-      if (solution && compute_real_residual)
-      {
-         // Store the matrix's executor
-         auto exec = matrix->get_executor();
-         // Compute the real residual vector by calling apply on the system
-         // First, compute res = A * x
-#if MFEM_GINKGO_VERSION < 10600
-         matrix->apply(gko::lend(solution), gko::lend(res));
-#else
-         matrix->apply(solution, res);
-#endif
-         // Now do res = res - b, depending on which vector/oper type
-         // Check if b is a Ginkgo vector or wrapped MFEM Vector
-         if (dynamic_cast<const VectorWrapper*>(b))
-         {
-            const VectorWrapper *b_cast = gko::as<const VectorWrapper>(b);
-            // Copy the MFEM Vector stored in b
-            VectorWrapper *res_cast = gko::as<VectorWrapper>(res);
-            res_cast->get_mfem_vec_ref() -= b_cast->get_mfem_vec_const_ref();
-         }
-         else
-         {
-            // Create a scalar containing the value -1.0
-            auto neg_one = gko::initialize<gko_dense>({-1.0}, exec);
-#if MFEM_GINKGO_VERSION < 10600
-            res->add_scaled(gko::lend(neg_one), gko::lend(b));
-#else
-            res->add_scaled(neg_one, b);
-#endif
-         }
-
-         // Compute the norm of the residual vector and add it to the
-         // `residual_norms` vector
-#if MFEM_GINKGO_VERSION < 10600
-         residual_norms.push_back(compute_norm(gko::lend(res)));
-#else
-         residual_norms.push_back(compute_norm(res));
-#endif
-      }
-      else
-      {
-         // If the solver shares an implicit or recurrent residual norm, log its value
-         if (implicit_sq_residual_norm)
-         {
-            auto dense_norm = gko::as<gko_dense>(implicit_sq_residual_norm);
-            // Add the norm to the `residual_norms` vector
-            residual_norms.push_back(get_norm(dense_norm));
-            // Otherwise, use the recurrent residual vector
-         }
-         else if (residual_norm)
-         {
-            auto dense_norm = gko::as<gko_dense>(residual_norm);
-            // Add the norm to the `residual_norms` vector
-            residual_norms.push_back(get_norm(dense_norm));
-            // Otherwise, use the recurrent residual vector
-         }
-         else
-         {
-            auto dense_residual = gko::as<gko_dense>(residual);
-            // Compute the residual vector's norm
-#if MFEM_GINKGO_VERSION < 10600
-            auto norm = compute_norm(gko::lend(dense_residual));
-#else
-            auto norm = compute_norm(dense_residual);
-#endif
-            // Add the computed norm to the `residual_norms` vector
-            residual_norms.push_back(norm);
-         }
-      }
-      // Add the current iteration number to the `iterations` vector
-      iterations.push_back(iteration);
-   }
-
    // Pointer to the system matrix
    const gko::LinOp *matrix;
    // Pointer to the right hand sides
@@ -515,40 +469,18 @@ public:
    /**
     * Constructor.
     * Takes an @p GinkgoExecType argument and creates an Executor.
-    * In Ginkgo, GPU Executors must have an associated host Executor.
-    * This routine will select a CPU Executor based on the OpenMP support
-    * for Ginkgo.
     */
    GinkgoExecutor(ExecType exec_type);
 
-   /**
-    * Constructor.
-    * Takes an @p GinkgoExecType argument and creates an Executor.
-    * In Ginkgo, GPU Executors must have an associated host Executor.
-    * This routine allows for explicite setting of the CPU Executor
-    * for GPU backends.
-    */
-   GinkgoExecutor(ExecType exec_type, ExecType host_exec_type);
-
    /**
     * Constructor.
     * Takes an MFEM @p Device object and creates an Executor
     * that "matches" (e.g., if MFEM is using the CPU, Ginkgo
-    * will choose the Reference or OmpExecutor based on MFEM's
-    * configuration and Ginkgo's capabilities; if MFEM is using
-    * CUDA, Ginkgo will choose the CudaExecutor with a default
-    * CPU Executor based on Ginkgo's OpenMP support).
+    * will choose the OmpExecutor; if MFEM is using CUDA,
+    * Ginkgo will choose the CudaExecutor).
     */
    GinkgoExecutor(Device &mfem_device);
 
-   /**
-    * Constructor.
-    * Takes an MFEM @p Device object and creates an Executor
-    * that "matches", but allows the user to specify the host
-    * Executor for GPU backends.
-    */
-   GinkgoExecutor(Device &mfem_device, ExecType host_exec_type);
-
    /**
     * Destructor.
     */

linalg/hypre.cpp

@@ -121,22 +121,27 @@ bool CanShallowCopy(const Memory<T> &src, MemoryClass mc)
 inline void HypreParVector::_SetDataAndSize_()
 {
    hypre_Vector *x_loc = hypre_ParVectorLocalVector(x);
-#if !defined(HYPRE_USING_GPU)
-   SetDataAndSize(hypre_VectorData(x_loc),
-                  internal::to_int(hypre_VectorSize(x_loc)));
-#else
-   size = internal::to_int(hypre_VectorSize(x_loc));
-   MemoryType mt = (hypre_VectorMemoryLocation(x_loc) == HYPRE_MEMORY_HOST
-                    ? MemoryType::HOST : GetHypreMemoryType());
-   if (hypre_VectorData(x_loc) != NULL)
+#if defined(HYPRE_USING_GPU)
+   if (HypreUsingGPU())
    {
-      data.Wrap(hypre_VectorData(x_loc), size, mt, false);
+      size = internal::to_int(hypre_VectorSize(x_loc));
+      MemoryType mt = (hypre_VectorMemoryLocation(x_loc) == HYPRE_MEMORY_HOST
+                       ? MemoryType::HOST : GetHypreMemoryType());
+      if (hypre_VectorData(x_loc) != NULL)
+      {
+         data.Wrap(hypre_VectorData(x_loc), size, mt, false);
+      }
+      else
+      {
+         data.Reset();
+      }
    }
    else
-   {
-      data.Reset();
-   }
 #endif
+   {
+      SetDataAndSize(hypre_VectorData(x_loc),
+                     internal::to_int(hypre_VectorSize(x_loc)));
+   }
 }
 
 HypreParVector::HypreParVector(MPI_Comm comm, HYPRE_BigInt glob_size,
@@ -264,6 +269,7 @@ void HypreParVector::WrapHypreParVector(hypre_ParVector *y, bool owner)
 Vector * HypreParVector::GlobalVector() const
 {
    hypre_Vector *hv = hypre_ParVectorToVectorAll(*this);
+   if (hv == nullptr) { return nullptr; }
    Vector *v = new Vector(hv->data, internal::to_int(hv->size));
    v->MakeDataOwner();
    hypre_SeqVectorSetDataOwner(hv,0);
@@ -315,7 +321,7 @@ void HypreParVector::HypreRead() const
    hypre_VectorData(x_loc) =
       const_cast<double*>(data.Read(GetHypreMemoryClass(), size));
 #ifdef HYPRE_USING_GPU
-   hypre_VectorMemoryLocation(x_loc) = HYPRE_MEMORY_DEVICE;
+   hypre_VectorMemoryLocation(x_loc) = mfem::GetHypreMemoryLocation();
 #endif
 }
 
@@ -324,7 +330,7 @@ void HypreParVector::HypreReadWrite()
    hypre_Vector *x_loc = hypre_ParVectorLocalVector(x);
    hypre_VectorData(x_loc) = data.ReadWrite(GetHypreMemoryClass(), size);
 #ifdef HYPRE_USING_GPU
-   hypre_VectorMemoryLocation(x_loc) = HYPRE_MEMORY_DEVICE;
+   hypre_VectorMemoryLocation(x_loc) = mfem::GetHypreMemoryLocation();
 #endif
 }
 
@@ -333,7 +339,7 @@ void HypreParVector::HypreWrite()
    hypre_Vector *x_loc = hypre_ParVectorLocalVector(x);
    hypre_VectorData(x_loc) = data.Write(GetHypreMemoryClass(), size);
 #ifdef HYPRE_USING_GPU
-   hypre_VectorMemoryLocation(x_loc) = HYPRE_MEMORY_DEVICE;
+   hypre_VectorMemoryLocation(x_loc) = mfem::GetHypreMemoryLocation();
 #endif
 }
 
@@ -347,7 +353,7 @@ void HypreParVector::WrapMemoryRead(const Memory<double> &mem)
    hypre_VectorData(x_loc) =
       const_cast<double*>(mem.Read(GetHypreMemoryClass(), size));
 #ifdef HYPRE_USING_GPU
-   hypre_VectorMemoryLocation(x_loc) = HYPRE_MEMORY_DEVICE;
+   hypre_VectorMemoryLocation(x_loc) = mfem::GetHypreMemoryLocation();
 #endif
    data.MakeAlias(mem, 0, size);
 }
@@ -361,7 +367,7 @@ void HypreParVector::WrapMemoryReadWrite(Memory<double> &mem)
    hypre_Vector *x_loc = hypre_ParVectorLocalVector(x);
    hypre_VectorData(x_loc) = mem.ReadWrite(GetHypreMemoryClass(), size);
 #ifdef HYPRE_USING_GPU
-   hypre_VectorMemoryLocation(x_loc) = HYPRE_MEMORY_DEVICE;
+   hypre_VectorMemoryLocation(x_loc) = mfem::GetHypreMemoryLocation();
 #endif
    data.MakeAlias(mem, 0, size);
 }
@@ -375,7 +381,7 @@ void HypreParVector::WrapMemoryWrite(Memory<double> &mem)
    hypre_Vector *x_loc = hypre_ParVectorLocalVector(x);
    hypre_VectorData(x_loc) = mem.Write(GetHypreMemoryClass(), size);
 #ifdef HYPRE_USING_GPU
-   hypre_VectorMemoryLocation(x_loc) = HYPRE_MEMORY_DEVICE;
+   hypre_VectorMemoryLocation(x_loc) = mfem::GetHypreMemoryLocation();
 #endif
    data.MakeAlias(mem, 0, size);
 }
@@ -535,6 +541,10 @@ void HypreParMatrix::Init()
 
 void HypreParMatrix::Read(MemoryClass mc) const
 {
+   if (GetHypreMemoryLocation() == HYPRE_MEMORY_HOST)
+   {
+      mc = Device::GetHostMemoryClass();
+   }
    hypre_CSRMatrix *diag = hypre_ParCSRMatrixDiag(A);
    hypre_CSRMatrix *offd = hypre_ParCSRMatrixOffd(A);
    const int num_rows = NumRows();
@@ -547,8 +557,8 @@ void HypreParMatrix::Read(MemoryClass mc) const
    offd->j = const_cast<HYPRE_Int*>(mem_offd.J.Read(mc, offd_nnz));
    offd->data = const_cast<double*>(mem_offd.data.Read(mc, offd_nnz));
 #if MFEM_HYPRE_VERSION >= 21800
-   decltype(diag->memory_location) ml =
-      (mc != GetHypreMemoryClass() ? HYPRE_MEMORY_HOST : HYPRE_MEMORY_DEVICE);
+   decltype(diag->memory_location) ml = (mc == MemoryClass::HOST) ?
+                                        HYPRE_MEMORY_HOST : GetHypreMemoryLocation();
    diag->memory_location = ml;
    offd->memory_location = ml;
 #endif
@@ -556,6 +566,10 @@ void HypreParMatrix::Read(MemoryClass mc) const
 
 void HypreParMatrix::ReadWrite(MemoryClass mc)
 {
+   if (GetHypreMemoryLocation() == HYPRE_MEMORY_HOST)
+   {
+      mc = Device::GetHostMemoryClass();
+   }
    hypre_CSRMatrix *diag = hypre_ParCSRMatrixDiag(A);
    hypre_CSRMatrix *offd = hypre_ParCSRMatrixOffd(A);
    const int num_rows = NumRows();
@@ -568,8 +582,8 @@ void HypreParMatrix::ReadWrite(MemoryClass mc)
    offd->j = mem_offd.J.ReadWrite(mc, offd_nnz);
    offd->data = mem_offd.data.ReadWrite(mc, offd_nnz);
 #if MFEM_HYPRE_VERSION >= 21800
-   decltype(diag->memory_location) ml =
-      (mc != GetHypreMemoryClass() ? HYPRE_MEMORY_HOST : HYPRE_MEMORY_DEVICE);
+   decltype(diag->memory_location) ml = (mc == MemoryClass::HOST) ?
+                                        HYPRE_MEMORY_HOST : GetHypreMemoryLocation();
    diag->memory_location = ml;
    offd->memory_location = ml;
 #endif
@@ -577,6 +591,10 @@ void HypreParMatrix::ReadWrite(MemoryClass mc)
 
 void HypreParMatrix::Write(MemoryClass mc, bool set_diag, bool set_offd)
 {
+   if (GetHypreMemoryLocation() == HYPRE_MEMORY_HOST)
+   {
+      mc = Device::GetHostMemoryClass();
+   }
    hypre_CSRMatrix *diag = hypre_ParCSRMatrixDiag(A);
    hypre_CSRMatrix *offd = hypre_ParCSRMatrixOffd(A);
    if (set_diag)
@@ -592,8 +610,8 @@ void HypreParMatrix::Write(MemoryClass mc, bool set_diag, bool set_offd)
       offd->data = mem_offd.data.Write(mc, mem_offd.data.Capacity());
    }
 #if MFEM_HYPRE_VERSION >= 21800
-   decltype(diag->memory_location) ml =
-      (mc != GetHypreMemoryClass() ? HYPRE_MEMORY_HOST : HYPRE_MEMORY_DEVICE);
+   decltype(diag->memory_location) ml = (mc == MemoryClass::HOST) ?
+                                        HYPRE_MEMORY_HOST : GetHypreMemoryLocation();
    if (set_diag) { diag->memory_location = ml; }
    if (set_offd) { offd->memory_location = ml; }
 #endif
@@ -759,7 +777,7 @@ signed char HypreParMatrix::HypreCsrToMem(hypre_CSRMatrix *h_mat,
       h_mat->data = mem.data.ReadWrite(hypre_mc, nnz);
       h_mat->owns_data = 0;
 #if MFEM_HYPRE_VERSION >= 21800
-      h_mat->memory_location = HYPRE_MEMORY_DEVICE;
+      h_mat->memory_location = mfem::GetHypreMemoryLocation();
 #endif
       return 3;
    }
@@ -1360,8 +1378,11 @@ hypre_ParCSRMatrix* HypreParMatrix::StealData()
    MFEM_ASSERT(ParCSROwner, "");
    hypre_ParCSRMatrix *R = A;
 #ifdef HYPRE_USING_GPU
-   if (diagOwner == -1) { HostReadWrite(); }
-   else { HypreReadWrite(); }
+   if (HypreUsingGPU())
+   {
+      if (diagOwner == -1) { HostReadWrite(); }
+      else { HypreReadWrite(); }
+   }
 #endif
    ParCSROwner = false;
    Destroy();
@@ -1709,7 +1730,10 @@ void HypreParMatrix::EnsureMultTranspose() const
 #if (MFEM_HYPRE_VERSION == 22500 && HYPRE_DEVELOP_NUMBER >= 1) || \
     (MFEM_HYPRE_VERSION > 22500)
 #ifdef HYPRE_USING_GPU
-   hypre_ParCSRMatrixLocalTranspose(A);
+   if (HypreUsingGPU())
+   {
+      hypre_ParCSRMatrixLocalTranspose(A);
+   }
 #endif
 #endif
 }
@@ -1719,15 +1743,18 @@ void HypreParMatrix::ResetTranspose() const
 #if (MFEM_HYPRE_VERSION == 22500 && HYPRE_DEVELOP_NUMBER >= 1) || \
     (MFEM_HYPRE_VERSION > 22500)
 #ifdef HYPRE_USING_GPU
-   if (A->diagT)
+   if (HypreUsingGPU())
    {
-      hypre_CSRMatrixDestroy(A->diagT);
-      A->diagT = NULL;
-   }
-   if (A->offdT)
-   {
-      hypre_CSRMatrixDestroy(A->offdT);
-      A->offdT = NULL;
+      if (A->diagT)
+      {
+         hypre_CSRMatrixDestroy(A->diagT);
+         A->diagT = NULL;
+      }
+      if (A->offdT)
+      {
+         hypre_CSRMatrixDestroy(A->offdT);
+         A->offdT = NULL;
+      }
    }
 #endif
 #endif
@@ -2416,8 +2443,10 @@ void HypreParMatrix::EliminateBC(const Array<int> &ess_dofs,
       }
 
       // Which of the local rows are to be eliminated?
-      MFEM_HYPRE_FORALL(i, diag_nrows, eliminate_row[i] = 0; );
-      MFEM_HYPRE_FORALL(i, n_ess_dofs, eliminate_row[ess_dofs_d[i]] = 1; );
+      //MFEM_HYPRE_FORALL(i, diag_nrows, eliminate_row[i] = 0; ); // TODO TMS
+      mfem::forall_switch(HypreUsingGPU(), diag_nrows, [=] MFEM_HOST_DEVICE (int i) { eliminate_row[i] = 0; });
+      //MFEM_HYPRE_FORALL(i, n_ess_dofs, eliminate_row[ess_dofs_d[i]] = 1; ); // TODO TMS
+      mfem::forall_switch(HypreUsingGPU(), n_ess_dofs, [=] MFEM_HOST_DEVICE (int i) { eliminate_row[ess_dofs_d[i]] = 1; });
 
       // Use a matvec communication pattern to find (in eliminate_col) which of
       // the local offd columns are to be eliminated
@@ -2428,26 +2457,37 @@ void HypreParMatrix::EliminateBC(const Array<int> &ess_dofs,
 
       HYPRE_Int *send_map_elmts;
 #if defined(HYPRE_USING_GPU)
-      hypre_ParCSRCommPkgCopySendMapElmtsToDevice(comm_pkg);
-      send_map_elmts = hypre_ParCSRCommPkgDeviceSendMapElmts(comm_pkg);
-#else
-      send_map_elmts = hypre_ParCSRCommPkgSendMapElmts(comm_pkg);
+      if (HypreUsingGPU())
+      {
+         hypre_ParCSRCommPkgCopySendMapElmtsToDevice(comm_pkg);
+         send_map_elmts = hypre_ParCSRCommPkgDeviceSendMapElmts(comm_pkg);
+      }
+      else
 #endif
-      MFEM_HYPRE_FORALL(i, int_buf_sz,
+      {
+         send_map_elmts = hypre_ParCSRCommPkgSendMapElmts(comm_pkg);
+      }
+      //MFEM_HYPRE_FORALL(i, int_buf_sz,
+      mfem::forall_switch(HypreUsingGPU(), int_buf_sz, [=] MFEM_HOST_DEVICE (int i)
       {
          int k = send_map_elmts[i];
          int_buf_data[i] = eliminate_row[k];
       });
 
 #if defined(HYPRE_USING_GPU)
-      // Try to use device-aware MPI for the communication if available
-      comm_handle = hypre_ParCSRCommHandleCreate_v2(
-                       11, comm_pkg, HYPRE_MEMORY_DEVICE, int_buf_data,
-                       HYPRE_MEMORY_DEVICE, eliminate_col);
-#else
-      comm_handle = hypre_ParCSRCommHandleCreate(
-                       11, comm_pkg, int_buf_data, eliminate_col );
+      if (HypreUsingGPU())
+      {
+         // Try to use device-aware MPI for the communication if available
+         comm_handle = hypre_ParCSRCommHandleCreate_v2(
+                          11, comm_pkg, HYPRE_MEMORY_DEVICE, int_buf_data,
+                          HYPRE_MEMORY_DEVICE, eliminate_col);
+      }
+      else
 #endif
+      {
+         comm_handle = hypre_ParCSRCommHandleCreate(
+                          11, comm_pkg, int_buf_data, eliminate_col );
+      }
    }
 
    // Eliminate rows and columns in the diagonal block
@@ -2456,7 +2496,8 @@ void HypreParMatrix::EliminateBC(const Array<int> &ess_dofs,
       const auto J = diag->j;
       auto data = diag->data;
 
-      MFEM_HYPRE_FORALL(i, n_ess_dofs,
+      //MFEM_HYPRE_FORALL(i, n_ess_dofs,
+      mfem::forall_switch(HypreUsingGPU(), n_ess_dofs, [=] MFEM_HOST_DEVICE (int i)
       {
          const int idof = ess_dofs_d[i];
          for (int j=I[idof]; j<I[idof+1]; ++j)
@@ -2494,7 +2535,8 @@ void HypreParMatrix::EliminateBC(const Array<int> &ess_dofs,
    {
       const auto I = offd->i;
       auto data = offd->data;
-      MFEM_HYPRE_FORALL(i, n_ess_dofs,
+      //MFEM_HYPRE_FORALL(i, n_ess_dofs,
+      mfem::forall_switch(HypreUsingGPU(), n_ess_dofs, [=] MFEM_HOST_DEVICE (int i)
       {
          const int idof = ess_dofs_d[i];
          for (int j=I[idof]; j<I[idof+1]; ++j)
@@ -2515,7 +2557,8 @@ void HypreParMatrix::EliminateBC(const Array<int> &ess_dofs,
       const auto I = offd->i;
       const auto J = offd->j;
       auto data = offd->data;
-      MFEM_HYPRE_FORALL(i, nrows_offd,
+      //MFEM_HYPRE_FORALL(i, nrows_offd,
+      mfem::forall_switch(HypreUsingGPU(), nrows_offd, [=] MFEM_HOST_DEVICE (int i)
       {
          for (int j=I[i]; j<I[i+1]; ++j)
          {
@@ -2691,7 +2734,7 @@ void HypreParMatrix::Destroy()
    if (A == NULL) { return; }
 
 #ifdef HYPRE_USING_GPU
-   if (ParCSROwner && (diagOwner < 0 || offdOwner < 0))
+   if (HypreUsingGPU() && ParCSROwner && (diagOwner < 0 || offdOwner < 0))
    {
       // Put the "host" or "hypre" pointers in {i,j,data} of A->{diag,offd}, so
       // that they can be destroyed by hypre when hypre_ParCSRMatrixDestroy(A)
@@ -2855,10 +2898,15 @@ HypreParMatrix * ParMult(const HypreParMatrix *A, const HypreParMatrix *B,
 {
    hypre_ParCSRMatrix * ab;
 #ifdef HYPRE_USING_GPU
-   ab = hypre_ParCSRMatMat(*A, *B);
-#else
-   ab = hypre_ParMatmul(*A,*B);
+   if (HypreUsingGPU())
+   {
+      ab = hypre_ParCSRMatMat(*A, *B);
+   }
+   else
 #endif
+   {
+      ab = hypre_ParMatmul(*A,*B);
+   }
    hypre_ParCSRMatrixSetNumNonzeros(ab);
 
    if (!hypre_ParCSRMatrixCommPkg(ab)) { hypre_MatvecCommPkgCreate(ab); }
@@ -2882,6 +2930,7 @@ HypreParMatrix * RAP(const HypreParMatrix *A, const HypreParMatrix *P)
    //        in ex28p.
    // Quick fix: add a diagonal matrix with 0 diagonal.
    // Maybe use hypre_CSRMatrixCheckDiagFirst to see if we need the fix.
+   if (HypreUsingGPU())
    {
       hypre_ParCSRMatrix *Q = hypre_ParCSRMatMat(*A,*P);
       const bool keepTranspose = false;
@@ -2891,25 +2940,27 @@ HypreParMatrix * RAP(const HypreParMatrix *A, const HypreParMatrix *P)
       // alternative:
       // hypre_ParCSRMatrixRAPKT
    }
-#else
-#if MFEM_HYPRE_VERSION <= 22200
-   HYPRE_Int P_owns_its_col_starts =
-      hypre_ParCSRMatrixOwnsColStarts((hypre_ParCSRMatrix*)(*P));
+   else
 #endif
-
-   hypre_BoomerAMGBuildCoarseOperator(*P,*A,*P,&rap);
-
-#if MFEM_HYPRE_VERSION <= 22200
-   /* Warning: hypre_BoomerAMGBuildCoarseOperator steals the col_starts
-      from P (even if it does not own them)! */
-   hypre_ParCSRMatrixSetRowStartsOwner(rap,0);
-   hypre_ParCSRMatrixSetColStartsOwner(rap,0);
-   if (P_owns_its_col_starts)
    {
-      hypre_ParCSRMatrixSetColStartsOwner(*P, 1);
+#if MFEM_HYPRE_VERSION <= 22200
+      HYPRE_Int P_owns_its_col_starts =
+         hypre_ParCSRMatrixOwnsColStarts((hypre_ParCSRMatrix*)(*P));
+#endif
+
+      hypre_BoomerAMGBuildCoarseOperator(*P,*A,*P,&rap);
+
+#if MFEM_HYPRE_VERSION <= 22200
+      /* Warning: hypre_BoomerAMGBuildCoarseOperator steals the col_starts
+         from P (even if it does not own them)! */
+      hypre_ParCSRMatrixSetRowStartsOwner(rap,0);
+      hypre_ParCSRMatrixSetColStartsOwner(rap,0);
+      if (P_owns_its_col_starts)
+      {
+         hypre_ParCSRMatrixSetColStartsOwner(*P, 1);
+      }
+#endif
    }
-#endif
-#endif
 
    hypre_ParCSRMatrixSetNumNonzeros(rap);
    // hypre_MatvecCommPkgCreate(rap);
@@ -2923,36 +2974,39 @@ HypreParMatrix * RAP(const HypreParMatrix * Rt, const HypreParMatrix *A,
    hypre_ParCSRMatrix * rap;
 
 #ifdef HYPRE_USING_GPU
+   if (HypreUsingGPU())
    {
       hypre_ParCSRMatrix *Q = hypre_ParCSRMatMat(*A,*P);
       rap = hypre_ParCSRTMatMat(*Rt,Q);
       hypre_ParCSRMatrixDestroy(Q);
    }
-#else
+   else
+#endif
+   {
 #if MFEM_HYPRE_VERSION <= 22200
-   HYPRE_Int P_owns_its_col_starts =
-      hypre_ParCSRMatrixOwnsColStarts((hypre_ParCSRMatrix*)(*P));
-   HYPRE_Int Rt_owns_its_col_starts =
-      hypre_ParCSRMatrixOwnsColStarts((hypre_ParCSRMatrix*)(*Rt));
+      HYPRE_Int P_owns_its_col_starts =
+         hypre_ParCSRMatrixOwnsColStarts((hypre_ParCSRMatrix*)(*P));
+      HYPRE_Int Rt_owns_its_col_starts =
+         hypre_ParCSRMatrixOwnsColStarts((hypre_ParCSRMatrix*)(*Rt));
 #endif
 
-   hypre_BoomerAMGBuildCoarseOperator(*Rt,*A,*P,&rap);
+      hypre_BoomerAMGBuildCoarseOperator(*Rt,*A,*P,&rap);
 
 #if MFEM_HYPRE_VERSION <= 22200
-   /* Warning: hypre_BoomerAMGBuildCoarseOperator steals the col_starts
-      from Rt and P (even if they do not own them)! */
-   hypre_ParCSRMatrixSetRowStartsOwner(rap,0);
-   hypre_ParCSRMatrixSetColStartsOwner(rap,0);
-   if (P_owns_its_col_starts)
-   {
-      hypre_ParCSRMatrixSetColStartsOwner(*P, 1);
-   }
-   if (Rt_owns_its_col_starts)
-   {
-      hypre_ParCSRMatrixSetColStartsOwner(*Rt, 1);
-   }
-#endif
+      /* Warning: hypre_BoomerAMGBuildCoarseOperator steals the col_starts
+         from Rt and P (even if they do not own them)! */
+      hypre_ParCSRMatrixSetRowStartsOwner(rap,0);
+      hypre_ParCSRMatrixSetColStartsOwner(rap,0);
+      if (P_owns_its_col_starts)
+      {
+         hypre_ParCSRMatrixSetColStartsOwner(*P, 1);
+      }
+      if (Rt_owns_its_col_starts)
+      {
+         hypre_ParCSRMatrixSetColStartsOwner(*Rt, 1);
+      }
 #endif
+   }
 
    hypre_ParCSRMatrixSetNumNonzeros(rap);
    // hypre_MatvecCommPkgCreate(rap);
@@ -3382,7 +3436,7 @@ int ParCSRRelax_FIR(hypre_ParCSRMatrix *A, // matrix to relax with
 
 HypreSmoother::HypreSmoother() : Solver()
 {
-   type = default_type;
+   type = DefaultType();
    relax_times = 1;
    relax_weight = 1.0;
    omega = 1.0;
@@ -3407,7 +3461,7 @@ HypreSmoother::HypreSmoother(const HypreParMatrix &A_, int type_,
                              double omega_, int poly_order_,
                              double poly_fraction_, int eig_est_cg_iter_)
 {
-   type = type_;
+   type = (type_ == -1) ? DefaultType() : type_;
    relax_times = relax_times_;
    relax_weight = relax_weight_;
    omega = omega_;
@@ -3522,17 +3576,22 @@ void HypreSmoother::SetOperator(const Operator &op)
    if (l1_norms && pos_l1_norms)
    {
 #if defined(HYPRE_USING_GPU)
-      double *d_l1_norms = l1_norms;  // avoid *this capture
-      MFEM_GPU_FORALL(i, height,
+      if (HypreUsingGPU())
       {
-         d_l1_norms[i] = std::abs(d_l1_norms[i]);
-      });
-#else
-      for (int i = 0; i < height; i++)
-      {
-         l1_norms[i] = std::abs(l1_norms[i]);
+         double *d_l1_norms = l1_norms;  // avoid *this capture
+         MFEM_GPU_FORALL(i, height,
+         {
+            d_l1_norms[i] = std::abs(d_l1_norms[i]);
+         });
       }
+      else
 #endif
+      {
+         for (int i = 0; i < height; i++)
+         {
+            l1_norms[i] = std::abs(l1_norms[i]);
+         }
+      }
    }
 
    if (type == 16)
@@ -4808,41 +4867,49 @@ HypreBoomerAMG::HypreBoomerAMG(const HypreParMatrix &A) : HypreSolver(&A)
 
 void HypreBoomerAMG::SetDefaultOptions()
 {
-#if !defined(HYPRE_USING_GPU)
-   // AMG coarsening options:
-   int coarsen_type = 10;   // 10 = HMIS, 8 = PMIS, 6 = Falgout, 0 = CLJP
-   int agg_levels   = 1;    // number of aggressive coarsening levels
-   double theta     = 0.25; // strength threshold: 0.25, 0.5, 0.8
-
    // AMG interpolation options:
-   int interp_type  = 6;    // 6 = extended+i, 0 = classical
-   int Pmax         = 4;    // max number of elements per row in P
+   int coarsen_type, agg_levels, interp_type, Pmax, relax_type, relax_sweeps,
+       print_level, max_levels;
+   double theta;
 
-   // AMG relaxation options:
-   int relax_type   = 8;    // 8 = l1-GS, 6 = symm. GS, 3 = GS, 18 = l1-Jacobi
-   int relax_sweeps = 1;    // relaxation sweeps on each level
+   if (!HypreUsingGPU())
+   {
+      // AMG coarsening options:
+      coarsen_type = 10;   // 10 = HMIS, 8 = PMIS, 6 = Falgout, 0 = CLJP
+      agg_levels   = 1;    // number of aggressive coarsening levels
+      theta     = 0.25; // strength threshold: 0.25, 0.5, 0.8
 
-   // Additional options:
-   int print_level  = 1;    // print AMG iterations? 1 = no, 2 = yes
-   int max_levels   = 25;   // max number of levels in AMG hierarchy
-#else
-   // AMG coarsening options:
-   int coarsen_type = 8;    // 10 = HMIS, 8 = PMIS, 6 = Falgout, 0 = CLJP
-   int agg_levels   = 0;    // number of aggressive coarsening levels
-   double theta     = 0.25; // strength threshold: 0.25, 0.5, 0.8
+      // AMG interpolation options:
+      interp_type  = 6;    // 6 = extended+i, 0 = classical
+      Pmax         = 4;    // max number of elements per row in P
 
-   // AMG interpolation options:
-   int interp_type  = 6;    // 6 = extended+i, or 18 = extended+e
-   int Pmax         = 4;    // max number of elements per row in P
+      // AMG relaxation options:
+      relax_type   = 8;    // 8 = l1-GS, 6 = symm. GS, 3 = GS, 18 = l1-Jacobi
+      relax_sweeps = 1;    // relaxation sweeps on each level
 
-   // AMG relaxation options:
-   int relax_type   = 18;   // 18 = l1-Jacobi, or 16 = Chebyshev
-   int relax_sweeps = 1;    // relaxation sweeps on each level
+      // Additional options:
+      print_level  = 1;    // print AMG iterations? 1 = no, 2 = yes
+      max_levels   = 25;   // max number of levels in AMG hierarchy
+   }
+   else
+   {
+      // AMG coarsening options:
+      coarsen_type = 8;    // 10 = HMIS, 8 = PMIS, 6 = Falgout, 0 = CLJP
+      agg_levels   = 0;    // number of aggressive coarsening levels
+      theta     = 0.25; // strength threshold: 0.25, 0.5, 0.8
 
-   // Additional options:
-   int print_level  = 1;    // print AMG iterations? 1 = no, 2 = yes
-   int max_levels   = 25;   // max number of levels in AMG hierarchy
-#endif
+      // AMG interpolation options:
+      interp_type  = 6;    // 6 = extended+i, or 18 = extended+e
+      Pmax         = 4;    // max number of elements per row in P
+
+      // AMG relaxation options:
+      relax_type   = 18;   // 18 = l1-Jacobi, or 16 = Chebyshev
+      relax_sweeps = 1;    // relaxation sweeps on each level
+
+      // Additional options:
+      print_level  = 1;    // print AMG iterations? 1 = no, 2 = yes
+      max_levels   = 25;   // max number of levels in AMG hierarchy
+   }
 
    HYPRE_BoomerAMGSetCoarsenType(amg_precond, coarsen_type);
    HYPRE_BoomerAMGSetAggNumLevels(amg_precond, agg_levels);
@@ -4979,14 +5046,20 @@ void HypreBoomerAMG::SetSystemsOptions(int dim, bool order_bynodes)
       // After the addition of hypre_IntArray, mapping is assumed
       // to be a device pointer. Previously, it was assumed to be
       // a host pointer.
+      HYPRE_Int *mapping = nullptr;
 #if defined(hypre_IntArrayData) && defined(HYPRE_USING_GPU)
-      HYPRE_Int *mapping = mfem_hypre_CTAlloc(HYPRE_Int, height);
-      hypre_TMemcpy(mapping, h_mapping, HYPRE_Int, height,
-                    HYPRE_MEMORY_DEVICE, HYPRE_MEMORY_HOST);
-      mfem_hypre_TFree_host(h_mapping);
-#else
-      HYPRE_Int *mapping = h_mapping;
+      if (HypreUsingGPU())
+      {
+         mapping = mfem_hypre_CTAlloc(HYPRE_Int, height);
+         hypre_TMemcpy(mapping, h_mapping, HYPRE_Int, height,
+                       HYPRE_MEMORY_DEVICE, HYPRE_MEMORY_HOST);
+         mfem_hypre_TFree_host(h_mapping);
+      }
+      else
 #endif
+      {
+         mapping = h_mapping;
+      }
 
       // hypre actually deletes the mapping pointer in HYPRE_BoomerAMGDestroy,
       // so we don't need to track it
@@ -5078,7 +5151,10 @@ void HypreBoomerAMG::RecomputeRBMs()
 void HypreBoomerAMG::SetElasticityOptions(ParFiniteElementSpace *fespace_)
 {
 #ifdef HYPRE_USING_GPU
-   MFEM_ABORT("this method is not supported in hypre built with GPU support");
+   if (HypreUsingGPU())
+   {
+      MFEM_ABORT("this method is not supported in hypre built with GPU support");
+   }
 #endif
 
    // Save the finite element space to support multiple calls to SetOperator()
@@ -5314,23 +5390,14 @@ void HypreAMS::MakeSolver(int sdim, int cycle_type)
    int rlx_sweeps       = 1;
    double rlx_weight    = 1.0;
    double rlx_omega     = 1.0;
-#if !defined(HYPRE_USING_GPU)
-   int amg_coarsen_type = 10;
-   int amg_agg_levels   = 1;
-   int amg_rlx_type     = 8;
-   int rlx_type         = 2;
+   const bool hypre_gpu = HypreUsingGPU();
+   int amg_coarsen_type = hypre_gpu ? 8 : 10;
+   int amg_agg_levels   = hypre_gpu ? 0 : 1;
+   int amg_rlx_type     = hypre_gpu ? 18 : 8;
+   int rlx_type         = hypre_gpu ? 1: 2;
    double theta         = 0.25;
    int amg_interp_type  = 6;
    int amg_Pmax         = 4;
-#else
-   int amg_coarsen_type = 8;
-   int amg_agg_levels   = 0;
-   int amg_rlx_type     = 18;
-   int rlx_type         = 1;
-   double theta         = 0.25;
-   int amg_interp_type  = 6;
-   int amg_Pmax         = 4;
-#endif
 
    space_dim = sdim;
    ams_cycle_type = cycle_type;
@@ -5691,23 +5758,14 @@ void HypreADS::MakeSolver()
    int rlx_sweeps       = 1;
    double rlx_weight    = 1.0;
    double rlx_omega     = 1.0;
-#if !defined(HYPRE_USING_GPU)
-   int rlx_type         = 2;
-   int amg_coarsen_type = 10;
-   int amg_agg_levels   = 1;
-   int amg_rlx_type     = 8;
+   const bool hypre_gpu = HypreUsingGPU();
+   int rlx_type         = hypre_gpu ? 1 : 2;
+   int amg_coarsen_type = hypre_gpu ? 8 : 10;
+   int amg_agg_levels   = hypre_gpu ? 0 : 1;
+   int amg_rlx_type     = hypre_gpu ? 18 : 8;
    double theta         = 0.25;
    int amg_interp_type  = 6;
    int amg_Pmax         = 4;
-#else
-   int rlx_type         = 1;
-   int amg_coarsen_type = 8;
-   int amg_agg_levels   = 0;
-   int amg_rlx_type     = 18;
-   double theta         = 0.25;
-   int amg_interp_type  = 6;
-   int amg_Pmax         = 4;
-#endif
 
    HYPRE_ADSCreate(&ads);
 

linalg/hypre.hpp

@@ -46,15 +46,15 @@
 // MFEM_HYPRE_FORALL is a macro similar to mfem::forall, but it executes on the
 // device that hypre was configured with (no matter what device was selected
 // in MFEM's runtime configuration).
-#if defined(HYPRE_USING_CUDA)
-#define MFEM_HYPRE_FORALL(i, N,...) CuWrap1D(N, [=] MFEM_DEVICE      \
-                                       (int i) {__VA_ARGS__})
-#elif defined(HYPRE_USING_HIP)
-#define MFEM_HYPRE_FORALL(i, N,...) HipWrap1D(N, [=] MFEM_DEVICE     \
-                                        (int i) {__VA_ARGS__})
-#else
-#define MFEM_HYPRE_FORALL(i, N,...) for (int i = 0; i < N; i++) { __VA_ARGS__ }
-#endif
+//#if defined(HYPRE_USING_CUDA)
+//#define MFEM_HYPRE_FORALL(i, N,...) CuWrap1D(N, [=] MFEM_DEVICE      \
+//                                       (int i) {__VA_ARGS__})
+//#elif defined(HYPRE_USING_HIP)
+//#define MFEM_HYPRE_FORALL(i, N,...) HipWrap1D(N, [=] MFEM_DEVICE     \
+//                                        (int i) {__VA_ARGS__})
+//#else
+//#define MFEM_HYPRE_FORALL(i, N,...) for (int i = 0; i < N; i++) { __VA_ARGS__ }
+//#endif
 
 #include "sparsemat.hpp"
 #include "hypre_parcsr.hpp"
@@ -134,14 +134,15 @@ inline int to_int(HYPRE_Int i)
 
 
 /// The MemoryClass used by Hypre objects.
-inline constexpr MemoryClass GetHypreMemoryClass()
+inline MemoryClass GetHypreMemoryClass()
 {
 #if !defined(HYPRE_USING_GPU)
    return MemoryClass::HOST;
 #elif defined(HYPRE_USING_UNIFIED_MEMORY)
    return MemoryClass::MANAGED;
 #else
-   return MemoryClass::DEVICE;
+   return (GetHypreMemoryLocation() == HYPRE_MEMORY_DEVICE) ? MemoryClass::DEVICE :
+          MemoryClass::HOST;
 #endif
 }
 
@@ -153,7 +154,19 @@ inline MemoryType GetHypreMemoryType()
 #elif defined(HYPRE_USING_UNIFIED_MEMORY)
    return MemoryType::MANAGED;
 #else
-   return MemoryType::DEVICE;
+   return (GetHypreMemoryLocation() == HYPRE_MEMORY_DEVICE) ? MemoryType::DEVICE :
+          Device::GetHostMemoryType();
+#endif
+}
+
+inline bool HypreUsingGPU()
+{
+#ifdef HYPRE_USING_GPU
+   HYPRE_MemoryLocation loc;
+   HYPRE_GetMemoryLocation(&loc);
+   return loc == HYPRE_MEMORY_DEVICE;
+#else
+   return false;
 #endif
 }
 
@@ -1026,6 +1039,7 @@ protected:
 
 public:
    /** Hypre smoother types:
+      -1    = Undefined, replaced with DefaultType() in constructors
        0    = Jacobi
        1    = l1-scaled Jacobi
        2    = l1-scaled block Gauss-Seidel/SSOR
@@ -1036,18 +1050,18 @@ public:
        16   = Chebyshev
        1001 = Taubin polynomial smoother
        1002 = FIR polynomial smoother. */
-   enum Type { Jacobi = 0, l1Jacobi = 1, l1GS = 2, l1GStr = 4, lumpedJacobi = 5,
+   enum Type { Undefined = -1, Jacobi = 0, l1Jacobi = 1, l1GS = 2, l1GStr = 4, lumpedJacobi = 5,
                GS = 6, OPFS = 10, Chebyshev = 16, Taubin = 1001, FIR = 1002
              };
-#if !defined(HYPRE_USING_GPU)
-   static constexpr Type default_type = l1GS;
-#else
-   static constexpr Type default_type = l1Jacobi;
-#endif
+
+   Type DefaultType()
+   {
+      return HypreUsingGPU() ? l1Jacobi : l1GS;
+   }
 
    HypreSmoother();
 
-   HypreSmoother(const HypreParMatrix &A_, int type = default_type,
+   HypreSmoother(const HypreParMatrix &A_, int type = Undefined,
                  int relax_times = 1, double relax_weight = 1.0,
                  double omega = 1.0, int poly_order = 2,
                  double poly_fraction = .3, int eig_est_cg_iter = 10);

linalg/hypre_parcsr.hpp

@@ -29,6 +29,34 @@ typedef HYPRE_Int HYPRE_BigInt;
 #define HYPRE_MPI_BIG_INT HYPRE_MPI_INT
 #endif
 
+
+namespace mfem
+{
+#if MFEM_HYPRE_VERSION >= 22600
+inline HYPRE_MemoryLocation GetHypreMemoryLocation()
+{
+   HYPRE_MemoryLocation loc;
+   HYPRE_GetMemoryLocation(&loc);
+   return loc;
+}
+inline HYPRE_ExecutionPolicy GetHypreExecutionPolicy()
+{
+   HYPRE_ExecutionPolicy pol;
+   HYPRE_GetExecutionPolicy(&pol);
+   return pol;
+}
+#else
+inline HYPRE_MemoryLocation GetHypreMemoryLocation()
+{
+#ifdef HYPRE_USING_GPU
+   return HYPRE_MEMORY_DEVICE;
+#else
+   return HYPRE_MEMORY_HOST;
+#endif // HYPRE_USING_GPU
+}
+#endif
+};
+
 // Define macro wrappers for hypre_TAlloc, hypre_CTAlloc and hypre_TFree:
 // mfem_hypre_TAlloc, mfem_hypre_CTAlloc, and mfem_hypre_TFree, respectively.
 // Note: these macros are used in hypre.cpp, hypre_parcsr.cpp, and perhaps
@@ -46,10 +74,10 @@ typedef HYPRE_Int HYPRE_BigInt;
 #else // MFEM_HYPRE_VERSION >= 21400
 
 #define mfem_hypre_TAlloc(type, size) \
-   hypre_TAlloc(type, size, HYPRE_MEMORY_DEVICE)
+   hypre_TAlloc(type, size, mfem::GetHypreMemoryLocation())
 #define mfem_hypre_CTAlloc(type, size) \
-   hypre_CTAlloc(type, size, HYPRE_MEMORY_DEVICE)
-#define mfem_hypre_TFree(ptr) hypre_TFree(ptr, HYPRE_MEMORY_DEVICE)
+   hypre_CTAlloc(type, size, mfem::GetHypreMemoryLocation())
+#define mfem_hypre_TFree(ptr) hypre_TFree(ptr, mfem::GetHypreMemoryLocation())
 
 #define mfem_hypre_TAlloc_host(type, size) \
    hypre_TAlloc(type, size, HYPRE_MEMORY_HOST)

mesh/mesh.cpp

@@ -6271,111 +6271,6 @@ int Mesh::CheckBdrElementOrientation(bool fix_it)
    return wo;
 }
 
-IntegrationPoint Mesh::TransformBdrElementToFace(Geometry::Type geom, int o,
-                                                 const IntegrationPoint &ip)
-{
-   IntegrationPoint fip = ip;
-   if (geom == Geometry::POINT)
-   {
-      return fip;
-   }
-   else if (geom == Geometry::SEGMENT)
-   {
-      MFEM_ASSERT(o >= 0 && o < 2, "Invalid orientation for Geometry::SEGMENT!");
-      if (o == 0)
-      {
-         fip.x = ip.x;
-      }
-      else if (o == 1)
-      {
-         fip.x = 1.0 - ip.x;
-      }
-   }
-   else if (geom == Geometry::TRIANGLE)
-   {
-      MFEM_ASSERT(o >= 0 && o < 6, "Invalid orientation for Geometry::TRIANGLE!");
-      if (o == 0)  // 0, 1, 2
-      {
-         fip.x = ip.x;
-         fip.y = ip.y;
-      }
-      else if (o == 5)  // 0, 2, 1
-      {
-         fip.x = ip.y;
-         fip.y = ip.x;
-      }
-      else if (o == 2)  // 1, 2, 0
-      {
-         fip.x = 1.0 - ip.x - ip.y;
-         fip.y = ip.x;
-      }
-      else if (o == 1)  // 1, 0, 2
-      {
-         fip.x = 1.0 - ip.x - ip.y;
-         fip.y = ip.y;
-      }
-      else if (o == 4)  // 2, 0, 1
-      {
-         fip.x = ip.y;
-         fip.y = 1.0 - ip.x - ip.y;
-      }
-      else if (o == 3)  // 2, 1, 0
-      {
-         fip.x = ip.x;
-         fip.y = 1.0 - ip.x - ip.y;
-      }
-   }
-   else if (geom == Geometry::SQUARE)
-   {
-      MFEM_ASSERT(o >= 0 && o < 8, "Invalid orientation for Geometry::SQUARE!");
-      if (o == 0)  // 0, 1, 2, 3
-      {
-         fip.x = ip.x;
-         fip.y = ip.y;
-      }
-      else if (o == 1)  // 0, 3, 2, 1
-      {
-         fip.x = ip.y;
-         fip.y = ip.x;
-      }
-      else if (o == 2)  // 1, 2, 3, 0
-      {
-         fip.x = ip.y;
-         fip.y = 1.0 - ip.x;
-      }
-      else if (o == 3)  // 1, 0, 3, 2
-      {
-         fip.x = 1.0 - ip.x;
-         fip.y = ip.y;
-      }
-      else if (o == 4)  // 2, 3, 0, 1
-      {
-         fip.x = 1.0 - ip.x;
-         fip.y = 1.0 - ip.y;
-      }
-      else if (o == 5)  // 2, 1, 0, 3
-      {
-         fip.x = 1.0 - ip.y;
-         fip.y = 1.0 - ip.x;
-      }
-      else if (o == 6)  // 3, 0, 1, 2
-      {
-         fip.x = 1.0 - ip.y;
-         fip.y = ip.x;
-      }
-      else if (o == 7)  // 3, 2, 1, 0
-      {
-         fip.x = ip.x;
-         fip.y = 1.0 - ip.y;
-      }
-   }
-   else
-   {
-      MFEM_ABORT("Unsupported face geometry for TransformBdrElementToFace!");
-   }
-   return fip;
-}
-
 int Mesh::GetNumGeometries(int dim) const
 {
    MFEM_ASSERT(0 <= dim && dim <= Dim, "invalid dim: " << dim);
@@ -6669,20 +6564,24 @@ Array<int> Mesh::FindFaceNeighbors(const int elem) const
 
 void Mesh::GetBdrElementFace(int i, int *f, int *o) const
 {
-   *f = GetBdrElementEdgeIndex(i);
+   const int *bv, *fv;
 
-   const int *fv = (Dim > 1) ? faces[*f]->GetVertices() : NULL;
-   const int *bv = boundary[i]->GetVertices();
+   *f = be_to_face[i];
+   bv = boundary[i]->GetVertices();
+   fv = faces[be_to_face[i]]->GetVertices();
 
    // find the orientation of the bdr. elem. w.r.t.
    // the corresponding face element (that's the base)
-   switch (GetBdrElementGeometry(i))
+   switch (GetBdrElementType(i))
    {
-      case Geometry::POINT:    *o = 0; break;
-      case Geometry::SEGMENT:  *o = (fv[0] == bv[0]) ? 0 : 1; break;
-      case Geometry::TRIANGLE: *o = GetTriOrientation(fv, bv); break;
-      case Geometry::SQUARE:   *o = GetQuadOrientation(fv, bv); break;
-      default: MFEM_ABORT("invalid geometry");
+      case Element::TRIANGLE:
+         *o = GetTriOrientation(fv, bv);
+         break;
+      case Element::QUADRILATERAL:
+         *o = GetQuadOrientation(fv, bv);
+         break;
+      default:
+         MFEM_ABORT("invalid geometry");
    }
 }
 
@@ -7284,7 +7183,7 @@ STable3D *Mesh::GetFacesTable()
             break;
          }
          default:
-            MFEM_ABORT("Unexpected type of Element: " << GetElementType(i));
+            MFEM_ABORT("Unexpected type of Element.");
       }
    }
    return faces_tbl;

mesh/mesh.hpp

@@ -322,38 +322,9 @@ protected:
    void ReadNURBSMesh(std::istream &input, int &curved, int &read_gf);
    void ReadInlineMesh(std::istream &input, bool generate_edges = false);
    void ReadGmshMesh(std::istream &input, int &curved, int &read_gf);
-
    /* Note NetCDF (optional library) is used for reading cubit files */
 #ifdef MFEM_USE_NETCDF
-
-   /// @brief Load a mesh from a Genesis file.
    void ReadCubit(const char *filename, int &curved, int &read_gf);
-
-   /// @brief The final step in constructing the mesh from a Genesis file. This
-   /// is only called if the mesh order == 2 (determined internally from the
-   /// cubit element type).
-   void FinalizeCubitSecondOrderMesh(const int cubit_element_type,
-                                     const int num_element_blocks,
-                                     const int num_nodes_per_element,
-                                     const int *start_of_block,
-                                     const double *coordx,
-                                     const double *coordy,
-                                     const double *coordz,
-                                     const int **element_blocks);
-
-   /// @brief Returns a pointer to a new mfem::Element based on the provided
-   /// cubit element type. This is used internally to create the mesh elements
-   /// from a Genesis file.
-   Element *CreateCubitElement(const int cubit_element_type,
-                               const int *vertex_ids,
-                               const int block_id);
-
-   /// @brief Returns a pointer to a new mfem::Element based on the provided
-   /// cubit face type. This is used internally to create the boundary elements
-   /// from a Genesis file.
-   Element *CreateCubitBoundaryElement(const int cubit_face_type,
-                                       const int *vertex_ids,
-                                       const int sideset_id) const;
 #endif
 
    /// Determine the mesh generator bitmask #meshgen, see MeshGenerator().
@@ -1359,12 +1330,7 @@ public:
        element @a elem, including @a elem. */
    Array<int> FindFaceNeighbors(const int elem) const;
 
-   /** Return the index and the orientation of the vertex of bdr element i. (1D)
-       Return the index and the orientation of the edge of bdr element i. (2D)
-       Return the index and the orientation of the face of bdr element i. (3D)
-
-       In 2D, the returned edge orientation is 0 or 1, not +/-1 as returned by
-       GetElementEdges/GetBdrElementEdges. */
+   /// Return the index and the orientation of the face of bdr element i. (3D)
    void GetBdrElementFace(int i, int *f, int *o) const;
 
    /** Return the vertex index of boundary element i. (1D)
@@ -1443,14 +1409,6 @@ public:
    /// should not be deleted by the caller.
    static FiniteElement *GetTransformationFEforElementType(Element::Type);
 
-   /** @brief For the vertex (1D), edge (2D), or face (3D) of a boundary element
-       with the orientation @a o, return the transformation of the boundary
-       element integration point @ ip to the face element. In 2D, the
-       the orientation is 0 or 1 as returned by GetBdrElementFace, not +/-1.
-       Supports both internal and external boundaries. */
-   static IntegrationPoint TransformBdrElementToFace(Geometry::Type geom, int o,
-                                                     const IntegrationPoint &ip);
-
    /// @anchor mfem_Mesh_elem_trans
    /// @name Access the coordinate transformation for individual elements
    ///

mesh/ncmesh.cpp

@@ -3021,13 +3021,13 @@ void NCMesh::TraverseTetEdge(int vn0, int vn1, const Point &p0, const Point &p1,
    if (nd.HasEdge())
    {
       // check if the edge is already a master in 'edge_list'
-      const auto eid_and_type = edge_list.GetMeshIdAndType(nd.edge_index);
-      if (eid_and_type.type == NCList::MeshIdType::MASTER
-          || eid_and_type.type == NCList::MeshIdType::CONFORMING)
+      int type;
+      const MeshId &eid = edge_list.LookUp(nd.edge_index, &type);
+      if (type == 1)
       {
          // in this case we need to add an edge-face constraint, because the
-         // non-slave edge is really a (face-)slave itself.
-         const MeshId &eid = *eid_and_type.id;
+         // master edge is really a (face-)slave itself
+
          face_list.slaves.Append(
             Slave(-1 - eid.index, eid.element, eid.local, Geometry::TRIANGLE));
 
@@ -3048,10 +3048,9 @@ void NCMesh::TraverseTetEdge(int vn0, int vn1, const Point &p0, const Point &p1,
    TraverseTetEdge(mid, vn1, pmid, p1, matrix_map);
 }
 
-NCMesh::TriFaceTraverseResults NCMesh::TraverseTriFace(int vn0, int vn1,
-                                                       int vn2,
-                                                       const PointMatrix& pm, int level,
-                                                       MatrixMap &matrix_map)
+bool NCMesh::TraverseTriFace(int vn0, int vn1, int vn2,
+                             const PointMatrix& pm, int level,
+                             MatrixMap &matrix_map)
 {
    if (level > 0)
    {
@@ -3070,7 +3069,7 @@ NCMesh::TriFaceTraverseResults NCMesh::TraverseTriFace(int vn0, int vn1,
          sl.local = ReorderFacePointMat(vn0, vn1, vn2, -1, elem, pm, pm_r);
          sl.matrix = matrix_map.GetIndex(pm_r);
 
-         return {true, elements[elem].rank != MyRank};
+         return true;
       }
    }
 
@@ -3078,7 +3077,7 @@ NCMesh::TriFaceTraverseResults NCMesh::TraverseTriFace(int vn0, int vn1,
    if (TriFaceSplit(vn0, vn1, vn2, mid))
    {
       Point pmid0(pm(0), pm(1)), pmid1(pm(1), pm(2)), pmid2(pm(2), pm(0));
-      TriFaceTraverseResults b[4];
+      bool b[4];
 
       b[0] = TraverseTriFace(vn0, mid[0], mid[2],
                              PointMatrix(pm(0), pmid0, pmid2),
@@ -3096,21 +3095,16 @@ NCMesh::TriFaceTraverseResults NCMesh::TraverseTriFace(int vn0, int vn1,
                              PointMatrix(pmid1, pmid2, pmid0),
                              level+1, matrix_map);
 
-      // Traverse possible tet edges constrained by the master face. This needs to occur if
-      // none of these first NC level faces are split further, OR if they are on different
-      // processors. The different processor constraint is needed in the case of local
-      // elements constrained by this face via the edge alone. Cannot know this a priori, so
-      // just constrain any edge attached to two neighbors.
-      if (HaveTets() && (!b[3].unsplit || b[3].ghost_neighbor))
+      // traverse possible tet edges constrained by the master face
+      if (HaveTets() && !b[3])
       {
-         // If the faces have no further splits, so would not be captured by normal face
-         // relations, add possible edge constraints.
-         if (!b[1].unsplit || b[1].ghost_neighbor) { TraverseTetEdge(mid[0],mid[1], pmid0,pmid1, matrix_map); }
-         if (!b[2].unsplit || b[2].ghost_neighbor) { TraverseTetEdge(mid[1],mid[2], pmid1,pmid2, matrix_map); }
-         if (!b[0].unsplit || b[0].ghost_neighbor) { TraverseTetEdge(mid[2],mid[0], pmid2,pmid0, matrix_map); }
+         if (!b[1]) { TraverseTetEdge(mid[0],mid[1], pmid0,pmid1, matrix_map); }
+         if (!b[2]) { TraverseTetEdge(mid[1],mid[2], pmid1,pmid2, matrix_map); }
+         if (!b[0]) { TraverseTetEdge(mid[2],mid[0], pmid2,pmid0, matrix_map); }
       }
    }
-   return {false, false};
+
+   return false;
 }
 
 void NCMesh::BuildFaceList()
@@ -3408,80 +3402,77 @@ void NCMesh::NCList::Clear()
       point_matrices[i].DeleteAll();
    }
 
-   inv_index.clear();
+   inv_index.DeleteAll();
 }
 
-NCMesh::NCList::MeshIdAndType
-NCMesh::NCList::GetMeshIdAndType(int index) const
+long NCMesh::NCList::TotalSize() const
 {
-   BuildIndex();
-   const auto it = inv_index.find(index);
-   auto ft = it != inv_index.end() ? it->second.first : MeshIdType::UNRECOGNIZED;
-   switch (ft)
+   return conforming.Size() + masters.Size() + slaves.Size();
+}
+
+const NCMesh::MeshId& NCMesh::NCList::LookUp(int index, int *type) const
+{
+   if (!inv_index.Size())
    {
-      case MeshIdType::CONFORMING:
-         return {&conforming[it->second.second], it->second.first};
-      case MeshIdType::MASTER:
-         return {&masters[it->second.second], it->second.first};
-      case MeshIdType::SLAVE:
-         return {&slaves[it->second.second], it->second.first};
-      case MeshIdType::UNRECOGNIZED:
-      default:
-         return {nullptr, MeshIdType::UNRECOGNIZED};
-   }
-}
-
-NCMesh::NCList::MeshIdType
-NCMesh::NCList::GetMeshIdType(int index) const
-{
-   BuildIndex();
-   auto it = inv_index.find(index);
-   return (it != inv_index.end()) ? it->second.first : MeshIdType::UNRECOGNIZED;
-}
-
-bool
-NCMesh::NCList::CheckMeshIdType(int index, MeshIdType ft) const
-{
-   return GetMeshIdType(index) == ft;
-}
-
-void
-NCMesh::NCList::BuildIndex() const
-{
-   if (inv_index.size() == 0)
-   {
-      auto index_compare = [](const MeshId &a, const MeshId &b) { return a.index < b.index; };
-      auto max_conforming = std::max_element(conforming.begin(), conforming.end(),
-                                             index_compare);
-      auto max_master = std::max_element(masters.begin(), masters.end(),
-                                         index_compare);
-      auto max_slave = std::max_element(slaves.begin(), slaves.end(), index_compare);
-
-      int max_conforming_index = max_conforming != nullptr ? max_conforming->index :
-                                 -1;
-      int max_master_index = max_master != nullptr ? max_master->index : -1;
-      int max_slave_index = max_slave != nullptr ? max_slave->index : -1;
-
-      inv_index.reserve(std::max({max_conforming_index, max_master_index, max_slave_index}));
+      int max_index = -1;
       for (int i = 0; i < conforming.Size(); i++)
       {
-         inv_index.emplace(conforming[i].index, std::make_pair(MeshIdType::CONFORMING,
-                                                               i));
+         max_index = std::max(conforming[i].index, max_index);
       }
       for (int i = 0; i < masters.Size(); i++)
       {
-         inv_index.emplace(masters[i].index, std::make_pair(MeshIdType::MASTER, i));
+         max_index = std::max(masters[i].index, max_index);
       }
       for (int i = 0; i < slaves.Size(); i++)
       {
-         inv_index.emplace(slaves[i].index, std::make_pair(MeshIdType::SLAVE, i));
+         if (slaves[i].index < 0) { continue; }
+         max_index = std::max(slaves[i].index, max_index);
+      }
+
+      inv_index.SetSize(max_index + 1);
+      inv_index = -1;
+
+      for (int i = 0; i < conforming.Size(); i++)
+      {
+         inv_index[conforming[i].index] = (i << 2);
+      }
+      for (int i = 0; i < masters.Size(); i++)
+      {
+         inv_index[masters[i].index] = (i << 2) + 1;
+      }
+      for (int i = 0; i < slaves.Size(); i++)
+      {
+         if (slaves[i].index < 0) { continue; }
+         inv_index[slaves[i].index] = (i << 2) + 2;
       }
    }
 
-   MFEM_ASSERT(inv_index.size() > 0,
-               "Empty inverse index, member lists must be populated before BuildIndex is called!");
+   MFEM_ASSERT(index >= 0 && index < inv_index.Size(), "");
+   int key = inv_index[index];
+
+   if (!type)
+   {
+      MFEM_VERIFY(key >= 0, "index " << index << " not found.");
+   }
+   else // return entity type if requested, don't abort when not found
+   {
+      *type = (key >= 0) ? (key & 0x3) : -1;
+
+      static MeshId invalid;
+      if (*type < 0) { return invalid; } // not found
+   }
+
+   // return found entity MeshId
+   switch (key & 0x3)
+   {
+      case 0: return conforming[key >> 2];
+      case 1: return masters[key >> 2];
+      case 2: return slaves[key >> 2];
+      default: MFEM_ABORT("internal error"); return conforming[0];
+   }
 }
 
+
 //// Neighbors /////////////////////////////////////////////////////////////////
 
 void NCMesh::CollectEdgeVertices(int v0, int v1, Array<int> &indices)

mesh/ncmesh.hpp

@@ -24,7 +24,6 @@
 #include <vector>
 #include <map>
 #include <iostream>
-#include <unordered_map>
 
 namespace mfem
 {
@@ -224,73 +223,30 @@ public:
          , master(-1), matrix(0), edge_flags(0) {}
    };
 
-
    /// Lists all edges/faces in the nonconforming mesh.
    struct NCList
    {
-      Array<MeshId> conforming; ///< All MeshIds corresponding to conformal faces
-      Array<Master> masters; ///< All MeshIds corresponding to master faces
-      Array<Slave> slaves; ///< All MeshIds corresponding to slave faces
+      Array<MeshId> conforming;
+      Array<Master> masters;
+      Array<Slave> slaves;
 
       /// List of unique point matrices for each slave geometry.
       Array<DenseMatrix*> point_matrices[Geometry::NumGeom];
 
+      /// Return the point matrix oriented according to the master and slave edges
       void OrientedPointMatrix(const Slave &slave,
                                DenseMatrix &oriented_matrix) const;
 
-      /// Particular MeshId type, used for allowing static casting to the
-      /// appropriate child type after searching the NCList. UNRECOGNIZED
-      /// denotes that an instance is not known within the NCList, meaning that
-      /// it does not play a part in NC mechanics. This can be because the index
-      /// did not exist in the original Mesh, or because the entry is a boundary
-      /// face, whose NC status is always conforming.
-      enum class MeshIdType : char {CONFORMING, MASTER, SLAVE, UNRECOGNIZED};
-
-      /// Helper storing a reference to a MeshId type, and the face type it can
-      /// be cast to
-      struct MeshIdAndType
-      {
-         const MeshId * const id; ///< Pointer to a possible MeshId, nullptr if not found
-         /// MeshIdType corresponding to the MeshId. UNRECOGNIZED if unfound.
-         const MeshIdType type;
-      };
-      /// Return a mesh id and type for a given nc index.
-      MeshIdAndType GetMeshIdAndType(int index) const;
-
-      /// Return a face type for a given nc index.
-      MeshIdType GetMeshIdType(int index) const;
-
-      /// Given an index, check if this is a certain face type.
-      bool CheckMeshIdType(int index, MeshIdType type) const;
-
-      /// Erase the contents of the conforming, master and slave arrays.
       void Clear();
-      /// Whether the NCList is empty.
-      bool Empty() const
-      {
-         return conforming.Size() == 0
-                && masters.Size() == 0
-                && slaves.Size() == 0;
-      }
-      /// The total size of the component arrays in the NCList.
-      long TotalSize() const
-      {
-         return conforming.Size() + masters.Size() + slaves.Size();
-      }
-      /// The memory usage of the three public arrays. Does not account for the
-      /// inverse index.
+      bool Empty() const { return !conforming.Size() && !masters.Size(); }
+      long TotalSize() const;
       long MemoryUsage() const;
+
+      const MeshId& LookUp(int index, int *type = NULL) const;
+
       ~NCList() { Clear(); }
    private:
-      // Check for existence or construct the inv_index list map if necessary.
-      // const because only modifies the mutable member inv_index.
-      void BuildIndex() const;
-
-      /// A lazily constructed map from index to MeshId. Built whenever
-      /// GetMeshIdAndType, GetMeshIdType or CheckMeshIdType is called for the
-      /// first time. The MeshIdType is stored with, to enable casting to Slave
-      /// or Master elements appropriately.
-      mutable std::unordered_map<int, std::pair<MeshIdType, int>> inv_index;
+      mutable Array<int> inv_index;
    };
 
    /// Return the current list of conforming and nonconforming faces.
@@ -771,14 +727,9 @@ protected: // implementation
    void TraverseQuadFace(int vn0, int vn1, int vn2, int vn3,
                          const PointMatrix& pm, int level, Face* eface[4],
                          MatrixMap &matrix_map);
-   struct TriFaceTraverseResults
-   {
-      bool unsplit; ///< Whether this face has no further splits.
-      bool ghost_neighbor; ///< Whether the face neighbor is a ghost.
-   };
-   TriFaceTraverseResults TraverseTriFace(int vn0, int vn1, int vn2,
-                                          const PointMatrix& pm, int level,
-                                          MatrixMap &matrix_map);
+   bool TraverseTriFace(int vn0, int vn1, int vn2,
+                        const PointMatrix& pm, int level,
+                        MatrixMap &matrix_map);
    void TraverseTetEdge(int vn0, int vn1, const Point &p0, const Point &p1,
                         MatrixMap &matrix_map);
    void TraverseEdge(int vn0, int vn1, double t0, double t1, int flags,

mesh/nurbs.cpp

@@ -39,7 +39,6 @@ KnotVector::KnotVector(int Order_, int NCP)
    Order = Order_;
    NumOfControlPoints = NCP;
    knot.SetSize(NumOfControlPoints + Order + 1);
-   NumOfElements = 0;
 
    knot = -1.;
 }
@@ -130,35 +129,27 @@ void KnotVector::Print(std::ostream &os) const
    knot.Print(os, knot.Size());
 }
 
+
 void KnotVector::PrintFunctions(std::ostream &os, int samples) const
 {
-   MFEM_VERIFY(GetNE(), "Elements not counted. Use GetElements().");
-
    Vector shape(Order+1);
 
    double x, dx = 1.0/double (samples - 1);
 
-   /* @a cnt is a counter including elements between repeated knots if
-      present. This is required for usage of CalcShape. */
-   int cnt = 0;
-
-   for (int e = 0; e < GetNE(); e++, cnt++)
+   for (int i = 0; i <GetNE() ; i++)
    {
-      // Avoid printing shapes between repeated knots
-      if (!isElement(cnt)) { e--; continue; }
-
       for (int j = 0; j <samples; j++)
       {
          x =j*dx;
-         os<< x + e;
+         os<< x + i;
 
-         CalcShape ( shape, cnt, x);
+         CalcShape ( shape, i, x);
          for (int d = 0; d < Order+1; d++) { os<<"\t"<<shape[d]; }
 
-         CalcDShape ( shape, cnt, x);
+         CalcDShape ( shape, i, x);
          for (int d = 0; d < Order+1; d++) { os<<"\t"<<shape[d]; }
 
-         CalcD2Shape ( shape, cnt, x);
+         CalcD2Shape ( shape, i, x);
          for (int d = 0; d < Order+1; d++) { os<<"\t"<<shape[d]; }
          os<<endl;
       }

mesh/nurbs.hpp

@@ -88,10 +88,6 @@ public:
 
    void Print(std::ostream &out) const;
 
-   /** Prints the non-zero shape functions and their first and second
-       derivatives associated with the KnotVector per element. Use GetElements()
-       to count the elements before using this function. @a samples is the
-       number of samples of the shape functions per element.*/
    void PrintFunctions(std::ostream &out, int samples=11) const;
 
    /// Destroys KnotVector

mesh/pmesh.cpp

@@ -34,6 +34,7 @@ ParMesh::ParMesh(const ParMesh &pmesh, bool copy_nodes)
      group_sedge(pmesh.group_sedge),
      group_stria(pmesh.group_stria),
      group_squad(pmesh.group_squad),
+     face_nbr_el_to_face(NULL),
      glob_elem_offset(-1),
      glob_offset_sequence(-1),
      gtopo(pmesh.gtopo)
@@ -105,7 +106,8 @@ ParMesh& ParMesh::operator=(ParMesh &&mesh)
 
 ParMesh::ParMesh(MPI_Comm comm, Mesh &mesh, int *partitioning_,
                  int part_method)
-   : glob_elem_offset(-1)
+   : face_nbr_el_to_face(NULL)
+   , glob_elem_offset(-1)
    , glob_offset_sequence(-1)
    , gtopo(comm)
 {
@@ -852,6 +854,7 @@ ParMesh::ParMesh(const ParNCMesh &pncmesh)
    : MyComm(pncmesh.MyComm)
    , NRanks(pncmesh.NRanks)
    , MyRank(pncmesh.MyRank)
+   , face_nbr_el_to_face(NULL)
    , glob_elem_offset(-1)
    , glob_offset_sequence(-1)
    , gtopo(MyComm)
@@ -903,7 +906,7 @@ void ParMesh::FinalizeParTopo()
    sface_lface.SetSize(nst + shared_quads.Size());
    if (sface_lface.Size())
    {
-      auto faces_tbl = std::unique_ptr<STable3D>(GetFacesTable());
+      STable3D *faces_tbl = GetFacesTable();
       for (int st = 0; st < nst; st++)
       {
          const int *v = shared_trias[st].v;
@@ -914,11 +917,13 @@ void ParMesh::FinalizeParTopo()
          const int *v = shared_quads[sq].v;
          sface_lface[nst+sq] = (*faces_tbl)(v[0], v[1], v[2], v[3]);
       }
+      delete faces_tbl;
    }
 }
 
 ParMesh::ParMesh(MPI_Comm comm, istream &input, bool refine)
-   : glob_elem_offset(-1)
+   : face_nbr_el_to_face(NULL)
+   , glob_elem_offset(-1)
    , glob_offset_sequence(-1)
    , gtopo(comm)
 {
@@ -1129,7 +1134,7 @@ void ParMesh::MakeRefined_(ParMesh &orig_mesh, int ref_factor, int ref_type)
    MyComm = orig_mesh.GetComm();
    NRanks = orig_mesh.GetNRanks();
    MyRank = orig_mesh.GetMyRank();
-   face_nbr_el_to_face = nullptr;
+   face_nbr_el_to_face = NULL;
    glob_elem_offset = -1;
    glob_offset_sequence = -1;
    gtopo = orig_mesh.gtopo;
@@ -2118,7 +2123,7 @@ void ParMesh::ExchangeFaceNbrData()
 
    if (Nonconforming())
    {
-      // with ParNCMesh we can set up face neighbors mostly without communication
+      // with ParNCMesh we can set up face neighbors without communication
       pncmesh->GetFaceNeighbors(*this);
       have_face_nbr_data = true;
 
@@ -2187,7 +2192,7 @@ void ParMesh::ExchangeFaceNbrData()
 
    if (Dim == 3)
    {
-      BuildFaceNbrElementToFaceTable();
+      GetFaceNbrElementToFaceTable();
    }
 
    if (del_tables) { delete gr_sface; }
@@ -2452,7 +2457,8 @@ void ParMesh::ExchangeFaceNbrData(Table *gr_sface, int *s2l_face)
 
    // convert the element data into face_nbr_elements
    face_nbr_elements.SetSize(face_nbr_elements_offset[num_face_nbrs]);
-   face_nbr_el_ori.reset(new Table(face_nbr_elements_offset[num_face_nbrs], 6));
+   face_nbr_el_ori.Clear();
+   face_nbr_el_ori.SetSize(face_nbr_elements_offset[num_face_nbrs], 6);
    while (true)
    {
       int fn;
@@ -2483,7 +2489,7 @@ void ParMesh::ExchangeFaceNbrData(Table *gr_sface, int *s2l_face)
          if (Dim == 3)
          {
             int nf = el->GetNFaces();
-            int * fn_ori = face_nbr_el_ori->GetRow(elem_off);
+            int * fn_ori = face_nbr_el_ori.GetRow(elem_off);
             for (int j = 0; j < nf; j++)
             {
                fn_ori[j] = recv_elemdata[j];
@@ -2493,7 +2499,7 @@ void ParMesh::ExchangeFaceNbrData(Table *gr_sface, int *s2l_face)
          face_nbr_elements[elem_off++] = el;
       }
    }
-   face_nbr_el_ori->Finalize();
+   face_nbr_el_ori.Finalize();
 
    MPI_Waitall(num_face_nbrs, send_requests, statuses);
 
@@ -2715,118 +2721,191 @@ STable3D *ParMesh::GetSharedFacesTable()
    return sfaces_tbl;
 }
 
-template <int N>
-void
-ParMesh::AddTriFaces(const Array<int> &elem_vertices,
-                     const std::unique_ptr<STable3D> &faces,
-                     const std::unique_ptr<STable3D> &shared_faces,
-                     int elem, int start, int end, const int fverts[][N])
+STable3D *ParMesh::GetFaceNbrElementToFaceTable(int ret_ftbl)
 {
-   for (int i = start; i < end; ++i)
+   int i, *v;
+   STable3D * faces_tbl = GetFacesTable();
+   STable3D * sfaces_tbl = GetSharedFacesTable();
+
+   if (face_nbr_el_to_face != NULL)
    {
-      // Reference face vertices.
-      const auto fv = fverts[i];
-      // Element specific face vertices.
-      const Vert3 elem_fv(elem_vertices[fv[0]], elem_vertices[fv[1]],
-                          elem_vertices[fv[2]]);
-
-      // Check amongst the faces of elements local to this rank for this set of vertices
-      const int lf = faces->Index(elem_fv.v[0], elem_fv.v[1], elem_fv.v[2]);
-
-      // If the face wasn't found amonst processor local elements, search the
-      // ghosts for this set of vertices.
-      const int sf = lf < 0 ? shared_faces->Index(elem_fv.v[0], elem_fv.v[1],
-                                                  elem_fv.v[2]) : -1;
-      // If find local face -> use that
-      //    else if find shared face -> shift and use that
-      //       else no face found -> set to -1
-      const int face_to_add = lf < 0 ? (sf >= 0 ? sf + NumOfFaces : -1) : lf;
-
-      MFEM_ASSERT(sf >= 0 ||
-                  lf >= 0, "Face must be from a local or a face neighbor element");
-
-      // Add this discovered face to the list of faces of this face neighbor element
-      face_nbr_el_to_face->Push(elem, face_to_add);
+      delete face_nbr_el_to_face;
    }
-}
-
-void ParMesh::BuildFaceNbrElementToFaceTable()
-{
-   const auto faces = std::unique_ptr<STable3D>(GetFacesTable());
-   const auto shared_faces = std::unique_ptr<STable3D>(GetSharedFacesTable());
-
-   face_nbr_el_to_face.reset(new Table(face_nbr_elements.Size(), 6));
-
-   Array<int> v;
-
-   // Helper for adding quadrilateral faces.
-   auto add_quad_faces = [&faces, &shared_faces, &v, this]
-                         (int elem, int start, int end, const int fverts[][4])
+   face_nbr_el_to_face = new Table(face_nbr_elements.Size(), 6);
+   for (i = 0; i < face_nbr_elements.Size(); i++)
    {
-      for (int i = start; i < end; ++i)
-      {
-         const int * const fv = fverts[i];
-         int k = 0;
-         int max = v[fv[0]];
-
-         if (max < v[fv[1]]) { max = v[fv[1]], k = 1; }
-         if (max < v[fv[2]]) { max = v[fv[2]], k = 2; }
-         if (max < v[fv[3]]) { k = 3; }
-
-         int v0 = -1, v1 = -1, v2 = -1;
-         switch (k)
-         {
-            case 0:
-               v0 = v[fv[1]]; v1 = v[fv[2]]; v2 = v[fv[3]];
-               break;
-            case 1:
-               v0 = v[fv[0]]; v1 = v[fv[2]]; v2 = v[fv[3]];
-               break;
-            case 2:
-               v0 = v[fv[0]]; v1 = v[fv[1]]; v2 = v[fv[3]];
-               break;
-            case 3:
-               v0 = v[fv[0]]; v1 = v[fv[1]]; v2 = v[fv[2]];
-               break;
-         }
-         int lf = faces->Index(v0, v1, v2);
-         if (lf < 0)
-         {
-            lf = shared_faces->Index(v0, v1, v2);
-            if (lf >= 0)
-            {
-               lf += NumOfFaces;
-            }
-         }
-         face_nbr_el_to_face->Push(elem, lf);
-      }
-   };
-
-   for (int i = 0; i < face_nbr_elements.Size(); i++)
-   {
-      face_nbr_elements[i]->GetVertices(v);
+      v = face_nbr_elements[i]->GetVertices();
       switch (face_nbr_elements[i]->GetType())
       {
          case Element::TETRAHEDRON:
          {
-            AddTriFaces(v, faces, shared_faces, i, 0, 4, tet_t::FaceVert);
+            for (int j = 0; j < 4; j++)
+            {
+               const int *fv = tet_t::FaceVert[j];
+               int lf = faces_tbl->Index(v[fv[0]], v[fv[1]], v[fv[2]]);
+               if (lf < 0)
+               {
+                  lf = sfaces_tbl->Index(v[fv[0]], v[fv[1]], v[fv[2]]);
+                  if (lf >= 0)
+                  {
+                     lf += NumOfFaces;
+                  }
+               }
+               face_nbr_el_to_face->Push(i, lf);
+            }
             break;
          }
          case Element::WEDGE:
          {
-            AddTriFaces(v, faces, shared_faces, i, 0, 2, pri_t::FaceVert);
-            add_quad_faces(i, 2, 5, pri_t::FaceVert);
+            for (int j = 0; j < 2; j++)
+            {
+               const int *fv = pri_t::FaceVert[j];
+               int lf = faces_tbl->Index(v[fv[0]], v[fv[1]], v[fv[2]]);
+               if (lf < 0)
+               {
+                  lf = sfaces_tbl->Index(v[fv[0]], v[fv[1]], v[fv[2]]);
+                  if (lf >= 0)
+                  {
+                     lf += NumOfFaces;
+                  }
+               }
+               face_nbr_el_to_face->Push(i, lf);
+            }
+            for (int j = 2; j < 5; j++)
+            {
+               const int *fv = pri_t::FaceVert[j];
+               int k = 0;
+               int max = v[fv[0]];
+
+               if (max < v[fv[1]]) { max = v[fv[1]], k = 1; }
+               if (max < v[fv[2]]) { max = v[fv[2]], k = 2; }
+               if (max < v[fv[3]]) { k = 3; }
+
+               int v0 = -1, v1 = -1, v2 = -1;
+               switch (k)
+               {
+                  case 0:
+                     v0 = v[fv[1]]; v1 = v[fv[2]]; v2 = v[fv[3]];
+                     break;
+                  case 1:
+                     v0 = v[fv[0]]; v1 = v[fv[2]]; v2 = v[fv[3]];
+                     break;
+                  case 2:
+                     v0 = v[fv[0]]; v1 = v[fv[1]]; v2 = v[fv[3]];
+                     break;
+                  case 3:
+                     v0 = v[fv[0]]; v1 = v[fv[1]]; v2 = v[fv[2]];
+                     break;
+               }
+               int lf = faces_tbl->Index(v0, v1, v2);
+               if (lf < 0)
+               {
+                  lf = sfaces_tbl->Index(v0, v1, v2);
+                  if (lf >= 0)
+                  {
+                     lf += NumOfFaces;
+                  }
+               }
+               face_nbr_el_to_face->Push(i, lf);
+            }
             break;
          }
          case Element::PYRAMID:
          {
-            add_quad_faces(i, 0, 1, pyr_t::FaceVert);
-            AddTriFaces(v, faces, shared_faces, i, 1, 5, pyr_t::FaceVert);
+            for (int j = 0; j < 1; j++)
+            {
+               const int *fv = pyr_t::FaceVert[j];
+               int k = 0;
+               int max = v[fv[0]];
+
+               if (max < v[fv[1]]) { max = v[fv[1]], k = 1; }
+               if (max < v[fv[2]]) { max = v[fv[2]], k = 2; }
+               if (max < v[fv[3]]) { k = 3; }
+
+               int v0 = -1, v1 = -1, v2 = -1;
+               switch (k)
+               {
+                  case 0:
+                     v0 = v[fv[1]]; v1 = v[fv[2]]; v2 = v[fv[3]];
+                     break;
+                  case 1:
+                     v0 = v[fv[0]]; v1 = v[fv[2]]; v2 = v[fv[3]];
+                     break;
+                  case 2:
+                     v0 = v[fv[0]]; v1 = v[fv[1]]; v2 = v[fv[3]];
+                     break;
+                  case 3:
+                     v0 = v[fv[0]]; v1 = v[fv[1]]; v2 = v[fv[2]];
+                     break;
+               }
+               int lf = faces_tbl->Index(v0, v1, v2);
+               if (lf < 0)
+               {
+                  lf = sfaces_tbl->Index(v0, v1, v2);
+                  if (lf >= 0)
+                  {
+                     lf += NumOfFaces;
+                  }
+               }
+               face_nbr_el_to_face->Push(i, lf);
+            }
+            for (int j = 1; j < 5; j++)
+            {
+               const int *fv = pyr_t::FaceVert[j];
+               int lf = faces_tbl->Index(v[fv[0]], v[fv[1]], v[fv[2]]);
+               if (lf < 0)
+               {
+                  lf = sfaces_tbl->Index(v[fv[0]], v[fv[1]], v[fv[2]]);
+                  if (lf >= 0)
+                  {
+                     lf += NumOfFaces;
+                  }
+               }
+               face_nbr_el_to_face->Push(i, lf);
+            }
             break;
          }
          case Element::HEXAHEDRON:
          {
-            add_quad_faces(i, 0, 6, hex_t::FaceVert);
+            // find the face by the vertices with the smallest 3 numbers
+            // z = 0, y = 0, x = 1, y = 1, x = 0, z = 1
+            for (int j = 0; j < 6; j++)
+            {
+               const int *fv = hex_t::FaceVert[j];
+               int k = 0;
+               int max = v[fv[0]];
+
+               if (max < v[fv[1]]) { max = v[fv[1]], k = 1; }
+               if (max < v[fv[2]]) { max = v[fv[2]], k = 2; }
+               if (max < v[fv[3]]) { k = 3; }
+
+               int v0 = -1, v1 = -1, v2 = -1;
+               switch (k)
+               {
+                  case 0:
+                     v0 = v[fv[1]]; v1 = v[fv[2]]; v2 = v[fv[3]];
+                     break;
+                  case 1:
+                     v0 = v[fv[0]]; v1 = v[fv[2]]; v2 = v[fv[3]];
+                     break;
+                  case 2:
+                     v0 = v[fv[0]]; v1 = v[fv[1]]; v2 = v[fv[3]];
+                     break;
+                  case 3:
+                     v0 = v[fv[0]]; v1 = v[fv[1]]; v2 = v[fv[2]];
+                     break;
+               }
+               int lf = faces_tbl->Index(v0, v1, v2);
+               if (lf < 0)
+               {
+                  lf = sfaces_tbl->Index(v0, v1, v2);
+                  if (lf >= 0)
+                  {
+                     lf += NumOfFaces;
+                  }
+               }
+               face_nbr_el_to_face->Push(i, lf);
+            }
             break;
          }
          default:
@@ -2834,6 +2913,14 @@ void ParMesh::BuildFaceNbrElementToFaceTable()
       }
    }
    face_nbr_el_to_face->Finalize();
+
+   delete sfaces_tbl;
+   if (ret_ftbl)
+   {
+      return faces_tbl;
+   }
+   delete faces_tbl;
+   return NULL;
 }
 
 int ParMesh::GetFaceNbrRank(int fn) const
@@ -2854,28 +2941,33 @@ int ParMesh::GetFaceNbrRank(int fn) const
 }
 
 void
-ParMesh::GetFaceNbrElementFaces(int i, Array<int> &faces,
-                                Array<int> &orientations) const
+ParMesh::GetFaceNbrElementFaces(int i, Array<int> &fcs, Array<int> &cor) const
 {
+   int n, j;
    int el_nbr = i - GetNE();
-   if (face_nbr_el_to_face != nullptr && el_nbr < face_nbr_el_to_face->Size())
+   if (face_nbr_el_to_face)
    {
-      face_nbr_el_to_face->GetRow(el_nbr, faces);
+      face_nbr_el_to_face->GetRow(el_nbr, fcs);
    }
    else
    {
       MFEM_ABORT("ParMesh::GetFaceNbrElementFaces(...) : "
-                 "face_nbr_el_to_face not generated correctly.");
+                 "face_nbr_el_to_face not generated.");
    }
-
-   if (face_nbr_el_ori != nullptr && el_nbr < face_nbr_el_ori->Size())
+   if (el_nbr < face_nbr_el_ori.Size())
    {
-      face_nbr_el_ori->GetRow(el_nbr, orientations);
+      const int * row = face_nbr_el_ori.GetRow(el_nbr);
+      n = fcs.Size();
+      cor.SetSize(n);
+      for (j=0; j<n; j++)
+      {
+         cor[j] = row[j];
+      }
    }
    else
    {
       MFEM_ABORT("ParMesh::GetFaceNbrElementFaces(...) : "
-                 "face_nbr_el_ori not generated correctly.");
+                 "face_nbr_el_to_face not generated.");
    }
 }
 
@@ -3135,7 +3227,7 @@ int ParMesh::GetSharedFace(int sface) const
    {
       MFEM_ASSERT(Dim > 1, "");
       const NCMesh::NCList &shared = pncmesh->GetSharedList(Dim-1);
-      int csize = shared.conforming.Size();
+      int csize = (int) shared.conforming.Size();
       return sface < csize
              ? shared.conforming[sface].index
              : shared.slaves[sface - csize].index;
@@ -4487,7 +4579,7 @@ void ParMesh::UniformRefinement3D()
 
    DSTable v_to_v(NumOfVertices);
    GetVertexToVertexTable(v_to_v);
-   auto faces_tbl = std::unique_ptr<STable3D>(GetFacesTable());
+   STable3D *faces_tbl = GetFacesTable();
 
    // call Mesh::UniformRefinement3D_base so that it won't update the nodes
    Array<int> f2qf;
@@ -4502,6 +4594,7 @@ void ParMesh::UniformRefinement3D()
    // update the groups
    UniformRefineGroups3D(old_nv, old_nedges, v_to_v, *faces_tbl,
                          f2qf.Size() ? &f2qf : NULL);
+   delete faces_tbl;
 
    UpdateNodes();
 }
@@ -6600,8 +6693,6 @@ void ParMesh::Swap(ParMesh &other)
    mfem::Swap(face_nbr_vertices, other.face_nbr_vertices);
    mfem::Swap(send_face_nbr_elements, other.send_face_nbr_elements);
    mfem::Swap(send_face_nbr_vertices, other.send_face_nbr_vertices);
-   std::swap(face_nbr_el_ori, other.face_nbr_el_ori);
-   std::swap(face_nbr_el_to_face, other.face_nbr_el_to_face);
 
    // Nodes, NCMesh, and NURBSExtension are taken care of by Mesh::Swap
    mfem::Swap(pncmesh, other.pncmesh);
@@ -6622,7 +6713,8 @@ void ParMesh::Destroy()
    }
    shared_edges.DeleteAll();
 
-   face_nbr_el_to_face = nullptr;
+   delete face_nbr_el_to_face;
+   face_nbr_el_to_face = NULL;
 }
 
 ParMesh::~ParMesh()

mesh/pmesh.hpp

@@ -77,11 +77,8 @@ protected:
    // sface ids: all triangles first, then all quads
    Array<int> sface_lface;
 
-   /// Table that maps from face neighbor element number, to the face numbers of
-   /// that element.
-   std::unique_ptr<Table> face_nbr_el_to_face;
-   /// orientations for each face (from nbr processor)
-   std::unique_ptr<Table> face_nbr_el_ori;
+   Table *face_nbr_el_to_face;
+   Table  face_nbr_el_ori; // orientations for each face (from nbr processor)
 
    IsoparametricTransformation FaceNbrTransformation;
 
@@ -116,32 +113,7 @@ protected:
    bool DecodeFaceSplittings(HashTable<Hashed2> &v_to_v, const int *v,
                              const Array<unsigned> &codes, int &pos);
 
-   // Given a completed FacesTable and SharedFacesTable, construct a table that
-   // maps from face neighbor element number, to the set of faces of that
-   // element. Store the resulting data in the member variable
-   // face_nbr_el_to_face. If the mesh is nonconforming, this also builds the
-   // the face_nbr_el_ori variable from the faces_info.
-   void BuildFaceNbrElementToFaceTable();
-
-   /**
-    * @brief Helper function for adding triangle face neighbor element to face
-    * table entries. Have to use a template here rather than lambda capture
-    * because the FaceVert entries in Geometry have inner size of 3 for tets and
-    * 4 for everything else.
-    *
-    * @tparam N Inner dimension on the fvert variable, 3 for tet, 4 otherwise
-    * @param[in] v Set of vertices for this element
-    * @param[in] faces Table of faces interior to this rank
-    * @param[in] shared_faces Table of faces shared by this rank and another
-    * @param[in] elem The face neighbor element
-    * @param[in] start Starting index into fverts
-    * @param[in] end End index into fverts
-    * @param[in] fverts Array of face vertices for this particular geometry.
-    */
-   template <int N>
-   void AddTriFaces(const Array<int> &v, const std::unique_ptr<STable3D> &faces,
-                    const std::unique_ptr<STable3D> &shared_faces,
-                    int elem, int start, int end, const int fverts[][N]);
+   STable3D *GetFaceNbrElementToFaceTable(int ret_ftbl = 0);
 
    void GetFaceNbrElementTransformation(
       int i, IsoparametricTransformation *ElTr);
@@ -315,7 +287,7 @@ protected:
 
 public:
    /// Default constructor. Create an empty @a ParMesh.
-   ParMesh() : MyComm(0), NRanks(0), MyRank(-1),
+   ParMesh() : MyComm(0), NRanks(0), MyRank(-1), face_nbr_el_to_face(NULL),
       glob_elem_offset(-1), glob_offset_sequence(-1),
       have_face_nbr_data(false), pncmesh(NULL) { }
 
@@ -493,8 +465,7 @@ public:
    int GetFaceNbrRank(int fn) const;
 
    /** Similar to Mesh::GetElementFaces */
-   void GetFaceNbrElementFaces(int i, Array<int> &faces,
-                               Array<int> &orientation) const;
+   void GetFaceNbrElementFaces(int i, Array<int> &fcs, Array<int> &cor) const;
 
    /** Similar to Mesh::GetFaceToElementTable with added face-neighbor elements
        with indices offset by the local number of elements. */

mesh/pncmesh.cpp

@@ -16,12 +16,10 @@
 #include "mesh_headers.hpp"
 #include "pncmesh.hpp"
 #include "../general/binaryio.hpp"
-#include "../general/communication.hpp"
 
 #include <numeric> // std::accumulate
 #include <map>
 #include <climits> // INT_MIN, INT_MAX
-#include <array>
 
 namespace mfem
 {
@@ -889,7 +887,6 @@ void ParNCMesh::GetFaceNeighbors(ParMesh &pmesh)
 
    Array<Element*> fnbr;
    Array<Connection> send_elems;
-   std::map<int, std::vector<int>> recv_elems;
 
    // Counts the number of slave faces of a master. This may be larger than the
    // number of shared slaves if there exist degenerate slave-faces from face-edge constraints.
@@ -904,11 +901,6 @@ void ParNCMesh::GetFaceNeighbors(ParMesh &pmesh)
    fnbr.Reserve(bound);
    send_elems.Reserve(bound);
 
-   // If there are face neighbor elements with triangular faces, the
-   // `face_nbr_el_ori` structure will need to be built. This requires
-   // communication so we attempt to avoid it by checking first.
-   bool face_nbr_w_tri_faces = false;
-
    // go over all shared faces and collect face neighbor elements
    for (int i = 0; i < shared.conforming.Size(); i++)
    {
@@ -922,12 +914,8 @@ void ParNCMesh::GetFaceNeighbors(ParMesh &pmesh)
       if (e[0]->rank == MyRank) { std::swap(e[0], e[1]); }
       MFEM_ASSERT(e[0]->rank != MyRank && e[1]->rank == MyRank, "");
 
-      face_nbr_w_tri_faces |= !Geometry::IsTensorProduct(Geometry::Type(e[0]->geom));
-      face_nbr_w_tri_faces |= !Geometry::IsTensorProduct(Geometry::Type(e[1]->geom));
-
       fnbr.Append(e[0]);
       send_elems.Append(Connection(e[0]->rank, e[1]->index));
-      recv_elems[e[0]->rank].push_back(e[0]->index);
    }
 
    for (int i = 0; i < shared.masters.Size(); i++)
@@ -936,7 +924,7 @@ void ParNCMesh::GetFaceNeighbors(ParMesh &pmesh)
       for (int j = mf.slaves_begin; j < mf.slaves_end; j++)
       {
          const Slave &sf = full_list.slaves[j];
-         if (sf.element < 0 || sf.index < 0) { continue; }
+         if (sf.element < 0) { continue; }
 
          MFEM_ASSERT(mf.element >= 0, "");
          Element* e[2] = { &elements[mf.element], &elements[sf.element] };
@@ -950,12 +938,8 @@ void ParNCMesh::GetFaceNeighbors(ParMesh &pmesh)
          }
          if (loc0) { std::swap(e[0], e[1]); }
 
-         face_nbr_w_tri_faces |= !Geometry::IsTensorProduct(Geometry::Type(e[0]->geom));
-         face_nbr_w_tri_faces |= !Geometry::IsTensorProduct(Geometry::Type(e[1]->geom));
-
          fnbr.Append(e[0]);
          send_elems.Append(Connection(e[0]->rank, e[1]->index));
-         recv_elems[e[0]->rank].push_back(e[0]->index);
       }
    }
 
@@ -1038,13 +1022,6 @@ void ParNCMesh::GetFaceNeighbors(ParMesh &pmesh)
    send_elems.Sort();
    send_elems.Unique();
 
-   for (auto &kv : recv_elems)
-   {
-      std::sort(kv.second.begin(), kv.second.end());
-      kv.second.erase(std::unique(kv.second.begin(), kv.second.end()),
-                      kv.second.end());
-   }
-
    for (int i = 0, last_rank = -1; i < send_elems.Size(); i++)
    {
       Connection &c = send_elems[i];
@@ -1198,117 +1175,6 @@ void ParNCMesh::GetFaceNeighbors(ParMesh &pmesh)
       }
    }
 
-
-   // In 3D some extra orientation data structures can be needed.
-   if (Dim == 3)
-   {
-      // Populates face_nbr_el_to_face, always needed.
-      pmesh.BuildFaceNbrElementToFaceTable();
-
-      if (face_nbr_w_tri_faces)
-      {
-         // There are face neighbor elements with triangular faces, need to
-         // perform communication to ensure the orientation is valid.
-         using RankToOrientation = std::map<int, std::vector<std::array<int, 6>>>;
-         constexpr std::array<int, 6> unset_ori{{-1,-1,-1,-1,-1,-1}};
-         const int rank = pmesh.GetMyRank();
-
-         // Loop over send elems, compute the orientation and place in the
-         // buffer to send to each processor. Note elements are
-         // lexicographically sorted with rank and element number, and this
-         // ordering holds across processors.
-         RankToOrientation send_rank_to_face_neighbor_orientations;
-         Array<int> orientations, faces;
-
-         // send_elems goes from rank of the receiving processor, to the index
-         // of the face neighbor element on this processor.
-         for (const auto &se : send_elems)
-         {
-            const auto &true_rank = pmesh.face_nbr_group[se.from];
-            pmesh.GetElementFaces(se.to, faces, orientations);
-
-            // Place a new entry of unset orientations
-            send_rank_to_face_neighbor_orientations[true_rank].emplace_back(unset_ori);
-
-            // Copy the entries, any unset faces will remain -1.
-            std::copy(orientations.begin(), orientations.end(),
-                      send_rank_to_face_neighbor_orientations[true_rank].back().begin());
-         }
-
-         // Initialize the receive buffers and resize to match the expected
-         // number of elements coming in. The copy ensures the appropriate rank
-         // pairings are in place, and for a purely conformal interface, the
-         // resize is a no-op.
-         auto recv_rank_to_face_neighbor_orientations =
-            send_rank_to_face_neighbor_orientations;
-         for (auto &kv : recv_rank_to_face_neighbor_orientations)
-         {
-            kv.second.resize(recv_elems[kv.first].size());
-         }
-
-         // For asynchronous send/recv, will use arrays of requests to monitor the
-         // status of the connections.
-         std::vector<MPI_Request> send_requests, recv_requests;
-         std::vector<MPI_Status> status(nranks);
-
-         // NOTE: This is CRITICAL, to ensure the addresses of these requests
-         // do not change between the send/recv and the wait.
-         send_requests.reserve(nranks);
-         recv_requests.reserve(nranks);
-
-         // Shared face communication is bidirectional -> any rank to whom
-         // orientations must be sent, will need to send orientations back. The
-         // orientation data is contiguous because std::array<int,6> is an
-         // aggregate. Loop over each communication pairing, and dispatch the
-         // buffer loaded with  all the orientation data.
-         for (const auto &kv : send_rank_to_face_neighbor_orientations)
-         {
-            send_requests.emplace_back(); // instantiate a request for tracking.
-
-            // low rank sends on low, high rank sends on high.
-            const int send_tag = (rank < kv.first)
-                                 ? std::min(rank, kv.first)
-                                 : std::max(rank, kv.first);
-            MPI_Isend(&kv.second[0][0], int(kv.second.size() * 6),
-                      MPI_INT, kv.first, send_tag, pmesh.MyComm, &send_requests.back());
-         }
-
-         // Loop over the communication pairing again, and receive the
-         // symmetric buffer from the other processor.
-         for (auto &kv : recv_rank_to_face_neighbor_orientations)
-         {
-            recv_requests.emplace_back(); // instantiate a request for tracking
-
-            // low rank receives on high, high rank receives on low.
-            const int recv_tag = (rank < kv.first)
-                                 ? std::max(rank, kv.first)
-                                 : std::min(rank, kv.first);
-            MPI_Irecv(&kv.second[0][0], int(kv.second.size() * 6),
-                      MPI_INT, kv.first, recv_tag, pmesh.MyComm, &recv_requests.back());
-         }
-
-         // Wait until all receive buffers are full before beginning to process.
-         MPI_Waitall(int(recv_requests.size()), recv_requests.data(), status.data());
-
-         pmesh.face_nbr_el_ori.reset(new Table(pmesh.face_nbr_elements.Size(), 6));
-         int elem = 0;
-         for (const auto &kv : recv_rank_to_face_neighbor_orientations)
-         {
-            // All elements associated to this face-neighbor rank
-            for (const auto &eo : kv.second)
-            {
-               std::copy(eo.begin(), eo.end(), pmesh.face_nbr_el_ori->GetRow(elem));
-               ++elem;
-            }
-         }
-         pmesh.face_nbr_el_ori->Finalize();
-
-         // Must wait for all send buffers to be released before the scope closes.
-         MPI_Waitall(int(send_requests.size()), send_requests.data(), status.data());
-      }
-   }
-
-
    // NOTE: this function skips ParMesh::send_face_nbr_vertices and
    // ParMesh::face_nbr_vertices_offset, these are not used outside of ParMesh
 }

mesh/pncmesh.hpp

@@ -108,7 +108,9 @@ public:
        passed. */
    void Rebalance(const Array<int> *custom_partition = NULL);
 
+
    // interface for ParFiniteElementSpace
+
    int GetNElements() const { return NElements; }
 
    int GetNGhostVertices() const { return NGhostVertices; }
@@ -139,8 +141,8 @@ public:
       return (index < NFaces) ? face_orient[index] : 0;
    }
 
-   using GroupId = short;
-   using CommGroup = std::vector<int>;
+   typedef short GroupId;
+   typedef std::vector<int> CommGroup;
 
    /// Return vertex/edge/face ('entity' == 0/1/2, resp.) owner.
    GroupId GetEntityOwnerId(int entity, int index)

miniapps/gslib/findpts.cpp

@@ -347,7 +347,8 @@ int main (int argc, char *argv[])
       {
          if ((double) rand() / RAND_MAX < 0.5)
          {
-            sc_fes.SetElementOrder(e, order + 1);
+            int element_order = sc_fes.GetElementOrder(e);
+            sc_fes.SetElementOrder(e, element_order + 1);
          }
       }
       sc_fes.Update(false);

miniapps/hdiv-linear-solver/discrete_divergence.cpp

@@ -85,7 +85,8 @@ void EliminateColumns(HypreParMatrix &D, const Array<int> &ess_dofs)
       const auto I = diag->i;
       const auto J = diag->j;
       auto data = diag->data;
-      MFEM_HYPRE_FORALL(i, nrows_diag,
+      //MFEM_HYPRE_FORALL(i, nrows_diag,
+      mfem::forall_switch(HypreUsingGPU(), nrows_diag, [=] MFEM_HOST_DEVICE (int i)
       {
          for (int jj=I[i]; jj<I[i+1]; ++jj)
          {
@@ -109,7 +110,8 @@ void EliminateColumns(HypreParMatrix &D, const Array<int> &ess_dofs)
       const auto I = offd->i;
       const auto J = offd->j;
       auto data = offd->data;
-      MFEM_HYPRE_FORALL(i, nrows_offd,
+      //MFEM_HYPRE_FORALL(i, nrows_offd,
+      mfem::forall_switch(HypreUsingGPU(), nrows_offd, [=] MFEM_HOST_DEVICE (int i)
       {
          for (int jj=I[i]; jj<I[i+1]; ++jj)
          {

miniapps/nurbs/CMakeLists.txt

@@ -17,10 +17,6 @@ add_mfem_miniapp(nurbs_curveint
   MAIN nurbs_curveint.cpp
   LIBRARIES mfem)
 
-add_mfem_miniapp(nurbs_printfunc
-  MAIN nurbs_printfunc.cpp
-  LIBRARIES mfem)
-
 add_mfem_miniapp(nurbs_patch_ex1
   MAIN nurbs_patch_ex1.cpp
   LIBRARIES mfem)
@@ -71,9 +67,6 @@ if (MFEM_ENABLE_TESTING)
   add_test(NAME nurbs_ex1_weak_patch_format_r1_ser
     COMMAND $<TARGET_FILE:nurbs_ex1> -no-vis
     -m ${PROJECT_SOURCE_DIR}/data/square-disc-nurbs-patch.mesh -o 2 --weak-bc -r 1)
-    
-  add_test(NAME nurbs_printfunc
-    COMMAND $<TARGET_FILE:nurbs_printfunc>)
 
   if (MFEM_USE_LAPACK)
     add_test(NAME nurbs_patch_ex1_o4_r2_iro8_patcha_ser

miniapps/nurbs/makefile

@@ -21,7 +21,7 @@ CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
 MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
-SEQ_MINIAPPS = nurbs_ex1 nurbs_patch_ex1 nurbs_curveint nurbs_printfunc
+SEQ_MINIAPPS = nurbs_ex1 nurbs_patch_ex1 nurbs_curveint
 PAR_MINIAPPS = nurbs_ex1p nurbs_ex11p
 ifeq ($(MFEM_USE_MPI),NO)
    MINIAPPS = $(SEQ_MINIAPPS)
@@ -109,9 +109,6 @@ nurbs_curveint-test-seq: nurbs_curveint
 	@$(call mfem-test,$<,, NURBS miniapp,$(CI_ARGS_1))
 	@$(call mfem-test,$<,, NURBS miniapp,$(CI_ARGS_2))
 
-nurbs_printfunc-test-seq: nurbs_printfunc
-	@$(call mfem-test,$<,, NURBS miniapp)
-
 EX1P_ARGS_1 :=
 EX1P_ARGS_2 := -m ../../data/pipe-nurbs-2d.mesh -o 2 -no-ibp
 EX1P_ARGS_3 := -m ../../data/ball-nurbs.mesh -o 2 --weak-bc -r 0

miniapps/nurbs/nurbs_printfunc.cpp

@@ -1,52 +0,0 @@
-//                          MFEM NURBS knot vector example
-//
-// Compile with: make nurbs_curveint
-//
-// Sample runs:  nurbs_curveint
-//
-// Description:  This example code demonstrates the use of MFEM to define a
-//               simple KnotVector and print its corresponding shape functions.
-
-#include <iostream>
-#include "mfem.hpp"
-
-using namespace std;
-using namespace mfem;
-
-int main(int argc, char *argv[])
-{
-   OptionsParser args(argc, argv);
-   bool visualization;
-
-   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
-                  "--no-visualization",
-                  "Enable or disable GLVis visualization. Dummy option to allow testing.");
-   args.Parse();
-   if (!args.Good())
-   {
-      args.PrintUsage(cout);
-      return 1;
-   }
-
-   KnotVector kv(2, 7);
-
-   kv[0] = 0;
-   kv[1] = 0;
-   kv[2] = 0;
-   kv[3] = 0.25;
-   kv[4] = 0.5;
-   kv[5] = 0.5; // Repeated knot
-   kv[6] = 0.75;
-   kv[7] = 1;
-   kv[8] = 1;
-   kv[9] = 1;
-
-   cout << "Printing knotvector:" << endl;
-   kv.Print(cout);
-
-   // Count number of elements, required for printing of shapes
-   kv.GetElements();
-
-   cout << "\nPrinting shapefunctions:" << endl;
-   kv.PrintFunctions(cout);
-}

tests/unit/fem/test_gslib.cpp

@@ -1,372 +0,0 @@
-// Copyright (c) 2010-2023, 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 "unit_tests.hpp"
-#include "mfem.hpp"
-
-using namespace mfem;
-#ifdef MFEM_USE_GSLIB
-namespace gslib_test
-{
-
-int func_order;
-
-// Scalar function to project
-double scalar_func(const Vector &x)
-{
-   const int dim = x.Size();
-   double res = 0.0;
-   for (int d = 0; d < dim; d++) { res += std::pow(x(d), func_order); }
-   return res;
-}
-
-void F_exact(const Vector &p, Vector &F)
-{
-   F(0) = scalar_func(p);
-   for (int i = 1; i < F.Size(); i++) { F(i) = (i+1)*F(0); }
-}
-
-enum class Space { H1, L2 };
-
-TEST_CASE("GSLIBInterpolate", "[GSLIBInterpolate][GSLIB]")
-{
-   auto space               = GENERATE(Space::H1, Space::L2);
-   auto simplex             = GENERATE(true, false);
-   int dim                  = GENERATE(2, 3);
-   func_order               = GENERATE(1, 2);
-   int mesh_order           = GENERATE(1, 2);
-   int mesh_node_ordering   = GENERATE(0, 1);
-   int point_ordering       = GENERATE(0, 1);
-   int ncomp                = GENERATE(1, 2);
-   int gf_ordering          = GENERATE(0, 1);
-   bool href                = GENERATE(true, false);
-   bool pref                = GENERATE(true, false);
-
-   int ne = 4;
-
-   CAPTURE(space, simplex, dim, func_order, mesh_order, mesh_node_ordering,
-           point_ordering, ncomp, gf_ordering, href, pref);
-
-   if (ncomp == 1 && gf_ordering == 1)
-   {
-      return;
-   }
-
-   Mesh mesh;
-   if (dim == 2)
-   {
-      Element::Type type = simplex ? Element::TRIANGLE : Element::QUADRILATERAL;
-      mesh = Mesh::MakeCartesian2D(ne, ne, type, 1, 1.0, 1.0);
-   }
-   else
-   {
-      Element::Type type = simplex ? Element::TETRAHEDRON : Element::HEXAHEDRON;
-      mesh = Mesh::MakeCartesian3D(ne, ne, ne, type, 1.0, 1.0, 1.0);
-   }
-
-   if (href || pref) { mesh.EnsureNCMesh(); }
-   if (href) { mesh.RandomRefinement(0.5); }
-
-   // Set Mesh NodalFESpace
-   H1_FECollection fecm(mesh_order, dim);
-   FiniteElementSpace fespacem(&mesh, &fecm, dim, mesh_node_ordering);
-   mesh.SetNodalFESpace(&fespacem);
-
-   // Set GridFunction to be interpolated
-   FiniteElementCollection *c_fec = nullptr;
-
-   switch (space)
-   {
-      case Space::H1:
-         c_fec = new H1_FECollection(func_order, dim);
-         break;
-      case Space::L2:
-         c_fec = new L2_FECollection(func_order, dim);
-         break;
-   }
-
-   FiniteElementSpace c_fespace =
-      FiniteElementSpace(&mesh, c_fec, ncomp, gf_ordering);
-   GridFunction field_vals(&c_fespace);
-
-   VectorFunctionCoefficient F(ncomp, F_exact);
-   field_vals.ProjectCoefficient(F);
-
-   // Generate points in the domain
-   Vector pos_min, pos_max;
-   mesh.GetBoundingBox(pos_min, pos_max, mesh_order);
-   const int pts_cnt_1D = 5;
-   int pts_cnt = pow(pts_cnt_1D, dim);
-   Vector vxyz(pts_cnt * dim);
-   NodalTensorFiniteElement *el = NULL;
-   if (dim == 2)
-   {
-      el = new L2_QuadrilateralElement(pts_cnt_1D-1,BasisType::ClosedUniform);
-   }
-   else
-   {
-      el = new L2_HexahedronElement(pts_cnt_1D - 1, BasisType::ClosedUniform);
-   }
-   const IntegrationRule &ir = el->GetNodes();
-   for (int i = 0; i < ir.GetNPoints(); i++)
-   {
-      const IntegrationPoint &ip = ir.IntPoint(i);
-      if (point_ordering == Ordering::byNODES)
-      {
-         vxyz(i)           = pos_min(0) + ip.x * (pos_max(0)-pos_min(0));
-         vxyz(pts_cnt + i) = pos_min(1) + ip.y * (pos_max(1)-pos_min(1));
-         if (dim == 3)
-         {
-            vxyz(2*pts_cnt + i) = pos_min(2) + ip.z * (pos_max(2)-pos_min(2));
-         }
-      }
-      else
-      {
-         vxyz(i*dim + 0) = pos_min(0) + ip.x * (pos_max(0)-pos_min(0));
-         vxyz(i*dim + 1) = pos_min(1) + ip.y * (pos_max(1)-pos_min(1));
-         if (dim == 3)
-         {
-            vxyz(i*dim + 2) = pos_min(2) + ip.z * (pos_max(2)-pos_min(2));
-         }
-      }
-   }
-   delete el;
-
-   // Find and interpolate FE Function values
-   Vector interp_vals(pts_cnt*ncomp);
-   FindPointsGSLIB finder;
-   finder.Setup(mesh);
-   finder.SetL2AvgType(FindPointsGSLIB::NONE);
-   finder.Interpolate(vxyz, field_vals, interp_vals, point_ordering);
-   Array<unsigned int> code_out    = finder.GetCode();
-   Vector dist_p_out = finder.GetDist();
-
-   int not_found = 0;
-   double err = 0.0, max_err = 0.0, max_dist = 0.0;
-   Vector pos(dim);
-
-   for (int i = 0; i < pts_cnt; i++)
-   {
-      max_dist = std::max(max_dist, dist_p_out(i));
-      for (int d = 0; d < dim; d++)
-      {
-         pos(d) = point_ordering == Ordering::byNODES ?
-                  vxyz(d*pts_cnt + i) :
-                  vxyz(i*dim + d);
-      }
-      Vector exact_val(ncomp);
-      F_exact(pos, exact_val);
-      for (int j = 0; j < ncomp; j++)
-      {
-         if (code_out[i] < 2)
-         {
-            err = gf_ordering == Ordering::byNODES ?
-                  fabs(exact_val(j) - interp_vals[i + j*pts_cnt]) :
-                  fabs(exact_val(j) - interp_vals[i*ncomp + j]);
-            max_err  = std::max(max_err, err);
-         }
-         else
-         {
-            if (j == 0) { not_found++; }
-         }
-      }
-   }
-
-   REQUIRE(max_err < 1e-12);
-   REQUIRE(max_dist < 1e-10);
-   REQUIRE(not_found == 0);
-
-   finder.FreeData();
-   delete c_fec;
-}
-
-// Generates meshes with different element types, followed by points at
-// element faces and interior, and finally checks to see if these points are
-// correctly detected at element boundary or not.
-TEST_CASE("GSLIBFindAtElementBoundary",
-          "[GSLIBFindAtElementBoundary][GSLIB]")
-{
-   int dim  = GENERATE(2, 3);
-   CAPTURE(dim);
-   int nex = 4;
-   int mesh_order = 4;
-   int l2_order = 4;
-
-   int netype = dim == 2 ? 2 : 4; // 2 element types in 2D, 4 in 3D.
-   int estart = dim == 2 ? 2 : 4; // starts at index 2 in 2D, 4 in 3D
-
-   for (int et = estart; et < estart+netype; et++)
-   {
-      // H1 - order 1, L2 - order 0 for pyramids
-      if (et == 7)
-      {
-         mesh_order = 1;
-         l2_order = 0;
-      }
-      Mesh mesh;
-      if (dim == 2)
-      {
-         mesh = Mesh::MakeCartesian2D(nex, nex, (Element::Type)et);
-      }
-      else
-      {
-         mesh = Mesh::MakeCartesian3D(nex, nex, nex, (Element::Type)et);
-      }
-
-      mesh.SetCurvature(mesh_order);
-      const FiniteElementSpace *n_fespace = mesh.GetNodalFESpace();
-      const GridFunction *nodes = mesh.GetNodes();
-
-      Array<double> xyz;
-
-      // Generate points on each element's face/edge
-      for (int e = 0; e < mesh.GetNE(); e++)
-      {
-         Array<int> faces,ori;
-         if (dim == 2)
-         {
-            mesh.GetElementEdges(e, faces, ori);
-         }
-         else
-         {
-            mesh.GetElementFaces(e, faces, ori);
-         }
-
-         for (int f = 0; f < faces.Size(); f++)
-         {
-            const FiniteElement *fe = n_fespace->GetFaceElement(faces[f]);
-            const IntegrationRule ir = fe->GetNodes();
-
-            DenseMatrix vals;
-            DenseMatrix tr;
-            nodes->GetFaceVectorValues(faces[f], 0, ir, vals, tr);
-            xyz.Append(vals.GetData(), vals.Height()*vals.Width());
-         }
-      }
-
-      int nptface = xyz.Size()/dim;
-
-      // Generate points inside each element
-      FiniteElementCollection *l2_fec = new L2_FECollection(l2_order, dim);
-      FiniteElementSpace l2_fespace =
-         FiniteElementSpace(&mesh, l2_fec, 1);
-      DenseMatrix vals;
-      DenseMatrix tr;
-      for (int e = 0; e < mesh.GetNE(); e++)
-      {
-         const FiniteElement *fe = l2_fespace.GetFE(e);
-         const IntegrationRule ir = fe->GetNodes();
-
-         nodes->GetVectorValues(e, ir, vals, tr);
-         xyz.Append(vals.GetData(), vals.Height()*vals.Width());
-      }
-
-      Vector xyzv(xyz.GetData(), xyz.Size());
-      int npt = xyzv.Size()/dim;
-
-      FindPointsGSLIB finder;
-      finder.Setup(mesh);
-      finder.FindPoints(xyzv, Ordering::byVDIM);
-      Array<unsigned int> code_out    = finder.GetCode();
-      unsigned int cmin = 5,
-                   cmax = 0;
-      for (int i = 0; i < nptface; i++)
-      {
-         cmin = std::min(code_out[i], cmin);
-         cmax = std::max(code_out[i], cmax);
-      }
-      REQUIRE((cmin == 1 && cmax == 1)); // should be found on element boundary
-
-      cmin = 5;
-      cmax = 0;
-      for (int i = nptface; i < npt; i++)
-      {
-         cmin = std::min(code_out[i], cmin);
-         cmax = std::max(code_out[i], cmax);
-      }
-      REQUIRE((cmin == 0 && cmax == 0)); // should be found inside element
-      delete l2_fec;
-   }
-}
-
-// Generate a 4x4 Quad/Hex Mesh and interpolate point in the center of domain
-// at element boundary. This tests L2 projection with and without averaging.
-TEST_CASE("GSLIBInterpolateL2ElementBoundary",
-          "[GSLIBInterpolateL2ElementBoundary][GSLIB]")
-{
-   int dim  = GENERATE(2, 3);
-   CAPTURE(dim);
-
-   int nex = 4;
-   int mesh_order = 2;
-   Mesh mesh;
-   if (dim == 2)
-   {
-      mesh = Mesh::MakeCartesian2D(nex, nex, Element::QUADRILATERAL);
-   }
-   else
-   {
-      mesh = Mesh::MakeCartesian3D(nex, nex, nex, Element::HEXAHEDRON);
-   }
-
-   mesh.SetCurvature(mesh_order);
-
-   // Set GridFunction to be interpolated
-   int func_order = 3;
-   FiniteElementCollection *c_fec = new L2_FECollection(func_order, dim);
-   FiniteElementSpace c_fespace =
-      FiniteElementSpace(&mesh, c_fec, 1);
-   GridFunction field_vals(&c_fespace);
-   Array<int> dofs;
-   double leftval = 1.0;
-   double rightval = 3.0;
-   for (int e = 0; e < mesh.GetNE(); e++)
-   {
-      Vector center(dim);
-      mesh.GetElementCenter(e, center);
-      double val_to_set = center(0) < 0.5 ? leftval : rightval;
-      c_fespace.GetElementDofs(e, dofs);
-      Vector vals(dofs.Size());
-      vals = val_to_set;
-      field_vals.SetSubVector(dofs, vals);
-   }
-
-   int npt = 1;
-   Vector xyz(npt*dim);
-   xyz = 0.0;
-   xyz(0) = 0.5;
-
-   // Find and interpolate FE Function values
-   Vector interp_vals(npt);
-   FindPointsGSLIB finder;
-   finder.Setup(mesh);
-   finder.SetL2AvgType(FindPointsGSLIB::NONE);
-   finder.Interpolate(xyz, field_vals, interp_vals, 1);
-   Array<unsigned int> code_out    = finder.GetCode();
-
-   // This point should have been found on element border. But the interpolated
-   // value will come from either of the elements that share this edge/face.
-   REQUIRE(code_out[0] == 1);
-   REQUIRE((interp_vals(0) == MFEM_Approx(leftval) ||
-            interp_vals(0) == MFEM_Approx(rightval)));
-
-   // Interpolated value should now be average of solution coming from
-   // adjacent elements.
-   finder.SetL2AvgType(FindPointsGSLIB::ARITHMETIC);
-   finder.Interpolate(xyz, field_vals, interp_vals, 1);
-   REQUIRE(interp_vals(0) == MFEM_Approx(0.5*(leftval+rightval)));
-
-   finder.FreeData();
-   delete c_fec;
-}
-
-} //namespace_gslib
-#endif

tests/unit/fem/test_pa_coeff.cpp

@@ -265,248 +265,332 @@ TEST_CASE("H1 PA Coefficient", "[PartialAssembly][Coefficient]")
 TEST_CASE("Hcurl/Hdiv PA Coefficient",
           "[CUDA][PartialAssembly][Coefficient]")
 {
-   const bool all_tests = launch_all_non_regression_tests;
-   enum MixedSpaces {Hcurl, Hdiv, HcurlHdiv, HdivHcurl, NumSpaceTypes};
-   // coeff_type: 0 - ConstantCoefficient
-   //             1 - FunctionCoefficient
-   //             2 - VectorFunctionCoefficient
-   //             3 - SymmetricMatrixFunctionCoefficient
-   //             4 - MatrixFunctionCoefficient
-
-   dimension = GENERATE(2, 3);
-   const int order = all_tests ? GENERATE(1, 2, 3) : GENERATE(1, 2);
-   const int coeff_type = GENERATE(0, 1, 2, 3, 4); // see comment above
-   const MixedSpaces space_type = GENERATE(Hcurl, Hdiv, HcurlHdiv, HdivHcurl);
-   CAPTURE(space_type, dimension, coeff_type, order);
-
-   const int ne = 2;
-   Mesh mesh = MakeCartesianNonaligned(dimension, ne);
-
-   std::unique_ptr<Coefficient> coeff;
-   std::unique_ptr<Coefficient> coeff2;
-   std::unique_ptr<VectorCoefficient> vcoeff;
-   std::unique_ptr<MatrixCoefficient> mcoeff;
-
-   if (coeff_type == 0)
+   for (dimension = 2; dimension < 4; ++dimension)
    {
-      coeff.reset(new ConstantCoefficient(12.34));
-      coeff2.reset(new ConstantCoefficient(12.34));
-   }
-   else if (coeff_type == 1)
-   {
-      coeff.reset(new FunctionCoefficient(&coeffFunction));
-      coeff2.reset(new FunctionCoefficient(&linearFunction));
-   }
-   else if (coeff_type == 2)
-   {
-      vcoeff.reset(new VectorFunctionCoefficient(dimension, &vectorCoeffFunction));
-      coeff2.reset(new FunctionCoefficient(&linearFunction));
-   }
-   else if (coeff_type == 3)
-   {
-      mcoeff.reset(new SymmetricMatrixFunctionCoefficient(dimension,
-                                                          &symmetricMatrixCoeffFunction));
-      coeff2.reset(new FunctionCoefficient(&linearFunction));
-   }
-   else if (coeff_type == 4)
-   {
-      mcoeff.reset(new MatrixFunctionCoefficient(dimension,
-                                                 &asymmetricMatrixCoeffFunction));
-      coeff2.reset(new FunctionCoefficient(&linearFunction));
-   }
+      const int ne = 3;
+      Mesh mesh = MakeCartesianNonaligned(dimension, ne);
 
-   std::unique_ptr<FiniteElementCollection> fec;
-   if (space_type == Hcurl || space_type == HcurlHdiv)
-   {
-      fec.reset(new ND_FECollection(order, dimension));
-   }
-   else if (space_type == HdivHcurl)
-   {
-      fec.reset(new RT_FECollection(order - 1, dimension));
-   }
-   else
-   {
-      fec.reset(new RT_FECollection(order, dimension));
-   }
-
-   FiniteElementSpace fes(&mesh, fec.get());
-
-   // Set essential boundary conditions on the entire boundary.
-   Array<int> ess_tdof_list;
-   fes.GetBoundaryTrueDofs(ess_tdof_list);
-
-   Vector xin(fes.GetTrueVSize());
-   xin.Randomize(1);
-
-   Vector y_fa, y_pa;
-
-   if (space_type == HcurlHdiv || space_type == HdivHcurl)
-   {
-      std::unique_ptr<FiniteElementCollection> fec_test;
-      if (space_type == HcurlHdiv)
+      for (int coeffType = 0; coeffType < 5; ++coeffType)
       {
-         fec_test.reset(new RT_FECollection(order - 1, dimension));
-      }
-      else
-      {
-         fec_test.reset(new ND_FECollection(order, dimension));
-      }
-
-      FiniteElementSpace fes_test(&mesh, fec_test.get());
-
-      MixedBilinearForm pa_form(&fes, &fes_test);
-      pa_form.SetAssemblyLevel(AssemblyLevel::PARTIAL);
-      MixedBilinearForm fa_form(&fes, &fes_test);
-
-      const int ndof_test = fes_test.GetTrueVSize();
-      y_fa.SetSize(ndof_test);
-      y_pa.SetSize(ndof_test);
-
-      if (mcoeff)
-      {
-         pa_form.AddDomainIntegrator(new VectorFEMassIntegrator(*mcoeff));
-         fa_form.AddDomainIntegrator(new VectorFEMassIntegrator(*mcoeff));
-      }
-      else if (vcoeff)
-      {
-         pa_form.AddDomainIntegrator(new VectorFEMassIntegrator(*vcoeff));
-         fa_form.AddDomainIntegrator(new VectorFEMassIntegrator(*vcoeff));
-      }
-      else
-      {
-         pa_form.AddDomainIntegrator(new VectorFEMassIntegrator(*coeff));
-         fa_form.AddDomainIntegrator(new VectorFEMassIntegrator(*coeff));
-      }
-
-      if (dimension == 3)
-      {
-         if (vcoeff)
+         Coefficient* coeff = nullptr;
+         Coefficient* coeff2 = nullptr;
+         VectorCoefficient* vcoeff = nullptr;
+         MatrixCoefficient* mcoeff = nullptr;
+         if (coeffType == 0)
          {
-            if (space_type == HcurlHdiv)
+            coeff = new ConstantCoefficient(12.34);
+            coeff2 = new ConstantCoefficient(12.34);
+         }
+         else if (coeffType == 1)
+         {
+            coeff = new FunctionCoefficient(&coeffFunction);
+            coeff2 = new FunctionCoefficient(&linearFunction);
+         }
+         else if (coeffType == 2)
+         {
+            vcoeff = new VectorFunctionCoefficient(dimension, &vectorCoeffFunction);
+            coeff2 = new FunctionCoefficient(&linearFunction);
+         }
+         else if (coeffType == 3)
+         {
+            mcoeff = new SymmetricMatrixFunctionCoefficient(dimension,
+                                                            &symmetricMatrixCoeffFunction);
+            coeff2 = new FunctionCoefficient(&linearFunction);
+         }
+         else if (coeffType == 4)
+         {
+            mcoeff = new MatrixFunctionCoefficient(dimension,
+                                                   &asymmetricMatrixCoeffFunction);
+            coeff2 = new FunctionCoefficient(&linearFunction);
+         }
+
+         enum MixedSpaces {Hcurl, Hdiv, HcurlHdiv, HdivHcurl, NumSpaceTypes};
+
+         for (int spaceType = 0; spaceType < NumSpaceTypes; ++spaceType)
+         {
+            const int numIntegrators =
+               (spaceType >= HcurlHdiv) ? 1 : ((coeffType == 2) ? 2 : 3);
+
+            for (int integrator = 0; integrator < numIntegrators; ++integrator)
             {
-               pa_form.AddDomainIntegrator(new MixedVectorCurlIntegrator(*vcoeff));
-               fa_form.AddDomainIntegrator(new MixedVectorCurlIntegrator(*vcoeff));
-            }
-            else
-            {
-               pa_form.AddDomainIntegrator(new MixedVectorWeakCurlIntegrator(*vcoeff));
-               fa_form.AddDomainIntegrator(new MixedVectorWeakCurlIntegrator(*vcoeff));
-            }
-         }
-         else
-         {
-            if (space_type == HcurlHdiv)
-            {
-               pa_form.AddDomainIntegrator(new MixedVectorCurlIntegrator(*coeff));
-               fa_form.AddDomainIntegrator(new MixedVectorCurlIntegrator(*coeff));
-            }
-            else
-            {
-               pa_form.AddDomainIntegrator(new MixedVectorWeakCurlIntegrator(*coeff));
-               fa_form.AddDomainIntegrator(new MixedVectorWeakCurlIntegrator(*coeff));
+               for (int order = 1; order < 4; ++order)
+               {
+                  CAPTURE(spaceType, dimension, coeffType, integrator, order);
+
+                  FiniteElementCollection* fec = nullptr;
+                  if (spaceType == Hcurl || spaceType == HcurlHdiv)
+                  {
+                     fec = (FiniteElementCollection*) new ND_FECollection(order, dimension);
+                  }
+                  else if (spaceType == HdivHcurl)
+                  {
+                     fec = (FiniteElementCollection*) new RT_FECollection(order - 1, dimension);
+                  }
+                  else
+                  {
+                     fec = (FiniteElementCollection*) new RT_FECollection(order, dimension);
+                  }
+
+                  FiniteElementSpace fespace(&mesh, fec);
+
+                  // Set essential boundary conditions on the entire boundary.
+                  Array<int> tdof_ess(fespace.GetVSize());
+                  tdof_ess = 0;
+
+                  for (int i=0; i<mesh.GetNBE(); ++i)
+                  {
+                     Array<int> dofs;
+                     fespace.GetBdrElementDofs(i, dofs);
+                     for (int j=0; j<dofs.Size(); ++j)
+                     {
+                        const int dof_j = (dofs[j] >= 0) ? dofs[j] : -1 - dofs[j];
+                        tdof_ess[dof_j] = 1;
+                     }
+                  }
+
+                  int num_ess = 0;
+                  for (int i=0; i<fespace.GetVSize(); ++i)
+                  {
+                     if (tdof_ess[i] == 1)
+                     {
+                        num_ess++;
+                     }
+                  }
+
+                  Array<int> ess_tdof_list(num_ess);
+                  num_ess = 0;
+                  for (int i=0; i<fespace.GetVSize(); ++i)
+                  {
+                     if (tdof_ess[i] == 1)
+                     {
+                        ess_tdof_list[num_ess] = i;
+                        num_ess++;
+                     }
+                  }
+
+                  Vector xin(fespace.GetTrueVSize());
+                  xin.Randomize(1);
+
+                  Vector y_mat, y_assembly, y_pa;
+
+                  if (spaceType >= HcurlHdiv)
+                  {
+                     FiniteElementCollection* fecTest = nullptr;
+                     if (spaceType == HcurlHdiv)
+                     {
+                        fecTest = (FiniteElementCollection*) new RT_FECollection(order - 1, dimension);
+                     }
+                     else
+                     {
+                        fecTest = (FiniteElementCollection*) new ND_FECollection(order, dimension);
+                     }
+
+                     FiniteElementSpace fespaceTest(&mesh, fecTest);
+
+                     MixedBilinearForm *paform = new MixedBilinearForm(&fespace, &fespaceTest);
+                     paform->SetAssemblyLevel(AssemblyLevel::PARTIAL);
+                     MixedBilinearForm *assemblyform = new MixedBilinearForm(&fespace, &fespaceTest);
+
+                     const int testSize = fespaceTest.GetTrueVSize();
+                     y_mat.SetSize(testSize);
+                     y_mat = 0.0;
+                     y_assembly.SetSize(testSize);
+                     y_assembly = 0.0;
+                     y_pa.SetSize(testSize);
+                     y_pa = 0.0;
+
+                     if (coeffType >= 3)
+                     {
+                        paform->AddDomainIntegrator(new VectorFEMassIntegrator(*mcoeff));
+                        assemblyform->AddDomainIntegrator(new VectorFEMassIntegrator(*mcoeff));
+                     }
+                     else if (coeffType == 2)
+                     {
+                        paform->AddDomainIntegrator(new VectorFEMassIntegrator(*vcoeff));
+                        assemblyform->AddDomainIntegrator(new VectorFEMassIntegrator(*vcoeff));
+                     }
+                     else
+                     {
+                        paform->AddDomainIntegrator(new VectorFEMassIntegrator(*coeff));
+                        assemblyform->AddDomainIntegrator(new VectorFEMassIntegrator(*coeff));
+                     }
+
+                     if (dimension == 3 && (spaceType == HcurlHdiv || spaceType == HdivHcurl))
+                     {
+                        if (coeffType == 2)
+                        {
+                           if (spaceType == HcurlHdiv)
+                           {
+                              paform->AddDomainIntegrator(new MixedVectorCurlIntegrator(*vcoeff));
+                              assemblyform->AddDomainIntegrator(new MixedVectorCurlIntegrator(*vcoeff));
+                           }
+                           else
+                           {
+                              paform->AddDomainIntegrator(new MixedVectorWeakCurlIntegrator(*vcoeff));
+                              assemblyform->AddDomainIntegrator(new MixedVectorWeakCurlIntegrator(*vcoeff));
+                           }
+                        }
+                        else if (coeffType < 2)
+                        {
+                           if (spaceType == HcurlHdiv)
+                           {
+                              paform->AddDomainIntegrator(new MixedVectorCurlIntegrator(*coeff));
+                              assemblyform->AddDomainIntegrator(new MixedVectorCurlIntegrator(*coeff));
+                           }
+                           else
+                           {
+                              paform->AddDomainIntegrator(new MixedVectorWeakCurlIntegrator(*coeff));
+                              assemblyform->AddDomainIntegrator(new MixedVectorWeakCurlIntegrator(*coeff));
+                           }
+                        }
+                     }
+
+                     Array<int> empty_ess; // empty
+
+                     paform->Assemble();
+                     OperatorHandle paopr;
+                     paform->FormRectangularSystemMatrix(ess_tdof_list, empty_ess, paopr);
+
+                     assemblyform->Assemble();
+                     assemblyform->Finalize();
+
+                     OperatorPtr A_explicit;
+                     assemblyform->FormRectangularSystemMatrix(ess_tdof_list, empty_ess, A_explicit);
+
+                     paopr->Mult(xin, y_pa);
+                     assemblyform->Mult(xin, y_assembly);
+                     A_explicit->Mult(xin, y_mat);
+
+                     // Test the transpose
+                     if ((spaceType == HcurlHdiv || spaceType == HdivHcurl) &&
+                         dimension == 3)
+                     {
+                        Vector u(testSize);
+                        u.Randomize();
+
+                        Vector v_mat(fespace.GetTrueVSize());
+                        v_mat = 0.0;
+                        Vector v_assembly(fespace.GetTrueVSize());
+                        v_assembly = 0.0;
+                        Vector v_pa(fespace.GetTrueVSize());
+                        v_pa = 0.0;
+
+                        const SparseMatrix& A_spmat = assemblyform->SpMat();
+                        paopr->MultTranspose(u, v_pa);
+                        assemblyform->MultTranspose(u, v_assembly);
+                        A_spmat.MultTranspose(u, v_mat);
+
+                        v_pa -= v_mat;
+                        double pa_error = v_pa.Norml2();
+                        REQUIRE(pa_error < 1.e-12);
+
+                        v_assembly -= v_mat;
+                        double assembly_error = v_assembly.Norml2();
+                        REQUIRE(assembly_error < 1.e-12);
+                     }
+
+                     delete paform;
+                     delete assemblyform;
+                     delete fecTest;
+                  }
+                  else
+                  {
+                     BilinearForm *paform = new BilinearForm(&fespace);
+                     paform->SetAssemblyLevel(AssemblyLevel::PARTIAL);
+                     BilinearForm *assemblyform = new BilinearForm(&fespace);
+
+                     y_mat.SetSize(xin.Size());
+                     y_mat = 0.0;
+                     y_assembly.SetSize(xin.Size());
+                     y_assembly = 0.0;
+                     y_pa.SetSize(xin.Size());
+                     y_pa = 0.0;
+
+                     if (integrator < 2)
+                     {
+                        if (coeffType >= 3)
+                        {
+                           paform->AddDomainIntegrator(new VectorFEMassIntegrator(*mcoeff));
+                           assemblyform->AddDomainIntegrator(new VectorFEMassIntegrator(*mcoeff));
+                        }
+                        else if (coeffType == 2)
+                        {
+                           paform->AddDomainIntegrator(new VectorFEMassIntegrator(*vcoeff));
+                           assemblyform->AddDomainIntegrator(new VectorFEMassIntegrator(*vcoeff));
+
+                        }
+                        else
+                        {
+                           paform->AddDomainIntegrator(new VectorFEMassIntegrator(*coeff));
+                           assemblyform->AddDomainIntegrator(new VectorFEMassIntegrator(*coeff));
+                        }
+                     }
+                     if (integrator > 0)
+                     {
+                        if (spaceType == Hcurl)
+                        {
+                           const FiniteElement *fel = fespace.GetFE(0);
+                           const IntegrationRule *intRule = &MassIntegrator::GetRule(*fel, *fel,
+                                                                                     *mesh.GetElementTransformation(0));
+
+                           if (coeffType >= 3 && dimension == 3)
+                           {
+                              paform->AddDomainIntegrator(new CurlCurlIntegrator(*mcoeff, intRule));
+                              assemblyform->AddDomainIntegrator(new CurlCurlIntegrator(*mcoeff, intRule));
+                           }
+                           else if (coeffType == 2 && dimension == 3)
+                           {
+                              paform->AddDomainIntegrator(new CurlCurlIntegrator(*vcoeff, intRule));
+                              assemblyform->AddDomainIntegrator(new CurlCurlIntegrator(*vcoeff, intRule));
+                           }
+                           else
+                           {
+                              paform->AddDomainIntegrator(new CurlCurlIntegrator(*coeff2));
+                              assemblyform->AddDomainIntegrator(new CurlCurlIntegrator(*coeff2));
+                           }
+                        }
+                        else
+                        {
+                           paform->AddDomainIntegrator(new DivDivIntegrator(*coeff2));
+                           assemblyform->AddDomainIntegrator(new DivDivIntegrator(*coeff2));
+                        }
+                     }
+                     paform->Assemble();
+                     OperatorHandle paopr;
+                     paform->FormSystemMatrix(ess_tdof_list, paopr);
+
+                     assemblyform->SetDiagonalPolicy(Matrix::DIAG_ONE);
+                     assemblyform->Assemble();
+                     OperatorPtr A_explicit;
+                     assemblyform->FormSystemMatrix(ess_tdof_list, A_explicit);
+
+                     paopr->Mult(xin, y_pa);
+                     assemblyform->Mult(xin, y_assembly);
+                     A_explicit->Mult(xin, y_mat);
+
+                     delete paform;
+                     delete assemblyform;
+                  }
+
+                  y_pa -= y_mat;
+                  double pa_error = y_pa.Norml2();
+                  REQUIRE(pa_error < 1.e-10);
+
+                  y_assembly -= y_mat;
+                  double assembly_error = y_assembly.Norml2();
+                  REQUIRE(assembly_error < 1.e-12);
+
+                  delete fec;
+               }
             }
          }
+
+         delete coeff;
+         delete coeff2;
+         delete vcoeff;
+         delete mcoeff;
       }
-
-      Array<int> empty_ess; // empty
-
-      OperatorHandle pa_op;
-      pa_form.Assemble();
-      pa_form.FormRectangularSystemMatrix(ess_tdof_list, empty_ess, pa_op);
-
-      OperatorPtr fa_op;
-      fa_form.Assemble();
-      fa_form.Finalize();
-      fa_form.FormRectangularSystemMatrix(ess_tdof_list, empty_ess, fa_op);
-
-      // Test the transpose
-      if (dimension == 3)
-      {
-         Vector u(ndof_test);
-         u.Randomize();
-
-         Vector v_pa(fes.GetTrueVSize());
-         Vector v_fa(fes.GetTrueVSize());
-
-         pa_op->MultTranspose(u, v_pa);
-         fa_op->MultTranspose(u, v_fa);
-
-         v_pa -= v_fa;
-         REQUIRE(v_pa.Norml2() == MFEM_Approx(0.0));
-      }
-
-      pa_op->Mult(xin, y_pa);
-      fa_op->Mult(xin, y_fa);
    }
-   else
-   {
-      BilinearForm pa_form(&fes);
-      pa_form.SetAssemblyLevel(AssemblyLevel::PARTIAL);
-      BilinearForm fa_form(&fes);
-
-      y_fa.SetSize(xin.Size());
-      y_pa.SetSize(xin.Size());
-
-      if (mcoeff)
-      {
-         pa_form.AddDomainIntegrator(new VectorFEMassIntegrator(*mcoeff));
-         fa_form.AddDomainIntegrator(new VectorFEMassIntegrator(*mcoeff));
-      }
-      else if (vcoeff)
-      {
-         pa_form.AddDomainIntegrator(new VectorFEMassIntegrator(*vcoeff));
-         fa_form.AddDomainIntegrator(new VectorFEMassIntegrator(*vcoeff));
-
-      }
-      else
-      {
-         pa_form.AddDomainIntegrator(new VectorFEMassIntegrator(*coeff));
-         fa_form.AddDomainIntegrator(new VectorFEMassIntegrator(*coeff));
-      }
-
-      if (space_type == Hcurl)
-      {
-         const FiniteElement *fel = fes.GetFE(0);
-         const IntegrationRule &ir =
-            MassIntegrator::GetRule(*fel, *fel, *mesh.GetElementTransformation(0));
-
-         if (coeff_type >= 3 && dimension == 3)
-         {
-            pa_form.AddDomainIntegrator(new CurlCurlIntegrator(*mcoeff, &ir));
-            fa_form.AddDomainIntegrator(new CurlCurlIntegrator(*mcoeff, &ir));
-         }
-         else if (coeff_type == 2 && dimension == 3)
-         {
-            pa_form.AddDomainIntegrator(new CurlCurlIntegrator(*vcoeff, &ir));
-            fa_form.AddDomainIntegrator(new CurlCurlIntegrator(*vcoeff, &ir));
-         }
-         else
-         {
-            pa_form.AddDomainIntegrator(new CurlCurlIntegrator(*coeff2));
-            fa_form.AddDomainIntegrator(new CurlCurlIntegrator(*coeff2));
-         }
-      }
-      else // space_type == Hdiv
-      {
-         pa_form.AddDomainIntegrator(new DivDivIntegrator(*coeff2));
-         fa_form.AddDomainIntegrator(new DivDivIntegrator(*coeff2));
-      }
-
-      OperatorHandle pa_op;
-      pa_form.Assemble();
-      pa_form.FormSystemMatrix(ess_tdof_list, pa_op);
-
-      OperatorPtr fa_op;
-      fa_form.SetDiagonalPolicy(Matrix::DIAG_ONE);
-      fa_form.Assemble();
-      fa_form.FormSystemMatrix(ess_tdof_list, fa_op);
-
-      pa_op->Mult(xin, y_pa);
-      fa_op->Mult(xin, y_fa);
-   }
-
-   y_pa -= y_fa;
-   REQUIRE(y_pa.Norml2() == MFEM_Approx(0.0, 1e-10));
 }
 
 TEST_CASE("Hcurl/Hdiv Mixed PA Coefficient",

tests/unit/fem/test_pa_kernels.cpp

@@ -385,17 +385,17 @@ void test_pa_convection(const std::string &meshname, int order, int prob,
    }
    int dim = mesh.Dimension();
 
-   std::unique_ptr<FiniteElementCollection> fec;
+   FiniteElementCollection *fec;
    if (prob)
    {
       auto basis = prob==3 ? BasisType::Positive : BasisType::GaussLobatto;
-      fec.reset(new L2_FECollection(order, dim, basis));
+      fec = new L2_FECollection(order, dim, basis);
    }
    else
    {
-      fec.reset(new H1_FECollection(order, dim));
+      fec = new H1_FECollection(order, dim);
    }
-   FiniteElementSpace fespace(&mesh, fec.get());
+   FiniteElementSpace fespace(&mesh, fec);
 
    L2_FECollection vel_fec(order, dim, BasisType::GaussLobatto);
    FiniteElementSpace vel_fespace(&mesh, &vel_fec, dim);
@@ -405,21 +405,21 @@ void test_pa_convection(const std::string &meshname, int order, int prob,
    BilinearForm k_pa(&fespace);
    BilinearForm k_fa(&fespace);
 
-   std::unique_ptr<VectorCoefficient> vel_coeff;
-   std::unique_ptr<Coefficient> rho;
+   VectorCoefficient *vel_coeff;
+   Coefficient *rho;
 
    // prob: 0: CG, 1: DG continuous coeff, 2: DG discontinuous coeff
    if (prob >= 2)
    {
       vel_gf.Randomize(1);
-      vel_coeff.reset(new VectorGridFunctionCoefficient(&vel_gf));
+      vel_coeff = new VectorGridFunctionCoefficient(&vel_gf);
       rho_gf.Randomize(1);
-      rho.reset(new GridFunctionCoefficient(&rho_gf));
+      rho = new GridFunctionCoefficient(&rho_gf);
    }
    else
    {
-      vel_coeff.reset(new VectorFunctionCoefficient(dim, velocity_function));
-      rho.reset(new ConstantCoefficient(1.0));
+      vel_coeff = new VectorFunctionCoefficient(dim, velocity_function);
+      rho = new ConstantCoefficient(1.0);
    }
 
 
@@ -451,6 +451,10 @@ void test_pa_convection(const std::string &meshname, int order, int prob,
    y_pa -= y_fa;
 
    REQUIRE(y_pa.Norml2() < 1.e-12);
+
+   delete vel_coeff;
+   delete rho;
+   delete fec;
 }
 
 // Basic unit tests for convection
@@ -464,7 +468,7 @@ TEST_CASE("PA Convection", "[PartialAssembly], [CUDA]")
    auto prob = GENERATE(0, 1, 2, 3);
    auto order = GENERATE(2);
    // refinement > 0 => Non-conforming mesh
-   auto refinement = GENERATE(0, 1);
+   auto refinement = GENERATE(0,1);
 
    SECTION("2D")
    {

tests/unit/mesh/test_face_orientations.cpp

@@ -10,10 +10,10 @@
 // CONTRIBUTING.md for details.
 
 #include "mfem.hpp"
-#include "unit_tests.hpp"
-
 using namespace mfem;
 
+#include "unit_tests.hpp"
+
 class TestMesh : public Mesh
 {
 public:
@@ -117,97 +117,3 @@ TEST_CASE("Face Orientation", "[FaceOrientation]")
       }
    }
 }
-
-template <Geometry::Type geom_t>
-constexpr Geometry::Type GetFaceType();
-
-template <>
-constexpr Geometry::Type GetFaceType<Geometry::SEGMENT>()
-{
-   return Geometry::POINT;
-}
-
-template <>
-constexpr Geometry::Type GetFaceType<Geometry::TRIANGLE>()
-{
-   return Geometry::SEGMENT;
-}
-
-template <>
-constexpr Geometry::Type GetFaceType<Geometry::SQUARE>()
-{
-   return Geometry::SEGMENT;
-}
-
-template <>
-constexpr Geometry::Type GetFaceType<Geometry::TETRAHEDRON>()
-{
-   return Geometry::TRIANGLE;
-}
-
-template <>
-constexpr Geometry::Type GetFaceType<Geometry::CUBE>()
-{
-   return Geometry::SQUARE;
-}
-
-TEMPLATE_TEST_CASE_SIG("Boundary Element Face Orientation", "[FaceOrientation]",
-                       ((Geometry::Type geom_t), geom_t),
-                       Geometry::SEGMENT, Geometry::TRIANGLE, Geometry::SQUARE,
-                       Geometry::TETRAHEDRON, Geometry::CUBE)
-{
-   constexpr auto face_t = GetFaceType<geom_t>();
-   using face_t_consts = Geometry::Constants<face_t>;
-
-   Mesh mesh;
-   constexpr int n1d = 1;
-   switch (geom_t)
-   {
-      case Geometry::SEGMENT:
-         mesh = Mesh::MakeCartesian1D(n1d, Element::SEGMENT);
-         break;
-      case Geometry::TRIANGLE:
-         mesh = Mesh::MakeCartesian2D(n1d, n1d, Element::TRIANGLE);
-         break;
-      case Geometry::SQUARE:
-         mesh = Mesh::MakeCartesian2D(n1d, n1d, Element::QUADRILATERAL);
-         break;
-      case Geometry::TETRAHEDRON:
-         mesh = Mesh::MakeCartesian3D(n1d, n1d, n1d, Element::TETRAHEDRON);
-         break;
-      case Geometry::CUBE:
-         mesh = Mesh::MakeCartesian3D(n1d, n1d, n1d, Element::HEXAHEDRON);
-         break;
-      default:
-         MFEM_ABORT("");
-   }
-   Element *be0 = mesh.GetBdrElement(0);
-   MFEM_VERIFY(be0->GetGeometryType() == face_t, "");
-   int f, o;
-   mesh.GetBdrElementFace(0, &f, &o);
-   const Element *face = mesh.GetFace(f);
-   int *be0_v = be0->GetVertices();
-   const int *face_v = face->GetVertices();
-   for (o = 0; o < face_t_consts::NumOrient; o++)
-   {
-      const int *face_perm = face_t_consts::Orient[o];
-      for (int i = 0; i < face_t_consts::NumVert; i++)
-      {
-         be0_v[i] = face_v[face_perm[i]];
-      }
-      IsoparametricTransformation bdr_tr, face_tr;
-      mesh.GetBdrElementTransformation(0, &bdr_tr);
-      mesh.GetFaceTransformation(f, &face_tr);
-      IntegrationPoint bdr_ip;
-      bdr_ip.Set3(0.1, 0.3, 0.0);
-      int inv_o;
-      mesh.GetBdrElementFace(0, &f, &inv_o);
-      MFEM_VERIFY(inv_o == face_t_consts::InvOrient[o], "");
-      IntegrationPoint face_ip = Mesh::TransformBdrElementToFace(
-                                    be0->GetGeometryType(), inv_o, bdr_ip);
-      Vector bdr_pt, face_pt;
-      bdr_tr.Transform(bdr_ip, bdr_pt);
-      face_tr.Transform(face_ip, face_pt);
-      REQUIRE(bdr_pt.DistanceTo(face_pt) == MFEM_Approx(0.0));
-   }
-}

tests/unit/mesh/test_fms.cpp

@@ -10,10 +10,10 @@
 // CONTRIBUTING.md for details.
 
 #include "mfem.hpp"
-#include "unit_tests.hpp"
-
 using namespace mfem;
 
+#include "unit_tests.hpp"
+
 #ifdef MFEM_USE_FMS
 TEST_CASE("Load FMS datacollection", "[FMS]")
 {

tests/unit/mesh/test_mesh.cpp

@@ -10,10 +10,10 @@
 // CONTRIBUTING.md for details.
 
 #include "mfem.hpp"
-#include "unit_tests.hpp"
-
 using namespace mfem;
 
+#include "unit_tests.hpp"
+
 TEST_CASE("Element-wise construction", "[Mesh]")
 {
    SECTION("Quadrilateral")

tests/unit/mesh/test_ncmesh.cpp

@@ -12,7 +12,6 @@
 #include "mfem.hpp"
 #include "unit_tests.hpp"
 
-#include <array>
 namespace mfem
 {
 
@@ -300,8 +299,8 @@ TEST_CASE("pNCMesh PA diagonal",  "[Parallel], [NCMesh]")
 
 // Given a parallel and a serial mesh, perform an L2 projection and check the
 // solutions match exactly.
-std::array<double, 2> CheckL2Projection(ParMesh& pmesh, Mesh& smesh, int order,
-                                        std::function<double(Vector const&)> exact_soln)
+void CheckL2Projection(ParMesh& pmesh, Mesh& smesh, int order,
+                       std::function<double(Vector const&)> exact_soln)
 {
    REQUIRE(pmesh.GetGlobalNE() == smesh.GetNE());
    REQUIRE(pmesh.Dimension() == smesh.Dimension());
@@ -369,6 +368,7 @@ std::array<double, 2> CheckL2Projection(ParMesh& pmesh, Mesh& smesh, int order,
       x = 0.0;
 
       double pnorm = x.ComputeL2Error(rhs_coef);
+
       b.AddDomainIntegrator(new DomainLFIntegrator(rhs_coef));
       b.Assemble();
 
@@ -393,11 +393,26 @@ std::array<double, 2> CheckL2Projection(ParMesh& pmesh, Mesh& smesh, int order,
       return x.ComputeL2Error(rhs_coef) / pnorm;
    }();
 
-   return {serror, perror};
+   constexpr double test_tol = 1e-9;
+   CHECK(std::abs(serror - perror) < test_tol);
 };
 
-TEST_CASE("EdgeFaceConstraint",  "[Parallel], [NCMesh]")
+
+TEST_CASE("FaceEdgeConstraint",  "[Parallel], [NCMesh]")
 {
+   constexpr int refining_rank = 0;
+   auto smesh = Mesh("../../data/ref-tetrahedron.mesh");
+
+   REQUIRE(smesh.GetNE() == 1);
+   {
+      // Start the test with two tetrahedra attached by triangle.
+      auto single_edge_refine = Array<Refinement>(1);
+      single_edge_refine[0].index = 0;
+      single_edge_refine[0].ref_type = Refinement::X;
+
+      smesh.GeneralRefinement(single_edge_refine, 0); // conformal
+   }
+
    auto exact_soln = [](const Vector& x)
    {
       // sin(|| x - d ||^2) -> non polynomial but very smooth.
@@ -407,239 +422,103 @@ TEST_CASE("EdgeFaceConstraint",  "[Parallel], [NCMesh]")
       return std::sin(d * d);
    };
 
-   SECTION("ReferenceTet")
+   REQUIRE(smesh.GetNE() == 2);
+   smesh.EnsureNCMesh(true);
+   smesh.Finalize();
+
+   auto partition = std::unique_ptr<int[]>(new int[smesh.GetNE()]);
+   partition[0] = 0;
+   partition[1] = Mpi::WorldSize() > 1 ? 1 : 0;
+
+   auto pmesh = ParMesh(MPI_COMM_WORLD, smesh, partition.get());
+
+   // Construct the NC refined mesh in parallel and serial. Once constructed a
+   // global L2 projected solution should match exactly on each.
+   Array<int> refines, serial_refines(1);
+   if (Mpi::WorldRank() == refining_rank)
    {
-      constexpr int refining_rank = 0;
-      auto smesh = Mesh("../../data/ref-tetrahedron.mesh");
+      refines.Append(0);
+   }
 
-      REQUIRE(smesh.GetNE() == 1);
-      {
-         // Start the test with two tetrahedra attached by triangle.
-         auto single_edge_refine = Array<Refinement>(1);
-         single_edge_refine[0].index = 0;
-         single_edge_refine[0].ref_type = Refinement::X;
+   // Must be called on all ranks as it uses MPI calls internally.
+   // All ranks will use the global element number dictated by rank 0 though.
+   serial_refines[0] = pmesh.GetGlobalElementNum(0);
+   MPI_Bcast(&serial_refines[0], 1, MPI_INT, refining_rank, MPI_COMM_WORLD);
 
-         smesh.GeneralRefinement(single_edge_refine, 0); // conformal
-      }
+   // Rank 0 refines the parallel mesh, all ranks refine the serial mesh
+   smesh.GeneralRefinement(serial_refines, 1); // nonconformal
+   pmesh.GeneralRefinement(refines, 1); // nonconformal
 
+   REQUIRE(pmesh.GetGlobalNE() == 8 + 1);
+   REQUIRE(smesh.GetNE() == 8 + 1);
 
-      REQUIRE(smesh.GetNE() == 2);
-      smesh.EnsureNCMesh(true);
-      smesh.Finalize();
+   // Each pair of indices here represents sequential element indices to refine.
+   // First the i element is refined, then in the resulting mesh the j element is
+   // refined. These pairs were arrived at by looping over all possible i,j pairs and
+   // checking for the addition of a face-edge constraint.
+   std::vector<std::pair<int,int>> indices{{2,13}, {3,13}, {6,2}, {6,3}};
 
-      auto partition = std::unique_ptr<int[]>(new int[smesh.GetNE()]);
-      partition[0] = 0;
-      partition[1] = Mpi::WorldSize() > 1 ? 1 : 0;
-
-      auto pmesh = ParMesh(MPI_COMM_WORLD, smesh, partition.get());
-
-      // Construct the NC refined mesh in parallel and serial. Once constructed a
-      // global L2 projected solution should match exactly on each.
-      Array<int> refines, serial_refines(1);
+   // Rank 0 has all but one element in the parallel mesh. The remaining element
+   // is owned by another processor if the number of ranks is greater than one.
+   for (const auto &ij : indices)
+   {
+      int i = ij.first;
+      int j = ij.second;
       if (Mpi::WorldRank() == refining_rank)
       {
-         refines.Append(0);
+         refines[0] = i;
       }
+      // Inform all ranks of the serial mesh
+      serial_refines[0] = pmesh.GetGlobalElementNum(i);
+      MPI_Bcast(&serial_refines[0], 1, MPI_INT, 0, MPI_COMM_WORLD);
 
-      // Must be called on all ranks as it uses MPI calls internally.
-      // All ranks will use the global element number dictated by rank 0 though.
-      serial_refines[0] = pmesh.GetGlobalElementNum(0);
-      MPI_Bcast(&serial_refines[0], 1, MPI_INT, refining_rank, MPI_COMM_WORLD);
+      ParMesh tmp(pmesh);
+      tmp.GeneralRefinement(refines);
 
-      // Rank 0 refines the parallel mesh, all ranks refine the serial mesh
-      smesh.GeneralRefinement(serial_refines, 1); // nonconformal
-      pmesh.GeneralRefinement(refines, 1); // nonconformal
+      REQUIRE(tmp.GetGlobalNE() == 1 + 8 - 1 + 8); // 16 elements
 
-      REQUIRE(pmesh.GetGlobalNE() == 8 + 1);
-      REQUIRE(smesh.GetNE() == 8 + 1);
+      Mesh stmp(smesh);
+      stmp.GeneralRefinement(serial_refines);
+      REQUIRE(stmp.GetNE() == 1 + 8 - 1 + 8); // 16 elements
 
-      // Each pair of indices here represents sequential element indices to refine.
-      // First the i element is refined, then in the resulting mesh the j element is
-      // refined. These pairs were arrived at by looping over all possible i,j pairs and
-      // checking for the addition of a face-edge constraint.
-      std::vector<std::pair<int,int>> indices{{2,13}, {3,13}, {6,2}, {6,3}};
-
-      // Rank 0 has all but one element in the parallel mesh. The remaining element
-      // is owned by another processor if the number of ranks is greater than one.
-      for (const auto &ij : indices)
+      if (Mpi::WorldRank() == refining_rank)
       {
-         int i = ij.first;
-         int j = ij.second;
-         if (Mpi::WorldRank() == refining_rank)
-         {
-            refines[0] = i;
-         }
-         // Inform all ranks of the serial mesh
-         serial_refines[0] = pmesh.GetGlobalElementNum(i);
-         MPI_Bcast(&serial_refines[0], 1, MPI_INT, 0, MPI_COMM_WORLD);
-
-         ParMesh tmp(pmesh);
-         tmp.GeneralRefinement(refines);
-
-         REQUIRE(tmp.GetGlobalNE() == 1 + 8 - 1 + 8); // 16 elements
-
-         Mesh stmp(smesh);
-         stmp.GeneralRefinement(serial_refines);
-         REQUIRE(stmp.GetNE() == 1 + 8 - 1 + 8); // 16 elements
-
-         if (Mpi::WorldRank() == refining_rank)
-         {
-            refines[0] = j;
-         }
-         // Inform all ranks of the serial mesh
-         serial_refines[0] = tmp.GetGlobalElementNum(j);
-         MPI_Bcast(&serial_refines[0], 1, MPI_INT, 0, MPI_COMM_WORLD);
-
-         ParMesh ttmp(tmp);
-         ttmp.GeneralRefinement(refines);
-
-         REQUIRE(ttmp.GetGlobalNE() == 1 + 8 - 1 + 8 - 1 + 8); // 23 elements
-
-         Mesh sttmp(stmp);
-         sttmp.GeneralRefinement(serial_refines);
-         REQUIRE(sttmp.GetNE() == 1 + 8 - 1 + 8 - 1 + 8); // 23 elements
-
-         // Loop over interior faces, fill and check face transform on the serial.
-         for (int iface = 0; iface < sttmp.GetNumFaces(); ++iface)
-         {
-            const auto face_transform = sttmp.GetFaceElementTransformations(iface);
-            CHECK(face_transform->CheckConsistency(0) < 1e-12);
-         }
-
-         for (int iface = 0; iface < ttmp.GetNumFacesWithGhost(); ++iface)
-         {
-            const auto face_transform = ttmp.GetFaceElementTransformations(iface);
-            CHECK(face_transform->CheckConsistency(0) < 1e-12);
-         }
-
-         // Use P4 to ensure there's a few fully interior DOF.
-         {
-            auto error = CheckL2Projection(ttmp, sttmp, 4, exact_soln);
-            double constexpr tol = 1e-9;
-            CHECK(std::abs(error[1] - error[0]) < tol);
-         }
-         ttmp.ExchangeFaceNbrData();
-         ttmp.Rebalance();
-         {
-            auto error = CheckL2Projection(ttmp, sttmp, 4, exact_soln);
-            double constexpr tol = 1e-9;
-            CHECK(std::abs(error[1] - error[0]) < tol);
-         }
+         refines[0] = j;
       }
+      // Inform all ranks of the serial mesh
+      serial_refines[0] = tmp.GetGlobalElementNum(j);
+      MPI_Bcast(&serial_refines[0], 1, MPI_INT, 0, MPI_COMM_WORLD);
+
+      ParMesh ttmp(tmp);
+      ttmp.GeneralRefinement(refines);
+
+      REQUIRE(ttmp.GetGlobalNE() == 1 + 8 - 1 + 8 - 1 + 8); // 23 elements
+
+      Mesh sttmp(stmp);
+      sttmp.GeneralRefinement(serial_refines);
+      REQUIRE(sttmp.GetNE() == 1 + 8 - 1 + 8 - 1 + 8); // 23 elements
+
+      // Loop over interior faces, fill and check face transform on the serial.
+      for (int iface = 0; iface < sttmp.GetNumFaces(); ++iface)
+      {
+         const auto face_transform = sttmp.GetFaceElementTransformations(iface);
+         CHECK(face_transform->CheckConsistency(0) < 1e-12);
+      }
+
+      for (int iface = 0; iface < ttmp.GetNumFacesWithGhost(); ++iface)
+      {
+         const auto face_transform = ttmp.GetFaceElementTransformations(iface);
+         CHECK(face_transform->CheckConsistency(0) < 1e-12);
+      }
+
+      // Use P4 to ensure there's a few fully interior DOF.
+      CheckL2Projection(ttmp, sttmp, 4, exact_soln);
+
+      ttmp.ExchangeFaceNbrData();
+      ttmp.Rebalance();
+
+      CheckL2Projection(ttmp, sttmp, 4, exact_soln);
    }
-
-   auto CheckSerialParallelH1Equivalence = [](Mesh &smesh)
-   {
-      constexpr int dim = 3;
-      constexpr int order = 2;
-      H1_FECollection nd_fec(order, dim);
-      FiniteElementSpace fes(&smesh, &nd_fec);
-      const auto serial_ntdof = fes.GetTrueVSize();
-
-      ParMesh mesh(MPI_COMM_WORLD, smesh);
-      ParFiniteElementSpace pfes(&mesh, &nd_fec);
-      const auto parallel_ntdof = pfes.GlobalTrueVSize();
-
-      // If nc constraints have been observed correctly, the number of true dof in
-      // parallel should match the number of true dof in serial. If the number of
-      // parallel dofs is greater, then a slave constraint has not been fully labeled.
-      CHECK(serial_ntdof == parallel_ntdof);
-   };
-
-   auto CheckSerialParallelNDEquivalence = [](Mesh &smesh)
-   {
-      constexpr int dim = 3;
-      constexpr int order = 1;
-      ND_FECollection nd_fec(order, dim);
-      FiniteElementSpace fes(&smesh, &nd_fec);
-      const auto serial_ntdof = fes.GetTrueVSize();
-
-      ParMesh mesh(MPI_COMM_WORLD, smesh);
-      ParFiniteElementSpace pfes(&mesh, &nd_fec);
-      const auto parallel_ntdof = pfes.GlobalTrueVSize();
-
-      // If nc constraints have been observed correctly, the number of true dof in
-      // parallel should match the number of true dof in serial. If the number of
-      // parallel dofs is greater, then a slave constraint has not been fully labeled.
-      CHECK(serial_ntdof == parallel_ntdof);
-   };
-
-   SECTION("LevelTwoRefinement")
-   {
-      Mesh smesh("../../data/ref-tetrahedron.mesh");
-      Array<Refinement> aniso_ref(1);
-      aniso_ref[0].index = 0;
-      aniso_ref[0].ref_type = Refinement::X;
-      smesh.GeneralRefinement(aniso_ref);
-      smesh.UniformRefinement();
-      smesh.EnsureNCMesh(true);
-      Array<int> el_to_refine(1);
-
-      for (int n = 0; n < smesh.GetNE(); n++)
-      {
-         Mesh smesh2(smesh);
-         el_to_refine[0] = n;
-         smesh2.GeneralRefinement(el_to_refine);
-         for (int m = 0; m < smesh2.GetNE(); m++)
-         {
-            Mesh smesh3(smesh2);
-            el_to_refine[0] = m;
-            smesh3.GeneralRefinement(el_to_refine);
-            CAPTURE(n,m);
-            CheckSerialParallelNDEquivalence(smesh3);
-            CheckSerialParallelH1Equivalence(smesh3);
-         }
-      }
-   }
-
-   SECTION("EdgeCasePartition")
-   {
-      Mesh smesh("../../data/ref-tetrahedron.mesh");
-      smesh.UniformRefinement();
-      smesh.EnsureNCMesh(true);
-      Array<int> el_to_refine(1);
-
-      el_to_refine[0] = 0;
-      smesh.GeneralRefinement(el_to_refine);
-
-      // This particular partition was found by brute force search. The default rebalancing
-      // can in rare cases produce similar local patterns, particularly for highly adapted meshes.
-      auto partition = std::unique_ptr<int[]>(new int[smesh.GetNE()]);
-      if (Mpi::WorldSize() > 1)
-      {
-         auto bad_partition = std::vector<int> {0, 0, 1, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0};
-         std::copy(bad_partition.begin(), bad_partition.end(), partition.get());
-      }
-      else
-      {
-         for (int i = 0; i < smesh.GetNE(); i++)
-         {
-            partition[i] = 0;
-         }
-      }
-      ParMesh pmesh(MPI_COMM_WORLD, smesh, partition.get());
-
-      {
-         constexpr int dim = 3;
-         constexpr int order = 1;
-         ND_FECollection nd_fec(order, dim);
-         FiniteElementSpace fes(&smesh, &nd_fec);
-         const auto serial_ntdof = fes.GetTrueVSize();
-         ParFiniteElementSpace pfes(&pmesh, &nd_fec);
-         pfes.ExchangeFaceNbrData();
-         const auto parallel_ntdof = pfes.GlobalTrueVSize();
-         CHECK(serial_ntdof == parallel_ntdof);
-      }
-
-      for (int order = 1; order <= 4; order++)
-      {
-         CAPTURE(order);
-         auto error = CheckL2Projection(pmesh, smesh, order, exact_soln);
-         double constexpr tol = 1e-9;
-         CHECK(std::abs(error[1] - error[0]) < tol);
-      }
-   }
-
 } // test case
 
 Mesh CylinderMesh(Geometry::Type el_type, bool quadratic, int variant = 0)
@@ -837,8 +716,7 @@ TEST_CASE("P2Q1PureTetHexPri",  "[Parallel], [NCMesh]")
       auto pmesh = ParMesh(MPI_COMM_WORLD, smesh);
 
       // P2 ensures there are triangles without dofs
-      auto error = CheckL2Projection(pmesh, smesh, 2, exact_soln);
-      CHECK(std::abs(error[1] - error[0]) < 1e-9);
+      CheckL2Projection(pmesh, smesh, 2, exact_soln);
    }
 } // test case
 
@@ -881,153 +759,11 @@ TEST_CASE("PNQ2PureTetHexPri",  "[Parallel], [NCMesh]")
 
       for (int p = 1; p < 3; ++p)
       {
-         auto error = CheckL2Projection(pmesh, smesh, p, exact_soln);
-         CHECK(std::abs(error[1] - error[0]) < 1e-9);
+         CheckL2Projection(pmesh, smesh, p, exact_soln);
       }
    }
 } // test case
 
-/**
- * @brief Test GetVectorValue on face neighbor elements for nonconformal meshes
- *
- * @param smesh The serial mesh to start from
- * @param nc_level Depth of refinement on processor boundaries
- * @param skip Refine every "skip" processor boundary element
- * @param use_ND Whether to use Nedelec elements (which are sensitive to orientation)
- */
-void TestVectorValueInVolume(Mesh &smesh, int nc_level, int skip, bool use_ND)
-{
-   auto vector_exact_soln = [](const Vector& x, Vector& v)
-   {
-      Vector d(3);
-      d[0] = -0.5; d[1] = -1; d[2] = -2; // arbitrary
-      v = (d -= x);
-   };
-
-   smesh.Finalize();
-   smesh.EnsureNCMesh(true);
-
-   auto pmesh = ParMesh(MPI_COMM_WORLD, smesh);
-
-   // Apply refinement on face neighbors to achieve a given nc level mismatch.
-   for (int i = 0; i < nc_level; ++i)
-   {
-      // To refine the face neighbors, need to know where they are.
-      pmesh.ExchangeFaceNbrData();
-      Array<int> elem_to_refine;
-      // Refine only on odd ranks.
-      if ((Mpi::WorldRank() + 1) % 2 == 0)
-      {
-         // Refine a subset of all shared faces. Using a subset helps to
-         // mix in conformal faces with nonconformal faces.
-         for (int n = 0; n < pmesh.GetNSharedFaces(); ++n)
-         {
-            if (n % skip != 0) { continue; }
-            const int local_face = pmesh.GetSharedFace(n);
-            const auto &face_info = pmesh.GetFaceInformation(local_face);
-            REQUIRE(face_info.IsShared());
-            REQUIRE(face_info.element[1].location == Mesh::ElementLocation::FaceNbr);
-            elem_to_refine.Append(face_info.element[0].index);
-         }
-      }
-      pmesh.GeneralRefinement(elem_to_refine);
-   }
-
-   // Do not rebalance again! The test is also checking for nc refinements
-   // along the processor boundary.
-
-   // Create a grid function of the mesh coordinates
-   pmesh.ExchangeFaceNbrData();
-   pmesh.EnsureNodes();
-   REQUIRE(pmesh.OwnsNodes());
-   GridFunction * const coords = pmesh.GetNodes();
-   dynamic_cast<ParGridFunction *>(pmesh.GetNodes())->ExchangeFaceNbrData();
-
-   // Project the linear function onto the mesh. Quadratic ND tetrahedral
-   // elements are the first to require face orientations.
-   const int order = 2, dim = 3;
-   std::unique_ptr<FiniteElementCollection> fec;
-   if (use_ND)
-   {
-      fec = std::unique_ptr<ND_FECollection>(new ND_FECollection(order, dim));
-   }
-   else
-   {
-      fec = std::unique_ptr<RT_FECollection>(new RT_FECollection(order, dim));
-   }
-   ParFiniteElementSpace pnd_fes(&pmesh, fec.get());
-
-   ParGridFunction psol(&pnd_fes);
-
-   VectorFunctionCoefficient func(3, vector_exact_soln);
-   psol.ProjectCoefficient(func);
-   psol.ExchangeFaceNbrData();
-
-   mfem::Vector value(3), exact(3), position(3);
-   const IntegrationRule &ir = mfem::IntRules.Get(Geometry::Type::TETRAHEDRON,
-                                                  order + 1);
-
-   // Check that non-ghost elements match up on the serial and parallel spaces.
-   for (int n = 0; n < pmesh.GetNE(); ++n)
-   {
-      constexpr double tol = 1e-12;
-      for (const auto &ip : ir)
-      {
-         coords->GetVectorValue(n, ip, position);
-         psol.GetVectorValue(n, ip, value);
-
-         vector_exact_soln(position, exact);
-
-         REQUIRE(value.Size() == exact.Size());
-         CHECK((value -= exact).Normlinf() < tol);
-      }
-   }
-
-   // Loop over face neighbor elements and check the vector values match in the
-   // face neighbor elements.
-   for (int n = 0; n < pmesh.GetNSharedFaces(); ++n)
-   {
-      const int local_face = pmesh.GetSharedFace(n);
-      const auto &face_info = pmesh.GetFaceInformation(local_face);
-      REQUIRE(face_info.IsShared());
-      REQUIRE(face_info.element[1].location == Mesh::ElementLocation::FaceNbr);
-
-      auto &T = *pmesh.GetFaceNbrElementTransformation(face_info.element[1].index);
-
-      constexpr double tol = 1e-12;
-      for (const auto &ip : ir)
-      {
-         T.SetIntPoint(&ip);
-         coords->GetVectorValue(T, ip, position);
-         psol.GetVectorValue(T, ip, value);
-
-         vector_exact_soln(position, exact);
-
-         REQUIRE(value.Size() == exact.Size());
-         CHECK((value -= exact).Normlinf() < tol);
-      }
-   }
-}
-
-TEST_CASE("GetVectorValueInFaceNeighborElement", "[Parallel], [NCMesh]")
-{
-   // The aim of this test is to verify the correct behaviour of the
-   // GetVectorValue method when called on face neighbor elements in a non
-   // conforming mesh.
-   auto smesh = Mesh("../../data/beam-tet.mesh");
-
-   for (int nc_level : {0,1,2,3})
-   {
-      for (int skip : {1,2})
-      {
-         for (bool use_ND : {false, true})
-         {
-            TestVectorValueInVolume(smesh, nc_level, skip, use_ND);
-         }
-      }
-   }
-}
-
 #endif // MFEM_USE_MPI
 
 } // namespace mfem