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)

.gitignore

@@ -52,13 +52,6 @@ examples/ex1[0-9]p
 examples/ex2[0-9]
 examples/ex2[0-9]p
 
-examples/soln-*
-examples/reference-*
-examples/*.dat
-examples/*.gif
-examples/*.mp4
-examples/*.gf
-
 examples/refined.mesh
 examples/displaced.mesh
 examples/mesh.*

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/combine-figures

@@ -1,30 +0,0 @@
-#!/bin/bash
-
-# Use imagemagick to combine the figures into a composite view
-
-for i in {0..19}; do
-    ii=`printf %04d $i`
-
-    # label reference density
-    convert reference-$ii.png -gravity south -annotate -50+50 "Reference Fine" ,reference-$ii.png
-
-    # label optimal density
-    convert opt-soln-$ii.png -gravity south -annotate -50+50 "Optimal AMR" ,opt-soln-$ii.png
-
-    # label greedy density
-    convert greedy-soln-$ii.png -gravity south -annotate -50+50 "Max Current AMR" ,greedy-soln-$ii.png
-
-    # stack optimal on top of error plot
-    convert ,opt-soln-$ii.png opt-err-$i.png -append ,opt-soln-werr-$ii.png
-
-    # stack greedy on top of error plot
-    convert ,greedy-soln-$ii.png greedy-err-$i.png -append ,greedy-soln-werr-$ii.png
-    
-    # lay the 3 side by side
-    convert ,reference-$ii.png ,opt-soln-werr-$ii.png ,greedy-soln-werr-$ii.png +append ,ref-sbs-soln-$ii.1.png
-
-    # Label the composite image with timestep
-    convert ,ref-sbs-soln-$ii.1.png -gravity northwest -annotate +20+20 "Timestep $i" ref-sbs-soln-$ii.png
-done
-
-rm ,*

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/ex18.cpp

@@ -42,7 +42,6 @@
 #include <fstream>
 #include <sstream>
 #include <iostream>
-#include <set>
 
 // Classes FE_Evolution, RiemannSolver, DomainIntegrator and FaceIntegrator
 // shared between the serial and parallel version of the example.
@@ -51,177 +50,14 @@
 // Choice for the problem setup. See InitialCondition in ex18.hpp.
 int problem;
 
-// The dim and hence number of equations is taken from the mesh.
-int dim;
-int num_equation;
-
 // Equation constant parameters.
+const int num_equation = 4;
 const double specific_heat_ratio = 1.4;
 const double gas_constant = 1.0;
 
 // Maximum characteristic speed (updated by integrators)
 double max_char_speed;
 
-void amr_update(FiniteElementSpace& fes,
-                FiniteElementSpace& dfes,
-                FiniteElementSpace& vfes,
-                BlockVector& u_block,
-                GridFunction& rho,
-                GridFunction& rho_u,
-                GridFunction& rho_e,
-                GridFunction& sol)
-{
-   fes.Update();
-   dfes.Update();
-   vfes.Update();
-
-   rho.Update();
-   rho_u.Update();
-   rho_e.Update();
-
-   Array<int> offsets(num_equation + 1);
-   for (int k = 0; k <= num_equation; k++) {
-      offsets[k] = k * vfes.GetNDofs();
-   }
-   u_block.Update(offsets);
-
-   Vector& sub0 = u_block.GetBlock(0);
-   Vector& sub1 = u_block.GetBlock(1);
-   Vector& sub2 = u_block.GetBlock(2);
-
-   sub0 = rho;
-   sub1 = rho_u;
-   sub2 = rho_e;
-
-   sol.MakeRef(&vfes, u_block.GetData());
-
-   rho.MakeRef(&fes, u_block.GetData() + offsets[0]);
-   rho_u.MakeRef(&dfes, u_block.GetData() + offsets[1]);
-   rho_e.MakeRef(&fes, u_block.GetData() + offsets[2]);
-}
-
-double estimate_dt(double cfl, Mesh& mesh,
-                   GridFunction& sol, NonlinearForm& A, int order)
-{
-   // Determine the minimum element size.
-   double hmin = 0.0;
-   hmin = mesh.GetElementSize(0, 1);
-   for (int i = 1; i < mesh.GetNE(); i++)
-   {
-      hmin = min(mesh.GetElementSize(i, 1), hmin);
-   }
-
-   double dt = -0.01;
-   // Find a safe dt, using a temporary vector. Calling Mult() computes the
-   // maximum char speed at all quadrature points on all faces.
-   Vector z(A.Width());
-   max_char_speed = 0.;
-   A.Mult(sol, z);
-   dt = cfl * hmin / max_char_speed / (2*order+1);
-
-   return dt;
-}
-
-void show_maps(Mesh& mesh, vector<int>& coarse_map, vector<int>& fine_map)
-{
-   for (int i = 0; i < mesh.GetNE(); i++) {
-      int depth = mesh.ncmesh->GetElementDepth(i);
-      if (depth == 0) {
-         printf("| %2d  ",i);
-      }
-      else {
-         printf("|%2d|%2d",i,i+1);
-         i++;
-      }
-   }
-   printf("| mesh i\n");
-
-   for (size_t i = 0; i < coarse_map.size(); i++) {
-      printf("| %2d  ",coarse_map[i]);
-   }
-   printf("| coarse map\n");
-
-   for (size_t i = 0; i < coarse_map.size(); i++) {
-      printf("| %2d  ",fine_map[i]);
-   }
-   printf("| fine map\n");
-
-   for (size_t i = 0; i < coarse_map.size(); i++) {
-      printf("| %2lu  ",i);
-   }
-   printf("| ref idx\n");
-
-   // some validation checks
-   for (size_t i = 0; i < coarse_map.size(); i++) {
-      //printf("checking reference index %d\n",i);
-      if (coarse_map[i] > -1) assert(  fine_map[i] == -1);
-      if (fine_map[i]   > -1) assert(coarse_map[i] == -1);
-      if (coarse_map[i] > -1) assert(mesh.ncmesh->GetElementDepth(coarse_map[i]) == 0);
-      if (fine_map[i]   > -1) assert(mesh.ncmesh->GetElementDepth(fine_map[i])   == 1);
-      if (fine_map[i]   > -1) assert(mesh.ncmesh->GetElementDepth(fine_map[i]+1) == 1);
-   }
-}
-
-void compute_reference_errors(int ti,
-                              FiniteElementCollection& fec,
-                              Mesh& mesh,
-                              GridFunction& den,
-                              const vector<int>& coarse_map,
-                              const vector<int>& fine_map,
-                              vector<double>& errors)
-{
-   // Create a temp copy of the mesh for refinement to reference resolution
-   Mesh tmp_mesh(mesh);
-   FiniteElementSpace tmp_fes(&tmp_mesh, &fec);
-   GridFunction tmp_den(&tmp_fes);
-   tmp_den = den; // copy data
-
-   // refine mesh and soln to reference resolution
-   Array<int> refs;
-   for (int i = 0; i < tmp_mesh.GetNE(); i++) {
-      if (tmp_mesh.ncmesh->GetElementDepth(i) == 0) {
-         refs.Append(i);
-      }
-   }
-   tmp_mesh.GeneralRefinement(refs);
-   tmp_den.FESpace()->Update();
-   tmp_den.Update();
-
-   // read in reference mesh
-   ostringstream fn;
-   fn << "reference-" << ti << ".mesh";
-   Mesh ref_mesh(fn.str().c_str());
-
-   // read in reference density
-   ostringstream sol_name;
-   sol_name << "reference-" << ti << "-" << 0 << ".gf";
-   ifstream gf_ifs(sol_name.str());
-   GridFunction ref_den(&ref_mesh, gf_ifs);
-
-   // integrate L2 errors
-   GridFunctionCoefficient ref_den_coeff(&ref_den);
-   Vector err(ref_mesh.GetNE());
-   tmp_den.ComputeElementL2Errors(ref_den_coeff, err);
-
-   {
-      ostringstream sol_name;
-      sol_name << "refined-current-rho-" << ti << ".gf";
-      ofstream sol_ofs(sol_name.str().c_str());
-      sol_ofs << tmp_den;
-   }
-
-
-   // sum L2 errors to reference mesh. The reference errors are always
-   // on a fully refined mesh.
-   int k = 0;
-   for (size_t i = 0; i < coarse_map.size(); i++) {
-      double e0 = err[k];
-      double e1 = err[k+1];
-      errors[i] = sqrt(e0*e0+e1*e1);
-      k += 2;
-   }
-}
-
 int main(int argc, char *argv[])
 {
    // 1. Parse command-line options.
@@ -233,15 +69,8 @@ int main(int argc, char *argv[])
    double t_final = 2.0;
    double dt = -0.01;
    double cfl = 0.3;
-   bool visualization = false;
+   bool visualization = true;
    int vis_steps = 50;
-   Array<int> rseq;
-   int nseq = 0;
-   int regrid_period = 1;
-   int greedy_refine = 0;
-   int lookahead_refine = 0;
-   int output_cycle_soln = 0;
-   int output_cycle_errors = 0;
 
    int precision = 8;
    cout.precision(precision);
@@ -269,18 +98,6 @@ int main(int argc, char *argv[])
                   "Enable or disable GLVis visualization.");
    args.AddOption(&vis_steps, "-vs", "--visualization-steps",
                   "Visualize every n-th timestep.");
-   args.AddOption(&rseq, "-rseq", "--refine-sequence",
-                  "Element sequence to refine.");
-   args.AddOption(&regrid_period, "-rp", "--regrid-period",
-                  "Timesteps per regrid.");
-   args.AddOption(&greedy_refine, "-gr", "--greedy-refine",
-                  "Use greedy refinement strategy.");
-   args.AddOption(&lookahead_refine, "-lr", "--lookahead-refine",
-                  "Use lookahead refinement strategy.");
-   args.AddOption(&output_cycle_soln, "-cs", "--output-cycle-solutions",
-                  "Output per-cycle solutions.");
-   args.AddOption(&output_cycle_errors, "-ce", "--output-cycle-errors",
-                  "Output per-cycle errors.");
 
    args.Parse();
    if (!args.Good())
@@ -290,15 +107,12 @@ int main(int argc, char *argv[])
    }
    args.PrintOptions(cout);
 
-   printf("dt = %f\n",dt);
-
    // 2. Read the mesh from the given mesh file. This example requires a 2D
    //    periodic mesh, such as ../data/periodic-square.mesh.
    Mesh mesh(mesh_file, 1, 1);
-   mesh.EnsureNCMesh();
    const int dim = mesh.Dimension();
 
-   num_equation = 2+dim;
+   MFEM_ASSERT(dim == 2, "Need a two-dimensional mesh for the problem definition");
 
    // 3. Define the ODE solver used for time integration. Several explicit
    //    Runge-Kutta methods are available.
@@ -348,9 +162,7 @@ int main(int argc, char *argv[])
    BlockVector u_block(offsets);
 
    // Momentum grid function on dfes for visualization.
-   GridFunction rho(&fes, u_block.GetData() + offsets[0]);
-   GridFunction rho_u(&dfes, u_block.GetData() + offsets[1]);
-   GridFunction rho_e(&fes, u_block.GetData() + offsets[2]);
+   GridFunction mom(&dfes, u_block.GetData() + offsets[1]);
 
    // Initialize the state.
    VectorFunctionCoefficient u0(num_equation, InitialCondition);
@@ -372,7 +184,6 @@ int main(int argc, char *argv[])
          sol_ofs.precision(precision);
          sol_ofs << uk;
       }
-
    }
 
    // 7. Set up the nonlinear form corresponding to the DG discretization of the
@@ -388,7 +199,7 @@ int main(int argc, char *argv[])
    // 8. Define the time-dependent evolution operator describing the ODE
    //    right-hand side, and perform time-integration (looping over the time
    //    iterations, ti, with a time-step dt).
-   FE_Evolution euler(vfes, A, Aflux);
+   FE_Evolution euler(vfes, A, Aflux.SpMat());
 
    // Visualize the density
    socketstream sout;
@@ -408,7 +219,7 @@ int main(int argc, char *argv[])
       else
       {
          sout.precision(precision);
-         sout << "solution\n" << mesh << rho;
+         sout << "solution\n" << mesh << mom;
          sout << "pause\n";
          sout << flush;
          cout << "GLVis visualization paused."
@@ -416,6 +227,17 @@ int main(int argc, char *argv[])
       }
    }
 
+   // Determine the minimum element size.
+   double hmin = 0.0;
+   if (cfl > 0)
+   {
+      hmin = mesh.GetElementSize(0, 1);
+      for (int i = 1; i < mesh.GetNE(); i++)
+      {
+         hmin = min(mesh.GetElementSize(i, 1), hmin);
+      }
+   }
+
    // Start the timer.
    tic_toc.Clear();
    tic_toc.Start();
@@ -424,348 +246,27 @@ int main(int argc, char *argv[])
    euler.SetTime(t);
    ode_solver->Init(euler);
 
-   // map from base coarse numbering to current coarse numbering.
-   // -1 entries for elements that are currently refined.
-   vector<int> coarse_map(mesh.GetNE());
-   for (size_t i = 0; i < coarse_map.size(); ++i) {
-      coarse_map[i] = i;
+   if (cfl > 0)
+   {
+      // Find a safe dt, using a temporary vector. Calling Mult() computes the
+      // maximum char speed at all quadrature points on all faces.
+      Vector z(A.Width());
+      max_char_speed = 0.;
+      A.Mult(sol, z);
+      dt = cfl * hmin / max_char_speed / (2*order+1);
    }
 
-   // map from base coarse numbering to current fine numbering. The
-   // second element is always this index+1.
-   // -1 entries for elements that are currently coarse.
-   vector<int> fine_map(mesh.GetNE());
-   for (size_t i = 0; i < fine_map.size(); ++i) {
-      fine_map[i] = -1;
-   }
-
-   //show_maps(mesh,coarse_map,fine_map);
-
-   // vector<int> coarse_map(mesh.GetNE());
-   // for (int i = 0; i < coarse_map.size(); ++i) {
-   //    coarse_map[i] = i;
-   // }
-
-   // the elements that are currently refined
-   set<int> cur_ref_set; // in base numbering
-
    // Integrate in time.
    bool done = false;
    for (int ti = 0; !done; )
    {
-
-      if (greedy_refine && !(ti % regrid_period )) {
-         vector<double> ref_errors(coarse_map.size());
-         compute_reference_errors(ti, fec, mesh, rho,
-                                  coarse_map, fine_map, ref_errors);
-         vector<double>::iterator it;
-         it = std::max_element(ref_errors.begin(), ref_errors.end());
-         int ref_max = std::distance(ref_errors.begin(), it);
-         std::cout << "max at: " << std::distance(ref_errors.begin(), it) << '\n';
-         rseq.SetSize(nseq+1);
-         rseq[nseq] = ref_max;
-      }
-
-      if (lookahead_refine && !(ti % regrid_period )) {
-
-         MFEM_ASSERT(cfl < 0, "need fixed time step for lookahead.");
-
-         double dt_real = min(dt, t_final - t);
-
-         // advance a copy of the solution
-         GridFunction tmp_sol(sol);
-
-         ode_solver->Step(tmp_sol, t, dt_real);
-
-         GridFunction tmp_rho(&fes, tmp_sol);
-
-         vector<double> ref_errors(coarse_map.size());
-         compute_reference_errors(ti+1, fec, mesh, tmp_rho,
-                                  coarse_map, fine_map, ref_errors);
-
-         vector<double>::iterator it;
-         it = std::max_element(ref_errors.begin(), ref_errors.end());
-         int ref_max = std::distance(ref_errors.begin(), it);
-         std::cout << "max at: " << std::distance(ref_errors.begin(), it) << '\n';
-         rseq.SetSize(nseq+1);
-         rseq[nseq] = ref_max;
-      }
-
-      // adapt mesh to new element(s)
-      if (rseq.Size() && !(ti % regrid_period )) {
-
-         //printf("*** begin refinement ***\n");
-         set<int> new_ref_set; // in base numbering
-
-         // use base numbering
-         int ne = coarse_map.size();
-         int el1 = rseq[nseq++];
-         int el2 = el1-1;
-         int el3 = el1+1;
-         el2 = (el2 + ne) % ne;
-         el3 = (el3 + ne) % ne;
-         new_ref_set.insert(el1);
-         new_ref_set.insert(el2);
-         new_ref_set.insert(el3);
-
-         // The coarsen set is elements from cur_ref_set not in new_ref_set
-         set<int> coarsen_set;
-         std::set_difference(cur_ref_set.begin(), cur_ref_set.end(),
-                             new_ref_set.begin(), new_ref_set.end(),
-                             std::inserter(coarsen_set, coarsen_set.begin()));
-
-         // The refine set is elements from new_ref_set not in cur_ref_set
-         set<int> refine_set;
-         std::set_difference(new_ref_set.begin(), new_ref_set.end(),
-                             cur_ref_set.begin(), cur_ref_set.end(),
-                             std::inserter(refine_set, refine_set.begin()));
-
-         set<int>::iterator it;
-
-         // printf("cur_ref set\n");
-         //
-         // for (it = cur_ref_set.begin(); it != cur_ref_set.end(); ++it) {
-         //    const int& i = *it;
-         //    printf("%d ",i);
-         // }
-         // printf("\n");
-
-         // printf("new_ref set\n");
-         // for (it = new_ref_set.begin(); it != new_ref_set.end(); ++it) {
-         //    const int& i = *it;
-         //    printf("%d ",i);
-         // }
-         // printf("\n");
-
-         // printf("refine set: ");
-         // for (it = refine_set.begin(); it != refine_set.end(); ++it) {
-         //    const int& i = *it;
-         //    printf("%d ",i);
-         // }
-         // printf("\n");
-
-         // printf("coarsen set: ");
-         // for (it = coarsen_set.begin(); it != coarsen_set.end(); ++it) {
-         //    const int& i = *it;
-         //    printf("%d ",i);
-         // }
-         // printf("\n");
-
-
-         // Perform any new refinements. Translate into current numbering.
-         Array<int> els;
-         for (it = refine_set.begin(); it != refine_set.end(); ++it) {
-            const int& i = *it;
-            assert(coarse_map[i] >= 0);
-            //printf("adding %d -> %d to be refined\n",i,coarse_map[i]);
-            els.Append(coarse_map[i]);
-         }
-
-         if (els.Size()) {
-            //printf("  ** start refinements **\n");
-            mesh.GeneralRefinement(els);
-
-            const CoarseFineTransformations &cft = mesh.GetRefinementTransforms();
-            //printf("updating coarse_map and fine_map after refinement\n");
-            Table c2f;
-            cft.GetCoarseToFineMap(mesh, c2f);
-            Array<int> row;
-            for (size_t i = 0; i < coarse_map.size(); ++i) {
-
-               int old = coarse_map[i];
-               if (old > -1) {
-                  // was coarse
-                  c2f.GetRow(old,row);
-                  if (row.Size() == 1) {
-                     // was coarse, still coarse
-                     coarse_map[i] = row[0];
-                  }
-                  else {
-                     // was coarse, now fine
-                     coarse_map[i] = -1;
-                     fine_map[i] = row[0];
-                  }
-               }
-               else {
-                  // was fine, still fine
-                  old = fine_map[i];
-                  assert(old != -1);
-                  c2f.GetRow(old,row);
-                  assert(row.Size() == 1);
-                  fine_map[i] = row[0];
-               }
-            }
-
-            // printf("new coarse_map is:\n");
-            // for (size_t i = 0; i < coarse_map.size(); i++) {
-            //    printf("%lu -> %d\n",i,coarse_map[i]);
-            // }
-            // printf("new fine_map is:\n");
-            // for (size_t i = 0; i < fine_map.size(); i++) {
-            //    printf("%lu -> %d\n",i,fine_map[i]);
-            // }
-
-            amr_update(fes, dfes, vfes, u_block, rho, rho_u, rho_e, sol);
-            A.Update();
-            Aflux.Update();
-            Aflux.Assemble();
-            euler.Update();
-            ode_solver->Init(euler);
-            //printf("  ** done refinements **\n");
-
-            //printf("maps after refinement\n");
-            //show_maps(mesh,coarse_map,fine_map);
-         }
-
-
-         cur_ref_set = new_ref_set;
-
-         // Perform any new derefinements
-         if (coarsen_set.size()) {
-
-            //printf("  ** start derefinements **\n");
-
-            Array<double> mock_error(mesh.GetNE());
-            mock_error = 1.0;
-
-            // We only mark one of the fine elements, but that's fine
-            // because we can set the threshold low enough to always
-            // derefine.
-            for (it = coarsen_set.begin(); it != coarsen_set.end(); ++it) {
-               const int& i = *it;
-               //printf("coarsening ref idx %d\n",i);
-               assert(coarse_map[i] == -1);
-               int i1 = fine_map[i];
-               //printf("which is fine element %d (and +1)\n",i1);
-               assert(i1 >= 0);
-               mock_error[i1] = 0.0;
-               //printf("setting mock error to 0.0 in %d\n",i1);
-            }
-            mesh.DerefineByError(mock_error, 2.0);
-            amr_update(fes, dfes, vfes, u_block, rho, rho_u, rho_e, sol);
-            A.Update();
-            Aflux.Update();
-            Aflux.Assemble();
-            euler.Update();
-            ode_solver->Init(euler);
-
-            //printf("updating coarse_map and fine_map after derefinement\n");
-
-            const CoarseFineTransformations& cft = mesh.ncmesh->GetDerefinementTransforms();
-
-            Table c2f;
-            cft.GetCoarseToFineMap(mesh, c2f);
-            //c2f.Print();
-
-            Array<int> row;
-            map<int,int> old2new;
-            for (int i = 0; i < c2f.Size(); ++i) {
-               c2f.GetRow(i,row);
-               for (int j = 0; j < row.Size(); j++) {
-                  old2new[row[j]] = i;
-               }
-            }
-
-            // map<int,int>::iterator it;
-            // printf("old2new\n");
-            // for (it = old2new.begin(); it != old2new.end(); ++it) {
-            //    printf("%d -> %d\n",it->first,it->second);
-            // }
-
-            for (size_t i = 0; i < coarse_map.size(); i++) {
-               //printf("ref i = %lu\n",i);
-               int jfine = fine_map[i];
-               int jcoarse = coarse_map[i];
-               //printf("mesh j old grid (fine) = %d\n",jfine);
-               if (jfine > -1) {
-                  int newj = old2new[jfine];
-                  //printf("mesh j new grid (coarse) = %d\n",newj);
-                  c2f.GetRow(newj,row);
-                  if (row.Size() > 1) {
-                     //printf("  was fine, now coarse\n");
-                     // was fine, now coarse
-                     coarse_map[i] = newj;
-                     fine_map[i] = -1;
-                  }
-                  else {
-                     //printf("  was fine, still fine\n");
-                     // was fine, still fine
-                     fine_map[i] = newj;
-                  }
-               }
-               //printf("mesh j old grid coarse = %d\n",jcoarse);
-               if (jcoarse > -1) {
-                  int newj = old2new[jcoarse];
-                  coarse_map[i] = newj;
-               }
-            }
-
-            // printf("new coarse_map is:\n");
-            // for (size_t i = 0; i < coarse_map.size(); i++) {
-            //    printf("%lu -> %d\n",i,coarse_map[i]);
-            // }
-            // printf("new fine_map is:\n");
-            // for (size_t i = 0; i < fine_map.size(); i++) {
-            //    printf("%lu -> %d\n",i,fine_map[i]);
-            // }
-
-            //printf("  ** done derefinements **\n");
-
-            //printf("maps after derefinement\n");
-            //show_maps(mesh,coarse_map,fine_map);
-
-         }
-
-      }
-
-      // Output the current solution.
-      if (output_cycle_soln) {
-
-         ostringstream fn;
-         fn << "soln-" << ti << ".mesh";
-         ofstream mesh_ofs(fn.str());
-         mesh_ofs.precision(precision);
-         mesh_ofs << mesh;
-
-         for (int k = 0; k < num_equation; k++)
-         {
-            GridFunction uk(&fes, u_block.GetBlock(k));
-            ostringstream sol_name;
-            sol_name << "soln-" << ti << "-" << k << ".gf";
-            ofstream sol_ofs(sol_name.str().c_str());
-            sol_ofs.precision(precision);
-            sol_ofs << uk;
-         }
-      }
-
-      // Output the current error on the reference mesh
-      if (output_cycle_errors) {
-         vector<double> ref_errors(coarse_map.size());
-         compute_reference_errors(ti, fec, mesh, rho,
-                                  coarse_map, fine_map, ref_errors);
-         ostringstream fn;
-         fn << "err-" << ti << ".dat";
-         ofstream err_ofs(fn.str());
-         vector<double>::iterator it;
-         for (it = ref_errors.begin(); it != ref_errors.end(); ++it) {
-            int i = std::distance(ref_errors.begin(), it);
-            double err = *it;
-            // "mark" refinements with negation. We'll use this to vis
-            // the refined regions on the error plots.
-            if (coarse_map[i] < 0) err *= -1;
-            err_ofs << i << " " << err << endl;
-         }
-      }
-
-
-      if (cfl > 0) {
-         dt = estimate_dt(cfl, mesh, sol, A, order);
-      }
-
       double dt_real = min(dt, t_final - t);
 
       ode_solver->Step(sol, t, dt_real);
-
+      if (cfl > 0)
+      {
+         dt = cfl * hmin / max_char_speed / (2*order+1);
+      }
       ti++;
 
       done = (t >= t_final - 1e-8*dt);
@@ -774,10 +275,9 @@ int main(int argc, char *argv[])
          cout << "time step: " << ti << ", time: " << t << endl;
          if (visualization)
          {
-            sout << "solution\n" << mesh << rho_u << flush;
+            sout << "solution\n" << mesh << mom << flush;
          }
       }
-
    }
 
    tic_toc.Stop();
@@ -785,57 +285,21 @@ int main(int argc, char *argv[])
 
    // 9. Save the final solution. This output can be viewed later using GLVis:
    //    "glvis -m vortex.mesh -g vortex-1-final.gf".
+   for (int k = 0; k < num_equation; k++)
    {
-      ofstream mesh_ofs("vortex-final.mesh");
-      mesh_ofs.precision(precision);
-      mesh_ofs << mesh;
-
-
-      for (int k = 0; k < num_equation; k++)
-      {
-         GridFunction uk(&fes, u_block.GetBlock(k));
-         ostringstream sol_name;
-         sol_name << "vortex-" << k << "-final.gf";
-         ofstream sol_ofs(sol_name.str().c_str());
-         sol_ofs.precision(precision);
-         sol_ofs << uk;
-      }
+      GridFunction uk(&fes, u_block.GetBlock(k));
+      ostringstream sol_name;
+      sol_name << "vortex-" << k << "-final.gf";
+      ofstream sol_ofs(sol_name.str().c_str());
+      sol_ofs.precision(precision);
+      sol_ofs << uk;
    }
 
    // 10. Compute the L2 solution error summed for all components.
-   if (rseq.Size())
+   if (t_final == 2.0)
    {
-      // const CoarseFineTransformations& cft = mesh.GetRefinementTransforms();
-      // Table c2f;
-      // cft.GetCoarseToFineMap(mesh, c2f);
-
-      // Refine to the reference mesh by refining every coarse el that
-      // wasn't refined.
-      Array<int> refs;
-      for (int i = 0; i < mesh.GetNE(); i++) {
-         if (mesh.ncmesh->GetElementDepth(i) == 0) {
-            refs.Append(i);
-         }
-      }
-      mesh.GeneralRefinement(refs);
-
-      amr_update(fes, dfes, vfes, u_block, rho, rho_u, rho_e, sol);
-
-      Mesh mesh_ref("reference.mesh");
-      ifstream rho_ifs("reference-rho.gf");
-      ifstream rho_u_ifs("reference-rho-u.gf");
-      ifstream rho_e_ifs("reference-rho-e.gf");
-      GridFunction rho_ref(&mesh_ref, rho_ifs);
-      GridFunction rho_u_ref(&mesh_ref, rho_u_ifs);
-      GridFunction rho_e_ref(&mesh_ref, rho_e_ifs);
-      rho_ref.Print();
-      rho.Print();
-
-      GridFunctionCoefficient rho_ref_coeff(&rho_ref);
-      const double err = rho.ComputeL2Error(rho_ref_coeff);
-      printf("final L2 error %e\n",err);
-      // const double error = sol.ComputeLpError(2, u0);
-      // cout << "Solution error: " << error << endl;
+      const double error = sol.ComputeLpError(2, u0);
+      cout << "Solution error: " << error << endl;
    }
 
    // Free the used memory.

examples/ex18.hpp

@@ -11,8 +11,7 @@ extern int problem;
 // Maximum characteristic speed (updated by integrators)
 extern double max_char_speed;
 
-extern int num_equation;
-
+extern const int num_equation;
 extern const double specific_heat_ratio;
 extern const double gas_constant;
 
@@ -25,8 +24,7 @@ private:
 
    FiniteElementSpace &vfes;
    Operator &A;
-   MixedBilinearForm& Aflux;
-   //SparseMatrix &Aflux;
+   SparseMatrix &Aflux;
    DenseTensor Me_inv;
 
    mutable Vector state;
@@ -38,13 +36,10 @@ private:
 
 public:
    FE_Evolution(FiniteElementSpace &vfes_,
-                Operator &A_,
-                MixedBilinearForm &Aflux_);
+                Operator &A_, SparseMatrix &Aflux_);
 
    virtual void Mult(const Vector &x, Vector &y) const;
 
-   void Update();
-
    virtual ~FE_Evolution() { }
 };
 
@@ -105,7 +100,7 @@ public:
 
 // Implementation of class FE_Evolution
 FE_Evolution::FE_Evolution(FiniteElementSpace &vfes_,
-                           Operator &A_, MixedBilinearForm &Aflux_)
+                           Operator &A_, SparseMatrix &Aflux_)
    : TimeDependentOperator(A_.Height()),
      dim(vfes_.GetFE(0)->GetDim()),
      vfes(vfes_),
@@ -128,32 +123,6 @@ FE_Evolution::FE_Evolution(FiniteElementSpace &vfes_,
       inv.Factor();
       inv.GetInverseMatrix(Me_inv(i));
    }
-   //Aflux.Assemble();
-}
-
-void FE_Evolution::Update()
-{
-   height = width = A.Height();
-
-   Me_inv.SetSize(vfes.GetFE(0)->GetDof(),
-                  vfes.GetFE(0)->GetDof(),
-                  vfes.GetNE());
-
-   Aflux.Assemble();
-
-   // Standard local assembly and inversion for energy mass matrices.
-   const int dof = vfes.GetFE(0)->GetDof();
-   DenseMatrix Me(dof);
-   DenseMatrixInverse inv(&Me);
-   MassIntegrator mi;
-   for (int i = 0; i < vfes.GetNE(); i++)
-   {
-      mi.AssembleElementMatrix(*vfes.GetFE(i), *vfes.GetElementTransformation(i), Me);
-      inv.Factor();
-      inv.GetInverseMatrix(Me_inv(i));
-   }
-
-   z.SetSize(A.Height());
 }
 
 void FE_Evolution::Mult(const Vector &x, Vector &y) const
@@ -416,8 +385,6 @@ void FaceIntegrator::AssembleFaceVector(const FiniteElement &el1,
                                         FaceElementTransformations &Tr,
                                         const Vector &elfun, Vector &elvect)
 {
-  int dim = el1.GetDim();
-
    // Compute the term <F.n(u),[w]> on the interior faces.
    const int dof1 = el1.GetDof();
    const int dof2 = el2.GetDof();
@@ -457,25 +424,16 @@ void FaceIntegrator::AssembleFaceVector(const FiniteElement &el1,
 
       Tr.SetAllIntPoints(&ip); // set face and element int. points
 
-      const IntegrationPoint &eip1 = Tr.GetElement1IntPoint();
-      const IntegrationPoint &eip2 = Tr.GetElement2IntPoint();
-
       // Calculate basis functions on both elements at the face
-      el1.CalcShape(eip1, shape1);
-      el2.CalcShape(eip2, shape2);
+      el1.CalcShape(Tr.GetElement1IntPoint(), shape1);
+      el2.CalcShape(Tr.GetElement2IntPoint(), shape2);
 
       // Interpolate elfun at the point
       elfun1_mat.MultTranspose(shape1, funval1);
       elfun2_mat.MultTranspose(shape2, funval2);
 
       // Get the normal vector and the flux on the face
-      if (1 == dim) {
-	nor(0) = 2*eip1.x - 1.0;
-      }
-      else {
-	CalcOrtho(Tr.Jacobian(), nor);
-      }
-
+      CalcOrtho(Tr.Jacobian(), nor);
       const double mcs = rsolver.Eval(funval1, funval2, nor, fluxN);
 
       // Update max char speed
@@ -546,147 +504,62 @@ bool StateIsPhysical(const Vector &state, const int dim)
 // Initial condition
 void InitialCondition(const Vector &x, Vector &y)
 {
-  int dim = x.Size();
+   MFEM_ASSERT(x.Size() == 2, "");
 
-  if (dim == 1) {
+   double radius = 0, Minf = 0, beta = 0;
+   if (problem == 1)
+   {
+      // "Fast vortex"
+      radius = 0.2;
+      Minf = 0.5;
+      beta = 1. / 5.;
+   }
+   else if (problem == 2)
+   {
+      // "Slow vortex"
+      radius = 0.2;
+      Minf = 0.05;
+      beta = 1. / 50.;
+   }
+   else
+   {
+      mfem_error("Cannot recognize problem."
+                 "Options are: 1 - fast vortex, 2 - slow vortex");
+   }
 
-     double rho;
-     double u;
-     double e;
+   const double xc = 0.0, yc = 0.0;
 
-     if (problem == 1) {
+   // Nice units
+   const double vel_inf = 1.;
+   const double den_inf = 1.;
 
-        double xc = x[0]-0.5;
-        rho = 1.0;
-        u = 0.0;
-        e = 1.0+0.1*exp(-xc*xc);
+   // Derive remainder of background state from this and Minf
+   const double pres_inf = (den_inf / specific_heat_ratio) * (vel_inf / Minf) *
+                           (vel_inf / Minf);
+   const double temp_inf = pres_inf / (den_inf * gas_constant);
 
-     }
-     else if (problem == 2) {
+   double r2rad = 0.0;
+   r2rad += (x(0) - xc) * (x(0) - xc);
+   r2rad += (x(1) - yc) * (x(1) - yc);
+   r2rad /= (radius * radius);
 
-        u = 0.0;
-        double gamma = 1.4;
-        double gm1 = gamma -1.0;
+   const double shrinv1 = 1.0 / (specific_heat_ratio - 1.);
 
-        double rhoL = 1.0;
-        double pL = 1.0;
+   const double velX = vel_inf * (1 - beta * (x(1) - yc) / radius * exp(
+                                     -0.5 * r2rad));
+   const double velY = vel_inf * beta * (x(0) - xc) / radius * exp(-0.5 * r2rad);
+   const double vel2 = velX * velX + velY * velY;
 
-        double rhoR = 0.125;
-        double pR = 0.1;
+   const double specific_heat = gas_constant * specific_heat_ratio * shrinv1;
+   const double temp = temp_inf - 0.5 * (vel_inf * beta) *
+                       (vel_inf * beta) / specific_heat * exp(-r2rad);
 
-        if (x[0] < 0.5) {
-           rho = rhoL;
-           e = 1./gm1*pL/rhoL;
-        }
-        else {
-           rho = rhoR;
-           e = 1./gm1*pR/rhoR;
-        }
-     }
-     else if (problem == 3) {
-
-        u = 0.0;
-        double gamma = 1.4;
-        double gm1 = gamma -1.0;
-
-        double rhoL = 1.0;
-        double pL = 1.0;
-
-        double rhoR = 0.125;
-        double pR = 0.1;
-
-        rho = rhoR;
-        e = 1./gm1*pR/rhoR;
-        if (x[0] > 0.3 && x[0] < 0.4) {
-           rho = rhoL;
-           e = 1./gm1*pL/rhoL;
-        }
-        if (x[0] > 0.5 && x[0] < 1.0) {
-           double xc = 0.75;
-           double x0 = x[0]-xc;
-           double a = 10.0;
-           double x1 = a*x0;
-           rho = rhoR +0.5*rhoL*exp(-x1*x1);
-           e = 1./gm1*pR/rhoL;
-        }
-     }
-     else {
-        printf("unknown problem!\n");
-        exit(1);
-     }
-
-    // if (x(0) < 0.5) {
-    //   rho = 1.0;
-    //   e = 1.05*(;
-    // }
-
-    int neq = 0;
-    y(neq++) = rho;
-    y(neq++) = rho * u;
-    y(neq++) = rho * e;
-  }
-  else {
-
-    MFEM_ASSERT(x.Size() == 2, "");
-
-    double radius = 0, Minf = 0, beta = 0;
-    if (problem == 1)
-      {
-	// "Fast vortex"
-	radius = 0.2;
-	Minf = 0.5;
-	beta = 1. / 5.;
-      }
-    else if (problem == 2)
-      {
-	// "Slow vortex"
-	radius = 0.2;
-	Minf = 0.05;
-	beta = 1. / 50.;
-      }
-    else
-      {
-	mfem_error("Cannot recognize problem."
-		   "Options are: 1 - fast vortex, 2 - slow vortex");
-      }
-
-    const double xc = 0.0, yc = 0.0;
-
-    // Nice units
-    const double vel_inf = 1.;
-    const double den_inf = 1.;
-
-    // Derive remainder of background state from this and Minf
-    const double pres_inf = (den_inf / specific_heat_ratio) * (vel_inf / Minf) *
-      (vel_inf / Minf);
-    const double temp_inf = pres_inf / (den_inf * gas_constant);
-
-    double r2rad = 0.0;
-    r2rad += (x(0) - xc) * (x(0) - xc);
-    r2rad += (x(1) - yc) * (x(1) - yc);
-    r2rad /= (radius * radius);
-
-    const double shrinv1 = 1.0 / (specific_heat_ratio - 1.);
-
-    const double velX = vel_inf * (1 - beta * (x(1) - yc) / radius * exp(
-									 -0.5 * r2rad));
-    const double velY = vel_inf * beta * (x(0) - xc) / radius * exp(-0.5 * r2rad);
-    const double vel2 = velX * velX + velY * velY;
-
-    const double specific_heat = gas_constant * specific_heat_ratio * shrinv1;
-    const double temp = temp_inf - 0.5 * (vel_inf * beta) *
-      (vel_inf * beta) / specific_heat * exp(-r2rad);
-
-    const double den = den_inf * pow(temp/temp_inf, shrinv1);
-    const double pres = den * gas_constant * temp;
-    const double energy = shrinv1 * pres / den + 0.5 * vel2;
-
-    int neq = 0;
-    y(neq++) = den;
-    y(neq++) = den * velX;
-    y(neq++) = den * velY;
-    y(neq++) = den * energy;
-
-  } // 2d
+   const double den = den_inf * pow(temp/temp_inf, shrinv1);
+   const double pres = den * gas_constant * temp;
+   const double energy = shrinv1 * pres / den + 0.5 * vel2;
 
+   y(0) = den;
+   y(1) = den * velX;
+   y(2) = den * velY;
+   y(3) = den * energy;
 }

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:

examples/generate-glvis-script

@@ -1,48 +0,0 @@
-#!/bin/bash
-
-if [ $# -ne 2 ]; then
-    echo "Usage: generate-glvis-script <ncyc> <output-file-prefix>"
-fi
-
-ncyc=$1
-prefix=$2
-comp=0
-
-ncycm1=$(( $ncyc - 1 ))
-
-cat <<EOF > glvis-script
-window 0 0 400 400
-
-# Initial solution
-solution soln-$comp.mesh soln-0-$comp.gf
-
-# Setup the GLVis scene. Executed after pressing the space bar.
-{
-   perspective off
-   viewcenter 0.2 0
-   valuerange 0.0 1.0
-   keys ca
-   autoscale off
-}
-
-# Take multiple screenshots. Executed after pressing space bar.
-{
-EOF
-
-for i in `seq 0 $ncycm1`; do
-    if [ ! -f soln-$i.mesh ]; then
-	echo "The file soln-$i.mesh doesn't exist. Stopping."
-	exit 1
-    fi
-    if [ ! -f soln-$i-$comp.gf ]; then
-	echo "The file soln-$i-$comp.gf doesn't exist. Stopping."
-	exit 1
-    fi
-    ii=`printf %04d $i`
-    echo "solution soln-$i.mesh soln-$i-$comp.gf" >> glvis-script
-    #echo "valuerange 0.0 1.0" >> glvis-script
-    echo "screenshot $prefix-soln-$ii.png" >> glvis-script
-done
-echo "}" >> glvis-script
-
-echo "% glvis -run glvis-script"

examples/generate-screenshot-cmds

@@ -1,11 +0,0 @@
-#!/bin/bash
-
-rm -f tmp
-for i in {0..19}; do
-    ii=`printf %04d $i`
-    echo "solution soln-$i.mesh soln-$i-0.gf" >> tmp
-    echo "valuerange 0.0 1.0" >> tmp
-    echo "screenshot greedy-soln-$ii.png" >> tmp
-done
-
-echo "output in ./tmp"

examples/make-err-plots.py

@@ -1,32 +0,0 @@
-
-# Take err-N.dat files and create error plots for each cycle.
-# Negative values are interpreted as refined regions and plotted in
-# red.
-
-import math
-import matplotlib.pyplot as plt
-
-for i in range(20):
-    print(f"{i}")
-    with open(f"err-{i}.dat") as f:
-        lines = f.readlines()
-        x = [float(line.split()[0]) for line in lines]
-        y = [float(line.split()[1]) for line in lines]
-
-        # errors in refined zones are marked as negative by convention
-        # from ex18. in order to detect "negative zero", we use the
-        # copysign method, which distinguishes negative zero from
-        # zero.
-        r = [(x,-y) for x,y in zip(x,y) if math.copysign(1, y) == -1.0]
-        xr,yr = zip(*r)
-        print(xr)
-        print(yr)
-        plt.scatter(x,y)
-        plt.scatter(xr,yr,c='r')
-        plt.title(f"Timestep {i}")
-        plt.xlabel("reference el")
-        plt.ylabel("L2 error")
-        plt.ylim(0,0.05)
-        #plt.show()
-        plt.savefig(f"opt-err-{i}.png",dpi=60)
-        plt.clf()

examples/make-glvis-figs

@@ -1,19 +0,0 @@
-window 0 0 400 400
-
-# Initial solution
-solution soln-0.mesh soln-0-0.gf
-
-# Setup the GLVis scene. Executed after pressing the space bar.
-{
-   perspective off
-#   view 0 0
-   viewcenter 0.2 0
-   valuerange 0.0 1.0
-#   zoom 1.5
-   keys ca
-}
-
-# Take multiple screenshots. Executed after pressing the space bar.
-{
-#insert the screenshot commands here
-}

examples/make-periodic-mesh

@@ -1,63 +0,0 @@
-#!/bin/bash
-
-if [ $# -ne 1 ]; then
-    echo "Usage: ./make-periodic-mesh <NE>"
-fi
-
-NE=$1
-NEm1=$(( $NE - 1 ))
-NEm2=$(( $NE - 2 ))
-
-filename="periodic-segment-$NE.mesh"
-
-cat <<ENDHERE > $filename
-MFEM mesh v1.0
-
-#
-# MFEM Geometry Types (see mesh/geom.hpp):
-#
-# POINT       = 0
-# SEGMENT     = 1
-# TRIANGLE    = 2
-# SQUARE      = 3
-# TETRAHEDRON = 4
-# CUBE        = 5
-#
-
-dimension
-1
-
-# format: <attribute> <geometry type> <vertex 0> <vertex 1> ...
-elements
-ENDHERE
-
-echo "$NE" >> $filename
-for i in `seq 0 $NEm2`; do
-    echo "1 1 $i $((i+1))" >> $filename
-done
-echo "1 1 $((i+1)) 0" >> $filename
-
-cat <<EOF >> $filename
-
-boundary
-0
-
-vertices
-$NE
-
-nodes
-FiniteElementSpace
-FiniteElementCollection: L2_T1_1D_P1
-VDim: 1
-Ordering: 1
- 
-EOF
-
-dx=`echo "1.0 / $NE" | bc -l`
-for i in `seq 0 $NEm1`; do
-    p1=`echo "$i * $dx" | bc -l`
-    p2=`echo "$((i+1)) * $dx" | bc -l`
-    echo "$p1 $p2" >> $filename
-done
-
-

examples/periodic-segment-16.mesh

@@ -1,65 +0,0 @@
-MFEM mesh v1.0
-
-#
-# MFEM Geometry Types (see mesh/geom.hpp):
-#
-# POINT       = 0
-# SEGMENT     = 1
-# TRIANGLE    = 2
-# SQUARE      = 3
-# TETRAHEDRON = 4
-# CUBE        = 5
-#
-
-dimension
-1
-
-# format: <attribute> <geometry type> <vertex 0> <vertex 1> ...
-elements
-16
-1 1 0 1
-1 1 1 2
-1 1 2 3
-1 1 3 4
-1 1 4 5
-1 1 5 6
-1 1 6 7
-1 1 7 8
-1 1 8 9
-1 1 9 10
-1 1 10 11
-1 1 11 12
-1 1 12 13
-1 1 13 14
-1 1 14 15
-1 1 15 0
-
-
-boundary
-0
-
-vertices
-16
-
-nodes
-FiniteElementSpace
-FiniteElementCollection: L2_T1_1D_P1
-VDim: 1
-Ordering: 1
-
-0.0000 0.0625
-0.0625 0.1250
-0.1250 0.1875
-0.1875 0.2500
-0.2500 0.3125
-0.3125 0.3750
-0.3750 0.4375
-0.4375 0.5000
-0.5000 0.5625
-0.5625 0.6250
-0.6250 0.6875
-0.6875 0.7500
-0.7500 0.8125
-0.8125 0.8750
-0.8750 0.9375
-0.9375 1.0000

examples/periodic-segment-8.mesh

@@ -1,50 +0,0 @@
-MFEM mesh v1.0
-
-#
-# MFEM Geometry Types (see mesh/geom.hpp):
-#
-# POINT       = 0
-# SEGMENT     = 1
-# TRIANGLE    = 2
-# SQUARE      = 3
-# TETRAHEDRON = 4
-# CUBE        = 5
-#
-
-dimension
-1
-
-# format: <attribute> <geometry type> <vertex 0> <vertex 1> ...
-elements
-8
-1 1 0 1
-1 1 1 2
-1 1 2 3
-1 1 3 4
-1 1 4 5
-1 1 5 6
-1 1 6 7
-1 1 7 0
-
-boundary
-0
-
-vertices
-8
-
-nodes
-FiniteElementSpace
-FiniteElementCollection: L2_T1_1D_P1
-VDim: 1
-Ordering: 1
-
-0.000 0.125
-0.125 0.250
-0.250 0.375
-0.375 0.500
-0.500 0.625
-0.625 0.750
-0.750 0.875
-0.875 1.000
-
-

examples/run-search

@@ -1,43 +0,0 @@
-#!/bin/bash
-
-# sod 5 regrids
-options="-dt 0.01 --cfl-number -1 " # fixed timestep
-options+="--regrid-period 4 "       # regrid period
-options+="--problem 2 "             # Sod problem
-options+="--order 0 "               # 0th order DG
-options+="--refine 0 "              # no initial uniform refinement
-options+="--t-final 0.20 "
-
-# create reference solution on 16 element mesh
-./ex18 $options -m periodic-segment-20.mesh
-./save-reference
-
-ne=10
-nem1=$((ne-1))
-
-echo "" > err-results.dat
-for i1 in $(seq 0 $nem1); do
-    for i2 in $(seq 0 $nem1); do
-	for i3 in $(seq 0 $nem1); do
-	    for i4 in $(seq 0 $nem1); do
-		for i5 in $(seq 0 $nem1); do
-		    dir="seq-$i1$i2$i3$i4$i5"
-		    mkdir $dir
-		    cd $dir
-		    ln -s ../reference.mesh .
-		    ln -s ../reference-rho.gf .
-		    ln -s ../reference-rho-u.gf .
-		    ln -s ../reference-rho-e.gf .
-		    echo "running $i1$i2$i3$i4$i5"
-		    seq="$i1 $i2 $i3 $i4 $i5"
-		    cmd="../ex18 -rseq \"$seq\" $options -m ../periodic-segment-10.mesh | grep err | awk \"{print \$2}\""
-		    echo "$cmd" > runcode
-		    err=$( source runcode )
-		    cd ..
-		    echo "$i1 $i2 $i3 $i4 $i5 $err" >> err-results.dat
-		done
-	    done
-	done
-    done
-done
-

examples/save-reference

@@ -1,16 +0,0 @@
-#!/bin/bash
-
-# The final time is the reference solution
-mv vortex.mesh reference.mesh
-mv vortex-0-final.gf reference-rho.gf
-mv vortex-1-final.gf reference-rho-u.gf
-mv vortex-2-final.gf reference-rho-e.gf
-
-# Save the timesteps as well for errors
-for i in {0..19}
-do
-    mv soln-$i.mesh reference-$i.mesh
-    mv soln-$i-0.gf reference-$i-0.gf
-    mv soln-$i-1.gf reference-$i-1.gf
-    mv soln-$i-2.gf reference-$i-2.gf
-done

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);
    }
@@ -3495,6 +3623,8 @@ void Mesh::SetMeshGen()
             mesh_geoms |= (1 << Geometry::TETRAHEDRON);
          case Element::TRIANGLE:
             mesh_geoms |= (1 << Geometry::TRIANGLE);
+         case Element::SEGMENT:
+            mesh_geoms |= (1 << Geometry::SEGMENT);
          case Element::POINT:
             mesh_geoms |= (1 << Geometry::POINT);
             meshgen |= 1;
@@ -3506,8 +3636,6 @@ void Mesh::SetMeshGen()
             mesh_geoms |= (1 << Geometry::SQUARE);
             mesh_geoms |= (1 << Geometry::SEGMENT);
             mesh_geoms |= (1 << Geometry::POINT);
-         case Element::SEGMENT:
-            mesh_geoms |= (1 << Geometry::SEGMENT);
             meshgen |= 2;
             break;
 
@@ -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;
 }
@@ -4253,19 +4253,6 @@ void NCMesh::GetPointMatrix(Geometry::Type geom, const char* ref_path,
             pm = PointMatrix(mid12, mid20, mid01);
          }
       }
-      else if (geom == Geometry::SEGMENT)
-      {
-         Point mid01(pm(0), pm(1));
-
-         if (child == 0)
-         {
-            pm = PointMatrix(pm(0), mid01);
-         }
-         else if (child == 1)
-         {
-            pm = PointMatrix(mid01, pm(1));
-         }
-      }
    }
 
    // write the points to the matrix

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