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
 ======================================

CMakeLists.txt

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

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

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

examples/ex25p.cpp

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

examples/spectra/ex11.cpp

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

examples/spectra/makefile

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

fem/bilininteg_convection_ea.cpp

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

fem/bilininteg_convection_pa.cpp

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

fem/bilininteg_dgtrace_ea.cpp

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

fem/bilininteg_dgtrace_pa.cpp

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

fem/bilininteg_diffusion_ea.cpp

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

fem/bilininteg_diffusion_pa.cpp

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

fem/bilininteg_divergence.cpp

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

fem/bilininteg_gradient.cpp

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

fem/bilininteg_hcurl.cpp

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

fem/bilininteg_hdiv.cpp

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

fem/bilininteg_mass_ea.cpp

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

fem/bilininteg_mass_pa.cpp

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

fem/bilininteg_vecdiffusion.cpp

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

fem/bilininteg_vecmass.cpp

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

fem/nonlininteg_vectorconvection.cpp

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

fem/qinterp/det.cpp

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

fem/qinterp/eval.hpp

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

fem/qinterp/grad.hpp

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

fem/quadinterpolator.cpp

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

fem/tmop.cpp

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

fem/tmop_amr.cpp

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

general/error.hpp

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

general/isockstream.cpp

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

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/kernels.hpp

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

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

linalg/sundials.cpp

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

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

@@ -2850,14 +2850,14 @@ void Mesh::Make3D(int nx, int ny, int nz, Element::Type type,
          y = sfc[3*k + 1];
          z = sfc[3*k + 2];
 
-         ind[0] = VTX(x, y, z  );
-         ind[1] = VTX(x+1, y, z  );
+         ind[0] = VTX(x  , y  , z  );
+         ind[1] = VTX(x+1, y  , z  );
          ind[2] = VTX(x+1, y+1, z  );
-         ind[3] = VTX(x, y+1, z  );
-         ind[4] = VTX(x, y, z+1);
-         ind[5] = VTX(x+1, y, z+1);
+         ind[3] = VTX(x  , y+1, z  );
+         ind[4] = VTX(x  , y  , z+1);
+         ind[5] = VTX(x+1, y  , z+1);
          ind[6] = VTX(x+1, y+1, z+1);
-         ind[7] = VTX(x, y+1, z+1);
+         ind[7] = VTX(x  , y+1, z+1);
 
          AddHex(ind, 1);
       }
@@ -2870,12 +2870,12 @@ void Mesh::Make3D(int nx, int ny, int nz, Element::Type type,
          {
             for (x = 0; x < nx; x++)
             {
-               ind[0] = VTX(x, y, z  );
-               ind[1] = VTX(x+1, y, z  );
+               ind[0] = VTX(x  , y  , z  );
+               ind[1] = VTX(x+1, y  , z  );
                ind[2] = VTX(x+1, y+1, z  );
-               ind[3] = VTX(x, y+1, z  );
-               ind[4] = VTX(x, y, z+1);
-               ind[5] = VTX(x+1, y, z+1);
+               ind[3] = VTX(x  , y+1, z  );
+               ind[4] = VTX(x  , y  , z+1);
+               ind[5] = VTX(x+1, y  , z+1);
                ind[6] = VTX(x+1, y+1, z+1);
                ind[7] = VTX(  x, y+1, z+1);
                if (type == Element::TETRAHEDRON)
@@ -2906,10 +2906,10 @@ void Mesh::Make3D(int nx, int ny, int nz, Element::Type type,
    {
       for (x = 0; x < nx; x++)
       {
-         ind[0] = VTX(x, y, 0);
-         ind[1] = VTX(x, y+1, 0);
+         ind[0] = VTX(x  , y  , 0);
+         ind[1] = VTX(x  , y+1, 0);
          ind[2] = VTX(x+1, y+1, 0);
-         ind[3] = VTX(x+1, y, 0);
+         ind[3] = VTX(x+1, y  , 0);
          if (type == Element::TETRAHEDRON)
          {
             AddBdrQuadAsTriangles(ind, 1);
@@ -2929,10 +2929,10 @@ void Mesh::Make3D(int nx, int ny, int nz, Element::Type type,
    {
       for (x = 0; x < nx; x++)
       {
-         ind[0] = VTX(x, y, nz);
-         ind[1] = VTX(x+1, y, nz);
+         ind[0] = VTX(x  , y  , nz);
+         ind[1] = VTX(x+1, y  , nz);
          ind[2] = VTX(x+1, y+1, nz);
-         ind[3] = VTX(x, y+1, nz);
+         ind[3] = VTX(x  , y+1, nz);
          if (type == Element::TETRAHEDRON)
          {
             AddBdrQuadAsTriangles(ind, 6);
@@ -2952,10 +2952,10 @@ void Mesh::Make3D(int nx, int ny, int nz, Element::Type type,
    {
       for (y = 0; y < ny; y++)
       {
-         ind[0] = VTX(0, y, z  );
-         ind[1] = VTX(0, y, z+1);
-         ind[2] = VTX(0, y+1, z+1);
-         ind[3] = VTX(0, y+1, z  );
+         ind[0] = VTX(0  , y  , z  );
+         ind[1] = VTX(0  , y  , z+1);
+         ind[2] = VTX(0  , y+1, z+1);
+         ind[3] = VTX(0  , y+1, z  );
          if (type == Element::TETRAHEDRON)
          {
             AddBdrQuadAsTriangles(ind, 5);
@@ -2971,10 +2971,10 @@ void Mesh::Make3D(int nx, int ny, int nz, Element::Type type,
    {
       for (y = 0; y < ny; y++)
       {
-         ind[0] = VTX(nx, y, z  );
+         ind[0] = VTX(nx, y  , z  );
          ind[1] = VTX(nx, y+1, z  );
          ind[2] = VTX(nx, y+1, z+1);
-         ind[3] = VTX(nx, y, z+1);
+         ind[3] = VTX(nx, y  , z+1);
          if (type == Element::TETRAHEDRON)
          {
             AddBdrQuadAsTriangles(ind, 3);
@@ -2990,10 +2990,10 @@ void Mesh::Make3D(int nx, int ny, int nz, Element::Type type,
    {
       for (z = 0; z < nz; z++)
       {
-         ind[0] = VTX(x, 0, z  );
+         ind[0] = VTX(x  , 0, z  );
          ind[1] = VTX(x+1, 0, z  );
          ind[2] = VTX(x+1, 0, z+1);
-         ind[3] = VTX(x, 0, z+1);
+         ind[3] = VTX(x  , 0, z+1);
          if (type == Element::TETRAHEDRON)
          {
             AddBdrQuadAsTriangles(ind, 2);
@@ -3009,8 +3009,8 @@ void Mesh::Make3D(int nx, int ny, int nz, Element::Type type,
    {
       for (z = 0; z < nz; z++)
       {
-         ind[0] = VTX(x, ny, z  );
-         ind[1] = VTX(x, ny, z+1);
+         ind[0] = VTX(x  , ny, z  );
+         ind[1] = VTX(x  , ny, z+1);
          ind[2] = VTX(x+1, ny, z+1);
          ind[3] = VTX(x+1, ny, z  );
          if (type == Element::TETRAHEDRON)

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;

mesh/ncmesh.cpp

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

miniapps/meshing/minimal-surface.cpp

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

miniapps/shifted/dist_solver.cpp

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

tests/convergence/prates.cpp

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

tests/convergence/rates.cpp

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

tests/unit/fem/test_doftrans.cpp

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