mod

.github/actions/sanitize/mfem/action.yml

@@ -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:

.github/actions/sanitize/restore/action.yml

@@ -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:

.github/actions/sanitize/setup/action.yml

@@ -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}}

.github/workflows/build-container.yml

@@ -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

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

@@ -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

.github/workflows/codeql-analysis.yml

@@ -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

.github/workflows/mfem-analysis.yml

@@ -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

.github/workflows/repo-check.yml

@@ -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
 

.github/workflows/sanitize-build-hypre.yml

@@ -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

.github/workflows/sanitize-build-libcxx.yml

@@ -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}}

.github/workflows/sanitize-build-lsan.yml

@@ -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

.github/workflows/sanitize-build-metis.yml

@@ -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

.github/workflows/sanitize-tests.yml

@@ -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:

.gitignore

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

.gitlab/jobs/dane.yml

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

.gitlab/jobs/matrix.yml

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

.gitlab/scripts/baseline

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

CHANGELOG

@@ -11,20 +11,8 @@
 Version 4.9.1 (development)
 ===========================
 
-- Policy for AI-assisted contribution added to CONTRIBUTING.md
-
 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.
-
 - Improved the gridfunction projection routines. Projections work for Scalar,
   Vector and VectorFE, also NURBS versions. Optionally different types of
   projections can be selected, default behaviour has not changed.

CONTRIBUTING.md

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

config/cmake/MFEMConfig.cmake.in

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

config/cmake/modules/FindUMPIRE.cmake

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

examples/CMakeLists.txt

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

examples/amgx/makefile

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

examples/caliper/makefile

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

examples/ex12p.cpp

@@ -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

examples/ex22p.cpp

@@ -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)
    {

examples/ex37.cpp

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

examples/ex37.hpp

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

examples/ex37p.cpp

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

examples/ginkgo/makefile

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

examples/makefile

@@ -71,7 +71,6 @@ endif
 
 SUBDIRS_ALL = $(addsuffix /all,$(SUBDIRS))
 SUBDIRS_TEST = $(addsuffix /test,$(SUBDIRS))
-SUBDIRS_TEST_NOCLEAN = $(addsuffix /test-noclean,$(SUBDIRS))
 SUBDIRS_CLEAN = $(addsuffix /clean,$(SUBDIRS))
 SUBDIRS_TPRINT = $(addsuffix /test-print,$(SUBDIRS))
 
@@ -88,9 +87,8 @@ SUBDIRS_TPRINT = $(addsuffix /test-print,$(SUBDIRS))
 
 all: $(EXAMPLES) $(SUBDIRS_ALL)
 
-.PHONY: $(SUBDIRS_ALL) $(SUBDIRS_TEST) $(SUBDIRS_TEST_NOCLEAN) \
-   $(SUBDIRS_CLEAN) $(SUBDIRS_TPRINT)
-$(SUBDIRS_ALL) $(SUBDIRS_TEST) $(SUBDIRS_TEST_NOCLEAN) $(SUBDIRS_CLEAN):
+.PHONY: $(SUBDIRS_ALL) $(SUBDIRS_TEST) $(SUBDIRS_CLEAN) $(SUBDIRS_TPRINT)
+$(SUBDIRS_ALL) $(SUBDIRS_TEST) $(SUBDIRS_CLEAN):
 	$(MAKE) -C $(@D) $(@F)
 $(SUBDIRS_TPRINT):
 	@$(MAKE) -C $(@D) $(@F)
@@ -109,7 +107,6 @@ endif
 MFEM_TESTS = EXAMPLES
 include $(MFEM_TEST_MK)
 test: $(SUBDIRS_TEST)
-test-noclean: $(SUBDIRS_TEST_NOCLEAN)
 test-print: $(SUBDIRS_TPRINT)
 
 # Testing: Parallel vs. serial runs
@@ -160,8 +157,6 @@ ex37-test-seq: ex37
 	@$(call mfem-test,$<,, Serial example,-mi 3)
 ex37p-test-par: ex37p
 	@$(call mfem-test,$<, $(RUN_MPI), Parallel example,-mi 3)
-ex39-test-seq: ex39
-	@$(call mfem-test,$<,, Serial example,-m ../data/compass.mesh)
 ex41-test-seq: ex41
 	@$(call mfem-test,$<,, Serial example,-tf 1.0)
 ex41p-test-par: ex41p

fem/bilinearform.cpp

@@ -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);
             }
          }

fem/bilininteg.hpp

@@ -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:

fem/bounds.cpp

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

fem/bounds.hpp

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

fem/coefficient.hpp

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

fem/datacollection.cpp

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

fem/dfem/doperator.hpp

@@ -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);
    }
 }

fem/dfem/integrate.hpp

@@ -533,4 +533,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

fem/dfem/interpolate.hpp

@@ -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])
       {

fem/dfem/qfunction_apply.hpp

@@ -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);
       }
    }
 }

fem/dfem/qfunction_transform.hpp

@@ -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

fem/dfem/util.hpp

@@ -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,

fem/dgmassinv_kernels.hpp

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

fem/estimators.hpp

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

fem/fe/fe_base.hpp

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

fem/fe/fe_nd.hpp

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

fem/fe/fe_rt.hpp

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

fem/fespace.cpp

@@ -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
 {

fem/fespace.hpp

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

fem/gridfunc.cpp

@@ -345,6 +345,27 @@ void GridFunction::ComputeFlux(BilinearFormIntegrator &blfi,
    }
 }
 
+int GridFunction::VectorDim() const
+{
+   const FiniteElement *fe = fes->GetTypicalFE();
+   if (!fe || fe->GetRangeType() == FiniteElement::SCALAR)
+   {
+      return fes->GetVDim();
+   }
+   return fes->GetVDim()*std::max(fes->GetMesh()->SpaceDimension(),
+                                  fe->GetRangeDim());
+}
+
+int GridFunction::CurlDim() const
+{
+   const FiniteElement *fe = fes->GetTypicalFE();
+   if (!fe || fe->GetRangeType() == FiniteElement::SCALAR)
+   {
+      return 2 * fes->GetMesh()->SpaceDimension() - 3;
+   }
+   return fes->GetVDim()*fe->GetCurlDim();
+}
+
 void GridFunction::GetTrueDofs(Vector &tv) const
 {
    const SparseMatrix *R = fes->GetRestrictionMatrix();
@@ -2029,18 +2050,6 @@ void GridFunction::AccumulateAndCountBdrValues(
    Coefficient *coeff[], VectorCoefficient *vcoeff, const Array<int> &attr,
    Array<int> &values_counter)
 {
-   if (vcoeff)
-   {
-      MFEM_VERIFY(fes->GetVDim() == vcoeff->GetVDim(),
-                  "vcoeff vdim != fes VDim");
-      MFEM_VERIFY(fes->GetTypicalBE()->GetMapType() == FiniteElement::VALUE &&
-                  fes->GetTypicalBE()->GetRangeType() ==
-                  FiniteElement::SCALAR,
-                  "Can only call ProjectBdrCoefficient on scalar value-type "
-                  "boundary elements. "
-                  "Did you intended to call ProjectBdrCoefficientNormal or "
-                  "ProjectBdrCoefficientTangent for vector finite elements?");
-   }
    Array<int> vdofs;
    Vector vc;
 
@@ -2193,9 +2202,6 @@ void GridFunction::AccumulateAndCountBdrTangentValues(
    VectorCoefficient &vcoeff, const Array<int> &bdr_attr,
    Array<int> &values_counter)
 {
-   MFEM_VERIFY(fes->GetTypicalBE()->GetPhysRangeDim(
-                  fes->GetMesh()->SpaceDimension()) == vcoeff.GetVDim(),
-               "vcoeff vdim != PhysRangeDim");
    const FiniteElement *fe;
    ElementTransformation *T;
    Array<int> dofs;
@@ -2349,9 +2355,6 @@ void GridFunction::ProjectDeltaCoefficient(DeltaCoefficient &delta_coeff,
 
 void GridFunction::ProjectCoefficient(Coefficient &coeff, ProjectType type)
 {
-   MFEM_VERIFY(
-      VectorDim() == 1,
-      "Cannot project scalar Coefficient onto vector GridFunction");
    DeltaCoefficient *delta_c = dynamic_cast<DeltaCoefficient *>(&coeff);
    DofTransformation doftrans;
    Array<int> vdofs;
@@ -2627,7 +2630,6 @@ void GridFunction::ProjectCoefficient(
 void GridFunction::ProjectCoefficient(VectorCoefficient &vcoeff,
                                       ProjectType type)
 {
-   MFEM_VERIFY(VectorDim() == vcoeff.GetVDim(), "vcoeff vdim != VectorDim()");
    Array<int> vdofs;
    Vector vals;
    DofTransformation doftrans;
@@ -2943,7 +2945,6 @@ void GridFunction::ProjectCoefficientElementL2(VectorCoefficient &vcoeff)
 void GridFunction::ProjectCoefficient(
    VectorCoefficient &vcoeff, Array<int> &dofs)
 {
-   MFEM_VERIFY(VectorDim() == vcoeff.GetVDim(), "vcoeff vdim != VectorDim()");
    int el = -1;
    ElementTransformation *T = NULL;
    const FiniteElement *fe = NULL;
@@ -2973,7 +2974,6 @@ void GridFunction::ProjectCoefficient(
 
 void GridFunction::ProjectCoefficient(VectorCoefficient &vcoeff, int attribute)
 {
-   MFEM_VERIFY(VectorDim() == vcoeff.GetVDim(), "vcoeff vdim != VectorDim()");
    int i;
    Array<int> vdofs;
    Vector vals;
@@ -3030,14 +3030,9 @@ void GridFunction::ProjectCoefficient(Coefficient *coeff[])
    }
 }
 
-void GridFunction::ProjectDiscCoefficient(
-   std::variant<Coefficient*, VectorCoefficient*> coeff, Array<int> &dof_attr)
+void GridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff,
+                                          Array<int> &dof_attr)
 {
-   std::visit([&](auto* c)
-   {
-      MFEM_VERIFY(VectorDim() == c->GetVDim(), "coeff vdim != VectorDim()");
-   }, coeff);
-
    Array<int> vdofs;
    Vector vals;
 
@@ -3051,10 +3046,7 @@ void GridFunction::ProjectDiscCoefficient(
    {
       fes->GetElementVDofs(i, vdofs);
       vals.SetSize(vdofs.Size());
-      std::visit([&](auto* c)
-      {
-         fes->GetFE(i)->Project(*c, *fes->GetElementTransformation(i), vals);
-      }, coeff);
+      fes->GetFE(i)->Project(coeff, *fes->GetElementTransformation(i), vals);
 
       // the values in shared dofs are determined from the element with maximal
       // attribute
@@ -3070,15 +3062,17 @@ void GridFunction::ProjectDiscCoefficient(
    }
 }
 
+void GridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff)
+{
+   Array<int> dof_attr;
+   ProjectDiscCoefficient(coeff, dof_attr);
+}
+
 void GridFunction::ProjectDiscCoefficient(Coefficient &coeff, AvgType type)
 {
    // Harmonic  (x1 ... xn) = [ (1/x1 + ... + 1/xn) / n ]^-1.
    // Arithmetic(x1 ... xn) = (x1 + ... + xn) / n.
 
-   MFEM_VERIFY(
-      VectorDim() == 1,
-      "Cannot project a scalar coefficient onto a vector GridFunction");
-
    Array<int> zones_per_vdof;
    AccumulateAndCountZones(coeff, type, zones_per_vdof);
 
@@ -3088,7 +3082,6 @@ void GridFunction::ProjectDiscCoefficient(Coefficient &coeff, AvgType type)
 void GridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff,
                                           AvgType type)
 {
-   MFEM_VERIFY(VectorDim() == coeff.GetVDim(), "coeff vdim != VectorDim()");
    Array<int> zones_per_vdof;
    AccumulateAndCountZones(coeff, type, zones_per_vdof);
 
@@ -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);
@@ -5252,30 +5242,6 @@ void GridFunction::GetElementBounds(const PLBound &plb,
                                     Vector &lower, Vector &upper,
                                     const int vdim) const
 {
-   if (UseDevice() && Device::Allows(Backend::DEVICE_MASK) &&
-       plb.GetBasisType() != BasisType::Positive &&
-       UsesTensorBasis(*fes))
-   {
-      const FiniteElement &fe = *fes->GetTypicalFE();
-      const int rdim = fe.GetDim();
-      const int fes_dim = fes->GetVDim();
-      const int nel = fes->GetNE();
-      const int nd = fe.GetDof();
-
-      Vector e_vec(nd*fes_dim*nel, Device::GetDeviceMemoryType());
-      e_vec.UseDevice(true);
-      const ElementRestrictionOperator *elem_restr =
-         fes->GetElementRestriction(ElementDofOrdering::LEXICOGRAPHIC);
-      MFEM_VERIFY(elem_restr != nullptr,
-                  "Element restriction is required for device bounds.");
-      elem_restr->Mult(*this, e_vec);
-
-      plb.GetElementBoundsKernel(rdim, fes_dim, e_vec, lower, upper, vdim);
-      lower.HostRead();
-      upper.HostRead();
-      return;
-   }
-
    int nel  = fes->GetNE();
    int fes_dim  = fes->GetVDim();
    lower.SetSize(nel*(vdim > 0 ? 1 :fes_dim));
@@ -5792,4 +5758,4 @@ std::pair<real_t, real_t> GridFunction::EstimateFunctionMaximum(
    return std::make_pair(global_max_lower, global_max_upper);
 }
 
-}
+}

fem/gridfunc.hpp

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

fem/gslib.cpp

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

fem/gslib.hpp

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

fem/gslib/findpts_local_2.cpp

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

fem/gslib/findpts_local_3.cpp

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

fem/integ/bilininteg_hdiv_ea.cpp

@@ -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;

fem/integ/bilininteg_mass_kernels.hpp

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

fem/intrules.hpp

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

fem/linearform.cpp

@@ -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);
             }
          }

fem/lininteg.hpp

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

fem/pgridfunc.cpp

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

fem/pgridfunc.hpp

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

general/arrays_by_name.hpp

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

general/device.cpp

@@ -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();

general/forall.hpp

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

general/optparser.cpp

@@ -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);

general/optparser.hpp

@@ -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)

general/text.hpp

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

linalg/hypre.cpp

@@ -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();

linalg/ode.cpp

@@ -284,7 +284,7 @@ void RK4Solver::Init(TimeDependentOperator &f_)
    int n = f->Width();
    y.SetSize(n, mem_type);
    k.SetSize(n, mem_type);
-   z.SetSize(n, mem_type);
+   z.SetSize(n, mem_type);   
 }
 
 void RK4Solver::Step(Vector &x, real_t &t, real_t &dt)
@@ -316,6 +316,67 @@ void RK4Solver::Step(Vector &x, real_t &t, real_t &dt)
    t += dt;
 }
 
+void RK4Solver::SetSolution(const Vector& x, real_t ts)
+{
+   hist.clear(); //we need only the last time step at time ts
+   StepRecord rec{ts,x};
+   hist.push_back(rec);
+}
+
+void RK4Solver::AdjointStep(Vector &lambda, real_t &t, real_t &dt)
+{
+   MFEM_VERIFY(!hist.empty(), "No forward history.");
+   MFEM_VERIFY(adj_mode == AdjointMode::Discrete, "Wrong adjoint mode.");
+
+   //Restore the solution history, i.e., perform a forward step.
+   auto rec = hist.back();
+   hist.pop_back();
+   const real_t tn = rec.t;
+   const real_t h  = dt;
+   const Vector &xn = rec.x;
+   // Recompute forward stage states and stage RHS:
+   f->SetTime(tn);
+   f->Mult(xn, k1); // k1
+   add(xn, h/2.0, k1, y2); // y2
+
+   f->SetTime(tn+h/2.0);
+   f->Mult(y2, k); // k2
+   add(xn, h/2.0, k, y3); // y3
+
+   f->Mult(y3, k); // k3
+   add(xn, h, k, y4); // y4
+
+
+   //set the initial values for a1,a2,a3,a4
+   a1.Set(h/6.0,lambda);
+   a2.Set(h/3.0,lambda);
+   a3.Set(h/3.0,lambda);
+   a4.Set(h/6.0,lambda);
+   
+   //stage 4
+   f->SetTime(tn+h);
+   f->JacobianMultTranspose(y4, a4, k);
+   lambda.Add(1.0,k);
+   a3.Add(h,k);
+
+   //stage 3
+   f->SetTime(tn+h/2.0);
+   f->JacobianMultTranspose(y3, a3, k);
+   lambda.Add(1.0,k);
+   a2.Add(h/2.0,k);
+
+   //stage 2
+   f->SetTime(tn+h/2.0);
+   f->JacobianMultTranspose(y2, a2, k);
+   lambda.Add(1.0,k);
+   a1.Add(h/2.0,k);
+
+   //stage 1
+   f->SetTime(tn);
+   f->JacobianMultTranspose(xn, a1, k);
+   lambda.Add(1.0,k);
+}
+
 ExplicitRKSolver::ExplicitRKSolver(int s_, const real_t *a_, const real_t *b_,
                                    const real_t *c_)
 {

linalg/ode.hpp

@@ -190,6 +190,78 @@ public:
       while (t < tf) { Step(x, t, dt); }
    }
 
+   /// @brief The enum defines the supported adjoint modes.
+   enum class AdjointMode { None, Continuous, Discrete };
+
+   AdjointMode adj_mode= AdjointMode::None;
+
+   /// Sets the solution at time ts. The solution will be used by the adjoint step
+   /// to recompute the internal stages/steps between ts and the time t used in the 
+   /// AdjointStep operator.
+   /// SetSolution(x,ts); and then calling AdjointStep(lambda, t, dt), will
+   /// first restore the solution/parametrization between  ts and t (if required by 
+   /// the time integration scheme) and then execute the adjoint step. 
+   virtual void SetSolution(const Vector& x, real_t ts){} 
+
+   
+   /// Return primal state x(t) for t in the recorded time interval.
+   /// Default false => solver does not provide primal interpolation access.
+   virtual bool GetForwardSolution(real_t t, Vector &x) const 
+   { 
+      return false; 
+   }
+
+   /// Perform one backward step for the adjoint variable lambda.
+   /// dt > 0 requests stepping from t to t - dt.
+   /// The solver may adjust dt (adaptive / grid snapping) and returns dt used.
+   /** @param[in,out] lambda   Approximate solution.
+       @param[in,out] t   Time associated with the approximate solution @a x.
+       @param[in,out] dt  Time step size.
+   */
+   virtual void AdjointStep(Vector &lambda, real_t &t, real_t &dt)
+   {
+      MFEM_ABORT("Adjoint stepping is not supported by this ODESolver.");
+   }
+
+   /// Perform adjoint integration from time @a t [in] to time @a tf [in].
+   /** @param[in,out] lambda   Approximate adjoint solution.
+       @param[in,out] t        Time associated with the approximate solution @a lambda.
+       @param[in,out] dt  Time step size.
+       @param[in]     t0  Requested final time. (t0<t)
+   */
+   virtual void AdjointRun(Vector &lambda, real_t &t, real_t &dt, real_t t0)
+   {
+      MFEM_VERIFY(t0 <= t, "AdjointRun expects t0 <= t.");
+      while (t > t0) { 
+         real_t d = std::min(dt, t - t0); 
+         AdjointStep(lambda, t, d); 
+         dt = d; 
+      }
+   }
+
+   /// Enable adjoint functionality and tell the solver 
+   /// what to record during forward stepping. 
+   /// Default: abort (adjoint mode is not supported).
+   virtual void EnableAdjoint(AdjointMode mod)
+   {
+      if(AdjointMode::None!=mod)
+      {
+         MFEM_ABORT("Adjoint is not supported by the ODESolver.");
+      }
+   }
+
+   /// Return true of the ODESolver supports the AdjointMode
+   /// and false if not.
+   virtual bool SupportsAdjoint(AdjointMode mod) const
+   { 
+      if(AdjointMode::None==mod)
+      {
+         return true;
+      }
+      //else
+      return false;
+   }
+
    /// Returns how many State vectors the ODE requires
    virtual int GetStateSize() { return 0; };
 
@@ -309,10 +381,73 @@ class RK4Solver : public ODESolver
 private:
    Vector y, k, z;
 
+   //storage for the adjoint solver
+   /// @brief stores the state x at time t
+   struct StepRecord{
+      real_t t;
+      Vector x;
+   };
+
+   /// @brief stores several steps
+   std::vector<StepRecord> hist;
+
+   Vector k1,y2,y3,y4; //RK4 stages
+   Vector a1,a2,a3,a4; //adjoint vectors
+
 public:
    void Init(TimeDependentOperator &f_) override;
 
    void Step(Vector &x, real_t &t, real_t &dt) override;
+
+   /// Return true of the ODESolver supports the AdjointMode and false
+   /// if not.
+   virtual bool SupportsAdjoint(AdjointMode mod) const override
+   { 
+      if(AdjointMode::Discrete==mod)
+      {
+         return true;
+      }
+      //else
+      return false;
+   }
+
+   /// Enable adjoint functionality and tell the solver 
+   /// what to record during forward stepping. 
+   /// Default: abort (adjoint mode is not supported).
+   virtual void EnableAdjoint(AdjointMode mod) override
+   {
+      if(AdjointMode::Continuous==mod)
+      {
+         MFEM_ABORT("Continuous adjoint not supported by the RK4Solver.");
+      }
+      adj_mode=mod;
+
+      // allocate the storage memeory
+      if(AdjointMode::Discrete==mod){
+         int n = f->Width();
+         k1.SetSize(n, mem_type);
+         y2.SetSize(n, mem_type);
+         y3.SetSize(n, mem_type);
+         y4.SetSize(n, mem_type);
+
+         a1.SetSize(n, mem_type);
+         a2.SetSize(n, mem_type);
+         a3.SetSize(n, mem_type);
+         a4.SetSize(n, mem_type);
+      }
+   }
+
+
+   void AdjointStep(Vector &lam, real_t &t, real_t &dt) override;
+
+   //Two posible ways to restore the histroy
+   //1) Use Step(Vector &x, real_t &t, real_t &dt) in discrete adjoint mode and record the internal states
+   //2) Use additional Method SetState(const Vector &x, real_t ts);
+   //   Calling AdjointStep(lam,t,dt) should restore the history between ts and t, and then execute the 
+   //   adjoint step in reverse to time t-dt with the recorded discrete states between ts and t. 
+   // Here we will use SetSolution
+   virtual void SetSolution(const Vector& x, real_t ts) override;
+
 };
 
 

linalg/operator.hpp

@@ -640,6 +640,15 @@ public:
    virtual int SUNMassMult(const Vector &x, Vector &v);
 
    virtual ~TimeDependentOperator() { }
+
+   /// y = (df/dx(x,t))^T * w
+   /// Default uses GetGradient(x).MultTranspose(w,y).
+   virtual void JacobianMultTranspose(const Vector &x,
+                                      const Vector &w,
+                                      Vector &y) const
+   {
+      GetGradient(x).MultTranspose(w, y);
+   }
 };
 
 

linalg/petsc.cpp

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

makefile

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

mesh/mesh.cpp

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

mesh/mesh.hpp

@@ -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

mesh/mesh_operators.cpp

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

mesh/mesh_operators.hpp

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

mesh/mesh_readers.cpp

@@ -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))
    {

mesh/pmesh.cpp

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

mesh/submesh/psubmesh.hpp

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

mesh/submesh/submesh.hpp

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

miniapps/CMakeLists.txt

@@ -31,6 +31,8 @@ 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)

miniapps/electromagnetics/CMakeLists.txt

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

miniapps/electromagnetics/makefile

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

miniapps/gslib/CMakeLists.txt

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

miniapps/gslib/makefile

@@ -83,16 +83,11 @@ include $(MFEM_TEST_MK)
 
 # Testing: Parallel vs. serial runs
 RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
-ifeq ($(MFEM_USE_MPI),YES)
-   RUN_MPI_SER = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) 1
-else
-   RUN_MPI_SER =
-endif
 TEST_NAME := GSLIB miniapp
 %-test-par: %
 	@$(call mfem-test,$<, $(RUN_MPI), $(TEST_NAME))
 %-test-seq: %
-	@$(call mfem-test,$<, $(RUN_MPI_SER), $(TEST_NAME))
+	@$(call mfem-test,$<,, $(TEST_NAME))
 
 # Testing: "test" target and mfem-test* variables are defined in config/test.mk
 

miniapps/hdiv-linear-solver/change_basis.cpp

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

miniapps/mtop/CMakeLists.txt

@@ -26,5 +26,10 @@ if(MFEM_USE_MPI AND MFEM_USE_DOUBLE)
 					 MAIN mtop_test_iso_elasticity.cpp
 					 ${PAR_MTOP_COMMON_FILES}
 					 LIBRARIES mfem)
+	add_mfem_miniapp(tst_rk4_adj
+				 	MAIN tst_rk4_adj.cpp
+					${PAR_MTOP_COMMON_FILES}
+                                         LIBRARIES mfem)
+
 
 endif (MFEM_USE_MPI AND MFEM_USE_DOUBLE)

miniapps/mtop/anisotropic/CMakeLists.txt

@@ -0,0 +1,32 @@
+# Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+if(MFEM_USE_MPI)
+
+list(APPEND AN_SOURCES linear_anisotropic_elasticity.hpp)
+list(APPEND AN_HEADERS linear_anisotropic_elasticity.cpp)
+
+convert_filenames_to_full_paths(AN_SOURCES)
+convert_filenames_to_full_paths(AN_HEADERS)
+
+set(EX_COMMON_FILES
+	EXTRA_SOURCES ${AN_SOURCES}
+	EXTRA_HEADERS ${AN_HEADERS})
+
+add_mfem_miniapp(mtop_test_anisotropic_le
+                      MAIN test_anisotropic_le.cpp
+		      ${EX_COMMON_FILES}
+                      LIBRARIES mfem)
+
+
+endif (MFEM_USE_MPI)
+
+

miniapps/mtop/chpt/CMakeLists.txt

@@ -0,0 +1,65 @@
+# 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.
+
+list(APPEND CHPT_SOURCES )
+list(APPEND CHPT_HEADERS )
+
+convert_filenames_to_full_paths(CHPT_SOURCES)
+convert_filenames_to_full_paths(CHPT_HEADERS)
+
+set(CHPT_COMMON_FILES
+	EXTRA_SOURCES ${CHPT_SOURCES}
+	EXTRA_HEADERS ${CHPT_HEADERS})
+
+add_mfem_miniapp(mtop_test_scalar_chpt
+                      MAIN tst_scalar_field.cpp
+		      ${CHPT_COMMON_FILES}
+                      LIBRARIES mfem)
+
+add_mfem_miniapp(mtop_test_vector_chpt
+                      MAIN tst_vector_field.cpp
+                      ${CHPT_COMMON_FILES}
+                      LIBRARIES mfem)
+	      
+add_mfem_miniapp(mtop_test_fixedslot
+		 MAIN tst_fixedslot_scalar_vector.cpp
+		 ${CHPT_COMMON_FILES}
+                 LIBRARIES mfem)
+
+add_mfem_miniapp(mtop_test_segmentio
+	         MAIN tst_segmentio_scalar_vector_mfem.cpp
+		 ${CHPT_COMMON_FILES}
+                 LIBRARIES mfem)
+
+add_mfem_miniapp(mtop_test_file
+		 MAIN tst_fileio_scalar_vector_mfem.cpp
+		 ${CHPT_COMMON_FILES}
+                 LIBRARIES mfem)
+
+
+add_mfem_miniapp(mtop_test_revolve
+                 MAIN revolve.cpp
+                 ${CHPT_COMMON_FILES}
+                 LIBRARIES mfem)
+
+add_mfem_miniapp(mtop_test_dyn
+                 MAIN tst_dynamic_chpt.cpp
+                 ${CHPT_COMMON_FILES}
+                 LIBRARIES mfem)
+
+add_mfem_miniapp(mtop_test_rev_fixedslot
+                 MAIN tst_revolve_chpt_fixedslot.cpp
+                 ${CHPT_COMMON_FILES}
+                 LIBRARIES mfem)
+
+
+
+

miniapps/mtop/chpt/dynamic_checkpointing.hpp

@@ -0,0 +1,309 @@
+#ifndef MFEM_DYNAMIC_CHECKPOINTING_HPP
+#define MFEM_DYNAMIC_CHECKPOINTING_HPP
+
+#include "mfem.hpp"
+
+#include <map>
+#include <memory>
+#include <limits>
+#include <vector>
+#include <utility>
+
+namespace mfem
+{
+
+template <typename Snapshot>
+class InMemoryCheckpointStorage
+{
+public:
+   using Handle = int;
+
+   Handle InvalidHandle() const { return -1; }
+   bool IsValid(const Handle &h) const { return h >= 0; }
+
+   Handle Store(Snapshot &&snap)
+   {
+      Handle id = InvalidHandle();
+      if (!free_.empty())
+      {
+         id = free_.back();
+         free_.pop_back();
+
+         if (slots_[id])
+         {
+            // reuse existing allocation
+            *slots_[id] = std::move(snap);
+         }
+         else
+         {
+            slots_[id].reset(new Snapshot(std::move(snap)));
+         }
+      }
+      else
+      {
+         id = (Handle) slots_.size();
+         slots_.push_back(std::unique_ptr<Snapshot>(new Snapshot(std::move(snap))));
+      }
+      return id;
+   }
+
+
+   template <typename Func>
+   void Read(const Handle &h, Func &&f) const
+   {
+      MFEM_VERIFY(IsValid(h), "InMemoryCheckpointStorage: invalid handle.");
+      MFEM_VERIFY(h < (Handle) slots_.size(), "InMemoryCheckpointStorage: handle out of range.");
+      MFEM_VERIFY(slots_[h].get() != nullptr, "InMemoryCheckpointStorage: empty slot.");
+      f(*slots_[h]);
+   }
+
+   void Erase(Handle &h)
+   {
+      if (!IsValid(h)) { h = InvalidHandle(); return; }
+      MFEM_VERIFY(h < (Handle) slots_.size(), "InMemoryCheckpointStorage: handle out of range.");
+      slots_[h].reset();
+      free_.push_back(h);
+      h = InvalidHandle();
+   }
+
+private:
+   std::vector<std::unique_ptr<Snapshot>> slots_;
+   std::vector<Handle> free_;
+};
+
+
+/**
+ * Dynamic checkpointing manager (Wang–Moin–Iaccarino 2009), with pluggable storage.
+ *
+ * Snapshot: checkpointed object (often the primal State itself).
+ * Storage : must provide:
+ *   using Handle = ...
+ *   Handle InvalidHandle() const;
+ *   bool   IsValid(const Handle&) const;
+ *   Handle Store(Snapshot&&);
+ *   template<class F> void Read(const Handle&, F&&) const; // calls f(const Snapshot&)
+ *   void Erase(Handle&);
+ */
+template <typename Snapshot,
+          typename Storage = InMemoryCheckpointStorage<Snapshot>>
+class DynamicCheckpointing
+{
+public:
+   using Step = long long;
+   using Handle = typename Storage::Handle;
+
+   struct Checkpoint
+   {
+      int level = 0;
+      Handle h; // InvalidHandle => placeholder
+      Checkpoint() = default;
+      Checkpoint(int lvl, const Handle &inv) : level(lvl), h(inv) {}
+   };
+
+   explicit DynamicCheckpointing(int s)
+      : s_(s),
+        owned_storage_(new Storage()),
+        storage_(owned_storage_.get())
+   {
+      MFEM_VERIFY(s_ > 0, "DynamicCheckpointing: s must be > 0.");
+      Reset();
+   }
+
+   DynamicCheckpointing(int s, Storage &external_storage)
+      : s_(s),
+        owned_storage_(nullptr),
+        storage_(&external_storage)
+   {
+      MFEM_VERIFY(s_ > 0, "DynamicCheckpointing: s must be > 0.");
+      Reset();
+   }
+
+   ~DynamicCheckpointing() { ReleaseAll(); }
+
+   void Reset()
+   {
+      ReleaseAll();
+      cps_.emplace(Step(0), Checkpoint(InfLevel(), storage_->InvalidHandle()));
+   }
+
+   struct CheckpointInfo
+   {
+      Step step;
+      int level;
+      bool stored;
+   };
+
+   std::vector<CheckpointInfo> GetCheckpointInfo() const
+   {
+      std::vector<CheckpointInfo> out;
+      out.reserve(cps_.size());
+      for (const auto &kv : cps_)
+      {
+         out.push_back({kv.first, kv.second.level, storage_->IsValid(kv.second.h)});
+      }
+      return out;
+   }
+
+   Step GetMaxStep() const
+   {
+      MFEM_VERIFY(!cps_.empty(), "DynamicCheckpointing: checkpoint map is empty.");
+      return cps_.rbegin()->first;
+   }
+
+   template <typename State, typename PrimalStepFn, typename MakeSnapshotFn>
+   void ForwardStep(const Step i,
+                    State &u_i_inout,
+                    PrimalStepFn &&primal_step,
+                    MakeSnapshotFn &&make_snapshot)
+   {
+      AllocateCheckpointForNextStep(i);
+
+      auto it = cps_.find(i);
+      if (it != cps_.end())
+      {
+         storage_->Erase(it->second.h);
+         Snapshot snap = make_snapshot(u_i_inout);
+         it->second.h = storage_->Store(std::move(snap));
+      }
+
+      primal_step(u_i_inout, i);
+   }
+
+   template <typename State,
+             typename AdjState,
+             typename PrimalStepFn,
+             typename AdjointStepFn,
+             typename MakeSnapshotFn,
+             typename RestoreSnapshotFn>
+   void BackwardStep(const Step i,
+                     AdjState &q_ip1_inout,
+                     State &u_work_inout,
+                     PrimalStepFn &&primal_step,
+                     AdjointStepFn &&adjoint_step,
+                     MakeSnapshotFn &&make_snapshot,
+                     RestoreSnapshotFn &&restore_snapshot)
+   {
+      // remove placeholder at i+1
+      const Step ph = i + 1;
+      auto it_ph = cps_.find(ph);
+
+      MFEM_VERIFY(it_ph != cps_.end(),
+                  "DynamicCheckpointing: expected checkpoint at i+1 before BackwardStep.");
+      storage_->Erase(it_ph->second.h);
+      cps_.erase(it_ph);
+
+      MFEM_ASSERT(GetMaxStep() <= i,
+                  "DynamicCheckpointing: found a checkpoint beyond current adjoint step.");
+
+      auto restore_from_handle = [&](Handle &h)
+      {
+         storage_->Read(h, [&](const Snapshot &snap)
+         {
+            restore_snapshot(snap, u_work_inout);
+         });
+         storage_->Erase(h); // retrieved => placeholder (Algorithm 4 semantics)
+      };
+
+      if (GetMaxStep() == i)
+      {
+         Handle h = TakeHandleMakePlaceholder(i);
+         restore_from_handle(h);
+      }
+      else
+      {
+         const Step k = GetMaxStep();
+         Handle hk = TakeHandleMakePlaceholder(k);
+         restore_from_handle(hk);
+
+         for (Step t = k; t < i; ++t)
+         {
+            ForwardStep(t, u_work_inout, primal_step, make_snapshot);
+         }
+      }
+
+      adjoint_step(q_ip1_inout, u_work_inout, i);
+   }
+
+private:
+   int s_ = 0;
+   std::unique_ptr<Storage> owned_storage_;
+   Storage *storage_ = nullptr;
+   std::map<Step, Checkpoint> cps_;
+
+   static int InfLevel() { return std::numeric_limits<int>::max(); }
+
+   void ReleaseAll()
+   {
+      if (!storage_) { cps_.clear(); return; }
+      for (auto &kv : cps_) { storage_->Erase(kv.second.h); }
+      cps_.clear();
+   }
+
+   bool FindDispensableLargestStep(Step &out_step) const
+   {
+      int max_level_seen = std::numeric_limits<int>::min();
+      for (auto it = cps_.rbegin(); it != cps_.rend(); ++it)
+      {
+         const Step step = it->first;
+         const int lvl = it->second.level;
+         if (max_level_seen > lvl) { out_step = step; return true; }
+         max_level_seen = (lvl > max_level_seen) ? lvl : max_level_seen;
+      }
+      return false;
+   }
+
+   void AllocateCheckpointForNextStep(const Step i)
+   {
+      const Step new_step = i + 1;
+      MFEM_VERIFY(cps_.find(new_step) == cps_.end(),
+                  "DynamicCheckpointing: checkpoint at i+1 already exists.");
+
+      const Handle inv = storage_->InvalidHandle();
+
+      // allow growth to s+1 entries (incl. placeholder)
+      if ((int)cps_.size() <= s_)
+      {
+         cps_.emplace(new_step, Checkpoint(0, inv));
+         return;
+      }
+
+      Step disp = -1;
+      if (FindDispensableLargestStep(disp))
+      {
+         auto it = cps_.find(disp);
+         MFEM_ASSERT(it != cps_.end(), "Internal error: dispensable checkpoint not found.");
+         storage_->Erase(it->second.h);
+         cps_.erase(it);
+         cps_.emplace(new_step, Checkpoint(0, inv));
+         return;
+      }
+
+      auto it_i = cps_.find(i);
+      MFEM_VERIFY(it_i != cps_.end(),
+                  "DynamicCheckpointing: promotion expected checkpoint at step i but none found.");
+      MFEM_VERIFY(i != 0, "DynamicCheckpointing: attempted to remove step 0 checkpoint.");
+
+      const int l = it_i->second.level;
+      storage_->Erase(it_i->second.h);
+      cps_.erase(it_i);
+      cps_.emplace(new_step, Checkpoint(l + 1, inv));
+   }
+
+   Handle TakeHandleMakePlaceholder(const Step i)
+   {
+      auto it = cps_.find(i);
+      MFEM_VERIFY(it != cps_.end(),
+                  "DynamicCheckpointing: TakeHandle requested a non-existent checkpoint.");
+      MFEM_VERIFY(storage_->IsValid(it->second.h),
+                  "DynamicCheckpointing: TakeHandle requested a checkpoint with no snapshot.");
+
+      Handle h = std::move(it->second.h);
+      it->second.h = storage_->InvalidHandle();
+      return h;
+   }
+};
+
+} // namespace mfem
+
+#endif // MFEM_DYNAMIC_CHECKPOINTING_HPP
+

miniapps/mtop/chpt/file_checkpoint_storage.hpp

@@ -0,0 +1,291 @@
+#ifndef MFEM_FILE_CHECKPOINT_STORAGE_HPP
+#define MFEM_FILE_CHECKPOINT_STORAGE_HPP
+
+#include "mfem.hpp"
+
+#include <cstdint>
+#include <cstdio>    // std::remove, std::rename
+#include <cstring>
+#include <string>
+#include <vector>
+#include <sstream>
+#include <iomanip>
+#include <fstream>
+#include <type_traits>
+
+#if __has_include(<filesystem>)
+  #include <filesystem>
+  namespace mfem_fs = std::filesystem;
+  #define MFEM_HAVE_FILESYSTEM 1
+#else
+  #define MFEM_HAVE_FILESYSTEM 0
+#endif
+
+namespace mfem
+{
+
+/**
+ * @brief Default binary serializer for checkpoint snapshots.
+ *
+ * Contract:
+ *   - static void     Write(std::ostream&, const Snapshot&)
+ *   - static Snapshot Read (std::istream&)
+ *
+ * Default implementation supports trivially-copyable POD types.
+ * Specialization for mfem::Vector is provided below.
+ *
+ * Portability note:
+ *   This binary format is NOT portable across endianness or differing sizeof(real_t).
+ *   If you need portability, provide a custom SnapshotIO (e.g., text, XDR, HDF5).
+ */
+template <typename Snapshot, typename Enable = void>
+struct DefaultCheckpointBinaryIO
+{
+   static void Write(std::ostream &, const Snapshot &)
+   {
+      static_assert(sizeof(Snapshot) == 0,
+                    "DefaultCheckpointBinaryIO: no implementation for this Snapshot type. "
+                    "Provide a custom SnapshotIO with Write/Read.");
+   }
+
+   static Snapshot Read(std::istream &)
+   {
+      static_assert(sizeof(Snapshot) == 0,
+                    "DefaultCheckpointBinaryIO: no implementation for this Snapshot type. "
+                    "Provide a custom SnapshotIO with Write/Read.");
+      return Snapshot();
+   }
+};
+
+// POD / trivially-copyable types (double, int, structs of POD, etc.)
+template <typename Snapshot>
+struct DefaultCheckpointBinaryIO<
+   Snapshot,
+   typename std::enable_if<std::is_trivially_copyable<Snapshot>::value>::type>
+{
+   static void Write(std::ostream &os, const Snapshot &x)
+   {
+      os.write(reinterpret_cast<const char*>(&x), sizeof(Snapshot));
+      MFEM_VERIFY(os.good(), "DefaultCheckpointBinaryIO: failed to write POD snapshot.");
+   }
+
+   static Snapshot Read(std::istream &is)
+   {
+      Snapshot x;
+      is.read(reinterpret_cast<char*>(&x), sizeof(Snapshot));
+      MFEM_VERIFY(is.good(), "DefaultCheckpointBinaryIO: failed to read POD snapshot.");
+      return x;
+   }
+};
+
+// Specialization for mfem::Vector
+template <>
+struct DefaultCheckpointBinaryIO<mfem::Vector, void>
+{
+   static void Write(std::ostream &os, const mfem::Vector &v)
+   {
+      const std::int64_t n = (std::int64_t) v.Size();
+      os.write(reinterpret_cast<const char*>(&n), sizeof(n));
+      MFEM_VERIFY(os.good(), "VectorBinaryIO: failed to write vector size.");
+
+      if (n > 0)
+      {
+         const mfem::real_t *data = v.GetData();
+         os.write(reinterpret_cast<const char*>(data),
+                  (std::streamsize)(n * (std::int64_t)sizeof(mfem::real_t)));
+         MFEM_VERIFY(os.good(), "VectorBinaryIO: failed to write vector data.");
+      }
+   }
+
+   static mfem::Vector Read(std::istream &is)
+   {
+      std::int64_t n = 0;
+      is.read(reinterpret_cast<char*>(&n), sizeof(n));
+      MFEM_VERIFY(is.good(), "VectorBinaryIO: failed to read vector size.");
+      MFEM_VERIFY(n >= 0, "VectorBinaryIO: invalid negative vector size.");
+
+      mfem::Vector v((int)n);
+      if (n > 0)
+      {
+         mfem::real_t *data = v.GetData();
+         is.read(reinterpret_cast<char*>(data),
+                 (std::streamsize)(n * (std::int64_t)sizeof(mfem::real_t)));
+         MFEM_VERIFY(is.good(), "VectorBinaryIO: failed to read vector data.");
+      }
+      return v;
+   }
+};
+
+
+/**
+ * @brief File-based checkpoint storage backend.
+ *
+ * Stores each snapshot in a separate file:
+ *   <directory>/<prefix><id><extension>
+ *
+ * Handle is a monotonically-increasing integer id (reused if keep_files==false).
+ *
+ * Template parameters:
+ *  - Snapshot   : stored snapshot type
+ *  - SnapshotIO : serializer with static Write/Read methods (see DefaultCheckpointBinaryIO)
+ *
+ * Threading:
+ *  - Not thread-safe. Typical adjoint/checkpointing usage is single-threaded control flow.
+ */
+template <typename Snapshot,
+          typename SnapshotIO = DefaultCheckpointBinaryIO<Snapshot>>
+class FileCheckpointStorage
+{
+public:
+   using Handle = std::int64_t;
+
+   /**
+    * @param directory   directory where checkpoint files live (created if create_dir==true)
+    * @param prefix      filename prefix (e.g. "ckpt_")
+    * @param extension   filename extension (e.g. ".bin")
+    * @param create_dir  create directory if missing (requires <filesystem>)
+    * @param keep_files  if true, Erase() will NOT delete files (useful for debugging),
+    *                    and ids are NOT reused.
+    */
+   FileCheckpointStorage(const std::string &directory,
+                         const std::string &prefix = "ckpt_",
+                         const std::string &extension = ".bin",
+                         bool create_dir = true,
+                         bool keep_files = false)
+      : dir_(directory),
+        prefix_(prefix),
+        ext_(extension),
+        keep_files_(keep_files)
+   {
+      MFEM_VERIFY(!dir_.empty(), "FileCheckpointStorage: directory must be non-empty.");
+      MFEM_VERIFY(!prefix_.empty(), "FileCheckpointStorage: prefix must be non-empty.");
+      MFEM_VERIFY(!ext_.empty(), "FileCheckpointStorage: extension must be non-empty.");
+
+      if (create_dir)
+      {
+#if MFEM_HAVE_FILESYSTEM
+         std::error_code ec;
+         mfem_fs::create_directories(mfem_fs::path(dir_), ec);
+         MFEM_VERIFY(!ec, "FileCheckpointStorage: failed to create directory: " << dir_);
+#else
+         MFEM_ABORT("FileCheckpointStorage: create_dir=true requires <filesystem> support.");
+#endif
+      }
+   }
+
+   Handle InvalidHandle() const { return (Handle)-1; }
+   bool IsValid(const Handle &h) const { return h >= 0; }
+
+   /**
+    * @brief Store snapshot to file, return handle.
+    *
+    * Uses atomic-ish pattern:
+    *  - write to "<path>.tmp"
+    *  - rename to "<path>"
+    */
+   Handle Store(Snapshot &&snap)
+   {
+      const Handle id = AllocateId_();
+      const std::string path = Path_(id);
+      const std::string tmp  = path + ".tmp";
+
+      {
+         std::ofstream os(tmp, std::ios::binary | std::ios::trunc);
+         MFEM_VERIFY(os.is_open(), "FileCheckpointStorage: failed to open for write: " << tmp);
+
+         // Write payload
+         SnapshotIO::Write(os, snap);
+         MFEM_VERIFY(os.good(), "FileCheckpointStorage: write failed for: " << tmp);
+      }
+
+      // Rename tmp -> final
+      const int rc = std::rename(tmp.c_str(), path.c_str());
+      MFEM_VERIFY(rc == 0, "FileCheckpointStorage: rename failed: " << tmp << " -> " << path);
+
+      return id;
+   }
+
+   /**
+    * @brief Read snapshot from file and pass it to callback f(const Snapshot&).
+    *
+    * The reference passed to f is valid only during the call.
+    */
+   template <typename Func>
+   void Read(const Handle &h, Func &&f) const
+   {
+      MFEM_VERIFY(IsValid(h), "FileCheckpointStorage: Read called with invalid handle.");
+      const std::string path = Path_(h);
+
+      std::ifstream is(path, std::ios::binary);
+      MFEM_VERIFY(is.is_open(), "FileCheckpointStorage: failed to open for read: " << path);
+
+      Snapshot snap = SnapshotIO::Read(is);
+      MFEM_VERIFY(is.good(), "FileCheckpointStorage: read failed for: " << path);
+
+      f(snap);
+   }
+
+   /**
+    * @brief Erase snapshot (delete file unless keep_files==true), set handle invalid.
+    *
+    * Id reuse policy:
+    *  - if keep_files_ == false: deleted ids are reused to avoid creating huge numbers of files
+    *  - if keep_files_ == true : ids are not reused (avoid overwriting old debug files)
+    */
+   void Erase(Handle &h)
+   {
+      if (!IsValid(h)) { h = InvalidHandle(); return; }
+
+      if (!keep_files_)
+      {
+         const std::string path = Path_(h);
+         const int rc = std::remove(path.c_str());
+         MFEM_VERIFY(rc == 0, "FileCheckpointStorage: failed to remove file: " << path);
+
+         free_.push_back(h);
+      }
+
+      h = InvalidHandle();
+   }
+
+private:
+   std::string dir_;
+   std::string prefix_;
+   std::string ext_;
+   bool keep_files_ = false;
+
+   mutable Handle next_id_ = 0;
+   mutable std::vector<Handle> free_;
+
+   Handle AllocateId_()
+   {
+      if (!keep_files_ && !free_.empty())
+      {
+         const Handle id = free_.back();
+         free_.pop_back();
+         return id;
+      }
+      return next_id_++;
+   }
+
+   std::string Path_(Handle id) const
+   {
+      MFEM_VERIFY(IsValid(id), "FileCheckpointStorage: Path_ called with invalid id.");
+
+      std::ostringstream oss;
+      oss << prefix_ << std::setw(12) << std::setfill('0') << id << ext_;
+
+#if MFEM_HAVE_FILESYSTEM
+      mfem_fs::path p = mfem_fs::path(dir_) / mfem_fs::path(oss.str());
+      return p.string();
+#else
+      // Fallback: simple concatenation; assumes dir_ ends without trailing slash if needed.
+      return dir_ + "/" + oss.str();
+#endif
+   }
+};
+
+} // namespace mfem
+
+#endif // MFEM_FILE_CHECKPOINT_STORAGE_HPP
+

miniapps/mtop/chpt/fixed_slot_checkpoint_storage.hpp

@@ -0,0 +1,385 @@
+#ifndef MFEM_FIXED_SLOT_CHECKPOINT_STORAGE_HPP
+#define MFEM_FIXED_SLOT_CHECKPOINT_STORAGE_HPP
+
+#include "mfem.hpp"
+
+#include <cstddef>
+#include <cstdint>
+#include <cstring>
+#include <string>
+#include <vector>
+#include <fstream>
+#include <limits>
+#include <type_traits>
+
+namespace mfem
+{
+
+// ------------------------------------------------------------
+// Packer 1: trivially-copyable snapshots (memcpy)
+// ------------------------------------------------------------
+template <typename Snapshot>
+class TrivialFixedPacker
+{
+public:
+   TrivialFixedPacker()
+   {
+      static_assert(std::is_trivially_copyable<Snapshot>::value,
+                    "TrivialFixedPacker requires Snapshot to be trivially copyable.");
+   }
+
+   std::size_t SlotBytes() const { return sizeof(Snapshot); }
+
+   void Pack(const Snapshot &snap, void *dst) const
+   {
+      std::memcpy(dst, &snap, sizeof(Snapshot));
+   }
+
+   void Unpack(const void *src, Snapshot &out) const
+   {
+      std::memcpy(&out, src, sizeof(Snapshot));
+   }
+};
+
+
+// ------------------------------------------------------------
+// Packer 2: mfem::Vector of fixed length n
+// (Serialized size is fixed: n*sizeof(real_t))
+// ------------------------------------------------------------
+class FixedVectorPacker
+{
+public:
+   explicit FixedVectorPacker(int n) : n_(n)
+   {
+      MFEM_VERIFY(n_ > 0, "FixedVectorPacker: n must be > 0.");
+   }
+
+   std::size_t SlotBytes() const
+   {
+      return (std::size_t)n_ * sizeof(mfem::real_t);
+   }
+
+   void Pack(const mfem::Vector &v, void *dst) const
+   {
+      MFEM_VERIFY(v.Size() == n_, "FixedVectorPacker: vector size mismatch.");
+      std::memcpy(dst, v.GetData(), SlotBytes());
+   }
+
+   void Unpack(const void *src, mfem::Vector &out) const
+   {
+      out.SetSize(n_);
+      std::memcpy(out.GetData(), src, SlotBytes());
+   }
+
+   int Size() const { return n_; }
+
+private:
+   int n_ = 0;
+};
+
+
+// ------------------------------------------------------------
+// Fixed-slot MEMORY storage: one big RAM block
+// ------------------------------------------------------------
+template <typename Snapshot, typename Packer = TrivialFixedPacker<Snapshot>>
+class FixedSlotMemoryCheckpointStorage
+{
+public:
+   using Handle = int;
+
+   FixedSlotMemoryCheckpointStorage(int max_slots, const Packer &packer = Packer())
+      : max_slots_(max_slots), packer_(packer)
+   {
+      MFEM_VERIFY(max_slots_ > 0, "FixedSlotMemoryCheckpointStorage: max_slots must be > 0.");
+
+      slot_bytes_ = packer_.SlotBytes();
+      MFEM_VERIFY(slot_bytes_ > 0, "FixedSlotMemoryCheckpointStorage: SlotBytes must be > 0.");
+
+      // Single contiguous block
+      bytes_.resize((std::size_t)max_slots_ * slot_bytes_);
+
+      in_use_.assign(max_slots_, 0);
+      free_.reserve(max_slots_);
+      for (int i = 0; i < max_slots_; ++i) { free_.push_back(i); }
+   }
+
+   Handle InvalidHandle() const { return -1; }
+   bool IsValid(const Handle &h) const { return h >= 0; }
+
+   int MaxSlots() const { return max_slots_; }
+   std::size_t SlotBytes() const { return slot_bytes_; }
+
+   Handle Store(Snapshot &&snap)
+   {
+      MFEM_VERIFY(!free_.empty(),
+                  "FixedSlotMemoryCheckpointStorage: out of slots (increase max_slots).");
+
+      const int slot = free_.back();
+      free_.pop_back();
+      in_use_[slot] = 1;
+
+      void *dst = SlotPtr_(slot);
+      packer_.Pack(snap, dst);
+
+      return slot;
+   }
+
+   template <typename Func>
+   void Read(const Handle &h, Func &&f) const
+   {
+      MFEM_VERIFY(IsValid(h), "FixedSlotMemoryCheckpointStorage: Read invalid handle.");
+      MFEM_VERIFY(h < max_slots_, "FixedSlotMemoryCheckpointStorage: Read handle out of range.");
+      MFEM_VERIFY(in_use_[h] == 1, "FixedSlotMemoryCheckpointStorage: Read from free slot.");
+
+      Snapshot tmp;
+      const void *src = SlotPtrConst_(h);
+      packer_.Unpack(src, tmp);
+
+      f(tmp);
+   }
+
+   void Erase(Handle &h)
+   {
+      if (!IsValid(h)) { h = InvalidHandle(); return; }
+
+      MFEM_VERIFY(h < max_slots_, "FixedSlotMemoryCheckpointStorage: Erase handle out of range.");
+      MFEM_VERIFY(in_use_[h] == 1, "FixedSlotMemoryCheckpointStorage: double-free / invalid erase.");
+
+      in_use_[h] = 0;
+      free_.push_back(h);
+      h = InvalidHandle();
+   }
+
+   /// Optional: return all slots to the free list (does not zero memory).
+   void Reset()
+   {
+      free_.clear();
+      for (int i = 0; i < max_slots_; ++i) { in_use_[i] = 0; free_.push_back(i); }
+   }
+
+private:
+   int max_slots_ = 0;
+   std::size_t slot_bytes_ = 0;
+   Packer packer_;
+
+   std::vector<unsigned char> bytes_; // single block
+   std::vector<unsigned char> in_use_;
+   std::vector<int> free_;
+
+   void *SlotPtr_(int slot)
+   {
+      return (void*)(&bytes_[(std::size_t)slot * slot_bytes_]);
+   }
+
+   const void *SlotPtrConst_(int slot) const
+   {
+      return (const void*)(&bytes_[(std::size_t)slot * slot_bytes_]);
+   }
+};
+
+
+// ------------------------------------------------------------
+// Fixed-slot FILE storage: one single pre-sized file
+// ------------------------------------------------------------
+template <typename Snapshot, typename Packer = TrivialFixedPacker<Snapshot>>
+class FixedSlotFileCheckpointStorage
+{
+public:
+   using Handle = int;
+
+   struct Header
+   {
+      char magic[8];              // "MFCKPTFS"
+      std::uint64_t version;      // 1
+      std::uint64_t slot_bytes;
+      std::uint64_t max_slots;
+      std::uint64_t reserved[4];  // future use / padding
+   };
+
+   FixedSlotFileCheckpointStorage(const std::string &path,
+                                 int max_slots,
+                                 const Packer &packer = Packer(),
+                                 bool truncate = true,
+                                 bool flush_on_store = false)
+      : path_(path),
+        max_slots_(max_slots),
+        packer_(packer),
+        flush_on_store_(flush_on_store)
+   {
+      MFEM_VERIFY(!path_.empty(), "FixedSlotFileCheckpointStorage: empty file path.");
+      MFEM_VERIFY(max_slots_ > 0, "FixedSlotFileCheckpointStorage: max_slots must be > 0.");
+
+      slot_bytes_ = packer_.SlotBytes();
+      MFEM_VERIFY(slot_bytes_ > 0, "FixedSlotFileCheckpointStorage: SlotBytes must be > 0.");
+
+      Open_(truncate);
+
+      in_use_.assign(max_slots_, 0);
+      free_.reserve(max_slots_);
+      for (int i = 0; i < max_slots_; ++i) { free_.push_back(i); }
+
+      scratch_.resize(slot_bytes_);
+   }
+
+   ~FixedSlotFileCheckpointStorage()
+   {
+      if (file_.is_open()) { file_.close(); }
+   }
+
+   Handle InvalidHandle() const { return -1; }
+   bool IsValid(const Handle &h) const { return h >= 0; }
+
+   int MaxSlots() const { return max_slots_; }
+   std::size_t SlotBytes() const { return slot_bytes_; }
+   const std::string &Path() const { return path_; }
+
+   Handle Store(Snapshot &&snap)
+   {
+      MFEM_VERIFY(!free_.empty(),
+                  "FixedSlotFileCheckpointStorage: out of slots (increase max_slots).");
+
+      const int slot = free_.back();
+      free_.pop_back();
+      in_use_[slot] = 1;
+
+      // Pack into scratch buffer then write into fixed slot offset
+      packer_.Pack(snap, scratch_.data());
+
+      const std::uint64_t off = SlotOffset_(slot);
+      file_.seekp((std::streamoff)off, std::ios::beg);
+      MFEM_VERIFY(file_.good(), "FixedSlotFileCheckpointStorage: seekp failed.");
+
+      file_.write(reinterpret_cast<const char*>(scratch_.data()),
+                  (std::streamsize)slot_bytes_);
+      MFEM_VERIFY(file_.good(), "FixedSlotFileCheckpointStorage: write failed.");
+
+      if (flush_on_store_) { file_.flush(); }
+
+      return slot;
+   }
+
+   template <typename Func>
+   void Read(const Handle &h, Func &&f) const
+   {
+      MFEM_VERIFY(IsValid(h), "FixedSlotFileCheckpointStorage: Read invalid handle.");
+      MFEM_VERIFY(h < max_slots_, "FixedSlotFileCheckpointStorage: Read handle out of range.");
+      MFEM_VERIFY(in_use_[h] == 1, "FixedSlotFileCheckpointStorage: Read from free slot.");
+
+      const std::uint64_t off = SlotOffset_(h);
+      file_.seekg((std::streamoff)off, std::ios::beg);
+      MFEM_VERIFY(file_.good(), "FixedSlotFileCheckpointStorage: seekg failed.");
+
+      file_.read(reinterpret_cast<char*>(scratch_.data()),
+                 (std::streamsize)slot_bytes_);
+      MFEM_VERIFY(file_.good(), "FixedSlotFileCheckpointStorage: read failed.");
+
+      Snapshot tmp;
+      packer_.Unpack(scratch_.data(), tmp);
+      f(tmp);
+   }
+
+   void Erase(Handle &h)
+   {
+      if (!IsValid(h)) { h = InvalidHandle(); return; }
+
+      MFEM_VERIFY(h < max_slots_, "FixedSlotFileCheckpointStorage: Erase handle out of range.");
+      MFEM_VERIFY(in_use_[h] == 1, "FixedSlotFileCheckpointStorage: double-free / invalid erase.");
+
+      // No file deletion; just return slot to free list.
+      in_use_[h] = 0;
+      free_.push_back(h);
+      h = InvalidHandle();
+   }
+
+   /// Optional: return all slots to free list (file contents remain).
+   void Reset()
+   {
+      free_.clear();
+      for (int i = 0; i < max_slots_; ++i) { in_use_[i] = 0; free_.push_back(i); }
+   }
+
+private:
+   std::string path_;
+   int max_slots_ = 0;
+   std::size_t slot_bytes_ = 0;
+   Packer packer_;
+   bool flush_on_store_ = false;
+
+   // mutable because Read() is const but needs to seek/read
+   mutable std::fstream file_;
+   mutable std::vector<unsigned char> scratch_;
+
+   std::vector<unsigned char> in_use_;
+   std::vector<int> free_;
+
+   static Header MakeHeader_(std::uint64_t slot_bytes, std::uint64_t max_slots)
+   {
+      Header h;
+      std::memset(&h, 0, sizeof(h));
+      h.magic[0] = 'M'; h.magic[1] = 'F'; h.magic[2] = 'C'; h.magic[3] = 'K';
+      h.magic[4] = 'P'; h.magic[5] = 'T'; h.magic[6] = 'F'; h.magic[7] = 'S';
+      h.version = 1;
+      h.slot_bytes = slot_bytes;
+      h.max_slots  = max_slots;
+      return h;
+   }
+
+   void Open_(bool truncate)
+   {
+      const std::ios::openmode mode =
+         std::ios::binary | std::ios::in | std::ios::out | (truncate ? std::ios::trunc : (std::ios::openmode)0);
+
+      file_.open(path_.c_str(), mode);
+      MFEM_VERIFY(file_.is_open(), "FixedSlotFileCheckpointStorage: failed to open file.");
+
+      const Header expected = MakeHeader_((std::uint64_t)slot_bytes_, (std::uint64_t)max_slots_);
+
+      if (truncate)
+      {
+         // Write header
+         file_.seekp(0, std::ios::beg);
+         file_.write(reinterpret_cast<const char*>(&expected), sizeof(expected));
+         MFEM_VERIFY(file_.good(), "FixedSlotFileCheckpointStorage: header write failed.");
+
+         // Pre-size file to: header + max_slots*slot_bytes
+         const std::uint64_t total = (std::uint64_t)sizeof(Header)
+                                   + (std::uint64_t)max_slots_ * (std::uint64_t)slot_bytes_;
+
+         MFEM_VERIFY(total > 0, "FixedSlotFileCheckpointStorage: invalid total file size.");
+         file_.seekp((std::streamoff)(total - 1), std::ios::beg);
+         MFEM_VERIFY(file_.good(), "FixedSlotFileCheckpointStorage: seekp for resize failed.");
+
+         const char zero = 0;
+         file_.write(&zero, 1);
+         MFEM_VERIFY(file_.good(), "FixedSlotFileCheckpointStorage: resize write failed.");
+         file_.flush();
+      }
+      else
+      {
+         // Validate existing header
+         Header got;
+         file_.seekg(0, std::ios::beg);
+         file_.read(reinterpret_cast<char*>(&got), sizeof(got));
+         MFEM_VERIFY(file_.good(), "FixedSlotFileCheckpointStorage: header read failed.");
+
+         MFEM_VERIFY(std::memcmp(got.magic, expected.magic, 8) == 0,
+                     "FixedSlotFileCheckpointStorage: magic mismatch.");
+         MFEM_VERIFY(got.version == expected.version,
+                     "FixedSlotFileCheckpointStorage: version mismatch.");
+         MFEM_VERIFY(got.slot_bytes == expected.slot_bytes,
+                     "FixedSlotFileCheckpointStorage: slot_bytes mismatch.");
+         MFEM_VERIFY(got.max_slots  == expected.max_slots,
+                     "FixedSlotFileCheckpointStorage: max_slots mismatch.");
+      }
+   }
+
+   std::uint64_t SlotOffset_(int slot) const
+   {
+      return (std::uint64_t)sizeof(Header) + (std::uint64_t)slot * (std::uint64_t)slot_bytes_;
+   }
+};
+
+} // namespace mfem
+
+#endif // MFEM_FIXED_SLOT_CHECKPOINT_STORAGE_HPP
+

miniapps/mtop/chpt/revolve.cpp

@@ -0,0 +1,694 @@
+#include "mfem.hpp"
+
+#include <cstdint>
+#include <cstring>
+#include <fstream>
+#include <string>
+#include <vector>
+#include <algorithm>
+
+using mfem::out;
+
+// -------------------------------
+// Fixed-slot storage: Memory
+// -------------------------------
+class FixedSlotMemoryStorage
+{
+private:
+   int max_slots_ = 0;
+   size_t slot_bytes_ = 0;
+   std::vector<uint8_t> data_;
+
+public:
+   FixedSlotMemoryStorage(int max_slots, size_t slot_bytes)
+      : max_slots_(max_slots), slot_bytes_(slot_bytes),
+        data_(size_t(max_slots)*slot_bytes, 0)
+   {
+      MFEM_VERIFY(max_slots_ > 0, "max_slots must be > 0");
+      MFEM_VERIFY(slot_bytes_ > 0, "slot_bytes must be > 0");
+   }
+
+   int MaxSlots() const { return max_slots_; }
+   size_t SlotBytes() const { return slot_bytes_; }
+
+   void Save(int slot, const uint8_t *src, size_t bytes)
+   {
+      MFEM_VERIFY(0 <= slot && slot < max_slots_, "Save: slot out of range");
+      MFEM_VERIFY(bytes == slot_bytes_, "Save: bytes mismatch");
+      MFEM_VERIFY(src != nullptr, "Save: src is null");
+      std::memcpy(data_.data() + size_t(slot)*slot_bytes_, src, slot_bytes_);
+   }
+
+   void Load(int slot, uint8_t *dst, size_t bytes) const
+   {
+      MFEM_VERIFY(0 <= slot && slot < max_slots_, "Load: slot out of range");
+      MFEM_VERIFY(bytes == slot_bytes_, "Load: bytes mismatch");
+      MFEM_VERIFY(dst != nullptr, "Load: dst is null");
+      std::memcpy(dst, data_.data() + size_t(slot)*slot_bytes_, slot_bytes_);
+   }
+};
+
+// -------------------------------
+// Fixed-slot storage: Single file
+// -------------------------------
+class FixedSlotFileStorage
+{
+private:
+   int max_slots_ = 0;
+   size_t slot_bytes_ = 0;
+   std::string filename_;
+   mutable std::fstream file_;
+
+   std::streamoff Offset(int slot) const
+   {
+      return std::streamoff(size_t(slot) * slot_bytes_);
+   }
+
+public:
+   FixedSlotFileStorage(const std::string &filename,
+                        int max_slots,
+                        size_t slot_bytes)
+      : max_slots_(max_slots), slot_bytes_(slot_bytes), filename_(filename)
+   {
+      MFEM_VERIFY(max_slots_ > 0, "max_slots must be > 0");
+      MFEM_VERIFY(slot_bytes_ > 0, "slot_bytes must be > 0");
+      MFEM_VERIFY(!filename_.empty(), "filename must not be empty");
+
+      // Create/truncate file and size it.
+      {
+         std::ofstream ofs(filename_, std::ios::binary | std::ios::trunc);
+         MFEM_VERIFY(ofs.good(), "Failed to create checkpoint file");
+         const size_t total = size_t(max_slots_) * slot_bytes_;
+         if (total > 0)
+         {
+            ofs.seekp(std::streamoff(total - 1));
+            char zero = 0;
+            ofs.write(&zero, 1);
+         }
+      }
+
+      file_.open(filename_, std::ios::binary | std::ios::in | std::ios::out);
+      MFEM_VERIFY(file_.good(), "Failed to open checkpoint file");
+   }
+
+   ~FixedSlotFileStorage()
+   {
+      if (file_.is_open()) { file_.close(); }
+   }
+
+   int MaxSlots() const { return max_slots_; }
+   size_t SlotBytes() const { return slot_bytes_; }
+
+   void Save(int slot, const uint8_t *src, size_t bytes)
+   {
+      MFEM_VERIFY(0 <= slot && slot < max_slots_, "Save: slot out of range");
+      MFEM_VERIFY(bytes == slot_bytes_, "Save: bytes mismatch");
+      MFEM_VERIFY(src != nullptr, "Save: src is null");
+
+      file_.seekp(Offset(slot));
+      MFEM_VERIFY(file_.good(), "Save: seekp failed");
+      file_.write(reinterpret_cast<const char*>(src), std::streamsize(slot_bytes_));
+      MFEM_VERIFY(file_.good(), "Save: write failed");
+      file_.flush();
+   }
+
+   void Load(int slot, uint8_t *dst, size_t bytes) const
+   {
+      MFEM_VERIFY(0 <= slot && slot < max_slots_, "Load: slot out of range");
+      MFEM_VERIFY(bytes == slot_bytes_, "Load: bytes mismatch");
+      MFEM_VERIFY(dst != nullptr, "Load: dst is null");
+
+      file_.seekg(Offset(slot));
+      MFEM_VERIFY(file_.good(), "Load: seekg failed");
+      file_.read(reinterpret_cast<char*>(dst), std::streamsize(slot_bytes_));
+      MFEM_VERIFY(file_.good(), "Load: read failed");
+   }
+};
+
+// -------------------------------
+// REVOLVE controller (actions)
+// (Transcribed from revolve.c used in ADOL-C; Algorithm 799.)
+// -------------------------------
+enum class RevolveAction
+{
+   advance,
+   takeshot,
+   restore,
+   firsturn,
+   youturn,
+   terminate
+};
+
+class RevolveController
+{
+private:
+   int snaps_ = 0;
+   int check_ = -1;
+   int capo_  = 0;
+   int fine_  = 0;
+
+   int turn_ = 0;
+   int oldfine_ = 0;
+   int oldsnaps_ = 0;
+
+   // Stack of checkpoint times, indexed by 'check_'.
+   std::vector<int> ch_;
+
+public:
+   RevolveController() = default;
+
+   RevolveController(int snaps, int capo0, int fine0)
+      : snaps_(snaps), check_(-1), capo_(capo0), fine_(fine0),
+        turn_(0), oldfine_(fine0), oldsnaps_(snaps), ch_(snaps, 0)
+   {
+      MFEM_VERIFY(snaps_ > 0, "REVOLVE snaps must be > 0");
+      MFEM_VERIFY(capo_ <= fine_, "REVOLVE: capo must be <= fine");
+
+      // Match revolve.c initialization behavior.
+      if (check_ == -1 && capo_ < fine_)
+      {
+         turn_ = 0;
+         ch_[0] = capo_ - 1;
+      }
+   }
+
+   int Snaps() const { return snaps_; }
+   int Check() const { return check_; }
+   int Capo() const  { return capo_; }
+   int Fine() const  { return fine_; }
+   const std::vector<int>& CheckpointTimes() const { return ch_; }
+
+   RevolveAction Next()
+   {
+      MFEM_VERIFY(!(check_ < -1), "REVOLVE: check < -1");
+      MFEM_VERIFY(!(capo_ > fine_), "REVOLVE: capo > fine");
+
+      if ((check_ == -1) && (capo_ < fine_))
+      {
+         turn_ = 0;
+         ch_[0] = capo_ - 1;
+      }
+
+      const int diff = fine_ - capo_;
+      switch (diff)
+      {
+         case 0:
+         {
+            // Terminate or restore to next checkpoint on the stack.
+            if (check_ == -1 || capo_ == ch_[0])
+            {
+               check_ -= 1; // mirror revolve.c behavior
+               return RevolveAction::terminate;
+            }
+            else
+            {
+               capo_ = ch_[check_];
+               oldfine_ = fine_;
+               return RevolveAction::restore;
+            }
+         }
+         case 1:
+         {
+            // One adjoint step available.
+            fine_ -= 1;
+            if (check_ >= 0 && ch_[check_] == capo_) { check_ -= 1; }
+
+            if (turn_ == 0)
+            {
+               turn_ = 1;
+               oldfine_ = fine_;
+               return RevolveAction::firsturn;
+            }
+            else
+            {
+               oldfine_ = fine_;
+               return RevolveAction::youturn;
+            }
+         }
+         default:
+         {
+            // diff > 1
+            if (check_ == -1 || ch_[check_] != capo_)
+            {
+               // Take a new checkpoint at current capo.
+               check_ += 1;
+               MFEM_VERIFY(check_ + 1 <= snaps_, "REVOLVE: exceeded snaps");
+               ch_[check_] = capo_;
+               oldfine_ = fine_;
+               return RevolveAction::takeshot;
+            }
+            else
+            {
+               // Advance capo forward within (capo, fine).
+               // This follows the binomial logic in revolve.c.
+               MFEM_VERIFY(!((oldfine_ < fine_) && (snaps_ == check_ + 1)),
+                           "REVOLVE: fine increased unexpectedly with full stack");
+
+               const int oldcapo = capo_;
+               const int ds = snaps_ - check_;
+               MFEM_VERIFY(ds >= 1, "REVOLVE: ds < 1");
+
+               int reps = 0;
+               long long range = 1;
+               while (range < (fine_ - capo_))
+               {
+                  reps += 1;
+                  // range = range*(reps+ds)/reps  (integer arithmetic)
+                  range = range * (reps + ds) / reps;
+               }
+               MFEM_VERIFY(reps >= 1, "REVOLVE: reps < 1");
+
+               // Binomial helper values (integer)
+               const long long bino1 = range * reps / (ds + reps);
+               const long long bino2 = (ds > 1) ? (bino1 * ds / (ds + reps - 1)) : 1;
+               const long long bino3 =
+                  (ds == 1) ? 0 :
+                  (ds > 2)  ? (bino2 * (ds - 1) / (ds + reps - 2)) : 1;
+               const long long bino4 = bino2 * (reps - 1) / ds;
+               const long long bino5 =
+                  (ds < 3) ? 0 :
+                  (ds > 3) ? (bino3 * (ds - 2) / reps) : 1;
+
+               // Kowarz "new version": keep l^ as small as possible
+               const long long bino6 = bino1 * ds / reps;
+
+               const long long gap = fine_ - capo_;
+               if (gap <= bino1 + bino3)
+               {
+                  capo_ += int(bino4);
+               }
+               else if (gap < bino1 + bino2)
+               {
+                  capo_ = fine_ - int(bino2 + bino3);
+               }
+               else if (gap <= bino1 + bino2 + bino5)
+               {
+                  capo_ += int(bino1 - bino3);
+               }
+               else
+               {
+                  capo_ = fine_ - int(bino6);
+               }
+
+               if (capo_ == oldcapo) { capo_ = oldcapo + 1; }
+
+               oldfine_ = fine_;
+               return RevolveAction::advance;
+            }
+         }
+      }
+   }
+};
+
+// -------------------------------
+// Fixed-step REVOLVE checkpointing
+// -------------------------------
+template <typename StorageT>
+class FixedStepRevolveCheckpointing
+{
+public:
+   struct Shot
+   {
+      int time = 0;   // state index
+      int slot = 0;   // checkpoint slot index
+   };
+
+private:
+   int num_steps_ = 0;        // total number of primal steps (0..num_steps)
+   int num_checkpoints_ = 0;  // REVOLVE "snaps"
+   size_t snapshot_bytes_ = 0;
+
+   StorageT *storage_ = nullptr;
+
+   std::vector<Shot> forward_shots_;
+   int forward_shot_cursor_ = 0;
+
+   // Controller state at the beginning of the reverse sweep (pre-firsturn).
+   RevolveController ctrl_init_;
+   RevolveController ctrl_;
+
+   // Two scratch buffers for (de)serialization.
+   std::vector<uint8_t> io_buf_;
+   std::vector<uint8_t> prefinal_buf_;
+   bool prefinal_valid_ = false;
+
+   // Reverse sweep bookkeeping.
+   bool reverse_started_ = false;
+   int u_work_time_ = -1;
+
+public:
+   FixedStepRevolveCheckpointing(int num_steps,
+                                int num_checkpoints,
+                                size_t snapshot_bytes,
+                                StorageT &storage)
+      : num_steps_(num_steps),
+        num_checkpoints_(num_checkpoints),
+        snapshot_bytes_(snapshot_bytes),
+        storage_(&storage),
+        io_buf_(snapshot_bytes, 0),
+        prefinal_buf_(snapshot_bytes, 0)
+   {
+      MFEM_VERIFY(num_steps_ >= 0, "num_steps must be >= 0");
+      MFEM_VERIFY(num_checkpoints_ > 0, "num_checkpoints must be > 0");
+      MFEM_VERIFY(snapshot_bytes_ > 0, "snapshot_bytes must be > 0");
+
+      MFEM_VERIFY(storage_->MaxSlots() == num_checkpoints_,
+                  "Storage MaxSlots() must match num_checkpoints");
+      MFEM_VERIFY(storage_->SlotBytes() == snapshot_bytes_,
+                  "Storage SlotBytes() must match snapshot_bytes");
+
+      BuildForwardPlanAndInitialControllerState();
+      Reset();
+   }
+
+   void Reset()
+   {
+      forward_shot_cursor_ = 0;
+      prefinal_valid_ = false;
+      reverse_started_ = false;
+      u_work_time_ = -1;
+      ctrl_ = ctrl_init_;
+   }
+
+   int NumSteps() const { return num_steps_; }
+   int NumCheckpoints() const { return num_checkpoints_; }
+   size_t SnapshotBytes() const { return snapshot_bytes_; }
+
+   // ForwardStep: called for i=0..num_steps-1
+   template <typename State, typename PrimalStep, typename MakeSnapshot>
+   void ForwardStep(int i,
+                    State &u,
+                    PrimalStep &&primal_step,
+                    MakeSnapshot &&make_snapshot)
+   {
+      MFEM_VERIFY(0 <= i && i < num_steps_, "ForwardStep: i out of range");
+
+      // Take any planned shots at time i (before advancing).
+      while (forward_shot_cursor_ < (int)forward_shots_.size() &&
+             forward_shots_[forward_shot_cursor_].time == i)
+      {
+         const int slot = forward_shots_[forward_shot_cursor_].slot;
+         make_snapshot(u, io_buf_.data(), snapshot_bytes_);
+         storage_->Save(slot, io_buf_.data(), snapshot_bytes_);
+         forward_shot_cursor_++;
+      }
+
+      // Cache u_{num_steps-1} so reverse can start there (REVOLVE expects capo=num_steps-1).
+      if (i == num_steps_ - 1)
+      {
+         make_snapshot(u, prefinal_buf_.data(), snapshot_bytes_);
+         prefinal_valid_ = true;
+      }
+
+      // Advance one step.
+      primal_step(i, u);
+   }
+
+   // BackwardStep: called for i=num_steps-1..0
+   template <typename State, typename AdjointState,
+             typename PrimalStep, typename AdjointStep,
+             typename MakeSnapshot, typename RestoreSnapshot>
+   void BackwardStep(int i,
+                     AdjointState &lambda,
+                     State &u_work,
+                     PrimalStep &&primal_step,
+                     AdjointStep &&adjoint_step,
+                     MakeSnapshot &&make_snapshot,
+                     RestoreSnapshot &&restore_snapshot)
+   {
+      MFEM_VERIFY(0 <= i && i < num_steps_, "BackwardStep: i out of range");
+
+      if (!reverse_started_)
+      {
+         MFEM_VERIFY(prefinal_valid_ || num_steps_ == 0,
+                     "Reverse started but prefinal state was not captured. "
+                     "Did you run the forward loop through i=num_steps-1?");
+
+         ctrl_ = ctrl_init_;
+         reverse_started_ = true;
+
+         if (num_steps_ > 0)
+         {
+            // Restore u_{num_steps-1} into u_work.
+            restore_snapshot(u_work, prefinal_buf_.data(), snapshot_bytes_);
+            u_work_time_ = num_steps_ - 1;
+         }
+      }
+
+      // Execute controller actions until we perform exactly one adjoint step.
+      for (;;)
+      {
+         RevolveAction act = ctrl_.Next();
+
+         switch (act)
+         {
+            case RevolveAction::takeshot:
+            {
+               const int slot = ctrl_.Check();
+               MFEM_VERIFY(u_work_time_ == ctrl_.Capo(),
+                           "takeshot: u_work_time must equal capo");
+               make_snapshot(u_work, io_buf_.data(), snapshot_bytes_);
+               storage_->Save(slot, io_buf_.data(), snapshot_bytes_);
+               break;
+            }
+            case RevolveAction::restore:
+            {
+               const int slot = ctrl_.Check();
+               storage_->Load(slot, io_buf_.data(), snapshot_bytes_);
+               restore_snapshot(u_work, io_buf_.data(), snapshot_bytes_);
+               u_work_time_ = ctrl_.Capo();
+               break;
+            }
+            case RevolveAction::advance:
+            {
+               const int target = ctrl_.Capo();
+               MFEM_VERIFY(u_work_time_ >= 0, "advance: u_work_time not initialized");
+               MFEM_VERIFY(target >= u_work_time_, "advance: target < current time");
+
+               for (int t = u_work_time_; t < target; ++t)
+               {
+                  primal_step(t, u_work);
+               }
+               u_work_time_ = target;
+               break;
+            }
+            case RevolveAction::firsturn:
+            case RevolveAction::youturn:
+            {
+               // After firsturn/youturn, ctrl_.Fine() has been decremented and equals ctrl_.Capo().
+               const int step = ctrl_.Fine();
+               MFEM_VERIFY(step == ctrl_.Capo(), "youturn: fine != capo");
+               MFEM_VERIFY(step == u_work_time_, "youturn: u_work_time != step");
+               MFEM_VERIFY(step == i, "BackwardStep called with i that doesn't match REVOLVE schedule");
+
+               // One adjoint step.
+               adjoint_step(step, u_work, lambda);
+               return;
+            }
+            case RevolveAction::terminate:
+            {
+               MFEM_ABORT("REVOLVE terminated early: BackwardStep called after completion?");
+               break;
+            }
+         }
+      }
+   }
+
+private:
+   void BuildForwardPlanAndInitialControllerState()
+   {
+      forward_shots_.clear();
+
+      RevolveController sim(num_checkpoints_, /*capo=*/0, /*fine=*/num_steps_);
+
+      // Forward-plan phase: only TAKESHOT/ADVANCE should occur while fine-capo>1.
+      while (sim.Fine() - sim.Capo() > 1)
+      {
+         RevolveAction a = sim.Next();
+         if (a == RevolveAction::takeshot)
+         {
+            forward_shots_.push_back({sim.Capo(), sim.Check()});
+         }
+         else if (a == RevolveAction::advance)
+         {
+            // nothing to record; capo moved forward inside sim
+         }
+         else
+         {
+            MFEM_ABORT("Unexpected REVOLVE action during forward planning phase");
+         }
+      }
+
+      // At this point, sim is in the pre-firsturn state (capo = num_steps-1, fine = num_steps)
+      // for num_steps>=1. For num_steps<=1, fine-capo<=1 from the start.
+      ctrl_init_ = sim;
+   }
+};
+
+// ============================================================
+// Mini example 1: scalar state (double)
+// u_{n+1} = factor * u_n
+// J = 0.5 (u_N - target)^2  => lambda_N = (u_N - target)
+// lambda_n = factor * lambda_{n+1}
+// ============================================================
+static void RunScalarExample(int Nsteps, int Ncheck, double factor)
+{
+   out << "\n--- Scalar REVOLVE example ---\n";
+   out << "Nsteps=" << Nsteps << ", Ncheckpoints=" << Ncheck << ", factor=" << factor << "\n";
+
+   const size_t snap_bytes = sizeof(double);
+   FixedSlotMemoryStorage storage(Ncheck, snap_bytes);
+   FixedStepRevolveCheckpointing<FixedSlotMemoryStorage> ckpt(Nsteps, Ncheck, snap_bytes, storage);
+
+   auto make_snapshot = [](const double &u, uint8_t *outb, size_t bytes)
+   {
+      MFEM_VERIFY(bytes == sizeof(double), "scalar snapshot size mismatch");
+      std::memcpy(outb, &u, sizeof(double));
+   };
+   auto restore_snapshot = [](double &u, const uint8_t *inb, size_t bytes)
+   {
+      MFEM_VERIFY(bytes == sizeof(double), "scalar snapshot size mismatch");
+      std::memcpy(&u, inb, sizeof(double));
+   };
+
+   auto primal_step = [factor](int /*step*/, double &u)
+   {
+      u *= factor;
+   };
+   auto adjoint_step = [factor](int /*step*/, const double & /*u_step*/, double &lambda)
+   {
+      lambda *= factor;
+   };
+
+   const double target = 2.0;
+
+   // Forward
+   double u = 1.0;
+   for (int i = 0; i < Nsteps; ++i)
+   {
+      ckpt.ForwardStep(i, u, primal_step, make_snapshot);
+   }
+   const double uN = u;
+
+   // Init adjoint at final state
+   double lambda = (uN - target);
+
+   // Reverse
+   double u_work = 0.0; // will be overwritten by restore_snapshot
+   for (int i = Nsteps - 1; i >= 0; --i)
+   {
+      ckpt.BackwardStep(i, lambda, u_work,
+                        primal_step, adjoint_step,
+                        make_snapshot, restore_snapshot);
+   }
+
+   // Analytic lambda_0 = factor^Nsteps * (uN - target)
+   double factorN = 1.0;
+   for (int k = 0; k < Nsteps; ++k) { factorN *= factor; }
+   const double lambda0_exact = factorN * (uN - target);
+
+   out << "uN = " << uN << "\n";
+   out << "lambda0 (computed) = " << lambda << "\n";
+   out << "lambda0 (exact)    = " << lambda0_exact << "\n";
+}
+
+// ============================================================
+// Mini example 2: mfem::Vector state
+// u_{n+1} = factor * u_n (elementwise scalar multiply)
+// J = 0.5 ||u_N - target||^2  => lambda_N = (u_N - target)
+// lambda_n = factor * lambda_{n+1}
+// ============================================================
+static void RunVectorExample(int Nsteps, int Ncheck, int dim, double factor)
+{
+   out << "\n--- mfem::Vector REVOLVE example ---\n";
+   out << "Nsteps=" << Nsteps << ", Ncheckpoints=" << Ncheck
+       << ", dim=" << dim << ", factor=" << factor << "\n";
+
+   MFEM_VERIFY(dim > 0, "dim must be > 0");
+
+   const size_t snap_bytes = sizeof(double) * size_t(dim);
+   FixedSlotMemoryStorage storage(Ncheck, snap_bytes);
+   FixedStepRevolveCheckpointing<FixedSlotMemoryStorage> ckpt(Nsteps, Ncheck, snap_bytes, storage);
+
+   auto make_snapshot = [](const mfem::Vector &u, uint8_t *outb, size_t bytes)
+   {
+      MFEM_VERIFY(bytes == sizeof(double) * size_t(u.Size()), "Vector snapshot size mismatch");
+      std::memcpy(outb, u.GetData(), bytes);
+   };
+   auto restore_snapshot = [](mfem::Vector &u, const uint8_t *inb, size_t bytes)
+   {
+      MFEM_VERIFY(bytes == sizeof(double) * size_t(u.Size()), "Vector snapshot size mismatch");
+      std::memcpy(u.GetData(), inb, bytes);
+   };
+
+   auto primal_step = [factor](int /*step*/, mfem::Vector &u)
+   {
+      u *= factor;
+   };
+   auto adjoint_step = [factor](int /*step*/, const mfem::Vector & /*u_step*/, mfem::Vector &lambda)
+   {
+      lambda *= factor;
+   };
+
+   mfem::Vector u(dim);
+   for (int k = 0; k < dim; ++k) { u[k] = 1.0 + 0.1*k; }
+
+   mfem::Vector target(dim);
+   target = 2.0;
+
+   // Forward
+   for (int i = 0; i < Nsteps; ++i)
+   {
+      ckpt.ForwardStep(i, u, primal_step, make_snapshot);
+   }
+   mfem::Vector uN(u); // copy for reporting
+
+   // Init adjoint: lambda_N = uN - target
+   mfem::Vector lambda(dim);
+   lambda = u;
+   lambda -= target;
+
+   // Reverse
+   mfem::Vector u_work(dim);
+   u_work = 0.0;
+   for (int i = Nsteps - 1; i >= 0; --i)
+   {
+      ckpt.BackwardStep(i, lambda, u_work,
+                        primal_step, adjoint_step,
+                        make_snapshot, restore_snapshot);
+   }
+
+   out << "||uN||_2      = " << uN.Norml2() << "\n";
+   out << "||lambda0||_2 = " << lambda.Norml2() << "\n";
+}
+
+int main(int argc, char *argv[])
+{
+   mfem::Device device("cpu");
+   device.Print();
+
+   int Nsteps = 20;
+   int Ncheck = 3;
+   int dim = 5;
+   double factor = 1.05;
+
+   mfem::OptionsParser args(argc, argv);
+   args.AddOption(&Nsteps, "-n", "--num-steps", "Number of primal steps.");
+   args.AddOption(&Ncheck, "-s", "--num-checkpoints", "Number of checkpoints (snaps).");
+   args.AddOption(&dim, "-d", "--dim", "Vector dimension for the mfem::Vector example.");
+   args.AddOption(&factor, "-f", "--factor", "Scalar factor in the toy update u_{n+1}=f*u_n.");
+   args.Parse();
+
+   if (!args.Good())
+   {
+      args.PrintUsage(out);
+      return 1;
+   }
+   args.PrintOptions(out);
+
+   RunScalarExample(Nsteps, Ncheck, factor);
+   RunVectorExample(Nsteps, Ncheck, dim, factor);
+
+   return 0;
+}
+