Cleanup pass

Merge branch 'master' into checkpointing_examples_gpu
Cleanup
2026-03-01 15:30:39 -08:00 · 2026-03-01 14:38:23 -08:00 · 2026-03-01 14:34:24 -08:00 · 2026-02-10 12:43:39 -08:00 · 2026-02-09 00:23:04 -08:00 · 2026-02-08 22:07:05 -08:00
150 changed files with 12006 additions and 6877 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)
 ===========================

@@ -652,8 +652,6 @@ foreach(TPL IN LISTS MFEM_TPLS)
  endif()
 endforeach(TPL)

-# reverse to remove the first instance of entries in TPL_LIBRARIES
-# so later duplicates are kept (for dependency ordering)
 list(REVERSE TPL_LIBRARIES)
 list(REMOVE_DUPLICATES TPL_LIBRARIES)
 list(REVERSE TPL_LIBRARIES)
@@ -1017,7 +1015,5 @@ install(DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/data
 # Create 'config.mk' from 'config.mk.in' for the build and install locations and
 # define install rules for 'config.mk' and 'test.mk'
 #-------------------------------------------------------------------------------
-if (MFEM_USE_CUDA OR MFEM_USE_HIP)
-  option(MFEM_EXPORT_GPU_CONFIG "Export config.mk for GPU-enabled downstream packages" ON)
-endif()
+
 mfem_export_mk_files()
@@ -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")
@@ -701,6 +701,7 @@ endfunction(mfem_find_library)
 # Extract compile and link options needed by the given target.
 #
 function(mfem_get_target_options Target CompileOptsVar LinkOptsVar)
+
  if (NOT TARGET ${Target})
    return()
  endif()
@@ -798,12 +799,7 @@ function(mfem_get_target_options Target CompileOptsVar LinkOptsVar)
          # message(STATUS "Lib = ${Lib}")
          # Filter-out generator expressions
          if (NOT ("${Lib}" MATCHES "^\\$"))
-            if(NOT ("${Lib}" STREQUAL "dl"))
-              list(APPEND LinkOpts "${Lib}")
-            else()
-              # for some reason libdl doesn't include the "-l"
-              list(APPEND LinkOpts "-ldl")
-            endif()
+            list(APPEND LinkOpts "${Lib}")
          endif()
        else()
          mfem_get_target_options(${Lib} COpts LOpts)
@@ -892,18 +888,9 @@ function(mfem_export_mk_files)
      set(${var} NO)
    endif()
  endforeach()
-  if (MFEM_USE_CUDA AND MFEM_EXPORT_GPU_CONFIG)
-    set(MFEM_CXX ${CMAKE_CUDA_COMPILER})
-    if(MFEM_CUDA_COMPILER_IS_NVCC)
-      set(MFEM_HOST_CXX ${CMAKE_CUDA_HOST_COMPILER})
-    else()
-      set(MFEM_HOST_CXX ${CMAKE_CXX_COMPILER})
-    endif()
-  else()
-    # mfem doesn't use enable_language(HIP)
-    set(MFEM_CXX ${CMAKE_CXX_COMPILER})
-    set(MFEM_HOST_CXX ${CMAKE_CXX_COMPILER})
-  endif()
+  # TODO: Add support for MFEM_USE_CUDA=YES
+  set(MFEM_CXX ${CMAKE_CXX_COMPILER})
+  set(MFEM_HOST_CXX ${MFEM_CXX})
  set(MFEM_CPPFLAGS "")
  get_target_property(cxx_std mfem CXX_STANDARD)
  # For now, we ignore the setting of the CXX_EXTENSIONS property. If this
@@ -913,50 +900,6 @@ function(mfem_export_mk_files)
  string(STRIP
         "${cxx_std_flag} ${CMAKE_CXX_FLAGS_${BUILD_TYPE}} ${CMAKE_CXX_FLAGS}"
         MFEM_CXXFLAGS)
-  if(MFEM_EXPORT_GPU_CONFIG)
-    if (MFEM_USE_CUDA)
-      set(MFEM_CXXFLAGS "${MFEM_CXXFLAGS} ${CMAKE_CUDA_FLAGS}")
-      if (MFEM_CUDA_COMPILER_IS_NVCC)
-        set(MFEM_CXXFLAGS "-x=cu ${MFEM_CXXFLAGS} -ccbin ${CMAKE_CXX_COMPILER} --forward-unknown-to-host-compiler")
-        # The following intentionally hides CUDA deprecation warnings
-        foreach(ENTRY IN LISTS CUDAToolkit_INCLUDE_DIRS)
-          set(MFEM_CXXFLAGS "${MFEM_CXXFLAGS} -isystem ${ENTRY}")
-        endforeach()
-        if (CMAKE_VERSION VERSION_GREATER_EQUAL 3.18.0)
-          # architecture flags not part of CMAKE_CUDA_FLAGS
-          if ("all" STREQUAL "${CMAKE_CUDA_ARCHITECTURES}"
-              OR "native" STREQUAL "${CMAKE_CUDA_ARCHITECTURES}"
-              OR "all-major" STREQUAL "${CMAKE_CUDA_ARCHITECTURES}")
-            set(MFEM_CXXFLAGS "${MFEM_CXXFLAGS} -arch=${CMAKE_CUDA_ARCHITECTURES}")
-          else()
-            foreach (ENTRY IN LISTS CMAKE_CUDA_ARCHITECTURES)
-              set(MFEM_CXXFLAGS
-                "${MFEM_CXXFLAGS} -gencode arch=compute_${ENTRY},code=sm_${ENTRY}")
-            endforeach()
-          endif()
-        endif()
-      else()
-        set(MFEM_CXXFLAGS "${MFEM_CXXFLAGS} -xcuda --cuda-path=${CUDAToolkit_LIBRARY_ROOT}")
-        if (CMAKE_VERSION VERSION_GREATER_EQUAL 3.18.0)
-          # architecture flags not part of CMAKE_CUDA_FLAGS
-          if ("all" STREQUAL "${CMAKE_CUDA_ARCHITECTURES}"
-              OR "native" STREQUAL "${CMAKE_CUDA_ARCHITECTURES}"
-              OR "all-major" STREQUAL "${CMAKE_CUDA_ARCHITECTURES}")
-            # TODO: not supported
-          else()
-            foreach(ENTRY IN LISTS CMAKE_CUDA_ARCHITECTURES)
-              set(MFEM_CXXFLAGS "-cuda-gpu-arch=sm_${ENTRY} ${MFEM_CXXFLAGS}")
-            endforeach()
-          endif()
-        endif()
-      endif()
-    elseif (MFEM_USE_HIP)
-      set(MFEM_CXXFLAGS "${MFEM_CXXFLAGS} -xhip")
-      foreach(ENTRY IN LISTS CMAKE_HIP_ARCHITECTURES)
-        set(MFEM_CXXFLAGS "--offload-arch=${ENTRY} ${MFEM_CXXFLAGS}")
-      endforeach()
-    endif()
-  endif()
  set(MFEM_TPLFLAGS "")
  foreach(dir ${TPL_INCLUDE_DIRS})
    set(MFEM_TPLFLAGS "${MFEM_TPLFLAGS} -I${dir}")
@@ -987,9 +930,6 @@ function(mfem_export_mk_files)
    set(MFEM_SHARED NO)
    set(MFEM_STATIC YES)
  endif()
-  if (MFEM_USE_CUDA)
-    set(MFEM_EXT_LIBS "${MFEM_EXT_LIBS} -lcudart")
-  endif()
  set(MFEM_BUILD_TAG "${CMAKE_SYSTEM}")
  set(MFEM_PREFIX "${CMAKE_INSTALL_PREFIX}")
  # For the next 4 variables, these are the values for the build-tree version of
@@ -998,15 +938,8 @@ function(mfem_export_mk_files)
  set(MFEM_LIB_DIR "${PROJECT_BINARY_DIR}")
  set(MFEM_TEST_MK "${PROJECT_SOURCE_DIR}/config/test.mk")
  set(MFEM_CONFIG_EXTRA "MFEM_BUILD_DIR ?= ${PROJECT_BINARY_DIR}")
-  if (MFEM_USE_CUDA AND MFEM_EXPORT_GPU_CONFIG)
-    if (MFEM_CUDA_COMPILER_IS_NVCC)
-      set(MFEM_XLINKER "-Xlinker=")
-    else()
-      set(MFEM_XLINKER "${CMAKE_CUDA_LINKER_WRAPPER_FLAG}")
-    endif()
-  else()
-    set(MFEM_XLINKER "${CMAKE_CXX_LINKER_WRAPPER_FLAG}")
-  endif()
+  # TODO: CUDA/HIP support:
+  set(MFEM_XLINKER "${CMAKE_CXX_LINKER_WRAPPER_FLAG}")
  set(MFEM_MPIEXEC ${MPIEXEC})
  if (NOT MFEM_MPIEXEC)
    set(MFEM_MPIEXEC "mpirun")
@@ -1054,21 +987,16 @@ function(mfem_export_mk_files)
    # handle interfaces (e.g., SCOREC::apf)
    if ("${lib}" MATCHES "SCOREC::.*" OR "${lib}" MATCHES "Ginkgo::.*" OR "${lib}" MATCHES "ParMoonolith::.*")
    elseif (TARGET "${lib}")
-      mfem_get_target_options(${lib} CompileOpts2 LinkOpts2)
-      # remove generator expressions
-      string(GENEX_STRIP "${CompileOpts2}" CompileOpts)
-      string(GENEX_STRIP "${LinkOpts2}" LinkOpts)
+      mfem_get_target_options(${lib} CompileOpts LinkOpts)
      # Removing duplicates may lead to issues:
      # list(REMOVE_DUPLICATES CompileOpts)
      # list(REMOVE_DUPLICATES LinkOpts)
-      # message(WARNING "${lib}[LinkOpts]: ${LinkOpts}")      
-      # message(WARNING "${lib}[CompileOpts]: ${CompileOpts}")
-      foreach(LOpt IN LISTS LinkOpts)
-        set(MFEM_EXT_LIBS "${MFEM_EXT_LIBS} ${LOpt}")
-      endforeach()
-      foreach(COpt IN LISTS CompileOpts)
-        set(MFEM_TPLFLAGS "${MFEM_TPLFLAGS} ${COpt}")
-      endforeach()
+      string(REPLACE ";" " " COpts "${CompileOpts}")
+      string(REPLACE ";" " " LOpts "${LinkOpts}")
+      # message(STATUS "${lib}[COpts]: '${COpts}'")
+      # message(STATUS "${lib}[LOpts]: '${LOpts}'")
+      set(MFEM_TPLFLAGS "${MFEM_TPLFLAGS} ${COpts}")
+      set(MFEM_EXT_LIBS "${MFEM_EXT_LIBS} ${LOpts}")
      # message(FATAL_ERROR "***** interface lib found ... exiting *****")
      # handle static and shared libs
    elseif ("${suffix}" STREQUAL "${CMAKE_SHARED_LIBRARY_SUFFIX}")
@@ -1076,7 +1004,7 @@ function(mfem_export_mk_files)
      get_filename_component(fullLibName ${lib} NAME_WE)
      string(REGEX REPLACE "^lib" "" libname ${fullLibName})
      set(MFEM_EXT_LIBS
-        "${MFEM_EXT_LIBS} ${shared_link_flag}${dir} -L${dir} -l${libname}")
+          "${MFEM_EXT_LIBS} ${shared_link_flag}${dir} -L${dir} -l${libname}")
    else()
      set(MFEM_EXT_LIBS "${MFEM_EXT_LIBS} ${lib}")
    endif()
@@ -1085,7 +1013,7 @@ function(mfem_export_mk_files)
  # Create the build-tree version of 'config.mk'
  configure_file(
    "${PROJECT_SOURCE_DIR}/config/config.mk.in"
-    "${PROJECT_BINARY_DIR}/config/config.mk" @ONLY)
+    "${PROJECT_BINARY_DIR}/config/config.mk")
  # Copy 'test.mk' from the source-tree to the build-tree
  configure_file(
    "${PROJECT_SOURCE_DIR}/config/test.mk"
@@ -1103,7 +1031,7 @@ function(mfem_export_mk_files)
  # Create the install-tree version of 'config.mk'
  configure_file(
    "${PROJECT_SOURCE_DIR}/config/config.mk.in"
-    "${PROJECT_BINARY_DIR}/config/config-install.mk" @ONLY)
+    "${PROJECT_BINARY_DIR}/config/config-install.mk")

  # Install rules for 'config.mk' and 'test.mk'
  install(FILES ${PROJECT_SOURCE_DIR}/config/test.mk
@@ -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
@@ -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
@@ -82,25 +82,6 @@ public:
   /// underlying #fes
   int VectorDim() const;

-   /// Copy assignment. Only the data of the base class Vector is copied.
-   /** It is assumed that this object and @a rhs use FiniteElementSpace%s that
-       have the same size.
-
-       @note Defining this method overwrites the implicitly defined copy
-       assignment operator. */
-   ComplexGridFunction &operator=(const ComplexGridFunction &rhs)
-   { return operator=((const Vector &)rhs); }
-
-   /// Copy the data from @a v.
-   /** The size of @a v must be equal to double of the size of the associated
-       FiniteElementSpace #fes. */
-   ComplexGridFunction &operator=(const Vector &v)
-   {
-      MFEM_ASSERT(fes && v.Size() == 2*fes->GetVSize(), "");
-      Vector::operator=(v);
-      return *this;
-   }
-
   /// Assign constant values to the ComplexGridFunction data.
   ComplexGridFunction &operator=(const std::complex<real_t> & value)
   { *gfr = value.real(); *gfi = value.imag(); return *this; }
@@ -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";
@@ -84,29 +84,33 @@ public:
      const std::vector<derivative_action_t> &derivative_actions,
      const FieldDescriptor &direction,
      const int &daction_l_size,
-      const std::vector<derivative_action_t> &derivative_actions_transpose,
+      const int &derivative_action_tr_l_size,
+      const std::vector<derivative_action_t> &derivative_tr_actions,
      const FieldDescriptor &transpose_direction,
      const int &daction_transpose_l_size,
      const std::vector<Vector *> &solutions_l,
      const std::vector<Vector *> &parameters_l,
      const restriction_callback_t &restriction_callback,
      const std::function<void(Vector &, Vector &)> &prolongation_transpose,
+      const std::function<void(Vector &, Vector &)> &tr_prolongation_transpose,
      const std::vector<assemble_derivative_sparsematrix_callback_t>
      &assemble_derivative_sparsematrix_callbacks,
-      const std::vector<assemble_derivative_hypreparmatrix_callback_t>
-      &assemble_derivative_hypreparmatrix_callbacks) :
+      const assemble_derivative_hypreparmatrix_callback_t
+      &assemble_derivative_hypreparmatrix_callback) :
      Operator(height, width),
      derivative_actions(derivative_actions),
      direction(direction),
      daction_l(daction_l_size),
      daction_l_size(daction_l_size),
-      derivative_actions_transpose(derivative_actions_transpose),
+      derivative_action_tr_l_size(derivative_action_tr_l_size),
+      derivative_tr_actions(derivative_tr_actions),
      transpose_direction(transpose_direction),
      prolongation_transpose(prolongation_transpose),
+      tr_prolongation_transpose(tr_prolongation_transpose),
      assemble_derivative_sparsematrix_callbacks(
         assemble_derivative_sparsematrix_callbacks),
-      assemble_derivative_hypreparmatrix_callbacks(
-         assemble_derivative_hypreparmatrix_callbacks)
+      assemble_derivative_hypreparmatrix_callback(
+         assemble_derivative_hypreparmatrix_callback)
   {
      std::vector<Vector> s_l(solutions_l.size());
      for (size_t i = 0; i < s_l.size(); i++)
@@ -156,18 +160,18 @@ public:
   /// direction_t on T-dofs.
   void MultTranspose(const Vector &direction_t, Vector &result_t) const override
   {
-      MFEM_ASSERT(!derivative_actions_transpose.empty(),
+      MFEM_ASSERT(!derivative_tr_actions.empty(),
                  "derivative can't be used to be multiplied in transpose mode");

-      daction_l.SetSize(width);
+      daction_l.SetSize(derivative_action_tr_l_size);
      daction_l = 0.0;

      prolongation(transpose_direction, direction_t, direction_l);
-      for (const auto &f : derivative_actions_transpose)
+      for (const auto &f : derivative_tr_actions)
      {
         f(fields_e, direction_l, daction_l);
      }
-      prolongation_transpose(daction_l, result_t);
+      tr_prolongation_transpose(daction_l, result_t);
   };

   /// @brief Assemble the derivative operator into a SparseMatrix.
@@ -183,6 +187,10 @@ public:
      {
         f(fields_e, A);
      }
+
+      // SparseMatrix A is finalized after all callbacks have contributed to
+      // it.
+      A->Finalize();
   }

   /// @brief Assemble the derivative operator into a HypreParMatrix.
@@ -191,13 +199,7 @@ public:
   /// be an uninitialized object.
   void Assemble(HypreParMatrix *&A)
   {
-      MFEM_ASSERT(!assemble_derivative_hypreparmatrix_callbacks.empty(),
-                  "derivative can't be assembled into a HypreParMatrix");
-
-      for (const auto &f : assemble_derivative_hypreparmatrix_callbacks)
-      {
-         f(fields_e, A);
-      }
+      assemble_derivative_hypreparmatrix_callback(fields_e, A);
   }

 private:
@@ -212,10 +214,12 @@ private:

   const int daction_l_size;

+   const int derivative_action_tr_l_size;
+
   /// Transpose Derivative action callbacks. Depending on the requested
   /// derivatives in DifferentiableOperator the callbacks represent certain
   /// combinations of actions of derivatives of the forward operator.
-   std::vector<derivative_action_t> derivative_actions_transpose;
+   std::vector<derivative_action_t> derivative_tr_actions;

   FieldDescriptor transpose_direction;

@@ -225,13 +229,15 @@ private:

   std::function<void(Vector &, Vector &)> prolongation_transpose;

+   std::function<void(Vector &, Vector &)> tr_prolongation_transpose;
+
   /// Callbacks that assemble derivatives into a SparseMatrix.
   std::vector<assemble_derivative_sparsematrix_callback_t>
   assemble_derivative_sparsematrix_callbacks;

   /// Callbacks that assemble derivatives into a HypreParMatrix.
-   std::vector<assemble_derivative_hypreparmatrix_callback_t>
-   assemble_derivative_hypreparmatrix_callbacks;
+   assemble_derivative_hypreparmatrix_callback_t
+   assemble_derivative_hypreparmatrix_callback;
 };

 /// Class representing a differentiable operator which acts on solution and
@@ -457,7 +463,10 @@ public:
         dir_l = s_l[derivative_idx];
      }

-      derivative_setup_callbacks[derivative_id][0](fields_e, dir_l);
+      for (size_t i = 0; i < derivative_setup_callbacks[derivative_id].size(); i++)
+      {
+         derivative_setup_callbacks[derivative_id][i](fields_e, dir_l);
+      }

      return std::make_shared<DerivativeOperator>(
                height,
@@ -465,15 +474,17 @@ public:
                derivative_action_callbacks[derivative_id],
                fields[derivative_idx],
                residual_l.Size(),
-                daction_transpose_callbacks[derivative_id],
+                derivative_action_tr_l_size[derivative_id],
+                derivative_action_tr_callbacks[derivative_id],
                fields[test_space_field_idx],
                GetVSize(fields[test_space_field_idx]),
                sol_l,
                par_l,
                restriction_callback,
                prolongation_transpose,
+                derivative_tr_prolongation_transpose[derivative_id],
                assemble_derivative_sparsematrix_callbacks[derivative_id],
-                assemble_derivative_hypreparmatrix_callbacks[derivative_id]);
+                assemble_derivative_hypreparmatrix_callback[derivative_id]);
   }

 private:
@@ -486,13 +497,15 @@ private:
   std::map<size_t,
       std::vector<derivative_action_t>> derivative_action_callbacks;
   std::map<size_t,
-       std::vector<derivative_action_t>> daction_transpose_callbacks;
+       std::vector<derivative_action_t>> derivative_action_tr_callbacks;
+   std::map<size_t,
+       std::function<void(Vector &, Vector &)>> derivative_tr_prolongation_transpose;
+   std::map<size_t, int> derivative_action_tr_l_size;
   std::map<size_t,
       std::vector<assemble_derivative_sparsematrix_callback_t>>
       assemble_derivative_sparsematrix_callbacks;
-   std::map<size_t,
-       std::vector<assemble_derivative_hypreparmatrix_callback_t>>
-       assemble_derivative_hypreparmatrix_callbacks;
+   std::map<size_t, assemble_derivative_hypreparmatrix_callback_t>
+   assemble_derivative_hypreparmatrix_callback;

   std::vector<FieldDescriptor> solutions;
   std::vector<FieldDescriptor> parameters;
@@ -510,7 +523,7 @@ private:
   std::function<void(Vector &, Vector &)> output_restriction_transpose;
   restriction_callback_t restriction_callback;

-   std::map<size_t, Vector> derivative_qp_caches;
+   std::map<size_t, std::vector<Vector>> derivative_qp_caches;

   std::map<size_t, size_t> assembled_vector_sizes;

@@ -769,9 +782,10 @@ void DifferentiableOperator::AddIntegrator(
   auto input_size_on_qp =
      get_input_size_on_qp(inputs, std::make_index_sequence<num_inputs> {});

+   // printf("calculate shmem action info\n");
   auto action_shmem_info =
      get_shmem_info<entity_t, num_fields, num_inputs, num_outputs>
-      (input_dtq_maps, output_dtq_maps, fields, num_entities, inputs, num_qp,
+      (input_dtq_maps, output_dtq_maps, fields, num_entities, num_qp,
       input_size_on_qp, residual_size_on_qp, element_dof_ordering);

   Vector shmem_cache(action_shmem_info.total_size);
@@ -894,9 +908,10 @@ void DifferentiableOperator::AddIntegrator(
         const int da_size_on_qp =
            GetSizeOnQP<entity_t>(output_fop, fields[test_space_field_idx]);

+         // printf("calculate shmem derivative action info\n");
         auto shmem_info =
            get_shmem_info<entity_t, num_fields, num_inputs, num_outputs>(
-               input_dtq_maps, output_dtq_maps, fields, num_entities, inputs,
+               input_dtq_maps, output_dtq_maps, fields, num_entities,
               num_qp, input_size_on_qp, residual_size_on_qp,
               element_dof_ordering, d_field_idx);

@@ -968,9 +983,12 @@ void DifferentiableOperator::AddIntegrator(
         // Quadrature point local derivative cache for each element, with data
         // layout:
         // [test_vdim, test_op_dim, trial_vdim, trial_op_dim, qp, num_entities].
-         derivative_qp_caches[derivative_id] = Vector(test_vdim * test_op_dim *
-                                                      trial_vdim *
-                                                      total_trial_op_dim * num_qp * num_entities);
+         derivative_qp_caches[derivative_id].push_back(
+            Vector(test_vdim * test_op_dim * trial_vdim * total_trial_op_dim * num_qp *
+                   num_entities));
+
+         const int cache_index = this->derivative_qp_caches[derivative_id].size() - 1;
+
         // Create local references for MSVC lambda capture compatibility
         auto& fields_ref = this->fields;
         auto& derivative_qp_caches_ref = this->derivative_qp_caches[derivative_id];
@@ -1012,6 +1030,8 @@ void DifferentiableOperator::AddIntegrator(
               trial_vdim,
               inputs_trial_op_dim,

+               qpdc_idx = cache_index,
+
               // capture by ref:
               &qpdc_mem = derivative_qp_caches_ref
            ](std::vector<Vector> &f_e, const Vector &dir_l) mutable
@@ -1024,7 +1044,7 @@ void DifferentiableOperator::AddIntegrator(
                                               shmem_info.direction_size,
                                               num_entities);

-            auto qpdc = Reshape(qpdc_mem.ReadWrite(), test_vdim, test_op_dim,
+            auto qpdc = Reshape(qpdc_mem[qpdc_idx].ReadWrite(), test_vdim, test_op_dim,
                                trial_vdim, total_trial_op_dim, num_qp, num_entities);

            auto itod = Reshape(inputs_trial_op_dim.Read(), num_inputs);
@@ -1095,6 +1115,7 @@ void DifferentiableOperator::AddIntegrator(
               inputs_trial_op_dim,
               total_trial_op_dim,
               trial_vdim,
+               qpdc_idx = cache_index,
               // capture by ref:
               &qpdc_mem = derivative_qp_caches_ref,
               &or_transpose
@@ -1112,7 +1133,7 @@ void DifferentiableOperator::AddIntegrator(
                                               shmem_info.direction_size,
                                               num_entities);

-            auto qpdc = Reshape(qpdc_mem.Read(), test_vdim, test_op_dim,
+            auto qpdc = Reshape(qpdc_mem[qpdc_idx].Read(), test_vdim, test_op_dim,
                                trial_vdim, total_trial_op_dim, num_qp, num_entities);

            auto itod = Reshape(inputs_trial_op_dim.Read(), num_inputs);
@@ -1157,6 +1178,154 @@ void DifferentiableOperator::AddIntegrator(
            or_transpose(derivative_action_e, der_action_l);
         });

+         // This prevents Sum/Identity. These are invalid
+         // as input FieldOperators anyways.
+         constexpr auto dummy_fop = Value<0> {};
+
+         auto [input_rt,
+               input_e_sz] = get_restriction_transpose<entity_t>
+                             (fields[d_field_idx],
+                              element_dof_ordering, dummy_fop);
+
+         const auto input_restriction_transpose = input_rt;
+
+         derivative_tr_prolongation_transpose[derivative_id] =
+            get_prolongation_transpose(
+               fields[d_field_idx], dummy_fop, mesh.GetComm());
+
+         const auto d_tr_field_idx = test_space_field_idx;
+         const auto direction_tr = fields[d_tr_field_idx];
+
+         auto output_size_on_qp =
+            get_input_size_on_qp(outputs, std::make_index_sequence<num_outputs> {});
+
+         const int residual_tr_size_on_qp = trial_vdim * total_trial_op_dim;
+
+         auto shmem_tr_info =
+            get_shmem_info<entity_t, num_fields, num_outputs, num_inputs>(
+               output_dtq_maps, input_dtq_maps, fields, num_entities,
+               num_qp, output_size_on_qp, residual_tr_size_on_qp,
+               element_dof_ordering, test_space_field_idx);
+
+         // print_shared_memory_info(shmem_tr_info);
+
+         // TODO: this is a hack to extend the shared memory with a known
+         // offset for a temp variable
+         Vector shmem_tr_cache(shmem_tr_info.total_size + residual_tr_size_on_qp *
+                               num_qp);
+
+         Vector direction_tr_e(get_restriction<entity_t>(
+                                  fields[test_space_field_idx],
+                                  element_dof_ordering)->Height());
+
+         derivative_action_tr_l_size[derivative_id] =
+            get_restriction<entity_t>(fields[d_field_idx],
+                                      element_dof_ordering)->Width();
+
+         Vector derivative_action_tr_e(input_e_sz);
+         derivative_action_tr_e = 0.0;
+
+         derivative_action_tr_callbacks[derivative_id].push_back(
+            [
+               // capture by copy:
+               dimension,             // int
+               num_entities,          // int
+               num_trial_dof,         // int
+               num_qp,                // int
+               q1d,                   // int
+               test_vdim,             // int (= output_fop.vdim)
+               test_op_dim,           // int (derived from output_fop)
+               inputs,                // mfem::future::tuple
+               outputs,               // mfem::future::tuple
+               attributes,            // Array<int>
+               ir_weights,            // DeviceTensor
+               use_sum_factorization, // bool
+               input_dtq_maps,        // std::array<DofToQuadMap, num_fields>
+               output_dtq_maps,       // std::array<DofToQuadMap, num_fields>
+               // output_fop,            // class derived from FieldOperator
+               thread_blocks,         // ThreadBlocks
+               shmem_tr_cache,           // Vector (local)
+               shmem_tr_info,            // SharedMemoryInfo
+               // TODO: make this Array<int> a member of the DifferentiableOperator
+               //       and capture it by ref.
+               elem_attributes,       // Array<int>
+
+               input_is_dependent,
+               direction_tr,             // FieldDescriptor
+               direction_tr_e,           // Vector
+               derivative_action_tr_e,   // Vector
+               element_dof_ordering,  // ElementDofOrdering
+               inputs_trial_op_dim,
+               total_trial_op_dim,
+               trial_vdim,
+               input_restriction_transpose,
+               qpdc_idx = cache_index,
+               // capture by ref:
+               &qpdc_mem = derivative_qp_caches_ref
+            ](
+               std::vector<Vector> &f_e, const Vector &dir_tr_l,
+               Vector &derivative_action_tr_l) mutable
+         {
+            restriction<entity_t>(direction_tr, dir_tr_l, direction_tr_e,
+                                  element_dof_ordering);
+            auto ye = Reshape(derivative_action_tr_e.ReadWrite(), num_trial_dof,
+                              trial_vdim, num_entities);
+            auto wrapped_fields_e = wrap_fields(f_e, shmem_tr_info.field_sizes,
+                                                num_entities);
+            auto wrapped_direction_e = Reshape(direction_tr_e.ReadWrite(),
+                                               shmem_tr_info.direction_size,
+                                               num_entities);
+
+            auto qpdc = Reshape(qpdc_mem[qpdc_idx].Read(), test_vdim, test_op_dim,
+                                trial_vdim, total_trial_op_dim, num_qp, num_entities);
+
+            auto itod = Reshape(inputs_trial_op_dim.Read(), num_inputs);
+
+            const bool has_attr = attributes.Size() > 0;
+            const auto d_attr = attributes.Read();
+            const auto d_elem_attr = elem_attributes->Read();
+
+            derivative_action_tr_e = 0.0;
+            forall([=] MFEM_HOST_DEVICE (int e, real_t *shmem)
+            {
+               if (has_attr && !d_attr[d_elem_attr[e] - 1]) { return; }
+
+               auto [output_dtq_shmem, input_dtq_shmem, fields_shmem,
+                                       direction_shmem, input_shmem,
+                                       shadow_shmem_, residual_shmem,
+                                       scratch_shmem] =
+                        unpack_shmem(shmem, shmem_tr_info, output_dtq_maps, input_dtq_maps,
+                                     wrapped_fields_e, wrapped_direction_e, num_qp, e);
+               auto &shadow_shmem = shadow_shmem_;
+
+               std::array<bool, num_outputs> all_true{true};
+               map_direction_to_quadrature_data_conditional(
+                  shadow_shmem, direction_shmem, output_dtq_shmem, outputs,
+                  ir_weights, scratch_shmem, all_true, dimension,
+                  use_sum_factorization);
+
+               auto fhat = Reshape(&residual_shmem(0, 0), trial_vdim,
+                                   total_trial_op_dim, num_qp);
+
+               auto qpdce = Reshape(&qpdc(0, 0, 0, 0, 0, e), test_vdim, test_op_dim,
+                                    trial_vdim, total_trial_op_dim, num_qp);
+
+               constexpr bool transpose = true;
+               apply_qpdc(fhat, shadow_shmem, qpdce, itod, q1d, dimension,
+                          use_sum_factorization, transpose);
+
+               auto y = Reshape(&ye(0, 0, e), num_trial_dof, trial_vdim);
+               auto fi_shmem = Reshape(shmem + shmem_tr_info.total_size, trial_vdim,
+                                       total_trial_op_dim, num_qp);
+
+               map_quadrature_data_to_fields_conditional(
+                  y, fhat, inputs, itod, input_dtq_shmem, scratch_shmem, fi_shmem,
+                  input_is_dependent, dimension, use_sum_factorization);
+            }, num_entities, thread_blocks, shmem_tr_info.total_size,
+            shmem_tr_cache.ReadWrite());
+            input_restriction_transpose(derivative_action_tr_e, derivative_action_tr_l);
+         });
+
         assemble_derivative_sparsematrix_callbacks[derivative_id].push_back(
            [
               // capture by copy:
@@ -1190,7 +1359,7 @@ void DifferentiableOperator::AddIntegrator(
               inputs_trial_op_dim,
               Ae_mem,
               output_to_field,
-
+               qpdc_idx = cache_index,
               // capture by ref:
               &qpdc_mem = derivative_qp_caches_ref,
               &fields = fields_ref
@@ -1202,7 +1371,7 @@ void DifferentiableOperator::AddIntegrator(
                                               shmem_info.direction_size,
                                               num_entities);

-            auto qpdc = Reshape(qpdc_mem.Read(), test_vdim, test_op_dim,
+            auto qpdc = Reshape(qpdc_mem[qpdc_idx].Read(), test_vdim, test_op_dim,
                                trial_vdim, total_trial_op_dim, num_qp, num_entities);

            auto itod = Reshape(inputs_trial_op_dim.Read(), num_inputs);
@@ -1250,7 +1419,10 @@ void DifferentiableOperator::AddIntegrator(
            auto test_fes = *std::get_if<const ParFiniteElementSpace *>
                            (&fields[output_to_field[0]].data);

-            A = new SparseMatrix(test_fes->GetVSize(), trial_fes->GetVSize());
+            if (A == nullptr)
+            {
+               A = new SparseMatrix(test_fes->GetVSize(), trial_fes->GetVSize());
+            }

            auto tmp = Reshape(Ae_mem.HostReadWrite(), num_test_dof * test_vdim,
                               num_trial_dof * trial_vdim, num_entities);
@@ -1321,14 +1493,15 @@ void DifferentiableOperator::AddIntegrator(
                  A->AddSubMatrix(test_vdofs, trial_vdofs, Aee, 1);
               }
            }
-            A->Finalize();
+            // Don't finalize here since multiple callbacks might contribute to the same matrix
+            // A->Finalize() will be called after all callbacks have contributed
         });

         // Create local references for MSVC lambda capture compatibility
         auto& assemble_derivative_sparsematrix_callbacks_ref =
            this->assemble_derivative_sparsematrix_callbacks[derivative_id];

-         assemble_derivative_hypreparmatrix_callbacks[derivative_id].push_back(
+         assemble_derivative_hypreparmatrix_callback[derivative_id] =
            [
               input_is_dependent,
               input_to_field,
@@ -1342,6 +1515,7 @@ void DifferentiableOperator::AddIntegrator(
            {
               f(f_e, spmat);
            }
+            spmat->Finalize();

            if (spmat == nullptr)
            {
@@ -1395,7 +1569,7 @@ void DifferentiableOperator::AddIntegrator(
                       trial_fes->Dof_TrueDof_Matrix());
            }
            delete spmat;
-         });
+         };
      }, derivative_ids);
   }
 }
@@ -90,8 +90,8 @@ void map_quadrature_data_to_fields_impl(
   }
   else
   {
-      MFEM_ABORT_KERNEL("quadrature data mapping to field is not implemented"
-                        " for this field descriptor");
+      MFEM_ABORT("quadrature data mapping to field is not implemented for"
+                 " this field descriptor");
   }
 }

@@ -169,9 +169,8 @@ void map_quadrature_data_to_fields_tensor_impl_1d(
   }
   else
   {
-      MFEM_ABORT_KERNEL("quadrature data mapping to field is not implemented"
-                        "for this field descriptor with sum factorization on"
-                        " tensor product elements");
+      MFEM_ABORT("quadrature data mapping to field is not implemented for"
+                 " this field descriptor with sum factorization on tensor product elements");
   }
 }

@@ -307,9 +306,8 @@ void map_quadrature_data_to_fields_tensor_impl_2d(
   }
   else
   {
-      MFEM_ABORT_KERNEL("quadrature data mapping to field is not implemented"
-                        " for this field descriptor with sum factorization on"
-                        " tensor product elements");
+      MFEM_ABORT("quadrature data mapping to field is not implemented for"
+                 " this field descriptor with sum factorization on tensor product elements");
   }
 }

@@ -494,9 +492,8 @@ void map_quadrature_data_to_fields_tensor_impl_3d(
   }
   else
   {
-      MFEM_ABORT_KERNEL("quadrature data mapping to field is not implemented"
-                        " for this field descriptor with sum factorization on"
-                        " tensor product elements");
+      MFEM_ABORT("quadrature data mapping to field is not implemented for"
+                 " this field descriptor with sum factorization on tensor product elements");
   }
 }

@@ -533,4 +530,65 @@ void map_quadrature_data_to_fields(
   }
 }

+template <size_t N, typename field_operator_ts>
+MFEM_HOST_DEVICE
+void map_quadrature_data_to_fields_conditional(
+   DeviceTensor<2, real_t> &y,
+   const DeviceTensor<3, real_t> &f,
+   const field_operator_ts &fops,
+   const DeviceTensor<1, const real_t> &op_dims,
+   const std::array<DofToQuadMap, N> &dtqmaps,
+   std::array<DeviceTensor<1>, 6> &scratch_mem,
+   const DeviceTensor<3> &fi_shmem,
+   const std::array<bool, N> &conditions,
+   const int &dimension,
+   const bool &use_sum_factorization)
+{
+   int offset = 0;
+   for_constexpr<N>([&](auto i)
+   {
+      if (conditions[i])
+      {
+         [[maybe_unused]] const auto [K, unused, M] = f.GetShape();
+         const int L = static_cast<int>(op_dims(static_cast<size_t>(i)));
+         auto fi = Reshape(&fi_shmem(0, 0, 0), K, L, M);
+         for (int k = 0; k < K; k++)
+         {
+            for (int l = 0; l < L; l++)
+            {
+               for (int m = 0; m < M; m++)
+               {
+                  fi(k, l, m) = f(k, l + offset, m);
+               }
+            }
+         }
+
+         if (use_sum_factorization)
+         {
+            if (dimension == 1)
+            {
+               map_quadrature_data_to_fields_tensor_impl_1d(
+                  y, fi, get<i>(fops), dtqmaps[i], scratch_mem);
+            }
+            else if (dimension == 2)
+            {
+               map_quadrature_data_to_fields_tensor_impl_2d(
+                  y, fi, get<i>(fops), dtqmaps[i], scratch_mem);
+            }
+            else if (dimension == 3)
+            {
+               map_quadrature_data_to_fields_tensor_impl_3d(
+                  y, fi, get<i>(fops), dtqmaps[i], scratch_mem);
+            }
+            else { MFEM_ABORT_KERNEL("dimension not supported"); }
+         }
+         else
+         {
+            map_quadrature_data_to_fields_impl(y, fi, get<i>(fops), dtqmaps[i]);
+         }
+         offset += L;
+      }
+   });
+}
+
 } // namespace mfem::future
@@ -505,13 +505,13 @@ void map_field_to_quadrature_data(
   }
 }

-template <typename field_operator_ts, size_t num_inputs, size_t num_fields>
+template <typename field_operator_ts, size_t N, size_t M>
 MFEM_HOST_DEVICE inline
 void map_fields_to_quadrature_data(
-   std::array<DeviceTensor<2>, num_inputs> &fields_qp,
-   const std::array<DeviceTensor<1>, num_fields> &fields_e,
-   const std::array<DofToQuadMap, num_inputs> &dtqmaps,
-   const std::array<size_t, num_inputs> &input_to_field,
+   std::array<DeviceTensor<2>, N> &fields_qp,
+   const std::array<DeviceTensor<1>, M> &fields_e,
+   const std::array<DofToQuadMap, N> &dtqmaps,
+   const std::array<size_t, N> &input_to_field,
   const field_operator_ts &fops,
   const DeviceTensor<1, const real_t> &integration_weights,
   const std::array<DeviceTensor<1>, 6> &scratch_mem,
@@ -523,7 +523,7 @@ void map_fields_to_quadrature_data(
   // attached to them and we create a dummy field which is not accessed
   // inside the functions it is passed to.
   const auto dummy_field_weight = DeviceTensor<1>(nullptr, 0);
-   for_constexpr<num_inputs>([&](auto i)
+   for_constexpr<N>([&](auto i)
   {
      const DeviceTensor<1> &field_e =
         (input_to_field[i] == SIZE_MAX) ? dummy_field_weight :
@@ -549,12 +549,7 @@ void map_fields_to_quadrature_data(
               fields_qp[i], dtqmaps[i], field_e, get<i>(fops),
               integration_weights, scratch_mem);
         }
-         else
-         {
-#if !(defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP))
-            MFEM_ABORT("unsupported dimension");
-#endif
-         }
+         else { MFEM_ABORT_KERNEL("unsupported dimension"); }
      }
      else
      {
@@ -627,20 +622,20 @@ void map_fields_to_quadrature_data_conditional(
   });
 }

-template <size_t num_inputs, typename field_operator_ts>
+template <size_t N, typename field_operator_ts>
 MFEM_HOST_DEVICE
 void map_direction_to_quadrature_data_conditional(
-   std::array<DeviceTensor<2>, num_inputs> &directions_qp,
+   std::array<DeviceTensor<2>, N> &directions_qp,
   const DeviceTensor<1> &direction_e,
-   const std::array<DofToQuadMap, num_inputs> &dtqmaps,
+   const std::array<DofToQuadMap, N> &dtqmaps,
   field_operator_ts fops,
   const DeviceTensor<1, const real_t> &integration_weights,
   const std::array<DeviceTensor<1>, 6> &scratch_mem,
-   const std::array<bool, num_inputs> &conditions,
+   const std::array<bool, N> &conditions,
   const int &dimension,
   const bool &use_sum_factorization)
 {
-   for_constexpr<num_inputs>([&](auto i)
+   for_constexpr<N>([&](auto i)
   {
      if (conditions[i])
      {
@@ -82,7 +82,7 @@ protected:
 };

 /// @brief Uniform parameter space
-class UniformParameterSpace : public ParameterSpace
+class UniformParameterSpace final : public ParameterSpace
 {
 public:
   /// @brief Constructor for a uniform parameter space
@@ -379,54 +379,102 @@ namespace detail
 /// @param shadow_shmem the shadow shared memory.
 /// @param qpdc the quadrature point data cache holding the resulting
 /// Jacobians on each quadrature point.
-/// @param itod inputs trial operator dimension.
-/// If input is dependent the value corresponds to the spatial dimension, otherwise
-/// a zero indicates non-dependence on the variable.
+/// @param op_dims operator dimensions.
+/// If an operator is dependent, the value corresponds to the spatial dimension.
+/// Otherwise a zero indicates indepence on the variable.
 /// @param q the current quadrature point index.
-template <size_t num_fields>
+/// @param transpose switch to use transpose action.
+template <size_t N>
 MFEM_HOST_DEVICE inline
 void apply_qpdc(
   DeviceTensor<3> &fhat,
-   const std::array<DeviceTensor<2>, num_fields> &shadow_shmem,
+   const std::array<DeviceTensor<2>, N> &shadow_shmem,
   const DeviceTensor<5, const real_t> &qpdc,
-   const DeviceTensor<1, const real_t> &itod,
-   const int &q)
+   const DeviceTensor<1, const real_t> &op_dims,
+   const int &q,
+   bool transpose)
 {
+   const size_t num_ops = op_dims.GetShape()[0];
+
   const int test_vdim = qpdc.GetShape()[0];
   const int test_op_dim = qpdc.GetShape()[1];
   const int trial_vdim = qpdc.GetShape()[2];
-   const int num_qp = qpdc.GetShape()[4];
-   const size_t num_inputs = itod.GetShape()[0];
+   const int total_trial_op_dim = qpdc.GetShape()[3];

-   for (int i = 0; i < test_vdim; i++)
+   const int num_qp = qpdc.GetShape()[4];
+
+   if (transpose)
   {
-      for (int k = 0; k < test_op_dim; k++)
+      for (int j = 0; j < trial_vdim; j++)
      {
-         real_t sum = 0.0;
-         int m_offset = 0;
-         for (size_t s = 0; s < num_inputs; s++)
+         for (int m = 0; m < total_trial_op_dim; m++)
         {
-            const int trial_op_dim = static_cast<int>(itod(s));
-            if (trial_op_dim == 0)
-            {
-               continue;
-            }
-            const auto d_qp =
-               Reshape(&(shadow_shmem[s])[0], trial_vdim, trial_op_dim, num_qp);
-            for (int j = 0; j < trial_vdim; j++)
-            {
-               for (int m = 0; m < trial_op_dim; m++)
-               {
-                  sum += qpdc(i, k, j, m + m_offset, q) * d_qp(j, m, q);
-               }
-            }
-            m_offset += trial_op_dim;
+            fhat(j, m, q) = 0.0;
+         }
+      }
+
+      // Since we don't support more than output space right now
+      // shadow_shmem will always be of size 1.
+      constexpr int shadow_idx_tr = 0;
+      auto d_qp = Reshape(&(shadow_shmem[shadow_idx_tr])[0], test_vdim, test_op_dim,
+                          num_qp);
+
+      int m_offset = 0;
+      for (size_t s = 0; s < num_ops; s++)
+      {
+         const int trial_op_dim = static_cast<int>(op_dims(s));
+         if (trial_op_dim == 0) { continue; }
+
+         for (int j = 0; j < trial_vdim; j++)
+         {
+            for (int m = 0; m < trial_op_dim; m++)
+            {
+               real_t sum = 0.0;
+               for (int i = 0; i < test_vdim; i++)
+               {
+                  for (int k = 0; k < test_op_dim; k++)
+                  {
+                     const real_t contrib = qpdc(i, k, j, m + m_offset, q) * d_qp(i, k, q);
+                     sum += contrib;
+                  }
+               }
+               fhat(j, m + m_offset, q) += sum;
+            }
+         }
+         m_offset += trial_op_dim;
+      }
+   }
+   else
+   {
+      for (int i = 0; i < test_vdim; i++)
+      {
+         for (int k = 0; k < test_op_dim; k++)
+         {
+            real_t sum = 0.0;
+            int m_offset = 0;
+            for (size_t s = 0; s < num_ops; s++)
+            {
+               const int trial_op_dim = static_cast<int>(op_dims(s));
+               if (trial_op_dim == 0) { continue; }
+
+               const auto d_qp =
+                  Reshape(&(shadow_shmem[s])[0], trial_vdim, trial_op_dim, num_qp);
+               for (int j = 0; j < trial_vdim; j++)
+               {
+                  for (int m = 0; m < trial_op_dim; m++)
+                  {
+                     sum += qpdc(i, k, j, m + m_offset, q) * d_qp(j, m, q);
+                  }
+               }
+               m_offset += trial_op_dim;
+            }
+            fhat(i, k, q) = sum;
         }
-         fhat(i, k, q) = sum;
      }
   }
 }
-}
+
+} // namespace detail

 /// @brief Apply the quadrature point data cache (qpdc) to a vector
 /// (usually a direction).
@@ -445,16 +493,18 @@ void apply_qpdc(
 /// @param q1d number of quadrature points in 1D.
 /// @param dimension spatial dimension.
 /// @param use_sum_factorization whether to use sum factorization.
-template <size_t num_fields>
+/// @param T switch to use transpose application.
+template <size_t N>
 MFEM_HOST_DEVICE inline
 void apply_qpdc(
   DeviceTensor<3> &fhat,
-   const std::array<DeviceTensor<2>, num_fields> &shadow_shmem,
+   const std::array<DeviceTensor<2>, N> &shadow_shmem,
   const DeviceTensor<5, const real_t> &qpdc,
   const DeviceTensor<1, const real_t> &itod,
   const int &q1d,
   const int &dimension,
-   const bool &use_sum_factorization)
+   const bool &use_sum_factorization,
+   const bool T = false)
 {
   if (use_sum_factorization)
   {
@@ -462,7 +512,7 @@ void apply_qpdc(
      {
         MFEM_FOREACH_THREAD_DIRECT(q, x, q1d)
         {
-            detail::apply_qpdc(fhat, shadow_shmem, qpdc, itod, q);
+            detail::apply_qpdc(fhat, shadow_shmem, qpdc, itod, q, T);
         }
      }
      else if (dimension == 2)
@@ -472,7 +522,7 @@ void apply_qpdc(
            MFEM_FOREACH_THREAD_DIRECT(qy, y, q1d)
            {
               const int q = qx + q1d * qy;
-               detail::apply_qpdc(fhat, shadow_shmem, qpdc, itod, q);
+               detail::apply_qpdc(fhat, shadow_shmem, qpdc, itod, q, T);
            }
         }
      }
@@ -485,7 +535,7 @@ void apply_qpdc(
               MFEM_FOREACH_THREAD_DIRECT(qz, z, q1d)
               {
                  const int q = qx + q1d * (qy + q1d * qz);
-                  detail::apply_qpdc(fhat, shadow_shmem, qpdc, itod, q);
+                  detail::apply_qpdc(fhat, shadow_shmem, qpdc, itod, q, T);
               }
            }
         }
@@ -500,7 +550,7 @@ void apply_qpdc(
      const int num_qp = qpdc.GetShape()[4];
      MFEM_FOREACH_THREAD_DIRECT(q, x, num_qp)
      {
-         detail::apply_qpdc(fhat, shadow_shmem, qpdc, itod, q);
+         detail::apply_qpdc(fhat, shadow_shmem, qpdc, itod, q, T);
      }
   }
 }
@@ -243,6 +243,35 @@ void process_qf_arg(
   }
 }

+template <typename T, int n>
+MFEM_HOST_DEVICE inline
+void process_qf_arg(
+   const DeviceTensor<1, T> &u,
+   const DeviceTensor<1, T> &v,
+   tensor<T, n> &arg)
+{
+   for (int i = 0; i < n; i++)
+   {
+      arg(i) = u(i);
+   }
+}
+
+template <typename T, int n, int m>
+MFEM_HOST_DEVICE inline
+void process_qf_arg(
+   const DeviceTensor<1, T> &u,
+   const DeviceTensor<1, T> &v,
+   tensor<T, n, m> &arg)
+{
+   for (int i = 0; i < m; i++)
+   {
+      for (int j = 0; j < n; j++)
+      {
+         arg(j, i) = u((i * n) + j);
+      }
+   }
+}
+
 template <typename arg_type>
 MFEM_HOST_DEVICE inline
 void process_qf_arg(const DeviceTensor<2> &u, arg_type &arg, int qp)
@@ -327,20 +356,4 @@ void process_qf_result(
   }
 }

-template <typename T, int n, int m>
-MFEM_HOST_DEVICE inline
-void process_qf_arg(
-   const DeviceTensor<1, T> &u,
-   const DeviceTensor<1, T> &v,
-   tensor<T, n, m> &arg)
-{
-   for (int i = 0; i < m; i++)
-   {
-      for (int j = 0; j < n; j++)
-      {
-         arg(j, i) = u((i * n) + j);
-      }
-   }
-}
-
 } // namespace mfem::future
@@ -1202,7 +1202,14 @@ std::function<void(const Vector&, Vector&)> get_prolongation_transpose(
   const Operator *P = get_prolongation(f);
   auto PT = [=](const Vector &r_local, Vector &y)
   {
-      P->MultTranspose(r_local, y);
+      if (P)
+      {
+         P->MultTranspose(r_local, y);
+      }
+      else
+      {
+         y = r_local;
+      }
   };
   return PT;
 }
@@ -1580,14 +1587,13 @@ struct SharedMemoryInfo
   std::array<int, 6> temp_sizes;
 };

-template <typename entity_t, std::size_t num_fields, std::size_t num_inputs, std::size_t num_outputs, typename input_t>
+template <typename entity_t, std::size_t num_fields, std::size_t num_inputs, std::size_t num_outputs>
 SharedMemoryInfo<num_fields, num_inputs, num_outputs>
 get_shmem_info(
   const std::array<DofToQuadMap, num_inputs> &input_dtq_maps,
   const std::array<DofToQuadMap, num_outputs> &output_dtq_maps,
   const std::vector<FieldDescriptor> &fields,
   const int &num_entities,
-   const input_t &inputs,
   const int &num_qp,
   const std::vector<int> &input_size_on_qp,
   const int &residual_size_on_qp,
@@ -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
@@ -1044,50 +1044,9 @@ void VectorFiniteElement::SetDerivMembers()
   switch (map_type)
   {
      case H_DIV:
-         switch (dim)
-         {
-            case 3: // div: 3D H_DIV -> 3D INTEGRAL
-               deriv_type = DIV;
-               deriv_range_type = SCALAR;
-               deriv_map_type = INTEGRAL;
-               break;
-            case 2: // div: 2D H_DIV -> 2D INTEGRAL
-               deriv_type = DIV;
-               deriv_range_type = SCALAR;
-               deriv_map_type = INTEGRAL;
-               break;
-            default:
-               MFEM_ABORT("Invalid dimension, Dim = " << dim);
-         }
-         break;
-      case H_DIV_R2D:
-         switch (dim)
-         {
-            case 2: // div: 2D H_DIV_R2D -> 2D INTEGRAL
-               deriv_type = DIV;
-               deriv_range_type = SCALAR;
-               deriv_map_type = INTEGRAL;
-               break;
-            case 1: // div: 1D H_DIV_R2D -> 1D INTEGRAL
-               deriv_type = DIV;
-               deriv_range_type = SCALAR;
-               deriv_map_type = INTEGRAL;
-               break;
-            default:
-               MFEM_ABORT("Invalid dimension, Dim = " << dim);
-         }
-         break;
-      case H_DIV_R1D:
-         switch (dim)
-         {
-            case 1: // div: 1D H_DIV_R1D -> 1D INTEGRAL
-               deriv_type = DIV;
-               deriv_range_type = SCALAR;
-               deriv_map_type = INTEGRAL;
-               break;
-            default:
-               MFEM_ABORT("Invalid dimension, Dim = " << dim);
-         }
+         deriv_type = DIV;
+         deriv_range_type = SCALAR;
+         deriv_map_type = INTEGRAL;
         break;
      case H_CURL:
         switch (dim)
@@ -1105,49 +1064,13 @@ void VectorFiniteElement::SetDerivMembers()
               break;
            case 1:
               deriv_type = NONE;
-               deriv_range_type = UNKNOWN_RANGE_TYPE;
-               deriv_map_type = UNKNOWN_MAP_TYPE;
+               deriv_range_type = SCALAR;
+               deriv_map_type = INTEGRAL;
               break;
            default:
               MFEM_ABORT("Invalid dimension, Dim = " << dim);
         }
         break;
-      case H_CURL_R2D:
-         switch (dim)
-         {
-            case 2:
-               // curl: 2D H_CURL_R2D -> H_DIV_R2D
-               deriv_type = CURL;
-               deriv_range_type = VECTOR;
-               deriv_map_type = H_DIV_R2D;
-               break;
-            case 1:
-               // curl: 1D H_CURL_R2D -> H_DIV_R2D
-               deriv_type = CURL;
-               deriv_range_type = VECTOR;
-               deriv_map_type = H_DIV_R2D;
-               break;
-            default:
-               MFEM_ABORT("Invalid dimension, Dim = " << dim);
-         }
-         break;
-      case H_CURL_R1D:
-         switch (dim)
-         {
-            case 1:
-               // curl: 1D H_CURL_R1D -> H_DIV_R1D
-               deriv_type = CURL;
-               deriv_range_type = VECTOR;
-               deriv_map_type = H_DIV_R1D;
-               break;
-            case 0:
-               deriv_type = NONE;
-               deriv_range_type = UNKNOWN_RANGE_TYPE;
-               deriv_map_type = UNKNOWN_MAP_TYPE;
-            default:
-               MFEM_ABORT("Invalid dimension, Dim = " << dim);
-         }
-         break;
      default:
         MFEM_ABORT("Invalid MapType = " << map_type);
   }
@@ -295,20 +295,10 @@ public:
                          $ u(x) = (1/w) \hat u(\hat x) $ */
      H_DIV,     /**< For vector fields; preserves surface integrals of the
                          normal component $ u(x) = (J/w) \hat u(\hat x) $ */
-      H_CURL,    /**< For vector fields; preserves line integrals of the
+      H_CURL     /**< For vector fields; preserves line integrals of the
                          tangential component
                          $ u(x) = J^{-t} \hat u(\hat x) $ (square J),
                          $ u(x) = J(J^t J)^{-1} \hat u(\hat x) $ (general J) */
-      H_DIV_R2D, /**< For 3-component vector fields in 2D; equivalent to a
-                          direct sum of an H_DIV basis and an INTEGRAL basis */
-      H_CURL_R2D,/**< For 3-component vector fields in 2D; equivalent to a
-                          direct sum of an H_CURL basis and a VALUE basis */
-      H_DIV_R1D, /**< For 3-component vector fields in 1D; equivalent to a
-                          direct sum of a VALUE basis and a pair of INTEGRAL
-                          bases */
-      H_CURL_R1D /**< For 3-component vector fields in 1D; equivalent to a
-                          direct sum of an INTEGRAL basis and a pair of VALUE
-                          bases */
   };

   /** @brief Enumeration for DerivType: defines which derivative method
@@ -340,28 +330,12 @@ public:
   int GetDim() const { return dim; }

   /** @brief Returns the vector dimension for vector-valued finite elements,
-       which is also the dimension of the interpolation operation and the
-       width of the DenseMatrix argument in
-       CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape). */
+       which is also the dimension of the interpolation operation. */
   int GetRangeDim() const { return vdim; }

-   /** @brief Returns the vector dimension, in physical space, for
-       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; }
-
-   /** Returns the dimension of the curl for vector-valued finite elements,
-       which is also the width of the DenseMatrix argument in
-       CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape). */
+   /// Returns the dimension of the curl for vector-valued finite elements.
   int GetCurlDim() const { return cdim; }

-   /** Returns the dimension, in physical space, of the curl for vector-valued
-       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; }
-
   /// Returns the Geometry::Type of the reference element.
   Geometry::Type GetGeomType() const { return geom_type; }

@@ -1016,8 +990,6 @@ protected:
 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; }
 };

 /// @brief Class for computing 1D special polynomials and their associated basis
@@ -2531,7 +2531,7 @@ void ND_FuentesPyramidElement::calcCurlBasis(const int p,

 ND_R1D_PointElement::ND_R1D_PointElement(int p)
   : VectorFiniteElement(1, Geometry::POINT, 2, p,
-                         H_CURL_R1D, FunctionSpace::Pk)
+                         H_CURL, FunctionSpace::Pk)
 {
   // VectorFiniteElement::SetDerivMembers doesn't support 0D H_CURL elements
   // so we mimic a 1D element and then correct the dimension here.
@@ -2562,7 +2562,7 @@ ND_R1D_SegmentElement::ND_R1D_SegmentElement(const int p,
                                             const int cb_type,
                                             const int ob_type)
   : VectorFiniteElement(1, Geometry::SEGMENT, 3 * p + 2, p,
-                         H_CURL_R1D, FunctionSpace::Pk),
+                         H_CURL, FunctionSpace::Pk),
     dof2tk(dof),
     cbasis1d(poly1d.GetBasis(p, VerifyClosed(cb_type))),
     obasis1d(poly1d.GetBasis(p - 1, VerifyOpen(ob_type)))
@@ -2839,7 +2839,7 @@ ND_R2D_SegmentElement::ND_R2D_SegmentElement(const int p,
                                             const int cb_type,
                                             const int ob_type)
   : VectorFiniteElement(1, Geometry::SEGMENT, 2 * p + 1, p,
-                         H_CURL_R2D, FunctionSpace::Pk),
+                         H_CURL, FunctionSpace::Pk),
     dof2tk(dof),
     cbasis1d(poly1d.GetBasis(p, VerifyClosed(cb_type))),
     obasis1d(poly1d.GetBasis(p - 1, VerifyOpen(ob_type)))
@@ -3023,7 +3023,7 @@ void ND_R2D_SegmentElement::Project(VectorCoefficient &vc,
 ND_R2D_FiniteElement::ND_R2D_FiniteElement(int p, Geometry::Type G, int Do,
                                           const real_t *tk_fe)
   : VectorFiniteElement(2, G, Do, p,
-                         H_CURL_R2D, FunctionSpace::Pk),
+                         H_CURL, FunctionSpace::Pk),
     tk(tk_fe),
     dof_map(dof),
     dof2tk(dof)
@@ -663,9 +663,6 @@ 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; }
-
   void CalcVShape(const IntegrationPoint &ip,
                   DenseMatrix &shape) const override;

@@ -708,9 +705,6 @@ private:
                           DenseMatrix &I) const;

 public:
-   int GetPhysRangeDim(int space_dim) const override { return 3; }
-   int GetPhysCurlDim(int space_dim) const override { return 3; }
-
   using FiniteElement::CalcVShape;
   using FiniteElement::CalcPhysCurlShape;

@@ -2006,7 +2006,7 @@ RT_R1D_SegmentElement::RT_R1D_SegmentElement(const int p,
                                             const int cb_type,
                                             const int ob_type)
   : VectorFiniteElement(1, Geometry::SEGMENT, 3 * p + 4, p + 1,
-                         H_DIV_R1D, FunctionSpace::Pk),
+                         H_DIV, FunctionSpace::Pk),
     dof2nk(dof),
     cbasis1d(poly1d.GetBasis(p + 1, VerifyClosed(cb_type))),
     obasis1d(poly1d.GetBasis(p, VerifyOpen(ob_type)))
@@ -2281,7 +2281,7 @@ const real_t RT_R2D_SegmentElement::nk[2] = { 0.,1.};
 RT_R2D_SegmentElement::RT_R2D_SegmentElement(const int p,
                                             const int ob_type)
   : VectorFiniteElement(1, Geometry::SEGMENT, p + 1, p + 1,
-                         H_DIV_R2D, FunctionSpace::Pk),
+                         H_DIV, FunctionSpace::Pk),
     dof2nk(dof),
     obasis1d(poly1d.GetBasis(p, VerifyOpen(ob_type)))
 {
@@ -2392,7 +2392,7 @@ void RT_R2D_SegmentElement::LocalInterpolation(const VectorFiniteElement &cfe,
 RT_R2D_FiniteElement::RT_R2D_FiniteElement(int p, Geometry::Type G, int Do,
                                           const real_t *nk_fe)
   : VectorFiniteElement(2, G, Do, p + 1,
-                         H_DIV_R2D, FunctionSpace::Pk),
+                         H_DIV, FunctionSpace::Pk),
     nk(nk_fe),
     dof_map(dof),
     dof2nk(dof)
@@ -510,9 +510,6 @@ 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; }
-
   void CalcVShape(const IntegrationPoint &ip,
                   DenseMatrix &shape) const override;

@@ -550,9 +547,6 @@ private:
                           DenseMatrix &I) const;

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

   void CalcVShape(ElementTransformation &Trans,
@@ -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
@@ -321,17 +321,12 @@ void ParL2FaceRestriction::DoubleValuedConformingMult(
   const int vd = vdim;
   const bool t = byvdim;
   const int threshold = ndofs;
-   const int nsdofs = pfes.GetFaceNbrVSize() / vd;
+   const int nsdofs = pfes.GetFaceNbrVSize();
   auto d_indices1 = scatter_indices1.Read();
   auto d_indices2 = scatter_indices2.Read();
   auto d_x = Reshape(x.Read(), t?vd:ndofs, t?ndofs:vd);
-   const int ne_shared = nsdofs / elem_dofs;
-   const int nedof = elem_dofs;
-   // Note: the shape of face_nbr_data, as determined by
-   // ParFiniteElementSpace::ExchangeFaceNbrData, is (elem_dofs, vdim,
-   // ne_shared), independent of the ordering (byNODES or byVDIM) of the finite
-   // element space.
-   auto d_x_shared = Reshape(face_nbr_data.Read(), elem_dofs, vd, ne_shared);
+   auto d_x_shared = Reshape(face_nbr_data.Read(),
+                             t?vd:nsdofs, t?nsdofs:vd);
   auto d_y = Reshape(y.Write(), nface_dofs, vd, 2, nf);
   mfem::forall(nfdofs, [=] MFEM_HOST_DEVICE (int i)
   {
@@ -351,9 +346,8 @@ void ParL2FaceRestriction::DoubleValuedConformingMult(
         }
         else if (idx2>=threshold) // shared boundary
         {
-            const int e_shared = (idx2 - threshold) / nedof;
-            const int i_shared = (idx2 - threshold) % nedof;
-            d_y(dof, c, 1, face) = d_x_shared(i_shared,c,e_shared);
+            d_y(dof, c, 1, face) = d_x_shared(t?c:(idx2-threshold),
+                                              t?(idx2-threshold):c);
         }
         else // true boundary
         {
@@ -1398,17 +1398,20 @@ void L2FaceRestriction::PermuteAndSetSharedFaceDofsScatterIndices2(
   const int dim = fes.GetMesh()->Dimension();
   const int dof1d = fes.GetTypicalFE()->GetOrder()+1;
   fes.GetTypicalFE()->GetFaceMap(face_id2, face_map);
+   Array<int> face_nbr_dofs;
+   const ParFiniteElementSpace &pfes =
+      static_cast<const ParFiniteElementSpace&>(this->fes);
+   pfes.GetFaceNbrElementVDofs(elem_index, face_nbr_dofs);

   for (int face_dof_elem1 = 0; face_dof_elem1 < face_dofs; ++face_dof_elem1)
   {
      const int face_dof_elem2 = PermuteFaceL2(dim, face_id1, face_id2,
                                               orientation, dof1d, face_dof_elem1);
      const int volume_dof_elem2 = face_map[face_dof_elem2];
-      // Encode the volume DOF index and element index
-      const int global_dof_elem2 = elem_index*elem_dofs + volume_dof_elem2;
+      const int global_dof_elem2 = face_nbr_dofs[volume_dof_elem2];
      const int restriction_dof_elem2 = face_dofs*face_index + face_dof_elem1;
      // Trick to differentiate dof location inter/shared
-      scatter_indices2[restriction_dof_elem2] = ndofs + global_dof_elem2;
+      scatter_indices2[restriction_dof_elem2] = ndofs+global_dof_elem2;
   }
 #endif
 }
@@ -14,7 +14,6 @@

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

 #include <iostream>
 #include <map>
@@ -248,8 +247,7 @@ inline void ArraysByName<T>::Print(std::ostream &os, int width) const
   os << data.size() << '\n';
   for (auto const &it : data)
   {
-      // Note: The method Load() can read any string formatted with std::quoted.
-      os << std::quoted(it.first) << '\n' << it.second.Size() << '\n';
+      os << '"' << it.first << '"' << '\n' << it.second.Size() << '\n';
      it.second.Print(os, width > 0 ? width : it.second.Size());
   }
 }
@@ -260,36 +258,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();
@@ -22,6 +22,14 @@
 //#define _WIN32
 //#define _aligned_malloc(s,a) malloc(s)

+#ifdef NVTX_DEBUG_HPP
+#undef NVTX_COLOR
+#define NVTX_COLOR ::nvtx::kGold
+#include NVTX_DEBUG_HPP
+#else
+#define dbg(...)
+#endif
+
 #ifndef _WIN32
 #include <unistd.h>
 #include <signal.h>
@@ -760,7 +768,12 @@ private:
      switch (mt)
      {
         case MT::HOST_DEBUG:
-            if (GetEnv("MFEM_MMU_STD")) { return new StdHostMemorySpace(); }
+            if (GetEnv("MFEM_MMU_STD"))
+            {
+               dbg("Using STD memory space for debug device!");
+               return new StdHostMemorySpace();
+            }
+            dbg("Using MMU memory space for debug device!");
            return new MmuHostMemorySpace();
 #ifdef MFEM_USE_UMPIRE
         case MT::HOST_UMPIRE:
@@ -200,6 +200,10 @@ void OptionsParser::Parse()
                  isValid = isValidAsInt(argv[i]);
                  *(int *)(options[j].var_ptr) = atoi(argv[i++]);
                  break;
+	       case LONG:
+		  isValid = isValidAsInt(argv[i]);
+                  *(long long *)(options[j].var_ptr) = atoi(argv[i++]);
+                  break;
               case DOUBLE:
                  isValid = isValidAsDouble(argv[i]);
                  *(real_t *)(options[j].var_ptr) = atof(argv[i++]);
@@ -278,6 +282,10 @@ void OptionsParser::WriteValue(const Option &opt, std::ostream &os)
      case INT:
         os << *(int *)(opt.var_ptr);
         break;
+      
+      case LONG:
+         os << *(long long *)(opt.var_ptr);
+         break;	 

      case DOUBLE:
         os << *(real_t *)(opt.var_ptr);
@@ -31,7 +31,7 @@ class Vector;
 class OptionsParser
 {
 public:
-   enum OptionType { INT, DOUBLE, STRING, STD_STRING, ENABLE, DISABLE, ARRAY, VECTOR };
+   enum OptionType { INT, DOUBLE, STRING, STD_STRING, ENABLE, DISABLE, ARRAY, VECTOR , LONG };

 private:
   struct Option
@@ -98,6 +98,14 @@ public:
                            required));
   }

+   /// Add a long integer option and set 'var' to receive the value.
+   void AddOption(long long *var, const char *short_name, const char *long_name,
+                  const char *description, bool required = false)
+   {
+      options.Append(Option(LONG, var, short_name, long_name, description,
+                            required));
+   }
+
   /// Add a double option and set 'var' to receive the value.
   void AddOption(real_t *var, const char *short_name, const char *long_name,
                  const char *description, bool required = false)
@@ -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
@@ -317,9 +317,6 @@ void HypreParVector::WrapHypreParVector(hypre_ParVector *y, bool owner)

 Vector * HypreParVector::GlobalVector() const
 {
-   MFEM_VERIFY(size > 0,
-               "GlobalVector method can only be called on vectors wherein each "
-               "process owns one or more entries");
   hypre_Vector *hv = hypre_ParVectorToVectorAll(*this);
   Vector *v = new Vector(hv->data, internal::to_int(hv->size));
   v->MakeDataOwner();
@@ -6556,7 +6553,7 @@ HypreLOBPCG::SetPreconditioner(Solver & precond)
 }

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

@@ -6603,7 +6600,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 +6621,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 +6863,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 +6921,7 @@ HypreAME::createDummyVectors() const
   }
 }

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

-Vector **
+HypreParVector **
 HypreAME::StealEigenvectors()
 {
   if ( eigenvectors == NULL )
@@ -6962,11 +6941,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
 {
@@ -38,13 +38,6 @@
 #if PETSC_VERSION_LT(3,19,0)
 #define PETSC_SUCCESS 0
 #endif
-#if PETSC_VERSION_LT(3,23,0)
-#define PetscContainerSetCtxDestroy(A,B) PetscContainerSetUserDestroy(A,B)
-typedef PetscErrorCode (PetscCtxDestroyFn)(void**);
-#endif
-#if PETSC_VERSION_LT(3,24,0)
-typedef PetscErrorCode KSPMonitorFn(KSP,PetscInt,PetscReal,void*);
-#endif

 #include <fstream>
 #include <iomanip>
@@ -84,17 +77,13 @@ static PetscErrorCode __mfem_mat_shell_apply_transpose(Mat,Vec,Vec);
 static PetscErrorCode __mfem_mat_shell_destroy(Mat);
 static PetscErrorCode __mfem_mat_shell_copy(Mat,Mat,MatStructure);
 #if PETSC_VERSION_LT(3,23,0)
-typedef void *PetscCtxRt;
-#elif PETSC_VERSION_LT(3,25,0)
-typedef void **PetscCtxRt;
-#endif
-static PetscErrorCode __mfem_array_container_destroy(PetscCtxRt);
-static PetscErrorCode __mfem_matarray_container_destroy(PetscCtxRt);
-#if PETSC_VERSION_LT(3,23,0)
-static PetscErrorCode __mfem_monitor_ctx_destroy(void**);
+static PetscErrorCode __mfem_array_container_destroy(void*);
+static PetscErrorCode __mfem_matarray_container_destroy(void *);
 #else
-static PetscErrorCode __mfem_monitor_ctx_destroy(PetscCtxRt);
+static PetscErrorCode __mfem_array_container_destroy(void**);
+static PetscErrorCode __mfem_matarray_container_destroy(void**);
 #endif
+static PetscErrorCode __mfem_monitor_ctx_destroy(void**);

 // auxiliary functions
 static PetscErrorCode Convert_Array_IS(MPI_Comm,bool,const mfem::Array<int>*,
@@ -1328,7 +1317,11 @@ BlockDiagonalConstructor(MPI_Comm comm,

         ierr = PetscContainerCreate(comm,&c); CCHKERRQ(comm,ierr);
         ierr = PetscContainerSetPointer(c,ptrs[i]); CCHKERRQ(comm,ierr);
+#if PETSC_VERSION_LT(3,23,0)
+         ierr = PetscContainerSetUserDestroy(c,__mfem_array_container_destroy);
+#else
         ierr = PetscContainerSetCtxDestroy(c,__mfem_array_container_destroy);
+#endif
         CCHKERRQ(comm,ierr);
         ierr = PetscObjectCompose((PetscObject)A,names[i],(PetscObject)c);
         CCHKERRQ(comm,ierr);
@@ -1655,7 +1648,11 @@ void PetscParMatrix::ConvertOperator(MPI_Comm comm, const Operator &op, Mat* A,
         PetscContainer c;
         ierr = PetscContainerCreate(comm,&c); CCHKERRQ(comm,ierr);
         ierr = PetscContainerSetPointer(c,vmatsl2l); PCHKERRQ(c,ierr);
+#if PETSC_VERSION_LT(3,23,0)
+         ierr = PetscContainerSetUserDestroy(c,__mfem_matarray_container_destroy);
+#else
         ierr = PetscContainerSetCtxDestroy(c,__mfem_matarray_container_destroy);
+#endif
         PCHKERRQ(c,ierr);
         ierr = PetscObjectCompose((PetscObject)(*A),"_MatIS_PtAP_l2l",(PetscObject)c);
         PCHKERRQ((*A),ierr);
@@ -1751,7 +1748,11 @@ void PetscParMatrix::ConvertOperator(MPI_Comm comm, const Operator &op, Mat* A,

            ierr = PetscContainerCreate(PETSC_COMM_SELF,&c); PCHKERRQ(B,ierr);
            ierr = PetscContainerSetPointer(c,ptrs[i]); PCHKERRQ(B,ierr);
+#if PETSC_VERSION_LT(3,23,0)
+            ierr = PetscContainerSetUserDestroy(c,__mfem_array_container_destroy);
+#else
            ierr = PetscContainerSetCtxDestroy(c,__mfem_array_container_destroy);
+#endif
            PCHKERRQ(B,ierr);
            ierr = PetscObjectCompose((PetscObject)(B),names[i],(PetscObject)c);
            PCHKERRQ(B,ierr);
@@ -2197,7 +2198,11 @@ PetscParMatrix * RAP(PetscParMatrix *Rt, PetscParMatrix *A, PetscParMatrix *P)
         ierr = PetscContainerCreate(PetscObjectComm((PetscObject)B),&c);
         PCHKERRQ(B,ierr);
         ierr = PetscContainerSetPointer(c,vmatsl2l); PCHKERRQ(c,ierr);
+#if PETSC_VERSION_LT(3,23,0)
+         ierr = PetscContainerSetUserDestroy(c,__mfem_matarray_container_destroy);
+#else
         ierr = PetscContainerSetCtxDestroy(c,__mfem_matarray_container_destroy);
+#endif
         PCHKERRQ(c,ierr);
         ierr = PetscObjectCompose((PetscObject)B,"_MatIS_PtAP_l2l",(PetscObject)c);
         PCHKERRQ(B,ierr);
@@ -2480,6 +2485,7 @@ void PetscSolver::SetMaxIter(int max_iter)

 void PetscSolver::SetPrintLevel(int plev)
 {
+   typedef PetscErrorCode (*myPetscFunc)(void**);
   PetscViewerAndFormat *vf = NULL;
   PetscViewer viewer = PETSC_VIEWER_STDOUT_(PetscObjectComm(obj));

@@ -2492,6 +2498,7 @@ void PetscSolver::SetPrintLevel(int plev)
   {
      // there are many other options, see the function KSPSetFromOptions() in
      // src/ksp/ksp/interface/itcl.c
+      typedef PetscErrorCode (*myMonitor)(KSP,PetscInt,PetscReal,void*);
      KSP ksp = (KSP)obj;
      if (plev >= 0)
      {
@@ -2500,29 +2507,29 @@ void PetscSolver::SetPrintLevel(int plev)
      if (plev == 1)
      {
 #if PETSC_VERSION_LT(3,15,0)
-         ierr = KSPMonitorSet(ksp,(KSPMonitorFn *)KSPMonitorDefault,vf,
+         ierr = KSPMonitorSet(ksp,(myMonitor)KSPMonitorDefault,vf,
 #else
-         ierr = KSPMonitorSet(ksp,(KSPMonitorFn *)KSPMonitorResidual,vf,
+         ierr = KSPMonitorSet(ksp,(myMonitor)KSPMonitorResidual,vf,
 #endif
-                              (PetscCtxDestroyFn *)PetscViewerAndFormatDestroy);
+                              (myPetscFunc)PetscViewerAndFormatDestroy);
         PCHKERRQ(ksp,ierr);
      }
      else if (plev > 1)
      {
         ierr = KSPSetComputeSingularValues(ksp,PETSC_TRUE); PCHKERRQ(ksp,ierr);
-         ierr = KSPMonitorSet(ksp,(KSPMonitorFn *)KSPMonitorSingularValue,vf,
-                              (PetscCtxDestroyFn *)PetscViewerAndFormatDestroy);
+         ierr = KSPMonitorSet(ksp,(myMonitor)KSPMonitorSingularValue,vf,
+                              (myPetscFunc)PetscViewerAndFormatDestroy);
         PCHKERRQ(ksp,ierr);
         if (plev > 2)
         {
            ierr = PetscViewerAndFormatCreate(viewer,PETSC_VIEWER_DEFAULT,&vf);
            PCHKERRQ(viewer,ierr);
 #if PETSC_VERSION_LT(3,15,0)
-            ierr = KSPMonitorSet(ksp,(KSPMonitorFn *)KSPMonitorTrueResidualNorm,vf,
+            ierr = KSPMonitorSet(ksp,(myMonitor)KSPMonitorTrueResidualNorm,vf,
 #else
-            ierr = KSPMonitorSet(ksp,(KSPMonitorFn *)KSPMonitorTrueResidual,vf,
+            ierr = KSPMonitorSet(ksp,(myMonitor)KSPMonitorTrueResidual,vf,
 #endif
-                                 (PetscCtxDestroyFn *)PetscViewerAndFormatDestroy);
+                                 (myPetscFunc)PetscViewerAndFormatDestroy);
            PCHKERRQ(ksp,ierr);
         }
      }
@@ -2538,7 +2545,7 @@ void PetscSolver::SetPrintLevel(int plev)
      if (plev > 0)
      {
         ierr = SNESMonitorSet(snes,(myMonitor)SNESMonitorDefault,vf,
-                               (PetscCtxDestroyFn *)PetscViewerAndFormatDestroy);
+                               (myPetscFunc)PetscViewerAndFormatDestroy);
         PCHKERRQ(snes,ierr);
      }
   }
@@ -4156,31 +4163,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(); }
 }
@@ -5333,27 +5329,21 @@ static PetscErrorCode __mfem_pc_shell_destroy(PC pc)
   PetscFunctionReturn(PETSC_SUCCESS);
 }

-static PetscErrorCode __mfem_array_container_destroy(PetscCtxRt ptr)
+#if PETSC_VERSION_LT(3,23,0)
+
+static PetscErrorCode __mfem_array_container_destroy(void *ptr)
 {
   PetscErrorCode ierr;

   PetscFunctionBeginUser;
-#if PETSC_VERSION_LT(3,23,0)
   ierr = PetscFree(ptr); CHKERRQ(ierr);
-#else
-   ierr = PetscFree(*(void**)ptr); CHKERRQ(ierr);
-#endif
   PetscFunctionReturn(PETSC_SUCCESS);
 }

-static PetscErrorCode __mfem_matarray_container_destroy(PetscCtxRt ptr)
+static PetscErrorCode __mfem_matarray_container_destroy(void *ptr)
 {
-#if PETSC_VERSION_LT(3,23,0)
   mfem::Array<Mat> *a = (mfem::Array<Mat>*)ptr;
-#else
-   mfem::Array<Mat> *a = *(mfem::Array<Mat>**)ptr;
-#endif
-   PetscErrorCode ierr;
+   PetscErrorCode   ierr;

   PetscFunctionBeginUser;
   for (int i=0; i<a->Size(); i++)
@@ -5366,16 +5356,41 @@ static PetscErrorCode __mfem_matarray_container_destroy(PetscCtxRt ptr)
   PetscFunctionReturn(PETSC_SUCCESS);
 }

-#if PETSC_VERSION_LT(3,23,0)
-static PetscErrorCode __mfem_monitor_ctx_destroy(void **ctx)
 #else
-static PetscErrorCode __mfem_monitor_ctx_destroy(PetscCtxRt ctx)
-#endif
+
+static PetscErrorCode __mfem_array_container_destroy(void **ptr)
 {
   PetscErrorCode ierr;

   PetscFunctionBeginUser;
-   ierr = PetscFree(*(void**)ctx); CHKERRQ(ierr);
+   ierr = PetscFree(*ptr); CHKERRQ(ierr);
+   PetscFunctionReturn(PETSC_SUCCESS);
+}
+
+static PetscErrorCode __mfem_matarray_container_destroy(void **ptr)
+{
+   mfem::Array<Mat> *a = (mfem::Array<Mat>*)*ptr;
+   PetscErrorCode   ierr;
+
+   PetscFunctionBeginUser;
+   for (int i=0; i<a->Size(); i++)
+   {
+      Mat M = (*a)[i];
+      MPI_Comm comm = PetscObjectComm((PetscObject)M);
+      ierr = MatDestroy(&M); CCHKERRQ(comm,ierr);
+   }
+   delete a;
+   PetscFunctionReturn(PETSC_SUCCESS);
+}
+
+#endif
+
+static PetscErrorCode __mfem_monitor_ctx_destroy(void **ctx)
+{
+   PetscErrorCode  ierr;
+
+   PetscFunctionBeginUser;
+   ierr = PetscFree(*ctx); CHKERRQ(ierr);
   PetscFunctionReturn(PETSC_SUCCESS);
 }

@@ -5620,7 +5635,11 @@ static PetscErrorCode MatConvert_hypreParCSR_AIJ(hypre_ParCSRMatrix* hA,Mat* pA)

      ierr = PetscContainerCreate(comm,&c); CHKERRQ(ierr);
      ierr = PetscContainerSetPointer(c,ptrs[i]); CHKERRQ(ierr);
+#if PETSC_VERSION_LT(3,23,0)
+      ierr = PetscContainerSetUserDestroy(c,__mfem_array_container_destroy);
+#else
      ierr = PetscContainerSetCtxDestroy(c,__mfem_array_container_destroy);
+#endif
      CHKERRQ(ierr);
      ierr = PetscObjectCompose((PetscObject)(*pA),names[i],(PetscObject)c);
      CHKERRQ(ierr);
@@ -5714,7 +5733,11 @@ static PetscErrorCode MatConvert_hypreParCSR_IS(hypre_ParCSRMatrix* hA,Mat* pA)

      ierr = PetscContainerCreate(PETSC_COMM_SELF,&c); CHKERRQ(ierr);
      ierr = PetscContainerSetPointer(c,ptrs[i]); CHKERRQ(ierr);
+#if PETSC_VERSION_LT(3,23,0)
+      ierr = PetscContainerSetUserDestroy(c,__mfem_array_container_destroy);
+#else
      ierr = PetscContainerSetCtxDestroy(c,__mfem_array_container_destroy);
+#endif
      CHKERRQ(ierr);
      ierr = PetscObjectCompose((PetscObject)lA,names[i],(PetscObject)c);
      CHKERRQ(ierr);
@@ -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
 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)

 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);
@@ -31,11 +31,12 @@ add_subdirectory(hdiv-linear-solver)
 add_subdirectory(hooke)
 add_subdirectory(meshing)
 add_subdirectory(mtop)
+add_subdirectory(mtop/chpt)
+add_subdirectory(mtop/examples)
 add_subdirectory(multidomain)
 add_subdirectory(nurbs)
 add_subdirectory(parelag)
 add_subdirectory(performance)
-add_subdirectory(plasma)
 add_subdirectory(shifted)
 add_subdirectory(solvers)
 add_subdirectory(spde)
@@ -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

--- a/Show More
+++ b/Show More
Author	SHA1	Message	Date
camierjs	9bef1cc28f	Cleanup pass	2026-03-01 15:30:39 -08:00
camierjs	6b45d9a69b	Merge branch 'master' into checkpointing_examples_gpu	2026-03-01 14:38:23 -08:00
camierjs	db32ae9807	Cleanup	2026-03-01 14:34:24 -08:00
camierjs	545a257ab5	LE style, compile fix	2026-02-10 12:43:39 -08:00
blaz	4b4aff736d	adjoint struct	2026-02-09 00:23:04 -08:00
blaz	727339e0ae	check for objective gradients	2026-02-08 22:07:05 -08:00
blaz	d8d4d5b874	Merge branch 'dfem-transpose' into checkpointing_examples	2026-02-08 20:24:10 -08:00
blaz	5a7a017a4e	objective integradn	2026-02-08 20:21:33 -08:00
blaz	224aff5ac9	forward part of the adjoint computations	2026-02-07 21:09:30 -08:00
blaz	e42762365f	adjoint initial version	2026-02-07 20:09:35 -08:00
blaz	765fb4a608	added spectral matrix preconditioner	2026-02-07 19:38:59 -08:00
John Camier	591156be2f	Merge branch 'master' into dfem-transpose	2026-01-30 18:42:29 -08:00
blaz	8e2478f964	Checkpointing examples	2026-01-25 19:45:26 -08:00
blaz	06c485bb3e	updates	2026-01-21 11:03:43 -08:00
bslazarov	1cd27605f0	memory transfer to Device for q-functions	2026-01-20 14:43:06 -08:00
blaz	547fd02911	anisotrpic elasticity structure	2026-01-20 09:57:29 -08:00
blaz	aa92498033	explicit linear elasticity	2026-01-20 09:51:15 -08:00
blaz	024576620a	initial REVOLVE	2026-01-16 11:29:25 -08:00
blaz	5b670bda3f	file and segement storage	2026-01-15 22:59:32 -08:00
blaz	6b0e135bdb	rename files	2026-01-15 22:18:33 -08:00
blaz	204b9df2d6	checkpinting into file - multiple files - one file per snapshot	2026-01-15 22:16:31 -08:00
blaz	8b49379f87	fixed slot in-memory or file storage for checkpointing	2026-01-15 22:15:09 -08:00
blaz	d32a00fb2f	CMake list	2026-01-15 19:57:10 -08:00
blaz	59c9ca70e4	initial version	2026-01-15 19:49:02 -08:00
camierjs	9646ce4d1a	Fix dFEM linearized sections leaks	2026-01-12 17:10:40 -08:00
Julian Andrej	f46a81bf3a	clearer comment	2026-01-12 13:09:39 -08:00
Julian Andrej	460e5ba483	correct sparsematrix test	2026-01-12 13:09:04 -08:00
Julian Andrej	34593dccac	properly initialize matrix pointer	2026-01-09 16:37:16 -08:00
Julian Andrej	248ab78e77	missing dual types	2026-01-09 16:04:01 -08:00
Julian Andrej	6b9bdb11df	fix hypreparmatrix assemble with multiple integrators	2026-01-09 13:18:37 -08:00
Julian Andrej	505d2e5123	fix multiple integrators	2026-01-09 12:59:25 -08:00
Julian Andrej	47cff96e2e	documentation	2026-01-08 10:57:34 -08:00
Julian Andrej	010f455382	documentation	2026-01-06 08:51:46 -08:00
Julian Andrej	acebc322b8	cleanup	2026-01-06 08:43:02 -08:00
Julian Andrej	25ec78df9a	make tests more strict	2026-01-05 16:21:58 -08:00
Julian Andrej	8c78387136	remove periodic tests for transposes	2025-12-30 13:40:59 -08:00
Julian Andrej	8684d9d5d1	working	2025-12-30 12:41:29 -08:00
Julian Andrej	2375953e95	debugging state	2025-12-29 13:35:24 -08:00
Julian Andrej	ecbb196551	bug	2025-12-23 08:27:01 -08:00
Julian Andrej	90ca758a98	again	2025-12-23 07:40:55 -08:00
Julian Andrej	93d674a7a2	argument order	2025-12-23 07:31:27 -08:00
Julian Andrej	03ec3f739b	correct size in parallel for ldof vectors	2025-12-23 07:27:36 -08:00
Julian Andrej	a3eb9c8e6b	doxygen	2025-12-22 16:30:39 -08:00
Julian Andrej	30f3e123cd	more warnings	2025-12-22 16:27:42 -08:00
Julian Andrej	4d4c2325f9	warnings	2025-12-22 16:25:39 -08:00
Julian Andrej	db84abedfc	make dual numbers work	2025-12-22 16:21:48 -08:00
Julian Andrej	2609e3af6e	transpose action using cache	2025-12-22 09:05:50 -08:00