add optional MemoryType to the SparseMatrix copy-ctor

Cherry-picking changes on top of our current hash.
There are things that hurt performance in laster mfem@master commits, not sure what yet. This includes: - feature/artv3/cusparse-Spmv - feature/tomstitt/temp-mem-type - patches to tmop.cpp and tmop_tools.cpp to address memory issues - patch to mem_manager/device to use mfem's default allcator instead of the host umpire one
2020-07-24 09:13:21 -07:00 · 2020-07-22 15:51:04 -07:00 · 2020-06-23 12:11:33 -07:00 · 2020-06-22 14:42:48 -07:00 · 2020-06-22 14:29:28 -07:00 · 2020-06-22 12:47:31 -07:00
33 changed files with 1673 additions and 290 deletions
@@ -125,6 +125,7 @@ MFEM_USE_GINKGO        = NO
 MFEM_USE_GNUTLS        = NO
 MFEM_USE_NETCDF        = NO
 MFEM_USE_PETSC         = NO
+MFEM_USE_SLEPC         = NO
 MFEM_USE_MPFR          = NO
 MFEM_USE_SIDRE         = NO
 MFEM_USE_CONDUIT       = NO
@@ -276,6 +277,20 @@ ifeq ($(PETSC_FOUND),YES)
      -L$(abspath $(PETSC_DIR))/lib -lpetsc $(PETSC_LIB)
 endif

+SLEPC_DIR := $(MFEM_DIR)/../slepc
+SLEPC_VARS := $(SLEPC_DIR)/lib/slepc/conf/slepc_variables
+SLEPC_FOUND := $(if $(wildcard $(SLEPC_VARS)),YES,)
+SLEPC_INC_VAR = SLEPC_INCLUDE
+SLEPC_LIB_VAR = SLEPC_EXTERNAL_LIB
+ifeq ($(SLEPC_FOUND),YES)
+   SLEPC_OPT := $(shell sed -n "s/$(SLEPC_INC_VAR) *= *//p" $(SLEPC_VARS))
+   # Some additional external libraries might be defined in this file
+   -include ${SLEPC_DIR}/${PETSC_ARCH}/lib/slepc/conf/slepcvariables
+   SLEPC_LIB := $(shell sed -n "s/$(SLEPC_LIB_VAR) *= *//p" $(SLEPC_VARS))
+   SLEPC_LIB := -Wl,-rpath,$(abspath $(SLEPC_DIR))/$(PETSC_ARCH)/lib\
+      -L$(abspath $(SLEPC_DIR))/$(PETSC_ARCH)/lib -lslepc $(SLEPC_LIB)
+endif
+
 # MPFR library configuration
 MPFR_OPT =
 MPFR_LIB = -lmpfr
@@ -324,9 +339,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 =
@@ -16,6 +16,7 @@
 //               mpirun -np 4 ex4p -m ../data/periodic-square.mesh -no-bc
 //               mpirun -np 4 ex4p -m ../data/periodic-cube.mesh -no-bc
 //               mpirun -np 4 ex4p -m ../data/amr-quad.mesh
+//               mpirun -np 3 ex4p -m ../data/amr-quad.mesh -o 2 -hb
 //               mpirun -np 4 ex4p -m ../data/amr-hex.mesh -o 2 -sc
 //               mpirun -np 4 ex4p -m ../data/amr-hex.mesh -o 2 -hb
 //               mpirun -np 4 ex4p -m ../data/star-surf.mesh -o 3 -hb
@@ -16,6 +16,7 @@
 //    mpirun -np 4 ex9p -m ../data/disc-nurbs.mesh -p 2 -rp 1 -dt 0.005 -tf 9
 //    mpirun -np 4 ex9p -m ../data/periodic-square.mesh -p 3 -rp 2 -dt 0.0025 -tf 9 -vs 20
 //    mpirun -np 4 ex9p -m ../data/periodic-cube.mesh -p 0 -o 2 -rp 1 -dt 0.01 -tf 8
+//    mpirun -np 3 ex9p -m ../data/amr-hex.mesh -p 1 -rs 1 -rp 0 -dt 0.005 -tf 0.5
 //
 // Device sample runs:
 //    mpirun -np 4 ex9p -pa
@@ -176,42 +176,41 @@ static void PADiffusionSetup3D(const int Q1D,
   auto J = Reshape(j.Read(), NQ, 3, 3, NE);
   auto C = const_c ? Reshape(c.Read(), 1, 1) : Reshape(c.Read(), NQ, NE);
   auto D = Reshape(d.Write(), NQ, 6, NE);
-   MFEM_FORALL(e, NE,
+   MFEM_FORALL(eq, NE*NQ,
   {
-      for (int q = 0; q < NQ; ++q)
-      {
-         const double J11 = J(q,0,0,e);
-         const double J21 = J(q,1,0,e);
-         const double J31 = J(q,2,0,e);
-         const double J12 = J(q,0,1,e);
-         const double J22 = J(q,1,1,e);
-         const double J32 = J(q,2,1,e);
-         const double J13 = J(q,0,2,e);
-         const double J23 = J(q,1,2,e);
-         const double J33 = J(q,2,2,e);
-         const double detJ = J11 * (J22 * J33 - J32 * J23) -
-         /* */               J21 * (J12 * J33 - J32 * J13) +
-         /* */               J31 * (J12 * J23 - J22 * J13);
-         const double coeff = const_c ? C(0,0) : C(q,e);
-         const double c_detJ = W[q] * coeff / detJ;
-         // adj(J)
-         const double A11 = (J22 * J33) - (J23 * J32);
-         const double A12 = (J32 * J13) - (J12 * J33);
-         const double A13 = (J12 * J23) - (J22 * J13);
-         const double A21 = (J31 * J23) - (J21 * J33);
-         const double A22 = (J11 * J33) - (J13 * J31);
-         const double A23 = (J21 * J13) - (J11 * J23);
-         const double A31 = (J21 * J32) - (J31 * J22);
-         const double A32 = (J31 * J12) - (J11 * J32);
-         const double A33 = (J11 * J22) - (J12 * J21);
-         // detJ J^{-1} J^{-T} = (1/detJ) adj(J) adj(J)^T
-         D(q,0,e) = c_detJ * (A11*A11 + A12*A12 + A13*A13); // 1,1
-         D(q,1,e) = c_detJ * (A11*A21 + A12*A22 + A13*A23); // 2,1
-         D(q,2,e) = c_detJ * (A11*A31 + A12*A32 + A13*A33); // 3,1
-         D(q,3,e) = c_detJ * (A21*A21 + A22*A22 + A23*A23); // 2,2
-         D(q,4,e) = c_detJ * (A21*A31 + A22*A32 + A23*A33); // 3,2
-         D(q,5,e) = c_detJ * (A31*A31 + A32*A32 + A33*A33); // 3,3
-      }
+      const int e = eq / NQ;
+      const int q = eq % NQ;
+      const double J11 = J(q,0,0,e);
+      const double J21 = J(q,1,0,e);
+      const double J31 = J(q,2,0,e);
+      const double J12 = J(q,0,1,e);
+      const double J22 = J(q,1,1,e);
+      const double J32 = J(q,2,1,e);
+      const double J13 = J(q,0,2,e);
+      const double J23 = J(q,1,2,e);
+      const double J33 = J(q,2,2,e);
+      const double detJ = J11 * (J22 * J33 - J32 * J23) -
+      /* */               J21 * (J12 * J33 - J32 * J13) +
+      /* */               J31 * (J12 * J23 - J22 * J13);
+      const double coeff = const_c ? C(0,0) : C(q,e);
+      const double c_detJ = W[q] * coeff / detJ;
+      // adj(J)
+      const double A11 = (J22 * J33) - (J23 * J32);
+      const double A12 = (J32 * J13) - (J12 * J33);
+      const double A13 = (J12 * J23) - (J22 * J13);
+      const double A21 = (J31 * J23) - (J21 * J33);
+      const double A22 = (J11 * J33) - (J13 * J31);
+      const double A23 = (J21 * J13) - (J11 * J23);
+      const double A31 = (J21 * J32) - (J31 * J22);
+      const double A32 = (J31 * J12) - (J11 * J32);
+      const double A33 = (J11 * J22) - (J12 * J21);
+      // detJ J^{-1} J^{-T} = (1/detJ) adj(J) adj(J)^T
+      D(q,0,e) = c_detJ * (A11*A11 + A12*A12 + A13*A13); // 1,1
+      D(q,1,e) = c_detJ * (A11*A21 + A12*A22 + A13*A23); // 2,1
+      D(q,2,e) = c_detJ * (A11*A31 + A12*A32 + A13*A33); // 3,1
+      D(q,3,e) = c_detJ * (A21*A21 + A22*A22 + A23*A23); // 2,2
+      D(q,4,e) = c_detJ * (A21*A31 + A22*A32 + A23*A33); // 3,2
+      D(q,5,e) = c_detJ * (A31*A31 + A32*A32 + A33*A33); // 3,3
   });
 }

@@ -25,6 +25,7 @@ namespace mfem

 void MassIntegrator::SetupPA(const FiniteElementSpace &fes, const bool force)
 {
+
   // Assuming the same element type
   fespace = &fes;
   Mesh *mesh = fes.GetMesh();
@@ -51,21 +52,30 @@ void MassIntegrator::SetupPA(const FiniteElementSpace &fes, const bool force)
   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) = 1.0;
+   }
+   else if (QuadratureCoefficient* cQ = dynamic_cast<QuadratureCoefficient*>(Q))
+   {
+      coeff = cQ->Data();
+      own_coeff = false;
   }
   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);
@@ -80,11 +90,11 @@ void MassIntegrator::SetupPA(const FiniteElementSpace &fes, const bool force)
   {
      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,
      {
@@ -104,11 +114,11 @@ void MassIntegrator::SetupPA(const FiniteElementSpace &fes, const bool force)
   {
      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,3,3,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,
      {
@@ -125,6 +135,8 @@ void MassIntegrator::SetupPA(const FiniteElementSpace &fes, const bool force)
         }
      });
   }
+
+   if (own_coeff) { delete coeff; }
 }

 void MassIntegrator::AssemblePA(const FiniteElementSpace &fes)
@@ -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)
 {
@@ -416,13 +424,43 @@ double DeterminantCoefficient::Eval(ElementTransformation &T,
   return ma.Det();
 }

-VectorSumCoefficient::VectorSumCoefficient(VectorCoefficient &A,
-                                           VectorCoefficient &B,
-                                           double _alpha, double _beta)
-   : VectorCoefficient(A.GetVDim()), a(&A), b(&B), alpha(_alpha), beta(_beta),
-     va(A.GetVDim())
+VectorSumCoefficient::VectorSumCoefficient(int dim)
+   : VectorCoefficient(dim),
+     ACoef(NULL), BCoef(NULL),
+     A(dim), B(dim),
+     alphaCoef(NULL), betaCoef(NULL),
+     alpha(1.0), beta(1.0)
 {
-   MFEM_ASSERT(A.GetVDim() == B.GetVDim(),
+   A = 0.0; B = 0.0;
+}
+
+VectorSumCoefficient::VectorSumCoefficient(VectorCoefficient &_A,
+                                           VectorCoefficient &_B,
+                                           double _alpha, double _beta)
+   : VectorCoefficient(_A.GetVDim()),
+     ACoef(&_A), BCoef(&_B),
+     A(_A.GetVDim()), B(_A.GetVDim()),
+     alphaCoef(NULL), betaCoef(NULL),
+     alpha(_alpha), beta(_beta)
+{
+   MFEM_ASSERT(_A.GetVDim() == _B.GetVDim(),
+               "VectorSumCoefficient:  "
+               "Arguments must have the same dimension.");
+}
+
+VectorSumCoefficient::VectorSumCoefficient(VectorCoefficient &_A,
+                                           VectorCoefficient &_B,
+                                           Coefficient &_alpha,
+                                           Coefficient &_beta)
+   : VectorCoefficient(_A.GetVDim()),
+     ACoef(&_A), BCoef(&_B),
+     A(_A.GetVDim()),
+     B(_A.GetVDim()),
+     alphaCoef(&_alpha),
+     betaCoef(&_beta),
+     alpha(0.0), beta(0.0)
+{
+   MFEM_ASSERT(_A.GetVDim() == _B.GetVDim(),
               "VectorSumCoefficient:  "
               "Arguments must have the same dimension.");
 }
@@ -430,26 +468,47 @@ VectorSumCoefficient::VectorSumCoefficient(VectorCoefficient &A,
 void VectorSumCoefficient::Eval(Vector &V, ElementTransformation &T,
                                const IntegrationPoint &ip)
 {
-   b->Eval(V, T, ip);
-   if ( beta != 1.0 ) { V *= beta; }
-   a->Eval(va, T, ip);
-   V.Add(alpha, va);
+   V.SetSize(A.Size());
+   if (    ACoef) { ACoef->Eval(A, T, ip); }
+   if (    BCoef) { BCoef->Eval(B, T, ip); }
+   if (alphaCoef) { alpha = alphaCoef->Eval(T, ip); }
+   if ( betaCoef) { beta  = betaCoef->Eval(T, ip); }
+   add(alpha, A, beta, B, V);
 }

+ScalarVectorProductCoefficient::ScalarVectorProductCoefficient(
+   double A,
+   VectorCoefficient &B)
+   : VectorCoefficient(B.GetVDim()), aConst(A), a(NULL), b(&B)
+{}
+
 ScalarVectorProductCoefficient::ScalarVectorProductCoefficient(
   Coefficient &A,
   VectorCoefficient &B)
-   : VectorCoefficient(B.GetVDim()), a(&A), b(&B)
+   : VectorCoefficient(B.GetVDim()), aConst(0.0), a(&A), b(&B)
 {}

 void ScalarVectorProductCoefficient::Eval(Vector &V, ElementTransformation &T,
                                          const IntegrationPoint &ip)
 {
-   double sa = a->Eval(T, ip);
+   double sa = (a == NULL) ? aConst : a->Eval(T, ip);
   b->Eval(V, T, ip);
   V *= sa;
 }

+NormalizedVectorCoefficient::NormalizedVectorCoefficient(VectorCoefficient &A,
+                                                         double _tol)
+   : VectorCoefficient(A.GetVDim()), a(&A), tol(_tol)
+{}
+
+void NormalizedVectorCoefficient::Eval(Vector &V, ElementTransformation &T,
+                                       const IntegrationPoint &ip)
+{
+   a->Eval(V, T, ip);
+   double nv = V.Norml2();
+   V *= (nv > tol) ? (1.0/nv) : 0.0;
+}
+
 VectorCrossProductCoefficient::VectorCrossProductCoefficient(
   VectorCoefficient &A,
   VectorCoefficient &B)
@@ -471,17 +530,18 @@ void VectorCrossProductCoefficient::Eval(Vector &V, ElementTransformation &T,
   V[2] = va[0] * vb[1] - va[1] * vb[0];
 }

-MatVecCoefficient::MatVecCoefficient(MatrixCoefficient &A,
-                                     VectorCoefficient &B)
+MatrixVectorProductCoefficient::MatrixVectorProductCoefficient(
+   MatrixCoefficient &A, VectorCoefficient &B)
   : VectorCoefficient(A.GetHeight()), a(&A), b(&B),
     ma(A.GetHeight(), A.GetWidth()), vb(B.GetVDim())
 {
   MFEM_ASSERT(A.GetWidth() == B.GetVDim(),
-               "MatVecCoefficient:  Arguments have incompatible dimensions.");
+               "MatrixVectorProductCoefficient:  "
+               "Arguments have incompatible dimensions.");
 }

-void MatVecCoefficient::Eval(Vector &V, ElementTransformation &T,
-                             const IntegrationPoint &ip)
+void MatrixVectorProductCoefficient::Eval(Vector &V, ElementTransformation &T,
+                                          const IntegrationPoint &ip)
 {
   a->Eval(ma, T, ip);
   b->Eval(vb, T, ip);
@@ -517,17 +577,23 @@ void MatrixSumCoefficient::Eval(DenseMatrix &M, ElementTransformation &T,
   M.Add(alpha, ma);
 }

+ScalarMatrixProductCoefficient::ScalarMatrixProductCoefficient(
+   double A,
+   MatrixCoefficient &B)
+   : MatrixCoefficient(B.GetHeight(), B.GetWidth()), aConst(A), a(NULL), b(&B)
+{}
+
 ScalarMatrixProductCoefficient::ScalarMatrixProductCoefficient(
   Coefficient &A,
   MatrixCoefficient &B)
-   : MatrixCoefficient(B.GetHeight(), B.GetWidth()), a(&A), b(&B)
+   : MatrixCoefficient(B.GetHeight(), B.GetWidth()), aConst(0.0), a(&A), b(&B)
 {}

 void ScalarMatrixProductCoefficient::Eval(DenseMatrix &M,
                                          ElementTransformation &T,
                                          const IntegrationPoint &ip)
 {
-   double sa = a->Eval(T, ip);
+   double sa = (a == NULL) ? aConst : a->Eval(T, ip);
   b->Eval(M, T, ip);
   M *= sa;
 }
@@ -581,6 +647,30 @@ void OuterProductCoefficient::Eval(DenseMatrix &M, ElementTransformation &T,
   }
 }

+CrossCrossCoefficient::CrossCrossCoefficient(Coefficient &A,
+                                             VectorCoefficient &K)
+   : MatrixCoefficient(K.GetVDim(), K.GetVDim()), aConst(0.0), a(&A), k(&K),
+     vk(K.GetVDim())
+{}
+
+void CrossCrossCoefficient::Eval(DenseMatrix &M, ElementTransformation &T,
+                                 const IntegrationPoint &ip)
+{
+   k->Eval(vk, T, ip);
+   M.SetSize(vk.Size(), vk.Size());
+   M = 0.0;
+   double k2 = vk*vk;
+   for (int i=0; i<vk.Size(); i++)
+   {
+      M(i, i) = k2;
+      for (int j=0; j<vk.Size(); j++)
+      {
+         M(i, j) -= vk[i] * vk[j];
+      }
+   }
+   M *= ((a == NULL ) ? aConst : a->Eval(T, ip) );
+}
+
 double LpNormLoop(double p, Coefficient &coeff, Mesh &mesh,
                  const IntegrationRule *irs[])
 {
@@ -84,6 +84,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
@@ -852,12 +879,14 @@ public:
                     const IntegrationPoint &ip);
 };

-/// Coefficients based on sums and products of other coefficients
-
-/// Scalar coefficient defined as the sum of two scalar coefficients
+/// Coefficients based on sums, products, or other functions of coefficients.
+///@{
+/** Scalar coefficient defined as the linear combination of two scalar
+    coefficients or a scalar and a scalar coefficient */
 class SumCoefficient : public Coefficient
 {
 private:
+   double aConst;
   Coefficient * a;
   Coefficient * b;

@@ -865,33 +894,141 @@ private:
   double beta;

 public:
-   /// Construct with the two coefficients.  Result is _alpha * A + _beta * B.
+   /// Constructor with one coefficient.  Result is _alpha * A + _beta * B
+   SumCoefficient(double A, Coefficient &B,
+                  double _alpha = 1.0, double _beta = 1.0)
+      : aConst(A), a(NULL), b(&B), alpha(_alpha), beta(_beta) { }
+
+   /// Constructor with two coefficients.  Result is _alpha * A + _beta * B.
   SumCoefficient(Coefficient &A, Coefficient &B,
                  double _alpha = 1.0, double _beta = 1.0)
-      : a(&A), b(&B), alpha(_alpha), beta(_beta) { }
+      : aConst(0.0), a(&A), b(&B), alpha(_alpha), beta(_beta) { }
+
+   /// Reset the first term in the linear combination as a constant
+   void SetAConst(double A) { a = NULL; aConst = A; }
+   /// Return the first term in the linear combination
+   double GetAConst() const { return aConst; }
+
+   /// Reset the first term in the linear combination
+   void SetACoef(Coefficient &A) { a = &A; }
+   /// Return the first term in the linear combination
+   Coefficient * GetACoef() const { return a; }
+
+   /// Reset the second term in the linear combination
+   void SetBCoef(Coefficient &B) { b = &B; }
+   /// Return the second term in the linear combination
+   Coefficient * GetBCoef() const { return b; }
+
+   /// Reset the factor in front of the first term in the linear combination
+   void SetAlpha(double _alpha) { alpha = _alpha; }
+   /// Return the factor in front of the first term in the linear combination
+   double GetAlpha() const { return alpha; }
+
+   /// Reset the factor in front of the second term in the linear combination
+   void SetBeta(double _beta) { beta = _beta; }
+   /// Return the factor in front of the second term in the linear combination
+   double GetBeta() const { return beta; }

   /// Evaluate the coefficient at @a ip.
   virtual double Eval(ElementTransformation &T,
                       const IntegrationPoint &ip)
-   { return alpha * a->Eval(T, ip) + beta * b->Eval(T, ip); }
+   {
+      return alpha * ((a == NULL ) ? aConst : a->Eval(T, ip) )
+             + beta * b->Eval(T, ip);
+   }
 };

-/// Scalar coefficient defined as the product of two scalar coefficients
+/** Scalar coefficient defined as the product of two scalar coefficients or
+    a scalar and a scalar coefficient. */
 class ProductCoefficient : public Coefficient
 {
 private:
+   double aConst;
   Coefficient * a;
   Coefficient * b;

 public:
-   /// Construct with the two coefficients.  Result is A * B.
+   /// Constructor with one coefficient.  Result is A * B.
+   ProductCoefficient(double A, Coefficient &B)
+      : aConst(A), a(NULL), b(&B) { }
+
+   /// Constructor with two coefficients.  Result is A * B.
   ProductCoefficient(Coefficient &A, Coefficient &B)
-      : a(&A), b(&B) { }
+      : aConst(0.0), a(&A), b(&B) { }
+
+   /// Reset the first term in the product as a constant
+   void SetAConst(double A) { a = NULL; aConst = A; }
+   /// Return the first term in the product
+   double GetAConst() const { return aConst; }
+
+   /// Reset the first term in the product
+   void SetACoef(Coefficient &A) { a = &A; }
+   /// Return the first term in the product
+   Coefficient * GetACoef() const { return a; }
+
+   /// Reset the second term in the product
+   void SetBCoef(Coefficient &B) { b = &B; }
+   /// Return the second term in the product
+   Coefficient * GetBCoef() const { return b; }

   /// Evaluate the coefficient at @a ip.
   virtual double Eval(ElementTransformation &T,
                       const IntegrationPoint &ip)
-   { return a->Eval(T, ip) * b->Eval(T, ip); }
+   { return ((a == NULL ) ? aConst : a->Eval(T, ip) ) * b->Eval(T, ip); }
+};
+
+/** Scalar coefficient defined as the ratio of two scalars where one or both
+    scalars are scalar coefficients. */
+class RatioCoefficient : public Coefficient
+{
+private:
+   double aConst;
+   double bConst;
+   Coefficient * a;
+   Coefficient * b;
+
+public:
+   /** Initialize a coefficient which returns A / B where @a A is a
+       constant and @a B is a scalar coefficient */
+   RatioCoefficient(double A, Coefficient &B)
+      : aConst(A), bConst(1.0), a(NULL), b(&B) { }
+   /** Initialize a coefficient which returns A / B where @a A and @a B are both
+       scalar coefficients */
+   RatioCoefficient(Coefficient &A, Coefficient &B)
+      : aConst(0.0), bConst(1.0), a(&A), b(&B) { }
+   /** Initialize a coefficient which returns A / B where @a A is a
+       scalar coefficient and @a B is a constant */
+   RatioCoefficient(Coefficient &A, double B)
+      : aConst(0.0), bConst(B), a(&A), b(NULL) { }
+
+   /// Reset the numerator in the ratio as a constant
+   void SetAConst(double A) { a = NULL; aConst = A; }
+   /// Return the numerator of the ratio
+   double GetAConst() const { return aConst; }
+
+   /// Reset the denominator in the ratio as a constant
+   void SetBConst(double B) { b = NULL; bConst = B; }
+   /// Return the denominator of the ratio
+   double GetBConst() const { return bConst; }
+
+   /// Reset the numerator in the ratio
+   void SetACoef(Coefficient &A) { a = &A; }
+   /// Return the numerator of the ratio
+   Coefficient * GetACoef() const { return a; }
+
+   /// Reset the denominator in the ratio
+   void SetBCoef(Coefficient &B) { b = &B; }
+   /// Return the denominator of the ratio
+   Coefficient * GetBCoef() const { return b; }
+
+   /// Evaluate the coefficient
+   virtual double Eval(ElementTransformation &T,
+                       const IntegrationPoint &ip)
+   {
+      double den = (b == NULL ) ? bConst : b->Eval(T, ip);
+      MFEM_ASSERT(den != 0.0, "Division by zero in RatioCoefficient");
+      return ((a == NULL ) ? aConst : a->Eval(T, ip) ) / den;
+   }
 };

 /// Scalar coefficient defined as a scalar raised to a power
@@ -907,6 +1044,16 @@ public:
   PowerCoefficient(Coefficient &A, double _p)
      : a(&A), p(_p) { }

+   /// Reset the base coefficient
+   void SetACoef(Coefficient &A) { a = &A; }
+   /// Return the base coefficient
+   Coefficient * GetACoef() const { return a; }
+
+   /// Reset the exponent
+   void SetExponent(double _p) { p = _p; }
+   /// Return the exponent
+   double GetExponent() const { return p; }
+
   /// Evaluate the coefficient at @a ip.
   virtual double Eval(ElementTransformation &T,
                       const IntegrationPoint &ip)
@@ -927,6 +1074,16 @@ public:
   /// Construct with the two vector coefficients.  Result is \f$ A \cdot B \f$.
   InnerProductCoefficient(VectorCoefficient &A, VectorCoefficient &B);

+   /// Reset the first vector in the inner product
+   void SetACoef(VectorCoefficient &A) { a = &A; }
+   /// Return the first vector coefficient in the inner product
+   VectorCoefficient * GetACoef() const { return a; }
+
+   /// Reset the second vector in the inner product
+   void SetBCoef(VectorCoefficient &B) { b = &B; }
+   /// Return the second vector coefficient in the inner product
+   VectorCoefficient * GetBCoef() const { return b; }
+
   /// Evaluate the coefficient at @a ip.
   virtual double Eval(ElementTransformation &T,
                       const IntegrationPoint &ip);
@@ -943,9 +1100,19 @@ private:
   mutable Vector vb;

 public:
-   /// Construct with the two vector coefficients.  Result is \f$ A_x B_y - A_y * B_x; \f$.
+   /// Constructor with two vector coefficients.  Result is \f$ A_x B_y - A_y * B_x; \f$.
   VectorRotProductCoefficient(VectorCoefficient &A, VectorCoefficient &B);

+   /// Reset the first vector in the product
+   void SetACoef(VectorCoefficient &A) { a = &A; }
+   /// Return the first vector of the product
+   VectorCoefficient * GetACoef() const { return a; }
+
+   /// Reset the second vector in the product
+   void SetBCoef(VectorCoefficient &B) { b = &B; }
+   /// Return the second vector of the product
+   VectorCoefficient * GetBCoef() const { return b; }
+
   /// Evaluate the coefficient at @a ip.
   virtual double Eval(ElementTransformation &T,
                       const IntegrationPoint &ip);
@@ -963,17 +1130,28 @@ public:
   /// Construct with the matrix.
   DeterminantCoefficient(MatrixCoefficient &A);

+   /// Reset the matrix coefficient
+   void SetACoef(MatrixCoefficient &A) { a = &A; }
+   /// Return the matrix coefficient
+   MatrixCoefficient * GetACoef() const { return a; }
+
   /// Evaluate the determinant coefficient at @a ip.
   virtual double Eval(ElementTransformation &T,
                       const IntegrationPoint &ip);
 };

-/// Vector coefficient defined as the sum of two vector coefficients
+/// Vector coefficient defined as the linear combination of two vectors
 class VectorSumCoefficient : public VectorCoefficient
 {
 private:
-   VectorCoefficient * a;
-   VectorCoefficient * b;
+   VectorCoefficient * ACoef;
+   VectorCoefficient * BCoef;
+
+   Vector A;
+   Vector B;
+
+   Coefficient * alphaCoef;
+   Coefficient * betaCoef;

   double alpha;
   double beta;
@@ -981,10 +1159,60 @@ private:
   mutable Vector va;

 public:
-   /// Construct with the two vector coefficients.  Result is _alpha * A + _beta * B.
+   /** Constructor with no coefficients.
+       To be used with the various "Set" methods */
+   VectorSumCoefficient(int dim);
+
+   /** Constructor with two vector coefficients.
+       Result is _alpha * A + _beta * B */
   VectorSumCoefficient(VectorCoefficient &A, VectorCoefficient &B,
                        double _alpha = 1.0, double _beta = 1.0);

+   /** Constructor with scalar coefficients.
+       Result is _alpha * _A + _beta * _B */
+   VectorSumCoefficient(VectorCoefficient &_A, VectorCoefficient &_B,
+                        Coefficient &_alpha, Coefficient &_beta);
+
+   /// Reset the first vector coefficient
+   void SetACoef(VectorCoefficient &A) { ACoef = &A; }
+   /// Return the first vector coefficient
+   VectorCoefficient * GetACoef() const { return ACoef; }
+
+   /// Reset the second vector coefficient
+   void SetBCoef(VectorCoefficient &B) { BCoef = &B; }
+   /// Return the second vector coefficient
+   VectorCoefficient * GetBCoef() const { return BCoef; }
+
+   /// Reset the factor in front of the first vector coefficient
+   void SetAlphaCoef(Coefficient &A) { alphaCoef = &A; }
+   /// Return the factor in front of the first vector coefficient
+   Coefficient * GetAlphaCoef() const { return alphaCoef; }
+
+   /// Reset the factor in front of the second vector coefficient
+   void SetBetaCoef(Coefficient &B) { betaCoef = &B; }
+   /// Return the factor in front of the second vector coefficient
+   Coefficient * GetBetaCoef() const { return betaCoef; }
+
+   /// Reset the first vector as a constant
+   void SetA(const Vector &_A) { A = _A; ACoef = NULL; }
+   /// Return the first vector constant
+   const Vector & GetA() const { return A; }
+
+   /// Reset the second vector as a constant
+   void SetB(const Vector &_B) { B = _B; BCoef = NULL; }
+   /// Return the second vector constant
+   const Vector & GetB() const { return B; }
+
+   /// Reset the factor in front of the first vector coefficient as a constant
+   void SetAlpha(double _alpha) { alpha = _alpha; alphaCoef = NULL; }
+   /// Return the factor in front of the first vector coefficient
+   double GetAlpha() const { return alpha; }
+
+   /// Reset the factor in front of the second vector coefficient as a constant
+   void SetBeta(double _beta) { beta = _beta; betaCoef = NULL; }
+   /// Return the factor in front of the second vector coefficient
+   double GetBeta() const { return beta; }
+
   /// Evaluate the coefficient at @a ip.
   virtual void Eval(Vector &V, ElementTransformation &T,
                     const IntegrationPoint &ip);
@@ -995,13 +1223,60 @@ public:
 class ScalarVectorProductCoefficient : public VectorCoefficient
 {
 private:
+   double aConst;
   Coefficient * a;
   VectorCoefficient * b;

 public:
-   /// Construct with the two coefficients.  Result is A * B.
+   /// Constructor with constant and vector coefficient.  Result is A * B.
+   ScalarVectorProductCoefficient(double A, VectorCoefficient &B);
+
+   /// Constructor with two coefficients.  Result is A * B.
   ScalarVectorProductCoefficient(Coefficient &A, VectorCoefficient &B);

+   /// Reset the scalar factor as a constant
+   void SetAConst(double A) { a = NULL; aConst = A; }
+   /// Return the scalar factor
+   double GetAConst() const { return aConst; }
+
+   /// Reset the scalar factor
+   void SetACoef(Coefficient &A) { a = &A; }
+   /// Return the scalar factor
+   Coefficient * GetACoef() const { return a; }
+
+   /// Reset the vector factor
+   void SetBCoef(VectorCoefficient &B) { b = &B; }
+   /// Return the vector factor
+   VectorCoefficient * GetBCoef() const { return b; }
+
+   /// Evaluate the coefficient at @a ip.
+   virtual void Eval(Vector &V, ElementTransformation &T,
+                     const IntegrationPoint &ip);
+   using VectorCoefficient::Eval;
+};
+
+/// Vector coefficient defined as a normalized vector field (returns v/|v|)
+class NormalizedVectorCoefficient : public VectorCoefficient
+{
+private:
+   VectorCoefficient * a;
+
+   double tol;
+
+public:
+   /** @brief Return a vector normalized to a length of one
+
+       This class evaluates the vector coefficient @a A and, if |A| > @a tol,
+       returns the normalized vector A / |A|.  If |A| <= @a tol, the zero
+       vector is returned.
+   */
+   NormalizedVectorCoefficient(VectorCoefficient &A, double tol = 1e-6);
+
+   /// Reset the vector coefficient
+   void SetACoef(VectorCoefficient &A) { a = &A; }
+   /// Return the vector coefficient
+   VectorCoefficient * GetACoef() const { return a; }
+
   /// Evaluate the coefficient at @a ip.
   virtual void Eval(Vector &V, ElementTransformation &T,
                     const IntegrationPoint &ip);
@@ -1022,6 +1297,16 @@ public:
   /// Construct with the two coefficients.  Result is A x B.
   VectorCrossProductCoefficient(VectorCoefficient &A, VectorCoefficient &B);

+   /// Reset the first term in the product
+   void SetACoef(VectorCoefficient &A) { a = &A; }
+   /// Return the first term in the product
+   VectorCoefficient * GetACoef() const { return a; }
+
+   /// Reset the second term in the product
+   void SetBCoef(VectorCoefficient &B) { b = &B; }
+   /// Return the second term in the product
+   VectorCoefficient * GetBCoef() const { return b; }
+
   /// Evaluate the coefficient at @a ip.
   virtual void Eval(Vector &V, ElementTransformation &T,
                     const IntegrationPoint &ip);
@@ -1030,7 +1315,7 @@ public:

 /** @brief Vector coefficient defined as a product of a matrix coefficient and
    a vector coefficient. */
-class MatVecCoefficient : public VectorCoefficient
+class MatrixVectorProductCoefficient : public VectorCoefficient
 {
 private:
   MatrixCoefficient * a;
@@ -1040,8 +1325,18 @@ private:
   mutable Vector vb;

 public:
-   /// Construct with the two coefficients.  Result is A*B.
-   MatVecCoefficient(MatrixCoefficient &A, VectorCoefficient &B);
+   /// Constructor with two coefficients.  Result is A*B.
+   MatrixVectorProductCoefficient(MatrixCoefficient &A, VectorCoefficient &B);
+
+   /// Reset the matrix coefficient
+   void SetACoef(MatrixCoefficient &A) { a = &A; }
+   /// Return the matrix coefficient
+   MatrixCoefficient * GetACoef() const { return a; }
+
+   /// Reset the vector coefficient
+   void SetBCoef(VectorCoefficient &B) { b = &B; }
+   /// Return the vector coefficient
+   VectorCoefficient * GetBCoef() const { return b; }

   /// Evaluate the vector coefficient at @a ip.
   virtual void Eval(Vector &V, ElementTransformation &T,
@@ -1049,6 +1344,9 @@ public:
   using VectorCoefficient::Eval;
 };

+/// Convenient alias for the MatrixVectorProductCoefficient
+typedef MatrixVectorProductCoefficient MatVecCoefficient;
+
 /// Constant matrix coefficient defined as the identity of dimension d
 class IdentityMatrixCoefficient : public MatrixCoefficient
 {
@@ -1065,7 +1363,7 @@ public:
                     const IntegrationPoint &ip);
 };

-/// Matrix coefficient defined as the sum of two matrix coefficients.
+/// Matrix coefficient defined as the linear combination of two matrices
 class MatrixSumCoefficient : public MatrixCoefficient
 {
 private:
@@ -1082,6 +1380,26 @@ public:
   MatrixSumCoefficient(MatrixCoefficient &A, MatrixCoefficient &B,
                        double _alpha = 1.0, double _beta = 1.0);

+   /// Reset the first matrix coefficient
+   void SetACoef(MatrixCoefficient &A) { a = &A; }
+   /// Return the first matrix coefficient
+   MatrixCoefficient * GetACoef() const { return a; }
+
+   /// Reset the second matrix coefficient
+   void SetBCoef(MatrixCoefficient &B) { b = &B; }
+   /// Return the second matrix coefficient
+   MatrixCoefficient * GetBCoef() const { return b; }
+
+   /// Reset the factor in front of the first matrix coefficient
+   void SetAlpha(double _alpha) { alpha = _alpha; }
+   /// Return the factor in front of the first matrix coefficient
+   double GetAlpha() const { return alpha; }
+
+   /// Reset the factor in front of the second matrix coefficient
+   void SetBeta(double _beta) { beta = _beta; }
+   /// Return the factor in front of the second matrix coefficient
+   double GetBeta() const { return beta; }
+
   /// Evaluate the matrix coefficient at @a ip.
   virtual void Eval(DenseMatrix &M, ElementTransformation &T,
                     const IntegrationPoint &ip);
@@ -1092,13 +1410,32 @@ public:
 class ScalarMatrixProductCoefficient : public MatrixCoefficient
 {
 private:
+   double aConst;
   Coefficient * a;
   MatrixCoefficient * b;

 public:
-   /// Construct with the two coefficients.  Result is A*B.
+   /// Constructor with one coefficient.  Result is A*B.
+   ScalarMatrixProductCoefficient(double A, MatrixCoefficient &B);
+
+   /// Constructor with two coefficients.  Result is A*B.
   ScalarMatrixProductCoefficient(Coefficient &A, MatrixCoefficient &B);

+   /// Reset the scalar factor as a constant
+   void SetAConst(double A) { a = NULL; aConst = A; }
+   /// Return the scalar factor
+   double GetAConst() const { return aConst; }
+
+   /// Reset the scalar factor
+   void SetACoef(Coefficient &A) { a = &A; }
+   /// Return the scalar factor
+   Coefficient * GetACoef() const { return a; }
+
+   /// Reset the matrix factor
+   void SetBCoef(MatrixCoefficient &B) { b = &B; }
+   /// Return the matrix factor
+   MatrixCoefficient * GetBCoef() const { return b; }
+
   /// Evaluate the matrix coefficient at @a ip.
   virtual void Eval(DenseMatrix &M, ElementTransformation &T,
                     const IntegrationPoint &ip);
@@ -1114,6 +1451,11 @@ public:
   /// Construct with the matrix coefficient.  Result is \f$ A^T \f$.
   TransposeMatrixCoefficient(MatrixCoefficient &A);

+   /// Reset the matrix coefficient
+   void SetACoef(MatrixCoefficient &A) { a = &A; }
+   /// Return the matrix coefficient
+   MatrixCoefficient * GetACoef() const { return a; }
+
   /// Evaluate the matrix coefficient at @a ip.
   virtual void Eval(DenseMatrix &M, ElementTransformation &T,
                     const IntegrationPoint &ip);
@@ -1129,6 +1471,11 @@ public:
   /// Construct with the matrix coefficient.  Result is \f$ A^{-1} \f$.
   InverseMatrixCoefficient(MatrixCoefficient &A);

+   /// Reset the matrix coefficient
+   void SetACoef(MatrixCoefficient &A) { a = &A; }
+   /// Return the matrix coefficient
+   MatrixCoefficient * GetACoef() const { return a; }
+
   /// Evaluate the matrix coefficient at @a ip.
   virtual void Eval(DenseMatrix &M, ElementTransformation &T,
                     const IntegrationPoint &ip);
@@ -1148,11 +1495,61 @@ public:
   /// Construct with two vector coefficients.  Result is \f$ A B^T \f$.
   OuterProductCoefficient(VectorCoefficient &A, VectorCoefficient &B);

+   /// Reset the first vector in the outer product
+   void SetACoef(VectorCoefficient &A) { a = &A; }
+   /// Return the first vector coefficient in the outer product
+   VectorCoefficient * GetACoef() const { return a; }
+
+   /// Reset the second vector in the outer product
+   void SetBCoef(VectorCoefficient &B) { b = &B; }
+   /// Return the second vector coefficient in the outer product
+   VectorCoefficient * GetBCoef() const { return b; }
+
   /// Evaluate the matrix coefficient at @a ip.
   virtual void Eval(DenseMatrix &M, ElementTransformation &T,
                     const IntegrationPoint &ip);
 };

+/** @brief Matrix coefficient defined as -a k x k x, for a vector k and scalar a
+
+    This coefficient returns \f$a * (|k|^2 I - k \otimes k)\f$, where I is
+    the identity matrix and \f$\otimes\f$ indicates the outer product.  This
+    can be evaluated for vectors of any dimension but in three
+    dimensions it corresponds to computing the cross product with k twice.
+*/
+class CrossCrossCoefficient : public MatrixCoefficient
+{
+private:
+   double aConst;
+   Coefficient * a;
+   VectorCoefficient * k;
+
+   mutable Vector vk;
+
+public:
+   CrossCrossCoefficient(double A, VectorCoefficient &K);
+   CrossCrossCoefficient(Coefficient &A, VectorCoefficient &K);
+
+   /// Reset the scalar factor as a constant
+   void SetAConst(double A) { a = NULL; aConst = A; }
+   /// Return the scalar factor
+   double GetAConst() const { return aConst; }
+
+   /// Reset the scalar factor
+   void SetACoef(Coefficient &A) { a = &A; }
+   /// Return the scalar factor
+   Coefficient * GetACoef() const { return a; }
+
+   /// Reset the vector factor
+   void SetKCoef(VectorCoefficient &K) { k = &K; }
+   /// Return the vector factor
+   VectorCoefficient * GetKCoef() const { return k; }
+
+   /// Evaluate the matrix coefficient at @a ip.
+   virtual void Eval(DenseMatrix &M, ElementTransformation &T,
+                     const IntegrationPoint &ip);
+};
+///@}

 class QuadratureFunction;

@@ -624,4 +624,77 @@ void FaceElementTransformations::Transform(const DenseMatrix &matrix,
   IsoparametricTransformation::Transform(matrix, result);
 }

+double FaceElementTransformations::CheckConsistency(int print_level,
+                                                    std::ostream &out)
+{
+   // Check that the face vertices are mapped to the same physical location
+   // when using the following three transformations:
+   // - the face transformation, *this
+   // - Loc1 + Elem1
+   // - Loc2 + Elem2, if present.
+
+   const bool have_face = (mask & 16);
+   const bool have_el1 = (mask & 1) && (mask & 4);
+   const bool have_el2 = (mask & 2) && (mask & 8) && (Elem2No >= 0);
+   if (int(have_face) + int(have_el1) + int(have_el2) < 2)
+   {
+      // need at least two different transformations to perform a check
+      return 0.0;
+   }
+
+   const IntegrationRule &v_ir = *Geometries.GetVertices(GetGeometryType());
+
+   double max_dist = 0.0;
+   Vector dist(v_ir.GetNPoints());
+   DenseMatrix coords_base, coords_el;
+   IntegrationRule v_eir(v_ir.GetNPoints());
+   if (have_face)
+   {
+      Transform(v_ir, coords_base);
+      if (print_level > 0)
+      {
+         out << "\nface vertex coordinates (from face transform):\n"
+             << "----------------------------------------------\n";
+         coords_base.PrintT(out, coords_base.Height());
+      }
+   }
+   if (have_el1)
+   {
+      Loc1.Transform(v_ir, v_eir);
+      Elem1->Transform(v_eir, coords_el);
+      if (print_level > 0)
+      {
+         out << "\nface vertex coordinates (from element 1 transform):\n"
+             << "---------------------------------------------------\n";
+         coords_el.PrintT(out, coords_el.Height());
+      }
+      if (have_face)
+      {
+         coords_el -= coords_base;
+         coords_el.Norm2(dist);
+         max_dist = std::max(max_dist, dist.Normlinf());
+      }
+      else
+      {
+         coords_base = coords_el;
+      }
+   }
+   if (have_el2)
+   {
+      Loc2.Transform(v_ir, v_eir);
+      Elem2->Transform(v_eir, coords_el);
+      if (print_level > 0)
+      {
+         out << "\nface vertex coordinates (from element 2 transform):\n"
+             << "---------------------------------------------------\n";
+         coords_el.PrintT(out, coords_el.Height());
+      }
+      coords_el -= coords_base;
+      coords_el.Norm2(dist);
+      max_dist = std::max(max_dist, dist.Normlinf());
+   }
+
+   return max_dist;
+}
+
 }
@@ -492,6 +492,26 @@ public:
   ElementTransformation & GetElement2Transformation();
   IntegrationPointTransformation & GetIntPoint1Transformation();
   IntegrationPointTransformation & GetIntPoint2Transformation();
+
+   /** @brief Check for self-consistency: compares the result of mapping the
+       reference face vertices to physical coordinates using the three
+       transformations: face, element 1, and element 2.
+
+       @param[in] print_level  If set to a positive number, print the physical
+                               coordinates of the face vertices computed through
+                               all available transformations: face, element 1,
+                               and/or element 2.
+       @param[in,out] out      The output stream to use for printing.
+
+       @returns A maximal distance between physical coordinates of face vertices
+                that should coincide. A successful check should return a small
+                number relative to the mesh extents. If less than 2 of the three
+                transformations are set, returns 0.
+
+       @warning This check will generally fail on periodic boundary faces.
+   */
+   double CheckConsistency(int print_level = 0,
+                           std::ostream &out = mfem::out);
 };

 /**                Elem1(Loc1(x)) = Face(x) = Elem2(Loc2(x))
@@ -1031,13 +1031,13 @@ int GridFunction::GetFaceVectorValues(
   }
   if (di == 0)
   {
-      Transf = fes->GetMesh()->GetFaceElementTransformations(i, 4);
+      Transf = fes->GetMesh()->GetFaceElementTransformations(i, 5);
      Transf->Loc1.Transform(ir, eir);
      GetVectorValues(*Transf->Elem1, eir, vals, &tr);
   }
   else
   {
-      Transf = fes->GetMesh()->GetFaceElementTransformations(i, 8);
+      Transf = fes->GetMesh()->GetFaceElementTransformations(i, 10);
      Transf->Loc2.Transform(ir, eir);
      GetVectorValues(*Transf->Elem2, eir, vals, &tr);
   }
@@ -63,6 +63,53 @@ void DomainLFIntegrator::AssembleDeltaElementVect(
   elvect *= delta->EvalDelta(Trans, Trans.GetIntPoint());
 }

+void DomainLFGradIntegrator::AssembleRHSElementVect(
+   const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
+{
+   int dof = el.GetDof();
+   int spaceDim = Tr.GetSpaceDim();
+
+   dshape.SetSize(dof, spaceDim);
+
+   elvect.SetSize(dof);
+   elvect = 0.0;
+
+   const IntegrationRule *ir = IntRule;
+   if (ir == NULL)
+   {
+      int intorder = 2 * el.GetOrder();
+      ir = &IntRules.Get(el.GetGeomType(), intorder);
+   }
+
+   for (int i = 0; i < ir->GetNPoints(); i++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(i);
+
+      Tr.SetIntPoint(&ip);
+      el.CalcPhysDShape(Tr, dshape);
+
+      Q.Eval(Qvec, Tr, ip);
+      Qvec *= ip.weight * Tr.Weight();
+
+      dshape.AddMult(Qvec, elvect);
+   }
+}
+
+void DomainLFGradIntegrator::AssembleDeltaElementVect(
+   const FiniteElement &fe, ElementTransformation &Trans, Vector &elvect)
+{
+   MFEM_ASSERT(vec_delta != NULL,"coefficient must be VectorDeltaCoefficient");
+   int dof = fe.GetDof();
+   int spaceDim = Trans.GetSpaceDim();
+
+   dshape.SetSize(dof, spaceDim);
+   fe.CalcPhysDShape(Trans, dshape);
+
+   vec_delta->EvalDelta(Qvec, Trans, Trans.GetIntPoint());
+
+   elvect.SetSize(dof);
+   dshape.Mult(Qvec, elvect);
+}

 void BoundaryLFIntegrator::AssembleRHSElementVect(
   const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
@@ -255,7 +302,6 @@ void VectorDomainLFIntegrator::AssembleDeltaElementVect(
   MultVWt(shape, Qvec, elvec_as_mat);
 }

-
 void VectorBoundaryLFIntegrator::AssembleRHSElementVect(
   const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
 {
@@ -332,7 +378,6 @@ void VectorBoundaryLFIntegrator::AssembleRHSElementVect(
   }
 }

-
 void VectorFEDomainLFIntegrator::AssembleRHSElementVect(
   const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
 {
@@ -362,7 +407,6 @@ void VectorFEDomainLFIntegrator::AssembleRHSElementVect(

      QF.Eval (vec, Tr, ip);
      vec *= ip.weight * Tr.Weight();
-
      vshape.AddMult (vec, elvect);
   }
 }
@@ -383,6 +427,125 @@ void VectorFEDomainLFIntegrator::AssembleDeltaElementVect(
   vshape.Mult(vec, elvect);
 }

+void VectorFEDomainLFCurlIntegrator::AssembleRHSElementVect(
+   const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
+{
+   int dof = el.GetDof();
+   int spaceDim = Tr.GetSpaceDim();
+   int n=(spaceDim == 3)? spaceDim : 1;
+   curlshape.SetSize(dof,n);
+   vec.SetSize(n);
+
+   elvect.SetSize(dof);
+   elvect = 0.0;
+
+   const IntegrationRule *ir = IntRule;
+   if (ir == NULL)
+   {
+      int intorder = 2*el.GetOrder();
+      ir = &IntRules.Get(el.GetGeomType(), intorder);
+   }
+
+   for (int i = 0; i < ir->GetNPoints(); i++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(i);
+
+      Tr.SetIntPoint (&ip);
+      el.CalcPhysCurlShape(Tr, curlshape);
+
+      switch (spaceDim)
+      {
+         case 3:
+            MFEM_VERIFY(QF, "VectorFunctionCoefficient not provided");
+            QF->Eval(vec, Tr, ip);
+            break;
+         case 2:
+            MFEM_VERIFY(Q, "FunctionCoefficient (Scalar) not provided");
+            vec[0] = Q->Eval(Tr, ip);
+            break;
+         default:
+            break; // This should be unreachable
+      }
+      vec *= ip.weight * Tr.Weight();
+      curlshape.AddMult (vec, elvect);
+   }
+}
+
+void VectorFEDomainLFCurlIntegrator::AssembleDeltaElementVect(
+   const FiniteElement &fe, ElementTransformation &Trans, Vector &elvect)
+{
+   int spaceDim = Trans.GetSpaceDim();
+   switch (spaceDim)
+   {
+      case 3:
+         MFEM_ASSERT(vec_delta != NULL,
+                     "coefficient must be VectorDeltaCoefficient");
+         break;
+      case 2:
+         MFEM_ASSERT(delta != NULL,
+                     "coefficient must be DeltaCoefficient");
+         break;
+      default:
+         break; // This should be unreachable
+   }
+   int dof = fe.GetDof();
+   int n=(spaceDim == 3)? spaceDim : 1;
+   curlshape.SetSize(dof, n);
+   elvect.SetSize(dof);
+   fe.CalcPhysCurlShape(Trans, curlshape);
+
+   switch (spaceDim)
+   {
+      case 3:
+         vec_delta->EvalDelta(vec, Trans, Trans.GetIntPoint());
+         curlshape.Mult(vec, elvect);
+         break;
+      case 2:
+         curlshape.GetColumn(0,elvect);
+         elvect *= delta->EvalDelta(Trans, Trans.GetIntPoint());
+         break;
+      default:
+         break; // This should be unreachable
+   }
+}
+
+void VectorFEDomainLFDivIntegrator::AssembleRHSElementVect(
+   const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
+{
+   int dof = el.GetDof();
+
+   divshape.SetSize(dof);       // vector of size dof
+   elvect.SetSize(dof);
+   elvect = 0.0;
+
+   const IntegrationRule *ir = IntRule;
+   if (ir == NULL)
+   {
+      int intorder = 2 * el.GetOrder();
+      ir = &IntRules.Get(el.GetGeomType(), intorder);
+   }
+
+   for (int i = 0; i < ir->GetNPoints(); i++)
+   {
+      const IntegrationPoint &ip = ir->IntPoint(i);
+
+      Tr.SetIntPoint (&ip);
+      double val = Tr.Weight() * Q.Eval(Tr, ip);
+      el.CalcPhysDivShape(Tr, divshape);
+
+      add(elvect, ip.weight * val, divshape, elvect);
+   }
+}
+
+void VectorFEDomainLFDivIntegrator::AssembleDeltaElementVect(
+   const FiniteElement &fe, ElementTransformation &Trans, Vector &elvect)
+{
+   MFEM_ASSERT(delta != NULL, "coefficient must be DeltaCoefficient");
+   elvect.SetSize(fe.GetDof());
+   fe.CalcPhysDivShape(Trans, elvect);
+   elvect *= delta->EvalDelta(Trans, Trans.GetIntPoint());
+}
+
 void VectorBoundaryFluxLFIntegrator::AssembleRHSElementVect(
   const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
 {
@@ -448,7 +611,6 @@ void VectorFEBoundaryFluxLFIntegrator::AssembleRHSElementVect(
   }
 }

-
 void VectorFEBoundaryTangentLFIntegrator::AssembleRHSElementVect(
   const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
 {
@@ -483,7 +645,6 @@ void VectorFEBoundaryTangentLFIntegrator::AssembleRHSElementVect(
   }
 }

-
 void BoundaryFlowIntegrator::AssembleRHSElementVect(
   const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
 {
@@ -548,7 +709,6 @@ void BoundaryFlowIntegrator::AssembleRHSElementVect(
   }
 }

-
 void DGDirichletLFIntegrator::AssembleRHSElementVect(
   const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
 {
@@ -637,7 +797,6 @@ void DGDirichletLFIntegrator::AssembleRHSElementVect(
   }
 }

-
 void DGElasticityDirichletLFIntegrator::AssembleRHSElementVect(
   const FiniteElement &el, ElementTransformation &Tr, Vector &elvect)
 {
@@ -119,6 +119,33 @@ public:
   using LinearFormIntegrator::AssembleRHSElementVect;
 };

+/// Class for domain integrator L(v) := (f, grad v)
+class DomainLFGradIntegrator : public DeltaLFIntegrator
+{
+private:
+   Vector shape, Qvec;
+   VectorCoefficient &Q;
+   DenseMatrix dshape;
+
+public:
+   /// Constructs the domain integrator (Q, grad v)
+   DomainLFGradIntegrator(VectorCoefficient &QF)
+      : DeltaLFIntegrator(QF), Q(QF) { }
+
+   /** Given a particular Finite Element and a transformation (Tr)
+       computes the element right hand side element vector, elvect. */
+   virtual void AssembleRHSElementVect(const FiniteElement &el,
+                                       ElementTransformation &Tr,
+                                       Vector &elvect);
+
+   virtual void AssembleDeltaElementVect(const FiniteElement &fe,
+                                         ElementTransformation &Trans,
+                                         Vector &elvect);
+
+   using LinearFormIntegrator::AssembleRHSElementVect;
+};
+
+
 /// Class for boundary integration L(v) := (g, v)
 class BoundaryLFIntegrator : public LinearFormIntegrator
 {
@@ -252,6 +279,56 @@ public:
   using LinearFormIntegrator::AssembleRHSElementVect;
 };

+/// \f$ (Q, curl v)_{\Omega} \f$ for Nedelec Elements)
+class VectorFEDomainLFCurlIntegrator : public DeltaLFIntegrator
+{
+private:
+   VectorCoefficient *QF=nullptr;
+   Coefficient *Q=nullptr;
+   DenseMatrix curlshape;
+   Vector vec;
+
+public:
+   /// Constructs the domain integrator (Q, curl v)
+   VectorFEDomainLFCurlIntegrator(VectorCoefficient &F)
+      : DeltaLFIntegrator(F), QF(&F) { }
+   VectorFEDomainLFCurlIntegrator(Coefficient &F)
+      : DeltaLFIntegrator(F), Q(&F) { }
+
+   virtual void AssembleRHSElementVect(const FiniteElement &el,
+                                       ElementTransformation &Tr,
+                                       Vector &elvect);
+
+   virtual void AssembleDeltaElementVect(const FiniteElement &fe,
+                                         ElementTransformation &Trans,
+                                         Vector &elvect);
+
+   using LinearFormIntegrator::AssembleRHSElementVect;
+};
+
+/// \f$ (Q, div v)_{\Omega} \f$ for RT Elements)
+class VectorFEDomainLFDivIntegrator : public DeltaLFIntegrator
+{
+private:
+   Vector divshape;
+   Coefficient &Q;
+public:
+   /// Constructs the domain integrator (Q, div v)
+   VectorFEDomainLFDivIntegrator(Coefficient &QF)
+      : DeltaLFIntegrator(QF), Q(QF) { }
+
+   /** Given a particular Finite Element and a transformation (Tr)
+       computes the element right hand side element vector, elvect. */
+   virtual void AssembleRHSElementVect(const FiniteElement &el,
+                                       ElementTransformation &Tr,
+                                       Vector &elvect);
+
+   virtual void AssembleDeltaElementVect(const FiniteElement &fe,
+                                         ElementTransformation &Trans,
+                                         Vector &elvect);
+
+   using LinearFormIntegrator::AssembleRHSElementVect;
+};

 /** \f$ (f, v \cdot n)_{\partial\Omega} \f$ for vector test function
    v=(v1,...,vn) where all vi are in the same scalar FE space and f is a
@@ -283,7 +360,7 @@ class VectorFEBoundaryFluxLFIntegrator : public LinearFormIntegrator
 private:
   Coefficient *F;
   Vector shape;
-   int oa, ob; // these contol the quadrature order, see DomainLFIntegrator
+   int oa, ob; // these control the quadrature order, see DomainLFIntegrator

 public:
   VectorFEBoundaryFluxLFIntegrator(int a = 1, int b = -1)
@@ -1172,7 +1172,7 @@ void ParFiniteElementSpace::GetFaceNbrFaceVDofs(int i, Array<int> &vdofs) const
 {
   // Works for NC mesh where 'i' is an index returned by
   // ParMesh::GetSharedFace() such that i >= Mesh::GetNumFaces(), i.e. 'i' is
-   // the index of a ghost.
+   // the index of a ghost face.
   MFEM_ASSERT(Nonconforming() && i >= pmesh->GetNumFaces(), "");
   int el1, el2, inf1, inf2;
   pmesh->GetFaceElements(i, &el1, &el2);
@@ -1212,11 +1212,14 @@ const FiniteElement *ParFiniteElementSpace::GetFaceNbrFE(int i) const

 const FiniteElement *ParFiniteElementSpace::GetFaceNbrFaceFE(int i) const
 {
+   // Works for NC mesh where 'i' is an index returned by
+   // ParMesh::GetSharedFace() such that i >= Mesh::GetNumFaces(), i.e. 'i' is
+   // the index of a ghost face.
   // Works in tandem with GetFaceNbrFaceVDofs() defined above.
+
   MFEM_ASSERT(Nonconforming() && !NURBSext, "");
-   Geometry::Type geom = (pmesh->Dimension() == 2) ?
-                         Geometry::SEGMENT : Geometry::SQUARE;
-   return fec->FiniteElementForGeometry(geom);
+   Geometry::Type face_geom = pmesh->GetFaceGeometryType(i);
+   return fec->FiniteElementForGeometry(face_geom);
 }

 void ParFiniteElementSpace::Lose_Dof_TrueDof_Matrix()
@@ -214,7 +214,7 @@ void ParGridFunction::ExchangeFaceNbrData()
   ParMesh *pmesh = pfes->GetParMesh();

   face_nbr_data.SetSize(pfes->GetFaceNbrVSize());
-   Vector send_data(pfes->send_face_nbr_ldof.Size_of_connections());
+   send_data.SetSize(pfes->send_face_nbr_ldof.Size_of_connections());

   int *send_offset = pfes->send_face_nbr_ldof.GetI();
   const int *d_send_ldof = mfem::Read(pfes->send_face_nbr_ldof.GetJMemory(),
@@ -38,6 +38,11 @@ protected:
       initialized by ExchangeFaceNbrData(). */
   Vector face_nbr_data;

+   /** @brief Vector used as an MPI buffer to send face-neighbor data
+        in ExchangeFaceNbrData() to neighboring processors. */
+   //TODO: Use temporary memory to avoid CUDA malloc allocation cost.
+   Vector send_data;
+
   void ProjectBdrCoefficient(Coefficient *coeff[], VectorCoefficient *vcoeff,
                              Array<int> &attr);

@@ -150,6 +150,7 @@ void QuadratureInterpolator::Eval3D(
   const int nq = maps.nqpt;
   const int ND = T_ND ? T_ND : nd;
   const int NQ = T_NQ ? T_NQ : nq;
+   const int NMAX = NQ > ND ? NQ : ND;
   const int VDIM = T_VDIM ? T_VDIM : vdim;
   MFEM_VERIFY(ND <= MAX_ND3D, "");
   MFEM_VERIFY(NQ <= MAX_NQ3D, "");
@@ -160,22 +161,24 @@ void QuadratureInterpolator::Eval3D(
   auto val = Reshape(q_val.Write(), NQ, VDIM, NE);
   auto der = Reshape(q_der.Write(), NQ, VDIM, 3, NE);
   auto det = Reshape(q_det.Write(), NQ, NE);
-   MFEM_FORALL(e, NE,
+   MFEM_FORALL_2D(e, NE, NMAX, 1, 1,
   {
      const int ND = T_ND ? T_ND : nd;
      const int NQ = T_NQ ? T_NQ : nq;
      const int VDIM = T_VDIM ? T_VDIM : vdim;
      constexpr int max_ND = T_ND ? T_ND : MAX_ND3D;
      constexpr int max_VDIM = T_VDIM ? T_VDIM : MAX_VDIM3D;
-      double s_E[max_VDIM*max_ND];
-      for (int d = 0; d < ND; d++)
+      MFEM_SHARED double s_E[max_VDIM*max_ND];
+      MFEM_FOREACH_THREAD(d, x, ND)
      {
         for (int c = 0; c < VDIM; c++)
         {
            s_E[c+d*VDIM] = E(d,c,e);
         }
      }
-      for (int q = 0; q < NQ; ++q)
+      MFEM_SYNC_THREAD;
+
+      MFEM_FOREACH_THREAD(q, x, NQ)
      {
         if (eval_flags & VALUES)
         {
@@ -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);
@@ -1058,6 +1058,7 @@ void DiscreteAdaptTC::SetDiscreteTargetBase(const GridFunction &tspec_)
   // make a copy of tspec->tspec_temp, increase its size, and
   // copy data from tspec_temp -> tspec, then add new entries
   Vector tspec_temp = tspec;
+   tspec.UseDevice(true);
   tspec.SetSize(ncomp*dof_cnt);

   for (int i = 0; i < tspec_temp.Size(); i++)
@@ -1214,6 +1215,7 @@ void DiscreteAdaptTC::ComputeElementTargets(int e_id, const FiniteElement &fe,
         Array<int> dofs;
         DenseMatrix D_rho(dim), Q_phi(dim), R_theta(dim);
         tspec_fesv->GetElementVDofs(e_id, dofs);
+         tspec.UseDevice(true);
         tspec.GetSubVector(dofs, tspec_vals);

         for (int i = 0; i < ir.GetNPoints(); i++)
@@ -33,10 +33,11 @@ void AdvectorCG::ComputeAtNewPosition(const Vector &new_nodes,
   const int pnt_cnt = new_field.Size()/ncomp;

   new_field = field0;
-
+   new_field.HostReadWrite();
+   Vector new_field_temp;
   for (int i = 0; i < ncomp; i++)
   {
-      Vector new_field_temp(new_field.GetData()+i*pnt_cnt, pnt_cnt);
+      new_field_temp.MakeRef(new_field, i*pnt_cnt, pnt_cnt);
      ComputeAtNewPositionScalar(new_nodes, new_field_temp);
   }

@@ -94,6 +95,7 @@ void AdvectorCG::ComputeAtNewPositionScalar(const Vector &new_nodes,
   double v_max = 0.0;
   const int s = new_field.Size();

+   u.HostReadWrite();
   for (int i = 0; i < s; i++)
   {
      double vel = 0.;
@@ -149,6 +151,7 @@ void AdvectorCG::ComputeAtNewPositionScalar(const Vector &new_nodes,
 #endif

   // Trim the overshoots and undershoots.
+   new_field.HostReadWrite();
   for (int i = 0; i < s; i++)
   {
      if (new_field(i) < glob_minv) { new_field(i) = glob_minv; }
@@ -61,15 +61,11 @@ Device Device::device_singleton;
 bool Device::device_env = false;
 bool Device::mem_host_env = false;
 bool Device::mem_device_env = false;
+#ifdef MFEM_USE_UMPIRE
+bool Device::use_umpire = true;
+#endif

-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)
   {
@@ -135,7 +131,7 @@ Device::Device() : mode(Device::SEQUENTIAL),
      {
         MFEM_ABORT("Unknown memory backend!");
      }
-      mm.Configure(host_mem_type, device_mem_type);
+      mm.Configure(host_mem_type, device_mem_type, device_mem_type);
   }

   if (getenv("MFEM_DEVICE"))
@@ -165,6 +161,8 @@ Device::~Device()
   Get().host_mem_class = MemoryClass::HOST;
   Get().device_mem_type = MemoryType::HOST;
   Get().device_mem_class = MemoryClass::HOST;
+   Get().device_temp_mem_type = MemoryType::HOST;
+   Get().device_temp_mem_class = MemoryClass::HOST;
 }

 void Device::Configure(const std::string &device, const int dev)
@@ -260,6 +258,10 @@ void Device::Print(std::ostream &out)
   if (Device::Allows(Backend::DEVICE_MASK))
   {
      out << ',' << MemoryTypeName[static_cast<int>(device_mem_type)];
+      if (device_temp_mem_type != device_mem_type)
+      {
+         out << ',' << MemoryTypeName[static_cast<int>(device_temp_mem_type)];
+      }
   }
   out << std::endl;
 }
@@ -272,7 +274,8 @@ 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; }
+   // TODO TMS: temporary
+   //if (!mem_host_env && use_umpire) { host_mem_type = MemoryType::HOST_UMPIRE; }
 #endif

   // Enable the device memory type
@@ -296,11 +299,16 @@ void Device::UpdateMemoryTypeAndClass()
         }
         else
         {
-#ifndef MFEM_USE_UMPIRE
-            device_mem_type = MemoryType::DEVICE;
-#else
-            device_mem_type = MemoryType::DEVICE_UMPIRE;
+#ifdef MFEM_USE_UMPIRE
+            if (use_umpire)
+            {
+               device_mem_type = MemoryType::DEVICE_UMPIRE;
+            }
+            else
 #endif
+            {
+               device_mem_type = MemoryType::DEVICE;
+            }
         }
      }
      device_mem_class = MemoryClass::DEVICE;
@@ -320,8 +328,20 @@ void Device::UpdateMemoryTypeAndClass()
      device_mem_type = MemoryType::DEVICE_DEBUG;
   }

+   // Setup device_temp_mem_{type,class}
+   switch (device_mem_type)
+   {
+      case MemoryType::DEVICE_UMPIRE:
+         device_temp_mem_type = MemoryType::DEVICE_TEMP_UMPIRE;
+         break;
+      default:
+         device_temp_mem_type = device_mem_type;
+         break;
+   }
+   device_temp_mem_class = device_mem_class;
+
   // Update the memory manager with the new settings
-   mm.Configure(host_mem_type, device_mem_type);
+   mm.Configure(host_mem_type, device_mem_type, device_temp_mem_type);
 }

 void Device::Enable()
@@ -119,20 +119,26 @@ private:

   static bool device_env, mem_host_env, mem_device_env;
   static Device device_singleton;
+#ifdef MFEM_USE_UMPIRE
+   static bool use_umpire;
+#endif

-   MODES mode;
+   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.
+   unsigned long backends{Backend::CPU}; ///< Bitwise-OR of all configured backends.
   /// 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
+   MemoryType device_mem_type{MemoryType::HOST};    ///< Current Device MemoryType
+   MemoryClass device_mem_class{MemoryClass::HOST};  ///< Current Device MemoryClass
+
+   MemoryType device_temp_mem_type{MemoryType::HOST};    ///< Current Device MemoryType
+   MemoryClass device_temp_mem_class{MemoryClass::HOST};  ///< Current Device MemoryClass

   char *device_option = NULL;
   Device(Device const&);
@@ -173,14 +179,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.
@@ -260,10 +258,24 @@ public:
   /** @deprecated Use GetDeviceMemoryClass() instead. */
   static inline MemoryClass GetMemoryClass() { return Get().device_mem_class; }

+   /** @brief Get the current Device Temporary MemoryType. This is the MemoryType used by
+       MFEM classes when allocating temporary memory to be used with device kernels.
+   */
+   static inline MemoryType GetDeviceTempMemoryType() { return Get().device_temp_mem_type; }
+
+   /** @brief Get the current Device Temporary MemoryClass. This is the MemoryClass used
+       by MFEM device kernels when they need to access temporary Memory objects. */
+   static inline MemoryClass GetDeviceTempMemoryClass() { return Get().device_temp_mem_class; }
+
   static void SetGPUAwareMPI(const bool force = true)
   { Get().mpi_gpu_aware = force; }

   static bool GetGPUAwareMPI() { return Get().mpi_gpu_aware; }
+
+#ifdef MFEM_USE_UMPIRE
+   static bool UseUmpire() { return Get().use_umpire; }
+   static void UseUmpire(bool use) { Get().use_umpire = use; }
+#endif
 };


@@ -67,15 +67,19 @@ MemoryType MemoryManager::GetDualMemoryType_(MemoryType mt)
 {
   switch (mt)
   {
-      case MemoryType::HOST:           return MemoryType::DEVICE;
+      // TODO TMS: temporary
+      case MemoryType::HOST:           return MemoryType::DEVICE_UMPIRE;
      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::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;
+      //case MemoryType::DEVICE_UMPIRE:  return MemoryType::HOST_UMPIRE;
+      case MemoryType::DEVICE_UMPIRE:  return MemoryType::HOST;
+      //case MemoryType::DEVICE_TEMP_UMPIRE:  return MemoryType::HOST_UMPIRE;
+      case MemoryType::DEVICE_TEMP_UMPIRE:  return MemoryType::HOST;
      default: mfem_error("Unknown memory type!");
   }
   MFEM_VERIFY(false,"");
@@ -88,6 +92,9 @@ static void MFEM_VERIFY_TYPES(const MemoryType h_mt, const MemoryType d_mt)
   MFEM_ASSERT(IsDeviceMemory(d_mt),"");
   const bool sync =
      (h_mt == MemoryType::HOST_UMPIRE && d_mt == MemoryType::DEVICE_UMPIRE) ||
+      (h_mt == MemoryType::HOST_UMPIRE && d_mt == MemoryType::DEVICE_TEMP_UMPIRE) ||
+      (h_mt == MemoryType::HOST && d_mt == MemoryType::DEVICE_UMPIRE) ||
+      (h_mt == MemoryType::HOST && d_mt == MemoryType::DEVICE_TEMP_UMPIRE) ||
      (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) ||
@@ -461,48 +468,96 @@ public:
 #ifndef MFEM_USE_UMPIRE
 class UmpireHostMemorySpace : public NoHostMemorySpace { };
 class UmpireDeviceMemorySpace : public NoDeviceMemorySpace { };
+class UmpireDeviceTempMemorySpace : public NoDeviceMemorySpace { };
 #else
+
+// TODO TMS: replace with um.hasAllocatorId(int) when it exists
+bool UmpireHasId(const umpire::ResourceManager & rm, int id)
+{
+   const auto & ids = rm.getAllocatorIds();
+   return std::find(ids.begin(), ids.end(), id) != ids.end();
+}
+
 /// The Umpire host memory space
 class UmpireHostMemorySpace : public HostMemorySpace
 {
 private:
-   const char *name;
   umpire::ResourceManager &rm;
   umpire::Allocator h_allocator;
-   umpire::strategy::AllocationStrategy *strat;
+   bool owns_allocator{false};
 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()) { }
+   // TODO: this only releases unused memory
+   ~UmpireHostMemorySpace() { if (owns_allocator) { h_allocator.release(); } }
+   UmpireHostMemorySpace(): HostMemorySpace(),
+      rm(umpire::ResourceManager::getInstance())
+   {
+      const int id = MemoryManager::GetUmpireHostAllocatorId();
+      if (!UmpireHasId(rm, id))
+      {
+         h_allocator = rm.makeAllocator<umpire::strategy::DynamicPool>("MFEM_HOST",
+                                                                       rm.getAllocator("HOST"));
+         owns_allocator = true;
+      }
+      else
+      {
+         h_allocator = rm.getAllocator(id);
+      }
+      MemoryManager::SetUmpireHostAllocatorId(id);
+   }
   void Alloc(void **ptr, size_t bytes) { *ptr = h_allocator.allocate(bytes); }
   void Dealloc(void *ptr) { h_allocator.deallocate(ptr); }
   void Insert(void *ptr, size_t bytes)
-   { rm.registerAllocation(ptr, {ptr, bytes, strat}); }
+   { mfem_error("UmpireHostMemorySpace::Insert is unsupported"); }
 };

 /// The Umpire device memory space
 #ifdef MFEM_USE_CUDA
-class UmpireDeviceMemorySpace : public DeviceMemorySpace
+class UmpireDeviceMemorySpaceImpl : public DeviceMemorySpace
 {
+public:
+   enum class AllocatorType { TEMPORARY, PERMANENT };
 private:
-   const char *name;
   umpire::ResourceManager &rm;
   umpire::Allocator d_allocator;
+   bool owns_allocator{false};
+
+   int SetupAllocator(int possible_id, const char * allocator_name)
+   {
+      if (!UmpireHasId(rm, possible_id))
+      {
+         d_allocator = rm.makeAllocator<umpire::strategy::DynamicPool>(allocator_name,
+                                                                       rm.getAllocator("DEVICE"));
+         owns_allocator = true;
+      }
+      else
+      {
+         d_allocator = rm.getAllocator(possible_id);
+      }
+
+      return d_allocator.getId();
+   }
 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"))) { }
+   // TODO: this only releases unused memory
+   ~UmpireDeviceMemorySpaceImpl() { if (owns_allocator) { d_allocator.release(); } }
+   UmpireDeviceMemorySpaceImpl(AllocatorType t): DeviceMemorySpace(),
+      rm(umpire::ResourceManager::getInstance())
+   {
+      switch (t)
+      {
+         case AllocatorType::PERMANENT:
+            MemoryManager::SetUmpireDeviceAllocatorId(SetupAllocator(
+                                                         MemoryManager::GetUmpireDeviceAllocatorId(),
+                                                         "MFEM_DEVICE"));
+            break;
+         case AllocatorType::TEMPORARY:
+            MemoryManager::SetUmpireDeviceTempAllocatorId(SetupAllocator(
+                                                             MemoryManager::GetUmpireDeviceTempAllocatorId(),
+                                                             "MFEM_DEVICE_TEMPORARY"));
+            break;
+         default:
+            mfem_error("Unknown Umpire AllocatorType");
+      }
+   }
   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)
@@ -536,8 +591,23 @@ public:
      //rm.copy(dst, const_cast<void*>(src), bytes); return dst;
   }
 };
+
+class UmpireDeviceMemorySpace : public UmpireDeviceMemorySpaceImpl
+{
+public:
+   UmpireDeviceMemorySpace() : UmpireDeviceMemorySpaceImpl(
+         AllocatorType::PERMANENT) {}
+};
+
+class UmpireDeviceTempMemorySpace : public UmpireDeviceMemorySpaceImpl
+{
+public:
+   UmpireDeviceTempMemorySpace() : UmpireDeviceMemorySpaceImpl(
+         AllocatorType::TEMPORARY) {}
+};
 #else
 class UmpireDeviceMemorySpace : public NoDeviceMemorySpace { };
+class UmpireDeviceTempMemorySpace : public NoDeviceMemorySpace { };
 #endif // MFEM_USE_CUDA
 #endif // MFEM_USE_UMPIRE

@@ -568,7 +638,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
@@ -610,8 +680,12 @@ 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();
+         case MT::HOST_UMPIRE: return new UmpireHostMemorySpace();
+         default: MFEM_ABORT("Unknown host memory controller!");
+      }
      return nullptr;
   }

@@ -620,6 +694,7 @@ private:
      switch (mt)
      {
         case MT::DEVICE_UMPIRE: return new UmpireDeviceMemorySpace();
+         case MT::DEVICE_TEMP_UMPIRE: return new UmpireDeviceTempMemorySpace();
         case MT::DEVICE_DEBUG: return new MmuDeviceMemorySpace();
         case MT::DEVICE:
         {
@@ -760,7 +835,7 @@ 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);
@@ -783,6 +858,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_TEMP_UMPIRE ||
                     d_mt == MemoryType::MANAGED,"");
         return true;
      }
@@ -1262,22 +1338,15 @@ MemoryManager::MemoryManager() { Init(); }
 MemoryManager::~MemoryManager() { if (exists) { Destroy(); } }

 void MemoryManager::Configure(const MemoryType host_mt,
-                              const MemoryType device_mt)
+                              const MemoryType device_mt,
+                              const MemoryType device_tmt)
 {
   Init();
   host_mem_type = host_mt;
   device_mem_type = device_mt;
+   device_temp_mem_type = device_tmt;
 }

-#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!");
@@ -1381,12 +1450,14 @@ MemoryManager mm;
 bool MemoryManager::exists = false;

 #ifdef MFEM_USE_UMPIRE
-const char* MemoryManager::h_umpire_name = "HOST";
-const char* MemoryManager::d_umpire_name = "DEVICE";
+int MemoryManager::h_umpire_id = -1;
+int MemoryManager::d_umpire_id = -1;
+int MemoryManager::d_umpire_temp_id = -1;
 #endif

 MemoryType MemoryManager::host_mem_type = MemoryType::HOST;
 MemoryType MemoryManager::device_mem_type = MemoryType::HOST;
+MemoryType MemoryManager::device_temp_mem_type = MemoryType::HOST;

 const char *MemoryTypeName[MemoryTypeSize] =
 {
@@ -1403,11 +1474,14 @@ const char *MemoryTypeName[MemoryTypeSize] =
 #endif
   "device-debug",
 #if defined(MFEM_USE_CUDA)
-   "cuda-umpire"
+   "cuda-umpire",
+   "cuda-umpire-temp"
 #elif defined(MFEM_USE_HIP)
-   "hip-umpire"
+   "hip-umpire",
+   "hip-umpire-temp"
 #else
-   "device-umpire"
+   "device-umpire",
+   "device-umpire-temp"
 #endif
 };

@@ -27,18 +27,19 @@ namespace mfem
 /// Memory types supported by MFEM.
 enum class MemoryType
 {
-   HOST,           ///< Host memory; using new[] and delete[]
-   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
-   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
+   HOST,                ///< Host memory; using new[] and delete[]
+   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
+   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
+   DEVICE_TEMP_UMPIRE,  ///< Temporary Device memory; using Umpire
+   SIZE                 ///< Number of host and device memory types
 };

 /// Static casts to 'int' and sizes of some useful memory types.
@@ -61,7 +62,7 @@ enum class MemoryClass
                                 HOST_UMPIRE, 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_TEMP_UMPIRE, MANAGED }
   MANAGED  ///< Memory types: { MANAGED }
 };

@@ -450,16 +451,16 @@ private:

   template <std::size_t align_bytes, bool dummy = true> struct Alloc
   {
+#if __cplusplus < 201703L
      static inline T *New(std::size_t)
      {
-#if __cplusplus < 201703L
         // Generate an error in debug mode
         MFEM_ASSERT(false, "overaligned type cannot use MemoryType::HOST");
         return nullptr;
-#else
-         return new T[size];
-#endif
      }
+#else
+      static inline T *New(std::size_t size) { return new T[size]; }
+#endif
   };

 #if __cplusplus < 201703L
@@ -489,6 +490,9 @@ private:
   /// Device memory type set during the Setup.
   static MemoryType device_mem_type;

+   /// Device temporary memory type set during the Setup.
+   static MemoryType device_temp_mem_type;
+
   /// Allow to detect if a global memory manager instance exists.
   static bool exists;

@@ -497,8 +501,9 @@ private:

   /// Host and device allocator names for Umpire.
 #ifdef MFEM_USE_UMPIRE
-   static const char *h_umpire_name;
-   static const char *d_umpire_name;
+   static int h_umpire_id;
+   static int d_umpire_id;
+   static int d_umpire_temp_id;
 #endif

 private: // Static methods used by the Memory<T> class
@@ -621,15 +626,21 @@ public:
   /// Initialize the memory manager.
   void Init();

-   /// Configure the Memory manager with given default host and device types
+   /// Configure the Memory manager with given default host, device, and device temporary types
   /// This method will be called when configuring a device.
-   void Configure(const MemoryType h_mt, const MemoryType d_mt);
+   void Configure(const MemoryType h_mt, const MemoryType d_mt,
+                  const MemoryType d_tmt);

 #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 and device Umpire allocator ids
+   static void SetUmpireHostAllocatorId(int h_id) { h_umpire_id = h_id; }
+   static void SetUmpireDeviceAllocatorId(int d_id) { d_umpire_id = d_id; }
+   static void SetUmpireDeviceTempAllocatorId(int d_id) { d_umpire_temp_id = d_id; }
+
+   /// Get the host and device Umpire allocator ids
+   static int GetUmpireHostAllocatorId() { return h_umpire_id; }
+   static int GetUmpireDeviceAllocatorId() { return d_umpire_id; }
+   static int GetUmpireDeviceTempAllocatorId() { return d_umpire_temp_id; }
 #endif

   /// Free all the device memories
@@ -654,6 +665,7 @@ public:

   static MemoryType GetHostMemoryType() { return host_mem_type; }
   static MemoryType GetDeviceMemoryType() { return device_mem_type; }
+   static MemoryType GetDeviceTempMemoryType() { return device_temp_mem_type; }
 };


@@ -3595,4 +3595,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
@@ -884,6 +884,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

@@ -28,6 +28,25 @@ 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
+}
+
 SparseMatrix::SparseMatrix(int nrows, int ncols)
   : AbstractSparseMatrix(nrows, (ncols >= 0) ? ncols : nrows),
     Rows(new RowNode *[nrows]),
@@ -50,6 +69,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 +88,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 +121,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 +144,11 @@ 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 +156,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 +168,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;
@@ -184,6 +211,8 @@ SparseMatrix::SparseMatrix(const SparseMatrix &mat, bool copy_graph)
   ColPtrNode = NULL;
   At = NULL;
   isSorted = mat.isSorted;
+
+   InitCuSparse();
 }

 SparseMatrix::SparseMatrix(const Vector &v)
@@ -211,6 +240,8 @@ SparseMatrix::SparseMatrix(const Vector &v)
      J[r] = r;
      A[r] = v[r];
   }
+
+   InitCuSparse();
 }

 SparseMatrix& SparseMatrix::operator=(const SparseMatrix &rhs)
@@ -250,6 +281,16 @@ void SparseMatrix::SetEmpty()
   NodesMem = NULL;
 #endif
   isSorted = false;
+
+#ifdef MFEM_USE_CUDA
+   if (initBuffers)
+   {
+      cusparseDestroySpMat(matA_descr);
+      cusparseDestroyDnVec(vecX_descr);
+      cusparseDestroyDnVec(vecY_descr);
+      initBuffers = false;
+   }
+#endif
 }

 int SparseMatrix::RowSize(const int i) const
@@ -494,24 +535,29 @@ void SparseMatrix::GetDiag(Vector & d) const

   d.SetSize(height);

-   int j, end;
-   for (int i = 0; i < height; i++)
-   {
+   auto I = this->ReadI();
+   auto J = this->ReadJ();
+   auto A = this->ReadData();
+   auto dd = d.Write();

-      end = I[i+1];
-      for (j = I[i]; j < end; j++)
+   MFEM_FORALL(i, height,
+   {
+      const int begin = I[i];
+      const int end = I[i+1];
+      int j;
+      for (j = begin; j < end; j++)
      {
         if (J[j] == i)
         {
-            d[i] = A[j];
+            dd[i] = A[j];
            break;
         }
      }
      if (j == end)
      {
-         d[i] = 0.;
+         dd[i] = 0.;
      }
-   }
+   });
 }

 /// Produces a DenseMatrix from a SparseMatrix
@@ -587,16 +633,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();
@@ -2145,31 +2247,46 @@ void SparseMatrix::DiagScale(const Vector &b, Vector &x, double sc) const
 {
   MFEM_VERIFY(Finalized(), "Matrix must be finalized.");

+   const int nnz = J.Capacity();
+
+   const bool use_dev = b.UseDevice() || x.UseDevice();
+
+   auto bp = b.Read(use_dev);
+   auto xp = x.Write(use_dev);
+
+   auto Ap = Read(A, nnz);
+   auto Ip = Read(I, height+1);
+   auto Jp = Read(J, nnz);
+
   bool scale = (sc != 1.0);
-   for (int i = 0, j = 0; i < height; i++)
+   MFEM_FORALL(i, height,
   {
-      int end = I[i+1];
-      for ( ; true; j++)
+      int end = Ip[i+1];
+      for (int j = Ip[i]; true; j++)
      {
-         MFEM_VERIFY(j != end, "Couldn't find diagonal in row. i = " << i
-                     << ", j = " << j
-                     << ", I[i+1] = " << end );
-         if (J[j] == i)
+         if (j == end)
         {
-            MFEM_VERIFY(std::abs(A[j]) > 0.0, "Diagonal " << j << " must be nonzero");
+            //MFEM_ABORT_KERNEL("Diagonal not found in SparseMatrix::DiagScale");
+         }
+         if (Jp[j] == i)
+         {
+            if (!(std::abs(Ap[j]) > 0.0))
+            {
+               //MFEM_ABORT_KERNEL("Zero diagonal in SparseMatrix::DiagScale");
+            }
+
            if (scale)
            {
-               x(i) = sc * b(i) / A[j];
+               xp[i] = sc * bp[i] / Ap[j];
            }
            else
            {
-               x(i) = b(i) / A[j];
+               xp[i] = bp[i] / Ap[j];
            }
            break;
         }
      }
-      j = end;
-   }
+   });
   return;
 }

@@ -2749,6 +2866,10 @@ void SparseMatrix::Print(std::ostream & out, int _width) const
      return;
   }

+   // HostRead forces synchronization
+   HostReadI();
+   HostReadJ();
+   HostReadData();
   for (i = 0; i < height; i++)
   {
      out << "[row " << i << "]\n";
@@ -2938,6 +3059,16 @@ void SparseMatrix::Destroy()
   delete NodesMem;
 #endif
   delete At;
+
+#ifdef MFEM_USE_CUDA
+   if (initBuffers)
+   {
+      cusparseDestroySpMat(matA_descr);
+      cusparseDestroyDnVec(vecX_descr);
+      cusparseDestroyDnVec(vecY_descr);
+      initBuffers = false;
+   }
+#endif
 }

 int SparseMatrix::ActualWidth() const
@@ -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,33 @@ 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 +143,14 @@ 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);
@@ -573,11 +606,28 @@ 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; }
 };

+inline std::ostream& operator<<(std::ostream& os, SparseMatrix const& mat)
+{
+   mat.Print(os);
+   return os;
+}
+
 /// Applies f() to each element of the matrix (after it is finalized).
 void SparseMatrixFunction(SparseMatrix &S, double (*f)(double));

@@ -911,20 +911,27 @@ FaceElementTransformations *Mesh::GetFaceElementTransformations(int FaceNo,
      // NC meshes: prepend slave edge/face transformation to Loc2
      if (Nonconforming() && IsSlaveFace(face_info))
      {
-         ApplyLocalSlaveTransformation(FaceElemTr.Loc2.Transf, face_info);
-
-         if (face_type == Element::SEGMENT)
-         {
-            // flip Loc2 to match Loc1 and Face
-            DenseMatrix &pm = FaceElemTr.Loc2.Transf.GetPointMat();
-            std::swap(pm(0,0), pm(0,1));
-            std::swap(pm(1,0), pm(1,1));
-         }
+         ApplyLocalSlaveTransformation(FaceElemTr, face_info, false);
      }
   }

   FaceElemTr.SetConfigurationMask(mask);

+   // This check can be useful for internal debugging, however it will fail on
+   // periodic boundary faces, so we keep it disabled in general.
+#if 0
+#ifdef MFEM_DEBUG
+   double dist = FaceElemTr.CheckConsistency();
+   if (dist >= 1e-12)
+   {
+      mfem::out << "\nInternal error: face id = " << FaceNo
+                << ", dist = " << dist << '\n';
+      FaceElemTr.CheckConsistency(1); // print coordinates
+      MFEM_ABORT("internal error");
+   }
+#endif
+#endif
+
   return &FaceElemTr;
 }

@@ -933,8 +940,8 @@ bool Mesh::IsSlaveFace(const FaceInfo &fi) const
   return fi.NCFace >= 0 && nc_faces_info[fi.NCFace].Slave;
 }

-void Mesh::ApplyLocalSlaveTransformation(IsoparametricTransformation &transf,
-                                         const FaceInfo &fi)
+void Mesh::ApplyLocalSlaveTransformation(FaceElementTransformations &FT,
+                                         const FaceInfo &fi, bool is_ghost)
 {
 #ifdef MFEM_THREAD_SAFE
   DenseMatrix composition;
@@ -942,8 +949,34 @@ void Mesh::ApplyLocalSlaveTransformation(IsoparametricTransformation &transf,
   static DenseMatrix composition;
 #endif
   MFEM_ASSERT(fi.NCFace >= 0, "");
-   transf.Transform(*nc_faces_info[fi.NCFace].PointMatrix, composition);
-   transf.SetPointMat(composition);
+   MFEM_ASSERT(nc_faces_info[fi.NCFace].Slave, "internal error");
+   if (!is_ghost)
+   {
+      // side 1 -> child side, side 2 -> parent side
+      IsoparametricTransformation &LT = FT.Loc2.Transf;
+      LT.Transform(*nc_faces_info[fi.NCFace].PointMatrix, composition);
+      // In 2D, we need to flip the point matrix since it is aligned with the
+      // parent side.
+      if (Dim == 2)
+      {
+         // swap points (columns) 0 and 1
+         std::swap(composition(0,0), composition(0,1));
+         std::swap(composition(1,0), composition(1,1));
+      }
+      LT.SetPointMat(composition);
+   }
+   else // is_ghost == true
+   {
+      // side 1 -> parent side, side 2 -> child side
+      IsoparametricTransformation &LT = FT.Loc1.Transf;
+      LT.Transform(*nc_faces_info[fi.NCFace].PointMatrix, composition);
+      // In 2D, there is no need to flip the point matrix since it is already
+      // aligned with the parent side, see also ParNCMesh::GetFaceNeighbors.
+      // In 3D the point matrix was flipped during construction in
+      // ParNCMesh::GetFaceNeighbors and due to that it is already aligned with
+      // the parent side.
+      LT.SetPointMat(composition);
+   }
 }

 FaceElementTransformations *Mesh::GetBdrFaceTransformations(int BdrElemNo)
@@ -988,7 +1021,21 @@ void Mesh::GetFaceInfos(int Face, int *Inf1, int *Inf2) const

 Geometry::Type Mesh::GetFaceGeometryType(int Face) const
 {
-   return (Dim == 1) ? Geometry::POINT : faces[Face]->GetGeometryType();
+   switch (Dim)
+   {
+      case 1: return Geometry::POINT;
+      case 2: return Geometry::SEGMENT;
+      case 3:
+         if (Face < NumOfFaces) // local (non-ghost) face
+         {
+            return faces[Face]->GetGeometryType();
+         }
+         // ghost face
+         const int nc_face_id = faces_info[Face].NCFace;
+         MFEM_ASSERT(nc_face_id >= 0, "parent ghost faces are not supported");
+         return faces[nc_faces_info[nc_face_id].MasterFace]->GetGeometryType();
+   }
+   return Geometry::INVALID;
 }

 Element::Type Mesh::GetFaceElementType(int Face) const
@@ -5261,7 +5308,10 @@ void Mesh::GenerateNCFaceInfo()

      slave_fi.Elem2No = master_fi.Elem1No;
      slave_fi.Elem2Inf = 64 * master_nc.MasterFace; // get lf no. stored above
-      // NOTE: orientation part of Elem2Inf is encoded in the point matrix
+      // NOTE: In 3D, the orientation part of Elem2Inf is encoded in the point
+      //       matrix. In 2D, the point matrix has the orientation of the parent
+      //       edge, so its columns need to be flipped when applying it, see
+      //       ApplyLocalSlaveTransformation.
   }
 }

@@ -10474,17 +10524,17 @@ GeometricFactors::GeometricFactors(const Mesh *mesh, const IntegrationRule &ir,
   unsigned eval_flags = 0;
   if (flags & GeometricFactors::COORDINATES)
   {
-      X.SetSize(vdim*NQ*NE);
+      X.SetSize(vdim*NQ*NE, Device::GetDeviceTempMemoryType());
      eval_flags |= QuadratureInterpolator::VALUES;
   }
   if (flags & GeometricFactors::JACOBIANS)
   {
-      J.SetSize(dim*vdim*NQ*NE);
+      J.SetSize(dim*vdim*NQ*NE, Device::GetDeviceTempMemoryType());
      eval_flags |= QuadratureInterpolator::DERIVATIVES;
   }
   if (flags & GeometricFactors::DETERMINANTS)
   {
-      detJ.SetSize(NQ*NE);
+      detJ.SetSize(NQ*NE, Device::GetDeviceTempMemoryType());
      eval_flags |= QuadratureInterpolator::DETERMINANTS;
   }

@@ -10494,7 +10544,7 @@ GeometricFactors::GeometricFactors(const Mesh *mesh, const IntegrationRule &ir,
   qi->SetOutputLayout(QVectorLayout::byNODES);
   if (elem_restr)
   {
-      Vector Enodes(vdim*ND*NE);
+      Vector Enodes(vdim*ND*NE, Device::GetDeviceTempMemoryType());
      elem_restr->Mult(*nodes, Enodes);
      qi->Mult(Enodes, eval_flags, X, J, detJ);
   }
@@ -131,11 +131,16 @@ protected:
   //   face. Elem2No is < 0 and -1-Elem2No is the index of the ghost
   //   face-neighbor element that generated this slave ghost face. In this
   //   case, Elem2Inf >= 0.
+   // Relevant methods: GenerateFaces(), GenerateNCFaceInfo(),
+   //                   ParNCMesh::GetFaceNeighbors(),
+   //                   ParMesh::ExchangeFaceNbrData()

   struct NCFaceInfo
   {
      bool Slave; // true if this is a slave face, false if master face
      int MasterFace; // if Slave, this is the index of the master face
+      // If not Slave, 'MasterFace' is the local face index of this master face
+      // as a face in the unique adjacent element.
      const DenseMatrix* PointMatrix; // if Slave, position within master face
      // (NOTE: PointMatrix points to a matrix owned by NCMesh.)

@@ -363,8 +368,9 @@ protected:

   /** Used in GetFaceElementTransformations to account for the fact that a
       slave face occupies only a portion of its master face. */
-   void ApplyLocalSlaveTransformation(IsoparametricTransformation &transf,
-                                      const FaceInfo &fi);
+   void ApplyLocalSlaveTransformation(FaceElementTransformations &FT,
+                                      const FaceInfo &fi, bool is_ghost);
+
   bool IsSlaveFace(const FaceInfo &fi) const;

   /// Returns the orientation of "test" relative to "base"
@@ -2363,16 +2363,16 @@ Table *ParMesh::GetFaceToAllElementTable() const
   return face_elem;
 }

-ElementTransformation* ParMesh::GetGhostFaceTransformation(
+void ParMesh::GetGhostFaceTransformation(
   FaceElementTransformations* FETr, Element::Type face_type,
   Geometry::Type face_geom)
 {
   // calculate composition of FETr->Loc1 and FETr->Elem1
-   DenseMatrix &face_pm = FaceTransformation.GetPointMat();
+   DenseMatrix &face_pm = FETr->GetPointMat();
   if (Nodes == NULL)
   {
      FETr->Elem1->Transform(FETr->Loc1.Transf.GetPointMat(), face_pm);
-      FaceTransformation.SetFE(GetTransformationFEforElementType(face_type));
+      FETr->SetFE(GetTransformationFEforElementType(face_type));
   }
   else
   {
@@ -2388,15 +2388,23 @@ ElementTransformation* ParMesh::GetGhostFaceTransformation(
      FETr->Loc1.Transform(face_el->GetNodes(), eir);
      Nodes->GetVectorValues(*FETr->Elem1, eir, face_pm);
 #endif
-      FaceTransformation.SetFE(face_el);
+      FETr->SetFE(face_el);
   }
-   return &FaceTransformation;
 }

 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];

@@ -2453,22 +2461,10 @@ GetSharedFaceTransformations(int sf, bool fill2)
   // adjust Loc1 or Loc2 of the master face if this is a slave face
   if (is_slave)
   {
-      // is a ghost slave? -> master not a ghost -> choose Elem1 local transf
-      // not a ghost slave? -> master is a ghost -> choose Elem2 local transf
-      IsoparametricTransformation &loctr =
-         is_ghost ? FaceElemTr.Loc1.Transf : FaceElemTr.Loc2.Transf;
-
      if (is_ghost || fill2)
      {
-         ApplyLocalSlaveTransformation(loctr, face_info);
-      }
-
-      if (face_type == Element::SEGMENT && fill2)
-      {
-         // fix slave orientation in 2D: flip Loc2 to match Loc1 and Face
-         DenseMatrix &pm = FaceElemTr.Loc2.Transf.GetPointMat();
-         std::swap(pm(0,0), pm(0,1));
-         std::swap(pm(1,0), pm(1,1));
+         // is_ghost -> modify side 1, otherwise -> modify side 2:
+         ApplyLocalSlaveTransformation(FaceElemTr, face_info, is_ghost);
      }
   }

@@ -2478,6 +2474,23 @@ GetSharedFaceTransformations(int sf, bool fill2)
      GetGhostFaceTransformation(&FaceElemTr, face_type, face_geom);
   }

+   FaceElemTr.SetConfigurationMask(fill2 ? 31 : 21);
+
+   // This check can be useful for internal debugging, however it will fail on
+   // periodic boundary faces, so we keep it disabled in general.
+#if 0
+#ifdef MFEM_DEBUG
+   double dist = FaceElemTr.CheckConsistency();
+   if (dist >= 1e-12)
+   {
+      mfem::out << "\nInternal error: face id = " << FaceNo
+                << ", dist = " << dist << ", rank = " << MyRank << '\n';
+      FaceElemTr.CheckConsistency(1); // print coordinates
+      MFEM_ABORT("internal error");
+   }
+#endif
+#endif
+
   return &FaceElemTr;
 }

@@ -107,7 +107,7 @@ protected:
   void GetFaceNbrElementTransformation(
      int i, IsoparametricTransformation *ElTr);

-   ElementTransformation* GetGhostFaceTransformation(
+   void GetGhostFaceTransformation(
      FaceElementTransformations* FETr, Element::Type face_type,
      Geometry::Type face_geom);

@@ -291,9 +291,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);

   /// Return the number of shared faces (3D), edges (2D), vertices (1D)
   int GetNSharedFaces() const;
@@ -1263,6 +1263,8 @@ void ParNCMesh::GetFaceNeighbors(ParMesh &pmesh)
            const DenseMatrix* pm = &sf.point_matrix;
            if (!sloc && Dim == 3)
            {
+               // TODO: does this handle triangle faces correctly?
+
               // ghost slave in 3D needs flipping orientation
               DenseMatrix* pm2 = new DenseMatrix(*pm);
               std::swap((*pm2)(0,1), (*pm2)(0,3));
@@ -1282,6 +1284,14 @@ void ParNCMesh::GetFaceNeighbors(ParMesh &pmesh)
               // processor, but on the other it is the element containing the
               // master face. In the latter case we need to flip the pm.
            }
+            else if (!sloc && Dim == 2)
+            {
+               fi.Elem2Inf ^= 1; // set orientation to 1
+               // The point matrix (used to define "side 1" which is the same as
+               // "parent side" in this case) does not require a flip since it
+               // is aligned with the parent side, so NO flip is performed in
+               // Mesh::ApplyLocalSlaveTransformation.
+            }

            MFEM_ASSERT(fi.NCFace < 0, "");
            fi.NCFace = pmesh.nc_faces_info.Size();
@@ -32,6 +32,16 @@ public:
      : attr(attr), coeff(coeff)
   {}

+   VelDirichletBC_T(VelDirichletBC_T &&obj)
+   {
+      // Deep copy the attribute array
+      this->attr = obj.attr;
+
+      // Move the coefficient pointer
+      this->coeff = obj.coeff;
+      obj.coeff = nullptr;
+   }
+
   ~VelDirichletBC_T() { delete coeff; }

   Array<int> attr;
@@ -46,6 +56,16 @@ public:
      : attr(attr), coeff(coeff)
   {}

+   PresDirichletBC_T(PresDirichletBC_T &&obj)
+   {
+      // Deep copy the attribute array
+      this->attr = obj.attr;
+
+      // Move the coefficient pointer
+      this->coeff = obj.coeff;
+      obj.coeff = nullptr;
+   }
+
   ~PresDirichletBC_T() { delete coeff; }

   Array<int> attr;
@@ -60,6 +80,16 @@ public:
      : attr(attr), coeff(coeff)
   {}

+   AccelTerm_T(AccelTerm_T &&obj)
+   {
+      // Deep copy the attribute array
+      this->attr = obj.attr;
+
+      // Move the coefficient pointer
+      this->coeff = obj.coeff;
+      obj.coeff = nullptr;
+   }
+
   ~AccelTerm_T() { delete coeff; }

   Array<int> attr;
Author	SHA1	Message	Date
Tom Stitt	0e8625cf89	add optional MemoryType to the SparseMatrix copy-ctor	2020-07-24 09:13:21 -07:00
Tom Stitt	3f7a432e74	Cherry-picking changes on top of our current hash. There are things that hurt performance in laster mfem@master commits, not sure what yet. This includes: - feature/artv3/cusparse-Spmv - feature/tomstitt/temp-mem-type - patches to tmop.cpp and tmop_tools.cpp to address memory issues - patch to mem_manager/device to use mfem's default allcator instead of the host umpire one	2020-07-22 15:51:04 -07:00
Tom Stitt	a57fd02a4c	Use flattened FORALL instead of FORALL_2D.	2020-06-23 12:11:33 -07:00
Tom Stitt	2ae97ff2da	Converts PADiffusionSetup3D and QuadratureInterpolator::Eval3D kernels from 1 element per thread to 1 qpt/dof per thread for better performance when offloading (there are not enough units of work with 1 element/thread)	2020-06-22 14:42:48 -07:00
Arturo Vargas	2b712207c6	fix merge conflicts	2020-06-22 14:29:28 -07:00
Tzanio Kolev	62e95fe7bc	Merge pull request #1564 from mfem/nc-amr-hybridization-fix Fix a bug with hybridization on 2D NC meshes	2020-06-22 12:47:31 -07:00
Veselin Dobrev	c62681aaf8	Merge branch 'master' into nc-amr-hybridization-fix	2020-06-19 21:07:07 -07:00
Veselin Dobrev	a25138e611	Add an optional output stream parameter to the new method FaceElementTransformations::CheckConsistency and improved its documentation. Extended Mesh::GetFaceGeometryType to support ghost child faces and use this to generalize ParFiniteElementSpace::GetFaceNbrFaceFE to support all face types. Cleanup some old and debug code.	2020-06-19 20:43:47 -07:00
Tzanio Kolev	cc3a4475d2	Merge pull request #1460 from mfem/LF-grad-curl-div Additional LF-integrators	2020-06-18 18:52:17 -07:00
Tzanio	8e00ae03e4	Removed convergence/bae.cpp, this will be handled in another PR.	2020-06-18 18:49:12 -07:00
Arturo Vargas	2930c1477f	promove temp vector class member	2020-06-18 13:35:11 -07:00
Tzanio Kolev	bf3423dcba	Merge pull request #1497 from mfem/coef-dev Adding accessor methods to compound coefficient classes [coef-dev]	2020-06-18 10:55:33 -07:00
Veselin Dobrev	dd23ccddb3	In ex4p, added a sample run that tests hybridization in 2D when there are ghost shared faces using amr-quad.mesh on 3 processors. In ex9p, added a sample run on amr-hex.mesh on 3 processors so have a test case with ghost shared faces in 3D.	2020-06-17 22:39:16 -07:00
Tzanio Kolev	3c5dc275a4	Merge pull request #1539 from mfem/navier-bookkeeping-fix Navier miniapp bookkeeping fix	2020-06-17 16:37:23 -07:00
Veselin Dobrev	6604ba702c	WIP: debugging and bugfix for the issue with hybridization on 2D nonconforming meshes, see issue #1105 on github.	2020-06-16 21:09:48 -07:00
Julian Andrej	395e75f69a	Merge branch 'master' into navier-bookkeeping-fix	2020-06-12 11:21:38 -07:00
Julian Andrej	1875c34055	Add move constructor to custom types	2020-06-12 11:20:07 -07:00
Stowell, Mark L	f58decb421	Adding Doxygen comments for the new accessor methods	2020-06-09 13:51:46 -07:00
Stowell, Mark L	635ac55e77	Creating Doxygen group for compound coefficient	2020-06-09 09:41:44 -07:00
Stowell, Mark L	6beb0decd7	Merge remote-tracking branch 'origin/master' into coef-dev # Conflicts: # fem/coefficient.hpp	2020-06-09 09:41:13 -07:00
Tzanio Kolev	253aebb734	Merge branch 'master' into LF-grad-curl-div	2020-06-09 07:28:20 -07:00
Stowell, Mark L	7de36b40bb	Wrapping math formulae in \f$	2020-06-08 12:27:59 -07:00
Stowell, Mark L	bf54fc6b3f	Expanding name of MatVecCoefficient to conform to standard set by other product coefficients	2020-06-07 13:59:16 -07:00
Stowell, Mark L	fb3b943cf0	Adding/improving comments	2020-06-07 13:53:24 -07:00
psocratis	df3d0b51c2	renamed example name	2020-06-02 15:54:36 -07:00
psocratis	4c4aeaeec3	renamed example	2020-06-02 15:51:28 -07:00
psocratis	cadb9778a4	Added the legal statement at the top. Corrected the mesh path	2020-05-28 15:14:43 -07:00
Socratis	32bdcb6ef2	make style	2020-05-27 19:01:16 -07:00
Socratis	e74d501ac3	fixed comments in example	2020-05-27 19:00:11 -07:00
Socratis	fcd2f0ceb2	Merge branch 'master' into LF-grad-curl-div	2020-05-27 18:44:36 -07:00
Socratis	c0879aaf76	fixed valgrind leaks	2020-05-27 18:44:06 -07:00
Stowell, Mark L	e024ac0f9c	Adding accessor methods to compound coefficient classes (and adding a couple new ones)	2020-05-20 15:43:17 -07:00
Arturo Vargas	b3ee631aa6	tmop_tools	2020-05-14 14:33:32 -07:00
Arturo Vargas	daf2fdecec	use tmop_tools from master	2020-05-14 14:24:52 -07:00
Arturo Vargas	d3a0d0a181	Merge branch 'master' into feature/apps/mesh-helper-tmop-ea	2020-05-14 14:19:13 -07:00
Arturo Vargas	07853b9c62	add missing host reads to tmop	2020-05-09 21:41:26 -07:00
Socratis	20089dd62b	valgrind checks passed	2020-05-05 15:36:16 -07:00
psocratis	a8ebdf16bd	Added convergence tests	2020-05-05 15:19:27 -07:00
psocratis	99372de933	Added parallel test	2020-05-04 20:35:27 -07:00
psocratis	70d84e441d	make style	2020-05-04 19:56:38 -07:00
psocratis	f6bc0f8988	Added test for the newly added LF integrators	2020-05-04 19:55:00 -07:00
psocratis	ef54519a29	Added LF integrators (Q, grard V), (Q, curl V) and (Q, div V)	2020-05-04 19:54:24 -07:00
Arturo Vargas	49a31c0cf7	Merge branch 'yohann/mass-emat' into feature/apps/mesh-helper-tmop-ea	2020-05-01 10:17:42 -07:00
Arturo Vargas	ba9b251007	Merge branch 'master' into feature/apps/mesh-helper-tmop-ea	2020-05-01 10:13:46 -07:00
Arturo Vargas	5ce2fa9ab9	fix free	2020-04-30 18:06:29 -07:00
Arturo Vargas	78c93de6ce	Merge branch 'gpu/artv3/batch-lu' into feature/apps/mesh-helper-tmop-ea	2020-04-30 13:09:53 -07:00
Arturo Vargas	3d4aa157cb	testing super code	2020-04-28 16:54:23 -07:00
Tomov	a7f1c177c5	Fixed a parallel communication bug in tmop's AdvectorCG.	2020-03-30 22:29:51 -07:00
Arturo Vargas	2d4e3cf77e	minor helper functions	2020-03-24 14:11:37 -07:00