style

Updates CHANGELOG
More fixes for building contact mini-app with make and cmake
2026-04-14 08:54:23 -07:00 · 2026-04-13 21:52:37 -07:00 · 2026-04-13 21:47:15 -07:00 · 2026-04-13 18:54:17 -07:00 · 2026-04-13 12:40:24 -07:00 · 2026-04-13 12:37:18 -07:00
240 changed files with 16794 additions and 5443 deletions
@@ -25,7 +25,7 @@ runs:
  steps:
    - uses: ./.github/actions/sanitize/config

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

@@ -17,11 +17,11 @@ jobs:
    runs-on: ubuntu-latest
    name: 2.19.0
    steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v6
      - uses: ./.github/actions/sanitize/config
      - name: Cache
        id: cache
-        uses: actions/cache@v4
+        uses: actions/cache@v5
        with:
          path: ${{env.HYPRE_DIR}}
          key: ${{runner.os}}-ompi-build-${{env.HYPRE_DIR}}-int32-fp64-v2.5
@@ -27,13 +27,13 @@ jobs:
            llvm_use_sanitizer: "Undefined"
    name: ${{matrix.sanitizer}}
    steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v6
      - uses: ./.github/actions/sanitize/config
        with:
          NO_FLAGS: true
      - name: Cache
        id: cache
-        uses: actions/cache@v4
+        uses: actions/cache@v5
        with:
          path: ${{env.LLVM_DIR}}
          key: build-libcxx-${{env.LLVM_VER}}-${{matrix.sanitizer}}
@@ -17,11 +17,11 @@ jobs:
    runs-on: ubuntu-latest
    name: lsan.supp
    steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v6
      - uses: ./.github/actions/sanitize/config
      - name: Cache
        id: cache
-        uses: actions/cache@v4
+        uses: actions/cache@v5
        with:
          path: ${{env.LSAN_DIR}}
          key: build-lsan-suppression-file
@@ -17,11 +17,11 @@ jobs:
    runs-on: ubuntu-latest
    name: 4.0.3
    steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v6
      - uses: ./.github/actions/sanitize/config
      - name: Cache
        id: cache
-        uses: actions/cache@v4
+        uses: actions/cache@v5
        with:
          path: ${{env.METIS_DIR}}
          key: ${{runner.os}}-build-${{env.METIS_DIR}}-v2.5
@@ -28,7 +28,7 @@ jobs:
  build:
    runs-on: ubuntu-latest
    steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v6
      - uses: ./.github/actions/sanitize/mfem
        with:
          par: ${{inputs.par}}
@@ -40,7 +40,7 @@ jobs:
    env:
      ex: ${{inputs.par && 'ex1p' || 'ex1'}}
    steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v6
      - uses: ./.github/actions/sanitize/restore
        id: restore
        with:
@@ -58,7 +58,7 @@ jobs:
    env:
      exclude: ${{inputs.par && '-E "_ser"' || ''}}
    steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v6
      - uses: ./.github/actions/sanitize/restore
        id: restore
        with:
@@ -82,7 +82,7 @@ jobs:
    env:
      exclude: ${{inputs.par && '-E "_ser"' || ''}}
    steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v6
      - uses: ./.github/actions/sanitize/restore
        id: restore
        with:
@@ -107,7 +107,7 @@ jobs:
      run: ${{inputs.par && '-R "_cpu_np"' || ''}}
      exclude: ${{inputs.par && '"unit_tests|debug"' || '"^unit_tests$|debug"'}}
    steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v6
      - uses: ./.github/actions/sanitize/restore
        id: restore
        with:
@@ -131,7 +131,7 @@ jobs:
    env:
      unit_tests: ${{inputs.par && 'punit_tests' || 'unit_tests'}}
    steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v6
      - uses: ./.github/actions/sanitize/restore
        id: restore
        with:
@@ -165,7 +165,7 @@ jobs:
      unit_tests: ${{inputs.par && 'punit_tests' || 'unit_tests'}}
      np: ${{inputs.par && '_np=2' || ''}}
    steps:
-      - uses: actions/checkout@v4
+      - uses: actions/checkout@v6
      - uses: ./.github/actions/sanitize/restore
        id: restore
        with:
@@ -313,6 +313,8 @@ miniapps/nurbs/nurbs_ex1
 miniapps/nurbs/nurbs_ex1p
 miniapps/nurbs/nurbs_ex3
 miniapps/nurbs/nurbs_ex5
+miniapps/nurbs/nurbs_ex10
+miniapps/nurbs/nurbs_ex10p
 miniapps/nurbs/nurbs_ex11p
 miniapps/nurbs/nurbs_ex24
 miniapps/nurbs/nurbs_solenoidal
@@ -338,7 +340,14 @@ miniapps/nurbs/nurbs_naca_cmesh
 miniapps/nurbs/naca-cmesh.mesh
 miniapps/nurbs/glvis_naca-cmesh.mesh
 miniapps/nurbs/Naca_cmesh
+miniapps/nurbs/nurbs_mesh_info
+miniapps/nurbs/k*_*.dat
 miniapps/nurbs/*-Surface.mesh
+miniapps/nurbs/*.mesh
+miniapps/nurbs/*.sol
+miniapps/nurbs/deformed.*
+miniapps/nurbs/elastic_energy.*
+miniapps/nurbs/velocity.*

 miniapps/performance/ex1
 miniapps/performance/ex1p
@@ -360,6 +369,7 @@ miniapps/shifted/lsf_integral
 miniapps/tools/display-basis
 miniapps/tools/load-dc
 miniapps/tools/convert-dc
+miniapps/tools/compare-dc
 miniapps/tools/gridfunction-bounds
 miniapps/tools/lor-transfer
 miniapps/tools/plor-transfer
@@ -8,9 +8,48 @@
                               https://mfem.org


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

+Discretization improvements
+---------------------------
+- Improved the gridfunction projection routines. Projections work for Scalar,
+  Vector and VectorFE, also NURBS versions. Optionally different types of
+  projections can be selected, default behaviour has not changed.
+
+- Added methods to estimate function extremum using piecewise linear bounds +
+  recursive subdivision.
+
+Meshing improvements
+--------------------
+- Improved support for 1D NURBS meshes with variable order, including using
+  the patches construct for 1D NURBS meshes.
+
+New and updated examples and miniapps
+-------------------------------------
+- Electromagnetics/lorentz miniapp has been updated to leverage the ParticleSet
+  capability.
+
+Miscellaneous
+-------------
+- Removes the SidreDataCollection class from MFEM in favor of the 
+  MFEMSidreDataCollection class in the Axom library (https://github.com/llnl/axom).
+

 Version 4.9, released on Dec 11, 2025
 =====================================
@@ -95,6 +134,23 @@ Linear and nonlinear solvers
  Filtering (AMGF), providing robust preconditioning for linear systems arising
  in constrained optimization problems such as frictionless contact.

+  Added 'GetResiduals' and 'GetFinalAbsResidualNorm' to 'HyprePCG',
+  'HypreGMRES', and 'HypreFGMRES' to get 'r' and '|r|_p'. Note that the latter
+  computes '|r|_p' from 'r' instead of returning a cached value like the
+  relative 'GetFinalResidualNorm'. These require Hypre >= 2.15.0.
+
+  Changed the default solver parameters for 'HyprePCG' to 'tol=1e-6' and
+  'max_iter=1000'. This matches the default parameters in Hypre 3.0.
+
+  Added various helper functions for querying/modifying Hypre solvers:
+  'HypreSmoother::GetType', 'HypreSmoother::GetSOROptions',
+  'HypreSmoother::GetPolyOptions', 'HypreSmoother::GetWindowParameters',
+  'HypreSmoother::IsOperatorSymmetric', 'HyprePCG::GetTol',
+  'HyprePCG::GetAbsTol', 'HyprePCG::GetMaxIter', 'HyprePCG::SetUseTwoNorm',
+  'HypreGMRES::GetTol', 'HypreGMRES::GetAbsTol', 'HypreGMRES::GetMaxIter',
+  'HypreGMRES::GetKDim', 'HypreFGMRES::GetTol', 'HypreFGMRES::GetMaxIter',
+  'HypreFGMRES::GetKDim', and 'HypreBoomerAMG::GetMaxIter'.
+
 GPU computing
 -------------
 - Added the 'gpu', 'raja-gpu', and 'ceed-gpu' backend aliases/shortcuts which
@@ -75,12 +75,10 @@ set(XSDK_ENABLE_Fortran OFF)

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

-# Axom/Sidre
-if (MFEM_USE_SIDRE)
-  find_package(Axom REQUIRED Axom)
-endif()
-
 # PUMI
 if (MFEM_USE_PUMI)
  # If PUMI_DIR was specified, only link to that directory,
@@ -629,7 +622,7 @@ find_package(Threads REQUIRED)
 #    integers, the METIS header (with 32-bit indices, as used by mfem) needs to
 #    be before SuiteSparse.
 set(MFEM_TPLS OPENMP HYPRE LAPACK BLAS SuperLUDist STRUMPACK METIS SuiteSparse
-    SUNDIALS PETSC SLEPC MUMPS AXOM FMS CONDUIT Ginkgo GNUTLS GSLIB HDF5
+    SUNDIALS PETSC SLEPC MUMPS FMS CONDUIT Ginkgo GNUTLS GSLIB HDF5
    NETCDF MPFR PUMI HIOP POSIXCLOCKS MFEMBacktrace ZLIB OCCA CEED RAJA UMPIRE
    ADIOS2 MKL_CPARDISO MKL_PARDISO AMGX MAGMA CUSPARSE CUBLAS CALIPER CODIPACK
    BENCHMARK PARELAG TRIBOL MPI_CXX HIP HIPBLAS HIPSPARSE MOONOLITH BLITZ
@@ -652,6 +645,8 @@ foreach(TPL IN LISTS MFEM_TPLS)
  endif()
 endforeach(TPL)

+# reverse to remove the first instance of entries in TPL_LIBRARIES
+# so later duplicates are kept (for dependency ordering)
 list(REVERSE TPL_LIBRARIES)
 list(REMOVE_DUPLICATES TPL_LIBRARIES)
 list(REVERSE TPL_LIBRARIES)
@@ -723,6 +718,7 @@ set(MFEM_INSTALL_DIR ${CMAKE_INSTALL_PREFIX})

 # Declaring the library
 mfem_add_library(mfem ${SOURCES} ${HEADERS} ${MASTER_HEADERS})
+target_compile_features(mfem PUBLIC cxx_std_${CMAKE_CXX_STANDARD})
 # message(STATUS "TPL_LIBRARIES = ${TPL_LIBRARIES}")
 target_link_libraries(mfem PUBLIC ${TPL_LIBRARIES} ${TPL_TARGETS})
 if (TPL_TARGETS)
@@ -869,11 +865,12 @@ add_dependencies(exec
 #  - https://cmake.org/Bug/view.php?id=8438

 # Add a target to copy the mfem data directory to the build directory
-add_custom_command(OUTPUT data_is_copied
-  COMMAND ${CMAKE_COMMAND} -E copy_directory ${PROJECT_SOURCE_DIR}/data data
-  COMMAND ${CMAKE_COMMAND} -E touch data_is_copied
-  COMMENT "Copying the data directory ...")
-add_custom_target(copy_data DEPENDS data_is_copied)
+# Implementable as a single copy_directory_if_different command w/ CMake >= 3.26
+file(GLOB DATA_FILES CONFIGURE_DEPENDS ${PROJECT_SOURCE_DIR}/data/*)
+add_custom_target(copy_data
+  COMMAND ${CMAKE_COMMAND} -E make_directory data
+  COMMAND ${CMAKE_COMMAND} -E copy_if_different ${DATA_FILES} data
+  COMMENT "Syncing the data directory ...")
 # Add 'copy_data' as a prerequisite for all executables, if the source and the
 # build directories are not the same.
 if (NOT ("${PROJECT_SOURCE_DIR}" STREQUAL "${PROJECT_BINARY_DIR}"))
@@ -1005,9 +1002,15 @@ install(FILES
 install(EXPORT ${PROJECT_NAME_UC}Targets
    DESTINATION ${INSTALL_CMAKE_DIR})

+# Install the data directory if present, i.e. if the copy_data target is built
+install(DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/data
+  DESTINATION ${MFEM_INSTALL_DIR} OPTIONAL)
+
 #-------------------------------------------------------------------------------
 # Create 'config.mk' from 'config.mk.in' for the build and install locations and
 # define install rules for 'config.mk' and 'test.mk'
 #-------------------------------------------------------------------------------
-
+if (MFEM_USE_CUDA OR MFEM_USE_HIP)
+  option(MFEM_EXPORT_GPU_CONFIG "Export config.mk for GPU-enabled downstream packages" ON)
+endif()
 mfem_export_mk_files()
@@ -452,13 +452,6 @@ MFEM_USE_MPFR = YES/NO
   quadrature rules. When enabled, this option uses the MPFR_* library options,
   see below.

-MFEM_USE_SIDRE = YES/NO
-   Sidre is a component of LLNL's axom project, https://github.com/LLNL/axom,
-   that provides an HDF5-based file format for visualization or restart
-   capability following the Conduit (https://github.com/LLNL/conduit) mesh
-   blueprint specification. When enabled, this option requires installation of
-   HDF5 (see also MFEM_USE_NETCDF), Conduit and LLNL's axom project.
-
 MFEM_USE_SIMD = YES/NO
   Enables the high performance templated classes to use architecture dependent
   SIMD intrinsics instead of the generic implementation of class AutoSIMD in
@@ -725,7 +718,9 @@ The specific libraries and their options are:
  URL: https://ginkgo-project.github.io
  Options: GINKGO_OPT, GINKGO_LIB, GINKGO_DIR, GINKGO_BUILD_TYPE (Release or
           Debug).
-  Versions: Ginkgo >= 1.9.0.
+  Versions: Ginkgo >= 1.9.0. When building Ginkgo with distributed support, a
+            recent version of the "develop" branch is required (1.11 as defined
+            in include/ginkgo/config.hpp).

 - AmgX (optional), used when MFEM_USE_AMGX = YES.
  URL: https://github.com/NVIDIA/AMGX
@@ -776,14 +771,6 @@ The specific libraries and their options are:
  Options: SLEPC_OPT, SLEPC_LIB.
  Versions: SLEPc >= 3.8.0.

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

-IF (DEFINED TPL_ENABLE_SIDRE)
-  SET(MFEM_USE_SIDRE ${TPL_ENABLE_SIDRE} CACHE BOOL "Enable Axom/Sidre usage" FORCE)
-ENDIF()
-
 IF (DEFINED TPL_ENABLE_FMS)
  SET(MFEM_USE_FMS ${TPL_ENABLE_FMS} CACHE BOOL "Enable FMS usage" FORCE)
 ENDIF()
@@ -46,7 +46,6 @@ set(MFEM_USE_NETCDF @MFEM_USE_NETCDF@)
 set(MFEM_USE_PETSC @MFEM_USE_PETSC@)
 set(MFEM_USE_SLEPC @MFEM_USE_SLEPC@)
 set(MFEM_USE_MPFR @MFEM_USE_MPFR@)
-set(MFEM_USE_SIDRE @MFEM_USE_SIDRE@)
 set(MFEM_USE_FMS @MFEM_USE_FMS@)
 set(MFEM_USE_CONDUIT @MFEM_USE_CONDUIT@)
 set(MFEM_USE_PUMI @MFEM_USE_PUMI@)
@@ -120,9 +120,6 @@
 // Enable secure socket streams based on the GNUTLS library.
 #cmakedefine MFEM_USE_GNUTLS

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

@@ -0,0 +1,24 @@
+# Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+# Defines the following variables:
+#   - ADIAK_FOUND
+#   - ADIAK_LIBRARIES
+#   - ADIAK_INCLUDE_DIRS
+
+include(MfemCmakeUtilities)
+
+mfem_find_package(Adiak ADIAK ADIAK_DIR
+      "include" "adiak.h"
+      "lib" "adiak"
+      "Paths to headers required by Adiak."
+      "Libraries required by Adiak.")
+
@@ -13,6 +13,9 @@
 #   - AXOM_FOUND
 #   - AXOM_LIBRARIES
 #   - AXOM_INCLUDE_DIRS
+#
+# MFEM itself does not depend on Axom, however Tribol does. This module exists
+# to support MFEM's Tribol integration (e.g. the contact miniapp).

 include(MfemCmakeUtilities)
 # Note: components are enabled based on the find_package() parameters.
@@ -0,0 +1,36 @@
+# Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+# Defines the following variables:
+#   - CAMP_FOUND
+#   - CAMP_LIBRARIES
+#   - CAMP_INCLUDE_DIRS
+
+include(MfemCmakeUtilities)
+
+mfem_find_package(CAMP CAMP CAMP_DIR
+      "include" "camp/camp.hpp"
+      "lib" "camp"
+      "Paths to headers required by CAMP."
+      "Libraries required by CAMP.")
+
+# RAJA commonly lists "camp" in INTERFACE_LINK_LIBRARIES. If there is no CMake
+# target named "camp", CMake treats it as a bare library name (-lcamp).
+if (CAMP_FOUND AND NOT TARGET camp)
+   list(GET CAMP_LIBRARIES 0 _camp_lib0)
+   add_library(camp UNKNOWN IMPORTED)
+   set_target_properties(camp PROPERTIES
+      IMPORTED_LOCATION "${_camp_lib0}"
+      INTERFACE_INCLUDE_DIRECTORIES "${CAMP_INCLUDE_DIRS}")
+   set(CAMP_LIBRARIES "camp" CACHE STRING "CAMP imported target." FORCE)
+   unset(_camp_lib0)
+endif()
+
@@ -21,6 +21,21 @@ mfem_find_package(Caliper CALIPER CALIPER_DIR
      "Paths to headers required by Caliper."
      "Libraries required by Caliper.")

+# Some downstream CMake packages (notably RAJA) may list "caliper" in their
+# INTERFACE_LINK_LIBRARIES. If there is no CMake target named "caliper", CMake
+# treats it as a bare library name and will pass -lcaliper to the linker.
+# Create a minimal imported target when we only located the library by path.
+if (CALIPER_FOUND AND NOT TARGET caliper)
+   list(GET CALIPER_LIBRARIES 0 _caliper_lib0)
+   add_library(caliper UNKNOWN IMPORTED)
+   set_target_properties(caliper PROPERTIES
+      IMPORTED_LOCATION "${_caliper_lib0}"
+      INTERFACE_INCLUDE_DIRECTORIES "${CALIPER_INCLUDE_DIRS}")
+   # Prefer linking via the target.
+   set(CALIPER_LIBRARIES "caliper" CACHE STRING "Caliper imported target." FORCE)
+   unset(_caliper_lib0)
+endif()
+
 # Append adiak path/lib if the user provided ADIAK_DIR
 if(ADIAK_DIR AND EXISTS ${ADIAK_DIR})
    find_package(adiak NO_DEFAULT_PATH REQUIRED PATHS ${ADIAK_DIR}/lib/cmake/adiak ${ADIAK_DIR})
@@ -701,7 +701,6 @@ endfunction(mfem_find_library)
 # Extract compile and link options needed by the given target.
 #
 function(mfem_get_target_options Target CompileOptsVar LinkOptsVar)
-
  if (NOT TARGET ${Target})
    return()
  endif()
@@ -799,7 +798,12 @@ function(mfem_get_target_options Target CompileOptsVar LinkOptsVar)
          # message(STATUS "Lib = ${Lib}")
          # Filter-out generator expressions
          if (NOT ("${Lib}" MATCHES "^\\$"))
-            list(APPEND LinkOpts "${Lib}")
+            if(NOT ("${Lib}" STREQUAL "dl"))
+              list(APPEND LinkOpts "${Lib}")
+            else()
+              # for some reason libdl doesn't include the "-l"
+              list(APPEND LinkOpts "-ldl")
+            endif()
          endif()
        else()
          mfem_get_target_options(${Lib} COpts LOpts)
@@ -874,7 +878,7 @@ function(mfem_export_mk_files)
      MFEM_USE_SUITESPARSE MFEM_USE_SUPERLU MFEM_USE_SUPERLU5 MFEM_USE_MUMPS
      MFEM_USE_STRUMPACK MFEM_USE_GINKGO MFEM_USE_AMGX MFEM_USE_MAGMA
      MFEM_USE_GNUTLS MFEM_USE_NETCDF MFEM_USE_PETSC MFEM_USE_SLEPC
-      MFEM_USE_MPFR MFEM_USE_SIDRE MFEM_USE_FMS MFEM_USE_CONDUIT MFEM_USE_PUMI
+      MFEM_USE_MPFR MFEM_USE_FMS MFEM_USE_CONDUIT MFEM_USE_PUMI
      MFEM_USE_HIOP MFEM_USE_GSLIB MFEM_USE_CUDA MFEM_USE_HIP MFEM_USE_RAJA
      MFEM_USE_OCCA MFEM_USE_CEED MFEM_USE_CALIPER MFEM_USE_UMPIRE MFEM_USE_SIMD
      MFEM_USE_ADIOS2 MFEM_USE_MKL_CPARDISO MFEM_USE_MKL_PARDISO
@@ -888,9 +892,18 @@ function(mfem_export_mk_files)
      set(${var} NO)
    endif()
  endforeach()
-  # TODO: Add support for MFEM_USE_CUDA=YES
-  set(MFEM_CXX ${CMAKE_CXX_COMPILER})
-  set(MFEM_HOST_CXX ${MFEM_CXX})
+  if (MFEM_USE_CUDA AND MFEM_EXPORT_GPU_CONFIG)
+    set(MFEM_CXX ${CMAKE_CUDA_COMPILER})
+    if(MFEM_CUDA_COMPILER_IS_NVCC)
+      set(MFEM_HOST_CXX ${CMAKE_CUDA_HOST_COMPILER})
+    else()
+      set(MFEM_HOST_CXX ${CMAKE_CXX_COMPILER})
+    endif()
+  else()
+    # mfem doesn't use enable_language(HIP)
+    set(MFEM_CXX ${CMAKE_CXX_COMPILER})
+    set(MFEM_HOST_CXX ${CMAKE_CXX_COMPILER})
+  endif()
  set(MFEM_CPPFLAGS "")
  get_target_property(cxx_std mfem CXX_STANDARD)
  # For now, we ignore the setting of the CXX_EXTENSIONS property. If this
@@ -900,6 +913,50 @@ function(mfem_export_mk_files)
  string(STRIP
         "${cxx_std_flag} ${CMAKE_CXX_FLAGS_${BUILD_TYPE}} ${CMAKE_CXX_FLAGS}"
         MFEM_CXXFLAGS)
+  if(MFEM_EXPORT_GPU_CONFIG)
+    if (MFEM_USE_CUDA)
+      set(MFEM_CXXFLAGS "${MFEM_CXXFLAGS} ${CMAKE_CUDA_FLAGS}")
+      if (MFEM_CUDA_COMPILER_IS_NVCC)
+        set(MFEM_CXXFLAGS "-x=cu ${MFEM_CXXFLAGS} -ccbin ${CMAKE_CXX_COMPILER} --forward-unknown-to-host-compiler")
+        # The following intentionally hides CUDA deprecation warnings
+        foreach(ENTRY IN LISTS CUDAToolkit_INCLUDE_DIRS)
+          set(MFEM_CXXFLAGS "${MFEM_CXXFLAGS} -isystem ${ENTRY}")
+        endforeach()
+        if (CMAKE_VERSION VERSION_GREATER_EQUAL 3.18.0)
+          # architecture flags not part of CMAKE_CUDA_FLAGS
+          if ("all" STREQUAL "${CMAKE_CUDA_ARCHITECTURES}"
+              OR "native" STREQUAL "${CMAKE_CUDA_ARCHITECTURES}"
+              OR "all-major" STREQUAL "${CMAKE_CUDA_ARCHITECTURES}")
+            set(MFEM_CXXFLAGS "${MFEM_CXXFLAGS} -arch=${CMAKE_CUDA_ARCHITECTURES}")
+          else()
+            foreach (ENTRY IN LISTS CMAKE_CUDA_ARCHITECTURES)
+              set(MFEM_CXXFLAGS
+                "${MFEM_CXXFLAGS} -gencode arch=compute_${ENTRY},code=sm_${ENTRY}")
+            endforeach()
+          endif()
+        endif()
+      else()
+        set(MFEM_CXXFLAGS "${MFEM_CXXFLAGS} -xcuda --cuda-path=${CUDAToolkit_LIBRARY_ROOT}")
+        if (CMAKE_VERSION VERSION_GREATER_EQUAL 3.18.0)
+          # architecture flags not part of CMAKE_CUDA_FLAGS
+          if ("all" STREQUAL "${CMAKE_CUDA_ARCHITECTURES}"
+              OR "native" STREQUAL "${CMAKE_CUDA_ARCHITECTURES}"
+              OR "all-major" STREQUAL "${CMAKE_CUDA_ARCHITECTURES}")
+            # TODO: not supported
+          else()
+            foreach(ENTRY IN LISTS CMAKE_CUDA_ARCHITECTURES)
+              set(MFEM_CXXFLAGS "-cuda-gpu-arch=sm_${ENTRY} ${MFEM_CXXFLAGS}")
+            endforeach()
+          endif()
+        endif()
+      endif()
+    elseif (MFEM_USE_HIP)
+      set(MFEM_CXXFLAGS "${MFEM_CXXFLAGS} -xhip")
+      foreach(ENTRY IN LISTS CMAKE_HIP_ARCHITECTURES)
+        set(MFEM_CXXFLAGS "--offload-arch=${ENTRY} ${MFEM_CXXFLAGS}")
+      endforeach()
+    endif()
+  endif()
  set(MFEM_TPLFLAGS "")
  foreach(dir ${TPL_INCLUDE_DIRS})
    set(MFEM_TPLFLAGS "${MFEM_TPLFLAGS} -I${dir}")
@@ -930,6 +987,9 @@ function(mfem_export_mk_files)
    set(MFEM_SHARED NO)
    set(MFEM_STATIC YES)
  endif()
+  if (MFEM_USE_CUDA)
+    set(MFEM_EXT_LIBS "${MFEM_EXT_LIBS} -lcudart")
+  endif()
  set(MFEM_BUILD_TAG "${CMAKE_SYSTEM}")
  set(MFEM_PREFIX "${CMAKE_INSTALL_PREFIX}")
  # For the next 4 variables, these are the values for the build-tree version of
@@ -938,8 +998,15 @@ function(mfem_export_mk_files)
  set(MFEM_LIB_DIR "${PROJECT_BINARY_DIR}")
  set(MFEM_TEST_MK "${PROJECT_SOURCE_DIR}/config/test.mk")
  set(MFEM_CONFIG_EXTRA "MFEM_BUILD_DIR ?= ${PROJECT_BINARY_DIR}")
-  # TODO: CUDA/HIP support:
-  set(MFEM_XLINKER "${CMAKE_CXX_LINKER_WRAPPER_FLAG}")
+  if (MFEM_USE_CUDA AND MFEM_EXPORT_GPU_CONFIG)
+    if (MFEM_CUDA_COMPILER_IS_NVCC)
+      set(MFEM_XLINKER "-Xlinker=")
+    else()
+      set(MFEM_XLINKER "${CMAKE_CUDA_LINKER_WRAPPER_FLAG}")
+    endif()
+  else()
+    set(MFEM_XLINKER "${CMAKE_CXX_LINKER_WRAPPER_FLAG}")
+  endif()
  set(MFEM_MPIEXEC ${MPIEXEC})
  if (NOT MFEM_MPIEXEC)
    set(MFEM_MPIEXEC "mpirun")
@@ -987,16 +1054,21 @@ function(mfem_export_mk_files)
    # handle interfaces (e.g., SCOREC::apf)
    if ("${lib}" MATCHES "SCOREC::.*" OR "${lib}" MATCHES "Ginkgo::.*" OR "${lib}" MATCHES "ParMoonolith::.*")
    elseif (TARGET "${lib}")
-      mfem_get_target_options(${lib} CompileOpts LinkOpts)
+      mfem_get_target_options(${lib} CompileOpts2 LinkOpts2)
+      # remove generator expressions
+      string(GENEX_STRIP "${CompileOpts2}" CompileOpts)
+      string(GENEX_STRIP "${LinkOpts2}" LinkOpts)
      # Removing duplicates may lead to issues:
      # list(REMOVE_DUPLICATES CompileOpts)
      # list(REMOVE_DUPLICATES LinkOpts)
-      string(REPLACE ";" " " COpts "${CompileOpts}")
-      string(REPLACE ";" " " LOpts "${LinkOpts}")
-      # message(STATUS "${lib}[COpts]: '${COpts}'")
-      # message(STATUS "${lib}[LOpts]: '${LOpts}'")
-      set(MFEM_TPLFLAGS "${MFEM_TPLFLAGS} ${COpts}")
-      set(MFEM_EXT_LIBS "${MFEM_EXT_LIBS} ${LOpts}")
+      # message(WARNING "${lib}[LinkOpts]: ${LinkOpts}")      
+      # message(WARNING "${lib}[CompileOpts]: ${CompileOpts}")
+      foreach(LOpt IN LISTS LinkOpts)
+        set(MFEM_EXT_LIBS "${MFEM_EXT_LIBS} ${LOpt}")
+      endforeach()
+      foreach(COpt IN LISTS CompileOpts)
+        set(MFEM_TPLFLAGS "${MFEM_TPLFLAGS} ${COpt}")
+      endforeach()
      # message(FATAL_ERROR "***** interface lib found ... exiting *****")
      # handle static and shared libs
    elseif ("${suffix}" STREQUAL "${CMAKE_SHARED_LIBRARY_SUFFIX}")
@@ -1004,7 +1076,7 @@ function(mfem_export_mk_files)
      get_filename_component(fullLibName ${lib} NAME_WE)
      string(REGEX REPLACE "^lib" "" libname ${fullLibName})
      set(MFEM_EXT_LIBS
-          "${MFEM_EXT_LIBS} ${shared_link_flag}${dir} -L${dir} -l${libname}")
+        "${MFEM_EXT_LIBS} ${shared_link_flag}${dir} -L${dir} -l${libname}")
    else()
      set(MFEM_EXT_LIBS "${MFEM_EXT_LIBS} ${lib}")
    endif()
@@ -1013,7 +1085,7 @@ function(mfem_export_mk_files)
  # Create the build-tree version of 'config.mk'
  configure_file(
    "${PROJECT_SOURCE_DIR}/config/config.mk.in"
-    "${PROJECT_BINARY_DIR}/config/config.mk")
+    "${PROJECT_BINARY_DIR}/config/config.mk" @ONLY)
  # Copy 'test.mk' from the source-tree to the build-tree
  configure_file(
    "${PROJECT_SOURCE_DIR}/config/test.mk"
@@ -1031,7 +1103,7 @@ function(mfem_export_mk_files)
  # Create the install-tree version of 'config.mk'
  configure_file(
    "${PROJECT_SOURCE_DIR}/config/config.mk.in"
-    "${PROJECT_BINARY_DIR}/config/config-install.mk")
+    "${PROJECT_BINARY_DIR}/config/config-install.mk" @ONLY)

  # Install rules for 'config.mk' and 'test.mk'
  install(FILES ${PROJECT_SOURCE_DIR}/config/test.mk
@@ -120,9 +120,6 @@
 // Enable secure socket streams based on the GNUTLS library.
 // #define MFEM_USE_GNUTLS

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

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

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

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

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

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

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

@@ -18,6 +18,7 @@
 # Some choices below are based on the OS type:
 NOTMAC := $(subst Darwin,,$(shell uname -s))

+ASTYLE_BIN = astyle
 ETAGS_BIN = $(shell command -v etags 2> /dev/null)
 EGREP_BIN = $(shell command -v egrep 2> /dev/null)

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

 # METIS library configuration
 ifeq ($(MFEM_USE_SUPERLU)$(MFEM_USE_STRUMPACK)$(MFEM_USE_MUMPS),NONONO)
+   # MFEM_USE_METIS_5: when the user supplies METIS_DIR, try to auto-detect
+   # METIS 5 installs that follow the common <prefix>/{include,lib,lib64} layout.
+   ifeq ($(MFEM_USE_METIS_5),NO)
+      ifneq ($(wildcard $(METIS_DIR)/include/metis.h),)
+         ifneq ($(wildcard $(METIS_DIR)/lib/libmetis.* $(METIS_DIR)/lib64/libmetis.*),)
+            MFEM_USE_METIS_5 = YES
+         endif
+      endif
+   endif
   ifeq ($(MFEM_USE_METIS_5),NO)
     METIS_DIR = @MFEM_DIR@/../metis-4.0
     METIS_OPT =
@@ -407,7 +416,7 @@ AMGX_LIB = -L$(AMGX_DIR)/lib -lamgx -lcusparse -lcusolver -lcublas -lnvToolsExt
 # MAGMA library configuration
 MAGMA_DIR = @MFEM_DIR@/../magma
 MAGMA_OPT = -I$(MAGMA_DIR)/include
-MAGMA_LIB = -L$(MAGMA_DIR)/lib -l:libmagma.a -lcublas -lcusparse $(LAPACK_LIB)
+MAGMA_LIB = -L$(MAGMA_DIR)/lib -l:libmagma.a $(LAPACK_LIB)

 # GnuTLS library configuration
 GNUTLS_OPT =
@@ -486,17 +495,6 @@ ifneq (,$(wildcard $(CONDUIT_HDF5_HEADER)))
                  -lhdf5 $(ZLIB_LIB)
 endif

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

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

 # Enzyme configuration
 ENZYME_DIR = @MFEM_DIR@/../enzyme
@@ -0,0 +1,86 @@
+MFEM NURBS mesh v1.0
+
+dimension
+1
+
+# Four segments with different NURBS orders, described via patches.
+elements
+4
+1 1 0 1
+2 1 2 3
+3 1 4 5
+4 1 6 7
+
+boundary
+0
+
+edges
+4
+0 0 1
+1 2 3
+2 4 5
+3 6 7
+
+vertices
+8
+
+patches
+
+# Patch 0: linear (order 1, 3 spans)
+knotvectors
+1
+1 4  0 0  .4 .6 1 1
+
+dimension
+2
+
+controlpoints
+0.0  0.0  1.0
+0.6  0.4  1.0
+0.4  0.6  1.0
+1.0  1.0  1.0
+
+# Patch 1: quadratic (order 2, 2 spans)
+knotvectors
+1
+2 4  0 0 0  .5 1 1 1
+
+dimension
+2
+
+controlpoints
+1.0  0.0  1.0
+1.9  0.0  1.21
+2.0  0.9  1.22
+2.0  1.0  1.0
+
+# Patch 2: cubic (order 3, 3 spans)
+knotvectors
+1
+3 6  0 0 0 0 .33 .66 1 1 1 1
+
+dimension
+2
+
+controlpoints
+2.0   0.0  1.0
+2.1   0.2  1.31
+3.5   0.4  1.32
+2.5   0.6  1.33
+2.9   1.0  1.34
+3.0   1.0  1.0
+
+# Patch 3: quartic (order 4, 1 span)
+knotvectors
+1
+4 5  0 0 0 0 0 1 1 1 1 1
+
+dimension
+2
+
+controlpoints
+3.0    0.0  1.0
+3.45   0.5  1.41
+3.50   1.0  1.42
+3.75   0.8  1.43
+4.0    0.0  1.0
@@ -0,0 +1,79 @@
+MFEM NURBS mesh v1.0
+
+#
+# MFEM Geometry Types (see fem/geom.hpp):
+#
+# SEGMENT     = 1
+# SQUARE      = 3
+# CUBE        = 5
+#
+
+dimension
+1
+
+# Three segments with different NURBS orders, described via patches.
+elements
+3
+1 1 0 1
+2 1 2 3
+3 1 4 5
+
+boundary
+6
+1 0 0
+1 0 1
+1 0 2
+1 0 3
+1 0 4
+1 0 5
+
+edges
+3
+0 0 1
+1 2 3
+2 4 5
+
+vertices
+6
+
+patches
+
+# Patch 0: linear (order 1, 2 control points)
+knotvectors
+1
+1 2  0 0  1 1
+
+dimension
+2
+
+controlpoints
+0.0  0.0  1.0
+1.0  1.0  1.0
+
+# Patch 1: quadratic (order 2, 3 control points)
+knotvectors
+1
+2 3  0 0 0  1 1 1
+
+dimension
+2
+
+controlpoints
+1.0  0.0   1.0
+1.02 1.02  1.2
+2.0  1.0   1.0
+
+# Patch 2: cubic (order 3, 4 control points)
+knotvectors
+1
+3 4  0 0 0 0  1 1 1 1
+
+dimension
+2
+
+controlpoints
+2.0   0.0   1.0
+2.03  0.83  1.31
+2.33  1.03  1.32
+3.0   1.0   1.0
+
@@ -0,0 +1,72 @@
+MFEM NURBS mesh v1.0
+
+#
+# MFEM Geometry Types (see fem/geom.hpp):
+#
+# SEGMENT     = 1
+# SQUARE      = 3
+# CUBE        = 5
+#
+
+dimension
+1
+
+elements
+3
+1 1 0 1
+2 1 2 3
+3 1 4 5
+
+boundary
+6
+1 0 0
+1 0 1
+1 0 2
+1 0 3
+1 0 4
+1 0 5
+
+edges
+3
+0 0 1
+1 2 3
+2 4 5
+
+vertices
+6
+
+# Edge 0: linear (order 1, 2 control points)
+# Edge 1: quadratic (order 2, 3 control points)
+# Edge 2: cubic (order 3, 4 control points)
+knotvectors
+3
+1 2  0 0  1 1
+2 3  0 0 0  1 1 1
+3 4  0 0 0 0  1 1 1 1
+
+# One weight per control point, in the same order as the control points; (2 + 3 + 4) = 9 weights total
+weights
+1
+1
+1
+1
+1
+1
+1.2
+1.31
+1.32
+
+FiniteElementSpace
+FiniteElementCollection: NURBS
+VDim: 2
+Ordering: 1
+
+0.0  0.0
+1.0  1.0
+1.0  0.0
+2.0  1.0
+2.0  0.0
+3.0  1.0
+1.02 1.02
+2.03 0.83
+2.33 1.03
@@ -0,0 +1,79 @@
+MFEM NURBS mesh v1.0
+
+#
+# MFEM Geometry Types (see fem/geom.hpp):
+#
+# SEGMENT     = 1
+# SQUARE      = 3
+# CUBE        = 5
+#
+
+dimension
+1
+
+# Three segments with different NURBS orders, described via patches.
+elements
+3
+1 1 0 1
+2 1 2 3
+3 1 4 5
+
+boundary
+6
+1 0 0
+1 0 1
+1 0 2
+1 0 3
+1 0 4
+1 0 5
+
+edges
+3
+0 0 1
+1 2 3
+2 4 5
+
+vertices
+6
+
+patches
+
+# Patch 0: linear (order 1, 2 control points)
+knotvectors
+1
+1 2  0 0  1 1
+
+dimension
+3
+
+controlpoints
+0.0  0.0  0.01 1.0
+1.0  1.0  1.01 1.0
+
+# Patch 1: quadratic (order 2, 3 control points)
+knotvectors
+1
+2 3  0 0 0  1 1 1
+
+dimension
+3
+
+controlpoints
+1.0  0.0   0.02 1.0
+1.02 1.02  0.52 1.2
+2.0  1.0   1.02 1.0
+
+# Patch 2: cubic (order 3, 4 control points)
+knotvectors
+1
+3 4  0 0 0 0  1 1 1 1
+
+dimension
+3
+
+controlpoints
+2.0   0.0   0.03 1.0
+2.03  0.83  0.33 1.31
+2.33  1.03  0.63 1.32
+3.0   1.0   1.03 1.0
+
@@ -0,0 +1,72 @@
+MFEM NURBS mesh v1.0
+
+#
+# MFEM Geometry Types (see fem/geom.hpp):
+#
+# SEGMENT     = 1
+# SQUARE      = 3
+# CUBE        = 5
+#
+
+dimension
+1
+
+elements
+3
+1 1 0 1
+2 1 2 3
+3 1 4 5
+
+boundary
+6
+1 0 0
+1 0 1
+1 0 2
+1 0 3
+1 0 4
+1 0 5
+
+edges
+3
+0 0 1
+1 2 3
+2 4 5
+
+vertices
+6
+
+# Edge 0: linear (order 1, 2 control points)
+# Edge 1: quadratic (order 2, 3 control points)
+# Edge 2: cubic (order 3, 4 control points)
+knotvectors
+3
+1 2  0 0  1 1
+2 3  0 0 0  1 1 1
+3 4  0 0 0 0  1 1 1 1
+
+# One weight per control point, in the same order as the control points; (2 + 3 + 4) = 9 weights total
+weights
+1
+1
+1
+1
+1
+1
+1.2
+1.31
+1.32
+
+FiniteElementSpace
+FiniteElementCollection: NURBS
+VDim: 3
+Ordering: 1
+
+0.0  0.0   0.01
+1.0  1.0   1.01
+1.0  0.0   0.02
+2.0  1.0   1.02
+2.0  0.0   0.03
+3.0  1.0   1.03
+1.02 1.02  0.52
+2.03 0.83  0.33
+2.33 1.03  0.63
@@ -117,6 +117,8 @@ namespace mfem {
 * - <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
+ * - <a class="el" href="ex41_8cpp_source.html">Example 41</a>: DG/CG IMEX time-dependent advection-diffusion
+ * - <a class="el" href="ex41p_8cpp_source.html">Example 41p</a>: parallel DG/CG IMEX time-dependent advection-diffusion
 *
 * <H4>AmgX Examples</H4>
 * - Variants of Examples
@@ -188,6 +190,8 @@ namespace mfem {
 *  <a class="el" href="nurbs__ex1p_8cpp_source.html">1p</a>,
 *  <a class="el" href="nurbs__ex3_8cpp_source.html">3</a>,
 *  <a class="el" href="nurbs__ex5_8cpp_source.html">5</a>,
+ *  <a class="el" href="nurbs__ex10_8cpp_source.html">10</a>,
+ *  <a class="el" href="nurbs__ex10p_8cpp_source.html">10p</a>,
 *  <a class="el" href="nurbs__ex11p_8cpp_source.html">11p</a>, and
 *  <a class="el" href="nurbs__ex24_8cpp_source.html">24</a>,
 *   demonstrating howto perform NURBS-based Isogeometric Analysis.
@@ -196,6 +200,7 @@ namespace mfem {
 * - <a class="el" href="nurbs__curveint_8cpp_source.html">NURBS Interpolation</a>: NURBS interpolation of given geometry
 * - <a class="el" href="nurbs__naca__cmesh_8cpp_source.html">NURBS NACA Mesher</a>: generate NURBS based mesh around a NACA foil
 * - <a class="el" href="nurbs__printfunc_8cpp_source.html">NURBS Printer</a>: print the NURBS-basis
+ * - <a class="el" href="nurbs__mesh_info_8cpp_source.html">NURBS Mesh info</a>: print the info of a NURBS mesh
 *
 * <H3>Miniapps</H3>
 * - <a class="el" href="volta_8cpp_source.html">Volta</a>: simple electrostatics simulation code
@@ -234,7 +239,8 @@ namespace mfem {
 * - <a class="el" href="miniapps_2performance_2ex1_8cpp_source.html">HPC Example 1</a>: high-performance nodal H1 FEM for the Poisson problem
 * - <a class="el" href="miniapps_2performance_2ex1p_8cpp_source.html">HPC Example 1p</a>: high-performance parallel nodal H1 FEM for the Poisson 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="contact-patch-test_8cpp_source.html">Tribol</a>: mortar contact patch test for elasticity
+ * - <a class="el" href="contact_8cpp_source.html">Contact</a>: Frictionless contact examples using <a class="el" href="classmfem_1_1IPSolver.html#details">IP optimization</a> and the <a class="el" href="classmfem_1_1AMGFSolver.html#details">AMGF solver</a>
 * - <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
@@ -5,9 +5,9 @@
 // Sample runs:
 //    mpirun -np 4 ex12p -m ../data/beam-tri.mesh
 //    mpirun -np 4 ex12p -m ../data/beam-quad.mesh
-//    mpirun -np 4 ex12p -m ../data/beam-tet.mesh -s 462 -n 10 -o 2 -elast
+//    mpirun -np 4 ex12p -m ../data/beam-tet.mesh -s 464 -n 10 -o 2 -elast
 //    mpirun -np 4 ex12p -m ../data/beam-hex.mesh -s 3878
-//    mpirun -np 4 ex12p -m ../data/beam-wedge.mesh -s 81
+//    mpirun -np 4 ex12p -m ../data/beam-wedge.mesh -s 82
 //    mpirun -np 4 ex12p -m ../data/beam-tri.mesh -s 3877 -o 2 -sys
 //    mpirun -np 4 ex12p -m ../data/beam-quad.mesh -s 4544 -n 6 -o 3 -elast
 //    mpirun -np 4 ex12p -m ../data/beam-quad-nurbs.mesh
@@ -105,6 +105,7 @@ int main(int argc, char *argv[])
   bool visualization = true;
   bool visit = false;
   int vis_steps = 5;
+   bool solve_implicit_state = false;

   int precision = 8;
   cout.precision(precision);
@@ -126,6 +127,9 @@ int main(int argc, char *argv[])
                  "Alpha coefficient.");
   args.AddOption(&kappa, "-k", "--kappa",
                  "Kappa coefficient offset.");
+   args.AddOption(&solve_implicit_state, "-imp-state", "--implicit-state",
+                  "-imp-slope", "--implicit-slope",
+                  "Implicitly solve for stage state or slope.");
   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                  "--no-visualization",
                  "Enable or disable GLVis visualization.");
@@ -179,6 +183,11 @@ int main(int argc, char *argv[])

   // 7. Initialize the conduction operator and the visualization.
   ConductionOperator oper(fespace, alpha, kappa, u);
+   using ImplicitVariableType = ConductionOperator::ImplicitVariableType;
+   ImplicitVariableType imp_var = solve_implicit_state ?
+                                  ImplicitVariableType::STATE
+                                  : ImplicitVariableType::SLOPE;
+   oper.SetImplicitVariableType(imp_var);

   u_gf.SetFromTrueDofs(u);
   {
@@ -316,11 +325,14 @@ void ConductionOperator::Mult(const Vector &u, Vector &du_dt) const
 }

 void ConductionOperator::ImplicitSolve(const real_t dt,
-                                       const Vector &u, Vector &du_dt)
+                                       const Vector &u, Vector &k)
 {
   // Solve the equation:
-   //    du_dt = M^{-1}*[-K(u + dt*du_dt)]
-   // for du_dt, where K is linearized by using u from the previous timestep
+   //    M*k = -K(u + dt*k) for k = du/dt, if solving for stage-slope
+   // or
+   //    M*k = -dt*K(k) + M*u for k = u_s, if solving for stage-state
+   // where K is linearized by using u from the previous timestep, and
+   // the stage-state and slope relation: du/dt = (u_s - u)/dt.
   if (!T)
   {
      T = Add(1.0, Mmat, dt, Kmat);
@@ -328,9 +340,20 @@ void ConductionOperator::ImplicitSolve(const real_t dt,
      T_solver.SetOperator(*T);
   }
   MFEM_VERIFY(dt == current_dt, ""); // SDIRK methods use the same dt
-   Kmat.Mult(u, z);
-   z.Neg();
-   T_solver.Mult(z, du_dt);
+
+   // Construct current right-hand side for stage state vs. slope solve
+   if (ImplicitVarTypeIsState())
+   {
+      // k, on return, is the stage value u_s
+      Mmat.Mult(u, z);
+   }
+   else
+   {
+      // k, on return, is the stage slope du/dt
+      Kmat.Mult(u, z);
+      z.Neg();
+   }
+   T_solver.Mult(z, k);
 }

 void ConductionOperator::SetParameters(const Vector &u)
@@ -115,6 +115,7 @@ int main(int argc, char *argv[])
   bool visit = false;
   int vis_steps = 5;
   bool adios2 = false;
+   bool solve_implicit_state = false;

   int precision = 8;
   cout.precision(precision);
@@ -138,6 +139,9 @@ int main(int argc, char *argv[])
                  "Alpha coefficient.");
   args.AddOption(&kappa, "-k", "--kappa",
                  "Kappa coefficient offset.");
+   args.AddOption(&solve_implicit_state, "-imp-state", "--implicit-state",
+                  "-imp-slope", "--implicit-slope",
+                  "Implicitly solve for stage state or slope.");
   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                  "--no-visualization",
                  "Enable or disable GLVis visualization.");
@@ -212,6 +216,11 @@ int main(int argc, char *argv[])

   // 9. Initialize the conduction operator and the VisIt visualization.
   ConductionOperator oper(fespace, alpha, kappa, u);
+   using ImplicitVariableType = ConductionOperator::ImplicitVariableType;
+   ImplicitVariableType imp_var = solve_implicit_state ?
+                                  ImplicitVariableType::STATE
+                                  : ImplicitVariableType::SLOPE;
+   oper.SetImplicitVariableType(imp_var);

   u_gf.SetFromTrueDofs(u);
   {
@@ -407,11 +416,14 @@ void ConductionOperator::Mult(const Vector &u, Vector &du_dt) const
 }

 void ConductionOperator::ImplicitSolve(const real_t dt,
-                                       const Vector &u, Vector &du_dt)
+                                       const Vector &u, Vector &k)
 {
   // Solve the equation:
-   //    du_dt = M^{-1}*[-K(u + dt*du_dt)]
-   // for du_dt, where K is linearized by using u from the previous timestep
+   //    M*k = -K(u + dt*k) for k = du/dt, if solving for stage-slope
+   // or
+   //    M*k = -dt*K(k) + M*u for k = u_s, if solving for stage-state
+   // where K is linearized by using u from the previous timestep, and
+   // the stage-state and slope relation: du/dt = (u_s - u)/dt.
   if (!T)
   {
      T = Add(1.0, Mmat, dt, Kmat);
@@ -419,9 +431,20 @@ void ConductionOperator::ImplicitSolve(const real_t dt,
      T_solver.SetOperator(*T);
   }
   MFEM_VERIFY(dt == current_dt, ""); // SDIRK methods use the same dt
-   Kmat.Mult(u, z);
-   z.Neg();
-   T_solver.Mult(z, du_dt);
+
+   // Construct current right-hand side for stage state vs. slope solve
+   if (ImplicitVarTypeIsState())
+   {
+      // k, on return, is the stage value u
+      Mmat.Mult(u, z);
+   }
+   else
+   {
+      // k, on return, is the stage slope du/dt
+      Kmat.Mult(u, z);
+      z.Neg();
+   }
+   T_solver.Mult(z, k);
 }

 void ConductionOperator::SetParameters(const Vector &u)
@@ -119,7 +119,7 @@ int main(int argc, char *argv[])
   }

   LinearForm b(&fespace);
-   b.AddDomainIntegrator(new VectorBoundaryLFIntegrator(f));
+   b.AddBoundaryIntegrator(new VectorBoundaryLFIntegrator(f));

   // 6. Set up the bilinear form a(.,.) on the finite element space
   //    corresponding to the linear elasticity integrator with piece-wise
@@ -140,7 +140,7 @@ int main(int argc, char *argv[])
   }

   ParLinearForm b(&fespace);
-   b.AddDomainIntegrator(new VectorBoundaryLFIntegrator(f));
+   b.AddBoundaryIntegrator(new VectorBoundaryLFIntegrator(f));

   // 6. Set up the bilinear form a(.,.) on the finite element space
   //    corresponding to the linear elasticity integrator with piece-wise
@@ -302,15 +302,21 @@ int main(int argc, char *argv[])
      char vishost[] = "localhost";
      int  visport   = 19916;
      socketstream sol_sock_r(vishost, visport);
-      socketstream sol_sock_i(vishost, visport);
      sol_sock_r << "parallel " << num_procs << " " << myid << "\n";
-      sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
      sol_sock_r.precision(8);
-      sol_sock_i.precision(8);
      sol_sock_r << "solution\n" << *pmesh << u_exact->real()
                 << "window_title 'Exact: Real Part'" << flush;
+      // Make sure all ranks have sent their real solution before initiating
+      // another set of GLVis connections (one from each rank):
+      MPI_Barrier(pmesh->GetComm());
+      socketstream sol_sock_i(vishost, visport);
+      sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock_i.precision(8);
      sol_sock_i << "solution\n" << *pmesh << u_exact->imag()
                 << "window_title 'Exact: Imaginary Part'" << flush;
+      // Make sure all ranks have sent their imaginary solution before initiating
+      // another set of GLVis connections (one from each rank):
+      MPI_Barrier(pmesh->GetComm());
   }

   // 11. Set up the parallel sesquilinear form a(.,.) on the finite element
@@ -534,15 +540,21 @@ int main(int argc, char *argv[])
      char vishost[] = "localhost";
      int  visport   = 19916;
      socketstream sol_sock_r(vishost, visport);
-      socketstream sol_sock_i(vishost, visport);
      sol_sock_r << "parallel " << num_procs << " " << myid << "\n";
-      sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
      sol_sock_r.precision(8);
-      sol_sock_i.precision(8);
      sol_sock_r << "solution\n" << *pmesh << u.real()
                 << "window_title 'Solution: Real Part'" << flush;
+      // Make sure all ranks have sent their real solution before initiating
+      // another set of GLVis connections (one from each rank):
+      MPI_Barrier(pmesh->GetComm());
+      socketstream sol_sock_i(vishost, visport);
+      sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock_i.precision(8);
      sol_sock_i << "solution\n" << *pmesh << u.imag()
                 << "window_title 'Solution: Imaginary Part'" << flush;
+      // Make sure all ranks have sent their imaginary solution before initiating
+      // another set of GLVis connections (one from each rank):
+      MPI_Barrier(pmesh->GetComm());
   }
   if (visualization && exact_sol)
   {
@@ -551,15 +563,21 @@ int main(int argc, char *argv[])
      char vishost[] = "localhost";
      int  visport   = 19916;
      socketstream sol_sock_r(vishost, visport);
-      socketstream sol_sock_i(vishost, visport);
      sol_sock_r << "parallel " << num_procs << " " << myid << "\n";
-      sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
      sol_sock_r.precision(8);
-      sol_sock_i.precision(8);
      sol_sock_r << "solution\n" << *pmesh << u_exact->real()
                 << "window_title 'Error: Real Part'" << flush;
+      // Make sure all ranks have sent their real solution before initiating
+      // another set of GLVis connections (one from each rank):
+      MPI_Barrier(pmesh->GetComm());
+      socketstream sol_sock_i(vishost, visport);
+      sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock_i.precision(8);
      sol_sock_i << "solution\n" << *pmesh << u_exact->imag()
                 << "window_title 'Error: Imaginary Part'" << flush;
+      // Make sure all ranks have sent their imaginary solution before initiating
+      // another set of GLVis connections (one from each rank):
+      MPI_Barrier(pmesh->GetComm());
   }
   if (visualization)
   {
@@ -9,6 +9,7 @@
 //               ex4 -m ../data/beam-hex.mesh -o 2 -pa
 //               ex4 -m ../data/escher.mesh
 //               ex4 -m ../data/fichera.mesh -o 2 -hb
+//               ex4 -m ../data/fichera.mesh -o 2 -hb -ea
 //               ex4 -m ../data/fichera-q2.vtk
 //               ex4 -m ../data/fichera-q3.mesh -o 2 -sc
 //               ex4 -m ../data/square-disc-nurbs.mesh
@@ -18,6 +19,7 @@
 //               ex4 -m ../data/amr-quad.mesh
 //               ex4 -m ../data/amr-hex.mesh
 //               ex4 -m ../data/amr-hex.mesh -o 2 -hb
+//               ex4 -m ../data/amr-hex.mesh -o 2 -hb -ea
 //               ex4 -m ../data/fichera-amr.mesh -o 2 -sc
 //               ex4 -m ../data/ref-prism.mesh -o 1
 //               ex4 -m ../data/octahedron.mesh -o 1
@@ -25,6 +27,8 @@
 //
 // Device sample runs:
 //               ex4 -m ../data/star.mesh -pa -d cuda
+//               ex4 -m ../data/star.mesh -hb -ea -d cuda
+//               ex4 -m ../data/amr-quad.mesh -hb -ea -d cuda
 //               ex4 -m ../data/star.mesh -pa -d raja-cuda
 //               ex4 -m ../data/star.mesh -pa -d raja-omp
 //               ex4 -m ../data/beam-hex.mesh -pa -d cuda
@@ -193,7 +197,7 @@ int main(int argc, char *argv[])
   cout << "Size of linear system: " << A->Height() << endl;

   // 11. Solve the linear system A X = B.
-   if (!pa)
+   if (!pa && (!ea || hybridization))
   {
 #ifndef MFEM_USE_SUITESPARSE
      // Use a simple symmetric Gauss-Seidel preconditioner with PCG.
@@ -433,16 +433,18 @@ int main(int argc, char *argv[])
   u.ProjectCoefficient(*u0);

   // Create data collection for solution output: either VisItDataCollection for
-   // ascii data files, or SidreDataCollection for binary data files.
+   // ascii data files, or ConduitDataCollection for binary data files.
   DataCollection *dc = NULL;
   if (visit)
   {
      if (binary)
      {
-#ifdef MFEM_USE_SIDRE
-         dc = new SidreDataCollection("Example41", &mesh);
+#ifdef MFEM_USE_CONDUIT
+         auto conduit_dc = new ConduitDataCollection("Example41", &mesh);
+         conduit_dc->SetProtocol("hdf5");
+         dc = conduit_dc;
 #else
-         MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
+         MFEM_ABORT("Must build with MFEM_USE_CONDUIT=YES for binary output.");
 #endif
      }
      else
@@ -518,10 +518,12 @@ int main(int argc, char *argv[])
   {
      if (binary)
      {
-#ifdef MFEM_USE_SIDRE
-         dc = new SidreDataCollection("Example41-Parallel", pmesh);
+#ifdef MFEM_USE_CONDUIT
+         auto conduit_dc = new ConduitDataCollection("Example41-Parallel", pmesh);
+         conduit_dc->SetProtocol("hdf5");
+         dc = conduit_dc;
 #else
-         MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
+         MFEM_ABORT("Must build with MFEM_USE_CONDUIT=YES for binary output.");
 #endif
      }
      else
@@ -9,6 +9,7 @@
 //               mpirun -np 4 ex4p -m ../data/beam-hex.mesh -o 2 -pa
 //               mpirun -np 4 ex4p -m ../data/escher.mesh -o 2 -sc
 //               mpirun -np 4 ex4p -m ../data/fichera.mesh -o 2 -hb
+//               mpirun -np 4 ex4p -m ../data/fichera.mesh -o 2 -hb -ea
 //               mpirun -np 4 ex4p -m ../data/fichera-q2.vtk
 //               mpirun -np 4 ex4p -m ../data/fichera-q3.mesh -o 2 -sc
 //               mpirun -np 4 ex4p -m ../data/square-disc-nurbs.mesh -o 3
@@ -17,14 +18,18 @@
 //               mpirun -np 4 ex4p -m ../data/periodic-cube.mesh -no-bc
 //               mpirun -np 4 ex4p -m ../data/amr-quad.mesh
 //               mpirun -np 3 ex4p -m ../data/amr-quad.mesh -o 2 -hb
+//               mpirun -np 3 ex4p -m ../data/amr-quad.mesh -o 2 -hb -ea
 //               mpirun -np 4 ex4p -m ../data/amr-hex.mesh -o 2 -sc
 //               mpirun -np 4 ex4p -m ../data/amr-hex.mesh -o 2 -hb
+//               mpirun -np 4 ex4p -m ../data/amr-hex.mesh -o 2 -hb -ea
 //               mpirun -np 4 ex4p -m ../data/ref-prism.mesh -o 1
 //               mpirun -np 4 ex4p -m ../data/octahedron.mesh -o 1
 //               mpirun -np 4 ex4p -m ../data/star-surf.mesh -o 3 -hb
 //
 // Device sample runs:
 //               mpirun -np 4 ex4p -m ../data/star.mesh -pa -d cuda
+//               mpirun -np 4 ex4p -m ../data/star.mesh -ea -hb -d cuda
+//               mpirun -np 4 ex4p -m ../data/amr-hex.mesh -ea -hb -d cuda
 //               mpirun -np 4 ex4p -m ../data/star.mesh -pa -d raja-cuda
 //               mpirun -np 4 ex4p -m ../data/star.mesh -pa -d raja-omp
 //               mpirun -np 4 ex4p -m ../data/beam-hex.mesh -pa -d cuda
@@ -230,7 +235,7 @@ int main(int argc, char *argv[])
   pcg->SetMaxIter(2000);
   pcg->SetPrintLevel(1);
   if (hybridization) { prec = new HypreBoomerAMG(*A.As<HypreParMatrix>()); }
-   else if (pa) { prec = new OperatorJacobiSmoother(*a, ess_tdof_list); }
+   else if (pa || ea) { prec = new OperatorJacobiSmoother(*a, ess_tdof_list); }
   else
   {
      ParFiniteElementSpace *prec_fespace =
@@ -160,6 +160,7 @@ int main(int argc, char *argv[])
   bool paraview = false;
   bool binary = false;
   int vis_steps = 5;
+   bool solve_implicit_state = false;

   int precision = 8;
   cout.precision(precision);
@@ -187,6 +188,9 @@ int main(int argc, char *argv[])
                  "Final time; start time is 0.");
   args.AddOption(&dt, "-dt", "--time-step",
                  "Time step.");
+   args.AddOption(&solve_implicit_state, "-imp-state", "--implicit-state",
+                  "-imp-slope", "--implicit-slope",
+                  "Implicitly solve for stage state or slope.");
   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                  "--no-visualization",
                  "Enable or disable GLVis visualization.");
@@ -301,16 +305,18 @@ int main(int argc, char *argv[])
   }

   // Create data collection for solution output: either VisItDataCollection for
-   // ascii data files, or SidreDataCollection for binary data files.
+   // ascii data files, or ConduitDataCollection for binary data files.
   DataCollection *dc = NULL;
   if (visit)
   {
      if (binary)
      {
-#ifdef MFEM_USE_SIDRE
-         dc = new SidreDataCollection("Example9", &mesh);
+#ifdef MFEM_USE_CONDUIT
+         auto conduit_dc = new ConduitDataCollection("Example9", &mesh);
+         conduit_dc->SetProtocol("hdf5");
+         dc = conduit_dc;
 #else
-         MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
+         MFEM_ABORT("Must build with MFEM_USE_CONDUIT=YES for binary output.");
 #endif
      }
      else
@@ -366,6 +372,11 @@ int main(int argc, char *argv[])
   //    right-hand side, and perform time-integration (looping over the time
   //    iterations, ti, with a time-step dt).
   FE_Evolution adv(m, k, b);
+   using ImplicitVariableType = FE_Evolution::ImplicitVariableType;
+   ImplicitVariableType imp_var = solve_implicit_state ?
+                                  ImplicitVariableType::STATE
+                                  : ImplicitVariableType::SLOPE;
+   adv.SetImplicitVariableType(imp_var);

   real_t t = 0.0;
   adv.SetTime(t);
@@ -459,7 +470,17 @@ void FE_Evolution::ImplicitSolve(const real_t dt, const Vector &x, Vector &k)
 {
   MFEM_VERIFY(dg_solver != NULL,
               "Implicit time integration is not supported with partial assembly");
-   K.Mult(x, z);
+   // Construct current right-hand side for stage state vs. slope solve
+   if (ImplicitVarTypeIsState())
+   {
+      // k, on return, is the stage value u
+      M.Mult(x, z);
+   }
+   else
+   {
+      // k, on return, is the stage slope du/dt
+      K.Mult(x, z);
+   }
   z += b;
   dg_solver->SetTimeStep(dt);
   dg_solver->Mult(z, k);
@@ -257,6 +257,7 @@ int main(int argc, char *argv[])
   bool adios2 = false;
   bool binary = false;
   int vis_steps = 5;
+   bool solve_implicit_state = false;
 #if MFEM_HYPRE_VERSION >= 21800
   PrecType prec_type = PrecType::AIR;
 #else
@@ -290,6 +291,9 @@ int main(int argc, char *argv[])
                  "Final time; start time is 0.");
   args.AddOption(&dt, "-dt", "--time-step",
                  "Time step.");
+   args.AddOption(&solve_implicit_state, "-imp-state", "--implicit-state",
+                  "-imp-slope", "--implicit-slope",
+                  "Implicitly solve for stage state or slope.");
   args.AddOption((int *)&prec_type, "-pt", "--prec-type", "Preconditioner for "
                  "implicit solves. 0 for ILU, 1 for pAIR-AMG.");
   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
@@ -437,16 +441,18 @@ int main(int argc, char *argv[])
   }

   // Create data collection for solution output: either VisItDataCollection for
-   // ascii data files, or SidreDataCollection for binary data files.
+   // ascii data files, or ConduitDataCollection for binary data files.
   DataCollection *dc = NULL;
   if (visit)
   {
      if (binary)
      {
-#ifdef MFEM_USE_SIDRE
-         dc = new SidreDataCollection("Example9-Parallel", pmesh);
+#ifdef MFEM_USE_CONDUIT
+         auto conduit_dc = new ConduitDataCollection("Example9-Parallel", pmesh);
+         conduit_dc->SetProtocol("hdf5");
+         dc = conduit_dc;
 #else
-         MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
+         MFEM_ABORT("Must build with MFEM_USE_CONDUIT=YES for binary output.");
 #endif
      }
      else
@@ -536,6 +542,11 @@ int main(int argc, char *argv[])
   //     right-hand side, and perform time-integration (looping over the time
   //     iterations, ti, with a time-step dt).
   FE_Evolution adv(*m, *k, *B, prec_type);
+   using ImplicitVariableType = FE_Evolution::ImplicitVariableType;
+   ImplicitVariableType imp_var = solve_implicit_state ?
+                                  ImplicitVariableType::STATE
+                                  : ImplicitVariableType::SLOPE;
+   adv.SetImplicitVariableType(imp_var);

   real_t t = 0.0;
   adv.SetTime(t);
@@ -676,7 +687,17 @@ FE_Evolution::FE_Evolution(ParBilinearForm &M_, ParBilinearForm &K_,
 //    (M - dt*K) d = K*u + b
 void FE_Evolution::ImplicitSolve(const real_t dt, const Vector &x, Vector &k)
 {
-   K->Mult(x, z);
+   // Construct current right-hand side for stage state vs. slope solve
+   if (ImplicitVarTypeIsState())
+   {
+      // k, on return, is the stage value u
+      M->Mult(x, z);
+   }
+   else
+   {
+      // k, on return, is the stage slope du/dt
+      K->Mult(x, z);
+   }
   z += b;
   dg_solver->SetTimeStep(dt);
   dg_solver->Mult(z, k);
@@ -14,6 +14,12 @@ list(APPEND GINKGO_EXAMPLES_SRCS
  ex1.cpp
 )

+if (MFEM_USE_MPI AND GINKGO_BUILD_MPI)
+    list(APPEND GINKGO_EXAMPLES_SRCS
+    ex1p.cpp
+)
+endif()
+
 # Include the source directory where mfem.hpp and mfem-performance.hpp are.
 include_directories(BEFORE ${PROJECT_BINARY_DIR})

@@ -207,7 +207,7 @@ int main(int argc, char *argv[])
            Ginkgo::IcPreconditioner ginkgo_precond(exec, "paric", 30);
            Ginkgo::CGSolver ginkgo_solver(exec, ginkgo_precond);
            ginkgo_solver.SetPrintLevel(print_lvl);
-            ginkgo_solver.SetRelTol(1e-12);
+            ginkgo_solver.SetRelTol(sqrt(1e-12));
            ginkgo_solver.SetAbsTol(0.0);
            ginkgo_solver.SetMaxIter(400);
            ginkgo_solver.SetOperator(*(A.Ptr()));
@@ -225,7 +225,7 @@ int main(int argc, char *argv[])
            Ginkgo::MFEMPreconditioner gko_M(exec, M);
            Ginkgo::CGSolver ginkgo_solver(exec, gko_M);
            ginkgo_solver.SetPrintLevel(print_lvl);
-            ginkgo_solver.SetRelTol(1e-12);
+            ginkgo_solver.SetRelTol(sqrt(1e-12));
            ginkgo_solver.SetAbsTol(0.0);
            ginkgo_solver.SetMaxIter(400);
            ginkgo_solver.SetOperator(*(A.Ptr()));
@@ -283,7 +283,7 @@ int main(int argc, char *argv[])
               Ginkgo::MFEMPreconditioner gko_M(exec, M);
               Ginkgo::CGSolver ginkgo_solver(exec, gko_M);
               ginkgo_solver.SetPrintLevel(print_lvl);
-               ginkgo_solver.SetRelTol(1e-12);
+               ginkgo_solver.SetRelTol(sqrt(1e-12));
               ginkgo_solver.SetAbsTol(0.0);
               ginkgo_solver.SetMaxIter(400);
               ginkgo_solver.SetOperator(*(A.Ptr()));
@@ -0,0 +1,436 @@
+//                       MFEM Example 1 - Parallel Version
+//                             GINKGO Modification
+//
+// Compile with: make ex1p
+//
+// Sample runs:  mpirun -np 4 ex1p -m ../data/square-disc.mesh
+//               mpirun -np 4 ex1p -m ../data/star.mesh
+//               mpirun -np 4 ex1p -m ../data/star-mixed.mesh
+//               mpirun -np 4 ex1p -m ../data/escher.mesh
+//               mpirun -np 4 ex1p -m ../data/fichera.mesh
+//               mpirun -np 4 ex1p -m ../data/fichera-mixed.mesh
+//               mpirun -np 4 ex1p -m ../data/toroid-wedge.mesh
+//               mpirun -np 4 ex1p -m ../data/octahedron.mesh -o 1
+//               mpirun -np 4 ex1p -m ../data/periodic-annulus-sector.msh
+//               mpirun -np 4 ex1p -m ../data/periodic-torus-sector.msh
+//               mpirun -np 4 ex1p -m ../data/square-disc-p2.vtk -o 2
+//               mpirun -np 4 ex1p -m ../data/square-disc-p3.mesh -o 3
+//               mpirun -np 4 ex1p -m ../data/square-disc-nurbs.mesh -o -1
+//               mpirun -np 4 ex1p -m ../data/star-mixed-p2.mesh -o 2
+//               mpirun -np 4 ex1p -m ../data/disc-nurbs.mesh -o -1
+//               mpirun -np 4 ex1p -m ../data/pipe-nurbs.mesh -o -1
+//               mpirun -np 4 ex1p -m ../data/ball-nurbs.mesh -o 2
+//               mpirun -np 4 ex1p -m ../data/fichera-mixed-p2.mesh -o 2
+//               mpirun -np 4 ex1p -m ../data/star-surf.mesh
+//               mpirun -np 4 ex1p -m ../data/square-disc-surf.mesh
+//               mpirun -np 4 ex1p -m ../data/inline-segment.mesh
+//               mpirun -np 4 ex1p -m ../data/amr-quad.mesh
+//               mpirun -np 4 ex1p -m ../data/amr-hex.mesh
+//               mpirun -np 4 ex1p -m ../data/mobius-strip.mesh
+//               mpirun -np 4 ex1p -m ../data/mobius-strip.mesh -o -1 -sc
+//
+// Device sample runs:
+//               mpirun -np 4 ex1p -pa -d cuda
+//               mpirun -np 4 ex1p -fa -d cuda
+//               mpirun -np 4 ex1p -pa -d occa-cuda
+//               mpirun -np 4 ex1p -pa -d raja-omp
+//               mpirun -np 4 ex1p -pa -d ceed-cpu
+//               mpirun -np 4 ex1p -pa -d ceed-cpu -o 4 -a
+//               mpirun -np 4 ex1p -pa -d ceed-cpu -m ../data/square-mixed.mesh
+//               mpirun -np 4 ex1p -pa -d ceed-cpu -m ../data/fichera-mixed.mesh
+//             * mpirun -np 4 ex1p -pa -d ceed-cuda
+//             * mpirun -np 4 ex1p -pa -d ceed-hip
+//               mpirun -np 4 ex1p -pa -d ceed-cuda:/gpu/cuda/shared
+//               mpirun -np 4 ex1p -pa -d ceed-cuda:/gpu/cuda/shared -m ../data/square-mixed.mesh
+//               mpirun -np 4 ex1p -pa -d ceed-cuda:/gpu/cuda/shared -m ../data/fichera-mixed.mesh
+//               mpirun -np 4 ex1p -m ../data/beam-tet.mesh -pa -d ceed-cpu
+//
+// Description:  This example code demonstrates the use of MFEM to define a
+//               simple finite element discretization of the Laplace problem
+//               -Delta u = 1 with homogeneous Dirichlet boundary conditions.
+//               Specifically, we discretize using a FE space of the specified
+//               order, or if order < 1 using an isoparametric/isogeometric
+//               space (i.e. quadratic for quadratic curvilinear mesh, NURBS for
+//               NURBS mesh, etc.)
+//
+//               The example highlights the use of mesh refinement, finite
+//               element grid functions, as well as linear and bilinear forms
+//               corresponding to the left-hand side and right-hand side of the
+//               discrete linear system. We also cover the explicit elimination
+//               of essential boundary conditions, static condensation, and the
+//               optional connection to the GLVis tool for visualization.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+#ifndef MFEM_USE_GINKGO
+#error This example requires that MFEM is built with MFEM_USE_GINKGO=YES
+#endif
+
+using namespace std;
+using namespace mfem;
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI and HYPRE.
+   Mpi::Init();
+   int num_procs = Mpi::WorldSize();
+   int myid = Mpi::WorldRank();
+   Hypre::Init();
+
+   // 2. Parse command-line options.
+   const char *mesh_file = "../../data/star.mesh";
+   int order = 1;
+   bool static_cond = false;
+   bool pa = false;
+   bool fa = false;
+   const char *device_config = "cpu";
+   bool visualization = true;
+   int solver_config = 0;
+   int print_lvl = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) or -1 for"
+                  " isoparametric space.");
+   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
+                  "--no-static-condensation", "Enable static condensation.");
+   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
+                  "--no-partial-assembly", "Enable Partial Assembly.");
+   args.AddOption(&fa, "-fa", "--full-assembly", "-no-fa",
+                  "--no-full-assembly", "Enable Full Assembly.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&solver_config, "-s", "--solver-config",
+                  "Solver and preconditioner combination: \n\t"
+                  "   0 - Ginkgo solver and Ginkgo preconditioner, \n\t"
+                  "   1 - Ginkgo solver and MFEM preconditioner, \n\t"
+                  "   2 - MFEM solver and Ginkgo preconditioner, \n\t"
+                  "   3 - MFEM solver and MFEM preconditioner.");
+   args.AddOption(&print_lvl, "-pl", "--print-level",
+                  "Print level for iterative solver (1 prints every iteration).");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   // 3. Enable hardware devices such as GPUs, and programming models such as
+   //    CUDA, OCCA, RAJA and OpenMP based on command line options.
+   Device device(device_config);
+   device.SetGPUAwareMPI(true);
+   if (myid == 0) { device.Print(); }
+
+   // 4. Read the (serial) mesh from the given mesh file on all processors.  We
+   //    can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
+   //    and volume meshes with the same code.
+   Mesh mesh(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+
+   // 5. Refine the serial mesh on all processors to increase the resolution. In
+   //    this example we do 'ref_levels' of uniform refinement. We choose
+   //    'ref_levels' to be the largest number that gives a final mesh with no
+   //    more than 10,000 elements.
+   {
+      int ref_levels =
+         (int)floor(log(10000./mesh.GetNE())/log(2.)/dim);
+      for (int l = 0; l < ref_levels; l++)
+      {
+         mesh.UniformRefinement();
+      }
+   }
+
+   // 6. Define a parallel mesh by a partitioning of the serial mesh. Refine
+   //    this mesh further in parallel to increase the resolution. Once the
+   //    parallel mesh is defined, the serial mesh can be deleted.
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
+   {
+      int par_ref_levels = 2;
+      for (int l = 0; l < par_ref_levels; l++)
+      {
+         pmesh.UniformRefinement();
+      }
+   }
+
+   // 7. Define a parallel finite element space on the parallel mesh. Here we
+   //    use continuous Lagrange finite elements of the specified order. If
+   //    order < 1, we instead use an isoparametric/isogeometric space.
+   FiniteElementCollection *fec;
+   bool delete_fec;
+   if (order > 0)
+   {
+      fec = new H1_FECollection(order, dim);
+      delete_fec = true;
+   }
+   else if (pmesh.GetNodes())
+   {
+      fec = pmesh.GetNodes()->OwnFEC();
+      delete_fec = false;
+      if (myid == 0)
+      {
+         cout << "Using isoparametric FEs: " << fec->Name() << endl;
+      }
+   }
+   else
+   {
+      fec = new H1_FECollection(order = 1, dim);
+      delete_fec = true;
+   }
+   ParFiniteElementSpace fespace(&pmesh, fec);
+   HYPRE_BigInt size = fespace.GlobalTrueVSize();
+   if (myid == 0)
+   {
+      cout << "Number of finite element unknowns: " << size << endl;
+   }
+
+   // 8. Determine the list of true (i.e. parallel conforming) essential
+   //    boundary dofs. In this example, the boundary conditions are defined
+   //    by marking all the boundary attributes from the mesh as essential
+   //    (Dirichlet) and converting them to a list of true dofs.
+   Array<int> ess_tdof_list;
+   if (pmesh.bdr_attributes.Size())
+   {
+      Array<int> ess_bdr(pmesh.bdr_attributes.Max());
+      ess_bdr = 1;
+      fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   // 9. Set up the parallel linear form b(.) which corresponds to the
+   //    right-hand side of the FEM linear system, which in this case is
+   //    (1,phi_i) where phi_i are the basis functions in fespace.
+   ParLinearForm b(&fespace);
+   ConstantCoefficient one(1.0);
+   b.AddDomainIntegrator(new DomainLFIntegrator(one));
+   b.Assemble();
+
+   // 10. Define the solution vector x as a parallel finite element grid
+   //     function corresponding to fespace. Initialize x with initial guess of
+   //     zero, which satisfies the boundary conditions.
+   ParGridFunction x(&fespace);
+   x = 0.0;
+
+   // 11. Set up the parallel bilinear form a(.,.) on the finite element space
+   //     corresponding to the Laplacian operator -Delta, by adding the
+   //     Diffusion domain integrator.
+   ParBilinearForm a(&fespace);
+   if (pa) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
+   if (fa)
+   {
+      a.SetAssemblyLevel(AssemblyLevel::FULL);
+      // Sort the matrix column indices when running on GPU or with OpenMP (i.e.
+      // when Device::IsEnabled() returns true). This makes the results
+      // bit-for-bit deterministic at the cost of somewhat longer run time.
+      a.EnableSparseMatrixSorting(Device::IsEnabled());
+   }
+   a.AddDomainIntegrator(new DiffusionIntegrator(one));
+
+   // 12. Assemble the parallel bilinear form and the corresponding linear
+   //     system, applying any necessary transformations such as: parallel
+   //     assembly, eliminating boundary conditions, applying conforming
+   //     constraints for non-conforming AMR, static condensation, etc.
+   if (static_cond) { a.EnableStaticCondensation(); }
+   a.Assemble();
+
+   OperatorPtr A;
+   Vector B, X;
+   a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
+
+   // 13. Solve the linear system A X = B.
+   if (!pa)
+   {
+      switch (solver_config)
+      {
+         // Solve the linear system with CG + Schwarz (with IC) from Ginkgo
+         case 0:
+         {
+            if (myid == 0) { cout << "Using Ginkgo solver + preconditioner...\n"; }
+            Ginkgo::GinkgoExecutor exec(device);
+            Ginkgo::IcPreconditioner local_solver(exec, "exact");
+            Ginkgo::SchwarzPreconditioner gko_M(exec, MPI_COMM_WORLD, local_solver);
+            Ginkgo::CGSolver ginkgo_solver(exec, MPI_COMM_WORLD, gko_M);
+            ginkgo_solver.SetPrintLevel(print_lvl);
+            ginkgo_solver.SetRelTol(sqrt(1e-12));
+            ginkgo_solver.SetAbsTol(0.0);
+            ginkgo_solver.SetMaxIter(400);
+            ginkgo_solver.SetOperator(*(A.Ptr()));
+            ginkgo_solver.Mult(B, X);
+            break;
+         }
+
+         // Solve the linear system with CG from Ginkgo + MFEM preconditioner
+         case 1:
+         {
+            if (myid == 0) { cout << "Using Ginkgo solver + MFEM preconditioner...\n"; }
+            Ginkgo::GinkgoExecutor exec(device);
+            //Create MFEM preconditioner and wrap it for Ginkgo's use.
+            HypreBoomerAMG M((HypreParMatrix&)(*A));
+            Ginkgo::MFEMPreconditioner gko_M(exec, M, MPI_COMM_WORLD);
+            Ginkgo::CGSolver ginkgo_solver(exec, MPI_COMM_WORLD, gko_M);
+            ginkgo_solver.SetPrintLevel(print_lvl);
+            ginkgo_solver.SetRelTol(sqrt(1e-12));
+            ginkgo_solver.SetAbsTol(0.0);
+            ginkgo_solver.SetMaxIter(400);
+            ginkgo_solver.SetOperator(*(A.Ptr()));
+            ginkgo_solver.Mult(B, X);
+            break;
+         }
+
+         // Ginkgo Schwarz preconditioner (local ParIC) +  MFEM CG solver
+         case 2:
+         {
+            if (myid == 0) { cout << "Using MFEM solver + Ginkgo preconditioner...\n"; }
+            Ginkgo::GinkgoExecutor exec(device);
+            Ginkgo::IcPreconditioner local_M(exec, "exact");
+            Ginkgo::SchwarzPreconditioner M(exec, MPI_COMM_WORLD, local_M);
+            M.SetOperator(*(A.Ptr()));  // Generate the preconditioner for the matrix A.
+            CGSolver cg(MPI_COMM_WORLD);
+            cg.SetRelTol(sqrt(1e-12));
+            cg.SetMaxIter(400);
+            cg.SetPrintLevel(1);
+            cg.SetPreconditioner(M);
+            cg.SetOperator(*A);
+            cg.Mult(B, X);
+            break;
+         }
+
+         // MFEM solver + MFEM preconditioner
+         case 3:
+         {
+            if (myid == 0) { cout << "Using MFEM solver + MFEM preconditioner...\n"; }
+            HypreBoomerAMG M((HypreParMatrix&)(*A));
+            CGSolver cg(MPI_COMM_WORLD);
+            cg.SetRelTol(sqrt(1e-12));
+            cg.SetMaxIter(400);
+            cg.SetPrintLevel(1);
+            cg.SetPreconditioner(M);
+            cg.SetOperator(*A);
+            cg.Mult(B, X);
+            break;
+         }
+      } // End switch on solver_config
+   }
+   // Partial assembly mode.  Cannot use Ginkgo preconditioners, but can use Ginkgo
+   // solvers.
+   else
+   {
+      if (UsesTensorBasis(fespace))
+      {
+         // Use Jacobi preconditioning in partial assembly mode.
+         OperatorJacobiSmoother M(a, ess_tdof_list);
+         switch (solver_config)
+         {
+            case 0:
+            {
+               if (myid == 0) { cout << "Using Ginkgo solver + preconditioner...\n"; }
+               MFEM_ABORT("Cannot use Ginkgo preconditioner in partial assembly mode.\n"
+                          "            Try -s 1 to test Ginkgo solver with an MFEM preconditioner.");
+               break;
+            }
+
+            // Use Ginkgo solver with MFEM preconditioner
+            case 1:
+            {
+               if (myid == 0) { cout << "Using Ginkgo solver + MFEM preconditioner...\n"; }
+               Ginkgo::GinkgoExecutor exec(device);
+               // Wrap MFEM preconditioner for Ginkgo's use.
+               Ginkgo::MFEMPreconditioner gko_M(exec, M, MPI_COMM_WORLD);
+               Ginkgo::CGSolver ginkgo_solver(exec, MPI_COMM_WORLD, gko_M);
+               ginkgo_solver.SetPrintLevel(print_lvl);
+               ginkgo_solver.SetRelTol(sqrt(1e-12));
+               ginkgo_solver.SetAbsTol(0.0);
+               ginkgo_solver.SetMaxIter(400);
+               ginkgo_solver.SetOperator(*(A.Ptr()));
+               ginkgo_solver.Mult(B, X);
+               break;
+            }
+
+            // No Ginkgo preconditioners work with matrix-free; error
+            case 2:
+            {
+               if (myid == 0) { cout << "Using Ginkgo solver + preconditioner...\n"; }
+               MFEM_ABORT("Cannot use Ginkgo preconditioner in partial assembly mode.\n"
+                          "            Try -s 1 to test Ginkgo solver with an MFEM preconditioner.");
+               break;
+            }
+
+            // Use MFEM solver and preconditioner
+            case 3:
+            {
+               if (myid == 0) { cout << "Using MFEM solver + MFEM preconditioner...\n"; }
+               CGSolver cg(MPI_COMM_WORLD);
+               cg.SetRelTol(sqrt(1e-12));
+               cg.SetMaxIter(400);
+               cg.SetPrintLevel(1);
+               cg.SetPreconditioner(M);
+               cg.SetOperator(*A);
+               cg.Mult(B, X);
+               break;
+            }
+         } // End switch on solver_config
+      }
+      else // CG with no preconditioning
+      {
+         if (myid == 0) { cout << "Using MFEM solver + no preconditioner...\n"; }
+         CGSolver cg(MPI_COMM_WORLD);
+         cg.SetRelTol(sqrt(1e-12));
+         cg.SetMaxIter(400);
+         cg.SetPrintLevel(1);
+         cg.SetOperator(*A);
+         cg.Mult(B, X);
+      }
+   }
+
+   // 14. Recover the parallel grid function corresponding to X. This is the
+   //     local finite element solution on each processor.
+   a.RecoverFEMSolution(X, b, x);
+
+   // 15. Save the refined mesh and the solution in parallel. This output can
+   //     be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
+   {
+      ostringstream mesh_name, sol_name;
+      mesh_name << "mesh." << setfill('0') << setw(6) << myid;
+      sol_name << "sol." << setfill('0') << setw(6) << myid;
+
+      ofstream mesh_ofs(mesh_name.str().c_str());
+      mesh_ofs.precision(8);
+      pmesh.Print(mesh_ofs);
+
+      ofstream sol_ofs(sol_name.str().c_str());
+      sol_ofs.precision(8);
+      x.Save(sol_ofs);
+   }
+
+   // 16. Send the solution by socket to a GLVis server.
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << pmesh << x << flush;
+   }
+
+   // 17. Free the used memory.
+   if (delete_fec)
+   {
+      delete fec;
+   }
+
+   return 0;
+}
@@ -20,9 +20,8 @@ CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
 MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)

-# Currently there are only serial Ginkgo examples
 SEQ_EXAMPLES = ex1
-PAR_EXAMPLES =
+PAR_EXAMPLES = ex1p
 ifeq ($(MFEM_USE_MPI),NO)
   EXAMPLES = $(SEQ_EXAMPLES)
 else
@@ -354,16 +354,18 @@ int main(int argc, char *argv[])
   }

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

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

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

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

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

-if (MFEM_USE_SIDRE)
-  list(APPEND SRCS sidredatacollection.cpp)
-  list(APPEND HDRS sidredatacollection.hpp)
-endif()
-
 if (MFEM_USE_CONDUIT)
  list(APPEND SRCS conduitdatacollection.cpp)
  list(APPEND HDRS conduitdatacollection.hpp)
@@ -729,7 +729,8 @@ void BilinearForm::Assemble(int skip_zeros)
         tr = mesh -> GetBdrFaceTransformations (i);
         if (tr != NULL)
         {
-            fes -> GetElementVDofs (tr -> Elem1No, vdofs);
+            mfem::DofTransformation doftrans;
+            fes -> GetElementVDofs (tr -> Elem1No, vdofs, doftrans);
            fe1 = fes -> GetFE (tr -> Elem1No);
            // The fe2 object is really a dummy and not used on the boundaries,
            // but we can't dereference a NULL pointer, and we don't want to
@@ -743,6 +744,7 @@ void BilinearForm::Assemble(int skip_zeros)

               boundary_face_integs[k] -> AssembleFaceMatrix (*fe1, *fe2, *tr,
                                                              elemmat);
+               doftrans.TransformDual(elemmat);
               mat -> AddSubMatrix (vdofs, vdofs, elemmat, skip_zeros);
            }
         }
@@ -825,14 +827,46 @@ void BilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list, Vector &x,
                                    Vector &b, OperatorHandle &A, Vector &X,
                                    Vector &B, int copy_interior)
 {
+   const SparseMatrix *P = fes->GetConformingProlongation();
+   const SparseMatrix *R = fes->GetConformingRestriction();
   if (ext)
   {
      if (hybridization)
      {
         FormSystemMatrix(ess_tdof_list, A);
-         ConstrainedOperator A_constrained(this, ess_tdof_list);
-         A_constrained.EliminateRHS(x, b);
-         hybridization->ReduceRHS(b, B);
+
+         std::unique_ptr<ConstrainedOperator> A_constrained([&]()
+         {
+            Operator *op;
+            Operator::FormSystemOperator(ess_tdof_list, op);
+            return dynamic_cast<ConstrainedOperator*>(op);
+         }());
+         MFEM_ASSERT(A_constrained != nullptr, "");
+
+         Vector conf_b, conf_x;
+         if (P)
+         {
+            // Nonconforming
+            conf_b.SetSize(P->Width());
+            conf_x.SetSize(P->Width());
+            P->MultTranspose(b, conf_b);
+            R->Mult(x, conf_x);
+         }
+         else
+         {
+            // Conforming
+            conf_b.MakeRef(b, 0, b.Size());
+            conf_x.MakeRef(x, 0, x.Size());
+         }
+
+         A_constrained->EliminateRHS(conf_x, conf_b);
+
+         if (P)
+         {
+            R->MultTranspose(conf_b, b); // store eliminated rhs in b
+         }
+
+         hybridization->ReduceRHS(conf_b, B);
         X.SetSize(B.Size());
         X = 0.0;
      }
@@ -842,7 +876,6 @@ void BilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list, Vector &x,
      }
      return;
   }
-   const SparseMatrix *P = fes->GetConformingProlongation();
   FormSystemMatrix(ess_tdof_list, A);

   // Transform the system and perform the elimination in B, based on the
@@ -878,7 +911,6 @@ void BilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list, Vector &x,
      if (hybridization)
      {
         // Reduction to the Lagrange multipliers system
-         const SparseMatrix *R = fes->GetConformingRestriction();
         Vector conf_b(P->Width()), conf_x(P->Width());
         P->MultTranspose(b, conf_b);
         R->Mult(x, conf_x);
@@ -891,7 +923,6 @@ void BilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list, Vector &x,
      else
      {
         // Variational restriction with P
-         const SparseMatrix *R = fes->GetConformingRestriction();
         B.SetSize(P->Width());
         P->MultTranspose(b, B);
         X.SetSize(R->Height());
@@ -1694,6 +1725,7 @@ void MixedBilinearForm::Assemble(int skip_zeros)
         }
      }

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


 /** Integrator for $(-Q u, \nabla v)$ for Nedelec ($u$) and $H^1$ ($v$) elements.
-    This is equivalent to a weak divergence of the $H(curl$ basis functions. */
+    This is equivalent to a weak divergence of the $H(curl)$ basis functions. */
 class VectorFEWeakDivergenceIntegrator: public BilinearFormIntegrator
 {
 protected:
@@ -39,8 +39,8 @@ void PLBound::Setup(const int nb_i, const int ncp_i,
   b_type = b_type_i;
   cp_type = cp_type_i;
   tol = tol_i;
-   lbound.SetSize(nb, ncp);
-   ubound.SetSize(nb, ncp);
+   lbound.SetSize(ncp, nb);
+   ubound.SetSize(ncp, nb);
   nodes.SetSize(nb);
   weights.SetSize(nb);
   control_points.SetSize(ncp);
@@ -125,21 +125,25 @@ void PLBound::Setup(const int nb_i, const int ncp_i,
      {
         if (j == 0)
         {
-            lbound(i, j) = bv(i);
-            ubound(i, j) = bv(i);
+            lbound(j,i) = bv(i);
+            ubound(j,i) = bv(i);
         }
         else if (j == ncp-1)
         {
-            lbound(i, j) = bv(i);
-            ubound(i, j) = bv(i);
+            lbound(j,i) = bv(i);
+            ubound(j,i) = bv(i);
         }
         else
         {
            vals(0)  = bv(i);
            vals(1) =  bmv(i) +  dm*bdmv(i);
            vals(2) =  bpv(i) +  dp*bdpv(i);
-            lbound(i, j) = vals.Min()-tol; // tolerance for good measure
-            ubound(i, j) = vals.Max()+tol; // tolerance for good measure
+            lbound(j,i) = vals.Min()-tol; // tolerance for good measure
+            ubound(j,i) = vals.Max()+tol; // tolerance for good measure
+            if (b_type == 2)
+            {
+               lbound(j,i) = std::max(lbound(j,i),0_r);
+            }
         }
      }
   }
@@ -273,8 +277,7 @@ void PLBound::Get1DBounds(const Vector &coeff, Vector &intmin,
   intmax.SetSize(ncp);
   intmin = 0.0;
   intmax = 0.0;
-   Vector coeffm(nb);
-   coeffm = 0.0;
+   Vector coeffm;

   real_t a0 = 0.0;
   real_t a1 = 0.0;
@@ -302,6 +305,8 @@ void PLBound::Get1DBounds(const Vector &coeff, Vector &intmin,
   // compute L2 projection for linear bases: a0 + a1*x
   if (proj)
   {
+      coeffm.SetSize(nb);
+      coeffm = 0.0;
      for (int i = 0; i < nb; i++)
      {
         x = 2.0*nodes_int(i)-1;
@@ -342,8 +347,8 @@ void PLBound::Get1DBounds(const Vector &coeff, Vector &intmin,
      real_t c = coeffm(i);
      for (int j = 0; j < ncp; j++)
      {
-         intmin(j) += min(lbound(i,j)*c, ubound(i,j)*c);
-         intmax(j) += max(lbound(i,j)*c, ubound(i,j)*c);
+         intmin(j) += min(lbound(j,i)*c, ubound(j,i)*c);
+         intmax(j) += max(lbound(j,i)*c, ubound(j,i)*c);
      }
   }
 }
@@ -474,10 +479,10 @@ void PLBound::Get2DBounds(const Vector &coeff, Vector &intmin,
         real_t w1 = intmaxT(id2++);
         for (int k = 0; k < ncp; k++) // kth row
         {
-            vals(0) = w0*lbound(j,k);
-            vals(1) = w0*ubound(j,k);
-            vals(2) = w1*lbound(j,k);
-            vals(3) = w1*ubound(j,k);
+            vals(0) = w0*lbound(k,j);
+            vals(1) = w0*ubound(k,j);
+            vals(2) = w1*lbound(k,j);
+            vals(3) = w1*ubound(k,j);
            intmin(k*ncp+i) += vals.Min();
            intmax(k*ncp+i) += vals.Max();
         }
@@ -553,17 +558,17 @@ void PLBound::Get3DBounds(const Vector &coeff, Vector &intmin,
            for (int i = 0; i < nb; i++)
            {
               x = 2.0*nodes(i)-1; // x-coordinate
-               minBounds(i) -= a0V(j) + a1V(j)*x;
-               maxBounds(i) -= a0V(j) + a1V(j)*x;
+               minNodalVals(i) -= a0V(j) + a1V(j)*x;
+               maxNodalVals(i) -= a0V(j) + a1V(j)*x;
            }
            // Compute Bernstein coefficients
            LUFactors lu(basisMatLU.GetData(), lu_ip.GetData());
-            lu.Solve(nb, 1, minBounds.GetData());
-            lu.Solve(nb, 1, maxBounds.GetData());
+            lu.Solve(nb, 1, minNodalVals.GetData());
+            lu.Solve(nb, 1, maxNodalVals.GetData());
            for (int i = 0; i < nb; i++)
            {
-               intminT(i*ncp2+j) = minBounds(i);
-               intmaxT(i*ncp2+j) = maxBounds(i);
+               intminT(i*ncp2+j) = minNodalVals(i);
+               intmaxT(i*ncp2+j) = maxNodalVals(i);
            }
         }
      }
@@ -617,10 +622,10 @@ void PLBound::Get3DBounds(const Vector &coeff, Vector &intmin,
         real_t w1 = intmaxT(id2++);
         for (int k = 0; k < ncp; k++) // kth slice
         {
-            vals(0) = w0*lbound(j,k);
-            vals(1) = w0*ubound(j,k);
-            vals(2) = w1*lbound(j,k);
-            vals(3) = w1*ubound(j,k);
+            vals(0) = w0*lbound(k,j);
+            vals(1) = w0*ubound(k,j);
+            vals(2) = w1*lbound(k,j);
+            vals(3) = w1*ubound(k,j);
            intmin(k*ncp2+i) += vals.Min();
            intmax(k*ncp2+i) += vals.Max();
         }
@@ -653,7 +658,8 @@ void PLBound::SetupBernsteinBasisMat(DenseMatrix &basisMat,
                                     Vector &nodesBern) const
 {
   const int nbern = nodesBern.Size();
-   L2_SegmentElement el(nbern-1, 2); // we use L2 to leverage lexicographic order
+   L2_SegmentElement el(nbern-1, 2);
+   // we use L2 to leverage lexicographic order
   Array<int> ordering = el.GetLexicographicOrdering();
   basisMat.SetSize(nbern, nbern);
   Vector shape(nbern);
@@ -666,6 +672,39 @@ void PLBound::SetupBernsteinBasisMat(DenseMatrix &basisMat,
   }
 }

+DenseMatrix PLBound::GetBoundingMatrix(int dim, bool is_lower) const
+{
+   if (dim > 1)
+   {
+      const int ncpd = static_cast<int>(std::pow(ncp, dim));
+      const int nbd = static_cast<int>(std::pow(nb, dim));
+      DenseMatrix boundND(ncpd, nbd);
+      Vector phimin, phimax, col;
+      Vector coeffs(nbd);
+      coeffs = 0.0;
+      for (int j = 0; j < nbd; j++)
+      {
+         coeffs(j) = 1.0;
+         boundND.GetColumnReference(j, col);
+         GetNDBounds(dim, coeffs, phimin, phimax);
+         col = is_lower ? phimin : phimax;
+         coeffs(j) = 0.0;
+      }
+      return boundND;
+   }
+   return is_lower ? lbound : ubound;
+}
+
+DenseMatrix PLBound::GetLowerBoundMatrix(int dim) const
+{
+   return GetBoundingMatrix(dim, true);
+}
+
+DenseMatrix PLBound::GetUpperBoundMatrix(int dim) const
+{
+   return GetBoundingMatrix(dim, false);
+}
+
 constexpr int PLBound::min_ncp_gl_x[2][11];
 constexpr int PLBound::min_ncp_gll_x[2][11];
 constexpr int PLBound::min_ncp_pos_x[2][11];
@@ -716,4 +755,4 @@ void PLBound::Print(std::ostream &outp) const
   ubound.Print(outp);
 }

-}
+}
@@ -19,14 +19,18 @@ namespace mfem
 {

 /** @name Piecewise linear bounds of bases
-    \brief Piecewise linear bounds of bases can be used to compute bounds on the grid function in each element. The bounds for the bases are constructed based on the following parameters:
+    \brief Piecewise linear bounds of bases can be used to compute bounds on
+    the grid function in each element. The bounds for the bases are constructed
+    based on the following parameters:

    (i) @b nb: number of bases/nodes in 1D (i.e. polynomial order+1),

-    (ii) @b b_type: bases type, 0 - Lagrange interpolants on Gauss-Legendre nodes, 1 - Lagrange interpolants on Gauss-Lobatto-Legendre nodes, and
+    (ii) @b b_type: bases type, 0 - Lagrange interpolants on Gauss-Legendre
+    nodes, 1 - Lagrange interpolants on Gauss-Lobatto-Legendre nodes, and
    2 - Positive/Bernstein bases on uniformly distributed nodes,

-    (iii) @b ncp: number of control points used to construct the piecewise linear bounds
+    (iii) @b ncp: number of control points used to construct the piecewise
+    linear bounds

    (iv) @b cp_type: control point distribution. 0 - GL + end-points,
                    1 - Chebyshev.
@@ -35,7 +39,9 @@ namespace mfem

    If the user does not specify @b ncp and @b cp_type, the minimum value of
    @b ncp is used that would bound the bases for the @b cp_type. We default
-    to @b cp_type = 0 as it requires fewer number of points to bound the bases. Typically, @b ncp = 2 @b nb is sufficient to get fairly compact bounds, and increasing @b ncp results in tighter bounds.
+    to @b cp_type = 0 as it requires fewer number of points to bound the bases.
+    Typically, @b ncp = 2 @b nb is sufficient to get fairly compact bounds, and
+    increasing @b ncp results in tighter bounds.

    Finally, only tensor-product elements are currently supported.

@@ -54,7 +60,7 @@ private:
   bool proj = true; // Use linear projection to compute bounds.
   real_t tol = 0.0; // offset bounds to avoid round-off errors
   Vector nodes, weights, control_points;
-   DenseMatrix lbound, ubound; // nb x ncp matrices with bounds of all bases
+   DenseMatrix lbound, ubound; // ncp x nb matrices with bounds of all bases
   // Some auxillary storage for computing the bounds with Bernstein
   DenseMatrix basisMatNodes; // Bernstein bases at equispaced nodes
   DenseMatrix basisMatInt;   // Bernstein bases at GLL nodes
@@ -80,6 +86,9 @@ private:
      {3,5,8,9,11,12,13,13,14,15,16}
   };

+   /// Helper function to extract lower or upper bounding matrix
+   DenseMatrix GetBoundingMatrix(int dim, bool is_lower) const;
+
 public:
   // Constructor
   PLBound(const int nb_i, const int ncp_i, const int b_type_i,
@@ -92,40 +101,82 @@ public:
   PLBound(const FiniteElementSpace *fes,
           const int ncp_i = -1, const int cp_type_i = 0);

-   // Get minimum number of control points needed to bound the given bases
+   /// Get minimum number of control points needed to bound the given bases
   int GetMinimumPointsForGivenBases(int nb_i, int b_type_i,
                                     int cp_type_i) const;

-   // Print information about the bounds
+   /// Print information about the bounds
   void Print(std::ostream &outp = mfem::out) const;

-   // Enable (default) or disable linear projection before bounding.
-   // This projection increases the computational cost but results in tighter
-   // bounds.
+   /** @brief Enable (default) or disable linear projection before bounding.
+    *
+    *  @details This projection increases the computational cost but results in
+    *  tighter bounds.
+    */
   void SetProjectionFlagForBounding(bool proj_) { proj = proj_; }

-   /// Compute piecewise linear bounds for the lexicographically-ordered
-   /// coefficients in @a coeff in 1D/2D/3D.
+   /** @brief Compute piecewise linear bounds for the lexicographically-ordered
+    *  nodal coefficients in @a coeff in 1D/2D/3D.
+    *
+    *  @param[in] rdim    The spatial dimension of the element (1, 2, or 3).
+    *  @param[in] coeff   The vector of lexicographically-ordered coefficients.
+    *                     Should be of size nb^rdim, where nb is the number of
+    *                     bases/nodes in 1D. These coefficients must correspond
+    *                     to the bases type and number of bases, used in the
+    *                     constructor of PLBound.
+    *
+    *  @param[out] intmin The vector of minimum bound for all control points.
+    *  @param[out] intmax The vector of maximum bound for all control points.
+    *                     Both intmin and intmax are of size ncp^rdim, where
+    *                     ncp is the number of control points in 1D, and are
+    *                     ordered lexicographically.
+    */
   void GetNDBounds(const int rdim, const Vector &coeff,
                    Vector &intmin, Vector &intmax) const;

   /// Get number of control points used to compute the bounds.
   int GetNControlPoints() const { return ncp; }
+
+   /// Get 1D control point locations (lexicographic order) in [0,1].
+   const Vector &GetControlPoints() const { return control_points; }
+
+   /** @brief Get lower and upper bounding matrix (ncp^dim x nb^dim)
+    *
+    *  @details The matrices can be used to compute the bounds at control points
+    *           by a simple matrix-vector product with the
+    *           lexicographically-ordered nodal coefficients.
+    *           The resulting output is also lexicographically-ordered.
+    *
+    *  @note These matrices do not account for the linear projection step that
+    *        is optionally done in GetNDBounds before bounding the function.
+    */
+   ///@{
+   DenseMatrix GetLowerBoundMatrix(int dim = 1) const;
+   DenseMatrix GetUpperBoundMatrix(int dim = 1) const;
+   ///@}
+
 private:
-   /// Compute piecewise linear bounds for the lexicographically-ordered
-   /// coefficients in @a coeff in 1D.
+   /** @brief Compute piecewise linear bounds for the lexicographically-ordered
+    *  nodal coefficients in @a coeff in 1D.
+    *  See GetNDBounds for details of the input and output parameters.
+    */
   void Get1DBounds(const Vector &coeff, Vector &intmin, Vector &intmax) const;

-   /// Compute piecewise linear bounds for the lexicographically-ordered
-   /// coefficients in @a coeff in 2D.
+   /** @brief Compute piecewise linear bounds for the lexicographically-ordered
+    *  nodal coefficients in @a coeff in 2D.
+    *  See GetNDBounds for details of the input and output parameters.
+    */
   void Get2DBounds(const Vector &coeff, Vector &intmin, Vector &intmax) const;

-   /// Compute piecewise linear bounds for the lexicographically-ordered
-   /// coefficients in @a coeff in 3D.
+   /** @brief Compute piecewise linear bounds for the lexicographically-ordered
+    *  nodal coefficients in @a coeff in 3D.
+    *  See GetNDBounds for details of the input and output parameters.
+    */
   void Get3DBounds(const Vector &coeff, Vector &intmin, Vector &intmax) const;

-   /// Setup matrix used to compute values at given 1D locations in [0,1]
-   /// for Bernstein bases.
+   /** @brief Setup matrix used to compute values at given 1D locations in [0,1]
+    *  for Bernstein bases.
+    */
   void SetupBernsteinBasisMat(DenseMatrix &basisMat, Vector &nodesBern) const;

   void Setup(const int nb_i, const int ncp_i, const int b_type_i,
@@ -1302,6 +1302,73 @@ real_t TraceCoefficient::Eval(ElementTransformation &T,
   return ma.Trace();
 }

+VectorComponentCoefficient::VectorComponentCoefficient(VectorCoefficient &A,
+                                                       int c)
+   : a(&A), va(A.GetVDim())
+{
+   SetComponent(c);
+}
+
+void VectorComponentCoefficient::SetComponent(int c)
+{
+   MFEM_ASSERT(c < a->GetVDim() && c >= 0,
+               "VectorComponentCoefficient:  "
+               "Index not in range.");
+
+   component = c;
+}
+
+void VectorComponentCoefficient::SetTime(real_t t)
+{
+   if (a) { a->SetTime(t); }
+   this->Coefficient::SetTime(t);
+}
+
+real_t VectorComponentCoefficient::Eval(ElementTransformation &T,
+                                        const IntegrationPoint &ip)
+{
+   a->Eval(va, T, ip);
+   return va[component];
+}
+
+MatrixComponentCoefficient::MatrixComponentCoefficient(MatrixCoefficient &A,
+                                                       int ri, int ci)
+   : a(&A), ma(A.GetHeight(), A.GetWidth())
+{
+   SetRowIndex(ri);
+   SetColumnIndex(ci);
+}
+
+void MatrixComponentCoefficient::SetRowIndex(int ri)
+{
+   MFEM_ASSERT(ri < a->GetHeight() && ri >= 0,
+               "MatrixComponentCoefficient:  "
+               "Row index not in range.");
+
+   row_idx = ri;
+}
+
+void MatrixComponentCoefficient::SetColumnIndex(int ci)
+{
+   MFEM_ASSERT(ci < a->GetWidth() && ci >= 0,
+               "MatrixComponentCoefficient:  "
+               "Column index not in range.");
+   col_idx = ci;
+}
+
+void MatrixComponentCoefficient::SetTime(real_t t)
+{
+   if (a) { a->SetTime(t); }
+   this->Coefficient::SetTime(t);
+}
+
+real_t MatrixComponentCoefficient::Eval(ElementTransformation &T,
+                                        const IntegrationPoint &ip)
+{
+   a->Eval(ma, T, ip);
+   return ma(row_idx,col_idx);
+}
+
 VectorSumCoefficient::VectorSumCoefficient(int dim)
   : VectorCoefficient(dim),
     ACoef(NULL), BCoef(NULL),
@@ -114,11 +114,10 @@ public:
   /// Construct the constant coefficient using a vector of constants.
   /** @a c should be a vector defined by attributes, so for region with
       attribute @a i @a c[i-1] is the coefficient in that region */
-   PWConstCoefficient(Vector &c)
-   { constants.SetSize(c.Size()); constants=c; }
+   PWConstCoefficient(const Vector &c) { UpdateConstants(c); }

   /// Update the constants with vector @a c.
-   void UpdateConstants(Vector &c) { constants.SetSize(c.Size()); constants=c; }
+   void UpdateConstants(const Vector &c) { constants = c; }

   /// Return a reference to the i-th constant
   real_t &operator()(int i) { return constants(i-1); }
@@ -1332,8 +1331,8 @@ public:
   /// Get the coefficient located at (i,j) in the matrix.
   Coefficient* GetCoeff (int i, int j) { return Coeff[i*width+j]; }

-   /** @brief Set the coefficient located at (i,j) in the matrix.  By default by
-       default this will take ownership of the Coefficient passed in, but this
+   /** @brief Set the coefficient located at (i,j) in the matrix.  By default
+       this will take ownership of the Coefficient passed in, but this
       can be overridden with the @a own parameter. */
   void Set(int i, int j, Coefficient * c, bool own=true);

@@ -1873,6 +1872,85 @@ public:
               const IntegrationPoint &ip) override;
 };

+/// Scalar coefficient defined as component of a vector coefficient
+class VectorComponentCoefficient : public Coefficient
+{
+private:
+   VectorCoefficient *a = nullptr;
+
+   mutable Vector va;
+   int component;
+
+public:
+   /// Construct with a vector coefficient.
+   VectorComponentCoefficient(VectorCoefficient &A)
+      : a(&A), va(A.GetVDim()), component(0) {};
+
+   VectorComponentCoefficient(VectorCoefficient &A, int c);
+
+   /// Set the time for internally stored coefficients
+   void SetTime(real_t t) override;
+
+   /// Reset the vector coefficient
+   void SetACoef(VectorCoefficient &A) { a = &A; }
+
+   /// Return the vector coefficient
+   VectorCoefficient * GetACoef() const { return a; }
+
+   /// Set the component
+   void SetComponent(int c);
+
+   /// Return the component
+   int GetComponent() const { return component; }
+
+   /// Evaluate the trace coefficient at @a ip.
+   real_t Eval(ElementTransformation &T,
+               const IntegrationPoint &ip) override;
+};
+
+/// Scalar coefficient defined as component of a matrix coefficient
+class MatrixComponentCoefficient : public Coefficient
+{
+private:
+   MatrixCoefficient *a = nullptr;
+
+   mutable DenseMatrix ma;
+   int row_idx,col_idx;
+
+public:
+   MatrixComponentCoefficient(MatrixCoefficient &A)
+      : a(&A), ma(A.GetHeight(), A.GetWidth()), row_idx(0), col_idx(0) {};
+
+   /// Construct with the matrix coefficient.
+   MatrixComponentCoefficient(MatrixCoefficient &A, int ri, int ci);
+
+   /// Set the time for internally stored coefficients
+   void SetTime(real_t t) override;
+
+   /// Reset the matrix coefficient
+   void SetACoef(MatrixCoefficient &A) { a = &A; }
+
+   /// Return the matrix coefficient
+   MatrixCoefficient * GetACoef() const { return a; }
+
+   /// Reset the index
+   void SetRowIndex(int ri);
+
+   /// Return the index
+   int GetRowIndex() const { return row_idx; }
+
+   /// Reset the index
+   void SetColumnIndex(int ci);
+
+   /// Return the index
+   int GetColumnIndex() const { return col_idx; }
+
+
+   /// Evaluate the trace coefficient at @a ip.
+   real_t Eval(ElementTransformation &T,
+               const IntegrationPoint &ip) override;
+};
+
 /// Vector coefficient defined as the linear combination of two vectors
 class VectorSumCoefficient : public VectorCoefficient
 {
@@ -82,6 +82,25 @@ public:
   /// underlying #fes
   int VectorDim() const;

+   /// Copy assignment. Only the data of the base class Vector is copied.
+   /** It is assumed that this object and @a rhs use FiniteElementSpace%s that
+       have the same size.
+
+       @note Defining this method overwrites the implicitly defined copy
+       assignment operator. */
+   ComplexGridFunction &operator=(const ComplexGridFunction &rhs)
+   { return operator=((const Vector &)rhs); }
+
+   /// Copy the data from @a v.
+   /** The size of @a v must be equal to double of the size of the associated
+       FiniteElementSpace #fes. */
+   ComplexGridFunction &operator=(const Vector &v)
+   {
+      MFEM_ASSERT(fes && v.Size() == 2*fes->GetVSize(), "");
+      Vector::operator=(v);
+      return *this;
+   }
+
   /// Assign constant values to the ComplexGridFunction data.
   ComplexGridFunction &operator=(const std::complex<real_t> & value)
   { *gfr = value.real(); *gfi = value.imag(); return *this; }
@@ -1453,8 +1453,6 @@ ConduitDataCollection::LoadMeshAndFields(int domain_id,
 std::string
 ConduitDataCollection::ElementTypeToShapeName(Element::Type element_type)
 {
-   // Adapted from SidreDataCollection
-
   // Note -- the mapping from Element::Type to string is based on
   //   enum Element::Type { POINT, SEGMENT, TRIANGLE, QUADRILATERAL,
   //                        TETRAHEDRON, HEXAHEDRON };
@@ -34,10 +34,10 @@ namespace mfem
    - HDF5 library, https://support.hdfgroup.org/HDF5

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

    This class also provides public static methods that convert between MFEM
    Meshes and GridFunctions and Conduit Mesh Blueprint descriptions.
@@ -90,8 +90,8 @@ void map_quadrature_data_to_fields_impl(
   }
   else
   {
-      MFEM_ABORT("quadrature data mapping to field is not implemented for"
-                 " this field descriptor");
+      MFEM_ABORT_KERNEL("quadrature data mapping to field is not implemented"
+                        " for this field descriptor");
   }
 }

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

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

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

@@ -10,7 +10,7 @@
 // CONTRIBUTING.md for details.
 #pragma once

-// This is serac's tuple implementation
+// This is smith's tuple implementation

 #include <ostream>
 #include "../../config/config.hpp"
@@ -12,7 +12,6 @@
 #include "dgmassinv.hpp"
 #include "bilinearform.hpp"
 #include "dgmassinv_kernels.hpp"
-#include "../general/forall.hpp"

 namespace mfem
 {
@@ -119,151 +118,6 @@ void DGMassInverse::Update()

 DGMassInverse::~DGMassInverse() = default;

-template<int DIM, int D1D, int Q1D>
-void DGMassInverse::DGMassCGIteration(const Vector &b_, Vector &u_) const
-{
-   using namespace internal; // host/device kernel functions
-
-   const int NE = fes.GetNE();
-   const int d1d = m->dofs1D;
-   const int q1d = m->quad1D;
-
-   const int ND = static_cast<int>(pow(d1d, DIM));
-
-   const auto B = m->maps->B.Read();
-   const auto Bt = m->maps->Bt.Read();
-   const auto pa_data = m->pa_data.Read();
-   const auto dinv = diag_inv.Read();
-   auto r = r_.Write();
-   auto d = d_.Write();
-   auto z = z_.Write();
-   auto u = u_.ReadWrite();
-
-   const real_t RELTOL = rel_tol;
-   const real_t ABSTOL = abs_tol;
-   const int MAXIT = max_iter;
-   const bool IT_MODE = iterative_mode;
-   const bool CHANGE_BASIS = (d2q != nullptr);
-
-   // b is the right-hand side (if no change of basis, this just points to the
-   // incoming RHS vector, if we have to change basis, this points to the
-   // internal b2 vector where we put the transformed RHS)
-   const real_t *b;
-   // the following are non-null if we have to change basis
-   real_t *b2 = nullptr; // non-const access to b2
-   const real_t *b_orig = nullptr; // RHS vector in "original" basis
-   const real_t *d2q_B = nullptr; // matrix to transform initial guess
-   const real_t *q2d_B = nullptr; // matrix to transform solution
-   const real_t *q2d_Bt = nullptr; // matrix to transform RHS
-   if (CHANGE_BASIS)
-   {
-      d2q_B = d2q->B.Read();
-      q2d_B = B_.Read();
-      q2d_Bt = Bt_.Read();
-
-      b2 = b2_.Write();
-      b_orig = b_.Read();
-      b = b2;
-   }
-   else
-   {
-      b = b_.Read();
-   }
-
-   static constexpr int NB = Q1D ? Q1D : 1; // block size
-
-   mfem::forall_2D(NE, NB, NB, [=] MFEM_HOST_DEVICE (int e)
-   {
-      // Perform change of basis if needed
-      if (CHANGE_BASIS)
-      {
-         // Transform RHS
-         DGMassBasis<DIM,D1D>(e, NE, q2d_Bt, b_orig, b2, d1d);
-         if (IT_MODE)
-         {
-            // Transform initial guess
-            DGMassBasis<DIM,D1D>(e, NE, d2q_B, u, u, d1d);
-         }
-      }
-
-      const int tid = MFEM_THREAD_ID(x) + NB*MFEM_THREAD_ID(y);
-
-      // Compute first residual
-      if (IT_MODE)
-      {
-         DGMassApply<DIM,D1D,Q1D>(e, NE, B, Bt, pa_data, u, r, d1d, q1d);
-         DGMassAxpy(e, NE, ND, 1.0, b, -1.0, r, r); // r = b - r
-      }
-      else
-      {
-         // if not in iterative mode, use zero initial guess
-         const int BX = MFEM_THREAD_SIZE(x);
-         const int BY = MFEM_THREAD_SIZE(y);
-         const int bxy = BX*BY;
-         const auto B = ConstDeviceMatrix(b, ND, NE);
-         auto U = DeviceMatrix(u, ND, NE);
-         auto R = DeviceMatrix(r, ND, NE);
-         for (int i = tid; i < ND; i += bxy)
-         {
-            U(i, e) = 0.0;
-            R(i, e) = B(i, e);
-         }
-         MFEM_SYNC_THREAD;
-      }
-
-      DGMassPreconditioner(e, NE, ND, dinv, r, z);
-      DGMassAxpy(e, NE, ND, 1.0, z, 0.0, z, d); // d = z
-
-      real_t nom = DGMassDot<NB>(e, NE, ND, d, r);
-      if (nom < 0.0) { return; /* Not positive definite */ }
-      real_t r0 = fmax(nom*RELTOL*RELTOL, ABSTOL*ABSTOL);
-      if (nom <= r0) { return; /* Converged */ }
-
-      DGMassApply<DIM,D1D,Q1D>(e, NE, B, Bt, pa_data, d, z, d1d, q1d);
-      real_t den = DGMassDot<NB>(e, NE, ND, z, d);
-      if (den <= 0.0)
-      {
-         DGMassDot<NB>(e, NE, ND, d, d);
-         // d2 > 0 => not positive definite
-         if (den == 0.0) { return; }
-      }
-
-      // start iteration
-      int i = 1;
-      while (true)
-      {
-         const real_t alpha = nom/den;
-         DGMassAxpy(e, NE, ND, 1.0, u, alpha, d, u); // u = u + alpha*d
-         DGMassAxpy(e, NE, ND, 1.0, r, -alpha, z, r); // r = r - alpha*A*d
-
-         DGMassPreconditioner(e, NE, ND, dinv, r, z);
-
-         real_t betanom = DGMassDot<NB>(e, NE, ND, r, z);
-         if (betanom < 0.0) { return; /* Not positive definite */ }
-         if (betanom <= r0) { break; /* Converged */ }
-
-         if (++i > MAXIT) { break; }
-
-         const real_t beta = betanom/nom;
-         DGMassAxpy(e, NE, ND, 1.0, z, beta, d, d); // d = z + beta*d
-         DGMassApply<DIM,D1D,Q1D>(e, NE, B, Bt, pa_data, d, z, d1d, q1d); // z = A d
-         den = DGMassDot<NB>(e, NE, ND, d, z);
-         if (den <= 0.0)
-         {
-            DGMassDot<NB>(e, NE, ND, d, d);
-            // d2 > 0 => not positive definite
-            if (den == 0.0) { break; }
-         }
-         nom = betanom;
-      }
-
-      if (CHANGE_BASIS)
-      {
-         DGMassBasis<DIM,D1D>(e, NE, q2d_B, u, u, d1d);
-      }
-   });
-}
-
 void DGMassInverse::Mult(const Vector &Mu, Vector &u) const
 {
   // Dispatch to templated version based on dim, d1d, and q1d.
@@ -306,23 +160,4 @@ DGMassInvKernels::DGMassInvKernels()
   k::Specialization<3,6,7>::Add();
 }

-/// @cond Suppress_Doxygen_warnings
-
-template <int DIM, int D1D, int Q1D>
-DGMassInverse::CGKernelType DGMassInverse::CGKernels::Kernel()
-{
-   return &DGMassInverse::DGMassCGIteration<DIM,D1D,Q1D>;
-}
-
-DGMassInverse::CGKernelType DGMassInverse::CGKernels::Fallback(
-   int dim, int, int)
-{
-   if (dim == 1) { return &DGMassInverse::DGMassCGIteration<1>; }
-   else if (dim == 2) { return &DGMassInverse::DGMassCGIteration<2>; }
-   else if (dim == 3) { return &DGMassInverse::DGMassCGIteration<3>; }
-   else { MFEM_ABORT("Unsupported dimension."); }
-}
-
-/// @endcond
-
 } // namespace mfem
@@ -15,6 +15,7 @@
 #include "../linalg/kernels.hpp"
 #include "kernels.hpp"
 #include "integ/bilininteg_mass_kernels.hpp"
+#include "dgmassinv.hpp"

 namespace mfem
 {
@@ -333,6 +334,170 @@ void DGMassBasis(const int e,

 } // namespace internal

+template<int DIM, int D1D, int Q1D>
+void DGMassInverse::DGMassCGIteration(const Vector &b_, Vector &u_) const
+{
+   using namespace internal; // host/device kernel functions
+
+   const int NE = fes.GetNE();
+   const int d1d = m->dofs1D;
+   const int q1d = m->quad1D;
+
+   const int ND = static_cast<int>(pow(d1d, DIM));
+
+   const auto B = m->maps->B.Read();
+   const auto Bt = m->maps->Bt.Read();
+   const auto pa_data = m->pa_data.Read();
+   const auto dinv = diag_inv.Read();
+   auto r = r_.Write();
+   auto d = d_.Write();
+   auto z = z_.Write();
+   auto u = u_.ReadWrite();
+
+   const real_t RELTOL = rel_tol;
+   const real_t ABSTOL = abs_tol;
+   const int MAXIT = max_iter;
+   const bool IT_MODE = iterative_mode;
+   const bool CHANGE_BASIS = (d2q != nullptr);
+
+   // b is the right-hand side (if no change of basis, this just points to the
+   // incoming RHS vector, if we have to change basis, this points to the
+   // internal b2 vector where we put the transformed RHS)
+   const real_t *b;
+   // the following are non-null if we have to change basis
+   real_t *b2 = nullptr; // non-const access to b2
+   const real_t *b_orig = nullptr; // RHS vector in "original" basis
+   const real_t *d2q_B = nullptr; // matrix to transform initial guess
+   const real_t *q2d_B = nullptr; // matrix to transform solution
+   const real_t *q2d_Bt = nullptr; // matrix to transform RHS
+   if (CHANGE_BASIS)
+   {
+      d2q_B = d2q->B.Read();
+      q2d_B = B_.Read();
+      q2d_Bt = Bt_.Read();
+
+      b2 = b2_.Write();
+      b_orig = b_.Read();
+      b = b2;
+   }
+   else
+   {
+      b = b_.Read();
+   }
+
+   static constexpr int NB = Q1D ? Q1D : 1; // block size
+
+   mfem::forall_2D<NB*NB>(NE, NB, NB, [=] MFEM_HOST_DEVICE (int e)
+   {
+      // Perform change of basis if needed
+      if (CHANGE_BASIS)
+      {
+         // Transform RHS
+         DGMassBasis<DIM,D1D>(e, NE, q2d_Bt, b_orig, b2, d1d);
+         if (IT_MODE)
+         {
+            // Transform initial guess
+            DGMassBasis<DIM,D1D>(e, NE, d2q_B, u, u, d1d);
+         }
+      }
+
+      const int tid = MFEM_THREAD_ID(x) + NB*MFEM_THREAD_ID(y);
+
+      // Compute first residual
+      if (IT_MODE)
+      {
+         DGMassApply<DIM,D1D,Q1D>(e, NE, B, Bt, pa_data, u, r, d1d, q1d);
+         DGMassAxpy(e, NE, ND, 1.0, b, -1.0, r, r); // r = b - r
+      }
+      else
+      {
+         // if not in iterative mode, use zero initial guess
+         const int BX = MFEM_THREAD_SIZE(x);
+         const int BY = MFEM_THREAD_SIZE(y);
+         const int bxy = BX*BY;
+         const auto B = ConstDeviceMatrix(b, ND, NE);
+         auto U = DeviceMatrix(u, ND, NE);
+         auto R = DeviceMatrix(r, ND, NE);
+         for (int i = tid; i < ND; i += bxy)
+         {
+            U(i, e) = 0.0;
+            R(i, e) = B(i, e);
+         }
+         MFEM_SYNC_THREAD;
+      }
+
+      DGMassPreconditioner(e, NE, ND, dinv, r, z);
+      DGMassAxpy(e, NE, ND, 1.0, z, 0.0, z, d); // d = z
+
+      real_t nom = DGMassDot<NB>(e, NE, ND, d, r);
+      if (nom < 0.0) { return; /* Not positive definite */ }
+      real_t r0 = fmax(nom*RELTOL*RELTOL, ABSTOL*ABSTOL);
+      if (nom <= r0) { return; /* Converged */ }
+
+      DGMassApply<DIM,D1D,Q1D>(e, NE, B, Bt, pa_data, d, z, d1d, q1d);
+      real_t den = DGMassDot<NB>(e, NE, ND, z, d);
+      if (den <= 0.0)
+      {
+         DGMassDot<NB>(e, NE, ND, d, d);
+         // d2 > 0 => not positive definite
+         if (den == 0.0) { return; }
+      }
+
+      // start iteration
+      int i = 1;
+      while (true)
+      {
+         const real_t alpha = nom/den;
+         DGMassAxpy(e, NE, ND, 1.0, u, alpha, d, u); // u = u + alpha*d
+         DGMassAxpy(e, NE, ND, 1.0, r, -alpha, z, r); // r = r - alpha*A*d
+
+         DGMassPreconditioner(e, NE, ND, dinv, r, z);
+
+         real_t betanom = DGMassDot<NB>(e, NE, ND, r, z);
+         if (betanom < 0.0) { return; /* Not positive definite */ }
+         if (betanom <= r0) { break; /* Converged */ }
+
+         if (++i > MAXIT) { break; }
+
+         const real_t beta = betanom/nom;
+         DGMassAxpy(e, NE, ND, 1.0, z, beta, d, d); // d = z + beta*d
+         DGMassApply<DIM,D1D,Q1D>(e, NE, B, Bt, pa_data, d, z, d1d, q1d); // z = A d
+         den = DGMassDot<NB>(e, NE, ND, d, z);
+         if (den <= 0.0)
+         {
+            DGMassDot<NB>(e, NE, ND, d, d);
+            // d2 > 0 => not positive definite
+            if (den == 0.0) { break; }
+         }
+         nom = betanom;
+      }
+
+      if (CHANGE_BASIS)
+      {
+         DGMassBasis<DIM,D1D>(e, NE, q2d_B, u, u, d1d);
+      }
+   });
+}
+
+/// @cond Suppress_Doxygen_warnings
+
+template <int DIM, int D1D, int Q1D>
+inline DGMassInverse::CGKernelType DGMassInverse::CGKernels::Kernel()
+{
+   return &DGMassInverse::DGMassCGIteration<DIM,D1D,Q1D>;
+}
+
+inline DGMassInverse::CGKernelType DGMassInverse::CGKernels::Fallback(
+   int dim, int, int)
+{
+   if (dim == 1) { return &DGMassInverse::DGMassCGIteration<1>; }
+   else if (dim == 2) { return &DGMassInverse::DGMassCGIteration<2>; }
+   else if (dim == 3) { return &DGMassInverse::DGMassCGIteration<3>; }
+   else { MFEM_ABORT("Unsupported dimension."); }
+}
+
+/// @endcond
+
 } // namespace mfem

 #endif
@@ -69,9 +69,9 @@ inline int ToLexOrdering2D(const int face_id, const int size1d, const int i)
 }

 /// @brief Given a face DOF index on a shared face, ordered lexicographically
-/// relative to element the element (where the local face is face_id), and
-/// return the corresponding face DOF index ordered lexicographically relative
-/// to the face itself.
+/// relative to the element (where the local face is face_id), return the
+/// corresponding face DOF index ordered lexicographically relative to the face
+/// itself.
 MFEM_HOST_DEVICE
 inline int PermuteFace2D(const int face_id, const int orientation,
                         const int size1d, const int index)
@@ -231,7 +231,7 @@ void FiniteElement::CalcPhysLaplacian(ElementTransformation &Trans,
   {
      for (int nd = 0; nd < dof; nd++)
      {
-         Laplacian[nd] = hess(nd,0) + hess(nd,4) + hess(nd,5);
+         Laplacian[nd] = hess(nd,0) + hess(nd,3) + hess(nd,5);
      }
   }
   else if (dim == 2)
@@ -268,11 +268,9 @@ void FiniteElement::CalcPhysLinLaplacian(ElementTransformation &Trans,
      scale[0] =   Gij(0,0);
      scale[1] = 2*Gij(0,1);
      scale[2] = 2*Gij(0,2);
-
-      scale[3] = 2*Gij(1,2);
-      scale[4] =   Gij(2,2);
-
-      scale[5] =   Gij(1,1);
+      scale[3] =   Gij(1,1);
+      scale[4] = 2*Gij(1,2);
+      scale[5] =   Gij(2,2);
   }
   else if (dim == 2)
   {
@@ -309,12 +307,12 @@ void  FiniteElement::CalcPhysHessian(ElementTransformation &Trans,
      map[2] = 2;

      map[3] = 1;
-      map[4] = 5;
-      map[5] = 3;
+      map[4] = 3;
+      map[5] = 4;

      map[6] = 2;
-      map[7] = 3;
-      map[8] = 4;
+      map[7] = 4;
+      map[8] = 5;
   }
   else if (dim == 2)
   {
@@ -382,11 +380,7 @@ const DofToQuad &FiniteElement::GetDofToQuad(const IntegrationRule &ir,
   #pragma omp critical (DofToQuad)
 #endif
   {
-      for (int i = 0; i < dof2quad_array.Size(); i++)
-      {
-         d2q = dof2quad_array[i];
-         if (d2q->IntRule != &ir || d2q->mode != mode) { d2q = nullptr; }
-      }
+      d2q = DofToQuad::SearchArray(dof2quad_array, ir, mode);
      if (!d2q)
      {
 #ifdef MFEM_THREAD_SAFE
@@ -661,58 +655,67 @@ void ScalarFiniteElement::ScalarLocalL2Restriction(
 void NodalFiniteElement::CreateLexicographicFullMap(const IntegrationRule &ir)
 const
 {
+   // Get the FULL version of the map. This call contains omp critical region,
+   // so it is done before the critical region below.
+   auto &d2q = GetDofToQuad(ir, DofToQuad::FULL);

 #if defined(MFEM_THREAD_SAFE) && defined(MFEM_USE_OPENMP)
   #pragma omp critical (DofToQuad)
 #endif
   {
-      // Get the FULL version of the map.
-      auto &d2q = GetDofToQuad(ir, DofToQuad::FULL);
-      //Undo the native ordering which is what FiniteElement::GetDofToQuad returns.
-      auto *d2q_new = new DofToQuad(d2q);
-      d2q_new->mode = DofToQuad::LEXICOGRAPHIC_FULL;
-      const int nqpt = ir.GetNPoints();
-
-      const int b_dim = (range_type == VECTOR) ? dim : 1;
-
-      for (int i = 0; i < nqpt; i++)
+      // Do not run if the new Dof2Quad is already present, e.g. added in a
+      // previous call or added by another omp thread.
+      if (DofToQuad::SearchArray(dof2quad_array, ir,
+                                 DofToQuad::LEXICOGRAPHIC_FULL) == nullptr)
      {
-         for (int d = 0; d < b_dim; d++)
+         // Undo the native ordering which is what FiniteElement::GetDofToQuad
+         // returns.
+         auto *d2q_new = new DofToQuad(d2q);
+         d2q_new->mode = DofToQuad::LEXICOGRAPHIC_FULL;
+         const int nqpt = ir.GetNPoints();
+
+         const int b_dim = (range_type == VECTOR) ? dim : 1;
+
+         for (int i = 0; i < nqpt; i++)
         {
-            for (int j = 0; j < dof; j++)
+            for (int d = 0; d < b_dim; d++)
            {
-               const double val = d2q.B[i + nqpt*(d+b_dim*lex_ordering[j])];
-               d2q_new->B[i+nqpt*(d+b_dim*j)] = val;
-               d2q_new->Bt[j+dof*(i+nqpt*d)] = val;
+               for (int j = 0; j < dof; j++)
+               {
+                  const double val = d2q.B[i + nqpt*(d+b_dim*lex_ordering[j])];
+                  d2q_new->B[i+nqpt*(d+b_dim*j)] = val;
+                  d2q_new->Bt[j+dof*(i+nqpt*d)] = val;
+               }
            }
         }
-      }

-      const int g_dim = [this]()
-      {
-         switch (deriv_type)
+         const int g_dim = [this]()
         {
-            case GRAD: return dim;
-            case DIV: return 1;
-            case CURL: return cdim;
-            default: return 0;
-         }
-      }();
-
-      for (int i = 0; i < nqpt; i++)
-      {
-         for (int d = 0; d < g_dim; d++)
-         {
-            for (int j = 0; j < dof; j++)
+            switch (deriv_type)
            {
-               const double val = d2q.G[i + nqpt*(d+g_dim*lex_ordering[j])];
-               d2q_new->G[i+nqpt*(d+g_dim*j)] = val;
-               d2q_new->Gt[j+dof*(i+nqpt*d)] = val;
+               case GRAD: return dim;
+               case DIV: return 1;
+               case CURL: return cdim;
+               default: return 0;
+            }
+         }();
+
+         for (int i = 0; i < nqpt; i++)
+         {
+            for (int d = 0; d < g_dim; d++)
+            {
+               for (int j = 0; j < dof; j++)
+               {
+                  const double val = d2q.G[i + nqpt*(d+g_dim*lex_ordering[j])];
+                  d2q_new->G[i+nqpt*(d+g_dim*j)] = val;
+                  d2q_new->Gt[j+dof*(i+nqpt*d)] = val;
+               }
            }
         }
+
+         dof2quad_array.Append(d2q_new);
      }

-      dof2quad_array.Append(d2q_new);
   }
 }

@@ -724,13 +727,7 @@ const DofToQuad &NodalFiniteElement::GetDofToQuad(const IntegrationRule &ir,
   #pragma omp critical (DofToQuad)
 #endif
   {
-      //Should make this loop a function of FiniteElement
-      for (int i = 0; i < dof2quad_array.Size(); i++)
-      {
-         d2q = dof2quad_array[i];
-         if (d2q->IntRule == &ir && d2q->mode == mode) { break; }
-         d2q = nullptr;
-      }
+      d2q = DofToQuad::SearchArray(dof2quad_array, ir, mode);
   }
   if (d2q) { return *d2q; }
   if (mode != DofToQuad::LEXICOGRAPHIC_FULL)
@@ -1047,9 +1044,50 @@ void VectorFiniteElement::SetDerivMembers()
   switch (map_type)
   {
      case H_DIV:
-         deriv_type = DIV;
-         deriv_range_type = SCALAR;
-         deriv_map_type = INTEGRAL;
+         switch (dim)
+         {
+            case 3: // div: 3D H_DIV -> 3D INTEGRAL
+               deriv_type = DIV;
+               deriv_range_type = SCALAR;
+               deriv_map_type = INTEGRAL;
+               break;
+            case 2: // div: 2D H_DIV -> 2D INTEGRAL
+               deriv_type = DIV;
+               deriv_range_type = SCALAR;
+               deriv_map_type = INTEGRAL;
+               break;
+            default:
+               MFEM_ABORT("Invalid dimension, Dim = " << dim);
+         }
+         break;
+      case H_DIV_R2D:
+         switch (dim)
+         {
+            case 2: // div: 2D H_DIV_R2D -> 2D INTEGRAL
+               deriv_type = DIV;
+               deriv_range_type = SCALAR;
+               deriv_map_type = INTEGRAL;
+               break;
+            case 1: // div: 1D H_DIV_R2D -> 1D INTEGRAL
+               deriv_type = DIV;
+               deriv_range_type = SCALAR;
+               deriv_map_type = INTEGRAL;
+               break;
+            default:
+               MFEM_ABORT("Invalid dimension, Dim = " << dim);
+         }
+         break;
+      case H_DIV_R1D:
+         switch (dim)
+         {
+            case 1: // div: 1D H_DIV_R1D -> 1D INTEGRAL
+               deriv_type = DIV;
+               deriv_range_type = SCALAR;
+               deriv_map_type = INTEGRAL;
+               break;
+            default:
+               MFEM_ABORT("Invalid dimension, Dim = " << dim);
+         }
         break;
      case H_CURL:
         switch (dim)
@@ -1067,13 +1105,49 @@ void VectorFiniteElement::SetDerivMembers()
               break;
            case 1:
               deriv_type = NONE;
-               deriv_range_type = SCALAR;
-               deriv_map_type = INTEGRAL;
+               deriv_range_type = UNKNOWN_RANGE_TYPE;
+               deriv_map_type = UNKNOWN_MAP_TYPE;
               break;
            default:
               MFEM_ABORT("Invalid dimension, Dim = " << dim);
         }
         break;
+      case H_CURL_R2D:
+         switch (dim)
+         {
+            case 2:
+               // curl: 2D H_CURL_R2D -> H_DIV_R2D
+               deriv_type = CURL;
+               deriv_range_type = VECTOR;
+               deriv_map_type = H_DIV_R2D;
+               break;
+            case 1:
+               // curl: 1D H_CURL_R2D -> H_DIV_R2D
+               deriv_type = CURL;
+               deriv_range_type = VECTOR;
+               deriv_map_type = H_DIV_R2D;
+               break;
+            default:
+               MFEM_ABORT("Invalid dimension, Dim = " << dim);
+         }
+         break;
+      case H_CURL_R1D:
+         switch (dim)
+         {
+            case 1:
+               // curl: 1D H_CURL_R1D -> H_DIV_R1D
+               deriv_type = CURL;
+               deriv_range_type = VECTOR;
+               deriv_map_type = H_DIV_R1D;
+               break;
+            case 0:
+               deriv_type = NONE;
+               deriv_range_type = UNKNOWN_RANGE_TYPE;
+               deriv_map_type = UNKNOWN_MAP_TYPE;
+            default:
+               MFEM_ABORT("Invalid dimension, Dim = " << dim);
+         }
+         break;
      default:
         MFEM_ABORT("Invalid MapType = " << map_type);
   }
@@ -2631,15 +2705,7 @@ const DofToQuad &TensorBasisElement::GetTensorDofToQuad(
   #pragma omp critical (DofToQuad)
 #endif
   {
-      for (int i = 0; i < dof2quad_array.Size(); i++)
-      {
-         auto* d2q_ = dof2quad_array[i];
-         if (d2q_->IntRule == &ir && d2q_->mode == mode)
-         {
-            d2q = d2q_;
-            break;
-         }
-      }
+      d2q = DofToQuad::SearchArray(dof2quad_array, ir, mode);
      if (!d2q)
      {
         d2q = new DofToQuad;
@@ -44,7 +44,7 @@ public:
      NumBasisTypes   = 9   /**< Keep track of maximum types to prevent
                                 hard-coding */
   };
-   /** @brief If the input does not represents a valid BasisType, abort with an
+   /** @brief If the input does not represent a valid BasisType, abort with an
       error; otherwise return the input. */
   static int Check(int b_type)
   {
@@ -52,7 +52,7 @@ public:
                  "unknown BasisType: " << b_type);
      return b_type;
   }
-   /** @brief If the input does not represents a valid nodal BasisType, abort
+   /** @brief If the input does not represent a valid nodal BasisType, abort
       with an error; otherwise return the input. */
   static int CheckNodal(int b_type)
   {
@@ -222,6 +222,12 @@ public:

   /// Returns absolute value of the maps
   DofToQuad Abs() const;
+
+   /// Auxiliary function for searching DofToQuad arrays.
+   static inline DofToQuad *SearchArray(
+      const Array<DofToQuad*> &dof2quad_array,
+      const IntegrationRule &ir,
+      DofToQuad::Mode mode);
 };

 /// Describes the function space on each element
@@ -289,10 +295,20 @@ public:
                          $ u(x) = (1/w) \hat u(\hat x) $ */
      H_DIV,     /**< For vector fields; preserves surface integrals of the
                          normal component $ u(x) = (J/w) \hat u(\hat x) $ */
-      H_CURL     /**< For vector fields; preserves line integrals of the
+      H_CURL,    /**< For vector fields; preserves line integrals of the
                          tangential component
                          $ u(x) = J^{-t} \hat u(\hat x) $ (square J),
                          $ u(x) = J(J^t J)^{-1} \hat u(\hat x) $ (general J) */
+      H_DIV_R2D, /**< For 3-component vector fields in 2D; equivalent to a
+                          direct sum of an H_DIV basis and an INTEGRAL basis */
+      H_CURL_R2D,/**< For 3-component vector fields in 2D; equivalent to a
+                          direct sum of an H_CURL basis and a VALUE basis */
+      H_DIV_R1D, /**< For 3-component vector fields in 1D; equivalent to a
+                          direct sum of a VALUE basis and a pair of INTEGRAL
+                          bases */
+      H_CURL_R1D /**< For 3-component vector fields in 1D; equivalent to a
+                          direct sum of an INTEGRAL basis and a pair of VALUE
+                          bases */
   };

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

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

-   /// Returns the dimension of the curl for vector-valued finite elements.
+   /** @brief Returns the vector dimension, in physical space, for
+       vector-valued finite elements, which is also the width of the
+       DenseMatrix argument in
+       CalcPhysVShape(ElementTransformation &Trans, DenseMatrix &shape). */
+   virtual int GetPhysRangeDim(int /* space_dim */) const { return vdim; }
+
+   /** Returns the dimension of the curl for vector-valued finite elements,
+       which is also the width of the DenseMatrix argument in
+       CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape). */
   int GetCurlDim() const { return cdim; }

+   /** Returns the dimension, in physical space, of the curl for vector-valued
+       finite elements, which is also the width of the DenseMatrix argument in
+       CalcPhysCurlShape(ElementTransformation &Trans, DenseMatrix &curl_shape).
+   */
+   virtual int GetPhysCurlDim(int /* space_dim */) const { return cdim; }
+
   /// Returns the Geometry::Type of the reference element.
   Geometry::Type GetGeomType() const { return geom_type; }

@@ -407,6 +439,7 @@ public:
   /** Each row of the result DenseMatrix @a Hessian contains upper triangular
       part of the Hessian of one shape function.
       The order in 2D is {u_xx, u_xy, u_yy}.
+       The order in 3D is {u_xx, u_xy, u_xz, u_yy, u_yz, u_zz}.
       The size (#dof x (#dim (#dim+1)/2) of @a Hessian must be set in advance.*/
   virtual void CalcHessian(const IntegrationPoint &ip,
                            DenseMatrix &Hessian) const;
@@ -983,6 +1016,8 @@ protected:
 public:
   VectorFiniteElement(int D, Geometry::Type G, int Do, int O, int M,
                       int F = FunctionSpace::Pk);
+
+   int GetPhysRangeDim(int space_dim) const override { return space_dim; }
 };

 /// @brief Class for computing 1D special polynomials and their associated basis
@@ -1120,7 +1155,7 @@ public:
      return GetPoints(p, btype, on_device);
   }

-   /// Get coordinates of a closed (GaussLegendre) set of points if degree @a p
+   /// Get coordinates of a closed (GaussLobatto) set of points if degree @a p
   const real_t *ClosedPoints(const int p,
                              const int btype = BasisType::GaussLobatto,
                              bool on_device = false)
@@ -1376,6 +1411,21 @@ public:
 void InvertLinearTrans(ElementTransformation &trans,
                       const IntegrationPoint &pt, Vector &x);

+
+// static inline method
+inline DofToQuad *DofToQuad::SearchArray(
+   const Array<DofToQuad*> &dof2quad_array,
+   const IntegrationRule &ir,
+   DofToQuad::Mode mode)
+{
+   for (int i = 0; i < dof2quad_array.Size(); i++)
+   {
+      DofToQuad *d2q = dof2quad_array[i];
+      if (d2q->IntRule == &ir && d2q->mode == mode) { return d2q; }
+   }
+   return nullptr;
+}
+
 } // namespace mfem

 #endif
@@ -60,6 +60,12 @@ void Linear1DFiniteElement::CalcDShape(const IntegrationPoint &ip,
   dshape(1,0) =  1.;
 }

+void Linear1DFiniteElement::CalcHessian(const IntegrationPoint &ip,
+                                        DenseMatrix &h) const
+{
+   h = 0.0;
+}
+
 Linear2DFiniteElement::Linear2DFiniteElement()
   : NodalFiniteElement(2, Geometry::TRIANGLE, 3, 1)
 {
@@ -87,6 +93,11 @@ void Linear2DFiniteElement::CalcDShape(const IntegrationPoint &ip,
   dshape(2,0) =  0.; dshape(2,1) =  1.;
 }

+void Linear2DFiniteElement::CalcHessian(const IntegrationPoint &ip,
+                                        DenseMatrix &h) const
+{
+   h = 0.0;
+}

 BiLinear2DFiniteElement::BiLinear2DFiniteElement()
   : NodalFiniteElement(2, Geometry::SQUARE, 4, 1, FunctionSpace::Qk)
@@ -1256,6 +1267,12 @@ void Linear3DFiniteElement::CalcDShape(const IntegrationPoint &ip,
   }
 }

+void Linear3DFiniteElement::CalcHessian(const IntegrationPoint &ip,
+                                        DenseMatrix &h) const
+{
+   h = 0.0;
+}
+
 void Linear3DFiniteElement::GetFaceDofs (int face, int **dofs, int *ndofs)
 const
 {
@@ -1632,6 +1649,37 @@ void TriLinear3DFiniteElement::CalcDShape(const IntegrationPoint &ip,
   dshape(7,2) =   ox *  y;
 }

+void TriLinear3DFiniteElement::CalcHessian(const IntegrationPoint &ip,
+                                           DenseMatrix &h) const
+{
+   real_t x = ip.x, y = ip.y, z = ip.z;
+   real_t ox = 1.-x, oy = 1.-y, oz = 1.-z;
+
+   h(0,0) = 0.;   h(0,1) =  oz;   h(0,2) = oy;
+   h(0,3) = 0.;   h(0,4) =  ox;   h(0,5) = 0.;
+
+   h(1,0) = 0.;   h(1,1) = -oz;   h(1,2) = -oy;
+   h(1,3) = 0.;   h(1,4) =  x;    h(1,5) = 0.;
+
+   h(2,0) = 0.;   h(2,1) =  oz;   h(2,2) = -y;
+   h(2,3) = 0.;   h(2,4) =  -x;   h(2,5) = 0.;
+
+   h(3,0) = 0.;   h(3,1) = -oz;   h(3,2) = y;
+   h(3,3) = 0.;   h(3,4) = -ox;   h(3,5) = 0.;
+
+   h(4,0) = 0.;   h(4,1) =  z;    h(4,2) = -oy;
+   h(4,3) = 0.;   h(4,4) = -ox;   h(4,5) = 0.;
+
+   h(5,0) = 0.;   h(5,1) = -z;    h(5,2) = oy;
+   h(5,3) = 0.;   h(5,4) = -x;    h(5,5) = 0.;
+
+   h(6,0) = 0.;   h(6,1) =  z;    h(6,2) = y;
+   h(6,3) = 0.;   h(6,4) =  x;    h(6,5) = 0.;
+
+   h(7,0) = 0.;   h(7,1) = -z;    h(7,2) = -y;
+   h(7,3) = 0.;   h(7,4) = ox;    h(7,5) = 0.;
+}
+

 P0SegmentFiniteElement::P0SegmentFiniteElement(int Ord)
   : NodalFiniteElement(1, Geometry::SEGMENT, 1, Ord)   // default Ord = 0
@@ -50,6 +50,8 @@ public:
       contains the derivative of one shape function */
   void CalcDShape(const IntegrationPoint &ip,
                   DenseMatrix &dshape) const override;
+   void CalcHessian(const IntegrationPoint &ip,
+                    DenseMatrix &h) const override;
 };

 /// A 2D linear element on triangle with nodes at the vertices of the triangle
@@ -70,6 +72,8 @@ public:
       so that each row contains the derivatives of one shape function */
   void CalcDShape(const IntegrationPoint &ip,
                   DenseMatrix &dshape) const override;
+   void CalcHessian(const IntegrationPoint &ip,
+                    DenseMatrix &h) const override;
   void ProjectDelta(int vertex, Vector &dofs) const override
   { dofs = 0.0; dofs(vertex) = 1.0; }
 };
@@ -404,6 +408,9 @@ public:
   void CalcDShape(const IntegrationPoint &ip,
                   DenseMatrix &dshape) const override;

+   void CalcHessian(const IntegrationPoint &ip,
+                    DenseMatrix &h) const override;
+
   void ProjectDelta(int vertex, Vector &dofs) const override
   { dofs = 0.0; dofs(vertex) = 1.0; }

@@ -445,7 +452,8 @@ public:
       so that each row contains the derivatives of one shape function */
   void CalcDShape(const IntegrationPoint &ip,
                   DenseMatrix &dshape) const override;
-
+   void CalcHessian(const IntegrationPoint &ip,
+                    DenseMatrix &h) const override;
   void ProjectDelta(int vertex, Vector &dofs) const override
   { dofs = 0.0; dofs(vertex) = 1.0; }
 };
@@ -589,7 +589,7 @@ void H1_TriangleElement::CalcHessian(const IntegrationPoint &ip,
   Vector   shape_x(p + 1),   shape_y(p + 1),   shape_l(p + 1);
   Vector  dshape_x(p + 1),  dshape_y(p + 1),  dshape_l(p + 1);
   Vector ddshape_x(p + 1), ddshape_y(p + 1), ddshape_l(p + 1);
-   DenseMatrix ddu(dof, dim);
+   DenseMatrix ddu(dof, (dim*(dim+1))/2);
 #endif

   poly1d.CalcBasis(p, ip.x, shape_x, dshape_x, ddshape_x);
@@ -2531,7 +2531,7 @@ void ND_FuentesPyramidElement::calcCurlBasis(const int p,

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

+   int GetPhysRangeDim(int space_dim) const override { return 2; }
+   int GetPhysCurlDim(int space_dim) const override { return 1; }
+
   void CalcVShape(const IntegrationPoint &ip,
                   DenseMatrix &shape) const override;

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

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

@@ -84,6 +84,46 @@ void NURBS1DFiniteElement::CalcHessian (const IntegrationPoint &ip,
   add(1.0, hess, (-d2sum + 2*dsum*dsum*sum)*sum*sum, shape_x, hess);
 }

+void NURBS1DFiniteElement::Project(Coefficient &coeff,
+                                   ElementTransformation &Trans,
+                                   Vector &dofs) const
+{
+   IntegrationPoint ip;
+
+   for (int i = 0; i <= order; i++)
+   {
+      real_t kx = kv[0]->GetBotella(ijk[0] + i);
+      if (!kv[0]->inSpan(kx, ijk[0]+order)) { continue; }
+      ip.x = kv[0]->GetRefPoint(kx, ijk[0]+order);
+
+      Trans.SetIntPoint(&ip);
+      dofs(i) = coeff.Eval(Trans, ip);
+   }
+}
+
+void NURBS1DFiniteElement::Project(VectorCoefficient &vc,
+                                   ElementTransformation &Trans,
+                                   Vector &dofs) const
+{
+   MFEM_ASSERT(dofs.Size() == vc.GetVDim()*dof, "");
+   Vector x(vc.GetVDim());
+   IntegrationPoint ip;
+
+   for (int i = 0; i <= order; i++)
+   {
+      real_t kx = kv[0]->GetBotella(ijk[0] + i);
+      if (!kv[0]->inSpan(kx, ijk[0]+order)) { continue; }
+      ip.x = kv[0]->GetRefPoint(kx, ijk[0]+order);
+
+      Trans.SetIntPoint(&ip);
+      vc.Eval(x, Trans, ip);
+      for (int j = 0; j < x.Size(); j++)
+      {
+         dofs(dof*j+i) = x(j);
+      }
+   }
+}
+

 void NURBS2DFiniteElement::SetOrder() const
 {
@@ -215,6 +255,63 @@ void NURBS2DFiniteElement::CalcHessian (const IntegrationPoint &ip,
   }
 }

+void NURBS2DFiniteElement::Project(Coefficient &coeff,
+                                   ElementTransformation &Trans,
+                                   Vector &dofs) const
+{
+   IntegrationPoint ip;
+   for (int o = 0, j = 0; j <= orders[1]; j++)
+   {
+      real_t ky = kv[1]->GetBotella(ijk[1] + j);
+      if (!kv[1]->inSpan(ky, ijk[1]+orders[1]))
+      {
+         o += orders[0] + 1;
+         continue;
+      }
+      ip.y = kv[1]->GetRefPoint(ky, ijk[1]+orders[1]);
+      for (int i = 0; i <= orders[0]; i++, o++)
+      {
+         real_t kx = kv[0]->GetBotella(ijk[0] + i);
+         if (!kv[0]->inSpan(kx, ijk[0]+orders[0])) { continue; }
+         ip.x = kv[0]->GetRefPoint(kx, ijk[0]+orders[0]);
+
+         Trans.SetIntPoint(&ip);
+         dofs(o) = coeff.Eval(Trans, ip);
+      }
+   }
+}
+
+void NURBS2DFiniteElement::Project(VectorCoefficient &vc,
+                                   ElementTransformation &Trans,
+                                   Vector &dofs) const
+{
+   MFEM_ASSERT(dofs.Size() == vc.GetVDim()*dof, "");
+   Vector x(vc.GetVDim());
+   IntegrationPoint ip;
+   for (int o = 0, j = 0; j <= orders[1]; j++)
+   {
+      real_t ky = kv[1]->GetBotella(ijk[1] + j);
+      if (!kv[1]->inSpan(ky, ijk[1]+orders[1]))
+      {
+         o += orders[0] + 1;
+         continue;
+      }
+      ip.y = kv[1]->GetRefPoint(ky, ijk[1]+orders[1]);
+      for (int i = 0; i <= orders[0]; i++, o++)
+      {
+         real_t kx = kv[0]->GetBotella(ijk[0] + i);
+         if (!kv[0]->inSpan(kx, ijk[0]+orders[0])) { continue; }
+         ip.x = kv[0]->GetRefPoint(kx, ijk[0]+orders[0]);
+
+         Trans.SetIntPoint(&ip);
+         vc.Eval(x, Trans, ip);
+         for (int v = 0; v < x.Size(); v++)
+         {
+            dofs(dof*v+o) = x(v);
+         }
+      }
+   }
+}

 void NURBS3DFiniteElement::SetOrder() const
 {
@@ -348,11 +445,10 @@ void NURBS3DFiniteElement::CalcHessian (const IntegrationPoint &ip,
            d2sum[0] += ( hessian(o,0) = d2sx*sy*sz*weights(o) );
            d2sum[1] += ( hessian(o,1) = dsx*dsy*sz*weights(o) );
            d2sum[2] += ( hessian(o,2) = dsx*sy*dsz*weights(o) );
+            d2sum[3] += ( hessian(o,3) = sx*d2sy*sz*weights(o) );
+            d2sum[4] += ( hessian(o,4) = sx*dsy*dsz*weights(o) );
+            d2sum[5] += ( hessian(o,5) = sx*sy*d2sz*weights(o) );

-            d2sum[3] += ( hessian(o,3) = sx*dsy*dsz*weights(o) );
-
-            d2sum[4] += ( hessian(o,4) = sx*sy*d2sz*weights(o) );
-            d2sum[5] += ( hessian(o,5) = sx*d2sy*sz*weights(o) );
         }
      }
   }
@@ -401,6 +497,85 @@ void NURBS3DFiniteElement::CalcHessian (const IntegrationPoint &ip,
   }
 }

+void NURBS3DFiniteElement::Project(Coefficient &coeff,
+                                   ElementTransformation &Trans,
+                                   Vector &dofs) const
+{
+   IntegrationPoint ip;
+
+   for (int o = 0, k = 0; k <= orders[2]; k++)
+   {
+      real_t kz = kv[2]->GetBotella(ijk[2] + k);
+      if (!kv[2]->inSpan(kz, ijk[2]+orders[2]))
+      {
+         o += (orders[0] + 1)*(orders[1] + 1);
+         continue;
+      }
+      ip.z = kv[2]->GetRefPoint(kz, ijk[2]+orders[2]);
+      for (int j = 0; j <= orders[1]; j++)
+      {
+         real_t ky = kv[1]->GetBotella(ijk[1] + j);
+         if (!kv[1]->inSpan(ky, ijk[1]+orders[1]))
+         {
+            o += orders[0] + 1;
+            continue;
+         }
+         ip.y = kv[1]->GetRefPoint(ky, ijk[1]+orders[1]);
+         for (int i = 0; i <= orders[0]; i++, o++)
+         {
+            real_t kx = kv[0]->GetBotella(ijk[0] + i);
+            if (!kv[0]->inSpan(kx, ijk[0]+orders[0])) { continue; }
+            ip.x = kv[0]->GetRefPoint(kx, ijk[0]+orders[0]);
+
+            Trans.SetIntPoint(&ip);
+            dofs(o) = coeff.Eval(Trans, ip);
+         }
+      }
+   }
+}
+
+void NURBS3DFiniteElement::Project(VectorCoefficient &vc,
+                                   ElementTransformation &Trans,
+                                   Vector &dofs) const
+{
+   MFEM_ASSERT(dofs.Size() == vc.GetVDim()*dof, "");
+   Vector x(vc.GetVDim());
+   IntegrationPoint ip;
+
+   for (int o = 0, k = 0; k <= orders[2]; k++)
+   {
+      real_t kz = kv[2]->GetBotella(ijk[2] + k);
+      if (!kv[2]->inSpan(kz, ijk[2]+orders[2]))
+      {
+         o += (orders[0] + 1)*(orders[1] + 1);
+         continue;
+      }
+      ip.z = kv[2]->GetRefPoint(kz, ijk[2]+orders[2]);
+      for (int j = 0; j <= orders[1]; j++)
+      {
+         real_t ky = kv[1]->GetBotella(ijk[1] + j);
+         if (!kv[1]->inSpan(ky, ijk[1]+orders[1]))
+         {
+            o += orders[0] + 1;
+            continue;
+         }
+         ip.y = kv[1]->GetRefPoint(ky, ijk[1]+orders[1]);
+         for (int i = 0; i <= orders[0]; i++, o++)
+         {
+            real_t kx = kv[0]->GetBotella(ijk[0] + i);
+            if (!kv[0]->inSpan(kx, ijk[0]+orders[0])) { continue; }
+            ip.x = kv[0]->GetRefPoint(kx, ijk[0]+orders[0]);
+
+            Trans.SetIntPoint(&ip);
+            vc.Eval(x, Trans, ip);
+            for (int v = 0; v < x.Size(); v++)
+            {
+               dofs(dof*v+o) = x(v);
+            }
+         }
+      }
+   }
+}

 void NURBS_HDiv2DFiniteElement::SetOrder() const
 {
@@ -517,6 +692,63 @@ void NURBS_HDiv2DFiniteElement::CalcDivShape(const IntegrationPoint &ip,
   }
 }

+void NURBS_HDiv2DFiniteElement::Project(VectorCoefficient &vc,
+                                        ElementTransformation &Trans,
+                                        Vector &dofs) const
+{
+   MFEM_ASSERT(dofs.Size() == dof, "");
+   MFEM_ASSERT(vc.GetVDim() == 2, "");
+   Vector x(2), mx(2);
+   IntegrationPoint ip;
+   int o = 0;
+
+   for (int j = 0; j <= orders[1]; j++)
+   {
+      real_t ky = kv[1]->GetBotella(ijk[1] + j);
+      if (!kv[1]->inSpan(ky, ijk[1]+orders[1]))
+      {
+         o += orders[0] + 2;
+         continue;
+      }
+      ip.y = kv[1]->GetRefPoint(ky, ijk[1]+orders[1]);
+      for (int i = 0; i <= orders[0]+1; i++, o++)
+      {
+         real_t kx = kv1[0]->GetBotella(ijk[0] + i);
+         if (!kv1[0]->inSpan(kx, ijk[0]+orders[0]+1)) { continue; }
+         ip.x = kv1[0]->GetRefPoint(kx, ijk[0]+orders[0]+1);
+
+         Trans.SetIntPoint(&ip);
+         vc.Eval(x, Trans, ip);
+
+         Trans.AdjugateJacobian().Mult(x,mx);
+         dofs(o) = mx(0);
+      }
+   }
+
+   for (int j = 0; j <= orders[1]+1; j++)
+   {
+      real_t ky = kv1[1]->GetBotella(ijk[1] + j);
+      if (!kv1[1]->inSpan(ky, ijk[1]+orders[1]+1))
+      {
+         o += orders[0] + 1;
+         continue;
+      }
+      ip.y = kv1[1]->GetRefPoint(ky, ijk[1]+orders[1]+1);
+      for (int i = 0; i <= orders[0]; i++, o++)
+      {
+         real_t kx = kv[0]->GetBotella(ijk[0] + i);
+         if (!kv[0]->inSpan(kx, ijk[0]+orders[0])) { continue; }
+         ip.x = kv[0]->GetRefPoint(kx, ijk[0]+orders[0]);
+
+         Trans.SetIntPoint(&ip);
+         vc.Eval(x, Trans, ip);
+
+         Trans.AdjugateJacobian().Mult(x,mx);
+         dofs(o) = mx(1);
+      }
+   }
+}
+
 NURBS_HDiv2DFiniteElement::~NURBS_HDiv2DFiniteElement()
 {
   if (kv1[0]) { delete kv1[0]; }
@@ -696,6 +928,120 @@ void NURBS_HDiv3DFiniteElement::CalcDivShape(const IntegrationPoint &ip,
   }
 }

+
+void NURBS_HDiv3DFiniteElement::Project(VectorCoefficient &vc,
+                                        ElementTransformation &Trans,
+                                        Vector &dofs) const
+{
+   MFEM_ASSERT(dofs.Size() == dof, "");
+   MFEM_ASSERT(vc.GetVDim() == 3, "");
+   Vector x(2), mx(3);
+   IntegrationPoint ip;
+
+   int o = 0;
+
+   for (int k = 0; k <= orders[2]; k++)
+   {
+      real_t kz = kv[2]->GetBotella(ijk[2] + k);
+      if (!kv[2]->inSpan(kz, ijk[2]+orders[2]))
+      {
+         o += (orders[0] + 2)*(orders[1] + 1);
+         continue;
+      }
+      ip.z = kv[2]->GetRefPoint(kz, ijk[2]+orders[2]);
+      for (int j = 0; j <= orders[1]; j++)
+      {
+         real_t ky = kv[1]->GetBotella(ijk[1] + j);
+         if (!kv[1]->inSpan(ky, ijk[1]+orders[1]))
+         {
+            o += orders[0] + 2;
+            continue;
+         }
+         ip.y = kv[1]->GetRefPoint(ky, ijk[1]+orders[1]);
+         for (int i = 0; i <= orders[0]+1; i++, o++)
+         {
+            real_t kx = kv1[0]->GetBotella(ijk[0] + i);
+            if (!kv1[0]->inSpan(kx, ijk[0]+orders[0]+1)) { continue; }
+            ip.x = kv1[0]->GetRefPoint(kx, ijk[0]+orders[0]+1);
+
+            Trans.SetIntPoint(&ip);
+            vc.Eval(x, Trans, ip);
+
+            Trans.AdjugateJacobian().Mult(x,mx);
+            dofs(o) = mx(0);
+         }
+      }
+   }
+
+   for (int k = 0; k <= orders[2]; k++)
+   {
+      real_t kz = kv[2]->GetBotella(ijk[2] + k);
+      if (!kv[2]->inSpan(kz, ijk[2]+orders[2]))
+      {
+         o += (orders[0] + 1)*(orders[1] + 2);
+         continue;
+      }
+      ip.z = kv[2]->GetRefPoint(kz, ijk[2]+orders[2]);
+      for (int j = 0; j <= orders[1]+1; j++)
+      {
+         real_t ky = kv1[1]->GetBotella(ijk[1] + j);
+         if (!kv1[1]->inSpan(ky, ijk[1]+orders[1]+1))
+         {
+            o += orders[0] + 1;
+            continue;
+         }
+         ip.y = kv1[1]->GetRefPoint(ky, ijk[1]+orders[1]+1);
+         for (int i = 0; i <= orders[0]; i++, o++)
+         {
+            real_t kx = kv[0]->GetBotella(ijk[0] + i);
+            if (!kv[0]->inSpan(kx, ijk[0]+orders[0])) { continue; }
+            ip.x = kv[0]->GetRefPoint(kx, ijk[0]+orders[0]);
+
+            Trans.SetIntPoint(&ip);
+            vc.Eval(x, Trans, ip);
+
+            Trans.AdjugateJacobian().Mult(x,mx);
+            dofs(o) = mx(1);
+         }
+      }
+   }
+
+   for (int k = 0; k <= orders[2]+1; k++)
+   {
+      real_t kz = kv1[2]->GetBotella(ijk[2] + k);
+      if (!kv1[2]->inSpan(kz, ijk[2]+orders[2]+1))
+      {
+         o += (orders[0] + 1)*(orders[1] + 1);
+         continue;
+      }
+      ip.z = kv1[2]->GetRefPoint(kz, ijk[2]+orders[2]+1);
+      for (int j = 0; j <= orders[1]; j++)
+      {
+         real_t ky = kv[1]->GetBotella(ijk[1] + j);
+         if (!kv[1]->inSpan(ky, ijk[1]+orders[1]))
+         {
+            o += orders[0] + 1;
+            continue;
+         }
+         ip.y = kv[1]->GetRefPoint(ky, ijk[1]+orders[1]);
+         for (int i = 0; i <= orders[0]; i++, o++)
+         {
+            real_t kx = kv[0]->GetBotella(ijk[0] + i);
+            if (!kv[0]->inSpan(kx, ijk[0]+orders[0])) { continue; }
+            ip.x = kv[0]->GetRefPoint(kx, ijk[0]+orders[0]);
+
+            Trans.SetIntPoint(&ip);
+            vc.Eval(x, Trans, ip);
+
+            Trans.AdjugateJacobian().Mult(x,mx);
+            dofs(o) = mx(2);
+         }
+      }
+   }
+
+}
+
+
 NURBS_HDiv3DFiniteElement::~NURBS_HDiv3DFiniteElement()
 {
   if (kv1[0]) { delete kv1[0]; }
@@ -817,13 +1163,68 @@ void NURBS_HCurl2DFiniteElement::CalcCurlShape(const IntegrationPoint &ip,
   }
 }

+void NURBS_HCurl2DFiniteElement::Project(VectorCoefficient &vc,
+                                         ElementTransformation &Trans,
+                                         Vector &dofs) const
+{
+   MFEM_ASSERT(dofs.Size() == dof, "");
+   MFEM_ASSERT(vc.GetVDim() == 2, "");
+   Vector x(2), xm(2);
+   IntegrationPoint ip;
+   int i, j, o;
+   for (o = 0, j = 0; j <= orders[1]+1; j++)
+   {
+      real_t ky = kv1[1]->GetBotella(ijk[1] + j);
+      if (!kv1[1]->inSpan(ky, ijk[1]+orders[1]+1))
+      {
+         o += orders[0] + 1;
+         continue;
+      }
+      ip.y = kv1[1]->GetRefPoint(ky, ijk[1]+orders[1]+1);
+      for (i = 0; i <= orders[0]; i++, o++)
+      {
+         real_t kx = kv[0]->GetBotella(ijk[0] + i);
+         if (!kv[0]->inSpan(kx, ijk[0]+orders[0])) { continue; }
+         ip.x = kv[0]->GetRefPoint(kx, ijk[0]+orders[0]);
+
+         Trans.SetIntPoint(&ip);
+         vc.Eval(x, Trans, ip);
+
+         Trans.Jacobian().MultTranspose(x,xm);
+         dofs(o) = xm(0);
+      }
+   }
+
+   for (j = 0; j <= orders[1]; j++)
+   {
+      real_t ky = kv[1]->GetBotella(ijk[1] + j);
+      if (!kv[1]->inSpan(ky, ijk[1]+orders[1]))
+      {
+         o += orders[0] + 2;
+         continue;
+      }
+      ip.y = kv[1]->GetRefPoint(ky, ijk[1]+orders[1]);
+      for (i = 0; i <= orders[0]+1; i++, o++)
+      {
+         real_t kx = kv1[0]->GetBotella(ijk[0] + i);
+         if (!kv1[0]->inSpan(kx, ijk[0]+orders[0]+1)) { continue; }
+         ip.x = kv1[0]->GetRefPoint(kx, ijk[0]+orders[0]+1);
+
+         Trans.SetIntPoint(&ip);
+         vc.Eval(x, Trans, ip);
+
+         Trans.Jacobian().MultTranspose(x,xm);
+         dofs(o) = xm(1);
+      }
+   }
+}
+
 NURBS_HCurl2DFiniteElement::~NURBS_HCurl2DFiniteElement()
 {
   if (kv1[0]) { delete kv1[0]; }
   if (kv1[1]) { delete kv1[1]; }
 }

-
 void NURBS_HCurl3DFiniteElement::SetOrder() const
 {
   orders[0] = kv[0]->GetOrder();
@@ -1003,11 +1404,124 @@ void NURBS_HCurl3DFiniteElement::CalcCurlShape(const IntegrationPoint &ip,
            curl_shape(o,0) = shape1_x(i)*dsy1_sz;
            curl_shape(o,1) = -dshape1_x(i)*sy1_sz;
            curl_shape(o,2) = 0.0;
+
         }
      }
   }
 }

+void NURBS_HCurl3DFiniteElement::Project(VectorCoefficient &vc,
+                                         ElementTransformation &Trans,
+                                         Vector &dofs) const
+{
+   MFEM_ASSERT(dofs.Size() == dof, "");
+   MFEM_ASSERT(vc.GetVDim() == 3, "");
+   Vector x(3), xm(3);
+   IntegrationPoint ip;
+
+   int o = 0;
+   for (int k = 0; k <= orders[2]+1; k++)
+   {
+      real_t kz = kv1[2]->GetBotella(ijk[2] + k);
+      if (!kv1[2]->inSpan(kz, ijk[2]+orders[2]+1))
+      {
+         o += (orders[0] + 1)*(orders[1] + 2);
+         continue;
+      }
+      ip.z = kv1[2]->GetRefPoint(kz, ijk[2]+orders[2]+1);
+      for (int j = 0; j <= orders[1]+1; j++)
+      {
+         real_t ky = kv1[1]->GetBotella(ijk[1] + j);
+         if (!kv1[1]->inSpan(ky, ijk[1]+orders[1]+1))
+         {
+            o += orders[0] + 1;
+            continue;
+         }
+         ip.y = kv1[1]->GetRefPoint(ky, ijk[1]+orders[1]+1);
+         for (int i = 0; i <= orders[0]; i++, o++)
+         {
+            real_t kx = kv[0]->GetBotella(ijk[0] + i);
+            if (!kv[0]->inSpan(kx, ijk[0]+orders[0])) { continue; }
+            ip.x = kv[0]->GetRefPoint(kx, ijk[0]+orders[0]);
+
+            Trans.SetIntPoint(&ip);
+            vc.Eval(x, Trans, ip);
+
+            Trans.Jacobian().MultTranspose(x,xm);
+            dofs(o) = xm(0);
+         }
+      }
+   }
+
+   for (int k = 0; k <= orders[2]+1; k++)
+   {
+      real_t kz = kv1[2]->GetBotella(ijk[2] + k);
+      if (!kv1[2]->inSpan(kz, ijk[2]+orders[2]+1))
+      {
+         o += (orders[0] + 2)*(orders[1] + 1);
+         continue;
+      }
+      ip.z = kv1[2]->GetRefPoint(kz, ijk[2]+orders[2]+1);
+      for (int j = 0; j <= orders[1]; j++)
+      {
+         real_t ky = kv[1]->GetBotella(ijk[1] + j);
+         if (!kv[1]->inSpan(ky, ijk[1]+orders[1]))
+         {
+            o += orders[0] + 2;
+            continue;
+         }
+         ip.y = kv[1]->GetRefPoint(ky, ijk[1]+orders[1]);
+         for (int i = 0; i <= orders[0]+1; i++, o++)
+         {
+            real_t kx = kv1[0]->GetBotella(ijk[0] + i);
+            if (!kv1[0]->inSpan(kx, ijk[0]+orders[0]+1)) { continue; }
+            ip.x = kv1[0]->GetRefPoint(kx, ijk[0]+orders[0]+1);
+
+            Trans.SetIntPoint(&ip);
+            vc.Eval(x, Trans, ip);
+
+            Trans.Jacobian().MultTranspose(x,xm);
+            dofs(o) = xm(1);
+         }
+      }
+   }
+
+   for (int k = 0; k <= orders[2]; k++)
+   {
+      real_t kz = kv[2]->GetBotella(ijk[2] + k);
+      if (!kv[2]->inSpan(kz, ijk[2]+orders[2]))
+      {
+         o += (orders[0] + 2)*(orders[1] + 2);
+         continue;
+      }
+      ip.z = kv[2]->GetRefPoint(kz, ijk[2]+orders[2]);
+      for (int j = 0; j <= orders[1]+1; j++)
+      {
+         real_t ky = kv1[1]->GetBotella(ijk[1] + j);
+         if (!kv1[1]->inSpan(ky, ijk[1]+orders[1]+1))
+         {
+            o += orders[0] + 2;
+            continue;
+         }
+         ip.y = kv1[1]->GetRefPoint(ky, ijk[1]+orders[1]+1);
+         for (int i = 0; i <= orders[0]+1; i++, o++)
+         {
+            real_t kx = kv1[0]->GetBotella(ijk[0] + i);
+            if (!kv1[0]->inSpan(kx, ijk[0]+orders[0]+1)) { continue; }
+            ip.x = kv1[0]->GetRefPoint(kx, ijk[0]+orders[0]+1);
+
+            Trans.SetIntPoint(&ip);
+            vc.Eval(x, Trans, ip);
+
+            Trans.Jacobian().MultTranspose(x,xm);
+            dofs(o) = xm(2);
+         }
+      }
+   }
+
+}
+
+
 NURBS_HCurl3DFiniteElement::~NURBS_HCurl3DFiniteElement()
 {
   if (kv1[0]) { delete kv1[0]; }
@@ -86,6 +86,18 @@ public:
                   DenseMatrix &dshape) const override;
   void CalcHessian (const IntegrationPoint &ip,
                     DenseMatrix &hessian) const override;
+
+   using FiniteElement::Project;
+
+   /** Evaluate the dofs that are defined on this element.
+       Dofs that can not be evaluated will remain unmodified. */
+   void Project(Coefficient &coeff,
+                ElementTransformation &Trans, Vector &dofs) const override;
+
+   /** Evaluate the dofs that are defined on this element.
+       Dofs that can not be evaluated will remain unmodified. */
+   void Project(VectorCoefficient &vcoeff,
+                ElementTransformation &Trans, Vector &dofs) const override;
 };

 /// An arbitrary order 2D NURBS element on a square
@@ -121,6 +133,18 @@ public:
                   DenseMatrix &dshape) const override;
   void CalcHessian (const IntegrationPoint &ip,
                     DenseMatrix &hessian) const override;
+
+   using FiniteElement::Project;
+
+   /** Evaluate the dofs that are defined on this element.
+       Dofs that can not be evaluated will remain unmodified. */
+   void Project(Coefficient &coeff,
+                ElementTransformation &Trans, Vector &dofs) const override;
+
+   /** Evaluate the dofs that are defined on this element.
+       Dofs that can not be evaluated will remain unmodified. */
+   void Project(VectorCoefficient &vcoeff,
+                ElementTransformation &Trans, Vector &dofs) const override;
 };

 /// An arbitrary order 3D NURBS element on a cube
@@ -161,6 +185,18 @@ public:
                   DenseMatrix &dshape) const override;
   void CalcHessian (const IntegrationPoint &ip,
                     DenseMatrix &hessian) const override;
+
+   using FiniteElement::Project;
+
+   /** Evaluate the dofs that are defined on this element.
+       Dofs that can not be evaluated will remain unmodified. */
+   void Project(Coefficient &coeff,
+                ElementTransformation &Trans, Vector &dofs) const override;
+
+   /** Evaluate the dofs that are defined on this element.
+       Dofs that can not be evaluated will remain unmodified. */
+   void Project(VectorCoefficient &vcoeff,
+                ElementTransformation &Trans, Vector &dofs) const override;
 };


@@ -242,6 +278,13 @@ public:
   void CalcDivShape(const IntegrationPoint &ip,
                     Vector &divshape) const override;

+   using FiniteElement::Project;
+
+   /** Evaluate the dofs that are defined on this element.
+       Dofs that can not be evaluated will remain unmodified. */
+   void Project(VectorCoefficient &vcoeff,
+                ElementTransformation &Trans, Vector &dofs) const override;
+
   ~NURBS_HDiv2DFiniteElement();
 };

@@ -336,6 +379,13 @@ public:
   void CalcDivShape(const IntegrationPoint &ip,
                     Vector &divshape) const override;

+   using FiniteElement::Project;
+
+   /** Evaluate the dofs that are defined on this element.
+       Dofs that can not be evaluated will remain unmodified. */
+   void Project(VectorCoefficient &vcoeff,
+                ElementTransformation &Trans, Vector &dofs) const override;
+
   ~NURBS_HDiv3DFiniteElement();
 };

@@ -415,6 +465,13 @@ public:
   void CalcCurlShape(const IntegrationPoint &ip,
                      DenseMatrix &curl_shape) const override;

+   using FiniteElement::Project;
+
+   /** Evaluate the dofs that are defined on this element.
+       Dofs that can not be evaluated will remain unmodified. */
+   void Project(VectorCoefficient &vcoeff,
+                ElementTransformation &Trans, Vector &dofs) const override;
+
   ~NURBS_HCurl2DFiniteElement();
 };

@@ -506,6 +563,13 @@ public:
   void CalcCurlShape(const IntegrationPoint &ip,
                      DenseMatrix &curl_shape) const override;

+   using FiniteElement::Project;
+
+   /** Evaluate the dofs that are defined on this element.
+       Dofs that can not be evaluated will remain unmodified. */
+   void Project(VectorCoefficient &vcoeff,
+                ElementTransformation &Trans, Vector &dofs) const override;
+
   ~NURBS_HCurl3DFiniteElement();
 };

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

+   int GetPhysRangeDim(int space_dim) const override { return 2; }
+   int GetPhysCurlDim(int space_dim) const override { return 0; }
+
   void CalcVShape(const IntegrationPoint &ip,
                   DenseMatrix &shape) const override;

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

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

   void CalcVShape(ElementTransformation &Trans,
@@ -111,36 +111,36 @@ public:
   | :------: | :---: | :---: | :-------: | :-----: | :---: |
   | H1_[DIM]_[ORDER] | H1 | * | 1 | VALUE | H1 nodal elements |
   | H1@[BTYPE]_[DIM]_[ORDER] | H1 | * | * | VALUE | H1 nodal elements |
-   | H1Pos_[DIM]_[ORDER] | H1 | * | 1 | VALUE | H1 nodal elements |
+   | H1Pos_[DIM]_[ORDER] | H1 | * | 2 | VALUE | H1 nodal elements |
   | H1Pos_Trace_[DIM]_[ORDER] | H^{1/2} | * | 2 | VALUE | H^{1/2}-conforming trace elements for H1 defined on the interface between mesh elements (faces,edges,vertices) |
   | H1_Trace_[DIM]_[ORDER] | H^{1/2} | * | 1 | VALUE | H^{1/2}-conforming trace elements for H1 defined on the interface between mesh elements (faces,edges,vertices) |
-   | H1_Trace@[BTYPE]_[DIM]_[ORDER] | H^{1/2} | * | 1 | VALUE | H^{1/2}-conforming trace elements for H1 defined on the interface between mesh elements (faces,edges,vertices) |
+   | H1_Trace@[BTYPE]_[DIM]_[ORDER] | H^{1/2} | * | * | VALUE | H^{1/2}-conforming trace elements for H1 defined on the interface between mesh elements (faces,edges,vertices) |
   | ND_[DIM]_[ORDER] | H(curl) | * | 1 / 0 | H_CURL | Nedelec vector elements |
   | ND@[CBTYPE][OBTYPE]_[DIM]_[ORDER] | H(curl) | * | * / * | H_CURL | Nedelec vector elements |
-   | ND_Trace_[DIM]_[ORDER] | H^{1/2} | * | 1 / 0  | H_CURL | H^{1/2}-conforming trace elements for H(curl) defined on the interface between mesh elements (faces) |
-   | ND_Trace@[CBTYPE][OBTYPE]_[DIM]_[ORDER] | H^{1/2} | * | 1 / 0 | H_CURL | H^{1/2}-conforming trace elements for H(curl) defined on the interface between mesh elements (faces) |
+   | ND_Trace_[DIM]_[ORDER] | H^{1/2} | * | 1 / 0  | H_CURL | H^{1/2}-conforming trace elements for H(curl) defined on the interface between mesh elements (faces,edges) |
+   | ND_Trace@[CBTYPE][OBTYPE]_[DIM]_[ORDER] | H^{1/2} | * | * / * | H_CURL | H^{1/2}-conforming trace elements for H(curl) defined on the interface between mesh elements (faces,edges) |
   | ND_R1D_[DIM]_[ORDER] | H(curl) | * | 1 / 0 | H_CURL | 3D H(curl)-conforming Nedelec vector elements in 1D. |
   | ND_R1D@[CBTYPE][OBTYPE]_[DIM]_[ORDER] | H(curl) | * | * / * | H_CURL | 3D H(curl)-conforming Nedelec vector elements in 1D. |
   | ND_R2D_[DIM]_[ORDER] | H(curl) | * | 1 / 0 | H_CURL | 3D H(curl)-conforming Nedelec vector elements in 2D. |
   | ND_R2D@[CBTYPE][OBTYPE]_[DIM]_[ORDER] | H(curl) | * | * / * | H_CURL | 3D H(curl)-conforming Nedelec vector elements in 2D. |
   | RT_[DIM]_[ORDER] | H(div) | * | 1 / 0 | H_DIV | Raviart-Thomas vector elements |
   | RT@[CBTYPE][OBTYPE]_[DIM]_[ORDER] | H(div) | * | * / * | H_DIV | Raviart-Thomas vector elements |
-   | RT_Trace_[DIM]_[ORDER] | H^{1/2} | * | 1 / 0 | INTEGRAL | H^{1/2}-conforming trace elements for H(div) defined on the interface between mesh elements (faces) |
-   | RT_ValTrace_[DIM]_[ORDER] | H^{1/2} | * | 1 / 0 | VALUE | H^{1/2}-conforming trace elements for H(div) defined on the interface between mesh elements (faces) |
-   | RT_Trace@[BTYPE]_[DIM]_[ORDER] | H^{1/2} | * | 1 / 0 | INTEGRAL | H^{1/2}-conforming trace elements for H(div) defined on the interface between mesh elements (faces) |
-   | RT_ValTrace@[BTYPE]_[DIM]_[ORDER] |  H^{1/2} | * | 1 / 0 | VALUE | H^{1/2}-conforming trace elements for H(div) defined on the interface between mesh elements (faces) |
+   | RT_Trace_[DIM]_[ORDER] | H^{1/2} | * | 0 | INTEGRAL | H^{1/2}-conforming trace elements for H(div) defined on the interface between mesh elements (faces) |
+   | RT_ValTrace_[DIM]_[ORDER] | H^{1/2} | * | 0 | VALUE | H^{1/2}-conforming trace elements for H(div) defined on the interface between mesh elements (faces) |
+   | RT_Trace@[BTYPE]_[DIM]_[ORDER] | H^{1/2} | * | * | INTEGRAL | H^{1/2}-conforming trace elements for H(div) defined on the interface between mesh elements (faces) |
+   | RT_ValTrace@[BTYPE]_[DIM]_[ORDER] |  H^{1/2} | * | * | VALUE | H^{1/2}-conforming trace elements for H(div) defined on the interface between mesh elements (faces) |
   | RT_R1D_[DIM]_[ORDER] | H(div) | * | 1 / 0 | H_DIV | 3D H(div)-conforming Raviart-Thomas vector elements in 1D. |
   | RT_R1D@[CBTYPE][OBTYPE]_[DIM]_[ORDER] | H(div) | * | * / * | H_DIV | 3D H(div)-conforming Raviart-Thomas vector elements in 1D. |
   | RT_R2D_[DIM]_[ORDER] | H(div) | * | 1 / 0 | H_DIV | 3D H(div)-conforming Raviart-Thomas vector elements in 2D. |
   | RT_R2D@[CBTYPE][OBTYPE]_[DIM]_[ORDER] | H(div) | * | * / * | H_DIV | 3D H(div)-conforming Raviart-Thomas vector elements in 2D. |
   | L2_[DIM]_[ORDER] | L2 | * | 0 | VALUE | Discontinuous L2 elements |
-   | L2_T[BTYPE]_[DIM]_[ORDER] | L2 | * | 0 | VALUE | Discontinuous L2 elements |
+   | L2_T[BTYPE]_[DIM]_[ORDER] | L2 | * | * | VALUE | Discontinuous L2 elements |
   | L2Int_[DIM]_[ORDER] | L2 | * | 0 | INTEGRAL | Discontinuous L2 elements |
-   | L2Int_T[BTYPE]_[DIM]_[ORDER] | L2 | * | 0 | INTEGRAL | Discontinuous L2 elements |
+   | L2Int_T[BTYPE]_[DIM]_[ORDER] | L2 | * | * | INTEGRAL | Discontinuous L2 elements |
   | DG_Iface_[DIM]_[ORDER] | - | * | 0 | VALUE | Discontinuous elements on the interface between mesh elements (faces) |
-   | DG_Iface@[BTYPE]_[DIM]_[ORDER] | - | * | 0 | VALUE | Discontinuous elements on the interface between mesh elements (faces) |
+   | DG_Iface@[BTYPE]_[DIM]_[ORDER] | - | * | * | VALUE | Discontinuous elements on the interface between mesh elements (faces) |
   | DG_IntIface_[DIM]_[ORDER] | - | * | 0 | INTEGRAL | Discontinuous elements on the interface between mesh elements (faces) |
-   | DG_IntIface@[BTYPE]_[DIM]_[ORDER] | - | * | 0 | INTEGRAL | Discontinuous elements on the interface between mesh elements (faces) |
+   | DG_IntIface@[BTYPE]_[DIM]_[ORDER] | - | * | * | INTEGRAL | Discontinuous elements on the interface between mesh elements (faces) |
   | NURBS[ORDER] | - | * | - | VALUE | Non-Uniform Rational B-Splines (NURBS) elements |
   | LinearNonConf3D | - | 1 | 1 | VALUE | Piecewise-linear nonconforming finite elements in 3D |
   | CrouzeixRaviart | - | - | - | - | Crouzeix-Raviart nonconforming elements in 2D |
@@ -172,7 +172,7 @@ public:
   | :------: | :--------: |
   | [DIM]    | Dimension of the elements (1D, 2D, 3D) |
   | [ORDER]  | Approximation order of the elements (P0, P1, P2, ...) |
-   | [BTYPE]  | BasisType of the element (0-GaussLegendre, 1 - GaussLobatto, 2-Bernstein, 3-OpenUniform, 4-CloseUniform, 5-OpenHalfUniform) |
+   | [BTYPE]  | BasisType of the element (0-GaussLegendre, 1-GaussLobatto, 2-Bernstein, 3-OpenUniform, 4-CloseUniform, 5-OpenHalfUniform 6-Serendipity 7-ClosedGL 8-IntegratedGLL) |
   | [OBTYPE] | Open BasisType of the element for elements which have both types |
   | [CBTYPE] | Closed BasisType of the element for elements which have both types |

@@ -62,10 +62,6 @@
 #include "pnonlinearform.hpp"
 #endif

-#ifdef MFEM_USE_SIDRE
-#include "sidredatacollection.hpp"
-#endif
-
 #ifdef MFEM_USE_CONDUIT
 #include "conduitdatacollection.hpp"
 #endif
@@ -282,14 +282,7 @@ int FiniteElementSpace::DofToVDof(int dof, int vd, int ndofs_) const
 void FiniteElementSpace::AdjustVDofs(Array<int> &vdofs)
 {
   int n = vdofs.Size(), *vdof = vdofs;
-   for (int i = 0; i < n; i++)
-   {
-      int j;
-      if ((j = vdof[i]) < 0)
-      {
-         vdof[i] = -1-j;
-      }
-   }
+   for (int i = 0; i < n; i++) { vdof[i] = UnsignIndex(vdof[i]); }
 }

 void FiniteElementSpace::GetElementVDofs(int i, Array<int> &vdofs,
@@ -483,13 +476,14 @@ void FiniteElementSpace::ReorderElementToDofTable()
   for (int k = 0, dof_counter = 0; k < nnz; k++)
   {
      const int sdof = J[k]; // signed dof
-      const int dof = (sdof < 0) ? -1-sdof : sdof;
+      const int dof = UnsignIndex(sdof);
      int new_dof = dof_marker[dof];
      if (new_dof < 0)
      {
         dof_marker[dof] = new_dof = dof_counter++;
      }
-      J[k] = (sdof < 0) ? -1-new_dof : new_dof; // preserve the sign of sdof
+      // Preserve the sign of sdof
+      J[k] = (sdof < 0) ? FlipIndexSign(new_dof) : new_dof;
   }
 }

@@ -547,7 +541,7 @@ void MarkDofs(const Array<int> &dofs, Array<int> &mark_array)
 {
   for (auto d : dofs)
   {
-      mark_array[d >= 0 ? d : -1 - d] = -1;
+      mark_array[UnsignIndex(d)] = -1;
   }
 }

@@ -931,7 +925,7 @@ void FiniteElementSpace::AddDependencies(
            if (std::abs(coef) > 1e-12)
            {
               const int mdof = master_dofs[j];
-               if (mdof != sdof && mdof != (-1-sdof))
+               if (mdof != sdof && mdof != FlipIndexSign(sdof))
               {
                  deps.Add(sdof, mdof, coef);
               }
@@ -1024,7 +1018,7 @@ int FiniteElementSpace::GetDegenerateFaceDofs(int index, Array<int> &dofs,
   // FiniteElementSpace::AddDependencies.

   Array<int> edof;
-   int order = GetEdgeDofs(-1 - index, edof, variant);
+   int order = GetEdgeDofs(FlipIndexSign(index), edof, variant);

   int nv = fec->DofForGeometry(Geometry::POINT);
   int ne = fec->DofForGeometry(Geometry::SEGMENT);
@@ -1516,36 +1510,76 @@ const FaceRestriction *FiniteElementSpace::GetFaceRestriction(
   const bool is_dg_space = IsDGSpace();
   const L2FaceValues m = (is_dg_space && mul==L2FaceValues::DoubleValued) ?
                          L2FaceValues::DoubleValued : L2FaceValues::SingleValued;
-   key_face key = std::make_tuple(is_dg_space, f_ordering, type, m);
+   auto key = std::make_tuple(is_dg_space, f_ordering, type, m);
   auto itr = L2F.find(key);
   if (itr != L2F.end())
   {
-      return itr->second;
+      return itr->second.get();
   }
   else
   {
-      FaceRestriction *res;
+      std::unique_ptr<FaceRestriction> res;
      if (is_dg_space)
      {
         if (Conforming())
         {
-            res = new L2FaceRestriction(*this, f_ordering, type, m);
+            res.reset(new L2FaceRestriction(*this, f_ordering, type, m));
         }
         else
         {
-            res = new NCL2FaceRestriction(*this, f_ordering, type, m);
+            res.reset(new NCL2FaceRestriction(*this, f_ordering, type, m));
         }
      }
      else if (dynamic_cast<const DG_Interface_FECollection*>(fec))
      {
-         res = new L2InterfaceFaceRestriction(*this, f_ordering, type);
+         res.reset(new L2InterfaceFaceRestriction(*this, f_ordering, type));
      }
      else
      {
-         res = new ConformingFaceRestriction(*this, f_ordering, type);
+         res.reset(new ConformingFaceRestriction(*this, f_ordering, type));
      }
-      L2F[key] = res;
-      return res;
+      return L2F.emplace(key, std::move(res)).first->second.get();
+   }
+}
+
+const InterpolationManager &FiniteElementSpace::GetInterpolationManager(
+   ElementDofOrdering f_ordering, FaceType type) const
+{
+   const auto key = make_tuple(f_ordering, type);
+
+   auto it = interpolations.find(key);
+   if (it != interpolations.end())
+   {
+      return *it->second;
+   }
+   else
+   {
+      auto interp = make_unique<InterpolationManager>(*this, f_ordering, type);
+
+      int face_idx = 0;
+      for (int f = 0; f < mesh->GetNumFacesWithGhost(); ++f)
+      {
+         Mesh::FaceInformation face = mesh->GetFaceInformation(f);
+         if (!face.IsOfFaceType(type) || face.IsNonconformingCoarse())
+         {
+            continue;
+         }
+         if (face.IsConforming() || face.IsBoundary())
+         {
+            interp->RegisterFaceConformingInterpolation(face, face_idx);
+         }
+         else
+         {
+            interp->RegisterFaceCoarseToFineInterpolation(face, face_idx);
+         }
+         ++face_idx;
+      }
+
+      // Transform the interpolation matrix map into contiguous memory.
+      interp->LinearizeInterpolatorMapIntoVector();
+      interp->InitializeNCInterpConfig();
+
+      return *interpolations.emplace(key, std::move(interp)).first->second;
   }
 }

@@ -1670,8 +1704,8 @@ SparseMatrix *FiniteElementSpace::RefinementMatrix_main(

         for (int i = 0; i < fine_ldof; i++)
         {
-            int r = DofToVDof(dofs[i], vd);
-            int m = (r >= 0) ? r : (-1 - r);
+            const int r = DofToVDof(dofs[i], vd);
+            const int m = UnsignIndex(r);

            if (!mark[m])
            {
@@ -1732,7 +1766,7 @@ SparseMatrix *FiniteElementSpace::VariableOrderRefinementMatrix(
         for (int i = 0; i < fine_ldof; i++)
         {
            const int r = DofToVDof(dofs[i], vd);
-            int m = (r >= 0) ? r : (-1 - r);
+            const int m = UnsignIndex(r);

            if (!mark[m])
            {
@@ -2442,8 +2476,8 @@ SparseMatrix* FiniteElementSpace::DerefinementMatrix(int old_ndofs,
         {
            if (!std::isfinite(lR(i, 0))) { continue; }

-            int r = DofToVDof(dofs[i], vd);
-            int m = (r >= 0) ? r : (-1 - r);
+            const int r = DofToVDof(dofs[i], vd);
+            const int m = UnsignIndex(r);

            if (is_dg || !mark[m])
            {
@@ -3161,7 +3195,7 @@ void FiniteElementSpace::CalcEdgeFaceVarOrders(
            else
            {
               // degenerate face (i.e., edge-face constraint)
-               slave_orders |= edge_orders[-1 - slave.index];
+               slave_orders |= edge_orders[FlipIndexSign(slave.index)];
            }
         }

@@ -3900,6 +3934,16 @@ const FiniteElement *FiniteElementSpace::GetBE(int i) const
   return BE;
 }

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

+const FiniteElement *FiniteElementSpace::GetTypicalFaceElement() const
+{
+   return fec->FiniteElementForGeometry(mesh->GetTypicalFaceGeometry());
+}
+
 const FiniteElement *FiniteElementSpace::GetEdgeElement(int i,
                                                        int variant) const
 {
@@ -3969,11 +4018,8 @@ void FiniteElementSpace::Destroy()
      delete E2Q_array[i];
   }
   E2Q_array.SetSize(0);
-   for (auto &x : L2F)
-   {
-      delete x.second;
-   }
   L2F.clear();
+   interpolations.clear();
   for (int i = 0; i < E2IFQ_array.Size(); i++)
   {
      delete E2IFQ_array[i];
@@ -13,6 +13,7 @@
 #define MFEM_FESPACE

 #include "../config/config.hpp"
+#include "../general/hash_util.hpp"
 #include "../linalg/ordering.hpp"
 #include "../linalg/sparsemat.hpp"
 #include "../mesh/mesh.hpp"
@@ -320,18 +321,11 @@ protected:
   mutable OperatorHandle L2E_nat, L2E_lex;
   /// The face restriction operators, see GetFaceRestriction().
   using key_face = std::tuple<bool, ElementDofOrdering, FaceType, L2FaceValues>;
-   struct key_hash
-   {
-      std::size_t operator()(const key_face& k) const
-      {
-         return std::get<0>(k)
-                + 2 * (int)std::get<1>(k)
-                + 4 * (int)std::get<2>(k)
-                + 8 * (int)std::get<3>(k);
-      }
-   };
-   using map_L2F = std::unordered_map<const key_face,FaceRestriction*,key_hash>;
-   mutable map_L2F L2F;
+   mutable std::unordered_map<key_face,std::unique_ptr<FaceRestriction>,
+           TupleHasher> L2F;
+
+   mutable std::unordered_map<std::tuple<ElementDofOrdering,FaceType>,
+           std::unique_ptr<InterpolationManager>, TupleHasher> interpolations;

   mutable Array<QuadratureInterpolator*> E2Q_array;
   mutable Array<FaceQuadratureInterpolator*> E2IFQ_array;
@@ -751,6 +745,9 @@ public:
      ElementDofOrdering f_ordering, FaceType,
      L2FaceValues mul = L2FaceValues::DoubleValued) const;

+   const InterpolationManager &GetInterpolationManager(
+      ElementDofOrdering f_ordering, FaceType type) const;
+
   /** @brief Return a QuadratureInterpolator that interpolates E-vectors to
       quadrature point values and/or derivatives (Q-vectors). */
   /** An E-vector represents the element-wise discontinuous version of the FE
@@ -842,7 +839,7 @@ public:
       Note: For vector-valued elements, the results pads up the range dimension
       to the spatial dimension. E.g., consider a stack of 5 vector-valued
       elements each representing 2D vectors, living in a 3 dimensional space.
-       Then this fucntion would give 15, not 10.
+       Then this function would give 15, not 10.
       */
   int GetVectorDim() const;

@@ -1153,7 +1150,7 @@ public:

   /// Helper to return the DOF associated with a sign encoded DOF
   static inline int DecodeDof(int dof)
-   { return (dof >= 0) ? dof : (-1 - dof); }
+   { return UnsignIndex(dof); }

   /// Helper to determine the DOF and sign of a sign encoded DOF
   static inline int DecodeDof(int dof, real_t& sign)
@@ -1326,12 +1323,24 @@ public:
        associated with i'th boundary face in the mesh object. */
   const FiniteElement *GetBE(int i) const;

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

+   /// @brief Return a "typical" face element.
+   ///
+   /// This can be used in situations where the local mesh partition may be
+   /// empty.
+   const FiniteElement *GetTypicalFaceElement() const;
+
   /** @brief Returns pointer to the FiniteElement in the FiniteElementCollection
        associated with i'th edge in the mesh object. */
   const FiniteElement *GetEdgeElement(int i, int variant = 0) const;
@@ -27,6 +27,24 @@
 namespace mfem
 {

+/** This enumerated type describes the three main projection types:
+    - ELEMENT, assigns the degree of freedom per element, as specified in the
+      specific element
+    - GLOBAL_L2, solves a global L2 projection
+    - ELEMENT_L2, solves a element level L2 projection. Inter element
+      connectivity is dealt with similar as in:
+      Bezier-Projection : A unified approach for local projection and
+      quadrature-free refinement and coarsening of NURBS and T-splines with
+      particular application to isogeometric design and analysis
+      [CMAME (284) 2015 pg 55-105]
+    - DEFAULT, for NURBS spaces this is ELEMENT_L2, while for all other spaces
+      this ELEMENT.
+   Note 1: ELEMENT_L2 also works for non NURBS elements
+   Note 2: For NURBS elements the ELEMENT projection gives results without
+   over and undershoots. However, the gradient near the boundary does not
+   converge.*/
+enum class ProjectType { DEFAULT, ELEMENT, GLOBAL_L2, ELEMENT_L2 };
+
 /// Class for grid function - Vector with associated FE space.
 class GridFunction : public Vector
 {
@@ -66,13 +84,17 @@ protected:
       degree of freedom. */
   void ProjectDiscCoefficient(VectorCoefficient &coeff, Array<int> &dof_attr);

+   /** Helper function for ProjectCoefficientElementL2 */
+   void ProjectCoefficientElementL2_(Coefficient &coeff, Vector &sol, Vector &Va);
+   void ProjectCoefficientElementL2_(VectorCoefficient &vcoeff, Vector &sol,
+                                     Vector &Va);
+
   /// Loading helper.
   void LegacyNCReorder();

   void Destroy();

 public:
-
   GridFunction() { fes = NULL; fec_owned = NULL; fes_sequence = 0; UseDevice(true); }

   /// Copy constructor. The internal true-dof vector #t_vec is not copied.
@@ -84,6 +106,10 @@ public:
   GridFunction(FiniteElementSpace *f) : Vector(f->GetVSize())
   { fes = f; fec_owned = NULL; fes_sequence = f->GetSequence(); UseDevice(true); }

+   /// Same as above but specify the memory type
+   GridFunction(FiniteElementSpace *f, MemoryType mt) : Vector(f->GetVSize(), mt)
+   { fes = f; fec_owned = NULL; fes_sequence = f->GetSequence(); UseDevice(true); }
+
   /// Construct a GridFunction using previously allocated array @a data.
   /** The GridFunction does not assume ownership of @a data which is assumed to
       be of size at least `f->GetVSize()`. Similar to the Vector constructor
@@ -124,11 +150,13 @@ public:

   FiniteElementCollection *OwnFEC() { return fec_owned; }

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

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

   /// Read only access to the (optional) internal true-dof Vector.
   const Vector &GetTrueVector() const
@@ -420,9 +448,30 @@ public:
   /** @brief Project @a coeff Coefficient to @a this GridFunction. The
       projection computation depends on the choice of the FiniteElementSpace
       #fes. Note that this is usually interpolation at the degrees of freedom
-       in each element (not L2 projection). For NURBS spaces these degrees of
-       freedom are not available and L2 projection is resorted to as fallback. */
-   virtual void ProjectCoefficient(Coefficient &coeff);
+       in each element (not L2 projection). For elements without a projection
+       member function one could use ProjectCoefficientGlobalL2 instead.
+       NOTE: For parallel simulations with NURBS elements some dofs might
+       not be defined, if the evaluation point does not reside on this rank.
+       If that is the case it is defined on another rank, and the issue is
+       rectified with the appropriate communication, see in ParGridFunction.
+       */
+   virtual void ProjectCoefficient(Coefficient &coeff,
+                                   ProjectType type = ProjectType::DEFAULT);
+
+   /** @brief Project @a coeff Coefficient to @a this GridFunction. The
+       projection is a global L2 projection. This routine can be used a
+       fallback for elements without a projection member function.*/
+   virtual void ProjectCoefficientGlobalL2(Coefficient &coeff,
+                                           real_t rtol = 1e-12,
+                                           int iter = 1000);
+
+   /** @brief Project @a coeff Coefficient to @a this GridFunction. The
+       projection is an element local L2 projection, with an appropriate
+       weighting for Dofs that are shared between elements. Inspired on
+       Bezier-Projection [CMAME (284) 2015 pg 55-105]
+       This routine can be used a fallback for elements without a projection
+       member function.*/
+   virtual void ProjectCoefficientElementL2(Coefficient &coeff);

   /** @brief Project @a coeff Coefficient to @a this GridFunction, using one
       element for each degree of freedom in @a dofs and nodal interpolation on
@@ -432,9 +481,26 @@ public:
   /** @brief Project @a vcoeff VectorCoefficient to @a this GridFunction. The
       projection computation depends on the choice of the FiniteElementSpace
       #fes. Note that this is usually interpolation at the degrees of freedom
-       in each element (not L2 projection). For NURBS spaces these degrees of
-       freedom are not available and L2 projection is resorted to as fallback. */
-   void ProjectCoefficient(VectorCoefficient &vcoeff);
+       in each element (not L2 projection). For elements without a projection
+       member function one could use ProjectCoefficientGlobalL2 instead.
+       NOTE: For parallel simulations with NURBS elements some dofs might
+       not be defined, if the evaluation point does not reside on this rank.
+       If that is the case it is defined on another rank, and the issue is
+       rectified with the appropriate communication, see in ParGridFunction.*/
+   virtual void ProjectCoefficient(VectorCoefficient &vcoeff,
+                                   ProjectType type = ProjectType::DEFAULT);
+
+   /** @brief Project @a coeff Coefficient to @a this GridFunction. The
+       projection is a global L2 projection. This routine can be used a
+       fallback for elements without a projection member function.*/
+   virtual void ProjectCoefficientGlobalL2(VectorCoefficient &vcoeff,
+                                           real_t rtol = 1e-12,
+                                           int iter = 1000);
+
+   /** @brief Project @a coeff Coefficient to @a this GridFunction. The
+       projection is a global L2 projection. This routine can be used a
+       fallback for elements without a projection member function.*/
+   virtual void ProjectCoefficientElementL2(VectorCoefficient &vcoeff);

   /** @brief Project @a vcoeff VectorCoefficient to @a this GridFunction, using
       one element for each degree of freedom in @a dofs and nodal interpolation
@@ -468,6 +534,9 @@ public:
   std::unique_ptr<GridFunction> ProlongateToMaxOrder() const;

 protected:
+   void ProjectBdrCoefficientNormal(Coefficient *coeff, VectorCoefficient *vcoeff,
+                                    const Array<int> &attr);
+
   /** @brief Accumulates (depending on @a type) the values of @a coeff at all
       shared vdofs and counts in how many zones each vdof appears. */
   void AccumulateAndCountZones(Coefficient &coeff, AvgType type,
@@ -500,6 +569,70 @@ protected:
   /// P-refinement version of Update().
   void UpdatePRef();

+   /** @brief Estimate the minimum value of the GridFunction in element @a elem
+    *  if it is below a certain @a min_threshold.
+    *
+    *  @details For a given element \p elem and grid function component \p vdim
+    *  an estimate of the function minimum is the minimum of the piecewise
+    *  linear lower bound obtained using the given PLBound object. The actual
+    *  minimum is between [minimum lower bound, minimum upper bound]. We
+    *  improve the estimate of the function minimum by recursively
+    *  subdividing the interval with the lowest lower bound, and computing
+    *  bounds on the sub-intervals.
+    *  This process continues until (i) the maximum recursion depth is reached
+    *  or (ii) the difference between the minimum upper bound and minimum lower
+    *  bound is less than a certain tolerance (\p tol * [initial maximum
+    *  upper bound - initial minimum lower bound]).
+    *  The function also terminates if the lowest minima estimate is found
+    *  to be above the given threshold \p min_threshold. This is useful when
+    *  we are interested in computing the global minimum of the function
+    *  over all elements. In this case we can reject elements where the lowest
+    *  bound is above the current global minimum. In case the function
+    *  minimum on the element is below the global minimum, we update
+    *  \p min_threshold.
+    *
+    *  We return a pair of values that bracket the actual minimum, i.e.
+    *  [min_lower_bound, min_upper_bound].
+    */
+   std::pair<real_t,real_t> EstimateFunctionMinimum(const int elem,
+                                                    const PLBound &plb,
+                                                    const int vdim,
+                                                    const int max_depth,
+                                                    const real_t tol,
+                                                    real_t &min_threshold)const;
+
+   /** @brief Estimate the maximum value of the GridFunction in element @a elem
+    *  if it is below a certain @a max_threshold.
+    *
+    *  @details For a given element \p elem and grid function component \p vdim
+    *  an estimate of the function maximum is the maximum of the piecewise
+    *  linear upper bound obtained using the given PLBound object. The actual
+    *  maximum is between [maximum lower bound, maximum upper bound]. We
+    *  improve the estimate of the function maximum by recursively
+    *  subdividing the interval with the highest upper bound, and computing
+    *  bounds on the sub-intervals.
+    *  This process continues until (i) the maximum recursion depth is reached
+    *  or (ii) the difference between the maximum upper bound and maximum lower
+    *  bound is less than a certain tolerance (\p tol * [initial maximum
+    *  upper bound - initial maximum lower bound]).
+    *  The function also terminates if the highest maxima estimate is found
+    *  to be below the given threshold \p max_threshold. This is useful when
+    *  we are interested in computing the global maximum of the function
+    *  over all elements. In this case we can reject elements where the upper
+    *  bound is below the current global maximum. In case the function
+    *  maximum on the element is above the global maximum, we update
+    *  \p max_threshold.
+    *
+    *  We return a pair of values that bracket the actual maximum, i.e.
+    *  [max_lower_bound, max_upper_bound].
+    */
+   std::pair<real_t,real_t> EstimateFunctionMaximum(const int elem,
+                                                    const PLBound &plb,
+                                                    const int vdim,
+                                                    const int max_depth,
+                                                    const real_t tol,
+                                                    real_t &max_threshold)const;
+
 public:
   /** @brief For each vdof, counts how many elements contain the vdof,
       as containment is determined by FiniteElementSpace::GetElementVDofs(). */
@@ -528,15 +661,26 @@ public:
   virtual void ProjectBdrCoefficient(Coefficient *coeff[],
                                      const Array<int> &attr);

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

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

@@ -1598,21 +1742,21 @@ public:
    */
   ///@{
   /// Computes the \ref PLBound for the gridfunction with number of control
-   /// points based on @a ref_factor, and returns the overall bounds for each
-   /// vdim (across all elements) in @b lower and @b upper. We also return the
+   /// points based on \p ref_factor, and returns the overall bounds for each
+   /// vdim (across all elements) in \p lower and \p upper. We also return the
   /// PLBound object used to compute the bounds.
-   /// We compute the bounds for each vdim if @a vdim < 1.
+   /// We compute the bounds for each vdim if \p vdim < 1.
   /// Note: For most cases, this method/interface will be sufficient.
   virtual PLBound GetBounds(Vector &lower, Vector &upper,
                             const int ref_factor=1, const int vdim=-1) const;

   /// Computes the \ref PLBound for the gridfunction with number of control
-   /// points based on @a ref_factor, and returns the bounds for each element
-   /// ordered byVDim:
+   /// points based on \p ref_factor, and returns the bounds for each element
+   /// ordered byNodes:
   /// lower_{0,0}, lower_{1,0}, ..., lower_{ne-1,0},
   /// lower_{0,1}, ..., lower_{ne-1,vdim-1}. We also return the
   /// PLBound object used to compute the bounds.
-   /// We compute the bounds for each vdim if @a vdim < 1.
+   /// We compute the bounds for each vdim if \p vdim < 1.
   PLBound GetElementBounds(Vector &lower, Vector &upper,
                            const int ref_factor=1, const int vdim=-1) const;

@@ -1623,6 +1767,18 @@ public:
                                        Vector &lower, Vector &upper,
                                        const int vdim = -1) const;

+   /** @brief Gets the bounds on given reference range inside an element.
+    *
+    *  @details @a ref_range is a vector of size 2*dim that specifies the
+    *  lower and upper limits in each dimension of the reference element.
+    *  For example, in 2D, ref_range = [rmin, smin, rmax, smax].
+    */
+   void GetElementBoundsAtControlPoints(const int elem, const PLBound &plb,
+                                        const Vector &ref_range,
+                                        const int vdim,
+                                        Vector &lower, Vector &upper,
+                                        Vector &control_pos) const;
+
   /// Compute bounds on the grid function for the given element.
   /// The bounds are stored in @b lower and @b upper.
   void GetElementBounds(const int elem, const PLBound &plb,
@@ -1630,11 +1786,45 @@ public:
                         const int vdim = -1) const;

   /// Compute bounds on the grid function for all the elements. The bounds
-   /// are returned in @b lower and @b upper, ordered byVDim:
+   /// are returned in @b lower and @b upper, ordered byNodes:
   /// lower_{0,0}, lower_{1,0}, ..., lower_{ne-1,0},
   /// lower_{0,1}, ..., lower_{ne-1,vdim-1}
   void GetElementBounds(const PLBound &plb, Vector &lower, Vector &upper,
                         const int vdim=-1) const;
+
+   /** @brief Estimate the minimum value of the GridFunction in element @a elem.
+    *
+    *  @details See the protected version of EstimateFunctionMinimum for
+    *  details.
+    */
+   std::pair<real_t, real_t> EstimateFunctionMinimum(const int elem,
+                                                     const PLBound &plb,
+                                                     const int vdim,
+                                                     const int max_depth,
+                                                     const real_t tol) const;
+
+   /** @brief Estimate the minimum value of the GridFunction in element @a elem.
+    *
+    *  @details See the protected version of EstimateFunctionMaximum for
+    *  details.
+    */
+   std::pair<real_t, real_t> EstimateFunctionMaximum(const int elem,
+                                                     const PLBound &plb,
+                                                     const int vdim,
+                                                     const int max_depth,
+                                                     const real_t tol) const;
+
+   /** @brief Estimate the GridFunction minimum across all elements. */
+   virtual std::pair<real_t,real_t> EstimateFunctionMinimum(const int vdim,
+                                                            const PLBound &plb,
+                                                            const int max_depth,
+                                                            const real_t tol) const;
+
+   /** @brief Estimate the GridFunction maximum across all elements. */
+   virtual std::pair<real_t,real_t> EstimateFunctionMaximum(const int vdim,
+                                                            const PLBound &plb,
+                                                            const int max_depth,
+                                                            const real_t tol) const;
   ///@}

   /// Destroys grid function.
@@ -1740,7 +1930,7 @@ real_t ComputeElementLpDistance(real_t p, int i,
                                GridFunction& gf1, GridFunction& gf2);


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

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

@@ -234,7 +234,7 @@ void FindPointsGSLIB::Setup(Mesh &m, const double bb_t, const double newt_tol,
 }

 void FindPointsGSLIB::FindPoints(const Vector &point_pos,
-                                 int point_pos_ordering)
+                                 const int point_pos_ordering)
 {
   MFEM_VERIFY(setupflag, "Use FindPointsGSLIB::Setup before finding points.");
   bool dev_mode = (point_pos.UseDevice() && Device::IsEnabled());
@@ -482,7 +482,7 @@ void FindPointsGSLIB::SetupDevice()
 }

 void FindPointsGSLIB::FindPointsOnDevice(const Vector &point_pos,
-                                         int point_pos_ordering)
+                                         const int point_pos_ordering)
 {
   if (!DEV.setup_device)
   {
@@ -505,13 +505,13 @@ void FindPointsGSLIB::FindPointsOnDevice(const Vector &point_pos,

   if (dim == 2)
   {
-      FindPointsLocal2(point_pos, point_pos_ordering, gsl_code, gsl_elem, gsl_ref,
-                       gsl_dist, points_cnt);
+      FindPointsLocal2(point_pos, point_pos_ordering, gsl_code, gsl_elem,
+                       gsl_ref, gsl_dist, points_cnt);
   }
   else
   {
-      FindPointsLocal3(point_pos, point_pos_ordering, gsl_code, gsl_elem, gsl_ref,
-                       gsl_dist, points_cnt);
+      FindPointsLocal3(point_pos, point_pos_ordering, gsl_code, gsl_elem,
+                       gsl_ref, gsl_dist, points_cnt);
   }

   // Sync from device to host
@@ -1085,7 +1085,7 @@ void FindPointsGSLIB::InterpolateOnDevice(const Vector &field_in_evec,
 #else
 void FindPointsGSLIB::SetupDevice() {};
 void FindPointsGSLIB::FindPointsOnDevice(const Vector &point_pos,
-                                         int point_pos_ordering) {};
+                                         const int point_pos_ordering) {};
 void FindPointsGSLIB::InterpolateOnDevice(const Vector &field_in_evec,
                                          Vector &field_out,
                                          const int nel, const int ncomp,
@@ -1094,7 +1094,8 @@ void FindPointsGSLIB::InterpolateOnDevice(const Vector &field_in_evec,
 #endif

 void FindPointsGSLIB::FindPoints(Mesh &m, const Vector &point_pos,
-                                 int point_pos_ordering, const double bb_t,
+                                 const int point_pos_ordering,
+                                 const double bb_t,
                                 const double newt_tol, const int npt_max)
 {
   if (!setupflag || (mesh != &m) )
@@ -1105,16 +1106,28 @@ void FindPointsGSLIB::FindPoints(Mesh &m, const Vector &point_pos,
 }

 void FindPointsGSLIB::Interpolate(const Vector &point_pos,
-                                  const GridFunction &field_in, Vector &field_out,
-                                  int point_pos_ordering)
+                                  const GridFunction &field_in,
+                                  Vector &field_out,
+                                  const int point_pos_ordering)
 {
   FindPoints(point_pos, point_pos_ordering);
   Interpolate(field_in, field_out);
 }

+void FindPointsGSLIB::Interpolate(const Vector &point_pos,
+                                  const GridFunction &field_in,
+                                  Vector &field_out,
+                                  const int point_pos_ordering,
+                                  const int field_out_ordering)
+{
+   FindPoints(point_pos, point_pos_ordering);
+   Interpolate(field_in, field_out, field_out_ordering);
+}
+
 void FindPointsGSLIB::Interpolate(Mesh &m, const Vector &point_pos,
-                                  const GridFunction &field_in, Vector &field_out,
-                                  int point_pos_ordering)
+                                  const GridFunction &field_in,
+                                  Vector &field_out,
+                                  const int point_pos_ordering)
 {
   FindPoints(m, point_pos, point_pos_ordering);
   Interpolate(field_in, field_out);
@@ -1470,7 +1483,7 @@ void FindPointsGSLIB::SetupSplitMeshesAndIntegrationRules(const int order)
 }

 void FindPointsGSLIB::GetNodalValues(const GridFunction *gf_in,
-                                     Vector &node_vals)
+                                     Vector &node_vals) const
 {
   const GridFunction *nodes     = gf_in;
   const FiniteElementSpace *fes = nodes->FESpace();
@@ -1758,6 +1771,13 @@ void FindPointsGSLIB::MapRefPosAndElemIndices()

 void FindPointsGSLIB::Interpolate(const GridFunction &field_in,
                                  Vector &field_out)
+{
+   Interpolate(field_in, field_out, field_in.FESpace()->GetOrdering());
+}
+
+void FindPointsGSLIB::Interpolate(const GridFunction &field_in,
+                                  Vector &field_out,
+                                  const int field_out_ordering)
 {
   const int  gf_order   = field_in.FESpace()->GetMaxElementOrder(),
              mesh_order = mesh->GetNodalFESpace()->GetMaxElementOrder();
@@ -1800,7 +1820,7 @@ void FindPointsGSLIB::Interpolate(const GridFunction &field_in,
      const int maxOrder = field_in.FESpace()->GetMaxElementOrder();

      InterpolateOnDevice(node_vals, field_out, NE_split_total, ncomp,
-                          maxOrder+1, field_in.FESpace()->GetOrdering());
+                          maxOrder+1, field_out_ordering);
      return;
 #endif
   }
@@ -1812,12 +1832,13 @@ void FindPointsGSLIB::Interpolate(const GridFunction &field_in,
       field_in.FESpace()->IsVariableOrder() ==
       mesh->GetNodalFESpace()->IsVariableOrder())
   {
-      InterpolateH1(field_in, field_out);
+      InterpolateH1(field_in, field_out, field_out_ordering);
      return;
   }
   else
   {
-      InterpolateGeneral(field_in, field_out);
+      InterpolateGeneral(field_in, field_out,
+                         field_out_ordering);
      if (!fec_l2 || avgtype == AvgType::NONE) { return; }
   }

@@ -1861,11 +1882,11 @@ void FindPointsGSLIB::Interpolate(const GridFunction &field_in,

      if (gf_order_h1 == mesh_order) // basis is GaussLobatto by default
      {
-         InterpolateH1(field_in_h1, field_out_l2);
+         InterpolateH1(field_in_h1, field_out_l2, field_out_ordering);
      }
      else
      {
-         InterpolateGeneral(field_in_h1, field_out_l2);
+         InterpolateGeneral(field_in_h1, field_out_l2, field_out_ordering);
      }

      // Copy interpolated values for the points on element border
@@ -1873,7 +1894,7 @@ void FindPointsGSLIB::Interpolate(const GridFunction &field_in,
      {
         for (int i = 0; i < indl2.Size(); i++)
         {
-            int idx = field_in_h1.FESpace()->GetOrdering() == Ordering::byNODES?
+            int idx = field_out_ordering == Ordering::byNODES?
                      indl2[i] + j*points_cnt:
                      indl2[i]*ncomp + j;
            field_out(idx) = field_out_l2(idx);
@@ -1883,7 +1904,8 @@ void FindPointsGSLIB::Interpolate(const GridFunction &field_in,
 }

 void FindPointsGSLIB::InterpolateH1(const GridFunction &field_in,
-                                    Vector &field_out)
+                                    Vector &field_out,
+                                    const int field_out_ordering)
 {
   FiniteElementSpace ind_fes(mesh, field_in.FESpace()->FEColl());
   if (field_in.FESpace()->IsVariableOrder())
@@ -1913,7 +1935,8 @@ void FindPointsGSLIB::InterpolateH1(const GridFunction &field_in,
                dataptrout = i*points_cnt;
      if (field_in.FESpace()->GetOrdering() == Ordering::byNODES)
      {
-         field_in_scalar.NewDataAndSize(field_in.GetData()+dataptrin, points_fld);
+         field_in_scalar.NewDataAndSize(field_in.GetData()+dataptrin,
+                                        points_fld);
      }
      else
      {
@@ -1945,7 +1968,7 @@ void FindPointsGSLIB::InterpolateH1(const GridFunction &field_in,
                        (gslib::findpts_data_3 *)this->fdataD);
      }
   }
-   if (field_in.FESpace()->GetOrdering() == Ordering::byVDIM)
+   if (field_out_ordering == Ordering::byVDIM)
   {
      Vector field_out_temp = field_out;
      for (int i = 0; i < ncomp; i++)
@@ -1959,7 +1982,8 @@ void FindPointsGSLIB::InterpolateH1(const GridFunction &field_in,
 }

 void FindPointsGSLIB::InterpolateGeneral(const GridFunction &field_in,
-                                         Vector &field_out)
+                                         Vector &field_out,
+                                         const int field_out_ordering)
 {
   int ncomp   = field_in.VectorDim(),
       nptorig = points_cnt,
@@ -1979,7 +2003,7 @@ void FindPointsGSLIB::InterpolateGeneral(const GridFunction &field_in,
         if (dim == 3) { ip.z = gsl_mfem_ref(index*dim + 2); }
         Vector localval(ncomp);
         field_in.GetVectorValue(gsl_mfem_elem[index], ip, localval);
-         if (field_in.FESpace()->GetOrdering() == Ordering::byNODES)
+         if (field_out_ordering == Ordering::byNODES)
         {
            for (int i = 0; i < ncomp; i++)
            {
@@ -2014,7 +2038,10 @@ void FindPointsGSLIB::InterpolateGeneral(const GridFunction &field_in,
      for (int index = 0; index < npt; index++)
      {
         if (gsl_code[index] == 2) { continue; }
-         for (int d = 0; d < dim; ++d) { pt->r[d]= gsl_mfem_ref(index*dim + d); }
+         for (int d = 0; d < dim; ++d)
+         {
+            pt->r[d]= gsl_mfem_ref(index*dim + d);
+         }
         pt->index = index;
         pt->proc  = gsl_proc[index];
         pt->el    = gsl_mfem_elem[index];
@@ -2104,7 +2131,7 @@ void FindPointsGSLIB::InterpolateGeneral(const GridFunction &field_in,
            sdpt = (struct send_pt *)sendpt->ptr;
            for (int index = 0; index < static_cast<int>(sendpt->n); index++)
            {
-               int idx = field_in.FESpace()->GetOrdering() == Ordering::byNODES ?
+               int idx = field_out_ordering == Ordering::byNODES ?
                         sdpt->index + j*nptorig :
                         sdpt->index*ncomp + j;
               field_out(idx) = sdpt->ival;
@@ -2246,7 +2273,7 @@ void FindPointsGSLIB::DistributeInterpolatedValues(const Vector &int_vals,
   }
 }

-void FindPointsGSLIB::GetAxisAlignedBoundingBoxes(Vector &aabb)
+void FindPointsGSLIB::GetAxisAlignedBoundingBoxes(Vector &aabb) const
 {
   MFEM_VERIFY(setupflag, "Call FindPointsGSLIB::Setup method first");
   auto *findptsData3 = (gslib::findpts_data_3 *)this->fdataD;
@@ -2317,7 +2344,7 @@ void FindPointsGSLIB::GetAxisAlignedBoundingBoxes(Vector &aabb)
 }

 void FindPointsGSLIB::GetOrientedBoundingBoxes(DenseTensor &obbA, Vector &obbC,
-                                               Vector &obbV)
+                                               Vector &obbV) const
 {
   MFEM_VERIFY(setupflag, "Call FindPointsGSLIB::Setup method first");
   auto *findptsData3 = (gslib::findpts_data_3 *)this->fdataD;
@@ -2502,8 +2529,8 @@ void OversetFindPointsGSLIB::Setup(Mesh &m, const int meshid,
 }

 void OversetFindPointsGSLIB::FindPoints(const Vector &point_pos,
-                                        Array<unsigned int> &point_id,
-                                        int point_pos_ordering)
+                                        const Array<unsigned int> &point_id,
+                                        const int point_pos_ordering)
 {
   MFEM_VERIFY(setupflag, "Use OversetFindPointsGSLIB::Setup before "
               "finding points.");
@@ -2582,10 +2609,10 @@ void OversetFindPointsGSLIB::FindPoints(const Vector &point_pos,
 }

 void OversetFindPointsGSLIB::Interpolate(const Vector &point_pos,
-                                         Array<unsigned int> &point_id,
+                                         const Array<unsigned int> &point_id,
                                         const GridFunction &field_in,
                                         Vector &field_out,
-                                         int point_pos_ordering)
+                                         const int point_pos_ordering)
 {
   FindPoints(point_pos, point_id, point_pos_ordering);
   Interpolate(field_in, field_out);
@@ -119,11 +119,13 @@ protected:
   } DEV;

   /// Use GSLIB for communication and interpolation
-   virtual void InterpolateH1(const GridFunction &field_in, Vector &field_out);
+   virtual void InterpolateH1(const GridFunction &field_in, Vector &field_out,
+                              const int field_out_ordering);
   /// Uses GSLIB Crystal Router for communication followed by MFEM's
   /// interpolation functions
   virtual void InterpolateGeneral(const GridFunction &field_in,
-                                   Vector &field_out);
+                                   Vector &field_out,
+                                   const int field_out_ordering);

   /// Since GSLIB is designed to work with quads/hexes, we split every
   /// triangle/tet/prism/pyramid element into quads/hexes.
@@ -140,7 +142,7 @@ protected:
   virtual void SetupSplitMeshesAndIntegrationRules(const int order);

   /// Get GridFunction value at the points expected by GSLIB.
-   virtual void GetNodalValues(const GridFunction *gf_in, Vector &node_vals);
+   virtual void GetNodalValues(const GridFunction *gf_in, Vector &node_vals) const;

   /// Map {r,s,t} coordinates from [-1,1] to [0,1] for MFEM. For simplices,
   /// find the original element number (that was split into micro quads/hexes)
@@ -182,7 +184,7 @@ protected:
       These positions can be ordered byNodes: (XXX...,YYY...,ZZZ) or
       byVDim: (XYZ,XYZ,....XYZ) specified by @a point_pos_ordering. */
   void FindPointsOnDevice(const Vector &point_pos,
-                           int point_pos_ordering = Ordering::byNODES);
+                           const int point_pos_ordering = Ordering::byNODES);

   /** Interpolation of field values at prescribed reference space positions.
       @param[in] field_in_evec E-vector of grid function to be interpolated.
@@ -253,10 +255,15 @@ public:
       #gsl_dist        Distance between the sought and the found point
                        in physical space. */
   void FindPoints(const Vector &point_pos,
-                   int point_pos_ordering = Ordering::byNODES);
+                   const int point_pos_ordering = Ordering::byNODES);
+   /// Convenience function when point positions are in a ParticleVector
+   void FindPoints(const ParticleVector &point_pos)
+   {
+      FindPoints(point_pos, point_pos.GetOrdering());
+   }
   /// Setup FindPoints and search positions
   void FindPoints(Mesh &m, const Vector &point_pos,
-                   int point_pos_ordering = Ordering::byNODES,
+                   const int point_pos_ordering = Ordering::byNODES,
                   const double bb_t = 0.1, const double newt_tol = 1.0e-12,
                   const int npt_max = 256);

@@ -266,20 +273,28 @@ public:
                              \p field_in is in H1 and in the same space as the
                              mesh that was given to Setup().
       @param[out] field_out  Interpolated values. For points that are not found
-                              the value is set to #default_interp_value. */
+                              the value is set to #default_interp_value.
+                              The output ordering is determined from field_in.*/
   virtual void Interpolate(const GridFunction &field_in, Vector &field_out);
+   /// Interpolation of field values, with output ordering specification.
+   virtual void Interpolate(const GridFunction &field_in, Vector &field_out,
+                            const int field_out_ordering);
   /** Search positions and interpolate. The ordering (byNODES or byVDIM) of
       the output values in \p field_out corresponds to the ordering used
       in the input GridFunction \p field_in. */
   void Interpolate(const Vector &point_pos, const GridFunction &field_in,
                    Vector &field_out,
-                    int point_pos_ordering = Ordering::byNODES);
+                    const int point_pos_ordering = Ordering::byNODES);
+   /// Search positions and interpolate with given point and output ordering.
+   void Interpolate(const Vector &point_pos, const GridFunction &field_in,
+                    Vector &field_out, const int point_pos_ordering,
+                    const int field_out_ordering);
   /** Setup FindPoints, search positions and interpolate. The ordering (byNODES
       or byVDIM) of the output values in \p field_out corresponds to the
       ordering used in the input GridFunction \p field_in. */
   void Interpolate(Mesh &m, const Vector &point_pos,
                    const GridFunction &field_in, Vector &field_out,
-                    int point_pos_ordering = Ordering::byNODES);
+                    const int point_pos_ordering = Ordering::byNODES);

   /// Average type to be used for L2 functions in-case a point is located at
   /// an element boundary where the function might be multi-valued.
@@ -376,7 +391,7 @@ public:
   /// The size of the returned vector is (nel x nverts x dim), where nel is the
   /// number of elements (after splitting for simplcies), nverts is number of
   /// vertices (4 in 2D, 8 in 3D), and dim is the spatial dimension.
-   void GetAxisAlignedBoundingBoxes(Vector &aabb);
+   void GetAxisAlignedBoundingBoxes(Vector &aabb) const;

   /// Return the oriented bounding boxes (OBB) computed during \ref Setup.
   /// Each OBB is represented using the inverse transformation (A^{-1}) and
@@ -386,7 +401,8 @@ public:
   /// size (dim x dim x nel), and the OBB centers are returned in \p obbC,
   /// a vector of size (nel x dim). The vertices of the OBBs are returned in
   /// \p obbV, a vector of size (nel x nverts x dim) .
-   void GetOrientedBoundingBoxes(DenseTensor &obbA, Vector &obbC, Vector &obbV);
+   void GetOrientedBoundingBoxes(DenseTensor &obbA, Vector &obbC,
+                                 Vector &obbV) const;
 };

 /** \brief OversetFindPointsGSLIB enables use of findpts for arbitrary number of
@@ -446,13 +462,14 @@ public:
                                       byNodes: (XXX...,YYY...,ZZZ) or
                                       byVDim: (XYZ,XYZ,....XYZ) */
   void FindPoints(const Vector &point_pos,
-                   Array<unsigned int> &point_id,
-                   int point_pos_ordering = Ordering::byNODES);
+                   const Array<unsigned int> &point_id,
+                   const int point_pos_ordering = Ordering::byNODES);

   /** Search positions and interpolate */
-   void Interpolate(const Vector &point_pos, Array<unsigned int> &point_id,
+   void Interpolate(const Vector &point_pos,
+                    const Array<unsigned int> &point_id,
                    const GridFunction &field_in, Vector &field_out,
-                    int point_pos_ordering = Ordering::byNODES);
+                    const int point_pos_ordering = Ordering::byNODES);
   using FindPointsGSLIB::Interpolate;
 };

@@ -789,7 +789,6 @@ void Hybridization::ComputeH()
   }
   else
   {
-      // TODO: add ones on the diagonal of zero rows
      V->Finalize();
      Array<HYPRE_BigInt> V_J(V->NumNonZeroElems());
      MFEM_ASSERT(c_pfes, "");
@@ -823,6 +822,13 @@ void Hybridization::ComputeH()
      MFEM_VERIFY(pH.Type() != Operator::PETSC_MATIS, "To be implemented");
      pH.MakePtAP(plpH, pP);
      delete lpH;
+
+      HypreParMatrix *hH = pH.As<HypreParMatrix>();
+      MFEM_ASSERT(hH, "");
+
+      SparseMatrix H_diag;
+      hH->GetDiag(H_diag);
+      H_diag.SetDiagIdentity();
   }
 #endif
 }
@@ -14,8 +14,11 @@

 #include "../config/config.hpp"
 #include "../general/array.hpp"
+#include "../linalg/operator.hpp"
 #include "../linalg/vector.hpp"

+#include <memory>
+
 namespace mfem
 {

@@ -45,15 +48,30 @@ protected:
   Array<int> hat_dof_gather_map;
   Array<DofType> hat_dof_marker;

-   Array<int> el_to_face;
-   Array<int> face_to_el;
+   Array<int> el_to_face; ///< Element to face connectivity.
+   Array<int> el_face_offsets; ///< Per-element offsets into @a el_to_face.
+   Array<int> face_to_el; ///< Face-to-element connectivity.
+   Array<int> face_face_offsets; ///< Face-to-face offsets.
+
+   int n_el_face; ///< Total number of element-to-face connections.
+   int n_face_face; ///< Total number of face-to-face connections.
+
   Vector Ct_mat; ///< Constraint matrix (transposed) stored element-wise.

+   /// @name For parallel non-conforming meshes
+   ///@{
+   std::unique_ptr<Operator> P_pc; ///< Partially conforming prolongation.
+   std::unique_ptr<Operator> P_nbr; ///< Face-neighbor prolongation.
+   ///@}
+
   Array<int> idofs, bdofs;

   Vector Ahat, Ahat_ii, Ahat_ib, Ahat_bi, Ahat_bb;
   Array<int> Ahat_ii_piv, Ahat_bb_piv;

+   /// Return the (partially) conforming prolongation on the constraint space.
+   const Operator &GetProlongation() const;
+
 public:
   /// Construct the constraint matrix.
   void ConstructC();
@@ -1004,13 +1004,16 @@ inline void SmemPADiffusionApply3D(const int NE,
   const int max_d1d = T_D1D ? T_D1D : DeviceDofQuadLimits::Get().MAX_D1D;
   MFEM_VERIFY(D1D <= max_d1d, "");
   MFEM_VERIFY(Q1D <= max_q1d, "");
-   auto b = Reshape(b_.Read(), Q1D, D1D);
-   auto g = Reshape(g_.Read(), Q1D, D1D);
-   auto d = Reshape(d_.Read(), Q1D, Q1D, Q1D, symmetric ? 6 : 9, NE);
-   auto x = Reshape(x_.Read(), D1D, D1D, D1D, NE);
+   const auto b = Reshape(b_.Read(), Q1D, D1D);
+   const auto g = Reshape(g_.Read(), Q1D, D1D);
+   const auto d = Reshape(d_.Read(), Q1D, Q1D, Q1D, symmetric ? 6 : 9, NE);
+   const auto x = Reshape(x_.Read(), D1D, D1D, D1D, NE);
   auto y = Reshape(y_.ReadWrite(), D1D, D1D, D1D, NE);
   MFEM_VERIFY(D1D <= Q1D, "THREAD_DIRECT requires D1D <= Q1D");
-   mfem::forall_3D(NE, Q1D, Q1D, Q1D, [=] MFEM_HOST_DEVICE (int e)
+
+   mfem::forall_3D<T_Q1D*T_Q1D*T_Q1D>(NE,
+                                      Q1D, Q1D, Q1D,
+                                      [=] MFEM_HOST_DEVICE (int e)
   {
      const int D1D = T_D1D ? T_D1D : d1d;
      const int Q1D = T_Q1D ? T_Q1D : q1d;
@@ -197,15 +197,21 @@ static void EAHdivAssemble3D(const int NE,
      // Assemble (one row per thread)
      MFEM_FOREACH_THREAD(idx_i, x, NDOF)
      {
+         // NOTE: due to an llvm backend bug, usage of the modulus operator
+         // has been removed from this foreach section.
         const int ic = idx_i / NDOF_C;
-         const int idx_ii = idx_i % NDOF_C;
+         const int idx_ii = idx_i - ic * NDOF_C; // idx_i % NDOF_C

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

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

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

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

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

            const real_t (&Bj1)[MQ1][MD1] = (jc == 0) ? r_Bc : r_Bo;
            const real_t (&Bj2)[MQ1][MD1] = (jc == 1) ? r_Bc : r_Bo;
@@ -1133,11 +1133,11 @@ inline void SmemPAMassApply3D(const int NE,
   const int max_d1d = T_D1D ? T_D1D : DeviceDofQuadLimits::Get().MAX_D1D;
   MFEM_VERIFY(D1D <= max_d1d, "");
   MFEM_VERIFY(Q1D <= max_q1d, "");
-   auto b = b_.Read();
-   auto d = d_.Read();
-   auto x = x_.Read();
+   const auto b = b_.Read();
+   const auto d = d_.Read();
+   const auto x = x_.Read();
   auto y = y_.ReadWrite();
-   mfem::forall_2D(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE (int e)
+   mfem::forall_2D<T_Q1D*T_Q1D>(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE (int e)
   {
      internal::SmemPAMassApply3D_Element<T_D1D,T_Q1D>(e, NE, b, d, x, y, d1d, q1d);
   });
@@ -1156,8 +1156,8 @@ inline void EAMassAssemble1D(const int NE,
   const int Q1D = T_Q1D ? T_Q1D : q1d;
   MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().MAX_D1D, "");
   MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().MAX_Q1D, "");
-   auto B = Reshape(basis.Read(), Q1D, D1D);
-   auto D = Reshape(padata.Read(), Q1D, NE);
+   const auto B = Reshape(basis.Read(), Q1D, D1D);
+   const auto D = Reshape(padata.Read(), Q1D, NE);
   auto M = Reshape(add ? eadata.ReadWrite() : eadata.Write(), D1D, D1D, NE);
   mfem::forall_2D(NE, D1D, D1D, [=] MFEM_HOST_DEVICE (int e)
   {
@@ -28,7 +28,7 @@ void NormalTraceJumpIntegrator::AssembleEAInteriorFaces(
   const FaceType ftype = FaceType::Interior;
   const int nf = mesh.GetNFbyType(ftype);

-   const Geometry::Type geom = mesh.GetFaceGeometry(0);
+   const Geometry::Type geom = mesh.GetTypicalFaceGeometry();
   const int trial_order = trial_fes.GetMaxElementOrder();
   const int test_order = test_fes.GetMaxElementOrder();
   const int qorder = test_order + trial_order - 1;
@@ -47,7 +47,7 @@ void NormalTraceJumpIntegrator::AssembleEAInteriorFaces(
      });
   }

-   const FiniteElement &trial_face_el = *trial_fes.GetFaceElement(0);
+   const FiniteElement &trial_face_el = *trial_fes.GetTypicalTraceElement();
   const auto maps = &trial_face_el.GetDofToQuad(ir, DofToQuad::TENSOR);
   const int ndof_face = trial_face_el.GetDof();

@@ -72,7 +72,7 @@ void NormalTraceJumpIntegrator::AssembleEAInteriorFaces(
      MFEM_ABORT("Unknown kernel.");
   }

-   const FiniteElement &test_el = *test_fes.GetFE(0);
+   const FiniteElement &test_el = *test_fes.GetTypicalFE();
   const int n_faces_per_el = 2*dim; // assuming tensor product
   // Get all the local face maps (mapping from lexicographic face index to
   // lexicographic volume index, depending on the local face index).
@@ -90,10 +90,10 @@ void NormalTraceJumpIntegrator::AssembleEAInteriorFaces(
   Array<int> face_info(nf * 4);
   {
      int fidx = 0;
-      for (int f = 0; f < mesh.GetNumFaces(); ++f)
+      for (int f = 0; f < mesh.GetNumFacesWithGhost(); ++f)
      {
         Mesh::FaceInformation finfo = mesh.GetFaceInformation(f);
-         if (!finfo.IsInterior()) { continue; }
+         if (!finfo.IsInterior() || finfo.IsNonconformingCoarse()) { continue; }
         face_info[0 + fidx*4] = finfo.element[0].local_face_id;
         face_info[1 + fidx*4] = finfo.element[0].orientation;
         face_info[2 + fidx*4] = finfo.element[1].local_face_id;
@@ -114,7 +114,7 @@ void NormalTraceJumpIntegrator::AssembleEAInteriorFaces(
   else
   {
      d_emat = emat.Write();
-      mfem::forall(emat.Size(), [=] MFEM_HOST_DEVICE (int i) { d_emat[i] = 0.0; });
+      emat = 0.0; // Will execute on device, since Write() sets the device flag
   }

   const auto face_mats = Reshape(mass_emat.Read(), ndof_face, ndof_face, nf);
@@ -133,26 +133,104 @@ void NormalTraceJumpIntegrator::AssembleEAInteriorFaces(
      }
   };

-   mfem::forall_3D(nf, ndof_face, ndof_face, 2, [=] MFEM_HOST_DEVICE (int f)
+   auto permute_face_2 = [=] MFEM_HOST_DEVICE(int local_face_1, int local_face_2,
+                                              int orient, int size1d, int index)
   {
-      MFEM_FOREACH_THREAD(el_i, z, 2)
+      if (dim == 2)
      {
-         const int lf_i = d_face_info(0, el_i, f);
-         const int orient = d_face_info(1, el_i, f);
-         // Loop over face indices in "native ordering"
-         MFEM_FOREACH_THREAD(i_lex, x, ndof_face)
+         return internal::PermuteFace2D(local_face_1, local_face_2, orient,
+                                        size1d, index);
+      }
+      else // dim == 3
+      {
+         return internal::PermuteFace3D(local_face_1, local_face_2, orient,
+                                        size1d, index);
+      }
+   };
+
+   if (mesh.Conforming())
+   {
+      mfem::forall_3D(nf, ndof_face, ndof_face, 2, [=] MFEM_HOST_DEVICE (int f)
+      {
+         MFEM_FOREACH_THREAD(el_i, z, 2)
         {
-            // Convert to lexicographic relative to the face itself
-            const int i_face = permute_face(lf_i, orient, d1d, i_lex);
-            // Convert from lexicographic face DOF to volume DOF
-            const int i = d_face_maps(i_lex, lf_i);
-            MFEM_FOREACH_THREAD(j, y, ndof_face)
+            const int lf_i = d_face_info(0, el_i, f);
+            const int orient = d_face_info(1, el_i, f);
+            // Loop over face indices in "native ordering"
+            MFEM_FOREACH_THREAD(i_lex, x, ndof_face)
            {
-               el_mats(i, j, el_i, f) += face_mats(i_face, j, f);
+               // Convert to lexicographic relative to the face itself
+               const int i_face = permute_face(lf_i, orient, d1d, i_lex);
+               // Convert from lexicographic face DOF to volume DOF
+               const int i = d_face_maps(i_lex, lf_i);
+               MFEM_FOREACH_THREAD(j, y, ndof_face)
+               {
+                  el_mats(i, j, el_i, f) += face_mats(i_face, j, f);
+               }
            }
         }
-      }
-   });
+      });
+   }
+   else
+   {
+      const InterpolationManager &interp =
+         test_fes.GetInterpolationManager(ElementDofOrdering::LEXICOGRAPHIC, ftype);
+
+      auto interp_configs = interp.GetFaceInterpConfig().Read();
+      const int nc_size = interp.GetNumInterpolators();
+      auto d_interp = Reshape(interp.GetInterpolators().Read(),
+                              ndof_face, ndof_face, nc_size);
+
+      mfem::forall(nf, [=] MFEM_HOST_DEVICE (int f)
+      {
+         const InterpConfig conf = interp_configs[f];
+         const int master_side = conf.master_side;
+         const int interp_index = conf.index;
+
+         const int lf_0 = d_face_info(0, 0, f);
+
+         for (int el_i = 0; el_i < 2; ++el_i)
+         {
+            const int lf_i = d_face_info(0, el_i, f);
+            const int orient = d_face_info(1, el_i, f);
+
+            for (int j = 0; j < ndof_face; j++)
+            {
+               for (int i_lex = 0; i_lex < ndof_face; i_lex++)
+               {
+                  real_t val = 0.0;
+                  if (conf.is_non_conforming && el_i == master_side)
+                  {
+                     // Interpolate from el_i (coarse element) to the fine face.
+                     // The mapping is given by d_interp, which uses indices
+                     // relative to element 0.
+
+                     // i0 is lexicographic relative to element 0
+                     const int i0 = permute_face_2(lf_i, lf_0, orient, d1d, i_lex);
+
+                     // k0 is lexicographic relative to element 0
+                     for (int k0 = 0; k0 < ndof_face; k0++)
+                     {
+                        // k is relative to the face itself
+                        const int k = permute_face(lf_0, orient, d1d, k0);
+                        val += d_interp(k0, i0, interp_index)
+                               * face_mats(k, j, f);
+                     }
+                  }
+                  else
+                  {
+                     // Convert to lexicographic relative to the face itself
+                     const int i_face = permute_face(lf_i, orient, d1d, i_lex);
+                     val = face_mats(i_face, j, f);
+                  }
+                  // Convert from lexicographic face DOF to volume DOF
+                  const int i = d_face_maps(i_lex, lf_i);
+                  el_mats(i, j, el_i, f) += val;
+               }
+            }
+         }
+      });
+   }
 }

 }
@@ -54,7 +54,7 @@ void SmemPAVectorDiffusionApply2D(const int NE,
   const auto XE = Reshape(x.Read(), D1D, D1D, SDIM, NE);
   auto YE = Reshape(y.ReadWrite(), D1D, D1D, SDIM, NE);

-   mfem::forall_2D(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE(int e)
+   mfem::forall_2D<T_Q1D*T_Q1D>(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE(int e)
   {
      constexpr int MD1 = T_D1D > 0 ? SetMaxOf(T_D1D) : DofQuadLimits::MAX_T1D;
      constexpr int MQ1 = T_Q1D > 0 ? SetMaxOf(T_Q1D) : DofQuadLimits::MAX_T1D;
@@ -120,7 +120,7 @@ void SmemPAVectorDiffusionApply3D(const int NE,
   const auto XE = Reshape(x.Read(), D1D, D1D, D1D, SDIM, NE);
   auto YE = Reshape(y.ReadWrite(), D1D, D1D, D1D, SDIM, NE);

-   mfem::forall_2D(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE(int e)
+   mfem::forall_2D<T_Q1D*T_Q1D>(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE(int e)
   {
      constexpr int MD1 = T_D1D > 0 ? SetMaxOf(T_D1D) : DofQuadLimits::MAX_T1D;
      constexpr int MQ1 = T_Q1D > 0 ? SetMaxOf(T_Q1D) : DofQuadLimits::MAX_T1D;
@@ -171,15 +171,15 @@ template<int DIM, int T_SDIM, int T_D1D, int T_Q1D>
 VectorDiffusionIntegrator::ApplyKernelType
 VectorDiffusionIntegrator::ApplyPAKernels::Kernel()
 {
-   if (DIM == 2)
+   if constexpr (DIM == 2)
   {
      return internal::SmemPAVectorDiffusionApply2D<T_SDIM, T_D1D, T_Q1D>;
   }
-   else if (DIM == 3)
+   else if constexpr (DIM == 3)
   {
      return internal::SmemPAVectorDiffusionApply3D<T_SDIM, T_D1D, T_Q1D>;
   }
-   else { MFEM_ABORT("Unsupported kernel"); }
+   MFEM_ABORT("Unsupported kernel");
 }

 inline VectorDiffusionIntegrator::ApplyKernelType
@@ -51,7 +51,7 @@ void SmemPAVectorMassApply2D(const int NE,
   const auto X = Reshape(x.Read(), D1D, D1D, VDIM, NE);
   auto Y = Reshape(y.ReadWrite(), D1D, D1D, VDIM, NE);

-   mfem::forall_2D(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE(int e)
+   mfem::forall_2D<T_Q1D*T_Q1D>(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE(int e)
   {
      constexpr int MD1 = T_D1D > 0 ? SetMaxOf(T_D1D) : DofQuadLimits::MAX_T1D;
      constexpr int MQ1 = T_Q1D > 0 ? SetMaxOf(T_Q1D) : DofQuadLimits::MAX_T1D;
@@ -119,7 +119,7 @@ void SmemPAVectorMassApply3D(const int NE,
   const auto X = Reshape(x.Read(), D1D, D1D, D1D, VDIM, NE);
   auto Y = Reshape(y.ReadWrite(), D1D, D1D, D1D, VDIM, NE);

-   mfem::forall_2D(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE(int e)
+   mfem::forall_2D<T_Q1D*T_Q1D>(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE(int e)
   {
      constexpr int MD1 = T_D1D > 0 ? SetMaxOf(T_D1D) : DofQuadLimits::MAX_T1D;
      constexpr int MQ1 = T_Q1D > 0 ? SetMaxOf(T_Q1D) : DofQuadLimits::MAX_T1D;
@@ -182,15 +182,15 @@ template<int DIM, int T_D1D, int T_Q1D>
 VectorMassIntegrator::VectorMassAddMultPAType
 VectorMassIntegrator::VectorMassAddMultPA::Kernel()
 {
-   if (DIM == 2)
+   if constexpr (DIM == 2)
   {
      return internal::SmemPAVectorMassApply2D<T_D1D,T_Q1D>;
   }
-   else if (DIM == 3)
+   else if constexpr (DIM == 3)
   {
      return internal::SmemPAVectorMassApply3D<T_D1D, T_Q1D>;
   }
-   else { MFEM_ABORT("Unsupported kernel"); }
+   MFEM_ABORT("Unsupported kernel");
 }

 inline VectorMassIntegrator::VectorMassAddMultPAType
@@ -301,18 +301,14 @@ template <int DIM, int T_D1D, int T_Q1D>
 DomainLFIntegrator::AssembleKernelType
 DomainLFIntegrator::AssembleKernels::Kernel()
 {
-   switch (DIM)
-   {
-      case 1:
-         return DLFEvalAssemble1D<T_D1D, T_Q1D>;
-      case 2:
-         return DLFEvalAssemble2D<T_D1D, T_Q1D>;
-      case 3:
-         return DLFEvalAssemble3D<T_D1D, T_Q1D>;
-   }
+   if constexpr (DIM == 1) { return DLFEvalAssemble1D<T_D1D, T_Q1D>; }
+   if constexpr (DIM == 2) { return DLFEvalAssemble2D<T_D1D, T_Q1D>; }
+   if constexpr (DIM == 3) { return DLFEvalAssemble3D<T_D1D, T_Q1D>; }
   MFEM_ABORT("");
 }
+
 /// \endcond DO_NOT_DOCUMENT

 } // namespace mfem
-#endif
+
+#endif // MFEM_LININTEG_DOMAIN_KERNELS_HPP
--- a/Show More
+++ b/Show More