make style

2. Add _USE_MATH_DEFINES to CmakeLists.txt (without this change VS doesn't "see" this preprocessor in config.hpp)
3. Remove asm definition - asm is not allowed in VS anymore
4. Convert "const double Epsilon" to constexpr (fixes CUDA compilation error)

CHANGELOG

@@ -12,6 +12,13 @@ Version 4.3.1 (development)
 ===========================
 - Added support for hr-adaptivity using TMOP-based error estimator.
 
+- Adding lowest order Nedelec and Raviart-Thomas basis functions on wedge
+  shaped elements.
+
+- Added initial support for meshes with pyramidal elements, including several
+  pyramidal meshes in the data/ directory and support for the lowest order H1,
+  Nedelec, Raviart-Thomas, and L2 basis functions on pyramids.
+
 - Updated the hypre interface according to changes in hypre-2.22.1. The ADS
   solver is now fully working on GPUs.
 

CMakeLists.txt

@@ -175,6 +175,10 @@ else()
   set(MFEM_DEBUG OFF)
 endif()
 
+if (WIN32)
+	add_definitions(-D_USE_MATH_DEFINES)
+endif()
+
 # MPI -> hypre; PETSc (optional)
 if (MFEM_USE_MPI)
   find_package(MPI REQUIRED)

data/inline-pyramid.mesh

@@ -0,0 +1,9 @@
+MFEM INLINE mesh v1.0
+
+type = pyramid
+nx = 4
+ny = 4
+nz = 4
+sx = 1.0
+sy = 1.0
+sz = 1.0

data/octahedron.mesh

@@ -0,0 +1,43 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+# PRISM       = 6
+# PYRAMID     = 7
+#
+
+dimension
+3
+
+elements
+2
+1 7 4 3 2 1 0
+1 7 1 2 3 4 5
+
+boundary
+8
+1 2 0 2 1
+2 2 0 3 2
+3 2 0 4 3
+4 2 0 1 4
+5 2 1 2 5
+6 2 2 3 5
+7 2 3 4 5
+8 2 4 1 5
+
+vertices
+6
+3
+0 0 -1
+1 0 0
+0 1 0
+-1 0 0
+0 -1 0
+0 0 1

data/ref-pyramid.mesh

@@ -0,0 +1,38 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see mesh/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+# PRISM       = 6
+# PYRAMID     = 7
+#
+
+dimension
+3
+
+elements
+1
+1 7 0 1 2 3 4
+
+boundary
+5
+1 3 3 2 1 0
+2 2 0 1 4
+3 2 1 2 4
+4 2 2 3 4
+5 2 3 0 4
+
+vertices
+5
+3
+0 0 0
+1 0 0
+1 1 0
+0 1 0
+0 0 1

examples/ex1.cpp

@@ -9,6 +9,7 @@
 //               ex1 -m ../data/fichera.mesh
 //               ex1 -m ../data/fichera-mixed.mesh
 //               ex1 -m ../data/toroid-wedge.mesh
+//               ex1 -m ../data/octahedron.mesh -o 1
 //               ex1 -m ../data/periodic-annulus-sector.msh
 //               ex1 -m ../data/periodic-torus-sector.msh
 //               ex1 -m ../data/square-disc-p2.vtk -o 2

examples/ex1p.cpp

@@ -9,6 +9,7 @@
 //               mpirun -np 4 ex1p -m ../data/fichera.mesh
 //               mpirun -np 4 ex1p -m ../data/fichera-mixed.mesh
 //               mpirun -np 4 ex1p -m ../data/toroid-wedge.mesh
+//               mpirun -np 4 ex1p -m ../data/octahedron.mesh -o 1
 //               mpirun -np 4 ex1p -m ../data/periodic-annulus-sector.msh
 //               mpirun -np 4 ex1p -m ../data/periodic-torus-sector.msh
 //               mpirun -np 4 ex1p -m ../data/square-disc-p2.vtk -o 2

examples/ex22.cpp

@@ -13,6 +13,8 @@
 //               ex22 -m ../data/inline-hex.mesh -o 2 -p 1
 //               ex22 -m ../data/inline-hex.mesh -o 2 -p 2
 //               ex22 -m ../data/inline-hex.mesh -o 2 -p 2 -pa
+//               ex22 -m ../data/inline-wedge.mesh -o 1
+//               ex22 -m ../data/inline-pyramid.mesh -o 1
 //               ex22 -m ../data/star.mesh -r 1 -o 2 -sigma 10.0
 //
 // Device sample runs:

examples/ex22p.cpp

@@ -13,6 +13,8 @@
 //               mpirun -np 4 ex22p -m ../data/inline-hex.mesh -o 2 -p 1
 //               mpirun -np 4 ex22p -m ../data/inline-hex.mesh -o 2 -p 2
 //               mpirun -np 4 ex22p -m ../data/inline-hex.mesh -o 1 -p 2 -pa
+//               mpirun -np 4 ex22p -m ../data/inline-wedge.mesh -o 1
+//               mpirun -np 4 ex22p -m ../data/inline-pyramid.mesh -o 1
 //               mpirun -np 4 ex22p -m ../data/star.mesh -o 2 -sigma 10.0
 //
 // Device sample runs:

examples/ex25.cpp

@@ -92,7 +92,7 @@ class PMLDiagMatrixCoefficient : public VectorCoefficient
 {
 private:
    CartesianPML * pml = nullptr;
-   void (*Function)(const Vector &, CartesianPML * , Vector &);
+   void (*Function)(const Vector &, CartesianPML *, Vector &);
 public:
    PMLDiagMatrixCoefficient(int dim, void(*F)(const Vector &, CartesianPML *,
                                               Vector &),

examples/ex25p.cpp

@@ -92,7 +92,7 @@ class PMLDiagMatrixCoefficient : public VectorCoefficient
 {
 private:
    CartesianPML * pml = nullptr;
-   void (*Function)(const Vector &, CartesianPML * , Vector &);
+   void (*Function)(const Vector &, CartesianPML *, Vector &);
 public:
    PMLDiagMatrixCoefficient(int dim, void(*F)(const Vector &, CartesianPML *,
                                               Vector &),

examples/ex3.cpp

@@ -16,6 +16,8 @@
 //               ex3 -m ../data/beam-hex-nurbs.mesh
 //               ex3 -m ../data/amr-hex.mesh
 //               ex3 -m ../data/fichera-amr.mesh
+//               ex3 -m ../data/ref-prism.mesh -o 1
+//               ex3 -m ../data/octahedron.mesh -o 1
 //               ex3 -m ../data/star-surf.mesh -o 1
 //               ex3 -m ../data/mobius-strip.mesh -f 0.1
 //               ex3 -m ../data/klein-bottle.mesh -f 0.1

examples/ex3p.cpp

@@ -16,6 +16,8 @@
 //               mpirun -np 4 ex3p -m ../data/beam-hex-nurbs.mesh
 //               mpirun -np 4 ex3p -m ../data/amr-quad.mesh -o 2
 //               mpirun -np 4 ex3p -m ../data/amr-hex.mesh
+//               mpirun -np 4 ex3p -m ../data/ref-prism.mesh -o 1
+//               mpirun -np 4 ex3p -m ../data/octahedron.mesh -o 1
 //               mpirun -np 4 ex3p -m ../data/star-surf.mesh -o 2
 //               mpirun -np 4 ex3p -m ../data/mobius-strip.mesh -o 2 -f 0.1
 //               mpirun -np 4 ex3p -m ../data/klein-bottle.mesh -o 2 -f 0.1

examples/ex4.cpp

@@ -19,6 +19,8 @@
 //               ex4 -m ../data/amr-hex.mesh
 //               ex4 -m ../data/amr-hex.mesh -o 2 -hb
 //               ex4 -m ../data/fichera-amr.mesh -o 2 -sc
+//               ex4 -m ../data/ref-prism.mesh -o 1
+//               ex4 -m ../data/octahedron.mesh -o 1
 //               ex4 -m ../data/star-surf.mesh -o 1
 //
 // Device sample runs:

examples/ex4p.cpp

@@ -19,6 +19,8 @@
 //               mpirun -np 3 ex4p -m ../data/amr-quad.mesh -o 2 -hb
 //               mpirun -np 4 ex4p -m ../data/amr-hex.mesh -o 2 -sc
 //               mpirun -np 4 ex4p -m ../data/amr-hex.mesh -o 2 -hb
+//               mpirun -np 4 ex4p -m ../data/ref-prism.mesh -o 1
+//               mpirun -np 4 ex4p -m ../data/octahedron.mesh -o 1
 //               mpirun -np 4 ex4p -m ../data/star-surf.mesh -o 3 -hb
 //
 // Device sample runs:

fem/bilininteg_convection_ea.cpp

@@ -17,14 +17,14 @@ namespace mfem
 {
 
 template<int T_D1D = 0, int T_Q1D = 0>
-static void EAConvectionAssemble1D(const int NE,
-                                   const Array<double> &b,
-                                   const Array<double> &g,
-                                   const Vector &padata,
-                                   Vector &eadata,
-                                   const bool add,
-                                   const int d1d = 0,
-                                   const int q1d = 0)
+void EAConvectionAssemble1D(const int NE,
+                            const Array<double> &b,
+                            const Array<double> &g,
+                            const Vector &padata,
+                            Vector &eadata,
+                            const bool add,
+                            const int d1d = 0,
+                            const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -69,14 +69,14 @@ static void EAConvectionAssemble1D(const int NE,
 }
 
 template<int T_D1D = 0, int T_Q1D = 0>
-static void EAConvectionAssemble2D(const int NE,
-                                   const Array<double> &b,
-                                   const Array<double> &g,
-                                   const Vector &padata,
-                                   Vector &eadata,
-                                   const bool add,
-                                   const int d1d = 0,
-                                   const int q1d = 0)
+void EAConvectionAssemble2D(const int NE,
+                            const Array<double> &b,
+                            const Array<double> &g,
+                            const Vector &padata,
+                            Vector &eadata,
+                            const bool add,
+                            const int d1d = 0,
+                            const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -146,14 +146,14 @@ static void EAConvectionAssemble2D(const int NE,
 }
 
 template<int T_D1D = 0, int T_Q1D = 0>
-static void EAConvectionAssemble3D(const int NE,
-                                   const Array<double> &b,
-                                   const Array<double> &g,
-                                   const Vector &padata,
-                                   Vector &eadata,
-                                   const bool add,
-                                   const int d1d = 0,
-                                   const int q1d = 0)
+void EAConvectionAssemble3D(const int NE,
+                            const Array<double> &b,
+                            const Array<double> &g,
+                            const Vector &padata,
+                            Vector &eadata,
+                            const bool add,
+                            const int d1d = 0,
+                            const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;

fem/bilininteg_convection_pa.cpp

@@ -21,13 +21,13 @@ namespace mfem
 // PA Convection Integrator
 
 // PA Convection Assemble 2D kernel
-static void PAConvectionSetup2D(const int NQ,
-                                const int NE,
-                                const Array<double> &w,
-                                const Vector &j,
-                                const Vector &vel,
-                                const double alpha,
-                                Vector &op)
+void PAConvectionSetup2D(const int NQ,
+                         const int NE,
+                         const Array<double> &w,
+                         const Vector &j,
+                         const Vector &vel,
+                         const double alpha,
+                         Vector &op)
 {
    constexpr int DIM = 2;
 
@@ -60,13 +60,13 @@ static void PAConvectionSetup2D(const int NQ,
 }
 
 // PA Convection Assemble 3D kernel
-static void PAConvectionSetup3D(const int NQ,
-                                const int NE,
-                                const Array<double> &w,
-                                const Vector &j,
-                                const Vector &vel,
-                                const double alpha,
-                                Vector &op)
+void PAConvectionSetup3D(const int NQ,
+                         const int NE,
+                         const Array<double> &w,
+                         const Vector &j,
+                         const Vector &vel,
+                         const double alpha,
+                         Vector &op)
 {
    constexpr int DIM = 3;
    constexpr int SDIM = DIM;
@@ -135,7 +135,7 @@ static void PAConvectionSetup(const int dim,
 }
 
 // PA Convection Apply 2D kernel
-template<int T_D1D = 0, int T_Q1D = 0> static
+template<int T_D1D = 0, int T_Q1D = 0>
 void PAConvectionApply2D(const int ne,
                          const Array<double> &b,
                          const Array<double> &g,
@@ -254,7 +254,7 @@ void PAConvectionApply2D(const int ne,
 }
 
 // Optimized PA Convection Apply 2D kernel
-template<int T_D1D = 0, int T_Q1D = 0, int T_NBZ = 0> static
+template<int T_D1D = 0, int T_Q1D = 0, int T_NBZ = 0>
 void SmemPAConvectionApply2D(const int ne,
                              const Array<double> &b,
                              const Array<double> &g,
@@ -382,7 +382,7 @@ void SmemPAConvectionApply2D(const int ne,
 }
 
 // PA Convection Apply 3D kernel
-template<int T_D1D = 0, int T_Q1D = 0> static
+template<int T_D1D = 0, int T_Q1D = 0>
 void PAConvectionApply3D(const int ne,
                          const Array<double> &b,
                          const Array<double> &g,
@@ -563,7 +563,7 @@ void PAConvectionApply3D(const int ne,
 }
 
 // Optimized PA Convection Apply 3D kernel
-template<int T_D1D = 0, int T_Q1D = 0> static
+template<int T_D1D = 0, int T_Q1D = 0>
 void SmemPAConvectionApply3D(const int ne,
                              const Array<double> &b,
                              const Array<double> &g,

fem/bilininteg_dgtrace_ea.cpp

@@ -16,12 +16,12 @@
 namespace mfem
 {
 
-static void EADGTraceAssemble1DInt(const int NF,
-                                   const Array<double> &basis,
-                                   const Vector &padata,
-                                   Vector &eadata_int,
-                                   Vector &eadata_ext,
-                                   const bool add)
+void EADGTraceAssemble1DInt(const int NF,
+                            const Array<double> &basis,
+                            const Vector &padata,
+                            Vector &eadata_int,
+                            Vector &eadata_ext,
+                            const bool add)
 {
    auto D = Reshape(padata.Read(), 2, 2, NF);
    auto A_int = Reshape(eadata_int.ReadWrite(), 2, NF);
@@ -50,11 +50,11 @@ static void EADGTraceAssemble1DInt(const int NF,
    });
 }
 
-static void EADGTraceAssemble1DBdr(const int NF,
-                                   const Array<double> &basis,
-                                   const Vector &padata,
-                                   Vector &eadata_bdr,
-                                   const bool add)
+void EADGTraceAssemble1DBdr(const int NF,
+                            const Array<double> &basis,
+                            const Vector &padata,
+                            Vector &eadata_bdr,
+                            const bool add)
 {
    auto D = Reshape(padata.Read(), 2, 2, NF);
    auto A_bdr = Reshape(eadata_bdr.ReadWrite(), NF);
@@ -72,14 +72,14 @@ static void EADGTraceAssemble1DBdr(const int NF,
 }
 
 template<int T_D1D = 0, int T_Q1D = 0>
-static void EADGTraceAssemble2DInt(const int NF,
-                                   const Array<double> &basis,
-                                   const Vector &padata,
-                                   Vector &eadata_int,
-                                   Vector &eadata_ext,
-                                   const bool add,
-                                   const int d1d = 0,
-                                   const int q1d = 0)
+void EADGTraceAssemble2DInt(const int NF,
+                            const Array<double> &basis,
+                            const Vector &padata,
+                            Vector &eadata_int,
+                            Vector &eadata_ext,
+                            const bool add,
+                            const int d1d = 0,
+                            const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -128,13 +128,13 @@ static void EADGTraceAssemble2DInt(const int NF,
 }
 
 template<int T_D1D = 0, int T_Q1D = 0>
-static void EADGTraceAssemble2DBdr(const int NF,
-                                   const Array<double> &basis,
-                                   const Vector &padata,
-                                   Vector &eadata_bdr,
-                                   const bool add,
-                                   const int d1d = 0,
-                                   const int q1d = 0)
+void EADGTraceAssemble2DBdr(const int NF,
+                            const Array<double> &basis,
+                            const Vector &padata,
+                            Vector &eadata_bdr,
+                            const bool add,
+                            const int d1d = 0,
+                            const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -170,14 +170,14 @@ static void EADGTraceAssemble2DBdr(const int NF,
 }
 
 template<int T_D1D = 0, int T_Q1D = 0>
-static void EADGTraceAssemble3DInt(const int NF,
-                                   const Array<double> &basis,
-                                   const Vector &padata,
-                                   Vector &eadata_int,
-                                   Vector &eadata_ext,
-                                   const bool add,
-                                   const int d1d = 0,
-                                   const int q1d = 0)
+void EADGTraceAssemble3DInt(const int NF,
+                            const Array<double> &basis,
+                            const Vector &padata,
+                            Vector &eadata_int,
+                            Vector &eadata_ext,
+                            const bool add,
+                            const int d1d = 0,
+                            const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -268,13 +268,13 @@ static void EADGTraceAssemble3DInt(const int NF,
 }
 
 template<int T_D1D = 0, int T_Q1D = 0>
-static void EADGTraceAssemble3DBdr(const int NF,
-                                   const Array<double> &basis,
-                                   const Vector &padata,
-                                   Vector &eadata_bdr,
-                                   const bool add,
-                                   const int d1d = 0,
-                                   const int q1d = 0)
+void EADGTraceAssemble3DBdr(const int NF,
+                            const Array<double> &basis,
+                            const Vector &padata,
+                            Vector &eadata_bdr,
+                            const bool add,
+                            const int d1d = 0,
+                            const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;

fem/bilininteg_dgtrace_pa.cpp

@@ -19,16 +19,16 @@ using namespace std;
 namespace mfem
 {
 // PA DG Trace Integrator
-static void PADGTraceSetup2D(const int Q1D,
-                             const int NF,
-                             const Array<double> &w,
-                             const Vector &det,
-                             const Vector &nor,
-                             const Vector &rho,
-                             const Vector &vel,
-                             const double alpha,
-                             const double beta,
-                             Vector &op)
+void PADGTraceSetup2D(const int Q1D,
+                      const int NF,
+                      const Array<double> &w,
+                      const Vector &det,
+                      const Vector &nor,
+                      const Vector &rho,
+                      const Vector &vel,
+                      const double alpha,
+                      const double beta,
+                      Vector &op)
 {
    const int VDIM = 2;
 
@@ -61,16 +61,16 @@ static void PADGTraceSetup2D(const int Q1D,
    });
 }
 
-static void PADGTraceSetup3D(const int Q1D,
-                             const int NF,
-                             const Array<double> &w,
-                             const Vector &det,
-                             const Vector &nor,
-                             const Vector &rho,
-                             const Vector &vel,
-                             const double alpha,
-                             const double beta,
-                             Vector &op)
+void PADGTraceSetup3D(const int Q1D,
+                      const int NF,
+                      const Array<double> &w,
+                      const Vector &det,
+                      const Vector &nor,
+                      const Vector &rho,
+                      const Vector &vel,
+                      const double alpha,
+                      const double beta,
+                      Vector &op)
 {
    const int VDIM = 3;
 
@@ -301,7 +301,7 @@ void DGTraceIntegrator::AssemblePABoundaryFaces(const FiniteElementSpace& fes)
 }
 
 // PA DGTrace Apply 2D kernel for Gauss-Lobatto/Bernstein
-template<int T_D1D = 0, int T_Q1D = 0> static
+template<int T_D1D = 0, int T_Q1D = 0>
 void PADGTraceApply2D(const int NF,
                       const Array<double> &b,
                       const Array<double> &bt,
@@ -392,7 +392,7 @@ void PADGTraceApply2D(const int NF,
 }
 
 // PA DGTrace Apply 3D kernel for Gauss-Lobatto/Bernstein
-template<int T_D1D = 0, int T_Q1D = 0> static
+template<int T_D1D = 0, int T_Q1D = 0>
 void PADGTraceApply3D(const int NF,
                       const Array<double> &b,
                       const Array<double> &bt,
@@ -537,7 +537,7 @@ void PADGTraceApply3D(const int NF,
 }
 
 // Optimized PA DGTrace Apply 3D kernel for Gauss-Lobatto/Bernstein
-template<int T_D1D = 0, int T_Q1D = 0, int T_NBZ = 0> static
+template<int T_D1D = 0, int T_Q1D = 0, int T_NBZ = 0>
 void SmemPADGTraceApply3D(const int NF,
                           const Array<double> &b,
                           const Array<double> &bt,
@@ -701,7 +701,7 @@ static void PADGTraceApply(const int dim,
 }
 
 // PA DGTrace Apply 2D kernel for Gauss-Lobatto/Bernstein
-template<int T_D1D = 0, int T_Q1D = 0> static
+template<int T_D1D = 0, int T_Q1D = 0>
 void PADGTraceApplyTranspose2D(const int NF,
                                const Array<double> &b,
                                const Array<double> &bt,
@@ -797,7 +797,7 @@ void PADGTraceApplyTranspose2D(const int NF,
 }
 
 // PA DGTrace Apply Transpose 3D kernel for Gauss-Lobatto/Bernstein
-template<int T_D1D = 0, int T_Q1D = 0> static
+template<int T_D1D = 0, int T_Q1D = 0>
 void PADGTraceApplyTranspose3D(const int NF,
                                const Array<double> &b,
                                const Array<double> &bt,
@@ -953,7 +953,7 @@ void PADGTraceApplyTranspose3D(const int NF,
 }
 
 // Optimized PA DGTrace Apply Transpose 3D kernel for Gauss-Lobatto/Bernstein
-template<int T_D1D = 0, int T_Q1D = 0, int T_NBZ = 0> static
+template<int T_D1D = 0, int T_Q1D = 0, int T_NBZ = 0>
 void SmemPADGTraceApplyTranspose3D(const int NF,
                                    const Array<double> &b,
                                    const Array<double> &bt,

fem/bilininteg_diffusion_ea.cpp

@@ -17,14 +17,14 @@ namespace mfem
 {
 
 template<int T_D1D = 0, int T_Q1D = 0>
-static void EADiffusionAssemble1D(const int NE,
-                                  const Array<double> &b,
-                                  const Array<double> &g,
-                                  const Vector &padata,
-                                  Vector &eadata,
-                                  const bool add,
-                                  const int d1d = 0,
-                                  const int q1d = 0)
+void EADiffusionAssemble1D(const int NE,
+                           const Array<double> &b,
+                           const Array<double> &g,
+                           const Vector &padata,
+                           Vector &eadata,
+                           const bool add,
+                           const int d1d = 0,
+                           const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -68,14 +68,14 @@ static void EADiffusionAssemble1D(const int NE,
 }
 
 template<int T_D1D = 0, int T_Q1D = 0>
-static void EADiffusionAssemble2D(const int NE,
-                                  const Array<double> &b,
-                                  const Array<double> &g,
-                                  const Vector &padata,
-                                  Vector &eadata,
-                                  const bool add,
-                                  const int d1d = 0,
-                                  const int q1d = 0)
+void EADiffusionAssemble2D(const int NE,
+                           const Array<double> &b,
+                           const Array<double> &g,
+                           const Vector &padata,
+                           Vector &eadata,
+                           const bool add,
+                           const int d1d = 0,
+                           const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -145,14 +145,14 @@ static void EADiffusionAssemble2D(const int NE,
 }
 
 template<int T_D1D = 0, int T_Q1D = 0>
-static void EADiffusionAssemble3D(const int NE,
-                                  const Array<double> &b,
-                                  const Array<double> &g,
-                                  const Vector &padata,
-                                  Vector &eadata,
-                                  const bool add,
-                                  const int d1d = 0,
-                                  const int q1d = 0)
+void EADiffusionAssemble3D(const int NE,
+                           const Array<double> &b,
+                           const Array<double> &g,
+                           const Vector &padata,
+                           Vector &eadata,
+                           const bool add,
+                           const int d1d = 0,
+                           const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;

fem/bilininteg_diffusion_pa.cpp

@@ -496,14 +496,14 @@ void DiffusionIntegrator::AssemblePA(const FiniteElementSpace &fes)
 }
 
 template<int T_D1D = 0, int T_Q1D = 0>
-static void PADiffusionDiagonal2D(const int NE,
-                                  const bool symmetric,
-                                  const Array<double> &b,
-                                  const Array<double> &g,
-                                  const Vector &d,
-                                  Vector &y,
-                                  const int d1d = 0,
-                                  const int q1d = 0)
+void PADiffusionDiagonal2D(const int NE,
+                           const bool symmetric,
+                           const Array<double> &b,
+                           const Array<double> &g,
+                           const Vector &d,
+                           Vector &y,
+                           const int d1d = 0,
+                           const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -562,14 +562,14 @@ static void PADiffusionDiagonal2D(const int NE,
 
 // Shared memory PA Diffusion Diagonal 2D kernel
 template<int T_D1D = 0, int T_Q1D = 0, int T_NBZ = 0>
-static void SmemPADiffusionDiagonal2D(const int NE,
-                                      const bool symmetric,
-                                      const Array<double> &b_,
-                                      const Array<double> &g_,
-                                      const Vector &d_,
-                                      Vector &y_,
-                                      const int d1d = 0,
-                                      const int q1d = 0)
+void SmemPADiffusionDiagonal2D(const int NE,
+                               const bool symmetric,
+                               const Array<double> &b_,
+                               const Array<double> &g_,
+                               const Vector &d_,
+                               Vector &y_,
+                               const int d1d = 0,
+                               const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -656,14 +656,14 @@ static void SmemPADiffusionDiagonal2D(const int NE,
 }
 
 template<int T_D1D = 0, int T_Q1D = 0>
-static void PADiffusionDiagonal3D(const int NE,
-                                  const bool symmetric,
-                                  const Array<double> &b,
-                                  const Array<double> &g,
-                                  const Vector &d,
-                                  Vector &y,
-                                  const int d1d = 0,
-                                  const int q1d = 0)
+void PADiffusionDiagonal3D(const int NE,
+                           const bool symmetric,
+                           const Array<double> &b,
+                           const Array<double> &g,
+                           const Vector &d,
+                           Vector &y,
+                           const int d1d = 0,
+                           const int q1d = 0)
 {
    constexpr int DIM = 3;
    const int D1D = T_D1D ? T_D1D : d1d;
@@ -757,14 +757,14 @@ static void PADiffusionDiagonal3D(const int NE,
 
 // Shared memory PA Diffusion Diagonal 3D kernel
 template<int T_D1D = 0, int T_Q1D = 0>
-static void SmemPADiffusionDiagonal3D(const int NE,
-                                      const bool symmetric,
-                                      const Array<double> &b_,
-                                      const Array<double> &g_,
-                                      const Vector &d_,
-                                      Vector &y_,
-                                      const int d1d = 0,
-                                      const int q1d = 0)
+void SmemPADiffusionDiagonal3D(const int NE,
+                               const bool symmetric,
+                               const Array<double> &b_,
+                               const Array<double> &g_,
+                               const Vector &d_,
+                               Vector &y_,
+                               const int d1d = 0,
+                               const int q1d = 0)
 {
    constexpr int DIM = 3;
    const int D1D = T_D1D ? T_D1D : d1d;
@@ -1034,17 +1034,17 @@ static void OccaPADiffusionApply3D(const int D1D,
 
 // PA Diffusion Apply 2D kernel
 template<int T_D1D = 0, int T_Q1D = 0>
-static void PADiffusionApply2D(const int NE,
-                               const bool symmetric,
-                               const Array<double> &b_,
-                               const Array<double> &g_,
-                               const Array<double> &bt_,
-                               const Array<double> &gt_,
-                               const Vector &d_,
-                               const Vector &x_,
-                               Vector &y_,
-                               const int d1d = 0,
-                               const int q1d = 0)
+void PADiffusionApply2D(const int NE,
+                        const bool symmetric,
+                        const Array<double> &b_,
+                        const Array<double> &g_,
+                        const Array<double> &bt_,
+                        const Array<double> &gt_,
+                        const Vector &d_,
+                        const Vector &x_,
+                        Vector &y_,
+                        const int d1d = 0,
+                        const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -1156,15 +1156,15 @@ static void PADiffusionApply2D(const int NE,
 
 // Shared memory PA Diffusion Apply 2D kernel
 template<int T_D1D = 0, int T_Q1D = 0, int T_NBZ = 0>
-static void SmemPADiffusionApply2D(const int NE,
-                                   const bool symmetric,
-                                   const Array<double> &b_,
-                                   const Array<double> &g_,
-                                   const Vector &d_,
-                                   const Vector &x_,
-                                   Vector &y_,
-                                   const int d1d = 0,
-                                   const int q1d = 0)
+void SmemPADiffusionApply2D(const int NE,
+                            const bool symmetric,
+                            const Array<double> &b_,
+                            const Array<double> &g_,
+                            const Vector &d_,
+                            const Vector &x_,
+                            Vector &y_,
+                            const int d1d = 0,
+                            const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -1314,16 +1314,16 @@ static void SmemPADiffusionApply2D(const int NE,
 
 // PA Diffusion Apply 3D kernel
 template<int T_D1D = 0, int T_Q1D = 0>
-static void PADiffusionApply3D(const int NE,
-                               const bool symmetric,
-                               const Array<double> &b,
-                               const Array<double> &g,
-                               const Array<double> &bt,
-                               const Array<double> &gt,
-                               const Vector &d_,
-                               const Vector &x_,
-                               Vector &y_,
-                               int d1d = 0, int q1d = 0)
+void PADiffusionApply3D(const int NE,
+                        const bool symmetric,
+                        const Array<double> &b,
+                        const Array<double> &g,
+                        const Array<double> &bt,
+                        const Array<double> &gt,
+                        const Vector &d_,
+                        const Vector &x_,
+                        Vector &y_,
+                        int d1d = 0, int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -1533,15 +1533,15 @@ static MFEM_HOST_DEVICE inline double sign(const int q, const int d)
 }
 
 template<int T_D1D = 0, int T_Q1D = 0>
-static void SmemPADiffusionApply3D(const int NE,
-                                   const bool symmetric,
-                                   const Array<double> &b_,
-                                   const Array<double> &g_,
-                                   const Vector &d_,
-                                   const Vector &x_,
-                                   Vector &y_,
-                                   const int d1d = 0,
-                                   const int q1d = 0)
+void SmemPADiffusionApply3D(const int NE,
+                            const bool symmetric,
+                            const Array<double> &b_,
+                            const Array<double> &g_,
+                            const Vector &d_,
+                            const Vector &x_,
+                            Vector &y_,
+                            const int d1d = 0,
+                            const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;

fem/bilininteg_divergence.cpp

@@ -21,12 +21,12 @@ namespace mfem
 // PA Divergence Integrator
 
 // PA Divergence Assemble 2D kernel
-static void PADivergenceSetup2D(const int Q1D,
-                                const int NE,
-                                const Array<double> &w,
-                                const Vector &j,
-                                const double COEFF,
-                                Vector &op)
+void PADivergenceSetup2D(const int Q1D,
+                         const int NE,
+                         const Array<double> &w,
+                         const Vector &j,
+                         const double COEFF,
+                         Vector &op)
 {
    const int NQ = Q1D*Q1D;
    auto W = w.Read();
@@ -51,12 +51,12 @@ static void PADivergenceSetup2D(const int Q1D,
 }
 
 // PA Divergence Assemble 3D kernel
-static void PADivergenceSetup3D(const int Q1D,
-                                const int NE,
-                                const Array<double> &w,
-                                const Vector &j,
-                                const double COEFF,
-                                Vector &op)
+void PADivergenceSetup3D(const int Q1D,
+                         const int NE,
+                         const Array<double> &w,
+                         const Vector &j,
+                         const double COEFF,
+                         Vector &op)
 {
    const int NQ = Q1D*Q1D*Q1D;
    auto W = w.Read();
@@ -160,16 +160,16 @@ void VectorDivergenceIntegrator::AssemblePA(const FiniteElementSpace &trial_fes,
 
 // PA Divergence Apply 2D kernel
 template<const int T_TR_D1D = 0, const int T_TE_D1D = 0, const int T_Q1D = 0>
-static void PADivergenceApply2D(const int NE,
-                                const Array<double> &b,
-                                const Array<double> &g,
-                                const Array<double> &bt,
-                                const Vector &op_,
-                                const Vector &x_,
-                                Vector &y_,
-                                const int tr_d1d = 0,
-                                const int te_d1d = 0,
-                                const int q1d = 0)
+void PADivergenceApply2D(const int NE,
+                         const Array<double> &b,
+                         const Array<double> &g,
+                         const Array<double> &bt,
+                         const Vector &op_,
+                         const Vector &x_,
+                         Vector &y_,
+                         const int tr_d1d = 0,
+                         const int te_d1d = 0,
+                         const int q1d = 0)
 {
    const int TR_D1D = T_TR_D1D ? T_TR_D1D : tr_d1d;
    const int TE_D1D = T_TE_D1D ? T_TE_D1D : te_d1d;
@@ -281,16 +281,16 @@ static void PADivergenceApply2D(const int NE,
 // Shared memory PA Divergence Apply 2D kernel
 template<const int T_TR_D1D = 0, const int T_TE_D1D = 0, const int T_Q1D = 0,
          const int T_NBZ = 0>
-static void SmemPADivergenceApply2D(const int NE,
-                                    const Array<double> &b_,
-                                    const Array<double> &g_,
-                                    const Array<double> &bt_,
-                                    const Vector &op_,
-                                    const Vector &x_,
-                                    Vector &y_,
-                                    const int tr_d1d = 0,
-                                    const int te_d1d = 0,
-                                    const int q1d = 0)
+void SmemPADivergenceApply2D(const int NE,
+                             const Array<double> &b_,
+                             const Array<double> &g_,
+                             const Array<double> &bt_,
+                             const Vector &op_,
+                             const Vector &x_,
+                             Vector &y_,
+                             const int tr_d1d = 0,
+                             const int te_d1d = 0,
+                             const int q1d = 0)
 {
    // TODO
    MFEM_ASSERT(false, "SHARED MEM NOT PROGRAMMED YET");
@@ -298,16 +298,16 @@ static void SmemPADivergenceApply2D(const int NE,
 
 // PA Divergence Apply 2D kernel transpose
 template<const int T_TR_D1D = 0, const int T_TE_D1D = 0, const int T_Q1D = 0>
-static void PADivergenceApplyTranspose2D(const int NE,
-                                         const Array<double> &bt,
-                                         const Array<double> &gt,
-                                         const Array<double> &b,
-                                         const Vector &op_,
-                                         const Vector &x_,
-                                         Vector &y_,
-                                         const int tr_d1d = 0,
-                                         const int te_d1d = 0,
-                                         const int q1d = 0)
+void PADivergenceApplyTranspose2D(const int NE,
+                                  const Array<double> &bt,
+                                  const Array<double> &gt,
+                                  const Array<double> &b,
+                                  const Vector &op_,
+                                  const Vector &x_,
+                                  Vector &y_,
+                                  const int tr_d1d = 0,
+                                  const int te_d1d = 0,
+                                  const int q1d = 0)
 {
    const int TR_D1D = T_TR_D1D ? T_TR_D1D : tr_d1d;
    const int TE_D1D = T_TE_D1D ? T_TE_D1D : te_d1d;
@@ -414,16 +414,16 @@ static void PADivergenceApplyTranspose2D(const int NE,
 
 // PA Vector Divergence Apply 3D kernel
 template<const int T_TR_D1D = 0, const int T_TE_D1D = 0, const int T_Q1D = 0>
-static void PADivergenceApply3D(const int NE,
-                                const Array<double> &b,
-                                const Array<double> &g,
-                                const Array<double> &bt,
-                                const Vector &op_,
-                                const Vector &x_,
-                                Vector &y_,
-                                int tr_d1d = 0,
-                                int te_d1d = 0,
-                                int q1d = 0)
+void PADivergenceApply3D(const int NE,
+                         const Array<double> &b,
+                         const Array<double> &g,
+                         const Array<double> &bt,
+                         const Vector &op_,
+                         const Vector &x_,
+                         Vector &y_,
+                         int tr_d1d = 0,
+                         int te_d1d = 0,
+                         int q1d = 0)
 {
    const int TR_D1D = T_TR_D1D ? T_TR_D1D : tr_d1d;
    const int TE_D1D = T_TE_D1D ? T_TE_D1D : te_d1d;
@@ -597,16 +597,16 @@ static void PADivergenceApply3D(const int NE,
 
 // PA Vector Divergence Apply 3D kernel
 template<const int T_TR_D1D = 0, const int T_TE_D1D = 0, const int T_Q1D = 0>
-static void PADivergenceApplyTranspose3D(const int NE,
-                                         const Array<double> &bt,
-                                         const Array<double> &gt,
-                                         const Array<double> &b,
-                                         const Vector &op_,
-                                         const Vector &x_,
-                                         Vector &y_,
-                                         int tr_d1d = 0,
-                                         int te_d1d = 0,
-                                         int q1d = 0)
+void PADivergenceApplyTranspose3D(const int NE,
+                                  const Array<double> &bt,
+                                  const Array<double> &gt,
+                                  const Array<double> &b,
+                                  const Vector &op_,
+                                  const Vector &x_,
+                                  Vector &y_,
+                                  int tr_d1d = 0,
+                                  int te_d1d = 0,
+                                  int q1d = 0)
 {
    const int TR_D1D = T_TR_D1D ? T_TR_D1D : tr_d1d;
    const int TE_D1D = T_TE_D1D ? T_TE_D1D : te_d1d;
@@ -775,16 +775,16 @@ static void PADivergenceApplyTranspose3D(const int NE,
 
 // Shared memory PA Vector Divergence Apply 3D kernel
 template<const int T_TR_D1D = 0, const int T_TE_D1D = 0, const int T_Q1D = 0>
-static void SmemPADivergenceApply3D(const int NE,
-                                    const Array<double> &b_,
-                                    const Array<double> &g_,
-                                    const Array<double> &bt_,
-                                    const Vector &q_,
-                                    const Vector &x_,
-                                    Vector &y_,
-                                    const int tr_d1d = 0,
-                                    const int te_d1d = 0,
-                                    const int q1d = 0)
+void SmemPADivergenceApply3D(const int NE,
+                             const Array<double> &b_,
+                             const Array<double> &g_,
+                             const Array<double> &bt_,
+                             const Vector &q_,
+                             const Vector &x_,
+                             Vector &y_,
+                             const int tr_d1d = 0,
+                             const int te_d1d = 0,
+                             const int q1d = 0)
 {
    const int TR_D1D = T_TR_D1D ? T_TR_D1D : tr_d1d;
    const int TE_D1D = T_TE_D1D ? T_TE_D1D : te_d1d;

fem/bilininteg_gradient.cpp

@@ -70,12 +70,12 @@ namespace mfem
    the \b MFEM_SHARED keyword for local arrays. */
 
 // PA Gradient Assemble 2D kernel
-static void PAGradientSetup2D(const int Q1D,
-                              const int NE,
-                              const Array<double> &w,
-                              const Vector &j,
-                              const Vector &c,
-                              Vector &op)
+void PAGradientSetup2D(const int Q1D,
+                       const int NE,
+                       const Array<double> &w,
+                       const Vector &j,
+                       const Vector &c,
+                       Vector &op)
 {
    const int NQ = Q1D*Q1D;
    auto W = w.Read();
@@ -105,12 +105,12 @@ static void PAGradientSetup2D(const int Q1D,
 }
 
 // PA Gradient Assemble 3D kernel
-static void PAGradientSetup3D(const int Q1D,
-                              const int NE,
-                              const Array<double> &w,
-                              const Vector &j,
-                              const Vector &c,
-                              Vector &op)
+void PAGradientSetup3D(const int Q1D,
+                       const int NE,
+                       const Array<double> &w,
+                       const Vector &j,
+                       const Vector &c,
+                       Vector &op)
 {
    const int NQ = Q1D*Q1D*Q1D;
    auto W = w.Read();
@@ -254,16 +254,16 @@ void GradientIntegrator::AssemblePA(const FiniteElementSpace &trial_fes,
 
 // PA Gradient Apply 2D kernel
 template<int T_TR_D1D = 0, int T_TE_D1D = 0, int T_Q1D = 0>
-static void PAGradientApply2D(const int NE,
-                              const Array<double> &b,
-                              const Array<double> &g,
-                              const Array<double> &bt,
-                              const Vector &op_,
-                              const Vector &x_,
-                              Vector &y_,
-                              const int tr_d1d = 0,
-                              const int te_d1d = 0,
-                              const int q1d = 0)
+void PAGradientApply2D(const int NE,
+                       const Array<double> &b,
+                       const Array<double> &g,
+                       const Array<double> &bt,
+                       const Vector &op_,
+                       const Vector &x_,
+                       Vector &y_,
+                       const int tr_d1d = 0,
+                       const int te_d1d = 0,
+                       const int q1d = 0)
 {
    const int TR_D1D = T_TR_D1D ? T_TR_D1D : tr_d1d;
    const int TE_D1D = T_TE_D1D ? T_TE_D1D : te_d1d;
@@ -384,16 +384,16 @@ static void PAGradientApplyTranspose2D(const int NE,
 
 // PA Gradient Apply 3D kernel
 template<const int T_TR_D1D = 0, const int T_TE_D1D = 0, const int T_Q1D = 0>
-static void PAGradientApply3D(const int NE,
-                              const Array<double> &b,
-                              const Array<double> &g,
-                              const Array<double> &bt,
-                              const Vector &op_,
-                              const Vector &x_,
-                              Vector &y_,
-                              int tr_d1d = 0,
-                              int te_d1d = 0,
-                              int q1d = 0)
+void PAGradientApply3D(const int NE,
+                       const Array<double> &b,
+                       const Array<double> &g,
+                       const Array<double> &bt,
+                       const Vector &op_,
+                       const Vector &x_,
+                       Vector &y_,
+                       int tr_d1d = 0,
+                       int te_d1d = 0,
+                       int q1d = 0)
 {
    const int TR_D1D = T_TR_D1D ? T_TR_D1D : tr_d1d;
    const int TE_D1D = T_TE_D1D ? T_TE_D1D : te_d1d;
@@ -579,16 +579,16 @@ static void PAGradientApplyTranspose3D(const int NE,
 
 // Shared memory PA Gradient Apply 3D kernel
 template<const int T_TR_D1D = 0, const int T_TE_D1D = 0, const int T_Q1D = 0>
-static void SmemPAGradientApply3D(const int NE,
-                                  const Array<double> &b_,
-                                  const Array<double> &g_,
-                                  const Array<double> &bt_,
-                                  const Vector &d_,
-                                  const Vector &x_,
-                                  Vector &y_,
-                                  const int tr_d1d = 0,
-                                  const int te_d1d = 0,
-                                  const int q1d = 0)
+void SmemPAGradientApply3D(const int NE,
+                           const Array<double> &b_,
+                           const Array<double> &g_,
+                           const Array<double> &bt_,
+                           const Vector &d_,
+                           const Vector &x_,
+                           Vector &y_,
+                           const int tr_d1d = 0,
+                           const int te_d1d = 0,
+                           const int q1d = 0)
 {
    const int TR_D1D = T_TR_D1D ? T_TR_D1D : tr_d1d;
    const int TE_D1D = T_TE_D1D ? T_TE_D1D : te_d1d;

fem/bilininteg_hcurl.cpp

@@ -791,12 +791,12 @@ void SmemPAHcurlMassApply3D(const int D1D,
 }
 
 // PA H(curl) curl-curl assemble 2D kernel
-static void PACurlCurlSetup2D(const int Q1D,
-                              const int NE,
-                              const Array<double> &w,
-                              const Vector &j,
-                              Vector &coeff,
-                              Vector &op)
+void PACurlCurlSetup2D(const int Q1D,
+                       const int NE,
+                       const Array<double> &w,
+                       const Vector &j,
+                       Vector &coeff,
+                       Vector &op)
 {
    const int NQ = Q1D*Q1D;
    auto W = w.Read();
@@ -818,13 +818,13 @@ static void PACurlCurlSetup2D(const int Q1D,
 }
 
 // PA H(curl) curl-curl assemble 3D kernel
-static void PACurlCurlSetup3D(const int Q1D,
-                              const int coeffDim,
-                              const int NE,
-                              const Array<double> &w,
-                              const Vector &j,
-                              Vector &coeff,
-                              Vector &op)
+void PACurlCurlSetup3D(const int Q1D,
+                       const int coeffDim,
+                       const int NE,
+                       const Array<double> &w,
+                       const Vector &j,
+                       Vector &coeff,
+                       Vector &op)
 {
    const int NQ = Q1D*Q1D*Q1D;
    const bool symmetric = (coeffDim != 9);
@@ -1045,16 +1045,16 @@ void CurlCurlIntegrator::AssemblePA(const FiniteElementSpace &fes)
    }
 }
 
-static void PACurlCurlApply2D(const int D1D,
-                              const int Q1D,
-                              const int NE,
-                              const Array<double> &bo,
-                              const Array<double> &bot,
-                              const Array<double> &gc,
-                              const Array<double> &gct,
-                              const Vector &pa_data,
-                              const Vector &x,
-                              Vector &y)
+void PACurlCurlApply2D(const int D1D,
+                       const int Q1D,
+                       const int NE,
+                       const Array<double> &bo,
+                       const Array<double> &bot,
+                       const Array<double> &gc,
+                       const Array<double> &gct,
+                       const Vector &pa_data,
+                       const Vector &x,
+                       Vector &y)
 {
    constexpr static int VDIM = 2;
    constexpr static int MAX_D1D = HCURL_MAX_D1D;
@@ -1166,19 +1166,19 @@ static void PACurlCurlApply2D(const int D1D,
 }
 
 template<int MAX_D1D = HCURL_MAX_D1D, int MAX_Q1D = HCURL_MAX_Q1D>
-static void PACurlCurlApply3D(const int D1D,
-                              const int Q1D,
-                              const bool symmetric,
-                              const int NE,
-                              const Array<double> &bo,
-                              const Array<double> &bc,
-                              const Array<double> &bot,
-                              const Array<double> &bct,
-                              const Array<double> &gc,
-                              const Array<double> &gct,
-                              const Vector &pa_data,
-                              const Vector &x,
-                              Vector &y)
+void PACurlCurlApply3D(const int D1D,
+                       const int Q1D,
+                       const bool symmetric,
+                       const int NE,
+                       const Array<double> &bo,
+                       const Array<double> &bc,
+                       const Array<double> &bot,
+                       const Array<double> &bct,
+                       const Array<double> &gc,
+                       const Array<double> &gct,
+                       const Vector &pa_data,
+                       const Vector &x,
+                       Vector &y)
 {
    MFEM_VERIFY(D1D <= MAX_D1D, "Error: D1D > MAX_D1D");
    MFEM_VERIFY(Q1D <= MAX_Q1D, "Error: Q1D > MAX_Q1D");
@@ -1677,19 +1677,19 @@ static void PACurlCurlApply3D(const int D1D,
 }
 
 template<int MAX_D1D = HCURL_MAX_D1D, int MAX_Q1D = HCURL_MAX_Q1D>
-static void SmemPACurlCurlApply3D(const int D1D,
-                                  const int Q1D,
-                                  const bool symmetric,
-                                  const int NE,
-                                  const Array<double> &bo,
-                                  const Array<double> &bc,
-                                  const Array<double> &bot,
-                                  const Array<double> &bct,
-                                  const Array<double> &gc,
-                                  const Array<double> &gct,
-                                  const Vector &pa_data,
-                                  const Vector &x,
-                                  Vector &y)
+void SmemPACurlCurlApply3D(const int D1D,
+                           const int Q1D,
+                           const bool symmetric,
+                           const int NE,
+                           const Array<double> &bo,
+                           const Array<double> &bc,
+                           const Array<double> &bot,
+                           const Array<double> &bct,
+                           const Array<double> &gc,
+                           const Array<double> &gct,
+                           const Vector &pa_data,
+                           const Vector &x,
+                           Vector &y)
 {
    MFEM_VERIFY(D1D <= MAX_D1D, "Error: D1D > MAX_D1D");
    MFEM_VERIFY(Q1D <= MAX_Q1D, "Error: Q1D > MAX_Q1D");
@@ -2032,13 +2032,13 @@ void CurlCurlIntegrator::AddMultPA(const Vector &x, Vector &y) const
    }
 }
 
-static void PACurlCurlAssembleDiagonal2D(const int D1D,
-                                         const int Q1D,
-                                         const int NE,
-                                         const Array<double> &bo,
-                                         const Array<double> &gc,
-                                         const Vector &pa_data,
-                                         Vector &diag)
+void PACurlCurlAssembleDiagonal2D(const int D1D,
+                                  const int Q1D,
+                                  const int NE,
+                                  const Array<double> &bo,
+                                  const Array<double> &gc,
+                                  const Vector &pa_data,
+                                  Vector &diag)
 {
    constexpr static int VDIM = 2;
    constexpr static int MAX_Q1D = HCURL_MAX_Q1D;
@@ -2087,16 +2087,16 @@ static void PACurlCurlAssembleDiagonal2D(const int D1D,
 }
 
 template<int MAX_D1D = HCURL_MAX_D1D, int MAX_Q1D = HCURL_MAX_Q1D>
-static void PACurlCurlAssembleDiagonal3D(const int D1D,
-                                         const int Q1D,
-                                         const bool symmetric,
-                                         const int NE,
-                                         const Array<double> &bo,
-                                         const Array<double> &bc,
-                                         const Array<double> &go,
-                                         const Array<double> &gc,
-                                         const Vector &pa_data,
-                                         Vector &diag)
+void PACurlCurlAssembleDiagonal3D(const int D1D,
+                                  const int Q1D,
+                                  const bool symmetric,
+                                  const int NE,
+                                  const Array<double> &bo,
+                                  const Array<double> &bc,
+                                  const Array<double> &go,
+                                  const Array<double> &gc,
+                                  const Vector &pa_data,
+                                  Vector &diag)
 {
    constexpr static int VDIM = 3;
    MFEM_VERIFY(D1D <= MAX_D1D, "Error: D1D > MAX_D1D");
@@ -2273,16 +2273,16 @@ static void PACurlCurlAssembleDiagonal3D(const int D1D,
 }
 
 template<int MAX_D1D = HCURL_MAX_D1D, int MAX_Q1D = HCURL_MAX_Q1D>
-static void SmemPACurlCurlAssembleDiagonal3D(const int D1D,
-                                             const int Q1D,
-                                             const bool symmetric,
-                                             const int NE,
-                                             const Array<double> &bo,
-                                             const Array<double> &bc,
-                                             const Array<double> &go,
-                                             const Array<double> &gc,
-                                             const Vector &pa_data,
-                                             Vector &diag)
+void SmemPACurlCurlAssembleDiagonal3D(const int D1D,
+                                      const int Q1D,
+                                      const bool symmetric,
+                                      const int NE,
+                                      const Array<double> &bo,
+                                      const Array<double> &bc,
+                                      const Array<double> &go,
+                                      const Array<double> &gc,
+                                      const Vector &pa_data,
+                                      Vector &diag)
 {
    MFEM_VERIFY(D1D <= MAX_D1D, "Error: D1D > MAX_D1D");
    MFEM_VERIFY(Q1D <= MAX_Q1D, "Error: Q1D > MAX_Q1D");
@@ -2955,18 +2955,18 @@ void MixedVectorCurlIntegrator::AssemblePA(const FiniteElementSpace &trial_fes,
 // Apply to x corresponding to DOF's in H(curl) (trial), whose curl is
 // integrated against H(curl) test functions corresponding to y.
 template<int MAX_D1D = HCURL_MAX_D1D, int MAX_Q1D = HCURL_MAX_Q1D>
-static void PAHcurlL2Apply3D(const int D1D,
-                             const int Q1D,
-                             const int coeffDim,
-                             const int NE,
-                             const Array<double> &bo,
-                             const Array<double> &bc,
-                             const Array<double> &bot,
-                             const Array<double> &bct,
-                             const Array<double> &gc,
-                             const Vector &pa_data,
-                             const Vector &x,
-                             Vector &y)
+void PAHcurlL2Apply3D(const int D1D,
+                      const int Q1D,
+                      const int coeffDim,
+                      const int NE,
+                      const Array<double> &bo,
+                      const Array<double> &bc,
+                      const Array<double> &bot,
+                      const Array<double> &bct,
+                      const Array<double> &gc,
+                      const Vector &pa_data,
+                      const Vector &x,
+                      Vector &y)
 {
    MFEM_VERIFY(D1D <= MAX_D1D, "Error: D1D > MAX_D1D");
    MFEM_VERIFY(Q1D <= MAX_Q1D, "Error: Q1D > MAX_Q1D");
@@ -3297,16 +3297,16 @@ static void PAHcurlL2Apply3D(const int D1D,
 // Apply to x corresponding to DOF's in H(curl) (trial), whose curl is
 // integrated against H(curl) test functions corresponding to y.
 template<int MAX_D1D = HCURL_MAX_D1D, int MAX_Q1D = HCURL_MAX_Q1D>
-static void SmemPAHcurlL2Apply3D(const int D1D,
-                                 const int Q1D,
-                                 const int coeffDim,
-                                 const int NE,
-                                 const Array<double> &bo,
-                                 const Array<double> &bc,
-                                 const Array<double> &gc,
-                                 const Vector &pa_data,
-                                 const Vector &x,
-                                 Vector &y)
+void SmemPAHcurlL2Apply3D(const int D1D,
+                          const int Q1D,
+                          const int coeffDim,
+                          const int NE,
+                          const Array<double> &bo,
+                          const Array<double> &bc,
+                          const Array<double> &gc,
+                          const Vector &pa_data,
+                          const Vector &x,
+                          Vector &y)
 {
    MFEM_VERIFY(D1D <= MAX_D1D, "Error: D1D > MAX_D1D");
    MFEM_VERIFY(Q1D <= MAX_Q1D, "Error: Q1D > MAX_Q1D");
@@ -3585,18 +3585,18 @@ static void SmemPAHcurlL2Apply3D(const int D1D,
 // Apply to x corresponding to DOF's in H(curl) (trial), whose curl is
 // integrated against H(div) test functions corresponding to y.
 template<int MAX_D1D = HCURL_MAX_D1D, int MAX_Q1D = HCURL_MAX_Q1D>
-static void PAHcurlHdivApply3D(const int D1D,
-                               const int D1Dtest,
-                               const int Q1D,
-                               const int NE,
-                               const Array<double> &bo,
-                               const Array<double> &bc,
-                               const Array<double> &bot,
-                               const Array<double> &bct,
-                               const Array<double> &gc,
-                               const Vector &pa_data,
-                               const Vector &x,
-                               Vector &y)
+void PAHcurlHdivApply3D(const int D1D,
+                        const int D1Dtest,
+                        const int Q1D,
+                        const int NE,
+                        const Array<double> &bo,
+                        const Array<double> &bc,
+                        const Array<double> &bot,
+                        const Array<double> &bct,
+                        const Array<double> &gc,
+                        const Vector &pa_data,
+                        const Vector &x,
+                        Vector &y)
 {
    MFEM_VERIFY(D1D <= MAX_D1D, "Error: D1D > MAX_D1D");
    MFEM_VERIFY(Q1D <= MAX_Q1D, "Error: Q1D > MAX_Q1D");
@@ -4071,18 +4071,18 @@ void MixedVectorWeakCurlIntegrator::AssemblePA(const FiniteElementSpace
 // Apply to x corresponding to DOF's in H(curl) (trial), integrated against curl
 // of H(curl) test functions corresponding to y.
 template<int MAX_D1D = HCURL_MAX_D1D, int MAX_Q1D = HCURL_MAX_Q1D>
-static void PAHcurlL2Apply3DTranspose(const int D1D,
-                                      const int Q1D,
-                                      const int coeffDim,
-                                      const int NE,
-                                      const Array<double> &bo,
-                                      const Array<double> &bc,
-                                      const Array<double> &bot,
-                                      const Array<double> &bct,
-                                      const Array<double> &gct,
-                                      const Vector &pa_data,
-                                      const Vector &x,
-                                      Vector &y)
+void PAHcurlL2Apply3DTranspose(const int D1D,
+                               const int Q1D,
+                               const int coeffDim,
+                               const int NE,
+                               const Array<double> &bo,
+                               const Array<double> &bc,
+                               const Array<double> &bot,
+                               const Array<double> &bct,
+                               const Array<double> &gct,
+                               const Vector &pa_data,
+                               const Vector &x,
+                               Vector &y)
 {
    // See PAHcurlL2Apply3D for comments.
 
@@ -4413,16 +4413,16 @@ static void PAHcurlL2Apply3DTranspose(const int D1D,
 }
 
 template<int MAX_D1D = HCURL_MAX_D1D, int MAX_Q1D = HCURL_MAX_Q1D>
-static void SmemPAHcurlL2Apply3DTranspose(const int D1D,
-                                          const int Q1D,
-                                          const int coeffDim,
-                                          const int NE,
-                                          const Array<double> &bo,
-                                          const Array<double> &bc,
-                                          const Array<double> &gc,
-                                          const Vector &pa_data,
-                                          const Vector &x,
-                                          Vector &y)
+void SmemPAHcurlL2Apply3DTranspose(const int D1D,
+                                   const int Q1D,
+                                   const int coeffDim,
+                                   const int NE,
+                                   const Array<double> &bo,
+                                   const Array<double> &bc,
+                                   const Array<double> &gc,
+                                   const Vector &pa_data,
+                                   const Vector &x,
+                                   Vector &y)
 {
    MFEM_VERIFY(D1D <= MAX_D1D, "Error: D1D > MAX_D1D");
    MFEM_VERIFY(Q1D <= MAX_Q1D, "Error: Q1D > MAX_Q1D");
@@ -4675,13 +4675,13 @@ void MixedVectorWeakCurlIntegrator::AddMultPA(const Vector &x, Vector &y) const
 // Apply to x corresponding to DOFs in H^1 (domain) the (topological) gradient
 // to get a dof in H(curl) (range). You can think of the range as the "test" space
 // and the domain as the "trial" space, but there's no integration.
-static void PAHcurlApplyGradient2D(const int c_dofs1D,
-                                   const int o_dofs1D,
-                                   const int NE,
-                                   const Array<double> &B_,
-                                   const Array<double> &G_,
-                                   const Vector &x_,
-                                   Vector &y_)
+void PAHcurlApplyGradient2D(const int c_dofs1D,
+                            const int o_dofs1D,
+                            const int NE,
+                            const Array<double> &B_,
+                            const Array<double> &G_,
+                            const Vector &x_,
+                            Vector &y_)
 {
    auto B = Reshape(B_.Read(), c_dofs1D, c_dofs1D);
    auto G = Reshape(G_.Read(), o_dofs1D, c_dofs1D);
@@ -4753,12 +4753,12 @@ static void PAHcurlApplyGradient2D(const int c_dofs1D,
 }
 
 // Specialization of PAHcurlApplyGradient2D to the case where B is identity
-static void PAHcurlApplyGradient2DBId(const int c_dofs1D,
-                                      const int o_dofs1D,
-                                      const int NE,
-                                      const Array<double> &G_,
-                                      const Vector &x_,
-                                      Vector &y_)
+void PAHcurlApplyGradient2DBId(const int c_dofs1D,
+                               const int o_dofs1D,
+                               const int NE,
+                               const Array<double> &G_,
+                               const Vector &x_,
+                               Vector &y_)
 {
    auto G = Reshape(G_.Read(), o_dofs1D, c_dofs1D);
 
@@ -4822,7 +4822,7 @@ static void PAHcurlApplyGradient2DBId(const int c_dofs1D,
    });
 }
 
-static void PAHcurlApplyGradientTranspose2D(
+void PAHcurlApplyGradientTranspose2D(
    const int c_dofs1D, const int o_dofs1D, const int NE,
    const Array<double> &B_, const Array<double> &G_,
    const Vector &x_, Vector &y_)
@@ -4898,7 +4898,7 @@ static void PAHcurlApplyGradientTranspose2D(
 
 // Specialization of PAHcurlApplyGradientTranspose2D to the case where
 // B is identity
-static void PAHcurlApplyGradientTranspose2DBId(
+void PAHcurlApplyGradientTranspose2DBId(
    const int c_dofs1D, const int o_dofs1D, const int NE,
    const Array<double> &G_,
    const Vector &x_, Vector &y_)
@@ -4965,13 +4965,13 @@ static void PAHcurlApplyGradientTranspose2DBId(
    });
 }
 
-static void PAHcurlApplyGradient3D(const int c_dofs1D,
-                                   const int o_dofs1D,
-                                   const int NE,
-                                   const Array<double> &B_,
-                                   const Array<double> &G_,
-                                   const Vector &x_,
-                                   Vector &y_)
+void PAHcurlApplyGradient3D(const int c_dofs1D,
+                            const int o_dofs1D,
+                            const int NE,
+                            const Array<double> &B_,
+                            const Array<double> &G_,
+                            const Vector &x_,
+                            Vector &y_)
 {
    auto B = Reshape(B_.Read(), c_dofs1D, c_dofs1D);
    auto G = Reshape(G_.Read(), o_dofs1D, c_dofs1D);
@@ -5154,12 +5154,12 @@ static void PAHcurlApplyGradient3D(const int c_dofs1D,
 }
 
 // Specialization of PAHcurlApplyGradient3D to the case where
-static void PAHcurlApplyGradient3DBId(const int c_dofs1D,
-                                      const int o_dofs1D,
-                                      const int NE,
-                                      const Array<double> &G_,
-                                      const Vector &x_,
-                                      Vector &y_)
+void PAHcurlApplyGradient3DBId(const int c_dofs1D,
+                               const int o_dofs1D,
+                               const int NE,
+                               const Array<double> &G_,
+                               const Vector &x_,
+                               Vector &y_)
 {
    auto G = Reshape(G_.Read(), o_dofs1D, c_dofs1D);
 
@@ -5322,7 +5322,7 @@ static void PAHcurlApplyGradient3DBId(const int c_dofs1D,
    });
 }
 
-static void PAHcurlApplyGradientTranspose3D(
+void PAHcurlApplyGradientTranspose3D(
    const int c_dofs1D, const int o_dofs1D, const int NE,
    const Array<double> &B_, const Array<double> &G_,
    const Vector &x_, Vector &y_)
@@ -5507,7 +5507,7 @@ static void PAHcurlApplyGradientTranspose3D(
 }
 
 // Specialization of PAHcurlApplyGradientTranspose3D to the case where
-static void PAHcurlApplyGradientTranspose3DBId(
+void PAHcurlApplyGradientTranspose3DBId(
    const int c_dofs1D, const int o_dofs1D, const int NE,
    const Array<double> &G_,
    const Vector &x_, Vector &y_)
@@ -5789,14 +5789,14 @@ void GradientInterpolator::AddMultTransposePA(const Vector &x, Vector &y) const
    }
 }
 
-static void PAHcurlVecH1IdentityApply3D(const int c_dofs1D,
-                                        const int o_dofs1D,
-                                        const int NE,
-                                        const Array<double> &Bclosed,
-                                        const Array<double> &Bopen,
-                                        const Vector &pa_data,
-                                        const Vector &x_,
-                                        Vector &y_)
+void PAHcurlVecH1IdentityApply3D(const int c_dofs1D,
+                                 const int o_dofs1D,
+                                 const int NE,
+                                 const Array<double> &Bclosed,
+                                 const Array<double> &Bopen,
+                                 const Vector &pa_data,
+                                 const Vector &x_,
+                                 Vector &y_)
 {
    auto Bc = Reshape(Bclosed.Read(), c_dofs1D, c_dofs1D);
    auto Bo = Reshape(Bopen.Read(), o_dofs1D, c_dofs1D);
@@ -6002,14 +6002,14 @@ static void PAHcurlVecH1IdentityApply3D(const int c_dofs1D,
    });
 }
 
-static void PAHcurlVecH1IdentityApplyTranspose3D(const int c_dofs1D,
-                                                 const int o_dofs1D,
-                                                 const int NE,
-                                                 const Array<double> &Bclosed,
-                                                 const Array<double> &Bopen,
-                                                 const Vector &pa_data,
-                                                 const Vector &x_,
-                                                 Vector &y_)
+void PAHcurlVecH1IdentityApplyTranspose3D(const int c_dofs1D,
+                                          const int o_dofs1D,
+                                          const int NE,
+                                          const Array<double> &Bclosed,
+                                          const Array<double> &Bopen,
+                                          const Vector &pa_data,
+                                          const Vector &x_,
+                                          Vector &y_)
 {
    auto Bc = Reshape(Bclosed.Read(), c_dofs1D, c_dofs1D);
    auto Bo = Reshape(Bopen.Read(), o_dofs1D, c_dofs1D);
@@ -6228,14 +6228,14 @@ static void PAHcurlVecH1IdentityApplyTranspose3D(const int c_dofs1D,
    });
 }
 
-static void PAHcurlVecH1IdentityApply2D(const int c_dofs1D,
-                                        const int o_dofs1D,
-                                        const int NE,
-                                        const Array<double> &Bclosed,
-                                        const Array<double> &Bopen,
-                                        const Vector &pa_data,
-                                        const Vector &x_,
-                                        Vector &y_)
+void PAHcurlVecH1IdentityApply2D(const int c_dofs1D,
+                                 const int o_dofs1D,
+                                 const int NE,
+                                 const Array<double> &Bclosed,
+                                 const Array<double> &Bopen,
+                                 const Vector &pa_data,
+                                 const Vector &x_,
+                                 Vector &y_)
 {
    auto Bc = Reshape(Bclosed.Read(), c_dofs1D, c_dofs1D);
    auto Bo = Reshape(Bopen.Read(), o_dofs1D, c_dofs1D);
@@ -6327,14 +6327,14 @@ static void PAHcurlVecH1IdentityApply2D(const int c_dofs1D,
    });
 }
 
-static void PAHcurlVecH1IdentityApplyTranspose2D(const int c_dofs1D,
-                                                 const int o_dofs1D,
-                                                 const int NE,
-                                                 const Array<double> &Bclosed,
-                                                 const Array<double> &Bopen,
-                                                 const Vector &pa_data,
-                                                 const Vector &x_,
-                                                 Vector &y_)
+void PAHcurlVecH1IdentityApplyTranspose2D(const int c_dofs1D,
+                                          const int o_dofs1D,
+                                          const int NE,
+                                          const Array<double> &Bclosed,
+                                          const Array<double> &Bopen,
+                                          const Vector &pa_data,
+                                          const Vector &x_,
+                                          Vector &y_)
 {
    auto Bc = Reshape(Bclosed.Read(), c_dofs1D, c_dofs1D);
    auto Bo = Reshape(Bopen.Read(), o_dofs1D, c_dofs1D);

fem/bilininteg_hdiv.cpp

@@ -539,12 +539,12 @@ void PAHdivMassApply3D(const int D1D,
 
 // PA H(div) div-div assemble 2D kernel
 // NOTE: this is identical to PACurlCurlSetup3D
-static void PADivDivSetup2D(const int Q1D,
-                            const int NE,
-                            const Array<double> &w,
-                            const Vector &j,
-                            Vector &coeff_,
-                            Vector &op)
+void PADivDivSetup2D(const int Q1D,
+                     const int NE,
+                     const Array<double> &w,
+                     const Vector &j,
+                     Vector &coeff_,
+                     Vector &op)
 {
    const int NQ = Q1D*Q1D;
    auto W = w.Read();
@@ -565,12 +565,12 @@ static void PADivDivSetup2D(const int Q1D,
    });
 }
 
-static void PADivDivSetup3D(const int Q1D,
-                            const int NE,
-                            const Array<double> &w,
-                            const Vector &j,
-                            Vector &coeff_,
-                            Vector &op)
+void PADivDivSetup3D(const int Q1D,
+                     const int NE,
+                     const Array<double> &w,
+                     const Vector &j,
+                     Vector &coeff_,
+                     Vector &op)
 {
    const int NQ = Q1D*Q1D*Q1D;
    auto W = w.Read();
@@ -599,16 +599,16 @@ static void PADivDivSetup3D(const int Q1D,
    });
 }
 
-static void PADivDivApply2D(const int D1D,
-                            const int Q1D,
-                            const int NE,
-                            const Array<double> &Bo_,
-                            const Array<double> &Gc_,
-                            const Array<double> &Bot_,
-                            const Array<double> &Gct_,
-                            const Vector &op_,
-                            const Vector &x_,
-                            Vector &y_)
+void PADivDivApply2D(const int D1D,
+                     const int Q1D,
+                     const int NE,
+                     const Array<double> &Bo_,
+                     const Array<double> &Gc_,
+                     const Array<double> &Bot_,
+                     const Array<double> &Gct_,
+                     const Vector &op_,
+                     const Vector &x_,
+                     Vector &y_)
 {
    constexpr static int VDIM = 2;
    constexpr static int MAX_D1D = HDIV_MAX_D1D;
@@ -718,16 +718,16 @@ static void PADivDivApply2D(const int D1D,
    }); // end of element loop
 }
 
-static void PADivDivApply3D(const int D1D,
-                            const int Q1D,
-                            const int NE,
-                            const Array<double> &Bo_,
-                            const Array<double> &Gc_,
-                            const Array<double> &Bot_,
-                            const Array<double> &Gct_,
-                            const Vector &op_,
-                            const Vector &x_,
-                            Vector &y_)
+void PADivDivApply3D(const int D1D,
+                     const int Q1D,
+                     const int NE,
+                     const Array<double> &Bo_,
+                     const Array<double> &Gc_,
+                     const Array<double> &Bot_,
+                     const Array<double> &Gct_,
+                     const Vector &op_,
+                     const Vector &x_,
+                     Vector &y_)
 {
    MFEM_VERIFY(D1D <= HDIV_MAX_D1D, "Error: D1D > HDIV_MAX_D1D");
    MFEM_VERIFY(Q1D <= HDIV_MAX_Q1D, "Error: Q1D > HDIV_MAX_Q1D");
@@ -967,13 +967,13 @@ void DivDivIntegrator::AddMultPA(const Vector &x, Vector &y) const
    }
 }
 
-static void PADivDivAssembleDiagonal2D(const int D1D,
-                                       const int Q1D,
-                                       const int NE,
-                                       const Array<double> &Bo_,
-                                       const Array<double> &Gc_,
-                                       const Vector &op_,
-                                       Vector &diag_)
+void PADivDivAssembleDiagonal2D(const int D1D,
+                                const int Q1D,
+                                const int NE,
+                                const Array<double> &Bo_,
+                                const Array<double> &Gc_,
+                                const Vector &op_,
+                                Vector &diag_)
 {
    constexpr static int VDIM = 2;
    constexpr static int MAX_Q1D = HDIV_MAX_Q1D;
@@ -1023,13 +1023,13 @@ static void PADivDivAssembleDiagonal2D(const int D1D,
    });
 }
 
-static void PADivDivAssembleDiagonal3D(const int D1D,
-                                       const int Q1D,
-                                       const int NE,
-                                       const Array<double> &Bo_,
-                                       const Array<double> &Gc_,
-                                       const Vector &op_,
-                                       Vector &diag_)
+void PADivDivAssembleDiagonal3D(const int D1D,
+                                const int Q1D,
+                                const int NE,
+                                const Array<double> &Bo_,
+                                const Array<double> &Gc_,
+                                const Vector &op_,
+                                Vector &diag_)
 {
    MFEM_VERIFY(D1D <= HDIV_MAX_D1D, "Error: D1D > HDIV_MAX_D1D");
    MFEM_VERIFY(Q1D <= HDIV_MAX_Q1D, "Error: Q1D > HDIV_MAX_Q1D");
@@ -1104,11 +1104,11 @@ void DivDivIntegrator::AssembleDiagonalPA(Vector& diag)
 }
 
 // PA H(div)-L2 (div u, p) assemble 2D kernel
-static void PADivL2Setup2D(const int Q1D,
-                           const int NE,
-                           const Array<double> &w,
-                           Vector &coeff_,
-                           Vector &op)
+void PADivL2Setup2D(const int Q1D,
+                    const int NE,
+                    const Array<double> &w,
+                    Vector &coeff_,
+                    Vector &op)
 {
    const int NQ = Q1D*Q1D;
    auto W = w.Read();
@@ -1123,11 +1123,11 @@ static void PADivL2Setup2D(const int Q1D,
    });
 }
 
-static void PADivL2Setup3D(const int Q1D,
-                           const int NE,
-                           const Array<double> &w,
-                           Vector &coeff_,
-                           Vector &op)
+void PADivL2Setup3D(const int Q1D,
+                    const int NE,
+                    const Array<double> &w,
+                    Vector &coeff_,
+                    Vector &op)
 {
    const int NQ = Q1D*Q1D*Q1D;
    auto W = w.Read();
@@ -1225,16 +1225,16 @@ VectorFEDivergenceIntegrator::AssemblePA(const FiniteElementSpace &trial_fes,
 
 // Apply to x corresponding to DOF's in H(div) (trial), whose divergence is
 // integrated against L_2 test functions corresponding to y.
-static void PAHdivL2Apply3D(const int D1D,
-                            const int Q1D,
-                            const int L2D1D,
-                            const int NE,
-                            const Array<double> &Bo_,
-                            const Array<double> &Gc_,
-                            const Array<double> &L2Bot_,
-                            const Vector &op_,
-                            const Vector &x_,
-                            Vector &y_)
+void PAHdivL2Apply3D(const int D1D,
+                     const int Q1D,
+                     const int L2D1D,
+                     const int NE,
+                     const Array<double> &Bo_,
+                     const Array<double> &Gc_,
+                     const Array<double> &L2Bot_,
+                     const Vector &op_,
+                     const Vector &x_,
+                     Vector &y_)
 {
    MFEM_VERIFY(D1D <= HDIV_MAX_D1D, "Error: D1D > HDIV_MAX_D1D");
    MFEM_VERIFY(Q1D <= HDIV_MAX_Q1D, "Error: Q1D > HDIV_MAX_Q1D");
@@ -1388,16 +1388,16 @@ static void PAHdivL2Apply3D(const int D1D,
 
 // Apply to x corresponding to DOF's in H(div) (trial), whose divergence is
 // integrated against L_2 test functions corresponding to y.
-static void PAHdivL2Apply2D(const int D1D,
-                            const int Q1D,
-                            const int L2D1D,
-                            const int NE,
-                            const Array<double> &Bo_,
-                            const Array<double> &Gc_,
-                            const Array<double> &L2Bot_,
-                            const Vector &op_,
-                            const Vector &x_,
-                            Vector &y_)
+void PAHdivL2Apply2D(const int D1D,
+                     const int Q1D,
+                     const int L2D1D,
+                     const int NE,
+                     const Array<double> &Bo_,
+                     const Array<double> &Gc_,
+                     const Array<double> &L2Bot_,
+                     const Vector &op_,
+                     const Vector &x_,
+                     Vector &y_)
 {
    constexpr static int VDIM = 2;
    constexpr static int MAX_D1D = HDIV_MAX_D1D;
@@ -1494,16 +1494,16 @@ static void PAHdivL2Apply2D(const int D1D,
    }); // end of element loop
 }
 
-static void PAHdivL2ApplyTranspose3D(const int D1D,
-                                     const int Q1D,
-                                     const int L2D1D,
-                                     const int NE,
-                                     const Array<double> &L2Bo_,
-                                     const Array<double> &Gct_,
-                                     const Array<double> &Bot_,
-                                     const Vector &op_,
-                                     const Vector &x_,
-                                     Vector &y_)
+void PAHdivL2ApplyTranspose3D(const int D1D,
+                              const int Q1D,
+                              const int L2D1D,
+                              const int NE,
+                              const Array<double> &L2Bo_,
+                              const Array<double> &Gct_,
+                              const Array<double> &Bot_,
+                              const Vector &op_,
+                              const Vector &x_,
+                              Vector &y_)
 {
    MFEM_VERIFY(D1D <= HDIV_MAX_D1D, "Error: D1D > HDIV_MAX_D1D");
    MFEM_VERIFY(Q1D <= HDIV_MAX_Q1D, "Error: Q1D > HDIV_MAX_Q1D");
@@ -1656,16 +1656,16 @@ static void PAHdivL2ApplyTranspose3D(const int D1D,
    }); // end of element loop
 }
 
-static void PAHdivL2ApplyTranspose2D(const int D1D,
-                                     const int Q1D,
-                                     const int L2D1D,
-                                     const int NE,
-                                     const Array<double> &L2Bo_,
-                                     const Array<double> &Gct_,
-                                     const Array<double> &Bot_,
-                                     const Vector &op_,
-                                     const Vector &x_,
-                                     Vector &y_)
+void PAHdivL2ApplyTranspose2D(const int D1D,
+                              const int Q1D,
+                              const int L2D1D,
+                              const int NE,
+                              const Array<double> &L2Bo_,
+                              const Array<double> &Gct_,
+                              const Array<double> &Bot_,
+                              const Vector &op_,
+                              const Vector &x_,
+                              Vector &y_)
 {
    constexpr static int VDIM = 2;
    constexpr static int MAX_D1D = HDIV_MAX_D1D;
@@ -1791,16 +1791,16 @@ void VectorFEDivergenceIntegrator::AddMultTransposePA(const Vector &x,
    }
 }
 
-static void PAHdivL2AssembleDiagonal_ADAt_3D(const int D1D,
-                                             const int Q1D,
-                                             const int L2D1D,
-                                             const int NE,
-                                             const Array<double> &L2Bo_,
-                                             const Array<double> &Gct_,
-                                             const Array<double> &Bot_,
-                                             const Vector &op_,
-                                             const Vector &D_,
-                                             Vector &diag_)
+void PAHdivL2AssembleDiagonal_ADAt_3D(const int D1D,
+                                      const int Q1D,
+                                      const int L2D1D,
+                                      const int NE,
+                                      const Array<double> &L2Bo_,
+                                      const Array<double> &Gct_,
+                                      const Array<double> &Bot_,
+                                      const Vector &op_,
+                                      const Vector &D_,
+                                      Vector &diag_)
 {
    MFEM_VERIFY(D1D <= HDIV_MAX_D1D, "Error: D1D > HDIV_MAX_D1D");
    MFEM_VERIFY(Q1D <= HDIV_MAX_Q1D, "Error: Q1D > HDIV_MAX_Q1D");
@@ -1916,16 +1916,16 @@ static void PAHdivL2AssembleDiagonal_ADAt_3D(const int D1D,
    }); // end of element loop
 }
 
-static void PAHdivL2AssembleDiagonal_ADAt_2D(const int D1D,
-                                             const int Q1D,
-                                             const int L2D1D,
-                                             const int NE,
-                                             const Array<double> &L2Bo_,
-                                             const Array<double> &Gct_,
-                                             const Array<double> &Bot_,
-                                             const Vector &op_,
-                                             const Vector &D_,
-                                             Vector &diag_)
+void PAHdivL2AssembleDiagonal_ADAt_2D(const int D1D,
+                                      const int Q1D,
+                                      const int L2D1D,
+                                      const int NE,
+                                      const Array<double> &L2Bo_,
+                                      const Array<double> &Gct_,
+                                      const Array<double> &Bot_,
+                                      const Vector &op_,
+                                      const Vector &D_,
+                                      Vector &diag_)
 {
    constexpr static int VDIM = 2;
 

fem/bilininteg_mass_ea.cpp

@@ -17,13 +17,13 @@ namespace mfem
 {
 
 template<int T_D1D = 0, int T_Q1D = 0>
-static void EAMassAssemble1D(const int NE,
-                             const Array<double> &basis,
-                             const Vector &padata,
-                             Vector &eadata,
-                             const bool add,
-                             const int d1d = 0,
-                             const int q1d = 0)
+void EAMassAssemble1D(const int NE,
+                      const Array<double> &basis,
+                      const Vector &padata,
+                      Vector &eadata,
+                      const bool add,
+                      const int d1d = 0,
+                      const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -67,13 +67,13 @@ static void EAMassAssemble1D(const int NE,
 }
 
 template<int T_D1D = 0, int T_Q1D = 0>
-static void EAMassAssemble2D(const int NE,
-                             const Array<double> &basis,
-                             const Vector &padata,
-                             Vector &eadata,
-                             const bool add,
-                             const int d1d = 0,
-                             const int q1d = 0)
+void EAMassAssemble2D(const int NE,
+                      const Array<double> &basis,
+                      const Vector &padata,
+                      Vector &eadata,
+                      const bool add,
+                      const int d1d = 0,
+                      const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -139,13 +139,13 @@ static void EAMassAssemble2D(const int NE,
 }
 
 template<int T_D1D = 0, int T_Q1D = 0>
-static void EAMassAssemble3D(const int NE,
-                             const Array<double> &basis,
-                             const Vector &padata,
-                             Vector &eadata,
-                             const bool add,
-                             const int d1d = 0,
-                             const int q1d = 0)
+void EAMassAssemble3D(const int NE,
+                      const Array<double> &basis,
+                      const Vector &padata,
+                      Vector &eadata,
+                      const bool add,
+                      const int d1d = 0,
+                      const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;

fem/bilininteg_mass_pa.cpp

@@ -155,12 +155,12 @@ void MassIntegrator::AssemblePA(const FiniteElementSpace &fes)
 }
 
 template<int T_D1D = 0, int T_Q1D = 0>
-static void PAMassAssembleDiagonal2D(const int NE,
-                                     const Array<double> &b,
-                                     const Vector &d,
-                                     Vector &y,
-                                     const int d1d = 0,
-                                     const int q1d = 0)
+void PAMassAssembleDiagonal2D(const int NE,
+                              const Array<double> &b,
+                              const Vector &d,
+                              Vector &y,
+                              const int d1d = 0,
+                              const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -201,12 +201,12 @@ static void PAMassAssembleDiagonal2D(const int NE,
 }
 
 template<int T_D1D = 0, int T_Q1D = 0, int T_NBZ = 0>
-static void SmemPAMassAssembleDiagonal2D(const int NE,
-                                         const Array<double> &b_,
-                                         const Vector &d_,
-                                         Vector &y_,
-                                         const int d1d = 0,
-                                         const int q1d = 0)
+void SmemPAMassAssembleDiagonal2D(const int NE,
+                                  const Array<double> &b_,
+                                  const Vector &d_,
+                                  Vector &y_,
+                                  const int d1d = 0,
+                                  const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -267,12 +267,12 @@ static void SmemPAMassAssembleDiagonal2D(const int NE,
 }
 
 template<int T_D1D = 0, int T_Q1D = 0>
-static void PAMassAssembleDiagonal3D(const int NE,
-                                     const Array<double> &b,
-                                     const Vector &d,
-                                     Vector &y,
-                                     const int d1d = 0,
-                                     const int q1d = 0)
+void PAMassAssembleDiagonal3D(const int NE,
+                              const Array<double> &b,
+                              const Vector &d,
+                              Vector &y,
+                              const int d1d = 0,
+                              const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -336,12 +336,12 @@ static void PAMassAssembleDiagonal3D(const int NE,
 }
 
 template<int T_D1D = 0, int T_Q1D = 0>
-static void SmemPAMassAssembleDiagonal3D(const int NE,
-                                         const Array<double> &b_,
-                                         const Vector &d_,
-                                         Vector &y_,
-                                         const int d1d = 0,
-                                         const int q1d = 0)
+void SmemPAMassAssembleDiagonal3D(const int NE,
+                                  const Array<double> &b_,
+                                  const Vector &d_,
+                                  Vector &y_,
+                                  const int d1d = 0,
+                                  const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -569,14 +569,14 @@ static void OccaPAMassApply3D(const int D1D,
 #endif // MFEM_USE_OCCA
 
 template<int T_D1D = 0, int T_Q1D = 0>
-static void PAMassApply2D(const int NE,
-                          const Array<double> &b_,
-                          const Array<double> &bt_,
-                          const Vector &d_,
-                          const Vector &x_,
-                          Vector &y_,
-                          const int d1d = 0,
-                          const int q1d = 0)
+void PAMassApply2D(const int NE,
+                   const Array<double> &b_,
+                   const Array<double> &bt_,
+                   const Vector &d_,
+                   const Vector &x_,
+                   Vector &y_,
+                   const int d1d = 0,
+                   const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -661,14 +661,14 @@ static void PAMassApply2D(const int NE,
 }
 
 template<int T_D1D = 0, int T_Q1D = 0, int T_NBZ = 0>
-static void SmemPAMassApply2D(const int NE,
-                              const Array<double> &b_,
-                              const Array<double> &bt_,
-                              const Vector &d_,
-                              const Vector &x_,
-                              Vector &y_,
-                              const int d1d = 0,
-                              const int q1d = 0)
+void SmemPAMassApply2D(const int NE,
+                       const Array<double> &b_,
+                       const Array<double> &bt_,
+                       const Vector &d_,
+                       const Vector &x_,
+                       Vector &y_,
+                       const int d1d = 0,
+                       const int q1d = 0)
 {
    MFEM_CONTRACT_VAR(bt_);
    const int D1D = T_D1D ? T_D1D : d1d;
@@ -784,14 +784,14 @@ static void SmemPAMassApply2D(const int NE,
 }
 
 template<int T_D1D = 0, int T_Q1D = 0>
-static void PAMassApply3D(const int NE,
-                          const Array<double> &b_,
-                          const Array<double> &bt_,
-                          const Vector &d_,
-                          const Vector &x_,
-                          Vector &y_,
-                          const int d1d = 0,
-                          const int q1d = 0)
+void PAMassApply3D(const int NE,
+                   const Array<double> &b_,
+                   const Array<double> &bt_,
+                   const Vector &d_,
+                   const Vector &x_,
+                   Vector &y_,
+                   const int d1d = 0,
+                   const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -925,14 +925,14 @@ static void PAMassApply3D(const int NE,
 }
 
 template<int T_D1D = 0, int T_Q1D = 0>
-static void SmemPAMassApply3D(const int NE,
-                              const Array<double> &b_,
-                              const Array<double> &bt_,
-                              const Vector &d_,
-                              const Vector &x_,
-                              Vector &y_,
-                              const int d1d = 0,
-                              const int q1d = 0)
+void SmemPAMassApply3D(const int NE,
+                       const Array<double> &b_,
+                       const Array<double> &bt_,
+                       const Vector &d_,
+                       const Vector &x_,
+                       Vector &y_,
+                       const int d1d = 0,
+                       const int q1d = 0)
 {
    MFEM_CONTRACT_VAR(bt_);
    const int D1D = T_D1D ? T_D1D : d1d;

fem/bilininteg_vecdiffusion.cpp

@@ -22,12 +22,12 @@ namespace mfem
 // PA Vector Diffusion Integrator
 
 // PA Diffusion Assemble 2D kernel
-static void PAVectorDiffusionSetup2D(const int Q1D,
-                                     const int NE,
-                                     const Array<double> &w,
-                                     const Vector &j,
-                                     const Vector &c,
-                                     Vector &op)
+void PAVectorDiffusionSetup2D(const int Q1D,
+                              const int NE,
+                              const Array<double> &w,
+                              const Vector &j,
+                              const Vector &c,
+                              Vector &op)
 {
    const int NQ = Q1D*Q1D;
    auto W = w.Read();
@@ -59,12 +59,12 @@ static void PAVectorDiffusionSetup2D(const int Q1D,
 }
 
 // PA Diffusion Assemble 3D kernel
-static void PAVectorDiffusionSetup3D(const int Q1D,
-                                     const int NE,
-                                     const Array<double> &w,
-                                     const Vector &j,
-                                     const Vector &c,
-                                     Vector &op)
+void PAVectorDiffusionSetup3D(const int Q1D,
+                              const int NE,
+                              const Array<double> &w,
+                              const Vector &j,
+                              const Vector &c,
+                              Vector &op)
 {
    const int NQ = Q1D*Q1D*Q1D;
    auto W = w.Read();
@@ -251,7 +251,7 @@ void VectorDiffusionIntegrator::AssemblePA(const FiniteElementSpace &fes)
 }
 
 // PA Diffusion Apply 2D kernel
-template<int T_D1D = 0, int T_Q1D = 0, int T_VDIM = 0> static
+template<int T_D1D = 0, int T_Q1D = 0, int T_VDIM = 0>
 void PAVectorDiffusionApply2D(const int NE,
                               const Array<double> &b,
                               const Array<double> &g,
@@ -374,7 +374,7 @@ void PAVectorDiffusionApply2D(const int NE,
 
 // PA Diffusion Apply 3D kernel
 template<const int T_D1D = 0,
-         const int T_Q1D = 0> static
+         const int T_Q1D = 0>
 void PAVectorDiffusionApply3D(const int NE,
                               const Array<double> &b,
                               const Array<double> &g,
@@ -606,13 +606,13 @@ void VectorDiffusionIntegrator::AddMultPA(const Vector &x, Vector &y) const
 }
 
 template<int T_D1D = 0, int T_Q1D = 0>
-static void PAVectorDiffusionDiagonal2D(const int NE,
-                                        const Array<double> &b,
-                                        const Array<double> &g,
-                                        const Vector &d,
-                                        Vector &y,
-                                        const int d1d = 0,
-                                        const int q1d = 0)
+void PAVectorDiffusionDiagonal2D(const int NE,
+                                 const Array<double> &b,
+                                 const Array<double> &g,
+                                 const Vector &d,
+                                 Vector &y,
+                                 const int d1d = 0,
+                                 const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -673,13 +673,13 @@ static void PAVectorDiffusionDiagonal2D(const int NE,
 }
 
 template<int T_D1D = 0, int T_Q1D = 0>
-static void PAVectorDiffusionDiagonal3D(const int NE,
-                                        const Array<double> &b,
-                                        const Array<double> &g,
-                                        const Vector &d,
-                                        Vector &y,
-                                        const int d1d = 0,
-                                        const int q1d = 0)
+void PAVectorDiffusionDiagonal3D(const int NE,
+                                 const Array<double> &b,
+                                 const Array<double> &g,
+                                 const Vector &d,
+                                 Vector &y,
+                                 const int d1d = 0,
+                                 const int q1d = 0)
 {
    constexpr int DIM = 3;
    const int D1D = T_D1D ? T_D1D : d1d;

fem/bilininteg_vecmass.cpp

@@ -104,14 +104,14 @@ void VectorMassIntegrator::AssemblePA(const FiniteElementSpace &fes)
 
 template<const int T_D1D = 0,
          const int T_Q1D = 0>
-static void PAVectorMassApply2D(const int NE,
-                                const Array<double> &B_,
-                                const Array<double> &Bt_,
-                                const Vector &op_,
-                                const Vector &x_,
-                                Vector &y_,
-                                const int d1d = 0,
-                                const int q1d = 0)
+void PAVectorMassApply2D(const int NE,
+                         const Array<double> &B_,
+                         const Array<double> &Bt_,
+                         const Vector &op_,
+                         const Vector &x_,
+                         Vector &y_,
+                         const int d1d = 0,
+                         const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -201,14 +201,14 @@ static void PAVectorMassApply2D(const int NE,
 
 template<const int T_D1D = 0,
          const int T_Q1D = 0>
-static void PAVectorMassApply3D(const int NE,
-                                const Array<double> &B_,
-                                const Array<double> &Bt_,
-                                const Vector &op_,
-                                const Vector &x_,
-                                Vector &y_,
-                                const int d1d = 0,
-                                const int q1d = 0)
+void PAVectorMassApply3D(const int NE,
+                         const Array<double> &B_,
+                         const Array<double> &Bt_,
+                         const Vector &op_,
+                         const Vector &x_,
+                         Vector &y_,
+                         const int d1d = 0,
+                         const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -379,13 +379,13 @@ void VectorMassIntegrator::AddMultPA(const Vector &x, Vector &y) const
 }
 
 template<const int T_D1D = 0, const int T_Q1D = 0>
-static void PAVectorMassAssembleDiagonal2D(const int NE,
-                                           const Array<double> &B_,
-                                           const Array<double> &Bt_,
-                                           const Vector &op_,
-                                           Vector &diag_,
-                                           const int d1d = 0,
-                                           const int q1d = 0)
+void PAVectorMassAssembleDiagonal2D(const int NE,
+                                    const Array<double> &B_,
+                                    const Array<double> &Bt_,
+                                    const Vector &op_,
+                                    Vector &diag_,
+                                    const int d1d = 0,
+                                    const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -431,13 +431,13 @@ static void PAVectorMassAssembleDiagonal2D(const int NE,
 }
 
 template<const int T_D1D = 0, const int T_Q1D = 0>
-static void PAVectorMassAssembleDiagonal3D(const int NE,
-                                           const Array<double> &B_,
-                                           const Array<double> &Bt_,
-                                           const Vector &op_,
-                                           Vector &diag_,
-                                           const int d1d = 0,
-                                           const int q1d = 0)
+void PAVectorMassAssembleDiagonal3D(const int NE,
+                                    const Array<double> &B_,
+                                    const Array<double> &Bt_,
+                                    const Vector &op_,
+                                    Vector &diag_,
+                                    const int d1d = 0,
+                                    const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;

fem/eltrans.cpp

@@ -380,6 +380,7 @@ void IsoparametricTransformation::SetIdentityTransformation(
       case Geometry::TETRAHEDRON : FElem = &TetrahedronFE; break;
       case Geometry::CUBE :        FElem = &HexahedronFE; break;
       case Geometry::PRISM :       FElem = &WedgeFE; break;
+      case Geometry::PYRAMID :     FElem = &PyramidFE; break;
       default:
          MFEM_ABORT("unknown Geometry::Type!");
    }

fem/fe.hpp

@@ -1313,6 +1313,64 @@ public:
                            DenseMatrix &dshape) const;
 };
 
+/// A linear element defined on a triangular prism
+class LinearWedgeFiniteElement : public NodalFiniteElement
+{
+public:
+   /// Construct the LinearWedgeFiniteElement
+   LinearWedgeFiniteElement();
+
+   /** @brief virtual function which evaluates the values of all
+       shape functions at a given point ip and stores
+       them in the vector shape of dimension Dof (4) */
+   virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
+
+   /** @brief virtual function which evaluates the values of all
+       partial derivatives of all shape functions at a given
+       point ip and stores them in the matrix dshape (Dof x Dim) (4 x 3)
+       so that each row contains the derivatives of one shape function */
+   virtual void CalcDShape(const IntegrationPoint &ip,
+                           DenseMatrix &dshape) const;
+
+   virtual void ProjectDelta(int vertex, Vector &dofs) const
+   { dofs = 0.0; dofs(vertex) = 1.0; }
+
+   /** @brief Get the dofs associated with the given @a face.
+       @a *dofs is set to an internal array of the local dofc on the
+       face, while *ndofs is set to the number of dofs on that face.
+   */
+   virtual void GetFaceDofs(int face, int **dofs, int *ndofs) const;
+};
+
+/// A linear element defined on a square pyramid
+class LinearPyramidFiniteElement : public NodalFiniteElement
+{
+public:
+   /// Construct the LinearPyramidFiniteElement
+   LinearPyramidFiniteElement();
+
+   /** @brief virtual function which evaluates the values of all
+       shape functions at a given point ip and stores
+       them in the vector shape of dimension Dof (4) */
+   virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
+
+   /** @brief virtual function which evaluates the values of all
+       partial derivatives of all shape functions at a given
+       point ip and stores them in the matrix dshape (Dof x Dim) (4 x 3)
+       so that each row contains the derivatives of one shape function */
+   virtual void CalcDShape(const IntegrationPoint &ip,
+                           DenseMatrix &dshape) const;
+
+   virtual void ProjectDelta(int vertex, Vector &dofs) const
+   { dofs = 0.0; dofs(vertex) = 1.0; }
+
+   /** @brief Get the dofs associated with the given @a face.
+       @a *dofs is set to an internal array of the local dofc on the
+       face, while *ndofs is set to the number of dofs on that face.
+   */
+   virtual void GetFaceDofs(int face, int **dofs, int *ndofs) const;
+};
+
 /// A 2D constant element on a triangle
 class P0TriangleFiniteElement : public NodalFiniteElement
 {
@@ -1690,6 +1748,32 @@ public:
    { dofs(0) = 1.0; }
 };
 
+/// A 3D constant element on a wedge
+class P0WdgFiniteElement : public NodalFiniteElement
+{
+public:
+   /// Construct the P0WdgFiniteElement
+   P0WdgFiniteElement ();
+   virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
+   virtual void CalcDShape(const IntegrationPoint &ip,
+                           DenseMatrix &dshape) const;
+   virtual void ProjectDelta(int vertex, Vector &dofs) const
+   { dofs(0) = 1.0; }
+};
+
+/// A 3D constant element on a pyramid
+class P0PyrFiniteElement : public NodalFiniteElement
+{
+public:
+   /// Construct the P0PyrFiniteElement
+   P0PyrFiniteElement ();
+   virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const;
+   virtual void CalcDShape(const IntegrationPoint &ip,
+                           DenseMatrix &dshape) const;
+   virtual void ProjectDelta(int vertex, Vector &dofs) const
+   { dofs(0) = 1.0; }
+};
+
 /** @brief Tensor products of 1D Lagrange1DFiniteElement
     (only degree 2 is functional) */
 class LagrangeHexFiniteElement : public NodalFiniteElement
@@ -1828,6 +1912,10 @@ public:
    using FiniteElement::Project;
    virtual void Project (VectorCoefficient &vc,
                          ElementTransformation &Trans, Vector &dofs) const;
+
+   virtual void ProjectGrad(const FiniteElement &fe,
+                            ElementTransformation &Trans,
+                            DenseMatrix &grad) const;
 };
 
 
@@ -1852,6 +1940,66 @@ public:
    using FiniteElement::Project;
    virtual void Project (VectorCoefficient &vc,
                          ElementTransformation &Trans, Vector &dofs) const;
+
+   virtual void ProjectGrad(const FiniteElement &fe,
+                            ElementTransformation &Trans,
+                            DenseMatrix &grad) const;
+};
+
+
+/// A 3D 1st order Nedelec element on a wedge
+class Nedelec1WdgFiniteElement : public VectorFiniteElement
+{
+private:
+   static const double tk[9][3];
+
+public:
+   /// Construct the Nedelec1WdgFiniteElement
+   Nedelec1WdgFiniteElement();
+   virtual void CalcVShape(const IntegrationPoint &ip,
+                           DenseMatrix &shape) const;
+   virtual void CalcVShape(ElementTransformation &Trans,
+                           DenseMatrix &shape) const
+   { CalcVShape_ND(Trans, shape); }
+   virtual void CalcCurlShape(const IntegrationPoint &ip,
+                              DenseMatrix &curl_shape) const;
+   virtual void GetLocalInterpolation (ElementTransformation &Trans,
+                                       DenseMatrix &I) const;
+   using FiniteElement::Project;
+   virtual void Project (VectorCoefficient &vc,
+                         ElementTransformation &Trans, Vector &dofs) const;
+
+   virtual void ProjectGrad(const FiniteElement &fe,
+                            ElementTransformation &Trans,
+                            DenseMatrix &grad) const;
+};
+
+
+/// A 3D 1st order Nedelec element on a pyramid
+class Nedelec1PyrFiniteElement : public VectorFiniteElement
+{
+private:
+   static const double tk[8][3];
+
+public:
+   /// Construct the Nedelec1PyrFiniteElement
+   Nedelec1PyrFiniteElement();
+   virtual void CalcVShape(const IntegrationPoint &ip,
+                           DenseMatrix &shape) const;
+   virtual void CalcVShape(ElementTransformation &Trans,
+                           DenseMatrix &shape) const
+   { CalcVShape_ND(Trans, shape); }
+   virtual void CalcCurlShape(const IntegrationPoint &ip,
+                              DenseMatrix &curl_shape) const;
+   virtual void GetLocalInterpolation (ElementTransformation &Trans,
+                                       DenseMatrix &I) const;
+   using FiniteElement::Project;
+   virtual void Project (VectorCoefficient &vc,
+                         ElementTransformation &Trans, Vector &dofs) const;
+
+   virtual void ProjectGrad(const FiniteElement &fe,
+                            ElementTransformation &Trans,
+                            DenseMatrix &grad) const;
 };
 
 
@@ -1945,6 +2093,77 @@ public:
 };
 
 
+/// A 3D 0th order Raviert-Thomas element on a wedge
+class RT0WdgFiniteElement : public VectorFiniteElement
+{
+private:
+   static const double nk[5][3];
+
+public:
+   /// Construct the RT0WdgFiniteElement
+   RT0WdgFiniteElement();
+
+   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;
+
+   using FiniteElement::Project;
+
+   virtual void Project (VectorCoefficient &vc,
+                         ElementTransformation &Trans, Vector &dofs) const;
+
+   virtual void ProjectCurl(const FiniteElement &fe,
+                            ElementTransformation &Trans,
+                            DenseMatrix &curl) const;
+};
+
+
+/// A 3D 0th order Raviert-Thomas element on a pyramid
+class RT0PyrFiniteElement : public VectorFiniteElement
+{
+private:
+   static const double nk[5][3];
+
+   // If true match RT0TetFiniteElement rather than RT_TetrahedronElement(0)
+   bool rt0;
+
+public:
+   /// Construct the RT0PyrFiniteElement
+   RT0PyrFiniteElement(bool rt0tets = true);
+
+   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;
+
+   using FiniteElement::Project;
+
+   virtual void Project (VectorCoefficient &vc,
+                         ElementTransformation &Trans, Vector &dofs) const;
+
+   virtual void ProjectCurl(const FiniteElement &fe,
+                            ElementTransformation &Trans,
+                            DenseMatrix &curl) const;
+};
+
+
 class RotTriLinearHexFiniteElement : public NodalFiniteElement
 {
 public:

fem/fe_coll.cpp

@@ -33,6 +33,9 @@ int FiniteElementCollection::HasFaceDofs(Geometry::Type geom, int p) const
       case Geometry::PRISM:
          return max(GetNumDof(Geometry::TRIANGLE, p),
                     GetNumDof(Geometry::SQUARE, p));
+      case Geometry::PYRAMID:
+         return max(GetNumDof(Geometry::TRIANGLE, p),
+                    GetNumDof(Geometry::SQUARE, p));
       default:
          MFEM_ABORT("unknown geometry type");
    }
@@ -574,6 +577,7 @@ LinearFECollection::FiniteElementForGeometry(Geometry::Type GeomType) const
       case Geometry::TETRAHEDRON: return &TetrahedronFE;
       case Geometry::CUBE:        return &ParallelepipedFE;
       case Geometry::PRISM:       return &WedgeFE;
+      case Geometry::PYRAMID:     return &PyramidFE;
       default:
          mfem_error ("LinearFECollection: unknown geometry type.");
    }
@@ -591,6 +595,7 @@ int LinearFECollection::DofForGeometry(Geometry::Type GeomType) const
       case Geometry::TETRAHEDRON: return 0;
       case Geometry::CUBE:        return 0;
       case Geometry::PRISM:       return 0;
+      case Geometry::PYRAMID:     return 0;
       default:
          mfem_error ("LinearFECollection: unknown geometry type.");
    }
@@ -1240,6 +1245,7 @@ Const3DFECollection::FiniteElementForGeometry(Geometry::Type GeomType) const
       case Geometry::TETRAHEDRON: return &TetrahedronFE;
       case Geometry::CUBE:        return &ParallelepipedFE;
       case Geometry::PRISM:       return &WedgeFE;
+      case Geometry::PYRAMID:     return &PyramidFE;
       default:
          mfem_error ("Const3DFECollection: unknown geometry type.");
    }
@@ -1257,6 +1263,7 @@ int Const3DFECollection::DofForGeometry(Geometry::Type GeomType) const
       case Geometry::TETRAHEDRON: return 1;
       case Geometry::CUBE:        return 1;
       case Geometry::PRISM:       return 1;
+      case Geometry::PYRAMID:     return 1;
       default:
          mfem_error ("Const3DFECollection: unknown geometry type.");
    }
@@ -1277,6 +1284,8 @@ LinearDiscont3DFECollection::FiniteElementForGeometry(
    switch (GeomType)
    {
       case Geometry::TETRAHEDRON: return &TetrahedronFE;
+      case Geometry::PYRAMID:     return &PyramidFE;
+      case Geometry::PRISM:       return &WedgeFE;
       case Geometry::CUBE:        return &ParallelepipedFE;
       default:
          mfem_error ("LinearDiscont3DFECollection: unknown geometry type.");
@@ -1293,6 +1302,8 @@ int LinearDiscont3DFECollection::DofForGeometry(Geometry::Type GeomType) const
       case Geometry::TRIANGLE:    return 0;
       case Geometry::SQUARE:      return 0;
       case Geometry::TETRAHEDRON: return 4;
+      case Geometry::PYRAMID:     return 5;
+      case Geometry::PRISM:       return 6;
       case Geometry::CUBE:        return 8;
       default:
          mfem_error ("LinearDiscont3DFECollection: unknown geometry type.");
@@ -1394,6 +1405,8 @@ ND1_3DFECollection::FiniteElementForGeometry(Geometry::Type GeomType) const
    {
       case Geometry::CUBE:        return &HexahedronFE;
       case Geometry::TETRAHEDRON: return &TetrahedronFE;
+      case Geometry::PRISM:       return &WedgeFE;
+      case Geometry::PYRAMID:     return &PyramidFE;
       default:
          mfem_error ("ND1_3DFECollection: unknown geometry type.");
    }
@@ -1410,6 +1423,8 @@ int ND1_3DFECollection::DofForGeometry(Geometry::Type GeomType) const
       case Geometry::SQUARE:      return 0;
       case Geometry::TETRAHEDRON: return 0;
       case Geometry::CUBE:        return 0;
+      case Geometry::PRISM:       return 0;
+      case Geometry::PYRAMID:     return 0;
       default:
          mfem_error ("ND1_3DFECollection: unknown geometry type.");
    }
@@ -1439,6 +1454,8 @@ RT0_3DFECollection::FiniteElementForGeometry(Geometry::Type GeomType) const
       case Geometry::SQUARE:      return &QuadrilateralFE;
       case Geometry::CUBE:        return &HexahedronFE;
       case Geometry::TETRAHEDRON: return &TetrahedronFE;
+      case Geometry::PRISM:       return &WedgeFE;
+      case Geometry::PYRAMID:     return &PyramidFE;
       default:
          mfem_error ("RT0_3DFECollection: unknown geometry type.");
    }
@@ -1455,6 +1472,8 @@ int RT0_3DFECollection::DofForGeometry(Geometry::Type GeomType) const
       case Geometry::SQUARE:      return 1;
       case Geometry::TETRAHEDRON: return 0;
       case Geometry::CUBE:        return 0;
+      case Geometry::PRISM:       return 0;
+      case Geometry::PYRAMID:     return 0;
       default:
          mfem_error ("RT0_3DFECollection: unknown geometry type.");
    }
@@ -1730,6 +1749,7 @@ H1_FECollection::H1_FECollection(const int p, const int dim, const int btype)
          H1_dof[Geometry::TETRAHEDRON] = (TriDof*pm3)/3;
          H1_dof[Geometry::CUBE] = QuadDof*pm1;
          H1_dof[Geometry::PRISM] = TriDof*pm1;
+         H1_dof[Geometry::PYRAMID] = 0;
          if (b_type == BasisType::Positive)
          {
             H1_Elements[Geometry::TETRAHEDRON] = new H1Pos_TetrahedronElement(p);
@@ -1743,6 +1763,7 @@ H1_FECollection::H1_FECollection(const int p, const int dim, const int btype)
             H1_Elements[Geometry::CUBE] = new H1_HexahedronElement(p, btype);
             H1_Elements[Geometry::PRISM] = new H1_WedgeElement(p, btype);
          }
+         H1_Elements[Geometry::PYRAMID] = new LinearPyramidFiniteElement;
 
          const int &TetDof = H1_dof[Geometry::TETRAHEDRON];
          TetDofOrd[0] = new int[24*TetDof];
@@ -1837,6 +1858,21 @@ H1_FECollection::H1_FECollection(const int p, const int dim, const int btype)
    }
 }
 
+const FiniteElement *
+H1_FECollection::FiniteElementForGeometry(Geometry::Type GeomType) const
+{
+   if (GeomType != Geometry::PYRAMID || this->GetOrder() == 1)
+   {
+      return H1_Elements[GeomType];
+   }
+   else
+   {
+      MFEM_ABORT("H1 Pyramid basis functions are not yet supported "
+                 "for order > 1.");
+      return NULL;
+   }
+}
+
 const int *H1_FECollection::DofOrderForOrientation(Geometry::Type GeomType,
                                                    int Or) const
 {
@@ -2076,9 +2112,12 @@ L2_FECollection::L2_FECollection(const int p, const int dim, const int btype,
          L2_Elements[Geometry::CUBE] = new L2_HexahedronElement(p, btype);
          L2_Elements[Geometry::PRISM] = new L2_WedgeElement(p, btype);
       }
+      L2_Elements[Geometry::PYRAMID] = new P0PyrFiniteElement;
+
       L2_Elements[Geometry::TETRAHEDRON]->SetMapType(map_type);
       L2_Elements[Geometry::CUBE]->SetMapType(map_type);
       L2_Elements[Geometry::PRISM]->SetMapType(map_type);
+      L2_Elements[Geometry::PYRAMID]->SetMapType(map_type);
       // Trace element use the default Gauss-Legendre nodal points for positive basis
       if (b_type == BasisType::Positive)
       {
@@ -2199,6 +2238,21 @@ L2_FECollection::L2_FECollection(const int p, const int dim, const int btype,
    }
 }
 
+const FiniteElement *
+L2_FECollection::FiniteElementForGeometry(Geometry::Type GeomType) const
+{
+   if (GeomType != Geometry::PYRAMID || this->GetOrder() == 0)
+   {
+      return L2_Elements[GeomType];
+   }
+   else
+   {
+      MFEM_ABORT("L2 Pyramid basis functions are not yet supported "
+                 "for order > 0.");
+      return NULL;
+   }
+}
+
 const int *L2_FECollection::DofOrderForOrientation(Geometry::Type GeomType,
                                                    int Or) const
 {
@@ -2290,6 +2344,12 @@ RT_FECollection::RT_FECollection(const int order, const int dim,
 
       RT_Elements[Geometry::CUBE] = new RT_HexahedronElement(p, cb_type, ob_type);
       RT_dof[Geometry::CUBE] = 3*p*pp1*pp1;
+
+      RT_Elements[Geometry::PRISM] = new RT0WdgFiniteElement;
+      RT_dof[Geometry::PRISM] = 0;
+
+      RT_Elements[Geometry::PYRAMID] = new RT0PyrFiniteElement(false);
+      RT_dof[Geometry::PYRAMID] = 0;
    }
    else
    {
@@ -2433,6 +2493,22 @@ void RT_FECollection::InitFaces(const int p, const int dim,
    }
 }
 
+const FiniteElement *
+RT_FECollection::FiniteElementForGeometry(Geometry::Type GeomType) const
+{
+   if ((GeomType != Geometry::PRISM && GeomType != Geometry::PYRAMID) ||
+       this->GetOrder() == 1)
+   {
+      return RT_Elements[GeomType];
+   }
+   else
+   {
+      MFEM_ABORT("RT Wedge and Pyramid basis functions are not yet supported "
+                 "for order > 0.");
+      return NULL;
+   }
+}
+
 const int *RT_FECollection::DofOrderForOrientation(Geometry::Type GeomType,
                                                    int Or) const
 {
@@ -2693,6 +2769,28 @@ ND_FECollection::ND_FECollection(const int p, const int dim,
       // TODO: cb_type and ob_type for tets
       ND_Elements[Geometry::TETRAHEDRON] = new ND_TetrahedronElement(p);
       ND_dof[Geometry::TETRAHEDRON] = p*pm1*pm2/2;
+
+      ND_Elements[Geometry::PRISM] = new Nedelec1WdgFiniteElement;
+      ND_dof[Geometry::PRISM] = 0;
+
+      ND_Elements[Geometry::PYRAMID] = new Nedelec1PyrFiniteElement;
+      ND_dof[Geometry::PYRAMID] = 0;
+   }
+}
+
+const FiniteElement *
+ND_FECollection::FiniteElementForGeometry(Geometry::Type GeomType) const
+{
+   if ((GeomType != Geometry::PRISM && GeomType != Geometry::PYRAMID) ||
+       this->GetOrder() == 1)
+   {
+      return ND_Elements[GeomType];
+   }
+   else
+   {
+      MFEM_ABORT("ND Wedge and Pyramid basis functions are not yet supported "
+                 "for order > 1.");
+      return NULL;
    }
 }
 

fem/fe_coll.hpp

@@ -228,8 +228,7 @@ public:
                             const int btype = BasisType::GaussLobatto);
 
    virtual const FiniteElement *FiniteElementForGeometry(
-      Geometry::Type GeomType) const
-   { return H1_Elements[GeomType]; }
+      Geometry::Type GeomType) const;
    virtual int DofForGeometry(Geometry::Type GeomType) const
    { return H1_dof[GeomType]; }
    virtual const int *DofOrderForOrientation(Geometry::Type GeomType,
@@ -302,10 +301,7 @@ public:
                    const int map_type = FiniteElement::VALUE);
 
    virtual const FiniteElement *FiniteElementForGeometry(
-      Geometry::Type GeomType) const
-   {
-      return L2_Elements[GeomType];
-   }
+      Geometry::Type GeomType) const;
    virtual int DofForGeometry(Geometry::Type GeomType) const
    {
       if (L2_Elements[GeomType])
@@ -371,8 +367,7 @@ public:
                    const int ob_type = BasisType::GaussLegendre);
 
    virtual const FiniteElement *FiniteElementForGeometry(
-      Geometry::Type GeomType) const
-   { return RT_Elements[GeomType]; }
+      Geometry::Type GeomType) const;
    virtual int DofForGeometry(Geometry::Type GeomType) const
    { return RT_dof[GeomType]; }
    virtual const int *DofOrderForOrientation(Geometry::Type GeomType,
@@ -430,8 +425,7 @@ public:
                    const int ob_type = BasisType::GaussLegendre);
 
    virtual const FiniteElement *
-   FiniteElementForGeometry(Geometry::Type GeomType) const
-   { return ND_Elements[GeomType]; }
+   FiniteElementForGeometry(Geometry::Type GeomType) const;
 
    virtual int DofForGeometry(Geometry::Type GeomType) const
    { return ND_dof[GeomType]; }
@@ -529,9 +523,10 @@ private:
    const BiLinear2DFiniteElement QuadrilateralFE;
    const Linear3DFiniteElement TetrahedronFE;
    const TriLinear3DFiniteElement ParallelepipedFE;
-   const H1_WedgeElement WedgeFE;
+   const LinearWedgeFiniteElement WedgeFE;
+   const LinearPyramidFiniteElement PyramidFE;
 public:
-   LinearFECollection() : FiniteElementCollection(1), WedgeFE(1) { }
+   LinearFECollection() : FiniteElementCollection(1) { }
 
    virtual const FiniteElement *
    FiniteElementForGeometry(Geometry::Type GeomType) const;
@@ -936,10 +931,11 @@ class Const3DFECollection : public FiniteElementCollection
 private:
    const P0TetFiniteElement TetrahedronFE;
    const P0HexFiniteElement ParallelepipedFE;
-   const L2_WedgeElement WedgeFE;
+   const P0WdgFiniteElement WedgeFE;
+   const P0PyrFiniteElement PyramidFE;
 
 public:
-   Const3DFECollection() : FiniteElementCollection(0), WedgeFE(0) { }
+   Const3DFECollection() : FiniteElementCollection(0) { }
 
    virtual const FiniteElement *
    FiniteElementForGeometry(Geometry::Type GeomType) const;
@@ -960,6 +956,8 @@ class LinearDiscont3DFECollection : public FiniteElementCollection
 {
 private:
    const Linear3DFiniteElement TetrahedronFE;
+   const LinearPyramidFiniteElement PyramidFE;
+   const LinearWedgeFiniteElement WedgeFE;
    const TriLinear3DFiniteElement ParallelepipedFE;
 
 public:
@@ -1036,6 +1034,8 @@ class ND1_3DFECollection : public FiniteElementCollection
 private:
    const Nedelec1HexFiniteElement HexahedronFE;
    const Nedelec1TetFiniteElement TetrahedronFE;
+   const Nedelec1WdgFiniteElement WedgeFE;
+   const Nedelec1PyrFiniteElement PyramidFE;
 
 public:
    ND1_3DFECollection() : FiniteElementCollection(1) { }
@@ -1061,6 +1061,8 @@ private:
    const P0QuadFiniteElement QuadrilateralFE;
    const RT0HexFiniteElement HexahedronFE;
    const RT0TetFiniteElement TetrahedronFE;
+   const RT0WdgFiniteElement WedgeFE;
+   const RT0PyrFiniteElement PyramidFE;
 public:
    RT0_3DFECollection() : FiniteElementCollection(1) { }
 

fem/geom.cpp

@@ -11,15 +11,19 @@
 
 #include "fem.hpp"
 #include "../mesh/wedge.hpp"
+#include "../mesh/pyramid.hpp"
 
 namespace mfem
 {
 
 const char *Geometry::Name[NumGeom] =
-{ "Point", "Segment", "Triangle", "Square", "Tetrahedron", "Cube", "Prism" };
+{
+   "Point", "Segment", "Triangle", "Square", "Tetrahedron", "Cube", "Prism",
+   "Pyramid"
+};
 
 const double Geometry::Volume[NumGeom] =
-{ 1.0, 1.0, 0.5, 1.0, 1./6, 1.0, 0.5 };
+{ 1.0, 1.0, 0.5, 1.0, 1./6, 1.0, 0.5, 1./3 };
 
 Geometry::Geometry()
 {
@@ -139,6 +143,28 @@ Geometry::Geometry()
    GeomVert[6]->IntPoint(5).y = 1.0;
    GeomVert[6]->IntPoint(5).z = 1.0;
 
+   // Vertices for Geometry::PYRAMID
+   GeomVert[7] = new IntegrationRule(5);
+   GeomVert[7]->IntPoint(0).x = 0.0;
+   GeomVert[7]->IntPoint(0).y = 0.0;
+   GeomVert[7]->IntPoint(0).z = 0.0;
+
+   GeomVert[7]->IntPoint(1).x = 1.0;
+   GeomVert[7]->IntPoint(1).y = 0.0;
+   GeomVert[7]->IntPoint(1).z = 0.0;
+
+   GeomVert[7]->IntPoint(2).x = 1.0;
+   GeomVert[7]->IntPoint(2).y = 1.0;
+   GeomVert[7]->IntPoint(2).z = 0.0;
+
+   GeomVert[7]->IntPoint(3).x = 0.0;
+   GeomVert[7]->IntPoint(3).y = 1.0;
+   GeomVert[7]->IntPoint(3).z = 0.0;
+
+   GeomVert[7]->IntPoint(4).x = 0.0;
+   GeomVert[7]->IntPoint(4).y = 0.0;
+   GeomVert[7]->IntPoint(4).z = 1.0;
+
    GeomCenter[POINT].x = 0.0;
    GeomCenter[POINT].y = 0.0;
    GeomCenter[POINT].z = 0.0;
@@ -167,6 +193,10 @@ Geometry::Geometry()
    GeomCenter[PRISM].y = 1.0 / 3.0;
    GeomCenter[PRISM].z = 0.5;
 
+   GeomCenter[PYRAMID].x = 0.375;
+   GeomCenter[PYRAMID].y = 0.375;
+   GeomCenter[PYRAMID].z = 0.25;
+
    GeomToPerfGeomJac[POINT]       = NULL;
    GeomToPerfGeomJac[SEGMENT]     = new DenseMatrix(1);
    GeomToPerfGeomJac[TRIANGLE]    = new DenseMatrix(2);
@@ -174,6 +204,7 @@ Geometry::Geometry()
    GeomToPerfGeomJac[TETRAHEDRON] = new DenseMatrix(3);
    GeomToPerfGeomJac[CUBE]        = new DenseMatrix(3);
    GeomToPerfGeomJac[PRISM]       = new DenseMatrix(3);
+   GeomToPerfGeomJac[PYRAMID]     = new DenseMatrix(3);
 
    PerfGeomToGeomJac[POINT]       = NULL;
    PerfGeomToGeomJac[SEGMENT]     = NULL;
@@ -182,6 +213,7 @@ Geometry::Geometry()
    PerfGeomToGeomJac[TETRAHEDRON] = new DenseMatrix(3);
    PerfGeomToGeomJac[CUBE]        = NULL;
    PerfGeomToGeomJac[PRISM]       = new DenseMatrix(3);
+   PerfGeomToGeomJac[PYRAMID]     = new DenseMatrix(3);
 
    GeomToPerfGeomJac[SEGMENT]->Diag(1.0, 1);
    {
@@ -210,6 +242,14 @@ Geometry::Geometry()
       *GeomToPerfGeomJac[PRISM] = pri_T.Jacobian();
       CalcInverse(pri_T.Jacobian(), *PerfGeomToGeomJac[PRISM]);
    }
+   {
+      IsoparametricTransformation pyr_T;
+      pyr_T.SetFE(&PyramidFE);
+      GetPerfPointMat (PYRAMID, pyr_T.GetPointMat());
+      pyr_T.SetIntPoint(&GeomCenter[PYRAMID]);
+      *GeomToPerfGeomJac[PYRAMID] = pyr_T.Jacobian();
+      CalcInverse(pyr_T.Jacobian(), *PerfGeomToGeomJac[PYRAMID]);
+   }
 }
 
 Geometry::~Geometry()
@@ -233,6 +273,7 @@ const IntegrationRule * Geometry::GetVertices(int GeomType)
       case Geometry::TETRAHEDRON: return GeomVert[4];
       case Geometry::CUBE:        return GeomVert[5];
       case Geometry::PRISM:       return GeomVert[6];
+      case Geometry::PYRAMID:     return GeomVert[7];
       default:
          mfem_error ("Geometry::GetVertices(...)");
    }
@@ -310,6 +351,25 @@ void Geometry::GetRandomPoint(int GeomType, IntegrationPoint &ip)
             ip.y = 1.0 - ip.y;
          }
          break;
+      case Geometry::PYRAMID:
+         ip.x = double(rand()) / RAND_MAX;
+         ip.y = double(rand()) / RAND_MAX;
+         ip.z = double(rand()) / RAND_MAX;
+         if (ip.x + ip.z > 1.0 && ip.y < ip.x)
+         {
+            double x = ip.x;
+            ip.x = ip.y;
+            ip.y = 1.0 - ip.z;
+            ip.z = 1.0 - x;
+         }
+         else if (ip.y + ip.z > 1.0)
+         {
+            double z = ip.z;
+            ip.z = 1.0 - ip.y;
+            ip.y = ip.x;
+            ip.x = 1.0 - z;
+         }
+         break;
       default:
          MFEM_ABORT("Unknown type of reference element!");
    }
@@ -371,6 +431,10 @@ bool Geometry::CheckPoint(int GeomType, const IntegrationPoint &ip)
          if (ip.x < 0.0 || ip.y < 0.0 || ip.x+ip.y > 1.0 ||
              ip.z < 0.0 || ip.z > 1.0) { return false; }
          break;
+      case Geometry::PYRAMID:
+         if (ip.x < 0.0 || ip.y < 0.0 || ip.x+ip.z > 1.0 || ip.y+ip.z > 1.0 ||
+             ip.z < 0.0 || ip.z > 1.0) { return false; }
+         break;
       default:
          MFEM_ABORT("Unknown type of reference element!");
    }
@@ -441,6 +505,17 @@ bool Geometry::CheckPoint(int GeomType, const IntegrationPoint &ip, double eps)
             return false;
          }
          break;
+      case Geometry::PYRAMID:
+         if (internal::FuzzyLT(ip.x, 0.0, eps)
+             || internal::FuzzyLT(ip.y, 0.0, eps)
+             || internal::FuzzyGT(ip.x+ip.z, 1.0, eps)
+             || internal::FuzzyGT(ip.y+ip.z, 1.0, eps)
+             || internal::FuzzyLT(ip.z, 0.0, eps)
+             || internal::FuzzyGT(ip.z, 1.0, eps) )
+         {
+            return false;
+         }
+         break;
       default:
          MFEM_ABORT("Unknown type of reference element!");
    }
@@ -555,6 +630,16 @@ bool Geometry::ProjectPoint(int GeomType, const IntegrationPoint &beg,
          double lbeg[5] = { beg.x, beg.y, beg.z, 1.0-beg.x-beg.y, 1.0-beg.z };
          return internal::IntersectSegment<5,3>(lbeg, lend, end);
       }
+      case Geometry::PYRAMID:
+      {
+         double lend[6] = { end.x, end.y, end.z,
+                            1.0-end.x-end.z, 1.0-end.y-end.z, 1.0-end.z
+                          };
+         double lbeg[6] = { beg.x, beg.y, beg.z,
+                            1.0-beg.x-beg.z, 1.0-beg.y-beg.z, 1.0-beg.z
+                          };
+         return internal::IntersectSegment<6,3>(lbeg, lend, end);
+      }
       default:
          MFEM_ABORT("Unknown type of reference element!");
    }
@@ -652,6 +737,43 @@ bool Geometry::ProjectPoint(int GeomType, IntegrationPoint &ip)
          return in_tri && in_z;
       }
 
+      case PYRAMID:
+      {
+         if (ip.x < 0.0)
+         {
+            ip.x = 0.0;
+            internal::ProjectTriangle(ip.y, ip.z);
+            return false;
+         }
+         if (ip.y < 0.0)
+         {
+            ip.y = 0.0;
+            internal::ProjectTriangle(ip.x, ip.z);
+            return false;
+         }
+         if (ip.z < 0.0)
+         {
+            ip.z = 0.0;
+            if (ip.x > 1.0) { ip.x = 1.0; }
+            if (ip.y > 1.0) { ip.y = 1.0; }
+            return false;
+         }
+         if (ip.x >= ip.y)
+         {
+            bool in_y = true;
+            bool in_tri = internal::ProjectTriangle(ip.x, ip.z);
+            if (ip.y > ip.z) { in_y = false; ip.y = ip.z; }
+            return in_tri && in_y;
+         }
+         else
+         {
+            bool in_x = true;
+            bool in_tri = internal::ProjectTriangle(ip.y, ip.z);
+            if (ip.x > ip.z) { in_x = false; ip.x = ip.z; }
+            return in_tri && in_x;
+         }
+      }
+
       default:
          MFEM_ABORT("Reference element type is not supported!");
    }
@@ -726,6 +848,17 @@ void Geometry::GetPerfPointMat(int GeomType, DenseMatrix &pm)
       }
       break;
 
+      case Geometry::PYRAMID:
+      {
+         pm.SetSize (3, 5);
+         pm(0,0) = 0.0;  pm(1,0) = 0.0;  pm(2,0) = 0.0;
+         pm(0,1) = 1.0;  pm(1,1) = 0.0;  pm(2,1) = 0.0;
+         pm(0,2) = 1.0;  pm(1,2) = 1.0;  pm(2,2) = 0.0;
+         pm(0,3) = 0.0;  pm(1,3) = 1.0;  pm(2,3) = 0.0;
+         pm(0,4) = 0.5;  pm(1,4) = 0.5;  pm(2,4) = 0.7071067811865475;
+      }
+      break;
+
       default:
          mfem_error ("Geometry::GetPerfPointMat (...)");
    }
@@ -744,13 +877,13 @@ void Geometry::JacToPerfJac(int GeomType, const DenseMatrix &J,
    }
 }
 
-const int Geometry::NumBdrArray[NumGeom] = { 0, 2, 3, 4, 4, 6, 5 };
-const int Geometry::Dimension[NumGeom] = { 0, 1, 2, 2, 3, 3, 3 };
+const int Geometry::NumBdrArray[NumGeom] = { 0, 2, 3, 4, 4, 6, 5, 5 };
+const int Geometry::Dimension[NumGeom] = { 0, 1, 2, 2, 3, 3, 3, 3 };
 const int Geometry::DimStart[MaxDim+2] =
 { POINT, SEGMENT, TRIANGLE, TETRAHEDRON, NUM_GEOMETRIES };
-const int Geometry::NumVerts[NumGeom] = { 1, 2, 3, 4, 4, 8, 6 };
-const int Geometry::NumEdges[NumGeom] = { 0, 1, 3, 4, 6, 12, 9 };
-const int Geometry::NumFaces[NumGeom] = { 0, 0, 1, 1, 4, 6, 5 };
+const int Geometry::NumVerts[NumGeom] = { 1, 2, 3, 4, 4, 8, 6, 5 };
+const int Geometry::NumEdges[NumGeom] = { 0, 1, 3, 4, 6, 12, 9, 8 };
+const int Geometry::NumFaces[NumGeom] = { 0, 0, 1, 1, 4, 6, 5, 5 };
 
 const int Geometry::
 Constants<Geometry::POINT>::Orient[1][1] = {{0}};
@@ -897,6 +1030,30 @@ Constants<Geometry::PRISM>::VertToVert::J[9][2] =
    {5, 4}                   // 4,5:4
 };
 
+const int Geometry::
+Constants<Geometry::PYRAMID>::Edges[8][2] =
+{{0, 1}, {1, 2}, {3, 2}, {0, 3}, {0, 4}, {1, 4}, {2, 4}, {3, 4}};
+const int Geometry::
+Constants<Geometry::PYRAMID>::FaceTypes[5] =
+{
+   Geometry::SQUARE,
+   Geometry::TRIANGLE, Geometry::TRIANGLE,
+   Geometry::TRIANGLE, Geometry::TRIANGLE
+};
+const int Geometry::
+Constants<Geometry::PYRAMID>::FaceVert[5][4] =
+{{3, 2, 1, 0}, {0, 1, 4, -1}, {1, 2, 4, -1}, {2, 3, 4, -1}, {3, 0, 4, -1}};
+const int Geometry::
+Constants<Geometry::PYRAMID>::VertToVert::I[5] = {0, 3, 5, 7, 8};
+const int Geometry::
+Constants<Geometry::PYRAMID>::VertToVert::J[8][2] =
+{
+   {1, 0}, {3, 3}, {4, 4}, // 0,1:0   0,3:3  0,4:4
+   {2, 1}, {4, 5},         // 1,2:1   1,4:5
+   {3,-3}, {4, 6},         // 2,3:-3  2,4:6
+   {4, 7}                  // 3,4:7
+};
+
 
 GeometryRefiner::GeometryRefiner()
 {
@@ -1262,6 +1419,104 @@ RefinedGeometry * GeometryRefiner::Refine(Geometry::Type Geom,
          return RG;
       }
 
+      case Geometry::PYRAMID:
+      {
+         const int n = Times;
+         RG = new RefinedGeometry ((n+1)*(n+2)*(2*n+3)/6,
+                                   5*n*(2*n-1)*(2*n+1)/3, 0);
+         RG->Times = Times;
+         RG->ETimes = ETimes;
+         RG->Type = type;
+         // enumerate and define the vertices
+         m = 0;
+         for (k = 0; k <= n; k++)
+         {
+            const double *cpij =
+               poly1d.GetPoints(Times - k, BasisType::GetNodalBasis(type));
+            for (j = 0; j <= n - k; j++)
+               for (i = 0; i <= n - k; i++)
+               {
+                  IntegrationPoint &ip = RG->RefPts.IntPoint(m);
+                  if (type == 0)
+                  {
+                     ip.x = (n > k) ? (double(i) / (n - k)) : 0.0;
+                     ip.y = (n > k) ? (double(j) / (n - k)) : 0.0;
+                     ip.z = double(k) / n;
+                  }
+                  else
+                  {
+                     ip.x = cpij[i] * (1.0 - cp[k]);
+                     ip.y = cpij[j] * (1.0 - cp[k]);
+                     ip.z = cp[k];
+                  }
+                  m++;
+               }
+         }
+         if (m != (n+1)*(n+2)*(2*n+3)/6)
+         {
+            mfem_error("GeometryRefiner::Refine() for PYRAMID #1");
+         }
+         // elements
+         Array<int> &G = RG->RefGeoms;
+         m = 0;
+         for (k = 0; k < n; k++)
+         {
+            int lk = k * (k * (2 * k - 6 * n - 9) + 6 * n * (n + 3) + 13) / 6;
+            int lkp1 = (k + 1) *
+                       (k * (2 * k - 6 * n -5) + 6 * n * (n + 2) + 6) / 6;
+            for (j = 0; j < n - k; j++)
+            {
+               for (i = 0; i < n - k; i++)
+               {
+                  G[m++] = lk + j * (n - k + 1) + i;
+                  G[m++] = lk + j * (n - k + 1) + i + 1;
+                  G[m++] = lk + (j + 1) * (n - k + 1) + i + 1;
+                  G[m++] = lk + (j + 1) * (n - k + 1) + i;
+                  G[m++] = lkp1 + j * (n - k) + i;
+               }
+            }
+            for (j = 0; j < n - k - 1; j++)
+            {
+               for (i = 0; i < n - k - 1; i++)
+               {
+                  G[m++] = lkp1 + j * (n - k) + i;
+                  G[m++] = lkp1 + (j + 1) * (n - k) + i;
+                  G[m++] = lkp1 + (j + 1) * (n - k) + i + 1;
+                  G[m++] = lkp1 + j * (n - k) + i + 1;
+                  G[m++] = lk + (j + 1) * (n - k + 1) + i + 1;
+               }
+            }
+            for (j = 0; j < n - k; j++)
+            {
+               for (i = 0; i < n - k - 1; i++)
+               {
+                  G[m++] = lk + j * (n - k + 1) + i + 1;
+                  G[m++] = lk + (j + 1) * (n - k + 1) + i + 1;
+                  G[m++] = lkp1 + j * (n - k) + i;
+                  G[m++] = lkp1 + j * (n - k) + i + 1;
+                  G[m++] = -1;
+               }
+            }
+            for (j = 0; j < n - k - 1; j++)
+            {
+               for (i = 0; i < n - k; i++)
+               {
+                  G[m++] = lk + (j + 1) * (n - k + 1) + i;
+                  G[m++] = lk + (j + 1) * (n - k + 1) + i + 1;
+                  G[m++] = lkp1 + (j + 1) * (n - k) + i;
+                  G[m++] = lkp1 + j * (n - k) + i;
+                  G[m++] = -1;
+               }
+            }
+         }
+         if (m != 5*n*(2*n-1)*(2*n+1)/3)
+         {
+            mfem_error("GeometryRefiner::Refine() for PYRAMID #2");
+         }
+         RGeom[Geometry::PYRAMID].Append(RG);
+         return RG;
+      }
+
       case Geometry::PRISM:
       {
          const int n = Times;

fem/geom.hpp

@@ -27,6 +27,7 @@ namespace mfem
     Geometry::TETRAHEDRON - w/ vert. (0,0,0),(1,0,0),(0,1,0),(0,0,1)
     Geometry::CUBE - the unit cube
     Geometry::PRISM - w/ vert. (0,0,0),(1,0,0),(0,1,0),(0,0,1),(1,0,1),(0,1,1)
+    Geometry::PYRAMID - w/ vert. (0,0,0),(1,0,0),(1,1,0),(0,1,0),(0,0,1)
 */
 class Geometry
 {
@@ -34,7 +35,7 @@ public:
    enum Type
    {
       INVALID = -1,
-      POINT = 0, SEGMENT, TRIANGLE, SQUARE, TETRAHEDRON, CUBE, PRISM,
+      POINT = 0, SEGMENT, TRIANGLE, SQUARE, TETRAHEDRON, CUBE, PRISM, PYRAMID,
       NUM_GEOMETRIES
    };
 
@@ -251,7 +252,26 @@ template <> struct Geometry::Constants<Geometry::PRISM>
    };
 };
 
-// Defined in fe.cpp to ensure construction after 'mfem::WedgeFE'.
+template <> struct Geometry::Constants<Geometry::PYRAMID>
+{
+   static const int Dimension = 3;
+   static const int NumVert = 5;
+   static const int NumEdges = 8;
+   static const int Edges[NumEdges][2];
+   static const int NumFaces = 5;
+   static const int FaceTypes[NumFaces];
+   static const int MaxFaceVert = 4;
+   static const int FaceVert[NumFaces][MaxFaceVert];
+   // Upper-triangular part of the local vertex-to-vertex graph.
+   struct VertToVert
+   {
+      static const int I[NumVert];
+      static const int J[NumEdges][2]; // {end,edge_idx}
+   };
+};
+
+// Defined in fe.cpp to ensure construction after 'mfem::TriangleFE' and
+// `mfem::TetrahedronFE`.
 extern Geometry Geometries;
 
 

fem/intrules.cpp

@@ -910,6 +910,9 @@ IntegrationRules::IntegrationRules(int Ref, int type_):
    TetrahedronIntRules.SetSize(32, h_mt);
    TetrahedronIntRules = NULL;
 
+   PyramidIntRules.SetSize(32, h_mt);
+   PyramidIntRules = NULL;
+
    PrismIntRules.SetSize(32, h_mt);
    PrismIntRules = NULL;
 
@@ -930,6 +933,7 @@ const IntegrationRule &IntegrationRules::Get(int GeomType, int Order)
       case Geometry::TETRAHEDRON: ir_array = &TetrahedronIntRules; break;
       case Geometry::CUBE:        ir_array = &CubeIntRules; break;
       case Geometry::PRISM:       ir_array = &PrismIntRules; break;
+      case Geometry::PYRAMID:     ir_array = &PyramidIntRules; break;
       default:
          mfem_error("IntegrationRules::Get(...) : Unknown geometry type!");
          ir_array = NULL;
@@ -976,6 +980,7 @@ void IntegrationRules::Set(int GeomType, int Order, IntegrationRule &IntRule)
       case Geometry::TETRAHEDRON: ir_array = &TetrahedronIntRules; break;
       case Geometry::CUBE:        ir_array = &CubeIntRules; break;
       case Geometry::PRISM:       ir_array = &PrismIntRules; break;
+      case Geometry::PYRAMID:     ir_array = &PyramidIntRules; break;
       default:
          mfem_error("IntegrationRules::Set(...) : Unknown geometry type!");
          ir_array = NULL;
@@ -1019,6 +1024,7 @@ IntegrationRules::~IntegrationRules()
    DeleteIntRuleArray(TetrahedronIntRules);
    DeleteIntRuleArray(CubeIntRules);
    DeleteIntRuleArray(PrismIntRules);
+   DeleteIntRuleArray(PyramidIntRules);
 }
 
 
@@ -1041,6 +1047,8 @@ IntegrationRule *IntegrationRules::GenerateIntegrationRule(int GeomType,
          return CubeIntegrationRule(Order);
       case Geometry::PRISM:
          return PrismIntegrationRule(Order);
+      case Geometry::PYRAMID:
+         return PyramidIntegrationRule(Order);
       default:
          mfem_error("IntegrationRules::Set(...) : Unknown geometry type!");
          return NULL;
@@ -1648,6 +1656,30 @@ IntegrationRule *IntegrationRules::TetrahedronIntegrationRule(int Order)
    }
 }
 
+// Integration rules for reference pyramid
+IntegrationRule *IntegrationRules::PyramidIntegrationRule(int Order)
+{
+   // This is a simple integration rule adapted from an integration
+   // rule for a cube which seems to be adequate for now. When we
+   // implement high order finite elements for pyramids we should
+   // revisit this and see if we can improve upon it.
+   const IntegrationRule &irc = Get(Geometry::CUBE, Order);
+   int npts = irc.GetNPoints();
+   AllocIntRule(PyramidIntRules, Order);
+   PyramidIntRules[Order] = new IntegrationRule(npts);
+
+   for (int k=0; k<npts; k++)
+   {
+      const IntegrationPoint & ipc = irc.IntPoint(k);
+      IntegrationPoint & ipp = PyramidIntRules[Order]->IntPoint(k);
+      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;
+   }
+   return PyramidIntRules[Order];
+}
+
 // Integration rules for reference prism
 IntegrationRule *IntegrationRules::PrismIntegrationRule(int Order)
 {

fem/intrules.hpp

@@ -323,6 +323,7 @@ private:
    Array<IntegrationRule *> TriangleIntRules;
    Array<IntegrationRule *> SquareIntRules;
    Array<IntegrationRule *> TetrahedronIntRules;
+   Array<IntegrationRule *> PyramidIntRules;
    Array<IntegrationRule *> PrismIntRules;
    Array<IntegrationRule *> CubeIntRules;
 
@@ -351,6 +352,7 @@ private:
    IntegrationRule *TriangleIntegrationRule(int Order);
    IntegrationRule *SquareIntegrationRule(int Order);
    IntegrationRule *TetrahedronIntegrationRule(int Order);
+   IntegrationRule *PyramidIntegrationRule(int Order);
    IntegrationRule *PrismIntegrationRule(int Order);
    IntegrationRule *CubeIntegrationRule(int Order);
 

fem/nonlininteg_vectorconvection.cpp

@@ -116,15 +116,15 @@ void VectorConvectionNLFIntegrator::AssemblePA(const FiniteElementSpace &fes)
 
 // PA Convection NL 2D kernel
 template<int T_D1D = 0, int T_Q1D = 0>
-static void PAConvectionNLApply2D(const int NE,
-                                  const Array<double> &b,
-                                  const Array<double> &g,
-                                  const Array<double> &bt,
-                                  const Vector &q_,
-                                  const Vector &x_,
-                                  Vector &y_,
-                                  const int d1d = 0,
-                                  const int q1d = 0)
+void PAConvectionNLApply2D(const int NE,
+                           const Array<double> &b,
+                           const Array<double> &g,
+                           const Array<double> &bt,
+                           const Vector &q_,
+                           const Vector &x_,
+                           Vector &y_,
+                           const int d1d = 0,
+                           const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -252,15 +252,15 @@ static void PAConvectionNLApply2D(const int NE,
 
 // PA Convection NL 3D kernel
 template<int T_D1D = 0, int T_Q1D = 0>
-static void PAConvectionNLApply3D(const int NE,
-                                  const Array<double> &b,
-                                  const Array<double> &g,
-                                  const Array<double> &bt,
-                                  const Vector &q_,
-                                  const Vector &x_,
-                                  Vector &y_,
-                                  const int d1d = 0,
-                                  const int q1d = 0)
+void PAConvectionNLApply3D(const int NE,
+                           const Array<double> &b,
+                           const Array<double> &g,
+                           const Array<double> &bt,
+                           const Vector &q_,
+                           const Vector &x_,
+                           Vector &y_,
+                           const int d1d = 0,
+                           const int q1d = 0)
 {
    constexpr int VDIM = 3;
    const int D1D = T_D1D ? T_D1D : d1d;
@@ -558,14 +558,14 @@ static void PAConvectionNLApply3D(const int NE,
 }
 
 template<int T_D1D = 0, int T_Q1D = 0, int T_MAX_D1D =0, int T_MAX_Q1D =0>
-static void SmemPAConvectionNLApply3D(const int NE,
-                                      const Array<double> &b_,
-                                      const Array<double> &g_,
-                                      const Vector &d_,
-                                      const Vector &x_,
-                                      Vector &y_,
-                                      const int d1d = 0,
-                                      const int q1d = 0)
+void SmemPAConvectionNLApply3D(const int NE,
+                               const Array<double> &b_,
+                               const Array<double> &g_,
+                               const Vector &d_,
+                               const Vector &x_,
+                               Vector &y_,
+                               const int d1d = 0,
+                               const int q1d = 0)
 {
    constexpr int VDIM = 3;
    const int D1D = T_D1D ? T_D1D : d1d;

fem/qinterp/det.cpp

@@ -27,14 +27,14 @@ namespace quadrature_interpolator
 {
 
 template<int T_D1D = 0, int T_Q1D = 0, int MAX_D1D = 0, int MAX_Q1D = 0>
-static void Det2D(const int NE,
-                  const double *b,
-                  const double *g,
-                  const double *x,
-                  double *y,
-                  const int vdim = 1,
-                  const int d1d = 0,
-                  const int q1d = 0)
+void Det2D(const int NE,
+           const double *b,
+           const double *g,
+           const double *x,
+           double *y,
+           const int vdim = 1,
+           const int d1d = 0,
+           const int q1d = 0)
 {
    constexpr int DIM = 2;
    static constexpr int NBZ = 1;
@@ -79,15 +79,15 @@ static void Det2D(const int NE,
 
 template<int T_D1D = 0, int T_Q1D = 0, int MAX_D1D = 0, int MAX_Q1D = 0,
          bool SMEM = true>
-static void Det3D(const int NE,
-                  const double *b,
-                  const double *g,
-                  const double *x,
-                  double *y,
-                  const int vdim = 1,
-                  const int d1d = 0,
-                  const int q1d = 0,
-                  Vector *d_buff = nullptr) // used only with SMEM = false
+void Det3D(const int NE,
+           const double *b,
+           const double *g,
+           const double *x,
+           double *y,
+           const int vdim = 1,
+           const int d1d = 0,
+           const int q1d = 0,
+           Vector *d_buff = nullptr) // used only with SMEM = false
 {
    constexpr int DIM = 3;
    static constexpr int MQ1 = T_Q1D ? T_Q1D : MAX_Q1D;

fem/qinterp/eval.hpp

@@ -31,13 +31,13 @@ namespace quadrature_interpolator
 template<QVectorLayout Q_LAYOUT,
          int T_VDIM = 0, int T_D1D = 0, int T_Q1D = 0,
          int T_NBZ = 1, int MAX_D1D = 0, int MAX_Q1D = 0>
-static void Values2D(const int NE,
-                     const double *b_,
-                     const double *x_,
-                     double *y_,
-                     const int vdim = 0,
-                     const int d1d = 0,
-                     const int q1d = 0)
+void Values2D(const int NE,
+              const double *b_,
+              const double *x_,
+              double *y_,
+              const int vdim = 0,
+              const int d1d = 0,
+              const int q1d = 0)
 {
    static constexpr int NBZ = T_NBZ ? T_NBZ : 1;
 
@@ -95,13 +95,13 @@ static void Values2D(const int NE,
 template<QVectorLayout Q_LAYOUT,
          int T_VDIM = 0, int T_D1D = 0, int T_Q1D = 0,
          int MAX_D1D = 0, int MAX_Q1D = 0>
-static void Values3D(const int NE,
-                     const double *b_,
-                     const double *x_,
-                     double *y_,
-                     const int vdim = 0,
-                     const int d1d = 0,
-                     const int q1d = 0)
+void Values3D(const int NE,
+              const double *b_,
+              const double *x_,
+              double *y_,
+              const int vdim = 0,
+              const int d1d = 0,
+              const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;

fem/qinterp/grad.hpp

@@ -31,15 +31,15 @@ namespace quadrature_interpolator
 template<QVectorLayout Q_LAYOUT, bool GRAD_PHYS,
          int T_VDIM = 0, int T_D1D = 0, int T_Q1D = 0,
          int T_NBZ = 1, int MAX_D1D = 0, int MAX_Q1D = 0>
-static void Derivatives2D(const int NE,
-                          const double *b_,
-                          const double *g_,
-                          const double *j_,
-                          const double *x_,
-                          double *y_,
-                          const int vdim = 0,
-                          const int d1d = 0,
-                          const int q1d = 0)
+void Derivatives2D(const int NE,
+                   const double *b_,
+                   const double *g_,
+                   const double *j_,
+                   const double *x_,
+                   double *y_,
+                   const int vdim = 0,
+                   const int d1d = 0,
+                   const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -139,15 +139,15 @@ static void Derivatives2D(const int NE,
 template<QVectorLayout Q_LAYOUT, bool GRAD_PHYS,
          int T_VDIM = 0, int T_D1D = 0, int T_Q1D = 0,
          int MAX_D1D = 0, int MAX_Q1D = 0>
-static void Derivatives3D(const int NE,
-                          const double *b_,
-                          const double *g_,
-                          const double *j_,
-                          const double *x_,
-                          double *y_,
-                          const int vdim = 0,
-                          const int d1d = 0,
-                          const int q1d = 0)
+void Derivatives3D(const int NE,
+                   const double *b_,
+                   const double *g_,
+                   const double *j_,
+                   const double *x_,
+                   double *y_,
+                   const int vdim = 0,
+                   const int d1d = 0,
+                   const int q1d = 0)
 {
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;

fem/quadinterpolator.cpp

@@ -61,16 +61,16 @@ namespace quadrature_interpolator
 // * assumes 'e_vec' is using ElementDofOrdering::NATIVE,
 // * assumes 'maps.mode == FULL'.
 template<const int T_VDIM, const int T_ND, const int T_NQ>
-static void Eval2D(const int NE,
-                   const int vdim,
-                   const QVectorLayout q_layout,
-                   const GeometricFactors *geom,
-                   const DofToQuad &maps,
-                   const Vector &e_vec,
-                   Vector &q_val,
-                   Vector &q_der,
-                   Vector &q_det,
-                   const int eval_flags)
+void Eval2D(const int NE,
+            const int vdim,
+            const QVectorLayout q_layout,
+            const GeometricFactors *geom,
+            const DofToQuad &maps,
+            const Vector &e_vec,
+            Vector &q_val,
+            Vector &q_der,
+            Vector &q_det,
+            const int eval_flags)
 {
    using QI = QuadratureInterpolator;
 
@@ -209,16 +209,16 @@ static void Eval2D(const int NE,
 // * assumes 'e_vec' is using ElementDofOrdering::NATIVE,
 // * assumes 'maps.mode == FULL'.
 template<const int T_VDIM, const int T_ND, const int T_NQ>
-static void Eval3D(const int NE,
-                   const int vdim,
-                   const QVectorLayout q_layout,
-                   const GeometricFactors *geom,
-                   const DofToQuad &maps,
-                   const Vector &e_vec,
-                   Vector &q_val,
-                   Vector &q_der,
-                   Vector &q_det,
-                   const int eval_flags)
+void Eval3D(const int NE,
+            const int vdim,
+            const QVectorLayout q_layout,
+            const GeometricFactors *geom,
+            const DofToQuad &maps,
+            const Vector &e_vec,
+            Vector &q_val,
+            Vector &q_der,
+            Vector &q_det,
+            const int eval_flags)
 {
    using QI = QuadratureInterpolator;
 

fem/tmop.cpp

@@ -1306,7 +1306,7 @@ namespace internal
 // MFEM_FORALL-based copy kernel -- used by protected methods below.
 // Needed as a workaround for the nvcc restriction that methods with MFEM_FORALL
 // in them must to be public.
-static inline void device_copy(double *d_dest, const double *d_src, int size)
+inline void device_copy(double *d_dest, const double *d_src, int size)
 {
    MFEM_FORALL(i, size, d_dest[i] = d_src[i];);
 }

fem/tmop_amr.cpp

@@ -234,7 +234,7 @@ void TMOPRefinerEstimator::SetTriIntRules()
 
    // Reftype = 0 // original element
    const int Nvert = 3, NEsplit = 1;
-   Mesh meshsplit(2, Nvert, NEsplit, 0 ,2);
+   Mesh meshsplit(2, Nvert, NEsplit, 0,2);
    const double tri_v[3][2] =
    {
       {0, 0}, {1, 0}, {0, 1}

general/error.hpp

@@ -164,7 +164,13 @@ __device__ void abort_msg(T & msg)
 #endif
 
 // Abort inside a device kernel
-#if defined(__CUDA_ARCH__)
+#if defined(__CUDA_ARCH__) && defined(_WIN32)
+#define MFEM_ABORT_KERNEL(msg) \
+   {                           \
+      printf(msg);             \
+      __debugbreak();          \
+   }
+#elif defined(__CUDA_ARCH__)
 #define MFEM_ABORT_KERNEL(msg) \
    {                           \
       printf(msg);             \

general/isockstream.cpp

@@ -79,7 +79,13 @@ int isockstream::establish()
    int on=1;
    setsockopt(port, SOL_SOCKET, SO_REUSEADDR, (char *)(&on), sizeof(on));
 
-   if (bind(port,(const sockaddr*)&sa,(socklen_t)sizeof(struct sockaddr_in)) < 0)
+   if (bind(
+#ifdef _WIN32
+          (SOCKET)port
+#else
+          port
+#endif
+          ,(const sockaddr*)&sa,(socklen_t)sizeof(struct sockaddr_in)) < 0)
    {
       mfem::err << "isockstream::establish(): bind() failed!" << endl;
       close(port);

linalg/kernels.hpp

@@ -378,7 +378,7 @@ void Swap(T &a, T &b)
    b = tmp;
 }
 
-const double Epsilon = std::numeric_limits<double>::epsilon();
+constexpr double Epsilon = std::numeric_limits<double>::epsilon();
 
 /// Utility function used in CalcSingularvalue<3>.
 MFEM_HOST_DEVICE static inline

linalg/sundials.cpp

@@ -1677,7 +1677,7 @@ int KINSolver::GradientMult(N_Vector v, N_Vector Jv, N_Vector u,
 
 // Wrapper for evaluating linear systems J u = b
 int KINSolver::LinSysSetup(N_Vector u, N_Vector, SUNMatrix J,
-                           void *, N_Vector , N_Vector )
+                           void *, N_Vector, N_Vector )
 {
    const SundialsNVector mfem_u(u);
    KINSolver *self = static_cast<KINSolver*>(GET_CONTENT(J));

mesh/CMakeLists.txt

@@ -19,6 +19,7 @@ set(SRCS
   ncmesh.cpp
   nurbs.cpp
   point.cpp
+  pyramid.cpp
   quadrilateral.cpp
   segment.cpp
   tetrahedron.cpp
@@ -38,6 +39,7 @@ set(HDRS
   ncmesh.hpp
   nurbs.hpp
   point.hpp
+  pyramid.hpp
   quadrilateral.hpp
   segment.hpp
   tetrahedron.hpp

mesh/element.hpp

@@ -39,7 +39,7 @@ public:
 
    /// Constants for the classes derived from Element.
    enum Type { POINT, SEGMENT, TRIANGLE, QUADRILATERAL,
-               TETRAHEDRON, HEXAHEDRON, WEDGE
+               TETRAHEDRON, HEXAHEDRON, WEDGE, PYRAMID
              };
 
    /// Default element constructor.

mesh/mesh.cpp

@@ -335,6 +335,7 @@ FiniteElement *Mesh::GetTransformationFEforElementType(Element::Type ElemType)
       case Element::TETRAHEDRON :    return &TetrahedronFE;
       case Element::HEXAHEDRON :     return &HexahedronFE;
       case Element::WEDGE :          return &WedgeFE;
+      case Element::PYRAMID :        return &PyramidFE;
       default:
          MFEM_ABORT("Unknown element type \"" << ElemType << "\"");
          break;
@@ -735,6 +736,31 @@ void Mesh::GetLocalTriToWdgTransformation(
    }
 }
 
+void Mesh::GetLocalTriToPyrTransformation(
+   IsoparametricTransformation &Transf, int i)
+{
+   DenseMatrix &locpm = Transf.GetPointMat();
+
+   Transf.SetFE(&TriangleFE);
+   //  (i/64) is the local face no. in the pyr
+   MFEM_VERIFY(i >= 64, "Local face index " << i/64
+               << " is not a triangular face of a pyramid.");
+   const int *pv = pyr_t::FaceVert[i/64];
+   //  (i%64) is the orientation of the pyramid face
+   //         w.r.t. the face element
+   const int *to = tri_t::Orient[i%64];
+   const IntegrationRule *PyrVert =
+      Geometries.GetVertices(Geometry::PYRAMID);
+   locpm.SetSize(3, 3);
+   for (int j = 0; j < 3; j++)
+   {
+      const IntegrationPoint &vert = PyrVert->IntPoint(pv[to[j]]);
+      locpm(0, j) = vert.x;
+      locpm(1, j) = vert.y;
+      locpm(2, j) = vert.z;
+   }
+}
+
 void Mesh::GetLocalQuadToHexTransformation(
    IsoparametricTransformation &Transf, int i)
 {
@@ -781,6 +807,29 @@ void Mesh::GetLocalQuadToWdgTransformation(
    }
 }
 
+void Mesh::GetLocalQuadToPyrTransformation(
+   IsoparametricTransformation &Transf, int i)
+{
+   DenseMatrix &locpm = Transf.GetPointMat();
+
+   Transf.SetFE(&QuadrilateralFE);
+   //  (i/64) is the local face no. in the pyr
+   MFEM_VERIFY(i < 64, "Local face index " << i/64
+               << " is not a quadrilateral face of a pyramid.");
+   const int *pv = pyr_t::FaceVert[i/64];
+   //  (i%64) is the orientation of the quad
+   const int *qo = quad_t::Orient[i%64];
+   const IntegrationRule *PyrVert = Geometries.GetVertices(Geometry::PYRAMID);
+   locpm.SetSize(3, 4);
+   for (int j = 0; j < 4; j++)
+   {
+      const IntegrationPoint &vert = PyrVert->IntPoint(pv[qo[j]]);
+      locpm(0, j) = vert.x;
+      locpm(1, j) = vert.y;
+      locpm(2, j) = vert.z;
+   }
+}
+
 const GeometricFactors* Mesh::GetGeometricFactors(const IntegrationRule& ir,
                                                   const int flags,
                                                   MemoryType d_mt)
@@ -862,10 +911,19 @@ void Mesh::GetLocalFaceTransformation(
          {
             GetLocalTriToTetTransformation(Transf, info);
          }
+         else if (elem_type == Element::WEDGE)
+         {
+            GetLocalTriToWdgTransformation(Transf, info);
+         }
+         else if (elem_type == Element::PYRAMID)
+         {
+            GetLocalTriToPyrTransformation(Transf, info);
+         }
          else
          {
-            MFEM_ASSERT(elem_type == Element::WEDGE, "");
-            GetLocalTriToWdgTransformation(Transf, info);
+            MFEM_ABORT("Mesh::GetLocalFaceTransformation not defined for "
+                       "face type " << face_type
+                       << " and element type " << elem_type << "\n");
          }
          break;
 
@@ -874,10 +932,19 @@ void Mesh::GetLocalFaceTransformation(
          {
             GetLocalQuadToHexTransformation(Transf, info);
          }
+         else if (elem_type == Element::WEDGE)
+         {
+            GetLocalQuadToWdgTransformation(Transf, info);
+         }
+         else if (elem_type == Element::PYRAMID)
+         {
+            GetLocalQuadToPyrTransformation(Transf, info);
+         }
          else
          {
-            MFEM_ASSERT(elem_type == Element::WEDGE, "");
-            GetLocalQuadToWdgTransformation(Transf, info);
+            MFEM_ABORT("Mesh::GetLocalFaceTransformation not defined for "
+                       "face type " << face_type
+                       << " and element type " << elem_type << "\n");
          }
          break;
    }
@@ -1370,6 +1437,20 @@ int Mesh::AddWedge(const int *vi, int attr)
    return NumOfElements++;
 }
 
+int Mesh::AddPyramid(int v1, int v2, int v3, int v4, int v5, int attr)
+{
+   CheckEnlarge(elements, NumOfElements);
+   elements[NumOfElements] = new Pyramid(v1, v2, v3, v4, v5, attr);
+   return NumOfElements++;
+}
+
+int Mesh::AddPyramid(const int *vi, int attr)
+{
+   CheckEnlarge(elements, NumOfElements);
+   elements[NumOfElements] = new Pyramid(vi, attr);
+   return NumOfElements++;
+}
+
 int Mesh::AddHex(int v1, int v2, int v3, int v4, int v5, int v6, int v7, int v8,
                  int attr)
 {
@@ -1423,6 +1504,25 @@ void Mesh::AddHexAsWedges(const int *vi, int attr)
    }
 }
 
+void Mesh::AddHexAsPyramids(const int *vi, int attr)
+{
+   static const int hex_to_pyr[6][5] =
+   {
+      { 0, 1, 2, 3, 8 }, { 0, 4, 5, 1, 8 }, { 1, 5, 6, 2, 8 },
+      { 2, 6, 7, 3, 8 }, { 3, 7, 4, 0, 8 }, { 7, 6, 5, 4, 8 }
+   };
+   int ti[5];
+
+   for (int i = 0; i < 6; i++)
+   {
+      for (int j = 0; j < 5; j++)
+      {
+         ti[j] = vi[hex_to_pyr[i][j]];
+      }
+      AddPyramid(ti, attr);
+   }
+}
+
 int Mesh::AddElement(Element *elem)
 {
    CheckEnlarge(elements, NumOfElements);
@@ -2692,11 +2792,16 @@ void Mesh::Make3D(int nx, int ny, int nz, Element::Type type,
       NElem *= 2;
       NBdrElem += 2*nx*ny;
    }
+   else if (type == Element::PYRAMID)
+   {
+      NElem *= 6;
+      NVert += nx * ny * nz;
+   }
 
    InitMesh(3, 3, NVert, NElem, NBdrElem);
 
    double coord[3];
-   int ind[8];
+   int ind[9];
 
    // Sets vertices and the corresponding coordinates
    for (z = 0; z <= nz; z++)
@@ -2712,8 +2817,25 @@ void Mesh::Make3D(int nx, int ny, int nz, Element::Type type,
          }
       }
    }
+   if (type == Element::PYRAMID)
+   {
+      for (z = 0; z < nz; z++)
+      {
+         coord[2] = (((double) z + 0.5) / nz) * sz;
+         for (y = 0; y < ny; y++)
+         {
+            coord[1] = (((double) y + 0.5 ) / ny) * sy;
+            for (x = 0; x < nx; x++)
+            {
+               coord[0] = (((double) x + 0.5 ) / nx) * sx;
+               AddVertex(coord);
+            }
+         }
+      }
+   }
 
 #define VTX(XC, YC, ZC) ((XC)+((YC)+(ZC)*(ny+1))*(nx+1))
+#define VTXP(XC, YC, ZC) ((nx+1)*(ny+1)*(nz+1)+(XC)+((YC)+(ZC)*ny)*nx)
 
    // Sets elements and the corresponding indices of vertices
    if (sfc_ordering && type == Element::HEXAHEDRON)
@@ -2728,14 +2850,14 @@ void Mesh::Make3D(int nx, int ny, int nz, Element::Type type,
          y = sfc[3*k + 1];
          z = sfc[3*k + 2];
 
-         ind[0] = VTX(x  , y  , z  );
-         ind[1] = VTX(x+1, y  , z  );
+         ind[0] = VTX(x, y, z  );
+         ind[1] = VTX(x+1, y, z  );
          ind[2] = VTX(x+1, y+1, z  );
-         ind[3] = VTX(x  , y+1, z  );
-         ind[4] = VTX(x  , y  , z+1);
-         ind[5] = VTX(x+1, y  , z+1);
+         ind[3] = VTX(x, y+1, z  );
+         ind[4] = VTX(x, y, z+1);
+         ind[5] = VTX(x+1, y, z+1);
          ind[6] = VTX(x+1, y+1, z+1);
-         ind[7] = VTX(x  , y+1, z+1);
+         ind[7] = VTX(x, y+1, z+1);
 
          AddHex(ind, 1);
       }
@@ -2748,12 +2870,12 @@ void Mesh::Make3D(int nx, int ny, int nz, Element::Type type,
          {
             for (x = 0; x < nx; x++)
             {
-               ind[0] = VTX(x  , y  , z  );
-               ind[1] = VTX(x+1, y  , z  );
+               ind[0] = VTX(x, y, z  );
+               ind[1] = VTX(x+1, y, z  );
                ind[2] = VTX(x+1, y+1, z  );
-               ind[3] = VTX(x  , y+1, z  );
-               ind[4] = VTX(x  , y  , z+1);
-               ind[5] = VTX(x+1, y  , z+1);
+               ind[3] = VTX(x, y+1, z  );
+               ind[4] = VTX(x, y, z+1);
+               ind[5] = VTX(x+1, y, z+1);
                ind[6] = VTX(x+1, y+1, z+1);
                ind[7] = VTX(  x, y+1, z+1);
                if (type == Element::TETRAHEDRON)
@@ -2764,6 +2886,11 @@ void Mesh::Make3D(int nx, int ny, int nz, Element::Type type,
                {
                   AddHexAsWedges(ind, 1);
                }
+               else if (type == Element::PYRAMID)
+               {
+                  ind[8] = VTXP(  x, y, z);
+                  AddHexAsPyramids(ind, 1);
+               }
                else
                {
                   AddHex(ind, 1);
@@ -2779,10 +2906,10 @@ void Mesh::Make3D(int nx, int ny, int nz, Element::Type type,
    {
       for (x = 0; x < nx; x++)
       {
-         ind[0] = VTX(x  , y  , 0);
-         ind[1] = VTX(x  , y+1, 0);
+         ind[0] = VTX(x, y, 0);
+         ind[1] = VTX(x, y+1, 0);
          ind[2] = VTX(x+1, y+1, 0);
-         ind[3] = VTX(x+1, y  , 0);
+         ind[3] = VTX(x+1, y, 0);
          if (type == Element::TETRAHEDRON)
          {
             AddBdrQuadAsTriangles(ind, 1);
@@ -2802,10 +2929,10 @@ void Mesh::Make3D(int nx, int ny, int nz, Element::Type type,
    {
       for (x = 0; x < nx; x++)
       {
-         ind[0] = VTX(x  , y  , nz);
-         ind[1] = VTX(x+1, y  , nz);
+         ind[0] = VTX(x, y, nz);
+         ind[1] = VTX(x+1, y, nz);
          ind[2] = VTX(x+1, y+1, nz);
-         ind[3] = VTX(x  , y+1, nz);
+         ind[3] = VTX(x, y+1, nz);
          if (type == Element::TETRAHEDRON)
          {
             AddBdrQuadAsTriangles(ind, 6);
@@ -2825,10 +2952,10 @@ void Mesh::Make3D(int nx, int ny, int nz, Element::Type type,
    {
       for (y = 0; y < ny; y++)
       {
-         ind[0] = VTX(0  , y  , z  );
-         ind[1] = VTX(0  , y  , z+1);
-         ind[2] = VTX(0  , y+1, z+1);
-         ind[3] = VTX(0  , y+1, z  );
+         ind[0] = VTX(0, y, z  );
+         ind[1] = VTX(0, y, z+1);
+         ind[2] = VTX(0, y+1, z+1);
+         ind[3] = VTX(0, y+1, z  );
          if (type == Element::TETRAHEDRON)
          {
             AddBdrQuadAsTriangles(ind, 5);
@@ -2844,10 +2971,10 @@ void Mesh::Make3D(int nx, int ny, int nz, Element::Type type,
    {
       for (y = 0; y < ny; y++)
       {
-         ind[0] = VTX(nx, y  , z  );
+         ind[0] = VTX(nx, y, z  );
          ind[1] = VTX(nx, y+1, z  );
          ind[2] = VTX(nx, y+1, z+1);
-         ind[3] = VTX(nx, y  , z+1);
+         ind[3] = VTX(nx, y, z+1);
          if (type == Element::TETRAHEDRON)
          {
             AddBdrQuadAsTriangles(ind, 3);
@@ -2863,10 +2990,10 @@ void Mesh::Make3D(int nx, int ny, int nz, Element::Type type,
    {
       for (z = 0; z < nz; z++)
       {
-         ind[0] = VTX(x  , 0, z  );
+         ind[0] = VTX(x, 0, z  );
          ind[1] = VTX(x+1, 0, z  );
          ind[2] = VTX(x+1, 0, z+1);
-         ind[3] = VTX(x  , 0, z+1);
+         ind[3] = VTX(x, 0, z+1);
          if (type == Element::TETRAHEDRON)
          {
             AddBdrQuadAsTriangles(ind, 2);
@@ -2882,8 +3009,8 @@ void Mesh::Make3D(int nx, int ny, int nz, Element::Type type,
    {
       for (z = 0; z < nz; z++)
       {
-         ind[0] = VTX(x  , ny, z  );
-         ind[1] = VTX(x  , ny, z+1);
+         ind[0] = VTX(x, ny, z  );
+         ind[1] = VTX(x, ny, z+1);
          ind[2] = VTX(x+1, ny, z+1);
          ind[3] = VTX(x+1, ny, z  );
          if (type == Element::TETRAHEDRON)
@@ -3428,6 +3555,7 @@ Element *Mesh::NewElement(int geom)
 #endif
       case Geometry::CUBE:      return (new Hexahedron);
       case Geometry::PRISM:     return (new Wedge);
+      case Geometry::PYRAMID:   return (new Pyramid);
       default:
          MFEM_ABORT("invalid Geometry::Type, geom = " << geom);
    }
@@ -3520,6 +3648,15 @@ void Mesh::SetMeshGen()
             meshgen |= 4;
             break;
 
+         case Element::PYRAMID:
+            mesh_geoms |= (1 << Geometry::PYRAMID);
+            mesh_geoms |= (1 << Geometry::SQUARE);
+            mesh_geoms |= (1 << Geometry::TRIANGLE);
+            mesh_geoms |= (1 << Geometry::SEGMENT);
+            mesh_geoms |= (1 << Geometry::POINT);
+            meshgen |= 8;
+            break;
+
          default:
             MFEM_ABORT("invalid element type: " << type);
             break;
@@ -5060,6 +5197,19 @@ int Mesh::CheckElementOrientation(bool fix_it)
                }
                break;
 
+            case Element::PYRAMID:
+               // only check the Jacobian at the center of the element
+               GetElementJacobian(i, J);
+               if (J.Det() < 0.0)
+               {
+                  wo++;
+                  if (fix_it)
+                  {
+                     // how?
+                  }
+               }
+               break;
+
             case Element::HEXAHEDRON:
                // only check the Jacobian at the center of the element
                GetElementJacobian(i, J);
@@ -6142,6 +6292,22 @@ void Mesh::GenerateFaces()
                }
                break;
             }
+            case Element::PYRAMID:
+            {
+               for (int j = 0; j < 1; j++)
+               {
+                  const int *fv = pyr_t::FaceVert[j];
+                  AddQuadFaceElement(j, ef[j], i,
+                                     v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
+               }
+               for (int j = 1; j < 5; j++)
+               {
+                  const int *fv = pyr_t::FaceVert[j];
+                  AddTriangleFaceElement(j, ef[j], i,
+                                         v[fv[0]], v[fv[1]], v[fv[2]]);
+               }
+               break;
+            }
             case Element::HEXAHEDRON:
             {
                for (int j = 0; j < 6; j++)
@@ -6234,6 +6400,20 @@ STable3D *Mesh::GetFacesTable()
             }
             break;
          }
+         case Element::PYRAMID:
+         {
+            for (int j = 0; j < 1; j++)
+            {
+               const int *fv = pyr_t::FaceVert[j];
+               faces_tbl->Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
+            }
+            for (int j = 1; j < 5; j++)
+            {
+               const int *fv = pyr_t::FaceVert[j];
+               faces_tbl->Push(v[fv[0]], v[fv[1]], v[fv[2]]);
+            }
+            break;
+         }
          case Element::WEDGE:
          {
             for (int j = 0; j < 2; j++)
@@ -6308,6 +6488,22 @@ STable3D *Mesh::GetElementToFaceTable(int ret_ftbl)
             }
             break;
          }
+         case Element::PYRAMID:
+         {
+            for (int j = 0; j < 1; j++)
+            {
+               const int *fv = pyr_t::FaceVert[j];
+               el_to_face->Push(
+                  i, faces_tbl->Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]));
+            }
+            for (int j = 1; j < 5; j++)
+            {
+               const int *fv = pyr_t::FaceVert[j];
+               el_to_face->Push(
+                  i, faces_tbl->Push(v[fv[0]], v[fv[1]], v[fv[2]]));
+            }
+            break;
+         }
          case Element::HEXAHEDRON:
          {
             // find the face by the vertices with the smallest 3 numbers
@@ -7644,6 +7840,18 @@ void Mesh::UniformRefinement3D_base(Array<int> *f2qf_ptr, DSTable *v_to_v_p,
       }
    }
 
+   int pyr_counter = 0;
+   if (HasGeometry(Geometry::PYRAMID))
+   {
+      for (int i = 0; i < elements.Size(); i++)
+      {
+         if (elements[i]->GetType() == Element::PYRAMID)
+         {
+            pyr_counter++;
+         }
+      }
+   }
+
    // Map from edge-index to vertex-index, needed for ReorientTetMesh() for
    // parallel meshes.
    // Note: with the removal of ReorientTetMesh() this may no longer
@@ -7703,7 +7911,7 @@ void Mesh::UniformRefinement3D_base(Array<int> *f2qf_ptr, DSTable *v_to_v_p,
    Array<Element*> new_boundary;
 
    vertices.SetSize(oelem + hex_counter);
-   new_elements.SetSize(8 * NumOfElements);
+   new_elements.SetSize(8 * NumOfElements + 2 * pyr_counter);
    CoarseFineTr.embeddings.SetSize(new_elements.Size());
 
    hex_counter = 0;
@@ -7969,6 +8177,73 @@ void Mesh::UniformRefinement3D_base(Array<int> *f2qf_ptr, DSTable *v_to_v_p,
          }
          break;
 
+         case Element::PYRAMID:
+         {
+            const int *f = el_to_face->GetRow(i);
+            // pyr_counter++;
+
+            for (int fi = 0; fi < 1; fi++)
+            {
+               for (int k = 0; k < 4; k++)
+               {
+                  vv[k] = v[pyr_t::FaceVert[fi][k]];
+               }
+               AverageVertices(vv, 4, oface + f2qf[f[fi]]);
+            }
+
+            for (int ei = 0; ei < 8; ei++)
+            {
+               for (int k = 0; k < 2; k++)
+               {
+                  vv[k] = v[pyr_t::Edges[ei][k]];
+               }
+               AverageVertices(vv, 2, oedge+e[ei]);
+            }
+
+            const int qf0 = f2qf[f[0]];
+
+            new_elements[j++] =
+               new Pyramid(v[0], oedge+e[0], oface+qf0,
+                           oedge+e[3], oedge+e[4], attr);
+
+            new_elements[j++] =
+               new Pyramid(oedge+e[0], v[1], oedge+e[1],
+                           oface+qf0, oedge+e[5], attr);
+
+            new_elements[j++] =
+               new Pyramid(oface+qf0, oedge+e[1], v[2],
+                           oedge+e[2], oedge+e[6], attr);
+
+            new_elements[j++] =
+               new Pyramid(oedge+e[3], oface+qf0, oedge+e[2],
+                           v[3], oedge+e[7], attr);
+
+            new_elements[j++] =
+               new Pyramid(oedge+e[4], oedge+e[5], oedge+e[6],
+                           oedge+e[7], v[4], attr);
+
+            new_elements[j++] =
+               new Pyramid(oedge+e[7], oedge+e[6], oedge+e[5],
+                           oedge+e[4], oface+qf0, attr);
+
+            new_elements[j++] =
+               new Tetrahedron(oedge+e[0], oedge+e[4], oedge+e[5],
+                               oface+qf0, attr);
+
+            new_elements[j++] =
+               new Tetrahedron(oedge+e[1], oedge+e[5], oedge+e[6],
+                               oface+qf0, attr);
+
+            new_elements[j++] =
+               new Tetrahedron(oedge+e[2], oedge+e[6], oedge+e[7],
+                               oface+qf0, attr);
+
+            new_elements[j++] =
+               new Tetrahedron(oedge+e[3], oedge+e[7], oedge+e[4],
+                               oface+qf0, attr);
+         }
+         break;
+
          case Element::HEXAHEDRON:
          {
             const int *f = el_to_face->GetRow(i);
@@ -8100,7 +8375,7 @@ void Mesh::UniformRefinement3D_base(Array<int> *f2qf_ptr, DSTable *v_to_v_p,
    }
    mfem::Swap(boundary, new_boundary);
 
-   static const double A = 0.0, B = 0.5, C = 1.0;
+   static const double A = 0.0, B = 0.5, C = 1.0, D = -1.0;
    static double tet_children[3*4*16] =
    {
       A,A,A, B,A,A, A,B,A, A,A,B,
@@ -8126,6 +8401,19 @@ void Mesh::UniformRefinement3D_base(Array<int> *f2qf_ptr, DSTable *v_to_v_p,
       A,A,B, A,B,B, B,A,B, B,B,A,
       A,A,B, B,A,B, B,A,A, B,B,A
    };
+   static double pyr_children[3*5*10] =
+   {
+      A,A,A, B,A,A, B,B,A, A,B,A, A,A,B,
+      B,A,A, C,A,A, C,B,A, B,B,A, B,A,B,
+      B,B,A, C,B,A, C,C,A, B,C,A, B,B,B,
+      A,B,A, B,B,A, B,C,A, A,C,A, A,B,B,
+      A,A,B, B,A,B, B,B,B, A,B,B, A,A,C,
+      A,B,B, B,B,B, B,A,B, A,A,B, B,B,A,
+      B,A,A, A,A,B, B,A,B, B,B,A, D,D,D,
+      C,B,A, B,A,B, B,B,B, B,B,A, D,D,D,
+      B,C,A, B,B,B, A,B,B, B,B,A, D,D,D,
+      A,B,A, A,B,B, A,A,B, B,B,A, D,D,D
+   };
    static double pri_children[3*6*8] =
    {
       A,A,A, B,A,A, A,B,A, A,A,B, B,A,B, A,B,B,
@@ -8151,6 +8439,8 @@ void Mesh::UniformRefinement3D_base(Array<int> *f2qf_ptr, DSTable *v_to_v_p,
 
    CoarseFineTr.point_matrices[Geometry::TETRAHEDRON]
    .UseExternalData(tet_children, 3, 4, 16);
+   CoarseFineTr.point_matrices[Geometry::PYRAMID]
+   .UseExternalData(pyr_children, 3, 5, 10);
    CoarseFineTr.point_matrices[Geometry::PRISM]
    .UseExternalData(pri_children, 3, 6, 8);
    CoarseFineTr.point_matrices[Geometry::CUBE]
@@ -8167,7 +8457,7 @@ void Mesh::UniformRefinement3D_base(Array<int> *f2qf_ptr, DSTable *v_to_v_p,
    }
 
    NumOfVertices    = vertices.Size();
-   NumOfElements    = 8 * NumOfElements;
+   NumOfElements    = 8 * NumOfElements + 2 * pyr_counter;
    NumOfBdrElements = 4 * NumOfBdrElements;
 
    GetElementToFaceTable();
@@ -8751,7 +9041,7 @@ void Mesh::GeneralRefinement(const Array<Refinement> &refinements,
    else if (nonconforming < 0)
    {
       // determine if nonconforming refinement is suitable
-      if ((meshgen & 2) || (meshgen & 4))
+      if ((meshgen & 2) || (meshgen & 4) || (meshgen & 8))
       {
          nonconforming = 1; // tensor product elements and wedges
       }
@@ -9535,6 +9825,7 @@ void Mesh::Printer(std::ostream &out, std::string section_delimiter) const
        "# TETRAHEDRON = 4\n"
        "# CUBE        = 5\n"
        "# PRISM       = 6\n"
+       "# PYRAMID     = 7\n"
        "#\n";
 
    out << "\ndimension\n" << Dim;

mesh/mesh.hpp

@@ -177,6 +177,7 @@ protected:
    int own_nodes;
 
    static const int vtk_quadratic_tet[10];
+   static const int vtk_quadratic_pyramid[13];
    static const int vtk_quadratic_wedge[18];
    static const int vtk_quadratic_hex[27];
 
@@ -195,6 +196,7 @@ public:
    typedef Geometry::Constants<Geometry::TETRAHEDRON> tet_t;
    typedef Geometry::Constants<Geometry::CUBE>        hex_t;
    typedef Geometry::Constants<Geometry::PRISM>       pri_t;
+   typedef Geometry::Constants<Geometry::PYRAMID>     pyr_t;
 
    enum Operation { NONE, REFINE, DEREFINE, REBALANCE };
 
@@ -373,11 +375,17 @@ protected:
    void GetLocalTriToWdgTransformation (IsoparametricTransformation &loc,
                                         int i);
    /// Used in GetFaceElementTransformations (...)
+   void GetLocalTriToPyrTransformation (IsoparametricTransformation &loc,
+                                        int i);
+   /// Used in GetFaceElementTransformations (...)
    void GetLocalQuadToHexTransformation (IsoparametricTransformation &loc,
                                          int i);
    /// Used in GetFaceElementTransformations (...)
    void GetLocalQuadToWdgTransformation (IsoparametricTransformation &loc,
                                          int i);
+   /// Used in GetFaceElementTransformations (...)
+   void GetLocalQuadToPyrTransformation (IsoparametricTransformation &loc,
+                                         int i);
 
    /** Used in GetFaceElementTransformations to account for the fact that a
        slave face occupies only a portion of its master face. */
@@ -657,11 +665,15 @@ public:
    int AddWedge(int v1, int v2, int v3, int v4, int v5, int v6, int attr = 1);
    int AddWedge(const int *vi, int attr = 1);
 
+   int AddPyramid(int v1, int v2, int v3, int v4, int v5, int attr = 1);
+   int AddPyramid(const int *vi, int attr = 1);
+
    int AddHex(int v1, int v2, int v3, int v4, int v5, int v6, int v7, int v8,
               int attr = 1);
    int AddHex(const int *vi, int attr = 1);
    void AddHexAsTets(const int *vi, int attr = 1);
    void AddHexAsWedges(const int *vi, int attr = 1);
+   void AddHexAsPyramids(const int *vi, int attr = 1);
 
    /// The parameter @a elem should be allocated using the NewElement() method
    int AddElement(Element *elem);
@@ -829,10 +841,16 @@ public:
 
    /** @brief Get the mesh generator/type.
 
+       The purpose of this is to be able to quickly tell what type of elements
+       one has in the mesh. Examination of this bitmask along with knowledge
+       of the mesh dimension can be used to identify which element types are
+       present.
+
        @return A bitmask:
        - bit 0 - simplices are present in the mesh (triangles, tets),
        - bit 1 - tensor product elements are present in the mesh (quads, hexes),
        - bit 2 - the mesh has wedge elements.
+       - bit 3 - the mesh has pyramid elements.
 
        In parallel, the result takes into account elements on all processors.
    */

mesh/mesh_headers.hpp

@@ -27,6 +27,7 @@
 #include "mesh_operators.hpp"
 #include "nurbs.hpp"
 #include "wedge.hpp"
+#include "pyramid.hpp"
 
 #ifdef MFEM_USE_MESQUITE
 #include "mesquite.hpp"

mesh/mesh_readers.cpp

@@ -352,6 +352,11 @@ void Mesh::ReadTrueGridMesh(std::istream &input)
 const int Mesh::vtk_quadratic_tet[10] =
 { 0, 1, 2, 3, 4, 7, 5, 6, 8, 9 };
 
+// see Pyramid::edges & Mesh::GenerateFaces
+// https://www.vtk.org/doc/nightly/html/classvtkBiQuadraticQuadraticWedge.html
+const int Mesh::vtk_quadratic_pyramid[13] =
+{ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12};
+
 // see Wedge::edges & Mesh::GenerateFaces
 // https://www.vtk.org/doc/nightly/html/classvtkBiQuadraticQuadraticWedge.html
 const int Mesh::vtk_quadratic_wedge[18] =
@@ -545,6 +550,8 @@ void Mesh::CreateVTKMesh(const Vector &points, const Array<int> &cell_data,
                   vtk_mfem = vtk_quadratic_hex; break;
                case Geometry::PRISM:
                   vtk_mfem = vtk_quadratic_wedge; break;
+               case Geometry::PYRAMID:
+                  vtk_mfem = vtk_quadratic_pyramid; break;
                default:
                   vtk_mfem = NULL; // suppress a warning
                   break;
@@ -1393,6 +1400,10 @@ void Mesh::ReadInlineMesh(std::istream &input, bool generate_edges)
          {
             type = Element::WEDGE;
          }
+         else if (eltype == "pyramid")
+         {
+            type = Element::PYRAMID;
+         }
          else if (eltype == "tet")
          {
             type = Element::TETRAHEDRON;
@@ -1447,7 +1458,7 @@ void Mesh::ReadInlineMesh(std::istream &input, bool generate_edges)
       Make2D(nx, ny, type, sx, sy, generate_edges, true);
    }
    else if (type == Element::TETRAHEDRON || type == Element::WEDGE ||
-            type == Element::HEXAHEDRON)
+            type == Element::HEXAHEDRON  || type == Element::PYRAMID)
    {
       MFEM_VERIFY(nx > 0 && ny > 0 && nz > 0 &&
                   sx > 0.0 && sy > 0.0 && sz > 0.0,

mesh/ncmesh.cpp

@@ -612,7 +612,7 @@ int NCMesh::NewSegment(int n0, int n1, int attr, int vattr1, int vattr2)
    // get (degenerate) faces and assign face attributes
    int v0 = el.node[0], v1 = el.node[1];
    faces.Get(v0, v0, v0, v0)->attribute = vattr1;
-   faces.Get(v1, v1, v1 ,v1)->attribute = vattr2;
+   faces.Get(v1, v1, v1,v1)->attribute = vattr2;
 
    return new_id;
 }

mesh/pyramid.cpp

@@ -0,0 +1,64 @@
+// 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.
+
+// Implementation of class Pyramid
+
+#include "mesh_headers.hpp"
+
+namespace mfem
+{
+
+Pyramid::Pyramid(const int *ind, int attr)
+   : Element(Geometry::PYRAMID)
+{
+   attribute = attr;
+   for (int i = 0; i < 5; i++)
+   {
+      indices[i] = ind[i];
+   }
+}
+
+Pyramid::Pyramid(int ind1, int ind2, int ind3, int ind4, int ind5, int attr)
+   : Element(Geometry::PYRAMID)
+{
+   attribute  = attr;
+   indices[0] = ind1;
+   indices[1] = ind2;
+   indices[2] = ind3;
+   indices[3] = ind4;
+   indices[4] = ind5;
+}
+
+void Pyramid::SetVertices(const int *ind)
+{
+   for (int i = 0; i < 5; i++)
+   {
+      indices[i] = ind[i];
+   }
+}
+
+void Pyramid::GetVertices(Array<int> &v) const
+{
+   v.SetSize(5);
+   for (int i = 0; i < 5; i++)
+   {
+      v[i] = indices[i];
+   }
+}
+
+int Pyramid::GetNFaces(int &nFaceVertices) const
+{
+   MFEM_ABORT("this method is not valid for Pyramid elements");
+   nFaceVertices = 4;
+   return 5;
+}
+
+}

mesh/pyramid.hpp

@@ -0,0 +1,78 @@
+// 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.
+
+#ifndef MFEM_PYRAMID
+#define MFEM_PYRAMID
+
+#include "../config/config.hpp"
+#include "element.hpp"
+
+namespace mfem
+{
+
+/// Data type Pyramid element
+class Pyramid : public Element
+{
+protected:
+   int indices[5];
+
+public:
+   typedef Geometry::Constants<Geometry::PYRAMID> geom_t;
+
+   Pyramid() : Element(Geometry::PYRAMID) { }
+
+   /// Constructs pyramid by specifying the indices and the attribute.
+   Pyramid(const int *ind, int attr = 1);
+
+   /// Constructs pyramid by specifying the indices and the attribute.
+   Pyramid(int ind1, int ind2, int ind3, int ind4, int ind5,
+           int attr = 1);
+
+   /// Return element's type.
+   virtual Type GetType() const { return Element::PYRAMID; }
+
+   /// Set the vertices according to the given input.
+   virtual void SetVertices(const int *ind);
+
+   /// Returns the indices of the element's  vertices.
+   virtual void GetVertices(Array<int> &v) const;
+
+   virtual int *GetVertices() { return indices; }
+
+   virtual int GetNVertices() const { return 5; }
+
+   virtual int GetNEdges() const { return 8; }
+
+   virtual const int *GetEdgeVertices(int ei) const
+   { return geom_t::Edges[ei]; }
+
+   /// @deprecated Use GetNFaces(void) and GetNFaceVertices(int) instead.
+   MFEM_DEPRECATED virtual int GetNFaces(int &nFaceVertices) const;
+
+   virtual int GetNFaces() const { return 5; }
+
+   virtual int GetNFaceVertices(int fi) const
+   { return ( ( fi < 1 ) ? 4 : 3); }
+
+   virtual const int *GetFaceVertices(int fi) const
+   { return geom_t::FaceVert[fi]; }
+
+   virtual Element *Duplicate(Mesh *m) const
+   { return new Pyramid(indices, attribute); }
+
+   virtual ~Pyramid() { }
+};
+
+extern class LinearPyramidFiniteElement PyramidFE;
+
+}
+
+#endif

mesh/wedge.hpp

@@ -71,8 +71,7 @@ public:
    virtual ~Wedge() { }
 };
 
-// Defined in fe.cpp to ensure construction after 'mfem::poly1d'.
-extern class H1_WedgeElement WedgeFE;
+extern class LinearWedgeFiniteElement WedgeFE;
 
 }
 

miniapps/meshing/minimal-surface.cpp

@@ -1159,8 +1159,8 @@ static double u0(const Vector &x) { return sin(3.0 * PI * (x[1] + x[0])); }
 
 enum {NORM, AREA};
 
-static double qf(const int order, const int ker, Mesh &m,
-                 FiniteElementSpace &fes, GridFunction &u)
+double qf(const int order, const int ker, Mesh &m,
+          FiniteElementSpace &fes, GridFunction &u)
 {
    const Geometry::Type type = m.GetElementBaseGeometry(0);
    const IntegrationRule &ir(IntRules.Get(type, order));

miniapps/shifted/dist_solver.cpp

@@ -483,11 +483,11 @@ void ScreenedPoisson::AssembleElementVector(const FiniteElement &el,
       pval=shapef*elfun;
       if (fval>0.0)
       {
-         elvect.Add( -w , shapef);
+         elvect.Add( -w, shapef);
       }
       else if (fval<0.0)
       {
-         elvect.Add(  w , shapef);
+         elvect.Add(  w, shapef);
       }
    }
 }
@@ -523,7 +523,7 @@ void ScreenedPoisson::AssembleElementGrad(const FiniteElement &el,
       el.CalcPhysDShape(trans, B);
       el.CalcPhysShape(trans,shapef);
 
-      AddMult_a_VVt(w , shapef, elmat);
+      AddMult_a_VVt(w, shapef, elmat);
       AddMult_a_AAt(w * diffcoef, B, elmat);
    }
 }
@@ -683,11 +683,11 @@ void PUMPLaplacian::AssembleElementVector(const FiniteElement &el,
       // add the external load -1 if tval > 0.0; 1 if tval < 0.0;
       if (tval>0.0)
       {
-         elvect.Add( -w*fval , shapef);
+         elvect.Add( -w*fval, shapef);
       }
       else  if (tval<0.0)
       {
-         elvect.Add(  w*fval , shapef);
+         elvect.Add(  w*fval, shapef);
       }
    }
 }

tests/convergence/prates.cpp

@@ -117,7 +117,7 @@ int main(int argc, char *argv[])
       MPI_Finalize();
       return 1;
    }
-   if (prob >3 || prob <0) prob = 0; // default problem = H1
+   if (prob >3 || prob <0) { prob = 0; } // default problem = H1
    if (prob == 3)
    {
       if (kappa < 0)
@@ -221,7 +221,7 @@ int main(int argc, char *argv[])
          gradu = new VectorFunctionCoefficient(dim,gradu_exact);
          b.AddDomainIntegrator(new DomainLFIntegrator(*f));
          b.AddBdrFaceIntegrator(
-      new DGDirichletLFIntegrator(*scalar_u, one, sigma, kappa));
+            new DGDirichletLFIntegrator(*scalar_u, one, sigma, kappa));
          a.AddDomainIntegrator(new DiffusionIntegrator(one));
          a.AddInteriorFaceIntegrator(new DGDiffusionIntegrator(one, sigma, kappa));
          a.AddBdrFaceIntegrator(new DGDiffusionIntegrator(one, sigma, kappa));
@@ -291,8 +291,8 @@ int main(int argc, char *argv[])
 
       x = *X;
       JumpScaling js(1.0, jump_scaling_type == 2 ? JumpScaling::P_SQUARED_OVER_H
-                        : jump_scaling_type == 1 ? JumpScaling::ONE_OVER_H
-                        : JumpScaling::CONSTANT);
+                     : jump_scaling_type == 1 ? JumpScaling::ONE_OVER_H
+                     : JumpScaling::CONSTANT);
       switch (prob)
       {
          case 0: rates.AddH1GridFunction(&x,scalar_u,gradu); break;
@@ -305,7 +305,7 @@ int main(int argc, char *argv[])
       delete B;
       delete A;
 
-      if (l==pr) break;
+      if (l==pr) { break; }
 
       pmesh->UniformRefinement();
       fespace->Update();

tests/convergence/rates.cpp

@@ -104,7 +104,7 @@ int main(int argc, char *argv[])
       args.PrintUsage(cout);
       return 1;
    }
-   if (prob >3 || prob <0) prob = 0; // default problem = H1
+   if (prob >3 || prob <0) { prob = 0; } // default problem = H1
    if (prob == 3)
    {
       if (kappa < 0)
@@ -194,7 +194,7 @@ int main(int argc, char *argv[])
          gradu = new VectorFunctionCoefficient(dim,gradu_exact);
          b.AddDomainIntegrator(new DomainLFIntegrator(*f));
          b.AddBdrFaceIntegrator(
-      new DGDirichletLFIntegrator(*scalar_u, one, sigma, kappa));
+            new DGDirichletLFIntegrator(*scalar_u, one, sigma, kappa));
          a.AddDomainIntegrator(new DiffusionIntegrator(one));
          a.AddInteriorFaceIntegrator(new DGDiffusionIntegrator(one, sigma, kappa));
          a.AddBdrFaceIntegrator(new DGDiffusionIntegrator(one, sigma, kappa));
@@ -224,8 +224,8 @@ int main(int argc, char *argv[])
       }
 
       JumpScaling js(1.0, jump_scaling_type == 2 ? JumpScaling::P_SQUARED_OVER_H
-                        : jump_scaling_type == 1 ? JumpScaling::ONE_OVER_H
-                        : JumpScaling::CONSTANT);
+                     : jump_scaling_type == 1 ? JumpScaling::ONE_OVER_H
+                     : JumpScaling::CONSTANT);
 
       switch (prob)
       {
@@ -235,7 +235,7 @@ int main(int argc, char *argv[])
          case 3: rates.AddL2GridFunction(&x,scalar_u,gradu,&one,js); break;
       }
 
-      if (l==sr) break;
+      if (l==sr) { break; }
 
       mesh->UniformRefinement();
       fespace->Update();

tests/unit/fem/test_doftrans.cpp

@@ -98,7 +98,7 @@ TEST_CASE("DoF Transformation Classes",
          double uAv = A.InnerProduct(v, u);
 
          REQUIRE(fabs(uAv -  At.InnerProduct(vt, u )) < tol * fabs(uAv));
-         REQUIRE(fabs(uAv -  tA.InnerProduct(v , ut)) < tol * fabs(uAv));
+         REQUIRE(fabs(uAv -  tA.InnerProduct(v, ut)) < tol * fabs(uAv));
          REQUIRE(fabs(uAv - tAt.InnerProduct(vt, ut)) < tol * fabs(uAv));
       }
       SECTION("Inner product of a primal vector and a dual vector")
@@ -119,7 +119,7 @@ TEST_CASE("DoF Transformation Classes",
          double fAv = A.InnerProduct(v, f);
 
          REQUIRE(fabs(fAv -  At.InnerProduct(vt, f )) < tol * fabs(fAv));
-         REQUIRE(fabs(fAv -  tA.InnerProduct(v , ft)) < tol * fabs(fAv));
+         REQUIRE(fabs(fAv -  tA.InnerProduct(v, ft)) < tol * fabs(fAv));
          REQUIRE(fabs(fAv - tAt.InnerProduct(vt, ft)) < tol * fabs(fAv));
       }
    }
@@ -185,9 +185,9 @@ TEST_CASE("DoF Transformation Functions",
 
       double fAv = A.InnerProduct(v, f);
 
-      REQUIRE(fabs(fAv - nAn.InnerProduct(v , f )) < tol * fabs(fAv));
+      REQUIRE(fabs(fAv - nAn.InnerProduct(v, f )) < tol * fabs(fAv));
       REQUIRE(fabs(fAv -  At.InnerProduct(vt, f )) < tol * fabs(fAv));
-      REQUIRE(fabs(fAv -  tA.InnerProduct(v , ft)) < tol * fabs(fAv));
+      REQUIRE(fabs(fAv -  tA.InnerProduct(v, ft)) < tol * fabs(fAv));
       REQUIRE(fabs(fAv - tAt.InnerProduct(vt, ft)) < tol * fabs(fAv));
    }
    SECTION("TransformDual")
@@ -217,9 +217,9 @@ TEST_CASE("DoF Transformation Functions",
 
       double uAv = A.InnerProduct(v, u);
 
-      REQUIRE(fabs(uAv - nAn.InnerProduct(v , u )) < tol * fabs(uAv));
+      REQUIRE(fabs(uAv - nAn.InnerProduct(v, u )) < tol * fabs(uAv));
       REQUIRE(fabs(uAv -  At.InnerProduct(vt, u )) < tol * fabs(uAv));
-      REQUIRE(fabs(uAv -  tA.InnerProduct(v , ut)) < tol * fabs(uAv));
+      REQUIRE(fabs(uAv -  tA.InnerProduct(v, ut)) < tol * fabs(uAv));
       REQUIRE(fabs(uAv - tAt.InnerProduct(vt, ut)) < tol * fabs(uAv));
    }
 }