Rule of zero for GridFunction and ParGridFunction

Use shared_ptr to handle the possibly owned finite element collection and space.

.gitignore

@@ -266,6 +266,7 @@ 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,6 +11,15 @@
 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
 ===========================================
@@ -31,6 +40,7 @@ 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,16 +331,30 @@ 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 (currently not needed)
+# Ginkgo library configuration
 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_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)))
+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)))
 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)
@@ -349,7 +363,8 @@ ifeq ($(GINKGO_BUILD_TYPE),Debug)
   endif
 else
 endif
-GINKGO_LIB = $(XLINKER)-rpath,$(GINKGO_LIB_DIR) -L$(GINKGO_LIB_DIR) $(GINKGO_LINK)
+GINKGO_LIB = $(XLINKER)-rpath,$(GINKGO_LINK_LIB_DIR) -L$(GINKGO_LINK_LIB_DIR)\
+ $(GINKGO_LINK)
 
 # AmgX library configuration
 AMGX_DIR = @MFEM_DIR@/../amgx

examples/ex13p.cpp

@@ -5,6 +5,7 @@
 // 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
@@ -54,6 +55,7 @@ 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";
 
@@ -69,6 +71,9 @@ 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.");
@@ -98,6 +103,10 @@ 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,6 +13,7 @@
 //               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/ex3.cpp

@@ -63,6 +63,7 @@ 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;
 
@@ -77,6 +78,9 @@ 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",
@@ -102,6 +106,11 @@ 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/ex36.cpp

@@ -194,7 +194,7 @@ int main(int argc, char *argv[])
    u_gf.ProjectCoefficient(IC_coef);
    u_old_gf = u_gf;
 
-   // 9. Initialize the slack variable ψₕ = exp(uₕ)
+   // 9. Initialize the slack variable ψₕ = ln(uₕ)
    LogarithmGridFunctionCoefficient ln_u(u_gf, obstacle);
    psi_gf.ProjectCoefficient(ln_u);
    psi_old_gf = psi_gf;

examples/ex36p.cpp

@@ -226,7 +226,7 @@ int main(int argc, char *argv[])
    u_gf.ProjectCoefficient(IC_coef);
    u_old_gf = u_gf;
 
-   // 9. Initialize the slack variable ψₕ = exp(uₕ)
+   // 9. Initialize the slack variable ψₕ = ln(uₕ)
    LogarithmGridFunctionCoefficient ln_u(u_gf, obstacle);
    psi_gf.ProjectCoefficient(ln_u);
    psi_old_gf = psi_gf;

examples/ex3p.cpp

@@ -5,6 +5,7 @@
 // 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
@@ -70,6 +71,7 @@ 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
@@ -87,6 +89,9 @@ 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",
@@ -124,6 +129,11 @@ 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,7 +450,8 @@ int main(int argc, char *argv[])
 
 // Implementation of class FE_Evolution
 FE_Evolution::FE_Evolution(BilinearForm &M_, BilinearForm &K_, const Vector &b_)
-   : TimeDependentOperator(M_.Height()), M(M_), K(K_), b(b_), z(M_.Height())
+   : TimeDependentOperator(M_.FESpace()->GetTrueVSize()),
+     M(M_), K(K_), b(b_), z(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_.Height()), b(b_),
+   : TimeDependentOperator(M_.ParFESpace()->GetTrueVSize()), b(b_),
      M_solver(M_.ParFESpace()->GetComm()),
-     z(M_.Height())
+     z(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 BlockMatrix class, as
+//               The example demonstrates the use of the BlockOperator 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_.Height(), 0.0,
+   : TimeDependentOperator(M_.ParFESpace()->GetTrueVSize(), 0.0,
                            M_in_lhs ? TimeDependentOperator::IMPLICIT
                            : TimeDependentOperator::EXPLICIT),
-     b(b_), comm(M_.ParFESpace()->GetComm()), M_solver(comm), z(M_.Height()),
+     b(b_), comm(M_.ParFESpace()->GetComm()), M_solver(comm), z(height),
      iJacobian(NULL), rJacobian(NULL)
 {
    MAlev = M_.GetAssemblyLevel();

examples/sundials/ex9.cpp

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

fem/bilininteg.cpp

@@ -2456,7 +2456,7 @@ void VectorFEMassIntegrator::AssembleElementMatrix(
 {
    int dof = el.GetDof();
    int spaceDim = Trans.GetSpaceDim();
-   int vdim = std::max(spaceDim, el.GetVDim());
+   int vdim = std::max(spaceDim, el.GetRangeDim());
 
    double w;
 
@@ -2524,7 +2524,7 @@ void VectorFEMassIntegrator::AssembleElementMatrix2(
    {
       // assume test_fe is scalar FE and trial_fe is vector FE
       int spaceDim = Trans.GetSpaceDim();
-      int vdim = std::max(spaceDim, trial_fe.GetVDim());
+      int vdim = std::max(spaceDim, trial_fe.GetRangeDim());
       int trial_dof = trial_fe.GetDof();
       int test_dof = test_fe.GetDof();
       double w;
@@ -2622,8 +2622,8 @@ void VectorFEMassIntegrator::AssembleElementMatrix2(
    {
       // assume both test_fe and trial_fe are vector FE
       int spaceDim = Trans.GetSpaceDim();
-      int trial_vdim = std::max(spaceDim, trial_fe.GetVDim());
-      int test_vdim = std::max(spaceDim, test_fe.GetVDim());
+      int trial_vdim = std::max(spaceDim, trial_fe.GetRangeDim());
+      int test_vdim = std::max(spaceDim, test_fe.GetRangeDim());
       int trial_dof = trial_fe.GetDof();
       int test_dof = test_fe.GetDof();
       double w;

fem/bilininteg.hpp

@@ -580,7 +580,7 @@ protected:
 
 
    inline virtual int GetTestVDim(const FiniteElement & test_fe)
-   { return std::max(space_dim, test_fe.GetVDim()); }
+   { return std::max(space_dim, test_fe.GetRangeDim()); }
 
    inline virtual void CalcTestShape(const FiniteElement & test_fe,
                                      ElementTransformation &Trans,
@@ -588,7 +588,7 @@ protected:
    { test_fe.CalcVShape(Trans, shape); }
 
    inline virtual int GetTrialVDim(const FiniteElement & trial_fe)
-   { return std::max(space_dim, trial_fe.GetVDim()); }
+   { return std::max(space_dim, trial_fe.GetRangeDim()); }
 
    inline virtual void CalcTrialShape(const FiniteElement & trial_fe,
                                       ElementTransformation &Trans,
@@ -678,7 +678,7 @@ protected:
 
 
    inline virtual int GetVDim(const FiniteElement & vector_fe)
-   { return std::max(space_dim, vector_fe.GetVDim()); }
+   { return std::max(space_dim, vector_fe.GetRangeDim()); }
 
    inline virtual void CalcVShape(const FiniteElement & vector_fe,
                                   ElementTransformation &Trans,
@@ -1105,7 +1105,7 @@ public:
       const FiniteElement & trial_fe,
       const FiniteElement & test_fe) const
    {
-      return (trial_fe.GetVDim() == 3 &&
+      return (trial_fe.GetRangeDim() == 3 &&
               trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR &&
               test_fe.GetRangeType()  == mfem::FiniteElement::SCALAR &&
               test_fe.GetDerivType()  == mfem::FiniteElement::GRAD );
@@ -1288,8 +1288,8 @@ public:
       const FiniteElement & trial_fe,
       const FiniteElement & test_fe) const
    {
-      return (trial_fe.GetCurlDim() == 3 && trial_fe.GetVDim() == 3 &&
-              test_fe.GetCurlDim() == 3 && test_fe.GetVDim() == 3 &&
+      return (trial_fe.GetCurlDim() == 3 && trial_fe.GetRangeDim() == 3 &&
+              test_fe.GetCurlDim() == 3 && test_fe.GetRangeDim() == 3 &&
               trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR &&
               trial_fe.GetDerivType() == mfem::FiniteElement::CURL &&
               test_fe.GetRangeType()  == mfem::FiniteElement::VECTOR &&
@@ -1419,7 +1419,7 @@ public:
       const FiniteElement & trial_fe,
       const FiniteElement & test_fe) const
    {
-      return (trial_fe.GetVDim() == 3 && test_fe.GetCurlDim() == 3 &&
+      return (trial_fe.GetRangeDim() == 3 && test_fe.GetCurlDim() == 3 &&
               trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR &&
               test_fe.GetRangeType()  == mfem::FiniteElement::VECTOR &&
               test_fe.GetDerivType()  == mfem::FiniteElement::CURL );
@@ -1489,7 +1489,7 @@ public:
       const FiniteElement & trial_fe,
       const FiniteElement & test_fe) const
    {
-      return (test_fe.GetVDim() == 3 &&
+      return (test_fe.GetRangeDim() == 3 &&
               trial_fe.GetRangeType() == mfem::FiniteElement::SCALAR &&
               trial_fe.GetDerivType() == mfem::FiniteElement::GRAD &&
               test_fe.GetRangeType()  == mfem::FiniteElement::VECTOR );
@@ -1529,7 +1529,7 @@ public:
       const FiniteElement & trial_fe,
       const FiniteElement & test_fe) const
    {
-      return (trial_fe.GetCurlDim() == 3 && test_fe.GetVDim() == 3 &&
+      return (trial_fe.GetCurlDim() == 3 && test_fe.GetRangeDim() == 3 &&
               trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR &&
               trial_fe.GetDerivType() == mfem::FiniteElement::CURL   &&
               test_fe.GetRangeType()  == mfem::FiniteElement::VECTOR );
@@ -1900,7 +1900,7 @@ protected:
       const FiniteElement & trial_fe,
       const FiniteElement & test_fe) const
    {
-      return (trial_fe.GetCurlDim() == 3 && test_fe.GetVDim() == 3 &&
+      return (trial_fe.GetCurlDim() == 3 && test_fe.GetRangeDim() == 3 &&
               trial_fe.GetDerivType() == mfem::FiniteElement::CURL  &&
               test_fe.GetRangeType()  == mfem::FiniteElement::VECTOR );
    }
@@ -1959,7 +1959,7 @@ protected:
       const FiniteElement & trial_fe,
       const FiniteElement & test_fe) const
    {
-      return (trial_fe.GetVDim() == 3 && test_fe.GetCurlDim() == 3 &&
+      return (trial_fe.GetRangeDim() == 3 && test_fe.GetCurlDim() == 3 &&
               trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR &&
               test_fe.GetDerivType()  == mfem::FiniteElement::CURL );
    }

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()));
    }
 }
 
@@ -1592,6 +1592,10 @@ void VectorQuadratureFunctionCoefficient::Eval(Vector &V,
    QuadF.HostRead();
 
    const int el_idx = QuadF.GetSpace()->GetEntityIndex(T);
+   // Handle the case of "interior boundary elements" and FaceQuadratureSpace
+   // with FaceType::Boundary.
+   if (el_idx < 0) { V = 0.0; return; }
+
    const int ip_idx = QuadF.GetSpace()->GetPermutedIndex(el_idx, ip.index);
 
    if (index == 0 && vdim == QuadF.GetVDim())
@@ -1629,6 +1633,9 @@ double QuadratureFunctionCoefficient::Eval(ElementTransformation &T,
    QuadF.HostRead();
    Vector temp(1);
    const int el_idx = QuadF.GetSpace()->GetEntityIndex(T);
+   // Handle the case of "interior boundary elements" and FaceQuadratureSpace
+   // with FaceType::Boundary.
+   if (el_idx < 0) { return 0.0; }
    const int ip_idx = QuadF.GetSpace()->GetPermutedIndex(el_idx, ip.index);
    QuadF.GetValues(el_idx, ip_idx, temp);
    return temp[0];

fem/coefficient.hpp

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

fem/eltrans.cpp

@@ -492,6 +492,7 @@ 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/fe/fe_base.cpp

@@ -807,7 +807,7 @@ void NodalFiniteElement::Project(
    else
    {
       DenseMatrix vshape(fe.GetDof(), std::max(Trans.GetSpaceDim(),
-                                               fe.GetVDim()));
+                                               fe.GetRangeDim()));
 
       I.SetSize(vshape.Width()*dof, fe.GetDof());
       for (int k = 0; k < dof; k++)

fem/fe/fe_base.hpp

@@ -307,19 +307,20 @@ public:
    FiniteElement(int D, Geometry::Type G, int Do, int O,
                  int F = FunctionSpace::Pk);
 
-   /// Returns the reference space dimension for the finite element
+   /// Returns the reference space dimension for the finite element.
    int GetDim() const { return dim; }
 
-   /// Returns the vector dimension for vector-valued finite elements
-   int GetVDim() const { return vdim; }
+   /** @brief Returns the vector dimension for vector-valued finite elements,
+       which is also the dimension of the interpolation operatrion. */
+   int GetRangeDim() const { return vdim; }
 
-   /// Returns the dimension of the curl for vector-valued finite elements
+   /// Returns the dimension of the curl for vector-valued finite elements.
    int GetCurlDim() const { return cdim; }
 
-   /// Returns the Geometry::Type of the reference element
+   /// Returns the Geometry::Type of the reference element.
    Geometry::Type GetGeomType() const { return geom_type; }
 
-   /// Returns the number of degrees of freedom in the finite element
+   /// Returns the number of degrees of freedom in the finite element.
    int GetDof() const { return dof; }
 
    /** @brief Returns the order of the finite element. In the case of

fem/fe/fe_nd.cpp

@@ -1852,7 +1852,7 @@ void ND_R1D_SegmentElement::Project(const FiniteElement &fe,
    else
    {
       double vk[Geometry::MaxDim];
-      DenseMatrix vshape(fe.GetDof(), fe.GetVDim());
+      DenseMatrix vshape(fe.GetDof(), fe.GetRangeDim());
 
       double * tk_ptr = const_cast<double*>(tk);
 
@@ -2293,7 +2293,7 @@ void ND_R2D_FiniteElement::Project(const FiniteElement &fe,
    else
    {
       double vk[Geometry::MaxDim];
-      DenseMatrix vshape(fe.GetDof(), fe.GetVDim());
+      DenseMatrix vshape(fe.GetDof(), fe.GetRangeDim());
 
       double * tk_ptr = const_cast<double*>(tk);
 

fem/fe/fe_rt.cpp

@@ -1486,7 +1486,7 @@ void RT_R1D_SegmentElement::Project(const FiniteElement &fe,
    else
    {
       double vk[Geometry::MaxDim];
-      DenseMatrix vshape(fe.GetDof(), fe.GetVDim());
+      DenseMatrix vshape(fe.GetDof(), fe.GetRangeDim());
 
       double * nk_ptr = const_cast<double*>(nk);
 
@@ -1523,7 +1523,7 @@ void RT_R1D_SegmentElement::ProjectCurl(const FiniteElement &fe,
                                         ElementTransformation &Trans,
                                         DenseMatrix &curl) const
 {
-   DenseMatrix curl_shape(fe.GetDof(), fe.GetVDim());
+   DenseMatrix curl_shape(fe.GetDof(), fe.GetRangeDim());
    Vector curl_k(fe.GetDof());
 
    double * nk_ptr = const_cast<double*>(nk);
@@ -1849,7 +1849,7 @@ void RT_R2D_FiniteElement::Project(const FiniteElement &fe,
    else
    {
       double vk[Geometry::MaxDim];
-      DenseMatrix vshape(fe.GetDof(), fe.GetVDim());
+      DenseMatrix vshape(fe.GetDof(), fe.GetRangeDim());
 
       double * nk_ptr = const_cast<double*>(nk);
 
@@ -1888,7 +1888,7 @@ void RT_R2D_FiniteElement::ProjectCurl(const FiniteElement &fe,
                                        ElementTransformation &Trans,
                                        DenseMatrix &curl) const
 {
-   DenseMatrix curl_shape(fe.GetDof(), fe.GetVDim());
+   DenseMatrix curl_shape(fe.GetDof(), fe.GetRangeDim());
    Vector curl_k(fe.GetDof());
 
    double * nk_ptr = const_cast<double*>(nk);

fem/fe_coll.cpp

@@ -87,6 +87,16 @@ int FiniteElementCollection::GetDerivMapType(int dim) const
    return FiniteElement::UNKNOWN_MAP_TYPE;
 }
 
+int FiniteElementCollection::GetRangeDim(int dim) const
+{
+   const FiniteElement *fe = FiniteElementForDim(dim);
+   if (fe)
+   {
+      return fe->GetRangeDim();
+   }
+   return 0;
+}
+
 int FiniteElementCollection::HasFaceDofs(Geometry::Type geom, int p) const
 {
    switch (geom)

fem/fespace.cpp

@@ -2428,6 +2428,7 @@ 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.
       }
    }
 
@@ -2446,6 +2447,7 @@ 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.
       }
    }
 
@@ -2656,7 +2658,6 @@ 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); }
          }
       }
@@ -2667,7 +2668,6 @@ int FiniteElementSpace::MakeDofTable(int ent_dim,
 
    // build the table
    entity_dofs.MakeFromList(num_ent+1, list);
-
    return total_dofs;
 }
 
@@ -2996,7 +2996,14 @@ 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, "");
+      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();
+      }());
    }
    else
    {

fem/gridfunc.cpp

@@ -38,8 +38,9 @@ GridFunction::GridFunction(Mesh *m, std::istream &input)
    // Grid functions are stored on the device
    UseDevice(true);
 
-   fes = new FiniteElementSpace;
-   fec = fes->Load(m, input);
+   owned_fes.reset(new FiniteElementSpace);
+   fes = owned_fes.get();
+   fec.reset(fes->Load(m, input));
 
    skip_comment_lines(input, '#');
    istream::int_type next_char = input.peek();
@@ -81,10 +82,11 @@ GridFunction::GridFunction(Mesh *m, GridFunction *gf_array[], int num_pieces)
    int vdim, ordering;
 
    fes = gf_array[0]->FESpace();
-   fec = FiniteElementCollection::New(fes->FEColl()->Name());
+   fec.reset(FiniteElementCollection::New(fes->FEColl()->Name()));
    vdim = fes->GetVDim();
    ordering = fes->GetOrdering();
-   fes = new FiniteElementSpace(m, fec, vdim, ordering);
+   owned_fes.reset(new FiniteElementSpace(m, fec.get(), vdim, ordering));
+   fes = owned_fes.get();
    SetSize(fes->GetVSize());
 
    if (m->NURBSext)
@@ -153,12 +155,9 @@ GridFunction::GridFunction(Mesh *m, GridFunction *gf_array[], int num_pieces)
 
 void GridFunction::Destroy()
 {
-   if (fec)
-   {
-      delete fes;
-      delete fec;
-      fec = NULL;
-   }
+   owned_fes.reset();
+   fec.reset();
+   fes = nullptr;
 }
 
 void GridFunction::Update()
@@ -340,7 +339,7 @@ int GridFunction::VectorDim() const
       return fes->GetVDim();
    }
    return fes->GetVDim()*std::max(fes->GetMesh()->SpaceDimension(),
-                                  fe->GetVDim());
+                                  fe->GetRangeDim());
 }
 
 int GridFunction::CurlDim() const
@@ -720,56 +719,6 @@ 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
@@ -834,18 +783,15 @@ 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 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);
+            int f, o;
+            fes->GetMesh()->GetBdrElementFace(T.ElementNo, &f, &o);
+            IntegrationPoint fip =
+               Mesh::TransformBdrElementToFace(FET->GetGeometryType(), o, ip);
 
             // Compute and set the point in element 1 from fip
             FET->SetAllIntPoints(&fip);
@@ -973,18 +919,15 @@ 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 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);
+            int f, o;
+            fes->GetMesh()->GetBdrElementFace(T.ElementNo, &f, &o);
+            IntegrationPoint fip =
+               Mesh::TransformBdrElementToFace(FET->GetGeometryType(), o, ip);
 
             // Compute and set the point in element 1 from fip
             FET->SetAllIntPoints(&fip);
@@ -997,6 +940,8 @@ 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
@@ -1041,7 +986,7 @@ void GridFunction::GetVectorValue(ElementTransformation &T,
    else
    {
       int spaceDim = fes->GetMesh()->SpaceDimension();
-      int vdim = std::max(spaceDim, fe->GetVDim());
+      int vdim = std::max(spaceDim, fe->GetRangeDim());
       DenseMatrix vshape(dof, vdim);
       fe->CalcVShape(T, vshape);
       val.SetSize(vdim);
@@ -1093,7 +1038,7 @@ void GridFunction::GetVectorValues(ElementTransformation &T,
    else
    {
       int spaceDim = fes->GetMesh()->SpaceDimension();
-      int vdim = std::max(spaceDim, FElem->GetVDim());
+      int vdim = std::max(spaceDim, FElem->GetRangeDim());
       DenseMatrix vshape(dof, vdim);
 
       vals.SetSize(vdim, nip);
@@ -1115,11 +1060,10 @@ int GridFunction::GetFaceVectorValues(
    int i, int side, const IntegrationRule &ir,
    DenseMatrix &vals, DenseMatrix &tr) const
 {
-   int n, di;
+   int di;
    FaceElementTransformations *Transf;
 
-   n = ir.GetNPoints();
-   IntegrationRule eir(n);  // ---
+   IntegrationRule eir(ir.GetNPoints());  // ---
    Transf = fes->GetMesh()->GetFaceElementTransformations(i, 0);
    if (side == 2)
    {
@@ -1141,12 +1085,14 @@ 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);
    }
@@ -1504,17 +1450,13 @@ 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 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);
+         int f, o;
+         fes->GetMesh()->GetBdrElementFace(T.ElementNo, &f, &o);
+         IntegrationPoint fip =
+            Mesh::TransformBdrElementToFace(FET->GetGeometryType(), o,
+                                            T.GetIntPoint());
 
          // Compute and set the point in element 1 from fip
          FET->SetAllIntPoints(&fip);
@@ -1601,17 +1543,13 @@ 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 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);
+         int f, o;
+         fes->GetMesh()->GetBdrElementFace(T.ElementNo, &f, &o);
+         IntegrationPoint fip =
+            Mesh::TransformBdrElementToFace(FET->GetGeometryType(), o,
+                                            T.GetIntPoint());
 
          // Compute and set the point in element 1 from fip
          FET->SetAllIntPoints(&fip);
@@ -1670,17 +1608,13 @@ 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 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);
+         int f, o;
+         fes->GetMesh()->GetBdrElementFace(T.ElementNo, &f, &o);
+         IntegrationPoint fip =
+            Mesh::TransformBdrElementToFace(FET->GetGeometryType(), o,
+                                            T.GetIntPoint());
 
          // Compute and set the point in element 1 from fip
          FET->SetAllIntPoints(&fip);
@@ -1756,17 +1690,13 @@ 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 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);
+         int f, o;
+         fes->GetMesh()->GetBdrElementFace(T.ElementNo, &f, &o);
+         IntegrationPoint fip =
+            Mesh::TransformBdrElementToFace(FET->GetGeometryType(), o,
+                                            T.GetIntPoint());
 
          // Compute and set the point in element 1 from fip
          FET->SetAllIntPoints(&fip);

fem/gridfunc.hpp

@@ -20,6 +20,7 @@
 #include "../general/adios2stream.hpp"
 #endif
 #include <limits>
+#include <memory>
 #include <ostream>
 #include <string>
 
@@ -30,14 +31,13 @@ namespace mfem
 class GridFunction : public Vector
 {
 protected:
-   /// FE space on which the grid function lives. Owned if #fec is not NULL.
+   /// FE space on which the grid function lives.
    FiniteElementSpace *fes;
 
    /** @brief Used when the grid function is read from a file. It can also be
-       set explicitly, see MakeOwner().
-
-       If not NULL, this pointer is owned by the GridFunction. */
-   FiniteElementCollection *fec;
+       set explicitly, see MakeOwner(). */
+   std::shared_ptr<FiniteElementCollection> fec;
+   std::shared_ptr<FiniteElementSpace> owned_fes;
 
    long fes_sequence; // see FiniteElementSpace::sequence, Mesh::sequence
 
@@ -74,11 +74,6 @@ 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); }
@@ -107,21 +102,12 @@ 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 = fec_; }
+   void MakeOwner(FiniteElementCollection *fec_) { fec.reset(fec_); }
 
-   FiniteElementCollection *OwnFEC() { return fec; }
+   FiniteElementCollection *OwnFEC() { return fec.get(); }
 
    int VectorDim() const;
    int CurlDim() const;
@@ -754,9 +740,6 @@ 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,6 +10,7 @@
 // CONTRIBUTING.md for details.
 
 #include "gslib.hpp"
+#include "geom.hpp"
 
 #ifdef MFEM_USE_GSLIB
 
@@ -238,7 +239,8 @@ 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.
+   // both simplices and quads. Also sets code to 1 for points found on element
+   // faces/edges.
    MapRefPosAndElemIndices();
 }
 
@@ -681,6 +683,9 @@ 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;
 
@@ -694,7 +699,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; };
+   struct out_pt { double r[3]; uint index, el, proc, code; };
    struct out_pt *pt;
    array_init(struct out_pt, outpt, nptsend);
    outpt->n=nptsend;
@@ -712,12 +717,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;
@@ -731,7 +736,13 @@ void FindPointsGSLIB::MapRefPosAndElemIndices()
       const Geometry::Type gt = fe->GetGeomType();
       pt->el = mesh_elem;
 
-      if (gt == Geometry::SQUARE || gt == Geometry::CUBE) { ++pt; continue; }
+      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;
+      }
       else if (gt == Geometry::TRIANGLE)
       {
          gf_rst_map_temp = gf_rst_map[0];
@@ -758,6 +769,10 @@ 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;
    }
 
@@ -774,6 +789,7 @@ void FindPointsGSLIB::MapRefPosAndElemIndices()
       {
          gsl_mfem_ref(d + pt->index*dim) = pt->r[d];
       }
+      gsl_code[pt->index] = pt->code;
       ++pt;
    }
    array_free(outpt);
@@ -784,12 +800,22 @@ 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) { continue; }
+         if (gt == Geometry::SQUARE || gt == Geometry::CUBE)
+         {
+            gsl_code[index] = Geometry::CheckPoint(gt, ip, -btol) ? 0 : 1;
+            continue;
+         }
          else if (gt == Geometry::TRIANGLE)
          {
             gf_rst_map_temp = gf_rst_map[0];
@@ -808,11 +834,12 @@ void FindPointsGSLIB::MapRefPosAndElemIndices()
          }
 
          int local_elem = split_element_index[elem];
-         IntegrationPoint ip;
-         Vector mfem_ref(gsl_mfem_ref.GetData()+index*dim, dim);
+         gf_rst_map_temp->GetVectorValue(local_elem, ip, mfem_ref);
+
+         // Check if the point is on element boundary
          ip.Set2(mfem_ref.GetData());
          if (dim == 3) { ip.z = mfem_ref(2); }
-         gf_rst_map_temp->GetVectorValue(local_elem, ip, mfem_ref);
+         gsl_code[index]  = Geometry::CheckPoint(gt, ip, -btol) ? 0 : 1;
       }
    }
 }

fem/integ/bilininteg_mass_pa.cpp

@@ -128,7 +128,7 @@ void MassIntegrator::AssemblePABoundary(const FiniteElementSpace &fes)
 
    int map_type = el.GetMapType();
    dim = el.GetDim(); // Dimension of the boundary element, *not* the mesh
-   ne = fes.GetMesh()->GetNBE();
+   ne = fes.GetMesh()->GetNFbyType(FaceType::Boundary);
    nq = ir->GetNPoints();
    face_geom = mesh->GetFaceGeometricFactors(*ir, GeometricFactors::DETERMINANTS,
                                              FaceType::Boundary, mt);

fem/intrules.cpp

@@ -1782,6 +1782,11 @@ IntegrationRule *IntegrationRules::PrismIntegrationRule(int Order)
    AllocIntRule(PrismIntRules, Order);
    PrismIntRules[Order] = new IntegrationRule(nt * ns);
    PrismIntRules[Order]->SetOrder(std::min(irt.GetOrder(), irs.GetOrder()));
+   while (Order < std::min(irt.GetOrder(), irs.GetOrder()))
+   {
+      AllocIntRule(PrismIntRules, ++Order);
+      PrismIntRules[Order] = PrismIntRules[Order-1];
+   }
 
    for (int ks=0; ks<ns; ks++)
    {

fem/linearform.cpp

@@ -27,15 +27,19 @@ LinearForm::LinearForm(FiniteElementSpace *f, LinearForm *lf)
    // Linear forms are stored on the device
    UseDevice(true);
 
-   // Copy the pointers to the integrators
+   // Copy the pointers to the integrators and the corresponding marker arrays
    domain_integs = lf->domain_integs;
+   domain_integs_marker = lf->domain_integs_marker;
 
    domain_delta_integs = lf->domain_delta_integs;
 
    boundary_integs = lf->boundary_integs;
+   boundary_integs_marker = lf->boundary_integs_marker;
 
    boundary_face_integs = lf->boundary_face_integs;
    boundary_face_integs_marker = lf->boundary_face_integs_marker;
+
+   interior_face_integs = lf->interior_face_integs;
 }
 
 void LinearForm::AddDomainIntegrator(LinearFormIntegrator *lfi)

fem/lininteg.cpp

@@ -395,8 +395,8 @@ void VectorBoundaryLFIntegrator::AssembleRHSElementVect(
    {
       const IntegrationPoint &ip = ir->IntPoint(i);
 
-      Q.Eval(vec, Tr, ip);
       Tr.SetIntPoint (&ip);
+      Q.Eval(vec, Tr, ip);
       vec *= Tr.Weight() * ip.weight;
       el.CalcShape(ip, shape);
       for (int k = 0; k < vdim; k++)
@@ -455,7 +455,7 @@ void VectorFEDomainLFIntegrator::AssembleRHSElementVect(
 {
    int dof = el.GetDof();
    int spaceDim = Tr.GetSpaceDim();
-   int vdim = std::max(spaceDim, el.GetVDim());
+   int vdim = std::max(spaceDim, el.GetRangeDim());
 
    vshape.SetSize(dof,vdim);
    vec.SetSize(vdim);
@@ -656,7 +656,7 @@ void VectorFEBoundaryTangentLFIntegrator::AssembleRHSElementVect(
 {
    int dof = el.GetDof();
    int dim = el.GetDim();
-   int vdim = el.GetVDim();
+   int vdim = el.GetRangeDim();
    DenseMatrix vshape(dof, vdim);
    Vector f_loc(3);
    Vector f_hat(2);

fem/pfespace.cpp

@@ -1946,9 +1946,8 @@ struct PMatrixRow
    void AddRow(const PMatrixRow &other, double coef)
    {
       elems.reserve(elems.size() + other.elems.size());
-      for (unsigned i = 0; i < other.elems.size(); i++)
+      for (const PMatrixElement &oei : other.elems)
       {
-         const PMatrixElement &oei = other.elems[i];
          elems.push_back(
             PMatrixElement(oei.column, oei.stride, coef * oei.value));
       }
@@ -2058,7 +2057,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).LookUp(ri.index);
+      const MeshId &id = *pncmesh->GetNCList(ri.entity).GetMeshIdAndType(ri.index).id;
       ent_ids[ri.entity].Append(id);
       row_idx[ri.entity].Append(i);
       group_ids[ri.entity].Append(ri.group);
@@ -2132,25 +2131,33 @@ void NeighborRowMessage::Decode(int rank)
    rows.clear();
    rows.reserve(nrows);
 
-   // read rows
+   // read rows ent = {0,1,2} means vertex, edge and face entity
    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
-         const int *ind = NULL;
+         // 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;
          if (ent == 1)
          {
+            // edge NC orientation is element defined.
             int eo = pncmesh->GetEdgeNCOrientation(id);
             ind = fec->DofOrderForOrientation(Geometry::SEGMENT, eo);
          }
          else if (ent == 2)
          {
-            Geometry::Type geom = pncmesh->GetFaceGeometry(id.index);
+            geom = pncmesh->GetFaceGeometry(id.index);
             int fo = pncmesh->GetFaceOrientation(id.index);
             ind = fec->DofOrderForOrientation(geom, fo);
          }
@@ -2165,13 +2172,14 @@ 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);
 
@@ -2181,6 +2189,69 @@ 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();
+         }
       }
    }
 }
@@ -2308,12 +2379,6 @@ 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.
-   virtual void UpdateMeshPointer(Mesh *new_mesh);
+   void UpdateMeshPointer(Mesh *new_mesh) override;
 
    /// 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.
-   virtual void CopyProlongationAndRestriction(const FiniteElementSpace &fes,
-                                               const Array<int> *perm);
+   void CopyProlongationAndRestriction(const FiniteElementSpace &fes,
+                                       const Array<int> *perm) override;
 
 public:
    // Face-neighbor data
@@ -282,32 +282,32 @@ public:
    { return Dof_TrueDof_Matrix()->GetGlobalNumCols(); }
 
    /// Return the number of local vector true dofs.
-   virtual int GetTrueVSize() const { return ltdof_size; }
+   int GetTrueVSize() const override { return ltdof_size; }
 
    /// Returns indexes of degrees of freedom in array dofs for i'th element.
-   virtual DofTransformation *GetElementDofs(int i, Array<int> &dofs) const;
+   DofTransformation *GetElementDofs(int i, Array<int> &dofs) const override;
 
    /// Returns indexes of degrees of freedom for i'th boundary element.
-   virtual DofTransformation *GetBdrElementDofs(int i, Array<int> &dofs) const;
+   DofTransformation *GetBdrElementDofs(int i, Array<int> &dofs) const override;
 
    /** Returns the indexes of the degrees of freedom for i'th face
        including the dofs for the edges and the vertices of the face. */
-   virtual int GetFaceDofs(int i, Array<int> &dofs, int variant = 0) const;
+   int GetFaceDofs(int i, Array<int> &dofs, int variant = 0) const override;
 
    /** 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. */
-   virtual const FiniteElement *GetFE(int i) const;
+   const FiniteElement *GetFE(int i) const override;
 
    /** 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 */
-   virtual const FaceRestriction *GetFaceRestriction(
+   const FaceRestriction *GetFaceRestriction(
       ElementDofOrdering f_ordering, FaceType type,
-      L2FaceValues mul = L2FaceValues::DoubleValued) const;
+      L2FaceValues mul = L2FaceValues::DoubleValued) const override;
 
    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
-   virtual void GetEssentialVDofs(const Array<int> &bdr_attr_is_ess,
-                                  Array<int> &ess_dofs,
-                                  int component = -1) const;
+   void GetEssentialVDofs(const Array<int> &bdr_attr_is_ess,
+                          Array<int> &ess_dofs,
+                          int component = -1) const override;
 
    /** Get a list of essential true dofs, ess_tdof_list, corresponding to the
        boundary attributes marked in the array bdr_attr_is_ess. */
-   virtual void GetEssentialTrueDofs(const Array<int> &bdr_attr_is_ess,
-                                     Array<int> &ess_tdof_list,
-                                     int component = -1);
+   void GetEssentialTrueDofs(const Array<int> &bdr_attr_is_ess,
+                             Array<int> &ess_tdof_list,
+                             int component = -1) override;
 
    /** 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;
 
-   virtual const Operator *GetProlongationMatrix() const;
+   const Operator *GetProlongationMatrix() const override;
    /** Get an Operator that performs the action of GetRestrictionMatrix(),
        but potentially with a non-assembled optimized matrix-free
        implementation. */
-   virtual const Operator *GetRestrictionOperator() const;
+   const Operator *GetRestrictionOperator() const override;
    /// Get the R matrix which restricts a local dof vector to true dof vector.
-   virtual const SparseMatrix *GetRestrictionMatrix() const
+   const SparseMatrix *GetRestrictionMatrix() const override
    { 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.
-   virtual void GetTrueTransferOperator(const FiniteElementSpace &coarse_fes,
-                                        OperatorHandle &T) const;
+   void GetTrueTransferOperator(const FiniteElementSpace &coarse_fes,
+                                OperatorHandle &T) const override;
 
    /** Reflect changes in the mesh. Calculate one of the refinement/derefinement
        /rebalance matrices, unless want_transform is false. */
-   virtual void Update(bool want_transform = true);
+   void Update(bool want_transform = true) override;
 
    /// Free ParGridFunction transformation matrix (if any), to save memory.
-   virtual void UpdatesFinished()
+   void UpdatesFinished() override
    {
       FiniteElementSpace::UpdatesFinished();
       old_dof_offsets.DeleteAll();
@@ -442,9 +442,9 @@ public:
 
    const GroupCommunicator &GetGroupCommunicator() const;
 
-   virtual void Mult(const Vector &x, Vector &y) const;
+   void Mult(const Vector &x, Vector &y) const override;
 
-   virtual void MultTranspose(const Vector &x, Vector &y) const;
+   void MultTranspose(const Vector &x, Vector &y) const override;
 };
 
 /// Auxiliary device class used by ParFiniteElementSpace.
@@ -493,9 +493,9 @@ public:
 
    virtual ~DeviceConformingProlongationOperator();
 
-   virtual void Mult(const Vector &x, Vector &y) const;
+   void Mult(const Vector &x, Vector &y) const override;
 
-   virtual void MultTranspose(const Vector &x, Vector &y) const;
+   void MultTranspose(const Vector &x, Vector &y) const override;
 };
 
 }

fem/pgridfunc.cpp

@@ -39,9 +39,11 @@ ParGridFunction::ParGridFunction(ParMesh *pmesh, const GridFunction *gf,
 {
    const FiniteElementSpace *glob_fes = gf->FESpace();
    // duplicate the FiniteElementCollection from 'gf'
-   fec = FiniteElementCollection::New(glob_fes->FEColl()->Name());
+   fec.reset(FiniteElementCollection::New(glob_fes->FEColl()->Name()));
    // create a local ParFiniteElementSpace from the global one:
-   fes = pfes = new ParFiniteElementSpace(pmesh, glob_fes, partitioning, fec);
+   fes = pfes = new ParFiniteElementSpace(pmesh, glob_fes, partitioning,
+                                          fec.get());
+   owned_fes.reset(pfes);
    SetSize(pfes->GetVSize());
 
    if (partitioning)
@@ -81,10 +83,10 @@ ParGridFunction::ParGridFunction(ParMesh *pmesh, std::istream &input)
    : GridFunction(pmesh, input)
 {
    // Convert the FiniteElementSpace, fes, to a ParFiniteElementSpace:
-   pfes = new ParFiniteElementSpace(pmesh, fec, fes->GetVDim(),
+   pfes = new ParFiniteElementSpace(pmesh, fec.get(), fes->GetVDim(),
                                     fes->GetOrdering());
-   delete fes;
    fes = pfes;
+   owned_fes.reset(pfes);
 }
 
 void ParGridFunction::Update()
@@ -469,17 +471,16 @@ void ParGridFunction::GetVectorValue(ElementTransformation &T,
    }
 
    Array<int> vdofs;
-   DofTransformation * doftrans = pfes->GetFaceNbrElementVDofs(nbr_el_no,
-                                                               vdofs);
-   const FiniteElement *fe = pfes->GetFaceNbrFE(nbr_el_no);
-
-   int dof = fe->GetDof();
+   DofTransformation * doftrans = pfes->GetFaceNbrElementVDofs(nbr_el_no, vdofs);
    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);
@@ -501,7 +502,7 @@ void ParGridFunction::GetVectorValue(ElementTransformation &T,
    else
    {
       int spaceDim = pfes->GetMesh()->SpaceDimension();
-      int vdim = std::max(spaceDim, fe->GetVDim());
+      int vdim = std::max(spaceDim, fe->GetRangeDim());
       DenseMatrix vshape(dof, vdim);
       fe->CalcVShape(T, vshape);
       val.SetSize(vdim);

fem/pgridfunc.hpp

@@ -49,10 +49,6 @@ 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.
@@ -93,15 +89,6 @@ 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; }
@@ -212,7 +199,7 @@ public:
    virtual double GetValue(int i, const IntegrationPoint &ip,
                            int vdim = 1) const;
    double GetValue(ElementTransformation &T)
-   { return GetValue(T.ElementNo, T.GetIntPoint()); }
+   { return GetValue(T, T.GetIntPoint()); }
 
    // Redefine to handle the case when T describes a face-neighbor element
    virtual double GetValue(ElementTransformation &T, const IntegrationPoint &ip,
@@ -459,8 +446,6 @@ public:
 
    /// Merge the local grid functions
    void SaveAsOne(std::ostream &out = mfem::out) const;
-
-   virtual ~ParGridFunction() { }
 };
 
 

fem/qspace.cpp

@@ -164,13 +164,20 @@ int FaceQuadratureSpace::GetPermutedIndex(int idx, int iq) const
 
 int FaceQuadratureSpace::GetEntityIndex(const ElementTransformation &T) const
 {
+   auto get_face_index = [this](const int idx)
+   {
+      const auto it = face_indices_inv.find(idx);
+      if (it == face_indices_inv.end()) { return -1; }
+      else { return it->second; }
+   };
+
    switch (T.ElementType)
    {
       case ElementTransformation::FACE:
-         return face_indices_inv.at(T.ElementNo);
+         return get_face_index(T.ElementNo);
       case ElementTransformation::BDR_ELEMENT:
       case ElementTransformation::BDR_FACE:
-         return face_indices_inv.at(mesh.GetBdrElementEdgeIndex(T.ElementNo));
+         return get_face_index(mesh.GetBdrElementEdgeIndex(T.ElementNo));
       default:
          MFEM_ABORT("Invalid element type.");
          return -1;

fem/qspace.hpp

@@ -87,7 +87,8 @@ public:
    ///
    /// For a QuadratureSpace defined on elements, this just returns the element
    /// index. For FaceQuadratureSpace, the returned index depends on the chosen
-   /// FaceType.
+   /// FaceType. If the entity is not found (for example, if @a T represents an
+   /// interior face, and the space has FaceType::Boundary) then -1 is returned.
    virtual int GetEntityIndex(const ElementTransformation &T) const = 0;
 
    /// Write the QuadratureSpace to the stream @a out.

general/kdtree.hpp

@@ -17,6 +17,7 @@
 #include <fstream>
 #include <iostream>
 #include <cmath>
+#include <cstdint>
 #include <tuple>
 
 namespace mfem

general/table.cpp

@@ -218,7 +218,8 @@ 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);
+   MFEM_ASSERT(i >=0 &&
+               i<size, "Index out of bounds.  i = " << i << " size " << size);
 
    for (int k = I[i], end = I[i+1]; k < end; k++)
    {

linalg/ginkgo.cpp

@@ -21,6 +21,7 @@
 #include <iomanip>
 #include <algorithm>
 #include <cmath>
+#include <cstring>
 
 namespace mfem
 {
@@ -28,8 +29,14 @@ 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:
@@ -49,13 +56,23 @@ GinkgoExecutor::GinkgoExecutor(ExecType exec_type)
 #ifdef MFEM_USE_CUDA
             int current_device = 0;
             MFEM_GPU_CHECK(cudaGetDevice(&current_device));
-            executor = gko::CudaExecutor::create(current_device,
-                                                 gko::OmpExecutor::create());
+            if (gko_with_omp_support)
+            {
+               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:
@@ -65,24 +82,111 @@ GinkgoExecutor::GinkgoExecutor(ExecType exec_type)
 #ifdef MFEM_USE_HIP
             int current_device = 0;
             MFEM_GPU_CHECK(hipGetDevice(&current_device));
-            executor = gko::HipExecutor::create(current_device,
-                                                gko::OmpExecutor::create());
+            if (gko_with_omp_support)
+            {
+               executor = gko::HipExecutor::create(current_device,
+                                                   gko::OmpExecutor::create());
+            }
+            else
+            {
+               executor = gko::HipExecutor::create(current_device,
+                                                   gko::ReferenceExecutor::create());
+            }
 #endif
          }
          else
-            mfem::err << "gko::HipExecutor::get_num_devices() did not report "
-                      << "any valid devices" << std::endl;
+         {
+            MFEM_ABORT("gko::HipExecutor::get_num_devices() did not report "
+                       "any valid devices.");
+         }
          break;
       }
       default:
-         mfem::err << "Invalid ExecType specified" << std::endl;
+         MFEM_ABORT("Invalid ExecType specified");
    }
 }
 
+// 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)
 {
-
-   // Pick "best match" Executor based on MFEM device configuration.
+   gko::version_info gko_version = gko::version_info::get();
+   bool gko_with_omp_support = (strcmp(gko_version.omp_version.tag,
+                                       "not compiled") != 0);
    if (mfem_device.Allows(Backend::CUDA_MASK))
    {
       if (gko::CudaExecutor::get_num_devices() > 0)
@@ -90,13 +194,23 @@ GinkgoExecutor::GinkgoExecutor(Device &mfem_device)
 #ifdef MFEM_USE_CUDA
          int current_device = 0;
          MFEM_GPU_CHECK(cudaGetDevice(&current_device));
-         executor = gko::CudaExecutor::create(current_device,
-                                              gko::OmpExecutor::create());
+         if (gko_with_omp_support)
+         {
+            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))
    {
@@ -105,16 +219,123 @@ GinkgoExecutor::GinkgoExecutor(Device &mfem_device)
 #ifdef MFEM_USE_HIP
          int current_device = 0;
          MFEM_GPU_CHECK(hipGetDevice(&current_device));
-         executor = gko::HipExecutor::create(current_device, gko::OmpExecutor::create());
+         if (gko_with_omp_support)
+         {
+            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
    {
-      executor = gko::OmpExecutor::create();
+      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();
+      }
    }
 }
 
@@ -195,7 +416,7 @@ const
                            gko::log::Logger::criterion_check_completed_mask);
 #endif
    residual_logger = std::make_shared<ResidualLogger<>>(executor,
-                                                        gko::lend(system_oper),b);
+                                                        system_oper.get(),b);
 
 }
 
@@ -334,7 +555,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::lend(gko_x));
+   initialize_ginkgo_log(gko_x.get());
 
    MFEM_VERIFY(convergence_logger, "convergence logger not initialized" );
    if (print_level==1)
@@ -350,7 +571,11 @@ 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();
@@ -467,6 +692,10 @@ 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);
 }
@@ -878,7 +1107,11 @@ 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)
@@ -917,6 +1150,10 @@ 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,9 +265,13 @@ 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.
-   b->compute_norm2(lend(b_norm));
+#if MFEM_GINKGO_VERSION < 10600
+   b->compute_norm2(gko::lend(b_norm));
+#else
+   b->compute_norm2(b_norm);
+#endif
    // Use the other utility function to return the norm contained in `b_norm``
-   return std::pow(get_norm(lend(b_norm)),2);
+   return std::pow(get_norm(b_norm.get()),2);
 }
 
 /**
@@ -319,84 +323,43 @@ struct ResidualLogger : gko::log::Logger
 
    using gko_dense = gko::matrix::Dense<ValueType>;
 
-   // 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
-   {
-      // 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);
-   }
-
-   // Version for solver that doesn't log implicit res norm
+   // 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
    {
-      on_iteration_complete(op, iteration, residual, solution, residual_norm,
-                            nullptr);
+      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,
+                              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 MFEM_GINKGO_VERSION > 10500
+   // Ginkgo 1.6 and newer
+   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
+   {
+      iteration_complete_core(iteration, residual, solution, residual_norm,
+                              implicit_sq_residual_norm);
+   }
+#endif
 
    // Construct the logger and store the system matrix and b vectors
    ResidualLogger(std::shared_ptr<const gko::Executor> exec,
@@ -428,6 +391,89 @@ 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
@@ -469,18 +515,40 @@ 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 OmpExecutor; if MFEM is using CUDA,
-    * Ginkgo will choose the CudaExecutor).
+    * 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).
     */
    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

@@ -5269,7 +5269,6 @@ HypreAMS::HypreAMS(const HypreParMatrix &A, HypreParMatrix *G_,
    MFEM_ASSERT(y != NULL, "");
    int sdim = (z == NULL) ? 2 : 3;
    int cycle_type = 13;
-
    MakeSolver(sdim, cycle_type);
 
    HYPRE_ParVector pz = z ? static_cast<HYPRE_ParVector>(*z) : NULL;
@@ -5277,6 +5276,39 @@ HypreAMS::HypreAMS(const HypreParMatrix &A, HypreParMatrix *G_,
    HYPRE_AMSSetDiscreteGradient(ams, *G);
 }
 
+void HypreAMS::Init(ParFiniteElementSpace *edge_fespace)
+{
+   ParMesh *pmesh = edge_fespace->GetParMesh();
+   int dim = pmesh->Dimension();
+   int sdim = pmesh->SpaceDimension();
+   int cycle_type = 13;
+
+   const FiniteElementCollection *edge_fec = edge_fespace->FEColl();
+   bool trace_space = dynamic_cast<const ND_Trace_FECollection *>(edge_fec);
+   bool rt_trace_space = dynamic_cast<const RT_Trace_FECollection *>(edge_fec);
+   trace_space = trace_space || rt_trace_space;
+
+   ND_Trace_FECollection *nd_tr_fec = NULL;
+   if (rt_trace_space)
+   {
+      MFEM_VERIFY(!edge_fespace->IsVariableOrder(),
+                  "HypreAMS does not support variable order spaces");
+      nd_tr_fec = new ND_Trace_FECollection(edge_fec->GetOrder(), dim);
+      edge_fespace = new ParFiniteElementSpace(pmesh, nd_tr_fec);
+   }
+
+   int vdim = edge_fespace->FEColl()->GetRangeDim(dim - trace_space);
+
+   MakeSolver(std::max(sdim, vdim), cycle_type);
+   MakeGradientAndInterpolation(edge_fespace, cycle_type);
+
+   if (rt_trace_space)
+   {
+      delete edge_fespace;
+      delete nd_tr_fec;
+   }
+}
+
 void HypreAMS::MakeSolver(int sdim, int cycle_type)
 {
    int rlx_sweeps       = 1;
@@ -5310,7 +5342,7 @@ void HypreAMS::MakeSolver(int sdim, int cycle_type)
    HYPRE_AMSSetCycleType(ams, cycle_type);
    HYPRE_AMSSetPrintLevel(ams, 1);
 
-   // set additional AMS options
+   // Set additional AMS options
    HYPRE_AMSSetSmoothingOptions(ams, rlx_type, rlx_sweeps, rlx_weight, rlx_omega);
    HYPRE_AMSSetAlphaAMGOptions(ams, amg_coarsen_type, amg_agg_levels, amg_rlx_type,
                                theta, amg_interp_type, amg_Pmax);
@@ -5330,27 +5362,20 @@ void HypreAMS::MakeSolver(int sdim, int cycle_type)
 void HypreAMS::MakeGradientAndInterpolation(
    ParFiniteElementSpace *edge_fespace, int cycle_type)
 {
-   int dim = edge_fespace->GetMesh()->Dimension();
-   int sdim = edge_fespace->GetMesh()->SpaceDimension();
    const FiniteElementCollection *edge_fec = edge_fespace->FEColl();
-
-   bool trace_space, rt_trace_space;
-   ND_Trace_FECollection *nd_tr_fec = NULL;
-   trace_space = dynamic_cast<const ND_Trace_FECollection*>(edge_fec);
-   rt_trace_space = dynamic_cast<const RT_Trace_FECollection*>(edge_fec);
-   trace_space = trace_space || rt_trace_space;
-
-   MFEM_VERIFY(!edge_fespace->IsVariableOrder(), "");
-   int p = edge_fec->GetOrder();
+   bool trace_space = dynamic_cast<const ND_Trace_FECollection *>(edge_fec);
 
    ParMesh *pmesh = edge_fespace->GetParMesh();
-   if (rt_trace_space)
-   {
-      nd_tr_fec = new ND_Trace_FECollection(p, dim);
-      edge_fespace = new ParFiniteElementSpace(pmesh, nd_tr_fec);
-   }
+   int dim = pmesh->Dimension();
+   int sdim = pmesh->SpaceDimension();
+   int vdim = edge_fespace->FEColl()->GetRangeDim(dim - trace_space);
 
-   // define the nodal linear finite element space associated with edge_fespace
+   // For dim = 1, ND_FECollection::GetOrder() returns p - 1
+   MFEM_VERIFY(!edge_fespace->IsVariableOrder(),
+               "HypreAMS does not support variable order spaces");
+   int p = edge_fec->GetOrder() + (dim - trace_space == 1 ? 1 : 0);
+
+   // Define the nodal linear finite element space associated with edge_fespace
    FiniteElementCollection *vert_fec;
    if (trace_space)
    {
@@ -5363,8 +5388,27 @@ void HypreAMS::MakeGradientAndInterpolation(
    ParFiniteElementSpace *vert_fespace = new ParFiniteElementSpace(pmesh,
                                                                    vert_fec);
 
-   // generate and set the vertex coordinates
-   if (p == 1 && pmesh->GetNodes() == NULL)
+   // generate and set the discrete gradient
+   ParDiscreteLinearOperator *grad;
+   grad = new ParDiscreteLinearOperator(vert_fespace, edge_fespace);
+   if (trace_space)
+   {
+      grad->AddTraceFaceInterpolator(new GradientInterpolator);
+   }
+   else
+   {
+      grad->AddDomainInterpolator(new GradientInterpolator);
+   }
+   grad->Assemble();
+   grad->Finalize();
+   G = grad->ParallelAssemble();
+   HYPRE_AMSSetDiscreteGradient(ams, *G);
+   delete grad;
+
+   // generate and set the vertex coordinates or Nedelec interpolation matrices
+   x = y = z = NULL;
+   Pi = Pix = Piy = Piz = NULL;
+   if (p == 1 && pmesh->GetNodes() == NULL && vdim <= sdim)
    {
       ParGridFunction x_coord(vert_fespace);
       ParGridFunction y_coord(vert_fespace);
@@ -5400,34 +5444,9 @@ void HypreAMS::MakeGradientAndInterpolation(
    }
    else
    {
-      x = NULL;
-      y = NULL;
-      z = NULL;
-   }
-
-   // generate and set the discrete gradient
-   ParDiscreteLinearOperator *grad;
-   grad = new ParDiscreteLinearOperator(vert_fespace, edge_fespace);
-   if (trace_space)
-   {
-      grad->AddTraceFaceInterpolator(new GradientInterpolator);
-   }
-   else
-   {
-      grad->AddDomainInterpolator(new GradientInterpolator);
-   }
-   grad->Assemble();
-   grad->Finalize();
-   G = grad->ParallelAssemble();
-   HYPRE_AMSSetDiscreteGradient(ams, *G);
-   delete grad;
-
-   // generate and set the Nedelec interpolation matrices
-   Pi = Pix = Piy = Piz = NULL;
-   if (p > 1 || pmesh->GetNodes() != NULL)
-   {
-      ParFiniteElementSpace *vert_fespace_d
-         = new ParFiniteElementSpace(pmesh, vert_fec, sdim, Ordering::byVDIM);
+      ParFiniteElementSpace *vert_fespace_d =
+         new ParFiniteElementSpace(pmesh, vert_fec, std::max(sdim, vdim),
+                                   Ordering::byVDIM);
 
       ParDiscreteLinearOperator *id_ND;
       id_ND = new ParDiscreteLinearOperator(vert_fespace_d, edge_fespace);
@@ -5451,8 +5470,8 @@ void HypreAMS::MakeGradientAndInterpolation(
          Array2D<HypreParMatrix *> Pi_blocks;
          id_ND->GetParBlocks(Pi_blocks);
          Pix = Pi_blocks(0,0);
-         if (sdim >= 2) { Piy = Pi_blocks(0,1); }
-         if (sdim == 3) { Piz = Pi_blocks(0,2); }
+         if (std::max(sdim, vdim) >= 2) { Piy = Pi_blocks(0,1); }
+         if (std::max(sdim, vdim) == 3) { Piz = Pi_blocks(0,2); }
       }
 
       delete id_ND;
@@ -5468,20 +5487,6 @@ void HypreAMS::MakeGradientAndInterpolation(
 
    delete vert_fespace;
    delete vert_fec;
-
-   if (rt_trace_space)
-   {
-      delete edge_fespace;
-      delete nd_tr_fec;
-   }
-}
-
-void HypreAMS::Init(ParFiniteElementSpace *edge_fespace)
-{
-   int cycle_type = 13;
-   int sdim = edge_fespace->GetMesh()->SpaceDimension();
-   MakeSolver(sdim, cycle_type);
-   MakeGradientAndInterpolation(edge_fespace, cycle_type);
 }
 
 void HypreAMS::ResetAMSPrecond()

linalg/hypre.hpp

@@ -1755,6 +1755,10 @@ private:
    void MakeGradientAndInterpolation(ParFiniteElementSpace *edge_fespace,
                                      int cycle_type);
 
+   // Recreates another AMS solver with the same options when SetOperator is
+   // called multiple times.
+   void ResetAMSPrecond();
+
    /// The underlying hypre solver object
    HYPRE_Solver ams;
    /// Vertex coordinates
@@ -1773,10 +1777,6 @@ private:
    /// Flag set if `SetPrintLevel` is called, needed in `ResetAMSPrecond`
    int print_level = 1;
 
-   // Recreates another AMS solver with the same options when SetOperator is
-   // called multiple times.
-   void ResetAMSPrecond();
-
 public:
    /// @brief Construct the AMS solver on the given edge finite element space.
    ///

linalg/sundials.cpp

@@ -1672,7 +1672,7 @@ void ARKStepSolver::UseMFEMMassLinearSolver(int tdep)
    MFEM_VERIFY(M, "error in SUNMatNewEmpty()");
 
    M->content      = this;
-   M->ops->getid   = SUNMatGetID;
+   M->ops->getid   = MatGetID;
    M->ops->matvec  = ARKStepSolver::MassMult1;
    M->ops->destroy = MatDestroy;
 

mesh/mesh.cpp

@@ -6271,6 +6271,111 @@ 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);
@@ -6564,24 +6669,20 @@ Array<int> Mesh::FindFaceNeighbors(const int elem) const
 
 void Mesh::GetBdrElementFace(int i, int *f, int *o) const
 {
-   const int *bv, *fv;
+   *f = GetBdrElementEdgeIndex(i);
 
-   *f = be_to_face[i];
-   bv = boundary[i]->GetVertices();
-   fv = faces[be_to_face[i]]->GetVertices();
+   const int *fv = (Dim > 1) ? faces[*f]->GetVertices() : NULL;
+   const int *bv = boundary[i]->GetVertices();
 
    // find the orientation of the bdr. elem. w.r.t.
    // the corresponding face element (that's the base)
-   switch (GetBdrElementType(i))
+   switch (GetBdrElementGeometry(i))
    {
-      case Element::TRIANGLE:
-         *o = GetTriOrientation(fv, bv);
-         break;
-      case Element::QUADRILATERAL:
-         *o = GetQuadOrientation(fv, bv);
-         break;
-      default:
-         MFEM_ABORT("invalid geometry");
+      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");
    }
 }
 
@@ -7183,7 +7284,7 @@ STable3D *Mesh::GetFacesTable()
             break;
          }
          default:
-            MFEM_ABORT("Unexpected type of Element.");
+            MFEM_ABORT("Unexpected type of Element: " << GetElementType(i));
       }
    }
    return faces_tbl;
@@ -9558,7 +9659,8 @@ void Mesh::NonconformingRefinement(const Array<Refinement> &refinements,
 double Mesh::AggregateError(const Array<double> &elem_error,
                             const int *fine, int nfine, int op)
 {
-   double error = elem_error[fine[0]];
+   double error = (op == 3) ? std::pow(elem_error[fine[0]],
+                                       2.0) : elem_error[fine[0]];
 
    for (int i = 1; i < nfine; i++)
    {
@@ -9570,9 +9672,11 @@ double Mesh::AggregateError(const Array<double> &elem_error,
          case 0: error = std::min(error, err_fine); break;
          case 1: error += err_fine; break;
          case 2: error = std::max(error, err_fine); break;
+         case 3: error += std::pow(err_fine, 2.0); break;
+         default: MFEM_ABORT("Invalid operation.");
       }
    }
-   return error;
+   return (op == 3) ? std::sqrt(error) : error;
 }
 
 bool Mesh::NonconformingDerefinement(Array<double> &elem_error,

mesh/mesh.hpp

@@ -322,9 +322,38 @@ 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().
@@ -1330,7 +1359,12 @@ public:
        element @a elem, including @a elem. */
    Array<int> FindFaceNeighbors(const int elem) const;
 
-   /// Return the index and the orientation of the face of bdr element i. (3D)
+   /** 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. */
    void GetBdrElementFace(int i, int *f, int *o) const;
 
    /** Return the vertex index of boundary element i. (1D)
@@ -1409,6 +1443,14 @@ 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'
-      int type;
-      const MeshId &eid = edge_list.LookUp(nd.edge_index, &type);
-      if (type == 1)
+      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)
       {
          // in this case we need to add an edge-face constraint, because the
-         // master edge is really a (face-)slave itself
-
+         // non-slave edge is really a (face-)slave itself.
+         const MeshId &eid = *eid_and_type.id;
          face_list.slaves.Append(
             Slave(-1 - eid.index, eid.element, eid.local, Geometry::TRIANGLE));
 
@@ -3048,9 +3048,10 @@ void NCMesh::TraverseTetEdge(int vn0, int vn1, const Point &p0, const Point &p1,
    TraverseTetEdge(mid, vn1, pmid, p1, matrix_map);
 }
 
-bool NCMesh::TraverseTriFace(int vn0, int vn1, int vn2,
-                             const PointMatrix& pm, int level,
-                             MatrixMap &matrix_map)
+NCMesh::TriFaceTraverseResults NCMesh::TraverseTriFace(int vn0, int vn1,
+                                                       int vn2,
+                                                       const PointMatrix& pm, int level,
+                                                       MatrixMap &matrix_map)
 {
    if (level > 0)
    {
@@ -3069,7 +3070,7 @@ bool NCMesh::TraverseTriFace(int vn0, int vn1, int vn2,
          sl.local = ReorderFacePointMat(vn0, vn1, vn2, -1, elem, pm, pm_r);
          sl.matrix = matrix_map.GetIndex(pm_r);
 
-         return true;
+         return {true, elements[elem].rank != MyRank};
       }
    }
 
@@ -3077,7 +3078,7 @@ bool NCMesh::TraverseTriFace(int vn0, int vn1, int vn2,
    if (TriFaceSplit(vn0, vn1, vn2, mid))
    {
       Point pmid0(pm(0), pm(1)), pmid1(pm(1), pm(2)), pmid2(pm(2), pm(0));
-      bool b[4];
+      TriFaceTraverseResults b[4];
 
       b[0] = TraverseTriFace(vn0, mid[0], mid[2],
                              PointMatrix(pm(0), pmid0, pmid2),
@@ -3095,16 +3096,21 @@ bool NCMesh::TraverseTriFace(int vn0, int vn1, int vn2,
                              PointMatrix(pmid1, pmid2, pmid0),
                              level+1, matrix_map);
 
-      // traverse possible tet edges constrained by the master face
-      if (HaveTets() && !b[3])
+      // 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))
       {
-         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); }
+         // 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); }
       }
    }
-
-   return false;
+   return {false, false};
 }
 
 void NCMesh::BuildFaceList()
@@ -3402,77 +3408,80 @@ void NCMesh::NCList::Clear()
       point_matrices[i].DeleteAll();
    }
 
-   inv_index.DeleteAll();
+   inv_index.clear();
 }
 
-long NCMesh::NCList::TotalSize() const
+NCMesh::NCList::MeshIdAndType
+NCMesh::NCList::GetMeshIdAndType(int index) const
 {
-   return conforming.Size() + masters.Size() + slaves.Size();
-}
-
-const NCMesh::MeshId& NCMesh::NCList::LookUp(int index, int *type) const
-{
-   if (!inv_index.Size())
+   BuildIndex();
+   const auto it = inv_index.find(index);
+   auto ft = it != inv_index.end() ? it->second.first : MeshIdType::UNRECOGNIZED;
+   switch (ft)
    {
-      int max_index = -1;
+      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}));
       for (int i = 0; i < conforming.Size(); i++)
       {
-         max_index = std::max(conforming[i].index, max_index);
+         inv_index.emplace(conforming[i].index, std::make_pair(MeshIdType::CONFORMING,
+                                                               i));
       }
       for (int i = 0; i < masters.Size(); i++)
       {
-         max_index = std::max(masters[i].index, max_index);
+         inv_index.emplace(masters[i].index, std::make_pair(MeshIdType::MASTER, i));
       }
       for (int i = 0; i < slaves.Size(); 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;
+         inv_index.emplace(slaves[i].index, std::make_pair(MeshIdType::SLAVE, i));
       }
    }
 
-   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];
-   }
+   MFEM_ASSERT(inv_index.size() > 0,
+               "Empty inverse index, member lists must be populated before BuildIndex is called!");
 }
 
-
 //// Neighbors /////////////////////////////////////////////////////////////////
 
 void NCMesh::CollectEdgeVertices(int v0, int v1, Array<int> &indices)

mesh/ncmesh.hpp

@@ -24,6 +24,7 @@
 #include <vector>
 #include <map>
 #include <iostream>
+#include <unordered_map>
 
 namespace mfem
 {
@@ -223,30 +224,73 @@ public:
          , master(-1), matrix(0), edge_flags(0) {}
    };
 
+
    /// Lists all edges/faces in the nonconforming mesh.
    struct NCList
    {
-      Array<MeshId> conforming;
-      Array<Master> masters;
-      Array<Slave> slaves;
+      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
 
       /// 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();
-      bool Empty() const { return !conforming.Size() && !masters.Size(); }
-      long TotalSize() const;
+      /// 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.
       long MemoryUsage() const;
-
-      const MeshId& LookUp(int index, int *type = NULL) const;
-
       ~NCList() { Clear(); }
    private:
-      mutable Array<int> inv_index;
+      // 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;
    };
 
    /// Return the current list of conforming and nonconforming faces.
@@ -727,9 +771,14 @@ protected: // implementation
    void TraverseQuadFace(int vn0, int vn1, int vn2, int vn3,
                          const PointMatrix& pm, int level, Face* eface[4],
                          MatrixMap &matrix_map);
-   bool TraverseTriFace(int vn0, int vn1, int vn2,
-                        const PointMatrix& pm, int level,
-                        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);
    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,6 +39,7 @@ KnotVector::KnotVector(int Order_, int NCP)
    Order = Order_;
    NumOfControlPoints = NCP;
    knot.SetSize(NumOfControlPoints + Order + 1);
+   NumOfElements = 0;
 
    knot = -1.;
 }
@@ -129,27 +130,35 @@ 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);
 
-   for (int i = 0; i <GetNE() ; i++)
+   /* @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++)
    {
+      // Avoid printing shapes between repeated knots
+      if (!isElement(cnt)) { e--; continue; }
+
       for (int j = 0; j <samples; j++)
       {
          x =j*dx;
-         os<< x + i;
+         os<< x + e;
 
-         CalcShape ( shape, i, x);
+         CalcShape ( shape, cnt, x);
          for (int d = 0; d < Order+1; d++) { os<<"\t"<<shape[d]; }
 
-         CalcDShape ( shape, i, x);
+         CalcDShape ( shape, cnt, x);
          for (int d = 0; d < Order+1; d++) { os<<"\t"<<shape[d]; }
 
-         CalcD2Shape ( shape, i, x);
+         CalcD2Shape ( shape, cnt, x);
          for (int d = 0; d < Order+1; d++) { os<<"\t"<<shape[d]; }
          os<<endl;
       }

mesh/nurbs.hpp

@@ -88,6 +88,10 @@ 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,7 +34,6 @@ 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)
@@ -106,8 +105,7 @@ ParMesh& ParMesh::operator=(ParMesh &&mesh)
 
 ParMesh::ParMesh(MPI_Comm comm, Mesh &mesh, int *partitioning_,
                  int part_method)
-   : face_nbr_el_to_face(NULL)
-   , glob_elem_offset(-1)
+   : glob_elem_offset(-1)
    , glob_offset_sequence(-1)
    , gtopo(comm)
 {
@@ -854,7 +852,6 @@ 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)
@@ -906,7 +903,7 @@ void ParMesh::FinalizeParTopo()
    sface_lface.SetSize(nst + shared_quads.Size());
    if (sface_lface.Size())
    {
-      STable3D *faces_tbl = GetFacesTable();
+      auto faces_tbl = std::unique_ptr<STable3D>(GetFacesTable());
       for (int st = 0; st < nst; st++)
       {
          const int *v = shared_trias[st].v;
@@ -917,13 +914,11 @@ 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)
-   : face_nbr_el_to_face(NULL)
-   , glob_elem_offset(-1)
+   : glob_elem_offset(-1)
    , glob_offset_sequence(-1)
    , gtopo(comm)
 {
@@ -1134,7 +1129,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 = NULL;
+   face_nbr_el_to_face = nullptr;
    glob_elem_offset = -1;
    glob_offset_sequence = -1;
    gtopo = orig_mesh.gtopo;
@@ -2123,7 +2118,7 @@ void ParMesh::ExchangeFaceNbrData()
 
    if (Nonconforming())
    {
-      // with ParNCMesh we can set up face neighbors without communication
+      // with ParNCMesh we can set up face neighbors mostly without communication
       pncmesh->GetFaceNeighbors(*this);
       have_face_nbr_data = true;
 
@@ -2192,7 +2187,7 @@ void ParMesh::ExchangeFaceNbrData()
 
    if (Dim == 3)
    {
-      GetFaceNbrElementToFaceTable();
+      BuildFaceNbrElementToFaceTable();
    }
 
    if (del_tables) { delete gr_sface; }
@@ -2457,8 +2452,7 @@ 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.Clear();
-   face_nbr_el_ori.SetSize(face_nbr_elements_offset[num_face_nbrs], 6);
+   face_nbr_el_ori.reset(new Table(face_nbr_elements_offset[num_face_nbrs], 6));
    while (true)
    {
       int fn;
@@ -2489,7 +2483,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];
@@ -2499,7 +2493,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);
 
@@ -2721,191 +2715,118 @@ STable3D *ParMesh::GetSharedFacesTable()
    return sfaces_tbl;
 }
 
-STable3D *ParMesh::GetFaceNbrElementToFaceTable(int ret_ftbl)
+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])
 {
-   int i, *v;
-   STable3D * faces_tbl = GetFacesTable();
-   STable3D * sfaces_tbl = GetSharedFacesTable();
+   for (int i = start; i < end; ++i)
+   {
+      // 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]]);
 
-   if (face_nbr_el_to_face != NULL)
-   {
-      delete face_nbr_el_to_face;
+      // 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);
    }
-   face_nbr_el_to_face = new Table(face_nbr_elements.Size(), 6);
-   for (i = 0; i < face_nbr_elements.Size(); i++)
+}
+
+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])
    {
-      v = face_nbr_elements[i]->GetVertices();
+      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);
       switch (face_nbr_elements[i]->GetType())
       {
          case Element::TETRAHEDRON:
          {
-            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);
-            }
+            AddTriFaces(v, faces, shared_faces, i, 0, 4, tet_t::FaceVert);
             break;
          }
          case Element::WEDGE:
          {
-            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);
-            }
+            AddTriFaces(v, faces, shared_faces, i, 0, 2, pri_t::FaceVert);
+            add_quad_faces(i, 2, 5, pri_t::FaceVert);
             break;
          }
          case Element::PYRAMID:
          {
-            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);
-            }
+            add_quad_faces(i, 0, 1, pyr_t::FaceVert);
+            AddTriFaces(v, faces, shared_faces, i, 1, 5, pyr_t::FaceVert);
             break;
          }
          case Element::HEXAHEDRON:
          {
-            // 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);
-            }
+            add_quad_faces(i, 0, 6, hex_t::FaceVert);
             break;
          }
          default:
@@ -2913,14 +2834,6 @@ STable3D *ParMesh::GetFaceNbrElementToFaceTable(int ret_ftbl)
       }
    }
    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
@@ -2941,33 +2854,28 @@ int ParMesh::GetFaceNbrRank(int fn) const
 }
 
 void
-ParMesh::GetFaceNbrElementFaces(int i, Array<int> &fcs, Array<int> &cor) const
+ParMesh::GetFaceNbrElementFaces(int i, Array<int> &faces,
+                                Array<int> &orientations) const
 {
-   int n, j;
    int el_nbr = i - GetNE();
-   if (face_nbr_el_to_face)
+   if (face_nbr_el_to_face != nullptr && el_nbr < face_nbr_el_to_face->Size())
    {
-      face_nbr_el_to_face->GetRow(el_nbr, fcs);
+      face_nbr_el_to_face->GetRow(el_nbr, faces);
    }
    else
    {
       MFEM_ABORT("ParMesh::GetFaceNbrElementFaces(...) : "
-                 "face_nbr_el_to_face not generated.");
+                 "face_nbr_el_to_face not generated correctly.");
    }
-   if (el_nbr < face_nbr_el_ori.Size())
+
+   if (face_nbr_el_ori != nullptr && el_nbr < face_nbr_el_ori->Size())
    {
-      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];
-      }
+      face_nbr_el_ori->GetRow(el_nbr, orientations);
    }
    else
    {
       MFEM_ABORT("ParMesh::GetFaceNbrElementFaces(...) : "
-                 "face_nbr_el_to_face not generated.");
+                 "face_nbr_el_ori not generated correctly.");
    }
 }
 
@@ -3227,7 +3135,7 @@ int ParMesh::GetSharedFace(int sface) const
    {
       MFEM_ASSERT(Dim > 1, "");
       const NCMesh::NCList &shared = pncmesh->GetSharedList(Dim-1);
-      int csize = (int) shared.conforming.Size();
+      int csize = shared.conforming.Size();
       return sface < csize
              ? shared.conforming[sface].index
              : shared.slaves[sface - csize].index;
@@ -4579,7 +4487,7 @@ void ParMesh::UniformRefinement3D()
 
    DSTable v_to_v(NumOfVertices);
    GetVertexToVertexTable(v_to_v);
-   STable3D *faces_tbl = GetFacesTable();
+   auto faces_tbl = std::unique_ptr<STable3D>(GetFacesTable());
 
    // call Mesh::UniformRefinement3D_base so that it won't update the nodes
    Array<int> f2qf;
@@ -4594,7 +4502,6 @@ void ParMesh::UniformRefinement3D()
    // update the groups
    UniformRefineGroups3D(old_nv, old_nedges, v_to_v, *faces_tbl,
                          f2qf.Size() ? &f2qf : NULL);
-   delete faces_tbl;
 
    UpdateNodes();
 }
@@ -6693,6 +6600,8 @@ 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);
@@ -6713,8 +6622,7 @@ void ParMesh::Destroy()
    }
    shared_edges.DeleteAll();
 
-   delete face_nbr_el_to_face;
-   face_nbr_el_to_face = NULL;
+   face_nbr_el_to_face = nullptr;
 }
 
 ParMesh::~ParMesh()

mesh/pmesh.hpp

@@ -77,8 +77,11 @@ protected:
    // sface ids: all triangles first, then all quads
    Array<int> sface_lface;
 
-   Table *face_nbr_el_to_face;
-   Table  face_nbr_el_ori; // orientations for each face (from nbr processor)
+   /// 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;
 
    IsoparametricTransformation FaceNbrTransformation;
 
@@ -113,7 +116,32 @@ protected:
    bool DecodeFaceSplittings(HashTable<Hashed2> &v_to_v, const int *v,
                              const Array<unsigned> &codes, int &pos);
 
-   STable3D *GetFaceNbrElementToFaceTable(int ret_ftbl = 0);
+   // 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]);
 
    void GetFaceNbrElementTransformation(
       int i, IsoparametricTransformation *ElTr);
@@ -287,7 +315,7 @@ protected:
 
 public:
    /// Default constructor. Create an empty @a ParMesh.
-   ParMesh() : MyComm(0), NRanks(0), MyRank(-1), face_nbr_el_to_face(NULL),
+   ParMesh() : MyComm(0), NRanks(0), MyRank(-1),
       glob_elem_offset(-1), glob_offset_sequence(-1),
       have_face_nbr_data(false), pncmesh(NULL) { }
 
@@ -465,7 +493,8 @@ public:
    int GetFaceNbrRank(int fn) const;
 
    /** Similar to Mesh::GetElementFaces */
-   void GetFaceNbrElementFaces(int i, Array<int> &fcs, Array<int> &cor) const;
+   void GetFaceNbrElementFaces(int i, Array<int> &faces,
+                               Array<int> &orientation) const;
 
    /** Similar to Mesh::GetFaceToElementTable with added face-neighbor elements
        with indices offset by the local number of elements. */

mesh/pncmesh.cpp

@@ -16,10 +16,12 @@
 #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
 {
@@ -887,6 +889,7 @@ 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.
@@ -901,6 +904,11 @@ 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++)
    {
@@ -914,8 +922,12 @@ 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++)
@@ -924,7 +936,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) { continue; }
+         if (sf.element < 0 || sf.index < 0) { continue; }
 
          MFEM_ASSERT(mf.element >= 0, "");
          Element* e[2] = { &elements[mf.element], &elements[sf.element] };
@@ -938,8 +950,12 @@ 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);
       }
    }
 
@@ -1022,6 +1038,13 @@ 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];
@@ -1175,6 +1198,117 @@ 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,9 +108,7 @@ public:
        passed. */
    void Rebalance(const Array<int> *custom_partition = NULL);
 
-
    // interface for ParFiniteElementSpace
-
    int GetNElements() const { return NElements; }
 
    int GetNGhostVertices() const { return NGhostVertices; }
@@ -141,8 +139,8 @@ public:
       return (index < NFaces) ? face_orient[index] : 0;
    }
 
-   typedef short GroupId;
-   typedef std::vector<int> CommGroup;
+   using GroupId = short;
+   using CommGroup = std::vector<int>;
 
    /// Return vertex/edge/face ('entity' == 0/1/2, resp.) owner.
    GroupId GetEntityOwnerId(int entity, int index)

miniapps/gslib/findpts.cpp

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

miniapps/nurbs/CMakeLists.txt

@@ -17,6 +17,10 @@ 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)
@@ -67,6 +71,9 @@ 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
+SEQ_MINIAPPS = nurbs_ex1 nurbs_patch_ex1 nurbs_curveint nurbs_printfunc
 PAR_MINIAPPS = nurbs_ex1p nurbs_ex11p
 ifeq ($(MFEM_USE_MPI),NO)
    MINIAPPS = $(SEQ_MINIAPPS)
@@ -109,6 +109,9 @@ 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

@@ -0,0 +1,52 @@
+//                          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

@@ -0,0 +1,372 @@
+// 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,332 +265,248 @@ TEST_CASE("H1 PA Coefficient", "[PartialAssembly][Coefficient]")
 TEST_CASE("Hcurl/Hdiv PA Coefficient",
           "[CUDA][PartialAssembly][Coefficient]")
 {
-   for (dimension = 2; dimension < 4; ++dimension)
+   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)
    {
-      const int ne = 3;
-      Mesh mesh = MakeCartesianNonaligned(dimension, ne);
+      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));
+   }
 
-      for (int coeffType = 0; coeffType < 5; ++coeffType)
+   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)
       {
-         Coefficient* coeff = nullptr;
-         Coefficient* coeff2 = nullptr;
-         VectorCoefficient* vcoeff = nullptr;
-         MatrixCoefficient* mcoeff = nullptr;
-         if (coeffType == 0)
-         {
-            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);
-         }
+         fec_test.reset(new RT_FECollection(order - 1, dimension));
+      }
+      else
+      {
+         fec_test.reset(new ND_FECollection(order, dimension));
+      }
 
-         enum MixedSpaces {Hcurl, Hdiv, HcurlHdiv, HdivHcurl, NumSpaceTypes};
+      FiniteElementSpace fes_test(&mesh, fec_test.get());
 
-         for (int spaceType = 0; spaceType < NumSpaceTypes; ++spaceType)
+      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)
          {
-            const int numIntegrators =
-               (spaceType >= HcurlHdiv) ? 1 : ((coeffType == 2) ? 2 : 3);
-
-            for (int integrator = 0; integrator < numIntegrators; ++integrator)
+            if (space_type == HcurlHdiv)
             {
-               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;
-               }
+               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));
             }
          }
-
-         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();
 
-   FiniteElementCollection *fec;
+   std::unique_ptr<FiniteElementCollection> fec;
    if (prob)
    {
       auto basis = prob==3 ? BasisType::Positive : BasisType::GaussLobatto;
-      fec = new L2_FECollection(order, dim, basis);
+      fec.reset(new L2_FECollection(order, dim, basis));
    }
    else
    {
-      fec = new H1_FECollection(order, dim);
+      fec.reset(new H1_FECollection(order, dim));
    }
-   FiniteElementSpace fespace(&mesh, fec);
+   FiniteElementSpace fespace(&mesh, fec.get());
 
    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);
 
-   VectorCoefficient *vel_coeff;
-   Coefficient *rho;
+   std::unique_ptr<VectorCoefficient> vel_coeff;
+   std::unique_ptr<Coefficient> rho;
 
    // prob: 0: CG, 1: DG continuous coeff, 2: DG discontinuous coeff
    if (prob >= 2)
    {
       vel_gf.Randomize(1);
-      vel_coeff = new VectorGridFunctionCoefficient(&vel_gf);
+      vel_coeff.reset(new VectorGridFunctionCoefficient(&vel_gf));
       rho_gf.Randomize(1);
-      rho = new GridFunctionCoefficient(&rho_gf);
+      rho.reset(new GridFunctionCoefficient(&rho_gf));
    }
    else
    {
-      vel_coeff = new VectorFunctionCoefficient(dim, velocity_function);
-      rho = new ConstantCoefficient(1.0);
+      vel_coeff.reset(new VectorFunctionCoefficient(dim, velocity_function));
+      rho.reset(new ConstantCoefficient(1.0));
    }
 
 
@@ -451,10 +451,6 @@ 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
@@ -468,7 +464,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"
-using namespace mfem;
-
 #include "unit_tests.hpp"
 
+using namespace mfem;
+
 class TestMesh : public Mesh
 {
 public:
@@ -117,3 +117,97 @@ 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"
-using namespace mfem;
-
 #include "unit_tests.hpp"
 
+using namespace mfem;
+
 #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"
-using namespace mfem;
-
 #include "unit_tests.hpp"
 
+using namespace mfem;
+
 TEST_CASE("Element-wise construction", "[Mesh]")
 {
    SECTION("Quadrilateral")

tests/unit/mesh/test_ncmesh.cpp

@@ -12,6 +12,7 @@
 #include "mfem.hpp"
 #include "unit_tests.hpp"
 
+#include <array>
 namespace mfem
 {
 
@@ -299,8 +300,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.
-void CheckL2Projection(ParMesh& pmesh, Mesh& smesh, int order,
-                       std::function<double(Vector const&)> exact_soln)
+std::array<double, 2> CheckL2Projection(ParMesh& pmesh, Mesh& smesh, int order,
+                                        std::function<double(Vector const&)> exact_soln)
 {
    REQUIRE(pmesh.GetGlobalNE() == smesh.GetNE());
    REQUIRE(pmesh.Dimension() == smesh.Dimension());
@@ -368,7 +369,6 @@ void 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,26 +393,11 @@ void CheckL2Projection(ParMesh& pmesh, Mesh& smesh, int order,
       return x.ComputeL2Error(rhs_coef) / pnorm;
    }();
 
-   constexpr double test_tol = 1e-9;
-   CHECK(std::abs(serror - perror) < test_tol);
+   return {serror, perror};
 };
 
-
-TEST_CASE("FaceEdgeConstraint",  "[Parallel], [NCMesh]")
+TEST_CASE("EdgeFaceConstraint",  "[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.
@@ -422,103 +407,239 @@ TEST_CASE("FaceEdgeConstraint",  "[Parallel], [NCMesh]")
       return std::sin(d * d);
    };
 
-   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)
+   SECTION("ReferenceTet")
    {
-      refines.Append(0);
-   }
+      constexpr int refining_rank = 0;
+      auto smesh = Mesh("../../data/ref-tetrahedron.mesh");
 
-   // 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);
+      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;
 
-   // Rank 0 refines the parallel mesh, all ranks refine the serial mesh
-   smesh.GeneralRefinement(serial_refines, 1); // nonconformal
-   pmesh.GeneralRefinement(refines, 1); // nonconformal
+         smesh.GeneralRefinement(single_edge_refine, 0); // conformal
+      }
 
-   REQUIRE(pmesh.GetGlobalNE() == 8 + 1);
-   REQUIRE(smesh.GetNE() == 8 + 1);
 
-   // 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}};
+      REQUIRE(smesh.GetNE() == 2);
+      smesh.EnsureNCMesh(true);
+      smesh.Finalize();
 
-   // 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;
+      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)
       {
-         refines[0] = i;
+         refines.Append(0);
       }
-      // 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);
+      // 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);
 
-      REQUIRE(tmp.GetGlobalNE() == 1 + 8 - 1 + 8); // 16 elements
+      // Rank 0 refines the parallel mesh, all ranks refine the serial mesh
+      smesh.GeneralRefinement(serial_refines, 1); // nonconformal
+      pmesh.GeneralRefinement(refines, 1); // nonconformal
 
-      Mesh stmp(smesh);
-      stmp.GeneralRefinement(serial_refines);
-      REQUIRE(stmp.GetNE() == 1 + 8 - 1 + 8); // 16 elements
+      REQUIRE(pmesh.GetGlobalNE() == 8 + 1);
+      REQUIRE(smesh.GetNE() == 8 + 1);
 
-      if (Mpi::WorldRank() == refining_rank)
+      // 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)
       {
-         refines[0] = j;
+         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);
+         }
       }
-      // 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)
@@ -716,7 +837,8 @@ TEST_CASE("P2Q1PureTetHexPri",  "[Parallel], [NCMesh]")
       auto pmesh = ParMesh(MPI_COMM_WORLD, smesh);
 
       // P2 ensures there are triangles without dofs
-      CheckL2Projection(pmesh, smesh, 2, exact_soln);
+      auto error = CheckL2Projection(pmesh, smesh, 2, exact_soln);
+      CHECK(std::abs(error[1] - error[0]) < 1e-9);
    }
 } // test case
 
@@ -759,11 +881,153 @@ TEST_CASE("PNQ2PureTetHexPri",  "[Parallel], [NCMesh]")
 
       for (int p = 1; p < 3; ++p)
       {
-         CheckL2Projection(pmesh, smesh, p, exact_soln);
+         auto error = CheckL2Projection(pmesh, smesh, p, exact_soln);
+         CHECK(std::abs(error[1] - error[0]) < 1e-9);
       }
    }
 } // 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