Fix signed char issue in Device::GetUUID

older hypre versions.

.github/workflows/builds-and-tests.yml

@@ -295,8 +295,7 @@ jobs:
           export HOMEBREW_NO_INSTALL_CLEANUP=1
           brew update
           brew install enzyme
-          ENZYME_LLVM=$(brew info enzyme | sed -n 's/^Required.*:.*\(llvm[^ ]*\).*/\1/p')
-          echo "ENZYME_LLVM=$ENZYME_LLVM"
+          ENZYME_LLVM=$(brew info enzyme | sed -n 's/^Required:.*\(llvm[^ ]*\).*/\1/p')
           LLVM_PREFIX=$(brew --prefix $ENZYME_LLVM)
           echo "LLVM_PREFIX=$LLVM_PREFIX" >> $GITHUB_ENV
           echo "OMPI_CC=$LLVM_PREFIX/bin/clang" >> $GITHUB_ENV

.gitignore

@@ -443,10 +443,6 @@ miniapps/diag-smoothers/mg-abs-l1-jacobi
 miniapps/contact/contact
 miniapps/contact/ParaView
 
-miniapps/plasma/pic/electrostatic-*
-!miniapps/plasma/pic/electrostatic-*.cpp
-miniapps/plasma/pic/*.csv
-
 # Unit test binary and outputs
 tests/unit/output_meshes
 tests/unit/unit_tests

config/cmake/MFEMConfig.cmake.in

@@ -109,10 +109,6 @@ if (MFEM_USE_RAJA)
    find_dependency(RAJA)
 endif()
 
-if (MFEM_USE_UMPIRE)
-   find_dependency(umpire)
-endif()
-
 if (NOT TARGET mfem)
   include(${CMAKE_CURRENT_LIST_DIR}/MFEMTargets.cmake)
 endif (NOT TARGET mfem)

config/cmake/modules/FindUMPIRE.cmake

@@ -14,12 +14,12 @@
 #   - UMPIRE_LIBRARIES
 #   - UMPIRE_INCLUDE_DIRS
 
-if (NOT umpire_ROOT AND UMPIRE_DIR)
-  set(umpire_ROOT ${UMPIRE_DIR})
+if (NOT umpire_DIR AND UMPIRE_DIR)
+  set(umpire_DIR ${UMPIRE_DIR}/lib/cmake/umpire)
 endif()
 message(STATUS "Looking for UMPIRE ...")
 message(STATUS "   in UMPIRE_DIR = ${UMPIRE_DIR}")
-message(STATUS "      umpire_ROOT = ${umpire_ROOT}")
+message(STATUS "      umpire_DIR = ${umpire_DIR}")
 find_package(umpire CONFIG)
 set(UMPIRE_FOUND ${umpire_FOUND})
 set(UMPIRE_LIBRARIES "umpire")

config/config.hpp.in

@@ -97,9 +97,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 SuperLU_DIST library.
 // #define MFEM_USE_SUPERLU
 // #define MFEM_USE_SUPERLU5

config/config.mk.in

@@ -32,7 +32,6 @@ MFEM_USE_MEMALLOC      = @MFEM_USE_MEMALLOC@
 MFEM_TIMER_TYPE        = @MFEM_TIMER_TYPE@
 MFEM_USE_SUNDIALS      = @MFEM_USE_SUNDIALS@
 MFEM_USE_SUITESPARSE   = @MFEM_USE_SUITESPARSE@
-MFEM_USE_ARPACK        = @MFEM_USE_ARPACK@
 MFEM_USE_SUPERLU       = @MFEM_USE_SUPERLU@
 MFEM_USE_SUPERLU5      = @MFEM_USE_SUPERLU5@
 MFEM_USE_MUMPS         = @MFEM_USE_MUMPS@

config/defaults.mk

@@ -178,7 +178,6 @@ MFEM_USE_ALGOIM        = NO
 MFEM_USE_UMPIRE        = NO
 MFEM_USE_SIMD          = NO
 MFEM_USE_ADIOS2        = NO
-MFEM_USE_ARPACK        = NO
 MFEM_USE_MKL_CPARDISO  = NO
 MFEM_USE_MKL_PARDISO   = NO
 MFEM_USE_MOONOLITH     = NO
@@ -428,14 +427,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
-ifeq ($(MFEM_USE_MPI),YES)
-ARPACK_LIB = -L$(ARPACK_DIR) -lparpack -larpack
-else
-ARPACK_LIB = -L$(ARPACK_DIR) -larpack
-endif
-
 # 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)

examples/CMakeLists.txt

@@ -49,13 +49,6 @@ list(APPEND ALL_EXE_SRCS
   ex41.cpp
   )
 
-if (MFEM_USE_ARPACK)
-list(APPEND ALL_EXE_SRCS
-  ex11.pp
-  ex13.pp
-)
-endif()
-
 if (MFEM_USE_MPI)
   list(APPEND ALL_EXE_SRCS
     ex0p.cpp

examples/ex11.cpp

@@ -1,298 +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;
-
-#ifdef MFEM_USE_ARPACK
-
-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;
-   bool arp_solver = true;
-
-   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();
-
-   Solver       * solver = NULL;
-
-#ifndef MFEM_USE_SUITESPARSE
-   // 6. 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());
-
-   // 7. Define and configure the ARPACK eigensolver
-   SymGenEigensolver * eig_solver = NULL;
-
-   if (arp_solver)
-   {
-      // ArPackSymGen * arpack = new ArPackSymGen();
-      ArPackSAUPD * arpack = new ArPackSAUPD();
-      arpack->SetMode(2);
-      arpack->SetNumModes(nev);
-      arpack->SetMaxIter(400);
-      arpack->SetTol(1e-8);
-      arpack->SetPrintLevel(2);
-      arpack->SetSolver(*solver);
-
-      eig_solver = arpack;
-   }
-
-   eig_solver->SetOperators(*a, *m);
-
-   // 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;
-   eig_solver->Solve();
-   eig_solver->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 Vector to GridFunction
-         x = eig_solver->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 Vector to GridFunction
-         x = eig_solver->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 eig_solver;
-   delete solver;
-   delete m;
-   delete a;
-
-   delete fespace;
-   if (order > 0)
-   {
-      delete fec;
-   }
-   delete mesh;
-
-   return 0;
-}
-
-#endif // MFEM_USE_ARPACK

examples/ex11p.cpp

@@ -72,8 +72,6 @@ int main(int argc, char *argv[])
    int seed = 75;
    bool slu_solver  = false;
    bool sp_solver = false;
-   bool lob_solver = true;
-   bool arp_solver = false;
    bool cpardiso_solver = false;
    bool visualization = 1;
 
@@ -99,10 +97,6 @@ int main(int argc, char *argv[])
    args.AddOption(&sp_solver, "-sp", "--strumpack", "-no-sp",
                   "--no-strumpack", "Use the STRUMPACK Solver.");
 #endif
-#ifdef MFEM_USE_ARPACK
-   args.AddOption(&arp_solver, "-arp", "--arpack", "-no-arp",
-                  "--no-arpack", "Use the Parallel ARPACK Solver.");
-#endif
 #ifdef MFEM_USE_MKL_CPARDISO
    args.AddOption(&cpardiso_solver, "-cpardiso", "--cpardiso", "-no-cpardiso",
                   "--no-cpardiso", "Use the MKL CPardiso Solver.");
@@ -119,11 +113,6 @@ int main(int argc, char *argv[])
               << "         Defaulting to SuperLU." << endl;
       sp_solver = false;
    }
-   if (arp_solver)
-   {
-      lob_solver = false;
-   }
-
    // The command line options are also passed to the STRUMPACK
    // solver. So do not exit if some options are not recognized.
    if (!sp_solver)
@@ -254,119 +243,70 @@ int main(int argc, char *argv[])
    // 8. Define and configure the LOBPCG eigensolver and the BoomerAMG
    //    preconditioner for A to be used within the solver. Set the matrices
    //    which define the generalized eigenproblem A x = lambda M x.
-   Solver * solver = NULL;
    Solver * precond = NULL;
    if (!slu_solver && !sp_solver && !cpardiso_solver)
    {
       HypreBoomerAMG * amg = new HypreBoomerAMG(*A);
       amg->SetPrintLevel(0);
       precond = amg;
-
-      if (arp_solver)
-      {
-         HyprePCG * pcg = new HyprePCG(*A);
-         pcg->SetTol(1e-12);
-         pcg->SetPreconditioner(*amg);
-         solver = pcg;
-      }
    }
-#ifdef MFEM_USE_SUPERLU
-   else if (slu_solver)
+   else
    {
-      SuperLUSolver * superlu = new SuperLUSolver(MPI_COMM_WORLD);
-      superlu->SetPrintStatistics(false);
-      superlu->SetSymmetricPattern(true);
-      superlu->SetColumnPermutation(superlu::PARMETIS);
-      superlu->SetOperator(*Arow);
-
-      if (arp_solver)
-      {
-         solver = superlu;
-      }
-      else
+#ifdef MFEM_USE_SUPERLU
+      if (slu_solver)
       {
+         SuperLUSolver * superlu = new SuperLUSolver(MPI_COMM_WORLD);
+         superlu->SetPrintStatistics(false);
+         superlu->SetSymmetricPattern(true);
+         superlu->SetColumnPermutation(superlu::PARMETIS);
+         superlu->SetOperator(*Arow);
          precond = superlu;
       }
-   }
 #endif
 #ifdef MFEM_USE_STRUMPACK
-   else if (sp_solver)
-   {
-      STRUMPACKSolver * strumpack = new STRUMPACKSolver(argc, argv,
-                                                        MPI_COMM_WORLD);
-      strumpack->SetPrintFactorStatistics(true);
-      strumpack->SetPrintSolveStatistics(false);
-      strumpack->SetKrylovSolver(strumpack::KrylovSolver::DIRECT);
-      strumpack->SetReorderingStrategy(strumpack::ReorderingStrategy::METIS);
-      strumpack->SetMatching(strumpack::MatchingJob::NONE);
-      strumpack->SetCompression(strumpack::CompressionType::NONE);
-      strumpack->SetOperator(*Arow);
-      strumpack->SetFromCommandLine();
-      if (arp_solver)
-      {
-         solver = strumpack;
-      }
-      else
+      if (sp_solver)
       {
+         STRUMPACKSolver * strumpack = new STRUMPACKSolver(MPI_COMM_WORLD, argc, argv);
+         strumpack->SetPrintFactorStatistics(true);
+         strumpack->SetPrintSolveStatistics(false);
+         strumpack->SetKrylovSolver(strumpack::KrylovSolver::DIRECT);
+         strumpack->SetReorderingStrategy(strumpack::ReorderingStrategy::METIS);
+         strumpack->SetMatching(strumpack::MatchingJob::NONE);
+         strumpack->SetCompression(strumpack::CompressionType::NONE);
+         strumpack->SetOperator(*Arow);
+         strumpack->SetFromCommandLine();
          precond = strumpack;
       }
-   }
 #endif
 #ifdef MFEM_USE_MKL_CPARDISO
-   else if (cpardiso_solver)
-   {
-      auto cpardiso = new CPardisoSolver(A->GetComm());
-      cpardiso->SetMatrixType(CPardisoSolver::MatType::REAL_STRUCTURE_SYMMETRIC);
-      cpardiso->SetPrintLevel(1);
-      cpardiso->SetOperator(*A);
-      if (arp_solver)
-      {
-         solver = cpardiso;
-      }
-      else
+      if (cpardiso_solver)
       {
+         auto cpardiso = new CPardisoSolver(A->GetComm());
+         cpardiso->SetMatrixType(CPardisoSolver::MatType::REAL_STRUCTURE_SYMMETRIC);
+         cpardiso->SetPrintLevel(1);
+         cpardiso->SetOperator(*A);
          precond = cpardiso;
       }
-   }
 #endif
-
-   SymGenEigensolver * eig_solver = NULL;
-
-   if (lob_solver)
-   {
-      HypreLOBPCG * lobpcg = new HypreLOBPCG(MPI_COMM_WORLD);
-      lobpcg->SetNumModes(nev);
-      lobpcg->SetRandomSeed(seed);
-      lobpcg->SetPreconditioner(*precond);
-      lobpcg->SetMaxIter(200);
-      lobpcg->SetTol(1e-8);
-      lobpcg->SetPrecondUsageMode(1);
-      lobpcg->SetPrintLevel(1);
-
-      eig_solver = lobpcg;
    }
-#ifdef MFEM_USE_ARPACK
-   else if (arp_solver)
-   {
-      ArPackPSAUPD * arpack = new ArPackPSAUPD(MPI_COMM_WORLD);
-      arpack->SetNumModes(nev);
-      arpack->SetMaxIter(400);
-      arpack->SetTol(1e-8);
-      arpack->SetMode(3);
-      arpack->SetPrintLevel(2);
-      arpack->SetSolver(*solver);
 
-      eig_solver = arpack;
-   }
-#endif
-   eig_solver->SetOperators(*A, *M);
+   HypreLOBPCG * lobpcg = new HypreLOBPCG(MPI_COMM_WORLD);
+   lobpcg->SetNumModes(nev);
+   lobpcg->SetRandomSeed(seed);
+   lobpcg->SetPreconditioner(*precond);
+   lobpcg->SetMaxIter(200);
+   lobpcg->SetTol(1e-8);
+   lobpcg->SetPrecondUsageMode(1);
+   lobpcg->SetPrintLevel(1);
+   lobpcg->SetMassMatrix(*M);
+   lobpcg->SetOperator(*A);
 
    // 9. 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<real_t> eigenvalues;
-   eig_solver->Solve();
-   eig_solver->GetEigenvalues(eigenvalues);
+   lobpcg->Solve();
+   lobpcg->GetEigenvalues(eigenvalues);
    ParGridFunction x(fespace);
 
    // 10. Save the refined mesh and the modes in parallel. This output can be
@@ -381,8 +321,8 @@ int main(int argc, char *argv[])
 
       for (int i=0; i<nev; i++)
       {
-         // convert eigenvector from Vector to ParGridFunction
-         x.Distribute(eig_solver->GetEigenvector(i));
+         // convert eigenvector from HypreParVector to ParGridFunction
+         x = lobpcg->GetEigenvector(i);
 
          mode_name << "mode_" << setfill('0') << setw(2) << i << "."
                    << setfill('0') << setw(6) << myid;
@@ -410,8 +350,8 @@ int main(int argc, char *argv[])
                  << ", Lambda = " << eigenvalues[i] << endl;
          }
 
-         // convert eigenvector from Vector to ParGridFunction
-         x.Distribute(eig_solver->GetEigenvector(i));
+         // convert eigenvector from HypreParVector to ParGridFunction
+         x = lobpcg->GetEigenvector(i);
 
          mode_sock << "parallel " << num_procs << " " << myid << "\n"
                    << "solution\n" << *pmesh << x << flush
@@ -435,8 +375,7 @@ int main(int argc, char *argv[])
    }
 
    // 12. Free the used memory.
-   delete eig_solver;
-   delete solver;
+   delete lobpcg;
    delete precond;
    delete M;
    delete A;

examples/ex11p_arpack.cpp

@@ -1,381 +0,0 @@
-//                       MFEM Example 11 - Parallel Version
-//
-// Compile with: make ex11p
-//
-// Sample runs:  mpirun -np 4 ex11p -m ../data/square-disc.mesh
-//               mpirun -np 4 ex11p -m ../data/star.mesh
-//               mpirun -np 4 ex11p -m ../data/escher.mesh
-//               mpirun -np 4 ex11p -m ../data/fichera.mesh
-//               mpirun -np 4 ex11p -m ../data/square-disc-p2.vtk -o 2
-//               mpirun -np 4 ex11p -m ../data/square-disc-p3.mesh -o 3
-//               mpirun -np 4 ex11p -m ../data/square-disc-nurbs.mesh -o -1
-//               mpirun -np 4 ex11p -m ../data/disc-nurbs.mesh -o -1 -n 20
-//               mpirun -np 4 ex11p -m ../data/pipe-nurbs.mesh -o -1
-//               mpirun -np 4 ex11p -m ../data/ball-nurbs.mesh -o 2
-//               mpirun -np 4 ex11p -m ../data/star-surf.mesh
-//               mpirun -np 4 ex11p -m ../data/square-disc-surf.mesh
-//               mpirun -np 4 ex11p -m ../data/inline-segment.mesh
-//               mpirun -np 4 ex11p -m ../data/amr-quad.mesh
-//               mpirun -np 4 ex11p -m ../data/amr-hex.mesh
-//               mpirun -np 4 ex11p -m ../data/mobius-strip.mesh -n 8
-//               mpirun -np 4 ex11p -m ../data/klein-bottle.mesh -n 10
-//
-// 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 LOBPCG and ARPACK
-//               eigenvalue solvers together with the BoomerAMG preconditioner
-//               in HYPRE, as well as optionally the SuperLU parallel direct
-//               solver. 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. Initialize MPI.
-   int num_procs, myid;
-   MPI_Init(&argc, &argv);
-   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
-   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
-
-   // 2. Parse command-line options.
-   const char *mesh_file = "../data/star.mesh";
-   int ser_ref_levels = 2;
-   int par_ref_levels = 1;
-   int order = 1;
-   int nev = 5;
-   bool slu_solver = false;
-   bool use_arpack = false;
-   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(&par_ref_levels, "-rp", "--refine-parallel",
-                  "Number of times to refine the mesh uniformly in parallel.");
-   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.");
-#ifdef MFEM_USE_SUPERLU
-   args.AddOption(&slu_solver, "-slu", "--superlu", "-no-slu",
-                  "--no-superlu", "Use the SuperLU Solver.");
-#endif
-#ifdef MFEM_USE_ARPACK
-   args.AddOption(&use_arpack, "-arpack", "--use-arpack", "-no-arpack",
-                  "--no-arpack",
-                  "Enable or disable the use of ARPACK.");
-#endif
-   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
-                  "--no-visualization",
-                  "Enable or disable GLVis visualization.");
-   args.Parse();
-   if (!args.Good())
-   {
-      if (myid == 0)
-      {
-         args.PrintUsage(cout);
-      }
-      MPI_Finalize();
-      return 1;
-   }
-   if (myid == 0)
-   {
-      args.PrintOptions(cout);
-   }
-
-   // 3. 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)
-   {
-      if (myid == 0)
-      {
-         cerr << "\nCan not open mesh file: " << mesh_file << '\n' << endl;
-      }
-      MPI_Finalize();
-      return 2;
-   }
-   mesh = new Mesh(imesh, 1, 1);
-   imesh.close();
-   int dim = mesh->Dimension();
-
-   // 4. 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();
-   }
-
-   // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
-   //    this mesh further in parallel to increase the resolution (1 time by
-   //    default, or specified on the command line with -rp). Once the parallel
-   //    mesh is defined, the serial mesh can be deleted.
-   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
-   delete mesh;
-   for (int lev = 0; lev < par_ref_levels; lev++)
-   {
-      pmesh->UniformRefinement();
-   }
-
-   // 6. Define a parallel finite element space on the parallel 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 (pmesh->GetNodes())
-   {
-      fec = pmesh->GetNodes()->OwnFEC();
-   }
-   else
-   {
-      fec = new H1_FECollection(order = 1, dim);
-   }
-   ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
-   HYPRE_Int size = fespace->GlobalTrueVSize();
-   if (myid == 0)
-   {
-      cout << "Number of unknowns: " << size << endl;
-   }
-
-   // 7. 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 (pmesh->bdr_attributes.Size())
-   {
-      ess_bdr.SetSize(pmesh->bdr_attributes.Max());
-      ess_bdr = 1;
-   }
-
-   ParBilinearForm *a = new ParBilinearForm(fespace);
-   a->AddDomainIntegrator(new DiffusionIntegrator(one));
-   if (pmesh->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();
-   a->EliminateEssentialBCDiag(ess_bdr, 1.0);
-   a->Finalize();
-
-   ParBilinearForm *m = new ParBilinearForm(fespace);
-   m->AddDomainIntegrator(new MassIntegrator(one));
-   m->Assemble();
-   // shift the eigenvalue corresponding to eliminated dofs to a large value
-   m->EliminateEssentialBCDiag(ess_bdr, numeric_limits<double>::min());
-   m->Finalize();
-
-   HypreParMatrix *A = a->ParallelAssemble();
-   HypreParMatrix *M = m->ParallelAssemble();
-
-#ifdef MFEM_USE_SUPERLU
-   Operator * Arow = NULL;
-   if (slu_solver)
-   {
-      Arow = new SuperLURowLocMatrix(*A);
-   }
-#endif
-
-   delete a;
-   delete m;
-
-   // 8. Define and configure the LOBPCG eigensolver and the BoomerAMG
-   //    preconditioner for A to be used within the solver. Set the matrices
-   //    which define the generalized eigenproblem A x = lambda M x.
-   Eigensolver * esolver = NULL;
-   Solver      *  solver = NULL;
-   Solver      * precond = NULL;
-
-   if (!slu_solver)
-   {
-      HypreBoomerAMG * amg = new HypreBoomerAMG(*A);
-      amg->SetPrintLevel(0);
-      precond = amg;
-
-#ifdef MFEM_USE_ARPACK
-      if ( use_arpack )
-      {
-         HyprePCG * pcg = new HyprePCG(*A);
-         pcg->SetTol(1e-12);
-         pcg->SetMaxIter(200);
-         pcg->SetPreconditioner(*amg);
-         pcg->SetPrintLevel(0);
-
-         solver = pcg;
-      }
-#endif
-   }
-#ifdef MFEM_USE_SUPERLU
-   else
-   {
-      SuperLUSolver * superlu = new SuperLUSolver(MPI_COMM_WORLD);
-      superlu->SetPrintStatistics(false);
-      superlu->SetSymmetricPattern(true);
-      superlu->SetColumnPermutation(superlu::PARMETIS);
-      superlu->SetOperator(*Arow);
-
-      solver  = use_arpack?superlu:NULL;
-      precond = use_arpack?NULL:superlu;
-   }
-#endif
-
-   if ( use_arpack )
-   {
-      ParArPackSym * arpack = new ParArPackSym(MPI_COMM_WORLD);
-      arpack->SetMode(3);
-      arpack->SetPrintLevel(2);
-      arpack->SetSolver(*solver);
-
-      esolver = arpack;
-   }
-   else
-   {
-      HypreLOBPCG * lobpcg = new HypreLOBPCG(MPI_COMM_WORLD);
-      lobpcg->SetPreconditioner(*precond);
-      lobpcg->SetPrecondUsageMode(1);
-      lobpcg->SetPrintLevel(1);
-
-      esolver = lobpcg;
-   }
-
-   esolver->SetNumModes(nev);
-   esolver->SetMaxIter(100);
-   esolver->SetTol(1e-8);
-
-   esolver->SetMassMatrix(*M);
-   esolver->SetOperator(*A);
-
-   // 9. 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;
-   esolver->Solve();
-   esolver->GetEigenvalues(eigenvalues);
-
-   if ( myid == 0 && use_arpack )
-   {
-      cout << endl;
-      for (int i=0; i<eigenvalues.Size(); i++)
-      {
-         cout << "Eigenvalue lambda " << eigenvalues[i] << endl;
-      }
-      cout << endl;
-   }
-
-   ParGridFunction x(fespace);
-
-   // 10. 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 << "mesh." << setfill('0') << setw(6) << myid;
-
-      ofstream mesh_ofs(mesh_name.str().c_str());
-      mesh_ofs.precision(8);
-      pmesh->Print(mesh_ofs);
-
-      for (int i=0; i<nev; i++)
-      {
-         // convert eigenvector from HypreParVector to ParGridFunction
-         x.Distribute(esolver->GetEigenvector(i));
-
-         mode_name << "mode_" << setfill('0') << setw(2) << i << "."
-                   << setfill('0') << setw(6) << myid;
-
-         ofstream mode_ofs(mode_name.str().c_str());
-         mode_ofs.precision(8);
-         x.Save(mode_ofs);
-         mode_name.str("");
-      }
-   }
-
-   // 11. 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++)
-      {
-         if ( myid == 0 )
-         {
-            cout << "Eigenmode " << i+1 << '/' << nev
-                 << ", Lambda = " << eigenvalues[i] << endl;
-         }
-
-         // convert eigenvector from HypreParVector to ParGridFunction
-         x.Distribute(esolver->GetEigenvector(i));
-
-         mode_sock << "parallel " << num_procs << " " << myid << "\n"
-                   << "solution\n" << *pmesh << x << flush
-                   << "window_title 'Eigenmode " << i+1 << '/' << nev
-                   << ", Lambda = " << eigenvalues[i] << "'" << endl;
-
-         char c;
-         if (myid == 0)
-         {
-            cout << "press (q)uit or (c)ontinue --> " << flush;
-            cin >> c;
-         }
-         MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
-
-         if (c != 'c')
-         {
-            break;
-         }
-      }
-      mode_sock.close();
-   }
-
-   // 12. Free the used memory.
-   delete esolver;
-   delete solver;
-   delete precond;
-   delete M;
-   delete A;
-
-   delete fespace;
-   if (order > 0)
-   {
-      delete fec;
-   }
-   delete pmesh;
-
-   MPI_Finalize();
-
-   return 0;
-}

examples/ex12p.cpp

@@ -276,8 +276,8 @@ int main(int argc, char *argv[])
 
       for (int i=0; i<nev; i++)
       {
-         // convert eigenvector from Vector to ParGridFunction
-         x.Distribute(lobpcg->GetEigenvector(i));
+         // convert eigenvector from HypreParVector to ParGridFunction
+         x = lobpcg->GetEigenvector(i);
 
          mode_name << "mode_" << setfill('0') << setw(2) << i << "."
                    << setfill('0') << setw(6) << myid;
@@ -303,7 +303,7 @@ int main(int argc, char *argv[])
       pmesh->Print(adios2output);
       for (int i=0; i<nev; i++)
       {
-         x.Distribute(lobpcg->GetEigenvector(i));
+         x = lobpcg->GetEigenvector(i);
          // x is a temporary that must be saved immediately
          x.Save(adios2output, "mode_" + std::to_string(i));
       }
@@ -326,8 +326,8 @@ int main(int argc, char *argv[])
                  << ", Lambda = " << eigenvalues[i] << endl;
          }
 
-         // convert eigenvector from Vector to ParGridFunction
-         x.Distribute(lobpcg->GetEigenvector(i));
+         // convert eigenvector from HypreParVector to ParGridFunction
+         x = lobpcg->GetEigenvector(i);
 
          mode_sock << "parallel " << num_procs << " " << myid << "\n"
                    << "solution\n" << *pmesh << x << flush

examples/ex13.cpp

@@ -1,282 +0,0 @@
-//                       MFEM Example 13
-//
-// Compile with: make ex3p
-//
-// Sample runs:  ex13 -m ../data/star.mesh -s 5
-//               ex13 -m ../data/square-disc.mesh -o 2 -n 4 // minres fails to conv.
-//               ex13 -m ../data/beam-hex.mesh
-//               ex13 -m ../data/square-disc.mesh -rs 1 -s 26
-//               ex13 -m ../data/square-disc-nurbs.mesh -rs 3 -s 26
-//               ex13 -m ../data/amr-quad.mesh -o 2 // minres fails to conv.
-//               ex13 -m ../data/mobius-strip.mesh -n 8
-//
-// Description:  This example code solves a simple 3D electromagnetic
-//               eigenmode problem corresponding to the second order
-//               Maxwell equation curl curl E = lambda E with boundary
-//               condition E x n = 0. We discretize with Nedelec finite
-//               elements.
-//
-//               The example demonstrates the use of H(curl) finite element
-//               spaces with the curl-curl and the (vector finite element) mass
-//               bilinear form, as well as the use of the ARPACK eigenmode
-//               solver for symmetric matrices using the shift-invert mode.
-//
-
-#include "mfem.hpp"
-#include <fstream>
-#include <iostream>
-
-using namespace std;
-using namespace mfem;
-
-#ifdef MFEM_USE_ARPACK
-
-int main(int argc, char *argv[])
-{
-   // 1. Parse command-line options.
-   const char *mesh_file = "../data/beam-tet.mesh";
-   int order = 1;
-   int nev = 5;
-   int sr = 2;
-   double sigma = 11.0;
-   bool visualization = 1;
-   bool arp_solver = true;
-
-   OptionsParser args(argc, argv);
-   args.AddOption(&mesh_file, "-m", "--mesh",
-                  "Mesh file to use.");
-   args.AddOption(&order, "-o", "--order",
-                  "Finite element order (polynomial degree).");
-   args.AddOption(&nev, "-n", "--num-eigs",
-                  "Number of desired eigenmodes.");
-   args.AddOption(&sr, "-rs", "--refine-serial",
-                  "Number of times to refine the mesh uniformly in serial.");
-   args.AddOption(&sigma, "-s", "--shift",
-                  "Average of the desired eigenvalue range.");
-   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 mesh from the given mesh file.  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 mesh to increase the resolution. In this example we do
-   //    'ref_levels' of uniform refinement.
-   {
-      int ref_levels = sr;
-      for (int l = 0; l < ref_levels; l++)
-      {
-         mesh->UniformRefinement();
-      }
-   }
-
-   // 4. Define a finite element space on the mesh. Here we use the lowest
-   //    order Nedelec finite elements, but we can easily switch
-   //    to higher-order spaces by changing the value of p.
-   FiniteElementCollection *fec = new ND_FECollection(order, dim);
-   FiniteElementSpace *fespace = new FiniteElementSpace(mesh, fec);
-   int size = fespace->GetVSize();
-
-   cout << "Number of unknowns: " << size << endl;
-   cout << "Number of boundary attributes: " << mesh->bdr_attributes.Max()
-        << endl;
-
-   // 5. Set up the parallel bilinear form corresponding to the EM diffusion
-   //    operator curl muinv curl - sigma I, by adding the curl-curl and the
-   //    mass domain integrators and finally imposing homogeneous Dirichlet
-   //    boundary conditions. The boundary conditions are implemented by
-   //    marking all the boundary attributes from the mesh as essential
-   //    (Dirichlet). After serial and parallel assembly we extract the
-   //    parallel matrices A and M.
-   Coefficient *muinv    = new ConstantCoefficient(1.0);
-   Coefficient *negSigma = new ConstantCoefficient(-sigma);
-
-   BilinearForm *a = new BilinearForm(fespace);
-   a->AddDomainIntegrator(new CurlCurlIntegrator(*muinv));
-   a->AddDomainIntegrator(new VectorFEMassIntegrator(*negSigma));
-   a->Assemble();
-   Array<int> ess_bdr(mesh->bdr_attributes.Max());
-   ess_bdr = 1;
-   a->EliminateEssentialBC(ess_bdr);
-   a->Finalize();
-
-   BilinearForm *m = new BilinearForm(fespace);
-   m->AddDomainIntegrator(new VectorFEMassIntegrator());
-   m->Assemble();
-   m->EliminateEssentialBCDiag(ess_bdr, sqrt(numeric_limits<double>::min()));
-   m->Finalize();
-
-   // 6. Define a parallel grid function to approximate each of the
-   // eigenmodes returned by the solver.  Use this as a template to
-   // create a special multi-vector object needed by the eigensolver
-   // which is then initialized with random values.
-   GridFunction x(fespace);
-   x = 0.0;
-
-   // 7. Define and configure the GMRES
-   //    solver to be used within the eigensolver.
-   Solver    * solver = NULL;
-   if ( false )
-   {
-      GMRESSolver *  gmres = new GMRESSolver();
-
-      gmres->SetOperator(*a);
-      gmres->SetRelTol(1e-8);
-      gmres->SetMaxIter(1000);
-      gmres->SetPrintLevel(0);
-      solver = gmres;
-   }
-   else
-   {
-#ifndef MFEM_USE_SUITESPARSE
-      cout << "Building MINRESSolver" << endl;
-      MINRESSolver * minres = new MINRESSolver();
-
-      minres->SetRelTol(1e-12);
-      minres->SetMaxIter(1000);
-      minres->SetPrintLevel(0);
-      solver = minres;
-#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(a->SpMat());
-
-   // 7. Define and configure the ARPACK eigensolver
-   SymGenEigensolver * eig_solver = NULL;
-   if (arp_solver)
-   {
-      ArPackSAUPD * arpack = new ArPackSAUPD();
-
-      arpack->SetNumModes(nev);
-      arpack->SetMaxIter(400);
-      arpack->SetTol(1e-8);
-      arpack->SetShift(sigma);
-      arpack->SetMode(3);
-      arpack->SetPrintLevel(2);
-      arpack->SetSolver(*solver);
-
-      eig_solver = arpack;
-   }
-
-   eig_solver->SetOperators(*a, *m);
-
-   // Obtain the eigenvalues and eigenvectors
-   Array<double> eigenvalues(nev);
-   eigenvalues = -1.0;
-
-   // arpack->Solve(eigenvalues, *eigenvectors);
-   eig_solver->Solve();
-
-   eig_solver->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<min(nev,eigenvalues.Size()); i++)
-   {
-      cout << "Eigenvalue lambda   " << eigenvalues[i] << endl;
-   }
-   cout.precision(old_prec);
-   cout.flags(old_fmt);
-   cout << endl;
-
-   VisItDataCollection visit_dc("Example13", mesh);
-   GridFunction ** mode = new GridFunction*[min(nev,eigenvalues.Size())];
-   for (int i=0; i<min(nev,eigenvalues.Size()); i++)
-   {
-      mode[i] = new GridFunction(fespace);
-      *mode[i] = eig_solver->GetEigenvector(i);
-
-      ostringstream modeName;
-      modeName << "mode_" << setfill('0') << setw(2) << i;
-      visit_dc.RegisterField(modeName.str().c_str(),mode[i]);
-   }
-   visit_dc.Save();
-
-   // 8. Save the refined mesh and the modes. This output can
-   //    be viewed later using GLVis: "glvis -m mesh -g mode".
-   {
-      ofstream mesh_ofs("refined.mesh");
-      mesh_ofs.precision(8);
-      mesh->Print(mesh_ofs);
-
-      for (int i=0; i<min(nev,eigenvalues.Size()); i++)
-      {
-         x = eig_solver->GetEigenvector(i);
-
-         ostringstream modeName;
-         modeName << "mode_" << setfill('0') << setw(2) << i;
-
-         ofstream mode_ofs(modeName.str().c_str());
-         mode_ofs.precision(8);
-         x.Save(mode_ofs);
-         modeName.str("");
-      }
-   }
-
-   // 9. 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<min(nev,eigenvalues.Size()); i++)
-      {
-         x = eig_solver->GetEigenvector(i);
-
-         mode_sock << "solution\n" << *mesh << x << flush;
-
-         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 a;
-   delete m;
-   delete negSigma;
-   delete muinv;
-   delete eig_solver;
-   delete solver;
-   // delete X;
-   delete fespace;
-   delete fec;
-   delete mesh;
-
-   return 0;
-}
-
-#endif // MFEM_USE_ARPACK

examples/ex13p.cpp

@@ -215,8 +215,8 @@ int main(int argc, char *argv[])
 
       for (int i=0; i<nev; i++)
       {
-         // convert eigenvector from Vector to ParGridFunction
-         x.Distribute(ame->GetEigenvector(i));
+         // convert eigenvector from HypreParVector to ParGridFunction
+         x = ame->GetEigenvector(i);
 
          mode_name << "mode_" << setfill('0') << setw(2) << i << "."
                    << setfill('0') << setw(6) << myid;
@@ -244,8 +244,8 @@ int main(int argc, char *argv[])
                  << ", Lambda = " << eigenvalues[i] << endl;
          }
 
-         // convert eigenvector from Vector to ParGridFunction
-         x.Distribute(ame->GetEigenvector(i));
+         // convert eigenvector from HypreParVector to ParGridFunction
+         x = ame->GetEigenvector(i);
 
          mode_sock << "parallel " << num_procs << " " << myid << "\n"
                    << "solution\n" << *pmesh << x << flush

examples/ex32p.cpp

@@ -228,7 +228,7 @@ int main(int argc, char *argv[])
       for (int i=0; i<nev; i++)
       {
          // convert eigenvector from HypreParVector to ParGridFunction
-         x.Distribute(ame->GetEigenvector(i));
+         x = ame->GetEigenvector(i);
          curl.Mult(x, dx);
 
          mode_name << "mode_" << setfill('0') << setw(2) << i << "."
@@ -295,7 +295,7 @@ int main(int argc, char *argv[])
             }
 
             // convert eigenvector from HypreParVector to ParGridFunction
-            x.Distribute(ame->GetEigenvector(i));
+            x = ame->GetEigenvector(i);
             curl.Mult(x, dx);
 
             {
@@ -469,7 +469,7 @@ int main(int argc, char *argv[])
             }
 
             // convert eigenvector from HypreParVector to ParGridFunction
-            x.Distribute(ame->GetEigenvector(i));
+            x = ame->GetEigenvector(i);
             curl.Mult(x, dx);
 
             {
@@ -599,7 +599,7 @@ int main(int argc, char *argv[])
             }
 
             // convert eigenvector from HypreParVector to ParGridFunction
-            x.Distribute(ame->GetEigenvector(i));
+            x = ame->GetEigenvector(i);
             curl.Mult(x, dx);
 
             mode_sock << "parallel " << num_procs << " " << myid << "\n"

examples/ex35p.cpp

@@ -658,7 +658,7 @@ void ScalarWaveGuide(int mode, ParGridFunction &x)
    lobpcg.SetOperator(*A);
    lobpcg.Solve();
 
-   x.Distribute(lobpcg.GetEigenvector(mode));
+   x = lobpcg.GetEigenvector(mode);
 
    delete A;
    delete M;
@@ -714,7 +714,7 @@ void VectorWaveGuide(int mode, ParGridFunction &x)
    ame.SetOperator(*A);
    ame.Solve();
 
-   x.Distribute(ame.GetEigenvector(mode));
+   x = ame.GetEigenvector(mode);
 
    delete A;
    delete M;
@@ -780,7 +780,7 @@ void PseudoScalarWaveGuide(int mode, ParGridFunction &x_l2)
    lobpcg.SetOperator(*A);
    lobpcg.Solve();
 
-   x.Distribute(lobpcg.GetEigenvector(mode));
+   x = lobpcg.GetEigenvector(mode);
 
    x_l2.ProjectCoefficient(xCoef);
 

examples/ex36.cpp

@@ -97,13 +97,7 @@ int main(int argc, char *argv[])
    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);
+   args.ParseCheck();
 
    // 2. Read the mesh from the mesh file.
    const char *mesh_file = "../data/disc-nurbs.mesh";
@@ -128,7 +122,7 @@ int main(int argc, char *argv[])
    *nodes /= scale;
 
    // 4. Define the necessary finite element spaces on the mesh.
-   H1_FECollection H1fec(order+1, dim);
+   H1Bubble_FECollection H1fec(order, order - 1, dim);
    FiniteElementSpace H1fes(&mesh, &H1fec);
 
    L2_FECollection L2fec(order-1, dim);

examples/ex36p.cpp

@@ -103,19 +103,7 @@ int main(int argc, char *argv[])
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
-   args.Parse();
-   if (!args.Good())
-   {
-      if (myid == 0)
-      {
-         args.PrintUsage(cout);
-      }
-      return 1;
-   }
-   if (myid == 0)
-   {
-      args.PrintOptions(cout);
-   }
+   args.ParseCheck();
 
    // 2. Read the mesh from the mesh file.
    const char *mesh_file = "../data/disc-nurbs.mesh";
@@ -143,7 +131,7 @@ int main(int argc, char *argv[])
    mesh.Clear();
 
    // 4. Define the necessary finite element spaces on the mesh.
-   H1_FECollection H1fec(order+1, dim);
+   H1Bubble_FECollection H1fec(order, order - 1, dim);
    ParFiniteElementSpace H1fes(&pmesh, &H1fec);
 
    L2_FECollection L2fec(order-1, dim);

examples/ex37.cpp

@@ -5,8 +5,8 @@
 // Sample runs:
 //     ex37 -alpha 10
 //     ex37 -alpha 10 -pv
-//     ex37 -lambda 0.1 -mu 0.1 -growth 1
-//     ex37 -o 2 -alpha 10.0 -mi 50 -vf 0.4 -ntol 1e-5 -growth 1.5
+//     ex37 -lambda 0.1 -mu 0.1
+//     ex37 -o 2 -alpha 5.0 -mi 50 -vf 0.4 -ntol 1e-5
 //     ex37 -r 6 -o 1 -alpha 25.0 -epsilon 0.02 -mi 50 -ntol 1e-5
 //
 // Description: This example code demonstrates the use of MFEM to solve a
@@ -55,6 +55,53 @@
 using namespace std;
 using namespace mfem;
 
+/**
+ * @brief Bregman projection of ρ = sigmoid(ψ) onto the subspace
+ *        ∫_Ω ρ dx = θ vol(Ω) as follows:
+ *
+ *        1. Compute the root of the R → R function
+ *            f(c) = ∫_Ω sigmoid(ψ + c) dx - θ vol(Ω)
+ *        2. Set ψ ← ψ + c.
+ *
+ * @param psi a GridFunction to be updated
+ * @param target_volume θ vol(Ω)
+ * @param tol Newton iteration tolerance
+ * @param max_its Newton maximum iteration number
+ * @return real_t Final volume, ∫_Ω sigmoid(ψ)
+ */
+real_t proj(GridFunction &psi, real_t target_volume, real_t tol=1e-12,
+            int max_its=10)
+{
+   MappedGridFunctionCoefficient sigmoid_psi(&psi, sigmoid);
+   MappedGridFunctionCoefficient der_sigmoid_psi(&psi, der_sigmoid);
+
+   LinearForm int_sigmoid_psi(psi.FESpace());
+   int_sigmoid_psi.AddDomainIntegrator(new DomainLFIntegrator(sigmoid_psi));
+   LinearForm int_der_sigmoid_psi(psi.FESpace());
+   int_der_sigmoid_psi.AddDomainIntegrator(new DomainLFIntegrator(
+                                              der_sigmoid_psi));
+   bool done = false;
+   for (int k=0; k<max_its; k++) // Newton iteration
+   {
+      int_sigmoid_psi.Assemble(); // Recompute f(c) with updated ψ
+      const real_t f = int_sigmoid_psi.Sum() - target_volume;
+
+      int_der_sigmoid_psi.Assemble(); // Recompute df(c) with updated ψ
+      const real_t df = int_der_sigmoid_psi.Sum();
+
+      const real_t dc = -f/df;
+      psi += dc;
+      if (abs(dc) < tol) { done = true; break; }
+   }
+   if (!done)
+   {
+      mfem_warning("Projection reached maximum iteration without converging. "
+                   "Result may not be accurate.");
+   }
+   int_sigmoid_psi.Assemble();
+   return int_sigmoid_psi.Sum();
+}
+
 /*
  * ---------------------------------------------------------------
  *                      ALGORITHM PREAMBLE
@@ -133,11 +180,10 @@ int main(int argc, char *argv[])
    int ref_levels = 5;
    int order = 2;
    real_t alpha = 1.0;
-   real_t growth = 2;
    real_t epsilon = 0.01;
    real_t vol_fraction = 0.5;
    int max_it = 1e3;
-   real_t itol = 1e-2;
+   real_t itol = 1e-1;
    real_t ntol = 1e-4;
    real_t rho_min = 1e-6;
    real_t lambda = 1.0;
@@ -152,8 +198,6 @@ int main(int argc, char *argv[])
                   "Order (degree) of the finite elements.");
    args.AddOption(&alpha, "-alpha", "--alpha-step-length",
                   "Step length for gradient descent.");
-   args.AddOption(&growth, "-growth", "--alpha-growth-rate",
-                  "Growth rate of step length for gradient descent.");
    args.AddOption(&epsilon, "-epsilon", "--epsilon-thickness",
                   "Length scale for ρ.");
    args.AddOption(&max_it, "-mi", "--max-it",
@@ -288,7 +332,6 @@ int main(int argc, char *argv[])
    }
    FilterSolver->SetEssentialBoundary(ess_bdr_filter);
    FilterSolver->SetupFEM();
-   FilterSolver->AssembleDiffusionBilinear();
 
    BilinearForm mass(&control_fes);
    mass.AddDomainIntegrator(new InverseIntegrator(new MassIntegrator(one)));
@@ -342,7 +385,7 @@ int main(int argc, char *argv[])
    // 11. Iterate:
    for (int k = 1; k <= max_it; k++)
    {
-      if (k > 1) { alpha = std::pow((real_t) k,growth); }
+      if (k > 1) { alpha *= ((real_t) k) / ((real_t) k-1); }
 
       mfem::out << "\nStep = " << k << std::endl;
 
@@ -379,9 +422,7 @@ int main(int argc, char *argv[])
 
       // Step 5 - Update design variable ψ ← proj(ψ - αG)
       psi.Add(-alpha, grad);
-      GridFunction alpha_grad(grad);
-      alpha_grad *= alpha;
-      const real_t material_volume = proj(psi, alpha_grad, target_volume);
+      const real_t material_volume = proj(psi, target_volume);
 
       // Compute ||ρ - ρ_old|| in control fes.
       real_t norm_increment = zerogf.ComputeL1Error(succ_diff_rho);

examples/ex37.hpp

@@ -137,7 +137,7 @@ public:
         exponent(exponent_), rho_min(rho_min_)
    {
       MFEM_ASSERT(rho_min_ >= 0.0, "rho_min must be >= 0");
-      MFEM_ASSERT(rho_min_ < 1.0,  "rho_min must be < 1");
+      MFEM_ASSERT(rho_min_ < 1.0,  "rho_min must be > 1");
       MFEM_ASSERT(u, "displacement field is not set");
       MFEM_ASSERT(rho_filter, "density field is not set");
    }
@@ -231,12 +231,9 @@ private:
    FiniteElementCollection * fec = nullptr;
    FiniteElementSpace * fes = nullptr;
    Array<int> ess_bdr;
-   Array<int> ess_tdof_list;
    Array<int> neumann_bdr;
    GridFunction * u = nullptr;
    LinearForm * b = nullptr;
-   BilinearForm * a = nullptr;
-   OperatorPtr A;
    bool parallel;
 #ifdef MFEM_USE_MPI
    ParMesh * pmesh = nullptr;
@@ -270,8 +267,6 @@ public:
    void ResetFEM();
    void SetupFEM();
 
-   void UpdateEssentialTDofs();
-   void AssembleDiffusionBilinear(bool update_ess_tdofs=true);
    void Solve();
    GridFunction * GetFEMSolution();
    LinearForm * GetLinearForm() {return b;}
@@ -376,130 +371,6 @@ public:
 
 };
 
-/**
- * @brief Bregman projection of ρ = sigmoid(ψ) onto the subspace
- *        ∫_Ω ρ dx = θ vol(Ω) as follows:
- *
- *        1. Compute the root of the R → R function
- *            f(c) = ∫_Ω sigmoid(ψ + c) dx - θ vol(Ω)
- *           using the Illinois method
- *        2. Set ψ ← ψ + c.
- *
- * @param psi a GridFunction to be updated
- * @param alpha_grad alpha multiplied by gradient
- * @param target_volume θ vol(Ω)
- * @param tol Illinois iteration tolerance
- * @param max_its Illinois maximum iteration number
- * @return real_t Final volume (∫_Ω sigmoid(ψ) dx)
- */
-real_t proj(GridFunction &psi, GridFunction &alpha_grad, real_t target_volume,
-            real_t tol = 1e-12, int max_its = 100)
-{
-#ifdef MFEM_USE_MPI
-   FiniteElementSpace *fes = psi.FESpace();
-   ParFiniteElementSpace *pfes = dynamic_cast<ParFiniteElementSpace*>(fes);
-#endif
-   ConstantCoefficient zero_cf(0.0);
-   real_t a = -alpha_grad.ComputeMaxError(zero_cf);
-   real_t b = -a;
-   real_t y = 0.0;
-
-   MappedGridFunctionCoefficient sigmoid_psi(
-   &psi, [&y](const real_t x) { return sigmoid(x + y); });
-   std::unique_ptr<LinearForm> int_sigmoid_psi;
-#ifdef MFEM_USE_MPI
-   ParGridFunction *par_psi = dynamic_cast<ParGridFunction *>(&psi);
-   if (par_psi)
-   {
-      int_sigmoid_psi.reset(new ParLinearForm(par_psi->ParFESpace()));
-   }
-   else
-   {
-      int_sigmoid_psi.reset(new LinearForm(psi.FESpace()));
-   }
-#else
-   int_sigmoid_psi.reset(new LinearForm(psi.FESpace()));
-#endif
-   int_sigmoid_psi->AddDomainIntegrator(new DomainLFIntegrator(sigmoid_psi));
-
-   y = a;
-   int_sigmoid_psi->Assemble();
-   real_t f_a = int_sigmoid_psi->Sum();  // f_a := f(a) + θ vol(Ω)
-
-   y = b;
-   int_sigmoid_psi->Assemble();
-   real_t f_b = int_sigmoid_psi->Sum();  // f_b := f(b) + θ vol(Ω)
-#ifdef MFEM_USE_MPI
-   if (pfes)
-   {
-      MPI_Allreduce(MPI_IN_PLACE, &f_a, 1, MPITypeMap<real_t>::mpi_type,
-                    MPI_SUM, MPI_COMM_WORLD);
-      MPI_Allreduce(MPI_IN_PLACE, &f_b, 1, MPITypeMap<real_t>::mpi_type,
-                    MPI_SUM, MPI_COMM_WORLD);
-   }
-#endif
-   f_a -= target_volume; // f_a := f(a)
-   f_b -= target_volume; // f_b := f(b)
-   real_t c = 0.0;
-   real_t f_c = 0.0;
-   int side = 0;
-
-   bool done = false;
-   for (int k=0; k < max_its; k++)
-   {
-      c = (f_a * b - f_b * a) / (f_a - f_b);
-
-      if (abs(b - a) < tol * abs(b + a)) { done = true; break; }
-
-      y = c;
-      int_sigmoid_psi->Assemble();
-      f_c = int_sigmoid_psi->Sum(); // f_c := f(c) + θ vol(Ω)
-#ifdef MFEM_USE_MPI
-      if (pfes)
-      {
-         MPI_Allreduce(MPI_IN_PLACE, &f_c, 1, MPITypeMap<real_t>::mpi_type,
-                       MPI_SUM, MPI_COMM_WORLD);
-      }
-#endif
-      f_c -= target_volume; // f_c := f(c)
-
-      if (f_c * f_b > 0)
-      {
-         b = c;
-         f_b = f_c;
-         if (side == -1) { f_a /= 2.0; }
-         side = -1;
-      }
-      else if (f_c * f_a > 0)
-      {
-         a = c;
-         f_a = f_c;
-         if (side == 1) { f_b /= 2.0; }
-         side = 1;
-      }
-      else
-      {
-         done = true; break;
-      }
-   }
-   if (!done)
-   {
-      mfem_warning("Projection reached maximum iteration without converging. "
-                   "Result may not be accurate.");
-   }
-   y = 0.0;
-   psi += c;
-   int_sigmoid_psi->Assemble();
-   real_t material_volume = int_sigmoid_psi->Sum();
-#ifdef MFEM_USE_MPI
-   if (pfes)
-   {
-      MPI_Allreduce(MPI_IN_PLACE, &material_volume, 1,
-                    MPITypeMap<real_t>::mpi_type, MPI_SUM, MPI_COMM_WORLD);
-   }
-#endif
-   return material_volume;
-}
 
 // Poisson solver
 
@@ -551,8 +422,12 @@ void DiffusionSolver::SetupFEM()
    }
 }
 
-void DiffusionSolver::UpdateEssentialTDofs()
+void DiffusionSolver::Solve()
 {
+   OperatorPtr A;
+   Vector B, X;
+   Array<int> ess_tdof_list;
+
 #ifdef MFEM_USE_MPI
    if (parallel)
    {
@@ -565,39 +440,7 @@ void DiffusionSolver::UpdateEssentialTDofs()
 #else
    fes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list);
 #endif
-}
-
-void DiffusionSolver::AssembleDiffusionBilinear(bool update_ess_tdofs)
-{
-   if (update_ess_tdofs)
-   {
-      UpdateEssentialTDofs();
-   }
-#ifdef MFEM_USE_MPI
-   if (parallel)
-   {
-      a = new ParBilinearForm(pfes);
-   }
-   else
-   {
-      a = new BilinearForm(fes);
-   }
-#else
-   a = new BilinearForm(fes);
-#endif
-   a->AddDomainIntegrator(new DiffusionIntegrator(*diffcf));
-   if (masscf)
-   {
-      a->AddDomainIntegrator(new MassIntegrator(*masscf));
-   }
-   a->Assemble();
-   a->FormSystemMatrix(ess_tdof_list, A);
-}
-
-void DiffusionSolver::Solve()
-{
-   Vector B, X;
-
+   *u=0.0;
    if (b)
    {
       delete b;
@@ -632,33 +475,31 @@ void DiffusionSolver::Solve()
 
    b->Assemble();
 
-   *u=0.0;
-   if (essbdr_cf)
-   {
-      u->ProjectBdrCoefficient(*essbdr_cf,ess_bdr);
-   }
+   BilinearForm * a = nullptr;
 
 #ifdef MFEM_USE_MPI
    if (parallel)
    {
-      X.SetSize(pfes->TrueVSize());
-      B.SetSize(pfes->TrueVSize());
-      dynamic_cast<ParGridFunction*>(u)->ParallelAssemble(X);
-      dynamic_cast<ParLinearForm*>(b)->ParallelAssemble(B);
-      dynamic_cast<ParBilinearForm*>(a)->ParallelEliminateTDofsInRHS(
-         ess_tdof_list, X, B);
+      a = new ParBilinearForm(pfes);
    }
    else
    {
-      X.NewDataAndSize(u->GetData(), u->Size());
-      B.NewDataAndSize(b->GetData(), b->Size());
-      a->EliminateVDofsInRHS(ess_tdof_list, X, B);
+      a = new BilinearForm(fes);
    }
 #else
-   X.NewDataAndSize(u->GetData(), u->Size());
-   B.NewDataAndSize(b->GetData(), b->Size());
-   a->EliminateVDofsInRHS(ess_tdof_list, X, B);
+   a = new BilinearForm(fes);
 #endif
+   a->AddDomainIntegrator(new DiffusionIntegrator(*diffcf));
+   if (masscf)
+   {
+      a->AddDomainIntegrator(new MassIntegrator(*masscf));
+   }
+   a->Assemble();
+   if (essbdr_cf)
+   {
+      u->ProjectBdrCoefficient(*essbdr_cf,ess_bdr);
+   }
+   a->FormLinearSystem(ess_tdof_list, *u, *b, A, X, B);
 
    CGSolver * cg = nullptr;
    Solver * M = nullptr;
@@ -687,6 +528,7 @@ void DiffusionSolver::Solve()
    delete M;
    delete cg;
    a->RecoverFEMSolution(X, *b, *u);
+   delete a;
 }
 
 GridFunction * DiffusionSolver::GetFEMSolution()
@@ -718,8 +560,6 @@ DiffusionSolver::~DiffusionSolver()
 #endif
    delete fec; fec = nullptr;
    delete b;
-   A.Clear();
-   delete a;
 }
 
 

examples/ex37p.cpp

@@ -4,8 +4,8 @@
 //
 // Sample runs:
 //    mpirun -np 4 ex37p -alpha 10 -pv
-//    mpirun -np 4 ex37p -lambda 0.1 -mu 0.1 -growth 1
-//    mpirun -np 4 ex37p -o 2 -alpha 10.0 -mi 50 -vf 0.4 -ntol 1e-5 -growth 1.5
+//    mpirun -np 4 ex37p -lambda 0.1 -mu 0.1
+//    mpirun -np 4 ex37p -o 2 -alpha 5.0 -mi 50 -vf 0.4 -ntol 1e-5
 //    mpirun -np 4 ex37p -r 6 -o 2 -alpha 10.0 -epsilon 0.02 -mi 50 -ntol 1e-5
 //
 // Description: This example code demonstrates the use of MFEM to solve a
@@ -54,6 +54,61 @@
 using namespace std;
 using namespace mfem;
 
+/**
+ * @brief Bregman projection of ρ = sigmoid(ψ) onto the subspace
+ *        ∫_Ω ρ dx = θ vol(Ω) as follows:
+ *
+ *        1. Compute the root of the R → R function
+ *            f(c) = ∫_Ω sigmoid(ψ + c) dx - θ vol(Ω)
+ *        2. Set ψ ← ψ + c.
+ *
+ * @param psi a GridFunction to be updated
+ * @param target_volume θ vol(Ω)
+ * @param tol Newton iteration tolerance
+ * @param max_its Newton maximum iteration number
+ * @return real_t Final volume, ∫_Ω sigmoid(ψ)
+ */
+real_t proj(ParGridFunction &psi, real_t target_volume, real_t tol=1e-12,
+            int max_its=10)
+{
+   MappedGridFunctionCoefficient sigmoid_psi(&psi, sigmoid);
+   MappedGridFunctionCoefficient der_sigmoid_psi(&psi, der_sigmoid);
+
+   ParLinearForm int_sigmoid_psi(psi.ParFESpace());
+   int_sigmoid_psi.AddDomainIntegrator(new DomainLFIntegrator(sigmoid_psi));
+   ParLinearForm int_der_sigmoid_psi(psi.ParFESpace());
+   int_der_sigmoid_psi.AddDomainIntegrator(new DomainLFIntegrator(
+                                              der_sigmoid_psi));
+   bool done = false;
+   for (int k=0; k<max_its; k++) // Newton iteration
+   {
+      int_sigmoid_psi.Assemble(); // Recompute f(c) with updated ψ
+      real_t f = int_sigmoid_psi.Sum();
+      MPI_Allreduce(MPI_IN_PLACE, &f, 1, MPITypeMap<real_t>::mpi_type,
+                    MPI_SUM, MPI_COMM_WORLD);
+      f -= target_volume;
+
+      int_der_sigmoid_psi.Assemble(); // Recompute df(c) with updated ψ
+      real_t df = int_der_sigmoid_psi.Sum();
+      MPI_Allreduce(MPI_IN_PLACE, &df, 1, MPITypeMap<real_t>::mpi_type,
+                    MPI_SUM, MPI_COMM_WORLD);
+
+      const real_t dc = -f/df;
+      psi += dc;
+      if (abs(dc) < tol) { done = true; break; }
+   }
+   if (!done)
+   {
+      mfem_warning("Projection reached maximum iteration without converging. "
+                   "Result may not be accurate.");
+   }
+   int_sigmoid_psi.Assemble();
+   real_t material_volume = int_sigmoid_psi.Sum();
+   MPI_Allreduce(MPI_IN_PLACE, &material_volume, 1,
+                 MPITypeMap<real_t>::mpi_type, MPI_SUM, MPI_COMM_WORLD);
+   return material_volume;
+}
+
 /*
  * ---------------------------------------------------------------
  *                      ALGORITHM PREAMBLE
@@ -138,11 +193,10 @@ int main(int argc, char *argv[])
    int ref_levels = 5;
    int order = 2;
    real_t alpha = 1.0;
-   real_t growth = 2;
    real_t epsilon = 0.01;
    real_t vol_fraction = 0.5;
    int max_it = 1e3;
-   real_t itol = 1e-2;
+   real_t itol = 1e-1;
    real_t ntol = 1e-4;
    real_t rho_min = 1e-6;
    real_t lambda = 1.0;
@@ -157,8 +211,6 @@ int main(int argc, char *argv[])
                   "Order (degree) of the finite elements.");
    args.AddOption(&alpha, "-alpha", "--alpha-step-length",
                   "Step length for gradient descent.");
-   args.AddOption(&growth, "-growth", "--alpha-growth-rate",
-                  "Growth rate of step length for gradient descent.");
    args.AddOption(&epsilon, "-epsilon", "--epsilon-thickness",
                   "Length scale for ρ.");
    args.AddOption(&max_it, "-mi", "--max-it",
@@ -307,7 +359,6 @@ int main(int argc, char *argv[])
    }
    FilterSolver->SetEssentialBoundary(ess_bdr_filter);
    FilterSolver->SetupFEM();
-   FilterSolver->AssembleDiffusionBilinear();
 
    ParBilinearForm mass(&control_fes);
    mass.AddDomainIntegrator(new InverseIntegrator(new MassIntegrator(one)));
@@ -361,7 +412,7 @@ int main(int argc, char *argv[])
    // 11. Iterate:
    for (int k = 1; k <= max_it; k++)
    {
-      if (k > 1) { alpha = std::pow((real_t) k,growth); }
+      if (k > 1) { alpha *= ((real_t) k) / ((real_t) k-1); }
 
       if (myid == 0)
       {
@@ -401,9 +452,7 @@ int main(int argc, char *argv[])
 
       // Step 5 - Update design variable ψ ← proj(ψ - αG)
       psi.Add(-alpha, grad);
-      ParGridFunction alpha_grad(grad);
-      alpha_grad *= alpha;
-      const real_t material_volume = proj(psi, alpha_grad, target_volume);
+      const real_t material_volume = proj(psi, target_volume);
 
       // Compute ||ρ - ρ_old|| in control fes.
       real_t norm_increment = zerogf.ComputeL1Error(succ_diff_rho);

examples/makefile

@@ -31,9 +31,6 @@ SEQ_DEVICE_EXAMPLES = ex1 ex3 ex4 ex5 ex6 ex9 ex14 ex22 ex24 ex25 ex26 ex34
 PAR_DEVICE_EXAMPLES = ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex9p ex13p ex14p \
    ex22p ex24p ex25p ex26p ex34p ex35p
 
-ifeq ($(MFEM_USE_ARPACK),YES)
-   SEQ_EXAMPLES += ex11 ex13
-endif
 ifeq ($(MFEM_USE_LAPACK),YES)
    SEQ_EXAMPLES += ex38
 endif
@@ -160,8 +157,6 @@ ex37-test-seq: ex37
 	@$(call mfem-test,$<,, Serial example,-mi 3)
 ex37p-test-par: ex37p
 	@$(call mfem-test,$<, $(RUN_MPI), Parallel example,-mi 3)
-ex39-test-seq: ex39
-	@$(call mfem-test,$<,, Serial example,-m ../data/compass.mesh)
 ex41-test-seq: ex41
 	@$(call mfem-test,$<,, Serial example,-tf 1.0)
 ex41p-test-par: ex41p

fem/CMakeLists.txt

@@ -73,6 +73,7 @@ set(SRCS
   fe/fe_base.cpp
   fe/fe_fixed_order.cpp
   fe/fe_h1.cpp
+  fe/fe_h1_bubble.cpp
   fe/fe_l2.cpp
   fe/fe_nd.cpp
   fe/fe_nurbs.cpp
@@ -133,7 +134,7 @@ set(SRCS
   tmop/assemble/diag2.cpp
   tmop/assemble/grad2_limit.cpp
   tmop/assemble/grad2.cpp
-  tmop/assemble/diag3_limit.cpp  
+  tmop/assemble/diag3_limit.cpp
   tmop/assemble/diag3.cpp
   tmop/assemble/grad3_limit.cpp
   tmop/assemble/grad3.cpp
@@ -221,6 +222,7 @@ set(HDRS
   fe/fe_base.hpp
   fe/fe_fixed_order.hpp
   fe/fe_h1.hpp
+  fe/fe_h1_bubble.hpp
   fe/fe_l2.hpp
   fe/fe_nd.hpp
   fe/fe_nurbs.hpp

fem/coefficient.hpp

@@ -52,9 +52,6 @@ public:
    /// Get the time for time dependent coefficients
    real_t GetTime() { return time; }
 
-   /// Returns dimension of the vector.
-   int GetVDim() { return 1; }
-
    /** @brief Evaluate the coefficient in the element described by @a T at the
        point @a ip. */
    /** @note When this method is called, the caller must make sure that the

fem/datacollection.cpp

@@ -492,8 +492,6 @@ void VisItDataCollection::SaveRootFile()
                            to_padded_string(cycle, pad_digits_cycle) +
                            ".mfem_root";
    std::ofstream root_file(root_name);
-   MFEM_VERIFY(root_file.is_open(),
-               "Failed to open ofstream " << root_name);
    root_file << GetVisItRootString();
    if (!root_file)
    {
@@ -979,10 +977,7 @@ void ParaViewDataCollection::Save()
    // Save the local part of the mesh and grid functions fields to the local
    // VTU file. Also save coefficient fields.
    {
-      std::string os_str = vtu_prefix + GenerateVTUFileName("proc", myid);
-      std::ofstream os(os_str);
-      MFEM_VERIFY(os.is_open(),
-                  "Failed to open ofstream " << os_str);
+      std::ofstream os(vtu_prefix + GenerateVTUFileName("proc", myid));
       os.precision(precision);
       SaveDataVTU(os, levels_of_detail);
    }
@@ -994,10 +989,7 @@ void ParaViewDataCollection::Save()
                   "QuadratureFunction output is not supported for "
                   "ParaViewDataCollection on domain boundary!");
       const std::string &field_name = qfield.first;
-      std::string os_str = vtu_prefix + GenerateVTUFileName(field_name, myid);
-      std::ofstream os(os_str);
-      MFEM_VERIFY(os.is_open(),
-                  "Failed to open ofstream " << os_str);
+      std::ofstream os(vtu_prefix + GenerateVTUFileName(field_name, myid));
       qfield.second->SaveVTU(os, pv_data_format, GetCompressionLevel(), field_name);
    }
 
@@ -1008,10 +1000,7 @@ void ParaViewDataCollection::Save()
    {
       // Create the main PVTU file
       {
-         std::string os_str = vtu_prefix + GeneratePVTUFileName("data");
-         std::ofstream pvtu_out(os_str);
-         MFEM_VERIFY(pvtu_out.is_open(),
-                     "Failed to open ofstream " << os_str);
+         std::ofstream pvtu_out(vtu_prefix + GeneratePVTUFileName("data"));
          WritePVTUHeader(pvtu_out);
 
          // Grid function fields and coefficient fields
@@ -1066,10 +1055,8 @@ void ParaViewDataCollection::Save()
          const std::string &q_field_name = q_field.first;
          std::string q_fname = GeneratePVTUPath() + "/"
                                + GeneratePVTUFileName(q_field_name);
-         std::string os_str = col_path + "/" + q_fname;
-         std::ofstream pvtu_out(os_str);
-         MFEM_VERIFY(pvtu_out.is_open(),
-                     "Failed to open ofstream " << os_str);
+
+         std::ofstream pvtu_out(col_path + "/" + q_fname);
          WritePVTUHeader(pvtu_out);
          int vec_dim = q_field.second->GetVDim();
          pvtu_out << "<PPointData>\n";

fem/estimators.hpp

@@ -320,8 +320,8 @@ public:
     error estimation procedure where the flux averaging is replaced by a global
     L2 projection (requiring a mass matrix solve).
 
-    The required BilinearFormIntegrator must implement the method
-    ComputeElementFlux().
+    The required BilinearFormIntegrator must implement the methods
+    ComputeElementFlux() and ComputeFluxEnergy().
 
     Implemented for the parallel case only.
  */
@@ -357,8 +357,8 @@ protected:
 
 public:
    /** @brief Construct a new L2ZienkiewiczZhuEstimator object.
-       @param integ    This BilinearFormIntegrator must implement the method
-                       ComputeElementFlux().
+       @param integ    This BilinearFormIntegrator must implement the methods
+                       ComputeElementFlux() and ComputeFluxEnergy().
        @param sol      The solution field whose error is to be estimated.
        @param flux_fes The L2ZienkiewiczZhuEstimator assumes ownership of this
                        FiniteElementSpace and will call its Update() method when
@@ -382,8 +382,8 @@ public:
    { }
 
    /** @brief Construct a new L2ZienkiewiczZhuEstimator object.
-       @param integ    This BilinearFormIntegrator must implement the method
-                       ComputeElementFlux().
+       @param integ    This BilinearFormIntegrator must implement the methods
+                       ComputeElementFlux() and ComputeFluxEnergy().
        @param sol      The solution field whose error is to be estimated.
        @param flux_fes The L2ZienkiewiczZhuEstimator does NOT assume ownership
                        of this FiniteElementSpace; will call its Update() method

fem/fe.hpp

@@ -20,6 +20,7 @@
 #include "fe/fe_base.hpp"
 #include "fe/fe_fixed_order.hpp"
 #include "fe/fe_h1.hpp"
+#include "fe/fe_h1_bubble.hpp"
 #include "fe/fe_nd.hpp"
 #include "fe/fe_rt.hpp"
 #include "fe/fe_l2.hpp"

fem/fe/fe_base.hpp

@@ -349,7 +349,7 @@ public:
        vector-valued finite elements, which is also the width of the
        DenseMatrix argument in
        CalcPhysVShape(ElementTransformation &Trans, DenseMatrix &shape). */
-   virtual int GetPhysRangeDim(int /* space_dim */) const { return vdim; }
+   int GetPhysRangeDim(int /* space_dim */) const { return vdim; }
 
    /** Returns the dimension of the curl for vector-valued finite elements,
        which is also the width of the DenseMatrix argument in
@@ -360,7 +360,7 @@ public:
        finite elements, which is also the width of the DenseMatrix argument in
        CalcPhysCurlShape(ElementTransformation &Trans, DenseMatrix &curl_shape).
    */
-   virtual int GetPhysCurlDim(int /* space_dim */) const { return cdim; }
+   int GetPhysCurlDim(int /* space_dim */) const { return cdim; }
 
    /// Returns the Geometry::Type of the reference element.
    Geometry::Type GetGeomType() const { return geom_type; }
@@ -1017,7 +1017,7 @@ public:
    VectorFiniteElement(int D, Geometry::Type G, int Do, int O, int M,
                        int F = FunctionSpace::Pk);
 
-   int GetPhysRangeDim(int space_dim) const override { return space_dim; }
+   int GetPhysRangeDim(int space_dim) const { return space_dim; }
 };
 
 /// @brief Class for computing 1D special polynomials and their associated basis

fem/fe/fe_h1_bubble.cpp

@@ -0,0 +1,973 @@
+// Copyright (c) 2010-2025, 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.
+
+// H1 Finite Element classes
+
+#include "fe_h1_bubble.hpp"
+
+namespace mfem
+{
+
+using namespace std;
+
+H1Bubble_TriangleElement::H1Bubble_TriangleElement(int p, int q, int btype)
+   : NodalFiniteElement(2, Geometry::TRIANGLE, 3*p + ((q+1)*(q+2))/2,
+                        max(p, 3 + q), FunctionSpace::Pk),
+     base_order(p), bubble_order(q)
+{
+   const real_t *cp = poly1d.ClosedPoints(p, VerifyNodal(VerifyClosed(btype)));
+   const real_t *cp2 = poly1d.ClosedPoints(
+                          q + 3, VerifyNodal(VerifyClosed(btype)));
+
+   const int n1d = max(p + 1, q + 1);
+   const int npq = ((p+1)*(p+2))/2 + ((q+1)*(q+2))/2;
+
+#ifndef MFEM_THREAD_SAFE
+   shape_x.SetSize(n1d);
+   shape_y.SetSize(n1d);
+   shape_l.SetSize(n1d);
+   dshape_x.SetSize(n1d);
+   dshape_y.SetSize(n1d);
+   dshape_l.SetSize(n1d);
+   u.SetSize(npq);
+   du.SetSize(npq, dim);
+#endif
+
+   // vertices
+   Nodes.IntPoint(0).Set2(cp[0], cp[0]);
+   Nodes.IntPoint(1).Set2(cp[p], cp[0]);
+   Nodes.IntPoint(2).Set2(cp[0], cp[p]);
+
+   // edges
+   int o = 3;
+   for (int i = 1; i < p; i++)
+   {
+      Nodes.IntPoint(o++).Set2(cp[i], cp[0]);
+   }
+   for (int i = 1; i < p; i++)
+   {
+      Nodes.IntPoint(o++).Set2(cp[p-i], cp[i]);
+   }
+   for (int i = 1; i < p; i++)
+   {
+      Nodes.IntPoint(o++).Set2(cp[0], cp[p-i]);
+   }
+
+   // Interior P_{q+3} nodes
+   for (int j = 1; j < q + 3; j++)
+   {
+      for (int i = 1; i + j < q + 3; i++)
+      {
+         const real_t w = cp2[i] + cp2[j] + cp2[q+3-i-j];
+         Nodes.IntPoint(o++).Set2(cp2[i]/w, cp2[j]/w);
+      }
+   }
+
+#ifdef MFEM_THREAD_SAFE
+   Vector shape_x(n1d), shape_y(n1d), shape_l(n1d);
+#endif
+
+   DenseMatrix Tt(dof, npq);
+   for (int k = 0; k < dof; ++k)
+   {
+      const IntegrationPoint &ip = Nodes.IntPoint(k);
+      poly1d.CalcBasis(p, ip.x, shape_x);
+      poly1d.CalcBasis(p, ip.y, shape_y);
+      poly1d.CalcBasis(p, 1. - ip.x - ip.y, shape_l);
+      o = 0;
+      for (int j = 0; j <= p; j++)
+      {
+         for (int i = 0; i + j <= p; i++)
+         {
+            Tt(k, o++) = shape_x[i]*shape_y[j]*shape_l[p-i-j];
+         }
+      }
+
+      poly1d.CalcBasis(q, ip.x, shape_x);
+      poly1d.CalcBasis(q, ip.y, shape_y);
+      poly1d.CalcBasis(q, 1. - ip.x - ip.y, shape_l);
+      const real_t b_T = ip.x * ip.y * (1 - ip.x - ip.y);
+      for (int j = 0; j <= q; j++)
+      {
+         for (int i = 0; i + j <= q; i++)
+         {
+            Tt(k, o++) = b_T*shape_x[i]*shape_y[j]*shape_l[q-i-j];
+         }
+      }
+   }
+
+   // Compute left inverse of T (given Tt = T^T).
+   DenseMatrix TtT(dof, dof);
+   MultAAt(Tt, TtT);
+
+   DenseMatrixInverse TtT_inv(TtT);
+   T_pinv.SetSize(dof, dof);
+   TtT_inv.Mult(Tt, T_pinv);
+}
+
+void H1Bubble_TriangleElement::CalcShape(const IntegrationPoint &ip,
+                                         Vector &shape) const
+{
+   const int p = base_order;
+   const int q = bubble_order;
+
+#ifdef MFEM_THREAD_SAFE
+   const int n1d = max(p + 1, q + 1);
+   const int npq = ((p+1)*(p+2))/2 + ((q+1)*(q+2))/2;
+   Vector shape_x(n1d), shape_y(n1d), shape_l(n1d), u(npq);
+#endif
+
+   poly1d.CalcBasis(p, ip.x, shape_x);
+   poly1d.CalcBasis(p, ip.y, shape_y);
+   poly1d.CalcBasis(p, 1. - ip.x - ip.y, shape_l);
+
+   int o = 0;
+   for (int j = 0; j <= p; j++)
+   {
+      for (int i = 0; i + j <= p; i++)
+      {
+         u(o++) = shape_x[i]*shape_y[j]*shape_l[p-i-j];
+      }
+   }
+
+   poly1d.CalcBasis(q, ip.x, shape_x);
+   poly1d.CalcBasis(q, ip.y, shape_y);
+   poly1d.CalcBasis(q, 1. - ip.x - ip.y, shape_l);
+   const real_t b_T = ip.x * ip.y * (1 - ip.x - ip.y);
+
+   for (int j = 0; j <= q; j++)
+   {
+      for (int i = 0; i + j <= q; i++)
+      {
+         u(o++) = b_T*shape_x[i]*shape_y[j]*shape_l[q-i-j];
+      }
+   }
+
+   T_pinv.Mult(u, shape);
+}
+
+void H1Bubble_TriangleElement::CalcDShape(const IntegrationPoint &ip,
+                                          DenseMatrix &dshape) const
+{
+   const int p = base_order;
+   const int q = bubble_order;
+
+#ifdef MFEM_THREAD_SAFE
+   const int n1d = max(p + 1, q + 1);
+   const int npq = ((p+1)*(p+2))/2 + ((q+1)*(q+2))/2;
+   Vector  shape_x(n1d),  shape_y(n1d),  shape_l(n1d);
+   Vector dshape_x(n1d), dshape_y(n1d), dshape_l(n1d);
+   DenseMatrix du(npq, dim);
+#endif
+
+   const real_t lambda = 1.0 - ip.x - ip.y;
+
+   poly1d.CalcBasis(p, ip.x, shape_x, dshape_x);
+   poly1d.CalcBasis(p, ip.y, shape_y, dshape_y);
+   poly1d.CalcBasis(p, lambda, shape_l, dshape_l);
+
+   int o = 0;
+   for (int j = 0; j <= p; j++)
+   {
+      for (int i = 0; i + j <= p; i++)
+      {
+         int k = p - i - j;
+         du(o,0) = (dshape_x[i]*shape_l[k] - shape_x[i]*dshape_l[k])*shape_y[j];
+         du(o,1) = (dshape_y[j]* shape_l[k] - shape_y[j]*dshape_l[k])*shape_x[i];
+         o++;
+      }
+   }
+
+   poly1d.CalcBasis(q, ip.x, shape_x, dshape_x);
+   poly1d.CalcBasis(q, ip.y, shape_y, dshape_y);
+   poly1d.CalcBasis(q, lambda, shape_l, dshape_l);
+   const real_t b_T = ip.x * ip.y * lambda;
+   const real_t dxb_T = ip.y * (lambda - ip.x);
+   const real_t dyb_T = ip.x * (lambda - ip.y);
+
+   for (int j = 0; j <= q; j++)
+   {
+      for (int i = 0; i + j <= q; i++)
+      {
+         int k = q - i - j;
+         du(o,0) = shape_y[j]*(dxb_T*shape_x[i]*shape_l[k]
+                               + b_T*dshape_x[i]*shape_l[k]
+                               - b_T*shape_x[i]*dshape_l[k]);
+         du(o,1) = shape_x[i]*(dyb_T*shape_y[j]*shape_l[k]
+                               + b_T*dshape_y[j]*shape_l[k]
+                               - b_T*shape_y[j]*dshape_l[k]);
+         o++;
+      }
+   }
+
+   Mult(T_pinv, du, dshape);
+}
+
+H1Bubble_QuadrilateralElement::H1Bubble_QuadrilateralElement(
+   int p, int q, int btype)
+   : NodalFiniteElement(2, Geometry::SQUARE, 4*p + (q+1)*(q+1),
+                        max(p, 2 + q), FunctionSpace::Qk),
+     base_order(p), bubble_order(q)
+{
+   const real_t *cp = poly1d.ClosedPoints(p, VerifyNodal(VerifyClosed(btype)));
+   const real_t *cp2 = poly1d.ClosedPoints(
+                          q + 2, VerifyNodal(VerifyClosed(btype)));
+
+   const int n1d = max(p + 1, q + 1);
+   const int npq = (p+1)*(p+1) + (q+1)*(q+1);
+
+#ifndef MFEM_THREAD_SAFE
+   shape_x.SetSize(n1d);
+   shape_y.SetSize(n1d);
+   dshape_x.SetSize(n1d);
+   dshape_y.SetSize(n1d);
+
+   u.SetSize(npq);
+   du.SetSize(npq, dim);
+#endif
+
+   // vertices
+   Nodes.IntPoint(0).Set2(cp[0], cp[0]);
+   Nodes.IntPoint(1).Set2(cp[p], cp[0]);
+   Nodes.IntPoint(2).Set2(cp[p], cp[p]);
+   Nodes.IntPoint(3).Set2(cp[0], cp[p]);
+
+   // edges
+   int o = 4;
+   for (int i = 1; i < p; i++)
+   {
+      Nodes.IntPoint(o++).Set2(cp[i], cp[0]);
+   }
+   for (int i = 1; i < p; i++)
+   {
+      Nodes.IntPoint(o++).Set2(cp[p], cp[i]);
+   }
+   for (int i = 1; i < p; i++)
+   {
+      Nodes.IntPoint(o++).Set2(cp[p-i], cp[p]);
+   }
+   for (int i = 1; i < p; i++)
+   {
+      Nodes.IntPoint(o++).Set2(cp[0], cp[p-i]);
+   }
+
+   // interior P_{q+2} nodes
+   for (int j = 1; j < q+2; j++)
+   {
+      for (int i = 1; i < q+2; i++)
+      {
+         Nodes.IntPoint(o++).Set2(cp2[i], cp2[j]);
+      }
+   }
+
+#ifdef MFEM_THREAD_SAFE
+   Vector shape_x(n1d), shape_y(n1d);
+#endif
+
+   DenseMatrix Tt(dof, npq);
+   for (int k = 0; k < dof; ++k)
+   {
+      const IntegrationPoint &ip = Nodes.IntPoint(k);
+      poly1d.CalcBasis(p, ip.x, shape_x);
+      poly1d.CalcBasis(p, ip.y, shape_y);
+
+      o = 0;
+      for (int j = 0; j <= p; j++)
+      {
+         for (int i = 0; i <= p; i++)
+         {
+            Tt(k, o++) = shape_x[i]*shape_y[j];
+         }
+      }
+
+      poly1d.CalcBasis(q, ip.x, shape_x);
+      poly1d.CalcBasis(q, ip.y, shape_y);
+      const real_t b_T = ip.x*(1.0 - ip.x)*ip.y*(1.0 - ip.y);
+      for (int j = 0; j <= q; j++)
+      {
+         for (int i = 0; i <= q; i++)
+         {
+            Tt(k, o++) = b_T*shape_x[i]*shape_y[j];
+         }
+      }
+   }
+
+   // Compute left inverse of T (given Tt = T^T).
+   DenseMatrix TtT(dof, dof);
+   MultAAt(Tt, TtT);
+
+   DenseMatrixInverse TtT_inv(TtT);
+   T_pinv.SetSize(dof, dof);
+   TtT_inv.Mult(Tt, T_pinv);
+}
+
+void H1Bubble_QuadrilateralElement::CalcShape(const IntegrationPoint &ip,
+                                              Vector &shape) const
+{
+   const int p = base_order;
+   const int q = bubble_order;
+
+#ifdef MFEM_THREAD_SAFE
+   const int n1d = max(p + 1, q + 1);
+   const int npq = (p+1)*(p+1) + (q+1)*(q+1);
+   Vector shape_x(n1d), shape_y(n1d), u(npq);
+#endif
+
+   poly1d.CalcBasis(p, ip.x, shape_x);
+   poly1d.CalcBasis(p, ip.y, shape_y);
+
+   int o = 0;
+   for (int j = 0; j <= p; j++)
+   {
+      for (int i = 0; i <= p; i++)
+      {
+         u(o++) = shape_x[i]*shape_y[j];
+      }
+   }
+
+   poly1d.CalcBasis(q, ip.x, shape_x);
+   poly1d.CalcBasis(q, ip.y, shape_y);
+   const real_t b_T = ip.x*(1.0 - ip.x)*ip.y*(1.0 - ip.y);
+
+   for (int j = 0; j <= q; j++)
+   {
+      for (int i = 0; i <= q; i++)
+      {
+         u(o++) = b_T*shape_x[i]*shape_y[j];
+      }
+   }
+
+   T_pinv.Mult(u, shape);
+}
+
+void H1Bubble_QuadrilateralElement::CalcDShape(const IntegrationPoint &ip,
+                                               DenseMatrix &dshape) const
+{
+   const int p = base_order;
+   const int q = bubble_order;
+
+#ifdef MFEM_THREAD_SAFE
+   const int n1d = max(p + 1, q + 1);
+   const int npq = (p+1)*(p+1) + (q+1)*(q+1);
+   Vector shape_x(n1d), shape_y(n1d), dshape_x(n1d), dshape_y(n1d);
+   DenseMatrix du(npq, dim);
+#endif
+
+   poly1d.CalcBasis(p, ip.x, shape_x, dshape_x);
+   poly1d.CalcBasis(p, ip.y, shape_y, dshape_y);
+
+   int o = 0;
+   for (int j = 0; j <= p; j++)
+   {
+      for (int i = 0; i <= p; i++)
+      {
+         du(o,0) = dshape_x[i]*shape_y[j];
+         du(o,1) = shape_x[i]*dshape_y[j];
+         o += 1;
+      }
+   }
+
+   poly1d.CalcBasis(q, ip.x, shape_x, dshape_x);
+   poly1d.CalcBasis(q, ip.y, shape_y, dshape_y);
+   const real_t b_T = ip.x*(1.0 - ip.x)*ip.y*(1.0 - ip.y);
+   const real_t dxb_T = (1.0 - 2*ip.x)*ip.y*(1.0 - ip.y);
+   const real_t dyb_T = ip.x*(1.0 - ip.x)*(1.0 - 2*ip.y);
+
+   for (int j = 0; j <= q; j++)
+   {
+      for (int i = 0; i <= q; i++)
+      {
+         du(o,0) = (dxb_T*shape_x[i] + b_T*dshape_x[i])*shape_y[j];
+         du(o,1) = (dyb_T*shape_y[j] + b_T*dshape_y[j])*shape_x[i];
+         o += 1;
+      }
+   }
+
+   Mult(T_pinv, du, dshape);
+}
+
+H1Bubble_TetrahedronElement::H1Bubble_TetrahedronElement(
+   int p, int q, int btype)
+   : NodalFiniteElement(3, Geometry::TETRAHEDRON,
+                        2*(p*p + 1) + ((q+1)*(q+2)*(q+3))/6,
+                        max(p, 4 + q), FunctionSpace::Pk),
+     base_order(p), bubble_order(q)
+{
+   const real_t *cp = poly1d.ClosedPoints(p, VerifyNodal(VerifyClosed(btype)));
+   const real_t *cp2 = poly1d.ClosedPoints(
+                          q + 4, VerifyNodal(VerifyClosed(btype)));
+
+   const int n1d = max(p+1, q+1);
+   const int npq = ((p+1)*(p+2)*(p+3))/6 + ((q+1)*(q+2)*(q+3))/6;
+
+#ifndef MFEM_THREAD_SAFE
+   shape_x.SetSize(n1d);
+   shape_y.SetSize(n1d);
+   shape_z.SetSize(n1d);
+   shape_l.SetSize(n1d);
+   dshape_x.SetSize(n1d);
+   dshape_y.SetSize(n1d);
+   dshape_z.SetSize(n1d);
+   dshape_l.SetSize(n1d);
+   u.SetSize(npq);
+   du.SetSize(npq, dim);
+#else
+   Vector shape_x(n1d), shape_y(n1d), shape_z(n1d), shape_l(n1d);
+#endif
+
+   // vertices
+   Nodes.IntPoint(0).Set3(cp[0], cp[0], cp[0]);
+   Nodes.IntPoint(1).Set3(cp[p], cp[0], cp[0]);
+   Nodes.IntPoint(2).Set3(cp[0], cp[p], cp[0]);
+   Nodes.IntPoint(3).Set3(cp[0], cp[0], cp[p]);
+
+   // edges (see Tetrahedron::edges in mesh/tetrahedron.cpp)
+   int o = 4;
+   for (int i = 1; i < p; i++)  // (0,1)
+   {
+      Nodes.IntPoint(o++).Set3(cp[i], cp[0], cp[0]);
+   }
+   for (int i = 1; i < p; i++)  // (0,2)
+   {
+      Nodes.IntPoint(o++).Set3(cp[0], cp[i], cp[0]);
+   }
+   for (int i = 1; i < p; i++)  // (0,3)
+   {
+      Nodes.IntPoint(o++).Set3(cp[0], cp[0], cp[i]);
+   }
+   for (int i = 1; i < p; i++)  // (1,2)
+   {
+      Nodes.IntPoint(o++).Set3(cp[p-i], cp[i], cp[0]);
+   }
+   for (int i = 1; i < p; i++)  // (1,3)
+   {
+      Nodes.IntPoint(o++).Set3(cp[p-i], cp[0], cp[i]);
+   }
+   for (int i = 1; i < p; i++)  // (2,3)
+   {
+      Nodes.IntPoint(o++).Set3(cp[0], cp[p-i], cp[i]);
+   }
+
+   // faces (see Mesh::GenerateFaces in mesh/mesh.cpp)
+   for (int j = 1; j < p; j++)
+   {
+      for (int i = 1; i + j < p; i++)  // (1,2,3)
+      {
+         real_t w = cp[i] + cp[j] + cp[p-i-j];
+         Nodes.IntPoint(o++).Set3(cp[p-i-j]/w, cp[i]/w, cp[j]/w);
+      }
+   }
+   for (int j = 1; j < p; j++)
+   {
+      for (int i = 1; i + j < p; i++)  // (0,3,2)
+      {
+         real_t w = cp[i] + cp[j] + cp[p-i-j];
+         Nodes.IntPoint(o++).Set3(cp[0], cp[j]/w, cp[i]/w);
+      }
+   }
+   for (int j = 1; j < p; j++)
+   {
+      for (int i = 1; i + j < p; i++)  // (0,1,3)
+      {
+         real_t w = cp[i] + cp[j] + cp[p-i-j];
+         Nodes.IntPoint(o++).Set3(cp[i]/w, cp[0], cp[j]/w);
+      }
+   }
+   for (int j = 1; j < p; j++)
+   {
+      for (int i = 1; i + j < p; i++)  // (0,2,1)
+      {
+         real_t w = cp[i] + cp[j] + cp[p-i-j];
+         Nodes.IntPoint(o++).Set3(cp[j]/w, cp[i]/w, cp[0]);
+      }
+   }
+
+   // Interior P_{q+4} nodes
+   for (int k = 1; k < q + 4; k++)
+   {
+      for (int j = 1; j + k < q + 4; j++)
+      {
+         for (int i = 1; i + j + k < q + 4; i++)
+         {
+            real_t w = cp2[i] + cp2[j] + cp2[k] + cp2[q+4-i-j-k];
+            Nodes.IntPoint(o++).Set3(cp2[i]/w, cp2[j]/w, cp2[k]/w);
+         }
+      }
+   }
+
+   DenseMatrix Tt(dof, npq);
+   for (int m = 0; m < dof; ++m)
+   {
+      const IntegrationPoint &ip = Nodes.IntPoint(m);
+      poly1d.CalcBasis(p, ip.x, shape_x);
+      poly1d.CalcBasis(p, ip.y, shape_y);
+      poly1d.CalcBasis(p, ip.z, shape_z);
+      poly1d.CalcBasis(p, 1. - ip.x - ip.y - ip.z, shape_l);
+
+      o = 0;
+      for (int k = 0; k <= p; k++)
+      {
+         for (int j = 0; j + k <= p; j++)
+         {
+            for (int i = 0; i + j + k <= p; i++)
+            {
+               Tt(m, o++) = shape_x[i]*shape_y[j]*shape_z[k]*shape_l[p-i-j-k];
+            }
+         }
+      }
+
+      poly1d.CalcBasis(q, ip.x, shape_x);
+      poly1d.CalcBasis(q, ip.y, shape_y);
+      poly1d.CalcBasis(q, ip.z, shape_z);
+      poly1d.CalcBasis(q, 1. - ip.x - ip.y - ip.z, shape_l);
+      const real_t b_T = ip.x * ip.y * ip.z * (1 - ip.x - ip.y - ip.z);
+
+      for (int k = 0; k <= q; k++)
+      {
+         for (int j = 0; j + k <= q; j++)
+         {
+            for (int i = 0; i + j + k <= q; i++)
+            {
+               Tt(m, o++) = b_T*shape_x[i]*shape_y[j]*shape_z[k]*shape_l[q-i-j-k];
+            }
+         }
+      }
+   }
+
+   // Compute left inverse of T (given Tt = T^T).
+   DenseMatrix TtT(dof, dof);
+   MultAAt(Tt, TtT);
+
+   DenseMatrixInverse TtT_inv(TtT);
+   T_pinv.SetSize(dof, dof);
+   TtT_inv.Mult(Tt, T_pinv);
+}
+
+void H1Bubble_TetrahedronElement::CalcShape(const IntegrationPoint &ip,
+                                            Vector &shape) const
+{
+   const int p = base_order;
+   const int q = bubble_order;
+
+#ifdef MFEM_THREAD_SAFE
+   const int n1d = max(p + 1, q + 1);
+   const int npq = ((p+1)*(p+2)*(p+3))/6 + ((q+1)*(q+2)*(q+3))/6;
+   Vector shape_x(n1d), shape_y(n1d), shape_z(n1d), shape_l(n1d), u(npq);
+#endif
+
+   poly1d.CalcBasis(p, ip.x, shape_x);
+   poly1d.CalcBasis(p, ip.y, shape_y);
+   poly1d.CalcBasis(p, ip.z, shape_z);
+   poly1d.CalcBasis(p, 1. - ip.x - ip.y - ip.z, shape_l);
+
+   int o = 0;
+   for (int k = 0; k <= p; k++)
+   {
+      for (int j = 0; j + k <= p; j++)
+      {
+         for (int i = 0; i + j + k <= p; i++)
+         {
+            u[o++] = shape_x[i]*shape_y[j]*shape_z[k]*shape_l[p-i-j-k];
+         }
+      }
+   }
+
+   poly1d.CalcBasis(q, ip.x, shape_x);
+   poly1d.CalcBasis(q, ip.y, shape_y);
+   poly1d.CalcBasis(q, ip.z, shape_z);
+   poly1d.CalcBasis(q, 1. - ip.x - ip.y - ip.z, shape_l);
+   const real_t b_T = ip.x * ip.y * ip.z * (1 - ip.x - ip.y - ip.z);
+
+   for (int k = 0; k <= q; k++)
+   {
+      for (int j = 0; j + k <= q; j++)
+      {
+         for (int i = 0; i + j + k <= q; i++)
+         {
+            u(o++) = b_T*shape_x[i]*shape_y[j]*shape_z[k]*shape_l[q-i-j-k];
+         }
+      }
+   }
+
+   T_pinv.Mult(u, shape);
+}
+
+void H1Bubble_TetrahedronElement::CalcDShape(const IntegrationPoint &ip,
+                                             DenseMatrix &dshape) const
+{
+   const int p = base_order;
+   const int q = bubble_order;
+
+#ifdef MFEM_THREAD_SAFE
+   const int n1d = max(p+1, q+1);
+   const int npq = ((p+1)*(p+2)*(p+3))/6 + ((q+1)*(q+2)*(q+3))/6;
+
+   Vector shape_x(n1d), shape_y(n1d), shape_z(n1d), shape_l(n1d);
+   Vector dshape_x(n1d), dshape_y(n1d), dshape_z(n1d), dshape_l(n1d);
+   DenseMatrix du(npq, dim);
+#endif
+
+   const real_t lambda = 1.0 - ip.x - ip.y - ip.z;
+
+   poly1d.CalcBasis(p, ip.x, shape_x, dshape_x);
+   poly1d.CalcBasis(p, ip.y, shape_y, dshape_y);
+   poly1d.CalcBasis(p, ip.z, shape_z, dshape_z);
+   poly1d.CalcBasis(p, lambda, shape_l, dshape_l);
+
+   int o = 0;
+   for (int k = 0; k <= p; k++)
+   {
+      for (int j = 0; j + k <= p; j++)
+      {
+         for (int i = 0; i + j + k <= p; i++)
+         {
+            int l = p - i - j - k;
+            du(o,0) = (dshape_x[i]*shape_l[l] - shape_x[i]*dshape_l[l])
+                      *shape_y[j]*shape_z[k];
+            du(o,1) = (dshape_y[j]*shape_l[l] - shape_y[j]*dshape_l[l])
+                      *shape_x[i]*shape_z[k];
+            du(o,2) = (dshape_z[k]*shape_l[l] - shape_z[k]*dshape_l[l])
+                      *shape_x[i]*shape_y[j];
+            o++;
+         }
+      }
+   }
+
+   poly1d.CalcBasis(q, ip.x, shape_x, dshape_x);
+   poly1d.CalcBasis(q, ip.y, shape_y, dshape_y);
+   poly1d.CalcBasis(q, ip.z, shape_z, dshape_z);
+   poly1d.CalcBasis(q, lambda, shape_l, dshape_l);
+   const real_t b_T = ip.x * ip.y * ip.z * (1 - ip.x - ip.y - ip.z);
+   const real_t dxb_T = ip.y * ip.z * (lambda - ip.x);
+   const real_t dyb_T = ip.x * ip.z * (lambda - ip.y);
+   const real_t dzb_T = ip.x * ip.y * (lambda - ip.z);
+
+   for (int k = 0; k <= q; k++)
+   {
+      for (int j = 0; j + k <= q; j++)
+      {
+         for (int i = 0; i + j + k <= q; i++)
+         {
+            int l = q - i - j - k;
+            du(o,0) = shape_y[j]*shape_z[k]*(dxb_T*shape_x[i]*shape_l[l]
+                                             + b_T*dshape_x[i]*shape_l[l]
+                                             - b_T*shape_x[i]*dshape_l[l]);
+            du(o,1) = shape_x[i]*shape_z[k]*(dyb_T*shape_y[j]*shape_l[l]
+                                             + b_T*dshape_y[j]*shape_l[l]
+                                             - b_T*shape_y[j]*dshape_l[l]);
+            du(o,2) = shape_x[i]*shape_y[j]*(dzb_T*shape_z[k]*shape_l[l]
+                                             + b_T*dshape_z[k]*shape_l[l]
+                                             - b_T*shape_z[k]*dshape_l[l]);
+            o++;
+         }
+      }
+   }
+
+   Mult(T_pinv, du, dshape);
+}
+
+H1Bubble_HexahedronElement::H1Bubble_HexahedronElement(
+   int p, int q, int btype)
+   : NodalFiniteElement(3, Geometry::CUBE, (2 + 6*p*p) + (q+1)*(q+1)*(q+1),
+                        max(p, 2 + q), FunctionSpace::Qk),
+     base_order(p), bubble_order(q)
+{
+   const real_t *cp = poly1d.ClosedPoints(p, VerifyNodal(VerifyClosed(btype)));
+   const real_t *cp2 = poly1d.ClosedPoints(
+                          q + 2, VerifyNodal(VerifyClosed(btype)));
+
+   const int n1d = max(p + 1, q + 1);
+   const int npq = (p+1)*(p+1)*(p+1) + (q+1)*(q+1)*(q+1);
+
+#ifndef MFEM_THREAD_SAFE
+   shape_x.SetSize(n1d);
+   shape_y.SetSize(n1d);
+   shape_z.SetSize(n1d);
+   dshape_x.SetSize(n1d);
+   dshape_y.SetSize(n1d);
+   dshape_z.SetSize(n1d);
+
+   u.SetSize(npq);
+   du.SetSize(npq, dim);
+#endif
+
+   // vertices
+   Nodes.IntPoint(0).Set3(cp[0], cp[0], cp[0]);
+   Nodes.IntPoint(1).Set3(cp[p], cp[0], cp[0]);
+   Nodes.IntPoint(2).Set3(cp[p], cp[p], cp[0]);
+   Nodes.IntPoint(3).Set3(cp[0], cp[p], cp[0]);
+
+   Nodes.IntPoint(4).Set3(cp[0], cp[0], cp[p]);
+   Nodes.IntPoint(5).Set3(cp[p], cp[0], cp[p]);
+   Nodes.IntPoint(6).Set3(cp[p], cp[p], cp[p]);
+   Nodes.IntPoint(7).Set3(cp[0], cp[p], cp[p]);
+
+   int o = 8;
+   for (int i = 1; i < p; i++)
+   {
+      Nodes.IntPoint(o++).Set3(cp[i], cp[0], cp[0]);   // (0,1)
+   }
+   for (int i = 1; i < p; i++)
+   {
+      Nodes.IntPoint(o++).Set3(cp[p], cp[i], cp[0]);   // (1,2)
+   }
+   for (int i = 1; i < p; i++)
+   {
+      Nodes.IntPoint(o++).Set3(cp[i], cp[p], cp[0]);   // (3,2)
+   }
+   for (int i = 1; i < p; i++)
+   {
+      Nodes.IntPoint(o++).Set3(cp[0], cp[i], cp[0]);   // (0,3)
+   }
+   for (int i = 1; i < p; i++)
+   {
+      Nodes.IntPoint(o++).Set3(cp[i], cp[0], cp[p]);   // (4,5)
+   }
+   for (int i = 1; i < p; i++)
+   {
+      Nodes.IntPoint(o++).Set3(cp[p], cp[i], cp[p]);   // (5,6)
+   }
+   for (int i = 1; i < p; i++)
+   {
+      Nodes.IntPoint(o++).Set3(cp[i], cp[p], cp[p]);   // (7,6)
+   }
+   for (int i = 1; i < p; i++)
+   {
+      Nodes.IntPoint(o++).Set3(cp[0], cp[i], cp[p]);   // (4,7)
+   }
+   for (int i = 1; i < p; i++)
+   {
+      Nodes.IntPoint(o++).Set3(cp[0], cp[0], cp[i]);   // (0,4)
+   }
+   for (int i = 1; i < p; i++)
+   {
+      Nodes.IntPoint(o++).Set3(cp[p], cp[0], cp[i]);   // (1,5)
+   }
+   for (int i = 1; i < p; i++)
+   {
+      Nodes.IntPoint(o++).Set3(cp[p], cp[p], cp[i]);   // (2,6)
+   }
+   for (int i = 1; i < p; i++)
+   {
+      Nodes.IntPoint(o++).Set3(cp[0], cp[p], cp[i]);   // (3,7)
+   }
+
+   // faces
+   for (int j = 1; j < p; j++)
+   {
+      for (int i = 1; i < p; i++)
+      {
+         Nodes.IntPoint(o++).Set3(cp[i], cp[p-j], cp[0]);   // (3,2,1,0)
+      }
+   }
+   for (int j = 1; j < p; j++)
+   {
+      for (int i = 1; i < p; i++)
+      {
+         Nodes.IntPoint(o++).Set3(cp[i], cp[0], cp[j]);   // (0,1,5,4)
+      }
+   }
+   for (int j = 1; j < p; j++)
+   {
+      for (int i = 1; i < p; i++)
+      {
+         Nodes.IntPoint(o++).Set3(cp[p], cp[i], cp[j]);   // (1,2,6,5)
+      }
+   }
+   for (int j = 1; j < p; j++)
+   {
+      for (int i = 1; i < p; i++)
+      {
+         Nodes.IntPoint(o++).Set3(cp[p-i], cp[p], cp[j]);   // (2,3,7,6)
+      }
+   }
+   for (int j = 1; j < p; j++)
+   {
+      for (int i = 1; i < p; i++)
+      {
+         Nodes.IntPoint(o++).Set3(cp[0], cp[p-i], cp[j]);   // (3,0,4,7)
+      }
+   }
+   for (int j = 1; j < p; j++)
+   {
+      for (int i = 1; i < p; i++)
+      {
+         Nodes.IntPoint(o++).Set3(cp[i], cp[j], cp[p]);   // (4,5,6,7)
+      }
+   }
+
+   // interior P_{q+2} nodes
+   for (int k = 1; k < q+2; k++)
+   {
+      for (int j = 1; j < q+2; j++)
+      {
+         for (int i = 1; i < q+2; i++)
+         {
+            Nodes.IntPoint(o++).Set3(cp2[i], cp2[j], cp2[k]);
+         }
+      }
+   }
+
+#ifdef MFEM_THREAD_SAFE
+   Vector shape_x(n1d), shape_y(n1d), shape_z(n1d);
+#endif
+
+   DenseMatrix Tt(dof, npq);
+   for (int m = 0; m < dof; ++m)
+   {
+      const IntegrationPoint &ip = Nodes.IntPoint(m);
+      poly1d.CalcBasis(p, ip.x, shape_x);
+      poly1d.CalcBasis(p, ip.y, shape_y);
+      poly1d.CalcBasis(p, ip.z, shape_z);
+
+      o = 0;
+      for (int k = 0; k <= p; k++)
+      {
+         for (int j = 0; j <= p; j++)
+         {
+            for (int i = 0; i <= p; i++)
+            {
+               Tt(m, o++) = shape_x[i]*shape_y[j]*shape_z[k];
+            }
+         }
+      }
+
+      poly1d.CalcBasis(q, ip.x, shape_x);
+      poly1d.CalcBasis(q, ip.y, shape_y);
+      poly1d.CalcBasis(q, ip.z, shape_z);
+      const real_t b_T = ip.x*(1.0 - ip.x)*ip.y*(1.0 - ip.y)*ip.z*(1.0 - ip.z);
+      for (int k = 0; k <= q; k++)
+      {
+         for (int j = 0; j <= q; j++)
+         {
+            for (int i = 0; i <= q; i++)
+            {
+               Tt(m, o++) = b_T*shape_x[i]*shape_y[j]*shape_z[k];
+            }
+         }
+      }
+   }
+
+   // Compute left inverse of T (given Tt = T^T).
+   DenseMatrix TtT(dof, dof);
+   MultAAt(Tt, TtT);
+
+   DenseMatrixInverse TtT_inv(TtT);
+   T_pinv.SetSize(dof, dof);
+   TtT_inv.Mult(Tt, T_pinv);
+}
+
+void H1Bubble_HexahedronElement::CalcShape(const IntegrationPoint &ip,
+                                           Vector &shape) const
+{
+   const int p = base_order;
+   const int q = bubble_order;
+
+#ifdef MFEM_THREAD_SAFE
+   const int n1d = max(p + 1, q + 1);
+   const int npq = (p+1)*(p+1)*(p+1) + (q+1)*(q+1)*(q+1);
+   Vector shape_x(n1d), shape_y(n1d), shape_z(n1d), u(npq);
+#endif
+
+   poly1d.CalcBasis(p, ip.x, shape_x);
+   poly1d.CalcBasis(p, ip.y, shape_y);
+   poly1d.CalcBasis(p, ip.z, shape_z);
+
+   int o = 0;
+   for (int k = 0; k <= p; k++)
+   {
+      for (int j = 0; j <= p; j++)
+      {
+         for (int i = 0; i <= p; i++)
+         {
+            u(o++) = shape_x[i]*shape_y[j]*shape_z[k];
+         }
+      }
+   }
+
+   poly1d.CalcBasis(q, ip.x, shape_x);
+   poly1d.CalcBasis(q, ip.y, shape_y);
+   poly1d.CalcBasis(q, ip.z, shape_z);
+   const real_t b_T = ip.x*(1.0 - ip.x)*ip.y*(1.0 - ip.y)*ip.z*(1.0 - ip.z);
+
+   for (int k = 0; k <= q; k++)
+   {
+      for (int j = 0; j <= q; j++)
+      {
+         for (int i = 0; i <= q; i++)
+         {
+            u(o++) = b_T*shape_x[i]*shape_y[j]*shape_z[k];
+         }
+      }
+   }
+
+   T_pinv.Mult(u, shape);
+}
+
+void H1Bubble_HexahedronElement::CalcDShape(const IntegrationPoint &ip,
+                                            DenseMatrix &dshape) const
+{
+   const int p = base_order;
+   const int q = bubble_order;
+
+#ifdef MFEM_THREAD_SAFE
+   const int n1d = max(p + 1, q + 1);
+   const int npq = (p+1)*(p+1)*(p+1) + (q+1)*(q+1)*(q+1);
+   Vector shape_x(n1d), shape_y(n1d), shape_z(n1d), dshape_x(n1d),
+          dshape_y(n1d), dshape_z(n1d);
+   DenseMatrix du(npq, dim);
+#endif
+
+   poly1d.CalcBasis(p, ip.x, shape_x, dshape_x);
+   poly1d.CalcBasis(p, ip.y, shape_y, dshape_y);
+   poly1d.CalcBasis(p, ip.z, shape_z, dshape_z);
+
+   int o = 0;
+   for (int k = 0; k <= p; k++)
+   {
+      for (int j = 0; j <= p; j++)
+      {
+         for (int i = 0; i <= p; i++)
+         {
+            du(o,0) = dshape_x[i]*shape_y[j]*shape_z[k];
+            du(o,1) = shape_x[i]*dshape_y[j]*shape_z[k];
+            du(o,2) = shape_x[i]*shape_y[j]*dshape_z[k];
+            o += 1;
+         }
+      }
+   }
+
+   poly1d.CalcBasis(q, ip.x, shape_x, dshape_x);
+   poly1d.CalcBasis(q, ip.y, shape_y, dshape_y);
+   poly1d.CalcBasis(q, ip.z, shape_z, dshape_z);
+   const real_t b_T = ip.x*(1.0 - ip.x)*ip.y*(1.0 - ip.y)*ip.z*(1.0 - ip.z);
+   const real_t dxb_T = (1.0 - 2*ip.x)*ip.y*(1.0 - ip.y)*ip.z*(1.0 - ip.z);
+   const real_t dyb_T = ip.x*(1.0 - ip.x)*(1.0 - 2*ip.y)*ip.z*(1.0 - ip.z);
+   const real_t dzb_T = ip.x*(1.0 - ip.x)*ip.y*(1.0 - ip.y)*(1.0 - 2*ip.z);
+
+   for (int k = 0; k <= q; k++)
+   {
+      for (int j = 0; j <= q; j++)
+      {
+         for (int i = 0; i <= q; i++)
+         {
+            du(o,0) = (dxb_T*shape_x[i] + b_T*dshape_x[i])*shape_y[j]*shape_z[k];
+            du(o,1) = (dyb_T*shape_y[j] + b_T*dshape_y[j])*shape_x[i]*shape_z[k];
+            du(o,2) = (dzb_T*shape_z[k] + b_T*dshape_z[k])*shape_x[i]*shape_y[j];
+            o += 1;
+         }
+      }
+   }
+
+   Mult(T_pinv, du, dshape);
+}
+
+}

fem/fe/fe_h1_bubble.hpp

@@ -0,0 +1,109 @@
+// Copyright (c) 2010-2025, 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_FE_H1_BUBBLE
+#define MFEM_FE_H1_BUBBLE
+
+#include "fe_base.hpp"
+
+namespace mfem
+{
+
+/// Arbitrary order H1 plus bubble elements in 2D on a triangle
+class H1Bubble_TriangleElement : public NodalFiniteElement
+{
+private:
+#ifndef MFEM_THREAD_SAFE
+   mutable Vector shape_x, shape_y, shape_l, dshape_x, dshape_y, dshape_l, u;
+   mutable DenseMatrix du;
+#endif
+   int base_order;
+   int bubble_order;
+   DenseMatrix T_pinv;
+
+public:
+   /// @brief Construct the triangular bubble element with degree-p polynomials,
+   /// enriched with cubic bubble times degree q polynomial.
+   H1Bubble_TriangleElement(int p, int q, int btype = BasisType::GaussLobatto);
+   void CalcShape(const IntegrationPoint &ip, Vector &shape) const override;
+   void CalcDShape(const IntegrationPoint &ip,
+                   DenseMatrix &dshape) const override;
+};
+
+/// Arbitrary order H1 plus bubble elements in 2D on a quadrilateral
+class H1Bubble_QuadrilateralElement : public NodalFiniteElement
+{
+private:
+#ifndef MFEM_THREAD_SAFE
+   mutable Vector shape_x, shape_y, dshape_x, dshape_y, u;
+   mutable DenseMatrix du;
+#endif
+   int base_order;
+   int bubble_order;
+   DenseMatrix T_pinv;
+
+public:
+   /// @brief Construct the quadrilateral bubble element with degree-p
+   /// polynomials, enriched with biquadratic bubble times degree q polynomial.
+   H1Bubble_QuadrilateralElement(
+      int p, int q, int btype = BasisType::GaussLobatto);
+   void CalcShape(const IntegrationPoint &ip, Vector &shape) const override;
+   void CalcDShape(const IntegrationPoint &ip,
+                   DenseMatrix &dshape) const override;
+};
+
+/// Arbitrary order H1 plus bubble elements in 3D on a tetrahedron
+class H1Bubble_TetrahedronElement : public NodalFiniteElement
+{
+private:
+#ifndef MFEM_THREAD_SAFE
+   mutable Vector shape_x, shape_y, shape_z, shape_l;
+   mutable Vector dshape_x, dshape_y, dshape_z, dshape_l, u;
+   mutable DenseMatrix du;
+#endif
+   int base_order;
+   int bubble_order;
+   DenseMatrix T_pinv;
+
+public:
+   /// @brief Construct the tetrahedral bubble element with degree-p
+   /// polynomials, enriched with quartic bubble times degree q polynomial.
+   H1Bubble_TetrahedronElement(int p, int q, int btype = BasisType::GaussLobatto);
+   void CalcShape(const IntegrationPoint &ip, Vector &shape) const override;
+   void CalcDShape(const IntegrationPoint &ip,
+                   DenseMatrix &dshape) const override;
+};
+
+/// Arbitrary order H1 plus bubble elements in 3D on a hexahedron
+class H1Bubble_HexahedronElement : public NodalFiniteElement
+{
+private:
+#ifndef MFEM_THREAD_SAFE
+   mutable Vector shape_x, shape_y, shape_z;
+   mutable Vector dshape_x, dshape_y, dshape_z, u;
+   mutable DenseMatrix du;
+#endif
+   int base_order;
+   int bubble_order;
+   DenseMatrix T_pinv;
+
+public:
+   /// @brief Construct the hexahedral bubble element with degree-p polynomials,
+   /// enriched with triquadratic bubble times degree q polynomial.
+   H1Bubble_HexahedronElement(int p, int q, int btype = BasisType::GaussLobatto);
+   void CalcShape(const IntegrationPoint &ip, Vector &shape) const override;
+   void CalcDShape(const IntegrationPoint &ip,
+                   DenseMatrix &dshape) const override;
+};
+
+} // namespace mfem
+
+#endif

fem/fe/fe_nd.hpp

@@ -663,8 +663,8 @@ public:
                          const int cb_type = BasisType::GaussLobatto,
                          const int ob_type = BasisType::GaussLegendre);
 
-   int GetPhysRangeDim(int space_dim) const override { return 2; }
-   int GetPhysCurlDim(int space_dim) const override { return 1; }
+   int GetPhysRangeDim(int space_dim) const { return 2; }
+   int GetPhysCurlDim(int space_dim) const { return 1; }
 
    void CalcVShape(const IntegrationPoint &ip,
                    DenseMatrix &shape) const override;
@@ -708,8 +708,8 @@ private:
                            DenseMatrix &I) const;
 
 public:
-   int GetPhysRangeDim(int space_dim) const override { return 3; }
-   int GetPhysCurlDim(int space_dim) const override { return 3; }
+   int GetPhysRangeDim(int space_dim) const { return 3; }
+   int GetPhysCurlDim(int space_dim) const { return 3; }
 
    using FiniteElement::CalcVShape;
    using FiniteElement::CalcPhysCurlShape;

fem/fe/fe_rt.hpp

@@ -510,8 +510,8 @@ public:
    RT_R2D_SegmentElement(const int p,
                          const int ob_type = BasisType::GaussLegendre);
 
-   int GetPhysRangeDim(int space_dim) const override { return 2; }
-   int GetPhysCurlDim(int space_dim) const override { return 0; }
+   int GetPhysRangeDim(int space_dim) const { return 2; }
+   int GetPhysCurlDim(int space_dim) const { return 0; }
 
    void CalcVShape(const IntegrationPoint &ip,
                    DenseMatrix &shape) const override;
@@ -550,8 +550,8 @@ private:
                            DenseMatrix &I) const;
 
 public:
-   int GetPhysRangeDim(int space_dim) const override { return 3; }
-   int GetPhysCurlDim(int space_dim) const override { return 0; }
+   int GetPhysRangeDim(int space_dim) const { return 3; }
+   int GetPhysCurlDim(int space_dim) const { return 0; }
 
    using FiniteElement::CalcVShape;
 

fem/fe_coll.cpp

@@ -243,11 +243,21 @@ FiniteElementCollection *FiniteElementCollection::New(const char *name)
    {
       fec = new H1Ser_FECollection(atoi(name + 10), atoi(name + 6));
    }
+   else if (!strncmp(name, "H1Bubble_", 9))
+   {
+      fec = new H1Bubble_FECollection(atoi(name + 13), atoi(name + 16),
+                                      atoi(name + 9));
+   }
    else if (!strncmp(name, "H1@", 3))
    {
       fec = new H1_FECollection(atoi(name + 9), atoi(name + 5),
                                 BasisType::GetType(name[3]));
    }
+   else if (!strncmp(name, "H1Bubble@", 9))
+   {
+      fec = new H1Bubble_FECollection(atoi(name + 15), atoi(name + 18),
+                                      atoi(name + 11), BasisType::GetType(name[9]));
+   }
    else if (!strncmp(name, "L2_T", 4))
       fec = new L2_FECollection(atoi(name + 10), atoi(name + 6),
                                 atoi(name + 4));
@@ -2122,6 +2132,173 @@ H1_FECollection::~H1_FECollection()
    }
 }
 
+static int GetBubbleSpaceOrder(int p, int q, int dim)
+{
+   switch (dim)
+   {
+      case 0: return 0;
+      case 1: return std::max(p, q + 2);
+      case 2: return std::max(p, q + 3);
+      case 3: return std::max(p, q + 4);
+   }
+   MFEM_ABORT("Unsupported dimension.");
+}
+
+H1Bubble_FECollection::H1Bubble_FECollection(const int p, const int q,
+                                             const int dim, const int btype)
+   : FiniteElementCollection(GetBubbleSpaceOrder(p, q, dim)),
+     dim(dim),
+     b_type(BasisType::Check(btype)),
+     h1_order(p),
+     bubble_order(q)
+{
+   MFEM_VERIFY(p >= 1, "H1Bubble_FECollection requires order >= 1.");
+   MFEM_VERIFY(dim >= 0 && dim <= 3, "Unsupported dimension.");
+
+   switch (btype)
+   {
+      case BasisType::GaussLobatto:
+      {
+         snprintf(fec_name, 32, "H1Bubble_%dD_P%d_P%d", dim, p, q);
+         break;
+      }
+      default:
+      {
+         const int pt_type = BasisType::GetQuadrature1D(btype);
+         MFEM_VERIFY(Quadrature1D::CheckClosed(pt_type) != Quadrature1D::Invalid,
+                     "unsupported BasisType: " << BasisType::Name(btype));
+         snprintf(fec_name, 32, "H1Bubble@%c_%dD_P%d_P%d",
+                  (int)BasisType::GetChar(btype), dim, p, q);
+      }
+   }
+
+   dofs[Geometry::POINT] = 1;
+   elements[Geometry::POINT] = make_unique<PointFiniteElement>();
+
+   if (dim >= 1)
+   {
+      dofs[Geometry::SEGMENT] = p - 1;
+      elements[Geometry::SEGMENT] = make_unique<H1_SegmentElement>(p, btype);
+   }
+
+   if (dim == 2)
+   {
+      dofs[Geometry::TRIANGLE] = ((q+1)*(q+2))/2;
+      dofs[Geometry::SQUARE] = (q+1)*(q+1);
+
+      elements[Geometry::TRIANGLE] =
+         make_unique<H1Bubble_TriangleElement>(p, q, btype);
+      elements[Geometry::SQUARE] =
+         make_unique<H1Bubble_QuadrilateralElement>(p, q, btype);
+   }
+
+   if (dim == 3)
+   {
+      dofs[Geometry::TRIANGLE] = ((p-1)*(p-2))/2;
+      dofs[Geometry::SQUARE] = (p-1)*(p-1);
+      dofs[Geometry::TETRAHEDRON] = ((q+1)*(q+2)*(q+3))/6;
+      dofs[Geometry::CUBE] = (q+1)*(q+1)*(q+1);
+
+      elements[Geometry::TRIANGLE] = make_unique<H1_TriangleElement>(p, btype);
+      elements[Geometry::SQUARE] = make_unique<H1_QuadrilateralElement>(p, btype);
+
+      elements[Geometry::TETRAHEDRON] =
+         make_unique<H1Bubble_TetrahedronElement>(p, q, btype);
+      elements[Geometry::CUBE] =
+         make_unique<H1Bubble_HexahedronElement>(p, q, btype);
+   }
+
+   // DOF orderings. Need only for lower-dimensional entities.
+   // Segment DOF orderings in 2D.
+   if (dim >= 2)
+   {
+      seg_dof_ord[0].resize(p - 1);
+      seg_dof_ord[1].resize(p - 1);
+      for (int i = 0; i < p - 1; i++)
+      {
+         seg_dof_ord[0][i] = i;
+         seg_dof_ord[1][i] = p - 2 - i;
+      }
+   }
+
+   // Face (triangle or quadrilateral) DOF orderings in 3D.
+   if (dim == 3)
+   {
+      const int n_tri_dof = dofs[Geometry::TRIANGLE];
+      for (int i = 0; i < 6; i++)
+      {
+         tri_dof_ord[i].resize(n_tri_dof);
+      }
+      // see Mesh::GetTriOrientation in mesh/mesh.cpp
+      const int pm1 = p - 1;
+      const int pm2 = p - 2;
+      for (int j = 0; j < pm2; j++)
+      {
+         for (int i = 0; i + j < pm2; i++)
+         {
+            int o = n_tri_dof - ((pm1 - j)*(pm2 - j))/2 + i;
+            int k = (p - 3) - j - i;
+            tri_dof_ord[0][o] = o;  // (0,1,2)
+            tri_dof_ord[1][o] = n_tri_dof - ((pm1-j)*(pm2-j))/2 + k;  // (1,0,2)
+            tri_dof_ord[2][o] = n_tri_dof - ((pm1-i)*(pm2-i))/2 + k;  // (2,0,1)
+            tri_dof_ord[3][o] = n_tri_dof - ((pm1-k)*(pm2-k))/2 + i;  // (2,1,0)
+            tri_dof_ord[4][o] = n_tri_dof - ((pm1-k)*(pm2-k))/2 + j;  // (1,2,0)
+            tri_dof_ord[5][o] = n_tri_dof - ((pm1-i)*(pm2-i))/2 + j;  // (0,2,1)
+         }
+      }
+
+      const int n_quad_dof = dofs[Geometry::SQUARE];
+      for (int i = 0; i < 8; i++)
+      {
+         quad_dof_ord[i].resize(n_quad_dof);
+      }
+      for (int j = 0; j < pm1; j++)
+      {
+         for (int i = 0; i < pm1; i++)
+         {
+            int o = i + j*pm1;
+            quad_dof_ord[0][o] = i + j*pm1;  // (0,1,2,3)
+            quad_dof_ord[1][o] = j + i*pm1;  // (0,3,2,1)
+            quad_dof_ord[2][o] = j + (pm2 - i)*pm1;  // (1,2,3,0)
+            quad_dof_ord[3][o] = (pm2 - i) + j*pm1;  // (1,0,3,2)
+            quad_dof_ord[4][o] = (pm2 - i) + (pm2 - j)*pm1;  // (2,3,0,1)
+            quad_dof_ord[5][o] = (pm2 - j) + (pm2 - i)*pm1;  // (2,1,0,3)
+            quad_dof_ord[6][o] = (pm2 - j) + i*pm1;  // (3,0,1,2)
+            quad_dof_ord[7][o] = i + (pm2 - j)*pm1;  // (3,2,1,0)
+         }
+      }
+   }
+}
+
+const FiniteElement *
+H1Bubble_FECollection::FiniteElementForGeometry(Geometry::Type GeomType) const
+{
+   return elements[GeomType].get();
+}
+
+const int *H1Bubble_FECollection::DofOrderForOrientation(
+   Geometry::Type GeomType, int Or) const
+{
+   if (GeomType == Geometry::SEGMENT)
+   {
+      return (Or > 0) ? seg_dof_ord[0].data() : seg_dof_ord[1].data();
+   }
+   else if (GeomType == Geometry::TRIANGLE)
+   {
+      return tri_dof_ord[Or%6].data();
+   }
+   else if (GeomType == Geometry::SQUARE)
+   {
+      return quad_dof_ord[Or%8].data();
+   }
+   return nullptr;
+}
+
+FiniteElementCollection *H1Bubble_FECollection::GetTraceCollection() const
+{
+   return (dim < 0) ? NULL : new H1_Trace_FECollection(h1_order, dim, b_type);
+}
+
 
 H1_Trace_FECollection::H1_Trace_FECollection(const int p, const int dim,
                                              const int btype)

fem/fe_coll.hpp

@@ -111,6 +111,8 @@ public:
    | :------: | :---: | :---: | :-------: | :-----: | :---: |
    | H1_[DIM]_[ORDER] | H1 | * | 1 | VALUE | H1 nodal elements |
    | H1@[BTYPE]_[DIM]_[ORDER] | H1 | * | * | VALUE | H1 nodal elements |
+   | H1Bubble_[DIM]_[ORDER]_[BUBBLE_ORDER] | H1 | * | 1 | VALUE | H1 nodal elements enriched with bubble functions |
+   | H1Bubble@[BTYPE]_[DIM]_[ORDER]_[BUBBLE_ORDER] | H1 | * | 1 | VALUE | H1 nodal elements enriched with bubble functions |
    | H1Pos_[DIM]_[ORDER] | H1 | * | 2 | VALUE | H1 nodal elements |
    | H1Pos_Trace_[DIM]_[ORDER] | H^{1/2} | * | 2 | VALUE | H^{1/2}-conforming trace elements for H1 defined on the interface between mesh elements (faces,edges,vertices) |
    | H1_Trace_[DIM]_[ORDER] | H^{1/2} | * | 1 | VALUE | H^{1/2}-conforming trace elements for H1 defined on the interface between mesh elements (faces,edges,vertices) |
@@ -317,6 +319,59 @@ public:
    virtual ~H1_FECollection();
 };
 
+/// @brief Arbitrary order $H^1$-conforming (continuous) finite elements
+/// enriched with bubble functions.
+///
+/// The bubble space consists of the standard $P_p$ or $Q_p$ space, enriched
+/// with bubble functions, which are degree-$q$ polynomials times $b$, where $b$
+/// is the lowest-order bubble function.
+///
+/// The traces are the same as the standard $H^1$ traces.
+class H1Bubble_FECollection : public FiniteElementCollection
+{
+protected:
+   int dim;
+   int b_type;
+   int h1_order;
+   int bubble_order;
+
+   char fec_name[32];
+   std::array<int, Geometry::NumGeom> dofs{}; // zero initialize
+   std::array<std::unique_ptr<FiniteElement>, Geometry::NumGeom> elements;
+
+   std::array<std::vector<int>, 2> seg_dof_ord;
+   std::array<std::vector<int>, 6> tri_dof_ord;
+   std::array<std::vector<int>, 8> quad_dof_ord;
+   std::array<std::vector<int>, 24> tet_dof_ord;
+
+public:
+   /// Construct the $H^1$ bubble collection consisting of degree-$p$
+   /// polynomials enriched with the bubble function times degree-$q$
+   /// polynomials.
+   explicit H1Bubble_FECollection(const int p, const int q, const int dim = 3,
+                                  const int btype = BasisType::GaussLobatto);
+
+   const FiniteElement *
+   FiniteElementForGeometry(Geometry::Type GeomType) const override;
+
+   int DofForGeometry(Geometry::Type GeomType) const override
+   { return dofs[GeomType]; }
+
+   const int *DofOrderForOrientation(Geometry::Type GeomType,
+                                     int Or) const override;
+
+   const char *Name() const override { return fec_name; }
+
+   int GetContType() const override { return CONTINUOUS; }
+
+   int GetBasisType() const { return b_type; }
+
+   FiniteElementCollection *GetTraceCollection() const override;
+
+   FiniteElementCollection *Clone(int p) const override
+   { return new H1Bubble_FECollection(p, bubble_order, dim, b_type); }
+};
+
 /** @brief Arbitrary order H1-conforming (continuous) finite elements with
     positive basis functions. */
 class H1Pos_FECollection : public H1_FECollection

fem/fespace.cpp

@@ -3877,9 +3877,12 @@ const FiniteElement *FiniteElementSpace::GetFE(int i) const
    else
    {
 #ifdef MFEM_DEBUG
-      // consistency check: fec->GetOrder() and FE->GetOrder() should return
-      // the same value (for standard, constant-order spaces)
-      if (!IsVariableOrder() && FE->GetDim() > 0)
+      // Consistency check: fec->GetOrder() and FE->GetOrder() should return
+      // the same value (for standard, constant-order spaces). Skip this check
+      // even for constant-order bubble spaces, since the bubble functions on
+      // different geometries have different orders.
+      if (!IsVariableOrder() && FE->GetDim() > 0 &&
+          dynamic_cast<const H1Bubble_FECollection*>(fec) == nullptr)
       {
          MFEM_ASSERT(FE->GetOrder() == fec->GetOrder(),
                      "internal error: " <<
@@ -3934,16 +3937,6 @@ const FiniteElement *FiniteElementSpace::GetBE(int i) const
    return BE;
 }
 
-const FiniteElement *FiniteElementSpace::GetTypicalBE() const
-{
-   if (mesh->GetNBE() > 0) { return GetBE(0); }
-
-   Geometry::Type geom = mesh->GetTypicalFaceGeometry();
-   const FiniteElement *be = fec->FiniteElementForGeometry(geom);
-   MFEM_VERIFY(be != nullptr, "Could not determine a typical BE!");
-   return be;
-}
-
 const FiniteElement *FiniteElementSpace::GetFaceElement(int i) const
 {
    MFEM_VERIFY(!IsVariableOrder(), "not implemented");
@@ -3974,11 +3967,6 @@ const FiniteElement *FiniteElementSpace::GetFaceElement(int i) const
    return fe;
 }
 
-const FiniteElement *FiniteElementSpace::GetTypicalFaceElement() const
-{
-   return fec->FiniteElementForGeometry(mesh->GetTypicalFaceGeometry());
-}
-
 const FiniteElement *FiniteElementSpace::GetEdgeElement(int i,
                                                         int variant) const
 {

fem/fespace.hpp

@@ -839,7 +839,7 @@ public:
        Note: For vector-valued elements, the results pads up the range dimension
        to the spatial dimension. E.g., consider a stack of 5 vector-valued
        elements each representing 2D vectors, living in a 3 dimensional space.
-       Then this function would give 15, not 10.
+       Then this fucntion would give 15, not 10.
        */
    int GetVectorDim() const;
 
@@ -1323,24 +1323,12 @@ public:
         associated with i'th boundary face in the mesh object. */
    const FiniteElement *GetBE(int i) const;
 
-   /// @brief Return a "typical" boundary element.
-   ///
-   /// This can be used in situations where the local mesh partition may be
-   /// empty.
-   const FiniteElement *GetTypicalBE() const;
-
    /** @brief Returns pointer to the FiniteElement in the FiniteElementCollection
         associated with i'th face in the mesh object.  Faces in this case refer
         to the MESHDIM-1 primitive so in 2D they are segments and in 1D they are
         points.*/
    const FiniteElement *GetFaceElement(int i) const;
 
-   /// @brief Return a "typical" face element.
-   ///
-   /// This can be used in situations where the local mesh partition may be
-   /// empty.
-   const FiniteElement *GetTypicalFaceElement() const;
-
    /** @brief Returns pointer to the FiniteElement in the FiniteElementCollection
         associated with i'th edge in the mesh object. */
    const FiniteElement *GetEdgeElement(int i, int variant = 0) const;

fem/gridfunc.cpp

@@ -345,6 +345,27 @@ void GridFunction::ComputeFlux(BilinearFormIntegrator &blfi,
    }
 }
 
+int GridFunction::VectorDim() const
+{
+   const FiniteElement *fe = fes->GetTypicalFE();
+   if (!fe || fe->GetRangeType() == FiniteElement::SCALAR)
+   {
+      return fes->GetVDim();
+   }
+   return fes->GetVDim()*std::max(fes->GetMesh()->SpaceDimension(),
+                                  fe->GetRangeDim());
+}
+
+int GridFunction::CurlDim() const
+{
+   const FiniteElement *fe = fes->GetTypicalFE();
+   if (!fe || fe->GetRangeType() == FiniteElement::SCALAR)
+   {
+      return 2 * fes->GetMesh()->SpaceDimension() - 3;
+   }
+   return fes->GetVDim()*fe->GetCurlDim();
+}
+
 void GridFunction::GetTrueDofs(Vector &tv) const
 {
    const SparseMatrix *R = fes->GetRestrictionMatrix();
@@ -2029,18 +2050,6 @@ void GridFunction::AccumulateAndCountBdrValues(
    Coefficient *coeff[], VectorCoefficient *vcoeff, const Array<int> &attr,
    Array<int> &values_counter)
 {
-   if (vcoeff)
-   {
-      MFEM_VERIFY(fes->GetVDim() == vcoeff->GetVDim(),
-                  "vcoeff vdim != fes VDim");
-      MFEM_VERIFY(fes->GetTypicalBE()->GetMapType() == FiniteElement::VALUE &&
-                  fes->GetTypicalBE()->GetRangeType() ==
-                  FiniteElement::SCALAR,
-                  "Can only call ProjectBdrCoefficient on scalar value-type "
-                  "boundary elements. "
-                  "Did you intended to call ProjectBdrCoefficientNormal or "
-                  "ProjectBdrCoefficientTangent for vector finite elements?");
-   }
    Array<int> vdofs;
    Vector vc;
 
@@ -2193,9 +2202,6 @@ void GridFunction::AccumulateAndCountBdrTangentValues(
    VectorCoefficient &vcoeff, const Array<int> &bdr_attr,
    Array<int> &values_counter)
 {
-   MFEM_VERIFY(fes->GetTypicalBE()->GetPhysRangeDim(
-                  fes->GetMesh()->SpaceDimension()) == vcoeff.GetVDim(),
-               "vcoeff vdim != PhysRangeDim");
    const FiniteElement *fe;
    ElementTransformation *T;
    Array<int> dofs;
@@ -2349,9 +2355,6 @@ void GridFunction::ProjectDeltaCoefficient(DeltaCoefficient &delta_coeff,
 
 void GridFunction::ProjectCoefficient(Coefficient &coeff, ProjectType type)
 {
-   MFEM_VERIFY(
-      VectorDim() == 1,
-      "Cannot project scalar Coefficient onto vector GridFunction");
    DeltaCoefficient *delta_c = dynamic_cast<DeltaCoefficient *>(&coeff);
    DofTransformation doftrans;
    Array<int> vdofs;
@@ -2627,7 +2630,6 @@ void GridFunction::ProjectCoefficient(
 void GridFunction::ProjectCoefficient(VectorCoefficient &vcoeff,
                                       ProjectType type)
 {
-   MFEM_VERIFY(VectorDim() == vcoeff.GetVDim(), "vcoeff vdim != VectorDim()");
    Array<int> vdofs;
    Vector vals;
    DofTransformation doftrans;
@@ -2943,7 +2945,6 @@ void GridFunction::ProjectCoefficientElementL2(VectorCoefficient &vcoeff)
 void GridFunction::ProjectCoefficient(
    VectorCoefficient &vcoeff, Array<int> &dofs)
 {
-   MFEM_VERIFY(VectorDim() == vcoeff.GetVDim(), "vcoeff vdim != VectorDim()");
    int el = -1;
    ElementTransformation *T = NULL;
    const FiniteElement *fe = NULL;
@@ -2973,7 +2974,6 @@ void GridFunction::ProjectCoefficient(
 
 void GridFunction::ProjectCoefficient(VectorCoefficient &vcoeff, int attribute)
 {
-   MFEM_VERIFY(VectorDim() == vcoeff.GetVDim(), "vcoeff vdim != VectorDim()");
    int i;
    Array<int> vdofs;
    Vector vals;
@@ -3030,14 +3030,9 @@ void GridFunction::ProjectCoefficient(Coefficient *coeff[])
    }
 }
 
-void GridFunction::ProjectDiscCoefficient(
-   std::variant<Coefficient*, VectorCoefficient*> coeff, Array<int> &dof_attr)
+void GridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff,
+                                          Array<int> &dof_attr)
 {
-   std::visit([&](auto* c)
-   {
-      MFEM_VERIFY(VectorDim() == c->GetVDim(), "coeff vdim != VectorDim()");
-   }, coeff);
-
    Array<int> vdofs;
    Vector vals;
 
@@ -3051,10 +3046,7 @@ void GridFunction::ProjectDiscCoefficient(
    {
       fes->GetElementVDofs(i, vdofs);
       vals.SetSize(vdofs.Size());
-      std::visit([&](auto* c)
-      {
-         fes->GetFE(i)->Project(*c, *fes->GetElementTransformation(i), vals);
-      }, coeff);
+      fes->GetFE(i)->Project(coeff, *fes->GetElementTransformation(i), vals);
 
       // the values in shared dofs are determined from the element with maximal
       // attribute
@@ -3070,15 +3062,17 @@ void GridFunction::ProjectDiscCoefficient(
    }
 }
 
+void GridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff)
+{
+   Array<int> dof_attr;
+   ProjectDiscCoefficient(coeff, dof_attr);
+}
+
 void GridFunction::ProjectDiscCoefficient(Coefficient &coeff, AvgType type)
 {
    // Harmonic  (x1 ... xn) = [ (1/x1 + ... + 1/xn) / n ]^-1.
    // Arithmetic(x1 ... xn) = (x1 + ... + xn) / n.
 
-   MFEM_VERIFY(
-      VectorDim() == 1,
-      "Cannot project a scalar coefficient onto a vector GridFunction");
-
    Array<int> zones_per_vdof;
    AccumulateAndCountZones(coeff, type, zones_per_vdof);
 
@@ -3088,7 +3082,6 @@ void GridFunction::ProjectDiscCoefficient(Coefficient &coeff, AvgType type)
 void GridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff,
                                           AvgType type)
 {
-   MFEM_VERIFY(VectorDim() == coeff.GetVDim(), "coeff vdim != VectorDim()");
    Array<int> zones_per_vdof;
    AccumulateAndCountZones(coeff, type, zones_per_vdof);
 
@@ -3146,16 +3139,12 @@ void GridFunction::ProjectBdrCoefficient(Coefficient *coeff[],
 void GridFunction::ProjectBdrCoefficientNormal(
    Coefficient *coeff, VectorCoefficient *vcoeff, const Array<int> &bdr_attr)
 {
-   MFEM_VERIFY(fes->GetVDim() == 1, "fespace VDim != 1");
-   MFEM_VERIFY(fes->GetTypicalBE()->GetRangeType() == FiniteElement::SCALAR &&
-               fes->GetTypicalBE()->GetMapType() == FiniteElement::INTEGRAL,
-               "Not an RT FE space!");
-   if (vcoeff)
+   if (fes->GetNBE() > 0)
    {
-      MFEM_VERIFY(vcoeff->GetVDim() == fes->GetMesh()->SpaceDimension(),
-                  "vcoeff vdim (" << vcoeff->GetVDim()
-                  << ") != SpaceDimension ("
-                  << fes->GetMesh()->SpaceDimension() << ")");
+      // TODO: Replace this by GetTypicalBdrElement() once implemented
+      const FiniteElement *be = fes->GetBE(0);
+      MFEM_VERIFY(be->GetRangeType() == FiniteElement::SCALAR &&
+                  be->GetMapType() == FiniteElement::INTEGRAL, "Not an RT FE space!");
    }
 
    // implementation for the case when the face dofs are scaled point
@@ -5768,4 +5757,4 @@ std::pair<real_t, real_t> GridFunction::EstimateFunctionMaximum(
    return std::make_pair(global_max_lower, global_max_upper);
 }
 
-}
+}

fem/gridfunc.hpp

@@ -23,7 +23,6 @@
 #include <limits>
 #include <ostream>
 #include <string>
-#include <variant>
 
 namespace mfem
 {
@@ -80,18 +79,10 @@ protected:
                         bool wcoef,
                         int subdomain);
 
-   /** @brief Project a discontinuous (vector) coefficient as a grid function on
-       a continuous finite element space. Return in dof_attr the maximal
-       attribute of the elements containing each degree of freedom. */
-   virtual void ProjectDiscCoefficient(
-      std::variant<Coefficient*, VectorCoefficient*> coeff, Array<int> &dof_attr);
-
-   /** @brief Project a discontinuous (vector) coefficient as a grid function on
-       a continuous finite element space. The values in shared dofs are
-       determined from the element with maximal attribute. */
-   virtual void ProjectDiscCoefficient(
-      std::variant<Coefficient*, VectorCoefficient*> coeff)
-   { Array<int> dof_attr; ProjectDiscCoefficient(coeff, dof_attr); };
+   /** Project a discontinuous vector coefficient in a continuous space and
+       return in dof_attr the maximal attribute of the elements containing each
+       degree of freedom. */
+   void ProjectDiscCoefficient(VectorCoefficient &coeff, Array<int> &dof_attr);
 
    /** Helper function for ProjectCoefficientElementL2 */
    void ProjectCoefficientElementL2_(Coefficient &coeff, Vector &sol, Vector &Va);
@@ -159,13 +150,11 @@ public:
 
    FiniteElementCollection *OwnFEC() { return fec_owned; }
 
-   /** @brief Shortcut for calling FiniteElementSpace::GetVectorDim() on the
-       underlying #fes */
-   int VectorDim() const { return fes->GetVectorDim(); }
+   /// Shortcut for calling FiniteElementSpace::GetVectorDim() on the underlying #fes
+   int VectorDim() const;
 
-   /** @brief Shortcut for calling FiniteElementSpace::GetCurlDim() on the
-       underlying #fes */
-   int CurlDim() const { return fes->GetCurlDim(); }
+   /// Shortcut for calling FiniteElementSpace::GetCurlDim() on the underlying #fes
+   int CurlDim() const;
 
    /// Read only access to the (optional) internal true-dof Vector.
    const Vector &GetTrueVector() const
@@ -524,17 +513,10 @@ public:
        but using an array of scalar coefficients for each component. */
    void ProjectCoefficient(Coefficient *coeff[]);
 
-   /** @brief Project a discontinuous coefficient as a grid function on
-       a continuous finite element space. The values in shared dofs are
-       determined from the element with maximal attribute. */
-   virtual void ProjectDiscCoefficient(Coefficient &coeff)
-   { ProjectDiscCoefficient(&coeff); }
-
    /** @brief Project a discontinuous vector coefficient as a grid function on
        a continuous finite element space. The values in shared dofs are
        determined from the element with maximal attribute. */
-   virtual void ProjectDiscCoefficient(VectorCoefficient &coeff)
-   { ProjectDiscCoefficient(&coeff); }
+   virtual void ProjectDiscCoefficient(VectorCoefficient &coeff);
 
    enum AvgType {ARITHMETIC, HARMONIC};
    /** @brief Projects a discontinuous coefficient so that the values in shared
@@ -1989,7 +1971,6 @@ public:
 
    void Eval(Vector &v, ElementTransformation &T,
              const IntegrationPoint &ip) override;
-   using VectorCoefficient::Eval;
 
    virtual ~VectorExtrudeCoefficient() { }
 };

fem/gslib.cpp

@@ -490,7 +490,7 @@ void FindPointsGSLIB::FindPointsOnDevice(const Vector &point_pos,
    }
    DEV.find_device = true;
 
-   const unsigned int id = gsl_comm->id, np = gsl_comm->np;
+   const int id = gsl_comm->id, np = gsl_comm->np;
 
    gsl_mfem_ref.SetSize(points_cnt * dim);
    gsl_mfem_elem.SetSize(points_cnt);
@@ -652,7 +652,7 @@ void FindPointsGSLIB::FindPointsOnDevice(const Vector &point_pos,
          {
             const int pp = hash_offset[i];
             /* don't send back to where it just came from */
-            if (static_cast<unsigned>(pp) == p->proc)
+            if (pp == p->proc)
             {
                continue;
             }
@@ -1068,7 +1068,7 @@ void FindPointsGSLIB::InterpolateOnDevice(const Vector &field_in_evec,
          sarray_transfer(struct evalOutPt_t, &outpt, proc, 1, cr);
 
          opt = (evalOutPt_t *)outpt.ptr;
-         for (size_t index = 0; index < outpt.n; index++)
+         for (int index = 0; index < outpt.n; index++)
          {
             int idx = ordering == Ordering::byNODES ?
                       opt->index + i*points_cnt :
@@ -1413,7 +1413,7 @@ void FindPointsGSLIB::SetupSplitMeshesAndIntegrationRules(const int order)
 {
    MFEM_VERIFY(mesh, "Setup FindPointsGSLIB with mesh first.");
    const int dof1D = order+1;
-   dim = mesh->Dimension();
+   const int dim = mesh->Dimension();
 
    SetupSplitMeshes();
    if (dim == 2)
@@ -2254,8 +2254,7 @@ void FindPointsGSLIB::DistributeInterpolatedValues(const Vector &int_vals,
       sarray_transfer(struct out_pt, outpt, proc, 1, cr);
 
       // Store received data
-      MFEM_VERIFY(outpt->n == static_cast<size_t>(points_cnt),
-                  "Incompatible size. Number of points "
+      MFEM_VERIFY(outpt->n == points_cnt, "Incompatible size. Number of points "
                   "received does not match the number of points originally "
                   "found using FindPoints.");
 

fem/gslib.hpp

@@ -202,19 +202,13 @@ protected:
                             const int dof1dsol, const int ordering);
 
 public:
-   /// Serial constructor
    FindPointsGSLIB();
-
-   /// Serial constructor + setup with given Mesh (see \ref Setup)
    FindPointsGSLIB(Mesh &mesh_in, const double bb_t = 0.1,
                    const double newt_tol = 1.0e-12,
                    const int npt_max = 256);
 
 #ifdef MFEM_USE_MPI
-   /// Constructor for ParMesh
    FindPointsGSLIB(MPI_Comm comm_);
-
-   /// Constructor + setup with given ParMesh (see \ref Setup)
    FindPointsGSLIB(ParMesh &mesh_in, const double bb_t = 0.1,
                    const double newt_tol = 1.0e-12,
                    const int npt_max = 256);

fem/gslib/findpts_local_2.cpp

@@ -254,7 +254,7 @@ get_edge(const double *elx[2], const double *wtend, int ei,
       edge.dxdn[d] = workspace + (2 + d) * pN; //dxdn and dydn at DOFs along edge
    }
 
-   if (static_cast<unsigned>(side_init) != (1u << ei))
+   if (side_init != (1u << ei))
    {
 #define ELX(d, j, k) elx[d][j + k * pN] // assumes lexicographic ordering
       for (int d = 0; d < 2; ++d)

fem/gslib/findpts_local_3.cpp

@@ -294,7 +294,7 @@ get_face(const double *elx[3], const double *wtend, int fi, double *workspace,
       face.dxdn[d] = workspace+(3+d)*p_Nfr;
    }
 
-   if (static_cast<unsigned>(side_init) != (1u << fi))
+   if (side_init != (1u << fi))
    {
       const int e_stride[3] = {1, pN, pN*pN};
 #define ELX(d, j, k, l) elx[d][j*e_stride[d1]+k*e_stride[d2]+l*e_stride[dn]]
@@ -342,7 +342,7 @@ get_edge(const double *elx[3], const double *wtend, int ei, double *workspace,
 
    if (jidx >= 3*pN) { return edge; }
 
-   if (static_cast<unsigned>(side_init) != (64u << ei))
+   if (side_init != (64u << ei))
    {
       const int e_stride[3] = {1, pN, pN*pN};
 #define ELX(d, j, k, l) elx[d][j*e_stride[de]+k*e_stride[dn1]+l*e_stride[dn2]]

fem/intrules.hpp

@@ -43,52 +43,56 @@ public:
       index = i;
    }
 
-   void Set3w(const real_t x1, const real_t x2, const real_t x3, const real_t w)
-   { x = x1; y = x2; z = x3; weight = w; }
-   void Set2w(const real_t x1, const real_t x2, const real_t w)
-   { x = x1; y = x2; weight = w; }
-   void Set1w(const real_t x1, const real_t w)
-   { x = x1; weight = w; }
-
-   void Set3w(const real_t *p) { Set3w(p[0], p[1], p[2], p[3]); }
-   void Set2w(const real_t *p) { Set2w(p[0], p[1], p[2]); }
-   void Set1w(const real_t *p) { Set1w(p[0], p[1]); }
-
-   void Set3(const real_t x1, const real_t x2, const real_t x3)
-   { x = x1; y = x2; z = x3; }
-   void Set2(const real_t x1, const real_t x2)
-   { x = x1; y = x2; }
-   void Set1(const real_t x1)
-   { x = x1; }
-
-   void Set3(const real_t *p) { Set3(p[0], p[1], p[2]); }
-   void Set2(const real_t *p) { Set2(p[0], p[1]); }
-   void Set1(const real_t *p) { Set1(p[0]); }
-
-   void Set(const real_t x1, const real_t x2, const real_t x3, const real_t w)
-   { Set3w(x1, x2, x3, w); }
-
    void Set(const real_t *p, const int dim)
    {
       MFEM_ASSERT(1 <= dim && dim <= 3, "invalid dim: " << dim);
-      switch (dim)
+      x = p[0];
+      if (dim > 1)
       {
-         case 3: Set3(p); break;
-         case 2: Set2(p); break;
-         case 1: Set1(p); break;
+         y = p[1];
+         if (dim > 2)
+         {
+            z = p[2];
+         }
       }
    }
 
    void Get(real_t *p, const int dim) const
    {
       MFEM_ASSERT(1 <= dim && dim <= 3, "invalid dim: " << dim);
-      switch (dim)
+      p[0] = x;
+      if (dim > 1)
       {
-         case 3: p[2] = z;
-         case 2: p[1] = y;
-         case 1: p[0] = x;
+         p[1] = y;
+         if (dim > 2)
+         {
+            p[2] = z;
+         }
       }
    }
+
+   void Set(const real_t x1, const real_t x2, const real_t x3, const real_t w)
+   { x = x1; y = x2; z = x3; weight = w; }
+
+   void Set3w(const real_t *p) { x = p[0]; y = p[1]; z = p[2]; weight = p[3]; }
+
+   void Set3(const real_t x1, const real_t x2, const real_t x3)
+   { x = x1; y = x2; z = x3; }
+
+   void Set3(const real_t *p) { x = p[0]; y = p[1]; z = p[2]; }
+
+   void Set2w(const real_t x1, const real_t x2, const real_t w)
+   { x = x1; y = x2; weight = w; }
+
+   void Set2w(const real_t *p) { x = p[0]; y = p[1]; weight = p[2]; }
+
+   void Set2(const real_t x1, const real_t x2) { x = x1; y = x2; }
+
+   void Set2(const real_t *p) { x = p[0]; y = p[1]; }
+
+   void Set1w(const real_t x1, const real_t w) { x = x1; weight = w; }
+
+   void Set1w(const real_t *p) { x = p[0]; weight = p[1]; }
 };
 
 /// Class for an integration rule - an Array of IntegrationPoint.

fem/lininteg.hpp

@@ -164,8 +164,8 @@ private:
 
 public:
    /// Constructs the domain integrator $ (Q, \nabla v) $
-   DomainLFGradIntegrator(VectorCoefficient &QF, const IntegrationRule *ir = NULL)
-      : DeltaLFIntegrator(QF, ir), Q(QF) { }
+   DomainLFGradIntegrator(VectorCoefficient &QF)
+      : DeltaLFIntegrator(QF), Q(QF) { }
 
    bool SupportsDevice() const override { return true; }
 

fem/pgridfunc.cpp

@@ -545,8 +545,6 @@ void ParGridFunction::GetElementDofValues(int el, Vector &dof_vals) const
 
 void ParGridFunction::ProjectCoefficient(Coefficient &coeff, ProjectType type)
 {
-   MFEM_VERIFY(VectorDim() == 1,
-               "Cannot project scalar coefficient onto vector ParGridFunction");
    DeltaCoefficient *delta_c = dynamic_cast<DeltaCoefficient *>(&coeff);
 
    if (delta_c == NULL)
@@ -717,8 +715,7 @@ void ParGridFunction::ProjectCoefficientElementL2(VectorCoefficient &vcoeff)
 }
 
 
-void ParGridFunction::ProjectDiscCoefficient(
-   std::variant<Coefficient*, VectorCoefficient*> coeff)
+void ParGridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff)
 {
    // local maximal element attribute for each dof
    Array<int> ldof_attr;
@@ -764,9 +761,6 @@ void ParGridFunction::ProjectDiscCoefficient(
 
 void ParGridFunction::ProjectDiscCoefficient(Coefficient &coeff, AvgType type)
 {
-   MFEM_VERIFY(
-      VectorDim() == 1,
-      "Cannot project scalar coefficient onto a vector ParGridFunction");
    // Harmonic  (x1 ... xn) = [ (1/x1 + ... + 1/xn) / n ]^-1.
    // Arithmetic(x1 ... xn) = (x1 + ... + xn) / n.
 
@@ -792,8 +786,6 @@ void ParGridFunction::ProjectDiscCoefficient(VectorCoefficient &vcoeff,
    // Harmonic  (x1 ... xn) = [ (1/x1 + ... + 1/xn) / n ]^-1.
    // Arithmetic(x1 ... xn) = (x1 + ... + xn) / n.
 
-   MFEM_VERIFY(VectorDim() == vcoeff.GetVDim(), "vcoeff vdim != VectorDim()");
-
    // Number of zones that contain a given dof.
    Array<int> zones_per_vdof;
    AccumulateAndCountZones(vcoeff, type, zones_per_vdof);
@@ -866,12 +858,6 @@ void ParGridFunction::ProjectBdrCoefficient(
 #endif
 }
 
-void ParGridFunction::ProjectBdrCoefficient(VectorCoefficient &vcoeff,
-                                            const Array<int> &attr)
-{
-   ProjectBdrCoefficient(NULL, &vcoeff, attr);
-}
-
 void ParGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient &vcoeff,
                                                    const Array<int> &bdr_attr)
 {

fem/pgridfunc.hpp

@@ -63,12 +63,6 @@ protected:
    void ProjectBdrCoefficient(Coefficient *coeff[], VectorCoefficient *vcoeff,
                               const Array<int> &attr);
 
-   /** @brief Project a discontinuous (vector) coefficient as a grid function on
-       a continuous finite element space. The values in shared dofs are
-       determined from the element with maximal attribute. */
-   virtual void ProjectDiscCoefficient(
-      std::variant<Coefficient*, VectorCoefficient*> coeff) override;
-
 public:
    ParGridFunction() { pfes = NULL; }
 
@@ -274,6 +268,11 @@ public:
                            ProjectType type = ProjectType::DEFAULT) override;
 
    using GridFunction::ProjectDiscCoefficient;
+   /** @brief Project a discontinuous vector coefficient as a grid function on
+       a continuous finite element space. The values in shared dofs are
+       determined from the element with maximal attribute. */
+   void ProjectDiscCoefficient(VectorCoefficient &coeff) override;
+
    void ProjectDiscCoefficient(Coefficient &coeff, AvgType type) override;
 
    void ProjectDiscCoefficient(VectorCoefficient &vcoeff, AvgType type) override;
@@ -281,7 +280,8 @@ public:
    using GridFunction::ProjectBdrCoefficient;
 
    void ProjectBdrCoefficient(VectorCoefficient &vcoeff,
-                              const Array<int> &attr) override;
+                              const Array<int> &attr) override
+   { ProjectBdrCoefficient(NULL, &vcoeff, attr); }
 
    void ProjectBdrCoefficient(Coefficient *coeff[],
                               const Array<int> &attr) override

general/arrays_by_name.hpp

@@ -14,7 +14,6 @@
 
 #include "../config/config.hpp"
 #include "array.hpp"
-#include "text.hpp"
 
 #include <iostream>
 #include <map>
@@ -248,8 +247,7 @@ inline void ArraysByName<T>::Print(std::ostream &os, int width) const
    os << data.size() << '\n';
    for (auto const &it : data)
    {
-      // Note: The method Load() can read any string formatted with std::quoted.
-      os << std::quoted(it.first) << '\n' << it.second.Size() << '\n';
+      os << '"' << it.first << '"' << '\n' << it.second.Size() << '\n';
       it.second.Print(os, width > 0 ? width : it.second.Size());
    }
 }
@@ -260,36 +258,40 @@ void ArraysByName<T>::Load(std::istream &in)
    int NumArrays;
    in >> NumArrays;
 
-   for (int i = 0; i < NumArrays; i++)
+   std::string ArrayLine, ArrayName;
+   for (int i=0; i < NumArrays; i++)
    {
       in >> std::ws;
-      // Read the name:
-      // - If the stream 'in' starts with " then parse it with the function
-      //   parse_quoted_string() from text.hpp. In this case, the name can be
-      //   empty. Note: this case allows for reading any string formatted using
-      //   std::quoted, e.g. as in the method Print().
-      // - If the name does not start with " then the name ends with the first
-      //   white space character (and the white space character is not included
-      //   in the name). Since white space characters are skipped before reading
-      //   the name, there will be at least one non-white-space character in the
-      //   name in this case.
-      std::string ArrayName;
-      if (in.peek() == '"')
+      getline(in, ArrayLine);
+
+      std::size_t q0 = ArrayLine.find('"');
+      std::size_t q1 = ArrayLine.rfind('"');
+
+      if (q0 != std::string::npos && q1 > q0)
       {
-         if (parse_quoted_string(ArrayName, in) != 0)
-         {
-            MFEM_ABORT("error parsing input!");
-         }
+         // Locate set name between first and last double quote
+         ArrayName = ArrayLine.substr(q0+1,q1-q0-1);
       }
       else
       {
-         in >> ArrayName;
-         MFEM_VERIFY(in.good(), "error parsing input!");
+         // If no double quotes found locate set name using white space
+         q1 = ArrayLine.find(' ');
+         ArrayName = ArrayLine.substr(0,q1-1);
+      }
+      if (q1+2 < ArrayLine.size())
+      {
+         // Read the remainder of the line which contains the array data
+         std::istringstream ArrayDataStream(ArrayLine.substr(q1+2,
+                                                             ArrayLine.size()));
+         data[ArrayName].Load(ArrayDataStream, 0);
+      }
+      else
+      {
+         // Read the array data starting on the next line
+         data[ArrayName].Load(in, 0);
       }
-
-      // Read the array
-      data[ArrayName].Load(in);
    }
+
 }
 
 }

general/text.hpp

@@ -50,48 +50,6 @@ inline void filter_dos(std::string &line)
    }
 }
 
-/** @brief Read a string formatted using std::quoted. Return nonzero on error.
-
-    The stream @a in must begin with @a delim. After clearing @a result and
-    extracting the opening @a delim, characters are extracted from @a in and
-    processed as follows:
-    - if the character is @a delim, return 0;
-    - if the character is different from @a escape, it is appended to @a result;
-    - if the character is @a escape, the next character from @a in is extracted
-      and if it is one of @a delim or @a escape, it is appended to @a result;
-      otherwise, both @a escape and the character after it are appended to
-      @a result; note that the latter case is not possible if the input was
-      formatted with std::quoted with the same @a delim and @a escape
-      characters.
-
-    If the stream @a in does not begin with @a delim, error code 1 is returned.
-    If reading the stream fails, error code 2 is returned. On success, zero is
-    returned and the closing @a delim character is the last character extracted
-    from @a in. */
-inline int parse_quoted_string(std::string &result, std::istream &in,
-                               char delim = '"', char escape = '\\')
-{
-   using tt = std::string::traits_type;  // std::char_traits<char>
-   auto equal = [](tt::int_type c1, tt::char_type c2) -> bool
-   {
-      return tt::eq_int_type(c1, tt::to_int_type(c2));
-   };
-   result.clear();
-   if (!equal(in.peek(), delim)) { return 1; }
-   in.get();  // extract delim
-   for (auto c = in.get(); !equal(c, delim); c = in.get())
-   {
-      if (equal(c, escape))
-      {
-         c = in.get();
-         if (!equal(c, escape) && !equal(c, delim)) { result += escape; }
-      }
-      if (!in) { return 2; }
-      result += tt::to_char_type(c);
-   }
-   return 0;
-}
-
 /// Convert an integer to a 0-padded string with the given number of @a digits
 inline std::string to_padded_string(int i, int digits)
 {

linalg/CMakeLists.txt

@@ -23,7 +23,6 @@ list(APPEND SRCS
   complex_operator.cpp
   constraints.cpp
   densemat.cpp
-  eigensolvers.cpp
   filteredsolver.cpp
   handle.cpp
   matrix.cpp
@@ -56,7 +55,6 @@ list(APPEND HDRS
   dinvariants.hpp
   dtensor.hpp
   dual.hpp
-  eigensolvers.hpp
   filteredsolver.hpp
   handle.hpp
   invariants.hpp
@@ -103,11 +101,6 @@ if (MFEM_USE_MPI)
   endif()
 endif()
 
-if (MFEM_USE_ARPACK)
-  list(APPEND SRCS arpack.cpp)
-  list(APPEND HDRS arpack.hpp)
-endif()
-
 if (MFEM_USE_SUNDIALS)
   list(APPEND SRCS sundials.cpp)
   list(APPEND HDRS sundials.hpp)

linalg/arpack.hpp

@@ -1,271 +0,0 @@
-// Copyright (c) 2010-2025, 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>
-
-#ifdef MFEM_USE_MPI
-#include <mpi.h>
-#include "hypre.hpp"
-#endif
-
-#include "operator.hpp"
-
-#define SSAUPD  ssaupd_
-#define SSEUPD  sseupd_
-#define DSAUPD  dsaupd_
-#define DSEUPD  dseupd_
-
-#ifdef MFEM_USE_MPI
-#define PSSAUPD pssaupd_
-#define PSSEUPD psseupd_
-#define PDSAUPD pdsaupd_
-#define PDSEUPD pdseupd_
-#endif
-
-extern "C" void SSAUPD(int *ido, char *bmat, int *n,
-                       char *which, int *nev, float *tol, float *resid,
-                       int *ncv, float *v, int *ldv,
-                       int *iparam, int *ipntr,
-                       float *workd, float *workl, int *lworkl, int *info);
-
-extern "C" void SSEUPD(int *, char *, int *, float *,
-                       float *, int *, float *, char *, int *, char *,
-                       int *, float *, float *, int *, float *,
-                       int *, int  *, int *, float *,
-                       float *, int *, int *);
-
-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 PSSAUPD(int *comm, int *ido, char *bmat, int *n,
-                        char *which, int *nev, float *tol, float *resid,
-                        int *ncv, float *v, int *ldv,
-                        int *iparam, int *ipntr,
-                        float *workd, float *workl, int *lworkl, int *info);
-
-extern "C" void PSSEUPD(int *comm, int *, char *, int *, float *,
-                        float *, int *, float *, char *, int *, char *,
-                        int *, float *, float *, int *, float *,
-                        int *, int  *, int *, float *,
-                        float *, int *, int *);
-
-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
-{
-
-/// Wrapper for the ARPACK routine SSAUPD or DSAUPD
-class ArPackSAUPD : public SymEigensolver, public SymGenEigensolver
-{
-public:
-
-   ArPackSAUPD();
-   virtual ~ArPackSAUPD();
-
-   /** 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(real_t tol) override         {      tol_ = tol;      }
-   inline void SetMaxIter(int max_iter) override   { max_iter_ = max_iter; }
-   inline void SetPrintLevel(int logging) override {  logging_ = logging;  }
-   inline void SetShift(real_t sigma)              {    sigma_ = sigma;    }
-   inline void SetNumModes(int num_eigs) override  {      nev_ = num_eigs; }
-
-   virtual void SetSolver(Solver & solver);
-   virtual void SetOperator(const Operator & A) override;
-   virtual void SetMassMatrix(const Operator & M);
-   virtual void SetOperators(const Operator & A, const Operator & B) override
-   { SetOperator(A); SetMassMatrix(B); }
-
-   void Solve() override;
-
-   virtual int GetNumConverged() const override { return iparam_[4]; }
-
-   /// Collect the converged eigenvalues
-   virtual void GetEigenvalues(Array<real_t> & eigenvalues) const override;
-
-   /// Extract a single eigenvector
-   virtual const Vector & GetEigenvector(unsigned int i) const override;
-
-   /// Transfer ownership of the converged eigenvectors
-   Vector ** StealEigenvectors() override;
-
-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
-
-   real_t tol_;     // relative accuracy bound for Ritz values
-   real_t sigma_;   // eigenvalue shift parameter
-
-   int    * select_;// workspace used during eigenvalue computation
-   real_t * dv_;    // Ritz values
-   real_t * v_;     // ncv Lanczos basis vectors
-   real_t * resid_; // residual vector
-   real_t * workd_; // work array for 3 vectors used in Arnoldi iteration
-   real_t * workl_; // work array
-
-   // Operators and Vectors needed outside of ARPACK
-   Solver   * solver_;
-   const Operator * A_;
-   const Operator * B_;
-
-   Vector * w_;
-   Vector * x_;
-   Vector * y_;
-   Vector * z_;
-
-   mutable Vector ** eigenvectors_;
-
-   std::string solverName_;
-
-   void reverseComm();
-
-   int reverseCommMode1();
-   int reverseCommMode2();
-   int reverseCommMode3();
-   int reverseCommMode4();
-   int reverseCommMode5();
-
-   virtual void prepareEigenvectors() const;
-
-   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 ArPackPSAUPD : public ArPackSAUPD
-{
-public:
-   ArPackPSAUPD(MPI_Comm comm);
-   virtual ~ArPackPSAUPD() {}
-
-   void SetOperator(const Operator & A);
-   void SetMassMatrix(const Operator & M);
-
-   /// Collect the converged eigenvalues
-   void GetEigenvalues(Array<real_t> & eigenvalues) const;
-
-   /// Extract a single eigenvector
-   const Vector & GetEigenvector(unsigned int i) const;
-
-   /// Transfer ownership of the converged eigenvectors
-   // HypreParVector ** StealEigenvectors();
-   Vector ** StealEigenvectors();
-
-protected:
-
-   void prepareEigenvectors() const;
-
-private:
-
-   MPI_Comm comm_;
-   MPI_Fint commf_; // Fortran style MPI communicator
-   int numProcs_;   // Number of processors
-
-   mutable 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/eigensolvers.cpp

@@ -1,20 +0,0 @@
-// Copyright (c) 2010-2025, 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 "eigensolvers.hpp"
-
-using namespace std;
-
-namespace mfem
-{
-
-};

linalg/eigensolvers.hpp

@@ -1,396 +0,0 @@
-// Copyright (c) 2010-2025, 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"
-
-namespace mfem
-{
-
-/// Abstract Eigenequation
-/// Defines the operator of the linear eigenvalue equation
-///    A x_i = lambda_i x_i
-/// Where A is a real-valued operator, the lambda_i are the eigenvalues,
-/// and x_i are the eigenvectors.
-class Eigenequation
-{
-protected:
-   Eigenequation() = default;
-
-public:
-   virtual ~Eigenequation() = default;
-
-   /// @brief Set the operator A of the eigenvalue equation
-   virtual void SetOperator(const Operator & A) = 0;
-};
-
-/// Abstract Complex-valued Eigenequation
-/// Defines the operator of the linear eigenvalue equation
-///    A x_i = lambda_i x_i
-/// Where A is a complex-valued operator, the lambda_i are the eigenvalues,
-/// and x_i are the eigenvectors.
-class ComplexEigenequation
-{
-protected:
-   ComplexEigenequation() = default;
-
-public:
-   virtual ~ComplexEigenequation() = default;
-
-   /// @brief Set the real and imaginary parts of the operator A
-   virtual void SetOperator(const Operator & Ar, const Operator & Ai) = 0;
-};
-
-/// Abstract Generalized Eigenequation
-/// Defines the operator of the linear eigenvalue equation
-///    A x_i = lambda_i B x_i
-/// Where A and B are real-valued operators, the lambda_i are the eigenvalues,
-/// and x_i are the eigenvectors.
-class GenEigenequation
-{
-protected:
-   GenEigenequation() = default;
-
-public:
-   virtual ~GenEigenequation() = default;
-
-   /// @brief Set the operators A and B of the generalized eigenvalue equation
-   virtual void SetOperators(const Operator & A, const Operator & B) = 0;
-};
-
-/// Abstract Complex-valued Generalized Eigenequation
-/// Defines the operator of the linear eigenvalue equation
-///    A x_i = lambda_i B x_i
-/// Where A and B are complex-valued operators, the lambda_i are the
-/// eigenvalues, and x_i are the eigenvectors.
-class ComplexGenEigenequation
-{
-protected:
-   ComplexGenEigenequation() = default;
-
-public:
-   virtual ~ComplexGenEigenequation() = default;
-
-   /// @brief Set the real and imaginary parts of the operators A and B
-   virtual void SetOperators(const Operator & Ar, const Operator & Ai,
-                             const Operator & Br, const Operator & Bi) = 0;
-};
-
-/// Abstract Eigensolver
-/// Computes eigenvalue/eigenvector pairs for the linear system
-///    A x_i = lambda_i x_i
-/// Where the lambda_i are the eigenvalues and x_i are the eigenvectors.
-class EigensolverBase
-{
-protected:
-   EigensolverBase() = default;
-
-public:
-   virtual ~EigensolverBase() = default;
-
-   /// @brief Stopping criteria based on numerical tolerance
-   ///
-   /// @note This may be defined differently by different solvers.
-   virtual void SetTol(real_t tol) = 0;
-
-   /// @brief Stopping criteria based on number of iterations required to
-   /// reach convergence.
-   ///
-   /// @note This may also be defined differently in different solvers.
-   virtual void SetMaxIter(int max_iter) = 0;
-
-   /// @brief Controls the type and amount of information printed to
-   ///        standard output.
-   virtual void SetPrintLevel(int logging) = 0;
-
-   /// @brief Set the number of desired eigenmodes to compute
-   virtual void SetNumModes(int num_eigs) = 0;
-
-   /// @brief Get the number of converged eigenmodes
-   virtual int GetNumConverged() const = 0;
-
-   /// @brief Perform the eigenvalue solve
-   virtual void Solve() = 0;
-};
-
-/// Symmetric Eigensolver
-/// If A^T = A the linear system must have real-valued eigenvalues
-/// and eigenvectors.
-class SymEigensolver : public EigensolverBase, public Eigenequation
-{
-protected:
-   SymEigensolver() = default;
-
-public:
-   virtual ~SymEigensolver() = default;
-
-   /// @brief Collect the converged eigenvalues
-   ///
-   /// The length of the array should equal the number of converged eigenvalues.
-   virtual void GetEigenvalues(Array<real_t> & eigenvalues) const = 0;
-
-   /// @brief Extract a single eigenvector
-   ///
-   /// The index i should be in the range [0, numConverged). The
-   virtual const Vector & GetEigenvector(unsigned int i) const = 0;
-
-   /// @brief Transfer ownership of the converged eigenvectors
-   ///
-   /// The array should contain numConverged vectors.
-   virtual Vector ** StealEigenvectors() = 0;
-};
-
-/// Symmetric Generalized Eigensolver
-/// If A^T = A and M^T = M the linear system must have real-valued eigenvalues
-/// and eigenvectors.
-class SymGenEigensolver : public EigensolverBase, public GenEigenequation
-{
-protected:
-   SymGenEigensolver() = default;
-
-public:
-   virtual ~SymGenEigensolver() = default;
-
-   /// @brief Collect the converged eigenvalues
-   ///
-   /// The length of the array should equal the number of converged eigenvalues.
-   virtual void GetEigenvalues(Array<real_t> & eigenvalues) const = 0;
-
-   /// @brief Extract a single eigenvector
-   ///
-   /// The index i should be in the range [0, numConverged). The
-   virtual const Vector & GetEigenvector(unsigned int i) const = 0;
-
-   /// @brief Transfer ownership of the converged eigenvectors
-   ///
-   /// The array should contain numConverged vectors.
-   virtual Vector ** StealEigenvectors() = 0;
-};
-
-/// Hermetian Eigensolver
-/// If A^H = A the linear system must have real-valued eigenvalues
-/// but may have complex-valued eigenvectors.
-class HermEigensolver : public EigensolverBase, public ComplexEigenequation
-{
-protected:
-   HermEigensolver() = default;
-
-public:
-   virtual ~HermEigensolver() = default;
-
-   /// @brief Collect the converged eigenvalues
-   ///
-   /// The length of the array should equal the number of converged eigenvalues.
-   virtual void GetEigenvalues(Array<real_t> & eigenvalues) const = 0;
-
-   /// @brief Extract a single eigenvector
-   ///
-   /// The index i should be in the range [0, 2*numConverged). The
-   /// vectors corresponding to even indices are the real parts of the
-   /// converged eigenvectors and the odd indices correspond to the
-   /// imaginary parts.
-   virtual const Vector & GetEigenvector(unsigned int i) const = 0;
-
-   /// @brief Transfer ownership of the converged eigenvectors
-   ///
-   /// The array should contain 2*numConverged vectors with the even
-   /// indices corresponding to the real parts of the converged
-   /// eigenvectors and the odd indices corresponding to the imaginary
-   /// parts.
-   virtual Vector ** StealEigenvectors() = 0;
-};
-
-/// Hermetian Generalized Eigensolver
-/// If A^H = A and M^H = M the linear system must have real-valued eigenvalues
-/// but may have complex-valued eigenvectors.
-class HermGenEigensolver :
-   public EigensolverBase, public ComplexGenEigenequation
-{
-protected:
-   HermGenEigensolver() = default;
-
-public:
-   virtual ~HermGenEigensolver() = default;
-
-   /// @brief Collect the converged eigenvalues
-   ///
-   /// The length of the array should equal the number of converged eigenvalues.
-   virtual void GetEigenvalues(Array<real_t> & eigenvalues) const = 0;
-
-   /// @brief Extract a single eigenvector
-   ///
-   /// The index i should be in the range [0, 2*numConverged). The
-   /// vectors corresponding to even indices are the real parts of the
-   /// converged eigenvectors and the odd indices correspond to the
-   /// imaginary parts.
-   virtual const Vector & GetEigenvector(unsigned int i) const = 0;
-
-   /// @brief Transfer ownership of the converged eigenvectors
-   ///
-   /// The array should contain 2*numConverged vectors with the even
-   /// indices corresponding to the real parts of the converged
-   /// eigenvectors and the odd indices corresponding to the imaginary
-   /// parts.
-   virtual Vector ** StealEigenvectors() = 0;
-};
-
-/// Non-Symmetric Eigensolver
-/// For general real-valued operators A the linear system must have
-/// eigenvalues and eigenvectors which form complex conjugate pairs.
-class NonSymEigensolver : public EigensolverBase, public Eigenequation
-{
-protected:
-   NonSymEigensolver() = default;
-
-public:
-   virtual ~NonSymEigensolver() = default;
-
-   /// @brief Collect the converged eigenvalues
-   ///
-   /// The length of the array should be the number of converged
-   /// eigenvalues. The complex-valued eigenvalues can be constructed
-   /// as: lambda_{2*j} = eig[2*j]+i*eig[2*j+1] and
-   /// lambda_{2*j+1} = eig[2*j]-i*eig[2*j+1]
-   /// With j in the range [0, numConverged/2)
-   virtual void GetEigenvalues(Array<real_t> & eig) const = 0;
-
-   /// @brief Extract a single eigenvector
-   ///
-   /// The index i should be in the range [0, numConverged). The
-   /// vectors corresponding to even indices are the real parts of the
-   /// converged eigenvectors and the odd indices correspond to the
-   /// imaginary parts. If needed, the complex conjugate pairs of
-   /// eigenvectors can be constructed in the same manner described
-   /// for the eigenvalues.
-   virtual const Vector & GetEigenvector(unsigned int i) const = 0;
-
-   /// @brief Transfer ownership of the converged eigenvectors
-   ///
-   /// The array should contain numConverged vectors with the even
-   /// indices corresponding to the real parts of the converged
-   /// eigenvectors and the odd indices corresponding to the imaginary
-   /// parts.
-   virtual Vector ** StealEigenvectors() = 0;
-};
-
-/// Non-Symmetric Eigensolver
-/// For general real-valued operators A and M the linear system must have
-/// eigenvalues and eigenvectors which form complex conjugate pairs.
-class NonSymGenEigensolver : public EigensolverBase, public GenEigenequation
-{
-protected:
-   NonSymGenEigensolver() = default;
-
-public:
-   virtual ~NonSymGenEigensolver() = default;
-
-   /// @brief Collect the converged eigenvalues
-   ///
-   /// The length of the array should be the number of converged
-   /// eigenvalues. The complex-valued eigenvalues can be constructed
-   /// as: lambda_{2*j} = eig[2*j]+i*eig[2*j+1] and
-   /// lambda_{2*j+1} = eig[2*j]-i*eig[2*j+1]
-   /// With j in the range [0, numConverged/2)
-   virtual void GetEigenvalues(Array<real_t> & eig) const = 0;
-
-   /// @brief Extract a single eigenvector
-   ///
-   /// The index i should be in the range [0, numConverged). The
-   /// vectors corresponding to even indices are the real parts of the
-   /// converged eigenvectors and the odd indices correspond to the
-   /// imaginary parts. If needed, the complex conjugate pairs of
-   /// eigenvectors can be constructed in the same manner described
-   /// for the eigenvalues.
-   virtual const Vector & GetEigenvector(unsigned int i) const = 0;
-
-   /// @brief Transfer ownership of the converged eigenvectors
-   ///
-   /// The array should contain numConverged vectors with the even
-   /// indices corresponding to the real parts of the converged
-   /// eigenvectors and the odd indices corresponding to the imaginary
-   /// parts.
-   virtual Vector ** StealEigenvectors() = 0;
-};
-
-/// Complex Eigensolver
-/// Can have arbitrary complex-valued eigenvalues and eigenvectors
-class ComplexEigensolver : public EigensolverBase, public ComplexEigenequation
-{
-protected:
-   ComplexEigensolver() = default;
-
-public:
-   virtual ~ComplexEigensolver() = default;
-
-   /// @brief Collect the converged eigenvalues
-   ///
-   /// The length of the array should be twice the number of converged
-   /// eigenvalues. The complex-valued eigenvalues can be constructed
-   /// as: lambda_j = eig[2*j]+i*eig[2*j+1]
-   virtual void GetEigenvalues(Array<real_t> & eig) const = 0;
-
-   /// @brief Extract a single eigenvector
-   ///
-   /// The index i should be in the range [0, 2*numConverged). The
-   /// vectors corresponding to even indices are the real parts of the
-   /// converged eigenvectors and the odd indices correspond to the
-   /// imaginary parts.
-   virtual const Vector & GetEigenvector(unsigned int i) const = 0;
-
-   /// @brief Transfer ownership of the converged eigenvectors
-   ///
-   /// The array should contain 2*numConverged vectors with the even
-   /// indices corresponding to the real parts of the converged
-   /// eigenvectors and the odd indices corresponding to the imaginary
-   /// parts.
-   virtual Vector ** StealEigenvectors() = 0;
-};
-
-/// Complex Generalized Eigensolver
-/// Can have arbitrary complex-valued eigenvalues and eigenvectors
-class ComplexGenEigensolver :
-   public EigensolverBase, public ComplexGenEigenequation
-{
-protected:
-   ComplexGenEigensolver() = default;
-
-public:
-   virtual ~ComplexGenEigensolver() = default;
-
-   /// @brief Collect the converged eigenvalues
-   ///
-   /// The length of the array should be twice the number of converged
-   /// eigenvalues. The complex-valued eigenvalues can be constructed
-   /// as: lambda_j = eig[2*j]+i*eig[2*j+1]
-   virtual void GetEigenvalues(Array<real_t> & eig) const = 0;
-
-   /// @brief Extract a single eigenvector
-   ///
-   /// The index i should be in the range [0, 2*numConverged). The
-   /// vectors corresponding to even indices are the real parts of the
-   /// converged eigenvectors and the odd indices correspond to the
-   /// imaginary parts.
-   virtual const Vector & GetEigenvector(unsigned int i) const = 0;
-
-   /// @brief Transfer ownership of the converged eigenvectors
-   ///
-   /// The array should contain 2*numConverged vectors with the even
-   /// indices corresponding to the real parts of the converged
-   /// eigenvectors and the odd indices corresponding to the imaginary
-   /// parts.
-   virtual Vector ** StealEigenvectors() = 0;
-};
-
-}
-
-#endif

linalg/hypre.cpp

@@ -6556,7 +6556,7 @@ HypreLOBPCG::SetPreconditioner(Solver & precond)
 }
 
 void
-HypreLOBPCG::SetOperator(const Operator & A)
+HypreLOBPCG::SetOperator(Operator & A)
 {
    HYPRE_BigInt locSize = A.Width();
 
@@ -6603,7 +6603,7 @@ HypreLOBPCG::SetOperator(const Operator & A)
 }
 
 void
-HypreLOBPCG::SetMassMatrix(const Operator & M)
+HypreLOBPCG::SetMassMatrix(Operator & M)
 {
    matvec_fn.MatvecCreate  = this->OperatorMatvecCreate;
    matvec_fn.Matvec        = this->OperatorMatvec;
@@ -6624,7 +6624,7 @@ HypreLOBPCG::GetEigenvalues(Array<real_t> & eigs) const
    }
 }
 
-const Vector &
+const HypreParVector &
 HypreLOBPCG::GetEigenvector(unsigned int i) const
 {
    return multi_vec->GetVector(i);
@@ -6866,24 +6866,6 @@ HypreAME::SetPreconditioner(HypreSolver & precond)
    ams_precond = &precond;
 }
 
-void
-HypreAME::SetOperators(const Operator & opA, const Operator & opB)
-{
-   const HypreParMatrix * A = dynamic_cast<const HypreParMatrix *>(&opA);
-   if (A == NULL)
-   {
-      mfem_error("HypreAME::SetOperator : first operator not HypreParMatrix!");
-   }
-   SetOperator(*A);
-
-   const HypreParMatrix * B = dynamic_cast<const HypreParMatrix *>(&opB);
-   if (B == NULL)
-   {
-      mfem_error("HypreAME::SetOperator : second operator not HypreParMatrix!");
-   }
-   SetMassMatrix(*B);
-}
-
 void
 HypreAME::SetOperator(const HypreParMatrix & A)
 {
@@ -6942,7 +6924,7 @@ HypreAME::createDummyVectors() const
    }
 }
 
-const Vector &
+const HypreParVector &
 HypreAME::GetEigenvector(unsigned int i) const
 {
    if ( eigenvectors == NULL )
@@ -6953,7 +6935,7 @@ HypreAME::GetEigenvector(unsigned int i) const
    return *eigenvectors[i];
 }
 
-Vector **
+HypreParVector **
 HypreAME::StealEigenvectors()
 {
    if ( eigenvectors == NULL )
@@ -6962,11 +6944,11 @@ HypreAME::StealEigenvectors()
    }
 
    // Set the local pointers to NULL so that they won't be deleted later
-   Vector ** vecs = (Vector**)eigenvectors;
+   HypreParVector ** vecs = eigenvectors;
    eigenvectors = NULL;
    multi_vec = NULL;
 
-   return (Vector**)vecs;
+   return vecs;
 }
 
 }

linalg/hypre.hpp

@@ -18,9 +18,7 @@
 
 #include "../general/globals.hpp"
 #include "sparsemat.hpp"
-#include "eigensolvers.hpp"
 #include "hypre_parcsr.hpp"
-#include "eigensolvers.hpp"
 #include <mpi.h>
 
 // Enable internal hypre timing routines
@@ -2148,7 +2146,7 @@ public:
     A. Knyazev, M. Argentati, I. Lashuk, and E. Ovtchinnikov, SISC, 29(5),
     2224-2239, 2007.
 */
-class HypreLOBPCG : public SymGenEigensolver
+class HypreLOBPCG
 {
 private:
    MPI_Comm comm;
@@ -2238,43 +2236,38 @@ public:
    HypreLOBPCG(MPI_Comm comm);
    ~HypreLOBPCG();
 
-   void SetTol(real_t tol) override;
+   void SetTol(real_t tol);
    // not implemented in HYPRE
    // real_t GetTol() const;
    void SetRelTol(real_t rel_tol);
    // not implemented in HYPRE
    // real_t GetRelTol() const;
-   void SetMaxIter(int max_iter) override;
+   void SetMaxIter(int max_iter);
    // not implemented in HYPRE
    // int GetMaxIter() const;
-   void SetPrintLevel(int logging) override;
-   void SetNumModes(int num_eigs) override { nev = num_eigs; }
+   void SetPrintLevel(int logging);
+   void SetNumModes(int num_eigs) { nev = num_eigs; }
    void SetPrecondUsageMode(int pcg_mode);
    void SetRandomSeed(int s) { seed = s; }
    void SetInitialVectors(int num_vecs, HypreParVector ** vecs);
 
    // The following four methods support general operators
    void SetPreconditioner(Solver & precond);
-   void SetOperators(const Operator & A, const Operator & B) override
-   { SetOperator(A); SetMassMatrix(B); }
-   void SetOperator(const Operator & A);
-   void SetMassMatrix(const Operator & M);
+   void SetOperator(Operator & A);
+   void SetMassMatrix(Operator & M);
    void SetSubSpaceProjector(Operator & proj) { subSpaceProj = &proj; }
 
    /// Solve the eigenproblem
-   void Solve() override;
-
-   int GetNumConverged() const override { return nev; }
+   void Solve();
 
    /// Collect the converged eigenvalues
-   void GetEigenvalues(Array<real_t> & eigenvalues) const override;
+   void GetEigenvalues(Array<real_t> & eigenvalues) const;
 
    /// Extract a single eigenvector
-   const Vector & GetEigenvector(unsigned int i) const override;
+   const HypreParVector & GetEigenvector(unsigned int i) const;
 
    /// Transfer ownership of the converged eigenvectors
-   Vector ** StealEigenvectors() override
-   { return (Vector**)multi_vec->StealVectors(); }
+   HypreParVector ** StealEigenvectors() { return multi_vec->StealVectors(); }
 };
 
 /** AME eigenvalue solver in hypre
@@ -2299,7 +2292,7 @@ public:
     mass matrix but it seems unlikely that this would be useful so it is not the
     default behavior.
 */
-class HypreAME : public SymGenEigensolver
+class HypreAME
 {
 private:
    int myid;
@@ -2328,31 +2321,28 @@ public:
    HypreAME(MPI_Comm comm);
    ~HypreAME();
 
-   void SetTol(real_t tol) override;
+   void SetTol(real_t tol);
    void SetRelTol(real_t rel_tol);
-   void SetMaxIter(int max_iter) override;
-   void SetPrintLevel(int logging) override;
-   void SetNumModes(int num_eigs) override;
+   void SetMaxIter(int max_iter);
+   void SetPrintLevel(int logging);
+   void SetNumModes(int num_eigs);
 
    // The following four methods support operators of type HypreParMatrix.
    void SetPreconditioner(HypreSolver & precond);
-   void SetOperators(const Operator & opA, const Operator & opB) override;
    void SetOperator(const HypreParMatrix & A);
    void SetMassMatrix(const HypreParMatrix & M);
 
    /// Solve the eigenproblem
-   void Solve() override;
-
-   int GetNumConverged() const override { return nev; }
+   void Solve();
 
    /// Collect the converged eigenvalues
-   void GetEigenvalues(Array<real_t> & eigenvalues) const override;
+   void GetEigenvalues(Array<real_t> & eigenvalues) const;
 
    /// Extract a single eigenvector
-   const Vector & GetEigenvector(unsigned int i) const override;
+   const HypreParVector & GetEigenvector(unsigned int i) const;
 
    /// Transfer ownership of the converged eigenvectors
-   Vector ** StealEigenvectors() override;
+   HypreParVector ** StealEigenvectors();
 };
 
 }

linalg/linalg.hpp

@@ -28,7 +28,6 @@
 #include "symmat.hpp"
 #include "ode.hpp"
 #include "solvers.hpp"
-#include "eigensolvers.hpp"
 #include "handle.hpp"
 #include "invariants.hpp"
 #include "constraints.hpp"
@@ -58,10 +57,6 @@
 #include "ginkgo.hpp"
 #endif
 
-#ifdef MFEM_USE_ARPACK
-#include "arpack.hpp"
-#endif
-
 #ifdef MFEM_USE_MKL_PARDISO
 #include "pardiso.hpp"
 #endif

linalg/operator.hpp

@@ -844,22 +844,6 @@ public:
 };
 
 
-/// Zero Operator N: x -> 0.
-class ZeroOperator : public Operator
-{
-public:
-   /// Create an zero operator of size @a n.
-   explicit ZeroOperator(int n) : Operator(n) { }
-
-   /// Operator application
-   void Mult(const Vector &x, Vector &y) const override
-   { y.SetSize(width); y = 0_r; }
-
-   /// Application of the transpose
-   void MultTranspose(const Vector &x, Vector &y) const override
-   { y.SetSize(width); y = 0_r; }
-};
-
 /// Identity Operator I: x -> x.
 class IdentityOperator : public Operator
 {

linalg/petsc.cpp

@@ -4156,31 +4156,20 @@ void PetscNonlinearSolver::SetUpdate(void (*update)(Operator *,int,
 void PetscNonlinearSolver::Mult(const Vector &b, Vector &x) const
 {
    SNES snes = (SNES)obj;
-   MPI_Comm comm = PetscObjectComm(obj);
 
-   // Reduction needed: some processes may have null local size while others don't,
-   // and VecPlaceArray (used by PlaceMemory) is a logically collective operation.
-   PetscBool b_nonempty = b.Size() ? PETSC_TRUE : PETSC_FALSE;
-#if PETSC_VERSION_LT(3,24,0)
-   mpiierr = MPI_Allreduce(MPI_IN_PLACE,&b_nonempty,1,MPIU_BOOL,MPI_LOR,comm);
-#else
-   mpiierr = MPI_Allreduce(MPI_IN_PLACE,&b_nonempty,1,MPI_C_BOOL,MPI_LOR,comm);
-#endif
-   CCHKERRQ(comm,mpiierr);
-
-   // Always create B with allocate=false so that PlaceMemory can be called on
-   // it regardless of whether b was empty on a previous call.
-   if (!B) { B = new PetscParVector(comm, *this, true, false); }
-   if (!X) { X = new PetscParVector(comm, *this, false, false); }
+   bool b_nonempty = b.Size();
+   if (!B) { B = new PetscParVector(PetscObjectComm(obj), *this, true); }
+   if (!X) { X = new PetscParVector(PetscObjectComm(obj), *this, false, false); }
    X->PlaceMemory(x.GetMemory(),iterative_mode);
    if (b_nonempty) { B->PlaceMemory(b.GetMemory()); }
+   else { *B = 0.0; }
 
    Customize();
 
    if (!iterative_mode) { *X = 0.; }
 
-   // Solve the system. Pass nullptr for b when empty (PETSc treats it as zero RHS).
-   ierr = SNESSolve(snes, b_nonempty ? B->x : nullptr, X->x); PCHKERRQ(snes, ierr);
+   // Solve the system.
+   ierr = SNESSolve(snes, B->x, X->x); PCHKERRQ(snes, ierr);
    X->ResetMemory();
    if (b_nonempty) { B->ResetMemory(); }
 }

makefile

@@ -126,11 +126,11 @@ EXAMPLE_TEST_DIRS := examples
 MINIAPP_SUBDIRS = common electromagnetics meshing performance tools \
  toys nurbs gslib adjoint solvers shifted mtop parelag tribol autodiff dfem \
  hooke multidomain dpg hdiv-linear-solver spde diag-smoothers contact \
- fluids/navier fluids/schrodinger-flow plasma plasma/pic
+ fluids/navier fluids/schrodinger-flow plasma
 MINIAPP_DIRS := $(addprefix miniapps/,$(MINIAPP_SUBDIRS))
 MINIAPP_TEST_DIRS := $(filter-out %/common,$(MINIAPP_DIRS))
 MINIAPP_USE_COMMON := $(addprefix miniapps/,electromagnetics meshing tools \
- toys shifted dpg diag-smoothers fluids/navier plasma plasma/pic)
+ toys shifted dpg diag-smoothers fluids/navier plasma)
 
 EM_DIRS = $(EXAMPLE_DIRS) $(MINIAPP_DIRS)
 
@@ -302,7 +302,7 @@ endif
 MFEM_REQ_LIB_DEPS = SUPERLU MUMPS METIS FMS CONDUIT SIDRE LAPACK SUNDIALS\
  SUITESPARSE STRUMPACK GINKGO GNUTLS HDF5 NETCDF SLEPC PETSC MPFR PUMI HIOP\
  GSLIB OCCA CEED RAJA UMPIRE MKL_CPARDISO MKL_PARDISO AMGX MAGMA CALIPER PARELAG\
- TRIBOL BENCHMARK MOONOLITH ALGOIM ARPACK
+ TRIBOL BENCHMARK MOONOLITH ALGOIM
 
 
 PETSC_ERROR_MSG = $(if $(PETSC_FOUND),,. PETSC config not found: $(PETSC_VARS))
@@ -371,8 +371,7 @@ MFEM_DEFINES = MFEM_VERSION MFEM_VERSION_STRING MFEM_GIT_STRING MFEM_USE_MPI\
  MFEM_USE_SIMD MFEM_USE_ADIOS2 MFEM_USE_MKL_CPARDISO MFEM_USE_MKL_PARDISO MFEM_USE_AMGX\
  MFEM_USE_MAGMA MFEM_USE_MUMPS MFEM_USE_ADFORWARD MFEM_USE_CODIPACK MFEM_USE_CALIPER\
  MFEM_USE_BENCHMARK MFEM_USE_PARELAG MFEM_USE_TRIBOL MFEM_USE_ALGOIM MFEM_USE_ENZYME\
- MFEM_SOURCE_DIR MFEM_INSTALL_DIR MFEM_SHARED_BUILD MFEM_USE_DOUBLE MFEM_USE_SINGLE\
-MFEM_USE_ARPACK
+ MFEM_SOURCE_DIR MFEM_INSTALL_DIR MFEM_SHARED_BUILD MFEM_USE_DOUBLE MFEM_USE_SINGLE
 
 # List of makefile variables that will be written to config.mk:
 MFEM_CONFIG_VARS = MFEM_CXX MFEM_HOST_CXX MFEM_CPPFLAGS MFEM_CXXFLAGS\
@@ -734,7 +733,6 @@ status info:
 	$(info MFEM_TIMER_TYPE        = $(MFEM_TIMER_TYPE))
 	$(info MFEM_USE_SUNDIALS      = $(MFEM_USE_SUNDIALS))
 	$(info MFEM_USE_SUITESPARSE   = $(MFEM_USE_SUITESPARSE))
-	$(info MFEM_USE_ARPACK        = $(MFEM_USE_ARPACK))
 	$(info MFEM_USE_SUPERLU       = $(MFEM_USE_SUPERLU))
 	$(info MFEM_USE_SUPERLU5      = $(MFEM_USE_SUPERLU5))
 	$(info MFEM_USE_MUMPS         = $(MFEM_USE_MUMPS))

mesh/mesh.cpp

@@ -1616,9 +1616,7 @@ Element::Type Mesh::GetFaceElementType(int Face) const
 
 Array<int> Mesh::GetFaceToBdrElMap() const
 {
-   Array<int> face_to_be(Dim == 1 ? NumOfVertices :
-                         Dim == 2 ? NumOfEdges :
-                         Dim == 3 ? NumOfFaces : 0);
+   Array<int> face_to_be(Dim == 2 ? NumOfEdges : NumOfFaces);
    face_to_be = -1;
    for (int i = 0; i < NumOfBdrElements; i++)
    {

mesh/mesh_operators.cpp

@@ -63,6 +63,7 @@ ThresholdRefiner::ThresholdRefiner(ErrorEstimator &est)
 
    threshold = 0.0;
    num_marked_elements = 0LL;
+   current_sequence = -1;
 
    non_conforming = -1;
    nc_limit = 0;
@@ -86,6 +87,7 @@ int ThresholdRefiner::MarkWithoutRefining(Mesh & mesh,
    threshold = 0.0;
    num_marked_elements = 0LL;
    refinements.SetSize(0);
+   current_sequence = mesh.GetSequence();
 
    const long long num_elements = mesh.GetGlobalNE();
    if (num_elements >= max_elements) { return STOP; }
@@ -147,6 +149,7 @@ int ThresholdRefiner::ApplyImpl(Mesh &mesh)
 void ThresholdRefiner::Reset()
 {
    estimator.Reset();
+   current_sequence = -1;
    num_marked_elements = 0LL;
    // marked_elements.SetSize(0); // not necessary
 }

mesh/mesh_operators.hpp

@@ -188,6 +188,7 @@ protected:
    long long num_marked_elements;
 
    Array<Refinement> marked_elements;
+   long current_sequence;
 
    int non_conforming;
    int nc_limit;

mesh/pmesh.cpp

@@ -5639,12 +5639,6 @@ Mesh ParMesh::GetSerialMesh(int save_rank) const
       }
    }
 
-   if (MyRank == save_rank)
-   {
-      attribute_sets.Copy(serialmesh.attribute_sets);
-      bdr_attribute_sets.Copy(serialmesh.bdr_attribute_sets);
-   }
-
    MPI_Barrier(MyComm);
    return serialmesh;
 }

mesh/submesh/psubmesh.hpp

@@ -227,29 +227,15 @@ public:
                                            const ParGridFunction &dst);
 
    /**
-   * @brief Check if Mesh @a m is a ParSubMesh.
+   * @brief Check if ParMesh @a m is a ParSubMesh.
    *
-   * @param m The input Mesh
+   * @param m The input ParMesh
    */
-   static bool IsParSubMesh(const Mesh *m)
+   static bool IsParSubMesh(const ParMesh *m)
    {
       return dynamic_cast<const ParSubMesh *>(m) != nullptr;
    }
 
-   /**
-   * @brief Check if Mesh @a sub is a ParSubMesh of Mesh @a parent.
-   *
-   * @param sub The potential submesh Mesh
-   * @param parent The potential parent Mesh
-   */
-   static bool IsParSubMesh(const Mesh* sub, const Mesh* parent)
-   {
-      while (IsParSubMesh(sub) &&
-             (sub = static_cast<const ParSubMesh *>(sub)->GetParent()) &&
-             sub != parent);
-      return sub == parent;
-   }
-
 private:
    ParSubMesh(const ParMesh &parent, SubMesh::From from,
               const Array<int> &attributes);

mesh/submesh/submesh.hpp

@@ -225,20 +225,6 @@ public:
       return dynamic_cast<const SubMesh *>(m) != nullptr;
    }
 
-   /**
-   * @brief Check if Mesh @a sub is a SubMesh of Mesh @a parent.
-   *
-   * @param sub The potential submesh Mesh
-   * @param parent The potential parent Mesh
-   */
-   static bool IsSubMesh(const Mesh* sub, const Mesh* parent)
-   {
-      while (IsSubMesh(sub) &&
-             (sub = static_cast<const SubMesh *>(sub)->GetParent()) &&
-             sub != parent);
-      return sub == parent;
-   }
-
 private:
    /// Private constructor
    SubMesh(const Mesh &parent, From from, const Array<int> &attributes);

miniapps/electromagnetics/CMakeLists.txt

@@ -43,39 +43,19 @@ endif()
   
   # Add the corresponding tests to the "test" target
   if (MFEM_ENABLE_TESTING)
-    add_test(NAME tesla_1_np=${MFEM_MPI_NP}
+    add_test(NAME tesla_np=4
       COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
       ${MPIEXEC_PREFLAGS}
       $<TARGET_FILE:tesla> -no-vis -maxit 2 -cr "0 0 -0.2 0 0 0.2 0.2 0.4 1"
       ${MPIEXEC_POSTFLAGS})
-    add_test(NAME tesla_2_np=${MFEM_MPI_NP}
+
+    add_test(NAME volta_np=4
       COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
       ${MPIEXEC_PREFLAGS}
-      $<TARGET_FILE:tesla>
-      -no-vis -maxit 2 -m ../../data/inline-hex.mesh -ubbc "0 0 1"
+      $<TARGET_FILE:volta> -no-vis -maxit 2 -dbcs 1 -dbcg -ds "0.0 0.0 0.0 0.2 8.0"
       ${MPIEXEC_POSTFLAGS})
 
-    add_test(NAME volta_1_np=${MFEM_MPI_NP}
-      COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
-      ${MPIEXEC_PREFLAGS}
-      $<TARGET_FILE:volta>
-      -no-vis -maxit 2 -dbcs 1 -dbcg -ds "0.0 0.0 0.0 0.2 8.0"
-      ${MPIEXEC_POSTFLAGS})
-    add_test(NAME volta_2_np=${MFEM_MPI_NP}
-      COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
-      ${MPIEXEC_PREFLAGS}
-      $<TARGET_FILE:volta>
-      -no-vis -maxit 2 -m ../../data/square-disc.mesh -dbcs "1 2 3 4 5 6 7 8"
-      -dbcv "0 0 0 0 1 1 1 1"
-      ${MPIEXEC_POSTFLAGS})
-    add_test(NAME volta_3_np=${MFEM_MPI_NP}
-      COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
-      ${MPIEXEC_PREFLAGS}
-      $<TARGET_FILE:volta>
-      -no-vis -maxit 2 -m ../../data/inline-hex.mesh -dbcs "1 6" -dbcv "0 1"
-      ${MPIEXEC_POSTFLAGS})
-
-    add_test(NAME joule_np=${MFEM_MPI_NP}
+    add_test(NAME joule_np=4
       COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
       ${MPIEXEC_PREFLAGS}
       $<TARGET_FILE:joule>
@@ -83,41 +63,12 @@ endif()
       ${MPIEXEC_POSTFLAGS})
 
     if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
-      add_test(NAME maxwell_np=${MFEM_MPI_NP}
+      add_test(NAME maxwell_np=4
         COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
         ${MPIEXEC_PREFLAGS}
         $<TARGET_FILE:maxwell>
         -no-vis -abcs "-1" -dp "-0.3 0.0 0.0 0.3 0.0 0.0 0.1 1 .5 .5"
         ${MPIEXEC_POSTFLAGS})
     endif()
-
-    if (MFEM_USE_GSLIB)
-      add_test(NAME lorentz_1_np=${MFEM_MPI_NP}
-        COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
-        ${MPIEXEC_PREFLAGS}
-        $<TARGET_FILE:lorentz>
-        -no-vis -er Volta-AMR-Parallel -ec 2 -npt 100 -xmin "0.0 0.0 0.0"
-        -xmax "1.0 1.0 1.0" -pmin "1 0 0" -pmax "1 0 0" -rdf 0 -vt 0 -nt 100
-        ${MPIEXEC_POSTFLAGS})
-      # Setup dependency on volta_3_np=<np>
-      set_tests_properties(volta_3_np=${MFEM_MPI_NP}
-        PROPERTIES FIXTURES_SETUP Volta3)
-      set_tests_properties(lorentz_1_np=${MFEM_MPI_NP}
-        PROPERTIES FIXTURES_REQUIRED Volta3)
-
-      add_test(NAME lorentz_2_np=${MFEM_MPI_NP}
-        COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
-        ${MPIEXEC_PREFLAGS}
-        $<TARGET_FILE:lorentz>
-        -no-vis -br Tesla-AMR-Parallel -bc 2 -npt 10 -xmin "0.0 0.0 0.0"
-        -xmax "1.0 1.0 1.0" -pmin "0 0.1 0.05" -pmax "0 0.4 0.1" -nt 1000 -rdf 0
-        -vt 0
-        ${MPIEXEC_POSTFLAGS})
-      # Setup dependency on tesla_2_np=<np>
-      set_tests_properties(tesla_2_np=${MFEM_MPI_NP}
-        PROPERTIES FIXTURES_SETUP Tesla2)
-      set_tests_properties(lorentz_2_np=${MFEM_MPI_NP}
-        PROPERTIES FIXTURES_REQUIRED Tesla2)
-    endif()
   endif()
 endif()

miniapps/electromagnetics/makefile

@@ -117,10 +117,10 @@ joule-test-par: joule
 lorentz-test-par: lorentz-test-1 lorentz-test-2
 lorentz-test-1: lorentz volta-test-3
 	@$(call mfem-test,$<, $(RUN_MPI), Electromagnetic miniapp,\
-	-er Volta-AMR-Parallel -ec 2 -npt 100 -xmin '0.0 0.0 0.0' -xmax '1.0 1.0 1.0' -pmin '1 0 0' -pmax '1 0 0' -rdf 0 -vt 0 -nt 100)
+	-er Volta-AMR-Parallel -ec 2 -npt 100 -xmin '0.0 0.0 0.0' -xmax '1.0 1.0 1.0' -pmin '1 0 0' -pmax '1 0 0' -rdf 0 -vt 0 -nt 100')
 lorentz-test-2: lorentz tesla-test-2
 	@$(call mfem-test,$<, $(RUN_MPI), Electromagnetic miniapp,\
-	-br Tesla-AMR-Parallel -bc 2 -npt 10 -xmin '0.0 0.0 0.0' -xmax '1.0 1.0 1.0' -pmin '0 0.1 0.05' -pmax '0 0.4 0.1' -nt 1000 -rdf 0 -vt 0)
+	-br Tesla-AMR-Parallel -bc 2 -br Tesla-AMR-Parallel -npt 10 -xmin '0.0 0.0 0.0' -xmax '1.0 1.0 1.0' -pmin '0 0.1 0.05' -pmax '0 0.4 0.1' -nt 1000 -rdf 0 -vt 0)
 
 # Testing: "test" target and mfem-test* variables are defined in config/test.mk
 

miniapps/hdiv-linear-solver/change_basis.cpp

@@ -22,7 +22,7 @@ void ComputeInverse(const Array<real_t> &A, Array<real_t> &Ainv)
 {
    Array<real_t> A2 = A;
    const int n2 = A.Size();
-   const int n = static_cast<int>(sqrt(n2));
+   const int n = static_cast<const int>(sqrt(n2));
    Array<int> ipiv(n);
    LUFactors lu(A2.GetData(), ipiv.GetData());
    lu.Factor(n);
@@ -58,7 +58,7 @@ void SubcellIntegrals(int n, const Poly_1D::Basis &basis, Array<real_t> &B)
 
 void Transpose(const Array<real_t> &B, Array<real_t> &Bt)
 {
-   const int n = static_cast<int>(sqrt(B.Size()));
+   const int n = static_cast<const int>(sqrt(B.Size()));
    Bt.SetSize(n*n);
    for (int i=0; i<n; ++i) for (int j=0; j<n; ++j) { Bt[i+j*n] = B[j+i*n]; }
 }

miniapps/nurbs/nurbs_ex11p.cpp

@@ -329,8 +329,8 @@ int main(int argc, char *argv[])
 
       for (int i=0; i<nev; i++)
       {
-         // convert eigenvector from Vector to ParGridFunction
-         x.Distribute(lobpcg->GetEigenvector(i));
+         // convert eigenvector from HypreParVector to ParGridFunction
+         x = lobpcg->GetEigenvector(i);
 
          mode_name << "mode_" << setfill('0') << setw(2) << i << "."
                    << setfill('0') << setw(6) << myid;
@@ -357,8 +357,8 @@ int main(int argc, char *argv[])
                  << ", Lambda = " << eigenvalues[i] << endl;
          }
 
-         // convert eigenvector from Vector to ParGridFunction
-         x.Distribute(lobpcg->GetEigenvector(i));
+         // convert eigenvector from HypreParVector to ParGridFunction
+         x = lobpcg->GetEigenvector(i);
 
          mode_sock << "parallel " << num_procs << " " << myid << "\n"
                    << "solution\n" << *pmesh << x << flush

miniapps/plasma/CMakeLists.txt

@@ -23,5 +23,3 @@ if (MFEM_USE_MPI)
       EXTRA_HEADERS ${PLASMA_COMMON_HEADERS})
 
 endif()
-
-add_subdirectory(pic)

miniapps/plasma/makefile

@@ -14,6 +14,9 @@ MFEM_DIR ?= ../..
 MFEM_BUILD_DIR ?= ../..
 SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/plasma/,)
 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)
@@ -26,8 +29,6 @@ else
    MINIAPPS = $(PAR_MINIAPPS) $(SEQ_MINIAPPS)
 endif
 
-PLASMA_SUBDIRS = pic
-
 .SUFFIXES:
 .SUFFIXES: .o .cpp .mk
 .PHONY: all lib-common clean clean-build clean-exec
@@ -46,12 +47,7 @@ COMMON_O=
 %: %.cpp
 %.o: %.cpp
 
-all: $(MINIAPPS) subdirs
-
-.PHONY: subdirs $(PLASMA_SUBDIRS)
-subdirs: $(PLASMA_SUBDIRS)
-$(PLASMA_SUBDIRS): lib-common
-	$(MAKE) -C $(BLD)$(@)
+all: $(MINIAPPS)
 
 # Rules for building the miniapps
 %: $(SRC)%.cpp $(COMMON_O) $(MFEM_LIB_FILE) $(CONFIG_MK) | lib-common
@@ -79,15 +75,11 @@ RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
 $(MFEM_LIB_FILE):
 	$(error The MFEM library is not built)
 
-ALL_CLEAN_SUBDIRS = $(addsuffix /clean,$(PLASMA_SUBDIRS))
-.PHONY: $(ALL_CLEAN_SUBDIRS)
-$(ALL_CLEAN_SUBDIRS):
-	$(MAKE) -C $(BLD)$(@D) $(@F)
-
 clean: clean-build clean-exec
 
-clean-build: $(addsuffix /clean,$(PLASMA_SUBDIRS))
+clean-build:
 	rm -f *.o *~ $(SEQ_MINIAPPS) $(PAR_MINIAPPS)
 	rm -rf *.dSYM *.TVD.*breakpoints
 
 clean-exec:
+

miniapps/plasma/pic/CMakeLists.txt

@@ -1,28 +0,0 @@
-# Copyright (c) 2010-2025, 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.
-
-if (MFEM_USE_MPI AND MFEM_USE_GSLIB)
-  add_mfem_miniapp(electrostatic-pic
-    MAIN electrostatic-pic.cpp
-    EXTRA_HEADERS ${MFEM_MINIAPPS_COMMON_HEADERS}
-    LIBRARIES mfem-common)
-
-  # Add the corresponding tests to the "test" target
-  if (MFEM_ENABLE_TESTING)
-    add_test(NAME electrostatic-pic_np=${MFEM_MPI_NP}
-      COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
-      ${MPIEXEC_PREFLAGS}
-      $<TARGET_FILE:electrostatic-pic> -rdi 2 -npt 40960 -k 0.2855993321 -a 0.05
-        -nt 200 -nx 16 -ny 16 -O 1 -q 0.01181640625 -m 0.01181640625 -oci 1000
-        -dt 0.1
-      ${MPIEXEC_POSTFLAGS})
-  endif()
-endif()

miniapps/plasma/pic/electrostatic-pic.cpp

@@ -1,788 +0,0 @@
-// Copyright (c) 2010-2025, 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.
-//
-//           -----------------------------------------------------
-//           Particle-In-Cell (PIC) Simulation (2D/3D)
-//           -----------------------------------------------------
-//
-// This miniapp performs a Particle-In-Cell simulation (supports 2D or 3D
-// spatial dimensions) of multiple charged particles subject to electric
-// field forces.
-//
-//                           dp/dt = q E
-//
-// The method used is explicit time integration with a leap-frog scheme.
-//
-// The electric field is computed from the particle charge distribution using
-// a Poisson solver. The particle trajectories are computed within a periodic
-// domain (2D or 3D).
-//
-// Solution process (per timestep, repeating steps 1-6):
-//   (1) Deposit charge from particles to grid via Dirac delta function
-//       to form the RHS of the Poisson equation
-//   (2) Solve Poisson equation (-Δφ = ρ - ρ_0) to compute potential φ, where
-//       ρ_0 is a constant neutralizing term that enforces global charge
-//       neutrality.
-//   (3) Compute electric field E = -∇φ from the potential
-//   (4) Interpolate E-field to particle positions
-//   (5) Push particles using leap-frog scheme (update momentum and position)
-//   (6) Redistribute particles across processors
-//
-// Compile with: make electrostatic-pic
-//
-// Sample runs:
-//
-//   2D2V Linear Landau damping test case (Ricketson & Hu, 2025):
-//      mpirun -n 4 ./electrostatic-pic -rdi 1 -npt 409600 -k 0.2855993321 -a 0.05 -nt 200 -nx 32 -ny 32 -O 1 -q 0.001181640625 -m 0.001181640625 -oci 1000 -dt 0.1
-//   3D3V Linear Landau damping test case (Zheng et al., 2025):
-//    * mpirun -n 128 ./electrostatic-pic -dim 3 -rdi 1 -npt 40960000 -k 0.5 -a 0.01 -nt 100 -nx 32 -ny 32 -nz 32 -O 1 -q 0.00004844730731 -m 0.00004844730731 -oci 1000 -dt 0.02 -no-vis
-
-#include "mfem.hpp"
-#include "../../../general/text.hpp"
-#include "../../common/fem_extras.hpp"
-#include "../../common/particles_extras.hpp"
-#include "../../common/pfem_extras.hpp"
-
-#include <ctime>
-#include <fstream>
-#include <iomanip>
-#include <iostream>
-#include <random>
-#include <string>
-#include <vector>
-
-#define EPSILON 1  // ε_0
-
-using namespace std;
-using namespace mfem;
-using namespace mfem::common;
-
-struct PICContext
-{
-   int dim = 2;     ///< Spatial dimension.
-   int order = 1;   ///< FE order for spatial discretization.
-   int nx = 100;    ///< Number of grid cells in x-direction.
-   int ny = 100;    ///< Number of grid cells in y-direction.
-   int nz = 100;    ///< Number of grid cells in z-direction.
-   real_t L = 1.0;  ///< Domain length.
-
-   int ordering = 1;  ///< Ordering of particles.
-   int npt = 1000;    ///< Number of particles.
-   real_t q = 1.0;    ///< Particle charge.
-   real_t m = 1.0;    ///< Particle mass.
-
-   real_t k = 1.0;      ///< Wave number (Landau damping init).
-   real_t alpha = 0.1;  ///< Perturbation amplitude (Landau damping init).
-
-   real_t dt = 1e-2;  ///< Time step size.
-
-   int nt = 1000;            ///< Number of time steps to run.
-   int redist_interval = 5;  ///< Redistribution and update E_gf interval.
-   int output_csv_interval = 1000;  ///< Interval for outputting CSV data files.
-
-   bool visualization = true;  ///< Enable visualization.
-   int visport = 19916;        ///< Port number for visualization server.
-   bool reproduce = true;      ///< Enable reproducible results.
-} ctx;
-
-/** This class implements explicit time integration for charged particles
-    in an electric field using ParticleSet. */
-class ParticleMover
-{
-public:
-   enum Fields
-   {
-      MASS,    // vdim = 1
-      CHARGE,  // vdim = 1
-      MOM,     // vdim = dim
-      EFIELD   // vdim = dim
-   };
-
-protected:
-   /// Pointers to E field GridFunctions
-   ParGridFunction* E_gf;
-
-   /// FindPointsGSLIB object for E field mesh
-   FindPointsGSLIB& E_finder;
-
-   /// ParticleSet of charged particles
-   std::unique_ptr<ParticleSet> charged_particles;
-
-   /// Temporary vectors for particle computation
-   mutable Vector pm_, pp_;
-
-public:
-   ParticleMover(MPI_Comm comm, ParGridFunction* E_gf_,
-                 FindPointsGSLIB& E_finder_, int num_particles,
-                 Ordering::Type pdata_ordering);
-
-   /// Initialize charged particles with given parameters
-   void InitializeChargedParticles(const real_t& k, const real_t& alpha,
-                                   real_t m, real_t q, real_t L,
-                                   bool reproduce = false);
-
-   /// Find Particles in mesh corresponding to E and field
-   void FindParticles();
-
-   /// Advance particles one time step using Boris algorithm
-   void Step(real_t& t, real_t dt, real_t L, bool first_step = false);
-
-   /// Redistribute particles across processors
-   void Redistribute();
-
-   /// Get reference to ParticleSet
-   ParticleSet& GetParticles() { return *charged_particles; }
-
-   /// Compute (global) kinetic energy from particles
-   /** Optionally, advance the particle momenta by time step @a dt. */
-   real_t ComputeKineticEnergy(real_t dt = 0.) const;
-};
-
-/** Field solver responsible for updating the electrostatic potential and field
-    from the particle charge density. Assembles and solves the periodic Poisson
-    problem, computes the electric field via a discrete gradient operator, and
-    provides utilities for field diagnostics (e.g. global field energy). */
-class FieldSolver
-{
-private:
-   real_t domain_volume;
-   real_t neutralizing_const;
-   ParLinearForm* precomputed_neutralizing_lf = nullptr;
-   bool precompute_neutralizing_const = false;
-   // Diffusion matrix
-   HypreParMatrix* diffusion_matrix;
-   // Gradient operator for computing E = -∇φ
-   ParDiscreteLinearOperator* grad_interpolator;
-   FindPointsGSLIB& E_finder;
-   ParLinearForm b;
-
-protected:
-   /** Compute neutralizing constant and initialize with the constant.
-       Returns a reference to the precomputed neutralizing ParLinearForm. */
-   const ParLinearForm& ComputeNeutralizingRHS(ParFiniteElementSpace* pfes,
-                                               const ParticleVector& Q,
-                                               MPI_Comm comm);
-
-   /** Deposit charge from particles into a ParLinearForm (RHS b).
-       b_i = sum_p q_p * φ_i(x_p) */
-   void DepositCharge(ParFiniteElementSpace* pfes, const ParticleVector& Q);
-
-public:
-   FieldSolver(ParFiniteElementSpace* phi_fes, ParFiniteElementSpace* E_fes,
-               FindPointsGSLIB& E_finder_,
-               bool precompute_neutralizing_const_ = false);
-
-   ~FieldSolver();
-
-   /** Update the phi_gf grid function from the particles.
-       Solve periodic Poisson: diffusion_matrix * phi = (rho - <rho>)
-       with zero-mean enforcement via OrthoSolver. */
-   void UpdatePhiGridFunction(ParticleSet& particles, ParGridFunction& phi_gf);
-
-   /** Update E_gf grid function from phi_gf grid function.
-       Compute the gradient: E = -∇φ. */
-   void UpdateEGridFunction(ParGridFunction& phi_gf, ParGridFunction& E_gf);
-
-   /// Compute (global) field energy: 0.5 * ∫ ||E||^2 dx
-   real_t ComputeFieldEnergy(const ParGridFunction& E_gf) const;
-};
-
-/// Prints the program's logo to the given output stream
-void display_banner(ostream& os);
-
-int main(int argc, char* argv[])
-{
-   Mpi::Init(argc, argv);
-   int num_ranks = Mpi::WorldSize();
-   int rank = Mpi::WorldRank();
-   Hypre::Init();
-
-   if (Mpi::Root()) { display_banner(cout); }
-
-   OptionsParser args(argc, argv);
-   args.AddOption(&ctx.dim, "-dim", "--dimension",
-                  "Spatial dimension (2 or 3)");
-   args.AddOption(&ctx.order, "-O", "--order",
-                  "Finite element polynomial degree");
-   args.AddOption(&ctx.nx, "-nx", "--num-x",
-                  "Number of elements in the x direction.");
-   args.AddOption(&ctx.ny, "-ny", "--num-y",
-                  "Number of elements in the y direction.");
-   args.AddOption(&ctx.nz, "-nz", "--num-z",
-                  "Number of elements in the z direction.");
-   args.AddOption(&ctx.q, "-q", "--charge", "Particle charge.");
-   args.AddOption(&ctx.m, "-m", "--mass", "Particle mass.");
-   args.AddOption(&ctx.dt, "-dt", "--time-step", "Time Step.");
-   args.AddOption(&ctx.nt, "-nt", "--num-timesteps", "Number of timesteps.");
-   args.AddOption(&ctx.npt, "-npt", "--num-particles",
-                  "Total number of particles.");
-   args.AddOption(&ctx.k, "-k", "--k", "Wave number for initial distribution.");
-   args.AddOption(&ctx.alpha, "-a", "--alpha",
-                  "Perturbation amplitude for initial distribution.");
-   args.AddOption(&ctx.ordering, "-o", "--ordering",
-                  "Ordering of particle data. 0 = byNODES, 1 = byVDIM.");
-   args.AddOption(&ctx.redist_interval, "-rdi", "--redist-interval",
-                  "Redistribution and update E_gf interval. Disabled if < 0.");
-   args.AddOption(&ctx.output_csv_interval, "-oci", "--output-csv-interval",
-                  "Output CSV interval. Disabled if < 0.");
-   args.AddOption(&ctx.visualization, "-vis", "--visualization", "-no-vis",
-                  "--no-visualization",
-                  "Enable or disable GLVis visualization.");
-   args.AddOption(&ctx.visport, "-p", "--send-port", "Socket for GLVis.");
-   args.AddOption(&ctx.reproduce, "-rep", "--reproduce", "-no-rep",
-                  "--no-reproduce",
-                  "Enable or disable reproducible random seed.");
-   args.Parse();
-   if (!args.Good())
-   {
-      if (Mpi::Root()) { args.PrintUsage(cout); }
-      return 1;
-   }
-   if (Mpi::Root()) { args.PrintOptions(cout); }
-
-   // Assert that dimension is 2 or 3
-   MFEM_VERIFY(ctx.dim == 2 || ctx.dim == 3,
-               "Dimension must be 2 or 3, got " << ctx.dim);
-   MFEM_VERIFY(ctx.alpha >= -1.0 && ctx.alpha < 1.0,
-               "Alpha should be in range [-1, 1).");
-   MFEM_VERIFY(ctx.k > 0.0,
-               "k must be nonzero for displacement initialization.");
-
-   ctx.L = 2.0 * M_PI / ctx.k;
-
-   // 1. make a Cartesian Mesh (2D or 3D)
-   Mesh serial_mesh;
-   std::vector<Vector> translations;
-
-   if (ctx.dim == 2)
-   {
-      serial_mesh = Mesh(Mesh::MakeCartesian2D(
-                            ctx.nx, ctx.ny, Element::QUADRILATERAL, false, ctx.L, ctx.L));
-      translations = {Vector({ctx.L, 0.0}), Vector({0.0, ctx.L})};
-   }
-   else  // ctx.dim == 3
-   {
-      serial_mesh = Mesh(Mesh::MakeCartesian3D(
-                            ctx.nx, ctx.ny, ctx.nz, Element::HEXAHEDRON, ctx.L, ctx.L, ctx.L));
-      translations = {Vector({ctx.L, 0.0, 0.0}), Vector({0.0, ctx.L, 0.0}),
-                      Vector({0.0, 0.0, ctx.L})
-                     };
-   }
-
-   Mesh periodic_mesh(Mesh::MakePeriodic(
-                         serial_mesh, serial_mesh.CreatePeriodicVertexMapping(translations)));
-   // 2. Partition and distribute the mesh
-   ParMesh mesh(MPI_COMM_WORLD, periodic_mesh);
-   serial_mesh.Clear();    // the serial mesh is no longer needed
-   periodic_mesh.Clear();  // the periodic mesh is no longer needed
-
-   // 3. Build the interpolator of E field
-   mesh.EnsureNodes();
-   FindPointsGSLIB E_finder(mesh);
-
-   // 4. Define finite element spaces on the parallel mesh
-   H1_FECollection phi_fec(ctx.order, ctx.dim);
-   ParFiniteElementSpace phi_fespace(&mesh, &phi_fec);
-   ND_FECollection E_fec(ctx.order, ctx.dim);
-   ParFiniteElementSpace E_fespace(&mesh, &E_fec);
-
-   // 5. Initialize the grid functions for the electric field and potential
-   ParGridFunction phi_gf(&phi_fespace);
-   ParGridFunction E_gf(&E_fespace);
-   phi_gf = 0.0;  // Initialize phi_gf to zero
-   E_gf = 0.0;    // Initialize E_gf to zero
-
-   // 6. Construct the field solver
-   FieldSolver field_solver(&phi_fespace, &E_fespace, E_finder, true);
-
-   // 7. Initialize ParticleMover
-   Ordering::Type ordering_type =
-      ctx.ordering == 0 ? Ordering::byNODES : Ordering::byVDIM;
-   int num_particles =
-      ctx.npt / num_ranks + (rank < (ctx.npt % num_ranks) ? 1 : 0);
-   ParticleMover particle_mover(MPI_COMM_WORLD, &E_gf, E_finder, num_particles,
-                                ordering_type);
-   particle_mover.InitializeChargedParticles(ctx.k, ctx.alpha, ctx.m, ctx.q,
-                                             ctx.L, ctx.reproduce);
-
-   // 8. Start the main loop
-   real_t t = 0;
-   real_t dt = ctx.dt;
-
-   mfem::StopWatch sw;
-   sw.Start();
-   for (int step = 1; step <= ctx.nt; step++)
-   {
-      // Step the FieldSolver
-      if (ctx.redist_interval > 0 &&
-          (step % ctx.redist_interval == 0 || step == 1) &&
-          particle_mover.GetParticles().GetGlobalNParticles() > 0)
-      {
-         // Redistribute
-         particle_mover.Redistribute();
-
-         // Update phi_gf from particles
-         field_solver.UpdatePhiGridFunction(particle_mover.GetParticles(),
-                                            phi_gf);
-         // Update E_gf from phi_gf
-         field_solver.UpdateEGridFunction(phi_gf, E_gf);
-
-         // Visualize fields if requested
-         if (ctx.visualization)
-         {
-            static socketstream vis_e, vis_phi;
-            common::VisualizeField(vis_e, "localhost", ctx.visport, E_gf,
-                                   "E_field", 0, 0, 500, 500);
-            common::VisualizeField(vis_phi, "localhost", ctx.visport, phi_gf,
-                                   "Potential", 500, 0, 500, 500);
-         }
-      }
-
-      // Step the ParticleMover
-      particle_mover.Step(t, dt, ctx.L, step == 1);
-      if (Mpi::Root())
-      {
-         mfem::out << "Step: " << step << " | Time: " << t;
-         mfem::out << " | Time per step: " << sw.RealTime() / step;
-         mfem::out << endl;
-      }
-      // Output particle data to CSV
-      if (ctx.output_csv_interval > 0 &&
-          (step % ctx.output_csv_interval == 0 || step == 1))
-      {
-         std::string csv_prefix = "PIC_Part_";
-         Array<int> field_idx{2}, tag_idx;
-         std::string file_name =
-            csv_prefix + mfem::to_padded_string(step, 6) + ".csv";
-         particle_mover.GetParticles().PrintCSV(file_name.c_str(), field_idx,
-                                                tag_idx);
-      }
-
-      if (ctx.redist_interval > 0 &&
-          (step % ctx.redist_interval == 0 || step == 1) &&
-          particle_mover.GetParticles().GetGlobalNParticles() > 0)
-      {
-         // Compute energies
-         // Note that particle momenta are a half time step ahead of the field
-         // after particle_mover.Step(). Therefore they are returned to the
-         // time level of the field for calculation of kinetic energy.
-         real_t kinetic_energy = particle_mover.ComputeKineticEnergy(-dt/2.);
-         real_t field_energy = field_solver.ComputeFieldEnergy(E_gf);
-
-         // Output energies
-         if (Mpi::Root())
-         {
-            cout << "Kinetic energy: " << kinetic_energy << "\t"
-                 << "Field energy: " << field_energy << "\t"
-                 << "Total energy: " << kinetic_energy + field_energy
-                 << endl;
-         }
-         // Write energies to a CSV file
-         if (Mpi::Root())
-         {
-            std::ofstream energy_file("energy.csv", std::ios::app);
-            energy_file << setprecision(10) << kinetic_energy << ","
-                        << field_energy << "," << kinetic_energy + field_energy
-                        << "\n";
-         }
-      }
-   }
-}
-
-ParticleMover::ParticleMover(MPI_Comm comm, ParGridFunction* E_gf_,
-                             FindPointsGSLIB& E_finder_, int num_particles,
-                             Ordering::Type pdata_ordering)
-   : E_gf(E_gf_), E_finder(E_finder_)
-{
-   MFEM_ASSERT(E_gf, "Must pass an E field to ParticleMover.");
-
-   int dim = E_gf->ParFESpace()->GetMesh()->SpaceDimension();
-
-   pm_.SetSize(dim);
-   pp_.SetSize(dim);
-
-   // Create particle set: 2 scalars of mass and charge,
-   // 2 vectors of size space dim for momentum and e field
-   Array<int> field_vdims({1, 1, dim, dim});
-   charged_particles = std::make_unique<ParticleSet>(
-                          comm, num_particles, dim, field_vdims, 1, pdata_ordering);
-}
-
-void ParticleMover::InitializeChargedParticles(const real_t& k,
-                                               const real_t& alpha, real_t m,
-                                               real_t q, real_t L,
-                                               bool reproduce)
-{
-   int rank;
-   MPI_Comm_rank(charged_particles->GetComm(), &rank);
-   // use time-based seed for randomness
-   std::mt19937 gen(
-      reproduce ? rank : (rank + static_cast<unsigned int>(time(nullptr))));
-   std::uniform_real_distribution<> real_dist(0.0, 1.0);
-   std::normal_distribution<> norm_dist(0.0, 1.0);
-
-   int dim = charged_particles->Coords().GetVDim();
-
-   ParticleVector& X = charged_particles->Coords();
-   ParticleVector& P = charged_particles->Field(ParticleMover::MOM);
-   ParticleVector& M = charged_particles->Field(ParticleMover::MASS);
-   ParticleVector& Q = charged_particles->Field(ParticleMover::CHARGE);
-
-   for (int i = 0; i < charged_particles->GetNParticles(); i++)
-   {
-      // Initialize momentum
-      for (int d = 0; d < dim; d++) { P(i, d) = m * norm_dist(gen); }
-
-      // Uniform positions (no accept-reject)
-      for (int d = 0; d < dim; d++) { X(i, d) = real_dist(gen) * L; }
-
-      // Displacement along x for perturbation ~ cos(k x)
-      for (int d = 0; d < dim; d++)
-      {
-         real_t x = X(i, d);
-         x -= (alpha / k) * std::sin(k * x);
-
-         // periodic wrap to [0, L)
-         x = std::fmod(x, L);
-         if (x < 0) { x += L; }
-
-         X(i, d) = x;
-      }
-
-      // Initialize mass + charge
-      M(i) = m;
-      Q(i) = q;
-   }
-   FindParticles();
-}
-
-void ParticleMover::FindParticles()
-{
-   E_finder.FindPoints(charged_particles->Coords());
-}
-
-void ParticleMover::Step(real_t& t, real_t dt, real_t L, bool first_step)
-{
-   // Update E field at particles
-   ParticleVector& E = charged_particles->Field(EFIELD);
-   E_finder.Interpolate(*E_gf, E, E.GetOrdering());
-
-   // Extract particle data
-   ParticleVector& X = charged_particles->Coords();
-   ParticleVector& P = charged_particles->Field(MOM);
-   ParticleVector& M = charged_particles->Field(MASS);
-   ParticleVector& Q = charged_particles->Field(CHARGE);
-
-   // Accelerate the particles by the electric field
-   const int npt = charged_particles->GetNParticles();
-   const int dim = X.GetVDim();
-
-   for (int particle = 0; particle < npt; ++particle)
-   {
-      for (int d = 0; d < dim; ++d)
-      {
-         P(particle, d) +=
-            (first_step ? dt / 2.0 : dt) * Q(particle) * E(particle, d);
-      }
-   }
-
-   // Periodic boundary: wrap coordinates to [0, L)
-   for (int particle = 0; particle < npt; ++particle)
-   {
-      for (int d = 0; d < dim; ++d)
-      {
-         X(particle, d) += dt / M(particle) * P(particle, d);
-         while (X(particle, d) > L) { X(particle, d) -= L; }
-         while (X(particle, d) < 0.0) { X(particle, d) += L; }
-      }
-   }
-
-   FindParticles();
-
-   // Update time
-   t += dt;
-}
-
-void ParticleMover::Redistribute()
-{
-   charged_particles->Redistribute(E_finder.GetProc());
-   FindParticles();
-}
-
-real_t ParticleMover::ComputeKineticEnergy(real_t dt) const
-{
-   const ParticleVector& P = charged_particles->Field(MOM);
-   const ParticleVector& M = charged_particles->Field(MASS);
-   const ParticleVector& Q = charged_particles->Field(CHARGE);
-   const ParticleVector& E = charged_particles->Field(EFIELD);
-
-   // Note the electric field is not reinterpolated here and the last
-   // update from Step() is used directly.
-
-   real_t kinetic_energy = 0.0;
-   for (int p = 0; p < charged_particles->GetNParticles(); ++p)
-   {
-      real_t p_square_p = 0.0;
-      for (int d = 0; d < P.GetVDim(); ++d)
-      {
-         const real_t P_m = P(p, d) + dt * Q(p) * E(p, d);
-         p_square_p += P_m * P_m;
-      }
-      kinetic_energy += 0.5 * p_square_p / M(p);
-   }
-
-   real_t global_kinetic_energy = 0.0;
-   MPI_Allreduce(&kinetic_energy, &global_kinetic_energy, 1, MPI_DOUBLE,
-                 MPI_SUM, charged_particles->GetComm());
-   return global_kinetic_energy;
-}
-
-FieldSolver::FieldSolver(ParFiniteElementSpace* phi_fes,
-                         ParFiniteElementSpace* E_fes,
-                         FindPointsGSLIB& E_finder_,
-                         bool precompute_neutralizing_const_)
-   : precompute_neutralizing_const(precompute_neutralizing_const_),
-     E_finder(E_finder_),
-     b(phi_fes)
-{
-   // compute domain volume
-   ParMesh* pmesh = phi_fes->GetParMesh();
-   real_t local_domain_volume = 0.0;
-   for (int i = 0; i < pmesh->GetNE(); i++)
-   {
-      local_domain_volume += pmesh->GetElementVolume(i);
-   }
-   MPI_Allreduce(&local_domain_volume, &domain_volume, 1, MPI_DOUBLE, MPI_SUM,
-                 phi_fes->GetParMesh()->GetComm());
-
-   {
-      // Par bilinear form for the gradgrad matrix
-      ParBilinearForm dm(phi_fes);
-      ConstantCoefficient epsilon(EPSILON);  // ε_0
-      dm.AddDomainIntegrator(
-         new DiffusionIntegrator(epsilon));  // ∫ ∇φ_i · ∇φ_j
-
-      dm.Assemble();
-      dm.Finalize();
-
-      diffusion_matrix = dm.ParallelAssemble();  // global gradgrad matrix
-   }
-
-   {
-      // Compute E = -∇φ using DiscreteLinearOperator
-      grad_interpolator = new ParDiscreteLinearOperator(phi_fes, E_fes);
-      grad_interpolator->AddDomainInterpolator(new GradientInterpolator);
-      grad_interpolator->Assemble();
-   }
-}
-
-FieldSolver::~FieldSolver()
-{
-   delete diffusion_matrix;
-   delete precomputed_neutralizing_lf;
-   delete grad_interpolator;
-}
-
-const ParLinearForm& FieldSolver::ComputeNeutralizingRHS(
-   ParFiniteElementSpace* pfes, const ParticleVector& Q, MPI_Comm comm)
-{
-   int npt = Q.Size();
-   // Get E_finder references
-   const Array<unsigned int>& code = E_finder.GetCode();
-
-   if (!precompute_neutralizing_const || precomputed_neutralizing_lf == nullptr)
-   {
-      // compute neutralizing constant
-      real_t local_sum = 0.0;
-      for (int p = 0; p < npt; ++p)
-      {
-         // Skip particles not successfully found
-         MFEM_ASSERT(code[p] != 2, "Particle " << p << " not found.");
-         local_sum += Q(p);
-      }
-
-      real_t global_sum = 0.0;
-      MPI_Allreduce(&local_sum, &global_sum, 1, MPI_DOUBLE, MPI_SUM, comm);
-
-      neutralizing_const = -global_sum / domain_volume;
-      if (Mpi::Root())
-      {
-         cout << "Total charge: " << global_sum
-              << ", Domain volume: " << domain_volume
-              << ", Neutralizing constant: " << neutralizing_const << endl;
-         if (precompute_neutralizing_const)
-         {
-            cout << "Further updates will use this precomputed neutralizing "
-                 "constant."
-                 << endl;
-         }
-      }
-      delete precomputed_neutralizing_lf;
-      precomputed_neutralizing_lf = new ParLinearForm(pfes);
-      *precomputed_neutralizing_lf = 0.0;
-      ConstantCoefficient neutralizing_coeff(neutralizing_const);
-      precomputed_neutralizing_lf->AddDomainIntegrator(
-         new DomainLFIntegrator(neutralizing_coeff));
-      precomputed_neutralizing_lf->Assemble();
-   }
-   return *precomputed_neutralizing_lf;
-}
-
-void FieldSolver::DepositCharge(ParFiniteElementSpace* pfes,
-                                const ParticleVector& Q)
-{
-   int npt = Q.Size();
-   ParMesh* pmesh = pfes->GetParMesh();
-   int dim = pmesh->SpaceDimension();
-   int curr_rank;
-   MPI_Comm_rank(pmesh->GetComm(), &curr_rank);
-
-   // Get E_finder references
-   // 0: inside, 1: boundary, 2: not found
-   const Array<unsigned int>& code = E_finder.GetCode();
-   const Array<unsigned int>& proc = E_finder.GetProc();  // owning MPI rank
-   const Array<unsigned int>& elem = E_finder.GetElem();  // local element id
-   const Vector& rref = E_finder.GetReferencePosition();  // (r,s,t) byVDIM
-
-   Array<int> dofs;
-
-   for (int p = 0; p < npt; ++p)
-   {
-      // Skip particles not successfully found
-      MFEM_ASSERT(code[p] != 2, "Particle " << p << " not found.");
-
-      // Assert particle is on the current rank
-      MFEM_ASSERT((int)proc[p] == curr_rank,
-                  "Particle " << p << " found in element owned by rank "
-                  << proc[p] << " but current rank is " << curr_rank
-                  << "." << endl
-                  << "You must call redistribute everytime before "
-                  "updating the density grid function.");
-      const int e = elem[p];
-
-      // Reference coordinates for this particle (r,s[,t]) with byVDIM layout
-      IntegrationPoint ip;
-      ip.Set(rref.GetData() + dim * p, dim);
-
-      const FiniteElement& fe = *pfes->GetFE(e);
-      const int ldofs = fe.GetDof();
-
-      Vector shape(ldofs);
-      fe.CalcShape(ip, shape);  // φ_i(x_p) in this element
-
-      pfes->GetElementDofs(e, dofs);  // local dof indices
-
-      const real_t q_p = Q(p);
-
-      // Add q_p * φ_i(x_p) to b_i
-      b.AddElementVector(dofs, q_p, shape);
-   }
-}
-
-void FieldSolver::UpdatePhiGridFunction(ParticleSet& particles,
-                                        ParGridFunction& phi_gf)
-{
-   // FE space / mesh
-   ParFiniteElementSpace* pfes = phi_gf.ParFESpace();
-
-   // Particle data: Q - charges (npt x 1)
-   ParticleVector& Q = particles.Field(ParticleMover::CHARGE);
-
-   // --------------------------------------------------------
-   // 1) Make RHS and pre-subtract averaged charge density for zero-mean RHS
-   // --------------------------------------------------------
-   MPI_Comm comm = pfes->GetComm();
-   b = ComputeNeutralizingRHS(pfes, Q, comm);
-
-   // --------------------------------------------------------
-   // 2) Deposit q_p * phi_i(x_p) into a ParLinearForm (RHS b)
-   //      b_i = sum_p q_p * φ_i(x_p)
-   // --------------------------------------------------------
-   DepositCharge(pfes, Q);
-
-   // Assemble to a global true-dof RHS vector compatible with MassMatrix
-   HypreParVector B(pfes);
-   b.ParallelAssemble(B);
-
-   // ------------------------------------------------------------------
-   // 3) Solve A * phi = B with zero-mean enforcement via OrthoSolver
-   // ------------------------------------------------------------------
-   phi_gf = 0.0;
-   HypreParVector Phi_true(pfes);
-   Phi_true = 0.0;
-
-   HyprePCG solver(diffusion_matrix->GetComm());
-   solver.SetOperator(*diffusion_matrix);
-   solver.SetTol(1e-12);
-   solver.SetMaxIter(200);
-   solver.SetPrintLevel(0);
-
-   HypreBoomerAMG prec(*diffusion_matrix);
-   prec.SetPrintLevel(0);
-   solver.SetPreconditioner(prec);
-
-   OrthoSolver ortho(comm);
-   ortho.SetSolver(solver);
-   ortho.Mult(B, Phi_true);
-
-   // Map true-dof solution back to the ParGridFunction
-   phi_gf.Distribute(Phi_true);
-}
-
-void FieldSolver::UpdateEGridFunction(ParGridFunction& phi_gf,
-                                      ParGridFunction& E_gf)
-{
-   // Compute ∇φ using precomputed gradient operator
-   grad_interpolator->Mult(phi_gf, E_gf);
-   // Scale by -1 to get E = -∇φ
-   E_gf.Neg();
-}
-
-real_t FieldSolver::ComputeFieldEnergy(const ParGridFunction& E_gf) const
-{
-   // ---- Field energy: 0.5 * ∫ ||E||^2 dx ----
-   const ParFiniteElementSpace* fes = E_gf.ParFESpace();
-   const ParMesh* pmesh = fes->GetParMesh();
-
-   const int order = fes->GetMaxElementOrder();
-   const int qorder = std::max(2, 2 * order + 1);
-
-   const IntegrationRule* irs[Geometry::NumGeom];
-   for (int g = 0; g < Geometry::NumGeom; g++)
-   {
-      irs[g] = &IntRules.Get(g, qorder);
-   }
-
-   real_t field_energy = 0.0;
-
-   Vector zero(pmesh->Dimension());
-   zero = 0.0;
-   VectorConstantCoefficient zero_vec(zero);
-
-   const real_t E_l2 = E_gf.ComputeL2Error(zero_vec, irs);
-   field_energy = 0.5 * EPSILON * E_l2 * E_l2;
-
-   return field_energy;
-}
-
-void display_banner(ostream& os)
-{
-   os << R"(
-      ██████╗░██╗░█████╗░
-      ██╔══██╗██║██╔══██╗
-      ██████╔╝██║██║░░╚═╝
-      ██╔═══╝░██║██║░░██╗
-      ██║░░░░░██║╚█████╔╝
-      ╚═╝░░░░░╚═╝░╚════╝░
-         )"
-      << endl
-      << flush;
-}

miniapps/plasma/pic/makefile

@@ -1,85 +0,0 @@
-# Copyright (c) 2010-2025, 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 ?= ../../..
-MFEM_INSTALL_DIR ?= ../../../mfem
-SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/plasma/pic/,)
-CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
-   $(wildcard $(MFEM_INSTALL_DIR)/share/mfem/config.mk))
-
-MFEM_LIB_FILE = mfem_is_not_built
--include $(CONFIG_MK)
-
-PAR_MINIAPPS =
-
-ifeq ($(MFEM_USE_GSLIB),YES)
-   PAR_MINIAPPS += electrostatic-pic
-endif
-
-ifeq ($(MFEM_USE_MPI),NO)
-   MINIAPPS =
-else
-   MINIAPPS = $(PAR_MINIAPPS)
-endif
-
-.SUFFIXES:
-.SUFFIXES: .o .cpp .mk
-.PHONY: all lib-common clean clean-build clean-exec
-.PRECIOUS: %.o
-
-COMMON_LIB = -L$(MFEM_BUILD_DIR)/miniapps/common -lmfem-common
-
-# If MFEM_SHARED is set, add the ../common rpath
-COMMON_LIB += $(if $(MFEM_SHARED:YES=),,\
-   $(MFEM_XLINKER)-rpath,$(abspath $(MFEM_BUILD_DIR)/miniapps/common))
-
-# Remove built-in rules
-%: %.cpp
-%.o: %.cpp
-
-all: $(MINIAPPS)
-
-# Rules for building the miniapps
-electrostatic-pic: electrostatic-pic.cpp $(MFEM_LIB_FILE) $(CONFIG_MK) | lib-common
-	$(MFEM_CXX) $(MFEM_FLAGS) -c $<
-	$(MFEM_CXX) $(MFEM_LINK_FLAGS) -o $@ $@.o $(COMMON_LIB) $(MFEM_LIBS)
-
-# Rule for building lib-common
-lib-common:
-	$(MAKE) -C $(MFEM_BUILD_DIR)/miniapps/common
-
-
-MFEM_TESTS = MINIAPPS
-include $(MFEM_TEST_MK)
-
-# Testing: "test" target and mfem-test* variables are defined in config/test.mk
-
-# Testing: Specific execution options
-RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
-electrostatic-pic-test-par: electrostatic-pic
-	@$(call mfem-test,$<, $(RUN_MPI), PIC miniapp,\
-		-rdi 2 -npt 40960 -k 0.2855993321 -a 0.05 -nt 200 -nx 16 -ny 16\
-		-O 1 -q 0.01181640625 -m 0.01181640625 -oci 1000 -dt 0.1)
-
-# 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_MINIAPPS) $(PAR_MINIAPPS)
-	rm -rf *.dSYM *.TVD.*breakpoints
-
-clean-exec:
-	@rm -rf electrostatic-pic_* *.csv energy.csv

tests/unit/CMakeLists.txt

@@ -71,12 +71,10 @@ set(UNIT_TESTS_SRCS
   linalg/test_ode2.cpp
   linalg/test_operator.cpp
   linalg/test_particlevector.cpp
-  linalg/test_petsc_nonlinear.cpp
   linalg/test_sparsesmoothers.cpp
   linalg/test_vector.cpp
   mesh/mesh_test_utils.cpp
   mesh/test_exodus_reader.cpp
-  mesh/test_mfem_mesh_reader.cpp
   mesh/test_exodus_writer.cpp
   mesh/test_face_orientations.cpp
   mesh/test_fms.cpp

tests/unit/data/compass-testing.mesh

@@ -1,118 +0,0 @@
-MFEM mesh v1.3
-
-#
-# MFEM Geometry Types (see mesh/geom.hpp):
-#
-# POINT       = 0
-# SEGMENT     = 1
-# TRIANGLE    = 2
-# SQUARE      = 3
-# TETRAHEDRON = 4
-# CUBE        = 5
-# PRISM       = 6
-#
-
-dimension
-2
-
-elements
-12
-10 2 7 0 1
-11 2 0 7 2
-12 2 9 0 2
-13 2 0 9 3
-14 2 11 0 3
-15 2 0 11 4
-16 2 5 0 4
-17 2 0 5 1
-9 3 1 5 6 7
-9 3 2 7 8 9
-9 3 3 9 10 11
-9 3 4 11 12 5
-
-attribute_sets
-16
-"Base" 1 9
-"E Even" 1 16
-"E Odd" 1 17
-"East"
-2
-16
-17
-"N Even" 1 10
-"N Odd" 1 11
-"North" 2 10     11
-"Rose" 8 10 11     12
-13    14 
-15    16     17
-"Rose Even" 4
-10
-12
-14
-16
-"Rose Odd"
-4
-11
-13
-15
-17
-"S Even" 1 14
-"S Odd" 1 15
-South 2
-14
-15
-"W Even" 1 12
-"W Odd" 1 13
-West 2 12 13
-
-boundary
-8
-1 1 5 6
-2 1 6 7
-3 1 7 8
-4 1 8 9
-5 1 9 10
-6 1 10 11
-7 1 11 12
-8 1 12 5
-
-bdr_attribute_sets
-13
-"Boundary" 8 1 2 3 4 5 6 7 8
-"ENE" 1 1
-"ESE" 1 8
-"Eastern Boundary" 2 1 8
-"NNE" 1 2
-"NNW" 1 3
-"Northern Boundary"
-2
-2
-3
-"SSE" 1 7
-"SSW" 1 6
-"Southern Boundary" 2
-6
-7
-"WNW" 1 4
-"WSW" 1 5
-"Western Boundary" 2 4
-5
-
-vertices
-13
-2
-0 0
-0.14142136 0.14142136
--0.14142136 0.14142136
--0.14142136 -0.14142136
-0.14142136 -0.14142136
-1 0
-0.70710678 0.70710678
-0 1
--0.70710678 0.70710678
--1 0
--0.70710678 -0.70710678
-0 -1
-0.70710678 -0.70710678
-
-mfem_mesh_end

tests/unit/fem/test_particleset.cpp

@@ -296,7 +296,7 @@ void TestRedistribute(Ordering::Type ordering)
       int wrong_proc_count = 0;
       for (int i = 0; i < procs.Size(); i++)
       {
-         if (static_cast<unsigned>(rank) != procs[i])
+         if (rank != procs[i])
          {
             wrong_proc_count++;
          }

tests/unit/fem/test_var_order.cpp

@@ -271,8 +271,6 @@ TEST_CASE("Variable Order FiniteElementSpace",
 
       const auto space_type = GENERATE(SpaceType::RT, SpaceType::ND);
       const int dim = GENERATE(2, 3);
-      CAPTURE(space_type);
-      CAPTURE(dim);
 
       Mesh mesh = MakeCartesianMesh(dim == 2 ? 4 : 2, dim);
       mesh.EnsureNCMesh();
@@ -700,14 +698,7 @@ static void TestSolveVec(FiniteElementSpace &fespace)
 
    GridFunction x(&fespace);
    x = 0.0;
-   if (x.FESpace()->GetTypicalBE()->GetRangeDim() == 0)
-   {
-      x.ProjectBdrCoefficientNormal(exsol, ess_attr);
-   }
-   else
-   {
-      x.ProjectBdrCoefficientTangent(exsol, ess_attr);
-   }
+   x.ProjectBdrCoefficient(exsol, ess_attr);
 
    // Assemble the linear form
    LinearForm lf(&fespace);
@@ -1091,14 +1082,7 @@ static void TestSolveParVec(ParFiniteElementSpace &fespace)
 
    ParGridFunction x(&fespace);
    x = 0.0;
-   if (x.FESpace()->GetTypicalBE()->GetRangeDim() == 0)
-   {
-      x.ProjectBdrCoefficientNormal(exsol, ess_attr);
-   }
-   else
-   {
-      x.ProjectBdrCoefficientTangent(exsol, ess_attr);
-   }
+   x.ProjectBdrCoefficient(exsol, ess_attr);
 
    // Assemble the linear form
    ParLinearForm lf(&fespace);

tests/unit/general/test_text.cpp

@@ -30,29 +30,3 @@ TEST_CASE("String Manipulation", "[General]")
       }
    }
 }
-
-TEST_CASE("Quoted String Input", "[General]")
-{
-   const auto test_strings =
-   {
-      "Test",
-      "Test with spaces",
-      "Test with \"quoted text\"",
-      "Test string ending with \\",
-      "\nTest with\tvarious white\v\rspace characters.",
-      "Test with some unicode characters: ∆, ∉, ∑, 🍎."
-   };
-
-   for (const auto c_str : test_strings)
-   {
-      CAPTURE(c_str);
-      const std::string str(c_str);
-      std::stringstream ss;
-      ss << std::quoted(str);
-
-      std::string read_str;
-      int error = parse_quoted_string(read_str, ss);
-      CHECK(error == 0);
-      CHECK(read_str == str);
-   }
-}

tests/unit/linalg/test_petsc_nonlinear.cpp

@@ -1,74 +0,0 @@
-// Copyright (c) 2010-2025, 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 "mfem.hpp"
-#include "unit_tests.hpp"
-
-using namespace mfem;
-
-#if defined(MFEM_USE_MPI) && defined(MFEM_USE_PETSC)
-
-namespace
-{
-struct PetscSession
-{
-   PetscSession() { MFEMInitializePetsc(); }
-   ~PetscSession() { MFEMFinalizePetsc(); }
-};
-
-class IdentityGradientOperator : public IdentityOperator
-{
-public:
-   IdentityGradientOperator() : IdentityOperator(1), _jac(1)
-   {
-      _jac.Add(0, 0, 1.0);
-      _jac.Finalize();
-   }
-
-   Operator &GetGradient(const Vector &) const override
-   {
-      return const_cast<SparseMatrix &>(_jac);
-   }
-
-private:
-   SparseMatrix _jac;
-};
-}
-
-TEST_CASE("PetscNonlinearSolver accepts non-empty rhs", "[Parallel][PETSc]")
-{
-   static PetscSession petsc_session;
-
-   IdentityGradientOperator oper;
-   PetscNonlinearSolver solver(MPI_COMM_WORLD, "nl_");
-   solver.SetRelTol(1.0e-12);
-   solver.SetAbsTol(1.0e-12);
-   solver.SetMaxIter(5);
-   solver.SetPrintLevel(0);
-   solver.SetJacobianType(Operator::PETSC_MATAIJ);
-   solver.SetOperator(oper);
-
-   Vector x(1);
-
-   Vector empty_rhs;
-   x = 0.0;
-   solver.Mult(empty_rhs, x);
-   REQUIRE(x(0) == MFEM_Approx(0.0));
-
-   Vector nonempty_rhs(1);
-   nonempty_rhs(0) = 2.5;
-   x = 0.0;
-   solver.Mult(nonempty_rhs, x);
-   REQUIRE(x.Size() == 1);
-   REQUIRE(x(0) == MFEM_Approx(nonempty_rhs(0)));
-}
-
-#endif

tests/unit/mesh/test_mfem_mesh_reader.cpp

@@ -1,108 +0,0 @@
-// Copyright (c) 2010-2025, 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 "mfem.hpp"
-#include "unit_tests.hpp"
-
-#include <algorithm>
-#include <string>
-#include <utility>
-#include <vector>
-
-using namespace mfem;
-
-TEST_CASE("MFEM Mesh Named Attributes", "[Mesh]")
-{
-   // Path relative to the directory tests/unit
-   Mesh mesh("data/compass-testing.mesh");
-
-   REQUIRE(mesh.Dimension() == 2);
-   REQUIRE(mesh.GetNE() == 12);
-   REQUIRE(mesh.GetNV() == 13);
-
-   REQUIRE(mesh.attribute_sets.attr_sets.Size() == 16);
-   REQUIRE(mesh.bdr_attribute_sets.attr_sets.Size() == 13);
-
-   std::vector<std::pair<std::string, std::vector<int>>> expected_attr_sets =
-   {
-      {"Base",  {9}},
-      {"E Even",  {16}},
-      {"E Odd",  {17}},
-      {"East",  {16, 17}},
-      {"N Even",  {10}},
-      {"N Odd",  {11}},
-      {"North",  {10, 11}},
-      {"Rose",  {10, 11, 12, 13, 14, 15, 16, 17}},
-      {"Rose Even",  {10, 12, 14, 16}},
-      {"Rose Odd",  {11, 13, 15, 17}},
-      {"S Even",  {14}},
-      {"S Odd",  {15}},
-      {"South",  {14, 15}},
-      {"W Even",  {12}},
-      {"W Odd",  {13}},
-      {"West",  {12, 13}}
-   };
-
-   for (auto const &attr_name_index_pair: expected_attr_sets )
-   {
-      REQUIRE(mesh.attribute_sets.AttributeSetExists(
-                 attr_name_index_pair.first));
-
-      auto const &attr_set = mesh.attribute_sets.GetAttributeSet(
-                                attr_name_index_pair.first);
-      auto const &expected_attr_set = attr_name_index_pair.second;
-
-      REQUIRE(static_cast<std::size_t>(attr_set.Size()) ==
-              expected_attr_set.size());
-
-      bool const elements_equal = std::equal(attr_set.begin(), attr_set.end(),
-                                             expected_attr_set.begin());
-
-      REQUIRE(elements_equal);
-   }
-
-   std::vector<std::pair<std::string, std::vector<int>>> expected_bdr_attr_sets
-   =
-   {
-      {"Boundary",  {1, 2, 3, 4, 5, 6, 7, 8}},
-      {"ENE",  { 1}},
-      {"ESE",  { 8}},
-      {"Eastern Boundary",  {1, 8}},
-      {"NNE",  { 2}},
-      {"NNW",  { 3}},
-      {"Northern Boundary", {2, 3}},
-      {"SSE",  { 7}},
-      {"SSW",  { 6}},
-      {"Southern Boundary", {6,7}},
-      {"WNW",  { 4}},
-      {"WSW",  { 5}},
-      {"Western Boundary", {4,5}}
-   };
-
-   for (auto const &attr_bdr_name_index_pair: expected_bdr_attr_sets )
-   {
-      REQUIRE(mesh.bdr_attribute_sets.AttributeSetExists(
-                 attr_bdr_name_index_pair.first));
-
-      auto const &bdr_attr_set = mesh.bdr_attribute_sets.GetAttributeSet(
-                                    attr_bdr_name_index_pair.first);
-      auto const &expected_bdr_attr_set = attr_bdr_name_index_pair.second;
-
-      REQUIRE(static_cast<std::size_t>(bdr_attr_set.Size()) ==
-              expected_bdr_attr_set.size());
-
-      bool const elements_equal = std::equal(bdr_attr_set.begin(),
-                                             bdr_attr_set.end(),
-                                             expected_bdr_attr_set.begin());
-
-      REQUIRE(elements_equal);
-   }
-}

tests/unit/mesh/test_psubmesh.cpp

@@ -486,14 +486,8 @@ void multidomain_test_3d(FECType fec_type)
          {
             cylinder_gf.ProjectCoefficient(vcoeff);
             outer_gf.ProjectCoefficient(vcoeff);
-            if (fec_type == FECType::RT)
-            {
-               outer_gf.ProjectBdrCoefficientNormal(vzerocoeff, outer_cyl_surf_marker);
-            }
-            else
-            {
-               outer_gf.ProjectBdrCoefficientTangent(vzerocoeff, outer_cyl_surf_marker);
-            }
+            outer_gf.ProjectBdrCoefficient(vzerocoeff,
+                                           outer_cyl_surf_marker);
             outer_gf_ex.ProjectCoefficient(vcoeff);
          }
          ParSubMesh::Transfer(cylinder_gf, outer_gf);
@@ -513,14 +507,8 @@ void multidomain_test_3d(FECType fec_type)
          {
             outer_gf.ProjectCoefficient(vcoeff);
             cylinder_gf.ProjectCoefficient(vcoeff);
-            if (fec_type == FECType::RT)
-            {
-               cylinder_gf.ProjectBdrCoefficientNormal(vzerocoeff, cylinder_cyl_surf_marker);
-            }
-            else
-            {
-               cylinder_gf.ProjectBdrCoefficientTangent(vzerocoeff, cylinder_cyl_surf_marker);
-            }
+            cylinder_gf.ProjectBdrCoefficient(vzerocoeff,
+                                              cylinder_cyl_surf_marker);
             cylinder_gf_ex.ProjectCoefficient(vcoeff);
          }
          ParSubMesh::Transfer(outer_gf, cylinder_gf);