make style

2. Add _USE_MATH_DEFINES to CmakeLists.txt (without this change VS doesn't "see" this preprocessor in config.hpp)
3. Remove asm definition - asm is not allowed in VS anymore
4. Convert "const double Epsilon" to constexpr (fixes CUDA compilation error)

.gitignore

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

CHANGELOG

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

CMakeLists.txt

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

config/config.hpp.in

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

config/config.mk.in

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

config/defaults.mk

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

data/sphere.geo

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

examples/arpack/ex11.cpp

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

examples/arpack/makefile

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

examples/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

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

examples/spectra/makefile

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

fem/bilininteg_convection_ea.cpp

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

fem/bilininteg_convection_pa.cpp

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

fem/bilininteg_dgtrace_ea.cpp

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

fem/bilininteg_dgtrace_pa.cpp

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

fem/bilininteg_diffusion_ea.cpp

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

fem/bilininteg_diffusion_pa.cpp

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

fem/bilininteg_divergence.cpp

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

fem/bilininteg_gradient.cpp

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

fem/bilininteg_hcurl.cpp

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

fem/bilininteg_hdiv.cpp

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

fem/bilininteg_mass_ea.cpp

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

fem/bilininteg_mass_pa.cpp

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

fem/bilininteg_vecdiffusion.cpp

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

fem/bilininteg_vecmass.cpp

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

fem/nonlininteg_vectorconvection.cpp

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

fem/qinterp/det.cpp

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

fem/qinterp/eval.hpp

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

fem/qinterp/grad.hpp

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

fem/quadinterpolator.cpp

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

fem/tmop.cpp

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

fem/tmop_amr.cpp

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

general/error.hpp

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

general/isockstream.cpp

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

linalg/CMakeLists.txt

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

linalg/arpack.hpp

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

linalg/eigen.hpp

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

linalg/eigensolver.cpp

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

linalg/eigensolver.hpp

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

linalg/kernels.hpp

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

linalg/linalg.hpp

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

linalg/spectra.cpp

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

linalg/spectra.hpp

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

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 ARPACK
+ OCCA CEED RAJA UMPIRE MKL_CPARDISO AMGX CALIPER
 
 PETSC_ERROR_MSG = $(if $(PETSC_FOUND),,. PETSC config not found: $(PETSC_VARS))
 SLEPC_ERROR_MSG = $(if $(SLEPC_FOUND),,. SLEPC config not found: $(SLEPC_VARS))
@@ -292,7 +292,7 @@ ifeq ($(MAKECMDGOALS),config)
 endif
 
 # List of MFEM dependencies, processed below
-MFEM_DEPENDENCIES = $(MFEM_REQ_LIB_DEPS) SPECTRA LIBUNWIND OPENMP CUDA HIP
+MFEM_DEPENDENCIES = $(MFEM_REQ_LIB_DEPS) LIBUNWIND OPENMP CUDA HIP
 
 # List of deprecated MFEM dependencies, processed below
 MFEM_LEGACY_DEPENDENCIES = OPENMP
@@ -340,7 +340,7 @@ MFEM_DEFINES = MFEM_VERSION MFEM_VERSION_STRING MFEM_GIT_STRING MFEM_USE_MPI\
  MFEM_USE_PUMI MFEM_USE_HIOP MFEM_USE_GSLIB MFEM_USE_CUDA MFEM_USE_HIP\
  MFEM_USE_OCCA MFEM_USE_CEED MFEM_USE_RAJA MFEM_USE_UMPIRE MFEM_USE_SIMD\
  MFEM_USE_ADIOS2 MFEM_USE_MKL_CPARDISO MFEM_USE_AMGX MFEM_USE_MUMPS\
- MFEM_USE_CALIPER MFEM_USE_ARPACK MFEM_USE_SPECTRA MFEM_SOURCE_DIR MFEM_INSTALL_DIR
+ MFEM_USE_CALIPER MFEM_SOURCE_DIR MFEM_INSTALL_DIR
 
 # List of makefile variables that will be written to config.mk:
 MFEM_CONFIG_VARS = MFEM_CXX MFEM_HOST_CXX MFEM_CPPFLAGS MFEM_CXXFLAGS\
@@ -652,8 +652,6 @@ status info:
 	$(info MFEM_USE_SUNDIALS      = $(MFEM_USE_SUNDIALS))
 	$(info MFEM_USE_MESQUITE      = $(MFEM_USE_MESQUITE))
 	$(info MFEM_USE_SUITESPARSE   = $(MFEM_USE_SUITESPARSE))
-	$(info MFEM_USE_ARPACK        = $(MFEM_USE_ARPACK))
-	$(info MFEM_USE_SPECTRA       = $(MFEM_USE_SPECTRA))
 	$(info MFEM_USE_SUPERLU       = $(MFEM_USE_SUPERLU))
 	$(info MFEM_USE_MUMPS         = $(MFEM_USE_MUMPS))
 	$(info MFEM_USE_STRUMPACK     = $(MFEM_USE_STRUMPACK))

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

mesh/ncmesh.cpp

@@ -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};
 
-static double qf(const int order, const int ker, Mesh &m,
-                 FiniteElementSpace &fes, GridFunction &u)
+double qf(const int order, const int ker, Mesh &m,
+          FiniteElementSpace &fes, GridFunction &u)
 {
    const Geometry::Type type = m.GetElementBaseGeometry(0);
    const IntegrationRule &ir(IntRules.Get(type, order));

miniapps/shifted/dist_solver.cpp

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

tests/convergence/prates.cpp

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

tests/convergence/rates.cpp

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

tests/unit/fem/test_doftrans.cpp

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