Added dual-transfer functionality

use temp mem in tmop+pa

config/defaults.mk

@@ -138,6 +138,7 @@ MFEM_USE_RAJA          = NO
 MFEM_USE_OCCA          = NO
 MFEM_USE_CEED          = NO
 MFEM_USE_UMPIRE        = NO
+MFEM_USE_CAMP          = NO
 MFEM_USE_SIMD          = NO
 MFEM_USE_ADIOS2        = NO
 
@@ -341,9 +342,9 @@ GSLIB_DIR = @MFEM_DIR@/../gslib/build
 GSLIB_OPT = -I$(GSLIB_DIR)/include
 GSLIB_LIB = -L$(GSLIB_DIR)/lib -lgs
 
-# CUDA library configuration (currently not needed)
+# CUDA library configuration
 CUDA_OPT =
-CUDA_LIB =
+CUDA_LIB = -lcusparse
 
 # HIP library configuration (currently not needed)
 HIP_OPT =
@@ -372,6 +373,11 @@ UMPIRE_DIR = @MFEM_DIR@/../umpire
 UMPIRE_OPT = -I$(UMPIRE_DIR)/include
 UMPIRE_LIB = -L$(UMPIRE_DIR)/lib -lumpire
 
+# CAMP library configuration
+CAMP_DIR = @MFEM_DIR@/../camp
+CAMP_OPT = -I$(CAMP_DIR)/include
+CAMP_LIB = -L$(CAMP_DIR)/lib
+
 # If YES, enable some informational messages
 VERBOSE = NO
 

fem/CMakeLists.txt

@@ -76,7 +76,9 @@ set(SRCS
   tmop_pa_h3s_c0.cpp
   tmop_pa_jp2.cpp
   tmop_pa_jp3.cpp
-  tmop_pa_jt3.cpp
+  tmop_pa_jt2_tc.cpp
+  tmop_pa_jt3_datc.cpp
+  tmop_pa_jt3_tc.cpp
   tmop_pa_p2.cpp
   tmop_pa_p2_c0.cpp
   tmop_pa_p3.cpp

fem/bilininteg_diffusion_pa.cpp

@@ -740,6 +740,7 @@ static void PADiffusionAssembleDiagonal(const int dim,
          case 0x23: return SmemPADiffusionDiagonal3D<2,3>(NE,B,G,D,Y);
          case 0x34: return SmemPADiffusionDiagonal3D<3,4>(NE,B,G,D,Y);
          case 0x45: return SmemPADiffusionDiagonal3D<4,5>(NE,B,G,D,Y);
+         case 0x46: return SmemPADiffusionDiagonal3D<4,6>(NE,B,G,D,Y);
          case 0x56: return SmemPADiffusionDiagonal3D<5,6>(NE,B,G,D,Y);
          case 0x67: return SmemPADiffusionDiagonal3D<6,7>(NE,B,G,D,Y);
          case 0x78: return SmemPADiffusionDiagonal3D<7,8>(NE,B,G,D,Y);
@@ -1700,6 +1701,7 @@ static void PADiffusionApply(const int dim,
       switch (ID)
       {
          case 0x23: return SmemPADiffusionApply3D<2,3>(NE,B,G,D,X,Y);
+         case 0x24: return SmemPADiffusionApply3D<2,4>(NE,B,G,D,X,Y);
          case 0x34: return SmemPADiffusionApply3D<3,4>(NE,B,G,D,X,Y);
          case 0x45: return SmemPADiffusionApply3D<4,5>(NE,B,G,D,X,Y);
          case 0x46: return SmemPADiffusionApply3D<4,6>(NE,B,G,D,X,Y);

fem/bilininteg_mass_pa.cpp

@@ -25,6 +25,7 @@ namespace mfem
 
 void MassIntegrator::SetupPA(const FiniteElementSpace &fes)
 {
+
    // Assuming the same element type
    fespace = &fes;
    Mesh *mesh = fes.GetMesh();
@@ -53,16 +54,24 @@ void MassIntegrator::SetupPA(const FiniteElementSpace &fes)
    dofs1D = maps->ndof;
    quad1D = maps->nqpt;
    pa_data.SetSize(ne*nq, Device::GetDeviceMemoryType());
-   Vector coeff;
+   Vector *coeff{nullptr};
+   bool own_coeff{true};
    if (Q == nullptr)
    {
-      coeff.SetSize(1);
-      coeff(0) = 1.0;
+      coeff = new Vector;
+      coeff->SetSize(1);
+      (*coeff)(0) = 1.0;
    }
    else if (ConstantCoefficient* cQ = dynamic_cast<ConstantCoefficient*>(Q))
    {
-      coeff.SetSize(1);
-      coeff(0) = cQ->constant;
+      coeff = new Vector;
+      coeff->SetSize(1);
+      (*coeff)(0) = cQ->constant;
+   }
+   else if (QuadratureCoefficient* cQ = dynamic_cast<QuadratureCoefficient*>(Q))
+   {
+      coeff = cQ->Data();
+      own_coeff = false;
    }
    else if (QuadratureFunctionCoefficient* cQ =
                dynamic_cast<QuadratureFunctionCoefficient*>(Q))
@@ -75,12 +84,13 @@ void MassIntegrator::SetupPA(const FiniteElementSpace &fes)
                   "IntegrationRule used within integrator and in"
                   " QuadratureFunction appear to be different");
       qFun.Read();
-      coeff.MakeRef(const_cast<QuadratureFunction &>(qFun),0);
+      coeff->MakeRef(const_cast<QuadratureFunction &>(qFun),0);
    }
    else
    {
-      coeff.SetSize(nq * ne);
-      auto C = Reshape(coeff.HostWrite(), nq, ne);
+      coeff = new Vector;
+      coeff->SetSize(nq * ne);
+      auto C = Reshape(coeff->HostWrite(), nq, ne);
       for (int e = 0; e < ne; ++e)
       {
          ElementTransformation& T = *fes.GetElementTransformation(e);
@@ -95,11 +105,11 @@ void MassIntegrator::SetupPA(const FiniteElementSpace &fes)
    {
       const int NE = ne;
       const int NQ = nq;
-      const bool const_c = coeff.Size() == 1;
+      const bool const_c = coeff->Size() == 1;
       auto w = ir->GetWeights().Read();
       auto J = Reshape(geom->J.Read(), NQ,2,2,NE);
       auto C =
-         const_c ? Reshape(coeff.Read(), 1,1) : Reshape(coeff.Read(), NQ,NE);
+         const_c ? Reshape(coeff->Read(), 1,1) : Reshape(coeff->Read(), NQ,NE);
       auto v = Reshape(pa_data.Write(), NQ, NE);
       MFEM_FORALL(e, NE,
       {
@@ -119,12 +129,12 @@ void MassIntegrator::SetupPA(const FiniteElementSpace &fes)
    {
       const int NE = ne;
       const int Q1D = quad1D;
-      const bool const_c = coeff.Size() == 1;
+      const bool const_c = coeff->Size() == 1;
       const auto W = Reshape(ir->GetWeights().Read(),Q1D,Q1D,Q1D);
       const auto J = Reshape(geom->J.Read(), Q1D,Q1D,Q1D,3,3,NE);
       const auto C = const_c ?
-                     Reshape(coeff.Read(), 1,1,1,1) :
-                     Reshape(coeff.Read(), Q1D,Q1D,Q1D,NE);
+                     Reshape(coeff->Read(), 1,1,1,1) :
+                     Reshape(coeff->Read(), Q1D,Q1D,Q1D,NE);
       auto V = Reshape(pa_data.Write(), Q1D,Q1D,Q1D,NE);
       MFEM_FORALL_3D(e, NE, Q1D, Q1D, Q1D,
       {
@@ -153,6 +163,8 @@ void MassIntegrator::SetupPA(const FiniteElementSpace &fes)
          }
       });
    }
+
+   if (own_coeff) { delete coeff; }
 }
 
 void MassIntegrator::AssemblePA(const FiniteElementSpace &fes)

fem/coefficient.cpp

@@ -12,6 +12,7 @@
 // Implementation of Coefficient class
 
 #include "fem.hpp"
+#include "../linalg/dtensor.hpp"
 
 #include <cmath>
 #include <limits>
@@ -21,6 +22,13 @@ namespace mfem
 
 using namespace std;
 
+double QuadratureCoefficient::Eval(ElementTransformation & T,
+                                   const IntegrationPoint & ip)
+{
+   auto coeff = mfem::Reshape(qData->HostRead(), nip, NE);
+   return coeff(ip.index, T.ElementNo);
+}
+
 double PWConstCoefficient::Eval(ElementTransformation & T,
                                 const IntegrationPoint & ip)
 {

fem/coefficient.hpp

@@ -87,6 +87,33 @@ public:
    { return (constant); }
 };
 
+
+/// class for quadrature coefficient
+class QuadratureCoefficient : public Coefficient
+{
+
+private:
+   const int nip;
+   const int NE;
+public:
+   Vector *qData{nullptr};
+
+   //Set external data
+   QuadratureCoefficient(Vector *Data, int in_nip, int in_NE)
+      : qData(Data), nip(in_nip), NE(in_NE)
+   { }
+
+   virtual double Eval(ElementTransformation &T,
+                       const IntegrationPoint &ip);
+
+   Vector *Data()
+   {
+      return qData;
+   }
+
+};
+
+/// class for piecewise constant coefficient
 /** @brief A piecewise constant coefficient with the constants keyed
     off the element attribute numbers. */
 class PWConstCoefficient : public Coefficient

fem/fespace.cpp

@@ -2865,6 +2865,8 @@ void L2ProjectionGridTransfer::L2Projection::Mult(
    {
       fes_ho.GetElementVDofs(iho, vdofs);
       x.GetSubVector(vdofs, xel_mat.GetData());
+
+      // SRW: R lives on L2ProjectionGridTransfer class 
       mfem::Mult(R(iho), xel_mat, yel_mat);
       // Place result correctly into the low-order vector
       for (int iref=0; iref<nref; ++iref)
@@ -2880,6 +2882,44 @@ void L2ProjectionGridTransfer::L2Projection::Mult(
    }
 }
 
+// TODO SRW: rewrite this routine; address the offsets member variable 
+// that we don't have access to
+void L2ProjectionGridTransfer::L2Projection::MultTranspose(
+   const Vector &x, Vector &y) const
+{
+   int vdim = fes_ho.GetVDim();
+   Array<int> vdofs;
+   DenseMatrix xel_mat, yel_mat;
+   y = 0.0;
+   for (int iho = 0; iho < fes_ho.GetNE(); ++iho)
+   {
+      int nref = ho2lor.RowSize(iho);
+      int ndof_ho = fes_ho.GetFE(iho)->GetDof();
+      int ndof_lor = fes_lor.GetFE(ho2lor.GetRow(iho)[0])->GetDof();
+      xel_mat.SetSize(ndof_lor*nref, vdim);
+      yel_mat.SetSize(ndof_ho, vdim);
+//      DenseMatrix R_iho(&R[offsets[iho]], ndof_lor*nref, ndof_ho);
+
+      // Extract the LOR DOFs
+      for (int iref=0; iref<nref; ++iref)
+      {
+         int ilor = ho2lor.GetRow(iho)[iref];
+         for (int vd=0; vd<vdim; ++vd)
+         {
+            fes_lor.GetElementDofs(ilor, vdofs);
+            fes_lor.DofsToVDofs(vd, vdofs);
+            x.GetSubVector(vdofs, &xel_mat(iref*ndof_lor, vd));
+         }
+      }
+      // Multiply locally by the transpose
+//      mfem::MultAtB(R_iho, xel_mat, yel_mat);
+      mfem::MultAtB(R(iho), xel_mat, yel_mat);
+      // Place the result in the HO vector
+      fes_ho.GetElementVDofs(iho, vdofs);
+      y.AddElementVector(vdofs, yel_mat.GetData());
+   }
+}
+
 void L2ProjectionGridTransfer::L2Projection::Prolongate(
    const Vector &x, Vector &y) const
 {
@@ -2901,6 +2941,8 @@ void L2ProjectionGridTransfer::L2Projection::Prolongate(
          }
       }
       // Locally prolongate
+
+      // SRW P lives on the L2ProjectionGridTransfer class
       mfem::Mult(P(iho), xel_mat, yel_mat);
       // Place the result in the HO vector
       fes_ho.GetElementVDofs(iho, vdofs);
@@ -2908,6 +2950,41 @@ void L2ProjectionGridTransfer::L2Projection::Prolongate(
    }
 }
 
+// TODO SRW: rewrite this routine; address the offsets member variable 
+// that we don't have access to
+void L2ProjectionGridTransfer::L2Projection::ProlongateTranspose(
+   const Vector &x, Vector &y) const
+{
+   int vdim = fes_ho.GetVDim();
+   Array<int> vdofs;
+   DenseMatrix xel_mat,yel_mat;
+   for (int iho = 0; iho < fes_ho.GetNE(); ++iho)
+   {
+      int nref = ho2lor.RowSize(iho);
+      int ndof_ho = fes_ho.GetFE(iho)->GetDof();
+      int ndof_lor = fes_lor.GetFE(ho2lor.GetRow(iho)[0])->GetDof();
+      xel_mat.SetSize(ndof_ho, vdim);
+      yel_mat.SetSize(ndof_lor*nref, vdim);
+//      DenseMatrix P_iho(&P[offsets[iho]], ndof_ho, ndof_lor*nref);
+
+      fes_ho.GetElementVDofs(iho, vdofs);
+      x.GetSubVector(vdofs, xel_mat.GetData());
+      mfem::MultAtB(P(iho), xel_mat, yel_mat);
+
+      // Place result correctly into the low-order vector
+      for (int iref = 0; iref < nref; ++iref)
+      {
+         int ilor = ho2lor.GetRow(iho)[iref];
+         for (int vd=0; vd<vdim; ++vd)
+         {
+            fes_lor.GetElementDofs(ilor, vdofs);
+            fes_lor.DofsToVDofs(vd, vdofs);
+            y.SetSubVector(vdofs, &yel_mat(iref*ndof_lor,vd));
+         }
+      }
+   }
+}
+
 const Operator &L2ProjectionGridTransfer::ForwardOperator()
 {
    if (!F) { F = new L2Projection(dom_fes, ran_fes); }

fem/fespace.hpp

@@ -943,9 +943,15 @@ protected:
                    const FiniteElementSpace &fes_lor_);
       /// Perform the L2 projection onto the LOR space
       virtual void Mult(const Vector &x, Vector &y) const;
+      /// Perform the transpose of L2 projection onto the LOR space, useful for
+      /// transferring dual fields.
+      virtual void MultTranspose(const Vector &x, Vector &y) const;
       /// Perform the mass conservative left-inverse prolongation operation.
       /// This functionality is also provided as an Operator by L2Prolongation.
       void Prolongate(const Vector &x, Vector &y) const;
+      /// Perform the transpose of the mass conservative left-inverse
+      /// prolongation operation, useful for transferring dual fields.
+      void ProlongateTranspose(const Vector &x, Vector &y) const;
       virtual ~L2Projection() { }
    };
 

fem/kernels.hpp

@@ -117,7 +117,7 @@ MFEM_HOST_DEVICE inline void LoadBGt(const int D1D, const int Q1D,
 
 /// Load 2D input scalar into shared memory
 template<int MD1, int NBZ>
-MFEM_HOST_DEVICE inline void LoadS(const int e, const int D1D,
+MFEM_HOST_DEVICE inline void LoadX(const int e, const int D1D,
                                    const DeviceTensor<3, const double> x,
                                    double sX[NBZ][MD1*MD1])
 {
@@ -593,6 +593,27 @@ MFEM_HOST_DEVICE inline void GradXt(const int D1D, const int Q1D,
 
 /// Load 3D scalar input vector into shared memory
 template<int MD1>
+MFEM_HOST_DEVICE inline void LoadX(const int e, const int D1D,
+                                   const DeviceTensor<4, const double> x,
+                                   double sm[MD1*MD1*MD1])
+{
+   DeviceCube X(sm, MD1, MD1, MD1);
+
+   MFEM_FOREACH_THREAD(dz,z,D1D)
+   {
+      MFEM_FOREACH_THREAD(dy,y,D1D)
+      {
+         MFEM_FOREACH_THREAD(dx,x,D1D)
+         {
+            X(dx,dy,dz) = x(dx,dy,dz,e);
+         }
+      }
+   }
+   MFEM_SYNC_THREAD;
+}
+
+/// Load 3D scalar input vector into shared memory, with comp
+template<int MD1>
 MFEM_HOST_DEVICE inline void LoadX(const int e, const int D1D, const int c,
                                    const DeviceTensor<5, const double> x,
                                    double sm[MD1*MD1*MD1])

fem/quadinterpolator_face.cpp

@@ -495,8 +495,9 @@ void FaceQuadratureInterpolator::Mult(
    }
 }
 
-void FaceQuadratureInterpolator::Values(
-   const Vector &e_vec, Vector &q_val) const
+
+void FaceQuadratureInterpolator::Values(const Vector &e_vec,
+                                        Vector &q_val) const
 {
    Vector q_der, q_det, q_nor;
    Mult(e_vec, VALUES, q_val, q_der, q_det, q_nor);

fem/tmop.hpp

@@ -681,6 +681,9 @@ public:
        nodes are used by all target types except IDEAL_SHAPE_UNIT_SIZE. */
    void SetNodes(const GridFunction &n) { nodes = &n; avg_volume = 0.0; }
 
+   /** @brief Get the nodes to be used in the target-matrix construction. */
+   const GridFunction *GetNodes() const { return nodes; }
+
    /// Used by target type IDEAL_SHAPE_EQUAL_SIZE. The default volume scale is 1.
    void SetVolumeScale(double vol_scale) { volume_scale = vol_scale; }
 
@@ -697,7 +700,13 @@ public:
                                       const Vector &elfun,
                                       DenseTensor &Jtr) const;
 
-   virtual bool ComputeElementTargetsPA(const IntegrationRule *ir,
+   template<int DIM>
+   bool ComputeElementTargetsPA(const FiniteElementSpace *fes,
+                                const IntegrationRule *ir,
+                                DenseTensor &Jtr,
+                                const Vector &xe = Vector()) const;
+   virtual bool ComputeElementTargetsPA(const FiniteElementSpace *fes,
+                                        const IntegrationRule *ir,
                                         DenseTensor &Jtr,
                                         const Vector &xe = Vector()) const;
 
@@ -746,7 +755,8 @@ public:
                                       const Vector &elfun,
                                       DenseTensor &Jtr) const;
 
-   virtual bool ComputeElementTargetsPA(const IntegrationRule *ir,
+   virtual bool ComputeElementTargetsPA(const FiniteElementSpace *fes,
+                                        const IntegrationRule *ir,
                                         DenseTensor &Jtr,
                                         const Vector &xe = Vector()) const;
 
@@ -901,7 +911,8 @@ public:
                                       const Vector &elfun,
                                       DenseTensor &Jtr) const;
 
-   virtual bool ComputeElementTargetsPA(const IntegrationRule *ir,
+   virtual bool ComputeElementTargetsPA(const FiniteElementSpace *fes,
+                                        const IntegrationRule *ir,
                                         DenseTensor &Jtr,
                                         const Vector &xe = Vector()) const;
 
@@ -989,9 +1000,11 @@ protected:
       mutable bool setup_Grad, setup_Jtr;
       mutable Vector E, O, W, X0, H, C0, LD, H0;
       const DofToQuad *maps;
+      const DofToQuad *maps_lim = nullptr;
       const GeometricFactors *geom;
       const FiniteElementSpace *fes;
       const Operator *R;
+      const IntegrationRule *ir;
    } PA;
 
    void ComputeNormalizationEnergies(const GridFunction &x,
@@ -1196,6 +1209,15 @@ public:
 
    /** @brief Flag to control if exact action of Integration is effected. */
    void SetExactActionFlag(bool flag_) { exact_action = flag_; }
+
+   void ReleaseTemporaryMemory()
+   {
+      if (PA.enabled)
+      {
+         PA.H.GetMemory().DeleteDevice();
+         PA.H0.GetMemory().DeleteDevice();
+      }
+   }
 };
 
 class TMOPComboIntegrator : public NonlinearFormIntegrator

fem/tmop_pa.cpp

@@ -43,20 +43,29 @@ void TMOP_Integrator::SetupGradPA(const Vector &xe) const
 // It is the case when EnableLimiting is called before the Setup => AssemblePA.
 void TMOP_Integrator::EnableLimitingPA(const GridFunction &n0)
 {
-   MFEM_ASSERT(PA.enabled, "EnableLimitingPA but PA is not enabled!");
+   MFEM_VERIFY(PA.enabled, "EnableLimitingPA but PA is not enabled!");
    const ElementDofOrdering ordering = ElementDofOrdering::LEXICOGRAPHIC;
-   PA.R = n0.FESpace()->GetElementRestriction(ordering);
 
    // Nodes0
-   PA.X0.SetSize(PA.R->Height(), Device::GetMemoryType());
+   const FiniteElementSpace *n0_fes = n0.FESpace();
+   const Operator *n0_R = n0_fes->GetElementRestriction(ordering);
+   PA.X0.SetSize(n0_R->Height(), Device::GetMemoryType());
    PA.X0.UseDevice(true);
-   PA.R->Mult(n0, PA.X0);
+   n0_R->Mult(n0, PA.X0);
+
+   // Get the 1D maps for the distance FE space.
+   const IntegrationRule &ir = *EnergyIntegrationRule(*n0.FESpace()->GetFE(0));
+   PA.maps_lim =
+      &lim_dist->FESpace()->GetFE(0)->GetDofToQuad(ir, DofToQuad::TENSOR);
 
    // lim_dist & lim_func checks
    MFEM_VERIFY(lim_dist, "No lim_dist!")
-   PA.LD.SetSize(PA.R->Height(), Device::GetMemoryType());
+   const FiniteElementSpace *ld_fes = lim_dist->FESpace();
+   const Operator *ld_R = ld_fes->GetElementRestriction(ordering);
+   MFEM_VERIFY(ld_R, "No lim_dist restriction operator found!");
+   PA.LD.SetSize(ld_R->Height(), Device::GetMemoryType());
    PA.LD.UseDevice(true);
-   PA.R->Mult(*lim_dist, PA.LD);
+   ld_R->Mult(*lim_dist, PA.LD);
 
    // Only TMOP_QuadraticLimiter is supported
    MFEM_VERIFY(lim_func, "No lim_func!")
@@ -64,6 +73,27 @@ void TMOP_Integrator::EnableLimitingPA(const GridFunction &n0)
                "Only TMOP_QuadraticLimiter is supported");
 }
 
+bool TargetConstructor::ComputeElementTargetsPA(const FiniteElementSpace *fes,
+                                                const IntegrationRule *ir,
+                                                DenseTensor &Jtr,
+                                                const Vector &xe) const
+{
+   MFEM_VERIFY(Jtr.SizeI() == Jtr.SizeJ() && Jtr.SizeI() > 1, "");
+   const int dim = Jtr.SizeI();
+   if (dim == 2) { return ComputeElementTargetsPA<2>(fes, ir, Jtr, xe); }
+   if (dim == 3) { return ComputeElementTargetsPA<3>(fes, ir, Jtr, xe); }
+   return false;
+}
+
+bool AnalyticAdaptTC::ComputeElementTargetsPA(const FiniteElementSpace *fes,
+                                              const IntegrationRule *ir,
+                                              DenseTensor &Jtr,
+                                              const Vector &xe) const
+{
+   return false;
+}
+
+
 // Code paths leading to ComputeElementTargets:
 // - GetElementEnergy(elfun) which is done through GetGridFunctionEnergyPA(x)
 // - AssembleElementVectorExact(elfun)
@@ -89,12 +119,19 @@ void TMOP_Integrator::ComputeElementTargetsPA(const Vector &xe) const
    PA.setup_Jtr = false;
    const FiniteElementSpace *fes = PA.fes;
    const IntegrationRule *ir = EnergyIntegrationRule(*fes->GetFE(0));
+   const TargetConstructor::TargetType &target_type = targetC->Type();
+   const DiscreteAdaptTC *discr_tc = GetDiscreteAdaptTC();
 
-   {
-      // Try to use the TargetConstructor ComputeElementTargetsPA
-      PA.setup_Jtr = targetC->ComputeElementTargetsPA(ir,PA.Jtr);
-      if (PA.setup_Jtr) { return; }
-   }
+   // Skip when TargetConstructor needs the nodes but have not been set
+   const bool use_nodes =
+      target_type == TargetConstructor::IDEAL_SHAPE_EQUAL_SIZE ||
+      target_type == TargetConstructor::IDEAL_SHAPE_GIVEN_SIZE ||
+      target_type == TargetConstructor::GIVEN_SHAPE_AND_SIZE;
+   if (targetC && !discr_tc && use_nodes && !targetC->GetNodes()) { return; }
+
+   // Try to use the TargetConstructor ComputeElementTargetsPA
+   PA.setup_Jtr = targetC->ComputeElementTargetsPA(fes, ir, PA.Jtr);
+   if (PA.setup_Jtr) { return; }
 
    // Defaulting to host version
    PA.Jtr.HostWrite();
@@ -104,15 +141,12 @@ void TMOP_Integrator::ComputeElementTargetsPA(const Vector &xe) const
    const int dim = PA.dim;
    DenseTensor &Jtr = PA.Jtr;
 
-   const TargetConstructor::TargetType &target_type = targetC->Type();
-
    Vector x;
    const bool useable_input_vector = xe.Size() > 0;
    const bool use_input_vector = target_type == TargetConstructor::GIVEN_FULL;
 
    if (use_input_vector && !useable_input_vector) { return; }
 
-   DiscreteAdaptTC *discr_tc = GetDiscreteAdaptTC();
    if (discr_tc && !discr_tc->GetTspecFesv()) { return; }
 
    if (use_input_vector)
@@ -148,7 +182,9 @@ void TMOP_Integrator::ComputeElementTargetsPA(const Vector &xe) const
 void TMOP_Integrator::AssemblePA(const FiniteElementSpace &fes)
 {
    PA.enabled = true;
-   const IntegrationRule *ir = EnergyIntegrationRule(*fes.GetFE(0));
+   MFEM_ASSERT(fes.GetMesh()->GetNE() > 0, "");
+   PA.ir = EnergyIntegrationRule(*fes.GetFE(0));
+   const IntegrationRule *ir = PA.ir;
    MFEM_ASSERT(fes.GetOrdering() == Ordering::byNODES,
                "PA Only supports Ordering::byNODES!");
 
@@ -171,12 +207,18 @@ void TMOP_Integrator::AssemblePA(const FiniteElementSpace &fes)
    PA.setup_Jtr = false;
    PA.setup_Grad = false;
 
+
+#ifdef MFEM_USE_UMPIRE
+   const MemoryType temp_type = Device::GetDeviceMemoryType() == MemoryType::DEVICE_UMPIRE
+      ? MemoryType::DEVICE_UMPIRE_2 : Device::GetDeviceMemoryType();
+#else
+   const MemoryType temp_type = Device::GetDeviceMemoryType();
+#endif
+
    // H for Grad
-   PA.H.UseDevice(true);
-   PA.H.SetSize(dim*dim * dim*dim * nq*ne, Device::GetDeviceMemoryType());
+   PA.H.SetSize(dim*dim * dim*dim * nq*ne, temp_type);
    // H0 for coeff0
-   PA.H0.UseDevice(true);
-   PA.H0.SetSize(dim * dim * nq*ne, Device::GetDeviceMemoryType());
+   PA.H0.SetSize(dim * dim * nq*ne, temp_type);
 
    // Restriction setup
    const ElementDofOrdering ordering = ElementDofOrdering::LEXICOGRAPHIC;

fem/tmop_pa_h2d_c0.cpp

@@ -76,6 +76,7 @@ MFEM_REGISTER_TMOP_KERNELS(void, AssembleDiagonalPA_Kernel_C0_2D,
                D(dx,dy,v,e) += d;
             }
          }
+         MFEM_SYNC_THREAD;
       }
    });
 }

fem/tmop_pa_h2s.cpp

@@ -236,8 +236,7 @@ void TMOP_Integrator::AssembleGradPA_2D(const Vector &X) const
    const int id = (D1D << 4 ) | Q1D;
    const double mn = metric_normal;
    const DenseTensor &J = PA.Jtr;
-   const IntegrationRule *ir = IntRule;
-   const Array<double> &W = ir->GetWeights();
+   const Array<double> &W = PA.ir->GetWeights();
    const Array<double> &B = PA.maps->B;
    const Array<double> &G = PA.maps->G;
    Vector &H = PA.H;

fem/tmop_pa_h2s_c0.cpp

@@ -26,6 +26,7 @@ MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_C0_2D,
                            const DenseTensor &j_,
                            const Array<double> &w_,
                            const Array<double> &b_,
+                           const Array<double> &bld_,
                            Vector &h0_,
                            const int d1d,
                            const int q1d)
@@ -39,10 +40,11 @@ MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_C0_2D,
    const auto C0 = const_c0 ?
                    Reshape(c0_.Read(), 1, 1, 1) :
                    Reshape(c0_.Read(), Q1D, Q1D, NE);
-   const auto LD = Reshape(lim_dist.Read(), D1D, D1D, DIM, NE);
+   const auto LD = Reshape(lim_dist.Read(), D1D, D1D, NE);
    const auto J = Reshape(j_.Read(), DIM, DIM, Q1D, Q1D, NE);
    const auto W = Reshape(w_.Read(), Q1D, Q1D);
    const auto b = Reshape(b_.Read(), Q1D, D1D);
+   const auto bld = Reshape(bld_.Read(), Q1D, D1D);
 
    auto H0 = Reshape(h0_.Write(), DIM, DIM, Q1D, Q1D, NE);
 
@@ -55,15 +57,19 @@ MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_C0_2D,
       constexpr int MD1 = T_D1D ? T_D1D : T_MAX;
 
       MFEM_SHARED double B[MQ1*MD1];
+      MFEM_SHARED double BLD[MQ1*MD1];
 
-      MFEM_SHARED double XY[2][NBZ][MD1*MD1];
-      MFEM_SHARED double DQ[2][NBZ][MD1*MQ1];
-      MFEM_SHARED double QQ[2][NBZ][MQ1*MQ1];
+      MFEM_SHARED double XY[NBZ][MD1*MD1];
+      MFEM_SHARED double DQ[NBZ][MD1*MQ1];
+      MFEM_SHARED double QQ[NBZ][MQ1*MQ1];
 
       kernels::LoadX<MD1,NBZ>(e,D1D,LD,XY);
+
       kernels::LoadB<MD1,MQ1>(D1D,Q1D,b,B);
-      kernels::EvalX<MD1,MQ1,NBZ>(D1D,Q1D,B,XY,DQ);
-      kernels::EvalY<MD1,MQ1,NBZ>(D1D,Q1D,B,DQ,QQ);
+      kernels::LoadB<MD1,MQ1>(D1D,Q1D,bld,BLD);
+
+      kernels::EvalX<MD1,MQ1,NBZ>(D1D,Q1D,BLD,XY,DQ);
+      kernels::EvalY<MD1,MQ1,NBZ>(D1D,Q1D,BLD,DQ,QQ);
 
       MFEM_FOREACH_THREAD(qy,y,Q1D)
       {
@@ -75,9 +81,9 @@ MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_C0_2D,
             const double coeff0 = const_c0 ? C0(0,0,0) : C0(qx,qy,e);
             const double weight_m = weight * lim_normal * coeff0;
 
-            double D[2];
+            double D;
             kernels::PullEval<MQ1,NBZ>(qx,qy,QQ,D);
-            const double dist = D[0]; // GetValues, default comp set to 0
+            const double dist = D; // GetValues, default comp set to 0
 
             // lim_func->Eval_d2(p1, p0, d_vals(q), grad_grad);
             // d2.Diag(1.0 / (dist * dist), x.Size());
@@ -107,13 +113,13 @@ void TMOP_Integrator::AssembleGradPA_C0_2D(const Vector &X) const
    const double ln = lim_normal;
    const Vector &LD = PA.LD;
    const DenseTensor &J = PA.Jtr;
-   const IntegrationRule *ir = IntRule;
-   const Array<double> &W = ir->GetWeights();
-   const Array<double> &B = PA.maps->B;
+   const Array<double> &W   = PA.ir->GetWeights();
+   const Array<double> &B   = PA.maps->B;
+   const Array<double> &BLD = PA.maps_lim->B;
    const Vector &C0 = PA.C0;
    Vector &H0 = PA.H0;
 
-   MFEM_LAUNCH_TMOP_KERNEL(SetupGradPA_C0_2D,id,ln,LD,C0,N,J,W,B,H0);
+   MFEM_LAUNCH_TMOP_KERNEL(SetupGradPA_C0_2D,id,ln,LD,C0,N,J,W,B,BLD,H0);
 }
 
 } // namespace mfem

fem/tmop_pa_h3s.cpp

@@ -270,8 +270,7 @@ void TMOP_Integrator::AssembleGradPA_3D(const Vector &X) const
    const int id = (D1D << 4 ) | Q1D;
    const double mn = metric_normal;
    const DenseTensor &J = PA.Jtr;
-   const IntegrationRule *ir = IntRule;
-   const Array<double> &W = ir->GetWeights();
+   const Array<double> &W = PA.ir->GetWeights();
    const Array<double> &B = PA.maps->B;
    const Array<double> &G = PA.maps->G;
    Vector &H = PA.H;

fem/tmop_pa_h3s_c0.cpp

@@ -26,6 +26,7 @@ MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_Kernel_C0_3D,
                            const DenseTensor &j_,
                            const Array<double> &w_,
                            const Array<double> &b_,
+                           const Array<double> &bld_,
                            Vector &h0_,
                            const int d1d,
                            const int q1d)
@@ -38,10 +39,11 @@ MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_Kernel_C0_3D,
    const auto C0 = const_c0 ?
                    Reshape(c0_.Read(), 1, 1, 1, 1) :
                    Reshape(c0_.Read(), Q1D, Q1D, Q1D, NE);
-   const auto LD = Reshape(lim_dist.Read(), D1D, D1D, D1D, DIM, NE);
+   const auto LD = Reshape(lim_dist.Read(), D1D, D1D, D1D, NE);
    const auto J = Reshape(j_.Read(), DIM, DIM, Q1D, Q1D, Q1D, NE);
    const auto W = Reshape(w_.Read(), Q1D, Q1D, Q1D);
    const auto b = Reshape(b_.Read(), Q1D, D1D);
+   const auto bld = Reshape(bld_.Read(), Q1D, D1D);
 
    auto H0 = Reshape(h0_.Write(), DIM, DIM, Q1D, Q1D, Q1D, NE);
 
@@ -54,16 +56,21 @@ MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_Kernel_C0_3D,
       constexpr int MD1 = T_D1D ? T_D1D : T_MAX;
 
       MFEM_SHARED double B[MQ1*MD1];
-      MFEM_SHARED double DDD[3][MD1*MD1*MD1];
-      MFEM_SHARED double DDQ[3][MD1*MD1*MQ1];
-      MFEM_SHARED double DQQ[3][MD1*MQ1*MQ1];
-      MFEM_SHARED double QQQ[3][MQ1*MQ1*MQ1];
+      MFEM_SHARED double BLD[MQ1*MD1];
+
+      MFEM_SHARED double DDD[MD1*MD1*MD1];
+      MFEM_SHARED double DDQ[MD1*MD1*MQ1];
+      MFEM_SHARED double DQQ[MD1*MQ1*MQ1];
+      MFEM_SHARED double QQQ[MQ1*MQ1*MQ1];
 
       kernels::LoadX<MD1>(e,D1D,LD,DDD);
+
       kernels::LoadB<MD1,MQ1>(D1D,Q1D,b,B);
-      kernels::EvalX<MD1,MQ1>(D1D,Q1D,B,DDD,DDQ);
-      kernels::EvalY<MD1,MQ1>(D1D,Q1D,B,DDQ,DQQ);
-      kernels::EvalZ<MD1,MQ1>(D1D,Q1D,B,DQQ,QQQ);
+      kernels::LoadB<MD1,MQ1>(D1D,Q1D,bld,BLD);
+
+      kernels::EvalX<MD1,MQ1>(D1D,Q1D,BLD,DDD,DDQ);
+      kernels::EvalY<MD1,MQ1>(D1D,Q1D,BLD,DDQ,DQQ);
+      kernels::EvalZ<MD1,MQ1>(D1D,Q1D,BLD,DQQ,QQQ);
 
       MFEM_FOREACH_THREAD(qz,z,Q1D)
       {
@@ -77,9 +84,9 @@ MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_Kernel_C0_3D,
                const double coeff0 = const_c0 ? C0(0,0,0,0) : C0(qx,qy,qz,e);
                const double weight_m = weight * lim_normal * coeff0;
 
-               double D[3];
+               double D;
                kernels::PullEval<MQ1>(qx,qy,qz,QQQ,D);
-               const double dist = D[0]; // GetValues, default comp set to 0
+               const double dist = D; // GetValues, default comp set to 0
 
                // lim_func->Eval_d2(p1, p0, d_vals(q), grad_grad);
                // d2.Diag(1.0 / (dist * dist), x.Size());
@@ -110,13 +117,13 @@ void TMOP_Integrator::AssembleGradPA_C0_3D(const Vector &X) const
    const double ln = lim_normal;
    const Vector &LD = PA.LD;
    const DenseTensor &J = PA.Jtr;
-   const IntegrationRule *ir = IntRule;
-   const Array<double> &W = ir->GetWeights();
+   const Array<double> &W = PA.ir->GetWeights();
    const Array<double> &B = PA.maps->B;
+   const Array<double> &BLD = PA.maps_lim->B;
    const Vector &C0 = PA.C0;
    Vector &H0 = PA.H0;
 
-   MFEM_LAUNCH_TMOP_KERNEL(SetupGradPA_Kernel_C0_3D,id,ln,LD,C0,N,J,W,B,H0);
+   MFEM_LAUNCH_TMOP_KERNEL(SetupGradPA_Kernel_C0_3D,id,ln,LD,C0,N,J,W,B,BLD,H0);
 }
 
 } // namespace mfem

fem/tmop_pa_jt2_tc.cpp

@@ -0,0 +1,155 @@
+// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "tmop.hpp"
+#include "tmop_pa.hpp"
+#include "gridfunc.hpp"
+#include "../general/forall.hpp"
+#include "../linalg/kernels.hpp"
+
+namespace mfem
+{
+
+MFEM_REGISTER_TMOP_KERNELS(bool, TC_IDEAL_SHAPE_UNIT_SIZE_2D_KERNEL,
+                           const int NE,
+                           const DenseMatrix w_, // copy
+                           DenseTensor &j_,
+                           const int d1d,
+                           const int q1d)
+{
+   constexpr int DIM = 2;
+   constexpr int NBZ = 1;
+
+   const int Q1D = T_Q1D ? T_Q1D : q1d;
+
+   const auto W = Reshape(w_.Read(), DIM,DIM);
+   auto J = Reshape(j_.Write(), DIM,DIM, Q1D,Q1D, NE);
+
+   MFEM_FORALL_2D(e, NE, Q1D, Q1D, NBZ,
+   {
+      const int Q1D = T_Q1D ? T_Q1D : q1d;
+      MFEM_FOREACH_THREAD(qy,y,Q1D)
+      {
+         MFEM_FOREACH_THREAD(qx,x,Q1D)
+         {
+            kernels::Set(DIM,DIM, 1.0, &W(0,0), &J(0,0,qx,qy,e));
+         }
+      }
+   });
+   return true;
+}
+
+MFEM_REGISTER_TMOP_KERNELS(bool, TC_IDEAL_SHAPE_GIVEN_SIZE_2D_KERNEL,
+                           const int NE,
+                           const Array<double> &b_,
+                           const Array<double> &g_,
+                           const DenseMatrix w_ideal_, // copy
+                           const Vector &x_,
+                           DenseTensor &j_,
+                           const int d1d,
+                           const int q1d)
+{
+   constexpr int DIM = 2;
+   constexpr int NBZ = 1;
+
+   const double detW = w_ideal_.Det();
+   const int D1D = T_D1D ? T_D1D : d1d;
+   const int Q1D = T_Q1D ? T_Q1D : q1d;
+
+   const auto b = Reshape(b_.Read(), Q1D, D1D);
+   const auto g = Reshape(g_.Read(), Q1D, D1D);
+   const auto Wideal = Reshape(w_ideal_.Read(), DIM,DIM);
+   const auto X = Reshape(x_.Read(), D1D, D1D, DIM, NE);
+   auto J = Reshape(j_.Write(), DIM,DIM, Q1D,Q1D, NE);
+
+   MFEM_FORALL_2D(e, NE, Q1D, Q1D, NBZ,
+   {
+      const int D1D = T_D1D ? T_D1D : d1d;
+      const int Q1D = T_Q1D ? T_Q1D : q1d;
+
+      constexpr int NBZ = 1;
+      constexpr int MQ1 = T_Q1D ? T_Q1D : T_MAX;
+      constexpr int MD1 = T_D1D ? T_D1D : T_MAX;
+
+      MFEM_SHARED double BG[2][MQ1*MD1];
+      MFEM_SHARED double XY[2][NBZ][MD1*MD1];
+      MFEM_SHARED double DQ[4][NBZ][MD1*MQ1];
+      MFEM_SHARED double QQ[4][NBZ][MQ1*MQ1];
+
+      kernels::LoadX<MD1,NBZ>(e,D1D,X,XY);
+      kernels::LoadBG<MD1,MQ1>(D1D,Q1D,b,g,BG);
+
+      kernels::GradX<MD1,MQ1,NBZ>(D1D,Q1D,BG,XY,DQ);
+      kernels::GradY<MD1,MQ1,NBZ>(D1D,Q1D,BG,DQ,QQ);
+
+      MFEM_FOREACH_THREAD(qy,y,Q1D)
+      {
+         MFEM_FOREACH_THREAD(qx,x,Q1D)
+         {
+            double Jtr[4];
+            const double *Wid = &Wideal(0,0);
+            kernels::PullGrad<MQ1,NBZ>(qx,qy,QQ,Jtr);
+            const double detJ = kernels::Det<2>(Jtr);
+            const double alpha = std::pow(detJ/detW,1./2);
+            kernels::Set(DIM,DIM,alpha,Wid,&J(0,0,qx,qy,e));
+         }
+      }
+   });
+   return true;
+}
+
+template<> bool
+TargetConstructor::ComputeElementTargetsPA<2>(const FiniteElementSpace *fes,
+                                              const IntegrationRule *ir,
+                                              DenseTensor &Jtr,
+                                              const Vector&) const
+{
+   MFEM_ASSERT(target_type == IDEAL_SHAPE_UNIT_SIZE || nodes != nullptr, "");
+   MFEM_VERIFY(fes->GetFE(0)->GetGeomType() == Geometry::SQUARE, "");
+   const DenseMatrix &W = Geometries.GetGeomToPerfGeomJac(Geometry::SQUARE);
+   const FiniteElement *fe = fes->GetFE(0);
+   const int NE = fes->GetMesh()->GetNE();
+   const DofToQuad::Mode mode = DofToQuad::TENSOR;
+   const DofToQuad &maps = fe->GetDofToQuad(*ir, mode);
+   const int D1D = maps.ndof;
+   const int Q1D = maps.nqpt;
+   const int id = (D1D << 4 ) | Q1D;
+
+   const Array<double> &B = maps.B;
+   const Array<double> &G = maps.G;
+
+   switch (target_type)
+   {
+      case IDEAL_SHAPE_UNIT_SIZE: // Jtr(i) = Wideal;
+      {
+         MFEM_LAUNCH_TMOP_KERNEL(TC_IDEAL_SHAPE_UNIT_SIZE_2D_KERNEL,
+                                 id,NE,W,Jtr);
+      }
+      case IDEAL_SHAPE_EQUAL_SIZE: return false;
+      case IDEAL_SHAPE_GIVEN_SIZE:
+      {
+         MFEM_VERIFY(nodes, "");
+         const ElementDofOrdering ordering = ElementDofOrdering::LEXICOGRAPHIC;
+         const Operator *R = fes->GetElementRestriction(ordering);
+         Vector X(R->Height(), Device::GetDeviceMemoryType());
+         X.UseDevice(true);
+         R->Mult(*nodes, X);
+         MFEM_ASSERT(nodes->FESpace()->GetVDim() == 2, "");
+         MFEM_LAUNCH_TMOP_KERNEL(TC_IDEAL_SHAPE_GIVEN_SIZE_2D_KERNEL,
+                                 id,NE,B,G,W,X,Jtr);
+      }
+      case GIVEN_SHAPE_AND_SIZE: return false;
+      default: return false;
+   }
+   return false;
+}
+
+} // namespace mfem

fem/tmop_pa_jt3_datc.cpp

@@ -17,20 +17,6 @@
 namespace mfem
 {
 
-bool TargetConstructor::ComputeElementTargetsPA(const IntegrationRule *ir,
-                                                DenseTensor &Jtr,
-                                                const Vector &xe) const
-{
-   return false;
-}
-
-bool AnalyticAdaptTC::ComputeElementTargetsPA(const IntegrationRule *ir,
-                                              DenseTensor &Jtr,
-                                              const Vector &xe) const
-{
-   return false;
-}
-
 MFEM_REGISTER_TMOP_KERNELS(bool, DatcSize,
                            const int NE,
                            const int ncomp,
@@ -127,7 +113,8 @@ MFEM_REGISTER_TMOP_KERNELS(bool, DatcSize,
 }
 
 // PA.Jtr Size = (dim, dim, PA.ne*PA.nq);
-bool DiscreteAdaptTC::ComputeElementTargetsPA(const IntegrationRule *ir,
+bool DiscreteAdaptTC::ComputeElementTargetsPA(const FiniteElementSpace *pa_fes,
+                                              const IntegrationRule *ir,
                                               DenseTensor &Jtr,
                                               const Vector &xe) const
 {

fem/tmop_pa_jt3_tc.cpp

@@ -0,0 +1,159 @@
+// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "tmop.hpp"
+#include "tmop_pa.hpp"
+#include "gridfunc.hpp"
+#include "../general/forall.hpp"
+#include "../linalg/kernels.hpp"
+
+namespace mfem
+{
+
+MFEM_REGISTER_TMOP_KERNELS(bool, TC_IDEAL_SHAPE_UNIT_SIZE_3D_KERNEL,
+                           const int NE,
+                           const DenseMatrix w_, // copy
+                           DenseTensor &j_,
+                           const int d1d,
+                           const int q1d)
+{
+   constexpr int DIM = 3;
+
+   const int Q1D = T_Q1D ? T_Q1D : q1d;
+
+   const auto W = Reshape(w_.Read(), DIM,DIM);
+   auto J = Reshape(j_.Write(), DIM,DIM, Q1D,Q1D,Q1D, NE);
+
+   MFEM_FORALL_3D(e, NE, Q1D, Q1D, Q1D,
+   {
+      const int Q1D = T_Q1D ? T_Q1D : q1d;
+      MFEM_FOREACH_THREAD(qy,y,Q1D)
+      {
+         MFEM_FOREACH_THREAD(qx,x,Q1D)
+         {
+            MFEM_FOREACH_THREAD(qz,z,Q1D)
+            {
+               kernels::Set(DIM,DIM, 1.0, &W(0,0), &J(0,0,qx,qy,qz,e));
+            }
+         }
+      }
+   });
+   return true;
+}
+
+MFEM_REGISTER_TMOP_KERNELS(bool, TC_IDEAL_SHAPE_GIVEN_SIZE_3D_KERNEL,
+                           const int NE,
+                           const Array<double> &b_,
+                           const Array<double> &g_,
+                           const DenseMatrix w_ideal_, // copy
+                           const Vector &x_,
+                           DenseTensor &j_,
+                           const int d1d,
+                           const int q1d)
+{
+   constexpr int DIM = 3;
+
+   const double detW = w_ideal_.Det();
+   const int D1D = T_D1D ? T_D1D : d1d;
+   const int Q1D = T_Q1D ? T_Q1D : q1d;
+
+   const auto b = Reshape(b_.Read(), Q1D, D1D);
+   const auto g = Reshape(g_.Read(), Q1D, D1D);
+   const auto Wideal = Reshape(w_ideal_.Read(), DIM,DIM);
+   const auto X = Reshape(x_.Read(), D1D, D1D, D1D, DIM, NE);
+   auto J = Reshape(j_.Write(), DIM,DIM, Q1D,Q1D,Q1D, NE);
+
+   MFEM_FORALL_3D(e, NE, Q1D, Q1D, Q1D,
+   {
+      const int D1D = T_D1D ? T_D1D : d1d;
+      const int Q1D = T_Q1D ? T_Q1D : q1d;
+
+      constexpr int MQ1 = T_Q1D ? T_Q1D : T_MAX;
+      constexpr int MD1 = T_D1D ? T_D1D : T_MAX;
+
+      MFEM_SHARED double BG[2][MQ1*MD1];
+      MFEM_SHARED double DDD[3][MD1*MD1*MD1];
+      MFEM_SHARED double DDQ[6][MD1*MD1*MQ1];
+      MFEM_SHARED double DQQ[9][MD1*MQ1*MQ1];
+      MFEM_SHARED double QQQ[9][MQ1*MQ1*MQ1];
+
+      kernels::LoadX<MD1>(e,D1D,X,DDD);
+      kernels::LoadBG<MD1,MQ1>(D1D,Q1D,b,g,BG);
+
+      kernels::GradX<MD1,MQ1>(D1D,Q1D,BG,DDD,DDQ);
+      kernels::GradY<MD1,MQ1>(D1D,Q1D,BG,DDQ,DQQ);
+      kernels::GradZ<MD1,MQ1>(D1D,Q1D,BG,DQQ,QQQ);
+
+      MFEM_FOREACH_THREAD(qz,z,Q1D)
+      {
+         MFEM_FOREACH_THREAD(qy,y,Q1D)
+         {
+            MFEM_FOREACH_THREAD(qx,x,Q1D)
+            {
+               double Jtr[9];
+               const double *Wid = &Wideal(0,0);
+               kernels::PullGrad<MQ1>(qx,qy,qz,QQQ,Jtr);
+               const double detJ = kernels::Det<3>(Jtr);
+               const double alpha = std::pow(detJ/detW,1./3);
+               kernels::Set(DIM,DIM,alpha,Wid,&J(0,0,qx,qy,qz,e));
+            }
+         }
+      }
+   });
+   return true;
+}
+
+template<> bool
+TargetConstructor::ComputeElementTargetsPA<3>(const FiniteElementSpace *fes,
+                                              const IntegrationRule *ir,
+                                              DenseTensor &Jtr,
+                                              const Vector&) const
+{
+   MFEM_ASSERT(target_type == IDEAL_SHAPE_UNIT_SIZE || nodes != nullptr, "");
+   MFEM_VERIFY(fes->GetFE(0)->GetGeomType() == Geometry::CUBE, "");
+   const DenseMatrix &W = Geometries.GetGeomToPerfGeomJac(Geometry::CUBE);
+   const FiniteElement *fe = fes->GetFE(0);
+   const int NE = fes->GetMesh()->GetNE();
+   const DofToQuad::Mode mode = DofToQuad::TENSOR;
+   const DofToQuad &maps = fe->GetDofToQuad(*ir, mode);
+   const Array<double> &B = maps.B;
+   const Array<double> &G = maps.G;
+   const int D1D = maps.ndof;
+   const int Q1D = maps.nqpt;
+   const int id = (D1D << 4 ) | Q1D;
+
+   switch (target_type)
+   {
+      case IDEAL_SHAPE_UNIT_SIZE: // Jtr(i) = Wideal;
+      {
+         MFEM_LAUNCH_TMOP_KERNEL(TC_IDEAL_SHAPE_UNIT_SIZE_3D_KERNEL,
+                                 id,NE,W,Jtr);
+      }
+      case IDEAL_SHAPE_EQUAL_SIZE: return false;
+      case IDEAL_SHAPE_GIVEN_SIZE:
+      {
+         MFEM_VERIFY(nodes, "");
+         const ElementDofOrdering ordering = ElementDofOrdering::LEXICOGRAPHIC;
+         const Operator *R = fes->GetElementRestriction(ordering);
+         Vector X(R->Height(), Device::GetDeviceMemoryType());
+         X.UseDevice(true);
+         R->Mult(*nodes, X);
+         MFEM_ASSERT(nodes->FESpace()->GetVDim() == 3, "");
+         MFEM_LAUNCH_TMOP_KERNEL(TC_IDEAL_SHAPE_GIVEN_SIZE_3D_KERNEL,
+                                 id,NE,B,G,W,X,Jtr);
+      }
+      case GIVEN_SHAPE_AND_SIZE: return false;
+      default: return false;
+   }
+   return false;
+}
+
+} // namespace mfem

fem/tmop_pa_p2.cpp

@@ -155,7 +155,7 @@ void TMOP_Integrator::AddMultPA_2D(const Vector &X, Vector &Y) const
    const int Q1D = PA.maps->nqpt;
    const int id = (D1D << 4 ) | Q1D;
    const DenseTensor &J = PA.Jtr;
-   const Array<double> &W = IntRule->GetWeights();
+   const Array<double> &W = PA.ir->GetWeights();
    const Array<double> &B = PA.maps->B;
    const Array<double> &G = PA.maps->G;
    const double mn = metric_normal;

fem/tmop_pa_p2_c0.cpp

@@ -26,6 +26,7 @@ MFEM_REGISTER_TMOP_KERNELS(void, AddMultPA_Kernel_C0_2D,
                            const DenseTensor &j_,
                            const Array<double> &w_,
                            const Array<double> &b_,
+                           const Array<double> &bld_,
                            const Vector &x0_,
                            const Vector &x1_,
                            Vector &y_,
@@ -43,9 +44,10 @@ MFEM_REGISTER_TMOP_KERNELS(void, AddMultPA_Kernel_C0_2D,
    const auto C0 = const_c0 ?
                    Reshape(c0_.Read(), 1, 1, 1) :
                    Reshape(c0_.Read(), Q1D, Q1D, NE);
-   const auto LD = Reshape(lim_dist.Read(), D1D, D1D, DIM, NE);
+   const auto LD = Reshape(lim_dist.Read(), D1D, D1D, NE);
    const auto J = Reshape(j_.Read(), DIM, DIM, Q1D, Q1D, NE);
    const auto b = Reshape(b_.Read(), Q1D, D1D);
+   const auto bld = Reshape(bld_.Read(), Q1D, D1D);
    const auto W = Reshape(w_.Read(), Q1D, Q1D);
    const auto X0 = Reshape(x0_.Read(), D1D, D1D, DIM, NE);
    const auto X1 = Reshape(x1_.Read(), D1D, D1D, DIM, NE);
@@ -61,10 +63,11 @@ MFEM_REGISTER_TMOP_KERNELS(void, AddMultPA_Kernel_C0_2D,
       constexpr int MD1 = T_D1D ? T_D1D : T_MAX;
 
       MFEM_SHARED double B[MQ1*MD1];
+      MFEM_SHARED double BLD[MQ1*MD1];
 
-      MFEM_SHARED double XY[2][NBZ][MD1*MD1];
-      MFEM_SHARED double DQ[2][NBZ][MD1*MQ1];
-      MFEM_SHARED double QQ[2][NBZ][MQ1*MQ1];
+      MFEM_SHARED double XY[NBZ][MD1*MD1];
+      MFEM_SHARED double DQ[NBZ][MD1*MQ1];
+      MFEM_SHARED double QQ[NBZ][MQ1*MQ1];
 
       MFEM_SHARED double XY0[2][NBZ][MD1*MD1];
       MFEM_SHARED double DQ0[2][NBZ][MD1*MQ1];
@@ -79,9 +82,10 @@ MFEM_REGISTER_TMOP_KERNELS(void, AddMultPA_Kernel_C0_2D,
       kernels::LoadX<MD1,NBZ>(e,D1D,X1,XY1);
 
       kernels::LoadB<MD1,MQ1>(D1D,Q1D,b,B);
+      kernels::LoadB<MD1,MQ1>(D1D,Q1D,bld,BLD);
 
-      kernels::EvalX<MD1,MQ1,NBZ>(D1D,Q1D,B,XY,DQ);
-      kernels::EvalY<MD1,MQ1,NBZ>(D1D,Q1D,B,DQ,QQ);
+      kernels::EvalX<MD1,MQ1,NBZ>(D1D,Q1D,BLD,XY,DQ);
+      kernels::EvalY<MD1,MQ1,NBZ>(D1D,Q1D,BLD,DQ,QQ);
 
       kernels::EvalX<MD1,MQ1,NBZ>(D1D,Q1D,B,XY0,DQ0);
       kernels::EvalY<MD1,MQ1,NBZ>(D1D,Q1D,B,DQ0,QQ0);
@@ -97,13 +101,13 @@ MFEM_REGISTER_TMOP_KERNELS(void, AddMultPA_Kernel_C0_2D,
             const double detJtr = kernels::Det<2>(Jtr);
             const double weight = W(qx,qy) * detJtr;
 
-            double ld[2], p0[2], p1[2];
+            double ld, p0[2], p1[2];
             const double coeff0 = const_c0 ? C0(0,0,0) : C0(qx,qy,e);
             kernels::PullEval<MQ1,NBZ>(qx,qy,QQ,ld);
             kernels::PullEval<MQ1,NBZ>(qx,qy,QQ0,p0);
             kernels::PullEval<MQ1,NBZ>(qx,qy,QQ1,p1);
 
-            const double dist = ld[0]; // GetValues, default comp set to 0
+            const double dist = ld; // GetValues, default comp set to 0
 
             double d1[2];
             // Eval_d1
@@ -132,12 +136,13 @@ void TMOP_Integrator::AddMultPA_C0_2D(const Vector &X, Vector &Y) const
    const double ln = lim_normal;
    const Vector &LD = PA.LD;
    const DenseTensor &J = PA.Jtr;
-   const Array<double> &W = IntRule->GetWeights();
-   const Array<double> &B = PA.maps->B;
+   const Array<double> &W   = PA.ir->GetWeights();
+   const Array<double> &B   = PA.maps->B;
+   const Array<double> &BLD = PA.maps_lim->B;
    const Vector &X0 = PA.X0;
    const Vector &C0 = PA.C0;
 
-   MFEM_LAUNCH_TMOP_KERNEL(AddMultPA_Kernel_C0_2D,id,ln,LD,C0,N,J,W,B,X0,X,Y);
+   MFEM_LAUNCH_TMOP_KERNEL(AddMultPA_Kernel_C0_2D,id,ln,LD,C0,N,J,W,B,BLD,X0,X,Y);
 }
 
 } // namespace mfem

fem/tmop_pa_p3.cpp

@@ -175,7 +175,7 @@ void TMOP_Integrator::AddMultPA_3D(const Vector &X, Vector &Y) const
    const int Q1D = PA.maps->nqpt;
    const int id = (D1D << 4 ) | Q1D;
    const DenseTensor &J = PA.Jtr;
-   const Array<double> &W = IntRule->GetWeights();
+   const Array<double> &W = PA.ir->GetWeights();
    const Array<double> &B = PA.maps->B;
    const Array<double> &G = PA.maps->G;
    const double mn = metric_normal;

fem/tmop_pa_p3_c0.cpp

@@ -26,6 +26,7 @@ MFEM_REGISTER_TMOP_KERNELS(void, AddMultPA_Kernel_C0_3D,
                            const DenseTensor &j_,
                            const Array<double> &w_,
                            const Array<double> &b_,
+                           const Array<double> &bld_,
                            const Vector &x0_,
                            const Vector &x1_,
                            Vector &y_,
@@ -41,9 +42,10 @@ MFEM_REGISTER_TMOP_KERNELS(void, AddMultPA_Kernel_C0_3D,
    const auto C0 = const_c0 ?
                    Reshape(c0_.Read(), 1, 1, 1, 1) :
                    Reshape(c0_.Read(), Q1D, Q1D, Q1D, NE);
-   const auto LD = Reshape(lim_dist.Read(), D1D, D1D, D1D, DIM, NE);
+   const auto LD = Reshape(lim_dist.Read(), D1D, D1D, D1D, NE);
    const auto J = Reshape(j_.Read(), DIM, DIM, Q1D, Q1D, Q1D, NE);
    const auto b = Reshape(b_.Read(), Q1D, D1D);
+   const auto bld = Reshape(bld_.Read(), Q1D, D1D);
    const auto W = Reshape(w_.Read(), Q1D, Q1D, Q1D);
    const auto X0 = Reshape(x0_.Read(), D1D, D1D, D1D, DIM, NE);
    const auto X1 = Reshape(x1_.Read(), D1D, D1D, D1D, DIM, NE);
@@ -58,11 +60,12 @@ MFEM_REGISTER_TMOP_KERNELS(void, AddMultPA_Kernel_C0_3D,
       constexpr int MD1 = T_D1D ? T_D1D : T_MAX;
 
       MFEM_SHARED double B[MQ1*MD1];
+      MFEM_SHARED double BLD[MQ1*MD1];
 
-      MFEM_SHARED double DDD[3][MD1*MD1*MD1];
-      MFEM_SHARED double DDQ[3][MD1*MD1*MQ1];
-      MFEM_SHARED double DQQ[3][MD1*MQ1*MQ1];
-      MFEM_SHARED double QQQ[3][MQ1*MQ1*MQ1];
+      MFEM_SHARED double DDD[MD1*MD1*MD1];
+      MFEM_SHARED double DDQ[MD1*MD1*MQ1];
+      MFEM_SHARED double DQQ[MD1*MQ1*MQ1];
+      MFEM_SHARED double QQQ[MQ1*MQ1*MQ1];
 
       MFEM_SHARED double DDD0[3][MD1*MD1*MD1];
       MFEM_SHARED double DDQ0[3][MD1*MD1*MQ1];
@@ -79,10 +82,11 @@ MFEM_REGISTER_TMOP_KERNELS(void, AddMultPA_Kernel_C0_3D,
       kernels::LoadX<MD1>(e,D1D,X1,DDD1);
 
       kernels::LoadB<MD1,MQ1>(D1D,Q1D,b,B);
+      kernels::LoadB<MD1,MQ1>(D1D,Q1D,bld,BLD);
 
-      kernels::EvalX<MD1,MQ1>(D1D,Q1D,B,DDD,DDQ);
-      kernels::EvalY<MD1,MQ1>(D1D,Q1D,B,DDQ,DQQ);
-      kernels::EvalZ<MD1,MQ1>(D1D,Q1D,B,DQQ,QQQ);
+      kernels::EvalX<MD1,MQ1>(D1D,Q1D,BLD,DDD,DDQ);
+      kernels::EvalY<MD1,MQ1>(D1D,Q1D,BLD,DDQ,DQQ);
+      kernels::EvalZ<MD1,MQ1>(D1D,Q1D,BLD,DQQ,QQQ);
 
       kernels::EvalX<MD1,MQ1>(D1D,Q1D,B,DDD0,DDQ0);
       kernels::EvalY<MD1,MQ1>(D1D,Q1D,B,DDQ0,DQQ0);
@@ -102,7 +106,7 @@ MFEM_REGISTER_TMOP_KERNELS(void, AddMultPA_Kernel_C0_3D,
                const double detJtr = kernels::Det<3>(Jtr);
                const double weight = W(qx,qy,qz) * detJtr;
 
-               double D[3], p0[3], p1[3];
+               double D, p0[3], p1[3];
                const double coeff0 = const_c0 ? C0(0,0,0,0) : C0(qx,qy,qz,e);
                kernels::PullEval<MQ1>(qx,qy,qz,QQQ,D);
                kernels::PullEval<MQ1>(qx,qy,qz,QQQ0,p0);
@@ -113,7 +117,7 @@ MFEM_REGISTER_TMOP_KERNELS(void, AddMultPA_Kernel_C0_3D,
                // subtract(1.0 / (dist * dist), x, x0, d1);
                // z = a * (x - y)
                // grad = a * (x - x0)
-               const double dist = D[0]; // GetValues, default comp set to 0
+               const double dist = D; // GetValues, default comp set to 0
                const double a = 1.0 / (dist * dist);
                const double w = weight * lim_normal * coeff0;
                kernels::Subtract<3>(w*a, p1, p0, d1);
@@ -138,12 +142,13 @@ void TMOP_Integrator::AddMultPA_C0_3D(const Vector &X, Vector &Y) const
    const double ln = lim_normal;
    const Vector &LD = PA.LD;
    const DenseTensor &J = PA.Jtr;
-   const Array<double> &W = IntRule->GetWeights();
+   const Array<double> &W = PA.ir->GetWeights();
    const Array<double> &B = PA.maps->B;
+   const Array<double> &BLD = PA.maps_lim->B;
    const Vector &X0 = PA.X0;
    const Vector &C0 = PA.C0;
 
-   MFEM_LAUNCH_TMOP_KERNEL(AddMultPA_Kernel_C0_3D,id,ln,LD,C0,N,J,W,B,X0,X,Y);
+   MFEM_LAUNCH_TMOP_KERNEL(AddMultPA_Kernel_C0_3D,id,ln,LD,C0,N,J,W,B,BLD,X0,X,Y);
 }
 
 } // namespace mfem

fem/tmop_pa_w2.cpp

@@ -143,8 +143,7 @@ double TMOP_Integrator::GetGridFunctionEnergyPA_2D(const Vector &X) const
    const int id = (D1D << 4 ) | Q1D;
    const double m = metric_normal;
    const DenseTensor &J = PA.Jtr;
-   const IntegrationRule *ir = IntRule;
-   const Array<double> &W = ir->GetWeights();
+   const Array<double> &W = PA.ir->GetWeights();
    const Array<double> &B = PA.maps->B;
    const Array<double> &G = PA.maps->G;
    const Vector &O = PA.O;

fem/tmop_pa_w2_c0.cpp

@@ -26,6 +26,7 @@ MFEM_REGISTER_TMOP_KERNELS(double, EnergyPA_C0_2D,
                            const DenseTensor &j_,
                            const Array<double> &w_,
                            const Array<double> &b_,
+                           const Array<double> &bld_,
                            const Vector &x0_,
                            const Vector &x1_,
                            const Vector &ones,
@@ -44,9 +45,10 @@ MFEM_REGISTER_TMOP_KERNELS(double, EnergyPA_C0_2D,
    const auto C0 = const_c0 ?
                    Reshape(c0_.Read(), 1, 1, 1) :
                    Reshape(c0_.Read(), Q1D, Q1D, NE);
-   const auto LD = Reshape(lim_dist.Read(), D1D, D1D, DIM, NE);
+   const auto LD = Reshape(lim_dist.Read(), D1D, D1D, NE);
    const auto J = Reshape(j_.Read(), DIM, DIM, Q1D, Q1D, NE);
    const auto b = Reshape(b_.Read(), Q1D, D1D);
+   const auto bld = Reshape(bld_.Read(), Q1D, D1D);
    const auto W = Reshape(w_.Read(), Q1D, Q1D);
    const auto X0 = Reshape(x0_.Read(), D1D, D1D, DIM, NE);
    const auto X1 = Reshape(x1_.Read(), D1D, D1D, DIM, NE);
@@ -62,10 +64,11 @@ MFEM_REGISTER_TMOP_KERNELS(double, EnergyPA_C0_2D,
       constexpr int MD1 = T_D1D ? T_D1D : T_MAX;
 
       MFEM_SHARED double B[MQ1*MD1];
+      MFEM_SHARED double BLD[MQ1*MD1];
 
-      MFEM_SHARED double XY[2][NBZ][MD1*MD1];
-      MFEM_SHARED double DQ[2][NBZ][MD1*MQ1];
-      MFEM_SHARED double QQ[2][NBZ][MQ1*MQ1];
+      MFEM_SHARED double XY[NBZ][MD1*MD1];
+      MFEM_SHARED double DQ[NBZ][MD1*MQ1];
+      MFEM_SHARED double QQ[NBZ][MQ1*MQ1];
 
       MFEM_SHARED double XY0[2][NBZ][MD1*MD1];
       MFEM_SHARED double DQ0[2][NBZ][MD1*MQ1];
@@ -80,9 +83,10 @@ MFEM_REGISTER_TMOP_KERNELS(double, EnergyPA_C0_2D,
       kernels::LoadX<MD1,NBZ>(e,D1D,X1,XY1);
 
       kernels::LoadB<MD1,MQ1>(D1D,Q1D,b,B);
+      kernels::LoadB<MD1,MQ1>(D1D,Q1D,bld,BLD);
 
-      kernels::EvalX<MD1,MQ1,NBZ>(D1D,Q1D,B,XY,DQ);
-      kernels::EvalY<MD1,MQ1,NBZ>(D1D,Q1D,B,DQ,QQ);
+      kernels::EvalX<MD1,MQ1,NBZ>(D1D,Q1D,BLD,XY,DQ);
+      kernels::EvalY<MD1,MQ1,NBZ>(D1D,Q1D,BLD,DQ,QQ);
 
       kernels::EvalX<MD1,MQ1,NBZ>(D1D,Q1D,B,XY0,DQ0);
       kernels::EvalY<MD1,MQ1,NBZ>(D1D,Q1D,B,DQ0,QQ0);
@@ -94,7 +98,7 @@ MFEM_REGISTER_TMOP_KERNELS(double, EnergyPA_C0_2D,
       {
          MFEM_FOREACH_THREAD(qx,x,Q1D)
          {
-            double ld[2], p0[2], p1[2];
+            double ld, p0[2], p1[2];
             const double *Jtr = &J(0,0,qx,qy,e);
             const double detJtr = kernels::Det<2>(Jtr);
             const double weight = W(qx,qy) * detJtr;
@@ -102,7 +106,7 @@ MFEM_REGISTER_TMOP_KERNELS(double, EnergyPA_C0_2D,
             kernels::PullEval<MQ1,NBZ>(qx,qy,QQ,ld);
             kernels::PullEval<MQ1,NBZ>(qx,qy,QQ0,p0);
             kernels::PullEval<MQ1,NBZ>(qx,qy,QQ1,p1);
-            const double dist = ld[0]; // GetValues, default comp set to 0
+            const double dist = ld; // GetValues, default comp set to 0
             const double id2 = 0.5 / (dist*dist);
             const double dsq = kernels::DistanceSquared<2>(p1,p0) * id2;
             E(qx,qy,e) = weight * lim_normal * dsq * coeff0;
@@ -121,15 +125,15 @@ double TMOP_Integrator::GetGridFunctionEnergyPA_C0_2D(const Vector &X) const
    const double ln = lim_normal;
    const Vector &LD = PA.LD;
    const DenseTensor &J = PA.Jtr;
-   const IntegrationRule *ir = IntRule;
-   const Array<double> &W = ir->GetWeights();
-   const Array<double> &B = PA.maps->B;
+   const Array<double> &W   = PA.ir->GetWeights();
+   const Array<double> &B   = PA.maps->B;
+   const Array<double> &BLD = PA.maps_lim->B;
    const Vector &X0 = PA.X0;
    const Vector &C0 = PA.C0;
    const Vector &O = PA.O;
    Vector &E = PA.E;
 
-   MFEM_LAUNCH_TMOP_KERNEL(EnergyPA_C0_2D,id,ln,LD,C0,N,J,W,B,X0,X,O,E);
+   MFEM_LAUNCH_TMOP_KERNEL(EnergyPA_C0_2D,id,ln,LD,C0,N,J,W,B,BLD,X0,X,O,E);
 }
 
 } // namespace mfem

fem/tmop_pa_w3.cpp

@@ -153,7 +153,7 @@ double TMOP_Integrator::GetGridFunctionEnergyPA_3D(const Vector &X) const
    const int id = (D1D << 4 ) | Q1D;
    const double mn = metric_normal;
    const DenseTensor &J = PA.Jtr;
-   const Array<double> &W = IntRule->GetWeights();
+   const Array<double> &W = PA.ir->GetWeights();
    const Array<double> &B = PA.maps->B;
    const Array<double> &G = PA.maps->G;
    const Vector &O = PA.O;

fem/tmop_pa_w3_c0.cpp

@@ -26,6 +26,7 @@ MFEM_REGISTER_TMOP_KERNELS(double, EnergyPA_C0_3D,
                            const DenseTensor &j_,
                            const Array<double> &w_,
                            const Array<double> &b_,
+                           const Array<double> &bld_,
                            const Vector &x0_,
                            const Vector &x1_,
                            const Vector &ones,
@@ -42,9 +43,10 @@ MFEM_REGISTER_TMOP_KERNELS(double, EnergyPA_C0_3D,
    const auto C0 = const_c0 ?
                    Reshape(c0_.Read(), 1, 1, 1, 1) :
                    Reshape(c0_.Read(), Q1D, Q1D, Q1D, NE);
-   const auto LD = Reshape(lim_dist.Read(), D1D, D1D, D1D, DIM, NE);
+   const auto LD = Reshape(lim_dist.Read(), D1D, D1D, D1D, NE);
    const auto J = Reshape(j_.Read(), DIM, DIM, Q1D, Q1D, Q1D, NE);
    const auto b = Reshape(b_.Read(), Q1D, D1D);
+   const auto bld = Reshape(bld_.Read(), Q1D, D1D);
    const auto W = Reshape(w_.Read(), Q1D, Q1D, Q1D);
    const auto X0 = Reshape(x0_.Read(), D1D, D1D, D1D, DIM, NE);
    const auto X1 = Reshape(x1_.Read(), D1D, D1D, D1D, DIM, NE);
@@ -59,11 +61,12 @@ MFEM_REGISTER_TMOP_KERNELS(double, EnergyPA_C0_3D,
       constexpr int MD1 = T_D1D ? T_D1D : T_MAX;
 
       MFEM_SHARED double B[MQ1*MD1];
+      MFEM_SHARED double BLD[MQ1*MD1];
 
-      MFEM_SHARED double DDD[3][MD1*MD1*MD1];
-      MFEM_SHARED double DDQ[3][MD1*MD1*MQ1];
-      MFEM_SHARED double DQQ[3][MD1*MQ1*MQ1];
-      MFEM_SHARED double QQQ[3][MQ1*MQ1*MQ1];
+      MFEM_SHARED double DDD[MD1*MD1*MD1];
+      MFEM_SHARED double DDQ[MD1*MD1*MQ1];
+      MFEM_SHARED double DQQ[MD1*MQ1*MQ1];
+      MFEM_SHARED double QQQ[MQ1*MQ1*MQ1];
 
       MFEM_SHARED double DDD0[3][MD1*MD1*MD1];
       MFEM_SHARED double DDQ0[3][MD1*MD1*MQ1];
@@ -80,10 +83,11 @@ MFEM_REGISTER_TMOP_KERNELS(double, EnergyPA_C0_3D,
       kernels::LoadX<MD1>(e,D1D,X1,DDD1);
 
       kernels::LoadB<MD1,MQ1>(D1D,Q1D,b,B);
+      kernels::LoadB<MD1,MQ1>(D1D,Q1D,bld,BLD);
 
-      kernels::EvalX<MD1,MQ1>(D1D,Q1D,B,DDD,DDQ);
-      kernels::EvalY<MD1,MQ1>(D1D,Q1D,B,DDQ,DQQ);
-      kernels::EvalZ<MD1,MQ1>(D1D,Q1D,B,DQQ,QQQ);
+      kernels::EvalX<MD1,MQ1>(D1D,Q1D,BLD,DDD,DDQ);
+      kernels::EvalY<MD1,MQ1>(D1D,Q1D,BLD,DDQ,DQQ);
+      kernels::EvalZ<MD1,MQ1>(D1D,Q1D,BLD,DQQ,QQQ);
 
       kernels::EvalX<MD1,MQ1>(D1D,Q1D,B,DDD0,DDQ0);
       kernels::EvalY<MD1,MQ1>(D1D,Q1D,B,DDQ0,DQQ0);
@@ -99,7 +103,7 @@ MFEM_REGISTER_TMOP_KERNELS(double, EnergyPA_C0_3D,
          {
             MFEM_FOREACH_THREAD(qx,x,Q1D)
             {
-               double D[3], p0[3], p1[3];
+               double D, p0[3], p1[3];
                const double *Jtr = &J(0,0,qx,qy,qz,e);
                const double detJtr = kernels::Det<3>(Jtr);
                const double weight = W(qx,qy,qz) * detJtr;
@@ -109,7 +113,7 @@ MFEM_REGISTER_TMOP_KERNELS(double, EnergyPA_C0_3D,
                kernels::PullEval<MQ1>(qx,qy,qz,QQQ0,p0);
                kernels::PullEval<MQ1>(qx,qy,qz,QQQ1,p1);
 
-               const double dist = D[0]; // GetValues, default comp set to 0
+               const double dist = D; // GetValues, default comp set to 0
                const double id2 = 0.5 / (dist*dist);
 
                const double dsq = kernels::DistanceSquared<3>(p1,p0) * id2;
@@ -130,14 +134,15 @@ double TMOP_Integrator::GetGridFunctionEnergyPA_C0_3D(const Vector &X) const
    const double ln = lim_normal;
    const Vector &LD = PA.LD;
    const DenseTensor &J = PA.Jtr;
-   const Array<double> &W = IntRule->GetWeights();
+   const Array<double> &W = PA.ir->GetWeights();
    const Array<double> &B = PA.maps->B;
+   const Array<double> &BLD = PA.maps_lim->B;
    const Vector &X0 = PA.X0;
    const Vector &C0 = PA.C0;
    const Vector &O = PA.O;
    Vector &E = PA.E;
 
-   MFEM_LAUNCH_TMOP_KERNEL(EnergyPA_C0_3D,id,ln,LD,C0,N,J,W,B,X0,X,O,E);
+   MFEM_LAUNCH_TMOP_KERNEL(EnergyPA_C0_3D,id,ln,LD,C0,N,J,W,B,BLD,X0,X,O,E);
 }
 
 } // namespace mfem

fem/tmop_tools.cpp

@@ -622,6 +622,7 @@ void TMOPNewtonSolver::ProcessNewState(const Vector &x) const
             Array<TMOP_Integrator *> ati = co->GetTMOPIntegrators();
             for (int j = 0; j < ati.Size(); j++)
             {
+               ati[j]->UpdateAfterMeshChange(x_loc);
                ati[j]->ComputeFDh(x_loc, *pfesc);
                UpdateDiscreteTC(*ati[j], x_loc);
             }
@@ -658,6 +659,7 @@ void TMOPNewtonSolver::ProcessNewState(const Vector &x) const
             Array<TMOP_Integrator *> ati = co->GetTMOPIntegrators();
             for (int j = 0; j < ati.Size(); j++)
             {
+               ati[j]->UpdateAfterMeshChange(x_loc);
                ati[j]->ComputeFDh(x_loc, *fesc);
                UpdateDiscreteTC(*ati[j], x_loc);
             }

general/cuda.cpp

@@ -1,4 +1,4 @@
-// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 // LICENSE and NOTICE for details. LLNL-CODE-806117.
 //
@@ -61,6 +61,21 @@ void* CuMallocManaged(void** dptr, size_t bytes)
    return *dptr;
 }
 
+void* CuMemAllocHostPinned(void** ptr, size_t bytes)
+{
+#ifdef MFEM_USE_CUDA
+#ifdef MFEM_TRACK_CUDA_MEM
+   mfem::out << "CuMemAllocHostPinned(): allocating " << bytes << " bytes ... "
+             << std::flush;
+#endif
+   MFEM_GPU_CHECK(cudaMallocHost(ptr, bytes));
+#ifdef MFEM_TRACK_CUDA_MEM
+   mfem::out << "done: " << *ptr << std::endl;
+#endif
+#endif
+   return *ptr;
+}
+
 void* CuMemFree(void *dptr)
 {
 #ifdef MFEM_USE_CUDA
@@ -76,6 +91,21 @@ void* CuMemFree(void *dptr)
    return dptr;
 }
 
+void* CuMemFreeHostPinned(void *ptr)
+{
+#ifdef MFEM_USE_CUDA
+#ifdef MFEM_TRACK_CUDA_MEM
+   mfem::out << "CuMemFreeHostPinned(): deallocating memory @ " << ptr << " ... "
+             << std::flush;
+#endif
+   MFEM_GPU_CHECK(cudaFreeHost(ptr));
+#ifdef MFEM_TRACK_CUDA_MEM
+   mfem::out << "done." << std::endl;
+#endif
+#endif
+   return ptr;
+}
+
 void* CuMemcpyHtoD(void* dst, const void* src, size_t bytes)
 {
 #ifdef MFEM_USE_CUDA

general/cuda.hpp

@@ -1,4 +1,4 @@
-// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 // LICENSE and NOTICE for details. LLNL-CODE-806117.
 //
@@ -43,7 +43,6 @@
 #if defined(MFEM_USE_CUDA) && defined(__CUDA_ARCH__)
 #define MFEM_SHARED __shared__
 #define MFEM_SYNC_THREAD __syncthreads()
-#define MFEM_BLOCK_ID(k) blockIdx.k
 #define MFEM_THREAD_ID(k) threadIdx.k
 #define MFEM_THREAD_SIZE(k) blockDim.k
 #define MFEM_FOREACH_THREAD(i,k,N) for(int i=threadIdx.k; i<N; i+=blockDim.k)
@@ -64,9 +63,15 @@ void* CuMemAlloc(void **d_ptr, size_t bytes);
 /// Allocates managed device memory
 void* CuMallocManaged(void **d_ptr, size_t bytes);
 
+/// Allocates page-locked (pinned) host memory
+void* CuMemAllocHostPinned(void **ptr, size_t bytes);
+
 /// Frees device memory and returns destination ptr.
 void* CuMemFree(void *d_ptr);
 
+/// Frees page-locked (pinned) host memory and returns destination ptr.
+void* CuMemFreeHostPinned(void *ptr);
+
 /// Copies memory from Host to Device and returns destination ptr.
 void* CuMemcpyHtoD(void *d_dst, const void *h_src, size_t bytes);
 

general/device.cpp

@@ -1,4 +1,4 @@
-// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 // LICENSE and NOTICE for details. LLNL-CODE-806117.
 //
@@ -12,9 +12,10 @@
 #include "forall.hpp"
 #include "occa.hpp"
 #ifdef MFEM_USE_CEED
-#include <ceed.h>
+#include "../fem/ceed/util.hpp"
 #endif
 
+#include <unordered_map>
 #include <string>
 #include <map>
 
@@ -33,13 +34,16 @@ occa::device occaDevice;
 
 #ifdef MFEM_USE_CEED
 Ceed ceed = NULL;
+
+ceed::BasisMap ceed_basis_map;
+ceed::RestrMap ceed_restr_map;
 #endif
 
 // Backends listed by priority, high to low:
 static const Backend::Id backend_list[Backend::NUM_BACKENDS] =
 {
    Backend::CEED_CUDA, Backend::OCCA_CUDA, Backend::RAJA_CUDA, Backend::CUDA,
-   Backend::HIP, Backend::DEBUG,
+   Backend::CEED_HIP, Backend::RAJA_HIP, Backend::HIP, Backend::DEBUG,
    Backend::OCCA_OMP, Backend::RAJA_OMP, Backend::OMP,
    Backend::CEED_CPU, Backend::OCCA_CPU, Backend::RAJA_CPU, Backend::CPU
 };
@@ -48,7 +52,7 @@ static const Backend::Id backend_list[Backend::NUM_BACKENDS] =
 static const char *backend_name[Backend::NUM_BACKENDS] =
 {
    "ceed-cuda", "occa-cuda", "raja-cuda", "cuda",
-   "hip", "debug",
+   "ceed-hip", "raja-hip", "hip", "debug",
    "occa-omp", "raja-omp", "omp",
    "ceed-cpu", "occa-cpu", "raja-cpu", "cpu"
 };
@@ -61,15 +65,9 @@ Device Device::device_singleton;
 bool Device::device_env = false;
 bool Device::mem_host_env = false;
 bool Device::mem_device_env = false;
+bool Device::mem_types_set = false;
 
-Device::Device() : mode(Device::SEQUENTIAL),
-   backends(Backend::CPU),
-   destroy_mm(false),
-   mpi_gpu_aware(false),
-   host_mem_type(MemoryType::HOST),
-   host_mem_class(MemoryClass::HOST),
-   device_mem_type(MemoryType::HOST),
-   device_mem_class(MemoryClass::HOST)
+Device::Device()
 {
    if (getenv("MFEM_MEMORY") && !mem_host_env && !mem_device_env)
    {
@@ -154,6 +152,18 @@ Device::~Device()
    {
       free(device_option);
 #ifdef MFEM_USE_CEED
+      // Destroy FES -> CeedBasis, CeedElemRestriction hash table contents
+      for (auto entry : internal::ceed_basis_map)
+      {
+         CeedBasisDestroy(&entry.second);
+      }
+      internal::ceed_basis_map.clear();
+      for (auto entry : internal::ceed_restr_map)
+      {
+         CeedElemRestrictionDestroy(&entry.second);
+      }
+      internal::ceed_restr_map.clear();
+      // Destroy Ceed context
       CeedDestroy(&internal::ceed);
 #endif
       mm.Destroy();
@@ -210,15 +220,24 @@ void Device::Configure(const std::string &device, const int dev)
       beg = end + 1;
    }
 
-   // OCCA_CUDA needs CUDA or RAJA_CUDA:
-   if (Allows(Backend::OCCA_CUDA) && !Allows(Backend::RAJA_CUDA))
+   // OCCA_CUDA and CEED_CUDA need CUDA or RAJA_CUDA:
+   if (Allows(Backend::OCCA_CUDA|Backend::CEED_CUDA) &&
+       !Allows(Backend::RAJA_CUDA))
    {
       Get().MarkBackend(Backend::CUDA);
    }
-   if (Allows(Backend::CEED_CUDA))
+   // CEED_HIP needs HIP:
+   if (Allows(Backend::CEED_HIP))
    {
-      Get().MarkBackend(Backend::CUDA);
+      Get().MarkBackend(Backend::HIP);
    }
+   // OCCA_OMP will use OMP or RAJA_OMP unless MFEM_USE_OPENMP=NO:
+#ifdef MFEM_USE_OPENMP
+   if (Allows(Backend::OCCA_OMP) && !Allows(Backend::RAJA_OMP))
+   {
+      Get().MarkBackend(Backend::OMP);
+   }
+#endif
 
    // Perform setup.
    Get().Setup(dev);
@@ -233,6 +252,30 @@ void Device::Configure(const std::string &device, const int dev)
    destroy_mm = true;
 }
 
+// static method
+void Device::SetMemoryTypes(MemoryType h_mt, MemoryType d_mt)
+{
+   // If the device and/or the MemoryTypes are configured through the
+   // environment (variables 'MFEM_DEVICE', 'MFEM_MEMORY'), ignore calls to this
+   // method.
+   if (mem_host_env || mem_device_env || device_env) { return; }
+
+   MFEM_VERIFY(!IsConfigured(), "the default MemoryTypes can only be set before"
+               " Device construction and configuration");
+   MFEM_VERIFY(IsHostMemory(h_mt),
+               "invalid host MemoryType, h_mt = " << (int)h_mt);
+   MFEM_VERIFY(IsDeviceMemory(d_mt) || d_mt == h_mt,
+               "invalid device MemoryType, d_mt = " << (int)d_mt
+               << " (h_mt = " << (int)h_mt << ')');
+
+   Get().host_mem_type = h_mt;
+   Get().device_mem_type = d_mt;
+   mem_types_set = true;
+
+   // h_mt and d_mt will be set as dual to each other during configuration by
+   // the call mm.Configure(...) in UpdateMemoryTypeAndClass()
+}
+
 void Device::Print(std::ostream &out)
 {
    out << "Device configuration: ";
@@ -272,7 +315,14 @@ void Device::UpdateMemoryTypeAndClass()
 
 #ifdef MFEM_USE_UMPIRE
    // If MFEM has been compiled with Umpire support, use it as the default
-   if (!mem_host_env) { host_mem_type = MemoryType::HOST_UMPIRE; }
+   if (!mem_host_env && !mem_types_set)
+   {
+      host_mem_type = MemoryType::HOST_UMPIRE;
+      if (!mem_device_env)
+      {
+         device_mem_type = MemoryType::HOST_UMPIRE;
+      }
+   }
 #endif
 
    // Enable the device memory type
@@ -294,7 +344,7 @@ void Device::UpdateMemoryTypeAndClass()
                   device_mem_type = MemoryType::DEVICE;
             }
          }
-         else
+         else if (!mem_types_set)
          {
 #ifndef MFEM_USE_UMPIRE
             device_mem_type = MemoryType::DEVICE;
@@ -320,6 +370,9 @@ void Device::UpdateMemoryTypeAndClass()
       device_mem_type = MemoryType::DEVICE_DEBUG;
    }
 
+   MFEM_VERIFY(!device || IsDeviceMemory(device_mem_type),
+               "invalid device memory configuration!");
+
    // Update the memory manager with the new settings
    mm.Configure(host_mem_type, device_mem_type);
 }
@@ -369,6 +422,8 @@ static void RajaDeviceSetup(const int dev, int &ngpu)
 {
 #ifdef MFEM_USE_CUDA
    if (ngpu <= 0) { DeviceSetup(dev, ngpu); }
+#elif defined(MFEM_USE_HIP)
+   HipDeviceSetup(dev, ngpu);
 #else
    MFEM_CONTRACT_VAR(dev);
    MFEM_CONTRACT_VAR(ngpu);
@@ -435,7 +490,8 @@ static void CeedDeviceSetup(const char* ceed_spec)
    CeedInit(ceed_spec, &internal::ceed);
    const char *ceed_backend;
    CeedGetResource(internal::ceed, &ceed_backend);
-   if (strcmp(ceed_spec, ceed_backend) && strcmp(ceed_spec, "/cpu/self"))
+   if (strcmp(ceed_spec, ceed_backend) && strcmp(ceed_spec, "/cpu/self") &&
+       strcmp(ceed_spec, "/gpu/hip"))
    {
       mfem::out << std::endl << "WARNING!!!\n"
                 "libCEED is not using the requested backend!!!\n"
@@ -473,12 +529,16 @@ void Device::Setup(const int device)
    MFEM_VERIFY(!Allows(Backend::CEED_MASK),
                "the CEED backends require MFEM built with MFEM_USE_CEED=YES");
 #else
-   MFEM_VERIFY(!Allows(Backend::CEED_CPU) || !Allows(Backend::CEED_CUDA),
+   int ceed_cpu  = Allows(Backend::CEED_CPU);
+   int ceed_cuda = Allows(Backend::CEED_CUDA);
+   int ceed_hip  = Allows(Backend::CEED_HIP);
+   MFEM_VERIFY(ceed_cpu + ceed_cuda + ceed_hip <= 1,
                "Only one CEED backend can be enabled at a time!");
 #endif
    if (Allows(Backend::CUDA)) { CudaDeviceSetup(dev, ngpu); }
    if (Allows(Backend::HIP)) { HipDeviceSetup(dev, ngpu); }
-   if (Allows(Backend::RAJA_CUDA)) { RajaDeviceSetup(dev, ngpu); }
+   if (Allows(Backend::RAJA_CUDA) || Allows(Backend::RAJA_HIP))
+   { RajaDeviceSetup(dev, ngpu); }
    // The check for MFEM_USE_OCCA is in the function OccaDeviceSetup().
    if (Allows(Backend::OCCA_MASK)) { OccaDeviceSetup(dev); }
    if (Allows(Backend::CEED_CPU))
@@ -504,6 +564,17 @@ void Device::Setup(const int device)
          CeedDeviceSetup(device_option);
       }
    }
+   if (Allows(Backend::CEED_HIP))
+   {
+      if (!device_option)
+      {
+         CeedDeviceSetup("/gpu/hip");
+      }
+      else
+      {
+         CeedDeviceSetup(device_option);
+      }
+   }
    if (Allows(Backend::DEBUG)) { ngpu = 1; }
 }
 

general/device.hpp

@@ -1,4 +1,4 @@
-// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 // LICENSE and NOTICE for details. LLNL-CODE-806117.
 //
@@ -46,26 +46,33 @@ struct Backend
       /** @brief [device] RAJA CUDA backend. Enabled when MFEM_USE_RAJA = YES
           and MFEM_USE_CUDA = YES. */
       RAJA_CUDA = 1 << 6,
+      /** @brief [device] RAJA HIP backend. Enabled when MFEM_USE_RAJA = YES
+          and MFEM_USE_HIP = YES. */
+      RAJA_HIP = 1 << 7,
       /** @brief [host] OCCA CPU backend: sequential execution on each MPI rank.
           Enabled when MFEM_USE_OCCA = YES. */
-      OCCA_CPU = 1 << 7,
+      OCCA_CPU = 1 << 8,
       /// [host] OCCA OpenMP backend. Enabled when MFEM_USE_OCCA = YES.
-      OCCA_OMP = 1 << 8,
+      OCCA_OMP = 1 << 9,
       /** @brief [device] OCCA CUDA backend. Enabled when MFEM_USE_OCCA = YES
           and MFEM_USE_CUDA = YES. */
-      OCCA_CUDA = 1 << 9,
+      OCCA_CUDA = 1 << 10,
       /** @brief [host] CEED CPU backend. GPU backends can still be used, but
           with expensive memory transfers. Enabled when MFEM_USE_CEED = YES. */
-      CEED_CPU  = 1 << 10,
+      CEED_CPU  = 1 << 11,
       /** @brief [device] CEED CUDA backend working together with the CUDA
           backend. Enabled when MFEM_USE_CEED = YES and MFEM_USE_CUDA = YES.
-          NOTE: The current default libCEED GPU backend is non-deterministic! */
-      CEED_CUDA = 1 << 11,
+          NOTE: The current default libCEED CUDA backend is non-deterministic! */
+      CEED_CUDA = 1 << 12,
+      /** @brief [device] CEED HIP backend working together with the HIP
+          backend. Enabled when MFEM_USE_CEED = YES and MFEM_USE_HIP = YES. */
+      CEED_HIP = 1 << 13,
       /** @brief [device] Debug backend: host memory is READ/WRITE protected
           while a device is in use. It allows to test the "device" code-path
           (using separate host/device memory pools and host <-> device
-          transfers) without any GPU hardware. */
-      DEBUG = 1 << 12
+          transfers) without any GPU hardware. As 'DEBUG' is sometimes used
+          as a macro, `_DEVICE` has been added to avoid conflicts. */
+      DEBUG = 1 << 14
    };
 
    /** @brief Additional useful constants. For example, the *_MASK constants can
@@ -73,23 +80,23 @@ struct Backend
    enum
    {
       /// Number of backends: from (1 << 0) to (1 << (NUM_BACKENDS-1)).
-      NUM_BACKENDS = 13,
+      NUM_BACKENDS = 15,
 
       /// Biwise-OR of all CPU backends
       CPU_MASK = CPU | RAJA_CPU | OCCA_CPU | CEED_CPU,
       /// Biwise-OR of all CUDA backends
       CUDA_MASK = CUDA | RAJA_CUDA | OCCA_CUDA | CEED_CUDA,
       /// Biwise-OR of all HIP backends
-      HIP_MASK = HIP,
+      HIP_MASK = HIP | RAJA_HIP | CEED_HIP,
       /// Biwise-OR of all OpenMP backends
       OMP_MASK = OMP | RAJA_OMP | OCCA_OMP,
       /// Bitwise-OR of all CEED backends
-      CEED_MASK = CEED_CPU | CEED_CUDA,
+      CEED_MASK = CEED_CPU | CEED_CUDA | CEED_HIP,
       /// Biwise-OR of all device backends
       DEVICE_MASK = CUDA_MASK | HIP_MASK | DEBUG,
 
       /// Biwise-OR of all RAJA backends
-      RAJA_MASK = RAJA_CPU | RAJA_OMP | RAJA_CUDA,
+      RAJA_MASK = RAJA_CPU | RAJA_OMP | RAJA_CUDA | RAJA_HIP,
       /// Biwise-OR of all OCCA backends
       OCCA_MASK = OCCA_CPU | OCCA_OMP | OCCA_CUDA
    };
@@ -117,22 +124,25 @@ private:
    friend class MemoryManager;
    enum MODES {SEQUENTIAL, ACCELERATED};
 
-   static bool device_env, mem_host_env, mem_device_env;
+   static bool device_env, mem_host_env, mem_device_env, mem_types_set;
    static Device device_singleton;
 
-   MODES mode;
-   int dev = 0; ///< Device ID of the configured device.
+   MODES mode = Device::SEQUENTIAL;
+   int dev = 0;   ///< Device ID of the configured device.
    int ngpu = -1; ///< Number of detected devices; -1: not initialized.
-   unsigned long backends; ///< Bitwise-OR of all configured backends.
+   /// Bitwise-OR of all configured backends.
+   unsigned long backends = Backend::CPU;
    /// Set to true during configuration, except in 'device_singleton'.
-   bool destroy_mm;
-   bool mpi_gpu_aware;
+   bool destroy_mm = false;
+   bool mpi_gpu_aware = false;
 
-   MemoryType host_mem_type;      ///< Current Host MemoryType
-   MemoryClass host_mem_class;    ///< Current Host MemoryClass
+   MemoryType host_mem_type = MemoryType::HOST;    ///< Current Host MemoryType
+   MemoryClass host_mem_class = MemoryClass::HOST; ///< Current Host MemoryClass
 
-   MemoryType device_mem_type;    ///< Current Device MemoryType
-   MemoryClass device_mem_class;  ///< Current Device MemoryClass
+   /// Current Device MemoryType
+   MemoryType device_mem_type = MemoryType::HOST;
+   /// Current Device MemoryClass
+   MemoryClass device_mem_class = MemoryClass::HOST;
 
    char *device_option = NULL;
    Device(Device const&);
@@ -173,14 +183,6 @@ public:
        @note This object should be destroyed after all other MFEM objects that
        use the Device are destroyed. */
    Device(const std::string &device, const int dev = 0)
-      : mode(Device::SEQUENTIAL),
-        backends(Backend::CPU),
-        destroy_mm(false),
-        mpi_gpu_aware(false),
-        host_mem_type(MemoryType::HOST),
-        host_mem_class(MemoryClass::HOST),
-        device_mem_type(MemoryType::HOST),
-        device_mem_class(MemoryClass::HOST)
    { Configure(device, dev); }
 
    /// Destructor.
@@ -193,25 +195,44 @@ public:
        * The available backends are described by the Backend class.
        * The string name of a backend is the lowercase version of the
          Backend::Id enumeration constant with '_' replaced by '-', e.g. the
-         string name of 'RAJA_CPU' is 'raja-cpu'.
+         string name of 'RAJA_CPU' is 'raja-cpu'. The string name of the debug
+         backend (Backend::Id 'DEBUG_DEVICE') is exceptionally set to 'debug'.
        * The 'cpu' backend is always enabled with lowest priority.
        * The current backend priority from highest to lowest is:
-         'ceed-cuda', 'occa-cuda', 'raja-cuda', 'cuda', 'hip', 'debug',
+         'ceed-cuda', 'occa-cuda', 'raja-cuda', 'cuda',
+         'ceed-hip', 'hip', 'debug',
          'occa-omp', 'raja-omp', 'omp',
          'ceed-cpu', 'occa-cpu', 'raja-cpu', 'cpu'.
        * Multiple backends can be configured at the same time.
        * Only one 'occa-*' backend can be configured at a time.
        * The backend 'occa-cuda' enables the 'cuda' backend unless 'raja-cuda'
          is already enabled.
+       * The backend 'occa-omp' enables the 'omp' backend (if MFEM was built
+         with MFEM_USE_OPENMP=YES) unless 'raja-omp' is already enabled.
+       * Only one 'ceed-*' backend can be configured at a time.
        * The backend 'ceed-cpu' delegates to a libCEED CPU backend the setup and
          evaluation of the operator.
        * The backend 'ceed-cuda' delegates to a libCEED CUDA backend the setup
-         and evaluation of the operator and enables the 'cuda' backend to avoid
-         transfer between host and device.
+         and evaluation of operators and enables the 'cuda' backend to avoid
+         transfers between host and device.
+       * The backend 'ceed-hip' delegates to a libCEED HIP backend the setup
+         and evaluation of operators and enables the 'hip' backend to avoid
+         transfers between host and device.
        * The 'debug' backend should not be combined with other device backends.
    */
    void Configure(const std::string &device, const int dev = 0);
 
+   /// Set the default host and device MemoryTypes, @a h_mt and @a d_mt.
+   /** The host and device MemoryTypes are also set to be dual to each other.
+
+       These two MemoryType%s are used by most MFEM classes when allocating
+       memory used on host and device, respectively.
+
+       This method can only be called before Device construction and
+       configuration, and the specified memory types must be compatible with
+       the subsequent Device configuration. */
+   static void SetMemoryTypes(MemoryType h_mt, MemoryType d_mt);
+
    /// Print the configuration of the MFEM virtual device object.
    void Print(std::ostream &out = mfem::out);
 

general/forall.hpp

@@ -23,8 +23,8 @@ namespace mfem
 {
 
 // Maximum size of dofs and quads in 1D.
-const int MAX_D1D = 14;
-const int MAX_Q1D = 14;
+const int MAX_D1D = 8;
+const int MAX_Q1D = 8;
 
 // MFEM pragma macros that can be used inside MFEM_FORALL macros.
 #define MFEM_PRAGMA(X) _Pragma(#X)
@@ -89,7 +89,11 @@ void OmpWrap(const int N, HBODY &&h_body)
 /// RAJA Cuda backend
 #if defined(MFEM_USE_RAJA) && defined(RAJA_ENABLE_CUDA)
 
+#if RAJA_VERSION_MAJOR == 0 && RAJA_VERSION_MINOR < 12
 using RAJA::statement::Segs;
+#else
+using RAJA::Segs;
+#endif
 
 template <const int BLOCKS = MFEM_CUDA_BLOCKS, typename DBODY>
 void RajaCudaWrap1D(const int N, DBODY &&d_body)
@@ -147,7 +151,12 @@ void RajaCudaWrap3D(const int N, DBODY &&d_body,
 /// RAJA OpenMP backend
 #if defined(MFEM_USE_RAJA) && defined(RAJA_ENABLE_OPENMP)
 
+#if RAJA_VERSION_MAJOR == 0 && RAJA_VERSION_MINOR < 12
 using RAJA::statement::Segs;
+#else
+using RAJA::Segs;
+#endif
+
 
 template <typename HBODY>
 void RajaOmpWrap(const int N, HBODY &&h_body)

general/hip.cpp

@@ -1,4 +1,4 @@
-// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 // LICENSE and NOTICE for details. LLNL-CODE-806117.
 //
@@ -61,6 +61,21 @@ void* HipMallocManaged(void** dptr, size_t bytes)
    return *dptr;
 }
 
+void* HipMemAllocHostPinned(void** ptr, size_t bytes)
+{
+#ifdef MFEM_USE_HIP
+#ifdef MFEM_TRACK_HIP_MEM
+   mfem::out << "HipMemAllocHostPinned(): allocating " << bytes << " bytes ... "
+             << std::flush;
+#endif
+   MFEM_GPU_CHECK(hipHostMalloc(ptr, bytes, hipHostMallocDefault));
+#ifdef MFEM_TRACK_HIP_MEM
+   mfem::out << "done: " << *ptr << std::endl;
+#endif
+#endif
+   return *ptr;
+}
+
 void* HipMemFree(void *dptr)
 {
 #ifdef MFEM_USE_HIP
@@ -76,6 +91,21 @@ void* HipMemFree(void *dptr)
    return dptr;
 }
 
+void* HipMemFreeHostPinned(void *ptr)
+{
+#ifdef MFEM_USE_HIP
+#ifdef MFEM_TRACK_HIP_MEM
+   mfem::out << "HipMemFreeHostPinned(): deallocating memory @ " << ptr << " ... "
+             << std::flush;
+#endif
+   MFEM_GPU_CHECK(hipHostFree(ptr));
+#ifdef MFEM_TRACK_HIP_MEM
+   mfem::out << "done." << std::endl;
+#endif
+#endif
+   return ptr;
+}
+
 void* HipMemcpyHtoD(void* dst, const void* src, size_t bytes)
 {
 #ifdef MFEM_USE_HIP
@@ -125,7 +155,7 @@ void* HipMemcpyDtoDAsync(void* dst, const void *src, size_t bytes)
 void* HipMemcpyDtoH(void *dst, const void *src, size_t bytes)
 {
 #ifdef MFEM_USE_HIP
-#ifdef MFEM_TRACK_HPI_MEM
+#ifdef MFEM_TRACK_HIP_MEM
    mfem::out << "HipMemcpyDtoH(): copying " << bytes << " bytes from "
              << src << " to " << dst << " ... " << std::flush;
 #endif

general/hip.hpp

@@ -1,4 +1,4 @@
-// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 // LICENSE and NOTICE for details. LLNL-CODE-806117.
 //
@@ -43,7 +43,6 @@
 #if defined(MFEM_USE_HIP) && defined(__HIP_DEVICE_COMPILE__)
 #define MFEM_SHARED __shared__
 #define MFEM_SYNC_THREAD __syncthreads()
-#define MFEM_BLOCK_ID(k) hipBlockIdx_ ##k
 #define MFEM_THREAD_ID(k) hipThreadIdx_ ##k
 #define MFEM_THREAD_SIZE(k) hipBlockDim_ ##k
 #define MFEM_FOREACH_THREAD(i,k,N) \
@@ -65,9 +64,15 @@ void* HipMemAlloc(void **d_ptr, size_t bytes);
 /// Allocates managed device memory
 void* HipMallocManaged(void **d_ptr, size_t bytes);
 
+/// Allocates page-locked (pinned) host memory
+void* HipMemAllocHostPinned(void **ptr, size_t bytes);
+
 /// Frees device memory
 void* HipMemFree(void *d_ptr);
 
+/// Frees page-locked (pinned) host memory and returns destination ptr.
+void* HipMemFreeHostPinned(void *ptr);
+
 /// Copies memory from Host to Device
 void* HipMemcpyHtoD(void *d_dst, const void *h_src, size_t bytes);
 

general/mem_manager.cpp

@@ -1,4 +1,4 @@
-// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 // LICENSE and NOTICE for details. LLNL-CODE-806117.
 //
@@ -62,38 +62,39 @@ MemoryType GetMemoryType(MemoryClass mc)
    return MemoryType::HOST;
 }
 
-// We want to keep this pairs, as it is checked in MFEM_VERIFY_TYPES
-MemoryType MemoryManager::GetDualMemoryType_(MemoryType mt)
-{
-   switch (mt)
-   {
-      case MemoryType::HOST:           return MemoryType::DEVICE;
-      case MemoryType::HOST_32:        return MemoryType::DEVICE;
-      case MemoryType::HOST_64:        return MemoryType::DEVICE;
-      case MemoryType::HOST_DEBUG:     return MemoryType::DEVICE_DEBUG;
-      case MemoryType::HOST_UMPIRE:    return MemoryType::DEVICE_UMPIRE;
-      case MemoryType::MANAGED:        return MemoryType::MANAGED;
-      case MemoryType::DEVICE:         return MemoryType::HOST;
-      case MemoryType::DEVICE_DEBUG:   return MemoryType::HOST_DEBUG;
-      case MemoryType::DEVICE_UMPIRE:  return MemoryType::HOST_UMPIRE;
-      default: mfem_error("Unknown memory type!");
-   }
-   MFEM_VERIFY(false,"");
-   return MemoryType::HOST;
-}
 
 static void MFEM_VERIFY_TYPES(const MemoryType h_mt, const MemoryType d_mt)
 {
-   MFEM_ASSERT(IsHostMemory(h_mt),"");
-   MFEM_ASSERT(IsDeviceMemory(d_mt),"");
+   MFEM_VERIFY(IsHostMemory(h_mt), "h_mt = " << (int)h_mt);
+   MFEM_VERIFY(IsDeviceMemory(d_mt) || d_mt == MemoryType::DEFAULT,
+               "d_mt = " << (int)d_mt);
+   // If h_mt == MemoryType::HOST_DEBUG, then d_mt == MemoryType::DEVICE_DEBUG
+   //                                      or d_mt == MemoryType::DEFAULT
+   MFEM_VERIFY(h_mt != MemoryType::HOST_DEBUG ||
+               d_mt == MemoryType::DEVICE_DEBUG ||
+               d_mt == MemoryType::DEFAULT,
+               "d_mt = " << MemoryTypeName[(int)d_mt]);
+   // If d_mt == MemoryType::DEVICE_DEBUG, then h_mt == MemoryType::HOST_DEBUG
+   MFEM_VERIFY(d_mt != MemoryType::DEVICE_DEBUG ||
+               h_mt == MemoryType::HOST_DEBUG,
+               "h_mt = " << MemoryTypeName[(int)h_mt]);
+#if 0
    const bool sync =
+      (h_mt == MemoryType::HOST_PINNED && d_mt == MemoryType::DEVICE) ||
+      (h_mt == MemoryType::HOST_PINNED && d_mt == MemoryType::DEVICE_UMPIRE) ||
+      (h_mt == MemoryType::HOST_PINNED && d_mt == MemoryType::DEVICE_UMPIRE_2) ||
+      (h_mt == MemoryType::HOST_UMPIRE && d_mt == MemoryType::DEVICE) ||
       (h_mt == MemoryType::HOST_UMPIRE && d_mt == MemoryType::DEVICE_UMPIRE) ||
+      (h_mt == MemoryType::HOST_UMPIRE && d_mt == MemoryType::DEVICE_UMPIRE_2) ||
       (h_mt == MemoryType::HOST_DEBUG && d_mt == MemoryType::DEVICE_DEBUG) ||
       (h_mt == MemoryType::MANAGED && d_mt == MemoryType::MANAGED) ||
       (h_mt == MemoryType::HOST_64 && d_mt == MemoryType::DEVICE) ||
       (h_mt == MemoryType::HOST_32 && d_mt == MemoryType::DEVICE) ||
-      (h_mt == MemoryType::HOST && d_mt == MemoryType::DEVICE);
+      (h_mt == MemoryType::HOST && d_mt == MemoryType::DEVICE) ||
+      (h_mt == MemoryType::HOST && d_mt == MemoryType::DEVICE_UMPIRE) ||
+      (h_mt == MemoryType::HOST && d_mt == MemoryType::DEVICE_UMPIRE_2);
    MFEM_VERIFY(sync, "");
+#endif
 }
 
 MemoryClass operator*(MemoryClass mc1, MemoryClass mc2)
@@ -135,9 +136,12 @@ struct Memory
    void *const h_ptr;
    void *d_ptr;
    const size_t bytes;
-   const MemoryType h_mt, d_mt;
+   const MemoryType h_mt;
+   MemoryType d_mt;
+   mutable bool h_rw, d_rw;
    Memory(void *p, size_t b, MemoryType h, MemoryType d):
-      h_ptr(p), d_ptr(nullptr), bytes(b), h_mt(h), d_mt(d) { }
+      h_ptr(p), d_ptr(nullptr), bytes(b), h_mt(h), d_mt(d),
+      h_rw(true), d_rw(true) { }
 };
 
 /// Alias class that holds the base memory region and the offset
@@ -173,8 +177,8 @@ public:
    virtual ~HostMemorySpace() { }
    virtual void Alloc(void **ptr, size_t bytes) { *ptr = std::malloc(bytes); }
    virtual void Dealloc(void *ptr) { std::free(ptr); }
-   virtual void Protect(const void*, size_t) { }
-   virtual void Unprotect(const void*, size_t) { }
+   virtual void Protect(const Memory&, size_t) { }
+   virtual void Unprotect(const Memory&, size_t) { }
    virtual void AliasProtect(const void*, size_t) { }
    virtual void AliasUnprotect(const void*, size_t) { }
 };
@@ -322,18 +326,16 @@ inline void MmuDealloc(void *ptr, const size_t bytes)
 /// MMU protection, through ::mprotect with no read/write accesses
 inline void MmuProtect(const void *ptr, const size_t bytes)
 {
-   static const bool mmu_protect_error = getenv("MFEM_MMU_PROTECT_ERROR");
    if (!::mprotect(const_cast<void*>(ptr), bytes, PROT_NONE)) { return; }
-   if (mmu_protect_error) { mfem_error("MMU protection (NONE) error"); }
+   mfem_error("MMU protection (NONE) error");
 }
 
 /// MMU un-protection, through ::mprotect with read/write accesses
 inline void MmuAllow(const void *ptr, const size_t bytes)
 {
    const int RW = PROT_READ | PROT_WRITE;
-   static const bool mmu_protect_error = getenv("MFEM_MMU_PROTECT_ERROR");
    if (!::mprotect(const_cast<void*>(ptr), bytes, RW)) { return; }
-   if (mmu_protect_error) { mfem_error("MMU protection (R/W) error"); }
+   mfem_error("MMU protection (R/W) error");
 }
 #else
 inline void MmuInit() { }
@@ -354,8 +356,10 @@ public:
    MmuHostMemorySpace(): HostMemorySpace() { MmuInit(); }
    void Alloc(void **ptr, size_t bytes) { MmuAlloc(ptr, bytes); }
    void Dealloc(void *ptr) { MmuDealloc(ptr, maps->memories.at(ptr).bytes); }
-   void Protect(const void *ptr, size_t bytes) { MmuProtect(ptr, bytes); }
-   void Unprotect(const void *ptr, size_t bytes) { MmuAllow(ptr, bytes); }
+   void Protect(const Memory& mem, size_t bytes)
+   { if (mem.h_rw) { mem.h_rw = false; MmuProtect(mem.h_ptr, bytes); } }
+   void Unprotect(const Memory &mem, size_t bytes)
+   { if (!mem.h_rw) { mem.h_rw = true; MmuAllow(mem.h_ptr, bytes); } }
    /// Aliases need to be restricted during protection
    void AliasProtect(const void *ptr, size_t bytes)
    { MmuProtect(MmuAddrR(ptr), MmuLengthR(ptr, bytes)); }
@@ -402,6 +406,31 @@ public:
    { return CuMemcpyDtoH(dst, src, bytes); }
 };
 
+/// The CUDA/HIP page-locked host memory space
+class HostPinnedMemorySpace: public HostMemorySpace
+{
+public:
+   HostPinnedMemorySpace(): HostMemorySpace() { }
+   void Alloc(void ** ptr, size_t bytes) override
+   {
+#ifdef MFEM_USE_CUDA
+      CuMemAllocHostPinned(ptr, bytes);
+#endif
+#ifdef MFEM_USE_HIP
+      HipMemAllocHostPinned(ptr, bytes);
+#endif
+   }
+   void Dealloc(void *ptr) override
+   {
+#ifdef MFEM_USE_CUDA
+      CuMemFreeHostPinned(ptr);
+#endif
+#ifdef MFEM_USE_HIP
+      HipMemFreeHostPinned(ptr);
+#endif
+   }
+};
+
 /// The HIP device memory space
 class HipDeviceMemorySpace: public DeviceMemorySpace
 {
@@ -444,8 +473,10 @@ public:
    MmuDeviceMemorySpace(): DeviceMemorySpace() { }
    void Alloc(Memory &m) { MmuAlloc(&m.d_ptr, m.bytes); }
    void Dealloc(Memory &m) { MmuDealloc(m.d_ptr, m.bytes); }
-   void Protect(const Memory &m) { MmuProtect(m.d_ptr, m.bytes); }
-   void Unprotect(const Memory &m) { MmuAllow(m.d_ptr, m.bytes); }
+   void Protect(const Memory &m)
+   { if (m.d_rw) { m.d_rw = false; MmuProtect(m.d_ptr, m.bytes); } }
+   void Unprotect(const Memory &m)
+   { if (!m.d_rw) { m.d_rw = true; MmuAllow(m.d_ptr, m.bytes); } }
    /// Aliases need to be restricted during protection
    void AliasProtect(const void *ptr, size_t bytes)
    { MmuProtect(MmuAddrR(ptr), MmuLengthR(ptr, bytes)); }
@@ -460,54 +491,64 @@ public:
    { return std::memcpy(dst, src, bytes); }
 };
 
-#ifndef MFEM_USE_UMPIRE
-class UmpireHostMemorySpace : public NoHostMemorySpace { };
-class UmpireDeviceMemorySpace : public NoDeviceMemorySpace { };
-#else
+#ifdef MFEM_USE_UMPIRE
+class UmpireMemorySpace
+{
+protected:
+   umpire::ResourceManager &rm;
+   umpire::Allocator allocator;
+   bool owns_allocator{false};
+
+public:
+   // TODO: this only releases unused memory
+   virtual ~UmpireMemorySpace() { if (owns_allocator) { allocator.release(); } }
+   UmpireMemorySpace(const char * name, const char * space)
+      : rm(umpire::ResourceManager::getInstance())
+   {
+      if (!rm.isAllocator(name))
+      {
+         allocator = rm.makeAllocator<umpire::strategy::DynamicPool>(
+                        name, rm.getAllocator(space));
+         owns_allocator = true;
+      }
+      else
+      {
+         allocator = rm.getAllocator(name);
+         owns_allocator = false;
+      }
+   }
+};
+
 /// The Umpire host memory space
-class UmpireHostMemorySpace : public HostMemorySpace
+class UmpireHostMemorySpace : public HostMemorySpace, public UmpireMemorySpace
 {
 private:
-   const char *name;
-   umpire::ResourceManager &rm;
-   umpire::Allocator h_allocator;
    umpire::strategy::AllocationStrategy *strat;
 public:
-   ~UmpireHostMemorySpace() { h_allocator.release(); }
-   UmpireHostMemorySpace():
-      HostMemorySpace(),
-      name(mm.GetUmpireAllocatorHostName()),
-      rm(umpire::ResourceManager::getInstance()),
-      h_allocator(rm.isAllocator(name)? rm.getAllocator(name):
-                  rm.makeAllocator<umpire::strategy::DynamicPool>
-                  (name, rm.getAllocator("HOST"))),
-      strat(h_allocator.getAllocationStrategy()) { }
-   void Alloc(void **ptr, size_t bytes) { *ptr = h_allocator.allocate(bytes); }
-   void Dealloc(void *ptr) { h_allocator.deallocate(ptr); }
+   UmpireHostMemorySpace(const char * name)
+      : HostMemorySpace(),
+        UmpireMemorySpace(name, "HOST"),
+        strat(allocator.getAllocationStrategy()) {}
+   void Alloc(void **ptr, size_t bytes) override
+   { *ptr = allocator.allocate(bytes); }
+   void Dealloc(void *ptr) override { allocator.deallocate(ptr); }
    void Insert(void *ptr, size_t bytes)
    { rm.registerAllocation(ptr, {ptr, bytes, strat}); }
 };
 
 /// The Umpire device memory space
-#ifdef MFEM_USE_CUDA
-class UmpireDeviceMemorySpace : public DeviceMemorySpace
+#if defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP)
+class UmpireDeviceMemorySpace : public DeviceMemorySpace,
+   public UmpireMemorySpace
 {
-private:
-   const char *name;
-   umpire::ResourceManager &rm;
-   umpire::Allocator d_allocator;
 public:
-   ~UmpireDeviceMemorySpace() { d_allocator.release(); }
-   UmpireDeviceMemorySpace():
-      DeviceMemorySpace(),
-      name(mm.GetUmpireAllocatorDeviceName()),
-      rm(umpire::ResourceManager::getInstance()),
-      d_allocator(rm.isAllocator(name)? rm.getAllocator(name):
-                  rm.makeAllocator<umpire::strategy::DynamicPool>
-                  (name, rm.getAllocator("DEVICE"))) { }
-   void Alloc(Memory &base) { base.d_ptr = d_allocator.allocate(base.bytes); }
-   void Dealloc(Memory &base) { d_allocator.deallocate(base.d_ptr); }
-   void *HtoD(void *dst, const void *src, size_t bytes)
+   UmpireDeviceMemorySpace(const char * name)
+      : DeviceMemorySpace(),
+        UmpireMemorySpace(name, "DEVICE") {}
+   void Alloc(Memory &base) override
+   { base.d_ptr = allocator.allocate(base.bytes); }
+   void Dealloc(Memory &base) override { rm.deallocate(base.d_ptr); }
+   void *HtoD(void *dst, const void *src, size_t bytes) override
    {
 #ifdef MFEM_USE_CUDA
       return CuMemcpyHtoD(dst, src, bytes);
@@ -517,7 +558,7 @@ public:
 #endif
       //rm.copy(dst, const_cast<void*>(src), bytes); return dst;
    }
-   void *DtoD(void* dst, const void* src, size_t bytes)
+   void *DtoD(void* dst, const void* src, size_t bytes) override
    {
 #ifdef MFEM_USE_CUDA
       return CuMemcpyDtoD(dst, src, bytes);
@@ -527,7 +568,7 @@ public:
 #endif
       //rm.copy(dst, const_cast<void*>(src), bytes); return dst;
    }
-   void *DtoH(void *dst, const void *src, size_t bytes)
+   void *DtoH(void *dst, const void *src, size_t bytes) override
    {
 #ifdef MFEM_USE_CUDA
       return CuMemcpyDtoH(dst, src, bytes);
@@ -539,8 +580,12 @@ public:
    }
 };
 #else
-class UmpireDeviceMemorySpace : public NoDeviceMemorySpace { };
-#endif // MFEM_USE_CUDA
+class UmpireDeviceMemorySpace : public NoDeviceMemorySpace
+{
+public:
+   UmpireDeviceMemorySpace(const char * /*unused*/) {}
+};
+#endif // MFEM_USE_CUDA || MFEM_USE_HIP
 #endif // MFEM_USE_UMPIRE
 
 /// Memory space controller class
@@ -570,7 +615,7 @@ public:
       host[static_cast<int>(MT::HOST_64)] = new Aligned64HostMemorySpace();
       // HOST_DEBUG is delayed, as it reroutes signals
       host[static_cast<int>(MT::HOST_DEBUG)] = nullptr;
-      host[static_cast<int>(MT::HOST_UMPIRE)] = new UmpireHostMemorySpace();
+      host[static_cast<int>(MT::HOST_UMPIRE)] = nullptr;
       host[static_cast<int>(MT::MANAGED)] = new UvmHostMemorySpace();
 
       // Filling the device memory backends, shifting with the device size
@@ -580,6 +625,7 @@ public:
       device[static_cast<int>(MemoryType::DEVICE)-shift] = nullptr;
       device[static_cast<int>(MT::DEVICE_DEBUG)-shift] = nullptr;
       device[static_cast<int>(MT::DEVICE_UMPIRE)-shift] = nullptr;
+      device[static_cast<int>(MT::DEVICE_UMPIRE_2)-shift] = nullptr;
    }
 
    HostMemorySpace* Host(const MemoryType mt)
@@ -612,8 +658,19 @@ public:
 private:
    HostMemorySpace* NewHostCtrl(const MemoryType mt)
    {
-      if (mt == MT::HOST_DEBUG) { return new MmuHostMemorySpace(); }
-      MFEM_ABORT("Unknown host memory controller!");
+      switch (mt)
+      {
+         case MT::HOST_DEBUG: return new MmuHostMemorySpace();
+#ifdef MFEM_USE_UMPIRE
+         case MT::HOST_UMPIRE:
+            return new UmpireHostMemorySpace(
+                      MemoryManager::GetUmpireHostAllocatorName());
+#else
+         case MT::HOST_UMPIRE: return new NoHostMemorySpace();
+#endif
+         case MT::HOST_PINNED: return new HostPinnedMemorySpace();
+         default: MFEM_ABORT("Unknown host memory controller!");
+      }
       return nullptr;
    }
 
@@ -621,7 +678,17 @@ private:
    {
       switch (mt)
       {
-         case MT::DEVICE_UMPIRE: return new UmpireDeviceMemorySpace();
+#ifdef MFEM_USE_UMPIRE
+         case MT::DEVICE_UMPIRE:
+            return new UmpireDeviceMemorySpace(
+                      MemoryManager::GetUmpireDeviceAllocatorName());
+         case MT::DEVICE_UMPIRE_2:
+            return new UmpireDeviceMemorySpace(
+                      MemoryManager::GetUmpireDevice2AllocatorName());
+#else
+         case MT::DEVICE_UMPIRE: return new NoDeviceMemorySpace();
+         case MT::DEVICE_UMPIRE_2: return new NoDeviceMemorySpace();
+#endif
          case MT::DEVICE_DEBUG: return new MmuDeviceMemorySpace();
          case MT::DEVICE:
          {
@@ -648,20 +715,50 @@ void *MemoryManager::New_(void *h_tmp, size_t bytes, MemoryType mt,
                           unsigned &flags)
 {
    MFEM_ASSERT(exists, "Internal error!");
-   MFEM_ASSERT(mt != MemoryType::HOST, "Internal error!");
-   const bool is_host_mem = IsHostMemory(mt);
-   const MemType dual_mt = GetDualMemoryType_(mt);
-   const MemType h_mt = is_host_mem ? mt : dual_mt;
-   const MemType d_mt = is_host_mem ? dual_mt : mt;
-   MFEM_VERIFY_TYPES(h_mt, d_mt);
-   void *h_ptr = h_tmp;
+   if (IsHostMemory(mt))
+   {
+      MFEM_ASSERT(mt != MemoryType::HOST && h_tmp == nullptr,
+                  "Internal error!");
+      // d_mt = MemoryType::DEFAULT means d_mt = GetDualMemoryType(h_mt),
+      // evaluated at the time when the device pointer is allocated, see
+      // GetDevicePtr() and GetAliasDevicePtr()
+      const MemoryType d_mt = MemoryType::DEFAULT;
+      // We rely on the next call using lazy dev alloc
+      return New_(h_tmp, bytes, mt, d_mt, Mem::VALID_HOST, flags);
+   }
+   else
+   {
+      const MemoryType h_mt = GetDualMemoryType(mt);
+      return New_(h_tmp, bytes, h_mt, mt, Mem::VALID_DEVICE, flags);
+   }
+}
+
+void *MemoryManager::New_(void *h_tmp, size_t bytes, MemoryType h_mt,
+                          MemoryType d_mt, unsigned valid_flags,
+                          unsigned &flags)
+{
+   MFEM_ASSERT(exists, "Internal error!");
+   MFEM_ASSERT(IsHostMemory(h_mt), "h_mt must be host type");
+   MFEM_ASSERT(IsDeviceMemory(d_mt) || d_mt == h_mt ||
+               d_mt == MemoryType::DEFAULT,
+               "d_mt must be device type, the same is h_mt, or DEFAULT");
+   MFEM_ASSERT((h_mt != MemoryType::HOST || h_tmp != nullptr) &&
+               (h_mt == MemoryType::HOST || h_tmp == nullptr),
+               "Internal error");
+   MFEM_ASSERT((valid_flags & ~(Mem::VALID_HOST | Mem::VALID_DEVICE)) == 0,
+               "Internal error");
+   void *h_ptr;
    if (h_tmp == nullptr) { ctrl->Host(h_mt)->Alloc(&h_ptr, bytes); }
-   flags = Mem::REGISTERED;
-   flags |= Mem::OWNS_INTERNAL | Mem::OWNS_HOST | Mem::OWNS_DEVICE;
-   flags |= is_host_mem ? Mem::VALID_HOST : Mem::VALID_DEVICE;
-   if (is_host_mem) { mm.Insert(h_ptr, bytes, h_mt, d_mt); }
-   else { mm.InsertDevice(nullptr, h_ptr, bytes, h_mt, d_mt); }
+   else { h_ptr = h_tmp; }
+   flags = Mem::REGISTERED | Mem::OWNS_INTERNAL | Mem::OWNS_HOST |
+           Mem::OWNS_DEVICE | valid_flags;
+   // The other New_() method relies on this lazy allocation behavior.
+   mm.Insert(h_ptr, bytes, h_mt, d_mt); // lazy dev alloc
+   // mm.InsertDevice(nullptr, h_ptr, bytes, h_mt, d_mt); // non-lazy dev alloc
+
+   // MFEM_VERIFY_TYPES(h_mt, mt); // done by mm.Insert() above
    CheckHostMemoryType_(h_mt, h_ptr);
+
    return h_ptr;
 }
 
@@ -673,9 +770,12 @@ void *MemoryManager::Register_(void *ptr, void *h_tmp, size_t bytes,
    MFEM_ASSERT(exists, "Internal error!");
    MFEM_ASSERT(!alias, "Cannot register an alias!");
    const bool is_host_mem = IsHostMemory(mt);
-   const MemType dual_mt = GetDualMemoryType_(mt);
-   const MemType h_mt = is_host_mem ? mt : dual_mt;
-   const MemType d_mt = is_host_mem ? dual_mt : mt;
+   const MemType h_mt = is_host_mem ? mt : GetDualMemoryType(mt);
+   const MemType d_mt = is_host_mem ? MemoryType::DEFAULT : mt;
+   // d_mt = MemoryType::DEFAULT means d_mt = GetDualMemoryType(h_mt),
+   // evaluated at the time when the device pointer is allocated, see
+   // GetDevicePtr() and GetAliasDevicePtr()
+
    MFEM_VERIFY_TYPES(h_mt, d_mt);
 
    if (ptr == nullptr && h_tmp == nullptr)
@@ -696,17 +796,42 @@ void *MemoryManager::Register_(void *ptr, void *h_tmp, size_t bytes,
    }
    else // DEVICE TYPES
    {
-      h_ptr = h_tmp;
-      if (own && h_tmp == nullptr) { ctrl->Host(h_mt)->Alloc(&h_ptr, bytes); }
+      MFEM_VERIFY(ptr, "cannot register NULL device pointer");
+      if (h_tmp == nullptr) { ctrl->Host(h_mt)->Alloc(&h_ptr, bytes); }
+      else { h_ptr = h_tmp; }
       mm.InsertDevice(ptr, h_ptr, bytes, h_mt, d_mt);
       flags = own ? flags | Mem::OWNS_DEVICE : flags & ~Mem::OWNS_DEVICE;
-      flags = own ? flags | Mem::OWNS_HOST   : flags & ~Mem::OWNS_HOST;
-      flags |= Mem::VALID_DEVICE;
+      flags |= (Mem::OWNS_HOST | Mem::VALID_DEVICE);
    }
    CheckHostMemoryType_(h_mt, h_ptr);
    return h_ptr;
 }
 
+void MemoryManager::Register_(void *h_ptr, void *d_ptr, size_t bytes,
+                              MemoryType h_mt, MemoryType d_mt,
+                              bool own, bool alias, unsigned &flags)
+{
+   MFEM_CONTRACT_VAR(alias);
+   MFEM_ASSERT(exists, "Internal error!");
+   MFEM_ASSERT(!alias, "Cannot register an alias!");
+   MFEM_VERIFY_TYPES(h_mt, d_mt);
+
+   if (h_ptr == nullptr && d_ptr == nullptr)
+   {
+      MFEM_VERIFY(bytes == 0, "internal error");
+      return;
+   }
+
+   flags |= Mem::REGISTERED | Mem::OWNS_INTERNAL;
+
+   mm.InsertDevice(d_ptr, h_ptr, bytes, h_mt, d_mt);
+   flags = (own ? flags | (Mem::OWNS_HOST | Mem::OWNS_DEVICE) :
+            flags & ~(Mem::OWNS_HOST | Mem::OWNS_DEVICE)) |
+           Mem::VALID_HOST;
+
+   CheckHostMemoryType_(h_mt, h_ptr);
+}
+
 void MemoryManager::Alias_(void *base_h_ptr, size_t offset, size_t bytes,
                            unsigned base_flags, unsigned &flags)
 {
@@ -716,6 +841,34 @@ void MemoryManager::Alias_(void *base_h_ptr, size_t offset, size_t bytes,
            ~(Mem::OWNS_HOST | Mem::OWNS_DEVICE);
 }
 
+void MemoryManager::SetDeviceMemoryType_(void *h_ptr, unsigned flags,
+                                         MemoryType d_mt)
+{
+   MFEM_VERIFY(h_ptr, "cannot set the device memory type: Memory is empty!");
+   if (!(flags & Mem::ALIAS))
+   {
+      auto mem_iter = maps->memories.find(h_ptr);
+      MFEM_VERIFY(mem_iter != maps->memories.end(), "internal error");
+      internal::Memory &mem = mem_iter->second;
+      if (mem.d_mt == d_mt) { return; }
+      MFEM_VERIFY(mem.d_ptr == nullptr, "cannot set the device memory type:"
+                  " device memory is allocated!");
+      mem.d_mt = d_mt;
+   }
+   else
+   {
+      auto alias_iter = maps->aliases.find(h_ptr);
+      MFEM_VERIFY(alias_iter != maps->aliases.end(), "internal error");
+      internal::Alias &alias = alias_iter->second;
+      internal::Memory &base_mem = *alias.mem;
+      if (base_mem.d_mt == d_mt) { return; }
+      MFEM_VERIFY(base_mem.d_ptr == nullptr,
+                  "cannot set the device memory type:"
+                  " alias' base device memory is allocated!");
+      base_mem.d_mt = d_mt;
+   }
+}
+
 MemoryType MemoryManager::Delete_(void *h_ptr, MemoryType mt, unsigned flags)
 {
    const bool alias = flags & Mem::ALIAS;
@@ -749,6 +902,16 @@ MemoryType MemoryManager::Delete_(void *h_ptr, MemoryType mt, unsigned flags)
    return mt;
 }
 
+void MemoryManager::DeleteDevice_(void *h_ptr, unsigned & flags)
+{
+   const bool owns_device = flags & Mem::OWNS_DEVICE;
+   if (owns_device)
+   {
+      mm.EraseDevice(h_ptr);
+      flags = (flags | Mem::VALID_HOST) & ~Mem::VALID_DEVICE;
+   }
+}
+
 bool MemoryManager::MemoryClassCheck_(MemoryClass mc, void *h_ptr,
                                       MemoryType h_mt, size_t bytes,
                                       unsigned flags)
@@ -762,11 +925,12 @@ bool MemoryManager::MemoryClassCheck_(MemoryClass mc, void *h_ptr,
    const bool known = mm.IsKnown(h_ptr);
    const bool alias = mm.IsAlias(h_ptr);
    const bool check = known || ((flags & Mem::ALIAS) && alias);
-   MFEM_VERIFY(check,"");
+   MFEM_VERIFY(check, "Unknown host pointer: " << h_ptr);
    const internal::Memory &mem =
       (flags & Mem::ALIAS) ?
       *maps->aliases.at(h_ptr).mem : maps->memories.at(h_ptr);
-   const MemoryType &d_mt = mem.d_mt;
+   MemoryType d_mt = mem.d_mt;
+   if (d_mt == MemoryType::DEFAULT) { d_mt = GetDualMemoryType(h_mt); }
    switch (mc)
    {
       case MemoryClass::HOST_32:
@@ -785,6 +949,7 @@ bool MemoryManager::MemoryClassCheck_(MemoryClass mc, void *h_ptr,
          MFEM_VERIFY(d_mt == MemoryType::DEVICE ||
                      d_mt == MemoryType::DEVICE_DEBUG ||
                      d_mt == MemoryType::DEVICE_UMPIRE ||
+                     d_mt == MemoryType::DEVICE_UMPIRE_2 ||
                      d_mt == MemoryType::MANAGED,"");
          return true;
       }
@@ -866,7 +1031,6 @@ void *MemoryManager::Write_(void *h_ptr, MemoryType h_mt, MemoryClass mc,
       if (flags & Mem::ALIAS)
       { return mm.GetAliasDevicePtr(h_ptr, bytes, false); }
       else { return mm.GetDevicePtr(h_ptr, bytes, false); }
-
    }
 }
 
@@ -971,11 +1135,8 @@ void MemoryManager::Copy_(void *dst_h_ptr, const void *src_h_ptr,
       {
          if (dst_h_ptr != src_d_ptr && bytes != 0)
          {
-            internal::Memory &dst_h_base = maps->memories.at(dst_h_ptr);
             internal::Memory &src_d_base = maps->memories.at(src_d_ptr);
-            MemoryType dst_h_mt = dst_h_base.h_mt;
             MemoryType src_d_mt = src_d_base.d_mt;
-            ctrl->Host(dst_h_mt)->Unprotect(dst_h_ptr, bytes);
             ctrl->Device(src_d_mt)->DtoH(dst_h_ptr, src_d_ptr, bytes);
          }
       }
@@ -1100,7 +1261,7 @@ void MemoryManager::Insert(void *h_ptr, size_t bytes,
 void MemoryManager::InsertDevice(void *d_ptr, void *h_ptr, size_t bytes,
                                  MemoryType h_mt, MemoryType d_mt)
 {
-   MFEM_VERIFY_TYPES(h_mt, d_mt);
+   // MFEM_VERIFY_TYPES(h_mt, d_mt); // done by Insert() below
    MFEM_ASSERT(h_ptr != NULL, "internal error");
    Insert(h_ptr, bytes, h_mt, d_mt);
    internal::Memory &mem = maps->memories.at(h_ptr);
@@ -1153,6 +1314,20 @@ void MemoryManager::Erase(void *h_ptr, bool free_dev_ptr)
    maps->memories.erase(mem_map_iter);
 }
 
+void MemoryManager::EraseDevice(void *h_ptr)
+{
+   if (!h_ptr) { return; }
+   auto mem_map_iter = maps->memories.find(h_ptr);
+   if (mem_map_iter == maps->memories.end()) { mfem_error("Unknown pointer!"); }
+   if (maps->aliases.find(h_ptr) != maps->aliases.end())
+   {
+      mfem_error("cannot delete aliased obj!");
+   }
+   internal::Memory &mem = mem_map_iter->second;
+   if (mem.d_ptr) { ctrl->Device(mem.d_mt)->Dealloc(mem);}
+   mem.d_ptr = nullptr;
+}
+
 void MemoryManager::EraseAlias(void *alias_ptr)
 {
    if (!alias_ptr) { return; }
@@ -1173,16 +1348,21 @@ void *MemoryManager::GetDevicePtr(const void *h_ptr, size_t bytes,
    }
    internal::Memory &mem = maps->memories.at(h_ptr);
    const MemoryType &h_mt = mem.h_mt;
-   const MemoryType &d_mt = mem.d_mt;
+   MemoryType &d_mt = mem.d_mt;
    MFEM_VERIFY_TYPES(h_mt, d_mt);
-   if (!mem.d_ptr) { ctrl->Device(d_mt)->Alloc(mem); }
+   if (!mem.d_ptr)
+   {
+      if (d_mt == MemoryType::DEFAULT) { d_mt = GetDualMemoryType(h_mt); }
+      ctrl->Device(d_mt)->Alloc(mem);
+   }
+   // Aliases might have done some protections
    ctrl->Device(d_mt)->Unprotect(mem);
    if (copy_data)
    {
       MFEM_ASSERT(bytes <= mem.bytes, "invalid copy size");
       ctrl->Device(d_mt)->HtoD(mem.d_ptr, h_ptr, bytes);
    }
-   ctrl->Host(h_mt)->Protect(h_ptr, bytes);
+   ctrl->Host(h_mt)->Protect(mem, bytes);
    return mem.d_ptr;
 }
 
@@ -1201,13 +1381,18 @@ void *MemoryManager::GetAliasDevicePtr(const void *alias_ptr, size_t bytes,
    const size_t offset = alias.offset;
    internal::Memory &mem = *alias.mem;
    const MemoryType &h_mt = mem.h_mt;
-   const MemoryType &d_mt = mem.d_mt;
+   MemoryType &d_mt = mem.d_mt;
    MFEM_VERIFY_TYPES(h_mt, d_mt);
-   if (!mem.d_ptr) { ctrl->Device(d_mt)->Alloc(mem); }
+   if (!mem.d_ptr)
+   {
+      if (d_mt == MemoryType::DEFAULT) { d_mt = GetDualMemoryType(h_mt); }
+      ctrl->Device(d_mt)->Alloc(mem);
+   }
    void *alias_h_ptr = static_cast<char*>(mem.h_ptr) + offset;
    void *alias_d_ptr = static_cast<char*>(mem.d_ptr) + offset;
    MFEM_ASSERT(alias_h_ptr == alias_ptr, "internal error");
    MFEM_ASSERT(bytes <= alias.bytes, "internal error");
+   mem.d_rw = false;
    ctrl->Device(d_mt)->AliasUnprotect(alias_d_ptr, bytes);
    ctrl->Host(h_mt)->AliasUnprotect(alias_ptr, bytes);
    if (copy) { ctrl->Device(d_mt)->HtoD(alias_d_ptr, alias_h_ptr, bytes); }
@@ -1223,8 +1408,8 @@ void *MemoryManager::GetHostPtr(const void *ptr, size_t bytes, bool copy)
    const MemoryType &h_mt = mem.h_mt;
    const MemoryType &d_mt = mem.d_mt;
    MFEM_VERIFY_TYPES(h_mt, d_mt);
-   ctrl->Host(h_mt)->Unprotect(mem.h_ptr, bytes);
    // Aliases might have done some protections
+   ctrl->Host(h_mt)->Unprotect(mem, bytes);
    if (mem.d_ptr) { ctrl->Device(d_mt)->Unprotect(mem); }
    if (copy && mem.d_ptr) { ctrl->Device(d_mt)->DtoH(mem.h_ptr, mem.d_ptr, bytes); }
    if (mem.d_ptr) { ctrl->Device(d_mt)->Protect(mem); }
@@ -1242,6 +1427,7 @@ void *MemoryManager::GetAliasHostPtr(const void *ptr, size_t bytes,
    void *alias_h_ptr = static_cast<char*>(mem->h_ptr) + alias.offset;
    void *alias_d_ptr = static_cast<char*>(mem->d_ptr) + alias.offset;
    MFEM_ASSERT(alias_h_ptr == ptr,  "internal error");
+   mem->h_rw = false;
    ctrl->Host(h_mt)->AliasUnprotect(alias_h_ptr, bytes);
    if (mem->d_ptr) { ctrl->Device(d_mt)->AliasUnprotect(alias_d_ptr, bytes); }
    if (copy_data && mem->d_ptr)
@@ -1263,23 +1449,48 @@ MemoryManager::MemoryManager() { Init(); }
 
 MemoryManager::~MemoryManager() { if (exists) { Destroy(); } }
 
+void MemoryManager::SetDualMemoryType(MemoryType mt, MemoryType dual_mt)
+{
+   MFEM_VERIFY(!configured, "changing the dual MemoryTypes is not allowed after"
+               " MemoryManager configuration!");
+   UpdateDualMemoryType(mt, dual_mt);
+}
+
+void MemoryManager::UpdateDualMemoryType(MemoryType mt, MemoryType dual_mt)
+{
+   MFEM_VERIFY((int)mt < MemoryTypeSize,
+               "invalid MemoryType, mt = " << (int)mt);
+   MFEM_VERIFY((int)dual_mt < MemoryTypeSize,
+               "invalid dual MemoryType, dual_mt = " << (int)dual_mt);
+
+   if ((IsHostMemory(mt) && IsDeviceMemory(dual_mt)) ||
+       (IsDeviceMemory(mt) && IsHostMemory(dual_mt)))
+   {
+      dual_map[(int)mt] = dual_mt;
+   }
+   else
+   {
+      // mt + dual_mt is not a pair of host + device types: this is only allowed
+      // when mt == dual_mt and mt is a host type; in this case we do not
+      // actually update the dual
+      MFEM_VERIFY(mt == dual_mt && IsHostMemory(mt),
+                  "invalid (mt, dual_mt) pair: ("
+                  << MemoryTypeName[(int)mt] << ", "
+                  << MemoryTypeName[(int)dual_mt] << ')');
+   }
+}
+
 void MemoryManager::Configure(const MemoryType host_mt,
                               const MemoryType device_mt)
 {
+   MemoryManager::UpdateDualMemoryType(host_mt, device_mt);
+   MemoryManager::UpdateDualMemoryType(device_mt, host_mt);
    Init();
    host_mem_type = host_mt;
    device_mem_type = device_mt;
+   configured = true;
 }
 
-#ifdef MFEM_USE_UMPIRE
-void MemoryManager::SetUmpireAllocatorNames(const char *h_name,
-                                            const char *d_name)
-{
-   h_umpire_name = h_name;
-   d_umpire_name = d_name;
-}
-#endif
-
 void MemoryManager::Destroy()
 {
    MFEM_VERIFY(exists, "MemoryManager has already been destroyed!");
@@ -1295,6 +1506,7 @@ void MemoryManager::Destroy()
    host_mem_type = MemoryType::HOST;
    device_mem_type = MemoryType::HOST;
    exists = false;
+   configured = false;
 }
 
 void MemoryManager::RegisterCheck(void *ptr)
@@ -1381,18 +1593,36 @@ void MemoryManager::CheckHostMemoryType_(MemoryType h_mt, void *h_ptr)
 MemoryManager mm;
 
 bool MemoryManager::exists = false;
-
-#ifdef MFEM_USE_UMPIRE
-const char* MemoryManager::h_umpire_name = "HOST";
-const char* MemoryManager::d_umpire_name = "DEVICE";
-#endif
+bool MemoryManager::configured = false;
 
 MemoryType MemoryManager::host_mem_type = MemoryType::HOST;
 MemoryType MemoryManager::device_mem_type = MemoryType::HOST;
 
+MemoryType MemoryManager::dual_map[MemoryTypeSize] =
+{
+   /* HOST            */  MemoryType::DEVICE,
+   /* HOST_32         */  MemoryType::DEVICE,
+   /* HOST_64         */  MemoryType::DEVICE,
+   /* HOST_DEBUG      */  MemoryType::DEVICE_DEBUG,
+   /* HOST_UMPIRE     */  MemoryType::DEVICE_UMPIRE,
+   /* HOST_PINNED     */  MemoryType::DEVICE,
+   /* MANAGED         */  MemoryType::MANAGED,
+   /* DEVICE          */  MemoryType::HOST,
+   /* DEVICE_DEBUG    */  MemoryType::HOST_DEBUG,
+   /* DEVICE_UMPIRE   */  MemoryType::HOST_UMPIRE,
+   /* DEVICE_UMPIRE_2 */  MemoryType::HOST_UMPIRE
+};
+
+#ifdef MFEM_USE_UMPIRE
+const char * MemoryManager::h_umpire_name = "MFEM_HOST";
+const char * MemoryManager::d_umpire_name = "MFEM_DEVICE";
+const char * MemoryManager::d_umpire_2_name = "MFEM_DEVICE_2";
+#endif
+
+
 const char *MemoryTypeName[MemoryTypeSize] =
 {
-   "host-std", "host-32", "host-64", "host-debug", "host-umpire",
+   "host-std", "host-32", "host-64", "host-debug", "host-umpire", "host-pinned",
 #if defined(MFEM_USE_CUDA)
    "cuda-uvm",
    "cuda",
@@ -1405,11 +1635,14 @@ const char *MemoryTypeName[MemoryTypeSize] =
 #endif
    "device-debug",
 #if defined(MFEM_USE_CUDA)
-   "cuda-umpire"
+   "cuda-umpire",
+   "cuda-umpire-2",
 #elif defined(MFEM_USE_HIP)
-   "hip-umpire"
+   "hip-umpire",
+   "hip-umpire-2",
 #else
-   "device-umpire"
+   "device-umpire",
+   "device-umpire-2",
 #endif
 };
 

general/mem_manager.hpp

@@ -1,4 +1,4 @@
-// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
 // LICENSE and NOTICE for details. LLNL-CODE-806117.
 //
@@ -31,14 +31,26 @@ enum class MemoryType
    HOST_32,        ///< Host memory; aligned at 32 bytes
    HOST_64,        ///< Host memory; aligned at 64 bytes
    HOST_DEBUG,     ///< Host memory; allocated from a "host-debug" pool
-   HOST_UMPIRE,    ///< Host memory; using Umpire
+   HOST_UMPIRE,    /**< Host memory; using an Umpire allocator which can be set
+                        with MemoryManager::SetUmpireHostAllocatorName */
+   HOST_PINNED,    ///< Host memory: pinned (page-locked)
    MANAGED,        /**< Managed memory; using CUDA or HIP *MallocManaged
                         and *Free */
    DEVICE,         ///< Device memory; using CUDA or HIP *Malloc and *Free
    DEVICE_DEBUG,   /**< Pseudo-device memory; allocated on host from a
                         "device-debug" pool */
-   DEVICE_UMPIRE,  ///< Device memory; using Umpire
-   SIZE            ///< Number of host and device memory types
+   DEVICE_UMPIRE,  /**< Device memory; using an Umpire allocator which can be
+                        set with MemoryManager::SetUmpireDeviceAllocatorName */
+   DEVICE_UMPIRE_2, /**< Device memory; using a second Umpire allocator settable
+                         with MemoryManager::SetUmpireDevice2AllocatorName */
+   SIZE,           ///< Number of host and device memory types
+
+   PRESERVE,       /**< Pseudo-MemoryType used as default value for MemoryType
+                        parameters to request preservation of existing
+                        MemoryType, e.g. in copy constructors. */
+   DEFAULT         /**< Pseudo-MemoryType used as default value for MemoryType
+                        parameters to request the use of the default host or
+                        device MemoryType. */
 };
 
 /// Static casts to 'int' and sizes of some useful memory types.
@@ -58,16 +70,22 @@ extern const char *MemoryTypeName[MemoryTypeSize];
 enum class MemoryClass
 {
    HOST,    /**< Memory types: { HOST, HOST_32, HOST_64, HOST_DEBUG,
-                                 HOST_UMPIRE, MANAGED } */
+                                 HOST_UMPIRE, HOST_PINNED, MANAGED } */
    HOST_32, ///< Memory types: { HOST_32, HOST_64, HOST_DEBUG }
    HOST_64, ///< Memory types: { HOST_64, HOST_DEBUG }
-   DEVICE,  ///< Memory types: { DEVICE, DEVICE_DEBUG, DEVICE_UMPIRE, MANAGED }
+   DEVICE,  /**< Memory types: { DEVICE, DEVICE_DEBUG, DEVICE_UMPIRE,
+                                 DEVICE_UMPIRE_2, MANAGED } */
    MANAGED  ///< Memory types: { MANAGED }
 };
 
 /// Return true if the given memory type is in MemoryClass::HOST.
 inline bool IsHostMemory(MemoryType mt) { return mt <= MemoryType::MANAGED; }
-inline bool IsDeviceMemory(MemoryType mt) { return mt >= MemoryType::MANAGED; }
+
+/// Return true if the given memory type is in MemoryClass::DEVICE
+inline bool IsDeviceMemory(MemoryType mt)
+{
+   return mt >= MemoryType::MANAGED && mt < MemoryType::SIZE;
+}
 
 /// Return a suitable MemoryType for a given MemoryClass.
 MemoryType GetMemoryType(MemoryClass mc);
@@ -177,6 +195,12 @@ public:
        MemoryType is still set as valid. */
    Memory(int size, MemoryType mt) { New(size, mt); }
 
+   /** @brief Allocate memory for @a size entries with the given host MemoryType
+       @a h_mt and device MemoryType @a d_mt. */
+   /** The newly allocated memory is not initialized. The host pointer is set as
+       valid. */
+   Memory(int size, MemoryType h_mt, MemoryType d_mt) { New(size, h_mt, d_mt); }
+
    /** @brief Wrap an externally allocated host pointer, @a ptr with the current
        host memory type returned by MemoryManager::GetHostMemoryType(). */
    /** The parameter @a own determines whether @a ptr will be deleted when the
@@ -261,9 +285,24 @@ public:
    /** The newly allocated memory is not initialized, however the given
        MemoryType is still set as valid.
 
+       When @a mt is a host type, the device MemoryType will be set later, if
+       requested, using the dual type of @a mt, see
+       MemoryManager::GetDualMemoryType().
+
+       When @a mt is a device type, the host MemoryType will be set immediately
+       to be the dual of @a mt, see MemoryManager::GetDualMemoryType().
+
        @note The current memory is NOT deleted by this method. */
    inline void New(int size, MemoryType mt);
 
+   /** @brief Allocate memory for @a size entries with the given host MemoryType
+       @a h_mt and device MemoryType @a d_mt. */
+   /** The newly allocated memory is not initialized. The host pointer is set as
+       valid.
+
+       @note The current memory is NOT deleted by this method. */
+   inline void New(int size, MemoryType h_mt, MemoryType d_mt);
+
    /** @brief Wrap an externally allocated host pointer, @a ptr with the current
        host memory type returned by MemoryManager::GetHostMemoryType(). */
    /** The parameter @a own determines whether @a ptr will be deleted when the
@@ -284,21 +323,18 @@ public:
        @note The current memory is NOT deleted by this method. */
    inline void Wrap(T *ptr, int size, MemoryType mt, bool own);
 
-   /** Wrap an externally pair of allocated pointers, @a h_ptr and @ d_ptr,
+   /** Wrap an externally pair of allocated pointers, @a h_ptr and @a d_ptr,
        of the given host MemoryType @a h_mt. */
    /** The new memory object will have the device MemoryType set as valid.
 
        The given @a h_ptr and @a d_ptr must be allocated appropriately for the
-       given host MemoryType and its associated device MemoryType:
-          - MANAGED => MANAGED,
-          - HOST_DEBUG => DEVICE_DEBUG,
-          - HOST_UMPIRE => DEVICE_UMPIRE,
-          - HOST, HOST_32, HOST_64 => DEVICE.
+       given host MemoryType and its dual device MemoryType as defined by
+       MemoryManager::GetDualMemoryType().
 
        The parameter @a own determines whether both @a h_ptr and @a d_ptr will
        be deleted when the method Delete() is called.
 
-       @note Ownership can also be controled by using the folowing methods:
+       @note Ownership can also be controlled by using the following methods:
          - ClearOwnerFlags,
          - SetHostPtrOwner,
          - SetDevicePtrOwner.
@@ -312,10 +348,29 @@ public:
        @note The current memory is NOT deleted by this method. */
    inline void MakeAlias(const Memory &base, int offset, int size);
 
+   /// Set the device MemoryType to be used by the Memory object.
+   /** If the specified @a d_mt is not a device MemoryType, i.e. not one of the
+       types in MemoryClass::DEVICE, then this method will return immediately.
+
+       If the device MemoryType has been previously set to a different type and
+       the actual device memory has been allocated, this method will trigger an
+       error. This method will not perform the actual device memory allocation,
+       however, the allocation may already exist if the MemoryType is the same
+       as the current one.
+
+       If the Memory is an alias Memory, the device MemoryType of its base will
+       be updated as described above. */
+   inline void SetDeviceMemoryType(MemoryType d_mt);
+
    /** @brief Delete the owned pointers. The Memory is not reset by this method,
        i.e. it will, generally, not be Empty() after this call. */
    inline void Delete();
 
+   /** @brief Delete the device pointer, if owned. If @a copy_to_host is true
+       and the data is valid only on device, move it to host before deleting.
+       Invalidates the device memory. */
+   inline void DeleteDevice(bool copy_to_host = true);
+
    /// Array subscript operator for host memory.
    inline T &operator[](int idx);
 
@@ -468,6 +523,12 @@ private:
       static inline T *New(std::size_t size) { return new T[size]; }
    };
 #endif
+
+   // Shortcut for Alloc<new_align_bytes>::New(size)
+   static inline T *NewHOST(std::size_t size)
+   {
+      return Alloc<new_align_bytes>::New(size);
+   }
 };
 
 
@@ -495,10 +556,22 @@ private:
    /// Return true if the global memory manager instance exists.
    static bool Exists() { return exists; }
 
+   /// Array defining the dual MemoryType for each MemoryType
+   /** The dual of a host MemoryType is a device MemoryType and vice versa: the
+       dual of a device MemoryType is a host MemoryType. */
+   static MemoryType dual_map[MemoryTypeSize];
+
+   /// Update the dual memory type of @a mt to be @a dual_mt.
+   static void UpdateDualMemoryType(MemoryType mt, MemoryType dual_mt);
+
+   /// True if Configure() was called.
+   static bool configured;
+
    /// Host and device allocator names for Umpire.
 #ifdef MFEM_USE_UMPIRE
-   static const char *h_umpire_name;
-   static const char *d_umpire_name;
+   static const char * h_umpire_name;
+   static const char * d_umpire_name;
+   static const char * d_umpire_2_name;
 #endif
 
 private: // Static methods used by the Memory<T> class
@@ -508,26 +581,37 @@ private: // Static methods used by the Memory<T> class
    /// memory type, e.g. CUDA (mt will not be HOST).
    static void *New_(void *h_tmp, size_t bytes, MemoryType mt, unsigned &flags);
 
+   static void *New_(void *h_tmp, size_t bytes, MemoryType h_mt,
+                     MemoryType d_mt, unsigned valid_flags, unsigned &flags);
+
    /// Register an external pointer of the given MemoryType.
    /// Return the host pointer.
    static void *Register_(void *ptr, void *h_ptr, size_t bytes, MemoryType mt,
                           bool own, bool alias, unsigned &flags);
 
+   /// Register a pair of external host and device pointers
+   static void Register_(void *h_ptr, void *d_ptr, size_t bytes,
+                         MemoryType h_mt, MemoryType d_mt,
+                         bool own, bool alias, unsigned &flags);
+
    /// Register an alias. Note: base_h_ptr may be an alias.
    static void Alias_(void *base_h_ptr, size_t offset, size_t bytes,
                       unsigned base_flags, unsigned &flags);
 
+   static void SetDeviceMemoryType_(void *h_ptr, unsigned flags,
+                                    MemoryType d_mt);
+
    /// Un-register and free memory identified by its host pointer. Returns the
    /// memory type of the host pointer.
    static MemoryType Delete_(void *h_ptr, MemoryType mt, unsigned flags);
 
+   /// Free device memory identified by its host pointer
+   static void DeleteDevice_(void *h_ptr, unsigned & flags);
+
    /// Check if the memory types given the memory class are valid
    static bool MemoryClassCheck_(MemoryClass mc, void *h_ptr,
                                  MemoryType h_mt, size_t bytes, unsigned flags);
 
-   /// Return the dual memory type of the given one.
-   static MemoryType GetDualMemoryType_(MemoryType mt);
-
    /// Return a pointer to the memory identified by the host pointer h_ptr for
    /// access with the given MemoryClass.
    static void *ReadWrite_(void *h_ptr, MemoryType h_mt, MemoryClass mc,
@@ -595,6 +679,9 @@ private:
    /// Erase an address from the memory map, as well as all its aliases
    void Erase(void *h_ptr, bool free_dev_ptr = true);
 
+   /// Erase device memory for a given host address
+   void EraseDevice(void *h_ptr);
+
    /// Erase an alias from the aliases map
    void EraseAlias(void *alias_ptr);
 
@@ -621,15 +708,59 @@ public:
    /// Initialize the memory manager.
    void Init();
 
-   /// Configure the Memory manager with given default host and device types
-   /// This method will be called when configuring a device.
+   /// Return the dual MemoryType of the given one, @a mt.
+   /** The default dual memory types are:
+
+       memory type     | dual type
+       --------------- | ---------
+       HOST            | DEVICE
+       HOST_32         | DEVICE
+       HOST_64         | DEVICE
+       HOST_DEBUG      | DEVICE_DEBUG
+       HOST_UMPIRE     | DEVICE_UMPIRE
+       HOST_PINNED     | DEVICE
+       MANAGED         | MANAGED
+       DEVICE          | HOST
+       DEVICE_DEBUG    | HOST_DEBUG
+       DEVICE_UMPIRE   | HOST_UMPIRE
+       DEVICE_UMPIRE_2 | HOST_UMPIRE
+
+       The dual types can be modified before device configuration using the
+       method SetDualMemoryType() or by calling Device::SetMemoryTypes(). */
+   static inline MemoryType GetDualMemoryType(MemoryType mt)
+   { return dual_map[(int)mt]; }
+
+   /// Set the dual memory type of @a mt to be @a dual_mt.
+   /** This method can only be called before configuration, i.e. before calling
+       Configure(), which is typically done during Device construction.
+
+       One of the types must be a host MemoryType and the other must be a device
+       MemoryType or both types must be the same host memory type. The latter
+       case is only allowed for convenience in setting up pure host execution,
+       so the actual dual is not updated. */
+   static void SetDualMemoryType(MemoryType mt, MemoryType dual_mt);
+
+   /** @brief Configure the Memory manager with given default host and device
+       types. This method will be called when configuring a device.
+
+       The host and device MemoryType%s, @a h_mt and @a d_mt, are set to be dual
+       to each other. */
    void Configure(const MemoryType h_mt, const MemoryType d_mt);
 
 #ifdef MFEM_USE_UMPIRE
-   /// Set the host and device UMpire allocator names
-   void SetUmpireAllocatorNames(const char *h_name, const char *d_name);
-   const char *GetUmpireAllocatorHostName() { return h_umpire_name; }
-   const char *GetUmpireAllocatorDeviceName() { return d_umpire_name; }
+   /// Set the host Umpire allocator name used with MemoryType::HOST_UMPIRE
+   static void SetUmpireHostAllocatorName(const char * h_name) { h_umpire_name = h_name; }
+   /// Set the device Umpire allocator name used with MemoryType::DEVICE_UMPIRE
+   static void SetUmpireDeviceAllocatorName(const char * d_name) { d_umpire_name = d_name; }
+   /// Set the device Umpire allocator name used with MemoryType::DEVICE_UMPIRE_2
+   static void SetUmpireDevice2AllocatorName(const char * d_name) { d_umpire_2_name = d_name; }
+
+   /// Get the host Umpire allocator name used with MemoryType::HOST_UMPIRE
+   static const char * GetUmpireHostAllocatorName() { return h_umpire_name; }
+   /// Get the device Umpire allocator name used with MemoryType::DEVICE_UMPIRE
+   static const char * GetUmpireDeviceAllocatorName() { return d_umpire_name; }
+   /// Get the device Umpire allocator name used with MemoryType::DEVICE_UMPIRE_2
+   static const char * GetUmpireDevice2AllocatorName() { return d_umpire_2_name; }
 #endif
 
    /// Free all the device memories
@@ -663,7 +794,7 @@ template <typename T>
 inline void Memory<T>::Reset()
 {
    h_ptr = NULL;
-   h_mt = MemoryManager::host_mem_type;
+   h_mt = MemoryManager::GetHostMemoryType();
    capacity = 0;
    flags = 0;
 }
@@ -682,8 +813,8 @@ inline void Memory<T>::New(int size)
 {
    capacity = size;
    flags = OWNS_HOST | VALID_HOST;
-   h_mt = MemoryManager::host_mem_type;
-   h_ptr = (h_mt == MemoryType::HOST) ? Alloc<new_align_bytes>::New(size) :
+   h_mt = MemoryManager::GetHostMemoryType();
+   h_ptr = (h_mt == MemoryType::HOST) ? NewHOST(size) :
            (T*)MemoryManager::New_(nullptr, size*sizeof(T), h_mt, flags);
 }
 
@@ -694,12 +825,21 @@ inline void Memory<T>::New(int size, MemoryType mt)
    const size_t bytes = size*sizeof(T);
    const bool mt_host = mt == MemoryType::HOST;
    if (mt_host) { flags = OWNS_HOST | VALID_HOST; }
-   h_mt = IsHostMemory(mt) ? mt : MemoryManager::GetDualMemoryType_(mt);
-   T *h_tmp = (h_mt == MemoryType::HOST) ?
-              Alloc<new_align_bytes>::New(size) : nullptr;
+   h_mt = IsHostMemory(mt) ? mt : MemoryManager::GetDualMemoryType(mt);
+   T *h_tmp = (h_mt == MemoryType::HOST) ? NewHOST(size) : nullptr;
    h_ptr = (mt_host) ? h_tmp : (T*)MemoryManager::New_(h_tmp, bytes, mt, flags);
 }
 
+template <typename T>
+inline void Memory<T>::New(int size, MemoryType h_mt, MemoryType d_mt)
+{
+   capacity = size;
+   const size_t bytes = size*sizeof(T);
+   this->h_mt = h_mt;
+   T *h_tmp = (h_mt == MemoryType::HOST) ? NewHOST(size) : nullptr;
+   h_ptr = (T*)MemoryManager::New_(h_tmp, bytes, h_mt, d_mt, VALID_HOST, flags);
+}
+
 template <typename T>
 inline void Memory<T>::Wrap(T *ptr, int size, bool own)
 {
@@ -707,7 +847,7 @@ inline void Memory<T>::Wrap(T *ptr, int size, bool own)
    capacity = size;
    const size_t bytes = size*sizeof(T);
    flags = (own ? OWNS_HOST : 0) | VALID_HOST;
-   h_mt = MemoryManager::host_mem_type;
+   h_mt = MemoryManager::GetHostMemoryType();
 #ifdef MFEM_DEBUG
    if (own && MemoryManager::Exists())
    { MFEM_VERIFY(h_mt == MemoryManager::GetHostMemoryType_(h_ptr),""); }
@@ -726,14 +866,14 @@ inline void Memory<T>::Wrap(T *ptr, int size, MemoryType mt, bool own)
       h_ptr = ptr;
       if (mt == MemoryType::HOST || !own)
       {
-         // Skip restration
+         // Skip registration
          flags = (own ? OWNS_HOST : 0) | VALID_HOST;
          return;
       }
    }
    else
    {
-      h_mt = MemoryManager::GetDualMemoryType_(mt);
+      h_mt = MemoryManager::GetDualMemoryType(mt);
       h_ptr = (h_mt == MemoryType::HOST) ? new T[size] : nullptr;
    }
    flags = 0;
@@ -750,8 +890,8 @@ inline void Memory<T>::Wrap(T *ptr, T *d_ptr, int size, MemoryType mt, bool own)
    capacity = size;
    MFEM_ASSERT(IsHostMemory(h_mt),"");
    const size_t bytes = size*sizeof(T);
-   const MemoryType d_mt = MemoryManager::GetDualMemoryType_(h_mt);
-   MemoryManager::Register_(d_ptr, h_ptr, bytes, d_mt, own, false, flags);
+   const MemoryType d_mt = MemoryManager::GetDualMemoryType(h_mt);
+   MemoryManager::Register_(h_ptr, d_ptr, bytes, h_mt, d_mt, own, false, flags);
 }
 
 template <typename T>
@@ -770,6 +910,18 @@ inline void Memory<T>::MakeAlias(const Memory &base, int offset, int size)
    }
 }
 
+template <typename T>
+inline void Memory<T>::SetDeviceMemoryType(MemoryType d_mt)
+{
+   if (!IsDeviceMemory(d_mt)) { return; }
+   if (!(flags & REGISTERED))
+   {
+      MemoryManager::Register_(h_ptr, nullptr, capacity*sizeof(T), h_mt,
+                               flags & OWNS_HOST, flags & ALIAS, flags);
+   }
+   MemoryManager::SetDeviceMemoryType_(h_ptr, flags, d_mt);
+}
+
 template <typename T>
 inline void Memory<T>::Delete()
 {
@@ -784,6 +936,16 @@ inline void Memory<T>::Delete()
    }
 }
 
+template <typename T>
+inline void Memory<T>::DeleteDevice(bool copy_to_host)
+{
+   if (flags & REGISTERED)
+   {
+      if (copy_to_host) { Read(MemoryClass::HOST, capacity); }
+      MemoryManager::DeleteDevice_((void*)h_ptr, flags);
+   }
+}
+
 template <typename T>
 inline T &Memory<T>::operator[](int idx)
 {
@@ -911,7 +1073,7 @@ inline void Memory<T>::CopyFrom(const Memory &src, int size)
    {
       if (h_ptr != src.h_ptr && size != 0)
       {
-         MFEM_ASSERT(h_ptr + size <= src || src + size <= h_ptr,
+         MFEM_ASSERT(h_ptr + size <= src.h_ptr || src.h_ptr + size <= h_ptr,
                      "data overlaps!");
          std::memcpy(h_ptr, src, size*sizeof(T));
       }

linalg/densemat.cpp

@@ -3586,4 +3586,114 @@ void BatchLUSolve(const DenseTensor &Mlu, const Array<int> &P, Vector &X)
 
 }
 
+void BatchLUFactor(Vector &Minv,const int m,const int NE, Array<int> &P)
+{
+   P.SetSize(m*NE);
+   auto data_all = mfem::Reshape(Minv.ReadWrite(), m, m, NE);
+   auto piv_all = mfem::Reshape(P.Write(), m, NE);
+
+   MFEM_FORALL(e, NE,
+   {
+
+      double *data = &data_all(0,0,e);
+      int *ipiv = &piv_all(0,e);
+      for (int i = 0; i < m; i++)
+      {
+
+         // pivoting
+         {
+            int piv = i;
+            double a = fabs(data[piv+i*m]);
+            for (int j = i+1; j < m; j++)
+            {
+               const double b = fabs(data[j+i*m]);
+               if (b > a)
+               {
+                  a = b;
+                  piv = j;
+               }
+            }
+            ipiv[i] = piv;
+            if (piv != i)
+            {
+               // swap rows i and piv in both L and U parts
+               for (int j = 0; j < m; j++)
+               {
+                  mfem::kernels::internal::Swap<double>(data[i+j*m], data[piv+j*m]);
+               }
+            }
+         }//pivot end
+
+         //Q: How to check for errors?
+         //if (abs(data[i + i*m]) <= TOL)
+         //{
+         //return false; // failed
+         //}
+
+         const double a_ii_inv = 1.0 / data[i+i*m];
+         for (int j = i+1; j < m; j++)
+         {
+            data[j+i*m] *= a_ii_inv;
+         }
+
+         for (int k = i+1; k < m; k++)
+         {
+            const double a_ik = data[i+k*m];
+            for (int j = i+1; j < m; j++)
+            {
+               data[j+k*m] -= a_ik * data[j+i*m];
+            }
+         }
+
+      }//m loop
+
+   });
+
+}
+
+void BatchLUSolve(Vector &Minv, int m, int NE,
+                  Array<int> &P, Vector &X)
+{
+
+   auto data_all = mfem::Reshape(Minv.Read(), m, m, NE);
+   auto piv_all = mfem::Reshape(P.Read(), m, NE);
+   auto x_all = mfem::Reshape(X.ReadWrite(), m, NE);
+
+   MFEM_FORALL(e, NE,
+   {
+
+      const double *data = &data_all(0,0,e);
+      const int *ipiv = &piv_all(0,e);
+      double *x = &x_all(0,e);
+
+      // X <- P X
+      for (int i = 0; i < m; i++)
+      {
+         mfem::kernels::internal::Swap<double>(x[i], x[ipiv[i]]);
+      }
+
+      // X <- L^{-1} X
+      for (int j = 0; j < m; j++)
+      {
+         const double x_j = x[j];
+         for (int i = j+1; i < m; i++)
+         {
+            x[i] -= data[i+j*m] * x_j;
+         }
+      }
+
+      // X <- U^{-1} X
+      for (int j = m-1; j >= 0; j--)
+      {
+         const double x_j = ( x[j] /= data[j+j*m] );
+         for (int i = 0; i < j; i++)
+         {
+            x[i] -= data[i+j*m] * x_j;
+         }
+      }
+   });
+
+}
+
+
 } // namespace mfem

linalg/densemat.hpp

@@ -746,6 +746,13 @@ public:
       tdata.New(i*j*k);
    }
 
+   DenseTensor(int i, int j, int k, MemoryType mt)
+      : Mk(NULL, i, j)
+   {
+      nk = k;
+      tdata.New(i*j*k, mt);
+   }
+
    /// Copy constructor: deep copy
    DenseTensor(const DenseTensor &other)
       : Mk(NULL, other.Mk.height, other.Mk.width), nk(other.nk)
@@ -768,9 +775,9 @@ public:
 
    int TotalSize() const { return SizeI()*SizeJ()*SizeK(); }
 
-   void SetSize(int i, int j, int k)
+   void SetSize(int i, int j, int k, MemoryType mt_ = MemoryType::SIZE)
    {
-      const MemoryType mt = tdata.GetMemoryType();
+      const MemoryType mt = mt_ == MemoryType::SIZE ? tdata.GetMemoryType() : mt_;
       tdata.Delete();
       Mk.UseExternalData(NULL, i, j);
       nk = k;
@@ -884,6 +891,10 @@ void BatchLUFactor(DenseTensor &Mlu, Array<int> &P, const double TOL = 0.0);
     dimension m x n. */
 void BatchLUSolve(const DenseTensor &Mlu, const Array<int> &P, Vector &X);
 
+void BatchLUFactor(Vector &Minv,int m,int NE, Array<int> &P);
+
+void BatchLUSolve(Vector &Minv, int m, int NE,
+                  Array<int> &P, Vector &X);
 
 // Inline methods
 

linalg/operator.cpp

@@ -597,7 +597,10 @@ double PowerMethod::EstimateLargestEigenvalue(Operator& opr, Vector& v0,
                                               int numSteps, double tolerance, int seed)
 {
    v1.SetSize(v0.Size());
-   v0.Randomize(seed);
+   if (seed != -1)
+   {
+      v0.Randomize(seed);
+   }
 
    double eigenvalue = 1.0;
 

linalg/operator.hpp

@@ -931,7 +931,8 @@ public:
        the eigenvector corresponding to the largest eigenvalue after convergence.
        The maximum number of iterations may set with \p numSteps, the relative
        tolerance with \p tolerance and the seed of the random initialization of
-       \p v0 with \p seed. */
+       \p v0 with \p seed, unless \p is -1, in which case \p v0 will be used
+       as-is. */
    double EstimateLargestEigenvalue(Operator& opr, Vector& v0,
                                     int numSteps = 10, double tolerance = 1e-8,
                                     int seed = 12345);

linalg/solvers.cpp

@@ -178,7 +178,11 @@ void OperatorJacobiSmoother::Setup(const Vector &diag)
    const double delta = damping;
    auto D = diag.Read();
    auto DI = dinv.Write();
-   MFEM_FORALL(i, N, DI[i] = delta / D[i]; );
+   MFEM_FORALL(i, N,
+   {
+      const double d_i = (D[i] < 0.0 && abs_values) ? -D[i] : D[i];
+      DI[i] = delta / d_i;
+   });
    auto I = ess_tdof_list.Read();
    MFEM_FORALL(i, ess_tdof_list.Size(), DI[I[i]] = delta; );
 }

linalg/solvers.hpp

@@ -139,6 +139,7 @@ public:
 
    void SetOperator(const Operator &op);
    void Setup(const Vector &diag);
+   void SetAbsValues(bool use_abs_values) { abs_values = use_abs_values; }
 
 private:
    const int N;
@@ -147,6 +148,7 @@ private:
    const Array<int> &ess_tdof_list;
    mutable Vector residual;
    const bool dynamic = false;
+   bool abs_values = false;
 
    const Operator *oper;
    const NonlinearForm *nlf = NULL;

linalg/sparsemat.cpp

@@ -28,6 +28,38 @@ namespace mfem
 
 using namespace std;
 
+#ifdef MFEM_USE_CUDA
+int SparseMatrix::SparseMatrixCount = 0;
+cusparseHandle_t SparseMatrix::handle;
+size_t SparseMatrix::bufferSize = 0;
+void * SparseMatrix::dBuffer = nullptr;
+#endif
+
+void SparseMatrix::InitCuSparse()
+{
+   /* Initialize CuSparse library */
+#ifdef MFEM_USE_CUDA
+   SparseMatrixCount++;
+   if (SparseMatrixCount == 1 && Device::Allows(Backend::CUDA_MASK))
+   {
+      cusparseCreate(&handle);
+   }
+#endif
+}
+
+void SparseMatrix::ClearCuSparseDescr()
+{
+#ifdef MFEM_USE_CUDA
+   if (initBuffers)
+   {
+      cusparseDestroySpMat(matA_descr);
+      cusparseDestroyDnVec(vecX_descr);
+      cusparseDestroyDnVec(vecY_descr);
+      initBuffers = false;
+   }
+#endif
+}
+
 SparseMatrix::SparseMatrix(int nrows, int ncols)
    : AbstractSparseMatrix(nrows, (ncols >= 0) ? ncols : nrows),
      Rows(new RowNode *[nrows]),
@@ -50,6 +82,8 @@ SparseMatrix::SparseMatrix(int nrows, int ncols)
 #ifdef MFEM_USE_MEMALLOC
    NodesMem = new RowNodeAlloc;
 #endif
+
+   InitCuSparse();
 }
 
 SparseMatrix::SparseMatrix(int *i, int *j, double *data, int m, int n)
@@ -67,6 +101,8 @@ SparseMatrix::SparseMatrix(int *i, int *j, double *data, int m, int n)
 #ifdef MFEM_USE_MEMALLOC
    NodesMem = NULL;
 #endif
+
+   InitCuSparse();
 }
 
 SparseMatrix::SparseMatrix(int *i, int *j, double *data, int m, int n,
@@ -98,6 +134,8 @@ SparseMatrix::SparseMatrix(int *i, int *j, double *data, int m, int n,
          A[i] = 0.0;
       }
    }
+
+   InitCuSparse();
 }
 
 SparseMatrix::SparseMatrix(int nrows, int ncols, int rowsize)
@@ -119,9 +157,12 @@ SparseMatrix::SparseMatrix(int nrows, int ncols, int rowsize)
    {
       I[i] = i * rowsize;
    }
+
+   InitCuSparse();
 }
 
-SparseMatrix::SparseMatrix(const SparseMatrix &mat, bool copy_graph)
+SparseMatrix::SparseMatrix(const SparseMatrix &mat, bool copy_graph,
+                           MemoryType mt)
    : AbstractSparseMatrix(mat.Height(), mat.Width())
 {
    if (mat.Finalized())
@@ -129,8 +170,8 @@ SparseMatrix::SparseMatrix(const SparseMatrix &mat, bool copy_graph)
       const int nnz = mat.I[height];
       if (copy_graph)
       {
-         I.New(height+1, mat.I.GetMemoryType());
-         J.New(nnz, mat.J.GetMemoryType());
+         I.New(height+1, mt == MemoryType::SIZE ? mat.I.GetMemoryType() : mt);
+         J.New(nnz, mt == MemoryType::SIZE ? mat.J.GetMemoryType() : mt);
          I.CopyFrom(mat.I, height+1);
          J.CopyFrom(mat.J, nnz);
       }
@@ -141,7 +182,7 @@ SparseMatrix::SparseMatrix(const SparseMatrix &mat, bool copy_graph)
          I.ClearOwnerFlags();
          J.ClearOwnerFlags();
       }
-      A.New(nnz, mat.A.GetMemoryType());
+      A.New(nnz, mt == MemoryType::SIZE ? mat.A.GetMemoryType() : mt);
       A.CopyFrom(mat.A, nnz);
 
       Rows = NULL;
@@ -174,9 +215,21 @@ SparseMatrix::SparseMatrix(const SparseMatrix &mat, bool copy_graph)
       }
 
       // We probably do not need to set the ownership flags here.
-      I.Reset(); I.SetHostPtrOwner(true);
-      J.Reset(); J.SetHostPtrOwner(true);
-      A.Reset(); A.SetHostPtrOwner(true);
+      if (mt != MemoryType::SIZE)
+      {
+         I.Reset(mt);
+         J.Reset(mt);
+         A.Reset(mt);
+      }
+      else
+      {
+         I.Reset();
+         J.Reset();
+         A.Reset();
+      }
+      I.SetHostPtrOwner(true);
+      J.SetHostPtrOwner(true);
+      A.SetHostPtrOwner(true);
    }
 
    current_row = -1;
@@ -184,6 +237,8 @@ SparseMatrix::SparseMatrix(const SparseMatrix &mat, bool copy_graph)
    ColPtrNode = NULL;
    At = NULL;
    isSorted = mat.isSorted;
+
+   InitCuSparse();
 }
 
 SparseMatrix::SparseMatrix(const Vector &v)
@@ -211,6 +266,8 @@ SparseMatrix::SparseMatrix(const Vector &v)
       J[r] = r;
       A[r] = v[r];
    }
+
+   InitCuSparse();
 }
 
 SparseMatrix& SparseMatrix::operator=(const SparseMatrix &rhs)
@@ -250,6 +307,8 @@ void SparseMatrix::SetEmpty()
    NodesMem = NULL;
 #endif
    isSorted = false;
+
+   ClearCuSparseDescr();
 }
 
 int SparseMatrix::RowSize(const int i) const
@@ -592,16 +651,72 @@ void SparseMatrix::AddMult(const Vector &x, Vector &y, const double a) const
    auto d_A = Read(A, nnz);
    auto d_x = x.Read();
    auto d_y = y.ReadWrite();
-   MFEM_FORALL(i, height,
+
+   //Skip if matrix has no non-zeros
+   if (nnz == 0) {return;}
+   if (Device::Allows(Backend::CUDA_MASK) && useCuSparse)
    {
-      double d = 0.0;
-      const int end = d_I[i+1];
-      for (int j = d_I[i]; j < end; j++)
+#ifdef MFEM_USE_CUDA
+      const double alpha = a;
+      const double beta  = 1.0;
+
+      //Setup descriptors
+      if (!initBuffers)
       {
-         d += d_A[j] * d_x[d_J[j]];
+         /* Setup matrix descriptor */
+         cusparseCreateCsr(&matA_descr,Height(), Width(), J.Capacity(),
+                           const_cast<int *>(d_I),
+                           const_cast<int *>(d_J), const_cast<double *>(d_A), CUSPARSE_INDEX_32I,
+                           CUSPARSE_INDEX_32I, CUSPARSE_INDEX_BASE_ZERO, CUDA_R_64F);
+
+         /*Create handles for input/output vectors */
+         cusparseCreateDnVec(&vecX_descr, x.Size(), const_cast<double *>(d_x),
+                             CUDA_R_64F);
+         cusparseCreateDnVec(&vecY_descr, y.Size(), d_y, CUDA_R_64F);
+
+         initBuffers = true;
       }
-      d_y[i] += a * d;
-   });
+
+      /*Allocate space for kernel. Buffer is shared between different sparsemats */
+      size_t newBufferSize = 0;
+      cusparseSpMV_bufferSize(handle, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha,
+                              matA_descr,
+                              vecX_descr, &beta, vecY_descr, CUDA_R_64F,
+                              CUSPARSE_CSRMV_ALG1, &newBufferSize);
+
+      //Check if need to resize
+      if (newBufferSize > bufferSize)
+      {
+         bufferSize = newBufferSize;
+         if (dBuffer != NULL) { CuMemFree(dBuffer); }
+         CuMemAlloc(&dBuffer, bufferSize);
+      }
+
+      //Update input/output vectors
+      cusparseDnVecSetValues(vecX_descr, const_cast<double *>(d_x));
+      cusparseDnVecSetValues(vecY_descr, d_y);
+
+      // Y = alpha A * X + beta * Y
+      cusparseSpMV(handle, CUSPARSE_OPERATION_NON_TRANSPOSE, &alpha, matA_descr,
+                   vecX_descr, &beta, vecY_descr, CUDA_R_64F, CUSPARSE_CSRMV_ALG1, dBuffer);
+#endif
+   }
+   else
+   {
+      //Native version
+      MFEM_FORALL(i, height,
+      {
+         double d = 0.0;
+         const int end = d_I[i+1];
+         for (int j = d_I[i]; j < end; j++)
+         {
+            d += d_A[j] * d_x[d_J[j]];
+         }
+         d_y[i] += a * d;
+      });
+
+   }
+
 #else
    const double *Ap = A, *xp = x.GetData();
    double *yp = y.GetData();
@@ -2956,6 +3071,10 @@ void SparseMatrix::Destroy()
    delete NodesMem;
 #endif
    delete At;
+
+#ifdef MFEM_USE_CUDA
+   ClearCuSparseDescr();
+#endif
 }
 
 int SparseMatrix::ActualWidth() const

linalg/sparsemat.hpp

@@ -21,6 +21,12 @@
 #include "../general/globals.hpp"
 #include "densemat.hpp"
 
+#ifdef MFEM_USE_CUDA
+#include <cusparse.h>
+#include <library_types.h>
+#include "../general/cuda.hpp"
+#endif
+
 namespace mfem
 {
 
@@ -80,9 +86,34 @@ protected:
    void Destroy();   // Delete all owned data
    void SetEmpty();  // Init all entries with empty values
 
+   bool useCuSparse{true}; //Use CuSparse if available
+
+   // Initialize CuSparse
+   void InitCuSparse();
+
+
+#ifdef MFEM_USE_CUDA
+   cusparseStatus_t status;
+   static cusparseHandle_t handle;
+   cusparseMatDescr_t descr=0;
+   static size_t bufferSize;
+   static void *dBuffer;
+   mutable bool initBuffers{false};
+
+   static int SparseMatrixCount;
+   mutable cusparseSpMatDescr_t matA_descr;
+   mutable cusparseDnVecDescr_t vecX_descr;
+   mutable cusparseDnVecDescr_t vecY_descr;
+#endif
+
 public:
    /// Create an empty SparseMatrix.
-   SparseMatrix() { SetEmpty(); }
+   SparseMatrix()
+   {
+      SetEmpty();
+
+      InitCuSparse();
+   }
 
    /** @brief Create a sparse matrix with flexible sparsity structure using a
        row-wise linked list (LIL) format. */
@@ -113,11 +144,15 @@ public:
    /** If @a mat is finalized and @a copy_graph is false, the #I and #J arrays
        will use a shallow copy (copy the pointers only) without transferring
        ownership. */
-   SparseMatrix(const SparseMatrix &mat, bool copy_graph = true);
+   SparseMatrix(const SparseMatrix &mat, bool copy_graph = true,
+                MemoryType mt = MemoryType::SIZE);
 
    /// Create a SparseMatrix with diagonal @a v, i.e. A = Diag(v)
    SparseMatrix(const Vector & v);
 
+   // Runtime option to use CuSparse
+   // Only valid when using a CUDA backend
+   void UseCuSparse(bool _useCuSparse = true) { useCuSparse = _useCuSparse;}
 
    /// Assignment operator: deep copy
    SparseMatrix& operator=(const SparseMatrix &rhs);
@@ -134,6 +169,9 @@ public:
    /// Clear the contents of the SparseMatrix.
    void Clear() { Destroy(); SetEmpty(); }
 
+   /// Release the cusparse descriptors.
+   void ClearCuSparseDescr();
+
    /// Check if the SparseMatrix is empty.
    bool Empty() const { return (A == NULL) && (Rows == NULL); }
 
@@ -573,7 +611,18 @@ public:
    void Swap(SparseMatrix &other);
 
    /// Destroys sparse matrix.
-   virtual ~SparseMatrix() { Destroy(); }
+   virtual ~SparseMatrix()
+   {
+      Destroy();
+#ifdef MFEM_USE_CUDA
+      if (handle && SparseMatrixCount==1 && Device::Allows(Backend::CUDA_MASK))
+      {
+         cusparseDestroy(handle);
+         CuMemFree(dBuffer);
+      }
+      SparseMatrixCount--;
+#endif
+   }
 
    Type GetType() const { return MFEM_SPARSEMAT; }
 };

makefile

@@ -261,7 +261,7 @@ endif
 # List of MFEM dependencies, that require the *_LIB variable to be non-empty
 MFEM_REQ_LIB_DEPS = SUPERLU METIS CONDUIT SIDRE LAPACK SUNDIALS MESQUITE\
  SUITESPARSE STRUMPACK GINKGO GNUTLS NETCDF PETSC SLEPC MPFR PUMI HIOP GSLIB\
- OCCA CEED RAJA UMPIRE
+ OCCA CEED RAJA UMPIRE CAMP
 PETSC_ERROR_MSG = $(if $(PETSC_FOUND),,. PETSC config not found: $(PETSC_VARS))
 SLEPC_ERROR_MSG = $(if $(SLEPC_FOUND),,. SLEPC config not found: $(SLEPC_VARS))
 
@@ -324,8 +324,8 @@ MFEM_DEFINES = MFEM_VERSION MFEM_VERSION_STRING MFEM_GIT_STRING MFEM_USE_MPI\
  MFEM_USE_SUPERLU MFEM_USE_STRUMPACK MFEM_USE_GNUTLS\
  MFEM_USE_NETCDF MFEM_USE_PETSC MFEM_USE_SLEPC MFEM_USE_MPFR MFEM_USE_SIDRE MFEM_USE_CONDUIT\
  MFEM_USE_PUMI MFEM_USE_HIOP MFEM_USE_GSLIB MFEM_USE_CUDA MFEM_USE_HIP\
- MFEM_USE_OCCA MFEM_USE_CEED MFEM_USE_RAJA MFEM_USE_UMPIRE MFEM_USE_SIMD\
- MFEM_USE_ADIOS2 MFEM_SOURCE_DIR MFEM_INSTALL_DIR
+ MFEM_USE_OCCA MFEM_USE_CEED MFEM_USE_RAJA MFEM_USE_UMPIRE MFEM_USE_CAMP\
+ MFEM_USE_SIMD MFEM_USE_ADIOS2 MFEM_SOURCE_DIR MFEM_INSTALL_DIR
 
 # List of makefile variables that will be written to config.mk:
 MFEM_CONFIG_VARS = MFEM_CXX MFEM_HOST_CXX MFEM_CPPFLAGS MFEM_CXXFLAGS\
@@ -645,6 +645,7 @@ status info:
 	$(info MFEM_USE_OCCA          = $(MFEM_USE_OCCA))
 	$(info MFEM_USE_CEED          = $(MFEM_USE_CEED))
 	$(info MFEM_USE_UMPIRE        = $(MFEM_USE_UMPIRE))
+	$(info MFEM_USE_CAMP          = $(MFEM_USE_CAMP))
 	$(info MFEM_USE_SIMD          = $(MFEM_USE_SIMD))
 	$(info MFEM_USE_ADIOS2        = $(MFEM_USE_ADIOS2))
 	$(info MFEM_CXX               = $(value MFEM_CXX))

mesh/mesh.cpp

@@ -755,7 +755,9 @@ void Mesh::GetLocalQuadToWdgTransformation(
 
 const GeometricFactors* Mesh::GetGeometricFactors(const IntegrationRule& ir,
                                                   const int flags,
-                                                  mfem::DofToQuad::Mode mode)
+                                                  mfem::DofToQuad::Mode mode,
+                                                  MemoryType d_mt
+                                               )
 {
    for (int i = 0; i < geom_factors.Size(); i++)
    {
@@ -767,7 +769,7 @@ const GeometricFactors* Mesh::GetGeometricFactors(const IntegrationRule& ir,
    }
 
    this->EnsureNodes();
-   GeometricFactors *gf = new GeometricFactors(this, ir, flags, mode);
+   GeometricFactors *gf = new GeometricFactors(this, ir, flags, mode, d_mt);
    geom_factors.Append(gf);
    return gf;
 }
@@ -5825,6 +5827,7 @@ int *Mesh::GeneratePartitioning(int nparts, int part_method)
    el_to_el = NULL;
 
    // Check for empty partitionings (a "feature" in METIS)
+   if (nparts > 1 && NumOfElements > nparts)
    {
       Array< Pair<int,int> > psize(nparts);
       int empty_parts;
@@ -10511,7 +10514,7 @@ int Mesh::FindPoints(DenseMatrix &point_mat, Array<int>& elem_ids,
 
 
 GeometricFactors::GeometricFactors(const Mesh *mesh, const IntegrationRule &ir,
-                                   int flags, DofToQuad::Mode mode)
+                                   int flags, DofToQuad::Mode mode, MemoryType d_mt)
 {
    this->mesh = mesh;
    IntRule = &ir;
@@ -10527,19 +10530,21 @@ GeometricFactors::GeometricFactors(const Mesh *mesh, const IntegrationRule &ir,
    const int NQ   = ir.GetNPoints();
 
    unsigned eval_flags = 0;
+   MemoryType my_d_mt = (d_mt != MemoryType::DEFAULT) ? d_mt :
+                         Device::GetDeviceMemoryType();
    if (flags & GeometricFactors::COORDINATES)
    {
-      X.SetSize(vdim*NQ*NE);
+      X.SetSize(vdim*NQ*NE, my_d_mt);
       eval_flags |= QuadratureInterpolator::VALUES;
    }
    if (flags & GeometricFactors::JACOBIANS)
    {
-      J.SetSize(dim*vdim*NQ*NE);
+      J.SetSize(dim*vdim*NQ*NE, my_d_mt);
       eval_flags |= QuadratureInterpolator::DERIVATIVES;
    }
    if (flags & GeometricFactors::DETERMINANTS)
    {
-      detJ.SetSize(NQ*NE);
+      detJ.SetSize(NQ*NE, my_d_mt);
       eval_flags |= QuadratureInterpolator::DETERMINANTS;
    }
 
@@ -10558,7 +10563,67 @@ GeometricFactors::GeometricFactors(const Mesh *mesh, const IntegrationRule &ir,
 
    if (elem_restr)
    {
-      Vector Enodes(vdim*ND*NE);
+      Vector Enodes(vdim*ND*NE, my_d_mt);
+      elem_restr->Mult(*nodes, Enodes);
+      qi->Mult(Enodes, eval_flags, X, J, detJ);
+   }
+   else
+   {
+      qi->Mult(*nodes, eval_flags, X, J, detJ);
+   }
+}
+
+GeometricFactors::GeometricFactors(const GridFunction *nodes_,
+                                   const IntegrationRule &ir,
+                                   int flags, DofToQuad::Mode mode)
+{
+   this->mesh = nullptr;
+   IntRule = &ir;
+   computed_factors = flags;
+
+   const GridFunction *nodes = nodes_;
+   const FiniteElementSpace *fespace = nodes->FESpace();
+   const FiniteElement *fe = fespace->GetFE(0);
+   const int dim  = fe->GetDim();
+   const int vdim = fespace->GetVDim();
+   const int NE   = fespace->GetNE();
+   const int ND   = fe->GetDof();
+   const int NQ   = ir.GetNPoints();
+
+   unsigned eval_flags = 0;
+   MemoryType d_mt = Device::GetDeviceMemoryType();
+   if (flags & GeometricFactors::COORDINATES)
+   {
+      X.SetSize(vdim*NQ*NE, d_mt);
+      eval_flags |= QuadratureInterpolator::VALUES;
+   }
+   if (flags & GeometricFactors::JACOBIANS)
+   {
+      J.SetSize(dim*vdim*NQ*NE, d_mt);
+      eval_flags |= QuadratureInterpolator::DERIVATIVES;
+   }
+   if (flags & GeometricFactors::DETERMINANTS)
+   {
+      detJ.SetSize(NQ*NE, d_mt);
+      eval_flags |= QuadratureInterpolator::DETERMINANTS;
+   }
+
+   const QuadratureInterpolator *qi =
+      fespace->GetQuadratureInterpolator(ir, mode);
+   const bool use_tensor_products = qi->UseTensorProducts();
+
+   // GeometricFactors arrays use a column-major layout
+   qi->SetOutputLayout(QVectorLayout::byNODES);
+
+   // Use LEXICOGRAPHIC ordering in case of tensor product evaluation
+   const ElementDofOrdering e_ordering = use_tensor_products ?
+                                         ElementDofOrdering::LEXICOGRAPHIC :
+                                         ElementDofOrdering::NATIVE;
+   const Operator *elem_restr = fespace->GetElementRestriction(e_ordering);
+
+   if (elem_restr)
+   {
+      Vector Enodes(vdim*ND*NE, d_mt);
       elem_restr->Mult(*nodes, Enodes);
       qi->Mult(Enodes, eval_flags, X, J, detJ);
    }

mesh/mesh.hpp

@@ -743,7 +743,8 @@ public:
        integration rule. */
    const GeometricFactors* GetGeometricFactors(const IntegrationRule& ir,
                                                const int flags,
-                                               DofToQuad::Mode = DofToQuad::FULL);
+                                               DofToQuad::Mode = DofToQuad::FULL,
+                                               MemoryType d_mt = MemoryType::DEFAULT);
 
    /** @brief Return the mesh geometric factors for the faces corresponding
         to the given integration rule. */
@@ -1374,6 +1375,11 @@ public:
    };
 
    GeometricFactors(const Mesh *mesh, const IntegrationRule &ir, int flags,
+                    DofToQuad::Mode = DofToQuad::FULL,
+                    MemoryType d_mt = MemoryType::DEFAULT);
+
+   GeometricFactors(const GridFunction *nodes_, const IntegrationRule &ir,
+                    int flags,
                     DofToQuad::Mode = DofToQuad::FULL);
 
    /// Mapped (physical) coordinates of all quadrature points.

mesh/pmesh.cpp

@@ -2394,9 +2394,18 @@ void ParMesh::GetGhostFaceTransformation(
 }
 
 FaceElementTransformations *ParMesh::
-GetSharedFaceTransformations(int sf, bool fill2)
+GetSharedFaceTransformations(int sf, bool fill2, bool direct)
 {
-   int FaceNo = GetSharedFace(sf);
+   //int FaceNo = GetSharedFace(sf);
+   int FaceNo;
+   if (direct)
+   {
+      FaceNo = sf;
+   }
+   else
+   {
+      FaceNo = GetSharedFace(sf);
+   }
 
    FaceInfo &face_info = faces_info[FaceNo];
 

mesh/pmesh.hpp

@@ -293,9 +293,9 @@ public:
 
    /** Get the FaceElementTransformations for the given shared face (edge 2D).
        In the returned object, 1 and 2 refer to the local and the neighbor
-       elements, respectively. */
+       elements, respectively. Use direct if sf is the face number */
    FaceElementTransformations *
-   GetSharedFaceTransformations(int sf, bool fill2 = true);
+   GetSharedFaceTransformations(int sf, bool fill2 = true, bool direct = false);
 
    ElementTransformation *
    GetFaceNbrElementTransformation(int i)

miniapps/meshing/mesh-optimizer.hpp

@@ -239,6 +239,15 @@ double adapt_lim_fun(const Vector &x)
    return val;
 }
 
+// Space-dependent limiting.
+double limit_func(const Vector &x)
+{
+   //if (x(0) < x(1) ) { return 1.0e-6; }
+   //else return 10000;
+
+   return 1.0;
+}
+
 void DiffuseField(GridFunction &field, int smooth_steps)
 {
    //Setup the Laplacian operator

miniapps/meshing/pmesh-optimizer.cpp

@@ -633,8 +633,13 @@ int main (int argc, char *argv[])
 
    // 14. Limit the node movement.
    // The limiting distances can be given by a general function of space.
-   ParGridFunction dist(pfespace);
-   dist = 1.0;
+   // TODO come up with proper example for space-dependent limiting that shows
+   //   differences between Bernstein and Gauss-Lobatto.
+   H1_FECollection lim_fec(mesh_poly_deg, dim, BasisType::Positive);
+   ParFiniteElementSpace lim_fes(pmesh, &lim_fec);
+   ParGridFunction dist(&lim_fes);
+   FunctionCoefficient lim_dist_coeff(limit_func);
+   dist.ProjectCoefficient(lim_dist_coeff);
    // The small_phys_size is relevant only with proper normalization.
    if (normalization) { dist = small_phys_size; }
    ConstantCoefficient lim_coeff(lim_const);
@@ -806,13 +811,17 @@ int main (int argc, char *argv[])
       minres->SetMaxIter(max_lin_iter);
       minres->SetRelTol(linsol_rtol);
       minres->SetAbsTol(0.0);
-      minres->SetPrintLevel(verbosity_level >= 2 ? 3 : -1);
+      if (verbosity_level > 2) { minres->SetPrintLevel(1); }
+      else { minres->SetPrintLevel(verbosity_level = 2 ? 3 : -1); }
       if (lin_solver == 3 || lin_solver == 4)
       {
          if (pa)
          {
             MFEM_VERIFY(lin_solver != 4, "PA l1-Jacobi is not implemented");
-            S_prec = new OperatorJacobiSmoother(a, a.GetEssentialTrueDofs());
+            OperatorJacobiSmoother *sm =
+               new OperatorJacobiSmoother(a, a.GetEssentialTrueDofs());
+            sm->SetAbsValues(true);
+            S_prec = sm;
          }
          else
          {