style

Updates CHANGELOG
More fixes for building contact mini-app with make and cmake
2026-04-14 08:54:23 -07:00 · 2026-04-13 21:52:37 -07:00 · 2026-04-13 21:47:15 -07:00 · 2026-04-13 18:54:17 -07:00 · 2026-04-13 12:40:24 -07:00 · 2026-04-13 12:37:18 -07:00
334 changed files with 8486 additions and 35116 deletions
@@ -94,16 +94,6 @@ inputs:
    description: If true, do not set any CXXFLAGS or LDFLAGS.
    default: false

-  # Unfortunately, "uses:" fields cannot have references to variables like
-  # ${{env.MFEM_ACTIONS_VERSION}}, so the branch/tag name has to be hard coded.
-  # Therefore, in the future, when updating the version of the
-  # mfem/github-actions to use, we'll have to replace:
-  # - all definitions of MFEM_ACTIONS_VERSION and
-  # - all "uses:" fields that refer to mfem/github-actions.
-  MFEM_ACTIONS_VERSION:
-    description: Version (branch or tag) of the mfem/github-actions to use.
-    default: v2.7
-
 runs:
  using: 'composite'
  steps:
@@ -128,7 +118,6 @@ runs:
        echo UBSAN_LDFLAGS=${{inputs.UBSAN_LDFLAGS}} >> $GITHUB_ENV
        echo MSAN_CXXFLAGS=${{inputs.MSAN_CXXFLAGS}} >> $GITHUB_ENV
        echo MSAN_LDFLAGS=${{inputs.MSAN_LDFLAGS}} >> $GITHUB_ENV
-        echo MFEM_ACTIONS_VERSION=${{inputs.MFEM_ACTIONS_VERSION}} >> $GITHUB_ENV
      shell: bash

    - name: Env (dir)
@@ -53,7 +53,7 @@ runs:
      run: echo CXXFLAGS=${{env.CXXFLAGS}} ${{env.UBSAN_CXXFLAGS}} >> $GITHUB_ENV
      shell: bash

-    - uses: mfem/github-actions/build-mfem@v2.7
+    - uses: mfem/github-actions/build-mfem@v2.5
      if: ${{steps.debug.outputs.cache-hit != 'true'}}
      env:
        CXXFLAGS: ${{env.CXXFLAGS}}
@@ -82,7 +82,7 @@ runs:
      run: find . -type f -name '*.o' -delete
      shell: bash

-    - uses: actions/upload-artifact@v7
+    - uses: actions/upload-artifact@v4
      with:
        name: build-${{inputs.par}}-${{inputs.sanitizer}}
        path: mfem/build
@@ -12,11 +12,6 @@
 name: 'Install MPI'
 description: 'Installs MPI and set up its environment variables'

-inputs:
-  NO_FLAGS:
-    description: If true, do not set any CXXFLAGS or LDFLAGS.
-    default: false
-
 runs:
  using: 'composite'
  steps:
@@ -32,7 +27,6 @@ runs:
      shell: bash

    - name: Env (bis)
-      if: ${{ inputs.NO_FLAGS != 'true' }}
      run: |
        echo CXXFLAGS=${{env.CXXFLAGS}} ${{env.MPI_INC}} >> $GITHUB_ENV
        echo LDFLAGS=${{env.LDFLAGS}} ${{env.MPI_LIB}} >> $GITHUB_ENV
@@ -49,7 +49,7 @@ runs:
        par: ${{inputs.par}}
        sanitizer: ${{inputs.sanitizer}}

-    - uses: actions/download-artifact@v8
+    - uses: actions/download-artifact@v4
      with:
        name: build-${{inputs.par}}-${{inputs.sanitizer}}
        path: mfem/build
@@ -37,14 +37,14 @@ runs:
      with:
        path: ${{env.HYPRE_DIR}}
        fail-on-cache-miss: true
-        key: ${{runner.os}}-ompi-build-${{env.HYPRE_DIR}}-int32-fp64-${{env.MFEM_ACTIONS_VERSION}}
+        key: ${{runner.os}}-ompi-build-${{env.HYPRE_DIR}}-int32-fp64-v2.5

    - uses: actions/cache/restore@v5  # Cache for Metis
      if: ${{inputs.par == 'true'}}
      with:
        path: ${{env.METIS_DIR}}
        fail-on-cache-miss: true
-        key: ${{runner.os}}-build-${{env.METIS_DIR}}-${{env.MFEM_ACTIONS_VERSION}}
+        key: ${{runner.os}}-build-${{env.METIS_DIR}}-v2.5

    - name: Hypre/Metis links
      if: ${{inputs.par == 'true'}}
@@ -1,42 +0,0 @@
-# MFEM Pull Request Review Agent Guide
-
-## Purpose and scope
-Review MFEM PRs for correctness, maintainability, performance, portability, test coverage, and MFEM consistency. Use the diff and PR context; reference source files, tests, and CI results when available. Follow `CONTRIBUTING.md`, especially Developer Guidelines, PR rules, checklist, and testing.
-
-## Critical review pillars
- Correctness and numerical behavior
- API and user-facing impact
- Performance implications
- Maintainability and portability
-
-## Review workflow
-1. Read the PR description, linked issues, and intended behavior.
-2. Inspect the diff before commenting.
-3. Identify affected MFEM components, examples, tests, build or docs changes, and downstream APIs.
-4. Analyze the code against the critical review pillars.
-5. Compare the change against nearby code and MFEM patterns; flag unmotivated deviations.
-6. Check whether tests and documentation were updated appropriately.
-7. Review CI results and suggest actions.
-8. Produce a structured review with prioritized findings.
-9. Always limit conclusions to available evidence.
-
-## MFEM-specific review checklist
- Component-aware scope: identify the touched subsystem (FEM, solvers, preconditioners, linear algebra, mesh, examples, miniapps, build, or docs) and assess its impact against the review pillars.
- Numerical and algorithmic behavior: assess issues in convergence, stability, tolerances, precision, iteration limits, and failure handling. If clear opportunities exist to improve the algorithmic approach, call them out with expected impact.
- API and user-facing impact: assess backward compatibility, user-visible behavior and default changes, migration impact, deprecations, and whether documentation clearly explains user-facing API changes.
- Data structure and memory semantics: assess ownership, lifetime, aliasing, container behavior, and device-host synchronization.
- Parallel and serial behavior: assess whether the change preserves equivalent semantics in serial and parallel modes where applicable; if logic is currently mode-specific, check whether extension to the other mode is straightforward (clear abstractions, no hard-wired assumptions), document constraints, and call out expected behavior differences explicitly.
- Backend and portability impact: assess likely cross-backend risks in CPU, CUDA, HIP, OCCA, RAJA, partial assembly, fallback paths, compiler compatibility, and platform assumptions.
- Build, dependency, and configuration impact: assess CMake or make changes, optional dependency behavior, and feature-flag interactions.
- Tests and docs alignment: check available regression or unit coverage evidence for changed behavior, and ensure docs are updated for new flags, APIs, options, or behavior changes.
- MFEM developer-guideline fit: keep code lean, simple, general, logically separated, and portable; suggest C++17 improvements when they clearly improve safety, clarity, or maintainability.
- New source files, examples, or miniapps: if a PR adds source/header files, verify they are properly wired into the relevant `makefile` and `CMakeLists.txt`, referenced in docs where applicable (including `doc/CodeDocumentation.dox`), and added to top-level `.gitignore` only when generated artifacts require it.
- Changelog: verify `CHANGELOG` is updated if the PR introduces significant new features or user-facing changes.
- MFEM conventions: use `real_t`; use `mfem::out`/`mfem::err` instead of `std::cout`/`std::cerr` in library code; flag large/binary files; if AI assistance is apparent but undisclosed, suggest following `CONTRIBUTING.md`.
- Edge cases: if the PR touches complex or error-prone areas, suggest additional tests for edge cases, failure modes, and parallel behavior.
-
-## Commenting guidelines
- Keep comments concise, actionable, and grounded in the diff.
- Focus on correctness, behavior changes, and user impact over style nits.
- Be professional, concise, collaborative, technically precise, and avoid unsupported assumptions.
-
@@ -13,7 +13,7 @@ Note that some of these scripts use the shared MFEM GitHub Actions from the exte

  <https://github.com/mfem/github-actions>

-For a particular action, e.g. `mfem/github-actions/build-mfem@v2.5`, the `v2.5` suffix denotes the branch (or tag) in the above from which the action is taken.
+For a particular action, e.g. `mfem/github-actions/build-mfem@v2.5`, the `v2.5` suffix denotes the branch in the above from which the action is taken.

 The current CI workflows are:

@@ -29,12 +29,16 @@ Runs a number of static repository-level sanity checks.

 - `branch-history` guards against accidental commits of large files using the `--history` option of the `config/githooks/pre-push` script.

+## `mfem-analysis.yml` (`build-analysis`)
+
+Checks if the code builds and satisfies minimal requirements.
+
+- `gitignore` builds hypre, METIS, and MFEM using `mfem/github-actions/build-hypre`, `mfem/github-actions/build-metis`, and `mfem/github-actions/build-mfem` and checks for correct `.gitignore` settings by running the `tests/scripts/gitignore` script.
+
 ## `builds-and-tests.yml`

 Runs a matrix of builds and tests runs with different compilers, OS, mfem/hypre settings, etc. Also processes and upload Codecov reports.

-One matrix job runs `tests/scripts/gitignore` after `make test-noclean` to check generated artifacts against `.gitignore`.
-
 Uses the following GitHub Actions from <https://github.com/mfem/github-actions>:

 - `mfem/github-actions/build-hypre`
@@ -40,7 +40,6 @@ env:
  METIS_ARCHIVE_MAC: metis-4.0.3-mac.tgz
  METIS_TOP_DIR: metis-4.0.3
  MFEM_TOP_DIR: mfem
-  MFEM_ACTIONS_VERSION: v2.7

 # Note for future improvements:
 #
@@ -111,7 +110,6 @@ jobs:
            build-system: make
            hypre-target: int64
            precision: fp64
-            gitignore-check: YES
          - os: ubuntu-latest
            target: opt
            codecov: NO
@@ -142,10 +140,6 @@ jobs:

    continue-on-error: ${{ matrix.enzyme && true || false }}

-    # Enable ccache for all jobs except Windows (would need sccache).
-    env:
-      USE_CCACHE: ${{ matrix.os != 'windows-latest' }}
-
    steps:
      # Fix 'No space left on device' errors for Ubuntu builds.
      - name: Run Actions Cleaner
@@ -176,6 +170,20 @@ jobs:
          env
        shell: bash

+      # For info on Xcode see:
+      #  - https://github.com/actions/runner-images/issues/12541
+      #  - https://github.com/actions/runner-images/blob/releases/macos-15-arm64/20250811/images/macos/macos-15-arm64-Readme.md#xcode
+      - name: Xcode version setup (MacOS)
+        if: matrix.os == 'macos-latest'
+        run: |
+          XCODE_PATH="/Applications/Xcode_16.4.app"
+          echo "> sudo xcode-select -s ${XCODE_PATH}"
+          sudo xcode-select -s ${XCODE_PATH}
+          echo "> g++ -v"
+          g++ -v
+          echo "> clang++ -v"
+          clang++ -v
+
      # Only get MPI if defined for the job.
      # TODO: It would be nice to have only one step, e.g. with a dedicated
      # action, but I (@adrienbernede) don't see how at the moment.
@@ -220,11 +228,11 @@ jobs:
        uses: actions/cache@v5
        with:
          path: ${{ env.HYPRE_TOP_DIR }}
-          key: ${{ runner.os }}-ompi-build-${{ env.HYPRE_TOP_DIR }}-${{ matrix.hypre-target }}-${{ matrix.precision }}-${{ env.MFEM_ACTIONS_VERSION }}
+          key: ${{ runner.os }}-ompi-build-${{ env.HYPRE_TOP_DIR }}-${{ matrix.hypre-target }}-${{ matrix.precision }}-v2.5

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

      - name: get hypre (Windows)
        if: matrix.mpi == 'par' && steps.hypre-cache.outputs.cache-hit != 'true' && matrix.os == 'windows-latest'
-        uses: mfem/github-actions/build-hypre@v2.7
+        uses: mfem/github-actions/build-hypre@v2.5
        with:
          archive: ${{ env.HYPRE_ARCHIVE }}
          dir: ${{ env.HYPRE_TOP_DIR }}
@@ -250,11 +258,11 @@ jobs:
        uses: actions/cache@v5
        with:
          path: ${{ env.METIS_TOP_DIR }}
-          key: ${{ runner.os }}-build-${{ env.METIS_TOP_DIR }}-${{ env.MFEM_ACTIONS_VERSION }}
+          key: ${{ runner.os }}-build-${{ env.METIS_TOP_DIR }}-v2.5

      - name: install metis
        if: matrix.mpi == 'par' && matrix.os != 'windows-latest' && steps.metis-cache.outputs.cache-hit != 'true'
-        uses: mfem/github-actions/build-metis@v2.7
+        uses: mfem/github-actions/build-metis@v2.5
        with:
          archive: ${{ matrix.os != 'macos-latest' && env.METIS_ARCHIVE || env.METIS_ARCHIVE_MAC }}
          dir: ${{ env.METIS_TOP_DIR }}
@@ -294,55 +302,9 @@ jobs:
          echo "OMPI_CC=$LLVM_PREFIX/bin/clang" >> $GITHUB_ENV
          echo "OMPI_CXX=$LLVM_PREFIX/bin/clang++" >> $GITHUB_ENV

-      # Restore the compiler cache (ccache). The key embeds the run id, so new
-      # runs save a fresh snapshot; the restore-keys prefix warm-starts from the
-      # most recent prior run (incl. the base branch for PRs).
-      - name: cache ccache
-        if: ${{ env.USE_CCACHE == 'true' }}
-        uses: actions/cache@v5
-        with:
-          path: .ccache
-          key: ccache-${{ matrix.os }}-${{ matrix.build-system }}-${{ matrix.target }}-${{ matrix.mpi }}-${{ matrix.hypre-target }}-${{ matrix.precision }}${{ matrix.enzyme && '-enzyme' || '' }}-${{ github.run_id }}
-          restore-keys: |
-            ccache-${{ matrix.os }}-${{ matrix.build-system }}-${{ matrix.target }}-${{ matrix.mpi }}-${{ matrix.hypre-target }}-${{ matrix.precision }}${{ matrix.enzyme && '-enzyme' || '' }}-
-
-      # Configure ccache and select how it is injected into the MFEM build:
-      # - make:  set CXX="ccache g++"; for MPI, OMPI_CXX="ccache g++" so mpicxx
-      #          runs ccache around g++ (not ccache around the mpicxx wrapper).
-      # - cmake: set CMAKE_<LANG>_COMPILER_LAUNCHER=ccache.
-      # - enzyme: wrap the brew clang++ via OMPI_CXX.
-      # The chosen options are passed through build-mfem's 'config-options'
-      # input (see the build step below).
-      - name: configure ccache
-        if: ${{ env.USE_CCACHE == 'true' }}
-        run: |
-          command -v ccache >/dev/null 2>&1 || {
-            if [[ "${{ runner.os }}" == "Linux" ]]; then
-              sudo apt-get update && sudo apt-get install -y ccache
-            else
-              brew install ccache
-            fi
-          }
-          echo "CCACHE_DIR=${{ github.workspace }}/.ccache" >> $GITHUB_ENV
-          echo "CCACHE_MAXSIZE=1G" >> $GITHUB_ENV
-          echo "CCACHE_COMPILERCHECK=content" >> $GITHUB_ENV
-          # Ignore header timestamps (restamped by each checkout) so direct mode hits.
-          echo "CCACHE_SLOPPINESS=include_file_mtime,include_file_ctime,time_macros" >> $GITHUB_ENV
-          # Hash absolute paths relative to the workspace.
-          echo "CCACHE_BASEDIR=${{ github.workspace }}" >> $GITHUB_ENV
-          if [[ "${{ matrix.enzyme }}" == "true" ]]; then
-            echo "OMPI_CXX=ccache $LLVM_PREFIX/bin/clang++" >> $GITHUB_ENV
-          elif [[ "${{ matrix.build-system }}" == "cmake" ]]; then
-            echo 'CCACHE_CONFIG_OPTS=-DCMAKE_CXX_COMPILER_LAUNCHER=ccache -DCMAKE_C_COMPILER_LAUNCHER=ccache' >> $GITHUB_ENV
-          else
-            echo "OMPI_CXX=ccache g++" >> $GITHUB_ENV
-            echo 'CCACHE_CONFIG_OPTS=CXX="ccache g++" MPICXX="mpicxx"' >> $GITHUB_ENV
-          fi
-        shell: bash
-
      # MFEM build and test
      - name: build
-        uses: mfem/github-actions/build-mfem@v2.7
+        uses: mfem/github-actions/build-mfem@v2.5
        env:
          VCPKG_DEFAULT_BINARY_CACHE: ${{ github.workspace }}/vcpkg_cache
        with:
@@ -355,14 +317,9 @@ jobs:
          metis-dir: ${{ env.METIS_TOP_DIR }}
          mfem-dir: ${{ env.MFEM_TOP_DIR }}
          precision: ${{ matrix.precision }}
-          config-options: ${{ matrix.config-opts }} ${{ env.CCACHE_CONFIG_OPTS }}
+          config-options: ${{ matrix.config-opts }}
          library-only: ${{ matrix.target == 'dbg' && matrix.os != 'ubuntu-latest' }}

-      - name: ccache stats
-        if: ${{ env.USE_CCACHE == 'true' }}
-        run: ccache -s
-        shell: bash
-
      # Run checks (and only checks) on debug targets
      - name: checks
        if: matrix.build-system == 'make' && matrix.target == 'dbg'
@@ -373,13 +330,7 @@ jobs:
      - name: tests
        if: matrix.build-system == 'make' && (matrix.target == 'opt' || matrix.os == 'ubuntu-latest')
        run: |
-          cd ${{ env.MFEM_TOP_DIR }}
-          if [[ "${{ matrix.gitignore-check }}" == "YES" ]]; then
-            make test-noclean
-          else
-            make test
-          fi
-        shell: bash
+          cd ${{ env.MFEM_TOP_DIR }} && make test

      - name: cmake checks
        if: matrix.build-system == 'cmake' && matrix.target == 'dbg'
@@ -424,16 +375,8 @@ jobs:
      # Code coverage (process and upload reports)
      - name: codecov
        if: matrix.codecov == 'YES'
-        uses: mfem/github-actions/upload-coverage@v2.7
+        uses: mfem/github-actions/upload-coverage@v2.5
        with:
-          name: ${{ matrix.os }}-${{ matrix.build-system }}-${{ matrix.target }}-${{ matrix.mpi }}-${{ matrix.hypre-target }}-${{ matrix.precision }}
+          name: ${{ matrix.os }}-${{ matrix.build-system }}-${{ matrix.target }}-${{ matrix.mpi }}-${{ matrix.hypre-target }}
          project_dir: ${{ env.MFEM_TOP_DIR }}
          directories: "fem general linalg mesh"
-        env:
-          CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }}
-
-      - name: gitignore
-        if: matrix.gitignore-check == 'YES'
-        run: |
-          cd ${{ env.MFEM_TOP_DIR }}/tests/scripts
-          ./runtest gitignore
@@ -1,42 +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.
---
-# A closed PR's caches can never be restored again, so delete them to free
-# space against the 10 GB per-repo cache limit.
-name: Cleanup PR caches
-
-on:
-  pull_request:
-    types: [closed]
-
-permissions:
-  actions: write
-
-jobs:
-  cleanup:
-    runs-on: ubuntu-latest
-    steps:
-      - name: Delete caches for the closed PR
-        env:
-          GH_TOKEN: ${{ secrets.GITHUB_TOKEN }}
-          GH_REPO: ${{ github.repository }}
-          PR_REF: refs/pull/${{ github.event.pull_request.number }}/merge
-        run: |
-          echo "Deleting caches for $PR_REF"
-          while :; do
-            ids=$(gh cache list --ref "$PR_REF" --limit 100 --json id --jq '.[].id')
-            [ -n "$ids" ] || break
-            echo "$ids" | while read -r id; do
-              [ -n "$id" ] || continue
-              echo "Deleting cache $id"
-              gh cache delete "$id" || echo "  (already gone)"
-            done
-          done
@@ -14,19 +14,9 @@ name: "Static Analysis"
 on:
  push:
    branches: ["master", "next"]
-    paths-ignore: &docs-only-paths
-      - "**/*.md"
-      - "doc/**"
-      - ".binder/**"
-      - "CITATION.cff"
-      - "LICENSE"
-      - "NOTICE"
-      - "CHANGELOG"
-      - "INSTALL"
  pull_request:
    # The branches below must be a subset of the branches above
    branches: ["master"]
-    paths-ignore: *docs-only-paths

 concurrency:
  group: ${{ github.workflow }}-${{ github.ref }}
@@ -0,0 +1,100 @@
+# 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.
+
+name: "Build Analysis"
+
+permissions:
+  actions: write
+
+on:
+  push:
+    branches:
+      - master
+      - next
+  pull_request:
+  workflow_dispatch:
+
+concurrency:
+  group: ${{ github.workflow }}-${{ github.ref }}
+  cancel-in-progress: true
+
+env:
+  HYPRE_ARCHIVE: v2.19.0.tar.gz
+  HYPRE_TOP_DIR: hypre-2.19.0
+  METIS_ARCHIVE: metis-4.0.3.tar.gz
+  METIS_TOP_DIR: metis-4.0.3
+  COVERAGE_ENV: mfem-coverage
+
+jobs:
+  gitignore:
+    runs-on: ubuntu-latest
+
+    steps:
+      - name: checkout MFEM
+        uses: actions/checkout@v6
+        with:
+          path: mfem
+
+      - name: Get MPI (Linux)
+        run: |
+          sudo apt-get install openmpi-bin libopenmpi-dev
+          export OMPI_MCA_rmaps_base_oversubscribe=1
+
+      - name: Cache Hypre Install
+        id: hypre-cache
+        uses: actions/cache@v5
+        with:
+          path: ${{ env.HYPRE_TOP_DIR }}
+          key: ${{ runner.os }}-ompi-build-${{ env.HYPRE_TOP_DIR }}-v2.5
+
+      - name: Get Hypre
+        if: steps.hypre-cache.outputs.cache-hit != 'true'
+        uses: mfem/github-actions/build-hypre@v2.5
+        with:
+          archive: ${{ env.HYPRE_ARCHIVE }}
+          dir: ${{ env.HYPRE_TOP_DIR }}
+          target: int32
+
+      - name: Cache Metis Install
+        id: metis-cache
+        uses: actions/cache@v5
+        with:
+          path: ${{ env.METIS_TOP_DIR }}
+          key: ${{ runner.os }}-build-${{ env.METIS_TOP_DIR }}-v2.5
+
+      - name: Install Metis
+        if: steps.metis-cache.outputs.cache-hit != 'true'
+        uses: mfem/github-actions/build-metis@v2.5
+        with:
+          archive: ${{ env.METIS_ARCHIVE }}
+          dir: ${{ env.METIS_TOP_DIR }}
+
+      # MFEM build and test
+      - name: build-mfem
+        uses: mfem/github-actions/build-mfem@v2.5
+        with:
+          os: ${{ runner.os }}
+          target: opt
+          codecov: NO
+          mpi: par
+          build-system: make
+          hypre-dir: ${{ env.HYPRE_TOP_DIR }}
+          metis-dir: ${{ env.METIS_TOP_DIR }}
+          mfem-dir: mfem
+
+      - name: test (no clean)
+        run: |
+          cd mfem && make test-noclean
+
+      - name: gitignore
+        run: |
+          cd mfem/tests/scripts
+          ./runtest gitignore
@@ -13,7 +13,6 @@ name: "Checks"

 permissions:
  actions: write
-  pull-requests: read

 on:
  push:
@@ -30,11 +29,6 @@ concurrency:
 # by checking if the workflow trigger is 'push' ("github.event_name == 'push'")
 # and if we are in a fork ("github.event.pull_request.head.repo.full_name !=
 # github.repository").
-#
-# The logic for the branch-history check is slightly different, since that check
-# also inspects the PR's labels to allow for overriding failures. In this case,
-# we run on all 'pull_request' triggers, but only run for 'push' triggers that
-# do not correspond to any open PRs.

 jobs:
  file-headers-check:
@@ -134,7 +128,10 @@ jobs:

  branch-history:
    if: |
-      github.ref != 'refs/heads/next' && github.ref != 'refs/heads/master'
+      github.ref != 'refs/heads/next' &&
+      github.ref != 'refs/heads/master' &&
+      (github.event_name == 'push' ||
+       github.event.pull_request.head.repo.full_name != github.repository)
    runs-on: ubuntu-latest
    steps:
      - name: checkout mfem
@@ -142,27 +139,7 @@ jobs:
        with:
          fetch-depth: 0

-      - name: check for pull request
-        id: check_pr
-        if: github.event_name == 'push'
-        env:
-          GH_TOKEN: ${{ github.token }}
-        run: |
-          pr_exists=$(gh pr list --repo "$GITHUB_REPOSITORY" \
-            --head "$GITHUB_REF_NAME" \
-            --state open \
-            --json number \
-            --jq 'length > 0')
-          echo "pr_exists=$pr_exists" >> "$GITHUB_OUTPUT"
-
      - name: branch-history
-        id: branch_history
-        if: |
-          (github.event_name == 'pull_request' ||
-           github.event_name == 'workflow_dispatch' ||
-           steps.check_pr.outputs.pr_exists == 'false')
-        continue-on-error: ${{ contains(github.event.pull_request.labels.*.name,
-                                        'branch-history-override') }}
        run: |
          # We override origin to make sure we point to the main repo.
          # This is to have consistent test results on PRs from forks.
@@ -170,20 +147,3 @@ jobs:
          git remote add origin https://github.com/mfem/mfem.git
          git checkout -b gh-actions-branch-history
          ./config/githooks/pre-push --history
-
-      - name: report branch-history override
-        if: steps.branch_history.outcome == 'failure'
-        run: |
-          echo "::warning::branch-history check failed, but the" \
-            "'branch-history-override' label is set."
-
-  codespell:
-    runs-on: ubuntu-latest
-    steps:
-    - name: Codespell with annotations
-      uses: codespell-project/actions-codespell@master
-      with:
-        check_filenames: true
-        check_hidden: true
-        ignore_words_list: allright,ba,equil,esy,fo,hda,lod,nd,ned,numer,ot,pres,ro,seh,shat,solfes,strat,tbe,te,warmup
-
@@ -19,22 +19,18 @@ jobs:
    steps:
      - uses: actions/checkout@v6
      - uses: ./.github/actions/sanitize/config
-        with:
-          NO_FLAGS: true
      - name: Cache
        id: cache
        uses: actions/cache@v5
        with:
          path: ${{env.HYPRE_DIR}}
-          key: ${{runner.os}}-ompi-build-${{env.HYPRE_DIR}}-int32-fp64-${{ env.MFEM_ACTIONS_VERSION }}
+          key: ${{runner.os}}-ompi-build-${{env.HYPRE_DIR}}-int32-fp64-v2.5
      - name: Setup
        if: steps.cache.outputs.cache-hit != 'true'
        uses: ./.github/actions/sanitize/mpi
-        with:
-          NO_FLAGS: true
      - name: Build
        if: steps.cache.outputs.cache-hit != 'true'
-        uses: mfem/github-actions/build-hypre@v2.7
+        uses: mfem/github-actions/build-hypre@v2.5
        with:
          archive: ${{env.HYPRE_TGZ}}
          dir: ${{env.HYPRE_DIR}}
@@ -19,22 +19,18 @@ jobs:
    steps:
      - uses: actions/checkout@v6
      - uses: ./.github/actions/sanitize/config
-        with:
-          NO_FLAGS: true
      - name: Cache
        id: cache
        uses: actions/cache@v5
        with:
          path: ${{env.METIS_DIR}}
-          key: ${{runner.os}}-build-${{env.METIS_DIR}}-${{env.MFEM_ACTIONS_VERSION}}
+          key: ${{runner.os}}-build-${{env.METIS_DIR}}-v2.5
      - name: Setup
        if: steps.cache.outputs.cache-hit != 'true'
        uses: ./.github/actions/sanitize/mpi
-        with:
-          NO_FLAGS: true
      - name: Build
        if: steps.cache.outputs.cache-hit != 'true'
-        uses: mfem/github-actions/build-metis@v2.7
+        uses: mfem/github-actions/build-metis@v2.5
        with:
          archive: ${{env.METIS_TGZ}}
          dir: ${{env.METIS_DIR}}
@@ -146,7 +146,7 @@ jobs:
        if: ${{steps.restore.outputs.cache-hit != 'true'}}
        working-directory: mfem/build/tests/unit
        run: find . -type f -name '*.o' -delete
-      - uses: actions/upload-artifact@v7
+      - uses: actions/upload-artifact@v4
        with:
          name: tests-${{inputs.par}}-${{inputs.sanitizer}}
          path: mfem/build/tests/unit/${{env.unit_tests}}
@@ -172,7 +172,7 @@ jobs:
          par: ${{inputs.par}}
          sanitizer: ${{inputs.sanitizer}}
          cache-path: mfem/build/tests/unit/${{env.unit_tests}}
-      - uses: actions/download-artifact@v8
+      - uses: actions/download-artifact@v4
        if: ${{steps.restore.outputs.cache-hit != 'true'}}
        with:
          name: tests-${{inputs.par}}-${{inputs.sanitizer}}
@@ -17,17 +17,7 @@ permissions:
 on:
  push:
    branches: ["master", "next"]
-    paths-ignore: &docs-only-paths
-      - "**/*.md"
-      - "doc/**"
-      - ".binder/**"
-      - "CITATION.cff"
-      - "LICENSE"
-      - "NOTICE"
-      - "CHANGELOG"
-      - "INSTALL"
  pull_request:
-    paths-ignore: *docs-only-paths
  workflow_dispatch:

 concurrency:
@@ -260,7 +260,6 @@ miniapps/meshing/polar-nc
 miniapps/meshing/mesh-quality
 miniapps/meshing/hpref
 miniapps/meshing/phpref
-miniapps/meshing/pref321
 miniapps/meshing/mobius-strip.mesh
 miniapps/meshing/klein-bottle.mesh
 miniapps/meshing/toroid-*.mesh
@@ -444,10 +443,6 @@ miniapps/diag-smoothers/mg-abs-l1-jacobi
 miniapps/contact/contact
 miniapps/contact/ParaView

-miniapps/plasma/pic/electrostatic-*
-!miniapps/plasma/pic/electrostatic-*.cpp
-miniapps/plasma/pic/*.csv
-
 # Unit test binary and outputs
 tests/unit/output_meshes
 tests/unit/unit_tests
@@ -102,14 +102,12 @@ report_baseline:
        mkdir -p ${MACHINE_NAME}
        rundir="${MACHINE_NAME}/$(date +%Y-%m-%d)-gitlab-${BASELINE_TEST}-${CI_COMMIT_REF_SLUG}"
        rundir=$(${CI_PROJECT_DIR}/.gitlab/scripts/safe_create_rundir ${rundir})
-        status=0
-        cp ${CI_PROJECT_DIR}/${ARTIFACTS_DIR}/* ${rundir} || { status=1; }
+        cp ${CI_PROJECT_DIR}/${ARTIFACTS_DIR}/* ${rundir}
        printf "%s\n" "" "Pipeline URL:" "$CI_PIPELINE_URL" \
          >> ${rundir}/pipeline.txt
        # We create an autotest-email.html file, because that's how we signal
        # that there was an error / diff (temporary).
-        if [[ $status -ne 0 ]] || \
-           [[ -f ${rundir}/${BASELINE_TEST}.err ]] || \
+        if [[ -f ${rundir}/${BASELINE_TEST}.err ]] || \
           [[ -f ${rundir}/${BASELINE_TEST}-${MACHINE_NAME}.diff ]]; then
          cp ${rundir}/pipeline.txt ${rundir}/autotest-email.html
        fi
@@ -85,8 +85,3 @@ opt_par_gcc_10_pumi:
  extends: .mfem_job_on_dane
  variables:
    SPEC: "%gcc@10.3.1 +pumi"
-
-opt_par_gcc_10_gslib:
-  extends: .mfem_job_on_dane
-  variables:
-    SPEC: "%gcc@10.3.1 +gslib"
@@ -63,8 +63,3 @@ opt_mpi_cuda_hypre_cuda_gcc:
  extends: .mfem_job_on_matrix
  variables:
    SPEC: "%gcc@10.3.1 +mpi +cuda cuda_arch=90 ^hypre+cuda"
-
-opt_mpi_cuda_gcc_gslib:
-  extends: .mfem_job_on_matrix
-  variables:
-    SPEC: "%gcc@10.3.1 +mpi +cuda +gslib cuda_arch=90 ^hypre+cuda"
@@ -32,9 +32,9 @@ mkdir _${BASELINE_TEST} && cd _${BASELINE_TEST}

 # run
 if [[ "${MACHINE_NAME}" == "dane" ]]; then
-  srun --nodes=1 -t 60 --exclusive --reservation=ci ../runtest ../../mfem "${BASELINE_TEST} ${TPLS_DIR}"
+  salloc --nodes=1 -t 60 --exclusive --reservation=ci ../runtest ../../mfem "${BASELINE_TEST} ${TPLS_DIR}"
 elif [[ ${MACHINE_NAME} == "corona" ]]; then
-  srun --nodes=1 -t 60 -p pbatch ../runtest ../../mfem "${BASELINE_TEST} ${TPLS_DIR}"
+  salloc --nodes=1 -t 60 -p pbatch ../runtest ../../mfem "${BASELINE_TEST} ${TPLS_DIR}"
 else
  echo "Unknown machine: MACHINE_NAME=$MACHINE_NAME"
  exit 1
@@ -8,115 +8,47 @@
                               https://mfem.org


-Version 4.9.1 (development)
-===========================
-
- Added policy for AI-assisted contribution to CONTRIBUTING.md.
+Version 4.10 (development)
+==========================

 Discretization improvements
 ---------------------------
- Improved FindPointsGSLIB surface mesh capability with support for simplices
-  and an option to specify axis-aligned bounding box padding for near-surface
-  point queries.
+- 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.

- Added GPU-enabled partial assembly for simplicial Bernstein H1 basis based on
-  ragged tensor algorithms (see DOI: 10.1137/11082539X) for mass and diffusion
-  integrators.
+Version 4.9.1 (development)
+===========================

- 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 and Vincent, DOI: 10.1016/j.camwa.2015.03.017
-  * Tet rules (d=1-13): Witherden and  Vincent (same as above)
-  * Tet rules (d=14-20): Chuluunbaatar et al., DOI: 10.1016/j.camwa.2022.08.016
-
- Added support for general 1D Gauss-Jacobi quadrature rules and Stroud conical
-  quadrature rules on triangles and tetrahedra.
-
- Improved the GridFunction projection routines. Projections work for Scalar,
+Discretization improvements
+---------------------------
+- Improved the gridfunction projection routines. Projections work for Scalar,
  Vector and VectorFE, also NURBS versions. Optionally different types of
-  projections can be selected, default behavior has not changed.
+  projections can be selected, default behaviour has not changed.

- Added GridFunction projection methods for trace spaces, i.e., project
-  coefficients on the mesh skeleton.
-
- Added methods to estimate function extremum using piecewise linear bounds plus
+- Added methods to estimate function extremum using piecewise linear bounds +
  recursive subdivision.

- Extend FindPointsGSLIB to support surface meshes.
-
- Added support for complex-valued mixed bilinear forms via the new classes
-  MixedSesquilinearForm and ParMixedSesquilinearForm, mirroring the existing
-  SesquilinearForm classes. Rectangular complex operators are now also
-  handled correctly by ComplexSparseMatrix::GetSystemMatrix and
-  ComplexHypreParMatrix::GetSystemMatrix, which previously assumed equal
-  trial and test spaces.
-
 Meshing improvements
 --------------------
- Added support for nonuniform anisotropic mesh refinement on parallel quad/hex
-  meshes with arbitrary spacing in each direction. This enables in particular
-  3:1 refinement in parallel, as demonstrated in the new meshing miniapp pref321.
-
- Added option to guarantee mesh validity during TMOP-based r-adaptivity, using
-  bounds on the determinant of the mesh transformation Jacobian.
-
- Added PA support for TMOP's adaptive limiting functionality. Multiple
-  GridFunctions and Coefficients can be combined to form a composite term.
-
 - Improved support for 1D NURBS meshes with variable order, including using
  the patches construct for 1D NURBS meshes.

- Added the option to include material interfaces (faces separating elements
-  with different element attributes) as additional boundary elements, for
-  parallel visualization, e.g. with GLVis. This is supported by both the Print
-  and PrintAsOne methods of ParMesh. See ParMesh::SetPrintInterfaces().
-
-Linear and nonlinear solvers
----------------------------
- Added support for trace spaces in PRefinementTransferOperator. This is used in
-  PRefinement multigrid methods for problems posed on trace spaces (see e.g. the
-  DPG miniapps).
-
-GPU computing
-------------
- Improved partial assembly for VectorDivergenceIntegrator with shared-memory
-  kernels, kernel registration, and transpose support.
-
- Improved partial-assembly diagonal kernels for VectorMassIntegrator (shared-
-  memory specializations) and ElasticityIntegrator (no scratch Q-vector).
-
- Added PA gradient and diagonal support for VectorConvectionNLFIntegrator
- (AssembleGradPA, AddMultGradPA, AssembleGradDiagonalPA).
-
- Added device assembly support for 3D H(curl) VectorFEDomainLFIntegrator.
-
- Added NVIDIA cuDSS library interface. Implementation examples have been
-  added to ex1 and ex1p. See https://developer.nvidia.com/cudss for more
-  details. Supported versions >= 0.6.0.
-
- Allow specifying GPU kernel launch bounds for native and RAJA GPU backends.
-
- Changed VectorFEMassIntegrator to use kernel specialization dispatch for
-  partial assembly.
-
- Added support for FiniteElement::MapType::INTEGRAL spaces to
-  QuadratureInterpolator.
-
 New and updated examples and miniapps
 -------------------------------------
- The Lorentz miniapp (in miniapps/electromagnetics) has been updated to
-  leverage the ParticleSet capability.
-
- Added (Complex)PRefinementMultigrid solver option in the DPG miniapps.
+- Electromagnetics/lorentz miniapp has been updated to leverage the ParticleSet
+  capability.

 Miscellaneous
 -------------
- Fixed signed DOF handling in ParGridFunction reading (read constructor) and
-  saving via SaveAsOne(). Simplified the process of applying the DOF signs by
-  using the new method ApplyDofSigns() in class ParFiniteElementSpace: the
-  method will return immediately if no sign flips are needed.
+- 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,12 +75,10 @@ set(XSDK_ENABLE_Fortran OFF)

 # Check if we need to enable C or Fortran.
 if (MFEM_USE_CONDUIT OR
-    MFEM_USE_SIDRE OR
    MFEM_USE_PETSC)
   # This seems to be needed by:
   #  * find_package(BLAS REQUIRED) and
-   #  * find_package(HDF5 REQUIRED) needed, in turn, by:
-   #    - find_package(AXOM REQUIRED)
+   #  * find_package(HDF5 REQUIRED) and
   #  * find_package(PETSc REQUIRED)
   set(XSDK_ENABLE_C ON)
 endif()
@@ -88,9 +86,18 @@ if (MFEM_USE_STRUMPACK OR MFEM_USE_MUMPS)
  # Just needed to find the MPI_Fortran libraries to link with
  set(XSDK_ENABLE_Fortran ON)
 endif()
-# RAJA requires C++20:
-if (MFEM_USE_RAJA AND ("${CMAKE_CXX_STANDARD}" LESS "20"))
-  set(CMAKE_CXX_STANDARD 20 CACHE STRING "C++ standard to use." FORCE)
+# Ginkgo requires C++17:
+if ((MFEM_USE_GINKGO) AND ("${CMAKE_CXX_STANDARD}" LESS "17"))
+  set(CMAKE_CXX_STANDARD 17 CACHE STRING "C++ standard to use." FORCE)
+# Google Benchmark, SUNDIALS, STRUMPACK, Tribol, RAJA and Umpire require C++14:
+elseif ((MFEM_USE_BENCHMARK OR
+         MFEM_USE_SUNDIALS OR
+         MFEM_USE_STRUMPACK OR
+         MFEM_USE_TRIBOL OR
+         MFEM_USE_RAJA OR
+         MFEM_USE_UMPIRE) AND
+    ("${CMAKE_CXX_STANDARD}" LESS "14"))
+  set(CMAKE_CXX_STANDARD 14 CACHE STRING "C++ standard to use." FORCE)
 endif()

 # Include xSDK default CMake file.
@@ -230,13 +237,6 @@ else()
  set(MFEM_DEBUG OFF)
 endif()

-# Shadow warnings for clang only; GCC's -Wshadow flags more.
-if (CMAKE_CXX_COMPILER_ID MATCHES "Clang")
-  set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -pedantic -Wall -Wshadow")
-elseif (CMAKE_CXX_COMPILER_ID STREQUAL "GNU")
-  set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} -pedantic -Wall")
-endif()
-
 # Shared build on Windows
 if (WIN32 AND BUILD_SHARED_LIBS)
  # CMAKE_WINDOWS_EXPORT_ALL_SYMBOLS works only with MSVC?
@@ -431,15 +431,6 @@ if (MFEM_USE_STRUMPACK)
  endif()
 endif()

-# cuDSS can only be enabled in CUDA
-if (MFEM_USE_CUDSS)
-  if (MFEM_USE_CUDA)
-    find_package(CUDSS REQUIRED)
-  else()
-    message(FATAL_ERROR " *** cuDSS requires that CUDA be enabled.")
-  endif()
-endif()
-
 # GnuTLS
 if (MFEM_USE_GNUTLS)
  find_package(_GnuTLS REQUIRED)
@@ -485,11 +476,6 @@ if (MFEM_USE_FMS)
    find_package(FMS REQUIRED fms)
 endif()

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

@@ -3,12 +3,12 @@
 </p>

 <p align="center">
-<a href="https://github.com/mfem/mfem/blob/master/LICENSE"><img alt="License" src="https://img.shields.io/badge/License-BSD-blue.svg"></a>
-<a href="https://github.com/mfem/mfem/releases/latest"><img alt="GitHub release" src="https://img.shields.io/github/v/release/mfem/mfem"></a>
-<a href="https://github.com/mfem/mfem/actions/workflows/repo-check.yml?query=branch%3Amaster"><img alt="Repo check" src="https://github.com/mfem/mfem/actions/workflows/repo-check.yml/badge.svg?branch=master"></a>
-<a href="https://github.com/mfem/mfem/actions/workflows/builds-and-tests.yml?query=branch%3Amaster"><img alt="Builds and Tests" src="https://github.com/mfem/mfem/actions/workflows/builds-and-tests.yml/badge.svg?branch=master"></a>
+<a href="https://github.com/mfem/mfem/blob/master/LICENSE"><img alt="License" src="https://img.shields.io/badge/License-BSD-brightgreen.svg"></a>
+<a href="https://github.com/mfem/mfem/actions?query=workflow%3Arepo-check+branch%3Amaster"><img alt="Repo check" src="https://github.com/mfem/mfem/actions/workflows/repo-check.yml/badge.svg?branch=master"></a>
+<a href="https://github.com/mfem/mfem/actions?query=workflow%3Abuild-analysis+branch%3Amaster"><img alt="Build Analysis" src="https://github.com/mfem/mfem/actions/workflows/mfem-analysis.yml/badge.svg?branch=master"></a>
+<a href="https://github.com/mfem/mfem/actions?query=workflow%3Abuilds-and-tests+branch%3Amaster"><img alt="Builds and Tests" src="https://github.com/mfem/mfem/actions/workflows/builds-and-tests.yml/badge.svg?branch=master"></a>
 <a href="https://ci.appveyor.com/project/mfem/mfem"><img alt="Build Status" src="https://ci.appveyor.com/api/projects/status/19non9sqm6msi2wy?svg=true"></a>
-<a href="https://docs.mfem.org/html/index.html"><img alt="Documentation" src="https://img.shields.io/badge/code-documented-brightgreen.svg"></a>
+<a href="https://docs.mfem.org/html/index.html"><img alt="Doxygen" src="https://img.shields.io/badge/code-documented-brightgreen.svg"></a>
 </p>


@@ -24,14 +24,6 @@ must be made under this license.
 Note also that MFEM has a [Code of Conduct](CODE_OF_CONDUCT.md). By participating
 in the MFEM community, you agree to abide by its rules.

-## AI Policy
- Use of AI code generation in MFEM is allowed but must be disclosed, e.g. by
-  selecting the `AI-assisted` label on the PR.
- By submitting a PR, the author acknowledges that they have reviewed and
-  understand the changes they are proposing.
- PR authors are still responsible for correctness, licensing, and attribution
-  of all changes.
-
 If you plan on contributing to MFEM, consider reviewing the
 [issue tracker](https://github.com/mfem/mfem/issues) first to check if a thread
 already exists for your desired feature or the bug you ran into. Use a pull
@@ -84,7 +76,7 @@ Origin](#developers-certificate-of-origin-11) at the end of this file.*
  follow the [MFEM PR Rules](#mfem-pr-rules).
 - When your contribution is fully working and ready to be reviewed, add
  the `ready-for-review` label.
- PRs are treated similarly to journal submission, with an "editor" assigning two
+- PRs are treated similarly to journal submission with an "editor" assigning two
  reviewers to evaluate the changes.
 - The reviewers have 3 weeks to evaluate the PR and work with the author to
  fix issues and implement improvements.
@@ -125,7 +117,7 @@ The MFEM source code has the following structure:
  │   ├── petsc
  │   ├── pumi
  │   ├── sundials
-  │   └── superlu
+  |   └── superlu
  ├── fem
  │   ├── ceed
  │   ├── dfem
@@ -137,6 +129,10 @@ The MFEM source code has the following structure:
  │   ├── moonolith
  │   ├── qinterp
  │   └── tmop
+  │   |   ├── assemble
+  │   |   ├── metrics
+  │   |   ├── mult
+  │   |   └── tools
  ├── general
  ├── linalg
  │   ├── batched
@@ -149,10 +145,11 @@ The MFEM source code has the following structure:
  │   ├── common
  │   ├── contact
  │   ├── dfem
-  │   ├── diag-smoothers
  │   ├── dpg
  │   ├── electromagnetics
  │   ├── fluids
+  │   │   ├── navier
+  │   │   └── schrodinger-flow
  │   ├── gslib
  │   ├── hdiv-linear-solver
  │   ├── hooke
@@ -162,7 +159,6 @@ The MFEM source code has the following structure:
  │   ├── nurbs
  │   ├── parelag
  │   ├── performance
-  │   ├── plasma
  │   ├── shifted
  │   ├── solvers
  │   ├── spde
@@ -193,15 +189,15 @@ respectively.

 - The main finite element classes are:
  + [`FiniteElement`](https://docs.mfem.org/html/classmfem_1_1FiniteElement.html)
-  + [`FiniteElementCollection`](https://docs.mfem.org/html/classmfem_1_1FiniteElementCollection.html)
+  + [`FiniteElementCollection`](https://docs.mfem.org/html/classmfem_1_1FiniteElement.html)
  + [`FiniteElementSpace`](https://docs.mfem.org/html/classmfem_1_1FiniteElementSpace.html)
  + [`GridFunction`](https://docs.mfem.org/html/classmfem_1_1GridFunction.html)
  + [`BilinearFormIntegrator`](https://docs.mfem.org/html/classmfem_1_1BilinearFormIntegrator.html) and [`LinearFormIntegrator`](https://docs.mfem.org/html/classmfem_1_1LinearFormIntegrator.html)
-  + [`LinearForm`](https://docs.mfem.org/html/classmfem_1_1LinearForm.html), [`BilinearForm`](https://docs.mfem.org/html/classmfem_1_1BilinearForm.html) and [`MixedBilinearForm`](https://docs.mfem.org/html/classmfem_1_1MixedBilinearForm.html)
+  + [`LinearForm`](https://docs.mfem.org/html/classmfem_1_1LinearFormIntegrator.html), [`BilinearForm`](https://docs.mfem.org/html/classmfem_1_1BilinearForm.html) and [`MixedBilinearForm`](https://docs.mfem.org/html/classmfem_1_1MixedBilinearForm.html)

 - The main linear algebra classes and sources are
  + [`Operator`](https://docs.mfem.org/html/classmfem_1_1Operator.html) and [`BilinearForm`](https://docs.mfem.org/html/classmfem_1_1BilinearForm.html)
-  + [`Vector`](https://docs.mfem.org/html/classmfem_1_1Vector.html) and [`LinearForm`](https://docs.mfem.org/html/classmfem_1_1LinearForm.html)
+  + [`Vector`](https://docs.mfem.org/html/classmfem_1_1BilinearForm.html) and [`LinearForm`](https://docs.mfem.org/html/classmfem_1_1LinearForm.html)
  + [`DenseMatrix`](https://docs.mfem.org/html/classmfem_1_1DenseMatrix.html) and [`SparseMatrix`](https://docs.mfem.org/html/classmfem_1_1SparseMatrix.html)
  + Sparse [smoothers](https://docs.mfem.org/html/sparsesmoothers_8hpp.html) and linear [solvers](https://docs.mfem.org/html/solvers_8hpp.html)

@@ -213,8 +209,8 @@ shared geometric entities between different tasks. The parallel source files
 have a `p` prefix, e.g. `pmesh.cpp` vs. the serial `mesh.cpp`.

 - The main parallel classes are
-  + [`ParMesh`](https://docs.mfem.org/html/classmfem_1_1ParMesh.html)
-  + [`ParNCMesh`](https://docs.mfem.org/html/classmfem_1_1ParNCMesh.html)
+  + [`ParMesh`](https://docs.mfem.org/html/solvers_8hpp.html)
+  + [`ParNCMesh`](https://docs.mfem.org/html/classmfem_1_1ParMesh.html)
  + [`ParFiniteElementSpace`](https://docs.mfem.org/html/classmfem_1_1ParFiniteElementSpace.html)
  + [`ParGridFunction`](https://docs.mfem.org/html/classmfem_1_1ParGridFunction.html)
  + [`ParBilinearForm`](https://docs.mfem.org/html/classmfem_1_1ParBilinearForm.html) and [`ParLinearForm`](https://docs.mfem.org/html/classmfem_1_1ParLinearForm.html)
@@ -224,14 +220,14 @@ have a `p` prefix, e.g. `pmesh.cpp` vs. the serial `mesh.cpp`.
 #### GPU and general device support

 GPU and multi-core CPU support is based on device kernels supporting different
-backends (CUDA, HIP, OCCA, RAJA, OpenMP, etc.) and an internal lightweight
+backends (CUDA, OCCA, RAJA, OpenMP, etc.) and an internal lightweight
 device/host memory manager.

 - The main device-relevant classes and sources are:
  + [`Device`](https://docs.mfem.org/html/device_8hpp.html)
  + [`MemoryManager`](https://docs.mfem.org/html/mem_manager_8hpp.html)
  + the [`mfem::forall`](https://docs.mfem.org/html/forall_8hpp.html) function
-  + the [`cuda.hpp`](https://docs.mfem.org/html/cuda_8hpp.html), [`hip.hpp`](https://docs.mfem.org/html/hip_8hpp.html) and [`occa.hpp`](https://docs.mfem.org/html/occa_8hpp.html) files
+  + the [`cuda.hpp`](https://docs.mfem.org/html/cuda_8hpp.html) and [`occa.hpp`](https://docs.mfem.org/html/occa_8hpp.html) files

 #### Utilities, building and documentation
 - The `general/` directory contains C++ classes that serve as utilities for
@@ -245,8 +241,8 @@ device/host memory manager.
 - `examples` and `miniapps` respectively gather simple and more fully-featured
  demonstrations of the usage on MFEM. They both rely on `data/` for the
  collection of meshes.
- The `tests/` directory contains a unit test suite, additional tests, and
-  benchmarks.
+- The `tests/` directory contains a unit test suite and will later contain more
+  tests that run example codes.

 See also the [code overview](https://mfem.org/code-overview/) section on the MFEM
 website.
@@ -280,8 +276,8 @@ Before you can start, you need a GitHub account, here are a few suggestions:
  the top of https://github.com/mfem.
 - Consider making your membership public by going to https://github.com/orgs/mfem/people
  and clicking on the organization visibility drop box next to your name.
- Project discussions and announcements will be posted at https://github.com/orgs/mfem/discussions,
-  tagging the `@mfem/everyone` team when appropriate.
+- Project discussions and announcements will be posted at
+  https://github.com/orgs/mfem/teams/everyone.

 #### Structure
 - The MFEM source code is in the [mfem](https://github.com/mfem/mfem)
@@ -341,12 +337,11 @@ Before you can start, you need a GitHub account, here are a few suggestions:
  - Well-designed simple code is frequently more general and powerful.
  - Lean code base is easier to understand by new collaborators.
  - New features should be added only if they are necessary or generally useful.
-  - Introduction of language constructs not currently used in MFEM should be
+  - Introduction of language constructions not currently used in MFEM should be
    justified and generally avoided (to maintain portability to various systems
    and compilers, including early access hardware).
-  - We prefer basic C++. Use C++17 features judiciously, prioritizing readability,
-    consistency with existing MFEM code, and portability to different systems,
-    compilers and device backends.
+  - We prefer basic C++ and the C++03 standard, to keep the code readable by
+    a large audience and to make sure it compiles anywhere.

 - *Keep the code general and reasonably efficient*
  - The main goal is fast prototyping for research and application development.
@@ -389,7 +384,7 @@ Before you can start, you need a GitHub account, here are a few suggestions:
 - When your branch is ready for other developers to review / comment on
  the code, create a pull request towards `mfem:master`.

- Pull requests typically have titles like:
+- Pull request typically have titles like:

     `Description [new-feature-dev]`

@@ -410,12 +405,12 @@ Before you can start, you need a GitHub account, here are a few suggestions:
 - Add a description, appropriate labels and assign yourself to the PR. The MFEM
  team will add reviewers as appropriate.

- List outstanding TODO items in the description.
+- List outstanding TODO items in the description, see PR #222 for an example.

 - When your contribution is fully working and ready to be reviewed, add
-  or request the `ready-for-review` label.
+  the `ready-for-review` label.

- PRs are treated similarly to journal submission, with an "editor" assigning
+- PRs are treated similarly to journal submission with an "editor" assigning
  two reviewers to evaluate the changes. The reviewers have 3 weeks to evaluate
  the PR and work with the author to implement improvements and fix issues.

@@ -441,7 +436,7 @@ Before you can start, you need a GitHub account, here are a few suggestions:
  checks in GitHub Actions enforce MFEM-specific rules which are explained in
  the error messages and the `tests/scripts` directory.

- Also note that the tests `branch-history` and `repo-check` found in GitHub
+- Also note that the tests `branch-history` and `repos-checks` found in GitHub
  Actions can be triggered automatically before each push using git hooks. See
  the [git hooks README](config/githooks/README.md) for a detailed explanation.

@@ -498,15 +493,15 @@ Everyone on the MFEM team can be asked to serve as a reviewer on a PR in their a

 3. To ensure the quality of the PR by making sure that the code adheres to the [Developer Guidelines](#developer-guidelines), e.g. all methods, data members, and functions have documentation, including data ownership and lifetime, new examples/miniapps have a corresponding PR in mfem/web, major features have `CHANGELOG` entries, etc.

-4. To seek help from the editors in case of difficulties.
+3. To seek help from the editors in case of difficulties.

-5. To complete the review in a timely manner: 3 weeks from assignment.
+4. To complete the review in a timely manner: 3 weeks from assignment.

-6. To test the PR thoroughly before merging in *next*. The PR author is also encouraged to perform testing and inform the reviewers about the results.
+5. To test the PR thoroughly before merging in *next*. The PR author is also encouraged to perform testing and inform the reviewers about the results.

-7. To monitor the PR impact on the testing in the *next* branch and alert the editors that the PR is ready for merging in *master*.
+6. To monitor the PR impact on the testing in the *next* branch and alert the editors that the PR is ready for merging in *master*.

-8. The review of bugfixes should be expedited proportional to their importance. The review window can be much less than three weeks in such cases.
+7. The review of bugfixes should be expedited proportional to their importance. The review window can be much less than three weeks in such cases.

 #### Responsibilities of Authors

@@ -532,30 +527,30 @@ Before a PR can be merged, it should satisfy the following:
 - [ ] Code builds.
 - [ ] Code passes `make style`.
 - [ ] Update `CHANGELOG`:
-   - [ ] Is this a new feature users need to be aware of? New or updated example or miniapp?
-   - [ ] Does it make sense to create a new section in the `CHANGELOG` to group with other related features?
+    - [ ] Is this a new feature users need to be aware of? New or updated example or miniapp?
+    - [ ] Does it make sense to create a new section in the `CHANGELOG` to group with other related features?
 - [ ] Update `INSTALL`:
-   - [ ] Has a new optional library been added? If so, what range of versions of this library are required? (*Make sure the external library is compatible with our BSD license, e.g. it is not licensed under GPL!*)
-   - [ ] Have the version ranges for any required or optional libraries changed?
-   - [ ] Does `make` or `cmake` have a new target?
-   - [ ] Did the requirements or the installation process change? *(rare)*
+    - [ ] Had a new optional library been added? If so, what range of versions of this library are required? (*Make sure the external library is compatible with our BSD license, e.g. it is not licensed under GPL!*)
+    - [ ] Have the version ranges for any required or optional libraries changed?
+    - [ ] Does `make` or `cmake` have a new target?
+    - [ ] Did the requirements or the installation process change? *(rare)*
 - [ ] Update continuous integration server configurations if necessary (e.g. with new version requirements for each of MFEM's dependencies)
-   - [ ] `.github`
-   - [ ] `.appveyor.yml`
+    - [ ] `.github`
+    - [ ] `.appveyor.yml`
 - [ ] Update `.gitignore`:
-   - [ ] Check if `make distclean; git status` shows any files that were generated from the source by the project (not an IDE) but we don't want to track in the repository.
-   - [ ] Add new patterns (just for the new files above) and re-run the above test.
+    - [ ] Check if `make distclean; git status` shows any files that were generated from the source by the project (not an IDE) but we don't want to track in the repository.
+    - [ ] Add new patterns (just for the new files above) and re-run the above test.
 - [ ] New examples:
-   - [ ] All sample runs at the top of the example source file work.
-   - [ ] Update `examples/makefile`:
+    - [ ] All sample runs at the top of the example source file work.
+    - [ ] Update `examples/makefile`:
      - [ ] Add the example code to the appropriate `SEQ_EXAMPLES` and `PAR_EXAMPLES` variables.
      - [ ] Add any files generated by it to the `clean` target.
      - [ ] Add the example binary and any files generated by it to the top-level `.gitignore` file.
-   - [ ] Update `examples/CMakeLists.txt`:
+    - [ ] Update `examples/CMakeLists.txt`:
      - [ ] Add the example code to the `ALL_EXE_SRCS` variable.
      - [ ] Make sure `THIS_TEST_OPTIONS` is set correctly for the new example.
   - [ ] List the new example in `doc/CodeDocumentation.dox`.
-   - [ ] If new examples directory (e.g. `examples/pumi`), list it in `doc/CodeDocumentation.conf.in`
+   - [ ] If new examples directory (e.g.`examples/pumi`), list it in `doc/CodeDocumentation.conf.in`
   - [ ] Companion pull request for documentation in [mfem/web](https://github.com/mfem/web) repo:
      - [ ] Update or add example-specific documentation, see e.g. the `src/examples.md`.
      - [ ] Add the description, labels and screenshots in `src/examples.md` and `src/img`.
@@ -572,13 +567,13 @@ Before a PR can be merged, it should satisfy the following:
      - [ ] Add/update the `CMakeLists.txt` file in the new miniapp directory.
      - [ ] Consider adding a new test for the new miniapp.
   - [ ] List the new miniapp in `doc/CodeDocumentation.dox`
-   - [ ] If new miniapps directory (e.g. `miniapps/nurbs`), add it to `MINIAPP_SUBDIRS` in the `makefile`.
-   - [ ] If new miniapps directory (e.g. `miniapps/nurbs`), list it in `doc/CodeDocumentation.conf.in`
+   - [ ] If new miniapps directory (e.g.`miniapps/nurbs`), add it to `MINIAPP_SUBDIRS` in the `makefile`.
+   - [ ] If new miniapps directory (e.g.`miniapps/nurbs`), list it in `doc/CodeDocumentation.conf.in`
   - [ ] Companion pull request for documentation in [mfem/web](https://github.com/mfem/web) repo:
-      - [ ] Update or add miniapp-specific documentation, see e.g. the `src/meshing.md` and `src/electromagnetics.md` files.
-      - [ ] Add the description, labels and screenshots in `src/examples.md` and `src/img`.
-      - [ ] The miniapps go at the end of the page, and are usually listed only under a specific "Application (PDE)" category.
-      - [ ] Add a short description of the miniapp in the "Extensive Examples" section of `features.md`.
+     - [ ] Update or add miniapp-specific documentation, see e.g. the `src/meshing.md` and `src/electromagnetics.md` files.
+     - [ ] Add the description, labels and screenshots in `src/examples.md` and `src/img`.
+     - [ ] The miniapps go at the end of the page, and are usually listed only under a specific "Application (PDE)" category.
+     - [ ] Add a short description of the miniapp in the "Extensive Examples" section of `features.md`.
 - [ ] New capability:
   - [ ] All new public, protected, and private classes, methods, data members, and functions have full Doxygen-style documentation in source comments. Documentation should include descriptions of member data, function arguments and return values, template parameters, and prerequisites for calling new functions.
   - [ ] Pointer arguments and return values must specify whether ownership is being transferred or lent with the call.
@@ -680,7 +675,7 @@ MFEM uses a `master`/`next`-branch workflow as described below:
 - [ ] Update URL shortlinks:
    - [ ] Create a shortlink at [http://bit.ly/](http://bit.ly/) for the release tarball, e.g. https://mfem.github.io/releases/mfem-3.1.tgz.
    - [ ] (LLNL only) Add and commit the new shortlink in the `links` and `links-mfem` files of the internal `mfem/downloads` repo.
-    - [ ] Add the new shortlinks to the MFEM package in `spack`.
+    - [ ] Add the new shortlinks to the MFEM packages in `spack`, `homebrew/science`, `VisIt`, etc.
 - [ ] Update website in `mfem/web` repo:
    - Update version and shortlinks in `src/index.md` and `src/download.md`.
    -  Use [cloc-1.62.pl](http://cloc.sourceforge.net/) and `ls -lh` to estimate the SLOC and the tarball size in `src/download.md`.
@@ -732,24 +727,22 @@ commit or push, see the [README](config/githooks/README.md) in the `config/githo
 directory.


-### GitHub Actions smoke tests
-
+### Linux and Mac smoke tests
 We use GitHub Actions to drive the default tests on the `master` and `next`
 branches. See the `.github/workflows` files and the logs at
 [https://github.com/mfem/mfem/actions](https://github.com/mfem/mfem/actions).

-GitHub Actions testing should be kept lightweight, as there is a time
-constraint on jobs. The current workflows cover Linux, macOS, and Windows
-configurations.
+Testing using GitHub Actions should be kept lightweight, as there is a time
+constraint on jobs. Two virtual machines are configured - Mac (OS X) and Linux.

 - Tests on the `master` branch are triggered whenever a PR is issued on this branch.
 - Tests on the `next` branch are currently scheduled to run each night.

-### Additional Windows smoke test

-We also use Appveyor to test building with the MS Visual C++ compiler in a Windows
-environment, as well as to test the CMake build. See the `.appveyor.yml` file
-and the build logs at
+### Windows smoke test
+We use Appveyor to test building with the MS Visual C++ compiler in a Windows
+environment, as well as to test the CMake build. See the `.appveyor` file and the
+build logs at
 [https://ci.appveyor.com/project/mfem/mfem](https://ci.appveyor.com/project/mfem/mfem).

 CMake is used to generate the MSVC Project files and drive the build.  A release
@@ -38,13 +38,14 @@ the option MFEM_USE_METIS.
 MFEM also includes support for devices such as GPUs, and programming models such
 as CUDA, HIP, OCCA, OpenMP and RAJA.

- Starting with version 4.9, MFEM requires a C++17 compiler.
+- Starting with version 4.0, MFEM requires a C++11 compiler. We recommend using
+  a newer compiler, e.g. GCC version 4.9 or higher.

 - CUDA support requires an NVIDIA GPU and an installation of the CUDA Toolkit
  https://developer.nvidia.com/cuda-toolkit

 - HIP support requires an AMD GPU and an installation of the ROCm software stack
-  https://rocm.docs.amd.com
+  https://rocmdocs.amd.com

 - OCCA support requires the OCCA library
  https://libocca.org
@@ -82,9 +83,9 @@ Serial build:
 Parallel build:
   (download hypre and METIS 4 from above URLs)
   (build METIS 4 in ../metis-4.0 relative to mfem/)
-   (for METIS 5, see https://mfem.org/building/#parallel-build-using-metis-5)
   (build hypre in ../hypre relative to mfem/)
   make parallel -j 4
+(For METIS 5, see https://mfem.org/building/#parallel-build-using-metis-5)

 CUDA build:
   make cuda -j 4
@@ -114,14 +115,14 @@ Serial build:
 Parallel build:
   (download hypre and METIS 4 from above URLs)
   (build METIS 4 in ../metis-4.0 relative to mfem/)
-   (for METIS 5, see https://mfem.org/building/#parallel-build-using-metis-5)
   (build hypre in ../hypre relative to mfem/)
   mkdir <mfem-build-dir> ; cd <mfem-build-dir>
   cmake <mfem-source-dir> -DMFEM_USE_MPI=YES
   make -j 4
+(For METIS 5, see https://mfem.org/building/#parallel-build-using-metis-5)

 Parallel build with fetching of hypre and METIS:
-   mkdir <mfem-build-dir> ; cd <mfem-build-dir>
+   mkdir <mfem-buil-dir> ; cd <mfem-build-dir>
   cmake <mfem-source-dir> -DMFEM_USE_MPI=YES -DMFEM_FETCH_TPLS=YES
   make -j 4

@@ -133,8 +134,7 @@ CUDA build:

 HIP build:
   mkdir <mfem-build-dir> ; cd <mfem-build-dir>
-   cmake <mfem-source-dir> -DMFEM_USE_HIP=YES -DHIP_ARCH=gfx942 \
-         -DCMAKE_CXX_COMPILER=amdclang++ -DCMAKE_HIP_COMPILER=amdclang++
+   cmake <mfem-source-dir> -DMFEM_USE_HIP=YES -DHIP_ARCH=gfx942 -DCMAKE_CXX_COMPILER=amdclang++ -DCMAKE_HIP_COMPILER=amdclang++
   make -j 4

 Example codes (serial/parallel, depending on the build):
@@ -269,7 +269,6 @@ Compilers:
   CXX      - C++ compiler, serial build
   MPICXX   - MPI C++ compiler, parallel build
   CUDA_CXX - The CUDA compiler, 'nvcc' or 'clang++'
-   HIP_CXX  - The HIP compiler, e.g. 'hipcc'

 Compiler options:
   OPTIM_FLAGS - Options for optimized build
@@ -396,11 +395,6 @@ MFEM_USE_STRUMPACK = YES/NO
   classes. When enabled, this option uses the STRUMPACK_* library options, see
   below.

-MFEM_USE_CUDSS = YES/NO
-   Enable MFEM functionality based on the cuDSS library. When using cuDSS, CUDA
-   support must be also enabled in MFEM, i.e. MFEM_USE_CUDA=YES must be set.
-   When enabled, this option uses the CUDSS_* library options, see below.
-
 MFEM_USE_GINKGO = YES/NO
   Enable MFEM functionality based on the Ginkgo library, which provides
   iterative linear solvers and preconditioners with OpenMP, CUDA backends, see
@@ -458,13 +452,6 @@ MFEM_USE_MPFR = YES/NO
   quadrature rules. When enabled, this option uses the MPFR_* library options,
   see below.

-MFEM_USE_SIDRE = YES/NO
-   Sidre is a component of LLNL's axom project, https://github.com/LLNL/axom,
-   that provides an HDF5-based file format for visualization or restart
-   capability following the Conduit (https://github.com/LLNL/conduit) mesh
-   blueprint specification. When enabled, this option requires installation of
-   HDF5 (see also MFEM_USE_NETCDF), Conduit and LLNL's axom project.
-
 MFEM_USE_SIMD = YES/NO
   Enables the high performance templated classes to use architecture dependent
   SIMD intrinsics instead of the generic implementation of class AutoSIMD in
@@ -560,13 +547,13 @@ MFEM_USE_RAJA = YES/NO
 MFEM_USE_OCCA = YES/NO
   Enables support for the OCCA library in MFEM. OCCA is an open-source library
   which aims to make it easy to program different types of devices (e.g. CPU,
-   GPU, FPGA) by providing a unified API for interacting with JIT-compiled
+   GPU, FPGA) by providing an unified API for interacting with JIT-compiled
   backends. In order to use the OCCA CUDA backend, CUDA support must be enabled
   in MFEM as well, i.e. MFEM_USE_CUDA=YES must be set.

 MFEM_USE_GSLIB = YES/NO
   Enables MFEM functionality based on the GSLIB library, and specifically its
-   FindPoints component, which provides robust algorithms to evaluate finite
+   FindPoints component, which provides a robust algorithms to evaluate finite
   element functions in a collection of points in physical space. When enabled,
   the user can use the GSLIB-FindPoints methods as shown in miniapps/gslib.

@@ -725,18 +712,9 @@ The specific libraries and their options are:
  Options: STRUMPACK_OPT, STRUMPACK_LIB.
  Versions: STRUMPACK >= 3.0.0.

- CUDSS (optional), used when MFEM_USE_CUDSS = YES. Note that CUDSS requires 
-  CUDA 12.x toolkit and the cuDSS libraries. The supported communication backend
-  is OpenMPI 4.x (default), and OpenMPI 4.x or a later version must be pre-built.
-  The source files in the cuDSS tarball provide guidance for developing custom 
-  MPI implementations.
-  URL: https://developer.nvidia.com/cudss
-       https://docs.nvidia.com/cuda/cudss/advanced_features.html#communication-layer-library-in-cudss
-  Options: CUDSS_OPT, CUDSS_LIB.
-  Versions: cuDSS >= 0.6.0.
-
- Ginkgo (optional), used when MFEM_USE_GINKGO = YES. Ginkgo may have additional
-  requirements and module-specific dependencies; see the webpage below.
+- Ginkgo (optional), used when MFEM_USE_GINKGO = YES. Note that Ginkgo needs a
+  C++ compiler that supports the C++-17 standard. For additional requirements
+  and dependencies of specific modules, see the Ginkgo webpage below.
  URL: https://ginkgo-project.github.io
  Options: GINKGO_OPT, GINKGO_LIB, GINKGO_DIR, GINKGO_BUILD_TYPE (Release or
           Debug).
@@ -793,14 +771,6 @@ The specific libraries and their options are:
  Options: SLEPC_OPT, SLEPC_LIB.
  Versions: SLEPc >= 3.8.0.

- Sidre (optional), part of LLNL's axom project, used when MFEM_USE_SIDRE = YES.
-  Starting with MFEM v4.1, Axom version 0.3.1 or later is required.
-  URL: https://github.com/LLNL/axom
-       https://github.com/LLNL/conduit (Conduit)
-       https://support.hdfgroup.org/HDF5 (HDF5)
-  Options: SIDRE_OPT, SIDRE_LIB.
-  Versions: Axom >= 0.3.1.
-
 - Conduit (optional), used when MFEM_USE_CONDUIT = YES. Conduit Mesh Blueprint
  support requires Conduit >= v0.3.1 and VisIt >= v2.13.1 to read the output.
  URL: https://github.com/LLNL/conduit (Conduit)
@@ -808,7 +778,7 @@ The specific libraries and their options are:
  Options: CONDUIT_OPT, CONDUIT_LIB.
  Versions: Conduit >= 0.3.1.

- ADIOS2 (optional), used when MFEM_USE_ADIOS2 = YES.
+- ADIOS2 (optional) used when MFEM_USE_ADIOS2 = YES.
  URL: https://adios2.readthedocs.io/
  Versions: ADIOS >= 2.5.0.

@@ -884,7 +854,7 @@ The specific libraries and their options are:
  Options: RAJA_DIR, RAJA_OPT, RAJA_LIB.
  Versions: RAJA >= 2022.10.3.

- Moonolith (optional), used when MFEM_USE_MOONOLITH = YES.
+- Moonolith (optional), use when MFEM_USE_MOONOLITH = YES.
  URL: https://bitbucket.org/zulianp/par_moonolith
  Options: MOONOLITH_DIR
  Versions: MOONOLITH >= 1.1.0.
@@ -972,7 +942,7 @@ CMAKE_BUILD_TYPE which can be set to standard values like "Debug", and "Release"
 To use a specific generator use the "-G <generator>" option of cmake:

   cmake <mfem-source-dir> -G "Xcode"
-   cmake <mfem-source-dir> -G "Visual Studio 17 2022"
+   cmake <mfem-source-dir> -G "Visual Studio 12 2013"
   cmake <mfem-source-dir> -G "MinGW Makefiles"

 With CMake it is possible to build MFEM as a shared library using the standard
@@ -1084,7 +1054,6 @@ MFEM_USE_OCCA
 MFEM_USE_CEED
 MFEM_USE_RAJA
 MFEM_USE_UMPIRE
-MFEM_USE_SIDRE
 MFEM_USE_MOONOLITH
 MFEM_USE_CALIPER
 MFEM_USE_FMS
@@ -1148,7 +1117,6 @@ The CMake build system adds auto-detection for the following packages/libraries:
 - OCCA
 - RAJA
 - UMPIRE
- - AXOM - Used when MFEM_USE_SIDRE is enabled
 - MOONOLITH
 - CALIPER
 - FMS
@@ -1217,7 +1185,7 @@ larger problems, there are two options:
 Specific options for HIP
 ========================
 MFEM expects the `ROCM_PATH` environment variable to be set to the path of the
-ROCm install, as well as having `$ROCM_PATH/bin` in `PATH`.
+ROCM install, as well as having `$ROCM_PATH/bin` in `PATH`.

 Specific options for RAJA+HIP+MPI
 =================================
@@ -28,7 +28,6 @@ license files. These software products and their licenses are as follows:
 * AmgXWrapper (linalg/amgxsolver.{hpp,cpp}) -- MIT license
 * Catch++ (tests/unit/catch.hpp) -- Boost 1.0 license
 * Gecko (general/gecko.{cpp,hpp}) -- BSD 3-clause license
-* gslib (fem/gslib.{cpp,hpp}, mesh/bb_grid_map.{cpp,hpp}) -- BSD 3-clause license
 * Picojson (fem/picojson.h) -- Custom 2-clause license
 * TinyXML2 (general/tinyxml2.{cpp,h}) -- zlib license
 * Zstr (general/zstr.hpp) -- MIT license
@@ -248,10 +248,6 @@ IF (DEFINED TPL_ENABLE_MPFR)
  SET(MFEM_USE_MPFR ${TPL_ENABLE_MPFR} CACHE BOOL "Enable MPFR usage." FORCE)
 ENDIF()

-IF (DEFINED TPL_ENABLE_SIDRE)
-  SET(MFEM_USE_SIDRE ${TPL_ENABLE_SIDRE} CACHE BOOL "Enable Axom/Sidre usage" FORCE)
-ENDIF()
-
 IF (DEFINED TPL_ENABLE_FMS)
  SET(MFEM_USE_FMS ${TPL_ENABLE_FMS} CACHE BOOL "Enable FMS usage" FORCE)
 ENDIF()
@@ -35,7 +35,6 @@ set(MFEM_USE_SUITESPARSE @MFEM_USE_SUITESPARSE@)
 set(MFEM_USE_SUPERLU @MFEM_USE_SUPERLU@)
 set(MFEM_USE_MUMPS @MFEM_USE_MUMPS@)
 set(MFEM_USE_STRUMPACK @MFEM_USE_STRUMPACK@)
-set(MFEM_USE_CUDSS @MFEM_USE_CUDSS@)
 set(MFEM_USE_GINKGO @MFEM_USE_GINKGO@)
 set(MFEM_USE_AMGX @MFEM_USE_AMGX@)
 set(MFEM_USE_MAGMA @MFEM_USE_MAGMA@)
@@ -47,7 +46,6 @@ set(MFEM_USE_NETCDF @MFEM_USE_NETCDF@)
 set(MFEM_USE_PETSC @MFEM_USE_PETSC@)
 set(MFEM_USE_SLEPC @MFEM_USE_SLEPC@)
 set(MFEM_USE_MPFR @MFEM_USE_MPFR@)
-set(MFEM_USE_SIDRE @MFEM_USE_SIDRE@)
 set(MFEM_USE_FMS @MFEM_USE_FMS@)
 set(MFEM_USE_CONDUIT @MFEM_USE_CONDUIT@)
 set(MFEM_USE_PUMI @MFEM_USE_PUMI@)
@@ -110,14 +108,6 @@ if (MFEM_USE_RAJA)
   find_dependency(RAJA)
 endif()

-if (MFEM_USE_CUDSS)
-  find_dependency(cudss)
-endif (MFEM_USE_CUDSS)
-
-if (MFEM_USE_UMPIRE)
-   find_dependency(umpire)
-endif()
-
 if (NOT TARGET mfem)
  include(${CMAKE_CURRENT_LIST_DIR}/MFEMTargets.cmake)
 endif (NOT TARGET mfem)
@@ -108,15 +108,6 @@
 // Enable MFEM functionality based on the STRUMPACK library.
 #cmakedefine MFEM_USE_STRUMPACK

-// Enable MFEM functionality based on the cuDSS library.
-#cmakedefine MFEM_USE_CUDSS
-
-// CUDSS communication layer library path
-#cmakedefine MFEM_CUDSS_COMM_LIB "@MFEM_CUDSS_COMM_LIB@"
-
-// CUDSS threading layer library path
-#cmakedefine MFEM_CUDSS_THREADING_LIB "@MFEM_CUDSS_THREADING_LIB@"
-
 // Enable functionality based on the Ginkgo library.
 #cmakedefine MFEM_USE_GINKGO

@@ -129,9 +120,6 @@
 // Enable secure socket streams based on the GNUTLS library.
 #cmakedefine MFEM_USE_GNUTLS

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

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

 include(MfemCmakeUtilities)
 # Note: components are enabled based on the find_package() parameters.
@@ -0,0 +1,36 @@
+# Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+# Defines the following variables:
+#   - CAMP_FOUND
+#   - CAMP_LIBRARIES
+#   - CAMP_INCLUDE_DIRS
+
+include(MfemCmakeUtilities)
+
+mfem_find_package(CAMP CAMP CAMP_DIR
+      "include" "camp/camp.hpp"
+      "lib" "camp"
+      "Paths to headers required by CAMP."
+      "Libraries required by CAMP.")
+
+# RAJA commonly lists "camp" in INTERFACE_LINK_LIBRARIES. If there is no CMake
+# target named "camp", CMake treats it as a bare library name (-lcamp).
+if (CAMP_FOUND AND NOT TARGET camp)
+   list(GET CAMP_LIBRARIES 0 _camp_lib0)
+   add_library(camp UNKNOWN IMPORTED)
+   set_target_properties(camp PROPERTIES
+      IMPORTED_LOCATION "${_camp_lib0}"
+      INTERFACE_INCLUDE_DIRECTORIES "${CAMP_INCLUDE_DIRS}")
+   set(CAMP_LIBRARIES "camp" CACHE STRING "CAMP imported target." FORCE)
+   unset(_camp_lib0)
+endif()
+
@@ -1,68 +0,0 @@
-if (NOT cudss_DIR AND CUDSS_DIR)
-  set(cudss_DIR ${CUDSS_DIR}/lib/cmake/cudss)
-endif()
-message(STATUS "Looking for CUDSS ...")
-message(STATUS "   in CUDSS_DIR = ${CUDSS_DIR}")
-message(STATUS "      cudss_DIR = ${cudss_DIR}")
-find_package(cudss)
-set(CUDSS_FOUND ${cudss_FOUND})
-set(CUDSS_LIBRARIES "cudss")
-if (CUDSS_FOUND)
-  message(STATUS
-    "Found CUDSS target: ${CUDSS_LIBRARIES} (version: ${cudss_VERSION})")
-else()
-  set(msg STATUS)
-  if (CUDSS_FIND_REQUIRED)
-    set(msg FATAL_ERROR)
-  endif()
-  message(${msg}
-    "CUDSS not found. Please set CUDSS_DIR to the install prefix.")
-endif()
-
-if(CUDSS_FOUND AND TARGET cudss)
-  get_target_property(CUDSS_LIBRARY_LOCATION cudss IMPORTED_LOCATION)
-  if(NOT CUDSS_LIBRARY_LOCATION)
-    get_target_property(CUDSS_LIBRARY_LOCATION cudss IMPORTED_LOCATION_RELEASE)
-  endif()
-  if(CUDSS_LIBRARY_LOCATION)
-    get_filename_component(CUDSS_LIBRARY_DIR "${CUDSS_LIBRARY_LOCATION}" DIRECTORY)
-  else()
-    message(WARNING "Could not determine the location of the cuDSS library.")
-  endif()
-else()
-    message(WARNING "cuDSS target not available; cannot determine library directory.")
-endif()
-
-# Set the full name of the cuDSS threading library if OpenMP is enabled.
-# The threading layer library (libcudss_mtlayer_gomp.so) is located under the
-# cuDSS library directory by default.
-if (MFEM_USE_OPENMP)
-  find_file(
-    CUDSS_THREADING_LIB
-    NAMES libcudss_mtlayer_gomp.so
-    PATHS ${CUDSS_LIBRARY_DIR}
-    NO_DEFAULT_PATH
-  )
-  if (NOT DEFINED MFEM_CUDSS_THREADING_LIB AND CUDSS_THREADING_LIB)
-    set(MFEM_CUDSS_THREADING_LIB "${CUDSS_THREADING_LIB}")
-  endif()
-  message(STATUS "CUDSS threading layer library: ${MFEM_CUDSS_THREADING_LIB}")
-endif()
-
-# Set the full name of the cuDSS communication library if MFEM use OpenMPI.
-# The communication layer library (libcudss_commlayer_mpi.so) is located under the
-# cuDSS library directory by default.
-# The communication layer library is used pre-built communication layers for OpenMPI 
-# by default.
-if (MFEM_USE_MPI)
-  find_file(
-    CUDSS_COMM_LIB
-    NAMES libcudss_commlayer_openmpi.so
-    PATHS ${CUDSS_LIBRARY_DIR}
-    NO_DEFAULT_PATH
-  )
-  if (NOT DEFINED MFEM_CUDSS_COMM_LIB AND CUDSS_COMM_LIB)
-    set(MFEM_CUDSS_COMM_LIB "${CUDSS_COMM_LIB}")
-  endif()
-  message(STATUS "CUDSS communication layer library: ${MFEM_CUDSS_COMM_LIB}")
-endif()
@@ -21,6 +21,21 @@ mfem_find_package(Caliper CALIPER CALIPER_DIR
      "Paths to headers required by Caliper."
      "Libraries required by Caliper.")

+# Some downstream CMake packages (notably RAJA) may list "caliper" in their
+# INTERFACE_LINK_LIBRARIES. If there is no CMake target named "caliper", CMake
+# treats it as a bare library name and will pass -lcaliper to the linker.
+# Create a minimal imported target when we only located the library by path.
+if (CALIPER_FOUND AND NOT TARGET caliper)
+   list(GET CALIPER_LIBRARIES 0 _caliper_lib0)
+   add_library(caliper UNKNOWN IMPORTED)
+   set_target_properties(caliper PROPERTIES
+      IMPORTED_LOCATION "${_caliper_lib0}"
+      INTERFACE_INCLUDE_DIRECTORIES "${CALIPER_INCLUDE_DIRS}")
+   # Prefer linking via the target.
+   set(CALIPER_LIBRARIES "caliper" CACHE STRING "Caliper imported target." FORCE)
+   unset(_caliper_lib0)
+endif()
+
 # Append adiak path/lib if the user provided ADIAK_DIR
 if(ADIAK_DIR AND EXISTS ${ADIAK_DIR})
    find_package(adiak NO_DEFAULT_PATH REQUIRED PATHS ${ADIAK_DIR}/lib/cmake/adiak ${ADIAK_DIR})
@@ -18,17 +18,19 @@

 if (MFEM_FETCH_GSLIB OR MFEM_FETCH_TPLS)
  enable_language(C)
-  set(GSLIB_FETCH_VERSION 1.0.9)
-  add_library(GSLIB STATIC IMPORTED)
-  # set options (technically flags because GSLIB does not use cmake)
  string(TOUPPER "${CMAKE_BUILD_TYPE}" BUILD_TYPE)
-  set(GSLIB_FLAGS "${CMAKE_C_FLAGS} ${CMAKE_C_FLAGS_${BUILD_TYPE}}")
-  if (BUILD_SHARED_LIBS)
-    set(GSLIB_FLAGS "${GSLIB_FLAGS} -fPIC")
+  set(GSLIB_FETCH_VERSION 1.0.9)
+  set(GSLIB_C_FLAGS ${CMAKE_C_FLAGS_${BUILD_TYPE}})
+  if (CMAKE_C_FLAGS)
+    set(GSLIB_C_FLAGS "${CMAKE_C_FLAGS} ${CMAKE_C_FLAGS_${BUILD_TYPE}}")
  endif()
+  if (BUILD_SHARED_LIBS)
+    set(GSLIB_C_FLAGS "${GSLIB_C_FLAGS} -fPIC")
+  endif()
+  add_library(GSLIB STATIC IMPORTED)
  # define external project and create future include directory so it is present
  # to pass CMake checks at end of MFEM configuration step
-  message(STATUS "Will fetch GSLIB ${GSLIB_FETCH_VERSION} to be built with ${GSLIB_FLAGS}")
+  message(STATUS "Will fetch GSLIB ${GSLIB_FETCH_VERSION} to be built with ${GSLIB_C_FLAGS}")
  set(PREFIX ${CMAKE_BINARY_DIR}/fetch/gslib)
  include(ExternalProject)
  ExternalProject_Add(gslib
@@ -38,7 +40,7 @@ if (MFEM_FETCH_GSLIB OR MFEM_FETCH_TPLS)
    UPDATE_DISCONNECTED TRUE
    PREFIX ${PREFIX}
    CONFIGURE_COMMAND ""
-    BUILD_COMMAND cd ${PREFIX}/src/gslib && $(MAKE) clean && $(MAKE) DESTDIR=${PREFIX} MPI=$<BOOL:${MFEM_USE_MPI}> "CFLAGS=${GSLIB_FLAGS}"
+    BUILD_COMMAND cd ${PREFIX}/src/gslib && $(MAKE) clean && $(MAKE) DESTDIR=${PREFIX} MPI=$<BOOL:${MFEM_USE_MPI}> "CFLAGS= ${GSLIB_C_FLAGS}"
    INSTALL_COMMAND "")
  file(MAKE_DIRECTORY ${PREFIX}/include)
  # set imported library target properties
@@ -44,9 +44,6 @@ if (MFEM_FETCH_HYPRE OR MFEM_FETCH_TPLS)
  # set options and associated dependencies
  set(HYPRE_CMAKE_OPTIONS "")
  list(APPEND HYPRE_CMAKE_OPTIONS -DCMAKE_BUILD_TYPE:STRING=${CMAKE_BUILD_TYPE})
-  if (BUILD_SHARED_LIBS)
-    list(APPEND HYPRE_CMAKE_OPTIONS -DCMAKE_POSITION_INDEPENDENT_CODE:BOOL=ON)
-  endif()
  # collect all HYPRE_ENABLE variables and pass them to hypre, assuming they are BOOL.
  get_cmake_property(all_vars VARIABLES)
  foreach(var ${all_vars})
@@ -98,6 +95,7 @@ if (MFEM_FETCH_HYPRE OR MFEM_FETCH_TPLS)
    UPDATE_DISCONNECTED TRUE
    SOURCE_SUBDIR src
    PREFIX ${HYPRE_INSTALL}
+    BUILD_COMMAND ${CMAKE_COMMAND} --build . -- -j${CMAKE_BUILD_PARALLEL_LEVEL}
    CMAKE_CACHE_ARGS -DCMAKE_INSTALL_PREFIX:PATH=${HYPRE_INSTALL} -DCMAKE_INSTALL_LIBDIR:PATH=lib ${HYPRE_CMAKE_OPTIONS})
  file(MAKE_DIRECTORY ${HYPRE_INSTALL}/include)
  # set imported library target properties
@@ -19,18 +19,10 @@
 #   - METIS_VERSION_5 (cache variable)

 if (MFEM_FETCH_METIS OR MFEM_FETCH_TPLS)
-  enable_language(C)
  set(METIS_FETCH_VERSION 4.0.3)
  add_library(METIS STATIC IMPORTED)
-  # set options (technically flags because METIS does not use cmake)
-  set(METIS_FLAGS "-Wno-implicit-int -Wno-incompatible-pointer-types")
-  string(TOUPPER "${CMAKE_BUILD_TYPE}" BUILD_TYPE)
-  set(METIS_FLAGS "${METIS_FLAGS} ${CMAKE_C_FLAGS} ${CMAKE_C_FLAGS_${BUILD_TYPE}}")
-  if (BUILD_SHARED_LIBS)
-    set(METIS_FLAGS "${METIS_FLAGS} -fPIC")
-  endif()
  # define external project
-  message(STATUS "Will fetch METIS ${METIS_FETCH_VERSION} to be built with ${METIS_FLAGS}")
+  message(STATUS "Will fetch METIS ${METIS_FETCH_VERSION} to be built with default options")
  set(PREFIX ${CMAKE_BINARY_DIR}/fetch/metis)
  include(ExternalProject)
  ExternalProject_Add(metis
@@ -40,7 +32,7 @@ if (MFEM_FETCH_METIS OR MFEM_FETCH_TPLS)
    UPDATE_DISCONNECTED TRUE
    PREFIX ${PREFIX}
    CONFIGURE_COMMAND tar -xzf ../metis/metis-${METIS_FETCH_VERSION}-mac.tgz --strip=1
-    BUILD_COMMAND $(MAKE) clean && $(MAKE) "OPTFLAGS=${METIS_FLAGS}"
+    BUILD_COMMAND $(MAKE) COPTIONS=-Wno-incompatible-pointer-types
    INSTALL_COMMAND mkdir -p ${PREFIX}/lib && cp libmetis.a ${PREFIX}/lib/)
  # set imported library target properties
  add_dependencies(METIS metis)
@@ -22,15 +22,15 @@ include(MfemCmakeUtilities)
 mfem_find_package(SuiteSparse SuiteSparse SuiteSparse_DIR "" "" "" ""
  "Paths to headers required by SuiteSparse."
  "Libraries required by SuiteSparse."
-  ADD_COMPONENT "UMFPACK" "include;include/suitesparse;suitesparse" umfpack.h "lib" umfpack
-  ADD_COMPONENT "KLU" "include;include/suitesparse;suitesparse" klu.h "lib" klu
-  ADD_COMPONENT "AMD" "include;include/suitesparse;suitesparse" amd.h "lib" amd
-  ADD_COMPONENT "BTF" "include;include/suitesparse;suitesparse" btf.h "lib" btf
-  ADD_COMPONENT "CHOLMOD" "include;include/suitesparse;suitesparse" cholmod.h "lib" cholmod
-  ADD_COMPONENT "COLAMD" "include;include/suitesparse;suitesparse" colamd.h "lib" colamd
-  ADD_COMPONENT "CAMD" "include;include/suitesparse;suitesparse" camd.h "lib" camd
-  ADD_COMPONENT "CCOLAMD" "include;include/suitesparse;suitesparse" ccolamd.h "lib" ccolamd
-  ADD_COMPONENT "config" "include;include/suitesparse;suitesparse" SuiteSparse_config.h "lib"
+  ADD_COMPONENT "UMFPACK" "include;suitesparse" umfpack.h "lib" umfpack
+  ADD_COMPONENT "KLU" "include;suitesparse" klu.h "lib" klu
+  ADD_COMPONENT "AMD" "include;suitesparse" amd.h "lib" amd
+  ADD_COMPONENT "BTF" "include;suitesparse" btf.h "lib" btf
+  ADD_COMPONENT "CHOLMOD" "include;suitesparse" cholmod.h "lib" cholmod
+  ADD_COMPONENT "COLAMD" "include;suitesparse" colamd.h "lib" colamd
+  ADD_COMPONENT "CAMD" "include;suitesparse" camd.h "lib" camd
+  ADD_COMPONENT "CCOLAMD" "include;suitesparse" ccolamd.h "lib" ccolamd
+  ADD_COMPONENT "config" "include;suitesparse" SuiteSparse_config.h "lib"
    suitesparseconfig)

 if (SuiteSparse_FOUND AND METIS_VERSION_5)
@@ -14,12 +14,12 @@
 #   - UMPIRE_LIBRARIES
 #   - UMPIRE_INCLUDE_DIRS

-if (NOT umpire_ROOT AND UMPIRE_DIR)
-  set(umpire_ROOT ${UMPIRE_DIR})
+if (NOT umpire_DIR AND UMPIRE_DIR)
+  set(umpire_DIR ${UMPIRE_DIR}/lib/cmake/umpire)
 endif()
 message(STATUS "Looking for UMPIRE ...")
 message(STATUS "   in UMPIRE_DIR = ${UMPIRE_DIR}")
-message(STATUS "      umpire_ROOT = ${umpire_ROOT}")
+message(STATUS "      umpire_DIR = ${umpire_DIR}")
 find_package(umpire CONFIG)
 set(UMPIRE_FOUND ${umpire_FOUND})
 set(UMPIRE_LIBRARIES "umpire")
@@ -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_SIDRE MFEM_USE_FMS MFEM_USE_CONDUIT MFEM_USE_PUMI
+      MFEM_USE_MPFR MFEM_USE_FMS MFEM_USE_CONDUIT MFEM_USE_PUMI
      MFEM_USE_HIOP MFEM_USE_GSLIB MFEM_USE_CUDA MFEM_USE_HIP MFEM_USE_RAJA
      MFEM_USE_OCCA MFEM_USE_CEED MFEM_USE_CALIPER MFEM_USE_UMPIRE MFEM_USE_SIMD
      MFEM_USE_ADIOS2 MFEM_USE_MKL_CPARDISO MFEM_USE_MKL_PARDISO
@@ -157,10 +157,4 @@ constexpr real_t operator""_r(unsigned long long v)
 #endif
 #endif // MFEM_USE_MPI not defined

-#ifndef MFEM_USE_CUDA
-#ifdef MFEM_USE_CUDSS
-#error Building with cuDSS (MFEM_USE_CUDSS=YES) requires CUDA (MFEM_USE_CUDA=YES)
-#endif
-#endif // MFEM_USE_CUDSS not defined
-
 #endif // MFEM_CONFIG_HPP
@@ -108,15 +108,6 @@
 // Enable MFEM functionality based on the STRUMPACK library.
 // #define MFEM_USE_STRUMPACK

-// Enable MFEM functionality based on the cuDSS library.
-// #define MFEM_USE_CUDSS
-
-// CUDSS communication layer library path
-// #define MFEM_CUDSS_COMM_LIB "@MFEM_CUDSS_COMM_LIB@"
-
-// CUDSS threading layer library path
-// #define MFEM_CUDSS_THREADING_LIB "@MFEM_CUDSS_THREADING_LIB@"
-
 // Enable MFEM features based on the Ginkgo library.
 // #define MFEM_USE_GINKGO

@@ -129,9 +120,6 @@
 // Enable secure socket streams based on the GNUTLS library.
 // #define MFEM_USE_GNUTLS

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

@@ -36,9 +36,6 @@ MFEM_USE_SUPERLU       = @MFEM_USE_SUPERLU@
 MFEM_USE_SUPERLU5      = @MFEM_USE_SUPERLU5@
 MFEM_USE_MUMPS         = @MFEM_USE_MUMPS@
 MFEM_USE_STRUMPACK     = @MFEM_USE_STRUMPACK@
-MFEM_USE_CUDSS         = @MFEM_USE_CUDSS@
-MFEM_CUDSS_COMM_LIB    = @MFEM_CUDSS_COMM_LIB@
-MFEM_CUDSS_THREADING_LIB = @MFEM_CUDSS_THREADING_LIB@
 MFEM_USE_GINKGO        = @MFEM_USE_GINKGO@
 MFEM_USE_AMGX          = @MFEM_USE_AMGX@
 MFEM_USE_MAGMA         = @MFEM_USE_MAGMA@
@@ -48,7 +45,6 @@ MFEM_USE_NETCDF        = @MFEM_USE_NETCDF@
 MFEM_USE_PETSC         = @MFEM_USE_PETSC@
 MFEM_USE_SLEPC         = @MFEM_USE_SLEPC@
 MFEM_USE_MPFR          = @MFEM_USE_MPFR@
-MFEM_USE_SIDRE         = @MFEM_USE_SIDRE@
 MFEM_USE_FMS           = @MFEM_USE_FMS@
 MFEM_USE_CONDUIT       = @MFEM_USE_CONDUIT@
 MFEM_USE_PUMI          = @MFEM_USE_PUMI@
@@ -38,7 +38,6 @@ option(MFEM_USE_SUPERLU "Enable SuperLU_DIST usage" OFF)
 option(MFEM_USE_SUPERLU5 "Use the old SuperLU_DIST 5.1 version" OFF)
 option(MFEM_USE_MUMPS "Enable MUMPS usage" OFF)
 option(MFEM_USE_STRUMPACK "Enable STRUMPACK usage" OFF)
-option(MFEM_USE_CUDSS "Enable cuDSS usage" OFF)
 option(MFEM_USE_GINKGO "Enable Ginkgo usage" OFF)
 option(MFEM_USE_AMGX "Enable AmgX usage" OFF)
 option(MFEM_USE_MAGMA "Enable MAGMA usage" OFF)
@@ -49,7 +48,6 @@ option(MFEM_USE_NETCDF "Enable NETCDF usage" OFF)
 option(MFEM_USE_PETSC "Enable PETSc support." OFF)
 option(MFEM_USE_SLEPC "Enable SLEPc support." OFF)
 option(MFEM_USE_MPFR "Enable MPFR usage." OFF)
-option(MFEM_USE_SIDRE "Enable Axom/Sidre usage" OFF)
 option(MFEM_USE_FMS "Enable FMS usage" OFF)
 option(MFEM_USE_CONDUIT "Enable Conduit usage" OFF)
 option(MFEM_USE_PUMI "Enable PUMI" OFF)
@@ -225,17 +223,8 @@ set(FMS_DIR "${MFEM_DIR}/../fms" CACHE PATH
 set(CONDUIT_DIR "${MFEM_DIR}/../conduit" CACHE PATH
    "Path to the Conduit library.")

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

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

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

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

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

@@ -27,13 +27,7 @@ MPICXX = mpicxx

 BASE_FLAGS  = -std=c++17
 OPTIM_FLAGS = -O3 $(BASE_FLAGS)
-
-# Shadow warnings for clang only; GCC's -Wshadow flags more.
-SHADOW_WARNING_FLAG = $(if $(findstring clang,\
- $(shell $(MFEM_HOST_CXX) --version 2>/dev/null)),-Wshadow,)
-WARNING_FLAGS = -pedantic -Wall $(SHADOW_WARNING_FLAG)
-
-DEBUG_FLAGS = $(strip -g $(addprefix $(XCOMPILER),$(WARNING_FLAGS)) $(BASE_FLAGS))
+DEBUG_FLAGS = -g $(XCOMPILER)-Wall $(BASE_FLAGS)

 # Prefixes for passing flags to the compiler and linker when using CXX or MPICXX
 CXX_XCOMPILER =
@@ -159,7 +153,6 @@ MFEM_USE_SUPERLU       = NO
 MFEM_USE_SUPERLU5      = NO
 MFEM_USE_MUMPS         = NO
 MFEM_USE_STRUMPACK     = NO
-MFEM_USE_CUDSS         = NO
 MFEM_USE_GINKGO        = NO
 MFEM_USE_AMGX          = NO
 MFEM_USE_MAGMA         = NO
@@ -169,7 +162,6 @@ MFEM_USE_NETCDF        = NO
 MFEM_USE_PETSC         = NO
 MFEM_USE_SLEPC         = NO
 MFEM_USE_MPFR          = NO
-MFEM_USE_SIDRE         = NO
 MFEM_USE_FMS           = NO
 MFEM_USE_CONDUIT       = NO
 MFEM_USE_PUMI          = NO
@@ -256,6 +248,15 @@ endif

 # METIS library configuration
 ifeq ($(MFEM_USE_SUPERLU)$(MFEM_USE_STRUMPACK)$(MFEM_USE_MUMPS),NONONO)
+   # MFEM_USE_METIS_5: when the user supplies METIS_DIR, try to auto-detect
+   # METIS 5 installs that follow the common <prefix>/{include,lib,lib64} layout.
+   ifeq ($(MFEM_USE_METIS_5),NO)
+      ifneq ($(wildcard $(METIS_DIR)/include/metis.h),)
+         ifneq ($(wildcard $(METIS_DIR)/lib/libmetis.* $(METIS_DIR)/lib64/libmetis.*),)
+            MFEM_USE_METIS_5 = YES
+         endif
+      endif
+   endif
   ifeq ($(MFEM_USE_METIS_5),NO)
     METIS_DIR = @MFEM_DIR@/../metis-4.0
     METIS_OPT =
@@ -375,19 +376,6 @@ STRUMPACK_OPT = -I$(STRUMPACK_DIR)/include $(SCOTCH_OPT)
 STRUMPACK_LIB = -L$(STRUMPACK_DIR)/lib -lstrumpack $(MPI_FORTRAN_LIB)\
 $(SCOTCH_LIB) $(SCALAPACK_LIB)

-# CUDSS library configuration
-CUDSS_DIR         = @MFEM_DIR@/../cudss
-CUDSS_INCLUDE_DIR = $(CUDSS_DIR)/include
-CUDSS_LIBRARY_DIR = $(CUDSS_DIR)/lib
-CUDSS_OPT         = -I$(CUDSS_INCLUDE_DIR)
-CUDSS_LIB         = \
- $(XLINKER)-rpath,$(CUDSS_LIBRARY_DIR) -L$(CUDSS_LIBRARY_DIR) -lcudss
-# The cuDSS communication and threading libraries. 
-MFEM_CUDSS_COMM_LIB = $(abspath $(wildcard $(or $(CUDSS_COMM_LIB),\
-   $(subst @MFEM_DIR@,$(MFEM_DIR), $(CUDSS_LIBRARY_DIR)/libcudss_commlayer_openmpi.so))))
-MFEM_CUDSS_THREADING_LIB = $(abspath $(wildcard $(or $(CUDSS_THREADING_LIB),\
-   $(subst @MFEM_DIR@,$(MFEM_DIR),$(CUDSS_LIBRARY_DIR)/libcudss_mtlayer_gomp.so))))
-
 # Ginkgo library configuration
 GINKGO_DIR = @MFEM_DIR@/../ginkgo/install
 GINKGO_SEARCH_DIR = $(subst @MFEM_DIR@,$(MFEM_DIR),$(GINKGO_DIR))
@@ -507,17 +495,6 @@ ifneq (,$(wildcard $(CONDUIT_HDF5_HEADER)))
                  -lhdf5 $(ZLIB_LIB)
 endif

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

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

 # Enzyme configuration
 ENZYME_DIR = @MFEM_DIR@/../enzyme
@@ -39,8 +39,3 @@ when a picture was added for documentation.
 If that is the case, make sure the failure is indeed justified, and rerun the
 push command with the `--no-verify` option. This will skip the hooks, allowing
 you to push those changes.
-
-The `branch-history` check is run automatically through GitHub Actions. If a
-branch is known to have a large number of changes that are legitimate, the
-check can be overridden by setting  the label 'branch-history-override' on the
-pull request.
@@ -1,38 +0,0 @@
-MFEM mesh v1.0
-
-#
-# MFEM Geometry Types (see fem/geom.hpp):
-#
-# POINT       = 0
-# SEGMENT     = 1
-# TRIANGLE    = 2
-# SQUARE      = 3
-# TETRAHEDRON = 4
-# CUBE        = 5
-# PRISM       = 6
-# PYRAMID     = 7
-
-dimension
-2
-
-elements
-2
-1 3 0 1 4 3
-1 2 1 2 4
-
-boundary
-5
-1 1 0 1
-1 1 1 2
-1 1 2 4
-1 1 4 3
-1 1 3 0
-
-vertices
-5
-2
-0 0
-1 0
-2 0
-0 1
-1 1
@@ -201,7 +201,6 @@ namespace mfem {
 * - <a class="el" href="nurbs__naca__cmesh_8cpp_source.html">NURBS NACA Mesher</a>: generate NURBS based mesh around a NACA foil
 * - <a class="el" href="nurbs__printfunc_8cpp_source.html">NURBS Printer</a>: print the NURBS-basis
 * - <a class="el" href="nurbs__mesh_info_8cpp_source.html">NURBS Mesh info</a>: print the info of a NURBS mesh
- * - <a class="el" href="nurbs__surface_8cpp_source.html">NURBS Surface</a>: interpolate a 3D Surface in a NURBS Patch
 *
 * <H3>Miniapps</H3>
 * - <a class="el" href="volta_8cpp_source.html">Volta</a>: simple electrostatics simulation code
@@ -246,9 +245,6 @@ namespace mfem {
 * - <a class="el" href="pdiffusion_8cpp_source.html">DPG Diffusion example</a>: DPG formulation for the diffusion problem
 * - <a class="el" href="pmaxwell_8cpp_source.html">DPG Maxwell example</a>: DPG formulation for the indefinite Maxwell problem
 * - <a class="el" href="lor__elast_8cpp_source.html">LOR Elasticity</a>: solve linear elasticity with LOR preconditioning on GPUs
- * - <a class="el" href="reflector_8cpp_source.html">Reflector Miniapp</a>: reflect a mesh about a plane
- * - <a class="el" href="ref321_8cpp_source.html">3:1 Refinement Miniapp</a>: perform 3:1 anisotropic mesh refinements
- * - <a class="el" href="pref321_8cpp_source.html">3:1 Refinement Miniapp</a>: parallel 3:1 anisotropic mesh refinements
 *
 * See also the <a class="el" href="https://mfem.org/examples/">examples documentation</a> online.
 */
@@ -215,7 +215,7 @@ if (MFEM_ENABLE_TESTING)
    add_test(NAME ex1p_ceed_np=${MFEM_MPI_NP}
      COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
      ${MPIEXEC_PREFLAGS}
-      $<TARGET_FILE:ex1p> "-no-vis" "-d" "ceed-cpu" "-pa" "-a"
+      $<TARGET_FILE:ex1p> "-no-vis" "-d ceed-cpu" "-pa" "-a"
      ${MPIEXEC_POSTFLAGS})
  endif()
 endif()
@@ -64,7 +64,7 @@ PARALLEL_NAME := Parallel AMGX example
 $(MFEM_LIB_FILE):
 	$(error The MFEM library is not build)

-clean: clean-build clean-exec
+clean: clean-build

 clean-build:
 	rm -f *.o *~ $(SEQ_EXAMPLES) $(PAR_EXAMPLES)
@@ -64,12 +64,12 @@ ex1p-test-par: ex1p
 $(MFEM_LIB_FILE):
 	$(error The MFEM library is not built)

-clean: clean-build clean-exec
+clean: clean-build clean-exec $(SUBDIRS_CLEAN)

 clean-build:
 	rm -f *.o *~ $(SEQ_EXAMPLES) $(PAR_EXAMPLES)
 	rm -rf *.dSYM *.TVD.*breakpoints

 clean-exec:
-	@rm -f refined.mesh mesh.*
-	@rm -f sol.*
+	@rm -f refined.mesh displaced.mesh mesh.* ex5.mesh
+	@rm -f sphere_refined.* sol.* sol_u.* sol_p.* sol_r.* sol_i.*
@@ -50,10 +50,6 @@
 //               ex1 -m ../data/beam-tet.mesh -pa -d ceed-cpu
 //               ex1 -m ../data/beam-tet.mesh -pa -d ceed-cuda:/gpu/cuda/ref
 //
-// Device simplices sample runs:
-//               ex1 -pa -d gpu -m ../data/inline-tet.mesh
-//               ex1 -pa -d gpu -m ../data/inline-tri.mesh
-//
 // Description:  This example code demonstrates the use of MFEM to define a
 //               simple finite element discretization of the Poisson problem
 //               -Delta u = 1 with homogeneous Dirichlet boundary conditions.
@@ -142,25 +138,25 @@ int main(int argc, char *argv[])
   }

   // 5. Define a finite element space on the mesh. Here we use continuous
-   //    Lagrange finite elements of the specified order.
-   //    - If order < 1, we instead use an isoparametric/isogeometric space.
-   //    - If the mesh is simplicial and partial assembly is requested,
-   //      we use the positive basis, which supports device execution.
+   //    Lagrange finite elements of the specified order. If order < 1, we
+   //    instead use an isoparametric/isogeometric space.
   FiniteElementCollection *fec;
-   auto basis_type = (pa && mesh.IsSimplexMesh()) ?
-                     BasisType::Positive : BasisType::GaussLobatto;
+   bool delete_fec;
   if (order > 0)
   {
-      fec = new H1_FECollection(order, dim, basis_type);
+      fec = new H1_FECollection(order, dim);
+      delete_fec = true;
   }
   else if (mesh.GetNodes())
   {
      fec = mesh.GetNodes()->OwnFEC();
+      delete_fec = false;
      cout << "Using isoparametric FEs: " << fec->Name() << endl;
   }
   else
   {
-      fec = new H1_FECollection(order = 1, dim, basis_type);
+      fec = new H1_FECollection(order = 1, dim);
+      delete_fec = true;
   }
   FiniteElementSpace fespace(&mesh, fec);
   cout << "Number of finite element unknowns: "
@@ -228,29 +224,17 @@ int main(int argc, char *argv[])
   // 11. Solve the linear system A X = B.
   if (!pa)
   {
-#ifdef MFEM_USE_CUDSS
-      if (Device::Allows(Backend::CUDA_MASK))
-      {
-         // Use cuDSS to solve the system.
-         CuDSSSolver cudss_solver;
-         cudss_solver.SetOperator(*A);
-         cudss_solver.Mult(B, X);
-      }
-      else
-#endif
-      {
 #ifndef MFEM_USE_SUITESPARSE
-         // Use a simple symmetric Gauss-Seidel preconditioner with PCG.
-         GSSmoother M((SparseMatrix&)(*A));
-         PCG(*A, M, B, X, 1, 200, 1e-12, 0.0);
+      // Use a simple symmetric Gauss-Seidel preconditioner with PCG.
+      GSSmoother M((SparseMatrix&)(*A));
+      PCG(*A, M, B, X, 1, 200, 1e-12, 0.0);
 #else
-         // If MFEM was compiled with SuiteSparse, use UMFPACK to solve the system.
-         UMFPackSolver umf_solver;
-         umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
-         umf_solver.SetOperator(*A);
-         umf_solver.Mult(B, X);
+      // If MFEM was compiled with SuiteSparse, use UMFPACK to solve the system.
+      UMFPackSolver umf_solver;
+      umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
+      umf_solver.SetOperator(*A);
+      umf_solver.Mult(B, X);
 #endif
-      }
   }
   else
   {
@@ -289,14 +273,17 @@ int main(int argc, char *argv[])
   if (visualization)
   {
      char vishost[] = "localhost";
-      int visport = 19916;
+      int  visport   = 19916;
      socketstream sol_sock(vishost, visport);
      sol_sock.precision(8);
      sol_sock << "solution\n" << mesh << x << flush;
   }

   // 15. Free the used memory.
-   if (order > 0) { delete fec; }
+   if (delete_fec)
+   {
+      delete fec;
+   }

   return 0;
 }
@@ -42,11 +42,7 @@
 //               mpirun -np 4 ex1p -pa -d ceed-cuda:/gpu/cuda/shared
 //               mpirun -np 4 ex1p -pa -d ceed-cuda:/gpu/cuda/shared -m ../data/square-mixed.mesh
 //               mpirun -np 4 ex1p -pa -d ceed-cuda:/gpu/cuda/shared -m ../data/fichera-mixed.mesh
-//               mpirun -np 4 ex1p -pa -d ceed-cpu -m ../data/beam-tet.mesh
-//
-// Device simplices sample runs:
-//               mpirun -np 4 ex1p -pa -d gpu -m ../data/inline-tet.mesh
-//               mpirun -np 4 ex1p -pa -d gpu -m ../data/inline-tri.mesh
+//               mpirun -np 4 ex1p -m ../data/beam-tet.mesh -pa -d ceed-cpu
 //
 // Description:  This example code demonstrates the use of MFEM to define a
 //               simple finite element discretization of the Poisson problem
@@ -87,9 +83,6 @@ int main(int argc, char *argv[])
   const char *device_config = "cpu";
   bool visualization = true;
   bool algebraic_ceed = false;
-#ifdef MFEM_USE_CUDSS
-   bool cudss_solver = false;
-#endif

   OptionsParser args(argc, argv);
   args.AddOption(&mesh_file, "-m", "--mesh",
@@ -109,10 +102,6 @@ int main(int argc, char *argv[])
   args.AddOption(&algebraic_ceed, "-a", "--algebraic",
                  "-no-a", "--no-algebraic",
                  "Use algebraic Ceed solver");
-#endif
-#ifdef MFEM_USE_CUDSS
-   args.AddOption(&cudss_solver, "-cudss", "--cudss-solver", "-no-cudss",
-                  "--no-cudss-solver", "Use the cuDSS Solver.");
 #endif
   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                  "--no-visualization",
@@ -169,20 +158,19 @@ int main(int argc, char *argv[])
   }

   // 7. Define a parallel finite element space on the parallel mesh. Here we
-   //    use continuous Lagrange finite elements of the specified order.
-   //    - If order < 1, we instead use an isoparametric/isogeometric space.
-   //    - If the mesh is simplicial and partial assembly is requested,
-   //      we use the positive basis, which supports device execution.
+   //    use continuous Lagrange finite elements of the specified order. If
+   //    order < 1, we instead use an isoparametric/isogeometric space.
   FiniteElementCollection *fec;
-   auto basis_type = (pa && pmesh.IsSimplexMesh()) ?
-                     BasisType::Positive : BasisType::GaussLobatto;
+   bool delete_fec;
   if (order > 0)
   {
-      fec = new H1_FECollection(order, dim, basis_type);
+      fec = new H1_FECollection(order, dim);
+      delete_fec = true;
   }
   else if (pmesh.GetNodes())
   {
      fec = pmesh.GetNodes()->OwnFEC();
+      delete_fec = false;
      if (myid == 0)
      {
         cout << "Using isoparametric FEs: " << fec->Name() << endl;
@@ -190,7 +178,8 @@ int main(int argc, char *argv[])
   }
   else
   {
-      fec = new H1_FECollection(order = 1, dim, basis_type);
+      fec = new H1_FECollection(order = 1, dim);
+      delete_fec = true;
   }
   ParFiniteElementSpace fespace(&pmesh, fec);
   HYPRE_BigInt size = fespace.GlobalTrueVSize();
@@ -259,51 +248,33 @@ int main(int argc, char *argv[])
   // 13. Solve the linear system A X = B.
   //     * With full assembly, use the BoomerAMG preconditioner from hypre.
   //     * With partial assembly, use Jacobi smoothing, for now.
-#ifdef MFEM_USE_CUDSS
-   if (!pa && (Device::Allows(Backend::CUDA_MASK) && cudss_solver))
+   Solver *prec = NULL;
+   if (pa)
   {
-      // Solve using a direct solver with cuDSS
-      CuDSSSolver cudss_solver(MPI_COMM_WORLD);
-      cudss_solver.SetMatrixSymType(
-         CuDSSSolver::SYMMETRIC_POSITIVE_DEFINITE);
-      cudss_solver.SetMatrixViewType(CuDSSSolver::UPPER);
-      cudss_solver.SetOperator(*A);
-      cudss_solver.Mult(B, X);
-   }
-   else
-#endif
-   {
-      Solver *prec = NULL;
-      if (pa)
+      if (UsesTensorBasis(fespace))
      {
-         if (UsesTensorBasis(fespace))
+         if (algebraic_ceed)
         {
-            if (algebraic_ceed)
-            {
-               prec = new ceed::AlgebraicSolver(a, ess_tdof_list);
-            }
-            else
-            {
-               prec = new OperatorJacobiSmoother(a, ess_tdof_list);
-            }
+            prec = new ceed::AlgebraicSolver(a, ess_tdof_list);
+         }
+         else
+         {
+            prec = new OperatorJacobiSmoother(a, ess_tdof_list);
         }
      }
-      else
-      {
-         prec = new HypreBoomerAMG;
-      }
-      CGSolver cg(MPI_COMM_WORLD);
-      cg.SetRelTol(1e-12);
-      cg.SetMaxIter(2000);
-      cg.SetPrintLevel(1);
-      if (prec)
-      {
-         cg.SetPreconditioner(*prec);
-      }
-      cg.SetOperator(*A);
-      cg.Mult(B, X);
-      delete prec;
   }
+   else
+   {
+      prec = new HypreBoomerAMG;
+   }
+   CGSolver cg(MPI_COMM_WORLD);
+   cg.SetRelTol(1e-12);
+   cg.SetMaxIter(2000);
+   cg.SetPrintLevel(1);
+   if (prec) { cg.SetPreconditioner(*prec); }
+   cg.SetOperator(*A);
+   cg.Mult(B, X);
+   delete prec;

   // 14. Recover the parallel grid function corresponding to X. This is the
   //     local finite element solution on each processor.
@@ -337,7 +308,10 @@ int main(int argc, char *argv[])
   }

   // 17. Free the used memory.
-   if (order > 0) { delete fec; }
+   if (delete_fec)
+   {
+      delete fec;
+   }

   return 0;
 }
@@ -177,7 +177,7 @@ int main(int argc, char *argv[])
   Array<int> ess_tdof_list(0);
   if (h1 && mesh->bdr_attributes.Size())
   {
-      // For a continuous basis the linear system must be modifed to enforce an
+      // For a continuous basis the linear system must be modified to enforce an
      // essential (Dirichlet) boundary condition. In the DG case this is not
      // necessary as the boundary condition will only be enforced weakly.
      fespace.GetEssentialTrueDofs(dbc_bdr, ess_tdof_list);
@@ -95,15 +95,6 @@ int main(int argc, char *argv[])
      args.PrintOptions(cout);
   }

-   if (amg_elast && !static_cond && reorder_space)
-   {
-      if (myid == 0)
-         cerr << "\nThe AMG elasticity solver requires ordering byVDIM! "
-              << "Ignoring the specified option -nodes/--by-nodes.\n"
-              << endl;
-      reorder_space = false;
-   }
-
   // 3. Enable hardware devices such as GPUs, and programming models such as
   //    CUDA, OCCA, RAJA and OpenMP based on command line options.
   Device device(device_config);
@@ -5,8 +5,8 @@
 // Sample runs:
 //     ex37 -alpha 10
 //     ex37 -alpha 10 -pv
-//     ex37 -lambda 0.1 -mu 0.1 -growth 1
-//     ex37 -o 2 -alpha 10.0 -mi 50 -vf 0.4 -ntol 1e-5 -growth 1.5
+//     ex37 -lambda 0.1 -mu 0.1
+//     ex37 -o 2 -alpha 5.0 -mi 50 -vf 0.4 -ntol 1e-5
 //     ex37 -r 6 -o 1 -alpha 25.0 -epsilon 0.02 -mi 50 -ntol 1e-5
 //
 // Description: This example code demonstrates the use of MFEM to solve a
@@ -55,6 +55,53 @@
 using namespace std;
 using namespace mfem;

+/**
+ * @brief Bregman projection of ρ = sigmoid(ψ) onto the subspace
+ *        ∫_Ω ρ dx = θ vol(Ω) as follows:
+ *
+ *        1. Compute the root of the R → R function
+ *            f(c) = ∫_Ω sigmoid(ψ + c) dx - θ vol(Ω)
+ *        2. Set ψ ← ψ + c.
+ *
+ * @param psi a GridFunction to be updated
+ * @param target_volume θ vol(Ω)
+ * @param tol Newton iteration tolerance
+ * @param max_its Newton maximum iteration number
+ * @return real_t Final volume, ∫_Ω sigmoid(ψ)
+ */
+real_t proj(GridFunction &psi, real_t target_volume, real_t tol=1e-12,
+            int max_its=10)
+{
+   MappedGridFunctionCoefficient sigmoid_psi(&psi, sigmoid);
+   MappedGridFunctionCoefficient der_sigmoid_psi(&psi, der_sigmoid);
+
+   LinearForm int_sigmoid_psi(psi.FESpace());
+   int_sigmoid_psi.AddDomainIntegrator(new DomainLFIntegrator(sigmoid_psi));
+   LinearForm int_der_sigmoid_psi(psi.FESpace());
+   int_der_sigmoid_psi.AddDomainIntegrator(new DomainLFIntegrator(
+                                              der_sigmoid_psi));
+   bool done = false;
+   for (int k=0; k<max_its; k++) // Newton iteration
+   {
+      int_sigmoid_psi.Assemble(); // Recompute f(c) with updated ψ
+      const real_t f = int_sigmoid_psi.Sum() - target_volume;
+
+      int_der_sigmoid_psi.Assemble(); // Recompute df(c) with updated ψ
+      const real_t df = int_der_sigmoid_psi.Sum();
+
+      const real_t dc = -f/df;
+      psi += dc;
+      if (abs(dc) < tol) { done = true; break; }
+   }
+   if (!done)
+   {
+      mfem_warning("Projection reached maximum iteration without converging. "
+                   "Result may not be accurate.");
+   }
+   int_sigmoid_psi.Assemble();
+   return int_sigmoid_psi.Sum();
+}
+
 /*
 * ---------------------------------------------------------------
 *                      ALGORITHM PREAMBLE
@@ -133,11 +180,10 @@ int main(int argc, char *argv[])
   int ref_levels = 5;
   int order = 2;
   real_t alpha = 1.0;
-   real_t growth = 2;
   real_t epsilon = 0.01;
   real_t vol_fraction = 0.5;
   int max_it = 1e3;
-   real_t itol = 1e-2;
+   real_t itol = 1e-1;
   real_t ntol = 1e-4;
   real_t rho_min = 1e-6;
   real_t lambda = 1.0;
@@ -152,8 +198,6 @@ int main(int argc, char *argv[])
                  "Order (degree) of the finite elements.");
   args.AddOption(&alpha, "-alpha", "--alpha-step-length",
                  "Step length for gradient descent.");
-   args.AddOption(&growth, "-growth", "--alpha-growth-rate",
-                  "Growth rate of step length for gradient descent.");
   args.AddOption(&epsilon, "-epsilon", "--epsilon-thickness",
                  "Length scale for ρ.");
   args.AddOption(&max_it, "-mi", "--max-it",
@@ -288,7 +332,6 @@ int main(int argc, char *argv[])
   }
   FilterSolver->SetEssentialBoundary(ess_bdr_filter);
   FilterSolver->SetupFEM();
-   FilterSolver->AssembleDiffusionBilinear();

   BilinearForm mass(&control_fes);
   mass.AddDomainIntegrator(new InverseIntegrator(new MassIntegrator(one)));
@@ -342,7 +385,7 @@ int main(int argc, char *argv[])
   // 11. Iterate:
   for (int k = 1; k <= max_it; k++)
   {
-      if (k > 1) { alpha = std::pow((real_t) k,growth); }
+      if (k > 1) { alpha *= ((real_t) k) / ((real_t) k-1); }

      mfem::out << "\nStep = " << k << std::endl;

@@ -379,9 +422,7 @@ int main(int argc, char *argv[])

      // Step 5 - Update design variable ψ ← proj(ψ - αG)
      psi.Add(-alpha, grad);
-      GridFunction alpha_grad(grad);
-      alpha_grad *= alpha;
-      const real_t material_volume = proj(psi, alpha_grad, target_volume);
+      const real_t material_volume = proj(psi, target_volume);

      // Compute ||ρ - ρ_old|| in control fes.
      real_t norm_increment = zerogf.ComputeL1Error(succ_diff_rho);
@@ -137,7 +137,7 @@ public:
        exponent(exponent_), rho_min(rho_min_)
   {
      MFEM_ASSERT(rho_min_ >= 0.0, "rho_min must be >= 0");
-      MFEM_ASSERT(rho_min_ < 1.0,  "rho_min must be < 1");
+      MFEM_ASSERT(rho_min_ < 1.0,  "rho_min must be > 1");
      MFEM_ASSERT(u, "displacement field is not set");
      MFEM_ASSERT(rho_filter, "density field is not set");
   }
@@ -231,12 +231,9 @@ private:
   FiniteElementCollection * fec = nullptr;
   FiniteElementSpace * fes = nullptr;
   Array<int> ess_bdr;
-   Array<int> ess_tdof_list;
   Array<int> neumann_bdr;
   GridFunction * u = nullptr;
   LinearForm * b = nullptr;
-   BilinearForm * a = nullptr;
-   OperatorPtr A;
   bool parallel;
 #ifdef MFEM_USE_MPI
   ParMesh * pmesh = nullptr;
@@ -270,8 +267,6 @@ public:
   void ResetFEM();
   void SetupFEM();

-   void UpdateEssentialTDofs();
-   void AssembleDiffusionBilinear(bool update_ess_tdofs=true);
   void Solve();
   GridFunction * GetFEMSolution();
   LinearForm * GetLinearForm() {return b;}
@@ -376,130 +371,6 @@ public:

 };

-/**
- * @brief Bregman projection of ρ = sigmoid(ψ) onto the subspace
- *        ∫_Ω ρ dx = θ vol(Ω) as follows:
- *
- *        1. Compute the root of the R → R function
- *            f(c) = ∫_Ω sigmoid(ψ + c) dx - θ vol(Ω)
- *           using the Illinois method
- *        2. Set ψ ← ψ + c.
- *
- * @param psi a GridFunction to be updated
- * @param alpha_grad alpha multiplied by gradient
- * @param target_volume θ vol(Ω)
- * @param tol Illinois iteration tolerance
- * @param max_its Illinois maximum iteration number
- * @return real_t Final volume (∫_Ω sigmoid(ψ) dx)
- */
-real_t proj(GridFunction &psi, GridFunction &alpha_grad, real_t target_volume,
-            real_t tol = 1e-12, int max_its = 100)
-{
-#ifdef MFEM_USE_MPI
-   FiniteElementSpace *fes = psi.FESpace();
-   ParFiniteElementSpace *pfes = dynamic_cast<ParFiniteElementSpace*>(fes);
-#endif
-   ConstantCoefficient zero_cf(0.0);
-   real_t a = -alpha_grad.ComputeMaxError(zero_cf);
-   real_t b = -a;
-   real_t y = 0.0;
-
-   MappedGridFunctionCoefficient sigmoid_psi(
-   &psi, [&y](const real_t x) { return sigmoid(x + y); });
-   std::unique_ptr<LinearForm> int_sigmoid_psi;
-#ifdef MFEM_USE_MPI
-   ParGridFunction *par_psi = dynamic_cast<ParGridFunction *>(&psi);
-   if (par_psi)
-   {
-      int_sigmoid_psi.reset(new ParLinearForm(par_psi->ParFESpace()));
-   }
-   else
-   {
-      int_sigmoid_psi.reset(new LinearForm(psi.FESpace()));
-   }
-#else
-   int_sigmoid_psi.reset(new LinearForm(psi.FESpace()));
-#endif
-   int_sigmoid_psi->AddDomainIntegrator(new DomainLFIntegrator(sigmoid_psi));
-
-   y = a;
-   int_sigmoid_psi->Assemble();
-   real_t f_a = int_sigmoid_psi->Sum();  // f_a := f(a) + θ vol(Ω)
-
-   y = b;
-   int_sigmoid_psi->Assemble();
-   real_t f_b = int_sigmoid_psi->Sum();  // f_b := f(b) + θ vol(Ω)
-#ifdef MFEM_USE_MPI
-   if (pfes)
-   {
-      MPI_Allreduce(MPI_IN_PLACE, &f_a, 1, MPITypeMap<real_t>::mpi_type,
-                    MPI_SUM, MPI_COMM_WORLD);
-      MPI_Allreduce(MPI_IN_PLACE, &f_b, 1, MPITypeMap<real_t>::mpi_type,
-                    MPI_SUM, MPI_COMM_WORLD);
-   }
-#endif
-   f_a -= target_volume; // f_a := f(a)
-   f_b -= target_volume; // f_b := f(b)
-   real_t c = 0.0;
-   real_t f_c = 0.0;
-   int side = 0;
-
-   bool done = false;
-   for (int k=0; k < max_its; k++)
-   {
-      c = (f_a * b - f_b * a) / (f_a - f_b);
-
-      if (abs(b - a) < tol * abs(b + a)) { done = true; break; }
-
-      y = c;
-      int_sigmoid_psi->Assemble();
-      f_c = int_sigmoid_psi->Sum(); // f_c := f(c) + θ vol(Ω)
-#ifdef MFEM_USE_MPI
-      if (pfes)
-      {
-         MPI_Allreduce(MPI_IN_PLACE, &f_c, 1, MPITypeMap<real_t>::mpi_type,
-                       MPI_SUM, MPI_COMM_WORLD);
-      }
-#endif
-      f_c -= target_volume; // f_c := f(c)
-
-      if (f_c * f_b > 0)
-      {
-         b = c;
-         f_b = f_c;
-         if (side == -1) { f_a /= 2.0; }
-         side = -1;
-      }
-      else if (f_c * f_a > 0)
-      {
-         a = c;
-         f_a = f_c;
-         if (side == 1) { f_b /= 2.0; }
-         side = 1;
-      }
-      else
-      {
-         done = true; break;
-      }
-   }
-   if (!done)
-   {
-      mfem_warning("Projection reached maximum iteration without converging. "
-                   "Result may not be accurate.");
-   }
-   y = 0.0;
-   psi += c;
-   int_sigmoid_psi->Assemble();
-   real_t material_volume = int_sigmoid_psi->Sum();
-#ifdef MFEM_USE_MPI
-   if (pfes)
-   {
-      MPI_Allreduce(MPI_IN_PLACE, &material_volume, 1,
-                    MPITypeMap<real_t>::mpi_type, MPI_SUM, MPI_COMM_WORLD);
-   }
-#endif
-   return material_volume;
-}

 // Poisson solver

@@ -551,8 +422,12 @@ void DiffusionSolver::SetupFEM()
   }
 }

-void DiffusionSolver::UpdateEssentialTDofs()
+void DiffusionSolver::Solve()
 {
+   OperatorPtr A;
+   Vector B, X;
+   Array<int> ess_tdof_list;
+
 #ifdef MFEM_USE_MPI
   if (parallel)
   {
@@ -565,39 +440,7 @@ void DiffusionSolver::UpdateEssentialTDofs()
 #else
   fes->GetEssentialTrueDofs(ess_bdr,ess_tdof_list);
 #endif
-}
-
-void DiffusionSolver::AssembleDiffusionBilinear(bool update_ess_tdofs)
-{
-   if (update_ess_tdofs)
-   {
-      UpdateEssentialTDofs();
-   }
-#ifdef MFEM_USE_MPI
-   if (parallel)
-   {
-      a = new ParBilinearForm(pfes);
-   }
-   else
-   {
-      a = new BilinearForm(fes);
-   }
-#else
-   a = new BilinearForm(fes);
-#endif
-   a->AddDomainIntegrator(new DiffusionIntegrator(*diffcf));
-   if (masscf)
-   {
-      a->AddDomainIntegrator(new MassIntegrator(*masscf));
-   }
-   a->Assemble();
-   a->FormSystemMatrix(ess_tdof_list, A);
-}
-
-void DiffusionSolver::Solve()
-{
-   Vector B, X;
-
+   *u=0.0;
   if (b)
   {
      delete b;
@@ -632,33 +475,31 @@ void DiffusionSolver::Solve()

   b->Assemble();

-   *u=0.0;
-   if (essbdr_cf)
-   {
-      u->ProjectBdrCoefficient(*essbdr_cf,ess_bdr);
-   }
+   BilinearForm * a = nullptr;

 #ifdef MFEM_USE_MPI
   if (parallel)
   {
-      X.SetSize(pfes->TrueVSize());
-      B.SetSize(pfes->TrueVSize());
-      dynamic_cast<ParGridFunction*>(u)->ParallelAssemble(X);
-      dynamic_cast<ParLinearForm*>(b)->ParallelAssemble(B);
-      dynamic_cast<ParBilinearForm*>(a)->ParallelEliminateTDofsInRHS(
-         ess_tdof_list, X, B);
+      a = new ParBilinearForm(pfes);
   }
   else
   {
-      X.NewDataAndSize(u->GetData(), u->Size());
-      B.NewDataAndSize(b->GetData(), b->Size());
-      a->EliminateVDofsInRHS(ess_tdof_list, X, B);
+      a = new BilinearForm(fes);
   }
 #else
-   X.NewDataAndSize(u->GetData(), u->Size());
-   B.NewDataAndSize(b->GetData(), b->Size());
-   a->EliminateVDofsInRHS(ess_tdof_list, X, B);
+   a = new BilinearForm(fes);
 #endif
+   a->AddDomainIntegrator(new DiffusionIntegrator(*diffcf));
+   if (masscf)
+   {
+      a->AddDomainIntegrator(new MassIntegrator(*masscf));
+   }
+   a->Assemble();
+   if (essbdr_cf)
+   {
+      u->ProjectBdrCoefficient(*essbdr_cf,ess_bdr);
+   }
+   a->FormLinearSystem(ess_tdof_list, *u, *b, A, X, B);

   CGSolver * cg = nullptr;
   Solver * M = nullptr;
@@ -687,6 +528,7 @@ void DiffusionSolver::Solve()
   delete M;
   delete cg;
   a->RecoverFEMSolution(X, *b, *u);
+   delete a;
 }

 GridFunction * DiffusionSolver::GetFEMSolution()
@@ -718,8 +560,6 @@ DiffusionSolver::~DiffusionSolver()
 #endif
   delete fec; fec = nullptr;
   delete b;
-   A.Clear();
-   delete a;
 }


@@ -4,8 +4,8 @@
 //
 // Sample runs:
 //    mpirun -np 4 ex37p -alpha 10 -pv
-//    mpirun -np 4 ex37p -lambda 0.1 -mu 0.1 -growth 1
-//    mpirun -np 4 ex37p -o 2 -alpha 10.0 -mi 50 -vf 0.4 -ntol 1e-5 -growth 1.5
+//    mpirun -np 4 ex37p -lambda 0.1 -mu 0.1
+//    mpirun -np 4 ex37p -o 2 -alpha 5.0 -mi 50 -vf 0.4 -ntol 1e-5
 //    mpirun -np 4 ex37p -r 6 -o 2 -alpha 10.0 -epsilon 0.02 -mi 50 -ntol 1e-5
 //
 // Description: This example code demonstrates the use of MFEM to solve a
@@ -54,6 +54,61 @@
 using namespace std;
 using namespace mfem;

+/**
+ * @brief Bregman projection of ρ = sigmoid(ψ) onto the subspace
+ *        ∫_Ω ρ dx = θ vol(Ω) as follows:
+ *
+ *        1. Compute the root of the R → R function
+ *            f(c) = ∫_Ω sigmoid(ψ + c) dx - θ vol(Ω)
+ *        2. Set ψ ← ψ + c.
+ *
+ * @param psi a GridFunction to be updated
+ * @param target_volume θ vol(Ω)
+ * @param tol Newton iteration tolerance
+ * @param max_its Newton maximum iteration number
+ * @return real_t Final volume, ∫_Ω sigmoid(ψ)
+ */
+real_t proj(ParGridFunction &psi, real_t target_volume, real_t tol=1e-12,
+            int max_its=10)
+{
+   MappedGridFunctionCoefficient sigmoid_psi(&psi, sigmoid);
+   MappedGridFunctionCoefficient der_sigmoid_psi(&psi, der_sigmoid);
+
+   ParLinearForm int_sigmoid_psi(psi.ParFESpace());
+   int_sigmoid_psi.AddDomainIntegrator(new DomainLFIntegrator(sigmoid_psi));
+   ParLinearForm int_der_sigmoid_psi(psi.ParFESpace());
+   int_der_sigmoid_psi.AddDomainIntegrator(new DomainLFIntegrator(
+                                              der_sigmoid_psi));
+   bool done = false;
+   for (int k=0; k<max_its; k++) // Newton iteration
+   {
+      int_sigmoid_psi.Assemble(); // Recompute f(c) with updated ψ
+      real_t f = int_sigmoid_psi.Sum();
+      MPI_Allreduce(MPI_IN_PLACE, &f, 1, MPITypeMap<real_t>::mpi_type,
+                    MPI_SUM, MPI_COMM_WORLD);
+      f -= target_volume;
+
+      int_der_sigmoid_psi.Assemble(); // Recompute df(c) with updated ψ
+      real_t df = int_der_sigmoid_psi.Sum();
+      MPI_Allreduce(MPI_IN_PLACE, &df, 1, MPITypeMap<real_t>::mpi_type,
+                    MPI_SUM, MPI_COMM_WORLD);
+
+      const real_t dc = -f/df;
+      psi += dc;
+      if (abs(dc) < tol) { done = true; break; }
+   }
+   if (!done)
+   {
+      mfem_warning("Projection reached maximum iteration without converging. "
+                   "Result may not be accurate.");
+   }
+   int_sigmoid_psi.Assemble();
+   real_t material_volume = int_sigmoid_psi.Sum();
+   MPI_Allreduce(MPI_IN_PLACE, &material_volume, 1,
+                 MPITypeMap<real_t>::mpi_type, MPI_SUM, MPI_COMM_WORLD);
+   return material_volume;
+}
+
 /*
 * ---------------------------------------------------------------
 *                      ALGORITHM PREAMBLE
@@ -138,11 +193,10 @@ int main(int argc, char *argv[])
   int ref_levels = 5;
   int order = 2;
   real_t alpha = 1.0;
-   real_t growth = 2;
   real_t epsilon = 0.01;
   real_t vol_fraction = 0.5;
   int max_it = 1e3;
-   real_t itol = 1e-2;
+   real_t itol = 1e-1;
   real_t ntol = 1e-4;
   real_t rho_min = 1e-6;
   real_t lambda = 1.0;
@@ -157,8 +211,6 @@ int main(int argc, char *argv[])
                  "Order (degree) of the finite elements.");
   args.AddOption(&alpha, "-alpha", "--alpha-step-length",
                  "Step length for gradient descent.");
-   args.AddOption(&growth, "-growth", "--alpha-growth-rate",
-                  "Growth rate of step length for gradient descent.");
   args.AddOption(&epsilon, "-epsilon", "--epsilon-thickness",
                  "Length scale for ρ.");
   args.AddOption(&max_it, "-mi", "--max-it",
@@ -307,7 +359,6 @@ int main(int argc, char *argv[])
   }
   FilterSolver->SetEssentialBoundary(ess_bdr_filter);
   FilterSolver->SetupFEM();
-   FilterSolver->AssembleDiffusionBilinear();

   ParBilinearForm mass(&control_fes);
   mass.AddDomainIntegrator(new InverseIntegrator(new MassIntegrator(one)));
@@ -361,7 +412,7 @@ int main(int argc, char *argv[])
   // 11. Iterate:
   for (int k = 1; k <= max_it; k++)
   {
-      if (k > 1) { alpha = std::pow((real_t) k,growth); }
+      if (k > 1) { alpha *= ((real_t) k) / ((real_t) k-1); }

      if (myid == 0)
      {
@@ -401,9 +452,7 @@ int main(int argc, char *argv[])

      // Step 5 - Update design variable ψ ← proj(ψ - αG)
      psi.Add(-alpha, grad);
-      ParGridFunction alpha_grad(grad);
-      alpha_grad *= alpha;
-      const real_t material_volume = proj(psi, alpha_grad, target_volume);
+      const real_t material_volume = proj(psi, target_volume);

      // Compute ||ρ - ρ_old|| in control fes.
      real_t norm_increment = zerogf.ComputeL1Error(succ_diff_rho);
@@ -433,16 +433,18 @@ int main(int argc, char *argv[])
   u.ProjectCoefficient(*u0);

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

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

   // Create data collection for solution output: either VisItDataCollection for
-   // ascii data files, or SidreDataCollection for binary data files.
+   // ascii data files, or ConduitDataCollection for binary data files.
   DataCollection *dc = NULL;
   if (visit)
   {
      if (binary)
      {
-#ifdef MFEM_USE_SIDRE
-         dc = new SidreDataCollection("Example9-Parallel", pmesh);
+#ifdef MFEM_USE_CONDUIT
+         auto conduit_dc = new ConduitDataCollection("Example9-Parallel", pmesh);
+         conduit_dc->SetProtocol("hdf5");
+         dc = conduit_dc;
 #else
-         MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
+         MFEM_ABORT("Must build with MFEM_USE_CONDUIT=YES for binary output.");
 #endif
      }
      else
@@ -76,4 +76,4 @@ clean-build:
 	rm -rf *.dSYM *.TVD.*breakpoints

 clean-exec:
-	@rm -f refined.mesh sol.gf mesh.* sol.*
+	@rm -f refined.mesh sol.gf
@@ -354,16 +354,18 @@ int main(int argc, char *argv[])
   }

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

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

 SUBDIRS_ALL = $(addsuffix /all,$(SUBDIRS))
 SUBDIRS_TEST = $(addsuffix /test,$(SUBDIRS))
-SUBDIRS_TEST_NOCLEAN = $(addsuffix /test-noclean,$(SUBDIRS))
 SUBDIRS_CLEAN = $(addsuffix /clean,$(SUBDIRS))
 SUBDIRS_TPRINT = $(addsuffix /test-print,$(SUBDIRS))

@@ -88,9 +87,8 @@ SUBDIRS_TPRINT = $(addsuffix /test-print,$(SUBDIRS))

 all: $(EXAMPLES) $(SUBDIRS_ALL)

-.PHONY: $(SUBDIRS_ALL) $(SUBDIRS_TEST) $(SUBDIRS_TEST_NOCLEAN) \
-   $(SUBDIRS_CLEAN) $(SUBDIRS_TPRINT)
-$(SUBDIRS_ALL) $(SUBDIRS_TEST) $(SUBDIRS_TEST_NOCLEAN) $(SUBDIRS_CLEAN):
+.PHONY: $(SUBDIRS_ALL) $(SUBDIRS_TEST) $(SUBDIRS_CLEAN) $(SUBDIRS_TPRINT)
+$(SUBDIRS_ALL) $(SUBDIRS_TEST) $(SUBDIRS_CLEAN):
 	$(MAKE) -C $(@D) $(@F)
 $(SUBDIRS_TPRINT):
 	@$(MAKE) -C $(@D) $(@F)
@@ -109,7 +107,6 @@ endif
 MFEM_TESTS = EXAMPLES
 include $(MFEM_TEST_MK)
 test: $(SUBDIRS_TEST)
-test-noclean: $(SUBDIRS_TEST_NOCLEAN)
 test-print: $(SUBDIRS_TPRINT)

 # Testing: Parallel vs. serial runs
@@ -160,8 +157,6 @@ ex37-test-seq: ex37
 	@$(call mfem-test,$<,, Serial example,-mi 3)
 ex37p-test-par: ex37p
 	@$(call mfem-test,$<, $(RUN_MPI), Parallel example,-mi 3)
-ex39-test-seq: ex39
-	@$(call mfem-test,$<,, Serial example,-m ../data/compass.mesh)
 ex41-test-seq: ex41
 	@$(call mfem-test,$<,, Serial example,-tf 1.0)
 ex41p-test-par: ex41p
@@ -368,16 +368,18 @@ int main(int argc, char *argv[])
   }

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

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

   // Create data collection for solution output: either VisItDataCollection for
-   // ascii data files, or SidreDataCollection for binary data files.
+   // ascii data files, or ConduitDataCollection for binary data files.
   DataCollection *dc = NULL;
   if (visit)
   {
      if (binary)
      {
-#ifdef MFEM_USE_SIDRE
-         dc = new SidreDataCollection("Example9-Parallel", pmesh);
+#ifdef MFEM_USE_CONDUIT
+         auto conduit_dc = new ConduitDataCollection("Example9-Parallel", pmesh);
+         conduit_dc->SetProtocol("hdf5");
+         dc = conduit_dc;
 #else
-         MFEM_ABORT("Must build with MFEM_USE_SIDRE=YES for binary output.");
+         MFEM_ABORT("Must build with MFEM_USE_CONDUIT=YES for binary output.");
 #endif
      }
      else
@@ -57,8 +57,6 @@ set(SRCS
  integ/lininteg_domain_grad.cpp
  integ/lininteg_domain_vectorfe.cpp
  integ/nonlininteg_vecconvection_pa.cpp
-  integ/nonlininteg_vecconvection_pa_diag.cpp
-  integ/nonlininteg_vecconvection_pa_grad.cpp
  integ/nonlininteg_vecconvection_mf.cpp
  coefficient.cpp
  complex_fem.cpp
@@ -135,7 +133,7 @@ set(SRCS
  tmop/assemble/diag2.cpp
  tmop/assemble/grad2_limit.cpp
  tmop/assemble/grad2.cpp
-  tmop/assemble/diag3_limit.cpp
+  tmop/assemble/diag3_limit.cpp  
  tmop/assemble/diag3.cpp
  tmop/assemble/grad3_limit.cpp
  tmop/assemble/grad3.cpp
@@ -173,12 +171,8 @@ set(SRCS
  tmop_tools.cpp
  tmop_amr.cpp
  gslib.cpp
-  gslib/findptsedge_local_2.cpp
-  gslib/findptsedge_local_3.cpp
-  gslib/findptssurf_local_3.cpp
  gslib/findpts_local_2.cpp
  gslib/findpts_local_3.cpp
-  gslib/interpolate_local_1.cpp
  gslib/interpolate_local_2.cpp
  gslib/interpolate_local_3.cpp
  transfer.cpp
@@ -197,20 +191,14 @@ set(HDRS
  integ/bilininteg_dgtrace_kernels.hpp
  integ/bilininteg_vecdiffusion_kernels.hpp
  integ/bilininteg_convection_kernels.hpp
-  integ/bilininteg_diffusion_pa_simplices.hpp
  integ/bilininteg_diffusion_kernels.hpp
  integ/bilininteg_elasticity_kernels.hpp
  integ/bilininteg_hcurl_kernels.hpp
  integ/bilininteg_hdiv_kernels.hpp
  integ/bilininteg_hcurlhdiv_kernels.hpp
  integ/bilininteg_mass_kernels.hpp
-  integ/bilininteg_mass_pa_simplices.hpp
  integ/bilininteg_vecdiffusion_pa.hpp
-  integ/bilininteg_vecdiv_pa.hpp
  integ/bilininteg_vecmass_pa.hpp
-  integ/nonlininteg_vecconvection_pa.hpp
-  integ/nonlininteg_vecconvection_pa_diag.hpp
-  integ/nonlininteg_vecconvection_pa_grad.hpp
  coefficient.hpp
  complex_fem.hpp
  convergence.hpp
@@ -317,7 +305,6 @@ set(HDRS
  tmop_tools.hpp
  tmop_amr.hpp
  gslib.hpp
-  gslib/gslib_kernel_helpers.hpp
  transfer.hpp
  hyperbolic.hpp
  integrator.hpp
@@ -325,11 +312,6 @@ set(HDRS
  particleset.hpp
  )

-if (MFEM_USE_SIDRE)
-  list(APPEND SRCS sidredatacollection.cpp)
-  list(APPEND HDRS sidredatacollection.hpp)
-endif()
-
 if (MFEM_USE_CONDUIT)
  list(APPEND SRCS conduitdatacollection.cpp)
  list(APPEND HDRS conduitdatacollection.hpp)
@@ -1345,8 +1345,7 @@ real_t DiffusionIntegrator::ComputeFluxEnergy
 }

 const IntegrationRule &DiffusionIntegrator::GetRule(
-   const FiniteElement &trial_fe, const FiniteElement &test_fe,
-   const bool stroud)
+   const FiniteElement &trial_fe, const FiniteElement &test_fe)
 {
   int order;
   if (trial_fe.Space() == FunctionSpace::Pk)
@@ -1363,15 +1362,7 @@ const IntegrationRule &DiffusionIntegrator::GetRule(
   {
      return RefinedIntRules.Get(trial_fe.GetGeomType(), order);
   }
-
-   if (stroud)
-   {
-      return StroudIntRules.Get(trial_fe.GetGeomType(), order);
-   }
-   else
-   {
-      return IntRules.Get(trial_fe.GetGeomType(), order);
-   }
+   return IntRules.Get(trial_fe.GetGeomType(), order);
 }

 MassIntegrator::MassIntegrator(const IntegrationRule *ir)
@@ -1458,8 +1449,7 @@ void MassIntegrator::AssembleElementMatrix2(

 const IntegrationRule &MassIntegrator::GetRule(const FiniteElement &trial_fe,
                                               const FiniteElement &test_fe,
-                                               const ElementTransformation &Trans,
-                                               const bool stroud)
+                                               const ElementTransformation &Trans)
 {
   // int order = trial_fe.GetOrder() + test_fe.GetOrder();
   const int order = trial_fe.GetOrder() + test_fe.GetOrder() + Trans.OrderW();
@@ -1468,15 +1458,7 @@ const IntegrationRule &MassIntegrator::GetRule(const FiniteElement &trial_fe,
   {
      return RefinedIntRules.Get(trial_fe.GetGeomType(), order);
   }
-
-   if (stroud)
-   {
-      return StroudIntRules.Get(trial_fe.GetGeomType(), order);
-   }
-   else
-   {
-      return IntRules.Get(trial_fe.GetGeomType(), order);
-   }
+   return IntRules.Get(trial_fe.GetGeomType(), order);
 }


@@ -2184,22 +2184,11 @@ public:
                                   const Vector&, const Vector&,
                                   Vector&, const int, const int);

-   using ApplySimplexKernelType = void(*)(const int, const bool, const Array<int>&,
-                                          const Array<int>&,
-                                          const Array<int>&, const Array<int>&, const Array<int>&,
-                                          const Array<real_t>&, const Array<real_t>&,
-                                          const Array<real_t>&, const Array<real_t>&,
-                                          const Array<real_t>&, const Array<real_t>&,
-                                          const Vector&, const Vector&,
-                                          Vector&, const int, const int);
-
   using DiagonalKernelType = void(*)(const int, const bool, const Array<real_t>&,
                                      const Array<real_t>&, const Vector&, Vector&,
                                      const int, const int);

   MFEM_REGISTER_KERNELS(ApplyPAKernels, ApplyKernelType, (int, int, int));
-   MFEM_REGISTER_KERNELS(ApplySimplexPAKernels, ApplySimplexKernelType, (int, int,
-                                                                         int));
   MFEM_REGISTER_KERNELS(DiagonalPAKernels, DiagonalKernelType, (int, int, int));
   struct Kernels { Kernels(); };

@@ -2352,8 +2341,7 @@ public:
   void AddMultPatchPA(const int patch, const Vector &x, Vector &y) const;

   static const IntegrationRule &GetRule(const FiniteElement &trial_fe,
-                                         const FiniteElement &test_fe,
-                                         const bool stroud = false);
+                                         const FiniteElement &test_fe);

   bool SupportsCeed() const override { return DeviceCanUseCeed(); }

@@ -2364,13 +2352,6 @@ public:
   {
      ApplyPAKernels::Specialization<DIM,D1D,Q1D>::Add();
      DiagonalPAKernels::Specialization<DIM,D1D,Q1D>::Add();
-      AddSimplexSpecialization<DIM,D1D,Q1D>();
-   }
-
-   template <int DIM, int D1D, int Q1D>
-   static void AddSimplexSpecialization()
-   {
-      ApplySimplexPAKernels::Specialization<DIM,D1D,Q1D>::Add();
   }
 protected:
   const IntegrationRule* GetDefaultIntegrationRule(
@@ -2407,22 +2388,11 @@ public:
                                   const Array<real_t>&, const Vector&,
                                   const Vector&, Vector&, const int, const int);

-   using ApplySimplexKernelType = void(*)(const int, const Array<int>&,
-                                          const Array<int>&,
-                                          const Array<int>&, const Array<int>&, const Array<int>&,
-                                          const Array<real_t>&, const Array<real_t>&,
-                                          const Array<real_t>&, const Array<real_t>&,
-                                          const Array<real_t>&, const Array<real_t>&,
-                                          const Vector&, const Vector&, Vector&,
-                                          const int, const int);
-
   using DiagonalKernelType =  void(*)(const int, const Array<real_t>&,
                                       const Vector&, Vector&, const int,
                                       const int);

   MFEM_REGISTER_KERNELS(ApplyPAKernels, ApplyKernelType, (int, int, int));
-   MFEM_REGISTER_KERNELS(ApplySimplexPAKernels, ApplySimplexKernelType, (int, int,
-                                                                         int));
   MFEM_REGISTER_KERNELS(DiagonalPAKernels, DiagonalKernelType, (int, int, int));
   struct Kernels { Kernels(); };

@@ -2471,8 +2441,7 @@ public:

   static const IntegrationRule &GetRule(const FiniteElement &trial_fe,
                                         const FiniteElement &test_fe,
-                                         const ElementTransformation &Trans,
-                                         const bool stroud = false);
+                                         const ElementTransformation &Trans);

   bool SupportsCeed() const override { return DeviceCanUseCeed(); }

@@ -2483,13 +2452,6 @@ public:
   {
      ApplyPAKernels::Specialization<DIM,D1D,Q1D>::Add();
      DiagonalPAKernels::Specialization<DIM,D1D,Q1D>::Add();
-      AddSimplexSpecialization<DIM,D1D,Q1D>();
-   }
-
-   template <int DIM, int D1D, int Q1D>
-   static void AddSimplexSpecialization()
-   {
-      ApplySimplexPAKernels::Specialization<DIM,D1D,Q1D>::Add();
   }

 protected:
@@ -2689,22 +2651,14 @@ public:
   void AddMultMF(const Vector &x, Vector &y) const override;
   bool SupportsCeed() const override { return DeviceCanUseCeed(); }

-   // PA AddMultPA kernels
   using VectorMassAddMultPAType =
      void(*)(const int,  const int,
              const Array<real_t>&, const Vector&,
              const Vector&, Vector&, const int, const int);
+
   MFEM_REGISTER_KERNELS(VectorMassAddMultPA,
                         VectorMassAddMultPAType,
                         (int, int, int));
-
-   // PA DiagonalPA kernels
-   using VectorMassAssembleDiagonalPAType =
-      void(*)(const int, const int, const int,
-              const real_t*, const real_t*, real_t*);
-   MFEM_REGISTER_KERNELS(VectorMassAssembleDiagonalPA,
-                         VectorMassAssembleDiagonalPAType,
-                         (int /*dim*/, int /*q1d*/));
 };


@@ -3003,10 +2957,11 @@ public:
    vector (diagonal matrix), or matrix), trial function $u$ is in $H(curl$ or
    $H(div)$, and test function $v$ is in $H(curl$, $H(div)$, or $v=(v_1,\dots,v_n)$, where
    $v_i$ are in $H^1$. */
-class VectorFEMassIntegrator : public BilinearFormIntegrator
+class VectorFEMassIntegrator: public BilinearFormIntegrator
 {
 private:
-   void Init(Coefficient *q, DiagonalMatrixCoefficient *dq, MatrixCoefficient *mq);
+   void Init(Coefficient *q, DiagonalMatrixCoefficient *dq, MatrixCoefficient *mq)
+   { Q = q; DQ = dq; MQ = mq; }

 #ifndef MFEM_THREAD_SAFE
   Vector shape;
@@ -3029,8 +2984,7 @@ protected:
   const DofToQuad *mapsOtest;     ///< Not owned. DOF-to-quad map, open.
   const DofToQuad *mapsCtest;     ///< Not owned. DOF-to-quad map, closed.
   const GeometricFactors *geom;   ///< Not owned
-   int dim, ne, nq, dofs1D, dofs1Dtest, quad1D;
-   FiniteElement::DerivType trial_fetype, test_fetype;
+   int dim, ne, nq, dofs1D, dofs1Dtest, quad1D, trial_fetype, test_fetype;
   bool symmetric = true; ///< False if using a nonsymmetric matrix coefficient

 public:
@@ -3061,29 +3015,6 @@ public:
                   const bool add) override;

   const Coefficient *GetCoefficient() const { return Q; }
-
-   using ApplyKernelType =
-      void (*)(const int NE, bool symmetric, const bool scalar_coeff,
-               const Array<real_t> &trialBO, const Array<real_t> &trialBC,
-               const Array<real_t> &testBOt, const Array<real_t> &testBCt,
-               const Vector &pa_data, const Vector &x, Vector &y,
-               const int triald1d, const int testd1d, const int q1d);
-
-   /// parameters: trial_fetype, test_fetype, ndims, trial_d1d, test_d1d, q1d
-   MFEM_REGISTER_KERNELS(ApplyPAKernels, ApplyKernelType,
-                         (FiniteElement::DerivType, FiniteElement::DerivType,
-                          int, int, int, int));
-
-   struct Kernels { Kernels(); };
-
-   template <FiniteElement::DerivType TrialType,
-             FiniteElement::DerivType TestType, int DIM, int TRIAL_D1D,
-             int TEST_D1D, int Q1D>
-   static void AddSpecialization()
-   {
-      ApplyPAKernels::Specialization<TrialType, TestType, DIM, TRIAL_D1D,
-                     TEST_D1D, Q1D>::Add();
-   }
 };

 /** Integrator for $(Q \nabla \cdot u, v)$ where $u=(u_1,\cdots,u_n)$ and all $u_i$ are in the same
@@ -3129,24 +3060,6 @@ public:
   void AddMultPA(const Vector &x, Vector &y) const override;
   void AddMultTransposePA(const Vector &x, Vector &y) const override;

-   using VectorDivergenceAddMultPAType =
-      void (*)(const int ne,
-               const Array<real_t> &b, const Array<real_t> &g, const Array<real_t> &bt,
-               const Vector &op, const Vector &x, Vector &y,
-               const int tr_d1d, const int te_d1d, const int q1d);
-   MFEM_REGISTER_KERNELS(VectorDivergenceAddMultPA,
-                         VectorDivergenceAddMultPAType,
-                         (int, int, int, int));
-
-   using VectorDivergenceAddMultTransposePAType =
-      void (*)(const int ne,
-               const Array<real_t> &bt, const Array<real_t> &gt, const Array<real_t> &b,
-               const Vector &q, const Vector &x, Vector &y,
-               const int tr_d1d, const int te_d1d, const int q1d);
-   MFEM_REGISTER_KERNELS(VectorDivergenceAddMultTransposePA,
-                         VectorDivergenceAddMultTransposePAType,
-                         (int, int, int, int));
-
   static const IntegrationRule &GetRule(const FiniteElement &trial_fe,
                                         const FiniteElement &test_fe,
                                         const ElementTransformation &Trans);
@@ -54,8 +54,6 @@ void Coefficient::Project(QuadratureFunction &qf)
   QuadratureSpaceBase &qspace = *qf.GetSpace();
   const int ne = qspace.GetNE();
   Vector values;
-   // GetValues makes a reference, but we need it to be valid on Host
-   qf.HostWrite();
   for (int iel = 0; iel < ne; ++iel)
   {
      qf.GetValues(iel, values);
@@ -329,8 +327,6 @@ void VectorCoefficient::Project(QuadratureFunction &qf)
   const int ne = qspace.GetNE();
   DenseMatrix values;
   Vector col;
-   // GetValues makes a reference, but we need it to be valid on Host
-   qf.HostWrite();
   for (int iel = 0; iel < ne; ++iel)
   {
      qf.GetValues(iel, values);
@@ -699,8 +695,6 @@ void MatrixCoefficient::Project(QuadratureFunction &qf, bool transpose)
   QuadratureSpaceBase &qspace = *qf.GetSpace();
   const int ne = qspace.GetNE();
   DenseMatrix values, matrix;
-   // GetValues makes a reference, but we need it to be valid on Host
-   qf.HostWrite();
   for (int iel = 0; iel < ne; ++iel)
   {
      qf.GetValues(iel, values);
@@ -52,9 +52,6 @@ public:
   /// Get the time for time dependent coefficients
   real_t GetTime() { return time; }

-   /// Returns dimension of the vector.
-   int GetVDim() { return 1; }
-
   /** @brief Evaluate the coefficient in the element described by @a T at the
       point @a ip. */
   /** @note When this method is called, the caller must make sure that the
@@ -166,75 +166,6 @@ public:
      return sqrt(err_r * err_r + err_i * err_i);
   }

-   /// @brief Returns Max|u_ex - u_h| error for complex-valued H1 or L2 elements
-   ///
-   /// Compute the $L_\infty$ error across the entire domain.
-   ///
-   /// @param[in] exsolr  Coefficient object reproducing the real part of the
-   ///                    anticipated values of the scalar field, Re(u_ex).
-   /// @param[in] exsoli  Coefficient object reproducing the imaginary part of
-   ///                    the anticipated values of the scalar field, Im(u_ex).
-   /// @param[in] irs        Optional pointer to an array of custom integration
-   ///                       rules e.g. higher order than the default rules. If
-   ///                       present the array will be indexed by
-   ///                       Geometry::Type.
-   ///
-   /// @note Uses ComputeLpError internally. See the ComputeLpError
-   ///       documentation for generalizations of this error computation.
-   ///
-   /// @note If an array of integration rules is provided through @a irs, be
-   ///       sure to include valid rules for each element type that may occur
-   ///       in the list of elements.
-   ///
-   virtual real_t ComputeMaxError(Coefficient &exsolr,
-                                  Coefficient &exsoli,
-                                  const IntegrationRule *irs[] = NULL) const
-   {
-      return ComputeLpError(infinity(), exsolr, exsoli, NULL, irs);
-   }
-
-   /// @brief Returns ||u_ex - u_h||_Lp for complex-valued H1 or L2 elements
-   ///
-   /// Computes:
-   ///    $$(\sum_{elems} \int_{elem} w \, |u_{ex} - u_h|^p)^{1/p}$$
-   /// Where:
-   ///    $$|u_{ex} - u_h| = \sqrt{Re(u_{ex} - u_h)^2 + Im(u_{ex} - u_h)^2}$$
-   ///
-   /// @param[in] p       Real value indicating the exponent of the $L^p$ norm.
-   ///                    To avoid domain errors p should have a positive value,
-   ///                    either finite or infinite.
-   /// @param[in] exsolr  Coefficient object reproducing the real part of the
-   ///                    anticipated values of the scalar field, Re(u_ex).
-   /// @param[in] exsoli  Coefficient object reproducing the imaginary part of
-   ///                    the anticipated values of the scalar field, Im(u_ex).
-   /// @param[in] weight  Optional pointer to a Coefficient object reproducing
-   ///                    a weighting function, w.
-   /// @param[in] irs     Optional pointer to an array of custom integration
-   ///                    rules e.g. higher order than the default rules. If
-   ///                    present the array will be indexed by Geometry::Type.
-   /// @param[in] elems   Optional pointer to a marker array, with a length
-   ///                    equal to the number of local elements, indicating
-   ///                    which elements to integrate over. Only those elements
-   ///                    corresponding to non-zero entries in @a elems will
-   ///                    contribute to the computed L2 error.
-   ///
-   /// @note If an array of integration rules is provided through @a irs, be
-   ///       sure to include valid rules for each element type that may occur
-   ///       in the list of elements.
-   ///
-   /// @note Quadratures with negative weights (as in some simplex integration
-   ///       rules in MFEM) can produce negative integrals even with
-   ///       non-negative integrands. To avoid returning negative errors this
-   ///       function uses the absolute values of the element-wise integrals.
-   ///       This may lead to results which are not entirely consistent with
-   ///       such integration rules.
-   virtual real_t ComputeLpError(const real_t p,
-                                 Coefficient &exsolr,
-                                 Coefficient &exsoli,
-                                 Coefficient *weight = NULL,
-                                 const IntegrationRule *irs[] = NULL,
-                                 const Array<int> *elems = NULL) const;
-
   /// Save the ComplexGridFunction to an output stream.
   virtual void Save(std::ostream &out) const;

@@ -505,186 +436,6 @@ public:
   virtual ~SesquilinearForm();
 };

-/** Class for a mixed sesquilinear form
-
-    A mixed sesquilinear form is a generalization of a mixed bilinear form to
-    complex-valued fields. Mixed sesquilinear forms are linear in the second
-    argument but the first argument involves a complex conjugate in the sense
-    that:
-
-                a(alpha u, beta v) = conj(alpha) beta a(u, v)
-
-    The @a convention argument in the class's constructor is documented in the
-    mfem::ComplexOperator class found in linalg/complex_operator.hpp.
-
-    When supplying integrators to the MixedSesquilinearForm either the real or
-    imaginary integrator can be NULL. This indicates that the corresponding
-    portion of the complex-valued material coefficient is equal to zero.
-*/
-class MixedSesquilinearForm
-{
-private:
-   ComplexOperator::Convention conv;
-
-   MixedBilinearForm * mblfr;
-   MixedBilinearForm * mblfi;
-
-   /* These methods check if the real/imag parts of the sesqulinear form are not
-      empty */
-   bool RealInteg();
-   bool ImagInteg();
-
-public:
-   MixedSesquilinearForm(
-      FiniteElementSpace * trial_fes,
-      FiniteElementSpace * test_fes,
-      ComplexOperator::Convention convention = ComplexOperator::HERMITIAN);
-
-   /** @brief Create a MixedSesquilinearForm on the given trial and test
-       FiniteElementSpaces, using the same integrators as the
-       MixedBilinearForms @a bfr and @a bfi.
-
-       The FiniteElementSpace pointers are not owned by the newly constructed
-       object.
-
-       The integrators are copied as pointers and they are not owned by the
-       newly constructed MixedSesquilinearForm. */
-   MixedSesquilinearForm(
-      FiniteElementSpace * trial_fes,
-      FiniteElementSpace * test_fes,
-      MixedBilinearForm * bfr,
-      MixedBilinearForm * bfi,
-      ComplexOperator::Convention convention = ComplexOperator::HERMITIAN);
-
-   ComplexOperator::Convention GetConvention() const { return conv; }
-   void SetConvention(const ComplexOperator::Convention & convention) { conv = convention; }
-
-   /// Set the desired assembly level.
-   /** Valid choices are:
-
-       - AssemblyLevel::LEGACY (default)
-       - AssemblyLevel::FULL
-       - AssemblyLevel::PARTIAL
-       - AssemblyLevel::ELEMENT
-       - AssemblyLevel::NONE
-
-       This method must be called before assembly. */
-   void SetAssemblyLevel(AssemblyLevel assembly_level)
-   {
-      mblfr->SetAssemblyLevel(assembly_level);
-      mblfi->SetAssemblyLevel(assembly_level);
-   }
-
-   MixedBilinearForm & real() { return *mblfr; }
-   MixedBilinearForm & imag() { return *mblfi; }
-   const MixedBilinearForm & real() const { return *mblfr; }
-   const MixedBilinearForm & imag() const { return *mblfi; }
-
-   /// Adds new Domain Integrator.
-   void AddDomainIntegrator(BilinearFormIntegrator * bfi_real,
-                            BilinearFormIntegrator * bfi_imag);
-
-   /// Adds new Domain Integrator, restricted to specific attributes.
-   void AddDomainIntegrator(BilinearFormIntegrator * bfi_real,
-                            BilinearFormIntegrator * bfi_imag,
-                            Array<int> & elem_marker);
-
-   /// Adds new Boundary Integrator.
-   void AddBoundaryIntegrator(BilinearFormIntegrator * bfi_real,
-                              BilinearFormIntegrator * bfi_imag);
-
-   /** @brief Adds new boundary Integrator, restricted to specific boundary
-       attributes.
-
-       Assumes ownership of @a bfi.
-
-       The mfem::array @a bdr_marker is stored internally as a pointer to the given
-       mfem::Array<int> object. */
-   void AddBoundaryIntegrator(BilinearFormIntegrator * bfi_real,
-                              BilinearFormIntegrator * bfi_imag,
-                              Array<int> & bdr_marker);
-
-   /// Adds new interior Face Integrator. Assumes ownership of @a bfi.
-   void AddInteriorFaceIntegrator(BilinearFormIntegrator * bfi_real,
-                                  BilinearFormIntegrator * bfi_imag);
-
-   /// Adds new boundary Face Integrator. Assumes ownership of @a bfi.
-   void AddBdrFaceIntegrator(BilinearFormIntegrator * bfi_real,
-                             BilinearFormIntegrator * bfi_imag);
-
-   /** @brief Adds new boundary Face Integrator, restricted to specific boundary
-       attributes.
-
-       Assumes ownership of @a bfi.
-
-       The mfem::array @a bdr_marker is stored internally as a pointer to the given
-       mfem::Array<int> object. */
-   void AddBdrFaceIntegrator(BilinearFormIntegrator * bfi_real,
-                             BilinearFormIntegrator * bfi_imag,
-                             Array<int> & bdr_marker);
-
-   /** @brief Add a trace face integrator. Assumes ownership of @a bfi.
-
-       This type of integrator assembles terms over all faces of the mesh using
-       the face FE from the trial space and the two adjacent volume FEs from
-       the test space. */
-   void AddTraceFaceIntegrator(BilinearFormIntegrator * bfi_real,
-                               BilinearFormIntegrator * bfi_imag);
-
-   /// Adds a boundary trace face integrator. Assumes ownership of @a bfi.
-   void AddBdrTraceFaceIntegrator(BilinearFormIntegrator * bfi_real,
-                                  BilinearFormIntegrator * bfi_imag);
-
-   /// Adds a boundary trace face integrator. Assumes ownership of @a bfi.
-   void AddBdrTraceFaceIntegrator(BilinearFormIntegrator * bfi_real,
-                                  BilinearFormIntegrator * bfi_imag,
-                                  Array<int> &bdr_marker);
-
-   /// Assemble the local matrix
-   void Assemble(int skip_zeros = 1);
-
-   /// Finalizes the matrix initialization.
-   void Finalize(int skip_zeros = 1);
-
-   /// Updates the internal mixed forms with the new finite element space.
-   virtual void Update();
-
-   /** @brief Return a ComplexSparseMatrix wrapping the local (L-dof) real
-       and imaginary matrices of the form.
-
-       The returned wrapper has to be deleted by the caller, but it does not
-       own the wrapped real and imaginary matrices, which remain owned by
-       this form. */
-   ComplexSparseMatrix *AssembleComplexSparseMatrix();
-
-   /// Return the trial FE space associated with the MixedSesquilinearForm.
-   FiniteElementSpace *TrialFESpace() { return mblfr->TrialFESpace(); }
-
-   /// Read-only access to the associated trial FiniteElementSpace.
-   const FiniteElementSpace *TrialFESpace() const { return mblfr->TrialFESpace(); }
-
-   /// Return the test FE space associated with the MixedSesquilinearForm.
-   FiniteElementSpace *TestFESpace() { return mblfr->TestFESpace(); }
-
-   /// Read-only access to the associated test FiniteElementSpace.
-   const FiniteElementSpace *TestFESpace() const { return mblfr->TestFESpace(); }
-
-
-   void FormRectangularLinearSystem(const Array<int> & ess_trial_tdof_list,
-                                    const Array<int> & ess_test_tdof_list,
-                                    Vector & x,
-                                    Vector & b,
-                                    OperatorHandle & A,
-                                    Vector & X,
-                                    Vector & B);
-
-   void FormRectangularSystemMatrix(const Array<int> & ess_trial_tdof_list,
-                                    const Array<int> & ess_test_tdof_list,
-                                    OperatorHandle & A);
-
-   virtual ~MixedSesquilinearForm();
-};
-
 #ifdef MFEM_USE_MPI

 /// Class for parallel complex-valued grid function - real + imaginary part
@@ -1101,169 +852,6 @@ public:
   virtual ~ParSesquilinearForm();
 };

-/** Class for a parallel mixed sesquilinear form
-
-    A mixed sesquilinear form is a generalization of a mixed bilinear form to
-    complex-valued fields. Mixed sesquilinear forms are linear in the second
-    argument but the first argument involves a complex conjugate in the sense
-    that:
-
-                a(alpha u, beta v) = conj(alpha) beta a(u, v)
-
-    The @a convention argument in the class's constructor is documented in the
-    mfem::ComplexOperator class found in linalg/complex_operator.hpp.
-
-    When supplying integrators to the ParMixedSesquilinearForm either the real
-    or imaginary integrator can be NULL. This indicates that the corresponding
-    portion of the complex-valued material coefficient is equal to zero.
-*/
-class ParMixedSesquilinearForm
-{
-private:
-   ComplexOperator::Convention conv;
-
-   ParMixedBilinearForm * pmblfr;
-   ParMixedBilinearForm * pmblfi;
-
-   /* These methods check if the real/imag parts of the sesqulinear form are
-    not empty */
-   bool RealInteg();
-   bool ImagInteg();
-
-public:
-   ParMixedSesquilinearForm(
-      ParFiniteElementSpace * trial_fes,
-      ParFiniteElementSpace * test_fes,
-      ComplexOperator::Convention convention = ComplexOperator::HERMITIAN);
-
-   /** @brief Create a ParMixedSesquilinearForm on the given trial and test
-       ParFiniteElementSpaces, using the same integrators as the
-       ParMixedBilinearForms @a pbfr and @a pbfi.
-
-       The ParFiniteElementSpace pointers are not owned by the newly
-       constructed object.
-
-       The integrators are copied as pointers and they are not owned by the
-       newly constructed ParMixedSesquilinearForm. */
-   ParMixedSesquilinearForm(
-      ParFiniteElementSpace * trial_fes,
-      ParFiniteElementSpace * test_fes,
-      ParMixedBilinearForm * pbfr,
-      ParMixedBilinearForm * pbfi,
-      ComplexOperator::Convention convention = ComplexOperator::HERMITIAN);
-
-   ComplexOperator::Convention GetConvention() const { return conv; }
-   void SetConvention(const ComplexOperator::Convention & convention) { conv = convention; }
-
-   /// Set the desired assembly level.
-   /** Valid choices are:
-
-       - AssemblyLevel::LEGACY (default)
-       - AssemblyLevel::FULL
-       - AssemblyLevel::PARTIAL
-       - AssemblyLevel::ELEMENT
-       - AssemblyLevel::NONE
-
-       This method must be called before assembly. */
-   void SetAssemblyLevel(AssemblyLevel assembly_level)
-   {
-      pmblfr->SetAssemblyLevel(assembly_level);
-      pmblfi->SetAssemblyLevel(assembly_level);
-   }
-
-   ParMixedBilinearForm & real() { return *pmblfr; }
-   ParMixedBilinearForm & imag() { return *pmblfi; }
-   const ParMixedBilinearForm & real() const { return *pmblfr; }
-   const ParMixedBilinearForm & imag() const { return *pmblfi; }
-
-   /// Adds new Domain Integrator.
-   void AddDomainIntegrator(BilinearFormIntegrator * bfi_real,
-                            BilinearFormIntegrator * bfi_imag);
-
-   /// Adds new Domain Integrator, restricted to specific attributes.
-   void AddDomainIntegrator(BilinearFormIntegrator * bfi_real,
-                            BilinearFormIntegrator * bfi_imag,
-                            Array<int> & elem_marker);
-
-   /// Adds new Boundary Integrator.
-   void AddBoundaryIntegrator(BilinearFormIntegrator * bfi_real,
-                              BilinearFormIntegrator * bfi_imag);
-
-   /** @brief Adds new boundary Integrator, restricted to specific boundary
-       attributes.
-
-       Assumes ownership of @a bfi.
-
-       The mfem::array @a bdr_marker is stored internally as a pointer to the given
-       mfem::Array<int> object. */
-   void AddBoundaryIntegrator(BilinearFormIntegrator * bfi_real,
-                              BilinearFormIntegrator * bfi_imag,
-                              Array<int> & bdr_marker);
-
-   /// Adds new interior Face Integrator. Assumes ownership of @a bfi.
-   void AddInteriorFaceIntegrator(BilinearFormIntegrator * bfi_real,
-                                  BilinearFormIntegrator * bfi_imag);
-
-   /// Adds new boundary Face Integrator. Assumes ownership of @a bfi.
-   void AddBdrFaceIntegrator(BilinearFormIntegrator * bfi_real,
-                             BilinearFormIntegrator * bfi_imag);
-
-   /** @brief Adds new boundary Face Integrator, restricted to specific boundary
-       attributes.
-
-       Assumes ownership of @a bfi.
-
-       The mfem::array @a bdr_marker is stored internally as a pointer to the given
-       mfem::Array<int> object. */
-   void AddBdrFaceIntegrator(BilinearFormIntegrator * bfi_real,
-                             BilinearFormIntegrator * bfi_imag,
-                             Array<int> & bdr_marker);
-
-   /** @brief Add a trace face integrator. Assumes ownership of @a bfi.
-
-      This type of integrator assembles terms over all faces of the mesh using
-      the face FE from the trial space and the two adjacent volume FEs from
-      the test space. */
-   void AddTraceFaceIntegrator(BilinearFormIntegrator * bfi_real,
-                               BilinearFormIntegrator * bfi_imag);
-
-   /// Adds a boundary trace face integrator. Assumes ownership of @a bfi.
-   void AddBdrTraceFaceIntegrator(BilinearFormIntegrator * bfi_real,
-                                  BilinearFormIntegrator * bfi_imag);
-
-   /// Adds a boundary trace face integrator. Assumes ownership of @a bfi.
-   void AddBdrTraceFaceIntegrator(BilinearFormIntegrator * bfi_real,
-                                  BilinearFormIntegrator * bfi_imag,
-                                  Array<int> &bdr_marker);
-
-   /// Assemble the local matrix
-   void Assemble(int skip_zeros = 1);
-
-   /// Finalizes the matrix initialization.
-   void Finalize(int skip_zeros = 1);
-
-   /// Updates the internal mixed forms with the new finite element space.
-   virtual void Update();
-
-   /// Returns the matrix assembled on the true dofs, i.e. P^t A P.
-   /** The returned matrix has to be deleted by the caller. */
-   ComplexHypreParMatrix * ParallelAssemble();
-
-   void FormRectangularLinearSystem(const Array<int> & ess_trial_tdof_list,
-                                    const Array<int> & ess_test_tdof_list,
-                                    Vector & x,
-                                    Vector & b,
-                                    OperatorHandle & A,
-                                    Vector & X,
-                                    Vector & B);
-
-   void FormRectangularSystemMatrix(const Array<int> & ess_trial_tdof_list,
-                                    const Array<int> & ess_test_tdof_list,
-                                    OperatorHandle & A);
-
-   virtual ~ParMixedSesquilinearForm();
-};
-
 #endif // MFEM_USE_MPI

 }
@@ -114,10 +114,6 @@ void ConduitDataCollection::Save()
                                         n_mesh["fields"][name]);
   }

-   // TODO: in parallel, we need to call ParFiniteElementSpace::ApplyDofSigns
-   // for all ParGridFunction objects before and after saving, see
-   // ParGridFunction::Save.
-
   // save mesh data
   SaveMeshAndFields(myid,
                     n_mesh,
@@ -1457,8 +1453,6 @@ ConduitDataCollection::LoadMeshAndFields(int domain_id,
 std::string
 ConduitDataCollection::ElementTypeToShapeName(Element::Type element_type)
 {
-   // Adapted from SidreDataCollection
-
   // Note -- the mapping from Element::Type to string is based on
   //   enum Element::Type { POINT, SEGMENT, TRIANGLE, QUADRILATERAL,
   //                        TETRAHEDRON, HEXAHEDRON };
@@ -34,10 +34,10 @@ namespace mfem
    - HDF5 library, https://support.hdfgroup.org/HDF5

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

    This class also provides public static methods that convert between MFEM
    Meshes and GridFunctions and Conduit Mesh Blueprint descriptions.
@@ -492,8 +492,6 @@ void VisItDataCollection::SaveRootFile()
                           to_padded_string(cycle, pad_digits_cycle) +
                           ".mfem_root";
   std::ofstream root_file(root_name);
-   MFEM_VERIFY(root_file.is_open(),
-               "Failed to open ofstream " << root_name);
   root_file << GetVisItRootString();
   if (!root_file)
   {
@@ -809,7 +807,7 @@ ParaViewDataCollectionBase::ParaViewDataCollectionBase(

 void ParaViewDataCollectionBase::SetLevelsOfDetail(int levels_of_detail_)
 {
-   levels_of_detail = std::max(levels_of_detail_, 1);
+   levels_of_detail = levels_of_detail_;
 }

 void ParaViewDataCollectionBase::SetHighOrderOutput(bool high_order_output_)
@@ -979,10 +977,7 @@ void ParaViewDataCollection::Save()
   // Save the local part of the mesh and grid functions fields to the local
   // VTU file. Also save coefficient fields.
   {
-      std::string os_str = vtu_prefix + GenerateVTUFileName("proc", myid);
-      std::ofstream os(os_str);
-      MFEM_VERIFY(os.is_open(),
-                  "Failed to open ofstream " << os_str);
+      std::ofstream os(vtu_prefix + GenerateVTUFileName("proc", myid));
      os.precision(precision);
      SaveDataVTU(os, levels_of_detail);
   }
@@ -994,10 +989,7 @@ void ParaViewDataCollection::Save()
                  "QuadratureFunction output is not supported for "
                  "ParaViewDataCollection on domain boundary!");
      const std::string &field_name = qfield.first;
-      std::string os_str = vtu_prefix + GenerateVTUFileName(field_name, myid);
-      std::ofstream os(os_str);
-      MFEM_VERIFY(os.is_open(),
-                  "Failed to open ofstream " << os_str);
+      std::ofstream os(vtu_prefix + GenerateVTUFileName(field_name, myid));
      qfield.second->SaveVTU(os, pv_data_format, GetCompressionLevel(), field_name);
   }

@@ -1008,10 +1000,7 @@ void ParaViewDataCollection::Save()
   {
      // Create the main PVTU file
      {
-         std::string os_str = vtu_prefix + GeneratePVTUFileName("data");
-         std::ofstream pvtu_out(os_str);
-         MFEM_VERIFY(pvtu_out.is_open(),
-                     "Failed to open ofstream " << os_str);
+         std::ofstream pvtu_out(vtu_prefix + GeneratePVTUFileName("data"));
         WritePVTUHeader(pvtu_out);

         // Grid function fields and coefficient fields
@@ -1066,10 +1055,8 @@ void ParaViewDataCollection::Save()
         const std::string &q_field_name = q_field.first;
         std::string q_fname = GeneratePVTUPath() + "/"
                               + GeneratePVTUFileName(q_field_name);
-         std::string os_str = col_path + "/" + q_fname;
-         std::ofstream pvtu_out(os_str);
-         MFEM_VERIFY(pvtu_out.is_open(),
-                     "Failed to open ofstream " << os_str);
+
+         std::ofstream pvtu_out(col_path + "/" + q_fname);
         WritePVTUHeader(pvtu_out);
         int vec_dim = q_field.second->GetVDim();
         pvtu_out << "<PPointData>\n";
@@ -1181,14 +1168,12 @@ void ParaViewDataCollection::SaveGFieldVTU(std::ostream &os, int ref_,
   DenseMatrix vval, pmat;
   std::vector<char> buf;
   int vec_dim = it->second->VectorDim();
-   int map_type = it->second->FESpace()->GetTypicalFE()->GetMapType();
   os << "<DataArray type=\"" << GetDataTypeString()
      << "\" Name=\"" << it->first
      << "\" NumberOfComponents=\"" << vec_dim << "\" "
      << VTKComponentLabels(vec_dim) << " "
      << "format=\"" << GetDataFormatString() << "\" >" << '\n';
-   if (vec_dim == 1 && (map_type == FiniteElement::VALUE ||
-                        map_type == FiniteElement::INTEGRAL))
+   if (vec_dim == 1)
   {
      for (int i = 0; i < mesh->GetNE(); i++)
      {
@@ -57,7 +57,7 @@ void DGMassApply(const int e,
      }
      else if (DIM == 3)
      {
-         SmemPAMassApply3D_Element<TD1D,TQ1D,NBZ,ACCUM>(e, NE, B, pa_data, x, y);
+         SmemPAMassApply3D_Element<TD1D,TQ1D,ACCUM>(e, NE, B, pa_data, x, y);
      }
      else
      {
@@ -320,8 +320,8 @@ public:
    error estimation procedure where the flux averaging is replaced by a global
    L2 projection (requiring a mass matrix solve).

-    The required BilinearFormIntegrator must implement the method
-    ComputeElementFlux().
+    The required BilinearFormIntegrator must implement the methods
+    ComputeElementFlux() and ComputeFluxEnergy().

    Implemented for the parallel case only.
 */
@@ -357,8 +357,8 @@ protected:

 public:
   /** @brief Construct a new L2ZienkiewiczZhuEstimator object.
-       @param integ    This BilinearFormIntegrator must implement the method
-                       ComputeElementFlux().
+       @param integ    This BilinearFormIntegrator must implement the methods
+                       ComputeElementFlux() and ComputeFluxEnergy().
       @param sol      The solution field whose error is to be estimated.
       @param flux_fes The L2ZienkiewiczZhuEstimator assumes ownership of this
                       FiniteElementSpace and will call its Update() method when
@@ -382,8 +382,8 @@ public:
   { }

   /** @brief Construct a new L2ZienkiewiczZhuEstimator object.
-       @param integ    This BilinearFormIntegrator must implement the method
-                       ComputeElementFlux().
+       @param integ    This BilinearFormIntegrator must implement the methods
+                       ComputeElementFlux() and ComputeFluxEnergy().
       @param sol      The solution field whose error is to be estimated.
       @param flux_fes The L2ZienkiewiczZhuEstimator does NOT assume ownership
                       of this FiniteElementSpace; will call its Update() method
@@ -167,15 +167,7 @@ public:
      /** @brief Full multidimensional representation which does not use tensor
          product structure. The ordering of the degrees of freedom is the
          same as TENSOR, but the sizes of B and G are the same as FULL.*/
-      LEXICOGRAPHIC_FULL,
-
-      /** @brief Ragged tensor product representation using 1D matrices/tensors
-          with dimensions using 1D number of quadrature points and ragged tensor degrees of
-          freedom. */
-      /** Used only for partial assembly of the H1 positive basis. The
-          size of B is d1d x qnpt x dim. Since different Gauss-Jacobi quadrature rules
-          are employed in each dimension, we need to store dim arrays. */
-      RAGGED_TENSOR
+      LEXICOGRAPHIC_FULL
   };

   /// Describes the contents of the #B, #Bt, #G, and #Gt arrays, see #Mode.
@@ -236,39 +228,6 @@ public:
      const Array<DofToQuad*> &dof2quad_array,
      const IntegrationRule &ir,
      DofToQuad::Mode mode);
-
-   virtual ~DofToQuad() = default;
-};
-
-/** @brief Structure representing the matrices/tensors needed to evaluate (in
-    reference space) the values, gradients, divergences, or curls of a positive
-    FiniteElement on simplices at the quadrature points of Stroud conical quadrature. */
-class RaggedDofToQuad : public DofToQuad
-{
-public:
-   /** @brief Special basis function structures for positive (Bernstein) basis with
-      partial assembly. The storage layout of Ba1 is ndof x nqpt for scalar elements.
-      The storage layout of Ba2 is ndof x ndof x nqpt. In particular, we have
-            Ba2(iqpt, a1, a2) = B^{p-a1}_{a2}(x_{iqpt}). */
-   Array<real_t> Ba1, Ba2, Ba3;
-   Array<real_t> Ba1t, Ba2t, Ba3t;
-
-   /** @brief Special structures for gradients of positive basis with partial assembly.
-      The gradient arrays exploit properties of the Bernstein basis which allow grad(B^p_alpha)
-      to be expressed as the sum of products of B^{p-1}_alpha and the barycentric coordinates.
-      Thus, Ga1 and Ga2 simply contain the ragged tensor product components of B^{p-1}_alpha */
-   Array<real_t> Ga1, Ga2, Ga3;
-   Array<real_t> Ga1t, Ga2t, Ga3t;
-
-   /** @brief Mapping from the Bernstein multi-index (a_1, ..., a_d) to the lexicographic
-       dof index. */
-   Array<int> lex_map;
-
-   Array<int> forward_map2d_diff, forward_map3d_diff;
-   Array<int> inverse_map2d_diff, inverse_map3d_diff;
-
-   Array<int> forward_map2d_mass, forward_map3d_mass;
-   Array<int> inverse_map2d_mass, inverse_map3d_mass;
 };

 /// Describes the function space on each element
@@ -307,12 +307,12 @@ public:

   /** @brief virtual function which evaluates the values of all
       shape functions at a given point ip and stores
-       them in the vector shape of dimension Dof (6) */
+       them in the vector shape of dimension Dof (4) */
   void CalcShape(const IntegrationPoint &ip, Vector &shape) const override;

   /** @brief virtual function which evaluates the values of all
       partial derivatives of all shape functions at a given
-       point ip and stores them in the matrix dshape (Dof x Dim) (6 x 3)
+       point ip and stores them in the matrix dshape (Dof x Dim) (4 x 3)
       so that each row contains the derivatives of one shape function */
   void CalcDShape(const IntegrationPoint &ip,
                   DenseMatrix &dshape) const override;
@@ -336,12 +336,12 @@ public:

   /** @brief virtual function which evaluates the values of all
       shape functions at a given point ip and stores
-       them in the vector shape of dimension Dof (5) */
+       them in the vector shape of dimension Dof (4) */
   void CalcShape(const IntegrationPoint &ip, Vector &shape) const override;

   /** @brief virtual function which evaluates the values of all
       partial derivatives of all shape functions at a given
-       point ip and stores them in the matrix dshape (Dof x Dim) (5 x 3)
+       point ip and stores them in the matrix dshape (Dof x Dim) (4 x 3)
       so that each row contains the derivatives of one shape function */
   void CalcDShape(const IntegrationPoint &ip,
                   DenseMatrix &dshape) const override;
@@ -1757,45 +1757,22 @@ H1_BergotPyramidElement::H1_BergotPyramidElement(const int p, const int btype)
      real_t y = (ip.z < 1.0) ? (ip.y / (1.0 - ip.z)) : 0.0;
      real_t z = ip.z;

-      if (std::abs(z - 1.0) < apex_tol)
-      {
-         // Compute the limit of the basis functions as z->1 with x and y on the
-         // line between the center of the base and the apex
-         o = 0;
-         for (int i = 0; i <= p; i++)
-            for (int j = 0; j <= p; j++)
-            {
-               int maxij = std::max(i, j);
-               for (int k = 0; k <= p - maxij; k++)
-                  if (i == 0 && j == 0)
-                  {
-                     T(o++, m) = ((k + 3.) * k + 2.) / 2.;
-                  }
-                  else
-                  {
-                     T(o++, m) = 0.;
-                  }
-            }
-      }
-      else
-      {
-         poly1d.CalcLegendre(p, x, shape_x.GetData());
-         poly1d.CalcLegendre(p, y, shape_y.GetData());
+      poly1d.CalcLegendre(p, x, shape_x.GetData());
+      poly1d.CalcLegendre(p, y, shape_y.GetData());

-         o = 0;
-         for (int i = 0; i <= p; i++)
+      o = 0;
+      for (int i = 0; i <= p; i++)
+      {
+         for (int j = 0; j <= p; j++)
         {
-            for (int j = 0; j <= p; j++)
-            {
-               int maxij = std::max(i, j);
-               FuentesPyramid::CalcScaledJacobi(p-maxij, 2.0 * (maxij + 1.0),
-                                                z, 1.0, shape_z);
+            int maxij = std::max(i, j);
+            FuentesPyramid::CalcScaledJacobi(p-maxij, 2.0 * (maxij + 1.0),
+                                             z, 1.0, shape_z);

-               for (int k = 0; k <= p - maxij; k++)
-               {
-                  T(o++, m) = shape_x(i) * shape_y(j) * shape_z(k) *
-                              pow(1.0 - ip.z, maxij);
-               }
+            for (int k = 0; k <= p - maxij; k++)
+            {
+               T(o++, m) = shape_x(i) * shape_y(j) * shape_z(k) *
+                           pow(1.0 - ip.z, maxij);
            }
         }
      }
@@ -1816,44 +1793,25 @@ void H1_BergotPyramidElement::CalcShape(const IntegrationPoint &ip,
   Vector u(dof);
 #endif

-   const real_t x = (ip.z < 1.0) ? (ip.x / (1.0 - ip.z)) : 0.0;
-   const real_t y = (ip.z < 1.0) ? (ip.y / (1.0 - ip.z)) : 0.0;
-   const real_t z = ip.z;
+   real_t x = (ip.z < 1.0) ? (ip.x / (1.0 - ip.z)) : 0.0;
+   real_t y = (ip.z < 1.0) ? (ip.y / (1.0 - ip.z)) : 0.0;
+   real_t z = ip.z;

-   if (std::abs(z - 1.0) < apex_tol)
-   {
-      // Compute the limit of the basis functions as z->1 with x and y on the
-      // line between the center of the base and the apex
-      u = 0.;
-      int o = 0;
-      for (int i = 0; i <= p; i++)
-         for (int j = 0; j <= p; j++)
-         {
-            int maxij = std::max(i, j);
-            for (int k = 0; k <= p - maxij; k++, o++)
-               if (i == 0 && j == 0)
-               {
-                  u(o) = ((k + 3.) * k + 2.) / 2.;
-               }
-         }
-   }
-   else
-   {
-      poly1d.CalcLegendre(p, x, shape_x.GetData());
-      poly1d.CalcLegendre(p, y, shape_y.GetData());
+   poly1d.CalcLegendre(p, x, shape_x.GetData());
+   poly1d.CalcLegendre(p, y, shape_y.GetData());
+
+   int o = 0;
+   for (int i = 0; i <= p; i++)
+      for (int j = 0; j <= p; j++)
+      {
+         int maxij = std::max(i, j);
+         FuentesPyramid::CalcScaledJacobi(p-maxij, 2.0 * (maxij + 1.0), z, 1.0,
+                                          shape_z);
+         for (int k = 0; k <= p - maxij; k++)
+            u[o++] = shape_x(i) * shape_y(j) * shape_z(k) *
+                     pow(1.0 - ip.z, maxij);
+      }

-      int o = 0;
-      for (int i = 0; i <= p; i++)
-         for (int j = 0; j <= p; j++)
-         {
-            int maxij = std::max(i, j);
-            FuentesPyramid::CalcScaledJacobi(p-maxij, 2.0 * (maxij + 1.0),
-                                             z, 1.0, shape_z);
-            for (int k = 0; k <= p - maxij; k++)
-               u[o++] = shape_x(i) * shape_y(j) * shape_z(k) *
-                        pow(1.0 - ip.z, maxij);
-         }
-   }
   Ti.Mult(u, shape);
 }

@@ -1872,68 +1830,37 @@ void H1_BergotPyramidElement::CalcDShape(const IntegrationPoint &ip,
   Vector dshape_z(order+1);
   Vector dshape_z_dt(order+1);
 #endif
-   const real_t x = (ip.z < 1.0) ? (ip.x / (1.0 - ip.z)) : 0.0;
-   const real_t y = (ip.z < 1.0) ? (ip.y / (1.0 - ip.z)) : 0.0;
-   const real_t z = ip.z;
+   real_t x = (ip.z < 1.0) ? (ip.x / (1.0 - ip.z)) : 0.0;
+   real_t y = (ip.z < 1.0) ? (ip.y / (1.0 - ip.z)) : 0.0;
+   real_t z = ip.z;

-   if (std::abs(z - 1.0) < apex_tol)
-   {
-      // Compute the limit of the gradients of the basis functions as
-      // z->1 with x and y on the line between the center of the base and the
-      // apex
-      du = 0.;
-      int o = 0;
-      for (int i = 0; i <= p; i++)
-         for (int j = 0; j <= p; j++)
+   poly1d.CalcLegendre(p, x, shape_x.GetData(), dshape_x.GetData());
+   poly1d.CalcLegendre(p, y, shape_y.GetData(), dshape_y.GetData());
+
+   int o = 0;
+   for (int i = 0; i <= p; i++)
+      for (int j = 0; j <= p; j++)
+      {
+         int maxij = std::max(i, j);
+         FuentesPyramid::CalcScaledJacobi(p-maxij, 2.0 * (maxij + 1.0), z, 1.0,
+                                          shape_z, dshape_z, dshape_z_dt);
+
+         for (int k = 0; k <= p - maxij; k++, o++)
         {
-            int maxij = std::max(i, j);
-            for (int k = 0; k <= p - maxij; k++, o++)
-            {
-               if (i == 0 && j == 0)
-               {
-                  du(o,2) = (((k + 6.) * k + 11.) * k + 6.) * k / 6.;
-               }
-               else if (i == 1 && j == 0)
-               {
-                  du(o,0) = ((((k + 10.) * k + 35.) * k + 50.) * k + 24.) / 24.;
-               }
-               else if (i == 0 && j == 1)
-               {
-                  du(o,1) = ((((k + 10.) * k + 35.) * k + 50.) * k + 24.) / 24.;
-               }
-            }
+            du(o,0) = dshape_x(i) * shape_y(j) * shape_z(k) *
+                      pow(1.0 - ip.z, maxij - 1);
+            du(o,1) = shape_x(i) * dshape_y(j) * shape_z(k) *
+                      pow(1.0 - ip.z, maxij - 1);
+            du(o,2) = shape_x(i) * shape_y(j) * dshape_z(k) *
+                      pow(1.0 - ip.z, maxij) +
+                      (ip.x * dshape_x(i) * shape_y(j) +
+                       ip.y * shape_x(i) * dshape_y(j)) *
+                      shape_z(k) * pow(1.0 - ip.z, maxij - 2) -
+                      maxij * shape_x(i) * shape_y(j) * shape_z(k) *
+                      pow(1.0 - ip.z, maxij - 1);
         }
-   }
-   else
-   {
-      poly1d.CalcLegendre(p, x, shape_x.GetData(), dshape_x.GetData());
-      poly1d.CalcLegendre(p, y, shape_y.GetData(), dshape_y.GetData());
+      }

-      int o = 0;
-      for (int i = 0; i <= p; i++)
-         for (int j = 0; j <= p; j++)
-         {
-            int maxij = std::max(i, j);
-            FuentesPyramid::CalcScaledJacobi(p-maxij, 2.0 * (maxij + 1.0),
-                                             z, 1.0,
-                                             shape_z, dshape_z, dshape_z_dt);
-
-            for (int k = 0; k <= p - maxij; k++, o++)
-            {
-               du(o,0) = dshape_x(i) * shape_y(j) * shape_z(k) *
-                         pow(1.0 - ip.z, maxij - 1);
-               du(o,1) = shape_x(i) * dshape_y(j) * shape_z(k) *
-                         pow(1.0 - ip.z, maxij - 1);
-               du(o,2) = shape_x(i) * shape_y(j) * dshape_z(k) *
-                         pow(1.0 - ip.z, maxij) +
-                         (ip.x * dshape_x(i) * shape_y(j) +
-                          ip.y * shape_x(i) * dshape_y(j)) *
-                         shape_z(k) * pow(1.0 - ip.z, maxij - 2) -
-                         maxij * shape_x(i) * shape_y(j) * shape_z(k) *
-                         (maxij > 0 ? pow(1.0 - ip.z, maxij - 1) : 0.0);
-            }
-         }
-   }
   Ti.Mult(du, dshape);
 }

@@ -208,8 +208,6 @@ private:
 #endif
   DenseMatrixInverse Ti;

-   static constexpr real_t apex_tol = 1e-8;
-
 public:
   H1_BergotPyramidElement(const int p,
                           const int btype = BasisType::GaussLobatto);
@@ -1106,16 +1106,9 @@ L2_BergotPyramidElement::L2_BergotPyramidElement(const int p, const int btype)
         {
            const real_t wik = op[i] + op[k] + op[p-i-k];
            const real_t w = wik * wjk * op[p-k];
-            if (std::abs(w) < apex_tol)
-            {
-               Nodes.IntPoint(o++).Set3(0.,0.,1.);
-            }
-            else
-            {
-               Nodes.IntPoint(o++).Set3(op[i] * (op[j] + op[p-j-k]) / w,
-                                        op[j] * (op[i] + op[p-i-k]) / w,
-                                        op[k] * op[p-k] / w);
-            }
+            Nodes.IntPoint(o++).Set3(op[i] * (op[j] + op[p-j-k]) / w,
+                                     op[j] * (op[j] + op[p-j-k]) / w,
+                                     op[k] * op[p-k] / w);
         }
      }

@@ -1132,45 +1125,22 @@ L2_BergotPyramidElement::L2_BergotPyramidElement(const int p, const int btype)
      const real_t y = (ip.z < 1.0) ? (ip.y / (1.0 - ip.z)) : 0.0;
      const real_t z = ip.z;

-      if (std::abs(z - 1.0) < apex_tol)
-      {
-         // Compute the limit of the basis functions as z->1 with x and y on the
-         // line between the center of the base and the apex
-         o = 0;
-         for (int i = 0; i <= p; i++)
-            for (int j = 0; j <= p; j++)
-            {
-               int maxij = std::max(i, j);
-               for (int k = 0; k <= p - maxij; k++)
-                  if (i == 0 && j == 0)
-                  {
-                     T(o++, m) = ((k + 3.) * k + 2.) / 2.;
-                  }
-                  else
-                  {
-                     T(o++, m) = 0.;
-                  }
-            }
-      }
-      else
-      {
-         poly1d.CalcLegendre(p, x, shape_x.GetData());
-         poly1d.CalcLegendre(p, y, shape_y.GetData());
+      poly1d.CalcLegendre(p, x, shape_x.GetData());
+      poly1d.CalcLegendre(p, y, shape_y.GetData());

-         o = 0;
-         for (int i = 0; i <= p; i++)
+      o = 0;
+      for (int i = 0; i <= p; i++)
+      {
+         for (int j = 0; j <= p; j++)
         {
-            for (int j = 0; j <= p; j++)
-            {
-               int maxij = std::max(i, j);
-               FuentesPyramid::CalcScaledJacobi(p-maxij, 2.0 * (maxij + 1.0),
-                                                z, 1.0, shape_z);
+            int maxij = std::max(i, j);
+            FuentesPyramid::CalcScaledJacobi(p-maxij, 2.0 * (maxij + 1.0),
+                                             z, 1.0, shape_z);

-               for (int k = 0; k <= p - maxij; k++)
-               {
-                  T(o++, m) = shape_x(i) * shape_y(j) * shape_z(k) *
-                              pow(1.0 - ip.z, maxij);
-               }
+            for (int k = 0; k <= p - maxij; k++)
+            {
+               T(o++, m) = shape_x(i) * shape_y(j) * shape_z(k) *
+                           pow(1.0 - ip.z, maxij);
            }
         }
      }
@@ -1195,41 +1165,26 @@ void L2_BergotPyramidElement::CalcShape(const IntegrationPoint &ip,
   const real_t y = (ip.z < 1.0) ? (ip.y / (1.0 - ip.z)) : 0.0;
   const real_t z = ip.z;

-   if (std::abs(z - 1.0) < apex_tol)
-   {
-      // Compute the limit of the basis functions as z->1 with x and y on the
-      // line between the center of the base and the apex
-      u = 0.;
-      int o = 0;
-      for (int i = 0; i <= p; i++)
-         for (int j = 0; j <= p; j++)
-         {
-            int maxij = std::max(i, j);
-            for (int k = 0; k <= p - maxij; k++, o++)
-               if (i == 0 && j == 0)
-               {
-                  u(o) = ((k + 3.) * k + 2.) / 2.;
-               }
-         }
-   }
-   else
-   {
-      poly1d.CalcLegendre(p, x, shape_x.GetData());
-      poly1d.CalcLegendre(p, y, shape_y.GetData());
+   poly1d.CalcLegendre(p, x, shape_x.GetData());
+   poly1d.CalcLegendre(p, y, shape_y.GetData());

-      int o = 0;
-      for (int i = 0; i <= p; i++)
-         for (int j = 0; j <= p; j++)
-         {
-            int maxij = std::max(i, j);
-            FuentesPyramid::CalcScaledJacobi(p-maxij, 2.0 * (maxij + 1.0),
-                                             z, 1.0, shape_z);
+   int o = 0;
+   for (int i = 0; i <= p; i++)
+   {
+      for (int j = 0; j <= p; j++)
+      {
+         int maxij = std::max(i, j);
+         FuentesPyramid::CalcScaledJacobi(p-maxij, 2.0 * (maxij + 1.0), z, 1.0,
+                                          shape_z);

-            for (int k = 0; k <= p - maxij; k++)
-               u[o++] = shape_x(i) * shape_y(j) * shape_z(k) *
-                        pow(1.0 - ip.z, maxij);
+         for (int k = 0; k <= p - maxij; k++)
+         {
+            u[o++] = shape_x(i) * shape_y(j) * shape_z(k) *
+                     pow(1.0 - ip.z, maxij);
         }
+      }
   }
+
   Ti.Mult(u, shape);
 }

@@ -1253,64 +1208,35 @@ void L2_BergotPyramidElement::CalcDShape(const IntegrationPoint &ip,
   const real_t y = (ip.z < 1.0) ? (ip.y / (1.0 - ip.z)) : 0.0;
   const real_t z = ip.z;

-   if (std::abs(z - 1.0) < apex_tol)
-   {
-      // Compute the limit of the gradients of the basis functions as
-      // z->1 with x and y on the line between the center of the base and the
-      // apex
-      du = 0.;
-      int o = 0;
-      for (int i = 0; i <= p; i++)
-         for (int j = 0; j <= p; j++)
-         {
-            int maxij = std::max(i, j);
-            for (int k = 0; k <= p - maxij; k++, o++)
-            {
-               if (i == 0 && j == 0)
-               {
-                  du(o,2) = (((k + 6.) * k + 11.) * k + 6.) * k / 6.;
-               }
-               else if (i == 1 && j == 0)
-               {
-                  du(o,0) = ((((k + 10.) * k + 35.) * k + 50.) * k + 24.) / 24.;
-               }
-               else if (i == 0 && j == 1)
-               {
-                  du(o,1) = ((((k + 10.) * k + 35.) * k + 50.) * k + 24.) / 24.;
-               }
-            }
-         }
-   }
-   else
-   {
-      Poly_1D::CalcLegendre(p, x, shape_x.GetData(), dshape_x.GetData());
-      Poly_1D::CalcLegendre(p, y, shape_y.GetData(), dshape_y.GetData());
+   Poly_1D::CalcLegendre(p, x, shape_x.GetData(), dshape_x.GetData());
+   Poly_1D::CalcLegendre(p, y, shape_y.GetData(), dshape_y.GetData());

-      int o = 0;
-      for (int i = 0; i <= p; i++)
-         for (int j = 0; j <= p; j++)
-         {
-            int maxij = std::max(i, j);
-            FuentesPyramid::CalcScaledJacobi(p-maxij, 2.0 * (maxij + 1.0),
-                                             z, 1.0,
-                                             shape_z, dshape_z, dshape_z_dt);
+   int o = 0;
+   for (int i = 0; i <= p; i++)
+   {
+      for (int j = 0; j <= p; j++)
+      {
+         int maxij = std::max(i, j);
+         FuentesPyramid::CalcScaledJacobi(p-maxij, 2.0 * (maxij + 1.0), z, 1.0,
+                                          shape_z, dshape_z, dshape_z_dt);

-            for (int k = 0; k <= p - maxij; k++, o++)
-            {
-               du(o,0) = dshape_x(i) * shape_y(j) * shape_z(k) *
-                         pow(1.0 - ip.z, maxij - 1);
-               du(o,1) = shape_x(i) * dshape_y(j) * shape_z(k) *
-                         pow(1.0 - ip.z, maxij - 1);
-               du(o,2) = shape_x(i) * shape_y(j) * dshape_z(k) *
-                         pow(1.0 - ip.z, maxij) +
-                         (ip.x * dshape_x(i) * shape_y(j) +
-                          ip.y * shape_x(i) * dshape_y(j)) *
-                         shape_z(k) * pow(1.0 - ip.z, maxij - 2) -
-                         maxij * shape_x(i) * shape_y(j) * shape_z(k) *
-                         (maxij > 0 ? pow(1.0 - ip.z, maxij - 1) : 0.0);
-            }
+         for (int k = 0; k <= p - maxij; k++, o++)
+         {
+            du(o,0) = dshape_x(i) * shape_y(j) * shape_z(k) *
+                      pow(1.0 - ip.z, maxij - 1);
+            du(o,1) = shape_x(i) * dshape_y(j) * shape_z(k) *
+                      pow(1.0 - ip.z, maxij - 1);
+            du(o,2) = shape_x(i) * shape_y(j) * dshape_z(k) *
+                      pow(1.0 - ip.z, maxij) +
+                      (ip.x * dshape_x(i) * shape_y(j) +
+                       ip.y * shape_x(i) * dshape_y(j)) *
+                      shape_z(k) * pow(1.0 - ip.z, maxij - 2) -
+                      ((maxij > 0) ? (maxij * shape_x(i) * shape_y(j) * shape_z(k) *
+                                      pow(1.0 - ip.z, maxij - 1)) : 0.0);
         }
+      }
   }
+
   Ti.Mult(du, dshape);
 }

@@ -225,8 +225,6 @@ private:
 #endif
   DenseMatrixInverse Ti;

-   static constexpr real_t apex_tol = 1e-8;
-
 public:
   /// Construct the L2_PyramidElement of order @a p and BasisType @a btype
   L2_BergotPyramidElement(const int p,
@@ -1282,49 +1282,12 @@ ND_SegmentElement::ND_SegmentElement(const int p, const int ob_type)
   }
 }

-void ND_SegmentElement::CalcShape(const IntegrationPoint &ip,
-                                  Vector &shape) const
-{
-   if (obasis1d.IsIntegratedType()) { obasis1d.ScaleIntegrated(false); }
-   obasis1d.Eval(ip.x, shape);
-}
-
 void ND_SegmentElement::CalcVShape(const IntegrationPoint &ip,
                                   DenseMatrix &shape) const
 {
   Vector vshape(shape.Data(), dof);

-   CalcShape(ip, vshape);
-}
-
-void ND_SegmentElement::ProjectIntegrated(VectorCoefficient &vc,
-                                          ElementTransformation &Trans,
-                                          Vector &dofs) const
-{
-   MFEM_ASSERT(obasis1d.IsIntegratedType(), "Not integrated type");
-   real_t vk[Geometry::MaxDim];
-   Vector xk(vk, vc.GetVDim());
-
-   const real_t *cp = poly1d.ClosedPoints(dof, BasisType::GaussLobatto);
-   const IntegrationRule &ir = IntRules.Get(Geometry::SEGMENT, dof);
-   IntegrationPoint ip;
-
-   for (int i = 0; i < dof; i++)
-   {
-      const real_t h = cp[i+1] - cp[i];
-      real_t val = 0.0;
-
-      for (int q = 0; q < ir.GetNPoints(); q++)
-      {
-         const IntegrationPoint &ip1d = ir.IntPoint(q);
-         ip.x = cp[i] + h*ip1d.x;
-         Trans.SetIntPoint(&ip);
-         vc.Eval(xk, Trans, ip);
-         val += ip1d.weight*Trans.Jacobian().InnerProduct(tk, vk);
-      }
-
-      dofs(i) = val*h;
-   }
+   obasis1d.Eval(ip.x, vshape);
 }

 const real_t ND_WedgeElement::tk[15] =
@@ -303,7 +303,8 @@ public:
   /** @brief Construct the ND_SegmentElement of order @a p and open
       BasisType @a ob_type */
   ND_SegmentElement(const int p, const int ob_type = BasisType::GaussLegendre);
-   void CalcShape(const IntegrationPoint &ip, Vector &shape) const override;
+   void CalcShape(const IntegrationPoint &ip, Vector &shape) const override
+   { obasis1d.Eval(ip.x, shape); }
   void CalcVShape(const IntegrationPoint &ip,
                   DenseMatrix &shape) const override;
   void CalcVShape(ElementTransformation &Trans,
@@ -324,10 +325,7 @@ public:
   using FiniteElement::Project;
   void Project(VectorCoefficient &vc,
                ElementTransformation &Trans, Vector &dofs) const override
-   {
-      if (obasis1d.IsIntegratedType()) { ProjectIntegrated(vc, Trans, dofs); }
-      else { Project_ND(tk, dof2tk, vc, Trans, dofs); }
-   }
+   { Project_ND(tk, dof2tk, vc, Trans, dofs); }
   void ProjectMatrixCoefficient(MatrixCoefficient &mc,
                                 ElementTransformation &T,
                                 Vector &dofs) const override
@@ -340,11 +338,6 @@ public:
                    ElementTransformation &Trans,
                    DenseMatrix &grad) const override
   { ProjectGrad_ND(tk, dof2tk, fe, Trans, grad); }
-
-protected:
-   void ProjectIntegrated(VectorCoefficient &vc,
-                          ElementTransformation &Trans,
-                          Vector &dofs) const;
 };

 class ND_WedgeElement : public VectorFiniteElement
@@ -557,101 +557,6 @@ H1Pos_TriangleElement::H1Pos_TriangleElement(const int p)
      }
 }

-const DofToQuad &H1Pos_TriangleElement::GetRaggedTensorDofToQuad(
-   const FiniteElement &fe, const IntegrationRule &ir,
-   DofToQuad::Mode mode,
-   Array<DofToQuad*> &dof2quad_array)
-{
-   DofToQuad *d2q = nullptr;
-   MFEM_VERIFY(mode == DofToQuad::RAGGED_TENSOR, "invalid mode requested");
-
-#if defined(MFEM_THREAD_SAFE) && defined(MFEM_USE_OPENMP)
-   #pragma omp critical (DofToQuad)
-#endif
-   {
-      for (int i = 0; i < dof2quad_array.Size(); i++)
-      {
-         d2q = dof2quad_array[i];
-         if (d2q->IntRule != &ir || d2q->mode != mode) { d2q = nullptr; }
-      }
-      if (!d2q)
-      {
-         d2q = new RaggedDofToQuad;
-         const int ndof = fe.GetOrder() + 1; // verify
-         const int nqpt = (int)floor(pow(ir.GetNPoints(), 1.0/fe.GetDim()) + 0.5);
-         d2q->FE = &fe;
-         d2q->IntRule = &ir;
-         d2q->mode = mode;
-         d2q->ndof = ndof;
-         d2q->nqpt = nqpt;
-
-         RaggedDofToQuad *rd2q = static_cast<RaggedDofToQuad*>(d2q);
-         rd2q->Ba1.SetSize(nqpt*ndof);
-         // second component of ragged tensor basis, technically dof*(dof-1)/2 entries
-         rd2q->Ba2.SetSize((int)nqpt*ndof*ndof);
-         rd2q->Ba1t.SetSize(nqpt*ndof);
-         rd2q->Ba2t.SetSize((int)nqpt*ndof*ndof);
-         // stores first component of ragged tensor basis with order p-1, for gradients only
-         rd2q->Ga1.SetSize(nqpt*(ndof -1));
-         // stores second component of ragged tensor basis with order p-1
-         rd2q->Ga2.SetSize(nqpt*(ndof-1)*(ndof -1));
-         rd2q->Ga1t.SetSize(nqpt*(ndof -1));
-         rd2q->Ga2t.SetSize(nqpt*(ndof-1)*(ndof -1));
-         rd2q->lex_map.SetSize(ndof * ndof);
-         Vector shape_a1(ndof), shape_a2(ndof * ndof);
-         Vector shape_Ga1(ndof-1), shape_Ga2((ndof-1) * (ndof-1));
-         for (int i = 0; i < nqpt; i++)
-         {
-            // The first 'nqpt' points in the first dimension 'ir' have the same x-coordinates as those
-            // of the 1D rule (ie. (2,0) Gauss-Jacobi rule). The first 'nqpt' points in the second dimension
-            // 'ir' have the same y-coordinates as those of the 1D rule for second dimension (i.e. (1,0)
-            // Gauss-Jacobi rule). Additionally, the Bernstein PA algorithms expect evaluation of the
-            // component 1D bases at the Stroud nodes pulled back to the unit square, so perform the pullback
-            // on the fly.
-            const real_t x = ir.IntPoint(i).x;
-            const real_t y = ir.IntPoint(nqpt*i).y / (1.0 - ir.IntPoint(nqpt*i).x);
-
-            Poly_1D::CalcBernstein(ndof-1, x, shape_a1);
-            Poly_1D::CalcBernstein(ndof-2, x, shape_Ga1);
-            for (int j = 0; j < ndof; j++)
-            {
-               rd2q->Ba1t[i+nqpt*j] = rd2q->Ba1[j+ndof*i] = shape_a1(j);
-               if (j < ndof-1)
-               {
-                  rd2q->Ga1t[i+nqpt*j] = rd2q->Ga1[j+(ndof-1)*i] = shape_Ga1(j);
-                  Poly_1D::CalcBernstein(ndof-2-j, y, shape_Ga2);
-               }
-
-               Poly_1D::CalcBernstein(ndof-1-j, y, shape_a2);
-               for (int k = 0; k < ndof-j; k++)
-               {
-                  rd2q->Ba2t[i + nqpt*(j + ndof*k)] = rd2q->Ba2[k + ndof*(j + ndof*i)] = shape_a2(
-                                                                                            k);
-                  if (j < ndof-1 && k < ndof-j-1)
-                  {
-                     rd2q->Ga2t[i + nqpt*(j + (ndof-1)*k)] = rd2q->Ga2[k + (ndof-1)*(j +
-                                                                                     (ndof-1)*i)] = shape_Ga2(k);
-                  }
-               }
-            }
-         }
-
-         // stores the mapping from 2D Bernstein multi-index (i,j,p-i-j) to the
-         // lexicographic DOF ordering
-         for (int i = 0; i < ndof; i++)
-         {
-            for (int j = 0; j < ndof-i; j++)
-            {
-               int idx = ((2 * (ndof-1) + 3) - j) * j / 2 + i;
-               rd2q->lex_map[j + ndof*i] = idx;
-            }
-         }
-         dof2quad_array.Append(d2q);
-      }
-   }
-   return *d2q;
-}
-
 // static method
 void H1Pos_TriangleElement::CalcShape(
   const int p, const real_t l1, const real_t l2, real_t *shape)
@@ -844,213 +749,6 @@ H1Pos_TetrahedronElement::H1Pos_TetrahedronElement(const int p)
         }
 }

-const DofToQuad &H1Pos_TetrahedronElement::GetRaggedTensorDofToQuad(
-   const FiniteElement &fe, const IntegrationRule &ir,
-   DofToQuad::Mode mode,
-   Array<DofToQuad*> &dof2quad_array)
-{
-   DofToQuad *d2q = nullptr;
-   MFEM_VERIFY(mode == DofToQuad::RAGGED_TENSOR, "invalid mode requested");
-
-#if defined(MFEM_THREAD_SAFE) && defined(MFEM_USE_OPENMP)
-   #pragma omp critical (DofToQuad)
-#endif
-   {
-      for (int i = 0; i < dof2quad_array.Size(); i++)
-      {
-         d2q = dof2quad_array[i];
-         if (d2q->IntRule != &ir || d2q->mode != mode) { d2q = nullptr; }
-      }
-      if (!d2q)
-      {
-         d2q = new RaggedDofToQuad;
-         const int ndof = fe.GetOrder() + 1; // verify
-         const int nqpt = (int)floor(pow(ir.GetNPoints(), 1.0/fe.GetDim()) + 0.5);
-         const int basis_dim2d = ndof*(ndof+1) / 2;
-         const int basis_dim3d = ndof*(ndof+1)*(ndof+2) / 6;
-         const int basis_dim2d_diff = (ndof-1)*(ndof) / 2;
-         const int basis_dim3d_diff = (ndof-1)*(ndof)*(ndof+1) / 6;
-         d2q->FE = &fe;
-         d2q->IntRule = &ir;
-         d2q->mode = mode;
-         d2q->ndof = ndof;
-         d2q->nqpt = nqpt;
-
-         RaggedDofToQuad *rd2q = static_cast<RaggedDofToQuad*>(d2q);
-         rd2q->Ba1.SetSize(nqpt * ndof);
-         // second component of ragged tensor basis, technically dof*(dof-1)/2 entries
-         rd2q->Ba2.SetSize(nqpt * basis_dim2d);
-         // third component of ragged tensor basis, technically dof*(dof-1)/2 entries
-         rd2q->Ba3.SetSize(nqpt * basis_dim3d);
-         rd2q->Ba1t.SetSize(nqpt * ndof);
-         rd2q->Ba2t.SetSize(nqpt * basis_dim2d);
-         rd2q->Ba3t.SetSize(nqpt * basis_dim3d);
-         // stores first component of ragged tensor basis with order p-1, for gradients only
-         rd2q->Ga1.SetSize(nqpt * (ndof-1));
-         // stores second component of ragged tensor basis with order p-1
-         rd2q->Ga2.SetSize(nqpt * basis_dim2d_diff);
-         // stores third component of ragged tensor basis with order p-1
-         rd2q->Ga3.SetSize(nqpt * basis_dim3d_diff);
-         rd2q->Ga1t.SetSize(nqpt * (ndof-1));
-         rd2q->Ga2t.SetSize(nqpt * basis_dim2d_diff);
-         rd2q->Ga3t.SetSize(nqpt * basis_dim3d_diff);
-         rd2q->lex_map.SetSize(ndof * ndof * ndof);
-
-         rd2q->forward_map2d_diff.SetSize((ndof-1) * (ndof-1));
-         rd2q->forward_map3d_diff.SetSize((ndof-1) * (ndof-1) * (ndof-1));
-         rd2q->inverse_map2d_diff.SetSize(2 * basis_dim2d_diff);
-         rd2q->inverse_map3d_diff.SetSize(3 * basis_dim3d_diff);
-
-         rd2q->forward_map2d_mass.SetSize(ndof * ndof);
-         rd2q->forward_map3d_mass.SetSize(ndof * ndof * ndof);
-         rd2q->inverse_map2d_mass.SetSize(2 * basis_dim2d);
-         rd2q->inverse_map3d_mass.SetSize(2 * basis_dim3d);
-
-         // forward and inverse maps for multi-index to collpased 1d index for diffusion, can combine
-         // these four loops, but need four idx's and clause for shorter diff loops
-         int idx = 0;
-         for (int i = 0; i < ndof-1; i++)
-         {
-            for (int j = 0; j < ndof-i-1; j++)
-            {
-               rd2q->forward_map2d_diff[j + (ndof-1)*i] = idx;
-               rd2q->inverse_map2d_diff[2*idx] = i;
-               rd2q->inverse_map2d_diff[1 + 2*idx] = j;
-               idx++;
-            }
-         }
-
-         idx = 0;
-         for (int k = 0; k < ndof-1; k++)
-         {
-            for (int j = 0; j < ndof-k-1; j++)
-            {
-               for (int i = 0; i < ndof-k-j-1; i++)
-               {
-                  rd2q->forward_map3d_diff[k + (ndof-1)*(j + (ndof-1)*i)] = idx;
-                  rd2q->inverse_map3d_diff[3*idx] = i;
-                  rd2q->inverse_map3d_diff[1 + 3*idx] = j;
-                  rd2q->inverse_map3d_diff[2 + 3*idx] = k;
-                  idx++;
-               }
-            }
-         }
-
-         // forward and inverse maps for multi-index to collpased 1d index for mass
-         idx = 0;
-         for (int j = 0; j < ndof; j++)
-         {
-            for (int i = 0; i < ndof-j; i++)
-            {
-               rd2q->forward_map2d_mass[j + ndof*i] = idx;
-               rd2q->inverse_map2d_mass[2*idx] = i;
-               rd2q->inverse_map2d_mass[1 + 2*idx] = j;
-               idx++;
-            }
-         }
-
-         idx = 0;
-         for (int k = 0; k < ndof; k++)
-         {
-            for (int j = 0; j < ndof-k; j++)
-            {
-               for (int i = 0; i < ndof-k-j; i++)
-               {
-                  rd2q->forward_map3d_mass[k + ndof*(j + ndof*i)] = idx;
-                  rd2q->inverse_map3d_mass[2*idx] = i;
-                  rd2q->inverse_map3d_mass[1 + 2*idx] = j;
-                  // d2q->inverse_map3d_mass[2 + 3*idx] = k;
-                  idx++;
-               }
-            }
-         }
-
-         Vector shape_a1(ndof), shape_a2(ndof * ndof), shape_a3(ndof * ndof * ndof);
-         Vector shape_Ga1(ndof-1), shape_Ga2(ndof-1), shape_Ga3(ndof-1);
-         for (int i = 0; i < nqpt; i++)
-         {
-            // The first 'nqpt' points in the first dimension 'ir' have the same x-coordinates as those
-            // of the 1D rule (ie. (2,0) Gauss-Jacobi rule). The first 'nqpt' points in the second dimension
-            // 'ir' have the same y-coordinates as those of the 1D rule for second dimension (i.e. (1,0)
-            // Gauss-Jacobi rule). The first 'nqpt' points in the third dimension have the same z-coordinates
-            // as those of the 1D rule for the third dimension (i.e. Gauss-Legendre rule). Additionally,
-            // the Bernstein PA algorithms expect evaluation of the component 1D bases at the Stroud nodes
-            // pulled back to the unit cube, so perform the pullback on the fly.
-            const real_t x = ir.IntPoint(i).x;
-            const real_t y = ir.IntPoint(nqpt*i).y / (1.0 - ir.IntPoint(nqpt*i).x);
-            const real_t z = ir.IntPoint(nqpt*nqpt*i).z / (1.0 - ir.IntPoint(
-                                                              nqpt*nqpt*i).x - ir.IntPoint(nqpt*nqpt*i).y);
-            Poly_1D::CalcBernstein(ndof-1, x, shape_a1);
-            Poly_1D::CalcBernstein(ndof-2, x, shape_Ga1);
-            for (int j = 0; j < ndof; j++)
-            {
-               rd2q->Ba1t[i+nqpt*j] = rd2q->Ba1[j+ndof*i] = shape_a1(j);
-               if (j < ndof-1)
-               {
-                  rd2q->Ga1t[i+nqpt*j] = rd2q->Ga1[j+(ndof-1)*i] = shape_Ga1(j);
-                  Poly_1D::CalcBernstein(ndof-2-j, y, shape_Ga2);
-               }
-
-               Poly_1D::CalcBernstein(ndof-1-j, y, shape_a2);
-               for (int k = 0; k < ndof-j; k++)
-               {
-                  const int a_2d_mass = rd2q->forward_map2d_mass[k + ndof*j];
-                  rd2q->Ba2t[i + nqpt*a_2d_mass] = rd2q->Ba2[a_2d_mass + basis_dim2d*i] =
-                                                      shape_a2(
-                                                         k);
-                  if (j < ndof-1 && k < ndof-j-1)
-                  {
-                     const int a_2d_diff = rd2q->forward_map2d_diff[k + (ndof-1)*j];
-                     rd2q->Ga2t[i + nqpt*a_2d_diff] = rd2q->Ga2[a_2d_diff + basis_dim2d_diff*i] =
-                                                         shape_Ga2(k);
-                     Poly_1D::CalcBernstein(ndof-2-j-k, z, shape_Ga3);
-                  }
-
-                  Poly_1D::CalcBernstein(ndof-1-j-k, z, shape_a3);
-                  for (int m = 0; m < ndof-j-k; m++)
-                  {
-                     const int a_3d_mass = rd2q->forward_map3d_mass[m + ndof*(k + ndof*j)];
-                     rd2q->Ba3t[i + nqpt*a_3d_mass] = rd2q->Ba3[a_3d_mass + basis_dim3d*i] =
-                                                         shape_a3(
-                                                            m);
-                     if (j < ndof-1 && k < ndof-j-1 && m < ndof-j-k-1)
-                     {
-                        // // collapsed 1D access
-                        // d2q->Ga3[i + nqpt*(m + d2q->offset3d[k + (ndof-1)*j])] = shape_Ga3(m);
-                        // collapsed 1D access with forward mapping
-                        const int a_3d_diff = rd2q->forward_map3d_diff[m + (ndof-1)*(k + (ndof-1)*j)];
-                        rd2q->Ga3t[i + nqpt*a_3d_diff] = rd2q->Ga3[a_3d_diff + basis_dim3d_diff*i] =
-                                                            shape_Ga3(m);
-                     }
-                  }
-               }
-            }
-         }
-
-         // stores the mapping from 3D Bernstein multi-index (i,j,k,p-i-j-k) to the
-         // lexicographic DOF ordering
-         int p = ndof - 1;
-         for (int i = 0; i < ndof; i++)
-         {
-            for (int j = 0; j < ndof-i; j++)
-            {
-               for (int k = 0; k < ndof-i-j; k++)
-               {
-                  int dof = (p+1)*(p+2)*(p+3) / 6;
-                  int tet = (p-k)*(p-k+1)*(p-k+2) / 6;
-                  int tri = (p+1-k-j)*(p+2-k-j)/2;
-                  int multi_idx = dof - tet - tri + i;
-                  rd2q->lex_map[k + ndof*(j + ndof*i)] = multi_idx;
-               }
-            }
-         }
-
-         dof2quad_array.Append(d2q);
-      }
-   }
-   return *d2q;
-}
-
 // static method
 void H1Pos_TetrahedronElement::CalcShape(
   const int p, const real_t l1, const real_t l2, const real_t l3,
@@ -191,21 +191,6 @@ public:
   /// Construct the H1Pos_TriangleElement of order @a p
   H1Pos_TriangleElement(const int p);

-   const DofToQuad &GetDofToQuad(const IntegrationRule &ir,
-                                 DofToQuad::Mode mode) const override
-   {
-      return (mode == DofToQuad::RAGGED_TENSOR) ?
-             GetRaggedTensorDofToQuad(*this, ir, mode, dof2quad_array) :
-             FiniteElement::GetDofToQuad(ir, mode);
-   }
-
-   static const DofToQuad &GetRaggedTensorDofToQuad(
-      const FiniteElement &fe, const IntegrationRule &ir,
-      DofToQuad::Mode mode,
-      Array<DofToQuad*> &dof2quad_array);
-
-   const Array<int> &GetDofMap() const { return dof_map; }
-
   // The size of shape is (p+1)(p+2)/2 (dof).
   static void CalcShape(const int p, const real_t x, const real_t y,
                         real_t *shape);
@@ -235,21 +220,6 @@ public:
   /// Construct the H1Pos_TetrahedronElement of order @a p
   H1Pos_TetrahedronElement(const int p);

-   const DofToQuad &GetDofToQuad(const IntegrationRule &ir,
-                                 DofToQuad::Mode mode) const override
-   {
-      return (mode == DofToQuad::RAGGED_TENSOR) ?
-             GetRaggedTensorDofToQuad(*this, ir, mode, dof2quad_array) :
-             FiniteElement::GetDofToQuad(ir, mode);
-   }
-
-   static const DofToQuad &GetRaggedTensorDofToQuad(
-      const FiniteElement &fe, const IntegrationRule &ir,
-      DofToQuad::Mode mode,
-      Array<DofToQuad*> &dof2quad_array);
-
-   const Array<int> &GetDofMap() const { return dof_map; }
-
   // The size of shape is (p+1)(p+2)(p+3)/6 (dof).
   static void CalcShape(const int p, const real_t x, const real_t y,
                         const real_t z, real_t *shape);
@@ -17,12 +17,6 @@
 namespace mfem
 {

-struct ScalarPyramid
-{
-   // Default basis type for H1 and L2 pyramids
-   static inline int DefaultType = 1; // Bergot(0) or Fuentes(1)
-};
-
 /** Base class for arbitrary order basis functions on pyramid-shaped elements

  This base class provides a common class to store temporary vectors,
@@ -228,19 +228,7 @@ FiniteElementCollection *FiniteElementCollection::New(const char *name)
   }
   else if (!strncmp(name, "H1_", 3))
   {
-      // Parse pyramid basis type if included in the name
-      const char *pyr = strstr(name, "Pyr");
-      if (pyr == NULL)
-      {
-         // Use default pyramid type elements
-         fec = new H1_FECollection(atoi(name + 7), atoi(name + 3));
-      }
-      else
-      {
-         // Use specific pyramid type elements
-         fec = new H1_FECollection(atoi(name + 7), atoi(name + 3),
-                                   BasisType::GaussLobatto, atoi(pyr + 3));
-      }
+      fec = new H1_FECollection(atoi(name + 7), atoi(name + 3));
   }
   else if (!strncmp(name, "H1Pos_Trace_", 12))
   {
@@ -257,44 +245,26 @@ FiniteElementCollection *FiniteElementCollection::New(const char *name)
   }
   else if (!strncmp(name, "H1@", 3))
   {
-      // Parse pyramid basis type if included in the name
-      const char *pyr = strstr(name, "Pyr");
-      if (pyr == NULL)
-      {
-         // Use default pyramid type elements
-         fec = new H1_FECollection(atoi(name + 9), atoi(name + 5),
-                                   BasisType::GetType(name[3]));
-      }
-      else
-      {
-         // Use specific pyramid type elements
-         fec = new H1_FECollection(atoi(name + 9), atoi(name + 5),
-                                   BasisType::GetType(name[3]),
-                                   atoi(pyr + 3));
-      }
+      fec = new H1_FECollection(atoi(name + 9), atoi(name + 5),
+                                BasisType::GetType(name[3]));
   }
-   else if (!strncmp(name, "L2", 2))
+   else if (!strncmp(name, "L2_T", 4))
+      fec = new L2_FECollection(atoi(name + 10), atoi(name + 6),
+                                atoi(name + 4));
+   else if (!strncmp(name, "L2_", 3))
   {
-      // Parse Map Type
-      const int mtype = strstr(name, "Int") == NULL ?
-                        FiniteElement::VALUE : FiniteElement::INTEGRAL;
-
-      // Parse the base order
-      const int p = atoi(strstr(name, "_P") + 2);
-
-      // Parse the mesh dimension
-      const int dim = atoi(strstr(name, "D") - 1);
-
-      // Parse basis type if specified
-      const char *t = strstr(name, "_T");
-      const int btype = t == NULL ? BasisType::GaussLegendre : atoi(t + 2);
-
-      // Parse the pyramid type if specified
-      const char *pyr = strstr(name, "Pyr");
-      const int ptype = pyr == NULL ? 1 : atoi(pyr + 3);
-
-      // Create collection
-      fec = new L2_FECollection(p, dim, btype, mtype, ptype);
+      fec = new L2_FECollection(atoi(name + 7), atoi(name + 3));
+   }
+   else if (!strncmp(name, "L2Int_T", 7))
+   {
+      fec = new L2_FECollection(atoi(name + 13), atoi(name + 9),
+                                atoi(name + 7), FiniteElement::INTEGRAL);
+   }
+   else if (!strncmp(name, "L2Int_", 6))
+   {
+      fec = new L2_FECollection(atoi(name + 10), atoi(name + 6),
+                                BasisType::GaussLegendre,
+                                FiniteElement::INTEGRAL);
   }
   else if (!strncmp(name, "RT_Trace_", 9))
   {
@@ -1739,10 +1709,9 @@ const int *RT1_3DFECollection::DofOrderForOrientation(Geometry::Type GeomType,


 H1_FECollection::H1_FECollection(const int p, const int dim, const int btype,
-                                 const int pyr_type)
+                                 const int pyrtype)
   : FiniteElementCollection(p)
   , dim(dim)
-   , p_type(pyr_type)
 {
   MFEM_VERIFY(p >= 1, "H1_FECollection requires order >= 1.");
   MFEM_VERIFY(dim >= 0 && dim <= 3, "H1_FECollection requires 0 <= dim <= 3.");
@@ -1755,14 +1724,7 @@ H1_FECollection::H1_FECollection(const int p, const int dim, const int btype,
   {
      case BasisType::GaussLobatto:
      {
-         if (pyr_type == ScalarPyramid::DefaultType)
-         {
-            snprintf(h1_name, 32, "H1_%dD_P%d", dim, p);
-         }
-         else
-         {
-            snprintf(h1_name, 32, "H1_%dD_P%d_Pyr%d", dim, p, pyr_type);
-         }
+         snprintf(h1_name, 32, "H1_%dD_P%d", dim, p);
         break;
      }
      case BasisType::Positive:
@@ -1948,11 +1910,11 @@ H1_FECollection::H1_FECollection(const int p, const int dim, const int btype,
         H1_dof[Geometry::TETRAHEDRON] = (TriDof*pm3)/3;
         H1_dof[Geometry::CUBE] = QuadDof*pm1;
         H1_dof[Geometry::PRISM] = TriDof*pm1;
-         if (pyr_type == 0 || b_type == BasisType::Positive)
+         if (pyrtype == 0 || b_type == BasisType::Positive)
         {
            H1_dof[Geometry::PYRAMID] = pm2*pm1*(2*p-3)/6; // Bergot (JSC)
         }
-         else if (pyr_type == 1)
+         else if (pyrtype == 1)
         {
            H1_dof[Geometry::PYRAMID] = pm1*pm1*pm1; // Fuentes
         }
@@ -1973,15 +1935,13 @@ H1_FECollection::H1_FECollection(const int p, const int dim, const int btype,
               new H1_TetrahedronElement(p, btype);
            H1_Elements[Geometry::CUBE] = new H1_HexahedronElement(p, btype);
            H1_Elements[Geometry::PRISM] = new H1_WedgeElement(p, btype);
-            if (pyr_type == 0)
+            if (pyrtype == 0)
            {
-               H1_Elements[Geometry::PYRAMID] =
-                  new H1_BergotPyramidElement(p, btype);
+               H1_Elements[Geometry::PYRAMID] = new H1_BergotPyramidElement(p, btype);
            }
            else
            {
-               H1_Elements[Geometry::PYRAMID] =
-                  new H1_FuentesPyramidElement(p, btype);
+               H1_Elements[Geometry::PYRAMID] = new H1_FuentesPyramidElement(p, btype);
            }
         }

@@ -2188,7 +2148,6 @@ L2_FECollection::L2_FECollection(const int p, const int dim, const int btype,
   : FiniteElementCollection(p)
   , dim(dim)
   , m_type(map_type)
-   , p_type(pyr_type)
 {
   MFEM_VERIFY(p >= 0, "L2_FECollection requires order >= 0.");

@@ -2204,25 +2163,10 @@ L2_FECollection::L2_FECollection(const int p, const int dim, const int btype,
   switch (btype)
   {
      case BasisType::GaussLegendre:
-         if (pyr_type == ScalarPyramid::DefaultType)
-         {
-            snprintf(d_name, 32, "%s_%dD_P%d", prefix, dim, p);
-         }
-         else
-         {
-            snprintf(d_name, 32, "%s_%dD_P%d_Pyr%d", prefix, dim, p, pyr_type);
-         }
+         snprintf(d_name, 32, "%s_%dD_P%d", prefix, dim, p);
         break;
      default:
-         if (pyr_type == ScalarPyramid::DefaultType)
-         {
-            snprintf(d_name, 32, "%s_T%d_%dD_P%d", prefix, btype, dim, p);
-         }
-         else
-         {
-            snprintf(d_name, 32, "%s_T%d_%dD_P%d_Pyr%d",
-                     prefix, btype, dim, p, pyr_type);
-         }
+         snprintf(d_name, 32, "%s_T%d_%dD_P%d", prefix, btype, dim, p);
   }

   for (int g = 0; g < Geometry::NumGeom; g++)
@@ -2341,13 +2285,11 @@ L2_FECollection::L2_FECollection(const int p, const int dim, const int btype,
         L2_Elements[Geometry::PRISM] = new L2_WedgeElement(p, btype);
         if (pyr_type == 0)
         {
-            L2_Elements[Geometry::PYRAMID] =
-               new L2_BergotPyramidElement(p, btype);
+            L2_Elements[Geometry::PYRAMID] = new L2_BergotPyramidElement(p, btype);
         }
         else
         {
-            L2_Elements[Geometry::PYRAMID] =
-               new L2_FuentesPyramidElement(p, btype);
+            L2_Elements[Geometry::PYRAMID] = new L2_FuentesPyramidElement(p, btype);
         }
      }

@@ -100,10 +100,6 @@ public:
      return FiniteElementForGeometry(GeomType);
   }

-   /** @brief Returns a collection of the trace elements.
-
-       @note The collection is owned by the caller and is NOT deleted in the
-       destructor. */
   virtual FiniteElementCollection *GetTraceCollection() const;

   virtual ~FiniteElementCollection();
@@ -254,14 +250,6 @@ public:
       its GetOrder() method. */
   virtual FiniteElementCollection *Clone(int p) const;

-   /** @brief Return the order parameter used to construct this collection.
-    *  This differs from GetOrder() depending on the collection type. */
-   virtual int GetConstructorOrder() const
-   {
-      MFEM_ABORT("Collection " << Name() << " does not support GetConstructorOrder");
-      return -1;
-   }
-
 protected:
   const int base_p; ///< Order as returned by GetOrder().

@@ -290,7 +278,7 @@ protected:
 class H1_FECollection : public FiniteElementCollection
 {
 protected:
-   int dim, b_type, p_type;
+   int dim, b_type;
   char h1_name[32];
   FiniteElement *H1_Elements[Geometry::NumGeom];
   int H1_dof[Geometry::NumGeom];
@@ -299,7 +287,7 @@ protected:
 public:
   explicit H1_FECollection(const int p, const int dim = 3,
                            const int btype = BasisType::GaussLobatto,
-                            const int pyr_type = ScalarPyramid::DefaultType);
+                            const int pyrtype = 1);

   const FiniteElement *
   FiniteElementForGeometry(Geometry::Type GeomType) const override;
@@ -324,10 +312,7 @@ public:
   const int *GetDofMap(Geometry::Type GeomType, int p) const;

   FiniteElementCollection *Clone(int p) const override
-   { return new H1_FECollection(p, dim, b_type, p_type); }
-
-   int GetConstructorOrder() const override
-   { return base_p; }
+   { return new H1_FECollection(p, dim, b_type); }

   virtual ~H1_FECollection();
 };
@@ -358,10 +343,6 @@ class H1_Trace_FECollection : public H1_FECollection
 public:
   H1_Trace_FECollection(const int p, const int dim,
                         const int btype = BasisType::GaussLobatto);
-
-   FiniteElementCollection *Clone(int p) const override
-   { return new H1_Trace_FECollection(p, dim+1, b_type); }
-
 };

 /// Arbitrary order "L2-conforming" discontinuous finite elements.
@@ -371,7 +352,6 @@ private:
   int dim;
   int b_type; // BasisType
   int m_type; // map type
-   int p_type; // Pyramid type (0 -> Bergot, 1 -> Fuentes)
   char d_name[32];
   ScalarFiniteElement *L2_Elements[Geometry::NumGeom];
   ScalarFiniteElement *Tr_Elements[Geometry::NumGeom];
@@ -384,7 +364,7 @@ public:
   L2_FECollection(const int p, const int dim,
                   const int btype = BasisType::GaussLegendre,
                   const int map_type = FiniteElement::VALUE,
-                   const int pyr_type = ScalarPyramid::DefaultType);
+                   const int pyrtype = 1);

   const FiniteElement *
   FiniteElementForGeometry(Geometry::Type GeomType) const override;
@@ -414,10 +394,7 @@ public:
   int GetBasisType() const { return b_type; }

   FiniteElementCollection *Clone(int p) const override
-   { return new L2_FECollection(p, dim, b_type, m_type, p_type); }
-
-   int GetConstructorOrder() const override
-   { return base_p; }
+   { return new L2_FECollection(p, dim, b_type, m_type); }

   virtual ~L2_FECollection();
 };
@@ -479,9 +456,6 @@ public:
   FiniteElementCollection *Clone(int p) const override
   { return new RT_FECollection(p, dim, cb_type, ob_type); }

-   int GetConstructorOrder() const override
-   { return base_p-1; }
-
   virtual ~RT_FECollection();
 };

@@ -562,9 +536,6 @@ public:
   FiniteElementCollection *Clone(int p) const override
   { return new ND_FECollection(p, dim, cb_type, ob_type); }

-   int GetConstructorOrder() const override
-   { return dim>1 ? base_p : base_p+1; }
-
   virtual ~ND_FECollection();
 };

@@ -577,9 +548,6 @@ public:
   ND_Trace_FECollection(const int p, const int dim,
                         const int cb_type = BasisType::GaussLobatto,
                         const int ob_type = BasisType::GaussLegendre);
-
-   FiniteElementCollection *Clone(int p) const override
-   { return new ND_Trace_FECollection(p, dim+1, cb_type, ob_type); }
 };

 /// Arbitrary order 3D H(curl)-conforming Nedelec finite elements in 1D.
@@ -62,10 +62,6 @@
 #include "pnonlinearform.hpp"
 #endif

-#ifdef MFEM_USE_SIDRE
-#include "sidredatacollection.hpp"
-#endif
-
 #ifdef MFEM_USE_CONDUIT
 #include "conduitdatacollection.hpp"
 #endif
@@ -4631,8 +4631,9 @@ FiniteElementCollection *FiniteElementSpace::Load(Mesh *m, std::istream &input)

 ElementDofOrdering GetEVectorOrdering(const FiniteElementSpace& fes)
 {
-   return (UsesTensorBasis(fes) || fes.UsesRaggedTensorBasis()) ?
+   return UsesTensorBasis(fes)?
          ElementDofOrdering::LEXICOGRAPHIC:
          ElementDofOrdering::NATIVE;
 }
+
 } // namespace mfem
@@ -1514,18 +1514,6 @@ public:
      return dynamic_cast<const L2_FECollection*>(fec) != NULL;
   }

-   /// @brief Return true if the mesh contains only one topology, the elements are
-   /// all triangles or tetrahedrons, and the elements are ragged tensor elements
-   /// i.e. Bernstein/positive basis.
-   bool UsesRaggedTensorBasis() const
-   {
-      bool simplex = this->GetMesh()->IsSimplexMesh();
-      bool positive =
-         dynamic_cast<const mfem::H1Pos_TriangleElement *>(this->GetTypicalFE()) ||
-         dynamic_cast<const mfem::H1Pos_TetrahedronElement *>(this->GetTypicalFE());
-      return simplex && positive;
-   }
-
   /** In variable-order spaces on nonconforming (NC) meshes, this function
       controls whether strict conformity is enforced in cases where coarse
       edges/faces have higher polynomial order than their fine NC neighbors.
@@ -2256,104 +2256,6 @@ void GridFunction::AccumulateAndCountBdrTangentValues(
   }
 }

-void GridFunction::AccumulateAndCountTraceValues(
-   Coefficient *coeff[], VectorCoefficient *vcoeff,
-   Array<int> &values_counter)
-{
-   if (vcoeff)
-   {
-      MFEM_VERIFY(fes->GetVDim() == vcoeff->GetVDim(),
-                  "vcoeff vdim != fes VDim");
-      MFEM_VERIFY(fes->GetTypicalTraceElement()->GetMapType() ==
-                  FiniteElement::VALUE &&
-                  fes->GetTypicalTraceElement()->GetRangeType() ==
-                  FiniteElement::SCALAR,
-                  "Can only call ProjectTraceCoefficient on scalar value-type "
-                  "trace elements. "
-                  "Use ProjectTraceCoefficientNormal for RT and "
-                  "ProjectTraceCoefficientTangent for ND finite elements.");
-   }
-
-   Array<int> vdofs;
-   Vector vc;
-
-   values_counter.SetSize(Size());
-   values_counter = 0;
-
-   const int vdim = fes->GetVDim();
-   HostReadWrite();
-
-   for (int i = 0; i < fes->GetMesh()->GetNumFaces(); i++)
-   {
-
-      const FiniteElement *fe = fes->GetFaceElement(i);
-      const int fdof = fe->GetDof();
-      ElementTransformation *transf = fes->GetMesh()->GetFaceTransformation(i);
-      const IntegrationRule &ir = fe->GetNodes();
-      fes->GetFaceVDofs(i, vdofs);
-
-      for (int j = 0; j < fdof; j++)
-      {
-         const IntegrationPoint &ip = ir.IntPoint(j);
-         transf->SetIntPoint(&ip);
-         if (vcoeff) { vcoeff->Eval(vc, *transf, ip); }
-         for (int d = 0; d < vdim; d++)
-         {
-            if (!vcoeff && !coeff[d]) { continue; }
-
-            real_t val = vcoeff ? vc(d) : coeff[d]->Eval(*transf, ip);
-            int ind = vdofs[fdof*d+j];
-            if ( ind < 0 )
-            {
-               val = -val, ind = -1-ind;
-            }
-            if (++values_counter[ind] == 1)
-            {
-               (*this)(ind) = val;
-            }
-            else
-            {
-               (*this)(ind) += val;
-            }
-         }
-      }
-   }
-}
-
-void GridFunction::AccumulateAndCountTraceTangentValues(
-   VectorCoefficient &vcoeff, Array<int> &values_counter)
-{
-   MFEM_VERIFY(fes->GetVDim() == 1, "fespace VDim != 1");
-   MFEM_VERIFY(fes->GetTypicalTraceElement()
-               ->GetRangeType() == FiniteElement::VECTOR &&
-               fes->GetTypicalTraceElement()
-               ->GetMapType() == FiniteElement::H_CURL,
-               "Not an ND FE space!");
-   MFEM_VERIFY(fes->GetTypicalTraceElement()->GetPhysRangeDim(
-                  fes->GetMesh()->SpaceDimension()) == vcoeff.GetVDim(),
-               "vcoeff vdim != PhysRangeDim");
-
-   const FiniteElement *fe;
-   ElementTransformation *T;
-   Array<int> dofs;
-   Vector lvec;
-
-   values_counter.SetSize(Size());
-   values_counter = 0;
-
-   HostReadWrite();
-
-   for (int i = 0; i < fes->GetMesh()->GetNumFaces(); i++)
-   {
-      fe = fes->GetFaceElement(i);
-      T = fes->GetMesh()->GetFaceTransformation(i);
-      fes->GetFaceVDofs(i, dofs);
-      lvec.SetSize(fe->GetDof());
-      fe->Project(vcoeff, *T, lvec);
-      accumulate_dofs(dofs, lvec, *this, values_counter);
-   }
-}
-
 void GridFunction::ComputeMeans(AvgType type, Array<int> &zones_per_vdof)
 {
   switch (type)
@@ -2796,74 +2698,6 @@ void GridFunction::ProjectCoefficient(VectorCoefficient &vcoeff,
   }
 }

-void GridFunction::ProjectTraceCoefficient(Coefficient *coeff[])
-{
-   Array<int> values_counter;
-   AccumulateAndCountTraceValues(coeff, NULL, values_counter);
-   ComputeMeans(ARITHMETIC, values_counter);
-}
-
-void GridFunction::ProjectTraceCoefficient(Coefficient &coeff)
-{
-   MFEM_VERIFY(FESpace()->GetVDim() == 1, "ProjectTraceCoefficient(Coefficient&)"
-               "is only valid for scalar GridFunction");
-   Coefficient *coeff_p = &coeff;
-   ProjectTraceCoefficient(&coeff_p);
-}
-
-void GridFunction::ProjectTraceCoefficient(VectorCoefficient &vcoeff)
-{
-   MFEM_VERIFY(FESpace()->GetVDim() == vcoeff.GetVDim(),
-               "Incompatible vcoeff vdim and fes vdim");
-   Array<int> values_counter;
-   AccumulateAndCountTraceValues(NULL, &vcoeff, values_counter);
-   ComputeMeans(ARITHMETIC, values_counter);
-}
-
-void GridFunction::ProjectTraceCoefficientNormal(VectorCoefficient &vcoeff)
-{
-   MFEM_VERIFY(fes->GetVDim() == 1, "fespace VDim != 1");
-   MFEM_VERIFY(fes->GetTypicalTraceElement()->GetRangeType() ==
-               FiniteElement::SCALAR &&
-               fes->GetTypicalTraceElement()->GetMapType() ==
-               FiniteElement::INTEGRAL, "Not an RT FE space!");
-   MFEM_VERIFY(vcoeff.GetVDim() == fes->GetMesh()->SpaceDimension(),
-               "vcoeff vdim (" << vcoeff.GetVDim()
-               << ") != SpaceDimension ("
-               << fes->GetMesh()->SpaceDimension() << ")");
-
-   const FiniteElement *fe;
-   ElementTransformation *T;
-   Array<int> dofs;
-   int dim = vcoeff.GetVDim();
-   Vector vc(dim), nor(dim), lvec;
-
-   for (int i = 0; i < fes->GetMesh()->GetNumFaces(); i++)
-   {
-      fe = fes->GetFaceElement(i);
-      T = fes->GetMesh()->GetFaceTransformation(i);
-      const IntegrationRule &ir = fe->GetNodes();
-      lvec.SetSize(fe->GetDof());
-      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);
-         lvec(j) = (vc * nor);
-      }
-      fes->GetFaceVDofs(i, dofs);
-      SetSubVector(dofs, lvec);
-   }
-}
-
-void GridFunction::ProjectTraceCoefficientTangent(VectorCoefficient &vcoeff)
-{
-   Array<int> values_counter;
-   AccumulateAndCountTraceTangentValues(vcoeff, values_counter);
-   ComputeMeans(ARITHMETIC, values_counter);
-}
-
 void GridFunction::ProjectCoefficientGlobalL2(VectorCoefficient &vcoeff,
                                              real_t rtol, int iter)
 {
@@ -3196,14 +3030,10 @@ void GridFunction::ProjectCoefficient(Coefficient *coeff[])
   }
 }

-void GridFunction::ProjectDiscCoefficient(
-   std::variant<Coefficient*, VectorCoefficient*> coeff, Array<int> &dof_attr)
+void GridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff,
+                                          Array<int> &dof_attr)
 {
-   std::visit([&](auto* c)
-   {
-      MFEM_VERIFY(VectorDim() == c->GetVDim(), "coeff vdim != VectorDim()");
-   }, coeff);
-
+   MFEM_VERIFY(VectorDim() == coeff.GetVDim(), "coeff vdim != VectorDim()");
   Array<int> vdofs;
   Vector vals;

@@ -3217,10 +3047,7 @@ void GridFunction::ProjectDiscCoefficient(
   {
      fes->GetElementVDofs(i, vdofs);
      vals.SetSize(vdofs.Size());
-      std::visit([&](auto* c)
-      {
-         fes->GetFE(i)->Project(*c, *fes->GetElementTransformation(i), vals);
-      }, coeff);
+      fes->GetFE(i)->Project(coeff, *fes->GetElementTransformation(i), vals);

      // the values in shared dofs are determined from the element with maximal
      // attribute
@@ -3236,6 +3063,13 @@ void GridFunction::ProjectDiscCoefficient(
   }
 }

+void GridFunction::ProjectDiscCoefficient(VectorCoefficient &coeff)
+{
+   MFEM_VERIFY(VectorDim() == coeff.GetVDim(), "coeff vdim != VectorDim()");
+   Array<int> dof_attr;
+   ProjectDiscCoefficient(coeff, dof_attr);
+}
+
 void GridFunction::ProjectDiscCoefficient(Coefficient &coeff, AvgType type)
 {
   // Harmonic  (x1 ... xn) = [ (1/x1 + ... + 1/xn) / n ]^-1.
@@ -5452,7 +5286,6 @@ PLBound GridFunction::GetBounds(Vector &lower, Vector &upper,
 {
   int max_order = fes->GetMaxElementOrder();
   PLBound plb(fes, ref_factor*(max_order+1));
-
   Vector lel, uel;
   GetElementBounds(plb, lel, uel, vdim);

@@ -23,7 +23,6 @@
 #include <limits>
 #include <ostream>
 #include <string>
-#include <variant>

 namespace mfem
 {
@@ -80,18 +79,10 @@ protected:
                        bool wcoef,
                        int subdomain);

-   /** @brief Project a discontinuous (vector) coefficient as a grid function on
-       a continuous finite element space. Return in dof_attr the maximal
-       attribute of the elements containing each degree of freedom. */
-   virtual void ProjectDiscCoefficient(
-      std::variant<Coefficient*, VectorCoefficient*> coeff, Array<int> &dof_attr);
-
-   /** @brief Project a discontinuous (vector) coefficient as a grid function on
-       a continuous finite element space. The values in shared dofs are
-       determined from the element with maximal attribute. */
-   virtual void ProjectDiscCoefficient(
-      std::variant<Coefficient*, VectorCoefficient*> coeff)
-   { Array<int> dof_attr; ProjectDiscCoefficient(coeff, dof_attr); };
+   /** Project a discontinuous vector coefficient in a continuous space and
+       return in dof_attr the maximal attribute of the elements containing each
+       degree of freedom. */
+   void ProjectDiscCoefficient(VectorCoefficient &coeff, Array<int> &dof_attr);

   /** Helper function for ProjectCoefficientElementL2 */
   void ProjectCoefficientElementL2_(Coefficient &coeff, Vector &sol, Vector &Va);
@@ -524,17 +515,10 @@ public:
       but using an array of scalar coefficients for each component. */
   void ProjectCoefficient(Coefficient *coeff[]);

-   /** @brief Project a discontinuous coefficient as a grid function on
-       a continuous finite element space. The values in shared dofs are
-       determined from the element with maximal attribute. */
-   virtual void ProjectDiscCoefficient(Coefficient &coeff)
-   { ProjectDiscCoefficient(&coeff); }
-
   /** @brief Project a discontinuous vector coefficient as a grid function on
       a continuous finite element space. The values in shared dofs are
       determined from the element with maximal attribute. */
-   virtual void ProjectDiscCoefficient(VectorCoefficient &coeff)
-   { ProjectDiscCoefficient(&coeff); }
+   virtual void ProjectDiscCoefficient(VectorCoefficient &coeff);

   enum AvgType {ARITHMETIC, HARMONIC};
   /** @brief Projects a discontinuous coefficient so that the values in shared
@@ -578,13 +562,6 @@ protected:
                                           const Array<int> &bdr_attr,
                                           Array<int> &values_counter);

-   void AccumulateAndCountTraceValues(Coefficient *coeff[],
-                                      VectorCoefficient *vcoeff,
-                                      Array<int> &values_counter);
-
-   void AccumulateAndCountTraceTangentValues(VectorCoefficient &vcoeff,
-                                             Array<int> &values_counter);
-
   // Complete the computation of averages; called e.g. after
   // AccumulateAndCountZones().
   void ComputeMeans(AvgType type, Array<int> &zones_per_vdof);
@@ -670,23 +647,6 @@ public:
      ProjectBdrCoefficient(&coeff_p, attr);
   }

-   /// Project a Coefficient on a GridFunction defined on H1 trace space
-   void ProjectTraceCoefficient(Coefficient *coeff[]);
-   void ProjectTraceCoefficient(Coefficient &coeff);
-
-   /** @brief Project a VectorCoefficient @a vcoeff on a GridFunction
-       defined on a Vector H1 trace space. Note that this also works
-       for a scalar H1 trace space, where only the first component of
-       @a vcoeff is used. */
-   void ProjectTraceCoefficient(VectorCoefficient &vcoeff);
-   /** @brief Project a VectorCoefficient on a GridFunction
-       defined on an RT trace space */
-   void ProjectTraceCoefficientNormal(VectorCoefficient &vcoeff);
-   /** @brief Project a VectorCoefficient on a GridFunction
-       defined on an ND trace space */
-   void ProjectTraceCoefficientTangent(VectorCoefficient &vcoeff);
-
-
   /** @brief Project a VectorCoefficient on the GridFunction, modifying only
       DOFs on the boundary associated with the boundary attributes marked in
       the @a attr array. */
@@ -1791,8 +1751,8 @@ public:
                             const int ref_factor=1, const int vdim=-1) const;

   /// Computes the \ref PLBound for the gridfunction with number of control
-   /// points based on @a ref_factor, and returns the bounds for each element
-   /// ordered byNODES:
+   /// points based on \p ref_factor, and returns the bounds for each element
+   /// ordered byNodes:
   /// lower_{0,0}, lower_{1,0}, ..., lower_{ne-1,0},
   /// lower_{0,1}, ..., lower_{ne-1,vdim-1}. We also return the
   /// PLBound object used to compute the bounds.
@@ -1826,7 +1786,7 @@ public:
                         const int vdim = -1) const;

   /// Compute bounds on the grid function for all the elements. The bounds
-   /// are returned in @b lower and @b upper, ordered byNODES:
+   /// are returned in @b lower and @b upper, ordered byNodes:
   /// lower_{0,0}, lower_{1,0}, ..., lower_{ne-1,0},
   /// lower_{0,1}, ..., lower_{ne-1,vdim-1}
   void GetElementBounds(const PLBound &plb, Vector &lower, Vector &upper,
--- a/Show More
+++ b/Show More
Author	SHA1	Message	Date
Kenneth Weiss	b86fb5308a	style	2026-04-14 08:54:23 -07:00
Kenneth Weiss	4463d921f6	Updates CHANGELOG	2026-04-13 21:52:37 -07:00
Kenneth Weiss	20c657559d	More fixes for building contact mini-app with make and cmake	2026-04-13 21:47:15 -07:00
Kenneth Weiss	7a999efb0d	yAdds back FindAxom in support of contact mini-app	2026-04-13 18:54:17 -07:00
Kenneth Weiss	442ac18b65	Removes references to Axom in INSTALL	2026-04-13 12:40:24 -07:00
Kenneth Weiss	5097a44411	Removes Axom from make-based and cmake-based build systems	2026-04-13 12:37:18 -07:00
Kenneth Weiss	e84fef4191	Removes SidreDataCollection from MFEM This has been superseded by the MFEMSidreDataCollection in Axom.	2026-04-13 12:34:02 -07:00