minor

AI policy

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

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

.gitignore

@@ -443,6 +443,10 @@ 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

.gitlab/jobs/dane.yml

@@ -85,3 +85,8 @@ opt_par_gcc_10_pumi:
   extends: .mfem_job_on_dane
   variables:
     SPEC: "%gcc@10.3.1 +pumi"
+
+opt_par_gcc_10_gslib:
+  extends: .mfem_job_on_dane
+  variables:
+    SPEC: "%gcc@10.3.1 +gslib"

.gitlab/jobs/matrix.yml

@@ -63,3 +63,8 @@ opt_mpi_cuda_hypre_cuda_gcc:
   extends: .mfem_job_on_matrix
   variables:
     SPEC: "%gcc@10.3.1 +mpi +cuda cuda_arch=90 ^hypre+cuda"
+
+opt_mpi_cuda_gcc_gslib:
+  extends: .mfem_job_on_matrix
+  variables:
+    SPEC: "%gcc@10.3.1 +mpi +cuda +gslib cuda_arch=90 ^hypre+cuda"

.gitlab/scripts/baseline

@@ -32,9 +32,9 @@ mkdir _${BASELINE_TEST} && cd _${BASELINE_TEST}
 
 # run
 if [[ "${MACHINE_NAME}" == "dane" ]]; then
-  salloc --nodes=1 -t 60 --exclusive --reservation=ci ../runtest ../../mfem "${BASELINE_TEST} ${TPLS_DIR}"
+  srun --nodes=1 -t 60 --exclusive --reservation=ci ../runtest ../../mfem "${BASELINE_TEST} ${TPLS_DIR}"
 elif [[ ${MACHINE_NAME} == "corona" ]]; then
-  salloc --nodes=1 -t 60 -p pbatch ../runtest ../../mfem "${BASELINE_TEST} ${TPLS_DIR}"
+  srun --nodes=1 -t 60 -p pbatch ../runtest ../../mfem "${BASELINE_TEST} ${TPLS_DIR}"
 else
   echo "Unknown machine: MACHINE_NAME=$MACHINE_NAME"
   exit 1

CHANGELOG

@@ -8,8 +8,10 @@
                                https://mfem.org
 
 
-Version 4.10 (development)
-==========================
+Version 4.9.1 (development)
+===========================
+
+- Policy for AI-assisted contribution added to CONTRIBUTING.md
 
 Discretization improvements
 ---------------------------
@@ -23,12 +25,6 @@ Discretization improvements
   Tet rules (d=14-20): Chuluunbaatar et al., Comput. Math. Appl. 124:89-97,
   2022.
 
-
-Version 4.9.1 (development)
-===========================
-
-Discretization improvements
----------------------------
 - Improved the gridfunction projection routines. Projections work for Scalar,
   Vector and VectorFE, also NURBS versions. Optionally different types of
   projections can be selected, default behaviour has not changed.

CONTRIBUTING.md

@@ -24,6 +24,14 @@ must be made under this license.
 Note also that MFEM has a [Code of Conduct](CODE_OF_CONDUCT.md). By participating
 in the MFEM community, you agree to abide by its rules.
 
+## AI Policy
+- Use of AI code generation in MFEM is allowed but must be disclosed, e.g. by
+  selecting the `AI-assisted` label on the PR.
+- By submitting a PR, the author acknowledges that they have reviewed and
+  understand the changes they are proposing.
+- PR authors are still responsible for correctness, licensing, and attribution
+  of all changes.
+
 If you plan on contributing to MFEM, consider reviewing the
 [issue tracker](https://github.com/mfem/mfem/issues) first to check if a thread
 already exists for your desired feature or the bug you ran into. Use a pull

config/cmake/MFEMConfig.cmake.in

@@ -109,6 +109,10 @@ 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_DIR AND UMPIRE_DIR)
-  set(umpire_DIR ${UMPIRE_DIR}/lib/cmake/umpire)
+if (NOT umpire_ROOT AND UMPIRE_DIR)
+  set(umpire_ROOT ${UMPIRE_DIR})
 endif()
 message(STATUS "Looking for UMPIRE ...")
 message(STATUS "   in UMPIRE_DIR = ${UMPIRE_DIR}")
-message(STATUS "      umpire_DIR = ${umpire_DIR}")
+message(STATUS "      umpire_ROOT = ${umpire_ROOT}")
 find_package(umpire CONFIG)
 set(UMPIRE_FOUND ${umpire_FOUND})
 set(UMPIRE_LIBRARIES "umpire")

examples/CMakeLists.txt

@@ -215,7 +215,7 @@ if (MFEM_ENABLE_TESTING)
     add_test(NAME ex1p_ceed_np=${MFEM_MPI_NP}
       COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
       ${MPIEXEC_PREFLAGS}
-      $<TARGET_FILE:ex1p> "-no-vis" "-d ceed-cpu" "-pa" "-a"
+      $<TARGET_FILE:ex1p> "-no-vis" "-d" "ceed-cpu" "-pa" "-a"
       ${MPIEXEC_POSTFLAGS})
   endif()
 endif()

examples/amgx/makefile

@@ -64,7 +64,7 @@ PARALLEL_NAME := Parallel AMGX example
 $(MFEM_LIB_FILE):
 	$(error The MFEM library is not build)
 
-clean: clean-build
+clean: clean-build clean-exec
 
 clean-build:
 	rm -f *.o *~ $(SEQ_EXAMPLES) $(PAR_EXAMPLES)

examples/caliper/makefile

@@ -64,12 +64,12 @@ ex1p-test-par: ex1p
 $(MFEM_LIB_FILE):
 	$(error The MFEM library is not built)
 
-clean: clean-build clean-exec $(SUBDIRS_CLEAN)
+clean: clean-build clean-exec
 
 clean-build:
 	rm -f *.o *~ $(SEQ_EXAMPLES) $(PAR_EXAMPLES)
 	rm -rf *.dSYM *.TVD.*breakpoints
 
 clean-exec:
-	@rm -f refined.mesh displaced.mesh mesh.* ex5.mesh
-	@rm -f sphere_refined.* sol.* sol_u.* sol_p.* sol_r.* sol_i.*
+	@rm -f refined.mesh mesh.*
+	@rm -f sol.*

examples/ex36.cpp

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

examples/ex36p.cpp

@@ -103,7 +103,19 @@ int main(int argc, char *argv[])
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
-   args.ParseCheck();
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
 
    // 2. Read the mesh from the mesh file.
    const char *mesh_file = "../data/disc-nurbs.mesh";
@@ -131,7 +143,7 @@ int main(int argc, char *argv[])
    mesh.Clear();
 
    // 4. Define the necessary finite element spaces on the mesh.
-   H1Bubble_FECollection H1fec(order, order - 1, dim);
+   H1_FECollection H1fec(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
-//     ex37 -o 2 -alpha 5.0 -mi 50 -vf 0.4 -ntol 1e-5
+//     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 -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,53 +55,6 @@
 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
@@ -180,10 +133,11 @@ 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-1;
+   real_t itol = 1e-2;
    real_t ntol = 1e-4;
    real_t rho_min = 1e-6;
    real_t lambda = 1.0;
@@ -198,6 +152,8 @@ 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",
@@ -332,6 +288,7 @@ 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)));
@@ -385,7 +342,7 @@ int main(int argc, char *argv[])
    // 11. Iterate:
    for (int k = 1; k <= max_it; k++)
    {
-      if (k > 1) { alpha *= ((real_t) k) / ((real_t) k-1); }
+      if (k > 1) { alpha = std::pow((real_t) k,growth); }
 
       mfem::out << "\nStep = " << k << std::endl;
 
@@ -422,7 +379,9 @@ int main(int argc, char *argv[])
 
       // Step 5 - Update design variable ψ ← proj(ψ - αG)
       psi.Add(-alpha, grad);
-      const real_t material_volume = proj(psi, target_volume);
+      GridFunction alpha_grad(grad);
+      alpha_grad *= alpha;
+      const real_t material_volume = proj(psi, alpha_grad, 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,9 +231,12 @@ 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;
@@ -267,6 +270,8 @@ public:
    void ResetFEM();
    void SetupFEM();
 
+   void UpdateEssentialTDofs();
+   void AssembleDiffusionBilinear(bool update_ess_tdofs=true);
    void Solve();
    GridFunction * GetFEMSolution();
    LinearForm * GetLinearForm() {return b;}
@@ -371,6 +376,130 @@ 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
 
@@ -422,12 +551,8 @@ void DiffusionSolver::SetupFEM()
    }
 }
 
-void DiffusionSolver::Solve()
+void DiffusionSolver::UpdateEssentialTDofs()
 {
-   OperatorPtr A;
-   Vector B, X;
-   Array<int> ess_tdof_list;
-
 #ifdef MFEM_USE_MPI
    if (parallel)
    {
@@ -440,7 +565,39 @@ void DiffusionSolver::Solve()
 #else
    fes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list);
 #endif
-   *u=0.0;
+}
+
+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;
+
    if (b)
    {
       delete b;
@@ -475,31 +632,33 @@ void DiffusionSolver::Solve()
 
    b->Assemble();
 
-   BilinearForm * a = nullptr;
-
-#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();
+   *u=0.0;
    if (essbdr_cf)
    {
       u->ProjectBdrCoefficient(*essbdr_cf,ess_bdr);
    }
-   a->FormLinearSystem(ess_tdof_list, *u, *b, A, X, B);
+
+#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);
+   }
+   else
+   {
+      X.NewDataAndSize(u->GetData(), u->Size());
+      B.NewDataAndSize(b->GetData(), b->Size());
+      a->EliminateVDofsInRHS(ess_tdof_list, X, B);
+   }
+#else
+   X.NewDataAndSize(u->GetData(), u->Size());
+   B.NewDataAndSize(b->GetData(), b->Size());
+   a->EliminateVDofsInRHS(ess_tdof_list, X, B);
+#endif
 
    CGSolver * cg = nullptr;
    Solver * M = nullptr;
@@ -528,7 +687,6 @@ void DiffusionSolver::Solve()
    delete M;
    delete cg;
    a->RecoverFEMSolution(X, *b, *u);
-   delete a;
 }
 
 GridFunction * DiffusionSolver::GetFEMSolution()
@@ -560,6 +718,8 @@ 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
-//    mpirun -np 4 ex37p -o 2 -alpha 5.0 -mi 50 -vf 0.4 -ntol 1e-5
+//    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 -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,61 +54,6 @@
 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
@@ -193,10 +138,11 @@ 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-1;
+   real_t itol = 1e-2;
    real_t ntol = 1e-4;
    real_t rho_min = 1e-6;
    real_t lambda = 1.0;
@@ -211,6 +157,8 @@ 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",
@@ -359,6 +307,7 @@ 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)));
@@ -412,7 +361,7 @@ int main(int argc, char *argv[])
    // 11. Iterate:
    for (int k = 1; k <= max_it; k++)
    {
-      if (k > 1) { alpha *= ((real_t) k) / ((real_t) k-1); }
+      if (k > 1) { alpha = std::pow((real_t) k,growth); }
 
       if (myid == 0)
       {
@@ -452,7 +401,9 @@ int main(int argc, char *argv[])
 
       // Step 5 - Update design variable ψ ← proj(ψ - αG)
       psi.Add(-alpha, grad);
-      const real_t material_volume = proj(psi, target_volume);
+      ParGridFunction alpha_grad(grad);
+      alpha_grad *= alpha;
+      const real_t material_volume = proj(psi, alpha_grad, target_volume);
 
       // Compute ||ρ - ρ_old|| in control fes.
       real_t norm_increment = zerogf.ComputeL1Error(succ_diff_rho);

examples/ginkgo/makefile

@@ -76,4 +76,4 @@ clean-build:
 	rm -rf *.dSYM *.TVD.*breakpoints
 
 clean-exec:
-	@rm -f refined.mesh sol.gf
+	@rm -f refined.mesh sol.gf mesh.* sol.*

examples/makefile

@@ -71,6 +71,7 @@ endif
 
 SUBDIRS_ALL = $(addsuffix /all,$(SUBDIRS))
 SUBDIRS_TEST = $(addsuffix /test,$(SUBDIRS))
+SUBDIRS_TEST_NOCLEAN = $(addsuffix /test-noclean,$(SUBDIRS))
 SUBDIRS_CLEAN = $(addsuffix /clean,$(SUBDIRS))
 SUBDIRS_TPRINT = $(addsuffix /test-print,$(SUBDIRS))
 
@@ -87,8 +88,9 @@ SUBDIRS_TPRINT = $(addsuffix /test-print,$(SUBDIRS))
 
 all: $(EXAMPLES) $(SUBDIRS_ALL)
 
-.PHONY: $(SUBDIRS_ALL) $(SUBDIRS_TEST) $(SUBDIRS_CLEAN) $(SUBDIRS_TPRINT)
-$(SUBDIRS_ALL) $(SUBDIRS_TEST) $(SUBDIRS_CLEAN):
+.PHONY: $(SUBDIRS_ALL) $(SUBDIRS_TEST) $(SUBDIRS_TEST_NOCLEAN) \
+   $(SUBDIRS_CLEAN) $(SUBDIRS_TPRINT)
+$(SUBDIRS_ALL) $(SUBDIRS_TEST) $(SUBDIRS_TEST_NOCLEAN) $(SUBDIRS_CLEAN):
 	$(MAKE) -C $(@D) $(@F)
 $(SUBDIRS_TPRINT):
 	@$(MAKE) -C $(@D) $(@F)
@@ -107,6 +109,7 @@ endif
 MFEM_TESTS = EXAMPLES
 include $(MFEM_TEST_MK)
 test: $(SUBDIRS_TEST)
+test-noclean: $(SUBDIRS_TEST_NOCLEAN)
 test-print: $(SUBDIRS_TPRINT)
 
 # Testing: Parallel vs. serial runs
@@ -157,6 +160,8 @@ 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,7 +73,6 @@ 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
@@ -134,7 +133,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
@@ -222,7 +221,6 @@ 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/bounds.cpp

@@ -41,9 +41,14 @@ void PLBound::Setup(const int nb_i, const int ncp_i,
    tol = tol_i;
    lbound.SetSize(ncp, nb);
    ubound.SetSize(ncp, nb);
+   lbound_t.SetSize(nb, ncp);
+   ubound_t.SetSize(nb, ncp);
    nodes.SetSize(nb);
    weights.SetSize(nb);
    control_points.SetSize(ncp);
+   xhat.SetSize(nb);
+   what.SetSize(nb);
+   cphat.SetSize(ncp);
 
    auto scalenodes = [](const Vector &in, const real_t a, const real_t b) -> Vector
    {
@@ -90,6 +95,10 @@ void PLBound::Setup(const int nb_i, const int ncp_i,
       MFEM_ABORT("Unsupported interval points. Use [0,1].\n");
    }
    control_points = scalenodes(control_points, 0.0, 1.0); // rescale to [0,1]
+   for (int i = 0; i < ncp; i++)
+   {
+      cphat(i) = 2.0*control_points(i) - 1.0;
+   }
 
    Poly_1D::Basis &basis1d(poly1d.GetBasis(nb-1, b_type));
 
@@ -145,6 +154,8 @@ void PLBound::Setup(const int nb_i, const int ncp_i,
                lbound(j,i) = std::max(lbound(j,i),0_r);
             }
          }
+         lbound_t(i,j) = lbound(j,i);
+         ubound_t(i,j) = ubound(j,i);
       }
    }
 
@@ -176,6 +187,11 @@ void PLBound::Setup(const int nb_i, const int ncp_i,
          nodes(i) = irule.IntPoint(i).x;
       }
    }
+   for (int i = 0; i < nb; i++)
+   {
+      xhat(i) = 2.0*nodes(i) - 1.0;
+      what(i) = 2.0*weights(i);
+   }
 
    if (b_type == 2)
    {
@@ -755,4 +771,4 @@ void PLBound::Print(std::ostream &outp) const
    ubound.Print(outp);
 }
 
-}
+}

fem/bounds.hpp

@@ -13,6 +13,7 @@
 #define MFEM_BOUNDS
 
 #include "../config/config.hpp"
+#include "../general/forall.hpp"
 #include "fespace.hpp"
 
 namespace mfem
@@ -60,7 +61,9 @@ private:
    bool proj = true; // Use linear projection to compute bounds.
    real_t tol = 0.0; // offset bounds to avoid round-off errors
    Vector nodes, weights, control_points;
+   Vector xhat, what, cphat;
    DenseMatrix lbound, ubound; // ncp x nb matrices with bounds of all bases
+   DenseMatrix lbound_t, ubound_t; // nb x ncp transposes for device kernel
    // Some auxillary storage for computing the bounds with Bernstein
    DenseMatrix basisMatNodes; // Bernstein bases at equispaced nodes
    DenseMatrix basisMatInt;   // Bernstein bases at GLL nodes
@@ -113,7 +116,10 @@ public:
     *  @details This projection increases the computational cost but results in
     *  tighter bounds.
     */
-   void SetProjectionFlagForBounding(bool proj_) { proj = proj_; }
+   void SetProjectionFlagForBounding(bool proj_)
+   {
+      proj = proj_;
+   }
 
    /** @brief Compute piecewise linear bounds for the lexicographically-ordered
     *  nodal coefficients in @a coeff in 1D/2D/3D.
@@ -137,9 +143,23 @@ public:
    /// Get number of control points used to compute the bounds.
    int GetNControlPoints() const { return ncp; }
 
+   /// Get the underlying 1D basis type.
+   int GetBasisType() const { return b_type; }
+
    /// Get 1D control point locations (lexicographic order) in [0,1].
    const Vector &GetControlPoints() const { return control_points; }
 
+   /** @brief Compute element-wise bounds from a lexicographic E-vector.
+    *
+    *  @details The expected layout of @a e_vec is `ND x VDIM x NE`, where
+    *  `ND = nb^rdim`, `VDIM = fes_vdim`, and `NE` is the number of elements.
+    *  The output layout matches GridFunction::GetElementBounds:
+    *  `NE x active_vdim`, with the element index varying fastest.
+    */
+   void GetElementBoundsKernel(const int rdim, const int fes_vdim,
+                               const Vector &e_vec, Vector &lower,
+                               Vector &upper, const int vdim = 0) const;
+
    /** @brief Get lower and upper bounding matrix (ncp^dim x nb^dim)
     *
     *  @details The matrices can be used to compute the bounds at control points
@@ -183,6 +203,600 @@ private:
               const int cp_type_i, const real_t tol_i);
 };
 
+namespace internal
+{
+
+struct PLBoundDeviceData
+{
+   int nb;
+   int ncp;
+   const real_t *xhat;
+   const real_t *what;
+   const real_t *cphat;
+   const real_t *lbound;
+   const real_t *ubound;
+};
+
+template<int T_NB = 0, bool T_PROJ = true>
+inline void GetElementBoundsKernel1D(const PLBoundDeviceData &data,
+                                     const int fes_vdim,
+                                     const int ne,
+                                     const Vector &e_vec,
+                                     Vector &lower,
+                                     Vector &upper,
+                                     const int comp0,
+                                     const int ncomp)
+{
+   constexpr int GENERIC_MAX_ND = 32;
+   constexpr int MAX_ND = T_NB ? T_NB : GENERIC_MAX_ND;
+   constexpr int BLOCK_X = 2*MAX_ND;
+
+   const int nd = T_NB ? T_NB : data.nb;
+   MFEM_VERIFY(nd <= MAX_ND,
+               "Device element bounds kernel supports up to 32 "
+               "1D degrees of freedom.");
+
+   const auto E = Reshape(e_vec.Read(), nd, fes_vdim, ne);
+   auto L = Reshape(lower.Write(), ne, ncomp);
+   auto U = Reshape(upper.Write(), ne, ncomp);
+
+   mfem::forall_2D<BLOCK_X>(ne*ncomp, BLOCK_X, 1,
+                            [=] MFEM_HOST_DEVICE (int ec)
+   {
+      const int e = ec % ne;
+      const int c = ec / ne;
+      const int vc = comp0 + c;
+      const real_t *coeff = &E(0, vc, e);
+      const int tid = MFEM_THREAD_ID(x);
+
+      MFEM_SHARED real_t sproj[MAX_ND];
+      MFEM_SHARED real_t ssum0[MAX_ND];
+      MFEM_SHARED real_t ssum1[MAX_ND];
+      MFEM_SHARED real_t smin[BLOCK_X];
+      MFEM_SHARED real_t smax[BLOCK_X];
+      MFEM_SHARED real_t sa0;
+      MFEM_SHARED real_t sa1;
+
+      MFEM_FOREACH_THREAD(i, x, nd)
+      {
+         if constexpr (T_PROJ)
+         {
+            const real_t x = data.xhat[i];
+            const real_t w = data.what[i];
+            ssum0[i] = 0.5*coeff[i]*w;
+            ssum1[i] = 1.5*coeff[i]*w*x;
+         }
+         else
+         {
+            ssum0[i] = 0.0;
+            ssum1[i] = 0.0;
+         }
+      }
+      MFEM_SYNC_THREAD;
+
+      MFEM_FOREACH_THREAD(ii, x, 1)
+      {
+         sa0 = 0.0;
+         sa1 = 0.0;
+         for (int i = 0; i < nd; i++)
+         {
+            sa0 += ssum0[i];
+            sa1 += ssum1[i];
+         }
+      }
+      MFEM_SYNC_THREAD;
+
+      MFEM_FOREACH_THREAD(i, x, nd)
+      {
+         if constexpr (T_PROJ)
+         {
+            const real_t x = data.xhat[i];
+            sproj[i] = coeff[i] - sa0 - sa1*x;
+         }
+         else
+         {
+            sproj[i] = coeff[i];
+         }
+      }
+      MFEM_SYNC_THREAD;
+
+      real_t lower_local = HUGE_VAL;
+      real_t upper_local = -HUGE_VAL;
+      MFEM_FOREACH_THREAD(j, x, data.ncp)
+      {
+         real_t lo = 0.0;
+         real_t hi = 0.0;
+         if constexpr (T_PROJ)
+         {
+            const real_t xcp = data.cphat[j];
+            lo = sa0 + sa1*xcp;
+            hi = lo;
+         }
+
+         for (int i = 0; i < nd; i++)
+         {
+            const real_t val = sproj[i];
+            const real_t lv = data.lbound[j + i*data.ncp]*val;
+            const real_t uv = data.ubound[j + i*data.ncp]*val;
+            lo += lv < uv ? lv : uv;
+            hi += lv > uv ? lv : uv;
+         }
+         lower_local = lower_local < lo ? lower_local : lo;
+         upper_local = upper_local > hi ? upper_local : hi;
+      }
+
+      smin[tid] = lower_local;
+      smax[tid] = upper_local;
+      MFEM_SYNC_THREAD;
+
+      MFEM_FOREACH_THREAD(ii, x, 1)
+      {
+         real_t lower_ec = smin[0];
+         real_t upper_ec = smax[0];
+         const int nthreads = MFEM_THREAD_SIZE(x);
+         const int nactive = data.ncp < nthreads ? data.ncp : nthreads;
+         for (int t = 1; t < nactive; t++)
+         {
+            lower_ec = lower_ec < smin[t] ? lower_ec : smin[t];
+            upper_ec = upper_ec > smax[t] ? upper_ec : smax[t];
+         }
+         L(e, c) = lower_ec;
+         U(e, c) = upper_ec;
+      }
+   });
+}
+
+template<int T_NB = 0, int T_NCP = 0, bool T_PROJ = true>
+inline void GetElementBoundsKernel2D(const PLBoundDeviceData &data,
+                                     const int fes_vdim,
+                                     const int ne,
+                                     const Vector &e_vec,
+                                     Vector &lower,
+                                     Vector &upper,
+                                     const int comp0,
+                                     const int ncomp)
+{
+   constexpr int DEFAULT_MAX_NB = 8;
+   constexpr int DEFAULT_MAX_CP = 3*DEFAULT_MAX_NB;
+   constexpr int MAX_NB = T_NB ? T_NB : DEFAULT_MAX_NB;
+   constexpr int MAX_CP = T_NCP ? T_NCP : DEFAULT_MAX_CP;
+   constexpr int MAX_THREADS = MAX_CP*MAX_CP;
+
+   const int nb = data.nb;
+   const int ncp = data.ncp;
+   const int nd = nb*nb;
+   MFEM_VERIFY(nb <= MAX_NB,
+               "Device 2D element bounds kernel exceeds its compile-time "
+               "1D degree bound.");
+   MFEM_VERIFY(ncp <= MAX_CP,
+               "Device 2D element bounds kernel exceeds its compile-time "
+               "control-point bound.");
+   MFEM_VERIFY(ncp*ncp <= MAX_THREADS,
+               "Device 2D element bounds kernel exceeds its compile-time "
+               "thread-block bound.");
+
+   const auto E = Reshape(e_vec.Read(), nd, fes_vdim, ne);
+   auto L = Reshape(lower.Write(), ne, ncomp);
+   auto U = Reshape(upper.Write(), ne, ncomp);
+
+   mfem::forall_2D<MAX_THREADS>(ne*ncomp, ncp, ncp,
+                                [=] MFEM_HOST_DEVICE (int ec)
+   {
+      const int e = ec % ne;
+      const int c = ec / ne;
+      const int vc = comp0 + c;
+      const real_t *coeff = &E(0, vc, e);
+      const int tx = MFEM_THREAD_ID(x);
+      const int ty = MFEM_THREAD_ID(y);
+
+      MFEM_SHARED real_t sproj[MAX_NB*MAX_NB];
+      MFEM_SHARED real_t srow_min[MAX_NB*MAX_CP];
+      MFEM_SHARED real_t srow_max[MAX_NB*MAX_CP];
+      MFEM_SHARED real_t srow_a0[MAX_NB];
+      MFEM_SHARED real_t srow_a1[MAX_NB];
+      MFEM_SHARED real_t sa0[MAX_CP];
+      MFEM_SHARED real_t sa1[MAX_CP];
+      MFEM_SHARED real_t smin[MAX_THREADS];
+      MFEM_SHARED real_t smax[MAX_THREADS];
+
+      // Stage 1a: for each nodal row, form the per-node contributions to the
+      // row-wise linear fit used by the first 1D bounding solve.
+      MFEM_FOREACH_THREAD(jrow, y, nb)
+      {
+         const real_t *row_coeff = coeff + jrow*nb;
+         const int row_ncp_off = jrow*MAX_CP;
+         MFEM_FOREACH_THREAD(i, x, nb)
+         {
+            if constexpr (T_PROJ)
+            {
+               const real_t x = data.xhat[i];
+               const real_t w = data.what[i];
+               srow_min[row_ncp_off + i] = 0.5*row_coeff[i]*w;
+               srow_max[row_ncp_off + i] = 1.5*row_coeff[i]*w*x;
+            }
+            else
+            {
+               srow_min[row_ncp_off + i] = 0.0;
+               srow_max[row_ncp_off + i] = 0.0;
+            }
+         }
+      }
+      MFEM_SYNC_THREAD;
+
+      // Stage 1b: reduce the row-wise projection coefficients a0/a1.
+      if constexpr (T_PROJ)
+      {
+         MFEM_FOREACH_THREAD(jrow, y, nb)
+         {
+            const int row_ncp_off = jrow*MAX_CP;
+            real_t a0 = 0.0;
+            real_t a1 = 0.0;
+            MFEM_FOREACH_THREAD(ii, x, 1)
+            {
+               for (int i = 0; i < nb; i++)
+               {
+                  a0 += srow_min[row_ncp_off + i];
+                  a1 += srow_max[row_ncp_off + i];
+               }
+               srow_a0[jrow] = a0;
+               srow_a1[jrow] = a1;
+            }
+         }
+         MFEM_SYNC_THREAD;
+      }
+
+      // Stage 1c: subtract the row-wise linear fit once and cache the
+      // projected row coefficients for reuse across all x-control points.
+      MFEM_FOREACH_THREAD(jrow, y, nb)
+      {
+         const real_t *row_coeff = coeff + jrow*nb;
+         MFEM_FOREACH_THREAD(i, x, nb)
+         {
+            if constexpr (T_PROJ)
+            {
+               const real_t x = data.xhat[i];
+               sproj[jrow*MAX_NB + i] = row_coeff[i]
+                                        - srow_a0[jrow] - srow_a1[jrow]*x;
+            }
+            else
+            {
+               sproj[jrow*MAX_NB + i] = row_coeff[i];
+            }
+         }
+      }
+      MFEM_SYNC_THREAD;
+
+      // Stage 1d: solve the first 1D bounding problem along each nodal row and
+      // store bounds at every x-direction control point.
+      MFEM_FOREACH_THREAD(icp, x, ncp)
+      {
+         MFEM_FOREACH_THREAD(jrow, y, nb)
+         {
+            const int row_cp_off = jrow*ncp;
+            real_t lo = 0.0;
+            real_t hi = 0.0;
+            if constexpr (T_PROJ)
+            {
+               const real_t xcp = data.cphat[icp];
+               lo = srow_a0[jrow] + srow_a1[jrow]*xcp;
+               hi = lo;
+            }
+            for (int i = 0; i < nb; i++)
+            {
+               const real_t val = sproj[jrow*MAX_NB + i];
+               const real_t lv = data.lbound[icp + i*data.ncp]*val;
+               const real_t uv = data.ubound[icp + i*data.ncp]*val;
+               lo += lv < uv ? lv : uv;
+               hi += lv > uv ? lv : uv;
+            }
+            srow_min[row_cp_off + icp] = lo;
+            srow_max[row_cp_off + icp] = hi;
+         }
+      }
+      MFEM_SYNC_THREAD;
+
+      // Stage 2a: from the row bounds, form the per-row contributions to the
+      // second 1D projection solve in the y-direction.
+      MFEM_FOREACH_THREAD(icp, x, ncp)
+      {
+         MFEM_FOREACH_THREAD(jrow, y, nb)
+         {
+            const int row_cp_off = jrow*ncp;
+            if constexpr (T_PROJ)
+            {
+               const real_t x = data.xhat[jrow];
+               const real_t w = data.what[jrow];
+               const real_t t = 0.5*(srow_min[row_cp_off + icp] +
+                                     srow_max[row_cp_off + icp]);
+               smin[row_cp_off + icp] = 0.5*t*w;
+               smax[row_cp_off + icp] = 1.5*t*w*x;
+            }
+            else
+            {
+               smin[row_cp_off + icp] = 0.0;
+               smax[row_cp_off + icp] = 0.0;
+            }
+         }
+      }
+      MFEM_SYNC_THREAD;
+
+      // Stage 2b: reduce the y-direction projection coefficients for each
+      // x-control-point column.
+      MFEM_FOREACH_THREAD(jj, y, 1)
+      {
+         MFEM_FOREACH_THREAD(icp, x, ncp)
+         {
+            real_t a0 = 0.0;
+            real_t a1 = 0.0;
+            for (int jrow = 0; jrow < nb; jrow++)
+            {
+               a0 += smin[jrow*ncp + icp];
+               a1 += smax[jrow*ncp + icp];
+            }
+            sa0[icp] = a0;
+            sa1[icp] = a1;
+         }
+      }
+      MFEM_SYNC_THREAD;
+
+      // Stage 2c: subtract the y-direction linear fit from the intermediate
+      // row bounds so the final tensor-product bound uses the perturbation.
+      if constexpr (T_PROJ)
+      {
+         MFEM_FOREACH_THREAD(icp, x, ncp)
+         {
+            MFEM_FOREACH_THREAD(jrow, y, nb)
+            {
+               const int row_cp_off = jrow*ncp;
+               const real_t x = data.xhat[jrow];
+               const real_t t = sa0[icp] + sa1[icp]*x;
+               srow_min[row_cp_off + icp] -= t;
+               srow_max[row_cp_off + icp] -= t;
+            }
+         }
+      }
+      MFEM_SYNC_THREAD;
+
+      // Stage 3: each thread now owns one 2D control point (icp, kcp) and
+      // accumulates its final lower/upper bound from the row-bound data.
+      MFEM_FOREACH_THREAD(icp, x, ncp)
+      {
+         MFEM_FOREACH_THREAD(kcp, y, ncp)
+         {
+            real_t lo = 0.0;
+            real_t hi = 0.0;
+            if constexpr (T_PROJ)
+            {
+               const real_t xcp = data.cphat[kcp];
+               lo = sa0[icp] + sa1[icp]*xcp;
+               hi = lo;
+            }
+            for (int jrow = 0; jrow < nb; jrow++)
+            {
+               const real_t w0 = srow_min[jrow*ncp + icp];
+               const real_t w1 = srow_max[jrow*ncp + icp];
+               const real_t lb = data.lbound[kcp + jrow*data.ncp];
+               const real_t ub = data.ubound[kcp + jrow*data.ncp];
+               const real_t v0 = lb*w0;
+               const real_t v1 = ub*w0;
+               const real_t v2 = lb*w1;
+               const real_t v3 = ub*w1;
+               real_t vlo = v0 < v1 ? v0 : v1;
+               real_t vhi = v0 > v1 ? v0 : v1;
+               vlo = vlo < v2 ? vlo : v2;
+               vlo = vlo < v3 ? vlo : v3;
+               vhi = vhi > v2 ? vhi : v2;
+               vhi = vhi > v3 ? vhi : v3;
+               lo += vlo;
+               hi += vhi;
+            }
+            const int slot = kcp*ncp + icp;
+            smin[slot] = lo;
+            smax[slot] = hi;
+         }
+      }
+      MFEM_SYNC_THREAD;
+
+      const int lane = ty*ncp + tx;
+      const int nactive = ncp*ncp;
+      const int nthreads = MFEM_THREAD_SIZE(x)*MFEM_THREAD_SIZE(y);
+
+      // Reduce all 2D control-point bounds to one lower/upper pair per
+      // (element, component).
+      if (nthreads == 1)
+      {
+         if (tx == 0 && ty == 0)
+         {
+            real_t lower_ec = smin[0];
+            real_t upper_ec = smax[0];
+            for (int t = 1; t < nactive; t++)
+            {
+               lower_ec = lower_ec < smin[t] ? lower_ec : smin[t];
+               upper_ec = upper_ec > smax[t] ? upper_ec : smax[t];
+            }
+            L(e, c) = lower_ec;
+            U(e, c) = upper_ec;
+         }
+      }
+      else
+      {
+         for (int stride = (nactive + 1)/2; stride > 0;
+              stride = (stride + 1)/2)
+         {
+            if (lane < stride && lane + stride < nactive)
+            {
+               smin[lane] = smin[lane] < smin[lane + stride] ?
+                            smin[lane] : smin[lane + stride];
+               smax[lane] = smax[lane] > smax[lane + stride] ?
+                            smax[lane] : smax[lane + stride];
+            }
+            MFEM_SYNC_THREAD;
+            if (stride == 1) { break; }
+         }
+
+         if (lane == 0)
+         {
+            L(e, c) = smin[0];
+            U(e, c) = smax[0];
+         }
+      }
+   });
+}
+
+} // namespace internal
+
+inline void PLBound::GetElementBoundsKernel(const int rdim, const int fes_vdim,
+                                            const Vector &e_vec,
+                                            Vector &lower, Vector &upper,
+                                            const int vdim) const
+{
+   MFEM_VERIFY(b_type != BasisType::Positive,
+               "Bernstein device bounds are not implemented.");
+   if (rdim == 3)
+   {
+      MFEM_ABORT("Device element bounds kernel currently only supports 1D/2D.");
+   }
+   MFEM_VERIFY(rdim == 1 || rdim == 2, "Invalid element dimension.");
+   MFEM_VERIFY(vdim >= -1 && vdim <= fes_vdim, "Invalid vector component.");
+   const int nd = static_cast<int>(std::pow(nb, rdim));
+   const int ne = e_vec.Size()/(nd*fes_vdim);
+   const int ncomp = (vdim > 0) ? 1 : fes_vdim;
+
+   lower.SetSize(ne*ncomp, e_vec);
+   upper.SetSize(ne*ncomp, e_vec);
+   lower.UseDevice(true);
+   upper.UseDevice(true);
+
+   if (!proj)
+   {
+      MFEM_ABORT("Device element bounds kernel currently requires projection "
+                 "enabled.");
+   }
+
+   const real_t *dxhat = xhat.Read();
+   const real_t *dwhat = what.Read();
+   const real_t *dcphat = cphat.Read();
+   const real_t *dlbound = lbound.Read();
+   const real_t *dubound = ubound.Read();
+
+   internal::PLBoundDeviceData data
+   {
+      nb,
+      ncp,
+      dxhat,
+      dwhat,
+      dcphat,
+      dlbound,
+      dubound
+   };
+
+   const int comp0 = (vdim > 0) ? (vdim - 1) : 0;
+
+   if (rdim == 1)
+   {
+      switch (nb)
+      {
+         case 2: return internal::GetElementBoundsKernel1D<2, true>(data, fes_vdim, ne,
+                                                                     e_vec, lower, upper,
+                                                                     comp0, ncomp);
+         case 3: return internal::GetElementBoundsKernel1D<3, true>(data, fes_vdim, ne,
+                                                                     e_vec, lower, upper,
+                                                                     comp0, ncomp);
+         case 4: return internal::GetElementBoundsKernel1D<4, true>(data, fes_vdim, ne,
+                                                                     e_vec, lower, upper,
+                                                                     comp0, ncomp);
+         case 5: return internal::GetElementBoundsKernel1D<5, true>(data, fes_vdim, ne,
+                                                                     e_vec, lower, upper,
+                                                                     comp0, ncomp);
+         case 6: return internal::GetElementBoundsKernel1D<6, true>(data, fes_vdim, ne,
+                                                                     e_vec, lower, upper,
+                                                                     comp0, ncomp);
+         case 7: return internal::GetElementBoundsKernel1D<7, true>(data, fes_vdim, ne,
+                                                                     e_vec, lower, upper,
+                                                                     comp0, ncomp);
+         case 8: return internal::GetElementBoundsKernel1D<8, true>(data, fes_vdim, ne,
+                                                                     e_vec, lower, upper,
+                                                                     comp0, ncomp);
+         case 9: return internal::GetElementBoundsKernel1D<9, true>(data, fes_vdim, ne,
+                                                                     e_vec, lower, upper,
+                                                                     comp0, ncomp);
+         case 10: return internal::GetElementBoundsKernel1D<10, true>(data, fes_vdim, ne,
+                                                                       e_vec, lower, upper,
+                                                                       comp0, ncomp);
+         default: return internal::GetElementBoundsKernel1D<0, true>(data, fes_vdim, ne,
+                                                                      e_vec, lower, upper,
+                                                                      comp0, ncomp);
+      }
+   }
+#define MFEM_PLBOUND_2D_DISPATCH(NB, NCP) \
+   return internal::GetElementBoundsKernel2D<NB, NCP, true>(data, fes_vdim, ne, \
+                                                            e_vec, lower, upper, \
+                                                            comp0, ncomp)
+   switch (nb)
+   {
+      case 2:
+         switch (ncp)
+         {
+            case 4: MFEM_PLBOUND_2D_DISPATCH(2, 4);
+            case 6: MFEM_PLBOUND_2D_DISPATCH(2, 6);
+            case 8: MFEM_PLBOUND_2D_DISPATCH(2, 8);
+         }
+         break;
+      case 3:
+         switch (ncp)
+         {
+            case 6: MFEM_PLBOUND_2D_DISPATCH(3, 6);
+            case 9: MFEM_PLBOUND_2D_DISPATCH(3, 9);
+            case 12: MFEM_PLBOUND_2D_DISPATCH(3, 12);
+         }
+         break;
+      case 4:
+         switch (ncp)
+         {
+            case 8: MFEM_PLBOUND_2D_DISPATCH(4, 8);
+            case 12: MFEM_PLBOUND_2D_DISPATCH(4, 12);
+            case 16: MFEM_PLBOUND_2D_DISPATCH(4, 16);
+         }
+         break;
+      case 5:
+         switch (ncp)
+         {
+            case 10: MFEM_PLBOUND_2D_DISPATCH(5, 10);
+            case 15: MFEM_PLBOUND_2D_DISPATCH(5, 15);
+            case 20: MFEM_PLBOUND_2D_DISPATCH(5, 20);
+         }
+         break;
+      case 6:
+         switch (ncp)
+         {
+            case 12: MFEM_PLBOUND_2D_DISPATCH(6, 12);
+            case 18: MFEM_PLBOUND_2D_DISPATCH(6, 18);
+            case 24: MFEM_PLBOUND_2D_DISPATCH(6, 24);
+         }
+         break;
+      case 7:
+         switch (ncp)
+         {
+            case 14: MFEM_PLBOUND_2D_DISPATCH(7, 14);
+            case 21: MFEM_PLBOUND_2D_DISPATCH(7, 21);
+            case 28: MFEM_PLBOUND_2D_DISPATCH(7, 28);
+         }
+         break;
+      case 8:
+         switch (ncp)
+         {
+            case 16: MFEM_PLBOUND_2D_DISPATCH(8, 16);
+            case 24: MFEM_PLBOUND_2D_DISPATCH(8, 24);
+            case 32: MFEM_PLBOUND_2D_DISPATCH(8, 32);
+         }
+         break;
+   }
+#undef MFEM_PLBOUND_2D_DISPATCH
+   return internal::GetElementBoundsKernel2D<0, 0, true>(data, fes_vdim, ne,
+                                                         e_vec, lower, upper,
+                                                         comp0, ncomp);
+}
+
 } // namespace mfem
 
 #endif // MFEM_BOUNDS

fem/coefficient.hpp

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

fem/dgmassinv_kernels.hpp

@@ -57,7 +57,7 @@ void DGMassApply(const int e,
       }
       else if (DIM == 3)
       {
-         SmemPAMassApply3D_Element<TD1D,TQ1D,ACCUM>(e, NE, B, pa_data, x, y);
+         SmemPAMassApply3D_Element<TD1D,TQ1D,NBZ,ACCUM>(e, NE, B, pa_data, x, y);
       }
       else
       {

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 methods
-    ComputeElementFlux() and ComputeFluxEnergy().
+    The required BilinearFormIntegrator must implement the method
+    ComputeElementFlux().
 
     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 methods
-                       ComputeElementFlux() and ComputeFluxEnergy().
+       @param integ    This BilinearFormIntegrator must implement the method
+                       ComputeElementFlux().
        @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 methods
-                       ComputeElementFlux() and ComputeFluxEnergy().
+       @param integ    This BilinearFormIntegrator must implement the method
+                       ComputeElementFlux().
        @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,7 +20,6 @@
 #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). */
-   int GetPhysRangeDim(int /* space_dim */) const { return vdim; }
+   virtual 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).
    */
-   int GetPhysCurlDim(int /* space_dim */) const { return cdim; }
+   virtual 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 { return space_dim; }
+   int GetPhysRangeDim(int space_dim) const override { return space_dim; }
 };
 
 /// @brief Class for computing 1D special polynomials and their associated basis

fem/fe/fe_h1_bubble.cpp

@@ -1,973 +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.
-
-// 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

@@ -1,109 +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_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 { return 2; }
-   int GetPhysCurlDim(int space_dim) const { return 1; }
+   int GetPhysRangeDim(int space_dim) const override { return 2; }
+   int GetPhysCurlDim(int space_dim) const override { 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 { return 3; }
-   int GetPhysCurlDim(int space_dim) const { return 3; }
+   int GetPhysRangeDim(int space_dim) const override { return 3; }
+   int GetPhysCurlDim(int space_dim) const override { 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 { return 2; }
-   int GetPhysCurlDim(int space_dim) const { return 0; }
+   int GetPhysRangeDim(int space_dim) const override { return 2; }
+   int GetPhysCurlDim(int space_dim) const override { 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 { return 3; }
-   int GetPhysCurlDim(int space_dim) const { return 0; }
+   int GetPhysRangeDim(int space_dim) const override { return 3; }
+   int GetPhysCurlDim(int space_dim) const override { return 0; }
 
    using FiniteElement::CalcVShape;
 

fem/fe_coll.cpp

@@ -243,21 +243,11 @@ 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));
@@ -2132,173 +2122,6 @@ 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,8 +111,6 @@ 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) |
@@ -319,59 +317,6 @@ 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,12 +3877,9 @@ 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). 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)
+      // consistency check: fec->GetOrder() and FE->GetOrder() should return
+      // the same value (for standard, constant-order spaces)
+      if (!IsVariableOrder() && FE->GetDim() > 0)
       {
          MFEM_ASSERT(FE->GetOrder() == fec->GetOrder(),
                      "internal error: " <<
@@ -3937,6 +3934,16 @@ 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");
@@ -3967,6 +3974,11 @@ 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 fucntion would give 15, not 10.
+       Then this function would give 15, not 10.
        */
    int GetVectorDim() const;
 
@@ -1323,12 +1323,24 @@ 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,27 +345,6 @@ 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();
@@ -2050,6 +2029,18 @@ 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;
 
@@ -2202,6 +2193,9 @@ 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;
@@ -2355,6 +2349,9 @@ 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;
@@ -2630,6 +2627,7 @@ 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;
@@ -2945,6 +2943,7 @@ 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;
@@ -2974,6 +2973,7 @@ 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,9 +3030,14 @@ void GridFunction::ProjectCoefficient(Coefficient *coeff[])
    }
 }
 
-void GridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff,
-                                          Array<int> &dof_attr)
+void GridFunction::ProjectDiscCoefficient(
+   std::variant<Coefficient*, VectorCoefficient*> coeff, Array<int> &dof_attr)
 {
+   std::visit([&](auto* c)
+   {
+      MFEM_VERIFY(VectorDim() == c->GetVDim(), "coeff vdim != VectorDim()");
+   }, coeff);
+
    Array<int> vdofs;
    Vector vals;
 
@@ -3046,7 +3051,10 @@ void GridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff,
    {
       fes->GetElementVDofs(i, vdofs);
       vals.SetSize(vdofs.Size());
-      fes->GetFE(i)->Project(coeff, *fes->GetElementTransformation(i), vals);
+      std::visit([&](auto* c)
+      {
+         fes->GetFE(i)->Project(*c, *fes->GetElementTransformation(i), vals);
+      }, coeff);
 
       // the values in shared dofs are determined from the element with maximal
       // attribute
@@ -3062,17 +3070,15 @@ void GridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff,
    }
 }
 
-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);
 
@@ -3082,6 +3088,7 @@ 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);
 
@@ -3139,12 +3146,16 @@ void GridFunction::ProjectBdrCoefficient(Coefficient *coeff[],
 void GridFunction::ProjectBdrCoefficientNormal(
    Coefficient *coeff, VectorCoefficient *vcoeff, const Array<int> &bdr_attr)
 {
-   if (fes->GetNBE() > 0)
+   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)
    {
-      // 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!");
+      MFEM_VERIFY(vcoeff->GetVDim() == fes->GetMesh()->SpaceDimension(),
+                  "vcoeff vdim (" << vcoeff->GetVDim()
+                  << ") != SpaceDimension ("
+                  << fes->GetMesh()->SpaceDimension() << ")");
    }
 
    // implementation for the case when the face dofs are scaled point
@@ -5241,6 +5252,30 @@ void GridFunction::GetElementBounds(const PLBound &plb,
                                     Vector &lower, Vector &upper,
                                     const int vdim) const
 {
+   if (UseDevice() && Device::Allows(Backend::DEVICE_MASK) &&
+       plb.GetBasisType() != BasisType::Positive &&
+       UsesTensorBasis(*fes))
+   {
+      const FiniteElement &fe = *fes->GetTypicalFE();
+      const int rdim = fe.GetDim();
+      const int fes_dim = fes->GetVDim();
+      const int nel = fes->GetNE();
+      const int nd = fe.GetDof();
+
+      Vector e_vec(nd*fes_dim*nel, Device::GetDeviceMemoryType());
+      e_vec.UseDevice(true);
+      const ElementRestrictionOperator *elem_restr =
+         fes->GetElementRestriction(ElementDofOrdering::LEXICOGRAPHIC);
+      MFEM_VERIFY(elem_restr != nullptr,
+                  "Element restriction is required for device bounds.");
+      elem_restr->Mult(*this, e_vec);
+
+      plb.GetElementBoundsKernel(rdim, fes_dim, e_vec, lower, upper, vdim);
+      lower.HostRead();
+      upper.HostRead();
+      return;
+   }
+
    int nel  = fes->GetNE();
    int fes_dim  = fes->GetVDim();
    lower.SetSize(nel*(vdim > 0 ? 1 :fes_dim));
@@ -5757,4 +5792,4 @@ std::pair<real_t, real_t> GridFunction::EstimateFunctionMaximum(
    return std::make_pair(global_max_lower, global_max_upper);
 }
 
-}
+}

fem/gridfunc.hpp

@@ -23,6 +23,7 @@
 #include <limits>
 #include <ostream>
 #include <string>
+#include <variant>
 
 namespace mfem
 {
@@ -79,10 +80,18 @@ protected:
                         bool wcoef,
                         int subdomain);
 
-   /** 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);
+   /** @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); };
 
    /** Helper function for ProjectCoefficientElementL2 */
    void ProjectCoefficientElementL2_(Coefficient &coeff, Vector &sol, Vector &Va);
@@ -150,11 +159,13 @@ public:
 
    FiniteElementCollection *OwnFEC() { return fec_owned; }
 
-   /// Shortcut for calling FiniteElementSpace::GetVectorDim() on the underlying #fes
-   int VectorDim() const;
+   /** @brief Shortcut for calling FiniteElementSpace::GetVectorDim() on the
+       underlying #fes */
+   int VectorDim() const { return fes->GetVectorDim(); }
 
-   /// Shortcut for calling FiniteElementSpace::GetCurlDim() on the underlying #fes
-   int CurlDim() const;
+   /** @brief Shortcut for calling FiniteElementSpace::GetCurlDim() on the
+       underlying #fes */
+   int CurlDim() const { return fes->GetCurlDim(); }
 
    /// Read only access to the (optional) internal true-dof Vector.
    const Vector &GetTrueVector() const
@@ -513,10 +524,17 @@ 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);
+   virtual void ProjectDiscCoefficient(VectorCoefficient &coeff)
+   { ProjectDiscCoefficient(&coeff); }
 
    enum AvgType {ARITHMETIC, HARMONIC};
    /** @brief Projects a discontinuous coefficient so that the values in shared
@@ -1971,6 +1989,7 @@ public:
 
    void Eval(Vector &v, ElementTransformation &T,
              const IntegrationPoint &ip) override;
+   using VectorCoefficient::Eval;
 
    virtual ~VectorExtrudeCoefficient() { }
 };

fem/gslib.cpp

@@ -106,9 +106,7 @@ FindPointsGSLIB::FindPointsGSLIB()
    gsl_comm = new gslib::comm;
    cr       = new gslib::crystal;
 #ifdef MFEM_USE_MPI
-   int initialized = 0;
-   MPI_Initialized(&initialized);
-   if (!initialized) { MPI_Init(NULL, NULL); }
+   if (!Mpi::IsInitialized()) { Mpi::Init(); }
    MPI_Comm comm = MPI_COMM_WORLD;
    comm_init(gsl_comm, comm);
 #else
@@ -490,7 +488,7 @@ void FindPointsGSLIB::FindPointsOnDevice(const Vector &point_pos,
    }
    DEV.find_device = true;
 
-   const int id = gsl_comm->id, np = gsl_comm->np;
+   const unsigned int id = gsl_comm->id, np = gsl_comm->np;
 
    gsl_mfem_ref.SetSize(points_cnt * dim);
    gsl_mfem_elem.SetSize(points_cnt);
@@ -652,7 +650,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 (pp == p->proc)
+            if (static_cast<unsigned>(pp) == p->proc)
             {
                continue;
             }
@@ -1068,7 +1066,7 @@ void FindPointsGSLIB::InterpolateOnDevice(const Vector &field_in_evec,
          sarray_transfer(struct evalOutPt_t, &outpt, proc, 1, cr);
 
          opt = (evalOutPt_t *)outpt.ptr;
-         for (int index = 0; index < outpt.n; index++)
+         for (size_t index = 0; index < outpt.n; index++)
          {
             int idx = ordering == Ordering::byNODES ?
                       opt->index + i*points_cnt :
@@ -1413,7 +1411,7 @@ void FindPointsGSLIB::SetupSplitMeshesAndIntegrationRules(const int order)
 {
    MFEM_VERIFY(mesh, "Setup FindPointsGSLIB with mesh first.");
    const int dof1D = order+1;
-   const int dim = mesh->Dimension();
+   dim = mesh->Dimension();
 
    SetupSplitMeshes();
    if (dim == 2)
@@ -2254,7 +2252,8 @@ void FindPointsGSLIB::DistributeInterpolatedValues(const Vector &int_vals,
       sarray_transfer(struct out_pt, outpt, proc, 1, cr);
 
       // Store received data
-      MFEM_VERIFY(outpt->n == points_cnt, "Incompatible size. Number of points "
+      MFEM_VERIFY(outpt->n == static_cast<size_t>(points_cnt),
+                  "Incompatible size. Number of points "
                   "received does not match the number of points originally "
                   "found using FindPoints.");
 
@@ -2623,9 +2622,7 @@ GSOPGSLIB::GSOPGSLIB(Array<long long> &ids)
    gsl_comm = new gslib::comm;
    cr       = new gslib::crystal;
 #ifdef MFEM_USE_MPI
-   int initialized;
-   MPI_Initialized(&initialized);
-   if (!initialized) { MPI_Init(NULL, NULL); }
+   if (!Mpi::IsInitialized()) { Mpi::Init(); }
    MPI_Comm comm = MPI_COMM_WORLD;
    comm_init(gsl_comm, comm);
 #else

fem/gslib.hpp

@@ -202,13 +202,19 @@ 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 (side_init != (1u << ei))
+   if (static_cast<unsigned>(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 (side_init != (1u << fi))
+   if (static_cast<unsigned>(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 (side_init != (64u << ei))
+   if (static_cast<unsigned>(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/integ/bilininteg_mass_kernels.hpp

@@ -181,6 +181,12 @@ constexpr int NBZ(int D1D)
 {
    return ipow(2, D(D1D) >= 0 ? D(D1D) : 0);
 }
+constexpr int NBZ3D(int MDQ)
+{
+   return MDQ > 0 ? std::min<int>(
+             (128 + MDQ * MDQ * MDQ - 1) / (MDQ * MDQ * MDQ), 64)
+          : 1;
+}
 }
 
 // Shared memory PA Mass Diagonal 2D kernel
@@ -804,19 +810,23 @@ void PAMassApply3D_Element(const int e,
    }
 }
 
-template<int T_D1D, int T_Q1D, bool ACCUMULATE = true>
-MFEM_HOST_DEVICE inline
-void SmemPAMassApply3D_Element(const int e,
-                               const int NE,
-                               const real_t *b_,
-                               const real_t *d_,
-                               const real_t *x_,
-                               real_t *y_,
-                               const int d1d = 0,
-                               const int q1d = 0)
+template <int T_D1D, int T_Q1D, int TBATCH, bool ACCUMULATE = true>
+MFEM_HOST_DEVICE inline void
+SmemPAMassApply3D_Element(const int e, const int NE, const real_t *b_,
+                          const real_t *d_, const real_t *x_, real_t *y_,
+                          int d1d = 0, int q1d = 0)
 {
-   constexpr int D1D = T_D1D ? T_D1D : d1d;
-   constexpr int Q1D = T_Q1D ? T_Q1D : q1d;
+   static_assert(TBATCH > 0, "TBATCH must be positive");
+#if defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__)
+   constexpr int tbatch = TBATCH;
+   const int tidz = MFEM_THREAD_ID(z);
+#else
+   // host always batch size 1
+   constexpr int tbatch = 1;
+   constexpr int tidz = 0;
+#endif
+   const int D1D = T_D1D ? T_D1D : d1d;
+   const int Q1D = T_Q1D ? T_Q1D : q1d;
    constexpr int MQ1 = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
    constexpr int MD1 = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
    constexpr int MDQ = (MQ1 > MD1) ? MQ1 : MD1;
@@ -829,33 +839,37 @@ void SmemPAMassApply3D_Element(const int e,
    MFEM_SHARED real_t sDQ[MQ1*MD1];
    real_t (*B)[MD1] = (real_t (*)[MD1]) sDQ;
    real_t (*Bt)[MQ1] = (real_t (*)[MQ1]) sDQ;
-   MFEM_SHARED real_t sm0[MDQ*MDQ*MDQ];
-   MFEM_SHARED real_t sm1[MDQ*MDQ*MDQ];
-   real_t (*X)[MD1][MD1]   = (real_t (*)[MD1][MD1]) sm0;
-   real_t (*DDQ)[MD1][MQ1] = (real_t (*)[MD1][MQ1]) sm1;
-   real_t (*DQQ)[MQ1][MQ1] = (real_t (*)[MQ1][MQ1]) sm0;
-   real_t (*QQQ)[MQ1][MQ1] = (real_t (*)[MQ1][MQ1]) sm1;
-   real_t (*QQD)[MQ1][MD1] = (real_t (*)[MQ1][MD1]) sm0;
-   real_t (*QDD)[MD1][MD1] = (real_t (*)[MD1][MD1]) sm1;
-   MFEM_FOREACH_THREAD(dy,y,D1D)
+   MFEM_SHARED real_t sm0[tbatch][MDQ*MDQ*MDQ];
+   MFEM_SHARED real_t sm1[tbatch][MDQ*MDQ*MDQ];
+   real_t (*X)[MD1][MD1]   = (real_t (*)[MD1][MD1]) (sm0+tidz);
+   real_t (*DDQ)[MD1][MQ1] = (real_t (*)[MD1][MQ1]) (sm1+tidz);
+   real_t (*DQQ)[MQ1][MQ1] = (real_t (*)[MQ1][MQ1]) (sm0+tidz);
+   real_t (*QQQ)[MQ1][MQ1] = (real_t (*)[MQ1][MQ1]) (sm1+tidz);
+   real_t (*QQD)[MQ1][MD1] = (real_t (*)[MQ1][MD1]) (sm0+tidz);
+   real_t (*QDD)[MD1][MD1] = (real_t (*)[MD1][MD1]) (sm1+tidz);
+   MFEM_FOREACH_THREAD(dy, y, D1D)
    {
-      MFEM_FOREACH_THREAD(dx,x,D1D)
+      MFEM_FOREACH_THREAD(dx, x, D1D)
       {
          MFEM_UNROLL(MD1)
          for (int dz = 0; dz < D1D; ++dz)
          {
-            X[dz][dy][dx] = x(dx,dy,dz,e);
+            X[dz][dy][dx] = x(dx, dy, dz, e);
          }
       }
-      MFEM_FOREACH_THREAD(dx,x,Q1D)
+      MFEM_FOREACH_THREAD(dx, x, Q1D) { B[dx][dy] = b(dx, dy); }
+   }
+   if (tidz == 0)
+   {
+      MFEM_FOREACH_THREAD(dy, y, D1D)
       {
-         B[dx][dy] = b(dx,dy);
+         MFEM_FOREACH_THREAD(dx, x, Q1D) { B[dx][dy] = b(dx, dy); }
       }
    }
    MFEM_SYNC_THREAD;
-   MFEM_FOREACH_THREAD(dy,y,D1D)
+   MFEM_FOREACH_THREAD(dy, y, D1D)
    {
-      MFEM_FOREACH_THREAD(qx,x,Q1D)
+      MFEM_FOREACH_THREAD(qx, x, Q1D)
       {
          real_t u[D1D];
          MFEM_UNROLL(MD1)
@@ -880,9 +894,9 @@ void SmemPAMassApply3D_Element(const int e,
       }
    }
    MFEM_SYNC_THREAD;
-   MFEM_FOREACH_THREAD(qy,y,Q1D)
+   MFEM_FOREACH_THREAD(qy, y, Q1D)
    {
-      MFEM_FOREACH_THREAD(qx,x,Q1D)
+      MFEM_FOREACH_THREAD(qx, x, Q1D)
       {
          real_t u[D1D];
          MFEM_UNROLL(MD1)
@@ -907,9 +921,9 @@ void SmemPAMassApply3D_Element(const int e,
       }
    }
    MFEM_SYNC_THREAD;
-   MFEM_FOREACH_THREAD(qy,y,Q1D)
+   MFEM_FOREACH_THREAD(qy, y, Q1D)
    {
-      MFEM_FOREACH_THREAD(qx,x,Q1D)
+      MFEM_FOREACH_THREAD(qx, x, Q1D)
       {
          real_t u[Q1D];
          MFEM_UNROLL(MQ1)
@@ -929,22 +943,22 @@ void SmemPAMassApply3D_Element(const int e,
          MFEM_UNROLL(MQ1)
          for (int qz = 0; qz < Q1D; qz++)
          {
-            QQQ[qz][qy][qx] = u[qz] * d(qx,qy,qz,e);
+            QQQ[qz][qy][qx] = u[qz] * d(qx, qy, qz, e);
          }
       }
    }
    MFEM_SYNC_THREAD;
-   MFEM_FOREACH_THREAD(di,y,D1D)
+   if (tidz == 0)
    {
-      MFEM_FOREACH_THREAD(q,x,Q1D)
+      MFEM_FOREACH_THREAD(di, y, D1D)
       {
-         Bt[di][q] = b(q,di);
+         MFEM_FOREACH_THREAD(q, x, Q1D) { Bt[di][q] = b(q, di); }
       }
    }
    MFEM_SYNC_THREAD;
-   MFEM_FOREACH_THREAD(qy,y,Q1D)
+   MFEM_FOREACH_THREAD(qy, y, Q1D)
    {
-      MFEM_FOREACH_THREAD(dx,x,D1D)
+      MFEM_FOREACH_THREAD(dx, x, D1D)
       {
          real_t u[Q1D];
          MFEM_UNROLL(MQ1)
@@ -969,9 +983,9 @@ void SmemPAMassApply3D_Element(const int e,
       }
    }
    MFEM_SYNC_THREAD;
-   MFEM_FOREACH_THREAD(dy,y,D1D)
+   MFEM_FOREACH_THREAD(dy, y, D1D)
    {
-      MFEM_FOREACH_THREAD(dx,x,D1D)
+      MFEM_FOREACH_THREAD(dx, x, D1D)
       {
          real_t u[Q1D];
          MFEM_UNROLL(MQ1)
@@ -996,9 +1010,9 @@ void SmemPAMassApply3D_Element(const int e,
       }
    }
    MFEM_SYNC_THREAD;
-   MFEM_FOREACH_THREAD(dy,y,D1D)
+   MFEM_FOREACH_THREAD(dy, y, D1D)
    {
-      MFEM_FOREACH_THREAD(dx,x,D1D)
+      MFEM_FOREACH_THREAD(dx, x, D1D)
       {
          real_t u[D1D];
          MFEM_UNROLL(MD1)
@@ -1020,11 +1034,11 @@ void SmemPAMassApply3D_Element(const int e,
          {
             if (ACCUMULATE)
             {
-               y(dx,dy,dz,e) += u[dz];
+               y(dx, dy, dz, e) += u[dz];
             }
             else
             {
-               y(dx,dy,dz,e) = u[dz];
+               y(dx, dy, dz, e) = u[dz];
             }
          }
       }
@@ -1115,8 +1129,8 @@ inline void PAMassApply3D(const int NE,
    });
 }
 
-// Shared memory PA Mass Apply 2D kernel
-template<int T_D1D = 0, int T_Q1D = 0>
+// Shared memory PA Mass Apply 3D kernel
+template<int T_D1D = 0, int T_Q1D = 0, int TBATCH=1>
 inline void SmemPAMassApply3D(const int NE,
                               const Array<real_t> &b_,
                               const Array<real_t> &bt_,
@@ -1126,6 +1140,9 @@ inline void SmemPAMassApply3D(const int NE,
                               const int d1d = 0,
                               const int q1d = 0)
 {
+   static_assert(T_D1D > 0, "T_D1D must be positive");
+   static_assert(T_Q1D > 0, "T_Q1D must be positive");
+   static_assert(TBATCH > 0, "TBATCH must be positive");
    MFEM_CONTRACT_VAR(bt_);
    const int D1D = T_D1D ? T_D1D : d1d;
    const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -1137,9 +1154,11 @@ inline void SmemPAMassApply3D(const int NE,
    const auto d = d_.Read();
    const auto x = x_.Read();
    auto y = y_.ReadWrite();
-   mfem::forall_2D<T_Q1D*T_Q1D>(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE (int e)
+   mfem::forall_2D_batch<T_Q1D * T_Q1D * TBATCH>(NE, Q1D, Q1D, TBATCH,
+                                                 [=] MFEM_HOST_DEVICE(int e)
    {
-      internal::SmemPAMassApply3D_Element<T_D1D,T_Q1D>(e, NE, b, d, x, y, d1d, q1d);
+      internal::SmemPAMassApply3D_Element<T_D1D, T_Q1D, TBATCH>(e, NE, b, d, x,
+                                                                y, d1d, q1d);
    });
 }
 
@@ -1394,7 +1413,16 @@ ApplyKernelType MassIntegrator::ApplyPAKernels::Kernel()
 {
    if constexpr (DIM == 1) { return internal::PAMassApply1D; }
    else if constexpr (DIM == 2) { return internal::SmemPAMassApply2D<T_D1D,T_Q1D>; }
-   else if constexpr (DIM == 3) { return internal::SmemPAMassApply3D<T_D1D, T_Q1D>; }
+   else if constexpr (DIM == 3)
+   {
+      constexpr int MDQ = T_D1D >= T_Q1D ? T_D1D : T_Q1D;
+      // max 64 threads in z limit in cuda and hip
+      if constexpr (MDQ > 0)
+      {
+         return internal::SmemPAMassApply3D<T_D1D, T_Q1D,
+                internal::mass::NBZ3D(MDQ)>;
+      }
+   }
    MFEM_ABORT("");
 }
 

fem/intrules.hpp

@@ -43,56 +43,52 @@ 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);
-      x = p[0];
-      if (dim > 1)
+      switch (dim)
       {
-         y = p[1];
-         if (dim > 2)
-         {
-            z = p[2];
-         }
+         case 3: Set3(p); break;
+         case 2: Set2(p); break;
+         case 1: Set1(p); break;
       }
    }
 
    void Get(real_t *p, const int dim) const
    {
       MFEM_ASSERT(1 <= dim && dim <= 3, "invalid dim: " << dim);
-      p[0] = x;
-      if (dim > 1)
+      switch (dim)
       {
-         p[1] = y;
-         if (dim > 2)
-         {
-            p[2] = z;
-         }
+         case 3: p[2] = z;
+         case 2: p[1] = y;
+         case 1: p[0] = x;
       }
    }
-
-   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)
-      : DeltaLFIntegrator(QF), Q(QF) { }
+   DomainLFGradIntegrator(VectorCoefficient &QF, const IntegrationRule *ir = NULL)
+      : DeltaLFIntegrator(QF, ir), Q(QF) { }
 
    bool SupportsDevice() const override { return true; }
 

fem/pgridfunc.cpp

@@ -545,6 +545,8 @@ 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)
@@ -715,7 +717,8 @@ void ParGridFunction::ProjectCoefficientElementL2(VectorCoefficient &vcoeff)
 }
 
 
-void ParGridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff)
+void ParGridFunction::ProjectDiscCoefficient(
+   std::variant<Coefficient*, VectorCoefficient*> coeff)
 {
    // local maximal element attribute for each dof
    Array<int> ldof_attr;
@@ -761,6 +764,9 @@ void ParGridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff)
 
 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.
 
@@ -786,6 +792,8 @@ 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);
@@ -858,6 +866,12 @@ 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,6 +63,12 @@ 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; }
 
@@ -268,11 +274,6 @@ 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;
@@ -280,8 +281,7 @@ public:
    using GridFunction::ProjectBdrCoefficient;
 
    void ProjectBdrCoefficient(VectorCoefficient &vcoeff,
-                              const Array<int> &attr) override
-   { ProjectBdrCoefficient(NULL, &vcoeff, attr); }
+                              const Array<int> &attr) override;
 
    void ProjectBdrCoefficient(Coefficient *coeff[],
                               const Array<int> &attr) override

general/arrays_by_name.hpp

@@ -14,6 +14,7 @@
 
 #include "../config/config.hpp"
 #include "array.hpp"
+#include "text.hpp"
 
 #include <iostream>
 #include <map>
@@ -247,7 +248,8 @@ inline void ArraysByName<T>::Print(std::ostream &os, int width) const
    os << data.size() << '\n';
    for (auto const &it : data)
    {
-      os << '"' << it.first << '"' << '\n' << it.second.Size() << '\n';
+      // Note: The method Load() can read any string formatted with std::quoted.
+      os << std::quoted(it.first) << '\n' << it.second.Size() << '\n';
       it.second.Print(os, width > 0 ? width : it.second.Size());
    }
 }
@@ -258,40 +260,36 @@ void ArraysByName<T>::Load(std::istream &in)
    int NumArrays;
    in >> NumArrays;
 
-   std::string ArrayLine, ArrayName;
-   for (int i=0; i < NumArrays; i++)
+   for (int i = 0; i < NumArrays; i++)
    {
       in >> std::ws;
-      getline(in, ArrayLine);
-
-      std::size_t q0 = ArrayLine.find('"');
-      std::size_t q1 = ArrayLine.rfind('"');
-
-      if (q0 != std::string::npos && q1 > q0)
+      // 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() == '"')
       {
-         // Locate set name between first and last double quote
-         ArrayName = ArrayLine.substr(q0+1,q1-q0-1);
+         if (parse_quoted_string(ArrayName, in) != 0)
+         {
+            MFEM_ABORT("error parsing input!");
+         }
       }
       else
       {
-         // 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);
+         in >> ArrayName;
+         MFEM_VERIFY(in.good(), "error parsing input!");
       }
+
+      // Read the array
+      data[ArrayName].Load(in);
    }
-
 }
 
 }

general/forall.hpp

@@ -1090,6 +1090,12 @@ inline void forall_2D_batch(int N, int X, int Y, int BZ, lambda &&body)
    ForallWrap<2>(true, N, body, X, Y, BZ);
 }
 
+template<int MAX_THREADS_PER_BLOCK, typename lambda>
+inline void forall_2D_batch(int N, int X, int Y, int BZ, lambda &&body)
+{
+   ForallWrap<2, MAX_THREADS_PER_BLOCK>(true, N, body, X, Y, BZ);
+}
+
 template<typename lambda>
 inline void forall_3D(int N, int X, int Y, int Z, lambda &&body)
 {

general/text.hpp

@@ -50,6 +50,48 @@ 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/petsc.cpp

@@ -4156,20 +4156,31 @@ 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);
 
-   bool b_nonempty = b.Size();
-   if (!B) { B = new PetscParVector(PetscObjectComm(obj), *this, true); }
-   if (!X) { X = new PetscParVector(PetscObjectComm(obj), *this, false, false); }
+   // 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); }
    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.
-   ierr = SNESSolve(snes, B->x, X->x); PCHKERRQ(snes, ierr);
+   // 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);
    X->ResetMemory();
    if (b_nonempty) { B->ResetMemory(); }
 }

makefile

@@ -123,15 +123,20 @@ EXAMPLE_SUBDIRS = amgx caliper ginkgo hiop petsc pumi sundials superlu moonolith
 EXAMPLE_DIRS := examples $(addprefix examples/,$(EXAMPLE_SUBDIRS))
 EXAMPLE_TEST_DIRS := examples
 
-MINIAPP_SUBDIRS = common electromagnetics meshing performance tools \
+MINIAPP_ALL_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
+ fluids/navier fluids/schrodinger-flow plasma plasma/pic
+MINIAPP_RECURSIVE_SUBDIRS = plasma/pic
+MINIAPP_SUBDIRS := $(filter-out \
+   $(MINIAPP_RECURSIVE_SUBDIRS),$(MINIAPP_ALL_SUBDIRS))
+MINIAPP_ALL_DIRS := $(addprefix miniapps/,$(MINIAPP_ALL_SUBDIRS))
 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)
+ toys gslib shifted dpg diag-smoothers fluids/navier plasma plasma/pic)
 
+EM_ALL_DIRS = $(EXAMPLE_DIRS) $(MINIAPP_ALL_DIRS)
 EM_DIRS = $(EXAMPLE_DIRS) $(MINIAPP_DIRS)
 
 TEST_SUBDIRS = unit
@@ -146,7 +151,7 @@ MFEM_BUILD_DIR ?= .
 BUILD_DIR := $(MFEM_BUILD_DIR)
 BUILD_REAL_DIR := $(abspath $(BUILD_DIR))
 ifneq ($(BUILD_REAL_DIR),$(MFEM_REAL_DIR))
-   BUILD_SUBDIRS = $(DIRS) config $(EM_DIRS) doc $(TEST_DIRS)
+   BUILD_SUBDIRS = $(DIRS) config $(EM_ALL_DIRS) doc $(TEST_DIRS)
    CONFIG_FILE_DEF = -DMFEM_CONFIG_FILE='"$(BUILD_REAL_DIR)/config/_config.hpp"'
    BLD := $(if $(BUILD_REAL_DIR:$(CURDIR)=),$(BUILD_DIR)/,)
    $(if $(word 2,$(BLD)),$(error Spaces in BLD = "$(BLD)" are not supported))
@@ -483,10 +488,10 @@ $(OBJECT_FILES): $(BLD)%.o: $(SRC)%.cpp $(CONFIG_MK)
 
 all: examples miniapps $(TEST_DIRS)
 
-.PHONY: miniapps $(EM_DIRS) $(TEST_DIRS)
+.PHONY: miniapps $(EM_ALL_DIRS) $(TEST_DIRS)
 miniapps: $(MINIAPP_DIRS)
 $(MINIAPP_USE_COMMON): miniapps/common
-$(EM_DIRS) $(TEST_DIRS): lib
+$(EM_ALL_DIRS) $(TEST_DIRS): lib
 	$(MAKE) -C $(BLD)$(@)
 
 .PHONY: doc
@@ -694,7 +699,7 @@ local-config:
 .PHONY: build-config
 build-config:
 	for d in $(BUILD_SUBDIRS); do mkdir -p $(BLD)$${d}; done
-	for dir in "" $(addsuffix /,config $(EM_DIRS) doc $(TEST_DIRS)); do \
+	for dir in "" $(addsuffix /,config $(EM_ALL_DIRS) doc $(TEST_DIRS)); do\
 	   printf "# Auto-generated file.\n%s\n%s\n" \
 	      "MFEM_DIR = $(MFEM_REAL_DIR)" \
 	      "include \$$(MFEM_DIR)/$${dir}makefile" \
@@ -796,13 +801,15 @@ status info:
 
 ASTYLE = $(ASTYLE_BIN) --options=$(SRC)config/mfem.astylerc
 ASTYLE_VER = "Artistic Style Version 3.1"
-FORMAT_FILES = $(foreach dir,$(DIRS) $(EM_DIRS) config,$(dir)/*.?pp)
+FORMAT_FILES = $(foreach dir,$(DIRS) $(EM_ALL_DIRS) config,$(dir)/*.?pp)
 TESTS_SUBDIRS = unit benchmarks convergence mem_manager par-mesh-format
-UNIT_TESTS_SUBDIRS = general linalg mesh fem miniapps ceed enzyme
-MINIAPPS_SUBDIRS = dpg/util hooke/operators hooke/preconditioners hooke/materials hooke/kernels
+UNIT_TESTS_SUBDIRS = general linalg mesh fem miniapps ceed enzyme dfem
+MINIAPPS_SUBDIRS = dpg/util hooke/operators hooke/preconditioners \
+   hooke/materials hooke/kernels
 FORMAT_FILES += $(foreach dir,$(TESTS_SUBDIRS),tests/$(dir)/*.?pp)
 FORMAT_FILES += $(foreach dir,$(UNIT_TESTS_SUBDIRS),tests/unit/$(dir)/*.?pp)
 FORMAT_FILES += $(foreach dir,$(MINIAPPS_SUBDIRS),miniapps/$(dir)/*.?pp)
+FORMAT_FILES += config/cmake/config.hpp.in config/config.hpp.in mfem*.hpp
 FORMAT_EXCLUDE = general/tinyxml2.cpp tests/unit/catch.hpp
 FORMAT_LIST = $(filter-out $(FORMAT_EXCLUDE),$(wildcard $(FORMAT_FILES)))
 
@@ -833,14 +840,29 @@ mfem_check_command = \
 # Verify the C++ code styling in MFEM and check that std::cout and std::cerr are
 # not used in the library (use mfem::out and mfem::err instead).
 style:
-	@echo "Applying C++ code style..."
 	@astyle_version="$$($(ASTYLE_BIN) --version)";\
 	 if [ "$$astyle_version" != $(ASTYLE_VER) ]; then\
 	    printf "%s\n" "Invalid astyle version: '$$astyle_version'"\
 	           "Please use: '"$(ASTYLE_VER)"'";\
 	    exit 1;\
 	 fi
-	@err_code=0;\
+	@err_code=0; \
+	if command -v git 2>&1 > /dev/null && [ -d $(MFEM_DIR)/.git ]; then \
+	   echo "Checking if all git files are selected for formatting ..."; \
+	   ls -1 $(FORMAT_FILES) | sort > format-files-make.txt; \
+	   git -C $(MFEM_DIR) ls-files '*.[ch]pp*' | sort \
+	      > format-files-git.txt; \
+	   cat format-files-make.txt format-files-git.txt | sort | uniq \
+	      > format-files-make-plus-git.txt; \
+	   rm -f format-files-git.txt; \
+	   $(call mfem_check_command,\
+	      diff format-files-make.txt format-files-make-plus-git.txt | \
+	      grep "^> ",\
+	      "All git files are selected for formatting",\
+	      "The above git files are NOT selected for formatting"); \
+	   rm -f format-files-make.txt format-files-make-plus-git.txt; \
+	fi; \
+	echo "Applying C++ code style...";\
 	$(call mfem_check_command,\
 	    $(ASTYLE) $(FORMAT_LIST) | grep Formatted,\
 	    "No source files were changed",\

mesh/mesh.cpp

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

mesh/mesh_operators.cpp

@@ -63,7 +63,6 @@ ThresholdRefiner::ThresholdRefiner(ErrorEstimator &est)
 
    threshold = 0.0;
    num_marked_elements = 0LL;
-   current_sequence = -1;
 
    non_conforming = -1;
    nc_limit = 0;
@@ -87,7 +86,6 @@ 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; }
@@ -149,7 +147,6 @@ 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,7 +188,6 @@ protected:
    long long num_marked_elements;
 
    Array<Refinement> marked_elements;
-   long current_sequence;
 
    int non_conforming;
    int nc_limit;

mesh/pmesh.cpp

@@ -5639,6 +5639,12 @@ 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,15 +227,29 @@ public:
                                            const ParGridFunction &dst);
 
    /**
-   * @brief Check if ParMesh @a m is a ParSubMesh.
+   * @brief Check if Mesh @a m is a ParSubMesh.
    *
-   * @param m The input ParMesh
+   * @param m The input Mesh
    */
-   static bool IsParSubMesh(const ParMesh *m)
+   static bool IsParSubMesh(const Mesh *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,6 +225,20 @@ 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,19 +43,39 @@ endif()
   
   # Add the corresponding tests to the "test" target
   if (MFEM_ENABLE_TESTING)
-    add_test(NAME tesla_np=4
+    add_test(NAME tesla_1_np=${MFEM_MPI_NP}
       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 volta_np=4
+    add_test(NAME tesla_2_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"
+      $<TARGET_FILE:tesla>
+      -no-vis -maxit 2 -m ../../data/inline-hex.mesh -ubbc "0 0 1"
       ${MPIEXEC_POSTFLAGS})
 
-    add_test(NAME joule_np=4
+    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}
       COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
       ${MPIEXEC_PREFLAGS}
       $<TARGET_FILE:joule>
@@ -63,12 +83,41 @@ endif()
       ${MPIEXEC_POSTFLAGS})
 
     if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
-      add_test(NAME maxwell_np=4
+      add_test(NAME maxwell_np=${MFEM_MPI_NP}
         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 -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)
+	-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)
 
 # Testing: "test" target and mfem-test* variables are defined in config/test.mk
 

miniapps/gslib/CMakeLists.txt

@@ -46,8 +46,16 @@ if (MFEM_USE_GSLIB)
 
   if (MFEM_ENABLE_TESTING)
     foreach (test "schwarz_ex1" "field-diff" "findpts" "field-interp")
-      add_test(NAME ${test}
-        COMMAND $<TARGET_FILE:${test}> -no-vis)
+      if (MFEM_USE_MPI)
+        add_test(NAME ${test}
+          COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} 1
+          ${MPIEXEC_PREFLAGS}
+          $<TARGET_FILE:${test}> -no-vis
+          ${MPIEXEC_POSTFLAGS})
+      else()
+        add_test(NAME ${test}
+          COMMAND $<TARGET_FILE:${test}> -no-vis)
+      endif()
     endforeach()
   endif()
 

miniapps/gslib/makefile

@@ -83,11 +83,16 @@ include $(MFEM_TEST_MK)
 
 # Testing: Parallel vs. serial runs
 RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
+ifeq ($(MFEM_USE_MPI),YES)
+   RUN_MPI_SER = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) 1
+else
+   RUN_MPI_SER =
+endif
 TEST_NAME := GSLIB miniapp
 %-test-par: %
 	@$(call mfem-test,$<, $(RUN_MPI), $(TEST_NAME))
 %-test-seq: %
-	@$(call mfem-test,$<,, $(TEST_NAME))
+	@$(call mfem-test,$<, $(RUN_MPI_SER), $(TEST_NAME))
 
 # 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<const int>(sqrt(n2));
+   const int n = static_cast<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<const int>(sqrt(B.Size()));
+   const int n = static_cast<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/plasma/CMakeLists.txt

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

miniapps/plasma/makefile

@@ -14,9 +14,6 @@ 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)
@@ -29,6 +26,14 @@ else
    MINIAPPS = $(PAR_MINIAPPS) $(SEQ_MINIAPPS)
 endif
 
+PLASMA_SUBDIRS = pic
+
+SUBDIRS_ALL = $(addsuffix /all,$(PLASMA_SUBDIRS))
+SUBDIRS_TEST = $(addsuffix /test,$(PLASMA_SUBDIRS))
+SUBDIRS_TEST_NOCLEAN = $(addsuffix /test-noclean,$(PLASMA_SUBDIRS))
+SUBDIRS_CLEAN = $(addsuffix /clean,$(PLASMA_SUBDIRS))
+SUBDIRS_TPRINT = $(addsuffix /test-print,$(PLASMA_SUBDIRS))
+
 .SUFFIXES:
 .SUFFIXES: .o .cpp .mk
 .PHONY: all lib-common clean clean-build clean-exec
@@ -38,26 +43,24 @@ 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=),,\
-   $(if $(MFEM_USE_CUDA:YES=),$(CXX_XLINKER),$(CUDA_XLINKER))-rpath,$(abspath\
-   $(MFEM_BUILD_DIR)/miniapps/common))
-
-COMMON_O=
+   $(MFEM_XLINKER)-rpath,$(abspath $(MFEM_BUILD_DIR)/miniapps/common))
 
 # Remove built-in rules
 %: %.cpp
 %.o: %.cpp
 
-all: $(MINIAPPS)
+all: $(MINIAPPS) $(SUBDIRS_ALL) 
+
+.PHONY: $(SUBDIRS_ALL) $(SUBDIRS_TEST) $(SUBDIRS_TEST_NOCLEAN) \
+   $(SUBDIRS_CLEAN) $(SUBDIRS_TPRINT)
+$(SUBDIRS_ALL) $(SUBDIRS_TEST) $(SUBDIRS_TEST_NOCLEAN) $(SUBDIRS_CLEAN):
+	$(MAKE) -C $(@D) $(@F)
+$(SUBDIRS_TPRINT):
+	@$(MAKE) -C $(@D) $(@F)
 
 # Rules for building the miniapps
-%: $(SRC)%.cpp $(COMMON_O) $(MFEM_LIB_FILE) $(CONFIG_MK) | lib-common
-	$(MFEM_CXX) $(MFEM_LINK_FLAGS) $< -o $@ $(COMMON_O) $(COMMON_LIB) \
-        $(MFEM_LIBS)
-
-# Rules for compiling miniapp dependencies
-$(COMMON_O) $(addsuffix _solver.o,$(MINIAPPS)): \
-%.o: $(SRC)%.cpp $(SRC)%.hpp $(CONFIG_MK)
-	$(MFEM_CXX) $(MFEM_FLAGS) -c $(<) -o $(@)
+%: $(SRC)%.cpp $(MFEM_LIB_FILE) $(CONFIG_MK) | lib-common
+	$(MFEM_CXX) $(MFEM_LINK_FLAGS) $< -o $@ $(COMMON_LIB) $(MFEM_LIBS)
 
 # Rule for building lib-common
 lib-common:
@@ -65,6 +68,9 @@ lib-common:
 
 MFEM_TESTS = MINIAPPS
 include $(MFEM_TEST_MK)
+test: $(SUBDIRS_TEST)
+test-noclean: $(SUBDIRS_TEST_NOCLEAN)
+test-print: $(SUBDIRS_TPRINT)
 
 # Testing: Specific execution options
 RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
@@ -75,11 +81,10 @@ RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
 $(MFEM_LIB_FILE):
 	$(error The MFEM library is not built)
 
-clean: clean-build clean-exec
+clean: clean-build clean-exec $(SUBDIRS_CLEAN)
 
 clean-build:
 	rm -f *.o *~ $(SEQ_MINIAPPS) $(PAR_MINIAPPS)
 	rm -rf *.dSYM *.TVD.*breakpoints
 
 clean-exec:
-

miniapps/plasma/pic/CMakeLists.txt

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

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

@@ -0,0 +1,86 @@
+# 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)
+
+SEQ_MINIAPPS =
+PAR_MINIAPPS =
+
+ifeq ($(MFEM_USE_GSLIB),YES)
+   PAR_MINIAPPS += electrostatic-pic
+endif
+
+ifeq ($(MFEM_USE_MPI),NO)
+   MINIAPPS = $(SEQ_MINIAPPS)
+else
+   MINIAPPS = $(PAR_MINIAPPS) $(SEQ_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
+%: $(SRC)%.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

miniapps/tools/CMakeLists.txt

@@ -61,6 +61,10 @@ if (MFEM_USE_MPI)
       LIBRARIES mfem-common)
       add_dependencies(gridfunction-bounds copy_miniapps_tools_data)
 
+    add_mfem_miniapp(random-gridfunction-bounds
+      MAIN random-gridfunction-bounds.cpp
+      LIBRARIES mfem)
+
     add_mfem_miniapp(plor-transfer
       MAIN plor-transfer.cpp LIBRARIES mfem)
 

miniapps/tools/makefile

@@ -23,7 +23,8 @@ MFEM_LIB_FILE = mfem_is_not_built
 SEQ_MINIAPPS = display-basis load-dc convert-dc get-values lor-transfer \
 	tmop-check-metric tmop-metric-magnitude compare-dc
 
-PAR_MINIAPPS = nodal-transfer plor-transfer gridfunction-bounds
+PAR_MINIAPPS = nodal-transfer plor-transfer gridfunction-bounds \
+	random-gridfunction-bounds
 
 ifeq ($(MFEM_USE_MPI),NO)
    MINIAPPS = $(SEQ_MINIAPPS)
@@ -79,7 +80,7 @@ RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
 # Do not test: display-basis, load-dc, convert-dc, get-values, lor-transfer, plor-transfer
 NO_TEST_APPS = display-basis load-dc convert-dc get-values lor-transfer \
 	plor-transfer tmop-check-metric tmop-metric-magnitude gridfunction-bounds \
-	compare-dc
+	random-gridfunction-bounds compare-dc
 $(foreach app,$(NO_TEST_APPS),$(app)-test-seq $(app)-test-par):
 	@true
 

miniapps/tools/random-gridfunction-bounds.cpp

@@ -0,0 +1,263 @@
+// 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.
+//
+//    ---------------------------------------------------------------------
+//     Compute bounds of a random grid function on a generated tensor mesh
+//    ---------------------------------------------------------------------
+//
+// This miniapp generates a 1D segment mesh or 2D quad mesh, builds a random
+// discontinuous grid function, computes element-wise piecewise linear bounds,
+// and visualizes the input field together with the lower and upper bounds.
+//
+// Compile with: make random-gridfunction-bounds
+//
+// Sample runs:
+//   mpirun -np 4 random-gridfunction-bounds
+//   mpirun -np 4 random-gridfunction-bounds -nx 64 -o 6 -ref 3 -d hip
+
+#include "mfem.hpp"
+
+#include <algorithm>
+#include <type_traits>
+
+using namespace mfem;
+using namespace std;
+
+void VisualizeField(ParMesh &pmesh, ParGridFunction &input,
+                    char *title, int pos_x, int pos_y);
+
+int main(int argc, char *argv[])
+{
+   Mpi::Init(argc, argv);
+   Hypre::Init();
+
+   int dim = 2;
+   int nx = 16;
+   int order = 4;
+   int num_comp = 2;
+   int ref = 2;
+   int niter = 1000;
+   int seed = 12345;
+   bool kernel_only = true;
+   bool visualization = false;
+   const char *device_config = "cpu";
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&dim, "-dim", "--dimension",
+                  "Dimension of the generated tensor-product mesh (1 or 2).");
+   args.AddOption(&nx, "-nx", "--num-elements",
+                  "Number of elements in each mesh direction.");
+   args.AddOption(&order, "-o", "--order",
+                  "Polynomial degree of the random discontinuous field.");
+   args.AddOption(&num_comp, "-nc", "--num-components",
+                  "Number of vector components in the ParFiniteElementSpace.");
+   args.AddOption(&ref, "-ref", "--piecewise-linear-ref-factor",
+                  "Scaling factor for the resolution of the piecewise linear "
+                  "bounds. If less than 2, the resolution is picked "
+                  "automatically.");
+   args.AddOption(&niter, "-ni", "--num-iters",
+                  "Number of times to evaluate the bounds.");
+   args.AddOption(&seed, "-rs", "--random-seed",
+                  "Random seed used to initialize the field.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&kernel_only, "-ko", "--kernel-only",
+                  "-no-ko", "--no-kernel-only",
+                  "Run only PLBound::GetElementBoundsKernel on a prebuilt "
+                  "element E-vector.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.ParseCheck();
+
+   MFEM_VERIFY(dim == 1 || dim == 2, "dim must be 1 or 2.");
+   MFEM_VERIFY(nx > 0, "nx must be positive.");
+   MFEM_VERIFY(order >= 0, "order must be non-negative.");
+   MFEM_VERIFY(num_comp > 0, "num_comp must be positive.");
+   MFEM_VERIFY(niter > 0, "niter must be positive.");
+
+   Device device(device_config);
+   if (Mpi::Root()) { device.Print(); }
+
+   Mesh mesh = (dim == 1) ?
+               Mesh::MakeCartesian1D(nx, 1.0) :
+               Mesh::MakeCartesian2D(nx, nx, Element::QUADRILATERAL, true,
+                                     1.0, 1.0);
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+
+   const int mesh_dim = pmesh.Dimension();
+   L2_FECollection fec(order, mesh_dim, BasisType::GaussLobatto);
+   ParFiniteElementSpace fes(&pmesh, &fec, num_comp, Ordering::byNODES);
+   ParGridFunction input(&fes);
+   input.Randomize(seed + Mpi::WorldRank());
+   input.UseDevice(true);
+
+   L2_FECollection fec_pc(0, mesh_dim);
+   ParFiniteElementSpace fes_pc(&pmesh, &fec_pc, num_comp, Ordering::byNODES);
+   ParGridFunction lowerb(&fes_pc), upperb(&fes_pc);
+   Vector lower_vec, upper_vec;
+
+   PLBound plb(&fes, ref*(fes.GetMaxElementOrder() + 1));
+   if (kernel_only)
+   {
+      const FiniteElement &fe = *fes.GetTypicalFE();
+      const int rdim = fe.GetDim();
+      const int nd = fe.GetDof();
+      const int fes_dim = fes.GetVDim();
+      Vector e_vec(nd*fes_dim*fes.GetNE(), Device::GetDeviceMemoryType());
+      e_vec.UseDevice(true);
+
+      const ElementRestrictionOperator *elem_restr =
+         fes.GetElementRestriction(ElementDofOrdering::LEXICOGRAPHIC);
+      MFEM_VERIFY(elem_restr != nullptr,
+                  "Element restriction is required for kernel-only mode.");
+      elem_restr->Mult(input, e_vec);
+
+      for (int i = 0; i < niter; i++)
+      {
+         plb.GetElementBoundsKernel(rdim, fes_dim, e_vec, lower_vec, upper_vec);
+      }
+   }
+   else
+   {
+      for (int i = 0; i < niter; i++)
+      {
+         input.GetElementBounds(plb, lower_vec, upper_vec);
+      }
+   }
+
+   const real_t *lower_data = lower_vec.HostRead();
+   const real_t *upper_data = upper_vec.HostRead();
+
+   // Build a host reference from the lexicographic E-vector and the scalar
+   // PLBound::GetNDBounds path to avoid re-entering the device dispatch.
+   const bool use_dev = input.UseDevice();
+   PLBound plb_host(&fes, ref*(fes.GetMaxElementOrder() + 1));
+   Vector lower_ref, upper_ref;
+   const FiniteElement &fe = *fes.GetTypicalFE();
+   const int rdim = fe.GetDim();
+   const int nd = fe.GetDof();
+   const int nel = fes.GetNE();
+   const int fes_dim = fes.GetVDim();
+   Vector e_vec_ref(nd*fes_dim*nel);
+   lower_ref.SetSize(nel*fes_dim);
+   upper_ref.SetSize(nel*fes_dim);
+   const ElementRestrictionOperator *elem_restr =
+      fes.GetElementRestriction(ElementDofOrdering::LEXICOGRAPHIC);
+   MFEM_VERIFY(elem_restr != nullptr,
+               "Element restriction is required for host reference.");
+   input.UseDevice(false);
+   input.HostRead();
+   elem_restr->Mult(input, e_vec_ref);
+   input.UseDevice(use_dev);
+   const real_t *e_ref_data = e_vec_ref.HostRead();
+
+   for (int d = 0; d < fes_dim; d++)
+   {
+      for (int e = 0; e < nel; e++)
+      {
+         Vector coeff(nd);
+         for (int i = 0; i < nd; i++)
+         {
+            coeff(i) = e_ref_data[i + nd*(d + fes_dim*e)];
+         }
+         Vector lower_c, upper_c;
+         plb_host.GetNDBounds(rdim, coeff, lower_c, upper_c);
+         lower_ref(e + d*nel) = lower_c.Min();
+         upper_ref(e + d*nel) = upper_c.Max();
+      }
+   }
+   const real_t *lower_ref_data = lower_ref.HostRead();
+   const real_t *upper_ref_data = upper_ref.HostRead();
+
+   MFEM_VERIFY(lower_vec.Size() == lower_ref.Size() &&
+               upper_vec.Size() == upper_ref.Size(),
+               "Reference element-bound vectors have inconsistent sizes.");
+
+   real_t lower_diff = 0.0;
+   real_t upper_diff = 0.0;
+   for (int i = 0; i < lower_vec.Size(); i++)
+   {
+      lower_diff = std::max(lower_diff,
+                            std::abs(lower_data[i] - lower_ref_data[i]));
+   }
+   for (int i = 0; i < upper_vec.Size(); i++)
+   {
+      upper_diff = std::max(upper_diff,
+                            std::abs(upper_data[i] - upper_ref_data[i]));
+   }
+   MPI_Allreduce(MPI_IN_PLACE, &lower_diff, 1, MPITypeMap<real_t>::mpi_type,
+                 MPI_MAX, pmesh.GetComm());
+   MPI_Allreduce(MPI_IN_PLACE, &upper_diff, 1, MPITypeMap<real_t>::mpi_type,
+                 MPI_MAX, pmesh.GetComm());
+
+   const real_t verify_tol = std::is_same<real_t, float>::value ?
+                             real_t(1.0e-5) : real_t(1.0e-12);
+   MFEM_VERIFY(lower_diff <= verify_tol && upper_diff <= verify_tol,
+               "Device element bounds do not match host reference.");
+
+   lowerb = lower_vec;
+   upperb = upper_vec;
+
+   real_t lower_min = lowerb.Min();
+   real_t upper_max = upperb.Max();
+   MPI_Allreduce(MPI_IN_PLACE, &lower_min, 1, MPITypeMap<real_t>::mpi_type,
+                 MPI_MIN, pmesh.GetComm());
+   MPI_Allreduce(MPI_IN_PLACE, &upper_max, 1, MPITypeMap<real_t>::mpi_type,
+                 MPI_MAX, pmesh.GetComm());
+
+   if (Mpi::Root())
+   {
+      cout << "dim: " << mesh_dim << '\n'
+           << "nx: " << nx << '\n'
+           << "order: " << order << '\n'
+           << "num components: " << num_comp << '\n'
+           << "PL bound control-point factor: " << ref << '\n'
+           << "iterations: " << niter << '\n'
+           << "kernel-only mode: " << (kernel_only ? "yes" : "no") << '\n'
+           << "host/device lower max diff: " << lower_diff << '\n'
+           << "host/device upper max diff: " << upper_diff << '\n'
+           << "global lower bound minimum: " << lower_min << '\n'
+           << "global upper bound maximum: " << upper_max << endl;
+   }
+
+   if (visualization)
+   {
+      char title1[] = "Random input gridfunction";
+      char title2[] = "Element-wise lower bound";
+      char title3[] = "Element-wise upper bound";
+      VisualizeField(pmesh, input, title1, 0, 0);
+      VisualizeField(pmesh, lowerb, title2, 450, 0);
+      VisualizeField(pmesh, upperb, title3, 900, 0);
+   }
+
+   return 0;
+}
+
+void VisualizeField(ParMesh &pmesh, ParGridFunction &input,
+                    char *title, int pos_x, int pos_y)
+{
+   socketstream sock;
+   if (pmesh.GetMyRank() == 0)
+   {
+      sock.open("localhost", 19916);
+      sock << "solution\n";
+   }
+   pmesh.PrintAsOne(sock);
+   input.SaveAsOne(sock);
+   if (pmesh.GetMyRank() == 0)
+   {
+      sock << "window_title '" << title << "'\n"
+           << "window_geometry "
+           << pos_x << " " << pos_y << " " << 400 << " " << 400 << "\n"
+           << "keys jRmclApppppppppppp//]]]]]]]]" << endl;
+   }
+}

tests/unit/CMakeLists.txt

@@ -71,10 +71,12 @@ 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

@@ -0,0 +1,118 @@
+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 (rank != procs[i])
+         if (static_cast<unsigned>(rank) != procs[i])
          {
             wrong_proc_count++;
          }

tests/unit/fem/test_var_order.cpp

@@ -271,6 +271,8 @@ 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();
@@ -698,7 +700,14 @@ static void TestSolveVec(FiniteElementSpace &fespace)
 
    GridFunction x(&fespace);
    x = 0.0;
-   x.ProjectBdrCoefficient(exsol, ess_attr);
+   if (x.FESpace()->GetTypicalBE()->GetRangeDim() == 0)
+   {
+      x.ProjectBdrCoefficientNormal(exsol, ess_attr);
+   }
+   else
+   {
+      x.ProjectBdrCoefficientTangent(exsol, ess_attr);
+   }
 
    // Assemble the linear form
    LinearForm lf(&fespace);
@@ -1082,7 +1091,14 @@ static void TestSolveParVec(ParFiniteElementSpace &fespace)
 
    ParGridFunction x(&fespace);
    x = 0.0;
-   x.ProjectBdrCoefficient(exsol, ess_attr);
+   if (x.FESpace()->GetTypicalBE()->GetRangeDim() == 0)
+   {
+      x.ProjectBdrCoefficientNormal(exsol, ess_attr);
+   }
+   else
+   {
+      x.ProjectBdrCoefficientTangent(exsol, ess_attr);
+   }
 
    // Assemble the linear form
    ParLinearForm lf(&fespace);

tests/unit/general/test_text.cpp

@@ -30,3 +30,29 @@ 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

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

@@ -0,0 +1,108 @@
+// 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,8 +486,14 @@ void multidomain_test_3d(FECType fec_type)
          {
             cylinder_gf.ProjectCoefficient(vcoeff);
             outer_gf.ProjectCoefficient(vcoeff);
-            outer_gf.ProjectBdrCoefficient(vzerocoeff,
-                                           outer_cyl_surf_marker);
+            if (fec_type == FECType::RT)
+            {
+               outer_gf.ProjectBdrCoefficientNormal(vzerocoeff, outer_cyl_surf_marker);
+            }
+            else
+            {
+               outer_gf.ProjectBdrCoefficientTangent(vzerocoeff, outer_cyl_surf_marker);
+            }
             outer_gf_ex.ProjectCoefficient(vcoeff);
          }
          ParSubMesh::Transfer(cylinder_gf, outer_gf);
@@ -507,8 +513,14 @@ void multidomain_test_3d(FECType fec_type)
          {
             outer_gf.ProjectCoefficient(vcoeff);
             cylinder_gf.ProjectCoefficient(vcoeff);
-            cylinder_gf.ProjectBdrCoefficient(vzerocoeff,
-                                              cylinder_cyl_surf_marker);
+            if (fec_type == FECType::RT)
+            {
+               cylinder_gf.ProjectBdrCoefficientNormal(vzerocoeff, cylinder_cyl_surf_marker);
+            }
+            else
+            {
+               cylinder_gf.ProjectBdrCoefficientTangent(vzerocoeff, cylinder_cyl_surf_marker);
+            }
             cylinder_gf_ex.ProjectCoefficient(vcoeff);
          }
          ParSubMesh::Transfer(outer_gf, cylinder_gf);