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

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

.gitignore

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

CHANGELOG

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

config/config.hpp.in

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

config/config.mk.in

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

config/defaults.mk

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

data/sphere.geo

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

examples/arpack/ex11.cpp

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

examples/arpack/makefile

@@ -0,0 +1,69 @@
+# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+# Use the MFEM build directory
+MFEM_DIR ?= ../..
+MFEM_BUILD_DIR ?= ../..
+SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/examples/arpack/,)
+CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
+# Use the MFEM install directory
+# MFEM_INSTALL_DIR = ../../mfem
+# CONFIG_MK = $(MFEM_INSTALL_DIR)/share/mfem/config.mk
+
+MFEM_LIB_FILE = mfem_is_not_built
+-include $(CONFIG_MK)
+
+SEQ_EXAMPLES = ex11
+PAR_EXAMPLES = 
+ifeq ($(MFEM_USE_MPI),NO)
+   EXAMPLES = $(SEQ_EXAMPLES)
+else
+   EXAMPLES = $(PAR_EXAMPLES)
+endif
+RC_FILES = $(patsubst $(SRC)%,%,$(wildcard $(SRC)rc_*))
+
+.SUFFIXES:
+.SUFFIXES: .o .cpp .mk
+.PHONY: all clean clean-build clean-exec
+
+# Remove built-in rule
+%: %.cpp
+
+# Replace the default implicit rule for *.cpp files
+%: $(SRC)%.cpp $(MFEM_LIB_FILE) $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ $(MFEM_LIBS)
+
+all: $(EXAMPLES)
+
+# Examples depend on their corresponding rc_* files:
+make-rc-rule = $(1): | $(filter rc_$(1)%,$(RC_FILES))
+$(foreach ex,$(EXAMPLES),$(eval $(call make-rc-rule,$(ex))))
+
+# Rules to copy the rc_* files when building out-of-source:
+ifneq ($(SRC),)
+$(RC_FILES): %: $(SRC)%
+	cp -pf $(<) .
+endif
+
+# Generate an error message if the MFEM library is not built and exit
+$(MFEM_LIB_FILE):
+	$(error The MFEM library is not built)
+
+clean: clean-build clean-exec
+
+clean-build:
+	rm -f *.o *~ $(SEQ_EXAMPLES) $(PAR_EXAMPLES)
+	rm -rf *.dSYM *.TVD.*breakpoints
+
+clean-exec:
+	@rm -rf mesh.* sol.* sol_p.* sol_u.* Example5*
+	@rm -f ex9-mesh.* ex9-init.* ex9-final.* Example9*
+	@rm -f deformed.* velocity.* elastic_energy.*

examples/spectra/ex11.cpp

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

examples/spectra/makefile

@@ -0,0 +1,67 @@
+# Copyright (c) 2010-2021, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+# Use the MFEM build directory
+MFEM_DIR ?= ../..
+MFEM_BUILD_DIR ?= ../..
+SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/examples/spectra/,)
+CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
+# Use the MFEM install directory
+# MFEM_INSTALL_DIR = ../../mfem
+# CONFIG_MK = $(MFEM_INSTALL_DIR)/share/mfem/config.mk
+
+MFEM_LIB_FILE = mfem_is_not_built
+-include $(CONFIG_MK)
+
+SEQ_EXAMPLES = ex11
+PAR_EXAMPLES = 
+ifeq ($(MFEM_USE_MPI),NO)
+   EXAMPLES = $(SEQ_EXAMPLES)
+else
+   EXAMPLES = $(PAR_EXAMPLES)
+endif
+RC_FILES = $(patsubst $(SRC)%,%,$(wildcard $(SRC)rc_*))
+
+.SUFFIXES:
+.SUFFIXES: .o .cpp .mk
+.PHONY: all clean clean-build clean-exec
+
+# Remove built-in rule
+%: %.cpp
+
+# Replace the default implicit rule for *.cpp files
+%: $(SRC)%.cpp $(MFEM_LIB_FILE) $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ $(MFEM_LIBS)
+
+all: $(EXAMPLES)
+
+# Examples depend on their corresponding rc_* files:
+make-rc-rule = $(1): | $(filter rc_$(1)%,$(RC_FILES))
+$(foreach ex,$(EXAMPLES),$(eval $(call make-rc-rule,$(ex))))
+
+# Rules to copy the rc_* files when building out-of-source:
+ifneq ($(SRC),)
+$(RC_FILES): %: $(SRC)%
+	cp -pf $(<) .
+endif
+
+# Generate an error message if the MFEM library is not built and exit
+$(MFEM_LIB_FILE):
+	$(error The MFEM library is not built)
+
+clean: clean-build clean-exec
+
+clean-build:
+	rm -f *.o *~ $(SEQ_EXAMPLES) $(PAR_EXAMPLES)
+	rm -rf *.dSYM *.TVD.*breakpoints
+
+clean-exec:
+	@rm -rf *.mesh mode_*

linalg/CMakeLists.txt

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

linalg/arpack.hpp

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

linalg/eigen.hpp

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

linalg/eigensolver.cpp

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

linalg/eigensolver.hpp

@@ -0,0 +1,53 @@
+// Copyright (c) 2010-2020, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_EIGENSOLVERS
+#define MFEM_EIGENSOLVERS
+
+#include "vector.hpp"
+#include "operator.hpp"
+
+namespace mfem
+{
+
+/// Abstract Eigensolver
+class Eigensolver
+{
+public:
+
+   Eigensolver();
+   virtual ~Eigensolver() {}
+
+   virtual void SetTol(double tol) = 0;
+   virtual void SetMaxIter(int max_iter) = 0;
+   virtual void SetPrintLevel(int logging) = 0;
+   virtual void SetNumModes(int num_eigs) = 0;
+
+   virtual void SetOperator(Operator & A) = 0;
+   virtual void SetMassMatrix(Operator & M) = 0;
+
+   /// Perform the eigenvalue solve
+   virtual void Solve() = 0;
+
+   /// Collect the converged eigenvalues
+   virtual void GetEigenvalues(Array<double> & eigenvalues) = 0;
+
+   /// Extract a single eigenvector
+   virtual Vector & GetEigenvector(unsigned int i) = 0;
+
+   /// Transfer ownership of the converged eigenvectors
+   virtual Vector ** StealEigenvectors() = 0;
+
+};
+
+}
+
+#endif

linalg/linalg.hpp

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

linalg/spectra.cpp

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

linalg/spectra.hpp

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

makefile

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

mesh/mesh_readers.cpp

@@ -1895,6 +1895,9 @@ void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
          ho_wdg[2] = wdg18; ho_wdg[3] = wdg40;
          ho_pyr[2] = pyr14; ho_pyr[3] = pyr30;
 
+         bool has_nonpositive_phys_domain = false;
+         bool has_positive_phys_domain = false;
+
          if (binary)
          {
             int n_elem_part = 0; // partial sum of elements that are read
@@ -1945,17 +1948,19 @@ void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
                      vert_indices[vi] = it->second;
                   }
 
-                  // non-positive attributes are not allowed in MFEM
+                  // Non-positive attributes are not allowed in MFEM. However,
+                  // by default, Gmsh sets the physical domain of all elements
+                  // to zero. In the case that all elements have physical domain
+                  // zero, we will given them attribute 1. If only some elements
+                  // have physical domain zero, we will throw an error.
                   if (phys_domain <= 0)
                   {
-                     MFEM_ABORT("Non-positive element attribute in Gmsh mesh!\n"
-                                "By default Gmsh sets element tags (attributes)"
-                                " to '0' but MFEM requires that they be"
-                                " positive integers.\n"
-                                "Use \"Physical Curve\", \"Physical Surface\","
-                                " or \"Physical Volume\" to set tags/attributes"
-                                " for all curves, surfaces, or volumes in your"
-                                " Gmsh geometry to values which are >= 1.");
+                     has_nonpositive_phys_domain = true;
+                     phys_domain = 1;
+                  }
+                  else
+                  {
+                     has_positive_phys_domain = true;
                   }
 
                   // initialize the mesh element
@@ -2172,17 +2177,19 @@ void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
                   vert_indices[vi] = it->second;
                }
 
-               // non-positive attributes are not allowed in MFEM
+               // Non-positive attributes are not allowed in MFEM. However,
+               // by default, Gmsh sets the physical domain of all elements
+               // to zero. In the case that all elements have physical domain
+               // zero, we will given them attribute 1. If only some elements
+               // have physical domain zero, we will throw an error.
                if (phys_domain <= 0)
                {
-                  MFEM_ABORT("Non-positive element attribute in Gmsh mesh!\n"
-                             "By default Gmsh sets element tags (attributes)"
-                             " to '0' but MFEM requires that they be"
-                             " positive integers.\n"
-                             "Use \"Physical Curve\", \"Physical Surface\","
-                             " or \"Physical Volume\" to set tags/attributes"
-                             " for all curves, surfaces, or volumes in your"
-                             " Gmsh geometry to values which are >= 1.");
+                  has_nonpositive_phys_domain = true;
+                  phys_domain = 1;
+               }
+               else
+               {
+                  has_positive_phys_domain = true;
                }
 
                // initialize the mesh element
@@ -2367,6 +2374,24 @@ void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
             } // el (all elements)
          } // if ASCII
 
+         if (has_positive_phys_domain && has_nonpositive_phys_domain)
+         {
+            MFEM_ABORT("Non-positive element attribute in Gmsh mesh!\n"
+                       "By default Gmsh sets element tags (attributes)"
+                       " to '0' but MFEM requires that they be"
+                       " positive integers.\n"
+                       "Use \"Physical Curve\", \"Physical Surface\","
+                       " or \"Physical Volume\" to set tags/attributes"
+                       " for all curves, surfaces, or volumes in your"
+                       " Gmsh geometry to values which are >= 1.");
+         }
+         else if (has_nonpositive_phys_domain)
+         {
+            mfem::out << "\nGmsh reader: all element attributes were zero.\n"
+                      << "MFEM only supports positive element attributes.\n"
+                      << "Setting element attributes to 1.\n\n";
+         }
+
          if (!elements_3D.empty())
          {
             Dim = 3;