small updates prior to passing the code along for comparison

updating contact and parproblems files
cleanup on finite difference test code. Making it easier for others to use
2024-07-29 16:34:06 -07:00 · 2024-07-15 10:47:54 -07:00 · 2024-07-15 08:08:00 -07:00 · 2024-07-14 20:43:48 -07:00 · 2024-06-17 13:15:59 -07:00 · 2024-06-01 14:03:22 -07:00
197 changed files with 67091 additions and 1594 deletions
@@ -33,6 +33,7 @@ env:
  HYPRE_ARCHIVE: v2.19.0.tar.gz
  HYPRE_TOP_DIR: hypre-2.19.0
  METIS_ARCHIVE: metis-4.0.3.tar.gz
+  METIS_ARCHIVE_MAC: metis-4.0.3-mac.tgz
  METIS_TOP_DIR: metis-4.0.3
  MFEM_TOP_DIR: mfem

@@ -52,6 +53,7 @@ jobs:
        mpi: [seq, par]
        build-system: [make, cmake]
        hypre-target: [int32]
+        precision: [fp64]
        exclude:
          - os: ubuntu-latest
            build-system: cmake
@@ -75,6 +77,8 @@ jobs:
          - os: ubuntu-latest
            target: dbg
            config-opts: 'CPPFLAGS+=-Og'
+          - os: macos-latest
+            codecov: NO
          - os: windows-latest
            codecov: NO
          - os: windows-latest
@@ -87,6 +91,7 @@ jobs:
            mpi: par
            build-system: cmake
            hypre-target: int32
+            precision: fp64
            # This option can be set to pass additional configuration options to
            # the MFEM configuration command.
            # config-opts: '-DCMAKE_VERBOSE_MAKEFILE=ON'
@@ -96,7 +101,15 @@ jobs:
            mpi: par
            build-system: make
            hypre-target: int64
-    name: ${{ matrix.os }}-${{ matrix.build-system }}-${{ matrix.target }}-${{ matrix.mpi }}-${{ matrix.hypre-target }}
+            precision: fp64
+          - os: ubuntu-latest
+            target: opt
+            codecov: NO
+            mpi: par
+            build-system: make
+            hypre-target: int32
+            precision: fp32
+    name: ${{ matrix.os }}-${{ matrix.build-system }}-${{ matrix.target }}-${{ matrix.mpi }}-${{ matrix.hypre-target }}-${{ matrix.precision }}

    runs-on: ${{ matrix.os }}

@@ -126,6 +139,17 @@ jobs:
          # Fetch the complete history for codecov to access commits ID
          fetch-depth: 0

+      - name: Xcode version setup (MacOS)
+        if: matrix.os == 'macos-latest'
+        run: |
+          XCODE_PATH="/Applications/Xcode_15.3.app"
+          echo "> sudo xcode-select -s ${XCODE_PATH}"
+          sudo xcode-select -s ${XCODE_PATH}
+          echo "> g++ -v"
+          g++ -v
+          echo "> clang++ -v"
+          clang++ -v
+
      # Only get MPI if defined for the job.
      # TODO: It would be nice to have only one step, e.g. with a dedicated
      # action, but I (@adrienbernede) don't see how at the moment.
@@ -169,25 +193,27 @@ jobs:
        uses: actions/cache@v4
        with:
          path: ${{ env.HYPRE_TOP_DIR }}
-          key: ${{ runner.os }}-build-${{ env.HYPRE_TOP_DIR }}-${{ matrix.hypre-target }}-v2.2
+          key: ${{ runner.os }}-build-${{ env.HYPRE_TOP_DIR }}-${{ matrix.hypre-target }}-${{ matrix.precision }}-v2.5

      - name: get hypre
        if: matrix.mpi == 'par' && steps.hypre-cache.outputs.cache-hit != 'true' && matrix.os != 'windows-latest'
-        uses: mfem/github-actions/build-hypre@v2.4
+        uses: mfem/github-actions/build-hypre@v2.5
        with:
          archive: ${{ env.HYPRE_ARCHIVE }}
          dir: ${{ env.HYPRE_TOP_DIR }}
          target: ${{ matrix.hypre-target }}
          build-system: make
+          precision: ${{ matrix.precision }}

      - name: get hypre (Windows)
        if: matrix.mpi == 'par' && steps.hypre-cache.outputs.cache-hit != 'true' && matrix.os == 'windows-latest'
-        uses: mfem/github-actions/build-hypre@v2.4
+        uses: mfem/github-actions/build-hypre@v2.5
        with:
          archive: ${{ env.HYPRE_ARCHIVE }}
          dir: ${{ env.HYPRE_TOP_DIR }}
          target: ${{ matrix.hypre-target }}
          build-system: cmake
+          precision: ${{ matrix.precision }}

      # Get Metis through cache, or build it.
      # Install will only run on cache miss.
@@ -197,13 +223,13 @@ jobs:
        uses: actions/cache@v4
        with:
          path: ${{ env.METIS_TOP_DIR }}
-          key: ${{ runner.os }}-build-${{ env.METIS_TOP_DIR }}-v2.2
+          key: ${{ runner.os }}-build-${{ env.METIS_TOP_DIR }}-v2.5

      - name: install metis
        if: matrix.mpi == 'par' && matrix.os != 'windows-latest' && steps.metis-cache.outputs.cache-hit != 'true'
-        uses: mfem/github-actions/build-metis@v2.4
+        uses: mfem/github-actions/build-metis@v2.5
        with:
-          archive: ${{ env.METIS_ARCHIVE }}
+          archive: ${{ matrix.os != 'macos-latest' && env.METIS_ARCHIVE || env.METIS_ARCHIVE_MAC }}
          dir: ${{ env.METIS_TOP_DIR }}

      - name: cache vcpkg (Windows)
@@ -228,7 +254,7 @@ jobs:

      # MFEM build and test
      - name: build
-        uses: mfem/github-actions/build-mfem@v2.4
+        uses: mfem/github-actions/build-mfem@v2.5
        env:
          VCPKG_DEFAULT_BINARY_CACHE: ${{ github.workspace }}/vcpkg_cache
        with:
@@ -240,6 +266,7 @@ jobs:
          hypre-dir: ${{ env.HYPRE_TOP_DIR }}
          metis-dir: ${{ env.METIS_TOP_DIR }}
          mfem-dir: ${{ env.MFEM_TOP_DIR }}
+          precision: ${{ matrix.precision }}
          config-options: ${{ matrix.config-opts }}
          library-only: ${{ matrix.target == 'dbg' && matrix.os != 'ubuntu-latest' }}

@@ -282,7 +309,7 @@ jobs:
      # Code coverage (process and upload reports)
      - name: codecov
        if: matrix.codecov == 'YES'
-        uses: mfem/github-actions/upload-coverage@v2.4
+        uses: mfem/github-actions/upload-coverage@v2.5
        with:
          name: ${{ matrix.os }}-${{ matrix.build-system }}-${{ matrix.target }}-${{ matrix.mpi }}-${{ matrix.hypre-target }}
          project_dir: ${{ env.MFEM_TOP_DIR }}
@@ -53,11 +53,11 @@ jobs:
        uses: actions/cache@v4
        with:
          path: ${{ env.HYPRE_TOP_DIR }}
-          key: ${{ runner.os }}-build-${{ env.HYPRE_TOP_DIR }}-v2.2
+          key: ${{ runner.os }}-build-${{ env.HYPRE_TOP_DIR }}-v2.5

      - name: Get Hypre
        if: steps.hypre-cache.outputs.cache-hit != 'true'
-        uses: mfem/github-actions/build-hypre@v2.4
+        uses: mfem/github-actions/build-hypre@v2.5
        with:
          archive: ${{ env.HYPRE_ARCHIVE }}
          dir: ${{ env.HYPRE_TOP_DIR }}
@@ -68,18 +68,18 @@ jobs:
        uses: actions/cache@v4
        with:
          path: ${{ env.METIS_TOP_DIR }}
-          key: ${{ runner.os }}-build-${{ env.METIS_TOP_DIR }}-v2.2
+          key: ${{ runner.os }}-build-${{ env.METIS_TOP_DIR }}-v2.5

      - name: Install Metis
        if: steps.metis-cache.outputs.cache-hit != 'true'
-        uses: mfem/github-actions/build-metis@v2.4
+        uses: mfem/github-actions/build-metis@v2.5
        with:
          archive: ${{ env.METIS_ARCHIVE }}
          dir: ${{ env.METIS_TOP_DIR }}

      # MFEM build and test
      - name: build-mfem
-        uses: mfem/github-actions/build-mfem@v2.4
+        uses: mfem/github-actions/build-mfem@v2.5
        with:
          os: ${{ runner.os }}
          target: opt
@@ -44,7 +44,7 @@ jobs:
          path: mfem

      - name: MFEM Build
-        uses: mfem/github-actions/build-mfem@v2.4
+        uses: mfem/github-actions/build-mfem@v2.5
        with:
          os: ${{ runner.os }}
          target: opt
@@ -57,6 +57,8 @@ examples/ex2[0-9]
 examples/ex2[0-9]p
 examples/ex3[0-9]
 examples/ex3[0-9]p
+examples/ex4[0-9]
+examples/ex4[0-9]p

 examples/refined.mesh
 examples/displaced.mesh
@@ -232,7 +234,7 @@ miniapps/meshing/mobius-strip.mesh
 miniapps/meshing/klein-bottle.mesh
 miniapps/meshing/toroid-*.mesh
 miniapps/meshing/twist-*.mesh
-miniapps/meshing/mesh-explorer.mesh
+miniapps/meshing/mesh-explorer.mesh*
 miniapps/meshing/partitioning.txt
 miniapps/meshing/mesh-explorer-visit*
 miniapps/meshing/mesh-explorer-paraview/
@@ -369,6 +371,8 @@ miniapps/dpg/ParaView
 miniapps/spde/generate_random_field
 miniapps/spde/ParaView

+miniapps/tribol/contact-patch-test
+
 # Unit test binary and outputs
 tests/unit/output_meshes
 tests/unit/unit_tests
@@ -13,6 +13,9 @@
 # at Lawrence Livermore National Laboratory (LLNL). This entire pipeline is
 # LLNL-specific!

+include:
+  - project: 'lc-templates/id_tokens'
+    file: 'id_tokens.yml'

 # The pipeline is divided into stages. Usually, jobs in a given stage wait for
 # the preceding stages to complete before to start. However, we sometimes use
@@ -9,6 +9,10 @@
 # terms of the BSD-3 license. We welcome feedback and contributions, see file
 # CONTRIBUTING.md for details.

+include:
+  - project: 'lc-templates/id_tokens'
+    file: 'id_tokens.yml'
+
 # We define the following GitLab pipeline variables:
 variables:

@@ -35,9 +35,8 @@ variables:
    - when: on_success

 # Lassen uses a different job scheduler (spectrum lsf) that does not allow
-# pre-allocation the same way slurm does.  We use pdebug queue on lassen
-# to speed-up the allocation.  However this would not be scalable to
-# multiple builds.
+# pre-allocation the same way slurm does.  We use the pci queue on lassen
+# to speed-up the allocation.
 .build_and_test_on_lassen:
  extends: [.on_lassen]
  stage: build_and_test
@@ -45,5 +44,5 @@ variables:
    - echo ${MFEM_DATA_DIR}
    - echo ${SPEC}
    # Next script uses 'THREADS': leaving it empty --> it uses 'make all -j'
-    - lalloc 1 -W 45 -q pdebug --atsdisable tests/gitlab/build_and_test --spec "${SPEC}" --data-dir "${MFEM_DATA_DIR}" --data
+    - lalloc 1 -W 45 -q pci --atsdisable tests/gitlab/build_and_test --spec "${SPEC}" --data-dir "${MFEM_DATA_DIR}" --data
  needs: [setup]
@@ -52,4 +52,4 @@ variables:
    - echo ${JOBID}
    - echo ${MFEM_DATA_DIR}
    - echo ${SPEC}
-    - srun $( [[ -n "${JOBID}" ]] && echo "--jobid=${JOBID}" ) -t 45 -N 1 tests/gitlab/build_and_test --spec "${SPEC}" --data-dir "${MFEM_DATA_DIR}" --data
+    - srun $( [[ -n "${JOBID}" ]] && echo "--jobid=${JOBID}" ) --reservation=ci -t 45 -N 1 tests/gitlab/build_and_test --spec "${SPEC}" --data-dir "${MFEM_DATA_DIR}" --data
@@ -14,14 +14,14 @@ stages:
  - build_and_test
  - report

-opt_mpi_cuda_xl_16_1_1_12:
+opt_mpi_cuda_gcc:
  variables:
-    SPEC: "%xl@16.1.1.12 +mpi +cuda cuda_arch=70"
+    SPEC: "%gcc@8.3.1 +mpi +cuda cuda_arch=70"
  extends: .build_and_test_on_lassen

-opt_mpi_cuda_hypre_cuda_xl:
+opt_mpi_cuda_hypre_cuda_gcc:
  variables:
-    SPEC: "%xl@16.1.1.12 +mpi +cuda cuda_arch=70 ^hypre+cuda~shared cuda_arch=70"
+    SPEC: "%gcc@8.3.1 +mpi +cuda cuda_arch=70 ^hypre+cuda~shared cuda_arch=70"
  extends: .build_and_test_on_lassen

 # Jobs report
@@ -32,11 +32,11 @@ mkdir _${BASELINE_TEST} && cd _${BASELINE_TEST}

 # run
 if [[ "${MACHINE_NAME}" == "quartz" || "${MACHINE_NAME}" == "ruby" ]]; then
-  salloc --nodes=1 -p pdebug ../runtest ../../mfem "${BASELINE_TEST} ${TPLS_DIR}"
+  salloc --nodes=1 --reservation=ci ../runtest ../../mfem "${BASELINE_TEST} ${TPLS_DIR}"
 elif [[ ${MACHINE_NAME} == "corona" ]]; then
  salloc --nodes=1 -t 60 -p pbatch ../runtest ../../mfem "${BASELINE_TEST} ${TPLS_DIR}"
 elif [[ ${MACHINE_NAME} == "lassen" ]]; then
-  lalloc 1 -q pdebug ../runtest ../../mfem "${BASELINE_TEST} ${TPLS_DIR}"
+  lalloc 1 -q pci ../runtest ../../mfem "${BASELINE_TEST} ${TPLS_DIR}"
 else
  echo "Unknown machine: MACHINE_NAME=$MACHINE_NAME"
  exit 1
@@ -8,74 +8,110 @@
                               https://mfem.org


-Version 4.6.1 (development)
+Version 4.7.1 (development)
 ===========================

+- Added an MFEM example for the eikonal equation. This new solver is based on
+  the proximal Galerkin method introduced by Keith and Surowiec.
+
+
+Version 4.7, released on May 7, 2024
+====================================
+
+- Added support for single precision (with corresponding hypre build). The MFEM
+  floating point type was generalized from `double` to `real_t`. For details see
+  https://github.com/orgs/mfem/discussions/4207.
+
 Meshing improvements
 --------------------
+- Added the capability to partition (big) serial meshes in serial code, see the
+  new classes MeshPartitioner and MeshPart. This capability is also exposed as a
+  menu option in the mesh-explorer miniapp in miniapps/meshing.
+
+- Added named attribute sets and basic supporting methods to the Mesh class as a
+  convenient means of referring to sets of domain or boundary attribute numbers.
+  See the new Example 39/39p and data/compass.mesh.
+
 - Introduced formulas for refinement of patches in NURBS meshes. Refinement by
  arbitrary integer factors is also enabled, e.g. in the mesh-explorer miniapp.
  NURBS coarsening and knot removal are also introduced.

-Discretization improvements
---------------------------
- Introduced support for higher order non conformal Nedelec elements on
-  simplices in ParMesh.
- Introduced support for internal boundary elements in nonconformal adapted
-  meshes.
-
- Added functionality for construction of cut-surface and cut-volume
-  IntegrationRules through a moment-fitting approach. The cut is specified by
-  the zero level set of a Coefficient. See fem/intrules_cut.hpp and Example 38.
-
- Added a new nonlinear integrator, `HyperbolicFormIntegrator`. This implements 
-  both element-wise weak divergence and face-wise numerical flux for a general
-  system of hyperbolic conservation laws. To use this integrator for a specific
-  flux function, users can define a derived class of `FluxFunction`. Currently,
-  advection, Burgers', shallow-water, Euler equations (see, Example 18) are
-  available.
-
-GPU support
----------------------------
- Added support for full assembly on simplices.
- Added functionality for BilinearFormIntegrators to use kernels that work for both
-  tensor and unstructured elements.
- Added partial assembly for linear elasticity. Does not use sum factorization for now.
-
-New and updated examples and miniapps
-------------------------------------
- Added a new block solver in miniapp/solvers for the Darcy problem.
-  The new solver is based on a Bramble-Pasciak preconditioning. User can
-  use and implement their own preconditioner for the mass matrix.
-
- Added miniapp to demonstrate new elasticity integrator and unstructured element GPU support,
-  and a block diagonal preconditioner using low order refinement. Allows comparison with
-  currently existing legacy mode integrator. See miniapps/solvers/lor_elast.
-
-Miscellaneous
-------------
- Added support for single and double precision, with corresponding hypre build.
-  Generalized the floating point type from `double` to `real_t`. For more
-  details see https://github.com/orgs/mfem/discussions/4207.
+- Added support for internal boundary elements in nonconforming meshes.

 - The ReadCubit Genesis mesh importer has been rewritten to improve readability.

+Discretization improvements
+---------------------------
+- Added a new nonlinear integrator, `HyperbolicFormIntegrator` that implements
+  both element-wise weak divergence and face-wise numerical flux for a general
+  system of hyperbolic conservation laws. To use the integrator for a specific
+  flux function, users can define a derived class of `FluxFunction`. Currently,
+  advection, Burgers, shallow-water and Euler equations (see Example 18/18p) are
+  available.
+
+- Added a capability to construct cut-surface and cut-volume IntegrationRules
+  through a moment-fitting approach. The cut is specified by the zero level set
+  of a Coefficient. See fem/intrules_cut.hpp and the new Example 38.
+
+- Introduced support for high-order nonconforming Nedelec elements on simplices.
+
+GPU computing
+-------------
+- Added partial assembly and GPU support for the DG diffusion integrator.
+
+- Efficient GPU-accelerated LOR assembly is now supported on surface meshes.
+
+- Added functionality to automatically configure hypre's compute policy to match
+  MFEM's compute policy when hypre is built with GPU support. Requires version
+  hypre-2.31.0 or later.
+
+- Added support for full assembly on simplices.
+
+- Added partial assembly for linear elasticity (no sum factorization for now).
+
+- Added functionality for BilinearFormIntegrators to use kernels that work for
+  both tensor and unstructured elements.
+
+- The RAJA backend will use `seq_exec` for serial loop execution when RAJA
+  v2023.06.00 and beyond is detected as `loop_exec` is deprecated.
+
+- API change: The macro MFEM_HYPRE_FORALL (from hypre.hpp) which was intended
+  for internal use, has been removed and replaced by the function template
+  mfem::hypre_forall in general/forall.hpp.
+
+New and updated examples and miniapps
+-------------------------------------
+- Added a new miniapp illustrating elastic contact based on the Tribol library,
+  (https://github.com/LLNL/Tribol). See miniapps/tribol.
+
+- Added a miniapp to demonstrate low order refined (LOR) block preconditioning
+  for linear elasticity on GPUs. See miniapps/solvers/lor_elast.
+
+- Added a new block solver in miniapp/solvers for the Darcy problem. The new
+  solver is based on a Bramble-Pasciak preconditioning. User can use and
+  implement their own preconditioner for the mass matrix.
+
+- Added a small miniapp for printing the shape functions of a KnotVector. See
+  miniapps/nurbs/nurbs_printfunc.cpp.
+
+- Added two new example codes: 38 and 39/39p described above. Substantially
+  updated Example 18/18p.
+
+Miscellaneous
+-------------
 - Updated the Doxygen documentation style, which now requires Doxygen version
  1.9.8 or later. See the doc/ directory.

- Improved thread safety for global variables in the library, for example
-  IntegrationRules IntRules, RefinedIntRules, GeometryRefiner
-  GlobGeometryRefiner, and FiniteElement::dof2quad_array.
+- Improved thread safety for global variables in the library, e.g. for IntRules,
+  RefinedIntRules, GlobGeometryRefiner, and FiniteElement::dof2quad_array.

 - PETSc integration now generally requires PETSc version 3.21 or later, though
  depending on the functionality older versions may still work.

- RAJA backend will use seq_exec for serial loop execution when RAJA
-  v2023.06.00 and beyond is detected as loop_exec is deprecated.
+- Various other simplifications, extensions, and bugfixes in the code.
+
+- Added GSLIB-based gather-scatter operator.

- Adding named attribute sets and basic supporting methods to the Mesh class as
-  a convenient means of referring to sets of domain or boundary attribute
-  numbers. Also adding related serial and parallel examples which illustrate.

 Version 4.6, released on September 27, 2023
 ===========================================
@@ -96,7 +132,6 @@ Meshing improvements
  * The edge to knot map for NURBS meshes can be determined automatically. It is
    no longer needed to specify this in the NURBS mesh.
  * Added curve interpolation method for NURBS.
-  * Added new small miniapp for printing of shape functions of a KnotVector
  * See miniapps/nurbs for example meshes and miniapps.

 Discretization improvements
@@ -143,8 +178,6 @@ Linear and nonlinear solvers

 - Added HIP support to the PETSc and SUNDIALS interfaces.

- Efficient GPU-accelerated LOR assembly now supports surface meshes.
-
 New and updated examples and miniapps
 -------------------------------------
 - Added a new H(div) solver miniapp demonstrating the use of a matrix-free
@@ -58,7 +58,7 @@ project(mfem NONE)
 # Current version of MFEM, see also `makefile`.
 #   mfem_VERSION = (string)
 #   MFEM_VERSION = (int)   [automatically derived from mfem_VERSION]
-set(${PROJECT_NAME}_VERSION 4.6.1)
+set(${PROJECT_NAME}_VERSION 4.7.1)

 # Prohibit in-source build
 if (${PROJECT_SOURCE_DIR} STREQUAL ${PROJECT_BINARY_DIR})
@@ -87,10 +87,11 @@ if (MFEM_USE_STRUMPACK OR MFEM_USE_MUMPS)
  # Just needed to find the MPI_Fortran libraries to link with
  set(XSDK_ENABLE_Fortran ON)
 endif()
-# SUNDIALS, STRUMPACK, Ginkgo, RAJA and Umpire require C++14:
+# SUNDIALS, STRUMPACK, Ginkgo, Tribol, RAJA and Umpire require C++14:
 if ((MFEM_USE_SUNDIALS OR
     MFEM_USE_STRUMPACK OR
     MFEM_USE_GINKGO OR
+     MFEM_USE_TRIBOL OR
     MFEM_USE_RAJA OR
     MFEM_USE_UMPIRE) AND
    ("${CMAKE_CXX_STANDARD}" LESS "14"))
@@ -503,6 +504,15 @@ if (MFEM_USE_PARELAG)
   find_package(PARELAG REQUIRED)
 endif()

+# Tribol
+if (MFEM_USE_TRIBOL)
+  if (MFEM_USE_MPI)
+    find_package(Tribol REQUIRED tribol redecomp)
+  else()
+    message(FATAL_ERROR " *** Tribol requires that MPI be enabled.")
+  endif()
+endif()
+
 # Enzyme
 if (MFEM_USE_ENZYME)
   find_package(ENZYME REQUIRED)
@@ -548,7 +558,7 @@ set(MFEM_TPLS OPENMP HYPRE LAPACK BLAS SuperLUDist STRUMPACK METIS SuiteSparse
    SUNDIALS PETSC SLEPC MUMPS AXOM FMS CONDUIT Ginkgo GNUTLS GSLIB
    NETCDF MPFR PUMI HIOP POSIXCLOCKS MFEMBacktrace ZLIB OCCA CEED RAJA UMPIRE
    ADIOS2 CUSPARSE MKL_CPARDISO MKL_PARDISO AMGX CALIPER CODIPACK
-    BENCHMARK PARELAG MPI_CXX HIP HIPSPARSE MOONOLITH BLITZ ALGOIM ENZYME)
+    BENCHMARK PARELAG TRIBOL MPI_CXX HIP HIPSPARSE MOONOLITH BLITZ ALGOIM ENZYME)

 # Add all *_FOUND libraries in the variable TPL_LIBRARIES.
 set(TPL_LIBRARIES "")
@@ -151,7 +151,8 @@ The MFEM source code has the following structure:
  │   ├── solvers
  │   ├── spde
  │   ├── tools
-  │   └── toys
+  │   ├── toys
+  │   └── tribol
  └── tests
      ├── benchmarks
      ├── convergence
@@ -75,6 +75,8 @@ and miniapps. See https://glvis.org and https://mfem.org/building.

 Quick start with GNU make
 =========================
+See also: https://mfem.org/building
+
 Serial build:
   make serial -j 4

@@ -83,6 +85,7 @@ Parallel build:
   (build METIS 4 in ../metis-4.0 relative to mfem/)
   (build hypre in ../hypre relative to mfem/)
   make parallel -j 4
+(For METIS 5, see https://mfem.org/building/#parallel-build-using-metis-5)

 CUDA build:
   make cuda -j 4
@@ -116,6 +119,7 @@ Parallel build:
   mkdir <mfem-build-dir> ; cd <mfem-build-dir>
   cmake <mfem-source-dir> -DMFEM_USE_MPI=YES
   make -j 4
+(For METIS 5, see https://mfem.org/building/#parallel-build-using-metis-5)

 CUDA build:
   (this build requires CMake 3.8 or newer)
@@ -571,6 +575,11 @@ MFEM_USE_PARELAG = YES/NO
   use ParELAG. In fact, ParELAG is dependent on MFEM. Therefore, this option
   currently only concerns the miniapps.

+MFEM_USE_TRIBOL = YES/NO
+   Enables the miniapps that use the Tribol library. MFEM does not currently
+   use Tribol. In fact, Tribol is dependent on MFEM. Therefore, this option
+   currently only concerns the miniapps.
+
 MFEM_USE_ENZYME = YES/NO
   Enables automatic differentiation support through the LLVM plugin Enzyme.
   This requires the compiler to be set to clang (>=14.0.0). We also advise to
@@ -607,9 +616,13 @@ The specific libraries and their options are:
            HYPRE >= 2.20.0   (HYPRE built with '--enable-mixedint')
            HYPRE >= 2.22.1   (HYPRE built with CUDA)
            HYPRE >= 2.23.0   (HYPRE built with HIP)
+            HYPRE >= 2.31.0   (runtime selectable HYPRE execution on CPU/GPU)

 - METIS, used when MFEM_USE_METIS = YES. If using METIS 5, set
-  MFEM_USE_METIS_5 = YES (default is to use METIS 4).
+  MFEM_USE_METIS_5 = YES (default is to use METIS 4). For building instructions,
+  see the following:
+    - METIS 4.0.3: https://mfem.org/building/#parallel-mpi-version-of-mfem
+    - METIS 5.1.0: https://mfem.org/building/#parallel-build-using-metis-5
  URL: https://github.com/mfem/tpls (MFEM mirror, see above)
  Options: METIS_OPT, METIS_LIB.
  Versions: METIS 4.0.3 or 5.1.0.
@@ -857,6 +870,10 @@ The specific libraries and their options are:
  URL: https://github.com/LLNL/parelag
  Options: PARELAG_DIR, PARELAG_OPT, PARELAG_LIB.

+- Tribol, used when MFEM_USE_TRIBOL = YES.
+  URL: https://github.com/LLNL/Tribol
+  Options: TRIBOL_DIR, TRIBOL_OPT, TRIBOL_LIB.
+
 - Enzyme, used when MFEM_USE_ENZYME = YES. Requires LLVM/Clang >= 14.0.0.
  URL: https://github.com/EnzymeAD/Enzyme
  Options: ENZYME_DIR, ENZYME_OPT, ENZYME_LIB.
@@ -1001,6 +1018,7 @@ MFEM_USE_CALIPER
 MFEM_USE_FMS
 MFEM_USE_BENCHMARK
 MFEM_USE_PARELAG
+MFEM_USE_TRIBOL
 MFEM_USE_ENZYME

 The following options are CMake specific:
@@ -287,3 +287,7 @@ ENDIF()
 IF (DEFINED TPL_ENABLE_PARELAG)
  SET(MFEM_USE_PARELAG ${TPL_ENABLE_PARELAG} CACHE BOOL "Enable ParELAG" FORCE)
 ENDIF()
+
+IF (DEFINED TPL_ENABLE_TRIBOL)
+  SET(MFEM_USE_TRIBOL ${TPL_ENABLE_TRIBOL} CACHE BOOL "Enable Tribol" FORCE)
+ENDIF()
@@ -64,6 +64,7 @@ set(MFEM_USE_CALIPER @MFEM_USE_CALIPER@)
 set(MFEM_USE_ALGOIM @MFEM_USE_ALGOIM@)
 set(MFEM_USE_BENCHMARK @MFEM_USE_BENCHMARK@)
 set(MFEM_USE_PARELAG @MFEM_USE_PARELAG@)
+set(MFEM_USE_TRIBOL @MFEM_USE_TRIBOL@)
 set(MFEM_USE_ENZYME @MFEM_USE_ENZYME@)

 set(MFEM_CXX_COMPILER "@CMAKE_CXX_COMPILER@")
@@ -18,4 +18,13 @@ include(MfemCmakeUtilities)
 # Note: components are enabled based on the find_package() parameters.
 mfem_find_package(Axom AXOM AXOM_DIR "include" "" "lib" ""
  "Paths to headers required by Axom." "Libraries required by Axom."
-  ADD_COMPONENT Axom "include" axom/config.hpp "lib" axom)
+  ADD_COMPONENT core "include" axom/core.hpp "lib" axom_core
+  ADD_COMPONENT inlet "include" axom/inlet.hpp "lib" axom_inlet
+  ADD_COMPONENT klee "include" axom/klee.hpp "lib" axom_klee
+  ADD_COMPONENT lumberjack "include" axom/lumberjack.hpp "lib" axom_lumberjack
+  ADD_COMPONENT mint "include" axom/mint.hpp "lib" axom_mint
+  ADD_COMPONENT multimat "include" axom/multimat.hpp "lib" axom_multimat
+  ADD_COMPONENT quest "include" axom/quest.hpp "lib" axom_quest
+  ADD_COMPONENT sidre "include" axom/sidre.hpp "lib" axom_sidre
+  ADD_COMPONENT slam "include" axom/slam.hpp "lib" axom_slam
+  ADD_COMPONENT slic "include" axom/slic.hpp "lib" axom_slic)
@@ -36,7 +36,11 @@ include(MfemCmakeUtilities)
 mfem_find_package(Conduit CONDUIT CONDUIT_DIR
  "include;include/conduit" conduit.hpp "lib" conduit
  "Paths to headers required by Conduit." "Libraries required by Conduit."
+  ADD_COMPONENT blueprint
+    "include;include/conduit" conduit_blueprint.hpp "lib" conduit_blueprint
+  ADD_COMPONENT blueprint_mpi
+    "include;include/conduit" conduit_blueprint_mpi.hpp "lib" conduit_blueprint_mpi
  ADD_COMPONENT relay
    "include;include/conduit" conduit_relay.hpp "lib" conduit_relay
-  ADD_COMPONENT blueprint
-    "include;include/conduit" conduit_blueprint.hpp "lib" conduit_blueprint)
+  ADD_COMPONENT relay_mpi
+    "include;include/conduit" conduit_relay_mpi.hpp "lib" conduit_relay_mpi)
@@ -16,12 +16,22 @@
 #   - MUMPS_VERSION

 include(MfemCmakeUtilities)
+
+# Toggle which precision of MUMPS to use depending on the precision of MFEM.
+if (MFEM_USE_DOUBLE)
+  set(_mumps_header dmumps_c.h)
+  set(_mumps_lib dmumps)
+elseif(MFEM_USE_SINGLE)
+  set(_mumps_header smumps_c.h)
+  set(_mumps_lib smumps)
+endif()
+
 mfem_find_package(MUMPS MUMPS MUMPS_DIR
-  "include" dmumps_c.h "lib" dmumps
+  "include" ${_mumps_header} "lib" ${_mumps_lib}
  "Paths to headers required by MUMPS."
  "Libraries required by MUMPS."
-  ADD_COMPONENT mumps_common "include" dmumps_c.h "lib" mumps_common
-  ADD_COMPONENT pord "include" dmumps_c.h "lib" pord)
+  ADD_COMPONENT mumps_common "include" ${_mumps_header} "lib" mumps_common
+  ADD_COMPONENT pord "include" ${_mumps_header} "lib" pord)

 if (MUMPS_FOUND AND (NOT MUMPS_VERSION))
  try_run(MUMPS_VERSION_RUN_RESULT MUMPS_VERSION_COMPILE_RESULT
@@ -0,0 +1,22 @@
+# Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+# Defines the following variables:
+#   - TRIBOL_FOUND
+#   - TRIBOL_LIBRARIES
+#   - TRIBOL_INCLUDE_DIRS
+
+include(MfemCmakeUtilities)
+# Note: components are enabled based on the find_package() parameters.
+mfem_find_package(Tribol TRIBOL TRIBOL_DIR "include" tribol/config.hpp "lib" tribol
+  "Paths to headers required by Tribol." "Libraries required by Tribol."
+  ADD_COMPONENT redecomp
+    "include" redecomp/redecomp.hpp "lib" redecomp)
@@ -852,8 +852,8 @@ function(mfem_export_mk_files)
      MFEM_USE_CUDA MFEM_USE_HIP MFEM_USE_RAJA MFEM_USE_OCCA MFEM_USE_CEED
      MFEM_USE_CALIPER MFEM_USE_UMPIRE MFEM_USE_SIMD MFEM_USE_ADIOS2
      MFEM_USE_MKL_CPARDISO MFEM_USE_MKL_PARDISO MFEM_USE_ADFORWARD
-      MFEM_USE_CODIPACK MFEM_USE_BENCHMARK MFEM_USE_PARELAG MFEM_USE_MOONOLITH
-      MFEM_USE_ALGOIM MFEM_USE_ENZYME)
+      MFEM_USE_CODIPACK MFEM_USE_BENCHMARK MFEM_USE_PARELAG MFEM_USE_TRIBOL
+      MFEM_USE_MOONOLITH MFEM_USE_ALGOIM MFEM_USE_ENZYME)
  foreach(var ${CONFIG_MK_BOOL_VARS})
    if (${var})
      set(${var} YES)
@@ -120,6 +120,15 @@ constexpr real_t operator""_r(unsigned long long v)

 // Check dependencies:

+// Define MFEM_MPI_REAL_T to be the appropriate MPI real type
+#ifdef MFEM_USE_MPI
+#ifdef MFEM_USE_SINGLE
+#define MFEM_MPI_REAL_T MPI_FLOAT
+#elif defined MFEM_USE_DOUBLE
+#define MFEM_MPI_REAL_T MPI_DOUBLE
+#endif
+#endif
+
 // Options that require MPI
 #ifndef MFEM_USE_MPI
 #ifdef MFEM_USE_SUPERLU
@@ -65,6 +65,7 @@ MFEM_USE_ADFORWARD     = @MFEM_USE_ADFORWARD@
 MFEM_USE_CODIPACK      = @MFEM_USE_CODIPACK@
 MFEM_USE_BENCHMARK     = @MFEM_USE_BENCHMARK@
 MFEM_USE_PARELAG       = @MFEM_USE_PARELAG@
+MFEM_USE_TRIBOL        = @MFEM_USE_TRIBOL@
 MFEM_USE_ENZYME        = @MFEM_USE_ENZYME@

 # Compiler, compile options, and link options
@@ -67,6 +67,7 @@ option(MFEM_USE_ADFORWARD "Enable forward mode for AD" OFF)
 option(MFEM_USE_CODIPACK "Enable automatic differentiation (AD) using CoDiPack" OFF)
 option(MFEM_USE_BENCHMARK "Enable Google Benchmark" OFF)
 option(MFEM_USE_PARELAG "Enable ParELAG" OFF)
+option(MFEM_USE_TRIBOL "Enable Tribol" OFF)
 option(MFEM_USE_ENZYME "Enable Enzyme" OFF)

 # Optional overrides for autodetected MPIEXEC and MPIEXEC_NUMPROC_FLAG
@@ -212,8 +213,15 @@ set(CONDUIT_DIR "${MFEM_DIR}/../conduit" CACHE PATH

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

 set(PUMI_DIR "${MFEM_DIR}/../pumi-2.1.0" CACHE STRING
    "Directory where PUMI is installed")
@@ -250,6 +258,10 @@ set(PARELAG_INCLUDE_DIRS "${PARELAG_DIR}/src;${PARELAG_DIR}/build/src" CACHE
 set(PARELAG_LIBRARIES "${PARELAG_DIR}/build/src/libParELAG.a" CACHE STRING
    "The ParELAG library.")

+set(TRIBOL_DIR "${MFEM_DIR}/../tribol" CACHE PATH "Path to Tribol")
+set(Tribol_REQUIRED_PACKAGES "Axom/core/mint/slam/slic" CACHE STRING 
+    "Additional packages required by Tribol")
+
 set(BLAS_INCLUDE_DIRS "" CACHE STRING "Path to BLAS headers.")
 set(BLAS_LIBRARIES "" CACHE STRING "The BLAS library.")
 set(LAPACK_INCLUDE_DIRS "" CACHE STRING "Path to LAPACK headers.")
@@ -167,8 +167,21 @@ MFEM_USE_ADFORWARD     = NO
 MFEM_USE_CODIPACK      = NO
 MFEM_USE_BENCHMARK     = NO
 MFEM_USE_PARELAG       = NO
+MFEM_USE_TRIBOL        = NO
 MFEM_USE_ENZYME        = NO

+# Process MFEM_PRECISION -> MFEM_USE_SINGLE, MFEM_USE_DOUBLE
+ifneq ($(filter double Double DOUBLE,$(MFEM_PRECISION)),)
+   MFEM_USE_DOUBLE = YES
+   MFEM_USE_SINGLE = NO
+else ifneq ($(filter single Single SINGLE,$(MFEM_PRECISION)),)
+   MFEM_USE_DOUBLE = NO
+   MFEM_USE_SINGLE = YES
+else ifeq ($(MAKECMDGOALS),config)
+   $(error Invalid floating-point precision: \
+     MFEM_PRECISION = $(MFEM_PRECISION))
+endif
+
 # MPI library compile and link flags
 # These settings are used only when building MFEM with MPI + HIP
 ifeq ($(MFEM_USE_MPI)$(MFEM_USE_HIP),YESYES)
@@ -318,13 +331,13 @@ MPI_FORTRAN_LIB = -lmpifort
 # MUMPS library configuration
 MUMPS_DIR = @MFEM_DIR@/../MUMPS_5.5.0
 MUMPS_OPT = -I$(MUMPS_DIR)/include
-MUMPS_LIB = $(XLINKER)-rpath,$(MUMPS_DIR)/lib -L$(MUMPS_DIR)/lib \
- -lmumps_common -lpord $(SCALAPACK_LIB) $(LAPACK_LIB) $(MPI_FORTRAN_LIB)
+MUMPS_LIB = $(XLINKER)-rpath,$(MUMPS_DIR)/lib -L$(MUMPS_DIR)/lib
 ifeq ($(MFEM_USE_SINGLE),YES)
   MUMPS_LIB += -lsmumps
 else
   MUMPS_LIB += -ldmumps
 endif
+MUMPS_LIB += -lmumps_common -lpord $(SCALAPACK_LIB) $(LAPACK_LIB) $(MPI_FORTRAN_LIB)

 # STRUMPACK library configuration
 STRUMPACK_DIR = @MFEM_DIR@/../STRUMPACK-build
@@ -375,7 +388,7 @@ GINKGO_LIB = $(XLINKER)-rpath,$(GINKGO_LINK_LIB_DIR) -L$(GINKGO_LINK_LIB_DIR)\
 # AmgX library configuration
 AMGX_DIR = @MFEM_DIR@/../amgx
 AMGX_OPT = -I$(AMGX_DIR)/include
-AMGX_LIB = -lcusparse -lcusolver -lcublas -lnvToolsExt -L$(AMGX_DIR)/lib -lamgx
+AMGX_LIB = -L$(AMGX_DIR)/lib -lamgx -lcusparse -lcusolver -lcublas -lnvToolsExt

 # GnuTLS library configuration
 GNUTLS_OPT =
@@ -576,6 +589,16 @@ PARELAG_DIR = @MFEM_DIR@/../parelag
 PARELAG_OPT = -I$(PARELAG_DIR)/src -I$(PARELAG_DIR)/build/src
 PARELAG_LIB = -L$(PARELAG_DIR)/build/src -lParELAG

+# Tribol library configuration
+ifeq ($(MFEM_USE_TRIBOL),YES)
+   BASE_FLAGS = -std=c++14
+endif
+AXOM_DIR = @MFEM_DIR@/../axom
+TRIBOL_DIR = @MFEM_DIR@/../tribol
+TRIBOL_OPT = -I$(TRIBOL_DIR)/include -I$(AXOM_DIR)/include
+TRIBOL_LIB = -L$(TRIBOL_DIR)/lib -ltribol -lredecomp -L$(AXOM_DIR)/lib -laxom_mint\
+   -laxom_slam -laxom_slic -laxom_core
+
 # Enzyme configuration

 # If you want to enable automatic differentiation at compile time, use the
@@ -110,4 +110,4 @@ config-mk:

 clean:
 	rm -f $(CONFIG_HPP) $(CONFIG_MK) sample-runs-build.log
-	rm -f $(GHV) $(GHV).out $(GMV) $(GMV).out
+	rm -f $(GHV) $(GHV).out $(GMV) $(GMV).out *.dSYM
@@ -92,4 +92,5 @@ vertices
 -0.70710678 -0.70710678
 0 -1
 0.70710678 -0.70710678
+
 mfem_mesh_end
@@ -48,7 +48,7 @@ PROJECT_NAME           = MFEM
 # could be handy for archiving the generated documentation or if some version
 # control system is used.

-PROJECT_NUMBER         = v4.6.1
+PROJECT_NUMBER         = v4.7.1

 # Using the PROJECT_BRIEF tag one can provide an optional one line description
 # for a project that appears at the top of each page and should give viewer a
@@ -987,6 +987,7 @@ INPUT                  = @MFEM_SOURCE_DIR@/doc/CodeDocumentation.dox \
                         @MFEM_SOURCE_DIR@/miniapps/solvers \
                         @MFEM_SOURCE_DIR@/miniapps/tools \
                         @MFEM_SOURCE_DIR@/miniapps/toys \
+                         @MFEM_SOURCE_DIR@/miniapps/tribol \
                         @MFEM_SOURCE_DIR@/miniapps/spde \
                         @MFEM_SOURCE_DIR@/miniapps/dpg \
                         @MFEM_SOURCE_DIR@/miniapps/dpg/util
@@ -110,9 +110,13 @@ namespace mfem {
 * - <a class="el" href="ex35p_8cpp_source.html">Example 35p</a>: parallel multi-domain damped harmonic oscillators
 * - <a class="el" href="ex36_8cpp_source.html">Example 36</a>: Proximal Galerkin FEM for the obstacle problem
 * - <a class="el" href="ex36p_8cpp_source.html">Example 36p</a>: parallel Proximal Galerkin FEM for the obstacle problem
- * - <a class="el" href="ex37_8cpp_source.html">Example 37</a>: Topology optimization
+ * - <a class="el" href="ex37_8cpp_source.html">Example 37</a>: topology optimization
 * - <a class="el" href="ex37p_8cpp_source.html">Example 37p</a>: parallel topology optimization
 * - <a class="el" href="ex38_8cpp_source.html">Example 38</a>: cut-surface and cut-volume integration
+ * - <a class="el" href="ex39_8cpp_source.html">Example 39</a>: named mesh attributes
+ * - <a class="el" href="ex39p_8cpp_source.html">Example 39p</a>: parallel named mesh attributes
+ * - <a class="el" href="ex40_8cpp_source.html">Example 40</a>: eikonal equation
+ * - <a class="el" href="ex40p_8cpp_source.html">Example 40p</a>: parallel eikonal equation
 *
 * <H4>AmgX Examples</H4>
 * - Variants of Examples
@@ -214,6 +218,8 @@ namespace mfem {
 * - <a class="el" href="miniapps_2performance_2ex1_8cpp_source.html">HPC Example 1</a>: high-performance nodal H1 FEM for the Laplace problem
 * - <a class="el" href="miniapps_2performance_2ex1p_8cpp_source.html">HPC Example 1p</a>: high-performance parallel nodal H1 FEM for the Laplace problem
 * - <a class="el" href="generate__random__field_8cpp_source.html">SPDE Solvers</a>: SPDE solver random field generation
+ * - <a class="el" href="contact-patch-test_8cpp_source.html">Contact</a>: mortar contact patch test for elasticity
+ * - <a class="el" href="multidomain_8cpp_source.html">Multidomain miniapp</a>: Multidomain and Submesh demonstration miniapp
 * - <a class="el" href="pdiffusion_8cpp_source.html">DPG Diffusion example</a>: DPG formulation for the diffusion problem
 * - <a class="el" href="pmaxwell_8cpp_source.html">DPG Maxwell example</a>: DPG formulation for the indefinite Maxwell problem
 * - <a class="el" href="lor__elast_8cpp_source.html">LOR Elasticity</a>: solve linear elasticity with LOR preconditioning on GPUs
@@ -46,7 +46,7 @@ class DoxygenAwesomeDarkModeToggle extends HTMLElement {
            DoxygenAwesomeDarkModeToggle.onSystemPreferenceChanged()
        })
        // Update the color scheme when the tab is made visible again.
-        // It is possible that the appearance was changed in another tab 
+        // It is possible that the appearance was changed in another tab
        // while this tab was in the background.
        document.addEventListener("visibilitychange", visibilityState => {
            if (document.visibilityState === 'visible') {
@@ -97,7 +97,7 @@ class DoxygenAwesomeDarkModeToggle extends HTMLElement {
     * @returns `true` for dark-mode, `false` for light-mode user preference
     */
    static get userPreference() {
-        return (!DoxygenAwesomeDarkModeToggle.systemPreference && localStorage.getItem(DoxygenAwesomeDarkModeToggle.prefersDarkModeInLightModeKey)) || 
+        return (!DoxygenAwesomeDarkModeToggle.systemPreference && localStorage.getItem(DoxygenAwesomeDarkModeToggle.prefersDarkModeInLightModeKey)) ||
        (DoxygenAwesomeDarkModeToggle.systemPreference && !localStorage.getItem(DoxygenAwesomeDarkModeToggle.prefersLightModeInDarkModeKey))
    }

@@ -45,6 +45,7 @@ list(APPEND ALL_EXE_SRCS
  ex37.cpp
  ex38.cpp
  ex39.cpp
+  ex40.cpp
  )

 if (MFEM_USE_MPI)
@@ -87,6 +88,7 @@ if (MFEM_USE_MPI)
    ex36p.cpp
    ex37p.cpp
    ex39p.cpp
+    ex40p.cpp
  )
 endif()

@@ -146,10 +148,10 @@ if (MFEM_ENABLE_TESTING)
  # Add CUDA/HIP tests.
  set(DEVICE_EXAMPLES
      # serial examples with device support:
-      ex1 ex3 ex4 ex5 ex6 ex9 ex22 ex24 ex25 ex26 ex34
+      ex1 ex3 ex4 ex5 ex6 ex9 ex14 ex22 ex24 ex25 ex26 ex34
      # parallel examples with device support:
-      ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex9p ex13p ex22p ex24p ex25p ex26p
-      ex34p ex35p)
+      ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex9p ex13p ex14p ex22p ex24p ex25p
+      ex26p ex34p ex35p)
  set(MFEM_TEST_DEVICE)
  if (MFEM_USE_CUDA)
    set(MFEM_TEST_DEVICE "cuda")
@@ -159,6 +161,11 @@ if (MFEM_ENABLE_TESTING)
  if (MFEM_TEST_DEVICE)
    foreach(TEST_NAME ${DEVICE_EXAMPLES})
      set(THIS_TEST_OPTIONS "-no-vis" "-d" "${MFEM_TEST_DEVICE}")
+      if (${TEST_NAME} MATCHES "ex14p")
+        list(APPEND THIS_TEST_OPTIONS "-rs" "2" "-rp" "0" "-pa")
+      elseif (${TEST_NAME} MATCHES "ex14")
+        list(APPEND THIS_TEST_OPTIONS "-r" "2" "-pa")
+      endif()
      if (NOT (${TEST_NAME} MATCHES ".*p$"))
        add_test(NAME ${TEST_NAME}_${MFEM_TEST_DEVICE}_ser
          COMMAND ${TEST_NAME} ${THIS_TEST_OPTIONS})
@@ -646,10 +646,7 @@ real_t HyperelasticOperator::ElasticEnergy(const ParGridFunction &x) const

 real_t HyperelasticOperator::KineticEnergy(const ParGridFunction &v) const
 {
-   real_t loc_energy = 0.5*M.InnerProduct(v, v);
-   real_t energy;
-   MPI_Allreduce(&loc_energy, &energy, 1, MPITypeMap<real_t>::mpi_type,
-                 MPI_SUM, fespace.GetComm());
+   real_t energy = 0.5*M.ParInnerProduct(v, v);
   return energy;
 }

@@ -18,10 +18,12 @@
 //               ex14 -m ../data/amr-quad.mesh -r 3
 //               ex14 -m ../data/amr-hex.mesh
 //               ex14 -m ../data/fichera-amr.mesh
+//               ex14 -pa -r 1 -o 3
+//               ex14 -pa -r 1 -o 3 -m ../data/fichera.mesh
 //
 // Device sample runs:
-//               ex14 -pa -d cuda -o 3
-//               ex14 -pa -d cuda -o 3 -m ../data/fichera.mesh
+//               ex14 -pa -r 2 -d cuda -o 3
+//               ex14 -pa -r 2 -d cuda -o 3 -m ../data/fichera.mesh
 //
 // Description:  This example code demonstrates the use of MFEM to define a
 //               discontinuous Galerkin (DG) finite element discretization of
@@ -119,7 +121,8 @@ int main(int argc, char *argv[])

   // 5. Define a finite element space on the mesh. Here we use discontinuous
   //    finite elements of the specified order >= 0.
-   DG_FECollection fec(order, dim, BasisType::GaussLobatto);
+   const auto bt = pa ? BasisType::GaussLobatto : BasisType::GaussLegendre;
+   DG_FECollection fec(order, dim, bt);
   FiniteElementSpace fespace(&mesh, &fec);
   cout << "Number of unknowns: " << fespace.GetVSize() << endl;

@@ -160,21 +163,15 @@ int main(int argc, char *argv[])

   // 9. Define a simple symmetric Gauss-Seidel preconditioner and use it to
   //    solve the system Ax=b with PCG in the symmetric case, and GMRES in the
-   //    non-symmetric one. (Note that tolerances are squared: 1e-24 corresponds
-   //    to a relative tolerance of 1e-12).
+   //    non-symmetric one. (Note that tolerances are squared: 1e-12 corresponds
+   //    to a relative tolerance of 1e-6).
   //
   //    If MFEM was compiled with SuiteSparse, use UMFPACK to solve the system.
   if (pa)
   {
-      const Operator &A = a;
-      if (sigma == -1.0)
-      {
-         CG(A, b, x, 1, 500, 1e-24, 0.0);
-      }
-      else
-      {
-         MFEM_ABORT("The case of PA with sigma != -1 is not yet supported.");
-      }
+      MFEM_VERIFY(sigma == -1.0,
+                  "The case of PA with sigma != -1 is not yet supported.");
+      CG(a, b, x, 1, 500, 1e-12, 0.0);
   }
   else
   {
@@ -183,11 +180,11 @@ int main(int argc, char *argv[])
      GSSmoother M(A);
      if (sigma == -1.0)
      {
-         PCG(A, M, b, x, 1, 500, 1e-24, 0.0);
+         PCG(A, M, b, x, 1, 500, 1e-12, 0.0);
      }
      else
      {
-         GMRES(A, M, b, x, 1, 500, 10, 1e-24, 0.0);
+         GMRES(A, M, b, x, 1, 500, 10, 1e-12, 0.0);
      }
 #else
      UMFPackSolver umf_solver;
@@ -17,10 +17,12 @@
 //               mpirun -np 4 ex14p -m ../data/inline-segment.mesh -rs 5
 //               mpirun -np 4 ex14p -m ../data/amr-quad.mesh -rs 3
 //               mpirun -np 4 ex14p -m ../data/amr-hex.mesh
+//               mpirun -np 4 ex14p -pa -rs 1 -rp 0 -o 3
+//               mpirun -np 4 ex14p -pa -rs 1 -rp 0 -m ../data/fichera.mesh -o 3
 //
 // Device sample runs:
-//               mpirun -np 4 ex14p -pa -d cuda -o 3
-//               mpirun -np 4 ex14p -pa -d cuda -m ../data/fichera.mesh -o 3
+//               mpirun -np 4 ex14p -pa -rs 2 -rp 0 -d cuda -o 3
+//               mpirun -np 4 ex14p -pa -rs 2 -rp 0 -d cuda -m ../data/fichera.mesh -o 3
 //
 // Description:  This example code demonstrates the use of MFEM to define a
 //               discontinuous Galerkin (DG) finite element discretization of
@@ -173,7 +175,8 @@ int main(int argc, char *argv[])

   // 6. Define a parallel finite element space on the parallel mesh. Here we
   //    use discontinuous finite elements of the specified order >= 0.
-   DG_FECollection fec(order, dim, BasisType::GaussLobatto);
+   const auto bt = pa ? BasisType::GaussLobatto : BasisType::GaussLegendre;
+   DG_FECollection fec(order, dim, bt);
   ParFiniteElementSpace fespace(&pmesh, &fec);
   HYPRE_BigInt size = fespace.GlobalTrueVSize();
   if (Mpi::Root())
@@ -39,8 +39,8 @@ private:
   // Base Nonlinear Form
   std::unique_ptr<NonlinearForm> nonlinearForm;
   // element-wise inverse mass matrix
-   std::vector<DenseMatrix> invmass; // local scalar inverse mass.
-   std::vector<DenseMatrix> weakdiv; // local weakdivergence. Trial space is ByDim.
+   std::vector<DenseMatrix> invmass; // local scalar inverse mass
+   std::vector<DenseMatrix> weakdiv; // local weak divergence (trial space ByDim)
   // global maximum characteristic speed. Updated by form integrators
   mutable real_t max_char_speed;
   // auxiliary variable used in Mult
@@ -169,9 +169,9 @@ void DGHyperbolicConservationLaws::Mult(const Vector &x, Vector &y) const
 {
   // 0. Reset wavespeed computation before operator application.
   formIntegrator->ResetMaxCharSpeed();
-   // 1. Apply Nonlinear form to obtain an axiliary result
+   // 1. Apply Nonlinear form to obtain an auxiliary result
   //         z = - <F̂(u_h,n), [[v]]>_e
-   //    If weak-divergencee is not preassembled, we also have weak-divergence
+   //    If weak-divergence is not preassembled, we also have weak-divergence
   //         z = - <F̂(u_h,n), [[v]]>_e + (F(u_h), ∇v)
   nonlinearForm->Mult(x, z);
   if (!weakdiv.empty()) // if weak divergence is pre-assembled
@@ -19,8 +19,11 @@
 //    ex33 -m ../data/amr-quad.mesh -ver -alpha 2.6 -o 2 -r 2
 //    ex33 -m ../data/inline-hex.mesh -ver -alpha 0.3 -o 2 -r 1
 //
-//  Note: the analytic solution to this problem is u = ∏_{i=0}^{dim-1} sin(π x_i)
-//        for all alpha.
+//  Note: The manufactured solution used in this problem is
+//
+//            u = ∏_{i=0}^{dim-1} sin(π x_i) ,
+//
+//        regardless of the value of alpha.
 //
 // Description:
 //
@@ -114,7 +117,8 @@ int main(int argc, char *argv[])
                  "Enable or disable GLVis visualization.");
   args.AddOption(&verification, "-ver", "--verification", "-no-ver",
                  "--no-verification",
-                  "Use sinusoidal function (f) for analytic comparison.");
+                  "Use sinusoidal function (f) for manufactured "
+                  "solution test.");
   args.Parse();
   if (!args.Good())
   {
@@ -163,7 +167,7 @@ int main(int argc, char *argv[])
   // 5. Define a finite element space on the mesh.
   H1_FECollection fec(order, dim);
   FiniteElementSpace fespace(&mesh, &fec);
-   cout << "Number of finite element unknowns: "
+   cout << "Number of degrees of freedom: "
        << fespace.GetTrueVSize() << endl;

   // 6. Determine the list of true (i.e. conforming) essential boundary dofs.
@@ -379,29 +383,29 @@ int main(int argc, char *argv[])
      FunctionCoefficient sol(solution);
      real_t l2_error = u.ComputeL2Error(sol);

-      string analytic_solution,expected_mesh;
+      string manufactured_solution,expected_mesh;
      switch (dim)
      {
         case 1:
-            analytic_solution = "sin(π x)";
+            manufactured_solution = "sin(π x)";
            expected_mesh = "inline_segment.mesh";
            break;
         case 2:
-            analytic_solution = "sin(π x) sin(π y)";
+            manufactured_solution = "sin(π x) sin(π y)";
            expected_mesh = "inline_quad.mesh";
            break;
         default:
-            analytic_solution = "sin(π x) sin(π y) sin(π z)";
+            manufactured_solution = "sin(π x) sin(π y) sin(π z)";
            expected_mesh = "inline_hex.mesh";
            break;
      }

      mfem::out << "\n" << string(80,'=')
                << "\n\nSolution Verification in "<< dim << "D \n\n"
-                << "Analytic solution : " << analytic_solution << "\n"
-                << "Expected mesh     : " << expected_mesh <<"\n"
-                << "Your mesh         : " << mesh_file << "\n"
-                << "L2 error          : " << l2_error << "\n\n"
+                << "Manufactured solution : " << manufactured_solution << "\n"
+                << "Expected mesh         : " << expected_mesh <<"\n"
+                << "Your mesh             : " << mesh_file << "\n"
+                << "L2 error              : " << l2_error << "\n\n"
                << string(80,'=') << endl;
   }

@@ -131,7 +131,7 @@ void RationalApproximation_AAA(const Vector &val, const Vector &pt,
      }

 #ifdef MFEM_USE_LAPACK
-      DenseMatrixSVD svd(Am,false,true);
+      DenseMatrixSVD svd(Am,'N','A');
      svd.Eval(Am);
      DenseMatrix &v = svd.RightSingularvectors();
      v.GetRow(k,w);
@@ -346,7 +346,7 @@ void ComputePartialFractionApproximation(real_t & alpha,
   }
   else
   {
-      if (abs(alpha - 0.5) > eps && print_warning)
+      if (abs(alpha - 0.5) > eps)
      {
         alpha = 0.5;
      }
@@ -368,6 +368,8 @@ void ComputePartialFractionApproximation(real_t & alpha,


   return;
+#else
+   MFEM_CONTRACT_VAR(print_warning);
 #endif

   Vector x(npoints);
@@ -19,8 +19,11 @@
 //    mpirun -np 4 ex33p -m ../data/amr-quad.mesh -ver -alpha 2.6 -o 2 -r 2
 //    mpirun -np 4 ex33p -m ../data/inline-hex.mesh -ver -alpha 0.3 -o 2 -r 1

-//  Note: the analytic solution to this problem is u = ∏_{i=0}^{dim-1} sin(π x_i)
-//        for all alpha.
+//  Note: The manufactured solution used in this problem is
+//
+//            u = ∏_{i=0}^{dim-1} sin(π x_i) ,
+//
+//        regardless of the value of alpha.
 //
 // Description:
 //
@@ -120,7 +123,8 @@ int main(int argc, char *argv[])
                  "Enable or disable GLVis visualization.");
   args.AddOption(&verification, "-ver", "--verification", "-no-ver",
                  "--no-verification",
-                  "Use sinusoidal function (f) for analytic comparison.");
+                  "Use sinusoidal function (f) for manufactured "
+                  "solution test.");
   args.Parse();
   if (!args.Good())
   {
@@ -180,10 +184,11 @@ int main(int argc, char *argv[])
   // 5. Define a finite element space on the mesh.
   H1_FECollection fec(order, dim);
   ParFiniteElementSpace fespace(&pmesh, &fec);
+   HYPRE_BigInt size = fespace.GlobalTrueVSize();
   if (Mpi::Root())
   {
-      cout << "Number of finite element unknowns: "
-           << fespace.GetTrueVSize() << endl;
+      cout << "Number of degrees of freedom: "
+           << size << endl;
   }

   // 6. Determine the list of true (i.e. conforming) essential boundary dofs.
@@ -223,7 +228,7 @@ int main(int argc, char *argv[])
   if (verification)
   {
      // This statement is only relevant for the verification of the code. It
-      // uses a different f such that an analytic solution is known and easy
+      // uses a different f such that an manufactured solution is known and easy
      // to compare with the numerical one. The FPDE becomes:
      // (-Δ)^α u = (2\pi ^2)^α sin(\pi x) sin(\pi y) on [0,1]^2
      // -> u(x,y) = sin(\pi x) sin(\pi y)
@@ -415,29 +420,29 @@ int main(int argc, char *argv[])

      if (Mpi::Root())
      {
-         string analytic_solution,expected_mesh;
+         string manufactured_solution,expected_mesh;
         switch (dim)
         {
            case 1:
-               analytic_solution = "sin(π x)";
+               manufactured_solution = "sin(π x)";
               expected_mesh = "inline_segment.mesh";
               break;
            case 2:
-               analytic_solution = "sin(π x) sin(π y)";
+               manufactured_solution = "sin(π x) sin(π y)";
               expected_mesh = "inline_quad.mesh";
               break;
            default:
-               analytic_solution = "sin(π x) sin(π y) sin(π z)";
+               manufactured_solution = "sin(π x) sin(π y) sin(π z)";
               expected_mesh = "inline_hex.mesh";
               break;
         }

         mfem::out << "\n" << string(80,'=')
                   << "\n\nSolution Verification in "<< dim << "D \n\n"
-                   << "Analytic solution : " << analytic_solution << "\n"
-                   << "Expected mesh     : " << expected_mesh <<"\n"
-                   << "Your mesh         : " << mesh_file << "\n"
-                   << "L2 error          : " << l2_error << "\n\n"
+                   << "Manufactured solution : " << manufactured_solution << "\n"
+                   << "Expected mesh         : " << expected_mesh <<"\n"
+                   << "Your mesh             : " << mesh_file << "\n"
+                   << "L2 error              : " << l2_error << "\n\n"
                   << string(80,'=') << endl;
      }
   }
@@ -199,7 +199,6 @@ public:
      {
         mesh->GetElementTransformation(elem, &Tr);
         MFIRs.GetSurfaceIntegrationRule(Tr, ir);
-         Vector w;
         MFIRs.GetSurfaceWeights(Tr, ir, w);
         SurfaceWeights.SetCol(elem, w);

@@ -0,0 +1,374 @@
+//                                MFEM Example 40
+//
+// Compile with: make ex40
+//
+// Sample runs: ex40 -step 10 -gr 2.0
+//              ex40 -step 10 -gr 2.0 -o 3 -r 1
+//              ex40 -step 10 -gr 2.0 -r 4 -m ../data/l-shape.mesh
+//              ex40 -step 10 -gr 2.0 -r 2 -m ../data/fichera.mesh
+//
+// Description: This example code demonstrates how to use MFEM to solve the
+//              eikonal equation,
+//
+//                      |∇𝑢| = 1 in Ω,  𝑢 = g on ∂Ω.
+//
+//              The solution of this problem coincides with the unique optimum of
+//              the nonlinear program
+//
+//                   maximize ∫_Ω 𝑢 d𝑥 subject to |∇𝑢| ≤ 1,   𝑢 = g on Ω,      (⋆)
+//
+//              which is the foundation for method implemented below.
+//
+//              Following the proximal Galerkin methodology [1] (see also Example
+//              36), we construct a Legendre function for the unit ball
+//              𝐵₁ := {𝑥 ∈ Rⁿ | |𝑥| < 1}. Our choice is the Hellinger entropy,
+//
+//                    h(𝑥) = −( 1 − |𝑥|² )^{1/2},
+//
+//              although other choices are possible, each leading to a slightly
+//              different algorithm. We then adaptively regularize the optimization
+//              problem (⋆) with the Bregman divergence of the Hellinger entropy,
+//
+//                 maximize  ∫_Ω 𝑢 d𝑥 - αₖ⁻¹ Dₕ(∇𝑢,∇𝑢ₖ₋₁)  subject to  𝑢 = g on Ω.
+//
+//              This results in a sequence of functions ( 𝜓ₖ , 𝑢ₖ ),
+//
+//                      𝑢ₖ → 𝑢,    𝜓ₖ/|𝜓ₖ| → ∇𝑢    as k → \infty,
+//
+//              defined by the nonlinear saddle-point problems
+//
+//               Find 𝜓ₖ ∈ H(div,Ω) and 𝑢ₖ ∈ L²(Ω) such that
+//               ( Zₖ(𝜓ₖ) , τ ) + ( 𝑢ₖ , ∇⋅τ ) = ⟨ g , τ⋅n ⟩         ∀ τ ∈ H(div,Ω)
+//               ( ∇⋅𝜓ₖ , v )                 = ( ∇⋅𝜓ₖ₋₁ - 1 , v )   ∀ v ∈ L²(Ω)
+//
+//              where Zₖ(𝜓) := ∇h⁻¹(αₖ 𝜓) = 𝜓 / ( αₖ⁻² + |𝜓|² )^{1/2} and step size
+//              αₖ > 0. These saddle-point problems are solved using a damped Newton's
+//              method. This example assumes that g = 0 and allows the step size to
+//              grow geometrically, αₖ = α₀rᵏ, where r ≥ 1 is the growth rate.
+//
+//              [1] Keith, B. and Surowiec, T. (2023) Proximal Galerkin: A structure-
+//                  preserving finite element method for pointwise bound constraints.
+//                  arXiv:2307.12444 [math.NA]
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+class ZCoefficient : public VectorCoefficient
+{
+protected:
+   GridFunction *psi;
+   real_t alpha;
+
+public:
+   ZCoefficient(int vdim, GridFunction &psi_, real_t alpha_ = 1.0)
+      : VectorCoefficient(vdim), psi(&psi_), alpha(alpha_) { }
+
+   virtual void Eval(Vector &V, ElementTransformation &T,
+                     const IntegrationPoint &ip);
+   void SetAlpha(real_t alpha_) { alpha = alpha_; }
+};
+
+class DZCoefficient : public MatrixCoefficient
+{
+protected:
+   GridFunction *psi;
+   real_t alpha;
+
+public:
+   DZCoefficient(int height, GridFunction &psi_, real_t alpha_ = 1.0)
+      : MatrixCoefficient(height),  psi(&psi_), alpha(alpha_) { }
+
+   virtual void Eval(DenseMatrix &K, ElementTransformation &T,
+                     const IntegrationPoint &ip);
+   void SetAlpha(real_t alpha_) { alpha = alpha_; }
+};
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file = "../data/star.mesh";
+   int order = 1;
+   int max_it = 5;
+   int ref_levels = 3;
+   real_t alpha = 1.0;
+   real_t growth_rate = 1.0;
+   real_t newton_scaling = 0.9;
+   real_t tichonov = 1e-1;
+   real_t tol = 1e-4;
+   bool visualization = true;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&ref_levels, "-r", "--refs",
+                  "Number of h-refinements.");
+   args.AddOption(&max_it, "-mi", "--max-it",
+                  "Maximum number of iterations");
+   args.AddOption(&tol, "-tol", "--tol",
+                  "Stopping criteria based on the difference between"
+                  "successive solution updates");
+   args.AddOption(&alpha, "-step", "--step",
+                  "Initial size alpha");
+   args.AddOption(&growth_rate, "-gr", "--growth-rate",
+                  "Growth rate of the step size alpha");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   // 2. Read the mesh from the mesh file.
+   Mesh mesh(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+   int sdim = mesh.SpaceDimension();
+
+   MFEM_ASSERT(mesh.bdr_attributes.Size(),
+               "This example does not currently support meshes"
+               " without boundary attributes."
+              )
+
+   // 3. Postprocess the mesh.
+   // 3A. Refine the mesh to increase the resolution.
+   for (int l = 0; l < ref_levels; l++)
+   {
+      mesh.UniformRefinement();
+   }
+
+   // 3B. Interpolate the geometry after refinement to control geometry error.
+   // NOTE: Minimum second-order interpolation is used to improve the accuracy.
+   int curvature_order = max(order,2);
+   mesh.SetCurvature(curvature_order);
+
+   // 4. Define the necessary finite element spaces on the mesh.
+   RT_FECollection RTfec(order, dim);
+   FiniteElementSpace RTfes(&mesh, &RTfec);
+
+   L2_FECollection L2fec(order, dim);
+   FiniteElementSpace L2fes(&mesh, &L2fec);
+
+   cout << "Number of H(div) dofs: "
+        << RTfes.GetTrueVSize() << endl;
+   cout << "Number of L² dofs: "
+        << L2fes.GetTrueVSize() << endl;
+
+   // 5. Define the offsets for the block matrices
+   Array<int> offsets(3);
+   offsets[0] = 0;
+   offsets[1] = RTfes.GetVSize();
+   offsets[2] = L2fes.GetVSize();
+   offsets.PartialSum();
+
+   BlockVector x(offsets), rhs(offsets);
+   x = 0.0; rhs = 0.0;
+
+   // 6. Define the solution vectors as a finite element grid functions
+   //    corresponding to the fespaces.
+   GridFunction u_gf, delta_psi_gf;
+   delta_psi_gf.MakeRef(&RTfes,x,offsets[0]);
+   u_gf.MakeRef(&L2fes,x,offsets[1]);
+
+   GridFunction psi_old_gf(&RTfes);
+   GridFunction psi_gf(&RTfes);
+   GridFunction u_old_gf(&L2fes);
+
+   // 7. Define initial guesses for the solution variables.
+   delta_psi_gf = 0.0;
+   psi_gf = 0.0;
+   u_gf = 0.0;
+   psi_old_gf = psi_gf;
+   u_old_gf = u_gf;
+
+   // 8. Prepare for glvis output.
+   char vishost[] = "localhost";
+   int  visport   = 19916;
+   socketstream sol_sock;
+   if (visualization)
+   {
+      sol_sock.open(vishost,visport);
+      sol_sock.precision(8);
+   }
+
+   // 9. Coefficients to be used later.
+   ConstantCoefficient neg_one(-1.0);
+   ConstantCoefficient zero(0.0);
+   ConstantCoefficient tichonov_cf(tichonov);
+   ConstantCoefficient neg_tichonov_cf(-1.0*tichonov);
+   ZCoefficient Z(sdim, psi_gf, alpha);
+   DZCoefficient DZ(sdim, psi_gf, alpha);
+   ScalarVectorProductCoefficient neg_Z(-1.0, Z);
+   DivergenceGridFunctionCoefficient div_psi_cf(&psi_gf);
+   DivergenceGridFunctionCoefficient div_psi_old_cf(&psi_old_gf);
+   SumCoefficient psi_old_minus_psi(div_psi_old_cf, div_psi_cf, 1.0, -1.0);
+
+   // 10. Assemble constant matrices/vectors to avoid reassembly in the loop.
+   LinearForm b0, b1;
+   b0.MakeRef(&RTfes,rhs.GetBlock(0),0);
+   b1.MakeRef(&L2fes,rhs.GetBlock(1),0);
+
+   b0.AddDomainIntegrator(new VectorFEDomainLFIntegrator(neg_Z));
+   b1.AddDomainIntegrator(new DomainLFIntegrator(neg_one));
+   b1.AddDomainIntegrator(new DomainLFIntegrator(psi_old_minus_psi));
+
+   BilinearForm a00(&RTfes);
+   a00.AddDomainIntegrator(new VectorFEMassIntegrator(DZ));
+   a00.AddDomainIntegrator(new VectorFEMassIntegrator(tichonov_cf));
+
+   MixedBilinearForm a10(&RTfes,&L2fes);
+   a10.AddDomainIntegrator(new VectorFEDivergenceIntegrator());
+   a10.Assemble();
+   a10.Finalize();
+   SparseMatrix &A10 = a10.SpMat();
+   SparseMatrix *A01 = Transpose(A10);
+
+   BilinearForm a11(&L2fes);
+   a11.AddDomainIntegrator(new MassIntegrator(neg_tichonov_cf));
+   a11.Assemble();
+   a11.Finalize();
+   SparseMatrix &A11 = a11.SpMat();
+
+   // 11. Iterate.
+   int k;
+   int total_iterations = 0;
+   real_t increment_u = 0.1;
+   GridFunction u_tmp(&L2fes);
+   for (k = 0; k < max_it; k++)
+   {
+      u_tmp = u_old_gf;
+      Z.SetAlpha(alpha);
+      DZ.SetAlpha(alpha);
+
+      mfem::out << "\nOUTER ITERATION " << k+1 << endl;
+
+      int j;
+      for ( j = 0; j < 5; j++)
+      {
+         total_iterations++;
+
+         b0.Assemble();
+         b1.Assemble();
+
+         a00.Assemble(false);
+         a00.Finalize(false);
+         SparseMatrix &A00 = a00.SpMat();
+
+         // Construct Schur-complement preconditioner
+         Vector A00_diag(a00.Height());
+         A00.GetDiag(A00_diag);
+         A00_diag.Reciprocal();
+         SparseMatrix *S = Mult_AtDA(*A01, A00_diag);
+
+         BlockDiagonalPreconditioner prec(offsets);
+         prec.SetDiagonalBlock(0,new DSmoother(A00));
+#ifndef MFEM_USE_SUITESPARSE
+         prec.SetDiagonalBlock(1,new GSSmoother(*S));
+#else
+         prec.SetDiagonalBlock(1,new UMFPackSolver(*S));
+#endif
+         prec.owns_blocks = 1;
+
+         BlockOperator A(offsets);
+         A.SetBlock(0,0,&A00);
+         A.SetBlock(1,0,&A10);
+         A.SetBlock(0,1,A01);
+         A.SetBlock(1,1,&A11);
+
+         GMRES(A,prec,rhs,x,0,2000,500,1e-12,0.0);
+         delete S;
+
+         u_tmp -= u_gf;
+         real_t Newton_update_size = u_tmp.ComputeL2Error(zero);
+         u_tmp = u_gf;
+
+         // Damped Newton update
+         psi_gf.Add(newton_scaling, delta_psi_gf);
+         a00.Update();
+
+         if (visualization)
+         {
+            sol_sock << "solution\n" << mesh << u_gf << "window_title 'Discrete solution'"
+                     << flush;
+         }
+
+         mfem::out << "Newton_update_size = " << Newton_update_size << endl;
+
+         if (Newton_update_size < increment_u)
+         {
+            break;
+         }
+      }
+
+      u_tmp = u_gf;
+      u_tmp -= u_old_gf;
+      increment_u = u_tmp.ComputeL2Error(zero);
+
+      mfem::out << "Number of Newton iterations = " << j+1 << endl;
+      mfem::out << "Increment (|| uₕ - uₕ_prvs||) = " << increment_u << endl;
+
+      u_old_gf = u_gf;
+      psi_old_gf = psi_gf;
+
+      if (increment_u < tol || k == max_it-1)
+      {
+         break;
+      }
+
+      alpha *= max(growth_rate, 1_r);
+
+   }
+
+   mfem::out << "\n Outer iterations: " << k+1
+             << "\n Total iterations: " << total_iterations
+             << "\n Total dofs:       " << RTfes.GetTrueVSize() + L2fes.GetTrueVSize()
+             << endl;
+
+   delete A01;
+   return 0;
+}
+
+void ZCoefficient::Eval(Vector &V, ElementTransformation &T,
+                        const IntegrationPoint &ip)
+{
+   MFEM_ASSERT(psi != NULL, "grid function is not set");
+   MFEM_ASSERT(alpha > 0, "alpha is not positive");
+
+   Vector psi_vals(vdim);
+   psi->GetVectorValue(T, ip, psi_vals);
+   real_t norm = psi_vals.Norml2();
+   real_t phi = 1.0 / sqrt(1.0/(alpha*alpha) + norm*norm);
+
+   V = psi_vals;
+   V *= phi;
+}
+
+void DZCoefficient::Eval(DenseMatrix &K, ElementTransformation &T,
+                         const IntegrationPoint &ip)
+{
+   MFEM_ASSERT(psi != NULL, "grid function is not set");
+   MFEM_ASSERT(alpha > 0, "alpha is not positive");
+
+   Vector psi_vals(height);
+   psi->GetVectorValue(T, ip, psi_vals);
+   real_t norm = psi_vals.Norml2();
+   real_t phi = 1.0 / sqrt(1.0/(alpha*alpha) + norm*norm);
+
+   K = 0.0;
+   for (int i = 0; i < height; i++)
+   {
+      K(i,i) = phi;
+      for (int j = 0; j < height; j++)
+      {
+         K(i,j) -= psi_vals(i) * psi_vals(j) * pow(phi, 3);
+      }
+   }
+}
@@ -0,0 +1,436 @@
+//                       MFEM Example 40 - Parallel Version
+//
+// Compile with: make ex40p
+//
+// Sample runs: mpirun -np 4 ex40p -step 10 -gr 2.0
+//              mpirun -np 4 ex40p -step 10 -gr 2.0 -o 3 -r 1
+//              mpirun -np 4 ex40p -step 10 -gr 2.0 -r 4 -m ../data/l-shape.mesh
+//              mpirun -np 4 ex40p -step 10 -gr 2.0 -r 2 -m ../data/fichera.mesh
+//
+// Description: This example code demonstrates how to use MFEM to solve the
+//              eikonal equation,
+//
+//                      |∇𝑢| = 1 in Ω,  𝑢 = g on ∂Ω.
+//
+//              The solution of this problem coincides with the unique optimum of
+//              the nonlinear program
+//
+//                   maximize ∫_Ω 𝑢 d𝑥 subject to |∇𝑢| ≤ 1,   𝑢 = g on Ω,      (⋆)
+//
+//              which is the foundation for method implemented below.
+//
+//              Following the proximal Galerkin methodology [1] (see also Example
+//              36), we construct a Legendre function for the unit ball
+//              𝐵₁ := {𝑥 ∈ Rⁿ | |𝑥| < 1}. Our choice is the Hellinger entropy,
+//
+//                    h(𝑥) = −( 1 − |𝑥|² )^{1/2},
+//
+//              although other choices are possible, each leading to a slightly
+//              different algorithm. We then adaptively regularize the optimization
+//              problem (⋆) with the Bregman divergence of the Hellinger entropy,
+//
+//                 maximize  ∫_Ω 𝑢 d𝑥 - αₖ⁻¹ Dₕ(∇𝑢,∇𝑢ₖ₋₁)  subject to  𝑢 = g on Ω.
+//
+//              This results in a sequence of functions ( 𝜓ₖ , 𝑢ₖ ),
+//
+//                      𝑢ₖ → 𝑢,    𝜓ₖ/|𝜓ₖ| → ∇𝑢    as k → \infty,
+//
+//              defined by the nonlinear saddle-point problems
+//
+//               Find 𝜓ₖ ∈ H(div,Ω) and 𝑢ₖ ∈ L²(Ω) such that
+//               ( Zₖ(𝜓ₖ) , τ ) + ( 𝑢ₖ , ∇⋅τ ) = ⟨ g , τ⋅n ⟩         ∀ τ ∈ H(div,Ω)
+//               ( ∇⋅𝜓ₖ , v )                 = ( ∇⋅𝜓ₖ₋₁ - 1 , v )   ∀ v ∈ L²(Ω)
+//
+//              where Zₖ(𝜓) := ∇h⁻¹(αₖ 𝜓) = 𝜓 / ( αₖ⁻² + |𝜓|² )^{1/2} and step size
+//              αₖ > 0. These saddle-point problems are solved using a damped Newton's
+//              method. This example assumes that g = 0 and allows the step size to
+//              grow geometrically, αₖ = α₀rᵏ, where r ≥ 1 is the growth rate.
+//
+//              [1] Keith, B. and Surowiec, T. (2023) Proximal Galerkin: A structure-
+//                  preserving finite element method for pointwise bound constraints.
+//                  arXiv:2307.12444 [math.NA]
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+class ZCoefficient : public VectorCoefficient
+{
+protected:
+   ParGridFunction *psi;
+   real_t alpha;
+
+public:
+   ZCoefficient(int vdim, ParGridFunction &psi_, real_t alpha_ = 1.0)
+      : VectorCoefficient(vdim), psi(&psi_), alpha(alpha_) { }
+
+   virtual void Eval(Vector &V, ElementTransformation &T,
+                     const IntegrationPoint &ip);
+   void SetAlpha(real_t alpha_) { alpha = alpha_; }
+};
+
+class DZCoefficient : public MatrixCoefficient
+{
+protected:
+   ParGridFunction *psi;
+   real_t alpha;
+
+public:
+   DZCoefficient(int height, ParGridFunction &psi_, real_t alpha_ = 1.0)
+      : MatrixCoefficient(height),  psi(&psi_), alpha(alpha_) { }
+
+   virtual void Eval(DenseMatrix &K, ElementTransformation &T,
+                     const IntegrationPoint &ip);
+   void SetAlpha(real_t alpha_) { alpha = alpha_; }
+};
+
+int main(int argc, char *argv[])
+{
+   // 0. Initialize MPI and HYPRE.
+   Mpi::Init();
+   int num_procs = Mpi::WorldSize();
+   int myid = Mpi::WorldRank();
+   Hypre::Init();
+
+   // 1. Parse command-line options.
+   const char *mesh_file = "../data/star.mesh";
+   int order = 1;
+   int max_it = 5;
+   int ref_levels = 3;
+   real_t alpha = 1.0;
+   real_t growth_rate = 1.0;
+   real_t newton_scaling = 0.9;
+   real_t tichonov = 1e-1;
+   real_t tol = 1e-4;
+   bool visualization = true;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&ref_levels, "-r", "--refs",
+                  "Number of h-refinements.");
+   args.AddOption(&max_it, "-mi", "--max-it",
+                  "Maximum number of iterations");
+   args.AddOption(&tol, "-tol", "--tol",
+                  "Stopping criteria based on the difference between"
+                  "successive solution updates");
+   args.AddOption(&alpha, "-step", "--step",
+                  "Initial size alpha");
+   args.AddOption(&growth_rate, "-gr", "--growth-rate",
+                  "Growth rate of the step size alpha");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   // 2. Read the mesh from the mesh file.
+   Mesh mesh(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+   int sdim = mesh.SpaceDimension();
+
+   MFEM_ASSERT(mesh.bdr_attributes.Size(),
+               "This example does not currently support meshes"
+               " without boundary attributes."
+              )
+
+   // 3. Postprocess the mesh.
+   // 3A. Refine the mesh to increase the resolution.
+   for (int l = 0; l < ref_levels; l++)
+   {
+      mesh.UniformRefinement();
+   }
+
+   // 3B. Interpolate the geometry after refinement to control geometry error.
+   // NOTE: Minimum second-order interpolation is used to improve the accuracy.
+   int curvature_order = max(order,2);
+   mesh.SetCurvature(curvature_order);
+
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
+
+   // 4. Define the necessary finite element spaces on the mesh.
+   RT_FECollection RTfec(order, dim);
+   ParFiniteElementSpace RTfes(&pmesh, &RTfec);
+
+   L2_FECollection L2fec(order, dim);
+   ParFiniteElementSpace L2fes(&pmesh, &L2fec);
+
+   int num_dofs_RT = RTfes.GlobalTrueVSize();
+   int num_dofs_L2 = L2fes.GlobalTrueVSize();
+   if (myid == 0)
+   {
+      cout << "Number of H(div) dofs: "
+           << num_dofs_RT << endl;
+      cout << "Number of L² dofs: "
+           << num_dofs_L2 << endl;
+   }
+
+   // 5. Define the offsets for the block matrices
+   Array<int> offsets(3);
+   offsets[0] = 0;
+   offsets[1] = RTfes.GetVSize();
+   offsets[2] = L2fes.GetVSize();
+   offsets.PartialSum();
+
+   Array<int> toffsets(3);
+   toffsets[0] = 0;
+   toffsets[1] = RTfes.GetTrueVSize();
+   toffsets[2] = L2fes.GetTrueVSize();
+   toffsets.PartialSum();
+
+   BlockVector x(offsets), rhs(offsets);
+   x = 0.0; rhs = 0.0;
+
+   BlockVector tx(toffsets), trhs(toffsets);
+   tx = 0.0; trhs = 0.0;
+
+   // 6. Define the solution vectors as a finite element grid functions
+   //    corresponding to the fespaces.
+   ParGridFunction u_gf, delta_psi_gf;
+   delta_psi_gf.MakeRef(&RTfes,x,offsets[0]);
+   u_gf.MakeRef(&L2fes,x,offsets[1]);
+
+   ParGridFunction psi_old_gf(&RTfes);
+   ParGridFunction psi_gf(&RTfes);
+   ParGridFunction u_old_gf(&L2fes);
+
+   // 7. Define initial guesses for the solution variables.
+   delta_psi_gf = 0.0;
+   psi_gf = 0.0;
+   u_gf = 0.0;
+   psi_old_gf = psi_gf;
+   u_old_gf = u_gf;
+
+   // 8. Prepare for glvis output.
+   char vishost[] = "localhost";
+   int  visport   = 19916;
+   socketstream sol_sock;
+   if (visualization)
+   {
+      sol_sock.open(vishost,visport);
+      sol_sock.precision(8);
+   }
+
+   // 9. Coefficients to be used later.
+   ConstantCoefficient neg_one(-1.0);
+   ConstantCoefficient zero(0.0);
+   ConstantCoefficient tichonov_cf(tichonov);
+   ConstantCoefficient neg_tichonov_cf(-1.0*tichonov);
+   ZCoefficient Z(sdim, psi_gf, alpha);
+   DZCoefficient DZ(sdim, psi_gf, alpha);
+   ScalarVectorProductCoefficient neg_Z(-1.0, Z);
+   DivergenceGridFunctionCoefficient div_psi_cf(&psi_gf);
+   DivergenceGridFunctionCoefficient div_psi_old_cf(&psi_old_gf);
+   SumCoefficient psi_old_minus_psi(div_psi_old_cf, div_psi_cf, 1.0, -1.0);
+
+   // 10. Assemble constant matrices/vectors to avoid reassembly in the loop.
+   ParLinearForm b0, b1;
+   b0.MakeRef(&RTfes,rhs.GetBlock(0),0);
+   b1.MakeRef(&L2fes,rhs.GetBlock(1),0);
+
+   b0.AddDomainIntegrator(new VectorFEDomainLFIntegrator(neg_Z));
+   b1.AddDomainIntegrator(new DomainLFIntegrator(neg_one));
+   b1.AddDomainIntegrator(new DomainLFIntegrator(psi_old_minus_psi));
+
+   ParBilinearForm a00(&RTfes);
+   a00.AddDomainIntegrator(new VectorFEMassIntegrator(DZ));
+   a00.AddDomainIntegrator(new VectorFEMassIntegrator(tichonov_cf));
+
+   ParMixedBilinearForm a10(&RTfes,&L2fes);
+   a10.AddDomainIntegrator(new VectorFEDivergenceIntegrator());
+   a10.Assemble();
+   a10.Finalize();
+   HypreParMatrix *A10 = a10.ParallelAssemble();
+
+   HypreParMatrix *A01 = A10->Transpose();
+
+   ParBilinearForm a11(&L2fes);
+   a11.AddDomainIntegrator(new MassIntegrator(neg_tichonov_cf));
+   a11.Assemble();
+   a11.Finalize();
+   HypreParMatrix *A11 = a11.ParallelAssemble();
+
+   // 11. Iterate.
+   int k;
+   int total_iterations = 0;
+   real_t increment_u = 0.1;
+   ParGridFunction u_tmp(&L2fes);
+   for (k = 0; k < max_it; k++)
+   {
+      u_tmp = u_old_gf;
+      Z.SetAlpha(alpha);
+      DZ.SetAlpha(alpha);
+
+      if (myid == 0)
+      {
+         mfem::out << "\nOUTER ITERATION " << k+1 << endl;
+      }
+
+      int j;
+      for ( j = 0; j < 5; j++)
+      {
+         total_iterations++;
+
+         b0.Assemble();
+         b0.ParallelAssemble(trhs.GetBlock(0));
+
+         b1.Assemble();
+         b1.ParallelAssemble(trhs.GetBlock(1));
+
+         a00.Assemble(false);
+         a00.Finalize(false);
+         HypreParMatrix *A00 = a00.ParallelAssemble();
+
+         // Construct Schur-complement preconditioner
+         HypreParVector A00_diag(MPI_COMM_WORLD, A00->GetGlobalNumRows(),
+                                 A00->GetRowStarts());
+         A00->GetDiag(A00_diag);
+         HypreParMatrix S_tmp(*A01);
+         S_tmp.InvScaleRows(A00_diag);
+         HypreParMatrix *S = ParMult(A10, &S_tmp, true);
+
+         BlockDiagonalPreconditioner prec(toffsets);
+         HypreBoomerAMG P00(*A00);
+         P00.SetPrintLevel(0);
+         HypreBoomerAMG P11(*S);
+         P11.SetPrintLevel(0);
+         prec.SetDiagonalBlock(0,&P00);
+         prec.SetDiagonalBlock(1,&P11);
+
+         BlockOperator A(toffsets);
+         A.SetBlock(0,0,A00);
+         A.SetBlock(1,0,A10);
+         A.SetBlock(0,1,A01);
+         A.SetBlock(1,1,A11);
+
+         GMRESSolver gmres(MPI_COMM_WORLD);
+         gmres.SetPrintLevel(-1);
+         gmres.SetRelTol(1e-8);
+         gmres.SetMaxIter(2000);
+         gmres.SetKDim(500);
+         gmres.SetOperator(A);
+         gmres.SetPreconditioner(prec);
+         gmres.Mult(trhs,tx);
+         delete S;
+         delete A00;
+
+         delta_psi_gf.SetFromTrueDofs(tx.GetBlock(0));
+         u_gf.SetFromTrueDofs(tx.GetBlock(1));
+
+         u_tmp -= u_gf;
+         real_t Newton_update_size = u_tmp.ComputeL2Error(zero);
+         u_tmp = u_gf;
+
+         // Damped Newton update
+         psi_gf.Add(newton_scaling, delta_psi_gf);
+         a00.Update();
+
+         if (visualization)
+         {
+            sol_sock << "parallel " << num_procs << " " << myid << "\n";
+            sol_sock << "solution\n" << pmesh << u_gf << "window_title 'Discrete solution'"
+                     << flush;
+         }
+
+         if (myid == 0)
+         {
+            mfem::out << "Newton_update_size = " << Newton_update_size << endl;
+         }
+
+         if (Newton_update_size < increment_u)
+         {
+            break;
+         }
+      }
+
+      u_tmp = u_gf;
+      u_tmp -= u_old_gf;
+      increment_u = u_tmp.ComputeL2Error(zero);
+
+      if (myid == 0)
+      {
+         mfem::out << "Number of Newton iterations = " << j+1 << endl;
+         mfem::out << "Increment (|| uₕ - uₕ_prvs||) = " << increment_u << endl;
+      }
+
+      u_old_gf = u_gf;
+      psi_old_gf = psi_gf;
+
+      if (increment_u < tol || k == max_it-1)
+      {
+         break;
+      }
+
+      alpha *= max(growth_rate, 1_r);
+
+   }
+
+   // 12. Print stats.
+   if (myid == 0)
+   {
+      mfem::out << "\n Outer iterations: " << k+1
+                << "\n Total iterations: " << total_iterations
+                << "\n Total dofs:       " << RTfes.GetTrueVSize() + L2fes.GetTrueVSize()
+                << endl;
+   }
+
+   // 13. Free the used memory.
+   delete A01;
+   delete A10;
+   delete A11;
+   return 0;
+}
+
+void ZCoefficient::Eval(Vector &V, ElementTransformation &T,
+                        const IntegrationPoint &ip)
+{
+   MFEM_ASSERT(psi != NULL, "grid function is not set");
+   MFEM_ASSERT(alpha > 0, "alpha is not positive");
+
+   Vector psi_vals(vdim);
+   psi->GetVectorValue(T, ip, psi_vals);
+   real_t norm = psi_vals.Norml2();
+   real_t phi = 1.0 / sqrt(1.0/(alpha*alpha) + norm*norm);
+
+   V = psi_vals;
+   V *= phi;
+}
+
+void DZCoefficient::Eval(DenseMatrix &K, ElementTransformation &T,
+                         const IntegrationPoint &ip)
+{
+   MFEM_ASSERT(psi != NULL, "grid function is not set");
+   MFEM_ASSERT(alpha > 0, "alpha is not positive");
+
+   Vector psi_vals(height);
+   psi->GetVectorValue(T, ip, psi_vals);
+   real_t norm = psi_vals.Norml2();
+   real_t phi = 1.0 / sqrt(1.0/(alpha*alpha) + norm*norm);
+
+   K = 0.0;
+   for (int i = 0; i < height; i++)
+   {
+      K(i,i) = phi;
+      for (int j = 0; j < height; j++)
+      {
+         K(i,j) -= psi_vals(i) * psi_vals(j) * pow(phi, 3);
+      }
+   }
+}
@@ -23,14 +23,14 @@ MFEM_LIB_FILE = mfem_is_not_built

 SEQ_EXAMPLES = ex0 ex1 ex2 ex3 ex4 ex5 ex6 ex7 ex8 ex9 ex10 ex14 ex15 ex16 \
   ex17 ex18 ex19 ex20 ex21 ex22 ex23 ex24 ex25 ex26 ex27 ex28 ex29 ex30 \
-   ex31 ex33 ex34 ex36 ex37 ex38 ex39
+   ex31 ex33 ex34 ex36 ex37 ex38 ex39 ex40
 PAR_EXAMPLES = ex0p ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex8p ex9p ex10p ex11p \
   ex12p ex13p ex14p ex15p ex16p ex17p ex18p ex19p ex20p ex21p ex22p ex24p \
   ex25p ex26p ex27p ex28p ex29p ex30p ex31p ex32p ex33p ex34p ex35p ex36p \
-   ex37p ex39p
-SEQ_DEVICE_EXAMPLES = ex1 ex3 ex4 ex5 ex6 ex9 ex22 ex24 ex25 ex26 ex34
-PAR_DEVICE_EXAMPLES = ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex9p ex13p ex22p \
-   ex24p ex25p ex26p ex34p ex35p
+   ex37p ex39p ex40p
+SEQ_DEVICE_EXAMPLES = ex1 ex3 ex4 ex5 ex6 ex9 ex14 ex22 ex24 ex25 ex26 ex34
+PAR_DEVICE_EXAMPLES = ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex9p ex13p ex14p \
+   ex22p ex24p ex25p ex26p ex34p ex35p

 ifeq ($(MFEM_USE_LAPACK),YES)
   SEQ_EXAMPLES += ex38
@@ -138,6 +138,14 @@ ex10-test-seq: ex10
 	@$(call mfem-test,$<,, Serial example,-tf 5)
 ex10p-test-par: ex10p
 	@$(call mfem-test,$<, $(RUN_MPI), Parallel example,-tf 5)
+ex14-test-seq-cuda: ex14
+	@$(call mfem-test,$<,, Serial CUDA example,-r 2 -pa -d cuda)
+ex14p-test-par-cuda: ex14p
+	@$(call mfem-test,$<, $(RUN_MPI), Parallel CUDA example,-rs 2 -rp 0 -pa -d cuda)
+ex14-test-seq-hip: ex14
+	@$(call mfem-test,$<,, Serial HIP example,-r 2 -pa -d hip)
+ex14p-test-par-hip: ex14p
+	@$(call mfem-test,$<, $(RUN_MPI), Parallel HIP example,-rs 2 -rp 0 -pa -d hip)
 ex15-test-seq: ex15
 	@$(call mfem-test,$<,, Serial example,-e 1)
 ex15p-test-par: ex15p
@@ -1,3 +1,14 @@
+// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
 #include <algorithm>
 #include <assert.h>
 #include <cstdlib>
@@ -709,10 +709,7 @@ real_t HyperelasticOperator::ElasticEnergy(const ParGridFunction &x) const

 real_t HyperelasticOperator::KineticEnergy(const ParGridFunction &v) const
 {
-   real_t loc_energy = 0.5*M.InnerProduct(v, v);
-   real_t energy;
-   MPI_Allreduce(&loc_energy, &energy, 1, MPITypeMap<real_t>::mpi_type, MPI_SUM,
-                 fespace.GetComm());
+   real_t energy = 0.5*M.ParInnerProduct(v, v);
   return energy;
 }

@@ -856,10 +856,7 @@ double HyperelasticOperator::ElasticEnergy(const ParGridFunction &x) const

 double HyperelasticOperator::KineticEnergy(const ParGridFunction &v) const
 {
-   double loc_energy = 0.5*M.InnerProduct(v, v);
-   double energy;
-   MPI_Allreduce(&loc_energy, &energy, 1, MPI_DOUBLE, MPI_SUM,
-                 fespace.GetComm());
+   double energy = 0.5*M.ParInnerProduct(v, v);
   return energy;
 }

@@ -340,9 +340,9 @@ public:
       $ M^{-1} $ (currently returns NULL) */
   virtual MatrixInverse *Inverse() const;

-   /** @brief  Finalizes the matrix initialization if the ::AssemblyLevel is
+   /** @brief Finalizes the matrix initialization if the ::AssemblyLevel is
       AssemblyLevel::LEGACY.
-       THe matrix that gets finalized is different if you are using static
+       The matrix that gets finalized is different if you are using static
       condensation or hybridization.*/
   virtual void Finalize(int skip_zeros = 1);

@@ -643,7 +643,7 @@ public:
   void EliminateVDofs(const Array<int> &vdofs, const Vector &sol, Vector &rhs,
                       DiagonalPolicy dpolicy = DIAG_ONE);

-   /** @brief  Eliminate the given @a vdofs, storing the eliminated part
+   /** @brief Eliminate the given @a vdofs, storing the eliminated part
       internally in $ M_e $.

       This method works in conjunction with EliminateVDofsInRHS() and allows
@@ -826,7 +826,7 @@ public:
       $ M^{-1} $ (currently unimplemented and returns NULL)*/
   virtual MatrixInverse *Inverse() const;

-   /** @brief  Finalizes the matrix initialization if the ::AssemblyLevel is
+   /** @brief Finalizes the matrix initialization if the ::AssemblyLevel is
       AssemblyLevel::LEGACY.*/
   virtual void Finalize(int skip_zeros = 1);

@@ -589,10 +589,27 @@ void PABilinearFormExtension::Mult(const Vector &x, Vector &y) const
   const int iFISz = intFaceIntegrators.Size();
   if (int_face_restrict_lex && iFISz>0)
   {
-      int_face_restrict_lex->Mult(x, int_face_X);
+      // When assembling interior face integrators for DG spaces, we need to
+      // exchange the face-neighbor information. This happens inside member
+      // functions of the 'int_face_restrict_lex'. To avoid repeated calls to
+      // ParGridFunction::ExchangeFaceNbrData, if we have a parallel space
+      // with interior face integrators, we create a ParGridFunction that
+      // will be used to cache the face-neighbor data. x_dg should be passed
+      // to any restriction operator that may need to use face-neighbor data.
+      const Vector *x_dg = &x;
+#ifdef MFEM_USE_MPI
+      ParGridFunction x_pgf;
+      if (auto *pfes = dynamic_cast<ParFiniteElementSpace*>(a->FESpace()))
+      {
+         x_pgf.MakeRef(pfes, const_cast<Vector&>(x), 0);
+         x_dg = &x_pgf;
+      }
+#endif
+
+      int_face_restrict_lex->Mult(*x_dg, int_face_X);
      if (int_face_dXdn.Size() > 0)
      {
-         int_face_restrict_lex->NormalDerivativeMult(x, int_face_dXdn);
+         int_face_restrict_lex->NormalDerivativeMult(*x_dg, int_face_dXdn);
      }
      if (int_face_X.Size() > 0)
      {
@@ -1741,7 +1741,7 @@ public:
   { vector_fe.CalcPhysDShape(Trans, shape); }
 };

-/** Class for integrating the bilinear form $a(u,v) := (-\hat{V} \cdot  \nabla \cdot u, \nabla \cdot v)$ in 2D
+/** Class for integrating the bilinear form $a(u,v) := (-\hat{V} \cdot \nabla u, \nabla \cdot v)$ in 2D
    or 3D and where $\hat{V}$ is a vector coefficient, $u$ is in $H^1$ and $v$ is in $H(div)$. */
 class MixedGradDivIntegrator : public MixedScalarVectorIntegrator
 {
@@ -1780,7 +1780,7 @@ public:
   { scalar_fe.CalcPhysDivShape(Trans, shape); }
 };

-/** Class for integrating the bilinear form $a(u,v) := (-\hat{V} \nabla \cdot u, \nabla \cdot v)$ in 2D
+/** Class for integrating the bilinear form $a(u,v) := (-\hat{V} \nabla \cdot u, \nabla v)$ in 2D
    or 3D and where $\hat{V}$ is a vector coefficient, $u$ is in $H(div)$ and $v$ is in $H^1$. */
 class MixedDivGradIntegrator : public MixedScalarVectorIntegrator
 {
@@ -1820,7 +1820,7 @@ public:
   { scalar_fe.CalcPhysDivShape(Trans, shape); }
 };

-/** Class for integrating the bilinear form $a(u,v) := (-\hat{V} u, \nabla \cdot v)$ in 2D or 3D
+/** Class for integrating the bilinear form $a(u,v) := (-\hat{V} u, \nabla v)$ in 2D or 3D
    and where $\hat{V}$ is a vector coefficient, $u$ is in $H^1$ or $L_2$ and $v$ is in $H^1$. */
 class MixedScalarWeakDivergenceIntegrator : public MixedScalarVectorIntegrator
 {
@@ -807,6 +807,7 @@ void SymmetricMatrixCoefficient::ProjectSymmetric(QuadratureFunction &qf)

   QuadratureSpaceBase &qspace = *qf.GetSpace();
   const int ne = qspace.GetNE();
+   qf.HostWrite();
   DenseMatrix values;
   DenseSymmetricMatrix matrix;
   for (int iel = 0; iel < ne; ++iel)
@@ -818,7 +819,7 @@ void SymmetricMatrixCoefficient::ProjectSymmetric(QuadratureFunction &qf)
      {
         const IntegrationPoint &ip = ir[iq];
         T.SetIntPoint(&ip);
-         matrix.UseExternalData(&values(0, iq), vdim);
+         matrix.UseExternalData(&values(0, iq), height);
         Eval(matrix, T, ip);
      }
   }
@@ -828,13 +829,12 @@ void SymmetricMatrixCoefficient::ProjectSymmetric(QuadratureFunction &qf)
 void SymmetricMatrixCoefficient::Eval(DenseMatrix &K, ElementTransformation &T,
                                      const IntegrationPoint &ip)
 {
-   mat.SetSize(height);
-   Eval(mat, T, ip);
+   Eval(mat_aux, T, ip);
   for (int j = 0; j < width; ++j)
   {
      for (int i = 0; i < height; ++ i)
      {
-         K(i, j) = mat(i, j);
+         K(i, j) = mat_aux(i, j);
      }
   }
 }
@@ -924,6 +924,75 @@ void MatrixArrayCoefficient::Eval(DenseMatrix &K, ElementTransformation &T,
   }
 }

+MatrixArrayVectorCoefficient::MatrixArrayVectorCoefficient (int dim)
+   : MatrixCoefficient (dim)
+{
+   Coeff.SetSize(height);
+   ownCoeff.SetSize(height);
+   for (int i = 0; i < height; i++)
+   {
+      Coeff[i] = NULL;
+      ownCoeff[i] = true;
+   }
+}
+
+void MatrixArrayVectorCoefficient::SetTime(real_t t)
+{
+   for (int i=0; i < height; i++)
+   {
+      if (Coeff[i]) { Coeff[i]->SetTime(t); }
+   }
+   this->MatrixCoefficient::SetTime(t);
+}
+
+void MatrixArrayVectorCoefficient::Set(int i, VectorCoefficient * c, bool own)
+{
+   MFEM_ASSERT(i < height && i >= 0, "Row "
+               << i << " does not exist. " <<
+               "Matrix height = " << height << ".");
+   if (ownCoeff[i]) { delete Coeff[i]; }
+   Coeff[i] = c;
+   ownCoeff[i] = own;
+}
+
+MatrixArrayVectorCoefficient::~MatrixArrayVectorCoefficient ()
+{
+   for (int i=0; i < height; i++)
+   {
+      if (ownCoeff[i]) { delete Coeff[i]; }
+   }
+}
+
+void MatrixArrayVectorCoefficient::Eval(int i, Vector &V,
+                                        ElementTransformation &T,
+                                        const IntegrationPoint &ip)
+{
+   MFEM_ASSERT(i < height && i >= 0, "Row "
+               << i << " does not exist. " <<
+               "Matrix height = " << height << ".");
+   if (Coeff[i])
+   {
+      Coeff[i] -> Eval(V, T, ip);
+   }
+   else
+   {
+      V = 0.0;
+   }
+}
+
+void MatrixArrayVectorCoefficient::Eval(DenseMatrix &K,
+                                        ElementTransformation &T,
+                                        const IntegrationPoint &ip)
+{
+   K.SetSize(height, width);
+   Vector V(width);
+   for (int i = 0; i < height; i++)
+   {
+      this->Eval(i, V, T, ip);
+      K.SetRow(i, V);
+   }
+}
+
 void MatrixRestrictedCoefficient::SetTime(real_t t)
 {
   if (c) { c->SetTime(t); }
@@ -1041,6 +1110,27 @@ real_t DeterminantCoefficient::Eval(ElementTransformation &T,
   return ma.Det();
 }

+TraceCoefficient::TraceCoefficient(MatrixCoefficient &A)
+   : a(&A), ma(A.GetHeight(), A.GetWidth())
+{
+   MFEM_ASSERT(A.GetHeight() == A.GetWidth(),
+               "TraceCoefficient:  "
+               "Argument must be a square matrix.");
+}
+
+void TraceCoefficient::SetTime(real_t t)
+{
+   if (a) { a->SetTime(t); }
+   this->Coefficient::SetTime(t);
+}
+
+real_t TraceCoefficient::Eval(ElementTransformation &T,
+                              const IntegrationPoint &ip)
+{
+   a->Eval(ma, T, ip);
+   return ma.Trace();
+}
+
 VectorSumCoefficient::VectorSumCoefficient(int dim)
   : VectorCoefficient(dim),
     ACoef(NULL), BCoef(NULL),
@@ -1326,6 +1416,30 @@ void InverseMatrixCoefficient::Eval(DenseMatrix &M,
   M.Invert();
 }

+ExponentialMatrixCoefficient::ExponentialMatrixCoefficient(MatrixCoefficient &A)
+   : MatrixCoefficient(A.GetHeight(), A.GetWidth()), a(&A)
+{
+   MFEM_ASSERT(A.GetHeight() == A.GetWidth() && A.GetHeight() == 2,
+               "ExponentialMatrixCoefficient:  "
+               << "Argument must be a square 2x2 matrix."
+               << "  Height = " << A.GetHeight()
+               << ", Width = " << A.GetWidth());
+}
+
+void ExponentialMatrixCoefficient::SetTime(real_t t)
+{
+   if (a) { a->SetTime(t); }
+   this->MatrixCoefficient::SetTime(t);
+}
+
+void ExponentialMatrixCoefficient::Eval(DenseMatrix &M,
+                                        ElementTransformation &T,
+                                        const IntegrationPoint &ip)
+{
+   a->Eval(M, T, ip);
+   M.Exponential();
+}
+
 OuterProductCoefficient::OuterProductCoefficient(VectorCoefficient &A,
                                                 VectorCoefficient &B)
   : MatrixCoefficient(A.GetVDim(), B.GetVDim()), a(&A), b(&B),
@@ -1334,6 +1334,46 @@ public:
   virtual ~MatrixArrayCoefficient();
 };

+/** @brief Matrix coefficient defined row-wise by an array of vector
+    coefficients. Rows that are not set will evaluate to zero. The
+    matrix coefficient is stored as an array indexing the rows of
+    the matrix. */
+class MatrixArrayVectorCoefficient : public MatrixCoefficient
+{
+private:
+   Array<VectorCoefficient *> Coeff;
+   Array<bool> ownCoeff;
+
+public:
+   /** @brief Construct a coefficient matrix of dimensions @a dim * @a dim. The
+       actual coefficients still need to be added with Set(). */
+   explicit MatrixArrayVectorCoefficient (int dim);
+
+   /// Set the time for internally stored coefficients
+   void SetTime(real_t t) override;
+
+   /// Get the vector coefficient located at the i-th row of the matrix
+   VectorCoefficient* GetCoeff (int i) { return Coeff[i]; }
+
+   /** @brief Set the coefficient located at the i-th row of the matrix.
+       By this will take ownership of the Coefficient passed in, but this
+       can be overridden with the @a own parameter. */
+   void Set(int i, VectorCoefficient * c, bool own=true);
+
+   using MatrixCoefficient::Eval;
+
+   /// Evaluate coefficient located at the i-th row of the matrix using integration
+   /// point @a ip.
+   void Eval(int i, Vector &V, ElementTransformation &T,
+             const IntegrationPoint &ip);
+
+   /// Evaluate the matrix coefficient @a ip.
+   void Eval(DenseMatrix &K, ElementTransformation &T,
+             const IntegrationPoint &ip) override;
+
+   virtual ~MatrixArrayVectorCoefficient();
+};
+

 /** @brief Derived matrix coefficient that has the value of the parent matrix
    coefficient where it is active and is zero otherwise. */
@@ -1426,12 +1466,13 @@ public:
 class SymmetricMatrixCoefficient : public MatrixCoefficient
 {
 protected:
+
   /// Internal matrix used when evaluating this coefficient as a DenseMatrix.
-   DenseSymmetricMatrix mat;
+   mutable DenseSymmetricMatrix mat_aux;
 public:
   /// Construct a dim x dim matrix coefficient.
   explicit SymmetricMatrixCoefficient(int dimension)
-      : MatrixCoefficient(dimension, true) { }
+      : MatrixCoefficient(dimension, true), mat_aux(height) { }

   /// Get the size of the matrix.
   int GetSize() const { return height; }
@@ -1464,8 +1505,9 @@ public:
   virtual void Eval(DenseMatrix &K, ElementTransformation &T,
                     const IntegrationPoint &ip);

-   /// Return a reference to the constant matrix.
-   const DenseSymmetricMatrix& GetMatrix() { return mat; }
+
+   /// @deprecated Return a reference to the internal matrix used when evaluating this coefficient as a DenseMatrix.
+   MFEM_DEPRECATED const DenseSymmetricMatrix& GetMatrix() { return mat_aux; }

   virtual ~SymmetricMatrixCoefficient() { }
 };
@@ -1485,6 +1527,10 @@ public:
   /// Evaluate the matrix coefficient at @a ip.
   virtual void Eval(DenseSymmetricMatrix &M, ElementTransformation &T,
                     const IntegrationPoint &ip) { M = mat; }
+
+   /// Return a reference to the constant matrix.
+   const DenseSymmetricMatrix& GetMatrix() { return mat; }
+
 };


@@ -1761,6 +1807,31 @@ public:
                       const IntegrationPoint &ip);
 };

+/// Scalar coefficient defined as the trace of a matrix coefficient
+class TraceCoefficient : public Coefficient
+{
+private:
+   MatrixCoefficient * a;
+
+   mutable DenseMatrix ma;
+
+public:
+   /// Construct with the matrix.
+   TraceCoefficient(MatrixCoefficient &A);
+
+   /// Set the time for internally stored coefficients
+   void SetTime(real_t t);
+
+   /// Reset the matrix coefficient
+   void SetACoef(MatrixCoefficient &A) { a = &A; }
+   /// Return the matrix coefficient
+   MatrixCoefficient * GetACoef() const { return a; }
+
+   /// Evaluate the trace coefficient at @a ip.
+   virtual real_t Eval(ElementTransformation &T,
+                       const IntegrationPoint &ip);
+};
+
 /// Vector coefficient defined as the linear combination of two vectors
 class VectorSumCoefficient : public VectorCoefficient
 {
@@ -2112,7 +2183,7 @@ public:
                     const IntegrationPoint &ip);
 };

-/// Matrix coefficient defined as the transpose a matrix coefficient
+/// Matrix coefficient defined as the transpose of a matrix coefficient
 class TransposeMatrixCoefficient : public MatrixCoefficient
 {
 private:
@@ -2135,7 +2206,7 @@ public:
                     const IntegrationPoint &ip);
 };

-/// Matrix coefficient defined as the inverse a matrix coefficient.
+/// Matrix coefficient defined as the inverse of a matrix coefficient.
 class InverseMatrixCoefficient : public MatrixCoefficient
 {
 private:
@@ -2158,6 +2229,29 @@ public:
                     const IntegrationPoint &ip);
 };

+/// Matrix coefficient defined as the exponential of a matrix coefficient.
+class ExponentialMatrixCoefficient : public MatrixCoefficient
+{
+private:
+   MatrixCoefficient * a;
+
+public:
+   /// Construct the matrix coefficient.  Result is $ \exp(A) $.
+   ExponentialMatrixCoefficient(MatrixCoefficient &A);
+
+   /// Set the time for internally stored coefficients
+   void SetTime(real_t t);
+
+   /// Reset the matrix coefficient
+   void SetACoef(MatrixCoefficient &A) { a = &A; }
+   /// Return the matrix coefficient
+   MatrixCoefficient * GetACoef() const { return a; }
+
+   /// Evaluate the matrix coefficient at @a ip.
+   virtual void Eval(DenseMatrix &M, ElementTransformation &T,
+                     const IntegrationPoint &ip);
+};
+
 /// Matrix coefficient defined as the outer product of two vector coefficients.
 class OuterProductCoefficient : public MatrixCoefficient
 {
@@ -1243,25 +1243,16 @@ ParSesquilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list,
         HypreParMatrix * Ah;
         A_i.Get(Ah);
         hypre_ParCSRMatrix *Aih = *Ah;
-#if !defined(HYPRE_USING_GPU)
-         ess_tdof_list.HostRead();
-         for (int k = 0; k < n; k++)
-         {
-            const int j = ess_tdof_list[k];
-            Aih->diag->data[Aih->diag->i[j]] = 0.0;
-         }
-#else
         Ah->HypreReadWrite();
         const int *d_ess_tdof_list =
-            ess_tdof_list.GetMemory().Read(MemoryClass::DEVICE, n);
-         const int *d_diag_i = Aih->diag->i;
+            ess_tdof_list.GetMemory().Read(GetHypreMemoryClass(), n);
+         HYPRE_Int *d_diag_i = Aih->diag->i;
         real_t *d_diag_data = Aih->diag->data;
-         MFEM_GPU_FORALL(k, n,
+         mfem::hypre_forall(n, [=] MFEM_HOST_DEVICE (int k)
         {
            const int j = d_ess_tdof_list[k];
            d_diag_data[d_diag_i[j]] = 0.0;
         });
-#endif
      }
      else
      {
@@ -1,8 +1,18 @@
+// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
 #include "convergence.hpp"

 using namespace std;

-
 namespace mfem
 {

@@ -101,7 +101,7 @@ void DGMassInverse::SetRelTol(const real_t rel_tol_) { rel_tol = rel_tol_; }

 void DGMassInverse::SetAbsTol(const real_t abs_tol_) { abs_tol = abs_tol_; }

-void DGMassInverse::SetMaxIter(const real_t max_iter_) { max_iter = max_iter_; }
+void DGMassInverse::SetMaxIter(const int max_iter_) { max_iter = max_iter_; }

 void DGMassInverse::Update()
 {
@@ -137,7 +137,7 @@ void DGMassInverse::DGMassCGIteration(const Vector &b_, Vector &u_) const

   const real_t RELTOL = rel_tol;
   const real_t ABSTOL = abs_tol;
-   const real_t MAXIT = max_iter;
+   const int MAXIT = max_iter;
   const bool IT_MODE = iterative_mode;
   const bool CHANGE_BASIS = (d2q != nullptr);

@@ -96,7 +96,7 @@ public:
   /// Set the absolute tolerance.
   void SetAbsTol(const real_t abs_tol_);
   /// Set the maximum number of iterations.
-   void SetMaxIter(const real_t max_iter_);
+   void SetMaxIter(const int max_iter_);
   /// Recompute operator and preconditioner (when coefficient or mesh changes).
   void Update();

@@ -316,7 +316,7 @@ public:
   int GetDim() const { return dim; }

   /** @brief Returns the vector dimension for vector-valued finite elements,
-       which is also the dimension of the interpolation operatrion. */
+       which is also the dimension of the interpolation operation. */
   int GetRangeDim() const { return vdim; }

   /// Returns the dimension of the curl for vector-valued finite elements.
@@ -1321,9 +1321,9 @@ void GridFunction::ProjectVectorFieldOn(GridFunction &vec_field, int comp)
   }
 }

-void GridFunction::AccumulateAndCountDerivativeValues(int comp, int der_comp,
-                                                      GridFunction &der,
-                                                      Array<int> &zones_per_dof)
+void GridFunction::AccumulateAndCountDerivativeValues(
+   int comp, int der_comp, GridFunction &der,
+   Array<int> &zones_per_dof) const
 {
   FiniteElementSpace * der_fes = der.FESpace();
   ElementTransformation * transf;
@@ -1374,7 +1374,8 @@ void GridFunction::AccumulateAndCountDerivativeValues(int comp, int der_comp,
   }
 }

-void GridFunction::GetDerivative(int comp, int der_comp, GridFunction &der)
+void GridFunction::GetDerivative(int comp, int der_comp,
+                                 GridFunction &der) const
 {
   Array<int> overlap;
   AccumulateAndCountDerivativeValues(comp, der_comp, der, overlap);
@@ -2061,41 +2062,37 @@ void GridFunction::AccumulateAndCountBdrValues(
   Coefficient *coeff[], VectorCoefficient *vcoeff, const Array<int> &attr,
   Array<int> &values_counter)
 {
-   int i, j, fdof, d, ind, vdim;
-   real_t val;
-   const FiniteElement *fe;
-   ElementTransformation *transf;
   Array<int> vdofs;
   Vector vc;

   values_counter.SetSize(Size());
   values_counter = 0;

-   vdim = fes->GetVDim();
-
+   const int vdim = fes->GetVDim();
   HostReadWrite();

-   for (i = 0; i < fes->GetNBE(); i++)
+   for (int i = 0; i < fes->GetNBE(); i++)
   {
      if (attr[fes->GetBdrAttribute(i) - 1] == 0) { continue; }

-      fe = fes->GetBE(i);
-      fdof = fe->GetDof();
-      transf = fes->GetBdrElementTransformation(i);
+      const FiniteElement *fe = fes->GetBE(i);
+      const int fdof = fe->GetDof();
+      ElementTransformation *transf = fes->GetBdrElementTransformation(i);
      const IntegrationRule &ir = fe->GetNodes();
      fes->GetBdrElementVDofs(i, vdofs);

-      for (j = 0; j < fdof; j++)
+      for (int j = 0; j < fdof; j++)
      {
         const IntegrationPoint &ip = ir.IntPoint(j);
         transf->SetIntPoint(&ip);
         if (vcoeff) { vcoeff->Eval(vc, *transf, ip); }
-         for (d = 0; d < vdim; d++)
+         for (int d = 0; d < vdim; d++)
         {
            if (!vcoeff && !coeff[d]) { continue; }

-            val = vcoeff ? vc(d) : coeff[d]->Eval(*transf, ip);
-            if ( (ind = vdofs[fdof*d+j]) < 0 )
+            real_t val = vcoeff ? vc(d) : coeff[d]->Eval(*transf, ip);
+            int ind = vdofs[fdof*d+j];
+            if ( ind < 0 )
            {
               val = -val, ind = -1-ind;
            }
@@ -2117,10 +2114,11 @@ void GridFunction::AccumulateAndCountBdrValues(
   // iff A_ij != 0. It is sufficient to resolve just the first level of
   // dependency, since A is a projection matrix: A^n = A due to cR.cP = I.
   // Cases like these arise in 3D when boundary edges are constrained by
-   // (depend on) internal faces/elements. We use the virtual method
-   // GetBoundaryClosure from NCMesh to resolve the dependencies.
-
-   if (fes->Nonconforming() && fes->GetMesh()->Dimension() == 3)
+   // (depend on) internal faces/elements, or for internal boundaries in 2 or
+   // 3D. We use the virtual method GetBoundaryClosure from NCMesh to resolve
+   // the dependencies.
+   if (fes->Nonconforming() && (fes->GetMesh()->Dimension() == 2 ||
+                                fes->GetMesh()->Dimension() == 3))
   {
      Vector vals;
      Mesh *mesh = fes->GetMesh();
@@ -2128,26 +2126,19 @@ void GridFunction::AccumulateAndCountBdrValues(
      Array<int> bdr_edges, bdr_vertices, bdr_faces;
      ncmesh->GetBoundaryClosure(attr, bdr_vertices, bdr_edges, bdr_faces);

-      for (i = 0; i < bdr_edges.Size(); i++)
+      auto mark_dofs = [&](ElementTransformation &transf, const FiniteElement &fe)
      {
-         int edge = bdr_edges[i];
-         fes->GetEdgeVDofs(edge, vdofs);
-         if (vdofs.Size() == 0) { continue; }
-
-         transf = mesh->GetEdgeTransformation(edge);
-         transf->Attribute = -1; // TODO: set the boundary attribute
-         fe = fes->GetEdgeElement(edge);
         if (!vcoeff)
         {
-            vals.SetSize(fe->GetDof());
-            for (d = 0; d < vdim; d++)
+            vals.SetSize(fe.GetDof());
+            for (int d = 0; d < vdim; d++)
            {
               if (!coeff[d]) { continue; }

-               fe->Project(*coeff[d], *transf, vals);
+               fe.Project(*coeff[d], transf, vals);
               for (int k = 0; k < vals.Size(); k++)
               {
-                  ind = vdofs[d*vals.Size()+k];
+                  const int ind = vdofs[d*vals.Size()+k];
                  if (++values_counter[ind] == 1)
                  {
                     (*this)(ind) = vals(k);
@@ -2161,11 +2152,11 @@ void GridFunction::AccumulateAndCountBdrValues(
         }
         else // vcoeff != NULL
         {
-            vals.SetSize(vdim*fe->GetDof());
-            fe->Project(*vcoeff, *transf, vals);
+            vals.SetSize(vdim*fe.GetDof());
+            fe.Project(*vcoeff, transf, vals);
            for (int k = 0; k < vals.Size(); k++)
            {
-               ind = vdofs[k];
+               const int ind = vdofs[k];
               if (++values_counter[ind] == 1)
               {
                  (*this)(ind) = vals(k);
@@ -2176,6 +2167,26 @@ void GridFunction::AccumulateAndCountBdrValues(
               }
            }
         }
+      };
+
+      for (auto edge : bdr_edges)
+      {
+         fes->GetEdgeVDofs(edge, vdofs);
+         if (vdofs.Size() == 0) { continue; }
+
+         ElementTransformation *transf = mesh->GetEdgeTransformation(edge);
+         const FiniteElement *fe = fes->GetEdgeElement(edge);
+         mark_dofs(*transf, *fe);
+      }
+
+      for (auto face : bdr_faces)
+      {
+         fes->GetFaceVDofs(face, vdofs);
+         if (vdofs.Size() == 0) { continue; }
+
+         ElementTransformation *transf = mesh->GetFaceTransformation(face);
+         const FiniteElement *fe = fes->GetFaceElement(face);
+         mark_dofs(*transf, *fe);
      }
   }
 }
@@ -2228,26 +2239,37 @@ void GridFunction::AccumulateAndCountBdrTangentValues(
      accumulate_dofs(dofs, lvec, *this, values_counter);
   }

-   if (fes->Nonconforming() && fes->GetMesh()->Dimension() == 3)
+   if (fes->Nonconforming() && (fes->GetMesh()->Dimension() == 2 ||
+                                fes->GetMesh()->Dimension() == 3))
   {
      Mesh *mesh = fes->GetMesh();
      NCMesh *ncmesh = mesh->ncmesh;
      Array<int> bdr_edges, bdr_vertices, bdr_faces;
      ncmesh->GetBoundaryClosure(bdr_attr, bdr_vertices, bdr_edges, bdr_faces);

-      for (int i = 0; i < bdr_edges.Size(); i++)
+      for (auto edge : bdr_edges)
      {
-         int edge = bdr_edges[i];
         fes->GetEdgeDofs(edge, dofs);
         if (dofs.Size() == 0) { continue; }

         T = mesh->GetEdgeTransformation(edge);
-         T->Attribute = -1; // TODO: set the boundary attribute
         fe = fes->GetEdgeElement(edge);
         lvec.SetSize(fe->GetDof());
         fe->Project(vcoeff, *T, lvec);
         accumulate_dofs(dofs, lvec, *this, values_counter);
      }
+
+      for (auto face : bdr_faces)
+      {
+         fes->GetFaceDofs(face, dofs);
+         if (dofs.Size() == 0) { continue; }
+
+         T = mesh->GetFaceTransformation(face);
+         fe = fes->GetFaceElement(face);
+         lvec.SetSize(fe->GetDof());
+         fe->Project(vcoeff, *T, lvec);
+         accumulate_dofs(dofs, lvec, *this, values_counter);
+      }
   }
 }

@@ -321,7 +321,7 @@ public:
       @param[out] der       The resulting derivative (scalar function). The
                             FiniteElementSpace of this function must be set
                             before the call. */
-   void GetDerivative(int comp, int der_comp, GridFunction &der);
+   void GetDerivative(int comp, int der_comp, GridFunction &der) const;

   real_t GetDivergence(ElementTransformation &tr) const;

@@ -443,7 +443,7 @@ protected:
       GetDerivative() method; see its documentation. */
   void AccumulateAndCountDerivativeValues(int comp, int der_comp,
                                           GridFunction &der,
-                                           Array<int> &zones_per_dof);
+                                           Array<int> &zones_per_dof) const;

   void AccumulateAndCountBdrValues(Coefficient *coeff[],
                                    VectorCoefficient *vcoeff,
@@ -1352,6 +1352,85 @@ void OversetFindPointsGSLIB::Interpolate(const Vector &point_pos,
   Interpolate(field_in, field_out);
 }

+GSOPGSLIB::GSOPGSLIB(Array<long long> &ids)
+{
+   gsl_comm = new gslib::comm;
+   cr       = new gslib::crystal;
+#ifdef MFEM_USE_MPI
+   int initialized;
+   MPI_Initialized(&initialized);
+   if (!initialized) { MPI_Init(NULL, NULL); }
+   MPI_Comm comm = MPI_COMM_WORLD;
+   comm_init(gsl_comm, comm);
+#else
+   comm_init(gsl_comm, 0);
+#endif
+   crystal_init(cr, gsl_comm);
+   UpdateIdentifiers(ids);
+}
+
+#ifdef MFEM_USE_MPI
+GSOPGSLIB::GSOPGSLIB(MPI_Comm comm_, Array<long long> &ids)
+   : cr(NULL), gsl_comm(NULL)
+{
+   gsl_comm = new gslib::comm;
+   cr      = new gslib::crystal;
+   comm_init(gsl_comm, comm_);
+   crystal_init(cr, gsl_comm);
+   UpdateIdentifiers(ids);
+}
+#endif
+
+GSOPGSLIB::~GSOPGSLIB()
+{
+   crystal_free(cr);
+   gslib_gs_free(gsl_data);
+   comm_free(gsl_comm);
+   delete gsl_comm;
+   delete cr;
+}
+
+void GSOPGSLIB::UpdateIdentifiers(const Array<long long> &ids)
+{
+   long long minval = ids.Min();
+#ifdef MFEM_USE_MPI
+   MPI_Allreduce(MPI_IN_PLACE, &minval, 1, MPI_LONG_LONG_INT,
+                 MPI_MIN, gsl_comm->c);
+#endif
+   MFEM_VERIFY(minval >= 0, "Unique identifier cannot be negative.");
+   if (gsl_data != NULL) { gslib_gs_free(gsl_data); }
+   num_ids = ids.Size();
+   gsl_data = gslib_gs_setup(ids.GetData(),
+                             ids.Size(),
+                             gsl_comm, 0,
+                             gslib::gs_crystal_router, 0);
+}
+
+void GSOPGSLIB::GS(Vector &senddata, GSOp op)
+{
+   MFEM_VERIFY(senddata.Size() == num_ids,
+               "Incompatible setup and GOP operation.");
+   if (op == GSOp::ADD)
+   {
+      gslib_gs(senddata.GetData(),gslib::gs_double,gslib::gs_add,0,gsl_data,0);
+   }
+   else if (op == GSOp::MUL)
+   {
+      gslib_gs(senddata.GetData(),gslib::gs_double,gslib::gs_mul,0,gsl_data,0);
+   }
+   else if (op == GSOp::MAX)
+   {
+      gslib_gs(senddata.GetData(),gslib::gs_double,gslib::gs_max,0,gsl_data,0);
+   }
+   else if (op == GSOp::MIN)
+   {
+      gslib_gs(senddata.GetData(),gslib::gs_double,gslib::gs_min,0,gsl_data,0);
+   }
+   else
+   {
+      MFEM_ABORT("Invalid GSOp operation.");
+   }
+}

 } // namespace mfem

@@ -23,13 +23,16 @@ struct comm;
 struct findpts_data_2;
 struct findpts_data_3;
 struct crystal;
+struct gs_data;
 }

 namespace mfem
 {

 /** \brief FindPointsGSLIB can robustly evaluate a GridFunction on an arbitrary
- *  collection of points. There are three key functions in FindPointsGSLIB:
+ *  collection of points.
+ *
+ *  There are three key functions in FindPointsGSLIB:
 *
 *  1. Setup - constructs the internal data structures of gslib.
 *
@@ -226,6 +229,7 @@ public:

 /** \brief OversetFindPointsGSLIB enables use of findpts for arbitrary number of
    overlapping grids.
+
    The parameters in this class are the same as FindPointsGSLIB with the
    difference of additional inputs required to account for more than 1 mesh. */
 class OversetFindPointsGSLIB : public FindPointsGSLIB
@@ -290,6 +294,63 @@ public:
   using FindPointsGSLIB::Interpolate;
 };

+/** \brief  Class for gather-scatter (gs) operations on Vectors based on
+    corresponding global identifiers.
+
+    This functionality is useful for gs-ops on DOF values across processor
+    boundary, where the global identifier would be the corresponding true DOF
+    index. Operations currently supported are min, max, sum, and multiplication.
+    Note: identifier 0 does not participate in the gather-scatter operation and
+    a given identifier can be included multiple times on a given rank.
+    For example, consider a vector, v:
+    - v = [0.3, 0.4, 0.25, 0.7] on rank1,
+    - v = [0.6, 0.1] on rank 2,
+    - v = [-0.2, 0.3, 0.7, 0.] on rank 3.
+
+    Consider a corresponding Array<int>, a:
+    - a = [1, 2, 3, 1] on rank 1,
+    - a = [3, 2] on rank 2,
+    - a = [1, 2, 0, 3] on rank 3.
+
+    A gather-scatter "minimum" operation, done as follows:
+    GSOPGSLIB gs = GSOPGSLIB(MPI_COMM_WORLD, a);
+    gs.GS(v, GSOp::MIN);
+    would return into v:
+    - v = [-0.2, 0.1, 0., -0.2] on rank 1,
+    - v = [0., 0.1] on rank 2,
+    - v = [-0.2, 0.1, 0.7, 0.] on rank 3,
+    where the values have been compared across all processors based on the
+    integer identifier. */
+class GSOPGSLIB
+{
+protected:
+   struct gslib::crystal *cr;               // gslib's internal data
+   struct gslib::comm *gsl_comm;            // gslib's internal data
+   struct gslib::gs_data *gsl_data = NULL;
+   int num_ids;
+
+public:
+   GSOPGSLIB(Array<long long> &ids);
+
+#ifdef MFEM_USE_MPI
+   GSOPGSLIB(MPI_Comm comm_, Array<long long> &ids);
+#endif
+
+   virtual ~GSOPGSLIB();
+
+   /// Supported operation types. See class description.
+   enum GSOp {ADD, MUL, MIN, MAX};
+
+   /// Update the identifiers used for the gather-scatter operator.
+   /// Same @a ids get grouped together and id == 0 does not participate.
+   /// See class description.
+   void UpdateIdentifiers(const Array<long long> &ids);
+
+   /// Gather-Scatter operation on senddata. Must match length of unique
+   /// identifiers used in the constructor. See class description.
+   void GS(Vector &senddata, GSOp op);
+};
+
 } // namespace mfem

 #endif // MFEM_USE_GSLIB
@@ -18,7 +18,6 @@
 namespace mfem
 {

-
 void HyperbolicFormIntegrator::AssembleElementVector(const FiniteElement &el,
                                                     ElementTransformation &Tr,
                                                     const Vector &elfun,
@@ -29,7 +28,7 @@ void HyperbolicFormIntegrator::AssembleElementVector(const FiniteElement &el,
   const int dof = el.GetDof();

 #ifdef MFEM_THREAD_SAFE
-   // Local storages for element integration
+   // Local storage for element integration

   // shape function value at an integration point
   Vector shape(dof);
@@ -62,7 +61,7 @@ void HyperbolicFormIntegrator::AssembleElementVector(const FiniteElement &el,
      ir = &IntRules.Get(Tr.GetGeometryType(), order);
   }

-   // loop over interation points
+   // loop over integration points
   for (int i = 0; i < ir->GetNPoints(); i++)
   {
      const IntegrationPoint &ip = ir->IntPoint(i);
@@ -92,7 +91,7 @@ void HyperbolicFormIntegrator::AssembleFaceVector(
   const int dof2 = el2.GetDof();

 #ifdef MFEM_THREAD_SAFE
-   // Local storages for element integration
+   // Local storage for element integration

   // shape function value at an integration point - first elem
   Vector shape1(dof1);
@@ -122,7 +121,7 @@ void HyperbolicFormIntegrator::AssembleFaceVector(
   DenseMatrix elvect2_mat(elvect.GetData() + dof1 * num_equations, dof2,
                           num_equations);

-   // obtain integration rule. If integration is rule is given, then use it.
+   // Obtain integration rule. If integration is rule is given, then use it.
   // Otherwise, get (2*p + IntOrderOffset) order integration rule
   const IntegrationRule *ir = IntRule;
   if (!ir)
@@ -149,7 +148,7 @@ void HyperbolicFormIntegrator::AssembleFaceVector(
      if (nor.Size() == 1)  // if 1D, use 1 or -1.
      {
         // This assume the 1D integration point is in (0,1). This may not work
-         // if this chages.
+         // if this changes.
         nor(0) = (Tr.GetElement1IntPoint().x - 0.5) * 2.0;
      }
      else
@@ -18,43 +18,36 @@
 namespace mfem
 {

-//                     MFEM Hyperbolic Conservation Laws
+// This file contains general hyperbolic conservation element/face form
+// integrators. HyperbolicFormIntegrator and RiemannSolver are defined.
 //
-// Description:
+// HyperbolicFormIntegrator is a NonlinearFormIntegrator that implements
+// element weak divergence and interface flux
 //
-//    This file contains general hyperbolic conservation element/face form
-//    integrators.
+//     ∫_T F(u):∇v,   -∫_e F̂(u)⋅[[v]]
 //
-//    HyperbolicFormIntegrator and RiemannSolver are defined.
-//    HyperbolicFormIntegrator is a NonlinearFormIntegrator that implements
-//    element weak divergence and interface flux
+// Here, T is an element, e is an edge, and [[⋅]] is jump. This form integrator
+// is coupled with RiemannSolver that implements the numerical flux F̂. For
+// RiemannSolver, the Rusanov flux, also known as local Lax-Friedrichs flux, is
+// provided.
 //
-//       ∫_T F(u):∇v,   -∫_e F̂(u)⋅[[v]]
+// To implement a specific hyperbolic conservation laws, users can create
+// derived classes from FluxFunction with overloaded ComputeFlux. One can
+// optionally overload ComputeFluxDotN to avoid creating dense matrix when
+// computing normal flux. Several example equations are also defined including:
+// advection, Burgers', shallow water, and Euler equations. Users can control
+// the quadrature rule by either providing the integration rule, or integration
+// order offset. Integration will use 2*p + IntOrderOffset order quadrature
+// rule.
 //
-//    Here, T is an element, e is an edge, and [[⋅]] is jump. This form
-//    integrator is coupled with RiemannSolver that implements the numerical
-//    flux F̂. For RiemannSolver, the Rusanov flux, also known as local
-//    Lax-Friedrichs flux, is provided.
-//
-//    To implement a specific hyperbolic conservation laws, users can create
-//    derived classes from FluxFunction with overloaded ComputeFlux. One can
-//    optionally overload ComputeFluxDotN to avoid creating dense matrix when
-//    computing normal flux. Several example equations are also defined
-//    including: advection, Burgers', shallow water, and Euler equations. Users
-//    can control the quadrature rule by either providing the integration rule,
-//    or integration order offset. Integration will use 2*p + IntOrderOffset
-//    order quadrature rule.
-//
-//    At each call of HyperbolicFormIntegrator::AssembleElementVector
-//    HyperbolicFormIntegrator::AssembleFaceVector, the maximum characteristic
-//    speed will be updated. This will not be reinitialized automatically.
-//    To reinitialize, use HyperbolicFormIntegrator::ResetMaxCharSpeed. See,
-//    ex18.hpp.
-//
-//    Note: To avoid communication overhead, we update the maximum
-//    characteristic speed within each process. Use a proper MPI routine to
-//    gather the information.
+// At each call of HyperbolicFormIntegrator::AssembleElementVector
+// HyperbolicFormIntegrator::AssembleFaceVector, the maximum characteristic
+// speed will be updated. This will not be reinitialized automatically. To
+// reinitialize, use HyperbolicFormIntegrator::ResetMaxCharSpeed. See, ex18.hpp.
 //
+// Note: To avoid communication overhead, we update the maximum characteristic
+// speed within each MPI process only. Use the appropriate MPI routine to gather
+// the information.

 /**
 * @brief Abstract class for hyperbolic flux for a system of hyperbolic
@@ -88,7 +81,7 @@ public:
   virtual real_t ComputeFlux(const Vector &state, ElementTransformation &Tr,
                              DenseMatrix &flux) const = 0;
   /**
-    * @brief Compute normal flux. Optionally overloadded in the
+    * @brief Compute normal flux. Optionally overloaded in the
    * derived class to avoid creating full dense matrix for flux.
    *
    * @param[in] state state at the current integration point
@@ -168,13 +161,13 @@ protected:
 class HyperbolicFormIntegrator : public NonlinearFormIntegrator
 {
 private:
-   // The maximum characterstic speed, updated during element/face vector assembly
+   // The maximum characteristic speed, updated during element/face vector assembly
   real_t max_char_speed;
   const RiemannSolver &rsolver;   // Numerical flux that maps F(u±,x) to hat(F)
   const FluxFunction &fluxFunction;
   const int IntOrderOffset; // integration order offset, 2*p + IntOrderOffset.
 #ifndef MFEM_THREAD_SAFE
-   // Local storages for element integration
+   // Local storage for element integration
   Vector shape;              // shape function value at an integration point
   Vector state;              // state value at an integration point
   DenseMatrix flux;          // flux value at an integration point
@@ -307,7 +307,7 @@ static void PADGDiffusionSetupFaceInfo2D(const int nf, const Mesh &mesh,
   }
 }

-// Assigns to perm the permuation:
+// Assigns to perm the permutation:
 //    perm[0] <- normal component
 //    perm[1] <- first tangential component
 //    perm[2] <- second tangential component
@@ -563,7 +563,7 @@ void DiffusionIntegrator::AssemblePatchMatrix_fullQuadrature(
   cdofs.SetSize(maxw[0], maxw[1], maxw[2]);

   // Compute sparsity of the sparse matrix
-   smati = new int[ndof+1];
+   smati = Memory<int>(ndof+1);
   smati[0] = 0;

   for (int dof_j=0; dof_j<ndof; ++dof_j)
@@ -586,8 +586,8 @@ void DiffusionIntegrator::AssemblePatchMatrix_fullQuadrature(
      nnz += ndd;
   }

-   smatj = new int[nnz];
-   smata = new real_t[nnz];
+   smatj = Memory<int>(nnz);
+   smata = Memory<real_t>(nnz);

   for (int i=0; i<nnz; ++i)
   {
@@ -973,7 +973,7 @@ void DiffusionIntegrator::AssemblePatchMatrix_reducedQuadrature(
   cdofs.SetSize(maxw[0], maxw[1], maxw[2]);

   // Compute sparsity of the sparse matrix
-   smati = new int[ndof+1];
+   smati = Memory<int>(ndof+1);
   smati[0] = 0;

   for (int dof_j=0; dof_j<ndof; ++dof_j)
@@ -996,8 +996,8 @@ void DiffusionIntegrator::AssemblePatchMatrix_reducedQuadrature(
      nnz += ndd;
   }

-   smatj = new int[nnz];
-   smata = new real_t[nnz];
+   smatj = Memory<int>(nnz);
+   smata = Memory<real_t>(nnz);

   for (int i=0; i<nnz; ++i)
   {
@@ -157,7 +157,7 @@ void ElasticityAddMultPA_(const int nDofs, const FiniteElementSpace &fespace,
   static constexpr int aSize = aUpper-aLower;
   static constexpr bool isComponent = (i_block >= 0);

-   //Assuming all elements are the same
+   // Assuming all elements are the same
   const auto &ir = QVec.GetIntRule(0);
   const QuadratureInterpolator *E_To_Q_Map = fespace.GetQuadratureInterpolator(
                                                 ir);
@@ -180,7 +180,7 @@ void ElasticityAddMultPA_(const int nDofs, const FiniteElementSpace &fespace,
         auto invJ = inv(make_tensor<d, d>(
         [&](int i, int j) { return J(p, i, j, e); }));
         tensor<real_t, aSize, d> gradx;
-         //load grad(x) into gradx
+         // load grad(x) into gradx
         if (isComponent)
         {
            for (int i = 0; i < d; i++)
@@ -198,11 +198,11 @@ void ElasticityAddMultPA_(const int nDofs, const FiniteElementSpace &fespace,
               }
            }
         }
-         //compute divergence
+         // compute divergence
         real_t div = 0.;
         for (int i = aLower; i < aUpper; i++)
         {
-            //take size of gradx into account
+            // take size of gradx into account
            const int iIndex = isComponent ? 0 : i;
            div += gradx(iIndex,i);
         }
@@ -211,11 +211,11 @@ void ElasticityAddMultPA_(const int nDofs, const FiniteElementSpace &fespace,
         {
            for (int q = qLower; q < qUpper; q++)
            {
-               //compute contraction of 4*sym(grad(u))sym(grad(v)) term.
-               //this contraction could be made slightly cheaper using Voigt
-               //notation, but repeated entries are summed for simplicity.
+               // compute contraction of 4*sym(grad(u))sym(grad(v)) term.
+               // this contraction could be made slightly cheaper using Voigt
+               // notation, but repeated entries are summed for simplicity.
               real_t contraction = 0.;
-               //not sure how to combine cases
+               // not sure how to combine cases
               if (isComponent)
               {
                  for (int a = 0; a < d; a++)
@@ -276,7 +276,7 @@ void ElasticityAssembleDiagonalPA_(const int nDofs,
                                   const CoefficientVector &mu, const GeometricFactors &geom,
                                   const DofToQuad &maps, QuadratureFunction &QVec, Vector &diag)
 {
-   //Assuming all elements are the same
+   // Assuming all elements are the same
   const auto &ir = QVec.GetIntRule(0);
   static constexpr int d = dim;
   const int numPoints = ir.GetNPoints();
@@ -299,9 +299,9 @@ void ElasticityAssembleDiagonalPA_(const int nDofs,
            {
               for (int q = 0; q < d; q++)
               {
-                  //compute contraction of 4*sym(grad(u))sym(grad(v)) term.
-                  //this contraction could be made slightly cheaper using Voigt
-                  //notation, but repeated entries are summed for simplicity.
+                  // compute contraction of 4*sym(grad(u))sym(grad(v)) term.
+                  // this contraction could be made slightly cheaper using Voigt
+                  // notation, but repeated entries are summed for simplicity.
                  real_t contraction = 0.;
                  for (int a = 0; a < d; a++)
                  {
@@ -321,7 +321,7 @@ void ElasticityAssembleDiagonalPA_(const int nDofs,
      }
   });

-   //Reduce quadrature function to an E-Vector
+   // Reduce quadrature function to an E-Vector
   const auto QRead = Reshape(QVec.Read(), numPoints, d, d, d, numEls);
   auto diagDev = Reshape(diag.Write(), nDofs, d, numEls);
   const auto G = Reshape(maps.G.Read(), numPoints, d, nDofs);
@@ -348,7 +348,7 @@ void ElasticityAssembleDiagonalPA_(const int nDofs,
   });
 }

-//Templated implementation of ElasticityAssembleEA.
+// Templated implementation of ElasticityAssembleEA.
 template<int dim>
 void ElasticityAssembleEA_(const int i_block,
                           const int j_block,
@@ -360,7 +360,7 @@ void ElasticityAssembleEA_(const int i_block,
                           const DofToQuad &maps,
                           Vector &emat)
 {
-   //Assuming all elements are the same
+   // Assuming all elements are the same
   static constexpr int d = dim;
   const int numPoints = ir.GetNPoints();
   const int numEls = lambda.Size()/numPoints;
@@ -386,7 +386,7 @@ void ElasticityAssembleEA_(const int i_block,
               {
                  for (int m = 0; m < d; m++)
                  {
-                     //compute contraction of 4*sym(grad(u))sym(grad(v)) term.
+                     // compute contraction of 4*sym(grad(u))sym(grad(v)) term.
                     real_t contraction = 0.;
                     for (int a = 0; a < d; a++)
                     {
@@ -101,7 +101,7 @@ void MomentFittingIntRules::InitVolume(int order, Coefficient& levelset,
         }
      }

-      // assamble the matrix
+      // assemble the matrix
      DenseMatrix Mat(nBasisVolume, ir.GetNPoints());
      for (int ip = 0; ip < ir.GetNPoints(); ip++)
      {
@@ -118,7 +118,7 @@ void MomentFittingIntRules::InitVolume(int order, Coefficient& levelset,
         Mat.SetCol(ip, shape);
      }

-      // compute the svd for the matrix
+      // compute the SVD for the matrix
      VolumeSVD = new DenseMatrixSVD(Mat, 'A', 'A');
      VolumeSVD->Eval(Mat);
   }
@@ -1239,7 +1239,7 @@ void MomentFittingIntRules::OrthoBasis2D(const IntegrationPoint& ip,

   shape.SetSize(nBasis, 2);

-   // evaluate basis inthe point
+   // evaluate basis in the point
   DenseMatrix preshape(nBasis, 2);
   DivFreeBasis2D(ip, shape);

@@ -1597,6 +1597,6 @@ void MomentFittingIntRules::GetSurfaceWeights(ElementTransformation& Tr,
   }
 }

-#endif //MFEM_USE_LAPACK
+#endif // MFEM_USE_LAPACK

 }
@@ -291,7 +291,7 @@ void BatchedLOR_AMS::FormCoordinateVectors(const Vector &X_vert)
   const auto ltdof_ldof = HypreRead(R->GetMemoryJ());

   // Go from E-vector format directly to T-vector format
-   MFEM_HYPRE_FORALL(i, ntdofs,
+   mfem::hypre_forall(ntdofs, [=] MFEM_HOST_DEVICE (int i)
   {
      const int j = d_offsets[ltdof_ldof[i]];
      for (int c = 0; c < sdim; ++c)
@@ -269,13 +269,13 @@ void BatchedLOR_H1::Assemble3D()
               real_t vx[8], vy[8], vz[8];
               LORVertexCoordinates3D<ORDER>(X, iel_ho, kx, ky, kz, vx, vy, vz);

-               //MFEM_UNROLL(2)
+               // MFEM_UNROLL(2)
               for (int iqz=0; iqz<2; ++iqz)
               {
-                  //MFEM_UNROLL(2)
+                  // MFEM_UNROLL(2)
                  for (int iqy=0; iqy<2; ++iqy)
                  {
-                     //MFEM_UNROLL(2)
+                     // MFEM_UNROLL(2)
                     for (int iqx=0; iqx<2; ++iqx)
                     {
                        const real_t x = iqx;
@@ -307,21 +307,21 @@ void BatchedLOR_H1::Assemble3D()
                  }
               }

-               //MFEM_UNROLL(2)
+               // MFEM_UNROLL(2)
               for (int iqx=0; iqx<2; ++iqx)
               {
-                  //MFEM_UNROLL(2)
+                  // MFEM_UNROLL(2)
                  for (int jz=0; jz<2; ++jz)
                  {
                     // Note loop starts at iz=jz here, taking advantage of
                     // symmetries.
-                     //MFEM_UNROLL(2)
+                     // MFEM_UNROLL(2)
                     for (int iz=jz; iz<2; ++iz)
                     {
-                        //MFEM_UNROLL(2)
+                        // MFEM_UNROLL(2)
                        for (int iqy=0; iqy<2; ++iqy)
                        {
-                           //MFEM_UNROLL(2)
+                           // MFEM_UNROLL(2)
                           for (int iqz=0; iqz<2; ++iqz)
                           {
                              const real_t mq = const_mq ? MQ(0,0,0,0) : MQ(kx+iqx, ky+iqy, kz+iqz, iel_ho);
@@ -356,10 +356,10 @@ void BatchedLOR_H1::Assemble3D()
                              real_t wdetJ = Q(6,iqz,iqy,iqx);
                              mass_A(iqy,iz,jz,iqx) += mq*wdetJ*biz*bjz;
                           }
-                           //MFEM_UNROLL(2)
+                           // MFEM_UNROLL(2)
                           for (int jy=0; jy<2; ++jy)
                           {
-                              //MFEM_UNROLL(2)
+                              // MFEM_UNROLL(2)
                              for (int iy=0; iy<2; ++iy)
                              {
                                 const real_t biy = (iy == iqy) ? 1.0 : 0.0;
@@ -382,16 +382,16 @@ void BatchedLOR_H1::Assemble3D()
                              }
                           }
                        }
-                        //MFEM_UNROLL(2)
+                        // MFEM_UNROLL(2)
                        for (int jy=0; jy<2; ++jy)
                        {
-                           //MFEM_UNROLL(2)
+                           // MFEM_UNROLL(2)
                           for (int jx=0; jx<2; ++jx)
                           {
-                              //MFEM_UNROLL(2)
+                              // MFEM_UNROLL(2)
                              for (int iy=0; iy<2; ++iy)
                              {
-                                 //MFEM_UNROLL(2)
+                                 // MFEM_UNROLL(2)
                                 for (int ix=0; ix<2; ++ix)
                                 {
                                    const real_t bix = (ix == iqx) ? 1.0 : 0.0;
@@ -431,7 +431,7 @@ void BatchedLOR_H1::Assemble3D()
               // Assemble the local matrix into the macro-element sparse matrix
               // in a format similar to coordinate format. The (I,J) arrays
               // are implicit (not stored explicitly).
-               //MFEM_UNROLL(8)
+               // MFEM_UNROLL(8)
               for (int ii_loc=0; ii_loc<nv; ++ii_loc)
               {
                  const int ix = ii_loc%2;
@@ -319,24 +319,6 @@ void L2NormalDerivativeFaceRestriction::AddMultTranspose(
 template <int T_D1D>
 void L2NormalDerivativeFaceRestriction::Mult2D(const Vector &x, Vector &y) const
 {
-   int ne_shared = 0;
-   const real_t *face_nbr_data = nullptr;
-
-#ifdef MFEM_USE_MPI
-   std::unique_ptr<ParGridFunction> x_gf;
-   if (const auto *pfes = dynamic_cast<const ParFiniteElementSpace*>(&fes))
-   {
-      if (face_type == FaceType::Interior)
-      {
-         x_gf.reset(new ParGridFunction(const_cast<ParFiniteElementSpace*>(pfes),
-                                        const_cast<Vector&>(x), 0));
-         x_gf->ExchangeFaceNbrData();
-         face_nbr_data = x_gf->FaceNbrData().Read();
-         ne_shared = pfes->GetParMesh()->GetNFaceNeighborElements();
-      }
-   }
-#endif
-
   const int vd = fes.GetVDim();
   const bool t = fes.GetOrdering() == Ordering::byVDIM;
   const int num_elem = ne;
@@ -347,6 +329,9 @@ void L2NormalDerivativeFaceRestriction::Mult2D(const Vector &x, Vector &y) const
   const int q = maps.nqpt;
   const int d = maps.ndof;

+   Vector face_nbr_data = GetLVectorFaceNbrData(fes, x, face_type);
+   const int ne_shared = face_nbr_data.Size() / d / d / vd;
+
   MFEM_VERIFY(q == d, "");
   MFEM_VERIFY(T_D1D == d || T_D1D == 0, "");

@@ -360,7 +345,7 @@ void L2NormalDerivativeFaceRestriction::Mult2D(const Vector &x, Vector &y) const
   // if byvdim, d_x has shape (vdim, nddof, nddof, ne)
   // otherwise, d_x has shape (nddof, nddof, ne, vdim)
   const auto d_x = Reshape(x.Read(), t?vd:d, d, t?d:ne, t?ne:vd);
-   const auto d_x_shared = Reshape(face_nbr_data,
+   const auto d_x_shared = Reshape(face_nbr_data.Read(),
                                   t?vd:d, d, t?d:ne_shared, t?ne_shared:vd);
   auto d_y = Reshape(y.Write(), q, vd, 2, nf);

@@ -446,24 +431,6 @@ void L2NormalDerivativeFaceRestriction::Mult2D(const Vector &x, Vector &y) const
 template <int T_D1D>
 void L2NormalDerivativeFaceRestriction::Mult3D(const Vector &x, Vector &y) const
 {
-   int ne_shared = 0;
-   const real_t *face_nbr_data = nullptr;
-
-#ifdef MFEM_USE_MPI
-   std::unique_ptr<ParGridFunction> x_gf;
-   if (const auto *pfes = dynamic_cast<const ParFiniteElementSpace*>(&fes))
-   {
-      if (face_type == FaceType::Interior)
-      {
-         x_gf.reset(new ParGridFunction(const_cast<ParFiniteElementSpace*>(pfes),
-                                        const_cast<Vector&>(x), 0));
-         x_gf->ExchangeFaceNbrData();
-         face_nbr_data = x_gf->FaceNbrData().Read();
-         ne_shared = pfes->GetParMesh()->GetNFaceNeighborElements();
-      }
-   }
-#endif
-
   const int vd = fes.GetVDim();
   const bool t = fes.GetOrdering() == Ordering::byVDIM;
   const int num_elem = ne;
@@ -475,6 +442,9 @@ void L2NormalDerivativeFaceRestriction::Mult3D(const Vector &x, Vector &y) const
   const int d = maps.ndof;
   const int q2d = q * q;

+   Vector face_nbr_data = GetLVectorFaceNbrData(fes, x, face_type);
+   const int ne_shared = face_nbr_data.Size() / d / d / d / vd;
+
   MFEM_VERIFY(q == d, "");
   MFEM_VERIFY(T_D1D == d || T_D1D == 0, "");

@@ -485,7 +455,7 @@ void L2NormalDerivativeFaceRestriction::Mult3D(const Vector &x, Vector &y) const

   // t ? (vdim, d, d, d, ne) : (d, d, d, ne, vdim)
   const auto d_x = Reshape(x.Read(), t?vd:d, d, d, t?d:ne, t?ne:vd);
-   const auto d_x_shared = Reshape(face_nbr_data,
+   const auto d_x_shared = Reshape(face_nbr_data.Read(),
                                   t?vd:d, d, d, t?d:ne_shared, t?ne_shared:vd);
   auto d_y = Reshape(y.Write(), q2d, vd, 2, nf);

@@ -368,6 +368,72 @@ const
   y.Add(a, Ytmp);
 }

+real_t ParBilinearForm::ParInnerProduct(const ParGridFunction &x,
+                                        const ParGridFunction &y) const
+{
+   MFEM_ASSERT(mat != NULL, "local matrix must be assembled");
+
+   real_t loc = InnerProduct(x, y);
+   real_t glob = 0.;
+
+   MPI_Allreduce(&loc, &glob, 1, MPITypeMap<real_t>::mpi_type, MPI_SUM,
+                 pfes->GetComm());
+
+   return glob;
+}
+
+real_t ParBilinearForm::TrueInnerProduct(const ParGridFunction &x,
+                                         const ParGridFunction &y) const
+{
+   MFEM_ASSERT(x.ParFESpace() == pfes, "the parallel spaces must match");
+   MFEM_ASSERT(y.ParFESpace() == pfes, "the parallel spaces must match");
+
+   HypreParVector *x_p = x.ParallelProject();
+   HypreParVector *y_p = y.ParallelProject();
+
+   real_t res = TrueInnerProduct(*x_p, *y_p);
+
+   delete x_p;
+   delete y_p;
+
+   return res;
+}
+
+real_t ParBilinearForm::TrueInnerProduct(HypreParVector &x,
+                                         HypreParVector &y) const
+{
+   MFEM_VERIFY(p_mat.Ptr() != NULL, "parallel matrix must be assembled");
+
+   if (p_mat->GetType() != Operator::Hypre_ParCSR)
+   {
+      return TrueInnerProduct((const Vector&)x, (const Vector&)y);
+   }
+
+   HypreParVector *Ax = new HypreParVector(pfes);
+   HypreParMatrix *A = p_mat.As<HypreParMatrix>();
+
+   A->Mult(x, *Ax);
+
+   real_t res = mfem::InnerProduct(y, *Ax);
+
+   delete Ax;
+
+   return res;
+}
+
+real_t ParBilinearForm::TrueInnerProduct(const Vector &x,
+                                         const Vector &y) const
+{
+   MFEM_VERIFY(p_mat.Ptr() != NULL, "parallel matrix must be assembled");
+
+   Vector Ax(pfes->GetTrueVSize());
+   p_mat->Mult(x, Ax);
+
+   real_t res = mfem::InnerProduct(pfes->GetComm(), y, Ax);
+
+   return res;
+}
+
 void ParBilinearForm::FormLinearSystem(
   const Array<int> &ess_tdof_list, Vector &x, Vector &b,
   OperatorHandle &A, Vector &X, Vector &B, int copy_interior)
@@ -173,6 +173,37 @@ public:
       vectors on the true dofs. */
   void TrueAddMult(const Vector &x, Vector &y, const real_t a = 1.0) const;

+   /// Compute $ y^T M x $
+   /** @warning The calculation is performed on local dofs, assuming that
+       the local vectors are consistent with the prolongations of the true
+       vectors (see ParGridFunction::Distribute()). If this is not the case,
+       use TrueInnerProduct(const ParGridFunction &, const ParGridFunction &)
+       instead.
+       @note It is assumed that the local matrix is assembled and it has
+       not been replaced by the parallel matrix through FormSystemMatrix().
+       @see TrueInnerProduct(const ParGridFunction&, const ParGridFunction&) */
+   real_t ParInnerProduct(const ParGridFunction &x,
+                          const ParGridFunction &y) const;
+
+   /// Compute $ y^T M x $ on true dofs (grid function version)
+   /** @note The ParGridFunction%s are restricted to the true-vectors for
+       for calculation.
+       @note It is assumed that the parallel system matrix is assembled,
+       see FormSystemMatrix().
+       @see ParInnerProduct(const ParGridFunction&, const ParGridFunction&) */
+   real_t TrueInnerProduct(const ParGridFunction &x,
+                           const ParGridFunction &y) const;
+
+   /// Compute $ y^T M x $ on true dofs (Hypre vector version)
+   /** @note It is assumed that the parallel system matrix is assembled,
+       see FormSystemMatrix(). */
+   real_t TrueInnerProduct(HypreParVector &x, HypreParVector &y) const;
+
+   /// Compute $ y^T M x $ on true dofs (true-vector version)
+   /** @note It is assumed that the parallel system matrix is assembled,
+       see FormSystemMatrix(). */
+   real_t TrueInnerProduct(const Vector &x, const Vector &y) const;
+
   /// Return the parallel FE space associated with the ParBilinearForm.
   ParFiniteElementSpace *ParFESpace() const { return pfes; }

@@ -861,17 +861,17 @@ void ParFiniteElementSpace::Build_Dof_TrueDof_Matrix() const // matrix P
         }
      }

-      HYPRE_Int *i_diag = new HYPRE_Int[ldof+1];
-      HYPRE_Int *j_diag = new HYPRE_Int[ltdof];
-      real_t    *d_diag = new real_t[ltdof];
+      HYPRE_Int *i_diag = Memory<HYPRE_Int>(ldof+1);
+      HYPRE_Int *j_diag = Memory<HYPRE_Int>(ltdof);
+      real_t    *d_diag = Memory<real_t>(ltdof);
      int diag_counter;

-      HYPRE_Int *i_offd = new HYPRE_Int[ldof+1];
-      HYPRE_Int *j_offd = new HYPRE_Int[nnz_offd];
-      real_t    *d_offd = new real_t[nnz_offd];
+      HYPRE_Int *i_offd = Memory<HYPRE_Int>(ldof+1);
+      HYPRE_Int *j_offd = Memory<HYPRE_Int>(nnz_offd);
+      real_t    *d_offd = Memory<real_t>(nnz_offd);
      int offd_counter;

-      HYPRE_BigInt *cmap   = new HYPRE_BigInt[ldof-ltdof];
+      HYPRE_BigInt *cmap   = Memory<HYPRE_BigInt>(ldof-ltdof);

      HYPRE_BigInt *col_starts = GetTrueDofOffsets();
      HYPRE_BigInt *row_starts = GetDofOffsets();
@@ -249,6 +249,8 @@ void ParGridFunction::ExchangeFaceNbrData()
   auto send_data_ptr = mpi_gpu_aware ? send_data.Read() : send_data.HostRead();
   auto face_nbr_data_ptr = mpi_gpu_aware ? face_nbr_data.Write() :
                            face_nbr_data.HostWrite();
+   // Wait for the kernel to be done since it updates what's sent and it may be async
+   if (mpi_gpu_aware) { MFEM_STREAM_SYNC; }
   for (int fn = 0; fn < num_face_nbrs; fn++)
   {
      int nbr_rank = pmesh->GetFaceNbrRank(fn);
@@ -518,7 +520,7 @@ void ParGridFunction::CountElementsPerVDof(Array<int> &elem_per_vdof) const
 }

 void ParGridFunction::GetDerivative(int comp, int der_comp,
-                                    ParGridFunction &der)
+                                    ParGridFunction &der) const
 {
   Array<int> overlap;
   AccumulateAndCountDerivativeValues(comp, der_comp, der, overlap);
@@ -713,10 +715,10 @@ void ParGridFunction::ProjectBdrCoefficient(
         }
      }
   }
+   gcomm.Bcast<int>(values_counter.HostReadWrite());
   for (int i = 0; i < values_counter.Size(); i++)
   {
-      MFEM_ASSERT(pfes->GetLocalTDofNumber(i) == -1 ||
-                  bool(values_counter[i]) == bool(ess_vdofs_marker[i]),
+      MFEM_ASSERT(bool(values_counter[i]) == bool(ess_vdofs_marker[i]),
                  "internal error");
   }
 #endif
@@ -753,10 +755,10 @@ void ParGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient &vcoeff,
 #ifdef MFEM_DEBUG
   Array<int> ess_vdofs_marker;
   pfes->GetEssentialVDofs(bdr_attr, ess_vdofs_marker);
+   gcomm.Bcast<int>(values_counter.HostReadWrite());
   for (int i = 0; i < values_counter.Size(); i++)
   {
-      MFEM_ASSERT(pfes->GetLocalTDofNumber(i) == -1 ||
-                  bool(values_counter[i]) == bool(ess_vdofs_marker[i]),
+      MFEM_ASSERT(bool(values_counter[i]) == bool(ess_vdofs_marker[i]),
                  "internal error: " << pfes->GetLocalTDofNumber(i) << ' ' << bool(
                     values_counter[i]));
   }
@@ -231,7 +231,7 @@ public:
   void CountElementsPerVDof(Array<int> &elem_per_vdof) const override;

   /// Parallel version of GridFunction::GetDerivative(); see its documentation.
-   void GetDerivative(int comp, int der_comp, ParGridFunction &der);
+   void GetDerivative(int comp, int der_comp, ParGridFunction &der) const;

   /** Sets the output vector @a dof_vals to the values of the degrees of
       freedom of element @a el. If @a el is greater than or equal to the number
@@ -262,7 +262,6 @@ public:
                              const Array<int> &attr) override
   { ProjectBdrCoefficient(coeff, NULL, attr); }

-   // Only the values in the master are guaranteed to be correct!
   void ProjectBdrCoefficientTangent(VectorCoefficient &vcoeff,
                                     const Array<int> &bdr_attr) override;

@@ -309,12 +309,8 @@ void ParL2FaceRestriction::DoubleValuedConformingMult(
   MFEM_ASSERT(
      m == L2FaceValues::DoubleValued,
      "This method should be called when m == L2FaceValues::DoubleValued.");
-   ParGridFunction x_gf;
-   x_gf.MakeRef(const_cast<ParFiniteElementSpace*>(&pfes),
-                const_cast<Vector&>(x), 0);
-   // Face-neighbor information is only needed for interior faces. For boundary
-   // faces, no communication is required.
-   if (type == FaceType::Interior) { x_gf.ExchangeFaceNbrData(); }
+
+   Vector face_nbr_data = GetLVectorFaceNbrData(fes, x, type);

   // Early return only after calling ParGridFunction::ExchangeFaceNbrData,
   // otherwise MPI communication can hang.
@@ -329,7 +325,7 @@ void ParL2FaceRestriction::DoubleValuedConformingMult(
   auto d_indices1 = scatter_indices1.Read();
   auto d_indices2 = scatter_indices2.Read();
   auto d_x = Reshape(x.Read(), t?vd:ndofs, t?ndofs:vd);
-   auto d_x_shared = Reshape(x_gf.FaceNbrData().Read(),
+   auto d_x_shared = Reshape(face_nbr_data.Read(),
                             t?vd:nsdofs, t?nsdofs:vd);
   auto d_y = Reshape(y.Write(), nface_dofs, vd, 2, nf);
   mfem::forall(nfdofs, [=] MFEM_HOST_DEVICE (int i)
@@ -13,7 +13,7 @@
 #include "normal_deriv_restriction.hpp"
 #include "gridfunc.hpp"
 #include "fespace.hpp"
-#include "pfespace.hpp"
+#include "pgridfunc.hpp"
 #include "qspace.hpp"
 #include "fe/face_map_utils.hpp"
 #include "../general/forall.hpp"
@@ -2282,4 +2282,33 @@ void NCL2FaceRestriction::ComputeGatherIndices()
   gather_offsets[0] = 0;
 }

+Vector GetLVectorFaceNbrData(
+   const FiniteElementSpace &fes, const Vector &x, FaceType ftype)
+{
+#ifdef MFEM_USE_MPI
+   if (ftype == FaceType::Interior)
+   {
+      if (auto *pfes = const_cast<ParFiniteElementSpace*>
+                       (dynamic_cast<const ParFiniteElementSpace*>(&fes)))
+      {
+         if (auto *x_gf = const_cast<ParGridFunction*>
+                          (dynamic_cast<const ParGridFunction*>(&x)))
+         {
+            Vector &gf_face_nbr = x_gf->FaceNbrData();
+            if (gf_face_nbr.Size() == 0) { x_gf->ExchangeFaceNbrData(); }
+            gf_face_nbr.Read();
+            return Vector(gf_face_nbr, 0, gf_face_nbr.Size());
+         }
+         else
+         {
+            ParGridFunction gf(pfes, const_cast<Vector&>(x));
+            gf.ExchangeFaceNbrData();
+            return std::move(gf.FaceNbrData());
+         }
+      }
+   }
+#endif
+   return Vector();
+}
+
 } // namespace mfem
@@ -1094,6 +1094,20 @@ int PermuteFaceL2(const int dim, const int face_id1,
                  const int face_id2, const int orientation,
                  const int size1d, const int index);

+/// @brief Return the face-neighbor data given the L-vector @a x.
+///
+/// If the input vector @a x is a ParGridFunction with non-empty face-neighbor
+/// data, return an alias to ParGridFunction::FaceNbrData() (avoiding an
+/// unneeded call to ParGridFunction::ExchangeFaceNbrData).
+///
+/// Otherwise, create a temporary ParGridFunction, exchange the face-neighbor
+/// data, and return the resulting vector.
+///
+/// If @a fes is not a parallel space, or if @a ftype is not FaceType::Interior,
+/// return an empty vector.
+Vector GetLVectorFaceNbrData(
+   const FiniteElementSpace &fes, const Vector &x, FaceType ftype);
+
 }

 #endif // MFEM_RESTRICTION
@@ -503,7 +503,7 @@ public:

      Array<int> vdofs;
      const Array<int> *dof_map = sol_fe.GetDofMap();
-      const int *dof_map_ = dof_map->GetData();
+      const int *dof_map_ = (dof_map) ? dof_map->GetData() : NULL;
      DenseMatrix M_loc_perm(dofs*vdim,dofs*vdim); // initialized with zeros

      const int NE = mesh.GetNE();
@@ -98,7 +98,7 @@ protected:
      const int vsize = sizeof(vint_t)/sizeof(attrib[0][0]);
      for (int i = 0; i < NE; i++)
      {
-         for (int j = 0; j < vsize; i++)
+         for (int j = 0; j < vsize; j++)
         {
            attrib[i][j] = elements[el+j+i*vsize]->GetAttribute();
         }
@@ -2949,6 +2949,15 @@ void TMOP_Integrator::EnableSurfaceFitting(const GridFunction &s0,
   MFEM_VERIFY(surf_fit_pos == NULL,
               "Using both fitting approaches is not supported.");

+   const int dim = s0.FESpace()->GetMesh()->Dimension();
+   Mesh *mesh = s0.FESpace()->GetMesh();
+   MFEM_VERIFY(mesh->GetNodes()->Size() == dim*s0.Size(),
+               "Mesh and level-set polynomial order must be the same.");
+   const H1_FECollection *fec = dynamic_cast<const H1_FECollection *>
+                                (s0.FESpace()->FEColl());
+   MFEM_VERIFY(fec, "Only H1_FECollection is supported for the surface fitting "
+               "grid function.");
+
   delete surf_fit_gf;
   surf_fit_gf = new GridFunction(s0);
   surf_fit_gf->CountElementsPerVDof(surf_fit_dof_count);
@@ -2987,12 +2996,24 @@ void TMOP_Integrator::EnableSurfaceFitting(const GridFunction &pos,
 void TMOP_Integrator::EnableSurfaceFitting(const ParGridFunction &s0,
                                           const Array<bool> &smarker,
                                           Coefficient &coeff,
-                                           AdaptivityEvaluator &ae)
+                                           AdaptivityEvaluator &ae,
+                                           AdaptivityEvaluator *aegrad,
+                                           AdaptivityEvaluator *aehess)
 {
   // To have both we must duplicate the markers.
   MFEM_VERIFY(surf_fit_pos == NULL,
               "Using both fitting approaches is not supported.");

+   const int dim = s0.FESpace()->GetMesh()->Dimension();
+   ParMesh *pmesh = s0.ParFESpace()->GetParMesh();
+   MFEM_VERIFY(pmesh->GetNodes()->Size() == dim*s0.Size(),
+               "Mesh and level-set polynomial order must be the same.");
+   const H1_FECollection *fec = dynamic_cast<const H1_FECollection *>
+                                (s0.FESpace()->FEColl());
+   MFEM_VERIFY(fec, "Only H1_FECollection is supported for the surface fitting "
+               "grid function.");
+
+
   delete surf_fit_gf;
   surf_fit_gf = new GridFunction(s0);
   s0.CountElementsPerVDof(surf_fit_dof_count);
@@ -3000,11 +3021,80 @@ void TMOP_Integrator::EnableSurfaceFitting(const ParGridFunction &s0,
   surf_fit_coeff = &coeff;
   surf_fit_eval = &ae;

-   surf_fit_eval->SetParMetaInfo(*s0.ParFESpace()->GetParMesh(),
-                                 *s0.ParFESpace());
+   surf_fit_eval->SetParMetaInfo(*pmesh, *s0.ParFESpace());
   surf_fit_eval->SetInitialField
   (*surf_fit_gf->FESpace()->GetMesh()->GetNodes(), *surf_fit_gf);
-   surf_fit_gf_bg = false;
+
+   if (!aegrad) { return; }
+
+   MFEM_VERIFY(aehess, "AdaptivityEvaluator for Hessians must be provided too.");
+
+   ParFiniteElementSpace *fes = s0.ParFESpace();
+
+   // FE space for gradients.
+   delete surf_fit_grad;
+   H1_FECollection *fec_grad = new H1_FECollection(fec->GetOrder(), dim,
+                                                   fec->GetBasisType());
+   ParFiniteElementSpace *fes_grad = new ParFiniteElementSpace(pmesh, fec_grad,
+                                                               dim);
+   // Initial gradients.
+   surf_fit_grad = new GridFunction(fes_grad);
+   surf_fit_grad->MakeOwner(fec_grad);
+   for (int d = 0; d < dim; d++)
+   {
+      ParGridFunction surf_fit_grad_comp(fes, surf_fit_grad->GetData()+d*s0.Size());
+      s0.GetDerivative(1, d, surf_fit_grad_comp);
+   }
+   surf_fit_eval_grad = aegrad;
+   surf_fit_eval_grad->SetParMetaInfo(*pmesh, *fes_grad);
+   surf_fit_eval_grad->SetInitialField(*pmesh->GetNodes(), *surf_fit_grad);
+
+   // FE space for Hessians.
+   delete surf_fit_hess;
+   H1_FECollection *fec_hess = new H1_FECollection(fec->GetOrder(), dim,
+                                                   fec->GetBasisType());
+   ParFiniteElementSpace *fes_hess = new ParFiniteElementSpace(pmesh, fec_hess,
+                                                               dim*dim);
+   // Initial Hessians.
+   surf_fit_hess = new GridFunction(fes_hess);
+   surf_fit_hess->MakeOwner(fec_hess);
+   int id = 0;
+   for (int d = 0; d < dim; d++)
+   {
+      for (int idir = 0; idir < dim; idir++)
+      {
+         ParGridFunction surf_fit_grad_comp(fes,
+                                            surf_fit_grad->GetData()+d*s0.Size());
+         ParGridFunction surf_fit_hess_comp(fes,
+                                            surf_fit_hess->GetData()+id*s0.Size());
+         surf_fit_grad_comp.GetDerivative(1, idir, surf_fit_hess_comp);
+         id++;
+      }
+   }
+   surf_fit_eval_hess = aehess;
+   surf_fit_eval_hess->SetParMetaInfo(*pmesh, *fes_hess);
+   surf_fit_eval_hess->SetInitialField(*pmesh->GetNodes(), *surf_fit_hess);
+
+   // Store DOF indices that are marked for fitting. Used to reduce work for
+   // transferring information between source/background and current mesh.
+   surf_fit_marker_dof_index.SetSize(0);
+#ifdef MFEM_USE_GSLIB
+   if (dynamic_cast<InterpolatorFP *>(surf_fit_eval) &&
+       dynamic_cast<InterpolatorFP *>(surf_fit_eval_grad) &&
+       dynamic_cast<InterpolatorFP *>(surf_fit_eval_hess))
+   {
+      for (int i = 0; i < surf_fit_marker->Size(); i++)
+      {
+         if ((*surf_fit_marker)[i] == true)
+         {
+            surf_fit_marker_dof_index.Append(i);
+         }
+      }
+   }
+#endif
+
+   *surf_fit_grad = 0.0;
+   *surf_fit_hess = 0.0;
 }

 void TMOP_Integrator::EnableSurfaceFittingFromSource(
@@ -3022,16 +3112,17 @@ void TMOP_Integrator::EnableSurfaceFittingFromSource(
   // Setup for level set function
   delete surf_fit_gf;
   surf_fit_gf = new GridFunction(s0);
-   *surf_fit_gf = 0.0;
   surf_fit_marker = &smarker;
   surf_fit_coeff = &coeff;
   surf_fit_eval = &ae;

-   surf_fit_gf_bg = true;
   surf_fit_eval->SetParMetaInfo(*s_bg.ParFESpace()->GetParMesh(),
                                 *s_bg.ParFESpace());
   surf_fit_eval->SetInitialField
   (*s_bg.FESpace()->GetMesh()->GetNodes(), s_bg);
+   GridFunction *nodes = s0.FESpace()->GetMesh()->GetNodes();
+   surf_fit_eval->ComputeAtNewPosition(*nodes, *surf_fit_gf,
+                                       nodes->FESpace()->GetOrdering());

   // Setup for gradient on background mesh
   MFEM_VERIFY(s_bg_grad.ParFESpace()->GetOrdering() ==
@@ -3041,11 +3132,11 @@ void TMOP_Integrator::EnableSurfaceFittingFromSource(
   delete surf_fit_grad;
   surf_fit_grad = new GridFunction(s0_grad);
   *surf_fit_grad = 0.0;
-   surf_fit_eval_bg_grad = &age;
-   surf_fit_eval_bg_hess = &ahe;
-   surf_fit_eval_bg_grad->SetParMetaInfo(*s_bg_grad.ParFESpace()->GetParMesh(),
-                                         *s_bg_grad.ParFESpace());
-   surf_fit_eval_bg_grad->SetInitialField
+   surf_fit_eval_grad = &age;
+   surf_fit_eval_hess = &ahe;
+   surf_fit_eval_grad->SetParMetaInfo(*s_bg_grad.ParFESpace()->GetParMesh(),
+                                      *s_bg_grad.ParFESpace());
+   surf_fit_eval_grad->SetInitialField
   (*s_bg_grad.FESpace()->GetMesh()->GetNodes(), s_bg_grad);

   // Setup for Hessian on background mesh
@@ -3056,9 +3147,9 @@ void TMOP_Integrator::EnableSurfaceFittingFromSource(
   delete surf_fit_hess;
   surf_fit_hess = new GridFunction(s0_hess);
   *surf_fit_hess = 0.0;
-   surf_fit_eval_bg_hess->SetParMetaInfo(*s_bg_hess.ParFESpace()->GetParMesh(),
-                                         *s_bg_hess.ParFESpace());
-   surf_fit_eval_bg_hess->SetInitialField
+   surf_fit_eval_hess->SetParMetaInfo(*s_bg_hess.ParFESpace()->GetParMesh(),
+                                      *s_bg_hess.ParFESpace());
+   surf_fit_eval_hess->SetInitialField
   (*s_bg_hess.FESpace()->GetMesh()->GetNodes(), s_bg_hess);

   // Count number of zones that share each of the DOFs
@@ -3863,7 +3954,7 @@ void TMOP_Integrator::AssembleElemVecSurfFit(const FiniteElement &el_x,

   Vector sigma_e(dof_s);
   DenseMatrix surf_fit_grad_e(dof_s, dim);
-   if (surf_fit_gf || surf_fit_gf_bg)
+   if (surf_fit_gf)
   {
      surf_fit_gf->GetSubVector(vdofs, sigma_e);

@@ -3871,7 +3962,7 @@ void TMOP_Integrator::AssembleElemVecSurfFit(const FiniteElement &el_x,
      // The FE coefficients of the gradient go in surf_fit_grad_e.
      Vector grad_ptr(surf_fit_grad_e.GetData(), dof_s * dim);
      DenseMatrix grad_phys; // This will be (dof x dim, dof).
-      if (surf_fit_gf_bg)
+      if (surf_fit_grad)
      {
         surf_fit_grad->FESpace()->GetElementVDofs(el_id, dofs);
         surf_fit_grad->GetSubVector(dofs, grad_ptr);
@@ -3945,7 +4036,7 @@ void TMOP_Integrator::AssembleElemGradSurfFit(const FiniteElement &el_x,
   Vector sigma_e(dof_s);
   DenseMatrix surf_fit_grad_e(dof_s, dim);
   DenseMatrix surf_fit_hess_e(dof_s, dim*dim);
-   if (surf_fit_gf || surf_fit_gf_bg)
+   if (surf_fit_gf)
   {
      surf_fit_gf->GetSubVector(vdofs, sigma_e);

@@ -3953,7 +4044,7 @@ void TMOP_Integrator::AssembleElemGradSurfFit(const FiniteElement &el_x,
      // The FE coefficients of the gradient go in surf_fit_grad_e.
      Vector grad_ptr(surf_fit_grad_e.GetData(), dof_s * dim);
      DenseMatrix grad_phys; // This will be (dof x dim, dof).
-      if (surf_fit_gf_bg)
+      if (surf_fit_grad)
      {
         surf_fit_grad->FESpace()->GetElementVDofs(el_id, dofs);
         surf_fit_grad->GetSubVector(dofs, grad_ptr);
@@ -3967,7 +4058,7 @@ void TMOP_Integrator::AssembleElemGradSurfFit(const FiniteElement &el_x,
      // Project the Hessian of sigma in the same space.
      // The FE coefficients of the Hessian go in surf_fit_hess_e.
      Vector hess_ptr(surf_fit_hess_e.GetData(), dof_s*dim*dim);
-      if (surf_fit_gf_bg)
+      if (surf_fit_hess)
      {
         surf_fit_hess->FESpace()->GetElementVDofs(el_id, dofs);
         surf_fit_hess->GetSubVector(dofs, hess_ptr);
@@ -3994,7 +4085,7 @@ void TMOP_Integrator::AssembleElemGradSurfFit(const FiniteElement &el_x,
      Tpr.SetIntPoint(&ip);
      real_t w = surf_fit_normal * surf_fit_coeff->Eval(Tpr, ip);

-      if (surf_fit_gf || surf_fit_gf_bg)
+      if (surf_fit_gf)
      {
         Vector gg_ptr(surf_fit_hess_s.GetData(), dim * dim);
         surf_fit_hess_e.GetRow(s, gg_ptr);
@@ -4376,6 +4467,130 @@ void TMOP_Integrator::ComputeMinJac(const Vector &x,
   dx = detv_avg_min / dxscale;
 }

+void TMOP_Integrator::RemapSurfaceFittingLevelSetAtNodes(const Vector &new_x,
+                                                         int new_x_ordering)
+{
+   if (!surf_fit_gf) { return; }
+
+   if (surf_fit_marker_dof_index.Size())
+   {
+      // Interpolate information only at DOFs marked for fitting.
+      const int dim = surf_fit_gf->FESpace()->GetMesh()->Dimension();
+      const int cnt = surf_fit_marker_dof_index.Size();
+      const int total_cnt = new_x.Size()/dim;
+      Vector new_x_sorted(cnt*dim);
+      if (new_x_ordering == 0)
+      {
+         for (int d = 0; d < dim; d++)
+         {
+            for (int i = 0; i < cnt; i++)
+            {
+               int dof_index = surf_fit_marker_dof_index[i];
+               new_x_sorted(i + d*cnt) = new_x(dof_index + d*total_cnt);
+            }
+         }
+      }
+      else
+      {
+         for (int i = 0; i < cnt; i++)
+         {
+            int dof_index = surf_fit_marker_dof_index[i];
+            for (int d = 0; d < dim; d++)
+            {
+               new_x_sorted(d + i*dim) = new_x(d + dof_index*dim);
+            }
+         }
+      }
+
+      // Interpolate values of the LS.
+      Vector surf_fit_gf_int, surf_fit_grad_int, surf_fit_hess_int;
+      surf_fit_eval->ComputeAtNewPosition(new_x_sorted, surf_fit_gf_int,
+                                          new_x_ordering);
+      for (int i = 0; i < cnt; i++)
+      {
+         int dof_index = surf_fit_marker_dof_index[i];
+         (*surf_fit_gf)[dof_index] = surf_fit_gf_int(i);
+      }
+
+      // Interpolate gradients of the LS.
+      surf_fit_eval_grad->ComputeAtNewPosition(new_x_sorted, surf_fit_grad_int,
+                                               new_x_ordering);
+      // Assumes surf_fit_grad and surf_fit_gf share the same space
+      const int grad_dim = surf_fit_grad->VectorDim();
+      const int grad_cnt = surf_fit_grad->Size()/grad_dim;
+      if (surf_fit_grad->FESpace()->GetOrdering() == Ordering::byNODES)
+      {
+         for (int d = 0; d < grad_dim; d++)
+         {
+            for (int i = 0; i < cnt; i++)
+            {
+               int dof_index = surf_fit_marker_dof_index[i];
+               (*surf_fit_grad)[dof_index + d*grad_cnt] =
+                  surf_fit_grad_int(i + d*cnt);
+            }
+         }
+      }
+      else
+      {
+         for (int i = 0; i < cnt; i++)
+         {
+            int dof_index = surf_fit_marker_dof_index[i];
+            for (int d = 0; d < grad_dim; d++)
+            {
+               (*surf_fit_grad)[dof_index*grad_dim + d] =
+                  surf_fit_grad_int(i*grad_dim + d);
+            }
+         }
+      }
+
+      // Interpolate Hessians of the LS.
+      surf_fit_eval_hess->ComputeAtNewPosition(new_x_sorted, surf_fit_hess_int,
+                                               new_x_ordering);
+      // Assumes surf_fit_hess and surf_fit_gf share the same space
+      const int hess_dim = surf_fit_hess->VectorDim();
+      const int hess_cnt = surf_fit_hess->Size()/hess_dim;
+      if (surf_fit_hess->FESpace()->GetOrdering() == Ordering::byNODES)
+      {
+         for (int d = 0; d < hess_dim; d++)
+         {
+            for (int i = 0; i < cnt; i++)
+            {
+               int dof_index = surf_fit_marker_dof_index[i];
+               (*surf_fit_hess)[dof_index + d*hess_cnt] =
+                  surf_fit_hess_int(i + d*cnt);
+            }
+         }
+      }
+      else
+      {
+         for (int i = 0; i < cnt; i++)
+         {
+            int dof_index = surf_fit_marker_dof_index[i];
+            for (int d = 0; d < hess_dim; d++)
+            {
+               (*surf_fit_hess)[dof_index*hess_dim + d] =
+                  surf_fit_hess_int(i*hess_dim + d);
+            }
+         }
+      }
+
+   }
+   else
+   {
+      surf_fit_eval->ComputeAtNewPosition(new_x, *surf_fit_gf, new_x_ordering);
+      if (surf_fit_eval_grad)
+      {
+         surf_fit_eval_grad->ComputeAtNewPosition(new_x, *surf_fit_grad,
+                                                  new_x_ordering);
+      }
+      if (surf_fit_eval_hess)
+      {
+         surf_fit_eval_hess->ComputeAtNewPosition(new_x, *surf_fit_hess,
+                                                  new_x_ordering);
+      }
+   }
+}
+
 void TMOP_Integrator::
 UpdateAfterMeshPositionChange(const Vector &x_new,
                              const FiniteElementSpace &x_fes)
@@ -4406,112 +4621,11 @@ UpdateAfterMeshPositionChange(const Vector &x_new,
      adapt_lim_eval->ComputeAtNewPosition(x_new, *adapt_lim_gf, ordering);
   }

-   // Update surf_fit_gf if surface fitting is enabled.
+   // Update surf_fit_gf (and optionally its gradients) if surface
+   // fitting is enabled.
   if (surf_fit_gf)
   {
-      if (surf_fit_gf_bg)
-      {
-         // Interpolate information for only DOFs marked for fitting.
-         const int dim = surf_fit_gf->FESpace()->GetMesh()->Dimension();
-         const int cnt = surf_fit_marker_dof_index.Size();
-         const int total_cnt = x_new.Size()/dim;
-         Vector new_x_sorted(cnt*dim);
-         if (ordering == 0)
-         {
-            for (int d = 0; d < dim; d++)
-            {
-               for (int i = 0; i < cnt; i++)
-               {
-                  int dof_index = surf_fit_marker_dof_index[i];
-                  new_x_sorted(i + d*cnt) = x_new(dof_index + d*total_cnt);
-               }
-            }
-         }
-         else
-         {
-            for (int i = 0; i < cnt; i++)
-            {
-               int dof_index = surf_fit_marker_dof_index[i];
-               for (int d = 0; d < dim; d++)
-               {
-                  new_x_sorted(d + i*dim) = x_new(d + dof_index*dim);
-               }
-            }
-         }
-
-         Vector surf_fit_gf_int, surf_fit_grad_int, surf_fit_hess_int;
-         surf_fit_eval->ComputeAtNewPosition(
-            new_x_sorted, surf_fit_gf_int, ordering);
-         for (int i = 0; i < cnt; i++)
-         {
-            int dof_index = surf_fit_marker_dof_index[i];
-            (*surf_fit_gf)[dof_index] = surf_fit_gf_int(i);
-         }
-
-         surf_fit_eval_bg_grad->ComputeAtNewPosition(
-            new_x_sorted, surf_fit_grad_int, ordering);
-         // Assumes surf_fit_grad and surf_fit_gf share the same space
-         const int grad_dim = surf_fit_grad->VectorDim();
-         const int grad_cnt = surf_fit_grad->Size()/grad_dim;
-         if (surf_fit_grad->FESpace()->GetOrdering() == Ordering::byNODES)
-         {
-            for (int d = 0; d < grad_dim; d++)
-            {
-               for (int i = 0; i < cnt; i++)
-               {
-                  int dof_index = surf_fit_marker_dof_index[i];
-                  (*surf_fit_grad)[dof_index + d*grad_cnt] =
-                     surf_fit_grad_int(i + d*cnt);
-               }
-            }
-         }
-         else
-         {
-            for (int i = 0; i < cnt; i++)
-            {
-               int dof_index = surf_fit_marker_dof_index[i];
-               for (int d = 0; d < grad_dim; d++)
-               {
-                  (*surf_fit_grad)[dof_index*dim + d] =
-                     surf_fit_grad_int(i*dim + d);
-               }
-            }
-         }
-
-         surf_fit_eval_bg_hess->ComputeAtNewPosition(
-            new_x_sorted, surf_fit_hess_int, ordering);
-         // Assumes surf_fit_hess and surf_fit_gf share the same space
-         const int hess_dim = surf_fit_hess->VectorDim();
-         const int hess_cnt = surf_fit_hess->Size()/hess_dim;
-         if (surf_fit_hess->FESpace()->GetOrdering() == Ordering::byNODES)
-         {
-            for (int d = 0; d < hess_dim; d++)
-            {
-               for (int i = 0; i < cnt; i++)
-               {
-                  int dof_index = surf_fit_marker_dof_index[i];
-                  (*surf_fit_hess)[dof_index + d*hess_cnt] =
-                     surf_fit_hess_int(i + d*cnt);
-               }
-            }
-         }
-         else
-         {
-            for (int i = 0; i < cnt; i++)
-            {
-               int dof_index = surf_fit_marker_dof_index[i];
-               for (int d = 0; d < hess_dim; d++)
-               {
-                  (*surf_fit_hess)[dof_index*dim + d] =
-                     surf_fit_hess_int(i*dim + d);
-               }
-            }
-         }
-      }
-      else
-      {
-         surf_fit_eval->ComputeAtNewPosition(x_new, *surf_fit_gf, ordering);
-      }
+      RemapSurfaceFittingLevelSetAtNodes(x_new, ordering);
   }
 }

@@ -1784,12 +1784,11 @@ protected:
   // Fitting to given physical positions.
   TMOP_QuadraticLimiter *surf_fit_limiter; // Owned. Created internally.
   const GridFunction *surf_fit_pos;        // Not owned. Positions to fit.
-   real_t surf_fit_normal;
-   bool surf_fit_gf_bg;
-   GridFunction *surf_fit_grad, *surf_fit_hess;
-   AdaptivityEvaluator *surf_fit_eval_bg_grad, *surf_fit_eval_bg_hess;
-   Array<int> surf_fit_dof_count;
-   Array<int> surf_fit_marker_dof_index;
+   real_t surf_fit_normal;                  // Normalization factor.
+   GridFunction *surf_fit_grad, *surf_fit_hess; // Owned. Created internally.
+   AdaptivityEvaluator *surf_fit_eval_grad, *surf_fit_eval_hess; // Not owned.
+   Array<int> surf_fit_dof_count;            // Number of dofs per node.
+   Array<int> surf_fit_marker_dof_index;     // Indices of nodes to fit.

   DiscreteAdaptTC *discr_tc;

@@ -1985,6 +1984,10 @@ protected:
   real_t ComputeUntanglerMaxMuBarrier(const Vector &x,
                                       const FiniteElementSpace &fes);

+   // Remaps the internal surface fitting gridfunction object at provided
+   // locations.
+   void RemapSurfaceFittingLevelSetAtNodes(const Vector &new_x,
+                                           int new_x_ordering);
 public:
   /** @param[in] m    TMOP_QualityMetric for r-adaptivity (not owned).
       @param[in] tc   Target-matrix construction algorithm to use (not owned).
@@ -2000,9 +2003,8 @@ public:
        surf_fit_marker(NULL), surf_fit_coeff(NULL),
        surf_fit_gf(NULL), surf_fit_eval(NULL),
        surf_fit_limiter(NULL), surf_fit_pos(NULL),
-        surf_fit_normal(1.0),
-        surf_fit_gf_bg(false), surf_fit_grad(NULL), surf_fit_hess(NULL),
-        surf_fit_eval_bg_grad(NULL), surf_fit_eval_bg_hess(NULL),
+        surf_fit_normal(1.0), surf_fit_grad(NULL), surf_fit_hess(NULL),
+        surf_fit_eval_grad(NULL), surf_fit_eval_hess(NULL),
        discr_tc(dynamic_cast<DiscreteAdaptTC *>(tc)),
        fdflag(false), dxscale(1.0e3), fd_call_flag(false), exact_action(false)
   { PA.enabled = false; }
@@ -2103,9 +2105,15 @@ public:

 #ifdef MFEM_USE_MPI
   /// Parallel support for surface fitting to the zero level set of a function.
+   /// Here, we add two optional inputs: @a aegrad and @a aehess. When provided,
+   /// the first and second derivative of the input level set are computed on
+   /// the initial mesh, and @a aegrad and @a aehess are used to remap grad_s(x)
+   /// from grad_s0(x0) and hess_s(x) from hess_s0(x0), respectively.
   void EnableSurfaceFitting(const ParGridFunction &s0,
                             const Array<bool> &smarker, Coefficient &coeff,
-                             AdaptivityEvaluator &ae);
+                             AdaptivityEvaluator &ae,
+                             AdaptivityEvaluator *aegrad = NULL,
+                             AdaptivityEvaluator *aehess = NULL);

   /** @brief Fitting of certain DOFs in the current mesh to the zero level set
       of a function defined on another (finer) source mesh.
@@ -429,11 +429,13 @@ real_t TMOPNewtonSolver::ComputeScalingFactor(const Vector &x,
 #endif

   real_t scale = 1.0;
-   real_t avg_surf_fit_err, max_surf_fit_err = 0.0;
-   if (surf_fit_max_threshold > 0.0)
+   bool fitting = IsSurfaceFittingEnabled();
+   real_t init_fit_avg_err, init_fit_max_err = 0.0;
+   if (fitting && surf_fit_converge_error)
   {
-      GetSurfaceFittingError(x_out_loc, avg_surf_fit_err, max_surf_fit_err);
-      if (max_surf_fit_err < surf_fit_max_threshold)
+      GetSurfaceFittingError(x_out_loc, init_fit_avg_err, init_fit_max_err);
+      // Check for convergence
+      if (init_fit_max_err < surf_fit_max_err_limit)
      {
         if (print_options.iterations)
         {
@@ -444,11 +446,12 @@ real_t TMOPNewtonSolver::ComputeScalingFactor(const Vector &x,
         return scale;
      }
   }
-   if (adapt_inc_count >= max_adapt_inc_count)
+
+   if (surf_fit_adapt_count >= surf_fit_adapt_count_limit)
   {
      if (print_options.iterations)
      {
-         mfem::out << "TMOPNewtonSolver converged "
+         mfem::out << "TMOPNewtonSolver terminated "
                   "based on max number of times surface fitting weight can"
                   "be increased. \n";
      }
@@ -467,7 +470,7 @@ real_t TMOPNewtonSolver::ComputeScalingFactor(const Vector &x,
      // reference to detect deteriorations.
      MFEM_VERIFY(min_det_ptr != NULL, " Initial mesh was valid, but"
                  " intermediate mesh is invalid. Contact TMOP Developers.");
-      MFEM_VERIFY(min_detJ_threshold == 0.0,
+      MFEM_VERIFY(min_detJ_limit == 0.0,
                  "This setup is not supported. Contact TMOP Developers.");
      *min_det_ptr = untangle_factor * min_detT_in;
   }
@@ -478,6 +481,7 @@ real_t TMOPNewtonSolver::ComputeScalingFactor(const Vector &x,
   bool x_out_ok = false;
   real_t energy_out = 0.0, min_detT_out;
   const real_t norm_in = Norm(r);
+   real_t avg_fit_err, max_fit_err = 0.0;

   const real_t detJ_factor = (solver_type == 1) ? 0.25 : 0.5;
   compute_metric_quantile_flag = false;
@@ -488,6 +492,9 @@ real_t TMOPNewtonSolver::ComputeScalingFactor(const Vector &x,
   // Perform the line search.
   for (int i = 0; i < 12; i++)
   {
+      avg_fit_err = 0.0;
+      max_fit_err = 0.0;
+
      // Update the mesh and get the L-vector in x_out_loc.
      add(x, -scale, c, x_out);
      if (serial)
@@ -502,7 +509,7 @@ real_t TMOPNewtonSolver::ComputeScalingFactor(const Vector &x,

      // Check the changes in detJ.
      min_detT_out = ComputeMinDet(x_out_loc, *fes);
-      if (untangling == false && min_detT_out <= min_detJ_threshold)
+      if (untangling == false && min_detT_out <= min_detJ_limit)
      {
         // No untangling, and detJ got negative (or small) -- no good.
         if (print_options.iterations)
@@ -529,18 +536,19 @@ real_t TMOPNewtonSolver::ComputeScalingFactor(const Vector &x,
      // Check the changes in total energy.
      ProcessNewState(x_out);

-      real_t avg_fit_err, max_fit_err = 0.0;
-      if (surf_fit_max_threshold > 0.0)
+      // Ensure sufficient decrease in fitting error if we are trying to
+      // converge based on error.
+      if (fitting && surf_fit_converge_error)
      {
         GetSurfaceFittingError(x_out_loc, avg_fit_err, max_fit_err);
-      }
-      if (surf_fit_max_threshold > 0.0 && max_fit_err >= 1.2*max_surf_fit_err)
-      {
-         if (print_options.iterations)
+         if (max_fit_err >= 1.2*init_fit_max_err)
         {
-            mfem::out << "Scale = " << scale << " Surf fit err increased.\n";
+            if (print_options.iterations)
+            {
+               mfem::out << "Scale = " << scale << " Surf fit err increased.\n";
+            }
+            scale *= 0.5; continue;
         }
-         scale *= 0.5; continue;
      }

      if (serial)
@@ -614,7 +622,7 @@ real_t TMOPNewtonSolver::ComputeScalingFactor(const Vector &x,

   if (x_out_ok == false) { scale = 0.0; }

-   if (surf_fit_scale_factor > 0.0) { update_surf_fit_coeff = true; }
+   if (surf_fit_scale_factor > 0.0) { surf_fit_coeff_update = true; }
   compute_metric_quantile_flag = true;

   return scale;
@@ -657,7 +665,7 @@ void TMOPNewtonSolver::GetSurfaceFittingWeight(Array<real_t> &weights) const
   for (int i = 0; i < integs.Size(); i++)
   {
      ti = dynamic_cast<TMOP_Integrator *>(integs[i]);
-      if (ti)
+      if (ti && ti->IsSurfaceFittingEnabled())
      {
         weight = ti->GetSurfaceFittingWeight();
         weights.Append(weight);
@@ -668,8 +676,11 @@ void TMOPNewtonSolver::GetSurfaceFittingWeight(Array<real_t> &weights) const
         Array<TMOP_Integrator *> ati = co->GetTMOPIntegrators();
         for (int j = 0; j < ati.Size(); j++)
         {
-            weight = ati[j]->GetSurfaceFittingWeight();
-            weights.Append(weight);
+            if (ati[j]->IsSurfaceFittingEnabled())
+            {
+               weight = ati[j]->GetSurfaceFittingWeight();
+               weights.Append(weight);
+            }
         }
      }
   }
@@ -716,6 +727,39 @@ void TMOPNewtonSolver::GetSurfaceFittingError(const Vector &x_loc,
   }
 }

+bool TMOPNewtonSolver::IsSurfaceFittingEnabled() const
+{
+   const NonlinearForm *nlf = dynamic_cast<const NonlinearForm *>(oper);
+   const Array<NonlinearFormIntegrator*> &integs = *nlf->GetDNFI();
+   TMOP_Integrator *ti  = NULL;
+   TMOPComboIntegrator *co = NULL;
+
+   for (int i = 0; i < integs.Size(); i++)
+   {
+      ti = dynamic_cast<TMOP_Integrator *>(integs[i]);
+      if (ti)
+      {
+         if (ti->IsSurfaceFittingEnabled())
+         {
+            return true;
+         }
+      }
+      co = dynamic_cast<TMOPComboIntegrator *>(integs[i]);
+      if (co)
+      {
+         Array<TMOP_Integrator *> ati = co->GetTMOPIntegrators();
+         for (int j = 0; j < ati.Size(); j++)
+         {
+            if (ati[j]->IsSurfaceFittingEnabled())
+            {
+               return true;
+            }
+         }
+      }
+   }
+   return false;
+}
+
 void TMOPNewtonSolver::ProcessNewState(const Vector &x) const
 {
   const NonlinearForm *nlf = dynamic_cast<const NonlinearForm *>(oper);
@@ -801,38 +845,46 @@ void TMOPNewtonSolver::ProcessNewState(const Vector &x) const
   // adaptive surface fitting is enabled. The idea is to increase the
   // coefficient if the surface fitting error does not sufficiently
   // decrease between subsequent TMOPNewtonSolver iterations.
-   if (update_surf_fit_coeff)
+   if (surf_fit_coeff_update)
   {
      // Get surface fitting errors.
-      GetSurfaceFittingError(x_loc, surf_fit_err_avg, surf_fit_err_max);
+      GetSurfaceFittingError(x_loc, surf_fit_avg_err, surf_fit_max_err);
      // Get array with surface fitting weights.
-      Array<real_t> weights;
-      GetSurfaceFittingWeight(weights);
+      Array<real_t> fitweights;
+      GetSurfaceFittingWeight(fitweights);

      if (print_options.iterations)
      {
         mfem::out << "Avg/Max surface fitting error: " <<
-                   surf_fit_err_avg << " " <<
-                   surf_fit_err_max << "\n";
+                   surf_fit_avg_err << " " <<
+                   surf_fit_max_err << "\n";
         mfem::out << "Min/Max surface fitting weight: " <<
-                   weights.Min() << " " << weights.Max() << "\n";
+                   fitweights.Min() << " " << fitweights.Max() << "\n";
      }

-      real_t change_surf_fit_err = surf_fit_err_avg_prvs-surf_fit_err_avg;
-      real_t rel_change_surf_fit_err = change_surf_fit_err/surf_fit_err_avg_prvs;
+      real_t change_surf_fit_err = surf_fit_avg_err_prvs-surf_fit_avg_err;
+      real_t rel_change_surf_fit_err = change_surf_fit_err/surf_fit_avg_err_prvs;
+
      // Increase the surface fitting coefficient if the surface fitting error
-      // does not decrease sufficiently.
-      if (rel_change_surf_fit_err < surf_fit_rel_change_threshold)
+      // does not decrease sufficiently. If we are converging based on residual,
+      // also make sure we have not reached the maximum fitting weight and
+      // error threshold.
+      if (rel_change_surf_fit_err < surf_fit_err_rel_change_limit &&
+          (surf_fit_converge_error ||
+           (fitweights.Max() < surf_fit_weight_limit &&
+            surf_fit_max_err > surf_fit_max_err_limit)))
      {
-         UpdateSurfaceFittingWeight(surf_fit_scale_factor);
-         adapt_inc_count += 1;
+         real_t scale_factor = std::min(surf_fit_scale_factor,
+                                        surf_fit_weight_limit/fitweights.Max());
+         UpdateSurfaceFittingWeight(scale_factor);
+         surf_fit_adapt_count += 1;
      }
      else
      {
-         adapt_inc_count = 0;
+         surf_fit_adapt_count = 0;
      }
-      surf_fit_err_avg_prvs = surf_fit_err_avg;
-      update_surf_fit_coeff = false;
+      surf_fit_avg_err_prvs = surf_fit_avg_err;
+      surf_fit_coeff_update = false;
   }
 }

@@ -134,18 +134,20 @@ protected:
   int solver_type;
   bool parallel;

-   // Line search step is rejected if min(detJ) <= min_detJ_threshold.
-   real_t min_detJ_threshold = 0.0;
+   // Line search step is rejected if min(detJ) <= min_detJ_limit.
+   real_t min_detJ_limit = 0.0;

   // Surface fitting variables.
-   mutable real_t surf_fit_err_avg_prvs = 10000.0;
-   mutable real_t surf_fit_err_avg, surf_fit_err_max;
-   mutable bool update_surf_fit_coeff = false;
-   real_t surf_fit_max_threshold = -1.0;
-   real_t surf_fit_rel_change_threshold = 0.001;
+   mutable real_t surf_fit_avg_err_prvs = 10000.0;
+   mutable real_t surf_fit_avg_err, surf_fit_max_err;
+   mutable bool surf_fit_coeff_update = false;
+   real_t surf_fit_max_err_limit = -1.0;
+   real_t surf_fit_err_rel_change_limit = 0.001;
   real_t surf_fit_scale_factor = 0.0;
-   mutable int adapt_inc_count = 0;
-   mutable int max_adapt_inc_count = 10;
+   mutable int surf_fit_adapt_count = 0;
+   mutable int surf_fit_adapt_count_limit = 10;
+   mutable real_t surf_fit_weight_limit = 1e10;
+   bool surf_fit_converge_error = false;

   // Minimum determinant over the whole mesh. Used for mesh untangling.
   real_t *min_det_ptr = nullptr;
@@ -191,6 +193,9 @@ protected:
   void GetSurfaceFittingWeight(Array<real_t> &weights) const;
   ///@}

+   /// Check if surface fitting is enabled.
+   bool IsSurfaceFittingEnabled() const;
+
 public:
 #ifdef MFEM_USE_MPI
   TMOPNewtonSolver(MPI_Comm comm, const IntegrationRule &irule, int type = 0)
@@ -224,38 +229,94 @@ public:
   /// (ii) surface fitting weight.
   virtual void ProcessNewState(const Vector &x) const;

-   /** @name Methods for adaptive surface fitting weight. (Experimental) */
-   /// Enable/Disable adaptive surface fitting weight.
-   /// The weight is modified after each TMOPNewtonSolver iteration as:
-   /// w_{k+1} = w_{k} * @a surf_fit_scale_factor if relative change in
-   /// max surface fitting error < @a surf_fit_rel_change_threshold.
-   /// The solver terminates if the maximum surface fitting error does
-   /// not sufficiently decrease for @a max_adapt_inc_count consecutive
-   /// solver iterations or if the max error falls below @a surf_fit_max_threshold.
-   void EnableAdaptiveSurfaceFitting()
-   {
-      surf_fit_scale_factor = 10.0;
-      surf_fit_rel_change_threshold = 0.001;
-   }
+   /** @name Methods for adaptive surface fitting.
+       \brief These methods control the behavior of the weight and the
+       termination of the solver. (Experimental)
+
+       Adaptive fitting weight: The weight is modified after each
+       TMOPNewtonSolver iteration as:
+       w_{k+1} = w_{k} * \ref surf_fit_scale_factor if the relative
+       change in average fitting error < \ref surf_fit_err_rel_change_limit.
+       When converging based on the residual, we enforce the fitting weight
+       to be at-most \ref surf_fit_weight_limit, and increase it only if the
+       fitting error is below user prescribed threshold
+       (\ref surf_fit_max_err_limit).
+       See \ref SetAdaptiveSurfaceFittingScalingFactor and
+       \ref SetAdaptiveSurfaceFittingRelativeChangeThreshold.
+
+       Note that the solver stops if the maximum surface fitting error
+       does not sufficiently decrease for \ref surf_fit_adapt_count_limit (default 10)
+       consecutive increments of the fitting weight during weight adaptation.
+       This typically occurs when the mesh cannot align with the level-set
+       without degrading element quality.
+       See \ref SetMaxNumberofIncrementsForAdaptiveFitting.
+
+       Convergence criterion: There are two modes, residual- and error-based,
+       which can be toggled using \ref SetSurfaceFittingConvergenceBasedOnError.
+
+       (i) Residual based (default): Stop when the norm of the gradient of the
+       TMOP objective reaches the prescribed tolerance. This method is best used
+       with a reasonable value for \ref surf_fit_weight_limit when the
+       adaptive surface fitting scheme is used. See method
+       \ref SetSurfaceFittingWeightLimit.
+
+       (ii) Error based: Stop when the maximum fitting error
+       reaches the user-prescribed threshold, \ref surf_fit_max_err_limit.
+       In this case, \ref surf_fit_weight_limit is ignored during weight
+       adaptation.
+   */
+   ///@{
   void SetAdaptiveSurfaceFittingScalingFactor(real_t factor)
   {
+      MFEM_VERIFY(factor > 1.0, "Scaling factor must be greater than 1.");
      surf_fit_scale_factor = factor;
   }
   void SetAdaptiveSurfaceFittingRelativeChangeThreshold(real_t threshold)
   {
-      surf_fit_rel_change_threshold = threshold;
+      surf_fit_err_rel_change_limit = threshold;
   }
+   /// Used for stopping based on the number of consecutive failed weight
+   /// adaptation iterations.
+   // TODO: Rename to SetMaxNumberofIncrementsForAdaptiveSurfaceFitting
+   // in future.
   void SetMaxNumberofIncrementsForAdaptiveFitting(int count)
   {
-      max_adapt_inc_count = count;
+      surf_fit_adapt_count_limit = count;
   }
+   /// Used for error-based surface fitting termination.
   void SetTerminationWithMaxSurfaceFittingError(real_t max_error)
   {
-      surf_fit_max_threshold = max_error;
+      surf_fit_max_err_limit = max_error;
+      surf_fit_converge_error = true;
   }
+   /// Could be used with both error-based or residual-based convergence.
+   void SetSurfaceFittingMaxErrorLimit(real_t max_error)
+   {
+      surf_fit_max_err_limit = max_error;
+   }
+   /// Used for residual-based surface fitting termination.
+   void SetSurfaceFittingWeightLimit(real_t weight)
+   {
+      surf_fit_weight_limit = weight;
+   }
+   /// Toggle convergence based on residual or error.
+   void SetSurfaceFittingConvergenceBasedOnError(bool mode)
+   {
+      surf_fit_converge_error = mode;
+      if (surf_fit_converge_error)
+      {
+         MFEM_VERIFY(surf_fit_max_err_limit >= 0,
+                     "Fitting error based convergence requires the user to "
+                     "first set the error threshold."
+                     "See SetTerminationWithMaxSurfaceFittingError");
+      }
+   }
+   ///@}
+
+   /// Set minimum determinant enforced during line-search.
   void SetMinimumDeterminantThreshold(real_t threshold)
   {
-      min_detJ_threshold = threshold;
+      min_detJ_limit = threshold;
   }

   virtual void Mult(const Vector &b, Vector &x) const
@@ -124,7 +124,7 @@ T Array<T>::Sum()
 }

 template <class T>
-int Array<T>::IsSorted()
+int Array<T>::IsSorted() const
 {
   T val_prev = operator[](0), val;
   for (int i = 1; i < size; i++)
@@ -74,10 +74,14 @@ public:
   inline Array(int asize, MemoryType mt)
      : size(asize) { asize > 0 ? data.New(asize, mt) : data.Reset(mt); }

-   /** @brief Creates array using an externally allocated pointer @a data_ to
-       @a asize elements. The data pointer will not be deleted by Array. */
-   inline Array(T *data_, int asize)
-   { data.Wrap(data_, asize, false); size = asize; }
+   /** @brief Creates array using an externally allocated host pointer @a data_
+       to @a asize elements. If @a own_data is true, the array takes ownership
+       of the pointer.
+
+       When @a own_data is true, the pointer @a data_ must be allocated with
+       MemoryType given by MemoryManager::GetHostMemoryType(). */
+   inline Array(T *data_, int asize, bool own_data = false)
+   { data.Wrap(data_, asize, own_data); size = asize; }

   /// Copy constructor: deep copy from @a src
   /** This method supports source arrays using any MemoryType. */
@@ -205,7 +209,14 @@ public:
   inline void Copy(Array &copy) const;

   /// Make this Array a reference to a pointer.
-   inline void MakeRef(T *, int);
+   /** When @a own_data is true, the pointer @a data_ must be allocated with
+       MemoryType given by MemoryManager::GetHostMemoryType(). */
+   inline void MakeRef(T *data_, int size_, bool own_data = false);
+
+   /// Make this Array a reference to a pointer.
+   /** When @a own_data is true, the pointer @a data_ must be allocated with
+       MemoryType given by @a mt. */
+   inline void MakeRef(T *data_, int size, MemoryType mt, bool own_data);

   /// Make this Array a reference to 'master'.
   inline void MakeRef(const Array &master);
@@ -262,7 +273,7 @@ public:
   }

   /// Return 1 if the array is sorted from lowest to highest.  Otherwise return 0.
-   int IsSorted();
+   int IsSorted() const;

   /// Fill the entries of the array with the cumulative sum of the entries.
   void PartialSum();
@@ -868,11 +879,19 @@ inline void Array<T>::Copy(Array &copy) const
 }

 template <class T>
-inline void Array<T>::MakeRef(T *p, int s)
+inline void Array<T>::MakeRef(T *data_, int size_, bool own_data)
 {
   data.Delete();
-   data.Wrap(p, s, false);
-   size = s;
+   data.Wrap(data_, size_, own_data);
+   size = size_;
+}
+
+template <class T>
+inline void Array<T>::MakeRef(T *data_, int size_, MemoryType mt, bool own_data)
+{
+   data.Delete();
+   data.Wrap(data_, size_, mt, own_data);
+   size = size_;
 }

 template <class T>
@@ -288,4 +288,3 @@ void ArraysByName<T>::Load(std::istream &in)
 }

 #endif
-
@@ -275,7 +275,7 @@ void GroupTopology::Save(ostream &os) const
   os << "\ncommunication_groups\n";
   os << "number_of_groups " << NGroups() << "\n\n";

-   os << "# number of entities in each group, followed by group ids in group\n";
+   os << "# number of entities in each group, followed by ranks in group\n";
   for (int group_id = 0; group_id < NGroups(); ++group_id)
   {
      int group_size = GetGroupSize(group_id);
@@ -14,6 +14,9 @@
 #ifdef MFEM_USE_CEED
 #include "../fem/ceed/interface/util.hpp"
 #endif
+#ifdef MFEM_USE_MPI
+#include "../linalg/hypre.hpp"
+#endif

 #include <unordered_map>
 #include <string>
@@ -250,6 +253,10 @@ void Device::Configure(const std::string &device, const int device_id)

   // Only '*this' will call the MemoryManager::Destroy() method.
   destroy_mm = true;
+
+#ifdef MFEM_USE_MPI
+   Hypre::InitDevice();
+#endif
 }

 // static method
@@ -19,6 +19,9 @@
 #include "device.hpp"
 #include "mem_manager.hpp"
 #include "../linalg/dtensor.hpp"
+#ifdef MFEM_USE_MPI
+#include <_hypre_utilities.h>
+#endif

 namespace mfem
 {
@@ -780,6 +783,63 @@ inline void forall_3D_grid(int N, int X, int Y, int Z, int G, lambda &&body)
   ForallWrap<3>(true, N, body, X, Y, Z, G);
 }

+#ifdef MFEM_USE_MPI
+
+// Function mfem::hypre_forall_cpu() similar to mfem::forall, but it always
+// executes on the CPU using sequential or OpenMP-parallel execution based on
+// the hypre build time configuration.
+template<typename lambda>
+inline void hypre_forall_cpu(int N, lambda &&body)
+{
+#ifdef HYPRE_USING_OPENMP
+   #pragma omp parallel for HYPRE_SMP_SCHEDULE
+#endif
+   for (int i = 0; i < N; i++) { body(i); }
+}
+
+// Function mfem::hypre_forall_gpu() similar to mfem::forall, but it always
+// executes on the GPU device that hypre was configured with at build time.
+#if defined(HYPRE_USING_GPU)
+template<typename lambda>
+inline void hypre_forall_gpu(int N, lambda &&body)
+{
+#if defined(HYPRE_USING_CUDA)
+   CuWrap1D(N, body);
+#elif defined(HYPRE_USING_HIP)
+   HipWrap1D(N, body);
+#else
+#error Unknown HYPRE GPU backend!
+#endif
+}
+#endif
+
+// Function mfem::hypre_forall() similar to mfem::forall, but it executes on the
+// device, CPU or GPU, that hypre was configured with at build time (when the
+// HYPRE version is < 2.31.0) or at runtime (when HYPRE was configured with GPU
+// support at build time and HYPRE's version is >= 2.31.0). This selection is
+// generally independent of what device was selected in MFEM's runtime
+// configuration.
+template<typename lambda>
+inline void hypre_forall(int N, lambda &&body)
+{
+#if !defined(HYPRE_USING_GPU)
+   hypre_forall_cpu(N, body);
+#elif MFEM_HYPRE_VERSION < 23100
+   hypre_forall_gpu(N, body);
+#else // HYPRE_USING_GPU is defined and MFEM_HYPRE_VERSION >= 23100
+   if (!HypreUsingGPU())
+   {
+      hypre_forall_cpu(N, body);
+   }
+   else
+   {
+      hypre_forall_gpu(N, body);
+   }
+#endif
+}
+
+#endif // MFEM_USE_MPI
+
 } // namespace mfem

 #endif // MFEM_FORALL_HPP
@@ -1154,6 +1154,10 @@ void MemoryManager::Copy_(void *dst_h_ptr, const void *src_h_ptr,
   //  dest   d  | h2d   d2d   d2d
   //        hd  | h2h   d2d   d2d

+   MFEM_ASSERT(bytes != 0, "this method should not be called with bytes = 0");
+   MFEM_ASSERT(dst_h_ptr != nullptr, "invalid dst_h_ptr = nullptr");
+   MFEM_ASSERT(src_h_ptr != nullptr, "invalid src_h_ptr = nullptr");
+
   const bool dst_on_host =
      (dst_flags & Mem::VALID_HOST) &&
      (!(dst_flags & Mem::VALID_DEVICE) ||
@@ -1229,6 +1233,10 @@ void MemoryManager::Copy_(void *dst_h_ptr, const void *src_h_ptr,
 void MemoryManager::CopyToHost_(void *dest_h_ptr, const void *src_h_ptr,
                                size_t bytes, unsigned src_flags)
 {
+   MFEM_ASSERT(bytes != 0, "this method should not be called with bytes = 0");
+   MFEM_ASSERT(dest_h_ptr != nullptr, "invalid dest_h_ptr = nullptr");
+   MFEM_ASSERT(src_h_ptr != nullptr, "invalid src_h_ptr = nullptr");
+
   const bool src_on_host = src_flags & Mem::VALID_HOST;
   if (src_on_host)
   {
@@ -1255,6 +1263,10 @@ void MemoryManager::CopyToHost_(void *dest_h_ptr, const void *src_h_ptr,
 void MemoryManager::CopyFromHost_(void *dest_h_ptr, const void *src_h_ptr,
                                  size_t bytes, unsigned &dest_flags)
 {
+   MFEM_ASSERT(bytes != 0, "this method should not be called with bytes = 0");
+   MFEM_ASSERT(dest_h_ptr != nullptr, "invalid dest_h_ptr = nullptr");
+   MFEM_ASSERT(src_h_ptr != nullptr, "invalid src_h_ptr = nullptr");
+
   const bool dest_on_host = dest_flags & Mem::VALID_HOST;
   if (dest_on_host)
   {
@@ -18,8 +18,14 @@
 #include <cstring> // std::memcpy
 #include <type_traits> // std::is_const
 #include <cstddef> // std::max_align_t
+
 #ifdef MFEM_USE_MPI
-#include <HYPRE_config.h> // HYPRE_USING_GPU
+// Enable internal hypre timing routines
+#define HYPRE_TIMING
+#include <HYPRE_utilities.h> // for HYPRE_GetMemoryLocation() and others
+#if (21400 <= MFEM_HYPRE_VERSION) && (MFEM_HYPRE_VERSION < 21900)
+#include <_hypre_utilities.h> // for HYPRE_MEMORY_HOST and others
+#endif
 #endif

 namespace mfem
@@ -869,6 +875,45 @@ public:
 };


+#ifdef MFEM_USE_MPI
+
+#if MFEM_HYPRE_VERSION < 21400
+#define HYPRE_MEMORY_DEVICE (0)
+#define HYPRE_MEMORY_HOST   (1)
+#endif
+#if MFEM_HYPRE_VERSION < 21900
+typedef int HYPRE_MemoryLocation;
+#endif
+
+/// Return the configured HYPRE_MemoryLocation
+inline HYPRE_MemoryLocation GetHypreMemoryLocation()
+{
+#if !defined(HYPRE_USING_GPU)
+   return HYPRE_MEMORY_HOST;
+#elif MFEM_HYPRE_VERSION < 23100
+   return HYPRE_MEMORY_DEVICE;
+#else // HYPRE_USING_GPU is defined and MFEM_HYPRE_VERSION >= 23100
+   HYPRE_MemoryLocation loc;
+   HYPRE_GetMemoryLocation(&loc);
+   return loc;
+#endif
+}
+
+/// Return true if HYPRE is configured to use GPU
+inline bool HypreUsingGPU()
+{
+#if !defined(HYPRE_USING_GPU)
+   return false;
+#elif MFEM_HYPRE_VERSION < 23100
+   return true;
+#else // HYPRE_USING_GPU is defined and MFEM_HYPRE_VERSION >= 23100
+   return GetHypreMemoryLocation() != HYPRE_MEMORY_HOST;
+#endif
+}
+
+#endif // MFEM_USE_MPI
+
+
 // Inline methods

 template <typename T>
@@ -1004,10 +1049,12 @@ inline void Memory<T>::MakeAlias(const Memory &base, int offset, int size)
         // If the following condition is true then MemoryManager::Exists()
         // should also be true:
         IsDeviceMemory(MemoryManager::GetDeviceMemoryType())
-#else
-         // When HYPRE_USING_GPU is defined we always register the 'base' if
-         // the MemoryManager::Exists():
+#elif MFEM_HYPRE_VERSION < 23100
+         // When HYPRE_USING_GPU is defined and HYPRE < 2.31.0, we always
+         // register the 'base' if the MemoryManager::Exists():
         MemoryManager::Exists()
+#else // HYPRE_USING_GPU is defined and MFEM_HYPRE_VERSION >= 23100
+         MemoryManager::Exists() && HypreUsingGPU()
 #endif
      )
      {
@@ -1213,9 +1260,10 @@ template <typename T>
 inline void Memory<T>::CopyFrom(const Memory &src, int size)
 {
   MFEM_VERIFY(src.capacity>=size && capacity>=size, "Incorrect size");
+   if (size <= 0) { return; }
   if (!(flags & Registered) && !(src.flags & Registered))
   {
-      if (h_ptr != src.h_ptr && size != 0)
+      if (h_ptr != src.h_ptr)
      {
         MFEM_ASSERT(h_ptr + size <= src.h_ptr || src.h_ptr + size <= h_ptr,
                     "data overlaps!");
@@ -1233,9 +1281,10 @@ template <typename T>
 inline void Memory<T>::CopyFromHost(const T *src, int size)
 {
   MFEM_VERIFY(capacity>=size, "Incorrect size");
+   if (size <= 0) { return; }
   if (!(flags & Registered))
   {
-      if (h_ptr != src && size != 0)
+      if (h_ptr != src)
      {
         MFEM_ASSERT(h_ptr + size <= src || src + size <= h_ptr,
                     "data overlaps!");
@@ -1252,7 +1301,6 @@ inline void Memory<T>::CopyFromHost(const T *src, int size)
 template <typename T>
 inline void Memory<T>::CopyTo(Memory &dest, int size) const
 {
-   MFEM_VERIFY(capacity>=size, "Incorrect size");
   dest.CopyFrom(*this, size);
 }

@@ -1260,9 +1308,10 @@ template <typename T>
 inline void Memory<T>::CopyToHost(T *dest, int size) const
 {
   MFEM_VERIFY(capacity>=size, "Incorrect size");
+   if (size <= 0) { return; }
   if (!(flags & Registered))
   {
-      if (h_ptr != dest && size != 0)
+      if (h_ptr != dest)
      {
         MFEM_ASSERT(h_ptr + size <= dest || dest + size <= h_ptr,
                     "data overlaps!");
@@ -134,7 +134,7 @@ int socketbuf::open(const char hostname[], int port)
      {
         closesocket(socket_descriptor);
         socket_descriptor = -2;
-         return -1;
+         continue;
      }
 #endif

@@ -148,7 +148,7 @@ int socketbuf::open(const char hostname[], int port)
   }

   freeaddrinfo(res);
-   return 0;
+   return (socket_descriptor < 0) ? -1 : 0;
 }

 int socketbuf::close()
@@ -207,7 +207,12 @@ template <> inline void Swap<Table>(Table &a, Table &b)
 void Transpose (const Table &A, Table &At, int ncols_A_ = -1);
 Table * Transpose (const Table &A);

-///  Transpose an Array<int>
+///  @brief Transpose an Array<int>.
+///
+/// The array @a A represents a table where each row @a i has exactly one
+/// connection to the column (TYPE II) index specified by @a A[i].
+///
+/// @note The column (TYPE II) indices in each row of @a At will be sorted.
 void Transpose(const Array<int> &A, Table &At, int ncols_A_ = -1);

 ///  C = A * B  (as boolean matrices)
@@ -400,6 +400,9 @@ inline double StopWatch::SystTime()

 StopWatch::StopWatch() : M(new internal::StopWatch) { }

+StopWatch::StopWatch(const StopWatch &sw)
+   : M(new internal::StopWatch(*(sw.M))) { }
+
 void StopWatch::Clear()
 {
   M->Clear();
@@ -40,6 +40,7 @@ private:
 public:
   /// Creates a new (stopped) StopWatch object.
   StopWatch();
+   StopWatch(const StopWatch &);

   /// Clear the elapsed time on the stopwatch and restart it if it's running.
   void Clear();
@@ -23,6 +23,10 @@
 #include "amgxsolver.hpp"
 #ifdef MFEM_USE_AMGX

+#ifdef MFEM_USE_MPI
+#include "../general/communication.hpp"
+#endif
+
 namespace mfem
 {

@@ -81,7 +81,7 @@ void BlockOperator::Mult (const Vector & x, Vector & y) const
      tmp.SetSize(row_offsets[iRow+1] - row_offsets[iRow]);
      for (int jCol=0; jCol < nColBlocks; ++jCol)
      {
-         if (op(iRow,jCol))
+         if (op(iRow,jCol) && coef(iRow,jCol) != 0.)
         {
            op(iRow,jCol)->Mult(xblock.GetBlock(jCol), tmp);
            yblock.GetBlock(iRow).Add(coef(iRow,jCol), tmp);
@@ -112,7 +112,7 @@ void BlockOperator::MultTranspose (const Vector & x, Vector & y) const
      tmp.SetSize(col_offsets[iRow+1] - col_offsets[iRow]);
      for (int jCol=0; jCol < nRowBlocks; ++jCol)
      {
-         if (op(jCol,iRow))
+         if (op(jCol,iRow) && coef(jCol,iRow) != 0.)
         {
            op(jCol,iRow)->MultTranspose(xblock.GetBlock(jCol), tmp);
            yblock.GetBlock(iRow).Add(coef(jCol,iRow), tmp);
@@ -1,3 +1,14 @@
+// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
 #include "cpardiso.hpp"
 #include "hypre.hpp"
 #include <algorithm>
@@ -532,6 +532,69 @@ MatrixInverse *DenseMatrix::Inverse() const
   return new DenseMatrixInverse(*this);
 }

+void DenseMatrix::Exponential()
+{
+   MFEM_ASSERT(Height() == Width() && Height() <= 2,
+               "The matrix must be square and "
+               << "of size less than or equal to 2."
+               << "  Height() = " << Height()
+               << ", Width() = " << Width());
+
+   switch (Height())
+   {
+      case 1:
+      {
+         data[0] = std::exp(data[0]);
+         break;
+      }
+      case 2:
+      {
+         /// Formulas from Corollary 2.4 of doi:10.1109/9.233156
+         /// Note typo in the paper, in the prefactor in the equation under (i).
+         const real_t a = data[0];
+         const real_t b = data[1];
+         const real_t c = data[2];
+         const real_t d = data[3];
+         const real_t e = (a - d)*(a - d) + 4*b*c;
+         const real_t f = std::exp((a + d)/2.0);
+         const real_t g = std::sqrt(std::abs(e)) / 2.0;
+
+         if (e == 0)
+         {
+            data[0] = 1.0 + (a - d)/2.0;
+            data[3] = 1.0 - (a - d)/2.0;
+         }
+         else if (e > 0)
+         {
+            data[0] = std::cosh(g) + (a - d)/2 * std::sinh(g) / g;
+            data[1] = b * std::sinh(g) / g;
+            data[2] = c * std::sinh(g) / g;
+            data[3] = std::cosh(g) - (a - d)/2 * std::sinh(g) / g;
+         }
+         else
+         {
+            data[0] = std::cos(g) + (a - d)/2 * std::sin(g) / g;
+            data[1] = b * std::sin(g) / g;
+            data[2] = c * std::sin(g) / g;
+            data[3] = std::cos(g) - (a - d)/2 * std::sin(g) / g;
+         }
+         for (int i = 0; i < 4; i++)
+         {
+            data[i] *= f;
+         }
+         break;
+      }
+      case 3:
+      {
+         MFEM_ABORT("3x3 matrices are not currently supported");
+      }
+      default:
+      {
+         MFEM_ABORT("Only 1x1 and 2x2 matrices are currently supported");
+      }
+   }
+}
+
 real_t DenseMatrix::Det() const
 {
   MFEM_ASSERT(Height() == Width() && Height() > 0,
@@ -3217,6 +3280,93 @@ void MultAtB(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &AtB)
 #endif
 }

+void AddMultAtB(const DenseMatrix &A, const DenseMatrix &B,
+                DenseMatrix &AtB)
+{
+   MFEM_ASSERT(AtB.Height() == A.Width() && AtB.Width() == B.Width() &&
+               A.Height() == B.Height(), "incompatible dimensions");
+
+#ifdef MFEM_USE_LAPACK
+   static char transa = 'T', transb = 'N';
+   static real_t alpha = 1.0, beta = 1.0;
+   int m = A.Width(), n = B.Width(), k = A.Height();
+
+#ifdef MFEM_USE_SINGLE
+   sgemm_(&transa, &transb, &m, &n, &k, &alpha, A.Data(), &k,
+#elif defined MFEM_USE_DOUBLE
+   dgemm_(&transa, &transb, &m, &n, &k, &alpha, A.Data(), &k,
+#endif
+          B.Data(), &k, &beta, AtB.Data(), &m);
+#else
+   const int ah = A.Height();
+   const int aw = A.Width();
+   const int bw = B.Width();
+   const real_t *ad = A.Data();
+   const real_t *bd = B.Data();
+   real_t *cd = AtB.Data();
+
+   for (int j = 0; j < bw; j++)
+   {
+      const real_t *ap = ad;
+      for (int i = 0; i < aw; i++)
+      {
+         real_t d = 0.0;
+         for (int k = 0; k < ah; k++)
+         {
+            d += ap[k] * bd[k];
+         }
+         *(cd++) += d;
+         ap += ah;
+      }
+      bd += ah;
+   }
+#endif
+}
+
+void AddMult_a_AtB(real_t a, const DenseMatrix &A, const DenseMatrix &B,
+                   DenseMatrix &AtB)
+{
+   MFEM_ASSERT(AtB.Height() == A.Width() && AtB.Width() == B.Width() &&
+               A.Height() == B.Height(), "incompatible dimensions");
+
+#ifdef MFEM_USE_LAPACK
+   static char transa = 'T', transb = 'N';
+   real_t alpha = a;
+   static real_t beta = 1.0;
+   int m = A.Width(), n = B.Width(), k = A.Height();
+
+#ifdef MFEM_USE_SINGLE
+   sgemm_(&transa, &transb, &m, &n, &k, &alpha, A.Data(), &k,
+#elif defined MFEM_USE_DOUBLE
+   dgemm_(&transa, &transb, &m, &n, &k, &alpha, A.Data(), &k,
+#endif
+          B.Data(), &k, &beta, AtB.Data(), &m);
+#else
+   const int ah = A.Height();
+   const int aw = A.Width();
+   const int bw = B.Width();
+   const real_t *ad = A.Data();
+   const real_t *bd = B.Data();
+   real_t *cd = AtB.Data();
+
+   for (int j = 0; j < bw; j++)
+   {
+      const real_t *ap = ad;
+      for (int i = 0; i < aw; i++)
+      {
+         real_t d = 0.0;
+         for (int k = 0; k < ah; k++)
+         {
+            d += ap[k] * bd[k];
+         }
+         *(cd++) += a * d;
+         ap += ah;
+      }
+      bd += ah;
+   }
+#endif
+}
+
 void AddMult_a_AAt(real_t a, const DenseMatrix &A, DenseMatrix &AAt)
 {
   real_t d;
@@ -207,6 +207,10 @@ public:
   /// Replaces the current matrix with its square root inverse
   void SquareRootInverse();

+   /// Replaces the current matrix with its exponential
+   /// (currently only supports 2x2 matrices)
+   void Exponential();
+
   /// Calculates the determinant of the matrix
   /// (optimized for 2x2, 3x3, and 4x4 matrices)
   real_t Det() const;
@@ -580,6 +584,13 @@ void AddMult_a_ABt(real_t a, const DenseMatrix &A, const DenseMatrix &B,
 /// Multiply the transpose of a matrix A with a matrix B:   At*B
 void MultAtB(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &AtB);

+/// AtB += A^t * B
+void AddMultAtB(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &AtB);
+
+/// AtB += a * A^t * B
+void AddMult_a_AtB(real_t a, const DenseMatrix &A, const DenseMatrix &B,
+                   DenseMatrix &AtB);
+
 /// AAt += a * A * A^t
 void AddMult_a_AAt(real_t a, const DenseMatrix &A, DenseMatrix &AAt);

--- a/Show More
+++ b/Show More