Fix single precision build.

Templates for some reducers.
More templates.
2025-12-12 18:29:30 -08:00 · 2025-12-12 17:09:03 -08:00 · 2025-12-12 13:57:12 -08:00 · 2025-12-12 13:25:10 -08:00 · 2025-12-12 12:35:54 -08:00 · 2025-12-12 12:34:20 -08:00
385 changed files with 9971 additions and 43905 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)
@@ -25,7 +25,7 @@ runs:
  steps:
    - uses: ./.github/actions/sanitize/config

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

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

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

-    - uses: actions/cache/restore@v5  # Cache for Metis
+    - uses: actions/cache/restore@v4  # Cache for Metis
      if: ${{inputs.par == 'true'}}
      with:
        path: ${{env.METIS_DIR}}
        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'}}
      run: ln -s -f ${{env.HYPRE_DIR}} hypre && ln -s -f ${{env.METIS_DIR}} metis-4.0
      shell: bash

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

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

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

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

      # 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:
@@ -318,13 +329,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'
@@ -369,16 +374,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
@@ -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 }}
@@ -50,11 +40,11 @@ jobs:

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

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

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

      - name: Perform CodeQL Analysis
-        uses: github/codeql-action/analyze@v4
+        uses: github/codeql-action/analyze@v2
@@ -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@v4
+        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@v4
+        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@v4
+        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
@@ -38,7 +38,7 @@ jobs:
       github.event.pull_request.head.repo.full_name != github.repository)
    steps:
      - name: checkout mfem
-        uses: actions/checkout@v6
+        uses: actions/checkout@v4

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

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

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

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

 miniapps/performance/ex1
 miniapps/performance/ex1p
@@ -369,7 +360,6 @@ miniapps/shifted/lsf_integral
 miniapps/tools/display-basis
 miniapps/tools/load-dc
 miniapps/tools/convert-dc
-miniapps/tools/compare-dc
 miniapps/tools/gridfunction-bounds
 miniapps/tools/lor-transfer
 miniapps/tools/plor-transfer
@@ -443,10 +433,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
@@ -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
@@ -11,87 +11,6 @@
 Version 4.9.1 (development)
 ===========================

- Added policy for AI-assisted contribution to CONTRIBUTING.md.
-
-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.
-
- 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.
-
- 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,
-  Vector and VectorFE, also NURBS versions. Optionally different types of
-  projections can be selected, default behavior 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
-  recursive subdivision.
-
- Extend FindPointsGSLIB to support surface meshes.
-
-Meshing improvements
--------------------
- 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
-------------
- 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.
-
-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.
-
-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.
-

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

-  Added 'GetResiduals' and 'GetFinalAbsResidualNorm' to 'HyprePCG',
-  'HypreGMRES', and 'HypreFGMRES' to get 'r' and '|r|_p'. Note that the latter
-  computes '|r|_p' from 'r' instead of returning a cached value like the
-  relative 'GetFinalResidualNorm'. These require Hypre >= 2.15.0.
-
-  Changed the default solver parameters for 'HyprePCG' to 'tol=1e-6' and
-  'max_iter=1000'. This matches the default parameters in Hypre 3.0.
-
-  Added various helper functions for querying/modifying Hypre solvers:
-  'HypreSmoother::GetType', 'HypreSmoother::GetSOROptions',
-  'HypreSmoother::GetPolyOptions', 'HypreSmoother::GetWindowParameters',
-  'HypreSmoother::IsOperatorSymmetric', 'HyprePCG::GetTol',
-  'HyprePCG::GetAbsTol', 'HyprePCG::GetMaxIter', 'HyprePCG::SetUseTwoNorm',
-  'HypreGMRES::GetTol', 'HypreGMRES::GetAbsTol', 'HypreGMRES::GetMaxIter',
-  'HypreGMRES::GetKDim', 'HypreFGMRES::GetTol', 'HypreFGMRES::GetMaxIter',
-  'HypreFGMRES::GetKDim', and 'HypreBoomerAMG::GetMaxIter'.
-
 GPU computing
 -------------
 - Added the 'gpu', 'raja-gpu', and 'ceed-gpu' backend aliases/shortcuts which
@@ -433,15 +433,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)
@@ -640,7 +631,7 @@ find_package(Threads REQUIRED)
 set(MFEM_TPLS OPENMP HYPRE LAPACK BLAS SuperLUDist STRUMPACK METIS SuiteSparse
    SUNDIALS PETSC SLEPC MUMPS AXOM FMS CONDUIT Ginkgo GNUTLS GSLIB HDF5
    NETCDF MPFR PUMI HIOP POSIXCLOCKS MFEMBacktrace ZLIB OCCA CEED RAJA UMPIRE
-    ADIOS2 MKL_CPARDISO MKL_PARDISO AMGX MAGMA CUSPARSE CUBLAS 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)

@@ -661,8 +652,6 @@ foreach(TPL IN LISTS MFEM_TPLS)
  endif()
 endforeach(TPL)

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

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

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

-# Install the data directory if present, i.e. if the copy_data target is built
-install(DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/data
-  DESTINATION ${MFEM_INSTALL_DIR} OPTIONAL)
-
 #-------------------------------------------------------------------------------
 # Create 'config.mk' from 'config.mk.in' for the build and install locations and
 # define install rules for 'config.mk' and 'test.mk'
 #-------------------------------------------------------------------------------
-if (MFEM_USE_CUDA OR MFEM_USE_HIP)
-  option(MFEM_EXPORT_GPU_CONFIG "Export config.mk for GPU-enabled downstream packages" ON)
-endif()
+
 mfem_export_mk_files()
@@ -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
@@ -560,13 +554,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,24 +719,13 @@ 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).
-  Versions: Ginkgo >= 1.9.0. When building Ginkgo with distributed support, a
-            recent version of the "develop" branch is required (1.11 as defined
-            in include/ginkgo/config.hpp).
+  Versions: Ginkgo >= 1.9.0.

 - AmgX (optional), used when MFEM_USE_AMGX = YES.
  URL: https://github.com/NVIDIA/AMGX
@@ -808,7 +791,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 +867,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 +955,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
@@ -1217,7 +1200,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
@@ -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@)
@@ -110,14 +109,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

@@ -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()
@@ -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")
@@ -701,6 +701,7 @@ endfunction(mfem_find_library)
 # Extract compile and link options needed by the given target.
 #
 function(mfem_get_target_options Target CompileOptsVar LinkOptsVar)
+
  if (NOT TARGET ${Target})
    return()
  endif()
@@ -798,12 +799,7 @@ function(mfem_get_target_options Target CompileOptsVar LinkOptsVar)
          # message(STATUS "Lib = ${Lib}")
          # Filter-out generator expressions
          if (NOT ("${Lib}" MATCHES "^\\$"))
-            if(NOT ("${Lib}" STREQUAL "dl"))
-              list(APPEND LinkOpts "${Lib}")
-            else()
-              # for some reason libdl doesn't include the "-l"
-              list(APPEND LinkOpts "-ldl")
-            endif()
+            list(APPEND LinkOpts "${Lib}")
          endif()
        else()
          mfem_get_target_options(${Lib} COpts LOpts)
@@ -892,18 +888,9 @@ function(mfem_export_mk_files)
      set(${var} NO)
    endif()
  endforeach()
-  if (MFEM_USE_CUDA AND MFEM_EXPORT_GPU_CONFIG)
-    set(MFEM_CXX ${CMAKE_CUDA_COMPILER})
-    if(MFEM_CUDA_COMPILER_IS_NVCC)
-      set(MFEM_HOST_CXX ${CMAKE_CUDA_HOST_COMPILER})
-    else()
-      set(MFEM_HOST_CXX ${CMAKE_CXX_COMPILER})
-    endif()
-  else()
-    # mfem doesn't use enable_language(HIP)
-    set(MFEM_CXX ${CMAKE_CXX_COMPILER})
-    set(MFEM_HOST_CXX ${CMAKE_CXX_COMPILER})
-  endif()
+  # TODO: Add support for MFEM_USE_CUDA=YES
+  set(MFEM_CXX ${CMAKE_CXX_COMPILER})
+  set(MFEM_HOST_CXX ${MFEM_CXX})
  set(MFEM_CPPFLAGS "")
  get_target_property(cxx_std mfem CXX_STANDARD)
  # For now, we ignore the setting of the CXX_EXTENSIONS property. If this
@@ -913,50 +900,6 @@ function(mfem_export_mk_files)
  string(STRIP
         "${cxx_std_flag} ${CMAKE_CXX_FLAGS_${BUILD_TYPE}} ${CMAKE_CXX_FLAGS}"
         MFEM_CXXFLAGS)
-  if(MFEM_EXPORT_GPU_CONFIG)
-    if (MFEM_USE_CUDA)
-      set(MFEM_CXXFLAGS "${MFEM_CXXFLAGS} ${CMAKE_CUDA_FLAGS}")
-      if (MFEM_CUDA_COMPILER_IS_NVCC)
-        set(MFEM_CXXFLAGS "-x=cu ${MFEM_CXXFLAGS} -ccbin ${CMAKE_CXX_COMPILER} --forward-unknown-to-host-compiler")
-        # The following intentionally hides CUDA deprecation warnings
-        foreach(ENTRY IN LISTS CUDAToolkit_INCLUDE_DIRS)
-          set(MFEM_CXXFLAGS "${MFEM_CXXFLAGS} -isystem ${ENTRY}")
-        endforeach()
-        if (CMAKE_VERSION VERSION_GREATER_EQUAL 3.18.0)
-          # architecture flags not part of CMAKE_CUDA_FLAGS
-          if ("all" STREQUAL "${CMAKE_CUDA_ARCHITECTURES}"
-              OR "native" STREQUAL "${CMAKE_CUDA_ARCHITECTURES}"
-              OR "all-major" STREQUAL "${CMAKE_CUDA_ARCHITECTURES}")
-            set(MFEM_CXXFLAGS "${MFEM_CXXFLAGS} -arch=${CMAKE_CUDA_ARCHITECTURES}")
-          else()
-            foreach (ENTRY IN LISTS CMAKE_CUDA_ARCHITECTURES)
-              set(MFEM_CXXFLAGS
-                "${MFEM_CXXFLAGS} -gencode arch=compute_${ENTRY},code=sm_${ENTRY}")
-            endforeach()
-          endif()
-        endif()
-      else()
-        set(MFEM_CXXFLAGS "${MFEM_CXXFLAGS} -xcuda --cuda-path=${CUDAToolkit_LIBRARY_ROOT}")
-        if (CMAKE_VERSION VERSION_GREATER_EQUAL 3.18.0)
-          # architecture flags not part of CMAKE_CUDA_FLAGS
-          if ("all" STREQUAL "${CMAKE_CUDA_ARCHITECTURES}"
-              OR "native" STREQUAL "${CMAKE_CUDA_ARCHITECTURES}"
-              OR "all-major" STREQUAL "${CMAKE_CUDA_ARCHITECTURES}")
-            # TODO: not supported
-          else()
-            foreach(ENTRY IN LISTS CMAKE_CUDA_ARCHITECTURES)
-              set(MFEM_CXXFLAGS "-cuda-gpu-arch=sm_${ENTRY} ${MFEM_CXXFLAGS}")
-            endforeach()
-          endif()
-        endif()
-      endif()
-    elseif (MFEM_USE_HIP)
-      set(MFEM_CXXFLAGS "${MFEM_CXXFLAGS} -xhip")
-      foreach(ENTRY IN LISTS CMAKE_HIP_ARCHITECTURES)
-        set(MFEM_CXXFLAGS "--offload-arch=${ENTRY} ${MFEM_CXXFLAGS}")
-      endforeach()
-    endif()
-  endif()
  set(MFEM_TPLFLAGS "")
  foreach(dir ${TPL_INCLUDE_DIRS})
    set(MFEM_TPLFLAGS "${MFEM_TPLFLAGS} -I${dir}")
@@ -987,9 +930,6 @@ function(mfem_export_mk_files)
    set(MFEM_SHARED NO)
    set(MFEM_STATIC YES)
  endif()
-  if (MFEM_USE_CUDA)
-    set(MFEM_EXT_LIBS "${MFEM_EXT_LIBS} -lcudart")
-  endif()
  set(MFEM_BUILD_TAG "${CMAKE_SYSTEM}")
  set(MFEM_PREFIX "${CMAKE_INSTALL_PREFIX}")
  # For the next 4 variables, these are the values for the build-tree version of
@@ -998,15 +938,8 @@ function(mfem_export_mk_files)
  set(MFEM_LIB_DIR "${PROJECT_BINARY_DIR}")
  set(MFEM_TEST_MK "${PROJECT_SOURCE_DIR}/config/test.mk")
  set(MFEM_CONFIG_EXTRA "MFEM_BUILD_DIR ?= ${PROJECT_BINARY_DIR}")
-  if (MFEM_USE_CUDA AND MFEM_EXPORT_GPU_CONFIG)
-    if (MFEM_CUDA_COMPILER_IS_NVCC)
-      set(MFEM_XLINKER "-Xlinker=")
-    else()
-      set(MFEM_XLINKER "${CMAKE_CUDA_LINKER_WRAPPER_FLAG}")
-    endif()
-  else()
-    set(MFEM_XLINKER "${CMAKE_CXX_LINKER_WRAPPER_FLAG}")
-  endif()
+  # TODO: CUDA/HIP support:
+  set(MFEM_XLINKER "${CMAKE_CXX_LINKER_WRAPPER_FLAG}")
  set(MFEM_MPIEXEC ${MPIEXEC})
  if (NOT MFEM_MPIEXEC)
    set(MFEM_MPIEXEC "mpirun")
@@ -1054,21 +987,16 @@ function(mfem_export_mk_files)
    # handle interfaces (e.g., SCOREC::apf)
    if ("${lib}" MATCHES "SCOREC::.*" OR "${lib}" MATCHES "Ginkgo::.*" OR "${lib}" MATCHES "ParMoonolith::.*")
    elseif (TARGET "${lib}")
-      mfem_get_target_options(${lib} CompileOpts2 LinkOpts2)
-      # remove generator expressions
-      string(GENEX_STRIP "${CompileOpts2}" CompileOpts)
-      string(GENEX_STRIP "${LinkOpts2}" LinkOpts)
+      mfem_get_target_options(${lib} CompileOpts LinkOpts)
      # Removing duplicates may lead to issues:
      # list(REMOVE_DUPLICATES CompileOpts)
      # list(REMOVE_DUPLICATES LinkOpts)
-      # message(WARNING "${lib}[LinkOpts]: ${LinkOpts}")      
-      # message(WARNING "${lib}[CompileOpts]: ${CompileOpts}")
-      foreach(LOpt IN LISTS LinkOpts)
-        set(MFEM_EXT_LIBS "${MFEM_EXT_LIBS} ${LOpt}")
-      endforeach()
-      foreach(COpt IN LISTS CompileOpts)
-        set(MFEM_TPLFLAGS "${MFEM_TPLFLAGS} ${COpt}")
-      endforeach()
+      string(REPLACE ";" " " COpts "${CompileOpts}")
+      string(REPLACE ";" " " LOpts "${LinkOpts}")
+      # message(STATUS "${lib}[COpts]: '${COpts}'")
+      # message(STATUS "${lib}[LOpts]: '${LOpts}'")
+      set(MFEM_TPLFLAGS "${MFEM_TPLFLAGS} ${COpts}")
+      set(MFEM_EXT_LIBS "${MFEM_EXT_LIBS} ${LOpts}")
      # message(FATAL_ERROR "***** interface lib found ... exiting *****")
      # handle static and shared libs
    elseif ("${suffix}" STREQUAL "${CMAKE_SHARED_LIBRARY_SUFFIX}")
@@ -1076,7 +1004,7 @@ function(mfem_export_mk_files)
      get_filename_component(fullLibName ${lib} NAME_WE)
      string(REGEX REPLACE "^lib" "" libname ${fullLibName})
      set(MFEM_EXT_LIBS
-        "${MFEM_EXT_LIBS} ${shared_link_flag}${dir} -L${dir} -l${libname}")
+          "${MFEM_EXT_LIBS} ${shared_link_flag}${dir} -L${dir} -l${libname}")
    else()
      set(MFEM_EXT_LIBS "${MFEM_EXT_LIBS} ${lib}")
    endif()
@@ -1085,7 +1013,7 @@ function(mfem_export_mk_files)
  # Create the build-tree version of 'config.mk'
  configure_file(
    "${PROJECT_SOURCE_DIR}/config/config.mk.in"
-    "${PROJECT_BINARY_DIR}/config/config.mk" @ONLY)
+    "${PROJECT_BINARY_DIR}/config/config.mk")
  # Copy 'test.mk' from the source-tree to the build-tree
  configure_file(
    "${PROJECT_SOURCE_DIR}/config/test.mk"
@@ -1103,7 +1031,7 @@ function(mfem_export_mk_files)
  # Create the install-tree version of 'config.mk'
  configure_file(
    "${PROJECT_SOURCE_DIR}/config/config.mk.in"
-    "${PROJECT_BINARY_DIR}/config/config-install.mk" @ONLY)
+    "${PROJECT_BINARY_DIR}/config/config-install.mk")

  # Install rules for 'config.mk' and 'test.mk'
  install(FILES ${PROJECT_SOURCE_DIR}/config/test.mk
@@ -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

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

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

@@ -153,7 +152,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
@@ -369,19 +367,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))
@@ -422,7 +407,7 @@ AMGX_LIB = -L$(AMGX_DIR)/lib -lamgx -lcusparse -lcusolver -lcublas -lnvToolsExt
 # MAGMA library configuration
 MAGMA_DIR = @MFEM_DIR@/../magma
 MAGMA_OPT = -I$(MAGMA_DIR)/include
-MAGMA_LIB = -L$(MAGMA_DIR)/lib -l:libmagma.a $(LAPACK_LIB)
+MAGMA_LIB = -L$(MAGMA_DIR)/lib -l:libmagma.a -lcublas -lcusparse $(LAPACK_LIB)

 # GnuTLS library configuration
 GNUTLS_OPT =
@@ -635,7 +620,7 @@ 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\
+TRIBOL_LIB = -L$(TRIBOL_DIR)/lib -ltribol -lredecomp -L$(AXOM_DIR)/lib -laxom_mint\
   -laxom_slam -laxom_slic -laxom_core

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

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

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

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

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

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

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

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

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

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

 void ConductionOperator::SetParameters(const Vector &u)
@@ -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;
 }
@@ -119,7 +119,7 @@ int main(int argc, char *argv[])
   }

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

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

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

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

   // 11. Set up the parallel sesquilinear form a(.,.) on the finite element
@@ -540,21 +534,15 @@ int main(int argc, char *argv[])
      char vishost[] = "localhost";
      int  visport   = 19916;
      socketstream sol_sock_r(vishost, visport);
+      socketstream sol_sock_i(vishost, visport);
      sol_sock_r << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
      sol_sock_r.precision(8);
+      sol_sock_i.precision(8);
      sol_sock_r << "solution\n" << *pmesh << u.real()
                 << "window_title 'Solution: Real Part'" << flush;
-      // Make sure all ranks have sent their real solution before initiating
-      // another set of GLVis connections (one from each rank):
-      MPI_Barrier(pmesh->GetComm());
-      socketstream sol_sock_i(vishost, visport);
-      sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
-      sol_sock_i.precision(8);
      sol_sock_i << "solution\n" << *pmesh << u.imag()
                 << "window_title 'Solution: Imaginary Part'" << flush;
-      // Make sure all ranks have sent their imaginary solution before initiating
-      // another set of GLVis connections (one from each rank):
-      MPI_Barrier(pmesh->GetComm());
   }
   if (visualization && exact_sol)
   {
@@ -563,21 +551,15 @@ int main(int argc, char *argv[])
      char vishost[] = "localhost";
      int  visport   = 19916;
      socketstream sol_sock_r(vishost, visport);
+      socketstream sol_sock_i(vishost, visport);
      sol_sock_r << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
      sol_sock_r.precision(8);
+      sol_sock_i.precision(8);
      sol_sock_r << "solution\n" << *pmesh << u_exact->real()
                 << "window_title 'Error: Real Part'" << flush;
-      // Make sure all ranks have sent their real solution before initiating
-      // another set of GLVis connections (one from each rank):
-      MPI_Barrier(pmesh->GetComm());
-      socketstream sol_sock_i(vishost, visport);
-      sol_sock_i << "parallel " << num_procs << " " << myid << "\n";
-      sol_sock_i.precision(8);
      sol_sock_i << "solution\n" << *pmesh << u_exact->imag()
                 << "window_title 'Error: Imaginary Part'" << flush;
-      // Make sure all ranks have sent their imaginary solution before initiating
-      // another set of GLVis connections (one from each rank):
-      MPI_Barrier(pmesh->GetComm());
   }
   if (visualization)
   {
@@ -139,7 +139,7 @@ void WaveOperator::ImplicitSolve(const real_t fac0, const real_t fac1,
   // for d2udt2
   if (!T)
   {
-      T = Add(1.0, Mmat, fac0, Kmat);
+      T = Add((real_t)1.0, Mmat, fac0, Kmat);
      T_solver.SetOperator(*T);
   }
   K->FullMult(u, z);
@@ -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);
@@ -56,6 +56,51 @@ void f_exact(const Vector &, Vector &);
 real_t freq = 1.0, kappa;
 int dim;

+void SolveSingle(SparseMatrix &A, const Vector &B, Vector &X)
+{
+   VectorMP<float> Bs, Xs;
+   const real_t *data = A.GetData();
+
+   const int n = A.GetI()[A.NumRows()];
+   int *Icopy = new int[A.NumRows() + 1];
+   int *Jcopy = new int[n];
+
+   float *sdata = new float[n];
+   for (int i=0; i<n; ++i)
+   {
+      sdata[i] = data[i];
+      Jcopy[i] = A.GetJ()[i];
+   }
+
+   for (int i=0; i<A.NumRows() + 1; ++i)
+   {
+      Icopy[i] = A.GetI()[i];
+   }
+
+   SparseMatrixMP<float> As(Icopy, Jcopy, sdata, A.NumRows(), A.NumCols());
+
+   Bs.SetSize(B.Size());
+   Xs.SetSize(X.Size());
+
+   for (int i=0; i<B.Size(); ++i)
+   {
+      Bs[i] = B[i];
+   }
+
+   for (int i=0; i<X.Size(); ++i)
+   {
+      Xs[i] = X[i];
+   }
+
+   GSSmootherMP<float> Ms(As);
+   PCG<float>(As, Ms, Bs, Xs, 1, 500, 1e-12, 0.0);
+
+   for (int i=0; i<X.Size(); ++i)
+   {
+      X[i] = Xs[i];
+   }
+}
+
 int main(int argc, char *argv[])
 {
   // 1. Parse command-line options.
@@ -185,6 +230,7 @@ int main(int argc, char *argv[])

   cout << "Size of linear system: " << A->Height() << endl;

+   /*
   // 11. Solve the linear system A X = B.
   if (pa) // Jacobi preconditioning in partial assembly mode
   {
@@ -193,20 +239,23 @@ int main(int argc, char *argv[])
   }
   else
   {
-#ifndef MFEM_USE_SUITESPARSE
+   #ifndef MFEM_USE_SUITESPARSE
      // 11. Define a simple symmetric Gauss-Seidel preconditioner and use it to
      //     solve the system Ax=b with PCG.
      GSSmoother M((SparseMatrix&)(*A));
      PCG(*A, M, B, X, 1, 500, 1e-12, 0.0);
-#else
+   #else
      // 11. 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
+   #endif
   }
+   */
+
+   SolveSingle((SparseMatrix&)(*A), B, X);

   // 12. Recover the solution as a finite element grid function.
   a->RecoverFEMSolution(X, *b, x);
@@ -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);
@@ -297,6 +297,18 @@ int main(int argc, char *argv[])
      sol_sock << "solution\n" << *pmesh << x << flush;
   }

+   VectorMP<float> xf(x.Size());
+   for (int i=0; i<x.Size(); ++i)
+   {
+      xf[i] = x[i];
+   }
+
+   if (myid == 0)
+   {
+      cout << "Norm of x " << x.Norml2() << endl;
+      cout << "Norm of xf " << xf.Norml2() << endl;
+   }
+
   // 18. Free the used memory.
   delete a;
   delete sigma;
@@ -9,7 +9,6 @@
 //               ex4 -m ../data/beam-hex.mesh -o 2 -pa
 //               ex4 -m ../data/escher.mesh
 //               ex4 -m ../data/fichera.mesh -o 2 -hb
-//               ex4 -m ../data/fichera.mesh -o 2 -hb -ea
 //               ex4 -m ../data/fichera-q2.vtk
 //               ex4 -m ../data/fichera-q3.mesh -o 2 -sc
 //               ex4 -m ../data/square-disc-nurbs.mesh
@@ -19,7 +18,6 @@
 //               ex4 -m ../data/amr-quad.mesh
 //               ex4 -m ../data/amr-hex.mesh
 //               ex4 -m ../data/amr-hex.mesh -o 2 -hb
-//               ex4 -m ../data/amr-hex.mesh -o 2 -hb -ea
 //               ex4 -m ../data/fichera-amr.mesh -o 2 -sc
 //               ex4 -m ../data/ref-prism.mesh -o 1
 //               ex4 -m ../data/octahedron.mesh -o 1
@@ -27,8 +25,6 @@
 //
 // Device sample runs:
 //               ex4 -m ../data/star.mesh -pa -d cuda
-//               ex4 -m ../data/star.mesh -hb -ea -d cuda
-//               ex4 -m ../data/amr-quad.mesh -hb -ea -d cuda
 //               ex4 -m ../data/star.mesh -pa -d raja-cuda
 //               ex4 -m ../data/star.mesh -pa -d raja-omp
 //               ex4 -m ../data/beam-hex.mesh -pa -d cuda
@@ -197,7 +193,7 @@ int main(int argc, char *argv[])
   cout << "Size of linear system: " << A->Height() << endl;

   // 11. Solve the linear system A X = B.
-   if (!pa && (!ea || hybridization))
+   if (!pa)
   {
 #ifndef MFEM_USE_SUITESPARSE
      // Use a simple symmetric Gauss-Seidel preconditioner with PCG.
@@ -9,7 +9,6 @@
 //               mpirun -np 4 ex4p -m ../data/beam-hex.mesh -o 2 -pa
 //               mpirun -np 4 ex4p -m ../data/escher.mesh -o 2 -sc
 //               mpirun -np 4 ex4p -m ../data/fichera.mesh -o 2 -hb
-//               mpirun -np 4 ex4p -m ../data/fichera.mesh -o 2 -hb -ea
 //               mpirun -np 4 ex4p -m ../data/fichera-q2.vtk
 //               mpirun -np 4 ex4p -m ../data/fichera-q3.mesh -o 2 -sc
 //               mpirun -np 4 ex4p -m ../data/square-disc-nurbs.mesh -o 3
@@ -18,18 +17,14 @@
 //               mpirun -np 4 ex4p -m ../data/periodic-cube.mesh -no-bc
 //               mpirun -np 4 ex4p -m ../data/amr-quad.mesh
 //               mpirun -np 3 ex4p -m ../data/amr-quad.mesh -o 2 -hb
-//               mpirun -np 3 ex4p -m ../data/amr-quad.mesh -o 2 -hb -ea
 //               mpirun -np 4 ex4p -m ../data/amr-hex.mesh -o 2 -sc
 //               mpirun -np 4 ex4p -m ../data/amr-hex.mesh -o 2 -hb
-//               mpirun -np 4 ex4p -m ../data/amr-hex.mesh -o 2 -hb -ea
 //               mpirun -np 4 ex4p -m ../data/ref-prism.mesh -o 1
 //               mpirun -np 4 ex4p -m ../data/octahedron.mesh -o 1
 //               mpirun -np 4 ex4p -m ../data/star-surf.mesh -o 3 -hb
 //
 // Device sample runs:
 //               mpirun -np 4 ex4p -m ../data/star.mesh -pa -d cuda
-//               mpirun -np 4 ex4p -m ../data/star.mesh -ea -hb -d cuda
-//               mpirun -np 4 ex4p -m ../data/amr-hex.mesh -ea -hb -d cuda
 //               mpirun -np 4 ex4p -m ../data/star.mesh -pa -d raja-cuda
 //               mpirun -np 4 ex4p -m ../data/star.mesh -pa -d raja-omp
 //               mpirun -np 4 ex4p -m ../data/beam-hex.mesh -pa -d cuda
@@ -235,7 +230,7 @@ int main(int argc, char *argv[])
   pcg->SetMaxIter(2000);
   pcg->SetPrintLevel(1);
   if (hybridization) { prec = new HypreBoomerAMG(*A.As<HypreParMatrix>()); }
-   else if (pa || ea) { prec = new OperatorJacobiSmoother(*a, ess_tdof_list); }
+   else if (pa) { prec = new OperatorJacobiSmoother(*a, ess_tdof_list); }
   else
   {
      ParFiniteElementSpace *prec_fespace =
@@ -160,7 +160,6 @@ int main(int argc, char *argv[])
   bool paraview = false;
   bool binary = false;
   int vis_steps = 5;
-   bool solve_implicit_state = false;

   int precision = 8;
   cout.precision(precision);
@@ -188,9 +187,6 @@ int main(int argc, char *argv[])
                  "Final time; start time is 0.");
   args.AddOption(&dt, "-dt", "--time-step",
                  "Time step.");
-   args.AddOption(&solve_implicit_state, "-imp-state", "--implicit-state",
-                  "-imp-slope", "--implicit-slope",
-                  "Implicitly solve for stage state or slope.");
   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                  "--no-visualization",
                  "Enable or disable GLVis visualization.");
@@ -370,11 +366,6 @@ int main(int argc, char *argv[])
   //    right-hand side, and perform time-integration (looping over the time
   //    iterations, ti, with a time-step dt).
   FE_Evolution adv(m, k, b);
-   using ImplicitVariableType = FE_Evolution::ImplicitVariableType;
-   ImplicitVariableType imp_var = solve_implicit_state ?
-                                  ImplicitVariableType::STATE
-                                  : ImplicitVariableType::SLOPE;
-   adv.SetImplicitVariableType(imp_var);

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

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

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

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

+# Currently there are only serial Ginkgo examples
 SEQ_EXAMPLES = ex1
-PAR_EXAMPLES = ex1p
+PAR_EXAMPLES =
 ifeq ($(MFEM_USE_MPI),NO)
   EXAMPLES = $(SEQ_EXAMPLES)
 else
@@ -76,4 +77,4 @@ clean-build:
 	rm -rf *.dSYM *.TVD.*breakpoints

 clean-exec:
-	@rm -f refined.mesh sol.gf mesh.* sol.*
+	@rm -f refined.mesh sol.gf
@@ -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
@@ -133,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
@@ -171,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
@@ -195,14 +191,12 @@ 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_vecmass_pa.hpp
  coefficient.hpp
@@ -311,7 +305,6 @@ set(HDRS
  tmop_tools.hpp
  tmop_amr.hpp
  gslib.hpp
-  gslib/gslib_kernel_helpers.hpp
  transfer.hpp
  hyperbolic.hpp
  integrator.hpp
@@ -729,8 +729,7 @@ void BilinearForm::Assemble(int skip_zeros)
         tr = mesh -> GetBdrFaceTransformations (i);
         if (tr != NULL)
         {
-            mfem::DofTransformation doftrans;
-            fes -> GetElementVDofs (tr -> Elem1No, vdofs, doftrans);
+            fes -> GetElementVDofs (tr -> Elem1No, vdofs);
            fe1 = fes -> GetFE (tr -> Elem1No);
            // The fe2 object is really a dummy and not used on the boundaries,
            // but we can't dereference a NULL pointer, and we don't want to
@@ -744,7 +743,6 @@ void BilinearForm::Assemble(int skip_zeros)

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

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

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


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

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

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

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

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

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

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

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

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

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

    Finally, only tensor-product elements are currently supported.

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

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

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

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

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

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

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

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

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

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

   void Setup(const int nb_i, const int ncp_i, const int b_type_i,
@@ -33,7 +33,6 @@ class FiniteElement;
 class FiniteElementSpace;
 class ElementTransformation;
 class IntegrationRule;
-class Vector;

 /** @brief Function that determines if a CEED kernel should be used, based on
    the current mfem::Device configuration. */
@@ -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);
@@ -1308,73 +1302,6 @@ real_t TraceCoefficient::Eval(ElementTransformation &T,
   return ma.Trace();
 }

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

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

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

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

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

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

-real_t
-ComplexGridFunction::ComputeLpError(const real_t p,
-                                    Coefficient &exsolr,
-                                    Coefficient &exsoli,
-                                    Coefficient *weight,
-                                    const IntegrationRule *irs[],
-                                    const Array<int> *elems) const
-{
-   real_t error = 0.0;
-   const FiniteElement *fe;
-   ElementTransformation *T;
-   Vector valsr;
-   Vector valsi;
-
-   const GridFunction& gf_r = real();
-   const GridFunction& gf_i = imag();
-
-   for (int i = 0; i < fes->GetNE(); i++)
-   {
-      if (elems != NULL && (*elems)[i] == 0) { continue; }
-      fe = fes->GetFE(i);
-      const IntegrationRule *ir;
-      if (irs)
-      {
-         ir = irs[fe->GetGeomType()];
-      }
-      else
-      {
-         int intorder = 2*fe->GetOrder() + 3;
-         ir = &(IntRules.Get(fe->GetGeomType(), intorder));
-      }
-      real_t elem_error = 0.0;
-      gf_r.GetValues(i, *ir, valsr);
-      gf_i.GetValues(i, *ir, valsi);
-      T = fes->GetElementTransformation(i);
-      for (int j = 0; j < ir->GetNPoints(); j++)
-      {
-         const IntegrationPoint &ip = ir->IntPoint(j);
-         T->SetIntPoint(&ip);
-         real_t diffr = valsr(j) - exsolr.Eval(*T, ip);
-         real_t diffi = valsi(j) - exsoli.Eval(*T, ip);
-         real_t diff = hypot(diffr, diffi);
-         if (p < infinity())
-         {
-            diff = pow(diff, p);
-            if (weight)
-            {
-               diff *= weight->Eval(*T, ip);
-            }
-            elem_error += ip.weight * T->Weight() * diff;
-         }
-         else
-         {
-            if (weight)
-            {
-               diff *= weight->Eval(*T, ip);
-            }
-            error = std::max(error, diff);
-         }
-      }
-      if (p < infinity())
-      {
-         // negative quadrature weights may cause the error to be negative
-         error += fabs(elem_error);
-      }
-   }
-
-   if (p < infinity())
-   {
-      error = pow(error, 1./p);
-   }
-
-   return error;
-}
-
 void ComplexGridFunction::Save(std::ostream &os) const
 {
   os << "ComplexGridFunction\n";
@@ -905,9 +830,15 @@ ParComplexGridFunction::ParComplexGridFunction(ParMesh *m, std::istream &input)
      int vsize = pfes->GetVSize();
      Vector::Load(input, 2*vsize);

-      real_t *h_data = HostReadWrite();
-      pfes->ApplyDofSigns(h_data);
-      pfes->ApplyDofSigns(h_data + vsize);
+      real_t *data_  = const_cast<real_t*>(HostRead());
+      for (int i = 0; i < vsize; i++)
+      {
+         if (pfes->GetDofSign(i) < 0)
+         {
+            data_[i] = -data_[i];
+            data_[i+vsize] = -data_[i+vsize];
+         }
+      }


      // if the mesh is a legacy (v1.1) NC mesh, it has old vertex ordering
@@ -1120,14 +1051,15 @@ void ParComplexGridFunction::Save(std::ostream &os) const
   os << '\n';

   int vsize = pfes->GetVSize();
-   // We use const_cast + HostRead (instead of HostReadWrite) because we only
-   // need to change the host data temporarily and this way we do not invalidate
-   // the data if it is on device. If we use HostReadWrite here, later calls to
-   // Read or ReadWrite will need to copy the data from host to device. With the
-   // approach used here, the host-to-device copy is avoided.
-   real_t *h_data = const_cast<real_t*>(HostRead());
-   pfes->ApplyDofSigns(h_data);
-   pfes->ApplyDofSigns(h_data + vsize);
+   real_t *data_  = const_cast<real_t*>(HostRead());
+   for (int i = 0; i < vsize; i++)
+   {
+      if (pfes->GetDofSign(i) < 0)
+      {
+         data_[i] = -data_[i];
+         data_[i+vsize] = -data_[i+vsize];
+      }
+   }

   if (pfes->GetOrdering() == Ordering::byNODES)
   {
@@ -1138,8 +1070,14 @@ void ParComplexGridFunction::Save(std::ostream &os) const
      Vector::Print(os, pfes->GetVDim());
   }

-   pfes->ApplyDofSigns(h_data);
-   pfes->ApplyDofSigns(h_data + vsize);
+   for (int i = 0; i < vsize; i++)
+   {
+      if (pfes->GetDofSign(i) < 0)
+      {
+         data_[i] = -data_[i];
+         data_[i+vsize] = -data_[i+vsize];
+      }
+   }

   os.flush();
 }
@@ -82,25 +82,6 @@ public:
   /// underlying #fes
   int VectorDim() const;

-   /// Copy assignment. Only the data of the base class Vector is copied.
-   /** It is assumed that this object and @a rhs use FiniteElementSpace%s that
-       have the same size.
-
-       @note Defining this method overwrites the implicitly defined copy
-       assignment operator. */
-   ComplexGridFunction &operator=(const ComplexGridFunction &rhs)
-   { return operator=((const Vector &)rhs); }
-
-   /// Copy the data from @a v.
-   /** The size of @a v must be equal to double of the size of the associated
-       FiniteElementSpace #fes. */
-   ComplexGridFunction &operator=(const Vector &v)
-   {
-      MFEM_ASSERT(fes && v.Size() == 2*fes->GetVSize(), "");
-      Vector::operator=(v);
-      return *this;
-   }
-
   /// Assign constant values to the ComplexGridFunction data.
   ComplexGridFunction &operator=(const std::complex<real_t> & value)
   { *gfr = value.real(); *gfi = value.imag(); return *this; }
@@ -166,75 +147,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;

@@ -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,
@@ -492,8 +492,6 @@ void VisItDataCollection::SaveRootFile()
                           to_padded_string(cycle, pad_digits_cycle) +
                           ".mfem_root";
   std::ofstream root_file(root_name);
-   MFEM_VERIFY(root_file.is_open(),
-               "Failed to open ofstream " << root_name);
   root_file << GetVisItRootString();
   if (!root_file)
   {
@@ -979,10 +977,7 @@ void ParaViewDataCollection::Save()
   // Save the local part of the mesh and grid functions fields to the local
   // VTU file. Also save coefficient fields.
   {
-      std::string os_str = vtu_prefix + GenerateVTUFileName("proc", myid);
-      std::ofstream os(os_str);
-      MFEM_VERIFY(os.is_open(),
-                  "Failed to open ofstream " << os_str);
+      std::ofstream os(vtu_prefix + GenerateVTUFileName("proc", myid));
      os.precision(precision);
      SaveDataVTU(os, levels_of_detail);
   }
@@ -994,10 +989,7 @@ void ParaViewDataCollection::Save()
                  "QuadratureFunction output is not supported for "
                  "ParaViewDataCollection on domain boundary!");
      const std::string &field_name = qfield.first;
-      std::string os_str = vtu_prefix + GenerateVTUFileName(field_name, myid);
-      std::ofstream os(os_str);
-      MFEM_VERIFY(os.is_open(),
-                  "Failed to open ofstream " << os_str);
+      std::ofstream os(vtu_prefix + GenerateVTUFileName(field_name, myid));
      qfield.second->SaveVTU(os, pv_data_format, GetCompressionLevel(), field_name);
   }

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

         // Grid function fields and coefficient fields
@@ -1066,10 +1055,8 @@ void ParaViewDataCollection::Save()
         const std::string &q_field_name = q_field.first;
         std::string q_fname = GeneratePVTUPath() + "/"
                               + GeneratePVTUFileName(q_field_name);
-         std::string os_str = col_path + "/" + q_fname;
-         std::ofstream pvtu_out(os_str);
-         MFEM_VERIFY(pvtu_out.is_open(),
-                     "Failed to open ofstream " << os_str);
+
+         std::ofstream pvtu_out(col_path + "/" + q_fname);
         WritePVTUHeader(pvtu_out);
         int vec_dim = q_field.second->GetVDim();
         pvtu_out << "<PPointData>\n";
@@ -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++)
      {
@@ -25,35 +25,21 @@
 namespace mfem
 {

-/// Lightweight adaptor over an std::map from type K to type to V
-template<typename K, typename V,
-         typename = typename std::enable_if<std::is_default_constructible<V>::value>::type>
-class GenericFieldMap
+/// Lightweight adaptor over an std::map from strings to pointer to T
+template<typename T>
+class NamedFieldsMap
 {
-private:
-   static constexpr bool ValueIsPointer = std::is_pointer<V>::value;
-
 public:
-   typedef std::map<K, V> MapType;
+   typedef std::map<std::string, T*> MapType;
   typedef typename MapType::iterator iterator;
   typedef typename MapType::const_iterator const_iterator;

-   /// Register field @a field with name @a key
-   /// Only enabled if the template parameter V is not a pointer
-   template<typename = std::enable_if<!ValueIsPointer, bool>>
-   void Register(const K& key, V field)
+   /// Register field @a field with name @a fname
+   /** Replace existing field associated with @a fname (and optionally
+       delete associated pointer if @a own_data is true) */
+   void Register(const std::string& fname, T* field, bool own_data)
   {
-      field_map[key] = field;
-   }
-
-   /// Register field @a field with name @a key
-   /** Replace existing field associated with @a key (and optionally
-       delete associated pointer if @a own_data is true).
-       Only enabled if the template parameter V is a pointer*/
-   template<typename = std::enable_if<ValueIsPointer, bool>>
-   void Register(const K& key, V field, bool own_data)
-   {
-      V& ref = field_map[key];
+      T*& ref = field_map[fname];
      if (own_data)
      {
         delete ref; // if newly allocated -> ref is null -> OK
@@ -61,40 +47,23 @@ public:
      ref = field;
   }

-   /// Unregister association between field @a field and name @a key
-   /// Only enabled if the template parameter V is not a pointer
-   template<typename = std::enable_if<!ValueIsPointer, bool>>
-   void Deregister(const K& key)
+   /// Unregister association between field @a field and name @a fname
+   /** Optionally delete associated pointer if @a own_data is true */
+   void Deregister(const std::string& fname, bool own_data)
   {
-      iterator it = field_map.find(key);
-      if ( it != field_map.end() )
-      {
-         field_map.erase(it);
-      }
-   }
-
-   /// Unregister association between field @a field and name @a key
-   /** Optionally delete associated pointer if @a own_data is true.
-       Only enabled if the template parameter V is a pointer */
-   template<typename = std::enable_if<ValueIsPointer, bool>>
-   void Deregister(const K& key, bool own_data)
-   {
-      iterator it = field_map.find(key);
+      iterator it = field_map.find(fname);
      if ( it != field_map.end() )
      {
         if (own_data)
         {
            delete it->second;
-            it->second = nullptr;
         }
         field_map.erase(it);
      }
   }

   /// Clear all associations between names and fields
-   /** Delete associated pointers when @a own_data is true.
-       Only enabled if the template parameter V is a pointer */
-   template<typename = std::enable_if<ValueIsPointer, bool>>
+   /** Delete associated pointers when @a own_data is true */
   void DeleteData(bool own_data)
   {
      for (iterator it = field_map.begin(); it != field_map.end(); ++it)
@@ -107,37 +76,22 @@ public:
      }
   }

-   /// Predicate to check if a field is associated with name @a key
-   bool Has(const K& key) const
+   /// Predicate to check if a field is associated with name @a fname
+   bool Has(const std::string& fname) const
   {
-      return field_map.find(key) != field_map.end();
+      return field_map.find(fname) != field_map.end();
   }

-   /// Get a pointer to the field associated with name @a key
-   /** @return Field associated with @a key or NULL,
-       if value is pointer and key not found */
-   V Get(const K& key) const
+   /// Get a pointer to the field associated with name @a fname
+   /** @return Pointer to field associated with @a fname or NULL */
+   T* Get(const std::string& fname) const
   {
-      const_iterator it = field_map.find(key);
-      if (it != field_map.end())
-      {
-         return it->second;
-      }
-      else
-      {
-         if constexpr (ValueIsPointer)
-         {
-            return nullptr;
-         }
-         else
-         {
-            return V(); // Return default-constructed value for non-pointer types
-         }
-      }
+      const_iterator it = field_map.find(fname);
+      return it != field_map.end() ? it->second : NULL;
   }

   /// Returns a const reference to the underlying map
-   const MapType &GetMap() const { return field_map; }
+   const MapType& GetMap() const { return field_map; }

   /// Returns the number of registered fields
   int NumFields() const { return field_map.size(); }
@@ -152,24 +106,21 @@ public:
   /// Returns an end const iterator to the registered fields
   const_iterator end() const { return field_map.end(); }

-   /// Returns an iterator to the field @a key
-   iterator find(const K& key)
-   { return field_map.find(key); }
+   /// Returns an iterator to the field @a fname
+   iterator find(const std::string& fname)
+   { return field_map.find(fname); }

-   /// Returns a const iterator to the field @a key
-   const_iterator find(const K& key) const
-   { return field_map.find(key); }
+   /// Returns a const iterator to the field @a fname
+   const_iterator find(const std::string& fname) const
+   { return field_map.find(fname); }

-   /// Clears the map of registered fields
+   /// Clears the map of registered fields without reclaiming memory
   void clear() { field_map.clear(); }

 protected:
   MapType field_map;
 };

-/// Lightweight adaptor over an std::map from strings to pointer to T
-template<typename T>
-using NamedFieldsMap = GenericFieldMap<std::string, T*>;

 /** A class for collecting finite element data that is part of the same
    simulation. Currently, this class groups together grid functions (fields),
@@ -90,8 +90,8 @@ void map_quadrature_data_to_fields_impl(
   }
   else
   {
-      MFEM_ABORT_KERNEL("quadrature data mapping to field is not implemented"
-                        " for this field descriptor");
+      MFEM_ABORT("quadrature data mapping to field is not implemented for"
+                 " this field descriptor");
   }
 }

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

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

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

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

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

 #include <ostream>
 #include "../../config/config.hpp"
@@ -671,20 +671,6 @@ public:
                    MPI_COMM_WORLD);
   }

-   Operator& GetGradient(const Vector &x0) const override
-   {
-      x = x0;
-      f.UseDevice(x.UseDevice());
-      xpev.UseDevice(x.UseDevice());
-
-      op.Mult(x, f);
-      const real_t xnorm_local = x.Norml2();
-      MPI_Allreduce(&xnorm_local, &xnorm, 1, MPITypeMap<real_t>::mpi_type, MPI_SUM,
-                    MPI_COMM_WORLD);
-
-      return const_cast<FDJacobian&>(*this);
-   }
-
   void Mult(const Vector &v, Vector &y) const override
   {
      // See [1] for choice of eps.
@@ -739,11 +725,11 @@ public:

 private:
   const Operator &op;
-   mutable Vector x, f;
+   Vector x, f;
   mutable Vector xpev;
   real_t lambda = 1.0e-6;
   real_t fixed_eps;
-   mutable real_t xnorm;
+   real_t xnorm;
 };

 /// @brief Find the index of a field descriptor in a vector of field descriptors.
@@ -12,6 +12,7 @@
 #include "dgmassinv.hpp"
 #include "bilinearform.hpp"
 #include "dgmassinv_kernels.hpp"
+#include "../general/forall.hpp"

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

 DGMassInverse::~DGMassInverse() = default;

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

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

 namespace mfem
 {
@@ -57,7 +56,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
      {
@@ -334,170 +333,6 @@ void DGMassBasis(const int e,

 } // namespace internal

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

 #endif
@@ -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
@@ -69,9 +69,9 @@ inline int ToLexOrdering2D(const int face_id, const int size1d, const int i)
 }

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

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

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

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

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

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

   /// Returns absolute value of the maps
   DofToQuad Abs() const;
-
-   /// Auxiliary function for searching DofToQuad arrays.
-   static inline DofToQuad *SearchArray(
-      const Array<DofToQuad*> &dof2quad_array,
-      const IntegrationRule &ir,
-      DofToQuad::Mode mode);
-
-   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
@@ -336,20 +289,10 @@ public:
                          $ u(x) = (1/w) \hat u(\hat x) $ */
      H_DIV,     /**< For vector fields; preserves surface integrals of the
                          normal component $ u(x) = (J/w) \hat u(\hat x) $ */
-      H_CURL,    /**< For vector fields; preserves line integrals of the
+      H_CURL     /**< For vector fields; preserves line integrals of the
                          tangential component
                          $ u(x) = J^{-t} \hat u(\hat x) $ (square J),
                          $ u(x) = J(J^t J)^{-1} \hat u(\hat x) $ (general J) */
-      H_DIV_R2D, /**< For 3-component vector fields in 2D; equivalent to a
-                          direct sum of an H_DIV basis and an INTEGRAL basis */
-      H_CURL_R2D,/**< For 3-component vector fields in 2D; equivalent to a
-                          direct sum of an H_CURL basis and a VALUE basis */
-      H_DIV_R1D, /**< For 3-component vector fields in 1D; equivalent to a
-                          direct sum of a VALUE basis and a pair of INTEGRAL
-                          bases */
-      H_CURL_R1D /**< For 3-component vector fields in 1D; equivalent to a
-                          direct sum of an INTEGRAL basis and a pair of VALUE
-                          bases */
   };

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

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

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

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

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

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

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

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

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

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

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

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

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

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

 P0SegmentFiniteElement::P0SegmentFiniteElement(int Ord)
   : NodalFiniteElement(1, Geometry::SEGMENT, 1, Ord)   // default Ord = 0
--- a/Show More
+++ b/Show More
Author	SHA1	Message	Date
dylan-copeland	0cf995cfc7	Fix single precision build.	2025-12-12 18:29:30 -08:00
dylan-copeland	c9535ee6fc	Templates for some reducers.	2025-12-12 17:09:03 -08:00
dylan-copeland	03766e06be	More templates.	2025-12-12 13:57:12 -08:00
dylan-copeland	9a96ed3cb1	Adding some missing templates.	2025-12-12 13:25:10 -08:00
Dylan Copeland	16afb15588	Merge branch 'master' of github.com:mfem/mfem into mp	2025-12-12 12:35:54 -08:00
Dylan Copeland	e480d543eb	Merge branch 'master' of github.com:mfem/mfem into mp	2025-12-12 12:34:20 -08:00
Dylan Copeland	cc762f8f39	More CI build fixes.	2025-01-31 19:41:13 -08:00
dylan-copeland	e9980453cc	Fixes for single precision build.	2025-01-31 17:04:49 -08:00
dylan-copeland	816ac475e5	Fixes for -Wshadow.	2025-01-31 16:15:28 -08:00
Dylan Copeland	b4fb21d3fa	Fixes for clang.	2025-01-31 16:03:07 -08:00
Dylan Copeland	7c7dc1469e	Fixes for another compiler.	2025-01-31 14:14:40 -08:00
dylan-copeland	b7903bfc6b	Fixed mac build.	2025-01-31 10:54:16 -08:00
Dylan Copeland	ad2483d36e	Templates for data types in serial matrices.	2025-01-31 10:31:38 -08:00
Dylan Copeland	738100e50a	Template for vector data type.	2025-01-30 12:21:47 -08:00