style

Updates CHANGELOG
More fixes for building contact mini-app with make and cmake
2026-04-14 08:54:23 -07:00 · 2026-04-13 21:52:37 -07:00 · 2026-04-13 21:47:15 -07:00 · 2026-04-13 18:54:17 -07:00 · 2026-04-13 12:40:24 -07:00 · 2026-04-13 12:37:18 -07:00
91 changed files with 644 additions and 4706 deletions
@@ -443,10 +443,6 @@ miniapps/diag-smoothers/mg-abs-l1-jacobi
 miniapps/contact/contact
 miniapps/contact/ParaView

-miniapps/plasma/pic/electrostatic-*
-!miniapps/plasma/pic/electrostatic-*.cpp
-miniapps/plasma/pic/*.csv
-
 # Unit test binary and outputs
 tests/unit/output_meshes
 tests/unit/unit_tests
@@ -85,8 +85,3 @@ 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"
@@ -63,8 +63,3 @@ 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"
@@ -32,9 +32,9 @@ mkdir _${BASELINE_TEST} && cd _${BASELINE_TEST}

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


-Version 4.9.1 (development)
-===========================
-
- Policy for AI-assisted contribution added to CONTRIBUTING.md
+Version 4.10 (development)
+==========================

 Discretization improvements
 ---------------------------
@@ -25,6 +23,11 @@ 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.
@@ -42,6 +45,11 @@ New and updated examples and miniapps
 - Electromagnetics/lorentz miniapp has been updated to leverage the ParticleSet
  capability.

+Miscellaneous
+-------------
+- Removes the SidreDataCollection class from MFEM in favor of the 
+  MFEMSidreDataCollection class in the Axom library (https://github.com/llnl/axom).
+

 Version 4.9, released on Dec 11, 2025
 =====================================
@@ -75,12 +75,10 @@ set(XSDK_ENABLE_Fortran OFF)

 # Check if we need to enable C or Fortran.
 if (MFEM_USE_CONDUIT OR
-    MFEM_USE_SIDRE OR
    MFEM_USE_PETSC)
   # This seems to be needed by:
   #  * find_package(BLAS REQUIRED) and
-   #  * find_package(HDF5 REQUIRED) needed, in turn, by:
-   #    - find_package(AXOM REQUIRED)
+   #  * find_package(HDF5 REQUIRED) and
   #  * find_package(PETSc REQUIRED)
   set(XSDK_ENABLE_C ON)
 endif()
@@ -478,11 +476,6 @@ if (MFEM_USE_FMS)
    find_package(FMS REQUIRED fms)
 endif()

-# Axom/Sidre
-if (MFEM_USE_SIDRE)
-  find_package(Axom REQUIRED Axom)
-endif()
-
 # PUMI
 if (MFEM_USE_PUMI)
  # If PUMI_DIR was specified, only link to that directory,
@@ -629,7 +622,7 @@ find_package(Threads REQUIRED)
 #    integers, the METIS header (with 32-bit indices, as used by mfem) needs to
 #    be before SuiteSparse.
 set(MFEM_TPLS OPENMP HYPRE LAPACK BLAS SuperLUDist STRUMPACK METIS SuiteSparse
-    SUNDIALS PETSC SLEPC MUMPS AXOM FMS CONDUIT Ginkgo GNUTLS GSLIB HDF5
+    SUNDIALS PETSC SLEPC MUMPS FMS CONDUIT Ginkgo GNUTLS GSLIB HDF5
    NETCDF MPFR PUMI HIOP POSIXCLOCKS MFEMBacktrace ZLIB OCCA CEED RAJA UMPIRE
    ADIOS2 MKL_CPARDISO MKL_PARDISO AMGX MAGMA CUSPARSE CUBLAS CALIPER CODIPACK
    BENCHMARK PARELAG TRIBOL MPI_CXX HIP HIPBLAS HIPSPARSE MOONOLITH BLITZ
@@ -24,14 +24,6 @@ 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
@@ -452,13 +452,6 @@ MFEM_USE_MPFR = YES/NO
   quadrature rules. When enabled, this option uses the MPFR_* library options,
   see below.

-MFEM_USE_SIDRE = YES/NO
-   Sidre is a component of LLNL's axom project, https://github.com/LLNL/axom,
-   that provides an HDF5-based file format for visualization or restart
-   capability following the Conduit (https://github.com/LLNL/conduit) mesh
-   blueprint specification. When enabled, this option requires installation of
-   HDF5 (see also MFEM_USE_NETCDF), Conduit and LLNL's axom project.
-
 MFEM_USE_SIMD = YES/NO
   Enables the high performance templated classes to use architecture dependent
   SIMD intrinsics instead of the generic implementation of class AutoSIMD in
@@ -778,14 +771,6 @@ The specific libraries and their options are:
  Options: SLEPC_OPT, SLEPC_LIB.
  Versions: SLEPc >= 3.8.0.

- Sidre (optional), part of LLNL's axom project, used when MFEM_USE_SIDRE = YES.
-  Starting with MFEM v4.1, Axom version 0.3.1 or later is required.
-  URL: https://github.com/LLNL/axom
-       https://github.com/LLNL/conduit (Conduit)
-       https://support.hdfgroup.org/HDF5 (HDF5)
-  Options: SIDRE_OPT, SIDRE_LIB.
-  Versions: Axom >= 0.3.1.
-
 - Conduit (optional), used when MFEM_USE_CONDUIT = YES. Conduit Mesh Blueprint
  support requires Conduit >= v0.3.1 and VisIt >= v2.13.1 to read the output.
  URL: https://github.com/LLNL/conduit (Conduit)
@@ -1069,7 +1054,6 @@ MFEM_USE_OCCA
 MFEM_USE_CEED
 MFEM_USE_RAJA
 MFEM_USE_UMPIRE
-MFEM_USE_SIDRE
 MFEM_USE_MOONOLITH
 MFEM_USE_CALIPER
 MFEM_USE_FMS
@@ -1133,7 +1117,6 @@ The CMake build system adds auto-detection for the following packages/libraries:
 - OCCA
 - RAJA
 - UMPIRE
- - AXOM - Used when MFEM_USE_SIDRE is enabled
 - MOONOLITH
 - CALIPER
 - FMS
@@ -248,10 +248,6 @@ IF (DEFINED TPL_ENABLE_MPFR)
  SET(MFEM_USE_MPFR ${TPL_ENABLE_MPFR} CACHE BOOL "Enable MPFR usage." FORCE)
 ENDIF()

-IF (DEFINED TPL_ENABLE_SIDRE)
-  SET(MFEM_USE_SIDRE ${TPL_ENABLE_SIDRE} CACHE BOOL "Enable Axom/Sidre usage" FORCE)
-ENDIF()
-
 IF (DEFINED TPL_ENABLE_FMS)
  SET(MFEM_USE_FMS ${TPL_ENABLE_FMS} CACHE BOOL "Enable FMS usage" FORCE)
 ENDIF()
@@ -46,7 +46,6 @@ set(MFEM_USE_NETCDF @MFEM_USE_NETCDF@)
 set(MFEM_USE_PETSC @MFEM_USE_PETSC@)
 set(MFEM_USE_SLEPC @MFEM_USE_SLEPC@)
 set(MFEM_USE_MPFR @MFEM_USE_MPFR@)
-set(MFEM_USE_SIDRE @MFEM_USE_SIDRE@)
 set(MFEM_USE_FMS @MFEM_USE_FMS@)
 set(MFEM_USE_CONDUIT @MFEM_USE_CONDUIT@)
 set(MFEM_USE_PUMI @MFEM_USE_PUMI@)
@@ -109,10 +108,6 @@ if (MFEM_USE_RAJA)
   find_dependency(RAJA)
 endif()

-if (MFEM_USE_UMPIRE)
-   find_dependency(umpire)
-endif()
-
 if (NOT TARGET mfem)
  include(${CMAKE_CURRENT_LIST_DIR}/MFEMTargets.cmake)
 endif (NOT TARGET mfem)
@@ -120,9 +120,6 @@
 // Enable secure socket streams based on the GNUTLS library.
 #cmakedefine MFEM_USE_GNUTLS

-// Enable Sidre support.
-#cmakedefine MFEM_USE_SIDRE
-
 // Enable the use of SIMD in the high performance templated classes.
 #cmakedefine MFEM_USE_SIMD

@@ -0,0 +1,24 @@
+# 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.
+
+# Defines the following variables:
+#   - ADIAK_FOUND
+#   - ADIAK_LIBRARIES
+#   - ADIAK_INCLUDE_DIRS
+
+include(MfemCmakeUtilities)
+
+mfem_find_package(Adiak ADIAK ADIAK_DIR
+      "include" "adiak.h"
+      "lib" "adiak"
+      "Paths to headers required by Adiak."
+      "Libraries required by Adiak.")
+
@@ -13,6 +13,9 @@
 #   - AXOM_FOUND
 #   - AXOM_LIBRARIES
 #   - AXOM_INCLUDE_DIRS
+#
+# MFEM itself does not depend on Axom, however Tribol does. This module exists
+# to support MFEM's Tribol integration (e.g. the contact miniapp).

 include(MfemCmakeUtilities)
 # Note: components are enabled based on the find_package() parameters.
@@ -0,0 +1,36 @@
+# 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.
+
+# Defines the following variables:
+#   - CAMP_FOUND
+#   - CAMP_LIBRARIES
+#   - CAMP_INCLUDE_DIRS
+
+include(MfemCmakeUtilities)
+
+mfem_find_package(CAMP CAMP CAMP_DIR
+      "include" "camp/camp.hpp"
+      "lib" "camp"
+      "Paths to headers required by CAMP."
+      "Libraries required by CAMP.")
+
+# RAJA commonly lists "camp" in INTERFACE_LINK_LIBRARIES. If there is no CMake
+# target named "camp", CMake treats it as a bare library name (-lcamp).
+if (CAMP_FOUND AND NOT TARGET camp)
+   list(GET CAMP_LIBRARIES 0 _camp_lib0)
+   add_library(camp UNKNOWN IMPORTED)
+   set_target_properties(camp PROPERTIES
+      IMPORTED_LOCATION "${_camp_lib0}"
+      INTERFACE_INCLUDE_DIRECTORIES "${CAMP_INCLUDE_DIRS}")
+   set(CAMP_LIBRARIES "camp" CACHE STRING "CAMP imported target." FORCE)
+   unset(_camp_lib0)
+endif()
+
@@ -21,6 +21,21 @@ mfem_find_package(Caliper CALIPER CALIPER_DIR
      "Paths to headers required by Caliper."
      "Libraries required by Caliper.")

+# Some downstream CMake packages (notably RAJA) may list "caliper" in their
+# INTERFACE_LINK_LIBRARIES. If there is no CMake target named "caliper", CMake
+# treats it as a bare library name and will pass -lcaliper to the linker.
+# Create a minimal imported target when we only located the library by path.
+if (CALIPER_FOUND AND NOT TARGET caliper)
+   list(GET CALIPER_LIBRARIES 0 _caliper_lib0)
+   add_library(caliper UNKNOWN IMPORTED)
+   set_target_properties(caliper PROPERTIES
+      IMPORTED_LOCATION "${_caliper_lib0}"
+      INTERFACE_INCLUDE_DIRECTORIES "${CALIPER_INCLUDE_DIRS}")
+   # Prefer linking via the target.
+   set(CALIPER_LIBRARIES "caliper" CACHE STRING "Caliper imported target." FORCE)
+   unset(_caliper_lib0)
+endif()
+
 # Append adiak path/lib if the user provided ADIAK_DIR
 if(ADIAK_DIR AND EXISTS ${ADIAK_DIR})
    find_package(adiak NO_DEFAULT_PATH REQUIRED PATHS ${ADIAK_DIR}/lib/cmake/adiak ${ADIAK_DIR})
@@ -14,12 +14,12 @@
 #   - UMPIRE_LIBRARIES
 #   - UMPIRE_INCLUDE_DIRS

-if (NOT umpire_ROOT AND UMPIRE_DIR)
-  set(umpire_ROOT ${UMPIRE_DIR})
+if (NOT umpire_DIR AND UMPIRE_DIR)
+  set(umpire_DIR ${UMPIRE_DIR}/lib/cmake/umpire)
 endif()
 message(STATUS "Looking for UMPIRE ...")
 message(STATUS "   in UMPIRE_DIR = ${UMPIRE_DIR}")
-message(STATUS "      umpire_ROOT = ${umpire_ROOT}")
+message(STATUS "      umpire_DIR = ${umpire_DIR}")
 find_package(umpire CONFIG)
 set(UMPIRE_FOUND ${umpire_FOUND})
 set(UMPIRE_LIBRARIES "umpire")
@@ -878,7 +878,7 @@ function(mfem_export_mk_files)
      MFEM_USE_SUITESPARSE MFEM_USE_SUPERLU MFEM_USE_SUPERLU5 MFEM_USE_MUMPS
      MFEM_USE_STRUMPACK MFEM_USE_GINKGO MFEM_USE_AMGX MFEM_USE_MAGMA
      MFEM_USE_GNUTLS MFEM_USE_NETCDF MFEM_USE_PETSC MFEM_USE_SLEPC
-      MFEM_USE_MPFR MFEM_USE_SIDRE MFEM_USE_FMS MFEM_USE_CONDUIT MFEM_USE_PUMI
+      MFEM_USE_MPFR MFEM_USE_FMS MFEM_USE_CONDUIT MFEM_USE_PUMI
      MFEM_USE_HIOP MFEM_USE_GSLIB MFEM_USE_CUDA MFEM_USE_HIP MFEM_USE_RAJA
      MFEM_USE_OCCA MFEM_USE_CEED MFEM_USE_CALIPER MFEM_USE_UMPIRE MFEM_USE_SIMD
      MFEM_USE_ADIOS2 MFEM_USE_MKL_CPARDISO MFEM_USE_MKL_PARDISO
@@ -120,9 +120,6 @@
 // Enable secure socket streams based on the GNUTLS library.
 // #define MFEM_USE_GNUTLS

-// Enable Sidre support.
-// #define MFEM_USE_SIDRE
-
 // Enable the use of SIMD in the high performance templated classes.
 // #define MFEM_USE_SIMD

@@ -45,7 +45,6 @@ MFEM_USE_NETCDF        = @MFEM_USE_NETCDF@
 MFEM_USE_PETSC         = @MFEM_USE_PETSC@
 MFEM_USE_SLEPC         = @MFEM_USE_SLEPC@
 MFEM_USE_MPFR          = @MFEM_USE_MPFR@
-MFEM_USE_SIDRE         = @MFEM_USE_SIDRE@
 MFEM_USE_FMS           = @MFEM_USE_FMS@
 MFEM_USE_CONDUIT       = @MFEM_USE_CONDUIT@
 MFEM_USE_PUMI          = @MFEM_USE_PUMI@
@@ -48,7 +48,6 @@ option(MFEM_USE_NETCDF "Enable NETCDF usage" OFF)
 option(MFEM_USE_PETSC "Enable PETSc support." OFF)
 option(MFEM_USE_SLEPC "Enable SLEPc support." OFF)
 option(MFEM_USE_MPFR "Enable MPFR usage." OFF)
-option(MFEM_USE_SIDRE "Enable Axom/Sidre usage" OFF)
 option(MFEM_USE_FMS "Enable FMS usage" OFF)
 option(MFEM_USE_CONDUIT "Enable Conduit usage" OFF)
 option(MFEM_USE_PUMI "Enable PUMI" OFF)
@@ -224,17 +223,8 @@ set(FMS_DIR "${MFEM_DIR}/../fms" CACHE PATH
 set(CONDUIT_DIR "${MFEM_DIR}/../conduit" CACHE PATH
    "Path to the Conduit library.")

-set(AXOM_DIR "${MFEM_DIR}/../axom" CACHE PATH "Path to the Axom library.")
-# May need to add "Boost" as requirement.
-if (MFEM_USE_SIDRE)
-    if (MFEM_USE_MPI)
-        set(Axom_REQUIRED_PACKAGES "Conduit/blueprint/blueprint_mpi/relay/relay_mpi" CACHE STRING
-            "Additional packages required by Axom.")
-    elseif()
-        set(Axom_REQUIRED_PACKAGES "Conduit/blueprint/relay" CACHE STRING
-            "Additional packages required by Axom.")
-    endif()
-endif()
+set(AXOM_DIR "${MFEM_DIR}/../axom" CACHE PATH
+    "Path to the Axom library (required by Tribol for the contact mini-app).")

 set(PUMI_DIR "${MFEM_DIR}/../pumi-2.1.0" CACHE STRING
    "Directory where PUMI is installed")
@@ -252,6 +242,7 @@ set(MKL_PARDISO_DIR "" CACHE STRING "MKL installation path.")

 set(OCCA_DIR "${MFEM_DIR}/../occa" CACHE PATH "Path to OCCA")
 set(RAJA_DIR "${MFEM_DIR}/../raja" CACHE PATH "Path to RAJA")
+set(CAMP_DIR "${MFEM_DIR}/../camp" CACHE PATH "Path to CAMP (required by RAJA/Umpire)")
 set(CEED_DIR "${MFEM_DIR}/../libCEED" CACHE PATH "Path to libCEED")
 set(UMPIRE_DIR "${MFEM_DIR}/../umpire" CACHE PATH "Path to Umpire")
 set(CALIPER_DIR "${MFEM_DIR}/../caliper" CACHE PATH "Path to Caliper")
@@ -272,8 +263,12 @@ set(PARELAG_LIBRARIES "${PARELAG_DIR}/build/src/libParELAG.a" CACHE STRING
    "The ParELAG library.")

 set(TRIBOL_DIR "${MFEM_DIR}/../tribol" CACHE PATH "Path to Tribol")
-set(Tribol_REQUIRED_PACKAGES "Axom/core/mint/slam/slic" CACHE STRING
-    "Additional packages required by Tribol")
+# Tribol requires Axom. Many Tribol builds also enable optional TPLs like
+# RAJA/UMPIRE/Caliper, and may pull additional Axom components (e.g. quest,
+# lumberjack) via its exported targets.
+set(Tribol_REQUIRED_PACKAGES
+    "REQUIRED:;Axom/core/primal/mint/slam/slic/quest/lumberjack;OPTIONAL:;Adiak;CAMP;RAJA;UMPIRE;Caliper"
+    CACHE STRING "Additional packages required by Tribol")

 set(ENZYME_DIR "${MFEM_DIR}/../enzyme" CACHE PATH "Path to Enzyme")

@@ -162,7 +162,6 @@ MFEM_USE_NETCDF        = NO
 MFEM_USE_PETSC         = NO
 MFEM_USE_SLEPC         = NO
 MFEM_USE_MPFR          = NO
-MFEM_USE_SIDRE         = NO
 MFEM_USE_FMS           = NO
 MFEM_USE_CONDUIT       = NO
 MFEM_USE_PUMI          = NO
@@ -249,6 +248,15 @@ endif

 # METIS library configuration
 ifeq ($(MFEM_USE_SUPERLU)$(MFEM_USE_STRUMPACK)$(MFEM_USE_MUMPS),NONONO)
+   # MFEM_USE_METIS_5: when the user supplies METIS_DIR, try to auto-detect
+   # METIS 5 installs that follow the common <prefix>/{include,lib,lib64} layout.
+   ifeq ($(MFEM_USE_METIS_5),NO)
+      ifneq ($(wildcard $(METIS_DIR)/include/metis.h),)
+         ifneq ($(wildcard $(METIS_DIR)/lib/libmetis.* $(METIS_DIR)/lib64/libmetis.*),)
+            MFEM_USE_METIS_5 = YES
+         endif
+      endif
+   endif
   ifeq ($(MFEM_USE_METIS_5),NO)
     METIS_DIR = @MFEM_DIR@/../metis-4.0
     METIS_OPT =
@@ -487,17 +495,6 @@ ifneq (,$(wildcard $(CONDUIT_HDF5_HEADER)))
                  -lhdf5 $(ZLIB_LIB)
 endif

-# Sidre and required libraries configuration
-# Be sure to check the HDF5_DIR (set above) is correct
-SIDRE_DIR = @MFEM_DIR@/../axom
-SIDRE_OPT = -I$(SIDRE_DIR)/include -I$(CONDUIT_DIR)/include/conduit\
- -I$(HDF5_DIR)/include
-SIDRE_LIB = \
-   $(XLINKER)-rpath,$(SIDRE_DIR)/lib -L$(SIDRE_DIR)/lib \
-   $(XLINKER)-rpath,$(CONDUIT_DIR)/lib -L$(CONDUIT_DIR)/lib \
-   $(XLINKER)-rpath,$(HDF5_DIR)/lib -L$(HDF5_DIR)/lib \
-   -laxom -lconduit -lconduit_relay -lconduit_blueprint -lhdf5 $(ZLIB_LIB) -ldl
-
 # PUMI
 # Note that PUMI_DIR is needed -- it is used to check for gmi_sim.h
 PUMI_DIR = @MFEM_DIR@/../pumi-2.1.0
@@ -579,7 +576,13 @@ ifdef CUB_DIR
   RAJA_OPT += -I$(CUB_DIR)
 endif

+# CAMP library configuration (required by RAJA/Umpire for most installs)
 CAMP_LIB = -lcamp
+# If the common sibling layout exists, use it as a default (handles versioned
+# directories like camp-<hash>).
+ifneq ($(wildcard $(RAJA_DIR)/../camp*/include/camp/camp.hpp),)
+   CAMP_DIR ?= $(patsubst %/include/camp/camp.hpp,%,$(firstword $(wildcard $(RAJA_DIR)/../camp*/include/camp/camp.hpp)))
+endif
 ifdef CAMP_DIR
   RAJA_OPT += -I$(CAMP_DIR)/include
   CAMP_LIB = $(XLINKER)-rpath,$(CAMP_DIR)/lib -L$(CAMP_DIR)/lib -lcamp
@@ -589,7 +592,12 @@ RAJA_LIB = $(XLINKER)-rpath,$(RAJA_DIR)/lib -L$(RAJA_DIR)/lib -lRAJA $(CAMP_LIB)
 # UMPIRE library configuration
 UMPIRE_DIR = @MFEM_DIR@/../umpire
 UMPIRE_OPT = -I$(UMPIRE_DIR)/include $(if $(CAMP_DIR), -I$(CAMP_DIR)/include)
-UMPIRE_LIB = -L$(UMPIRE_DIR)/lib -L$(UMPIRE_DIR)/lib64 -lumpire $(CAMP_LIB)
+UMPIRE_LIB = -L$(UMPIRE_DIR)/lib -L$(UMPIRE_DIR)/lib64 -lumpire $(CAMP_LIB) -lpthread
+# If the common sibling layout exists, use it as a default (handles versioned
+# directories like fmt-<hash>).
+ifneq ($(wildcard $(UMPIRE_DIR)/../fmt*/include/fmt/format.h),)
+   FMT_DIR ?= $(patsubst %/include/fmt/format.h,%,$(firstword $(wildcard $(UMPIRE_DIR)/../fmt*/include/fmt/format.h)))
+endif
 ifdef FMT_DIR
   UMPIRE_OPT += -I$(FMT_DIR)/include
   UMPIRE_LIB += -L$(FMT_DIR)/lib -L$(FMT_DIR)/lib64 -lfmt
@@ -621,8 +629,63 @@ PARELAG_LIB = -L$(PARELAG_DIR)/build/src -lParELAG
 AXOM_DIR = @MFEM_DIR@/../axom
 TRIBOL_DIR = @MFEM_DIR@/../tribol
 TRIBOL_OPT = -I$(TRIBOL_DIR)/include -I$(AXOM_DIR)/include
-TRIBOL_LIB = -L$(TRIBOL_DIR)/lib -ltribol -lredecomp -L$(AXOM_DIR)/lib -laxom_mint\
-   -laxom_slam -laxom_slic -laxom_core
+# Tribol may be built with optional dependencies (e.g. RAJA/UMPIRE/CALIPER).
+# Add those options only when the corresponding headers/libraries exist.
+ifneq ($(wildcard $(RAJA_DIR)/include/RAJA/RAJA.hpp),)
+   TRIBOL_OPT += $(RAJA_OPT)
+endif
+ifneq ($(wildcard $(UMPIRE_DIR)/include/umpire/Umpire.hpp),)
+   TRIBOL_OPT += $(UMPIRE_OPT)
+endif
+ifneq ($(wildcard $(CALIPER_DIR)/include/caliper/cali.h),)
+   TRIBOL_OPT += $(CALIPER_OPT)
+endif
+
+TRIBOL_LIB = -L$(TRIBOL_DIR)/lib -L$(TRIBOL_DIR)/lib64
+ifneq ($(wildcard $(TRIBOL_DIR)/lib/libtribol.* $(TRIBOL_DIR)/lib64/libtribol.*),)
+   TRIBOL_LIB += -ltribol
+endif
+ifneq ($(wildcard $(TRIBOL_DIR)/lib/libtribol_shared.* $(TRIBOL_DIR)/lib64/libtribol_shared.*),)
+   TRIBOL_LIB += -ltribol_shared
+endif
+ifneq ($(wildcard $(TRIBOL_DIR)/lib/libredecomp.* $(TRIBOL_DIR)/lib64/libredecomp.*),)
+   TRIBOL_LIB += -lredecomp
+endif
+
+TRIBOL_LIB += -L$(AXOM_DIR)/lib -L$(AXOM_DIR)/lib64
+ifneq ($(wildcard $(AXOM_DIR)/lib/libaxom_quest.* $(AXOM_DIR)/lib64/libaxom_quest.*),)
+   TRIBOL_LIB += -laxom_quest
+endif
+ifneq ($(wildcard $(AXOM_DIR)/lib/libaxom_mint.* $(AXOM_DIR)/lib64/libaxom_mint.*),)
+   TRIBOL_LIB += -laxom_mint
+endif
+ifneq ($(wildcard $(AXOM_DIR)/lib/libaxom_slam.* $(AXOM_DIR)/lib64/libaxom_slam.*),)
+   TRIBOL_LIB += -laxom_slam
+endif
+ifneq ($(wildcard $(AXOM_DIR)/lib/libaxom_slic.* $(AXOM_DIR)/lib64/libaxom_slic.*),)
+   TRIBOL_LIB += -laxom_slic
+endif
+ifneq ($(wildcard $(AXOM_DIR)/lib/libaxom_lumberjack.* $(AXOM_DIR)/lib64/libaxom_lumberjack.*),)
+   TRIBOL_LIB += -laxom_lumberjack
+endif
+ifneq ($(wildcard $(AXOM_DIR)/lib/libaxom_core.* $(AXOM_DIR)/lib64/libaxom_core.*),)
+   TRIBOL_LIB += -laxom_core
+endif
+
+# Add common optional Tribol TPLs when their libraries are present.
+ifneq ($(wildcard $(ADIAK_DIR)/lib/libadiak.* $(ADIAK_DIR)/lib64/libadiak.*),)
+   TRIBOL_LIB += $(XLINKER)-rpath,$(ADIAK_DIR)/lib64 $(XLINKER)-rpath,$(ADIAK_DIR)/lib \
+      -L$(ADIAK_DIR)/lib64 -L$(ADIAK_DIR)/lib -ladiak -ldl
+endif
+ifneq ($(wildcard $(UMPIRE_DIR)/lib/libumpire.* $(UMPIRE_DIR)/lib64/libumpire.*),)
+   TRIBOL_LIB += $(UMPIRE_LIB)
+endif
+ifneq ($(wildcard $(RAJA_DIR)/lib/libRAJA.* $(RAJA_DIR)/lib64/libRAJA.*),)
+   TRIBOL_LIB += $(RAJA_LIB)
+endif
+ifneq ($(wildcard $(CALIPER_DIR)/lib/libcaliper.* $(CALIPER_DIR)/lib64/libcaliper.*),)
+   TRIBOL_LIB += $(CALIPER_LIB)
+endif

 # Enzyme configuration
 ENZYME_DIR = @MFEM_DIR@/../enzyme
@@ -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()
@@ -64,7 +64,7 @@ PARALLEL_NAME := Parallel AMGX example
 $(MFEM_LIB_FILE):
 	$(error The MFEM library is not build)

-clean: clean-build clean-exec
+clean: clean-build

 clean-build:
 	rm -f *.o *~ $(SEQ_EXAMPLES) $(PAR_EXAMPLES)
@@ -64,12 +64,12 @@ ex1p-test-par: ex1p
 $(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_EXAMPLES) $(PAR_EXAMPLES)
 	rm -rf *.dSYM *.TVD.*breakpoints

 clean-exec:
-	@rm -f refined.mesh mesh.*
-	@rm -f sol.*
+	@rm -f refined.mesh displaced.mesh mesh.* ex5.mesh
+	@rm -f sphere_refined.* sol.* sol_u.* sol_p.* sol_r.* sol_i.*
@@ -5,8 +5,8 @@
 // Sample runs:
 //     ex37 -alpha 10
 //     ex37 -alpha 10 -pv
-//     ex37 -lambda 0.1 -mu 0.1 -growth 1
-//     ex37 -o 2 -alpha 10.0 -mi 50 -vf 0.4 -ntol 1e-5 -growth 1.5
+//     ex37 -lambda 0.1 -mu 0.1
+//     ex37 -o 2 -alpha 5.0 -mi 50 -vf 0.4 -ntol 1e-5
 //     ex37 -r 6 -o 1 -alpha 25.0 -epsilon 0.02 -mi 50 -ntol 1e-5
 //
 // Description: This example code demonstrates the use of MFEM to solve a
@@ -55,6 +55,53 @@
 using namespace std;
 using namespace mfem;

+/**
+ * @brief Bregman projection of ρ = sigmoid(ψ) onto the subspace
+ *        ∫_Ω ρ dx = θ vol(Ω) as follows:
+ *
+ *        1. Compute the root of the R → R function
+ *            f(c) = ∫_Ω sigmoid(ψ + c) dx - θ vol(Ω)
+ *        2. Set ψ ← ψ + c.
+ *
+ * @param psi a GridFunction to be updated
+ * @param target_volume θ vol(Ω)
+ * @param tol Newton iteration tolerance
+ * @param max_its Newton maximum iteration number
+ * @return real_t Final volume, ∫_Ω sigmoid(ψ)
+ */
+real_t proj(GridFunction &psi, real_t target_volume, real_t tol=1e-12,
+            int max_its=10)
+{
+   MappedGridFunctionCoefficient sigmoid_psi(&psi, sigmoid);
+   MappedGridFunctionCoefficient der_sigmoid_psi(&psi, der_sigmoid);
+
+   LinearForm int_sigmoid_psi(psi.FESpace());
+   int_sigmoid_psi.AddDomainIntegrator(new DomainLFIntegrator(sigmoid_psi));
+   LinearForm int_der_sigmoid_psi(psi.FESpace());
+   int_der_sigmoid_psi.AddDomainIntegrator(new DomainLFIntegrator(
+                                              der_sigmoid_psi));
+   bool done = false;
+   for (int k=0; k<max_its; k++) // Newton iteration
+   {
+      int_sigmoid_psi.Assemble(); // Recompute f(c) with updated ψ
+      const real_t f = int_sigmoid_psi.Sum() - target_volume;
+
+      int_der_sigmoid_psi.Assemble(); // Recompute df(c) with updated ψ
+      const real_t df = int_der_sigmoid_psi.Sum();
+
+      const real_t dc = -f/df;
+      psi += dc;
+      if (abs(dc) < tol) { done = true; break; }
+   }
+   if (!done)
+   {
+      mfem_warning("Projection reached maximum iteration without converging. "
+                   "Result may not be accurate.");
+   }
+   int_sigmoid_psi.Assemble();
+   return int_sigmoid_psi.Sum();
+}
+
 /*
 * ---------------------------------------------------------------
 *                      ALGORITHM PREAMBLE
@@ -133,11 +180,10 @@ int main(int argc, char *argv[])
   int ref_levels = 5;
   int order = 2;
   real_t alpha = 1.0;
-   real_t growth = 2;
   real_t epsilon = 0.01;
   real_t vol_fraction = 0.5;
   int max_it = 1e3;
-   real_t itol = 1e-2;
+   real_t itol = 1e-1;
   real_t ntol = 1e-4;
   real_t rho_min = 1e-6;
   real_t lambda = 1.0;
@@ -152,8 +198,6 @@ int main(int argc, char *argv[])
                  "Order (degree) of the finite elements.");
   args.AddOption(&alpha, "-alpha", "--alpha-step-length",
                  "Step length for gradient descent.");
-   args.AddOption(&growth, "-growth", "--alpha-growth-rate",
-                  "Growth rate of step length for gradient descent.");
   args.AddOption(&epsilon, "-epsilon", "--epsilon-thickness",
                  "Length scale for ρ.");
   args.AddOption(&max_it, "-mi", "--max-it",
@@ -288,7 +332,6 @@ int main(int argc, char *argv[])
   }
   FilterSolver->SetEssentialBoundary(ess_bdr_filter);
   FilterSolver->SetupFEM();
-   FilterSolver->AssembleDiffusionBilinear();

   BilinearForm mass(&control_fes);
   mass.AddDomainIntegrator(new InverseIntegrator(new MassIntegrator(one)));
@@ -342,7 +385,7 @@ int main(int argc, char *argv[])
   // 11. Iterate:
   for (int k = 1; k <= max_it; k++)
   {
-      if (k > 1) { alpha = std::pow((real_t) k,growth); }
+      if (k > 1) { alpha *= ((real_t) k) / ((real_t) k-1); }

      mfem::out << "\nStep = " << k << std::endl;

@@ -379,9 +422,7 @@ int main(int argc, char *argv[])

      // Step 5 - Update design variable ψ ← proj(ψ - αG)
      psi.Add(-alpha, grad);
-      GridFunction alpha_grad(grad);
-      alpha_grad *= alpha;
-      const real_t material_volume = proj(psi, alpha_grad, target_volume);
+      const real_t material_volume = proj(psi, target_volume);

      // Compute ||ρ - ρ_old|| in control fes.
      real_t norm_increment = zerogf.ComputeL1Error(succ_diff_rho);
@@ -137,7 +137,7 @@ public:
        exponent(exponent_), rho_min(rho_min_)
   {
      MFEM_ASSERT(rho_min_ >= 0.0, "rho_min must be >= 0");
-      MFEM_ASSERT(rho_min_ < 1.0,  "rho_min must be < 1");
+      MFEM_ASSERT(rho_min_ < 1.0,  "rho_min must be > 1");
      MFEM_ASSERT(u, "displacement field is not set");
      MFEM_ASSERT(rho_filter, "density field is not set");
   }
@@ -231,12 +231,9 @@ private:
   FiniteElementCollection * fec = nullptr;
   FiniteElementSpace * fes = nullptr;
   Array<int> ess_bdr;
-   Array<int> ess_tdof_list;
   Array<int> neumann_bdr;
   GridFunction * u = nullptr;
   LinearForm * b = nullptr;
-   BilinearForm * a = nullptr;
-   OperatorPtr A;
   bool parallel;
 #ifdef MFEM_USE_MPI
   ParMesh * pmesh = nullptr;
@@ -270,8 +267,6 @@ public:
   void ResetFEM();
   void SetupFEM();

-   void UpdateEssentialTDofs();
-   void AssembleDiffusionBilinear(bool update_ess_tdofs=true);
   void Solve();
   GridFunction * GetFEMSolution();
   LinearForm * GetLinearForm() {return b;}
@@ -376,130 +371,6 @@ public:

 };

-/**
- * @brief Bregman projection of ρ = sigmoid(ψ) onto the subspace
- *        ∫_Ω ρ dx = θ vol(Ω) as follows:
- *
- *        1. Compute the root of the R → R function
- *            f(c) = ∫_Ω sigmoid(ψ + c) dx - θ vol(Ω)
- *           using the Illinois method
- *        2. Set ψ ← ψ + c.
- *
- * @param psi a GridFunction to be updated
- * @param alpha_grad alpha multiplied by gradient
- * @param target_volume θ vol(Ω)
- * @param tol Illinois iteration tolerance
- * @param max_its Illinois maximum iteration number
- * @return real_t Final volume (∫_Ω sigmoid(ψ) dx)
- */
-real_t proj(GridFunction &psi, GridFunction &alpha_grad, real_t target_volume,
-            real_t tol = 1e-12, int max_its = 100)
-{
-#ifdef MFEM_USE_MPI
-   FiniteElementSpace *fes = psi.FESpace();
-   ParFiniteElementSpace *pfes = dynamic_cast<ParFiniteElementSpace*>(fes);
-#endif
-   ConstantCoefficient zero_cf(0.0);
-   real_t a = -alpha_grad.ComputeMaxError(zero_cf);
-   real_t b = -a;
-   real_t y = 0.0;
-
-   MappedGridFunctionCoefficient sigmoid_psi(
-   &psi, [&y](const real_t x) { return sigmoid(x + y); });
-   std::unique_ptr<LinearForm> int_sigmoid_psi;
-#ifdef MFEM_USE_MPI
-   ParGridFunction *par_psi = dynamic_cast<ParGridFunction *>(&psi);
-   if (par_psi)
-   {
-      int_sigmoid_psi.reset(new ParLinearForm(par_psi->ParFESpace()));
-   }
-   else
-   {
-      int_sigmoid_psi.reset(new LinearForm(psi.FESpace()));
-   }
-#else
-   int_sigmoid_psi.reset(new LinearForm(psi.FESpace()));
-#endif
-   int_sigmoid_psi->AddDomainIntegrator(new DomainLFIntegrator(sigmoid_psi));
-
-   y = a;
-   int_sigmoid_psi->Assemble();
-   real_t f_a = int_sigmoid_psi->Sum();  // f_a := f(a) + θ vol(Ω)
-
-   y = b;
-   int_sigmoid_psi->Assemble();
-   real_t f_b = int_sigmoid_psi->Sum();  // f_b := f(b) + θ vol(Ω)
-#ifdef MFEM_USE_MPI
-   if (pfes)
-   {
-      MPI_Allreduce(MPI_IN_PLACE, &f_a, 1, MPITypeMap<real_t>::mpi_type,
-                    MPI_SUM, MPI_COMM_WORLD);
-      MPI_Allreduce(MPI_IN_PLACE, &f_b, 1, MPITypeMap<real_t>::mpi_type,
-                    MPI_SUM, MPI_COMM_WORLD);
-   }
-#endif
-   f_a -= target_volume; // f_a := f(a)
-   f_b -= target_volume; // f_b := f(b)
-   real_t c = 0.0;
-   real_t f_c = 0.0;
-   int side = 0;
-
-   bool done = false;
-   for (int k=0; k < max_its; k++)
-   {
-      c = (f_a * b - f_b * a) / (f_a - f_b);
-
-      if (abs(b - a) < tol * abs(b + a)) { done = true; break; }
-
-      y = c;
-      int_sigmoid_psi->Assemble();
-      f_c = int_sigmoid_psi->Sum(); // f_c := f(c) + θ vol(Ω)
-#ifdef MFEM_USE_MPI
-      if (pfes)
-      {
-         MPI_Allreduce(MPI_IN_PLACE, &f_c, 1, MPITypeMap<real_t>::mpi_type,
-                       MPI_SUM, MPI_COMM_WORLD);
-      }
-#endif
-      f_c -= target_volume; // f_c := f(c)
-
-      if (f_c * f_b > 0)
-      {
-         b = c;
-         f_b = f_c;
-         if (side == -1) { f_a /= 2.0; }
-         side = -1;
-      }
-      else if (f_c * f_a > 0)
-      {
-         a = c;
-         f_a = f_c;
-         if (side == 1) { f_b /= 2.0; }
-         side = 1;
-      }
-      else
-      {
-         done = true; break;
-      }
-   }
-   if (!done)
-   {
-      mfem_warning("Projection reached maximum iteration without converging. "
-                   "Result may not be accurate.");
-   }
-   y = 0.0;
-   psi += c;
-   int_sigmoid_psi->Assemble();
-   real_t material_volume = int_sigmoid_psi->Sum();
-#ifdef MFEM_USE_MPI
-   if (pfes)
-   {
-      MPI_Allreduce(MPI_IN_PLACE, &material_volume, 1,
-                    MPITypeMap<real_t>::mpi_type, MPI_SUM, MPI_COMM_WORLD);
-   }
-#endif
-   return material_volume;
-}

 // Poisson solver

@@ -551,8 +422,12 @@ void DiffusionSolver::SetupFEM()
   }
 }

-void DiffusionSolver::UpdateEssentialTDofs()
+void DiffusionSolver::Solve()
 {
+   OperatorPtr A;
+   Vector B, X;
+   Array<int> ess_tdof_list;
+
 #ifdef MFEM_USE_MPI
   if (parallel)
   {
@@ -565,39 +440,7 @@ void DiffusionSolver::UpdateEssentialTDofs()
 #else
   fes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list);
 #endif
-}
-
-void DiffusionSolver::AssembleDiffusionBilinear(bool update_ess_tdofs)
-{
-   if (update_ess_tdofs)
-   {
-      UpdateEssentialTDofs();
-   }
-#ifdef MFEM_USE_MPI
-   if (parallel)
-   {
-      a = new ParBilinearForm(pfes);
-   }
-   else
-   {
-      a = new BilinearForm(fes);
-   }
-#else
-   a = new BilinearForm(fes);
-#endif
-   a->AddDomainIntegrator(new DiffusionIntegrator(*diffcf));
-   if (masscf)
-   {
-      a->AddDomainIntegrator(new MassIntegrator(*masscf));
-   }
-   a->Assemble();
-   a->FormSystemMatrix(ess_tdof_list, A);
-}
-
-void DiffusionSolver::Solve()
-{
-   Vector B, X;
-
+   *u=0.0;
   if (b)
   {
      delete b;
@@ -632,33 +475,31 @@ void DiffusionSolver::Solve()

   b->Assemble();

-   *u=0.0;
-   if (essbdr_cf)
-   {
-      u->ProjectBdrCoefficient(*essbdr_cf,ess_bdr);
-   }
+   BilinearForm * a = nullptr;

 #ifdef MFEM_USE_MPI
   if (parallel)
   {
-      X.SetSize(pfes->TrueVSize());
-      B.SetSize(pfes->TrueVSize());
-      dynamic_cast<ParGridFunction*>(u)->ParallelAssemble(X);
-      dynamic_cast<ParLinearForm*>(b)->ParallelAssemble(B);
-      dynamic_cast<ParBilinearForm*>(a)->ParallelEliminateTDofsInRHS(
-         ess_tdof_list, X, B);
+      a = new ParBilinearForm(pfes);
   }
   else
   {
-      X.NewDataAndSize(u->GetData(), u->Size());
-      B.NewDataAndSize(b->GetData(), b->Size());
-      a->EliminateVDofsInRHS(ess_tdof_list, X, B);
+      a = new BilinearForm(fes);
   }
 #else
-   X.NewDataAndSize(u->GetData(), u->Size());
-   B.NewDataAndSize(b->GetData(), b->Size());
-   a->EliminateVDofsInRHS(ess_tdof_list, X, B);
+   a = new BilinearForm(fes);
 #endif
+   a->AddDomainIntegrator(new DiffusionIntegrator(*diffcf));
+   if (masscf)
+   {
+      a->AddDomainIntegrator(new MassIntegrator(*masscf));
+   }
+   a->Assemble();
+   if (essbdr_cf)
+   {
+      u->ProjectBdrCoefficient(*essbdr_cf,ess_bdr);
+   }
+   a->FormLinearSystem(ess_tdof_list, *u, *b, A, X, B);

   CGSolver * cg = nullptr;
   Solver * M = nullptr;
@@ -687,6 +528,7 @@ void DiffusionSolver::Solve()
   delete M;
   delete cg;
   a->RecoverFEMSolution(X, *b, *u);
+   delete a;
 }

 GridFunction * DiffusionSolver::GetFEMSolution()
@@ -718,8 +560,6 @@ DiffusionSolver::~DiffusionSolver()
 #endif
   delete fec; fec = nullptr;
   delete b;
-   A.Clear();
-   delete a;
 }


@@ -4,8 +4,8 @@
 //
 // Sample runs:
 //    mpirun -np 4 ex37p -alpha 10 -pv
-//    mpirun -np 4 ex37p -lambda 0.1 -mu 0.1 -growth 1
-//    mpirun -np 4 ex37p -o 2 -alpha 10.0 -mi 50 -vf 0.4 -ntol 1e-5 -growth 1.5
+//    mpirun -np 4 ex37p -lambda 0.1 -mu 0.1
+//    mpirun -np 4 ex37p -o 2 -alpha 5.0 -mi 50 -vf 0.4 -ntol 1e-5
 //    mpirun -np 4 ex37p -r 6 -o 2 -alpha 10.0 -epsilon 0.02 -mi 50 -ntol 1e-5
 //
 // Description: This example code demonstrates the use of MFEM to solve a
@@ -54,6 +54,61 @@
 using namespace std;
 using namespace mfem;

+/**
+ * @brief Bregman projection of ρ = sigmoid(ψ) onto the subspace
+ *        ∫_Ω ρ dx = θ vol(Ω) as follows:
+ *
+ *        1. Compute the root of the R → R function
+ *            f(c) = ∫_Ω sigmoid(ψ + c) dx - θ vol(Ω)
+ *        2. Set ψ ← ψ + c.
+ *
+ * @param psi a GridFunction to be updated
+ * @param target_volume θ vol(Ω)
+ * @param tol Newton iteration tolerance
+ * @param max_its Newton maximum iteration number
+ * @return real_t Final volume, ∫_Ω sigmoid(ψ)
+ */
+real_t proj(ParGridFunction &psi, real_t target_volume, real_t tol=1e-12,
+            int max_its=10)
+{
+   MappedGridFunctionCoefficient sigmoid_psi(&psi, sigmoid);
+   MappedGridFunctionCoefficient der_sigmoid_psi(&psi, der_sigmoid);
+
+   ParLinearForm int_sigmoid_psi(psi.ParFESpace());
+   int_sigmoid_psi.AddDomainIntegrator(new DomainLFIntegrator(sigmoid_psi));
+   ParLinearForm int_der_sigmoid_psi(psi.ParFESpace());
+   int_der_sigmoid_psi.AddDomainIntegrator(new DomainLFIntegrator(
+                                              der_sigmoid_psi));
+   bool done = false;
+   for (int k=0; k<max_its; k++) // Newton iteration
+   {
+      int_sigmoid_psi.Assemble(); // Recompute f(c) with updated ψ
+      real_t f = int_sigmoid_psi.Sum();
+      MPI_Allreduce(MPI_IN_PLACE, &f, 1, MPITypeMap<real_t>::mpi_type,
+                    MPI_SUM, MPI_COMM_WORLD);
+      f -= target_volume;
+
+      int_der_sigmoid_psi.Assemble(); // Recompute df(c) with updated ψ
+      real_t df = int_der_sigmoid_psi.Sum();
+      MPI_Allreduce(MPI_IN_PLACE, &df, 1, MPITypeMap<real_t>::mpi_type,
+                    MPI_SUM, MPI_COMM_WORLD);
+
+      const real_t dc = -f/df;
+      psi += dc;
+      if (abs(dc) < tol) { done = true; break; }
+   }
+   if (!done)
+   {
+      mfem_warning("Projection reached maximum iteration without converging. "
+                   "Result may not be accurate.");
+   }
+   int_sigmoid_psi.Assemble();
+   real_t material_volume = int_sigmoid_psi.Sum();
+   MPI_Allreduce(MPI_IN_PLACE, &material_volume, 1,
+                 MPITypeMap<real_t>::mpi_type, MPI_SUM, MPI_COMM_WORLD);
+   return material_volume;
+}
+
 /*
 * ---------------------------------------------------------------
 *                      ALGORITHM PREAMBLE
@@ -138,11 +193,10 @@ int main(int argc, char *argv[])
   int ref_levels = 5;
   int order = 2;
   real_t alpha = 1.0;
-   real_t growth = 2;
   real_t epsilon = 0.01;
   real_t vol_fraction = 0.5;
   int max_it = 1e3;
-   real_t itol = 1e-2;
+   real_t itol = 1e-1;
   real_t ntol = 1e-4;
   real_t rho_min = 1e-6;
   real_t lambda = 1.0;
@@ -157,8 +211,6 @@ int main(int argc, char *argv[])
                  "Order (degree) of the finite elements.");
   args.AddOption(&alpha, "-alpha", "--alpha-step-length",
                  "Step length for gradient descent.");
-   args.AddOption(&growth, "-growth", "--alpha-growth-rate",
-                  "Growth rate of step length for gradient descent.");
   args.AddOption(&epsilon, "-epsilon", "--epsilon-thickness",
                  "Length scale for ρ.");
   args.AddOption(&max_it, "-mi", "--max-it",
@@ -307,7 +359,6 @@ int main(int argc, char *argv[])
   }
   FilterSolver->SetEssentialBoundary(ess_bdr_filter);
   FilterSolver->SetupFEM();
-   FilterSolver->AssembleDiffusionBilinear();

   ParBilinearForm mass(&control_fes);
   mass.AddDomainIntegrator(new InverseIntegrator(new MassIntegrator(one)));
@@ -361,7 +412,7 @@ int main(int argc, char *argv[])
   // 11. Iterate:
   for (int k = 1; k <= max_it; k++)
   {
-      if (k > 1) { alpha = std::pow((real_t) k,growth); }
+      if (k > 1) { alpha *= ((real_t) k) / ((real_t) k-1); }

      if (myid == 0)
      {
@@ -401,9 +452,7 @@ int main(int argc, char *argv[])

      // Step 5 - Update design variable ψ ← proj(ψ - αG)
      psi.Add(-alpha, grad);
-      ParGridFunction alpha_grad(grad);
-      alpha_grad *= alpha;
-      const real_t material_volume = proj(psi, alpha_grad, target_volume);
+      const real_t material_volume = proj(psi, target_volume);

      // Compute ||ρ - ρ_old|| in control fes.
      real_t norm_increment = zerogf.ComputeL1Error(succ_diff_rho);
@@ -433,16 +433,18 @@ int main(int argc, char *argv[])
   u.ProjectCoefficient(*u0);

   // Create data collection for solution output: either VisItDataCollection for
-   // ascii data files, or SidreDataCollection for binary data files.
+   // ascii data files, or ConduitDataCollection for binary data files.
   DataCollection *dc = NULL;
   if (visit)
   {
      if (binary)
      {
-#ifdef MFEM_USE_SIDRE
-         dc = new SidreDataCollection("Example41", &mesh);
+#ifdef MFEM_USE_CONDUIT
+         auto conduit_dc = new ConduitDataCollection("Example41", &mesh);
+         conduit_dc->SetProtocol("hdf5");
+         dc = conduit_dc;
 #else
-         MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
+         MFEM_ABORT("Must build with MFEM_USE_CONDUIT=YES for binary output.");
 #endif
      }
      else
@@ -518,10 +518,12 @@ int main(int argc, char *argv[])
   {
      if (binary)
      {
-#ifdef MFEM_USE_SIDRE
-         dc = new SidreDataCollection("Example41-Parallel", pmesh);
+#ifdef MFEM_USE_CONDUIT
+         auto conduit_dc = new ConduitDataCollection("Example41-Parallel", pmesh);
+         conduit_dc->SetProtocol("hdf5");
+         dc = conduit_dc;
 #else
-         MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
+         MFEM_ABORT("Must build with MFEM_USE_CONDUIT=YES for binary output.");
 #endif
      }
      else
@@ -305,16 +305,18 @@ int main(int argc, char *argv[])
   }

   // Create data collection for solution output: either VisItDataCollection for
-   // ascii data files, or SidreDataCollection for binary data files.
+   // ascii data files, or ConduitDataCollection for binary data files.
   DataCollection *dc = NULL;
   if (visit)
   {
      if (binary)
      {
-#ifdef MFEM_USE_SIDRE
-         dc = new SidreDataCollection("Example9", &mesh);
+#ifdef MFEM_USE_CONDUIT
+         auto conduit_dc = new ConduitDataCollection("Example9", &mesh);
+         conduit_dc->SetProtocol("hdf5");
+         dc = conduit_dc;
 #else
-         MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
+         MFEM_ABORT("Must build with MFEM_USE_CONDUIT=YES for binary output.");
 #endif
      }
      else
@@ -441,16 +441,18 @@ int main(int argc, char *argv[])
   }

   // Create data collection for solution output: either VisItDataCollection for
-   // ascii data files, or SidreDataCollection for binary data files.
+   // ascii data files, or ConduitDataCollection for binary data files.
   DataCollection *dc = NULL;
   if (visit)
   {
      if (binary)
      {
-#ifdef MFEM_USE_SIDRE
-         dc = new SidreDataCollection("Example9-Parallel", pmesh);
+#ifdef MFEM_USE_CONDUIT
+         auto conduit_dc = new ConduitDataCollection("Example9-Parallel", pmesh);
+         conduit_dc->SetProtocol("hdf5");
+         dc = conduit_dc;
 #else
-         MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
+         MFEM_ABORT("Must build with MFEM_USE_CONDUIT=YES for binary output.");
 #endif
      }
      else
@@ -76,4 +76,4 @@ clean-build:
 	rm -rf *.dSYM *.TVD.*breakpoints

 clean-exec:
-	@rm -f refined.mesh sol.gf mesh.* sol.*
+	@rm -f refined.mesh sol.gf
@@ -354,16 +354,18 @@ int main(int argc, char *argv[])
   }

   // Create data collection for solution output: either VisItDataCollection for
-   // ascii data files, or SidreDataCollection for binary data files.
+   // ascii data files, or ConduitDataCollection for binary data files.
   DataCollection *dc = NULL;
   if (visit)
   {
      if (binary)
      {
-#ifdef MFEM_USE_SIDRE
-         dc = new SidreDataCollection("Example9", mesh);
+#ifdef MFEM_USE_CONDUIT
+         auto conduit_dc = new ConduitDataCollection("Example9", mesh);
+         conduit_dc->SetProtocol("hdf5");
+         dc = conduit_dc;
 #else
-         MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
+         MFEM_ABORT("Must build with MFEM_USE_CONDUIT=YES for binary output.");
 #endif
      }
      else
@@ -414,16 +414,18 @@ int main(int argc, char *argv[])
   }

   // Create data collection for solution output: either VisItDataCollection for
-   // ascii data files, or SidreDataCollection for binary data files.
+   // ascii data files, or ConduitDataCollection for binary data files.
   DataCollection *dc = NULL;
   if (visit)
   {
      if (binary)
      {
-#ifdef MFEM_USE_SIDRE
-         dc = new SidreDataCollection("Example9-Parallel", pmesh);
+#ifdef MFEM_USE_CONDUIT
+         auto conduit_dc = new ConduitDataCollection("Example9-Parallel", pmesh);
+         conduit_dc->SetProtocol("hdf5");
+         dc = conduit_dc;
 #else
-         MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
+         MFEM_ABORT("Must build with MFEM_USE_CONDUIT=YES for binary output.");
 #endif
      }
      else
@@ -71,7 +71,6 @@ 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))

@@ -88,9 +87,8 @@ SUBDIRS_TPRINT = $(addsuffix /test-print,$(SUBDIRS))

 all: $(EXAMPLES) $(SUBDIRS_ALL)

-.PHONY: $(SUBDIRS_ALL) $(SUBDIRS_TEST) $(SUBDIRS_TEST_NOCLEAN) \
-   $(SUBDIRS_CLEAN) $(SUBDIRS_TPRINT)
-$(SUBDIRS_ALL) $(SUBDIRS_TEST) $(SUBDIRS_TEST_NOCLEAN) $(SUBDIRS_CLEAN):
+.PHONY: $(SUBDIRS_ALL) $(SUBDIRS_TEST) $(SUBDIRS_CLEAN) $(SUBDIRS_TPRINT)
+$(SUBDIRS_ALL) $(SUBDIRS_TEST) $(SUBDIRS_CLEAN):
 	$(MAKE) -C $(@D) $(@F)
 $(SUBDIRS_TPRINT):
 	@$(MAKE) -C $(@D) $(@F)
@@ -109,7 +107,6 @@ 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
@@ -160,8 +157,6 @@ ex37-test-seq: ex37
 	@$(call mfem-test,$<,, Serial example,-mi 3)
 ex37p-test-par: ex37p
 	@$(call mfem-test,$<, $(RUN_MPI), Parallel example,-mi 3)
-ex39-test-seq: ex39
-	@$(call mfem-test,$<,, Serial example,-m ../data/compass.mesh)
 ex41-test-seq: ex41
 	@$(call mfem-test,$<,, Serial example,-tf 1.0)
 ex41p-test-par: ex41p
@@ -368,16 +368,18 @@ int main(int argc, char *argv[])
   }

   // Create data collection for solution output: either VisItDataCollection for
-   // ascii data files, or SidreDataCollection for binary data files.
+   // ascii data files, or ConduitDataCollection for binary data files.
   DataCollection *dc = NULL;
   if (visit)
   {
      if (binary)
      {
-#ifdef MFEM_USE_SIDRE
-         dc = new SidreDataCollection("Example9-Parallel", pmesh);
+#ifdef MFEM_USE_CONDUIT
+         auto conduit_dc = new ConduitDataCollection("Example9-Parallel", pmesh);
+         conduit_dc->SetProtocol("hdf5");
+         dc = conduit_dc;
 #else
-         MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
+         MFEM_ABORT("Must build with MFEM_USE_CONDUIT=YES for binary output.");
 #endif
      }
      else
@@ -316,16 +316,18 @@ int main(int argc, char *argv[])
   }

   // Create data collection for solution output: either VisItDataCollection for
-   // ascii data files, or SidreDataCollection for binary data files.
+   // ascii data files, or ConduitDataCollection for binary data files.
   DataCollection *dc = NULL;
   if (visit)
   {
      if (binary)
      {
-#ifdef MFEM_USE_SIDRE
-         dc = new SidreDataCollection("Example9", &mesh);
+#ifdef MFEM_USE_CONDUIT
+         auto conduit_dc = new ConduitDataCollection("Example9", &mesh);
+         conduit_dc->SetProtocol("hdf5");
+         dc = conduit_dc;
 #else
-         MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
+         MFEM_ABORT("Must build with MFEM_USE_CONDUIT=YES for binary output.");
 #endif
      }
      else
@@ -453,16 +453,18 @@ int main(int argc, char *argv[])
   }

   // Create data collection for solution output: either VisItDataCollection for
-   // ascii data files, or SidreDataCollection for binary data files.
+   // ascii data files, or ConduitDataCollection for binary data files.
   DataCollection *dc = NULL;
   if (visit)
   {
      if (binary)
      {
-#ifdef MFEM_USE_SIDRE
-         dc = new SidreDataCollection("Example9-Parallel", pmesh);
+#ifdef MFEM_USE_CONDUIT
+         auto conduit_dc = new ConduitDataCollection("Example9-Parallel", pmesh);
+         conduit_dc->SetProtocol("hdf5");
+         dc = conduit_dc;
 #else
-         MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
+         MFEM_ABORT("Must build with MFEM_USE_CONDUIT=YES for binary output.");
 #endif
      }
      else
@@ -312,11 +312,6 @@ set(HDRS
  particleset.hpp
  )

-if (MFEM_USE_SIDRE)
-  list(APPEND SRCS sidredatacollection.cpp)
-  list(APPEND HDRS sidredatacollection.hpp)
-endif()
-
 if (MFEM_USE_CONDUIT)
  list(APPEND SRCS conduitdatacollection.cpp)
  list(APPEND HDRS conduitdatacollection.hpp)
@@ -41,14 +41,9 @@ 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
   {
@@ -95,10 +90,6 @@ 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));

@@ -154,8 +145,6 @@ 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);
      }
   }

@@ -187,11 +176,6 @@ 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)
   {
@@ -771,4 +755,4 @@ void PLBound::Print(std::ostream &outp) const
   ubound.Print(outp);
 }

-}
+}
@@ -13,7 +13,6 @@
 #define MFEM_BOUNDS

 #include "../config/config.hpp"
-#include "../general/forall.hpp"
 #include "fespace.hpp"

 namespace mfem
@@ -61,9 +60,7 @@ 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
@@ -116,10 +113,7 @@ 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.
@@ -143,23 +137,9 @@ 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
@@ -203,600 +183,6 @@ 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
@@ -52,9 +52,6 @@ public:
   /// Get the time for time dependent coefficients
   real_t GetTime() { return time; }

-   /// Returns dimension of the vector.
-   int GetVDim() { return 1; }
-
   /** @brief Evaluate the coefficient in the element described by @a T at the
       point @a ip. */
   /** @note When this method is called, the caller must make sure that the
@@ -1453,8 +1453,6 @@ ConduitDataCollection::LoadMeshAndFields(int domain_id,
 std::string
 ConduitDataCollection::ElementTypeToShapeName(Element::Type element_type)
 {
-   // Adapted from SidreDataCollection
-
   // Note -- the mapping from Element::Type to string is based on
   //   enum Element::Type { POINT, SEGMENT, TRIANGLE, QUADRILATERAL,
   //                        TETRAHEDRON, HEXAHEDRON };
@@ -34,10 +34,10 @@ namespace mfem
    - HDF5 library, https://support.hdfgroup.org/HDF5

    @note The ConduitDataCollection only wraps the mfem objects to save them and
-    creates them on load, Conduit does not own any of the data. The
-    SidreDataCollection provides more features, for example the
-    SidreDataCollection allocates and will own the data backing the mfem objects
-    in the data collection.
+    creates them on load, Conduit does not own any of the data.
+    The MFEMSidreDataCollection in the Axom package (https://github.com/LLNL/axom)
+    derives from mfem::DataCollection and provides more features, for example
+    it allocates and will own the data backing the mfem objects in the data collection.

    This class also provides public static methods that convert between MFEM
    Meshes and GridFunctions and Conduit Mesh Blueprint descriptions.
@@ -492,8 +492,6 @@ void VisItDataCollection::SaveRootFile()
                           to_padded_string(cycle, pad_digits_cycle) +
                           ".mfem_root";
   std::ofstream root_file(root_name);
-   MFEM_VERIFY(root_file.is_open(),
-               "Failed to open ofstream " << root_name);
   root_file << GetVisItRootString();
   if (!root_file)
   {
@@ -979,10 +977,7 @@ void ParaViewDataCollection::Save()
   // Save the local part of the mesh and grid functions fields to the local
   // VTU file. Also save coefficient fields.
   {
-      std::string os_str = vtu_prefix + GenerateVTUFileName("proc", myid);
-      std::ofstream os(os_str);
-      MFEM_VERIFY(os.is_open(),
-                  "Failed to open ofstream " << os_str);
+      std::ofstream os(vtu_prefix + GenerateVTUFileName("proc", myid));
      os.precision(precision);
      SaveDataVTU(os, levels_of_detail);
   }
@@ -994,10 +989,7 @@ void ParaViewDataCollection::Save()
                  "QuadratureFunction output is not supported for "
                  "ParaViewDataCollection on domain boundary!");
      const std::string &field_name = qfield.first;
-      std::string os_str = vtu_prefix + GenerateVTUFileName(field_name, myid);
-      std::ofstream os(os_str);
-      MFEM_VERIFY(os.is_open(),
-                  "Failed to open ofstream " << os_str);
+      std::ofstream os(vtu_prefix + GenerateVTUFileName(field_name, myid));
      qfield.second->SaveVTU(os, pv_data_format, GetCompressionLevel(), field_name);
   }

@@ -1008,10 +1000,7 @@ void ParaViewDataCollection::Save()
   {
      // Create the main PVTU file
      {
-         std::string os_str = vtu_prefix + GeneratePVTUFileName("data");
-         std::ofstream pvtu_out(os_str);
-         MFEM_VERIFY(pvtu_out.is_open(),
-                     "Failed to open ofstream " << os_str);
+         std::ofstream pvtu_out(vtu_prefix + GeneratePVTUFileName("data"));
         WritePVTUHeader(pvtu_out);

         // Grid function fields and coefficient fields
@@ -1066,10 +1055,8 @@ void ParaViewDataCollection::Save()
         const std::string &q_field_name = q_field.first;
         std::string q_fname = GeneratePVTUPath() + "/"
                               + GeneratePVTUFileName(q_field_name);
-         std::string os_str = col_path + "/" + q_fname;
-         std::ofstream pvtu_out(os_str);
-         MFEM_VERIFY(pvtu_out.is_open(),
-                     "Failed to open ofstream " << os_str);
+
+         std::ofstream pvtu_out(col_path + "/" + q_fname);
         WritePVTUHeader(pvtu_out);
         int vec_dim = q_field.second->GetVDim();
         pvtu_out << "<PPointData>\n";
@@ -57,7 +57,7 @@ void DGMassApply(const int e,
      }
      else if (DIM == 3)
      {
-         SmemPAMassApply3D_Element<TD1D,TQ1D,NBZ,ACCUM>(e, NE, B, pa_data, x, y);
+         SmemPAMassApply3D_Element<TD1D,TQ1D,ACCUM>(e, NE, B, pa_data, x, y);
      }
      else
      {
@@ -320,8 +320,8 @@ public:
    error estimation procedure where the flux averaging is replaced by a global
    L2 projection (requiring a mass matrix solve).

-    The required BilinearFormIntegrator must implement the method
-    ComputeElementFlux().
+    The required BilinearFormIntegrator must implement the methods
+    ComputeElementFlux() and ComputeFluxEnergy().

    Implemented for the parallel case only.
 */
@@ -357,8 +357,8 @@ protected:

 public:
   /** @brief Construct a new L2ZienkiewiczZhuEstimator object.
-       @param integ    This BilinearFormIntegrator must implement the method
-                       ComputeElementFlux().
+       @param integ    This BilinearFormIntegrator must implement the methods
+                       ComputeElementFlux() and ComputeFluxEnergy().
       @param sol      The solution field whose error is to be estimated.
       @param flux_fes The L2ZienkiewiczZhuEstimator assumes ownership of this
                       FiniteElementSpace and will call its Update() method when
@@ -382,8 +382,8 @@ public:
   { }

   /** @brief Construct a new L2ZienkiewiczZhuEstimator object.
-       @param integ    This BilinearFormIntegrator must implement the method
-                       ComputeElementFlux().
+       @param integ    This BilinearFormIntegrator must implement the methods
+                       ComputeElementFlux() and ComputeFluxEnergy().
       @param sol      The solution field whose error is to be estimated.
       @param flux_fes The L2ZienkiewiczZhuEstimator does NOT assume ownership
                       of this FiniteElementSpace; will call its Update() method
@@ -62,10 +62,6 @@
 #include "pnonlinearform.hpp"
 #endif

-#ifdef MFEM_USE_SIDRE
-#include "sidredatacollection.hpp"
-#endif
-
 #ifdef MFEM_USE_CONDUIT
 #include "conduitdatacollection.hpp"
 #endif
@@ -3030,14 +3030,10 @@ void GridFunction::ProjectCoefficient(Coefficient *coeff[])
   }
 }

-void GridFunction::ProjectDiscCoefficient(
-   std::variant<Coefficient*, VectorCoefficient*> coeff, Array<int> &dof_attr)
+void GridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff,
+                                          Array<int> &dof_attr)
 {
-   std::visit([&](auto* c)
-   {
-      MFEM_VERIFY(VectorDim() == c->GetVDim(), "coeff vdim != VectorDim()");
-   }, coeff);
-
+   MFEM_VERIFY(VectorDim() == coeff.GetVDim(), "coeff vdim != VectorDim()");
   Array<int> vdofs;
   Vector vals;

@@ -3051,10 +3047,7 @@ void GridFunction::ProjectDiscCoefficient(
   {
      fes->GetElementVDofs(i, vdofs);
      vals.SetSize(vdofs.Size());
-      std::visit([&](auto* c)
-      {
-         fes->GetFE(i)->Project(*c, *fes->GetElementTransformation(i), vals);
-      }, coeff);
+      fes->GetFE(i)->Project(coeff, *fes->GetElementTransformation(i), vals);

      // the values in shared dofs are determined from the element with maximal
      // attribute
@@ -3070,6 +3063,13 @@ void GridFunction::ProjectDiscCoefficient(
   }
 }

+void GridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff)
+{
+   MFEM_VERIFY(VectorDim() == coeff.GetVDim(), "coeff vdim != VectorDim()");
+   Array<int> dof_attr;
+   ProjectDiscCoefficient(coeff, dof_attr);
+}
+
 void GridFunction::ProjectDiscCoefficient(Coefficient &coeff, AvgType type)
 {
   // Harmonic  (x1 ... xn) = [ (1/x1 + ... + 1/xn) / n ]^-1.
@@ -5252,30 +5252,6 @@ 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));
@@ -23,7 +23,6 @@
 #include <limits>
 #include <ostream>
 #include <string>
-#include <variant>

 namespace mfem
 {
@@ -80,18 +79,10 @@ protected:
                        bool wcoef,
                        int subdomain);

-   /** @brief Project a discontinuous (vector) coefficient as a grid function on
-       a continuous finite element space. Return in dof_attr the maximal
-       attribute of the elements containing each degree of freedom. */
-   virtual void ProjectDiscCoefficient(
-      std::variant<Coefficient*, VectorCoefficient*> coeff, Array<int> &dof_attr);
-
-   /** @brief Project a discontinuous (vector) coefficient as a grid function on
-       a continuous finite element space. The values in shared dofs are
-       determined from the element with maximal attribute. */
-   virtual void ProjectDiscCoefficient(
-      std::variant<Coefficient*, VectorCoefficient*> coeff)
-   { Array<int> dof_attr; ProjectDiscCoefficient(coeff, dof_attr); };
+   /** Project a discontinuous vector coefficient in a continuous space and
+       return in dof_attr the maximal attribute of the elements containing each
+       degree of freedom. */
+   void ProjectDiscCoefficient(VectorCoefficient &coeff, Array<int> &dof_attr);

   /** Helper function for ProjectCoefficientElementL2 */
   void ProjectCoefficientElementL2_(Coefficient &coeff, Vector &sol, Vector &Va);
@@ -524,17 +515,10 @@ public:
       but using an array of scalar coefficients for each component. */
   void ProjectCoefficient(Coefficient *coeff[]);

-   /** @brief Project a discontinuous coefficient as a grid function on
-       a continuous finite element space. The values in shared dofs are
-       determined from the element with maximal attribute. */
-   virtual void ProjectDiscCoefficient(Coefficient &coeff)
-   { ProjectDiscCoefficient(&coeff); }
-
   /** @brief Project a discontinuous vector coefficient as a grid function on
       a continuous finite element space. The values in shared dofs are
       determined from the element with maximal attribute. */
-   virtual void ProjectDiscCoefficient(VectorCoefficient &coeff)
-   { ProjectDiscCoefficient(&coeff); }
+   virtual void ProjectDiscCoefficient(VectorCoefficient &coeff);

   enum AvgType {ARITHMETIC, HARMONIC};
   /** @brief Projects a discontinuous coefficient so that the values in shared
@@ -106,7 +106,9 @@ FindPointsGSLIB::FindPointsGSLIB()
   gsl_comm = new gslib::comm;
   cr       = new gslib::crystal;
 #ifdef MFEM_USE_MPI
-   if (!Mpi::IsInitialized()) { Mpi::Init(); }
+   int initialized = 0;
+   MPI_Initialized(&initialized);
+   if (!initialized) { MPI_Init(NULL, NULL); }
   MPI_Comm comm = MPI_COMM_WORLD;
   comm_init(gsl_comm, comm);
 #else
@@ -488,7 +490,7 @@ void FindPointsGSLIB::FindPointsOnDevice(const Vector &point_pos,
   }
   DEV.find_device = true;

-   const unsigned int id = gsl_comm->id, np = gsl_comm->np;
+   const int id = gsl_comm->id, np = gsl_comm->np;

   gsl_mfem_ref.SetSize(points_cnt * dim);
   gsl_mfem_elem.SetSize(points_cnt);
@@ -650,7 +652,7 @@ void FindPointsGSLIB::FindPointsOnDevice(const Vector &point_pos,
         {
            const int pp = hash_offset[i];
            /* don't send back to where it just came from */
-            if (static_cast<unsigned>(pp) == p->proc)
+            if (pp == p->proc)
            {
               continue;
            }
@@ -1066,7 +1068,7 @@ void FindPointsGSLIB::InterpolateOnDevice(const Vector &field_in_evec,
         sarray_transfer(struct evalOutPt_t, &outpt, proc, 1, cr);

         opt = (evalOutPt_t *)outpt.ptr;
-         for (size_t index = 0; index < outpt.n; index++)
+         for (int index = 0; index < outpt.n; index++)
         {
            int idx = ordering == Ordering::byNODES ?
                      opt->index + i*points_cnt :
@@ -1411,7 +1413,7 @@ void FindPointsGSLIB::SetupSplitMeshesAndIntegrationRules(const int order)
 {
   MFEM_VERIFY(mesh, "Setup FindPointsGSLIB with mesh first.");
   const int dof1D = order+1;
-   dim = mesh->Dimension();
+   const int dim = mesh->Dimension();

   SetupSplitMeshes();
   if (dim == 2)
@@ -2252,8 +2254,7 @@ void FindPointsGSLIB::DistributeInterpolatedValues(const Vector &int_vals,
      sarray_transfer(struct out_pt, outpt, proc, 1, cr);

      // Store received data
-      MFEM_VERIFY(outpt->n == static_cast<size_t>(points_cnt),
-                  "Incompatible size. Number of points "
+      MFEM_VERIFY(outpt->n == points_cnt, "Incompatible size. Number of points "
                  "received does not match the number of points originally "
                  "found using FindPoints.");

@@ -2622,7 +2623,9 @@ GSOPGSLIB::GSOPGSLIB(Array<long long> &ids)
   gsl_comm = new gslib::comm;
   cr       = new gslib::crystal;
 #ifdef MFEM_USE_MPI
-   if (!Mpi::IsInitialized()) { Mpi::Init(); }
+   int initialized;
+   MPI_Initialized(&initialized);
+   if (!initialized) { MPI_Init(NULL, NULL); }
   MPI_Comm comm = MPI_COMM_WORLD;
   comm_init(gsl_comm, comm);
 #else
@@ -202,19 +202,13 @@ protected:
                            const int dof1dsol, const int ordering);

 public:
-   /// Serial constructor
   FindPointsGSLIB();
-
-   /// Serial constructor + setup with given Mesh (see \ref Setup)
   FindPointsGSLIB(Mesh &mesh_in, const double bb_t = 0.1,
                   const double newt_tol = 1.0e-12,
                   const int npt_max = 256);

 #ifdef MFEM_USE_MPI
-   /// Constructor for ParMesh
   FindPointsGSLIB(MPI_Comm comm_);
-
-   /// Constructor + setup with given ParMesh (see \ref Setup)
   FindPointsGSLIB(ParMesh &mesh_in, const double bb_t = 0.1,
                   const double newt_tol = 1.0e-12,
                   const int npt_max = 256);
@@ -254,7 +254,7 @@ get_edge(const double *elx[2], const double *wtend, int ei,
      edge.dxdn[d] = workspace + (2 + d) * pN; //dxdn and dydn at DOFs along edge
   }

-   if (static_cast<unsigned>(side_init) != (1u << ei))
+   if (side_init != (1u << ei))
   {
 #define ELX(d, j, k) elx[d][j + k * pN] // assumes lexicographic ordering
      for (int d = 0; d < 2; ++d)
@@ -294,7 +294,7 @@ get_face(const double *elx[3], const double *wtend, int fi, double *workspace,
      face.dxdn[d] = workspace+(3+d)*p_Nfr;
   }

-   if (static_cast<unsigned>(side_init) != (1u << fi))
+   if (side_init != (1u << fi))
   {
      const int e_stride[3] = {1, pN, pN*pN};
 #define ELX(d, j, k, l) elx[d][j*e_stride[d1]+k*e_stride[d2]+l*e_stride[dn]]
@@ -342,7 +342,7 @@ get_edge(const double *elx[3], const double *wtend, int ei, double *workspace,

   if (jidx >= 3*pN) { return edge; }

-   if (static_cast<unsigned>(side_init) != (64u << ei))
+   if (side_init != (64u << ei))
   {
      const int e_stride[3] = {1, pN, pN*pN};
 #define ELX(d, j, k, l) elx[d][j*e_stride[de]+k*e_stride[dn1]+l*e_stride[dn2]]
@@ -181,12 +181,6 @@ 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
@@ -810,23 +804,19 @@ void PAMassApply3D_Element(const int e,
   }
 }

-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)
+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)
 {
-   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 D1D = T_D1D ? T_D1D : d1d;
+   constexpr 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;
@@ -839,37 +829,33 @@ SmemPAMassApply3D_Element(const int e, const int NE, const real_t *b_,
   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[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_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_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) { B[dx][dy] = b(dx, dy); }
-   }
-   if (tidz == 0)
-   {
-      MFEM_FOREACH_THREAD(dy, y, D1D)
+      MFEM_FOREACH_THREAD(dx,x,Q1D)
      {
-         MFEM_FOREACH_THREAD(dx, x, Q1D) { B[dx][dy] = b(dx, dy); }
+         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)
@@ -894,9 +880,9 @@ SmemPAMassApply3D_Element(const int e, const int NE, const real_t *b_,
      }
   }
   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)
@@ -921,9 +907,9 @@ SmemPAMassApply3D_Element(const int e, const int NE, const real_t *b_,
      }
   }
   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)
@@ -943,22 +929,22 @@ SmemPAMassApply3D_Element(const int e, const int NE, const real_t *b_,
         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;
-   if (tidz == 0)
+   MFEM_FOREACH_THREAD(di,y,D1D)
   {
-      MFEM_FOREACH_THREAD(di, y, D1D)
+      MFEM_FOREACH_THREAD(q,x,Q1D)
      {
-         MFEM_FOREACH_THREAD(q, x, Q1D) { Bt[di][q] = b(q, di); }
+         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)
@@ -983,9 +969,9 @@ SmemPAMassApply3D_Element(const int e, const int NE, const real_t *b_,
      }
   }
   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)
@@ -1010,9 +996,9 @@ SmemPAMassApply3D_Element(const int e, const int NE, const real_t *b_,
      }
   }
   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)
@@ -1034,11 +1020,11 @@ SmemPAMassApply3D_Element(const int e, const int NE, const real_t *b_,
         {
            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];
            }
         }
      }
@@ -1129,8 +1115,8 @@ inline void PAMassApply3D(const int NE,
   });
 }

-// Shared memory PA Mass Apply 3D kernel
-template<int T_D1D = 0, int T_Q1D = 0, int TBATCH=1>
+// Shared memory PA Mass Apply 2D kernel
+template<int T_D1D = 0, int T_Q1D = 0>
 inline void SmemPAMassApply3D(const int NE,
                              const Array<real_t> &b_,
                              const Array<real_t> &bt_,
@@ -1140,9 +1126,6 @@ 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;
@@ -1154,11 +1137,9 @@ inline void SmemPAMassApply3D(const int NE,
   const auto d = d_.Read();
   const auto x = x_.Read();
   auto y = y_.ReadWrite();
-   mfem::forall_2D_batch<T_Q1D * T_Q1D * TBATCH>(NE, Q1D, Q1D, TBATCH,
-                                                 [=] MFEM_HOST_DEVICE(int e)
+   mfem::forall_2D<T_Q1D*T_Q1D>(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE (int e)
   {
-      internal::SmemPAMassApply3D_Element<T_D1D, T_Q1D, TBATCH>(e, NE, b, d, x,
-                                                                y, d1d, q1d);
+      internal::SmemPAMassApply3D_Element<T_D1D,T_Q1D>(e, NE, b, d, x, y, d1d, q1d);
   });
 }

@@ -1413,16 +1394,7 @@ 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)
-   {
-      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)>;
-      }
-   }
+   else if constexpr (DIM == 3) { return internal::SmemPAMassApply3D<T_D1D, T_Q1D>; }
   MFEM_ABORT("");
 }

@@ -43,52 +43,56 @@ public:
      index = i;
   }

-   void Set3w(const real_t x1, const real_t x2, const real_t x3, const real_t w)
-   { x = x1; y = x2; z = x3; weight = w; }
-   void Set2w(const real_t x1, const real_t x2, const real_t w)
-   { x = x1; y = x2; weight = w; }
-   void Set1w(const real_t x1, const real_t w)
-   { x = x1; weight = w; }
-
-   void Set3w(const real_t *p) { Set3w(p[0], p[1], p[2], p[3]); }
-   void Set2w(const real_t *p) { Set2w(p[0], p[1], p[2]); }
-   void Set1w(const real_t *p) { Set1w(p[0], p[1]); }
-
-   void Set3(const real_t x1, const real_t x2, const real_t x3)
-   { x = x1; y = x2; z = x3; }
-   void Set2(const real_t x1, const real_t x2)
-   { x = x1; y = x2; }
-   void Set1(const real_t x1)
-   { x = x1; }
-
-   void Set3(const real_t *p) { Set3(p[0], p[1], p[2]); }
-   void Set2(const real_t *p) { Set2(p[0], p[1]); }
-   void Set1(const real_t *p) { Set1(p[0]); }
-
-   void Set(const real_t x1, const real_t x2, const real_t x3, const real_t w)
-   { Set3w(x1, x2, x3, w); }
-
   void Set(const real_t *p, const int dim)
   {
      MFEM_ASSERT(1 <= dim && dim <= 3, "invalid dim: " << dim);
-      switch (dim)
+      x = p[0];
+      if (dim > 1)
      {
-         case 3: Set3(p); break;
-         case 2: Set2(p); break;
-         case 1: Set1(p); break;
+         y = p[1];
+         if (dim > 2)
+         {
+            z = p[2];
+         }
      }
   }

   void Get(real_t *p, const int dim) const
   {
      MFEM_ASSERT(1 <= dim && dim <= 3, "invalid dim: " << dim);
-      switch (dim)
+      p[0] = x;
+      if (dim > 1)
      {
-         case 3: p[2] = z;
-         case 2: p[1] = y;
-         case 1: p[0] = x;
+         p[1] = y;
+         if (dim > 2)
+         {
+            p[2] = z;
+         }
      }
   }
+
+   void Set(const real_t x1, const real_t x2, const real_t x3, const real_t w)
+   { x = x1; y = x2; z = x3; weight = w; }
+
+   void Set3w(const real_t *p) { x = p[0]; y = p[1]; z = p[2]; weight = p[3]; }
+
+   void Set3(const real_t x1, const real_t x2, const real_t x3)
+   { x = x1; y = x2; z = x3; }
+
+   void Set3(const real_t *p) { x = p[0]; y = p[1]; z = p[2]; }
+
+   void Set2w(const real_t x1, const real_t x2, const real_t w)
+   { x = x1; y = x2; weight = w; }
+
+   void Set2w(const real_t *p) { x = p[0]; y = p[1]; weight = p[2]; }
+
+   void Set2(const real_t x1, const real_t x2) { x = x1; y = x2; }
+
+   void Set2(const real_t *p) { x = p[0]; y = p[1]; }
+
+   void Set1w(const real_t x1, const real_t w) { x = x1; weight = w; }
+
+   void Set1w(const real_t *p) { x = p[0]; weight = p[1]; }
 };

 /// Class for an integration rule - an Array of IntegrationPoint.
@@ -164,8 +164,8 @@ private:

 public:
   /// Constructs the domain integrator $ (Q, \nabla v) $
-   DomainLFGradIntegrator(VectorCoefficient &QF, const IntegrationRule *ir = NULL)
-      : DeltaLFIntegrator(QF, ir), Q(QF) { }
+   DomainLFGradIntegrator(VectorCoefficient &QF)
+      : DeltaLFIntegrator(QF), Q(QF) { }

   bool SupportsDevice() const override { return true; }

@@ -717,9 +717,9 @@ void ParGridFunction::ProjectCoefficientElementL2(VectorCoefficient &vcoeff)
 }


-void ParGridFunction::ProjectDiscCoefficient(
-   std::variant<Coefficient*, VectorCoefficient*> coeff)
+void ParGridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff)
 {
+   MFEM_VERIFY(VectorDim() == coeff.GetVDim(), "coeff vdim != VectorDim()");
   // local maximal element attribute for each dof
   Array<int> ldof_attr;

@@ -63,12 +63,6 @@ protected:
   void ProjectBdrCoefficient(Coefficient *coeff[], VectorCoefficient *vcoeff,
                              const Array<int> &attr);

-   /** @brief Project a discontinuous (vector) coefficient as a grid function on
-       a continuous finite element space. The values in shared dofs are
-       determined from the element with maximal attribute. */
-   virtual void ProjectDiscCoefficient(
-      std::variant<Coefficient*, VectorCoefficient*> coeff) override;
-
 public:
   ParGridFunction() { pfes = NULL; }

@@ -274,6 +268,11 @@ public:
                           ProjectType type = ProjectType::DEFAULT) override;

   using GridFunction::ProjectDiscCoefficient;
+   /** @brief Project a discontinuous vector coefficient as a grid function on
+       a continuous finite element space. The values in shared dofs are
+       determined from the element with maximal attribute. */
+   void ProjectDiscCoefficient(VectorCoefficient &coeff) override;
+
   void ProjectDiscCoefficient(Coefficient &coeff, AvgType type) override;

   void ProjectDiscCoefficient(VectorCoefficient &vcoeff, AvgType type) override;
@@ -1,539 +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_SIDREDATACOLLECTION
-#define MFEM_SIDREDATACOLLECTION
-
-#include "../config/config.hpp"
-
-#ifdef MFEM_USE_SIDRE
-
-#include "datacollection.hpp"
-
-// Ignore warnings from the axom/sidre header (GCC + Clang versions)
-#ifdef MFEM_HAVE_GCC_PRAGMA_DIAGNOSTIC
-# pragma GCC diagnostic push
-# if defined(__clang__)
-#  pragma GCC diagnostic ignored "-Wextra-semi"
-# else // real GCC?
-#  pragma GCC diagnostic ignored "-Wpedantic"
-# endif
-#endif
-#include <axom/sidre.hpp>
-#ifdef MFEM_HAVE_GCC_PRAGMA_DIAGNOSTIC
-# pragma GCC diagnostic pop
-#endif
-
-namespace mfem
-{
-
-/** @brief Data collection with Sidre routines following the Conduit mesh
-    blueprint specification. */
-/** SidreDataCollection provides an HDF5-based file format for visualization or
-    restart capability.  This functionality is aimed primarily at customers of
-    LLNL's axom project that run problems at extreme scales.
-
-    For more information, see:
-    - Sidre component of LLNL's axom project (to be open-sourced), http://goo.gl/cZyJdn
-    - LLNL conduit/blueprint library, https://github.com/LLNL/conduit
-    - HDF5 library, https://support.hdfgroup.org/HDF5
-
-    The layout created in the Sidre DataStore is: (`"──"` denote groups,
-    `"─•"` denote views, `"─>"` denote links, i.e. shallow-copy view)
-
-        <root>
-         ├── <collection-name>_global          (global group)
-         │    └── blueprint_index
-         │         └── <collection-name>       (bp_index group)
-         │              ├── state
-         │              │    ├─• cycle
-         │              │    ├─• time
-         │              │    └─• number_of_domains = <mesh-mpi-comm-size>
-         │              ├── coordsets
-         │              │    └── coords
-         │              │         ├─• path = "<bp-path>/coordsets/coords"
-         │              │         ├─• type ─> <bp-grp>/coordsets/coords/type = "explicit"
-         │              │         └─• coord_system = "x"|"xy"|"xyz"
-         │              ├── topologies
-         │              │    ├── mesh
-         │              │    │    ├─• path = "<bp-path>/topologies/mesh"
-         │              │    │    ├─• type              ─> <bp-grp>/topologies/mesh/type = "unstructured"
-         │              │    │    ├─• coordset          ─> <bp-grp>/topologies/mesh/coordset = "coords"
-         │              │    │    ├─• grid_function     ─> <bp-grp>/topologies/mesh/grid_function = "<nodes-field-name>"
-         │              │    │    └─• boundary_topology ─> <bp-grp>/topologies/mesh/boundary_topology = "boundary"
-         │              │    └── boundary
-         │              │         ├─• path = "<bp-path>/topologies/mesh"
-         │              │         ├─• type     ─> <bp-grp>/topologies/boundary/type = "unstructured"
-         │              │         └─• coordset ─> <bp-grp>/topologies/boundary/coordset = "coords"
-         │              └── fields
-         │                   ├── mesh_material_attribute
-         │                   │    ├─• path = "<bp-path>/fields/mesh_material_attribute"
-         │                   │    ├─• association ─> <bp-grp>/fields/mesh_material_attribute/association = "element"
-         │                   │    ├─• topology    ─> <bp-grp>/fields/mesh_material_attribute/topology = "mesh"
-         │                   │    └─• number_of_components = 1
-         │                   ├── boundary_material_attribute
-         │                   │    ├─• path = "<bp-path>/fields/boundary_material_attribute"
-         │                   │    ├─• association ─> <bp-grp>/fields/boundary_material_attribute/association = "element"
-         │                   │    ├─• topology    ─> <bp-grp>/fields/boundary_material_attribute/topology = "boundary"
-         │                   │    └─• number_of_components = 1
-         │                   ├── grid-function-1
-         │                   │    ├─• path = "<bp-path>/fields/grid-function-1"
-         │                   │    ├─• basis    ─> <bp-grp>/fields/grid-function-1/basis = "<fe-coll-name>"
-         │                   │    ├─• topology ─> <bp-grp>/fields/grid-function-1/topology = "mesh"
-         │                   │    └─• number_of_components = gf1->VectorDim()
-         │                   ├── grid-function-2
-         │                   │    ├─• path = "<bp-path>/fields/grid-function-2"
-         │                   │    ├─• basis    ─> <bp-grp>/fields/grid-function-2/basis = "<fe-coll-name>"
-         │                   │    ├─• topology ─> <bp-grp>/fields/grid-function-2/topology = "mesh"
-         │                   │    └─• number_of_components = gf2->VectorDim()
-         │                   ├── ...
-         │                  ...
-         └── <collection-name>                 (domain group)
-              ├── blueprint                    (blueprint group)
-              │    ├── state
-              │    │    ├─• cycle
-              │    │    ├─• time
-              │    │    ├─• domain = <mesh-mpi-rank>
-              │    │    └─• time_step
-              │    ├── coordsets
-              │    │    └── coords
-              │    │         ├─• type = "explicit"
-              │    │         └── values
-              │    │              ├─• x = view in <vertex-coords-buffer>/<ext-double-data>
-              │    │              ├─• y = view in <vertex-coords-buffer>/<ext-double-data>
-              │    │              └─• z = view in <vertex-coords-buffer>/<ext-double-data>
-              │    ├── topologies
-              │    │    ├── mesh
-              │    │    │    ├─• type = "unstructured"
-              │    │    │    ├── elements
-              │    │    │    │    ├─• shape = "points"|"lines"|...
-              │    │    │    │    └─• connectivity = <vert-idx-array>
-              │    │    │    ├─• coordset = "coords"
-              │    │    │    ├─• grid_function = "<nodes-field-name>"
-              │    │    │    └─• boundary_topology = "boundary"
-              │    │    └── boundary
-              │    │         ├─• type = "unstructured"
-              │    │         ├── elements
-              │    │         │    ├─• shape = "points"|"lines"|...
-              │    │         │    └─• connectivity = <vert-idx-array>
-              │    │         └─• coordset = "coords"
-              │    └── fields
-              │         ├── mesh_material_attribute
-              │         │    ├─• association = "element"
-              │         │    ├─• topology = "mesh"
-              │         │    └─• values = <attr-array>
-              │         ├── boundary_material_attribute
-              │         │    ├─• association = "element"
-              │         │    ├─• topology = "boundary"
-              │         │    └─• values = <attr-array>
-              │         ├── grid-function-1   (name can include path)
-              │         │    ├─• basis = "<fe-coll-name>"
-              │         │    ├─• topology = "mesh"
-              │         │    └─• values = <ext-double-array>/<named-buffer> (vdim == 1)
-              │         ├── grid-function-2   (name can include path)
-              │         │    ├─• basis = "<fe-coll-name>"
-              │         │    ├─• topology = "mesh"
-              │         │    └── values   (vdim > 1)
-              │         │         ├─• x0 = view into <ext-double-array>/<named-buffer>
-              │         │         ├─• x1 = view into <ext-double-array>/<named-buffer>
-              │         │         └─• x2 = view into <ext-double-array>/<named-buffer>
-              │         ├── ...
-              │        ...
-              └── named_buffers                (named_buffers group)
-                   ├─• vertex_coords = <double-array>
-                   ├─• grid-function-1 = <double-array>
-                   ├─• grid-function-2 = <double-array>
-                  ...
-
-    @note blueprint_index is used both in serial and in parallel. In parallel,
-    only rank 0 will add entries to the blueprint index.
-
-    @note QuadratureFunction%s (q-fields) are not supported.
-
-    @note SidreDataCollection does not manage the FiniteElementSpace%s and
-    FiniteElementCollection%s associated with registered GridFunction%s.
-    Therefore, field registration is left to the user of SidreDataCollection and
-    there are no methods that automatically register GridFunction%s using just
-    the content of the Sidre DataStore. Such capabilities can be implemented in
-    a derived class, adding any desired object management routines.
-
-    @warning This class is still _experimental_, meaning that in future
-    releases, it may not be backward compatible, and the output files generated
-    by the current version may become unreadable.
-*/
-class SidreDataCollection : public DataCollection
-{
-public:
-   typedef NamedFieldsMap< Array<int> > AttributeFieldMap;
-   AttributeFieldMap attr_map;
-
-public:
-
-   /// Constructor that allocates and initializes a Sidre DataStore.
-   /**
-       @param[in] collection_name  Name of the collection used as a file name
-                                   when saving
-       @param[in] the_mesh         Mesh shared by all grid functions in the
-                                   collection (can be NULL)
-       @param[in] owns_mesh_data   Does the SidreDC own the mesh vertices?
-
-       With this constructor, the SidreDataCollection owns the allocated Sidre
-       DataStore.
-    */
-   SidreDataCollection(const std::string& collection_name,
-                       Mesh *the_mesh = NULL,
-                       bool owns_mesh_data = false);
-
-   /// Constructor that links to an external Sidre DataStore.
-   /** Specifically, the global and domain groups can be at arbitrary paths.
-
-       @param[in] collection_name  Name of the collection used as a file name
-                                   when saving
-       @param[in] bp_index_grp     Pointer to the blueprint index group in the
-                                   datastore, see the above schematic
-       @param[in] domain_grp       Pointer to the domain group in the datastore,
-                                   see the above schematic
-       @param[in] owns_mesh_data   Does the SidreDC own the mesh vertices?
-
-       With this constructor, the SidreDataCollection does not own the Sidre
-       DataStore.
-       @note No mesh or fields are read from the given Groups. The mesh has
-       to be set with SetMesh() and fields registered with RegisterField().
-    */
-   SidreDataCollection(const std::string& collection_name,
-                       axom::sidre::Group * bp_index_grp,
-                       axom::sidre::Group * domain_grp,
-                       bool owns_mesh_data = false);
-
-#ifdef MFEM_USE_MPI
-   /// Associate an MPI communicator with the collection.
-   /** If no mesh was associated with the collection, this method should be
-       called before using any of the Load() methods to read parallel data. */
-   void SetComm(MPI_Comm comm);
-#endif
-
-   /// Register a GridFunction in the Sidre DataStore.
-   /** This method is a shortcut for the call
-       `RegisterField(field_name, gf, field_name, 0)`.
-    */
-   virtual void RegisterField(const std::string &field_name, GridFunction *gf)
-   {
-      RegisterField(field_name, gf, field_name, 0);
-   }
-
-   /// Register a GridFunction in the Sidre DataStore.
-   /** The registration procedure is as follows:
-       - if (@a gf's data is NULL), allocate named buffer with the name
-         @a buffer_name with size _offset + gf->FESpace()->GetVSize()_ and use
-         its data (plus the given @a offset) to set @a gf's data;
-       - else, if (DataStore has a named buffer @a buffer_name), replace @a gf's
-         data array with that named buffer plus the given @a offset;
-       - else, use @a gf's data as external data associated with @a field_name
-         in the DataStore;
-       - register @a field_name in #field_map.
-
-       Both the @a field_name and @a buffer_name can contain a path prefix.
-       @note If @a field_name or @a buffer_name is empty, the method does
-       nothing.
-       @note If the GridFunction pointer @a gf or it's FiniteElementSpace
-       pointer are NULL, the method does nothing.
-    */
-   void RegisterField(const std::string &field_name, GridFunction *gf,
-                      const std::string &buffer_name,
-                      axom::sidre::IndexType offset);
-
-   /// Registers an attribute field in the Sidre DataStore
-   /** The registration process is similar to that of RegisterField()
-       The attribute field is associated with the elements of the mesh
-       when @a is_bdry is false, and with the boundary elements, when
-       @a is_bdry is true.
-       @sa RegisterField()  */
-   void RegisterAttributeField(const std::string& name, bool is_bdry);
-   void DeregisterAttributeField(const std::string& name);
-
-   /** Returns a pointer to the attribute field associated with
-       @a field_name, or NULL when there is no associated field */
-   Array<int>* GetAttributeField(const std::string& field_name) const
-   { return attr_map.Get(field_name); }
-
-   /** Checks if there is an attribute field associated with @a field_name */
-   bool HasAttributeField(const std::string& field_name) const
-   { return attr_map.Has(field_name); }
-
-   /** Checks if any rank in the mesh has boundary elements */
-   bool HasBoundaryMesh() const;
-
-   /// Set the name of the mesh nodes field.
-   /** This name will be used by SetMesh() to register the mesh nodes, if not
-       already registered. Also, this method should be called if the mesh nodes
-       GridFunction was or will be registered directly by the user. The default
-       value for the name is "mesh_nodes". */
-   void SetMeshNodesName(const std::string &nodes_name)
-   {
-      if (!nodes_name.empty()) { m_meshNodesGFName = nodes_name; }
-   }
-
-   /// De-register @a field_name from the SidreDataCollection.
-   /** The field is removed from the #field_map and the DataStore, including
-       deleting it from the named_buffers group, if allocated. */
-   virtual void DeregisterField(const std::string& field_name);
-
-   /// Delete all owned data.
-   virtual ~SidreDataCollection();
-
-   /// Set/change the mesh associated with the collection
-   /** Uses the field name "mesh_nodes" or the value set by SetMeshNodesName()
-       to register the mesh nodes GridFunction, if the mesh uses nodes. */
-   virtual void SetMesh(Mesh *new_mesh);
-
-#ifdef MFEM_USE_MPI
-   /// Set/change the mesh associated with the collection
-   /** Uses the field name "mesh_nodes" or the value set by SetMeshNodesName()
-       to register the mesh nodes GridFunction, if the mesh uses nodes. */
-   virtual void SetMesh(MPI_Comm comm, Mesh *new_mesh);
-#endif
-
-   /// Reset the domain and global datastore group pointers.
-   /** These are set in the constructor, but if a host code changes the
-       datastore contents ( such as wiping out the datastore and loading in new
-       contents from a file, i.e. a restart ) these pointers will need to be
-       reset to valid groups in the datastore.
-       @sa Load(const std::string &path, const std::string &protocol).
-    */
-   void SetGroupPointers(axom::sidre::Group * global_grp,
-                         axom::sidre::Group * domain_grp);
-
-   axom::sidre::Group * GetBPGroup() { return m_bp_grp; }
-   axom::sidre::Group * GetBPIndexGroup() { return m_bp_index_grp; }
-
-   /// Prepare the DataStore for writing
-   virtual void PrepareToSave();
-
-   /// Save the collection to file.
-   /** This method calls `Save(collection_name, "sidre_hdf5")`. */
-   virtual void Save();
-
-   /// Save the collection to @a filename.
-   /** The collection path prefix is prepended to the @a filename and the
-       current cycle is appended, if cycle >= 0. */
-   void Save(const std::string& filename, const std::string& protocol);
-
-   /// Load the Sidre DataStore from file.
-   /** No mesh or fields are read from the loaded DataStore.
-
-       If the data collection created the datastore, it knows the layout of
-       where the domain and global groups are, and can restore them after the
-       Load().
-
-       If, however, the data collection does not own the datastore (e.g. it did
-       not create the datastore), the host code must reset these pointers after
-       the load operation, using SetGroupPointers(), and also reset the state
-       variables, using UpdateStateFromDS().
-    */
-   void Load(const std::string& path, const std::string& protocol);
-
-   /// Load SidreDataCollection from file.
-   /** The used file path is based on the current prefix path, collection name,
-       and the given @a cycle_. The protocol is "sidre_hdf5".
-       @sa Load(const std::string &path, const std::string &protocol).
-    */
-   virtual void Load(int cycle_ = 0)
-   {
-      SetCycle(cycle_);
-      Load(get_file_path(name), "sidre_hdf5");
-   }
-
-   /// Load external data after registering externally owned fields.
-   void LoadExternalData(const std::string& path);
-
-   /** @brief Updates the DataCollection's cycle, time, and time-step variables
-       with the values from the data store. */
-   void UpdateStateFromDS();
-
-   /** @brief Updates the data store's cycle, time, and time-step variables with
-       the values from the SidreDataCollection. */
-   void UpdateStateToDS();
-
-   /** @name Methods for named buffer access and manipulation. */
-   ///@{
-
-   /** @brief Get a pointer to the sidre::View holding the named buffer for
-       @a buffer_name. */
-   /** If such named buffer is not allocated, the method returns NULL.
-       @note To access the underlying pointer, use View::getData().
-       @note To query the size of the buffer, use View::getNumElements().
-    */
-   axom::sidre::View *
-   GetNamedBuffer(const std::string& buffer_name) const
-   {
-      return named_buffers_grp()->hasView(buffer_name)
-             ? named_buffers_grp()->getView(buffer_name)
-             : NULL;
-   }
-
-   /// Return newly allocated or existing named buffer for @a buffer_name.
-   /** The buffer is stored in the named_buffers group. If the currently
-       allocated buffer size is smaller than @a sz, then the buffer is
-       reallocated with size @a sz, destroying its contents.
-       @note To access the underlying pointer, use View::getData().
-    */
-   axom::sidre::View *
-   AllocNamedBuffer(const std::string& buffer_name,
-                    axom::sidre::IndexType sz,
-                    axom::sidre::TypeID type =
-                       axom::sidre::DOUBLE_ID);
-
-   /// Deallocate the named buffer @a buffer_name.
-   void FreeNamedBuffer(const std::string& buffer_name)
-   { named_buffers_grp()->destroyViewAndData(buffer_name); }
-
-   ///@}
-
-private:
-   // Used if the Sidre data collection is providing the datastore itself.
-   const bool m_owns_datastore;
-
-   // TODO - Need to evaluate if this bool member can be combined with own_data
-   // in parent data collection class. m_owns_mesh_data indicates whether the
-   // Sidre dc owns the mesh element data and node positions gf. The DC base
-   // class own_data indicates if the dc owns the mesh object pointer itself and
-   // GF objects. Can we use one flag and just have DC own all objects vs none?
-   const bool m_owns_mesh_data;
-
-   // Name to be used for registering the mesh nodes in the SidreDataCollection.
-   // This name is used by SetMesh() and can be overwritten by the method
-   // SetMeshNodesName().
-   // Default value: "mesh_nodes".
-   std::string m_meshNodesGFName;
-
-   // If the data collection owns the datastore, it will store a pointer to it.
-   // Otherwise, this pointer is NULL.
-   axom::sidre::DataStore * m_datastore_ptr;
-
-protected:
-   axom::sidre::Group *named_buffers_grp() const;
-
-   axom::sidre::View *
-   alloc_view(axom::sidre::Group *grp,
-              const std::string &view_name);
-
-   axom::sidre::View *
-   alloc_view(axom::sidre::Group *grp,
-              const std::string &view_name,
-              const axom::sidre::DataType &dtype);
-
-   axom::sidre::Group *
-   alloc_group(axom::sidre::Group *grp,
-               const std::string &group_name);
-
-   // return the filename based on prefix_path, collection name and cycle.
-   std::string get_file_path(const std::string &filename) const;
-
-private:
-   // If the data collection does not own the datastore, it will need pointers
-   // to the blueprint and blueprint index group to use.
-   axom::sidre::Group * m_bp_grp;
-   axom::sidre::Group * m_bp_index_grp;
-
-   // This is stored for convenience.
-   axom::sidre::Group * m_named_bufs_grp;
-
-   // Private helper functions
-
-   void RegisterFieldInBPIndex(const std::string& field_name,
-                               GridFunction *gf);
-   void DeregisterFieldInBPIndex(const std::string & field_name);
-
-   void RegisterAttributeFieldInBPIndex(const std::string& attr_name);
-   void DeregisterAttributeFieldInBPIndex(const std::string& attr_name);
-
-   /** @brief Return a string with the conduit blueprint name for the given
-       Element::Type. */
-   std::string getElementName( Element::Type elementEnum );
-
-   /**
-    * \brief A private helper function to set up the views associated with the
-       data of a scalar valued grid function in the blueprint style.
-    * \pre gf is not null
-    * \note This function is expected to be called by RegisterField()
-    * \note Handles cases where hierarchy is already set up,
-    *      where the data was allocated by this data collection
-    *      and where the grid function data is external to Sidre
-    */
-   void addScalarBasedGridFunction(const std::string& field_name,
-                                   GridFunction* gf,
-                                   const std::string &buffer_name,
-                                   axom::sidre::IndexType offset);
-
-   /**
-    * \brief A private helper function to set up the views associated with the
-       data of a vector valued grid function in the blueprint style.
-    * \pre gf is not null
-    * \note This function is expected to be called by RegisterField()
-    * \note Handles cases where hierarchy is already set up,
-    *      where the data was allocated by this data collection
-    *      and where the grid function data is external to Sidre
-    */
-   void addVectorBasedGridFunction(const std::string& field_name,
-                                   GridFunction* gf,
-                                   const std::string &buffer_name,
-                                   axom::sidre::IndexType offset);
-
-   /** @brief A private helper function to set up the Views associated with
-       attribute field named @a field_name */
-   void addIntegerAttributeField(const std::string& field_name, bool is_bdry);
-
-   /// Sets up the four main mesh blueprint groups.
-   /**
-    * \param hasBP Indicates whether the blueprint has already been set up.
-    */
-   void createMeshBlueprintStubs(bool hasBP);
-
-   /// Sets up the mesh blueprint 'state' group.
-   /**
-    * \param hasBP Indicates whether the blueprint has already been set up.
-    */
-   void createMeshBlueprintState(bool hasBP);
-
-   /// Sets up the mesh blueprint 'coordsets' group.
-   /**
-    * \param hasBP Indicates whether the blueprint has already been set up.
-    */
-   void createMeshBlueprintCoordset(bool hasBP);
-
-   /// Sets up the mesh blueprint 'topologies' group.
-   /**
-    * This method is called from SetMesh().
-    * \param hasBP Indicates whether the blueprint has already been set up.
-    * \param mesh_name The name of the topology.
-    * \note Valid values for @a mesh_name are "mesh" and "boundary" and the
-            former has to be created with this method before the latter.
-    */
-   void createMeshBlueprintTopologies(bool hasBP, const std::string& mesh_name);
-
-#ifdef MFEM_USE_MPI
-   /// Sets up the mesh blueprint 'adjacencies' group.
-   /**
-    * \param hasBP Indicates whether the blueprint has already been set up.
-    * \note Only valid when using parallel meshes
-    */
-   void createMeshBlueprintAdjacencies(bool hasBP);
-#endif
-
-   /// Verifies that the contents of the mesh blueprint data is valid.
-   void verifyMeshBlueprint();
-};
-
-} // end namespace mfem
-
-#endif
-
-#endif
@@ -14,7 +14,6 @@

 #include "../config/config.hpp"
 #include "array.hpp"
-#include "text.hpp"

 #include <iostream>
 #include <map>
@@ -248,8 +247,7 @@ inline void ArraysByName<T>::Print(std::ostream &os, int width) const
   os << data.size() << '\n';
   for (auto const &it : data)
   {
-      // Note: The method Load() can read any string formatted with std::quoted.
-      os << std::quoted(it.first) << '\n' << it.second.Size() << '\n';
+      os << '"' << it.first << '"' << '\n' << it.second.Size() << '\n';
      it.second.Print(os, width > 0 ? width : it.second.Size());
   }
 }
@@ -260,36 +258,40 @@ void ArraysByName<T>::Load(std::istream &in)
   int NumArrays;
   in >> NumArrays;

-   for (int i = 0; i < NumArrays; i++)
+   std::string ArrayLine, ArrayName;
+   for (int i=0; i < NumArrays; i++)
   {
      in >> std::ws;
-      // Read the name:
-      // - If the stream 'in' starts with " then parse it with the function
-      //   parse_quoted_string() from text.hpp. In this case, the name can be
-      //   empty. Note: this case allows for reading any string formatted using
-      //   std::quoted, e.g. as in the method Print().
-      // - If the name does not start with " then the name ends with the first
-      //   white space character (and the white space character is not included
-      //   in the name). Since white space characters are skipped before reading
-      //   the name, there will be at least one non-white-space character in the
-      //   name in this case.
-      std::string ArrayName;
-      if (in.peek() == '"')
+      getline(in, ArrayLine);
+
+      std::size_t q0 = ArrayLine.find('"');
+      std::size_t q1 = ArrayLine.rfind('"');
+
+      if (q0 != std::string::npos && q1 > q0)
      {
-         if (parse_quoted_string(ArrayName, in) != 0)
-         {
-            MFEM_ABORT("error parsing input!");
-         }
+         // Locate set name between first and last double quote
+         ArrayName = ArrayLine.substr(q0+1,q1-q0-1);
      }
      else
      {
-         in >> ArrayName;
-         MFEM_VERIFY(in.good(), "error parsing input!");
+         // If no double quotes found locate set name using white space
+         q1 = ArrayLine.find(' ');
+         ArrayName = ArrayLine.substr(0,q1-1);
+      }
+      if (q1+2 < ArrayLine.size())
+      {
+         // Read the remainder of the line which contains the array data
+         std::istringstream ArrayDataStream(ArrayLine.substr(q1+2,
+                                                             ArrayLine.size()));
+         data[ArrayName].Load(ArrayDataStream, 0);
+      }
+      else
+      {
+         // Read the array data starting on the next line
+         data[ArrayName].Load(in, 0);
      }
-
-      // Read the array
-      data[ArrayName].Load(in);
   }
+
 }

 }
@@ -1090,12 +1090,6 @@ 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)
 {
@@ -50,48 +50,6 @@ inline void filter_dos(std::string &line)
   }
 }

-/** @brief Read a string formatted using std::quoted. Return nonzero on error.
-
-    The stream @a in must begin with @a delim. After clearing @a result and
-    extracting the opening @a delim, characters are extracted from @a in and
-    processed as follows:
-    - if the character is @a delim, return 0;
-    - if the character is different from @a escape, it is appended to @a result;
-    - if the character is @a escape, the next character from @a in is extracted
-      and if it is one of @a delim or @a escape, it is appended to @a result;
-      otherwise, both @a escape and the character after it are appended to
-      @a result; note that the latter case is not possible if the input was
-      formatted with std::quoted with the same @a delim and @a escape
-      characters.
-
-    If the stream @a in does not begin with @a delim, error code 1 is returned.
-    If reading the stream fails, error code 2 is returned. On success, zero is
-    returned and the closing @a delim character is the last character extracted
-    from @a in. */
-inline int parse_quoted_string(std::string &result, std::istream &in,
-                               char delim = '"', char escape = '\\')
-{
-   using tt = std::string::traits_type;  // std::char_traits<char>
-   auto equal = [](tt::int_type c1, tt::char_type c2) -> bool
-   {
-      return tt::eq_int_type(c1, tt::to_int_type(c2));
-   };
-   result.clear();
-   if (!equal(in.peek(), delim)) { return 1; }
-   in.get();  // extract delim
-   for (auto c = in.get(); !equal(c, delim); c = in.get())
-   {
-      if (equal(c, escape))
-      {
-         c = in.get();
-         if (!equal(c, escape) && !equal(c, delim)) { result += escape; }
-      }
-      if (!in) { return 2; }
-      result += tt::to_char_type(c);
-   }
-   return 0;
-}
-
 /// Convert an integer to a 0-padded string with the given number of @a digits
 inline std::string to_padded_string(int i, int digits)
 {
@@ -160,9 +160,6 @@ const char *GetConfigStr()
 #ifdef MFEM_USE_RAJA
      "MFEM_USE_RAJA\n"
 #endif
-#ifdef MFEM_USE_SIDRE
-      "MFEM_USE_SIDRE\n"
-#endif
 #ifdef MFEM_USE_SIMD
      "MFEM_USE_SIMD\n"
 #endif
@@ -123,20 +123,15 @@ EXAMPLE_SUBDIRS = amgx caliper ginkgo hiop petsc pumi sundials superlu moonolith
 EXAMPLE_DIRS := examples $(addprefix examples/,$(EXAMPLE_SUBDIRS))
 EXAMPLE_TEST_DIRS := examples

-MINIAPP_ALL_SUBDIRS = common electromagnetics meshing performance tools \
+MINIAPP_SUBDIRS = common electromagnetics meshing performance tools \
 toys nurbs gslib adjoint solvers shifted mtop parelag tribol autodiff dfem \
 hooke multidomain dpg hdiv-linear-solver spde diag-smoothers contact \
- fluids/navier fluids/schrodinger-flow plasma plasma/pic
-MINIAPP_RECURSIVE_SUBDIRS = plasma/pic
-MINIAPP_SUBDIRS := $(filter-out \
-   $(MINIAPP_RECURSIVE_SUBDIRS),$(MINIAPP_ALL_SUBDIRS))
-MINIAPP_ALL_DIRS := $(addprefix miniapps/,$(MINIAPP_ALL_SUBDIRS))
+ fluids/navier fluids/schrodinger-flow plasma
 MINIAPP_DIRS := $(addprefix miniapps/,$(MINIAPP_SUBDIRS))
 MINIAPP_TEST_DIRS := $(filter-out %/common,$(MINIAPP_DIRS))
 MINIAPP_USE_COMMON := $(addprefix miniapps/,electromagnetics meshing tools \
- toys gslib shifted dpg diag-smoothers fluids/navier plasma plasma/pic)
+ toys shifted dpg diag-smoothers fluids/navier plasma)

-EM_ALL_DIRS = $(EXAMPLE_DIRS) $(MINIAPP_ALL_DIRS)
 EM_DIRS = $(EXAMPLE_DIRS) $(MINIAPP_DIRS)

 TEST_SUBDIRS = unit
@@ -151,7 +146,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_ALL_DIRS) doc $(TEST_DIRS)
+   BUILD_SUBDIRS = $(DIRS) config $(EM_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))
@@ -304,7 +299,7 @@ ifeq ($(MFEM_USE_LEGACY_OPENMP),YES)
 endif

 # List of MFEM dependencies, that require the *_LIB variable to be non-empty
-MFEM_REQ_LIB_DEPS = SUPERLU MUMPS METIS FMS CONDUIT SIDRE LAPACK SUNDIALS\
+MFEM_REQ_LIB_DEPS = SUPERLU MUMPS METIS FMS CONDUIT LAPACK SUNDIALS\
 SUITESPARSE STRUMPACK GINKGO GNUTLS HDF5 NETCDF SLEPC PETSC MPFR PUMI HIOP\
 GSLIB OCCA CEED RAJA UMPIRE MKL_CPARDISO MKL_PARDISO AMGX MAGMA CALIPER PARELAG\
 TRIBOL BENCHMARK MOONOLITH ALGOIM
@@ -370,7 +365,7 @@ MFEM_DEFINES = MFEM_VERSION MFEM_VERSION_STRING MFEM_GIT_STRING MFEM_USE_MPI\
 MFEM_USE_LEGACY_OPENMP MFEM_USE_MEMALLOC MFEM_TIMER_TYPE MFEM_USE_SUNDIALS\
 MFEM_USE_SUITESPARSE MFEM_USE_GINKGO MFEM_USE_SUPERLU MFEM_USE_SUPERLU5\
 MFEM_USE_STRUMPACK MFEM_USE_GNUTLS MFEM_USE_HDF5 MFEM_USE_NETCDF MFEM_USE_PETSC\
- MFEM_USE_SLEPC MFEM_USE_MPFR MFEM_USE_SIDRE MFEM_USE_FMS MFEM_USE_CONDUIT\
+ MFEM_USE_SLEPC MFEM_USE_MPFR MFEM_USE_FMS MFEM_USE_CONDUIT\
 MFEM_USE_PUMI MFEM_USE_HIOP MFEM_USE_GSLIB MFEM_USE_CUDA MFEM_USE_HIP\
 MFEM_USE_OCCA MFEM_USE_MOONOLITH MFEM_USE_CEED MFEM_USE_RAJA MFEM_USE_UMPIRE\
 MFEM_USE_SIMD MFEM_USE_ADIOS2 MFEM_USE_MKL_CPARDISO MFEM_USE_MKL_PARDISO MFEM_USE_AMGX\
@@ -488,10 +483,10 @@ $(OBJECT_FILES): $(BLD)%.o: $(SRC)%.cpp $(CONFIG_MK)

 all: examples miniapps $(TEST_DIRS)

-.PHONY: miniapps $(EM_ALL_DIRS) $(TEST_DIRS)
+.PHONY: miniapps $(EM_DIRS) $(TEST_DIRS)
 miniapps: $(MINIAPP_DIRS)
 $(MINIAPP_USE_COMMON): miniapps/common
-$(EM_ALL_DIRS) $(TEST_DIRS): lib
+$(EM_DIRS) $(TEST_DIRS): lib
 	$(MAKE) -C $(BLD)$(@)

 .PHONY: doc
@@ -699,7 +694,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_ALL_DIRS) doc $(TEST_DIRS)); do\
+	for dir in "" $(addsuffix /,config $(EM_DIRS) doc $(TEST_DIRS)); do \
 	   printf "# Auto-generated file.\n%s\n%s\n" \
 	      "MFEM_DIR = $(MFEM_REAL_DIR)" \
 	      "include \$$(MFEM_DIR)/$${dir}makefile" \
@@ -751,7 +746,6 @@ status info:
 	$(info MFEM_USE_PETSC         = $(MFEM_USE_PETSC))
 	$(info MFEM_USE_SLEPC         = $(MFEM_USE_SLEPC))
 	$(info MFEM_USE_MPFR          = $(MFEM_USE_MPFR))
-	$(info MFEM_USE_SIDRE         = $(MFEM_USE_SIDRE))
 	$(info MFEM_USE_FMS           = $(MFEM_USE_FMS))
 	$(info MFEM_USE_CONDUIT       = $(MFEM_USE_CONDUIT))
 	$(info MFEM_USE_PUMI          = $(MFEM_USE_PUMI))
@@ -801,15 +795,13 @@ status info:

 ASTYLE = $(ASTYLE_BIN) --options=$(SRC)config/mfem.astylerc
 ASTYLE_VER = "Artistic Style Version 3.1"
-FORMAT_FILES = $(foreach dir,$(DIRS) $(EM_ALL_DIRS) config,$(dir)/*.?pp)
+FORMAT_FILES = $(foreach dir,$(DIRS) $(EM_DIRS) config,$(dir)/*.?pp)
 TESTS_SUBDIRS = unit benchmarks convergence mem_manager par-mesh-format
-UNIT_TESTS_SUBDIRS = general linalg mesh fem miniapps ceed enzyme dfem
-MINIAPPS_SUBDIRS = dpg/util hooke/operators hooke/preconditioners \
-   hooke/materials hooke/kernels
+UNIT_TESTS_SUBDIRS = general linalg mesh fem miniapps ceed enzyme
+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)))

@@ -840,29 +832,14 @@ 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; \
-	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...";\
+	@err_code=0;\
 	$(call mfem_check_command,\
 	    $(ASTYLE) $(FORMAT_LIST) | grep Formatted,\
 	    "No source files were changed",\
@@ -1616,9 +1616,7 @@ Element::Type Mesh::GetFaceElementType(int Face) const

 Array<int> Mesh::GetFaceToBdrElMap() const
 {
-   Array<int> face_to_be(Dim == 1 ? NumOfVertices :
-                         Dim == 2 ? NumOfEdges :
-                         Dim == 3 ? NumOfFaces : 0);
+   Array<int> face_to_be(Dim == 2 ? NumOfEdges : NumOfFaces);
   face_to_be = -1;
   for (int i = 0; i < NumOfBdrElements; i++)
   {
@@ -63,6 +63,7 @@ ThresholdRefiner::ThresholdRefiner(ErrorEstimator &est)

   threshold = 0.0;
   num_marked_elements = 0LL;
+   current_sequence = -1;

   non_conforming = -1;
   nc_limit = 0;
@@ -86,6 +87,7 @@ int ThresholdRefiner::MarkWithoutRefining(Mesh & mesh,
   threshold = 0.0;
   num_marked_elements = 0LL;
   refinements.SetSize(0);
+   current_sequence = mesh.GetSequence();

   const long long num_elements = mesh.GetGlobalNE();
   if (num_elements >= max_elements) { return STOP; }
@@ -147,6 +149,7 @@ int ThresholdRefiner::ApplyImpl(Mesh &mesh)
 void ThresholdRefiner::Reset()
 {
   estimator.Reset();
+   current_sequence = -1;
   num_marked_elements = 0LL;
   // marked_elements.SetSize(0); // not necessary
 }
@@ -188,6 +188,7 @@ protected:
   long long num_marked_elements;

   Array<Refinement> marked_elements;
+   long current_sequence;

   int non_conforming;
   int nc_limit;
@@ -227,29 +227,15 @@ public:
                                           const ParGridFunction &dst);

   /**
-   * @brief Check if Mesh @a m is a ParSubMesh.
+   * @brief Check if ParMesh @a m is a ParSubMesh.
   *
-   * @param m The input Mesh
+   * @param m The input ParMesh
   */
-   static bool IsParSubMesh(const Mesh *m)
+   static bool IsParSubMesh(const ParMesh *m)
   {
      return dynamic_cast<const ParSubMesh *>(m) != nullptr;
   }

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

-   /**
-   * @brief Check if Mesh @a sub is a SubMesh of Mesh @a parent.
-   *
-   * @param sub The potential submesh Mesh
-   * @param parent The potential parent Mesh
-   */
-   static bool IsSubMesh(const Mesh* sub, const Mesh* parent)
-   {
-      while (IsSubMesh(sub) &&
-             (sub = static_cast<const SubMesh *>(sub)->GetParent()) &&
-             sub != parent);
-      return sub == parent;
-   }
-
 private:
   /// Private constructor
   SubMesh(const Mesh &parent, From from, const Array<int> &attributes);
@@ -82,6 +82,43 @@ Although Tribol can be built automatically via **uberenv** and **Spack**,
 for this miniapp it is simpler to build **Axom** and **MFEM** manually and 
 point Tribol to them. The steps are as follows:

+### Using pre-built Tribol/Axom installs
+
+If you already have compatible installs of Tribol and Axom, point MFEM to the install prefixes.
+
+- Hypre install prefix: `<path/to/hypre>`
+- METIS install prefix: `<path/to/metis>`
+- Axom install prefix: `<path/to/axom>`
+- Tribol install prefix: `<path/to/tribol>`
+
+**MFEM make build (configure):**
+```bash
+make config MFEM_USE_MPI=YES MFEM_USE_METIS=YES MFEM_USE_TRIBOL=YES \
+  HYPRE_DIR=<path/to/hypre> METIS_DIR=<path/to/metis> \
+  AXOM_DIR=<path/to/axom> TRIBOL_DIR=<path/to/tribol> ADIAK_DIR=<path/to/adiak> CAMP_DIR=<path/to/camp> RAJA_DIR=<path/to/raja> \
+  UMPIRE_DIR=<path/to/umpire> FMT_DIR=<path/to/fmt> CALIPER_DIR=<path/to/caliper>
+```
+
+**MFEM CMake build (configure):**
+```bash
+cmake -S . -B <mfem-build-dir> -DMFEM_USE_MPI=YES -DMFEM_USE_METIS=YES -DMFEM_USE_TRIBOL=YES \
+  HYPRE_DIR=<path/to/hypre> METIS_DIR=<path/to/metis> \
+  AXOM_DIR=<path/to/axom> TRIBOL_DIR=<path/to/tribol> ADIAK_DIR=<path/to/adiak> CAMP_DIR=<path/to/camp> RAJA_DIR=<path/to/raja> \
+  UMPIRE_DIR=<path/to/umpire> FMT_DIR=<path/to/fmt> CALIPER_DIR=<path/to/caliper>
+```
+
+Note: RAJA/UMPIRE/CALIPER are optional for MFEM itself, but many Tribol builds
+enable them. If your Tribol install does not depend on them, you can omit the
+corresponding `*_DIR` entries above.
+
+Note: `FMT_DIR` only needs to be added for the make-based build (and only when
+the Umpire install uses `fmt`). If `FMT_DIR` is not set and a sibling `fmt-*`
+directory exists next to your `UMPIRE_DIR`, MFEM's make configuration will try
+to pick it up automatically.
+
+Note: when using pre-built Tribol/Axom, you typically need to use a compatible
+compiler/MPI wrapper (same C++ standard library ABI).
+
 ### Manual Build Steps

 1. Pull axom and tribol (starting from the mfem folder):
@@ -99,7 +136,7 @@ point Tribol to them. The steps are as follows:
   TRIBOL_DIR = @MFEM_DIR@/../tribol-repo/tribol
   TRIBOL_OPT = -I$(TRIBOL_DIR)/include -I$(AXOM_DIR)/include
   TRIBOL_LIB = -L$(TRIBOL_DIR)/lib -ltribol -lredecomp -L$(AXOM_DIR)/lib \
-   -laxom_mint -laxom_slam -laxom_slic -laxom_core
+   -laxom_quest -laxom_mint -laxom_slam -laxom_slic -laxom_lumberjack -laxom_core
   ```
 3. [**Axom:**](https://github.com/LLNL/axom.git) Starting from the MFEM root
   directory (we assume this directory is named mfem):
@@ -43,39 +43,19 @@ endif()
  
  # Add the corresponding tests to the "test" target
  if (MFEM_ENABLE_TESTING)
-    add_test(NAME tesla_1_np=${MFEM_MPI_NP}
+    add_test(NAME tesla_np=4
      COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
      ${MPIEXEC_PREFLAGS}
      $<TARGET_FILE:tesla> -no-vis -maxit 2 -cr "0 0 -0.2 0 0 0.2 0.2 0.4 1"
      ${MPIEXEC_POSTFLAGS})
-    add_test(NAME tesla_2_np=${MFEM_MPI_NP}
+
+    add_test(NAME volta_np=4
      COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
      ${MPIEXEC_PREFLAGS}
-      $<TARGET_FILE:tesla>
-      -no-vis -maxit 2 -m ../../data/inline-hex.mesh -ubbc "0 0 1"
+      $<TARGET_FILE:volta> -no-vis -maxit 2 -dbcs 1 -dbcg -ds "0.0 0.0 0.0 0.2 8.0"
      ${MPIEXEC_POSTFLAGS})

-    add_test(NAME volta_1_np=${MFEM_MPI_NP}
-      COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
-      ${MPIEXEC_PREFLAGS}
-      $<TARGET_FILE:volta>
-      -no-vis -maxit 2 -dbcs 1 -dbcg -ds "0.0 0.0 0.0 0.2 8.0"
-      ${MPIEXEC_POSTFLAGS})
-    add_test(NAME volta_2_np=${MFEM_MPI_NP}
-      COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
-      ${MPIEXEC_PREFLAGS}
-      $<TARGET_FILE:volta>
-      -no-vis -maxit 2 -m ../../data/square-disc.mesh -dbcs "1 2 3 4 5 6 7 8"
-      -dbcv "0 0 0 0 1 1 1 1"
-      ${MPIEXEC_POSTFLAGS})
-    add_test(NAME volta_3_np=${MFEM_MPI_NP}
-      COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
-      ${MPIEXEC_PREFLAGS}
-      $<TARGET_FILE:volta>
-      -no-vis -maxit 2 -m ../../data/inline-hex.mesh -dbcs "1 6" -dbcv "0 1"
-      ${MPIEXEC_POSTFLAGS})
-
-    add_test(NAME joule_np=${MFEM_MPI_NP}
+    add_test(NAME joule_np=4
      COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
      ${MPIEXEC_PREFLAGS}
      $<TARGET_FILE:joule>
@@ -83,41 +63,12 @@ endif()
      ${MPIEXEC_POSTFLAGS})

    if (MFEM_USE_DOUBLE) # otherwise returns MFEM_SKIP_RETURN_VALUE
-      add_test(NAME maxwell_np=${MFEM_MPI_NP}
+      add_test(NAME maxwell_np=4
        COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
        ${MPIEXEC_PREFLAGS}
        $<TARGET_FILE:maxwell>
        -no-vis -abcs "-1" -dp "-0.3 0.0 0.0 0.3 0.0 0.0 0.1 1 .5 .5"
        ${MPIEXEC_POSTFLAGS})
    endif()
-
-    if (MFEM_USE_GSLIB)
-      add_test(NAME lorentz_1_np=${MFEM_MPI_NP}
-        COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
-        ${MPIEXEC_PREFLAGS}
-        $<TARGET_FILE:lorentz>
-        -no-vis -er Volta-AMR-Parallel -ec 2 -npt 100 -xmin "0.0 0.0 0.0"
-        -xmax "1.0 1.0 1.0" -pmin "1 0 0" -pmax "1 0 0" -rdf 0 -vt 0 -nt 100
-        ${MPIEXEC_POSTFLAGS})
-      # Setup dependency on volta_3_np=<np>
-      set_tests_properties(volta_3_np=${MFEM_MPI_NP}
-        PROPERTIES FIXTURES_SETUP Volta3)
-      set_tests_properties(lorentz_1_np=${MFEM_MPI_NP}
-        PROPERTIES FIXTURES_REQUIRED Volta3)
-
-      add_test(NAME lorentz_2_np=${MFEM_MPI_NP}
-        COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
-        ${MPIEXEC_PREFLAGS}
-        $<TARGET_FILE:lorentz>
-        -no-vis -br Tesla-AMR-Parallel -bc 2 -npt 10 -xmin "0.0 0.0 0.0"
-        -xmax "1.0 1.0 1.0" -pmin "0 0.1 0.05" -pmax "0 0.4 0.1" -nt 1000 -rdf 0
-        -vt 0
-        ${MPIEXEC_POSTFLAGS})
-      # Setup dependency on tesla_2_np=<np>
-      set_tests_properties(tesla_2_np=${MFEM_MPI_NP}
-        PROPERTIES FIXTURES_SETUP Tesla2)
-      set_tests_properties(lorentz_2_np=${MFEM_MPI_NP}
-        PROPERTIES FIXTURES_REQUIRED Tesla2)
-    endif()
  endif()
 endif()
@@ -117,10 +117,10 @@ joule-test-par: joule
 lorentz-test-par: lorentz-test-1 lorentz-test-2
 lorentz-test-1: lorentz volta-test-3
 	@$(call mfem-test,$<, $(RUN_MPI), Electromagnetic miniapp,\
-	-er Volta-AMR-Parallel -ec 2 -npt 100 -xmin '0.0 0.0 0.0' -xmax '1.0 1.0 1.0' -pmin '1 0 0' -pmax '1 0 0' -rdf 0 -vt 0 -nt 100)
+	-er Volta-AMR-Parallel -ec 2 -npt 100 -xmin '0.0 0.0 0.0' -xmax '1.0 1.0 1.0' -pmin '1 0 0' -pmax '1 0 0' -rdf 0 -vt 0 -nt 100')
 lorentz-test-2: lorentz tesla-test-2
 	@$(call mfem-test,$<, $(RUN_MPI), Electromagnetic miniapp,\
-	-br Tesla-AMR-Parallel -bc 2 -npt 10 -xmin '0.0 0.0 0.0' -xmax '1.0 1.0 1.0' -pmin '0 0.1 0.05' -pmax '0 0.4 0.1' -nt 1000 -rdf 0 -vt 0)
+	-br Tesla-AMR-Parallel -bc 2 -br Tesla-AMR-Parallel -npt 10 -xmin '0.0 0.0 0.0' -xmax '1.0 1.0 1.0' -pmin '0 0.1 0.05' -pmax '0 0.4 0.1' -nt 1000 -rdf 0 -vt 0)

 # Testing: "test" target and mfem-test* variables are defined in config/test.mk

@@ -46,16 +46,8 @@ if (MFEM_USE_GSLIB)

  if (MFEM_ENABLE_TESTING)
    foreach (test "schwarz_ex1" "field-diff" "findpts" "field-interp")
-      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()
+      add_test(NAME ${test}
+        COMMAND $<TARGET_FILE:${test}> -no-vis)
    endforeach()
  endif()

@@ -83,16 +83,11 @@ 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,$<, $(RUN_MPI_SER), $(TEST_NAME))
+	@$(call mfem-test,$<,, $(TEST_NAME))

 # Testing: "test" target and mfem-test* variables are defined in config/test.mk

@@ -22,7 +22,7 @@ void ComputeInverse(const Array<real_t> &A, Array<real_t> &Ainv)
 {
   Array<real_t> A2 = A;
   const int n2 = A.Size();
-   const int n = static_cast<int>(sqrt(n2));
+   const int n = static_cast<const int>(sqrt(n2));
   Array<int> ipiv(n);
   LUFactors lu(A2.GetData(), ipiv.GetData());
   lu.Factor(n);
@@ -58,7 +58,7 @@ void SubcellIntegrals(int n, const Poly_1D::Basis &basis, Array<real_t> &B)

 void Transpose(const Array<real_t> &B, Array<real_t> &Bt)
 {
-   const int n = static_cast<int>(sqrt(B.Size()));
+   const int n = static_cast<const int>(sqrt(B.Size()));
   Bt.SetSize(n*n);
   for (int i=0; i<n; ++i) for (int j=0; j<n; ++j) { Bt[i+j*n] = B[j+i*n]; }
 }
@@ -23,5 +23,3 @@ if (MFEM_USE_MPI)
      EXTRA_HEADERS ${PLASMA_COMMON_HEADERS})

 endif()
-
-add_subdirectory(pic)
@@ -14,6 +14,9 @@ MFEM_DIR ?= ../..
 MFEM_BUILD_DIR ?= ../..
 SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/plasma/,)
 CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
+# Use the MFEM install directory
+# MFEM_INSTALL_DIR = ../../mfem
+# CONFIG_MK = $(MFEM_INSTALL_DIR)/share/mfem/config.mk

 MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
@@ -26,14 +29,6 @@ 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
@@ -43,24 +38,26 @@ 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))
+   $(if $(MFEM_USE_CUDA:YES=),$(CXX_XLINKER),$(CUDA_XLINKER))-rpath,$(abspath\
+   $(MFEM_BUILD_DIR)/miniapps/common))
+
+COMMON_O=

 # Remove built-in rules
 %: %.cpp
 %.o: %.cpp

-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)
+all: $(MINIAPPS)

 # Rules for building the miniapps
-%: $(SRC)%.cpp $(MFEM_LIB_FILE) $(CONFIG_MK) | lib-common
-	$(MFEM_CXX) $(MFEM_LINK_FLAGS) $< -o $@ $(COMMON_LIB) $(MFEM_LIBS)
+%: $(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 $(@)

 # Rule for building lib-common
 lib-common:
@@ -68,9 +65,6 @@ 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)
@@ -81,10 +75,11 @@ 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 $(SUBDIRS_CLEAN)
+clean: clean-build clean-exec

 clean-build:
 	rm -f *.o *~ $(SEQ_MINIAPPS) $(PAR_MINIAPPS)
 	rm -rf *.dSYM *.TVD.*breakpoints

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

@@ -18,7 +18,6 @@
 //
 // Currently supported data collection type options:
 //    visit:                VisItDataCollection (default)
-//    sidre or sidre_hdf5:  SidreDataCollection
 //    json:                 ConduitDataCollection w/ protocol json
 //    conduit_json:         ConduitDataCollection w/ protocol conduit_json
 //    conduit_bin:          ConduitDataCollection w/ protocol conduit_bin
@@ -52,14 +51,6 @@ DataCollection *create_data_collection(const std::string &dc_name,
      dc = new VisItDataCollection(MPI_COMM_WORLD, dc_name);
 #else
      dc = new VisItDataCollection(dc_name);
-#endif
-   }
-   else if ( dc_type == "sidre" || dc_type == "sidre_hdf5")
-   {
-#ifdef MFEM_USE_SIDRE
-      dc = new SidreDataCollection(dc_name);
-#else
-      MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for sidre support.");
 #endif
   }
   else if ( dc_type == "json" ||
@@ -140,7 +131,6 @@ int main(int argc, char *argv[])
   args.AddOption(&src_coll_type, "-st", "--source-type",
                  "Set the source data collection type. Options:\n"
                  "\t   visit:                VisItDataCollection (default)\n"
-                  "\t   sidre or sidre_hdf5:  SidreDataCollection\n"
                  "\t   json:                 ConduitDataCollection w/ protocol json\n"
                  "\t   conduit_json:         ConduitDataCollection w/ protocol conduit_json\n"
                  "\t   conduit_bin:          ConduitDataCollection w/ protocol conduit_bin\n"
@@ -152,7 +142,6 @@ int main(int argc, char *argv[])
   args.AddOption(&out_coll_type, "-ot", "--output-type",
                  "Set the output data collection type. Options:\n"
                  "\t   visit:                VisItDataCollection (default)\n"
-                  "\t   sidre or sidre_hdf5:  SidreDataCollection\n"
                  "\t   json:                 ConduitDataCollection w/ protocol json\n"
                  "\t   conduit_json:         ConduitDataCollection w/ protocol conduit_json\n"
                  "\t   conduit_bin:          ConduitDataCollection w/ protocol conduit_bin\n"
@@ -23,8 +23,7 @@ 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 \
-	random-gridfunction-bounds
+PAR_MINIAPPS = nodal-transfer plor-transfer gridfunction-bounds

 ifeq ($(MFEM_USE_MPI),NO)
   MINIAPPS = $(SEQ_MINIAPPS)
@@ -80,7 +79,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 \
-	random-gridfunction-bounds compare-dc
+	compare-dc
 $(foreach app,$(NO_TEST_APPS),$(app)-test-seq $(app)-test-par):
 	@true

@@ -1,263 +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.
-//
-//    ---------------------------------------------------------------------
-//     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;
-   }
-}
@@ -76,7 +76,6 @@ set(UNIT_TESTS_SRCS
  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
@@ -1,118 +0,0 @@
-MFEM mesh v1.3
-
-#
-# MFEM Geometry Types (see mesh/geom.hpp):
-#
-# POINT       = 0
-# SEGMENT     = 1
-# TRIANGLE    = 2
-# SQUARE      = 3
-# TETRAHEDRON = 4
-# CUBE        = 5
-# PRISM       = 6
-#
-
-dimension
-2
-
-elements
-12
-10 2 7 0 1
-11 2 0 7 2
-12 2 9 0 2
-13 2 0 9 3
-14 2 11 0 3
-15 2 0 11 4
-16 2 5 0 4
-17 2 0 5 1
-9 3 1 5 6 7
-9 3 2 7 8 9
-9 3 3 9 10 11
-9 3 4 11 12 5
-
-attribute_sets
-16
-"Base" 1 9
-"E Even" 1 16
-"E Odd" 1 17
-"East"
-2
-16
-17
-"N Even" 1 10
-"N Odd" 1 11
-"North" 2 10     11
-"Rose" 8 10 11     12
-13    14 
-15    16     17
-"Rose Even" 4
-10
-12
-14
-16
-"Rose Odd"
-4
-11
-13
-15
-17
-"S Even" 1 14
-"S Odd" 1 15
-South 2
-14
-15
-"W Even" 1 12
-"W Odd" 1 13
-West 2 12 13
-
-boundary
-8
-1 1 5 6
-2 1 6 7
-3 1 7 8
-4 1 8 9
-5 1 9 10
-6 1 10 11
-7 1 11 12
-8 1 12 5
-
-bdr_attribute_sets
-13
-"Boundary" 8 1 2 3 4 5 6 7 8
-"ENE" 1 1
-"ESE" 1 8
-"Eastern Boundary" 2 1 8
-"NNE" 1 2
-"NNW" 1 3
-"Northern Boundary"
-2
-2
-3
-"SSE" 1 7
-"SSW" 1 6
-"Southern Boundary" 2
-6
-7
-"WNW" 1 4
-"WSW" 1 5
-"Western Boundary" 2 4
-5
-
-vertices
-13
-2
-0 0
-0.14142136 0.14142136
-0.14142136 0.14142136
-0.14142136 -0.14142136
-0.14142136 -0.14142136
-1 0
-0.70710678 0.70710678
-0 1
-0.70710678 0.70710678
-1 0
-0.70710678 -0.70710678
-0 -1
-0.70710678 -0.70710678
-
-mfem_mesh_end
@@ -296,7 +296,7 @@ void TestRedistribute(Ordering::Type ordering)
      int wrong_proc_count = 0;
      for (int i = 0; i < procs.Size(); i++)
      {
-         if (static_cast<unsigned>(rank) != procs[i])
+         if (rank != procs[i])
         {
            wrong_proc_count++;
         }
@@ -30,29 +30,3 @@ TEST_CASE("String Manipulation", "[General]")
      }
   }
 }
-
-TEST_CASE("Quoted String Input", "[General]")
-{
-   const auto test_strings =
-   {
-      "Test",
-      "Test with spaces",
-      "Test with \"quoted text\"",
-      "Test string ending with \\",
-      "\nTest with\tvarious white\v\rspace characters.",
-      "Test with some unicode characters: ∆, ∉, ∑, 🍎."
-   };
-
-   for (const auto c_str : test_strings)
-   {
-      CAPTURE(c_str);
-      const std::string str(c_str);
-      std::stringstream ss;
-      ss << std::quoted(str);
-
-      std::string read_str;
-      int error = parse_quoted_string(read_str, ss);
-      CHECK(error == 0);
-      CHECK(read_str == str);
-   }
-}
@@ -1,108 +0,0 @@
-// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
-// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-// LICENSE and NOTICE for details. LLNL-CODE-806117.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability visit https://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the BSD-3 license. We welcome feedback and contributions, see file
-// CONTRIBUTING.md for details.
-
-#include "mfem.hpp"
-#include "unit_tests.hpp"
-
-#include <algorithm>
-#include <string>
-#include <utility>
-#include <vector>
-
-using namespace mfem;
-
-TEST_CASE("MFEM Mesh Named Attributes", "[Mesh]")
-{
-   // Path relative to the directory tests/unit
-   Mesh mesh("data/compass-testing.mesh");
-
-   REQUIRE(mesh.Dimension() == 2);
-   REQUIRE(mesh.GetNE() == 12);
-   REQUIRE(mesh.GetNV() == 13);
-
-   REQUIRE(mesh.attribute_sets.attr_sets.Size() == 16);
-   REQUIRE(mesh.bdr_attribute_sets.attr_sets.Size() == 13);
-
-   std::vector<std::pair<std::string, std::vector<int>>> expected_attr_sets =
-   {
-      {"Base",  {9}},
-      {"E Even",  {16}},
-      {"E Odd",  {17}},
-      {"East",  {16, 17}},
-      {"N Even",  {10}},
-      {"N Odd",  {11}},
-      {"North",  {10, 11}},
-      {"Rose",  {10, 11, 12, 13, 14, 15, 16, 17}},
-      {"Rose Even",  {10, 12, 14, 16}},
-      {"Rose Odd",  {11, 13, 15, 17}},
-      {"S Even",  {14}},
-      {"S Odd",  {15}},
-      {"South",  {14, 15}},
-      {"W Even",  {12}},
-      {"W Odd",  {13}},
-      {"West",  {12, 13}}
-   };
-
-   for (auto const &attr_name_index_pair: expected_attr_sets )
-   {
-      REQUIRE(mesh.attribute_sets.AttributeSetExists(
-                 attr_name_index_pair.first));
-
-      auto const &attr_set = mesh.attribute_sets.GetAttributeSet(
-                                attr_name_index_pair.first);
-      auto const &expected_attr_set = attr_name_index_pair.second;
-
-      REQUIRE(static_cast<std::size_t>(attr_set.Size()) ==
-              expected_attr_set.size());
-
-      bool const elements_equal = std::equal(attr_set.begin(), attr_set.end(),
-                                             expected_attr_set.begin());
-
-      REQUIRE(elements_equal);
-   }
-
-   std::vector<std::pair<std::string, std::vector<int>>> expected_bdr_attr_sets
-   =
-   {
-      {"Boundary",  {1, 2, 3, 4, 5, 6, 7, 8}},
-      {"ENE",  { 1}},
-      {"ESE",  { 8}},
-      {"Eastern Boundary",  {1, 8}},
-      {"NNE",  { 2}},
-      {"NNW",  { 3}},
-      {"Northern Boundary", {2, 3}},
-      {"SSE",  { 7}},
-      {"SSW",  { 6}},
-      {"Southern Boundary", {6,7}},
-      {"WNW",  { 4}},
-      {"WSW",  { 5}},
-      {"Western Boundary", {4,5}}
-   };
-
-   for (auto const &attr_bdr_name_index_pair: expected_bdr_attr_sets )
-   {
-      REQUIRE(mesh.bdr_attribute_sets.AttributeSetExists(
-                 attr_bdr_name_index_pair.first));
-
-      auto const &bdr_attr_set = mesh.bdr_attribute_sets.GetAttributeSet(
-                                    attr_bdr_name_index_pair.first);
-      auto const &expected_bdr_attr_set = attr_bdr_name_index_pair.second;
-
-      REQUIRE(static_cast<std::size_t>(bdr_attr_set.Size()) ==
-              expected_bdr_attr_set.size());
-
-      bool const elements_equal = std::equal(bdr_attr_set.begin(),
-                                             bdr_attr_set.end(),
-                                             expected_bdr_attr_set.begin());
-
-      REQUIRE(elements_equal);
-   }
-}
Author	SHA1	Message	Date
Kenneth Weiss	b86fb5308a	style	2026-04-14 08:54:23 -07:00
Kenneth Weiss	4463d921f6	Updates CHANGELOG	2026-04-13 21:52:37 -07:00
Kenneth Weiss	20c657559d	More fixes for building contact mini-app with make and cmake	2026-04-13 21:47:15 -07:00
Kenneth Weiss	7a999efb0d	yAdds back FindAxom in support of contact mini-app	2026-04-13 18:54:17 -07:00
Kenneth Weiss	442ac18b65	Removes references to Axom in INSTALL	2026-04-13 12:40:24 -07:00
Kenneth Weiss	5097a44411	Removes Axom from make-based and cmake-based build systems	2026-04-13 12:37:18 -07:00
Kenneth Weiss	e84fef4191	Removes SidreDataCollection from MFEM This has been superseded by the MFEMSidreDataCollection in Axom.	2026-04-13 12:34:02 -07:00