Rule of zero for GridFunction and ParGridFunction

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

.gitignore

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

CHANGELOG

@@ -87,8 +87,6 @@ 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,34 +641,16 @@ 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}.tmp"
+    file(WRITE "${PROJECT_BINARY_DIR}/${Header}"
 "// 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}.tmp"
+    file(WRITE "${PROJECT_BINARY_DIR}/InstallHeaders/${Header}"
 "// 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,7 +659,8 @@ 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.
+  The support for STRUMPACK was added in MFEM v3.3.2 and it requires STRUMPACK
+  2.0.0 or later.
   URL: http://portal.nersc.gov/project/sparse/strumpack
   Options: STRUMPACK_OPT, STRUMPACK_LIB.
   Versions: STRUMPACK >= 3.0.0.

config/defaults.cmake

@@ -157,8 +157,7 @@ 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"
-    "Scotch/ptscotch/ptscotcherr/scotch/scotcherr"
-    "ScaLAPACK" "LAPACK" "BLAS" CACHE STRING
+    "ScaLAPACK" "Scotch/ptscotch/ptscotcherr/scotch/scotcherr" 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,13 +262,12 @@ int main(int argc, char *argv[])
 #ifdef MFEM_USE_STRUMPACK
       if (sp_solver)
       {
-         STRUMPACKSolver * strumpack = new STRUMPACKSolver(MPI_COMM_WORLD, argc, argv);
+         STRUMPACKSolver * strumpack = new STRUMPACKSolver(argc, argv, MPI_COMM_WORLD);
          strumpack->SetPrintFactorStatistics(true);
          strumpack->SetPrintSolveStatistics(false);
          strumpack->SetKrylovSolver(strumpack::KrylovSolver::DIRECT);
          strumpack->SetReorderingStrategy(strumpack::ReorderingStrategy::METIS);
-         strumpack->SetMatching(strumpack::MatchingJob::NONE);
-         strumpack->SetCompression(strumpack::CompressionType::NONE);
+         strumpack->DisableMatching();
          strumpack->SetOperator(*Arow);
          strumpack->SetFromCommandLine();
          precond = strumpack;

examples/ex25p.cpp

@@ -170,7 +170,6 @@ 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";
@@ -201,11 +200,6 @@ 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",
@@ -215,14 +209,13 @@ 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 + strumpack_solver > 1)
+   if (slu_solver && mumps_solver)
    {
       if (myid == 0)
-         cout << "WARNING: More than one of SuperLU, MUMPS, and STRUMPACK have"
-              << " been selected, please choose only one." << endl
+         cout << "WARNING: Both SuperLU and MUMPS have been selected,"
+              << " please choose either one." << endl
               << "         Defaulting to SuperLU." << endl;
       mumps_solver = false;
-      strumpack_solver = false;
    }
 
    if (iprob > 4) { iprob = 4; }
@@ -481,24 +474,6 @@ 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)
    {
@@ -518,7 +493,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 && !strumpack_solver))
+   if (pa || (!slu_solver && !mumps_solver))
    {
       ConstantCoefficient absomeg(pow(omega, 2) * epsilon);
       RestrictedCoefficient restr_absomeg(absomeg,attr);

examples/petsc/ex11p.cpp

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

fem/CMakeLists.txt

@@ -96,6 +96,9 @@ 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,7 +101,6 @@ 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;
@@ -434,6 +433,7 @@ 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,9 +441,6 @@ 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++)
             {
@@ -451,8 +448,9 @@ void BilinearForm::Assemble(int skip_zeros)
                     (*(domain_integs_marker[k]))[elem_attr-1] == 1)
                    && !domain_integs[k]->Patchwise())
                {
-                  domain_integs[k]->AssembleElementMatrix(*fes->GetFE(i),
-                                                          *eltrans, elemmat);
+                  const FiniteElement &fe = *fes->GetFE(i);
+                  eltrans = fes->GetElementTransformation(i);
+                  domain_integs[k]->AssembleElementMatrix(fe, *eltrans, elemmat);
                   if (elmat.Size() == 0)
                   {
                      elmat = elemmat;
@@ -1224,14 +1222,11 @@ MixedBilinearForm::MixedBilinearForm (FiniteElementSpace *tr_fes,
 
    // Copy the pointers to the integrators
    domain_integs = mbf->domain_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_integs_marker = mbf->boundary_integs_marker;
    boundary_trace_face_integs_marker = mbf->boundary_trace_face_integs_marker;
 
    assembly = AssemblyLevel::LEGACY;
@@ -1354,14 +1349,6 @@ 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)
@@ -1396,7 +1383,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)
    {
@@ -1418,20 +1405,8 @@ 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);
@@ -1440,14 +1415,10 @@ void MixedBilinearForm::Assemble(int skip_zeros)
          elmat = 0.0;
          for (int k = 0; k < domain_integs.Size(); k++)
          {
-            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;
-            }
+            domain_integs[k] -> AssembleElementMatrix2 (*trial_fes -> GetFE(i),
+                                                        *test_fes  -> GetFE(i),
+                                                        *eltrans, elemmat);
+            elmat += elemmat;
          }
          if (ran_dof_trans || dom_dof_trans)
          {
@@ -1970,56 +1941,41 @@ void DiscreteLinearOperator::Assemble(int skip_zeros)
       return;
    }
 
-   ElementTransformation *eltrans;
+   Array<int> dom_vdofs, ran_vdofs;
+   ElementTransformation *T;
    DofTransformation * dom_dof_trans;
    DofTransformation * ran_dof_trans;
-   DenseMatrix elmat;
-
-   Mesh *mesh = test_fes->GetMesh();
+   const FiniteElement *dom_fe, *ran_fe;
+   DenseMatrix totelmat, elmat;
 
    if (mat == NULL)
    {
       mat = new SparseMatrix(height, width);
    }
 
-   if (domain_integs.Size())
+   if (domain_integs.Size() > 0)
    {
-      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);
+         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);
 
-         elmat.SetSize(test_vdofs.Size(), trial_vdofs.Size());
-         elmat = 0.0;
-         for (int k = 0; k < domain_integs.Size(); k++)
+         domain_integs[0]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T,
+                                                  totelmat);
+         for (int j = 1; j < domain_integs.Size(); j++)
          {
-            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;
-            }
+            domain_integs[j]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T,
+                                                     elmat);
+            totelmat += elmat;
          }
          if (ran_dof_trans || dom_dof_trans)
          {
-            TransformPrimal(ran_dof_trans, dom_dof_trans, elemmat);
+            TransformPrimal(ran_dof_trans, dom_dof_trans, totelmat);
          }
-         mat->SetSubMatrix(test_vdofs, trial_vdofs, elemmat, skip_zeros);
+         mat->SetSubMatrix(ran_vdofs, dom_vdofs, totelmat, skip_zeros);
       }
    }
 
@@ -2028,20 +1984,21 @@ void DiscreteLinearOperator::Assemble(int skip_zeros)
       const int nfaces = test_fes->GetMesh()->GetNumFaces();
       for (int i = 0; i < nfaces; i++)
       {
-         trial_fes->GetFaceVDofs(i, trial_vdofs);
-         test_fes->GetFaceVDofs(i, test_vdofs);
-         eltrans = test_fes->GetMesh()->GetFaceTransformation(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);
 
-         elmat.SetSize(test_vdofs.Size(), trial_vdofs.Size());
-         elmat = 0.0;
-         for (int k = 0; k < trace_face_integs.Size(); k++)
+         trace_face_integs[0]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T,
+                                                      totelmat);
+         for (int j = 1; j < trace_face_integs.Size(); j++)
          {
-            trace_face_integs[k]->AssembleElementMatrix2(*trial_fes->GetFaceElement(i),
-                                                         *test_fes->GetFaceElement(i),
-                                                         *eltrans, elemmat);
-            elmat += elemmat;
+            trace_face_integs[j]->AssembleElementMatrix2(*dom_fe, *ran_fe, *T,
+                                                         elmat);
+            totelmat += elmat;
          }
-         mat->SetSubMatrix(test_vdofs, trial_vdofs, elmat, skip_zeros);
+         mat->SetSubMatrix(ran_vdofs, dom_vdofs, totelmat, 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; ///< Entries are not owned.
+   Array<Array<int>*>             domain_integs_marker;
 
    /// Set of Boundary Integrators to be applied.
    Array<BilinearFormIntegrator*> boundary_integs;
@@ -722,13 +722,10 @@ protected:
 
    /// Domain integrators.
    Array<BilinearFormIntegrator*> domain_integs;
-   /// Entries are not owned.
-   Array<Array<int>*> domain_integs_marker;
 
    /// Boundary integrators.
    Array<BilinearFormIntegrator*> boundary_integs;
-   /// Entries are not owned.
-   Array<Array<int>*> boundary_integs_marker;
+   Array<Array<int>*> boundary_integs_marker; ///< Entries are not owned.
 
    /// Trace face (skeleton) integrators.
    Array<BilinearFormIntegrator*> trace_face_integs;
@@ -808,16 +805,12 @@ 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.
 
@@ -827,18 +820,14 @@ 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; }
@@ -1076,9 +1065,6 @@ 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)
@@ -1086,7 +1072,6 @@ 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.GetBdrElementFaceIndex(i);
+         const int f = mesh.GetBdrElementEdgeIndex(i);
          f_to_be[f] = i;
       }
       const int nf_bdr = trial_fes->GetNFbyType(FaceType::Boundary);

fem/doftrans.cpp

@@ -14,166 +14,175 @@
 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)
 {
-   // No action if both transformations are NULL
-   if (ran_dof_trans)
+   if (ran_dof_trans && dom_dof_trans)
+   {
+      ran_dof_trans->TransformPrimalCols(elmat);
+      dom_dof_trans->TransformDualRows(elmat);
+   }
+   else if (ran_dof_trans)
    {
       ran_dof_trans->TransformPrimalCols(elmat);
    }
-   if (dom_dof_trans)
+   else 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)
 {
-   // No action if both transformations are NULL
-   if (ran_dof_trans)
+   if (ran_dof_trans && dom_dof_trans)
+   {
+      ran_dof_trans->TransformDualCols(elmat);
+      dom_dof_trans->TransformDualRows(elmat);
+   }
+   else if (ran_dof_trans)
    {
       ran_dof_trans->TransformDualCols(elmat);
    }
-   if (dom_dof_trans)
+   else 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_DofTransformation::T_data[24] =
+const double ND_StatelessDofTransformation::T_data[24] =
 {
    1.0,  0.0,  0.0,  1.0,
    -1.0, -1.0,  0.0,  1.0,
@@ -183,11 +192,11 @@ const double ND_DofTransformation::T_data[24] =
    0.0,  1.0,  1.0,  0.0
 };
 
-const DenseTensor ND_DofTransformation
-::T(const_cast<double *>(ND_DofTransformation::T_data), 2, 2, 6);
+const DenseTensor ND_StatelessDofTransformation
+::T(const_cast<double*>(ND_StatelessDofTransformation::T_data), 2, 2, 6);
 
 // ordering (i0j0, i1j0, i0j1, i1j1), each row is a column major matrix
-const double ND_DofTransformation::TInv_data[24] =
+const double ND_StatelessDofTransformation::TInv_data[24] =
 {
    1.0,  0.0,  0.0,  1.0,
    -1.0, -1.0,  0.0,  1.0,
@@ -197,11 +206,12 @@ const double ND_DofTransformation::TInv_data[24] =
    0.0,  1.0,  1.0,  0.0
 };
 
-const DenseTensor ND_DofTransformation
-::TInv(const_cast<double *>(TInv_data), 2, 2, 6);
+const DenseTensor ND_StatelessDofTransformation
+::TInv(const_cast<double*>(TInv_data), 2, 2, 6);
 
-ND_DofTransformation::ND_DofTransformation(int size, int p, int num_edges,
-                                           int num_tri_faces)
+ND_StatelessDofTransformation::ND_StatelessDofTransformation(int size, int p,
+                                                             int num_edges,
+                                                             int num_tri_faces)
    : StatelessDofTransformation(size)
    , order(p)
    , nedofs(p)
@@ -211,19 +221,18 @@ ND_DofTransformation::ND_DofTransformation(int size, int p, int num_edges,
 {
 }
 
-void ND_DofTransformation::TransformPrimal(const Array<int> & Fo,
-                                           double *v) const
+void ND_StatelessDofTransformation::TransformPrimal(const Array<int> & Fo,
+                                                    double *v) const
 {
    // Return immediately when no face DoFs are present
-   if (IsIdentity()) { return; }
+   if (nfdofs < 2) { return; }
 
    MFEM_VERIFY(Fo.Size() >= nfaces,
                "Face orientation array is shorter than the number of faces in "
-               "ND_DofTransformation");
+               "ND_StatelessDofTransformation");
 
    double data[2];
    Vector v2(data, 2);
-   DenseMatrix T2;
 
    // Transform face DoFs
    for (int f=0; f<nfaces; f++)
@@ -231,25 +240,23 @@ void ND_DofTransformation::TransformPrimal(const Array<int> & Fo,
       for (int i=0; i<nfdofs/2; i++)
       {
          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]);
+         T(Fo[f]).Mult(v2, &v[nedges*nedofs + f*nfdofs + 2*i]);
       }
    }
 }
 
-void ND_DofTransformation::InvTransformPrimal(const Array<int> & Fo,
-                                              double *v) const
+void ND_StatelessDofTransformation::InvTransformPrimal(const Array<int> & Fo,
+                                                       double *v) const
 {
    // Return immediately when no face DoFs are present
-   if (IsIdentity()) { return; }
+   if (nfdofs < 2) { return; }
 
    MFEM_VERIFY(Fo.Size() >= nfaces,
                "Face orientation array is shorter than the number of faces in "
-               "ND_DofTransformation");
+               "ND_StatelessDofTransformation");
 
    double data[2];
    Vector v2(data, 2);
-   DenseMatrix T2Inv;
 
    // Transform face DoFs
    for (int f=0; f<nfaces; f++)
@@ -257,24 +264,23 @@ void ND_DofTransformation::InvTransformPrimal(const Array<int> & Fo,
       for (int i=0; i<nfdofs/2; i++)
       {
          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]);
+         TInv(Fo[f]).Mult(v2, &v[nedges*nedofs + f*nfdofs + 2*i]);
       }
    }
 }
 
-void ND_DofTransformation::TransformDual(const Array<int> & Fo, double *v) const
+void ND_StatelessDofTransformation::TransformDual(const Array<int> & Fo,
+                                                  double *v) const
 {
    // Return immediately when no face DoFs are present
-   if (IsIdentity()) { return; }
+   if (nfdofs < 2) { return; }
 
    MFEM_VERIFY(Fo.Size() >= nfaces,
                "Face orientation array is shorter than the number of faces in "
-               "ND_DofTransformation");
+               "ND_StatelessDofTransformation");
 
    double data[2];
    Vector v2(data, 2);
-   DenseMatrix T2Inv;
 
    // Transform face DoFs
    for (int f=0; f<nfaces; f++)
@@ -282,25 +288,23 @@ void ND_DofTransformation::TransformDual(const Array<int> & Fo, double *v) const
       for (int i=0; i<nfdofs/2; i++)
       {
          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]);
+         TInv(Fo[f]).MultTranspose(v2, &v[nedges*nedofs + f*nfdofs + 2*i]);
       }
    }
 }
 
-void ND_DofTransformation::InvTransformDual(const Array<int> & Fo,
-                                            double *v) const
+void ND_StatelessDofTransformation::InvTransformDual(const Array<int> & Fo,
+                                                     double *v) const
 {
    // Return immediately when no face DoFs are present
-   if (IsIdentity()) { return; }
+   if (nfdofs < 2) { return; }
 
    MFEM_VERIFY(Fo.Size() >= nfaces,
                "Face orientation array is shorter than the number of faces in "
-               "ND_DofTransformation");
+               "ND_StatelessDofTransformation");
 
    double data[2];
    Vector v2(data, 2);
-   DenseMatrix T2;
 
    // Transform face DoFs
    for (int f=0; f<nfaces; f++)
@@ -308,8 +312,7 @@ void ND_DofTransformation::InvTransformDual(const Array<int> & Fo,
       for (int i=0; i<nfdofs/2; i++)
       {
          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]);
+         T(Fo[f]).MultTranspose(v2, &v[nedges*nedofs + f*nfdofs + 2*i]);
       }
    }
 }

fem/doftrans.hpp

@@ -80,9 +80,6 @@ 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
@@ -118,8 +115,6 @@ 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
@@ -138,76 +133,35 @@ 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
+class DofTransformation : virtual public StatelessDofTransformation
 {
 protected:
-   Array<int> Fo_;
-   const StatelessDofTransformation * dof_trans_;
-   int vdim_;
-   int ordering_;
+   Array<int> Fo;
+
+   DofTransformation(int size)
+      : StatelessDofTransformation(size) {}
 
 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> & Fo)
-   { Fo_ = Fo; }
+   inline void SetFaceOrientations(const Array<int> & face_orientation)
+   { Fo = face_orientation; }
 
-   /// Return the face orientations for the current element
-   inline const Array<int> & GetFaceOrientations() const { return Fo_; }
+   inline const Array<int> & GetFaceOrientations() const { return Fo; }
 
-   /// 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(); }
+   using StatelessDofTransformation::TransformPrimal;
+   using StatelessDofTransformation::InvTransformPrimal;
+   using StatelessDofTransformation::TransformDual;
+   using StatelessDofTransformation::InvTransformDual;
 
    /** 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. */
-   void TransformPrimal(double *v) const;
+   inline void TransformPrimal(double *v) const
+   { TransformPrimal(Fo, v); }
    inline void TransformPrimal(Vector &v) const
    { TransformPrimal(v.GetData()); }
 
@@ -225,18 +179,21 @@ public:
        transform the vector obtained using GridFunction::GetSubVector before it
        can be used to compute a local interpolation.
    */
-   void InvTransformPrimal(double *v) const;
+   inline void InvTransformPrimal(double *v) const
+   { InvTransformPrimal(Fo, v); }
    inline void InvTransformPrimal(Vector &v) const
    { InvTransformPrimal(v.GetData()); }
 
    /** Transform dual DoFs as computed by a LinearFormIntegrator before summing
        into a LinearForm object. */
-   void TransformDual(double *v) const;
+   inline void TransformDual(double *v) const
+   { TransformDual(Fo, v); }
    inline void TransformDual(Vector &v) const
    { TransformDual(v.GetData()); }
 
    /** Inverse Transform dual DoFs */
-   void InvTransformDual(double *v) const;
+   inline void InvTransformDual(double *v) const
+   { InvTransformDual(Fo, v); }
    inline void InvTransformDual(Vector &v) const
    { InvTransformDual(v.GetData()); }
 
@@ -268,6 +225,8 @@ public:
          TransformDual(V.GetColumn(c));
       }
    }
+
+   virtual ~DofTransformation() = default;
 };
 
 /** Transform a matrix of DoFs entries from different finite element spaces as
@@ -286,6 +245,145 @@ 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.
 
@@ -298,7 +396,7 @@ void TransformDual(const DofTransformation *ran_dof_trans,
     be accessed as DenseMatrices using the GetFaceTransform() and
     GetFaceInverseTransform() methods.
 */
-class ND_DofTransformation : public StatelessDofTransformation
+class ND_StatelessDofTransformation : virtual public StatelessDofTransformation
 {
 private:
    static const double T_data[24];
@@ -312,7 +410,8 @@ protected:
    const int nedges; // number of edges per element
    const int nfaces; // number of triangular faces per element
 
-   ND_DofTransformation(int size, int order, int num_edges, int num_tri_faces);
+   ND_StatelessDofTransformation(int size, int order,
+                                 int num_edges, int num_tri_faces);
 
 public:
    // Return the 2x2 transformation operator for the given face orientation
@@ -322,41 +421,116 @@ public:
    static const DenseMatrix & GetFaceInverseTransform(int ori)
    { return TInv(ori); }
 
-   bool IsIdentity() const override { return nfdofs < 2; }
+   void TransformPrimal(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;
+   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;
 };
 
 /// Stateless DoF transformation implementation for the Nedelec basis on
 /// triangles
-class ND_TriDofTransformation : public ND_DofTransformation
+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
 {
 public:
    ND_TriDofTransformation(int order)
-      : ND_DofTransformation(order*(order + 2), order, 3, 1)
+      : 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)
    {}
 };
 
 /// DoF transformation implementation for the Nedelec basis on tetrahedra
-class ND_TetDofTransformation : public ND_DofTransformation
+class ND_TetDofTransformation : public DofTransformation,
+   public ND_TetStatelessDofTransformation
 {
 public:
    ND_TetDofTransformation(int order)
-      : ND_DofTransformation(order*(order + 2)*(order + 3)/2, order, 6, 4)
+      : 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)
    {}
 };
 
 /// DoF transformation implementation for the Nedelec basis on wedge elements
-class ND_WedgeDofTransformation : public ND_DofTransformation
+class ND_WedgeDofTransformation : public DofTransformation,
+   public ND_WedgeStatelessDofTransformation
 {
 public:
    ND_WedgeDofTransformation(int order)
-      : ND_DofTransformation(3 * order * ((order + 1) * (order + 2))/2,
-                             order, 9, 2)
+      : StatelessDofTransformation(3 * order * ((order + 1) * (order + 2))/2)
+      , DofTransformation(3 * order * ((order + 1) * (order + 2))/2)
+      , ND_WedgeStatelessDofTransformation(order)
    {}
+
+   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 const StatelessDofTransformation *GetDofTransformation() const
+   virtual 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;
 
-   ND_TetDofTransformation doftrans;
+   mutable ND_TetStatelessDofTransformation 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 const StatelessDofTransformation *GetDofTransformation() const
+   virtual 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;
 
-   ND_TriDofTransformation doftrans;
+   mutable ND_TriStatelessDofTransformation 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 const StatelessDofTransformation *GetDofTransformation() const
+   virtual 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;
 
-   ND_WedgeDofTransformation doftrans;
+   mutable ND_WedgeStatelessDofTransformation doftrans;
 
    H1_TriangleElement H1TriangleFE;
    ND_TriangleElement NDTriangleFE;
@@ -379,7 +379,7 @@ public:
                                   DenseMatrix &I) const
    { LocalInterpolation_ND(CheckVectorFE(fe), tk, dof2tk, Trans, I); }
 
-   virtual const StatelessDofTransformation *GetDofTransformation() const
+   virtual StatelessDofTransformation * GetDofTransformation() const
    { return &doftrans; }
 
    using FiniteElement::Project;

fem/fe_coll.cpp

@@ -2896,7 +2896,7 @@ ND_FECollection::FiniteElementForGeometry(Geometry::Type GeomType) const
    }
 }
 
-const StatelessDofTransformation *
+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 const StatelessDofTransformation *
+   virtual StatelessDofTransformation *
    DofTransformationForGeometry(Geometry::Type GeomType) const
    { return NULL; }
 
@@ -483,7 +483,7 @@ public:
    int DofForGeometry(Geometry::Type GeomType) const override
    { return ND_dof[GeomType]; }
 
-   const StatelessDofTransformation *
+   StatelessDofTransformation *
    DofTransformationForGeometry(Geometry::Type GeomType) const override;
 
    const int *DofOrderForOrientation(Geometry::Type GeomType,

fem/fespace.cpp

@@ -63,6 +63,7 @@ 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)
@@ -71,6 +72,7 @@ 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;
@@ -210,7 +212,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; }
@@ -262,7 +264,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++)
@@ -275,36 +277,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
 {
-   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);
+   DofTransformation * doftrans = GetElementDofs(i, vdofs);
    DofsToVDofs(vdofs);
-   doftrans.SetVDim(vdim, ordering);
+   if (vdim == 1 || doftrans == NULL)
+   {
+      return doftrans;
+   }
+   else
+   {
+      VDoFTrans.SetDofTransformation(*doftrans);
+      return &VDoFTrans;
+   }
 }
 
 DofTransformation *
 FiniteElementSpace::GetBdrElementVDofs(int i, Array<int> &vdofs) const
 {
-   DoFTrans.SetDofTransformation(NULL);
-   GetBdrElementVDofs(i, vdofs, DoFTrans);
-   return DoFTrans.GetDofTransformation() ? &DoFTrans : NULL;
+   DofTransformation * doftrans = GetBdrElementDofs(i, vdofs);
+   DofsToVDofs(vdofs);
+   if (vdim == 1 || doftrans == NULL)
+   {
+      return doftrans;
+   }
+   else
+   {
+      VDoFTrans.SetDofTransformation(*doftrans);
+      return &VDoFTrans;
+   }
 }
 
 void FiniteElementSpace::GetPatchVDofs(int i, Array<int> &vdofs) const
@@ -775,9 +777,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++)
    {
@@ -800,9 +802,11 @@ void FiniteElementSpace::AddDependencies(
    }
 }
 
-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
@@ -1529,12 +1533,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.");
@@ -1549,7 +1553,7 @@ FiniteElementSpace::RefinementOperator::RefinementOperator(
       fespace->GetLocalRefinementMatrices(elem_geoms[i], localP[elem_geoms[i]]);
    }
 
-   ConstructDoFTransArray();
+   ConstructDoFTrans();
 }
 
 FiniteElementSpace::RefinementOperator::RefinementOperator(
@@ -1574,58 +1578,59 @@ FiniteElementSpace::RefinementOperator::RefinementOperator(
       old_elem_fos = new Table(*coarse_fes->GetElementToFaceOrientationTable());
    }
 
-   ConstructDoFTransArray();
+   ConstructDoFTrans();
 }
 
 FiniteElementSpace::RefinementOperator::~RefinementOperator()
 {
    delete old_elem_dof;
    delete old_elem_fos;
-   for (int i=0; i<old_DoFTransArray.Size(); i++)
+   for (int i=0; i<old_DoFTrans.Size(); i++)
    {
-      delete old_DoFTransArray[i];
+      delete old_DoFTrans[i];
    }
 }
 
-void FiniteElementSpace::RefinementOperator::ConstructDoFTransArray()
+void FiniteElementSpace::RefinementOperator
+::ConstructDoFTrans()
 {
-   old_DoFTransArray.SetSize(Geometry::NUM_GEOMETRIES);
-   for (int i=0; i<old_DoFTransArray.Size(); i++)
+   old_DoFTrans.SetSize(Geometry::NUM_GEOMETRIES);
+   for (int i=0; i<old_DoFTrans.Size(); i++)
    {
-      old_DoFTransArray[i] = NULL;
+      old_DoFTrans[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_DoFTransArray[Geometry::TRIANGLE] =
+         old_DoFTrans[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_DoFTransArray[Geometry::TETRAHEDRON] =
+         old_DoFTrans[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_DoFTransArray[Geometry::PRISM] =
+         old_DoFTrans[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 =
@@ -1657,7 +1662,6 @@ void FiniteElementSpace::RefinementOperator::Mult(const Vector &x,
             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);
@@ -1666,30 +1670,40 @@ void FiniteElementSpace::RefinementOperator::Mult(const Vector &x,
       else
       {
          old_elem_fos->GetRow(emb.parent, old_Fo);
-         old_DoFTrans.SetDofTransformation(*old_DoFTransArray[geom]);
-         old_DoFTrans.SetFaceOrientations(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();
+         }
 
-         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.InvTransformPrimal(subX);
+            old_DoFTrans[geom]->InvTransformPrimal(subX);
             lP.Mult(subX, subY);
             doftrans->TransformPrimal(subY);
             y.SetSubVector(vdofs, subY);
          }
-         doftrans->SetVDim(rvdim, fespace->GetOrdering());
+
+         if (vdoftrans)
+         {
+            doftrans = new_doftrans;
+         }
       }
    }
 }
 
-void FiniteElementSpace::RefinementOperator::MultTranspose(const Vector &x,
-                                                           Vector &y) const
+void FiniteElementSpace::RefinementOperator
+::MultTranspose(const Vector &x, Vector &y) const
 {
    y = 0.0;
 
@@ -1713,7 +1727,7 @@ void FiniteElementSpace::RefinementOperator::MultTranspose(const Vector &x,
       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)
@@ -1728,6 +1742,7 @@ void FiniteElementSpace::RefinementOperator::MultTranspose(const Vector &x,
             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])])
@@ -1735,6 +1750,7 @@ void FiniteElementSpace::RefinementOperator::MultTranspose(const Vector &x,
                   subX[p] = 0.0;
                }
             }
+
             lP.MultTranspose(subX, subY);
             y.AddElementVector(c_vdofs, subY);
          }
@@ -1744,10 +1760,17 @@ void FiniteElementSpace::RefinementOperator::MultTranspose(const Vector &x,
          subYt.SetSize(lP.Width());
 
          old_elem_fos->GetRow(emb.parent, old_Fo);
-         old_DoFTrans.SetDofTransformation(*old_DoFTransArray[geom]);
-         old_DoFTrans.SetFaceOrientations(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();
+         }
 
-         doftrans->SetVDim();
          for (int vd = 0; vd < rvdim; vd++)
          {
             f_dofs.Copy(f_vdofs);
@@ -1764,11 +1787,16 @@ void FiniteElementSpace::RefinementOperator::MultTranspose(const Vector &x,
                   subX[p] = 0.0;
                }
             }
+
             lP.MultTranspose(subX, subYt);
-            old_DoFTrans.TransformDual(subYt);
+            old_DoFTrans[geom]->TransformDual(subYt);
             y.AddElementVector(c_vdofs, subYt);
          }
-         doftrans->SetVDim(rvdim, fespace->GetOrdering());
+
+         if (vdoftrans)
+         {
+            doftrans = new_doftrans;
+         }
       }
 
       for (int p = 0; p < f_dofs.Size(); ++p)
@@ -1996,8 +2024,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;
@@ -2199,7 +2227,7 @@ void FiniteElementSpace::Constructor(Mesh *mesh_, NURBSExtension *NURBSext_,
       R_transpose.reset();
       cP_is_set = false;
 
-      ConstructDoFTransArray();
+      ConstructDoFTrans();
    }
    else
    {
@@ -2211,39 +2239,40 @@ void FiniteElementSpace::Constructor(Mesh *mesh_, NURBSExtension *NURBSext_,
    BuildElementToDofTable();
 }
 
-void FiniteElementSpace::ConstructDoFTransArray()
+void FiniteElementSpace::ConstructDoFTrans()
 {
-   DestroyDoFTransArray();
+   DestroyDoFTrans();
 
-   DoFTransArray.SetSize(Geometry::NUM_GEOMETRIES);
-   for (int i=0; i<DoFTransArray.Size(); i++)
+   VDoFTrans.SetVDim(vdim);
+   DoFTrans.SetSize(Geometry::NUM_GEOMETRIES);
+   for (int i=0; i<DoFTrans.Size(); i++)
    {
-      DoFTransArray[i] = NULL;
+      DoFTrans[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)
       {
-         DoFTransArray[Geometry::TRIANGLE] =
+         DoFTrans[Geometry::TRIANGLE] =
             new ND_TriDofTransformation(nd_tri->GetOrder());
       }
 
-      const FiniteElement *nd_tet =
+      const FiniteElement * nd_tet =
          fec->FiniteElementForGeometry(Geometry::TETRAHEDRON);
       if (nd_tet)
       {
-         DoFTransArray[Geometry::TETRAHEDRON] =
+         DoFTrans[Geometry::TETRAHEDRON] =
             new ND_TetDofTransformation(nd_tet->GetOrder());
       }
 
-      const FiniteElement *nd_pri =
+      const FiniteElement * nd_pri =
          fec->FiniteElementForGeometry(Geometry::PRISM);
       if (nd_pri)
       {
-         DoFTransArray[Geometry::PRISM] =
+         DoFTrans[Geometry::PRISM] =
             new ND_WedgeDofTransformation(nd_pri->GetOrder());
       }
    }
@@ -2295,7 +2324,7 @@ void FiniteElementSpace::BuildNURBSFaceToDofTable() const
    face_to_be = -1;
    for (int b = 0; b < GetNBE(); b++)
    {
-      int f = mesh->GetBdrElementFaceIndex(b);
+      int f = mesh->GetBdrElementEdgeIndex(b);
       face_to_be[f] = b;
    }
 
@@ -2447,7 +2476,7 @@ void FiniteElementSpace::Construct()
 
    ndofs = nvdofs + nedofs + nfdofs + nbdofs;
 
-   ConstructDoFTransArray();
+   ConstructDoFTrans();
 
    // record the current mesh sequence number to detect refinement etc.
    mesh_sequence = mesh->GetSequence();
@@ -2472,8 +2501,9 @@ 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(), "");
@@ -2697,8 +2727,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.";
 
-void FiniteElementSpace::GetElementDofs(int elem, Array<int> &dofs,
-                                        DofTransformation &doftrans) const
+DofTransformation *
+FiniteElementSpace::GetElementDofs(int elem, Array<int> &dofs) const
 {
    MFEM_VERIFY(!orders_changed, msg_orders_changed);
 
@@ -2706,16 +2736,13 @@ void FiniteElementSpace::GetElementDofs(int elem, Array<int> &dofs,
    {
       elem_dof->GetRow(elem, dofs);
 
-      if (DoFTransArray[mesh->GetElementBaseGeometry(elem)])
+      if (DoFTrans[mesh->GetElementBaseGeometry(elem)])
       {
          Array<int> Fo;
          elem_fos -> GetRow (elem, Fo);
-         doftrans.SetDofTransformation(
-            *DoFTransArray[mesh->GetElementBaseGeometry(elem)]);
-         doftrans.SetFaceOrientations(Fo);
-         doftrans.SetVDim();
+         DoFTrans[mesh->GetElementBaseGeometry(elem)]->SetFaceOrientations(Fo);
       }
-      return;
+      return DoFTrans[mesh->GetElementBaseGeometry(elem)];
    }
 
    Array<int> V, E, Eo, F, Fo; // TODO: LocalArray
@@ -2739,12 +2766,10 @@ void FiniteElementSpace::GetElementDofs(int elem, Array<int> &dofs,
       {
          nfd += fec->GetNumDof(mesh->GetFaceGeometry(F[i]), order);
       }
-      if (DoFTransArray[mesh->GetElementBaseGeometry(elem)])
+      if (DoFTrans[mesh->GetElementBaseGeometry(elem)])
       {
-         doftrans.SetDofTransformation(
-            *DoFTransArray[mesh->GetElementBaseGeometry(elem)]);
-         doftrans.SetFaceOrientations(Fo);
-         doftrans.SetVDim();
+         DoFTrans[mesh->GetElementBaseGeometry(elem)]
+         -> SetFaceOrientations(Fo);
       }
    }
 
@@ -2803,18 +2828,54 @@ void FiniteElementSpace::GetElementDofs(int elem, Array<int> &dofs,
          dofs.Append(bbase + j);
       }
    }
+   return DoFTrans[mesh->GetElementBaseGeometry(elem)];
 }
 
-DofTransformation *FiniteElementSpace::GetElementDofs(int elem,
-                                                      Array<int> &dofs) const
+void FiniteElementSpace::GetPatchDofs(int patch, Array<int> &dofs) const
 {
-   DoFTrans.SetDofTransformation(NULL);
-   GetElementDofs(elem, dofs, DoFTrans);
-   return DoFTrans.GetDofTransformation() ? &DoFTrans : NULL;
+   MFEM_ASSERT(NURBSext,
+               "FiniteElementSpace::GetPatchDofs needs a NURBSExtension");
+   NURBSext->GetPatchDofs(patch, dofs);
 }
 
-void FiniteElementSpace::GetBdrElementDofs(int bel, Array<int> &dofs,
-                                           DofTransformation &doftrans) const
+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
 {
    MFEM_VERIFY(!orders_changed, msg_orders_changed);
 
@@ -2822,19 +2883,17 @@ void FiniteElementSpace::GetBdrElementDofs(int bel, Array<int> &dofs,
    {
       bdr_elem_dof->GetRow(bel, dofs);
 
-      if (DoFTransArray[mesh->GetBdrElementBaseGeometry(bel)])
+      if (DoFTrans[mesh->GetBdrElementBaseGeometry(bel)])
       {
          Array<int> Fo;
          bdr_elem_fos -> GetRow (bel, Fo);
-         doftrans.SetDofTransformation(
-            *DoFTransArray[mesh->GetBdrElementBaseGeometry(bel)]);
-         doftrans.SetFaceOrientations(Fo);
-         doftrans.SetVDim();
+         DoFTrans[mesh->GetBdrElementBaseGeometry(bel)]->
+         SetFaceOrientations(Fo);
       }
-      return;
+      return DoFTrans[mesh->GetBdrElementBaseGeometry(bel)];
    }
 
-   Array<int> V, E, Eo; // TODO: LocalArray
+   Array<int> V, E, Eo, Fo; // TODO: LocalArray
    int F, oF;
 
    int dim = mesh->Dimension();
@@ -2858,14 +2917,11 @@ void FiniteElementSpace::GetBdrElementDofs(int bel, Array<int> &dofs,
    {
       mesh->GetBdrElementFace(bel, &F, &oF);
 
-      if (DoFTransArray[mesh->GetBdrElementBaseGeometry(bel)])
+      if (DoFTrans[mesh->GetBdrElementBaseGeometry(bel)])
       {
-         mfem::Array<int> Fo(1);
-         Fo[0] = oF;
-         doftrans.SetDofTransformation(
-            *DoFTransArray[mesh->GetBdrElementBaseGeometry(bel)]);
-         doftrans.SetFaceOrientations(Fo);
-         doftrans.SetVDim();
+         Fo.Append(oF);
+         DoFTrans[mesh->GetBdrElementBaseGeometry(bel)]->
+         SetFaceOrientations(Fo);
       }
    }
 
@@ -2907,14 +2963,8 @@ void FiniteElementSpace::GetBdrElementDofs(int bel, Array<int> &dofs,
          dofs.Append(EncodeDof(nvdofs + nedofs + fbase, ind[j]));
       }
    }
-}
 
-DofTransformation *FiniteElementSpace::GetBdrElementDofs(int bel,
-                                                         Array<int> &dofs) const
-{
-   DoFTrans.SetDofTransformation(NULL);
-   GetBdrElementDofs(bel, dofs, DoFTrans);
-   return DoFTrans.GetDofTransformation() ? &DoFTrans : NULL;
+   return DoFTrans[mesh->GetBdrElementBaseGeometry(bel)];
 }
 
 int FiniteElementSpace::GetFaceDofs(int face, Array<int> &dofs,
@@ -3084,6 +3134,18 @@ 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");
@@ -3108,61 +3170,6 @@ 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();
@@ -3235,8 +3242,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);
 }
@@ -3276,7 +3283,7 @@ void FiniteElementSpace::Destroy()
    }
    E2BFQ_array.SetSize(0);
 
-   DestroyDoFTransArray();
+   DestroyDoFTrans();
 
    dof_elem_array.DeleteAll();
    dof_ldof_array.DeleteAll();
@@ -3294,18 +3301,19 @@ void FiniteElementSpace::Destroy()
       delete bdr_elem_dof;
       delete bdr_elem_fos;
       delete face_dof;
+
       delete [] bdofs;
    }
    ceed::RemoveBasisAndRestriction(this);
 }
 
-void FiniteElementSpace::DestroyDoFTransArray()
+void FiniteElementSpace::DestroyDoFTrans()
 {
-   for (int i = 0; i < DoFTransArray.Size(); i++)
+   for (int i = 0; i < DoFTrans.Size(); i++)
    {
-      delete DoFTransArray[i];
+      delete DoFTrans[i];
    }
-   DoFTransArray.SetSize(0);
+   DoFTrans.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<StatelessDofTransformation *> DoFTransArray;
-   mutable DofTransformation DoFTrans;
+   Array<DofTransformation*> DoFTrans;
+   mutable VDofTransformation VDoFTrans;
 
    /** 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 ConstructDoFTransArray();
-   void DestroyDoFTransArray();
+   void ConstructDoFTrans();
+   void DestroyDoFTrans();
 
    void BuildElementToDofTable() const;
    void BuildBdrElementToDofTable() const;
@@ -416,10 +416,10 @@ protected:
       Table* old_elem_dof; // Owned.
       Table* old_elem_fos; // Owned.
 
-      Array<StatelessDofTransformation*> old_DoFTransArray;
-      mutable DofTransformation old_DoFTrans;
+      Array<DofTransformation*> old_DoFTrans;
+      mutable VDofTransformation old_VDoFTrans;
 
-      void ConstructDoFTransArray();
+      void ConstructDoFTrans();
 
    public:
       /** Construct the operator based on the elem_dof table of the original
@@ -803,16 +803,7 @@ public:
    /// with triangular faces.
    ///
    /// @note The returned object should NOT be deleted by the caller.
-   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;
+   virtual DofTransformation *GetElementDofs(int elem, Array<int> &dofs) const;
 
    /// @brief Returns indices of degrees of freedom for boundary element 'bel'.
    /// The returned indices are offsets into an @ref ldof vector. See also
@@ -826,16 +817,13 @@ public:
    /// with triangular faces.
    ///
    /// @note The returned object should NOT be deleted by the caller.
-   DofTransformation *GetBdrElementDofs(int bel, Array<int> &dofs) const;
+   virtual DofTransformation *GetBdrElementDofs(int bel,
+                                                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 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 Returns the indices of the degrees of freedom for the specified
    /// face, including the DOFs for the edges and the vertices of the face.
@@ -882,13 +870,6 @@ 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.
    ///
@@ -901,6 +882,13 @@ 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.
    ///
@@ -909,11 +897,6 @@ 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
@@ -1040,15 +1023,6 @@ 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
@@ -1064,15 +1038,6 @@ 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->GetBdrElementFaceIndex(i));
+            bdr_eles[FMS_VERTEX].push_back(mmesh->GetBdrElementEdgeIndex(i));
             break;
          case Element::SEGMENT:
-            bdr_eles[FMS_EDGE].push_back(mmesh->GetBdrElementFaceIndex(i));
+            bdr_eles[FMS_EDGE].push_back(mmesh->GetBdrElementEdgeIndex(i));
             break;
          case Element::TRIANGLE:
-            bdr_eles[FMS_TRIANGLE].push_back(mmesh->GetBdrElementFaceIndex(i));
+            bdr_eles[FMS_TRIANGLE].push_back(mmesh->GetBdrElementEdgeIndex(i));
             break;
          case Element::QUADRILATERAL:
-            bdr_eles[FMS_QUADRILATERAL].push_back(mmesh->GetBdrElementFaceIndex(i));
+            bdr_eles[FMS_QUADRILATERAL].push_back(mmesh->GetBdrElementEdgeIndex(i));
             break;
          case Element::TETRAHEDRON:
-            bdr_eles[FMS_TETRAHEDRON].push_back(mmesh->GetBdrElementFaceIndex(i));
+            bdr_eles[FMS_TETRAHEDRON].push_back(mmesh->GetBdrElementEdgeIndex(i));
             break;
          case Element::HEXAHEDRON:
-            bdr_eles[FMS_HEXAHEDRON].push_back(mmesh->GetBdrElementFaceIndex(i));
+            bdr_eles[FMS_HEXAHEDRON].push_back(mmesh->GetBdrElementEdgeIndex(i));
             break;
          default:
             MFEM_WARNING("Unsupported boundary element " << betype << " at boundary index "

fem/gridfunc.cpp

@@ -38,8 +38,9 @@ GridFunction::GridFunction(Mesh *m, std::istream &input)
    // Grid functions are stored on the device
    UseDevice(true);
 
-   fes = new FiniteElementSpace;
-   fec = fes->Load(m, input);
+   owned_fes.reset(new FiniteElementSpace);
+   fes = owned_fes.get();
+   fec.reset(fes->Load(m, input));
 
    skip_comment_lines(input, '#');
    istream::int_type next_char = input.peek();
@@ -81,10 +82,11 @@ GridFunction::GridFunction(Mesh *m, GridFunction *gf_array[], int num_pieces)
    int vdim, ordering;
 
    fes = gf_array[0]->FESpace();
-   fec = FiniteElementCollection::New(fes->FEColl()->Name());
+   fec.reset(FiniteElementCollection::New(fes->FEColl()->Name()));
    vdim = fes->GetVDim();
    ordering = fes->GetOrdering();
-   fes = new FiniteElementSpace(m, fec, vdim, ordering);
+   owned_fes.reset(new FiniteElementSpace(m, fec.get(), vdim, ordering));
+   fes = owned_fes.get();
    SetSize(fes->GetVSize());
 
    if (m->NURBSext)
@@ -153,12 +155,9 @@ GridFunction::GridFunction(Mesh *m, GridFunction *gf_array[], int num_pieces)
 
 void GridFunction::Destroy()
 {
-   if (fec)
-   {
-      delete fes;
-      delete fec;
-      fec = NULL;
-   }
+   owned_fes.reset();
+   fec.reset();
+   fes = nullptr;
 }
 
 void GridFunction::Update()

fem/gridfunc.hpp

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

fem/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.GetBdrElementFaceIndex(be);
+         const int f = mesh.GetBdrElementEdgeIndex(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.GetBdrElementFaceIndex(i);
+         const int f = mesh.GetBdrElementEdgeIndex(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 (sdim,
+   // In this function, we need to convert X_vert (which has the shape (dim,
    // ndof_per_el, nel_ho)) to T-DOF format.
    //
-   // We place the results in the vector xyz_tvec, which has shape (ntdofs, sdim)
+   // We place the results in the vector xyz_tvec, which has shape (ntdofs, dim)
    // and then make the hypre vectors x, y, and z point to subvectors.
    //
-   // When the space dimension is 2, z is NULL.
+   // In 2D, 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 = vert_fes.GetMesh()->SpaceDimension();
+   const int sdim = dim;
    const int ntdofs = R->Height();
 
    const MemoryClass mc = GetHypreMemoryClass();
    bool dev = (mc == MemoryClass::DEVICE);
 
-   xyz_tvec = new Vector(ntdofs*sdim);
+   xyz_tvec = new Vector(ntdofs*dim);
 
-   auto xyz_tv = Reshape(HypreWrite(xyz_tvec->GetMemory()), ntdofs, sdim);
+   auto xyz_tv = Reshape(HypreWrite(xyz_tvec->GetMemory()), ntdofs, dim);
    const auto xyz_e =
-      Reshape(HypreRead(X_vert.GetMemory()), sdim, ndof_per_el, nel_ho);
+      Reshape(HypreRead(X_vert.GetMemory()), dim, 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 (sdim == 3)
+   if (dim == 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,7 +77,6 @@ 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;
@@ -95,7 +94,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(sdim*ndof_per_el*nel_ho);
+   X_vert.SetSize(dim*ndof_per_el*nel_ho);
    const QuadratureInterpolator *quad_interp =
       nodal_fes->GetQuadratureInterpolator(ir);
    quad_interp->SetOutputLayout(QVectorLayout::byVDIM);
@@ -381,49 +380,44 @@ 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();
 
-   Assemble_(kernel, dim, sdim, order);
+   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 << "!");
+      }
+   }
 }
 
 void BatchedLORAssembly::AssembleWithoutBC(BilinearForm &a, OperatorHandle &A)

fem/lor/lor_h1.cpp

@@ -9,6 +9,7 @@
 // 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"
@@ -16,7 +17,7 @@
 namespace mfem
 {
 
-template <int ORDER, int SDIM>
+template <int ORDER>
 void BatchedLOR_H1::Assemble2D()
 {
    const int nel_ho = fes_ho.GetNE();
@@ -73,8 +74,31 @@ void BatchedLOR_H1::Assemble2D()
 
             for (int i=0; i<sz_local_mat; ++i) { local_mat[i] = 0.0; }
 
-            SetupLORQuadData2D<ORDER,SDIM,false,false>(X, iel_ho, kx, ky, Q, false);
+            double vx[4], vy[4];
+            LORVertexCoordinates2D<ORDER>(X, iel_ho, kx, ky, vx, vy);
 
+            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)
@@ -495,4 +519,34 @@ 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_h1.hpp

@@ -22,22 +22,15 @@ namespace mfem
 class BatchedLOR_H1 : BatchedLORKernel
 {
 public:
-   template <int ORDER, int SDIM> void Assemble2D();
+   template <int ORDER> void Assemble2D();
    template <int ORDER> void Assemble3D();
    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);
-   }
+                 Array<int> &sparse_mapping_);
 };
 
 }
 
-#include "lor_h1_impl.hpp"
-
 #endif

fem/lor/lor_nd.cpp

@@ -9,6 +9,7 @@
 // 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"
@@ -16,7 +17,7 @@
 namespace mfem
 {
 
-template <int ORDER, int SDIM>
+template <int ORDER>
 void BatchedLOR_ND::Assemble2D()
 {
    const int nel_ho = fes_ho.GetNE();
@@ -82,8 +83,31 @@ 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; }
 
-            SetupLORQuadData2D<ORDER,SDIM,false,true>(X, iel_ho, kx, ky, Q, true);
+            double vx[4], vy[4];
+            LORVertexCoordinates2D<ORDER>(X, iel_ho, kx, ky, vx, vy);
 
+            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)
@@ -539,4 +563,34 @@ 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_nd.hpp

@@ -22,22 +22,15 @@ namespace mfem
 class BatchedLOR_ND : BatchedLORKernel
 {
 public:
-   template <int ORDER, int SDIM> void Assemble2D();
+   template <int ORDER> void Assemble2D();
    template <int ORDER> void Assemble3D();
    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);
-   }
+                 Array<int> &sparse_mapping_);
 };
 
 }
 
-#include "lor_nd_impl.hpp"
-
 #endif

fem/lor/lor_rt.cpp

@@ -9,6 +9,7 @@
 // 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"
@@ -16,7 +17,7 @@
 namespace mfem
 {
 
-template <int ORDER, int SDIM>
+template <int ORDER>
 void BatchedLOR_RT::Assemble2D()
 {
    const int nel_ho = fes_ho.GetNE();
@@ -78,8 +79,31 @@ 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; }
 
-            SetupLORQuadData2D<ORDER,SDIM,true,false>(X, iel_ho, kx, ky, Q, true);
+            double vx[4], vy[4];
+            LORVertexCoordinates2D<ORDER>(X, iel_ho, kx, ky, vx, vy);
 
+            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)
@@ -523,4 +547,34 @@ 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_rt.hpp

@@ -22,22 +22,15 @@ namespace mfem
 class BatchedLOR_RT : BatchedLORKernel
 {
 public:
-   template <int ORDER, int SDIM> void Assemble2D();
+   template <int ORDER> void Assemble2D();
    template <int ORDER> void Assemble3D();
    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);
-   }
+                 Array<int> &sparse_mapping_);
 };
 
 }
 
-#include "lor_rt_impl.hpp"
-
 #endif

fem/lor/lor_util.hpp

@@ -20,22 +20,11 @@
 namespace mfem
 {
 
-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>
+template <int ORDER>
 MFEM_HOST_DEVICE inline void LORVertexCoordinates2D(
-   const double *X, int iel_ho, int kx, int ky, double **v)
+   const double *X, int iel_ho, int kx, int ky, double vx[4], double vy[4])
 {
+   const int dim = 2;
    const int nd1d = ORDER + 1;
    const int nvert_per_el = nd1d*nd1d;
 
@@ -44,31 +33,23 @@ MFEM_HOST_DEVICE inline void LORVertexCoordinates2D(
    const int v2 = kx + 1 + nd1d*(ky + 1);
    const int v3 = kx + nd1d*(ky + 1);
 
-   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);
+   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);
 
    // Vertex coordinates
-   v[0][0] = X[e0 + 0];
-   v[1][0] = X[e0 + 1];
+   vx[0] = X[e0 + 0];
+   vy[0] = X[e0 + 1];
 
-   v[0][1] = X[e1 + 0];
-   v[1][1] = X[e1 + 1];
+   vx[1] = X[e1 + 0];
+   vy[1] = X[e1 + 1];
 
-   v[0][2] = X[e2 + 0];
-   v[1][2] = X[e2 + 1];
+   vx[2] = X[e2 + 0];
+   vy[2] = X[e2 + 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];
-   }
+   vx[3] = X[e3 + 0];
+   vy[3] = X[e3 + 1];
 }
 
 template <int ORDER>
@@ -131,80 +112,15 @@ 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, double **v, DeviceMatrix &J);
-
-template <> MFEM_HOST_DEVICE inline void Jacobian2D<2>(
-   const double x, const double y, double **v, DeviceMatrix &J)
+   const double x, const double y, const double vx[4], const double vy[4],
+   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(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];
 
-   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;
-         }
-      }
-   }
+   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];
 }
 
 MFEM_HOST_DEVICE inline void Jacobian3D(
@@ -264,6 +180,18 @@ 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,37 +97,12 @@ 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);
@@ -135,9 +110,6 @@ 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);
          }
       }
@@ -203,32 +175,8 @@ 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);
@@ -236,9 +184,6 @@ 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);
@@ -377,32 +322,8 @@ 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);
@@ -410,9 +331,6 @@ 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);
@@ -643,6 +561,13 @@ 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)
 {
@@ -667,7 +592,6 @@ 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)
@@ -677,33 +601,8 @@ 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)
          {
@@ -715,13 +614,9 @@ 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)
@@ -761,7 +656,6 @@ 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;
@@ -776,32 +670,8 @@ 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)
          {
@@ -813,9 +683,6 @@ 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);
 
@@ -831,6 +698,7 @@ void BlockNonlinearForm::MultBlocked(const BlockVector &bx,
 
    if (fnfi.Size())
    {
+      Mesh *mesh = fes[0]->GetMesh();
       FaceElementTransformations *tr;
 
       for (int i = 0; i < mesh->GetNumFaces(); ++i)
@@ -868,8 +736,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);
@@ -990,7 +858,6 @@ 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)
    {
@@ -1021,32 +888,8 @@ 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)
          {
@@ -1058,9 +901,6 @@ 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)
@@ -1083,6 +923,7 @@ 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)
       {
@@ -1119,6 +960,7 @@ 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,7 +37,6 @@ 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
@@ -109,12 +108,7 @@ public:
 
    /// Adds new Domain Integrator.
    void AddDomainIntegrator(NonlinearFormIntegrator *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); }
+   { dnfi.Append(nlfi); }
 
    /// Access all integrators added with AddDomainIntegrator().
    Array<NonlinearFormIntegrator*> *GetDNFI() { return &dnfi; }
@@ -233,14 +227,13 @@ 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). */
@@ -305,12 +298,7 @@ public:
 
    /// Adds new Domain Integrator.
    void AddDomainIntegrator(BlockNonlinearFormIntegrator *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); }
+   { dnfi.Append(nlfi); }
 
    /// Adds new Interior Face Integrator.
    void AddInteriorFaceIntegrator(BlockNonlinearFormIntegrator *nlfi)
@@ -323,8 +311,7 @@ public:
    /** @brief Adds new Boundary Face Integrator, restricted to specific boundary
        attributes. */
    void AddBdrFaceIntegrator(BlockNonlinearFormIntegrator *nlfi,
-                             Array<int> &bdr_marker)
-   { bfnfi.Append(nlfi); bfnfi_marker.Append(&bdr_marker); }
+                             Array<int> &bdr_marker);
 
    virtual void SetEssentialBC(const Array<Array<int> *>&bdr_attr_is_ess,
                                Array<Vector *> &rhs);

fem/pfespace.cpp

@@ -466,54 +466,53 @@ void ParFiniteElementSpace::ApplyLDofSigns(Table &el_dof) const
    ApplyLDofSigns(all_dofs);
 }
 
-void ParFiniteElementSpace::GetElementDofs(int i, Array<int> &dofs,
-                                           DofTransformation &doftrans) const
+DofTransformation *
+ParFiniteElementSpace::GetElementDofs(int i, Array<int> &dofs) const
 {
    if (elem_dof)
    {
       elem_dof->GetRow(i, dofs);
 
-      if (DoFTransArray[mesh->GetElementBaseGeometry(i)])
+      if (DoFTrans[mesh->GetElementBaseGeometry(i)])
       {
          Array<int> Fo;
          elem_fos->GetRow(i, Fo);
-         doftrans.SetDofTransformation(
-            *DoFTransArray[mesh->GetElementBaseGeometry(i)]);
-         doftrans.SetFaceOrientations(Fo);
-         doftrans.SetVDim();
+         DoFTrans[mesh->GetElementBaseGeometry(i)]->SetFaceOrientations(Fo);
+         return DoFTrans[mesh->GetElementBaseGeometry(i)];
       }
-      return;
+      return NULL;
    }
-   FiniteElementSpace::GetElementDofs(i, dofs, doftrans);
+   DofTransformation * doftrans = FiniteElementSpace::GetElementDofs(i, dofs);
    if (Conforming())
    {
       ApplyLDofSigns(dofs);
    }
+   return doftrans;
 }
 
-void ParFiniteElementSpace::GetBdrElementDofs(int i, Array<int> &dofs,
-                                              DofTransformation &doftrans) const
+DofTransformation *
+ParFiniteElementSpace::GetBdrElementDofs(int i, Array<int> &dofs) const
 {
    if (bdr_elem_dof)
    {
       bdr_elem_dof->GetRow(i, dofs);
 
-      if (DoFTransArray[mesh->GetBdrElementBaseGeometry(i)])
+      if (DoFTrans[mesh->GetBdrElementBaseGeometry(i)])
       {
          Array<int> Fo;
-         bdr_elem_fos->GetRow(i, Fo);
-         doftrans.SetDofTransformation(
-            *DoFTransArray[mesh->GetBdrElementBaseGeometry(i)]);
-         doftrans.SetFaceOrientations(Fo);
-         doftrans.SetVDim();
+         bdr_elem_fos -> GetRow (i, Fo);
+         DoFTrans[mesh->GetBdrElementBaseGeometry(i)]->SetFaceOrientations(Fo);
+         return DoFTrans[mesh->GetBdrElementBaseGeometry(i)];
       }
-      return;
+      return NULL;
    }
-   FiniteElementSpace::GetBdrElementDofs(i, dofs, doftrans);
+   DofTransformation * doftrans =
+      FiniteElementSpace::GetBdrElementDofs(i, dofs);
    if (Conforming())
    {
       ApplyLDofSigns(dofs);
    }
+   return doftrans;
 }
 
 int ParFiniteElementSpace::GetFaceDofs(int i, Array<int> &dofs,
@@ -940,8 +939,8 @@ void ParFiniteElementSpace::Build_Dof_TrueDof_Matrix() const // matrix P
          }
          else if (i_offd[i+1] == i_offd[i] + 2)
          {
-            const double *T =
-               ND_DofTransformation::GetFaceTransform(ltori[i]).GetData();
+            const double * T = ND_StatelessDofTransformation
+                               ::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++;
@@ -1455,28 +1454,29 @@ void ParFiniteElementSpace::ExchangeFaceNbrData()
    delete [] requests;
 }
 
-void ParFiniteElementSpace::GetFaceNbrElementVDofs(
-   int i, Array<int> &vdofs, DofTransformation &doftrans) const
-{
-   face_nbr_element_dof.GetRow(i, vdofs);
-
-   if (DoFTransArray[GetFaceNbrFE(i)->GetGeomType()])
-   {
-      Array<int> F, Fo;
-      pmesh->GetFaceNbrElementFaces(pmesh->GetNE() + i, F, Fo);
-      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;
+   face_nbr_element_dof.GetRow(i, vdofs);
+
+   DofTransformation *doftrans = NULL;
+   Geometry::Type geom = GetFaceNbrFE(i)->GetGeomType();
+   if (DoFTrans[geom])
+   {
+      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;
+   }
 }
 
 void ParFiniteElementSpace::GetFaceNbrFaceVDofs(int i, Array<int> &vdofs) const
@@ -2235,13 +2235,19 @@ 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();
 
-            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);
+            const auto T = [&fo]()
+            {
+               auto T = ND_StatelessDofTransformation::GetFaceTransform(fo);
+               T(0,0) -= 1;
+               T(1,1) -= 1;
+               return T;
+            }();
+
+            first_row.AddRow(initial_first_row, T(0,0));
+            first_row.AddRow(initial_second_row, T(0,1));
+            second_row.AddRow(initial_first_row, T(1,0));
+            second_row.AddRow(initial_second_row, T(1,1));
 
             first_row.Collapse();
             second_row.Collapse();

fem/pfespace.hpp

@@ -284,17 +284,11 @@ 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 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 in array dofs for i'th element.
+   DofTransformation *GetElementDofs(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 indexes of degrees of freedom for i'th boundary element.
+   DofTransformation *GetBdrElementDofs(int i, Array<int> &dofs) const override;
 
    /** Returns the indexes of the degrees of freedom for i'th face
        including the dofs for the edges and the vertices of the face. */
@@ -388,8 +382,6 @@ 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

@@ -39,9 +39,11 @@ ParGridFunction::ParGridFunction(ParMesh *pmesh, const GridFunction *gf,
 {
    const FiniteElementSpace *glob_fes = gf->FESpace();
    // duplicate the FiniteElementCollection from 'gf'
-   fec = FiniteElementCollection::New(glob_fes->FEColl()->Name());
+   fec.reset(FiniteElementCollection::New(glob_fes->FEColl()->Name()));
    // create a local ParFiniteElementSpace from the global one:
-   fes = pfes = new ParFiniteElementSpace(pmesh, glob_fes, partitioning, fec);
+   fes = pfes = new ParFiniteElementSpace(pmesh, glob_fes, partitioning,
+                                          fec.get());
+   owned_fes.reset(pfes);
    SetSize(pfes->GetVSize());
 
    if (partitioning)
@@ -81,10 +83,10 @@ ParGridFunction::ParGridFunction(ParMesh *pmesh, std::istream &input)
    : GridFunction(pmesh, input)
 {
    // Convert the FiniteElementSpace, fes, to a ParFiniteElementSpace:
-   pfes = new ParFiniteElementSpace(pmesh, fec, fes->GetVDim(),
+   pfes = new ParFiniteElementSpace(pmesh, fec.get(), fes->GetVDim(),
                                     fes->GetOrdering());
-   delete fes;
    fes = pfes;
+   owned_fes.reset(pfes);
 }
 
 void ParGridFunction::Update()
@@ -693,23 +695,7 @@ void ParGridFunction::ProjectBdrCoefficient(
 
 #ifdef MFEM_DEBUG
    Array<int> 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]);
-         }
-      }
-   }
+   pfes->GetEssentialVDofs(attr, ess_vdofs_marker);
    for (int i = 0; i < values_counter.Size(); i++)
    {
       MFEM_ASSERT(pfes->GetLocalTDofNumber(i) == -1 ||

fem/pgridfunc.hpp

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

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.GetBdrElementFaceIndex(T.ElementNo));
+         return get_face_index(mesh.GetBdrElementEdgeIndex(T.ElementNo));
       default:
          MFEM_ABORT("Invalid element type.");
          return -1;

fem/tmop.cpp

@@ -462,52 +462,21 @@ void TMOP_Metric_009::AssembleH(const DenseMatrix &Jpt,
    ie.Assemble_ddI1b(weight, A.GetData());
 }
 
-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();
-}
-
+// mu_14 = |T-I|^2
 double TMOP_Metric_014::EvalW(const DenseMatrix &Jpt) const
 {
-   // mu_14 = |J - I|^2 = I1[J-I].
-   DenseMatrix Mat(Jpt);
-   Mat(0,0) -= 1.0;
-   Mat(1,1) -= 1.0;
+   MFEM_VERIFY(Jtr != NULL,
+               "Requires a target Jacobian, use SetTargetJacobian().");
 
-   ie.SetJacobian(Mat.GetData());
-   return ie.Get_I1();
-}
+   DenseMatrix Id(2,2);
 
-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();
-}
+   Id(0,0) = 1; Id(0,1) = 0;
+   Id(1,0) = 0; Id(1,1) = 1;
 
-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());
+   DenseMatrix Mat(2,2);
+   Mat = Jpt;
+   Mat.Add(-1,Id);
+   return Mat.FNorm2();
 }
 
 double TMOP_Metric_022::EvalW(const DenseMatrix &Jpt) const

fem/tmop.hpp

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

general/communication.cpp

@@ -26,10 +26,6 @@
 #include "sort_pairs.hpp"
 #include "globals.hpp"
 
-#ifdef MFEM_USE_STRUMPACK
-#include <StrumpackConfig.hpp> // STRUMPACK_USE_PTSCOTCH, etc.
-#endif
-
 #include <iostream>
 #include <map>
 
@@ -38,14 +34,6 @@ 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,6 +22,7 @@
 #include "globals.hpp"
 #include <mpi.h>
 
+
 namespace mfem
 {
 
@@ -31,34 +32,10 @@ namespace mfem
 class Mpi
 {
 public:
-   /// 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();
-   }
+   /// 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); }
    /// Finalize MPI (if it has been initialized and not yet already finalized).
    static void Finalize()
    {
@@ -94,19 +71,20 @@ 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 the Mpi singleton.
-   static Mpi &Singleton()
+   /// Initialize MPI
+   static void Init_(int *argc, char ***argv)
    {
+      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 <winsock2.h>
-#include <ws2tcpip.h>
+#include <winsock.h>
 #ifdef _MSC_VER
 typedef int ssize_t;
 // Link with ws2_32.lib
@@ -51,66 +51,47 @@ int isockstream::establish()
 {
    // char myname[129];
    char   myname[] = "localhost";
-   int    sfd;
-   struct addrinfo hints, *res, *rp;
+   int    port;
+   struct sockaddr_in sa;
+   struct hostent *hp;
 
-   memset(&hints, 0, sizeof(hints));
-   hints.ai_family = AF_UNSPEC;
-   hints.ai_socktype = SOCK_STREAM;
-   hints.ai_protocol = 0;
+   memset(&sa, 0, sizeof(struct sockaddr_in));
+   // gethostname(myname, 128);
+   hp= gethostbyname(myname);
 
-   int s = getaddrinfo(myname, NULL, &hints, &res);
-   if (s != 0)
+   if (hp == NULL)
    {
-      mfem::err << "isockstream::establish(): getaddrinfo() failed!\n"
-                << "isockstream::establish(): getaddrinfo() returned: '"
+      mfem::err << "isockstream::establish(): gethostbyname() failed!\n"
+                << "isockstream::establish(): gethostname() returned: '"
                 << myname << "'" << endl;
       error = 1;
       return (-1);
    }
 
-   // loop the list of address structures returned by getaddrinfo()
-   for (rp = res; rp != NULL; rp = rp->ai_next)
+   sa.sin_family= hp->h_addrtype;
+   sa.sin_port= htons(portnum);
+
+   if ((port = socket(AF_INET, SOCK_STREAM, 0)) < 0)
    {
-      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";
+      mfem::err << "isockstream::establish(): socket() failed!" << endl;
+      error = 2;
       return (-1);
    }
 
-   freeaddrinfo(res);
-   listen(sfd, 4);
-   return (sfd);
+   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);
 }
 
 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>      // getaddrinfo
+#include <netdb.h>      // gethostbyname
 #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 <winsock2.h>
-#include <ws2tcpip.h>
+#include <winsock.h>
 #ifdef _MSC_VER
 typedef int ssize_t;
 // Link with ws2_32.lib
@@ -93,7 +93,8 @@ int socketbuf::attach(int sd)
 
 int socketbuf::open(const char hostname[], int port)
 {
-   struct addrinfo     hints, *res, *rp;
+   struct sockaddr_in  sa;
+   struct hostent     *hp;
 
    if (!wsInit_.Initialized())
    {
@@ -104,47 +105,42 @@ int socketbuf::open(const char hostname[], int port)
    setg(NULL, NULL, NULL);
    setp(obuf, obuf + buflen);
 
-   hints.ai_family = AF_UNSPEC;
-   hints.ai_socktype = SOCK_STREAM;
-   hints.ai_protocol = 0;
-
-   int s = getaddrinfo(hostname, NULL, &hints, &res);
-   if (s != 0)
+   hp = gethostbyname(hostname);
+   if (hp == NULL)
    {
       socket_descriptor = -3;
       return -1;
    }
-
-   for (rp = res; rp != NULL; rp = rp->ai_next)
+   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)
    {
-      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,
-                     &on, sizeof(on)) < 0)
-      {
-         closesocket(socket_descriptor);
-         socket_descriptor = -2;
-         return -1;
-      }
-#endif
-
-      if (connect(socket_descriptor, rp->ai_addr, rp->ai_addrlen) < 0)
-      {
-         closesocket(socket_descriptor);
-         socket_descriptor = -2;
-         continue;
-      }
+      return -1;
    }
 
-   freeaddrinfo(res);
+#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;
+   }
+#endif
+
+   if (connect(socket_descriptor,
+               (const struct sockaddr *)&sa, sizeof(sa)) < 0)
+   {
+      closesocket(socket_descriptor);
+      socket_descriptor = -2;
+      return -1;
+   }
    return 0;
 }
 

linalg/strumpack.cpp

@@ -16,471 +16,238 @@
 
 #include "strumpack.hpp"
 
+using namespace std;
+using namespace strumpack;
+
 namespace mfem
 {
 
 STRUMPACKRowLocMatrix::STRUMPACKRowLocMatrix(MPI_Comm comm,
-                                             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)
+                                             int num_loc_rows, int first_loc_row,
+                                             int glob_nrows, int glob_ncols,
+                                             int *I, int *J, double *data)
+   : comm_(comm), A_(NULL)
 {
    // Set mfem::Operator member data
    height = num_loc_rows;
    width  = 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);
+   // 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;
    dist[0] = 0;
-   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);
-
-#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
-   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
+   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;
 }
 
-STRUMPACKRowLocMatrix::STRUMPACKRowLocMatrix(const Operator &op,
-                                             bool sym_sparse)
+STRUMPACKRowLocMatrix::STRUMPACKRowLocMatrix(const HypreParMatrix & hypParMat)
+   : comm_(hypParMat.GetComm()),
+     A_(NULL)
 {
-   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);
+   hypre_ParCSRMatrix * parcsr_op =
+      (hypre_ParCSRMatrix *)const_cast<HypreParMatrix&>(hypParMat);
 
-   // 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;
+   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
-   HYPRE_BigInt    *Jptr   = csr_op->big_j;
-#else
-   HYPRE_Int       *Jptr   = csr_op->j;
+   // 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
-   double          *data   = csr_op->data;
 
-   HYPRE_BigInt fst_row = parcsr_op->first_row_index;
-   HYPRE_Int    m_loc   = csr_op->num_rows;
+   height = csr_op->num_rows;
+   width  = 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);
+   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;
-   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);
+   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;
 
-#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
+   // Everything has been copied or abducted so delete the structure
    hypre_CSRMatrixDestroy(csr_op);
 }
 
 STRUMPACKRowLocMatrix::~STRUMPACKRowLocMatrix()
 {
-   delete A_;
+   // Delete the struct
+   if ( A_ != NULL ) { 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)
+STRUMPACKSolver::STRUMPACKSolver( int argc, char* argv[], MPI_Comm comm )
+   : comm_(comm),
+     APtr_(NULL),
+     solver_(NULL)
 {
-   solver_ = new STRUMPACKSolverType(comm, argc, argv, false);
+   this->Init(argc, argv);
 }
 
-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)
+STRUMPACKSolver::STRUMPACKSolver( STRUMPACKRowLocMatrix & A )
+   : comm_(A.GetComm()),
+     APtr_(&A),
+     solver_(NULL)
 {
-   solver_ = new STRUMPACKSolverType(A.GetComm(), argc, argv, false);
-   SetOperator(A);
+   height = A.Height();
+   width  = A.Width();
+
+   this->Init(0, NULL);
 }
 
-template <typename STRUMPACKSolverType>
-STRUMPACKSolverBase<STRUMPACKSolverType>::
-~STRUMPACKSolverBase()
+STRUMPACKSolver::~STRUMPACKSolver()
 {
-   delete solver_;
+   if ( solver_ != NULL ) { delete solver_; }
 }
 
-template <typename STRUMPACKSolverType>
-void STRUMPACKSolverBase<STRUMPACKSolverType>::
-SetFromCommandLine()
+void STRUMPACKSolver::Init( int argc, char* argv[] )
 {
-   solver_->options().set_from_command_line();
+   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);
 }
 
-template <typename STRUMPACKSolverType>
-void STRUMPACKSolverBase<STRUMPACKSolverType>::
-SetPrintFactorStatistics(bool print_stat)
+void STRUMPACKSolver::SetFromCommandLine( )
+{
+   solver_->options().set_from_command_line( );
+}
+
+void STRUMPACKSolver::SetPrintFactorStatistics( bool print_stat )
 {
    factor_verbose_ = print_stat;
 }
 
-template <typename STRUMPACKSolverType>
-void STRUMPACKSolverBase<STRUMPACKSolverType>::
-SetPrintSolveStatistics(bool print_stat)
+void STRUMPACKSolver::SetPrintSolveStatistics( bool print_stat )
 {
    solve_verbose_ = print_stat;
 }
 
-template <typename STRUMPACKSolverType>
-void STRUMPACKSolverBase<STRUMPACKSolverType>
-::SetRelTol(double rtol)
+void STRUMPACKSolver::SetKrylovSolver( strumpack::KrylovSolver method )
 {
-   solver_->options().set_rel_tol(rtol);
+   solver_->options().set_Krylov_solver( method );
 }
 
-template <typename STRUMPACKSolverType>
-void STRUMPACKSolverBase<STRUMPACKSolverType>
-::SetAbsTol(double atol)
+void STRUMPACKSolver::SetReorderingStrategy( strumpack::ReorderingStrategy
+                                             method )
 {
-   solver_->options().set_abs_tol(atol);
+   solver_->options().set_reordering_method( method );
 }
 
-template <typename STRUMPACKSolverType>
-void STRUMPACKSolverBase<STRUMPACKSolverType>
-::SetMaxIter(int max_it)
+void STRUMPACKSolver::DisableMatching( )
 {
-   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);
+#if STRUMPACK_VERSION_MAJOR >= 3
+   solver_->options().set_matching( strumpack::MatchingJob::NONE );
 #else
-   switch (type)
+   solver_->options().set_mc64job( strumpack::MC64Job::NONE );
+#endif
+}
+
+void STRUMPACKSolver::EnableMatching( )
+{
+#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 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;
+      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("Invalid compression type for STRUMPACK version " <<
-                    STRUMPACK_VERSION_MAJOR << "!");
-         break;
+      {
+         MFEM_ABORT("STRUMPACK: 'factor()' error code = " << ret);
+      }
    }
-#endif
+   solver_->options().set_verbose( solve_verbose_ );
+   solver_->solve(xPtr, yPtr);
+
 }
 
-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)
+void STRUMPACKSolver::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!");
+   APtr_ = dynamic_cast<const STRUMPACKRowLocMatrix*>(&op);
+   if ( APtr_ == NULL )
+   {
+      mfem_error("STRUMPACKSolver::SetOperator : not STRUMPACKRowLocMatrix!");
+   }
+
+   solver_->set_matrix( *(APtr_->getA()) );
 
    // 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,14 +16,12 @@
 
 #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
 {
@@ -36,80 +34,63 @@ 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, HYPRE_BigInt first_loc_row,
-                         HYPRE_BigInt glob_nrows, HYPRE_BigInt glob_ncols,
-                         int *I, HYPRE_BigInt *J, double *data,
-                         bool sym_sparse = false);
+                         int num_loc_rows, int first_loc_row,
+                         int glob_nrows, int glob_ncols,
+                         int *I, int *J, double *data);
 
    /** 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 Operator &op, bool sym_sparse = false);
+   STRUMPACKRowLocMatrix(const HypreParMatrix & hypParMat);
 
    ~STRUMPACKRowLocMatrix();
 
    void Mult(const Vector &x, Vector &y) const
    {
-      MFEM_ABORT("STRUMPACKRowLocMatrix::Mult: Matrix vector products are not "
-                 "supported!");
+      mfem_error("STRUMPACKRowLocMatrix::Mult(...)\n"
+                 "  matrix vector products are not supported.");
    }
 
-   MPI_Comm GetComm() const { return A_->comm(); }
+   MPI_Comm GetComm() const { return comm_; }
 
-   strumpack::CSRMatrixMPI<double, HYPRE_BigInt> *GetA() const { return A_; }
+   strumpack::CSRMatrixMPI<double,int>* getA() const { return A_; }
 
 private:
-   strumpack::CSRMatrixMPI<double, HYPRE_BigInt> *A_;
-};
+   MPI_Comm   comm_;
+   strumpack::CSRMatrixMPI<double,int>* A_;
+
+}; // mfem::STRUMPACKRowLocMatrix
 
 /** 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
 */
-template <typename STRUMPACKSolverType>
-class STRUMPACKSolverBase : public Solver
+class STRUMPACKSolver : public mfem::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.
-   virtual ~STRUMPACKSolverBase();
+   ~STRUMPACKSolver( void );
 
    // Factor and solve the linear system y = Op^{-1} x.
-   void Mult(const Vector &x, Vector &y) const;
-   void ArrayMult(const Array<const Vector *> &X, Array<Vector *> &Y) const;
+   void Mult( const Vector & x, 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);
-
-   // 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();
+   void SetFromCommandLine( );
+   void SetPrintFactorStatistics( bool print_stat );
+   void SetPrintSolveStatistics( bool print_stat );
+   void SetRelTol( double rtol );
+   void SetAbsTol( double atol );
 
    /**
     * STRUMPACK is an (approximate) direct solver. It can be used as a direct
@@ -119,153 +100,70 @@ 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:
-    *    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
+    *    METIS, PARMETIS, SCOTCH, PTSCOTCH, RCM
     */
-   void SetReorderingStrategy(strumpack::ReorderingStrategy method);
+   void SetReorderingStrategy( strumpack::ReorderingStrategy method );
 
    /**
-    * Configure static pivoting for stability. The static pivoting in STRUMPACK
+    * Disable 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 SetMatching(strumpack::MatchingJob job);
+   void DisableMatching();
 
    /**
-    * 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.
+    * 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.
     */
-   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);
+   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();
 #endif
 
 private:
-   // Helper method for calling the STRUMPACK factoriation routine.
-   void FactorInternal() const;
+   void Init( int argc, char* argv[] );
 
 protected:
-   const STRUMPACKRowLocMatrix *APtr_;
-   STRUMPACKSolverType         *solver_;
+
+   MPI_Comm      comm_;
+   int           numProcs_;
+   int           myid_;
 
    bool factor_verbose_;
    bool solve_verbose_;
-   bool reorder_reuse_;
 
-   mutable Vector rhs_, sol_;
-   mutable int    nrhs_;
-};
+   const STRUMPACKRowLocMatrix * APtr_;
+   strumpack::StrumpackSparseSolverMPIDist<double,int> * solver_;
 
-class STRUMPACKSolver :
-   public STRUMPACKSolverBase<strumpack::
-   SparseSolverMPIDist<double, HYPRE_BigInt>>
-{
-public:
-   // Constructor with MPI_Comm parameter.
-   STRUMPACKSolver(MPI_Comm comm);
+}; // mfem::STRUMPACKSolver class
 
-   // 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
+} // mfem namespace
 
 #endif // MFEM_USE_MPI
 #endif // MFEM_USE_STRUMPACK

linalg/superlu.cpp

@@ -650,7 +650,6 @@ 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,6 +35,10 @@
 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()
@@ -440,7 +444,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;
+   int CheckFinite() const { return mfem::CheckFinite(HostRead(), size); }
 
    /// Destroys vector.
    virtual ~Vector();
@@ -490,8 +494,6 @@ 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;
@@ -502,11 +504,6 @@ 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,12 +384,6 @@ 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)
 {
@@ -434,6 +428,19 @@ 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);
@@ -494,12 +501,6 @@ 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();
@@ -1101,7 +1102,7 @@ void Mesh::ApplyLocalSlaveTransformation(FaceElementTransformations &FT,
 FaceElementTransformations *Mesh::GetBdrFaceTransformations(int BdrElemNo)
 {
    FaceElementTransformations *tr;
-   int fn = GetBdrElementFaceIndex(BdrElemNo);
+   int fn = GetBdrFace(BdrElemNo);
 
    // Check if the face is interior, shared, or nonconforming.
    if (FaceIsTrueInterior(fn) || faces_info[fn].NCFace >= 0)
@@ -1116,6 +1117,24 @@ 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;
@@ -1429,7 +1448,7 @@ Array<int> Mesh::GetFaceToBdrElMap() const
    face_to_be = -1;
    for (int i = 0; i < NumOfBdrElements; i++)
    {
-      face_to_be[GetBdrElementFaceIndex(i)] = i;
+      face_to_be[GetBdrElementEdgeIndex(i)] = i;
    }
    return face_to_be;
 }
@@ -1520,6 +1539,7 @@ void Mesh::Destroy()
    faces.DeleteAll();
    faces_info.DeleteAll();
    nc_faces_info.DeleteAll();
+   be_to_edge.DeleteAll();
    be_to_face.DeleteAll();
 
    // TODO:
@@ -1886,7 +1906,12 @@ int Mesh::AddBdrPoint(int v, int attr)
 
 void Mesh::GenerateBoundaryElements()
 {
-   for (int i = 0; i < boundary.Size(); i++)
+   int i, j;
+   Array<int> &be2face = (Dim == 2) ? be_to_edge : be_to_face;
+
+   // GenerateFaces();
+
+   for (i = 0; i < boundary.Size(); i++)
    {
       FreeElement(boundary[i]);
    }
@@ -1899,24 +1924,22 @@ void Mesh::GenerateBoundaryElements()
 
    // count the 'NumOfBdrElements'
    NumOfBdrElements = 0;
-   for (int i = 0; i < faces_info.Size(); i++)
+   for (i = 0; i < faces_info.Size(); i++)
    {
       if (faces_info[i].Elem2No < 0) { NumOfBdrElements++; }
    }
 
-   // Add the boundary elements
    boundary.SetSize(NumOfBdrElements);
-   be_to_face.SetSize(NumOfBdrElements);
-   for (int i = 0, j = 0; i < faces_info.Size(); i++)
+   be2face.SetSize(NumOfBdrElements);
+   for (j = i = 0; i < faces_info.Size(); i++)
    {
       if (faces_info[i].Elem2No < 0)
       {
          boundary[j] = faces[i]->Duplicate(this);
-         be_to_face[j++] = i;
+         be2face[j++] = i;
       }
    }
-
-   // Note: in 3D, 'bel_to_edge' is destroyed but it's not updated.
+   // In 3D, 'bel_to_edge' is destroyed but it's not updated.
 }
 
 void Mesh::FinalizeCheck()
@@ -1946,7 +1969,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);
+      NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
       GenerateFaces();
       CheckBdrElementOrientation();
    }
@@ -1974,7 +1997,7 @@ void Mesh::FinalizeQuadMesh(int generate_edges, int refine,
    if (generate_edges)
    {
       el_to_edge = new Table;
-      NumOfEdges = GetElementToEdgeTable(*el_to_edge);
+      NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
       GenerateFaces();
       CheckBdrElementOrientation();
    }
@@ -2276,7 +2299,8 @@ void Mesh::ReorderElements(const Array<int> &ordering, bool reorder_vertices)
    // - edge_vertex - no need to rebuild
    // - geom_factors - no need to rebuild
 
-   // - be_to_face
+   // - be_to_edge  - 2D only
+   // - be_to_face  - 3D only
 
    // - Nodes
 
@@ -2362,9 +2386,9 @@ void Mesh::ReorderElements(const Array<int> &ordering, bool reorder_vertices)
 
    if (Dim > 1)
    {
-      // generate el_to_edge, be_to_face (2D), bel_to_edge (3D)
+      // generate el_to_edge, be_to_edge (2D), bel_to_edge (3D)
       el_to_edge = new Table;
-      NumOfEdges = GetElementToEdgeTable(*el_to_edge);
+      NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
    }
    if (Dim > 2)
    {
@@ -2754,8 +2778,8 @@ void Mesh::DoNodeReorder(DSTable *old_v_to_v, Table *old_elem_vert)
    }
    if (el_to_edge)
    {
-      // update 'el_to_edge', 'be_to_face' (2D), 'bel_to_edge' (3D)
-      NumOfEdges = GetElementToEdgeTable(*el_to_edge);
+      // update 'el_to_edge', 'be_to_edge' (2D), 'bel_to_edge' (3D)
+      NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
       if (Dim == 2)
       {
          // update 'faces' and 'faces_info'
@@ -2832,7 +2856,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);
+      NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
    }
    else
    {
@@ -2867,7 +2891,7 @@ void Mesh::FinalizeWedgeMesh(int generate_edges, int refine,
    if (generate_edges == 1)
    {
       el_to_edge = new Table;
-      NumOfEdges = GetElementToEdgeTable(*el_to_edge);
+      NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
    }
    else
    {
@@ -2899,7 +2923,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);
+      NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
    }
    else
    {
@@ -2978,7 +3002,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);
+      NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
       if (Dim == 2)
       {
          GenerateFaces(); // 'Faces' in 2D refers to the edges
@@ -2998,17 +3022,8 @@ 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)
@@ -3540,7 +3555,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);
+      NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
       GenerateFaces();
       CheckBdrElementOrientation();
    }
@@ -3598,11 +3613,6 @@ 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);
 }
@@ -3656,6 +3666,9 @@ 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++)
@@ -5256,7 +5269,7 @@ void Mesh::UpdateNURBS()
 
    if (el_to_edge)
    {
-      NumOfEdges = GetElementToEdgeTable(*el_to_edge);
+      NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
    }
 
    if (el_to_face)
@@ -6156,15 +6169,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);
+         NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
          GenerateFaces(); // 'Faces' in 2D refers to the edges
       }
       for (int i = 0; i < NumOfBdrElements; i++)
       {
-         if (faces_info[be_to_face[i]].Elem2No < 0) // boundary face
+         if (faces_info[be_to_edge[i]].Elem2No < 0) // boundary face
          {
             int *bv = boundary[i]->GetVertices();
-            int *fv = faces[be_to_face[i]]->GetVertices();
+            int *fv = faces[be_to_edge[i]]->GetVertices();
             if (bv[0] != fv[0])
             {
                if (fix_it)
@@ -6415,7 +6428,7 @@ void Mesh::GetBdrElementEdges(int i, Array<int> &edges, Array<int> &cor) const
    {
       edges.SetSize(1);
       cor.SetSize(1);
-      edges[0] = be_to_face[i];
+      edges[0] = be_to_edge[i];
       const int *v = boundary[i]->GetVertices();
       cor[0] = (v[0] < v[1]) ? (1) : (-1);
    }
@@ -6656,7 +6669,7 @@ Array<int> Mesh::FindFaceNeighbors(const int elem) const
 
 void Mesh::GetBdrElementFace(int i, int *f, int *o) const
 {
-   *f = GetBdrElementFaceIndex(i);
+   *f = GetBdrElementEdgeIndex(i);
 
    const int *fv = (Dim > 1) ? faces[*f]->GetVertices() : NULL;
    const int *bv = boundary[i]->GetVertices();
@@ -6673,9 +6686,21 @@ 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 = GetBdrElementFaceIndex(bdr_el);
+   int fid = GetBdrElementEdgeIndex(bdr_el);
 
    const FaceInfo &fi = faces_info[fid];
    MFEM_ASSERT(fi.Elem1Inf % 64 == 0, "internal error"); // orientation == 0
@@ -6697,7 +6722,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 = GetBdrElementFaceIndex(bdr_el);
+   int fid = GetBdrElementEdgeIndex(bdr_el);
 
    const FaceInfo &fi = faces_info[fid];
    MFEM_ASSERT(fi.Elem1Inf % 64 == 0, "internal error"); // orientation == 0
@@ -6823,7 +6848,7 @@ void Mesh::GetVertexToVertexTable(DSTable &v_to_v) const
    }
 }
 
-int Mesh::GetElementToEdgeTable(Table &e_to_f)
+int Mesh::GetElementToEdgeTable(Table & e_to_f, Array<int> &be_to_f)
 {
    int i, NumberOfEdges;
 
@@ -6838,11 +6863,11 @@ int Mesh::GetElementToEdgeTable(Table &e_to_f)
    if (Dim == 2)
    {
       // Initialize the indices for the boundary elements.
-      be_to_face.SetSize(NumOfBdrElements);
+      be_to_f.SetSize(NumOfBdrElements);
       for (i = 0; i < NumOfBdrElements; i++)
       {
          const int *v = boundary[i]->GetVertices();
-         be_to_face[i] = v_to_v(v[0], v[1]);
+         be_to_f[i] = v_to_v(v[0], v[1]);
       }
    }
    else if (Dim == 3)
@@ -7445,7 +7470,7 @@ void Mesh::ReorientTetMesh()
       GenerateFaces();
       if (el_to_edge)
       {
-         NumOfEdges = GetElementToEdgeTable(*el_to_edge);
+         NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
       }
    }
    else
@@ -8480,7 +8505,7 @@ void Mesh::UniformRefinement2D_base(bool update_nodes)
    if (el_to_edge == NULL)
    {
       el_to_edge = new Table;
-      NumOfEdges = GetElementToEdgeTable(*el_to_edge);
+      NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
    }
 
    int quad_counter = 0;
@@ -8569,8 +8594,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_face[i], attr);
-      new_boundary[j++] = new Segment(oedge+be_to_face[i], v[1], attr);
+      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);
 
       FreeElement(boundary[i]);
    }
@@ -8610,7 +8635,7 @@ void Mesh::UniformRefinement2D_base(bool update_nodes)
    NumOfBdrElements = 2 * NumOfBdrElements;
    NumOfFaces       = 0;
 
-   NumOfEdges = GetElementToEdgeTable(*el_to_edge);
+   NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
    GenerateFaces();
 
    last_operation = Mesh::REFINE;
@@ -8640,7 +8665,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);
+      NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
    }
 
    if (el_to_face == NULL)
@@ -9331,7 +9356,7 @@ void Mesh::UniformRefinement3D_base(Array<int> *f2qf_ptr, DSTable *v_to_v_p,
    CheckBdrElementOrientation(false);
 #endif
 
-   NumOfEdges = GetElementToEdgeTable(*el_to_edge);
+   NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
 
    last_operation = Mesh::REFINE;
    sequence++;
@@ -9465,7 +9490,7 @@ void Mesh::LocalRefinement(const Array<int> &marked_el, int type)
 
       if (el_to_edge != NULL)
       {
-         NumOfEdges = GetElementToEdgeTable(*el_to_edge);
+         NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
          GenerateFaces();
       }
 
@@ -9561,7 +9586,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);
+         NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
       }
       if (el_to_face != NULL)
       {
@@ -9752,7 +9777,7 @@ void Mesh::InitFromNCMesh(const NCMesh &ncmesh_)
    if (Dim > 1)
    {
       el_to_edge = new Table;
-      NumOfEdges = GetElementToEdgeTable(*el_to_edge);
+      NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
    }
    if (Dim > 2)
    {
@@ -9799,6 +9824,7 @@ 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);
@@ -12253,9 +12279,9 @@ void Mesh::RemoveUnusedVertices()
    DeleteTables();
    if (Dim > 1)
    {
-      // generate el_to_edge, be_to_face (2D), bel_to_edge (3D)
+      // generate el_to_edge, be_to_edge (2D), bel_to_edge (3D)
       el_to_edge = new Table;
-      NumOfEdges = GetElementToEdgeTable(*el_to_edge);
+      NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
    }
    if (Dim > 2)
    {
@@ -12286,7 +12312,7 @@ void Mesh::RemoveInternalBoundaries()
    int new_bel_to_edge_nnz = 0;
    for (int i = 0; i < GetNBE(); i++)
    {
-      if (FaceIsInterior(GetBdrElementFaceIndex(i)))
+      if (FaceIsInterior(GetBdrElementEdgeIndex(i)))
       {
          FreeElement(boundary[i]);
       }
@@ -12303,24 +12329,32 @@ void Mesh::RemoveInternalBoundaries()
    if (num_bdr_elem == GetNBE()) { return; }
 
    Array<Element *> new_boundary(num_bdr_elem);
-   Array<int> new_be_to_face;
+   Array<int> new_be_to_edge, new_be_to_face;
    Table *new_bel_to_edge = NULL;
    new_boundary.SetSize(0);
-   new_be_to_face.Reserve(num_bdr_elem);
-   if (Dim == 3)
+   if (Dim == 2)
    {
+      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(GetBdrElementFaceIndex(i)))
+      if (!FaceIsInterior(GetBdrElementEdgeIndex(i)))
       {
          new_boundary.Append(boundary[i]);
-         int row = new_be_to_face.Size();
-         new_be_to_face.Append(be_to_face[i]);
-         if (Dim == 3)
+         if (Dim == 2)
          {
+            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);
@@ -12336,10 +12370,13 @@ void Mesh::RemoveInternalBoundaries()
    NumOfBdrElements = new_boundary.Size();
    mfem::Swap(boundary, new_boundary);
 
-   mfem::Swap(be_to_face, new_be_to_face);
-
-   if (Dim == 3)
+   if (Dim == 2)
    {
+      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_face; // faces = vertices (1D), edges (2D), faces (3D)
-
-   Table *bel_to_edge;    // for 3D only
+   Array<int> be_to_edge;  // for 2D
+   Table *bel_to_edge;     // for 3D
+   Array<int> be_to_face;
 
    // 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 &);
+   int GetElementToEdgeTable(Table &, Array<int> &);
 
    /// Used in GenerateFaces()
    void AddPointFaceElement(int lf, int gf, int el);
@@ -1367,6 +1367,11 @@ 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.
 
@@ -1387,19 +1392,8 @@ public:
        @sa GetBdrElementAdjacentElement() */
    void GetBdrElementAdjacentElement2(int bdr_el, int &el, int &info) 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); }
+   /// Return the local face index for the given boundary face.
+   int GetBdrFace(int BdrElemNo) const;
 
    /// @}
 

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);
+         NumOfEdges = Mesh::GetElementToEdgeTable(*el_to_edge, be_to_edge);
       }
 
       STable3D *faces_tbl = NULL;
@@ -197,19 +197,6 @@ 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
@@ -448,7 +435,7 @@ int ParMesh::BuildLocalBoundary(const Mesh& mesh, const int* partitioning,
 
       for (int i = 0; i < mesh.GetNBE(); i++)
       {
-         int edge = mesh.GetBdrElementFaceIndex(i);
+         int edge = mesh.GetBdrElementEdgeIndex(i);
          int el1 = edge_element->GetRow(edge)[0];
          if (partitioning[el1] == MyRank)
          {
@@ -464,7 +451,7 @@ int ParMesh::BuildLocalBoundary(const Mesh& mesh, const int* partitioning,
       boundary.SetSize(nbdry);
       for (int i = 0; i < mesh.GetNBE(); i++)
       {
-         int edge = mesh.GetBdrElementFaceIndex(i);
+         int edge = mesh.GetBdrElementEdgeIndex(i);
          int el1 = edge_element->GetRow(edge)[0];
          if (partitioning[el1] == MyRank)
          {
@@ -3333,7 +3320,7 @@ void ParMesh::ReorientTetMesh()
       GenerateFaces();
       if (el_to_edge)
       {
-         NumOfEdges = GetElementToEdgeTable(*el_to_edge);
+         NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
       }
    }
    else
@@ -3568,7 +3555,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);
+         NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
       }
    } //  'if (Dim == 3)'
 
@@ -3806,7 +3793,7 @@ void ParMesh::LocalRefinement(const Array<int> &marked_el, int type)
 
       if (el_to_edge != NULL)
       {
-         NumOfEdges = GetElementToEdgeTable(*el_to_edge);
+         NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
          GenerateFaces();
       }
    } //  'if (Dim == 2)'
@@ -3868,8 +3855,6 @@ void ParMesh::NonconformingRefinement(const Array<Refinement> &refinements,
                  "serial Mesh)");
    }
 
-   ResetLazyData();
-
    DeleteFaceNbrData();
 
    // NOTE: no check of !refinements.Size(), in parallel we would have to reduce
@@ -5301,7 +5286,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,6 +1028,9 @@ 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,13 +235,10 @@ ParSubMesh::ParSubMesh(const ParMesh &parent, SubMesh::From from,
 
    // Add boundaries
    {
-      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; }
-      }();
+      int num_of_faces_or_edges =
+         (Dim == 3) ? NumOfFaces :
+         ((Dim == 2) ? NumOfEdges : NumOfVertices);
+      Array<int> &be2face = (Dim == 2) ? be_to_edge : be_to_face;
 
       if (Dim == 3)
       {
@@ -252,7 +249,7 @@ ParSubMesh::ParSubMesh(const ParMesh &parent, SubMesh::From from,
       }
 
       NumOfBdrElements = 0;
-      for (int i = 0; i < num_codim_1; i++)
+      for (int i = 0; i < num_of_faces_or_edges; i++)
       {
          if (GetFaceInformation(i).IsBoundary())
          {
@@ -261,17 +258,14 @@ ParSubMesh::ParSubMesh(const ParMesh &parent, SubMesh::From from,
       }
 
       boundary.SetSize(NumOfBdrElements);
-      be_to_face.SetSize(NumOfBdrElements);
+      be2face.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_codim_1; i++)
+      for (int i = 0, j = 0; i < num_of_faces_or_edges; 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]] :
@@ -289,7 +283,7 @@ ParSubMesh::ParSubMesh(const ParMesh &parent, SubMesh::From from,
             {
                boundary[j]->SetAttribute(SubMesh::GENERATED_ATTRIBUTE);
             }
-            ++j;
+            be2face[j++] = i;
          }
       }
    }
@@ -329,7 +323,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);
+      NumOfEdges = GetElementToEdgeTable(*el_to_edge, be_to_edge);
    }
 
    SetAttributes();

mesh/submesh/ptransfermap.cpp

@@ -317,7 +317,8 @@ ParTransferMap::CorrectFaceOrientations(const ParFiniteElementSpace &fes,
 
    if (parent_face_ori.Size() == 0) { return; }
 
-   DofTransformation doftrans(fes.GetVDim(), fes.GetOrdering());
+   VDofTransformation vdoftrans(fes.GetVDim(),
+                                fes.GetOrdering());
 
    int dim = mesh->Dimension();
    bool face = (dim == 3);
@@ -331,13 +332,17 @@ 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);;
 
-      if (!fec->DofTransformationForGeometry(geom)) { continue; }
-      doftrans.SetDofTransformation(*fec->DofTransformationForGeometry(geom));
+      StatelessDofTransformation * doftrans =
+         fec->DofTransformationForGeometry(geom);
+
+      if (doftrans == NULL) { continue; }
+
+      vdoftrans.SetDofTransformation(*doftrans);
 
       Fo[0] = parent_face_ori[i];
-      doftrans.SetFaceOrientations(Fo);
+      vdoftrans.SetFaceOrientations(Fo);
 
       if (face)
       {
@@ -351,12 +356,12 @@ ParTransferMap::CorrectFaceOrientations(const ParFiniteElementSpace &fes,
       if (sub_to_parent_map)
       {
          src.GetSubVector(vdofs, face_vector);
-         doftrans.TransformPrimal(face_vector);
+         vdoftrans.TransformPrimal(face_vector);
       }
       else
       {
          dst.GetSubVector(vdofs, face_vector);
-         doftrans.InvTransformPrimal(face_vector);
+         vdoftrans.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_[GetBdrElementFaceIndex(i)]];
+         int pbeid = parent_face_to_be[parent_face_ids_[GetBdrFace(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_[GetBdrElementFaceIndex(i)]];
+            int pbeid = parent_face_to_be[parent_edge_ids_[GetBdrFace(i)]];
             if (pbeid != -1)
             {
                int attr = parent.GetBdrElement(pbeid)->GetAttribute();

mesh/submesh/transfermap.cpp

@@ -241,7 +241,8 @@ void TransferMap::CorrectFaceOrientations(const FiniteElementSpace &fes,
 
    if (parent_face_ori.Size() == 0) { return; }
 
-   DofTransformation doftrans(fes.GetVDim(), fes.GetOrdering());
+   VDofTransformation vdoftrans(fes.GetVDim(),
+                                fes.GetOrdering());
 
    int dim = mesh->Dimension();
    bool face = (dim == 3);
@@ -255,13 +256,17 @@ 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);;
 
-      if (!fec->DofTransformationForGeometry(geom)) { continue; }
-      doftrans.SetDofTransformation(*fec->DofTransformationForGeometry(geom));
+      StatelessDofTransformation * doftrans =
+         fec->DofTransformationForGeometry(geom);
+
+      if (doftrans == NULL) { continue; }
+
+      vdoftrans.SetDofTransformation(*doftrans);
 
       Fo[0] = parent_face_ori[i];
-      doftrans.SetFaceOrientations(Fo);
+      vdoftrans.SetFaceOrientations(Fo);
 
       if (face)
       {
@@ -275,12 +280,12 @@ void TransferMap::CorrectFaceOrientations(const FiniteElementSpace &fes,
       if (sub_to_parent_map)
       {
          src.GetSubVector(vdofs, face_vector);
-         doftrans.TransformPrimal(face_vector);
+         vdoftrans.TransformPrimal(face_vector);
       }
       else
       {
          dst.GetSubVector(vdofs, face_vector);
-         doftrans.InvTransformPrimal(face_vector);
+         vdoftrans.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
+//     mesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 14 -tid 4 -ni 100 -bnd -qt 1 -qo 8 -fd
 //   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 = GetBdrElementFaceIndex(i);
+         const int fn = GetBdrElementEdgeIndex(i);
          MFEM_VERIFY(!FaceIsTrueInterior(fn),"");
          Array<int> vertices;
          GetFaceVertices(fn, vertices);

miniapps/meshing/pmesh-fitting.cpp

@@ -359,6 +359,7 @@ 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);
@@ -431,6 +432,7 @@ 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);
@@ -628,7 +630,6 @@ 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
+//     mpirun -np 4 pmesh-optimizer -m square01.mesh -o 2 -rs 2 -mid 14 -tid 4 -ni 200 -bnd -qt 1 -qo 8 -fd
 //   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 = GetBdrElementFaceIndex(i);
+         const int fn = GetBdrElementEdgeIndex(i);
          MFEM_VERIFY(!FaceIsTrueInterior(fn),"");
          Array<int> vertices;
          GetFaceVertices(fn, vertices);

miniapps/multidomain/multidomain.cpp

@@ -322,19 +322,17 @@ int main(int argc, char *argv[])
 
    char vishost[] = "localhost";
    int  visport   = 19916;
-   socketstream cyl_sol_sock;
+   socketstream cyl_sol_sock(vishost, visport);
    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;
+   socketstream block_sol_sock(vishost, visport);
    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,13 +281,12 @@ int main(int argc, char *argv[])
 #ifdef MFEM_USE_STRUMPACK
       if (sp_solver)
       {
-         STRUMPACKSolver * strumpack = new STRUMPACKSolver(MPI_COMM_WORLD, argc, argv);
+         STRUMPACKSolver * strumpack = new STRUMPACKSolver(argc, argv, MPI_COMM_WORLD);
          strumpack->SetPrintFactorStatistics(true);
          strumpack->SetPrintSolveStatistics(false);
          strumpack->SetKrylovSolver(strumpack::KrylovSolver::DIRECT);
          strumpack->SetReorderingStrategy(strumpack::ReorderingStrategy::METIS);
-         strumpack->SetMatching(strumpack::MatchingJob::NONE);
-         strumpack->SetCompression(strumpack::CompressionType::NONE);
+         strumpack->DisableMatching();
          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.GetBdrElementFaceIndex(i), dofs);
+               pfes_sltn->GetFaceDofs(pmesh.GetBdrFace(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.GetBdrElementFaceIndex(i), dofs);
+               pfes_sltn->GetFaceDofs(pmesh.GetBdrFace(i), dofs);
                sface_dof_list.Append(dofs);
             }
          }

miniapps/solvers/lor_solvers.cpp

@@ -57,9 +57,6 @@
 //    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
@@ -113,15 +110,12 @@ int main(int argc, char *argv[])
    double kappa = (order+1)*(order+1); // Penalty used for DG discretizations
 
    Mesh mesh(mesh_file, 1, 1);
-   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.");
+   int dim = mesh.Dimension();
+   MFEM_VERIFY(dim == 2 || dim == 3, "Spatial dimension must be 2 or 3.");
    for (int l = 0; l < ref_levels; l++) { mesh.UniformRefinement(); }
 
    FunctionCoefficient f_coeff(f(1.0)), u_coeff(u);
-   VectorFunctionCoefficient f_vec_coeff(sdim, f_vec(RT)),
-                             u_vec_coeff(sdim, u_vec);
+   VectorFunctionCoefficient f_vec_coeff(dim, f_vec(RT)), u_vec_coeff(dim, u_vec);
 
    int b1 = BasisType::GaussLobatto, b2 = BasisType::IntegratedGLL;
    unique_ptr<FiniteElementCollection> fec;
@@ -155,9 +149,8 @@ int main(int argc, char *argv[])
       a.AddInteriorFaceIntegrator(new DGDiffusionIntegrator(-1.0, kappa));
       a.AddBdrFaceIntegrator(new DGDiffusionIntegrator(-1.0, kappa));
    }
-   // 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); }
+   // TODO: L2 diffusion not implemented with partial assembly
+   if (!L2) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
    a.Assemble();
 
    LinearForm b(&fes);
@@ -195,12 +188,9 @@ int main(int argc, char *argv[])
 
    a.RecoverFEMSolution(X, b, x);
 
-   if (sdim == dim)
-   {
-      double er =
-         (H1 || L2) ? x.ComputeL2Error(u_coeff) : x.ComputeL2Error(u_vec_coeff);
-      cout << "L2 error: " << er << endl;
-   }
+   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,9 +55,6 @@
 //    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
@@ -116,10 +113,8 @@ int main(int argc, char *argv[])
    double kappa = (order+1)*(order+1); // Penalty used for DG discretizations
 
    Mesh serial_mesh(mesh_file, 1, 1);
-   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.");
+   int dim = serial_mesh.Dimension();
+   MFEM_VERIFY(dim == 2 || dim == 3, "Spatial dimension must be 2 or 3.");
    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(); }
@@ -129,8 +124,7 @@ 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(sdim, f_vec(RT)),
-                             u_vec_coeff(sdim, u_vec);
+   VectorFunctionCoefficient f_vec_coeff(dim, f_vec(RT)), u_vec_coeff(dim, u_vec);
 
    int b1 = BasisType::GaussLobatto, b2 = BasisType::IntegratedGLL;
    unique_ptr<FiniteElementCollection> fec;
@@ -165,9 +159,8 @@ int main(int argc, char *argv[])
       a.AddInteriorFaceIntegrator(new DGDiffusionIntegrator(-1.0, kappa));
       a.AddBdrFaceIntegrator(new DGDiffusionIntegrator(-1.0, kappa));
    }
-   // 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); }
+   // TODO: L2 diffusion not implemented with partial assembly
+   if (!L2) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
    a.Assemble();
 
    ParLinearForm b(&fes);
@@ -188,6 +181,7 @@ 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)
    {
@@ -213,12 +207,9 @@ int main(int argc, char *argv[])
 
    a.RecoverFEMSolution(X, b, x);
 
-   if (sdim == dim)
-   {
-      double er =
-         (H1 || L2) ? x.ComputeL2Error(u_coeff) : x.ComputeL2Error(u_vec_coeff);
-      if (Mpi::Root()) { cout << "L2 error: " << er << endl; }
-   }
+   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,13 +33,8 @@ if (MFEM_USE_MPI)
     MAIN plor-transfer.cpp LIBRARIES mfem)
 endif()
 
-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)
+add_mfem_miniapp(check-tmop-metric
+  MAIN check-tmop-metric.cpp LIBRARIES mfem)
 
 if (MFEM_USE_MPI)
    add_mfem_miniapp(nodal-transfer

miniapps/tools/check-tmop-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 tmop-check-metric
+// Compile with: make check-tmop-metric
 //
-// Sample runs:  tmop-check-metric -mid 360
+// Sample runs:  check-tmop-metric -mid 360
 
 #include "mfem.hpp"
 #include <iostream>

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 \
-	tmop-check-metric tmop-metric-magnitude
+	check-tmop-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 tmop-check-metric tmop-metric-magnitude
+	plor-transfer check-tmop-metric tmop-metric-magnitude
 $(foreach app,$(NO_TEST_APPS),$(app)-test-seq $(app)-test-par):
 	@true
 

tests/unit/fem/test_doftrans.cpp

@@ -22,14 +22,13 @@ 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 Tnd(p);
+      ND_TetDofTransformation T(p);
 
       Array<int> ori(4);
       ori[0] = 1;
@@ -37,191 +36,102 @@ TEST_CASE("DoF Transformation Classes",
       ori[2] = 5;
       ori[3] = 1;
 
-      SECTION("VDim == 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")
       {
-         DofTransformation T(Tnd);
-         T.SetFaceOrientations(ori);
+         Vector w;
 
-         Vector u(T.Width());
-         Vector v(T.Width());
-         Vector f(T.Width());
-         Vector ut;
-         Vector vt;
-         Vector ft;
+         ut = u; T.TransformPrimal(ut);
+         w = ut; T.InvTransformPrimal(w);
 
-         u.Randomize(seed);
-         v.Randomize(seed+1);
-         f.Randomize(seed+2);
+         w -= u;
 
-         SECTION("Inverse DoF transformation")
+         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++)
          {
-            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));
+            A.GetColumnReference(i, Ac);
+            Ac.Randomize(seed+i);
          }
       }
-      SECTION("VDim > 1")
+
+      SECTION("Inner product of two primal vectors")
       {
-         Vector v(vdim * Tnd.Width());
-         Vector f(vdim * Tnd.Width());
-         Vector vt;
-         Vector ft;
+         // The matrix A in this case should be regarded as a BilinearForm.
+         DenseMatrix tA;
+         DenseMatrix At;
+         DenseMatrix tAt;
 
-         v.Randomize(seed);
-         f.Randomize(seed+1);
+         ut = u; T.TransformPrimal(ut);
+         vt = v; T.TransformPrimal(vt);
 
-         SECTION("Ordering == byNODES")
-         {
-            DofTransformation T(Tnd, vdim, Ordering::byNODES);
-            T.SetFaceOrientations(ori);
+         At = A; T.TransformDualRows(At);
+         tA = A; T.TransformDualCols(tA);
+         tAt = A; T.TransformDual(tAt);
 
-            SECTION("Inverse DoF transformation")
-            {
-               Vector w;
+         double uAv = A.InnerProduct(v, u);
 
-               vt = v; T.TransformPrimal(vt);
-               w = vt; T.InvTransformPrimal(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;
 
-               w -= v;
+         ft = f; T.TransformDual(ft);
+         vt = v; T.TransformPrimal(vt);
 
-               REQUIRE(w.Norml2() < tol * v.Norml2());
-            }
-            SECTION("Inverse Dual DoF transformation")
-            {
-               Vector w;
+         At = A; T.TransformDualRows(At);
+         tA = A; T.TransformPrimalCols(tA);
+         tAt = At; T.TransformPrimalCols(tAt);
 
-               vt = v; T.TransformDual(vt);
-               w = vt; T.InvTransformDual(w);
+         double fAv = A.InnerProduct(v, f);
 
-               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));
-            }
-         }
+         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));
       }
    }
 }
@@ -236,8 +146,8 @@ TEST_CASE("DoF Transformation Functions",
 
    double tol = 1e-13;
 
-   ND_TetDofTransformation Tndp(p);
-   ND_TetDofTransformation Tndq(q);
+   ND_TetDofTransformation Tp(p);
+   ND_TetDofTransformation Tq(q);
 
    Array<int> ori(4);
    ori[0] = 1;
@@ -245,7 +155,6 @@ TEST_CASE("DoF Transformation Functions",
    ori[2] = 5;
    ori[3] = 1;
 
-   DofTransformation Tp(Tndp), Tq(Tndq);
    Tp.SetFaceOrientations(ori);
    Tq.SetFaceOrientations(ori);
 
@@ -326,4 +235,153 @@ 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,7 +91,6 @@ void TestBatchedLOR()
 {
    const int order = 5;
    const auto mesh_fname = GENERATE(
-                              "../../data/star-surf.mesh",
                               "../../data/star-q3.mesh",
                               "../../data/fichera-q3.mesh"
                            );
@@ -125,11 +124,8 @@ void TestBatchedLOR()
    IntegrationRules irs(0, Quadrature1D::GaussLobatto);
    const IntegrationRule &ir = irs.Get(mesh.GetElementGeometry(0), 1);
    const GeometricFactors::FactorFlags dets = GeometricFactors::DETERMINANTS;
-   if (mesh.Dimension() == mesh.SpaceDimension())
-   {
-      REQUIRE(
-         lor.GetFESpace().GetMesh()->GetGeometricFactors(ir, dets)->detJ.Min() > 0.0);
-   }
+   REQUIRE(lor.GetFESpace().GetMesh()->GetGeometricFactors(ir, dets)->detJ.Min()
+           > 0.0);
 
    lor.LegacyAssembleSystem(a, ess_dofs);
    SparseMatrix A1 = lor.GetAssembledMatrix(); // deep copy
@@ -187,7 +183,6 @@ 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"
                            );
@@ -246,7 +241,6 @@ 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"
                      );
@@ -257,7 +251,6 @@ 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; }
@@ -295,18 +288,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 (sdim == 3) { z_coord(i) = coord[2]; }
+      if (dim == 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 (sdim == 3) { z.reset(z_coord.ParallelProject()); }
+   if (dim == 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 (sdim == 3)
+   if (dim == 3)
    {
       *z -= *batched_lor.GetZCoordinate();
       REQUIRE(z->Normlinf() == MFEM_Approx(0.0));

tests/unit/linalg/test_direct_solvers.cpp

@@ -26,9 +26,6 @@ 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)
 
@@ -106,7 +103,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)
       {
@@ -190,13 +187,13 @@ TEST_CASE("Parallel Direct Solvers", "[Parallel], [CUDA]")
 {
    int rank;
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
-   const int ne = 4;
+   const int ne = 2;
    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)
       {
@@ -315,39 +312,6 @@ 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]")
 {