Fix a small issue in miniapps/shifted/distance.cpp

In GridFunction::ComputeLpError, use the device implementation only
in supported cases.
2025-08-18 19:06:08 -07:00 · 2025-08-18 16:31:56 -07:00 · 2025-08-18 13:52:10 -07:00 · 2025-08-18 13:39:23 -07:00 · 2025-08-18 13:37:17 -07:00 · 2025-08-16 16:01:21 -07:00
98 changed files with 3412 additions and 586 deletions
@@ -289,10 +289,10 @@ jobs:
        run: |
          export HOMEBREW_NO_INSTALL_CLEANUP=1
          brew update
-          brew install llvm@19 enzyme
-          echo "LLVM_PREFIX=$(brew --prefix llvm@19)" >> $GITHUB_ENV
-          echo "OMPI_CC=$(brew --prefix llvm@19)/bin/clang" >> $GITHUB_ENV
-          echo "OMPI_CXX=$(brew --prefix llvm@19)/bin/clang++" >> $GITHUB_ENV
+          brew install llvm@20 enzyme
+          echo "LLVM_PREFIX=$(brew --prefix llvm@20)" >> $GITHUB_ENV
+          echo "OMPI_CC=$(brew --prefix llvm@20)/bin/clang" >> $GITHUB_ENV
+          echo "OMPI_CXX=$(brew --prefix llvm@20)/bin/clang++" >> $GITHUB_ENV

      # MFEM build and test
      - name: build
@@ -211,7 +211,7 @@ miniapps/electromagnetics/joule
 miniapps/electromagnetics/Volta-AMR*
 miniapps/electromagnetics/Tesla-AMR*
 miniapps/electromagnetics/Maxwell-Parallel*
-miniapps/electromagnetics/Joule_*
+miniapps/electromagnetics/Joule_[0-9]*

 miniapps/gslib/field-diff
 miniapps/gslib/field-interp
@@ -267,9 +267,9 @@ miniapps/meshing/bounding-box*
 miniapps/meshing/jacobian-determinant*

 miniapps/mtop/parheat
-miniapps/mtop/ParHeat*
+miniapps/mtop/ParHeat/*
 miniapps/mtop/seqheat
-miniapps/mtop/SeqHeat*
+miniapps/mtop/SeqHeat/*

 miniapps/autodiff/paradiff
 miniapps/autodiff/seqadiff
@@ -277,7 +277,7 @@ miniapps/autodiff/seqtest
 miniapps/autodiff/par_example
 miniapps/autodiff/seq_example
 miniapps/autodiff/seq_test
-miniapps/autodiff/Exampl*
+miniapps/autodiff/Example/*

 miniapps/navier/navier_mms
 miniapps/navier/navier_kovasznay
@@ -300,6 +300,7 @@ miniapps/nurbs/nurbs_solenoidal
 miniapps/nurbs/nurbs_printfunc
 miniapps/nurbs/nurbs_patch_ex1
 miniapps/nurbs/nurbs_curveint
+miniapps/nurbs/nurbs_surface
 miniapps/nurbs/refined.mesh
 miniapps/nurbs/mesh.*
 miniapps/nurbs/sol_?.gf
@@ -318,6 +319,7 @@ miniapps/nurbs/nurbs_naca_cmesh
 miniapps/nurbs/naca-cmesh.mesh
 miniapps/nurbs/glvis_naca-cmesh.mesh
 miniapps/nurbs/Naca_cmesh
+miniapps/nurbs/*-Surface.mesh

 miniapps/performance/ex1
 miniapps/performance/ex1p
@@ -22,7 +22,7 @@ include:
 # the "needs" keyword and express the DAG of jobs for more efficiency.
 # - We use setup and setup_baseline phases to download content outside of mfem
 #   directory.
-# - Allocate/Release is where ruby resource are allocated/released once for all.
+# - Allocate/Release is where Dane resource are allocated/released once for all.
 # - Build and Test is where we build and MFEM for multiple toolchains.
 # - Baseline_checks gathers baseline-type test suites execution
 # - Baseline_publish, only available on master, allows to update baseline
@@ -53,7 +53,7 @@ variables:
  AUTOTEST_COMMIT: "YES"

 # Trigger subpipelines:
-ruby-build-and-test:
+dane-build-and-test:
  stage: sub-pipelines
  variables:
    # Explicitly pass down values that we want to be able to set when triggering
@@ -61,10 +61,10 @@ ruby-build-and-test:
    AUTOTEST: "${AUTOTEST}"
    AUTOTEST_COMMIT: "${AUTOTEST_COMMIT}"
  trigger:
-    include: .gitlab/ruby-build-and-test.yml
+    include: .gitlab/dane-build-and-test.yml
    strategy: depend

-ruby-baseline:
+dane-baseline:
  stage: sub-pipelines
  variables:
    # Explicitly pass down values that we want to be able to set when triggering
@@ -73,7 +73,7 @@ ruby-baseline:
    AUTOTEST: "${AUTOTEST}"
    AUTOTEST_COMMIT: "${AUTOTEST_COMMIT}"
  trigger:
-    include: .gitlab/ruby-baseline.yml
+    include: .gitlab/dane-baseline.yml
    strategy: depend

 lassen-build-and-test:
@@ -24,7 +24,7 @@ and `test type`.

 Machines typically include:

-* Ruby: 2nd Gen Intel Xeon (Cascade Lake)
+* Dane: Intel Sapphire Rapids
 * Lassen: Power9 + Nvidia GPU
 * Corona: AMD GPU

@@ -76,13 +76,13 @@ with a spack spec of MFEM, within the limits permitted by the MFEM spack
 package.

 In any build-and-test sub-pipeline a job basically consists in defining the
-spack spec to use. Adding a job on ruby for example resumes to:
+spack spec to use. Adding a job on Dane for example resumes to:

 ```yaml
 <job_name>:
  variables:
    SPEC: "<spack_spec>"
-  extends: .build_and_test_on_ruby
+  extends: .build_and_test_on_dane
 ```

 The remaining and non trivial work is to make sure this spec is working. To
@@ -24,7 +24,7 @@ variables:
 # TODO: add a clean-up mechanism
  BUILD_ROOT: ${USER_CI_TOP_DIR}/${CI_PROJECT_NAME}-${MACHINE_NAME}-pipeline-${CI_PIPELINE_ID}

-# On LLNL's ruby, there is only one allocation shared among jobs in order to
+# On LLNL's Dane, there is only one allocation shared among jobs in order to
 # save time and resource. This allocation has to be uniquely named so that we
 # are sure to retrieve it.
  ALLOC_NAME: ${CI_PROJECT_NAME}_ci_${CI_PIPELINE_ID}
@@ -9,17 +9,17 @@
 # terms of the BSD-3 license. We welcome feedback and contributions, see file
 # CONTRIBUTING.md for details.

-# GitLab pipelines configurations for the Ruby machine at LLNL
+# GitLab pipelines configurations for the Dane machine at LLNL
 variables:
-  MACHINE_NAME: ruby
+  MACHINE_NAME: dane

-.on_ruby:
+.on_dane:
  tags:
    - shell
-    - ruby
+    - dane
  rules:
-    # Don't run ruby jobs if...
-    - if: '$CI_COMMIT_BRANCH =~ /_qnone/ || $ON_RUBY == "OFF"'
+    # Don't run dane jobs if...
+    - if: '$CI_COMMIT_BRANCH =~ /_qnone/ || $ON_DANE == "OFF"'
      when: never
    # Don't run autotest update if...
    - if: '$CI_JOB_NAME =~ /report/ && $AUTOTEST != "YES"'
@@ -40,16 +40,17 @@ variables:
    - when: on_success

 # Spack helped builds
-# Generic ruby build job, extending build script
-.build_and_test_on_ruby:
-  extends: [.on_ruby]
+# Generic dane build job, extending build script
+.build_and_test_on_dane:
+  extends: [.on_dane]
  stage: build_and_test
  script:
    # THREADS is used by 'tests/gitlab/build_and_test', run below
-    - export THREADS=16
+    # Dane has 224 threads/node and we run 7 separate jobs: 224=7*32
+    - export THREADS=28
    - echo ${ALLOC_NAME}
    - export JOBID=$(squeue -h --name=${ALLOC_NAME} --format=%A)
    - echo ${JOBID}
    - echo ${MFEM_DATA_DIR}
    - echo ${SPEC}
-    - srun $( [[ -n "${JOBID}" ]] && echo "--jobid=${JOBID}" ) --reservation=ci -t 45 -N 1 tests/gitlab/build_and_test --spec "${SPEC}" --data-dir "${MFEM_DATA_DIR}" --data
+    - srun $( [[ -n "${JOBID}" ]] && echo "--jobid=${JOBID}" ) --reservation=ci -t 60 -N 1 tests/gitlab/build_and_test --spec "${SPEC}" --data-dir "${MFEM_DATA_DIR}" --data
@@ -18,7 +18,7 @@
 setup_baseline:
  tags:
    - shell
-    - ruby
+    - dane
  stage: setup
  variables:
    GIT_STRATEGY: none
@@ -16,7 +16,7 @@
 setup:
  tags:
    - shell
-    - ruby
+    - dane
  stage: setup
  variables:
    GIT_STRATEGY: none
@@ -19,8 +19,8 @@ stages:
  - cleanup
  - baseline_publish

-baselinecheck_mfem_intel_ruby:
-  extends: [.on_ruby]
+baselinecheck_mfem_intel_dane:
+  extends: [.on_dane]
  stage: baseline_check
  variables:
    # TPLS_DIR is used in .gitlab/scripts/baseline to provide the tpls location
@@ -31,8 +31,8 @@ baselinecheck_mfem_intel_ruby:
  script:
    - echo ${BUILD_ROOT}
    - echo ${TPLS_DIR}
-    # Used by the tests in MFEM/tests:
-    - export MFEM_TEST_NP=48
+    # Used by the tests in MFEM/tests, dane has 224 threads/node:
+    - export MFEM_TEST_NP=192
    # The next script uses the following environment variables:
    # * BASELINE_TEST, SYS_TYPE, CI_PROJECT_DIR, ARTIFACTS_DIR,
    # * BUILD_ROOT, TPLS_DIR, MACHINE_NAME
@@ -44,7 +44,7 @@ baselinecheck_mfem_intel_ruby:
  allow_failure: true

 cleanup:
-  extends: .on_ruby
+  extends: .on_dane
  stage: cleanup
  variables:
    GIT_STRATEGY: none
@@ -53,7 +53,7 @@ cleanup:
    - rm -rf "${BUILD_ROOT}" || true

 report_baseline:
-  extends: [.on_ruby]
+  extends: [.on_dane]
  stage: baseline_report
  script:
    - echo ${MACHINE_NAME}
@@ -113,8 +113,8 @@ report_baseline:
        exit $err
      ) 9> autotest.lock

-baselinepublish_mfem_ruby:
-  extends: [.on_ruby]
+baselinepublish_mfem_dane:
+  extends: [.on_dane]
  stage: baseline_publish
  rules:
    # - if: '$CI_COMMIT_BRANCH == "master" || $REBASELINE == "YES"'
@@ -129,5 +129,5 @@ baselinepublish_mfem_ruby:

 include:
  - local: .gitlab/configs/common.yml
-  - local: .gitlab/configs/ruby-config.yml
+  - local: .gitlab/configs/dane-config.yml
  - local: .gitlab/configs/setup-baseline.yml
@@ -19,54 +19,54 @@ stages:
 allocate_resource:
  variables:
    GIT_STRATEGY: none
-  extends: .on_ruby
+  extends: .on_dane
  stage: allocate_resource
  script:
    - echo ${ALLOC_NAME}
    - salloc --exclusive --nodes=1 --reservation=ci --time=60 --no-shell --job-name=${ALLOC_NAME}
  timeout: 6h

-# GitLab jobs for the Ruby machine at LLNL
+# GitLab jobs for the Dane machine at LLNL
 debug_ser_gcc_10:
  variables:
    SPEC: "%gcc@10.3.1 +debug~mpi"
-  extends: .build_and_test_on_ruby
+  extends: .build_and_test_on_dane

 debug_par_gcc_10:
  variables:
    SPEC: "%gcc@10.3.1 +debug+mpi"
-  extends: .build_and_test_on_ruby
+  extends: .build_and_test_on_dane

 opt_ser_gcc_10:
  variables:
    SPEC: "%gcc@10.3.1 ~mpi"
-  extends: .build_and_test_on_ruby
+  extends: .build_and_test_on_dane

 opt_par_gcc_10:
  variables:
    SPEC: "%gcc@10.3.1"
-  extends: .build_and_test_on_ruby
+  extends: .build_and_test_on_dane

 opt_par_gcc_10_sundials:
  variables:
    SPEC: "%gcc@10.3.1 +sundials"
-  extends: .build_and_test_on_ruby
+  extends: .build_and_test_on_dane

 opt_par_gcc_10_petsc:
  variables:
    SPEC: "%gcc@10.3.1 +petsc ^petsc+mumps~superlu-dist"
-  extends: .build_and_test_on_ruby
+  extends: .build_and_test_on_dane

 opt_par_gcc_10_pumi:
  variables:
    SPEC: "%gcc@10.3.1 +pumi"
-  extends: .build_and_test_on_ruby
+  extends: .build_and_test_on_dane

 # Release
 release_resource:
  variables:
    GIT_STRATEGY: none
-  extends: .on_ruby
+  extends: .on_dane
  stage: release_resource_and_report
  script:
    - echo ${ALLOC_NAME}
@@ -78,17 +78,17 @@ release_resource:
 report_job_success:
  stage: release_resource_and_report
  extends:
-    - .on_ruby
+    - .on_dane
    - .report_job_success

 report_job_failure:
  stage: release_resource_and_report
  extends:
-    - .on_ruby
+    - .on_dane
    - .report_job_failure

 include:
  - local: .gitlab/configs/common.yml
-  - local: .gitlab/configs/ruby-config.yml
+  - local: .gitlab/configs/dane-config.yml
  - local: .gitlab/configs/setup-build-and-test.yml
  - local: .gitlab/configs/report-build-and-test.yml
@@ -14,7 +14,7 @@
 # locals
 glob_err=${BASELINE_TEST}.err
 base=${BASELINE_TEST}-${SYS_TYPE}
-if [[ "${MACHINE_NAME}" == "ruby" ]]; then
+if [[ "${MACHINE_NAME}" == "dane" ]]; then
  base="${BASELINE_TEST}-${MACHINE_NAME}"
 fi
 base_diff=${base}.diff
@@ -31,7 +31,7 @@ cd tests
 mkdir _${BASELINE_TEST} && cd _${BASELINE_TEST}

 # run
-if [[ "${MACHINE_NAME}" == "ruby" ]]; then
+if [[ "${MACHINE_NAME}" == "dane" ]]; then
  salloc --nodes=1 --exclusive --reservation=ci ../runtest ../../mfem "${BASELINE_TEST} ${TPLS_DIR}"
 elif [[ ${MACHINE_NAME} == "corona" ]]; then
  salloc --nodes=1 -t 60 -p pbatch ../runtest ../../mfem "${BASELINE_TEST} ${TPLS_DIR}"
@@ -11,7 +11,7 @@
 # terms of the BSD-3 license. We welcome feedback and contributions, see file
 # CONTRIBUTING.md for details.

-# There will be collision between corona and ruby baselines.
+# There will be collision between corona and dane baselines.
 # Once the corresponding files have been generated, we can switch to machine
 # specific ref.
 ARTIFACT_PATH=${CI_PROJECT_DIR}/${ARTIFACTS_DIR}/baseline-${SYS_TYPE}
@@ -21,7 +21,7 @@ PATCH_FILE=${ARTIFACT_PATH}.patch
 FULL_FILE=${ARTIFACT_PATH}.out
 DIFF_FILE=${ARTIFACT_PATH}.diff

-# There will be collision between corona and ruby baselines.
+# There will be collision between corona and dane baselines.
 # Once the corresponding files have been generated, we can switch to machine
 # specific ref.
 SAVED_NAME=baseline-${SYS_TYPE}.saved
@@ -29,9 +29,14 @@ Discretization improvements

 Meshing improvements
 --------------------
+
 - Added support for higher order meshes in Mesh::MakeSimplicial and
  ParMesh::MakeSimplicial.

+- Added a new miniapp for interpolating a surface grid of points in 3D using a
+  smooth NURBS surface, that can then be sampled at arbitrary resolution while
+  staying close to the original geometry. See miniapps/nurbs/nurbs_surface.
+
 GPU computing
 -------------
 - The function Vector::SetSubVector(const Array<int> &, const real_t) now
@@ -41,6 +46,11 @@ GPU computing
  where host execution is always needed (e.g. when the DOFs array is small).
 - Introduced MFEM_FOREACH_THREAD_DIRECT, which directly maps loop tasks to GPU
  threads, assigning one task per thread.
+- Implemented a GPU-accelerated matrix-free AMR derefinement `GridFunction`
+  update operator. This supports mixed geometry meshes and variable order
+  spaces, and is the default derefinement operator constructed by
+  `FiniteElementSpace::Update` and `ParFiniteElementSpace::Update`.
+  The operator requires `FiniteElementSpace::Nonconforming() == true`.

 New and updated examples and miniapps
 -------------------------------------
@@ -598,14 +598,20 @@ set(MFEM_TPLS OPENMP HYPRE LAPACK BLAS SuperLUDist STRUMPACK METIS SuiteSparse
    BENCHMARK PARELAG TRIBOL MPI_CXX HIP HIPBLAS HIPSPARSE MOONOLITH BLITZ
    ALGOIM ENZYME)

-# Add all *_FOUND libraries in the variable TPL_LIBRARIES.
+# Add all created targets and *_FOUND libraries in the variables TPL_TARGETS and
+# TPL_LIBRARIES, respectively.
+set(TPL_TARGETS)
 set(TPL_LIBRARIES "")
 set(TPL_INCLUDE_DIRS "")
 foreach(TPL IN LISTS MFEM_TPLS)
-  if (${TPL}_FOUND)
+  if (${TPL}_FOUND OR TARGET ${TPL})
    message(STATUS "MFEM: using package ${TPL}")
-    list(APPEND TPL_LIBRARIES ${${TPL}_LIBRARIES})
-    list(APPEND TPL_INCLUDE_DIRS ${${TPL}_INCLUDE_DIRS})
+    if (TARGET ${TPL})
+      list(APPEND TPL_TARGETS ${TPL})
+    else()
+      list(APPEND TPL_LIBRARIES ${${TPL}_LIBRARIES})
+      list(APPEND TPL_INCLUDE_DIRS ${${TPL}_INCLUDE_DIRS})
+    endif()
  endif()
 endforeach(TPL)
 list(REVERSE TPL_LIBRARIES)
@@ -680,7 +686,10 @@ set(MFEM_INSTALL_DIR ${CMAKE_INSTALL_PREFIX})
 # Declaring the library
 mfem_add_library(mfem ${SOURCES} ${HEADERS} ${MASTER_HEADERS})
 # message(STATUS "TPL_LIBRARIES = ${TPL_LIBRARIES}")
-target_link_libraries(mfem PUBLIC ${TPL_LIBRARIES})
+target_link_libraries(mfem PUBLIC ${TPL_LIBRARIES} ${TPL_TARGETS})
+if (TPL_TARGETS)
+  add_dependencies(mfem ${TPL_TARGETS})
+endif()
 if (MINGW)
  target_link_libraries(mfem PRIVATE ws2_32)
 endif()
@@ -121,6 +121,11 @@ Parallel build:
   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-buil-dir> ; cd <mfem-build-dir>
+   cmake <mfem-source-dir> -DMFEM_USE_MPI=YES -DFETCH_TPLS=YES
+   make -j 4
+
 CUDA build:
   (this build requires CMake 3.17 or newer)
   mkdir <mfem-build-dir> ; cd <mfem-build-dir>
@@ -663,6 +668,7 @@ The specific libraries and their options are:
 - OpenMP (optional), usually part of compiler, used when either MFEM_USE_OPENMP
  or MFEM_USE_LEGACY_OPENMP is set to YES.
  Options: OPENMP_OPT, OPENMP_LIB.
+  Versions: OpenMP >= 3.1 when MFEM_USE_OPENMP=YES.

 - High-resolution POSIX clocks: when using MFEM_TIMER_TYPE = 2, it may be
  necessary to link with a system library (e.g. librt.so).
@@ -842,6 +848,7 @@ The specific libraries and their options are:
 - HIP (optional), used when MFEM_USE_HIP = YES.
  URL: https://rocmdocs.amd.com
  Options: HIP_CXX, HIP_ARCH, HIP_OPT, HIP_LIB.
+  Versions: ROCm >= 5.6.1.

 - OCCA (optional), used when MFEM_USE_OCCA = YES.
  URL: https://libocca.org
@@ -1074,6 +1081,9 @@ The following options are CMake specific:
 MFEM_ENABLE_TESTING  - Enable the ctest framework for testing.
 MFEM_ENABLE_EXAMPLES - Build all of the examples by default.
 MFEM_ENABLE_MINIAPPS - Build all of the miniapps by default.
+FETCH_TPLS           - Enable fetching of all supported third-party libraries.
+HYPRE_FETCH          - Enable fetching of hypre.
+METIS_FETCH          - Enable fetching of metis.

 External libraries (CMake):
 ---------------------------
@@ -1135,6 +1145,12 @@ The following built-in CMake packages are also used:
      set the <LIBNAME>_LIBRARIES option directly; the configuration option
      <LIBNAME>_DIR is not supported.

+The MFEM CMake build system also provides fetching (automated building) for the
+packages/libraries listed below. Note that when fetching is enabled, any related
+auto-detection functionality is disabled.
+
+ - HYPRE
+ - METIS

 Building without GNU make or CMake
 ==================================
@@ -9,15 +9,18 @@
 # terms of the BSD-3 license. We welcome feedback and contributions, see file
 # CONTRIBUTING.md for details.

-# Defines the following variables:
+# Defines the following variables if fetching of TPLs is disabled (default):
 #   - HYPRE_FOUND
 #   - HYPRE_LIBRARIES
 #   - HYPRE_INCLUDE_DIRS
 #   - HYPRE_VERSION
 #   - HYPRE_USING_CUDA (internal)
 #   - HYPRE_USING_HIP (internal)
+# otherwise, the following are defined:
+#   - HYPRE (imported library target)
+#   - HYPRE_VERSION (cache variable)

-if (HYPRE_FOUND)
+if (HYPRE_FOUND OR TARGET HYPRE)
  if (HYPRE_USING_CUDA)
    find_package(CUDAToolkit REQUIRED)
  endif()
@@ -33,6 +36,55 @@ if (HYPRE_FOUND)
  endif()
 endif()

+if (HYPRE_FETCH OR FETCH_TPLS)
+  set(HYPRE_FETCH_VERSION 2.33.0)
+  add_library(HYPRE STATIC IMPORTED)
+  # set options and associated dependencies
+  set(CMAKE_OPTIONS)
+  list(APPEND CMAKE_OPTIONS -DCMAKE_BUILD_TYPE:STRING=${CMAKE_BUILD_TYPE})
+  if (MFEM_USE_CUDA)
+    list(APPEND CMAKE_OPTIONS -DHYPRE_WITH_CUDA:BOOL=ON)
+    find_package(CUDAToolkit REQUIRED)
+    target_link_libraries(HYPRE INTERFACE CUDA::cusparse CUDA::curand CUDA::cublas)
+  elseif (MFEM_USE_HIP)
+    list(APPEND CMAKE_OPTIONS -DHYPRE_WITH_HIP:BOOL=ON)
+    find_package(rocsparse REQUIRED)
+    find_package(rocrand REQUIRED)
+    target_link_libraries(HYPRE INTERFACE rocsparse rocrand)
+  endif()
+  if (MFEM_USE_SINGLE)
+    list(APPEND CMAKE_OPTIONS -DHYPRE_ENABLE_SINGLE:BOOL=ON)
+  endif()
+  # 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 HYPRE ${HYPRE_FETCH_VERSION} to be built with ${CMAKE_OPTIONS}")
+  set(PREFIX ${CMAKE_BINARY_DIR}/fetch/hypre)
+  include(ExternalProject)
+  ExternalProject_Add(hypre
+    GIT_REPOSITORY https://github.com/hypre-space/hypre.git
+    GIT_TAG v${HYPRE_FETCH_VERSION}
+    GIT_SHALLOW TRUE
+    UPDATE_DISCONNECTED TRUE
+    SOURCE_SUBDIR src
+    PREFIX ${PREFIX}
+    CMAKE_CACHE_ARGS -DCMAKE_INSTALL_PREFIX:PATH=${PREFIX} -DCMAKE_INSTALL_LIBDIR:PATH=lib ${CMAKE_OPTIONS})
+  file(MAKE_DIRECTORY ${PREFIX}/include)
+  # set imported library target properties
+  add_dependencies(HYPRE hypre)
+  set_target_properties(HYPRE PROPERTIES
+    IMPORTED_LOCATION ${PREFIX}/lib/libHYPRE.a
+    INTERFACE_INCLUDE_DIRECTORIES ${PREFIX}/include)
+  # convert HYPRE version to integer
+  string(REGEX MATCHALL "[0-9]+" HYPRE_SPLIT_VERSION ${HYPRE_FETCH_VERSION})
+  list(GET HYPRE_SPLIT_VERSION 0 HYPRE_MAJOR_VERSION)
+  list(GET HYPRE_SPLIT_VERSION 1 HYPRE_MINOR_VERSION)
+  list(GET HYPRE_SPLIT_VERSION 2 HYPRE_PATCH_VERSION)
+  math(EXPR HYPRE_VERSION "10000*${HYPRE_MAJOR_VERSION} + 100*${HYPRE_MINOR_VERSION} + ${HYPRE_PATCH_VERSION}")
+  # set cache variables that would otherwise be set after mfem_find_package call
+  set(HYPRE_VERSION ${HYPRE_VERSION} CACHE STRING "HYPRE version." FORCE)
+  return()
+endif()
+
 include(MfemCmakeUtilities)
 mfem_find_package(HYPRE HYPRE HYPRE_DIR "include" "HYPRE.h" "lib" "HYPRE"
  "Paths to headers required by HYPRE." "Libraries required by HYPRE."
@@ -9,10 +9,38 @@
 # terms of the BSD-3 license. We welcome feedback and contributions, see file
 # CONTRIBUTING.md for details.

-# Defines the following variables:
+# Defines the following variables if fetching of TPLs is disabled (default):
 #   - METIS_FOUND
 #   - METIS_LIBRARIES
 #   - METIS_INCLUDE_DIRS
+#   - METIS_VERSION_5
+# otherwise, the following are defined:
+#   - METIS (imported library target)
+#   - METIS_VERSION_5 (cache variable)
+
+if (METIS_FETCH OR FETCH_TPLS)
+  set(METIS_FETCH_VERSION 4.0.3)
+  add_library(METIS STATIC IMPORTED)
+  # define external project
+  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
+    GIT_REPOSITORY https://github.com/mfem/tpls
+    GIT_TAG b60352fbe9675d374b00828055e55be4584c7995 # tag from 1/16/25
+    GIT_SHALLOW TRUE
+    UPDATE_DISCONNECTED TRUE
+    PREFIX ${PREFIX}
+    CONFIGURE_COMMAND tar -xzf ../metis/metis-${METIS_FETCH_VERSION}-mac.tgz --strip=1
+    INSTALL_COMMAND mkdir -p ${PREFIX}/lib && cp libmetis.a ${PREFIX}/lib/)
+  # set imported library target properties
+  add_dependencies(METIS metis)
+  set_target_properties(METIS PROPERTIES
+    IMPORTED_LOCATION ${PREFIX}/lib/libmetis.a)
+  # set cache variables that would otherwise be set after mfem_find_package call
+  set(METIS_VERSION_5 FALSE CACHE BOOL "Is METIS version 5?")
+  return()
+endif()

 include(MfemCmakeUtilities)
 mfem_find_package(METIS METIS METIS_DIR "include;Lib" "metis.h"
@@ -23,11 +23,14 @@
 #include "_config.hpp"
 #endif

+#include <cstdint>
+#include <climits>
+
 namespace mfem
 {

 #if (defined(MFEM_USE_CUDA) && defined(__CUDACC__)) || \
-    (defined(MFEM_USE_HIP) && defined(__HIPCC__))
+    (defined(MFEM_USE_HIP) && defined(__HIP__))
 #define MFEM_HOST_DEVICE __host__ __device__
 #else
 #define MFEM_HOST_DEVICE
@@ -89,6 +89,12 @@ option(MFEM_ENABLE_EXAMPLES "Build all of the examples" OFF)
 option(MFEM_ENABLE_MINIAPPS "Build all of the miniapps" OFF)
 option(MFEM_ENABLE_BENCHMARKS "Build all of the benchmarks" OFF)

+# Allow a user to specify fetching of certain third-party libraries instead of
+# searching for existing installations.
+option(FETCH_TPLS "Enable fetching of all supported third-party libraries" OFF)
+option(HYPRE_FETCH "Enable fetching of hypre" OFF)
+option(METIS_FETCH "Enable fetching of METIS" OFF)
+
 # Setting CXX/MPICXX on the command line or in user.cmake will overwrite the
 # autodetected C++ compiler.
 # set(CXX g++)
@@ -82,6 +82,8 @@ set(SRCS
  fe/fe_ser.cpp
  fe_coll.cpp
  fespace.cpp
+  derefmat_op.cpp
+  pderefmat_op.cpp
  geom.cpp
  gridfunc.cpp
  hybridization.cpp
@@ -255,6 +255,8 @@ PABilinearFormExtension::PABilinearFormExtension(BilinearForm *form)

 void PABilinearFormExtension::SetupRestrictionOperators(const L2FaceValues m)
 {
+   MFEM_PERF_FUNCTION;
+
   if ( Device::Allows(Backend::CEED_MASK) ) { return; }
   ElementDofOrdering ordering = GetEVectorOrdering(*a->FESpace());
   elem_restrict = trial_fes->GetElementRestriction(ordering);
@@ -266,11 +268,7 @@ void PABilinearFormExtension::SetupRestrictionOperators(const L2FaceValues m)

      // Gather the attributes on the host from all the elements
      const Mesh &mesh = *trial_fes->GetMesh();
-      elem_attributes.SetSize(mesh.GetNE());
-      for (int i = 0; i < mesh.GetNE(); ++i)
-      {
-         elem_attributes[i] = mesh.GetAttribute(i);
-      }
+      elem_attributes = &mesh.GetElementAttributes();
   }

   // Construct face restriction operators only if the bilinear form has
@@ -329,50 +327,14 @@ void PABilinearFormExtension::SetupRestrictionOperators(const L2FaceValues m)
         bdr_face_dYdn.SetSize(bdr_face_restrict_lex->Height());
      }

-      const Mesh &mesh = *trial_fes->GetMesh();
-      // See LinearFormExtension::Update for explanation of f_to_be logic.
-      std::unordered_map<int,int> f_to_be;
-      for (int i = 0; i < mesh.GetNBE(); ++i)
-      {
-         const int f = mesh.GetBdrElementFaceIndex(i);
-         f_to_be[f] = i;
-      }
-      const int nf_bdr = trial_fes->GetNFbyType(FaceType::Boundary);
-      bdr_attributes.SetSize(nf_bdr);
-      int f_ind = 0;
-      int missing_bdr_elems = 0;
-      for (int f = 0; f < mesh.GetNumFaces(); ++f)
-      {
-         if (!mesh.GetFaceInformation(f).IsOfFaceType(FaceType::Boundary))
-         {
-            continue;
-         }
-         int attribute = 1; // default value
-         if (f_to_be.find(f) != f_to_be.end())
-         {
-            const int be = f_to_be[f];
-            attribute = mesh.GetBdrAttribute(be);
-         }
-         else
-         {
-            // If a boundary face does not correspond to the a boundary element,
-            // we assign it the default attribute of 1. We also generate a
-            // warning at runtime with the number of such missing elements.
-            ++missing_bdr_elems;
-         }
-         bdr_attributes[f_ind] = attribute;
-         ++f_ind;
-      }
-      if (missing_bdr_elems)
-      {
-         MFEM_WARNING("Missing " << missing_bdr_elems << " boundary elements "
-                      "for boundary faces.");
-      }
+      bdr_face_attributes = &trial_fes->GetMesh()->GetBdrFaceAttributes();
   }
 }

 void PABilinearFormExtension::Assemble()
 {
+   MFEM_PERF_FUNCTION;
+
   SetupRestrictionOperators(L2FaceValues::DoubleValued);

   Array<BilinearFormIntegrator*> &integrators = *a->GetDBFI();
@@ -429,7 +391,7 @@ void PABilinearFormExtension::AssembleDiagonal(Vector &y) const
         mfem::forall(ne, [=] MFEM_HOST_DEVICE (int e)
         {
            const int attr = d_attr[e];
-            if (d_m[attr - 1] == 0)
+            if (attr <= 0 || d_m[attr - 1] == 0)
            {
               for (int i = 0; i < nd; ++i)
               {
@@ -450,7 +412,7 @@ void PABilinearFormExtension::AssembleDiagonal(Vector &y) const
         for (int i = 0; i < iSz; ++i)
         {
            assemble_diagonal_with_markers(*integrators[i], elem_markers[i],
-                                           elem_attributes, localY);
+                                           *elem_attributes, localY);
         }
         const ElementRestriction* H1elem_restrict =
            dynamic_cast<const ElementRestriction*>(elem_restrict);
@@ -476,7 +438,7 @@ void PABilinearFormExtension::AssembleDiagonal(Vector &y) const
      for (int i = 0; i < iSz; ++i)
      {
         assemble_diagonal_with_markers(*integrators[i], elem_markers[i],
-                                        elem_attributes, y);
+                                        *elem_attributes, y);
      }
   }

@@ -489,7 +451,7 @@ void PABilinearFormExtension::AssembleDiagonal(Vector &y) const
      for (int i = 0; i < n_bdr_integs; ++i)
      {
         assemble_diagonal_with_markers(*bdr_integs[i], bdr_markers[i],
-                                        bdr_attributes, bdr_face_Y);
+                                        *bdr_face_attributes, bdr_face_Y);
      }
      bdr_face_restrict_lex->AddAbsMultTranspose(bdr_face_Y, y);
   }
@@ -529,6 +491,8 @@ void PABilinearFormExtension::FormLinearSystem(const Array<int> &ess_tdof_list,
 void PABilinearFormExtension::MultInternal(const Vector &x, Vector &y,
                                           const bool useAbs) const
 {
+   MFEM_PERF_FUNCTION;
+
   Array<BilinearFormIntegrator*> &integrators = *a->GetDBFI();

   const int iSz = integrators.Size();
@@ -588,7 +552,7 @@ void PABilinearFormExtension::MultInternal(const Vector &x, Vector &y,
         for (int i = 0; i < iSz; ++i)
         {
            AddMultWithMarkers(*integrators[i], localX, elem_markers[i],
-                               elem_attributes, false, localY, useAbs);
+                               *elem_attributes, false, localY, useAbs);
         }
         if (H1elem_restrict && useAbs)
         {
@@ -690,8 +654,8 @@ void PABilinearFormExtension::MultInternal(const Vector &x, Vector &y,
         }
         for (int i = 0; i < n_bdr_integs; ++i)
         {
-            AddMultWithMarkers(*bdr_integs[i], bdr_face_X, bdr_markers[i], bdr_attributes,
-                               false, bdr_face_Y);
+            AddMultWithMarkers(*bdr_integs[i], bdr_face_X, bdr_markers[i],
+                               *bdr_face_attributes, false, bdr_face_Y);
         }
         for (int i = 0; i < n_bdr_face_integs; ++i)
         {
@@ -699,12 +663,14 @@ void PABilinearFormExtension::MultInternal(const Vector &x, Vector &y,
            {
               AddMultNormalDerivativesWithMarkers(
                  *bdr_face_integs[i], bdr_face_X, bdr_face_dXdn,
-                  bdr_face_markers[i], bdr_attributes, bdr_face_Y, bdr_face_dYdn);
+                  bdr_face_markers[i], *bdr_face_attributes, bdr_face_Y,
+                  bdr_face_dYdn);
            }
            else
            {
-               AddMultWithMarkers(*bdr_face_integs[i], bdr_face_X, bdr_face_markers[i],
-                                  bdr_attributes, false, bdr_face_Y);
+               AddMultWithMarkers(*bdr_face_integs[i], bdr_face_X,
+                                  bdr_face_markers[i], *bdr_face_attributes, false,
+                                  bdr_face_Y);
            }
         }
         bdr_face_restrict_lex->AddMultTransposeInPlace(bdr_face_Y, y);
@@ -727,7 +693,7 @@ void PABilinearFormExtension::MultTranspose(const Vector &x, Vector &y) const
      localY = 0.0;
      for (int i = 0; i < iSz; ++i)
      {
-         AddMultWithMarkers(*integrators[i], localX, elem_markers[i], elem_attributes,
+         AddMultWithMarkers(*integrators[i], localX, elem_markers[i], *elem_attributes,
                            true, localY);
      }
      elem_restrict->MultTranspose(localY, y);
@@ -774,13 +740,14 @@ void PABilinearFormExtension::MultTranspose(const Vector &x, Vector &y) const
         bdr_face_Y = 0.0;
         for (int i = 0; i < n_bdr_integs; ++i)
         {
-            AddMultWithMarkers(*bdr_integs[i], bdr_face_X, bdr_markers[i], bdr_attributes,
-                               true, bdr_face_Y);
+            AddMultWithMarkers(*bdr_integs[i], bdr_face_X, bdr_markers[i],
+                               *bdr_face_attributes, true, bdr_face_Y);
         }
         for (int i = 0; i < n_bdr_face_integs; ++i)
         {
-            AddMultWithMarkers(*bdr_face_integs[i], bdr_face_X, bdr_face_markers[i],
-                               bdr_attributes, true, bdr_face_Y);
+            AddMultWithMarkers(*bdr_face_integs[i], bdr_face_X,
+                               bdr_face_markers[i], *bdr_face_attributes, true,
+                               bdr_face_Y);
         }
         bdr_face_restrict_lex->AddMultTransposeInPlace(bdr_face_Y, y);
      }
@@ -804,7 +771,7 @@ static void AddWithMarkers_(
   mfem::forall(ne, [=] MFEM_HOST_DEVICE (int e)
   {
      const int attr = d_attr[e];
-      if (d_m[attr - 1] == 0) { return; }
+      if (attr <= 0 || d_m[attr - 1] == 0) { return; }
      for (int i = 0; i < nd; ++i)
      {
         d_y(i, e) += d_x(i, e);
@@ -898,6 +865,7 @@ EABilinearFormExtension::EABilinearFormExtension(BilinearForm *form)

 void EABilinearFormExtension::Assemble()
 {
+   MFEM_PERF_FUNCTION;
   SetupRestrictionOperators(L2FaceValues::SingleValued);

   ne = trial_fes->GetMesh()->GetNE();
@@ -920,7 +888,8 @@ void EABilinearFormExtension::Assemble()
      {
         const int i = idx % sz;
         const int e = idx / sz;
-         const real_t val = d_m[d_a[e] - 1] ? d_ea_1(i, e) : 0.0;
+         const real_t val =
+            d_a[e] > 0 ? (d_m[d_a[e] - 1] ? d_ea_1(i, e) : 0) : 0;
         if (add)
         {
            d_ea_2(i, e) += val;
@@ -953,7 +922,7 @@ void EABilinearFormExtension::Assemble()
            ea_data_tmp.SetSize(ea_data.Size());
            integrators[i]->AssembleEA(*a->FESpace(), ea_data_tmp, false);
            add_with_markers(ea_data_tmp, ea_data, ne, *markers,
-                             elem_attributes, add);
+                             *elem_attributes, add);
         }
      }
   }
@@ -982,7 +951,7 @@ void EABilinearFormExtension::Assemble()
            ea_data_tmp.SetSize(ea_data_bdr.Size());
            bdr_integs[i]->AssembleEABoundary(*a->FESpace(), ea_data_tmp, add);
            add_with_markers(ea_data_tmp, ea_data_bdr, nf_bdr, *markers,
-                             bdr_attributes, add);
+                             *bdr_face_attributes, add);
         }
      }
   }
@@ -1031,7 +1000,7 @@ void EABilinearFormExtension::Assemble()
                                                        ea_data_tmp,
                                                        add);
            add_with_markers(ea_data_tmp, ea_data_bdr, nf_bdr, *markers,
-                             bdr_attributes, add);
+                             *bdr_face_attributes, add);
         }
      }
   }
@@ -1445,6 +1414,7 @@ FABilinearFormExtension::FABilinearFormExtension(BilinearForm *form)

 void FABilinearFormExtension::Assemble()
 {
+   MFEM_PERF_FUNCTION;
   EABilinearFormExtension::Assemble();
   FiniteElementSpace &fes = *a->FESpace();
   int width = fes.GetVSize();
@@ -69,7 +69,8 @@ class PABilinearFormExtension : public BilinearFormExtension
 protected:
   const FiniteElementSpace *trial_fes, *test_fes; // Not owned
   /// Attributes of all mesh elements.
-   Array<int> elem_attributes, bdr_attributes;
+   const Array<int> *elem_attributes; // Not owned
+   const Array<int> *bdr_face_attributes; // Not owned
   mutable Vector tmp_evec; // Work array
   mutable Vector localX, localY;
   mutable Vector int_face_X, int_face_Y;
@@ -812,7 +812,7 @@ protected:
      const FiniteElement & test_fe) const
   {
      return (trial_fe.GetDim() == 1 && test_fe.GetDim() == 1 &&
-              trial_fe.GetDerivType() == mfem::FiniteElement::GRAD  &&
+              trial_fe.GetDerivType() == mfem::FiniteElement::GRAD &&
              test_fe.GetRangeType()  == mfem::FiniteElement::SCALAR );
   }

@@ -884,7 +884,7 @@ protected:
      const FiniteElement & trial_fe,
      const FiniteElement & test_fe) const
   {
-      return (trial_fe.GetDerivType() == mfem::FiniteElement::DIV  &&
+      return (trial_fe.GetDerivType() == mfem::FiniteElement::DIV &&
              test_fe.GetRangeType()  == mfem::FiniteElement::SCALAR );
   }

@@ -919,7 +919,7 @@ protected:
      const FiniteElement & trial_fe,
      const FiniteElement & test_fe) const
   {
-      return (trial_fe.GetDerivType() == mfem::FiniteElement::DIV  &&
+      return (trial_fe.GetDerivType() == mfem::FiniteElement::DIV &&
              test_fe.GetRangeType()  == mfem::FiniteElement::VECTOR );
   }

@@ -1600,7 +1600,7 @@ public:
   {
      return (trial_fe.GetCurlDim() == 3 && test_fe.GetRangeDim() == 3 &&
              trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR &&
-              trial_fe.GetDerivType() == mfem::FiniteElement::CURL   &&
+              trial_fe.GetDerivType() == mfem::FiniteElement::CURL  &&
              test_fe.GetRangeType()  == mfem::FiniteElement::VECTOR );
   }

@@ -1635,7 +1635,7 @@ public:
   {
      return (trial_fe.GetDim() == 2 && test_fe.GetDim() == 2 &&
              trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR &&
-              trial_fe.GetDerivType() == mfem::FiniteElement::CURL   &&
+              trial_fe.GetDerivType() == mfem::FiniteElement::CURL  &&
              test_fe.GetRangeType()  == mfem::FiniteElement::VECTOR );
   }

@@ -1669,7 +1669,7 @@ public:
   {
      return (trial_fe.GetDim() == 2 && test_fe.GetDim() == 2 &&
              trial_fe.GetRangeType() == mfem::FiniteElement::SCALAR &&
-              trial_fe.GetDerivType() == mfem::FiniteElement::GRAD   &&
+              trial_fe.GetDerivType() == mfem::FiniteElement::GRAD  &&
              test_fe.GetRangeType()  == mfem::FiniteElement::SCALAR );
   }

@@ -1760,7 +1760,7 @@ public:
      const FiniteElement & test_fe) const
   {
      return (trial_fe.GetRangeType() == mfem::FiniteElement::SCALAR &&
-              trial_fe.GetDerivType() == mfem::FiniteElement::GRAD   &&
+              trial_fe.GetDerivType() == mfem::FiniteElement::GRAD  &&
              test_fe.GetRangeType()  == mfem::FiniteElement::SCALAR );
   }

@@ -1793,7 +1793,7 @@ public:
      const FiniteElement & test_fe) const
   {
      return (trial_fe.GetRangeType() == mfem::FiniteElement::SCALAR &&
-              trial_fe.GetDerivType() == mfem::FiniteElement::GRAD   &&
+              trial_fe.GetDerivType() == mfem::FiniteElement::GRAD  &&
              test_fe.GetRangeType()  == mfem::FiniteElement::VECTOR &&
              test_fe.GetDerivType()  == mfem::FiniteElement::DIV   );
   }
@@ -1832,7 +1832,7 @@ public:
      const FiniteElement & test_fe) const
   {
      return (trial_fe.GetRangeType() == mfem::FiniteElement::VECTOR &&
-              trial_fe.GetDerivType() == mfem::FiniteElement::DIV    &&
+              trial_fe.GetDerivType() == mfem::FiniteElement::DIV   &&
              test_fe.GetRangeType()  == mfem::FiniteElement::SCALAR &&
              test_fe.GetDerivType()  == mfem::FiniteElement::GRAD
             );
@@ -1973,7 +1973,7 @@ protected:
      const FiniteElement & test_fe) const override
   {
      return (trial_fe.GetCurlDim() == 3 && test_fe.GetRangeDim() == 3 &&
-              trial_fe.GetDerivType() == mfem::FiniteElement::CURL  &&
+              trial_fe.GetDerivType() == mfem::FiniteElement::CURL &&
              test_fe.GetRangeType()  == mfem::FiniteElement::VECTOR );
   }

@@ -50,6 +50,7 @@ ElementTransformation *RefinedToCoarse(

 void Coefficient::Project(QuadratureFunction &qf)
 {
+   MFEM_PERF_FUNCTION;
   QuadratureSpaceBase &qspace = *qf.GetSpace();
   const int ne = qspace.GetNE();
   Vector values;
@@ -912,7 +912,7 @@ ConduitDataCollection::GridFunctionToBlueprintField(mfem::GridFunction *gf,

   if (vdim == 1) // scalar case
   {
-      n_field["values"].set_external(gf->GetData(),
+      n_field["values"].set_external(const_cast<real_t *>(gf->HostRead()),
                                     ndofs);
   }
   else // vector case
@@ -925,18 +925,18 @@ ConduitDataCollection::GridFunctionToBlueprintField(mfem::GridFunction *gf,
      int vdim_stride  = (ordering == Ordering::byNODES ? ndofs : 1);

      index_t offset = 0;
-      index_t stride = sizeof(double) * entry_stride;
+      index_t stride = sizeof(real_t) * entry_stride;

      for (int d = 0;  d < vdim; d++)
      {
         std::ostringstream oss;
         oss << "v" << d;
         std::string comp_name = oss.str();
-         n_field["values"][comp_name].set_external(gf->GetData(),
+         n_field["values"][comp_name].set_external(const_cast<real_t *>(gf->HostRead()),
                                                   ndofs,
                                                   offset,
                                                   stride);
-         offset +=  sizeof(double) * vdim_stride;
+         offset +=  sizeof(real_t) * vdim_stride;
      }
   }

@@ -0,0 +1,266 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "derefmat_op.hpp"
+#include "fes_kernels.hpp"
+
+/// \cond DO_NOT_DOCUMENT
+namespace mfem
+{
+
+namespace internal
+{
+template <Ordering::Type Order, bool Atomic>
+static void DerefMultKernelImpl(const DerefineMatrixOp &op, const Vector &x,
+                                Vector &y)
+{
+   DerefineMatrixOpMultFunctor<Order, Atomic> func;
+   func.xptr = x.Read();
+   y.UseDevice();
+   y = 0.;
+   func.yptr = y.ReadWrite();
+   func.bsptr = op.block_storage.Read();
+   func.boptr = op.block_offsets.Read();
+   func.brptr = op.block_row_idcs_offsets.Read();
+   func.bcptr = op.block_col_idcs_offsets.Read();
+   func.rptr = op.row_idcs.Read();
+   func.cptr = op.col_idcs.Read();
+   func.vdims = op.fespace->GetVDim();
+   func.nblocks = op.block_offsets.Size();
+   func.width = op.Width() / func.vdims;
+   func.height = op.Height() / func.vdims;
+   func.Run(op.max_rows);
+}
+
+} // namespace internal
+
+DerefineMatrixOp::DerefineMatrixOp(FiniteElementSpace &fespace_, int old_ndofs,
+                                   const Table *old_elem_dof,
+                                   const Table *old_elem_fos)
+   : Operator(fespace_.GetVSize(), old_ndofs * fespace_.GetVDim()),
+     fespace(&fespace_)
+{
+   static Kernels kernels;
+   constexpr int max_team_size = 256;
+   /// TODO: Implement DofTransformation support
+
+   MFEM_VERIFY(fespace->Nonconforming(),
+               "Not implemented for conforming meshes.");
+   MFEM_VERIFY(old_ndofs, "Missing previous (finer) space.");
+   MFEM_VERIFY(fespace->GetNDofs() <= old_ndofs,
+               "Previous space is not finer.");
+
+   const CoarseFineTransformations &dtrans =
+      fespace->GetMesh()->ncmesh->GetDerefinementTransforms();
+
+   MFEM_ASSERT(dtrans.embeddings.Size() == old_elem_dof->Size(), "");
+
+   const bool is_dg = fespace->FEColl()->GetContType()
+                      == FiniteElementCollection::DISCONTINUOUS;
+   DenseMatrix localRVO; // for variable-order only
+
+   DenseTensor localR[Geometry::NumGeom];
+   int total_rows = 0;
+   int total_cols = 0;
+   block_offsets.SetSize(dtrans.embeddings.Size());
+   block_offsets.HostWrite();
+   if (fespace->IsVariableOrder())
+   {
+      // TODO: any potential for some compression here?
+      // determine storage size and offsets
+      block_offsets[0] = 0;
+      int total_size = 0;
+      for (int k = 0; k < dtrans.embeddings.Size(); ++k)
+      {
+         const Embedding &emb = dtrans.embeddings[k];
+         const FiniteElement *fe = fespace->GetFE(emb.parent);
+         const int ldof = fe->GetDof();
+         if (k + 1 < dtrans.embeddings.Size())
+         {
+            block_offsets[k + 1] = block_offsets[k] + ldof * ldof;
+         }
+         total_rows += ldof;
+         total_cols += ldof;
+         total_size += ldof * ldof;
+      }
+      block_storage.SetSize(total_size);
+   }
+   else
+   {
+      // compression scheme:
+      // block_offsets is the start of each block, potentially repeated
+      // only need to store localR for used shapes
+      Mesh::GeometryList elem_geoms(*fespace->GetMesh());
+
+      int geom_offsets[Geometry::NumGeom];
+      {
+         int size = 0;
+         for (int i = 0; i < elem_geoms.Size(); ++i)
+         {
+            fespace->GetLocalDerefinementMatrices(elem_geoms[i],
+                                                  localR[elem_geoms[i]]);
+            geom_offsets[elem_geoms[i]] = size;
+            size += localR[elem_geoms[i]].TotalSize();
+         }
+         block_storage.SetSize(size);
+         // copy blocks into block_storage
+         auto bs_ptr = block_storage.HostWrite();
+         for (int i = 0; i < elem_geoms.Size(); ++i)
+         {
+            std::copy(localR[elem_geoms[i]].Data(),
+                      localR[elem_geoms[i]].Data()
+                      + localR[elem_geoms[i]].TotalSize(),
+                      bs_ptr);
+            bs_ptr += localR[elem_geoms[i]].TotalSize();
+         }
+      }
+      for (int k = 0; k < dtrans.embeddings.Size(); ++k)
+      {
+         const Embedding &emb = dtrans.embeddings[k];
+         Geometry::Type geom =
+            fespace->GetMesh()->GetElementBaseGeometry(emb.parent);
+
+         auto size = localR[geom].SizeI() * localR[geom].SizeJ();
+         total_rows += localR[geom].SizeI();
+         total_cols += localR[geom].SizeJ();
+         // set block offsets and sizes
+         block_offsets[k] = geom_offsets[geom] + size * emb.matrix;
+      }
+   }
+   row_idcs.SetSize(total_rows);
+   row_idcs.HostWrite();
+   col_idcs.SetSize(total_cols);
+   col_idcs.HostWrite();
+   block_row_idcs_offsets.SetSize(dtrans.embeddings.Size() + 1);
+   block_row_idcs_offsets.HostWrite();
+   block_col_idcs_offsets.SetSize(dtrans.embeddings.Size() + 1);
+   block_col_idcs_offsets.HostWrite();
+   block_row_idcs_offsets[0] = 0;
+   block_col_idcs_offsets[0] = 0;
+
+   // compute index information
+   Array<int> dofs, old_dofs;
+   max_rows = 1;
+
+   {
+      Array<int> mark(fespace->GetNDofs());
+      mark = 0;
+      auto bs_ptr = block_storage.HostWrite();
+      int ridx = 0;
+      int cidx = 0;
+      int num_marked = 0;
+      for (int k = 0; k < dtrans.embeddings.Size(); k++)
+      {
+         const Embedding &emb = dtrans.embeddings[k];
+         Geometry::Type geom =
+            fespace->GetMesh()->GetElementBaseGeometry(emb.parent);
+
+         if (fespace->IsVariableOrder())
+         {
+            const FiniteElement *fe = fespace->GetFE(emb.parent);
+            const DenseTensor &pmats = dtrans.point_matrices[geom];
+            const int ldof = fe->GetDof();
+
+            IsoparametricTransformation isotr;
+            isotr.SetIdentityTransformation(geom);
+
+            localRVO.SetSize(ldof, ldof);
+            isotr.SetPointMat(pmats(emb.matrix));
+            // Local restriction is size ldofxldof assuming that the parent
+            // and child are of same polynomial order.
+            fe->GetLocalRestriction(isotr, localRVO);
+            // copy block
+            auto size = localRVO.Height() * localRVO.Width();
+            std::copy(localRVO.Data(), localRVO.Data() + size, bs_ptr);
+            bs_ptr += size;
+         }
+         DenseMatrix &lR =
+            fespace->IsVariableOrder() ? localRVO : localR[geom](emb.matrix);
+         block_row_idcs_offsets[k + 1] =
+            block_row_idcs_offsets[k] + lR.Height();
+         block_col_idcs_offsets[k + 1] = block_col_idcs_offsets[k] + lR.Width();
+         max_rows = std::max(lR.Height(), max_rows);
+         // index information
+         fespace->elem_dof->GetRow(emb.parent, dofs);
+         old_elem_dof->GetRow(k, old_dofs);
+         MFEM_VERIFY(old_dofs.Size() == dofs.Size(),
+                     "Parent and child must have same #dofs.");
+         for (int i = 0; i < lR.Height(); ++i, ++ridx)
+         {
+            if (!std::isfinite(lR(i, 0)))
+            {
+               row_idcs[ridx] = INT_MAX;
+               continue;
+            }
+            int r = dofs[i];
+            int m = (r >= 0) ? r : (-1 - r);
+            if (is_dg || !mark[m])
+            {
+               row_idcs[ridx] = r;
+               mark[m] = 1;
+               ++num_marked;
+            }
+            else
+            {
+               row_idcs[ridx] = INT_MAX;
+            }
+         }
+         for (int i = 0; i < lR.Width(); ++i, ++cidx)
+         {
+            col_idcs[cidx] = old_dofs[i];
+         }
+      }
+      if (!is_dg && !fespace->IsVariableOrder())
+      {
+         MFEM_VERIFY(num_marked * fespace->GetVDim() == Height(),
+                     "internal error: not all rows were set.");
+      }
+   }
+   // if not using GPU, set max_rows/max_cols to zero
+   if (Device::Allows(Backend::DEVICE_MASK))
+   {
+      max_rows = std::min(max_rows, max_team_size);
+   }
+   else
+   {
+      max_rows = 1;
+   }
+}
+
+void DerefineMatrixOp::Mult(const Vector &x, Vector &y) const
+{
+   const bool is_dg = fespace->FEColl()->GetContType()
+                      == FiniteElementCollection::DISCONTINUOUS;
+   // DG needs atomic summation
+   MultKernel::Run(fespace->GetOrdering(), is_dg, *this, x, y);
+}
+
+DerefineMatrixOp::Kernels::Kernels()
+{
+   MultKernel::Specialization<Ordering::byNODES, false>::Add();
+   MultKernel::Specialization<Ordering::byVDIM, false>::Add();
+   MultKernel::Specialization<Ordering::byNODES, true>::Add();
+   MultKernel::Specialization<Ordering::byVDIM, true>::Add();
+}
+
+template <Ordering::Type Order, bool Atomic>
+DerefineMatrixOp::MultKernelType DerefineMatrixOp::MultKernel::Kernel()
+{
+   return internal::DerefMultKernelImpl<Order, Atomic>;
+}
+
+DerefineMatrixOp::MultKernelType
+DerefineMatrixOp::MultKernel::Fallback(Ordering::Type, bool)
+{
+   MFEM_ABORT("invalid MultKernel parameters");
+}
+} // namespace mfem
+/// \endcond DO_NOT_DOCUMENT
@@ -0,0 +1,65 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_DEREFMAT_OP
+#define MFEM_DEREFMAT_OP
+
+#include "fespace.hpp"
+
+#include "kernel_dispatch.hpp"
+
+namespace mfem
+{
+
+/// \cond DO_NOT_DOCUMENT
+
+struct DerefineMatrixOp : public Operator
+{
+   FiniteElementSpace *fespace;
+   /// offsets into block_storage
+   Array<int> block_offsets;
+   /// offsets into row_idcs
+   Array<int> block_row_idcs_offsets;
+   /// offsets into col_idcs
+   Array<int> block_col_idcs_offsets;
+   /// mapping for row dofs, INT_MAX indicates the block row should be ignored.
+   /// negative means the row data should be negated.
+   Array<int> row_idcs;
+   /// mapping for col dofs, negative means the col data should be negated.
+   Array<int> col_idcs;
+   /// dense block matrices which can be reused to construct the full matrix
+   /// operation. These are stored contiguously and blocks have no restrictions
+   /// on shape (can be rectangle and differ from block to block).
+   Vector block_storage;
+   /// maximum height of any block in block_storage for GPU
+   /// parallelization, or 1 for CPU runs.
+   int max_rows;
+
+   using MultKernelType = void (*)(const DerefineMatrixOp &, const Vector &,
+                                   Vector &);
+   /// template args: ordering, atomic
+   MFEM_REGISTER_KERNELS(MultKernel, MultKernelType, (Ordering::Type, bool));
+
+   struct Kernels
+   {
+      Kernels();
+   };
+
+   void Mult(const Vector &x, Vector &y) const;
+
+   DerefineMatrixOp(FiniteElementSpace &fespace_, int old_ndofs,
+                    const Table *old_elem_dof, const Table *old_elem_fos);
+};
+
+/// \endcond DO_NOT_DOCUMENT
+
+} // namespace mfem
+#endif
@@ -241,6 +241,7 @@ public:
   {
      MFEM_ASSERT(!action_callbacks.empty(), "no integrators have been set");
      prolongation(solutions, solutions_t, solutions_l);
+      residual_l = 0.0;
      for (auto &action : action_callbacks)
      {
         action(solutions_l, parameters_l, residual_l);
@@ -327,8 +327,8 @@ void print_mpi_sync(const std::string& msg)
   // First gather string lengths
   size_t msg_len = msg.length();
   std::vector<size_t> lengths(nranks);
-   MPI_Gather(&msg_len, 1, MPI_INT,
-              lengths.data(), 1, MPI_INT,
+   MPI_Gather(&msg_len, 1, MPITypeMap<size_t>::mpi_type,
+              lengths.data(), 1, MPITypeMap<size_t>::mpi_type,
              0, MPI_COMM_WORLD);

   if (myrank == 0)
@@ -568,7 +568,7 @@ struct ThreadBlocks
   int z = 1;
 };

-#if (defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP))
+#if defined(MFEM_USE_CUDA_OR_HIP)
 template <typename func_t>
 __global__ void forall_kernel_shmem(func_t f, int n)
 {
@@ -591,7 +591,7 @@ void forall(func_t f,
   if (Device::Allows(Backend::CUDA_MASK) ||
       Device::Allows(Backend::HIP_MASK))
   {
-#if (defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP))
+#if defined(MFEM_USE_CUDA_OR_HIP)
      // int gridsize = (N + Z - 1) / Z;
      int num_bytes = num_shmem * sizeof(decltype(shmem));
      dim3 block_size(blocks.x, blocks.y, blocks.z);
@@ -987,7 +987,7 @@ get_restriction_transpose(
   {
      auto RT = [=](const Vector &v_e, Vector &v_l)
      {
-         v_l = v_e;
+         v_l += v_e;
      };
      return std::make_tuple(RT, 1);
   }
@@ -996,7 +996,7 @@ get_restriction_transpose(
      const Operator *R = get_restriction<entity_t>(f, o);
      std::function<void(const Vector&, Vector&)> RT = [=](const Vector &x, Vector &y)
      {
-         R->MultTranspose(x, y);
+         R->AddMultTranspose(x, y);
      };
      return std::make_tuple(RT, R->Height());
   }
@@ -0,0 +1,249 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_FES_KERNELS_HPP
+#define MFEM_FES_KERNELS_HPP
+
+#include "../general/forall.hpp"
+
+#include <climits>
+
+namespace mfem
+{
+/// \cond DO_NOT_DOCUMENT
+namespace internal
+{
+
+///
+/// Implements matrix-vector multiply $y = A x$ for a sparse matrix composed of
+/// a sum of smaller dense blocks. There is additional permutation/sign
+/// information associated with each block. The base class only implements
+/// helper routines such as computing block widths, index into x, index into y,
+/// and column in A given sub-block information.
+/// @sa DerefineMatrixOpMultFunctor
+///
+/// @tparam Order vdim ordering for x and y. Note that for Diag = false this is
+/// ignored for x as x has a special interleaved order.
+/// @tparam Base used for the curious recurring template pattern (CRTP) so the
+/// base class can access child class fields without virtual functions
+/// @tparam Diag true if this corresponds to the diagonal block (coarse element
+/// and fine element are on our rank), false otherwise (coarse element is on our
+/// rank, fine element is on a different rank).
+///
+template <Ordering::Type Order, class Base, bool Diag = true>
+struct DerefineMatrixOpFunctorBase;
+
+template <class Base>
+struct DerefineMatrixOpFunctorBase<Ordering::byNODES, Base, true>
+{
+   /// block column indices offsets
+   const int *bcptr;
+   /// column indices
+   const int *cptr;
+
+   int MFEM_HOST_DEVICE BlockWidth(int k) const
+   {
+      return bcptr[k + 1] - bcptr[k];
+   }
+
+   void MFEM_HOST_DEVICE Col(int j, int k, int &col, int &sign) const
+   {
+      col = cptr[bcptr[k] + j];
+      if (col < 0)
+      {
+         col = -1 - col;
+         sign = -sign;
+      }
+   }
+
+   int MFEM_HOST_DEVICE IndexX(int col, int vdim, int) const
+   {
+      return col + vdim * static_cast<const Base *>(this)->width;
+   }
+   int MFEM_HOST_DEVICE IndexY(int row, int vdim) const
+   {
+      return row + vdim * static_cast<const Base *>(this)->height;
+   }
+};
+
+template <class Base>
+struct DerefineMatrixOpFunctorBase<Ordering::byVDIM, Base, true>
+{
+   /// block column indices offsets
+   const int *bcptr;
+   /// column indices
+   const int *cptr;
+
+   int MFEM_HOST_DEVICE BlockWidth(int k) const
+   {
+      return bcptr[k + 1] - bcptr[k];
+   }
+
+   void MFEM_HOST_DEVICE Col(int j, int k, int &col, int &sign) const
+   {
+      col = cptr[bcptr[k] + j];
+      if (col < 0)
+      {
+         col = -1 - col;
+         sign = -sign;
+      }
+   }
+
+   int MFEM_HOST_DEVICE IndexX(int col, int vdim, int) const
+   {
+      return vdim + col * static_cast<const Base *>(this)->vdims;
+   }
+   int MFEM_HOST_DEVICE IndexY(int row, int vdim) const
+   {
+      return vdim + row * static_cast<const Base *>(this)->vdims;
+   }
+};
+
+template <class Base>
+struct DerefineMatrixOpFunctorBase<Ordering::byNODES, Base, false>
+{
+   /// receive segment offsets
+   const int *segptr;
+   /// receive segment index
+   const int *rsptr;
+   /// off-diagonal block column offsets
+   const int *coptr;
+   /// off-diagonal block widths
+   const int *bwptr;
+
+   int MFEM_HOST_DEVICE BlockWidth(int k) const { return bwptr[k]; }
+
+   void MFEM_HOST_DEVICE Col(int j, int k, int &col, int &sign) const
+   {
+      col = coptr[k] + j;
+   }
+
+   int MFEM_HOST_DEVICE IndexX(int col, int vdim, int k) const
+   {
+      int tmp = rsptr[k];
+      int segwidth = segptr[tmp + 1] - segptr[tmp];
+      return segptr[tmp] * static_cast<const Base *>(this)->vdims + col +
+             vdim * segwidth;
+   }
+   int MFEM_HOST_DEVICE IndexY(int row, int vdim) const
+   {
+      return row + vdim * static_cast<const Base *>(this)->height;
+   }
+};
+
+template <class Base>
+struct DerefineMatrixOpFunctorBase<Ordering::byVDIM, Base, false>
+{
+   /// receive segment offsets
+   const int *segptr;
+   /// receive segment index
+   const int *rsptr;
+   /// off-diagonal block column offsets
+   const int *coptr;
+   /// off-diagonal block widths
+   const int *bwptr;
+
+   int MFEM_HOST_DEVICE BlockWidth(int k) const { return bwptr[k]; }
+
+   void MFEM_HOST_DEVICE Col(int j, int k, int &col, int &sign) const
+   {
+      col = coptr[k] + j;
+   }
+
+   int MFEM_HOST_DEVICE IndexX(int col, int vdim, int k) const
+   {
+      int tmp = rsptr[k];
+      int segwidth = segptr[tmp + 1] - segptr[tmp];
+      return segptr[tmp] * static_cast<const Base *>(this)->vdims + col +
+             vdim * segwidth;
+   }
+   int MFEM_HOST_DEVICE IndexY(int row, int vdim) const
+   {
+      return vdim + row * static_cast<const Base *>(this)->vdims;
+   }
+};
+
+/// internally used to implement the derefinement operator Mult diagonal
+/// block
+template <Ordering::Type Order, bool Atomic, bool Diag = true>
+struct DerefineMatrixOpMultFunctor
+   : public DerefineMatrixOpFunctorBase<
+     Order, DerefineMatrixOpMultFunctor<Order, Atomic, Diag>, Diag>
+{
+   const real_t *xptr;
+   real_t *yptr;
+   /// block storage
+   const real_t *bsptr;
+   /// block offsets
+   const int *boptr;
+   /// block row index offsets
+   const int *brptr;
+   /// row indices
+   const int *rptr;
+
+   // number of blocks
+   int nblocks;
+   // number of components
+   int vdims;
+   /// overall operator height (for vdim = 1)
+   int height;
+   /// overall operator width (for vdim = 1)
+   int width;
+   void MFEM_HOST_DEVICE operator()(int kidx) const
+   {
+      int k = kidx % nblocks;
+      int vdim = kidx / nblocks;
+
+      int block_height = brptr[k + 1] - brptr[k];
+      int block_width = this->BlockWidth(k);
+      MFEM_FOREACH_THREAD(i, x, block_height)
+      {
+         int row = rptr[brptr[k] + i];
+         int rsign = 1;
+         if (row < 0)
+         {
+            row = -1 - row;
+            rsign = -1;
+         }
+         if (row < INT_MAX)
+         {
+            // row not marked as unused
+            real_t sum = 0;
+            for (int j = 0; j < block_width; ++j)
+            {
+               int col, sign = rsign;
+               this->Col(j, k, col, sign);
+               sum += sign * bsptr[boptr[k] + i + j * block_height] *
+                      xptr[this->IndexX(col, vdim, k)];
+            }
+#if defined(__CUDA_ARCH__) or defined(__HIP_DEVICE_COMPILE__)
+            if (Atomic)
+            {
+               atomicAdd(yptr + this->IndexY(row, vdim), sum);
+            }
+            else
+#endif
+            {
+               yptr[this->IndexY(row, vdim)] += sum;
+            }
+         }
+      }
+   }
+
+   /// N is the max block row size (doesn't have to be a power of 2)
+   void Run(int N) const { forall_2D(nblocks * vdims, N, 1, *this); }
+};
+
+} // namespace internal
+/// \endcond DO_NOT_DOCUMENT
+} // namespace mfem
+
+#endif
@@ -17,6 +17,9 @@
 #include "fem.hpp"
 #include "ceed/interface/util.hpp"

+#include "derefmat_op.hpp"
+
+#include <algorithm>
 #include <cmath>
 #include <cstdarg>

@@ -24,9 +27,9 @@ using namespace std;

 namespace mfem
 {
-
-template <> void Ordering::
-DofsToVDofs<Ordering::byNODES>(int ndofs, int vdim, Array<int> &dofs)
+template <>
+void Ordering::DofsToVDofs<Ordering::byNODES>(int ndofs, int vdim,
+                                              Array<int> &dofs)
 {
   // static method
   int size = dofs.Size();
@@ -40,8 +43,9 @@ DofsToVDofs<Ordering::byNODES>(int ndofs, int vdim, Array<int> &dofs)
   }
 }

-template <> void Ordering::
-DofsToVDofs<Ordering::byVDIM>(int ndofs, int vdim, Array<int> &dofs)
+template <>
+void Ordering::DofsToVDofs<Ordering::byVDIM>(int ndofs, int vdim,
+                                             Array<int> &dofs)
 {
   // static method
   int size = dofs.Size();
@@ -55,7 +59,6 @@ DofsToVDofs<Ordering::byVDIM>(int ndofs, int vdim, Array<int> &dofs)
   }
 }

-
 FiniteElementSpace::FiniteElementSpace()
   : mesh(NULL), fec(NULL), vdim(0), ordering(Ordering::byNODES),
     ndofs(0), nvdofs(0), nedofs(0), nfdofs(0), nbdofs(0),
@@ -98,7 +101,10 @@ FiniteElementSpace::FiniteElementSpace(const FiniteElementSpace &orig,
 FiniteElementSpace::FiniteElementSpace(Mesh *mesh,
                                       const FiniteElementCollection *fec,
                                       int vdim, int ordering)
-{ Constructor(mesh, NULL, fec, vdim, ordering); }
+{
+   MFEM_PERF_FUNCTION;
+   Constructor(mesh, NULL, fec, vdim, ordering);
+}

 FiniteElementSpace::FiniteElementSpace(Mesh *mesh, NURBSExtension *ext,
                                       const FiniteElementCollection *fec,
@@ -390,6 +396,8 @@ void FiniteElementSpace::BuildElementToDofTable() const
 {
   if (elem_dof) { return; }

+   MFEM_PERF_FUNCTION;
+
   // TODO: can we call GetElementDofs only once per element?
   Table *el_dof = new Table;
   Table *el_fos = (mesh->Dimension() > 2) ? (new Table) : NULL;
@@ -2745,6 +2753,8 @@ void FiniteElementSpace::BuildNURBSFaceToDofTable() const

 void FiniteElementSpace::Construct()
 {
+   MFEM_PERF_FUNCTION;
+
   // This method should be used only for non-NURBS spaces.
   MFEM_VERIFY(!NURBSext, "internal error");

@@ -4244,7 +4254,11 @@ void FiniteElementSpace::Update(bool want_transform)
         case Mesh::DEREFINE:
         {
            BuildConformingInterpolation();
+#if 0
            Th.Reset(DerefinementMatrix(old_ndofs, old_elem_dof, old_elem_fos));
+#else
+            Th.Reset(new DerefineMatrixOp(*this, old_ndofs, old_elem_dof, old_elem_fos));
+#endif
            if (IsVariableOrder())
            {
               if (cP && cR_hp)
@@ -113,7 +113,7 @@ class QuadratureSpace;
 class QuadratureInterpolator;
 class FaceQuadratureInterpolator;
 class PRefinementTransferOperator;
-
+struct DerefineMatrixOp;

 /** @brief Class FiniteElementSpace - responsible for providing FEM view of the
    mesh, mainly managing the set of degrees of freedom.
@@ -246,6 +246,7 @@ class FiniteElementSpace
   friend class PRefinementTransferOperator;
   friend void Mesh::Swap(Mesh &, bool);
   friend class LORBase;
+   friend struct DerefineMatrixOp;

 protected:
   /// The mesh that FE space lives on (not owned).
@@ -19,6 +19,7 @@
 #include "../mesh/nurbs.hpp"
 #include "../mesh/vtkhdf.hpp"
 #include "../general/text.hpp"
+#include "../general/reducers.hpp"

 #ifdef MFEM_USE_MPI
 #include "pfespace.hpp"
@@ -3326,8 +3327,126 @@ real_t GridFunction::ComputeLpError(const real_t p, Coefficient &exsol,
                                    const IntegrationRule *irs[],
                                    const Array<int> *elems) const
 {
+   MFEM_PERF_FUNCTION;
+
+   MFEM_VERIFY(fes->GetVDim() == 1, "invalid vector dimension!");
+
   real_t error = 0.0;
-   const FiniteElement *fe;
+
+   bool device_eval = true;
+   // TODO: check for cases that are not supported on device:
+   //       * mixed meshes
+   //       * meshes with non-tensor-product elements can have negative weights
+   //       * variable orders
+   //       * weight is not NULL
+   //       * elems is not NULL
+   //       * map type is not VALUE
+   //       * ...
+   Mesh *mesh = fes->GetMesh();
+   const FiniteElement *fe = fes->GetTypicalFE();
+   if (mesh->GetNumGeometries(mesh->Dimension()) > 1 ||
+       (mesh->Dimension() > 1 && mesh->MeshGenerator() != 2) ||
+       fes->IsVariableOrder() ||
+       weight != nullptr ||
+       elems != nullptr ||
+       fe->GetMapType() != FiniteElement::MapType::VALUE)
+   {
+      device_eval = false;
+   }
+   if (device_eval)
+   {
+      Geometry::Type geom = mesh->GetTypicalElementGeometry();
+      const IntegrationRule *ir_p;
+      if (irs)
+      {
+         ir_p = irs[geom];
+      }
+      else
+      {
+         int intorder = 2*fe->GetOrder() + 3; // <----------
+         ir_p = &(IntRules.Get(geom, intorder));
+      }
+      const IntegrationRule &ir = *ir_p;
+      QuadratureSpace qs(*mesh, ir);
+      CoefficientVector coeff(exsol, qs, CoefficientStorage::FULL);
+
+      const QVectorLayout ql = QVectorLayout::byNODES;
+      const MemoryType d_mt = MemoryType::DEFAULT;
+      Vector q_vals;
+      // TODO: make this a method
+      {
+         // const FiniteElement *fe = fes->GetTypicalFE();
+         const int vdim = fes->GetVDim();
+         const int NE   = fes->GetNE();
+         const int ND   = fe->GetDof();
+         const int NQ   = ir.GetNPoints();
+         MemoryType my_d_mt = (d_mt != MemoryType::DEFAULT) ? d_mt :
+                              Device::GetDeviceMemoryType();
+         // byNODES : NQPT x VDIM x NE
+         // byVDIM  : VDIM x NQPT x NE
+         q_vals.SetSize(vdim*NQ*NE, my_d_mt);
+         const QuadratureInterpolator &qi = *fes->GetQuadratureInterpolator(ir);
+         qi.SetOutputLayout(ql);
+         const bool use_tensor_products = UsesTensorBasis(*fes);
+         qi.DisableTensorProducts(!use_tensor_products);
+         const ElementDofOrdering e_ordering =
+            use_tensor_products ?
+            ElementDofOrdering::LEXICOGRAPHIC :
+            ElementDofOrdering::NATIVE;
+         const Operator *elem_restr = fes->GetElementRestriction(e_ordering);
+         if (fe->GetMapType() == FiniteElement::MapType::INTEGRAL)
+         {
+            // Pre-compute the geometric factors in order to set the desired
+            // MemoryType they use:
+            fes->GetMesh()->GetGeometricFactors(
+               ir, GeometricFactors::DETERMINANTS, my_d_mt);
+         }
+         if (elem_restr)
+         {
+            Vector f_e(vdim*ND*NE, my_d_mt);
+            elem_restr->Mult(*this, f_e);
+            qi.PhysValues(f_e, q_vals);
+         }
+         else
+         {
+            qi.PhysValues(*this, q_vals);
+         }
+      }
+
+      const real_t *exact_d = coeff.Read();
+      const real_t *gridf_d = q_vals.Read();
+      // FIXME: reuse the workspace vector from vector.cpp?
+      static Array<real_t> workspace;
+      if (p < infinity())
+      {
+         MemoryType my_d_mt = (d_mt != MemoryType::DEFAULT) ? d_mt :
+                              Device::GetDeviceMemoryType();
+         const GeometricFactors *geom_factors =
+            fes->GetMesh()->GetGeometricFactors(
+               ir, GeometricFactors::DETERMINANTS, my_d_mt);
+         const real_t *detJ_d = geom_factors->detJ.Read();
+         const real_t *w_d = ir.GetWeights().Read();
+         const int NQ = ir.GetNPoints();
+         mfem::reduce(q_vals.Size(), error,
+                      [=] MFEM_HOST_DEVICE(int i, real_t &r)
+         {
+            const real_t diff = fabs(exact_d[i] - gridf_d[i]);
+            r += w_d[i%NQ] * detJ_d[i] * pow(diff, p);
+         }, SumReducer<real_t> {}, true, workspace);
+         error = pow(error, 1./p);
+      }
+      else
+      {
+         mfem::reduce(q_vals.Size(), error,
+                      [=] MFEM_HOST_DEVICE(int i, real_t &r)
+         {
+            const real_t diff = fabs(exact_d[i] - gridf_d[i]);
+            r = fmax(r, diff);
+         }, MaxReducer<real_t> {}, true, workspace);
+      }
+      return error;
+   }
+
   ElementTransformation *T;
   Vector vals;

@@ -161,7 +161,8 @@ static void EADiffusionAssemble3D(const int NE,
   auto B = Reshape(b.Read(), Q1D, D1D);
   auto G = Reshape(g.Read(), Q1D, D1D);
   auto D = Reshape(padata.Read(), Q1D, Q1D, Q1D, 6, NE);
-   auto A = Reshape(eadata.ReadWrite(), D1D, D1D, D1D, D1D, D1D, D1D, NE);
+   auto A = Reshape(add ? eadata.ReadWrite() : eadata.Write(),
+                    D1D, D1D, D1D, D1D, D1D, D1D, NE);
   mfem::forall_3D(NE, D1D, D1D, D1D, [=] MFEM_HOST_DEVICE (int e)
   {
      const int D1D = T_D1D ? T_D1D : d1d;
@@ -246,58 +247,60 @@ void DiffusionIntegrator::AssembleEA(const FiniteElementSpace &fes,
                                     Vector &ea_data,
                                     const bool add)
 {
+   MFEM_PERF_FUNCTION;
   AssemblePA(fes);
   ne = fes.GetMesh()->GetNE();
   const Array<real_t> &B = maps->B;
   const Array<real_t> &G = maps->G;
+   decltype(&EADiffusionAssemble1D<>) kernel = nullptr;
   if (dim == 1)
   {
      switch ((dofs1D << 4 ) | quad1D)
      {
-         case 0x22: return EADiffusionAssemble1D<2,2>(ne,B,G,pa_data,ea_data,add);
-         case 0x33: return EADiffusionAssemble1D<3,3>(ne,B,G,pa_data,ea_data,add);
-         case 0x44: return EADiffusionAssemble1D<4,4>(ne,B,G,pa_data,ea_data,add);
-         case 0x55: return EADiffusionAssemble1D<5,5>(ne,B,G,pa_data,ea_data,add);
-         case 0x66: return EADiffusionAssemble1D<6,6>(ne,B,G,pa_data,ea_data,add);
-         case 0x77: return EADiffusionAssemble1D<7,7>(ne,B,G,pa_data,ea_data,add);
-         case 0x88: return EADiffusionAssemble1D<8,8>(ne,B,G,pa_data,ea_data,add);
-         case 0x99: return EADiffusionAssemble1D<9,9>(ne,B,G,pa_data,ea_data,add);
-         default:   return EADiffusionAssemble1D(ne,B,G,pa_data,ea_data,add,
-                                                    dofs1D,quad1D);
+         case 0x22: kernel = EADiffusionAssemble1D<2,2>;
+         case 0x33: kernel = EADiffusionAssemble1D<3,3>;
+         case 0x44: kernel = EADiffusionAssemble1D<4,4>;
+         case 0x55: kernel = EADiffusionAssemble1D<5,5>;
+         case 0x66: kernel = EADiffusionAssemble1D<6,6>;
+         case 0x77: kernel = EADiffusionAssemble1D<7,7>;
+         case 0x88: kernel = EADiffusionAssemble1D<8,8>;
+         case 0x99: kernel = EADiffusionAssemble1D<9,9>;
+         default:   kernel = EADiffusionAssemble1D<>;
      }
   }
   else if (dim == 2)
   {
      switch ((dofs1D << 4 ) | quad1D)
      {
-         case 0x22: return EADiffusionAssemble2D<2,2>(ne,B,G,pa_data,ea_data,add);
-         case 0x33: return EADiffusionAssemble2D<3,3>(ne,B,G,pa_data,ea_data,add);
-         case 0x44: return EADiffusionAssemble2D<4,4>(ne,B,G,pa_data,ea_data,add);
-         case 0x55: return EADiffusionAssemble2D<5,5>(ne,B,G,pa_data,ea_data,add);
-         case 0x66: return EADiffusionAssemble2D<6,6>(ne,B,G,pa_data,ea_data,add);
-         case 0x77: return EADiffusionAssemble2D<7,7>(ne,B,G,pa_data,ea_data,add);
-         case 0x88: return EADiffusionAssemble2D<8,8>(ne,B,G,pa_data,ea_data,add);
-         case 0x99: return EADiffusionAssemble2D<9,9>(ne,B,G,pa_data,ea_data,add);
-         default:   return EADiffusionAssemble2D(ne,B,G,pa_data,ea_data,add,
-                                                    dofs1D,quad1D);
+         case 0x22: kernel = EADiffusionAssemble2D<2,2>;
+         case 0x33: kernel = EADiffusionAssemble2D<3,3>;
+         case 0x44: kernel = EADiffusionAssemble2D<4,4>;
+         case 0x55: kernel = EADiffusionAssemble2D<5,5>;
+         case 0x66: kernel = EADiffusionAssemble2D<6,6>;
+         case 0x77: kernel = EADiffusionAssemble2D<7,7>;
+         case 0x88: kernel = EADiffusionAssemble2D<8,8>;
+         case 0x99: kernel = EADiffusionAssemble2D<9,9>;
+         default:   kernel = EADiffusionAssemble2D<>;
      }
   }
   else if (dim == 3)
   {
      switch ((dofs1D << 4 ) | quad1D)
      {
-         case 0x23: return EADiffusionAssemble3D<2,3>(ne,B,G,pa_data,ea_data,add);
-         case 0x34: return EADiffusionAssemble3D<3,4>(ne,B,G,pa_data,ea_data,add);
-         case 0x45: return EADiffusionAssemble3D<4,5>(ne,B,G,pa_data,ea_data,add);
-         case 0x56: return EADiffusionAssemble3D<5,6>(ne,B,G,pa_data,ea_data,add);
-         case 0x67: return EADiffusionAssemble3D<6,7>(ne,B,G,pa_data,ea_data,add);
-         case 0x78: return EADiffusionAssemble3D<7,8>(ne,B,G,pa_data,ea_data,add);
-         case 0x89: return EADiffusionAssemble3D<8,9>(ne,B,G,pa_data,ea_data,add);
-         default:   return EADiffusionAssemble3D(ne,B,G,pa_data,ea_data,add,
-                                                    dofs1D,quad1D);
+         case 0x23: kernel = EADiffusionAssemble3D<2,3>;
+         case 0x34: kernel = EADiffusionAssemble3D<3,4>;
+         case 0x45: kernel = EADiffusionAssemble3D<4,5>;
+         case 0x56: kernel = EADiffusionAssemble3D<5,6>;
+         case 0x67: kernel = EADiffusionAssemble3D<6,7>;
+         case 0x78: kernel = EADiffusionAssemble3D<7,8>;
+         case 0x89: kernel = EADiffusionAssemble3D<8,9>;
+         default:   kernel = EADiffusionAssemble3D<>;
      }
   }
-   MFEM_ABORT("Unknown kernel.");
+   MFEM_VERIFY(kernel != nullptr, "Unknown kernel.");
+   kernel(ne,B,G,pa_data,ea_data,add,dofs1D,quad1D);
+   // Free the PA data:
+   pa_data.Destroy();
 }

 }
@@ -39,6 +39,8 @@ void DiffusionIntegrator::AssembleDiagonalPA(Vector &diag)
 // PA Diffusion Apply kernel
 void DiffusionIntegrator::AddMultPA(const Vector &x, Vector &y) const
 {
+   MFEM_PERF_FUNCTION;
+
   if (DeviceCanUseCeed())
   {
      ceedOp->AddMult(x, y);
@@ -88,6 +90,8 @@ void DiffusionIntegrator::AddMultTransposePA(const Vector &x, Vector &y) const

 void DiffusionIntegrator::AssemblePA(const FiniteElementSpace &fes)
 {
+   MFEM_PERF_FUNCTION;
+
   const MemoryType mt = (pa_mt == MemoryType::DEFAULT) ?
                         Device::GetDeviceMemoryType() : pa_mt;
   // Assuming the same element type
@@ -23,6 +23,8 @@ namespace mfem

 void MassIntegrator::AssemblePA(const FiniteElementSpace &fes)
 {
+   MFEM_PERF_FUNCTION;
+
   const MemoryType mt = (pa_mt == MemoryType::DEFAULT) ?
                         Device::GetDeviceMemoryType() : pa_mt;

@@ -59,26 +61,23 @@ void MassIntegrator::AssemblePA(const FiniteElementSpace &fes)

   QuadratureSpace qs(*mesh, *ir);
   CoefficientVector coeff(Q, qs, CoefficientStorage::COMPRESSED);
-
-   const int NE = ne;
-   const int Q1D = quad1D;
-   const int NQ = static_cast<int>(std::pow(Q1D, dim));
-   const bool const_c = coeff.Size() == 1;
-   const bool by_val = map_type == FiniteElement::VALUE;
-   const auto W = Reshape(ir->GetWeights().Read(), NQ);
-   const auto J = Reshape(geom->detJ.Read(), NQ, NE);
-   const auto C = const_c ? Reshape(coeff.Read(), 1, 1) :
-                  Reshape(coeff.Read(), NQ,NE);
-   auto v = Reshape(pa_data.Write(), NQ, NE);
-   mfem::forall_2D(NE, NQ, 1, [=] MFEM_HOST_DEVICE (int e)
   {
-      MFEM_FOREACH_THREAD(i, x, NQ)
+      const int NE = ne;
+      const int NQ = nq;
+      const bool const_c = coeff.Size() == 1;
+      const bool by_val = map_type == FiniteElement::VALUE;
+      const auto W = Reshape(ir->GetWeights().Read(), NQ);
+      const auto J = Reshape(geom->detJ.Read(), NQ, NE);
+      const auto C =
+         const_c ? Reshape(coeff.Read(), 1, 1) : Reshape(coeff.Read(), NQ, NE);
+      auto v = Reshape(pa_data.Write(), NQ, NE);
+      mfem::forall(NQ, NE, [=] MFEM_HOST_DEVICE(int q, int e)
      {
-         const real_t detJ = J(i,e);
-         const real_t coeff = const_c ? C(0,0) : C(i,e);
-         v(i,e) =  W(i) * coeff * (by_val ? detJ : 1.0/detJ);
-      }
-   });
+         const real_t detJ = J(q, e);
+         const real_t coeff = const_c ? C(0, 0) : C(q, e);
+         v(q, e) = W(q) * coeff * (by_val ? detJ : 1.0 / detJ);
+      });
+   }
 }

 void MassIntegrator::AssemblePABoundary(const FiniteElementSpace &fes)
@@ -109,50 +108,22 @@ void MassIntegrator::AssemblePABoundary(const FiniteElementSpace &fes)
   CoefficientVector coeff(Q, qs, CoefficientStorage::COMPRESSED);

   const int NE = ne;
-   const int Q1D = quad1D;
+   const int NQ = nq;
   const bool const_c = coeff.Size() == 1;
   const bool by_val = map_type == FiniteElement::VALUE;
-   if (dim==1)
   {
-      const auto W = Reshape(ir->GetWeights().Read(), Q1D);
-      const auto J = Reshape(face_geom->detJ.Read(), Q1D, NE);
-      const auto C = const_c ? Reshape(coeff.Read(), 1, 1) :
-                     Reshape(coeff.Read(), Q1D, NE);
-      auto v = Reshape(pa_data.Write(), Q1D, NE);
-      mfem::forall_2D(NE, Q1D, 1, [=] MFEM_HOST_DEVICE (int e)
+      const auto W = Reshape(ir->GetWeights().Read(), NQ);
+      const auto J = Reshape(face_geom->detJ.Read(), NQ, NE);
+      const auto C = const_c ? Reshape(coeff.Read(), 1, 1)
+                     : Reshape(coeff.Read(), NQ, NE);
+      auto v = Reshape(pa_data.Write(), NQ, NE);
+      mfem::forall(NQ, NE, [=] MFEM_HOST_DEVICE(int q, int e)
      {
-         MFEM_FOREACH_THREAD(qx,x,Q1D)
-         {
-            const real_t detJ = J(qx,e);
-            const real_t coeff = const_c ? C(0,0) : C(qx,e);
-            v(qx,e) =  W(qx) * coeff * (by_val ? detJ : 1.0/detJ);
-         }
+         const real_t detJ = J(q, e);
+         const real_t coeff = const_c ? C(0, 0) : C(q, e);
+         v(q, e) = W(q) * coeff * (by_val ? detJ : 1.0 / detJ);
      });
   }
-   else if (dim==2)
-   {
-      const auto W = Reshape(ir->GetWeights().Read(), Q1D,Q1D);
-      const auto J = Reshape(face_geom->detJ.Read(), Q1D,Q1D,NE);
-      const auto C = const_c ? Reshape(coeff.Read(), 1,1,1) :
-                     Reshape(coeff.Read(), Q1D,Q1D,NE);
-      auto v = Reshape(pa_data.Write(), Q1D,Q1D, NE);
-      mfem::forall_2D(NE, Q1D, Q1D, [=] MFEM_HOST_DEVICE (int e)
-      {
-         MFEM_FOREACH_THREAD(qx,x,Q1D)
-         {
-            MFEM_FOREACH_THREAD(qy,y,Q1D)
-            {
-               const real_t detJ = J(qx,qy,e);
-               const real_t coeff = const_c ? C(0,0,0) : C(qx,qy,e);
-               v(qx,qy,e) =  W(qx,qy) * coeff * (by_val ? detJ : 1.0/detJ);
-            }
-         }
-      });
-   }
-   else
-   {
-      MFEM_ABORT("Not supported.");
-   }
 }

 void MassIntegrator::AssembleDiagonalPA(Vector &diag)
@@ -170,6 +141,8 @@ void MassIntegrator::AssembleDiagonalPA(Vector &diag)

 void MassIntegrator::AddMultPA(const Vector &x, Vector &y) const
 {
+   MFEM_PERF_FUNCTION;
+
   if (DeviceCanUseCeed())
   {
      ceedOp->AddMult(x, y);
@@ -242,6 +242,8 @@ void DomainLFIntegrator::AssembleDevice(const FiniteElementSpace &fes,
                                        const Array<int> &markers,
                                        Vector &b)
 {
+   MFEM_PERF_FUNCTION;
+
   const FiniteElement &fe = *fes.GetTypicalFE();
   const int qorder = oa * fe.GetOrder() + ob;
   const Geometry::Type gtype = fe.GetGeomType();
@@ -346,13 +346,13 @@ private:
      template<typename T>
      T operator() (const blitz::TinyVector<T,3>& x) const
      {
-         int el_order=el->GetOrder();
-         T u1[el_order+1];
-         T u2[el_order+1];
-         T u3[el_order+1];
-         TmplPoly_1D::CalcBernstein(el_order, x[0], u1);
-         TmplPoly_1D::CalcBernstein(el_order, x[1], u2);
-         TmplPoly_1D::CalcBernstein(el_order, x[2], u3);
+         const int el_order = el->GetOrder();
+         std::vector<T> u1(el_order+1);
+         std::vector<T> u2(el_order+1);
+         std::vector<T> u3(el_order+1);
+         TmplPoly_1D::CalcBernstein(el_order, x[0], u1.data());
+         TmplPoly_1D::CalcBernstein(el_order, x[1], u2.data());
+         TmplPoly_1D::CalcBernstein(el_order, x[2], u3.data());

         const Array<int>& dof_map=el->GetDofMap();

@@ -370,17 +370,17 @@ private:
      template<typename T>
      blitz::TinyVector<T,3> grad(const blitz::TinyVector<T,3>& x) const
      {
-         int el_order=el->GetOrder();
-         T u1[el_order+1];
-         T u2[el_order+1];
-         T u3[el_order+1];
-         T d1[el_order+1];
-         T d2[el_order+1];
-         T d3[el_order+1];
+         const int el_order = el->GetOrder();
+         std::vector<T> u1(el_order+1);
+         std::vector<T> u2(el_order+1);
+         std::vector<T> u3(el_order+1);
+         std::vector<T> d1(el_order+1);
+         std::vector<T> d2(el_order+1);
+         std::vector<T> d3(el_order+1);

-         TmplPoly_1D::CalcBernstein(el_order,x[0], u1, d1);
-         TmplPoly_1D::CalcBernstein(el_order,x[1], u2, d2);
-         TmplPoly_1D::CalcBernstein(el_order,x[2], u3, d3);
+         TmplPoly_1D::CalcBernstein(el_order,x[0], u1.data(), d1.data());
+         TmplPoly_1D::CalcBernstein(el_order,x[1], u2.data(), d2.data());
+         TmplPoly_1D::CalcBernstein(el_order,x[2], u3.data(), d3.data());

         blitz::TinyVector<T,3> res(T(0.0),T(0.0),T(0.0));

@@ -415,11 +415,11 @@ private:
      template<typename T>
      T operator() (const blitz::TinyVector<T,2>& x) const
      {
-         int el_order=el->GetOrder();
-         T u1[el_order+1];
-         T u2[el_order+1];
-         TmplPoly_1D::CalcBernstein(el_order, x[0], u1);
-         TmplPoly_1D::CalcBernstein(el_order, x[1], u2);
+         const int el_order = el->GetOrder();
+         std::vector<T> u1(el_order+1);
+         std::vector<T> u2(el_order+1);
+         TmplPoly_1D::CalcBernstein(el_order, x[0], u1.data());
+         TmplPoly_1D::CalcBernstein(el_order, x[1], u2.data());

         const Array<int>& dof_map=el->GetDofMap();

@@ -437,14 +437,14 @@ private:
      template<typename T>
      blitz::TinyVector<T,2> grad(const blitz::TinyVector<T,2>& x) const
      {
-         int el_order=el->GetOrder();
-         T u1[el_order+1];
-         T u2[el_order+1];
-         T d1[el_order+1];
-         T d2[el_order+1];
+         const int el_order = el->GetOrder();
+         std::vector<T> u1(el_order+1);
+         std::vector<T> u2(el_order+1);
+         std::vector<T> d1(el_order+1);
+         std::vector<T> d2(el_order+1);

-         TmplPoly_1D::CalcBernstein(el_order,x[0], u1, d1);
-         TmplPoly_1D::CalcBernstein(el_order,x[1], u2, d2);
+         TmplPoly_1D::CalcBernstein(el_order,x[0], u1.data(), d1.data());
+         TmplPoly_1D::CalcBernstein(el_order,x[1], u2.data(), d2.data());

         blitz::TinyVector<T,2> res(T(0.0),T(0.0));

@@ -161,6 +161,7 @@ bool LinearForm::SupportsDevice() const

 void LinearForm::UseFastAssembly(bool use_fa)
 {
+   MFEM_PERF_FUNCTION;
   fast_assembly = use_fa;

   if (fast_assembly && SupportsDevice() && !ext)
@@ -171,6 +172,8 @@ void LinearForm::UseFastAssembly(bool use_fa)

 void LinearForm::Assemble()
 {
+   MFEM_PERF_FUNCTION;
+
   Array<int> vdofs;
   ElementTransformation *eltrans;
   Vector elemvect;
@@ -15,10 +15,16 @@
 namespace mfem
 {

-LinearFormExtension::LinearFormExtension(LinearForm *lf): lf(lf) { Update(); }
+LinearFormExtension::LinearFormExtension(LinearForm *lf): lf(lf)
+{
+   MFEM_PERF_FUNCTION;
+   Update();
+}

 void LinearFormExtension::Assemble()
 {
+   MFEM_PERF_FUNCTION;
+
   const FiniteElementSpace &fes = *lf->FESpace();
   MFEM_VERIFY(lf->SupportsDevice(), "Not supported.");
   MFEM_VERIFY(lf->Size() == fes.GetVSize(), "LinearForm size does not "
@@ -51,7 +57,7 @@ void LinearFormExtension::Assemble()
      {
         // scan the attributes to set the markers to 0 or 1
         const int NE = fes.GetNE();
-         const auto attr = attributes.Read();
+         const auto attr = attributes->Read();
         const auto dimk = domain_integs_marker_k->Read();
         auto markers_w = markers.Write();
         mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
@@ -93,13 +99,14 @@ void LinearFormExtension::Assemble()
      else
      {
         // scan the attributes to set the markers to 0 or 1
-         const int NBE = bdr_attributes.Size();
-         const auto attr = bdr_attributes.Read();
+         const int NBE = bdr_face_attributes->Size();
+         const auto attr = bdr_face_attributes->Read();
         const auto attr_markers = boundary_integs_marker_k->Read();
         auto markers_w = bdr_markers.Write();
-         mfem::forall(NBE, [=] MFEM_HOST_DEVICE (int e)
+         mfem::forall(NBE, [=] MFEM_HOST_DEVICE(int e)
         {
-            markers_w[e] = attr_markers[attr[e]-1] == 1;
+            markers_w[e] =
+               attr[e] > 0 ? (attr_markers[attr[e] - 1] == 1) : false;
         });
      }

@@ -112,6 +119,8 @@ void LinearFormExtension::Assemble()

 void LinearFormExtension::Update()
 {
+   MFEM_PERF_FUNCTION;
+
   const FiniteElementSpace &fes = *lf->FESpace();
   const Mesh &mesh = *fes.GetMesh();
   constexpr ElementDofOrdering ordering = ElementDofOrdering::LEXICOGRAPHIC;
@@ -125,8 +134,7 @@ void LinearFormExtension::Update()
      //markers.UseDevice(true);

      // Gather the attributes on the host from all the elements
-      attributes.SetSize(NE);
-      for (int i = 0; i < NE; ++i) { attributes[i] = mesh.GetAttribute(i); }
+      attributes = &mesh.GetElementAttributes();

      elem_restrict_lex = fes.GetElementRestriction(ordering);
      MFEM_VERIFY(elem_restrict_lex, "Element restriction not available");
@@ -136,34 +144,12 @@ void LinearFormExtension::Update()

   if (lf->boundary_integs.Size() > 0)
   {
-      const int nf_bdr = fes.GetNFbyType(FaceType::Boundary);
+      bdr_face_attributes = &mesh.GetBdrFaceAttributes();
+
+      const int nf_bdr = bdr_face_attributes->Size();
      bdr_markers.SetSize(nf_bdr);
      // bdr_markers.UseDevice(true);

-      // The face restriction will give us "face E-vectors" on the boundary that
-      // are numbered in the order of the faces of mesh. This numbering will be
-      // different than the numbering of the boundary elements. We compute
-      // mappings so that the array `bdr_attributes[i]` gives the boundary
-      // attribute of the `i`th boundary face in the mesh face order.
-      std::unordered_map<int,int> f_to_be;
-      for (int i = 0; i < mesh.GetNBE(); ++i)
-      {
-         const int f = mesh.GetBdrElementFaceIndex(i);
-         f_to_be[f] = i;
-      }
-      MFEM_VERIFY(size_t(nf_bdr) == f_to_be.size(), "Incompatible sizes");
-      bdr_attributes.SetSize(nf_bdr);
-      int f_ind = 0;
-      for (int f = 0; f < mesh.GetNumFaces(); ++f)
-      {
-         if (f_to_be.find(f) != f_to_be.end())
-         {
-            const int be = f_to_be[f];
-            bdr_attributes[f_ind] = mesh.GetBdrAttribute(be);
-            ++f_ind;
-         }
-      }
-
      bdr_restrict_lex =
         dynamic_cast<const FaceRestriction*>(
            fes.GetFaceRestriction(ordering, FaceType::Boundary,
@@ -26,7 +26,8 @@ class LinearForm;
 class LinearFormExtension
 {
   /// Attributes of all mesh elements.
-   Array<int> attributes, bdr_attributes;
+   const Array<int> *attributes; // Not owned
+   const Array<int> *bdr_face_attributes; // Not owned

   /// Temporary markers for device kernels.
   Array<int> markers, bdr_markers;
@@ -365,6 +365,8 @@ FiniteElementSpace &LORBase::GetFESpace() const

 void LORBase::AssembleSystem(BilinearForm &a_ho, const Array<int> &ess_dofs)
 {
+   MFEM_PERF_FUNCTION;
+
   A.Clear();
   delete a;
   if (BatchedLORAssembly::FormIsSupported(a_ho))
@@ -360,6 +360,8 @@ void BatchedLORAssembly::FillJAndData(SparseMatrix &A) const

 void BatchedLORAssembly::SparseIJToCSR(OperatorHandle &A) const
 {
+   MFEM_PERF_FUNCTION;
+
   const int nvdof = fes_ho.GetVSize();

   // If A contains an existing SparseMatrix, reuse it (and try to reuse its
@@ -417,6 +419,8 @@ static void Assemble_(LOR_KERNEL &kernel, int dim, int sdim, int order)
 template <typename LOR_KERNEL>
 void BatchedLORAssembly::AssemblyKernel(BilinearForm &a)
 {
+   MFEM_PERF_FUNCTION;
+
   LOR_KERNEL kernel(a, fes_ho, X_vert, sparse_ij, sparse_mapping);

   const int dim = fes_ho.GetMesh()->Dimension();
@@ -184,6 +184,8 @@ void BatchedLOR_H1::Assemble2D()
 template <int ORDER>
 void BatchedLOR_H1::Assemble3D()
 {
+   MFEM_PERF_FUNCTION;
+
   const int nel_ho = fes_ho.GetNE();
   static constexpr int nv = 8;
   static constexpr int dim = 3;
@@ -10,6 +10,7 @@
 // CONTRIBUTING.md for details.

 #include "multigrid.hpp"
+#include "../general/annotation.hpp"

 namespace mfem
 {
@@ -135,6 +136,8 @@ void MultigridBase::Mult(const Vector& x, Vector& y) const
 void MultigridBase::ArrayMult(const Array<const Vector*>& X_,
                              Array<Vector*>& Y_) const
 {
+   MFEM_PERF_FUNCTION;
+
   MFEM_ASSERT(operators.Size() > 0,
               "Multigrid solver does not have operators set!");
   MFEM_ASSERT(X_.Size() == Y_.Size(),
@@ -158,6 +161,8 @@ void MultigridBase::ArrayMult(const Array<const Vector*>& X_,

 void MultigridBase::SmoothingStep(int level, bool zero, bool transpose) const
 {
+   MFEM_PERF_FUNCTION;
+
   // y = y + S (x - A y) or y = y + S^T (x - A y)

   // Note: 'zero' == true means that Y(level,*) are not initialized and we
@@ -201,6 +206,8 @@ void MultigridBase::SmoothingStep(int level, bool zero, bool transpose) const

 void MultigridBase::CoarseSolve(bool zero) const
 {
+   MFEM_PERF_FUNCTION;
+
   // See the comment about iterative mode in SmoothingStep()
   coarse_solver->iterative_mode = false;

@@ -124,6 +124,8 @@ void ParBilinearForm::pAllocMat()
 void ParBilinearForm::ParallelRAP(SparseMatrix &loc_A, OperatorHandle &A,
                                  bool steal_loc_A)
 {
+   MFEM_PERF_FUNCTION;
+
   ParFiniteElementSpace &pfespace = *ParFESpace();

   // Create a block diagonal parallel matrix
@@ -0,0 +1,591 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "pderefmat_op.hpp"
+
+#ifdef MFEM_USE_MPI
+
+#include "fes_kernels.hpp"
+/// \cond DO_NOT_DOCUMENT
+namespace mfem
+{
+namespace internal
+{
+template <Ordering::Type Order, bool Atomic>
+static void ParDerefMultKernelImpl(const ParDerefineMatrixOp &op,
+                                   const Vector &x, Vector &y)
+{
+   // pack sends
+   if (op.xghost_send.Size())
+   {
+      auto src = x.Read();
+      auto idcs = op.send_permutations.Read();
+      auto dst = Device::GetGPUAwareMPI() ? op.xghost_send.Write()
+                 : op.xghost_send.HostWrite();
+      auto vdims = op.fespace->GetVDim();
+      auto sptr = op.send_segment_idcs.Read();
+      auto lptr = op.send_segments.Read();
+      auto old_ndofs = x.Size() / vdims;
+
+      forall(op.send_permutations.Size(), [=] MFEM_HOST_DEVICE(int i)
+      {
+         int seg = sptr[i];
+         int width = lptr[seg + 1] - lptr[seg];
+         auto tdst = dst + i + lptr[seg] * vdims;
+         int sign = 1;
+         int col = idcs[i];
+         if (col < 0)
+         {
+            sign = -1;
+            col = -1 - col;
+         }
+         for (int vdim = 0; vdim < vdims; ++vdim)
+         {
+            tdst[vdim * width] =
+               sign
+               * src[Order == Ordering::byNODES ? (col + vdim * old_ndofs)
+                           : (col * vdims + vdim)];
+         }
+      });
+      // TODO: is this needed so we can send the packed data correctly?
+      // unclear for GPU-aware MPI, definitely required otherwise
+      MFEM_DEVICE_SYNC;
+   }
+   // initialize off-diagonal receive and send
+   op.requests.clear();
+   if (op.xghost_recv.Size())
+   {
+      auto vdims = op.fespace->GetVDim();
+      auto rcv = Device::GetGPUAwareMPI() ? op.xghost_recv.Write()
+                 : op.xghost_recv.HostWrite();
+      for (int i = 0; i < op.recv_ranks.Size(); ++i)
+      {
+         op.requests.emplace_back();
+         MPI_Irecv(rcv + op.recv_segments[i] * vdims,
+                   (op.recv_segments[i + 1] - op.recv_segments[i]) * vdims,
+                   MPITypeMap<real_t>::mpi_type, op.recv_ranks[i],
+                   MessageTag::DEREFINEMENT_MATRIX_CONSTRUCTION_DATA,
+                   op.fespace->GetComm(), &op.requests.back());
+      }
+   }
+   if (op.xghost_send.Size())
+   {
+      auto vdims = op.fespace->GetVDim();
+      // only is a GPU mem ptr if GPU-aware MPI is enabled
+      auto dst = Device::GetGPUAwareMPI() ? op.xghost_send.Write()
+                 : op.xghost_send.HostWrite();
+      for (int i = 0; i < op.send_ranks.Size(); ++i)
+      {
+         op.requests.emplace_back();
+         MPI_Isend(dst + op.send_segments[i] * vdims,
+                   (op.send_segments[i + 1] - op.send_segments[i]) * vdims,
+                   MPITypeMap<real_t>::mpi_type, op.send_ranks[i],
+                   MessageTag::DEREFINEMENT_MATRIX_CONSTRUCTION_DATA,
+                   op.fespace->GetComm(), &op.requests.back());
+      }
+   }
+   {
+      // diagonal
+      DerefineMatrixOpMultFunctor<Order, Atomic, true> func;
+      func.xptr = x.Read();
+      y.UseDevice();
+      y = 0.;
+      func.yptr = y.ReadWrite();
+      func.bsptr = op.block_storage.Read();
+      func.boptr = op.block_offsets.Read();
+      func.brptr = op.block_row_idcs_offsets.Read();
+      func.bcptr = op.block_col_idcs_offsets.Read();
+      func.rptr = op.row_idcs.Read();
+      func.cptr = op.col_idcs.Read();
+      func.vdims = op.fespace->GetVDim();
+      func.nblocks = op.block_offsets.Size();
+      func.width = op.Width() / func.vdims;
+      func.height = op.Height() / func.vdims;
+      func.Run(op.max_rows);
+   }
+   // wait for comm to finish, if any
+   if (op.requests.size())
+   {
+      MPI_Waitall(op.requests.size(), op.requests.data(), MPI_STATUSES_IGNORE);
+      if (op.xghost_recv.Size())
+      {
+         // off-diagonal kernel
+         DerefineMatrixOpMultFunctor<Order, Atomic, false> func;
+         // directly read from host-pinned memory if not using GPU-aware MPI
+         func.xptr = Device::GetGPUAwareMPI() ? op.xghost_recv.Read()
+                     : op.xghost_recv.HostRead();
+         func.yptr = y.ReadWrite();
+         func.bsptr = op.block_storage.Read();
+         func.boptr = op.off_diag_block_offsets.Read();
+         func.brptr = op.block_off_diag_row_idcs_offsets.Read();
+         func.rsptr = op.recv_segment_idcs.Read();
+         func.segptr = op.recv_segments.Read();
+         func.coptr = op.block_off_diag_col_offsets.Read();
+         func.bwptr = op.block_off_diag_widths.Read();
+         func.rptr = op.row_off_diag_idcs.Read();
+         func.vdims = op.fespace->GetVDim();
+         func.nblocks = op.off_diag_block_offsets.Size();
+         func.width = op.xghost_recv.Size() / func.vdims;
+         func.height = op.Height() / func.vdims;
+         func.Run(op.max_rows);
+      }
+   }
+}
+} // namespace internal
+
+template <Ordering::Type Order, bool Atomic>
+ParDerefineMatrixOp::MultKernelType ParDerefineMatrixOp::MultKernel::Kernel()
+{
+   return internal::ParDerefMultKernelImpl<Order, Atomic>;
+}
+
+ParDerefineMatrixOp::MultKernelType
+ParDerefineMatrixOp::MultKernel::Fallback(Ordering::Type, bool)
+{
+   MFEM_ABORT("invalid MultKernel parameters");
+}
+
+ParDerefineMatrixOp::Kernels::Kernels()
+{
+   MultKernel::Specialization<Ordering::byNODES, false>::Add();
+   MultKernel::Specialization<Ordering::byVDIM, false>::Add();
+   MultKernel::Specialization<Ordering::byNODES, true>::Add();
+   MultKernel::Specialization<Ordering::byVDIM, true>::Add();
+}
+
+void ParDerefineMatrixOp::Mult(const Vector &x, Vector &y) const
+{
+   const bool is_dg = fespace->FEColl()->GetContType()
+                      == FiniteElementCollection::DISCONTINUOUS;
+   // DG needs atomic summation
+   MultKernel::Run(fespace->GetOrdering(), is_dg, *this, x, y);
+   // use this to prevent xghost* from being re-purposed for subsequent Mult
+   // calls
+   MFEM_DEVICE_SYNC;
+}
+
+ParDerefineMatrixOp::ParDerefineMatrixOp(ParFiniteElementSpace &fespace_,
+                                         int old_ndofs,
+                                         const Table *old_elem_dof,
+                                         const Table *old_elem_fos)
+   : Operator(fespace_.GetVSize(), old_ndofs * fespace_.GetVDim()),
+     fespace(&fespace_)
+{
+   static Kernels kernels;
+   constexpr int max_team_size = 256;
+
+   const int NRanks = fespace->GetNRanks();
+
+   const int nrk = HYPRE_AssumedPartitionCheck() ? 2 : NRanks;
+
+   MFEM_VERIFY(fespace->Nonconforming(),
+               "Not implemented for conforming meshes.");
+   MFEM_VERIFY(fespace->old_dof_offsets[nrk],
+               "Missing previous (finer) space.");
+
+   const int MyRank = fespace->GetMyRank();
+   ParNCMesh *old_pncmesh = fespace->GetParMesh()->pncmesh;
+   const CoarseFineTransformations &dtrans =
+      old_pncmesh->GetDerefinementTransforms();
+   const Array<int> &old_ranks = old_pncmesh->GetDerefineOldRanks();
+
+   const bool is_dg = fespace->FEColl()->GetContType()
+                      == FiniteElementCollection::DISCONTINUOUS;
+   DenseMatrix localRVO; // for variable-order only
+
+   DenseTensor localR[Geometry::NumGeom];
+   int diag_rows = 0;
+   int off_diag_rows = 0;
+   int diag_cols = 0;
+
+   auto get_ldofs = [&](int k) -> int
+   {
+      const Embedding &emb = dtrans.embeddings[k];
+      if (fespace->IsVariableOrder())
+      {
+         const FiniteElement *fe = fespace->GetFE(emb.parent);
+         return fe->GetDof();
+      }
+      else
+      {
+         Geometry::Type geom =
+         fespace->GetParMesh()->GetElementBaseGeometry(emb.parent);
+         return fespace->FEColl()->FiniteElementForGeometry(geom)->GetDof();
+      }
+   };
+   Array<int> dofs, old_dofs;
+   max_rows = 1;
+   // first pass:
+   // - determine memory block lengths
+   // - identify dofs in x we need to send/receive
+   // don't need to send the indices, fine rank will re-arrange and sign
+   // change x before transmitting the ghost data
+
+   // key: coarse rank to send to
+   // value: old dofs to send (with sign)
+   std::map<int, std::vector<int>> to_send;
+   // key: fine rank
+   // value: indices into dtrans.embeddings
+   std::map<int, std::vector<int>> od_ks;
+   // key: fine rank
+   // value: recv segment length
+   std::map<int, int> od_seg_lens;
+   int send_len = 0;
+   int recv_len = 0;
+   // size of block_storage, if fespace->IsVariableOrder()
+   // otherwise unused
+   int total_size = 0;
+   int num_diagonal_blocks = 0;
+   int num_offdiagonal_blocks = 0;
+   for (int k = 0; k < dtrans.embeddings.Size(); ++k)
+   {
+      const Embedding &emb = dtrans.embeddings[k];
+      int fine_rank = old_ranks[k];
+      int coarse_rank = (emb.parent < 0) ? (-1 - emb.parent)
+                        : old_pncmesh->ElementRank(emb.parent);
+      if (coarse_rank != MyRank && fine_rank == MyRank)
+      {
+         // this rank needs to send data in x to course_rank
+         old_elem_dof->GetRow(k, old_dofs);
+         auto &tmp = to_send[coarse_rank];
+         send_len += old_dofs.Size();
+         for (int i = 0; i < old_dofs.Size(); ++i)
+         {
+            tmp.emplace_back(old_dofs[i]);
+         }
+      }
+      else if (coarse_rank == MyRank && fine_rank != MyRank)
+      {
+         // this rank needs to receive data in x from fine_rank
+         MFEM_ASSERT(emb.parent >= 0, "");
+         auto ldofs = get_ldofs(k);
+         off_diag_rows += ldofs;
+         recv_len += ldofs;
+         od_ks[fine_rank].emplace_back(k);
+         od_seg_lens[fine_rank] += ldofs;
+         ++num_offdiagonal_blocks;
+         if (fespace->IsVariableOrder())
+         {
+            total_size += ldofs * ldofs;
+         }
+      }
+      else if (coarse_rank == MyRank && fine_rank == MyRank)
+      {
+         MFEM_ASSERT(emb.parent >= 0, "");
+         // diagonal
+         ++num_diagonal_blocks;
+         auto ldofs = get_ldofs(k);
+         diag_rows += ldofs;
+         diag_cols += ldofs;
+         if (fespace->IsVariableOrder())
+         {
+            total_size += ldofs * ldofs;
+         }
+      }
+   }
+   send_segments.SetSize(to_send.size() + 1);
+   send_segments.HostWrite();
+   send_ranks.SetSize(to_send.size());
+   send_ranks.HostWrite();
+   {
+      int idx = 0;
+      send_segments[0] = 0;
+      for (auto &tmp : to_send)
+      {
+         send_ranks[idx] = tmp.first;
+         send_segments[idx + 1] = send_segments[idx] + tmp.second.size();
+         ++idx;
+      }
+   }
+   recv_segment_idcs.SetSize(off_diag_rows);
+   recv_segment_idcs.HostWrite();
+   recv_segments.SetSize(od_ks.size() + 1);
+   recv_segments.HostWrite();
+   recv_ranks.SetSize(od_ks.size());
+   recv_ranks.HostWrite();
+
+   // set sizes
+   row_idcs.SetSize(diag_rows);
+   row_idcs.HostWrite();
+   row_off_diag_idcs.SetSize(off_diag_rows);
+   row_off_diag_idcs.HostWrite();
+   col_idcs.SetSize(diag_cols);
+   col_idcs.HostWrite();
+   block_row_idcs_offsets.SetSize(num_diagonal_blocks + 1);
+   block_row_idcs_offsets.HostWrite();
+   block_col_idcs_offsets.SetSize(num_diagonal_blocks + 1);
+   block_col_idcs_offsets.HostWrite();
+   block_off_diag_row_idcs_offsets.SetSize(num_offdiagonal_blocks + 1);
+   block_off_diag_row_idcs_offsets.HostWrite();
+   block_off_diag_col_offsets.SetSize(num_offdiagonal_blocks);
+   block_off_diag_col_offsets.HostWrite();
+   block_off_diag_widths.SetSize(num_offdiagonal_blocks);
+   block_off_diag_widths.HostWrite();
+   pack_col_idcs.SetSize(send_len);
+   // memory manager doesn't appear to have a graceful fallback for
+   // HOST_PINNED if not built with CUDA or HIP
+#if defined(MFEM_USE_CUDA) or defined(MFEM_USE_HIP)
+   xghost_send.SetSize(send_len * fespace->GetVDim(),
+                       Device::GetGPUAwareMPI() ? MemoryType::DEFAULT
+                       : MemoryType::HOST_PINNED);
+   xghost_recv.SetSize(recv_len * fespace->GetVDim(),
+                       Device::GetGPUAwareMPI() ? MemoryType::DEFAULT
+                       : MemoryType::HOST_PINNED);
+#else
+   xghost_send.SetSize(send_len * fespace->GetVDim());
+   xghost_recv.SetSize(recv_len * fespace->GetVDim());
+#endif
+   send_permutations.SetSize(send_len);
+   send_segment_idcs.SetSize(send_len);
+   block_offsets.SetSize(num_diagonal_blocks);
+   block_offsets.HostWrite();
+   off_diag_block_offsets.SetSize(num_offdiagonal_blocks);
+   off_diag_block_offsets.HostWrite();
+   int geom_offsets[Geometry::NumGeom];
+   real_t *bs_ptr;
+
+   if (fespace->IsVariableOrder())
+   {
+      block_storage.SetSize(total_size);
+      bs_ptr = block_storage.HostWrite();
+      // compute block data later
+   }
+   else
+   {
+      // compression scheme:
+      // block_offsets is the start of each block, potentially repeated
+      // only need to store localR for used shapes
+      Mesh::GeometryList elem_geoms(*fespace->GetMesh());
+
+      int size = 0;
+      for (int i = 0; i < elem_geoms.Size(); ++i)
+      {
+         fespace->GetLocalDerefinementMatrices(elem_geoms[i],
+                                               localR[elem_geoms[i]]);
+         geom_offsets[elem_geoms[i]] = size;
+         size += localR[elem_geoms[i]].TotalSize();
+      }
+      block_storage.SetSize(size);
+      bs_ptr = block_storage.HostWrite();
+      // copy blocks into block_storage
+      for (int i = 0; i < elem_geoms.Size(); ++i)
+      {
+         std::copy(localR[elem_geoms[i]].Data(),
+                   localR[elem_geoms[i]].Data()
+                   + localR[elem_geoms[i]].TotalSize(),
+                   bs_ptr);
+         bs_ptr += localR[elem_geoms[i]].TotalSize();
+      }
+   }
+
+   // second pass:
+   // - initialize buffers
+
+   {
+      auto ptr = send_permutations.HostWrite();
+      auto ptr2 = send_segment_idcs.HostWrite();
+      int i = 0;
+      for (auto &v : to_send)
+      {
+         ptr = std::copy(v.second.begin(), v.second.end(), ptr);
+         for (size_t idx = 0; idx < v.second.size(); ++idx)
+         {
+            *ptr2 = i;
+            ++ptr2;
+         }
+         ++i;
+      }
+   }
+
+   block_row_idcs_offsets[0] = 0;
+   block_col_idcs_offsets[0] = 0;
+   block_off_diag_row_idcs_offsets[0] = 0;
+   Array<int> mark(fespace->GetNDofs());
+   mark = 0;
+   {
+      int idx = 0;
+      recv_segments[0] = 0;
+      for (auto &v : od_seg_lens)
+      {
+         recv_ranks[idx] = v.first;
+         recv_segments[idx + 1] = recv_segments[idx] + v.second;
+         ++idx;
+      }
+   }
+   // key: index into dtrans.embeddings
+   // value: off-diagonal block offset, od_ridx, seg id
+   std::unordered_map<int, std::array<int, 3>> ks_map;
+   {
+      int od_ridx = 0;
+      int seg_id = 0;
+      for (auto &v1 : od_ks)
+      {
+         for (auto k : v1.second)
+         {
+            auto &tmp = ks_map[k];
+            tmp[0] = ks_map.size() - 1;
+            tmp[1] = od_ridx;
+            tmp[2] = seg_id;
+            od_ridx += get_ldofs(k);
+         }
+         ++seg_id;
+      }
+   }
+   int diag_idx = 0;
+   int var_offset = 0;
+   int ridx = 0;
+   int cidx = 0;
+   // can't break this up into separate diagonals/off-diagonals loops because
+   // of mark
+   for (int k = 0; k < dtrans.embeddings.Size(); ++k)
+   {
+      const Embedding &emb = dtrans.embeddings[k];
+      if (emb.parent < 0)
+      {
+         continue;
+      }
+      int fine_rank = old_ranks[k];
+      int coarse_rank = (emb.parent < 0) ? (-1 - emb.parent)
+                        : old_pncmesh->ElementRank(emb.parent);
+      if (coarse_rank == MyRank)
+      {
+         // either diagonal or off-diagonal
+         Geometry::Type geom =
+            fespace->GetMesh()->GetElementBaseGeometry(emb.parent);
+         if (fespace->IsVariableOrder())
+         {
+            const FiniteElement *fe = fespace->GetFE(emb.parent);
+            const DenseTensor &pmats = dtrans.point_matrices[geom];
+            const int ldof = fe->GetDof();
+
+            IsoparametricTransformation isotr;
+            isotr.SetIdentityTransformation(geom);
+
+            localRVO.SetSize(ldof, ldof);
+            isotr.SetPointMat(pmats(emb.matrix));
+            // Local restriction is size ldofxldof assuming that the parent
+            // and child are of same polynomial order.
+            fe->GetLocalRestriction(isotr, localRVO);
+            // copy block
+            auto s = localRVO.Height() * localRVO.Width();
+            std::copy(localRVO.Data(), localRVO.Data() + s, bs_ptr);
+            bs_ptr += s;
+         }
+         DenseMatrix &lR =
+            fespace->IsVariableOrder() ? localRVO : localR[geom](emb.matrix);
+         max_rows = std::max(lR.Height(), max_rows);
+         auto size = lR.Height() * lR.Width();
+         fespace->elem_dof->GetRow(emb.parent, dofs);
+         if (fine_rank == MyRank)
+         {
+            // diagonal
+            old_elem_dof->GetRow(k, old_dofs);
+            MFEM_VERIFY(old_dofs.Size() == dofs.Size(),
+                        "Parent and child must have same #dofs.");
+            block_row_idcs_offsets[diag_idx + 1] =
+               block_row_idcs_offsets[diag_idx] + lR.Height();
+            block_col_idcs_offsets[diag_idx + 1] =
+               block_col_idcs_offsets[diag_idx] + lR.Width();
+
+            if (fespace->IsVariableOrder())
+            {
+               block_offsets[diag_idx] = var_offset;
+               var_offset += size;
+            }
+            else
+            {
+               block_offsets[diag_idx] = geom_offsets[geom] + size * emb.matrix;
+            }
+            for (int i = 0; i < lR.Height(); ++i, ++ridx)
+            {
+               if (!std::isfinite(lR(i, 0)))
+               {
+                  row_idcs[ridx] = INT_MAX;
+                  continue;
+               }
+               int r = dofs[i];
+               int m = (r >= 0) ? r : (-1 - r);
+               if (is_dg || !mark[m])
+               {
+                  row_idcs[ridx] = r;
+                  mark[m] = 1;
+               }
+               else
+               {
+                  row_idcs[ridx] = INT_MAX;
+               }
+            }
+            for (int i = 0; i < lR.Width(); ++i, ++cidx)
+            {
+               col_idcs[cidx] = old_dofs[i];
+            }
+            ++diag_idx;
+         }
+         else
+         {
+            // off-diagonal
+            auto &tmp = ks_map.at(k);
+            auto od_idx = tmp[0];
+            auto od_ridx = tmp[1];
+            block_off_diag_row_idcs_offsets[od_idx + 1] =
+               block_off_diag_row_idcs_offsets[od_idx] + lR.Height();
+            block_off_diag_col_offsets[od_idx] = od_ridx;
+            block_off_diag_widths[od_idx] = lR.Width();
+            recv_segment_idcs[od_idx] = tmp[2];
+
+            if (fespace->IsVariableOrder())
+            {
+               off_diag_block_offsets[od_idx] = var_offset;
+               var_offset += size;
+            }
+            else
+            {
+               off_diag_block_offsets[od_idx] =
+                  geom_offsets[geom] + size * emb.matrix;
+            }
+            for (int i = 0; i < lR.Height(); ++i, ++od_ridx)
+            {
+               if (!std::isfinite(lR(i, 0)))
+               {
+                  row_off_diag_idcs[od_ridx] = INT_MAX;
+                  continue;
+               }
+               int r = dofs[i];
+               int m = (r >= 0) ? r : (-1 - r);
+               if (is_dg || !mark[m])
+               {
+                  row_off_diag_idcs[od_ridx] = r;
+                  mark[m] = 1;
+               }
+               else
+               {
+                  row_off_diag_idcs[od_ridx] = INT_MAX;
+               }
+            }
+            ++od_idx;
+         }
+      }
+   }
+   // if not using GPU, set max_rows/max_cols to zero
+   if (Device::Allows(Backend::DEVICE_MASK))
+   {
+      max_rows = std::min(max_rows, max_team_size);
+   }
+   else
+   {
+      max_rows = 1;
+   }
+   requests.reserve(recv_ranks.Size() + send_ranks.Size());
+}
+} // namespace mfem
+/// \endcond DO_NOT_DOCUMENT
+
+#endif
@@ -0,0 +1,111 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_PDEREFMAT_OP
+#define MFEM_PDEREFMAT_OP
+
+#include "../config/config.hpp"
+
+#ifdef MFEM_USE_MPI
+
+#include "pfespace.hpp"
+
+#include "kernel_dispatch.hpp"
+
+#include <vector>
+
+namespace mfem
+{
+
+/// \cond DO_NOT_DOCUMENT
+struct ParDerefineMatrixOp : public Operator
+{
+   ParFiniteElementSpace *fespace;
+   /// offsets into block_storage for diagonal
+   Array<int> block_offsets;
+   /// offsets into row_idcs for diagonal
+   Array<int> block_row_idcs_offsets;
+   /// offsets into col_idcs for diagonal
+   Array<int> block_col_idcs_offsets;
+
+   /// offsets into block_storage for off-diagonal
+   Array<int> off_diag_block_offsets;
+   /// offsets into row_idcs for off-diagonal
+   Array<int> block_off_diag_row_idcs_offsets;
+   Array<int> block_off_diag_col_offsets;
+   Array<int> block_off_diag_widths;
+   /// mapping for row dofs, INT_MAX indicates the block row should be ignored.
+   /// negative means the row data should be negated.
+   /// only for diagonal blocks
+   Array<int> row_idcs;
+   /// mapping for col dofs, negative means the col data should be negated.
+   /// only for diagonal blocks
+   Array<int> col_idcs;
+
+   Array<int> pack_col_idcs;
+
+   /// mapping for row dofs, INT_MAX indicates the block row should be ignored.
+   /// negative means the row data should be negated.
+   /// only for off-diagonal blocks
+   Array<int> row_off_diag_idcs;
+   /// dense block matrices which can be reused to construct the full matrix
+   /// operation. These are stored contiguously and blocks have no restrictions
+   /// on shape (can be rectangle and differ from block to block).
+   /// This is only for the diagonal block.
+   Vector block_storage;
+   /// maximum height of any block in block_storage for GPU
+   /// parallelization, or 1 for CPU runs.
+   int max_rows;
+
+   /// quasi Ordering::byNODES, broken into sections by ranks we need to send
+   /// the data to
+   mutable Vector xghost_send;
+   /// quasi Ordering::byNODES, broken into sections by ranks we received
+   /// the data from
+   mutable Vector xghost_recv;
+   /// maps off-diagonal k to segment
+   Array<int> recv_segment_idcs;
+   /// cumulative count of dofs which will be received from other ranks
+   Array<int> recv_segments;
+   /// Source rank of each recv segment
+   Array<int> recv_ranks;
+   /// What send segment each entry in send_permutations corresponds to
+   Array<int> send_segment_idcs;
+   /// cumulative count of dofs which will be sent to other ranks
+   Array<int> send_segments;
+   /// Destination rank of each send segment
+   Array<int> send_ranks;
+   /// how to permute/sign change values from our local x to send to other ranks
+   Array<int> send_permutations;
+   /// internal buffer for MPI requests
+   mutable std::vector<MPI_Request> requests;
+
+   using MultKernelType = void (*)(const ParDerefineMatrixOp &, const Vector &,
+                                   Vector &);
+   /// template args: ordering, atomic
+   MFEM_REGISTER_KERNELS(MultKernel, MultKernelType, (Ordering::Type, bool));
+
+   struct Kernels
+   {
+      Kernels();
+   };
+
+   void Mult(const Vector &x, Vector &y) const;
+
+   ParDerefineMatrixOp(ParFiniteElementSpace &fespace_, int old_ndofs,
+                       const Table *old_elem_dof, const Table *old_elem_fos);
+};
+/// \endcond DO_NOT_DOCUMENT
+} // namespace mfem
+
+#endif
+
+#endif
@@ -22,12 +22,13 @@
 #include "../mesh/mesh_headers.hpp"
 #include "../general/binaryio.hpp"

+#include "pderefmat_op.hpp"
+
 #include <limits>
 #include <list>

 namespace mfem
 {
-
 ParFiniteElementSpace::ParFiniteElementSpace(
   const ParFiniteElementSpace &orig, ParMesh *pmesh,
   const FiniteElementCollection *fec)
@@ -62,9 +63,11 @@ ParFiniteElementSpace::ParFiniteElementSpace(

 ParFiniteElementSpace::ParFiniteElementSpace(
   ParMesh *pm, const FiniteElementCollection *f, int dim, int ordering)
-   : FiniteElementSpace(pm, f, dim, ordering)
+   : FiniteElementSpace((MFEM_PERF_BEGIN(_MFEM_FUNC_NAME), pm),
+                        f, dim, ordering)
 {
   ParInit(pm);
+   MFEM_PERF_END(_MFEM_FUNC_NAME);
 }

 ParFiniteElementSpace::ParFiniteElementSpace(
@@ -91,6 +94,7 @@ ParNURBSExtension *ParFiniteElementSpace::MakeLocalNURBSext(

 void ParFiniteElementSpace::ParInit(ParMesh *pm)
 {
+   MFEM_PERF_FUNCTION;
   pmesh = pm;
   pncmesh = nullptr;

@@ -180,6 +184,7 @@ void ParFiniteElementSpace::CommunicateGhostOrder()

 void ParFiniteElementSpace::Construct()
 {
+   MFEM_PERF_FUNCTION;
   if (NURBSext)
   {
      ConstructTrueNURBSDofs();
@@ -838,6 +843,8 @@ void ParFiniteElementSpace::Build_Dof_TrueDof_Matrix() const // matrix P

   if (P) { return; }

+   MFEM_PERF_FUNCTION;
+
   if (!nd_strias)
   {
      // Safe to assume 1-1 correspondence between shared dofs
@@ -1423,6 +1430,7 @@ const Operator *ParFiniteElementSpace::GetRestrictionOperator() const
      if (NRanks == 1)
      {
         R_transpose.reset(new IdentityOperator(GetTrueVSize()));
+         Rconf = new IdentityOperator(GetTrueVSize());
      }
      else
      {
@@ -1435,8 +1443,8 @@ const Operator *ParFiniteElementSpace::GetRestrictionOperator() const
            R_transpose.reset(
               new DeviceConformingProlongationOperator(*this, true));
         }
+         Rconf = new TransposeOperator(*R_transpose);
      }
-      Rconf = new TransposeOperator(*R_transpose);
      return Rconf;
   }
   else
@@ -4487,13 +4495,6 @@ ParFiniteElementSpace::RebalanceMatrix(int old_ndofs,
   return M;
 }

-
-struct DerefDofMessage
-{
-   std::vector<HYPRE_BigInt> dofs;
-   MPI_Request request;
-};
-
 HypreParMatrix*
 ParFiniteElementSpace::ParallelDerefinementMatrix(int old_ndofs,
                                                  const Table* old_elem_dof,
@@ -4536,7 +4537,13 @@ ParFiniteElementSpace::ParallelDerefinementMatrix(int old_ndofs,
      old_pncmesh->GetDerefinementTransforms();
   const Array<int> &old_ranks = old_pncmesh->GetDerefineOldRanks();

-   std::map<int, DerefDofMessage> messages;
+   // key: other rank
+   // value: send or recieve buffer
+   std::map<int, std::vector<HYPRE_BigInt>> to_send;
+   std::map<int, std::vector<HYPRE_BigInt>> to_recv;
+   // key: index into dtrans.embeddings
+   // value: [start, stop]
+   std::unordered_map<int, std::array<size_t, 2>> recv_messages;

   HYPRE_BigInt old_offset = HYPRE_AssumedPartitionCheck()
                             ? old_dof_offsets[0] : old_dof_offsets[MyRank];
@@ -4556,30 +4563,46 @@ ParFiniteElementSpace::ParallelDerefinementMatrix(int old_ndofs,
         old_elem_dof->GetRow(k, dofs);
         DofsToVDofs(dofs, old_ndofs);

-         DerefDofMessage &msg = messages[k];
-         msg.dofs.resize(dofs.Size());
+         std::vector<HYPRE_BigInt>& send_buf = to_send[coarse_rank];
+         auto pos = send_buf.size();
+         send_buf.resize(pos + dofs.Size());
         for (int i = 0; i < dofs.Size(); i++)
         {
-            msg.dofs[i] = old_offset + dofs[i];
+            send_buf[pos + i] = old_offset + dofs[i];
         }
-
-         MPI_Isend(&msg.dofs[0], static_cast<int>(msg.dofs.size()), HYPRE_MPI_BIG_INT,
-                   coarse_rank, 291, MyComm, &msg.request);
      }
      else if (coarse_rank == MyRank && fine_rank != MyRank)
      {
         MFEM_ASSERT(emb.parent >= 0, "");
         Geometry::Type geom = mesh->GetElementBaseGeometry(emb.parent);

-         DerefDofMessage &msg = messages[k];
-         msg.dofs.resize(ldof[geom]*vdim);
-
-         MPI_Irecv(&msg.dofs[0], ldof[geom]*vdim, HYPRE_MPI_BIG_INT,
-                   fine_rank, 291, MyComm, &msg.request);
+         std::vector<HYPRE_BigInt>& recv_buf = to_recv[fine_rank];
+         auto& msg = recv_messages[k];
+         msg[0] = recv_buf.size();
+         recv_buf.resize(recv_buf.size() + ldof[geom] * vdim);
+         msg[1] = recv_buf.size();
      }
-      // TODO: coalesce Isends/Irecvs to the same rank. Typically, on uniform
-      // derefinement, there should be just one send to MyRank-1 and one recv
-      // from MyRank+1
+   }
+
+   // assume embedding orders are consistent (i.e. what we expect to receive
+   // first from a given rank is sent first, etc.)
+   std::vector<MPI_Request> requests;
+   requests.reserve(to_send.size() + to_recv.size());
+   // enqueue recvs
+   for (auto &v : to_recv)
+   {
+      requests.emplace_back();
+      MPI_Irecv(v.second.data(), v.second.size(), HYPRE_MPI_BIG_INT, v.first,
+                MessageTag::DEREFINEMENT_MATRIX_CONSTRUCTION_DATA, MyComm,
+                &requests.back());
+   }
+   // enqueue sends
+   for (auto &v : to_send)
+   {
+      requests.emplace_back();
+      MPI_Isend(v.second.data(), v.second.size(), HYPRE_MPI_BIG_INT, v.first,
+                MessageTag::DEREFINEMENT_MATRIX_CONSTRUCTION_DATA, MyComm,
+                &requests.back());
   }

   DenseTensor localR[Geometry::NumGeom];
@@ -4637,10 +4660,7 @@ ParFiniteElementSpace::ParallelDerefinementMatrix(int old_ndofs,
   diag->Finalize();

   // wait for all sends/receives to complete
-   for (auto it = messages.begin(); it != messages.end(); ++it)
-   {
-      MPI_Wait(&it->second.request, MPI_STATUS_IGNORE);
-   }
+   MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);

   // create the off-diagonal part of the derefinement matrix
   SparseMatrix *offd = new SparseMatrix(ndofs*vdim, 1);
@@ -4661,13 +4681,14 @@ ParFiniteElementSpace::ParallelDerefinementMatrix(int old_ndofs,

         elem_dof->GetRow(emb.parent, dofs);

-         DerefDofMessage &msg = messages[k];
-         MFEM_ASSERT(msg.dofs.size(), "");
+         auto& odofs = to_recv.at(fine_rank);
+         auto &msg = recv_messages[k];
+         MFEM_ASSERT(msg[1] > msg[0], "");

         for (int vd = 0; vd < vdim; vd++)
         {
            MFEM_ASSERT(ldof[geom], "");
-            HYPRE_BigInt* remote_dofs = &msg.dofs[vd*ldof[geom]];
+            HYPRE_BigInt *remote_dofs = odofs.data() + msg[0] + vd * ldof[geom];

            for (int i = 0; i < lR.Height(); i++)
            {
@@ -4694,7 +4715,6 @@ ParFiniteElementSpace::ParallelDerefinementMatrix(int old_ndofs,
      }
   }

-   messages.clear();
   offd->Finalize(0);
   offd->SetWidth(static_cast<int>(col_map.size()));

@@ -4946,8 +4966,13 @@ void ParFiniteElementSpace::Update(bool want_transform)

         case Mesh::DEREFINE:
         {
+#if 0
            Th.Reset(ParallelDerefinementMatrix(old_ndofs, old_elem_dof,
                                                old_elem_fos));
+#else
+            Th.Reset(new ParDerefineMatrixOp(*this, old_ndofs, old_elem_dof,
+                                             old_elem_fos));
+#endif
            if (Nonconforming())
            {
               Th.SetOperatorOwner(false);
@@ -24,9 +24,12 @@
 namespace mfem
 {

+struct ParDerefineMatrixOp;
+
 /// Abstract parallel finite element space.
 class ParFiniteElementSpace : public FiniteElementSpace
 {
+   friend struct ParDerefineMatrixOp;
 private:
   /// MPI data.
   MPI_Comm MyComm;
@@ -45,6 +45,8 @@ void ParLinearForm::MakeRef(ParFiniteElementSpace *pf, Vector &v, int v_offset)

 void ParLinearForm::Assemble()
 {
+   MFEM_PERF_FUNCTION;
+
   LinearForm::Assemble();

   if (interior_face_integs.Size())
@@ -214,16 +214,16 @@ public:
 inline void QuadratureFunction::GetValues(
   int idx, Vector &values)
 {
-   const int s_offset = qspace->offsets[idx];
-   const int sl_size = qspace->offsets[idx+1] - s_offset;
+   const int s_offset = qspace->Offset(idx);
+   const int sl_size = qspace->Offset(idx + 1) - s_offset;
   values.MakeRef(*this, vdim*s_offset, vdim*sl_size);
 }

 inline void QuadratureFunction::GetValues(
   int idx, Vector &values) const
 {
-   const int s_offset = qspace->offsets[idx];
-   const int sl_size = qspace->offsets[idx+1] - s_offset;
+   const int s_offset = qspace->Offset(idx);
+   const int sl_size = qspace->Offset(idx + 1) - s_offset;
   values.SetSize(vdim*sl_size);
   values.HostWrite();
   const real_t *q = HostRead() + vdim*s_offset;
@@ -236,14 +236,14 @@ inline void QuadratureFunction::GetValues(
 inline void QuadratureFunction::GetValues(
   int idx, const int ip_num, Vector &values)
 {
-   const int s_offset = qspace->offsets[idx] * vdim + ip_num * vdim;
+   const int s_offset = qspace->Offset(idx) * vdim + ip_num * vdim;
   values.MakeRef(*this, s_offset, vdim);
 }

 inline void QuadratureFunction::GetValues(
   int idx, const int ip_num, Vector &values) const
 {
-   const int s_offset = qspace->offsets[idx] * vdim + ip_num * vdim;
+   const int s_offset = qspace->Offset(idx) * vdim + ip_num * vdim;
   values.SetSize(vdim);
   values.HostWrite();
   const real_t *q = HostRead() + s_offset;
@@ -256,8 +256,8 @@ inline void QuadratureFunction::GetValues(
 inline void QuadratureFunction::GetValues(
   int idx, DenseMatrix &values)
 {
-   const int s_offset = qspace->offsets[idx];
-   const int sl_size = qspace->offsets[idx+1] - s_offset;
+   const int s_offset = qspace->Offset(idx);
+   const int sl_size = qspace->Offset(idx + 1) - s_offset;
   // Make the values matrix memory an alias of the quadrature function memory
   Memory<real_t> &values_mem = values.GetMemory();
   values_mem.Delete();
@@ -268,8 +268,8 @@ inline void QuadratureFunction::GetValues(
 inline void QuadratureFunction::GetValues(
   int idx, DenseMatrix &values) const
 {
-   const int s_offset = qspace->offsets[idx];
-   const int sl_size = qspace->offsets[idx+1] - s_offset;
+   const int s_offset = qspace->Offset(idx);
+   const int sl_size = qspace->Offset(idx + 1) - s_offset;
   values.SetSize(vdim, sl_size);
   values.HostWrite();
   const real_t *q = HostRead() + vdim*s_offset;
@@ -96,17 +96,33 @@ void QuadratureSpaceBase::Integrate(VectorCoefficient &coeff,

 void QuadratureSpace::ConstructOffsets()
 {
+   MFEM_PERF_FUNCTION;
   const int num_elem = mesh.GetNE();
-   offsets.SetSize(num_elem + 1);
-   int offset = 0;
-   for (int i = 0; i < num_elem; i++)
+   ne = num_elem;
+
+   if (mesh.GetNumGeometries(mesh.Dimension()) == 1)
   {
-      offsets[i] = offset;
-      int geom = mesh.GetElementBaseGeometry(i);
-      MFEM_ASSERT(int_rule[geom] != NULL, "Missing integration rule.");
-      offset += int_rule[geom]->GetNPoints();
+      Array<Geometry::Type> geoms;
+      mesh.GetGeometries(mesh.Dimension(), geoms);
+      offsets.SetSize(1);
+      offsets.HostWrite();
+      offsets[0] = int_rule[geoms[0]]->GetNPoints();
+      size = num_elem * offsets[0];
+   }
+   else
+   {
+      offsets.SetSize(num_elem + 1);
+      int offset = 0;
+      for (int i = 0; i < num_elem; i++)
+      {
+         offsets[i] = offset;
+         const Geometry::Type geom = mesh.GetElementBaseGeometry(i);
+         MFEM_ASSERT(int_rule[geom] != NULL, "Missing integration rule.");
+         offset += int_rule[geom]->GetNPoints();
+      }
+      offsets[num_elem] = offset;
+      size = offsets.Last();
   }
-   offsets[num_elem] = size = offset;
 }

 void QuadratureSpace::Construct()
@@ -188,6 +204,7 @@ void FaceQuadratureSpace::ConstructOffsets()
 {
   face_indices.SetSize(num_faces);
   offsets.SetSize(num_faces + 1);
+   ne = num_faces;
   int offset = 0;
   int f_idx = 0;
   for (int i = 0; i < mesh.GetNumFacesWithGhost(); i++)
@@ -30,13 +30,17 @@ protected:
   Mesh &mesh; ///< The underlying mesh.
   int order; ///< The order of integration rule.
   int size; ///< Total number of quadrature points.
+   int ne; ///< Actual number of entities
   mutable Vector weights; ///< Integration weights.
   mutable long nodes_sequence = 0; ///< Nodes counter for cache invalidation.

-   /// @brief Entity quadrature point offset array, of size num_entities + 1.
+   /// @brief Entity quadrature point offset array.
+   ///
+   /// Supports a constant compression scheme for meshes which have a single
+   /// geometry type. When compressed, will have a single value. The true offset
+   /// can be computed as i * offsets[0], where i is the entity index. Otherwise
+   /// has size num_entities + 1.
   ///
-   /// The quadrature point values for entity i are stored in the indices between
-   /// offsets[i] and offsets[i+1].
   Array<int> offsets;
   /// The quadrature rules used for each geometry type.
   const IntegrationRule *int_rule[Geometry::NumGeom];
@@ -59,6 +63,24 @@ protected:
   void ConstructWeights() const;

 public:
+   /// @brief Gets the offset for a given entity @a idx.
+   ///
+   /// The quadrature point values for entity i are stored in the indices
+   /// between Offset(i) and Offset(i+1)
+   int Offset(int idx) const
+   {
+      return (offsets.Size() == 1) ? (idx * offsets[0]) : offsets[idx];
+   }
+
+   /// @brief Entity quadrature point offset array.
+   ///
+   /// Supports a constant compression scheme for meshes which have a single
+   /// geometry type. When compressed, will have a single value. The true offset
+   /// can be computed as i * offsets[0], where i is the entity index. Otherwise
+   /// has size num_entities + 1.
+   ///
+   const Array<int> &Offsets() const { return offsets; }
+
   /// Return the total number of quadrature points.
   int GetSize() const { return size; }

@@ -66,7 +88,7 @@ public:
   int GetOrder() const { return order; }

   /// Return the number of entities.
-   int GetNE() const { return offsets.Size() - 1; }
+   int GetNE() const { return ne; }

   /// Returns the mesh.
   inline Mesh *GetMesh() const { return &mesh; }
@@ -503,6 +503,7 @@ void QuadratureInterpolator::Mult(const Vector &e_vec,
                                  Vector &q_der,
                                  Vector &q_det) const
 {
+   MFEM_PERF_FUNCTION;
   using namespace internal::quadrature_interpolator;

   const int ne = fespace->GetNE();
@@ -25,7 +25,7 @@ namespace mfem

 ElementRestriction::ElementRestriction(const FiniteElementSpace &f,
                                       ElementDofOrdering e_ordering)
-   : fes(f),
+   : fes((MFEM_PERF_BEGIN(_MFEM_FUNC_NAME), f)),
     ne(fes.GetNE()),
     vdim(fes.GetVDim()),
     byvdim(fes.GetOrdering() == Ordering::byVDIM),
@@ -104,10 +104,13 @@ ElementRestriction::ElementRestriction(const FiniteElementSpace &f,
      offsets[i] = offsets[i - 1];
   }
   offsets[0] = 0;
+   MFEM_PERF_END(_MFEM_FUNC_NAME);
 }

 void ElementRestriction::Mult(const Vector& x, Vector& y) const
 {
+   MFEM_PERF_FUNCTION;
+
   // Assumes all elements have the same number of dofs
   const int nd = dof;
   const int vd = vdim;
@@ -152,6 +155,8 @@ void ElementRestriction::AbsMult(const Vector& x, Vector& y) const
 template <bool ADD>
 void ElementRestriction::TAddMultTranspose(const Vector& x, Vector& y) const
 {
+   MFEM_PERF_FUNCTION;
+
   // Assumes all elements have the same number of dofs
   const int nd = dof;
   const int vd = vdim;
@@ -2328,6 +2328,8 @@ static void SmemProlongation3D(const int NE,
                               const Vector& localL, Vector& localH,
                               const Array<real_t> &b, const Vector& mask)
 {
+   MFEM_PERF_FUNCTION;
+
   auto u_lo = Reshape(localL.Read(), DLO, DLO, DLO, NE);
   auto u_hi = Reshape(localH.Write(), DHI, DHI, DHI, NE);
   auto d_b = Reshape(b.Read(), DHI, DLO);
@@ -2523,6 +2525,8 @@ static void SmemProlongationTranspose3D(
   const int NE, const Vector& localH, Vector& localL,
   const Array<real_t>& bt, const Vector& mask)
 {
+   MFEM_PERF_FUNCTION;
+
   auto u_h = Reshape(localH.Read(), DHI, DHI, DHI, NE);
   auto u_l = Reshape(localL.Write(), DLO, DLO, DLO, NE);
   auto d_bt = Reshape(bt.Read(), DLO, DHI);
@@ -2662,6 +2666,7 @@ void ProlongationTranspose3D(const int NE, const int D1D, const int Q1D,
 void TensorProductPRefinementTransferOperator::Mult(const Vector& x,
                                                    Vector& y) const
 {
+   MFEM_PERF_FUNCTION;
   using namespace TransferKernels;

   if (lFESpace.GetMesh()->GetNE() == 0)
@@ -2704,6 +2709,7 @@ void TensorProductPRefinementTransferOperator::Mult(const Vector& x,
 void TensorProductPRefinementTransferOperator::MultTranspose(const Vector& x,
                                                             Vector& y) const
 {
+   MFEM_PERF_FUNCTION;
   using namespace TransferKernels;

   if (lFESpace.GetMesh()->GetNE() == 0)
@@ -2754,20 +2760,20 @@ TrueTransferOperator::TrueTransferOperator(const FiniteElementSpace& lFESpace_,

   P = lFESpace.GetProlongationMatrix();
   R = hFESpace.IsVariableOrder() ? hFESpace.GetHpRestrictionMatrix() :
-       hFESpace.GetRestrictionMatrix();
+       hFESpace.GetRestrictionOperator();

   // P and R can be both null
   // P can be null and R not null
   // If P is not null it is assumed that R is not null as well
   if (P) { MFEM_VERIFY(R, "Both P and R have to be not NULL") }

-   if (P)
+   if (!IsIdentityProlongation(P))
   {
      tmpL.SetSize(lFESpace_.GetVSize());
      tmpH.SetSize(hFESpace_.GetVSize());
   }
   // P can be null and R not null
-   else if (R)
+   else if (!IsIdentityProlongation(R))
   {
      tmpH.SetSize(hFESpace_.GetVSize());
   }
@@ -2780,13 +2786,14 @@ TrueTransferOperator::~TrueTransferOperator()

 void TrueTransferOperator::Mult(const Vector& x, Vector& y) const
 {
-   if (P)
+   MFEM_PERF_FUNCTION;
+   if (!IsIdentityProlongation(P))
   {
      P->Mult(x, tmpL);
      localTransferOperator->Mult(tmpL, tmpH);
      R->Mult(tmpH, y);
   }
-   else if (R)
+   else if (!IsIdentityProlongation(R))
   {
      localTransferOperator->Mult(x, tmpH);
      R->Mult(tmpH, y);
@@ -2799,13 +2806,14 @@ void TrueTransferOperator::Mult(const Vector& x, Vector& y) const

 void TrueTransferOperator::MultTranspose(const Vector& x, Vector& y) const
 {
-   if (P)
+   MFEM_PERF_FUNCTION;
+   if (!IsIdentityProlongation(P))
   {
      R->MultTranspose(x, tmpH);
      localTransferOperator->MultTranspose(tmpH, tmpL);
      P->MultTranspose(tmpL, y);
   }
-   else if (R)
+   else if (!IsIdentityProlongation(R))
   {
      R->MultTranspose(x, tmpH);
      localTransferOperator->MultTranspose(tmpH, y);
@@ -639,7 +639,7 @@ private:
   const FiniteElementSpace& lFESpace;
   const FiniteElementSpace& hFESpace;
   const Operator * P = nullptr;
-   const SparseMatrix * R = nullptr;
+   const Operator * R = nullptr;
   TransferOperator* localTransferOperator;
   mutable Vector tmpL;
   mutable Vector tmpH;
@@ -98,7 +98,7 @@ struct ScopeAnnotation
 #else // #ifdef MFEM_USE_CALIPER

 #define MFEM_PERF_FUNCTION
-#define MFEM_PERF_BEGIN(s)
+#define MFEM_PERF_BEGIN(s) ((void)(0))
 #define MFEM_PERF_END(s)
 #define MFEM_PERF_SCOPE(name)

@@ -14,7 +14,7 @@

 #include "../config/config.hpp"

-#ifdef MFEM_USE_CUDA
+#if defined(MFEM_USE_CUDA) && defined(__CUDACC__)
 #include <cusparse.h>
 #include <library_types.h>
 #include <cuda_runtime.h>
@@ -22,7 +22,7 @@
 #endif
 #include "cuda.hpp"

-#ifdef MFEM_USE_HIP
+#if defined(MFEM_USE_HIP) && defined(__HIP__)
 #include <hip/hip_runtime.h>
 #endif
 #include "hip.hpp"
@@ -43,7 +43,7 @@
 #endif
 #endif

-#if !(defined(MFEM_USE_CUDA) || defined(MFEM_USE_HIP))
+#if !defined(MFEM_USE_CUDA_OR_HIP)
 #define MFEM_DEVICE
 #define MFEM_HOST
 #define MFEM_LAMBDA
@@ -55,7 +55,7 @@
 #endif

 #if !((defined(MFEM_USE_CUDA) && defined(__CUDA_ARCH__)) || \
-      (defined(MFEM_USE_HIP)  && defined(__HIP_DEVICE_COMPILE__)))
+      (defined(MFEM_USE_HIP) && defined(__HIP_DEVICE_COMPILE__)))
 #define MFEM_SHARED
 #define MFEM_SYNC_THREAD
 #define MFEM_BLOCK_ID(k) 0
@@ -66,7 +66,7 @@
 #endif

 // 'double' and 'float' atomicAdd implementation for previous versions of CUDA
-#if defined(MFEM_USE_CUDA) && defined(__CUDA_ARCH__) && __CUDA_ARCH__ < 600
+#if defined(MFEM_USE_CUDA) && defined(__CUDA_ARCH__) && (__CUDA_ARCH__ < 600)
 MFEM_DEVICE inline mfem::real_t atomicAdd(mfem::real_t *add, mfem::real_t val)
 {
   unsigned long long int *ptr = (unsigned long long int *) add;
@@ -94,7 +94,7 @@ template <typename T>
 MFEM_HOST_DEVICE T AtomicAdd(T &add, const T val)
 {
 #if ((defined(MFEM_USE_CUDA) && defined(__CUDA_ARCH__)) || \
-     (defined(MFEM_USE_HIP)  && defined(__HIP_DEVICE_COMPILE__)))
+     (defined(MFEM_USE_HIP) && defined(__HIP_DEVICE_COMPILE__)))
   return atomicAdd(&add,val);
 #else
   T old = add;
@@ -347,13 +347,23 @@ void GroupTopology::Swap(GroupTopology &other)
   mfem::Swap(group_mgroup, other.group_mgroup);
 }

+/// \cond DO_NOT_DOCUMENT
 // Initialize the static mpi_type for the specializations of MPITypeMap:
+const MPI_Datatype MPITypeMap<bool>::mpi_type = MFEM_MPI_CXX_BOOL;
+const MPI_Datatype MPITypeMap<char>::mpi_type = MPI_CHAR;
+const MPI_Datatype MPITypeMap<unsigned char>::mpi_type = MPI_UNSIGNED_CHAR;
+const MPI_Datatype MPITypeMap<short>::mpi_type = MPI_SHORT;
+const MPI_Datatype MPITypeMap<unsigned short>::mpi_type = MPI_UNSIGNED_SHORT;
 const MPI_Datatype MPITypeMap<int>::mpi_type = MPI_INT;
+const MPI_Datatype MPITypeMap<unsigned int>::mpi_type = MPI_UNSIGNED;
+const MPI_Datatype MPITypeMap<long>::mpi_type = MPI_LONG;
+const MPI_Datatype MPITypeMap<unsigned long>::mpi_type = MPI_UNSIGNED_LONG;
+const MPI_Datatype MPITypeMap<long long>::mpi_type = MPI_LONG_LONG;
+const MPI_Datatype MPITypeMap<unsigned long long>::mpi_type =
+   MPI_UNSIGNED_LONG_LONG;
 const MPI_Datatype MPITypeMap<float>::mpi_type = MPI_FLOAT;
 const MPI_Datatype MPITypeMap<double>::mpi_type = MPI_DOUBLE;
-const MPI_Datatype MPITypeMap<int64_t>::mpi_type = MPI_INT64_T;
-const MPI_Datatype MPITypeMap<uint64_t>::mpi_type = MPI_UINT64_T;
-
+/// \endcond DO_NOT_DOCUMENT

 GroupCommunicator::GroupCommunicator(const GroupTopology &gt, Mode m)
   : gtopo(gt), mode(m)
@@ -421,16 +421,24 @@ public:
   ~GroupCommunicator();
 };

+/// General MPI message tags used by MFEM
+enum MessageTag
+{
+   DEREFINEMENT_MATRIX_CONSTRUCTION_DATA =
+      291, /// ParFiniteElementSpace ParallelDerefinementMatrix and
+   /// ParDerefineMatrixOp
+};
+
 enum VarMessageTag
 {
-   NEIGHBOR_ELEMENT_RANK_VM,  ///< NeighborElementRankMessage
-   NEIGHBOR_ORDER_VM,         ///< NeighborOrderMessage
-   NEIGHBOR_DEREFINEMENT_VM,  ///< NeighborDerefinementMessage
-   NEIGHBOR_REFINEMENT_VM,    ///< NeighborRefinementMessage
-   NEIGHBOR_PREFINEMENT_VM,   ///< NeighborPRefinementMessage
-   NEIGHBOR_ROW_VM,           ///< NeighborRowMessage
-   REBALANCE_VM,              ///< RebalanceMessage
-   REBALANCE_DOF_VM           ///< RebalanceDofMessage
+   NEIGHBOR_ELEMENT_RANK_VM, ///< NeighborElementRankMessage
+   NEIGHBOR_ORDER_VM,        ///< NeighborOrderMessage
+   NEIGHBOR_DEREFINEMENT_VM, ///< NeighborDerefinementMessage
+   NEIGHBOR_REFINEMENT_VM,   ///< NeighborRefinementMessage
+   NEIGHBOR_PREFINEMENT_VM,  ///< NeighborPRefinementMessage
+   NEIGHBOR_ROW_VM,          ///< NeighborRowMessage
+   REBALANCE_VM,             ///< RebalanceMessage
+   REBALANCE_DOF_VM,         ///< RebalanceDofMessage
 };

 /// \brief Variable-length MPI message containing unspecific binary data.
@@ -595,10 +603,50 @@ protected:
 template <typename Type> struct MPITypeMap;

 // Specializations of MPITypeMap; mpi_type initialized in communication.cpp:
+template<> struct MPITypeMap<bool>
+{
+   static MFEM_EXPORT const MPI_Datatype mpi_type;
+};
+template<> struct MPITypeMap<char>
+{
+   static MFEM_EXPORT const MPI_Datatype mpi_type;
+};
+template<> struct MPITypeMap<unsigned char>
+{
+   static MFEM_EXPORT const MPI_Datatype mpi_type;
+};
+template<> struct MPITypeMap<short>
+{
+   static MFEM_EXPORT const MPI_Datatype mpi_type;
+};
+template<> struct MPITypeMap<unsigned short>
+{
+   static MFEM_EXPORT const MPI_Datatype mpi_type;
+};
 template<> struct MPITypeMap<int>
 {
   static MFEM_EXPORT const MPI_Datatype mpi_type;
 };
+template<> struct MPITypeMap<unsigned int>
+{
+   static MFEM_EXPORT const MPI_Datatype mpi_type;
+};
+template<> struct MPITypeMap<long>
+{
+   static MFEM_EXPORT const MPI_Datatype mpi_type;
+};
+template<> struct MPITypeMap<unsigned long>
+{
+   static MFEM_EXPORT const MPI_Datatype mpi_type;
+};
+template<> struct MPITypeMap<long long>
+{
+   static MFEM_EXPORT const MPI_Datatype mpi_type;
+};
+template<> struct MPITypeMap<unsigned long long>
+{
+   static MFEM_EXPORT const MPI_Datatype mpi_type;
+};
 template<> struct MPITypeMap<double>
 {
   static MFEM_EXPORT const MPI_Datatype mpi_type;
@@ -607,15 +655,6 @@ template<> struct MPITypeMap<float>
 {
   static MFEM_EXPORT const MPI_Datatype mpi_type;
 };
-template<> struct MPITypeMap<int64_t>
-{
-   static MFEM_EXPORT const MPI_Datatype mpi_type;
-};
-template<> struct MPITypeMap<uint64_t>
-{
-   static MFEM_EXPORT const MPI_Datatype mpi_type;
-};
-

 /** Reorder MPI ranks to follow the Z-curve within the physical machine topology
    (provided that functions to query physical node coordinates are available).
@@ -24,7 +24,7 @@ void mfem_cuda_error(cudaError_t err, const char *expr, const char *func,
                     const char *file, int line)
 {
   mfem::err << "\n\nCUDA error: (" << expr << ") failed with error:\n --> "
-             << cudaGetErrorString(err)
+             << cudaGetErrorString(err) << " [code: " << (int)err << ']'
             << "\n ... in function: " << func
             << "\n ... in file: " << file << ':' << line << '\n';
   mfem_error();
@@ -18,7 +18,7 @@
 // CUDA block size used by MFEM.
 #define MFEM_CUDA_BLOCKS 256

-#ifdef MFEM_USE_CUDA
+#if defined(MFEM_USE_CUDA) && defined(__CUDACC__)
 #define MFEM_USE_CUDA_OR_HIP
 #define MFEM_DEVICE __device__
 #define MFEM_HOST __host__
@@ -37,10 +37,9 @@
                              __FILE__, __LINE__);                             \
    }                                                                          \
  } while (0)
-#endif // MFEM_USE_CUDA

 // Define the MFEM inner threading macros
-#if defined(MFEM_USE_CUDA) && defined(__CUDA_ARCH__)
+#if defined(__CUDA_ARCH__)
 #define MFEM_SHARED __shared__
 #define MFEM_SYNC_THREAD __syncthreads()
 #define MFEM_BLOCK_ID(k) blockIdx.k
@@ -48,12 +47,13 @@
 #define MFEM_THREAD_SIZE(k) blockDim.k
 #define MFEM_FOREACH_THREAD(i,k,N) for(int i=threadIdx.k; i<N; i+=blockDim.k)
 #define MFEM_FOREACH_THREAD_DIRECT(i,k,N) if(const int i=threadIdx.k; i<N)
-#endif
+#endif // defined(__CUDA_ARCH__)
+#endif // defined(MFEM_USE_CUDA) && defined(__CUDACC__)

 namespace mfem
 {

-#ifdef MFEM_USE_CUDA
+#if defined(MFEM_USE_CUDA) && defined(__CUDACC__)
 // Function used by the macro MFEM_GPU_CHECK.
 void mfem_cuda_error(cudaError_t err, const char *expr, const char *func,
                     const char *file, int line);
@@ -631,7 +631,7 @@ void Device::Setup(const std::string &device_option, const int device_id)
   if (Allows(Backend::DEBUG_DEVICE)) { ngpu = 1; }
 }

-MemoryType Device::QueryMemoryType(void *ptr)
+MemoryType Device::QueryMemoryType(const void* ptr)
 {
   // from HYPRE's hypre_GetPointerLocation
   MemoryType res = MemoryType::HOST;
@@ -297,9 +297,9 @@ public:
   /// Get the status of GPU-aware MPI flag.
   static bool GetGPUAwareMPI() { return Get().mpi_gpu_aware; }

-   /** @brief Query the device driver for what memory type a given @a ptr is
-       allocated with. */
-   static MemoryType QueryMemoryType(void *ptr);
+   /** Query the device driver for what memory type a given @a ptr is allocated
+    * with. */
+   static MemoryType QueryMemoryType(const void* ptr);

   /** @brief The number of hardware compute units/streaming multiprocessors
       available on a given compute device @a device_id. */
@@ -176,7 +176,7 @@ __device__ void abort_msg(T & msg)
      printf(__VA_ARGS__);     \
      asm("trap;");            \
   }
-#elif defined(MFEM_USE_HIP)
+#elif defined(__HIP_DEVICE_COMPILE__)
 #define MFEM_ABORT_KERNEL(...) \
   {                           \
      printf(__VA_ARGS__);     \
@@ -158,8 +158,8 @@ private:
 #define MFEM_PRAGMA(X) _Pragma(#X)

 // MFEM_UNROLL pragma macro that can be used inside MFEM_FORALL macros.
-#if defined(MFEM_USE_CUDA) && defined(__CUDA_ARCH__)
-#ifdef __NVCC__
+#if defined(MFEM_USE_CUDA) && defined(__CUDA_ARCH__) // Clang cuda or nvcc
+#ifdef __NVCC__ // nvcc specifically
 #define MFEM_UNROLL(N) MFEM_PRAGMA(unroll(N))
 #else // Assuming Clang CUDA
 #define MFEM_UNROLL(N) MFEM_PRAGMA(unroll N)
@@ -169,12 +169,12 @@ private:
 #endif

 // MFEM_GPU_FORALL: "parallel for" executed with CUDA or HIP based on the MFEM
-// build-time configuration (MFEM_USE_CUDA or MFEM_USE_HIP). If neither CUDA nor
-// HIP is enabled, this macro is a no-op.
-#if defined(MFEM_USE_CUDA)
+// build-time configuration (MFEM_USE_CUDA or MFEM_USE_HIP), and if compiling
+// with CUDA/HIP language. Otherwise, this macro is a no-op.
+#if defined(MFEM_USE_CUDA) && defined(__CUDACC__)
 #define MFEM_GPU_FORALL(i, N,...) CuWrap1D(N, [=] MFEM_DEVICE      \
                                       (int i) {__VA_ARGS__})
-#elif defined(MFEM_USE_HIP)
+#elif defined(MFEM_USE_HIP) && defined(__HIP__)
 #define MFEM_GPU_FORALL(i, N,...) HipWrap1D(N, [=] MFEM_DEVICE     \
                                        (int i) {__VA_ARGS__})
 #else
@@ -225,6 +225,52 @@ void OmpWrap(const int N, HBODY &&h_body)
 #endif
 }

+template <typename HBODY>
+void OmpWrap2D(const int Nx, const int Ny, HBODY &&h_body)
+{
+#ifdef MFEM_USE_OPENMP
+   // requires OpenMP 3.1
+   #pragma omp parallel for collapse(2)
+   for (int j = 0; j < Ny; j++)
+   {
+      for (int i = 0; i < Nx; i++)
+      {
+         h_body(i, j);
+      }
+   }
+#else
+   MFEM_CONTRACT_VAR(Nx);
+   MFEM_CONTRACT_VAR(Ny);
+   MFEM_CONTRACT_VAR(h_body);
+   MFEM_ABORT("OpenMP requested for MFEM but OpenMP is not enabled!");
+#endif
+}
+
+template <typename HBODY>
+void OmpWrap3D(const int Nx, const int Ny, const int Nz, HBODY &&h_body)
+{
+#ifdef MFEM_USE_OPENMP
+   // requires OpenMP 3.1
+   #pragma omp parallel for collapse(3)
+   for (int k = 0; k < Nz; k++)
+   {
+      for (int j = 0; j < Ny; j++)
+      {
+         for (int i = 0; i < Nx; i++)
+         {
+            h_body(i, j, k);
+         }
+      }
+   }
+#else
+   MFEM_CONTRACT_VAR(Nx);
+   MFEM_CONTRACT_VAR(Ny);
+   MFEM_CONTRACT_VAR(Nz);
+   MFEM_CONTRACT_VAR(h_body);
+   MFEM_ABORT("OpenMP requested for MFEM but OpenMP is not enabled!");
+#endif
+}
+

 /// RAJA Cuda and Hip backends
 #if defined(MFEM_USE_RAJA) && defined(RAJA_ENABLE_CUDA)
@@ -454,6 +500,42 @@ void RajaOmpWrap(const int N, HBODY &&h_body)
   RAJA::forall<RAJA::omp_parallel_for_exec>(RAJA::RangeSegment(0,N), h_body);
 }

+template <typename HBODY>
+void RajaOmpWrap2D(const int Nx, const int Ny, HBODY &&h_body)
+{
+   using omp_launch_policy = RAJA::LaunchPolicy<RAJA::omp_launch_t>;
+   using global_thread_xy = RAJA::LoopPolicy<RAJA::omp_for_exec>;
+   RAJA::RangeSegment xrange(0, Nx);
+   RAJA::RangeSegment yrange(0, Ny);
+   RAJA::launch<omp_launch_policy>(RAJA::ExecPlace::HOST, RAJA::LaunchParams(),
+                                   [=](RAJA::LaunchContext ctx)
+   {
+      // contiguous in x
+      RAJA::expt::loop<global_thread_xy>(ctx, xrange, yrange, [&](int i, int j)
+      {
+         h_body(i, j);
+      });
+   });
+}
+
+template <typename HBODY>
+void RajaOmpWrap3D(const int Nx, const int Ny, const int Nz, HBODY &&h_body)
+{
+   using omp_launch_policy = RAJA::LaunchPolicy<RAJA::omp_launch_t>;
+   using global_thread_xyz = RAJA::LoopPolicy<RAJA::omp_for_exec>;
+   RAJA::RangeSegment xrange(0, Nx);
+   RAJA::RangeSegment yrange(0, Ny);
+   RAJA::RangeSegment zrange(0, Nz);
+   RAJA::launch<omp_launch_policy>(RAJA::ExecPlace::HOST, RAJA::LaunchParams(),
+                                   [=](RAJA::LaunchContext ctx)
+   {
+      // contiguous in x
+      RAJA::expt::loop<global_thread_xyz>(ctx, xrange, yrange, zrange,
+                                          [&](int i, int j, int k)
+      { h_body(i, j, k); });
+   });
+}
+
 #endif


@@ -481,7 +563,7 @@ void RajaSeqWrap(const int N, HBODY &&h_body)


 /// CUDA backend
-#ifdef MFEM_USE_CUDA
+#if defined(MFEM_USE_CUDA) && defined(__CUDACC__)

 template <typename BODY> __global__ static
 void CuKernel1D(const int N, BODY body)
@@ -573,11 +655,11 @@ struct CuWrap<3>
   }
 };

-#endif // MFEM_USE_CUDA
+#endif // defined(MFEM_USE_CUDA) && defined(__CUDACC__)


 /// HIP backend
-#ifdef MFEM_USE_HIP
+#if defined(MFEM_USE_HIP) && defined(__HIP__)

 template <typename BODY> __global__ static
 void HipKernel1D(const int N, BODY body)
@@ -668,7 +750,7 @@ struct HipWrap<3>
   }
 };

-#endif // MFEM_USE_HIP
+#endif // defined(MFEM_USE_HIP) && defined(__HIP__)


 /// The forall kernel body wrapper
@@ -701,7 +783,7 @@ inline void ForallWrap(const bool use_dev, const int N,
   }
 #endif

-#ifdef MFEM_USE_CUDA
+#if defined(MFEM_USE_CUDA) && defined(__CUDACC__)
   // If Backend::CUDA is allowed, use it
   if (Device::Allows(Backend::CUDA))
   {
@@ -709,7 +791,7 @@ inline void ForallWrap(const bool use_dev, const int N,
   }
 #endif

-#ifdef MFEM_USE_HIP
+#if defined(MFEM_USE_HIP) && defined(__HIP__)
   // If Backend::HIP is allowed, use it
   if (Device::Allows(Backend::HIP))
   {
@@ -753,6 +835,83 @@ inline void ForallWrap(const bool use_dev, const int N, lambda &&body,
 template<typename lambda>
 inline void forall(int N, lambda &&body) { ForallWrap<1>(true, N, body); }

+template<typename lambda>
+inline void forall(int Nx, int Ny, lambda &&body)
+{
+   if (Device::Allows(Backend::DEVICE_MASK))
+   {
+      forall(Nx * Ny, [=] MFEM_HOST_DEVICE(int idx)
+      {
+         int j = idx / Nx;
+         int i = idx % Nx;
+         body(i, j);
+      });
+   }
+#if defined(MFEM_USE_RAJA) && defined(RAJA_ENABLE_OPENMP)
+   else if (Device::Allows(Backend::RAJA_OMP))
+   {
+      return RajaOmpWrap2D(Nx, Ny, body);
+   }
+#endif
+#ifdef MFEM_USE_OPENMP
+   else if (Device::Allows(Backend::OMP))
+   {
+      return OmpWrap2D(Nx, Ny, body);
+   }
+#endif
+   else
+   {
+      for (int j = 0; j < Ny; ++j)
+      {
+         for (int i = 0; i < Nx; ++i)
+         {
+            body(i, j);
+         }
+      }
+   }
+}
+
+template<typename lambda>
+inline void forall(int Nx, int Ny, int Nz, lambda &&body)
+{
+   if (Device::Allows(Backend::DEVICE_MASK))
+   {
+      forall(Nx * Ny * Nz, [=] MFEM_HOST_DEVICE(int idx)
+      {
+         int i = idx % Nx;
+         int j = idx / Nx;
+         int k = j / Ny;
+         j = j % Ny;
+         body(i, j, k);
+      });
+   }
+#if defined(MFEM_USE_RAJA) && defined(RAJA_ENABLE_OPENMP)
+   else if (Device::Allows(Backend::RAJA_OMP))
+   {
+      return RajaOmpWrap3D(Nx, Ny, Nz, body);
+   }
+#endif
+#ifdef MFEM_USE_OPENMP
+   else if (Device::Allows(Backend::OMP))
+   {
+      return OmpWrap3D(Nx, Ny, Nz, body);
+   }
+#endif
+   else
+   {
+      for (int k = 0; k < Nz; ++k)
+      {
+         for (int j = 0; j < Ny; ++j)
+         {
+            for (int i = 0; i < Nx; ++i)
+            {
+               body(i, j, k);
+            }
+         }
+      }
+   }
+}
+
 template<typename lambda>
 inline void forall_switch(bool use_dev, int N, lambda &&body)
 {
@@ -24,7 +24,7 @@ void mfem_hip_error(hipError_t err, const char *expr, const char *func,
                    const char *file, int line)
 {
   mfem::err << "\n\nHIP error: (" << expr << ") failed with error:\n --> "
-             << hipGetErrorString(err)
+             << hipGetErrorString(err) << " [code: " << (int)err << ']'
             << "\n ... in function: " << func
             << "\n ... in file: " << file << ':' << line << '\n';
   mfem_error();
@@ -18,7 +18,7 @@
 // HIP block size used by MFEM.
 #define MFEM_HIP_BLOCKS 256

-#ifdef MFEM_USE_HIP
+#if defined(MFEM_USE_HIP) && defined(__HIP__)
 #define MFEM_USE_CUDA_OR_HIP
 #define MFEM_DEVICE __device__
 #define MFEM_HOST __host__
@@ -37,10 +37,9 @@
                             __FILE__, __LINE__);                              \
    }                                                                          \
  } while (0)
-#endif // MFEM_USE_HIP

 // Define the MFEM inner threading macros
-#if defined(MFEM_USE_HIP) && defined(__HIP_DEVICE_COMPILE__)
+#if defined(__HIP_DEVICE_COMPILE__)
 #define MFEM_SHARED __shared__
 #define MFEM_SYNC_THREAD __syncthreads()
 #define MFEM_BLOCK_ID(k) hipBlockIdx_ ##k
@@ -50,7 +49,8 @@
   for(int i=hipThreadIdx_ ##k; i<N; i+=hipBlockDim_ ##k)
 #define MFEM_FOREACH_THREAD_DIRECT(i,k,N) \
   if(const int i=hipThreadIdx_ ##k; i<N)
-#endif
+#endif // defined(__HIP_DEVICE_COMPILE__)
+#endif // defined(MFEM_USE_HIP) && defined(__HIP__)

 namespace mfem
 {
@@ -513,10 +513,7 @@ public:
   void *HtoD(void *dst, const void *src, size_t bytes) override
   { return HipMemcpyHtoD(dst, src, bytes); }
   void *DtoD(void* dst, const void* src, size_t bytes) override
-   // Unlike cudaMemcpy(DtoD), hipMemcpy(DtoD) causes a host-side synchronization so
-   // instead we use hipMemcpyAsync to get similar behavior.
-   // for more info see: https://github.com/mfem/mfem/pull/2780
-   { return HipMemcpyDtoDAsync(dst, src, bytes); }
+   { return HipMemcpyDtoD(dst, src, bytes); }
   void *DtoH(void *dst, const void *src, size_t bytes) override
   { return HipMemcpyDtoH(dst, src, bytes); }
 };
@@ -658,10 +655,7 @@ public:
      return CuMemcpyDtoD(dst, src, bytes);
 #endif
 #ifdef MFEM_USE_HIP
-      // Unlike cudaMemcpy(DtoD), hipMemcpy(DtoD) causes a host-side synchronization so
-      // instead we use hipMemcpyAsync to get similar behavior.
-      // for more info see: https://github.com/mfem/mfem/pull/2780
-      return HipMemcpyDtoDAsync(dst, src, bytes);
+      return HipMemcpyDtoD(dst, src, bytes);
 #endif
      // rm.copy(dst, const_cast<void*>(src), bytes); return dst;
   }
@@ -897,6 +897,7 @@ inline HYPRE_MemoryLocation GetHypreMemoryLocation()
 #elif MFEM_HYPRE_VERSION < 23100
   return HYPRE_MEMORY_DEVICE;
 #else // HYPRE_USING_GPU is defined and MFEM_HYPRE_VERSION >= 23100
+   if (!HYPRE_Initialized()) { return HYPRE_MEMORY_HOST; }
   HYPRE_MemoryLocation loc;
   HYPRE_GetMemoryLocation(&loc);
   return loc;
@@ -1059,7 +1060,8 @@ inline void Memory<T>::MakeAlias(const Memory &base, int offset, int size)
         // register the 'base' if the MemoryManager::Exists():
         MemoryManager::Exists()
 #else // HYPRE_USING_GPU is defined and MFEM_HYPRE_VERSION >= 23100
-         MemoryManager::Exists() && HypreUsingGPU()
+         IsDeviceMemory(MemoryManager::GetDeviceMemoryType()) ||
+         (MemoryManager::Exists() && HypreUsingGPU())
 #endif
      )
      {
@@ -537,7 +537,7 @@ void reduce(int N, T &res, B &&body, const R &reducer, bool use_dev,
      return;
   }

-#if defined(MFEM_USE_HIP) || defined(MFEM_USE_CUDA)
+#if defined(MFEM_USE_CUDA_OR_HIP)
   if (use_dev &&
       mfem::Device::Allows(Backend::CUDA | Backend::HIP | Backend::RAJA_CUDA |
                            Backend::RAJA_HIP))
@@ -4405,4 +4405,32 @@ void BatchLUSolve(const DenseTensor &Mlu, const Array<int> &P, Vector &X)
   BatchedLinAlg::LUSolve(Mlu, P, X);
 }

+#ifdef MFEM_USE_LAPACK
+void BandedSolve(int KL, int KU, DenseMatrix &AB, DenseMatrix &B,
+                 Array<int> &ipiv)
+{
+   int LDAB = (2*KL) + KU + 1;
+   int N = AB.NumCols();
+   int NRHS = B.NumCols();
+   int info;
+   ipiv.SetSize(N);
+   MFEM_LAPACK_PREFIX(gbsv_)(&N, &KL, &KU, &NRHS, AB.GetData(), &LDAB,
+                             ipiv.GetData(), B.GetData(), &N, &info);
+   MFEM_ASSERT(info == 0, "BandedSolve failed in LAPACK");
+}
+
+void BandedFactorizedSolve(int KL, int KU, DenseMatrix &AB, DenseMatrix &B,
+                           bool transpose, Array<int> &ipiv)
+{
+   int LDAB = (2*KL) + KU + 1;
+   int N = AB.NumCols();
+   int NRHS = B.NumCols();
+   char trans = transpose ? 'T' : 'N';
+   int info;
+   MFEM_LAPACK_PREFIX(gbtrs_)(&trans, &N, &KL, &KU, &NRHS, AB.GetData(), &LDAB,
+                              ipiv.GetData(), B.GetData(), &N, &info);
+   MFEM_ASSERT(info == 0, "BandedFactorizedSolve failed in LAPACK");
+}
+#endif
+
 } // namespace mfem
@@ -1329,6 +1329,13 @@ void BatchLUFactor(DenseTensor &Mlu, Array<int> &P, const real_t TOL = 0.0);
    dimension m x n. */
 void BatchLUSolve(const DenseTensor &Mlu, const Array<int> &P, Vector &X);

+#ifdef MFEM_USE_LAPACK
+void BandedSolve(int KL, int KU, DenseMatrix &AB, DenseMatrix &B,
+                 Array<int> &ipiv);
+void BandedFactorizedSolve(int KL, int KU, DenseMatrix &AB, DenseMatrix &B,
+                           bool transpose, Array<int> &ipiv);
+#endif
+
 // Inline methods

 inline real_t &DenseMatrix::operator()(int i, int j)
@@ -1868,6 +1868,8 @@ HYPRE_Int HypreParMatrix::Mult(HypreParVector &x, HypreParVector &y,

 void HypreParMatrix::Mult(real_t a, const Vector &x, real_t b, Vector &y) const
 {
+   MFEM_PERF_FUNCTION;
+
   MFEM_ASSERT(x.Size() == Width(), "invalid x.Size() = " << x.Size()
               << ", expected size = " << Width());
   MFEM_ASSERT(y.Size() == Height(), "invalid y.Size() = " << y.Size()
@@ -1926,6 +1928,8 @@ void HypreParMatrix::Mult(real_t a, const Vector &x, real_t b, Vector &y) const
 void HypreParMatrix::MultTranspose(real_t a, const Vector &x,
                                   real_t b, Vector &y) const
 {
+   MFEM_PERF_FUNCTION;
+
   MFEM_ASSERT(x.Size() == Height(), "invalid x.Size() = " << x.Size()
               << ", expected size = " << Height());
   MFEM_ASSERT(y.Size() == Width(), "invalid y.Size() = " << y.Size()
@@ -4091,6 +4095,8 @@ void HypreSolver::Setup(const HypreParVector &b, HypreParVector &x) const
 {
   if (setup_called) { return; }

+   MFEM_PERF_FUNCTION;
+
   MFEM_VERIFY(A != NULL, "HypreParMatrix A is missing");

   HYPRE_Int err_flag = SetupFcn()(*this, *A, b, x);
@@ -4116,6 +4122,8 @@ void HypreSolver::Setup(const Vector &b, Vector &x) const

 void HypreSolver::Mult(const HypreParVector &b, HypreParVector &x) const
 {
+   MFEM_PERF_FUNCTION;
+
   HYPRE_Int err_flag;
   if (A == NULL)
   {
@@ -42,6 +42,13 @@ extern "C" void
 MFEM_LAPACK_PREFIX(getri_)(int *N, real_t *A, int *LDA, int *IPIV, real_t *WORK,
                           int *LWORK, int *INFO);
 extern "C" void
+MFEM_LAPACK_PREFIX(gbsv_)(int *, int *, int *, int *, real_t *, int *, int *,
+                          real_t *, int *, int *);
+extern "C" void
+MFEM_LAPACK_PREFIX(gbtrs_)(char *, int *, int *, int *, int *, real_t *, int *,
+                           int *, real_t *, int *, int *);
+
+extern "C" void
 MFEM_LAPACK_PREFIX(syevr_)(char *JOBZ, char *RANGE, char *UPLO, int *N,
                           real_t *A, int *LDA, real_t *VL, real_t *VU, int *IL,
                           int *IU, real_t *ABSTOL, int *M, real_t *W,
@@ -588,6 +588,8 @@ void ConstrainedOperator::EliminateRHS(const Vector &x, Vector &b) const
 void ConstrainedOperator::ConstrainedMult(const Vector &x, Vector &y,
                                          const bool transpose) const
 {
+   MFEM_PERF_FUNCTION;
+
   const int csz = constraint_list.Size();
   if (csz == 0)
   {
@@ -787,6 +789,8 @@ void RectangularConstrainedOperator::EliminateRHS(const Vector &x,

 void RectangularConstrainedOperator::Mult(const Vector &x, Vector &y) const
 {
+   MFEM_PERF_FUNCTION;
+
   const int trial_csz = trial_constraints.Size();
   const int test_csz = test_constraints.Size();
   if (trial_csz == 0)
@@ -822,6 +826,8 @@ void RectangularConstrainedOperator::Mult(const Vector &x, Vector &y) const
 void RectangularConstrainedOperator::MultTranspose(const Vector &x,
                                                   Vector &y) const
 {
+   MFEM_PERF_FUNCTION;
+
   const int trial_csz = trial_constraints.Size();
   const int test_csz = test_constraints.Size();
   if (test_csz == 0)
@@ -13,6 +13,7 @@
 #define MFEM_OPERATOR

 #include "vector.hpp"
+#include "../general/annotation.hpp"

 namespace mfem
 {
@@ -825,10 +826,12 @@ public:
   explicit IdentityOperator(int n) : Operator(n) { }

   /// Operator application
-   void Mult(const Vector &x, Vector &y) const override { y = x; }
+   void Mult(const Vector &x, Vector &y) const override
+   { MFEM_PERF_FUNCTION; y = x; }

   /// Application of the transpose
-   void MultTranspose(const Vector &x, Vector &y) const override { y = x; }
+   void MultTranspose(const Vector &x, Vector &y) const override
+   { MFEM_PERF_FUNCTION; y = x; }
 };

 /// Returns true if P is the identity prolongation, i.e. if it is either NULL or
@@ -55,6 +55,8 @@ IterativeSolver::IterativeSolver(MPI_Comm comm_)

 real_t IterativeSolver::Dot(const Vector &x, const Vector &y) const
 {
+   MFEM_PERF_FUNCTION;
+
 #ifndef MFEM_USE_MPI
   return (x * y);
 #else
@@ -373,7 +375,7 @@ OperatorChebyshevSmoother::OperatorChebyshevSmoother(const Operator &oper_,
                                                     real_t power_tolerance,
                                                     int power_seed)
 #endif
-   : Solver(d.Size()),
+   : Solver((MFEM_PERF_BEGIN(_MFEM_FUNC_NAME), d.Size())),
     order(order_),
     N(d.Size()),
     dinv(N),
@@ -400,6 +402,7 @@ OperatorChebyshevSmoother::OperatorChebyshevSmoother(const Operator &oper_,
                                                            power_seed);

   Setup();
+   MFEM_PERF_END(_MFEM_FUNC_NAME);
 }

 OperatorChebyshevSmoother::OperatorChebyshevSmoother(const Operator* oper_,
@@ -521,6 +524,8 @@ void OperatorChebyshevSmoother::SetOrder(int new_order)

 void OperatorChebyshevSmoother::Mult(const Vector& x, Vector &y) const
 {
+   MFEM_PERF_FUNCTION;
+
   if (iterative_mode)
   {
      MFEM_ABORT("Chebyshev smoother not implemented for iterative mode");
@@ -764,6 +764,8 @@ void SparseMatrix::Mult(const Vector &x, Vector &y) const

 void SparseMatrix::AddMult(const Vector &x, Vector &y, const real_t a) const
 {
+   MFEM_PERF_FUNCTION;
+
   MFEM_ASSERT(width == x.Size(), "Input vector size (" << x.Size()
               << ") must match matrix width (" << width << ")");
   MFEM_ASSERT(height == y.Size(), "Output vector size (" << y.Size()
@@ -964,6 +966,8 @@ void SparseMatrix::MultTranspose(const Vector &x, Vector &y) const
 void SparseMatrix::AddMultTranspose(const Vector &x, Vector &y,
                                    const real_t a) const
 {
+   MFEM_PERF_FUNCTION;
+
   MFEM_ASSERT(height == x.Size(), "Input vector size (" << x.Size()
               << ") must match matrix height (" << height << ")");
   MFEM_ASSERT(width == y.Size(), "Output vector size (" << y.Size()
@@ -117,7 +117,7 @@ Vector::Vector(const Vector &v)
   UseDevice(v.UseDevice());
 }

-Vector::Vector(Vector &&v)
+Vector::Vector(Vector &&v) : Vector()
 {
   *this = std::move(v);
 }
@@ -205,14 +205,16 @@ Vector &Vector::operator=(const Vector &v)
   data.CopyFrom(v.data, v.Size());
   UseDevice(v.UseDevice());
 #else
-   SetSize(v.Size());
   const bool vuse = v.UseDevice();
   const bool use_dev = UseDevice() || vuse;
+   if (use_dev) { MFEM_PERF_BEGIN(_MFEM_FUNC_NAME); }
+   SetSize(v.Size());
   v.UseDevice(use_dev);
   // keep 'data' where it is, unless 'use_dev' is true
   if (use_dev) { Write(); }
   data.CopyFrom(v.data, v.Size());
   v.UseDevice(vuse);
+   if (use_dev) { MFEM_PERF_END(_MFEM_FUNC_NAME); }
 #endif
   return *this;
 }
@@ -227,9 +229,11 @@ Vector &Vector::operator=(Vector &&v)
 Vector &Vector::operator=(real_t value)
 {
   const bool use_dev = UseDevice();
+   if (use_dev) { MFEM_PERF_BEGIN(_MFEM_FUNC_NAME); }
   const int N = size;
   auto y = Write(use_dev);
   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] = value; });
+   if (use_dev) { MFEM_PERF_END(_MFEM_FUNC_NAME); }
   return *this;
 }

@@ -290,10 +294,12 @@ Vector &Vector::operator-=(const Vector &v)
   MFEM_ASSERT(size == v.size, "incompatible Vectors!");

   const bool use_dev = UseDevice() || v.UseDevice();
+   if (use_dev) { MFEM_PERF_BEGIN(_MFEM_FUNC_NAME); }
   const int N = size;
   const auto x = v.Read(use_dev);
   auto y = ReadWrite(use_dev);
   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] -= x[i]; });
+   if (use_dev) { MFEM_PERF_END(_MFEM_FUNC_NAME); }
   return *this;
 }

@@ -311,10 +317,12 @@ Vector &Vector::operator+=(const Vector &v)
   MFEM_ASSERT(size == v.size, "incompatible Vectors!");

   const bool use_dev = UseDevice() || v.UseDevice();
+   if (use_dev) { MFEM_PERF_BEGIN(_MFEM_FUNC_NAME); }
   const int N = size;
   const auto x = v.Read(use_dev);
   auto y = ReadWrite(use_dev);
   mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] += x[i]; });
+   if (use_dev) { MFEM_PERF_END(_MFEM_FUNC_NAME); }
   return *this;
 }

@@ -326,9 +334,11 @@ Vector &Vector::Add(const real_t a, const Vector &Va)
   {
      const int N = size;
      const bool use_dev = UseDevice() || Va.UseDevice();
+      if (use_dev) { MFEM_PERF_BEGIN(_MFEM_FUNC_NAME); }
      const auto x = Va.Read(use_dev);
      auto y = ReadWrite(use_dev);
      mfem::forall_switch(use_dev, N, [=] MFEM_HOST_DEVICE (int i) { y[i] += a * x[i]; });
+      if (use_dev) { MFEM_PERF_END(_MFEM_FUNC_NAME); }
   }
   return *this;
 }
@@ -445,6 +455,7 @@ void add(const Vector &v1, real_t alpha, const Vector &v2, Vector &v)
   {
 #if !defined(MFEM_USE_LEGACY_OPENMP)
      const bool use_dev = v1.UseDevice() || v2.UseDevice() || v.UseDevice();
+      if (use_dev) { MFEM_PERF_BEGIN(_MFEM_FUNC_NAME); }
      const int N = v.size;
      // Note: get read access first, in case v is the same as v1/v2.
      const auto d_x = v1.Read(use_dev);
@@ -454,6 +465,7 @@ void add(const Vector &v1, real_t alpha, const Vector &v2, Vector &v)
      {
         d_z[i] = d_x[i] + alpha * d_y[i];
      });
+      if (use_dev) { MFEM_PERF_END(_MFEM_FUNC_NAME); }
 #else
      const real_t *v1p = v1.data, *v2p = v2.data;
      real_t *vp = v.data;
@@ -569,6 +581,7 @@ void subtract(const Vector &x, const Vector &y, Vector &z)

 #if !defined(MFEM_USE_LEGACY_OPENMP)
   const bool use_dev = x.UseDevice() || y.UseDevice() || z.UseDevice();
+   if (use_dev) { MFEM_PERF_BEGIN(_MFEM_FUNC_NAME); }
   const int N = x.size;
   // Note: get read access first, in case z is the same as x/y.
   const auto xd = x.Read(use_dev);
@@ -578,6 +591,7 @@ void subtract(const Vector &x, const Vector &y, Vector &z)
   {
      zd[i] = xd[i] - yd[i];
   });
+   if (use_dev) { MFEM_PERF_END(_MFEM_FUNC_NAME); }
 #else
   const real_t *xp = x.data;
   const real_t *yp = y.data;
@@ -32,6 +32,7 @@
 #include <cstring>
 #include <ctime>
 #include <functional>
+#include <set>
 #include <numeric>
 #include <unordered_map>
 #include <unordered_set>
@@ -921,8 +922,66 @@ const FaceGeometricFactors* Mesh::GetFaceGeometricFactors(
   return gf;
 }

+const Array<int>& Mesh::GetBdrFaceAttributes() const
+{
+   if (bdr_face_attrs_cache.Size() == 0)
+   {
+      std::unordered_map<int, int> f_to_be;
+      for (int i = 0; i < GetNBE(); ++i)
+      {
+         const int f = GetBdrElementFaceIndex(i);
+         f_to_be[f] = i;
+      }
+      const int nf_bdr = GetNFbyType(FaceType::Boundary);
+      // MFEM_VERIFY(size_t(nf_bdr) == f_to_be.size(), "Incompatible sizes");
+      bdr_face_attrs_cache.SetSize(nf_bdr);
+      int f_ind = 0;
+      const int nf = GetNumFaces();
+      for (int f = 0; f < nf; ++f)
+      {
+         if (!GetFaceInformation(f).IsOfFaceType(FaceType::Boundary))
+         {
+            continue;
+         }
+         int attribute = -1; // default value
+         auto iter = f_to_be.find(f);
+         if (iter != f_to_be.end())
+         {
+            const int be = iter->second;
+            attribute = GetBdrAttribute(be);
+         }
+         else
+         {
+            // If a boundary face does not correspond to the a boundary element,
+            // we assign it the default attribute of -1.
+         }
+         bdr_face_attrs_cache[f_ind] = attribute;
+         ++f_ind;
+      }
+   }
+   return bdr_face_attrs_cache;
+}
+
+const Array<int>& Mesh::GetElementAttributes() const
+{
+   if (elem_attrs_cache.Size() == 0)
+   {
+      // re-compute cache
+      elem_attrs_cache.SetSize(GetNE());
+      elem_attrs_cache.HostWrite();
+      for (int i = 0; i < GetNE(); ++i)
+      {
+         elem_attrs_cache[i] = GetAttribute(i);
+         MFEM_ASSERT(elem_attrs_cache[i] > 0,
+                     "Negative attribute on element " << i);
+      }
+   }
+   return elem_attrs_cache;
+}
+
 void Mesh::DeleteGeometricFactors()
 {
+   MFEM_PERF_FUNCTION;
   for (int i = 0; i < geom_factors.Size(); i++)
   {
      delete geom_factors[i];
@@ -1804,6 +1863,8 @@ void Mesh::Destroy()
   TetMemory.Clear();
 #endif

+   elem_attrs_cache.DeleteAll();
+   bdr_face_attrs_cache.DeleteAll();
   attributes.DeleteAll();
   bdr_attributes.DeleteAll();
 }
@@ -1817,36 +1878,48 @@ void Mesh::ResetLazyData()
   DeleteGeometricFactors();
   nbInteriorFaces = -1;
   nbBoundaryFaces = -1;
+   // set size to 0 so re-computations can potentially avoid a new allocation
+   bdr_face_attrs_cache.SetSize(0);
+   elem_attrs_cache.SetSize(0);
 }

-void Mesh::SetAttributes()
+void Mesh::SetAttributes(bool elem_attrs_changed, bool bdr_face_attrs_changed)
 {
-   Array<int> attribs;
+   if (bdr_face_attrs_changed)
+   {
+      bdr_face_attrs_cache.SetSize(0); // Invalidate the cache

-   attribs.SetSize(GetNBE());
-   for (int i = 0; i < attribs.Size(); i++)
-   {
-      attribs[i] = GetBdrAttribute(i);
-   }
-   attribs.Sort();
-   attribs.Unique();
-   attribs.Copy(bdr_attributes);
-   if (bdr_attributes.Size() > 0 && bdr_attributes[0] <= 0)
-   {
-      MFEM_WARNING("Non-positive attributes on the boundary!");
+      // Get sorted list of unique boundary element attributes
+      std::set<int> attribs;
+      for (int i = 0; i < GetNBE(); i++)
+      {
+         attribs.emplace(GetBdrAttribute(i));
+      }
+
+      bdr_attributes.SetSize(attribs.size());
+      bdr_attributes.HostWrite();
+      std::copy(attribs.begin(), attribs.end(), bdr_attributes.begin());
+      if (bdr_attributes.Size() > 0 && bdr_attributes[0] <= 0)
+      {
+         MFEM_WARNING("Non-positive attributes on the boundary!");
+      }
   }

-   attribs.SetSize(GetNE());
-   for (int i = 0; i < attribs.Size(); i++)
+   if (elem_attrs_changed)
   {
-      attribs[i] = GetAttribute(i);
-   }
-   attribs.Sort();
-   attribs.Unique();
-   attribs.Copy(attributes);
-   if (attributes.Size() > 0 && attributes[0] <= 0)
-   {
-      MFEM_WARNING("Non-positive attributes in the domain!");
+      // Re-compute the attributes cache
+      elem_attrs_cache.SetSize(0);
+      GetElementAttributes();
+      // Get sorted list of unique element attributes
+      std::set<int> attribs(elem_attrs_cache.begin(), elem_attrs_cache.end());
+      attributes.SetSize(attribs.size());
+      attributes.HostWrite();
+      std::copy(attribs.begin(), attribs.end(), attributes.begin());
+
+      if (attributes.Size() > 0 && attributes[0] <= 0)
+      {
+         MFEM_WARNING("Non-positive attributes in the domain!");
+      }
   }
 }

@@ -4454,6 +4527,10 @@ Mesh::Mesh(const Mesh &mesh, bool copy_nodes)
      Nodes = mesh.Nodes;
      own_nodes = 0;
   }
+
+   // copy attribute caches
+   elem_attrs_cache = mesh.elem_attrs_cache;
+   bdr_face_attrs_cache = mesh.bdr_face_attrs_cache;
 }

 Mesh::Mesh(Mesh &&mesh) : Mesh()
@@ -6503,6 +6580,7 @@ void XYZ_VectorFunction(const Vector &p, Vector &v)

 void Mesh::GetNodes(GridFunction &nodes) const
 {
+   MFEM_PERF_FUNCTION;
   if (Nodes == NULL || Nodes->FESpace() != nodes.FESpace())
   {
      const int newSpaceDim = nodes.FESpace()->GetVDim();
@@ -6523,6 +6601,7 @@ void Mesh::SetNodalFESpace(FiniteElementSpace *nfes)

 void Mesh::EnsureNodes()
 {
+   MFEM_PERF_FUNCTION;
   if (Nodes)
   {
      const FiniteElementCollection *fec = GetNodalFESpace()->FEColl();
@@ -6575,6 +6654,7 @@ const FiniteElementSpace *Mesh::GetNodalFESpace() const

 void Mesh::SetCurvature(int order, bool discont, int space_dim, int ordering)
 {
+   MFEM_PERF_FUNCTION;
   if (order <= 0)
   {
      delete Nodes;
@@ -7721,6 +7801,12 @@ void Mesh::GetBdrElementAdjacentElement2(
 void Mesh::SetAttribute(int i, int attr)
 {
  elements[i]->SetAttribute(attr);
+  if (elem_attrs_cache.Size() == GetNE())
+  {
+     // update the existing cache instead of deleting it
+     elem_attrs_cache.HostReadWrite();
+     elem_attrs_cache[i] = attr;
+  }
  if (ncmesh) ncmesh->SetAttribute(i, attr);
 }

@@ -10850,6 +10936,10 @@ void Mesh::Swap(Mesh& other, bool non_geometry)
      mfem::Swap(nodes_sequence, other.nodes_sequence);
      mfem::Swap(last_operation, other.last_operation);
   }
+
+   // copy attribute caches
+   mfem::Swap(elem_attrs_cache, other.elem_attrs_cache);
+   mfem::Swap(bdr_face_attrs_cache, other.bdr_face_attrs_cache);
 }

 void Mesh::GetElementData(const Array<Element*> &elem_array, int geom,
@@ -13257,7 +13347,7 @@ void Mesh::ScaleElements(real_t sf)
   delete [] vn;
 }

-void Mesh::Transform(void (*f)(const Vector&, Vector&))
+void Mesh::Transform(std::function<void(const Vector &, Vector&)> f)
 {
   // TODO: support for different new spaceDim.
   if (Nodes == NULL)
@@ -13270,7 +13360,7 @@ void Mesh::Transform(void (*f)(const Vector&, Vector&))
            vold(j) = vertices[i](j);
         }
         vnew.SetData(vertices[i]());
-         (*f)(vold, vnew);
+         f(vold, vnew);
      }
   }
   else
@@ -14649,7 +14739,7 @@ GeometricFactors::GeometricFactors(const GridFunction &nodes,
 void GeometricFactors::Compute(const GridFunction &nodes,
                               MemoryType d_mt)
 {
-
+   MFEM_PERF_FUNCTION;
   const FiniteElementSpace *fespace = nodes.FESpace();
   const FiniteElement *fe = fespace->GetTypicalFE();
   const int dim  = fe->GetDim();
@@ -107,6 +107,11 @@ protected:
   Array<Element *> boundary;
   Array<Element *> faces;

+   /// internal cache for element attributes
+   mutable Array<int> elem_attrs_cache;
+   /// internal cache for boundary element attributes
+   mutable Array<int> bdr_face_attrs_cache;
+
   /** @brief This structure stores the low level information necessary to
       interpret the configuration of elements on a specific face. This
       information can be accessed using methods like GetFaceElements(),
@@ -588,9 +593,10 @@ protected:
   void Loader(std::istream &input, int generate_edges = 0,
               std::string parse_tag = "");

-   /** If NURBS mesh, write NURBS format. If NCMesh, write mfem v1.1 format.
-       If section_delimiter is empty, write mfem v1.0 format. Otherwise, write
-       mfem v1.2 format with the given section_delimiter at the end.
+   /** @brief If NURBS mesh, write NURBS format. If NCMesh, write mfem v1.1
+       format. If section_delimiter is empty, write mfem v1.0 format. Otherwise,
+       write mfem v1.2 format with the given section_delimiter at the end.
+
       If @a comments is non-empty, it will be printed after the first line of
       the file, and each line should begin with '#'. */
   void Printer(std::ostream &os = mfem::out,
@@ -1121,13 +1127,14 @@ public:
       Mesh vertices or nodes are set. */
   virtual void Finalize(bool refine = false, bool fix_orientation = false);

-   /// @brief Determine the sets of unique attribute values in domain and
-   /// boundary elements.
+   /// @brief Determine the sets of unique attribute values in domain if @a
+   /// elem_attrs_changed and boundary elements if @a bdr_face_attrs_changed.
   ///
   /// Separately scan the domain and boundary elements to generate unique,
   /// sorted sets of the element attribute values present in the mesh and
   /// store these in the Mesh::attributes and Mesh::bdr_attributes arrays.
-   virtual void SetAttributes();
+   virtual void SetAttributes(bool elem_attrs_changed = true,
+                              bool bdr_face_attrs_changed = true);

   /// Check (and optionally attempt to fix) the orientation of the elements
   /** @param[in] fix_it  If `true`, attempt to fix the orientations of some
@@ -2253,7 +2260,7 @@ public:
   void ScaleSubdomains (real_t sf);
   void ScaleElements (real_t sf);

-   void Transform(void (*f)(const Vector&, Vector&));
+   void Transform(std::function<void(const Vector &, Vector&)> f);
   void Transform(VectorCoefficient &deformation);

   /** @brief This function should be called after the mesh node coordinates
@@ -2266,6 +2273,35 @@ public:
       method does not modify the nodes. */
   void NodesUpdated() { DeleteGeometricFactors(); }

+   /// @brief Returns the attributes for all elements in this mesh. The i'th
+   /// entry of the array is the attribute of the i'th element of the mesh.
+   ///
+   /// The returned array points to an internal object that may be invalidated
+   /// by mesh operations such as refinement or any element attributes are
+   /// modified. Since not all such modifications can be tracked by the Mesh
+   /// class (e.g. if a user calls GetElement() then changes the element
+   /// attribute directly), one needs to account for such changes by calling the
+   /// method SetAttributes().
+   const Array<int>& GetElementAttributes() const;
+
+   /// @brief Returns the attributes for all boundary elements in this mesh.
+   ///
+   /// The face restriction will give "face E-vectors" on the boundary that
+   /// are numbered in the order of the faces of mesh. This numbering will be
+   /// different than the numbering of the boundary elements. We compute
+   /// mappings so that the array `bdr_attributes[i]` gives the boundary
+   /// attribute of the `i`th boundary face in the mesh face order.
+   /// Attributes <= 0 indicate there is no boundary element and should be
+   /// skipped.
+   ///
+   /// The returned array points to an internal object that may be invalidated
+   /// by mesh operations such as refinement or any element attributes are
+   /// modified. Since not all such modifications can be tracked by the Mesh
+   /// class (e.g. if a user calls GetElement() then changes the element
+   /// attribute directly), one needs to account for such changes by calling the
+   /// method SetAttributes().
+   const Array<int>& GetBdrFaceAttributes() const;
+
   /// @}

   /// @anchor mfem_Mesh_gf_nodes
@@ -2482,10 +2518,12 @@ public:
   /// Print the mesh to the given stream using Netgen/Truegrid format.
   virtual void PrintXG(std::ostream &os = mfem::out) const;

-   /// Print the mesh to the given stream using the default MFEM mesh format.
-   /// \see mfem::ofgzstream() for on-the-fly compression of ascii outputs. If
-   /// @a comments is non-empty, it will be printed after the first line of the
-   /// file, and each line should begin with '#'.
+   /** @brief Print the mesh to the given stream using the default MFEM mesh
+       format.
+
+       \see mfem::ofgzstream() for on-the-fly compression of ascii outputs. If
+       @a comments is non-empty, it will be printed after the first line of the
+       file, and each line should begin with '#'. */
   virtual void Print(std::ostream &os = mfem::out,
                      const std::string &comments = "") const
   { Printer(os, "", comments); }
@@ -53,7 +53,7 @@ KnotVector::KnotVector(int order, int NCP)
 }

 KnotVector::KnotVector(int order, const Vector& intervals,
-                       const Array<int>& continuity )
+                       const Array<int>& continuity)
 {
   // NOTE: This may need to be generalized to support periodicity
   // in the future.
@@ -151,7 +151,7 @@ void KnotVector::UniformRefinement(Vector &newknots, int rf) const
      {
         for (int m = 1; m < rf; ++m)
         {
-            newknots(j) = m * h * (knot(i) + knot(i+1));
+            newknots(j) = ((1.0 - (m * h)) * knot(i)) + (m * h * knot(i+1));
            j++;
         }
      }
@@ -340,7 +340,7 @@ void KnotVector::PrintFunctions(std::ostream &os, int samples) const
   }
 }

-// Routine from "The NURBS book" - 2nd ed - Piegl and Tiller
+// Routine from "The NURBS Book" - 2nd ed - Piegl and Tiller
 // Algorithm A2.2 p. 70
 void KnotVector::CalcShape(Vector &shape, int i, real_t xi) const
 {
@@ -367,7 +367,7 @@ void KnotVector::CalcShape(Vector &shape, int i, real_t xi) const
   }
 }

-// Routine from "The NURBS book" - 2nd ed - Piegl and Tiller
+// Routine from "The NURBS Book" - 2nd ed - Piegl and Tiller
 // Algorithm A2.3 p. 72
 void KnotVector::CalcDShape(Vector &grad, int i, real_t xi) const
 {
@@ -425,7 +425,7 @@ void KnotVector::CalcDShape(Vector &grad, int i, real_t xi) const
   }
 }

-// Routine from "The NURBS book" - 2nd ed - Piegl and Tiller
+// Routine from "The NURBS Book" - 2nd ed - Piegl and Tiller
 // Algorithm A2.3 p. 72
 void KnotVector::CalcDnShape(Vector &gradn, int n, int i, real_t xi) const
 {
@@ -545,11 +545,11 @@ void KnotVector::FindMaxima(Array<int> &ks, Vector &xi, Vector &u) const
         int i = j - d;
         if (isElement(i))
         {
-            arg1 = 1e-16;
+            arg1 = std::numeric_limits<real_t>::epsilon() / 2_r;
            CalcShape(shape, i, arg1);
            max1 = shape[d];

-            arg2 = 1-(1e-16);
+            arg2 = 1_r - arg1;
            CalcShape(shape, i, arg2);
            max2 = shape[d];

@@ -587,9 +587,9 @@ void KnotVector::FindMaxima(Array<int> &ks, Vector &xi, Vector &u) const
   }
 }

-// Routine from "The NURBS book" - 2nd ed - Piegl and Tiller
+// Routine from "The NURBS Book" - 2nd ed - Piegl and Tiller
 // Algorithm A9.1 p. 369
-void KnotVector::FindInterpolant(Array<Vector*> &x)
+void KnotVector::FindInterpolant(Array<Vector*> &x, bool reuse_inverse)
 {
   int order = GetOrder();
   int ncp = GetNCP();
@@ -597,29 +597,93 @@ void KnotVector::FindInterpolant(Array<Vector*> &x)
   // Find interpolation points
   Vector xi_args, u_args;
   Array<int> i_args;
-   FindMaxima(i_args,xi_args, u_args);
+   FindMaxima(i_args, xi_args, u_args);

   // Assemble collocation matrix
-   Vector shape(order+1);
-   DenseMatrix A(ncp,ncp);
-   A = 0.0;
+#ifdef MFEM_USE_LAPACK
+   // If using LAPACK, we use banded matrix storage (order + 1 nonzeros per row).
+   // Find banded structure of matrix.
+   int KL = 0; // Number of subdiagonals
+   int KU = 0; // Number of superdiagonals
   for (int i = 0; i < ncp; i++)
   {
-      CalcShape(shape, i_args[i], xi_args[i]);
      for (int p = 0; p < order+1; p++)
      {
-         A(i,i_args[i] + p) =  shape[p];
+         const int col = i_args[i] + p;
+         if (col < i)
+         {
+            KL = std::max(KL, i - col);
+         }
+         else if (i < col)
+         {
+            KU = std::max(KU, col - i);
+         }
      }
   }

-   // Solve problems
-   A.Invert();
+   const int LDAB = (2*KL) + KU + 1;
+   const int N = ncp;
+
+   fact_AB.SetSize(LDAB, N);
+#else
+   // Without LAPACK, we store and invert a DenseMatrix (inefficient).
+   if (!reuse_inverse)
+   {
+      A_coll_inv.SetSize(ncp, ncp);
+      A_coll_inv = 0.0;
+   }
+#endif
+
+   Vector shape(order+1);
+
+   if (!reuse_inverse) // Set collocation matrix entries
+   {
+      for (int i = 0; i < ncp; i++)
+      {
+         CalcShape(shape, i_args[i], xi_args[i]);
+         for (int p = 0; p < order+1; p++)
+         {
+            const int j = i_args[i] + p;
+#ifdef MFEM_USE_LAPACK
+            fact_AB(KL+KU+i-j,j) = shape[p];
+#else
+            A_coll_inv(i,j) = shape[p];
+#endif
+         }
+      }
+   }
+
+   // Solve the system
+#ifdef MFEM_USE_LAPACK
+   const int NRHS = x.Size();
+   DenseMatrix B(N, NRHS);
+   for (int j=0; j<NRHS; ++j)
+   {
+      for (int i=0; i<N; ++i) { B(i, j) = (*x[j])[i]; }
+   }
+
+   if (reuse_inverse)
+   {
+      BandedFactorizedSolve(KL, KU, fact_AB, B, false, fact_ipiv);
+   }
+   else
+   {
+      BandedSolve(KL, KU, fact_AB, B, fact_ipiv);
+   }
+
+   for (int j=0; j<NRHS; ++j)
+   {
+      for (int i=0; i<N; ++i) { (*x[j])[i] = B(i, j); }
+   }
+#else
+   if (!reuse_inverse) { A_coll_inv.Invert(); }
   Vector tmp;
-   for (int i= 0; i < x.Size(); i++)
+   for (int i = 0; i < x.Size(); i++)
   {
      tmp = *x[i];
-      A.Mult(tmp,*x[i]);
+      A_coll_inv.Mult(tmp, *x[i]);
   }
+#endif
 }

 int KnotVector::findKnotSpan(real_t u) const
@@ -1421,7 +1485,7 @@ void NURBSPatch::DegreeElevate(int t)
   }
 }

-// Routine from "The NURBS book" - 2nd ed - Piegl and Tiller
+// Routine from "The NURBS Book" - 2nd ed - Piegl and Tiller
 void NURBSPatch::DegreeElevate(int dir, int t)
 {
   if (dir >= kv.Size() || dir < 0)
@@ -2385,7 +2449,7 @@ NURBSExtension::NURBSExtension(Mesh *mesh_array[], int num_pieces)
 }

 NURBSExtension::NURBSExtension(const Mesh *patch_topology,
-                               const Array<const NURBSPatch*> patches_)
+                               const Array<const NURBSPatch*> &patches_)
 {
   // Basic topology checks
   MFEM_VERIFY(patches_.Size() > 0, "Must have at least one patch");
@@ -143,8 +143,13 @@ public:
   /** @brief Global curve interpolation through the points @a x (overwritten).
       @a x is an array with the length of the spatial dimension containing
       vectors with spatial coordinates. The control points of the interpolated
-       curve are returned in @a x in the same form. */
-   void FindInterpolant(Array<Vector*> &x);
+       curve are returned in @a x in the same form.
+
+       The inverse of the collocation matrix, used in the interpolation, is
+       stored for repeated calls and used if @a reuse_inverse is true. Reuse is
+       valid only if this KnotVector has not changed since the initial call with
+       @a reuse_inverse false. */
+   void FindInterpolant(Array<Vector*> &x, bool reuse_inverse = false);

   /** Set @a diff, comprised of knots in @a kv not contained in this KnotVector.
       @a kv must be of the same order as this KnotVector. The current
@@ -202,6 +207,14 @@ public:
   /** Flag to indicate whether the KnotVector has been coarsened, which means
       it is ready for non-nested refinement. */
   bool coarse;
+
+#ifdef MFEM_USE_LAPACK
+   // Data for reusing banded matrix factorization in FindInterpolant().
+   DenseMatrix fact_AB; /// Banded matrix factorization
+   Array<int> fact_ipiv; /// Row pivot indices
+#else
+   DenseMatrix A_coll_inv; /// Collocation matrix inverse
+#endif
 };


@@ -285,7 +298,7 @@ public:
       includes the weight. The array of control point coordinates stores each
       point's coordinates contiguously, and points are ordered in a standard
       ijk grid ordering. */
-   NURBSPatch(Array<const KnotVector *> &kv_,  int dim_,
+   NURBSPatch(Array<const KnotVector *> &kv_, int dim_,
              const real_t* control_points);

   /// Constructor for a patch of dimension equal to the size of @a kv.
@@ -700,7 +713,8 @@ public:

   NURBSExtension(Mesh *mesh_array[], int num_pieces);

-   NURBSExtension(const Mesh *patch_topology, const Array<const NURBSPatch*> p);
+   NURBSExtension(const Mesh *patch_topology,
+                  const Array<const NURBSPatch*> &patches_);

   /// Copy assignment not supported.
   NURBSExtension& operator=(const NURBSExtension&) = delete;
@@ -1590,21 +1590,27 @@ void ParMesh::DistributeAttributes(Array<int> &attr)
   delete [] glb_attr_marker;
 }

-void ParMesh::SetAttributes()
+void ParMesh::SetAttributes(bool elem_attrs_changed, bool bdr_attrs_changed)
 {
   // Determine the attributes occurring in local interior and boundary elements
-   Mesh::SetAttributes();
+   Mesh::SetAttributes(elem_attrs_changed, bdr_attrs_changed);

-   DistributeAttributes(bdr_attributes);
-   if (bdr_attributes.Size() > 0 && bdr_attributes[0] <= 0)
+   if (bdr_attrs_changed)
   {
-      MFEM_WARNING("Non-positive boundary element attributes found!");
+      DistributeAttributes(bdr_attributes);
+      if (bdr_attributes.Size() > 0 && bdr_attributes[0] <= 0)
+      {
+         MFEM_WARNING("Non-positive boundary element attributes found!");
+      }
   }

-   DistributeAttributes(attributes);
-   if (attributes.Size() > 0 && attributes[0] <= 0)
+   if (elem_attrs_changed)
   {
-      MFEM_WARNING("Non-positive element attributes found!");
+      DistributeAttributes(attributes);
+      if (attributes.Size() > 0 && attributes[0] <= 0)
+      {
+         MFEM_WARNING("Non-positive element attributes found!");
+      }
   }
 }

@@ -2010,6 +2016,7 @@ void ParMesh::DeleteFaceNbrData()

 void ParMesh::SetCurvature(int order, bool discont, int space_dim, int ordering)
 {
+   MFEM_PERF_FUNCTION;
   DeleteFaceNbrData();
   space_dim = (space_dim == -1) ? spaceDim : space_dim;
   FiniteElementCollection* nfec;
@@ -394,7 +394,8 @@ public:

   void Finalize(bool refine = false, bool fix_orientation = false) override;

-   void SetAttributes() override;
+   void SetAttributes(bool elem_attrs_changed = true,
+                      bool bdr_attrs_changed = true) override;

   /// Checks if any rank in the mesh has boundary elements
   bool HasBoundaryElements() const override;
@@ -257,15 +257,7 @@ template <typename SubMeshT>
 void AddBoundaryElements(SubMeshT &mesh,
                         const std::unordered_map<int,int> &lface_to_boundary_attribute)
 {
-   mesh.Dimension();
-   const int num_codim_1 = [&mesh]()
-   {
-      auto Dim = mesh.Dimension();
-      if (Dim == 1) { return mesh.GetNV(); }
-      else if (Dim == 2) { return mesh.GetNEdges(); }
-      else if (Dim == 3) { return mesh.GetNFaces(); }
-      else { MFEM_ABORT("Invalid dimension."); return -1; }
-   }();
+   const int num_codim_1 = mesh.GetNumFaces();

   if (mesh.Dimension() == 3)
   {
@@ -175,6 +175,56 @@ real_t w_xx(int n, real_t x)

 #endif // CEED_SOLVER_BP_SOLUTION_OPTION

+using BPSFunctionType = real_t(*)(int n, const real_t *xyz);
+
+template <BPSFunctionType F>
+void ProjectBPSFunction(int n, QuadratureFunction &qf)
+{
+   MFEM_PERF_FUNCTION;
+   QuadratureSpaceBase &qs = *qf.GetSpace();
+   Mesh &mesh = *qs.GetMesh();
+   const IntegrationRule &ir = qs.GetIntRule(0);
+
+   auto *geom = mesh.GetGeometricFactors(ir, GeometricFactors::COORDINATES);
+
+   const int dim = qs.GetMesh()->Dimension();
+   const int nq = ir.Size();
+   const int N = qf.Size();
+
+   const real_t *d_x = geom->X.Read();
+   real_t *d_q = qf.Write();
+
+   mfem::forall(N, [=] MFEM_HOST_DEVICE (int ii)
+   {
+      const int i = ii / nq;
+      const int j = ii % nq;
+      real_t xvec[3];
+      for (int d = 0; d < dim; ++d)
+      {
+         xvec[d] = d_x[j + d*nq + i*dim*nq];
+      }
+      d_q[ii] = F(n, xvec);
+   });
+}
+
+MFEM_HOST_DEVICE inline
+real_t sol_1d(const int n, const real_t *xyz)
+{
+   return w(n, xyz[0]);
+}
+
+MFEM_HOST_DEVICE inline
+real_t sol_2d(const int n, const real_t *xyz)
+{
+   return w(n, xyz[0])*w(n, xyz[1]);
+}
+
+MFEM_HOST_DEVICE inline
+real_t sol_3d(const int n, const real_t *xyz)
+{
+   return w(n, xyz[0])*w(n, xyz[1])*w(n, xyz[2]);
+}
+
 struct ExactSolution : Coefficient
 {
   int dim, n;
@@ -198,6 +248,17 @@ struct ExactSolution : Coefficient
         return w(n, xyz[0])*w(n, xyz[1])*w(n, xyz[2]);
      }
   }
+
+   void Project(QuadratureFunction &qf) override
+   {
+      switch (dim)
+      {
+         case 1: ProjectBPSFunction<sol_1d>(n, qf); break;
+         case 2: ProjectBPSFunction<sol_2d>(n, qf); break;
+         case 3: ProjectBPSFunction<sol_3d>(n, qf); break;
+         default: MFEM_ABORT("Unsupported dimension.");
+      }
+   }
 };

 struct ExactGrad : VectorCoefficient
@@ -254,45 +315,14 @@ real_t rhs_3d(const int n, const real_t *xyz)
          - w(n, xyz[0])*w(n, xyz[1])*w_xx(n, xyz[2]);
 }

-using RHSFunctionType = real_t(*)(int dim, const real_t *xyz);
-
-template <RHSFunctionType F>
-void ProjectRHS_(int n, QuadratureFunction &qf)
-{
-   QuadratureSpaceBase &qs = *qf.GetSpace();
-   Mesh &mesh = *qs.GetMesh();
-   const IntegrationRule &ir = qs.GetIntRule(0);
-
-   auto *geom = mesh.GetGeometricFactors(ir, GeometricFactors::COORDINATES);
-
-   const int dim = qs.GetMesh()->Dimension();
-   const int nq = ir.Size();
-   const int N = qf.Size();
-
-   const real_t *d_x = geom->X.Read();
-   real_t *d_q = qf.Write();
-
-   mfem::forall(N, [=] MFEM_HOST_DEVICE (int ii)
-   {
-      const int i = ii / nq;
-      const int j = ii % nq;
-      real_t xvec[3];
-      for (int d = 0; d < dim; ++d)
-      {
-         xvec[d] = d_x[j + d*nq + i*dim*nq];
-      }
-      d_q[ii] = F(n, xvec);
-   });
-}
-
 void ProjectRHS(int n, QuadratureFunction &qf)
 {
   const int dim = qf.GetSpace()->GetMesh()->Dimension();
   switch (dim)
   {
-      case 1: ProjectRHS_<rhs_1d>(n, qf); break;
-      case 2: ProjectRHS_<rhs_2d>(n, qf); break;
-      case 3: ProjectRHS_<rhs_3d>(n, qf); break;
+      case 1: ProjectBPSFunction<rhs_1d>(n, qf); break;
+      case 2: ProjectBPSFunction<rhs_2d>(n, qf); break;
+      case 3: ProjectBPSFunction<rhs_3d>(n, qf); break;
      default: MFEM_ABORT("Unsupported dimension.");
   }
 }
@@ -3,14 +3,15 @@ ssep="----------------------------------------"
 # Enable GPU-aware MPI:
 # gpu_aware_mpi_env_cmd="env MPICH_GPU_SUPPORT_ENABLED=1"
 # gpu_aware_mpi="-g"
-# number of nodes, number of MPI ranks:
-nnodes=1
+# number of MPI ranks, number of ranks per node, number of nodes:
 np=1
-# dev="-d hip ${gpu_aware_mpi}"
+nrnode=4
+((nnodes = (np+nrnode-1)/nrnode))
+# dev="-d gpu ${gpu_aware_mpi}"
 eps="0.3"
 # mpirun_np="mpirun -np"
 mpirun_np="env MFEM_REPORT_KERNELS=1 mpirun -np"
-# mpirun_np="${gpu_aware_mpi_env_cmd} flux run --exclusive -N ${nnodes} -n"
+# mpirun_np="${gpu_aware_mpi_env_cmd} flux run -x -N ${nnodes} -n"
 # dry run:
 # mpirun_np="echo ${mpirun_np}"
 #    p-MG/LOR + FA-hypre, or diagonal (Jacobi smoother)
@@ -37,7 +38,7 @@ while ((np_ > 8)); do
 done
 ((mf = 2**mm))
 ((mff = 3*mf))
-echo " *** mf = ${mf}, mff = ${mff}"
+echo " *** np = ${np}, mf = ${mf}, mff = ${mff}"
 for mg in "${mg_set[@]}"; do
   echo "${bsep}"
   p=(${mg})
@@ -98,26 +99,26 @@ for mg in "${mg_set[@]}"; do
      echo "np = ${np}, p = ${p}, ndofs = ${ndofs}"
      if [[ "$prec_type" == "p-mg" ]]; then
         # p-MG
-         printf "$mpirun_np ${np} ./solver-bp ${dev}"
+         printf "$mpirun_np ${np} ./solver-bp ${dev} -nrn ${nrnode}"
         printf " -ey ${eps} -mg \"${mg}\" -cs 1 ${p_mg_opts}"
         printf " -nx ${nx} -rp ${rp}\n"
         echo "${ssep}"
-         $mpirun_np "${np}" ./solver-bp ${dev} \
+         $mpirun_np "${np}" ./solver-bp ${dev} -nrn ${nrnode} \
            -ey ${eps} -mg "${mg}" -cs 1 ${p_mg_opts} -nx "${nx}" -rp "${rp}"
      elif [[ "$prec_type" == "lor" ]]; then
         # LOR
-         printf "$mpirun_np ${np} ./solver-bp ${dev}"
+         printf "$mpirun_np ${np} ./solver-bp ${dev} -nrn ${nrnode}"
         printf " -ey ${eps} -mg \"${p}\" -cs 2 ${lor_opts}"
         printf " -nx ${nx} -rp ${rp}\n"
         echo "${ssep}"
-         $mpirun_np "${np}" ./solver-bp ${dev} \
+         $mpirun_np "${np}" ./solver-bp ${dev} -nrn ${nrnode} \
            -ey ${eps} -mg "${p}" -cs 2 ${lor_opts} -nx "${nx}" -rp "${rp}"
      elif [[ "$prec_type" == "diag" ]]; then
         # Diag
-         printf "$mpirun_np ${np} ./solver-bp ${dev}"
+         printf "$mpirun_np ${np} ./solver-bp ${dev} -nrn ${nrnode}"
         printf " -ey ${eps} -mg \"${p}\" -cs 0 -nx ${nx} -rp ${rp}\n"
         echo "${ssep}"
-         $mpirun_np "${np}" ./solver-bp ${dev} \
+         $mpirun_np "${np}" ./solver-bp ${dev} -nrn ${nrnode} \
            -ey ${eps} -mg "${p}" -cs 0 -nx "${nx}" -rp "${rp}"
      fi
   done
@@ -71,6 +71,8 @@ struct CGMonitor : IterativeSolverMonitor
   void MonitorResidual(int it, real_t norm, const Vector &r, bool final)
   override
   {
+      MFEM_PERF_FUNCTION;
+
      MFEM_CONTRACT_VAR(norm);
      // Avoid recomputing the norm if it was already computed -- this method
      // is called two times for the final iteration: once with final = false
@@ -139,11 +141,13 @@ int main(int argc, char *argv[])
   DiffusionIntegrator::AddSpecialization<3,4,4>();
   DiffusionIntegrator::AddSpecialization<3,5,5>();
   DiffusionIntegrator::AddSpecialization<3,6,6>();
+   DiffusionIntegrator::AddSpecialization<3,7,7>();

   Mpi::Init(argc, argv);
   Hypre::Init();

   const char *device_config = "cpu";
+   int nrnode = 4; // number of ranks per node, used for CSV output only
   bool gpu_aware_mpi = false;
   int nx = 6, ny = -1, nz = -1;
   int rhs_n = -1;
@@ -159,6 +163,9 @@ int main(int argc, char *argv[])
   OptionsParser args(argc, argv);
   args.AddOption(&device_config, "-d", "--device",
                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&nrnode, "-nrn", "--num-ranks-per-node",
+                  "Number of ranks per compute node. Used to compute the number"
+                  " of nodes written in CSV output.");
   args.AddOption(&gpu_aware_mpi, "-g", "--gpu-aware-mpi", "-no-g",
                  "--no-gpu-aware-mpi", "Enable GPU-aware MPI.");
   args.AddOption(&mg_spec, "-mg", "--multigrid-spec",
@@ -215,9 +222,15 @@ int main(int argc, char *argv[])
   report_env_vars();

   // Generate mesh
+   MFEM_PERF_BEGIN("CreateKershawMesh");
   ParMesh mesh_coarse = CreateKershawMesh(nx, ny, nz, epsy, epsz);
+   MFEM_PERF_END("CreateKershawMesh");
   const int dim = mesh_coarse.Dimension();
-   for (int i=0; i<ref_par; ++i) { mesh_coarse.UniformRefinement(); }
+   for (int i=0; i<ref_par; ++i)
+   {
+      MFEM_PERF_SCOPE("Mesh UniformRefinement");
+      mesh_coarse.UniformRefinement();
+   }

   int coarse_order = 0, order = 0, h_ref = ref_par;
   // Parse order specification
@@ -292,6 +305,7 @@ int main(int argc, char *argv[])
   }
 #endif

+   MFEM_PERF_BEGIN("Setup [hierarchy]");
   vector<unique_ptr<FiniteElementCollection>> fe_collections;
   fe_collections.emplace_back(new H1_FECollection(coarse_order, dim));
   ParFiniteElementSpace fes_coarse(&mesh_coarse, fe_collections.back().get());
@@ -310,6 +324,7 @@ int main(int argc, char *argv[])
         hierarchy.AddOrderRefinedLevel(fe_collections.back().get());
      }
   }
+   MFEM_PERF_END("Setup [hierarchy]");

   const int nlevels = hierarchy.GetNumLevels();
   if (Mpi::Root())
@@ -344,7 +359,9 @@ int main(int argc, char *argv[])

   ParFiniteElementSpace &fes = hierarchy.GetFinestFESpace();
   ParMesh &mesh = *fes.GetParMesh();
+   MFEM_PERF_BEGIN("ParMesh PrintInfo");
   mesh.PrintInfo(cout);
+   MFEM_PERF_END("ParMesh PrintInfo");
   HYPRE_Int ndof = fes.GlobalTrueVSize();
   if (Mpi::Root())
   {
@@ -363,6 +380,7 @@ int main(int argc, char *argv[])
   ConstantCoefficient coeff(1.0); // Diffusion coefficient
   // Set up RHS
   if (Mpi::Root()) { cout << "Assembling right-hand side..." << endl; }
+   MFEM_PERF_BEGIN("Setup [RHS]");
   RHS rhs_coeff(dim, rhs_n);
   ParLinearForm b(&fes);
   const int rhs_ir_inc = 2*q1d_inc+1;
@@ -370,12 +388,20 @@ int main(int argc, char *argv[])
   b.AddDomainIntegrator(new DomainLFIntegrator(rhs_coeff, 2, rhs_ir_inc));
   b.UseFastAssembly(true);
   b.Assemble();
+   MFEM_PERF_END("Setup [RHS]");
   if (Mpi::Root()) { cout << "Assembling right-hand side... Done." << endl; }

+   // Free device memory: the geometric facros computed so far are:
+   // * the coordinates, for the rhs coefficient evaluation, and
+   // * the detJ, for the DomainLFIntegrator.
+   // These are no-longer needed (?), so we can free the memory.
+   mesh.DeleteGeometricFactors();
+
   // make sure the GPU is done with any previous tasks:
   if (Device::Allows(Backend::DEVICE_MASK)) { MFEM_STREAM_SYNC; }
   // make sure all ranks are done with any previous tasks:
   MPI_Barrier(MPI_COMM_WORLD);
+   MFEM_PERF_BEGIN("Setup [DiffusionMultigrid]");
   tic();
   // Set up operators in the multigrid hierarchy
   DiffusionMultigrid MG(hierarchy, coeff, ess_bdr, coarse_solver, q1d_inc,
@@ -385,37 +411,43 @@ int main(int argc, char *argv[])
   if (Device::Allows(Backend::DEVICE_MASK)) { MFEM_STREAM_SYNC; }
   // make sure all ranks are done with all setup tasks:
   MPI_Barrier(MPI_COMM_WORLD);
-   const real_t t_setup = tic_toc.RealTime();
+   const double t_setup = tic_toc.RealTime();
+   MFEM_PERF_END("Setup [DiffusionMultigrid]");

   ParGridFunction x(&fes);
   x = 0.0;

   OperatorPtr A;
   Vector X, B;
+   MFEM_PERF_BEGIN("Setup [MG.FormFineLinearSystem]");
   MG.FormFineLinearSystem(x, b, A, X, B);
+   MFEM_PERF_END("Setup [MG.FormFineLinearSystem]");

   const real_t l2_tol = 1e-8;
   CGMonitor monitor(l2_tol);

   CGSolver cg(MPI_COMM_WORLD);
   cg.SetRelTol(0.0); // use the 'monitor' for convergence
-   cg.SetPrintLevel(1);
+   cg.SetPrintLevel(3);
   cg.SetOperator(*A);
   cg.SetPreconditioner(MG);
   cg.SetMonitor(monitor);
   // Run 2 CG iterations to ensure everything is allocated and initialized for
   // the full CG solve:
   if (Mpi::Root()) { cout << "Running 2 warm-up CG iterations ...\n"; }
+   MFEM_PERF_BEGIN("Warm-up");
   cg.SetMaxIter(2);
   {
      Vector X_save(X);
      cg.Mult(B, X);
      X = X_save;
   }
+   MFEM_PERF_END("Warm-up");
   if (coarse_solver.inner_sli &&
       ((coarse_solver.type == SolverConfig::FA_HYPRE /* && order > 1 */) ||
        coarse_solver.type == SolverConfig::LOR_HYPRE))
   {
+      MFEM_PERF_SCOPE("Auto-tuning");
      // timing data: (t-solve,sli-iter,cheby-order,pcg-iter)
      std::vector<std::tuple<double,int,int,int>> timings;
      Vector X_save(X);
@@ -427,6 +459,8 @@ int main(int argc, char *argv[])
         for (int cheby_order = 1; cheby_order <= smoothers_cheby_order;
              cheby_order++)
         {
+            MFEM_PERF_SCOPE(("Timing [" + to_string(sli_it) + "," +
+                             to_string(cheby_order) + "]").c_str());
            MG.SetSmoothersChebyshevOrder(cheby_order);

            if (Mpi::Root())
@@ -507,6 +541,7 @@ int main(int argc, char *argv[])
   if (Device::Allows(Backend::DEVICE_MASK)) { MFEM_STREAM_SYNC; }
   // make sure all ranks are done with any previous tasks:
   MPI_Barrier(MPI_COMM_WORLD);
+   MFEM_PERF_BEGIN("Final CG Solve");
   tic();
   cg.Mult(B, X);
   // make sure the GPU is done with all solve tasks:
@@ -514,6 +549,7 @@ int main(int argc, char *argv[])
   // make sure all ranks are done with all solve tasks:
   MPI_Barrier(MPI_COMM_WORLD);
   const double t_solve = tic_toc.RealTime();
+   MFEM_PERF_END("Final CG Solve");

   const int niter = cg.GetConverged() ? cg.GetNumIterations() : -1;

@@ -526,10 +562,12 @@ int main(int argc, char *argv[])

   MG.RecoverFineFEMSolution(X, b, x);

+   MFEM_PERF_BEGIN("Compute L2 Error");
   ExactSolution exact_coeff(dim, rhs_n);
   // ExactGrad exact_grad_coeff(dim, rhs_n);
   real_t L2_err = x.ComputeL2Error(exact_coeff);
   // real_t grad_err = x.ComputeGradError(&exact_grad_coeff);
+   MFEM_PERF_END("Compute L2 Error");
   if (Mpi::Root())
   {
      cout << "\nL2 Error: " << setprecision(10) << scientific
@@ -669,7 +707,7 @@ int main(int argc, char *argv[])
      string host_id = regex_replace(hostname ? hostname : "(unknown)",
                                     regex("[0-9]*$"), "");
      cout << ',' << host_id; // 3
-      cout << ',' << (fes.GetNRanks() + 7)/8; // 4 (assuming 8 ranks/node !!)
+      cout << ',' << (fes.GetNRanks() + (nrnode-1))/nrnode; // 4
      cout << ',' << fes.GetNRanks(); // 5
      cout << ',' << nx << ',' << ny << ',' << nz; // 6,7,8
      cout << ',' << order; // 9
@@ -691,7 +729,7 @@ int main(int argc, char *argv[])

 void report_hypre_gpu_status(bool gpu_aware_mpi_requested)
 {
-#ifdef HYPRE_WITH_GPU_AWARE_MPI
+#if defined(HYPRE_WITH_GPU_AWARE_MPI) || defined(HYPRE_USING_GPU_AWARE_MPI)
   bool hypre_gpu_aware_mpi = true;
 #else
   bool hypre_gpu_aware_mpi = false;
@@ -705,11 +743,7 @@ void report_hypre_gpu_status(bool gpu_aware_mpi_requested)
                  "GPU-aware MPI requested but HYPRE's GPU-aware MPI support"
                  " is not enabled");
      cout << "\nHYPRE GPU support: "
-#ifdef HYPRE_USING_GPU
-           << "enabled";
-#else
-           << "disabled";
-#endif
+           << (HypreUsingGPU() ? "enabled" : "disabled");
      cout << "\nHYPRE GPU-aware MPI support: "
           << (hypre_gpu_aware_mpi ? "enabled" : "disabled") << endl;
   }
@@ -80,6 +80,10 @@ add_mfem_miniapp(nurbs_solenoidal
  LIBRARIES mfem)
  add_dependencies(nurbs_solenoidal copy_miniapps_nurbs_data)

+add_mfem_miniapp(nurbs_surface
+  MAIN nurbs_surface.cpp
+  LIBRARIES mfem)
+
 if (MFEM_ENABLE_TESTING)
  add_test(NAME nurbs_ex1_1d_r1_o2_ser
    COMMAND $<TARGET_FILE:nurbs_ex1> -no-vis
@@ -247,6 +251,14 @@ if (MFEM_ENABLE_TESTING)
    COMMAND $<TARGET_FILE:nurbs_solenoidal> -no-vis
        -m ${PROJECT_SOURCE_DIR}/data/cube-nurbs.mesh -r 1 -o 2)

+  add_test(NAME nurbs_surface_10_10_10_10_ex1_o3_ser
+    COMMAND $<TARGET_FILE:nurbs_surface> -no-vis
+        -o 3 -nx 10 -ny 10 -fnx 10 -fny 10 -ex 1 -orig)
+
+  add_test(NAME nurbs_surface_10_10_40_40_ex1_o3_ser
+    COMMAND $<TARGET_FILE:nurbs_surface> -no-vis
+        -o 3 -nx 10 -ny 10 -fnx 40 -fny 14 -ex 1)
+
 endif()

 if (MFEM_USE_MPI)
@@ -21,7 +21,7 @@ MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)

 SEQ_MINIAPPS = nurbs_ex1 nurbs_patch_ex1 nurbs_ex3 nurbs_ex5 nurbs_ex24 \
-   nurbs_curveint nurbs_printfunc nurbs_solenoidal nurbs_naca_cmesh
+   nurbs_curveint nurbs_printfunc nurbs_solenoidal nurbs_naca_cmesh nurbs_surface
 PAR_MINIAPPS = nurbs_ex1p nurbs_ex11p
 ifeq ($(MFEM_USE_MPI),NO)
   MINIAPPS = $(SEQ_MINIAPPS)
@@ -158,6 +158,13 @@ nurbs_naca_cmesh-test-seq: nurbs_naca_cmesh
 nurbs_printfunc-test-seq: nurbs_printfunc
 	@$(call mfem-test,$<,, NURBS miniapp)

+SURF_ARGS_1 := -o 3 -nx 10 -ny 10 -fnx 10 -fny 10 -ex 1 -orig
+SURF_ARGS_2 := -o 3 -nx 10 -ny 10 -fnx 40 -fny 40 -ex 1
+
+nurbs_surface-test-seq: nurbs_surface
+	@$(call mfem-test,$<,, NURBS miniapp,$(SURF_ARGS_1))
+	@$(call mfem-test,$<,, NURBS miniapp,$(SURF_ARGS_2))
+
 EX1P_ARGS_1 :=
 EX1P_ARGS_2 := -m ../../data/pipe-nurbs-2d.mesh -o 2 -no-ibp
 EX1P_ARGS_3 := -m ../../data/ball-nurbs.mesh -o 2 --weak-bc -r 0
@@ -192,6 +199,6 @@ clean-build:

 clean-exec:
 	@rm -f refined.mesh sin-fit.mesh ex5.mesh exsol.mesh mesh.* sol.* mode_*
-	@rm -f naca-cmesh.mesh sol_?.gf
+	@rm -f naca-cmesh.mesh sol_?.gf *-Surface.mesh
 	@rm -rf Example1* Example3* Example5* Solenoidal_* ParaView
 	@rm -rf CurveInt Naca_cmesh glvis_naca-cmesh.mesh solution.dat
@@ -0,0 +1,655 @@
+// 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.
+//
+//          --------------------------------------------------------
+//          NURBS Surface: Interpolate a 3D Surface in a NURBS Patch
+//          --------------------------------------------------------
+//
+// Compile with: make nurbs_surface
+//
+// Sample runs:  nurbs_surface -o 3 -nx 10 -ny 10 -fnx 10 -fny 10 -ex 1 -orig
+//               nurbs_surface -o 3 -nx 10 -ny 10 -fnx 40 -fny 40 -ex 1
+//               nurbs_surface -o 3 -nx 20 -ny 20 -fnx 10 -fny 10 -ex 1
+//               nurbs_surface -o 3 -nx 20 -ny 20 -fnx 40 -fny 40 -ex 1 -j 0.5
+//               nurbs_surface -o 3 -nx 10 -ny 10 -fnx 10 -fny 10 -ex 2 -orig
+//               nurbs_surface -o 3 -nx 10 -ny 10 -fnx 40 -fny 40 -ex 2
+//               nurbs_surface -o 3 -nx 20 -ny 20 -fnx 10 -fny 10 -ex 2
+//               nurbs_surface -o 3 -nx 10 -ny 10 -fnx 10 -fny 10 -ex 3 -orig
+//               nurbs_surface -o 3 -nx 10 -ny 10 -fnx 40 -fny 40 -ex 3
+//               nurbs_surface -o 3 -nx 20 -ny 20 -fnx 10 -fny 10 -ex 3
+//               nurbs_surface -o 3 -nx 20 -ny 10 -fnx 20 -fny 10 -ex 4 -orig
+//             * nurbs_surface -o 3 -nx 20 -ny 10 -fnx 80 -fny 40 -ex 4
+//             * nurbs_surface -o 3 -nx 40 -ny 20 -fnx 20 -fny 10 -ex 4
+//             * nurbs_surface -o 3 -nx 100 -ny 100 -fnx 100 -fny 100 -ex 5 -orig
+//             * nurbs_surface -o 3 -nx 100 -ny 100 -fnx 400 -fny 400 -ex 5
+//             * nurbs_surface -o 3 -nx 200 -ny 200 -fnx 100 -fny 100 -ex 5
+//
+// Description:  This example demonstrates the use of MFEM to interpolate an
+//               input surface point grid in 3D using a NURBS surface. The NURBS
+//               surface can then be sampled to generate an output mesh of
+//               arbitrary resolution while staying close to the input geometry.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+// Example data for 3D point grid on surface, given by an analytic function.
+void SurfaceGridExample(int example, int nx, int ny, Array3D<real_t> &vertices,
+                        real_t jitter);
+
+// Write a linear surface mesh with given vertex positions in v.
+void WriteLinearMesh(int nx, int ny, const Array3D<real_t> &v,
+                     const std::string &basename, bool visualization = false,
+                     int x = 0, int y = 0, int w = 500, int h = 500);
+
+// Given an input grid of 3D points on a surface, this class computes a NURBS
+// surface of given order that interpolates the vertices of the input grid.
+class SurfaceInterpolator
+{
+public:
+   /// Constructor for a given 2D point grid size and NURBS order.
+   SurfaceInterpolator(int num_elem_x, int num_elem_y, int order);
+
+   /// Create a surface interpolating the 2D grid of 3D points in @a input3D.
+   void CreateSurface(const Array3D<real_t> &input3D);
+
+   /// Sample the surface with the given grid size, storing points in
+   /// @a output3D.
+   void SampleSurface(int num_elem_x, int num_elem_y, bool compareOriginal,
+                      Array3D<real_t> &output3D);
+
+   /** @brief Write the NURBS surface mesh to file, defined coordinate-wise by
+       the entries of @a cmesh. */
+   void WriteNURBSMesh(const std::string &basename, bool visualization = false,
+                       int x = 0, int y = 0, int w = 500, int h = 500);
+
+protected:
+   /** @brief Compute the NURBS mesh interpolating the given coordinate of the
+       grid of 3D points in @a input3D. */
+   void ComputeNURBS(int coordinate, const Array3D<real_t> &input3D);
+
+private:
+   int nx, ny; // Number of elements in two directions of the surface grid
+   int orderNURBS; // NURBS degree
+   real_t hx, hy, hz; // Grid size in reference space
+
+   Array3D<real_t> initial3D; // Initial grid of points
+
+   static constexpr int dim = 3;
+   Array<int> ncp; // Number of control points in each direction
+   Array<int> nks; // Number of knot-spans in each direction
+
+   std::vector<Vector> ugrid; // Parameter space [0,1]^2 grid point coordinates
+
+   std::vector<KnotVector> kv; // KnotVectors in each direction
+
+   std::unique_ptr<NURBSPatch> patch; // Pointer to the only patch in the mesh
+
+   Mesh mesh; // NURBS mesh representing the surface
+   std::vector<Mesh> cmesh; // NURBS meshes representing point components
+};
+
+
+int main(int argc, char *argv[])
+{
+   // Parse command-line options
+   int nx = 4;
+   int ny = 4;
+   int fnx = 40;
+   int fny = 40;
+   int order = 3;
+   int example = 1;
+   bool visualization = true;
+   bool compareOriginal = false;
+   real_t jitter = 0.0;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&example, "-ex", "--example",
+                  "Example data");
+   args.AddOption(&nx, "-nx", "--nx",
+                  "Number of elements in x");
+   args.AddOption(&ny, "-ny", "--ny",
+                  "Number of elements in y");
+   args.AddOption(&fnx, "-fnx", "--fnx",
+                  "Number of resampled elements in x");
+   args.AddOption(&fny, "-fny", "--fny",
+                  "Number of resampled elements in y");
+   args.AddOption(&order, "-o", "--order",
+                  "NURBS finite element order (polynomial degree)");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&compareOriginal, "-orig", "--compare-original", "-no-orig",
+                  "--no-compare-original",
+                  "Compare to the original mesh?");
+   args.AddOption(&jitter, "-j", "--jitter",
+                  "Relative jittering in (0,1) to add to the input point "
+                  "coordinates on a uniform nx x ny grid (0 by default)");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   args.PrintOptions(cout);
+
+   if (compareOriginal && (fnx != nx || fny != ny))
+   {
+      cout << "Comparing to the original mesh requires the same number of "
+           << "samples!\n";
+      return 1;
+   }
+
+   // Dimensions of the 3 surfaces (Input, NURBS, Output)
+   cout << "Input Surface:  " << nx << " x " << ny << " linear elements\n";
+   cout << "NURBS Surface:  " << nx + 1 - order << " x " << ny + 1 - order
+        << " knot elements of order " << order << "\n";
+   cout << "Output Surface: " << fnx << " x " << fny << " linear elements\n";
+
+   // Set the vertex coordinates of the initial linear mesh
+   constexpr int dim = 3;
+   Array3D<real_t> input3D(nx + 1, ny + 1, dim);
+   SurfaceGridExample(example, nx, ny, input3D, jitter);
+
+   // Create a NURBS surface for the given nx, ny and order parameters that
+   // interpolates the input vertex coordinates
+   SurfaceInterpolator surf(nx, ny, order);
+   surf.CreateSurface(input3D);
+
+   // Compute the vertex coordinates of the output linear mesh by sampling the
+   // values from the NURBS surface
+   Array3D<real_t> output3D(fnx + 1, fny + 1, dim);
+   surf.SampleSurface(fnx, fny, compareOriginal, output3D);
+
+   // Save and optionally visualize the 3 surfaces (Input, NURBS, Output)
+   WriteLinearMesh(nx, ny, input3D, "Input-Surface", visualization, 0, 0);
+   surf.WriteNURBSMesh("NURBS-Surface", visualization, 502, 0);
+   WriteLinearMesh(fnx, fny, output3D, "Output-Surface", visualization, 1004, 0);
+
+   return 0;
+}
+
+// f(x,y) = sin(2 * pi * x) * sin(2 * pi * y)
+void Function1(real_t u, real_t v, real_t &x, real_t &y, real_t &z)
+{
+   x = u;
+   y = v;
+   z = sin(2.0 * M_PI * u) * sin(2.0 * M_PI * v);
+}
+
+// Part of the parametric surface of a sphere, using spherical coordinates.
+void Function2(real_t u, real_t v, real_t &x, real_t &y, real_t &z)
+{
+   constexpr real_t r = 1.0;
+   constexpr real_t pi_4 = M_PI * 0.25;
+   constexpr real_t phi0 = -3*pi_4;
+   constexpr real_t phi1 = 3*pi_4;
+   constexpr real_t theta0 = pi_4;
+   constexpr real_t theta1 = 3 * pi_4;
+
+   const real_t phi = (phi0 * (1.0 - v)) + (phi1 * v);
+   const real_t theta = (theta0 * (1.0 - u)) + (theta1 * u);
+   x = r * sin(theta) * cos(phi);
+   y = r * sin(theta) * sin(phi);
+   z = r * cos(theta);
+}
+
+// Helicoid surface
+void Function3(real_t u, real_t v, real_t &x, real_t &y, real_t &z)
+{
+   x = u * cos(2.0 * M_PI * v);
+   y = u * sin(2.0 * M_PI * v);
+   z = v;
+}
+
+// Mobius strip
+void Function4(real_t u, real_t v, real_t &x, real_t &y, real_t &z)
+{
+   constexpr int twists = 1;
+   const real_t a = 1.0 + 0.5 * ((2.0 * v) - 1.0) * cos(2.0 * M_PI * twists * u);
+   x = a * cos(2.0 * M_PI * u);
+   y = a * sin(2.0 * M_PI * u);
+   z = 0.5 * (2.0 * v - 1.0) * sin(2.0 * M_PI * twists * u);
+}
+
+// Breather surface
+void Function5(real_t u, real_t v, real_t &x, real_t &y, real_t &z)
+{
+   const real_t m = 13.2 * ((2.0 * u) - 1.0);
+   const real_t n = 37.4 * ((2.0 * v) - 1.0);
+   constexpr real_t b = 0.4;
+   constexpr real_t r = 1.0 - (b*b);
+   const real_t w = sqrt(r);
+   const real_t denom = b * (pow(w*cosh(b*m),2) + pow(b*sin(w*n),2));
+   x = -m + (2*r*cosh(b*m)*sinh(b*m)) / denom;
+   y = (2*w*cosh(b*m)*(-(w*cos(n)*cos(w*n)) - sin(n)*sin(w*n))) / denom;
+   z = (2*w*cosh(b*m)*(-(w*sin(n)*cos(w*n)) + cos(n)*sin(w*n))) / denom;
+}
+
+void SurfaceFunction(int example, real_t u, real_t v,
+                     real_t &x, real_t &y, real_t &z)
+{
+   switch (example)
+   {
+      case 1:
+         Function1(u, v, x, y, z);
+         break;
+      case 2:
+         Function2(u, v, x, y, z);
+         break;
+      case 3:
+         Function3(u, v, x, y, z);
+         break;
+      case 4:
+         Function4(u, v, x, y, z);
+         break;
+      default:
+         Function5(u, v, x, y, z);
+   };
+}
+
+// Example data for 3D point grid on surface, given by an analytic function.
+void SurfaceExample(int example, const std::vector<Vector> &grid,
+                    Array3D<real_t> &v3D, real_t jitter)
+{
+   int seed = (int)time(0);
+   srand((unsigned)seed);
+
+   real_t h0 = grid[0][1]-grid[0][0], h1 = grid[1][1]-grid[1][0];
+   for (int i = 0; i < grid[0].Size(); i++)
+   {
+      for (int j = 0; j < grid[1].Size(); j++)
+      {
+         if (i != 0 && i != grid[0].Size()-1 && j != 0 && j != grid[1].Size()-1)
+         {
+            SurfaceFunction(example, grid[0][i] + rand_real()*h0*jitter,
+                            grid[1][j] + rand_real()*h1*jitter,
+                            v3D(i, j, 0), v3D(i, j, 1), v3D(i, j, 2));
+         }
+         else
+         {
+            SurfaceFunction(example, grid[0][i], grid[1][j],
+                            v3D(i, j, 0), v3D(i, j, 1), v3D(i, j, 2));
+         }
+      }
+   }
+}
+
+void SurfaceGridExample(int example, int nx, int ny, Array3D<real_t> &vertices,
+                        real_t jitter = 0)
+{
+   // Define a uniform grid of the reference parameter space [0,1]^2
+   std::vector<Vector> uniformGrid(2);
+   for (int i = 0; i < 2; ++i)
+   {
+      const int n = (i == 0) ? nx : ny;
+      const real_t h = 1.0 / n;
+      uniformGrid[i].SetSize(n + 1);
+      for (int j = 0; j <= n; ++j) { uniformGrid[i][j] = j * h; }
+   }
+
+   SurfaceExample(example, uniformGrid, vertices, jitter);
+}
+
+// Write a linear surface mesh with given vertex positions in v.
+void WriteLinearMesh(int nx, int ny, const Array3D<real_t> &v,
+                     const std::string &basename, bool visualization,
+                     int x, int y, int w, int h)
+{
+   const int nv = (nx + 1) * (ny + 1);
+   const int nelem = nx * ny;
+   constexpr int dim = 3; // Spatial dimension
+
+   Mesh lmesh(2, nv, nelem, 0, dim);
+   Vector vertex(dim);
+
+   for (int i = 0; i <= nx; ++i)
+   {
+      for (int j = 0; j <= ny; ++j)
+      {
+         for (int k = 0; k < dim; ++k) { vertex[k] = v(i, j, k); }
+         lmesh.AddVertex(vertex);
+      }
+   }
+
+   Array<int> verts(4);
+
+   auto vID = [&](int i, int j)
+   {
+      return j + (i * (ny + 1));
+   };
+
+   for (int i = 0; i < nx; ++i)
+   {
+      for (int j = 0; j < ny; ++j)
+      {
+         verts[0] = vID(i, j);
+         verts[1] = vID(i+1, j);
+         verts[2] = vID(i+1, j+1);
+         verts[3] = vID(i, j+1);
+
+         Element* el = lmesh.NewElement(Element::QUADRILATERAL);
+         el->SetVertices(verts);
+         lmesh.AddElement(el);
+      }
+   }
+
+   lmesh.FinalizeTopology();
+
+   ofstream mesh_ofs(basename + ".mesh");
+   mesh_ofs.precision(8);
+   lmesh.Print(mesh_ofs);
+
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      constexpr int visport = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock.precision(8);
+      sol_sock << "mesh\n" << lmesh
+               << "window_title '" << basename << "'"
+               << "window_geometry "
+               << x << " " << y << " " << w << " " << h << "\n"
+               << "keys PPPPPPPPAattttt******\n"
+               << flush;
+   }
+}
+
+
+// Compute error of interpolation with respect to an input grid of point data.
+void CheckError(const Array3D<real_t> &a, const Array3D<real_t> &b, int c,
+                int nx, int ny)
+{
+   real_t maxErr = 0.0;
+   for (int i = 0; i <= nx; ++i)
+   {
+      for (int j = 0; j <= ny; ++j)
+      {
+         const real_t err_ij = std::abs(a(i, j, c) - b(i, j, 2));
+         maxErr = std::max(maxErr, err_ij);
+      }
+   }
+
+   cout << "Max error: " << maxErr << " for coordinate " << c << endl;
+}
+
+
+// Sample a NURBS mesh to generate a first-order mesh.
+void SampleNURBS(bool uniform, int nx, int ny, const Mesh &mesh,
+                 const Array<int> &nks, const std::vector<Vector> &ugrid,
+                 Array3D<real_t> &vpos)
+{
+   const GridFunction *nodes = mesh.GetNodes();
+
+   const real_t hx = 1.0 / (real_t) nx;
+   const real_t hy = 1.0 / (real_t) ny;
+
+   const real_t hxks = 1.0 / (real_t) nks[0];
+   const real_t hyks = 1.0 / (real_t) nks[1];
+
+   Vector vertex;
+   IntegrationPoint ip;
+
+   ip.z = 1.0;
+   for (int i = 0; i <= nx; ++i)
+   {
+      const real_t xref = uniform ? i * hx : ugrid[0][i];
+      const int nurbsElem0 = std::min((int) (xref / hxks), nks[0] - 1);
+      const real_t ipx = (xref - (nurbsElem0 * hxks)) / hxks;
+      ip.x = ipx;
+
+      for (int j = 0; j <= ny; ++j)
+      {
+         const real_t yref = uniform ? j * hy : ugrid[1][j];
+         const int nurbsElem1 = std::min((int) (yref / hyks), nks[1] - 1);
+         const real_t ipy = (yref - (nurbsElem1 * hyks)) / hyks;
+         ip.y = ipy;
+
+         const int nurbsElem = nurbsElem0 + (nurbsElem1 * nks[0]);
+         nodes->GetVectorValue(nurbsElem, ip, vertex);
+
+         for (int k = 0; k < 3; ++k)
+         {
+            vpos(i, j, k) = vertex[k];
+         }
+      }
+   }
+}
+
+
+SurfaceInterpolator::SurfaceInterpolator(int num_elem_x, int num_elem_y,
+                                         int order) :
+   nx(num_elem_x), ny(num_elem_y), orderNURBS(order),
+   ncp(dim), nks(dim), ugrid(dim - 1)
+{
+   ncp[0] = nx + 1;
+   ncp[1] = ny + 1;
+   ncp[2] = order + 1;
+
+   for (int i = 0; i < dim; ++i)
+   {
+      nks[i] = ncp[i] - order;
+
+      Vector intervals(nks[i]);
+      Array<int> continuity(nks[i] + 1);
+
+      intervals = 1.0 / (real_t) nks[i];
+      continuity = order - 1;
+      continuity[0] = -1;
+      continuity[nks[i]] = -1;
+
+      kv.emplace_back(order, intervals, continuity);
+   }
+
+   patch.reset(new NURBSPatch(&kv[0], &kv[1], &kv[2], dim + 1));
+
+   hx = 1.0 / (real_t) (ncp[0] - 1);
+   hy = 1.0 / (real_t) (ncp[1] - 1);
+   hz = 1.0 / (real_t) (ncp[2] - 1);
+
+   Vector xi_args;
+   Array<int> i_args;
+   for (int i = 0; i < 2; ++i)
+   {
+      kv[i].FindMaxima(i_args, xi_args, ugrid[i]);
+   }
+}
+
+void SurfaceInterpolator::CreateSurface(const Array3D<real_t> &input3D)
+{
+   cmesh.clear();
+   for (int c = 0; c < dim; ++c) // Loop over coordinates
+   {
+      ComputeNURBS(c, input3D);
+      cmesh.emplace_back(mesh);
+   }
+
+   initial3D = input3D;
+}
+
+void SurfaceInterpolator::SampleSurface(int num_elem_x, int num_elem_y,
+                                        bool compareOriginal,
+                                        Array3D<real_t> &output3D)
+{
+   Array3D<real_t> vpos(num_elem_x + 1, num_elem_y + 1, dim);
+   for (int c = 0; c < dim; ++c) // Loop over coordinates
+   {
+      SampleNURBS(true, num_elem_x, num_elem_y, cmesh[c], nks, ugrid, vpos);
+
+      if (compareOriginal)
+      {
+         SampleNURBS(false, num_elem_x, num_elem_y, cmesh[c], nks, ugrid, vpos);
+         CheckError(initial3D, vpos, c, nx, ny);
+      }
+
+      for (int i = 0; i <= num_elem_x; ++i)
+      {
+         for (int j = 0; j <= num_elem_y; ++j)
+         {
+            output3D(i,j,c) = vpos(i,j,2);
+         }
+      }
+   }
+}
+
+void SurfaceInterpolator::ComputeNURBS(int coordinate,
+                                       const Array3D<real_t> &input3D)
+{
+   Array<Vector*> x;
+   for (int i = 0; i < dim; ++i) { x.Append(new Vector(ncp[0])); }
+
+   for (int k = 0; k < ncp[2]; ++k)
+   {
+      const real_t z = k * hz;
+
+      // For each horizontal slice (fixed k), interpolate a 2D surface by
+      // sweeping curve interpolations in each direction. See Algorithm A9.4 of
+      // "The NURBS Book" - 2nd ed - Piegl and Tiller.
+
+      // Resize for sweep in first direction
+      for (int i = 0; i < dim; ++i) { x[i]->SetSize(ncp[0]); }
+
+      // Sweep in the first direction
+      for (int j = 0; j < ncp[1]; ++j)
+      {
+         for (int i = 0; i < ncp[0]; i++)
+         {
+            (*x[0])[i] = ugrid[0][i];
+            (*x[1])[i] = ugrid[1][j];
+
+            const real_t s_ij = input3D(i, j, coordinate);
+            (*x[2])[i] = -1.0 + z + s_ij;
+         }
+
+         const bool reuse_factorization = j > 0;
+         kv[0].FindInterpolant(x, reuse_factorization);
+
+         for (int i = 0; i < ncp[0]; i++)
+         {
+            (*patch)(i,j,k,0) = (*x[0])[i];
+            (*patch)(i,j,k,1) = (*x[1])[i];
+            (*patch)(i,j,k,2) = (*x[2])[i];
+            (*patch)(i,j,k,3) = 1.0; // weight
+         }
+      }
+
+      // Resize for sweep in second direction
+      for (int i = 0; i < dim; ++i) { x[i]->SetSize(ncp[1]); }
+
+      // Do another sweep in the second direction
+      for (int i = 0; i < ncp[0]; i++)
+      {
+         for (int j = 0; j < ncp[1]; ++j)
+         {
+            (*x[0])[j] = (*patch)(i,j,k,0);
+            (*x[1])[j] = (*patch)(i,j,k,1);
+            (*x[2])[j] = (*patch)(i,j,k,2);
+         }
+
+         const bool reuse_factorization = i > 0;
+         kv[1].FindInterpolant(x, reuse_factorization);
+
+         for (int j = 0; j < ncp[1]; ++j)
+         {
+            (*patch)(i,j,k,0) = (*x[0])[j];
+            (*patch)(i,j,k,1) = (*x[1])[j];
+            (*patch)(i,j,k,2) = (*x[2])[j];
+         }
+      }
+   }
+
+   for (auto p : x) { delete p; }
+
+   Array<const NURBSPatch*> patches(1);
+   patches[0] = patch.get();
+   Mesh patch_topology = Mesh::MakeCartesian3D(1, 1, 1, Element::HEXAHEDRON);
+   NURBSExtension nurbsExt(&patch_topology, patches);
+
+   mesh = Mesh(nurbsExt);
+}
+
+void SurfaceInterpolator::WriteNURBSMesh(const std::string &basename,
+                                         bool visualization,
+                                         int x, int y, int w, int h)
+{
+   GridFunction *nodes = cmesh[0].GetNodes();
+   NURBSPatch patch2D(&kv[0], &kv[1], dim);
+   Array<const NURBSPatch*> patches(1);
+   patches[0] = &patch2D;
+   Mesh patch_topology = Mesh::MakeCartesian2D(1, 1, Element::QUADRILATERAL);
+   Array<int> dofs;
+   cmesh[0].NURBSext->GetPatchDofs(0, dofs);
+
+   MFEM_VERIFY(dofs.Size() == (nx + 1) * (ny + 1) * (orderNURBS + 1), "");
+
+   for (int j = 0; j < ncp[1]; ++j)
+   {
+      for (int i = 0; i < ncp[0]; i++)
+      {
+         const int dof = dofs[i + (ncp[0] * (j + (ncp[1] * orderNURBS)))];
+         for (int k = 0; k < 2; ++k) { patch2D(i,j,k) = (*nodes)[dim*dof + k]; }
+         patch2D(i,j,2) = 1.0; // weight
+      }
+   }
+
+   NURBSExtension nurbsExt(&patch_topology, patches);
+   Mesh mesh2D(nurbsExt);
+
+   FiniteElementCollection *fec = nodes->OwnFEC();
+   FiniteElementSpace fespace(&mesh2D, fec, dim, Ordering::byVDIM);
+   GridFunction nodes2D(&fespace);
+
+   const int n = mesh2D.GetNodes()->Size() / (dim - 1);
+   MFEM_VERIFY((dim - 1) * n == mesh2D.GetNodes()->Size(), "");
+   MFEM_VERIFY(dim * n == nodes2D.Size(), "");
+
+   Array<int> dofs2D;
+   mesh2D.NURBSext->GetPatchDofs(0, dofs2D);
+
+   for (int k = 0; k < dim; ++k)
+   {
+      const GridFunction &nodes_k = *cmesh[k].GetNodes();
+
+      for (int j = 0; j < ncp[1]; ++j)
+      {
+         for (int i = 0; i < ncp[0]; i++)
+         {
+            const int dof = dofs[i + (ncp[0] * (j + (ncp[1] * orderNURBS)))];
+            const int dof2D = dofs2D[i + (ncp[0] * j)];
+            nodes2D[(dim*dof2D) + k] = nodes_k[dim*dof + 2];
+         }
+      }
+   }
+
+   // Make mesh2D into a surface mesh with nodes given by nodes2D
+   mesh2D.NewNodes(nodes2D);
+
+   ofstream mesh_ofs(basename + ".mesh");
+   mesh_ofs.precision(8);
+   mesh2D.Print(mesh_ofs);
+
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      constexpr int visport = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock.precision(8);
+      sol_sock << "mesh\n" << mesh2D
+               << "window_title '" << basename << "'"
+               << "window_geometry "
+               << x << " " << y << " " << w << " " << h << "\n"
+               << "keys PPPPPPPPAattttt******\n"
+               << flush;
+   }
+}
@@ -105,8 +105,11 @@ MFEM_PERF_CXXFLAGS_xlc = -mcpu=native
 # - Clang extra options:
 ifeq ($(MFEM_MACHINE),riscv64)
  MFEM_PERF_CXXFLAGS_clang += -march=rv64gc
-else ifneq ($(MFEM_MACHINE),arm64)
-  # -march=native is unavailable on clang/ARM64 as of 05/2021: support could be added later.
+else ifneq (,$(findstring ppc,$(MFEM_MACHINE)))
+  MFEM_PERF_CXXFLAGS_clang += -mcpu=native -mtune=native
+else ifeq ($(MFEM_MACHINE),arm64)
+  MFEM_PERF_CXXFLAGS_clang += -mcpu=native -mtune=native
+else
  MFEM_PERF_CXXFLAGS_clang += -march=native
 endif
 MFEM_PERF_CXXFLAGS_clang += $(PEDANTIC_FLAG) -Wall
@@ -384,8 +384,10 @@ int main(int argc, char *argv[])
   dacol.Save();

   ConstantCoefficient zero(0.0);
+   Vector zero_vec(dim); zero_vec = 0_r;
+   VectorConstantCoefficient vzero(zero_vec);
   const real_t s_norm  = distance_s.ComputeL2Error(zero),
-                v_norm  = distance_v.ComputeL2Error(zero);
+                v_norm  = distance_v.ComputeL2Error(vzero);
   if (myid == 0)
   {
      cout << fixed << setprecision(10) << "Norms: "
@@ -235,7 +235,7 @@ void test_derefine_L2_element(int order, Element::Type el_type, int basis_type)
   }
 }

-TEST_CASE("AMR Coarsen L2 Element","[AMR][Coarsen]")
+TEST_CASE("AMR Coarsen L2 Element","[AMR][Coarsen][CUDA]")
 {
   std::vector<int> orders_1d{0,1,2,3};
   std::vector<int> orders_2d{0,1,2,3};
@@ -420,7 +420,7 @@ void stress_parallel_coarsen(int order, Element::Type el_type, int basis_type)
 #endif
 }

-TEST_CASE("Parallel AMR Coarsen Stress Test", "[AMR][Coarsen][Parallel]")
+TEST_CASE("Parallel AMR Coarsen Stress Test", "[AMR][Coarsen][Parallel][CUDA]")
 {
   std::vector<int> orders_1d{0,1,2,3};
   std::vector<int> orders_2d{0,1,2,3};
@@ -482,7 +482,7 @@ TEST_CASE("Parallel AMR Coarsen Stress Test", "[AMR][Coarsen][Parallel]")

 }

-TEST_CASE("ParDerefine", "[Parallel]")
+TEST_CASE("ParDerefine", "[Parallel][CUDA]")
 {
   for (dimension = 2; dimension <= 3; ++dimension)
   {