Merge branch 'eigen-dev' of github.com:mfem/mfem into eigen-dev

Allow reading Gmsh meshes where all elements have attribute zero [gmsh-zero-attributes]

.gitignore

@@ -145,6 +145,14 @@ examples/petsc/velocity.*
 examples/petsc/elastic_energy.*
 examples/petsc/mode_*
 
+examples/arpack/ex11
+examples/arpack/mode_*
+examples/arpack/ex11.mesh
+
+examples/spectra/ex11
+examples/spectra/mode_*
+examples/spectra/ex11.mesh
+
 examples/pumi/ex1
 examples/pumi/ex[126]p
 examples/pumi/refined.mesh
@@ -310,7 +318,13 @@ tests/convergence/prates
 tests/par-mesh-format/ex1p
 
 # VPATH builds
-build-*/*
+build-*/
+
+# User config
+user-*
+
+# VSCode
+.vscode 
 
 # PETSc automated build
 petsc-build/*

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.
 
@@ -20,6 +27,9 @@ Version 4.3.1 (development)
   functions on wedges and pyramids which are not amenable to reordering. The
   ReorientTetMesh method of the Mesh and ParMesh classes has been deprecated.
 
+- Gmsh meshes where all elements have zero physical tag (the default Gmsh
+  output format if no physical groups are defined) are now successfully loaded,
+  and elements are reassigned attribute number 1.
 
 Version 4.3, released on July 29, 2021
 ======================================

config/config.hpp.in

@@ -91,6 +91,12 @@
 // Enable MFEM functionality based on the SuiteSparse library.
 // #define MFEM_USE_SUITESPARSE
 
+// Enable MFEM functionality based on the ARPACK library.
+// #define MFEM_USE_ARPACK
+
+// Enable MFEM functionality based on the SPECTRA library.
+// #define MFEM_USE_SPECTRA
+
 // Enable MFEM functionality based on the SuperLU library.
 // #define MFEM_USE_SUPERLU
 // #define MFEM_USE_SUPERLU5

config/config.mk.in

@@ -31,6 +31,8 @@ MFEM_TIMER_TYPE        = @MFEM_TIMER_TYPE@
 MFEM_USE_SUNDIALS      = @MFEM_USE_SUNDIALS@
 MFEM_USE_MESQUITE      = @MFEM_USE_MESQUITE@
 MFEM_USE_SUITESPARSE   = @MFEM_USE_SUITESPARSE@
+MFEM_USE_ARPACK        = @MFEM_USE_ARPACK@
+MFEM_USE_SPECTRA       = @MFEM_USE_SPECTRA@
 MFEM_USE_SUPERLU       = @MFEM_USE_SUPERLU@
 MFEM_USE_SUPERLU5      = @MFEM_USE_SUPERLU5@
 MFEM_USE_MUMPS         = @MFEM_USE_MUMPS@

config/defaults.mk

@@ -151,6 +151,8 @@ MFEM_USE_UMPIRE        = NO
 MFEM_USE_SIMD          = NO
 MFEM_USE_ADIOS2        = NO
 MFEM_USE_MKL_CPARDISO  = NO
+MFEM_USE_ARPACK        = NO
+MFEM_USE_SPECTRA       = NO
 
 # MPI library compile and link flags
 # These settings are used only when building MFEM with MPI + HIP
@@ -328,6 +330,19 @@ NETCDF_LIB = $(XLINKER)-rpath,$(NETCDF_DIR)/lib -L$(NETCDF_DIR)/lib\
  $(XLINKER)-rpath,$(HDF5_DIR)/lib -L$(HDF5_DIR)/lib\
  -lnetcdf -lhdf5_hl -lhdf5 $(ZLIB_LIB)
 
+# ARPACK library configuration
+ARPACK_DIR = @MFEM_DIR@/../ARPACK
+ARPACK_OPT = -I$(ARPACK_DIR)
+ARPACK_LIB = -L$(ARPACK_DIR) -lparpack -larpack
+
+# EIGEN library configuration
+EIGEN_DIR = @MFEM_DIR@/../eigen
+EIGEN_OPT = -I$(EIGEN_DIR)
+
+# SPECTRA library configuration
+SPECTRA_DIR = @MFEM_DIR@/../spectra/include
+SPECTRA_OPT = -I$(SPECTRA_DIR) $(EIGEN_OPT)
+
 # PETSc library configuration (version greater or equal to 3.8 or the dev branch)
 PETSC_ARCH := arch-linux2-c-debug
 PETSC_DIR  := $(MFEM_DIR)/../petsc/$(PETSC_ARCH)

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

data/sphere.geo

@@ -0,0 +1,47 @@
+Mesh.Algorithm = 6;
+
+lc = 0.1;
+Point(1) = {0.0,0.0,0.0,lc};
+Point(2) = {1,0.0,0.0,lc};
+Point(3) = {0,1,0.0,lc};
+Circle(1) = {2,1,3};
+Point(4) = {-1,0,0.0,lc};
+Point(5) = {0,-1,0.0,lc};
+Circle(2) = {3,1,4};
+Circle(3) = {4,1,5};
+Circle(4) = {5,1,2};
+Point(6) = {0,0,-1,lc};
+Point(7) = {0,0,1,lc};
+Circle(5) = {3,1,6};
+Circle(6) = {6,1,5};
+Circle(7) = {5,1,7};
+Circle(8) = {7,1,3};
+Circle(9) = {2,1,7};
+Circle(10) = {7,1,4};
+Circle(11) = {4,1,6};
+Circle(12) = {6,1,2};
+Curve Loop(13) = {2,8,-10};
+Surface(14) = {13};
+Curve Loop(15) = {10,3,7};
+Surface(16) = {15};
+Curve Loop(17) = {-8,-9,1};
+Surface(18) = {17};
+Curve Loop(19) = {-11,-2,5};
+Surface(20) = {19};
+Curve Loop(21) = {-5,-12,-1};
+Surface(22) = {21};
+Curve Loop(23) = {-3,11,6};
+Surface(24) = {23};
+Curve Loop(25) = {-7,4,9};
+Surface(26) = {25};
+Curve Loop(27) = {-4,12,-6};
+Surface(28) = {27};
+Surface Loop(29) = {28,26,16,14,20,24,22,18};
+Volume(30) = {29};
+
+Physical Surface(1) = {28,26,16,14,20,24,22,18};
+Physical Volume(2) = 30;
+
+// Generate 2D mesh
+Mesh 2;
+Mesh.MshFileVersion = 2.2;

examples/arpack/ex11.cpp

@@ -0,0 +1,286 @@
+//                       MFEM Example 11 - Serial Version
+//
+// Compile with: make ex11
+//
+// Sample runs:  ex11 -m ../data/square-disc.mesh
+//               ex11 -m ../data/star.mesh
+//               ex11 -m ../data/star-mixed.mesh
+//               ex11 -m ../data/periodic-annulus-sector.msh
+//               ex11 -m ../data/square-disc-p2.vtk -o 2
+//               ex11 -m ../data/square-disc-p3.mesh -o 3
+//               ex11 -m ../data/square-disc-nurbs.mesh -o -1
+//               ex11 -m ../data/disc-nurbs.mesh -o -1 -n 20
+//               ex11 -m ../data/star-surf.mesh
+//               ex11 -m ../data/square-disc-surf.mesh
+//               ex11 -m ../data/inline-segment.mesh
+//               ex11 -m ../data/inline-quad.mesh
+//               ex11 -m ../data/inline-tri.mesh
+//               ex11 -m ../data/amr-quad.mesh
+//               ex11 -m ../data/amr-hex.mesh
+//               ex11 -m ../data/mobius-strip.mesh -n 8
+//
+// Description:  This example code demonstrates the use of MFEM to solve the
+//               eigenvalue problem -Delta u = lambda u with homogeneous
+//               Dirichlet boundary conditions.
+//
+//               We compute a number of the lowest eigenmodes by discretizing
+//               the Laplacian and Mass operators using a FE space of the
+//               specified order, or an isoparametric/isogeometric space if
+//               order < 1 (quadratic for quadratic curvilinear mesh, NURBS for
+//               NURBS mesh, etc.)
+//
+//               The example highlights the use of the ARPACK eigenvalue solver
+//               (regular inverse mode). Reusing a single GLVis visualization
+//               window for multiple eigenfunctions is also illustrated.
+//
+//               We recommend viewing Example 1 before viewing this example.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file = "../../data/star.mesh";
+   int ser_ref_levels = 3;
+   int order = 1;
+   int nev = 5;
+   double dbc_eig = 1e3;
+   bool visualization = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
+                  "Number of times to refine the mesh uniformly in serial.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&nev, "-n", "--num-eigs",
+                  "Number of desired eigenmodes.");
+   args.AddOption(&dbc_eig, "-d", "--dbc-eig",
+                  "Eigenvalues associated with Dirichlet BC "
+                  "(should be larger than the maximum desired eigenvalue).");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   // 2. Read the (serial) mesh from the given mesh file on all processors. We
+   //    can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
+   //    and volume meshes with the same code.
+   Mesh *mesh;
+   ifstream imesh(mesh_file);
+   if (!imesh)
+   {
+      cerr << "\nCan not open mesh file: " << mesh_file << '\n' << endl;
+      return 2;
+   }
+   mesh = new Mesh(imesh, 1, 1);
+   imesh.close();
+   int dim = mesh->Dimension();
+
+   // 3. Refine the serial mesh on all processors to increase the resolution. In
+   //    this example we do 'ref_levels' of uniform refinement (2 by default, or
+   //    specified on the command line with -rs).
+   for (int lev = 0; lev < ser_ref_levels; lev++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   // 4. Define a finite element space on the mesh. Here we
+   //    use continuous Lagrange finite elements of the specified order. If
+   //    order < 1, we instead use an isoparametric/isogeometric space.
+   FiniteElementCollection *fec;
+   if (order > 0)
+   {
+      fec = new H1_FECollection(order, dim);
+   }
+   else if (mesh->GetNodes())
+   {
+      fec = mesh->GetNodes()->OwnFEC();
+   }
+   else
+   {
+      fec = new H1_FECollection(order = 1, dim);
+   }
+   FiniteElementSpace *fespace = new FiniteElementSpace(mesh, fec);
+   int size = fespace->GetVSize();
+
+   cout << "Number of unknowns: " << size << endl;
+
+   // 5. Set up the parallel bilinear forms a(.,.) and m(.,.) on the finite
+   //    element space. The first corresponds to the Laplacian operator -Delta,
+   //    while the second is a simple mass matrix needed on the right hand side
+   //    of the generalized eigenvalue problem below. The boundary conditions
+   //    are implemented by elimination with special values on the diagonal to
+   //    shift the Dirichlet eigenvalues out of the computational range. After
+   //    serial and parallel assembly we extract the corresponding parallel
+   //    matrices A and M.
+   ConstantCoefficient one(1.0);
+   Array<int> ess_bdr;
+   if (mesh->bdr_attributes.Size())
+   {
+      ess_bdr.SetSize(mesh->bdr_attributes.Max());
+      ess_bdr = 1;
+   }
+
+   BilinearForm *a = new BilinearForm(fespace);
+   a->AddDomainIntegrator(new DiffusionIntegrator(one));
+   if (mesh->bdr_attributes.Size() == 0)
+   {
+      // Add a mass term if the mesh has no boundary, e.g. periodic mesh or
+      // closed surface.
+      a->AddDomainIntegrator(new MassIntegrator(one));
+   }
+   a->Assemble();
+   if (mesh->bdr_attributes.Size() != 0)
+   {
+      a->EliminateEssentialBCDiag(ess_bdr, dbc_eig);
+   }
+   a->Finalize();
+
+   BilinearForm *m = new BilinearForm(fespace);
+   m->AddDomainIntegrator(new MassIntegrator(one));
+   m->Assemble();
+   if (mesh->bdr_attributes.Size() != 0)
+   {
+      // shift the eigenvalue corresponding to eliminated dofs to a large value
+      m->EliminateEssentialBCDiag(ess_bdr, 1.0);
+   }
+   m->Finalize();
+
+   // 6. Define and configure the ARPACK eigensolver
+   ArPackSym * arpack = new ArPackSym();
+   Solver    * solver = NULL;
+
+#ifndef MFEM_USE_SUITESPARSE
+   // 7. Define a simple symmetric Gauss-Seidel preconditioner and use it to
+   //    solve the system A X = B with PCG.
+   cout << "Building CGSolver" << endl;
+   GSSmoother M(m->SpMat());
+   CGSolver * cg_solver = new CGSolver;
+   cg_solver->SetPreconditioner(M);
+   cg_solver->SetRelTol(1.0e-12);
+   solver = cg_solver;
+#else
+   // 7. If MFEM was compiled with SuiteSparse, use UMFPACK to solve the system.
+   cout << "Building UMFPackSolver" << endl;
+   UMFPackSolver * umf_solver = new UMFPackSolver;
+   umf_solver->Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
+   solver = umf_solver;
+#endif
+   solver->SetOperator(m->SpMat());
+
+   arpack->SetNumModes(nev);
+   arpack->SetMaxIter(400);
+   arpack->SetTol(1e-8);
+   arpack->SetMode(2);
+   arpack->SetPrintLevel(2);
+
+   arpack->SetOperator(*a);
+   arpack->SetMassMatrix(*m);
+   arpack->SetSolver(*solver);
+
+   // 8. Compute the eigenmodes and extract the array of eigenvalues. Define a
+   //    parallel grid function to represent each of the eigenmodes returned by
+   //    the solver.
+   Array<double> eigenvalues;
+   arpack->Solve();
+   arpack->GetEigenvalues(eigenvalues);
+
+   cout << endl;
+   std::ios::fmtflags old_fmt = cout.flags();
+   cout.setf(std::ios::scientific);
+   std::streamsize old_prec = cout.precision(14);
+   for (int i=0; i<nev; i++)
+   {
+      cout << "Eigenvalue lambda   " << eigenvalues[i] << endl;
+   }
+   cout.precision(old_prec);
+   cout.flags(old_fmt);
+   cout << endl;
+
+   GridFunction x(fespace);
+
+   // 9. Save the refined mesh and the modes in parallel. This output can be
+   //    viewed later using GLVis: "glvis -np <np> -m mesh -g mode".
+   {
+      ostringstream mesh_name, mode_name;
+      mesh_name << "ex11.mesh";
+
+      ofstream mesh_ofs(mesh_name.str().c_str());
+      mesh_ofs.precision(8);
+      mesh->Print(mesh_ofs);
+
+      for (int i=0; i<nev; i++)
+      {
+         // convert eigenvector from HypreParVector to ParGridFunction
+         x = arpack->GetEigenvector(i);
+
+         mode_name << "mode_" << setfill('0') << setw(2) << i;
+
+         ofstream mode_ofs(mode_name.str().c_str());
+         mode_ofs.precision(8);
+         x.Save(mode_ofs);
+         mode_name.str("");
+      }
+   }
+
+   // 10. Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream mode_sock(vishost, visport);
+      mode_sock.precision(8);
+
+      for (int i=0; i<nev; i++)
+      {
+         cout << "Eigenmode " << i+1 << '/' << nev
+              << ", Lambda = " << eigenvalues[i] << endl;
+
+         // convert eigenvector from HypreParVector to ParGridFunction
+         x = arpack->GetEigenvector(i);
+
+         mode_sock << "solution\n" << *mesh << x << flush
+                   << "window_title 'Eigenmode " << i+1 << '/' << nev
+                   << ", Lambda = " << eigenvalues[i] << "'" << endl;
+
+         char c;
+         cout << "press (q)uit or (c)ontinue --> " << flush;
+         cin >> c;
+
+         if (c != 'c')
+         {
+            break;
+         }
+      }
+      mode_sock.close();
+   }
+
+   // 11. Free the used memory.
+   delete arpack;
+   delete solver;
+   delete m;
+   delete a;
+
+   delete fespace;
+   if (order > 0)
+   {
+      delete fec;
+   }
+   delete mesh;
+
+   return 0;
+}

examples/arpack/makefile

@@ -0,0 +1,69 @@
+# 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.
+
+# Use the MFEM build directory
+MFEM_DIR ?= ../..
+MFEM_BUILD_DIR ?= ../..
+SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/examples/arpack/,)
+CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
+# Use the MFEM install directory
+# MFEM_INSTALL_DIR = ../../mfem
+# CONFIG_MK = $(MFEM_INSTALL_DIR)/share/mfem/config.mk
+
+MFEM_LIB_FILE = mfem_is_not_built
+-include $(CONFIG_MK)
+
+SEQ_EXAMPLES = ex11
+PAR_EXAMPLES = 
+ifeq ($(MFEM_USE_MPI),NO)
+   EXAMPLES = $(SEQ_EXAMPLES)
+else
+   EXAMPLES = $(PAR_EXAMPLES)
+endif
+RC_FILES = $(patsubst $(SRC)%,%,$(wildcard $(SRC)rc_*))
+
+.SUFFIXES:
+.SUFFIXES: .o .cpp .mk
+.PHONY: all clean clean-build clean-exec
+
+# Remove built-in rule
+%: %.cpp
+
+# Replace the default implicit rule for *.cpp files
+%: $(SRC)%.cpp $(MFEM_LIB_FILE) $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ $(MFEM_LIBS)
+
+all: $(EXAMPLES)
+
+# Examples depend on their corresponding rc_* files:
+make-rc-rule = $(1): | $(filter rc_$(1)%,$(RC_FILES))
+$(foreach ex,$(EXAMPLES),$(eval $(call make-rc-rule,$(ex))))
+
+# Rules to copy the rc_* files when building out-of-source:
+ifneq ($(SRC),)
+$(RC_FILES): %: $(SRC)%
+	cp -pf $(<) .
+endif
+
+# Generate an error message if the MFEM library is not built and exit
+$(MFEM_LIB_FILE):
+	$(error The MFEM library is not built)
+
+clean: clean-build clean-exec
+
+clean-build:
+	rm -f *.o *~ $(SEQ_EXAMPLES) $(PAR_EXAMPLES)
+	rm -rf *.dSYM *.TVD.*breakpoints
+
+clean-exec:
+	@rm -rf mesh.* sol.* sol_p.* sol_u.* Example5*
+	@rm -f ex9-mesh.* ex9-init.* ex9-final.* Example9*
+	@rm -f deformed.* velocity.* elastic_energy.*

examples/ex1.cpp

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

examples/ex1p.cpp

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

examples/ex22.cpp

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

examples/ex22p.cpp

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

examples/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/spectra/ex11.cpp

@@ -0,0 +1,250 @@
+//                       MFEM Example 11 - Serial Version
+//
+// Compile with: make ex11
+//
+// Sample runs:  ex11 -m ../data/square-disc.mesh
+//               ex11 -m ../data/star.mesh
+//               ex11 -m ../data/star-mixed.mesh
+//               ex11 -m ../data/periodic-annulus-sector.msh
+//               ex11 -m ../data/square-disc-p2.vtk -o 2
+//               ex11 -m ../data/square-disc-p3.mesh -o 3
+//               ex11 -m ../data/square-disc-nurbs.mesh -o -1
+//               ex11 -m ../data/disc-nurbs.mesh -o -1 -n 20
+//               ex11 -m ../data/star-surf.mesh
+//               ex11 -m ../data/square-disc-surf.mesh
+//               ex11 -m ../data/inline-segment.mesh
+//               ex11 -m ../data/inline-quad.mesh
+//               ex11 -m ../data/inline-tri.mesh
+//               ex11 -m ../data/amr-quad.mesh
+//               ex11 -m ../data/amr-hex.mesh
+//               ex11 -m ../data/mobius-strip.mesh -n 8
+//
+// Description:  This example code demonstrates the use of MFEM to solve the
+//               eigenvalue problem -Delta u = lambda u with homogeneous
+//               Dirichlet boundary conditions.
+//
+//               We compute a number of the lowest eigenmodes by discretizing
+//               the Laplacian and Mass operators using a FE space of the
+//               specified order, or an isoparametric/isogeometric space if
+//               order < 1 (quadratic for quadratic curvilinear mesh, NURBS for
+//               NURBS mesh, etc.)
+//
+//               The example highlights the use of the ARPACK eigenvalue solver
+//               (regular inverse mode). Reusing a single GLVis visualization
+//               window for multiple eigenfunctions is also illustrated.
+//
+//               We recommend viewing Example 1 before viewing this example.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file = "../../data/star.mesh";
+   int ser_ref_levels = 1;
+   int order = 1;
+   int nev = 5;
+   double dbc_eig = 1e3;
+   bool visualization = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
+                  "Number of times to refine the mesh uniformly in serial.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&nev, "-n", "--num-eigs",
+                  "Number of desired eigenmodes.");
+   args.AddOption(&dbc_eig, "-d", "--dbc-eig",
+                  "Eigenvalues associated with Dirichlet BC "
+                  "(should be larger than the maximum desired eigenvalue).");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   // 2. Read the (serial) mesh from the given mesh file on all processors. We
+   //    can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
+   //    and volume meshes with the same code.
+   Mesh *mesh;
+   ifstream imesh(mesh_file);
+   if (!imesh)
+   {
+      cerr << "\nCan not open mesh file: " << mesh_file << '\n' << endl;
+      return 2;
+   }
+   mesh = new Mesh(imesh, 1, 1);
+   imesh.close();
+   int dim = mesh->Dimension();
+
+   // 3. Refine the serial mesh on all processors to increase the resolution. In
+   //    this example we do 'ref_levels' of uniform refinement (2 by default, or
+   //    specified on the command line with -rs). 
+   for (int lev = 0; lev < ser_ref_levels; lev++)
+   {
+      mesh->UniformRefinement();
+   }
+
+   // 4. Define a finite element space on the mesh. Here we
+   //    use continuous Lagrange finite elements of the specified order. If
+   //    order < 1, we instead use an isoparametric/isogeometric space.
+   FiniteElementCollection *fec;
+   if (order > 0)
+   {
+      fec = new H1_FECollection(order, dim);
+   }
+   else if (mesh->GetNodes())
+   {
+      fec = mesh->GetNodes()->OwnFEC();
+   }
+   else
+   {
+      fec = new H1_FECollection(order = 1, dim);
+   }
+   FiniteElementSpace *fespace = new FiniteElementSpace(mesh, fec);
+   int size = fespace->GetVSize();
+
+   cout << "Number of unknowns: " << size << endl;
+
+   // 5. Set up the parallel bilinear forms a(.,.) and m(.,.) on the finite
+   //    element space. The first corresponds to the Laplacian operator -Delta,
+   //    while the second is a simple mass matrix needed on the right hand side
+   //    of the generalized eigenvalue problem below. The boundary conditions
+   //    are implemented by elimination with special values on the diagonal to
+   //    shift the Dirichlet eigenvalues out of the computational range. After
+   //    serial and parallel assembly we extract the corresponding parallel
+   //    matrices A and M.
+   ConstantCoefficient one(1.0);
+   Array<int> ess_bdr;
+   if (mesh->bdr_attributes.Size())
+   {
+      ess_bdr.SetSize(mesh->bdr_attributes.Max());
+      ess_bdr = 1;
+   }
+
+   BilinearForm *a = new BilinearForm(fespace);
+   a->AddDomainIntegrator(new DiffusionIntegrator(one));
+   if (mesh->bdr_attributes.Size() == 0)
+   {
+      // Add a mass term if the mesh has no boundary, e.g. periodic mesh or
+      // closed surface.
+      a->AddDomainIntegrator(new MassIntegrator(one));
+   }
+   a->Assemble();
+   if (mesh->bdr_attributes.Size() != 0)
+   {
+      a->EliminateEssentialBCDiag(ess_bdr, dbc_eig);
+   }
+   a->Finalize();
+
+   BilinearForm *m = new BilinearForm(fespace);
+   m->AddDomainIntegrator(new MassIntegrator(one));
+   m->Assemble();
+   if (mesh->bdr_attributes.Size() != 0)
+   {
+      // shift the eigenvalue corresponding to eliminated dofs to a large value
+      m->EliminateEssentialBCDiag(ess_bdr, 1.0);
+   }
+   m->Finalize();
+
+   // 6. Define and configure the SPECTRA eigensolver and solve problem
+   SpectraEigenSolver spectra;
+
+   spectra.SetNumModes(nev)
+      .SetKrylov(10)
+      .SetMaxIter(5000)
+      .SetTol(1e-5)
+      .SetOperators(*a, *m)
+      .Solve();
+
+   Eigen::VectorXd eigenvalues = spectra.GetEigenvalues(nev);
+
+   // 7. Define a grid function to represent each of the eigenmodes returned by the solver.
+   
+   GridFunction x(fespace);
+
+   // 8. Save the refined mesh and the modes in parallel. 
+   //    This output can be viewed later using GLVis: "glvis -np <np> -m mesh -g mode"
+   {
+      ostringstream mesh_name, mode_name;
+      mesh_name << "ex11.mesh";
+
+      ofstream mesh_ofs(mesh_name.str().c_str());
+      mesh_ofs.precision(8);
+      mesh->Print(mesh_ofs);
+
+      for (int i = 0; i < nev; i++) {
+         // conver Eigen Vector to MFEM Vector
+         Vector eigenvector = VectorConverter<double>::from(spectra.GetEigenvector(i));
+         
+         // convert eigenvector from Vector to GridFunction
+         x = eigenvector;
+
+         mode_name << "mode_" << setfill('0') << setw(2) << i;
+
+         ofstream mode_ofs(mode_name.str().c_str());
+         mode_ofs.precision(8);
+         x.Save(mode_ofs);
+         mode_name.str("");
+      }
+   }
+
+   // 10. Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream mode_sock(vishost, visport);
+      mode_sock.precision(8);
+
+      for (int i=0; i<nev; i++)
+      {
+         cout << "Eigenmode " << i+1 << '/' << nev
+              << ", Lambda = " << eigenvalues[i] << endl;
+
+         // convert eigenvector from HypreParVector to ParGridFunction
+         Vector eigenvector = VectorConverter<double>::from(spectra.GetEigenvector(i));
+         x = eigenvector;
+
+         mode_sock << "solution\n" << *mesh << x << flush
+                   << "window_title 'Eigenmode " << i+1 << '/' << nev
+                   << ", Lambda = " << eigenvalues[i] << "'" << endl;
+
+         char c;
+         cout << "press (q)uit or (c)ontinue --> " << flush;
+         cin >> c;
+
+         if (c != 'c')
+         {
+            break;
+         }
+      }
+      mode_sock.close();
+   }
+
+   // 10. Free the used memory.
+   delete m;
+   delete a;
+
+   delete fespace;
+   if (order > 0)
+   {
+      delete fec;
+   }
+   delete mesh;
+
+   return 0;
+}

examples/spectra/makefile

@@ -0,0 +1,67 @@
+# 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.
+
+# Use the MFEM build directory
+MFEM_DIR ?= ../..
+MFEM_BUILD_DIR ?= ../..
+SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/examples/spectra/,)
+CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
+# Use the MFEM install directory
+# MFEM_INSTALL_DIR = ../../mfem
+# CONFIG_MK = $(MFEM_INSTALL_DIR)/share/mfem/config.mk
+
+MFEM_LIB_FILE = mfem_is_not_built
+-include $(CONFIG_MK)
+
+SEQ_EXAMPLES = ex11
+PAR_EXAMPLES = 
+ifeq ($(MFEM_USE_MPI),NO)
+   EXAMPLES = $(SEQ_EXAMPLES)
+else
+   EXAMPLES = $(PAR_EXAMPLES)
+endif
+RC_FILES = $(patsubst $(SRC)%,%,$(wildcard $(SRC)rc_*))
+
+.SUFFIXES:
+.SUFFIXES: .o .cpp .mk
+.PHONY: all clean clean-build clean-exec
+
+# Remove built-in rule
+%: %.cpp
+
+# Replace the default implicit rule for *.cpp files
+%: $(SRC)%.cpp $(MFEM_LIB_FILE) $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ $(MFEM_LIBS)
+
+all: $(EXAMPLES)
+
+# Examples depend on their corresponding rc_* files:
+make-rc-rule = $(1): | $(filter rc_$(1)%,$(RC_FILES))
+$(foreach ex,$(EXAMPLES),$(eval $(call make-rc-rule,$(ex))))
+
+# Rules to copy the rc_* files when building out-of-source:
+ifneq ($(SRC),)
+$(RC_FILES): %: $(SRC)%
+	cp -pf $(<) .
+endif
+
+# Generate an error message if the MFEM library is not built and exit
+$(MFEM_LIB_FILE):
+	$(error The MFEM library is not built)
+
+clean: clean-build clean-exec
+
+clean-build:
+	rm -f *.o *~ $(SEQ_EXAMPLES) $(PAR_EXAMPLES)
+	rm -rf *.dSYM *.TVD.*breakpoints
+
+clean-exec:
+	@rm -rf *.mesh mode_*

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

linalg/CMakeLists.txt

@@ -78,6 +78,16 @@ if (MFEM_USE_MPI)
   endif()
 endif()
 
+if (MFEM_USE_ARPACK)
+  list(APPEND SRCS eigensolvers.cpp arpack.cpp)
+  list(APPEND HDRS eigensolvers.hpp arpack.hpp)
+endif()
+
+if (MFEM_USE_SPECTRA)
+  list(APPEND SRCS spectra.cpp)
+  list(APPEND HDRS eigen.hpp spectra.hpp)
+endif()
+
 if (MFEM_USE_SUNDIALS)
   list(APPEND SRCS sundials.cpp)
   list(APPEND HDRS sundials.hpp)

linalg/arpack.hpp

@@ -0,0 +1,240 @@
+// Copyright (c) 2010-2020, 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_ARPACK
+#define MFEM_ARPACK
+
+#include "../config/config.hpp"
+
+#ifdef MFEM_USE_ARPACK
+
+#include <string>
+
+using namespace std;
+
+#ifdef MFEM_USE_MPI
+#include <mpi.h>
+#include "hypre.hpp"
+#endif
+
+#include "operator.hpp"
+
+#define DSAUPD  dsaupd_
+#define DSEUPD  dseupd_
+
+#ifdef MFEM_USE_MPI
+#define PDSAUPD pdsaupd_
+#define PDSEUPD pdseupd_
+#endif
+
+extern "C" void DSAUPD(int *ido,char *bmat, int *n,
+                       char *which, int *nev,double *tol,double *resid,
+                       int *ncv,double *v, int *ldv,
+                       int *iparam, int *ipntr,
+                       double *workd, double *workl, int *lworkl, int *info);
+
+extern "C" void DSEUPD(int *, char *,int *, double *,
+                       double *,int *, double *,char *, int *, char *,
+                       int *,double *,double *,int *, double *,
+                       int *, int  *,int *, double *,
+                       double *,int *, int *);
+
+#ifdef MFEM_USE_MPI
+
+extern "C" void PDSAUPD(int *comm, int *ido,char *bmat, int *n,
+                        char *which, int *nev,double *tol,double *resid,
+                        int *ncv,double *v, int *ldv,
+                        int *iparam, int *ipntr,
+                        double *workd, double *workl, int *lworkl, int *info);
+
+extern "C" void PDSEUPD(int *comm, int *, char *,int *, double *,
+                        double *,int *, double *,char *, int *, char *,
+                        int *,double *,double *,int *, double *,
+                        int *, int  *,int *, double *,
+                        double *,int *, int *);
+
+#endif
+
+extern "C" {
+   void arpackgetcommdbg_(int *,int *,int *);
+   void arpacksetcommdbg_(int *,int *,int *);
+   void arpacksymdbg_(int *,int *,int *,int *,int *,int *,int *);
+   void arpacknonsymdbg_(int *,int *,int *,int *,int *,int *,int *);
+   void arpackcmplxdbg_(int *,int *,int *,int *,int *,int *,int *);
+}
+
+namespace mfem
+{
+
+class ArPackSym : public Eigensolver
+{
+public:
+
+   ArPackSym();
+   virtual ~ArPackSym();
+
+   /** ARPACK modes are described in section 3.5 of the ARPACK manual.
+       Mode 1: regular mode to solve A x = lambda x
+               No solver and no mass matrix are needed.
+       Mode 2: regular inverse mode to solve A x = lambda M x
+               Both A and M are needed and the solver should compute M^{-1}.
+       Mode 3: shift-invert mode to solve either A x = lambda x
+               or A x = lambda M x
+               Mass matrix is optional.  The solver should compute
+               (A-sigma I)^{-1} or (A-sigma M)^{-1}.  The shift parameter,
+               sigma, also needs to be set with SetShift().
+       Mode 4: Buckling mode to solve K x = lambda K_G x
+               K is set using SetMassMatrix(), K_G is set using SetOperator(),
+          and the solver should compute (K-sigma K_G)^{-1}.  The shift
+               parameter, sigma, also needs to be set with SetShift().
+       Mode 5: Cayley mode to solve A x = lambda M x
+               Both A and M are needed and the solver should compute
+               (A - sigma M)^{-1}.  The shift parameter, sigma, also needs
+               to be set with SetShift().
+    */
+   void SetMode(int mode);
+
+   inline void SetTol(double tol)         {      tol_ = tol;      }
+   inline void SetMaxIter(int max_iter)   { max_iter_ = max_iter; }
+   inline void SetPrintLevel(int logging) {  logging_ = logging;  }
+   inline void SetShift(double sigma)     {    sigma_ = sigma;    }
+   inline void SetNumModes(int num_eigs)  {      nev_ = num_eigs; }
+
+   virtual void SetSolver(Solver & solver);
+   virtual void SetOperator(Operator & A);
+   virtual void SetMassMatrix(Operator & M);
+
+   void Solve();
+
+   /// Collect the converged eigenvalues
+   virtual void GetEigenvalues(Array<double> & eigenvalues);
+
+   /// Extract a single eigenvector
+   virtual Vector & GetEigenvector(unsigned int i);
+
+   /// Transfer ownership of the converged eigenvectors
+   Vector ** StealEigenvectors();
+
+protected:
+
+   int myid_;       // Index of this processor
+   int max_iter_;
+   int logging_;
+
+   // The following variables are for ARPACK
+   int nloc_;       // number of items stored locally
+   int nev_;        // number of requested eigenvalues
+   int ncv_;        // number of ritz vectors
+   int rvec_;       // boolean to return eigenvectors as well
+   int mode_;       // 1 = standard, 2 = generalized, 3 = shift invert,
+   // 4 = buckling, 5 = Cayley
+   int lworkl_;     // length of lworkl_ work array
+   int iparam_[12]; // arpack parameters
+   int ipntr_[12];  // arpack pointers
+
+   char bmat_;      // I for standard problem, G for generalized
+   char which_[3];  // spectrum portion: LA, SA, LM, SM, BE
+   char hwmny_;     // DSEUPD: A for all eigenvalues, S for some
+
+   double tol_;     // relative accuracy bound for Ritz values
+   double sigma_;   // eigenvalue shift parameter
+
+   int    * select_;// workspace used during eigenvalue computation
+   double * dv_;    // Ritz values
+   double * v_;     // ncv Lanczos basis vectors
+   double * resid_; // residual vector
+   double * workd_; // work array for 3 vectors used in Arnoldi iteration
+   double * workl_; // work array
+
+   // Operators and Vectors needed outside of ARPACK
+   Solver   * solver_;
+   Operator * A_;
+   Operator * B_;
+
+   Vector * w_;
+   Vector * x_;
+   Vector * y_;
+   Vector * z_;
+
+   Vector ** eigenvectors_;
+
+   string solverName_;
+
+   void reverseComm();
+
+   int reverseCommMode1();
+   int reverseCommMode2();
+   int reverseCommMode3();
+   int reverseCommMode4();
+   int reverseCommMode5();
+
+   virtual void prepareEigenvectors();
+
+   void printErrors(const int & info, const int iparam[],
+                    const char & bmat, const int & n,
+                    const char which[],
+                    const int & nev, const int & ncv,
+                    const int & lworkl );
+
+private:
+
+   virtual int computeNlocf() { return nloc_; }
+   virtual int computeIter(int & ido);
+   virtual int computeEigs();
+
+};
+
+#ifdef MFEM_USE_MPI
+
+class ParArPackSym : public ArPackSym
+{
+public:
+   ParArPackSym(MPI_Comm comm);
+   virtual ~ParArPackSym() {}
+
+   void SetOperator(Operator & A);
+   void SetMassMatrix(Operator & M);
+
+   /// Collect the converged eigenvalues
+   void GetEigenvalues(Array<double> & eigenvalues);
+
+   /// Extract a single eigenvector
+   Vector & GetEigenvector(unsigned int i);
+
+   /// Transfer ownership of the converged eigenvectors
+   // HypreParVector ** StealEigenvectors();
+   Vector ** StealEigenvectors();
+
+protected:
+
+   void prepareEigenvectors();
+
+private:
+
+   MPI_Comm comm_;
+   MPI_Fint commf_; // Fortran style MPI communicator
+   int numProcs_;   // Number of processors
+
+   HYPRE_Int * part_; // parallel partitioning for eigenvectors
+
+   int computeNlocf();
+   int computeIter(int & ido);
+   int computeEigs();
+
+};
+
+#endif // MFEM_USE_MPI
+
+};
+
+#endif // MFEM_USE_ARPACK
+
+#endif // MFEM_ARPACK

linalg/eigen.hpp

@@ -0,0 +1,94 @@
+#ifndef MFEM_EIGEN_HPP
+#define MFEM_EIGEN_HPP
+
+#include <vector>
+#include <Eigen/Sparse>
+
+#include "vector.hpp"
+#include "sparsemat.hpp"
+#include "densemat.hpp"
+
+namespace mfem{
+    /** @brief Eigen template specialization for vector conversion */
+    template <typename T>
+    struct VectorConverter {
+        static Vector from(const Eigen::Matrix<T, Eigen::Dynamic, 1>& other)
+        {
+            Vector v(other.rows());
+
+            for (size_t i = 0; i < v.Size(); i++)
+                v(i) = other(i);
+
+            return std::move(v);
+        }
+
+        static Eigen::Matrix<T, Eigen::Dynamic, 1> to(const Vector& other)
+        {
+            Eigen::Matrix<T, Eigen::Dynamic, 1> v(other.Size());
+
+            for (size_t i = 0; i < v.Size(); i++)
+                v(i) = other(i);
+
+            return std::move(v);
+        }
+    };
+
+    /** @brief Eigen template specialization for dense matrix conversion */
+    template <typename T>
+    struct DenseMatrixConverter {
+        static DenseMatrix from(const Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic>& other)
+        {
+            DenseMatrix mat(other.rows(), other.cols());
+
+            for (size_t j = 0; j < mat.Width(); j++)
+                for (size_t i = 0; i < mat.Height(); i++)
+                    mat(i, j) = other(i, j);
+
+            return mat;
+        }
+
+        static Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> to(const DenseMatrix& other)
+        {
+            Eigen::Matrix<T, Eigen::Dynamic, Eigen::Dynamic> mat(other.Height(), other.Width());
+
+            for (size_t j = 0; j < mat.cols(); j++)
+                for (size_t i = 0; i < mat.rows(); i++)
+                    mat(i, j) = other(i, j);
+
+            return mat;
+        }
+    };
+
+    /** @brief Eigen template specialization for sparse matrix conversion */
+    template <class T>
+    struct SparseMatrixConverter {
+        static SparseMatrix from(const Eigen::SparseMatrix<T, Eigen::RowMajor>& other)
+        {
+            return SparseMatrix(other.outerIndexPtr(), other.innerIndexPtr(), other.valuePtr(), other.rows(), other.cols());
+        }
+
+        static Eigen::SparseMatrix<T, Eigen::RowMajor> to(const SparseMatrix& other)
+        {
+            // MFEM memory info
+            const int *I = other.GetI(), *J = other.GetJ();
+            const T* Data = other.GetData();
+
+            // Eigen triplet
+            std::vector<Eigen::Triplet<double>> tripletList;
+            tripletList.reserve(other.GetMemoryData().Capacity());
+
+            for (size_t i = 0; i < other.Size(); i++) {
+                for (size_t k = I[i], end = I[i + 1]; k < end; k++)
+                    tripletList.push_back(Eigen::Triplet<double>(i, J[k], Data[k]));
+            }
+
+            // Create Eigen sparse matrix
+            Eigen::SparseMatrix<T, Eigen::RowMajor> mat(other.Height(), other.Width());
+            mat.setFromTriplets(tripletList.begin(), tripletList.end());
+
+            return mat;
+        }
+    };   
+}
+
+#endif // MFEM_EIGEN_HPP

linalg/eigensolver.cpp

@@ -0,0 +1,23 @@
+// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "linalg.hpp"
+#include "eigensolver.hpp"
+
+using namespace std;
+
+namespace mfem
+{
+
+Eigensolver::Eigensolver()
+{}
+
+};

linalg/eigensolver.hpp

@@ -0,0 +1,53 @@
+// Copyright (c) 2010-2020, 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_EIGENSOLVERS
+#define MFEM_EIGENSOLVERS
+
+#include "vector.hpp"
+#include "operator.hpp"
+
+namespace mfem
+{
+
+/// Abstract Eigensolver
+class Eigensolver
+{
+public:
+
+   Eigensolver();
+   virtual ~Eigensolver() {}
+
+   virtual void SetTol(double tol) = 0;
+   virtual void SetMaxIter(int max_iter) = 0;
+   virtual void SetPrintLevel(int logging) = 0;
+   virtual void SetNumModes(int num_eigs) = 0;
+
+   virtual void SetOperator(Operator & A) = 0;
+   virtual void SetMassMatrix(Operator & M) = 0;
+
+   /// Perform the eigenvalue solve
+   virtual void Solve() = 0;
+
+   /// Collect the converged eigenvalues
+   virtual void GetEigenvalues(Array<double> & eigenvalues) = 0;
+
+   /// Extract a single eigenvector
+   virtual Vector & GetEigenvector(unsigned int i) = 0;
+
+   /// Transfer ownership of the converged eigenvectors
+   virtual Vector ** StealEigenvectors() = 0;
+
+};
+
+}
+
+#endif

linalg/linalg.hpp

@@ -48,6 +48,16 @@
 #include "ginkgo.hpp"
 #endif
 
+#ifdef MFEM_USE_ARPACK
+#include "eigensolver.hpp"
+#include "arpack.hpp"
+#endif
+
+#ifdef MFEM_USE_SPECTRA
+#include "eigen.hpp"
+#include "spectra.hpp"
+#endif
+
 #ifdef MFEM_USE_MPI
 #include "hypre_parcsr.hpp"
 #include "hypre.hpp"

linalg/spectra.cpp

@@ -0,0 +1,157 @@
+#include "spectra.hpp"
+
+#include "../fem/bilinearform.hpp"
+
+namespace mfem {
+    SpectraEigenSolver::SpectraEigenSolver()
+    {
+        // Init params
+        _nconv = 0;
+        _nev = 1;
+        _ncv = 1;
+        _max_iter = 1000;
+        _tol = 1e-3;
+    }
+
+    SpectraEigenSolver::~SpectraEigenSolver()
+    {
+        delete _A_s, _B_s, _S, _G;
+    }
+
+    /// Set dimension of Krylov subspace in the Lanczos method
+    SpectraEigenSolver& SpectraEigenSolver::SetKrylov(double ncv)
+    {
+        _ncv = ncv;
+
+        return *this;
+    }
+
+    /// Set solver tolerance
+    SpectraEigenSolver& SpectraEigenSolver::SetTol(double tol)
+    {
+        _tol = tol;
+
+        return *this;
+    }
+
+    /// Set maximum number of iterations
+    SpectraEigenSolver& SpectraEigenSolver::SetMaxIter(int max_iter)
+    {
+        _max_iter = max_iter;
+
+        return *this;
+    }
+
+    /// Set the number of required eigenmodes
+    SpectraEigenSolver& SpectraEigenSolver::SetNumModes(int nev)
+    {
+        _nev = nev;
+
+        return *this;
+    }
+
+    /// Set operator for standard eigenvalue problem (A*x = lambda*x)
+    SpectraEigenSolver& SpectraEigenSolver::SetOperator(const Operator& A)
+    {
+        // Set EIGEN operators
+        _A_e = SparseMatrixConverter<double>::to(static_cast<const BilinearForm&>(A).SpMat());
+
+        // Set SPECTRA operators
+        _A_s = new SparseSymMatProd<double>(_A_e);
+
+        return *this;
+    }
+
+    /// Set operator for generalized eigenvalue problem (A*x = lambda*B*x)
+    SpectraEigenSolver& SpectraEigenSolver::SetOperators(const Operator& A, const Operator& B)
+    {
+        // Set EIGEN operators
+        _A_e = SparseMatrixConverter<double>::to(static_cast<const BilinearForm&>(A).SpMat());
+        _B_e = SparseMatrixConverter<double>::to(static_cast<const BilinearForm&>(B).SpMat());
+
+        // Set SPECTRA operators
+        _A_s = new SparseSymMatProd<double>(_A_e);
+        _B_s = new SparseCholesky<double>(_B_e);
+
+        return *this;
+    }
+
+    /// Solve the eigenvalue problem for the specified number of eigenvalues
+    void SpectraEigenSolver::Solve()
+    {
+        // Set the dimension of the Krilov space equal to the number of requested eigenvalues if necessary
+        if (_ncv < _nev)
+            _ncv = _nev;
+
+        if (!_B_s) {
+            _S = new SymEigsSolver<SparseSymMatProd<double>>(*_A_s, _nev, _ncv);
+            _S->init();
+            _nconv = _S->compute(SortRule::SmallestMagn, _max_iter, _tol, SortRule::SmallestMagn);
+        }
+        else {
+            _G = new SymGEigsSolver<SparseSymMatProd<double>, SparseCholesky<double>, GEigsMode::Cholesky>(*_A_s, *_B_s, _nev, _ncv);
+            _G->init();
+            _nconv = _G->compute(SortRule::SmallestMagn, _max_iter, _tol, SortRule::SmallestMagn);
+        }
+    }
+
+    /// Get the number of converged eigenvalues
+    int SpectraEigenSolver::GetNumConverged()
+    {
+        return _nconv;
+    }
+
+    /// Get the corresponding eigenvalue
+    double SpectraEigenSolver::GetEigenvalue(unsigned int i) const
+    {
+        if (!_B_s) {
+            if (_S->info() == CompInfo::Successful && i < _nconv)
+                return _S->eigenvalues()[i];
+            else
+                return 0;
+        }
+        else {
+            if (_G->info() == CompInfo::Successful && i < _nconv)
+                return _G->eigenvalues()[i];
+            else
+                return 0;
+        }
+    }
+
+    Eigen::VectorXd SpectraEigenSolver::GetEigenvalues(unsigned int i) const
+    {
+        if (!_B_s) {
+            if (_S->info() == CompInfo::Successful && i < _nconv)
+                return _S->eigenvalues().segment(0, i);
+        }
+        else {
+            if (_G->info() == CompInfo::Successful && i < _nconv)
+                return _G->eigenvalues().segment(0, i);
+        }
+    }
+
+    /// Get the corresponding eigenvector
+    Eigen::VectorXd SpectraEigenSolver::GetEigenvector(unsigned int i) const
+    {
+        if (!_B_s) {
+            if (_S->info() == CompInfo::Successful && i < _nconv)
+                return _S->eigenvectors().col(i);
+        }
+        else {
+            if (_G->info() == CompInfo::Successful && i < _nconv)
+                return _G->eigenvectors().col(i);
+        }
+    }
+
+    Eigen::MatrixXd SpectraEigenSolver::GetEigenvectors(unsigned int i) const
+    {
+        if (!_B_s) {
+            if (_S->info() == CompInfo::Successful && i < _nconv)
+                return _S->eigenvectors().topRows(i);
+        }
+        else {
+            if (_G->info() == CompInfo::Successful && i < _nconv)
+                return _G->eigenvectors().topRows(i);
+        }
+    }
+} // namespace mfem

linalg/spectra.hpp

@@ -0,0 +1,77 @@
+#ifndef MFEM_SPECTRA_HPP
+#define MFEM_SPECTRA_HPP
+
+#include <Spectra/GenEigsSolver.h>
+#include <Spectra/MatOp/SparseCholesky.h>
+#include <Spectra/MatOp/SparseGenMatProd.h>
+#include <Spectra/SymEigsSolver.h>
+#include <Spectra/SymGEigsSolver.h>
+
+#include "eigen.hpp"
+
+using namespace Spectra;
+
+namespace mfem {
+    class SpectraEigenSolver {
+    public:
+        SpectraEigenSolver();
+
+        virtual ~SpectraEigenSolver();
+
+        /// Set dimension of Krylov subspace in the Lanczos method
+        SpectraEigenSolver& SetKrylov(double ncv);
+
+        /// Set solver tolerance
+        SpectraEigenSolver& SetTol(double tol);
+
+        /// Set maximum number of iterations
+        SpectraEigenSolver& SetMaxIter(int max_iter);
+
+        /// Set the number of required eigenmodes
+        SpectraEigenSolver& SetNumModes(int nev);
+
+        /// Set operator for standard eigenvalue problem (A*x = lambda*x)
+        SpectraEigenSolver& SetOperator(const Operator& A);
+
+        /// Set operator for generalized eigenvalue problem (A*x = lambda*B*x)
+        SpectraEigenSolver& SetOperators(const Operator& A, const Operator& B);
+
+        /// Solve the eigenvalue problem for the specified number of eigenvalues
+        void Solve();
+
+        /// Get the number of converged eigenvalues
+        int GetNumConverged();
+
+        /// Get the corresponding eigenvalue
+        double GetEigenvalue(unsigned int i) const;
+
+        Eigen::VectorXd GetEigenvalues(unsigned int i = 0) const;
+
+        /// Get the corresponding eigenvector
+        Eigen::VectorXd GetEigenvector(unsigned int i) const;
+
+        Eigen::MatrixXd GetEigenvectors(unsigned int i) const;
+
+    protected:
+        // Params
+        int _nconv, _nev, _ncv, _max_iter;
+        double _tol;
+
+        // EIGEN Operators
+        Eigen::SparseMatrix<double> _A_e, _B_e;
+
+        // Spectra Operators
+        SparseSymMatProd<double>* _A_s = nullptr;
+        SparseCholesky<double>* _B_s = nullptr;
+
+        // Eigenvalue solution based on Spectra
+        SymEigsSolver<SparseSymMatProd<double>>* _S = nullptr;
+        SymGEigsSolver<SparseSymMatProd<double>, SparseCholesky<double>, GEigsMode::Cholesky>* _G = nullptr;
+
+        // // Eigenvalue solution based on Eigen
+        // Eigen::SelfAdjointEigenSolver<Eigen::MatrixXd>* _S = nullptr;
+        // Eigen::GeneralizedSelfAdjointEigenSolver<Eigen::MatrixXd>* _G = nullptr;
+    };
+} // namespace mfem
+
+#endif // MFEM_SPECTRA_HPP

makefile

@@ -274,7 +274,7 @@ endif
 # List of MFEM dependencies, that require the *_LIB variable to be non-empty
 MFEM_REQ_LIB_DEPS = SUPERLU MUMPS METIS FMS CONDUIT SIDRE LAPACK SUNDIALS MESQUITE\
  SUITESPARSE STRUMPACK GINKGO GNUTLS NETCDF PETSC SLEPC MPFR PUMI HIOP GSLIB\
- OCCA CEED RAJA UMPIRE MKL_CPARDISO AMGX CALIPER
+ OCCA CEED RAJA UMPIRE MKL_CPARDISO AMGX CALIPER ARPACK
 
 PETSC_ERROR_MSG = $(if $(PETSC_FOUND),,. PETSC config not found: $(PETSC_VARS))
 SLEPC_ERROR_MSG = $(if $(SLEPC_FOUND),,. SLEPC config not found: $(SLEPC_VARS))
@@ -292,7 +292,7 @@ ifeq ($(MAKECMDGOALS),config)
 endif
 
 # List of MFEM dependencies, processed below
-MFEM_DEPENDENCIES = $(MFEM_REQ_LIB_DEPS) LIBUNWIND OPENMP CUDA HIP
+MFEM_DEPENDENCIES = $(MFEM_REQ_LIB_DEPS) SPECTRA LIBUNWIND OPENMP CUDA HIP
 
 # List of deprecated MFEM dependencies, processed below
 MFEM_LEGACY_DEPENDENCIES = OPENMP
@@ -340,7 +340,7 @@ MFEM_DEFINES = MFEM_VERSION MFEM_VERSION_STRING MFEM_GIT_STRING MFEM_USE_MPI\
  MFEM_USE_PUMI MFEM_USE_HIOP MFEM_USE_GSLIB MFEM_USE_CUDA MFEM_USE_HIP\
  MFEM_USE_OCCA MFEM_USE_CEED MFEM_USE_RAJA MFEM_USE_UMPIRE MFEM_USE_SIMD\
  MFEM_USE_ADIOS2 MFEM_USE_MKL_CPARDISO MFEM_USE_AMGX MFEM_USE_MUMPS\
- MFEM_USE_CALIPER MFEM_SOURCE_DIR MFEM_INSTALL_DIR
+ MFEM_USE_CALIPER MFEM_USE_ARPACK MFEM_USE_SPECTRA MFEM_SOURCE_DIR MFEM_INSTALL_DIR
 
 # List of makefile variables that will be written to config.mk:
 MFEM_CONFIG_VARS = MFEM_CXX MFEM_HOST_CXX MFEM_CPPFLAGS MFEM_CXXFLAGS\
@@ -652,6 +652,8 @@ status info:
 	$(info MFEM_USE_SUNDIALS      = $(MFEM_USE_SUNDIALS))
 	$(info MFEM_USE_MESQUITE      = $(MFEM_USE_MESQUITE))
 	$(info MFEM_USE_SUITESPARSE   = $(MFEM_USE_SUITESPARSE))
+	$(info MFEM_USE_ARPACK        = $(MFEM_USE_ARPACK))
+	$(info MFEM_USE_SPECTRA       = $(MFEM_USE_SPECTRA))
 	$(info MFEM_USE_SUPERLU       = $(MFEM_USE_SUPERLU))
 	$(info MFEM_USE_MUMPS         = $(MFEM_USE_MUMPS))
 	$(info MFEM_USE_STRUMPACK     = $(MFEM_USE_STRUMPACK))

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)
@@ -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);
@@ -3428,6 +3555,7 @@ Element *Mesh::NewElement(int geom)
 #endif
       case Geometry::CUBE:      return (new Hexahedron);
       case Geometry::PRISM:     return (new Wedge);
+      case Geometry::PYRAMID:   return (new Pyramid);
       default:
          MFEM_ABORT("invalid Geometry::Type, geom = " << geom);
    }
@@ -3520,6 +3648,15 @@ void Mesh::SetMeshGen()
             meshgen |= 4;
             break;
 
+         case Element::PYRAMID:
+            mesh_geoms |= (1 << Geometry::PYRAMID);
+            mesh_geoms |= (1 << Geometry::SQUARE);
+            mesh_geoms |= (1 << Geometry::TRIANGLE);
+            mesh_geoms |= (1 << Geometry::SEGMENT);
+            mesh_geoms |= (1 << Geometry::POINT);
+            meshgen |= 8;
+            break;
+
          default:
             MFEM_ABORT("invalid element type: " << type);
             break;
@@ -5060,6 +5197,19 @@ int Mesh::CheckElementOrientation(bool fix_it)
                }
                break;
 
+            case Element::PYRAMID:
+               // only check the Jacobian at the center of the element
+               GetElementJacobian(i, J);
+               if (J.Det() < 0.0)
+               {
+                  wo++;
+                  if (fix_it)
+                  {
+                     // how?
+                  }
+               }
+               break;
+
             case Element::HEXAHEDRON:
                // only check the Jacobian at the center of the element
                GetElementJacobian(i, J);
@@ -6142,6 +6292,22 @@ void Mesh::GenerateFaces()
                }
                break;
             }
+            case Element::PYRAMID:
+            {
+               for (int j = 0; j < 1; j++)
+               {
+                  const int *fv = pyr_t::FaceVert[j];
+                  AddQuadFaceElement(j, ef[j], i,
+                                     v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
+               }
+               for (int j = 1; j < 5; j++)
+               {
+                  const int *fv = pyr_t::FaceVert[j];
+                  AddTriangleFaceElement(j, ef[j], i,
+                                         v[fv[0]], v[fv[1]], v[fv[2]]);
+               }
+               break;
+            }
             case Element::HEXAHEDRON:
             {
                for (int j = 0; j < 6; j++)
@@ -6234,6 +6400,20 @@ STable3D *Mesh::GetFacesTable()
             }
             break;
          }
+         case Element::PYRAMID:
+         {
+            for (int j = 0; j < 1; j++)
+            {
+               const int *fv = pyr_t::FaceVert[j];
+               faces_tbl->Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
+            }
+            for (int j = 1; j < 5; j++)
+            {
+               const int *fv = pyr_t::FaceVert[j];
+               faces_tbl->Push(v[fv[0]], v[fv[1]], v[fv[2]]);
+            }
+            break;
+         }
          case Element::WEDGE:
          {
             for (int j = 0; j < 2; j++)
@@ -6308,6 +6488,22 @@ STable3D *Mesh::GetElementToFaceTable(int ret_ftbl)
             }
             break;
          }
+         case Element::PYRAMID:
+         {
+            for (int j = 0; j < 1; j++)
+            {
+               const int *fv = pyr_t::FaceVert[j];
+               el_to_face->Push(
+                  i, faces_tbl->Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]));
+            }
+            for (int j = 1; j < 5; j++)
+            {
+               const int *fv = pyr_t::FaceVert[j];
+               el_to_face->Push(
+                  i, faces_tbl->Push(v[fv[0]], v[fv[1]], v[fv[2]]));
+            }
+            break;
+         }
          case Element::HEXAHEDRON:
          {
             // find the face by the vertices with the smallest 3 numbers
@@ -7644,6 +7840,18 @@ void Mesh::UniformRefinement3D_base(Array<int> *f2qf_ptr, DSTable *v_to_v_p,
       }
    }
 
+   int pyr_counter = 0;
+   if (HasGeometry(Geometry::PYRAMID))
+   {
+      for (int i = 0; i < elements.Size(); i++)
+      {
+         if (elements[i]->GetType() == Element::PYRAMID)
+         {
+            pyr_counter++;
+         }
+      }
+   }
+
    // Map from edge-index to vertex-index, needed for ReorientTetMesh() for
    // parallel meshes.
    // Note: with the removal of ReorientTetMesh() this may no longer
@@ -7703,7 +7911,7 @@ void Mesh::UniformRefinement3D_base(Array<int> *f2qf_ptr, DSTable *v_to_v_p,
    Array<Element*> new_boundary;
 
    vertices.SetSize(oelem + hex_counter);
-   new_elements.SetSize(8 * NumOfElements);
+   new_elements.SetSize(8 * NumOfElements + 2 * pyr_counter);
    CoarseFineTr.embeddings.SetSize(new_elements.Size());
 
    hex_counter = 0;
@@ -7969,6 +8177,73 @@ void Mesh::UniformRefinement3D_base(Array<int> *f2qf_ptr, DSTable *v_to_v_p,
          }
          break;
 
+         case Element::PYRAMID:
+         {
+            const int *f = el_to_face->GetRow(i);
+            // pyr_counter++;
+
+            for (int fi = 0; fi < 1; fi++)
+            {
+               for (int k = 0; k < 4; k++)
+               {
+                  vv[k] = v[pyr_t::FaceVert[fi][k]];
+               }
+               AverageVertices(vv, 4, oface + f2qf[f[fi]]);
+            }
+
+            for (int ei = 0; ei < 8; ei++)
+            {
+               for (int k = 0; k < 2; k++)
+               {
+                  vv[k] = v[pyr_t::Edges[ei][k]];
+               }
+               AverageVertices(vv, 2, oedge+e[ei]);
+            }
+
+            const int qf0 = f2qf[f[0]];
+
+            new_elements[j++] =
+               new Pyramid(v[0], oedge+e[0], oface+qf0,
+                           oedge+e[3], oedge+e[4], attr);
+
+            new_elements[j++] =
+               new Pyramid(oedge+e[0], v[1], oedge+e[1],
+                           oface+qf0, oedge+e[5], attr);
+
+            new_elements[j++] =
+               new Pyramid(oface+qf0, oedge+e[1], v[2],
+                           oedge+e[2], oedge+e[6], attr);
+
+            new_elements[j++] =
+               new Pyramid(oedge+e[3], oface+qf0, oedge+e[2],
+                           v[3], oedge+e[7], attr);
+
+            new_elements[j++] =
+               new Pyramid(oedge+e[4], oedge+e[5], oedge+e[6],
+                           oedge+e[7], v[4], attr);
+
+            new_elements[j++] =
+               new Pyramid(oedge+e[7], oedge+e[6], oedge+e[5],
+                           oedge+e[4], oface+qf0, attr);
+
+            new_elements[j++] =
+               new Tetrahedron(oedge+e[0], oedge+e[4], oedge+e[5],
+                               oface+qf0, attr);
+
+            new_elements[j++] =
+               new Tetrahedron(oedge+e[1], oedge+e[5], oedge+e[6],
+                               oface+qf0, attr);
+
+            new_elements[j++] =
+               new Tetrahedron(oedge+e[2], oedge+e[6], oedge+e[7],
+                               oface+qf0, attr);
+
+            new_elements[j++] =
+               new Tetrahedron(oedge+e[3], oedge+e[7], oedge+e[4],
+                               oface+qf0, attr);
+         }
+         break;
+
          case Element::HEXAHEDRON:
          {
             const int *f = el_to_face->GetRow(i);
@@ -8100,7 +8375,7 @@ void Mesh::UniformRefinement3D_base(Array<int> *f2qf_ptr, DSTable *v_to_v_p,
    }
    mfem::Swap(boundary, new_boundary);
 
-   static const double A = 0.0, B = 0.5, C = 1.0;
+   static const double A = 0.0, B = 0.5, C = 1.0, D = -1.0;
    static double tet_children[3*4*16] =
    {
       A,A,A, B,A,A, A,B,A, A,A,B,
@@ -8126,6 +8401,19 @@ void Mesh::UniformRefinement3D_base(Array<int> *f2qf_ptr, DSTable *v_to_v_p,
       A,A,B, A,B,B, B,A,B, B,B,A,
       A,A,B, B,A,B, B,A,A, B,B,A
    };
+   static double pyr_children[3*5*10] =
+   {
+      A,A,A, B,A,A, B,B,A, A,B,A, A,A,B,
+      B,A,A, C,A,A, C,B,A, B,B,A, B,A,B,
+      B,B,A, C,B,A, C,C,A, B,C,A, B,B,B,
+      A,B,A, B,B,A, B,C,A, A,C,A, A,B,B,
+      A,A,B, B,A,B, B,B,B, A,B,B, A,A,C,
+      A,B,B, B,B,B, B,A,B, A,A,B, B,B,A,
+      B,A,A, A,A,B, B,A,B, B,B,A, D,D,D,
+      C,B,A, B,A,B, B,B,B, B,B,A, D,D,D,
+      B,C,A, B,B,B, A,B,B, B,B,A, D,D,D,
+      A,B,A, A,B,B, A,A,B, B,B,A, D,D,D
+   };
    static double pri_children[3*6*8] =
    {
       A,A,A, B,A,A, A,B,A, A,A,B, B,A,B, A,B,B,
@@ -8151,6 +8439,8 @@ void Mesh::UniformRefinement3D_base(Array<int> *f2qf_ptr, DSTable *v_to_v_p,
 
    CoarseFineTr.point_matrices[Geometry::TETRAHEDRON]
    .UseExternalData(tet_children, 3, 4, 16);
+   CoarseFineTr.point_matrices[Geometry::PYRAMID]
+   .UseExternalData(pyr_children, 3, 5, 10);
    CoarseFineTr.point_matrices[Geometry::PRISM]
    .UseExternalData(pri_children, 3, 6, 8);
    CoarseFineTr.point_matrices[Geometry::CUBE]
@@ -8167,7 +8457,7 @@ void Mesh::UniformRefinement3D_base(Array<int> *f2qf_ptr, DSTable *v_to_v_p,
    }
 
    NumOfVertices    = vertices.Size();
-   NumOfElements    = 8 * NumOfElements;
+   NumOfElements    = 8 * NumOfElements + 2 * pyr_counter;
    NumOfBdrElements = 4 * NumOfBdrElements;
 
    GetElementToFaceTable();
@@ -8751,7 +9041,7 @@ void Mesh::GeneralRefinement(const Array<Refinement> &refinements,
    else if (nonconforming < 0)
    {
       // determine if nonconforming refinement is suitable
-      if ((meshgen & 2) || (meshgen & 4))
+      if ((meshgen & 2) || (meshgen & 4) || (meshgen & 8))
       {
          nonconforming = 1; // tensor product elements and wedges
       }
@@ -9535,6 +9825,7 @@ void Mesh::Printer(std::ostream &out, std::string section_delimiter) const
        "# TETRAHEDRON = 4\n"
        "# CUBE        = 5\n"
        "# PRISM       = 6\n"
+       "# PYRAMID     = 7\n"
        "#\n";
 
    out << "\ndimension\n" << Dim;

mesh/mesh.hpp

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

mesh/mesh_headers.hpp

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

mesh/mesh_readers.cpp

@@ -352,6 +352,11 @@ void Mesh::ReadTrueGridMesh(std::istream &input)
 const int Mesh::vtk_quadratic_tet[10] =
 { 0, 1, 2, 3, 4, 7, 5, 6, 8, 9 };
 
+// see Pyramid::edges & Mesh::GenerateFaces
+// https://www.vtk.org/doc/nightly/html/classvtkBiQuadraticQuadraticWedge.html
+const int Mesh::vtk_quadratic_pyramid[13] =
+{ 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12};
+
 // see Wedge::edges & Mesh::GenerateFaces
 // https://www.vtk.org/doc/nightly/html/classvtkBiQuadraticQuadraticWedge.html
 const int Mesh::vtk_quadratic_wedge[18] =
@@ -545,6 +550,8 @@ void Mesh::CreateVTKMesh(const Vector &points, const Array<int> &cell_data,
                   vtk_mfem = vtk_quadratic_hex; break;
                case Geometry::PRISM:
                   vtk_mfem = vtk_quadratic_wedge; break;
+               case Geometry::PYRAMID:
+                  vtk_mfem = vtk_quadratic_pyramid; break;
                default:
                   vtk_mfem = NULL; // suppress a warning
                   break;
@@ -1393,6 +1400,10 @@ void Mesh::ReadInlineMesh(std::istream &input, bool generate_edges)
          {
             type = Element::WEDGE;
          }
+         else if (eltype == "pyramid")
+         {
+            type = Element::PYRAMID;
+         }
          else if (eltype == "tet")
          {
             type = Element::TETRAHEDRON;
@@ -1447,7 +1458,7 @@ void Mesh::ReadInlineMesh(std::istream &input, bool generate_edges)
       Make2D(nx, ny, type, sx, sy, generate_edges, true);
    }
    else if (type == Element::TETRAHEDRON || type == Element::WEDGE ||
-            type == Element::HEXAHEDRON)
+            type == Element::HEXAHEDRON  || type == Element::PYRAMID)
    {
       MFEM_VERIFY(nx > 0 && ny > 0 && nz > 0 &&
                   sx > 0.0 && sy > 0.0 && sz > 0.0,
@@ -1884,6 +1895,9 @@ void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
          ho_wdg[2] = wdg18; ho_wdg[3] = wdg40;
          ho_pyr[2] = pyr14; ho_pyr[3] = pyr30;
 
+         bool has_nonpositive_phys_domain = false;
+         bool has_positive_phys_domain = false;
+
          if (binary)
          {
             int n_elem_part = 0; // partial sum of elements that are read
@@ -1934,17 +1948,19 @@ void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
                      vert_indices[vi] = it->second;
                   }
 
-                  // non-positive attributes are not allowed in MFEM
+                  // Non-positive attributes are not allowed in MFEM. However,
+                  // by default, Gmsh sets the physical domain of all elements
+                  // to zero. In the case that all elements have physical domain
+                  // zero, we will given them attribute 1. If only some elements
+                  // have physical domain zero, we will throw an error.
                   if (phys_domain <= 0)
                   {
-                     MFEM_ABORT("Non-positive element attribute in Gmsh mesh!\n"
-                                "By default Gmsh sets element tags (attributes)"
-                                " to '0' but MFEM requires that they be"
-                                " positive integers.\n"
-                                "Use \"Physical Curve\", \"Physical Surface\","
-                                " or \"Physical Volume\" to set tags/attributes"
-                                " for all curves, surfaces, or volumes in your"
-                                " Gmsh geometry to values which are >= 1.");
+                     has_nonpositive_phys_domain = true;
+                     phys_domain = 1;
+                  }
+                  else
+                  {
+                     has_positive_phys_domain = true;
                   }
 
                   // initialize the mesh element
@@ -2161,17 +2177,19 @@ void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
                   vert_indices[vi] = it->second;
                }
 
-               // non-positive attributes are not allowed in MFEM
+               // Non-positive attributes are not allowed in MFEM. However,
+               // by default, Gmsh sets the physical domain of all elements
+               // to zero. In the case that all elements have physical domain
+               // zero, we will given them attribute 1. If only some elements
+               // have physical domain zero, we will throw an error.
                if (phys_domain <= 0)
                {
-                  MFEM_ABORT("Non-positive element attribute in Gmsh mesh!\n"
-                             "By default Gmsh sets element tags (attributes)"
-                             " to '0' but MFEM requires that they be"
-                             " positive integers.\n"
-                             "Use \"Physical Curve\", \"Physical Surface\","
-                             " or \"Physical Volume\" to set tags/attributes"
-                             " for all curves, surfaces, or volumes in your"
-                             " Gmsh geometry to values which are >= 1.");
+                  has_nonpositive_phys_domain = true;
+                  phys_domain = 1;
+               }
+               else
+               {
+                  has_positive_phys_domain = true;
                }
 
                // initialize the mesh element
@@ -2356,6 +2374,24 @@ void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
             } // el (all elements)
          } // if ASCII
 
+         if (has_positive_phys_domain && has_nonpositive_phys_domain)
+         {
+            MFEM_ABORT("Non-positive element attribute in Gmsh mesh!\n"
+                       "By default Gmsh sets element tags (attributes)"
+                       " to '0' but MFEM requires that they be"
+                       " positive integers.\n"
+                       "Use \"Physical Curve\", \"Physical Surface\","
+                       " or \"Physical Volume\" to set tags/attributes"
+                       " for all curves, surfaces, or volumes in your"
+                       " Gmsh geometry to values which are >= 1.");
+         }
+         else if (has_nonpositive_phys_domain)
+         {
+            mfem::out << "\nGmsh reader: all element attributes were zero.\n"
+                      << "MFEM only supports positive element attributes.\n"
+                      << "Setting element attributes to 1.\n\n";
+         }
+
          if (!elements_3D.empty())
          {
             Dim = 3;

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