Adding finite difference unit test for divergence of RT bases

functionality is not implemented.

.gitignore

@@ -220,6 +220,7 @@ miniapps/meshing/trimmer.mesh
 miniapps/meshing/optimized*
 miniapps/meshing/perturbed*
 miniapps/meshing/polar-nc.mesh
+miniapps/meshing/ext-mesh-mapping
 
 miniapps/mtop/parheat
 miniapps/mtop/ParHeat*

examples/test_pyramid_h1.cpp

@@ -0,0 +1,86 @@
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+double testFunc(const Vector &p)
+{
+   const double x = p[0];
+   const double y = p[1];
+   const double z = p[2];
+
+   return cos(M_PI * x) * cos(M_PI * y) * cos(M_PI * z);
+}
+
+int main(int argc, char *argv[])
+{
+   int e = (int)Element::PYRAMID;
+   int nx = 1;
+   int r = 3;
+   int o = 1;
+   // int pyrtype = 0;
+   double d = 0.0;
+   bool visualization = true;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&e, "-e", "--elem-type", "Element Type: [4,7]");
+   args.AddOption(&nx, "-n", "--n", "Num elems in 1D");
+   args.AddOption(&r, "-r", "--refine", "Number of refinements");
+   args.AddOption(&o, "-o", "--order", "Number of refinements");
+   // args.AddOption(&pyrtype, "-p", "--pyramid-type", "0-Bergot, 1-Fuentes");
+   args.AddOption(&d, "-d", "--deformation", "Mesh deformation [0,1)");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.ParseCheck();
+
+   FunctionCoefficient testCoef(testFunc);
+
+   H1_FECollection fec(o, 3, BasisType::GaussLobatto);
+
+   Vector errs(r+1); errs = -1.0;
+   Vector conv(r);
+   for (int i = 0; i <= r; i++)
+   {
+      int n = nx * pow(2, i);
+      Mesh mesh = Mesh::MakeCartesian3D(n,n,n,(Element::Type)e);
+
+      if (d > 0.0)
+      {
+         const double max = (double)(RAND_MAX) + 1.0;
+         const double h = 1.0 / n;
+
+         Vector disp(3*mesh.GetNV());
+         for (int j=0; j<disp.Size(); j++)
+         {
+            disp[j] = (2.0 * rand()/max - 1.0) * h * d;
+         }
+
+         mesh.MoveVertices(disp);
+      }
+
+      FiniteElementSpace fes(&mesh, &fec);
+
+      GridFunction x(&fes);
+      x.ProjectCoefficient(testCoef);
+      errs[i] = x.ComputeL2Error(testCoef);
+      cout << "DoFs / L2 Error / Conv: " << fes.GetNDofs() << " / " << errs[i];
+      if (i > 0)
+      {
+         conv[i-1] = errs[i-1] / errs[i];
+         cout << " / " << conv[i-1];
+      }
+      cout << endl;
+
+      if (visualization)
+      {
+         char vishost[] = "localhost";
+         int  visport   = 19916;
+         socketstream sol_sock(vishost, visport);
+         sol_sock.precision(8);
+         sol_sock << "solution\n" << mesh << x << flush;
+      }
+   }
+}

examples/test_pyramid_hcurl.cpp

@@ -0,0 +1,92 @@
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+void testFunc(const Vector &p, Vector &V)
+{
+   const double x = p[0];
+   const double y = p[1];
+   const double z = p[2];
+
+   V.SetSize(3);
+   V[0] = cos(M_PI * x) * cos(M_PI * y) * cos(M_PI * z);
+   V[1] = cos(M_PI * x) * cos(M_PI * y) * cos(M_PI * z);
+   V[2] = cos(M_PI * x) * cos(M_PI * y) * cos(M_PI * z);
+   V /= sqrt(3.0);
+}
+
+int main(int argc, char *argv[])
+{
+   int e = (int)Element::PYRAMID;
+   int nx = 1;
+   int r = 3;
+   int o = 1;
+   double d = 0.0;
+   bool visualization = true;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&e, "-e", "--elem-type", "Element Type: [4,7]");
+   args.AddOption(&nx, "-n", "--n", "Num elems in 1D");
+   args.AddOption(&r, "-r", "--refine", "Number of refinements");
+   args.AddOption(&o, "-o", "--order", "Number of refinements");
+   args.AddOption(&d, "-d", "--deformation", "Mesh deformation [0,1)");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.ParseCheck();
+
+   VectorFunctionCoefficient testCoef(3, testFunc);
+
+   ND1_3DFECollection fec1;
+   ND_FECollection fec2(o, 3);
+
+   FiniteElementCollection &fec = (o == 1) ?
+                                  (FiniteElementCollection&)fec1 : (FiniteElementCollection&)fec2;
+
+   Vector errs(r+1); errs = -1.0;
+   Vector conv(r);
+   for (int i = 0; i <= r; i++)
+   {
+      int n = nx * pow(2, i);
+      Mesh mesh = Mesh::MakeCartesian3D(n,n,n,(Element::Type)e);
+
+      if (d > 0.0)
+      {
+         const double max = (double)(RAND_MAX) + 1.0;
+         const double h = 1.0 / n;
+
+         Vector disp(3*mesh.GetNV());
+         for (int j=0; j<disp.Size(); j++)
+         {
+            disp[j] = (2.0 * rand()/max - 1.0) * h * d;
+         }
+
+         mesh.MoveVertices(disp);
+      }
+
+      FiniteElementSpace fes(&mesh, &fec);
+
+      GridFunction x(&fes);
+      x.ProjectCoefficient(testCoef);
+      errs[i] = x.ComputeL2Error(testCoef);
+      cout << "DoFs / L2 Error / Conv: " << fes.GetNDofs() << " / " << errs[i];
+      if (i > 0)
+      {
+         conv[i-1] = errs[i-1] / errs[i];
+         cout << " / " << conv[i-1];
+      }
+      cout << endl;
+
+      if (visualization)
+      {
+         char vishost[] = "localhost";
+         int  visport   = 19916;
+         socketstream sol_sock(vishost, visport);
+         sol_sock.precision(8);
+         sol_sock << "solution\n" << mesh << x << flush;
+      }
+   }
+}

examples/test_pyramid_hdiv.cpp

@@ -0,0 +1,92 @@
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+void testFunc(const Vector &p, Vector &V)
+{
+   const double x = p[0];
+   const double y = p[1];
+   const double z = p[2];
+
+   V.SetSize(3);
+   V[0] = cos(M_PI * x) * cos(M_PI * y) * cos(M_PI * z);
+   V[1] = cos(M_PI * x) * cos(M_PI * y) * cos(M_PI * z);
+   V[2] = cos(M_PI * x) * cos(M_PI * y) * cos(M_PI * z);
+   V /= sqrt(3.0);
+}
+
+int main(int argc, char *argv[])
+{
+   int e = (int)Element::PYRAMID;
+   int nx = 1;
+   int r = 3;
+   int o = 1;
+   double d = 0.0;
+   bool visualization = true;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&e, "-e", "--elem-type", "Element Type: [4,7]");
+   args.AddOption(&nx, "-n", "--n", "Num elems in 1D");
+   args.AddOption(&r, "-r", "--refine", "Number of refinements");
+   args.AddOption(&o, "-o", "--order", "Number of refinements");
+   args.AddOption(&d, "-d", "--deformation", "Mesh deformation [0,1)");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.ParseCheck();
+
+   VectorFunctionCoefficient testCoef(3, testFunc);
+
+   RT0_3DFECollection fec1;
+   RT_FECollection fec2(o,3);
+
+   FiniteElementCollection &fec = (o == 1) ?
+                                  (FiniteElementCollection&)fec1 : (FiniteElementCollection&)fec2;
+
+   Vector errs(r+1); errs = -1.0;
+   Vector conv(r);
+   for (int i = 0; i <= r; i++)
+   {
+      int n = nx * pow(2, i);
+      Mesh mesh = Mesh::MakeCartesian3D(n,n,n,(Element::Type)e);
+
+      if (d > 0.0)
+      {
+         const double max = (double)(RAND_MAX) + 1.0;
+         const double h = 1.0 / n;
+
+         Vector disp(3*mesh.GetNV());
+         for (int j=0; j<disp.Size(); j++)
+         {
+            disp[j] = (2.0 * rand()/max - 1.0) * h * d;
+         }
+
+         mesh.MoveVertices(disp);
+      }
+
+      FiniteElementSpace fes(&mesh, &fec);
+
+      GridFunction x(&fes);
+      x.ProjectCoefficient(testCoef);
+      errs[i] = x.ComputeL2Error(testCoef);
+      cout << "DoFs / L2 Error / Conv: " << fes.GetNDofs() << " / " << errs[i];
+      if (i > 0)
+      {
+         conv[i-1] = errs[i-1] / errs[i];
+         cout << " / " << conv[i-1];
+      }
+      cout << endl;
+
+      if (visualization)
+      {
+         char vishost[] = "localhost";
+         int  visport   = 19916;
+         socketstream sol_sock(vishost, visport);
+         sol_sock.precision(8);
+         sol_sock << "solution\n" << mesh << x << flush;
+      }
+   }
+}

fem/coefficient.cpp

@@ -669,6 +669,51 @@ void MatrixFunctionCoefficient::EvalSymmetric(Vector &K,
    }
 }
 
+MatrixGridFunctionCoefficient::MatrixGridFunctionCoefficient(
+   const GridFunction *gf, int h, int w) : MatrixCoefficient(h, w)
+{
+   GridFunc = gf;
+   if (gf)
+   {
+      MFEM_ASSERT(gf->VectorDim() == h*w, "GridFunction vector dimension"
+                  "must equate to the number of matrix entries (h*w)");
+   }
+}
+
+void MatrixGridFunctionCoefficient::SetGridFunction(const GridFunction *gf,
+                                                    int h, int w)
+{
+   GridFunc = gf;
+   height = h;
+   width = w;
+   if (gf)
+   {
+      MFEM_ASSERT(gf->VectorDim() == h*w, "GridFunction vector dimension"
+                  "must equate to the number of matrix entries (h*w)");
+   }
+}
+
+void MatrixGridFunctionCoefficient::Eval(DenseMatrix &K,
+                                         ElementTransformation &T,
+                                         const IntegrationPoint &ip)
+{
+   Vector V;
+   Mesh *gf_mesh = GridFunc->FESpace()->GetMesh();
+   if (T.mesh == gf_mesh)
+   {
+      GridFunc->GetVectorValue(T, ip, V);
+   }
+   else
+   {
+      IntegrationPoint coarse_ip;
+      ElementTransformation *coarse_T = RefinedToCoarse(*gf_mesh, T, ip, coarse_ip);
+      GridFunc->GetVectorValue(*coarse_T, coarse_ip, V);
+   }
+
+   K.SetSize(height, width);
+   K.Set(1.0, V.begin());
+}
+
 void SymmetricMatrixCoefficient::Eval(DenseMatrix &K, ElementTransformation &T,
                                       const IntegrationPoint &ip)
 {

fem/coefficient.hpp

@@ -1158,6 +1158,36 @@ public:
    virtual ~MatrixArrayCoefficient();
 };
 
+/** @brief A matrix coefficient defined using a long vector grid function that
+    will have it's vector components interpreted as an m x n matrix in column
+    major order.  */
+class MatrixGridFunctionCoefficient : public MatrixCoefficient
+{
+protected:
+   const GridFunction *GridFunc;
+
+public:
+   /** @brief Construct the coefficient with the vector grid function @a gf and set
+       the m x n matrix dimensions that the vector grid function will be used to
+       interpret the long vector as a matrix in column major order.  The vector
+       dimension of @a gf must equate to the product of m and n.*/
+   MatrixGridFunctionCoefficient(const GridFunction *gf, int h, int w);
+
+   /** @brief Set the vector grid function for this coefficient to @a gf and set
+       the m x n matrix dimensions that the vector grid function will be used to
+       interpret the long vector as a matrix in column major order.  The vector
+       dimension of @a gf must equate to the product of m and n.*/
+   void SetGridFunction(const GridFunction *gf, int h, int w);
+
+   ///  Returns a pointer to the grid function in this Coefficient
+   const GridFunction * GetGridFunction() const { return GridFunc; }
+
+   virtual void Eval(DenseMatrix &K, ElementTransformation &T,
+                     const IntegrationPoint &ip);
+
+   virtual ~MatrixGridFunctionCoefficient() { }
+};
+
 
 /** @brief Derived matrix coefficient that has the value of the parent matrix
     coefficient where it is active and is zero otherwise. */

fem/fe/fe_fixed_order.cpp

@@ -6171,6 +6171,346 @@ void RT0PyrFiniteElement::ProjectCurl(const FiniteElement &fe,
    }
 }
 
+const double RT1PyrFiniteElement::nk[15] =
+{0,0,-1, 0,-1,0, 1,0,1, 0,1,1, -1,0,0};
+
+RT1PyrFiniteElement::RT1PyrFiniteElement()
+   : VectorFiniteElement(3, Geometry::PYRAMID, 28, 2, H_DIV), dof2nk(dof)
+{
+   const int p = order - 1;
+
+   const double *iop = poly1d.OpenPoints(p);
+   const double *icp = poly1d.ClosedPoints(p + 1);
+   const double *bop = poly1d.OpenPoints(p);
+
+#ifndef MFEM_THREAD_SAFE
+   u.SetSize(dof, dim);
+   divu.SetSize(dof);
+#else
+   DenseMatrix u(dof, dim);
+#endif
+
+   int o = 0;
+   // quadrilateral face
+   for (int j = 0; j <= p; j++)
+      for (int i = 0; i <= p; i++)  // (3,2,1,0)
+      {
+         Nodes.IntPoint(o).Set3(bop[i], bop[p-j], 0.);
+         dof2nk[o++] = 0;
+      }
+   // triangular faces
+   for (int j = 0; j <= p; j++)
+      for (int i = 0; i + j <= p; i++)  // (0,1,4)
+      {
+         double w = bop[i] + bop[j] + bop[p-i-j];
+         Nodes.IntPoint(o).Set3(bop[i]/w, 0., bop[j]/w);
+         dof2nk[o++] = 1;
+      }
+   for (int j = 0; j <= p; j++)
+      for (int i = 0; i + j <= p; i++)  // (1,2,4)
+      {
+         double w = bop[i] + bop[j] + bop[p-i-j];
+         Nodes.IntPoint(o).Set3(1.-bop[j]/w, bop[i]/w, bop[j]/w);
+         dof2nk[o++] = 2;
+      }
+   for (int j = 0; j <= p; j++)
+      for (int i = p - j; i >= 0; i--)  // (2,3,4)
+      {
+         double w = bop[i] + bop[j] + bop[p-i-j];
+         Nodes.IntPoint(o).Set3(bop[i]/w, 1.0-bop[j]/w, bop[j]/w);
+         dof2nk[o++] = 3;
+      }
+   for (int j = 0; j <= p; j++)
+      for (int i = p - j; i >= 0; i--)  // (3,0,4)
+      {
+         double w = bop[i] + bop[j] + bop[p-i-j];
+         Nodes.IntPoint(o).Set3(0., bop[i]/w, bop[j]/w);
+         dof2nk[o++] = 4;
+      }
+
+   // interior
+   // x-components
+   for (int k = 0; k <= p; k++)
+      for (int j = 0; j <= p; j++)
+         for (int i = 1; i <= p; i++)
+         {
+            double w = 1.0 - iop[k];
+            Nodes.IntPoint(o).Set3(icp[i]*w, iop[j]*w, iop[k]);
+            dof2nk[o++] = 4;
+         }
+   // y-components
+   for (int k = 0; k <= p; k++)
+      for (int j = 1; j <= p; j++)
+         for (int i = 0; i <= p; i++)
+         {
+            double w = 1.0 - iop[k];
+            Nodes.IntPoint(o).Set3(iop[i]*w, icp[j]*w, iop[k]);
+            dof2nk[o++] = 1;
+         }
+   // z-components
+   for (int k = 1; k <= p; k++)
+      for (int j = 0; j <= p; j++)
+         for (int i = 0; i <= p; i++)
+         {
+            double w = 1.0 - icp[k];
+            Nodes.IntPoint(o).Set3(iop[i]*w, iop[j]*w, icp[k]);
+            dof2nk[o++] = 0;
+         }
+
+   DenseMatrix T(dof);
+
+   for (int m = 0; m < dof; m++)
+   {
+      const IntegrationPoint &ip = Nodes.IntPoint(m);
+      const Vector nm({nk[3*dof2nk[m]], nk[3*dof2nk[m]+1], nk[3*dof2nk[m]+2]});
+      calcBasis(ip, u);
+      u.Mult(nm, T.GetColumn(m));
+   }
+
+   Ti.Factor(T);
+}
+
+void RT1PyrFiniteElement::CalcVShape(const IntegrationPoint &ip,
+                                     DenseMatrix &shape) const
+{
+#ifdef MFEM_THREAD_SAFE
+   DenseMatrix u(dof, dim);
+#endif
+
+   calcBasis(ip, u);
+
+   Ti.Mult(u, shape);
+}
+
+void RT1PyrFiniteElement::CalcRawVShape(const IntegrationPoint &ip,
+                                        DenseMatrix &shape) const
+{
+   calcBasis(ip, shape);
+}
+
+void RT1PyrFiniteElement::CalcDivShape(const IntegrationPoint &ip,
+                                       Vector &divshape) const
+{
+#ifdef MFEM_THREAD_SAFE
+   Vector divu(dof);
+#endif
+
+   calcDivBasis(ip, divu);
+
+   Ti.Mult(divu, divshape);
+}
+
+void RT1PyrFiniteElement::CalcRawDivShape(const IntegrationPoint &ip,
+                                          Vector &divshape) const
+{
+   calcDivBasis(ip, divshape);
+}
+
+void RT1PyrFiniteElement::calcBasis(const IntegrationPoint &ip,
+                                    DenseMatrix &F) const
+{
+   const double x = ip.x;
+   const double y = ip.y;
+   const double z = ip.z, oz = 1.0 - z;
+   const double x2 = 2.0*ip.x, y2 = 2.0*ip.y, z2 = 2.0*ip.z;
+
+   const double tol = 1e-6;
+
+   F = 0.0;
+
+   if (oz <= tol)
+   {
+      // At the apex the basis functions are not single valued. The following
+      // values are computed in the limit x->(1-z)/2, y->(1-z)/2, z->1.
+      F( 4, 0) = -0.25; F( 4, 1) = -0.75; F( 4, 2) =  0.5;
+      F( 6, 0) = -0.5;  F( 6, 1) = -1.5;  F( 6, 2) =  1.0;
+      F( 7, 0) =  0.25; F( 7, 1) = -0.25; F( 7, 2) = -0.5;
+      F( 9, 0) =  0.5;  F( 9, 1) = -0.5;  F( 9, 2) = -1.0;
+      F(10, 0) =  0.75; F(10, 1) =  0.25; F(10, 2) = -0.5;
+      F(12, 0) =  1.5;  F(12, 1) =  0.5;  F(12, 2) = -1.0;
+      F(13, 0) =  0.25; F(13, 1) = -0.25; F(13, 2) =  0.5;
+      F(15, 0) =  0.5;  F(15, 1) = -0.5;  F(15, 2) =  1.0;
+      F(16, 0) =  0.0;  F(16, 1) = -0.5;  F(16, 2) =  0.0;
+      F(18, 0) =  0.5;  F(18, 1) =  0.0;  F(18, 2) =  0.0;
+      F(26, 0) = -0.5;  F(26, 1) =  0.0;  F(26, 2) =  0.0;
+      F(27, 0) =  0.0;  F(27, 1) =  0.5;  F(27, 2) =  0.0;
+   }
+   else
+   {
+      const double ozi = 1.0 / oz;
+
+      F( 0, 0) = -x;
+      F( 0, 1) = -y;
+      F( 0, 2) =  oz;
+
+      F( 1, 0) =  x * (oz - x2) * ozi;
+      F( 1, 1) =  y * (oz - x2) * ozi;
+      F( 1, 2) =  x2 - oz;
+
+
+      F( 2, 0) =  x * (oz - y2) * ozi;
+      F( 2, 1) =  y * (oz - y2) * ozi;
+      F( 2, 2) =  y2 - oz;
+
+      F( 3, 0) = -x * (x2 - oz) * ozi * (y2 - oz) * ozi;
+      F( 3, 1) = -y * (x2 - oz) * ozi * (y2 - oz) * ozi;
+      F( 3, 2) = (x2 - oz) * (y2 - oz) * ozi;
+
+      F( 4, 0) = -0.5 * x * z * ozi;
+      F( 4, 1) = -0.5 * (oz * (2.0 - y) - y) * ozi;
+      F( 4, 2) =  0.5 * z;
+
+      F( 5, 0) =  0.5 * x * z * (x2 - oz) * ozi * (y - oz) * ozi;
+      F( 5, 1) =  0.5 * (2.0 - y) * (x2 - oz) * (y - oz) * ozi;
+      F( 5, 2) = -0.5 * z * (x2 - oz) * (y - oz) * ozi;
+
+      F( 6, 0) = -0.5 * x * z * (2.0 - 3.0 * oz + y) * ozi;
+      F( 6, 1) = -(oz * ((3.0 * oz - y) * (0.5 * y - 1.0) + 2.0) - y) * ozi;
+      F( 6, 2) =  0.5 * z * (2.0 - 3.0 * oz + y);
+
+      F( 7, 0) =  0.5 * x * z * ozi;
+      F( 7, 1) = -0.5 * y * (2.0 - z) * ozi;
+      F( 7, 2) = -0.5 * z;
+
+      F( 8, 0) =  0.5 * x * y * z * (x2 - oz) * ozi * ozi;
+      F( 8, 1) = -0.5 * y * (1.0 + y) * (x2 - oz) * ozi;
+      F( 8, 2) = -0.5 * y * z * (x2 - oz) * ozi;
+
+      F( 9, 0) = -0.5 * x * z * (y - z2) * ozi;
+      F( 9, 1) =  0.5 * y * (oz * (2.0 * oz + y + 1.0) - 2.0) * ozi;
+      F( 9, 2) =  0.5 * z * (y - z2);
+
+      F(10, 0) = -0.5 * (oz * (x - 2.0) + x) * ozi;
+      F(10, 1) =  0.5 * y * z * ozi;
+      F(10, 2) = -0.5 * z;
+
+      F(11, 0) = -0.5 * (2.0 - x) * (x - oz) * (y2 - oz) * ozi;
+      F(11, 1) = -0.5 * y * z * (x - oz) * ozi * (y2 - oz) * ozi;
+      F(11, 2) =  0.5 * z * (x - oz) * (y2 - oz) * ozi;
+
+      F(12, 0) =  0.5 * (oz * ((3.0 * oz - x) * (x - 2.0) + 4.0) - x2) * ozi;
+      F(12, 1) = -0.5 * y * z * (3.0 * oz - x - 2.0) * ozi;
+      F(12, 2) =  0.5 * z * (3.0 * oz - x - 2.0);
+
+      F(13, 0) =  0.5 * x * (1.0 + oz) * ozi;
+      F(13, 1) = -0.5 * y * z * ozi;
+      F(13, 2) =  0.5 * z;
+
+      F(14, 0) =  0.5 * x * (1.0 + x) * (y2 - oz) * ozi;
+      F(14, 1) = -0.5 * x * y * z * (y2 - oz) * ozi * ozi;
+      F(14, 2) =  0.5 * x * z * (y2 - oz) * ozi;
+
+      F(15, 0) = -0.5 * x * (1.0 + x - (5.0 + x) * z + 2.0 * z * z) * ozi;
+      F(15, 1) = -0.5 * y * z * (2.0 * z - x) * ozi;
+      F(15, 2) = -0.5 * x * z + z * z;
+
+      F(16, 0) = -x * z * (y2 - oz) * ozi * ozi;
+      F(16, 1) = -y * (y + z2 - 1.0) * ozi;
+      F(16, 2) =  z * (y2 - oz) * ozi;
+
+      F(17, 0) = -x * z * ozi * (x2 - oz) * ozi * (y2 - oz) * ozi;
+      F(17, 1) = -y * (x2 - oz) * ozi * (y + z2 - 1.0) * ozi;
+      F(17, 2) =  z * (x2 - oz) * ozi * (y2 - oz) * ozi;
+
+      F(18, 0) =  x * (x + z2 - 1.0) * ozi;
+      F(18, 1) =  y * z * ozi * (x2 - oz) * ozi;
+      F(18, 2) = -z * (x2 - oz) * ozi;
+
+      F(19, 0) =  x * (y2 - oz) * ozi * (x + z2 - 1.0) * ozi;
+      F(19, 1) =  y * z * ozi * (x2 - oz) * ozi * (y2 - oz) * ozi;
+      F(19, 2) = -z * (x2 - oz) * ozi * (y2 - oz) * ozi;
+
+      F(20, 0) = -2.0 * x * ozi * (x - oz) * ozi * (y2 - oz) * ozi;
+      F(20, 1) =  2.0 * y * ozi * (x2 - oz) * ozi * (y - oz) * ozi;
+      F(20, 2) =  0.0;
+
+      F(21, 0) =  x * z;
+      F(21, 1) =  y * z;
+      F(21, 2) = -z * oz;
+
+      F(22, 0) =  x * z * (x2 - oz) * ozi;
+      F(22, 1) =  y * z * (x2 - oz) * ozi;
+      F(22, 2) = -z * (x2 - oz);
+
+      F(23, 0) =  x * z * (y2 - oz) * ozi;
+      F(23, 1) =  y * z * (y2 - oz) * ozi;
+      F(23, 2) = -z * (y2 - oz);
+
+      F(24, 0) =  x * z * (x2 - oz) * ozi * (y2 - oz) * ozi;
+      F(24, 1) =  y * z * (x2 - oz) * ozi * (y2 - oz) * ozi;
+      F(24, 2) = -z * (x2 - oz) * (y2 - oz) * ozi;
+
+      F(25, 0) = -x * ozi * x * ozi * z * (y2 - oz) * ozi;
+      F(25, 1) = -y * ozi * y * ozi * z * (x2 - oz) * ozi;
+      F(25, 2) =  z * (x2 - oz) * ozi * (y2 - oz) * ozi;
+
+      F(26, 0) = -x * z * ozi;
+      F(26, 1) = -y * ozi * z * (x2 - oz) * ozi;
+      F(26, 2) =  z * (x2 - oz) * ozi;
+
+      F(27, 0) =  x * ozi * z * (y2 - oz) * ozi;
+      F(27, 1) =  y * ozi * z;
+      F(27, 2) = -z * (y2 - oz) * ozi;
+   }
+}
+
+void RT1PyrFiniteElement::calcDivBasis(const IntegrationPoint &ip,
+                                       Vector &dF) const
+{
+   const double x = ip.x;
+   const double y = ip.y;
+   const double z = ip.z, oz = 1.0 - z;
+   const double x2 = 2.0*ip.x, y2 = 2.0*ip.y;
+
+   const double tol = 1e-6;
+
+   dF = 0.0;
+
+   if (oz <= tol)
+   {
+      // At the apex the divergence is not single valued. The following values
+      // are computed in the limit x->(1-z)/2, y->(1-z)/2, z->1.
+      dF( 0) = -3.0;
+      dF( 4) =  1.5;
+      dF( 6) =  3.75;
+      dF( 7) = -1.5;
+      dF( 9) = -3.75;
+      dF(10) = -1.5;
+      dF(12) = -3.75;
+      dF(13) =  1.5;
+      dF(15) =  3.75;
+      dF(21) =  3.0;
+   }
+   else
+   {
+      const double ozi = 1.0 / oz;
+
+      dF( 0) = -3.0;
+      dF( 1) =  3.0 * (oz - x2) * ozi;
+      dF( 2) =  3.0 * (oz - y2) * ozi;
+      dF( 3) = -3.0 * (oz - x2) * ozi * (oz - y2) * ozi;
+      dF( 4) =  1.5;
+      dF( 5) =  0.5 * (oz - x2) * ozi * (4.0 * oz * (y - oz) - y) * ozi;
+      dF( 6) =  0.5 * (4.0 * oz * (2.0 - 3.0 * oz + y) - y) * ozi;
+      dF( 7) = -1.5;
+      dF( 8) =  0.5 * (oz - x2) * ozi * (oz * (1.0 + 4.0 * y) - y) * ozi;
+      dF( 9) =  0.5 * (oz * (8.0 * oz + 4.0 * y - 7.0) - y) * ozi;
+      dF(10) = -1.5;
+      dF(11) =  0.5 * (4.0 * oz * (oz - x) + x) * ozi * (oz - y2) * ozi;
+      dF(12) =  0.5 * (4.0 * oz * (3.0 * oz - x - 2.0) + x) * ozi;
+      dF(13) =  1.5;
+      dF(14) = -0.5 * (oz * (1.0 + 4.0 * x) - x) * ozi * (oz - y2) * ozi;
+      dF(15) = -0.5 * (oz * (8.0 * oz + 4.0 * x - 7.0) - x) * ozi;
+      dF(21) = -(1.0 - 4.0 * z);
+      dF(22) =  (1.0 - 4.0 * z) * (oz - x2) * ozi;
+      dF(23) =  (1.0 - 4.0 * z) * (oz - y2) * ozi;
+      dF(24) = -(1.0 - 4.0 * z) * (oz - x2) * ozi * (oz - y2) * ozi;
+      dF(25) =  (oz - x2) * ozi * (oz - y2) * ozi;
+      dF(26) = -(oz - x2) * ozi;
+      dF(27) =  (oz - y2) * ozi;
+   }
+}
+
 RotTriLinearHexFiniteElement::RotTriLinearHexFiniteElement()
    : NodalFiniteElement(3, Geometry::CUBE, 6, 2, FunctionSpace::Qk)
 {

fem/fe/fe_fixed_order.hpp

@@ -1179,6 +1179,72 @@ public:
 };
 
 
+/// A 3D 1st order Raviert-Thomas element on a pyramid
+class RT1PyrFiniteElement : public VectorFiniteElement
+{
+private:
+   static const double nk[15];
+
+#ifndef MFEM_THREAD_SAFE
+   mutable DenseMatrix u;
+   mutable Vector      divu;
+#endif
+   Array<int> dof2nk;
+   DenseMatrixInverse Ti;
+
+   void calcBasis(const IntegrationPoint &ip,
+                  DenseMatrix &F) const;
+
+   void calcDivBasis(const IntegrationPoint &ip,
+                     Vector &dF) const;
+
+public:
+   /// Construct the RT0PyrFiniteElement
+   RT1PyrFiniteElement();
+
+   virtual void CalcVShape(const IntegrationPoint &ip,
+                           DenseMatrix &shape) const;
+
+   virtual void CalcVShape(ElementTransformation &Trans,
+                           DenseMatrix &shape) const
+   { CalcVShape_RT(Trans, shape); }
+
+   virtual void CalcDivShape(const IntegrationPoint &ip,
+                             Vector &divshape) const;
+
+   virtual void GetLocalInterpolation(ElementTransformation &Trans,
+                                      DenseMatrix &I) const
+   { LocalInterpolation_RT(*this, nk, dof2nk, Trans, I); }
+   virtual void GetLocalRestriction(ElementTransformation &Trans,
+                                    DenseMatrix &R) const
+   { LocalRestriction_RT(nk, dof2nk, Trans, R); }
+   virtual void GetTransferMatrix(const FiniteElement &fe,
+                                  ElementTransformation &Trans,
+                                  DenseMatrix &I) const
+   { LocalInterpolation_RT(CheckVectorFE(fe), nk, dof2nk, Trans, I); }
+   using FiniteElement::Project;
+   virtual void Project(VectorCoefficient &vc,
+                        ElementTransformation &Trans, Vector &dofs) const
+   { Project_RT(nk, dof2nk, vc, Trans, dofs); }
+   virtual void ProjectMatrixCoefficient(
+      MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const
+   { ProjectMatrixCoefficient_RT(nk, dof2nk, mc, T, dofs); }
+   virtual void Project(const FiniteElement &fe, ElementTransformation &Trans,
+                        DenseMatrix &I) const
+   { Project_RT(nk, dof2nk, fe, Trans, I); }
+   virtual void ProjectCurl(const FiniteElement &fe,
+                            ElementTransformation &Trans,
+                            DenseMatrix &curl) const
+   { ProjectCurl_RT(nk, dof2nk, fe, Trans, curl); }
+
+   void CalcRawVShape(const IntegrationPoint &ip,
+                      DenseMatrix &shape) const;
+
+   void CalcRawDivShape(const IntegrationPoint &ip,
+                        Vector &dshape) const;
+};
+
+
 class RotTriLinearHexFiniteElement : public NodalFiniteElement
 {
 public:

fem/fe_coll.cpp

@@ -2348,8 +2348,15 @@ RT_FECollection::RT_FECollection(const int order, const int dim,
       RT_Elements[Geometry::PRISM] = new RT_WedgeElement(p);
       RT_dof[Geometry::PRISM] = p*pp1*(3*p + 4)/2;
 
-      RT_Elements[Geometry::PYRAMID] = new RT0PyrFiniteElement(false);
-      RT_dof[Geometry::PYRAMID] = 0;
+      if (p == 0)
+      {
+         RT_Elements[Geometry::PYRAMID] = new RT0PyrFiniteElement(false);
+      }
+      else if (p == 1)
+      {
+         RT_Elements[Geometry::PYRAMID] = new RT1PyrFiniteElement();
+      }
+      RT_dof[Geometry::PYRAMID] = 3*p*pp1*pp1;
    }
    else
    {
@@ -2496,14 +2503,15 @@ void RT_FECollection::InitFaces(const int p, const int dim_,
 const FiniteElement *
 RT_FECollection::FiniteElementForGeometry(Geometry::Type GeomType) const
 {
-   if (GeomType != Geometry::PYRAMID || this->GetOrder() == 1)
+   if (GeomType != Geometry::PYRAMID ||
+       this->GetOrder() == 1 || this->GetOrder() == 2 )
    {
       return RT_Elements[GeomType];
    }
    else
    {
       MFEM_ABORT("RT Pyramid basis functions are not yet supported "
-                 "for order > 0.");
+                 "for order > 1.");
       return NULL;
    }
 }

fem/intrules.cpp

@@ -1676,7 +1676,7 @@ IntegrationRule *IntegrationRules::PyramidIntegrationRule(int Order)
       ipp.x = ipc.x * (1.0 - ipc.z);
       ipp.y = ipc.y * (1.0 - ipc.z);
       ipp.z = ipc.z;
-      ipp.weight = ipc.weight / 3.0;
+      ipp.weight = ipc.weight * pow(1.0 - ipc.z, 2);
    }
    return PyramidIntRules[Order];
 }

mesh/mesh.cpp

@@ -31,6 +31,7 @@
 #include <functional>
 #include <map>
 #include <set>
+#include <unordered_set>
 
 // Include the METIS header, if using version 5. If using METIS 4, the needed
 // declarations are inlined below, i.e. no header is needed.
@@ -1456,6 +1457,7 @@ void Mesh::InitTables()
 {
    el_to_edge =
       el_to_face = el_to_el = bel_to_edge = face_edge = edge_vertex = NULL;
+   vertex_to_el = vertex_to_face = face_to_vertex = vertex_to_edge = NULL;
 }
 
 void Mesh::SetEmpty()
@@ -1537,9 +1539,15 @@ void Mesh::Destroy()
 
 void Mesh::ResetLazyData()
 {
-   delete el_to_el;     el_to_el = NULL;
-   delete face_edge;    face_edge = NULL;
-   delete edge_vertex;  edge_vertex = NULL;
+   delete el_to_el;        el_to_el = NULL;
+   delete face_edge;       face_edge = NULL;
+   delete edge_vertex;     edge_vertex = NULL;
+   delete vertex_to_edge;  vertex_to_edge = NULL;
+   delete vertex_to_face;  vertex_to_face = NULL;
+   delete face_to_vertex;  face_to_vertex = NULL;
+   delete vertex_to_el;    vertex_to_el = NULL;
+
+
    DeleteGeometricFactors();
    nbInteriorFaces = -1;
    nbBoundaryFaces = -1;
@@ -3615,6 +3623,12 @@ Mesh::Mesh(const Mesh &mesh, bool copy_nodes)
    // Do NOT copy the face-to-edge Table, face_edge
    face_edge = NULL;
 
+   // Do NOT copy the vertex to edge face and element tables
+   vertex_to_edge = NULL;
+   vertex_to_face = NULL;
+   face_to_vertex = NULL;
+   vertex_to_el = NULL;
+
    // Copy the edge-to-vertex Table, edge_vertex
    edge_vertex = (mesh.edge_vertex) ? new Table(*mesh.edge_vertex) : NULL;
 
@@ -5972,6 +5986,7 @@ void Mesh::GetBdrElementEdges(int i, Array<int> &edges, Array<int> &cor) const
    }
 }
 
+
 void Mesh::GetFaceEdges(int i, Array<int> &edges, Array<int> &o) const
 {
    if (Dim == 2)
@@ -6065,13 +6080,72 @@ Table *Mesh::GetEdgeVertexTable() const
    return edge_vertex;
 }
 
-Table *Mesh::GetVertexToElementTable()
+
+Table *Mesh::GetVertexToEdgeTable() const
+{
+   if (vertex_to_edge) {return vertex_to_edge;}
+   if (!edge_vertex) {GetEdgeVertexTable();}
+   vertex_to_edge = Transpose(*edge_vertex);
+   vertex_to_edge->Finalize();
+   return vertex_to_edge;
+}
+
+Table *Mesh::GetFaceToVertexTable() const
 {
    int i, j, nv, *v;
 
-   Table *vert_elem = new Table;
+   if (face_to_vertex) {return face_to_vertex;}
+   face_to_vertex = new Table;
 
-   vert_elem->MakeI(NumOfVertices);
+   face_to_vertex->MakeI(NumOfFaces);
+
+   for (i = 0; i < NumOfFaces; i++)
+   {
+      nv = faces[i]->GetNVertices();
+      v  = faces[i]->GetVertices();
+      for (j = 0; j < nv; j++)
+      {
+         face_to_vertex->AddAColumnInRow(i);
+      }
+   }
+
+   face_to_vertex->MakeJ();
+
+   for (i = 0; i < NumOfFaces; i++)
+   {
+      nv = faces[i]->GetNVertices();
+      v  = faces[i]->GetVertices();
+      for (j = 0; j < nv; j++)
+      {
+         face_to_vertex->AddConnection(i, v[j]);
+      }
+   }
+
+   face_to_vertex->ShiftUpI();
+   face_to_vertex->Finalize();
+
+   return face_to_vertex;
+}
+
+
+Table *Mesh::GetVertexToFaceTable() const
+{
+   if (vertex_to_face) {return vertex_to_face;}
+   if (!face_to_vertex) {GetFaceToVertexTable();}
+   vertex_to_face = Transpose(*face_to_vertex);
+   vertex_to_face->Finalize();
+   return vertex_to_face;
+}
+
+
+Table *Mesh::GetVertexToElementTable() const
+{
+   int i, j, nv, *v;
+
+   if (vertex_to_el) {return vertex_to_el;}
+   vertex_to_el = new Table;
+
+   vertex_to_el->MakeI(NumOfVertices);
 
    for (i = 0; i < NumOfElements; i++)
    {
@@ -6079,11 +6153,11 @@ Table *Mesh::GetVertexToElementTable()
       v  = elements[i]->GetVertices();
       for (j = 0; j < nv; j++)
       {
-         vert_elem->AddAColumnInRow(v[j]);
+         vertex_to_el->AddAColumnInRow(v[j]);
       }
    }
 
-   vert_elem->MakeJ();
+   vertex_to_el->MakeJ();
 
    for (i = 0; i < NumOfElements; i++)
    {
@@ -6091,13 +6165,14 @@ Table *Mesh::GetVertexToElementTable()
       v  = elements[i]->GetVertices();
       for (j = 0; j < nv; j++)
       {
-         vert_elem->AddConnection(v[j], i);
+         vertex_to_el->AddConnection(v[j], i);
       }
    }
 
-   vert_elem->ShiftUpI();
+   vertex_to_el->ShiftUpI();
+   vertex_to_el->Finalize();
 
-   return vert_elem;
+   return vertex_to_el;
 }
 
 Table *Mesh::GetFaceToElementTable() const
@@ -6180,6 +6255,184 @@ void Mesh::GetBdrElementFace(int i, int *f, int *o) const
    }
 }
 
+void Mesh::ElemsWithVert(Array<int> &elems, int vi)
+{
+   if (!vertex_to_el) {GetVertexToElementTable();}
+   vertex_to_el->GetRow(vi, elems);
+}
+
+void Mesh::FacesWithVert(Array<int> &faces_, int vi)
+{
+   if (!vertex_to_face) {GetVertexToFaceTable();}
+   vertex_to_face->GetRow(vi, faces_);
+}
+
+void Mesh::EdgesWithVert(Array<int> &edges, int vi)
+{
+   if (!vertex_to_edge) {GetVertexToEdgeTable();}
+   vertex_to_edge->GetRow(vi, edges);
+}
+
+
+void Mesh::ElemsWithAllVerts(Array<int> &elems, const Array<int> &verts)
+{
+   if (!vertex_to_el) {GetVertexToElementTable();}
+   //elems.Reserve(verts.Size());
+
+   //Find all the elements touched by the vertices
+   std::set<int> touched_elems;
+   for (int i = 0; i < verts.Size(); ++i)
+   {
+      int row_sz = vertex_to_el->RowSize(verts[i]);
+      const int *row = vertex_to_el->GetRow(verts[i]);
+      for (int j = 0; j < row_sz; ++j)
+      {
+         touched_elems.insert(row[j]);
+      }
+   }
+
+   //Put the verts into a hashing set for fast finding
+   std::unordered_set<int> vert_set(verts.begin(), verts.end());
+
+   //Run through the touched elems and put the ones that fully covered in
+   int num_elems = 0;
+   elems.SetSize(touched_elems.size());
+   for (auto ei = touched_elems.begin(); ei != touched_elems.end(); ++ei)
+   {
+      Element *elem = elements[*ei];
+      bool elem_covered = true;
+      int *elem_verts = elem->GetVertices();
+      for (int j = 0; j < elem->GetNVertices(); ++j)
+      {
+         if (vert_set.count(elem_verts[j]) < 1)
+         {
+            elem_covered = false;
+            break;
+         }
+      }
+
+      if (elem_covered)
+      {
+         elems[num_elems] = *ei;
+         num_elems ++;
+      }
+   }
+   elems.SetSize(num_elems);
+}
+
+void Mesh::FacesWithAllVerts(Array<int> &faces_, const Array<int> &verts)
+{
+   if (!vertex_to_face) {GetVertexToFaceTable();}
+   if (!face_to_vertex) {GetFaceToVertexTable();}
+
+   //Find all the elements touched by the vertices
+   std::set<int> touched_faces;
+   for (int i = 0; i < verts.Size(); ++i)
+   {
+      int row_sz = vertex_to_face->RowSize(verts[i]);
+      const int *row = vertex_to_face->GetRow(verts[i]);
+      for (int j = 0; j < row_sz; ++j)
+      {
+         touched_faces.insert(row[j]);
+      }
+   }
+
+   //Put the verts into a hashing set for fast finding
+   std::unordered_set<int> vert_set(verts.begin(), verts.end());
+
+   //Run through the touched faces and put the ones that fully covered in
+   int num_faces = 0;
+   faces_.SetSize(touched_faces.size());
+   for (auto fi = touched_faces.begin(); fi != touched_faces.end(); ++fi)
+   {
+      int row_sz = face_to_vertex->RowSize(*fi);
+      const int *row = face_to_vertex->GetRow(*fi);
+      bool face_covered = true;
+      for (int j = 0; j < row_sz; ++j)
+      {
+         if (vert_set.count(row[j]) < 1)
+         {
+            face_covered = false;
+            break;
+         }
+      }
+
+      if (face_covered)
+      {
+         faces_[num_faces] = *fi;
+         num_faces ++;
+      }
+   }
+   faces_.SetSize(num_faces);
+}
+
+void Mesh::EdgesWithAllVerts(Array<int> &edges, const Array<int> &verts)
+{
+   if (!vertex_to_edge) {GetVertexToEdgeTable();}
+   if (!edge_vertex) {GetEdgeVertexTable();}
+
+   //Find all the elements touched by the vertices
+   std::set<int> touched_edges;
+   for (int i = 0; i < verts.Size(); ++i)
+   {
+      int row_sz = vertex_to_edge->RowSize(verts[i]);
+      const int *row = vertex_to_edge->GetRow(verts[i]);
+      for (int j = 0; j < row_sz; ++j)
+      {
+         touched_edges.insert(row[j]);
+      }
+   }
+
+   //Put the verts into a hashing set for fast finding
+   std::unordered_set<int> vert_set(verts.begin(), verts.end());
+
+   //Run through the touched faces and put the ones that fully covered in
+   int num_edges = 0;
+   edges.SetSize(touched_edges.size());
+   for (auto ei = touched_edges.begin(); ei != touched_edges.end(); ++ei)
+   {
+      int row_sz = edge_vertex->RowSize(*ei);
+      const int *row = edge_vertex->GetRow(*ei);
+      bool edge_covered = true;
+      for (int j = 0; j < row_sz; ++j)
+      {
+         if (vert_set.count(row[j]) < 1)
+         {
+            edge_covered = false;
+            break;
+         }
+      }
+
+      if (edge_covered)
+      {
+         edges[num_edges] = *ei;
+         num_edges ++;
+      }
+   }
+   edges.SetSize(num_edges);
+}
+
+
+void Mesh::EdgesInBdrElems(Array<int> &edges, const Array<int> &belems)
+{
+   std::set<int> touched_edges;
+   Array<int> be_edges;
+   Array<int> cor;
+   for (int bei = 0; bei < belems.Size(); ++bei)
+   {
+      GetBdrElementEdges(bei, be_edges, cor);
+      for (int edgei = 0; edgei < be_edges.Size(); ++edgei)
+      {
+         touched_edges.insert(be_edges[edgei]);
+      }
+   }
+
+   edges.SetSize(touched_edges.size());
+   std::copy(touched_edges.begin(), touched_edges.end(), edges.begin());
+}
+
+
+
 int Mesh::GetBdrElementEdgeIndex(int i) const
 {
    switch (Dim)

mesh/mesh.hpp

@@ -225,6 +225,11 @@ protected:
    Array<int> be_to_face;
    mutable Table *face_edge;
    mutable Table *edge_vertex;
+   mutable Table *vertex_to_el;
+   mutable Table *vertex_to_face;
+   mutable Table *face_to_vertex;
+   mutable Table *vertex_to_edge;
+
 
    IsoparametricTransformation Transformation, Transformation2;
    IsoparametricTransformation BdrTransformation;
@@ -1145,6 +1150,22 @@ public:
    /// Return the index and the orientation of the face of bdr element i. (3D)
    void GetBdrElementFace(int i, int *f, int *o) const;
 
+   /// Return the indices of the faces connected to vertex vi.
+   void FacesWithVert(Array<int> &faces, int vi);
+
+   /// Return the indices of the elements connected to vertex vi.
+   void ElemsWithVert(Array<int> &elems, int vi);
+
+   /// Return the indices of the edges conected to vertex vi.
+   void EdgesWithVert(Array<int> &edges, int vi);
+
+   void ElemsWithAllVerts(Array<int> &elems, const Array<int> &verts);
+   void FacesWithAllVerts(Array<int> &faces, const Array<int> &verts);
+   void EdgesWithAllVerts(Array<int> &edges, const Array<int> &verts);
+
+   /// Return the edges found in the list of boundary elements
+   void EdgesInBdrElems(Array<int> &edges, const Array<int> &belems);
+
    /** Return the vertex index of boundary element i. (1D)
        Return the edge index of boundary element i. (2D)
        Return the face index of boundary element i. (3D) */
@@ -1470,8 +1491,17 @@ public:
 
    const Table &ElementToEdgeTable() const;
 
-   ///  The returned Table must be destroyed by the caller
-   Table *GetVertexToElementTable();
+   /// Returns the vertex-to-element Table
+   Table *GetVertexToElementTable() const;
+
+   /// Returns the vertex-to-edge Table (3D)
+   Table *GetVertexToEdgeTable() const;
+
+   /// Returns the vertex-to-face Table (3D)
+   Table *GetVertexToFaceTable() const;
+
+   /// Returns the face_to_vertex Table (3D)
+   Table *GetFaceToVertexTable() const;
 
    /** Return the "face"-element Table. Here "face" refers to face (3D),
        edge (2D), or vertex (1D).

mesh/ncmesh.cpp

@@ -33,10 +33,11 @@ void NCMesh::GeomInfo::InitGeom(Geometry::Type geom)
    {
       case Geometry::CUBE: elem = new Hexahedron; break;
       case Geometry::PRISM: elem = new Wedge; break;
-      case Geometry::SQUARE: elem = new Quadrilateral; break;
-      case Geometry::SEGMENT: elem = new Segment; break;
-      case Geometry::TRIANGLE: elem = new Triangle; break;
       case Geometry::TETRAHEDRON: elem = new Tetrahedron; break;
+      case Geometry::PYRAMID: elem = new Pyramid; break;
+      case Geometry::SQUARE: elem = new Quadrilateral; break;
+      case Geometry::TRIANGLE: elem = new Triangle; break;
+      case Geometry::SEGMENT: elem = new Segment; break;
       default: MFEM_ABORT("unsupported geometry " << geom);
    }
 
@@ -97,7 +98,7 @@ static void CheckSupportedGeom(Geometry::Type geom)
    MFEM_VERIFY(geom == Geometry::SEGMENT ||
                geom == Geometry::TRIANGLE || geom == Geometry::SQUARE ||
                geom == Geometry::CUBE || geom == Geometry::PRISM ||
-               geom == Geometry::TETRAHEDRON,
+               geom == Geometry::PYRAMID || geom == Geometry::TETRAHEDRON,
                "Element type " << geom << " is not supported by NCMesh.");
 }
 
@@ -119,6 +120,13 @@ NCMesh::NCMesh(const Mesh *mesh)
       CheckSupportedGeom(geom);
       GI[geom].InitGeom(geom);
 
+      //If we have pyramids we will need tets after refinement
+      if (geom == Geometry::PYRAMID)
+      {
+         CheckSupportedGeom(Geometry::TETRAHEDRON);
+         GI[Geometry::TETRAHEDRON].InitGeom(Geometry::TETRAHEDRON);
+      }
+
       // create NCMesh::Element for this mfem::Element
       int root_id = AddElement(Element(geom, elem->GetAttribute()));
       MFEM_ASSERT(root_id == i, "");
@@ -454,7 +462,7 @@ NCMesh::Element::Element(Geometry::Type geom, int attr)
    : geom(geom), ref_type(0), tet_type(0), flag(0), index(-1)
    , rank(0), attribute(attr), parent(-1)
 {
-   for (int i = 0; i < 8; i++) { node[i] = -1; }
+   for (int i = 0; i < MaxElemNodes; i++) { node[i] = -1; }
 
    // NOTE: in 2D the 8-element node/child arrays are not optimal, however,
    // testing shows we would only save 17% of the total NCMesh memory if
@@ -476,7 +484,7 @@ int NCMesh::NewHexahedron(int n0, int n1, int n2, int n3,
    el.node[4] = n4, el.node[5] = n5, el.node[6] = n6, el.node[7] = n7;
 
    // get faces and assign face attributes
-   Face* f[6];
+   Face* f[MaxElemFaces];
    const GeomInfo &gi_hex = GI[Geometry::CUBE];
    for (int i = 0; i < gi_hex.nf; i++)
    {
@@ -549,6 +557,35 @@ int NCMesh::NewTetrahedron(int n0, int n1, int n2, int n3, int attr,
 
    return new_id;
 }
+int NCMesh::NewPyramid(int n0, int n1, int n2, int n3, int n4, int attr,
+                       int fattr0, int fattr1, int fattr2, int fattr3,
+                       int fattr4)
+{
+   // create new element, initialize nodes
+   int new_id = AddElement(Element(Geometry::PYRAMID, attr));
+   Element &el = elements[new_id];
+
+   el.node[0] = n0, el.node[1] = n1, el.node[2] = n2, el.node[3] = n3;
+   el.node[4] = n4;
+
+   // get faces and assign face attributes
+   Face* f[5];
+   const GeomInfo &gi_pyr = GI[Geometry::PYRAMID];
+   for (int i = 0; i < gi_pyr.nf; i++)
+   {
+      const int* fv = gi_pyr.faces[i];
+      f[i] = faces.Get(el.node[fv[0]], el.node[fv[1]],
+                       el.node[fv[2]], el.node[fv[3]]);
+   }
+
+   f[0]->attribute = fattr0;
+   f[1]->attribute = fattr1;
+   f[2]->attribute = fattr2;
+   f[3]->attribute = fattr3;
+   f[4]->attribute = fattr4;
+
+   return new_id;
+}
 
 int NCMesh::NewQuadrilateral(int n0, int n1, int n2, int n3,
                              int attr,
@@ -771,6 +808,30 @@ void NCMesh::CheckAnisoPrism(int vn1, int vn2, int vn3, int vn4,
    }
 }
 
+void NCMesh::CheckAnisoPyramid(int vn1, int vn2, int vn3, int vn4,
+                               const Refinement *refs, int nref)
+{
+   MeshId buf[4];
+   Array<MeshId> eid(buf, 4);
+   FindEdgeElements(vn1, vn2, vn3, vn4, eid);
+
+   // see if there is an element that has not been force-refined yet
+   for (int i = 0, j; i < eid.Size(); i++)
+   {
+      int elem = eid[i].element;
+      for (j = 0; j < nref; j++)
+      {
+         if (refs[j].index == elem) { break; }
+      }
+      if (j == nref) // elem not found in refs[]
+      {
+         // schedule prism refinement along Z axis
+         MFEM_ASSERT(elements[elem].Geom() == Geometry::PYRAMID, "");
+         ref_stack.Append(Refinement(elem, 4));
+      }
+   }
+}
+
 
 void NCMesh::CheckAnisoFace(int vn1, int vn2, int vn3, int vn4,
                             int mid12, int mid34, int level)
@@ -828,6 +889,21 @@ void NCMesh::CheckAnisoFace(int vn1, int vn2, int vn3, int vn4,
                CheckAnisoPrism(mid23, vn3, vn4, mid41, NULL, 0);
             }
          }
+         if (HavePyramids() && nodes[midf].HasEdge())
+         {
+            // Check if there is a pyramid with edge (mid23, mid41) that we may
+            // have missed in 'CheckAnisoFace', and force-refine it if present.
+
+            if (ref_stack.Size() > rs)
+            {
+               CheckAnisoPyramid(mid23, vn3, vn4, mid41,
+                                 &ref_stack[rs], ref_stack.Size() - rs);
+            }
+            else
+            {
+               CheckAnisoPyramid(mid23, vn3, vn4, mid41, NULL, 0);
+            }
+         }
 
          // perform the reparents all at once at the end
          if (level == 0)
@@ -895,7 +971,7 @@ void NCMesh::RefineElement(int elem, char ref_type)
       char remaining = ref_type & ~el.ref_type;
 
       // do the remaining splits on the children
-      for (int i = 0; i < 8; i++)
+      for (int i = 0; i < MaxElemChildren; i++)
       {
          if (el.child[i] >= 0) { RefineElement(el.child[i], remaining); }
       }
@@ -912,11 +988,11 @@ void NCMesh::RefineElement(int elem, char ref_type)
    int* no = el.node;
    int attr = el.attribute;
 
-   int child[8];
-   for (int i = 0; i < 8; i++) { child[i] = -1; }
+   int child[MaxElemChildren];
+   for (int i = 0; i < MaxElemChildren; i++) { child[i] = -1; }
 
    // get parent's face attributes
-   int fa[6];
+   int fa[MaxElemFaces];
    GeomInfo& gi = GI[el.Geom()];
    for (int i = 0; i < gi.nf; i++)
    {
@@ -1404,6 +1480,56 @@ void NCMesh::RefineElement(int elem, char ref_type)
                                    -1, -1, fa[3], -1);
       }
    }
+   else if (el.Geom() == Geometry::PYRAMID)
+   {
+      // Pyramid vertex numbering:
+      //
+      //    4
+      //     + \_                       Faces: 0 bottom (3,2,1,0)
+      //     |\\_ \_                           1 front (0, 1, 4)
+      //     ||  \_ \__                        2 right (1, 2, 4)
+      //     | \   \_   \__                    3 back (2, 3, 4)
+      //     |  +____\_ ____\                  4 left (3, 0, 4)
+      //     | /3      \_   2   Z  Y
+      //     |/          \ /    | /
+      //     +------------+     *--X
+      //    0              1
+
+      ref_type = Refinement::XYZ; // for consistence
+
+      int mid01 = GetMidEdgeNode(no[0], no[1]);
+      int mid12 = GetMidEdgeNode(no[1], no[2]);
+      int mid23 = GetMidEdgeNode(no[2], no[3]);
+      int mid03 = GetMidEdgeNode(no[0], no[3]);
+      int mid04 = GetMidEdgeNode(no[0], no[4]);
+      int mid14 = GetMidEdgeNode(no[1], no[4]);
+      int mid24 = GetMidEdgeNode(no[2], no[4]);
+      int mid34 = GetMidEdgeNode(no[3], no[4]);
+      int midf0 = GetMidFaceNode(mid23, mid12, mid01, mid03);
+
+      child[0] = NewPyramid(no[0], mid01, midf0, mid03, mid04,
+                            attr, fa[0], fa[1], -1, -1, fa[4]);
+      child[1] = NewTetrahedron(mid01, midf0, mid04, mid14,
+                                attr, -1, -1, -1, fa[1]);
+      child[2] = NewPyramid(mid01, no[1], mid12, midf0, mid14,
+                            attr, fa[0], fa[1], fa[2], -1, -1);
+      child[3] = NewTetrahedron(midf0, mid14, mid12, mid24,
+                                attr, -1, -1, fa[2], -1);
+      child[4] = NewPyramid(midf0, mid12, no[2], mid23, mid24,
+                            attr, fa[0], -1, fa[2], fa[3], -1);
+      child[5] = NewTetrahedron(midf0, mid23, mid34, mid24,
+                                attr, -1, -1, fa[3], -1);
+      child[6] = NewPyramid(mid03, midf0, mid23, no[3], mid34,
+                            attr, fa[0], -1, -1, fa[3], fa[4]);
+      child[7] = NewTetrahedron(mid03, mid04, midf0, mid34,
+                                attr, -1, fa[4], -1, -1);
+      child[8] = NewPyramid(mid24, mid14, mid04, mid34, midf0,
+                            attr, -1, -1, -1, -1, -1);
+      child[9] = NewPyramid(mid04, mid14, mid24, mid34, no[4],
+                            attr, -1, fa[1], fa[2], fa[3], fa[4]);
+
+      CheckIsoFace(no[3], no[2], no[1], no[0], mid23, mid12, mid01, mid03, midf0);
+   }
    else if (el.Geom() == Geometry::SQUARE)
    {
       ref_type &= 0x3; // ignore Z bit
@@ -1486,20 +1612,20 @@ void NCMesh::RefineElement(int elem, char ref_type)
    }
 
    // start using the nodes of the children, create edges & faces
-   for (int i = 0; i < 8 && child[i] >= 0; i++)
+   for (int i = 0; i < MaxElemChildren && child[i] >= 0; i++)
    {
       ReferenceElement(child[i]);
    }
 
-   int buf[6];
-   Array<int> parentFaces(buf, 6);
+   int buf[MaxElemFaces];
+   Array<int> parentFaces(buf, MaxElemFaces);
    parentFaces.SetSize(0);
 
    // sign off of all nodes of the parent, clean up unused nodes, but keep faces
    UnreferenceElement(elem, parentFaces);
 
    // register the children in their faces
-   for (int i = 0; i < 8 && child[i] >= 0; i++)
+   for (int i = 0; i < MaxElemChildren && child[i] >= 0; i++)
    {
       RegisterFaces(child[i]);
    }
@@ -1508,7 +1634,7 @@ void NCMesh::RefineElement(int elem, char ref_type)
    DeleteUnusedFaces(parentFaces);
 
    // make the children inherit our rank; set the parent element
-   for (int i = 0; i < 8 && child[i] >= 0; i++)
+   for (int i = 0; i < MaxElemChildren && child[i] >= 0; i++)
    {
       Element &ch = elements[child[i]];
       ch.rank = el.rank;
@@ -1586,6 +1712,12 @@ int NCMesh::RetrieveNode(const Element &el, int index)
          ch = el.child[ch];
          break;
 
+      case Geometry::PYRAMID:
+         ch = pyramid_deref_table[el.ref_type - 1][index];
+         MFEM_ASSERT(ch != -1, "");
+         ch = el.child[ch];
+         break;
+
       case Geometry::SQUARE:
          ch = el.child[quad_deref_table[el.ref_type - 1][index]];
          break;
@@ -1608,11 +1740,11 @@ void NCMesh::DerefineElement(int elem)
    Element &el = elements[elem];
    if (!el.ref_type) { return; }
 
-   int child[8];
+   int child[MaxElemChildren];
    std::memcpy(child, el.child, sizeof(child));
 
    // first make sure that all children are leaves, derefine them if not
-   for (int i = 0; i < 8 && child[i] >= 0; i++)
+   for (int i = 0; i < MaxElemChildren && child[i] >= 0; i++)
    {
       if (elements[child[i]].ref_type)
       {
@@ -1620,10 +1752,10 @@ void NCMesh::DerefineElement(int elem)
       }
    }
 
-   int faces_attribute[6];
+   int faces_attribute[MaxElemFaces];
    int ref_type_key = el.ref_type - 1;
 
-   for (int i = 0; i < 8; i++) { el.node[i] = -1; }
+   for (int i = 0; i < MaxElemNodes; i++) { el.node[i] = -1; }
 
    // retrieve original corner nodes and face attributes from the children
    if (el.Geom() == Geometry::CUBE)
@@ -1678,6 +1810,34 @@ void NCMesh::DerefineElement(int elem)
                               ->attribute;
       }
    }
+   else if (el.Geom() == Geometry::PYRAMID)
+   {
+      MFEM_ASSERT(pyramid_deref_table[ref_type_key][0] != -1,
+                  "invalid pyramid refinement");
+      constexpr int nb_pyramid_childs = 5;
+      for (int i = 0; i < nb_pyramid_childs; i++)
+      {
+         const int child_local_index = pyramid_deref_table[ref_type_key][i];
+         const int child_global_index = child[child_local_index];
+         Element &ch = elements[child_global_index];
+         el.node[i] = ch.node[i];
+      }
+      el.node[5] = el.node[6] = el.node[7] = -1;
+
+
+      constexpr int nb_pyramid_faces = 5;
+      for (int i = 0; i < nb_pyramid_faces; i++)
+      {
+         const int child_local_index = pyramid_deref_table[ref_type_key]
+                                       [i + nb_pyramid_childs];
+         const int child_global_index = child[child_local_index];
+         Element &ch = elements[child_global_index];
+         const int* fv = GI[el.Geom()].faces[i];
+         faces_attribute[i] = faces.Find(ch.node[fv[0]], ch.node[fv[1]],
+                                         ch.node[fv[2]], ch.node[fv[3]])
+                              ->attribute;
+      }
+   }
    else if (el.Geom() == Geometry::TETRAHEDRON)
    {
       for (int i = 0; i < 4; i++)
@@ -1745,13 +1905,13 @@ void NCMesh::DerefineElement(int elem)
    // sign in to all nodes
    ReferenceElement(elem);
 
-   int buf[8*6];
-   Array<int> childFaces(buf, 8*6);
+   int buf[MaxElemChildren*MaxElemFaces];
+   Array<int> childFaces(buf, MaxElemChildren*MaxElemFaces);
    childFaces.SetSize(0);
 
    // delete children, determine rank
    el.rank = std::numeric_limits<int>::max();
-   for (int i = 0; i < 8 && child[i] >= 0; i++)
+   for (int i = 0; i < MaxElemChildren && child[i] >= 0; i++)
    {
       el.rank = std::min(el.rank, elements[child[i]].rank);
       UnreferenceElement(child[i], childFaces);
@@ -1775,7 +1935,7 @@ void NCMesh::CollectDerefinements(int elem, Array<Connection> &list)
    if (!el.ref_type) { return; }
 
    int total = 0, ref = 0, ghost = 0;
-   for (int i = 0; i < 8 && el.child[i] >= 0; i++)
+   for (int i = 0; i < MaxElemChildren && el.child[i] >= 0; i++)
    {
       total++;
       Element &ch = elements[el.child[i]];
@@ -1787,7 +1947,7 @@ void NCMesh::CollectDerefinements(int elem, Array<Connection> &list)
    {
       // can be derefined, add to list
       int next_row = list.Size() ? (list.Last().from + 1) : 0;
-      for (int i = 0; i < 8 && el.child[i] >= 0; i++)
+      for (int i = 0; i < MaxElemChildren && el.child[i] >= 0; i++)
       {
          Element &ch = elements[el.child[i]];
          list.Append(Connection(next_row, ch.index));
@@ -1795,7 +1955,7 @@ void NCMesh::CollectDerefinements(int elem, Array<Connection> &list)
    }
    else
    {
-      for (int i = 0; i < 8 && el.child[i] >= 0; i++)
+      for (int i = 0; i < MaxElemChildren && el.child[i] >= 0; i++)
       {
          CollectDerefinements(el.child[i], list);
       }
@@ -1905,7 +2065,7 @@ void NCMesh::SetDerefMatrixCodes(int parent, Array<int> &fine_coarse)
 {
    // encode the ref_type and child number for GetDerefinementTransforms()
    Element &prn = elements[parent];
-   for (int i = 0; i < 8 && prn.child[i] >= 0; i++)
+   for (int i = 0; i < MaxElemChildren && prn.child[i] >= 0; i++)
    {
       Element &ch = elements[prn.child[i]];
       if (ch.index >= 0)
@@ -1977,7 +2137,7 @@ void NCMesh::CollectLeafElements(int elem, int state, Array<int> &ghosts,
       }
       else // no space filling curve tables yet for remaining cases
       {
-         for (int i = 0; i < 8; i++)
+         for (int i = 0; i < MaxElemChildren; i++)
          {
             if (el.child[i] >= 0)
             {
@@ -2231,7 +2391,8 @@ void NCMesh::InitRootState(int root_count)
       if (v_in < 0) { v_in = 0; }
 
       // determine which nodes are shared with the next element
-      bool shared[8] = { 0, 0, 0, 0, 0, 0, 0, 0 };
+      bool shared[MaxElemNodes];
+      for (int ni = 0; ni < MaxElemNodes; ++ni) { shared[ni] = 0; }
       if (i+1 < root_count)
       {
          Element &next = elements[i+1];
@@ -2265,6 +2426,7 @@ mfem::Element* NCMesh::NewMeshElement(int geom) const
    {
       case Geometry::CUBE: return new mfem::Hexahedron;
       case Geometry::PRISM: return new mfem::Wedge;
+      case Geometry::PYRAMID: return new mfem::Pyramid;
       case Geometry::TETRAHEDRON: return new mfem::Tetrahedron;
       case Geometry::SQUARE: return new mfem::Quadrilateral;
       case Geometry::TRIANGLE: return new mfem::Triangle;
@@ -2347,7 +2509,8 @@ void NCMesh::GetMeshComponents(Mesh &mesh) const
          if (face->Boundary())
          {
             if ((nc_elem.geom == Geometry::CUBE) ||
-                (nc_elem.geom == Geometry::PRISM && nfv == 4))
+                ((nc_elem.geom == Geometry::PRISM ||
+                  nc_elem.geom == Geometry::PYRAMID) && nfv == 4))
             {
                auto* quad = (Quadrilateral*) mesh.NewElement(Geometry::SQUARE);
                quad->SetAttribute(face->attribute);
@@ -2358,6 +2521,7 @@ void NCMesh::GetMeshComponents(Mesh &mesh) const
                mesh.boundary.Append(quad);
             }
             else if (nc_elem.geom == Geometry::PRISM ||
+                     nc_elem.geom == Geometry::PYRAMID ||
                      nc_elem.geom == Geometry::TETRAHEDRON)
             {
                MFEM_ASSERT(nfv == 3, "");
@@ -2603,7 +2767,7 @@ bool NCMesh::TriFaceSplit(int v1, int v2, int v3, int mid[3]) const
 
 int NCMesh::find_node(const Element &el, int node)
 {
-   for (int i = 0; i < 8; i++)
+   for (int i = 0; i < MaxElemNodes; i++)
    {
       if (el.node[i] == node) { return i; }
    }
@@ -3625,7 +3789,7 @@ void NCMesh::FindNeighbors(int elem, Array<int> &neighbors,
       }
       else
       {
-         for (int i = 0; i < 8 && el.child[i] >= 0; i++)
+         for (int i = 0; i < MaxElemChildren && el.child[i] >= 0; i++)
          {
             stack.Append(el.child[i]);
          }
@@ -3737,8 +3901,8 @@ int NCMesh::GetVertexRootCoord(int elem, RefCoord coord[3]) const
       MFEM_ASSERT(pa.ref_type, "internal error");
 
       int ch = 0;
-      while (ch < 8 && pa.child[ch] != elem) { ch++; }
-      MFEM_ASSERT(ch < 8, "internal error");
+      while (ch < MaxElemChildren && pa.child[ch] != elem) { ch++; }
+      MFEM_ASSERT(ch < MaxElemChildren, "internal error");
 
       MFEM_ASSERT(geom_parent[el.Geom()], "unsupported geometry");
       const RefTrf &tr = geom_parent[el.Geom()][(int) pa.ref_type][ch];
@@ -3768,6 +3932,11 @@ static bool RefPointInside(Geometry::Type geom, const RefCoord pt[3])
          return (pt[0] >= 0) && (pt[1] >= 0) && (pt[0] + pt[1] <= T_ONE) &&
                 (pt[2] >= 0) && (pt[2] <= T_ONE);
 
+      case Geometry::PYRAMID:
+         return (pt[0] >= 0) && (pt[1] >= 0) && (pt[2] >= 0.0) &&
+                (pt[0] + pt[2] <= T_ONE) && (pt[1] + pt[2] <= T_ONE) &&
+                (pt[2] <= T_ONE);
+
       default:
          MFEM_ABORT("unsupported geometry");
          return false;
@@ -3785,7 +3954,7 @@ void NCMesh::CollectIncidentElements(int elem, const RefCoord coord[3],
    }
 
    RefCoord tcoord[3];
-   for (int ch = 0; ch < 8 && el.child[ch] >= 0; ch++)
+   for (int ch = 0; ch < MaxElemChildren && el.child[ch] >= 0; ch++)
    {
       const RefTrf &tr = geom_child[el.Geom()][(int) el.ref_type][ch];
       tr.Apply(coord, tcoord);
@@ -3856,6 +4025,10 @@ NCMesh::PointMatrix NCMesh::pm_prism_identity(
    Point(0, 0, 0), Point(1, 0, 0), Point(0, 1, 0),
    Point(0, 0, 1), Point(1, 0, 1), Point(0, 1, 1)
 );
+NCMesh::PointMatrix NCMesh::pm_pyramid_identity(
+   Point(0, 0, 0), Point(1, 0, 0), Point(1, 1, 0),
+   Point(0, 1, 0), Point(0, 0, 1)
+);
 NCMesh::PointMatrix NCMesh::pm_hex_identity(
    Point(0, 0, 0), Point(1, 0, 0), Point(1, 1, 0), Point(0, 1, 0),
    Point(0, 0, 1), Point(1, 0, 1), Point(1, 1, 1), Point(0, 1, 1)
@@ -3870,6 +4043,7 @@ const NCMesh::PointMatrix& NCMesh::GetGeomIdentity(Geometry::Type geom)
       case Geometry::SQUARE:      return pm_quad_identity;
       case Geometry::TETRAHEDRON: return pm_tet_identity;
       case Geometry::PRISM:       return pm_prism_identity;
+      case Geometry::PYRAMID:     return pm_pyramid_identity;
       case Geometry::CUBE:        return pm_hex_identity;
       default:
          MFEM_ABORT("unsupported geometry " << geom);
@@ -4174,6 +4348,55 @@ void NCMesh::GetPointMatrix(Geometry::Type geom, const char* ref_path,
             }
          }
       }
+      else if (geom == Geometry::PYRAMID)
+      {
+         Point mid01(pm(0), pm(1)), mid23(pm(2), pm(3));
+         Point mid03(pm(0), pm(3)), mid12(pm(1), pm(2));
+         Point mid04(pm(0), pm(4)), mid14(pm(1), pm(4));
+         Point mid24(pm(2), pm(4)), mid34(pm(3), pm(4));
+         Point midf0(mid23, mid12, mid01, mid03);
+
+         if (child == 0)   //Pyramid
+         {
+            pm = PointMatrix(pm(0), mid01, midf0, mid03, mid04);
+         }
+         if (child == 1)   //Tet
+         {
+            pm = PointMatrix(mid01, midf0, mid04, mid14);
+         }
+         if (child == 2)   //Pyramid
+         {
+            pm = PointMatrix(mid01, pm(1), mid12, midf0, mid14);
+         }
+         if (child == 3)   //Tet
+         {
+            pm = PointMatrix(midf0, mid14, mid12, mid24);
+         }
+         if (child == 4)   //Pyramid
+         {
+            pm = PointMatrix(midf0, mid12, pm(2), mid23, mid24);
+         }
+         if (child == 5)   //Tet
+         {
+            pm = PointMatrix(midf0, mid23, mid34, mid24);
+         }
+         if (child == 6)   //Pyramid
+         {
+            pm = PointMatrix(mid03, midf0, mid23, pm(3), mid34);
+         }
+         if (child == 7)   //Tet
+         {
+            pm = PointMatrix(mid03, mid04, midf0, mid34);
+         }
+         if (child == 8)   //Pyramid
+         {
+            pm = PointMatrix(mid24, mid14, mid04, mid34, midf0);
+         }
+         if (child == 9)   //Pyramid
+         {
+            pm = PointMatrix(mid04, mid14, mid24, mid34, pm(4));
+         }
+      }
       else if (geom == Geometry::TETRAHEDRON)
       {
          Point mid01(pm(0), pm(1)), mid12(pm(1), pm(2)), mid02(pm(2), pm(0));
@@ -4333,7 +4556,7 @@ void NCMesh::TraverseRefinements(int elem, int coarse_index,
       ref_path.push_back(el.ref_type);
       ref_path.push_back(0);
 
-      for (int i = 0; i < 8; i++)
+      for (int i = 0; i < MaxElemChildren; i++)
       {
          if (el.child[i] >= 0)
          {
@@ -5091,14 +5314,14 @@ void NCMesh::CountSplits(int elem, int splits[3]) const
    const int* node = el.node;
    GeomInfo& gi = GI[el.Geom()];
 
-   int elevel[12];
+   int elevel[MaxElemEdges];
    for (int i = 0; i < gi.ne; i++)
    {
       const int* ev = gi.edges[i];
       elevel[i] = EdgeSplitLevel(node[ev[0]], node[ev[1]]);
    }
 
-   int flevel[6][2];
+   int flevel[MaxElemFaces][2];
    if (Dim >= 3)
    {
       for (int i = 0; i < gi.nf; i++)
@@ -5142,6 +5365,18 @@ void NCMesh::CountSplits(int elem, int splits[3]) const
       splits[2] = max6(flevel[2][1], flevel[3][1], flevel[4][1],
                        elevel[6], elevel[7], elevel[8]);
    }
+   else if (el.Geom() == Geometry::PYRAMID)
+   {
+      splits[0] = std::max(
+                     max6(flevel[0][0], flevel[1][0], 0,
+                          flevel[2][0], flevel[3][0], flevel[4][0]),
+                     max8(elevel[0], elevel[1], elevel[2],
+                          elevel[3], elevel[4], elevel[5],
+                          elevel[6], elevel[7]));
+
+      splits[1] = splits[0];
+      splits[2] = splits[0];
+   }
    else if (el.Geom() == Geometry::TETRAHEDRON)
    {
       splits[0] = std::max(
@@ -5404,7 +5639,8 @@ void NCMesh::Print(std::ostream &os) const
       "# SQUARE      = 3\n"
       "# TETRAHEDRON = 4\n"
       "# CUBE        = 5\n"
-      "# PRISM       = 6\n";
+      "# PRISM       = 6\n"
+      "# PYRAMID     = 7\n";
 
    os << "\ndimension\n" << Dim << "\n";
 
@@ -5425,7 +5661,7 @@ void NCMesh::Print(std::ostream &os) const
       if (el.parent == -2) { os << "-1\n"; continue; } // unused element
 
       os << int(el.geom) << " " << int(el.ref_type);
-      for (int j = 0; j < 8 && el.node[j] >= 0; j++)
+      for (int j = 0; j < MaxElemNodes && el.node[j] >= 0; j++)
       {
          os << " " << el.node[j];
       }
@@ -5482,7 +5718,7 @@ void NCMesh::InitRootElements()
       Element &el = elements[i];
       if (el.ref_type)
       {
-         for (int j = 0; j < 8 && el.child[j] >= 0; j++)
+         for (int j = 0; j < MaxElemChildren && el.child[j] >= 0; j++)
          {
             int child = el.child[j];
             MFEM_VERIFY(child < elements.Size(), "invalid mesh file: "
@@ -5704,7 +5940,7 @@ void NCMesh::CopyElements(int elem,
    Element &el = elements[elem];
    if (el.ref_type)
    {
-      for (int i = 0; i < 8 && el.child[i] >= 0; i++)
+      for (int i = 0; i < MaxElemChildren && el.child[i] >= 0; i++)
       {
          int old_id = el.child[i];
          // here we know 'free_element_ids' is empty
@@ -6073,7 +6309,7 @@ void NCMesh::DebugLeafOrder(std::ostream &os) const
       {
          double sum = 0.0;
          int count = 0;
-         for (int k = 0; k < 8; k++)
+         for (int k = 0; k < MaxElemNodes; k++)
          {
             if (elem->node[k] >= 0)
             {

mesh/ncmesh.hpp

@@ -153,6 +153,11 @@ public:
    int GetNFaces() const { return NFaces; }
    virtual int GetNGhostElements() const { return 0; }
 
+   /** NCMesh can change the vertex ordering after refinement, coarsening, or on creation
+    *  of the NCMesh object.  After update operation the Vertex ID Map contains the remapping
+    *  information. */
+   const Array<int> &GetVertexIDMap() {return vertex_nodeId;}
+
    /** Perform the given batch of refinements. Please note that in the presence
        of anisotropic splits additional refinements may be necessary to keep
        the mesh consistent. However, the function always performs at least the
@@ -282,7 +287,6 @@ public:
       }
    }
 
-
    // coarse/fine transforms
 
    /** Remember the current layer of leaf elements before the mesh is refined.
@@ -421,6 +425,15 @@ protected: // implementation
    int Geoms; ///< bit mask of element geometries present, see InitGeomFlags()
    bool Legacy; ///< true if the mesh was loaded from the legacy v1.1 format
 
+   static const int MaxElemNodes =
+      8;       ///< Size of the buffers for nodes of an element
+   static const int MaxElemEdges =
+      12;      ///< Size of the buffers for faces of an element
+   static const int MaxElemFaces =
+      6;       ///< Size of the buffers for faces of an element
+   static const int MaxElemChildren =
+      10;   ///< Size of the buffers for children of an element
+
    /** A Node can hold a vertex, an edge, or both. Elements directly point to
        their corner nodes, but edge nodes also exist and can be accessed using
        a hash-table given their two end-point node IDs. All nodes can be
@@ -482,8 +495,8 @@ protected: // implementation
       int attribute;
       union
       {
-         int node[8];  ///< element corners (if ref_type == 0)
-         int child[8]; ///< 2-8 children (if ref_type != 0)
+         int node[MaxElemNodes];  ///< element corners (if ref_type == 0)
+         int child[MaxElemChildren]; ///< 2-10 children (if ref_type != 0)
       };
       int parent; ///< parent element, -1 if this is a root element, -2 if free'd
 
@@ -586,6 +599,9 @@ protected: // implementation
    /// Return true if the mesh contains prism elements.
    bool HavePrisms() const { return Geoms & (1 << Geometry::PRISM); }
 
+   /// Return true if the mesh contains pyramid elements.
+   bool HavePyramids() const { return Geoms & (1 << Geometry::PYRAMID); }
+
    /// Return true if the mesh contains tetrahedral elements.
    bool HaveTets() const   { return Geoms & (1 << Geometry::TETRAHEDRON); }
 
@@ -642,6 +658,10 @@ protected: // implementation
    int NewTetrahedron(int n0, int n1, int n2, int n3, int attr,
                       int fattr0, int fattr1, int fattr2, int fattr3);
 
+   int NewPyramid(int n0, int n1, int n2, int n3, int n4, int attr,
+                  int fattr0, int fattr1, int fattr2, int fattr3,
+                  int fattr4);
+
    int NewQuadrilateral(int n0, int n1, int n2, int n3, int attr,
                         int eattr0, int eattr1, int eattr2, int eattr3);
 
@@ -665,6 +685,9 @@ protected: // implementation
    void CheckAnisoPrism(int vn1, int vn2, int vn3, int vn4,
                         const Refinement *refs, int nref);
 
+   void CheckAnisoPyramid(int vn1, int vn2, int vn3, int vn4,
+                          const Refinement *refs, int nref);
+
    void CheckAnisoFace(int vn1, int vn2, int vn3, int vn4,
                        int mid12, int mid34, int level = 0);
 
@@ -845,7 +868,7 @@ protected: // implementation
    struct PointMatrix
    {
       int np;
-      Point points[8];
+      Point points[MaxElemNodes];
 
       PointMatrix() : np(0) {}
 
@@ -858,6 +881,13 @@ protected: // implementation
       PointMatrix(const Point& p0, const Point& p1, const Point& p2, const Point& p3)
       { np = 4; points[0] = p0; points[1] = p1; points[2] = p2; points[3] = p3; }
 
+      PointMatrix(const Point& p0, const Point& p1, const Point& p2,
+                  const Point& p3, const Point& p4)
+      {
+         np = 5;
+         points[0] = p0; points[1] = p1; points[2] = p2;
+         points[3] = p3; points[4] = p4;
+      }
       PointMatrix(const Point& p0, const Point& p1, const Point& p2,
                   const Point& p3, const Point& p4, const Point& p5)
       {
@@ -887,6 +917,7 @@ protected: // implementation
    static PointMatrix pm_quad_identity;
    static PointMatrix pm_tet_identity;
    static PointMatrix pm_prism_identity;
+   static PointMatrix pm_pyramid_identity;
    static PointMatrix pm_hex_identity;
 
    static const PointMatrix& GetGeomIdentity(Geometry::Type geom);
@@ -976,9 +1007,9 @@ protected: // implementation
    struct GeomInfo
    {
       int nv, ne, nf;   // number of: vertices, edges, faces
-      int edges[12][2]; // edge vertices (up to 12 edges)
-      int faces[6][4];  // face vertices (up to 6 faces)
-      int nfv[6];       // number of face vertices
+      int edges[MaxElemEdges][2]; // edge vertices (up to 12 edges)
+      int faces[MaxElemFaces][4];  // face vertices (up to 6 faces)
+      int nfv[MaxElemFaces];       // number of face vertices
 
       bool initialized;
       GeomInfo() : initialized(false) {}

mesh/ncmesh_tables.hpp

@@ -49,6 +49,16 @@ const int prism_deref_table[7][6 + 5] =
    { 0, 1, 2, 4, 5, 6, /**/  0, 5, 0, 5, 0 }  // 7 - iso
 };
 
+const int pyramid_deref_table[7][5 + 5] =
+{
+   {-1,-1,-1,-1,-1, /**/ -1,-1,-1,-1,-1 }, // 1
+   {-1,-1,-1,-1,-1, /**/ -1,-1,-1,-1,-1 }, // 2
+   {-1,-1,-1,-1,-1, /**/ -1,-1,-1,-1,-1 }, // 3
+   {-1,-1,-1,-1,-1, /**/ -1,-1,-1,-1,-1 }, // 4
+   {-1,-1,-1,-1,-1, /**/ -1,-1,-1,-1,-1 }, // 5
+   {-1,-1,-1,-1,-1, /**/ -1,-1,-1,-1,-1 }, // 6
+   { 0, 2, 4, 6, 9, /**/  0, 9, 9, 9, 9 }  // 7 - iso
+};
 
 // child ordering tables
 
@@ -140,8 +150,17 @@ const RefCoord prism_corners[6][3] =
    {    0, T_ONE, T_ONE}
 };
 
+const RefCoord pyramid_corners[5][3] =
+{
+   {    0,     0,     0},
+   {T_ONE,     0,     0},
+   {T_ONE, T_ONE,     0},
+   {    0, T_ONE,     0},
+   {    0,     0, T_ONE}
+};
+
 typedef RefCoord RefPoint[3];
-const RefPoint* geom_corners[7] =
+const RefPoint* geom_corners[8] =
 {
    NULL, // point
    NULL, // segment
@@ -149,7 +168,8 @@ const RefPoint* geom_corners[7] =
    quad_corners,
    NULL, // tetrahedron
    hex_corners,
-   prism_corners
+   prism_corners,
+   pyramid_corners
 };
 
 // reference domain transform: 3 scales, 3 translations

miniapps/electromagnetics/bfield_remhos.hpp

@@ -0,0 +1,448 @@
+// Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights
+// reserved. See files LICENSE and NOTICE for details.
+//
+// This file is part of CEED, a collection of benchmarks, miniapps, software
+// libraries and APIs for efficient high-order finite element and spectral
+// element discretizations for exascale applications. For more information and
+// source code availability see http://github.com/ceed.
+//
+// The CEED research is supported by the Exascale Computing Project 17-SC-20-SC,
+// a collaborative effort of two U.S. Department of Energy organizations (Office
+// of Science and the National Nuclear Security Administration) responsible for
+// the planning and preparation of a capable exascale ecosystem, including
+// software, applications, hardware, advanced system engineering and early
+// testbed platforms, in support of the nation's exascale computing imperative.
+
+#ifndef MFEM_BFIELD_REMHOS
+#define MFEM_BFIELD_REMHOS
+
+#include "mfem.hpp"
+#include "general/forall.hpp"
+
+#define EMPTY_ZONE_TOL 1e-12
+
+namespace mfem
+{
+
+namespace electromagnetics
+{
+
+// Class storing information on dofs needed for the low order methods and FCT.
+class DofInfo
+{
+private:
+   // 0 is overlap, see ComputeOverlapBounds().
+   // 1 is sparcity, see ComputeMatrixSparcityBounds().
+   int bounds_type;
+   ParMesh *pmesh;
+   ParFiniteElementSpace &pfes;
+
+   // The min and max bounds are represented as CG functions of the same order
+   // as the solution, thus having 1:1 dof correspondence inside each element.
+   H1_FECollection fec_bounds;
+   ParFiniteElementSpace pfes_bounds;
+   ParGridFunction x_min, x_max;
+
+   // For each DOF on an element boundary, the global index of the DOF on the
+   // opposite site is computed and stored in a list. This is needed for lumping
+   // the flux contributions, as in the paper. Right now it works on 1D meshes,
+   // quad meshes in 2D and 3D meshes of ordered cubes.
+   // NOTE: The mesh is assumed to consist of segments, quads or hexes.
+   // NOTE: This approach will not work for meshes with hanging nodes.
+   void FillNeighborDofs();
+
+   // A list is filled to later access the correct element-global indices given
+   // the subcell number and subcell index.
+   // NOTE: The mesh is assumed to consist of segments, quads or hexes.
+   void FillSubcell2CellDof();
+
+   // Computes the admissible interval of values for each DG dof from the values
+   // of all elements that feature the dof at its physical location.
+   // A given DOF gets bounds from the elements it touches (in Gauss-Lobatto
+   // sense, i.e., a face dof touches two elements, vertex dofs can touch many).
+   void ComputeOverlapBounds(const Vector &el_min, const Vector &el_max,
+                             Vector &dof_min, Vector &dof_max,
+                             Array<bool> *active_el = NULL);
+
+   // A given DOF gets bounds from its own element and its face-neighbors.
+   void ComputeMatrixSparsityBounds(const Vector &el_min, const Vector &el_max,
+                                    Vector &dof_min, Vector &dof_max,
+                                    Array<bool> *active_el = NULL);
+
+public:
+   Vector xi_min, xi_max; // min/max values for each dof
+   Vector xe_min, xe_max; // min/max values for each element
+
+   DenseMatrix BdrDofs, Sub2Ind;
+   DenseTensor NbrDof;
+
+   int numBdrs, numFaceDofs, numSubcells, numDofsSubcell;
+
+   DofInfo(ParFiniteElementSpace &pfes_sltn, int btype = 0);
+
+   // Computes the admissible interval of values for each DG dof from the values
+   // of all elements that feature the dof at its physical location.
+   void ComputeBounds(const Vector &el_min, const Vector &el_max,
+                      Vector &dof_min, Vector &dof_max,
+                      Array<bool> *active_el = NULL)
+   {
+      if (bounds_type == 0)
+      {
+         ComputeOverlapBounds(el_min, el_max, dof_min, dof_max, active_el);
+      }
+      else if (bounds_type == 1)
+      {
+         ComputeMatrixSparsityBounds(el_min, el_max,
+                                     dof_min, dof_max, active_el);
+      }
+      else { MFEM_ABORT("Wrong option for bounds computation."); }
+   }
+
+   // Computes the min and max values of u over each element.
+   void ComputeElementsMinMax(const Vector &u,
+                              Vector &u_min, Vector &u_max,
+                              Array<bool> *active_el,
+                              Array<bool> *active_dof) const;
+};
+
+struct LowOrderMethod
+{
+   bool subcell_scheme;
+   FiniteElementSpace *SubFes0, *SubFes1;
+   Array <int> smap;
+   SparseMatrix D;
+   ParBilinearForm* pk;
+   VectorCoefficient* coef;
+   VectorCoefficient* subcellCoeff;
+   const IntegrationRule* irF;
+   BilinearFormIntegrator* VolumeTerms;
+};
+
+class SmoothnessIndicator
+{
+private:
+   const int type;
+   const double param;
+   H1_FECollection fec_sub;
+   ParFiniteElementSpace pfes_CG_sub;
+   ParFiniteElementSpace &pfes_DG;
+   SparseMatrix Mmat, LaplaceOp, *MassMixed;
+   BilinearFormIntegrator *MassInt;
+   Vector lumpedMH1;
+   DenseMatrix ShapeEval;
+
+   void ComputeVariationalMatrix(DofInfo &dof_info);
+   void ApproximateLaplacian(const Vector &x, ParGridFunction &y);
+   void ComputeFromSparsity(const SparseMatrix &K, const ParGridFunction &x,
+                            Vector &x_min, Vector &x_max);
+
+public:
+   SmoothnessIndicator(int type_id,
+                       ParMesh &subcell_mesh,
+                       ParFiniteElementSpace &pfes_DG_,
+                       ParGridFunction &u,
+                       DofInfo &dof_info);
+   ~SmoothnessIndicator();
+
+   void ComputeSmoothnessIndicator(const Vector &u, ParGridFunction &si_vals_u);
+   void UpdateBounds(int dof_id, double u_HO,
+                     const ParGridFunction &si_vals,
+                     double &u_min, double &u_max);
+
+   Vector DG2CG;
+};
+
+class Assembly
+{
+private:
+   const int exec_mode;
+   const GridFunction &inflow_gf;
+   mutable ParGridFunction x_gf;
+   BilinearFormIntegrator *VolumeTerms;
+   FiniteElementSpace *fes, *SubFes0, *SubFes1;
+   Mesh *subcell_mesh;
+
+public:
+   Assembly(DofInfo &_dofs, LowOrderMethod &inlom, const GridFunction &inflow,
+            ParFiniteElementSpace &pfes, ParMesh *submesh, int mode);
+
+   // Auxiliary member variables that need to be accessed during time-stepping.
+   DofInfo &dofs;
+
+   LowOrderMethod &lom;
+   // Data structures storing Galerkin contributions. These are updated for
+   // remap but remain constant for transport.
+   // bdrInt - eq (32).
+   // SubcellWeights - above eq (49).
+   DenseTensor bdrInt, SubcellWeights;
+
+   void ComputeFluxTerms(const int e_id, const int BdrID,
+                         FaceElementTransformations *Trans,
+                         LowOrderMethod &lom);
+
+   void ComputeSubcellWeights(const int k, const int m);
+
+   void LinearFluxLumping(const int k, const int nd,
+                          const int BdrID, const Vector &x,
+                          Vector &y, const Vector &x_nd,
+                          const Vector &alpha) const;
+   void NonlinFluxLumping(const int k, const int nd,
+                          const int BdrID, const Vector &x,
+                          Vector &y, const Vector &x_nd,
+                          const Vector &alpha) const;
+
+   const FiniteElementSpace * GetFes() {return fes;}
+
+   int GetExecMode() const { return exec_mode;}
+
+   Mesh *GetSubCellMesh() { return subcell_mesh;}
+};
+
+
+// Class for local assembly of M_L M_C^-1 K, where M_L and M_C are the lumped
+// and consistent mass matrices and K is the convection matrix. The spaces are
+// assumed to be L2 conforming.
+class PrecondConvectionIntegrator: public BilinearFormIntegrator
+{
+private:
+#ifndef MFEM_THREAD_SAFE
+   DenseMatrix dshape, adjJ, Q_ir;
+   Vector shape, vec2, BdFidxT;
+#endif
+   VectorCoefficient &Q;
+   double alpha;
+
+public:
+   PrecondConvectionIntegrator(VectorCoefficient &q, double a = 1.0)
+      : Q(q) { alpha = a; }
+   virtual void AssembleElementMatrix(const FiniteElement &,
+                                      ElementTransformation &,
+                                      DenseMatrix &);
+};
+
+// alpha (q . grad u, v)
+class MixedConvectionIntegrator : public BilinearFormIntegrator
+{
+private:
+#ifndef MFEM_THREAD_SAFE
+   DenseMatrix dshape, adjJ, Q_ir;
+   Vector shape, vec2, BdFidxT;
+#endif
+   VectorCoefficient &Q;
+   double alpha;
+
+public:
+   MixedConvectionIntegrator(VectorCoefficient &q, double a = 1.0)
+      : Q(q) { alpha = a; }
+   virtual void AssembleElementMatrix2(const FiniteElement &tr_el,
+                                       const FiniteElement &te_el,
+                                       ElementTransformation &Trans,
+                                       DenseMatrix &elmat);
+};
+
+
+// Low-Order Solver.
+class LOSolver
+{
+protected:
+   ParFiniteElementSpace &pfes;
+   double dt = -1.0; // usually not known at creation, updated later.
+
+public:
+   LOSolver(ParFiniteElementSpace &space) : pfes(space) { }
+
+   virtual ~LOSolver() { }
+
+   virtual void UpdateTimeStep(double dt_new) { dt = dt_new; }
+
+   virtual void CalcLOSolution(const Vector &u, Vector &du) const = 0;
+};
+
+class DiscreteUpwind : public LOSolver
+{
+protected:
+   const SparseMatrix &K;
+   mutable SparseMatrix D;
+   const Array<int> &K_smap;
+   const Vector &M_lumped;
+   Assembly &assembly;
+   const bool update_D;
+
+   void ComputeDiscreteUpwindMatrix() const;
+
+public:
+   DiscreteUpwind(ParFiniteElementSpace &space, const SparseMatrix &adv,
+                  const Array<int> &adv_smap, const Vector &Mlump,
+                  Assembly &asmbly, bool updateD);
+
+   virtual void CalcLOSolution(const Vector &u, Vector &du) const;
+};
+
+// High-Order Solver.
+// Conserve mass / provide high-order convergence / may violate the bounds.
+class HOSolver
+{
+protected:
+   ParFiniteElementSpace &pfes;
+
+public:
+   HOSolver(ParFiniteElementSpace &space) : pfes(space) { }
+
+   virtual ~HOSolver() { }
+
+   virtual void CalcHOSolution(const Vector &u, Vector &du) const = 0;
+};
+
+class LocalInverseHOSolver : public HOSolver
+{
+protected:
+   ParBilinearForm &M, &K;
+
+public:
+   LocalInverseHOSolver(ParFiniteElementSpace &space,
+                        ParBilinearForm &Mbf, ParBilinearForm &Kbf);
+
+   virtual void CalcHOSolution(const Vector &u, Vector &du) const;
+};
+
+// Monotone, High-order, Conservative Solver.
+class FCTSolver
+{
+protected:
+   ParFiniteElementSpace &pfes;
+   SmoothnessIndicator *smth_indicator;
+   double dt;
+   const bool needs_LO_input_for_products;
+
+   // Computes a compatible slope (piecewise constan = mass_us / mass_u).
+   // It could also update s_min and s_max, if required.
+   void CalcCompatibleLOProduct(const ParGridFunction &us,
+                                const Vector &m, const Vector &d_us_HO,
+                                Vector &s_min, Vector &s_max,
+                                const Vector &u_new,
+                                const Array<bool> &active_el,
+                                const Array<bool> &active_dofs,
+                                Vector &d_us_LO_new);
+   void ScaleProductBounds(const Vector &s_min, const Vector &s_max,
+                           const Vector &u_new, const Array<bool> &active_el,
+                           const Array<bool> &active_dofs,
+                           Vector &us_min, Vector &us_max);
+
+public:
+   FCTSolver(ParFiniteElementSpace &space,
+             SmoothnessIndicator *si, double dt_, bool needs_LO_prod)
+      : pfes(space), smth_indicator(si), dt(dt_),
+        needs_LO_input_for_products(needs_LO_prod) { }
+
+   virtual ~FCTSolver() { }
+
+   virtual void UpdateTimeStep(double dt_new) { dt = dt_new; }
+
+   bool NeedsLOProductInput() const { return needs_LO_input_for_products; }
+
+   // Calculate du that satisfies the following:
+   // bounds preservation: u_min_i <= u_i + dt du_i <= u_max_i,
+   // conservation:        sum m_i (u_i + dt du_ho_i) = sum m_i (u_i + dt du_i).
+   // Some methods utilize du_lo as a backup choice, as it satisfies the above.
+   virtual void CalcFCTSolution(const ParGridFunction &u, const Vector &m,
+                                const Vector &du_ho, const Vector &du_lo,
+                                const Vector &u_min, const Vector &u_max,
+                                Vector &du) const = 0;
+
+   // Used in the case of product remap.
+   // Given the input, calculates d_us, so that:
+   // bounds preservation: s_min_i <= (us_i + dt d_us_i) / u_new_i <= s_max_i,
+   // conservation: sum m_i (us_i + dt d_us_HO_i) = sum m_i (us_i + dt d_us_i).
+   virtual void CalcFCTProduct(const ParGridFunction &us, const Vector &m,
+                               const Vector &d_us_HO, const Vector &d_us_LO,
+                               Vector &s_min, Vector &s_max,
+                               const Vector &u_new,
+                               const Array<bool> &active_el,
+                               const Array<bool> &active_dofs, Vector &d_us)
+   {
+      MFEM_ABORT("Product remap is not implemented for the chosen solver");
+   }
+};
+
+class FluxBasedFCT : public FCTSolver
+{
+protected:
+   const SparseMatrix &K, &M;
+   const Array<int> &K_smap;
+
+   // Temporary computation objects.
+   mutable SparseMatrix flux_ij;
+   mutable ParGridFunction gp, gm;
+
+   const int iter_cnt;
+
+   void ComputeFluxMatrix(const ParGridFunction &u, const Vector &du_ho,
+                          SparseMatrix &flux_mat) const;
+   void AddFluxesAtDofs(const SparseMatrix &flux_mat,
+                        Vector &flux_pos, Vector &flux_neg) const;
+   void ComputeFluxCoefficients(const Vector &u, const Vector &du_lo,
+                                const Vector &m, const Vector &u_min, const Vector &u_max,
+                                Vector &coeff_pos, Vector &coeff_neg) const;
+   void UpdateSolutionAndFlux(const Vector &du_lo, const Vector &m,
+                              ParGridFunction &coeff_pos, ParGridFunction &coeff_neg,
+                              SparseMatrix &flux_mat, Vector &du) const;
+
+public:
+   FluxBasedFCT(ParFiniteElementSpace &space,
+                SmoothnessIndicator *si, double delta_t,
+                const SparseMatrix &adv_mat, const Array<int> &adv_smap,
+                const SparseMatrix &mass_mat, int fct_iterations = 1)
+      : FCTSolver(space, si, delta_t, true),
+        K(adv_mat), M(mass_mat), K_smap(adv_smap), flux_ij(adv_mat),
+        gp(&pfes), gm(&pfes), iter_cnt(fct_iterations) { }
+
+   virtual void CalcFCTSolution(const ParGridFunction &u, const Vector &m,
+                                const Vector &du_ho, const Vector &du_lo,
+                                const Vector &u_min, const Vector &u_max,
+                                Vector &du) const;
+
+   virtual void CalcFCTProduct(const ParGridFunction &us, const Vector &m,
+                               const Vector &d_us_HO, const Vector &d_us_LO,
+                               Vector &s_min, Vector &s_max,
+                               const Vector &u_new,
+                               const Array<bool> &active_el,
+                               const Array<bool> &active_dofs, Vector &d_us);
+};
+
+int GetLocalFaceDofIndex(int dim, int loc_face_id, int face_orient,
+                         int face_dof_id, int face_dof1D_cnt);
+
+void ExtractBdrDofs(int p, Geometry::Type gtype, DenseMatrix &dofs);
+
+void GetMinMax(const ParGridFunction &g, double &min, double &max);
+
+// Utility function to build a map to the offset of the symmetric entry in a
+// sparse matrix.
+Array<int> SparseMatrix_Build_smap(const SparseMatrix &A);
+
+// Given a matrix K, matrix D (initialized with same sparsity as K) is computed,
+// such that (K+D)_ij >= 0 for i != j.
+void ComputeDiscreteUpwindingMatrix(const SparseMatrix &K,
+                                    Array<int> smap, SparseMatrix& D);
+
+void VisualizeField(socketstream &sock, const char *vishost, int visport,
+                    ParGridFunction &gf, const char *title,
+                    int x, int y, int w, int h,
+                    const char *keys = NULL, bool vec = false);
+
+void ComputeBoolIndicators(int NE, const Vector &u,
+                           Array<bool> &ind_elem, Array<bool> &ind_dofs);
+
+void ComputeRatio(int NE, const Vector &u_s, const Vector &u,
+                  Vector &s, Array<bool> &bool_el, Array<bool> &bool_dof);
+
+void ZeroOutEmptyDofs(const Array<bool> &ind_elem,
+                      const Array<bool> &ind_dofs, Vector &u);
+
+
+} // namespace electromagnetics
+
+} // namespace mfem
+
+#endif // MFEM_REMHOS_TOOLS

miniapps/electromagnetics/bfieldadvect.cpp

@@ -0,0 +1,205 @@
+#include "bfieldadvect_solver.hpp"
+#include <random>
+
+using namespace std;
+using namespace mfem;
+using namespace mfem::electromagnetics;
+
+void BFieldFunc(const Vector &, Vector&);
+void PeturbBoxMesh(Mesh *mesh, double xlen, double ylen, double zlen,
+                   double window);
+
+int main(int argc, char *argv[])
+{
+   MPI_Session mpi(argc, argv);
+
+   // Parse command-line options.
+   const char *mesh_file = "../../data/toroid-hex.mesh";
+   int order = 1;
+   int serial_ref_levels = 0;
+   int parallel_ref_levels = 0;
+   bool visualization = false;
+   bool visit = true;
+
+   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(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&visit, "-visit", "--visit", "-no-visit",
+                  "--no-visualization",
+                  "Enable or disable VisIt visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (mpi.Root())
+      {
+         args.PrintUsage(cout);
+      }
+      return 1;
+   }
+   if (mpi.Root())
+   {
+      args.PrintOptions(cout);
+   }
+
+
+   Mesh *test_mesh = new Mesh("../../data/ref-pyramid.mesh", 1, 1);
+   std::cout << test_mesh->GetElementVolume(0) << std::endl;
+   Array<Refinement> ref(1);
+   ref[0].ref_type = Refinement::XYZ;
+   ref[0].index = 0;
+   test_mesh->GeneralRefinement(ref, 1);
+   double sum = 0.0;
+   for (int i = 0; i < 10; ++i)
+   {
+      sum += test_mesh->GetElementVolume(i);
+   }
+   std::cout << sum << std::endl;
+   Array<double> elem_error(10);
+   for (int i = 0; i < 10; ++i)
+   {
+      elem_error[i] = 0.0;
+   }
+   test_mesh->DerefineByError(elem_error, 1.0);
+   std::cout << test_mesh->GetElementVolume(0) << std::endl;
+
+   VisItDataCollection visit_dc_test("test", test_mesh);
+   visit_dc_test.SetCycle(0);
+   visit_dc_test.SetTime(0);
+   visit_dc_test.Save();
+
+
+
+   Mesh *mesh_old = new Mesh(20, 20, 5, Element::HEXAHEDRON, false, 4.0, 4.0, 1.0);
+   Mesh *mesh_new = new Mesh(*mesh_old, true);
+   double dl = 1.0/5.0;
+   PeturbBoxMesh(mesh_new, 4.0, 4.0, 1.0, 0.05*dl);
+
+
+   // 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_old->UniformRefinement();
+      mesh_new->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_old(MPI_COMM_WORLD, *mesh_old);
+   ParMesh pmesh_new(MPI_COMM_WORLD, *mesh_new);
+   delete mesh_old;
+   delete mesh_new;
+
+   // Refine this mesh in parallel to increase the resolution.
+   for (int l = 0; l < parallel_ref_levels; l++)
+   {
+      pmesh_old.UniformRefinement();
+      pmesh_new.UniformRefinement();
+   }
+
+   //Set up the pre and post advection fields on the relevant meshes/spaces
+   RT_ParFESpace *HDivFESpaceOld  = new RT_ParFESpace(&pmesh_old,order,
+                                                      pmesh_old.Dimension());
+   RT_ParFESpace *HDivFESpaceNew  = new RT_ParFESpace(&pmesh_new,order,
+                                                      pmesh_new.Dimension());
+   ParGridFunction *b = new ParGridFunction(HDivFESpaceOld);
+   ParGridFunction *b_new = new ParGridFunction(HDivFESpaceNew);
+   ParGridFunction *b_new_exact = new ParGridFunction(HDivFESpaceNew);
+
+   //Set the initial B value
+   *b = 0.0;
+   *b_new = 0.0;
+   *b_new_exact = 0.0;
+   VectorFunctionCoefficient BFieldCoef(3,BFieldFunc);
+   b->ProjectCoefficient(BFieldCoef);
+   b_new_exact->ProjectCoefficient(BFieldCoef);
+
+   BFieldAdvector advector(&pmesh_old, &pmesh_new, 1);
+   advector.Advect(b, b_new);
+   ParGridFunction *b_recon = advector.GetReconstructedB();
+   ParGridFunction *curl_b = advector.GetCurlB();
+   ParGridFunction *a = advector.GetA();
+   ParGridFunction *a_new = advector.GetANew();
+
+   Vector diff_b(*b_new_exact);
+   diff_b -= *b_new;    //diff = b_new_exact - b_new
+   std::cout << "Vector diff in B field on the new mesh:  " << diff_b.Normlinf() <<
+             std::endl;
+
+   Vector diff_a(*a);
+   diff_a -= *a_new;    //diff = b_new_exact - b_new
+   std::cout << "Vector diff in A field on the new mesh:  " << diff_a.Normlinf() <<
+             ", " << a->Normlinf() << std::endl;
+
+   // Handle the visit visualization
+   if (visit)
+   {
+      VisItDataCollection visit_dc_old("bfa-old", &pmesh_old);
+      visit_dc_old.RegisterField("B", b);
+      visit_dc_old.RegisterField("Curl_B", curl_b);
+      visit_dc_old.RegisterField("A", a);
+      visit_dc_old.RegisterField("B_recon", b_recon);
+      visit_dc_old.SetCycle(0);
+      visit_dc_old.SetTime(0);
+      visit_dc_old.Save();
+
+      VisItDataCollection visit_dc_new("bfa-new", &pmesh_new);
+      visit_dc_new.RegisterField("A_new", a_new);
+      visit_dc_new.RegisterField("B_new", b_new);
+      visit_dc_new.SetCycle(0);
+      visit_dc_new.SetTime(0);
+      visit_dc_new.Save();
+   }
+
+}
+
+
+void BFieldFunc(const Vector &x, Vector &B)
+{
+   B.SetSize(3);
+   B[0] =  x[1] - 2.0;
+   B[1] = -(x[0] - 2.0);
+   B[2] =  0.0;
+}
+
+
+void PeturbBoxMesh(Mesh *mesh, double xlen, double ylen, double zlen,
+                   double window)
+{
+   Vector displacements(3*mesh->GetNV());
+   displacements = 0.0;
+   std::random_device
+   rd;  // Will be used to obtain a seed for the random number engine
+   std::mt19937 gen(rd()); // Standard mersenne_twister_engine seeded with rd()
+   std::uniform_real_distribution<> r(-0.5*window, 0.5*window);
+   for (int vi = 0; vi < mesh->GetNV(); ++vi)
+   {
+      double *v = mesh->GetVertex(vi);
+      if (fabs(v[0]) > 1e-6 && fabs(v[0] - xlen) > 1e-6)
+      {
+         displacements[3*vi+0] = r(gen);
+      }
+
+      if (fabs(v[1]) > 1e-6 && fabs(v[1] - ylen) > 1e-6)
+      {
+         displacements[3*vi+1] = r(gen);
+      }
+
+      if (fabs(v[2]) > 1e-6 && fabs(v[2] - zlen) > 1e-6)
+      {
+         displacements[3*vi+2] = r(gen);
+      }
+   }
+   std::cout << "Displacement norm:  " << displacements.Norml2() << std::endl;
+   mesh->MoveVertices(displacements);
+}

miniapps/electromagnetics/bfieldadvect_solver.cpp

@@ -0,0 +1,866 @@
+// Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "bfieldadvect_solver.hpp"
+
+
+//notes
+//Look at DivergenceFreeProjector to pull th divergence out of a vecrot in Hcurl
+// Matrix-Vector Multiplication AddMult(x,y,val):  y = y + val*A*x, Mult(x, y): y = A*x
+
+using namespace std;
+
+namespace mfem
+{
+
+namespace electromagnetics
+{
+
+BFieldAdvector::BFieldAdvector(ParMesh *pmesh_old, ParMesh *pmesh_new,
+                               int order_) :
+   order(order_),
+   pmeshOld(nullptr),
+   pmeshNew(nullptr),
+   H1FESpaceOld(nullptr),
+   HCurlFESpaceOld(nullptr),
+   HDivFESpaceOld(nullptr),
+   L2FESpaceOld(nullptr),
+   H1FESpaceNew(nullptr),
+   HCurlFESpaceNew(nullptr),
+   HDivFESpaceNew(nullptr),
+   L2FESpaceNew(nullptr),
+   grad(nullptr),
+   curl_old(nullptr),
+   curl_new(nullptr),
+   weakCurl(nullptr),
+   WC(nullptr),
+   m1(nullptr),
+   curlCurl(nullptr),
+   divFreeProj(nullptr),
+   a(nullptr),
+   a_new(nullptr),
+   curl_b(nullptr),
+   clean_curl_b(nullptr),
+   recon_b(nullptr)
+{
+   myComm = pmesh_old->GetComm();
+   SetMesh(pmesh_old, pmesh_new);
+}
+
+
+void BFieldAdvector::SetMesh(ParMesh *pmesh_old, ParMesh *pmesh_new)
+{
+   CleanInternals();
+
+   pmeshOld = pmesh_old;
+   pmeshOld->EnsureNodes();
+   pmeshNew = pmesh_new;
+   pmeshNew->EnsureNodes();
+
+   //Set up the various spaces on the meshes
+   H1FESpaceOld    = new H1_ParFESpace(pmesh_old,order,pmesh_old->Dimension());
+   HCurlFESpaceOld = new ND_ParFESpace(pmesh_old,order,pmesh_old->Dimension());
+   HDivFESpaceOld  = new RT_ParFESpace(pmesh_old,order,pmesh_old->Dimension());
+   L2FESpaceOld    = new L2_ParFESpace(pmesh_old,order,pmesh_old->Dimension());
+   H1FESpaceNew    = new H1_ParFESpace(pmesh_new,order,pmesh_new->Dimension());
+   HCurlFESpaceNew = new ND_ParFESpace(pmesh_new,order,pmesh_new->Dimension());
+   HDivFESpaceNew  = new RT_ParFESpace(pmesh_new,order,pmesh_new->Dimension());
+   L2FESpaceNew    = new L2_ParFESpace(pmesh_new,order,pmesh_new->Dimension());
+
+   //Discrete Differential Operators
+   grad = new ParDiscreteGradOperator(H1FESpaceOld, HCurlFESpaceOld);
+   grad->Assemble();
+   grad->Finalize();
+   curl_old = new ParDiscreteCurlOperator(HCurlFESpaceOld, HDivFESpaceOld);
+   curl_old->Assemble();
+   curl_old->Finalize();
+   curl_new = new ParDiscreteCurlOperator(HCurlFESpaceOld, HDivFESpaceOld);
+   curl_new->Assemble();
+   curl_new->Finalize();
+
+   //Weak curl operator for taking the curl of B living in Hdiv
+   ConstantCoefficient oneCoef(1.0);
+   weakCurl = new ParMixedBilinearForm(HDivFESpaceOld, HCurlFESpaceOld);
+   weakCurl->AddDomainIntegrator(new VectorFECurlIntegrator(oneCoef));
+   weakCurl->Assemble();
+   weakCurl->Finalize();
+   WC = weakCurl->ParallelAssemble();
+
+   m1 = new ParBilinearForm(HCurlFESpaceOld);
+   m1->AddDomainIntegrator(new VectorFEMassIntegrator(oneCoef));
+   m1->Assemble();
+   m1->Finalize();
+
+   //CurlCurl operator
+   curlCurl  = new ParBilinearForm(HCurlFESpaceOld);
+   curlCurl->AddDomainIntegrator(new CurlCurlIntegrator(oneCoef));
+   curlCurl->Assemble();
+   curlCurl->Finalize();
+
+   //Projector to clean the divergence out of vectors in Hcurl
+   int irOrder = H1FESpaceOld->GetElementTransformation(0)->OrderW()+ 2 * order;
+   divFreeProj = new DivergenceFreeProjector(*H1FESpaceOld, *HCurlFESpaceOld,
+                                             irOrder, NULL, NULL, grad);
+
+   // Build internal grid functions on the spaces
+   a  = new ParGridFunction(
+      HCurlFESpaceOld);            //Vector potential A in HCurl
+   a_new = new ParGridFunction(
+      HCurlFESpaceNew);         //Vector potential A in Hcurl on the new mesh
+   curl_b = new ParGridFunction(
+      HCurlFESpaceOld);        //curl B in Hcurl from the weak curl
+   clean_curl_b = new ParGridFunction(HCurlFESpaceOld);  //B in Hcurl
+   recon_b = new ParGridFunction(HDivFESpaceOld);        //Reconstructed B from A
+}
+
+
+void BFieldAdvector::CleanInternals()
+{
+   if (H1FESpaceOld != nullptr) { delete H1FESpaceOld; }
+   if (HCurlFESpaceOld != nullptr) { delete HCurlFESpaceOld; }
+   if (HDivFESpaceOld != nullptr) { delete HDivFESpaceOld; }
+   if (L2FESpaceOld != nullptr) { delete L2FESpaceOld; }
+   if (H1FESpaceNew != nullptr) { delete H1FESpaceNew; }
+   if (HCurlFESpaceNew != nullptr) { delete HCurlFESpaceNew; }
+   if (HDivFESpaceNew != nullptr) { delete HDivFESpaceNew; }
+   if (L2FESpaceNew != nullptr) { delete L2FESpaceNew; }
+
+   if (grad != nullptr) { delete grad; }
+   if (curl_old != nullptr) { delete curl_old; }
+   if (curl_new != nullptr) { delete curl_new; }
+
+   if (weakCurl != nullptr) { delete weakCurl; }
+   if (divFreeProj != nullptr) { delete divFreeProj; }
+   if (curlCurl != nullptr) { delete curlCurl; }
+
+   if (a != nullptr) { delete a; }
+   if (a_new != nullptr) { delete a_new; }
+   if (curl_b != nullptr) { delete curl_b; }
+   if (clean_curl_b != nullptr) { delete clean_curl_b; }
+   if (recon_b != nullptr) { delete recon_b; }
+}
+
+
+void BFieldAdvector::Advect(ParGridFunction* b_old, ParGridFunction* b_new)
+{
+   ComputeA(b_old);
+   FindPtsInterpolateToTargetMesh(a, a_new, 1);
+   curl_new->Mult(*a_new, *b_new);
+}
+
+//Solve Curl Curl A = Curl B for A using AMS
+void BFieldAdvector::ComputeA(ParGridFunction* b)
+{
+   Array<int> ess_bdr;
+   ess_bdr.SetSize(pmeshOld->bdr_attributes.Max());
+   ess_bdr = 0;   // All outer surfaces
+   //ess_bdr[0] = 1;
+   Array<int> ess_bdr_tdofs;
+   HCurlFESpaceOld->GetEssentialTrueDofs(ess_bdr, ess_bdr_tdofs);
+
+   //Set up a linear form with a curl operator on B
+   VectorGridFunctionCoefficient b_coef(b);
+   ParLinearForm rhs(HCurlFESpaceOld);
+   rhs.AddDomainIntegrator(new VectorFEDomainLFCurlIntegrator(b_coef));
+   rhs.Assemble();
+
+   // Apply Dirichlet BCs to matrix and right hand side and otherwise
+   // prepare the linear system
+   HypreParMatrix M;
+   Vector A, RHS;
+   *a = 0;
+
+   //curlCurl->FormLinearSystem(ess_bdr_tdofs, *a, *clean_curl_b, M, A, RHS);
+   curlCurl->FormLinearSystem(ess_bdr_tdofs, *a, rhs, M, A, RHS);
+
+   // Define and apply a parallel PCG solver for M A = RHS with the AMS
+   // preconditioner from hypre.
+   HypreAMS ams(M, HCurlFESpaceOld);
+   ams.SetSingularProblem();
+
+   HyprePCG pcg(M);
+   pcg.SetTol(1e-12);
+   pcg.SetMaxIter(50);
+   pcg.SetPrintLevel(2);
+   pcg.SetPreconditioner(ams);
+   pcg.Mult(RHS, A);
+
+   // Extract the parallel grid function corresponding to the finite
+   // element approximation A. This is the local solution on each
+   // processor.
+   curlCurl->RecoverFEMSolution(A, rhs, *a);
+
+   //Compute the reconstructed b field for comparison
+   curl_old->Mult(*a, *recon_b);
+
+   //
+   Vector diff(*b);
+   diff -= *recon_b;    //diff = b - recon_b
+   std::cout << "L2 Error in reconstructed B field on old mesh:  " << diff.Norml2()
+             << std::endl;
+
+}
+
+
+//Given b in Hdiv compute the curl of b_ in Hcurl
+//and then clean any divergence out of it
+void BFieldAdvector::ComputeCleanCurlB(ParGridFunction* b)
+{
+   HypreParMatrix M1;
+   ParGridFunction rhs(HCurlFESpaceOld);
+   Vector RHS(HCurlFESpaceOld->GetTrueVSize());
+   Vector X(RHS.Size());
+   Vector P(RHS.Size());
+
+   Array<int> ess_tdof_list;
+   Array<int> ess_bdr;
+   if (pmeshOld->bdr_attributes.Size())
+   {
+      ess_bdr.SetSize(pmeshOld->bdr_attributes.Max());
+      ess_bdr = 1;
+      HCurlFESpaceOld->GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   b->GetTrueDofs(P);
+   WC->Mult(P,RHS);
+   rhs.SetFromTrueDofs(RHS);
+   *curl_b = 0.0;
+   m1->FormLinearSystem(ess_tdof_list, *curl_b, rhs, M1, X, RHS);
+
+   HypreDiagScale Jacobi(M1);
+   HyprePCG pcg(M1);
+   pcg.SetTol(1e-12);
+   pcg.SetMaxIter(1000);
+   pcg.SetPrintLevel(2);
+   pcg.SetPreconditioner(Jacobi);
+   X = 0.0;
+   pcg.Mult(RHS, X);
+
+   m1->RecoverFEMSolution(X, rhs, *curl_b);
+
+   divFreeProj->Mult(*curl_b, *clean_curl_b);
+}
+
+
+void BFieldAdvector::FindPtsInterpolateToTargetMesh(const ParGridFunction
+                                                    *old_gf, ParGridFunction *new_gf, int fieldtype)
+{
+   MFEM_ASSERT(fieldtype >= 0 &&
+               fieldtype <= 3, "Method expects a field type of 0, 1, 2, or 3");
+
+   int dim = pmeshOld->Dimension();
+   int vdim = old_gf->VectorDim();
+   int num_target_elem = pmeshOld->GetNE();
+   ParFiniteElementSpace *target_fes = new_gf->ParFESpace();
+
+   // Loop through the elements in case we have a mixed mesh
+   int num_target_pts = 0;
+   for (int e = 0; e < num_target_elem; ++e)
+   {
+      num_target_pts += target_fes->GetFE(e)->GetNodes().GetNPoints();
+   }
+
+   //Extract the target points from the nodes of the elements of the
+   //new mesh and then line them up in a vector V(x1,x2,...,y1,y2...,z1,z2...)
+   Vector vxyz(num_target_pts*dim);
+   int vxyz_pos = 0;
+   for (int e = 0; e < num_target_elem; e++)
+   {
+      const FiniteElement *fe = target_fes->GetFE(e);
+      const IntegrationRule ir = fe->GetNodes();
+      int elem_num_nodes = fe->GetNodes().GetNPoints();
+      ElementTransformation *trans = target_fes->GetElementTransformation(e);
+
+      DenseMatrix pos;
+      trans->Transform(ir, pos);
+      Vector rowx(vxyz.GetData() + vxyz_pos, elem_num_nodes),
+             rowy(vxyz.GetData() + num_target_pts + vxyz_pos, elem_num_nodes),
+             rowz;
+      if (dim == 3)
+      {
+         rowz.SetDataAndSize(vxyz.GetData() + 2*num_target_pts + vxyz_pos,
+                             elem_num_nodes);
+      }
+      pos.GetRow(0, rowx);
+      pos.GetRow(1, rowy);
+      if (dim == 3) { pos.GetRow(2, rowz); }
+      vxyz_pos += elem_num_nodes;
+   }
+
+   // Interpolate the values at the new_gf nodes
+   Vector interp_vals(num_target_pts*vdim);
+   FindPointsGSLIB finder(myComm);
+   finder.Setup(*pmeshOld);
+   finder.Interpolate(vxyz, *old_gf, interp_vals);
+
+
+   //I'll need to integrat this in and think about the differences between
+   //The H1/L2 versions and the ND/RT versions.
+   //It may be a good idea to have isH1, isL2, isND, and isRT methods in the FEC
+   // Project the interpolated values to the target FiniteElementSpace.
+   if (fieldtype == 0 || fieldtype == 3) // H1 or L2
+   {
+      if ((fieldtype == 0) || fieldtype == 3)
+      {
+         (*new_gf) = interp_vals;
+      }
+      else // H1 - but mesh order != GridFunction order
+      {
+         Array<int> vdofs;
+         //Vector vals;
+         int ivals_pos = 0;
+         for (int e = 0; e < num_target_elem; e++)
+         {
+            const FiniteElement *fe = target_fes->GetFE(e);
+            int elem_num_nodes = fe->GetNodes().GetNPoints();
+            Vector elem_dof_vals(elem_num_nodes*vdim);
+
+            target_fes->GetElementVDofs(e, vdofs);
+            //vals.SetSize(vdofs.Size());
+            for (int j = 0; j < elem_num_nodes; j++)
+            {
+               for (int d = 0; d < vdim; d++)
+               {
+                  // Arrange values byNodes
+                  elem_dof_vals(j+d*elem_num_nodes) = interp_vals(d*num_target_pts + ivals_pos +
+                                                                  j);
+               }
+            }
+            new_gf->SetSubVector(vdofs, elem_dof_vals);
+            ivals_pos += elem_num_nodes;
+         }
+      }
+   }
+   else // H(div) or H(curl)
+   {
+      std::cout << "Writing data into H(curl)/H(Div) dofs" <<std::endl;
+      Array<int> vdofs;
+      Vector vals;
+      int ivals_pos = 0;
+      for (int e = 0; e < num_target_elem; e++)
+      {
+         const FiniteElement *fe = target_fes->GetFE(e);
+         int elem_num_nodes = fe->GetNodes().GetNPoints();
+         Vector elem_dof_vals(elem_num_nodes*vdim);
+         target_fes->GetElementVDofs(e, vdofs);
+         vals.SetSize(vdofs.Size());
+         for (int j = 0; j < elem_num_nodes; j++)
+         {
+            for (int d = 0; d < vdim; d++)
+            {
+               // Arrange values byVDim
+               elem_dof_vals(j*vdim+d) = interp_vals(d*num_target_pts + ivals_pos + j);
+            }
+         }
+         fe->ProjectFromNodes(elem_dof_vals,
+                              *target_fes->GetElementTransformation(e),
+                              vals);
+         new_gf->SetSubVector(vdofs, vals);
+         ivals_pos += elem_num_nodes;
+      }
+   }
+}
+
+
+void BFieldAdvector::RemhosRemap(ParGridFunction* b_old, ParGridFunction* b_new)
+{
+#if 0
+   const int NE  = pmeshOld->GetNE();
+   const int mesh_order = pmeshOld->GetNodalFESpace()->GetMaxElementOrder();
+   const int dim = pmeshOld->Dimension();
+   const bool verify_bounds = false;
+   const bool forced_bounds = false;
+   const int problem_num = 11;
+   const int vis_steps = 10;
+   const int bounds_type = 0;
+   double t_final = 1.0;
+
+
+   FiniteElementCollection *mesh_fec = new H1_FECollection(mesh_order, dim,
+                                                           BasisType::GaussLobatto);
+
+   // Define the ODE solver used for time integration.
+   ODESolver *ode_solver = new RK3SSPSolver;
+
+
+   // Current mesh positions.
+   ParFiniteElementSpace mesh_pfes(pmeshOld, mesh_fec, dim);
+   ParGridFunction x(&mesh_pfes);
+   pmeshOld->SetNodalGridFunction(&x);
+
+   ParFiniteElementSpace mesh_pfes_new(pmeshNew, mesh_fec, dim);
+   ParGridFunction x_new(&mesh_pfes);
+   pmeshOld->SetNodalGridFunction(&x_new);
+
+   // Store initial mesh positions.
+   Vector x0(x.Size());
+   x0 = x;
+
+   // Initial time step estimate.
+   // Since we are talking about a local change lets do a fixed number of steps
+   // TODO:  Ask Vladimir about this assumption
+   double dt = 0.25;
+
+   // Mesh velocity.
+   // Obtain the mesh velocity by moving the mesh to the final
+   // mesh positions, and taking the displacement vector.
+   // The mesh motion resembles a time-dependent deformation, e.g., similar to
+   // a deformation that is obtained by a Lagrangian simulation.
+   GridFunction v_gf(x.FESpace());
+   VectorGridFunctionCoefficient v_mesh_coeff(&v_gf);
+   ParGridFunction v(&mesh_pfes);
+   v = x_new - x;
+   double t = 0.0;
+   while (t < t_final)
+   {
+      t += dt;
+      // Move the mesh nodes.
+      x.Add(std::min(dt, t_final-t), v);
+      //No need to update v here since the pseudo velocity doesn't change
+      //During this remap step
+      //TODO:  For higher order we will need a velocity function that can
+      //change over the course of the pseudo-timestep
+   }
+   add(x, -1.0, x0, v_gf); // Pseudotime velocity.
+   x = x0; // Return the mesh to the initial configuration.
+
+   const int btype = BasisType::Positive;
+   DG_FECollection fec(order, dim, btype);
+   ParFiniteElementSpace pfes(pmeshOld, &fec);
+
+   ParGridFunction inflow_gf(&pfes);
+   inflow_gf = 0.0;
+
+   // Set up the bilinear and linear forms corresponding to the DG
+   // discretization.
+   ParBilinearForm m(&pfes);
+   m.AddDomainIntegrator(new MassIntegrator);
+
+   ParBilinearForm M_HO(&pfes);
+   M_HO.AddDomainIntegrator(new MassIntegrator);
+
+   ParBilinearForm k(&pfes);
+   ParBilinearForm K_HO(&pfes);
+   k.AddDomainIntegrator(new ConvectionIntegrator(v_mesh_coeff));
+   K_HO.AddDomainIntegrator(new ConvectionIntegrator(v_mesh_coeff));
+
+   auto dgt_i = new DGTraceIntegrator(v_mesh_coeff, -1.0, -0.5);
+   auto dgt_b = new DGTraceIntegrator(v_mesh_coeff, -1.0, -0.5);
+   K_HO.AddInteriorFaceIntegrator(new TransposeIntegrator(dgt_i));
+   K_HO.AddBdrFaceIntegrator(new TransposeIntegrator(dgt_b));
+   K_HO.KeepNbrBlock(true);
+
+   K_HO.SetAssemblyLevel(AssemblyLevel::FULL);
+   M_HO.Assemble();
+   K_HO.Assemble(0);
+   M_HO.Finalize();
+   K_HO.Finalize(0);
+
+   // Compute the lumped mass matrix.
+   Vector lumpedM;
+   ParBilinearForm ml(&pfes);
+   ml.AddDomainIntegrator(new LumpedIntegrator(new MassIntegrator));
+   ml.Assemble();
+   ml.Finalize();
+   ml.SpMat().GetDiag(lumpedM);
+
+   m.Assemble();
+   m.Finalize();
+   int skip_zeros = 0;
+   k.Assemble(skip_zeros);
+   k.Finalize(skip_zeros);
+
+   // Store topological dof data.
+   DofInfo dofs(pfes, bounds_type);
+
+   // Precompute data required for high and low order schemes. This could be put
+   // into a separate routine. I am using a struct now because the various
+   // schemes require quite different information.
+   LowOrderMethod lom;
+   lom.subcell_scheme = false;
+
+   lom.pk = NULL;
+   lom.smap = SparseMatrix_Build_smap(k.SpMat());
+   lom.D = k.SpMat();
+   lom.coef = &v_mesh_coeff;
+
+   // Face integration rule.
+   const FaceElementTransformations *ft =
+      pmeshOld->GetFaceElementTransformations(0);
+   const int el_order = pfes.GetFE(0)->GetOrder();
+   int ft_order = ft->Elem1->OrderW() + 2 * el_order;
+   if (pfes.GetFE(0)->Space() == FunctionSpace::Pk) { ft_order++; }
+   lom.irF = &IntRules.Get(ft->FaceGeom, ft_order);
+
+   DG_FECollection fec0(0, dim, btype);
+   DG_FECollection fec1(1, dim, btype);
+
+   ParMesh *subcell_mesh = NULL;
+   lom.SubFes0 = NULL;
+   lom.SubFes1 = NULL;
+   FiniteElementCollection *fec_sub = NULL;
+   ParFiniteElementSpace *pfes_sub = NULL;;
+   ParGridFunction *xsub = NULL;
+   ParGridFunction v_sub_gf;
+   VectorGridFunctionCoefficient v_sub_coef;
+   Vector x0_sub;
+
+   if (order > 1)
+   {
+      // The mesh corresponding to Bezier subcells of order p is constructed.
+      // NOTE: The mesh is assumed to consist of quads or hexes.
+      MFEM_VERIFY(order > 1, "This code should not be entered for order = 1.");
+
+      // Get a uniformly refined mesh.
+      const int btype = BasisType::ClosedUniform;
+      subcell_mesh = new ParMesh(ParMesh::MakeRefined(*pmeshOld, order, btype));
+
+      // Check if the mesh is periodic.
+      const L2_FECollection *L2_coll = dynamic_cast<const L2_FECollection *>
+                                       (pmeshOld->GetNodes()->FESpace()->FEColl());
+      // Standard non-periodic mesh.
+      // Note that the fine mesh is always linear.
+      fec_sub = new H1_FECollection(1, dim, BasisType::ClosedUniform);
+      pfes_sub = new ParFiniteElementSpace(subcell_mesh, fec_sub, dim);
+      xsub = new ParGridFunction(pfes_sub);
+      subcell_mesh->SetCurvature(1);
+      subcell_mesh->SetNodalGridFunction(xsub);
+
+      lom.SubFes0 = new FiniteElementSpace(subcell_mesh, &fec0);
+      lom.SubFes1 = new FiniteElementSpace(subcell_mesh, &fec1);
+
+      // Submesh velocity.
+      v_sub_gf.SetSpace(pfes_sub);
+      v_sub_gf.ProjectCoefficient(v_mesh_coeff);
+
+      // Zero it out on boundaries (not moving boundaries).
+      Array<int> ess_bdr, ess_vdofs;
+      if (subcell_mesh->bdr_attributes.Size() > 0)
+      {
+         ess_bdr.SetSize(subcell_mesh->bdr_attributes.Max());
+      }
+      ess_bdr = 1;
+      xsub->ParFESpace()->GetEssentialVDofs(ess_bdr, ess_vdofs);
+      for (int i = 0; i < ess_vdofs.Size(); i++)
+      {
+         if (ess_vdofs[i] == -1) { v_sub_gf(i) = 0.0; }
+      }
+      v_sub_coef.SetGridFunction(&v_sub_gf);
+
+      // Store initial submesh positions.
+      x0_sub = *xsub;
+
+      lom.subcellCoeff = &v_sub_coef;
+      lom.VolumeTerms = new MixedConvectionIntegrator(v_sub_coef);
+   }
+   else { subcell_mesh = pmeshOld; }
+
+   Assembly asmbl(dofs, lom, inflow_gf, pfes, subcell_mesh, 1);
+
+   // Setup the initial conditions.
+   const int vsize = pfes.GetVSize();
+   Array<int> offset(2);      //2 because we are assuming product_sync = 0
+   for (int i = 0; i < offset.Size(); i++) { offset[i] = i*vsize; }
+   BlockVector S(offset, Device::GetMemoryType());
+   // Primary scalar field is u.
+   ParGridFunction u(&pfes);
+   u.MakeRef(&pfes, S, offset[0]);
+   u = *b_old;    //Set u to the pre remap b state
+   u.SyncAliasMemory(S);
+
+   //No Product sync
+
+   //No Smoothness indicator
+   SmoothnessIndicator *smth_indicator = NULL;
+
+   // Setup of the high-order solver
+   HOSolver *ho_solver = new LocalInverseHOSolver(pfes, M_HO, K_HO);
+
+   // Setup the low order solver
+   const bool time_dep = true;
+   Array<int> lo_smap = SparseMatrix_Build_smap(k.SpMat());
+   LOSolver *lo_solver = new DiscreteUpwind(pfes, k.SpMat(), lo_smap,
+                                            lumpedM, asmbl, time_dep);
+
+   // Setup of the FCT solver.
+   Array<int> K_HO_smap;
+   FCTSolver *fct_solver = NULL;
+   K_HO.SpMat().HostReadI();
+   K_HO.SpMat().HostReadJ();
+   K_HO.SpMat().HostReadData();
+   K_HO_smap = SparseMatrix_Build_smap(K_HO.SpMat());
+   const int fct_iterations = 1;
+   fct_solver = new FluxBasedFCT(pfes, smth_indicator, dt, K_HO.SpMat(),
+                                 K_HO_smap, M_HO.SpMat(), fct_iterations);
+
+   AdvectionOperator adv(S.Size(), m, ml, lumpedM, k, M_HO, K_HO,
+                         x, xsub, v_gf, v_sub_gf, asmbl, lom, dofs,
+                         ho_solver, lo_solver, fct_solver);
+
+   t = 0.0;
+   adv.SetTime(t);
+   ode_solver->Init(adv);
+
+   double umin, umax;
+   GetMinMax(u, umin, umax);
+
+   adv.SetRemapStartPos(x0, x0_sub);
+
+
+   ParGridFunction res = u;
+   double residual;
+   double s_min_glob = numeric_limits<double>::infinity(),
+          s_max_glob = -numeric_limits<double>::infinity();
+
+   // Time-integration (loop over the time iterations, ti, with a time-step dt).
+   bool done = false;
+   BlockVector Sold(S);
+   int ti_total = 0, ti = 0;
+   while (done == false)
+   {
+      double dt_real = min(dt, t_final - t);
+
+      // This also resets the time step estimate when automatic dt is on.
+      adv.SetDt(dt_real);
+      if (lo_solver)  { lo_solver->UpdateTimeStep(dt_real); }
+      if (fct_solver) { fct_solver->UpdateTimeStep(dt_real); }
+
+      Sold = S;
+      ode_solver->Step(S, t, dt_real);
+      ti++;
+      ti_total++;
+
+      //Time step control is fixed
+
+      // S has been modified, update the alias
+      u.SyncMemory(S);
+
+      // Monotonicity check for debug purposes mainly.
+      if (verify_bounds && forced_bounds && smth_indicator == NULL)
+      {
+         double umin_new, umax_new;
+         GetMinMax(u, umin_new, umax_new);
+         if (problem_num % 10 != 6 && problem_num % 10 != 7)
+         {
+            if (pmeshOld->GetMyRank() == 0)
+            {
+               MFEM_VERIFY(umin_new > umin - 1e-12,
+                           "Undershoot of " << umin - umin_new);
+               MFEM_VERIFY(umax_new < umax + 1e-12,
+                           "Overshoot of " << umax_new - umax);
+            }
+            umin = umin_new;
+            umax = umax_new;
+         }
+         else
+         {
+            if (pmeshOld->GetMyRank() == 0)
+            {
+               MFEM_VERIFY(umin_new > 0.0 - 1e-12,
+                           "Undershoot of " << 0.0 - umin_new);
+               MFEM_VERIFY(umax_new < 1.0 + 1e-12,
+                           "Overshoot of " << umax_new - 1.0);
+            }
+         }
+      }
+
+      x0.HostReadWrite(); v_sub_gf.HostReadWrite();
+      x.HostReadWrite();
+      add(x0, t, v_gf, x);
+      x0_sub.HostReadWrite(); v_sub_gf.HostReadWrite();
+      MFEM_VERIFY(xsub != NULL,
+                  "xsub == NULL/This code should not be entered for order = 1.");
+      xsub->HostReadWrite();
+      add(x0_sub, t, v_sub_gf, *xsub);
+
+      done = (t >= t_final - 1.e-8*dt);
+      if (done || ti % vis_steps == 0)
+      {
+         if (pmeshOld->GetMyRank() == 0)
+         {
+            std::cout << "time step: " << ti << ", time: " << t
+                      << ", dt: " << dt << ", residual: " << residual << std::endl;
+         }
+      }
+   }
+
+   std::cout << "Total time steps: " << ti_total
+             << " (" << ti_total-ti << " repeated)." << std::endl;
+
+   *b_new = u;
+#endif
+}
+
+
+
+
+AdvectionOperator::AdvectionOperator(int size, BilinearForm &Mbf_,
+                                     BilinearForm &_ml, Vector &_lumpedM,
+                                     ParBilinearForm &Kbf_,
+                                     ParBilinearForm &M_HO_, ParBilinearForm &K_HO_,
+                                     GridFunction &pos, GridFunction *sub_pos,
+                                     GridFunction &vel, GridFunction &sub_vel,
+                                     Assembly &_asmbl,
+                                     LowOrderMethod &_lom, DofInfo &_dofs,
+                                     HOSolver *hos, LOSolver *los, FCTSolver *fct) :
+   TimeDependentOperator(size), Mbf(Mbf_), ml(_ml), Kbf(Kbf_),
+   M_HO(M_HO_), K_HO(K_HO_),
+   lumpedM(_lumpedM),
+   start_mesh_pos(pos.Size()), start_submesh_pos(sub_vel.Size()),
+   mesh_pos(pos), submesh_pos(sub_pos),
+   mesh_vel(vel), submesh_vel(sub_vel),
+   x_gf(Kbf.ParFESpace()),
+   asmbl(_asmbl), lom(_lom), dofs(_dofs),
+   ho_solver(hos), lo_solver(los), fct_solver(fct)
+{
+   MFEM_VERIFY(fct_solver && ho_solver &&
+               lo_solver, "Bfield Remhos advector requires FCT, ho and lo solvers.");
+}
+
+void AdvectionOperator::Mult(const Vector &X, Vector &Y) const
+{
+   MFEM_VERIFY(ho_solver && lo_solver, "FCT requires HO and LO solvers.");
+
+   // Move the mesh positions.
+   const double t = GetTime();
+   add(start_mesh_pos, t, mesh_vel, mesh_pos);
+   if (submesh_pos)
+   {
+      add(start_submesh_pos, t, submesh_vel, *submesh_pos);
+   }
+   // Reset precomputed geometric data.
+   Mbf.FESpace()->GetMesh()->DeleteGeometricFactors();
+
+   // Reassemble on the new mesh. Element contributions.
+   // Currently needed to have the sparse matrices used by the LO methods.
+   Mbf.BilinearForm::operator=(0.0);
+   Mbf.Assemble();
+   Kbf.BilinearForm::operator=(0.0);
+   Kbf.Assemble(0);
+   ml.BilinearForm::operator=(0.0);
+   ml.Assemble();
+   lumpedM.HostReadWrite();
+   ml.SpMat().GetDiag(lumpedM);
+
+   M_HO.BilinearForm::operator=(0.0);
+   M_HO.Assemble();
+   K_HO.BilinearForm::operator=(0.0);
+   K_HO.Assemble(0);
+
+   if (lom.pk)
+   {
+      lom.pk->BilinearForm::operator=(0.0);
+      lom.pk->Assemble();
+   }
+
+   // Face contributions.
+   asmbl.bdrInt = 0.;
+   Mesh *mesh = M_HO.FESpace()->GetMesh();
+   const int dim = mesh->Dimension(), ne = mesh->GetNE();
+   Array<int> bdrs, orientation;
+   FaceElementTransformations *Trans;
+
+   for (int k = 0; k < ne; k++)
+   {
+      if (dim == 1)      { mesh->GetElementVertices(k, bdrs); }
+      else if (dim == 2) { mesh->GetElementEdges(k, bdrs, orientation); }
+      else if (dim == 3) { mesh->GetElementFaces(k, bdrs, orientation); }
+
+      for (int i = 0; i < dofs.numBdrs; i++)
+      {
+         Trans = mesh->GetFaceElementTransformations(bdrs[i]);
+         asmbl.ComputeFluxTerms(k, i, Trans, lom);
+      }
+   }
+
+   const int size = Kbf.ParFESpace()->GetVSize();
+   const int NE   = Kbf.ParFESpace()->GetNE();
+
+   // Needed because X and Y are allocated on the host by the ODESolver.
+   X.Read(); Y.Read();
+
+   Vector u, d_u;
+   Vector* xptr = const_cast<Vector*>(&X);
+   u.MakeRef(*xptr, 0, size);
+   d_u.MakeRef(Y, 0, size);
+   Vector du_HO(u.Size()), du_LO(u.Size());
+
+   x_gf = u;
+   x_gf.ExchangeFaceNbrData();
+
+   if (fct_solver)
+   {
+      MFEM_VERIFY(ho_solver && lo_solver, "FCT requires HO and LO solvers.");
+
+      lo_solver->CalcLOSolution(u, du_LO);
+      ho_solver->CalcHOSolution(u, du_HO);
+
+      dofs.ComputeElementsMinMax(u, dofs.xe_min, dofs.xe_max, NULL, NULL);
+      dofs.ComputeBounds(dofs.xe_min, dofs.xe_max, dofs.xi_min, dofs.xi_max);
+      fct_solver->CalcFCTSolution(x_gf, lumpedM, du_HO, du_LO,
+                                  dofs.xi_min, dofs.xi_max, d_u);
+   }
+
+   d_u.SyncAliasMemory(Y);
+
+   // Remap the product field, if there is a product field.
+   if (X.Size() > size)
+   {
+      Vector us, d_us;
+      us.MakeRef(*xptr, size, size);
+      d_us.MakeRef(Y, size, size);
+
+      x_gf = us;
+      x_gf.ExchangeFaceNbrData();
+
+      if (fct_solver)
+      {
+         MFEM_VERIFY(ho_solver && lo_solver, "FCT requires HO and LO solvers.");
+
+         Vector d_us_HO(us.Size()), d_us_LO;
+         if (fct_solver->NeedsLOProductInput())
+         {
+            d_us_LO.SetSize(us.Size());
+            lo_solver->CalcLOSolution(us, d_us_LO);
+         }
+         ho_solver->CalcHOSolution(us, d_us_HO);
+
+         // Compute the ratio s = us_old / u_old, and old active dofs.
+         Vector s(size);
+         Array<bool> s_bool_el, s_bool_dofs;
+         ComputeRatio(NE, us, u, s, s_bool_el, s_bool_dofs);
+
+         // Bounds for s, based on the old values (and old active dofs).
+         // This doesn't consider s values from the old inactive dofs, because
+         // there were no bounds restriction on them at the previous time step.
+         dofs.ComputeElementsMinMax(s, dofs.xe_min, dofs.xe_max,
+                                    &s_bool_el, &s_bool_dofs);
+         dofs.ComputeBounds(dofs.xe_min, dofs.xe_max,
+                            dofs.xi_min, dofs.xi_max, &s_bool_el);
+
+         // Evolve u and get the new active dofs.
+         Vector u_new(size);
+         add(1.0, u, dt, d_u, u_new);
+         Array<bool> s_bool_el_new, s_bool_dofs_new;
+         ComputeBoolIndicators(NE, u_new, s_bool_el_new, s_bool_dofs_new);
+
+         fct_solver->CalcFCTProduct(x_gf, lumpedM, d_us_HO, d_us_LO,
+                                    dofs.xi_min, dofs.xi_max,
+                                    u_new,
+                                    s_bool_el_new, s_bool_dofs_new, d_us);
+      }
+
+      d_us.SyncAliasMemory(Y);
+   }
+}
+
+} // namespace electromagnetics
+} // namespace mfem

miniapps/electromagnetics/bfieldadvect_solver.hpp

@@ -0,0 +1,130 @@
+// Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "../common/pfem_extras.hpp"
+#include "../common/mesh_extras.hpp"
+#include "bfield_remhos.hpp"
+
+namespace mfem
+{
+
+using common::H1_ParFESpace;
+using common::ND_ParFESpace;
+using common::RT_ParFESpace;
+using common::L2_ParFESpace;
+using common::ParDiscreteGradOperator;
+using common::ParDiscreteCurlOperator;
+using common::DivergenceFreeProjector;
+
+namespace electromagnetics
+{
+
+class BFieldAdvector
+{
+public:
+   BFieldAdvector(ParMesh *pmesh_old, ParMesh *pmesh_new, int order);
+   void SetMesh(ParMesh *pmesh_old, ParMesh *pmesh_new);
+   void SetMeshNodes(ParGridFunction *old_nodes, ParGridFunction *new_nodes);
+   void Advect(ParGridFunction* b_old, ParGridFunction* b_new);
+
+   ParGridFunction* GetVectorPotential() {return a;}
+   ParGridFunction* GetCurlB() {return curl_b;}
+   ParGridFunction* GetCleanCurlB() {return clean_curl_b;}
+   ParGridFunction* GetReconstructedB() {return recon_b;}
+   ParGridFunction* GetA() {return a;}
+   ParGridFunction* GetANew() {return a_new;}
+
+private:
+   void CleanInternals();
+   void ComputeCleanCurlB(ParGridFunction* b);
+   void ComputeA(ParGridFunction* b);
+
+   /// Given a grid function on the old and new meshes interpolate a field from the old to the new
+   /** The fieldtype variable goes from 0 forms to 3 forms in order 0-H1, 1-H(curl), 2-H(div), 3-L2
+   **/
+   void FindPtsInterpolateToTargetMesh(const ParGridFunction *old_gf,
+                                       ParGridFunction *new_gf, int fieldtype);
+   void RemhosRemap(ParGridFunction* b_old, ParGridFunction* b_new);
+
+   int order;
+   MPI_Comm myComm;
+   ParMesh *pmeshOld, *pmeshNew;
+   H1_ParFESpace *H1FESpaceOld, *H1FESpaceNew;
+   ND_ParFESpace *HCurlFESpaceOld, *HCurlFESpaceNew;
+   RT_ParFESpace *HDivFESpaceOld, *HDivFESpaceNew;
+   L2_ParFESpace *L2FESpaceOld, *L2FESpaceNew;
+   ParDiscreteGradOperator *grad;
+   ParDiscreteCurlOperator *curl_old;
+   ParDiscreteCurlOperator *curl_new;
+   ParMixedBilinearForm *weakCurl;
+   HypreParMatrix *WC;
+   ParBilinearForm *m1;
+   ParBilinearForm *curlCurl;
+   DivergenceFreeProjector *divFreeProj;
+   ParGridFunction *a;
+   ParGridFunction *a_new;
+   ParGridFunction *curl_b;
+   ParGridFunction *clean_curl_b;
+   ParGridFunction *recon_b;
+};
+
+
+class AdvectionOperator : public TimeDependentOperator
+{
+private:
+   BilinearForm &Mbf, &ml;
+   ParBilinearForm &Kbf;
+   ParBilinearForm &M_HO, &K_HO;
+   Vector &lumpedM;
+
+   Vector start_mesh_pos, start_submesh_pos;
+   GridFunction &mesh_pos, *submesh_pos, &mesh_vel, &submesh_vel;
+
+   mutable ParGridFunction x_gf;
+
+   double dt;
+   mutable double dt_est;
+   Assembly &asmbl;
+
+   LowOrderMethod &lom;
+   DofInfo &dofs;
+
+   HOSolver *ho_solver;
+   LOSolver *lo_solver;
+   FCTSolver *fct_solver;
+
+public:
+   AdvectionOperator(int size, BilinearForm &Mbf_, BilinearForm &_ml,
+                     Vector &_lumpedM,
+                     ParBilinearForm &Kbf_,
+                     ParBilinearForm &M_HO_, ParBilinearForm &K_HO_,
+                     GridFunction &pos, GridFunction *sub_pos,
+                     GridFunction &vel, GridFunction &sub_vel,
+                     Assembly &_asmbl, LowOrderMethod &_lom, DofInfo &_dofs,
+                     HOSolver *hos, LOSolver *los, FCTSolver *fct);
+
+   virtual void Mult(const Vector &x, Vector &y) const;
+
+   void SetDt(double dt_) { dt = dt_; dt_est = dt; }
+   double GetTimeStepEstimate() { return dt_est; }
+
+   void SetRemapStartPos(const Vector &m_pos, const Vector &sm_pos)
+   {
+      start_mesh_pos    = m_pos;
+      start_submesh_pos = sm_pos;
+   }
+
+   virtual ~AdvectionOperator() { }
+};
+
+} // namespace electromagnetics
+
+} // namespace mfem

miniapps/electromagnetics/makefile

@@ -26,7 +26,7 @@ MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
 SEQ_MINIAPPS =
-PAR_MINIAPPS = volta tesla maxwell joule
+PAR_MINIAPPS = volta tesla maxwell joule bfieldadvect
 ifeq ($(MFEM_USE_MPI),NO)
    MINIAPPS = $(SEQ_MINIAPPS)
 else
@@ -54,8 +54,9 @@ all: $(MINIAPPS)
 # Rules for building the miniapps
 %: $(SRC)%.cpp %_solver.o $(MFEM_LIB_FILE) $(CONFIG_MK) | lib-common
 	$(MFEM_CXX) $(MFEM_FLAGS) -c $(<)
+	$(MFEM_CXX) $(MFEM_FLAGS) -c bfield_remhos.cpp
 	$(MFEM_CXX) $(MFEM_LINK_FLAGS) -o $@ $@.o $@_solver.o $(COMMON_LIB) \
-	  $(MFEM_LIBS)
+	  $(MFEM_LIBS) bfield_remhos.o
 
 # Rules for compiling miniapp dependencies
 $(addsuffix _solver.o,$(MINIAPPS)): \

miniapps/meshing/ext-mesh-mapping.cpp

@@ -0,0 +1,268 @@
+// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+//
+//             ---------------------------------------------------
+//                       External Mesh Mapping Miniapp
+//             ---------------------------------------------------
+//
+// This miniapp starts with a serial non-conforming mesh in a dummy format and demonstrates
+// how to build up a corresponding non-conforming MFEM Mesh and then decompose it
+// into a parallel ParMesh.  As part of this process we will demonstrate how to obtain
+// and compose the vertex ID mappings from the various steps that can shuffle the
+// vertices.  In the end this will let us map between the vertex ID numbers in the
+// external dummy mesh and the parallel non-conforming mesh constructed in MFEM.  This
+// is the sort of thing you will have to do if you intend to add MFEM meshes to an existing
+// simulation code and need them to exist and share data with other kinds of meshes in that
+// code.  If you are starting a new MFEM code, it makes much more sense to do everything with
+// MFEM meshes.
+//
+// Compile with: make ext-mesh-mapping
+//
+// Sample runs:   mpirun -np 4 ext_mesh_mapping
+
+#include "mfem.hpp"
+#include "ext-mesh-mapping.hpp"
+
+using namespace mfem;
+
+Mesh *build_mfem_mesh(DummyMesh *dmesh);
+void create_pmesh_to_mesh_vmaps(ParMesh *pmesh, Mesh *mesh, Array<int> &vmap);
+
+void print_dmesh_verts(DummyMesh *dmesh);
+void print_mesh_verts(Mesh *mesh, const Array<int> &vmap = Array<int>());
+void print_pmesh_verts(ParMesh *pmesh, const Array<int> &vmap = Array<int>());
+
+int main(int argc, char *argv[])
+{
+   Mpi::Init(argc, argv);
+   Hypre::Init();
+
+   //Initilize our dummy mesh and display the coordinates of the vertices
+   DummyMesh *dmesh = new DummyMesh();
+   print_dmesh_verts(dmesh);
+
+   //Now build the mfem Mesh object with all of the elements in it on all processors and capture the vertex id mapping
+   //that occurs when non-conforming meshes are finalized.
+   Mesh * mesh = build_mfem_mesh(dmesh);
+   Array<int> mesh_to_dmesh_vmap;
+   const Array<int> vmap = mesh->ncmesh->GetVertexIDMap();
+   mesh_to_dmesh_vmap = vmap;
+   print_mesh_verts(mesh,
+                    mesh_to_dmesh_vmap);     //Print the vertices reordered using the vmap
+
+   //Now enable parallel, given the following partition of the elements.
+   //Note that we only have local vertex ids in the pmesh object.
+   int partition[5] = {0,0,1,2,3};
+   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh, partition);
+
+   //Now compute the mappings between the local vertex ids on each processor in pmesh
+   //and the global ids in mesh.
+   Array<int> pmesh_to_mesh_vmap;
+   create_pmesh_to_mesh_vmaps(pmesh, mesh, pmesh_to_mesh_vmap);
+
+   //Now compose the maps to create a final map between the local vertex numbering in the
+   //pmesh object on this processor and the original dmesh vertex numbering
+   Array<int> final_vmap(pmesh->GetNV());
+   for (int local_vi = 0; local_vi < pmesh->GetNV(); ++local_vi)
+   {
+      final_vmap[local_vi] = mesh_to_dmesh_vmap[pmesh_to_mesh_vmap[local_vi]];
+   }
+   print_pmesh_verts(pmesh,
+                     final_vmap);           //Print the vertices with global IDs from the final_vmap
+}
+
+
+// Build up an MFEM Mesh from the data in the Dummy Mesh.  Since we are
+// building the vertex and element lists in the same order as we found them
+// in the dmesh, the element id and vertex id mappings will be the identity maps.
+Mesh *build_mfem_mesh(DummyMesh *dmesh)
+{
+   //Initilize the the dimension and memory for the mesh
+   Mesh *mesh = new Mesh(2,    // The dimension of the mesh
+                         dmesh->num_vertices,
+                         dmesh->num_elements,
+                         dmesh->num_belements,
+                         2    // The dimension of the space the mesh lives in (different for surface meshes)
+                        );
+
+
+   //Add the vertices into the mesh in the same order as they are in the dmesh
+   for (int vid = 0; vid < dmesh->num_vertices; ++vid)
+   {
+      mesh->AddVertex(dmesh->V[vid].x, dmesh->V[vid].y);
+   }
+
+   //Add the elements into the mesh in the same order as they are in the dmesh
+   //The dmesh is purely quads, but if it had other types we would use
+   //the other mesh.Add* methods here as well.
+   for (int eid = 0; eid < dmesh->num_elements; ++eid)
+   {
+      std::array<int,4> vid = dmesh->E[eid].vertex_ids;
+
+      //The order of vertices in the MFEM quad elements is different from the
+      //order of the vertices in out DummyMesh format so we must permute them
+      //here.  See the ref-*.mesh files in mfem/data to establish to ordering
+      //of the vertices in the MFEM elements.
+      //  Dummy    MFEM
+      //  2--3     3--2
+      //  |  |     |  |
+      //  0--1     0--1
+      mesh->AddQuad(vid[0], vid[1], vid[3], vid[2]);
+   }
+
+   //Add the boundary elements into the mesh in the same order they are in the dmesh
+   for (int bid = 0; bid < dmesh->num_belements; ++bid)
+   {
+      std::array<int,2> vid = dmesh->B[bid].vertex_ids;
+      mesh->AddBdrSegment(vid[0], vid[1]);
+   }
+
+   //Finally add vertex parents to mark element 0 for anisotropic refinement
+   for (int vpi = 0; vpi < dmesh->num_vparents; ++vpi)
+   {
+      mesh->AddVertexParents(std::get<0>(dmesh->VP[vpi]),std::get<1>(dmesh->VP[vpi]),
+                             std::get<2>(dmesh->VP[vpi]));
+   }
+
+   //This will make the mesh usable
+   mesh->FinalizeMesh();
+
+   return mesh;
+}
+
+// We can use the local elements with local vertex ids defined in pmesh and the global
+// elements defined with global vertex ids in mesh to define a mapping from the
+// local vertex ids on each processor to their global vertex id numbers.
+void create_pmesh_to_mesh_vmaps(ParMesh *pmesh, Mesh *mesh, Array<int> &vmap)
+{
+   vmap.SetSize(pmesh->GetNV());
+   for (int local_eid = 0; local_eid < pmesh->GetNE(); ++local_eid)
+   {
+      int global_eid = int(pmesh->GetGlobalElementNum(
+                              local_eid));   //This comes back as a long long
+      Array<int> local_elem_verts, global_elem_verts;
+      pmesh->GetElement(local_eid)->GetVertices(local_elem_verts);
+      mesh->GetElement(global_eid)->GetVertices(global_elem_verts);
+      for (int vi = 0; vi < local_elem_verts.Size(); ++vi)
+      {
+         vmap[local_elem_verts[vi]] = global_elem_verts[vi];
+      }
+   }
+}
+
+void print_dmesh_verts(DummyMesh *dmesh)
+{
+   if (Mpi::Root())
+   {
+      std::cout << "6-----7-----8" << std::endl;
+      std::cout << "|     |     |" << std::endl;
+      std::cout << "|  3  |  4  |" << std::endl;
+      std::cout << "|     |     |" << std::endl;
+      std::cout << "3-----4-----5" << std::endl;
+      std::cout << "|  1  |     |" << std::endl;
+      std::cout << "9----10  2  |" << std::endl;
+      std::cout << "|  0  |     |" << std::endl;
+      std::cout << "0-----1-----2" << std::endl;
+
+      std::cout << "DummyMesh vertices:  "  << std::endl;
+      for (int vid = 0; vid < dmesh->num_vertices; ++vid)
+      {
+         std::cout << vid << ":  " << dmesh->V[vid].x << ", " << dmesh->V[vid].y <<
+                   std::endl;
+      }
+      std::cout << std::endl;
+   }
+}
+
+void print_mesh_verts(Mesh *mesh, const Array<int> &vmap)
+{
+   if (Mpi::Root())
+   {
+      std::cout << std::endl << "MFEM Mesh Vertices:  "  << std::endl;
+      if (vmap.Size() > 0)
+      {
+         std::cout << "(Remapped vertex ids)" << std::endl;
+      }
+      for (int vid = 0; vid < mesh->GetNV(); ++vid)
+      {
+         int id = vmap.Size() > 0 ? vmap[vid] : vid;
+         double *vertex = mesh->GetVertex(id);
+         std::cout << vid << ":  " <<  vertex[0] << ", " << vertex[1] << std::endl;
+      }
+      std::cout << std::endl;
+   }
+}
+
+
+void print_pmesh_verts(ParMesh *pmesh, const Array<int> &vmap)
+{
+   int my_rank = Mpi::WorldRank();
+   int num_rank = Mpi::WorldSize();
+
+   Array<int> num_verts(num_rank);
+   int my_num_verts = pmesh->GetNV();
+   MPI_Allgather(&my_num_verts, 1, MPI_INTEGER,
+                 num_verts.GetData(), 1, MPI_INTEGER, MPI_COMM_WORLD);
+
+   int max_num_verts = *std::max_element(num_verts.begin(), num_verts.end());
+   Array<int> id_data(max_num_verts);
+   Array<double> x_data(max_num_verts);
+   Array<double> y_data(max_num_verts);
+
+   //Send the data to rank 0
+   for (int vid = 0; vid < pmesh->GetNV(); ++vid)
+   {
+      int id = vmap.Size() > 0 ? vmap[vid] : vid;
+      double *vertex = pmesh->GetVertex(vid);
+      id_data[vid] = id;
+      x_data[vid] =  vertex[0];
+      y_data[vid] =  vertex[1];
+   }
+   if (my_rank != 0)
+   {
+      MPI_Send(id_data.GetData(), pmesh->GetNV(), MPI_INTEGER, 0, 0, MPI_COMM_WORLD);
+      MPI_Send(x_data.GetData(), pmesh->GetNV(), MPI_DOUBLE, 0, 1, MPI_COMM_WORLD);
+      MPI_Send(y_data.GetData(), pmesh->GetNV(), MPI_DOUBLE, 0, 2, MPI_COMM_WORLD);
+   }
+
+
+   if (my_rank == 0)
+   {
+      std::cout << "MFEM ParMesh Vertices on each processor:  "  << std::endl;
+      if (vmap.Size() > 0)
+      {
+         std::cout << "(Remapped vertex ids)" << std::endl;
+      }
+
+      for (int p = 0; p < num_rank; ++p)
+      {
+         if (p != 0)
+         {
+            MPI_Status status;
+            MPI_Recv(id_data.GetData(), num_verts[p], MPI_INTEGER, p, 0, MPI_COMM_WORLD,
+                     &status);
+            MPI_Recv(x_data.GetData(), num_verts[p], MPI_DOUBLE, p, 1, MPI_COMM_WORLD,
+                     &status);
+            MPI_Recv(y_data.GetData(), num_verts[p], MPI_DOUBLE, p, 2, MPI_COMM_WORLD,
+                     &status);
+         }
+
+         for (int vid = 0; vid < num_verts[p]; ++vid)
+         {
+            std::cout << "rank (" << p << ") id (" << id_data[vid] << "):  "
+                      << x_data[vid]  << ", " << y_data[vid] << std::endl;
+         }
+      }
+      std::cout << std::endl;
+   }
+
+   MPI_Barrier(MPI_COMM_WORLD);
+}

miniapps/meshing/ext-mesh-mapping.hpp

@@ -0,0 +1,118 @@
+// Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+// MFEM External Mesh Mapping Miniapp
+
+#include <array>
+
+// A very simple example of an external mesh class with a baby mesh in it.
+//The element and vertex ids are as follows:
+// 6-----7-----8
+// |     |     |
+// |  3  |  4  |
+// |     |     |
+// 3-----4-----5
+// |  1  |     |
+// 9----10  2  |
+// |  0  |     |
+// 0-----1-----2
+class DummyMesh
+{
+   class Vertex
+   {
+   public:
+      Vertex(): x(0.0), y(0.0) {}
+
+      void Set(double x_, double y_)
+      {
+         x = x_;
+         y = y_;
+      }
+
+      double x;
+      double y;
+   };
+
+   class Element
+   {
+   public:
+      Element() {}
+
+      void Set(int id0, int id1, int id2, int id3)
+      {
+         vertex_ids[0] = id0; vertex_ids[1] = id1;
+         vertex_ids[2] = id2; vertex_ids[3] = id3;
+      }
+
+      const int num_vertices = 4;
+      std::array<int,4> vertex_ids;
+   };
+
+   class BElement
+   {
+   public:
+      BElement() {}
+
+      void Set(int id0, int id1)
+      {
+         vertex_ids[0] = id0; vertex_ids[1] = id1;
+      }
+
+      const int num_vertices = 2;
+      std::array<int,2> vertex_ids;
+   };
+
+public:
+   DummyMesh() : num_vertices(11), num_elements(5), num_belements(9),
+      num_vparents(2)
+   {
+      //Vertices
+      V[6].Set(0.0, 2.0);  V[7].Set(1.0, 2.0);  V[8].Set(2.0, 2.0);
+      V[3].Set(0.0, 1.0);  V[4].Set(1.0, 1.0);  V[5].Set(2.0, 1.0);
+      V[0].Set(0.0, 0.0);  V[1].Set(1.0, 0.0);  V[2].Set(2.0, 0.0);
+
+      V[9].Set(0.0, 0.5);  V[10].Set(1.0, 0.5);
+
+      // Elements in this dummy format have their vertex indices listed
+      // lexographic order rather than going around the element as is normal
+      // in MFEM.
+      E[0].Set(0,1,9,10);
+      E[1].Set(9,10,3,4);
+      E[2].Set(1,2,4,5);
+      E[3].Set(3,4,6,7);
+      E[4].Set(4,5,7,8);
+
+      //Boundary Elements
+      //Bottom
+      B[0].Set(0,1);       B[1].Set(1,2);
+      //Top
+      B[2].Set(6,7);       B[3].Set(7,8);
+      //Left
+      B[4].Set(0,9);       B[5].Set(9,3);       B[6].Set(3,6);
+      //Right
+      B[7].Set(2,5);       B[8].Set(5,8);
+
+      //Set the vertex parents
+      //In the future replace this with element parents and demonstrate
+      //computation of vertex parents
+      VP[0] = std::make_tuple(9,0,3);
+      VP[1] = std::make_tuple(10,1,4);
+   }
+
+   const int num_vertices;
+   const int num_elements;
+   const int num_belements;
+   const int num_vparents;
+   Vertex V[11];                    //Mesh vertices
+   Element E[5];                    //Mesh elements
+   BElement B[9];                   //Mesh boundary elements
+   std::tuple<int,int,int> VP[2];   //Vertex parents (vid, parent1id, parent2id)
+};

miniapps/meshing/makefile

@@ -27,7 +27,7 @@ MFEM_LIB_FILE = mfem_is_not_built
 
 SEQ_MINIAPPS = mobius-strip klein-bottle toroid trimmer twist mesh-explorer\
 	shaper extruder mesh-optimizer minimal-surface polar-nc
-PAR_MINIAPPS = pmesh-optimizer pminimal-surface
+PAR_MINIAPPS = pmesh-optimizer pminimal-surface ext_mesh_mapping
 ifeq ($(MFEM_USE_MPI),NO)
    MINIAPPS = $(SEQ_MINIAPPS)
 else

tests/unit/fem/test_calcdivshape.cpp

@@ -100,6 +100,7 @@ TEST_CASE("CalcDivShape RT",
           "[RT_QuadrilateralElement]"
           "[RT_TetrahedronElement]"
           "[RT_WedgeElement]"
+          "[RT_PyramidElement]"
           "[RT_HexahedronElement]")
 {
    const int maxOrder = 5;
@@ -144,6 +145,23 @@ TEST_CASE("CalcDivShape RT",
       TestCalcDivShape(&fe, &T, resolution);
    }
 
+   SECTION("RT_PyramidElement")
+   {
+      IsoparametricTransformation T;
+      GetReferenceTransformation(Element::PYRAMID, T);
+
+      if (order == 1)
+      {
+         RT0PyrFiniteElement fe;
+         TestCalcDivShape(&fe, &T, resolution);
+      }
+      else if (order == 2)
+      {
+         RT1PyrFiniteElement fe;
+         TestCalcDivShape(&fe, &T, resolution);
+      }
+   }
+
    SECTION("RT_HexahedronElement")
    {
       IsoparametricTransformation T;
@@ -153,3 +171,184 @@ TEST_CASE("CalcDivShape RT",
       TestCalcDivShape(&fe, &T, resolution);
    }
 }
+
+/**
+ * Tests fe->CalcDivShape() over a set of IntegrationPoints
+ * chosen based on the order. Compares the computed derivatives against
+ * approximate derivatives computed using the secant method.
+ */
+void TestFDCalcDivShape(FiniteElement* fe, ElementTransformation * T, int order)
+{
+   int  dof = fe->GetDof();
+   int  dim = fe->GetDim();
+
+   DenseMatrix pshape(dof, dim);
+   DenseMatrix mshape(dof, dim);
+   Vector pcomp;
+   Vector mcomp;
+   Vector dshape(dof);
+   Vector fdcomp(dof), fdshape(dof);
+
+   // Optimal step size for central difference
+   double h = std::cbrt(std::numeric_limits<double>::epsilon());
+   double inv2h = 0.5 / h;
+
+   // Error in the finite difference approximation of the derivative of a
+   // Legendre polynomial: P_n'''(1) h^2 / 6. Because we use shifted and scaled
+   // Legendre polynomials we need to increase these estimates by 2^3. We also
+   // make use of the fact that the third derivatives of Legendre polynomials
+   // are bounded by +/- (n+1)(n+2)(n+3)(n+4)(n+5)(n+6)/48.
+   double err_est = (order + 1) * (order + 2) * (order + 3) *
+                    (order + 4) * (order + 5) * (order + 6) * h * h / 36.0;
+
+   bool pyr = fe->GetGeomType() == Geometry::PYRAMID;
+
+   const IntegrationRule *ir = &IntRules.Get(fe->GetGeomType(), 2*order+dim-1);
+
+   IntegrationPoint ptp;
+   IntegrationPoint ptm;
+
+   int npoints = ir->GetNPoints();
+   for (int i=0; i < npoints; ++i)
+   {
+      // Get the current integration point from the integration rule
+      IntegrationPoint pt = ir->IntPoint(i);
+      fe->CalcDivShape(pt, dshape);
+
+      CAPTURE(pt.x, pt.y, pt.z);
+
+      fdshape = 0.0;
+      for (int d=0; d<dim; d++)
+      {
+         // Compute shifted integration points
+         switch (d)
+         {
+            case 0:
+               ptm.x = pt.x - h; ptm.y = pt.y; ptm.z = pt.z;
+               ptp.x = pt.x + h; ptp.y = pt.y; ptp.z = pt.z;
+               break;
+            case 1:
+               ptm.x = pt.x; ptm.y = pt.y - h; ptm.z = pt.z;
+               ptp.x = pt.x; ptp.y = pt.y + h; ptp.z = pt.z;
+               break;
+            case 2:
+               ptm.x = pt.x; ptm.y = pt.y; ptm.z = pt.z - h;
+               ptp.x = pt.x; ptp.y = pt.y; ptp.z = pt.z + h;
+               break;
+            default:
+               ptm = pt;
+               ptp = pt;
+         }
+
+         // Compute shape functions at the shifted points
+         fe->CalcVShape(ptm, mshape);
+         fe->CalcVShape(ptp, pshape);
+
+         // Extract the component to be differentiated
+         mshape.GetColumnReference(d, mcomp);
+         pshape.GetColumnReference(d, pcomp);
+
+         // Compute approximate derivatives using the secant method
+         add(inv2h, pcomp, -inv2h, mcomp, fdcomp);
+
+         fdshape += fdcomp;
+      }
+
+      // Compute the difference between the computed derivative and its
+      // finite difference approximation
+      fdshape -= dshape;
+
+      // Due to the scaling of the Legendre polynomials, as the integration
+      // points approach the apex of a pyramid the derivatives in the x and y
+      // directions become infinite. Therefore, we need to scale the finite
+      // difference error estimate by the following z-dependent factor.
+      double pyr_fac = pyr ? std::pow(1.0/(1.0-pt.z), 3) : 1.0;
+
+      // Determine the maximum difference between the two derivative
+      // calculations
+      double max_err = fdshape.Normlinf();
+
+      // The first factor of dim is added to account for the product
+      // rule used in computing derivatives of our basis functions which are
+      // products of Legendre polynomials in the different coordinates. The
+      // second factor of dim is added to account for the sum of derivatives
+      // in each direction needed to form the divergence.
+      REQUIRE( max_err < dim * dim * pyr_fac * err_est );
+   }
+}
+
+TEST_CASE("CalcDivShape vs FD RT",
+          "[RT_TriangleElement]"
+          "[RT_QuadrilateralElement]"
+          "[RT_TetrahedronElement]"
+          "[RT_WedgeElement]"
+          "[RT_FuentesPyramidElement]"
+          "[RT_HexahedronElement]")
+{
+   const int maxOrder = 5;
+   auto order = GENERATE_COPY(range(1, maxOrder + 1));
+
+   CAPTURE(order);
+
+   SECTION("RT_TriangleElement")
+   {
+      IsoparametricTransformation T;
+      GetReferenceTransformation(Element::TRIANGLE, T);
+
+      RT_TriangleElement fe(order - 1);
+      TestFDCalcDivShape(&fe, &T, order);
+   }
+
+   SECTION("RT_QuadrilateralElement")
+   {
+      IsoparametricTransformation T;
+      GetReferenceTransformation(Element::QUADRILATERAL, T);
+
+      RT_QuadrilateralElement fe(order - 1);
+      TestFDCalcDivShape(&fe, &T, order);
+   }
+
+   SECTION("RT_TetrahedronElement")
+   {
+      IsoparametricTransformation T;
+      GetReferenceTransformation(Element::TETRAHEDRON, T);
+
+      RT_TetrahedronElement fe(order - 1);
+      TestFDCalcDivShape(&fe, &T, order);
+   }
+
+   SECTION("RT_WedgeElement")
+   {
+      IsoparametricTransformation T;
+      GetReferenceTransformation(Element::WEDGE, T);
+
+      RT_WedgeElement fe(order - 1);
+      TestFDCalcDivShape(&fe, &T, order);
+   }
+
+   SECTION("RT_PyramidElement")
+   {
+      IsoparametricTransformation T;
+      GetReferenceTransformation(Element::PYRAMID, T);
+
+      if (order == 1)
+      {
+         RT0PyrFiniteElement fe;
+         TestFDCalcDivShape(&fe, &T, order);
+      }
+      else if (order == 2)
+      {
+         RT1PyrFiniteElement fe;
+         TestFDCalcDivShape(&fe, &T, order);
+      }
+   }
+
+   SECTION("RT_HexahedronElement")
+   {
+      IsoparametricTransformation T;
+      GetReferenceTransformation(Element::HEXAHEDRON, T);
+
+      RT_HexahedronElement fe(order - 1);
+      TestFDCalcDivShape(&fe, &T, order);
+   }
+}

tests/unit/fem/test_calcvshape.cpp

@@ -138,6 +138,25 @@ void GetReferenceTransformation(const Element::Type ElemType,
          T.GetPointMat()(2, 5) = 1.0;
          T.SetFE(&WedgeFE);
          break;
+      case Element::PYRAMID :
+         T.GetPointMat().SetSize(3, 5);
+         T.GetPointMat()(0, 0) = 0.0;
+         T.GetPointMat()(1, 0) = 0.0;
+         T.GetPointMat()(2, 0) = 0.0;
+         T.GetPointMat()(0, 1) = 1.0;
+         T.GetPointMat()(1, 1) = 0.0;
+         T.GetPointMat()(2, 1) = 0.0;
+         T.GetPointMat()(0, 2) = 1.0;
+         T.GetPointMat()(1, 2) = 1.0;
+         T.GetPointMat()(2, 2) = 0.0;
+         T.GetPointMat()(0, 3) = 0.0;
+         T.GetPointMat()(1, 3) = 1.0;
+         T.GetPointMat()(2, 3) = 0.0;
+         T.GetPointMat()(0, 4) = 0.0;
+         T.GetPointMat()(1, 4) = 0.0;
+         T.GetPointMat()(2, 4) = 1.0;
+         T.SetFE(&PyramidFE);
+         break;
       default:
          MFEM_ABORT("Unknown element type \"" << ElemType << "\"");
          break;