add optional MemoryType to the SparseMatrix copy-ctor

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

config/defaults.mk

@@ -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 =

examples/ex4p.cpp

@@ -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

examples/ex9p.cpp

@@ -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

fem/bilininteg_diffusion_pa.cpp

@@ -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
    });
 }
 

fem/bilininteg_mass_pa.cpp

@@ -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)

fem/coefficient.cpp

@@ -12,6 +12,7 @@
 // Implementation of Coefficient class
 
 #include "fem.hpp"
+#include "../linalg/dtensor.hpp"
 
 #include <cmath>
 #include <limits>
@@ -21,6 +22,13 @@ namespace mfem
 
 using namespace std;
 
+double QuadratureCoefficient::Eval(ElementTransformation & T,
+                                   const IntegrationPoint & ip)
+{
+   auto coeff = mfem::Reshape(qData->HostRead(), nip, NE);
+   return coeff(ip.index, T.ElementNo);
+}
+
 double PWConstCoefficient::Eval(ElementTransformation & T,
                                 const IntegrationPoint & ip)
 {
@@ -285,22 +293,6 @@ void VectorRestrictedCoefficient::Eval(
    }
 }
 
-void UnitNormalCoefficient::Eval(Vector &V, ElementTransformation &T,
-                                 const IntegrationPoint &ip)
-{
-   V.SetSize(vdim);
-   V = 0.0;
-
-   const DenseMatrix & J = T.Jacobian();
-   if (J.Width() == J.Height() - 1)
-   {
-      CalcOrtho(J, V);
-      double norm = V.Norml2();
-      MFEM_ASSERT(norm > 0.0, "Length of normal vector is non-positive!");
-      V /= norm;
-   }
-}
-
 void MatrixFunctionCoefficient::Eval(DenseMatrix &K, ElementTransformation &T,
                                      const IntegrationPoint &ip)
 {
@@ -432,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.");
 }
@@ -446,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)
@@ -487,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);
@@ -533,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;
 }
@@ -597,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[])
 {

fem/coefficient.hpp

@@ -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
@@ -685,16 +712,6 @@ public:
                      const IntegrationRule &ir);
 };
 
-/// VectorCoefficient which computes unit normal vector on the mesh boundary
-class UnitNormalCoefficient : public VectorCoefficient
-{
-public:
-   UnitNormalCoefficient(int dim) : VectorCoefficient(dim) {}
-
-   virtual void Eval(Vector &V, ElementTransformation &T,
-                     const IntegrationPoint &ip);
-};
-
 
 /// Base class for Matrix Coefficients that optionally depend on time and space.
 class MatrixCoefficient
@@ -862,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;
 
@@ -875,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
@@ -917,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)
@@ -937,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);
@@ -953,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);
@@ -973,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;
@@ -991,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);
@@ -1005,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);
@@ -1032,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);
@@ -1040,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;
@@ -1050,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,
@@ -1059,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
 {
@@ -1075,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:
@@ -1092,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);
@@ -1102,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);
@@ -1124,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);
@@ -1139,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);
@@ -1158,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;
 

fem/eltrans.cpp

@@ -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;
+}
+
 }

fem/eltrans.hpp

@@ -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))

fem/gridfunc.cpp

@@ -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);
    }

fem/lininteg.cpp

@@ -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)
 {

fem/lininteg.hpp

@@ -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)

fem/pfespace.cpp

@@ -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()

fem/pgridfunc.cpp

@@ -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(),

fem/pgridfunc.hpp

@@ -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);
 

fem/quadinterpolator.cpp

@@ -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)
          {

fem/quadinterpolator_face.cpp

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

fem/tmop.cpp

@@ -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++)

fem/tmop_tools.cpp

@@ -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; }

general/device.cpp

@@ -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()

general/device.hpp

@@ -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
 };
 
 

general/mem_manager.cpp

@@ -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
 };
 

general/mem_manager.hpp

@@ -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; }
 };
 
 

linalg/densemat.cpp

@@ -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

linalg/densemat.hpp

@@ -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
 

linalg/sparsemat.cpp

@@ -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

linalg/sparsemat.hpp

@@ -21,6 +21,12 @@
 #include "../general/globals.hpp"
 #include "densemat.hpp"
 
+#ifdef MFEM_USE_CUDA
+#include <cusparse.h>
+#include <library_types.h>
+#include "../general/cuda.hpp"
+#endif
+
 namespace mfem
 {
 
@@ -80,9 +86,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));
 

mesh/mesh.cpp

@@ -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);
    }

mesh/mesh.hpp

@@ -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"

mesh/pmesh.cpp

@@ -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;
 }
 

mesh/pmesh.hpp

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

mesh/pncmesh.cpp

@@ -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();

miniapps/common/fem_extras.cpp

@@ -56,45 +56,9 @@ RT_FESpace::~RT_FESpace()
    delete FEC_;
 }
 
-L2_FESpace::L2_FESpace(Mesh *m, const int p, const int space_dim,
-                       int vdim, int order)
-   : FiniteElementSpace(m, new L2_FECollection(p,space_dim),vdim,order)
-{
-   FEC_ = this->FiniteElementSpace::fec;
-}
-
-L2_FESpace::~L2_FESpace()
-{
-   delete FEC_;
-}
-
-DiscreteInterpolationOperator::~DiscreteInterpolationOperator()
-{}
-
-DiscreteGradOperator::DiscreteGradOperator(FiniteElementSpace *dfes,
-                                           FiniteElementSpace *rfes)
-   : DiscreteInterpolationOperator(dfes, rfes)
-{
-   this->AddDomainInterpolator(new GradientInterpolator);
-}
-
-DiscreteCurlOperator::DiscreteCurlOperator(FiniteElementSpace *dfes,
-                                           FiniteElementSpace *rfes)
-   : DiscreteInterpolationOperator(dfes, rfes)
-{
-   this->AddDomainInterpolator(new CurlInterpolator);
-}
-
-DiscreteDivOperator::DiscreteDivOperator(FiniteElementSpace *dfes,
-                                         FiniteElementSpace *rfes)
-   : DiscreteInterpolationOperator(dfes, rfes)
-{
-   this->AddDomainInterpolator(new DivergenceInterpolator);
-}
-
 void VisualizeMesh(socketstream &sock, const char *vishost, int visport,
                    Mesh &mesh, const char *title,
-                   int x, int y, int w, int h, const char *keys, bool vec)
+                   int x, int y, int w, int h, const char * keys, bool vec)
 {
    bool newly_opened = false;
    int connection_failed;
@@ -129,7 +93,7 @@ void VisualizeMesh(socketstream &sock, const char *vishost, int visport,
 
 void VisualizeField(socketstream &sock, const char *vishost, int visport,
                     GridFunction &gf, const char *title,
-                    int x, int y, int w, int h, const char *keys, bool vec)
+                    int x, int y, int w, int h, const char * keys, bool vec)
 {
    Mesh &mesh = *gf.FESpace()->GetMesh();
 

miniapps/common/fem_extras.hpp

@@ -66,50 +66,6 @@ private:
 };
 
 
-/** The L2_FESpace class is a FiniteElementSpace which automatically
-    allocates and destroys its own FiniteElementCollection, in this
-    case an L2_FECollection object.
-*/
-class L2_FESpace : public FiniteElementSpace
-{
-public:
-   L2_FESpace(Mesh *m, const int p, const int space_dim,
-              int vdim = 1, int order = Ordering::byNODES);
-   ~L2_FESpace();
-private:
-   const FiniteElementCollection *FEC_;
-};
-
-class DiscreteInterpolationOperator : public DiscreteLinearOperator
-{
-public:
-   DiscreteInterpolationOperator(FiniteElementSpace *dfes,
-                                 FiniteElementSpace *rfes)
-      : DiscreteLinearOperator(dfes, rfes) {}
-   virtual ~DiscreteInterpolationOperator();
-};
-
-class DiscreteGradOperator : public DiscreteInterpolationOperator
-{
-public:
-   DiscreteGradOperator(FiniteElementSpace *dfes,
-                        FiniteElementSpace *rfes);
-};
-
-class DiscreteCurlOperator : public DiscreteInterpolationOperator
-{
-public:
-   DiscreteCurlOperator(FiniteElementSpace *dfes,
-                        FiniteElementSpace *rfes);
-};
-
-class DiscreteDivOperator : public DiscreteInterpolationOperator
-{
-public:
-   DiscreteDivOperator(FiniteElementSpace *dfes,
-                       FiniteElementSpace *rfes);
-};
-
 /// Visualize the given mesh object, using a GLVis server on the
 /// specified host and port. Set the visualization window title, and optionally,
 /// its geometry.

miniapps/electromagnetics/CMakeLists.txt

@@ -9,12 +9,6 @@
 # terms of the BSD-3 license. We welcome feedback and contributions, see file
 # CONTRIBUTING.md for details.
 
-add_mfem_miniapp(hertz_ser
-  MAIN hertz_ser.cpp
-  EXTRA_SOURCES hertz_ser_solver.cpp
-  EXTRA_HEADERS hertz_ser_solver.hpp ${MFEM_MINIAPPS_COMMON_HEADERS}
-  LIBRARIES mfem mfem_miniapps_common)
-
 if (MFEM_USE_MPI)
   add_mfem_miniapp(tesla
     MAIN tesla.cpp
@@ -28,12 +22,6 @@ if (MFEM_USE_MPI)
     volta_solver.hpp ${MFEM_MINIAPPS_COMMON_HEADERS}
     LIBRARIES mfem mfem-common)
 
- add_mfem_miniapp(hertz
-   MAIN hertz.cpp
-   EXTRA_SOURCES hertz_solver.cpp
-   EXTRA_HEADERS hertz_solver.hpp ${MFEM_MINIAPPS_COMMON_HEADERS}
-   LIBRARIES mfem mfem_miniapps_common)
-
  add_mfem_miniapp(joule
    MAIN joule.cpp
    joule_solver.cpp

miniapps/electromagnetics/hertz.cpp

@@ -1,672 +0,0 @@
-// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
-// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
-// reserved. See file COPYRIGHT for details.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability see http://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the GNU Lesser General Public License (as published by the Free
-// Software Foundation) version 2.1 dated February 1999.
-//
-//   -----------------------------------------------------------------------
-//   Hertz Miniapp:  Simple Frequency-Domain Electromagnetic Simulation Code
-//   -----------------------------------------------------------------------
-//
-//   Assumes that all sources and boundary conditions oscillate with the same
-//   frequency although not necessarily in phase with one another.  This
-//   assumptions implies that we can factor out the time dependence which we
-//   take to be of the form exp(i omega t).  With these assumptions we can
-//   write the Maxwell equations in the form:
-//
-//   i omega epsilon E = Curl mu^{-1} B - J - sigma E
-//   i omega B         = - Curl E
-//
-//   Which combine to yield:
-//
-//   Curl mu^{-1} Curl E - omega^2 epsilon E + i omega sigma E = - i omega J
-//
-//   We discretize this equation with H(Curl) a.k.a Nedelec basis
-//   functions.  The curl curl operator must be handled with
-//   integration by parts which yields a surface integral:
-//
-//   (W, Curl mu^{-1} Curl E) = (Curl W, mu^{-1} Curl E)
-//               + (W, n x (mu^{-1} Curl E))_{\Gamma}
-//
-//   or
-//
-//   (W, Curl mu^{-1} Curl E) = (Curl W, mu^{-1} Curl E)
-//               - i omega (W, n x H)_{\Gamma}
-//
-//   For plane waves
-//     omega B = - k x E
-//     omega D = k x H, assuming n x k = 0 => n x H = omega epsilon E / |k|
-//
-//   c = omega/|k|
-//
-//   (W, Curl mu^{-1} Curl E) = (Curl W, mu^{-1} Curl E)
-//               - i omega sqrt{epsilon/mu} (W, E)_{\Gamma}
-//
-//
-// Compile with: make hertz
-//
-// Sample runs:
-//
-//   By default the sources and fields are all zero
-//     mpirun -np 4 hertz
-//
-//   Current source in a metal sphere
-//     mpirun -np 4 hertz -m ../../data/ball-nurbs.mesh -rs 2
-//                        -dbcs '-1' -f 3e8 -herm
-//                        -do '-0.3 0.0 0.0 0.3 0.0 0.0 0.1 1 .5 .5'
-//
-//   Current source in a sphere with absorbing boundary conditions
-//     mpirun -np 4 hertz -m ../../data/ball-nurbs.mesh -rs 2
-//                        -abcs '-1' -f 3e8
-//                        -do '-0.3 0.0 0.0 0.3 0.0 0.0 0.1 1 .5 .5'
-//
-//   Current source in a metal sphere with dielectric and conducting materials
-//     mpirun -np 4 hertz -m ../../data/ball-nurbs.mesh -rs 2
-//                        -dbcs '-1' -f 3e8
-//                        -do '-0.3 0.0 0.0 0.3 0.0 0.0 0.1 1 .5 .5'
-//                        -cs '0.0 0.0 -0.5 .2 10'
-//                        -ds '0.0 0.0 0.5 .2 10'
-//
-//   Current source in a metal box
-//     mpirun -np 4 hertz -m ../../data/fichera.mesh -rs 3
-//                        -dbcs '-1' -f 3e8
-//                        -do '-0.5 -0.5 0.0 -0.5 -0.5 1.0 0.1 1 .5 1'
-//
-//   Current source with a mixture of absorbing and reflecting boundaries
-//     mpirun -np 4 hertz -m ../../data/fichera.mesh -rs 3
-//                        -do '-0.5 -0.5 0.0 -0.5 -0.5 1.0 0.1 1 .5 1'
-//                        -dbcs '4 8 19 21' -abcs '5 18' -f 3e8
-//
-
-#include "hertz_solver.hpp"
-#include <fstream>
-#include <iostream>
-
-using namespace std;
-using namespace mfem;
-using namespace mfem::electromagnetics;
-
-// Permittivity Functions
-Coefficient * SetupPermittivityCoefficient();
-
-static Vector pw_eps_(0);     // Piecewise permittivity values
-static Vector ds_params_(0);  // Center, Radius, and Permittivity
-//                               of dielectric sphere
-double dielectric_sphere(const Vector &);
-
-// Permeability Function
-Coefficient * SetupInvPermeabilityCoefficient();
-
-static Vector pw_mu_(0);      // Piecewise permeability values
-static Vector pw_mu_inv_(0);  // Piecewise inverse permeability values
-static Vector ms_params_(0);  // Center, Inner and Outer Radii, and
-//                               Permeability of magnetic shell
-double magnetic_shell(const Vector &);
-double magnetic_shell_inv(const Vector & x) { return 1.0/magnetic_shell(x); }
-
-// Conductivity Functions
-Coefficient * SetupConductivityCoefficient();
-
-static Vector pw_sigma_(0);   // Piecewise conductivity values
-static Vector cs_params_(0);  // Center, Radius, and Conductivity
-//                               of conductive sphere
-double conductive_sphere(const Vector &);
-
-// Impedance
-Coefficient * SetupAdmittanceCoefficient(const Mesh & mesh,
-                                         const Array<int> & abcs);
-
-static Vector pw_eta_(0);      // Piecewise impedance values
-static Vector pw_eta_inv_(0);  // Piecewise inverse impedance values
-
-// Current Density Function
-static Vector do_params_(0);  // Axis Start, Axis End, Rod Radius,
-//                               Total Current of Rod
-void dipole_oscillator(const Vector &x, Vector &j);
-void j_src(const Vector &x, Vector &j) { dipole_oscillator(x, j); }
-
-// Electric Field Boundary Condition: The following function returns zero but
-// any function could be used.
-void e_bc_r(const Vector &x, Vector &E);
-void e_bc_i(const Vector &x, Vector &E);
-
-static double freq_ = 1.0;
-
-// Prints the program's logo to the given output stream
-void display_banner(ostream & os);
-
-int main(int argc, char *argv[])
-{
-   MPI_Session mpi(argc, argv);
-
-   if ( mpi.Root() ) { display_banner(cout); }
-
-   // Parse command-line options.
-   const char *mesh_file = "../../data/ball-nurbs.mesh";
-   int order = 1;
-   int maxit = 100;
-   int serial_ref_levels = 0;
-   int parallel_ref_levels = 0;
-   int sol = 3;
-   int prec = 4;
-   bool herm_conv = false;
-   bool visualization = true;
-   bool visit = true;
-
-   Array<int> abcs;
-   Array<int> dbcs;
-
-   SolverOptions solOpts;
-   solOpts.maxIter = 1000;
-   solOpts.kDim = 50;
-   solOpts.printLvl = 1;
-   solOpts.relTol = 1e-4;
-   solOpts.euLvl = 1;
-
-   OptionsParser args(argc, argv);
-   args.AddOption(&mesh_file, "-m", "--mesh",
-                  "Mesh file to use.");
-   args.AddOption(&order, "-o", "--order",
-                  "Finite element order (polynomial degree).");
-   args.AddOption(&serial_ref_levels, "-rs", "--serial-ref-levels",
-                  "Number of serial refinement levels.");
-   args.AddOption(&parallel_ref_levels, "-rp", "--parallel-ref-levels",
-                  "Number of parallel refinement levels.");
-   args.AddOption(&freq_, "-f", "--frequency",
-                  "Frequency in Hertz (of course...)");
-   args.AddOption(&prec, "-pc", "--precond",
-                  "Preconditioner: 1 - Diagonal Scaling, 2 - ParaSails, "
-                  "3 - Euclid, 4 - AMS");
-   args.AddOption(&sol, "-s", "--solver",
-                  "Solver: 1 - GMRES, 2 - FGMRES, 3 - MINRES"
-#ifdef MFEM_USE_SUPERLU
-                  ", 4 - SuperLU"
-#endif
-#ifdef MFEM_USE_STRUMPACK
-                  ", 5 - STRUMPACK"
-#endif
-                 );
-   args.AddOption(&solOpts.maxIter, "-sol-it", "--solver-iterations",
-                  "Maximum number of solver iterations.");
-   args.AddOption(&solOpts.kDim, "-sol-k-dim", "--solver-krylov-dimension",
-                  "Krylov space dimension for GMRES and FGMRES.");
-   args.AddOption(&solOpts.relTol, "-sol-tol", "--solver-tolerance",
-                  "Relative tolerance for GMRES or FGMRES.");
-   args.AddOption(&solOpts.printLvl, "-sol-prnt-lvl", "--solver-print-level",
-                  "Logging level for solvers.");
-   args.AddOption(&solOpts.euLvl, "-eu-lvl", "--euclid-level",
-                  "Euclid factorization level for ILU(k).");
-   args.AddOption(&pw_eps_, "-pwe", "--piecewise-eps",
-                  "Piecewise values of Permittivity");
-   args.AddOption(&ds_params_, "-ds", "--dielectric-sphere-params",
-                  "Center, Radius, and Permittivity of Dielectric Sphere");
-   args.AddOption(&pw_mu_, "-pwm", "--piecewise-mu",
-                  "Piecewise values of Permeability");
-   args.AddOption(&ms_params_, "-ms", "--magnetic-shell-params",
-                  "Center, Inner Radius, Outer Radius, "
-                  "and Permeability of Magnetic Shell");
-   args.AddOption(&pw_sigma_, "-pws", "--piecewise-sigma",
-                  "Piecewise values of Conductivity");
-   args.AddOption(&cs_params_, "-cs", "--conductive-sphere-params",
-                  "Center, Radius, and Conductivity of Conductive Sphere");
-   args.AddOption(&pw_eta_, "-pwz", "--piecewise-eta",
-                  "Piecewise values of Impedance (one value per abc surface)");
-   args.AddOption(&do_params_, "-do", "--dipole-oscillator-params",
-                  "Axis End Points, Radius, and Amplitude");
-   args.AddOption(&abcs, "-abcs", "--absorbing-bc-surf",
-                  "Absorbing Boundary Condition Surfaces");
-   args.AddOption(&dbcs, "-dbcs", "--dirichlet-bc-surf",
-                  "Dirichlet Boundary Condition Surfaces");
-   args.AddOption(&maxit, "-maxit", "--max-amr-iterations",
-                  "Max number of iterations in the main AMR loop.");
-   args.AddOption(&herm_conv, "-herm", "--hermitian", "-no-herm",
-                  "--no-hermitian", "Use convention for Hermitian operators.");
-   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
-                  "--no-visualization",
-                  "Enable or disable GLVis visualization.");
-   args.AddOption(&visit, "-visit", "--visit", "-no-visit", "--no-visit",
-                  "Enable or disable VisIt visualization.");
-   args.Parse();
-   if (!args.Good())
-   {
-      if (mpi.Root())
-      {
-         args.PrintUsage(cout);
-      }
-      return 1;
-   }
-   if (mpi.Root())
-   {
-      args.PrintOptions(cout);
-   }
-
-   ComplexOperator::Convention conv =
-      herm_conv ? ComplexOperator::HERMITIAN : ComplexOperator::BLOCK_SYMMETRIC;
-
-   // Read the (serial) mesh from the given mesh file on all processors.  We
-   // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
-   // and volume meshes with the same code.
-   Mesh *mesh = new Mesh(mesh_file, 1, 1);
-
-   if (mpi.Root())
-   {
-      cout << "Starting initialization." << endl;
-   }
-
-   // Project a NURBS mesh to a piecewise-quadratic curved mesh
-   if (mesh->NURBSext)
-   {
-      mesh->UniformRefinement();
-      if (serial_ref_levels > 0) { serial_ref_levels--; }
-
-      mesh->SetCurvature(2);
-   }
-
-   // Ensure that quad and hex meshes are treated as non-conforming.
-   mesh->EnsureNCMesh();
-
-   // Refine the serial mesh on all processors to increase the resolution. In
-   // this example we do 'ref_levels' of uniform refinement.
-   for (int l = 0; l < serial_ref_levels; l++)
-   {
-      mesh->UniformRefinement();
-   }
-
-   // Define a parallel mesh by a partitioning of the serial mesh. Refine
-   // this mesh further in parallel to increase the resolution. Once the
-   // parallel mesh is defined, the serial mesh can be deleted.
-   ParMesh pmesh(MPI_COMM_WORLD, *mesh);
-   delete mesh;
-
-   // Refine this mesh in parallel to increase the resolution.
-   int par_ref_levels = parallel_ref_levels;
-   for (int l = 0; l < par_ref_levels; l++)
-   {
-      pmesh.UniformRefinement();
-   }
-
-   // Create a coefficient describing the dielectric permittivity
-   Coefficient * epsCoef = SetupPermittivityCoefficient();
-
-   // Create a coefficient describing the magnetic permeability
-   Coefficient * muInvCoef = SetupInvPermeabilityCoefficient();
-
-   // Create a coefficient describing the electrical conductivity
-   Coefficient * sigmaCoef = SetupConductivityCoefficient();
-
-   // Create a coefficient describing the surface admittance
-   Coefficient * etaInvCoef = SetupAdmittanceCoefficient(pmesh, abcs);
-
-   // Create the Magnetostatic solver
-   HertzSolver Hertz(pmesh, order, freq_,
-                     (HertzSolver::SolverType)sol, solOpts,
-                     (HertzSolver::PrecondType)prec,
-                     conv, *epsCoef, *muInvCoef, sigmaCoef, etaInvCoef,
-                     abcs, dbcs,
-                     e_bc_r, e_bc_i,
-                     (do_params_.Size() > 0 ) ? j_src : NULL, NULL
-                    );
-
-   // Initialize GLVis visualization
-   if (visualization)
-   {
-      Hertz.InitializeGLVis();
-   }
-
-   // Initialize VisIt visualization
-   VisItDataCollection visit_dc("Hertz-AMR-Parallel", &pmesh);
-
-   if ( visit )
-   {
-      Hertz.RegisterVisItFields(visit_dc);
-   }
-   if (mpi.Root()) { cout << "Initialization done." << endl; }
-
-   // The main AMR loop. In each iteration we solve the problem on the current
-   // mesh, visualize the solution, estimate the error on all elements, refine
-   // the worst elements and update all objects to work with the new mesh. We
-   // refine until the maximum number of dofs in the Nedelec finite element
-   // space reaches 10 million.
-   const int max_dofs = 10000000;
-   for (int it = 1; it <= maxit; it++)
-   {
-      if (mpi.Root())
-      {
-         cout << "\nAMR Iteration " << it << endl;
-      }
-
-      // Display the current number of DoFs in each finite element space
-      Hertz.PrintSizes();
-
-      // Assemble all forms
-      Hertz.Assemble();
-
-      // Solve the system and compute any auxiliary fields
-      Hertz.Solve();
-
-      // Determine the current size of the linear system
-      int prob_size = Hertz.GetProblemSize();
-
-      // Write fields to disk for VisIt
-      if ( visit )
-      {
-         Hertz.WriteVisItFields(it);
-      }
-
-      // Send the solution by socket to a GLVis server.
-      if (visualization)
-      {
-         Hertz.DisplayToGLVis();
-      }
-
-      if (mpi.Root())
-      {
-         cout << "AMR iteration " << it << " complete." << endl;
-      }
-
-      // Check stopping criteria
-      if (prob_size > max_dofs)
-      {
-         if (mpi.Root())
-         {
-            cout << "Reached maximum number of dofs, exiting..." << endl;
-         }
-         break;
-      }
-      if ( it == maxit )
-      {
-         break;
-      }
-
-      // Wait for user input. Ask every 10th iteration.
-      char c = 'c';
-      if (mpi.Root() && (it % 10 == 0))
-      {
-         cout << "press (q)uit or (c)ontinue --> " << flush;
-         cin >> c;
-      }
-      MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
-
-      if (c != 'c')
-      {
-         break;
-      }
-
-      // Estimate element errors using the Zienkiewicz-Zhu error estimator.
-      Vector errors(pmesh.GetNE());
-      Hertz.GetErrorEstimates(errors);
-
-      double local_max_err = errors.Max();
-      double global_max_err;
-      MPI_Allreduce(&local_max_err, &global_max_err, 1,
-                    MPI_DOUBLE, MPI_MAX, pmesh.GetComm());
-
-      // Refine the elements whose error is larger than a fraction of the
-      // maximum element error.
-      const double frac = 0.5;
-      double threshold = frac * global_max_err;
-      if (mpi.Root()) { cout << "Refining ..." << endl; }
-      pmesh.RefineByError(errors, threshold);
-
-      // Update the magnetostatic solver to reflect the new state of the mesh.
-      Hertz.Update();
-
-      if (pmesh.Nonconforming() && mpi.WorldSize() > 1 && false)
-      {
-         if (mpi.Root()) { cout << "Rebalancing ..." << endl; }
-         pmesh.Rebalance();
-
-         // Update again after rebalancing
-         Hertz.Update();
-      }
-   }
-
-   // Send the solution by socket to a GLVis server.
-   if (visualization)
-   {
-      Hertz.DisplayAnimationToGLVis();
-   }
-
-   delete epsCoef;
-   delete muInvCoef;
-   delete sigmaCoef;
-   delete etaInvCoef;
-
-   return 0;
-}
-
-// Print the Hertz ascii logo to the given ostream
-void display_banner(ostream & os)
-{
-   os << "     ____  ____              __           " << endl
-      << "    /   / /   / ____________/  |_________ " << endl
-      << "   /   /_/   /_/ __ \\_  __ \\   __\\___   / " << endl
-      << "  /   __    / \\  ___/|  | \\/|  |  /   _/  " << endl
-      << " /___/ /_  /   \\___  >__|   |__| /_____ \\ " << endl
-      << "         \\/        \\/                  \\/ " << endl << flush;
-}
-
-// The Permittivity is a required coefficient which may be defined in
-// various ways so we'll determine the appropriate coefficient type here.
-Coefficient *
-SetupPermittivityCoefficient()
-{
-   Coefficient * coef = NULL;
-
-   if ( ds_params_.Size() > 0 )
-   {
-      coef = new FunctionCoefficient(dielectric_sphere);
-   }
-   else if ( pw_eps_.Size() > 0 )
-   {
-      coef = new PWConstCoefficient(pw_eps_);
-   }
-   else
-   {
-      coef = new ConstantCoefficient(epsilon0_);
-   }
-
-   return coef;
-}
-
-// The Permeability is a required coefficient which may be defined in
-// various ways so we'll determine the appropriate coefficient type here.
-Coefficient *
-SetupInvPermeabilityCoefficient()
-{
-   Coefficient * coef = NULL;
-
-   if ( ms_params_.Size() > 0 )
-   {
-      coef = new FunctionCoefficient(magnetic_shell_inv);
-   }
-   else if ( pw_mu_.Size() > 0 )
-   {
-      pw_mu_inv_.SetSize(pw_mu_.Size());
-      for (int i = 0; i < pw_mu_.Size(); i++)
-      {
-         MFEM_ASSERT( pw_mu_[i] > 0.0, "permeability values must be positive" );
-         pw_mu_inv_[i] = 1.0/pw_mu_[i];
-      }
-      coef = new PWConstCoefficient(pw_mu_inv_);
-   }
-   else
-   {
-      coef = new ConstantCoefficient(1.0/mu0_);
-   }
-
-   return coef;
-}
-
-// The Conductivity is an optional coefficient which may be defined in
-// various ways so we'll determine the appropriate coefficient type here.
-Coefficient *
-SetupConductivityCoefficient()
-{
-   Coefficient * coef = NULL;
-
-   if ( cs_params_.Size() > 0 )
-   {
-      coef = new FunctionCoefficient(conductive_sphere);
-   }
-   else if ( pw_sigma_.Size() > 0 )
-   {
-      coef = new PWConstCoefficient(pw_sigma_);
-   }
-
-   return coef;
-}
-
-// The Admittance is an optional coefficient defined on boundary surfaces which
-// can be used in conjunction with absorbing boundary conditions.
-Coefficient *
-SetupAdmittanceCoefficient(const Mesh & mesh, const Array<int> & abcs)
-{
-   Coefficient * coef = NULL;
-
-   if ( pw_eta_.Size() > 0 )
-   {
-      MFEM_VERIFY(pw_eta_.Size() == abcs.Size(),
-                  "Each impedance value must be associated with exactly one "
-                  "absorbing boundary surface.");
-
-      pw_eta_inv_.SetSize(mesh.bdr_attributes.Size());
-
-      if ( abcs[0] == -1 )
-      {
-         pw_eta_inv_ = 1.0 / pw_eta_[0];
-      }
-      else
-      {
-         pw_eta_inv_ = 0.0;
-
-         for (int i=0; i<pw_eta_.Size(); i++)
-         {
-            pw_eta_inv_[abcs[i]-1] = 1.0 / pw_eta_[i];
-         }
-      }
-      coef = new PWConstCoefficient(pw_eta_inv_);
-   }
-
-   return coef;
-}
-
-// A sphere with constant permittivity.  The sphere has a radius,
-// center, and permittivity specified on the command line and stored
-// in ds_params_.
-double dielectric_sphere(const Vector &x)
-{
-   double r2 = 0.0;
-
-   for (int i=0; i<x.Size(); i++)
-   {
-      r2 += (x(i)-ds_params_(i))*(x(i)-ds_params_(i));
-   }
-
-   if ( sqrt(r2) <= ds_params_(x.Size()) )
-   {
-      return ds_params_(x.Size()+1) * epsilon0_;
-   }
-   return epsilon0_;
-}
-
-// A spherical shell with constant permeability.  The sphere has inner
-// and outer radii, center, and relative permeability specified on the
-// command line and stored in ms_params_.
-double magnetic_shell(const Vector &x)
-{
-   double r2 = 0.0;
-
-   for (int i = 0; i < x.Size(); i++)
-   {
-      r2 += (x(i) - ms_params_(i))*(x(i) - ms_params_(i));
-   }
-
-   if ( sqrt(r2) >= ms_params_(x.Size()) &&
-        sqrt(r2) <= ms_params_(x.Size()+1) )
-   {
-      return mu0_*ms_params_(x.Size()+2);
-   }
-   return mu0_;
-}
-
-// A sphere with constant conductivity.  The sphere has a radius,
-// center, and conductivity specified on the command line and stored
-// in ls_params_.
-double conductive_sphere(const Vector &x)
-{
-   double r2 = 0.0;
-
-   for (int i=0; i<x.Size(); i++)
-   {
-      r2 += (x(i)-cs_params_(i))*(x(i)-cs_params_(i));
-   }
-
-   if ( sqrt(r2) <= cs_params_(x.Size()) )
-   {
-      return cs_params_(x.Size()+1);
-   }
-   return 0.0;
-}
-
-// A cylindrical rod of current density.  The rod has two axis end
-// points, a radus, a current amplitude in Amperes.  All of these
-// parameters are stored in do_params_.
-void dipole_oscillator(const Vector &x, Vector &j)
-{
-   MFEM_ASSERT(x.Size() == 3, "current source requires 3D space.");
-
-   j.SetSize(x.Size());
-   j = 0.0;
-
-   Vector  v(x.Size());  // Normalized Axis vector
-   Vector xu(x.Size());  // x vector relative to the axis end-point
-
-   xu = x;
-
-   for (int i=0; i<x.Size(); i++)
-   {
-      xu[i] -= do_params_[i];
-      v[i]   = do_params_[x.Size()+i] - do_params_[i];
-   }
-
-   double h = v.Norml2();
-
-   if ( h == 0.0 )
-   {
-      return;
-   }
-   v /= h;
-
-   double r = do_params_[2*x.Size()+0];
-   double a = do_params_[2*x.Size()+1];
-
-   double xv = xu * v;
-
-   // Compute perpendicular vector from axis to x
-   xu.Add(-xv, v);
-
-   double xp = xu.Norml2();
-
-   if ( xv >= 0.0 && xv <= h && xp <= r )
-   {
-      j.Add(a, v);
-   }
-}
-
-void e_bc_r(const Vector &x, Vector &E)
-{
-   E.SetSize(3);
-   E = 0.0;
-}
-
-void e_bc_i(const Vector &x, Vector &E)
-{
-   E.SetSize(3);
-   E = 0.0;
-}

miniapps/electromagnetics/hertz_ser.cpp

@@ -1,638 +0,0 @@
-// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
-// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
-// reserved. See file COPYRIGHT for details.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability see http://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the GNU Lesser General Public License (as published by the Free
-// Software Foundation) version 2.1 dated February 1999.
-//
-//   -----------------------------------------------------------------------
-//   Hertz Miniapp:  Simple Frequency-Domain Electromagnetic Simulation Code
-//   -----------------------------------------------------------------------
-//
-//   Assumes that all sources and boundary conditions oscillate with the same
-//   frequency although not necessarily in phase with one another.  This
-//   assumptions implies that we can factor out the time dependence which we
-//   take to be of the form exp(i omega t).  With these assumptions we can
-//   write the Maxwell equations in the form:
-//
-//   i omega epsilon E = Curl mu^{-1} B - J - sigma E
-//   i omega B         = - Curl E
-//
-//   Which combine to yield:
-//
-//   Curl mu^{-1} Curl E - omega^2 epsilon E + i omega sigma E = - i omega J
-//
-//   We discretize this equation with H(Curl) a.k.a Nedelec basis
-//   functions.  The curl curl operator must be handled with
-//   integration by parts which yields a surface integral:
-//
-//   (W, Curl mu^{-1} Curl E) = (Curl W, mu^{-1} Curl E)
-//               + (W, n x (mu^{-1} Curl E))_{\Gamma}
-//
-//   or
-//
-//   (W, Curl mu^{-1} Curl E) = (Curl W, mu^{-1} Curl E)
-//               - i omega (W, n x H)_{\Gamma}
-//
-//   For plane waves
-//     omega B = - k x E
-//     omega D = k x H, assuming n x k = 0 => n x H = omega epsilon E / |k|
-//
-//   c = omega/|k|
-//
-//   (W, Curl mu^{-1} Curl E) = (Curl W, mu^{-1} Curl E)
-//               - i omega sqrt{epsilon/mu} (W, E)_{\Gamma}
-//
-//
-// Compile with: make hertz
-//
-// Sample runs:
-//
-//   By default the sources and fields are all zero
-//     hertz_ser
-//
-//   Current source in a sphere with absorbing boundary conditions
-//     hertz_ser -m ../../data/ball-nurbs.mesh -rs 2
-//               -abcs '-1'
-//               -do '-0.3 0.0 0.0 0.3 0.0 0.0 0.1 1 .5 .5'
-//
-//   Current source in a metal sphere with dielectric and conducting materials
-//     hertz_ser -m ../../data/ball-nurbs.mesh -rs 2
-//               -dbcs '-1'
-//               -do '-0.3 0.0 0.0 0.3 0.0 0.0 0.1 1 .5 .5'
-//               -cs '0.0 0.0 -0.5 .2 10'
-//               -ds '0.0 0.0 0.5 .2 10'
-//
-//   Current source in a metal box
-//     hertz_ser -m ../../data/fichera.mesh -rs 3
-//               -dbcs '-1'
-//               -do '-0.5 -0.5 0.0 -0.5 -0.5 1.0 0.1 1 .5 1'
-//
-//   Current source with a mixture of absorbing and reflecting boundaries
-//     hertz_ser -m ../../data/fichera.mesh -rs 3
-//               -do '-0.5 -0.5 0.0 -0.5 -0.5 1.0 0.1 1 .5 1'
-//               -dbcs '4 8 19 21' -abcs '5 18'
-//
-
-#include "hertz_ser_solver.hpp"
-#include <fstream>
-#include <iostream>
-
-using namespace std;
-using namespace mfem;
-using namespace mfem::electromagnetics;
-
-// Permittivity Functions
-Coefficient * SetupPermittivityCoefficient();
-
-static Vector pw_eps_(0);     // Piecewise permittivity values
-static Vector ds_params_(0);  // Center, Radius, and Permittivity
-//                               of dielectric sphere
-double dielectric_sphere(const Vector &);
-
-// Permeability Function
-Coefficient * SetupInvPermeabilityCoefficient();
-
-static Vector pw_mu_(0);      // Piecewise permeability values
-static Vector pw_mu_inv_(0);  // Piecewise inverse permeability values
-static Vector ms_params_(0);  // Center, Inner and Outer Radii, and
-//                               Permeability of magnetic shell
-double magnetic_shell(const Vector &);
-double magnetic_shell_inv(const Vector & x) { return 1.0/magnetic_shell(x); }
-
-// Conductivity Functions
-Coefficient * SetupConductivityCoefficient();
-
-static Vector pw_sigma_(0);   // Piecewise conductivity values
-static Vector cs_params_(0);  // Center, Radius, and Conductivity
-//                               of conductive sphere
-double conductive_sphere(const Vector &);
-
-// Impedance
-Coefficient * SetupAdmittanceCoefficient(const Mesh & mesh,
-                                         const Array<int> & abcs);
-
-static Vector pw_eta_(0);      // Piecewise impedance values
-static Vector pw_eta_inv_(0);  // Piecewise inverse impedance values
-
-// Current Density Function
-static Vector do_params_(0);  // Axis Start, Axis End, Rod Radius,
-//                               Total Current of Rod
-void dipole_oscillator(const Vector &x, Vector &j);
-void j_src(const Vector &x, Vector &j) { dipole_oscillator(x, j); }
-
-// Electric Field Boundary Condition: The following function returns zero but
-// any function could be used.
-void e_bc_r(const Vector &x, Vector &E);
-void e_bc_i(const Vector &x, Vector &E);
-
-static double freq_ = 1.0;
-
-// Prints the program's logo to the given output stream
-void display_banner(ostream & os);
-
-int main(int argc, char *argv[])
-{
-   display_banner(cout);
-
-   // Parse command-line options.
-   const char *mesh_file = "../../data/ball-nurbs.mesh";
-   int order = 1;
-   int maxit = 100;
-   int serial_ref_levels = 0;
-   int sol = 1;
-   bool herm_conv = false;
-   bool visualization = true;
-   bool visit = true;
-
-   Array<int> abcs;
-   Array<int> dbcs;
-
-   OptionsParser args(argc, argv);
-   args.AddOption(&mesh_file, "-m", "--mesh",
-                  "Mesh file to use.");
-   args.AddOption(&order, "-o", "--order",
-                  "Finite element order (polynomial degree).");
-   args.AddOption(&serial_ref_levels, "-rs", "--serial-ref-levels",
-                  "Number of serial refinement levels.");
-   args.AddOption(&freq_, "-f", "--frequency",
-                  "Frequency in Hertz (of course...)");
-   args.AddOption(&sol, "-s", "--solver",
-                  "Solver: 1 - GMRES");
-   args.AddOption(&pw_eps_, "-pwe", "--piecewise-eps",
-                  "Piecewise values of Permittivity");
-   args.AddOption(&ds_params_, "-ds", "--dielectric-sphere-params",
-                  "Center, Radius, and Permittivity of Dielectric Sphere");
-   args.AddOption(&pw_mu_, "-pwm", "--piecewise-mu",
-                  "Piecewise values of Permeability");
-   args.AddOption(&ms_params_, "-ms", "--magnetic-shell-params",
-                  "Center, Inner Radius, Outer Radius, "
-                  "and Permeability of Magnetic Shell");
-   args.AddOption(&pw_sigma_, "-pws", "--piecewise-sigma",
-                  "Piecewise values of Conductivity");
-   args.AddOption(&cs_params_, "-cs", "--conductive-sphere-params",
-                  "Center, Radius, and Conductivity of Conductive Sphere");
-   args.AddOption(&pw_eta_, "-pwz", "--piecewise-eta",
-                  "Piecewise values of Impedance (one value per abc surface)");
-   args.AddOption(&do_params_, "-do", "--dipole-oscillator-params",
-                  "Axis End Points, Radius, and Amplitude");
-   args.AddOption(&abcs, "-abcs", "--absorbing-bc-surf",
-                  "Absorbing Boundary Condition Surfaces");
-   args.AddOption(&dbcs, "-dbcs", "--dirichlet-bc-surf",
-                  "Dirichlet Boundary Condition Surfaces");
-   /*
-   args.AddOption(&dbcv, "-dbcv", "--dirichlet-bc-vals",
-                  "Dirichlet Boundary Condition Values");
-   args.AddOption(&dbcg, "-dbcg", "--dirichlet-bc-gradient",
-                  "-no-dbcg", "--no-dirichlet-bc-gradient",
-                  "Dirichlet Boundary Condition Gradient (phi = -z)");
-   args.AddOption(&nbcs, "-nbcs", "--neumann-bc-surf",
-                  "Neumann Boundary Condition Surfaces");
-   args.AddOption(&nbcv, "-nbcv", "--neumann-bc-vals",
-                  "Neumann Boundary Condition Values");
-   args.AddOption(&kbcs, "-kbcs", "--surface-current-bc",
-                  "Surfaces for the Surface Current (K) Boundary Condition");
-   args.AddOption(&vbcs, "-vbcs", "--voltage-bc-surf",
-                  "Voltage Boundary Condition Surfaces (to drive K)");
-   args.AddOption(&vbcv, "-vbcv", "--voltage-bc-vals",
-                  "Voltage Boundary Condition Values (to drive K)");
-   */
-   args.AddOption(&maxit, "-maxit", "--max-amr-iterations",
-                  "Max number of iterations in the main AMR loop.");
-   args.AddOption(&herm_conv, "-herm", "--hermitian", "-no-herm",
-                  "--no-hermitian", "Use convention for Hermitian operators.");
-   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
-                  "--no-visualization",
-                  "Enable or disable GLVis visualization.");
-   args.AddOption(&visit, "-visit", "--visit", "-no-visit", "--no-visit",
-                  "Enable or disable VisIt visualization.");
-   args.Parse();
-   if (!args.Good())
-   {
-      args.PrintUsage(cout);
-      return 1;
-   }
-   args.PrintOptions(cout);
-
-   ComplexOperator::Convention conv =
-      herm_conv ? ComplexOperator::HERMITIAN : ComplexOperator::BLOCK_SYMMETRIC;
-
-   // Read the (serial) mesh from the given mesh file on all processors.  We
-   // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
-   // and volume meshes with the same code.
-   Mesh *mesh = new Mesh(mesh_file, 1, 1);
-
-   cout << "Starting initialization." << endl;
-
-   // Project a NURBS mesh to a piecewise-quadratic curved mesh
-   if (mesh->NURBSext)
-   {
-      mesh->UniformRefinement();
-      if (serial_ref_levels > 0) { serial_ref_levels--; }
-
-      mesh->SetCurvature(2);
-   }
-
-   // Ensure that quad and hex meshes are treated as non-conforming.
-   mesh->EnsureNCMesh();
-
-   // Refine the serial mesh on all processors to increase the resolution. In
-   // this example we do 'ref_levels' of uniform refinement.
-   for (int l = 0; l < serial_ref_levels; l++)
-   {
-      mesh->UniformRefinement();
-   }
-
-   // If values for Voltage BCs were not set issue a warning and exit
-   /*
-   if ( ( vbcs.Size() > 0 && kbcs.Size() == 0 ) ||
-        ( kbcs.Size() > 0 && vbcs.Size() == 0 ) ||
-        ( vbcv.Size() < vbcs.Size() ) )
-   {
-      if ( mpi.Root() )
-      {
-         cout << "The surface current (K) boundary condition requires "
-              << "surface current boundary condition surfaces (with -kbcs), "
-              << "voltage boundary condition surface (with -vbcs), "
-              << "and voltage boundary condition values (with -vbcv)."
-              << endl;
-      }
-      return 3;
-   }
-   */
-   // Create a coefficient describing the dielectric permittivity
-   Coefficient * epsCoef = SetupPermittivityCoefficient();
-
-   // Create a coefficient describing the magnetic permeability
-   Coefficient * muInvCoef = SetupInvPermeabilityCoefficient();
-
-   // Create a coefficient describing the electrical conductivity
-   Coefficient * sigmaCoef = SetupConductivityCoefficient();
-
-   // Create a coefficient describing the surface admittance
-   Coefficient * etaInvCoef = SetupAdmittanceCoefficient(*mesh, abcs);
-
-   // Create the Magnetostatic solver
-   HertzSolver Hertz(*mesh, order, freq_, (HertzSolver::SolverType)sol,
-                     conv, *epsCoef, *muInvCoef, sigmaCoef, etaInvCoef,
-                     abcs, dbcs,
-                     e_bc_r, e_bc_i,
-                     (do_params_.Size() > 0 ) ? j_src : NULL, NULL
-                    );
-
-   //(b_uniform_.Size() > 0 ) ? a_bc_uniform  : NULL,
-   //(cr_params_.Size() > 0 ) ? current_ring  : NULL,
-   //(bm_params_.Size() > 0 ) ? bar_magnet    :
-   //(ha_params_.Size() > 0 ) ? halbach_array : NULL);
-
-   // Initialize GLVis visualization
-   if (visualization)
-   {
-      Hertz.InitializeGLVis();
-   }
-
-   // Initialize VisIt visualization
-   VisItDataCollection visit_dc("Hertz-AMR-Serial", mesh);
-
-   if ( visit )
-   {
-      Hertz.RegisterVisItFields(visit_dc);
-   }
-
-   // ZienkiewiczZhuEstimator estimator(*integ, x, flux_fespace);
-   ZienkiewiczZhuEstimator estimator(*Hertz.GetErrorInteg(),
-                                     *Hertz.GetErrorField(),
-                                     Hertz.GetErrorFluxFES());
-   ThresholdRefiner refiner(estimator);
-   refiner.SetTotalErrorFraction(0.5);
-
-   cout << "Initialization done." << endl;
-
-   // The main AMR loop. In each iteration we solve the problem on the current
-   // mesh, visualize the solution, estimate the error on all elements, refine
-   // the worst elements and update all objects to work with the new mesh. We
-   // refine until the maximum number of dofs in the Nedelec finite element
-   // space reaches 10 million.
-   const int max_dofs = 10000000;
-   for (int it = 1; it <= maxit; it++)
-   {
-      cout << "\nAMR Iteration " << it << endl;
-
-      // Display the current number of DoFs in each finite element space
-      Hertz.PrintSizes();
-
-      // Assemble all forms
-      Hertz.Assemble();
-
-      // Solve the system and compute any auxiliary fields
-      Hertz.Solve();
-
-      // Determine the current size of the linear system
-      int prob_size = Hertz.GetProblemSize();
-
-      // Write fields to disk for VisIt
-      if ( visit )
-      {
-         Hertz.WriteVisItFields(it);
-      }
-
-      // Send the solution by socket to a GLVis server.
-      if (visualization)
-      {
-         Hertz.DisplayToGLVis();
-      }
-
-      cout << "AMR iteration " << it << " complete." << endl;
-
-      // Check stopping criteria
-      if (prob_size > max_dofs)
-      {
-         cout << "Reached maximum number of dofs, exiting..." << endl;
-         break;
-      }
-      if ( it == maxit )
-      {
-         break;
-      }
-
-      // Wait for user input. Ask every 10th iteration.
-      char c = 'c';
-      if (it % 10 == 0)
-      {
-         cout << "press (q)uit or (c)ontinue --> " << flush;
-         cin >> c;
-      }
-
-      if (c != 'c')
-      {
-         break;
-      }
-
-      // Estimate element errors using the Zienkiewicz-Zhu error estimator.
-      refiner.Apply(*mesh);
-      if (refiner.Stop())
-      {
-         cout << "Stopping criterion satisfied. Stop." << endl;
-         break;
-      }
-      /*
-      Vector errors(mesh.GetNE());
-      Hertz.GetErrorEstimates(errors);
-
-      double global_max_err = errors.Max();
-
-      // Refine the elements whose error is larger than a fraction of the
-      // maximum element error.
-      const double frac = 0.5;
-      double threshold = frac * global_max_err;
-      cout << "Refining ..." << endl;
-      mesh.RefineByError(errors, threshold);
-      */
-
-      // Update the magnetostatic solver to reflect the new state of the mesh.
-      Hertz.Update();
-
-   }
-
-   delete epsCoef;
-   delete muInvCoef;
-   delete sigmaCoef;
-
-   return 0;
-}
-
-// Print the Hertz ascii logo to the given ostream
-void display_banner(ostream & os)
-{
-   os << "     ____  ____              __           " << endl
-      << "    /   / /   / ____________/  |_________ " << endl
-      << "   /   /_/   /_/ __ \\_  __ \\   __\\___   / " << endl
-      << "  /   __    / \\  ___/|  | \\/|  |  /   _/  " << endl
-      << " /___/ /_  /   \\___  >__|   |__| /_____ \\ " << endl
-      << "         \\/        \\/                  \\/ " << endl << flush;
-}
-
-// The Permittivity is a required coefficient which may be defined in
-// various ways so we'll determine the appropriate coefficient type here.
-Coefficient *
-SetupPermittivityCoefficient()
-{
-   Coefficient * coef = NULL;
-
-   if ( ds_params_.Size() > 0 )
-   {
-      coef = new FunctionCoefficient(dielectric_sphere);
-   }
-   else if ( pw_eps_.Size() > 0 )
-   {
-      coef = new PWConstCoefficient(pw_eps_);
-   }
-   else
-   {
-      coef = new ConstantCoefficient(epsilon0_);
-   }
-
-   return coef;
-}
-
-// The Permeability is a required coefficient which may be defined in
-// various ways so we'll determine the appropriate coefficient type here.
-Coefficient *
-SetupInvPermeabilityCoefficient()
-{
-   Coefficient * coef = NULL;
-
-   if ( ms_params_.Size() > 0 )
-   {
-      coef = new FunctionCoefficient(magnetic_shell_inv);
-   }
-   else if ( pw_mu_.Size() > 0 )
-   {
-      pw_mu_inv_.SetSize(pw_mu_.Size());
-      for (int i = 0; i < pw_mu_.Size(); i++)
-      {
-         MFEM_ASSERT( pw_mu_[i] > 0.0, "permeability values must be positive" );
-         pw_mu_inv_[i] = 1.0/pw_mu_[i];
-      }
-      coef = new PWConstCoefficient(pw_mu_inv_);
-   }
-   else
-   {
-      coef = new ConstantCoefficient(1.0/mu0_);
-   }
-
-   return coef;
-}
-
-// The Conductivity is an optional coefficient which may be defined in
-// various ways so we'll determine the appropriate coefficient type here.
-Coefficient *
-SetupConductivityCoefficient()
-{
-   Coefficient * coef = NULL;
-
-   if ( cs_params_.Size() > 0 )
-   {
-      coef = new FunctionCoefficient(conductive_sphere);
-   }
-   else if ( pw_sigma_.Size() > 0 )
-   {
-      coef = new PWConstCoefficient(pw_sigma_);
-   }
-
-   return coef;
-}
-
-// The Admittance is an optional coefficient defined on boundary surfaces which
-// can be used in conjunction with absorbing boundary conditions.
-Coefficient *
-SetupAdmittanceCoefficient(const Mesh & mesh, const Array<int> & abcs)
-{
-   Coefficient * coef = NULL;
-
-   if ( pw_eta_.Size() > 0 )
-   {
-      MFEM_VERIFY(pw_eta_.Size() == abcs.Size(),
-                  "Each impedance value must be associated with exactly one "
-                  "absorbing boundary surface.");
-
-      pw_eta_inv_.SetSize(mesh.bdr_attributes.Size());
-
-      if ( abcs[0] == -1 )
-      {
-         pw_eta_inv_ = 1.0 / pw_eta_[0];
-      }
-      else
-      {
-         pw_eta_inv_ = 0.0;
-
-         for (int i=0; i<pw_eta_.Size(); i++)
-         {
-            pw_eta_inv_[abcs[i]-1] = 1.0 / pw_eta_[i];
-         }
-      }
-      coef = new PWConstCoefficient(pw_eta_inv_);
-   }
-
-   return coef;
-}
-
-// A sphere with constant permittivity.  The sphere has a radius,
-// center, and permittivity specified on the command line and stored
-// in ds_params_.
-double dielectric_sphere(const Vector &x)
-{
-   double r2 = 0.0;
-
-   for (int i=0; i<x.Size(); i++)
-   {
-      r2 += (x(i)-ds_params_(i))*(x(i)-ds_params_(i));
-   }
-
-   if ( sqrt(r2) <= ds_params_(x.Size()) )
-   {
-      return ds_params_(x.Size()+1) * epsilon0_;
-   }
-   return epsilon0_;
-}
-
-// A spherical shell with constant permeability.  The sphere has inner
-// and outer radii, center, and relative permeability specified on the
-// command line and stored in ms_params_.
-double magnetic_shell(const Vector &x)
-{
-   double r2 = 0.0;
-
-   for (int i = 0; i < x.Size(); i++)
-   {
-      r2 += (x(i) - ms_params_(i))*(x(i) - ms_params_(i));
-   }
-
-   if ( sqrt(r2) >= ms_params_(x.Size()) &&
-        sqrt(r2) <= ms_params_(x.Size()+1) )
-   {
-      return mu0_*ms_params_(x.Size()+2);
-   }
-   return mu0_;
-}
-
-// A sphere with constant conductivity.  The sphere has a radius,
-// center, and conductivity specified on the command line and stored
-// in ls_params_.
-double conductive_sphere(const Vector &x)
-{
-   double r2 = 0.0;
-
-   for (int i=0; i<x.Size(); i++)
-   {
-      r2 += (x(i)-cs_params_(i))*(x(i)-cs_params_(i));
-   }
-
-   if ( sqrt(r2) <= cs_params_(x.Size()) )
-   {
-      return cs_params_(x.Size()+1);
-   }
-   return 0.0;
-}
-
-// A cylindrical rod of current density.  The rod has two axis end
-// points, a radus, a current amplitude in Amperes.  All of these
-// parameters are stored in do_params_.
-void dipole_oscillator(const Vector &x, Vector &j)
-{
-   MFEM_ASSERT(x.Size() == 3, "current source requires 3D space.");
-
-   j.SetSize(x.Size());
-   j = 0.0;
-
-   Vector  v(x.Size());  // Normalized Axis vector
-   Vector xu(x.Size());  // x vector relative to the axis end-point
-
-   xu = x;
-
-   for (int i=0; i<x.Size(); i++)
-   {
-      xu[i] -= do_params_[i];
-      v[i]   = do_params_[x.Size()+i] - do_params_[i];
-   }
-
-   double h = v.Norml2();
-
-   if ( h == 0.0 )
-   {
-      return;
-   }
-   v /= h;
-
-   double r = do_params_[2*x.Size()+0];
-   double a = do_params_[2*x.Size()+1];
-
-   double xv = xu * v;
-
-   // Compute perpendicular vector from axis to x
-   xu.Add(-xv, v);
-
-   double xp = xu.Norml2();
-
-   if ( xv >= 0.0 && xv <= h && xp <= r )
-   {
-      j.Add(a, v);
-   }
-}
-
-void e_bc_r(const Vector &x, Vector &E)
-{
-   E.SetSize(3);
-   E = 0.0;
-
-}
-
-void e_bc_i(const Vector &x, Vector &E)
-{
-   E.SetSize(3);
-   E = 0.0;
-}

miniapps/electromagnetics/hertz_ser_solver.hpp

@@ -1,248 +0,0 @@
-// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
-// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
-// reserved. See file COPYRIGHT for details.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability see http://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the GNU Lesser General Public License (as published by the Free
-// Software Foundation) version 2.1 dated February 1999.
-
-#ifndef MFEM_HERTZ_SOLVER
-#define MFEM_HERTZ_SOLVER
-
-#include "../common/fem_extras.hpp"
-
-#include <string>
-#include <map>
-
-namespace mfem
-{
-
-using common::H1_FESpace;
-using common::ND_FESpace;
-using common::RT_FESpace;
-using common::DiscreteGradOperator;
-using common::DiscreteCurlOperator;
-//using common::DivergenceFreeProjector;
-
-namespace electromagnetics
-{
-
-// Physical Constants
-// Permittivity of Free Space (units F/m)
-static const double epsilon0_ = 8.8541878176e-12;
-// Permeability of Free Space (units H/m)
-static const double mu0_ = 4.0e-7*M_PI;
-
-//class SurfaceCurrent;
-class HertzSolver
-{
-public:
-
-   enum SolverType
-   {
-      INVALID   = -1,
-      GMRES     =  1,
-      FGMRES    =  2,
-      MINRES    =  3,
-      SUPERLU   =  4,
-      STRUMPACK =  5
-   };
-
-   HertzSolver(Mesh & mesh, int order, double freq,
-               HertzSolver::SolverType s,
-               ComplexOperator::Convention conv,
-               Coefficient & epsCoef,
-               Coefficient & muInvCoef,
-               Coefficient * sigmaCoef,
-               Coefficient * etaInvCoef,
-               Array<int> & abcs,
-               Array<int> & dbcs,
-               void   (*e_r_bc )(const Vector&, Vector&),
-               void   (*e_i_bc )(const Vector&, Vector&),
-               void   (*j_r_src)(const Vector&, Vector&),
-               void   (*j_i_src)(const Vector&, Vector&));
-   ~HertzSolver();
-
-   int GetProblemSize();
-
-   void PrintSizes();
-
-   void Assemble();
-
-   void Update();
-
-   void Solve();
-
-   //void GetErrorEstimates(Vector & errors);
-
-   void RegisterVisItFields(VisItDataCollection & visit_dc);
-
-   void WriteVisItFields(int it = 0);
-
-   void InitializeGLVis();
-
-   void DisplayToGLVis();
-
-   BilinearFormIntegrator * GetErrorInteg() { return err_integ_; }
-
-   GridFunction * GetErrorField() { return &e_->real(); }
-
-   FiniteElementSpace * GetErrorFluxFES()
-   {
-      return new FiniteElementSpace(mesh_, err_flux_fec_);
-   }
-   // const ParGridFunction & GetVectorPotential() { return *a_; }
-
-private:
-
-   int order_;
-   int logging_;
-
-   SolverType sol_;
-
-   ComplexOperator::Convention conv_;
-
-   bool ownsEtaInv_;
-
-   double freq_;
-
-   Mesh * mesh_;
-
-   // H1_ParFESpace * H1FESpace_;
-   ND_FESpace * HCurlFESpace_;
-   // RT_ParFESpace * HDivFESpace_;
-
-   Array<int> blockTrueOffsets_;
-
-   // ParSesquilinearForm * a0_;
-   SesquilinearForm * a1_;
-   BilinearForm * b1_;
-
-   // ParGridFunction * e_r_;  // Real part of electric field (HCurl)
-   // ParGridFunction * e_i_;  // Imaginary part of electric field (HCurl)
-   ComplexGridFunction * e_;  // Complex electric field (HCurl)
-   ComplexGridFunction * j_;  // Complex current density (HCurl)
-   // ParGridFunction * j_i_;  // Imaginary part of current density (HCurl)
-
-   ComplexLinearForm   * jd_; // Dual of complex current density (HCurl)
-   // ParLinearForm   * jd_r_; // Dual of real part of current density (HCurl)
-   // ParLinearForm   * jd_i_; // Dual of imaginary part of current density (HCurl)
-
-   /*
-     ParBilinearForm * curlMuInvCurl_;
-     ParBilinearForm * hCurlMass_;
-     ParMixedBilinearForm * hDivHCurlMuInv_;
-     ParMixedBilinearForm * weakCurlMuInv_;
-    */
-   // ParDiscreteGradOperator * grad_;
-   // ParDiscreteCurlOperator * curl_;
-   /*
-    ParGridFunction * a_;  // Vector Potential (HCurl)
-    ParGridFunction * b_;  // Magnetic Flux (HDiv)
-    ParGridFunction * h_;  // Magnetic Field (HCurl)
-    ParGridFunction * jr_; // Raw Volumetric Current Density (HCurl)
-    ParGridFunction * j_;  // Volumetric Current Density (HCurl)
-    ParGridFunction * k_;  // Surface Current Density (HCurl)
-    ParGridFunction * m_;  // Magnetization (HDiv)
-    ParGridFunction * bd_; // Dual of B (HCurl)
-    ParGridFunction * jd_; // Dual of J, the rhs vector (HCurl)
-   */
-   // DivergenceFreeProjector * DivFreeProj_;
-   // SurfaceCurrent          * SurfCur_;
-
-   Coefficient       * epsCoef_;   // Dielectric Material Coefficient
-   Coefficient       * muInvCoef_; // Dia/Paramagnetic Material Coefficient
-   Coefficient       * sigmaCoef_; // Electrical Conductivity Coefficient
-   Coefficient       * etaInvCoef_; // Admittance Coefficient
-
-   Coefficient * omegaCoef_;     // omega expressed as a Coefficient
-   Coefficient * negOmegaCoef_;  // -omega expressed as a Coefficient
-   Coefficient * omega2Coef_;    // omega^2 expressed as a Coefficient
-   Coefficient * negOmega2Coef_; // -omega^2 expressed as a Coefficient
-   Coefficient * massCoef_;      // -omega^2 epsilon
-   Coefficient * posMassCoef_;   // omega^2 epsilon
-   Coefficient * lossCoef_;      // -omega sigma
-   // Coefficient * gainCoef_;    // omega sigma
-   Coefficient * abcCoef_;       // -omega eta^{-1}
-
-   // VectorCoefficient * aBCCoef_;   // Vector Potential BC Function
-   VectorCoefficient * jrCoef_;     // Volume Current Density Function
-   VectorCoefficient * jiCoef_;     // Volume Current Density Function
-   VectorCoefficient * erCoef_;     // Electric Field Boundary Condition
-   VectorCoefficient * eiCoef_;     // Electric Field Boundary Condition
-   // VectorCoefficient * mCoef_;     // Magnetization Vector Function
-
-   BilinearFormIntegrator * err_integ_;
-   FiniteElementCollection * err_flux_fec_;
-
-   // void   (*a_bc_ )(const Vector&, Vector&);
-   void   (*j_r_src_)(const Vector&, Vector&);
-   void   (*j_i_src_)(const Vector&, Vector&);
-   // void   (*m_src_)(const Vector&, Vector&);
-
-   // Array of 0's and 1's marking the location of absorbing surfaces
-   Array<int> abc_marker_;
-
-   // Array of 0's and 1's marking the location of Dirichlet boundaries
-   Array<int> dbc_marker_;
-   void   (*e_r_bc_)(const Vector&, Vector&);
-   void   (*e_i_bc_)(const Vector&, Vector&);
-
-   Array<int> * dbcs_;
-   Array<int> ess_bdr_;
-   Array<int> ess_bdr_tdofs_;
-   Array<int> non_k_bdr_;
-
-   VisItDataCollection * visit_dc_;
-
-   std::map<std::string,socketstream*> socks_;
-};
-/*
-class SurfaceCurrent
-{
-public:
- SurfaceCurrent(ParFiniteElementSpace & H1FESpace,
-                ParDiscreteGradOperator & Grad,
-                Array<int> & kbcs, Array<int> & vbcs, Vector & vbcv);
- ~SurfaceCurrent();
-
- void InitSolver() const;
-
- void ComputeSurfaceCurrent(ParGridFunction & k);
-
- void Update();
-
- ParGridFunction * GetPsi() { return psi_; }
-
-private:
- int myid_;
-
- ParFiniteElementSpace   * H1FESpace_;
- ParDiscreteGradOperator * grad_;
- Array<int>              * kbcs_;
- Array<int>              * vbcs_;
- Vector                  * vbcv_;
-
- ParBilinearForm * s0_;
- ParGridFunction * psi_;
- ParGridFunction * rhs_;
-
- HypreParMatrix  * S0_;
- mutable Vector Psi_;
- mutable Vector RHS_;
-
- mutable HypreBoomerAMG  * amg_;
- mutable HyprePCG        * pcg_;
-
- Array<int> ess_bdr_, ess_bdr_tdofs_;
- Array<int> non_k_bdr_;
-};
-*/
-} // namespace electromagnetics
-
-} // namespace mfem
-
-#endif // MFEM_HERTZ_SOLVER

miniapps/electromagnetics/hertz_solver.cpp

@@ -1,711 +0,0 @@
-// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
-// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
-// reserved. See file COPYRIGHT for details.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability see http://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the GNU Lesser General Public License (as published by the Free
-// Software Foundation) version 2.1 dated February 1999.
-
-#include "hertz_solver.hpp"
-
-#ifdef MFEM_USE_MPI
-
-using namespace std;
-namespace mfem
-{
-using namespace common;
-
-namespace electromagnetics
-{
-
-// Used for combining scalar coefficients
-double prodFunc(double a, double b) { return a * b; }
-
-HertzSolver::HertzSolver(ParMesh & pmesh, int order, double freq,
-                         HertzSolver::SolverType sol, SolverOptions & sOpts,
-                         HertzSolver::PrecondType prec,
-                         ComplexOperator::Convention conv,
-                         Coefficient & epsCoef,
-                         Coefficient & muInvCoef,
-                         Coefficient * sigmaCoef,
-                         Coefficient * etaInvCoef,
-                         Array<int> & abcs,
-                         Array<int> & dbcs,
-                         void (*e_r_bc )(const Vector&, Vector&),
-                         void (*e_i_bc )(const Vector&, Vector&),
-                         void (*j_r_src)(const Vector&, Vector&),
-                         void (*j_i_src)(const Vector&, Vector&))
-   : myid_(0),
-     num_procs_(1),
-     order_(order),
-     logging_(1),
-     sol_(sol),
-     solOpts_(sOpts),
-     prec_(prec),
-     conv_(conv),
-     ownsEtaInv_(etaInvCoef == NULL),
-     freq_(freq),
-     pmesh_(&pmesh),
-     HCurlFESpace_(NULL),
-     a1_(NULL),
-     b1_(NULL),
-     e_(NULL),
-     e_t_(NULL),
-     j_(NULL),
-     jd_(NULL),
-     epsCoef_(&epsCoef),
-     muInvCoef_(&muInvCoef),
-     sigmaCoef_(sigmaCoef),
-     etaInvCoef_(etaInvCoef),
-     omegaCoef_(new ConstantCoefficient(2.0 * M_PI * freq_)),
-     negOmegaCoef_(new ConstantCoefficient(-2.0 * M_PI * freq_)),
-     omega2Coef_(new ConstantCoefficient(pow(2.0 * M_PI * freq_, 2))),
-     negOmega2Coef_(new ConstantCoefficient(-pow(2.0 * M_PI * freq_, 2))),
-     massCoef_(NULL),
-     posMassCoef_(NULL),
-     lossCoef_(NULL),
-     abcCoef_(NULL),
-     posAbcCoef_(NULL),
-     jrCoef_(NULL),
-     jiCoef_(NULL),
-     erCoef_(NULL),
-     eiCoef_(NULL),
-     j_r_src_(j_r_src),
-     j_i_src_(j_i_src),
-     e_r_bc_(e_r_bc),
-     e_i_bc_(e_i_bc),
-     dbcs_(&dbcs),
-     visit_dc_(NULL)
-{
-   // Initialize MPI variables
-   MPI_Comm_size(pmesh_->GetComm(), &num_procs_);
-   MPI_Comm_rank(pmesh_->GetComm(), &myid_);
-
-   // Define compatible parallel finite element spaces on the parallel
-   // mesh. Here we use arbitrary order Nedelec finite elements.
-   HCurlFESpace_ = new ND_ParFESpace(pmesh_,order,pmesh_->Dimension());
-
-   // Set the size of the 2x2 block representation of the complex linear system
-   blockTrueOffsets_.SetSize(3);
-   blockTrueOffsets_[0] = 0;
-   blockTrueOffsets_[1] = HCurlFESpace_->TrueVSize();
-   blockTrueOffsets_[2] = HCurlFESpace_->TrueVSize();
-   blockTrueOffsets_.PartialSum();
-
-   // Setup Dirichlet BC
-   ess_bdr_.SetSize(pmesh.bdr_attributes.Max());
-   if ( dbcs_ != NULL )
-   {
-      if ( dbcs_->Size() == 1 && (*dbcs_)[0] == -1 )
-      {
-         ess_bdr_ = 1;
-      }
-      else
-      {
-         ess_bdr_ = 0;
-         for (int i=0; i<dbcs_->Size(); i++)
-         {
-            ess_bdr_[(*dbcs_)[i]-1] = 1;
-         }
-      }
-      HCurlFESpace_->GetEssentialTrueDofs(ess_bdr_, ess_bdr_tdofs_);
-
-      if (e_r_bc_)
-      {
-         erCoef_ = new VectorFunctionCoefficient(pmesh_->SpaceDimension(),
-                                                 e_r_bc_);
-         if (e_i_bc_ == NULL)
-         {
-            Vector e(3); e = 0.0;
-            eiCoef_ = new VectorConstantCoefficient(e);
-         }
-      }
-      if (e_i_bc_)
-      {
-         eiCoef_ = new VectorFunctionCoefficient(pmesh_->SpaceDimension(),
-                                                 e_i_bc_);
-         if (e_r_bc_ == NULL)
-         {
-            Vector e(3); e = 0.0;
-            erCoef_ = new VectorConstantCoefficient(e);
-         }
-      }
-   }
-
-   // Setup various coefficients
-   massCoef_ = new TransformedCoefficient(negOmega2Coef_, epsCoef_, prodFunc);
-   posMassCoef_ = new TransformedCoefficient(omega2Coef_, epsCoef_, prodFunc);
-   if ( sigmaCoef_ )
-   {
-      lossCoef_ = new TransformedCoefficient(omegaCoef_, sigmaCoef_, prodFunc);
-   }
-
-   // Impedance of free space for the Absorbing boundary condition
-   if ( abcs.Size() > 0 )
-   {
-      if ( myid_ == 0 && logging_ > 0 )
-      {
-         cout << "Creating Admittance Coefficient" << endl;
-      }
-
-      abc_marker_.SetSize(pmesh.bdr_attributes.Max());
-      if ( abcs.Size() == 1 && abcs[0] < 0 )
-      {
-         // Mark all boundaries as absorbing
-         abc_marker_ = 1;
-      }
-      else
-      {
-         // Mark select boundaries as absorbing
-         abc_marker_ = 0;
-         for (int i=0; i<abcs.Size(); i++)
-         {
-            abc_marker_[abcs[i]-1] = 1;
-         }
-      }
-      if ( etaInvCoef_ == NULL )
-      {
-         etaInvCoef_ = new ConstantCoefficient(sqrt(epsilon0_/mu0_));
-      }
-      abcCoef_ = new TransformedCoefficient(negOmegaCoef_, etaInvCoef_,
-                                            prodFunc);
-      posAbcCoef_ = new TransformedCoefficient(omegaCoef_, etaInvCoef_,
-                                               prodFunc);
-   }
-
-   // Volume Current Density
-   if ( j_r_src_ != NULL )
-   {
-      jrCoef_ = new VectorFunctionCoefficient(pmesh_->SpaceDimension(),
-                                              j_r_src_);
-   }
-   else
-   {
-      Vector j(3); j = 0.0;
-      jrCoef_ = new VectorConstantCoefficient(j);
-   }
-   if ( j_i_src_ != NULL )
-   {
-      jiCoef_ = new VectorFunctionCoefficient(pmesh_->SpaceDimension(),
-                                              j_i_src_);
-   }
-   else
-   {
-      Vector j(3); j = 0.0;
-      jiCoef_ = new VectorConstantCoefficient(j);
-   }
-
-   // Bilinear Forms
-   // Primary system operator
-   a1_ = new ParSesquilinearForm(HCurlFESpace_, conv_);
-   a1_->AddDomainIntegrator(new CurlCurlIntegrator(*muInvCoef_), NULL);
-   a1_->AddDomainIntegrator(new VectorFEMassIntegrator(*massCoef_), NULL);
-   if ( lossCoef_ )
-   {
-      a1_->AddDomainIntegrator(NULL, new VectorFEMassIntegrator(*lossCoef_));
-   }
-   if ( abcCoef_ )
-   {
-      a1_->AddBoundaryIntegrator(NULL, new VectorFEMassIntegrator(*abcCoef_),
-                                 abc_marker_);
-   }
-
-   // Operator used with the perconditioner
-   b1_ = new ParBilinearForm(HCurlFESpace_);
-   b1_->AddDomainIntegrator(new CurlCurlIntegrator(*muInvCoef_));
-   b1_->AddDomainIntegrator(new VectorFEMassIntegrator(*posMassCoef_));
-   if ( lossCoef_ )
-   {
-      b1_->AddDomainIntegrator(new VectorFEMassIntegrator(*lossCoef_));
-   }
-   if ( abcCoef_ )
-   {
-      b1_->AddBoundaryIntegrator(new VectorFEMassIntegrator(*posAbcCoef_),
-                                 abc_marker_);
-   }
-
-   // Build grid functions
-   // The solution vector is the Electric field
-   e_   = new ParComplexGridFunction(HCurlFESpace_);
-   e_t_ = new ParGridFunction(HCurlFESpace_);
-   if (erCoef_ && eiCoef_)
-   {
-      e_->ProjectCoefficient(*erCoef_, *eiCoef_);
-   }
-   else
-   {
-      *e_ = 0.0;
-   }
-
-   // A GridFunction to visualize the volumetric current density
-   j_ = new ParComplexGridFunction(HCurlFESpace_);
-   j_->ProjectCoefficient(*jrCoef_, *jiCoef_);
-
-   // A LineatForm representation of the current denisty for the RHS
-   jd_ = new ParComplexLinearForm(HCurlFESpace_, conv_);
-   jd_->AddDomainIntegrator(new VectorFEDomainLFIntegrator(*jrCoef_),
-                            new VectorFEDomainLFIntegrator(*jiCoef_));
-   jd_->real().Vector::operator=(0.0);
-   jd_->imag().Vector::operator=(0.0);
-}
-
-HertzSolver::~HertzSolver()
-{
-   delete jrCoef_;
-   delete jiCoef_;
-   delete erCoef_;
-   delete eiCoef_;
-   delete massCoef_;
-   delete posMassCoef_;
-   delete lossCoef_;
-   delete abcCoef_;
-   delete posAbcCoef_;
-   if ( ownsEtaInv_ ) { delete etaInvCoef_; }
-   delete omegaCoef_;
-   delete negOmegaCoef_;
-   delete omega2Coef_;
-   delete negOmega2Coef_;
-
-   delete e_;
-   delete e_t_;
-   delete j_;
-   delete jd_;
-
-   delete a1_;
-   delete b1_;
-
-   delete HCurlFESpace_;
-
-   map<string,socketstream*>::iterator mit;
-   for (mit=socks_.begin(); mit!=socks_.end(); mit++)
-   {
-      delete mit->second;
-   }
-}
-
-HYPRE_Int
-HertzSolver::GetProblemSize()
-{
-   return 2 * HCurlFESpace_->GlobalTrueVSize();
-}
-
-void
-HertzSolver::PrintSizes()
-{
-   HYPRE_Int size_nd = HCurlFESpace_->GlobalTrueVSize();
-   if (myid_ == 0)
-   {
-      cout << "Number of H(Curl) unknowns: " << size_nd << endl;
-   }
-}
-
-void
-HertzSolver::Assemble()
-{
-   if ( myid_ == 0 && logging_ > 0 ) { cout << "Assembling ..." << flush; }
-
-   a1_->Assemble();
-   a1_->Finalize();
-
-   b1_->Assemble();
-   b1_->Finalize();
-
-   jd_->Assemble();
-
-   if ( myid_ == 0 && logging_ > 0 ) { cout << " done." << endl; }
-}
-
-void
-HertzSolver::Update()
-{
-   if ( myid_ == 0 && logging_ > 0 ) { cout << "Updating ..." << endl; }
-
-   // Inform the spaces that the mesh has changed
-   HCurlFESpace_->Update();
-
-   if ( ess_bdr_.Size() > 0 )
-   {
-      HCurlFESpace_->GetEssentialTrueDofs(ess_bdr_, ess_bdr_tdofs_);
-   }
-
-   blockTrueOffsets_[0] = 0;
-   blockTrueOffsets_[1] = HCurlFESpace_->TrueVSize();
-   blockTrueOffsets_[2] = HCurlFESpace_->TrueVSize();
-   blockTrueOffsets_.PartialSum();
-
-   // Inform the grid functions that the space has changed.
-   e_->Update();
-   if (erCoef_ && eiCoef_)
-   {
-      e_->ProjectCoefficient(*erCoef_, *eiCoef_);
-   }
-
-   j_->Update();
-   j_->ProjectCoefficient(*jrCoef_, *jiCoef_);
-
-   jd_->Update();
-
-   // Inform the bilinear forms that the space has changed.
-   a1_->Update();
-   b1_->Update();
-}
-
-void
-HertzSolver::Solve()
-{
-   if ( myid_ == 0 && logging_ > 0 ) { cout << "Running solver ... " << endl; }
-
-   OperatorHandle A1;
-   Vector E, RHS;
-
-   a1_->FormLinearSystem(ess_bdr_tdofs_, *e_, *jd_, A1, E, RHS);
-
-   OperatorHandle PCOp;
-   b1_->FormSystemMatrix(ess_bdr_tdofs_, PCOp);
-
-   tic_toc.Clear();
-   tic_toc.Start();
-
-   Operator * pcr = NULL;
-   Operator * pci = NULL;
-   BlockDiagonalPreconditioner * BDP = NULL;
-
-   if (sol_ == FGMRES || sol_ == MINRES)
-   {
-      switch (prec_)
-      {
-         case INVALID_PC:
-            if ( myid_ == 0 && logging_ > 0 )
-            {
-               cout << "No Preconditioner Requested" << endl;
-            }
-            break;
-         case DIAG_SCALE:
-            if ( myid_ == 0 && logging_ > 0 )
-            {
-               cout << "Diagonal Scaling Preconditioner Requested" << endl;
-            }
-            pcr = new HypreDiagScale(dynamic_cast<HypreParMatrix&>(*PCOp.Ptr()));
-            break;
-         case PARASAILS:
-            if ( myid_ == 0 && logging_ > 0 )
-            {
-               cout << "ParaSails Preconditioner Requested" << endl;
-            }
-            pcr = new HypreParaSails(dynamic_cast<HypreParMatrix&>(*PCOp.Ptr()));
-            dynamic_cast<HypreParaSails*>(pcr)->SetSymmetry(1);
-            break;
-         case EUCLID:
-            if ( myid_ == 0 && logging_ > 0 )
-            {
-               cout << "Euclid Preconditioner Requested" << endl;
-            }
-            pcr = new HypreEuclid(dynamic_cast<HypreParMatrix&>(*PCOp.Ptr()));
-            if (solOpts_.euLvl != 1)
-            {
-               HypreSolver * pc = dynamic_cast<HypreSolver*>(pcr);
-               HYPRE_EuclidSetLevel(*pc, solOpts_.euLvl);
-            }
-            break;
-         case AMS:
-            if ( myid_ == 0 && logging_ > 0 )
-            {
-               cout << "AMS Preconditioner Requested" << endl;
-            }
-            pcr = new HypreAMS(dynamic_cast<HypreParMatrix&>(*PCOp.Ptr()),
-                               HCurlFESpace_);
-            break;
-         default:
-            MFEM_ABORT("Requested preconditioner is not available.");
-            break;
-      }
-      pci = pcr;
-
-      if (pcr)
-      {
-         BDP = new BlockDiagonalPreconditioner(blockTrueOffsets_);
-         BDP->SetDiagonalBlock(0, pcr);
-         BDP->SetDiagonalBlock(1, pci);
-         BDP->owns_blocks = 0;
-      }
-   }
-
-   switch (sol_)
-   {
-      case GMRES:
-      {
-         if ( myid_ == 0 && logging_ > 0 )
-         {
-            cout << "GMRES Solver Requested" << endl;
-         }
-         GMRESSolver gmres(HCurlFESpace_->GetComm());
-         gmres.SetOperator(*A1.Ptr());
-         gmres.SetRelTol(solOpts_.relTol);
-         gmres.SetMaxIter(solOpts_.maxIter);
-         gmres.SetKDim(solOpts_.kDim);
-         gmres.SetPrintLevel(solOpts_.printLvl);
-
-         gmres.Mult(RHS, E);
-      }
-      break;
-      case FGMRES:
-      {
-         if ( myid_ == 0 && logging_ > 0 )
-         {
-            cout << "FGMRES Solver Requested" << endl;
-         }
-         FGMRESSolver fgmres(HCurlFESpace_->GetComm());
-         if (BDP) { fgmres.SetPreconditioner(*BDP); }
-         fgmres.SetOperator(*A1.Ptr());
-         fgmres.SetRelTol(solOpts_.relTol);
-         fgmres.SetMaxIter(solOpts_.maxIter);
-         fgmres.SetKDim(solOpts_.kDim);
-         fgmres.SetPrintLevel(solOpts_.printLvl);
-
-         fgmres.Mult(RHS, E);
-      }
-      break;
-      case MINRES:
-      {
-         if ( myid_ == 0 && logging_ > 0 )
-         {
-            cout << "MINRES Solver Requested" << endl;
-         }
-         MINRESSolver minres(HCurlFESpace_->GetComm());
-         if (BDP) { minres.SetPreconditioner(*BDP); }
-         minres.SetOperator(*A1.Ptr());
-         minres.SetRelTol(solOpts_.relTol);
-         minres.SetMaxIter(solOpts_.maxIter);
-         minres.SetPrintLevel(solOpts_.printLvl);
-
-         minres.Mult(RHS, E);
-      }
-      break;
-#ifdef MFEM_USE_SUPERLU
-      case SUPERLU:
-      {
-         if ( myid_ == 0 && logging_ > 0 )
-         {
-            cout << "SuperLU Solver Requested" << endl;
-         }
-         ComplexHypreParMatrix * A1Z = A1.As<ComplexHypreParMatrix>();
-         HypreParMatrix * A1C = A1Z->GetSystemMatrix();
-         SuperLURowLocMatrix A_SuperLU(*A1C);
-         SuperLUSolver solver(MPI_COMM_WORLD);
-         solver.SetOperator(A_SuperLU);
-         solver.Mult(RHS, E);
-         delete A1C;
-      }
-      break;
-#endif
-#ifdef MFEM_USE_STRUMPACK
-      case STRUMPACK:
-      {
-         if ( myid_ == 0 && logging_ > 0 )
-         {
-            cout << "STRUMPACK Solver Requested" << endl;
-         }
-         ComplexHypreParMatrix * A1Z = A1.As<ComplexHypreParMatrix>();
-         HypreParMatrix * A1C = A1Z->GetSystemMatrix();
-         STRUMPACKRowLocMatrix A_STRUMPACK(*A1C);
-         STRUMPACKSolver solver(0, NULL, MPI_COMM_WORLD);
-         solver.SetPrintFactorStatistics(true);
-         solver.SetPrintSolveStatistics(false);
-         solver.SetKrylovSolver(strumpack::KrylovSolver::DIRECT);
-         solver.SetReorderingStrategy(strumpack::ReorderingStrategy::METIS);
-         solver.DisableMatching();
-         solver.SetOperator(A_STRUMPACK);
-         solver.SetFromCommandLine();
-         solver.Mult(RHS, E);
-         delete A1C;
-      }
-      break;
-#endif
-      default:
-         break;
-   };
-
-   tic_toc.Stop();
-
-   e_->Distribute(E);
-
-   delete BDP;
-   if (pci != pcr) { delete pci; }
-   delete pcr;
-
-   if ( myid_ == 0 && logging_ > 0 )
-   {
-      cout << " Solver done in " << tic_toc.RealTime() << " seconds." << endl;
-   }
-}
-
-void
-HertzSolver::GetErrorEstimates(Vector & errors)
-{
-   if ( myid_ == 0 && logging_ > 0 )
-   { cout << "Estimating Error ... " << flush; }
-
-   // Space for the discontinuous (original) flux
-   CurlCurlIntegrator flux_integrator(*muInvCoef_);
-   RT_FECollection flux_fec(order_-1, pmesh_->SpaceDimension());
-   ParFiniteElementSpace flux_fes(pmesh_, &flux_fec);
-
-   // Space for the smoothed (conforming) flux
-   double norm_p = 1;
-   ND_FECollection smooth_flux_fec(order_, pmesh_->Dimension());
-   ParFiniteElementSpace smooth_flux_fes(pmesh_, &smooth_flux_fec);
-
-   L2ZZErrorEstimator(flux_integrator, e_->real(),
-                      smooth_flux_fes, flux_fes, errors, norm_p);
-
-   if ( myid_ == 0 && logging_ > 0 ) { cout << "done." << endl; }
-}
-
-void
-HertzSolver::RegisterVisItFields(VisItDataCollection & visit_dc)
-{
-   visit_dc_ = &visit_dc;
-
-   visit_dc.RegisterField("Re(E)", &e_->real());
-   visit_dc.RegisterField("Im(E)", &e_->imag());
-
-   if ( j_ )
-   {
-      visit_dc.RegisterField("Re(J)", &j_->real());
-      visit_dc.RegisterField("Im(J)", &j_->imag());
-   }
-}
-
-void
-HertzSolver::WriteVisItFields(int it)
-{
-   if ( visit_dc_ )
-   {
-      if (myid_ == 0) { cout << "Writing VisIt files ..." << flush; }
-
-      if ( j_ )
-      {
-         j_->ProjectCoefficient(*jrCoef_, *jiCoef_);
-      }
-
-      HYPRE_Int prob_size = this->GetProblemSize();
-      visit_dc_->SetCycle(it);
-      visit_dc_->SetTime(prob_size);
-      visit_dc_->Save();
-
-      if (myid_ == 0) { cout << " done." << endl; }
-   }
-}
-
-void
-HertzSolver::InitializeGLVis()
-{
-   if ( myid_ == 0 ) { cout << "Opening GLVis sockets." << endl; }
-
-   socks_["Er"] = new socketstream;
-   socks_["Er"]->precision(8);
-
-   socks_["Ei"] = new socketstream;
-   socks_["Ei"]->precision(8);
-
-   if ( j_ )
-   {
-      socks_["Jr"] = new socketstream;
-      socks_["Jr"]->precision(8);
-
-      socks_["Ji"] = new socketstream;
-      socks_["Ji"]->precision(8);
-   }
-
-   if ( myid_ == 0 ) { cout << "GLVis sockets open." << endl; }
-}
-
-void
-HertzSolver::DisplayToGLVis()
-{
-   if (myid_ == 0) { cout << "Sending data to GLVis ..." << flush; }
-
-   char vishost[] = "localhost";
-   int  visport   = 19916;
-
-   int Wx = 0, Wy = 0; // window position
-   int Ww = 350, Wh = 350; // window size
-   int offx = Ww+10, offy = Wh+45; // window offsets
-
-   VisualizeField(*socks_["Er"], vishost, visport,
-                  e_->real(), "Electric Field, Re(E)", Wx, Wy, Ww, Wh);
-   Wx += offx;
-
-   VisualizeField(*socks_["Ei"], vishost, visport,
-                  e_->imag(), "Electric Field, Im(E)", Wx, Wy, Ww, Wh);
-
-   Wx = 0; Wy += offy; // next line
-
-   if ( j_ )
-   {
-      j_->ProjectCoefficient(*jrCoef_, *jiCoef_);
-
-      VisualizeField(*socks_["Jr"], vishost, visport,
-                     j_->real(), "Current Density, Re(J)", Wx, Wy, Ww, Wh);
-      Wx += offx;
-      VisualizeField(*socks_["Ji"], vishost, visport,
-                     j_->imag(), "Current Density, Im(J)", Wx, Wy, Ww, Wh);
-   }
-
-   Wx = 0; Wy += offy; // next line
-
-   if (myid_ == 0) { cout << " done." << endl; }
-}
-
-void
-HertzSolver::DisplayAnimationToGLVis()
-{
-   if (myid_ == 0) { cout << "Sending animation data to GLVis ..." << flush; }
-
-   Vector zeroVec(3); zeroVec = 0.0;
-   VectorConstantCoefficient zeroCoef(zeroVec);
-
-   double norm_r = e_->real().ComputeMaxError(zeroCoef);
-   double norm_i = e_->imag().ComputeMaxError(zeroCoef);
-
-   *e_t_ = e_->real();
-
-   char vishost[] = "localhost";
-   int  visport   = 19916;
-   socketstream sol_sock(vishost, visport);
-   sol_sock << "parallel " << num_procs_ << " " << myid_ << "\n";
-   sol_sock.precision(8);
-   sol_sock << "solution\n" << *pmesh_ << *e_t_
-            << "window_title 'Harmonic Solution (t = 0.0 T)'"
-            << "valuerange 0.0 " << max(norm_r, norm_i) << "\n"
-            << "autoscale off\n"
-            << "keys cvvv\n"
-            << "pause\n" << flush;
-   if (myid_ == 0)
-      cout << "GLVis visualization paused."
-           << " Press space (in the GLVis window) to resume it.\n";
-   int num_frames = 24;
-   int i = 0;
-   while (sol_sock)
-   {
-      double t = (double)(i % num_frames) / num_frames;
-      ostringstream oss;
-      oss << "Harmonic Solution (t = " << t << " T)";
-
-      add( cos( 2.0 * M_PI * t), e_->real(),
-           sin( 2.0 * M_PI * t), e_->imag(), *e_t_);
-      sol_sock << "parallel " << num_procs_ << " " << myid_ << "\n";
-      sol_sock << "solution\n" << *pmesh_ << *e_t_
-               << "window_title '" << oss.str() << "'" << flush;
-      i++;
-   }
-}
-
-} // namespace electromagnetics
-
-} // namespace mfem
-
-#endif // MFEM_USE_MPI

miniapps/electromagnetics/hertz_solver.hpp

@@ -1,194 +0,0 @@
-// Copyright (c) 2010, Lawrence Livermore National Security, LLC. Produced at
-// the Lawrence Livermore National Laboratory. LLNL-CODE-443211. All Rights
-// reserved. See file COPYRIGHT for details.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability see http://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the GNU Lesser General Public License (as published by the Free
-// Software Foundation) version 2.1 dated February 1999.
-
-#ifndef MFEM_HERTZ_SOLVER
-#define MFEM_HERTZ_SOLVER
-
-#include "../common/pfem_extras.hpp"
-
-#ifdef MFEM_USE_MPI
-
-#include <string>
-#include <map>
-
-namespace mfem
-{
-
-using common::H1_ParFESpace;
-using common::ND_ParFESpace;
-using common::RT_ParFESpace;
-using common::ParDiscreteGradOperator;
-using common::ParDiscreteCurlOperator;
-using common::DivergenceFreeProjector;
-
-namespace electromagnetics
-{
-
-// Physical Constants
-// Permittivity of Free Space (units F/m)
-static const double epsilon0_ = 8.8541878176e-12;
-// Permeability of Free Space (units H/m)
-static const double mu0_ = 4.0e-7*M_PI;
-
-// Solver options
-struct SolverOptions
-{
-   int maxIter;
-   int kDim;
-   int printLvl;
-   double relTol;
-
-   // Euclid Options
-   int euLvl;
-};
-
-//class SurfaceCurrent;
-class HertzSolver
-{
-public:
-
-   enum PrecondType
-   {
-      INVALID_PC  = -1,
-      DIAG_SCALE  =  1,
-      PARASAILS   =  2,
-      EUCLID      =  3,
-      AMS         =  4
-   };
-
-   enum SolverType
-   {
-      INVALID   = -1,
-      GMRES     =  1,
-      FGMRES    =  2,
-      MINRES    =  3,
-      SUPERLU   =  4,
-      STRUMPACK =  5
-   };
-
-   HertzSolver(ParMesh & pmesh, int order, double freq,
-               HertzSolver::SolverType s, SolverOptions & sOpts,
-               HertzSolver::PrecondType p,
-               ComplexOperator::Convention conv,
-               Coefficient & epsCoef,
-               Coefficient & muInvCoef,
-               Coefficient * sigmaCoef,
-               Coefficient * etaInvCoef,
-               Array<int> & abcs,
-               Array<int> & dbcs,
-               void (*e_r_bc )(const Vector&, Vector&),
-               void (*e_i_bc )(const Vector&, Vector&),
-               void (*j_r_src)(const Vector&, Vector&),
-               void (*j_i_src)(const Vector&, Vector&));
-   ~HertzSolver();
-
-   HYPRE_Int GetProblemSize();
-
-   void PrintSizes();
-
-   void Assemble();
-
-   void Update();
-
-   void Solve();
-
-   void GetErrorEstimates(Vector & errors);
-
-   void RegisterVisItFields(VisItDataCollection & visit_dc);
-
-   void WriteVisItFields(int it = 0);
-
-   void InitializeGLVis();
-
-   void DisplayToGLVis();
-
-   void DisplayAnimationToGLVis();
-
-private:
-
-   int myid_;
-   int num_procs_;
-   int order_;
-   int logging_;
-
-   SolverType sol_;
-   SolverOptions & solOpts_;
-   PrecondType prec_;
-
-   ComplexOperator::Convention conv_;
-
-   bool ownsEtaInv_;
-
-   double freq_;
-
-   ParMesh * pmesh_;
-
-   ND_ParFESpace * HCurlFESpace_;
-
-   Array<HYPRE_Int> blockTrueOffsets_;
-
-   ParSesquilinearForm * a1_;
-   ParBilinearForm * b1_;
-
-   ParComplexGridFunction * e_;   // Complex electric field (HCurl)
-   ParGridFunction        * e_t_; // Real electric field (HCurl)
-   ParComplexGridFunction * j_;   // Complex current density (HCurl)
-
-   ParComplexLinearForm   * jd_; // Dual of complex current density (HCurl)
-
-   Coefficient       * epsCoef_;    // Dielectric Material Coefficient
-   Coefficient       * muInvCoef_;  // Dia/Paramagnetic Material Coefficient
-   Coefficient       * sigmaCoef_;  // Electrical Conductivity Coefficient
-   Coefficient       * etaInvCoef_; // Admittance Coefficient
-
-   Coefficient * omegaCoef_;     // omega expressed as a Coefficient
-   Coefficient * negOmegaCoef_;  // -omega expressed as a Coefficient
-   Coefficient * omega2Coef_;    // omega^2 expressed as a Coefficient
-   Coefficient * negOmega2Coef_; // -omega^2 expressed as a Coefficient
-   Coefficient * massCoef_;      // -omega^2 epsilon
-   Coefficient * posMassCoef_;   // omega^2 epsilon
-   Coefficient * lossCoef_;      // -omega sigma
-   Coefficient * abcCoef_;       // -omega eta^{-1}
-   Coefficient * posAbcCoef_;    // omega eta^{-1}
-
-   VectorCoefficient * jrCoef_;     // Volume Current Density Function
-   VectorCoefficient * jiCoef_;     // Volume Current Density Function
-   VectorCoefficient * erCoef_;     // Electric Field Boundary Condition
-   VectorCoefficient * eiCoef_;     // Electric Field Boundary Condition
-
-   void (*j_r_src_)(const Vector&, Vector&);
-   void (*j_i_src_)(const Vector&, Vector&);
-
-   // Array of 0's and 1's marking the location of absorbing surfaces
-   Array<int> abc_marker_;
-
-   // Array of 0's and 1's marking the location of Dirichlet boundaries
-   Array<int> dbc_marker_;
-   void (*e_r_bc_)(const Vector&, Vector&);
-   void (*e_i_bc_)(const Vector&, Vector&);
-
-   Array<int> * dbcs_;
-   Array<int> ess_bdr_;
-   Array<int> ess_bdr_tdofs_;
-   Array<int> non_k_bdr_;
-
-   VisItDataCollection * visit_dc_;
-
-   std::map<std::string,socketstream*> socks_;
-};
-
-} // namespace electromagnetics
-
-} // namespace mfem
-
-#endif // MFEM_USE_MPI
-
-#endif // MFEM_HERTZ_SOLVER

miniapps/electromagnetics/makefile

@@ -25,8 +25,8 @@ include $(DEFAULTS_MK)
 MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
-SEQ_MINIAPPS = hertz_ser
-PAR_MINIAPPS = volta tesla maxwell hertz joule
+SEQ_MINIAPPS =
+PAR_MINIAPPS = volta tesla maxwell joule
 ifeq ($(MFEM_USE_MPI),NO)
    MINIAPPS = $(SEQ_MINIAPPS)
 else
@@ -118,6 +118,4 @@ clean-build:
 	rm -rf *.dSYM *.TVD.*breakpoints
 
 clean-exec:
-	@rm -rf Volta-AMR* Tesla-AMR* Maxwell-Parallel* Joule_* \
-	Hertz-AMR-Parallel* Hertz-AMR-Serial*
-
+	@rm -rf Volta-AMR* Tesla-AMR* Maxwell-Parallel* Joule_*

miniapps/electromagnetics/volta_solver.cpp

@@ -448,7 +448,7 @@ VoltaSolver::Solve()
    divEpsGrad_->RecoverFEMSolution(Phi, *rhod_, *phi_);
 
    // Compute the negative Gradient of the solution vector.  This is
-   // the electric field corresponding to the scalar potential
+   // the magnetic field corresponding to the scalar potential
    // represented by phi.
    grad_->Mult(*phi_, *e_); *e_ *= -1.0;
 

miniapps/navier/navier_solver.hpp

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