Remove problematic inline.

Wrong assert code in ParGF::ProjectBdrCoefficient

.gitignore

@@ -300,7 +300,7 @@ miniapps/tools/convert-dc
 miniapps/tools/lor-transfer
 miniapps/tools/plor-transfer
 miniapps/tools/get-values
-miniapps/tools/check-tmop-metric
+miniapps/tools/tmop-check-metric
 miniapps/tools/tmop-metric-magnitude
 miniapps/tools/nodal-transfer
 miniapps/tools/ParaView

CHANGELOG

@@ -87,6 +87,8 @@ Linear and nonlinear solvers
 
 - Added HIP support to the PETSc and SUNDIALS interfaces.
 
+- Efficient GPU-accelerated LOR assembly now supports surface meshes.
+
 New and updated examples and miniapps
 -------------------------------------
 - Added a new H(div) solver miniapp demonstrating the use of a matrix-free

CMakeLists.txt

@@ -641,16 +641,34 @@ if (NOT ("${PROJECT_SOURCE_DIR}" STREQUAL "${PROJECT_BINARY_DIR}"))
   foreach(Header mfem.hpp mfem-performance.hpp)
     message(STATUS
       "Writing substitute header --> \"${Header}\"")
-    file(WRITE "${PROJECT_BINARY_DIR}/${Header}"
+    file(WRITE "${PROJECT_BINARY_DIR}/${Header}.tmp"
 "// Auto-generated file.
 #define MFEM_CONFIG_FILE \"${PROJECT_BINARY_DIR}/config/_config.hpp\"
 #include \"${PROJECT_SOURCE_DIR}/${Header}\"
 ")
+
+    execute_process(COMMAND ${CMAKE_COMMAND} -E copy_if_different 
+      "${PROJECT_BINARY_DIR}/${Header}.tmp"
+      "${PROJECT_BINARY_DIR}/${Header}"
+    )
+    execute_process(COMMAND ${CMAKE_COMMAND} -E remove
+      "${PROJECT_BINARY_DIR}/${Header}.tmp"
+    )
+
     # This version will be installed in the top include directory:
-    file(WRITE "${PROJECT_BINARY_DIR}/InstallHeaders/${Header}"
+    file(WRITE "${PROJECT_BINARY_DIR}/InstallHeaders/${Header}.tmp"
 "// Auto-generated file.
 #include \"mfem/${Header}\"
 ")
+
+    execute_process(COMMAND ${CMAKE_COMMAND} -E copy_if_different 
+      "${PROJECT_BINARY_DIR}/InstallHeaders/${Header}.tmp"
+      "${PROJECT_BINARY_DIR}/InstallHeaders/${Header}"
+    )
+    execute_process(COMMAND ${CMAKE_COMMAND} -E remove
+      "${PROJECT_BINARY_DIR}/InstallHeaders/${Header}.tmp"
+    )
+
   endforeach()
 endif()
 

INSTALL

@@ -659,8 +659,7 @@ The specific libraries and their options are:
   requires the PT-Scotch and Scalapack libraries as well as ParMETIS, which
   includes METIS 5 in its distribution. Starting with STRUMPACK v2.2.0, ParMETIS
   and PT-Scotch are optional dependencies.
-  The support for STRUMPACK was added in MFEM v3.3.2 and it requires STRUMPACK
-  2.0.0 or later.
+  The support for STRUMPACK was added in MFEM v3.3.2.
   URL: http://portal.nersc.gov/project/sparse/strumpack
   Options: STRUMPACK_OPT, STRUMPACK_LIB.
   Versions: STRUMPACK >= 3.0.0.

config/defaults.cmake

@@ -157,7 +157,8 @@ set(STRUMPACK_DIR "${MFEM_DIR}/../STRUMPACK-build" CACHE PATH
 # STRUMPACK may also depend on "OpenMP", depending on how it was compiled.
 # Starting with v2.2.0 of STRUMPACK, ParMETIS and Scotch are optional.
 set(STRUMPACK_REQUIRED_PACKAGES "MPI" "MPI_Fortran" "ParMETIS" "METIS"
-    "ScaLAPACK" "Scotch/ptscotch/ptscotcherr/scotch/scotcherr" CACHE STRING
+    "Scotch/ptscotch/ptscotcherr/scotch/scotcherr"
+    "ScaLAPACK" "LAPACK" "BLAS" CACHE STRING
     "Additional packages required by STRUMPACK.")
 # If the MPI package does not find all required Fortran libraries:
 # set(STRUMPACK_REQUIRED_LIBRARIES "gfortran" "mpi_mpifh" CACHE STRING

examples/ex11p.cpp

@@ -262,12 +262,13 @@ int main(int argc, char *argv[])
 #ifdef MFEM_USE_STRUMPACK
       if (sp_solver)
       {
-         STRUMPACKSolver * strumpack = new STRUMPACKSolver(argc, argv, MPI_COMM_WORLD);
+         STRUMPACKSolver * strumpack = new STRUMPACKSolver(MPI_COMM_WORLD, argc, argv);
          strumpack->SetPrintFactorStatistics(true);
          strumpack->SetPrintSolveStatistics(false);
          strumpack->SetKrylovSolver(strumpack::KrylovSolver::DIRECT);
          strumpack->SetReorderingStrategy(strumpack::ReorderingStrategy::METIS);
-         strumpack->DisableMatching();
+         strumpack->SetMatching(strumpack::MatchingJob::NONE);
+         strumpack->SetCompression(strumpack::CompressionType::NONE);
          strumpack->SetOperator(*Arow);
          strumpack->SetFromCommandLine();
          precond = strumpack;

examples/ex25p.cpp

@@ -170,6 +170,7 @@ int main(int argc, char *argv[])
    bool herm_conv = true;
    bool slu_solver  = false;
    bool mumps_solver = false;
+   bool strumpack_solver = false;
    bool visualization = 1;
    bool pa = false;
    const char *device_config = "cpu";
@@ -200,6 +201,11 @@ int main(int argc, char *argv[])
 #ifdef MFEM_USE_MUMPS
    args.AddOption(&mumps_solver, "-mumps", "--mumps-solver", "-no-mumps",
                   "--no-mumps-solver", "Use the MUMPS Solver.");
+#endif
+#ifdef MFEM_USE_STRUMPACK
+   args.AddOption(&strumpack_solver, "-strumpack", "--strumpack-solver",
+                  "-no-strumpack", "--no-strumpack-solver",
+                  "Use the STRUMPACK Solver.");
 #endif
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
@@ -209,13 +215,14 @@ int main(int argc, char *argv[])
    args.AddOption(&device_config, "-d", "--device",
                   "Device configuration string, see Device::Configure().");
    args.Parse();
-   if (slu_solver && mumps_solver)
+   if (slu_solver + mumps_solver + strumpack_solver > 1)
    {
       if (myid == 0)
-         cout << "WARNING: Both SuperLU and MUMPS have been selected,"
-              << " please choose either one." << endl
+         cout << "WARNING: More than one of SuperLU, MUMPS, and STRUMPACK have"
+              << " been selected, please choose only one." << endl
               << "         Defaulting to SuperLU." << endl;
       mumps_solver = false;
+      strumpack_solver = false;
    }
 
    if (iprob > 4) { iprob = 4; }
@@ -474,6 +481,24 @@ int main(int argc, char *argv[])
       delete A;
    }
 #endif
+#ifdef MFEM_USE_STRUMPACK
+   if (!pa && strumpack_solver)
+   {
+      HypreParMatrix *A = Ah.As<ComplexHypreParMatrix>()->GetSystemMatrix();
+      STRUMPACKRowLocMatrix SA(*A);
+      STRUMPACKSolver strumpack(MPI_COMM_WORLD, argc, argv);
+      strumpack.SetPrintFactorStatistics(false);
+      strumpack.SetPrintSolveStatistics(false);
+      strumpack.SetKrylovSolver(strumpack::KrylovSolver::DIRECT);
+      strumpack.SetReorderingStrategy(strumpack::ReorderingStrategy::METIS);
+      strumpack.SetMatching(strumpack::MatchingJob::NONE);
+      strumpack.SetCompression(strumpack::CompressionType::NONE);
+      strumpack.SetFromCommandLine();
+      strumpack.SetOperator(SA);
+      strumpack.Mult(B, X);
+      delete A;
+   }
+#endif
 #ifdef MFEM_USE_MUMPS
    if (!pa && mumps_solver)
    {
@@ -493,7 +518,7 @@ int main(int argc, char *argv[])
    //
    //    In PML:   1/mu (abs(1/det(J) J^T J) Curl E, Curl F)
    //              + omega^2 * epsilon (abs(det(J) * (J^T J)^-1) * E, F)
-   if (pa || (!slu_solver && !mumps_solver))
+   if (pa || (!slu_solver && !mumps_solver && !strumpack_solver))
    {
       ConstantCoefficient absomeg(pow(omega, 2) * epsilon);
       RestrictedCoefficient restr_absomeg(absomeg,attr);

examples/petsc/ex11p.cpp

@@ -273,12 +273,13 @@ int main(int argc, char *argv[])
 #ifdef MFEM_USE_STRUMPACK
          if (sp_solver)
          {
-            STRUMPACKSolver * strumpack = new STRUMPACKSolver(argc, argv, MPI_COMM_WORLD);
+            STRUMPACKSolver * strumpack = new STRUMPACKSolver(MPI_COMM_WORLD, argc, argv);
             strumpack->SetPrintFactorStatistics(true);
             strumpack->SetPrintSolveStatistics(false);
             strumpack->SetKrylovSolver(strumpack::KrylovSolver::DIRECT);
             strumpack->SetReorderingStrategy(strumpack::ReorderingStrategy::METIS);
-            strumpack->DisableMatching();
+            strumpack->SetMatching(strumpack::MatchingJob::NONE);
+            strumpack->SetCompression(strumpack::CompressionType::NONE);
             strumpack->SetOperator(*Arow);
             strumpack->SetFromCommandLine();
             precond = strumpack;

fem/CMakeLists.txt

@@ -96,9 +96,6 @@ set(SRCS
   lor/lor_ads.cpp
   lor/lor_ams.cpp
   lor/lor_batched.cpp
-  lor/lor_h1.cpp
-  lor/lor_nd.cpp
-  lor/lor_rt.cpp
   multigrid.cpp
   nonlinearform.cpp
   nonlinearform_ext.cpp

fem/bilinearform.cpp

@@ -101,6 +101,7 @@ BilinearForm::BilinearForm (FiniteElementSpace * f, BilinearForm * bf, int ps)
 
    // Copy the pointers to the integrators
    domain_integs = bf->domain_integs;
+   domain_integs_marker = bf->domain_integs_marker;
 
    boundary_integs = bf->boundary_integs;
    boundary_integs_marker = bf->boundary_integs_marker;
@@ -433,7 +434,6 @@ void BilinearForm::Assemble(int skip_zeros)
       // Element-wise integration
       for (int i = 0; i < fes -> GetNE(); i++)
       {
-         doftrans = fes->GetElementVDofs(i, vdofs);
          if (element_matrices)
          {
             elmat_p = &(*element_matrices)(i);
@@ -441,6 +441,9 @@ void BilinearForm::Assemble(int skip_zeros)
          else
          {
             const int elem_attr = fes->GetMesh()->GetAttribute(i);
+            doftrans = fes->GetElementVDofs(i, vdofs);
+            eltrans = fes->GetElementTransformation(i);
+
             elmat.SetSize(0);
             for (int k = 0; k < domain_integs.Size(); k++)
             {
@@ -448,9 +451,8 @@ void BilinearForm::Assemble(int skip_zeros)
                     (*(domain_integs_marker[k]))[elem_attr-1] == 1)
                    && !domain_integs[k]->Patchwise())
                {
-                  const FiniteElement &fe = *fes->GetFE(i);
-                  eltrans = fes->GetElementTransformation(i);
-                  domain_integs[k]->AssembleElementMatrix(fe, *eltrans, elemmat);
+                  domain_integs[k]->AssembleElementMatrix(*fes->GetFE(i),
+                                                          *eltrans, elemmat);
                   if (elmat.Size() == 0)
                   {
                      elmat = elemmat;
@@ -1222,11 +1224,14 @@ MixedBilinearForm::MixedBilinearForm (FiniteElementSpace *tr_fes,
 
    // Copy the pointers to the integrators
    domain_integs = mbf->domain_integs;
-   boundary_integs = mbf->boundary_integs;
-   trace_face_integs = mbf->trace_face_integs;
-   boundary_trace_face_integs = mbf->boundary_trace_face_integs;
+   domain_integs_marker = mbf->domain_integs_marker;
 
+   boundary_integs = mbf->boundary_integs;
    boundary_integs_marker = mbf->boundary_integs_marker;
+
+   trace_face_integs = mbf->trace_face_integs;
+
+   boundary_trace_face_integs = mbf->boundary_trace_face_integs;
    boundary_trace_face_integs_marker = mbf->boundary_trace_face_integs_marker;
 
    assembly = AssemblyLevel::LEGACY;
@@ -1349,6 +1354,14 @@ void MixedBilinearForm::GetBlocks(Array2D<SparseMatrix *> &blocks) const
 void MixedBilinearForm::AddDomainIntegrator (BilinearFormIntegrator * bfi)
 {
    domain_integs.Append (bfi);
+   domain_integs_marker.Append(NULL); // NULL marker means apply everywhere
+}
+
+void MixedBilinearForm::AddDomainIntegrator (BilinearFormIntegrator * bfi,
+                                             Array<int> &elem_marker)
+{
+   domain_integs.Append (bfi);
+   domain_integs_marker.Append(&elem_marker);
 }
 
 void MixedBilinearForm::AddBoundaryIntegrator (BilinearFormIntegrator * bfi)
@@ -1383,7 +1396,7 @@ void MixedBilinearForm::AddBdrTraceFaceIntegrator(BilinearFormIntegrator *bfi,
    boundary_trace_face_integs_marker.Append(&bdr_marker);
 }
 
-void MixedBilinearForm::Assemble (int skip_zeros)
+void MixedBilinearForm::Assemble(int skip_zeros)
 {
    if (ext)
    {
@@ -1405,8 +1418,20 @@ void MixedBilinearForm::Assemble (int skip_zeros)
 
    if (domain_integs.Size())
    {
+      for (int k = 0; k < domain_integs.Size(); k++)
+      {
+         if (domain_integs_marker[k] != NULL)
+         {
+            MFEM_VERIFY(domain_integs_marker[k]->Size() ==
+                        (mesh->attributes.Size() ? mesh->attributes.Max() : 0),
+                        "invalid element marker for domain integrator #"
+                        << k << ", counting from zero");
+         }
+      }
+
       for (int i = 0; i < test_fes -> GetNE(); i++)
       {
+         const int elem_attr = mesh->GetAttribute(i);
          dom_dof_trans = trial_fes -> GetElementVDofs (i, trial_vdofs);
          ran_dof_trans = test_fes  -> GetElementVDofs (i, test_vdofs);
          eltrans = test_fes -> GetElementTransformation (i);
@@ -1415,10 +1440,14 @@ void MixedBilinearForm::Assemble (int skip_zeros)
          elmat = 0.0;
          for (int k = 0; k < domain_integs.Size(); k++)
          {
-            domain_integs[k] -> AssembleElementMatrix2 (*trial_fes -> GetFE(i),
-                                                        *test_fes  -> GetFE(i),
-                                                        *eltrans, elemmat);
-            elmat += elemmat;
+            if (domain_integs_marker[k] == NULL ||
+                (*(domain_integs_marker[k]))[elem_attr-1] == 1)
+            {
+               domain_integs[k] -> AssembleElementMatrix2 (*trial_fes -> GetFE(i),
+                                                           *test_fes  -> GetFE(i),
+                                                           *eltrans, elemmat);
+               elmat += elemmat;
+            }
          }
          if (ran_dof_trans || dom_dof_trans)
          {
@@ -1941,41 +1970,56 @@ void DiscreteLinearOperator::Assemble(int skip_zeros)
       return;
    }
 
-   Array<int> dom_vdofs, ran_vdofs;
-   ElementTransformation *T;
+   ElementTransformation *eltrans;
    DofTransformation * dom_dof_trans;
    DofTransformation * ran_dof_trans;
-   const FiniteElement *dom_fe, *ran_fe;
-   DenseMatrix totelmat, elmat;
+   DenseMatrix elmat;
+
+   Mesh *mesh = test_fes->GetMesh();
 
    if (mat == NULL)
    {
       mat = new SparseMatrix(height, width);
    }
 
-   if (domain_integs.Size() > 0)
+   if (domain_integs.Size())
    {
+      for (int k = 0; k < domain_integs.Size(); k++)
+      {
+         if (domain_integs_marker[k] != NULL)
+         {
+            MFEM_VERIFY(domain_integs_marker[k]->Size() ==
+                        (mesh->attributes.Size() ? mesh->attributes.Max() : 0),
+                        "invalid element marker for domain integrator #"
+                        << k << ", counting from zero");
+         }
+      }
+
       for (int i = 0; i < test_fes->GetNE(); i++)
       {
-         dom_dof_trans = trial_fes->GetElementVDofs(i, dom_vdofs);
-         ran_dof_trans = test_fes->GetElementVDofs(i, ran_vdofs);
-         T = test_fes->GetElementTransformation(i);
-         dom_fe = trial_fes->GetFE(i);
-         ran_fe = test_fes->GetFE(i);
+         const int elem_attr = mesh->GetAttribute(i);
+         dom_dof_trans = trial_fes->GetElementVDofs(i, trial_vdofs);
+         ran_dof_trans = test_fes->GetElementVDofs(i, test_vdofs);
+         eltrans = test_fes->GetElementTransformation(i);
 
-         domain_integs[0]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T,
-                                                  totelmat);
-         for (int j = 1; j < domain_integs.Size(); j++)
+         elmat.SetSize(test_vdofs.Size(), trial_vdofs.Size());
+         elmat = 0.0;
+         for (int k = 0; k < domain_integs.Size(); k++)
          {
-            domain_integs[j]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T,
-                                                     elmat);
-            totelmat += elmat;
+            if (domain_integs_marker[k] == NULL ||
+                (*(domain_integs_marker[k]))[elem_attr-1] == 1)
+            {
+               domain_integs[k]->AssembleElementMatrix2(*trial_fes->GetFE(i),
+                                                        *test_fes->GetFE(i),
+                                                        *eltrans, elemmat);
+               elmat += elemmat;
+            }
          }
          if (ran_dof_trans || dom_dof_trans)
          {
-            TransformPrimal(ran_dof_trans, dom_dof_trans, totelmat);
+            TransformPrimal(ran_dof_trans, dom_dof_trans, elemmat);
          }
-         mat->SetSubMatrix(ran_vdofs, dom_vdofs, totelmat, skip_zeros);
+         mat->SetSubMatrix(test_vdofs, trial_vdofs, elemmat, skip_zeros);
       }
    }
 
@@ -1984,21 +2028,20 @@ void DiscreteLinearOperator::Assemble(int skip_zeros)
       const int nfaces = test_fes->GetMesh()->GetNumFaces();
       for (int i = 0; i < nfaces; i++)
       {
-         trial_fes->GetFaceVDofs(i, dom_vdofs);
-         test_fes->GetFaceVDofs(i, ran_vdofs);
-         T = test_fes->GetMesh()->GetFaceTransformation(i);
-         dom_fe = trial_fes->GetFaceElement(i);
-         ran_fe = test_fes->GetFaceElement(i);
+         trial_fes->GetFaceVDofs(i, trial_vdofs);
+         test_fes->GetFaceVDofs(i, test_vdofs);
+         eltrans = test_fes->GetMesh()->GetFaceTransformation(i);
 
-         trace_face_integs[0]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T,
-                                                      totelmat);
-         for (int j = 1; j < trace_face_integs.Size(); j++)
+         elmat.SetSize(test_vdofs.Size(), trial_vdofs.Size());
+         elmat = 0.0;
+         for (int k = 0; k < trace_face_integs.Size(); k++)
          {
-            trace_face_integs[j]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T,
-                                                         elmat);
-            totelmat += elmat;
+            trace_face_integs[k]->AssembleElementMatrix2(*trial_fes->GetFaceElement(i),
+                                                         *test_fes->GetFaceElement(i),
+                                                         *eltrans, elemmat);
+            elmat += elemmat;
          }
-         mat->SetSubMatrix(ran_vdofs, dom_vdofs, totelmat, skip_zeros);
+         mat->SetSubMatrix(test_vdofs, trial_vdofs, elmat, skip_zeros);
       }
    }
 }

fem/bilinearform.hpp

@@ -100,7 +100,7 @@ protected:
    /// Includes all by default.
    /// 0 - ignore attribute
    /// 1 - include attribute
-   Array<Array<int>*>             domain_integs_marker;
+   Array<Array<int>*> domain_integs_marker; ///< Entries are not owned.
 
    /// Set of Boundary Integrators to be applied.
    Array<BilinearFormIntegrator*> boundary_integs;
@@ -722,10 +722,13 @@ protected:
 
    /// Domain integrators.
    Array<BilinearFormIntegrator*> domain_integs;
+   /// Entries are not owned.
+   Array<Array<int>*> domain_integs_marker;
 
    /// Boundary integrators.
    Array<BilinearFormIntegrator*> boundary_integs;
-   Array<Array<int>*> boundary_integs_marker; ///< Entries are not owned.
+   /// Entries are not owned.
+   Array<Array<int>*> boundary_integs_marker;
 
    /// Trace face (skeleton) integrators.
    Array<BilinearFormIntegrator*> trace_face_integs;
@@ -805,12 +808,16 @@ public:
    /// Adds a domain integrator. Assumes ownership of @a bfi.
    void AddDomainIntegrator(BilinearFormIntegrator *bfi);
 
+   /// Adds a domain integrator. Assumes ownership of @a bfi.
+   void AddDomainIntegrator(BilinearFormIntegrator *bfi,
+                            Array<int> &elem_marker);
+
    /// Adds a boundary integrator. Assumes ownership of @a bfi.
    void AddBoundaryIntegrator(BilinearFormIntegrator *bfi);
 
    /// Adds a boundary integrator. Assumes ownership of @a bfi.
-   void AddBoundaryIntegrator (BilinearFormIntegrator * bfi,
-                               Array<int> &bdr_marker);
+   void AddBoundaryIntegrator(BilinearFormIntegrator * bfi,
+                              Array<int> &bdr_marker);
 
    /** @brief Add a trace face integrator. Assumes ownership of @a bfi.
 
@@ -820,14 +827,18 @@ public:
    void AddTraceFaceIntegrator(BilinearFormIntegrator *bfi);
 
    /// Adds a boundary trace face integrator. Assumes ownership of @a bfi.
-   void AddBdrTraceFaceIntegrator (BilinearFormIntegrator * bfi);
+   void AddBdrTraceFaceIntegrator(BilinearFormIntegrator * bfi);
 
    /// Adds a boundary trace face integrator. Assumes ownership of @a bfi.
-   void AddBdrTraceFaceIntegrator (BilinearFormIntegrator * bfi,
-                                   Array<int> &bdr_marker);
+   void AddBdrTraceFaceIntegrator(BilinearFormIntegrator * bfi,
+                                  Array<int> &bdr_marker);
 
    /// Access all integrators added with AddDomainIntegrator().
    Array<BilinearFormIntegrator*> *GetDBFI() { return &domain_integs; }
+   /** @brief Access all domain markers added with AddDomainIntegrator().
+       If no marker was specified when the integrator was added, the
+       corresponding pointer (to Array<int>) will be NULL. */
+   Array<Array<int>*> *GetDBFI_Marker() { return &domain_integs_marker; }
 
    /// Access all integrators added with AddBoundaryIntegrator().
    Array<BilinearFormIntegrator*> *GetBBFI() { return &boundary_integs; }
@@ -1065,6 +1076,9 @@ public:
    /// Adds a domain interpolator. Assumes ownership of @a di.
    void AddDomainInterpolator(DiscreteInterpolator *di)
    { AddDomainIntegrator(di); }
+   void AddDomainInterpolator(DiscreteInterpolator *di,
+                              Array<int> &elem_marker)
+   { AddDomainIntegrator(di, elem_marker); }
 
    /// Adds a trace face interpolator. Assumes ownership of @a di.
    void AddTraceFaceInterpolator(DiscreteInterpolator *di)
@@ -1072,6 +1086,7 @@ public:
 
    /// Access all interpolators added with AddDomainInterpolator().
    Array<BilinearFormIntegrator*> *GetDI() { return &domain_integs; }
+   Array<Array<int>*> *GetDI_Marker() { return &domain_integs_marker; }
 
    /// Set the desired assembly level. The default is AssemblyLevel::FULL.
    /** This method must be called before assembly. */

fem/bilinearform_ext.cpp

@@ -303,7 +303,7 @@ void PABilinearFormExtension::SetupRestrictionOperators(const L2FaceValues m)
       std::unordered_map<int,int> f_to_be;
       for (int i = 0; i < mesh.GetNBE(); ++i)
       {
-         const int f = mesh.GetBdrElementEdgeIndex(i);
+         const int f = mesh.GetBdrElementFaceIndex(i);
          f_to_be[f] = i;
       }
       const int nf_bdr = trial_fes->GetNFbyType(FaceType::Boundary);

fem/doftrans.cpp

@@ -14,175 +14,166 @@
 namespace mfem
 {
 
+void DofTransformation::TransformPrimal(double *v) const
+{
+   MFEM_ASSERT(dof_trans_,
+               "DofTransformation has no local transformation, call "
+               "SetDofTransformation first!");
+   int size = dof_trans_->Size();
+
+   if (vdim_ == 1 || (Ordering::Type)ordering_ == Ordering::byNODES)
+   {
+      for (int i=0; i<vdim_; i++)
+      {
+         dof_trans_->TransformPrimal(Fo_, &v[i*size]);
+      }
+   }
+   else
+   {
+      Vector vec(size);
+      for (int i=0; i<vdim_; i++)
+      {
+         for (int j=0; j<size; j++)
+         {
+            vec(j) = v[j*vdim_+i];
+         }
+         dof_trans_->TransformPrimal(Fo_, vec);
+         for (int j=0; j<size; j++)
+         {
+            v[j*vdim_+i] = vec(j);
+         }
+      }
+   }
+}
+
+void DofTransformation::InvTransformPrimal(double *v) const
+{
+   MFEM_ASSERT(dof_trans_,
+               "DofTransformation has no local transformation, call "
+               "SetDofTransformation first!");
+   int size = dof_trans_->Height();
+
+   if (vdim_ == 1 || (Ordering::Type)ordering_ == Ordering::byNODES)
+   {
+      for (int i=0; i<vdim_; i++)
+      {
+         dof_trans_->InvTransformPrimal(Fo_, &v[i*size]);
+      }
+   }
+   else
+   {
+      Vector vec(size);
+      for (int i=0; i<vdim_; i++)
+      {
+         for (int j=0; j<size; j++)
+         {
+            vec(j) = v[j*vdim_+i];
+         }
+         dof_trans_->InvTransformPrimal(Fo_, vec);
+         for (int j=0; j<size; j++)
+         {
+            v[j*vdim_+i] = vec(j);
+         }
+      }
+   }
+}
+
+void DofTransformation::TransformDual(double *v) const
+{
+   MFEM_ASSERT(dof_trans_,
+               "DofTransformation has no local transformation, call "
+               "SetDofTransformation first!");
+   int size = dof_trans_->Size();
+
+   if (vdim_ == 1 || (Ordering::Type)ordering_ == Ordering::byNODES)
+   {
+      for (int i=0; i<vdim_; i++)
+      {
+         dof_trans_->TransformDual(Fo_, &v[i*size]);
+      }
+   }
+   else
+   {
+      Vector vec(size);
+      for (int i=0; i<vdim_; i++)
+      {
+         for (int j=0; j<size; j++)
+         {
+            vec(j) = v[j*vdim_+i];
+         }
+         dof_trans_->TransformDual(Fo_, vec);
+         for (int j=0; j<size; j++)
+         {
+            v[j*vdim_+i] = vec(j);
+         }
+      }
+   }
+}
+
+void DofTransformation::InvTransformDual(double *v) const
+{
+   MFEM_ASSERT(dof_trans_,
+               "DofTransformation has no local transformation, call "
+               "SetDofTransformation first!");
+   int size = dof_trans_->Size();
+
+   if (vdim_ == 1 || (Ordering::Type)ordering_ == Ordering::byNODES)
+   {
+      for (int i=0; i<vdim_; i++)
+      {
+         dof_trans_->InvTransformDual(Fo_, &v[i*size]);
+      }
+   }
+   else
+   {
+      Vector vec(size);
+      for (int i=0; i<vdim_; i++)
+      {
+         for (int j=0; j<size; j++)
+         {
+            vec(j) = v[j*vdim_+i];
+         }
+         dof_trans_->InvTransformDual(Fo_, vec);
+         for (int j=0; j<size; j++)
+         {
+            v[j*vdim_+i] = vec(j);
+         }
+      }
+   }
+}
+
 void TransformPrimal(const DofTransformation *ran_dof_trans,
                      const DofTransformation *dom_dof_trans,
                      DenseMatrix &elmat)
 {
-   if (ran_dof_trans && dom_dof_trans)
-   {
-      ran_dof_trans->TransformPrimalCols(elmat);
-      dom_dof_trans->TransformDualRows(elmat);
-   }
-   else if (ran_dof_trans)
+   // No action if both transformations are NULL
+   if (ran_dof_trans)
    {
       ran_dof_trans->TransformPrimalCols(elmat);
    }
-   else if (dom_dof_trans)
+   if (dom_dof_trans)
    {
       dom_dof_trans->TransformDualRows(elmat);
    }
-   else
-   {
-      // If both transformations are NULL this function should not be called
-   }
 }
 
 void TransformDual(const DofTransformation *ran_dof_trans,
                    const DofTransformation *dom_dof_trans,
                    DenseMatrix &elmat)
 {
-   if (ran_dof_trans && dom_dof_trans)
-   {
-      ran_dof_trans->TransformDualCols(elmat);
-      dom_dof_trans->TransformDualRows(elmat);
-   }
-   else if (ran_dof_trans)
+   // No action if both transformations are NULL
+   if (ran_dof_trans)
    {
       ran_dof_trans->TransformDualCols(elmat);
    }
-   else if (dom_dof_trans)
+   if (dom_dof_trans)
    {
       dom_dof_trans->TransformDualRows(elmat);
    }
-   else
-   {
-      // If both transformations are NULL this function should not be called
-   }
-}
-
-void StatelessVDofTransformation::TransformPrimal(const Array<int> & face_ori,
-                                                  double *v) const
-{
-   int size = sdoftrans_->Size();
-
-   if ((Ordering::Type)ordering_ == Ordering::byNODES || vdim_ == 1)
-   {
-      for (int i=0; i<vdim_; i++)
-      {
-         sdoftrans_->TransformPrimal(face_ori, &v[i*size]);
-      }
-   }
-   else
-   {
-      Vector vec(size);
-      for (int i=0; i<vdim_; i++)
-      {
-         for (int j=0; j<size; j++)
-         {
-            vec(j) = v[j*vdim_+i];
-         }
-         sdoftrans_->TransformPrimal(face_ori, vec);
-         for (int j=0; j<size; j++)
-         {
-            v[j*vdim_+i] = vec(j);
-         }
-      }
-   }
-}
-
-void StatelessVDofTransformation::InvTransformPrimal(
-   const Array<int> & face_ori,
-   double *v) const
-{
-   int size = sdoftrans_->Height();
-
-   if ((Ordering::Type)ordering_ == Ordering::byNODES)
-   {
-      for (int i=0; i<vdim_; i++)
-      {
-         sdoftrans_->InvTransformPrimal(face_ori, &v[i*size]);
-      }
-   }
-   else
-   {
-      Vector vec(size);
-      for (int i=0; i<vdim_; i++)
-      {
-         for (int j=0; j<size; j++)
-         {
-            vec(j) = v[j*vdim_+i];
-         }
-         sdoftrans_->InvTransformPrimal(face_ori, vec);
-         for (int j=0; j<size; j++)
-         {
-            v[j*vdim_+i] = vec(j);
-         }
-      }
-   }
-}
-
-void StatelessVDofTransformation::TransformDual(const Array<int> & face_ori,
-                                                double *v) const
-{
-   int size = sdoftrans_->Size();
-
-   if ((Ordering::Type)ordering_ == Ordering::byNODES)
-   {
-      for (int i=0; i<vdim_; i++)
-      {
-         sdoftrans_->TransformDual(face_ori, &v[i*size]);
-      }
-   }
-   else
-   {
-      Vector vec(size);
-      for (int i=0; i<vdim_; i++)
-      {
-         for (int j=0; j<size; j++)
-         {
-            vec(j) = v[j*vdim_+i];
-         }
-         sdoftrans_->TransformDual(face_ori, vec);
-         for (int j=0; j<size; j++)
-         {
-            v[j*vdim_+i] = vec(j);
-         }
-      }
-   }
-}
-
-void StatelessVDofTransformation::InvTransformDual(const Array<int> & face_ori,
-                                                   double *v) const
-{
-   int size = sdoftrans_->Size();
-
-   if ((Ordering::Type)ordering_ == Ordering::byNODES)
-   {
-      for (int i=0; i<vdim_; i++)
-      {
-         sdoftrans_->InvTransformDual(face_ori, &v[i*size]);
-      }
-   }
-   else
-   {
-      Vector vec(size);
-      for (int i=0; i<vdim_; i++)
-      {
-         for (int j=0; j<size; j++)
-         {
-            vec(j) = v[j*vdim_+i];
-         }
-         sdoftrans_->InvTransformDual(face_ori, vec);
-         for (int j=0; j<size; j++)
-         {
-            v[j*vdim_+i] = vec(j);
-         }
-      }
-   }
 }
 
 // ordering (i0j0, i1j0, i0j1, i1j1), each row is a column major matrix
-const double ND_StatelessDofTransformation::T_data[24] =
+const double ND_DofTransformation::T_data[24] =
 {
    1.0,  0.0,  0.0,  1.0,
    -1.0, -1.0,  0.0,  1.0,
@@ -192,11 +183,11 @@ const double ND_StatelessDofTransformation::T_data[24] =
    0.0,  1.0,  1.0,  0.0
 };
 
-const DenseTensor ND_StatelessDofTransformation
-::T(const_cast<double*>(ND_StatelessDofTransformation::T_data), 2, 2, 6);
+const DenseTensor ND_DofTransformation
+::T(const_cast<double *>(ND_DofTransformation::T_data), 2, 2, 6);
 
 // ordering (i0j0, i1j0, i0j1, i1j1), each row is a column major matrix
-const double ND_StatelessDofTransformation::TInv_data[24] =
+const double ND_DofTransformation::TInv_data[24] =
 {
    1.0,  0.0,  0.0,  1.0,
    -1.0, -1.0,  0.0,  1.0,
@@ -206,12 +197,11 @@ const double ND_StatelessDofTransformation::TInv_data[24] =
    0.0,  1.0,  1.0,  0.0
 };
 
-const DenseTensor ND_StatelessDofTransformation
-::TInv(const_cast<double*>(TInv_data), 2, 2, 6);
+const DenseTensor ND_DofTransformation
+::TInv(const_cast<double *>(TInv_data), 2, 2, 6);
 
-ND_StatelessDofTransformation::ND_StatelessDofTransformation(int size, int p,
-                                                             int num_edges,
-                                                             int num_tri_faces)
+ND_DofTransformation::ND_DofTransformation(int size, int p, int num_edges,
+                                           int num_tri_faces)
    : StatelessDofTransformation(size)
    , order(p)
    , nedofs(p)
@@ -221,18 +211,19 @@ ND_StatelessDofTransformation::ND_StatelessDofTransformation(int size, int p,
 {
 }
 
-void ND_StatelessDofTransformation::TransformPrimal(const Array<int> & Fo,
-                                                    double *v) const
+void ND_DofTransformation::TransformPrimal(const Array<int> & Fo,
+                                           double *v) const
 {
    // Return immediately when no face DoFs are present
-   if (nfdofs < 2) { return; }
+   if (IsIdentity()) { return; }
 
    MFEM_VERIFY(Fo.Size() >= nfaces,
                "Face orientation array is shorter than the number of faces in "
-               "ND_StatelessDofTransformation");
+               "ND_DofTransformation");
 
    double data[2];
    Vector v2(data, 2);
+   DenseMatrix T2;
 
    // Transform face DoFs
    for (int f=0; f<nfaces; f++)
@@ -240,23 +231,25 @@ void ND_StatelessDofTransformation::TransformPrimal(const Array<int> & Fo,
       for (int i=0; i<nfdofs/2; i++)
       {
          v2 = &v[nedges*nedofs + f*nfdofs + 2*i];
-         T(Fo[f]).Mult(v2, &v[nedges*nedofs + f*nfdofs + 2*i]);
+         T2.UseExternalData(const_cast<double *>(T.GetData(Fo[f])), 2, 2);
+         T2.Mult(v2, &v[nedges*nedofs + f*nfdofs + 2*i]);
       }
    }
 }
 
-void ND_StatelessDofTransformation::InvTransformPrimal(const Array<int> & Fo,
-                                                       double *v) const
+void ND_DofTransformation::InvTransformPrimal(const Array<int> & Fo,
+                                              double *v) const
 {
    // Return immediately when no face DoFs are present
-   if (nfdofs < 2) { return; }
+   if (IsIdentity()) { return; }
 
    MFEM_VERIFY(Fo.Size() >= nfaces,
                "Face orientation array is shorter than the number of faces in "
-               "ND_StatelessDofTransformation");
+               "ND_DofTransformation");
 
    double data[2];
    Vector v2(data, 2);
+   DenseMatrix T2Inv;
 
    // Transform face DoFs
    for (int f=0; f<nfaces; f++)
@@ -264,23 +257,24 @@ void ND_StatelessDofTransformation::InvTransformPrimal(const Array<int> & Fo,
       for (int i=0; i<nfdofs/2; i++)
       {
          v2 = &v[nedges*nedofs + f*nfdofs + 2*i];
-         TInv(Fo[f]).Mult(v2, &v[nedges*nedofs + f*nfdofs + 2*i]);
+         T2Inv.UseExternalData(const_cast<double *>(TInv.GetData(Fo[f])), 2, 2);
+         T2Inv.Mult(v2, &v[nedges*nedofs + f*nfdofs + 2*i]);
       }
    }
 }
 
-void ND_StatelessDofTransformation::TransformDual(const Array<int> & Fo,
-                                                  double *v) const
+void ND_DofTransformation::TransformDual(const Array<int> & Fo, double *v) const
 {
    // Return immediately when no face DoFs are present
-   if (nfdofs < 2) { return; }
+   if (IsIdentity()) { return; }
 
    MFEM_VERIFY(Fo.Size() >= nfaces,
                "Face orientation array is shorter than the number of faces in "
-               "ND_StatelessDofTransformation");
+               "ND_DofTransformation");
 
    double data[2];
    Vector v2(data, 2);
+   DenseMatrix T2Inv;
 
    // Transform face DoFs
    for (int f=0; f<nfaces; f++)
@@ -288,23 +282,25 @@ void ND_StatelessDofTransformation::TransformDual(const Array<int> & Fo,
       for (int i=0; i<nfdofs/2; i++)
       {
          v2 = &v[nedges*nedofs + f*nfdofs + 2*i];
-         TInv(Fo[f]).MultTranspose(v2, &v[nedges*nedofs + f*nfdofs + 2*i]);
+         T2Inv.UseExternalData(const_cast<double *>(TInv.GetData(Fo[f])), 2, 2);
+         T2Inv.MultTranspose(v2, &v[nedges*nedofs + f*nfdofs + 2*i]);
       }
    }
 }
 
-void ND_StatelessDofTransformation::InvTransformDual(const Array<int> & Fo,
-                                                     double *v) const
+void ND_DofTransformation::InvTransformDual(const Array<int> & Fo,
+                                            double *v) const
 {
    // Return immediately when no face DoFs are present
-   if (nfdofs < 2) { return; }
+   if (IsIdentity()) { return; }
 
    MFEM_VERIFY(Fo.Size() >= nfaces,
                "Face orientation array is shorter than the number of faces in "
-               "ND_StatelessDofTransformation");
+               "ND_DofTransformation");
 
    double data[2];
    Vector v2(data, 2);
+   DenseMatrix T2;
 
    // Transform face DoFs
    for (int f=0; f<nfaces; f++)
@@ -312,7 +308,8 @@ void ND_StatelessDofTransformation::InvTransformDual(const Array<int> & Fo,
       for (int i=0; i<nfdofs/2; i++)
       {
          v2 = &v[nedges*nedofs + f*nfdofs + 2*i];
-         T(Fo[f]).MultTranspose(v2, &v[nedges*nedofs + f*nfdofs + 2*i]);
+         T2.UseExternalData(const_cast<double *>(T.GetData(Fo[f])), 2, 2);
+         T2.MultTranspose(v2, &v[nedges*nedofs + f*nfdofs + 2*i]);
       }
    }
 }

fem/doftrans.hpp

@@ -80,6 +80,9 @@ public:
    inline int Width() const { return size_; }
    inline int NumCols() const { return size_; }
 
+   /// If the DofTransformation performs no transformation
+   virtual bool IsIdentity() const = 0;
+
    /** Transform local DoFs to align with the global DoFs. For example, this
        transformation can be used to map the local vector computed by
        FiniteElement::Project() to the transformed vector stored within a
@@ -115,6 +118,8 @@ public:
    inline void InvTransformDual(const Array<int> & face_orientation,
                                 Vector &v) const
    { InvTransformDual(face_orientation, v.GetData()); }
+
+   virtual ~StatelessDofTransformation() = default;
 };
 
 /** The DofTransformation class is an extension of the
@@ -133,35 +138,76 @@ public:
     transferring finite element degrees of freedom between different meshes.
     For examples of its use see the TransferMap used by the SubMesh class.
    */
-class DofTransformation : virtual public StatelessDofTransformation
+class DofTransformation
 {
 protected:
-   Array<int> Fo;
-
-   DofTransformation(int size)
-      : StatelessDofTransformation(size) {}
+   Array<int> Fo_;
+   const StatelessDofTransformation * dof_trans_;
+   int vdim_;
+   int ordering_;
 
 public:
+   /** @brief Default constructor which requires that SetDofTransformation be
+       called before use. */
+   DofTransformation(int vdim = 1, int ordering = 0)
+      : dof_trans_(NULL)
+      , vdim_(vdim)
+      , ordering_(ordering)
+   {}
+
+   /// Constructor with a known StatelessDofTransformation
+   DofTransformation(const StatelessDofTransformation & dof_trans,
+                     int vdim = 1, int ordering = 0)
+      : dof_trans_(&dof_trans)
+      , vdim_(vdim)
+      , ordering_(ordering)
+   {}
 
    /** @brief Configure the transformation using face orientations for the
        current element. */
    /// The face_orientation array can be obtained from Mesh::GetElementFaces.
-   inline void SetFaceOrientations(const Array<int> & face_orientation)
-   { Fo = face_orientation; }
+   inline void SetFaceOrientations(const Array<int> & Fo)
+   { Fo_ = Fo; }
 
-   inline const Array<int> & GetFaceOrientations() const { return Fo; }
+   /// Return the face orientations for the current element
+   inline const Array<int> & GetFaceOrientations() const { return Fo_; }
 
-   using StatelessDofTransformation::TransformPrimal;
-   using StatelessDofTransformation::InvTransformPrimal;
-   using StatelessDofTransformation::TransformDual;
-   using StatelessDofTransformation::InvTransformDual;
+   /// Set or change the nested StatelessDofTransformation object
+   inline void SetDofTransformation(const StatelessDofTransformation & dof_trans)
+   {
+      dof_trans_ = &dof_trans;
+   }
+   inline void SetDofTransformation(const StatelessDofTransformation * dof_trans)
+   {
+      dof_trans_ = dof_trans;
+   }
+
+   /// Return the nested StatelessDofTransformation object
+   inline const StatelessDofTransformation * GetDofTransformation() const
+   { return dof_trans_; }
+
+   /// Set or change the vdim and ordering parameter
+   inline void SetVDim(int vdim = 1, int ordering = 0)
+   {
+      vdim_ = vdim;
+      ordering_ = ordering;
+   }
+
+   /// Return the current vdim value
+   inline int GetVDim() const { return vdim_; }
+
+   inline int Size() const { return dof_trans_->Size(); }
+   inline int Height() const { return dof_trans_->Height(); }
+   inline int NumRows() const { return dof_trans_->NumRows(); }
+   inline int Width() const { return dof_trans_->Width(); }
+   inline int NumCols() const { return dof_trans_->NumCols(); }
+   inline bool IsIdentity() const { return dof_trans_->IsIdentity(); }
 
    /** Transform local DoFs to align with the global DoFs. For example, this
        transformation can be used to map the local vector computed by
        FiniteElement::Project() to the transformed vector stored within a
        GridFunction object. */
-   inline void TransformPrimal(double *v) const
-   { TransformPrimal(Fo, v); }
+   void TransformPrimal(double *v) const;
    inline void TransformPrimal(Vector &v) const
    { TransformPrimal(v.GetData()); }
 
@@ -179,21 +225,18 @@ public:
        transform the vector obtained using GridFunction::GetSubVector before it
        can be used to compute a local interpolation.
    */
-   inline void InvTransformPrimal(double *v) const
-   { InvTransformPrimal(Fo, v); }
+   void InvTransformPrimal(double *v) const;
    inline void InvTransformPrimal(Vector &v) const
    { InvTransformPrimal(v.GetData()); }
 
    /** Transform dual DoFs as computed by a LinearFormIntegrator before summing
        into a LinearForm object. */
-   inline void TransformDual(double *v) const
-   { TransformDual(Fo, v); }
+   void TransformDual(double *v) const;
    inline void TransformDual(Vector &v) const
    { TransformDual(v.GetData()); }
 
    /** Inverse Transform dual DoFs */
-   inline void InvTransformDual(double *v) const
-   { InvTransformDual(Fo, v); }
+   void InvTransformDual(double *v) const;
    inline void InvTransformDual(Vector &v) const
    { InvTransformDual(v.GetData()); }
 
@@ -225,8 +268,6 @@ public:
          TransformDual(V.GetColumn(c));
       }
    }
-
-   virtual ~DofTransformation() = default;
 };
 
 /** Transform a matrix of DoFs entries from different finite element spaces as
@@ -245,145 +286,6 @@ void TransformDual(const DofTransformation *ran_dof_trans,
                    const DofTransformation *dom_dof_trans,
                    DenseMatrix &elmat);
 
-/** The StatelessVDofTransformation class implements a nested transformation
-    where an arbitrary StatelessDofTransformation is replicated with a
-    vdim >= 1.
-*/
-class StatelessVDofTransformation : virtual public StatelessDofTransformation
-{
-protected:
-   int vdim_;
-   int ordering_;
-   StatelessDofTransformation * sdoftrans_;
-
-public:
-   /** @brief Default constructor which requires that SetDofTransformation be
-       called before use. */
-   StatelessVDofTransformation(int vdim = 1, int ordering = 0)
-      : StatelessDofTransformation(0)
-      , vdim_(vdim)
-      , ordering_(ordering)
-      , sdoftrans_(NULL)
-   {}
-
-   /// Constructor with a known StatelessDofTransformation
-   StatelessVDofTransformation(StatelessDofTransformation & doftrans,
-                               int vdim = 1,
-                               int ordering = 0)
-      : StatelessDofTransformation(vdim * doftrans.Size())
-      , vdim_(vdim)
-      , ordering_(ordering)
-      , sdoftrans_(&doftrans)
-   {}
-
-   /// Set or change the vdim parameter
-   inline void SetVDim(int vdim)
-   {
-      vdim_ = vdim;
-      if (sdoftrans_)
-      {
-         size_ = vdim_ * sdoftrans_->Size();
-      }
-   }
-
-   /// Return the current vdim value
-   inline int GetVDim() const { return vdim_; }
-
-   /// Set or change the nested StatelessDofTransformation object
-   inline void SetDofTransformation(StatelessDofTransformation & doftrans)
-   {
-      size_ = vdim_ * doftrans.Size();
-      sdoftrans_ = &doftrans;
-   }
-
-   /// Return the nested StatelessDofTransformation object
-   inline StatelessDofTransformation * GetDofTransformation() const
-   { return sdoftrans_; }
-
-   using StatelessDofTransformation::TransformPrimal;
-   using StatelessDofTransformation::InvTransformPrimal;
-   using StatelessDofTransformation::TransformDual;
-   using StatelessDofTransformation::InvTransformDual;
-
-   /** Specializations of these base class methods which account for the vdim
-       and ordering of the full set of DoFs.
-   */
-   void TransformPrimal(const Array<int> & face_ori, double *v) const;
-   void InvTransformPrimal(const Array<int> & face_ori, double *v) const;
-   void TransformDual(const Array<int> & face_ori, double *v) const;
-   void InvTransformDual(const Array<int> & face_ori, double *v) const;
-};
-
-/** The VDofTransformation class implements a nested transformation where an
-    arbitrary DofTransformation is replicated with a vdim >= 1.
-*/
-class VDofTransformation : public StatelessVDofTransformation,
-   public DofTransformation
-{
-protected:
-   DofTransformation * doftrans_;
-
-public:
-   /** @brief Default constructor which requires that SetDofTransformation be
-       called before use. */
-   VDofTransformation(int vdim = 1, int ordering = 0)
-      : StatelessDofTransformation(0)
-      , StatelessVDofTransformation(vdim, ordering)
-      , DofTransformation(0)
-      , doftrans_(NULL)
-   {}
-
-   /// Constructor with a known DofTransformation
-   /// @note The face orientations in @a doftrans will be copied into the
-   /// new VDofTransformation object.
-   VDofTransformation(DofTransformation & doftrans, int vdim = 1,
-                      int ordering = 0)
-      : StatelessDofTransformation(vdim * doftrans.Size())
-      , StatelessVDofTransformation(doftrans, vdim, ordering)
-      , DofTransformation(vdim * doftrans.Size())
-      , doftrans_(&doftrans)
-   {
-      DofTransformation::SetFaceOrientations(doftrans.GetFaceOrientations());
-   }
-
-   using StatelessVDofTransformation::SetDofTransformation;
-
-   /// Set or change the nested DofTransformation object
-   /// @note The face orientations in @a doftrans will be copied into the
-   /// VDofTransformation object.
-   void SetDofTransformation(DofTransformation & doftrans)
-   {
-      doftrans_ = &doftrans;
-      StatelessVDofTransformation::SetDofTransformation(doftrans);
-      DofTransformation::SetFaceOrientations(doftrans.GetFaceOrientations());
-   }
-
-   /// Return the nested DofTransformation object
-   inline DofTransformation * GetDofTransformation() const { return doftrans_; }
-
-   /// Set new face orientations in both the VDofTransformation and the
-   /// DofTransformation contained within (if there is one).
-   inline void SetFaceOrientations(const Array<int> & face_orientation)
-   {
-      DofTransformation::SetFaceOrientations(face_orientation);
-      if (doftrans_) { doftrans_->SetFaceOrientations(face_orientation); }
-   }
-
-   using DofTransformation::TransformPrimal;
-   using DofTransformation::InvTransformPrimal;
-   using DofTransformation::TransformDual;
-   using DofTransformation::InvTransformDual;
-
-   inline void TransformPrimal(double *v) const
-   { TransformPrimal(Fo, v); }
-   inline void InvTransformPrimal(double *v) const
-   { InvTransformPrimal(Fo, v); }
-   inline void TransformDual(double *v) const
-   { TransformDual(Fo, v); }
-   inline void InvTransformDual(double *v) const
-   { InvTransformDual(Fo, v); }
-};
-
 /** Abstract base class for high-order Nedelec spaces on elements with
     triangular faces.
 
@@ -396,7 +298,7 @@ public:
     be accessed as DenseMatrices using the GetFaceTransform() and
     GetFaceInverseTransform() methods.
 */
-class ND_StatelessDofTransformation : virtual public StatelessDofTransformation
+class ND_DofTransformation : public StatelessDofTransformation
 {
 private:
    static const double T_data[24];
@@ -410,8 +312,7 @@ protected:
    const int nedges; // number of edges per element
    const int nfaces; // number of triangular faces per element
 
-   ND_StatelessDofTransformation(int size, int order,
-                                 int num_edges, int num_tri_faces);
+   ND_DofTransformation(int size, int order, int num_edges, int num_tri_faces);
 
 public:
    // Return the 2x2 transformation operator for the given face orientation
@@ -421,116 +322,41 @@ public:
    static const DenseMatrix & GetFaceInverseTransform(int ori)
    { return TInv(ori); }
 
-   void TransformPrimal(const Array<int> & face_orientation,
-                        double *v) const;
+   bool IsIdentity() const override { return nfdofs < 2; }
 
-   void InvTransformPrimal(const Array<int> & face_orientation,
-                           double *v) const;
-
-   void TransformDual(const Array<int> & face_orientation,
-                      double *v) const;
-
-   void InvTransformDual(const Array<int> & face_orientation,
-                         double *v) const;
+   void TransformPrimal(const Array<int> & Fo, double *v) const override;
+   void InvTransformPrimal(const Array<int> & Fo, double *v) const override;
+   void TransformDual(const Array<int> & Fo, double *v) const override;
+   void InvTransformDual(const Array<int> & Fo, double *v) const override;
 };
 
 /// Stateless DoF transformation implementation for the Nedelec basis on
 /// triangles
-class ND_TriStatelessDofTransformation : public ND_StatelessDofTransformation
-{
-public:
-   ND_TriStatelessDofTransformation(int order)
-      : StatelessDofTransformation(order*(order + 2))
-      , ND_StatelessDofTransformation(order*(order + 2), order, 3, 1)
-   {}
-};
-
-/// DoF transformation implementation for the Nedelec basis on triangles
-class ND_TriDofTransformation : public DofTransformation,
-   public ND_TriStatelessDofTransformation
+class ND_TriDofTransformation : public ND_DofTransformation
 {
 public:
    ND_TriDofTransformation(int order)
-      : StatelessDofTransformation(order*(order + 2))
-      , DofTransformation(order*(order + 2))
-      , ND_TriStatelessDofTransformation(order)
-   {}
-
-   using DofTransformation::TransformPrimal;
-   using DofTransformation::InvTransformPrimal;
-   using DofTransformation::TransformDual;
-   using DofTransformation::InvTransformDual;
-
-   using ND_TriStatelessDofTransformation::TransformPrimal;
-   using ND_TriStatelessDofTransformation::InvTransformPrimal;
-   using ND_TriStatelessDofTransformation::TransformDual;
-   using ND_TriStatelessDofTransformation::InvTransformDual;
-};
-
-/// DoF transformation implementation for the Nedelec basis on tetrahedra
-class ND_TetStatelessDofTransformation : public ND_StatelessDofTransformation
-{
-public:
-   ND_TetStatelessDofTransformation(int order)
-      : StatelessDofTransformation(order*(order + 2)*(order + 3)/2)
-      , ND_StatelessDofTransformation(order*(order + 2)*(order + 3)/2, order,
-                                      6, 4)
+      : ND_DofTransformation(order*(order + 2), order, 3, 1)
    {}
 };
 
 /// DoF transformation implementation for the Nedelec basis on tetrahedra
-class ND_TetDofTransformation : public DofTransformation,
-   public ND_TetStatelessDofTransformation
+class ND_TetDofTransformation : public ND_DofTransformation
 {
 public:
    ND_TetDofTransformation(int order)
-      : StatelessDofTransformation(order*(order + 2)*(order + 3)/2)
-      , DofTransformation(order*(order + 2)*(order + 3)/2)
-      , ND_TetStatelessDofTransformation(order)
-   {}
-
-   using DofTransformation::TransformPrimal;
-   using DofTransformation::InvTransformPrimal;
-   using DofTransformation::TransformDual;
-   using DofTransformation::InvTransformDual;
-
-   using ND_TetStatelessDofTransformation::TransformPrimal;
-   using ND_TetStatelessDofTransformation::InvTransformPrimal;
-   using ND_TetStatelessDofTransformation::TransformDual;
-   using ND_TetStatelessDofTransformation::InvTransformDual;
-};
-
-/// DoF transformation implementation for the Nedelec basis on wedge elements
-class ND_WedgeStatelessDofTransformation : public ND_StatelessDofTransformation
-{
-public:
-   ND_WedgeStatelessDofTransformation(int order)
-      : StatelessDofTransformation(3 * order * ((order + 1) * (order + 2))/2)
-      , ND_StatelessDofTransformation(3 * order * ((order + 1) * (order + 2))/2,
-                                      order, 9, 2)
+      : ND_DofTransformation(order*(order + 2)*(order + 3)/2, order, 6, 4)
    {}
 };
 
 /// DoF transformation implementation for the Nedelec basis on wedge elements
-class ND_WedgeDofTransformation : public DofTransformation,
-   public ND_WedgeStatelessDofTransformation
+class ND_WedgeDofTransformation : public ND_DofTransformation
 {
 public:
    ND_WedgeDofTransformation(int order)
-      : StatelessDofTransformation(3 * order * ((order + 1) * (order + 2))/2)
-      , DofTransformation(3 * order * ((order + 1) * (order + 2))/2)
-      , ND_WedgeStatelessDofTransformation(order)
+      : ND_DofTransformation(3 * order * ((order + 1) * (order + 2))/2,
+                             order, 9, 2)
    {}
-
-   using DofTransformation::TransformPrimal;
-   using DofTransformation::InvTransformPrimal;
-   using DofTransformation::TransformDual;
-   using DofTransformation::InvTransformDual;
-
-   using ND_WedgeStatelessDofTransformation::TransformPrimal;
-   using ND_WedgeStatelessDofTransformation::InvTransformPrimal;
-   using ND_WedgeStatelessDofTransformation::TransformDual;
-   using ND_WedgeStatelessDofTransformation::InvTransformDual;
 };
 
 } // namespace mfem

fem/fe/fe_base.hpp

@@ -596,7 +596,7 @@ public:
    /** @brief Return a DoF transformation object for this particular type of
        basis.
    */
-   virtual StatelessDofTransformation * GetDofTransformation() const
+   virtual const StatelessDofTransformation *GetDofTransformation() const
    { return NULL; }
 
    /// Deconstruct the FiniteElement

fem/fe/fe_fixed_order.cpp

@@ -6031,7 +6031,7 @@ void RT0PyrFiniteElement::CalcVShape(const IntegrationPoint &ip,
    shape(1,2) = z;
 
    shape(2,0) = x * (2.0 - z) * ozi;
-   shape(2,1) = - y * z * ozi;;
+   shape(2,1) = - y * z * ozi;
    shape(2,2) = z;
 
    shape(3,0) = - x * z * ozi;

fem/fe/fe_nd.hpp

@@ -179,7 +179,7 @@ class ND_TetrahedronElement : public VectorFiniteElement
    Array<int> dof2tk;
    DenseMatrixInverse Ti;
 
-   mutable ND_TetStatelessDofTransformation doftrans;
+   ND_TetDofTransformation doftrans;
 
 public:
    /// Construct the ND_TetrahedronElement of order @a p
@@ -201,7 +201,7 @@ public:
                                   ElementTransformation &Trans,
                                   DenseMatrix &I) const
    { LocalInterpolation_ND(CheckVectorFE(fe), tk, dof2tk, Trans, I); }
-   virtual StatelessDofTransformation * GetDofTransformation() const
+   virtual const StatelessDofTransformation *GetDofTransformation() const
    { return &doftrans; }
    using FiniteElement::Project;
    virtual void Project(VectorCoefficient &vc,
@@ -242,7 +242,7 @@ class ND_TriangleElement : public VectorFiniteElement
    Array<int> dof2tk;
    DenseMatrixInverse Ti;
 
-   mutable ND_TriStatelessDofTransformation doftrans;
+   ND_TriDofTransformation doftrans;
 
 public:
    /// Construct the ND_TriangleElement of order @a p
@@ -264,7 +264,7 @@ public:
                                   ElementTransformation &Trans,
                                   DenseMatrix &I) const
    { LocalInterpolation_ND(CheckVectorFE(fe), tk, dof2tk, Trans, I); }
-   virtual StatelessDofTransformation * GetDofTransformation() const
+   virtual const StatelessDofTransformation *GetDofTransformation() const
    { return &doftrans; }
    using FiniteElement::Project;
    virtual void Project(VectorCoefficient &vc,
@@ -346,7 +346,7 @@ private:
 #endif
    Array<int> dof2tk, t_dof, s_dof;
 
-   mutable ND_WedgeStatelessDofTransformation doftrans;
+   ND_WedgeDofTransformation doftrans;
 
    H1_TriangleElement H1TriangleFE;
    ND_TriangleElement NDTriangleFE;
@@ -379,7 +379,7 @@ public:
                                   DenseMatrix &I) const
    { LocalInterpolation_ND(CheckVectorFE(fe), tk, dof2tk, Trans, I); }
 
-   virtual StatelessDofTransformation * GetDofTransformation() const
+   virtual const StatelessDofTransformation *GetDofTransformation() const
    { return &doftrans; }
 
    using FiniteElement::Project;

fem/fe_coll.cpp

@@ -2896,7 +2896,7 @@ ND_FECollection::FiniteElementForGeometry(Geometry::Type GeomType) const
    }
 }
 
-StatelessDofTransformation *
+const StatelessDofTransformation *
 ND_FECollection::DofTransformationForGeometry(Geometry::Type GeomType) const
 {
    if (!Geometry::IsTensorProduct(GeomType) && this->GetOrder() > 1)

fem/fe_coll.hpp

@@ -63,7 +63,7 @@ public:
    /** @brief Returns a DoF transformation object compatible with this basis
        and geometry type.
    */
-   virtual StatelessDofTransformation *
+   virtual const StatelessDofTransformation *
    DofTransformationForGeometry(Geometry::Type GeomType) const
    { return NULL; }
 
@@ -483,7 +483,7 @@ public:
    int DofForGeometry(Geometry::Type GeomType) const override
    { return ND_dof[GeomType]; }
 
-   StatelessDofTransformation *
+   const StatelessDofTransformation *
    DofTransformationForGeometry(Geometry::Type GeomType) const override;
 
    const int *DofOrderForOrientation(Geometry::Type GeomType,

fem/fespace.cpp

@@ -63,7 +63,6 @@ FiniteElementSpace::FiniteElementSpace()
      elem_dof(NULL), elem_fos(NULL), bdr_elem_dof(NULL), bdr_elem_fos(NULL),
      face_dof(NULL),
      NURBSext(NULL), own_ext(false),
-     DoFTrans(0), VDoFTrans(vdim, ordering),
      cP_is_set(false),
      Th(Operator::ANY_TYPE),
      sequence(0), mesh_sequence(0), orders_changed(false), relaxed_hp(false)
@@ -72,7 +71,6 @@ FiniteElementSpace::FiniteElementSpace()
 FiniteElementSpace::FiniteElementSpace(const FiniteElementSpace &orig,
                                        Mesh *mesh_,
                                        const FiniteElementCollection *fec_)
-   : VDoFTrans(orig.vdim, orig.ordering)
 {
    mesh_ = mesh_ ? mesh_ : orig.mesh;
    fec_ = fec_ ? fec_ : orig.fec;
@@ -212,7 +210,7 @@ void FiniteElementSpace::GetVDofs(int vd, Array<int>& dofs, int ndofs_) const
    }
 }
 
-void FiniteElementSpace::DofsToVDofs (Array<int> &dofs, int ndofs_) const
+void FiniteElementSpace::DofsToVDofs(Array<int> &dofs, int ndofs_) const
 {
    if (vdim == 1) { return; }
    if (ndofs_ < 0) { ndofs_ = this->ndofs; }
@@ -264,7 +262,7 @@ int FiniteElementSpace::DofToVDof(int dof, int vd, int ndofs_) const
 }
 
 // static function
-void FiniteElementSpace::AdjustVDofs (Array<int> &vdofs)
+void FiniteElementSpace::AdjustVDofs(Array<int> &vdofs)
 {
    int n = vdofs.Size(), *vdof = vdofs;
    for (int i = 0; i < n; i++)
@@ -277,36 +275,36 @@ void FiniteElementSpace::AdjustVDofs (Array<int> &vdofs)
    }
 }
 
+void FiniteElementSpace::GetElementVDofs(int i, Array<int> &vdofs,
+                                         DofTransformation &doftrans) const
+{
+   GetElementDofs(i, vdofs, doftrans);
+   DofsToVDofs(vdofs);
+   doftrans.SetVDim(vdim, ordering);
+}
+
 DofTransformation *
 FiniteElementSpace::GetElementVDofs(int i, Array<int> &vdofs) const
 {
-   DofTransformation * doftrans = GetElementDofs(i, vdofs);
+   DoFTrans.SetDofTransformation(NULL);
+   GetElementVDofs(i, vdofs, DoFTrans);
+   return DoFTrans.GetDofTransformation() ? &DoFTrans : NULL;
+}
+
+void FiniteElementSpace::GetBdrElementVDofs(int i, Array<int> &vdofs,
+                                            DofTransformation &doftrans) const
+{
+   GetBdrElementDofs(i, vdofs, doftrans);
    DofsToVDofs(vdofs);
-   if (vdim == 1 || doftrans == NULL)
-   {
-      return doftrans;
-   }
-   else
-   {
-      VDoFTrans.SetDofTransformation(*doftrans);
-      return &VDoFTrans;
-   }
+   doftrans.SetVDim(vdim, ordering);
 }
 
 DofTransformation *
 FiniteElementSpace::GetBdrElementVDofs(int i, Array<int> &vdofs) const
 {
-   DofTransformation * doftrans = GetBdrElementDofs(i, vdofs);
-   DofsToVDofs(vdofs);
-   if (vdim == 1 || doftrans == NULL)
-   {
-      return doftrans;
-   }
-   else
-   {
-      VDoFTrans.SetDofTransformation(*doftrans);
-      return &VDoFTrans;
-   }
+   DoFTrans.SetDofTransformation(NULL);
+   GetBdrElementVDofs(i, vdofs, DoFTrans);
+   return DoFTrans.GetDofTransformation() ? &DoFTrans : NULL;
 }
 
 void FiniteElementSpace::GetPatchVDofs(int i, Array<int> &vdofs) const
@@ -777,9 +775,9 @@ FiniteElementSpace::H2L_GlobalRestrictionMatrix (FiniteElementSpace *lfes)
    return R;
 }
 
-void FiniteElementSpace
-::AddDependencies(SparseMatrix& deps, Array<int>& master_dofs,
-                  Array<int>& slave_dofs, DenseMatrix& I, int skipfirst)
+void FiniteElementSpace::AddDependencies(
+   SparseMatrix& deps, Array<int>& master_dofs, Array<int>& slave_dofs,
+   DenseMatrix& I, int skipfirst)
 {
    for (int i = skipfirst; i < slave_dofs.Size(); i++)
    {
@@ -802,11 +800,9 @@ void FiniteElementSpace
    }
 }
 
-void FiniteElementSpace
-::AddEdgeFaceDependencies(SparseMatrix &deps, Array<int> &master_dofs,
-                          const FiniteElement *master_fe,
-                          Array<int> &slave_dofs, int slave_face,
-                          const DenseMatrix *pm) const
+void FiniteElementSpace::AddEdgeFaceDependencies(
+   SparseMatrix &deps, Array<int> &master_dofs, const FiniteElement *master_fe,
+   Array<int> &slave_dofs, int slave_face, const DenseMatrix *pm) const
 {
    // In variable-order spaces in 3D, we need to only constrain interior face
    // DOFs (this is done one level up), since edge dependencies can be more
@@ -1533,12 +1529,12 @@ SparseMatrix* FiniteElementSpace::RefinementMatrix(int old_ndofs,
                                 localP);
 }
 
-FiniteElementSpace::RefinementOperator::RefinementOperator
-(const FiniteElementSpace* fespace, Table* old_elem_dof, Table* old_elem_fos,
- int old_ndofs)
-   : fespace(fespace)
-   , old_elem_dof(old_elem_dof)
-   , old_elem_fos(old_elem_fos)
+FiniteElementSpace::RefinementOperator::RefinementOperator(
+   const FiniteElementSpace* fespace, Table* old_elem_dof, Table* old_elem_fos,
+   int old_ndofs)
+   : fespace(fespace),
+     old_elem_dof(old_elem_dof),
+     old_elem_fos(old_elem_fos)
 {
    MFEM_VERIFY(fespace->GetNE() >= old_elem_dof->Size(),
                "Previous mesh is not coarser.");
@@ -1553,7 +1549,7 @@ FiniteElementSpace::RefinementOperator::RefinementOperator
       fespace->GetLocalRefinementMatrices(elem_geoms[i], localP[elem_geoms[i]]);
    }
 
-   ConstructDoFTrans();
+   ConstructDoFTransArray();
 }
 
 FiniteElementSpace::RefinementOperator::RefinementOperator(
@@ -1578,59 +1574,58 @@ FiniteElementSpace::RefinementOperator::RefinementOperator(
       old_elem_fos = new Table(*coarse_fes->GetElementToFaceOrientationTable());
    }
 
-   ConstructDoFTrans();
+   ConstructDoFTransArray();
 }
 
 FiniteElementSpace::RefinementOperator::~RefinementOperator()
 {
    delete old_elem_dof;
    delete old_elem_fos;
-   for (int i=0; i<old_DoFTrans.Size(); i++)
+   for (int i=0; i<old_DoFTransArray.Size(); i++)
    {
-      delete old_DoFTrans[i];
+      delete old_DoFTransArray[i];
    }
 }
 
-void FiniteElementSpace::RefinementOperator
-::ConstructDoFTrans()
+void FiniteElementSpace::RefinementOperator::ConstructDoFTransArray()
 {
-   old_DoFTrans.SetSize(Geometry::NUM_GEOMETRIES);
-   for (int i=0; i<old_DoFTrans.Size(); i++)
+   old_DoFTransArray.SetSize(Geometry::NUM_GEOMETRIES);
+   for (int i=0; i<old_DoFTransArray.Size(); i++)
    {
-      old_DoFTrans[i] = NULL;
+      old_DoFTransArray[i] = NULL;
    }
 
    const FiniteElementCollection *fec_ref = fespace->FEColl();
    if (dynamic_cast<const ND_FECollection*>(fec_ref))
    {
-      const FiniteElement * nd_tri =
+      const FiniteElement *nd_tri =
          fec_ref->FiniteElementForGeometry(Geometry::TRIANGLE);
       if (nd_tri)
       {
-         old_DoFTrans[Geometry::TRIANGLE] =
+         old_DoFTransArray[Geometry::TRIANGLE] =
             new ND_TriDofTransformation(nd_tri->GetOrder());
       }
 
-      const FiniteElement * nd_tet =
+      const FiniteElement *nd_tet =
          fec_ref->FiniteElementForGeometry(Geometry::TETRAHEDRON);
       if (nd_tet)
       {
-         old_DoFTrans[Geometry::TETRAHEDRON] =
+         old_DoFTransArray[Geometry::TETRAHEDRON] =
             new ND_TetDofTransformation(nd_tet->GetOrder());
       }
 
-      const FiniteElement * nd_pri =
+      const FiniteElement *nd_pri =
          fec_ref->FiniteElementForGeometry(Geometry::PRISM);
       if (nd_pri)
       {
-         old_DoFTrans[Geometry::PRISM] =
+         old_DoFTransArray[Geometry::PRISM] =
             new ND_WedgeDofTransformation(nd_pri->GetOrder());
       }
    }
 }
 
-void FiniteElementSpace::RefinementOperator
-::Mult(const Vector &x, Vector &y) const
+void FiniteElementSpace::RefinementOperator::Mult(const Vector &x,
+                                                  Vector &y) const
 {
    Mesh* mesh_ref = fespace->GetMesh();
    const CoarseFineTransformations &trans_ref =
@@ -1662,6 +1657,7 @@ void FiniteElementSpace::RefinementOperator
             fespace->DofsToVDofs(vd, vdofs);
             old_dofs.Copy(old_vdofs);
             fespace->DofsToVDofs(vd, old_vdofs, old_ndofs);
+
             x.GetSubVector(old_vdofs, subX);
             lP.Mult(subX, subY);
             y.SetSubVector(vdofs, subY);
@@ -1670,40 +1666,30 @@ void FiniteElementSpace::RefinementOperator
       else
       {
          old_elem_fos->GetRow(emb.parent, old_Fo);
-         old_DoFTrans[geom]->SetFaceOrientations(old_Fo);
-
-         DofTransformation *new_doftrans = NULL;
-         VDofTransformation *vdoftrans =
-            dynamic_cast<VDofTransformation*>(doftrans);
-         if (vdoftrans)
-         {
-            new_doftrans = doftrans;
-            doftrans = vdoftrans->GetDofTransformation();
-         }
+         old_DoFTrans.SetDofTransformation(*old_DoFTransArray[geom]);
+         old_DoFTrans.SetFaceOrientations(old_Fo);
 
+         doftrans->SetVDim();
          for (int vd = 0; vd < rvdim; vd++)
          {
             dofs.Copy(vdofs);
             fespace->DofsToVDofs(vd, vdofs);
             old_dofs.Copy(old_vdofs);
             fespace->DofsToVDofs(vd, old_vdofs, old_ndofs);
+
             x.GetSubVector(old_vdofs, subX);
-            old_DoFTrans[geom]->InvTransformPrimal(subX);
+            old_DoFTrans.InvTransformPrimal(subX);
             lP.Mult(subX, subY);
             doftrans->TransformPrimal(subY);
             y.SetSubVector(vdofs, subY);
          }
-
-         if (vdoftrans)
-         {
-            doftrans = new_doftrans;
-         }
+         doftrans->SetVDim(rvdim, fespace->GetOrdering());
       }
    }
 }
 
-void FiniteElementSpace::RefinementOperator
-::MultTranspose(const Vector &x, Vector &y) const
+void FiniteElementSpace::RefinementOperator::MultTranspose(const Vector &x,
+                                                           Vector &y) const
 {
    y = 0.0;
 
@@ -1727,7 +1713,7 @@ void FiniteElementSpace::RefinementOperator
       const Geometry::Type geom = mesh_ref->GetElementBaseGeometry(k);
       const DenseMatrix &lP = localP[geom](emb.matrix);
 
-      DofTransformation * doftrans = fespace->GetElementDofs(k, f_dofs);
+      DofTransformation *doftrans = fespace->GetElementDofs(k, f_dofs);
       old_elem_dof->GetRow(emb.parent, c_dofs);
 
       if (!doftrans)
@@ -1742,7 +1728,6 @@ void FiniteElementSpace::RefinementOperator
             fespace->DofsToVDofs(vd, c_vdofs, old_ndofs);
 
             x.GetSubVector(f_vdofs, subX);
-
             for (int p = 0; p < f_dofs.Size(); ++p)
             {
                if (processed[DecodeDof(f_dofs[p])])
@@ -1750,7 +1735,6 @@ void FiniteElementSpace::RefinementOperator
                   subX[p] = 0.0;
                }
             }
-
             lP.MultTranspose(subX, subY);
             y.AddElementVector(c_vdofs, subY);
          }
@@ -1760,17 +1744,10 @@ void FiniteElementSpace::RefinementOperator
          subYt.SetSize(lP.Width());
 
          old_elem_fos->GetRow(emb.parent, old_Fo);
-         old_DoFTrans[geom]->SetFaceOrientations(old_Fo);
-
-         DofTransformation *new_doftrans = NULL;
-         VDofTransformation *vdoftrans =
-            dynamic_cast<VDofTransformation*>(doftrans);
-         if (vdoftrans)
-         {
-            new_doftrans = doftrans;
-            doftrans = vdoftrans->GetDofTransformation();
-         }
+         old_DoFTrans.SetDofTransformation(*old_DoFTransArray[geom]);
+         old_DoFTrans.SetFaceOrientations(old_Fo);
 
+         doftrans->SetVDim();
          for (int vd = 0; vd < rvdim; vd++)
          {
             f_dofs.Copy(f_vdofs);
@@ -1787,16 +1764,11 @@ void FiniteElementSpace::RefinementOperator
                   subX[p] = 0.0;
                }
             }
-
             lP.MultTranspose(subX, subYt);
-            old_DoFTrans[geom]->TransformDual(subYt);
+            old_DoFTrans.TransformDual(subYt);
             y.AddElementVector(c_vdofs, subYt);
          }
-
-         if (vdoftrans)
-         {
-            doftrans = new_doftrans;
-         }
+         doftrans->SetVDim(rvdim, fespace->GetOrdering());
       }
 
       for (int p = 0; p < f_dofs.Size(); ++p)
@@ -2024,8 +1996,8 @@ FiniteElementSpace::DerefinementOperator::~DerefinementOperator()
    delete coarse_elem_dof;
 }
 
-void FiniteElementSpace::DerefinementOperator
-::Mult(const Vector &x, Vector &y) const
+void FiniteElementSpace::DerefinementOperator::Mult(const Vector &x,
+                                                    Vector &y) const
 {
    Array<int> c_vdofs, f_vdofs;
    Vector loc_x, loc_y;
@@ -2227,7 +2199,7 @@ void FiniteElementSpace::Constructor(Mesh *mesh_, NURBSExtension *NURBSext_,
       R_transpose.reset();
       cP_is_set = false;
 
-      ConstructDoFTrans();
+      ConstructDoFTransArray();
    }
    else
    {
@@ -2239,40 +2211,39 @@ void FiniteElementSpace::Constructor(Mesh *mesh_, NURBSExtension *NURBSext_,
    BuildElementToDofTable();
 }
 
-void FiniteElementSpace::ConstructDoFTrans()
+void FiniteElementSpace::ConstructDoFTransArray()
 {
-   DestroyDoFTrans();
+   DestroyDoFTransArray();
 
-   VDoFTrans.SetVDim(vdim);
-   DoFTrans.SetSize(Geometry::NUM_GEOMETRIES);
-   for (int i=0; i<DoFTrans.Size(); i++)
+   DoFTransArray.SetSize(Geometry::NUM_GEOMETRIES);
+   for (int i=0; i<DoFTransArray.Size(); i++)
    {
-      DoFTrans[i] = NULL;
+      DoFTransArray[i] = NULL;
    }
    if (mesh->Dimension() < 3) { return; }
    if (dynamic_cast<const ND_FECollection*>(fec))
    {
-      const FiniteElement * nd_tri =
+      const FiniteElement *nd_tri =
          fec->FiniteElementForGeometry(Geometry::TRIANGLE);
       if (nd_tri)
       {
-         DoFTrans[Geometry::TRIANGLE] =
+         DoFTransArray[Geometry::TRIANGLE] =
             new ND_TriDofTransformation(nd_tri->GetOrder());
       }
 
-      const FiniteElement * nd_tet =
+      const FiniteElement *nd_tet =
          fec->FiniteElementForGeometry(Geometry::TETRAHEDRON);
       if (nd_tet)
       {
-         DoFTrans[Geometry::TETRAHEDRON] =
+         DoFTransArray[Geometry::TETRAHEDRON] =
             new ND_TetDofTransformation(nd_tet->GetOrder());
       }
 
-      const FiniteElement * nd_pri =
+      const FiniteElement *nd_pri =
          fec->FiniteElementForGeometry(Geometry::PRISM);
       if (nd_pri)
       {
-         DoFTrans[Geometry::PRISM] =
+         DoFTransArray[Geometry::PRISM] =
             new ND_WedgeDofTransformation(nd_pri->GetOrder());
       }
    }
@@ -2324,7 +2295,7 @@ void FiniteElementSpace::BuildNURBSFaceToDofTable() const
    face_to_be = -1;
    for (int b = 0; b < GetNBE(); b++)
    {
-      int f = mesh->GetBdrElementEdgeIndex(b);
+      int f = mesh->GetBdrElementFaceIndex(b);
       face_to_be[f] = b;
    }
 
@@ -2476,7 +2447,7 @@ void FiniteElementSpace::Construct()
 
    ndofs = nvdofs + nedofs + nfdofs + nbdofs;
 
-   ConstructDoFTrans();
+   ConstructDoFTransArray();
 
    // record the current mesh sequence number to detect refinement etc.
    mesh_sequence = mesh->GetSequence();
@@ -2501,9 +2472,8 @@ int FiniteElementSpace::MinOrder(VarOrderBits bits)
    return 0;
 }
 
-void FiniteElementSpace
-::CalcEdgeFaceVarOrders(Array<VarOrderBits> &edge_orders,
-                        Array<VarOrderBits> &face_orders) const
+void FiniteElementSpace::CalcEdgeFaceVarOrders(
+   Array<VarOrderBits> &edge_orders, Array<VarOrderBits> &face_orders) const
 {
    MFEM_ASSERT(IsVariableOrder(), "");
    MFEM_ASSERT(Nonconforming(), "");
@@ -2727,8 +2697,8 @@ int FiniteElementSpace::GetNVariants(int entity, int index) const
 static const char* msg_orders_changed =
    "Element orders changed, you need to Update() the space first.";
 
-DofTransformation *
-FiniteElementSpace::GetElementDofs(int elem, Array<int> &dofs) const
+void FiniteElementSpace::GetElementDofs(int elem, Array<int> &dofs,
+                                        DofTransformation &doftrans) const
 {
    MFEM_VERIFY(!orders_changed, msg_orders_changed);
 
@@ -2736,13 +2706,16 @@ FiniteElementSpace::GetElementDofs(int elem, Array<int> &dofs) const
    {
       elem_dof->GetRow(elem, dofs);
 
-      if (DoFTrans[mesh->GetElementBaseGeometry(elem)])
+      if (DoFTransArray[mesh->GetElementBaseGeometry(elem)])
       {
          Array<int> Fo;
          elem_fos -> GetRow (elem, Fo);
-         DoFTrans[mesh->GetElementBaseGeometry(elem)]->SetFaceOrientations(Fo);
+         doftrans.SetDofTransformation(
+            *DoFTransArray[mesh->GetElementBaseGeometry(elem)]);
+         doftrans.SetFaceOrientations(Fo);
+         doftrans.SetVDim();
       }
-      return DoFTrans[mesh->GetElementBaseGeometry(elem)];
+      return;
    }
 
    Array<int> V, E, Eo, F, Fo; // TODO: LocalArray
@@ -2766,10 +2739,12 @@ FiniteElementSpace::GetElementDofs(int elem, Array<int> &dofs) const
       {
          nfd += fec->GetNumDof(mesh->GetFaceGeometry(F[i]), order);
       }
-      if (DoFTrans[mesh->GetElementBaseGeometry(elem)])
+      if (DoFTransArray[mesh->GetElementBaseGeometry(elem)])
       {
-         DoFTrans[mesh->GetElementBaseGeometry(elem)]
-         -> SetFaceOrientations(Fo);
+         doftrans.SetDofTransformation(
+            *DoFTransArray[mesh->GetElementBaseGeometry(elem)]);
+         doftrans.SetFaceOrientations(Fo);
+         doftrans.SetVDim();
       }
    }
 
@@ -2828,54 +2803,18 @@ FiniteElementSpace::GetElementDofs(int elem, Array<int> &dofs) const
          dofs.Append(bbase + j);
       }
    }
-   return DoFTrans[mesh->GetElementBaseGeometry(elem)];
 }
 
-void FiniteElementSpace::GetPatchDofs(int patch, Array<int> &dofs) const
+DofTransformation *FiniteElementSpace::GetElementDofs(int elem,
+                                                      Array<int> &dofs) const
 {
-   MFEM_ASSERT(NURBSext,
-               "FiniteElementSpace::GetPatchDofs needs a NURBSExtension");
-   NURBSext->GetPatchDofs(patch, dofs);
+   DoFTrans.SetDofTransformation(NULL);
+   GetElementDofs(elem, dofs, DoFTrans);
+   return DoFTrans.GetDofTransformation() ? &DoFTrans : NULL;
 }
 
-const FiniteElement *FiniteElementSpace::GetFE(int i) const
-{
-   if (i < 0 || i >= mesh->GetNE())
-   {
-      if (mesh->GetNE() == 0)
-      {
-         MFEM_ABORT("Empty MPI partitions are not permitted!");
-      }
-      MFEM_ABORT("Invalid element id:" << i << "; minimum allowed:" << 0 <<
-                 ", maximum allowed:" << mesh->GetNE()-1);
-   }
-
-   const FiniteElement *FE =
-      fec->GetFE(mesh->GetElementGeometry(i), GetElementOrderImpl(i));
-
-   if (NURBSext)
-   {
-      NURBSext->LoadFE(i, FE);
-   }
-   else
-   {
-#ifdef MFEM_DEBUG
-      // consistency check: fec->GetOrder() and FE->GetOrder() should return
-      // the same value (for standard, constant-order spaces)
-      if (!IsVariableOrder() && FE->GetDim() > 0)
-      {
-         MFEM_ASSERT(FE->GetOrder() == fec->GetOrder(),
-                     "internal error: " <<
-                     FE->GetOrder() << " != " << fec->GetOrder());
-      }
-#endif
-   }
-
-   return FE;
-}
-
-DofTransformation *
-FiniteElementSpace::GetBdrElementDofs(int bel, Array<int> &dofs) const
+void FiniteElementSpace::GetBdrElementDofs(int bel, Array<int> &dofs,
+                                           DofTransformation &doftrans) const
 {
    MFEM_VERIFY(!orders_changed, msg_orders_changed);
 
@@ -2883,17 +2822,19 @@ FiniteElementSpace::GetBdrElementDofs(int bel, Array<int> &dofs) const
    {
       bdr_elem_dof->GetRow(bel, dofs);
 
-      if (DoFTrans[mesh->GetBdrElementBaseGeometry(bel)])
+      if (DoFTransArray[mesh->GetBdrElementBaseGeometry(bel)])
       {
          Array<int> Fo;
          bdr_elem_fos -> GetRow (bel, Fo);
-         DoFTrans[mesh->GetBdrElementBaseGeometry(bel)]->
-         SetFaceOrientations(Fo);
+         doftrans.SetDofTransformation(
+            *DoFTransArray[mesh->GetBdrElementBaseGeometry(bel)]);
+         doftrans.SetFaceOrientations(Fo);
+         doftrans.SetVDim();
       }
-      return DoFTrans[mesh->GetBdrElementBaseGeometry(bel)];
+      return;
    }
 
-   Array<int> V, E, Eo, Fo; // TODO: LocalArray
+   Array<int> V, E, Eo; // TODO: LocalArray
    int F, oF;
 
    int dim = mesh->Dimension();
@@ -2917,11 +2858,14 @@ FiniteElementSpace::GetBdrElementDofs(int bel, Array<int> &dofs) const
    {
       mesh->GetBdrElementFace(bel, &F, &oF);
 
-      if (DoFTrans[mesh->GetBdrElementBaseGeometry(bel)])
+      if (DoFTransArray[mesh->GetBdrElementBaseGeometry(bel)])
       {
-         Fo.Append(oF);
-         DoFTrans[mesh->GetBdrElementBaseGeometry(bel)]->
-         SetFaceOrientations(Fo);
+         mfem::Array<int> Fo(1);
+         Fo[0] = oF;
+         doftrans.SetDofTransformation(
+            *DoFTransArray[mesh->GetBdrElementBaseGeometry(bel)]);
+         doftrans.SetFaceOrientations(Fo);
+         doftrans.SetVDim();
       }
    }
 
@@ -2963,8 +2907,14 @@ FiniteElementSpace::GetBdrElementDofs(int bel, Array<int> &dofs) const
          dofs.Append(EncodeDof(nvdofs + nedofs + fbase, ind[j]));
       }
    }
+}
 
-   return DoFTrans[mesh->GetBdrElementBaseGeometry(bel)];
+DofTransformation *FiniteElementSpace::GetBdrElementDofs(int bel,
+                                                         Array<int> &dofs) const
+{
+   DoFTrans.SetDofTransformation(NULL);
+   GetBdrElementDofs(bel, dofs, DoFTrans);
+   return DoFTrans.GetDofTransformation() ? &DoFTrans : NULL;
 }
 
 int FiniteElementSpace::GetFaceDofs(int face, Array<int> &dofs,
@@ -3134,18 +3084,6 @@ int FiniteElementSpace::GetNumElementInteriorDofs(int i) const
                          GetElementOrderImpl(i));
 }
 
-void FiniteElementSpace::GetEdgeInteriorDofs(int i, Array<int> &dofs) const
-{
-   MFEM_VERIFY(!IsVariableOrder(), "not implemented");
-
-   int ne = fec->DofForGeometry(Geometry::SEGMENT);
-   dofs.SetSize (ne);
-   for (int j = 0, k = nvdofs+i*ne; j < ne; j++, k++)
-   {
-      dofs[j] = k;
-   }
-}
-
 void FiniteElementSpace::GetFaceInteriorDofs(int i, Array<int> &dofs) const
 {
    MFEM_VERIFY(!IsVariableOrder(), "not implemented");
@@ -3170,6 +3108,61 @@ void FiniteElementSpace::GetFaceInteriorDofs(int i, Array<int> &dofs) const
    }
 }
 
+void FiniteElementSpace::GetEdgeInteriorDofs(int i, Array<int> &dofs) const
+{
+   MFEM_VERIFY(!IsVariableOrder(), "not implemented");
+
+   int ne = fec->DofForGeometry(Geometry::SEGMENT);
+   dofs.SetSize (ne);
+   for (int j = 0, k = nvdofs+i*ne; j < ne; j++, k++)
+   {
+      dofs[j] = k;
+   }
+}
+
+void FiniteElementSpace::GetPatchDofs(int patch, Array<int> &dofs) const
+{
+   MFEM_ASSERT(NURBSext,
+               "FiniteElementSpace::GetPatchDofs needs a NURBSExtension");
+   NURBSext->GetPatchDofs(patch, dofs);
+}
+
+const FiniteElement *FiniteElementSpace::GetFE(int i) const
+{
+   if (i < 0 || i >= mesh->GetNE())
+   {
+      if (mesh->GetNE() == 0)
+      {
+         MFEM_ABORT("Empty MPI partitions are not permitted!");
+      }
+      MFEM_ABORT("Invalid element id:" << i << "; minimum allowed:" << 0 <<
+                 ", maximum allowed:" << mesh->GetNE()-1);
+   }
+
+   const FiniteElement *FE =
+      fec->GetFE(mesh->GetElementGeometry(i), GetElementOrderImpl(i));
+
+   if (NURBSext)
+   {
+      NURBSext->LoadFE(i, FE);
+   }
+   else
+   {
+#ifdef MFEM_DEBUG
+      // consistency check: fec->GetOrder() and FE->GetOrder() should return
+      // the same value (for standard, constant-order spaces)
+      if (!IsVariableOrder() && FE->GetDim() > 0)
+      {
+         MFEM_ASSERT(FE->GetOrder() == fec->GetOrder(),
+                     "internal error: " <<
+                     FE->GetOrder() << " != " << fec->GetOrder());
+      }
+#endif
+   }
+
+   return FE;
+}
+
 const FiniteElement *FiniteElementSpace::GetBE(int i) const
 {
    int order = fec->GetOrder();
@@ -3242,8 +3235,8 @@ const FiniteElement *FiniteElementSpace::GetEdgeElement(int i,
    return fec->GetFE(Geometry::SEGMENT, eo);
 }
 
-const FiniteElement *FiniteElementSpace
-::GetTraceElement(int i, Geometry::Type geom_type) const
+const FiniteElement *FiniteElementSpace::GetTraceElement(
+   int i, Geometry::Type geom_type) const
 {
    return fec->TraceFiniteElementForGeometry(geom_type);
 }
@@ -3283,7 +3276,7 @@ void FiniteElementSpace::Destroy()
    }
    E2BFQ_array.SetSize(0);
 
-   DestroyDoFTrans();
+   DestroyDoFTransArray();
 
    dof_elem_array.DeleteAll();
    dof_ldof_array.DeleteAll();
@@ -3301,19 +3294,18 @@ void FiniteElementSpace::Destroy()
       delete bdr_elem_dof;
       delete bdr_elem_fos;
       delete face_dof;
-
       delete [] bdofs;
    }
    ceed::RemoveBasisAndRestriction(this);
 }
 
-void FiniteElementSpace::DestroyDoFTrans()
+void FiniteElementSpace::DestroyDoFTransArray()
 {
-   for (int i = 0; i < DoFTrans.Size(); i++)
+   for (int i = 0; i < DoFTransArray.Size(); i++)
    {
-      delete DoFTrans[i];
+      delete DoFTransArray[i];
    }
-   DoFTrans.SetSize(0);
+   DoFTransArray.SetSize(0);
 }
 
 void FiniteElementSpace::GetTransferOperator(

fem/fespace.hpp

@@ -271,8 +271,8 @@ protected:
    int own_ext;
    mutable Array<int> face_to_be; // NURBS FE space only
 
-   Array<DofTransformation*> DoFTrans;
-   mutable VDofTransformation VDoFTrans;
+   Array<StatelessDofTransformation *> DoFTransArray;
+   mutable DofTransformation DoFTrans;
 
    /** Matrix representing the prolongation from the global conforming dofs to
        a set of intermediate partially conforming dofs, e.g. the dofs associated
@@ -328,8 +328,8 @@ protected:
    void Construct();
    void Destroy();
 
-   void ConstructDoFTrans();
-   void DestroyDoFTrans();
+   void ConstructDoFTransArray();
+   void DestroyDoFTransArray();
 
    void BuildElementToDofTable() const;
    void BuildBdrElementToDofTable() const;
@@ -416,10 +416,10 @@ protected:
       Table* old_elem_dof; // Owned.
       Table* old_elem_fos; // Owned.
 
-      Array<DofTransformation*> old_DoFTrans;
-      mutable VDofTransformation old_VDoFTrans;
+      Array<StatelessDofTransformation*> old_DoFTransArray;
+      mutable DofTransformation old_DoFTrans;
 
-      void ConstructDoFTrans();
+      void ConstructDoFTransArray();
 
    public:
       /** Construct the operator based on the elem_dof table of the original
@@ -803,7 +803,16 @@ public:
    /// with triangular faces.
    ///
    /// @note The returned object should NOT be deleted by the caller.
-   virtual DofTransformation *GetElementDofs(int elem, Array<int> &dofs) const;
+   DofTransformation *GetElementDofs(int elem, Array<int> &dofs) const;
+
+   /// @brief The same as GetElementDofs(), but with a user-allocated
+   /// DofTransformation object. @a doftrans must be allocated in advance and
+   /// will be owned by the caller. The user can use the
+   /// DofTransformation::GetDofTransformation method on the returned
+   /// @a doftrans object to detect if the DofTransformation should actually be
+   /// used.
+   virtual void GetElementDofs(int elem, Array<int> &dofs,
+                               DofTransformation &doftrans) const;
 
    /// @brief Returns indices of degrees of freedom for boundary element 'bel'.
    /// The returned indices are offsets into an @ref ldof vector. See also
@@ -817,13 +826,16 @@ public:
    /// with triangular faces.
    ///
    /// @note The returned object should NOT be deleted by the caller.
-   virtual DofTransformation *GetBdrElementDofs(int bel,
-                                                Array<int> &dofs) const;
+   DofTransformation *GetBdrElementDofs(int bel, Array<int> &dofs) const;
 
-   /** @brief Returns indices of degrees of freedom for NURBS patch index
-    @a patch. Cartesian ordering is used, for the tensor-product degrees of
-    freedom. */
-   void GetPatchDofs(int patch, Array<int> &dofs) const;
+   /// @brief The same as GetBdrElementDofs(), but with a user-allocated
+   /// DofTransformation object. @a doftrans must be allocated in advance and
+   /// will be owned by the caller. The user can use the
+   /// DofTransformation::GetDofTransformation method on the returned
+   /// @a doftrans object to detect if the DofTransformation should actually be
+   /// used.
+   virtual void GetBdrElementDofs(int bel, Array<int> &dofs,
+                                  DofTransformation &doftrans) const;
 
    /// @brief Returns the indices of the degrees of freedom for the specified
    /// face, including the DOFs for the edges and the vertices of the face.
@@ -870,6 +882,13 @@ public:
    /// GetElementInteriorVDofs().
    void GetElementInteriorDofs(int i, Array<int> &dofs) const;
 
+   /// @brief Returns the number of degrees of freedom associated with the
+   /// interior of the specified element.
+   ///
+   /// See GetElementInteriorDofs() for more information or to obtain the
+   /// relevant indices.
+   int GetNumElementInteriorDofs(int i) const;
+
    /// @brief Returns the indices of the degrees of freedom for the interior
    /// of the specified face.
    ///
@@ -882,13 +901,6 @@ public:
    /// GetFaceInteriorVDofs().
    void GetFaceInteriorDofs(int i, Array<int> &dofs) const;
 
-   /// @brief Returns the number of degrees of freedom associated with the
-   /// interior of the specified element.
-   ///
-   /// See GetElementInteriorDofs() for more information or to obtain the
-   /// relevant indices.
-   int GetNumElementInteriorDofs(int i) const;
-
    /// @brief Returns the indices of the degrees of freedom for the interior
    /// of the specified edge.
    ///
@@ -897,6 +909,11 @@ public:
    void GetEdgeInteriorDofs(int i, Array<int> &dofs) const;
    ///@}
 
+   /** @brief Returns indices of degrees of freedom for NURBS patch index
+    @a patch. Cartesian ordering is used, for the tensor-product degrees of
+    freedom. */
+   void GetPatchDofs(int patch, Array<int> &dofs) const;
+
    /// @anchor dof2vdof @name DoF To VDoF Conversion methods
    /// These methods convert between local dof and local vector dof using the
    /// appropriate relationship based on the Ordering::Type defined in this
@@ -1023,6 +1040,15 @@ public:
    /// @note The returned object should NOT be deleted by the caller.
    DofTransformation *GetElementVDofs(int i, Array<int> &vdofs) const;
 
+   /// @brief The same as GetElementVDofs(), but with a user-allocated
+   /// DofTransformation object. @a doftrans must be allocated in advance and
+   /// will be owned by the caller. The user can use the
+   /// DofTransformation::GetDofTransformation method on the returned
+   /// @a doftrans object to detect if the DofTransformation should actually be
+   /// used.
+   void GetElementVDofs(int i, Array<int> &vdofs,
+                        DofTransformation &doftrans) const;
+
    /// @brief Returns indices of degrees of freedom for @a i'th boundary
    /// element.
    /// The returned indices are offsets into an @ref ldof vector with @b vdim
@@ -1038,6 +1064,15 @@ public:
    /// @note The returned object should NOT be deleted by the caller.
    DofTransformation *GetBdrElementVDofs(int i, Array<int> &vdofs) const;
 
+   /// @brief The same as GetBdrElementVDofs(), but with a user-allocated
+   /// DofTransformation object. @a doftrans must be allocated in advance and
+   /// will be owned by the caller. The user can use the
+   /// DofTransformation::GetDofTransformation method on the returned
+   /// @a doftrans object to detect if the DofTransformation should actually be
+   /// used.
+   void GetBdrElementVDofs(int i, Array<int> &vdofs,
+                           DofTransformation &doftrans) const;
+
    /// Returns indices of degrees of freedom in @a vdofs for NURBS patch @a i.
    void GetPatchVDofs(int i, Array<int> &vdofs) const;
 

fem/fmsconvert.cpp

@@ -31,13 +31,13 @@ FmsBasisTypeToMfemBasis(FmsBasisType b)
    switch (b)
    {
       case FMS_NODAL_GAUSS_OPEN:
-         retval = mfem::BasisType::GaussLegendre;;
+         retval = mfem::BasisType::GaussLegendre;
          break;
       case FMS_NODAL_GAUSS_CLOSED:
-         retval = mfem::BasisType::GaussLobatto;;
+         retval = mfem::BasisType::GaussLobatto;
          break;
       case FMS_POSITIVE:
-         retval = mfem::BasisType::Positive;;
+         retval = mfem::BasisType::Positive;
          break;
       case FMS_NODAL_UNIFORM_OPEN:
          retval = mfem::BasisType::OpenUniform;
@@ -1812,22 +1812,22 @@ MeshToFmsMesh(const Mesh *mmesh, FmsMesh *fmesh, FmsComponent *volume)
       switch (betype)
       {
          case Element::POINT:
-            bdr_eles[FMS_VERTEX].push_back(mmesh->GetBdrElementEdgeIndex(i));
+            bdr_eles[FMS_VERTEX].push_back(mmesh->GetBdrElementFaceIndex(i));
             break;
          case Element::SEGMENT:
-            bdr_eles[FMS_EDGE].push_back(mmesh->GetBdrElementEdgeIndex(i));
+            bdr_eles[FMS_EDGE].push_back(mmesh->GetBdrElementFaceIndex(i));
             break;
          case Element::TRIANGLE:
-            bdr_eles[FMS_TRIANGLE].push_back(mmesh->GetBdrElementEdgeIndex(i));
+            bdr_eles[FMS_TRIANGLE].push_back(mmesh->GetBdrElementFaceIndex(i));
             break;
          case Element::QUADRILATERAL:
-            bdr_eles[FMS_QUADRILATERAL].push_back(mmesh->GetBdrElementEdgeIndex(i));
+            bdr_eles[FMS_QUADRILATERAL].push_back(mmesh->GetBdrElementFaceIndex(i));
             break;
          case Element::TETRAHEDRON:
-            bdr_eles[FMS_TETRAHEDRON].push_back(mmesh->GetBdrElementEdgeIndex(i));
+            bdr_eles[FMS_TETRAHEDRON].push_back(mmesh->GetBdrElementFaceIndex(i));
             break;
          case Element::HEXAHEDRON:
-            bdr_eles[FMS_HEXAHEDRON].push_back(mmesh->GetBdrElementEdgeIndex(i));
+            bdr_eles[FMS_HEXAHEDRON].push_back(mmesh->GetBdrElementFaceIndex(i));
             break;
          default:
             MFEM_WARNING("Unsupported boundary element " << betype << " at boundary index "

fem/linearform.cpp

@@ -137,7 +137,7 @@ bool LinearForm::SupportsDevice() const
       // Make sure every boundary element corresponds to a boundary face
       for (int be = 0; be < fes->GetNBE(); ++be)
       {
-         const int f = mesh.GetBdrElementEdgeIndex(be);
+         const int f = mesh.GetBdrElementFaceIndex(be);
          const auto face_info = mesh.GetFaceInformation(f);
          if (!face_info.IsBoundary())
          {

fem/linearform_ext.cpp

@@ -148,7 +148,7 @@ void LinearFormExtension::Update()
       std::unordered_map<int,int> f_to_be;
       for (int i = 0; i < mesh.GetNBE(); ++i)
       {
-         const int f = mesh.GetBdrElementEdgeIndex(i);
+         const int f = mesh.GetBdrElementFaceIndex(i);
          f_to_be[f] = i;
       }
       MFEM_VERIFY(size_t(nf_bdr) == f_to_be.size(), "Incompatible sizes");

fem/lor/lor_ams.cpp

@@ -257,13 +257,13 @@ void BatchedLOR_AMS::FormCoordinateVectors(const Vector &X_vert)
    // vertices of the LOR mesh. The vertex coordinates are already computed in
    // E-vector format and passed in in X_vert.
    //
-   // In this function, we need to convert X_vert (which has the shape (dim,
+   // In this function, we need to convert X_vert (which has the shape (sdim,
    // ndof_per_el, nel_ho)) to T-DOF format.
    //
-   // We place the results in the vector xyz_tvec, which has shape (ntdofs, dim)
+   // We place the results in the vector xyz_tvec, which has shape (ntdofs, sdim)
    // and then make the hypre vectors x, y, and z point to subvectors.
    //
-   // In 2D, z is NULL.
+   // When the space dimension is 2, z is NULL.
 
    // Create the H1 vertex space and get the element restriction
    ElementDofOrdering ordering = ElementDofOrdering::LEXICOGRAPHIC;
@@ -275,17 +275,17 @@ void BatchedLOR_AMS::FormCoordinateVectors(const Vector &X_vert)
    const int nel_ho = vert_fes.GetNE();
    const int ndp1 = order + 1;
    const int ndof_per_el = static_cast<int>(pow(ndp1, dim));
-   const int sdim = dim;
+   const int sdim = vert_fes.GetMesh()->SpaceDimension();
    const int ntdofs = R->Height();
 
    const MemoryClass mc = GetHypreMemoryClass();
    bool dev = (mc == MemoryClass::DEVICE);
 
-   xyz_tvec = new Vector(ntdofs*dim);
+   xyz_tvec = new Vector(ntdofs*sdim);
 
-   auto xyz_tv = Reshape(HypreWrite(xyz_tvec->GetMemory()), ntdofs, dim);
+   auto xyz_tv = Reshape(HypreWrite(xyz_tvec->GetMemory()), ntdofs, sdim);
    const auto xyz_e =
-      Reshape(HypreRead(X_vert.GetMemory()), dim, ndof_per_el, nel_ho);
+      Reshape(HypreRead(X_vert.GetMemory()), sdim, ndof_per_el, nel_ho);
    const auto d_offsets = HypreRead(el_restr->Offsets().GetMemory());
    const auto d_indices = HypreRead(el_restr->Indices().GetMemory());
    const auto ltdof_ldof = HypreRead(R->GetMemoryJ());
@@ -309,7 +309,7 @@ void BatchedLOR_AMS::FormCoordinateVectors(const Vector &X_vert)
    x = new HypreParVector(vert_fes.GetComm(), glob_size, d_x_ptr, cols, dev);
    double *d_y_ptr = xyz_tv + 1*ntdofs;
    y = new HypreParVector(vert_fes.GetComm(), glob_size, d_y_ptr, cols, dev);
-   if (dim == 3)
+   if (sdim == 3)
    {
       double *d_z_ptr = xyz_tv + 2*ntdofs;
       z = new HypreParVector(vert_fes.GetComm(), glob_size, d_z_ptr, cols, dev);

fem/lor/lor_batched.cpp

@@ -77,6 +77,7 @@ void BatchedLORAssembly::FormLORVertexCoordinates(FiniteElementSpace &fes_ho,
 
    // Get nodal points at the LOR vertices
    const int dim = mesh_ho.Dimension();
+   const int sdim = mesh_ho.SpaceDimension();
    const int nel_ho = mesh_ho.GetNE();
    const int order = fes_ho.GetMaxElementOrder();
    const int nd1d = order + 1;
@@ -94,7 +95,7 @@ void BatchedLORAssembly::FormLORVertexCoordinates(FiniteElementSpace &fes_ho,
    IntegrationRule ir = GetCollocatedIntRule(fes_ho);
 
    // Map from nodal E-vector to Q-vector at the LOR vertex points
-   X_vert.SetSize(dim*ndof_per_el*nel_ho);
+   X_vert.SetSize(sdim*ndof_per_el*nel_ho);
    const QuadratureInterpolator *quad_interp =
       nodal_fes->GetQuadratureInterpolator(ir);
    quad_interp->SetOutputLayout(QVectorLayout::byVDIM);
@@ -380,44 +381,49 @@ void BatchedLORAssembly::SparseIJToCSR(OperatorHandle &A) const
    FillJAndData(*A_mat);
 }
 
+template <int ORDER, int SDIM, typename LOR_KERNEL>
+static void Assemble_(LOR_KERNEL &kernel, int dim)
+{
+   if (dim == 2) { kernel.template Assemble2D<ORDER,SDIM>(); }
+   else if (dim == 3) { kernel.template Assemble3D<ORDER>(); }
+   else { MFEM_ABORT("Unsupported dimension"); }
+}
+
+template <int ORDER, typename LOR_KERNEL>
+static void Assemble_(LOR_KERNEL &kernel, int dim, int sdim)
+{
+   if (sdim == 2) { Assemble_<ORDER,2>(kernel, dim); }
+   else if (sdim == 3) { Assemble_<ORDER,3>(kernel, dim); }
+   else { MFEM_ABORT("Unsupported space dimension."); }
+}
+
+template <typename LOR_KERNEL>
+static void Assemble_(LOR_KERNEL &kernel, int dim, int sdim, int order)
+{
+   switch (order)
+   {
+      case 1: Assemble_<1>(kernel, dim, sdim); break;
+      case 2: Assemble_<2>(kernel, dim, sdim); break;
+      case 3: Assemble_<3>(kernel, dim, sdim); break;
+      case 4: Assemble_<4>(kernel, dim, sdim); break;
+      case 5: Assemble_<5>(kernel, dim, sdim); break;
+      case 6: Assemble_<6>(kernel, dim, sdim); break;
+      case 7: Assemble_<7>(kernel, dim, sdim); break;
+      case 8: Assemble_<8>(kernel, dim, sdim); break;
+      default: MFEM_ABORT("No kernel order " << order << "!");
+   }
+}
+
 template <typename LOR_KERNEL>
 void BatchedLORAssembly::AssemblyKernel(BilinearForm &a)
 {
    LOR_KERNEL kernel(a, fes_ho, X_vert, sparse_ij, sparse_mapping);
 
    const int dim = fes_ho.GetMesh()->Dimension();
+   const int sdim = fes_ho.GetMesh()->SpaceDimension();
    const int order = fes_ho.GetMaxElementOrder();
 
-   if (dim == 2)
-   {
-      switch (order)
-      {
-         case 1: kernel.template Assemble2D<1>(); break;
-         case 2: kernel.template Assemble2D<2>(); break;
-         case 3: kernel.template Assemble2D<3>(); break;
-         case 4: kernel.template Assemble2D<4>(); break;
-         case 5: kernel.template Assemble2D<5>(); break;
-         case 6: kernel.template Assemble2D<6>(); break;
-         case 7: kernel.template Assemble2D<7>(); break;
-         case 8: kernel.template Assemble2D<8>(); break;
-         default: MFEM_ABORT("No kernel order " << order << "!");
-      }
-   }
-   else if (dim == 3)
-   {
-      switch (order)
-      {
-         case 1: kernel.template Assemble3D<1>(); break;
-         case 2: kernel.template Assemble3D<2>(); break;
-         case 3: kernel.template Assemble3D<3>(); break;
-         case 4: kernel.template Assemble3D<4>(); break;
-         case 5: kernel.template Assemble3D<5>(); break;
-         case 6: kernel.template Assemble3D<6>(); break;
-         case 7: kernel.template Assemble3D<7>(); break;
-         case 8: kernel.template Assemble3D<8>(); break;
-         default: MFEM_ABORT("No kernel order " << order << "!");
-      }
-   }
+   Assemble_(kernel, dim, sdim, order);
 }
 
 void BatchedLORAssembly::AssembleWithoutBC(BilinearForm &a, OperatorHandle &A)

fem/lor/lor_h1.hpp

@@ -22,15 +22,22 @@ namespace mfem
 class BatchedLOR_H1 : BatchedLORKernel
 {
 public:
-   template <int ORDER> void Assemble2D();
+   template <int ORDER, int SDIM> void Assemble2D();
    template <int ORDER> void Assemble3D();
    BatchedLOR_H1(BilinearForm &a,
                  FiniteElementSpace &fes_ho_,
                  Vector &X_vert_,
                  Vector &sparse_ij_,
-                 Array<int> &sparse_mapping_);
+                 Array<int> &sparse_mapping_)
+      : BatchedLORKernel(fes_ho_, X_vert_, sparse_ij_, sparse_mapping_)
+   {
+      ProjectLORCoefficient<MassIntegrator>(a, c1);
+      ProjectLORCoefficient<DiffusionIntegrator>(a, c2);
+   }
 };
 
 }
 
+#include "lor_h1_impl.hpp"
+
 #endif

fem/lor/lor_h1_impl.hpp

@@ -9,7 +9,6 @@
 // terms of the BSD-3 license. We welcome feedback and contributions, see file
 // CONTRIBUTING.md for details.
 
-#include "lor_h1.hpp"
 #include "lor_util.hpp"
 #include "../../linalg/dtensor.hpp"
 #include "../../general/forall.hpp"
@@ -17,7 +16,7 @@
 namespace mfem
 {
 
-template <int ORDER>
+template <int ORDER, int SDIM>
 void BatchedLOR_H1::Assemble2D()
 {
    const int nel_ho = fes_ho.GetNE();
@@ -74,31 +73,8 @@ void BatchedLOR_H1::Assemble2D()
 
             for (int i=0; i<sz_local_mat; ++i) { local_mat[i] = 0.0; }
 
-            double vx[4], vy[4];
-            LORVertexCoordinates2D<ORDER>(X, iel_ho, kx, ky, vx, vy);
+            SetupLORQuadData2D<ORDER,SDIM,false,false>(X, iel_ho, kx, ky, Q, false);
 
-            for (int iqy=0; iqy<2; ++iqy)
-            {
-               for (int iqx=0; iqx<2; ++iqx)
-               {
-                  const double x = iqx;
-                  const double y = iqy;
-                  const double w = 1.0/4.0;
-
-                  double J_[2*2];
-                  DeviceTensor<2> J(J_, 2, 2);
-
-                  Jacobian2D(x, y, vx, vy, J);
-
-                  const double detJ = Det2D(J);
-                  const double w_detJ = w/detJ;
-
-                  Q(0,iqy,iqx) = w_detJ * (J(0,1)*J(0,1) + J(1,1)*J(1,1)); // 1,1
-                  Q(1,iqy,iqx) = -w_detJ * (J(0,1)*J(0,0) + J(1,1)*J(1,0)); // 1,2
-                  Q(2,iqy,iqx) = w_detJ * (J(0,0)*J(0,0) + J(1,0)*J(1,0)); // 2,2
-                  Q(3,iqy,iqx) = w*detJ;
-               }
-            }
             for (int iqx=0; iqx<2; ++iqx)
             {
                for (int iqy=0; iqy<2; ++iqy)
@@ -519,34 +495,4 @@ void BatchedLOR_H1::Assemble3D()
    }
 }
 
-// Explicit template instantiations
-template void BatchedLOR_H1::Assemble2D<1>();
-template void BatchedLOR_H1::Assemble2D<2>();
-template void BatchedLOR_H1::Assemble2D<3>();
-template void BatchedLOR_H1::Assemble2D<4>();
-template void BatchedLOR_H1::Assemble2D<5>();
-template void BatchedLOR_H1::Assemble2D<6>();
-template void BatchedLOR_H1::Assemble2D<7>();
-template void BatchedLOR_H1::Assemble2D<8>();
-
-template void BatchedLOR_H1::Assemble3D<1>();
-template void BatchedLOR_H1::Assemble3D<2>();
-template void BatchedLOR_H1::Assemble3D<3>();
-template void BatchedLOR_H1::Assemble3D<4>();
-template void BatchedLOR_H1::Assemble3D<5>();
-template void BatchedLOR_H1::Assemble3D<6>();
-template void BatchedLOR_H1::Assemble3D<7>();
-template void BatchedLOR_H1::Assemble3D<8>();
-
-BatchedLOR_H1::BatchedLOR_H1(BilinearForm &a,
-                             FiniteElementSpace &fes_ho_,
-                             Vector &X_vert_,
-                             Vector &sparse_ij_,
-                             Array<int> &sparse_mapping_)
-   : BatchedLORKernel(fes_ho_, X_vert_, sparse_ij_, sparse_mapping_)
-{
-   ProjectLORCoefficient<MassIntegrator>(a, c1);
-   ProjectLORCoefficient<DiffusionIntegrator>(a, c2);
-}
-
 } // namespace mfem

fem/lor/lor_nd.hpp

@@ -22,15 +22,22 @@ namespace mfem
 class BatchedLOR_ND : BatchedLORKernel
 {
 public:
-   template <int ORDER> void Assemble2D();
+   template <int ORDER, int SDIM> void Assemble2D();
    template <int ORDER> void Assemble3D();
    BatchedLOR_ND(BilinearForm &a,
                  FiniteElementSpace &fes_ho_,
                  Vector &X_vert_,
                  Vector &sparse_ij_,
-                 Array<int> &sparse_mapping_);
+                 Array<int> &sparse_mapping_)
+      : BatchedLORKernel(fes_ho_, X_vert_, sparse_ij_, sparse_mapping_)
+   {
+      ProjectLORCoefficient<VectorFEMassIntegrator>(a, c1);
+      ProjectLORCoefficient<CurlCurlIntegrator>(a, c2);
+   }
 };
 
 }
 
+#include "lor_nd_impl.hpp"
+
 #endif

fem/lor/lor_nd_impl.hpp

@@ -9,7 +9,6 @@
 // terms of the BSD-3 license. We welcome feedback and contributions, see file
 // CONTRIBUTING.md for details.
 
-#include "lor_nd.hpp"
 #include "lor_util.hpp"
 #include "../../linalg/dtensor.hpp"
 #include "../../general/forall.hpp"
@@ -17,7 +16,7 @@
 namespace mfem
 {
 
-template <int ORDER>
+template <int ORDER, int SDIM>
 void BatchedLOR_ND::Assemble2D()
 {
    const int nel_ho = fes_ho.GetNE();
@@ -83,31 +82,8 @@ void BatchedLOR_ND::Assemble2D()
             // local_mat is the local (dense) stiffness matrix
             for (int i=0; i<sz_local_mat; ++i) { local_mat[i] = 0.0; }
 
-            double vx[4], vy[4];
-            LORVertexCoordinates2D<ORDER>(X, iel_ho, kx, ky, vx, vy);
+            SetupLORQuadData2D<ORDER,SDIM,false,true>(X, iel_ho, kx, ky, Q, true);
 
-            for (int iqx=0; iqx<2; ++iqx)
-            {
-               for (int iqy=0; iqy<2; ++iqy)
-               {
-                  const double x = iqx;
-                  const double y = iqy;
-                  const double w = 1.0/4.0;
-
-                  double J_[2*2];
-                  DeviceTensor<2> J(J_, 2, 2);
-
-                  Jacobian2D(x, y, vx, vy, J);
-
-                  const double detJ = Det2D(J);
-                  const double w_detJ = w/detJ;
-
-                  Q(0,iqy,iqx) = w_detJ * (J(0,1)*J(0,1) + J(1,1)*J(1,1)); // 1,1
-                  Q(1,iqy,iqx) = -w_detJ * (J(0,1)*J(0,0) + J(1,1)*J(1,0)); // 1,2
-                  Q(2,iqy,iqx) = w_detJ * (J(0,0)*J(0,0) + J(1,0)*J(1,0)); // 2,2
-                  Q(3,iqy,iqx) = w_detJ;
-               }
-            }
             for (int iqx=0; iqx<2; ++iqx)
             {
                for (int iqy=0; iqy<2; ++iqy)
@@ -563,34 +539,4 @@ void BatchedLOR_ND::Assemble3D()
    }
 }
 
-// Explicit template instantiations
-template void BatchedLOR_ND::Assemble2D<1>();
-template void BatchedLOR_ND::Assemble2D<2>();
-template void BatchedLOR_ND::Assemble2D<3>();
-template void BatchedLOR_ND::Assemble2D<4>();
-template void BatchedLOR_ND::Assemble2D<5>();
-template void BatchedLOR_ND::Assemble2D<6>();
-template void BatchedLOR_ND::Assemble2D<7>();
-template void BatchedLOR_ND::Assemble2D<8>();
-
-template void BatchedLOR_ND::Assemble3D<1>();
-template void BatchedLOR_ND::Assemble3D<2>();
-template void BatchedLOR_ND::Assemble3D<3>();
-template void BatchedLOR_ND::Assemble3D<4>();
-template void BatchedLOR_ND::Assemble3D<5>();
-template void BatchedLOR_ND::Assemble3D<6>();
-template void BatchedLOR_ND::Assemble3D<7>();
-template void BatchedLOR_ND::Assemble3D<8>();
-
-BatchedLOR_ND::BatchedLOR_ND(BilinearForm &a,
-                             FiniteElementSpace &fes_ho_,
-                             Vector &X_vert_,
-                             Vector &sparse_ij_,
-                             Array<int> &sparse_mapping_)
-   : BatchedLORKernel(fes_ho_, X_vert_, sparse_ij_, sparse_mapping_)
-{
-   ProjectLORCoefficient<VectorFEMassIntegrator>(a, c1);
-   ProjectLORCoefficient<CurlCurlIntegrator>(a, c2);
-}
-
 } // namespace mfem

fem/lor/lor_rt.hpp

@@ -22,15 +22,22 @@ namespace mfem
 class BatchedLOR_RT : BatchedLORKernel
 {
 public:
-   template <int ORDER> void Assemble2D();
+   template <int ORDER, int SDIM> void Assemble2D();
    template <int ORDER> void Assemble3D();
    BatchedLOR_RT(BilinearForm &a,
                  FiniteElementSpace &fes_ho_,
                  Vector &X_vert_,
                  Vector &sparse_ij_,
-                 Array<int> &sparse_mapping_);
+                 Array<int> &sparse_mapping_)
+      : BatchedLORKernel(fes_ho_, X_vert_, sparse_ij_, sparse_mapping_)
+   {
+      ProjectLORCoefficient<VectorFEMassIntegrator>(a, c1);
+      ProjectLORCoefficient<DivDivIntegrator>(a, c2);
+   }
 };
 
 }
 
+#include "lor_rt_impl.hpp"
+
 #endif

fem/lor/lor_rt_impl.hpp

@@ -9,7 +9,6 @@
 // terms of the BSD-3 license. We welcome feedback and contributions, see file
 // CONTRIBUTING.md for details.
 
-#include "lor_rt.hpp"
 #include "lor_util.hpp"
 #include "../../linalg/dtensor.hpp"
 #include "../../general/forall.hpp"
@@ -17,7 +16,7 @@
 namespace mfem
 {
 
-template <int ORDER>
+template <int ORDER, int SDIM>
 void BatchedLOR_RT::Assemble2D()
 {
    const int nel_ho = fes_ho.GetNE();
@@ -79,31 +78,8 @@ void BatchedLOR_RT::Assemble2D()
             // local_mat is the local (dense) stiffness matrix
             for (int i=0; i<sz_local_mat; ++i) { local_mat[i] = 0.0; }
 
-            double vx[4], vy[4];
-            LORVertexCoordinates2D<ORDER>(X, iel_ho, kx, ky, vx, vy);
+            SetupLORQuadData2D<ORDER,SDIM,true,false>(X, iel_ho, kx, ky, Q, true);
 
-            for (int iqx=0; iqx<2; ++iqx)
-            {
-               for (int iqy=0; iqy<2; ++iqy)
-               {
-                  const double x = iqx;
-                  const double y = iqy;
-                  const double w = 1.0/4.0;
-
-                  double J_[2*2];
-                  DeviceTensor<2> J(J_, 2, 2);
-
-                  Jacobian2D(x, y, vx, vy, J);
-
-                  const double detJ = Det2D(J);
-                  const double w_detJ = w/detJ;
-
-                  Q(0,iqy,iqx) = w_detJ * (J(0,0)*J(0,0) + J(1,0)*J(1,0)); // 1,1
-                  Q(1,iqy,iqx) = w_detJ * (J(0,0)*J(0,1) + J(1,0)*J(1,1)); // 1,2
-                  Q(2,iqy,iqx) = w_detJ * (J(0,1)*J(0,1) + J(1,1)*J(1,1)); // 2,2
-                  Q(3,iqy,iqx) = w_detJ;
-               }
-            }
             for (int iqx=0; iqx<2; ++iqx)
             {
                for (int iqy=0; iqy<2; ++iqy)
@@ -547,34 +523,4 @@ void BatchedLOR_RT::Assemble3D()
    }
 }
 
-// Explicit template instantiations
-template void BatchedLOR_RT::Assemble2D<1>();
-template void BatchedLOR_RT::Assemble2D<2>();
-template void BatchedLOR_RT::Assemble2D<3>();
-template void BatchedLOR_RT::Assemble2D<4>();
-template void BatchedLOR_RT::Assemble2D<5>();
-template void BatchedLOR_RT::Assemble2D<6>();
-template void BatchedLOR_RT::Assemble2D<7>();
-template void BatchedLOR_RT::Assemble2D<8>();
-
-template void BatchedLOR_RT::Assemble3D<1>();
-template void BatchedLOR_RT::Assemble3D<2>();
-template void BatchedLOR_RT::Assemble3D<3>();
-template void BatchedLOR_RT::Assemble3D<4>();
-template void BatchedLOR_RT::Assemble3D<5>();
-template void BatchedLOR_RT::Assemble3D<6>();
-template void BatchedLOR_RT::Assemble3D<7>();
-template void BatchedLOR_RT::Assemble3D<8>();
-
-BatchedLOR_RT::BatchedLOR_RT(BilinearForm &a,
-                             FiniteElementSpace &fes_ho_,
-                             Vector &X_vert_,
-                             Vector &sparse_ij_,
-                             Array<int> &sparse_mapping_)
-   : BatchedLORKernel(fes_ho_, X_vert_, sparse_ij_, sparse_mapping_)
-{
-   ProjectLORCoefficient<VectorFEMassIntegrator>(a, c1);
-   ProjectLORCoefficient<DivDivIntegrator>(a, c2);
-}
-
 } // namespace mfem

fem/lor/lor_util.hpp

@@ -20,11 +20,22 @@
 namespace mfem
 {
 
-template <int ORDER>
-MFEM_HOST_DEVICE inline void LORVertexCoordinates2D(
-   const double *X, int iel_ho, int kx, int ky, double vx[4], double vy[4])
+MFEM_HOST_DEVICE inline double Det2D(DeviceMatrix &J)
+{
+   return J(0,0)*J(1,1) - J(1,0)*J(0,1);
+}
+
+MFEM_HOST_DEVICE inline double Det3D(DeviceMatrix &J)
+{
+   return J(0,0) * (J(1,1) * J(2,2) - J(2,1) * J(1,2)) -
+          J(1,0) * (J(0,1) * J(2,2) - J(2,1) * J(0,2)) +
+          J(2,0) * (J(0,1) * J(1,2) - J(1,1) * J(0,2));
+}
+
+template <int ORDER, int SDIM=2>
+MFEM_HOST_DEVICE inline void LORVertexCoordinates2D(
+   const double *X, int iel_ho, int kx, int ky, double **v)
 {
-   const int dim = 2;
    const int nd1d = ORDER + 1;
    const int nvert_per_el = nd1d*nd1d;
 
@@ -33,23 +44,31 @@ MFEM_HOST_DEVICE inline void LORVertexCoordinates2D(
    const int v2 = kx + 1 + nd1d*(ky + 1);
    const int v3 = kx + nd1d*(ky + 1);
 
-   const int e0 = dim*(v0 + nvert_per_el*iel_ho);
-   const int e1 = dim*(v1 + nvert_per_el*iel_ho);
-   const int e2 = dim*(v2 + nvert_per_el*iel_ho);
-   const int e3 = dim*(v3 + nvert_per_el*iel_ho);
+   const int e0 = SDIM*(v0 + nvert_per_el*iel_ho);
+   const int e1 = SDIM*(v1 + nvert_per_el*iel_ho);
+   const int e2 = SDIM*(v2 + nvert_per_el*iel_ho);
+   const int e3 = SDIM*(v3 + nvert_per_el*iel_ho);
 
    // Vertex coordinates
-   vx[0] = X[e0 + 0];
-   vy[0] = X[e0 + 1];
+   v[0][0] = X[e0 + 0];
+   v[1][0] = X[e0 + 1];
 
-   vx[1] = X[e1 + 0];
-   vy[1] = X[e1 + 1];
+   v[0][1] = X[e1 + 0];
+   v[1][1] = X[e1 + 1];
 
-   vx[2] = X[e2 + 0];
-   vy[2] = X[e2 + 1];
+   v[0][2] = X[e2 + 0];
+   v[1][2] = X[e2 + 1];
 
-   vx[3] = X[e3 + 0];
-   vy[3] = X[e3 + 1];
+   v[0][3] = X[e3 + 0];
+   v[1][3] = X[e3 + 1];
+
+   if (SDIM == 3)
+   {
+      v[2][0] = X[e0 + 2];
+      v[2][1] = X[e1 + 2];
+      v[2][2] = X[e2 + 2];
+      v[2][3] = X[e3 + 2];
+   }
 }
 
 template <int ORDER>
@@ -112,15 +131,80 @@ MFEM_HOST_DEVICE inline void LORVertexCoordinates3D(
    vz[7] = X[e7 + 2];
 }
 
+template <int SDIM=2>
 MFEM_HOST_DEVICE inline void Jacobian2D(
-   const double x, const double y, const double vx[4], const double vy[4],
-   DeviceMatrix &J)
-{
-   J(0,0) = -(1-y)*vx[0] + (1-y)*vx[1] + y*vx[2] - y*vx[3];
-   J(0,1) = -(1-x)*vx[0] - x*vx[1] + x*vx[2] + (1-x)*vx[3];
+   const double x, const double y, double **v, DeviceMatrix &J);
 
-   J(1,0) = -(1-y)*vy[0] + (1-y)*vy[1] + y*vy[2] - y*vy[3];
-   J(1,1) = -(1-x)*vy[0] - x*vy[1] + x*vy[2] + (1-x)*vy[3];
+template <> MFEM_HOST_DEVICE inline void Jacobian2D<2>(
+   const double x, const double y, double **v, DeviceMatrix &J)
+{
+   J(0,0) = -(1-y)*v[0][0] + (1-y)*v[0][1] + y*v[0][2] - y*v[0][3];
+   J(0,1) = -(1-x)*v[0][0] - x*v[0][1] + x*v[0][2] + (1-x)*v[0][3];
+
+   J(1,0) = -(1-y)*v[1][0] + (1-y)*v[1][1] + y*v[1][2] - y*v[1][3];
+   J(1,1) = -(1-x)*v[1][0] - x*v[1][1] + x*v[1][2] + (1-x)*v[1][3];
+}
+
+template <> MFEM_HOST_DEVICE inline void Jacobian2D<3>(
+   const double x, const double y, double **v, DeviceMatrix &J)
+{
+   J(0,0) = -(1-y)*v[0][0] + (1-y)*v[0][1] + y*v[0][2] - y*v[0][3];
+   J(0,1) = -(1-x)*v[0][0] - x*v[0][1] + x*v[0][2] + (1-x)*v[0][3];
+
+   J(1,0) = -(1-y)*v[1][0] + (1-y)*v[1][1] + y*v[1][2] - y*v[1][3];
+   J(1,1) = -(1-x)*v[1][0] - x*v[1][1] + x*v[1][2] + (1-x)*v[1][3];
+
+   J(2,0) = -(1-y)*v[2][0] + (1-y)*v[2][1] + y*v[2][2] - y*v[2][3];
+   J(2,1) = -(1-x)*v[2][0] - x*v[2][1] + x*v[2][2] + (1-x)*v[2][3];
+}
+
+template <int ORDER, int SDIM, bool RT, bool ND>
+MFEM_HOST_DEVICE inline void SetupLORQuadData2D(
+   const double *X, int iel_ho, int kx, int ky, DeviceTensor<3> &Q, bool piola)
+{
+   double vx[4], vy[4], vz[4];
+   double *v[] = {vx, vy, vz};
+   LORVertexCoordinates2D<ORDER,SDIM>(X, iel_ho, kx, ky, v);
+
+   for (int iqy=0; iqy<2; ++iqy)
+   {
+      for (int iqx=0; iqx<2; ++iqx)
+      {
+         const double x = iqx;
+         const double y = iqy;
+         const double w = 1.0/4.0;
+
+         double J_[SDIM*2];
+         DeviceTensor<2> J(J_, SDIM, 2);
+
+         Jacobian2D<SDIM>(x, y, v, J);
+
+         if (SDIM == 2)
+         {
+            const double detJ = Det2D(J);
+            const double w_detJ = w/detJ;
+            const double E = J(0,0)*J(0,0) + J(1,0)*J(1,0);
+            const double F = J(0,0)*J(0,1) + J(1,0)*J(1,1);
+            const double G = J(0,1)*J(0,1) + J(1,1)*J(1,1);
+            Q(0,iqy,iqx) = w_detJ * (RT ? E : G); // 1,1
+            Q(1,iqy,iqx) = w_detJ * (RT ? F : -F); // 1,2
+            Q(2,iqy,iqx) = w_detJ * (RT ? G : E); // 2,2
+            Q(3,iqy,iqx) = (ND || RT) ? w_detJ : w*detJ;
+         }
+         else
+         {
+            const double E = J(0,0)*J(0,0) + J(1,0)*J(1,0) + J(2,0)*J(2,0);
+            const double F = J(0,0)*J(0,1) + J(1,0)*J(1,1) + J(2,0)*J(2,1);
+            const double G = J(0,1)*J(0,1) + J(1,1)*J(1,1) + J(2,1)*J(2,1);
+            const double detJ = sqrt(E*G - F*F);
+            const double w_detJ = w/detJ;
+            Q(0,iqy,iqx) = w_detJ * (RT ? E : G); // 1,1
+            Q(1,iqy,iqx) = w_detJ * (RT ? F : -F); // 1,2
+            Q(2,iqy,iqx) = w_detJ * (RT ? G : E); // 2,2
+            Q(3,iqy,iqx) =  (ND || RT) ? w_detJ : w*detJ;
+         }
+      }
+   }
 }
 
 MFEM_HOST_DEVICE inline void Jacobian3D(
@@ -180,18 +264,6 @@ MFEM_HOST_DEVICE inline void Adjugate3D(const DeviceMatrix &J, DeviceMatrix &A)
    A(2,2) = (J(0,0) * J(1,1)) - (J(0,1) * J(1,0));
 }
 
-MFEM_HOST_DEVICE inline double Det2D(DeviceMatrix &J)
-{
-   return J(0,0)*J(1,1) - J(1,0)*J(0,1);
-}
-
-MFEM_HOST_DEVICE inline double Det3D(DeviceMatrix &J)
-{
-   return J(0,0) * (J(1,1) * J(2,2) - J(2,1) * J(1,2)) -
-          J(1,0) * (J(0,1) * J(2,2) - J(2,1) * J(0,2)) +
-          J(2,0) * (J(0,1) * J(1,2) - J(1,1) * J(0,2));
-}
-
 }
 
 #endif

fem/nonlinearform.cpp

@@ -97,12 +97,37 @@ double NonlinearForm::GetGridFunctionEnergy(const Vector &x) const
    const FiniteElement *fe;
    ElementTransformation *T;
    DofTransformation *doftrans;
+   Mesh *mesh = fes->GetMesh();
    double energy = 0.0;
 
    if (dnfi.Size())
    {
+      // Which attributes need to be processed?
+      Array<int> attr_marker(mesh->attributes.Size() ?
+                             mesh->attributes.Max() : 0);
+      attr_marker = 0;
+      for (int k = 0; k < dnfi.Size(); k++)
+      {
+         if (dnfi_marker[k] == NULL)
+         {
+            attr_marker = 1;
+            break;
+         }
+         Array<int> &marker = *dnfi_marker[k];
+         MFEM_ASSERT(marker.Size() == attr_marker.Size(),
+                     "invalid marker for domain integrator #"
+                     << k << ", counting from zero");
+         for (int i = 0; i < attr_marker.Size(); i++)
+         {
+            attr_marker[i] |= marker[i];
+         }
+      }
+
       for (int i = 0; i < fes->GetNE(); i++)
       {
+         const int attr = mesh->GetAttribute(i);
+         if (attr_marker[attr-1] == 0) { continue; }
+
          fe = fes->GetFE(i);
          doftrans = fes->GetElementVDofs(i, vdofs);
          T = fes->GetElementTransformation(i);
@@ -110,6 +135,9 @@ double NonlinearForm::GetGridFunctionEnergy(const Vector &x) const
          if (doftrans) {doftrans->InvTransformPrimal(el_x); }
          for (int k = 0; k < dnfi.Size(); k++)
          {
+            if (dnfi_marker[k] &&
+                (*dnfi_marker[k])[attr-1] == 0) { continue; }
+
             energy += dnfi[k]->GetElementEnergy(*fe, *T, el_x);
          }
       }
@@ -175,8 +203,32 @@ void NonlinearForm::Mult(const Vector &x, Vector &y) const
 
    if (dnfi.Size())
    {
+      // Which attributes need to be processed?
+      Array<int> attr_marker(mesh->attributes.Size() ?
+                             mesh->attributes.Max() : 0);
+      attr_marker = 0;
+      for (int k = 0; k < dnfi.Size(); k++)
+      {
+         if (dnfi_marker[k] == NULL)
+         {
+            attr_marker = 1;
+            break;
+         }
+         Array<int> &marker = *dnfi_marker[k];
+         MFEM_ASSERT(marker.Size() == attr_marker.Size(),
+                     "invalid marker for domain integrator #"
+                     << k << ", counting from zero");
+         for (int i = 0; i < attr_marker.Size(); i++)
+         {
+            attr_marker[i] |= marker[i];
+         }
+      }
+
       for (int i = 0; i < fes->GetNE(); i++)
       {
+         const int attr = mesh->GetAttribute(i);
+         if (attr_marker[attr-1] == 0) { continue; }
+
          fe = fes->GetFE(i);
          doftrans = fes->GetElementVDofs(i, vdofs);
          T = fes->GetElementTransformation(i);
@@ -184,6 +236,9 @@ void NonlinearForm::Mult(const Vector &x, Vector &y) const
          if (doftrans) {doftrans->InvTransformPrimal(el_x); }
          for (int k = 0; k < dnfi.Size(); k++)
          {
+            if (dnfi_marker[k] &&
+                (*dnfi_marker[k])[attr-1] == 0) { continue; }
+
             dnfi[k]->AssembleElementVector(*fe, *T, el_x, el_y);
             if (doftrans) {doftrans->TransformDual(el_y); }
             py.AddElementVector(vdofs, el_y);
@@ -322,8 +377,32 @@ Operator &NonlinearForm::GetGradient(const Vector &x) const
 
    if (dnfi.Size())
    {
+      // Which attributes need to be processed?
+      Array<int> attr_marker(mesh->attributes.Size() ?
+                             mesh->attributes.Max() : 0);
+      attr_marker = 0;
+      for (int k = 0; k < dnfi.Size(); k++)
+      {
+         if (dnfi_marker[k] == NULL)
+         {
+            attr_marker = 1;
+            break;
+         }
+         Array<int> &marker = *dnfi_marker[k];
+         MFEM_ASSERT(marker.Size() == attr_marker.Size(),
+                     "invalid marker for domain integrator #"
+                     << k << ", counting from zero");
+         for (int i = 0; i < attr_marker.Size(); i++)
+         {
+            attr_marker[i] |= marker[i];
+         }
+      }
+
       for (int i = 0; i < fes->GetNE(); i++)
       {
+         const int attr = mesh->GetAttribute(i);
+         if (attr_marker[attr-1] == 0) { continue; }
+
          fe = fes->GetFE(i);
          doftrans = fes->GetElementVDofs(i, vdofs);
          T = fes->GetElementTransformation(i);
@@ -331,6 +410,9 @@ Operator &NonlinearForm::GetGradient(const Vector &x) const
          if (doftrans) {doftrans->InvTransformPrimal(el_x); }
          for (int k = 0; k < dnfi.Size(); k++)
          {
+            if (dnfi_marker[k] &&
+                (*dnfi_marker[k])[attr-1] == 0) { continue; }
+
             dnfi[k]->AssembleElementGrad(*fe, *T, el_x, elmat);
             if (doftrans) { doftrans->TransformDual(elmat); }
             Grad->AddSubMatrix(vdofs, vdofs, elmat, skip_zeros);
@@ -561,13 +643,6 @@ BlockNonlinearForm::BlockNonlinearForm(Array<FiniteElementSpace *> &f) :
    SetSpaces(f);
 }
 
-void BlockNonlinearForm::AddBdrFaceIntegrator(BlockNonlinearFormIntegrator *nfi,
-                                              Array<int> &bdr_attr_marker)
-{
-   bfnfi.Append(nfi);
-   bfnfi_marker.Append(&bdr_attr_marker);
-}
-
 void BlockNonlinearForm::SetEssentialBC(
    const Array<Array<int> *> &bdr_attr_is_ess, Array<Vector *> &rhs)
 {
@@ -592,6 +667,7 @@ double BlockNonlinearForm::GetEnergyBlocked(const BlockVector &bx) const
    Array<const FiniteElement *> fe(fes.Size());
    ElementTransformation *T;
    DofTransformation *doftrans;
+   Mesh *mesh = fes[0]->GetMesh();
    double energy = 0.0;
 
    for (int i=0; i<fes.Size(); ++i)
@@ -601,8 +677,33 @@ double BlockNonlinearForm::GetEnergyBlocked(const BlockVector &bx) const
    }
 
    if (dnfi.Size())
+   {
+      // Which attributes need to be processed?
+      Array<int> attr_marker(mesh->attributes.Size() ?
+                             mesh->attributes.Max() : 0);
+      attr_marker = 0;
+      for (int k = 0; k < dnfi.Size(); k++)
+      {
+         if (dnfi_marker[k] == NULL)
+         {
+            attr_marker = 1;
+            break;
+         }
+         Array<int> &marker = *dnfi_marker[k];
+         MFEM_ASSERT(marker.Size() == attr_marker.Size(),
+                     "invalid marker for domain integrator #"
+                     << k << ", counting from zero");
+         for (int i = 0; i < attr_marker.Size(); i++)
+         {
+            attr_marker[i] |= marker[i];
+         }
+      }
+
       for (int i = 0; i < fes[0]->GetNE(); ++i)
       {
+         const int attr = mesh->GetAttribute(i);
+         if (attr_marker[attr-1] == 0) { continue; }
+
          T = fes[0]->GetElementTransformation(i);
          for (int s=0; s<fes.Size(); ++s)
          {
@@ -614,9 +715,13 @@ double BlockNonlinearForm::GetEnergyBlocked(const BlockVector &bx) const
 
          for (int k = 0; k < dnfi.Size(); ++k)
          {
+            if (dnfi_marker[k] &&
+                (*dnfi_marker[k])[attr-1] == 0) { continue; }
+
             energy += dnfi[k]->GetElementEnergy(fe, *T, el_x_const);
          }
       }
+   }
 
    // free the allocated memory
    for (int i = 0; i < fes.Size(); ++i)
@@ -656,6 +761,7 @@ void BlockNonlinearForm::MultBlocked(const BlockVector &bx,
    Array<const FiniteElement *> fe2(fes.Size());
    ElementTransformation *T;
    Array<DofTransformation *> doftrans(fes.Size()); doftrans = nullptr;
+   Mesh *mesh = fes[0]->GetMesh();
 
    by.UseDevice(true);
    by = 0.0;
@@ -670,8 +776,32 @@ void BlockNonlinearForm::MultBlocked(const BlockVector &bx,
 
    if (dnfi.Size())
    {
+      // Which attributes need to be processed?
+      Array<int> attr_marker(mesh->attributes.Size() ?
+                             mesh->attributes.Max() : 0);
+      attr_marker = 0;
+      for (int k = 0; k < dnfi.Size(); k++)
+      {
+         if (dnfi_marker[k] == NULL)
+         {
+            attr_marker = 1;
+            break;
+         }
+         Array<int> &marker = *dnfi_marker[k];
+         MFEM_ASSERT(marker.Size() == attr_marker.Size(),
+                     "invalid marker for domain integrator #"
+                     << k << ", counting from zero");
+         for (int i = 0; i < attr_marker.Size(); i++)
+         {
+            attr_marker[i] |= marker[i];
+         }
+      }
+
       for (int i = 0; i < fes[0]->GetNE(); ++i)
       {
+         const int attr = mesh->GetAttribute(i);
+         if (attr_marker[attr-1] == 0) { continue; }
+
          T = fes[0]->GetElementTransformation(i);
          for (int s = 0; s < fes.Size(); ++s)
          {
@@ -683,6 +813,9 @@ void BlockNonlinearForm::MultBlocked(const BlockVector &bx,
 
          for (int k = 0; k < dnfi.Size(); ++k)
          {
+            if (dnfi_marker[k] &&
+                (*dnfi_marker[k])[attr-1] == 0) { continue; }
+
             dnfi[k]->AssembleElementVector(fe, *T,
                                            el_x_const, el_y);
 
@@ -698,7 +831,6 @@ void BlockNonlinearForm::MultBlocked(const BlockVector &bx,
 
    if (fnfi.Size())
    {
-      Mesh *mesh = fes[0]->GetMesh();
       FaceElementTransformations *tr;
 
       for (int i = 0; i < mesh->GetNumFaces(); ++i)
@@ -736,8 +868,8 @@ void BlockNonlinearForm::MultBlocked(const BlockVector &bx,
 
    if (bfnfi.Size())
    {
-      Mesh *mesh = fes[0]->GetMesh();
       FaceElementTransformations *tr;
+
       // Which boundary attributes need to be processed?
       Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
                                  mesh->bdr_attributes.Max() : 0);
@@ -858,6 +990,7 @@ void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx) const
    Array<const FiniteElement *>fe2(fes.Size());
    ElementTransformation * T;
    Array<DofTransformation *> doftrans(fes.Size()); doftrans = nullptr;
+   Mesh *mesh = fes[0]->GetMesh();
 
    for (int i=0; i<fes.Size(); ++i)
    {
@@ -888,8 +1021,32 @@ void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx) const
 
    if (dnfi.Size())
    {
+      // Which attributes need to be processed?
+      Array<int> attr_marker(mesh->attributes.Size() ?
+                             mesh->attributes.Max() : 0);
+      attr_marker = 0;
+      for (int k = 0; k < dnfi.Size(); k++)
+      {
+         if (dnfi_marker[k] == NULL)
+         {
+            attr_marker = 1;
+            break;
+         }
+         Array<int> &marker = *dnfi_marker[k];
+         MFEM_ASSERT(marker.Size() == attr_marker.Size(),
+                     "invalid marker for domain integrator #"
+                     << k << ", counting from zero");
+         for (int i = 0; i < attr_marker.Size(); i++)
+         {
+            attr_marker[i] |= marker[i];
+         }
+      }
+
       for (int i = 0; i < fes[0]->GetNE(); ++i)
       {
+         const int attr = mesh->GetAttribute(i);
+         if (attr_marker[attr-1] == 0) { continue; }
+
          T = fes[0]->GetElementTransformation(i);
          for (int s = 0; s < fes.Size(); ++s)
          {
@@ -901,6 +1058,9 @@ void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx) const
 
          for (int k = 0; k < dnfi.Size(); ++k)
          {
+            if (dnfi_marker[k] &&
+                (*dnfi_marker[k])[attr-1] == 0) { continue; }
+
             dnfi[k]->AssembleElementGrad(fe, *T, el_x_const, elmats);
 
             for (int j=0; j<fes.Size(); ++j)
@@ -923,7 +1083,6 @@ void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx) const
    if (fnfi.Size())
    {
       FaceElementTransformations *tr;
-      Mesh *mesh = fes[0]->GetMesh();
 
       for (int i = 0; i < mesh->GetNumFaces(); ++i)
       {
@@ -960,7 +1119,6 @@ void BlockNonlinearForm::ComputeGradientBlocked(const BlockVector &bx) const
    if (bfnfi.Size())
    {
       FaceElementTransformations *tr;
-      Mesh *mesh = fes[0]->GetMesh();
 
       // Which boundary attributes need to be processed?
       Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?

fem/nonlinearform.hpp

@@ -37,6 +37,7 @@ protected:
 
    /// Set of Domain Integrators to be assembled (added).
    Array<NonlinearFormIntegrator*> dnfi; // owned
+   Array<Array<int>*>              dnfi_marker; // not owned
 
    /// Set of interior face Integrators to be assembled (added).
    Array<NonlinearFormIntegrator*> fnfi; // owned
@@ -108,7 +109,12 @@ public:
 
    /// Adds new Domain Integrator.
    void AddDomainIntegrator(NonlinearFormIntegrator *nlfi)
-   { dnfi.Append(nlfi); }
+   { dnfi.Append(nlfi); dnfi_marker.Append(NULL); }
+
+   /// Adds new Domain Integrator, restricted to specific attributes.
+   void AddDomainIntegrator(NonlinearFormIntegrator *nlfi,
+                            Array<int> &elem_marker)
+   { dnfi.Append(nlfi); dnfi_marker.Append(&elem_marker); }
 
    /// Access all integrators added with AddDomainIntegrator().
    Array<NonlinearFormIntegrator*> *GetDNFI() { return &dnfi; }
@@ -227,13 +233,14 @@ protected:
 
    /// Set of Domain Integrators to be assembled (added).
    Array<BlockNonlinearFormIntegrator*> dnfi;
+   Array<Array<int>*>                   dnfi_marker;
 
    /// Set of interior face Integrators to be assembled (added).
    Array<BlockNonlinearFormIntegrator*> fnfi;
 
    /// Set of Boundary Face Integrators to be assembled (added).
    Array<BlockNonlinearFormIntegrator*> bfnfi;
-   Array<Array<int>*>           bfnfi_marker;
+   Array<Array<int>*>                   bfnfi_marker;
 
    /** Auxiliary block-vectors for wrapping input and output vectors or holding
        GridFunction-like block-vector data (e.g. in parallel). */
@@ -298,7 +305,12 @@ public:
 
    /// Adds new Domain Integrator.
    void AddDomainIntegrator(BlockNonlinearFormIntegrator *nlfi)
-   { dnfi.Append(nlfi); }
+   { dnfi.Append(nlfi); dnfi_marker.Append(NULL); }
+
+   /// Adds new Domain Integrator, restricted to specific attributes.
+   void AddDomainIntegrator(BlockNonlinearFormIntegrator *nlfi,
+                            Array<int> &elem_marker)
+   { dnfi.Append(nlfi); dnfi_marker.Append(&elem_marker); }
 
    /// Adds new Interior Face Integrator.
    void AddInteriorFaceIntegrator(BlockNonlinearFormIntegrator *nlfi)
@@ -311,7 +323,8 @@ public:
    /** @brief Adds new Boundary Face Integrator, restricted to specific boundary
        attributes. */
    void AddBdrFaceIntegrator(BlockNonlinearFormIntegrator *nlfi,
-                             Array<int> &bdr_marker);
+                             Array<int> &bdr_marker)
+   { bfnfi.Append(nlfi); bfnfi_marker.Append(&bdr_marker); }
 
    virtual void SetEssentialBC(const Array<Array<int> *>&bdr_attr_is_ess,
                                Array<Vector *> &rhs);

fem/pfespace.cpp

@@ -466,53 +466,54 @@ void ParFiniteElementSpace::ApplyLDofSigns(Table &el_dof) const
    ApplyLDofSigns(all_dofs);
 }
 
-DofTransformation *
-ParFiniteElementSpace::GetElementDofs(int i, Array<int> &dofs) const
+void ParFiniteElementSpace::GetElementDofs(int i, Array<int> &dofs,
+                                           DofTransformation &doftrans) const
 {
    if (elem_dof)
    {
       elem_dof->GetRow(i, dofs);
 
-      if (DoFTrans[mesh->GetElementBaseGeometry(i)])
+      if (DoFTransArray[mesh->GetElementBaseGeometry(i)])
       {
          Array<int> Fo;
          elem_fos->GetRow(i, Fo);
-         DoFTrans[mesh->GetElementBaseGeometry(i)]->SetFaceOrientations(Fo);
-         return DoFTrans[mesh->GetElementBaseGeometry(i)];
+         doftrans.SetDofTransformation(
+            *DoFTransArray[mesh->GetElementBaseGeometry(i)]);
+         doftrans.SetFaceOrientations(Fo);
+         doftrans.SetVDim();
       }
-      return NULL;
+      return;
    }
-   DofTransformation * doftrans = FiniteElementSpace::GetElementDofs(i, dofs);
+   FiniteElementSpace::GetElementDofs(i, dofs, doftrans);
    if (Conforming())
    {
       ApplyLDofSigns(dofs);
    }
-   return doftrans;
 }
 
-DofTransformation *
-ParFiniteElementSpace::GetBdrElementDofs(int i, Array<int> &dofs) const
+void ParFiniteElementSpace::GetBdrElementDofs(int i, Array<int> &dofs,
+                                              DofTransformation &doftrans) const
 {
    if (bdr_elem_dof)
    {
       bdr_elem_dof->GetRow(i, dofs);
 
-      if (DoFTrans[mesh->GetBdrElementBaseGeometry(i)])
+      if (DoFTransArray[mesh->GetBdrElementBaseGeometry(i)])
       {
          Array<int> Fo;
-         bdr_elem_fos -> GetRow (i, Fo);
-         DoFTrans[mesh->GetBdrElementBaseGeometry(i)]->SetFaceOrientations(Fo);
-         return DoFTrans[mesh->GetBdrElementBaseGeometry(i)];
+         bdr_elem_fos->GetRow(i, Fo);
+         doftrans.SetDofTransformation(
+            *DoFTransArray[mesh->GetBdrElementBaseGeometry(i)]);
+         doftrans.SetFaceOrientations(Fo);
+         doftrans.SetVDim();
       }
-      return NULL;
+      return;
    }
-   DofTransformation * doftrans =
-      FiniteElementSpace::GetBdrElementDofs(i, dofs);
+   FiniteElementSpace::GetBdrElementDofs(i, dofs, doftrans);
    if (Conforming())
    {
       ApplyLDofSigns(dofs);
    }
-   return doftrans;
 }
 
 int ParFiniteElementSpace::GetFaceDofs(int i, Array<int> &dofs,
@@ -939,8 +940,8 @@ void ParFiniteElementSpace::Build_Dof_TrueDof_Matrix() const // matrix P
          }
          else if (i_offd[i+1] == i_offd[i] + 2)
          {
-            const double * T = ND_StatelessDofTransformation
-                               ::GetFaceTransform(ltori[i]).GetData();
+            const double *T =
+               ND_DofTransformation::GetFaceTransform(ltori[i]).GetData();
             j_offd[i_offd[i] + 1] = j_offd[i_offd[i]] + 1;
             d_offd[i_offd[i]] = T[0]; d_offd[i_offd[i] + 1] = T[2];
             i++;
@@ -1454,31 +1455,30 @@ void ParFiniteElementSpace::ExchangeFaceNbrData()
    delete [] requests;
 }
 
-DofTransformation *ParFiniteElementSpace::GetFaceNbrElementVDofs(
-   int i, Array<int> &vdofs) const
+void ParFiniteElementSpace::GetFaceNbrElementVDofs(
+   int i, Array<int> &vdofs, DofTransformation &doftrans) const
 {
    face_nbr_element_dof.GetRow(i, vdofs);
 
-   DofTransformation *doftrans = NULL;
-   Geometry::Type geom = GetFaceNbrFE(i)->GetGeomType();
-   if (DoFTrans[geom])
+   if (DoFTransArray[GetFaceNbrFE(i)->GetGeomType()])
    {
       Array<int> F, Fo;
       pmesh->GetFaceNbrElementFaces(pmesh->GetNE() + i, F, Fo);
-      doftrans = DoFTrans[geom];
-      doftrans->SetFaceOrientations(Fo);
-   }
-   if (vdim == 1 || doftrans == NULL)
-   {
-      return doftrans;
-   }
-   else
-   {
-      VDoFTrans.SetDofTransformation(*doftrans);
-      return &VDoFTrans;
+      doftrans.SetDofTransformation(
+         *DoFTransArray[GetFaceNbrFE(i)->GetGeomType()]);
+      doftrans.SetFaceOrientations(Fo);
+      doftrans.SetVDim(vdim, ordering);
    }
 }
 
+DofTransformation *ParFiniteElementSpace::GetFaceNbrElementVDofs(
+   int i, Array<int> &vdofs) const
+{
+   DoFTrans.SetDofTransformation(NULL);
+   GetFaceNbrElementVDofs(i, vdofs, DoFTrans);
+   return DoFTrans.GetDofTransformation() ? &DoFTrans : NULL;
+}
+
 void ParFiniteElementSpace::GetFaceNbrFaceVDofs(int i, Array<int> &vdofs) const
 {
    // Works for NC mesh where 'i' is an index returned by
@@ -2235,19 +2235,13 @@ void NeighborRowMessage::Decode(int rank)
 
             // This is the second "fundamental unit" used in the transformation.
             const auto initial_second_row = second_row;
+            const double *T =
+               ND_DofTransformation::GetFaceTransform(fo).GetData();
 
-            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.AddRow(initial_first_row, T[0] - 1.0);
+            first_row.AddRow(initial_second_row, T[2]);
+            second_row.AddRow(initial_first_row, T[1]);
+            second_row.AddRow(initial_second_row, T[3] - 1.0);
 
             first_row.Collapse();
             second_row.Collapse();

fem/pfespace.hpp

@@ -284,11 +284,17 @@ public:
    /// Return the number of local vector true dofs.
    int GetTrueVSize() const override { return ltdof_size; }
 
-   /// Returns indexes of degrees of freedom in array dofs for i'th element.
-   DofTransformation *GetElementDofs(int i, Array<int> &dofs) const override;
+   /// Returns indexes of degrees of freedom in array dofs for i'th element and
+   /// returns the DofTransformation data in a user-provided object.
+   using FiniteElementSpace::GetElementDofs;
+   void GetElementDofs(int i, Array<int> &dofs,
+                       DofTransformation &doftrans) const override;
 
-   /// Returns indexes of degrees of freedom for i'th boundary element.
-   DofTransformation *GetBdrElementDofs(int i, Array<int> &dofs) const override;
+   /// Returns indexes of degrees of freedom for i'th boundary element and
+   /// returns the DofTransformation data in a user-provided object.
+   using FiniteElementSpace::GetBdrElementDofs;
+   void GetBdrElementDofs(int i, Array<int> &dofs,
+                          DofTransformation &doftrans) 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. */
@@ -382,6 +388,8 @@ public:
    // Face-neighbor functions
    void ExchangeFaceNbrData();
    int GetFaceNbrVSize() const { return num_face_nbr_dofs; }
+   void GetFaceNbrElementVDofs(int i, Array<int> &vdofs,
+                               DofTransformation &doftrans) const;
    DofTransformation *GetFaceNbrElementVDofs(int i, Array<int> &vdofs) const;
    void GetFaceNbrFaceVDofs(int i, Array<int> &vdofs) const;
    const FiniteElement *GetFaceNbrFE(int i) const;

fem/pgridfunc.cpp

@@ -693,7 +693,23 @@ void ParGridFunction::ProjectBdrCoefficient(
 
 #ifdef MFEM_DEBUG
    Array<int> ess_vdofs_marker;
-   pfes->GetEssentialVDofs(attr, ess_vdofs_marker);
+   if (vcoeff) { pfes->GetEssentialVDofs(attr, ess_vdofs_marker); }
+   else
+   {
+      ess_vdofs_marker.SetSize(Size());
+      ess_vdofs_marker = 0;
+      for (int i = 0; i < fes->GetVDim(); i++)
+      {
+         if (!coeff[i]) { continue; }
+         Array<int> component_dof_marker;
+         pfes->GetEssentialVDofs(attr, component_dof_marker,i);
+         for (int j = 0; j<Size(); j++)
+         {
+            ess_vdofs_marker[j] = bool(ess_vdofs_marker[j]) ||
+                                  bool(component_dof_marker[j]);
+         }
+      }
+   }
    for (int i = 0; i < values_counter.Size(); i++)
    {
       MFEM_ASSERT(pfes->GetLocalTDofNumber(i) == -1 ||

fem/qspace.cpp

@@ -177,7 +177,7 @@ int FaceQuadratureSpace::GetEntityIndex(const ElementTransformation &T) const
          return get_face_index(T.ElementNo);
       case ElementTransformation::BDR_ELEMENT:
       case ElementTransformation::BDR_FACE:
-         return get_face_index(mesh.GetBdrElementEdgeIndex(T.ElementNo));
+         return get_face_index(mesh.GetBdrElementFaceIndex(T.ElementNo));
       default:
          MFEM_ABORT("Invalid element type.");
          return -1;

fem/tmop.cpp

@@ -462,21 +462,52 @@ void TMOP_Metric_009::AssembleH(const DenseMatrix &Jpt,
    ie.Assemble_ddI1b(weight, A.GetData());
 }
 
-// mu_14 = |T-I|^2
+double TMOP_Metric_014::EvalWMatrixForm(const DenseMatrix &Jpt) const
+{
+   // mu_14 = |J - I|^2.
+   DenseMatrix Mat(Jpt);
+   Mat(0,0) -= 1.0;
+   Mat(1,1) -= 1.0;
+   return Mat.FNorm2();
+}
+
 double TMOP_Metric_014::EvalW(const DenseMatrix &Jpt) const
 {
-   MFEM_VERIFY(Jtr != NULL,
-               "Requires a target Jacobian, use SetTargetJacobian().");
+   // mu_14 = |J - I|^2 = I1[J-I].
+   DenseMatrix Mat(Jpt);
+   Mat(0,0) -= 1.0;
+   Mat(1,1) -= 1.0;
 
-   DenseMatrix Id(2,2);
+   ie.SetJacobian(Mat.GetData());
+   return ie.Get_I1();
+}
 
-   Id(0,0) = 1; Id(0,1) = 0;
-   Id(1,0) = 0; Id(1,1) = 1;
+void TMOP_Metric_014::EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
+{
+   // P = dI1[J-I] d/dJ[J-I] = dI1[J-I].
+   DenseMatrix JptMinusId = Jpt;
+   for (int i = 0; i < Jpt.Size(); i++)
+   {
+      JptMinusId(i, i) -= 1.0;
+   }
+   ie.SetJacobian(JptMinusId.GetData());
+   P = ie.Get_dI1();
+}
 
-   DenseMatrix Mat(2,2);
-   Mat = Jpt;
-   Mat.Add(-1,Id);
-   return Mat.FNorm2();
+void TMOP_Metric_014::AssembleH(const DenseMatrix &Jpt,
+                                const DenseMatrix &DS,
+                                const double weight,
+                                DenseMatrix &A) const
+{
+   // dP = ddI1[J-I].
+   DenseMatrix JptMinusId = Jpt;
+   for (int i = 0; i < Jpt.Size(); i++)
+   {
+      JptMinusId(i, i) -= 1.0;
+   }
+   ie.SetJacobian(JptMinusId.GetData());
+   ie.SetDerivativeMatrix(DS.Height(), DS.GetData());
+   ie.Assemble_ddI1(weight, A.GetData());
 }
 
 double TMOP_Metric_022::EvalW(const DenseMatrix &Jpt) const

fem/tmop.hpp

@@ -373,16 +373,20 @@ public:
 /// 2D non-barrier Shape+Size+Orientation (VOS) metric (polyconvex).
 class TMOP_Metric_014 : public TMOP_QualityMetric
 {
+protected:
+   mutable InvariantsEvaluator2D<double> ie;
+
 public:
-   // W = |T-I|^2.
+   // W = |J - I|^2.
+   virtual double EvalWMatrixForm(const DenseMatrix &Jpt) const;
+
+   // W = I1[J-I].
    virtual double EvalW(const DenseMatrix &Jpt) const;
 
-   virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const
-   { MFEM_ABORT("Not implemented"); }
+   virtual void EvalP(const DenseMatrix &Jpt, DenseMatrix &P) const;
 
    virtual void AssembleH(const DenseMatrix &Jpt, const DenseMatrix &DS,
-                          const double weight, DenseMatrix &A) const
-   { MFEM_ABORT("Not implemented"); }
+                          const double weight, DenseMatrix &A) const;
 };
 
 /// 2D Shifted barrier form of shape metric (mu_2).

general/communication.cpp

@@ -26,6 +26,10 @@
 #include "sort_pairs.hpp"
 #include "globals.hpp"
 
+#ifdef MFEM_USE_STRUMPACK
+#include <StrumpackConfig.hpp> // STRUMPACK_USE_PTSCOTCH, etc.
+#endif
+
 #include <iostream>
 #include <map>
 
@@ -34,6 +38,14 @@ using namespace std;
 namespace mfem
 {
 
+#if defined(MFEM_USE_STRUMPACK) && \
+    (defined(STRUMPACK_USE_PTSCOTCH) || defined(STRUMPACK_USE_SLATE_SCALAPACK))
+int Mpi::default_thread_required = MPI_THREAD_MULTIPLE;
+#else
+int Mpi::default_thread_required = MPI_THREAD_SINGLE;
+#endif
+
+
 GroupTopology::GroupTopology(const GroupTopology &gt)
    : MyComm(gt.MyComm),
      group_lproc(gt.group_lproc)

general/communication.hpp

@@ -22,7 +22,6 @@
 #include "globals.hpp"
 #include <mpi.h>
 
-
 namespace mfem
 {
 
@@ -32,10 +31,34 @@ namespace mfem
 class Mpi
 {
 public:
-   /// Singleton creation with Mpi::Init();
-   static void Init() { Init_(NULL, NULL); }
-   /// Singleton creation with Mpi::Init(argc,argv);
-   static void Init(int &argc, char **&argv) { Init_(&argc, &argv); }
+   /// Singleton creation with Mpi::Init(argc, argv).
+   static void Init(int &argc, char **&argv,
+                    int required = default_thread_required,
+                    int *provided = nullptr)
+   { Init(&argc, &argv, required, provided); }
+   /// Singleton creation with Mpi::Init().
+   static void Init(int *argc = nullptr, char ***argv = nullptr,
+                    int required = default_thread_required,
+                    int *provided = nullptr)
+   {
+      MFEM_VERIFY(!IsInitialized(), "MPI already initialized!");
+      if (required == MPI_THREAD_SINGLE)
+      {
+         int mpi_err = MPI_Init(argc, argv);
+         MFEM_VERIFY(!mpi_err, "error in MPI_Init()!");
+         if (provided) { *provided = MPI_THREAD_SINGLE; }
+      }
+      else
+      {
+         int mpi_provided;
+         int mpi_err = MPI_Init_thread(argc, argv, required, &mpi_provided);
+         MFEM_VERIFY(!mpi_err, "error in MPI_Init()!");
+         if (provided) { *provided = mpi_provided; }
+      }
+      // The Mpi singleton object below needs to be created after MPI_Init() for
+      // some MPI implementations.
+      Singleton();
+   }
    /// Finalize MPI (if it has been initialized and not yet already finalized).
    static void Finalize()
    {
@@ -71,20 +94,19 @@ public:
    }
    /// Return true if the rank in MPI_COMM_WORLD is zero.
    static bool Root() { return WorldRank() == 0; }
+   /// Default level of thread support for MPI_Init_thread.
+   static MFEM_EXPORT int default_thread_required;
 private:
-   /// Initialize MPI
-   static void Init_(int *argc, char ***argv)
+   /// Initialize the Mpi singleton.
+   static Mpi &Singleton()
    {
-      MFEM_VERIFY(!IsInitialized(), "MPI already initialized!")
-      MPI_Init(argc, argv);
-      // The "mpi" object below needs to be created after MPI_Init() for some
-      // MPI implementations
       static Mpi mpi;
+      return mpi;
    }
-   /// Finalize MPI
+   /// Finalize MPI.
    ~Mpi() { Finalize(); }
-   /// Prevent direct construction of objects of this class
-   Mpi() { }
+   /// Prevent direct construction of objects of this class.
+   Mpi() {}
 };
 
 /** @brief A simple convenience class based on the Mpi singleton class above.

general/isockstream.cpp

@@ -16,13 +16,13 @@
 #include <cstdlib>
 #include <errno.h>
 #ifndef _WIN32
-#include <netinet/in.h>
 #include <netdb.h>
 #include <sys/types.h>
 #include <sys/socket.h>
 #include <unistd.h>
 #else
-#include <winsock.h>
+#include <winsock2.h>
+#include <ws2tcpip.h>
 #ifdef _MSC_VER
 typedef int ssize_t;
 // Link with ws2_32.lib
@@ -51,47 +51,66 @@ int isockstream::establish()
 {
    // char myname[129];
    char   myname[] = "localhost";
-   int    port;
-   struct sockaddr_in sa;
-   struct hostent *hp;
+   int    sfd;
+   struct addrinfo hints, *res, *rp;
 
-   memset(&sa, 0, sizeof(struct sockaddr_in));
-   // gethostname(myname, 128);
-   hp= gethostbyname(myname);
+   memset(&hints, 0, sizeof(hints));
+   hints.ai_family = AF_UNSPEC;
+   hints.ai_socktype = SOCK_STREAM;
+   hints.ai_protocol = 0;
 
-   if (hp == NULL)
+   int s = getaddrinfo(myname, NULL, &hints, &res);
+   if (s != 0)
    {
-      mfem::err << "isockstream::establish(): gethostbyname() failed!\n"
-                << "isockstream::establish(): gethostname() returned: '"
+      mfem::err << "isockstream::establish(): getaddrinfo() failed!\n"
+                << "isockstream::establish(): getaddrinfo() returned: '"
                 << myname << "'" << endl;
       error = 1;
       return (-1);
    }
 
-   sa.sin_family= hp->h_addrtype;
-   sa.sin_port= htons(portnum);
-
-   if ((port = socket(AF_INET, SOCK_STREAM, 0)) < 0)
+   // loop the list of address structures returned by getaddrinfo()
+   for (rp = res; rp != NULL; rp = rp->ai_next)
    {
-      mfem::err << "isockstream::establish(): socket() failed!" << endl;
-      error = 2;
+      if ((sfd = socket(rp->ai_family, rp->ai_socktype, rp->ai_protocol)) < 0)
+      {
+         mfem::err << "isockstream::establish(): socket() failed!" << endl;
+         error = 2;
+         return (-1);
+      }
+
+      int on = 1;
+      if (setsockopt(sfd, SOL_SOCKET, SO_REUSEADDR, (char *)&on, sizeof(on)) < 0)
+      {
+         mfem::err << "isockstream::establish(): setsockopt() failed!" << endl;
+         return (-1);
+      }
+
+#if defined(__APPLE__)
+      if (bind(sfd, (const struct sockaddr *)rp->ai_addr, rp->ai_addrlen) < 0)
+#else
+      if (bind(sfd, rp->ai_addr, rp->ai_addrlen) < 0)
+#endif
+      {
+         mfem::err << "isockstream::establish(): bind() failed!" << endl;
+         close(sfd);
+         error = 3;
+         continue;
+      }
+
+      break;
+   }
+
+   // No address succeeded
+   if (rp == NULL)
+   {
+      mfem::err << "Could not bind\n";
       return (-1);
    }
 
-   int on=1;
-   setsockopt(port, SOL_SOCKET, SO_REUSEADDR, (char *)(&on), sizeof(on));
-
-   if (bind(port,(const sockaddr*)&sa,(socklen_t)sizeof(struct sockaddr_in)) < 0)
-   {
-      mfem::err << "isockstream::establish(): bind() failed!" << endl;
-      close(port);
-      error = 3;
-      return (-1);
-   }
-
-   listen(port, 4);
-   error = 0;
-   return (port);
+   freeaddrinfo(res);
+   listen(sfd, 4);
+   return (sfd);
 }
 
 int isockstream::read_data(int s, char *buf, int n)

general/socketstream.cpp

@@ -19,15 +19,15 @@
 #include <cstring>      // memset, memcpy, strerror
 #include <cerrno>       // errno
 #ifndef _WIN32
-#include <netdb.h>      // gethostbyname
+#include <netdb.h>      // getaddrinfo
 #include <arpa/inet.h>  // htons
 #include <sys/types.h>  // socket, setsockopt, connect, recv, send
 #include <sys/socket.h> // socket, setsockopt, connect, recv, send
 #include <unistd.h>     // close
-#include <netinet/in.h> // sockaddr_in
 #define closesocket (::close)
 #else
-#include <winsock.h>
+#include <winsock2.h>
+#include <ws2tcpip.h>
 #ifdef _MSC_VER
 typedef int ssize_t;
 // Link with ws2_32.lib
@@ -93,8 +93,7 @@ int socketbuf::attach(int sd)
 
 int socketbuf::open(const char hostname[], int port)
 {
-   struct sockaddr_in  sa;
-   struct hostent     *hp;
+   struct addrinfo     hints, *res, *rp;
 
    if (!wsInit_.Initialized())
    {
@@ -105,42 +104,47 @@ int socketbuf::open(const char hostname[], int port)
    setg(NULL, NULL, NULL);
    setp(obuf, obuf + buflen);
 
-   hp = gethostbyname(hostname);
-   if (hp == NULL)
+   hints.ai_family = AF_UNSPEC;
+   hints.ai_socktype = SOCK_STREAM;
+   hints.ai_protocol = 0;
+
+   int s = getaddrinfo(hostname, NULL, &hints, &res);
+   if (s != 0)
    {
       socket_descriptor = -3;
       return -1;
    }
-   memset(&sa, 0, sizeof(sa));
-   memcpy((char *)&sa.sin_addr, hp->h_addr, hp->h_length);
-   sa.sin_family = hp->h_addrtype;
-   sa.sin_port = htons(port);
-   socket_descriptor = socket(hp->h_addrtype, SOCK_STREAM, 0);
-   if (socket_descriptor < 0)
+
+   for (rp = res; rp != NULL; rp = rp->ai_next)
    {
-      return -1;
-   }
+      socket_descriptor = socket(rp->ai_family, rp->ai_socktype, rp->ai_protocol);
+      if (socket_descriptor < 0)
+      {
+         continue;
+      }
 
 #if defined __APPLE__
-   // OS X does not support the MSG_NOSIGNAL option of send().
-   // Instead we can use the SO_NOSIGPIPE socket option.
-   int on = 1;
-   if (setsockopt(socket_descriptor, SOL_SOCKET, SO_NOSIGPIPE,
-                  (char *)(&on), sizeof(on)) < 0)
-   {
-      closesocket(socket_descriptor);
-      socket_descriptor = -2;
-      return -1;
-   }
+      // OS X does not support the MSG_NOSIGNAL option of send().
+      // Instead we can use the SO_NOSIGPIPE socket option.
+      int on = 1;
+      if (setsockopt(socket_descriptor, SOL_SOCKET, SO_NOSIGPIPE,
+                     &on, sizeof(on)) < 0)
+      {
+         closesocket(socket_descriptor);
+         socket_descriptor = -2;
+         return -1;
+      }
 #endif
 
-   if (connect(socket_descriptor,
-               (const struct sockaddr *)&sa, sizeof(sa)) < 0)
-   {
-      closesocket(socket_descriptor);
-      socket_descriptor = -2;
-      return -1;
+      if (connect(socket_descriptor, rp->ai_addr, rp->ai_addrlen) < 0)
+      {
+         closesocket(socket_descriptor);
+         socket_descriptor = -2;
+         continue;
+      }
    }
+
+   freeaddrinfo(res);
    return 0;
 }
 

linalg/strumpack.cpp

@@ -16,238 +16,471 @@
 
 #include "strumpack.hpp"
 
-using namespace std;
-using namespace strumpack;
-
 namespace mfem
 {
 
 STRUMPACKRowLocMatrix::STRUMPACKRowLocMatrix(MPI_Comm comm,
-                                             int num_loc_rows, int first_loc_row,
-                                             int glob_nrows, int glob_ncols,
-                                             int *I, int *J, double *data)
-   : comm_(comm), A_(NULL)
+                                             int num_loc_rows,
+                                             HYPRE_BigInt first_loc_row,
+                                             HYPRE_BigInt glob_nrows,
+                                             HYPRE_BigInt glob_ncols,
+                                             int *I, HYPRE_BigInt *J,
+                                             double *data, bool sym_sparse)
 {
    // Set mfem::Operator member data
    height = num_loc_rows;
    width  = num_loc_rows;
 
-   // Allocate STRUMPACK's CSRMatrixMPI
-   int nprocs, rank;
-   MPI_Comm_rank(comm_, &rank);
-   MPI_Comm_size(comm_, &nprocs);
-   int * dist = new int[nprocs + 1];
-   dist[rank + 1] = first_loc_row + num_loc_rows;
+   // Allocate STRUMPACK's CSRMatrixMPI (copies all inputs)
+   int rank, nprocs;
+   MPI_Comm_rank(comm, &rank);
+   MPI_Comm_size(comm, &nprocs);
+   Array<HYPRE_BigInt> dist(nprocs + 1);
    dist[0] = 0;
-   MPI_Allgather(MPI_IN_PLACE, 0, MPI_INT, dist + 1, 1, MPI_INT, comm_);
-   A_ = new CSRMatrixMPI<double,int>(num_loc_rows, I, J, data, dist, comm_, false);
-   delete[] dist;
-}
+   dist[rank + 1] = first_loc_row + (HYPRE_BigInt)num_loc_rows;
+   MPI_Allgather(MPI_IN_PLACE, 0, MPI_DATATYPE_NULL,
+                 dist.GetData() + 1, 1, HYPRE_MPI_BIG_INT, comm);
 
-STRUMPACKRowLocMatrix::STRUMPACKRowLocMatrix(const HypreParMatrix & hypParMat)
-   : comm_(hypParMat.GetComm()),
-     A_(NULL)
-{
-   // First cast the parameter to a hypre_ParCSRMatrix
-   hypre_ParCSRMatrix * parcsr_op =
-      (hypre_ParCSRMatrix *)const_cast<HypreParMatrix&>(hypParMat);
-
-   MFEM_ASSERT(parcsr_op != NULL,"STRUMPACK: const_cast failed in SetOperator");
-
-   // Create the CSRMatrixMPI A_ by borrowing the internal data from a
-   // hypre_CSRMatrix.
-   hypParMat.HostRead();
-   hypre_CSRMatrix * csr_op = hypre_MergeDiagAndOffd(parcsr_op);
-   hypParMat.HypreRead();
-   hypre_CSRMatrixSetDataOwner(csr_op,0);
-#if MFEM_HYPRE_VERSION >= 21600
-   // For now, this method assumes that HYPRE_Int is int. Also, csr_op->num_cols
-   // is of type HYPRE_Int, so if we want to check for big indices in
-   // csr_op->big_j, we'll have to check all entries and that check will only be
-   // necessary in HYPRE_MIXEDINT mode which is not supported at the moment.
-   hypre_CSRMatrixBigJtoJ(csr_op);
-#endif
-
-   height = csr_op->num_rows;
-   width  = csr_op->num_rows;
-
-   int nprocs, rank;
-   MPI_Comm_rank(comm_, &rank);
-   MPI_Comm_size(comm_, &nprocs);
-   int * dist = new int[nprocs + 1];
-   dist[rank + 1] = parcsr_op->first_row_index + csr_op->num_rows;
-   dist[0] = 0;
-   MPI_Allgather(MPI_IN_PLACE, 0, MPI_INT, dist + 1, 1, MPI_INT, comm_);
-   A_ = new CSRMatrixMPI<double,int>(csr_op->num_rows, csr_op->i, csr_op->j,
-                                     csr_op->data, dist, comm_, false);
-   delete[] dist;
-
-   // Everything has been copied or abducted so delete the structure
-   hypre_CSRMatrixDestroy(csr_op);
-}
-
-STRUMPACKRowLocMatrix::~STRUMPACKRowLocMatrix()
-{
-   // Delete the struct
-   if ( A_ != NULL ) { delete A_; }
-}
-
-STRUMPACKSolver::STRUMPACKSolver( int argc, char* argv[], MPI_Comm comm )
-   : comm_(comm),
-     APtr_(NULL),
-     solver_(NULL)
-{
-   this->Init(argc, argv);
-}
-
-STRUMPACKSolver::STRUMPACKSolver( STRUMPACKRowLocMatrix & A )
-   : comm_(A.GetComm()),
-     APtr_(&A),
-     solver_(NULL)
-{
-   height = A.Height();
-   width  = A.Width();
-
-   this->Init(0, NULL);
-}
-
-STRUMPACKSolver::~STRUMPACKSolver()
-{
-   if ( solver_ != NULL ) { delete solver_; }
-}
-
-void STRUMPACKSolver::Init( int argc, char* argv[] )
-{
-   MPI_Comm_size(comm_, &numProcs_);
-   MPI_Comm_rank(comm_, &myid_);
-
-   factor_verbose_ = false;
-   solve_verbose_ = false;
-
-   solver_ = new StrumpackSparseSolverMPIDist<double,int>(comm_, argc, argv,
-                                                          false);
-}
-
-void STRUMPACKSolver::SetFromCommandLine( )
-{
-   solver_->options().set_from_command_line( );
-}
-
-void STRUMPACKSolver::SetPrintFactorStatistics( bool print_stat )
-{
-   factor_verbose_ = print_stat;
-}
-
-void STRUMPACKSolver::SetPrintSolveStatistics( bool print_stat )
-{
-   solve_verbose_ = print_stat;
-}
-
-void STRUMPACKSolver::SetKrylovSolver( strumpack::KrylovSolver method )
-{
-   solver_->options().set_Krylov_solver( method );
-}
-
-void STRUMPACKSolver::SetReorderingStrategy( strumpack::ReorderingStrategy
-                                             method )
-{
-   solver_->options().set_reordering_method( method );
-}
-
-void STRUMPACKSolver::DisableMatching( )
-{
-#if STRUMPACK_VERSION_MAJOR >= 3
-   solver_->options().set_matching( strumpack::MatchingJob::NONE );
+#if !(defined(HYPRE_BIGINT) || defined(HYPRE_MIXEDINT))
+   A_ = new strumpack::CSRMatrixMPI<double, HYPRE_BigInt>(
+      (HYPRE_BigInt)num_loc_rows, I, J, data, dist.GetData(),
+      comm, sym_sparse);
 #else
-   solver_->options().set_mc64job( strumpack::MC64Job::NONE );
+   Array<HYPRE_BigInt> II(num_loc_rows+1);
+   for (int i = 0; i <= num_loc_rows; i++) { II[i] = (HYPRE_BigInt)I[i]; }
+   A_ = new strumpack::CSRMatrixMPI<double, HYPRE_BigInt>(
+      (HYPRE_BigInt)num_loc_rows, II.GetData(), J, data, dist.GetData(),
+      comm, sym_sparse);
 #endif
 }
 
-void STRUMPACKSolver::EnableMatching( )
+STRUMPACKRowLocMatrix::STRUMPACKRowLocMatrix(const Operator &op,
+                                             bool sym_sparse)
 {
-#if STRUMPACK_VERSION_MAJOR >= 3
-   solver_->options().set_matching
-   ( strumpack::MatchingJob::MAX_DIAGONAL_PRODUCT_SCALING );
-#else
-   solver_->options().set_mc64job
-   ( strumpack::MC64Job::MAX_DIAGONAL_PRODUCT_SCALING );
-#endif
-}
-
-#if STRUMPACK_VERSION_MAJOR >= 3
-void STRUMPACKSolver::EnableParallelMatching( )
-{
-   solver_->options().set_matching
-   ( strumpack::MatchingJob::COMBBLAS );
-}
-#endif
-
-void STRUMPACKSolver::SetRelTol( double rtol )
-{
-   solver_->options().set_rel_tol( rtol );
-}
-
-void STRUMPACKSolver::SetAbsTol( double atol )
-{
-   solver_->options().set_abs_tol( atol );
-}
-
-
-void STRUMPACKSolver::Mult( const Vector & x, Vector & y ) const
-{
-   MFEM_ASSERT(APtr_ != NULL,
-               "STRUMPACK Error: The operator must be set before"
-               " the system can be solved.");
-   MFEM_ASSERT(x.Size() == Width(), "invalid x.Size() = " << x.Size()
-               << ", expected size = " << Width());
-   MFEM_ASSERT(y.Size() == Height(), "invalid y.Size() = " << y.Size()
-               << ", expected size = " << Height());
-
-   double*  yPtr = y.HostWrite();
-   const double*  xPtr = x.HostRead();
-
-   solver_->options().set_verbose( factor_verbose_ );
-   ReturnCode ret = solver_->factor();
-   switch (ret)
-   {
-      case ReturnCode::SUCCESS: break;
-      case ReturnCode::MATRIX_NOT_SET:
-      {
-         MFEM_ABORT("STRUMPACK:  Matrix was not set!");
-      }
-      break;
-      case ReturnCode::REORDERING_ERROR:
-      {
-         MFEM_ABORT("STRUMPACK:  Matrix reordering failed!");
-      }
-      break;
-      default:
-      {
-         MFEM_ABORT("STRUMPACK: 'factor()' error code = " << ret);
-      }
-   }
-   solver_->options().set_verbose( solve_verbose_ );
-   solver_->solve(xPtr, yPtr);
-
-}
-
-void STRUMPACKSolver::SetOperator( const Operator & op )
-{
-   // Verify that we have a compatible operator
-   APtr_ = dynamic_cast<const STRUMPACKRowLocMatrix*>(&op);
-   if ( APtr_ == NULL )
-   {
-      mfem_error("STRUMPACKSolver::SetOperator : not STRUMPACKRowLocMatrix!");
-   }
-
-   solver_->set_matrix( *(APtr_->getA()) );
+   const HypreParMatrix *APtr = dynamic_cast<const HypreParMatrix *>(&op);
+   MFEM_VERIFY(APtr, "Not a compatible matrix type");
+   MPI_Comm comm = APtr->GetComm();
 
    // Set mfem::Operator member data
    height = op.Height();
    width  = op.Width();
 
+   // First cast the parameter to a hypre_ParCSRMatrix
+   hypre_ParCSRMatrix *parcsr_op =
+      (hypre_ParCSRMatrix *)const_cast<HypreParMatrix &>(*APtr);
+
+   // Create the CSRMatrixMPI A by taking the internal data from a
+   // hypre_CSRMatrix
+   APtr->HostRead();
+   hypre_CSRMatrix *csr_op = hypre_MergeDiagAndOffd(parcsr_op);
+   APtr->HypreRead();
+   HYPRE_Int       *Iptr   = csr_op->i;
+#if MFEM_HYPRE_VERSION >= 21600
+   HYPRE_BigInt    *Jptr   = csr_op->big_j;
+#else
+   HYPRE_Int       *Jptr   = csr_op->j;
+#endif
+   double          *data   = csr_op->data;
+
+   HYPRE_BigInt fst_row = parcsr_op->first_row_index;
+   HYPRE_Int    m_loc   = csr_op->num_rows;
+
+   // Allocate STRUMPACK's CSRMatrixMPI
+   int rank, nprocs;
+   MPI_Comm_rank(comm, &rank);
+   MPI_Comm_size(comm, &nprocs);
+   Array<HYPRE_BigInt> dist(nprocs + 1);
+   dist[0] = 0;
+   dist[rank + 1] = fst_row + (HYPRE_BigInt)m_loc;
+   MPI_Allgather(MPI_IN_PLACE, 0, MPI_DATATYPE_NULL,
+                 dist.GetData() + 1, 1, HYPRE_MPI_BIG_INT, comm);
+
+#if !defined(HYPRE_MIXEDINT)
+   A_ = new strumpack::CSRMatrixMPI<double, HYPRE_BigInt>(
+      (HYPRE_BigInt)m_loc, Iptr, Jptr, data, dist.GetData(),
+      comm, sym_sparse);
+#else
+   Array<HYPRE_BigInt> II(m_loc+1);
+   for (int i = 0; i <= m_loc; i++) { II[i] = (HYPRE_BigInt)Iptr[i]; }
+   A_ = new strumpack::CSRMatrixMPI<double, HYPRE_BigInt>(
+      (HYPRE_BigInt)m_loc, II.GetData(), Jptr, data, dist.GetData(),
+      comm, sym_sparse);
+#endif
+
+   // Everything has been copied so delete the structure
+   hypre_CSRMatrixDestroy(csr_op);
 }
 
+STRUMPACKRowLocMatrix::~STRUMPACKRowLocMatrix()
+{
+   delete A_;
+}
+
+template <typename STRUMPACKSolverType>
+STRUMPACKSolverBase<STRUMPACKSolverType>::
+STRUMPACKSolverBase(MPI_Comm comm, int argc, char *argv[])
+   : APtr_(NULL),
+     factor_verbose_(false),
+     solve_verbose_(false),
+     reorder_reuse_(false),
+     nrhs_(-1)
+{
+   solver_ = new STRUMPACKSolverType(comm, argc, argv, false);
+}
+
+template <typename STRUMPACKSolverType>
+STRUMPACKSolverBase<STRUMPACKSolverType>::
+STRUMPACKSolverBase(STRUMPACKRowLocMatrix &A, int argc, char *argv[])
+   : APtr_(&A),
+     factor_verbose_(false),
+     solve_verbose_(false),
+     reorder_reuse_(false),
+     nrhs_(-1)
+{
+   solver_ = new STRUMPACKSolverType(A.GetComm(), argc, argv, false);
+   SetOperator(A);
+}
+
+template <typename STRUMPACKSolverType>
+STRUMPACKSolverBase<STRUMPACKSolverType>::
+~STRUMPACKSolverBase()
+{
+   delete solver_;
+}
+
+template <typename STRUMPACKSolverType>
+void STRUMPACKSolverBase<STRUMPACKSolverType>::
+SetFromCommandLine()
+{
+   solver_->options().set_from_command_line();
+}
+
+template <typename STRUMPACKSolverType>
+void STRUMPACKSolverBase<STRUMPACKSolverType>::
+SetPrintFactorStatistics(bool print_stat)
+{
+   factor_verbose_ = print_stat;
+}
+
+template <typename STRUMPACKSolverType>
+void STRUMPACKSolverBase<STRUMPACKSolverType>::
+SetPrintSolveStatistics(bool print_stat)
+{
+   solve_verbose_ = print_stat;
+}
+
+template <typename STRUMPACKSolverType>
+void STRUMPACKSolverBase<STRUMPACKSolverType>
+::SetRelTol(double rtol)
+{
+   solver_->options().set_rel_tol(rtol);
+}
+
+template <typename STRUMPACKSolverType>
+void STRUMPACKSolverBase<STRUMPACKSolverType>
+::SetAbsTol(double atol)
+{
+   solver_->options().set_abs_tol(atol);
+}
+
+template <typename STRUMPACKSolverType>
+void STRUMPACKSolverBase<STRUMPACKSolverType>
+::SetMaxIter(int max_it)
+{
+   solver_->options().set_maxit(max_it);
+}
+
+template <typename STRUMPACKSolverType>
+void STRUMPACKSolverBase<STRUMPACKSolverType>
+::SetReorderingReuse(bool reuse)
+{
+   reorder_reuse_ = reuse;
+}
+
+template <typename STRUMPACKSolverType>
+void STRUMPACKSolverBase<STRUMPACKSolverType>
+::EnableGPU()
+{
+   solver_->options().enable_gpu();
+}
+
+template <typename STRUMPACKSolverType>
+void STRUMPACKSolverBase<STRUMPACKSolverType>
+::DisableGPU()
+{
+   solver_->options().disable_gpu();
+}
+
+template <typename STRUMPACKSolverType>
+void STRUMPACKSolverBase<STRUMPACKSolverType>::
+SetKrylovSolver(strumpack::KrylovSolver method)
+{
+   solver_->options().set_Krylov_solver(method);
+}
+
+template <typename STRUMPACKSolverType>
+void STRUMPACKSolverBase<STRUMPACKSolverType>::
+SetReorderingStrategy(strumpack::ReorderingStrategy method)
+{
+   solver_->options().set_reordering_method(method);
+}
+
+template <typename STRUMPACKSolverType>
+void STRUMPACKSolverBase<STRUMPACKSolverType>::
+SetMatching(strumpack::MatchingJob job)
+{
+   solver_->options().set_matching(job);
+}
+
+template <typename STRUMPACKSolverType>
+void STRUMPACKSolverBase<STRUMPACKSolverType>::
+SetCompression(strumpack::CompressionType type)
+{
+#if STRUMPACK_VERSION_MAJOR >= 5
+   solver_->options().set_compression(type);
+#else
+   switch (type)
+   {
+      case strumpack::NONE:
+         solver_->options().disable_BLR();
+         solver_->options().disable_HSS();
+         break;
+      case strumpack::BLR:
+         solver_->options().enable_BLR();
+         break;
+      case strumpack::HSS:
+         solver_->options().enable_HSS();
+         break;
+      default:
+         MFEM_ABORT("Invalid compression type for STRUMPACK version " <<
+                    STRUMPACK_VERSION_MAJOR << "!");
+         break;
+   }
+#endif
+}
+
+template <typename STRUMPACKSolverType>
+void STRUMPACKSolverBase<STRUMPACKSolverType>::
+SetCompressionRelTol(double rtol)
+{
+#if STRUMPACK_VERSION_MAJOR >= 5
+   solver_->options().set_compression_rel_tol(rtol);
+#else
+   solver_->options().BLR_options().set_rel_tol(rtol);
+   solver_->options().HSS_options().set_rel_tol(rtol);
+#endif
+}
+
+template <typename STRUMPACKSolverType>
+void STRUMPACKSolverBase<STRUMPACKSolverType>::
+SetCompressionAbsTol(double atol)
+{
+#if STRUMPACK_VERSION_MAJOR >= 5
+   solver_->options().set_compression_abs_tol(atol);
+#else
+   solver_->options().BLR_options().set_abs_tol(atol);
+   solver_->options().HSS_options().set_abs_tol(atol);
+#endif
+}
+
+#if STRUMPACK_VERSION_MAJOR >= 5
+template <typename STRUMPACKSolverType>
+void STRUMPACKSolverBase<STRUMPACKSolverType>::
+SetCompressionLossyPrecision(int precision)
+{
+   solver_->options().set_lossy_precision(precision);
+}
+
+template <typename STRUMPACKSolverType>
+void STRUMPACKSolverBase<STRUMPACKSolverType>::
+SetCompressionButterflyLevels(int levels)
+{
+   solver_->options().HODLR_options().set_butterfly_levels(levels);
+}
+#endif
+
+template <typename STRUMPACKSolverType>
+void STRUMPACKSolverBase<STRUMPACKSolverType>::
+SetOperator(const Operator &op)
+{
+   // Verify that we have a compatible operator
+   bool first_mat = !APtr_;
+   APtr_ = dynamic_cast<const STRUMPACKRowLocMatrix *>(&op);
+   MFEM_VERIFY(APtr_,
+               "STRUMPACK: Operator is not a STRUMPACKRowLocMatrix!");
+
+   // Set mfem::Operator member data
+   height = op.Height();
+   width  = op.Width();
+
+   if (first_mat || !reorder_reuse_)
+   {
+      solver_->set_matrix(*(APtr_->GetA()));
+   }
+   else
+   {
+      solver_->update_matrix_values(*(APtr_->GetA()));
+   }
+}
+
+template <typename STRUMPACKSolverType>
+void STRUMPACKSolverBase<STRUMPACKSolverType>::
+FactorInternal() const
+{
+   MFEM_ASSERT(APtr_,
+               "STRUMPACK: Operator must be set before the system can be "
+               "solved!");
+   solver_->options().set_verbose(factor_verbose_);
+   strumpack::ReturnCode ret = solver_->factor();
+   if (ret != strumpack::ReturnCode::SUCCESS)
+   {
+#if STRUMPACK_VERSION_MAJOR >= 7
+      MFEM_ABORT("STRUMPACK: Factor failed with return code " << ret << "!");
+#else
+      MFEM_ABORT("STRUMPACK: Factor failed!");
+#endif
+   }
+}
+
+template <typename STRUMPACKSolverType>
+void STRUMPACKSolverBase<STRUMPACKSolverType>::
+Mult(const Vector &x, Vector &y) const
+{
+   MFEM_ASSERT(x.Size() == Width(),
+               "STRUMPACK: Invalid x.Size() = " << x.Size() <<
+               ", expected size = " << Width() << "!");
+   MFEM_ASSERT(y.Size() == Height(),
+               "STRUMPACK: Invalid y.Size() = " << y.Size() <<
+               ", expected size = " << Height() << "!");
+
+   const double *xPtr = x.HostRead();
+   double *yPtr       = y.HostReadWrite();
+
+   FactorInternal();
+   solver_->options().set_verbose(solve_verbose_);
+   strumpack::ReturnCode ret = solver_->solve(xPtr, yPtr, false);
+   if (ret != strumpack::ReturnCode::SUCCESS)
+   {
+#if STRUMPACK_VERSION_MAJOR >= 7
+      MFEM_ABORT("STRUMPACK: Solve failed with return code " << ret << "!");
+#else
+      MFEM_ABORT("STRUMPACK: Solve failed!");
+#endif
+   }
+}
+
+template <typename STRUMPACKSolverType>
+void STRUMPACKSolverBase<STRUMPACKSolverType>::
+ArrayMult(const Array<const Vector *> &X, Array<Vector *> &Y) const
+{
+   MFEM_ASSERT(X.Size() == Y.Size(),
+               "Number of columns mismatch in STRUMPACK solve!");
+   if (X.Size() == 1)
+   {
+      nrhs_ = 1;
+      MFEM_ASSERT(X[0] && Y[0], "Missing Vector in STRUMPACK solve!");
+      Mult(*X[0], *Y[0]);
+      return;
+   }
+
+   // Multiple RHS case
+   int ldx = Height();
+   if (nrhs_ != X.Size())
+   {
+      rhs_.SetSize(X.Size() * ldx);
+      sol_.SetSize(X.Size() * ldx);
+      nrhs_ = X.Size();
+   }
+   for (int i = 0; i < nrhs_; i++)
+   {
+      MFEM_ASSERT(X[i] && X[i]->Size() == Width(),
+                  "STRUMPACK: Missing or invalid sized RHS Vector in solve!");
+      Vector s(rhs_, i * ldx, ldx);
+      s = *X[i];
+      rhs_.SyncMemory(s);  // Update flags for rhs_ if updated on device
+   }
+   const double *xPtr = rhs_.HostRead();
+   double *yPtr       = sol_.HostReadWrite();
+
+   FactorInternal();
+   solver_->options().set_verbose(solve_verbose_);
+   strumpack::ReturnCode ret = solver_->solve(nrhs_, xPtr, ldx, yPtr, ldx,
+                                              false);
+   if (ret != strumpack::ReturnCode::SUCCESS)
+   {
+#if STRUMPACK_VERSION_MAJOR >= 7
+      MFEM_ABORT("STRUMPACK: Solve failed with return code " << ret << "!");
+#else
+      MFEM_ABORT("STRUMPACK: Solve failed!");
+#endif
+   }
+
+   for (int i = 0; i < nrhs_; i++)
+   {
+      MFEM_ASSERT(Y[i] && Y[i]->Size() == Width(),
+                  "STRUMPACK: Missing or invalid sized solution Vector in solve!");
+      Vector s(sol_, i * ldx, ldx);
+      *Y[i] = s;
+   }
+}
+
+STRUMPACKSolver::
+STRUMPACKSolver(MPI_Comm comm)
+   : STRUMPACKSolverBase<strumpack::
+     SparseSolverMPIDist<double, HYPRE_BigInt>>
+     (comm, 0, NULL) {}
+
+STRUMPACKSolver::
+STRUMPACKSolver(STRUMPACKRowLocMatrix &A)
+   : STRUMPACKSolverBase<strumpack::
+     SparseSolverMPIDist<double, HYPRE_BigInt>>
+     (A, 0, NULL) {}
+
+STRUMPACKSolver::
+STRUMPACKSolver(MPI_Comm comm, int argc, char *argv[])
+   : STRUMPACKSolverBase<strumpack::
+     SparseSolverMPIDist<double, HYPRE_BigInt>>
+     (comm, argc, argv) {}
+
+STRUMPACKSolver::
+STRUMPACKSolver(STRUMPACKRowLocMatrix &A, int argc, char *argv[])
+   : STRUMPACKSolverBase<strumpack::
+     SparseSolverMPIDist<double, HYPRE_BigInt>>
+     (A, argc, argv) {}
+
+#if STRUMPACK_VERSION_MAJOR >= 7
+STRUMPACKMixedPrecisionSolver::
+STRUMPACKMixedPrecisionSolver(MPI_Comm comm)
+   : STRUMPACKSolverBase<strumpack::
+     SparseSolverMixedPrecisionMPIDist<float, double, HYPRE_BigInt>>
+     (comm, 0, NULL) {}
+
+STRUMPACKMixedPrecisionSolver::
+STRUMPACKMixedPrecisionSolver(STRUMPACKRowLocMatrix &A)
+   : STRUMPACKSolverBase<strumpack::
+     SparseSolverMixedPrecisionMPIDist<float, double, HYPRE_BigInt>>
+     (A, 0, NULL) {}
+
+STRUMPACKMixedPrecisionSolver::
+STRUMPACKMixedPrecisionSolver(MPI_Comm comm, int argc, char *argv[])
+   : STRUMPACKSolverBase<strumpack::
+     SparseSolverMixedPrecisionMPIDist<float, double, HYPRE_BigInt>>
+     (comm, argc, argv) {}
+
+STRUMPACKMixedPrecisionSolver::
+STRUMPACKMixedPrecisionSolver(STRUMPACKRowLocMatrix &A, int argc, char *argv[])
+   : STRUMPACKSolverBase<strumpack::
+     SparseSolverMixedPrecisionMPIDist<float, double, HYPRE_BigInt>>
+     (A, argc, argv) {}
+#endif
+
+template class STRUMPACKSolverBase<strumpack::
+                                   SparseSolverMPIDist<double, HYPRE_BigInt>>;
+#if STRUMPACK_VERSION_MAJOR >= 7
+template class STRUMPACKSolverBase<strumpack::
+                                   SparseSolverMixedPrecisionMPIDist<float, double, HYPRE_BigInt>>;
+#endif
+
 } // mfem namespace
 
 #endif // MFEM_USE_MPI

linalg/strumpack.hpp

@@ -16,12 +16,14 @@
 
 #ifdef MFEM_USE_STRUMPACK
 #ifdef MFEM_USE_MPI
+
 #include "operator.hpp"
 #include "hypre.hpp"
-
 #include <mpi.h>
 
+// STRUMPACK headers
 #include "StrumpackSparseSolverMPIDist.hpp"
+#include "StrumpackSparseSolverMixedPrecisionMPIDist.hpp"
 
 namespace mfem
 {
@@ -34,63 +36,80 @@ public:
        be of size (local) nrows by (global) glob_ncols. The new parallel matrix
        contains copies of all input arrays (so they can be deleted). */
    STRUMPACKRowLocMatrix(MPI_Comm comm,
-                         int num_loc_rows, int first_loc_row,
-                         int glob_nrows, int glob_ncols,
-                         int *I, int *J, double *data);
+                         int num_loc_rows, HYPRE_BigInt first_loc_row,
+                         HYPRE_BigInt glob_nrows, HYPRE_BigInt glob_ncols,
+                         int *I, HYPRE_BigInt *J, double *data,
+                         bool sym_sparse = false);
 
    /** Creates a copy of the parallel matrix hypParMat in STRUMPACK's RowLoc
        format. All data is copied so the original matrix may be deleted. */
-   STRUMPACKRowLocMatrix(const HypreParMatrix & hypParMat);
+   STRUMPACKRowLocMatrix(const Operator &op, bool sym_sparse = false);
 
    ~STRUMPACKRowLocMatrix();
 
    void Mult(const Vector &x, Vector &y) const
    {
-      mfem_error("STRUMPACKRowLocMatrix::Mult(...)\n"
-                 "  matrix vector products are not supported.");
+      MFEM_ABORT("STRUMPACKRowLocMatrix::Mult: Matrix vector products are not "
+                 "supported!");
    }
 
-   MPI_Comm GetComm() const { return comm_; }
+   MPI_Comm GetComm() const { return A_->comm(); }
 
-   strumpack::CSRMatrixMPI<double,int>* getA() const { return A_; }
+   strumpack::CSRMatrixMPI<double, HYPRE_BigInt> *GetA() const { return A_; }
 
 private:
-   MPI_Comm   comm_;
-   strumpack::CSRMatrixMPI<double,int>* A_;
-
-}; // mfem::STRUMPACKRowLocMatrix
+   strumpack::CSRMatrixMPI<double, HYPRE_BigInt> *A_;
+};
 
 /** The MFEM STRUMPACK Direct Solver class.
 
     The mfem::STRUMPACKSolver class uses the STRUMPACK library to perform LU
     factorization of a parallel sparse matrix. The solver is capable of handling
-    double precision types. See http://portal.nersc.gov/project/sparse/strumpack
+    double precision types. See
+    http://portal.nersc.gov/project/sparse/strumpack/.
 */
-class STRUMPACKSolver : public mfem::Solver
+template <typename STRUMPACKSolverType>
+class STRUMPACKSolverBase : public Solver
 {
+protected:
+   // Constructor with MPI_Comm parameter and command line arguments.
+   STRUMPACKSolverBase(MPI_Comm comm, int argc, char *argv[]);
+
+   // Constructor with STRUMPACK matrix object and command line arguments.
+   STRUMPACKSolverBase(STRUMPACKRowLocMatrix &A, int argc, char *argv[]);
+
 public:
-   // Constructor with MPI_Comm parameter.
-   STRUMPACKSolver( int argc, char* argv[], MPI_Comm comm );
-
-   // Constructor with STRUMPACK Matrix Object.
-   STRUMPACKSolver( STRUMPACKRowLocMatrix & A);
-
    // Default destructor.
-   ~STRUMPACKSolver( void );
+   virtual ~STRUMPACKSolverBase();
 
    // Factor and solve the linear system y = Op^{-1} x.
-   void Mult( const Vector & x, Vector & y ) const;
+   void Mult(const Vector &x, Vector &y) const;
+   void ArrayMult(const Array<const Vector *> &X, Array<Vector *> &Y) const;
 
    // Set the operator.
-   void SetOperator( const Operator & op );
+   void SetOperator(const Operator &op);
 
    // Set various solver options. Refer to STRUMPACK documentation for
    // details.
-   void SetFromCommandLine( );
-   void SetPrintFactorStatistics( bool print_stat );
-   void SetPrintSolveStatistics( bool print_stat );
-   void SetRelTol( double rtol );
-   void SetAbsTol( double atol );
+   void SetFromCommandLine();
+   void SetPrintFactorStatistics(bool print_stat);
+   void SetPrintSolveStatistics(bool print_stat);
+
+   // Set tolerances and iterations for iterative solvers. Compression
+   // tolerance is handled below.
+   void SetRelTol(double rtol);
+   void SetAbsTol(double atol);
+   void SetMaxIter(int max_it);
+
+   // Set the flag controlling reuse of the symbolic factorization for multiple
+   // operators. This method has to be called before repeated calls to
+   // SetOperator.
+   void SetReorderingReuse(bool reuse);
+
+   // Enable or not GPU off-loading available if STRUMPACK was compiled with CUDA. Note
+   // that input/output from MFEM to STRUMPACK is all still through host memory.
+   void EnableGPU();
+   void DisableGPU();
 
    /**
     * STRUMPACK is an (approximate) direct solver. It can be used as a direct
@@ -100,70 +119,153 @@ public:
     * used without preconditioner.
     *
     * Supported values are:
-    *    AUTO:           Use iterative refinement if no HSS compression is used,
-    *                    otherwise use GMRes.
-    *    DIRECT:         No outer iterative solver, just a single application of
-    *                    the multifrontal solver.
-    *    REFINE:         Iterative refinement.
-    *    PREC_GMRES:     Preconditioned GMRes.
-    *                    The preconditioner is the (approx) multifrontal solver.
-    *    GMRES:          UN-preconditioned GMRes. (for testing mainly)
-    *    PREC_BICGSTAB:  Preconditioned BiCGStab.
-    *                    The preconditioner is the (approx) multifrontal solver.
+    *    AUTO:           Use iterative refinement if no HSS compression is
+    *                    used, otherwise use GMRes
+    *    DIRECT:         No outer iterative solver, just a single application
+    *                    of the multifrontal solver
+    *    REFINE:         Iterative refinement
+    *    PREC_GMRES:     Preconditioned GMRes
+    *                    The preconditioner is the (approx) multifrontal solver
+    *    GMRES:          UN-preconditioned GMRes (for testing mainly)
+    *    PREC_BICGSTAB:  Preconditioned BiCGStab
+    *                    The preconditioner is the (approx) multifrontal solver
     *    BICGSTAB:       UN-preconditioned BiCGStab. (for testing mainly)
     */
-   void SetKrylovSolver( strumpack::KrylovSolver method );
+   void SetKrylovSolver(strumpack::KrylovSolver method);
 
    /**
     * Supported reorderings are:
-    *    METIS, PARMETIS, SCOTCH, PTSCOTCH, RCM
+    *    NATURAL:    Do not reorder the system
+    *    METIS:      Use Metis nested-dissection reordering (default)
+    *    PARMETIS:   Use ParMetis nested-dissection reordering
+    *    SCOTCH:     Use Scotch nested-dissection reordering
+    *    PTSCOTCH:   Use PT-Scotch nested-dissection reordering
+    *    RCM:        Use RCM reordering
+    *    GEOMETRIC:  A simple geometric nested dissection code that
+    *                only works for regular meshes
+    *    AMD:        Approximate minimum degree
+    *    MMD:        Multiple minimum degree
+    *    AND:        Nested dissection
+    *    MLF:        Minimum local fill
+    *    SPECTRAL:   Spectral nested dissection
     */
-   void SetReorderingStrategy( strumpack::ReorderingStrategy method );
+   void SetReorderingStrategy(strumpack::ReorderingStrategy method);
 
    /**
-    * Disable static pivoting for stability. The static pivoting in strumpack
+    * Configure static pivoting for stability. The static pivoting in STRUMPACK
     * permutes the sparse input matrix in order to get large (nonzero) elements
     * on the diagonal. If the input matrix is already diagonally dominant, this
     * reordering can be disabled.
+    *
+    * Supported matching algorithms are:
+    *    NONE:                          Don't do anything
+    *    MAX_CARDINALITY:               Maximum cardinality
+    *    MAX_SMALLEST_DIAGONAL:         Maximum smallest diagonal value
+    *    MAX_SMALLEST_DIAGONAL_2:       Same as MAX_SMALLEST_DIAGONAL
+    *                                   but different algorithm
+    *    MAX_DIAGONAL_SUM:              Maximum sum of diagonal values
+    *    MAX_DIAGONAL_PRODUCT_SCALING:  Maximum product of diagonal values
+    *                                   and row and column scaling (default)
+    *    COMBBLAS:                      Use AWPM from CombBLAS (only with
+    *                                   version >= 3)
     */
-   void DisableMatching();
+   void SetMatching(strumpack::MatchingJob job);
 
    /**
-    * Enable static pivoting for stability using the MC64 algorithm with
-    * job=5. Using a matching algorithm, this will permute the sparse input
-    * matrix in order to get nonzero elements (as large as possible) on the
-    * diagonal. And will also scale the rows and columns of the matrix.
+    * Enable support for rank-structured data formats, which can be used
+    * for compression within the sparse solver.
+    *
+    * Supported compression types are:
+    *    NONE:           No compression, purely direct solver (default)
+    *    HSS:            HSS compression of frontal matrices
+    *    BLR:            Block low-rank compression of fronts
+    *    HODLR:          Hierarchically Off-diagonal Low-Rank
+    *                    compression of frontal matrices
+    *    BLR_HODLR:      Block low-rank compression of medium
+    *                    fronts and Hierarchically Off-diagonal
+    *                    Low-Rank compression of large fronts
+    *    ZFP_BLR_HODLR:  ZFP compression for small fronts,
+    *                    Block low-rank compression of medium
+    *                    fronts and Hierarchically Off-diagonal
+    *                    Low-Rank compression of large fronts
+    *    LOSSLESS:       Lossless compression
+    *    LOSSY:          Lossy compression
+    *
+    * For versions of STRUMPACK < 5, we support only NONE, HSS, and BLR.
+    * BLR_HODLR and ZPR_BLR_HODLR are supported in STRUMPACK >= 6.
     */
-   void EnableMatching();
-
-#if STRUMPACK_VERSION_MAJOR >= 3
-   /**
-    * Use the AWPM (approximate weight perfect matching) algorithm from the
-    * Combinatorial BLAS library for static pivoting, i.e. getting large
-    * nonzeros on the diagonal. This requires that strumpack was compiled with
-    * support for Combinatorial BLAS.
-    */
-   void EnableParallelMatching();
+   void SetCompression(strumpack::CompressionType type);
+   void SetCompressionRelTol(double rtol);
+   void SetCompressionAbsTol(double atol);
+#if STRUMPACK_VERSION_MAJOR >= 5
+   void SetCompressionLossyPrecision(int precision);
+   void SetCompressionButterflyLevels(int levels);
 #endif
 
 private:
-   void Init( int argc, char* argv[] );
+   // Helper method for calling the STRUMPACK factoriation routine.
+   void FactorInternal() const;
 
 protected:
-
-   MPI_Comm      comm_;
-   int           numProcs_;
-   int           myid_;
+   const STRUMPACKRowLocMatrix *APtr_;
+   STRUMPACKSolverType         *solver_;
 
    bool factor_verbose_;
    bool solve_verbose_;
+   bool reorder_reuse_;
 
-   const STRUMPACKRowLocMatrix * APtr_;
-   strumpack::StrumpackSparseSolverMPIDist<double,int> * solver_;
+   mutable Vector rhs_, sol_;
+   mutable int    nrhs_;
+};
 
-}; // mfem::STRUMPACKSolver class
+class STRUMPACKSolver :
+   public STRUMPACKSolverBase<strumpack::
+   SparseSolverMPIDist<double, HYPRE_BigInt>>
+{
+public:
+   // Constructor with MPI_Comm parameter.
+   STRUMPACKSolver(MPI_Comm comm);
 
-} // mfem namespace
+   // Constructor with STRUMPACK matrix object.
+   STRUMPACKSolver(STRUMPACKRowLocMatrix &A);
+
+   // Constructor with MPI_Comm parameter and command line arguments.
+   STRUMPACKSolver(MPI_Comm comm, int argc, char *argv[]);
+   MFEM_DEPRECATED STRUMPACKSolver(int argc, char *argv[], MPI_Comm comm)
+      : STRUMPACKSolver(comm, argc, argv) {}
+
+   // Constructor with STRUMPACK matrix object and command line arguments.
+   STRUMPACKSolver(STRUMPACKRowLocMatrix &A, int argc, char *argv[]);
+
+   // Destructor.
+   ~STRUMPACKSolver() {}
+};
+
+#if STRUMPACK_VERSION_MAJOR >= 7
+class STRUMPACKMixedPrecisionSolver :
+   public STRUMPACKSolverBase<strumpack::
+   SparseSolverMixedPrecisionMPIDist<float, double, HYPRE_BigInt>>
+{
+public:
+   // Constructor with MPI_Comm parameter.
+   STRUMPACKMixedPrecisionSolver(MPI_Comm comm);
+
+   // Constructor with STRUMPACK matrix object.
+   STRUMPACKMixedPrecisionSolver(STRUMPACKRowLocMatrix &A);
+
+   // Constructor with MPI_Comm parameter and command line arguments.
+   STRUMPACKMixedPrecisionSolver(MPI_Comm comm, int argc, char *argv[]);
+
+   // Constructor with STRUMPACK matrix object and command line arguments.
+   STRUMPACKMixedPrecisionSolver(STRUMPACKRowLocMatrix &A,
+                                 int argc, char *argv[]);
+
+   // Destructor.
+   ~STRUMPACKMixedPrecisionSolver() {}
+};
+#endif
+
+} // namespace mfem
 
 #endif // MFEM_USE_MPI
 #endif // MFEM_USE_STRUMPACK

linalg/superlu.cpp

@@ -650,6 +650,7 @@ void SuperLUSolver::ArrayMult(const Array<const Vector *> &X,
          MFEM_ASSERT(X[i], "Missing Vector in SuperLUSolver::Mult!");
          Vector s(sol_, i * ldx, ldx);
          s = *X[i];
+         sol_.SyncMemory(s);  // Update flags for sol_ if updated on device
       }
    }
 

linalg/vector.hpp

@@ -35,10 +35,6 @@
 namespace mfem
 {
 
-/** Count the number of entries in an array of doubles for which isfinite
-    is false, i.e. the entry is a NaN or +/-Inf. */
-inline int CheckFinite(const double *v, const int n);
-
 /// Define a shortcut for std::numeric_limits<double>::infinity()
 #ifndef __CYGWIN__
 inline double infinity()
@@ -444,7 +440,7 @@ public:
 
    /** @brief Count the number of entries in the Vector for which isfinite
        is false, i.e. the entry is a NaN or +/-Inf. */
-   int CheckFinite() const { return mfem::CheckFinite(HostRead(), size); }
+   int CheckFinite() const;
 
    /// Destroys vector.
    virtual ~Vector();
@@ -494,6 +490,8 @@ inline bool IsFinite(const double &val)
 #endif
 }
 
+/** Count the number of entries in an array of doubles for which isfinite
+    is false, i.e. the entry is a NaN or +/-Inf. */
 inline int CheckFinite(const double *v, const int n)
 {
    int bad = 0;
@@ -504,6 +502,11 @@ inline int CheckFinite(const double *v, const int n)
    return bad;
 }
 
+inline int Vector::CheckFinite() const
+{
+   return mfem::CheckFinite(HostRead(), size);
+}
+
 inline Vector::Vector(int s)
 {
    MFEM_ASSERT(s>=0,"Unexpected negative size.");

mesh/mesh.cpp

@@ -384,6 +384,12 @@ void Mesh::GetElementTransformation(int i, IsoparametricTransformation *ElTr)
    }
 }
 
+ElementTransformation *Mesh::GetElementTransformation(int i)
+{
+   GetElementTransformation(i, &Transformation);
+   return &Transformation;
+}
+
 void Mesh::GetElementTransformation(int i, const Vector &nodes,
                                     IsoparametricTransformation *ElTr)
 {
@@ -428,19 +434,6 @@ void Mesh::GetElementTransformation(int i, const Vector &nodes,
    }
 }
 
-ElementTransformation *Mesh::GetElementTransformation(int i)
-{
-   GetElementTransformation(i, &Transformation);
-
-   return &Transformation;
-}
-
-ElementTransformation *Mesh::GetBdrElementTransformation(int i)
-{
-   GetBdrElementTransformation(i, &BdrTransformation);
-   return &BdrTransformation;
-}
-
 void Mesh::GetBdrElementTransformation(int i, IsoparametricTransformation* ElTr)
 {
    ElTr->Attribute = GetBdrAttribute(i);
@@ -501,6 +494,12 @@ void Mesh::GetBdrElementTransformation(int i, IsoparametricTransformation* ElTr)
    }
 }
 
+ElementTransformation *Mesh::GetBdrElementTransformation(int i)
+{
+   GetBdrElementTransformation(i, &BdrTransformation);
+   return &BdrTransformation;
+}
+
 void Mesh::GetFaceTransformation(int FaceNo, IsoparametricTransformation *FTr)
 {
    FTr->Attribute = (Dim == 1) ? 1 : faces[FaceNo]->GetAttribute();
@@ -1102,7 +1101,7 @@ void Mesh::ApplyLocalSlaveTransformation(FaceElementTransformations &FT,
 FaceElementTransformations *Mesh::GetBdrFaceTransformations(int BdrElemNo)
 {
    FaceElementTransformations *tr;
-   int fn = GetBdrFace(BdrElemNo);
+   int fn = GetBdrElementFaceIndex(BdrElemNo);
 
    // Check if the face is interior, shared, or nonconforming.
    if (FaceIsTrueInterior(fn) || faces_info[fn].NCFace >= 0)
@@ -1117,24 +1116,6 @@ FaceElementTransformations *Mesh::GetBdrFaceTransformations(int BdrElemNo)
    return tr;
 }
 
-int Mesh::GetBdrFace(int BdrElemNo) const
-{
-   int fn;
-   if (Dim == 3)
-   {
-      fn = be_to_face[BdrElemNo];
-   }
-   else if (Dim == 2)
-   {
-      fn = be_to_edge[BdrElemNo];
-   }
-   else
-   {
-      fn = boundary[BdrElemNo]->GetVertices()[0];
-   }
-   return fn;
-}
-
 Mesh::FaceInformation Mesh::GetFaceInformation(int f) const
 {
    FaceInformation face;
@@ -1448,7 +1429,7 @@ Array<int> Mesh::GetFaceToBdrElMap() const
    face_to_be = -1;
    for (int i = 0; i < NumOfBdrElements; i++)
    {
-      face_to_be[GetBdrElementEdgeIndex(i)] = i;
+      face_to_be[GetBdrElementFaceIndex(i)] = i;
    }
    return face_to_be;
 }
@@ -1539,7 +1520,6 @@ void Mesh::Destroy()
    faces.DeleteAll();
    faces_info.DeleteAll();
    nc_faces_info.DeleteAll();
-   be_to_edge.DeleteAll();
    be_to_face.DeleteAll();
 
    // TODO:
@@ -1906,12 +1886,7 @@ int Mesh::AddBdrPoint(int v, int attr)
 
 void Mesh::GenerateBoundaryElements()
 {
-   int i, j;
-   Array<int> &be2face = (Dim == 2) ? be_to_edge : be_to_face;
-
-   // GenerateFaces();
-
-   for (i = 0; i < boundary.Size(); i++)
+   for (int i = 0; i < boundary.Size(); i++)
    {
       FreeElement(boundary[i]);
    }
@@ -1924,22 +1899,24 @@ void Mesh::GenerateBoundaryElements()
 
    // count the 'NumOfBdrElements'
    NumOfBdrElements = 0;
-   for (i = 0; i < faces_info.Size(); i++)
+   for (int i = 0; i < faces_info.Size(); i++)
    {
       if (faces_info[i].Elem2No < 0) { NumOfBdrElements++; }
    }
 
+   // Add the boundary elements
    boundary.SetSize(NumOfBdrElements);
-   be2face.SetSize(NumOfBdrElements);
-   for (j = i = 0; i < faces_info.Size(); i++)
+   be_to_face.SetSize(NumOfBdrElements);
+   for (int i = 0, j = 0; i < faces_info.Size(); i++)
    {
       if (faces_info[i].Elem2No < 0)
       {
          boundary[j] = faces[i]->Duplicate(this);
-         be2face[j++] = i;
+         be_to_face[j++] = i;
       }
    }
-   // In 3D, 'bel_to_edge' is destroyed but it's not updated.
+
+   // Note: in 3D, 'bel_to_edge' is destroyed but it's not updated.
 }
 
 void Mesh::FinalizeCheck()
@@ -1969,7 +1946,7 @@ void Mesh::FinalizeTriMesh(int generate_edges, int refine, bool fix_orientation)
    if (generate_edges)
    {
       el_to_edge = new Table;
-      NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
+      NumOfEdges = GetElementToEdgeTable(*el_to_edge);
       GenerateFaces();
       CheckBdrElementOrientation();
    }
@@ -1997,7 +1974,7 @@ void Mesh::FinalizeQuadMesh(int generate_edges, int refine,
    if (generate_edges)
    {
       el_to_edge = new Table;
-      NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
+      NumOfEdges = GetElementToEdgeTable(*el_to_edge);
       GenerateFaces();
       CheckBdrElementOrientation();
    }
@@ -2299,8 +2276,7 @@ void Mesh::ReorderElements(const Array<int> &ordering, bool reorder_vertices)
    // - edge_vertex - no need to rebuild
    // - geom_factors - no need to rebuild
 
-   // - be_to_edge  - 2D only
-   // - be_to_face  - 3D only
+   // - be_to_face
 
    // - Nodes
 
@@ -2386,9 +2362,9 @@ void Mesh::ReorderElements(const Array<int> &ordering, bool reorder_vertices)
 
    if (Dim > 1)
    {
-      // generate el_to_edge, be_to_edge (2D), bel_to_edge (3D)
+      // generate el_to_edge, be_to_face (2D), bel_to_edge (3D)
       el_to_edge = new Table;
-      NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
+      NumOfEdges = GetElementToEdgeTable(*el_to_edge);
    }
    if (Dim > 2)
    {
@@ -2778,8 +2754,8 @@ void Mesh::DoNodeReorder(DSTable *old_v_to_v, Table *old_elem_vert)
    }
    if (el_to_edge)
    {
-      // update 'el_to_edge', 'be_to_edge' (2D), 'bel_to_edge' (3D)
-      NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
+      // update 'el_to_edge', 'be_to_face' (2D), 'bel_to_edge' (3D)
+      NumOfEdges = GetElementToEdgeTable(*el_to_edge);
       if (Dim == 2)
       {
          // update 'faces' and 'faces_info'
@@ -2856,7 +2832,7 @@ void Mesh::FinalizeTetMesh(int generate_edges, int refine, bool fix_orientation)
    if (generate_edges == 1)
    {
       el_to_edge = new Table;
-      NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
+      NumOfEdges = GetElementToEdgeTable(*el_to_edge);
    }
    else
    {
@@ -2891,7 +2867,7 @@ void Mesh::FinalizeWedgeMesh(int generate_edges, int refine,
    if (generate_edges == 1)
    {
       el_to_edge = new Table;
-      NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
+      NumOfEdges = GetElementToEdgeTable(*el_to_edge);
    }
    else
    {
@@ -2923,7 +2899,7 @@ void Mesh::FinalizeHexMesh(int generate_edges, int refine, bool fix_orientation)
    if (generate_edges)
    {
       el_to_edge = new Table;
-      NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
+      NumOfEdges = GetElementToEdgeTable(*el_to_edge);
    }
    else
    {
@@ -3002,7 +2978,7 @@ void Mesh::FinalizeTopology(bool generate_bdr)
    {
       // el_to_edge may already be allocated (P2 VTK meshes)
       if (!el_to_edge) { el_to_edge = new Table; }
-      NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
+      NumOfEdges = GetElementToEdgeTable(*el_to_edge);
       if (Dim == 2)
       {
          GenerateFaces(); // 'Faces' in 2D refers to the edges
@@ -3022,8 +2998,17 @@ void Mesh::FinalizeTopology(bool generate_bdr)
       GenerateFaces();
       if (NumOfBdrElements == 0 && generate_bdr)
       {
+         // be_to_face will be set inside GenerateBoundaryElements
          GenerateBoundaryElements();
       }
+      else
+      {
+         be_to_face.SetSize(NumOfBdrElements);
+         for (int i = 0; i < NumOfBdrElements; ++i)
+         {
+            be_to_face[i] = boundary[i]->GetVertices()[0];
+         }
+      }
    }
 
    if (ncmesh)
@@ -3555,7 +3540,7 @@ void Mesh::Make2D(int nx, int ny, Element::Type type,
    if (generate_edges == 1)
    {
       el_to_edge = new Table;
-      NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
+      NumOfEdges = GetElementToEdgeTable(*el_to_edge);
       GenerateFaces();
       CheckBdrElementOrientation();
    }
@@ -3613,6 +3598,11 @@ void Mesh::Make1D(int n, double sx)
    SetMeshGen();
    GenerateFaces();
 
+   // Set be_to_face
+   be_to_face.SetSize(2);
+   be_to_face[0] = 0;
+   be_to_face[1] = n;
+
    attributes.Append(1);
    bdr_attributes.Append(1); bdr_attributes.Append(2);
 }
@@ -3666,9 +3656,6 @@ Mesh::Mesh(const Mesh &mesh, bool copy_nodes)
    // Copy the boundary-to-edge Table, bel_to_edge (3D)
    bel_to_edge = (mesh.bel_to_edge) ? new Table(*mesh.bel_to_edge) : NULL;
 
-   // Copy the boundary-to-edge Array, be_to_edge (2D)
-   mesh.be_to_edge.Copy(be_to_edge);
-
    // Duplicate the faces and faces_info.
    faces.SetSize(mesh.faces.Size());
    for (int i = 0; i < faces.Size(); i++)
@@ -5269,7 +5256,7 @@ void Mesh::UpdateNURBS()
 
    if (el_to_edge)
    {
-      NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
+      NumOfEdges = GetElementToEdgeTable(*el_to_edge);
    }
 
    if (el_to_face)
@@ -6169,15 +6156,15 @@ int Mesh::CheckBdrElementOrientation(bool fix_it)
       if (el_to_edge == NULL) // edges were not generated
       {
          el_to_edge = new Table;
-         NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
+         NumOfEdges = GetElementToEdgeTable(*el_to_edge);
          GenerateFaces(); // 'Faces' in 2D refers to the edges
       }
       for (int i = 0; i < NumOfBdrElements; i++)
       {
-         if (faces_info[be_to_edge[i]].Elem2No < 0) // boundary face
+         if (faces_info[be_to_face[i]].Elem2No < 0) // boundary face
          {
             int *bv = boundary[i]->GetVertices();
-            int *fv = faces[be_to_edge[i]]->GetVertices();
+            int *fv = faces[be_to_face[i]]->GetVertices();
             if (bv[0] != fv[0])
             {
                if (fix_it)
@@ -6428,7 +6415,7 @@ void Mesh::GetBdrElementEdges(int i, Array<int> &edges, Array<int> &cor) const
    {
       edges.SetSize(1);
       cor.SetSize(1);
-      edges[0] = be_to_edge[i];
+      edges[0] = be_to_face[i];
       const int *v = boundary[i]->GetVertices();
       cor[0] = (v[0] < v[1]) ? (1) : (-1);
    }
@@ -6669,7 +6656,7 @@ Array<int> Mesh::FindFaceNeighbors(const int elem) const
 
 void Mesh::GetBdrElementFace(int i, int *f, int *o) const
 {
-   *f = GetBdrElementEdgeIndex(i);
+   *f = GetBdrElementFaceIndex(i);
 
    const int *fv = (Dim > 1) ? faces[*f]->GetVertices() : NULL;
    const int *bv = boundary[i]->GetVertices();
@@ -6686,21 +6673,9 @@ void Mesh::GetBdrElementFace(int i, int *f, int *o) const
    }
 }
 
-int Mesh::GetBdrElementEdgeIndex(int i) const
-{
-   switch (Dim)
-   {
-      case 1: return boundary[i]->GetVertices()[0];
-      case 2: return be_to_edge[i];
-      case 3: return be_to_face[i];
-      default: MFEM_ABORT("invalid dimension!");
-   }
-   return -1;
-}
-
 void Mesh::GetBdrElementAdjacentElement(int bdr_el, int &el, int &info) const
 {
-   int fid = GetBdrElementEdgeIndex(bdr_el);
+   int fid = GetBdrElementFaceIndex(bdr_el);
 
    const FaceInfo &fi = faces_info[fid];
    MFEM_ASSERT(fi.Elem1Inf % 64 == 0, "internal error"); // orientation == 0
@@ -6722,7 +6697,7 @@ void Mesh::GetBdrElementAdjacentElement(int bdr_el, int &el, int &info) const
 
 void Mesh::GetBdrElementAdjacentElement2(int bdr_el, int &el, int &info) const
 {
-   int fid = GetBdrElementEdgeIndex(bdr_el);
+   int fid = GetBdrElementFaceIndex(bdr_el);
 
    const FaceInfo &fi = faces_info[fid];
    MFEM_ASSERT(fi.Elem1Inf % 64 == 0, "internal error"); // orientation == 0
@@ -6848,7 +6823,7 @@ void Mesh::GetVertexToVertexTable(DSTable &v_to_v) const
    }
 }
 
-int Mesh::GetElementToEdgeTable(Table & e_to_f, Array<int> &be_to_f)
+int Mesh::GetElementToEdgeTable(Table &e_to_f)
 {
    int i, NumberOfEdges;
 
@@ -6863,11 +6838,11 @@ int Mesh::GetElementToEdgeTable(Table & e_to_f, Array<int> &be_to_f)
    if (Dim == 2)
    {
       // Initialize the indices for the boundary elements.
-      be_to_f.SetSize(NumOfBdrElements);
+      be_to_face.SetSize(NumOfBdrElements);
       for (i = 0; i < NumOfBdrElements; i++)
       {
          const int *v = boundary[i]->GetVertices();
-         be_to_f[i] = v_to_v(v[0], v[1]);
+         be_to_face[i] = v_to_v(v[0], v[1]);
       }
    }
    else if (Dim == 3)
@@ -7470,7 +7445,7 @@ void Mesh::ReorientTetMesh()
       GenerateFaces();
       if (el_to_edge)
       {
-         NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
+         NumOfEdges = GetElementToEdgeTable(*el_to_edge);
       }
    }
    else
@@ -8505,7 +8480,7 @@ void Mesh::UniformRefinement2D_base(bool update_nodes)
    if (el_to_edge == NULL)
    {
       el_to_edge = new Table;
-      NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
+      NumOfEdges = GetElementToEdgeTable(*el_to_edge);
    }
 
    int quad_counter = 0;
@@ -8594,8 +8569,8 @@ void Mesh::UniformRefinement2D_base(bool update_nodes)
       const int attr = boundary[i]->GetAttribute();
       int *v = boundary[i]->GetVertices();
 
-      new_boundary[j++] = new Segment(v[0], oedge+be_to_edge[i], attr);
-      new_boundary[j++] = new Segment(oedge+be_to_edge[i], v[1], attr);
+      new_boundary[j++] = new Segment(v[0], oedge+be_to_face[i], attr);
+      new_boundary[j++] = new Segment(oedge+be_to_face[i], v[1], attr);
 
       FreeElement(boundary[i]);
    }
@@ -8635,7 +8610,7 @@ void Mesh::UniformRefinement2D_base(bool update_nodes)
    NumOfBdrElements = 2 * NumOfBdrElements;
    NumOfFaces       = 0;
 
-   NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
+   NumOfEdges = GetElementToEdgeTable(*el_to_edge);
    GenerateFaces();
 
    last_operation = Mesh::REFINE;
@@ -8665,7 +8640,7 @@ void Mesh::UniformRefinement3D_base(Array<int> *f2qf_ptr, DSTable *v_to_v_p,
    if (el_to_edge == NULL)
    {
       el_to_edge = new Table;
-      NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
+      NumOfEdges = GetElementToEdgeTable(*el_to_edge);
    }
 
    if (el_to_face == NULL)
@@ -9356,7 +9331,7 @@ void Mesh::UniformRefinement3D_base(Array<int> *f2qf_ptr, DSTable *v_to_v_p,
    CheckBdrElementOrientation(false);
 #endif
 
-   NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
+   NumOfEdges = GetElementToEdgeTable(*el_to_edge);
 
    last_operation = Mesh::REFINE;
    sequence++;
@@ -9490,7 +9465,7 @@ void Mesh::LocalRefinement(const Array<int> &marked_el, int type)
 
       if (el_to_edge != NULL)
       {
-         NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
+         NumOfEdges = GetElementToEdgeTable(*el_to_edge);
          GenerateFaces();
       }
 
@@ -9586,7 +9561,7 @@ void Mesh::LocalRefinement(const Array<int> &marked_el, int type)
       // 5. Update element-to-edge and element-to-face relations.
       if (el_to_edge != NULL)
       {
-         NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
+         NumOfEdges = GetElementToEdgeTable(*el_to_edge);
       }
       if (el_to_face != NULL)
       {
@@ -9777,7 +9752,7 @@ void Mesh::InitFromNCMesh(const NCMesh &ncmesh_)
    if (Dim > 1)
    {
       el_to_edge = new Table;
-      NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
+      NumOfEdges = GetElementToEdgeTable(*el_to_edge);
    }
    if (Dim > 2)
    {
@@ -9824,7 +9799,6 @@ void Mesh::Swap(Mesh& other, bool non_geometry)
    mfem::Swap(el_to_edge, other.el_to_edge);
    mfem::Swap(el_to_face, other.el_to_face);
    mfem::Swap(el_to_el, other.el_to_el);
-   mfem::Swap(be_to_edge, other.be_to_edge);
    mfem::Swap(bel_to_edge, other.bel_to_edge);
    mfem::Swap(be_to_face, other.be_to_face);
    mfem::Swap(face_edge, other.face_edge);
@@ -12279,9 +12253,9 @@ void Mesh::RemoveUnusedVertices()
    DeleteTables();
    if (Dim > 1)
    {
-      // generate el_to_edge, be_to_edge (2D), bel_to_edge (3D)
+      // generate el_to_edge, be_to_face (2D), bel_to_edge (3D)
       el_to_edge = new Table;
-      NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
+      NumOfEdges = GetElementToEdgeTable(*el_to_edge);
    }
    if (Dim > 2)
    {
@@ -12312,7 +12286,7 @@ void Mesh::RemoveInternalBoundaries()
    int new_bel_to_edge_nnz = 0;
    for (int i = 0; i < GetNBE(); i++)
    {
-      if (FaceIsInterior(GetBdrElementEdgeIndex(i)))
+      if (FaceIsInterior(GetBdrElementFaceIndex(i)))
       {
          FreeElement(boundary[i]);
       }
@@ -12329,32 +12303,24 @@ void Mesh::RemoveInternalBoundaries()
    if (num_bdr_elem == GetNBE()) { return; }
 
    Array<Element *> new_boundary(num_bdr_elem);
-   Array<int> new_be_to_edge, new_be_to_face;
+   Array<int> new_be_to_face;
    Table *new_bel_to_edge = NULL;
    new_boundary.SetSize(0);
-   if (Dim == 2)
+   new_be_to_face.Reserve(num_bdr_elem);
+   if (Dim == 3)
    {
-      new_be_to_edge.Reserve(num_bdr_elem);
-   }
-   else if (Dim == 3)
-   {
-      new_be_to_face.Reserve(num_bdr_elem);
       new_bel_to_edge = new Table;
       new_bel_to_edge->SetDims(num_bdr_elem, new_bel_to_edge_nnz);
    }
    for (int i = 0; i < GetNBE(); i++)
    {
-      if (!FaceIsInterior(GetBdrElementEdgeIndex(i)))
+      if (!FaceIsInterior(GetBdrElementFaceIndex(i)))
       {
          new_boundary.Append(boundary[i]);
-         if (Dim == 2)
+         int row = new_be_to_face.Size();
+         new_be_to_face.Append(be_to_face[i]);
+         if (Dim == 3)
          {
-            new_be_to_edge.Append(be_to_edge[i]);
-         }
-         else if (Dim == 3)
-         {
-            int row = new_be_to_face.Size();
-            new_be_to_face.Append(be_to_face[i]);
             int *e = bel_to_edge->GetRow(i);
             int ne = bel_to_edge->RowSize(i);
             int *new_e = new_bel_to_edge->GetRow(row);
@@ -12370,13 +12336,10 @@ void Mesh::RemoveInternalBoundaries()
    NumOfBdrElements = new_boundary.Size();
    mfem::Swap(boundary, new_boundary);
 
-   if (Dim == 2)
+   mfem::Swap(be_to_face, new_be_to_face);
+
+   if (Dim == 3)
    {
-      mfem::Swap(be_to_edge, new_be_to_edge);
-   }
-   else if (Dim == 3)
-   {
-      mfem::Swap(be_to_face, new_be_to_face);
       delete bel_to_edge;
       bel_to_edge = new_bel_to_edge;
    }

mesh/mesh.hpp

@@ -220,9 +220,9 @@ protected:
    Table *el_to_edge;
    Table *el_to_face;
    Table *el_to_el;
-   Array<int> be_to_edge;  // for 2D
-   Table *bel_to_edge;     // for 3D
-   Array<int> be_to_face;
+   Array<int> be_to_face; // faces = vertices (1D), edges (2D), faces (3D)
+
+   Table *bel_to_edge;    // for 3D only
 
    // Note that the following tables are owned by this class and should not be
    // deleted by the caller. Of these three tables, only face_edge and
@@ -533,7 +533,7 @@ protected:
        nodes in the elements. For example, if T is the element to edge table
        T(i, 0) gives the index of edge in element i that connects vertex 0
        to vertex 1, etc. Returns the number of the edges. */
-   int GetElementToEdgeTable(Table &, Array<int> &);
+   int GetElementToEdgeTable(Table &);
 
    /// Used in GenerateFaces()
    void AddPointFaceElement(int lf, int gf, int el);
@@ -1367,11 +1367,6 @@ public:
        GetElementEdges/GetBdrElementEdges. */
    void GetBdrElementFace(int i, int *f, int *o) const;
 
-   /** Return the vertex index of boundary element i. (1D)
-       Return the edge index of boundary element i. (2D)
-       Return the face index of boundary element i. (3D) */
-   int GetBdrElementEdgeIndex(int i) const;
-
    /** @brief For the given boundary element, bdr_el, return its adjacent
        element and its info, i.e. 64*local_bdr_index+bdr_orientation.
 
@@ -1392,8 +1387,19 @@ public:
        @sa GetBdrElementAdjacentElement() */
    void GetBdrElementAdjacentElement2(int bdr_el, int &el, int &info) const;
 
-   /// Return the local face index for the given boundary face.
-   int GetBdrFace(int BdrElemNo) const;
+   /// @brief Return the local face (codimension-1) index for the given boundary
+   /// element index.
+   int GetBdrElementFaceIndex(int be_idx) const { return be_to_face[be_idx]; }
+
+   /// Deprecated in favor of GetBdrElementFaceIndex().
+   MFEM_DEPRECATED int GetBdrFace(int i) const { return GetBdrElementFaceIndex(i); }
+
+   /** Return the vertex index of boundary element i. (1D)
+       Return the edge index of boundary element i. (2D)
+       Return the face index of boundary element i. (3D)
+
+       Deprecated in favor of GetBdrElementFaceIndex(). */
+   MFEM_DEPRECATED int GetBdrElementEdgeIndex(int i) const { return GetBdrElementFaceIndex(i); }
 
    /// @}
 

mesh/pmesh.cpp

@@ -186,7 +186,7 @@ ParMesh::ParMesh(MPI_Comm comm, Mesh &mesh, int *partitioning_,
       if (Dim > 1)
       {
          el_to_edge = new Table;
-         NumOfEdges = Mesh::GetElementToEdgeTable(*el_to_edge, be_to_edge);
+         NumOfEdges = Mesh::GetElementToEdgeTable(*el_to_edge);
       }
 
       STable3D *faces_tbl = NULL;
@@ -197,6 +197,19 @@ ParMesh::ParMesh(MPI_Comm comm, Mesh &mesh, int *partitioning_,
 
       GenerateFaces();
 
+      // Make sure the be_to_face array is initialized.
+      // In 2D, it will be set in the above call to Mesh::GetElementToEdgeTable.
+      // In 3D, it will be set in GetElementToFaceTable.
+      // In 1D, we need to set it manually.
+      if (Dim == 1)
+      {
+         be_to_face.SetSize(NumOfBdrElements);
+         for (int i = 0; i < NumOfBdrElements; ++i)
+         {
+            be_to_face[i] = boundary[i]->GetVertices()[0];
+         }
+      }
+
       ListOfIntegerSets  groups;
       {
          // the first group is the local one
@@ -435,7 +448,7 @@ int ParMesh::BuildLocalBoundary(const Mesh& mesh, const int* partitioning,
 
       for (int i = 0; i < mesh.GetNBE(); i++)
       {
-         int edge = mesh.GetBdrElementEdgeIndex(i);
+         int edge = mesh.GetBdrElementFaceIndex(i);
          int el1 = edge_element->GetRow(edge)[0];
          if (partitioning[el1] == MyRank)
          {
@@ -451,7 +464,7 @@ int ParMesh::BuildLocalBoundary(const Mesh& mesh, const int* partitioning,
       boundary.SetSize(nbdry);
       for (int i = 0; i < mesh.GetNBE(); i++)
       {
-         int edge = mesh.GetBdrElementEdgeIndex(i);
+         int edge = mesh.GetBdrElementFaceIndex(i);
          int el1 = edge_element->GetRow(edge)[0];
          if (partitioning[el1] == MyRank)
          {
@@ -3320,7 +3333,7 @@ void ParMesh::ReorientTetMesh()
       GenerateFaces();
       if (el_to_edge)
       {
-         NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
+         NumOfEdges = GetElementToEdgeTable(*el_to_edge);
       }
    }
    else
@@ -3555,7 +3568,7 @@ void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
       // 6. Update element-to-edge relations.
       if (el_to_edge != NULL)
       {
-         NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
+         NumOfEdges = GetElementToEdgeTable(*el_to_edge);
       }
    } //  'if (Dim == 3)'
 
@@ -3793,7 +3806,7 @@ void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
 
       if (el_to_edge != NULL)
       {
-         NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
+         NumOfEdges = GetElementToEdgeTable(*el_to_edge);
          GenerateFaces();
       }
    } //  'if (Dim == 2)'
@@ -3855,6 +3868,8 @@ void ParMesh::NonconformingRefinement(const Array<Refinement> &refinements,
                  "serial Mesh)");
    }
 
+   ResetLazyData();
+
    DeleteFaceNbrData();
 
    // NOTE: no check of !refinements.Size(), in parallel we would have to reduce
@@ -5286,7 +5301,7 @@ Mesh ParMesh::GetSerialMesh(int save_rank) const
       for (int e = 0; e < NumOfElements; e++)
       {
          const int attr = elements[e]->GetAttribute();
-         const int geom_type = elements[e]->GetGeometryType();;
+         const int geom_type = elements[e]->GetGeometryType();
          ints.Append(attr);
          ints.Append(geom_type);
          pfespace_linear.GetElementDofs(e, dofs);

mesh/pumi.cpp

@@ -1028,9 +1028,6 @@ void ParPumiMesh::UpdateMesh(const ParMesh* AdaptedpMesh)
    bel_to_edge = (AdaptedpMesh->bel_to_edge) ?
                  new Table(*(AdaptedpMesh->bel_to_edge)) : NULL;
 
-   // Copy the boudary-to-edge Array, be_to_edge (2D)
-   AdaptedpMesh->be_to_edge.Copy(be_to_edge);
-
    // Duplicate the faces and faces_info.
    faces.SetSize(AdaptedpMesh->faces.Size());
    for (int i = 0; i < faces.Size(); i++)

mesh/submesh/psubmesh.cpp

@@ -235,10 +235,13 @@ ParSubMesh::ParSubMesh(const ParMesh &parent, SubMesh::From from,
 
    // Add boundaries
    {
-      int num_of_faces_or_edges =
-         (Dim == 3) ? NumOfFaces :
-         ((Dim == 2) ? NumOfEdges : NumOfVertices);
-      Array<int> &be2face = (Dim == 2) ? be_to_edge : be_to_face;
+      const int num_codim_1 = [this]()
+      {
+         if (Dim == 1) { return NumOfVertices; }
+         else if (Dim == 2) { return NumOfEdges; }
+         else if (Dim == 3) { return NumOfFaces; }
+         else { MFEM_ABORT("Invalid dimension."); return -1; }
+      }();
 
       if (Dim == 3)
       {
@@ -249,7 +252,7 @@ ParSubMesh::ParSubMesh(const ParMesh &parent, SubMesh::From from,
       }
 
       NumOfBdrElements = 0;
-      for (int i = 0; i < num_of_faces_or_edges; i++)
+      for (int i = 0; i < num_codim_1; i++)
       {
          if (GetFaceInformation(i).IsBoundary())
          {
@@ -258,14 +261,17 @@ ParSubMesh::ParSubMesh(const ParMesh &parent, SubMesh::From from,
       }
 
       boundary.SetSize(NumOfBdrElements);
-      be2face.SetSize(NumOfBdrElements);
+      be_to_face.SetSize(NumOfBdrElements);
       Array<int> parent_face_to_be = parent.GetFaceToBdrElMap();
       int max_bdr_attr = parent.bdr_attributes.Max();
-      for (int i = 0, j = 0; i < num_of_faces_or_edges; i++)
+
+      for (int i = 0, j = 0; i < num_codim_1; i++)
       {
          if (GetFaceInformation(i).IsBoundary())
          {
             boundary[j] = faces[i]->Duplicate(this);
+            be_to_face[j] = i;
+
             if (from == SubMesh::From::Domain && Dim >= 2)
             {
                int pbeid = Dim == 3 ? parent_face_to_be[parent_face_ids_[i]] :
@@ -283,7 +289,7 @@ ParSubMesh::ParSubMesh(const ParMesh &parent, SubMesh::From from,
             {
                boundary[j]->SetAttribute(SubMesh::GENERATED_ATTRIBUTE);
             }
-            be2face[j++] = i;
+            ++j;
          }
       }
    }
@@ -323,7 +329,7 @@ ParSubMesh::ParSubMesh(const ParMesh &parent, SubMesh::From from,
    if (Dim > 1)
    {
       if (!el_to_edge) { el_to_edge = new Table; }
-      NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
+      NumOfEdges = GetElementToEdgeTable(*el_to_edge);
    }
 
    SetAttributes();

mesh/submesh/ptransfermap.cpp

@@ -317,8 +317,7 @@ ParTransferMap::CorrectFaceOrientations(const ParFiniteElementSpace &fes,
 
    if (parent_face_ori.Size() == 0) { return; }
 
-   VDofTransformation vdoftrans(fes.GetVDim(),
-                                fes.GetOrdering());
+   DofTransformation doftrans(fes.GetVDim(), fes.GetOrdering());
 
    int dim = mesh->Dimension();
    bool face = (dim == 3);
@@ -332,17 +331,13 @@ ParTransferMap::CorrectFaceOrientations(const ParFiniteElementSpace &fes,
       if (parent_face_ori[i] == 0) { continue; }
 
       Geometry::Type geom = face ? mesh->GetFaceGeometry(i) :
-                            mesh->GetElementGeometry(i);;
+                            mesh->GetElementGeometry(i);
 
-      StatelessDofTransformation * doftrans =
-         fec->DofTransformationForGeometry(geom);
-
-      if (doftrans == NULL) { continue; }
-
-      vdoftrans.SetDofTransformation(*doftrans);
+      if (!fec->DofTransformationForGeometry(geom)) { continue; }
+      doftrans.SetDofTransformation(*fec->DofTransformationForGeometry(geom));
 
       Fo[0] = parent_face_ori[i];
-      vdoftrans.SetFaceOrientations(Fo);
+      doftrans.SetFaceOrientations(Fo);
 
       if (face)
       {
@@ -356,12 +351,12 @@ ParTransferMap::CorrectFaceOrientations(const ParFiniteElementSpace &fes,
       if (sub_to_parent_map)
       {
          src.GetSubVector(vdofs, face_vector);
-         vdoftrans.TransformPrimal(face_vector);
+         doftrans.TransformPrimal(face_vector);
       }
       else
       {
          dst.GetSubVector(vdofs, face_vector);
-         vdoftrans.InvTransformPrimal(face_vector);
+         doftrans.InvTransformPrimal(face_vector);
       }
 
       for (int j = 0; j < vdofs.Size(); j++)

mesh/submesh/submesh.cpp

@@ -65,7 +65,7 @@ SubMesh::SubMesh(const Mesh &parent, From from,
 
       for (int i = 0; i < NumOfBdrElements; i++)
       {
-         int pbeid = parent_face_to_be[parent_face_ids_[GetBdrFace(i)]];
+         int pbeid = parent_face_to_be[parent_face_ids_[GetBdrElementFaceIndex(i)]];
          if (pbeid != -1)
          {
             int attr = parent.GetBdrElement(pbeid)->GetAttribute();
@@ -117,7 +117,7 @@ SubMesh::SubMesh(const Mesh &parent, From from,
 
          for (int i = 0; i < NumOfBdrElements; i++)
          {
-            int pbeid = parent_face_to_be[parent_edge_ids_[GetBdrFace(i)]];
+            int pbeid = parent_face_to_be[parent_edge_ids_[GetBdrElementFaceIndex(i)]];
             if (pbeid != -1)
             {
                int attr = parent.GetBdrElement(pbeid)->GetAttribute();

mesh/submesh/transfermap.cpp

@@ -241,8 +241,7 @@ void TransferMap::CorrectFaceOrientations(const FiniteElementSpace &fes,
 
    if (parent_face_ori.Size() == 0) { return; }
 
-   VDofTransformation vdoftrans(fes.GetVDim(),
-                                fes.GetOrdering());
+   DofTransformation doftrans(fes.GetVDim(), fes.GetOrdering());
 
    int dim = mesh->Dimension();
    bool face = (dim == 3);
@@ -256,17 +255,13 @@ void TransferMap::CorrectFaceOrientations(const FiniteElementSpace &fes,
       if (parent_face_ori[i] == 0) { continue; }
 
       Geometry::Type geom = face ? mesh->GetFaceGeometry(i) :
-                            mesh->GetElementGeometry(i);;
+                            mesh->GetElementGeometry(i);
 
-      StatelessDofTransformation * doftrans =
-         fec->DofTransformationForGeometry(geom);
-
-      if (doftrans == NULL) { continue; }
-
-      vdoftrans.SetDofTransformation(*doftrans);
+      if (!fec->DofTransformationForGeometry(geom)) { continue; }
+      doftrans.SetDofTransformation(*fec->DofTransformationForGeometry(geom));
 
       Fo[0] = parent_face_ori[i];
-      vdoftrans.SetFaceOrientations(Fo);
+      doftrans.SetFaceOrientations(Fo);
 
       if (face)
       {
@@ -280,12 +275,12 @@ void TransferMap::CorrectFaceOrientations(const FiniteElementSpace &fes,
       if (sub_to_parent_map)
       {
          src.GetSubVector(vdofs, face_vector);
-         vdoftrans.TransformPrimal(face_vector);
+         doftrans.TransformPrimal(face_vector);
       }
       else
       {
          dst.GetSubVector(vdofs, face_vector);
-         vdoftrans.InvTransformPrimal(face_vector);
+         doftrans.InvTransformPrimal(face_vector);
       }
 
       for (int j = 0; j < vdofs.Size(); j++)

miniapps/dpg/pmaxwell.cpp

@@ -1403,7 +1403,7 @@ void maxwell_solution_curlcurl(const Vector & X,
    curlcurlE.resize(dim);
    for (int i = 0; i < dim; ++i)
    {
-      curlcurlE[i] = 0.0;;
+      curlcurlE[i] = 0.0;
    }
    switch (prob)
    {

miniapps/hdiv-linear-solver/grad_div.cpp

@@ -278,4 +278,4 @@ void SolveCG(Operator &A, Solver &P, const Vector &B, Vector &X)
       cout << "Done.\nIterations: " << cg.GetNumIterations()
            << "\nElapsed: " << tic_toc.RealTime() << endl;
    }
-};
+}

miniapps/meshing/mesh-optimizer.cpp

@@ -35,7 +35,7 @@
 //   Adapted analytic shape:
 //     mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 2 -tid 4 -ni 200 -bnd -qt 1 -qo 8
 //   Adapted analytic size+orientation:
-//     mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 14 -tid 4 -ni 100 -bnd -qt 1 -qo 8 -fd
+//     mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 14 -tid 4 -ni 100 -bnd -qt 1 -qo 8
 //   Adapted analytic shape+orientation:
 //     mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 85 -tid 4 -ni 100 -bnd -qt 1 -qo 8 -fd
 //

miniapps/meshing/minimal-surface.cpp

@@ -1012,7 +1012,7 @@ struct QuarterPeach: public Surface
       for (int i = 0; i < GetNBE(); i++)
       {
          Element *el = GetBdrElement(i);
-         const int fn = GetBdrElementEdgeIndex(i);
+         const int fn = GetBdrElementFaceIndex(i);
          MFEM_VERIFY(!FaceIsTrueInterior(fn),"");
          Array<int> vertices;
          GetFaceVertices(fn, vertices);

miniapps/meshing/pmesh-fitting.cpp

@@ -359,7 +359,6 @@ int main (int argc, char *argv[])
    }
    pmesh->ExchangeFaceNbrData();
 
-
    // Surface fitting.
    L2_FECollection mat_coll(0, dim);
    H1_FECollection surf_fit_fec(mesh_poly_deg, dim);
@@ -432,7 +431,6 @@ int main (int argc, char *argv[])
          }
          else { surf_fit_bg_gf0->ProjectCoefficient(*ls_coeff); }
 
-
          surf_fit_bg_grad_fes =
             new ParFiniteElementSpace(pmesh_surf_fit_bg, surf_fit_bg_fec, dim);
          surf_fit_bg_grad = new ParGridFunction(surf_fit_bg_grad_fes);
@@ -630,6 +628,7 @@ int main (int argc, char *argv[])
          }
       }
    }
+   pmesh->SetAttributes();
 
    // 13. Setup the final NonlinearForm (which defines the integral of interest,
    //     its first and second derivatives). Here we can use a combination of

miniapps/meshing/pmesh-optimizer.cpp

@@ -35,7 +35,7 @@
 //   Adapted analytic shape:
 //     mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 2 -tid 4 -ni 200 -bnd -qt 1 -qo 8
 //   Adapted analytic size+orientation:
-//     mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 14 -tid 4 -ni 200 -bnd -qt 1 -qo 8 -fd
+//     mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 14 -tid 4 -ni 200 -bnd -qt 1 -qo 8
 //   Adapted analytic shape+orientation:
 //     mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 85 -tid 4 -ni 100 -bnd -qt 1 -qo 8 -fd
 //

miniapps/meshing/pminimal-surface.cpp

@@ -1019,7 +1019,7 @@ struct QuarterPeach: public Surface
       for (int i = 0; i < GetNBE(); i++)
       {
          Element *el = GetBdrElement(i);
-         const int fn = GetBdrElementEdgeIndex(i);
+         const int fn = GetBdrElementFaceIndex(i);
          MFEM_VERIFY(!FaceIsTrueInterior(fn),"");
          Array<int> vertices;
          GetFaceVertices(fn, vertices);

miniapps/multidomain/multidomain.cpp

@@ -322,17 +322,19 @@ int main(int argc, char *argv[])
 
    char vishost[] = "localhost";
    int  visport   = 19916;
-   socketstream cyl_sol_sock(vishost, visport);
+   socketstream cyl_sol_sock;
    if (visualization)
    {
+      cyl_sol_sock.open(vishost, visport);
       cyl_sol_sock << "parallel " << num_procs << " " << myid << "\n";
       cyl_sol_sock.precision(8);
       cyl_sol_sock << "solution\n" << cylinder_submesh << temperature_cylinder_gf <<
                    "pause\n" << std::flush;
    }
-   socketstream block_sol_sock(vishost, visport);
+   socketstream block_sol_sock;
    if (visualization)
    {
+      block_sol_sock.open(vishost, visport);
       block_sol_sock << "parallel " << num_procs << " " << myid << "\n";
       block_sol_sock.precision(8);
       block_sol_sock << "solution\n" << block_submesh << temperature_block_gf <<

miniapps/nurbs/nurbs_ex1.cpp

@@ -49,8 +49,8 @@ public:
    Data(double x_, double val_) {x=x_; val=val_;};
 };
 
-inline bool operator==(const Data& d1,const Data& d2) { return (d1.x == d2.x); };
-inline bool operator <(const Data& d1,const Data& d2) { return (d1.x  < d2.x); };
+inline bool operator==(const Data& d1,const Data& d2) { return (d1.x == d2.x); }
+inline bool operator <(const Data& d1,const Data& d2) { return (d1.x  < d2.x); }
 
 /** Class for integrating the bilinear form a(u,v) := (Q Laplace u, v) where Q
     can be a scalar coefficient. */

miniapps/nurbs/nurbs_ex11p.cpp

@@ -281,12 +281,13 @@ int main(int argc, char *argv[])
 #ifdef MFEM_USE_STRUMPACK
       if (sp_solver)
       {
-         STRUMPACKSolver * strumpack = new STRUMPACKSolver(argc, argv, MPI_COMM_WORLD);
+         STRUMPACKSolver * strumpack = new STRUMPACKSolver(MPI_COMM_WORLD, argc, argv);
          strumpack->SetPrintFactorStatistics(true);
          strumpack->SetPrintSolveStatistics(false);
          strumpack->SetKrylovSolver(strumpack::KrylovSolver::DIRECT);
          strumpack->SetReorderingStrategy(strumpack::ReorderingStrategy::METIS);
-         strumpack->DisableMatching();
+         strumpack->SetMatching(strumpack::MatchingJob::NONE);
+         strumpack->SetCompression(strumpack::CompressionType::NONE);
          strumpack->SetOperator(*Arow);
          strumpack->SetFromCommandLine();
          precond = strumpack;

miniapps/shifted/extrapolate.cpp

@@ -111,7 +111,7 @@ int main(int argc, char *argv[])
 {
    // Initialize MPI and HYPRE.
    Mpi::Init();
-   int myid = Mpi::WorldRank();;
+   int myid = Mpi::WorldRank();
    Hypre::Init();
 
    // Parse command-line options.

miniapps/shifted/marking.cpp

@@ -160,7 +160,7 @@ void ShiftedFaceMarker::ListShiftedFaceDofs(const Array<int> &elem_marker,
          {
             if (elem_marker[tr->Elem1No] >= SBElementType::CUT)
             {
-               pfes_sltn->GetFaceDofs(pmesh.GetBdrFace(i), dofs);
+               pfes_sltn->GetFaceDofs(pmesh.GetBdrElementFaceIndex(i), dofs);
                sface_dof_list.Append(dofs);
             }
          }
@@ -343,7 +343,7 @@ void ShiftedFaceMarker::ListShiftedFaceDofs2(const Array<int> &elem_marker,
          {
             if (elem_marker[tr->Elem1No] >= SBElementType::CUT)
             {
-               pfes_sltn->GetFaceDofs(pmesh.GetBdrFace(i), dofs);
+               pfes_sltn->GetFaceDofs(pmesh.GetBdrElementFaceIndex(i), dofs);
                sface_dof_list.Append(dofs);
             }
          }

miniapps/solvers/lor_solvers.cpp

@@ -57,6 +57,9 @@
 //    lor_solvers -m ../../data/amr-quad.mesh -fe n
 //    lor_solvers -m ../../data/amr-quad.mesh -fe r
 //    lor_solvers -m ../../data/amr-quad.mesh -fe l
+//    lor_solvers -m ../../data/star-surf.mesh -fe h
+//    lor_solvers -m ../../data/star-surf.mesh -fe n
+//    lor_solvers -m ../../data/star-surf.mesh -fe r
 //
 // Device sample runs:
 //    lor_solvers -fe h -d cuda
@@ -110,12 +113,15 @@ int main(int argc, char *argv[])
    double kappa = (order+1)*(order+1); // Penalty used for DG discretizations
 
    Mesh mesh(mesh_file, 1, 1);
-   int dim = mesh.Dimension();
-   MFEM_VERIFY(dim == 2 || dim == 3, "Spatial dimension must be 2 or 3.");
+   const int dim = mesh.Dimension();
+   const int sdim = mesh.SpaceDimension();
+   MFEM_VERIFY(dim == 2 || dim == 3, "Mesh dimension must be 2 or 3.");
+   MFEM_VERIFY(!L2 || dim == sdim, "DG surface meshes not supported.");
    for (int l = 0; l < ref_levels; l++) { mesh.UniformRefinement(); }
 
    FunctionCoefficient f_coeff(f(1.0)), u_coeff(u);
-   VectorFunctionCoefficient f_vec_coeff(dim, f_vec(RT)), u_vec_coeff(dim, u_vec);
+   VectorFunctionCoefficient f_vec_coeff(sdim, f_vec(RT)),
+                             u_vec_coeff(sdim, u_vec);
 
    int b1 = BasisType::GaussLobatto, b2 = BasisType::IntegratedGLL;
    unique_ptr<FiniteElementCollection> fec;
@@ -149,8 +155,9 @@ int main(int argc, char *argv[])
       a.AddInteriorFaceIntegrator(new DGDiffusionIntegrator(-1.0, kappa));
       a.AddBdrFaceIntegrator(new DGDiffusionIntegrator(-1.0, kappa));
    }
-   // TODO: L2 diffusion not implemented with partial assembly
-   if (!L2) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
+   // Partial assembly not currently supported for DG or for surface meshes with
+   // vector finite elements (ND or RT).
+   if (!L2 && (H1 || sdim == dim)) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
    a.Assemble();
 
    LinearForm b(&fes);
@@ -188,9 +195,12 @@ int main(int argc, char *argv[])
 
    a.RecoverFEMSolution(X, b, x);
 
-   double er =
-      (H1 || L2) ? x.ComputeL2Error(u_coeff) : x.ComputeL2Error(u_vec_coeff);
-   cout << "L2 error: " << er << endl;
+   if (sdim == dim)
+   {
+      double er =
+         (H1 || L2) ? x.ComputeL2Error(u_coeff) : x.ComputeL2Error(u_vec_coeff);
+      cout << "L2 error: " << er << endl;
+   }
 
    if (visualization)
    {

miniapps/solvers/plor_solvers.cpp

@@ -55,6 +55,9 @@
 //    mpirun -np 4 plor_solvers -m ../../data/fichera.mesh -fe l
 //    mpirun -np 4 plor_solvers -m ../../data/amr-hex.mesh -fe h -rs 0 -o 2
 //    mpirun -np 4 plor_solvers -m ../../data/amr-hex.mesh -fe l -rs 0 -o 2
+//    mpirun -np 4 plor_solvers -m ../../data/star-surf.mesh -fe h
+//    mpirun -np 4 plor_solvers -m ../../data/star-surf.mesh -fe n
+//    mpirun -np 4 plor_solvers -m ../../data/star-surf.mesh -fe r
 //
 // Device sample runs:
 //    mpirun -np 4 plor_solvers -m ../../data/fichera.mesh -fe h -d cuda
@@ -113,8 +116,10 @@ int main(int argc, char *argv[])
    double kappa = (order+1)*(order+1); // Penalty used for DG discretizations
 
    Mesh serial_mesh(mesh_file, 1, 1);
-   int dim = serial_mesh.Dimension();
-   MFEM_VERIFY(dim == 2 || dim == 3, "Spatial dimension must be 2 or 3.");
+   const int dim = serial_mesh.Dimension();
+   const int sdim = serial_mesh.SpaceDimension();
+   MFEM_VERIFY(dim == 2 || dim == 3, "Mesh dimension must be 2 or 3.");
+   MFEM_VERIFY(!L2 || dim == sdim, "DG surface meshes not supported.");
    for (int l = 0; l < ser_ref_levels; l++) { serial_mesh.UniformRefinement(); }
    ParMesh mesh(MPI_COMM_WORLD, serial_mesh);
    for (int l = 0; l < par_ref_levels; l++) { mesh.UniformRefinement(); }
@@ -124,7 +129,8 @@ int main(int argc, char *argv[])
    { MFEM_ABORT("LOR AMS and ADS solvers are not supported with AMR meshes."); }
 
    FunctionCoefficient f_coeff(f(1.0)), u_coeff(u);
-   VectorFunctionCoefficient f_vec_coeff(dim, f_vec(RT)), u_vec_coeff(dim, u_vec);
+   VectorFunctionCoefficient f_vec_coeff(sdim, f_vec(RT)),
+                             u_vec_coeff(sdim, u_vec);
 
    int b1 = BasisType::GaussLobatto, b2 = BasisType::IntegratedGLL;
    unique_ptr<FiniteElementCollection> fec;
@@ -159,8 +165,9 @@ int main(int argc, char *argv[])
       a.AddInteriorFaceIntegrator(new DGDiffusionIntegrator(-1.0, kappa));
       a.AddBdrFaceIntegrator(new DGDiffusionIntegrator(-1.0, kappa));
    }
-   // TODO: L2 diffusion not implemented with partial assembly
-   if (!L2) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
+   // Partial assembly not currently supported for DG or for surface meshes with
+   // vector finite elements (ND or RT).
+   if (!L2 && (H1 || sdim == dim)) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
    a.Assemble();
 
    ParLinearForm b(&fes);
@@ -181,7 +188,6 @@ int main(int argc, char *argv[])
    OperatorHandle A;
    a.FormLinearSystem(ess_dofs, x, b, A, X, B);
 
-
    unique_ptr<Solver> solv_lor;
    if (H1 || L2)
    {
@@ -207,9 +213,12 @@ int main(int argc, char *argv[])
 
    a.RecoverFEMSolution(X, b, x);
 
-   double er =
-      (H1 || L2) ? x.ComputeL2Error(u_coeff) : x.ComputeL2Error(u_vec_coeff);
-   if (Mpi::Root()) { cout << "L2 error: " << er << endl; }
+   if (sdim == dim)
+   {
+      double er =
+         (H1 || L2) ? x.ComputeL2Error(u_coeff) : x.ComputeL2Error(u_vec_coeff);
+      if (Mpi::Root()) { cout << "L2 error: " << er << endl; }
+   }
 
    if (visualization)
    {

miniapps/spde/spde_solver.cpp

@@ -237,7 +237,7 @@ void Boundary::AddInhomogeneousDirichletBoundaryCondition(int boundary,
 void Boundary::SetRobinCoefficient(double coefficient)
 {
    robin_coefficient = coefficient;
-};
+}
 
 double IntegrateBC(const ParGridFunction &x, const Array<int> &bdr,
                    double alpha, double beta, double gamma, double &glb_err)
@@ -549,7 +549,7 @@ void SPDESolver::SetupRandomFieldGenerator(int seed)
       new WhiteGaussianNoiseDomainLFIntegrator(fespace_ptr_->GetComm(), seed);
    b_wn = new ParLinearForm(fespace_ptr_);
    b_wn->AddDomainIntegrator(integ);
-};
+}
 
 void SPDESolver::GenerateRandomField(ParGridFunction &x)
 {
@@ -565,7 +565,7 @@ void SPDESolver::GenerateRandomField(ParGridFunction &x)
 
    // Call back to solve to generate the random field
    Solve(*b_wn, x);
-};
+}
 
 double SPDESolver::ConstructNormalizationCoefficient(double nu, double l1,
                                                      double l2, double l3,

miniapps/tools/CMakeLists.txt

@@ -33,8 +33,13 @@ if (MFEM_USE_MPI)
     MAIN plor-transfer.cpp LIBRARIES mfem)
 endif()
 
-add_mfem_miniapp(check-tmop-metric
-  MAIN check-tmop-metric.cpp LIBRARIES mfem)
+add_mfem_miniapp(tmop-check-metric
+  MAIN tmop-check-metric.cpp LIBRARIES mfem)
+
+add_mfem_miniapp(tmop-metric-magnitude
+  MAIN tmop-metric-magnitude.cpp
+  ${MFEM_MINIAPPS_COMMON_HEADERS}
+  LIBRARIES mfem mfem-common)
 
 if (MFEM_USE_MPI)
    add_mfem_miniapp(nodal-transfer

miniapps/tools/display-basis.cpp

@@ -206,7 +206,7 @@ int main(int argc, char *argv[])
       {
          cout << endl;
          cout << "Element Type:          " << elemTypeStr(eType) << endl;
-         cout << "Basis Type:            " << basisTypeStr(bType) << endl;;
+         cout << "Basis Type:            " << basisTypeStr(bType) << endl;
          cout << "Basis function order:  " << bOrder << endl;
          cout << "Map Type:              " << mapTypeStr(mType) << endl;
       }

miniapps/tools/makefile

@@ -26,7 +26,7 @@ MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
 SEQ_MINIAPPS = display-basis load-dc convert-dc get-values lor-transfer \
-	check-tmop-metric tmop-metric-magnitude
+	tmop-check-metric tmop-metric-magnitude
 
 PAR_MINIAPPS = nodal-transfer plor-transfer
 
@@ -84,7 +84,7 @@ RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
 # Testing: Specific execution options
 # Do not test: display-basis, load-dc, convert-dc, get-values, lor-transfer, plor-transfer
 NO_TEST_APPS = display-basis load-dc convert-dc get-values lor-transfer \
-	plor-transfer check-tmop-metric tmop-metric-magnitude
+	plor-transfer tmop-check-metric tmop-metric-magnitude
 $(foreach app,$(NO_TEST_APPS),$(app)-test-seq $(app)-test-par):
 	@true
 

miniapps/tools/tmop-check-metric.cpp

@@ -16,9 +16,9 @@
 // This miniapp checks the evaluation, 1st, and 2nd derivatives of a TMOP
 // metric. Works only in serial.
 //
-// Compile with: make check-tmop-metric
+// Compile with: make tmop-check-metric
 //
-// Sample runs:  check-tmop-metric -mid 360
+// Sample runs:  tmop-check-metric -mid 360
 
 #include "mfem.hpp"
 #include <iostream>

tests/unit/fem/test_doftrans.cpp

@@ -22,13 +22,14 @@ TEST_CASE("DoF Transformation Classes",
           "[ND_TetDofTransformation]")
 {
    int p = 4;
+   int vdim = 3;
    int seed = 123;
 
    double tol = 1e-13;
 
    SECTION("Nedelec Tetrahedral Transformations")
    {
-      ND_TetDofTransformation T(p);
+      ND_TetDofTransformation Tnd(p);
 
       Array<int> ori(4);
       ori[0] = 1;
@@ -36,102 +37,191 @@ TEST_CASE("DoF Transformation Classes",
       ori[2] = 5;
       ori[3] = 1;
 
-      T.SetFaceOrientations(ori);
-
-      Vector u(T.Width());
-      Vector v(T.Width());
-      Vector f(T.Width());
-      Vector ut;
-      Vector vt;
-      Vector ft;
-
-      u.Randomize(seed);
-      v.Randomize(seed+1);
-      f.Randomize(seed+2);
-
-      SECTION("Inverse DoF transformation")
+      SECTION("VDim == 1")
       {
-         Vector w;
+         DofTransformation T(Tnd);
+         T.SetFaceOrientations(ori);
 
-         ut = u; T.TransformPrimal(ut);
-         w = ut; T.InvTransformPrimal(w);
+         Vector u(T.Width());
+         Vector v(T.Width());
+         Vector f(T.Width());
+         Vector ut;
+         Vector vt;
+         Vector ft;
 
-         w -= u;
+         u.Randomize(seed);
+         v.Randomize(seed+1);
+         f.Randomize(seed+2);
 
-         REQUIRE(w.Norml2() < tol * u.Norml2());
-      }
-      SECTION("Inverse Dual DoF transformation")
-      {
-         Vector w;
-
-         ut = u; T.TransformDual(ut);
-         w = ut; T.InvTransformDual(w);
-
-         w -= u;
-
-         REQUIRE(w.Norml2() < tol * u.Norml2());
-      }
-
-      SECTION("Inner product with linear form f(v)")
-      {
-         vt = v; T.TransformPrimal(vt);
-         ft = f; T.TransformDual(ft);
-
-         double fv = f * v;
-
-         REQUIRE(fabs(fv - ft * vt) < tol * fabs(fv));
-      }
-
-      DenseMatrix A(T.Width());
-      {
-         Vector Ac;
-         for (int i=0; i<A.Width(); i++)
+         SECTION("Inverse DoF transformation")
          {
-            A.GetColumnReference(i, Ac);
-            Ac.Randomize(seed+i);
+            Vector w;
+
+            ut = u; T.TransformPrimal(ut);
+            w = ut; T.InvTransformPrimal(w);
+
+            w -= u;
+
+            REQUIRE(w.Norml2() < tol * u.Norml2());
+         }
+         SECTION("Inverse Dual DoF transformation")
+         {
+            Vector w;
+
+            ut = u; T.TransformDual(ut);
+            w = ut; T.InvTransformDual(w);
+
+            w -= u;
+
+            REQUIRE(w.Norml2() < tol * u.Norml2());
+         }
+
+         SECTION("Inner product with linear form f(v)")
+         {
+            vt = v; T.TransformPrimal(vt);
+            ft = f; T.TransformDual(ft);
+
+            double fv = f * v;
+
+            REQUIRE(fabs(fv - ft * vt) < tol * fabs(fv));
+         }
+
+         DenseMatrix A(T.Width());
+         {
+            Vector Ac;
+            for (int i=0; i<A.Width(); i++)
+            {
+               A.GetColumnReference(i, Ac);
+               Ac.Randomize(seed+i);
+            }
+         }
+
+         SECTION("Inner product of two primal vectors")
+         {
+            // The matrix A in this case should be regarded as a BilinearForm.
+            DenseMatrix tA;
+            DenseMatrix At;
+            DenseMatrix tAt;
+
+            ut = u; T.TransformPrimal(ut);
+            vt = v; T.TransformPrimal(vt);
+
+            At = A; T.TransformDualRows(At);
+            tA = A; T.TransformDualCols(tA);
+            tAt = A; T.TransformDual(tAt);
+
+            double uAv = A.InnerProduct(v, u);
+
+            REQUIRE(fabs(uAv -  At.InnerProduct(vt, u )) < tol * fabs(uAv));
+            REQUIRE(fabs(uAv -  tA.InnerProduct(v, ut)) < tol * fabs(uAv));
+            REQUIRE(fabs(uAv - tAt.InnerProduct(vt, ut)) < tol * fabs(uAv));
+         }
+         SECTION("Inner product of a primal vector and a dual vector")
+         {
+            // The matrix A in this case should be regarded as a
+            // DiscreteLinearOperator.
+            DenseMatrix tA;
+            DenseMatrix At;
+            DenseMatrix tAt;
+
+            ft = f; T.TransformDual(ft);
+            vt = v; T.TransformPrimal(vt);
+
+            At = A; T.TransformDualRows(At);
+            tA = A; T.TransformPrimalCols(tA);
+            tAt = At; T.TransformPrimalCols(tAt);
+
+            double fAv = A.InnerProduct(v, f);
+
+            REQUIRE(fabs(fAv -  At.InnerProduct(vt, f )) < tol * fabs(fAv));
+            REQUIRE(fabs(fAv -  tA.InnerProduct(v, ft)) < tol * fabs(fAv));
+            REQUIRE(fabs(fAv - tAt.InnerProduct(vt, ft)) < tol * fabs(fAv));
          }
       }
-
-      SECTION("Inner product of two primal vectors")
+      SECTION("VDim > 1")
       {
-         // The matrix A in this case should be regarded as a BilinearForm.
-         DenseMatrix tA;
-         DenseMatrix At;
-         DenseMatrix tAt;
+         Vector v(vdim * Tnd.Width());
+         Vector f(vdim * Tnd.Width());
+         Vector vt;
+         Vector ft;
 
-         ut = u; T.TransformPrimal(ut);
-         vt = v; T.TransformPrimal(vt);
+         v.Randomize(seed);
+         f.Randomize(seed+1);
 
-         At = A; T.TransformDualRows(At);
-         tA = A; T.TransformDualCols(tA);
-         tAt = A; T.TransformDual(tAt);
+         SECTION("Ordering == byNODES")
+         {
+            DofTransformation T(Tnd, vdim, Ordering::byNODES);
+            T.SetFaceOrientations(ori);
 
-         double uAv = A.InnerProduct(v, u);
+            SECTION("Inverse DoF transformation")
+            {
+               Vector w;
 
-         REQUIRE(fabs(uAv -  At.InnerProduct(vt, u )) < tol * fabs(uAv));
-         REQUIRE(fabs(uAv -  tA.InnerProduct(v, ut)) < tol * fabs(uAv));
-         REQUIRE(fabs(uAv - tAt.InnerProduct(vt, ut)) < tol * fabs(uAv));
-      }
-      SECTION("Inner product of a primal vector and a dual vector")
-      {
-         // The matrix A in this case should be regarded as a
-         // DiscreteLinearOperator.
-         DenseMatrix tA;
-         DenseMatrix At;
-         DenseMatrix tAt;
+               vt = v; T.TransformPrimal(vt);
+               w = vt; T.InvTransformPrimal(w);
 
-         ft = f; T.TransformDual(ft);
-         vt = v; T.TransformPrimal(vt);
+               w -= v;
 
-         At = A; T.TransformDualRows(At);
-         tA = A; T.TransformPrimalCols(tA);
-         tAt = At; T.TransformPrimalCols(tAt);
+               REQUIRE(w.Norml2() < tol * v.Norml2());
+            }
+            SECTION("Inverse Dual DoF transformation")
+            {
+               Vector w;
 
-         double fAv = A.InnerProduct(v, f);
+               vt = v; T.TransformDual(vt);
+               w = vt; T.InvTransformDual(w);
 
-         REQUIRE(fabs(fAv -  At.InnerProduct(vt, f )) < tol * fabs(fAv));
-         REQUIRE(fabs(fAv -  tA.InnerProduct(v, ft)) < tol * fabs(fAv));
-         REQUIRE(fabs(fAv - tAt.InnerProduct(vt, ft)) < tol * fabs(fAv));
+               w -= v;
+
+               REQUIRE(w.Norml2() < tol * v.Norml2());
+            }
+            SECTION("Inner product with linear form f(v)")
+            {
+               vt = v; T.TransformPrimal(vt);
+               ft = f; T.TransformDual(ft);
+
+               double fv = f * v;
+
+               REQUIRE(fabs(fv - ft * vt) < tol * fabs(fv));
+            }
+         }
+         SECTION("Ordering == byVDIM")
+         {
+            DofTransformation T(Tnd, vdim, Ordering::byVDIM);
+            T.SetFaceOrientations(ori);
+
+            SECTION("Inverse DoF transformation")
+            {
+               Vector w;
+
+               vt = v; T.TransformPrimal(vt);
+               w = vt; T.InvTransformPrimal(w);
+
+               w -= v;
+
+               REQUIRE(w.Norml2() < tol * v.Norml2());
+            }
+            SECTION("Inverse Dual DoF transformation")
+            {
+               Vector w;
+
+               vt = v; T.TransformDual(vt);
+               w = vt; T.InvTransformDual(w);
+
+               w -= v;
+
+               REQUIRE(w.Norml2() < tol * v.Norml2());
+            }
+            SECTION("Inner product with linear form f(v)")
+            {
+               vt = v; T.TransformPrimal(vt);
+               ft = f; T.TransformDual(ft);
+
+               double fv = f * v;
+
+               REQUIRE(fabs(fv - ft * vt) < tol * fabs(fv));
+            }
+         }
       }
    }
 }
@@ -146,8 +236,8 @@ TEST_CASE("DoF Transformation Functions",
 
    double tol = 1e-13;
 
-   ND_TetDofTransformation Tp(p);
-   ND_TetDofTransformation Tq(q);
+   ND_TetDofTransformation Tndp(p);
+   ND_TetDofTransformation Tndq(q);
 
    Array<int> ori(4);
    ori[0] = 1;
@@ -155,6 +245,7 @@ TEST_CASE("DoF Transformation Functions",
    ori[2] = 5;
    ori[3] = 1;
 
+   DofTransformation Tp(Tndp), Tq(Tndq);
    Tp.SetFaceOrientations(ori);
    Tq.SetFaceOrientations(ori);
 
@@ -235,153 +326,4 @@ TEST_CASE("DoF Transformation Functions",
    }
 }
 
-TEST_CASE("VDoF Transformation Class",
-          "[DofTransformation]"
-          "[VDofTransformation]")
-{
-   int p = 4;
-   int vdim = 3;
-   int seed = 123;
-
-   double tol = 1e-13;
-
-   ND_TetDofTransformation Tnd(p);
-
-   Array<int> ori(4);
-   ori[0] = 1;
-   ori[1] = 3;
-   ori[2] = 5;
-   ori[3] = 1;
-
-   Tnd.SetFaceOrientations(ori);
-
-   SECTION("VDim == 1")
-   {
-      VDofTransformation T(Tnd);
-
-      Vector v(T.Width());
-      Vector f(T.Width());
-      Vector vt;
-      Vector ft;
-
-      v.Randomize(seed);
-      f.Randomize(seed+1);
-
-      SECTION("Inverse DoF transformation")
-      {
-         Vector w;
-
-         vt = v; T.TransformPrimal(vt);
-         w = vt; T.InvTransformPrimal(w);
-
-         w -= v;
-
-         REQUIRE(w.Norml2() < tol * v.Norml2());
-      }
-      SECTION("Inverse Dual DoF transformation")
-      {
-         Vector w;
-
-         vt = v; T.TransformDual(vt);
-         w = vt; T.InvTransformDual(w);
-
-         w -= v;
-
-         REQUIRE(w.Norml2() < tol * v.Norml2());
-      }
-      SECTION("Inner product with linear form f(v)")
-      {
-         vt = v; T.TransformPrimal(vt);
-         ft = f; T.TransformDual(ft);
-
-         double fv = f * v;
-
-         REQUIRE(fabs(fv - ft * vt) < tol * fabs(fv));
-      }
-   }
-   SECTION("VDim > 1")
-   {
-      Vector v(vdim * Tnd.Width());
-      Vector f(vdim * Tnd.Width());
-      Vector vt;
-      Vector ft;
-
-      v.Randomize(seed);
-      f.Randomize(seed+1);
-
-      SECTION("Ordering == byNODES")
-      {
-         VDofTransformation T(Tnd, vdim, Ordering::byNODES);
-
-         SECTION("Inverse DoF transformation")
-         {
-            Vector w;
-
-            vt = v; T.TransformPrimal(vt);
-            w = vt; T.InvTransformPrimal(w);
-
-            w -= v;
-
-            REQUIRE(w.Norml2() < tol * v.Norml2());
-         }
-         SECTION("Inverse Dual DoF transformation")
-         {
-            Vector w;
-
-            vt = v; T.TransformDual(vt);
-            w = vt; T.InvTransformDual(w);
-
-            w -= v;
-
-            REQUIRE(w.Norml2() < tol * v.Norml2());
-         }
-         SECTION("Inner product with linear form f(v)")
-         {
-            vt = v; T.TransformPrimal(vt);
-            ft = f; T.TransformDual(ft);
-
-            double fv = f * v;
-
-            REQUIRE(fabs(fv - ft * vt) < tol * fabs(fv));
-         }
-      }
-      SECTION("Ordering == byVDIM")
-      {
-         VDofTransformation T(Tnd, vdim, Ordering::byVDIM);
-
-         SECTION("Inverse DoF transformation")
-         {
-            Vector w;
-
-            vt = v; T.TransformPrimal(vt);
-            w = vt; T.InvTransformPrimal(w);
-
-            w -= v;
-
-            REQUIRE(w.Norml2() < tol * v.Norml2());
-         }
-         SECTION("Inverse Dual DoF transformation")
-         {
-            Vector w;
-
-            vt = v; T.TransformDual(vt);
-            w = vt; T.InvTransformDual(w);
-
-            w -= v;
-
-            REQUIRE(w.Norml2() < tol * v.Norml2());
-         }
-         SECTION("Inner product with linear form f(v)")
-         {
-            vt = v; T.TransformPrimal(vt);
-            ft = f; T.TransformDual(ft);
-
-            double fv = f * v;
-
-            REQUIRE(fabs(fv - ft * vt) < tol * fabs(fv));
-         }
-      }
-   }
-}
-
 } // namespace doftrans

tests/unit/fem/test_estimator.cpp

@@ -244,7 +244,7 @@ TEST_CASE("Least-squares ZZ estimator on 3D NCMesh", "[NCMesh]")
       {
          REQUIRE(local_errors(i) < 1e-10);
       }
-      REQUIRE(estimator.GetTotalError() < 1e-10);;
+      REQUIRE(estimator.GetTotalError() < 1e-10);
    }
 
    SECTION("Perfect Approximation Y")

tests/unit/fem/test_lor_batched.cpp

@@ -91,6 +91,7 @@ void TestBatchedLOR()
 {
    const int order = 5;
    const auto mesh_fname = GENERATE(
+                              "../../data/star-surf.mesh",
                               "../../data/star-q3.mesh",
                               "../../data/fichera-q3.mesh"
                            );
@@ -124,8 +125,11 @@ void TestBatchedLOR()
    IntegrationRules irs(0, Quadrature1D::GaussLobatto);
    const IntegrationRule &ir = irs.Get(mesh.GetElementGeometry(0), 1);
    const GeometricFactors::FactorFlags dets = GeometricFactors::DETERMINANTS;
-   REQUIRE(lor.GetFESpace().GetMesh()->GetGeometricFactors(ir, dets)->detJ.Min()
-           > 0.0);
+   if (mesh.Dimension() == mesh.SpaceDimension())
+   {
+      REQUIRE(
+         lor.GetFESpace().GetMesh()->GetGeometricFactors(ir, dets)->detJ.Min() > 0.0);
+   }
 
    lor.LegacyAssembleSystem(a, ess_dofs);
    SparseMatrix A1 = lor.GetAssembledMatrix(); // deep copy
@@ -183,6 +187,7 @@ void ParTestBatchedLOR()
    const bool all_tests = launch_all_non_regression_tests;
    const int order = !all_tests ? 5 : GENERATE(1,3,5);
    const auto mesh_fname = GENERATE(
+                              "../../data/star-surf.mesh",
                               "../../data/star-q3.mesh",
                               "../../data/fichera-q3.mesh"
                            );
@@ -241,6 +246,7 @@ TEST_CASE("LOR AMS", "[LOR][BatchedLOR][AMS][Parallel][CUDA]")
    enum SpaceType { ND, RT };
    auto space_type = GENERATE(ND, RT);
    auto mesh_fname = GENERATE(
+                        "../../data/star-surf.mesh",
                         "../../data/star-q3.mesh",
                         "../../data/fichera-q3.mesh"
                      );
@@ -251,6 +257,7 @@ TEST_CASE("LOR AMS", "[LOR][BatchedLOR][AMS][Parallel][CUDA]")
    serial_mesh.Clear();
 
    const int dim = mesh.Dimension();
+   const int sdim = mesh.SpaceDimension();
 
    // Only test RT spaces in 2D
    if (space_type == RT && dim == 3) { return; }
@@ -288,18 +295,18 @@ TEST_CASE("LOR AMS", "[LOR][BatchedLOR][AMS][Parallel][CUDA]")
       const double *coord = edge_fespace.GetMesh()->GetVertex(i);
       x_coord(i) = coord[0];
       y_coord(i) = coord[1];
-      if (dim == 3) { z_coord(i) = coord[2]; }
+      if (sdim == 3) { z_coord(i) = coord[2]; }
    }
    std::unique_ptr<HypreParVector> x(x_coord.ParallelProject());
    std::unique_ptr<HypreParVector> y(y_coord.ParallelProject());
    std::unique_ptr<HypreParVector> z;
-   if (dim == 3) { z.reset(z_coord.ParallelProject()); }
+   if (sdim == 3) { z.reset(z_coord.ParallelProject()); }
 
    *x -= *batched_lor.GetXCoordinate();
    REQUIRE(x->Normlinf() == MFEM_Approx(0.0));
    *y -= *batched_lor.GetYCoordinate();
    REQUIRE(y->Normlinf() == MFEM_Approx(0.0));
-   if (dim == 3)
+   if (sdim == 3)
    {
       *z -= *batched_lor.GetZCoordinate();
       REQUIRE(z->Normlinf() == MFEM_Approx(0.0));

tests/unit/linalg/test_direct_solvers.cpp

@@ -26,6 +26,9 @@ using namespace mfem;
 #ifdef MFEM_USE_SUPERLU
 #define DIRECT_SOLVE_PARALLEL
 #endif
+#ifdef MFEM_USE_STRUMPACK
+#define DIRECT_SOLVE_PARALLEL
+#endif
 
 #if defined(DIRECT_SOLVE_SERIAL) || defined(DIRECT_SOLVE_PARALLEL)
 
@@ -103,7 +106,7 @@ TEST_CASE("Serial Direct Solvers", "[CUDA]")
       Mesh mesh;
       if (dim == 1)
       {
-         mesh = Mesh::MakeCartesian1D(ne,  1.0);
+         mesh = Mesh::MakeCartesian1D(ne, 1.0);
       }
       else if (dim == 2)
       {
@@ -187,13 +190,13 @@ TEST_CASE("Parallel Direct Solvers", "[Parallel], [CUDA]")
 {
    int rank;
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
-   const int ne = 2;
+   const int ne = 4;
    for (int dim = 1; dim < 4; ++dim)
    {
       Mesh mesh;
       if (dim == 1)
       {
-         mesh = Mesh::MakeCartesian1D(ne,  1.0);
+         mesh = Mesh::MakeCartesian1D(ne, 1.0);
       }
       else if (dim == 2)
       {
@@ -312,6 +315,39 @@ TEST_CASE("Parallel Direct Solvers", "[Parallel], [CUDA]")
          REQUIRE(error < 1.e-12);
       }
 #endif
+#ifdef MFEM_USE_STRUMPACK
+      // Transform to monolithic HypreParMatrix
+      {
+         STRUMPACKRowLocMatrix SA(*A.As<HypreParMatrix>());
+         STRUMPACKSolver strumpack(MPI_COMM_WORLD);
+         strumpack.SetPrintFactorStatistics(false);
+         strumpack.SetPrintSolveStatistics(false);
+         strumpack.SetKrylovSolver(strumpack::KrylovSolver::DIRECT);
+         strumpack.SetReorderingStrategy(dim > 1 ? strumpack::ReorderingStrategy::METIS :
+                                         strumpack::ReorderingStrategy::NATURAL);
+         strumpack.SetOperator(SA);
+         strumpack.Mult(B, X);
+
+         Vector Y(X.Size());
+         A->Mult(X, Y);
+         Y -= B;
+         REQUIRE(Y.Norml2() < 1.e-12);
+
+         strumpack.ArrayMult(BB, XX);
+
+         for (int i = 0; i < XX.Size(); i++)
+         {
+            A->Mult(*XX[i], Y);
+            Y -= *BB[i];
+            REQUIRE(Y.Norml2() < 1.e-12);
+         }
+
+         a.RecoverFEMSolution(X, b, x);
+         VectorFunctionCoefficient grad(dim, gradexact);
+         double error = x.ComputeH1Error(&uex, &grad);
+         REQUIRE(error < 1.e-12);
+      }
+#endif
    }
 }
 

tests/unit/mesh/test_ncmesh.cpp

@@ -394,7 +394,7 @@ std::array<double, 2> CheckL2Projection(ParMesh& pmesh, Mesh& smesh, int order,
    }();
 
    return {serror, perror};
-};
+}
 
 TEST_CASE("EdgeFaceConstraint",  "[Parallel], [NCMesh]")
 {