Fix style error

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

@@ -45,11 +30,6 @@ 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
 =====================================
@@ -75,10 +75,12 @@ 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) and
+   #  * find_package(HDF5 REQUIRED) needed, in turn, by:
+   #    - find_package(AXOM REQUIRED)
   #  * find_package(PETSc REQUIRED)
   set(XSDK_ENABLE_C ON)
 endif()
@@ -476,6 +478,11 @@ 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,
@@ -622,7 +629,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 FMS CONDUIT Ginkgo GNUTLS GSLIB HDF5
+    SUNDIALS PETSC SLEPC MUMPS AXOM 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
@@ -452,6 +452,13 @@ 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
@@ -771,6 +778,14 @@ 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)
@@ -1054,6 +1069,7 @@ MFEM_USE_OCCA
 MFEM_USE_CEED
 MFEM_USE_RAJA
 MFEM_USE_UMPIRE
+MFEM_USE_SIDRE
 MFEM_USE_MOONOLITH
 MFEM_USE_CALIPER
 MFEM_USE_FMS
@@ -1117,6 +1133,7 @@ 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,6 +248,10 @@ 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,6 +46,7 @@ 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,6 +120,9 @@
 // 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

@@ -1,24 +0,0 @@
-# Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
-# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-# LICENSE and NOTICE for details. LLNL-CODE-806117.
-#
-# This file is part of the MFEM library. For more information and source code
-# availability visit https://mfem.org.
-#
-# MFEM is free software; you can redistribute it and/or modify it under the
-# terms of the BSD-3 license. We welcome feedback and contributions, see file
-# CONTRIBUTING.md for details.
-
-# 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,9 +13,6 @@
 #   - 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.
@@ -1,36 +0,0 @@
-# Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
-# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-# LICENSE and NOTICE for details. LLNL-CODE-806117.
-#
-# This file is part of the MFEM library. For more information and source code
-# availability visit https://mfem.org.
-#
-# MFEM is free software; you can redistribute it and/or modify it under the
-# terms of the BSD-3 license. We welcome feedback and contributions, see file
-# CONTRIBUTING.md for details.
-
-# 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,21 +21,6 @@ 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})
@@ -878,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_FMS MFEM_USE_CONDUIT MFEM_USE_PUMI
+      MFEM_USE_MPFR MFEM_USE_SIDRE 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
@@ -120,6 +120,9 @@
 // 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,6 +45,7 @@ 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,6 +48,7 @@ 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)
@@ -223,8 +224,17 @@ 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 (required by Tribol for the contact mini-app).")
+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(PUMI_DIR "${MFEM_DIR}/../pumi-2.1.0" CACHE STRING
    "Directory where PUMI is installed")
@@ -242,7 +252,6 @@ 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")
@@ -263,12 +272,8 @@ 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")
-# 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(Tribol_REQUIRED_PACKAGES "Axom/core/mint/slam/slic" CACHE STRING
+    "Additional packages required by Tribol")

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

@@ -162,6 +162,7 @@ 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,15 +249,6 @@ 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 =
@@ -495,6 +487,17 @@ 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
@@ -576,13 +579,7 @@ 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
@@ -592,12 +589,7 @@ 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) -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
+UMPIRE_LIB = -L$(UMPIRE_DIR)/lib -L$(UMPIRE_DIR)/lib64 -lumpire $(CAMP_LIB)
 ifdef FMT_DIR
   UMPIRE_OPT += -I$(FMT_DIR)/include
   UMPIRE_LIB += -L$(FMT_DIR)/lib -L$(FMT_DIR)/lib64 -lfmt
@@ -629,63 +621,8 @@ 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 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
+TRIBOL_LIB = -L$(TRIBOL_DIR)/lib -ltribol -lredecomp -L$(AXOM_DIR)/lib -laxom_mint\
+   -laxom_slam -laxom_slic -laxom_core

 # Enzyme configuration
 ENZYME_DIR = @MFEM_DIR@/../enzyme
@@ -5,9 +5,9 @@
 // Sample runs:
 //    mpirun -np 4 ex12p -m ../data/beam-tri.mesh
 //    mpirun -np 4 ex12p -m ../data/beam-quad.mesh
-//    mpirun -np 4 ex12p -m ../data/beam-tet.mesh -s 464 -n 10 -o 2 -elast
+//    mpirun -np 4 ex12p -m ../data/beam-tet.mesh -s 462 -n 10 -o 2 -elast
 //    mpirun -np 4 ex12p -m ../data/beam-hex.mesh -s 3878
-//    mpirun -np 4 ex12p -m ../data/beam-wedge.mesh -s 82
+//    mpirun -np 4 ex12p -m ../data/beam-wedge.mesh -s 81
 //    mpirun -np 4 ex12p -m ../data/beam-tri.mesh -s 3877 -o 2 -sys
 //    mpirun -np 4 ex12p -m ../data/beam-quad.mesh -s 4544 -n 6 -o 3 -elast
 //    mpirun -np 4 ex12p -m ../data/beam-quad-nurbs.mesh
@@ -302,21 +302,15 @@ int main(int argc, char *argv[])
      char vishost[] = "localhost";
      int  visport   = 19916;
      socketstream sol_sock_r(vishost, visport);
+      socketstream sol_sock_i(vishost, visport);
      sol_sock_r << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
      sol_sock_r.precision(8);
+      sol_sock_i.precision(8);
      sol_sock_r << "solution\n" << *pmesh << u_exact->real()
                 << "window_title 'Exact: Real Part'" << flush;
-      // Make sure all ranks have sent their real solution before initiating
-      // another set of GLVis connections (one from each rank):
-      MPI_Barrier(pmesh->GetComm());
-      socketstream sol_sock_i(vishost, visport);
-      sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
-      sol_sock_i.precision(8);
      sol_sock_i << "solution\n" << *pmesh << u_exact->imag()
                 << "window_title 'Exact: Imaginary Part'" << flush;
-      // Make sure all ranks have sent their imaginary solution before initiating
-      // another set of GLVis connections (one from each rank):
-      MPI_Barrier(pmesh->GetComm());
   }

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

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

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

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

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

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

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

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

   // Create data collection for solution output: either VisItDataCollection for
-   // ascii data files, or ConduitDataCollection for binary data files.
+   // ascii data files, or SidreDataCollection for binary data files.
   DataCollection *dc = NULL;
   if (visit)
   {
      if (binary)
      {
-#ifdef MFEM_USE_CONDUIT
-         auto conduit_dc = new ConduitDataCollection("Example9-Parallel", pmesh);
-         conduit_dc->SetProtocol("hdf5");
-         dc = conduit_dc;
+#ifdef MFEM_USE_SIDRE
+         dc = new SidreDataCollection("Example9-Parallel", pmesh);
 #else
-         MFEM_ABORT("Must build with MFEM_USE_CONDUIT=YES for binary output.");
+         MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
 #endif
      }
      else
@@ -213,6 +213,7 @@ set(HDRS
  dfem/qfunction_apply.hpp
  dfem/qfunction_transform.hpp
  dfem/tuple.hpp
+  dfem/univarsolvers.hpp
  dfem/util.hpp
  eltrans.hpp
  estimators.hpp
@@ -312,6 +313,11 @@ 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,8 +729,7 @@ void BilinearForm::Assemble(int skip_zeros)
         tr = mesh -> GetBdrFaceTransformations (i);
         if (tr != NULL)
         {
-            mfem::DofTransformation doftrans;
-            fes -> GetElementVDofs (tr -> Elem1No, vdofs, doftrans);
+            fes -> GetElementVDofs (tr -> Elem1No, vdofs);
            fe1 = fes -> GetFE (tr -> Elem1No);
            // The fe2 object is really a dummy and not used on the boundaries,
            // but we can't dereference a NULL pointer, and we don't want to
@@ -744,7 +743,6 @@ void BilinearForm::Assemble(int skip_zeros)

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

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


 /** Integrator for $(-Q u, \nabla v)$ for Nedelec ($u$) and $H^1$ ($v$) elements.
-    This is equivalent to a weak divergence of the $H(curl)$ basis functions. */
+    This is equivalent to a weak divergence of the $H(curl$ basis functions. */
 class VectorFEWeakDivergenceIntegrator: public BilinearFormIntegrator
 {
 protected:
@@ -1453,6 +1453,8 @@ 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 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.
+    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.

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

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

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

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

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

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

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

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

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

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

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

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

   using FiniteElement::CalcVShape;

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

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

 #ifdef MFEM_USE_MPI
 #include "pfespace.hpp"
@@ -62,6 +63,10 @@
 #include "pnonlinearform.hpp"
 #endif

+#ifdef MFEM_USE_SIDRE
+#include "sidredatacollection.hpp"
+#endif
+
 #ifdef MFEM_USE_CONDUIT
 #include "conduitdatacollection.hpp"
 #endif
@@ -3934,16 +3934,6 @@ const FiniteElement *FiniteElementSpace::GetBE(int i) const
   return BE;
 }

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

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

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

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

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

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

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

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

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

@@ -3065,7 +3064,6 @@ void GridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff,

 void GridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff)
 {
-   MFEM_VERIFY(VectorDim() == coeff.GetVDim(), "coeff vdim != VectorDim()");
   Array<int> dof_attr;
   ProjectDiscCoefficient(coeff, dof_attr);
 }
@@ -3075,10 +3073,6 @@ void GridFunction::ProjectDiscCoefficient(Coefficient &coeff, AvgType type)
   // Harmonic  (x1 ... xn) = [ (1/x1 + ... + 1/xn) / n ]^-1.
   // Arithmetic(x1 ... xn) = (x1 + ... + xn) / n.

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

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

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

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

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

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

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

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

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

-}
+}
@@ -150,13 +150,11 @@ public:

   FiniteElementCollection *OwnFEC() { return fec_owned; }

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

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

   /// Read only access to the (optional) internal true-dof Vector.
   const Vector &GetTrueVector() const
@@ -534,9 +532,6 @@ public:
   std::unique_ptr<GridFunction> ProlongateToMaxOrder() const;

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

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

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

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


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

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

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

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

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

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

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

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

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

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

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

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

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

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

   if (delta_c == NULL)
@@ -719,7 +717,6 @@ void ParGridFunction::ProjectCoefficientElementL2(VectorCoefficient &vcoeff)

 void ParGridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff)
 {
-   MFEM_VERIFY(VectorDim() == coeff.GetVDim(), "coeff vdim != VectorDim()");
   // local maximal element attribute for each dof
   Array<int> ldof_attr;

@@ -764,9 +761,6 @@ void ParGridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff)

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

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

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

-void ParGridFunction::ProjectBdrCoefficient(VectorCoefficient &vcoeff,
-                                            const Array<int> &attr)
-{
-   ProjectBdrCoefficient(NULL, &vcoeff, attr);
-}
-
 void ParGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient &vcoeff,
                                                   const Array<int> &bdr_attr)
 {
@@ -280,7 +280,8 @@ public:
   using GridFunction::ProjectBdrCoefficient;

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

   void ProjectBdrCoefficient(Coefficient *coeff[],
                              const Array<int> &attr) override
@@ -0,0 +1,539 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_SIDREDATACOLLECTION
+#define MFEM_SIDREDATACOLLECTION
+
+#include "../config/config.hpp"
+
+#ifdef MFEM_USE_SIDRE
+
+#include "datacollection.hpp"
+
+// Ignore warnings from the axom/sidre header (GCC + Clang versions)
+#ifdef MFEM_HAVE_GCC_PRAGMA_DIAGNOSTIC
+# pragma GCC diagnostic push
+# if defined(__clang__)
+#  pragma GCC diagnostic ignored "-Wextra-semi"
+# else // real GCC?
+#  pragma GCC diagnostic ignored "-Wpedantic"
+# endif
+#endif
+#include <axom/sidre.hpp>
+#ifdef MFEM_HAVE_GCC_PRAGMA_DIAGNOSTIC
+# pragma GCC diagnostic pop
+#endif
+
+namespace mfem
+{
+
+/** @brief Data collection with Sidre routines following the Conduit mesh
+    blueprint specification. */
+/** SidreDataCollection provides an HDF5-based file format for visualization or
+    restart capability.  This functionality is aimed primarily at customers of
+    LLNL's axom project that run problems at extreme scales.
+
+    For more information, see:
+    - Sidre component of LLNL's axom project (to be open-sourced), http://goo.gl/cZyJdn
+    - LLNL conduit/blueprint library, https://github.com/LLNL/conduit
+    - HDF5 library, https://support.hdfgroup.org/HDF5
+
+    The layout created in the Sidre DataStore is: (`"──"` denote groups,
+    `"─•"` denote views, `"─>"` denote links, i.e. shallow-copy view)
+
+        <root>
+         ├── <collection-name>_global          (global group)
+         │    └── blueprint_index
+         │         └── <collection-name>       (bp_index group)
+         │              ├── state
+         │              │    ├─• cycle
+         │              │    ├─• time
+         │              │    └─• number_of_domains = <mesh-mpi-comm-size>
+         │              ├── coordsets
+         │              │    └── coords
+         │              │         ├─• path = "<bp-path>/coordsets/coords"
+         │              │         ├─• type ─> <bp-grp>/coordsets/coords/type = "explicit"
+         │              │         └─• coord_system = "x"|"xy"|"xyz"
+         │              ├── topologies
+         │              │    ├── mesh
+         │              │    │    ├─• path = "<bp-path>/topologies/mesh"
+         │              │    │    ├─• type              ─> <bp-grp>/topologies/mesh/type = "unstructured"
+         │              │    │    ├─• coordset          ─> <bp-grp>/topologies/mesh/coordset = "coords"
+         │              │    │    ├─• grid_function     ─> <bp-grp>/topologies/mesh/grid_function = "<nodes-field-name>"
+         │              │    │    └─• boundary_topology ─> <bp-grp>/topologies/mesh/boundary_topology = "boundary"
+         │              │    └── boundary
+         │              │         ├─• path = "<bp-path>/topologies/mesh"
+         │              │         ├─• type     ─> <bp-grp>/topologies/boundary/type = "unstructured"
+         │              │         └─• coordset ─> <bp-grp>/topologies/boundary/coordset = "coords"
+         │              └── fields
+         │                   ├── mesh_material_attribute
+         │                   │    ├─• path = "<bp-path>/fields/mesh_material_attribute"
+         │                   │    ├─• association ─> <bp-grp>/fields/mesh_material_attribute/association = "element"
+         │                   │    ├─• topology    ─> <bp-grp>/fields/mesh_material_attribute/topology = "mesh"
+         │                   │    └─• number_of_components = 1
+         │                   ├── boundary_material_attribute
+         │                   │    ├─• path = "<bp-path>/fields/boundary_material_attribute"
+         │                   │    ├─• association ─> <bp-grp>/fields/boundary_material_attribute/association = "element"
+         │                   │    ├─• topology    ─> <bp-grp>/fields/boundary_material_attribute/topology = "boundary"
+         │                   │    └─• number_of_components = 1
+         │                   ├── grid-function-1
+         │                   │    ├─• path = "<bp-path>/fields/grid-function-1"
+         │                   │    ├─• basis    ─> <bp-grp>/fields/grid-function-1/basis = "<fe-coll-name>"
+         │                   │    ├─• topology ─> <bp-grp>/fields/grid-function-1/topology = "mesh"
+         │                   │    └─• number_of_components = gf1->VectorDim()
+         │                   ├── grid-function-2
+         │                   │    ├─• path = "<bp-path>/fields/grid-function-2"
+         │                   │    ├─• basis    ─> <bp-grp>/fields/grid-function-2/basis = "<fe-coll-name>"
+         │                   │    ├─• topology ─> <bp-grp>/fields/grid-function-2/topology = "mesh"
+         │                   │    └─• number_of_components = gf2->VectorDim()
+         │                   ├── ...
+         │                  ...
+         └── <collection-name>                 (domain group)
+              ├── blueprint                    (blueprint group)
+              │    ├── state
+              │    │    ├─• cycle
+              │    │    ├─• time
+              │    │    ├─• domain = <mesh-mpi-rank>
+              │    │    └─• time_step
+              │    ├── coordsets
+              │    │    └── coords
+              │    │         ├─• type = "explicit"
+              │    │         └── values
+              │    │              ├─• x = view in <vertex-coords-buffer>/<ext-double-data>
+              │    │              ├─• y = view in <vertex-coords-buffer>/<ext-double-data>
+              │    │              └─• z = view in <vertex-coords-buffer>/<ext-double-data>
+              │    ├── topologies
+              │    │    ├── mesh
+              │    │    │    ├─• type = "unstructured"
+              │    │    │    ├── elements
+              │    │    │    │    ├─• shape = "points"|"lines"|...
+              │    │    │    │    └─• connectivity = <vert-idx-array>
+              │    │    │    ├─• coordset = "coords"
+              │    │    │    ├─• grid_function = "<nodes-field-name>"
+              │    │    │    └─• boundary_topology = "boundary"
+              │    │    └── boundary
+              │    │         ├─• type = "unstructured"
+              │    │         ├── elements
+              │    │         │    ├─• shape = "points"|"lines"|...
+              │    │         │    └─• connectivity = <vert-idx-array>
+              │    │         └─• coordset = "coords"
+              │    └── fields
+              │         ├── mesh_material_attribute
+              │         │    ├─• association = "element"
+              │         │    ├─• topology = "mesh"
+              │         │    └─• values = <attr-array>
+              │         ├── boundary_material_attribute
+              │         │    ├─• association = "element"
+              │         │    ├─• topology = "boundary"
+              │         │    └─• values = <attr-array>
+              │         ├── grid-function-1   (name can include path)
+              │         │    ├─• basis = "<fe-coll-name>"
+              │         │    ├─• topology = "mesh"
+              │         │    └─• values = <ext-double-array>/<named-buffer> (vdim == 1)
+              │         ├── grid-function-2   (name can include path)
+              │         │    ├─• basis = "<fe-coll-name>"
+              │         │    ├─• topology = "mesh"
+              │         │    └── values   (vdim > 1)
+              │         │         ├─• x0 = view into <ext-double-array>/<named-buffer>
+              │         │         ├─• x1 = view into <ext-double-array>/<named-buffer>
+              │         │         └─• x2 = view into <ext-double-array>/<named-buffer>
+              │         ├── ...
+              │        ...
+              └── named_buffers                (named_buffers group)
+                   ├─• vertex_coords = <double-array>
+                   ├─• grid-function-1 = <double-array>
+                   ├─• grid-function-2 = <double-array>
+                  ...
+
+    @note blueprint_index is used both in serial and in parallel. In parallel,
+    only rank 0 will add entries to the blueprint index.
+
+    @note QuadratureFunction%s (q-fields) are not supported.
+
+    @note SidreDataCollection does not manage the FiniteElementSpace%s and
+    FiniteElementCollection%s associated with registered GridFunction%s.
+    Therefore, field registration is left to the user of SidreDataCollection and
+    there are no methods that automatically register GridFunction%s using just
+    the content of the Sidre DataStore. Such capabilities can be implemented in
+    a derived class, adding any desired object management routines.
+
+    @warning This class is still _experimental_, meaning that in future
+    releases, it may not be backward compatible, and the output files generated
+    by the current version may become unreadable.
+*/
+class SidreDataCollection : public DataCollection
+{
+public:
+   typedef NamedFieldsMap< Array<int> > AttributeFieldMap;
+   AttributeFieldMap attr_map;
+
+public:
+
+   /// Constructor that allocates and initializes a Sidre DataStore.
+   /**
+       @param[in] collection_name  Name of the collection used as a file name
+                                   when saving
+       @param[in] the_mesh         Mesh shared by all grid functions in the
+                                   collection (can be NULL)
+       @param[in] owns_mesh_data   Does the SidreDC own the mesh vertices?
+
+       With this constructor, the SidreDataCollection owns the allocated Sidre
+       DataStore.
+    */
+   SidreDataCollection(const std::string& collection_name,
+                       Mesh *the_mesh = NULL,
+                       bool owns_mesh_data = false);
+
+   /// Constructor that links to an external Sidre DataStore.
+   /** Specifically, the global and domain groups can be at arbitrary paths.
+
+       @param[in] collection_name  Name of the collection used as a file name
+                                   when saving
+       @param[in] bp_index_grp     Pointer to the blueprint index group in the
+                                   datastore, see the above schematic
+       @param[in] domain_grp       Pointer to the domain group in the datastore,
+                                   see the above schematic
+       @param[in] owns_mesh_data   Does the SidreDC own the mesh vertices?
+
+       With this constructor, the SidreDataCollection does not own the Sidre
+       DataStore.
+       @note No mesh or fields are read from the given Groups. The mesh has
+       to be set with SetMesh() and fields registered with RegisterField().
+    */
+   SidreDataCollection(const std::string& collection_name,
+                       axom::sidre::Group * bp_index_grp,
+                       axom::sidre::Group * domain_grp,
+                       bool owns_mesh_data = false);
+
+#ifdef MFEM_USE_MPI
+   /// Associate an MPI communicator with the collection.
+   /** If no mesh was associated with the collection, this method should be
+       called before using any of the Load() methods to read parallel data. */
+   void SetComm(MPI_Comm comm);
+#endif
+
+   /// Register a GridFunction in the Sidre DataStore.
+   /** This method is a shortcut for the call
+       `RegisterField(field_name, gf, field_name, 0)`.
+    */
+   virtual void RegisterField(const std::string &field_name, GridFunction *gf)
+   {
+      RegisterField(field_name, gf, field_name, 0);
+   }
+
+   /// Register a GridFunction in the Sidre DataStore.
+   /** The registration procedure is as follows:
+       - if (@a gf's data is NULL), allocate named buffer with the name
+         @a buffer_name with size _offset + gf->FESpace()->GetVSize()_ and use
+         its data (plus the given @a offset) to set @a gf's data;
+       - else, if (DataStore has a named buffer @a buffer_name), replace @a gf's
+         data array with that named buffer plus the given @a offset;
+       - else, use @a gf's data as external data associated with @a field_name
+         in the DataStore;
+       - register @a field_name in #field_map.
+
+       Both the @a field_name and @a buffer_name can contain a path prefix.
+       @note If @a field_name or @a buffer_name is empty, the method does
+       nothing.
+       @note If the GridFunction pointer @a gf or it's FiniteElementSpace
+       pointer are NULL, the method does nothing.
+    */
+   void RegisterField(const std::string &field_name, GridFunction *gf,
+                      const std::string &buffer_name,
+                      axom::sidre::IndexType offset);
+
+   /// Registers an attribute field in the Sidre DataStore
+   /** The registration process is similar to that of RegisterField()
+       The attribute field is associated with the elements of the mesh
+       when @a is_bdry is false, and with the boundary elements, when
+       @a is_bdry is true.
+       @sa RegisterField()  */
+   void RegisterAttributeField(const std::string& name, bool is_bdry);
+   void DeregisterAttributeField(const std::string& name);
+
+   /** Returns a pointer to the attribute field associated with
+       @a field_name, or NULL when there is no associated field */
+   Array<int>* GetAttributeField(const std::string& field_name) const
+   { return attr_map.Get(field_name); }
+
+   /** Checks if there is an attribute field associated with @a field_name */
+   bool HasAttributeField(const std::string& field_name) const
+   { return attr_map.Has(field_name); }
+
+   /** Checks if any rank in the mesh has boundary elements */
+   bool HasBoundaryMesh() const;
+
+   /// Set the name of the mesh nodes field.
+   /** This name will be used by SetMesh() to register the mesh nodes, if not
+       already registered. Also, this method should be called if the mesh nodes
+       GridFunction was or will be registered directly by the user. The default
+       value for the name is "mesh_nodes". */
+   void SetMeshNodesName(const std::string &nodes_name)
+   {
+      if (!nodes_name.empty()) { m_meshNodesGFName = nodes_name; }
+   }
+
+   /// De-register @a field_name from the SidreDataCollection.
+   /** The field is removed from the #field_map and the DataStore, including
+       deleting it from the named_buffers group, if allocated. */
+   virtual void DeregisterField(const std::string& field_name);
+
+   /// Delete all owned data.
+   virtual ~SidreDataCollection();
+
+   /// Set/change the mesh associated with the collection
+   /** Uses the field name "mesh_nodes" or the value set by SetMeshNodesName()
+       to register the mesh nodes GridFunction, if the mesh uses nodes. */
+   virtual void SetMesh(Mesh *new_mesh);
+
+#ifdef MFEM_USE_MPI
+   /// Set/change the mesh associated with the collection
+   /** Uses the field name "mesh_nodes" or the value set by SetMeshNodesName()
+       to register the mesh nodes GridFunction, if the mesh uses nodes. */
+   virtual void SetMesh(MPI_Comm comm, Mesh *new_mesh);
+#endif
+
+   /// Reset the domain and global datastore group pointers.
+   /** These are set in the constructor, but if a host code changes the
+       datastore contents ( such as wiping out the datastore and loading in new
+       contents from a file, i.e. a restart ) these pointers will need to be
+       reset to valid groups in the datastore.
+       @sa Load(const std::string &path, const std::string &protocol).
+    */
+   void SetGroupPointers(axom::sidre::Group * global_grp,
+                         axom::sidre::Group * domain_grp);
+
+   axom::sidre::Group * GetBPGroup() { return m_bp_grp; }
+   axom::sidre::Group * GetBPIndexGroup() { return m_bp_index_grp; }
+
+   /// Prepare the DataStore for writing
+   virtual void PrepareToSave();
+
+   /// Save the collection to file.
+   /** This method calls `Save(collection_name, "sidre_hdf5")`. */
+   virtual void Save();
+
+   /// Save the collection to @a filename.
+   /** The collection path prefix is prepended to the @a filename and the
+       current cycle is appended, if cycle >= 0. */
+   void Save(const std::string& filename, const std::string& protocol);
+
+   /// Load the Sidre DataStore from file.
+   /** No mesh or fields are read from the loaded DataStore.
+
+       If the data collection created the datastore, it knows the layout of
+       where the domain and global groups are, and can restore them after the
+       Load().
+
+       If, however, the data collection does not own the datastore (e.g. it did
+       not create the datastore), the host code must reset these pointers after
+       the load operation, using SetGroupPointers(), and also reset the state
+       variables, using UpdateStateFromDS().
+    */
+   void Load(const std::string& path, const std::string& protocol);
+
+   /// Load SidreDataCollection from file.
+   /** The used file path is based on the current prefix path, collection name,
+       and the given @a cycle_. The protocol is "sidre_hdf5".
+       @sa Load(const std::string &path, const std::string &protocol).
+    */
+   virtual void Load(int cycle_ = 0)
+   {
+      SetCycle(cycle_);
+      Load(get_file_path(name), "sidre_hdf5");
+   }
+
+   /// Load external data after registering externally owned fields.
+   void LoadExternalData(const std::string& path);
+
+   /** @brief Updates the DataCollection's cycle, time, and time-step variables
+       with the values from the data store. */
+   void UpdateStateFromDS();
+
+   /** @brief Updates the data store's cycle, time, and time-step variables with
+       the values from the SidreDataCollection. */
+   void UpdateStateToDS();
+
+   /** @name Methods for named buffer access and manipulation. */
+   ///@{
+
+   /** @brief Get a pointer to the sidre::View holding the named buffer for
+       @a buffer_name. */
+   /** If such named buffer is not allocated, the method returns NULL.
+       @note To access the underlying pointer, use View::getData().
+       @note To query the size of the buffer, use View::getNumElements().
+    */
+   axom::sidre::View *
+   GetNamedBuffer(const std::string& buffer_name) const
+   {
+      return named_buffers_grp()->hasView(buffer_name)
+             ? named_buffers_grp()->getView(buffer_name)
+             : NULL;
+   }
+
+   /// Return newly allocated or existing named buffer for @a buffer_name.
+   /** The buffer is stored in the named_buffers group. If the currently
+       allocated buffer size is smaller than @a sz, then the buffer is
+       reallocated with size @a sz, destroying its contents.
+       @note To access the underlying pointer, use View::getData().
+    */
+   axom::sidre::View *
+   AllocNamedBuffer(const std::string& buffer_name,
+                    axom::sidre::IndexType sz,
+                    axom::sidre::TypeID type =
+                       axom::sidre::DOUBLE_ID);
+
+   /// Deallocate the named buffer @a buffer_name.
+   void FreeNamedBuffer(const std::string& buffer_name)
+   { named_buffers_grp()->destroyViewAndData(buffer_name); }
+
+   ///@}
+
+private:
+   // Used if the Sidre data collection is providing the datastore itself.
+   const bool m_owns_datastore;
+
+   // TODO - Need to evaluate if this bool member can be combined with own_data
+   // in parent data collection class. m_owns_mesh_data indicates whether the
+   // Sidre dc owns the mesh element data and node positions gf. The DC base
+   // class own_data indicates if the dc owns the mesh object pointer itself and
+   // GF objects. Can we use one flag and just have DC own all objects vs none?
+   const bool m_owns_mesh_data;
+
+   // Name to be used for registering the mesh nodes in the SidreDataCollection.
+   // This name is used by SetMesh() and can be overwritten by the method
+   // SetMeshNodesName().
+   // Default value: "mesh_nodes".
+   std::string m_meshNodesGFName;
+
+   // If the data collection owns the datastore, it will store a pointer to it.
+   // Otherwise, this pointer is NULL.
+   axom::sidre::DataStore * m_datastore_ptr;
+
+protected:
+   axom::sidre::Group *named_buffers_grp() const;
+
+   axom::sidre::View *
+   alloc_view(axom::sidre::Group *grp,
+              const std::string &view_name);
+
+   axom::sidre::View *
+   alloc_view(axom::sidre::Group *grp,
+              const std::string &view_name,
+              const axom::sidre::DataType &dtype);
+
+   axom::sidre::Group *
+   alloc_group(axom::sidre::Group *grp,
+               const std::string &group_name);
+
+   // return the filename based on prefix_path, collection name and cycle.
+   std::string get_file_path(const std::string &filename) const;
+
+private:
+   // If the data collection does not own the datastore, it will need pointers
+   // to the blueprint and blueprint index group to use.
+   axom::sidre::Group * m_bp_grp;
+   axom::sidre::Group * m_bp_index_grp;
+
+   // This is stored for convenience.
+   axom::sidre::Group * m_named_bufs_grp;
+
+   // Private helper functions
+
+   void RegisterFieldInBPIndex(const std::string& field_name,
+                               GridFunction *gf);
+   void DeregisterFieldInBPIndex(const std::string & field_name);
+
+   void RegisterAttributeFieldInBPIndex(const std::string& attr_name);
+   void DeregisterAttributeFieldInBPIndex(const std::string& attr_name);
+
+   /** @brief Return a string with the conduit blueprint name for the given
+       Element::Type. */
+   std::string getElementName( Element::Type elementEnum );
+
+   /**
+    * \brief A private helper function to set up the views associated with the
+       data of a scalar valued grid function in the blueprint style.
+    * \pre gf is not null
+    * \note This function is expected to be called by RegisterField()
+    * \note Handles cases where hierarchy is already set up,
+    *      where the data was allocated by this data collection
+    *      and where the grid function data is external to Sidre
+    */
+   void addScalarBasedGridFunction(const std::string& field_name,
+                                   GridFunction* gf,
+                                   const std::string &buffer_name,
+                                   axom::sidre::IndexType offset);
+
+   /**
+    * \brief A private helper function to set up the views associated with the
+       data of a vector valued grid function in the blueprint style.
+    * \pre gf is not null
+    * \note This function is expected to be called by RegisterField()
+    * \note Handles cases where hierarchy is already set up,
+    *      where the data was allocated by this data collection
+    *      and where the grid function data is external to Sidre
+    */
+   void addVectorBasedGridFunction(const std::string& field_name,
+                                   GridFunction* gf,
+                                   const std::string &buffer_name,
+                                   axom::sidre::IndexType offset);
+
+   /** @brief A private helper function to set up the Views associated with
+       attribute field named @a field_name */
+   void addIntegerAttributeField(const std::string& field_name, bool is_bdry);
+
+   /// Sets up the four main mesh blueprint groups.
+   /**
+    * \param hasBP Indicates whether the blueprint has already been set up.
+    */
+   void createMeshBlueprintStubs(bool hasBP);
+
+   /// Sets up the mesh blueprint 'state' group.
+   /**
+    * \param hasBP Indicates whether the blueprint has already been set up.
+    */
+   void createMeshBlueprintState(bool hasBP);
+
+   /// Sets up the mesh blueprint 'coordsets' group.
+   /**
+    * \param hasBP Indicates whether the blueprint has already been set up.
+    */
+   void createMeshBlueprintCoordset(bool hasBP);
+
+   /// Sets up the mesh blueprint 'topologies' group.
+   /**
+    * This method is called from SetMesh().
+    * \param hasBP Indicates whether the blueprint has already been set up.
+    * \param mesh_name The name of the topology.
+    * \note Valid values for @a mesh_name are "mesh" and "boundary" and the
+            former has to be created with this method before the latter.
+    */
+   void createMeshBlueprintTopologies(bool hasBP, const std::string& mesh_name);
+
+#ifdef MFEM_USE_MPI
+   /// Sets up the mesh blueprint 'adjacencies' group.
+   /**
+    * \param hasBP Indicates whether the blueprint has already been set up.
+    * \note Only valid when using parallel meshes
+    */
+   void createMeshBlueprintAdjacencies(bool hasBP);
+#endif
+
+   /// Verifies that the contents of the mesh blueprint data is valid.
+   void verifyMeshBlueprint();
+};
+
+} // end namespace mfem
+
+#endif
+
+#endif
@@ -278,18 +278,9 @@ void ArraysByName<T>::Load(std::istream &in)
         q1 = ArrayLine.find(' ');
         ArrayName = ArrayLine.substr(0,q1-1);
      }
-      if (q1+2 < ArrayLine.size())
-      {
-         // Read the remainder of the line which contains the array data
-         std::istringstream ArrayDataStream(ArrayLine.substr(q1+2,
-                                                             ArrayLine.size()));
-         data[ArrayName].Load(ArrayDataStream, 0);
-      }
-      else
-      {
-         // Read the array data starting on the next line
-         data[ArrayName].Load(in, 0);
-      }
+
+      // Ignore the remainder of the line which may contain explanatory comments
+      data[ArrayName].Load(in, 0);
   }

 }
@@ -726,16 +726,16 @@ std::string Device::GetUUID(const int device_id)
   MFEM_GPU_CHECK(cudaGetDeviceProperties(&prop, device_id));
   for (int i = 0; i < 16; ++i)
   {
-      const unsigned b = static_cast<unsigned char>(prop.uuid.bytes[i]);
-      res << std::setfill('0') << std::setw(2) << std::hex << b;
+      res << std::setfill('0') << std::setw(2) << std::hex
+          << static_cast<unsigned>(prop.uuid.bytes[i]);
   }
 #elif defined(MFEM_USE_HIP)
   hipUUID uuid;
   MFEM_GPU_CHECK(hipDeviceGetUuid(&uuid, device_id));
   for (int i = 0; i < 16; ++i)
   {
-      const unsigned b = static_cast<unsigned char>(uuid.bytes[i]);
-      res << std::setfill('0') << std::setw(2) << std::hex << b;
+      res << std::setfill('0') << std::setw(2) << std::hex
+          << static_cast<unsigned>(uuid.bytes[i]);
   }
 #endif
   return res.str();
@@ -160,6 +160,9 @@ const char *GetConfigStr()
 #ifdef MFEM_USE_RAJA
      "MFEM_USE_RAJA\n"
 #endif
+#ifdef MFEM_USE_SIDRE
+      "MFEM_USE_SIDRE\n"
+#endif
 #ifdef MFEM_USE_SIMD
      "MFEM_USE_SIMD\n"
 #endif
@@ -4156,31 +4156,20 @@ void PetscNonlinearSolver::SetUpdate(void (*update)(Operator *,int,
 void PetscNonlinearSolver::Mult(const Vector &b, Vector &x) const
 {
   SNES snes = (SNES)obj;
-   MPI_Comm comm = PetscObjectComm(obj);

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

   Customize();

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

-   // Solve the system. Pass nullptr for b when empty (PETSc treats it as zero RHS).
-   ierr = SNESSolve(snes, b_nonempty ? B->x : nullptr, X->x); PCHKERRQ(snes, ierr);
+   // Solve the system.
+   ierr = SNESSolve(snes, B->x, X->x); PCHKERRQ(snes, ierr);
   X->ResetMemory();
   if (b_nonempty) { B->ResetMemory(); }
 }
@@ -299,7 +299,7 @@ ifeq ($(MFEM_USE_LEGACY_OPENMP),YES)
 endif

 # List of MFEM dependencies, that require the *_LIB variable to be non-empty
-MFEM_REQ_LIB_DEPS = SUPERLU MUMPS METIS FMS CONDUIT LAPACK SUNDIALS\
+MFEM_REQ_LIB_DEPS = SUPERLU MUMPS METIS FMS CONDUIT SIDRE LAPACK SUNDIALS\
 SUITESPARSE STRUMPACK GINKGO GNUTLS HDF5 NETCDF SLEPC PETSC MPFR PUMI HIOP\
 GSLIB OCCA CEED RAJA UMPIRE MKL_CPARDISO MKL_PARDISO AMGX MAGMA CALIPER PARELAG\
 TRIBOL BENCHMARK MOONOLITH ALGOIM
@@ -365,7 +365,7 @@ MFEM_DEFINES = MFEM_VERSION MFEM_VERSION_STRING MFEM_GIT_STRING MFEM_USE_MPI\
 MFEM_USE_LEGACY_OPENMP MFEM_USE_MEMALLOC MFEM_TIMER_TYPE MFEM_USE_SUNDIALS\
 MFEM_USE_SUITESPARSE MFEM_USE_GINKGO MFEM_USE_SUPERLU MFEM_USE_SUPERLU5\
 MFEM_USE_STRUMPACK MFEM_USE_GNUTLS MFEM_USE_HDF5 MFEM_USE_NETCDF MFEM_USE_PETSC\
- MFEM_USE_SLEPC MFEM_USE_MPFR MFEM_USE_FMS MFEM_USE_CONDUIT\
+ MFEM_USE_SLEPC MFEM_USE_MPFR MFEM_USE_SIDRE MFEM_USE_FMS MFEM_USE_CONDUIT\
 MFEM_USE_PUMI MFEM_USE_HIOP MFEM_USE_GSLIB MFEM_USE_CUDA MFEM_USE_HIP\
 MFEM_USE_OCCA MFEM_USE_MOONOLITH MFEM_USE_CEED MFEM_USE_RAJA MFEM_USE_UMPIRE\
 MFEM_USE_SIMD MFEM_USE_ADIOS2 MFEM_USE_MKL_CPARDISO MFEM_USE_MKL_PARDISO MFEM_USE_AMGX\
@@ -746,6 +746,7 @@ status info:
 	$(info MFEM_USE_PETSC         = $(MFEM_USE_PETSC))
 	$(info MFEM_USE_SLEPC         = $(MFEM_USE_SLEPC))
 	$(info MFEM_USE_MPFR          = $(MFEM_USE_MPFR))
+	$(info MFEM_USE_SIDRE         = $(MFEM_USE_SIDRE))
 	$(info MFEM_USE_FMS           = $(MFEM_USE_FMS))
 	$(info MFEM_USE_CONDUIT       = $(MFEM_USE_CONDUIT))
 	$(info MFEM_USE_PUMI          = $(MFEM_USE_PUMI))
@@ -3206,22 +3206,10 @@ public:


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

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

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

-   if (MyRank == save_rank)
-   {
-      attribute_sets.Copy(serialmesh.attribute_sets);
-      bdr_attribute_sets.Copy(serialmesh.bdr_attribute_sets);
-   }
-
   MPI_Barrier(MyComm);
   return serialmesh;
 }
@@ -82,43 +82,6 @@ Although Tribol can be built automatically via **uberenv** and **Spack**,
 for this miniapp it is simpler to build **Axom** and **MFEM** manually and 
 point Tribol to them. The steps are as follows:

-### Using pre-built Tribol/Axom installs
-
-If you already have compatible installs of Tribol and Axom, point MFEM to the install prefixes.
-
- Hypre install prefix: `<path/to/hypre>`
- METIS install prefix: `<path/to/metis>`
- Axom install prefix: `<path/to/axom>`
- Tribol install prefix: `<path/to/tribol>`
-
-**MFEM make build (configure):**
-```bash
-make config MFEM_USE_MPI=YES MFEM_USE_METIS=YES MFEM_USE_TRIBOL=YES \
-  HYPRE_DIR=<path/to/hypre> METIS_DIR=<path/to/metis> \
-  AXOM_DIR=<path/to/axom> TRIBOL_DIR=<path/to/tribol> ADIAK_DIR=<path/to/adiak> CAMP_DIR=<path/to/camp> RAJA_DIR=<path/to/raja> \
-  UMPIRE_DIR=<path/to/umpire> FMT_DIR=<path/to/fmt> CALIPER_DIR=<path/to/caliper>
-```
-
-**MFEM CMake build (configure):**
-```bash
-cmake -S . -B <mfem-build-dir> -DMFEM_USE_MPI=YES -DMFEM_USE_METIS=YES -DMFEM_USE_TRIBOL=YES \
-  HYPRE_DIR=<path/to/hypre> METIS_DIR=<path/to/metis> \
-  AXOM_DIR=<path/to/axom> TRIBOL_DIR=<path/to/tribol> ADIAK_DIR=<path/to/adiak> CAMP_DIR=<path/to/camp> RAJA_DIR=<path/to/raja> \
-  UMPIRE_DIR=<path/to/umpire> FMT_DIR=<path/to/fmt> CALIPER_DIR=<path/to/caliper>
-```
-
-Note: RAJA/UMPIRE/CALIPER are optional for MFEM itself, but many Tribol builds
-enable them. If your Tribol install does not depend on them, you can omit the
-corresponding `*_DIR` entries above.
-
-Note: `FMT_DIR` only needs to be added for the make-based build (and only when
-the Umpire install uses `fmt`). If `FMT_DIR` is not set and a sibling `fmt-*`
-directory exists next to your `UMPIRE_DIR`, MFEM's make configuration will try
-to pick it up automatically.
-
-Note: when using pre-built Tribol/Axom, you typically need to use a compatible
-compiler/MPI wrapper (same C++ standard library ABI).
-
 ### Manual Build Steps

 1. Pull axom and tribol (starting from the mfem folder):
@@ -136,7 +99,7 @@ compiler/MPI wrapper (same C++ standard library ABI).
   TRIBOL_DIR = @MFEM_DIR@/../tribol-repo/tribol
   TRIBOL_OPT = -I$(TRIBOL_DIR)/include -I$(AXOM_DIR)/include
   TRIBOL_LIB = -L$(TRIBOL_DIR)/lib -ltribol -lredecomp -L$(AXOM_DIR)/lib \
-   -laxom_quest -laxom_mint -laxom_slam -laxom_slic -laxom_lumberjack -laxom_core
+   -laxom_mint -laxom_slam -laxom_slic -laxom_core
   ```
 3. [**Axom:**](https://github.com/LLNL/axom.git) Starting from the MFEM root
   directory (we assume this directory is named mfem):
@@ -18,6 +18,7 @@
 //
 // Currently supported data collection type options:
 //    visit:                VisItDataCollection (default)
+//    sidre or sidre_hdf5:  SidreDataCollection
 //    json:                 ConduitDataCollection w/ protocol json
 //    conduit_json:         ConduitDataCollection w/ protocol conduit_json
 //    conduit_bin:          ConduitDataCollection w/ protocol conduit_bin
@@ -51,6 +52,14 @@ DataCollection *create_data_collection(const std::string &dc_name,
      dc = new VisItDataCollection(MPI_COMM_WORLD, dc_name);
 #else
      dc = new VisItDataCollection(dc_name);
+#endif
+   }
+   else if ( dc_type == "sidre" || dc_type == "sidre_hdf5")
+   {
+#ifdef MFEM_USE_SIDRE
+      dc = new SidreDataCollection(dc_name);
+#else
+      MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for sidre support.");
 #endif
   }
   else if ( dc_type == "json" ||
@@ -131,6 +140,7 @@ int main(int argc, char *argv[])
   args.AddOption(&src_coll_type, "-st", "--source-type",
                  "Set the source data collection type. Options:\n"
                  "\t   visit:                VisItDataCollection (default)\n"
+                  "\t   sidre or sidre_hdf5:  SidreDataCollection\n"
                  "\t   json:                 ConduitDataCollection w/ protocol json\n"
                  "\t   conduit_json:         ConduitDataCollection w/ protocol conduit_json\n"
                  "\t   conduit_bin:          ConduitDataCollection w/ protocol conduit_bin\n"
@@ -142,6 +152,7 @@ int main(int argc, char *argv[])
   args.AddOption(&out_coll_type, "-ot", "--output-type",
                  "Set the output data collection type. Options:\n"
                  "\t   visit:                VisItDataCollection (default)\n"
+                  "\t   sidre or sidre_hdf5:  SidreDataCollection\n"
                  "\t   json:                 ConduitDataCollection w/ protocol json\n"
                  "\t   conduit_json:         ConduitDataCollection w/ protocol conduit_json\n"
                  "\t   conduit_bin:          ConduitDataCollection w/ protocol conduit_bin\n"
@@ -39,6 +39,7 @@ set(UNIT_TESTS_SRCS
  dfem/test_divergence.cpp
  dfem/test_lvector_interface.cpp
  dfem/test_mass.cpp
+  dfem/test_univarsolver.cpp
  general/test_array.cpp
  general/test_scan.cpp
  general/test_arrays_by_name.cpp
@@ -71,7 +72,6 @@ set(UNIT_TESTS_SRCS
  linalg/test_ode2.cpp
  linalg/test_operator.cpp
  linalg/test_particlevector.cpp
-  linalg/test_petsc_nonlinear.cpp
  linalg/test_sparsesmoothers.cpp
  linalg/test_vector.cpp
  mesh/mesh_test_utils.cpp
@@ -0,0 +1,359 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include <algorithm>
+#include <cmath>
+
+#include "mfem.hpp"
+#include "unit_tests.hpp"
+
+#ifdef MFEM_USE_ENZYME
+
+using mfem::real_t;
+using namespace mfem::future;
+
+MFEM_HOST_DEVICE inline real_t FlowResistance(real_t eqps, real_t sigma_y, real_t n, real_t ep_0)
+{
+    return sigma_y*(1.0 + std::pow((eqps)/ep_0, n));
+}
+
+using J2PlasticityParameters = tuple<real_t, real_t, real_t, real_t, real_t, real_t>;
+
+// Residual function that is solved in the plasticity model.
+// Made a free function to facilitate Enzyme differentiation.
+real_t J2PlasticityResidual(real_t delta_eqps, J2PlasticityParameters p)
+{
+    auto [eqps, q, G, sigma_y, n, ep_0] = p;
+    return q - 3.0*G*delta_eqps - FlowResistance(eqps + delta_eqps, sigma_y, n, ep_0);
+}
+
+struct J2Plasticity {
+  static constexpr int dim = 3;  ///< spatial dimension
+  static constexpr int N_INTERNAL_STATES = 10;
+  static constexpr real_t tol = 1e-10;  ///< relative tolerance on residual mag to judge convergence of return map
+
+  real_t E;        ///< Young's modulus
+  real_t nu;       ///< Poisson's ratio
+  real_t sigma_y;  ///< Yield strength
+  real_t n;        ///< Hardening index
+  real_t ep_0;     ///< Reference plastic strain
+  
+  /// @brief variables required to characterize the hysteresis response
+  struct InternalState {
+    tensor<real_t, dim, dim> plastic_strain;
+    real_t accumulated_plastic_strain;
+  };
+
+  /// Internal state variables in a flattened array for storing in a global field
+  using PackedInternalState = mfem::future::tensor<real_t, N_INTERNAL_STATES>;
+
+  // Unflatten internal state variables
+  MFEM_HOST_DEVICE static inline InternalState unpack_internal_state(
+      const mfem::future::tensor<real_t, N_INTERNAL_STATES>& packed_state)
+  {
+    auto plastic_strain =
+        mfem::future::make_tensor<dim, dim>([&packed_state](int i, int j) { return packed_state[dim * i + j]; });
+    real_t accumulated_plastic_strain = packed_state[N_INTERNAL_STATES - 1];
+    return {plastic_strain, accumulated_plastic_strain};
+  }
+
+  // Flatten internal state variables (for repacking into global field)
+  MFEM_HOST_DEVICE static inline PackedInternalState pack_internal_state(
+      const mfem::future::tensor<real_t, dim, dim>& plastic_strain, real_t accumulated_plastic_strain)
+  {
+    PackedInternalState packed_state{};
+    for (int i = 0, ij = 0; i < dim; i++) {
+      for (int j = 0; j < dim; j++, ij++) {
+        packed_state[ij] = plastic_strain[i][j];
+      }
+    }
+    packed_state[N_INTERNAL_STATES - 1] = accumulated_plastic_strain;
+    return packed_state;
+  }
+
+  // Compute the new stress and the internal state variables
+  MFEM_HOST_DEVICE inline tuple<tensor<real_t, dim, dim>, PackedInternalState>
+  update(tensor<real_t, dim, dim> dudxi,
+         PackedInternalState internal_state,
+         tensor<real_t, dim, dim> J,
+         real_t w) const
+  {
+    auto invJ = inv(J);
+    const auto dudX = dudxi * invJ;
+    auto I = IdentityMatrix<dim>();
+    const real_t K = E / (3.0 * (1.0 - 2.0 * nu));
+    const real_t G = 0.5 * E / (1.0 + nu);
+
+    auto [plastic_strain, accumulated_plastic_strain] = unpack_internal_state(internal_state);
+    
+    auto el_strain = sym(dudX) - plastic_strain;
+    auto p = K * tr(el_strain);
+    auto s = 2.0 * G * dev(el_strain);
+    auto q = std::sqrt(1.5) * norm(s);
+    real_t denom = q > 0.0? q : 1.0;
+    auto Np = 1.5 * s / denom;
+
+    if (q > FlowResistance(accumulated_plastic_strain, sigma_y, n, ep_0)) {
+        real_t lb = 0.0;
+        real_t ub = (q - FlowResistance(accumulated_plastic_strain, sigma_y, n, ep_0))/(3*G);
+        SolverSettings settings{1e-10*sigma_y, 1e-10, {lb, ub}};
+        // Use the differentiable univariate root finder.
+        // This has custom derivatives, so it's ok to differentiate this enclosing function.
+        real_t delta_eqps = SolveNewtonBisection<J2PlasticityResidual>(
+            0.5*(lb + ub), make_tuple(accumulated_plastic_strain, q, G, sigma_y, n, ep_0), settings);
+        accumulated_plastic_strain += delta_eqps;
+        plastic_strain += delta_eqps * Np;
+        s -= 2.0 * G * delta_eqps * Np;
+    }
+    auto Q_new = pack_internal_state(plastic_strain, accumulated_plastic_strain);
+    auto stress = s + p * I;
+    const real_t dV = det(J)*w;
+    return {stress*transpose(invJ)*dV, Q_new};
+  }
+  
+  MFEM_HOST_DEVICE inline tensor<real_t, dim, dim>
+  stress(tensor<real_t, dim, dim> dudxi,
+         PackedInternalState internal_state,
+         tensor<real_t, dim, dim> J,
+         real_t w) const
+  {
+    auto [stress, internal_state_new] = update(dudxi, internal_state, J, w);
+    return stress;
+  }
+
+  MFEM_HOST_DEVICE inline PackedInternalState
+  internal_state_new(tensor<real_t, dim, dim> dudxi,
+                     PackedInternalState internal_state,
+                     tensor<real_t, dim, dim> J,
+                     real_t w) const
+  {
+    auto [stress, internal_state_new] = update(dudxi, internal_state, J, w);
+    return internal_state_new;
+  }
+};
+
+
+// Register the custom derivatives for the solver.
+// This needs to be done for every residual function that the solver is applied on,
+// since the SolveNewtonBisection_impl is a function template, and we need a real
+// function with an address to specify the custom derivative.
+
+// Forward mode
+__attribute__((used))
+void *  __enzyme_register_derivative_newton_bisection_on_j2[2] = {
+  (void*) mfem::internal::SolveNewtonBisection_impl<J2PlasticityResidual, J2PlasticityParameters>,
+  (void*) mfem::internal::SolveNewtonBisection_impl_fwddiff<J2PlasticityResidual, J2PlasticityParameters>
+};
+
+ // Reverse mode
+__attribute__((used))
+void* __enzyme_register_gradient_SolveNewtonBisectionJ2[3] = {
+  (void*)mfem::internal::SolveNewtonBisection_impl<J2PlasticityResidual, J2PlasticityParameters>,
+  (void*)mfem::internal::SolveNewtonBisection_impl_aug<J2PlasticityResidual, J2PlasticityParameters>,
+  (void*)mfem::internal::SolveNewtonBisection_impl_rev<J2PlasticityResidual, J2PlasticityParameters>
+};
+
+// Create free functions for Enzyme to differentiate in the tests
+// Return by value version
+tensor<real_t, 3, 3> ComputeStress(
+    J2Plasticity* material, tensor<real_t, 3, 3> dudxi, 
+    J2Plasticity::PackedInternalState Q, tensor<real_t, 3, 3> J, real_t w)
+{
+    return material->stress(dudxi, Q, J, w);
+}
+
+// Return by reference version
+void ComputeStressRef(const J2Plasticity* material, const tensor<real_t, 3, 3>& dudxi,
+                      const J2Plasticity::PackedInternalState& Q,
+                      const tensor<real_t, 3, 3>& J, real_t w,
+                      tensor<real_t, 3, 3>& sigma)
+{
+  sigma = material->stress(dudxi, Q, J, w);
+}
+
+template <int dim>
+real_t elementwise_max_norm(tensor<real_t, dim, dim> A) {
+  real_t maxval = 0;
+  for (int i = 0; i < dim; i++) {
+    for (int j = 0; j < dim; j++) {
+      maxval = std::max(std::abs(A[i][j]), maxval);
+    }
+  }
+  return maxval;
+}
+
+
+TEST_CASE("Univariate function solver in a qfunction", "[univar]")
+{
+    J2Plasticity material{70.0e3, 0.34, 240.0, 0.15, 1e-3};
+    tensor<real_t, 3, 3> H{{{0.947667  , 0.9785799 , 0.33229148},
+                            {0.46866846, 0.5698887 , 0.16550303},
+                            {0.3101946 , 0.68948054, 0.74676657}}};
+    J2Plasticity::PackedInternalState Q{};
+    const tensor<real_t, 3, 3> J = IdentityMatrix<3>();
+    const real_t w = 1.0;
+
+    SECTION("Correctness")
+    {
+      // Checks that stress after update is on the yield surface.
+      auto [stress, Q_new] = material.update(H, Q, IdentityMatrix<3>(), 1.0);
+      real_t mises = std::sqrt(1.5)*norm(dev(stress));
+      real_t eqps = Q_new[9];
+      // This test only makes sense if the displacement gradient is big enough to 
+      // cuase yielding.
+      REQUIRE(eqps > 1e-9);
+      real_t Y = FlowResistance(eqps, material.sigma_y, material.n, material.ep_0);
+      CHECK(mises == MFEM_Approx(Y, 0.0, 1e-8));
+    }
+
+    SECTION("JVP")
+    {
+      // Compare forward mode derivative to finite difference approximation
+
+      tensor<real_t, 3, 3> H_dot{{{1.0, 0.0 , 0.0},
+                                  {0.0, 0.0 , 0.0},
+                                  {0.0, 0.0 , 0.0}}};
+
+      // Enzyme directional derivative (uses custom derivative of solver)
+      auto sigma_dot = __enzyme_fwddiff<tensor<real_t, 3, 3>>((void*)ComputeStress, 
+                                                              enzyme_const, &material,
+                                                              enzyme_dup, H, H_dot,
+                                                              enzyme_const, Q,
+                                                              enzyme_const, J,
+                                                              enzyme_const, w);
+      // sigma_dot = ∂sigma / ∂H[0, 0]
+      REQUIRE(sigma_dot[0][0] > 0.0);
+
+      // Finite difference derivative approximation
+      constexpr int dim = 3;
+      real_t eps = 1e-5;
+      tensor<real_t, 3, 3> sigma = ComputeStress(&material, H, Q, J, w);
+      tensor<real_t, 3, 3> sigma_p = ComputeStress(&material, H + eps*H_dot, Q, J, w);
+      tensor<real_t, 3, 3> sigma_dot_h = (1.0/eps)*(sigma_p - sigma);
+
+      tensor<real_t, 3, 3> rel_error = sigma_dot - sigma_dot_h;
+      for (int i = 0; i < dim; i++) {
+        for (int j = 0; j < dim; j++) {
+          real_t denom = sigma[i][j] != 0? sigma[i][j] : 1.0;
+          rel_error[i][j] /= denom;
+        }
+      }
+
+      CHECK(elementwise_max_norm(rel_error) < 1e-5);
+    }
+
+    SECTION("VJP")
+    {
+      // compare reverse mode derivative to finite differences
+
+      tensor<real_t, 3, 3> sigma;
+      ComputeStressRef(&material, H, Q, J, w, sigma);
+      double epsilon = 1e-6;
+      tensor<real_t, 3, 3> dH{{{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}, {0.0, 0.0, 0.0}}};
+      auto H_p = H + epsilon*dH;
+      tensor<real_t, 3, 3> sigma_p;
+      ComputeStressRef(&material, H_p, Q, J, w, sigma_p);
+      auto sigma_dot_h = (sigma_p - sigma)/epsilon;
+      // Note: sigma_dot_h[i,j] = ∂sigma[i,j]/∂H[0,0]
+
+      tensor<real_t, 3, 3> sigma_bar{{{1.0, 0.0, 0.0},
+                                      {0.0, 0.0, 0.0},
+                                      {0.0, 0.0, 0.0}}};
+
+      tensor<real_t, 3, 3> H_bar{};
+      J2Plasticity::PackedInternalState Q_bar{};
+      tensor<real_t, 3, 3> J_bar{};
+      __enzyme_autodiff<void>(
+        (void*)ComputeStressRef,  enzyme_const, &material, enzyme_dup, &H, &H_bar,
+        enzyme_dup, &Q, &Q_bar, enzyme_dup, &J, &J_bar, enzyme_const, w,
+        enzyme_dup, &sigma, &sigma_bar);
+
+      // H_bar[ij] = ∂sigma[0,0]/∂H[i,j]
+      // For this model, we expect the major symmetries in the tangent operator.
+      // Hence H_bar \approx sigma_dot_h
+
+      const double abs_tol = 1e-12;
+      const double rel_tol = 5e-6;
+      for (int i = 0; i < 3; i++) {
+        for (int j = 0; j < 3; j++) {
+          CHECK(H_bar[i][j] == MFEM_Approx(sigma_dot_h[i][j], abs_tol, rel_tol));
+        }
+      }
+    }
+}
+
+
+real_t nthroot_res(real_t x, tuple<real_t, real_t> p)
+{
+    auto [index, radicand] = p;
+    return std::pow(x, index) - radicand;
+}
+
+__attribute__((used))
+void* __enzyme_register_gradient_solver[3] = {
+  (void*)mfem::internal::SolveNewtonBisection_impl<nthroot_res, tuple<real_t, real_t>>,
+  (void*)mfem::internal::SolveNewtonBisection_impl_aug<nthroot_res, tuple<real_t, real_t>>,
+  (void*)mfem::internal::SolveNewtonBisection_impl_rev<nthroot_res, tuple<real_t, real_t>>
+};
+
+TEST_CASE("Univariate solver reverse mode derivative", "[univar]")
+{
+  auto mysqrt = [](real_t x) -> real_t
+    {
+      real_t x0 = x;
+      real_t index = 2.0;
+      real_t ub = std::max(1.0, x);
+      SolverSettings settings{1e-12, 1e-12, {0, ub}};
+      return SolveNewtonBisection<nthroot_res>(x0, make_tuple(index, x), settings);
+    };
+
+  real_t x = 2.0;
+  real_t dydx = __enzyme_autodiff<real_t>((void*)+mysqrt, enzyme_out, x);
+  CHECK(dydx == MFEM_Approx(0.5/std::sqrt(2.0)));
+}
+
+TEST_CASE("Univariate function solver robustness", "[univar]")
+{
+  SolverSettings settings{1e-12, 1e-12};
+
+  SECTION("Simple case")
+  {
+      auto Nthroot = [&settings](real_t x, real_t n) {
+          real_t x0 = std::max(x, 1.0);
+          settings.bounds = {0.0, x0};
+          return SolveNewtonBisection<nthroot_res>(x0, make_tuple(n, x), settings);
+      };
+      real_t x = 8.0;
+      real_t y = Nthroot(x, 3.0);
+      CHECK(y == MFEM_Approx(2.0));
+  }
+
+    SECTION("Stiff problem")
+    {
+        auto f = [](real_t x, real_t p) { return std::pow(x, p) - 1.0; };
+        real_t x0 = 0.1;
+        real_t p = 50;
+        settings.bounds = {0.0, 5.1};
+        real_t x = SolveNewtonBisection<+f>(x0, p, settings);
+        CHECK(x == MFEM_Approx(1.0));
+    }
+
+    SECTION("Works where Newton diverges")
+    {
+        auto f = [](double x, int) { return std::atan(x); };
+        real_t x0 = 1.5;
+        settings.bounds = {0.0, 2.0};
+        real_t x = SolveNewtonBisection<+f>(x0, int{}, settings);
+        CHECK(std::abs(x) == MFEM_Approx(0.0));
+    }
+}
+
+#endif // MFEM_USE_ENZYME
@@ -105,39 +105,39 @@ TEST_CASE("Integration rule order initialization", "[IntegrationRules]")
   SECTION("Segment rule constructed by accessing square rule")
   {
      auto &quad5_ir = intrules.Get(Geometry::SQUARE, 5);
-      REQUIRE(quad5_ir.GetOrder() >= 5);
+      REQUIRE(quad5_ir.GetOrder() == 5);
      // The segment integration rule of order 5 is lazy constructed when we get
      // the square integration rule of order 5. Make sure its order was
      // properly set:
      auto &line5_ir = intrules.Get(Geometry::SEGMENT, 5);
-      REQUIRE(line5_ir.GetOrder() >= 5);
+      REQUIRE(line5_ir.GetOrder() == 5);
   }

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

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

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

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

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

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

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

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

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

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

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

   // Assemble the linear form
   ParLinearForm lf(&fespace);
@@ -200,39 +200,3 @@ TEST_CASE("ArraysByName Sort/Unique Methods", "[ArraysByName]")
      }
   }
 }
-
-TEST_CASE("ArraysByName Print/Load Methods", "[ArraysByName]")
-{
-   ArraysByName<int> abn;
-
-   FillArraysByName(abn);
-
-   // Print object to string using default format
-   std::ostringstream oss1;
-   abn.Print(oss1);
-
-   // Load new object from printed output
-   ArraysByName<int> abn_load1;
-   std::istringstream iss1(oss1.str());
-   abn_load1.Load(iss1);
-   REQUIRE(abn == abn_load1);
-
-   // Print object to string using one line per array
-   std::ostringstream oss2;
-   oss2 << abn.Size() << '\n';
-   for (auto a : abn)
-   {
-      oss2 << '"' << a.first << "\" " << a.second.Size();
-      for (auto d : a.second)
-      {
-         oss2 << ' ' << d;
-      }
-      oss2 << '\n';
-   }
-
-   // Load new object from printed output
-   ArraysByName<int> abn_load2;
-   std::istringstream iss2(oss2.str());
-   abn_load2.Load(iss2);
-   REQUIRE(abn == abn_load2);
-}
@@ -1,74 +0,0 @@
-// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
-// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
-// LICENSE and NOTICE for details. LLNL-CODE-806117.
-//
-// This file is part of the MFEM library. For more information and source code
-// availability visit https://mfem.org.
-//
-// MFEM is free software; you can redistribute it and/or modify it under the
-// terms of the BSD-3 license. We welcome feedback and contributions, see file
-// CONTRIBUTING.md for details.
-
-#include "mfem.hpp"
-#include "unit_tests.hpp"
-
-using namespace mfem;
-
-#if defined(MFEM_USE_MPI) && defined(MFEM_USE_PETSC)
-
-namespace
-{
-struct PetscSession
-{
-   PetscSession() { MFEMInitializePetsc(); }
-   ~PetscSession() { MFEMFinalizePetsc(); }
-};
-
-class IdentityGradientOperator : public IdentityOperator
-{
-public:
-   IdentityGradientOperator() : IdentityOperator(1), _jac(1)
-   {
-      _jac.Add(0, 0, 1.0);
-      _jac.Finalize();
-   }
-
-   Operator &GetGradient(const Vector &) const override
-   {
-      return const_cast<SparseMatrix &>(_jac);
-   }
-
-private:
-   SparseMatrix _jac;
-};
-}
-
-TEST_CASE("PetscNonlinearSolver accepts non-empty rhs", "[Parallel][PETSc]")
-{
-   static PetscSession petsc_session;
-
-   IdentityGradientOperator oper;
-   PetscNonlinearSolver solver(MPI_COMM_WORLD, "nl_");
-   solver.SetRelTol(1.0e-12);
-   solver.SetAbsTol(1.0e-12);
-   solver.SetMaxIter(5);
-   solver.SetPrintLevel(0);
-   solver.SetJacobianType(Operator::PETSC_MATAIJ);
-   solver.SetOperator(oper);
-
-   Vector x(1);
-
-   Vector empty_rhs;
-   x = 0.0;
-   solver.Mult(empty_rhs, x);
-   REQUIRE(x(0) == MFEM_Approx(0.0));
-
-   Vector nonempty_rhs(1);
-   nonempty_rhs(0) = 2.5;
-   x = 0.0;
-   solver.Mult(nonempty_rhs, x);
-   REQUIRE(x.Size() == 1);
-   REQUIRE(x(0) == MFEM_Approx(nonempty_rhs(0)));
-}
-
-#endif
@@ -486,14 +486,8 @@ void multidomain_test_3d(FECType fec_type)
         {
            cylinder_gf.ProjectCoefficient(vcoeff);
            outer_gf.ProjectCoefficient(vcoeff);
-            if (fec_type == FECType::RT)
-            {
-               outer_gf.ProjectBdrCoefficientNormal(vzerocoeff, outer_cyl_surf_marker);
-            }
-            else
-            {
-               outer_gf.ProjectBdrCoefficientTangent(vzerocoeff, outer_cyl_surf_marker);
-            }
+            outer_gf.ProjectBdrCoefficient(vzerocoeff,
+                                           outer_cyl_surf_marker);
            outer_gf_ex.ProjectCoefficient(vcoeff);
         }
         ParSubMesh::Transfer(cylinder_gf, outer_gf);
@@ -513,14 +507,8 @@ void multidomain_test_3d(FECType fec_type)
         {
            outer_gf.ProjectCoefficient(vcoeff);
            cylinder_gf.ProjectCoefficient(vcoeff);
-            if (fec_type == FECType::RT)
-            {
-               cylinder_gf.ProjectBdrCoefficientNormal(vzerocoeff, cylinder_cyl_surf_marker);
-            }
-            else
-            {
-               cylinder_gf.ProjectBdrCoefficientTangent(vzerocoeff, cylinder_cyl_surf_marker);
-            }
+            cylinder_gf.ProjectBdrCoefficient(vzerocoeff,
+                                              cylinder_cyl_surf_marker);
            cylinder_gf_ex.ProjectCoefficient(vcoeff);
         }
         ParSubMesh::Transfer(outer_gf, cylinder_gf);
Author	SHA1	Message	Date
Brandon Talamini	4b13e46bc5	Fix style error	2026-03-29 17:12:36 -07:00
Brandon Talamini	de114112b5	Remove unused variable from test	2026-03-27 21:08:43 -07:00
Brandon Talamini	ed1c21f656	Apply style to tensor.hpp	2026-03-27 21:01:16 -07:00
Brandon Talamini	811254989f	Get rid of designated initializers	2026-03-27 20:56:18 -07:00
Brandon Talamini	1f639a826f	Fix ODR problems caused by dropped inline keywords	2026-03-27 20:32:59 -07:00
Brandon Talamini	41e6f97bc7	Use relative paths for all includes	2026-03-27 20:02:46 -07:00
Brandon Talamini	7045b49973	Fix path to config.hpp	2026-03-27 19:57:40 -07:00
Brandon Talamini	91a121e427	Rename a badly named variable	2026-03-27 19:43:08 -07:00
Brandon Talamini	1f6d6ba82a	Style	2026-03-27 17:31:20 -07:00
Brandon Talamini	379ec4a393	Move test to dfem directory	2026-03-27 17:09:15 -07:00
Brandon Talamini	34895e545e	Cleanup through comments and renaming	2026-03-27 17:02:40 -07:00
Brandon Talamini	e0659eecd4	Guard tests for enzyme build only	2026-03-27 17:02:40 -07:00
Brandon Talamini	20b774792a	reorder things in test for clarity	2026-03-27 17:02:40 -07:00
Brandon Talamini	227bbe0283	Remove debug prints and add some comments	2026-03-27 17:02:40 -07:00
Brandon Talamini	88e8055d54	Remove some debug output and simplify test	2026-03-27 17:02:40 -07:00
Brandon Talamini	8e968658b5	Clean up tolerances and settings	2026-03-27 17:02:40 -07:00
Brandon Talamini	978e2bbd99	Fix Enzyme compilation error in Release mode Credit to William Moses.	2026-03-27 17:02:40 -07:00
Brandon Talamini	6465d872b5	Comment and move things around for clarity	2026-03-27 17:02:40 -07:00
Brandon Talamini	0439e51f6b	Make bounds for nth root functions better so they will work for any input	2026-03-27 17:02:40 -07:00
Brandon Talamini	e2917c9e15	Clamp initial guess with brackets, no warning if guess is moved	2026-03-27 17:02:40 -07:00
Brandon Talamini	f0ea358c36	Change lack of root bracketing handling from an assertion to a warning so the user is alerted in Release builds	2026-03-27 17:02:40 -07:00
Brandon Talamini	c26460c0ce	Change divergence error to an MFEM_ABORT macro	2026-03-27 17:02:40 -07:00
Brandon Talamini	9975f4c9db	Inline tensor operations so that Enzyme can analyze enzyme_const operations correctly	2026-03-27 17:02:40 -07:00
Brandon Talamini	b29d5856c2	Solve problem where setting bounds that depended on state caused incorrect derivatives	2026-03-27 17:02:40 -07:00
Brandon Talamini	ce54fdf7b4	Put in a robustness test for a case where standard Newton diverges	2026-03-27 17:02:40 -07:00
Brandon Talamini	7ca58c106a	Clean up tests: comment and use non-halting test conditions	2026-03-27 17:02:39 -07:00
Brandon Talamini	63056ab08b	Put in a real test condition for the VJP	2026-03-27 17:02:39 -07:00
Brandon Talamini	cd9456db2a	Get reverse mode working Derivatives become nans if solver upper bound depends on state. Need to debug this with Bill. But derivatives are now being computed, and appear correct.	2026-03-27 17:02:39 -07:00
Brandon Talamini	7a8fb3ca79	Fixing Enzyme problems with Bill Moses	2026-03-27 17:02:39 -07:00
Brandon Talamini	dc251d9518	Move solver from linalg to dfem	2026-03-27 17:02:39 -07:00
Brandon Talamini	2c118aaf85	Move implmentation inside separate namespace	2026-03-27 17:02:39 -07:00
Brandon Talamini	e0e32976ec	Fix custom revrse mode	2026-03-27 17:02:39 -07:00
Brandon Talamini	ac71a377e1	Make wrapper for residual use values and references so that it is more symmetric with residual signature	2026-03-27 17:02:39 -07:00
Brandon Talamini	ac9b4d3d38	Make a basic test of reverse mode and get it to pass	2026-03-27 17:02:39 -07:00
Brandon Talamini	6be47674d6	Implement custom reverse mode derivative of solver	2026-03-27 17:02:39 -07:00
Brandon Talamini	75b6f2ac19	Check enzyme derivative with finite difference	2026-03-27 17:02:39 -07:00
Brandon Talamini	76c805cbad	Put in a test that computes derivative of q-function with solver inside	2026-03-27 17:02:39 -07:00
Brandon Talamini	592f87591d	Fix bugs to make tests pass	2026-03-27 17:02:39 -07:00
Brandon Talamini	5decf944cb	Fix bugs in error checking Forgot to reverse inequalities when changing from ERROR_IF to ASSERT checks.	2026-03-27 17:02:39 -07:00
Brandon Talamini	6eccede39e	Add a unit test of new solver	2026-03-27 17:02:39 -07:00
Brandon Talamini	f8fd930618	Fix compilation errors	2026-03-27 17:02:39 -07:00
Brandon Talamini	c450622e6e	Initial draft of solver	2026-03-27 17:02:39 -07:00