add PointMatrix for Geometry::SEGMENT

.gitignore

@@ -145,14 +145,6 @@ 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
@@ -318,13 +310,7 @@ tests/convergence/prates
 tests/par-mesh-format/ex1p
 
 # VPATH builds
-build-*/
-
-# User config
-user-*
-
-# VSCode
-.vscode 
+build-*/*
 
 # PETSc automated build
 petsc-build/*

CHANGELOG

@@ -12,13 +12,6 @@ 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.
 
@@ -27,9 +20,6 @@ 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,12 +91,6 @@
 // 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,8 +31,6 @@ 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,8 +151,6 @@ 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
@@ -330,19 +328,6 @@ 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

@@ -1,9 +0,0 @@
-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

@@ -1,43 +0,0 @@
-MFEM mesh v1.0
-
-#
-# MFEM Geometry Types (see mesh/geom.hpp):
-#
-# POINT       = 0
-# SEGMENT     = 1
-# TRIANGLE    = 2
-# SQUARE      = 3
-# TETRAHEDRON = 4
-# CUBE        = 5
-# 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

@@ -1,38 +0,0 @@
-MFEM mesh v1.0
-
-#
-# MFEM Geometry Types (see mesh/geom.hpp):
-#
-# POINT       = 0
-# SEGMENT     = 1
-# TRIANGLE    = 2
-# SQUARE      = 3
-# TETRAHEDRON = 4
-# CUBE        = 5
-# 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

@@ -1,47 +0,0 @@
-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

@@ -1,286 +0,0 @@
-//                       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

@@ -1,69 +0,0 @@
-# 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,7 +9,6 @@
 //               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,7 +9,6 @@
 //               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,8 +13,6 @@
 //               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,8 +13,6 @@
 //               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,8 +16,6 @@
 //               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,8 +16,6 @@
 //               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,8 +19,6 @@
 //               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,8 +19,6 @@
 //               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

@@ -1,250 +0,0 @@
-//                       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

@@ -1,67 +0,0 @@
-# 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,7 +380,6 @@ 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,64 +1313,6 @@ 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
 {
@@ -1748,32 +1690,6 @@ 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
@@ -1912,10 +1828,6 @@ 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;
 };
 
 
@@ -1940,66 +1852,6 @@ 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;
 };
 
 
@@ -2093,77 +1945,6 @@ 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,9 +33,6 @@ 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");
    }
@@ -577,7 +574,6 @@ 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.");
    }
@@ -595,7 +591,6 @@ 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.");
    }
@@ -1245,7 +1240,6 @@ 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.");
    }
@@ -1263,7 +1257,6 @@ 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.");
    }
@@ -1284,8 +1277,6 @@ 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.");
@@ -1302,8 +1293,6 @@ 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.");
@@ -1405,8 +1394,6 @@ 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.");
    }
@@ -1423,8 +1410,6 @@ 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.");
    }
@@ -1454,8 +1439,6 @@ 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.");
    }
@@ -1472,8 +1455,6 @@ 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.");
    }
@@ -1749,7 +1730,6 @@ 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);
@@ -1763,7 +1743,6 @@ 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];
@@ -1858,21 +1837,6 @@ 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
 {
@@ -2112,12 +2076,9 @@ 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)
       {
@@ -2238,21 +2199,6 @@ 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
 {
@@ -2344,12 +2290,6 @@ 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
    {
@@ -2493,22 +2433,6 @@ 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
 {
@@ -2769,28 +2693,6 @@ 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,7 +228,8 @@ public:
                             const int btype = BasisType::GaussLobatto);
 
    virtual const FiniteElement *FiniteElementForGeometry(
-      Geometry::Type GeomType) const;
+      Geometry::Type GeomType) const
+   { return H1_Elements[GeomType]; }
    virtual int DofForGeometry(Geometry::Type GeomType) const
    { return H1_dof[GeomType]; }
    virtual const int *DofOrderForOrientation(Geometry::Type GeomType,
@@ -301,7 +302,10 @@ public:
                    const int map_type = FiniteElement::VALUE);
 
    virtual const FiniteElement *FiniteElementForGeometry(
-      Geometry::Type GeomType) const;
+      Geometry::Type GeomType) const
+   {
+      return L2_Elements[GeomType];
+   }
    virtual int DofForGeometry(Geometry::Type GeomType) const
    {
       if (L2_Elements[GeomType])
@@ -367,7 +371,8 @@ public:
                    const int ob_type = BasisType::GaussLegendre);
 
    virtual const FiniteElement *FiniteElementForGeometry(
-      Geometry::Type GeomType) const;
+      Geometry::Type GeomType) const
+   { return RT_Elements[GeomType]; }
    virtual int DofForGeometry(Geometry::Type GeomType) const
    { return RT_dof[GeomType]; }
    virtual const int *DofOrderForOrientation(Geometry::Type GeomType,
@@ -425,7 +430,8 @@ public:
                    const int ob_type = BasisType::GaussLegendre);
 
    virtual const FiniteElement *
-   FiniteElementForGeometry(Geometry::Type GeomType) const;
+   FiniteElementForGeometry(Geometry::Type GeomType) const
+   { return ND_Elements[GeomType]; }
 
    virtual int DofForGeometry(Geometry::Type GeomType) const
    { return ND_dof[GeomType]; }
@@ -523,10 +529,9 @@ private:
    const BiLinear2DFiniteElement QuadrilateralFE;
    const Linear3DFiniteElement TetrahedronFE;
    const TriLinear3DFiniteElement ParallelepipedFE;
-   const LinearWedgeFiniteElement WedgeFE;
-   const LinearPyramidFiniteElement PyramidFE;
+   const H1_WedgeElement WedgeFE;
 public:
-   LinearFECollection() : FiniteElementCollection(1) { }
+   LinearFECollection() : FiniteElementCollection(1), WedgeFE(1) { }
 
    virtual const FiniteElement *
    FiniteElementForGeometry(Geometry::Type GeomType) const;
@@ -931,11 +936,10 @@ class Const3DFECollection : public FiniteElementCollection
 private:
    const P0TetFiniteElement TetrahedronFE;
    const P0HexFiniteElement ParallelepipedFE;
-   const P0WdgFiniteElement WedgeFE;
-   const P0PyrFiniteElement PyramidFE;
+   const L2_WedgeElement WedgeFE;
 
 public:
-   Const3DFECollection() : FiniteElementCollection(0) { }
+   Const3DFECollection() : FiniteElementCollection(0), WedgeFE(0) { }
 
    virtual const FiniteElement *
    FiniteElementForGeometry(Geometry::Type GeomType) const;
@@ -956,8 +960,6 @@ class LinearDiscont3DFECollection : public FiniteElementCollection
 {
 private:
    const Linear3DFiniteElement TetrahedronFE;
-   const LinearPyramidFiniteElement PyramidFE;
-   const LinearWedgeFiniteElement WedgeFE;
    const TriLinear3DFiniteElement ParallelepipedFE;
 
 public:
@@ -1034,8 +1036,6 @@ class ND1_3DFECollection : public FiniteElementCollection
 private:
    const Nedelec1HexFiniteElement HexahedronFE;
    const Nedelec1TetFiniteElement TetrahedronFE;
-   const Nedelec1WdgFiniteElement WedgeFE;
-   const Nedelec1PyrFiniteElement PyramidFE;
 
 public:
    ND1_3DFECollection() : FiniteElementCollection(1) { }
@@ -1061,8 +1061,6 @@ 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,19 +11,15 @@
 
 #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",
-   "Pyramid"
-};
+{ "Point", "Segment", "Triangle", "Square", "Tetrahedron", "Cube", "Prism" };
 
 const double Geometry::Volume[NumGeom] =
-{ 1.0, 1.0, 0.5, 1.0, 1./6, 1.0, 0.5, 1./3 };
+{ 1.0, 1.0, 0.5, 1.0, 1./6, 1.0, 0.5 };
 
 Geometry::Geometry()
 {
@@ -143,28 +139,6 @@ 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;
@@ -193,10 +167,6 @@ 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);
@@ -204,7 +174,6 @@ 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;
@@ -213,7 +182,6 @@ 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);
    {
@@ -242,14 +210,6 @@ 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()
@@ -273,7 +233,6 @@ 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(...)");
    }
@@ -351,25 +310,6 @@ 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!");
    }
@@ -431,10 +371,6 @@ 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!");
    }
@@ -505,17 +441,6 @@ 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!");
    }
@@ -630,16 +555,6 @@ 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!");
    }
@@ -737,43 +652,6 @@ 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!");
    }
@@ -848,17 +726,6 @@ 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 (...)");
    }
@@ -877,13 +744,13 @@ void Geometry::JacToPerfJac(int GeomType, const DenseMatrix &J,
    }
 }
 
-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::NumBdrArray[NumGeom] = { 0, 2, 3, 4, 4, 6, 5 };
+const int Geometry::Dimension[NumGeom] = { 0, 1, 2, 2, 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, 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::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::
 Constants<Geometry::POINT>::Orient[1][1] = {{0}};
@@ -1030,30 +897,6 @@ 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()
 {
@@ -1419,104 +1262,6 @@ 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,7 +27,6 @@ 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
 {
@@ -35,7 +34,7 @@ public:
    enum Type
    {
       INVALID = -1,
-      POINT = 0, SEGMENT, TRIANGLE, SQUARE, TETRAHEDRON, CUBE, PRISM, PYRAMID,
+      POINT = 0, SEGMENT, TRIANGLE, SQUARE, TETRAHEDRON, CUBE, PRISM,
       NUM_GEOMETRIES
    };
 
@@ -252,26 +251,7 @@ template <> struct Geometry::Constants<Geometry::PRISM>
    };
 };
 
-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`.
+// Defined in fe.cpp to ensure construction after 'mfem::WedgeFE'.
 extern Geometry Geometries;
 
 

fem/intrules.cpp

@@ -910,9 +910,6 @@ 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;
 
@@ -933,7 +930,6 @@ 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;
@@ -980,7 +976,6 @@ 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;
@@ -1024,7 +1019,6 @@ IntegrationRules::~IntegrationRules()
    DeleteIntRuleArray(TetrahedronIntRules);
    DeleteIntRuleArray(CubeIntRules);
    DeleteIntRuleArray(PrismIntRules);
-   DeleteIntRuleArray(PyramidIntRules);
 }
 
 
@@ -1047,8 +1041,6 @@ 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;
@@ -1656,30 +1648,6 @@ 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,7 +323,6 @@ private:
    Array<IntegrationRule *> TriangleIntRules;
    Array<IntegrationRule *> SquareIntRules;
    Array<IntegrationRule *> TetrahedronIntRules;
-   Array<IntegrationRule *> PyramidIntRules;
    Array<IntegrationRule *> PrismIntRules;
    Array<IntegrationRule *> CubeIntRules;
 
@@ -352,7 +351,6 @@ 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,16 +78,6 @@ 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

@@ -1,240 +0,0 @@
-// 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

@@ -1,94 +0,0 @@
-#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

@@ -1,23 +0,0 @@
-// 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

@@ -1,53 +0,0 @@
-// 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,16 +48,6 @@
 #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

@@ -1,157 +0,0 @@
-#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

@@ -1,77 +0,0 @@
-#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 ARPACK
+ OCCA CEED RAJA UMPIRE MKL_CPARDISO AMGX CALIPER
 
 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) SPECTRA LIBUNWIND OPENMP CUDA HIP
+MFEM_DEPENDENCIES = $(MFEM_REQ_LIB_DEPS) 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_USE_ARPACK MFEM_USE_SPECTRA MFEM_SOURCE_DIR MFEM_INSTALL_DIR
+ MFEM_USE_CALIPER 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,8 +652,6 @@ 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,7 +19,6 @@ set(SRCS
   ncmesh.cpp
   nurbs.cpp
   point.cpp
-  pyramid.cpp
   quadrilateral.cpp
   segment.cpp
   tetrahedron.cpp
@@ -39,7 +38,6 @@ 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, PYRAMID
+               TETRAHEDRON, HEXAHEDRON, WEDGE
              };
 
    /// Default element constructor.

mesh/mesh.cpp

@@ -335,7 +335,6 @@ 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;
@@ -736,31 +735,6 @@ 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)
 {
@@ -807,29 +781,6 @@ 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)
@@ -911,19 +862,10 @@ 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_ABORT("Mesh::GetLocalFaceTransformation not defined for "
-                       "face type " << face_type
-                       << " and element type " << elem_type << "\n");
+            MFEM_ASSERT(elem_type == Element::WEDGE, "");
+            GetLocalTriToWdgTransformation(Transf, info);
          }
          break;
 
@@ -932,19 +874,10 @@ 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_ABORT("Mesh::GetLocalFaceTransformation not defined for "
-                       "face type " << face_type
-                       << " and element type " << elem_type << "\n");
+            MFEM_ASSERT(elem_type == Element::WEDGE, "");
+            GetLocalQuadToWdgTransformation(Transf, info);
          }
          break;
    }
@@ -1437,20 +1370,6 @@ 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)
 {
@@ -1504,25 +1423,6 @@ 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);
@@ -2792,16 +2692,11 @@ 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[9];
+   int ind[8];
 
    // Sets vertices and the corresponding coordinates
    for (z = 0; z <= nz; z++)
@@ -2817,25 +2712,8 @@ 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)
@@ -2886,11 +2764,6 @@ 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);
@@ -3555,7 +3428,6 @@ 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);
    }
@@ -3623,8 +3495,6 @@ void Mesh::SetMeshGen()
             mesh_geoms |= (1 << Geometry::TETRAHEDRON);
          case Element::TRIANGLE:
             mesh_geoms |= (1 << Geometry::TRIANGLE);
-         case Element::SEGMENT:
-            mesh_geoms |= (1 << Geometry::SEGMENT);
          case Element::POINT:
             mesh_geoms |= (1 << Geometry::POINT);
             meshgen |= 1;
@@ -3636,6 +3506,8 @@ void Mesh::SetMeshGen()
             mesh_geoms |= (1 << Geometry::SQUARE);
             mesh_geoms |= (1 << Geometry::SEGMENT);
             mesh_geoms |= (1 << Geometry::POINT);
+         case Element::SEGMENT:
+            mesh_geoms |= (1 << Geometry::SEGMENT);
             meshgen |= 2;
             break;
 
@@ -3648,15 +3520,6 @@ 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;
@@ -5197,19 +5060,6 @@ 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);
@@ -6292,22 +6142,6 @@ 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++)
@@ -6400,20 +6234,6 @@ 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++)
@@ -6488,22 +6308,6 @@ 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
@@ -7840,18 +7644,6 @@ 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
@@ -7911,7 +7703,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 + 2 * pyr_counter);
+   new_elements.SetSize(8 * NumOfElements);
    CoarseFineTr.embeddings.SetSize(new_elements.Size());
 
    hex_counter = 0;
@@ -8177,73 +7969,6 @@ 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);
@@ -8375,7 +8100,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, D = -1.0;
+   static const double A = 0.0, B = 0.5, C = 1.0;
    static double tet_children[3*4*16] =
    {
       A,A,A, B,A,A, A,B,A, A,A,B,
@@ -8401,19 +8126,6 @@ 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,
@@ -8439,8 +8151,6 @@ 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]
@@ -8457,7 +8167,7 @@ void Mesh::UniformRefinement3D_base(Array<int> *f2qf_ptr, DSTable *v_to_v_p,
    }
 
    NumOfVertices    = vertices.Size();
-   NumOfElements    = 8 * NumOfElements + 2 * pyr_counter;
+   NumOfElements    = 8 * NumOfElements;
    NumOfBdrElements = 4 * NumOfBdrElements;
 
    GetElementToFaceTable();
@@ -9041,7 +8751,7 @@ void Mesh::GeneralRefinement(const Array<Refinement> &refinements,
    else if (nonconforming < 0)
    {
       // determine if nonconforming refinement is suitable
-      if ((meshgen & 2) || (meshgen & 4) || (meshgen & 8))
+      if ((meshgen & 2) || (meshgen & 4))
       {
          nonconforming = 1; // tensor product elements and wedges
       }
@@ -9825,7 +9535,6 @@ 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,7 +177,6 @@ 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];
 
@@ -196,7 +195,6 @@ 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 };
 
@@ -375,17 +373,11 @@ 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. */
@@ -665,15 +657,11 @@ 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);
@@ -841,16 +829,10 @@ 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,7 +27,6 @@
 #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,11 +352,6 @@ 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] =
@@ -550,8 +545,6 @@ 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;
@@ -1400,10 +1393,6 @@ 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;
@@ -1458,7 +1447,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::PYRAMID)
+            type == Element::HEXAHEDRON)
    {
       MFEM_VERIFY(nx > 0 && ny > 0 && nz > 0 &&
                   sx > 0.0 && sy > 0.0 && sz > 0.0,
@@ -1895,9 +1884,6 @@ 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
@@ -1948,19 +1934,17 @@ void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
                      vert_indices[vi] = it->second;
                   }
 
-                  // 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.
+                  // non-positive attributes are not allowed in MFEM
                   if (phys_domain <= 0)
                   {
-                     has_nonpositive_phys_domain = true;
-                     phys_domain = 1;
-                  }
-                  else
-                  {
-                     has_positive_phys_domain = true;
+                     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.");
                   }
 
                   // initialize the mesh element
@@ -2177,19 +2161,17 @@ void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
                   vert_indices[vi] = it->second;
                }
 
-               // 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.
+               // non-positive attributes are not allowed in MFEM
                if (phys_domain <= 0)
                {
-                  has_nonpositive_phys_domain = true;
-                  phys_domain = 1;
-               }
-               else
-               {
-                  has_positive_phys_domain = true;
+                  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.");
                }
 
                // initialize the mesh element
@@ -2374,24 +2356,6 @@ 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/ncmesh.cpp

@@ -4253,6 +4253,19 @@ void NCMesh::GetPointMatrix(Geometry::Type geom, const char* ref_path,
             pm = PointMatrix(mid12, mid20, mid01);
          }
       }
+      else if (geom == Geometry::SEGMENT)
+      {
+         Point mid01(pm(0), pm(1));
+
+         if (child == 0)
+         {
+            pm = PointMatrix(pm(0), mid01);
+         }
+         else if (child == 1)
+         {
+            pm = PointMatrix(mid01, pm(1));
+         }
+      }
    }
 
    // write the points to the matrix

mesh/pyramid.cpp

@@ -1,64 +0,0 @@
-// 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

@@ -1,78 +0,0 @@
-// 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,7 +71,8 @@ public:
    virtual ~Wedge() { }
 };
 
-extern class LinearWedgeFiniteElement WedgeFE;
+// Defined in fe.cpp to ensure construction after 'mfem::poly1d'.
+extern class H1_WedgeElement WedgeFE;
 
 }