Fix style error

Remove unused variable from test
Apply style to tensor.hpp
2026-03-29 17:12:36 -07:00 · 2026-03-27 21:08:43 -07:00 · 2026-03-27 21:01:16 -07:00 · 2026-03-27 20:56:18 -07:00 · 2026-03-27 20:32:59 -07:00 · 2026-03-27 20:02:46 -07:00
104 changed files with 1813 additions and 6420 deletions
@@ -25,7 +25,7 @@ runs:
  steps:
    - uses: ./.github/actions/sanitize/config

-    - uses: actions/cache@v5
+    - uses: actions/cache@v4
      if: ${{env.DEBUG == 'true'}}
      id: debug
      with:
@@ -36,7 +36,7 @@ runs:
  steps:
    - uses: ./.github/actions/sanitize/config

-    - uses: actions/cache@v5
+    - uses: actions/cache@v4
      if: ${{env.DEBUG == 'true' && inputs.cache-skip != 'true'}}
      id: debug
      with:
@@ -23,7 +23,7 @@ inputs:
 runs:
  using: 'composite'
  steps:
-    - uses: actions/cache/restore@v5  # Cache for LLVM libcxx
+    - uses: actions/cache/restore@v4  # Cache for LLVM libcxx
      with:
        path: ${{env.LLVM_DIR}}
        fail-on-cache-miss: true
@@ -32,14 +32,14 @@ runs:
    - uses: ./.github/actions/sanitize/mpi
      if: ${{inputs.par == 'true'}}

-    - uses: actions/cache/restore@v5  # Cache for Hypre
+    - uses: actions/cache/restore@v4  # Cache for Hypre
      if: ${{inputs.par == 'true'}}
      with:
        path: ${{env.HYPRE_DIR}}
        fail-on-cache-miss: true
        key: ${{runner.os}}-ompi-build-${{env.HYPRE_DIR}}-int32-fp64-v2.5

-    - uses: actions/cache/restore@v5  # Cache for Metis
+    - uses: actions/cache/restore@v4  # Cache for Metis
      if: ${{inputs.par == 'true'}}
      with:
        path: ${{env.METIS_DIR}}
@@ -51,13 +51,13 @@ runs:
      run: ln -s -f ${{env.HYPRE_DIR}} hypre && ln -s -f ${{env.METIS_DIR}} metis-4.0
      shell: bash

-    - uses: actions/cache/restore@v5  # Cache for LSAN suppression file
+    - uses: actions/cache/restore@v4  # Cache for LSAN suppression file
      with:
        path: ${{env.LSAN_DIR}}
        fail-on-cache-miss: true
        key: build-lsan-suppression-file

-    - uses: actions/checkout@v6  # Checkout the repository
+    - uses: actions/checkout@v4  # Checkout the repository
      with:
        path: mfem
        # ref: ${{env.BRANCH}}
@@ -43,7 +43,7 @@ jobs:
          remove-docker-images: 'true'

      - name: Checkout
-        uses: actions/checkout@v6
+        uses: actions/checkout@v4

      # It's easier to reference named variables than indexes of the matrix
      - name: Set Environment
@@ -153,7 +153,7 @@ jobs:
      # /home/runner/work/mfem/mfem/mfem
      # Note: Done now to access "install-hypre" and "install-metis" actions.
      - name: checkout mfem
-        uses: actions/checkout@v6
+        uses: actions/checkout@v4
        with:
          path: ${{ env.MFEM_TOP_DIR }}
          # Fetch the complete history for codecov to access commits ID
@@ -225,7 +225,7 @@ jobs:
      - name: cache hypre
        id: hypre-cache
        if: matrix.mpi == 'par'
-        uses: actions/cache@v5
+        uses: actions/cache@v4
        with:
          path: ${{ env.HYPRE_TOP_DIR }}
          key: ${{ runner.os }}-ompi-build-${{ env.HYPRE_TOP_DIR }}-${{ matrix.hypre-target }}-${{ matrix.precision }}-v2.5
@@ -255,7 +255,7 @@ jobs:
      - name: cache metis
        id: metis-cache
        if: matrix.mpi == 'par' && matrix.os != 'windows-latest'
-        uses: actions/cache@v5
+        uses: actions/cache@v4
        with:
          path: ${{ env.METIS_TOP_DIR }}
          key: ${{ runner.os }}-build-${{ env.METIS_TOP_DIR }}-v2.5
@@ -270,7 +270,7 @@ jobs:
      - name: cache vcpkg (Windows)
        id: vcpkg-cache
        if: matrix.os == 'windows-latest'
-        uses: actions/cache@v5
+        uses: actions/cache@v4
        with:
          path: vcpkg_cache
          key: ${{ runner.os }}-${{ matrix.mpi }}-vcpkg-v1
@@ -295,8 +295,7 @@ jobs:
          export HOMEBREW_NO_INSTALL_CLEANUP=1
          brew update
          brew install enzyme
-          ENZYME_LLVM=$(brew info enzyme | sed -n 's/^Required.*:.*\(llvm[^ ]*\).*/\1/p')
-          echo "ENZYME_LLVM=$ENZYME_LLVM"
+          ENZYME_LLVM=$(brew info enzyme | sed -n 's/^Required:.*\(llvm[^ ]*\).*/\1/p')
          LLVM_PREFIX=$(brew --prefix $ENZYME_LLVM)
          echo "LLVM_PREFIX=$LLVM_PREFIX" >> $GITHUB_ENV
          echo "OMPI_CC=$LLVM_PREFIX/bin/clang" >> $GITHUB_ENV
@@ -40,11 +40,11 @@ jobs:

    steps:
      - name: Checkout repository
-        uses: actions/checkout@v6
+        uses: actions/checkout@v4

      # Initializes the CodeQL tools for scanning.
      - name: Initialize CodeQL
-        uses: github/codeql-action/init@v4
+        uses: github/codeql-action/init@v2
        with:
          languages: ${{ matrix.language }}
          # If you wish to specify custom queries, you can do so here or in a config file.
@@ -57,7 +57,7 @@ jobs:
      # Autobuild attempts to build any compiled languages  (C/C++, C#, or Java).
      # If this step fails, then you should remove it and run the build manually (see below)
      - name: Autobuild
-        uses: github/codeql-action/autobuild@v4
+        uses: github/codeql-action/autobuild@v2

      # ℹ️ Command-line programs to run using the OS shell.
      # 📚 See https://docs.github.com/en/actions/using-workflows/workflow-syntax-for-github-actions#jobsjob_idstepsrun
@@ -70,4 +70,4 @@ jobs:
      #   ./location_of_script_within_repo/buildscript.sh

      - name: Perform CodeQL Analysis
-        uses: github/codeql-action/analyze@v4
+        uses: github/codeql-action/analyze@v2
@@ -39,7 +39,7 @@ jobs:

    steps:
      - name: checkout MFEM
-        uses: actions/checkout@v6
+        uses: actions/checkout@v4
        with:
          path: mfem

@@ -50,7 +50,7 @@ jobs:

      - name: Cache Hypre Install
        id: hypre-cache
-        uses: actions/cache@v5
+        uses: actions/cache@v4
        with:
          path: ${{ env.HYPRE_TOP_DIR }}
          key: ${{ runner.os }}-ompi-build-${{ env.HYPRE_TOP_DIR }}-v2.5
@@ -65,7 +65,7 @@ jobs:

      - name: Cache Metis Install
        id: metis-cache
-        uses: actions/cache@v5
+        uses: actions/cache@v4
        with:
          path: ${{ env.METIS_TOP_DIR }}
          key: ${{ runner.os }}-build-${{ env.METIS_TOP_DIR }}-v2.5
@@ -38,7 +38,7 @@ jobs:
       github.event.pull_request.head.repo.full_name != github.repository)
    steps:
      - name: checkout mfem
-        uses: actions/checkout@v6
+        uses: actions/checkout@v4

      - name: copyright check
        id: copyright
@@ -93,7 +93,7 @@ jobs:
       github.event.pull_request.head.repo.full_name != github.repository)
    steps:
      - name: checkout mfem
-        uses: actions/checkout@v6
+        uses: actions/checkout@v4

      - name: get astyle
        run: |
@@ -110,7 +110,7 @@ jobs:
       github.event.pull_request.head.repo.full_name != github.repository)
    steps:
      - name: checkout mfem
-        uses: actions/checkout@v6
+        uses: actions/checkout@v4

      - name: get doxygen and graphviz
        run: |
@@ -135,7 +135,7 @@ jobs:
    runs-on: ubuntu-latest
    steps:
      - name: checkout mfem
-        uses: actions/checkout@v6
+        uses: actions/checkout@v4
        with:
          fetch-depth: 0

@@ -17,11 +17,11 @@ jobs:
    runs-on: ubuntu-latest
    name: 2.19.0
    steps:
-      - uses: actions/checkout@v6
+      - uses: actions/checkout@v4
      - uses: ./.github/actions/sanitize/config
      - name: Cache
        id: cache
-        uses: actions/cache@v5
+        uses: actions/cache@v4
        with:
          path: ${{env.HYPRE_DIR}}
          key: ${{runner.os}}-ompi-build-${{env.HYPRE_DIR}}-int32-fp64-v2.5
@@ -27,13 +27,13 @@ jobs:
            llvm_use_sanitizer: "Undefined"
    name: ${{matrix.sanitizer}}
    steps:
-      - uses: actions/checkout@v6
+      - uses: actions/checkout@v4
      - uses: ./.github/actions/sanitize/config
        with:
          NO_FLAGS: true
      - name: Cache
        id: cache
-        uses: actions/cache@v5
+        uses: actions/cache@v4
        with:
          path: ${{env.LLVM_DIR}}
          key: build-libcxx-${{env.LLVM_VER}}-${{matrix.sanitizer}}
@@ -17,11 +17,11 @@ jobs:
    runs-on: ubuntu-latest
    name: lsan.supp
    steps:
-      - uses: actions/checkout@v6
+      - uses: actions/checkout@v4
      - uses: ./.github/actions/sanitize/config
      - name: Cache
        id: cache
-        uses: actions/cache@v5
+        uses: actions/cache@v4
        with:
          path: ${{env.LSAN_DIR}}
          key: build-lsan-suppression-file
@@ -17,11 +17,11 @@ jobs:
    runs-on: ubuntu-latest
    name: 4.0.3
    steps:
-      - uses: actions/checkout@v6
+      - uses: actions/checkout@v4
      - uses: ./.github/actions/sanitize/config
      - name: Cache
        id: cache
-        uses: actions/cache@v5
+        uses: actions/cache@v4
        with:
          path: ${{env.METIS_DIR}}
          key: ${{runner.os}}-build-${{env.METIS_DIR}}-v2.5
@@ -28,7 +28,7 @@ jobs:
  build:
    runs-on: ubuntu-latest
    steps:
-      - uses: actions/checkout@v6
+      - uses: actions/checkout@v4
      - uses: ./.github/actions/sanitize/mfem
        with:
          par: ${{inputs.par}}
@@ -40,7 +40,7 @@ jobs:
    env:
      ex: ${{inputs.par && 'ex1p' || 'ex1'}}
    steps:
-      - uses: actions/checkout@v6
+      - uses: actions/checkout@v4
      - uses: ./.github/actions/sanitize/restore
        id: restore
        with:
@@ -58,7 +58,7 @@ jobs:
    env:
      exclude: ${{inputs.par && '-E "_ser"' || ''}}
    steps:
-      - uses: actions/checkout@v6
+      - uses: actions/checkout@v4
      - uses: ./.github/actions/sanitize/restore
        id: restore
        with:
@@ -82,7 +82,7 @@ jobs:
    env:
      exclude: ${{inputs.par && '-E "_ser"' || ''}}
    steps:
-      - uses: actions/checkout@v6
+      - uses: actions/checkout@v4
      - uses: ./.github/actions/sanitize/restore
        id: restore
        with:
@@ -107,7 +107,7 @@ jobs:
      run: ${{inputs.par && '-R "_cpu_np"' || ''}}
      exclude: ${{inputs.par && '"unit_tests|debug"' || '"^unit_tests$|debug"'}}
    steps:
-      - uses: actions/checkout@v6
+      - uses: actions/checkout@v4
      - uses: ./.github/actions/sanitize/restore
        id: restore
        with:
@@ -131,7 +131,7 @@ jobs:
    env:
      unit_tests: ${{inputs.par && 'punit_tests' || 'unit_tests'}}
    steps:
-      - uses: actions/checkout@v6
+      - uses: actions/checkout@v4
      - uses: ./.github/actions/sanitize/restore
        id: restore
        with:
@@ -165,7 +165,7 @@ jobs:
      unit_tests: ${{inputs.par && 'punit_tests' || 'unit_tests'}}
      np: ${{inputs.par && '_np=2' || ''}}
    steps:
-      - uses: actions/checkout@v6
+      - uses: actions/checkout@v4
      - uses: ./.github/actions/sanitize/restore
        id: restore
        with:
@@ -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
@@ -8,22 +8,6 @@
                               https://mfem.org


-Version 4.10 (development)
-==========================
-
-Discretization improvements
---------------------------
- Replaced legacy simplex quadrature rules with symmetric positive-weight
-  rules for triangles (orders 0-25) and tetrahedra (orders 0-20). These
-  rules guarantee all-positive weights and interior quadrature points,
-  improving numerical stability. Higher orders fall back to Grundmann-Moller.
-  Triangle rules: Witherden & Vincent, Comput. Math. Appl. 69(10):1232-1241,
-  2015.
-  Tet rules (d=1-13): Witherden & Vincent (ibid).
-  Tet rules (d=14-20): Chuluunbaatar et al., Comput. Math. Appl. 124:89-97,
-  2022.
-
-
 Version 4.9.1 (development)
 ===========================

@@ -109,10 +109,6 @@ if (MFEM_USE_RAJA)
   find_dependency(RAJA)
 endif()

-if (MFEM_USE_UMPIRE)
-   find_dependency(umpire)
-endif()
-
 if (NOT TARGET mfem)
  include(${CMAKE_CURRENT_LIST_DIR}/MFEMTargets.cmake)
 endif (NOT TARGET mfem)
@@ -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")
@@ -97,9 +97,6 @@
 // Enable MFEM functionality based on the SuiteSparse library.
 // #define MFEM_USE_SUITESPARSE

-// Enable MFEM functionality based on the ARPACK library.
-// #define MFEM_USE_ARPACK
-
 // Enable MFEM functionality based on the SuperLU_DIST library.
 // #define MFEM_USE_SUPERLU
 // #define MFEM_USE_SUPERLU5
@@ -32,7 +32,6 @@ MFEM_USE_MEMALLOC      = @MFEM_USE_MEMALLOC@
 MFEM_TIMER_TYPE        = @MFEM_TIMER_TYPE@
 MFEM_USE_SUNDIALS      = @MFEM_USE_SUNDIALS@
 MFEM_USE_SUITESPARSE   = @MFEM_USE_SUITESPARSE@
-MFEM_USE_ARPACK        = @MFEM_USE_ARPACK@
 MFEM_USE_SUPERLU       = @MFEM_USE_SUPERLU@
 MFEM_USE_SUPERLU5      = @MFEM_USE_SUPERLU5@
 MFEM_USE_MUMPS         = @MFEM_USE_MUMPS@
@@ -178,7 +178,6 @@ MFEM_USE_ALGOIM        = NO
 MFEM_USE_UMPIRE        = NO
 MFEM_USE_SIMD          = NO
 MFEM_USE_ADIOS2        = NO
-MFEM_USE_ARPACK        = NO
 MFEM_USE_MKL_CPARDISO  = NO
 MFEM_USE_MKL_PARDISO   = NO
 MFEM_USE_MOONOLITH     = NO
@@ -428,14 +427,6 @@ NETCDF_LIB = $(XLINKER)-rpath,$(NETCDF_DIR)/lib -L$(NETCDF_DIR)/lib\
 $(XLINKER)-rpath,$(HDF5_DIR)/lib -L$(HDF5_DIR)/lib\
 -lnetcdf -lhdf5_hl -lhdf5 $(ZLIB_LIB)

-# ARPACK library configuration
-ARPACK_DIR = @MFEM_DIR@/../ARPACK
-ifeq ($(MFEM_USE_MPI),YES)
-ARPACK_LIB = -L$(ARPACK_DIR) -lparpack -larpack
-else
-ARPACK_LIB = -L$(ARPACK_DIR) -larpack
-endif
-
 # PETSc library configuration (version greater or equal to 3.8 or the dev branch)
 PETSC_ARCH := arch-linux2-c-debug
 PETSC_DIR  := $(MFEM_DIR)/../petsc/$(PETSC_ARCH)
@@ -49,13 +49,6 @@ list(APPEND ALL_EXE_SRCS
  ex41.cpp
  )

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

@@ -99,10 +97,6 @@ int main(int argc, char *argv[])
   args.AddOption(&sp_solver, "-sp", "--strumpack", "-no-sp",
                  "--no-strumpack", "Use the STRUMPACK Solver.");
 #endif
-#ifdef MFEM_USE_ARPACK
-   args.AddOption(&arp_solver, "-arp", "--arpack", "-no-arp",
-                  "--no-arpack", "Use the Parallel ARPACK Solver.");
-#endif
 #ifdef MFEM_USE_MKL_CPARDISO
   args.AddOption(&cpardiso_solver, "-cpardiso", "--cpardiso", "-no-cpardiso",
                  "--no-cpardiso", "Use the MKL CPardiso Solver.");
@@ -119,11 +113,6 @@ int main(int argc, char *argv[])
              << "         Defaulting to SuperLU." << endl;
      sp_solver = false;
   }
-   if (arp_solver)
-   {
-      lob_solver = false;
-   }
-
   // The command line options are also passed to the STRUMPACK
   // solver. So do not exit if some options are not recognized.
   if (!sp_solver)
@@ -254,119 +243,70 @@ int main(int argc, char *argv[])
   // 8. Define and configure the LOBPCG eigensolver and the BoomerAMG
   //    preconditioner for A to be used within the solver. Set the matrices
   //    which define the generalized eigenproblem A x = lambda M x.
-   Solver * solver = NULL;
   Solver * precond = NULL;
   if (!slu_solver && !sp_solver && !cpardiso_solver)
   {
      HypreBoomerAMG * amg = new HypreBoomerAMG(*A);
      amg->SetPrintLevel(0);
      precond = amg;
-
-      if (arp_solver)
-      {
-         HyprePCG * pcg = new HyprePCG(*A);
-         pcg->SetTol(1e-12);
-         pcg->SetPreconditioner(*amg);
-         solver = pcg;
-      }
   }
-#ifdef MFEM_USE_SUPERLU
-   else if (slu_solver)
+   else
   {
-      SuperLUSolver * superlu = new SuperLUSolver(MPI_COMM_WORLD);
-      superlu->SetPrintStatistics(false);
-      superlu->SetSymmetricPattern(true);
-      superlu->SetColumnPermutation(superlu::PARMETIS);
-      superlu->SetOperator(*Arow);
-
-      if (arp_solver)
-      {
-         solver = superlu;
-      }
-      else
+#ifdef MFEM_USE_SUPERLU
+      if (slu_solver)
      {
+         SuperLUSolver * superlu = new SuperLUSolver(MPI_COMM_WORLD);
+         superlu->SetPrintStatistics(false);
+         superlu->SetSymmetricPattern(true);
+         superlu->SetColumnPermutation(superlu::PARMETIS);
+         superlu->SetOperator(*Arow);
         precond = superlu;
      }
-   }
 #endif
 #ifdef MFEM_USE_STRUMPACK
-   else if (sp_solver)
-   {
-      STRUMPACKSolver * strumpack = new STRUMPACKSolver(argc, argv,
-                                                        MPI_COMM_WORLD);
-      strumpack->SetPrintFactorStatistics(true);
-      strumpack->SetPrintSolveStatistics(false);
-      strumpack->SetKrylovSolver(strumpack::KrylovSolver::DIRECT);
-      strumpack->SetReorderingStrategy(strumpack::ReorderingStrategy::METIS);
-      strumpack->SetMatching(strumpack::MatchingJob::NONE);
-      strumpack->SetCompression(strumpack::CompressionType::NONE);
-      strumpack->SetOperator(*Arow);
-      strumpack->SetFromCommandLine();
-      if (arp_solver)
-      {
-         solver = strumpack;
-      }
-      else
+      if (sp_solver)
      {
+         STRUMPACKSolver * strumpack = new STRUMPACKSolver(MPI_COMM_WORLD, argc, argv);
+         strumpack->SetPrintFactorStatistics(true);
+         strumpack->SetPrintSolveStatistics(false);
+         strumpack->SetKrylovSolver(strumpack::KrylovSolver::DIRECT);
+         strumpack->SetReorderingStrategy(strumpack::ReorderingStrategy::METIS);
+         strumpack->SetMatching(strumpack::MatchingJob::NONE);
+         strumpack->SetCompression(strumpack::CompressionType::NONE);
+         strumpack->SetOperator(*Arow);
+         strumpack->SetFromCommandLine();
         precond = strumpack;
      }
-   }
 #endif
 #ifdef MFEM_USE_MKL_CPARDISO
-   else if (cpardiso_solver)
-   {
-      auto cpardiso = new CPardisoSolver(A->GetComm());
-      cpardiso->SetMatrixType(CPardisoSolver::MatType::REAL_STRUCTURE_SYMMETRIC);
-      cpardiso->SetPrintLevel(1);
-      cpardiso->SetOperator(*A);
-      if (arp_solver)
-      {
-         solver = cpardiso;
-      }
-      else
+      if (cpardiso_solver)
      {
+         auto cpardiso = new CPardisoSolver(A->GetComm());
+         cpardiso->SetMatrixType(CPardisoSolver::MatType::REAL_STRUCTURE_SYMMETRIC);
+         cpardiso->SetPrintLevel(1);
+         cpardiso->SetOperator(*A);
         precond = cpardiso;
      }
-   }
 #endif
-
-   SymGenEigensolver * eig_solver = NULL;
-
-   if (lob_solver)
-   {
-      HypreLOBPCG * lobpcg = new HypreLOBPCG(MPI_COMM_WORLD);
-      lobpcg->SetNumModes(nev);
-      lobpcg->SetRandomSeed(seed);
-      lobpcg->SetPreconditioner(*precond);
-      lobpcg->SetMaxIter(200);
-      lobpcg->SetTol(1e-8);
-      lobpcg->SetPrecondUsageMode(1);
-      lobpcg->SetPrintLevel(1);
-
-      eig_solver = lobpcg;
   }
-#ifdef MFEM_USE_ARPACK
-   else if (arp_solver)
-   {
-      ArPackPSAUPD * arpack = new ArPackPSAUPD(MPI_COMM_WORLD);
-      arpack->SetNumModes(nev);
-      arpack->SetMaxIter(400);
-      arpack->SetTol(1e-8);
-      arpack->SetMode(3);
-      arpack->SetPrintLevel(2);
-      arpack->SetSolver(*solver);

-      eig_solver = arpack;
-   }
-#endif
-   eig_solver->SetOperators(*A, *M);
+   HypreLOBPCG * lobpcg = new HypreLOBPCG(MPI_COMM_WORLD);
+   lobpcg->SetNumModes(nev);
+   lobpcg->SetRandomSeed(seed);
+   lobpcg->SetPreconditioner(*precond);
+   lobpcg->SetMaxIter(200);
+   lobpcg->SetTol(1e-8);
+   lobpcg->SetPrecondUsageMode(1);
+   lobpcg->SetPrintLevel(1);
+   lobpcg->SetMassMatrix(*M);
+   lobpcg->SetOperator(*A);

   // 9. Compute the eigenmodes and extract the array of eigenvalues. Define a
   //    parallel grid function to represent each of the eigenmodes returned by
   //    the solver.
   Array<real_t> eigenvalues;
-   eig_solver->Solve();
-   eig_solver->GetEigenvalues(eigenvalues);
+   lobpcg->Solve();
+   lobpcg->GetEigenvalues(eigenvalues);
   ParGridFunction x(fespace);

   // 10. Save the refined mesh and the modes in parallel. This output can be
@@ -381,8 +321,8 @@ int main(int argc, char *argv[])

      for (int i=0; i<nev; i++)
      {
-         // convert eigenvector from Vector to ParGridFunction
-         x.Distribute(eig_solver->GetEigenvector(i));
+         // convert eigenvector from HypreParVector to ParGridFunction
+         x = lobpcg->GetEigenvector(i);

         mode_name << "mode_" << setfill('0') << setw(2) << i << "."
                   << setfill('0') << setw(6) << myid;
@@ -410,8 +350,8 @@ int main(int argc, char *argv[])
                 << ", Lambda = " << eigenvalues[i] << endl;
         }

-         // convert eigenvector from Vector to ParGridFunction
-         x.Distribute(eig_solver->GetEigenvector(i));
+         // convert eigenvector from HypreParVector to ParGridFunction
+         x = lobpcg->GetEigenvector(i);

         mode_sock << "parallel " << num_procs << " " << myid << "\n"
                   << "solution\n" << *pmesh << x << flush
@@ -435,8 +375,7 @@ int main(int argc, char *argv[])
   }

   // 12. Free the used memory.
-   delete eig_solver;
-   delete solver;
+   delete lobpcg;
   delete precond;
   delete M;
   delete A;
@@ -1,381 +0,0 @@
-//                       MFEM Example 11 - Parallel Version
-//
-// Compile with: make ex11p
-//
-// Sample runs:  mpirun -np 4 ex11p -m ../data/square-disc.mesh
-//               mpirun -np 4 ex11p -m ../data/star.mesh
-//               mpirun -np 4 ex11p -m ../data/escher.mesh
-//               mpirun -np 4 ex11p -m ../data/fichera.mesh
-//               mpirun -np 4 ex11p -m ../data/square-disc-p2.vtk -o 2
-//               mpirun -np 4 ex11p -m ../data/square-disc-p3.mesh -o 3
-//               mpirun -np 4 ex11p -m ../data/square-disc-nurbs.mesh -o -1
-//               mpirun -np 4 ex11p -m ../data/disc-nurbs.mesh -o -1 -n 20
-//               mpirun -np 4 ex11p -m ../data/pipe-nurbs.mesh -o -1
-//               mpirun -np 4 ex11p -m ../data/ball-nurbs.mesh -o 2
-//               mpirun -np 4 ex11p -m ../data/star-surf.mesh
-//               mpirun -np 4 ex11p -m ../data/square-disc-surf.mesh
-//               mpirun -np 4 ex11p -m ../data/inline-segment.mesh
-//               mpirun -np 4 ex11p -m ../data/amr-quad.mesh
-//               mpirun -np 4 ex11p -m ../data/amr-hex.mesh
-//               mpirun -np 4 ex11p -m ../data/mobius-strip.mesh -n 8
-//               mpirun -np 4 ex11p -m ../data/klein-bottle.mesh -n 10
-//
-// Description:  This example code demonstrates the use of MFEM to solve the
-//               eigenvalue problem -Delta u = lambda u with homogeneous
-//               Dirichlet boundary conditions.
-//
-//               We compute a number of the lowest eigenmodes by discretizing
-//               the Laplacian and Mass operators using a FE space of the
-//               specified order, or an isoparametric/isogeometric space if
-//               order < 1 (quadratic for quadratic curvilinear mesh, NURBS for
-//               NURBS mesh, etc.)
-//
-//               The example highlights the use of the LOBPCG and ARPACK
-//               eigenvalue solvers together with the BoomerAMG preconditioner
-//               in HYPRE, as well as optionally the SuperLU parallel direct
-//               solver. Reusing a single GLVis visualization window for
-//               multiple eigenfunctions is also illustrated.
-//
-//               We recommend viewing Example 1 before viewing this example.
-
-#include "mfem.hpp"
-#include <fstream>
-#include <iostream>
-
-using namespace std;
-using namespace mfem;
-
-int main(int argc, char *argv[])
-{
-   // 1. Initialize MPI.
-   int num_procs, myid;
-   MPI_Init(&argc, &argv);
-   MPI_Comm_size(MPI_COMM_WORLD, &num_procs);
-   MPI_Comm_rank(MPI_COMM_WORLD, &myid);
-
-   // 2. Parse command-line options.
-   const char *mesh_file = "../data/star.mesh";
-   int ser_ref_levels = 2;
-   int par_ref_levels = 1;
-   int order = 1;
-   int nev = 5;
-   bool slu_solver = false;
-   bool use_arpack = false;
-   bool visualization = 1;
-
-   OptionsParser args(argc, argv);
-   args.AddOption(&mesh_file, "-m", "--mesh",
-                  "Mesh file to use.");
-   args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
-                  "Number of times to refine the mesh uniformly in serial.");
-   args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
-                  "Number of times to refine the mesh uniformly in parallel.");
-   args.AddOption(&order, "-o", "--order",
-                  "Finite element order (polynomial degree) or -1 for"
-                  " isoparametric space.");
-   args.AddOption(&nev, "-n", "--num-eigs",
-                  "Number of desired eigenmodes.");
-#ifdef MFEM_USE_SUPERLU
-   args.AddOption(&slu_solver, "-slu", "--superlu", "-no-slu",
-                  "--no-superlu", "Use the SuperLU Solver.");
-#endif
-#ifdef MFEM_USE_ARPACK
-   args.AddOption(&use_arpack, "-arpack", "--use-arpack", "-no-arpack",
-                  "--no-arpack",
-                  "Enable or disable the use of ARPACK.");
-#endif
-   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
-                  "--no-visualization",
-                  "Enable or disable GLVis visualization.");
-   args.Parse();
-   if (!args.Good())
-   {
-      if (myid == 0)
-      {
-         args.PrintUsage(cout);
-      }
-      MPI_Finalize();
-      return 1;
-   }
-   if (myid == 0)
-   {
-      args.PrintOptions(cout);
-   }
-
-   // 3. Read the (serial) mesh from the given mesh file on all processors. We
-   //    can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
-   //    and volume meshes with the same code.
-   Mesh *mesh;
-   ifstream imesh(mesh_file);
-   if (!imesh)
-   {
-      if (myid == 0)
-      {
-         cerr << "\nCan not open mesh file: " << mesh_file << '\n' << endl;
-      }
-      MPI_Finalize();
-      return 2;
-   }
-   mesh = new Mesh(imesh, 1, 1);
-   imesh.close();
-   int dim = mesh->Dimension();
-
-   // 4. Refine the serial mesh on all processors to increase the resolution. In
-   //    this example we do 'ref_levels' of uniform refinement (2 by default, or
-   //    specified on the command line with -rs).
-   for (int lev = 0; lev < ser_ref_levels; lev++)
-   {
-      mesh->UniformRefinement();
-   }
-
-   // 5. Define a parallel mesh by a partitioning of the serial mesh. Refine
-   //    this mesh further in parallel to increase the resolution (1 time by
-   //    default, or specified on the command line with -rp). Once the parallel
-   //    mesh is defined, the serial mesh can be deleted.
-   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
-   delete mesh;
-   for (int lev = 0; lev < par_ref_levels; lev++)
-   {
-      pmesh->UniformRefinement();
-   }
-
-   // 6. Define a parallel finite element space on the parallel mesh. Here we
-   //    use continuous Lagrange finite elements of the specified order. If
-   //    order < 1, we instead use an isoparametric/isogeometric space.
-   FiniteElementCollection *fec;
-   if (order > 0)
-   {
-      fec = new H1_FECollection(order, dim);
-   }
-   else if (pmesh->GetNodes())
-   {
-      fec = pmesh->GetNodes()->OwnFEC();
-   }
-   else
-   {
-      fec = new H1_FECollection(order = 1, dim);
-   }
-   ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
-   HYPRE_Int size = fespace->GlobalTrueVSize();
-   if (myid == 0)
-   {
-      cout << "Number of unknowns: " << size << endl;
-   }
-
-   // 7. Set up the parallel bilinear forms a(.,.) and m(.,.) on the finite
-   //    element space. The first corresponds to the Laplacian operator -Delta,
-   //    while the second is a simple mass matrix needed on the right hand side
-   //    of the generalized eigenvalue problem below. The boundary conditions
-   //    are implemented by elimination with special values on the diagonal to
-   //    shift the Dirichlet eigenvalues out of the computational range. After
-   //    serial and parallel assembly we extract the corresponding parallel
-   //    matrices A and M.
-   ConstantCoefficient one(1.0);
-   Array<int> ess_bdr;
-   if (pmesh->bdr_attributes.Size())
-   {
-      ess_bdr.SetSize(pmesh->bdr_attributes.Max());
-      ess_bdr = 1;
-   }
-
-   ParBilinearForm *a = new ParBilinearForm(fespace);
-   a->AddDomainIntegrator(new DiffusionIntegrator(one));
-   if (pmesh->bdr_attributes.Size() == 0)
-   {
-      // Add a mass term if the mesh has no boundary, e.g. periodic mesh or
-      // closed surface.
-      a->AddDomainIntegrator(new MassIntegrator(one));
-   }
-   a->Assemble();
-   a->EliminateEssentialBCDiag(ess_bdr, 1.0);
-   a->Finalize();
-
-   ParBilinearForm *m = new ParBilinearForm(fespace);
-   m->AddDomainIntegrator(new MassIntegrator(one));
-   m->Assemble();
-   // shift the eigenvalue corresponding to eliminated dofs to a large value
-   m->EliminateEssentialBCDiag(ess_bdr, numeric_limits<double>::min());
-   m->Finalize();
-
-   HypreParMatrix *A = a->ParallelAssemble();
-   HypreParMatrix *M = m->ParallelAssemble();
-
-#ifdef MFEM_USE_SUPERLU
-   Operator * Arow = NULL;
-   if (slu_solver)
-   {
-      Arow = new SuperLURowLocMatrix(*A);
-   }
-#endif
-
-   delete a;
-   delete m;
-
-   // 8. Define and configure the LOBPCG eigensolver and the BoomerAMG
-   //    preconditioner for A to be used within the solver. Set the matrices
-   //    which define the generalized eigenproblem A x = lambda M x.
-   Eigensolver * esolver = NULL;
-   Solver      *  solver = NULL;
-   Solver      * precond = NULL;
-
-   if (!slu_solver)
-   {
-      HypreBoomerAMG * amg = new HypreBoomerAMG(*A);
-      amg->SetPrintLevel(0);
-      precond = amg;
-
-#ifdef MFEM_USE_ARPACK
-      if ( use_arpack )
-      {
-         HyprePCG * pcg = new HyprePCG(*A);
-         pcg->SetTol(1e-12);
-         pcg->SetMaxIter(200);
-         pcg->SetPreconditioner(*amg);
-         pcg->SetPrintLevel(0);
-
-         solver = pcg;
-      }
-#endif
-   }
-#ifdef MFEM_USE_SUPERLU
-   else
-   {
-      SuperLUSolver * superlu = new SuperLUSolver(MPI_COMM_WORLD);
-      superlu->SetPrintStatistics(false);
-      superlu->SetSymmetricPattern(true);
-      superlu->SetColumnPermutation(superlu::PARMETIS);
-      superlu->SetOperator(*Arow);
-
-      solver  = use_arpack?superlu:NULL;
-      precond = use_arpack?NULL:superlu;
-   }
-#endif
-
-   if ( use_arpack )
-   {
-      ParArPackSym * arpack = new ParArPackSym(MPI_COMM_WORLD);
-      arpack->SetMode(3);
-      arpack->SetPrintLevel(2);
-      arpack->SetSolver(*solver);
-
-      esolver = arpack;
-   }
-   else
-   {
-      HypreLOBPCG * lobpcg = new HypreLOBPCG(MPI_COMM_WORLD);
-      lobpcg->SetPreconditioner(*precond);
-      lobpcg->SetPrecondUsageMode(1);
-      lobpcg->SetPrintLevel(1);
-
-      esolver = lobpcg;
-   }
-
-   esolver->SetNumModes(nev);
-   esolver->SetMaxIter(100);
-   esolver->SetTol(1e-8);
-
-   esolver->SetMassMatrix(*M);
-   esolver->SetOperator(*A);
-
-   // 9. Compute the eigenmodes and extract the array of eigenvalues. Define a
-   //    parallel grid function to represent each of the eigenmodes returned by
-   //    the solver.
-   Array<double> eigenvalues;
-   esolver->Solve();
-   esolver->GetEigenvalues(eigenvalues);
-
-   if ( myid == 0 && use_arpack )
-   {
-      cout << endl;
-      for (int i=0; i<eigenvalues.Size(); i++)
-      {
-         cout << "Eigenvalue lambda " << eigenvalues[i] << endl;
-      }
-      cout << endl;
-   }
-
-   ParGridFunction x(fespace);
-
-   // 10. Save the refined mesh and the modes in parallel. This output can be
-   //     viewed later using GLVis: "glvis -np <np> -m mesh -g mode".
-   {
-      ostringstream mesh_name, mode_name;
-      mesh_name << "mesh." << setfill('0') << setw(6) << myid;
-
-      ofstream mesh_ofs(mesh_name.str().c_str());
-      mesh_ofs.precision(8);
-      pmesh->Print(mesh_ofs);
-
-      for (int i=0; i<nev; i++)
-      {
-         // convert eigenvector from HypreParVector to ParGridFunction
-         x.Distribute(esolver->GetEigenvector(i));
-
-         mode_name << "mode_" << setfill('0') << setw(2) << i << "."
-                   << setfill('0') << setw(6) << myid;
-
-         ofstream mode_ofs(mode_name.str().c_str());
-         mode_ofs.precision(8);
-         x.Save(mode_ofs);
-         mode_name.str("");
-      }
-   }
-
-   // 11. Send the solution by socket to a GLVis server.
-   if (visualization)
-   {
-      char vishost[] = "localhost";
-      int  visport   = 19916;
-      socketstream mode_sock(vishost, visport);
-      mode_sock.precision(8);
-
-      for (int i=0; i<nev; i++)
-      {
-         if ( myid == 0 )
-         {
-            cout << "Eigenmode " << i+1 << '/' << nev
-                 << ", Lambda = " << eigenvalues[i] << endl;
-         }
-
-         // convert eigenvector from HypreParVector to ParGridFunction
-         x.Distribute(esolver->GetEigenvector(i));
-
-         mode_sock << "parallel " << num_procs << " " << myid << "\n"
-                   << "solution\n" << *pmesh << x << flush
-                   << "window_title 'Eigenmode " << i+1 << '/' << nev
-                   << ", Lambda = " << eigenvalues[i] << "'" << endl;
-
-         char c;
-         if (myid == 0)
-         {
-            cout << "press (q)uit or (c)ontinue --> " << flush;
-            cin >> c;
-         }
-         MPI_Bcast(&c, 1, MPI_CHAR, 0, MPI_COMM_WORLD);
-
-         if (c != 'c')
-         {
-            break;
-         }
-      }
-      mode_sock.close();
-   }
-
-   // 12. Free the used memory.
-   delete esolver;
-   delete solver;
-   delete precond;
-   delete M;
-   delete A;
-
-   delete fespace;
-   if (order > 0)
-   {
-      delete fec;
-   }
-   delete pmesh;
-
-   MPI_Finalize();
-
-   return 0;
-}
@@ -5,9 +5,9 @@
 // Sample runs:
 //    mpirun -np 4 ex12p -m ../data/beam-tri.mesh
 //    mpirun -np 4 ex12p -m ../data/beam-quad.mesh
-//    mpirun -np 4 ex12p -m ../data/beam-tet.mesh -s 464 -n 10 -o 2 -elast
+//    mpirun -np 4 ex12p -m ../data/beam-tet.mesh -s 462 -n 10 -o 2 -elast
 //    mpirun -np 4 ex12p -m ../data/beam-hex.mesh -s 3878
-//    mpirun -np 4 ex12p -m ../data/beam-wedge.mesh -s 82
+//    mpirun -np 4 ex12p -m ../data/beam-wedge.mesh -s 81
 //    mpirun -np 4 ex12p -m ../data/beam-tri.mesh -s 3877 -o 2 -sys
 //    mpirun -np 4 ex12p -m ../data/beam-quad.mesh -s 4544 -n 6 -o 3 -elast
 //    mpirun -np 4 ex12p -m ../data/beam-quad-nurbs.mesh
@@ -276,8 +276,8 @@ int main(int argc, char *argv[])

      for (int i=0; i<nev; i++)
      {
-         // convert eigenvector from Vector to ParGridFunction
-         x.Distribute(lobpcg->GetEigenvector(i));
+         // convert eigenvector from HypreParVector to ParGridFunction
+         x = lobpcg->GetEigenvector(i);

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

-         // convert eigenvector from Vector to ParGridFunction
-         x.Distribute(lobpcg->GetEigenvector(i));
+         // convert eigenvector from HypreParVector to ParGridFunction
+         x = lobpcg->GetEigenvector(i);

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

      for (int i=0; i<nev; i++)
      {
-         // convert eigenvector from Vector to ParGridFunction
-         x.Distribute(ame->GetEigenvector(i));
+         // convert eigenvector from HypreParVector to ParGridFunction
+         x = ame->GetEigenvector(i);

         mode_name << "mode_" << setfill('0') << setw(2) << i << "."
                   << setfill('0') << setw(6) << myid;
@@ -244,8 +244,8 @@ int main(int argc, char *argv[])
                 << ", Lambda = " << eigenvalues[i] << endl;
         }

-         // convert eigenvector from Vector to ParGridFunction
-         x.Distribute(ame->GetEigenvector(i));
+         // convert eigenvector from HypreParVector to ParGridFunction
+         x = ame->GetEigenvector(i);

         mode_sock << "parallel " << num_procs << " " << myid << "\n"
                   << "solution\n" << *pmesh << x << flush
@@ -302,21 +302,15 @@ int main(int argc, char *argv[])
      char vishost[] = "localhost";
      int  visport   = 19916;
      socketstream sol_sock_r(vishost, visport);
+      socketstream sol_sock_i(vishost, visport);
      sol_sock_r << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
      sol_sock_r.precision(8);
+      sol_sock_i.precision(8);
      sol_sock_r << "solution\n" << *pmesh << u_exact->real()
                 << "window_title 'Exact: Real Part'" << flush;
-      // Make sure all ranks have sent their real solution before initiating
-      // another set of GLVis connections (one from each rank):
-      MPI_Barrier(pmesh->GetComm());
-      socketstream sol_sock_i(vishost, visport);
-      sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
-      sol_sock_i.precision(8);
      sol_sock_i << "solution\n" << *pmesh << u_exact->imag()
                 << "window_title 'Exact: Imaginary Part'" << flush;
-      // Make sure all ranks have sent their imaginary solution before initiating
-      // another set of GLVis connections (one from each rank):
-      MPI_Barrier(pmesh->GetComm());
   }

   // 11. Set up the parallel sesquilinear form a(.,.) on the finite element
@@ -540,21 +534,15 @@ int main(int argc, char *argv[])
      char vishost[] = "localhost";
      int  visport   = 19916;
      socketstream sol_sock_r(vishost, visport);
+      socketstream sol_sock_i(vishost, visport);
      sol_sock_r << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
      sol_sock_r.precision(8);
+      sol_sock_i.precision(8);
      sol_sock_r << "solution\n" << *pmesh << u.real()
                 << "window_title 'Solution: Real Part'" << flush;
-      // Make sure all ranks have sent their real solution before initiating
-      // another set of GLVis connections (one from each rank):
-      MPI_Barrier(pmesh->GetComm());
-      socketstream sol_sock_i(vishost, visport);
-      sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
-      sol_sock_i.precision(8);
      sol_sock_i << "solution\n" << *pmesh << u.imag()
                 << "window_title 'Solution: Imaginary Part'" << flush;
-      // Make sure all ranks have sent their imaginary solution before initiating
-      // another set of GLVis connections (one from each rank):
-      MPI_Barrier(pmesh->GetComm());
   }
   if (visualization && exact_sol)
   {
@@ -563,21 +551,15 @@ int main(int argc, char *argv[])
      char vishost[] = "localhost";
      int  visport   = 19916;
      socketstream sol_sock_r(vishost, visport);
+      socketstream sol_sock_i(vishost, visport);
      sol_sock_r << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
      sol_sock_r.precision(8);
+      sol_sock_i.precision(8);
      sol_sock_r << "solution\n" << *pmesh << u_exact->real()
                 << "window_title 'Error: Real Part'" << flush;
-      // Make sure all ranks have sent their real solution before initiating
-      // another set of GLVis connections (one from each rank):
-      MPI_Barrier(pmesh->GetComm());
-      socketstream sol_sock_i(vishost, visport);
-      sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
-      sol_sock_i.precision(8);
      sol_sock_i << "solution\n" << *pmesh << u_exact->imag()
                 << "window_title 'Error: Imaginary Part'" << flush;
-      // Make sure all ranks have sent their imaginary solution before initiating
-      // another set of GLVis connections (one from each rank):
-      MPI_Barrier(pmesh->GetComm());
   }
   if (visualization)
   {
@@ -228,7 +228,7 @@ int main(int argc, char *argv[])
      for (int i=0; i<nev; i++)
      {
         // convert eigenvector from HypreParVector to ParGridFunction
-         x.Distribute(ame->GetEigenvector(i));
+         x = ame->GetEigenvector(i);
         curl.Mult(x, dx);

         mode_name << "mode_" << setfill('0') << setw(2) << i << "."
@@ -295,7 +295,7 @@ int main(int argc, char *argv[])
            }

            // convert eigenvector from HypreParVector to ParGridFunction
-            x.Distribute(ame->GetEigenvector(i));
+            x = ame->GetEigenvector(i);
            curl.Mult(x, dx);

            {
@@ -469,7 +469,7 @@ int main(int argc, char *argv[])
            }

            // convert eigenvector from HypreParVector to ParGridFunction
-            x.Distribute(ame->GetEigenvector(i));
+            x = ame->GetEigenvector(i);
            curl.Mult(x, dx);

            {
@@ -599,7 +599,7 @@ int main(int argc, char *argv[])
            }

            // convert eigenvector from HypreParVector to ParGridFunction
-            x.Distribute(ame->GetEigenvector(i));
+            x = ame->GetEigenvector(i);
            curl.Mult(x, dx);

            mode_sock << "parallel " << num_procs << " " << myid << "\n"
@@ -658,7 +658,7 @@ void ScalarWaveGuide(int mode, ParGridFunction &x)
   lobpcg.SetOperator(*A);
   lobpcg.Solve();

-   x.Distribute(lobpcg.GetEigenvector(mode));
+   x = lobpcg.GetEigenvector(mode);

   delete A;
   delete M;
@@ -714,7 +714,7 @@ void VectorWaveGuide(int mode, ParGridFunction &x)
   ame.SetOperator(*A);
   ame.Solve();

-   x.Distribute(ame.GetEigenvector(mode));
+   x = ame.GetEigenvector(mode);

   delete A;
   delete M;
@@ -780,7 +780,7 @@ void PseudoScalarWaveGuide(int mode, ParGridFunction &x_l2)
   lobpcg.SetOperator(*A);
   lobpcg.Solve();

-   x.Distribute(lobpcg.GetEigenvector(mode));
+   x = lobpcg.GetEigenvector(mode);

   x_l2.ProjectCoefficient(xCoef);

@@ -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);
@@ -31,9 +31,6 @@ SEQ_DEVICE_EXAMPLES = ex1 ex3 ex4 ex5 ex6 ex9 ex14 ex22 ex24 ex25 ex26 ex34
 PAR_DEVICE_EXAMPLES = ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex9p ex13p ex14p \
   ex22p ex24p ex25p ex26p ex34p ex35p

-ifeq ($(MFEM_USE_ARPACK),YES)
-   SEQ_EXAMPLES += ex11 ex13
-endif
 ifeq ($(MFEM_USE_LAPACK),YES)
   SEQ_EXAMPLES += ex38
 endif
@@ -160,8 +157,6 @@ ex37-test-seq: ex37
 	@$(call mfem-test,$<,, Serial example,-mi 3)
 ex37p-test-par: ex37p
 	@$(call mfem-test,$<, $(RUN_MPI), Parallel example,-mi 3)
-ex39-test-seq: ex39
-	@$(call mfem-test,$<,, Serial example,-m ../data/compass.mesh)
 ex41-test-seq: ex41
 	@$(call mfem-test,$<,, Serial example,-tf 1.0)
 ex41p-test-par: ex41p
@@ -213,6 +213,7 @@ set(HDRS
  dfem/qfunction_apply.hpp
  dfem/qfunction_transform.hpp
  dfem/tuple.hpp
+  dfem/univarsolvers.hpp
  dfem/util.hpp
  eltrans.hpp
  estimators.hpp
@@ -729,8 +729,7 @@ void BilinearForm::Assemble(int skip_zeros)
         tr = mesh -> GetBdrFaceTransformations (i);
         if (tr != NULL)
         {
-            mfem::DofTransformation doftrans;
-            fes -> GetElementVDofs (tr -> Elem1No, vdofs, doftrans);
+            fes -> GetElementVDofs (tr -> Elem1No, vdofs);
            fe1 = fes -> GetFE (tr -> Elem1No);
            // The fe2 object is really a dummy and not used on the boundaries,
            // but we can't dereference a NULL pointer, and we don't want to
@@ -744,7 +743,6 @@ void BilinearForm::Assemble(int skip_zeros)

               boundary_face_integs[k] -> AssembleFaceMatrix (*fe1, *fe2, *tr,
                                                              elemmat);
-               doftrans.TransformDual(elemmat);
               mat -> AddSubMatrix (vdofs, vdofs, elemmat, skip_zeros);
            }
         }
@@ -1725,7 +1723,6 @@ void MixedBilinearForm::Assemble(int skip_zeros)
         }
      }

-      DofTransformation dom_dof_trans, ran_dof_trans;
      for (int i = 0; i < trial_fes -> GetNBE(); i++)
      {
         const int bdr_attr = mesh->GetBdrAttribute(i);
@@ -1734,8 +1731,8 @@ void MixedBilinearForm::Assemble(int skip_zeros)
         ftr = mesh -> GetBdrFaceTransformations (i);
         if (ftr != NULL)
         {
-            trial_fes->GetElementVDofs(ftr->Elem1No, trial_vdofs, dom_dof_trans);
-            test_fes->GetElementVDofs(ftr->Elem1No, test_vdofs, ran_dof_trans);
+            trial_fes->GetElementVDofs(ftr->Elem1No, trial_vdofs);
+            test_fes->GetElementVDofs(ftr->Elem1No, test_vdofs);
            trial_fe1 = trial_fes->GetFE(ftr->Elem1No);
            test_fe1 = test_fes->GetFE(ftr->Elem1No);
            // The test_fe2 object is really a dummy and not used on the
@@ -1751,7 +1748,6 @@ void MixedBilinearForm::Assemble(int skip_zeros)
               boundary_face_integs[k]->AssembleFaceMatrix(*trial_fe1, *test_fe1, *trial_fe2,
                                                           *test_fe2,
                                                           *ftr, elemmat);
-               TransformDual(ran_dof_trans, dom_dof_trans, elemmat);
               mat->AddSubMatrix(test_vdofs, trial_vdofs, elemmat, skip_zeros);
            }
         }
@@ -2710,7 +2710,7 @@ public:


 /** Integrator for $(-Q u, \nabla v)$ for Nedelec ($u$) and $H^1$ ($v$) elements.
-    This is equivalent to a weak divergence of the $H(curl)$ basis functions. */
+    This is equivalent to a weak divergence of the $H(curl$ basis functions. */
 class VectorFEWeakDivergenceIntegrator: public BilinearFormIntegrator
 {
 protected:
@@ -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
@@ -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";
@@ -0,0 +1,263 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+/**
+ * @file univarsolvers.hpp
+ *
+ * @brief Solvers of functions of a single variable suitable for use in ∂FEM q-functions.
+ */
+
+#ifndef MFEM_UNIVARSOLVERS
+#define MFEM_UNIVARSOLVERS
+
+#include <cmath>
+#include <limits>
+
+#include "../../config/config.hpp"
+#include "../../general/enzyme.hpp"
+#include "../../general/error.hpp"
+
+#ifdef MFEM_USE_ENZYME
+
+// Currently needed to work around a bug in LLVM
+extern void __enzyme_double(void*, size_t);
+
+namespace mfem
+{
+
+namespace future
+{
+/// Representation of root search bounds
+struct Bounds
+{
+   real_t lower, upper;
+};
+
+
+/// Settings for univariate solver
+struct SolverSettings
+{
+   real_t residual_abs_tol; ///< Tolerance for convergence check on absolute value of residual
+   real_t residual_rel_tol; ///< Tolerance for convergence check on absolute value of current residual relative to absolute value of residual at initial guess
+   Bounds bounds; ///< Bounds on root
+   int max_iters = 50;
+};
+} // namespace future
+
+namespace internal
+{
+/// @cond DO_NOT_DOCUMENT
+
+using future::SolverSettings;
+
+// The noinline attribute is neccessary for Enzyme. If this function were to be
+//  inlined in the calling function, The custom derivative rules would not be
+// found (since the function they refer to would no longer exist).
+template <auto f, typename T>
+__attribute__((noinline))
+MFEM_HOST_DEVICE void SolveNewtonBisection_impl(const real_t* x0_ptr,
+                                                const T* p_ptr, const SolverSettings* settings_ptr, real_t* x_ptr)
+{
+   int max_iters = settings_ptr->max_iters;
+
+   const real_t& x0 = *x0_ptr;
+   const T& p = *p_ptr;
+   const SolverSettings& settings = *settings_ptr;
+   const real_t& left_bracket = settings.bounds.lower;
+   const real_t& right_bracket = settings.bounds.upper;
+   real_t& x = *x_ptr;
+   using std::abs;
+
+   auto fprime = [&p](real_t x)
+   {
+      real_t x_dot = 1.0;
+      return __enzyme_fwddiff<real_t>((void*)+f, enzyme_dup, x, x_dot, enzyme_const,
+                                      p);
+   };
+
+   real_t fl = f(left_bracket, p);
+   real_t fh = f(right_bracket, p);
+
+   // handle corner cases where one of the brackets is the root
+   if (abs(fl) < settings.residual_abs_tol)
+   {
+      x = left_bracket;
+      return;
+   }
+   else if (abs(fh) < settings.residual_abs_tol)
+   {
+      x = right_bracket;
+      return;
+   }
+
+   if (fl * fh > 0)
+   {
+      MFEM_WARNING("Root is not bracketed, solver may diverge.");
+   }
+
+   // clamp initial guess within root brackets
+   x = x0 > right_bracket? right_bracket : x0 < left_bracket? left_bracket : x0;
+
+   // Orient search so that f(xl) < 0
+   real_t xl = left_bracket;
+   real_t xh = right_bracket;
+   if (fl > 0.0)
+   {
+      xl = right_bracket;
+      xh = left_bracket;
+      real_t tmp = fl;
+      fl = fh;
+      fh = tmp;
+   }
+
+   real_t dx_old = abs(right_bracket - left_bracket);
+   real_t dx = dx_old;
+   x = x0;
+   real_t r = f(x, p);
+   real_t dr_dx = fprime(x);
+   real_t r0 = r;
+   for (int i = 0; i < max_iters; i++)
+   {
+      if ((((x - xh) * dr_dx - r)*((x - xl)*dr_dx - r) >= 0.0) ||
+          // Newton out of range
+          (std::abs(2.0*r) > std::abs(
+              dx_old*dr_dx)))              // Newton decreasing bracket slower than bisection
+      {
+         // Take bisection step
+         dx_old = dx;
+         dx = 0.5*(xh - xl);
+         real_t x_old = x;
+         x = xl + dx;
+         if (x == x_old) { return; }
+      }
+      else
+      {
+         // Take Newton step
+         dx_old = dx;
+         dx = -r/dr_dx;
+         real_t x_old = x;
+         x += dx;
+         if (x == x_old) { return; }
+      }
+
+      // update residual and jacobian
+      r = f(x, p);
+      dr_dx = fprime(x);
+
+      // Check convergence
+      if (abs(r) < settings.residual_rel_tol*r0 ||
+          abs(r) < settings.residual_abs_tol)
+      {
+         return;
+      }
+
+      // Update bracket
+      if (r < 0.0)
+      {
+         xl = x;
+         fl = r;
+      }
+      else
+      {
+         xh = x;
+         fh = r;
+      }
+   }
+   MFEM_ABORT("Univariate solve did not converge.");
+}
+
+template <auto f, typename T>
+void SolveNewtonBisection_impl_fwddiff(const real_t* x0,
+                                       const real_t* /* unused shadow */,
+                                       const T* p, const T* dp,
+                                       const SolverSettings* settings, const SolverSettings* /* unused shadow */,
+                                       real_t* x, real_t* dx)
+{
+   SolveNewtonBisection_impl<f>(x0, p, settings, x);
+   real_t dfdx = __enzyme_fwddiff<real_t>((void*)+f, enzyme_dup, *x, 1.0,
+                                          enzyme_const, *p);
+   real_t dfdp = __enzyme_fwddiff<real_t>((void*)+f, enzyme_const, *x, enzyme_dup,
+                                          *p, *dp);
+   *dx = -dfdp/dfdx;
+}
+
+
+template<auto f, typename T>
+void SolveNewtonBisection_impl_aug(const real_t* x0, real_t* x0_bar,
+                                   const T* p, T* p_bar,
+                                   const SolverSettings* settings, SolverSettings* settings_bar,
+                                   real_t* x, real_t* x_bar)
+{
+   SolveNewtonBisection_impl<f>(x0, p, settings, x);
+}
+
+// Change the residual function to return-by-reference so that there is a
+// slot to provide the downstream cotangent (ie the shadow for y)
+// in the reverse mode call.
+template<auto f, typename T>
+void rbr_wrapper(real_t x, T& p, real_t& y)
+{
+   y = f(x, p);
+}
+
+template<auto f, typename T>
+void SolveNewtonBisection_impl_rev(const real_t* x0, real_t* x0_bar,
+                                   const T* p, T* p_bar,
+                                   const SolverSettings* settings, SolverSettings* settings_bar,
+                                   real_t* x, real_t* x_bar)
+{
+   real_t drdx = __enzyme_fwddiff<real_t>((void*)+f, enzyme_dup, *x, 1.0,
+                                          enzyme_const, *p);
+   real_t lambda = -(*x_bar / drdx);
+   real_t r;
+   __enzyme_autodiff<void>((void*)rbr_wrapper<f, T>, enzyme_const, *x, enzyme_dup,
+                           p, p_bar, enzyme_dupnoneed, &r, &lambda);
+
+   // These are logically constants, the root has no sensitivity to these
+   *x0_bar = 0.0;
+   *settings_bar = SolverSettings{};
+}
+
+/// @endcond
+} // namespace internal
+
+
+namespace future
+{
+
+/**
+ * @brief Find the root of a univariate funtion
+ */
+template<auto f, typename T>
+MFEM_HOST_DEVICE __attribute__((always_inline)) real_t SolveNewtonBisection(
+   real_t x0, T p, SolverSettings settings)
+{
+   // We need to tell Enzyme how much memory in the SolverSettings object is
+   // used by active variables (in the sense of Enzyme activity analysis).
+   // Without this, it seems that a bug in LLVM causes this information to
+   // be lost during some optimization pass, and the Enzyme pass fails in
+   // Release builds.
+   // There are 4 real_t members in settings, which is what Enzyme will
+   // consider active.
+   // TODO: File an issue on Enzyme to remind Bill to fix this in LLVM.
+   __enzyme_double((void*)&settings, sizeof(real_t)*4);
+
+   real_t x;
+   internal::SolveNewtonBisection_impl<f>(&x0, &p, &settings, &x);
+   return x;
+}
+
+
+} // namespace future
+} // namespace mfem
+
+#endif // MFEM_USE_ENZYME
+#endif // MFEM_UNIVARSOLVERS
@@ -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
@@ -349,7 +349,7 @@ public:
       vector-valued finite elements, which is also the width of the
       DenseMatrix argument in
       CalcPhysVShape(ElementTransformation &Trans, DenseMatrix &shape). */
-   virtual int GetPhysRangeDim(int /* space_dim */) const { return vdim; }
+   int GetPhysRangeDim(int /* space_dim */) const { return vdim; }

   /** Returns the dimension of the curl for vector-valued finite elements,
       which is also the width of the DenseMatrix argument in
@@ -360,7 +360,7 @@ public:
       finite elements, which is also the width of the DenseMatrix argument in
       CalcPhysCurlShape(ElementTransformation &Trans, DenseMatrix &curl_shape).
   */
-   virtual int GetPhysCurlDim(int /* space_dim */) const { return cdim; }
+   int GetPhysCurlDim(int /* space_dim */) const { return cdim; }

   /// Returns the Geometry::Type of the reference element.
   Geometry::Type GetGeomType() const { return geom_type; }
@@ -1017,7 +1017,7 @@ public:
   VectorFiniteElement(int D, Geometry::Type G, int Do, int O, int M,
                       int F = FunctionSpace::Pk);

-   int GetPhysRangeDim(int space_dim) const override { return space_dim; }
+   int GetPhysRangeDim(int space_dim) const { return space_dim; }
 };

 /// @brief Class for computing 1D special polynomials and their associated basis
@@ -663,8 +663,8 @@ public:
                         const int cb_type = BasisType::GaussLobatto,
                         const int ob_type = BasisType::GaussLegendre);

-   int GetPhysRangeDim(int space_dim) const override { return 2; }
-   int GetPhysCurlDim(int space_dim) const override { return 1; }
+   int GetPhysRangeDim(int space_dim) const { return 2; }
+   int GetPhysCurlDim(int space_dim) const { return 1; }

   void CalcVShape(const IntegrationPoint &ip,
                   DenseMatrix &shape) const override;
@@ -708,8 +708,8 @@ private:
                           DenseMatrix &I) const;

 public:
-   int GetPhysRangeDim(int space_dim) const override { return 3; }
-   int GetPhysCurlDim(int space_dim) const override { return 3; }
+   int GetPhysRangeDim(int space_dim) const { return 3; }
+   int GetPhysCurlDim(int space_dim) const { return 3; }

   using FiniteElement::CalcVShape;
   using FiniteElement::CalcPhysCurlShape;
@@ -510,8 +510,8 @@ public:
   RT_R2D_SegmentElement(const int p,
                         const int ob_type = BasisType::GaussLegendre);

-   int GetPhysRangeDim(int space_dim) const override { return 2; }
-   int GetPhysCurlDim(int space_dim) const override { return 0; }
+   int GetPhysRangeDim(int space_dim) const { return 2; }
+   int GetPhysCurlDim(int space_dim) const { return 0; }

   void CalcVShape(const IntegrationPoint &ip,
                   DenseMatrix &shape) const override;
@@ -550,8 +550,8 @@ private:
                           DenseMatrix &I) const;

 public:
-   int GetPhysRangeDim(int space_dim) const override { return 3; }
-   int GetPhysCurlDim(int space_dim) const override { return 0; }
+   int GetPhysRangeDim(int space_dim) const { return 3; }
+   int GetPhysCurlDim(int space_dim) const { return 0; }

   using FiniteElement::CalcVShape;

@@ -53,6 +53,7 @@
 #include "particleset.hpp"

 #include "dfem/doperator.hpp"
+#include "dfem/univarsolvers.hpp"

 #ifdef MFEM_USE_MPI
 #include "pfespace.hpp"
@@ -3934,16 +3934,6 @@ const FiniteElement *FiniteElementSpace::GetBE(int i) const
   return BE;
 }

-const FiniteElement *FiniteElementSpace::GetTypicalBE() const
-{
-   if (mesh->GetNBE() > 0) { return GetBE(0); }
-
-   Geometry::Type geom = mesh->GetTypicalFaceGeometry();
-   const FiniteElement *be = fec->FiniteElementForGeometry(geom);
-   MFEM_VERIFY(be != nullptr, "Could not determine a typical BE!");
-   return be;
-}
-
 const FiniteElement *FiniteElementSpace::GetFaceElement(int i) const
 {
   MFEM_VERIFY(!IsVariableOrder(), "not implemented");
@@ -3974,11 +3964,6 @@ const FiniteElement *FiniteElementSpace::GetFaceElement(int i) const
   return fe;
 }

-const FiniteElement *FiniteElementSpace::GetTypicalFaceElement() const
-{
-   return fec->FiniteElementForGeometry(mesh->GetTypicalFaceGeometry());
-}
-
 const FiniteElement *FiniteElementSpace::GetEdgeElement(int i,
                                                        int variant) const
 {
@@ -839,7 +839,7 @@ public:
       Note: For vector-valued elements, the results pads up the range dimension
       to the spatial dimension. E.g., consider a stack of 5 vector-valued
       elements each representing 2D vectors, living in a 3 dimensional space.
-       Then this function would give 15, not 10.
+       Then this fucntion would give 15, not 10.
       */
   int GetVectorDim() const;

@@ -1323,24 +1323,12 @@ public:
        associated with i'th boundary face in the mesh object. */
   const FiniteElement *GetBE(int i) const;

-   /// @brief Return a "typical" boundary element.
-   ///
-   /// This can be used in situations where the local mesh partition may be
-   /// empty.
-   const FiniteElement *GetTypicalBE() const;
-
   /** @brief Returns pointer to the FiniteElement in the FiniteElementCollection
        associated with i'th face in the mesh object.  Faces in this case refer
        to the MESHDIM-1 primitive so in 2D they are segments and in 1D they are
        points.*/
   const FiniteElement *GetFaceElement(int i) const;

-   /// @brief Return a "typical" face element.
-   ///
-   /// This can be used in situations where the local mesh partition may be
-   /// empty.
-   const FiniteElement *GetTypicalFaceElement() const;
-
   /** @brief Returns pointer to the FiniteElement in the FiniteElementCollection
        associated with i'th edge in the mesh object. */
   const FiniteElement *GetEdgeElement(int i, int variant = 0) const;
@@ -345,6 +345,27 @@ void GridFunction::ComputeFlux(BilinearFormIntegrator &blfi,
   }
 }

+int GridFunction::VectorDim() const
+{
+   const FiniteElement *fe = fes->GetTypicalFE();
+   if (!fe || fe->GetRangeType() == FiniteElement::SCALAR)
+   {
+      return fes->GetVDim();
+   }
+   return fes->GetVDim()*std::max(fes->GetMesh()->SpaceDimension(),
+                                  fe->GetRangeDim());
+}
+
+int GridFunction::CurlDim() const
+{
+   const FiniteElement *fe = fes->GetTypicalFE();
+   if (!fe || fe->GetRangeType() == FiniteElement::SCALAR)
+   {
+      return 2 * fes->GetMesh()->SpaceDimension() - 3;
+   }
+   return fes->GetVDim()*fe->GetCurlDim();
+}
+
 void GridFunction::GetTrueDofs(Vector &tv) const
 {
   const SparseMatrix *R = fes->GetRestrictionMatrix();
@@ -2029,18 +2050,6 @@ void GridFunction::AccumulateAndCountBdrValues(
   Coefficient *coeff[], VectorCoefficient *vcoeff, const Array<int> &attr,
   Array<int> &values_counter)
 {
-   if (vcoeff)
-   {
-      MFEM_VERIFY(fes->GetVDim() == vcoeff->GetVDim(),
-                  "vcoeff vdim != fes VDim");
-      MFEM_VERIFY(fes->GetTypicalBE()->GetMapType() == FiniteElement::VALUE &&
-                  fes->GetTypicalBE()->GetRangeType() ==
-                  FiniteElement::SCALAR,
-                  "Can only call ProjectBdrCoefficient on scalar value-type "
-                  "boundary elements. "
-                  "Did you intended to call ProjectBdrCoefficientNormal or "
-                  "ProjectBdrCoefficientTangent for vector finite elements?");
-   }
   Array<int> vdofs;
   Vector vc;

@@ -2193,9 +2202,6 @@ void GridFunction::AccumulateAndCountBdrTangentValues(
   VectorCoefficient &vcoeff, const Array<int> &bdr_attr,
   Array<int> &values_counter)
 {
-   MFEM_VERIFY(fes->GetTypicalBE()->GetPhysRangeDim(
-                  fes->GetMesh()->SpaceDimension()) == vcoeff.GetVDim(),
-               "vcoeff vdim != PhysRangeDim");
   const FiniteElement *fe;
   ElementTransformation *T;
   Array<int> dofs;
@@ -2349,9 +2355,6 @@ void GridFunction::ProjectDeltaCoefficient(DeltaCoefficient &delta_coeff,

 void GridFunction::ProjectCoefficient(Coefficient &coeff, ProjectType type)
 {
-   MFEM_VERIFY(
-      VectorDim() == 1,
-      "Cannot project scalar Coefficient onto vector GridFunction");
   DeltaCoefficient *delta_c = dynamic_cast<DeltaCoefficient *>(&coeff);
   DofTransformation doftrans;
   Array<int> vdofs;
@@ -2627,7 +2630,6 @@ void GridFunction::ProjectCoefficient(
 void GridFunction::ProjectCoefficient(VectorCoefficient &vcoeff,
                                      ProjectType type)
 {
-   MFEM_VERIFY(VectorDim() == vcoeff.GetVDim(), "vcoeff vdim != VectorDim()");
   Array<int> vdofs;
   Vector vals;
   DofTransformation doftrans;
@@ -2943,7 +2945,6 @@ void GridFunction::ProjectCoefficientElementL2(VectorCoefficient &vcoeff)
 void GridFunction::ProjectCoefficient(
   VectorCoefficient &vcoeff, Array<int> &dofs)
 {
-   MFEM_VERIFY(VectorDim() == vcoeff.GetVDim(), "vcoeff vdim != VectorDim()");
   int el = -1;
   ElementTransformation *T = NULL;
   const FiniteElement *fe = NULL;
@@ -2973,7 +2974,6 @@ void GridFunction::ProjectCoefficient(

 void GridFunction::ProjectCoefficient(VectorCoefficient &vcoeff, int attribute)
 {
-   MFEM_VERIFY(VectorDim() == vcoeff.GetVDim(), "vcoeff vdim != VectorDim()");
   int i;
   Array<int> vdofs;
   Vector vals;
@@ -3030,14 +3030,9 @@ void GridFunction::ProjectCoefficient(Coefficient *coeff[])
   }
 }

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

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

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

+void GridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff)
+{
+   Array<int> dof_attr;
+   ProjectDiscCoefficient(coeff, dof_attr);
+}
+
 void GridFunction::ProjectDiscCoefficient(Coefficient &coeff, AvgType type)
 {
   // Harmonic  (x1 ... xn) = [ (1/x1 + ... + 1/xn) / n ]^-1.
   // Arithmetic(x1 ... xn) = (x1 + ... + xn) / n.

-   MFEM_VERIFY(
-      VectorDim() == 1,
-      "Cannot project a scalar coefficient onto a vector GridFunction");
-
   Array<int> zones_per_vdof;
   AccumulateAndCountZones(coeff, type, zones_per_vdof);

@@ -3088,7 +3082,6 @@ void GridFunction::ProjectDiscCoefficient(Coefficient &coeff, AvgType type)
 void GridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff,
                                          AvgType type)
 {
-   MFEM_VERIFY(VectorDim() == coeff.GetVDim(), "coeff vdim != VectorDim()");
   Array<int> zones_per_vdof;
   AccumulateAndCountZones(coeff, type, zones_per_vdof);

@@ -3144,33 +3137,52 @@ void GridFunction::ProjectBdrCoefficient(Coefficient *coeff[],
 }

 void GridFunction::ProjectBdrCoefficientNormal(
-   Coefficient *coeff, VectorCoefficient *vcoeff, const Array<int> &bdr_attr)
+   VectorCoefficient &vcoeff, const Array<int> &bdr_attr)
 {
-   MFEM_VERIFY(fes->GetVDim() == 1, "fespace VDim != 1");
-   MFEM_VERIFY(fes->GetTypicalBE()->GetRangeType() == FiniteElement::SCALAR &&
-               fes->GetTypicalBE()->GetMapType() == FiniteElement::INTEGRAL,
-               "Not an RT FE space!");
-   if (vcoeff)
-   {
-      MFEM_VERIFY(vcoeff->GetVDim() == fes->GetMesh()->SpaceDimension(),
-                  "vcoeff vdim (" << vcoeff->GetVDim()
-                  << ") != SpaceDimension ("
-                  << fes->GetMesh()->SpaceDimension() << ")");
-   }
+#if 0
+   // implementation for the case when the face dofs are integrals of the
+   // normal component.
+   const FiniteElement *fe;
+   ElementTransformation *T;
+   Array<int> dofs;
+   int dim = vcoeff.GetVDim();
+   Vector vc(dim), nor(dim), lvec, shape;

+   for (int i = 0; i < fes->GetNBE(); i++)
+   {
+      if (bdr_attr[fes->GetBdrAttribute(i)-1] == 0)
+      {
+         continue;
+      }
+      fe = fes->GetBE(i);
+      T = fes->GetBdrElementTransformation(i);
+      int intorder = 2*fe->GetOrder(); // !!!
+      const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(), intorder);
+      int nd = fe->GetDof();
+      lvec.SetSize(nd);
+      shape.SetSize(nd);
+      lvec = 0.0;
+      for (int j = 0; j < ir.GetNPoints(); j++)
+      {
+         const IntegrationPoint &ip = ir.IntPoint(j);
+         T->SetIntPoint(&ip);
+         vcoeff.Eval(vc, *T, ip);
+         CalcOrtho(T->Jacobian(), nor);
+         fe->CalcShape(ip, shape);
+         lvec.Add(ip.weight * (vc * nor), shape);
+      }
+      fes->GetBdrElementDofs(i, dofs);
+      SetSubVector(dofs, lvec);
+   }
+#else
   // implementation for the case when the face dofs are scaled point
   // values of the normal component.
   const FiniteElement *fe;
   ElementTransformation *T;
   Array<int> dofs;
-   Vector vc, nor, lvec;
+   int dim = vcoeff.GetVDim();
+   Vector vc(dim), nor(dim), lvec;
   DofTransformation doftrans;
-   if (vcoeff)
-   {
-      const int dim = vcoeff->GetVDim();
-      vc.SetSize(dim);
-      nor.SetSize(dim);
-   }

   for (int i = 0; i < fes->GetNBE(); i++)
   {
@@ -3186,22 +3198,15 @@ void GridFunction::ProjectBdrCoefficientNormal(
      {
         const IntegrationPoint &ip = ir.IntPoint(j);
         T->SetIntPoint(&ip);
-         if (coeff)
-         {
-            const real_t c = coeff->Eval(*T, ip);
-            lvec(j) = c * T->Weight();
-         }
-         else if (vcoeff)
-         {
-            vcoeff->Eval(vc, *T, ip);
-            CalcOrtho(T->Jacobian(), nor);
-            lvec(j) = (vc * nor);
-         }
+         vcoeff.Eval(vc, *T, ip);
+         CalcOrtho(T->Jacobian(), nor);
+         lvec(j) = (vc * nor);
      }
      fes->GetBdrElementDofs(i, dofs, doftrans);
      doftrans.TransformPrimal(lvec);
      SetSubVector(dofs, lvec);
   }
+#endif
 }

 void GridFunction::ProjectBdrCoefficientTangent(
@@ -5002,14 +5007,6 @@ real_t ExtrudeCoefficient::Eval(ElementTransformation &T,
   return sol_in.Eval(*T_in, ip);
 }

-void VectorExtrudeCoefficient::Eval(Vector &v, ElementTransformation &T,
-                                    const IntegrationPoint &ip)
-{
-   ElementTransformation *T_in =
-      mesh_in->GetElementTransformation(T.ElementNo / n);
-   T_in->SetIntPoint(&ip);
-   sol_in.Eval(v, *T_in, ip);
-}

 GridFunction *Extrude1DGridFunction(Mesh *mesh, Mesh *mesh2d,
                                    GridFunction *sol, const int ny)
@@ -5060,17 +5057,10 @@ GridFunction *Extrude1DGridFunction(Mesh *mesh, Mesh *mesh2d,
      return NULL;
   }
   FiniteElementSpace *solfes2d;
-   const int vdim = sol->FESpace()->GetVDim();
-   solfes2d = new FiniteElementSpace(mesh2d, solfec2d, vdim);
+   // assuming sol is scalar
+   solfes2d = new FiniteElementSpace(mesh2d, solfec2d);
   sol2d = new GridFunction(solfes2d);
   sol2d->MakeOwner(solfec2d);
-   if (vdim > 1)
-   {
-      VectorGridFunctionCoefficient vcsol(sol);
-      VectorExtrudeCoefficient vc2d(mesh, vcsol, ny);
-      sol2d->ProjectCoefficient(vc2d);
-   }
-   else
   {
      GridFunctionCoefficient csol(sol);
      ExtrudeCoefficient c2d(mesh, csol, ny);
@@ -5768,4 +5758,4 @@ std::pair<real_t, real_t> GridFunction::EstimateFunctionMaximum(
   return std::make_pair(global_max_lower, global_max_upper);
 }

-}
+}
@@ -23,7 +23,6 @@
 #include <limits>
 #include <ostream>
 #include <string>
-#include <variant>

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

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

   /** Helper function for ProjectCoefficientElementL2 */
   void ProjectCoefficientElementL2_(Coefficient &coeff, Vector &sol, Vector &Va);
@@ -159,13 +150,11 @@ public:

   FiniteElementCollection *OwnFEC() { return fec_owned; }

-   /** @brief Shortcut for calling FiniteElementSpace::GetVectorDim() on the
-       underlying #fes */
-   int VectorDim() const { return fes->GetVectorDim(); }
+   /// Shortcut for calling FiniteElementSpace::GetVectorDim() on the underlying #fes
+   int VectorDim() const;

-   /** @brief Shortcut for calling FiniteElementSpace::GetCurlDim() on the
-       underlying #fes */
-   int CurlDim() const { return fes->GetCurlDim(); }
+   /// Shortcut for calling FiniteElementSpace::GetCurlDim() on the underlying #fes
+   int CurlDim() const;

   /// Read only access to the (optional) internal true-dof Vector.
   const Vector &GetTrueVector() const
@@ -524,17 +513,10 @@ public:
       but using an array of scalar coefficients for each component. */
   void ProjectCoefficient(Coefficient *coeff[]);

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

   enum AvgType {ARITHMETIC, HARMONIC};
   /** @brief Projects a discontinuous coefficient so that the values in shared
@@ -550,9 +532,6 @@ public:
   std::unique_ptr<GridFunction> ProlongateToMaxOrder() const;

 protected:
-   void ProjectBdrCoefficientNormal(Coefficient *coeff, VectorCoefficient *vcoeff,
-                                    const Array<int> &attr);
-
   /** @brief Accumulates (depending on @a type) the values of @a coeff at all
       shared vdofs and counts in how many zones each vdof appears. */
   void AccumulateAndCountZones(Coefficient &coeff, AvgType type,
@@ -677,26 +656,15 @@ public:
   virtual void ProjectBdrCoefficient(Coefficient *coeff[],
                                      const Array<int> &attr);

-   /** @brief Project the normal component of the given VectorCoefficient on
-       the boundary. */
-   /** Only boundary attributes that are marked in @a bdr_attr are
-       projected. Assumes RT-type vector finite element GridFunction. */
+   /** Project the normal component of the given VectorCoefficient on
+       the boundary. Only boundary attributes that are marked in
+       'bdr_attr' are projected. Assumes RT-type VectorFE GridFunction. */
   void ProjectBdrCoefficientNormal(VectorCoefficient &vcoeff,
-                                    const Array<int> &bdr_attr)
-   { ProjectBdrCoefficientNormal(NULL, &vcoeff, bdr_attr); }
-
-   /** @brief Project the given Coefficient in the normal direction on the
-       boundary. */
-   /** Only boundary attributes that are marked in @a bdr_attr are projected.
-       Assumes RT-type vector finite element GridFunction. */
-   void ProjectBdrCoefficientNormal(Coefficient &coeff,
-                                    const Array<int> &bdr_attr)
-   { ProjectBdrCoefficientNormal(&coeff, NULL, bdr_attr); }
+                                    const Array<int> &bdr_attr);

   /** @brief Project the tangential components of the given VectorCoefficient
-       on the boundary. */
-   /** Only boundary attributes that are marked in @a bdr_attr
-       are projected. Assumes ND-type vector finite element GridFunction. */
+       on the boundary. Only boundary attributes that are marked in @a bdr_attr
+       are projected. Assumes ND-type VectorFE GridFunction. */
   virtual void ProjectBdrCoefficientTangent(VectorCoefficient &vcoeff,
                                             const Array<int> &bdr_attr);

@@ -1946,7 +1914,7 @@ real_t ComputeElementLpDistance(real_t p, int i,
                                GridFunction& gf1, GridFunction& gf2);


-/// Class used for extruding a scalar coefficient
+/// Class used for extruding scalar GridFunctions
 class ExtrudeCoefficient : public Coefficient
 {
 private:
@@ -1954,53 +1922,13 @@ private:
   Mesh *mesh_in;
   Coefficient &sol_in;
 public:
-   /// Constructs an instance of VectorExtrudeCoefficient
-   /**
-    * @param m      1D mesh
-    * @param s      1D vector coefficient
-    * @param n_     number of transverse elements of the extruded mesh
-    */
   ExtrudeCoefficient(Mesh *m, Coefficient &s, int n_)
-      : n(n_), mesh_in(m), sol_in(s)
-   { MFEM_VERIFY(n > 0, "Number of transverse elements must be positive!"); }
-
+      : n(n_), mesh_in(m), sol_in(s) { }
   real_t Eval(ElementTransformation &T, const IntegrationPoint &ip) override;
-
   virtual ~ExtrudeCoefficient() { }
 };

-/// Class used for extruding a vector coefficient
-class VectorExtrudeCoefficient : public VectorCoefficient
-{
-private:
-   int n;
-   Mesh *mesh_in;
-   VectorCoefficient &sol_in;
-public:
-   /// Constructs an instance of VectorExtrudeCoefficient
-   /**
-    * @param m      1D mesh
-    * @param s      1D vector coefficient
-    * @param n_     number of transverse elements of the extruded mesh
-    */
-   VectorExtrudeCoefficient(Mesh *m, VectorCoefficient &s, int n_)
-      : VectorCoefficient(s.GetVDim()), n(n_), mesh_in(m), sol_in(s)
-   { MFEM_VERIFY(n > 0, "Number of transverse elements must be positive!"); }
-
-   void Eval(Vector &v, ElementTransformation &T,
-             const IntegrationPoint &ip) override;
-   using VectorCoefficient::Eval;
-
-   virtual ~VectorExtrudeCoefficient() { }
-};
-
-/// Extrude a 1D GridFunction, after extruding the mesh with Extrude1D()
-/**
- * @param mesh      1D mesh
- * @param mesh2d    extruded mesh
- * @param sol       grid function
- * @param ny        number of transverse elements of the extruded mesh
- */
+/// Extrude a scalar 1D GridFunction, after extruding the mesh with Extrude1D.
 GridFunction *Extrude1DGridFunction(Mesh *mesh, Mesh *mesh2d,
                                    GridFunction *sol, const int ny);

@@ -490,7 +490,7 @@ void FindPointsGSLIB::FindPointsOnDevice(const Vector &point_pos,
   }
   DEV.find_device = true;

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

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

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

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

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

@@ -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]]
@@ -197,21 +197,15 @@ static void EAHdivAssemble3D(const int NE,
      // Assemble (one row per thread)
      MFEM_FOREACH_THREAD(idx_i, x, NDOF)
      {
-         // NOTE: due to an llvm backend bug, usage of the modulus operator
-         // has been removed from this foreach section.
         const int ic = idx_i / NDOF_C;
-         const int idx_ii = idx_i - ic * NDOF_C; // idx_i % NDOF_C
+         const int idx_ii = idx_i % NDOF_C;

         const int nx_i = (ic == 0) ? D1D : D1D-1;
         const int ny_i = (ic == 1) ? D1D : D1D-1;

-         const int qx_i = idx_ii / nx_i;
-         const int ix = idx_ii - qx_i * nx_i; // idx_ii % nx_i
-
-         const int qy_i = qx_i / ny_i;
-         const int iy = qx_i - qy_i * ny_i; // (idx_ii / nx_i) % ny_i
-
-         const int iz = qy_i; // (idx_ii / nx_i) / ny_i
+         const int ix = idx_ii % nx_i;
+         const int iy = (idx_ii / nx_i) % ny_i;
+         const int iz = (idx_ii / nx_i) / ny_i;

         const real_t (&Bi1)[MQ1][MD1] = (ic == 0) ? r_Bc : r_Bo;
         const real_t (&Bi2)[MQ1][MD1] = (ic == 1) ? r_Bc : r_Bo;
@@ -220,18 +214,14 @@ static void EAHdivAssemble3D(const int NE,
         for (int idx_j = 0; idx_j < NDOF; ++idx_j)
         {
            const int jc = idx_j / NDOF_C;
-            const int idx_jj = idx_j - jc * NDOF_C; // idx_j % NDOF_C
+            const int idx_jj = idx_j % NDOF_C;

            const int nx_j = (jc == 0) ? D1D : D1D-1;
            const int ny_j = (jc == 1) ? D1D : D1D-1;

-            const int qx_j = idx_jj / nx_j;
-            const int jx = idx_jj - qx_j * nx_j; // idx_jj % nx_j
-
-            const int qy_j = qx_j / ny_j;
-            const int jy = qx_j - qy_j * ny_j; // (idx_jj / nx_j) % ny_j
-
-            const int jz = qy_j; // (idx_jj / nx_j) / ny_j
+            const int jx = idx_jj % nx_j;
+            const int jy = (idx_jj / nx_j) % ny_j;
+            const int jz = (idx_jj / nx_j) / ny_j;

            const real_t (&Bj1)[MQ1][MD1] = (jc == 0) ? r_Bc : r_Bo;
            const real_t (&Bj2)[MQ1][MD1] = (jc == 1) ? r_Bc : r_Bo;
@@ -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.
@@ -121,6 +125,18 @@ private:
   void AddTriPoints3b(const int off, const real_t b, const real_t weight)
   { AddTriPoints3(off, (1. - b)/2., b, weight); }

+   void AddTriPoints3R(const int off, const real_t a, const real_t b,
+                       const real_t c, const real_t weight)
+   {
+      IntPoint(off + 0).Set2w(a, b, weight);
+      IntPoint(off + 1).Set2w(c, a, weight);
+      IntPoint(off + 2).Set2w(b, c, weight);
+   }
+
+   void AddTriPoints3R(const int off, const real_t a, const real_t b,
+                       const real_t weight)
+   { AddTriPoints3R(off, a, b, 1. - a - b, weight); }
+
   void AddTriPoints6(const int off, const real_t a, const real_t b,
                      const real_t c, const real_t weight)
   {
@@ -167,6 +183,14 @@ private:
      AddTetPoints3(off + 1, a, 1. - 3.*a, weight);
   }

+   // given b, add the permutations of (a,a,a,b), where 3*a + b = 1
+   void AddTetPoints4b(const int off, const real_t b, const real_t weight)
+   {
+      const real_t a = (1. - b)/3.;
+      IntPoint(off).Set(a, a, a, weight);
+      AddTetPoints3(off + 1, a, b, weight);
+   }
+
   // add the permutations of (a,a,b,b), 2*(a + b) = 1
   void AddTetPoints6(const int off, const real_t a, const real_t weight)
   {
@@ -185,37 +209,14 @@ private:
      AddTetPoints6(off + 6, a, bc, cb, weight);
   }

-   // add all 24 permutations of (a,b,c,d) where a+b+c+d = 1, all distinct
-   void AddTetPoints24(const int off, const real_t a, const real_t b,
-                       const real_t c, const real_t weight)
+   // given (b,c), add the permutations of (a,a,b,c), 2*a + b + c = 1
+   void AddTetPoints12bc(const int off, const real_t b, const real_t c,
+                         const real_t weight)
   {
-      const real_t d = 1. - a - b - c;
-      // all 24 permutations of 4 distinct barycentric coordinates
-      // permuting which coordinate goes to x, y, z (4th is 1-x-y-z)
-      IntPoint(off +  0).Set(a, b, c, weight);
-      IntPoint(off +  1).Set(a, b, d, weight);
-      IntPoint(off +  2).Set(a, c, b, weight);
-      IntPoint(off +  3).Set(a, c, d, weight);
-      IntPoint(off +  4).Set(a, d, b, weight);
-      IntPoint(off +  5).Set(a, d, c, weight);
-      IntPoint(off +  6).Set(b, a, c, weight);
-      IntPoint(off +  7).Set(b, a, d, weight);
-      IntPoint(off +  8).Set(b, c, a, weight);
-      IntPoint(off +  9).Set(b, c, d, weight);
-      IntPoint(off + 10).Set(b, d, a, weight);
-      IntPoint(off + 11).Set(b, d, c, weight);
-      IntPoint(off + 12).Set(c, a, b, weight);
-      IntPoint(off + 13).Set(c, a, d, weight);
-      IntPoint(off + 14).Set(c, b, a, weight);
-      IntPoint(off + 15).Set(c, b, d, weight);
-      IntPoint(off + 16).Set(c, d, a, weight);
-      IntPoint(off + 17).Set(c, d, b, weight);
-      IntPoint(off + 18).Set(d, a, b, weight);
-      IntPoint(off + 19).Set(d, a, c, weight);
-      IntPoint(off + 20).Set(d, b, a, weight);
-      IntPoint(off + 21).Set(d, b, c, weight);
-      IntPoint(off + 22).Set(d, c, a, weight);
-      IntPoint(off + 23).Set(d, c, b, weight);
+      const real_t a = (1. - b - c)/2.;
+      AddTetPoints3(off,     a, b, weight);
+      AddTetPoints3(off + 3, a, c, weight);
+      AddTetPoints6(off + 6, a, b, c, weight);
   }

 public:
@@ -297,8 +297,7 @@ void LinearForm::Assemble()
         tr = mesh->GetBdrFaceTransformations(i);
         if (tr != NULL)
         {
-            mfem::DofTransformation doftrans;
-            fes -> GetElementVDofs (tr -> Elem1No, vdofs, doftrans);
+            fes -> GetElementVDofs (tr -> Elem1No, vdofs);
            for (int k = 0; k < boundary_face_integs.Size(); k++)
            {
               if (boundary_face_integs_marker[k] &&
@@ -308,7 +307,6 @@ void LinearForm::Assemble()
               boundary_face_integs[k]->
               AssembleRHSElementVect(*fes->GetFE(tr->Elem1No),
                                      *tr, elemvect);
-               doftrans.TransformDual(elemvect);
               AddElementVector (vdofs, elemvect);
            }
         }
@@ -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; }

@@ -545,8 +545,6 @@ void ParGridFunction::GetElementDofValues(int el, Vector &dof_vals) const

 void ParGridFunction::ProjectCoefficient(Coefficient &coeff, ProjectType type)
 {
-   MFEM_VERIFY(VectorDim() == 1,
-               "Cannot project scalar coefficient onto vector ParGridFunction");
   DeltaCoefficient *delta_c = dynamic_cast<DeltaCoefficient *>(&coeff);

   if (delta_c == NULL)
@@ -717,8 +715,7 @@ void ParGridFunction::ProjectCoefficientElementL2(VectorCoefficient &vcoeff)
 }


-void ParGridFunction::ProjectDiscCoefficient(
-   std::variant<Coefficient*, VectorCoefficient*> coeff)
+void ParGridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff)
 {
   // local maximal element attribute for each dof
   Array<int> ldof_attr;
@@ -764,9 +761,6 @@ void ParGridFunction::ProjectDiscCoefficient(

 void ParGridFunction::ProjectDiscCoefficient(Coefficient &coeff, AvgType type)
 {
-   MFEM_VERIFY(
-      VectorDim() == 1,
-      "Cannot project scalar coefficient onto a vector ParGridFunction");
   // Harmonic  (x1 ... xn) = [ (1/x1 + ... + 1/xn) / n ]^-1.
   // Arithmetic(x1 ... xn) = (x1 + ... + xn) / n.

@@ -792,8 +786,6 @@ void ParGridFunction::ProjectDiscCoefficient(VectorCoefficient &vcoeff,
   // Harmonic  (x1 ... xn) = [ (1/x1 + ... + 1/xn) / n ]^-1.
   // Arithmetic(x1 ... xn) = (x1 + ... + xn) / n.

-   MFEM_VERIFY(VectorDim() == vcoeff.GetVDim(), "vcoeff vdim != VectorDim()");
-
   // Number of zones that contain a given dof.
   Array<int> zones_per_vdof;
   AccumulateAndCountZones(vcoeff, type, zones_per_vdof);
@@ -866,12 +858,6 @@ void ParGridFunction::ProjectBdrCoefficient(
 #endif
 }

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

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

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

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

   void ProjectDiscCoefficient(VectorCoefficient &vcoeff, AvgType type) override;
@@ -281,7 +280,8 @@ public:
   using GridFunction::ProjectBdrCoefficient;

   void ProjectBdrCoefficient(VectorCoefficient &vcoeff,
-                              const Array<int> &attr) override;
+                              const Array<int> &attr) override
+   { ProjectBdrCoefficient(NULL, &vcoeff, attr); }

   void ProjectBdrCoefficient(Coefficient *coeff[],
                              const Array<int> &attr) override
@@ -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,31 @@ 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);
      }

-      // Read the array
-      data[ArrayName].Load(in);
+      // Ignore the remainder of the line which may contain explanatory comments
+      data[ArrayName].Load(in, 0);
   }
+
 }

 }
@@ -726,16 +726,16 @@ std::string Device::GetUUID(const int device_id)
   MFEM_GPU_CHECK(cudaGetDeviceProperties(&prop, device_id));
   for (int i = 0; i < 16; ++i)
   {
-      const unsigned b = static_cast<unsigned char>(prop.uuid.bytes[i]);
-      res << std::setfill('0') << std::setw(2) << std::hex << b;
+      res << std::setfill('0') << std::setw(2) << std::hex
+          << static_cast<unsigned>(prop.uuid.bytes[i]);
   }
 #elif defined(MFEM_USE_HIP)
   hipUUID uuid;
   MFEM_GPU_CHECK(hipDeviceGetUuid(&uuid, device_id));
   for (int i = 0; i < 16; ++i)
   {
-      const unsigned b = static_cast<unsigned char>(uuid.bytes[i]);
-      res << std::setfill('0') << std::setw(2) << std::hex << b;
+      res << std::setfill('0') << std::setw(2) << std::hex
+          << static_cast<unsigned>(uuid.bytes[i]);
   }
 #endif
   return res.str();
@@ -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)
 {
@@ -23,7 +23,6 @@ list(APPEND SRCS
  complex_operator.cpp
  constraints.cpp
  densemat.cpp
-  eigensolvers.cpp
  filteredsolver.cpp
  handle.cpp
  matrix.cpp
@@ -56,7 +55,6 @@ list(APPEND HDRS
  dinvariants.hpp
  dtensor.hpp
  dual.hpp
-  eigensolvers.hpp
  filteredsolver.hpp
  handle.hpp
  invariants.hpp
@@ -103,11 +101,6 @@ if (MFEM_USE_MPI)
  endif()
 endif()

-if (MFEM_USE_ARPACK)
-  list(APPEND SRCS arpack.cpp)
-  list(APPEND HDRS arpack.hpp)
-endif()
-
 if (MFEM_USE_SUNDIALS)
  list(APPEND SRCS sundials.cpp)
  list(APPEND HDRS sundials.hpp)
@@ -1,271 +0,0 @@
-// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
-// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-// LICENSE and NOTICE for details. LLNL-CODE-806117.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability visit https://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the BSD-3 license. We welcome feedback and contributions, see file
-// CONTRIBUTING.md for details.
-
-#ifndef MFEM_ARPACK
-#define MFEM_ARPACK
-
-#include "../config/config.hpp"
-
-#ifdef MFEM_USE_ARPACK
-
-#include <string>
-
-#ifdef MFEM_USE_MPI
-#include <mpi.h>
-#include "hypre.hpp"
-#endif
-
-#include "operator.hpp"
-
-#define SSAUPD  ssaupd_
-#define SSEUPD  sseupd_
-#define DSAUPD  dsaupd_
-#define DSEUPD  dseupd_
-
-#ifdef MFEM_USE_MPI
-#define PSSAUPD pssaupd_
-#define PSSEUPD psseupd_
-#define PDSAUPD pdsaupd_
-#define PDSEUPD pdseupd_
-#endif
-
-extern "C" void SSAUPD(int *ido, char *bmat, int *n,
-                       char *which, int *nev, float *tol, float *resid,
-                       int *ncv, float *v, int *ldv,
-                       int *iparam, int *ipntr,
-                       float *workd, float *workl, int *lworkl, int *info);
-
-extern "C" void SSEUPD(int *, char *, int *, float *,
-                       float *, int *, float *, char *, int *, char *,
-                       int *, float *, float *, int *, float *,
-                       int *, int  *, int *, float *,
-                       float *, int *, int *);
-
-extern "C" void DSAUPD(int *ido, char *bmat, int *n,
-                       char *which, int *nev, double *tol, double *resid,
-                       int *ncv, double *v, int *ldv,
-                       int *iparam, int *ipntr,
-                       double *workd, double *workl, int *lworkl, int *info);
-
-extern "C" void DSEUPD(int *, char *, int *, double *,
-                       double *, int *, double *, char *, int *, char *,
-                       int *, double *, double *, int *, double *,
-                       int *, int  *, int *, double *,
-                       double *, int *, int *);
-
-#ifdef MFEM_USE_MPI
-
-extern "C" void PSSAUPD(int *comm, int *ido, char *bmat, int *n,
-                        char *which, int *nev, float *tol, float *resid,
-                        int *ncv, float *v, int *ldv,
-                        int *iparam, int *ipntr,
-                        float *workd, float *workl, int *lworkl, int *info);
-
-extern "C" void PSSEUPD(int *comm, int *, char *, int *, float *,
-                        float *, int *, float *, char *, int *, char *,
-                        int *, float *, float *, int *, float *,
-                        int *, int  *, int *, float *,
-                        float *, int *, int *);
-
-extern "C" void PDSAUPD(int *comm, int *ido, char *bmat, int *n,
-                        char *which, int *nev, double *tol, double *resid,
-                        int *ncv, double *v, int *ldv,
-                        int *iparam, int *ipntr,
-                        double *workd, double *workl, int *lworkl, int *info);
-
-extern "C" void PDSEUPD(int *comm, int *, char *, int *, double *,
-                        double *, int *, double *, char *, int *, char *,
-                        int *, double *, double *, int *, double *,
-                        int *, int  *, int *, double *,
-                        double *, int *, int *);
-
-#endif
-
-extern "C" {
-   void arpackgetcommdbg_(int *,int *,int *);
-   void arpacksetcommdbg_(int *,int *,int *);
-   void arpacksymdbg_(int *,int *,int *,int *,int *,int *,int *);
-   void arpacknonsymdbg_(int *,int *,int *,int *,int *,int *,int *);
-   void arpackcmplxdbg_(int *,int *,int *,int *,int *,int *,int *);
-}
-
-namespace mfem
-{
-
-/// Wrapper for the ARPACK routine SSAUPD or DSAUPD
-class ArPackSAUPD : public SymEigensolver, public SymGenEigensolver
-{
-public:
-
-   ArPackSAUPD();
-   virtual ~ArPackSAUPD();
-
-   /** ARPACK modes are described in section 3.5 of the ARPACK manual.
-       Mode 1: regular mode to solve A x = lambda x
-               No solver and no mass matrix are needed.
-       Mode 2: regular inverse mode to solve A x = lambda M x
-               Both A and M are needed and the solver should compute M^{-1}.
-       Mode 3: shift-invert mode to solve either A x = lambda x
-               or A x = lambda M x
-               Mass matrix is optional.  The solver should compute
-               (A-sigma I)^{-1} or (A-sigma M)^{-1}.  The shift parameter,
-               sigma, also needs to be set with SetShift().
-       Mode 4: Buckling mode to solve K x = lambda K_G x
-               K is set using SetMassMatrix(), K_G is set using SetOperator(),
-          and the solver should compute (K-sigma K_G)^{-1}.  The shift
-               parameter, sigma, also needs to be set with SetShift().
-       Mode 5: Cayley mode to solve A x = lambda M x
-               Both A and M are needed and the solver should compute
-               (A - sigma M)^{-1}.  The shift parameter, sigma, also needs
-               to be set with SetShift().
-    */
-   void SetMode(int mode);
-
-   inline void SetTol(real_t tol) override         {      tol_ = tol;      }
-   inline void SetMaxIter(int max_iter) override   { max_iter_ = max_iter; }
-   inline void SetPrintLevel(int logging) override {  logging_ = logging;  }
-   inline void SetShift(real_t sigma)              {    sigma_ = sigma;    }
-   inline void SetNumModes(int num_eigs) override  {      nev_ = num_eigs; }
-
-   virtual void SetSolver(Solver & solver);
-   virtual void SetOperator(const Operator & A) override;
-   virtual void SetMassMatrix(const Operator & M);
-   virtual void SetOperators(const Operator & A, const Operator & B) override
-   { SetOperator(A); SetMassMatrix(B); }
-
-   void Solve() override;
-
-   virtual int GetNumConverged() const override { return iparam_[4]; }
-
-   /// Collect the converged eigenvalues
-   virtual void GetEigenvalues(Array<real_t> & eigenvalues) const override;
-
-   /// Extract a single eigenvector
-   virtual const Vector & GetEigenvector(unsigned int i) const override;
-
-   /// Transfer ownership of the converged eigenvectors
-   Vector ** StealEigenvectors() override;
-
-protected:
-
-   int myid_;       // Index of this processor
-   int max_iter_;
-   int logging_;
-
-   // The following variables are for ARPACK
-   int nloc_;       // number of items stored locally
-   int nev_;        // number of requested eigenvalues
-   int ncv_;        // number of ritz vectors
-   int rvec_;       // boolean to return eigenvectors as well
-   int mode_;       // 1 = standard, 2 = generalized, 3 = shift invert,
-   // 4 = buckling, 5 = Cayley
-   int lworkl_;     // length of lworkl_ work array
-   int iparam_[12]; // arpack parameters
-   int ipntr_[12];  // arpack pointers
-
-   char bmat_;      // I for standard problem, G for generalized
-   char which_[3];  // spectrum portion: LA, SA, LM, SM, BE
-   char hwmny_;     // DSEUPD: A for all eigenvalues, S for some
-
-   real_t tol_;     // relative accuracy bound for Ritz values
-   real_t sigma_;   // eigenvalue shift parameter
-
-   int    * select_;// workspace used during eigenvalue computation
-   real_t * dv_;    // Ritz values
-   real_t * v_;     // ncv Lanczos basis vectors
-   real_t * resid_; // residual vector
-   real_t * workd_; // work array for 3 vectors used in Arnoldi iteration
-   real_t * workl_; // work array
-
-   // Operators and Vectors needed outside of ARPACK
-   Solver   * solver_;
-   const Operator * A_;
-   const Operator * B_;
-
-   Vector * w_;
-   Vector * x_;
-   Vector * y_;
-   Vector * z_;
-
-   mutable Vector ** eigenvectors_;
-
-   std::string solverName_;
-
-   void reverseComm();
-
-   int reverseCommMode1();
-   int reverseCommMode2();
-   int reverseCommMode3();
-   int reverseCommMode4();
-   int reverseCommMode5();
-
-   virtual void prepareEigenvectors() const;
-
-   void printErrors(const int & info, const int iparam[],
-                    const char & bmat, const int & n,
-                    const char which[],
-                    const int & nev, const int & ncv,
-                    const int & lworkl );
-
-private:
-
-   virtual int computeNlocf() { return nloc_; }
-   virtual int computeIter(int & ido);
-   virtual int computeEigs();
-
-};
-
-#ifdef MFEM_USE_MPI
-
-class ArPackPSAUPD : public ArPackSAUPD
-{
-public:
-   ArPackPSAUPD(MPI_Comm comm);
-   virtual ~ArPackPSAUPD() {}
-
-   void SetOperator(const Operator & A);
-   void SetMassMatrix(const Operator & M);
-
-   /// Collect the converged eigenvalues
-   void GetEigenvalues(Array<real_t> & eigenvalues) const;
-
-   /// Extract a single eigenvector
-   const Vector & GetEigenvector(unsigned int i) const;
-
-   /// Transfer ownership of the converged eigenvectors
-   // HypreParVector ** StealEigenvectors();
-   Vector ** StealEigenvectors();
-
-protected:
-
-   void prepareEigenvectors() const;
-
-private:
-
-   MPI_Comm comm_;
-   MPI_Fint commf_; // Fortran style MPI communicator
-   int numProcs_;   // Number of processors
-
-   mutable HYPRE_Int * part_; // parallel partitioning for eigenvectors
-
-   int computeNlocf();
-   int computeIter(int & ido);
-   int computeEigs();
-
-};
-
-#endif // MFEM_USE_MPI
-
-};
-
-#endif // MFEM_USE_ARPACK
-
-#endif // MFEM_ARPACK
@@ -1,20 +0,0 @@
-// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
-// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-// LICENSE and NOTICE for details. LLNL-CODE-806117.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability visit https://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the BSD-3 license. We welcome feedback and contributions, see file
-// CONTRIBUTING.md for details.
-
-#include "linalg.hpp"
-#include "eigensolvers.hpp"
-
-using namespace std;
-
-namespace mfem
-{
-
-};
@@ -1,396 +0,0 @@
-// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
-// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-// LICENSE and NOTICE for details. LLNL-CODE-806117.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability visit https://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the BSD-3 license. We welcome feedback and contributions, see file
-// CONTRIBUTING.md for details.
-
-#ifndef MFEM_EIGENSOLVERS
-#define MFEM_EIGENSOLVERS
-
-#include "vector.hpp"
-
-namespace mfem
-{
-
-/// Abstract Eigenequation
-/// Defines the operator of the linear eigenvalue equation
-///    A x_i = lambda_i x_i
-/// Where A is a real-valued operator, the lambda_i are the eigenvalues,
-/// and x_i are the eigenvectors.
-class Eigenequation
-{
-protected:
-   Eigenequation() = default;
-
-public:
-   virtual ~Eigenequation() = default;
-
-   /// @brief Set the operator A of the eigenvalue equation
-   virtual void SetOperator(const Operator & A) = 0;
-};
-
-/// Abstract Complex-valued Eigenequation
-/// Defines the operator of the linear eigenvalue equation
-///    A x_i = lambda_i x_i
-/// Where A is a complex-valued operator, the lambda_i are the eigenvalues,
-/// and x_i are the eigenvectors.
-class ComplexEigenequation
-{
-protected:
-   ComplexEigenequation() = default;
-
-public:
-   virtual ~ComplexEigenequation() = default;
-
-   /// @brief Set the real and imaginary parts of the operator A
-   virtual void SetOperator(const Operator & Ar, const Operator & Ai) = 0;
-};
-
-/// Abstract Generalized Eigenequation
-/// Defines the operator of the linear eigenvalue equation
-///    A x_i = lambda_i B x_i
-/// Where A and B are real-valued operators, the lambda_i are the eigenvalues,
-/// and x_i are the eigenvectors.
-class GenEigenequation
-{
-protected:
-   GenEigenequation() = default;
-
-public:
-   virtual ~GenEigenequation() = default;
-
-   /// @brief Set the operators A and B of the generalized eigenvalue equation
-   virtual void SetOperators(const Operator & A, const Operator & B) = 0;
-};
-
-/// Abstract Complex-valued Generalized Eigenequation
-/// Defines the operator of the linear eigenvalue equation
-///    A x_i = lambda_i B x_i
-/// Where A and B are complex-valued operators, the lambda_i are the
-/// eigenvalues, and x_i are the eigenvectors.
-class ComplexGenEigenequation
-{
-protected:
-   ComplexGenEigenequation() = default;
-
-public:
-   virtual ~ComplexGenEigenequation() = default;
-
-   /// @brief Set the real and imaginary parts of the operators A and B
-   virtual void SetOperators(const Operator & Ar, const Operator & Ai,
-                             const Operator & Br, const Operator & Bi) = 0;
-};
-
-/// Abstract Eigensolver
-/// Computes eigenvalue/eigenvector pairs for the linear system
-///    A x_i = lambda_i x_i
-/// Where the lambda_i are the eigenvalues and x_i are the eigenvectors.
-class EigensolverBase
-{
-protected:
-   EigensolverBase() = default;
-
-public:
-   virtual ~EigensolverBase() = default;
-
-   /// @brief Stopping criteria based on numerical tolerance
-   ///
-   /// @note This may be defined differently by different solvers.
-   virtual void SetTol(real_t tol) = 0;
-
-   /// @brief Stopping criteria based on number of iterations required to
-   /// reach convergence.
-   ///
-   /// @note This may also be defined differently in different solvers.
-   virtual void SetMaxIter(int max_iter) = 0;
-
-   /// @brief Controls the type and amount of information printed to
-   ///        standard output.
-   virtual void SetPrintLevel(int logging) = 0;
-
-   /// @brief Set the number of desired eigenmodes to compute
-   virtual void SetNumModes(int num_eigs) = 0;
-
-   /// @brief Get the number of converged eigenmodes
-   virtual int GetNumConverged() const = 0;
-
-   /// @brief Perform the eigenvalue solve
-   virtual void Solve() = 0;
-};
-
-/// Symmetric Eigensolver
-/// If A^T = A the linear system must have real-valued eigenvalues
-/// and eigenvectors.
-class SymEigensolver : public EigensolverBase, public Eigenequation
-{
-protected:
-   SymEigensolver() = default;
-
-public:
-   virtual ~SymEigensolver() = default;
-
-   /// @brief Collect the converged eigenvalues
-   ///
-   /// The length of the array should equal the number of converged eigenvalues.
-   virtual void GetEigenvalues(Array<real_t> & eigenvalues) const = 0;
-
-   /// @brief Extract a single eigenvector
-   ///
-   /// The index i should be in the range [0, numConverged). The
-   virtual const Vector & GetEigenvector(unsigned int i) const = 0;
-
-   /// @brief Transfer ownership of the converged eigenvectors
-   ///
-   /// The array should contain numConverged vectors.
-   virtual Vector ** StealEigenvectors() = 0;
-};
-
-/// Symmetric Generalized Eigensolver
-/// If A^T = A and M^T = M the linear system must have real-valued eigenvalues
-/// and eigenvectors.
-class SymGenEigensolver : public EigensolverBase, public GenEigenequation
-{
-protected:
-   SymGenEigensolver() = default;
-
-public:
-   virtual ~SymGenEigensolver() = default;
-
-   /// @brief Collect the converged eigenvalues
-   ///
-   /// The length of the array should equal the number of converged eigenvalues.
-   virtual void GetEigenvalues(Array<real_t> & eigenvalues) const = 0;
-
-   /// @brief Extract a single eigenvector
-   ///
-   /// The index i should be in the range [0, numConverged). The
-   virtual const Vector & GetEigenvector(unsigned int i) const = 0;
-
-   /// @brief Transfer ownership of the converged eigenvectors
-   ///
-   /// The array should contain numConverged vectors.
-   virtual Vector ** StealEigenvectors() = 0;
-};
-
-/// Hermetian Eigensolver
-/// If A^H = A the linear system must have real-valued eigenvalues
-/// but may have complex-valued eigenvectors.
-class HermEigensolver : public EigensolverBase, public ComplexEigenequation
-{
-protected:
-   HermEigensolver() = default;
-
-public:
-   virtual ~HermEigensolver() = default;
-
-   /// @brief Collect the converged eigenvalues
-   ///
-   /// The length of the array should equal the number of converged eigenvalues.
-   virtual void GetEigenvalues(Array<real_t> & eigenvalues) const = 0;
-
-   /// @brief Extract a single eigenvector
-   ///
-   /// The index i should be in the range [0, 2*numConverged). The
-   /// vectors corresponding to even indices are the real parts of the
-   /// converged eigenvectors and the odd indices correspond to the
-   /// imaginary parts.
-   virtual const Vector & GetEigenvector(unsigned int i) const = 0;
-
-   /// @brief Transfer ownership of the converged eigenvectors
-   ///
-   /// The array should contain 2*numConverged vectors with the even
-   /// indices corresponding to the real parts of the converged
-   /// eigenvectors and the odd indices corresponding to the imaginary
-   /// parts.
-   virtual Vector ** StealEigenvectors() = 0;
-};
-
-/// Hermetian Generalized Eigensolver
-/// If A^H = A and M^H = M the linear system must have real-valued eigenvalues
-/// but may have complex-valued eigenvectors.
-class HermGenEigensolver :
-   public EigensolverBase, public ComplexGenEigenequation
-{
-protected:
-   HermGenEigensolver() = default;
-
-public:
-   virtual ~HermGenEigensolver() = default;
-
-   /// @brief Collect the converged eigenvalues
-   ///
-   /// The length of the array should equal the number of converged eigenvalues.
-   virtual void GetEigenvalues(Array<real_t> & eigenvalues) const = 0;
-
-   /// @brief Extract a single eigenvector
-   ///
-   /// The index i should be in the range [0, 2*numConverged). The
-   /// vectors corresponding to even indices are the real parts of the
-   /// converged eigenvectors and the odd indices correspond to the
-   /// imaginary parts.
-   virtual const Vector & GetEigenvector(unsigned int i) const = 0;
-
-   /// @brief Transfer ownership of the converged eigenvectors
-   ///
-   /// The array should contain 2*numConverged vectors with the even
-   /// indices corresponding to the real parts of the converged
-   /// eigenvectors and the odd indices corresponding to the imaginary
-   /// parts.
-   virtual Vector ** StealEigenvectors() = 0;
-};
-
-/// Non-Symmetric Eigensolver
-/// For general real-valued operators A the linear system must have
-/// eigenvalues and eigenvectors which form complex conjugate pairs.
-class NonSymEigensolver : public EigensolverBase, public Eigenequation
-{
-protected:
-   NonSymEigensolver() = default;
-
-public:
-   virtual ~NonSymEigensolver() = default;
-
-   /// @brief Collect the converged eigenvalues
-   ///
-   /// The length of the array should be the number of converged
-   /// eigenvalues. The complex-valued eigenvalues can be constructed
-   /// as: lambda_{2*j} = eig[2*j]+i*eig[2*j+1] and
-   /// lambda_{2*j+1} = eig[2*j]-i*eig[2*j+1]
-   /// With j in the range [0, numConverged/2)
-   virtual void GetEigenvalues(Array<real_t> & eig) const = 0;
-
-   /// @brief Extract a single eigenvector
-   ///
-   /// The index i should be in the range [0, numConverged). The
-   /// vectors corresponding to even indices are the real parts of the
-   /// converged eigenvectors and the odd indices correspond to the
-   /// imaginary parts. If needed, the complex conjugate pairs of
-   /// eigenvectors can be constructed in the same manner described
-   /// for the eigenvalues.
-   virtual const Vector & GetEigenvector(unsigned int i) const = 0;
-
-   /// @brief Transfer ownership of the converged eigenvectors
-   ///
-   /// The array should contain numConverged vectors with the even
-   /// indices corresponding to the real parts of the converged
-   /// eigenvectors and the odd indices corresponding to the imaginary
-   /// parts.
-   virtual Vector ** StealEigenvectors() = 0;
-};
-
-/// Non-Symmetric Eigensolver
-/// For general real-valued operators A and M the linear system must have
-/// eigenvalues and eigenvectors which form complex conjugate pairs.
-class NonSymGenEigensolver : public EigensolverBase, public GenEigenequation
-{
-protected:
-   NonSymGenEigensolver() = default;
-
-public:
-   virtual ~NonSymGenEigensolver() = default;
-
-   /// @brief Collect the converged eigenvalues
-   ///
-   /// The length of the array should be the number of converged
-   /// eigenvalues. The complex-valued eigenvalues can be constructed
-   /// as: lambda_{2*j} = eig[2*j]+i*eig[2*j+1] and
-   /// lambda_{2*j+1} = eig[2*j]-i*eig[2*j+1]
-   /// With j in the range [0, numConverged/2)
-   virtual void GetEigenvalues(Array<real_t> & eig) const = 0;
-
-   /// @brief Extract a single eigenvector
-   ///
-   /// The index i should be in the range [0, numConverged). The
-   /// vectors corresponding to even indices are the real parts of the
-   /// converged eigenvectors and the odd indices correspond to the
-   /// imaginary parts. If needed, the complex conjugate pairs of
-   /// eigenvectors can be constructed in the same manner described
-   /// for the eigenvalues.
-   virtual const Vector & GetEigenvector(unsigned int i) const = 0;
-
-   /// @brief Transfer ownership of the converged eigenvectors
-   ///
-   /// The array should contain numConverged vectors with the even
-   /// indices corresponding to the real parts of the converged
-   /// eigenvectors and the odd indices corresponding to the imaginary
-   /// parts.
-   virtual Vector ** StealEigenvectors() = 0;
-};
-
-/// Complex Eigensolver
-/// Can have arbitrary complex-valued eigenvalues and eigenvectors
-class ComplexEigensolver : public EigensolverBase, public ComplexEigenequation
-{
-protected:
-   ComplexEigensolver() = default;
-
-public:
-   virtual ~ComplexEigensolver() = default;
-
-   /// @brief Collect the converged eigenvalues
-   ///
-   /// The length of the array should be twice the number of converged
-   /// eigenvalues. The complex-valued eigenvalues can be constructed
-   /// as: lambda_j = eig[2*j]+i*eig[2*j+1]
-   virtual void GetEigenvalues(Array<real_t> & eig) const = 0;
-
-   /// @brief Extract a single eigenvector
-   ///
-   /// The index i should be in the range [0, 2*numConverged). The
-   /// vectors corresponding to even indices are the real parts of the
-   /// converged eigenvectors and the odd indices correspond to the
-   /// imaginary parts.
-   virtual const Vector & GetEigenvector(unsigned int i) const = 0;
-
-   /// @brief Transfer ownership of the converged eigenvectors
-   ///
-   /// The array should contain 2*numConverged vectors with the even
-   /// indices corresponding to the real parts of the converged
-   /// eigenvectors and the odd indices corresponding to the imaginary
-   /// parts.
-   virtual Vector ** StealEigenvectors() = 0;
-};
-
-/// Complex Generalized Eigensolver
-/// Can have arbitrary complex-valued eigenvalues and eigenvectors
-class ComplexGenEigensolver :
-   public EigensolverBase, public ComplexGenEigenequation
-{
-protected:
-   ComplexGenEigensolver() = default;
-
-public:
-   virtual ~ComplexGenEigensolver() = default;
-
-   /// @brief Collect the converged eigenvalues
-   ///
-   /// The length of the array should be twice the number of converged
-   /// eigenvalues. The complex-valued eigenvalues can be constructed
-   /// as: lambda_j = eig[2*j]+i*eig[2*j+1]
-   virtual void GetEigenvalues(Array<real_t> & eig) const = 0;
-
-   /// @brief Extract a single eigenvector
-   ///
-   /// The index i should be in the range [0, 2*numConverged). The
-   /// vectors corresponding to even indices are the real parts of the
-   /// converged eigenvectors and the odd indices correspond to the
-   /// imaginary parts.
-   virtual const Vector & GetEigenvector(unsigned int i) const = 0;
-
-   /// @brief Transfer ownership of the converged eigenvectors
-   ///
-   /// The array should contain 2*numConverged vectors with the even
-   /// indices corresponding to the real parts of the converged
-   /// eigenvectors and the odd indices corresponding to the imaginary
-   /// parts.
-   virtual Vector ** StealEigenvectors() = 0;
-};
-
-}
-
-#endif
@@ -6556,7 +6556,7 @@ HypreLOBPCG::SetPreconditioner(Solver & precond)
 }

 void
-HypreLOBPCG::SetOperator(const Operator & A)
+HypreLOBPCG::SetOperator(Operator & A)
 {
   HYPRE_BigInt locSize = A.Width();

@@ -6603,7 +6603,7 @@ HypreLOBPCG::SetOperator(const Operator & A)
 }

 void
-HypreLOBPCG::SetMassMatrix(const Operator & M)
+HypreLOBPCG::SetMassMatrix(Operator & M)
 {
   matvec_fn.MatvecCreate  = this->OperatorMatvecCreate;
   matvec_fn.Matvec        = this->OperatorMatvec;
@@ -6624,7 +6624,7 @@ HypreLOBPCG::GetEigenvalues(Array<real_t> & eigs) const
   }
 }

-const Vector &
+const HypreParVector &
 HypreLOBPCG::GetEigenvector(unsigned int i) const
 {
   return multi_vec->GetVector(i);
@@ -6866,24 +6866,6 @@ HypreAME::SetPreconditioner(HypreSolver & precond)
   ams_precond = &precond;
 }

-void
-HypreAME::SetOperators(const Operator & opA, const Operator & opB)
-{
-   const HypreParMatrix * A = dynamic_cast<const HypreParMatrix *>(&opA);
-   if (A == NULL)
-   {
-      mfem_error("HypreAME::SetOperator : first operator not HypreParMatrix!");
-   }
-   SetOperator(*A);
-
-   const HypreParMatrix * B = dynamic_cast<const HypreParMatrix *>(&opB);
-   if (B == NULL)
-   {
-      mfem_error("HypreAME::SetOperator : second operator not HypreParMatrix!");
-   }
-   SetMassMatrix(*B);
-}
-
 void
 HypreAME::SetOperator(const HypreParMatrix & A)
 {
@@ -6942,7 +6924,7 @@ HypreAME::createDummyVectors() const
   }
 }

-const Vector &
+const HypreParVector &
 HypreAME::GetEigenvector(unsigned int i) const
 {
   if ( eigenvectors == NULL )
@@ -6953,7 +6935,7 @@ HypreAME::GetEigenvector(unsigned int i) const
   return *eigenvectors[i];
 }

-Vector **
+HypreParVector **
 HypreAME::StealEigenvectors()
 {
   if ( eigenvectors == NULL )
@@ -6962,11 +6944,11 @@ HypreAME::StealEigenvectors()
   }

   // Set the local pointers to NULL so that they won't be deleted later
-   Vector ** vecs = (Vector**)eigenvectors;
+   HypreParVector ** vecs = eigenvectors;
   eigenvectors = NULL;
   multi_vec = NULL;

-   return (Vector**)vecs;
+   return vecs;
 }

 }
@@ -18,9 +18,7 @@

 #include "../general/globals.hpp"
 #include "sparsemat.hpp"
-#include "eigensolvers.hpp"
 #include "hypre_parcsr.hpp"
-#include "eigensolvers.hpp"
 #include <mpi.h>

 // Enable internal hypre timing routines
@@ -2148,7 +2146,7 @@ public:
    A. Knyazev, M. Argentati, I. Lashuk, and E. Ovtchinnikov, SISC, 29(5),
    2224-2239, 2007.
 */
-class HypreLOBPCG : public SymGenEigensolver
+class HypreLOBPCG
 {
 private:
   MPI_Comm comm;
@@ -2238,43 +2236,38 @@ public:
   HypreLOBPCG(MPI_Comm comm);
   ~HypreLOBPCG();

-   void SetTol(real_t tol) override;
+   void SetTol(real_t tol);
   // not implemented in HYPRE
   // real_t GetTol() const;
   void SetRelTol(real_t rel_tol);
   // not implemented in HYPRE
   // real_t GetRelTol() const;
-   void SetMaxIter(int max_iter) override;
+   void SetMaxIter(int max_iter);
   // not implemented in HYPRE
   // int GetMaxIter() const;
-   void SetPrintLevel(int logging) override;
-   void SetNumModes(int num_eigs) override { nev = num_eigs; }
+   void SetPrintLevel(int logging);
+   void SetNumModes(int num_eigs) { nev = num_eigs; }
   void SetPrecondUsageMode(int pcg_mode);
   void SetRandomSeed(int s) { seed = s; }
   void SetInitialVectors(int num_vecs, HypreParVector ** vecs);

   // The following four methods support general operators
   void SetPreconditioner(Solver & precond);
-   void SetOperators(const Operator & A, const Operator & B) override
-   { SetOperator(A); SetMassMatrix(B); }
-   void SetOperator(const Operator & A);
-   void SetMassMatrix(const Operator & M);
+   void SetOperator(Operator & A);
+   void SetMassMatrix(Operator & M);
   void SetSubSpaceProjector(Operator & proj) { subSpaceProj = &proj; }

   /// Solve the eigenproblem
-   void Solve() override;
-
-   int GetNumConverged() const override { return nev; }
+   void Solve();

   /// Collect the converged eigenvalues
-   void GetEigenvalues(Array<real_t> & eigenvalues) const override;
+   void GetEigenvalues(Array<real_t> & eigenvalues) const;

   /// Extract a single eigenvector
-   const Vector & GetEigenvector(unsigned int i) const override;
+   const HypreParVector & GetEigenvector(unsigned int i) const;

   /// Transfer ownership of the converged eigenvectors
-   Vector ** StealEigenvectors() override
-   { return (Vector**)multi_vec->StealVectors(); }
+   HypreParVector ** StealEigenvectors() { return multi_vec->StealVectors(); }
 };

 /** AME eigenvalue solver in hypre
@@ -2299,7 +2292,7 @@ public:
    mass matrix but it seems unlikely that this would be useful so it is not the
    default behavior.
 */
-class HypreAME : public SymGenEigensolver
+class HypreAME
 {
 private:
   int myid;
@@ -2328,31 +2321,28 @@ public:
   HypreAME(MPI_Comm comm);
   ~HypreAME();

-   void SetTol(real_t tol) override;
+   void SetTol(real_t tol);
   void SetRelTol(real_t rel_tol);
-   void SetMaxIter(int max_iter) override;
-   void SetPrintLevel(int logging) override;
-   void SetNumModes(int num_eigs) override;
+   void SetMaxIter(int max_iter);
+   void SetPrintLevel(int logging);
+   void SetNumModes(int num_eigs);

   // The following four methods support operators of type HypreParMatrix.
   void SetPreconditioner(HypreSolver & precond);
-   void SetOperators(const Operator & opA, const Operator & opB) override;
   void SetOperator(const HypreParMatrix & A);
   void SetMassMatrix(const HypreParMatrix & M);

   /// Solve the eigenproblem
-   void Solve() override;
-
-   int GetNumConverged() const override { return nev; }
+   void Solve();

   /// Collect the converged eigenvalues
-   void GetEigenvalues(Array<real_t> & eigenvalues) const override;
+   void GetEigenvalues(Array<real_t> & eigenvalues) const;

   /// Extract a single eigenvector
-   const Vector & GetEigenvector(unsigned int i) const override;
+   const HypreParVector & GetEigenvector(unsigned int i) const;

   /// Transfer ownership of the converged eigenvectors
-   Vector ** StealEigenvectors() override;
+   HypreParVector ** StealEigenvectors();
 };

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

-#ifdef MFEM_USE_ARPACK
-#include "arpack.hpp"
-#endif
-
 #ifdef MFEM_USE_MKL_PARDISO
 #include "pardiso.hpp"
 #endif
@@ -844,22 +844,6 @@ public:
 };


-/// Zero Operator N: x -> 0.
-class ZeroOperator : public Operator
-{
-public:
-   /// Create an zero operator of size @a n.
-   explicit ZeroOperator(int n) : Operator(n) { }
-
-   /// Operator application
-   void Mult(const Vector &x, Vector &y) const override
-   { y.SetSize(width); y = 0_r; }
-
-   /// Application of the transpose
-   void MultTranspose(const Vector &x, Vector &y) const override
-   { y.SetSize(width); y = 0_r; }
-};
-
 /// Identity Operator I: x -> x.
 class IdentityOperator : public Operator
 {
@@ -4156,31 +4156,20 @@ void PetscNonlinearSolver::SetUpdate(void (*update)(Operator *,int,
 void PetscNonlinearSolver::Mult(const Vector &b, Vector &x) const
 {
   SNES snes = (SNES)obj;
-   MPI_Comm comm = PetscObjectComm(obj);

-   // Reduction needed: some processes may have null local size while others don't,
-   // and VecPlaceArray (used by PlaceMemory) is a logically collective operation.
-   PetscBool b_nonempty = b.Size() ? PETSC_TRUE : PETSC_FALSE;
-#if PETSC_VERSION_LT(3,24,0)
-   mpiierr = MPI_Allreduce(MPI_IN_PLACE,&b_nonempty,1,MPIU_BOOL,MPI_LOR,comm);
-#else
-   mpiierr = MPI_Allreduce(MPI_IN_PLACE,&b_nonempty,1,MPI_C_BOOL,MPI_LOR,comm);
-#endif
-   CCHKERRQ(comm,mpiierr);
-
-   // Always create B with allocate=false so that PlaceMemory can be called on
-   // it regardless of whether b was empty on a previous call.
-   if (!B) { B = new PetscParVector(comm, *this, true, false); }
-   if (!X) { X = new PetscParVector(comm, *this, false, false); }
+   bool b_nonempty = b.Size();
+   if (!B) { B = new PetscParVector(PetscObjectComm(obj), *this, true); }
+   if (!X) { X = new PetscParVector(PetscObjectComm(obj), *this, false, false); }
   X->PlaceMemory(x.GetMemory(),iterative_mode);
   if (b_nonempty) { B->PlaceMemory(b.GetMemory()); }
+   else { *B = 0.0; }

   Customize();

   if (!iterative_mode) { *X = 0.; }

-   // Solve the system. Pass nullptr for b when empty (PETSc treats it as zero RHS).
-   ierr = SNESSolve(snes, b_nonempty ? B->x : nullptr, X->x); PCHKERRQ(snes, ierr);
+   // Solve the system.
+   ierr = SNESSolve(snes, B->x, X->x); PCHKERRQ(snes, ierr);
   X->ResetMemory();
   if (b_nonempty) { B->ResetMemory(); }
 }
@@ -126,11 +126,11 @@ EXAMPLE_TEST_DIRS := examples
 MINIAPP_SUBDIRS = common electromagnetics meshing performance tools \
 toys nurbs gslib adjoint solvers shifted mtop parelag tribol autodiff dfem \
 hooke multidomain dpg hdiv-linear-solver spde diag-smoothers contact \
- fluids/navier fluids/schrodinger-flow plasma plasma/pic
+ fluids/navier fluids/schrodinger-flow plasma
 MINIAPP_DIRS := $(addprefix miniapps/,$(MINIAPP_SUBDIRS))
 MINIAPP_TEST_DIRS := $(filter-out %/common,$(MINIAPP_DIRS))
 MINIAPP_USE_COMMON := $(addprefix miniapps/,electromagnetics meshing tools \
- toys shifted dpg diag-smoothers fluids/navier plasma plasma/pic)
+ toys shifted dpg diag-smoothers fluids/navier plasma)

 EM_DIRS = $(EXAMPLE_DIRS) $(MINIAPP_DIRS)

@@ -302,7 +302,7 @@ endif
 MFEM_REQ_LIB_DEPS = SUPERLU MUMPS METIS FMS CONDUIT SIDRE LAPACK SUNDIALS\
 SUITESPARSE STRUMPACK GINKGO GNUTLS HDF5 NETCDF SLEPC PETSC MPFR PUMI HIOP\
 GSLIB OCCA CEED RAJA UMPIRE MKL_CPARDISO MKL_PARDISO AMGX MAGMA CALIPER PARELAG\
- TRIBOL BENCHMARK MOONOLITH ALGOIM ARPACK
+ TRIBOL BENCHMARK MOONOLITH ALGOIM


 PETSC_ERROR_MSG = $(if $(PETSC_FOUND),,. PETSC config not found: $(PETSC_VARS))
@@ -371,8 +371,7 @@ MFEM_DEFINES = MFEM_VERSION MFEM_VERSION_STRING MFEM_GIT_STRING MFEM_USE_MPI\
 MFEM_USE_SIMD MFEM_USE_ADIOS2 MFEM_USE_MKL_CPARDISO MFEM_USE_MKL_PARDISO MFEM_USE_AMGX\
 MFEM_USE_MAGMA MFEM_USE_MUMPS MFEM_USE_ADFORWARD MFEM_USE_CODIPACK MFEM_USE_CALIPER\
 MFEM_USE_BENCHMARK MFEM_USE_PARELAG MFEM_USE_TRIBOL MFEM_USE_ALGOIM MFEM_USE_ENZYME\
- MFEM_SOURCE_DIR MFEM_INSTALL_DIR MFEM_SHARED_BUILD MFEM_USE_DOUBLE MFEM_USE_SINGLE\
-MFEM_USE_ARPACK
+ MFEM_SOURCE_DIR MFEM_INSTALL_DIR MFEM_SHARED_BUILD MFEM_USE_DOUBLE MFEM_USE_SINGLE

 # List of makefile variables that will be written to config.mk:
 MFEM_CONFIG_VARS = MFEM_CXX MFEM_HOST_CXX MFEM_CPPFLAGS MFEM_CXXFLAGS\
@@ -734,7 +733,6 @@ status info:
 	$(info MFEM_TIMER_TYPE        = $(MFEM_TIMER_TYPE))
 	$(info MFEM_USE_SUNDIALS      = $(MFEM_USE_SUNDIALS))
 	$(info MFEM_USE_SUITESPARSE   = $(MFEM_USE_SUITESPARSE))
-	$(info MFEM_USE_ARPACK        = $(MFEM_USE_ARPACK))
 	$(info MFEM_USE_SUPERLU       = $(MFEM_USE_SUPERLU))
 	$(info MFEM_USE_SUPERLU5      = $(MFEM_USE_SUPERLU5))
 	$(info MFEM_USE_MUMPS         = $(MFEM_USE_MUMPS))
@@ -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++)
   {
@@ -3206,22 +3206,10 @@ public:


 /// Extrude a 1D mesh
-/**
- * @param mesh      1D mesh
- * @param ny        number of transverse elements of the extruded mesh
- * @param sy        physical size in the direction of extrusion
- * @param closed    if false, only the original boundaries are extruded,
- *                  otherwise boundaries are generated all around the domain
- */
 Mesh *Extrude1D(Mesh *mesh, const int ny, const real_t sy,
                const bool closed = false);

 /// Extrude a 2D mesh
-/**
- * @param mesh      2D mesh
- * @param nz        number of transverse elements of the extruded mesh
- * @param sz        physical size in the direction of extrusion
- */
 Mesh *Extrude2D(Mesh *mesh, const int nz, const real_t sz);

 /** @brief Constructs the smallest possible [0,1]^dim serial mesh that can be
@@ -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;
@@ -1516,15 +1516,12 @@ void Mesh::ReadInlineMesh(std::istream &input, bool generate_edges)
 void Mesh::ReadGmshMesh(std::istream &input, int &curved, int &read_gf)
 {
   string buff;
-   string version;
+   real_t version;
   int binary, dsize;
   input >> version >> binary >> dsize;
-   if (version != "2.2")
+   if (version < 2.2)
   {
-      MFEM_ABORT("Gmsh file version must be 2.2, found version "
-                 << version << ".\n"
-                 "To convert your mesh to the required format, use:\n"
-                 "  gmsh -format msh22 -save -o output.msh input.msh");
+      MFEM_ABORT("Gmsh file version < 2.2");
   }
   if (dsize != sizeof(double))
   {
@@ -5639,12 +5639,6 @@ Mesh ParMesh::GetSerialMesh(int save_rank) const
      }
   }

-   if (MyRank == save_rank)
-   {
-      attribute_sets.Copy(serialmesh.attribute_sets);
-      bdr_attribute_sets.Copy(serialmesh.bdr_attribute_sets);
-   }
-
   MPI_Barrier(MyComm);
   return serialmesh;
 }
@@ -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);
@@ -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

@@ -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]; }
 }
@@ -329,8 +329,8 @@ int main(int argc, char *argv[])

      for (int i=0; i<nev; i++)
      {
-         // convert eigenvector from Vector to ParGridFunction
-         x.Distribute(lobpcg->GetEigenvector(i));
+         // convert eigenvector from HypreParVector to ParGridFunction
+         x = lobpcg->GetEigenvector(i);

         mode_name << "mode_" << setfill('0') << setw(2) << i << "."
                   << setfill('0') << setw(6) << myid;
@@ -357,8 +357,8 @@ int main(int argc, char *argv[])
                 << ", Lambda = " << eigenvalues[i] << endl;
         }

-         // convert eigenvector from Vector to ParGridFunction
-         x.Distribute(lobpcg->GetEigenvector(i));
+         // convert eigenvector from HypreParVector to ParGridFunction
+         x = lobpcg->GetEigenvector(i);

         mode_sock << "parallel " << num_procs << " " << myid << "\n"
                   << "solution\n" << *pmesh << x << flush
@@ -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,8 +29,6 @@ else
   MINIAPPS = $(PAR_MINIAPPS) $(SEQ_MINIAPPS)
 endif

-PLASMA_SUBDIRS = pic
-
 .SUFFIXES:
 .SUFFIXES: .o .cpp .mk
 .PHONY: all lib-common clean clean-build clean-exec
@@ -46,12 +47,7 @@ COMMON_O=
 %: %.cpp
 %.o: %.cpp

-all: $(MINIAPPS) subdirs
-
-.PHONY: subdirs $(PLASMA_SUBDIRS)
-subdirs: $(PLASMA_SUBDIRS)
-$(PLASMA_SUBDIRS): lib-common
-	$(MAKE) -C $(BLD)$(@)
+all: $(MINIAPPS)

 # Rules for building the miniapps
 %: $(SRC)%.cpp $(COMMON_O) $(MFEM_LIB_FILE) $(CONFIG_MK) | lib-common
@@ -79,15 +75,11 @@ RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
 $(MFEM_LIB_FILE):
 	$(error The MFEM library is not built)

-ALL_CLEAN_SUBDIRS = $(addsuffix /clean,$(PLASMA_SUBDIRS))
-.PHONY: $(ALL_CLEAN_SUBDIRS)
-$(ALL_CLEAN_SUBDIRS):
-	$(MAKE) -C $(BLD)$(@D) $(@F)
-
 clean: clean-build clean-exec

-clean-build: $(addsuffix /clean,$(PLASMA_SUBDIRS))
+clean-build:
 	rm -f *.o *~ $(SEQ_MINIAPPS) $(PAR_MINIAPPS)
 	rm -rf *.dSYM *.TVD.*breakpoints

 clean-exec:
+
@@ -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,85 +0,0 @@
-# Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
-# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-# LICENSE and NOTICE for details. LLNL-CODE-806117.
-#
-# This file is part of the MFEM library. For more information and source code
-# availability visit https://mfem.org.
-#
-# MFEM is free software; you can redistribute it and/or modify it under the
-# terms of the BSD-3 license. We welcome feedback and contributions, see file
-# CONTRIBUTING.md for details.
-
-# Use the MFEM build directory
-MFEM_DIR ?= ../../..
-MFEM_BUILD_DIR ?= ../../..
-MFEM_INSTALL_DIR ?= ../../../mfem
-SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/plasma/pic/,)
-CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
-   $(wildcard $(MFEM_INSTALL_DIR)/share/mfem/config.mk))
-
-MFEM_LIB_FILE = mfem_is_not_built
-include $(CONFIG_MK)
-
-PAR_MINIAPPS =
-
-ifeq ($(MFEM_USE_GSLIB),YES)
-   PAR_MINIAPPS += electrostatic-pic
-endif
-
-ifeq ($(MFEM_USE_MPI),NO)
-   MINIAPPS =
-else
-   MINIAPPS = $(PAR_MINIAPPS)
-endif
-
-.SUFFIXES:
-.SUFFIXES: .o .cpp .mk
-.PHONY: all lib-common clean clean-build clean-exec
-.PRECIOUS: %.o
-
-COMMON_LIB = -L$(MFEM_BUILD_DIR)/miniapps/common -lmfem-common
-
-# If MFEM_SHARED is set, add the ../common rpath
-COMMON_LIB += $(if $(MFEM_SHARED:YES=),,\
-   $(MFEM_XLINKER)-rpath,$(abspath $(MFEM_BUILD_DIR)/miniapps/common))
-
-# Remove built-in rules
-%: %.cpp
-%.o: %.cpp
-
-all: $(MINIAPPS)
-
-# Rules for building the miniapps
-electrostatic-pic: electrostatic-pic.cpp $(MFEM_LIB_FILE) $(CONFIG_MK) | lib-common
-	$(MFEM_CXX) $(MFEM_FLAGS) -c $<
-	$(MFEM_CXX) $(MFEM_LINK_FLAGS) -o $@ $@.o $(COMMON_LIB) $(MFEM_LIBS)
-
-# Rule for building lib-common
-lib-common:
-	$(MAKE) -C $(MFEM_BUILD_DIR)/miniapps/common
-
-
-MFEM_TESTS = MINIAPPS
-include $(MFEM_TEST_MK)
-
-# Testing: "test" target and mfem-test* variables are defined in config/test.mk
-
-# Testing: Specific execution options
-RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
-electrostatic-pic-test-par: electrostatic-pic
-	@$(call mfem-test,$<, $(RUN_MPI), PIC miniapp,\
-		-rdi 2 -npt 40960 -k 0.2855993321 -a 0.05 -nt 200 -nx 16 -ny 16\
-		-O 1 -q 0.01181640625 -m 0.01181640625 -oci 1000 -dt 0.1)
-
-# Generate an error message if the MFEM library is not built and exit
-$(MFEM_LIB_FILE):
-	$(error The MFEM library is not built)
-
-clean: clean-build clean-exec
-
-clean-build:
-	rm -f *.o *~ $(SEQ_MINIAPPS) $(PAR_MINIAPPS)
-	rm -rf *.dSYM *.TVD.*breakpoints
-
-clean-exec:
-	@rm -rf electrostatic-pic_* *.csv energy.csv
@@ -39,6 +39,7 @@ set(UNIT_TESTS_SRCS
  dfem/test_divergence.cpp
  dfem/test_lvector_interface.cpp
  dfem/test_mass.cpp
+  dfem/test_univarsolver.cpp
  general/test_array.cpp
  general/test_scan.cpp
  general/test_arrays_by_name.cpp
@@ -71,12 +72,10 @@ set(UNIT_TESTS_SRCS
  linalg/test_ode2.cpp
  linalg/test_operator.cpp
  linalg/test_particlevector.cpp
-  linalg/test_petsc_nonlinear.cpp
  linalg/test_sparsesmoothers.cpp
  linalg/test_vector.cpp
  mesh/mesh_test_utils.cpp
  mesh/test_exodus_reader.cpp
-  mesh/test_mfem_mesh_reader.cpp
  mesh/test_exodus_writer.cpp
  mesh/test_face_orientations.cpp
  mesh/test_fms.cpp
@@ -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
@@ -0,0 +1,359 @@
+// 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 <algorithm>
+#include <cmath>
+
+#include "mfem.hpp"
+#include "unit_tests.hpp"
+
+#ifdef MFEM_USE_ENZYME
+
+using mfem::real_t;
+using namespace mfem::future;
+
+MFEM_HOST_DEVICE inline real_t FlowResistance(real_t eqps, real_t sigma_y, real_t n, real_t ep_0)
+{
+    return sigma_y*(1.0 + std::pow((eqps)/ep_0, n));
+}
+
+using J2PlasticityParameters = tuple<real_t, real_t, real_t, real_t, real_t, real_t>;
+
+// Residual function that is solved in the plasticity model.
+// Made a free function to facilitate Enzyme differentiation.
+real_t J2PlasticityResidual(real_t delta_eqps, J2PlasticityParameters p)
+{
+    auto [eqps, q, G, sigma_y, n, ep_0] = p;
+    return q - 3.0*G*delta_eqps - FlowResistance(eqps + delta_eqps, sigma_y, n, ep_0);
+}
+
+struct J2Plasticity {
+  static constexpr int dim = 3;  ///< spatial dimension
+  static constexpr int N_INTERNAL_STATES = 10;
+  static constexpr real_t tol = 1e-10;  ///< relative tolerance on residual mag to judge convergence of return map
+
+  real_t E;        ///< Young's modulus
+  real_t nu;       ///< Poisson's ratio
+  real_t sigma_y;  ///< Yield strength
+  real_t n;        ///< Hardening index
+  real_t ep_0;     ///< Reference plastic strain
+  
+  /// @brief variables required to characterize the hysteresis response
+  struct InternalState {
+    tensor<real_t, dim, dim> plastic_strain;
+    real_t accumulated_plastic_strain;
+  };
+
+  /// Internal state variables in a flattened array for storing in a global field
+  using PackedInternalState = mfem::future::tensor<real_t, N_INTERNAL_STATES>;
+
+  // Unflatten internal state variables
+  MFEM_HOST_DEVICE static inline InternalState unpack_internal_state(
+      const mfem::future::tensor<real_t, N_INTERNAL_STATES>& packed_state)
+  {
+    auto plastic_strain =
+        mfem::future::make_tensor<dim, dim>([&packed_state](int i, int j) { return packed_state[dim * i + j]; });
+    real_t accumulated_plastic_strain = packed_state[N_INTERNAL_STATES - 1];
+    return {plastic_strain, accumulated_plastic_strain};
+  }
+
+  // Flatten internal state variables (for repacking into global field)
+  MFEM_HOST_DEVICE static inline PackedInternalState pack_internal_state(
+      const mfem::future::tensor<real_t, dim, dim>& plastic_strain, real_t accumulated_plastic_strain)
+  {
+    PackedInternalState packed_state{};
+    for (int i = 0, ij = 0; i < dim; i++) {
+      for (int j = 0; j < dim; j++, ij++) {
+        packed_state[ij] = plastic_strain[i][j];
+      }
+    }
+    packed_state[N_INTERNAL_STATES - 1] = accumulated_plastic_strain;
+    return packed_state;
+  }
+
+  // Compute the new stress and the internal state variables
+  MFEM_HOST_DEVICE inline tuple<tensor<real_t, dim, dim>, PackedInternalState>
+  update(tensor<real_t, dim, dim> dudxi,
+         PackedInternalState internal_state,
+         tensor<real_t, dim, dim> J,
+         real_t w) const
+  {
+    auto invJ = inv(J);
+    const auto dudX = dudxi * invJ;
+    auto I = IdentityMatrix<dim>();
+    const real_t K = E / (3.0 * (1.0 - 2.0 * nu));
+    const real_t G = 0.5 * E / (1.0 + nu);
+
+    auto [plastic_strain, accumulated_plastic_strain] = unpack_internal_state(internal_state);
+    
+    auto el_strain = sym(dudX) - plastic_strain;
+    auto p = K * tr(el_strain);
+    auto s = 2.0 * G * dev(el_strain);
+    auto q = std::sqrt(1.5) * norm(s);
+    real_t denom = q > 0.0? q : 1.0;
+    auto Np = 1.5 * s / denom;
+
+    if (q > FlowResistance(accumulated_plastic_strain, sigma_y, n, ep_0)) {
+        real_t lb = 0.0;
+        real_t ub = (q - FlowResistance(accumulated_plastic_strain, sigma_y, n, ep_0))/(3*G);
+        SolverSettings settings{1e-10*sigma_y, 1e-10, {lb, ub}};
+        // Use the differentiable univariate root finder.
+        // This has custom derivatives, so it's ok to differentiate this enclosing function.
+        real_t delta_eqps = SolveNewtonBisection<J2PlasticityResidual>(
+            0.5*(lb + ub), make_tuple(accumulated_plastic_strain, q, G, sigma_y, n, ep_0), settings);
+        accumulated_plastic_strain += delta_eqps;
+        plastic_strain += delta_eqps * Np;
+        s -= 2.0 * G * delta_eqps * Np;
+    }
+    auto Q_new = pack_internal_state(plastic_strain, accumulated_plastic_strain);
+    auto stress = s + p * I;
+    const real_t dV = det(J)*w;
+    return {stress*transpose(invJ)*dV, Q_new};
+  }
+  
+  MFEM_HOST_DEVICE inline tensor<real_t, dim, dim>
+  stress(tensor<real_t, dim, dim> dudxi,
+         PackedInternalState internal_state,
+         tensor<real_t, dim, dim> J,
+         real_t w) const
+  {
+    auto [stress, internal_state_new] = update(dudxi, internal_state, J, w);
+    return stress;
+  }
+
+  MFEM_HOST_DEVICE inline PackedInternalState
+  internal_state_new(tensor<real_t, dim, dim> dudxi,
+                     PackedInternalState internal_state,
+                     tensor<real_t, dim, dim> J,
+                     real_t w) const
+  {
+    auto [stress, internal_state_new] = update(dudxi, internal_state, J, w);
+    return internal_state_new;
+  }
+};
+
+
+// Register the custom derivatives for the solver.
+// This needs to be done for every residual function that the solver is applied on,
+// since the SolveNewtonBisection_impl is a function template, and we need a real
+// function with an address to specify the custom derivative.
+
+// Forward mode
+__attribute__((used))
+void *  __enzyme_register_derivative_newton_bisection_on_j2[2] = {
+  (void*) mfem::internal::SolveNewtonBisection_impl<J2PlasticityResidual, J2PlasticityParameters>,
+  (void*) mfem::internal::SolveNewtonBisection_impl_fwddiff<J2PlasticityResidual, J2PlasticityParameters>
+};
+
+ // Reverse mode
+__attribute__((used))
+void* __enzyme_register_gradient_SolveNewtonBisectionJ2[3] = {
+  (void*)mfem::internal::SolveNewtonBisection_impl<J2PlasticityResidual, J2PlasticityParameters>,
+  (void*)mfem::internal::SolveNewtonBisection_impl_aug<J2PlasticityResidual, J2PlasticityParameters>,
+  (void*)mfem::internal::SolveNewtonBisection_impl_rev<J2PlasticityResidual, J2PlasticityParameters>
+};
+
+// Create free functions for Enzyme to differentiate in the tests
+// Return by value version
+tensor<real_t, 3, 3> ComputeStress(
+    J2Plasticity* material, tensor<real_t, 3, 3> dudxi, 
+    J2Plasticity::PackedInternalState Q, tensor<real_t, 3, 3> J, real_t w)
+{
+    return material->stress(dudxi, Q, J, w);
+}
+
+// Return by reference version
+void ComputeStressRef(const J2Plasticity* material, const tensor<real_t, 3, 3>& dudxi,
+                      const J2Plasticity::PackedInternalState& Q,
+                      const tensor<real_t, 3, 3>& J, real_t w,
+                      tensor<real_t, 3, 3>& sigma)
+{
+  sigma = material->stress(dudxi, Q, J, w);
+}
+
+template <int dim>
+real_t elementwise_max_norm(tensor<real_t, dim, dim> A) {
+  real_t maxval = 0;
+  for (int i = 0; i < dim; i++) {
+    for (int j = 0; j < dim; j++) {
+      maxval = std::max(std::abs(A[i][j]), maxval);
+    }
+  }
+  return maxval;
+}
+
+
+TEST_CASE("Univariate function solver in a qfunction", "[univar]")
+{
+    J2Plasticity material{70.0e3, 0.34, 240.0, 0.15, 1e-3};
+    tensor<real_t, 3, 3> H{{{0.947667  , 0.9785799 , 0.33229148},
+                            {0.46866846, 0.5698887 , 0.16550303},
+                            {0.3101946 , 0.68948054, 0.74676657}}};
+    J2Plasticity::PackedInternalState Q{};
+    const tensor<real_t, 3, 3> J = IdentityMatrix<3>();
+    const real_t w = 1.0;
+
+    SECTION("Correctness")
+    {
+      // Checks that stress after update is on the yield surface.
+      auto [stress, Q_new] = material.update(H, Q, IdentityMatrix<3>(), 1.0);
+      real_t mises = std::sqrt(1.5)*norm(dev(stress));
+      real_t eqps = Q_new[9];
+      // This test only makes sense if the displacement gradient is big enough to 
+      // cuase yielding.
+      REQUIRE(eqps > 1e-9);
+      real_t Y = FlowResistance(eqps, material.sigma_y, material.n, material.ep_0);
+      CHECK(mises == MFEM_Approx(Y, 0.0, 1e-8));
+    }
+
+    SECTION("JVP")
+    {
+      // Compare forward mode derivative to finite difference approximation
+
+      tensor<real_t, 3, 3> H_dot{{{1.0, 0.0 , 0.0},
+                                  {0.0, 0.0 , 0.0},
+                                  {0.0, 0.0 , 0.0}}};
+
+      // Enzyme directional derivative (uses custom derivative of solver)
+      auto sigma_dot = __enzyme_fwddiff<tensor<real_t, 3, 3>>((void*)ComputeStress, 
+                                                              enzyme_const, &material,
+                                                              enzyme_dup, H, H_dot,
+                                                              enzyme_const, Q,
+                                                              enzyme_const, J,
+                                                              enzyme_const, w);
+      // sigma_dot = ∂sigma / ∂H[0, 0]
+      REQUIRE(sigma_dot[0][0] > 0.0);
+
+      // Finite difference derivative approximation
+      constexpr int dim = 3;
+      real_t eps = 1e-5;
+      tensor<real_t, 3, 3> sigma = ComputeStress(&material, H, Q, J, w);
+      tensor<real_t, 3, 3> sigma_p = ComputeStress(&material, H + eps*H_dot, Q, J, w);
+      tensor<real_t, 3, 3> sigma_dot_h = (1.0/eps)*(sigma_p - sigma);
+
+      tensor<real_t, 3, 3> rel_error = sigma_dot - sigma_dot_h;
+      for (int i = 0; i < dim; i++) {
+        for (int j = 0; j < dim; j++) {
+          real_t denom = sigma[i][j] != 0? sigma[i][j] : 1.0;
+          rel_error[i][j] /= denom;
+        }
+      }
+
+      CHECK(elementwise_max_norm(rel_error) < 1e-5);
+    }
+
+    SECTION("VJP")
+    {
+      // compare reverse mode derivative to finite differences
+
+      tensor<real_t, 3, 3> sigma;
+      ComputeStressRef(&material, H, Q, J, w, sigma);
+      double epsilon = 1e-6;
+      tensor<real_t, 3, 3> dH{{{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}}};
+      auto H_p = H + epsilon*dH;
+      tensor<real_t, 3, 3> sigma_p;
+      ComputeStressRef(&material, H_p, Q, J, w, sigma_p);
+      auto sigma_dot_h = (sigma_p - sigma)/epsilon;
+      // Note: sigma_dot_h[i,j] = ∂sigma[i,j]/∂H[0,0]
+
+      tensor<real_t, 3, 3> sigma_bar{{{1.0, 0.0, 0.0},
+                                      {0.0, 0.0, 0.0},
+                                      {0.0, 0.0, 0.0}}};
+
+      tensor<real_t, 3, 3> H_bar{};
+      J2Plasticity::PackedInternalState Q_bar{};
+      tensor<real_t, 3, 3> J_bar{};
+      __enzyme_autodiff<void>(
+        (void*)ComputeStressRef,  enzyme_const, &material, enzyme_dup, &H, &H_bar,
+        enzyme_dup, &Q, &Q_bar, enzyme_dup, &J, &J_bar, enzyme_const, w,
+        enzyme_dup, &sigma, &sigma_bar);
+
+      // H_bar[ij] = ∂sigma[0,0]/∂H[i,j]
+      // For this model, we expect the major symmetries in the tangent operator.
+      // Hence H_bar \approx sigma_dot_h
+
+      const double abs_tol = 1e-12;
+      const double rel_tol = 5e-6;
+      for (int i = 0; i < 3; i++) {
+        for (int j = 0; j < 3; j++) {
+          CHECK(H_bar[i][j] == MFEM_Approx(sigma_dot_h[i][j], abs_tol, rel_tol));
+        }
+      }
+    }
+}
+
+
+real_t nthroot_res(real_t x, tuple<real_t, real_t> p)
+{
+    auto [index, radicand] = p;
+    return std::pow(x, index) - radicand;
+}
+
+__attribute__((used))
+void* __enzyme_register_gradient_solver[3] = {
+  (void*)mfem::internal::SolveNewtonBisection_impl<nthroot_res, tuple<real_t, real_t>>,
+  (void*)mfem::internal::SolveNewtonBisection_impl_aug<nthroot_res, tuple<real_t, real_t>>,
+  (void*)mfem::internal::SolveNewtonBisection_impl_rev<nthroot_res, tuple<real_t, real_t>>
+};
+
+TEST_CASE("Univariate solver reverse mode derivative", "[univar]")
+{
+  auto mysqrt = [](real_t x) -> real_t
+    {
+      real_t x0 = x;
+      real_t index = 2.0;
+      real_t ub = std::max(1.0, x);
+      SolverSettings settings{1e-12, 1e-12, {0, ub}};
+      return SolveNewtonBisection<nthroot_res>(x0, make_tuple(index, x), settings);
+    };
+
+  real_t x = 2.0;
+  real_t dydx = __enzyme_autodiff<real_t>((void*)+mysqrt, enzyme_out, x);
+  CHECK(dydx == MFEM_Approx(0.5/std::sqrt(2.0)));
+}
+
+TEST_CASE("Univariate function solver robustness", "[univar]")
+{
+  SolverSettings settings{1e-12, 1e-12};
+
+  SECTION("Simple case")
+  {
+      auto Nthroot = [&settings](real_t x, real_t n) {
+          real_t x0 = std::max(x, 1.0);
+          settings.bounds = {0.0, x0};
+          return SolveNewtonBisection<nthroot_res>(x0, make_tuple(n, x), settings);
+      };
+      real_t x = 8.0;
+      real_t y = Nthroot(x, 3.0);
+      CHECK(y == MFEM_Approx(2.0));
+  }
+
+    SECTION("Stiff problem")
+    {
+        auto f = [](real_t x, real_t p) { return std::pow(x, p) - 1.0; };
+        real_t x0 = 0.1;
+        real_t p = 50;
+        settings.bounds = {0.0, 5.1};
+        real_t x = SolveNewtonBisection<+f>(x0, p, settings);
+        CHECK(x == MFEM_Approx(1.0));
+    }
+
+    SECTION("Works where Newton diverges")
+    {
+        auto f = [](double x, int) { return std::atan(x); };
+        real_t x0 = 1.5;
+        settings.bounds = {0.0, 2.0};
+        real_t x = SolveNewtonBisection<+f>(x0, int{}, settings);
+        CHECK(std::abs(x) == MFEM_Approx(0.0));
+    }
+}
+
+#endif // MFEM_USE_ENZYME
@@ -105,39 +105,39 @@ TEST_CASE("Integration rule order initialization", "[IntegrationRules]")
   SECTION("Segment rule constructed by accessing square rule")
   {
      auto &quad5_ir = intrules.Get(Geometry::SQUARE, 5);
-      REQUIRE(quad5_ir.GetOrder() >= 5);
+      REQUIRE(quad5_ir.GetOrder() == 5);
      // The segment integration rule of order 5 is lazy constructed when we get
      // the square integration rule of order 5. Make sure its order was
      // properly set:
      auto &line5_ir = intrules.Get(Geometry::SEGMENT, 5);
-      REQUIRE(line5_ir.GetOrder() >= 5);
+      REQUIRE(line5_ir.GetOrder() == 5);
   }

   SECTION("Segment rule constructed by accessing cube rule")
   {
      auto &hex7_ir = intrules.Get(Geometry::CUBE, 7);
-      REQUIRE(hex7_ir.GetOrder() >= 7);
+      REQUIRE(hex7_ir.GetOrder() == 7);
      // The segment integration rule of order 7 is lazy constructed when we get
      // the cube integration rule of order 7. Make sure its order was properly
      // set:
      auto &line7_ir = intrules.Get(Geometry::SEGMENT, 7);
-      REQUIRE(line7_ir.GetOrder() >= 7);
+      REQUIRE(line7_ir.GetOrder() == 7);
   }

   SECTION("Segment and triangle rules constructed by accessing prism rule")
   {
      auto &prism3_ir = intrules.Get(Geometry::PRISM, 3);
-      REQUIRE(prism3_ir.GetOrder() >= 3);
+      REQUIRE(prism3_ir.GetOrder() == 3);
      // The segment integration rule of order 3 is lazy constructed when we get
      // the prism integration rule of order 3. Make sure its order was properly
      // set:
      auto &line3_ir = intrules.Get(Geometry::SEGMENT, 3);
-      REQUIRE(line3_ir.GetOrder() >= 3);
+      REQUIRE(line3_ir.GetOrder() == 3);
      // The triangle integration rule of order 3 is lazy constructed when we
      // get the prism integration rule of order 3. Make sure its order was
      // properly set:
      auto &tri3_ir = intrules.Get(Geometry::TRIANGLE, 3);
-      REQUIRE(tri3_ir.GetOrder() >= 3);
+      REQUIRE(tri3_ir.GetOrder() == 3);
   }
 }

@@ -271,43 +271,3 @@ TEST_CASE("Simplex integration rules", "[SimplexRules]")
      }
   }
 }
-
-
-// Monomial exactness is tested by [SimplexRules] above, which now uses
-// positive-weight rules by default. The tests below verify properties
-// specific to the positive-weight rules: weight positivity, stability,
-// and interior point placement.
-
-TEST_CASE("Simplex rule positivity", "[IntegrationRules]")
-{
-   IntegrationRules rules;
-
-   SECTION("triangle rules have all positive weights for orders 0-25")
-   {
-      for (int order = 0; order <= 25; order++)
-      {
-         const IntegrationRule &ir = rules.Get(Geometry::TRIANGLE, order);
-         for (int i = 0; i < ir.GetNPoints(); i++)
-         {
-            INFO("order=" << order << ", point=" << i);
-            REQUIRE(ir.IntPoint(i).weight > 0.0);
-         }
-      }
-   }
-
-   SECTION("tet rules have all positive weights for orders 0-20")
-   {
-      for (int order = 0; order <= 20; order++)
-      {
-         const IntegrationRule &ir =
-            rules.Get(Geometry::TETRAHEDRON, order);
-         for (int i = 0; i < ir.GetNPoints(); i++)
-         {
-            INFO("order=" << order << ", point=" << i);
-            REQUIRE(ir.IntPoint(i).weight > 0.0);
-         }
-      }
-   }
-}
-
-
@@ -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++;
         }
@@ -25,201 +25,15 @@ void Func_3D_lin(const Vector &x, Vector &v)
   v[2] = -2.572 * x[0] + 1.321 * x[1] + 3.234 * x[2];
 }

-TEST_CASE("3D ProjectBdrCoefficientNormal Vector",
-          "[GridFunction]"
-          "[VectorGridFunctionCoefficient]")
-{
-   const int n = 1;
-   const int dim = 3;
-   const int order = 1;
-
-   const double tol = 1e-6;
-
-   for (int type = (int)Element::TETRAHEDRON;
-        type <= (int)Element::HEXAHEDRON; type++)
-   {
-      Mesh mesh = Mesh::MakeCartesian3D(
-                     n, n, n, (Element::Type)type, 2.0, 3.0, 5.0);
-
-      VectorFunctionCoefficient funcCoef(dim, Func_3D_lin);
-
-      SECTION("3D GetVectorValue tests for element type " +
-              std::to_string(type))
-      {
-         RT_FECollection  rt_fec(order+1, dim);
-
-         FiniteElementSpace  rt_fespace(&mesh,  &rt_fec);
-
-         GridFunction  rt_x( &rt_fespace);
-
-         VectorGridFunctionCoefficient  rt_xCoef( &rt_x);
-
-         Array<int> bdr_marker(6);
-
-         Vector normal(dim);
-         Vector  f_val(dim);
-         Vector rt_val(dim);
-
-         for (int b = 1; b<=6; b++)
-         {
-            bdr_marker = 0;
-            bdr_marker[b-1] = 1;
-
-            rt_x = 0.0;
-            rt_x.ProjectBdrCoefficientNormal(funcCoef, bdr_marker);
-
-            for (int be = 0; be < mesh.GetNBE(); be++)
-            {
-               Element *e = mesh.GetBdrElement(be);
-               if (e->GetAttribute() != b) { continue; }
-
-               ElementTransformation *T = mesh.GetBdrElementTransformation(be);
-               const FiniteElement   *fe = rt_fespace.GetBE(be);
-               const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(),
-                                                        2*order + 2);
-
-               double rt_err = 0.0;
-
-               for (int j=0; j<ir.GetNPoints(); j++)
-               {
-                  const IntegrationPoint &ip = ir.IntPoint(j);
-                  T->SetIntPoint(&ip);
-
-                  CalcOrtho(T->Jacobian(), normal);
-
-                  funcCoef.Eval(f_val, *T, ip);
-                  rt_xCoef.Eval(rt_val, *T, ip);
-
-                  rt_val -= f_val;
-
-                  double rt_dist = rt_val * normal;
-
-                  rt_err += rt_dist;
-
-                  if (verbose_tests && rt_dist > tol)
-                  {
-                     mfem::out << be << ":" << j << " rt ("
-                               << f_val[0] << "," << f_val[1] << "," << f_val[2]
-                               << ") vs. ("
-                               << rt_val[0] << "," << rt_val[1] << ","
-                               << rt_val[2] << ") " << rt_dist << std::endl;
-                  }
-               }
-               rt_err  /= ir.GetNPoints();
-
-               REQUIRE( rt_err == MFEM_Approx(0.0));
-            }
-         }
-      }
-   }
-}
-
-TEST_CASE("3D ProjectBdrCoefficientNormal Scalar",
-          "[GridFunction]"
-          "[VectorGridFunctionCoefficient]")
-{
-   const int n = 1;
-   const int dim = 3;
-   const int order = 1;
-
-   const double tol = 1e-6;
-
-   const char bdrs_axis[] = {2, 1, 0, 1, 0, 2};
-   const char bdrs_sign[] = {-1, -1, +1, +1, -1, +1};
-
-   for (int type = (int)Element::TETRAHEDRON;
-        type <= (int)Element::HEXAHEDRON; type++)
-   {
-      Mesh mesh = Mesh::MakeCartesian3D(
-                     n, n, n, (Element::Type)type, 2.0, 3.0, 5.0);
-
-      VectorFunctionCoefficient funcCoef(dim, Func_3D_lin);
-
-      SECTION("3D GetVectorValue tests for element type " +
-              std::to_string(type))
-      {
-         RT_FECollection  rt_fec(order+1, dim);
-
-         FiniteElementSpace  rt_fespace(&mesh,  &rt_fec);
-
-         GridFunction  rt_x( &rt_fespace);
-
-         VectorGridFunctionCoefficient  rt_xCoef( &rt_x);
-
-         Array<int> bdr_marker(6);
-
-         Vector normal(dim);
-         Vector  f_val(dim);
-         Vector rt_val(dim);
-
-         for (int b = 1; b<=6; b++)
-         {
-            bdr_marker = 0;
-            bdr_marker[b-1] = 1;
-
-            rt_x = 0.0;
-
-            normal = 0.;
-            normal(bdrs_axis[b-1]) = (bdrs_sign[b-1] > 0)?(+1.):(-1.);
-            VectorConstantCoefficient normCoef(normal);
-            InnerProductCoefficient prodCoef(funcCoef, normCoef);
-            rt_x.ProjectBdrCoefficientNormal(prodCoef, bdr_marker);
-
-            for (int be = 0; be < mesh.GetNBE(); be++)
-            {
-               Element *e = mesh.GetBdrElement(be);
-               if (e->GetAttribute() != b) { continue; }
-
-               ElementTransformation *T = mesh.GetBdrElementTransformation(be);
-               const FiniteElement   *fe = rt_fespace.GetBE(be);
-               const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(),
-                                                        2*order + 2);
-
-               double rt_err = 0.0;
-
-               for (int j=0; j<ir.GetNPoints(); j++)
-               {
-                  const IntegrationPoint &ip = ir.IntPoint(j);
-                  T->SetIntPoint(&ip);
-
-                  CalcOrtho(T->Jacobian(), normal);
-
-                  funcCoef.Eval(f_val, *T, ip);
-                  rt_xCoef.Eval(rt_val, *T, ip);
-
-                  rt_val -= f_val;
-
-                  double rt_dist = rt_val * normal;
-
-                  rt_err += rt_dist;
-
-                  if (verbose_tests && rt_dist > tol)
-                  {
-                     mfem::out << be << ":" << j << " rt ("
-                               << f_val[0] << "," << f_val[1] << "," << f_val[2]
-                               << ") vs. ("
-                               << rt_val[0] << "," << rt_val[1] << ","
-                               << rt_val[2] << ") " << rt_dist << std::endl;
-                  }
-               }
-               rt_err  /= ir.GetNPoints();
-
-               REQUIRE( rt_err == MFEM_Approx(0.0));
-            }
-         }
-      }
-   }
-}
-
 TEST_CASE("3D ProjectBdrCoefficientTangent",
          "[GridFunction]"
          "[VectorGridFunctionCoefficient]")
 {
-   const int n = 1;
-   const int dim = 3;
-   const int order = 1;
+   int n = 1;
+   int dim = 3;
+   int order = 1;

-   const double tol = 1e-6;
+   double tol = 1e-6;

   for (int type = (int)Element::TETRAHEDRON;
        type <= (int)Element::HEXAHEDRON; type++)
@@ -271,8 +271,6 @@ TEST_CASE("Variable Order FiniteElementSpace",

      const auto space_type = GENERATE(SpaceType::RT, SpaceType::ND);
      const int dim = GENERATE(2, 3);
-      CAPTURE(space_type);
-      CAPTURE(dim);

      Mesh mesh = MakeCartesianMesh(dim == 2 ? 4 : 2, dim);
      mesh.EnsureNCMesh();
@@ -700,14 +698,7 @@ static void TestSolveVec(FiniteElementSpace &fespace)

   GridFunction x(&fespace);
   x = 0.0;
-   if (x.FESpace()->GetTypicalBE()->GetRangeDim() == 0)
-   {
-      x.ProjectBdrCoefficientNormal(exsol, ess_attr);
-   }
-   else
-   {
-      x.ProjectBdrCoefficientTangent(exsol, ess_attr);
-   }
+   x.ProjectBdrCoefficient(exsol, ess_attr);

   // Assemble the linear form
   LinearForm lf(&fespace);
@@ -1091,14 +1082,7 @@ static void TestSolveParVec(ParFiniteElementSpace &fespace)

   ParGridFunction x(&fespace);
   x = 0.0;
-   if (x.FESpace()->GetTypicalBE()->GetRangeDim() == 0)
-   {
-      x.ProjectBdrCoefficientNormal(exsol, ess_attr);
-   }
-   else
-   {
-      x.ProjectBdrCoefficientTangent(exsol, ess_attr);
-   }
+   x.ProjectBdrCoefficient(exsol, ess_attr);

   // Assemble the linear form
   ParLinearForm lf(&fespace);
@@ -200,39 +200,3 @@ TEST_CASE("ArraysByName Sort/Unique Methods", "[ArraysByName]")
      }
   }
 }
-
-TEST_CASE("ArraysByName Print/Load Methods", "[ArraysByName]")
-{
-   ArraysByName<int> abn;
-
-   FillArraysByName(abn);
-
-   // Print object to string using default format
-   std::ostringstream oss1;
-   abn.Print(oss1);
-
-   // Load new object from printed output
-   ArraysByName<int> abn_load1;
-   std::istringstream iss1(oss1.str());
-   abn_load1.Load(iss1);
-   REQUIRE(abn == abn_load1);
-
-   // Print object to string using one line per array
-   std::ostringstream oss2;
-   oss2 << abn.Size() << '\n';
-   for (auto a : abn)
-   {
-      oss2 << '"' << a.first << "\" " << a.second.Size();
-      for (auto d : a.second)
-      {
-         oss2 << ' ' << d;
-      }
-      oss2 << '\n';
-   }
-
-   // Load new object from printed output
-   ArraysByName<int> abn_load2;
-   std::istringstream iss2(oss2.str());
-   abn_load2.Load(iss2);
-   REQUIRE(abn == abn_load2);
-}
--- a/Show More
+++ b/Show More
Author	SHA1	Message	Date
Brandon Talamini	4b13e46bc5	Fix style error	2026-03-29 17:12:36 -07:00
Brandon Talamini	de114112b5	Remove unused variable from test	2026-03-27 21:08:43 -07:00
Brandon Talamini	ed1c21f656	Apply style to tensor.hpp	2026-03-27 21:01:16 -07:00
Brandon Talamini	811254989f	Get rid of designated initializers	2026-03-27 20:56:18 -07:00
Brandon Talamini	1f639a826f	Fix ODR problems caused by dropped inline keywords	2026-03-27 20:32:59 -07:00
Brandon Talamini	41e6f97bc7	Use relative paths for all includes	2026-03-27 20:02:46 -07:00
Brandon Talamini	7045b49973	Fix path to config.hpp	2026-03-27 19:57:40 -07:00
Brandon Talamini	91a121e427	Rename a badly named variable	2026-03-27 19:43:08 -07:00
Brandon Talamini	1f6d6ba82a	Style	2026-03-27 17:31:20 -07:00
Brandon Talamini	379ec4a393	Move test to dfem directory	2026-03-27 17:09:15 -07:00
Brandon Talamini	34895e545e	Cleanup through comments and renaming	2026-03-27 17:02:40 -07:00
Brandon Talamini	e0659eecd4	Guard tests for enzyme build only	2026-03-27 17:02:40 -07:00
Brandon Talamini	20b774792a	reorder things in test for clarity	2026-03-27 17:02:40 -07:00
Brandon Talamini	227bbe0283	Remove debug prints and add some comments	2026-03-27 17:02:40 -07:00
Brandon Talamini	88e8055d54	Remove some debug output and simplify test	2026-03-27 17:02:40 -07:00
Brandon Talamini	8e968658b5	Clean up tolerances and settings	2026-03-27 17:02:40 -07:00
Brandon Talamini	978e2bbd99	Fix Enzyme compilation error in Release mode Credit to William Moses.	2026-03-27 17:02:40 -07:00
Brandon Talamini	6465d872b5	Comment and move things around for clarity	2026-03-27 17:02:40 -07:00
Brandon Talamini	0439e51f6b	Make bounds for nth root functions better so they will work for any input	2026-03-27 17:02:40 -07:00
Brandon Talamini	e2917c9e15	Clamp initial guess with brackets, no warning if guess is moved	2026-03-27 17:02:40 -07:00
Brandon Talamini	f0ea358c36	Change lack of root bracketing handling from an assertion to a warning so the user is alerted in Release builds	2026-03-27 17:02:40 -07:00
Brandon Talamini	c26460c0ce	Change divergence error to an MFEM_ABORT macro	2026-03-27 17:02:40 -07:00
Brandon Talamini	9975f4c9db	Inline tensor operations so that Enzyme can analyze enzyme_const operations correctly	2026-03-27 17:02:40 -07:00
Brandon Talamini	b29d5856c2	Solve problem where setting bounds that depended on state caused incorrect derivatives	2026-03-27 17:02:40 -07:00
Brandon Talamini	ce54fdf7b4	Put in a robustness test for a case where standard Newton diverges	2026-03-27 17:02:40 -07:00
Brandon Talamini	7ca58c106a	Clean up tests: comment and use non-halting test conditions	2026-03-27 17:02:39 -07:00
Brandon Talamini	63056ab08b	Put in a real test condition for the VJP	2026-03-27 17:02:39 -07:00
Brandon Talamini	cd9456db2a	Get reverse mode working Derivatives become nans if solver upper bound depends on state. Need to debug this with Bill. But derivatives are now being computed, and appear correct.	2026-03-27 17:02:39 -07:00
Brandon Talamini	7a8fb3ca79	Fixing Enzyme problems with Bill Moses	2026-03-27 17:02:39 -07:00
Brandon Talamini	dc251d9518	Move solver from linalg to dfem	2026-03-27 17:02:39 -07:00
Brandon Talamini	2c118aaf85	Move implmentation inside separate namespace	2026-03-27 17:02:39 -07:00
Brandon Talamini	e0e32976ec	Fix custom revrse mode	2026-03-27 17:02:39 -07:00
Brandon Talamini	ac71a377e1	Make wrapper for residual use values and references so that it is more symmetric with residual signature	2026-03-27 17:02:39 -07:00
Brandon Talamini	ac9b4d3d38	Make a basic test of reverse mode and get it to pass	2026-03-27 17:02:39 -07:00
Brandon Talamini	6be47674d6	Implement custom reverse mode derivative of solver	2026-03-27 17:02:39 -07:00
Brandon Talamini	75b6f2ac19	Check enzyme derivative with finite difference	2026-03-27 17:02:39 -07:00
Brandon Talamini	76c805cbad	Put in a test that computes derivative of q-function with solver inside	2026-03-27 17:02:39 -07:00
Brandon Talamini	592f87591d	Fix bugs to make tests pass	2026-03-27 17:02:39 -07:00
Brandon Talamini	5decf944cb	Fix bugs in error checking Forgot to reverse inequalities when changing from ERROR_IF to ASSERT checks.	2026-03-27 17:02:39 -07:00
Brandon Talamini	6eccede39e	Add a unit test of new solver	2026-03-27 17:02:39 -07:00
Brandon Talamini	f8fd930618	Fix compilation errors	2026-03-27 17:02:39 -07:00
Brandon Talamini	c450622e6e	Initial draft of solver	2026-03-27 17:02:39 -07:00