Merge master in lor-gpu

LOR H1/ND/RT explicit specializations
Add lor_h1 batched kernel for order 1
2022-10-23 09:54:23 -07:00 · 2022-09-15 07:28:55 -07:00 · 2022-09-13 17:26:27 -07:00 · 2022-09-13 12:56:41 -07:00 · 2022-09-08 16:42:02 -07:00 · 2022-09-08 16:33:48 -07:00
148 changed files with 4245 additions and 2808 deletions
@@ -62,7 +62,7 @@ jobs:

      - name: GHCR Login
        if: (github.event_name != 'pull_request')
-        uses: docker/login-action@v2
+        uses: docker/login-action@v1
        with:
          registry: ghcr.io
          username: ${{ github.actor }}
@@ -94,7 +94,7 @@ jobs:
    # This external action allows to interrupt a workflow already running on
    # the same branch to save resource
    - name: Cancel Previous Runs
-      uses: styfle/cancel-workflow-action@0.11.0
+      uses: styfle/cancel-workflow-action@0.9.0
      with:
        access_token: ${{ github.token }}

@@ -102,7 +102,7 @@ jobs:
    # /home/runner/work/mfem/mfem/mfem
    # Note: Done now to access "install-hypre" and "install-metis" actions.
    - name: checkout mfem
-      uses: actions/checkout@v3
+      uses: actions/checkout@v2
      with:
        path: ${{ env.MFEM_TOP_DIR }}
        # Fetch the complete history for codecov to access commits ID
@@ -115,25 +115,25 @@ jobs:
      if: matrix.mpi == 'par' && matrix.os == 'ubuntu-latest'
      run: |
        sudo apt-get install mpich libmpich-dev
+        export MAKE_CXX_FLAG="MPICXX=mpic++"

    - name: get lcov (Linux)
      if: matrix.codecov == 'YES' && matrix.os == 'ubuntu-latest'
      run: |
        sudo apt-get install lcov

-    # Keep the following section in case we need it again in the future,
-    # see: https://github.com/mfem/mfem/pull/3385#discussion_r1058013032
-    # - name: Set up Homebrew
-    #   if: ( matrix.mpi == 'par' || matrix.codecov == 'YES' ) && matrix.os == 'macos-latest'
-    #   uses: Homebrew/actions/setup-homebrew@master
+    - name: Set up Homebrew
+      if: ( matrix.mpi == 'par' || matrix.codecov == 'YES' ) && matrix.os == 'macos-latest'
+      uses: Homebrew/actions/setup-homebrew@c4aafe8c4620bf08883dd4679c374f11e73329d3

    - name: get MPI (MacOS)
      if: matrix.mpi == 'par' && matrix.os == 'macos-latest'
      run: |
        export HOMEBREW_NO_INSTALL_CLEANUP=1
        brew install openmpi
+        export MAKE_CXX_FLAG="MPICXX=mpic++"

-    - name: get lcov (MacOS)
+    - name: get MPI (MacOS)
      if: matrix.codecov == 'YES' && matrix.os == 'macos-latest'
      run: |
        export HOMEBREW_NO_INSTALL_CLEANUP=1
@@ -141,14 +141,14 @@ jobs:

    - name: get MPI (Windows)
      if: matrix.mpi == 'par' && matrix.os == 'windows-latest'
-      uses: mpi4py/setup-mpi@v1.1.2
+      uses: mpi4py/setup-mpi@v1.0.3

    # Get Hypre through cache, or build it.
    # Install will only run on cache miss.
    - name: cache hypre
      id: hypre-cache
      if: matrix.mpi == 'par'
-      uses: actions/cache@v3
+      uses: actions/cache@v2
      with:
        path: ${{ env.HYPRE_TOP_DIR }}
        key: ${{ runner.os }}-build-${{ env.HYPRE_TOP_DIR }}-${{ matrix.hypre-target }}-v2.2
@@ -176,7 +176,7 @@ jobs:
    - name: cache metis
      id: metis-cache
      if: matrix.mpi == 'par' && matrix.os != 'windows-latest'
-      uses: actions/cache@v3
+      uses: actions/cache@v2
      with:
        path: ${{ env.METIS_TOP_DIR }}
        key: ${{ runner.os }}-build-${{ env.METIS_TOP_DIR }}-v2.2
@@ -190,13 +190,12 @@ jobs:

    - name: cache vcpkg (Windows)
      id: vcpkg-cache
-      if: matrix.os == 'windows-latest'
      uses: actions/cache@v3
      with:
        path: vcpkg_cache
        key: ${{ runner.os }}-${{ matrix.mpi }}-vcpkg-v1

-    - name: prepare vcpkg binary cache location (Windows)
+    - name: prepare binary cache location
      if: matrix.os == 'windows-latest' && steps.vcpkg-cache.outputs.cache-hit != 'true'
      run: |
        mkdir -p vcpkg_cache
@@ -231,7 +230,11 @@ jobs:
      run: |
        cd ${{ env.MFEM_TOP_DIR }}  && make check

-    # Note: 'tests' include the unit tests
+    - name: unit tests
+      if: matrix.build-system == 'make' && matrix.target == 'opt'
+      run: |
+        cd ${{ env.MFEM_TOP_DIR }} && make unittest
+
    - name: tests
      if: matrix.build-system == 'make' && matrix.target == 'opt'
      run: |
@@ -244,7 +247,7 @@ jobs:
        cd ${{ env.MFEM_TOP_DIR }} && cmake --build build --target check --config ${CTEST_CONFIG}
      shell: bash

-    - name: cmake unit tests (Ubuntu)
+    - name: cmake unit tests (Ubuntu 20.04)
      if: matrix.build-system == 'cmake' && matrix.target == 'opt' && matrix.os == 'ubuntu-latest'
      run: |
        CTEST_CONFIG="Release"
@@ -262,7 +265,7 @@ jobs:
    # Code coverage (process and upload reports)
    - name: codecov
      if: matrix.codecov == 'YES'
-      uses: mfem/github-actions/upload-coverage@v2.2
+      uses: mfem/github-actions/upload-coverage@v2.0
      with:
        name: ${{ matrix.os }}-${{ matrix.build-system }}-${{ matrix.target }}-${{ matrix.mpi }}-${{ matrix.hypre-target }}
        project_dir: ${{ env.MFEM_TOP_DIR }}
@@ -35,22 +35,23 @@ jobs:

    steps:
    - name: Cancel Previous Runs
-      uses: styfle/cancel-workflow-action@0.11.0
+      uses: styfle/cancel-workflow-action@0.9.0
      with:
        access_token: ${{ github.token }}

    - name: checkout MFEM
-      uses: actions/checkout@v3
+      uses: actions/checkout@v2
      with:
        path: mfem

    - name: Get MPI (Linux)
      run: |
        sudo apt-get install mpich libmpich-dev
+        export MAKE_CXX_FLAG="MPICXX=mpic++"

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

    - name: Cache Metis Install
      id: metis-cache
-      uses: actions/cache@v3
+      uses: actions/cache@v2
      with:
        path: ${{ env.METIS_TOP_DIR }}
        key: ${{ runner.os }}-build-${{ env.METIS_TOP_DIR }}-v2.2
@@ -34,12 +34,12 @@ jobs:
       github.event.pull_request.head.repo.full_name != github.repository)
    steps:
    - name: Cancel Previous Runs
-      uses: styfle/cancel-workflow-action@0.11.0
+      uses: styfle/cancel-workflow-action@0.9.0
      with:
        access_token: ${{ github.token }}

    - name: checkout mfem
-      uses: actions/checkout@v3
+      uses: actions/checkout@v2

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

    - name: get astyle
      run: |
@@ -101,14 +101,11 @@ jobs:
       github.event.pull_request.head.repo.full_name != github.repository)
    steps:
    - name: checkout mfem
-      uses: actions/checkout@v3
+      uses: actions/checkout@v2

    - name: get doxygen and graphviz
      run: |
        sudo apt-get install doxygen graphviz
-
-    - name: update doxygen config file
-      run: |
        cd doc
        doxygen -u CodeDocumentation.conf.in

@@ -126,7 +123,7 @@ jobs:
    runs-on: ubuntu-latest
    steps:
    - name: checkout mfem
-      uses: actions/checkout@v3
+      uses: actions/checkout@v2
      with:
        fetch-depth: 0

@@ -18,6 +18,11 @@ CMakeFiles/
 # Backup files
 *~

+*.sqlite
+*.nsys-rep
+*.qdstrm
+*.csv
+
 # Default install location
 /mfem/

@@ -329,6 +334,7 @@ tests/benchmarks/bench_ceed
 tests/benchmarks/bench_tmop
 tests/benchmarks/bench_vector
 tests/benchmarks/bench_virtuals
+tests/benchmarks/bench_lor

 # Test script output
 tests/scripts/*.err
@@ -8,51 +8,6 @@
                               https://mfem.org


-
-
-Version 4.5.1 (development)
-===========================
-
-Meshing improvements
--------------------
- Added support for pyramids in non-conforming meshes. Currently only isotropic
-  refinement is supported in this case.
-
-Discretization improvements
---------------------------
- TBD
-
-Linear and nonlinear solvers
----------------------------
- Added a fast normalization-based distance solver, see the Distance miniapp
-  in the miniapps/shifted/ directory.
-
-New and updated examples and miniapps
-------------------------------------
- TBD
-
-Integrations, testing and documentation
---------------------------------------
- Removed the support for the Mesquite toolkit. We recommend using MFEM's TMOP
-  functionality instead for mesh optimization. See the mesh-optimizer miniapp.
-
-Miscellaneous
-------------
- VisItDataCollection now correctly handles collection names containing
-  underscores.
-
- VisItDataCollection::SetPadDigits() no longer alters the number of digits
-  used to represent the MPI rank because VisIt seems to require 6 digits.
-  This parameter can still be explicitly overridden with
-  VisItDataCollection::SetPadDigitsRank().
-
-API changes
-----------
- The implicit cast methods of class Vector to 'double *' and 'const double *'
-  have been deprecated and generate deprecation warnings if used. They will be
-  removed in a future release.
-
-
 Version 4.5, released on October 22, 2022
 =========================================

@@ -10,9 +10,7 @@
 # CONTRIBUTING.md for details.

 # The variable CMAKE_CXX_STANDARD and related were introduced in CMake v3.1
-# Version 3.8 fixes the handling of CMAKE_CXX_STANDARD for try_compile.
-# Version 3.8 or newer is required for direct CUDA support.
-cmake_minimum_required(VERSION 3.8)
+cmake_minimum_required(VERSION 3.1)
 set(USER_CONFIG "${CMAKE_CURRENT_SOURCE_DIR}/config/user.cmake" CACHE PATH
  "Path to optional user configuration file.")

@@ -53,7 +51,7 @@ project(mfem NONE)
 # Current version of MFEM, see also `makefile`.
 #   mfem_VERSION = (string)
 #   MFEM_VERSION = (int)   [automatically derived from mfem_VERSION]
-set(${PROJECT_NAME}_VERSION 4.5.1)
+set(${PROJECT_NAME}_VERSION 4.5.0)

 # Prohibit in-source build
 if (${PROJECT_SOURCE_DIR} STREQUAL ${PROJECT_BINARY_DIR})
@@ -68,7 +66,8 @@ set(XSDK_ENABLE_C OFF)
 set(XSDK_ENABLE_Fortran OFF)

 # Check if we need to enable C or Fortran.
-if (MFEM_USE_CONDUIT OR
+if (CMAKE_VERSION VERSION_LESS 3.2 OR
+    MFEM_USE_CONDUIT OR
    MFEM_USE_SIDRE OR
    MFEM_USE_PETSC)
   # This seems to be needed by:
@@ -82,9 +81,11 @@ if (MFEM_USE_STRUMPACK)
  # Just needed to find the MPI_Fortran libraries to link with
  set(XSDK_ENABLE_Fortran ON)
 endif()
-# SUNDIALS, STRUMPACK, and Ginkgo require C++14:
-if ((MFEM_USE_SUNDIALS OR MFEM_USE_STRUMPACK OR MFEM_USE_GINKGO) AND
-    ("${CMAKE_CXX_STANDARD}" LESS "14"))
+# SUNDIALS >= 6.4.0 requires C++14:
+if (MFEM_USE_SUNDIALS AND ("${CMAKE_CXX_STANDARD}" LESS "14"))
+  set(CMAKE_CXX_STANDARD 14)
+endif()
+if (MFEM_USE_GINKGO AND ("${CMAKE_CXX_STANDARD}" LESS "14"))
  set(CMAKE_CXX_STANDARD 14)
 endif()

@@ -102,6 +103,8 @@ if (MFEM_USE_CUDA)
   if (MFEM_USE_HIP)
      message(FATAL_ERROR " *** MFEM_USE_HIP cannot be combined with MFEM_USE_CUDA.")
   endif()
+   # MFEM_USE_CUDA requires CMake 3.8 or newer (for direct CUDA support)
+   cmake_minimum_required(VERSION 3.8 FATAL_ERROR)
   # Use ${CMAKE_CXX_COMPILER} as the cuda host compiler.
   if (NOT CMAKE_CUDA_HOST_COMPILER)
      set(CMAKE_CUDA_HOST_COMPILER ${CMAKE_CXX_COMPILER})
@@ -196,26 +199,6 @@ if (MFEM_USE_HIP)
  find_package(HIPSPARSE REQUIRED)
 endif()

-# OpenMP
-if (MFEM_USE_OPENMP OR MFEM_USE_LEGACY_OPENMP)
-  if (NOT MFEM_THREAD_SAFE AND MFEM_USE_LEGACY_OPENMP)
-    message(FATAL_ERROR " *** MFEM_USE_LEGACY_OPENMP requires MFEM_THREAD_SAFE=ON.")
-  endif()
-  find_package(OpenMP REQUIRED)
-  set(OPENMP_LIBRARIES ${OpenMP_CXX_LIBRARIES})
-  if(APPLE)
-    # On macOS, the compiler needs additional help to find the <omp.h> header.
-    # See issue #2642 for more information.
-    set(OPENMP_INCLUDE_DIRS ${OpenMP_CXX_INCLUDE_DIRS})
-  endif(APPLE)
-  if (OPENMP_FOUND)
-    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
-    if (MFEM_USE_CUDA)
-      set(CMAKE_CUDA_FLAGS "${CMAKE_CUDA_FLAGS} -Xcompiler=${OpenMP_CXX_FLAGS}")
-    endif()
-  endif()
-endif()
-
 # MPI -> hypre; PETSc (optional)
 if (MFEM_USE_MPI)
  find_package(MPI REQUIRED)
@@ -281,6 +264,20 @@ if (MFEM_USE_LAPACK)
  find_package(LAPACK REQUIRED)
 endif()

+# OpenMP
+if (MFEM_USE_OPENMP OR MFEM_USE_LEGACY_OPENMP)
+  if (NOT MFEM_THREAD_SAFE AND MFEM_USE_LEGACY_OPENMP)
+    message(FATAL_ERROR " *** MFEM_USE_LEGACY_OPENMP requires MFEM_THREAD_SAFE=ON.")
+  endif()
+  find_package(OpenMP REQUIRED)
+  set(OPENMP_LIBRARIES ${OpenMP_CXX_LIBRARIES})
+  if(APPLE)
+    # On macOS, the compiler needs additional help to find the <omp.h> header.
+    # See issue #2642 for more information.
+    set(OPENMP_INCLUDE_DIRS ${OpenMP_CXX_INCLUDE_DIRS})
+  endif(APPLE)
+endif()
+
 # SuiteSparse (before SUNDIALS which may depend on KLU)
 if (MFEM_USE_SUITESPARSE)
  find_package(SuiteSparse REQUIRED
@@ -299,6 +296,11 @@ if (MFEM_USE_SUNDIALS)
  find_package(SUNDIALS REQUIRED ${SUNDIALS_COMPONENTS})
 endif()

+# Mesquite
+if (MFEM_USE_MESQUITE)
+  find_package(Mesquite REQUIRED)
+endif()
+
 # SuperLU_DIST can only be enabled in parallel
 if (MFEM_USE_SUPERLU)
  if (MFEM_USE_MPI)
@@ -486,10 +488,6 @@ endif()
 #   an ALIAS target is missing?
 # Call Stack (most recent call first):
 #   CMakeLists.txt:474 (mfem_add_library)
-#
-# NOTE: We need to figure out which TPL library adds the dependency on
-#       "Threads::Threads" and call the next line only when that TPL library is
-#       enabled.  -V. Dobrev
 find_package(Threads REQUIRED)

 # List all possible libraries in order of dependencies.
@@ -497,8 +495,8 @@ find_package(Threads REQUIRED)
 #    With newer versions of SuiteSparse which include METIS header using 64-bit
 #    integers, the METIS header (with 32-bit indices, as used by mfem) needs to
 #    be before SuiteSparse.
-set(MFEM_TPLS OPENMP HYPRE LAPACK BLAS SuperLUDist STRUMPACK METIS SuiteSparse
-    SUNDIALS PETSC SLEPC MUMPS AXOM FMS CONDUIT Ginkgo GNUTLS GSLIB
+set(MFEM_TPLS OPENMP HYPRE LAPACK BLAS SuperLUDist METIS SuiteSparse SUNDIALS
+    PETSC SLEPC MESQUITE MUMPS STRUMPACK AXOM FMS CONDUIT Ginkgo GNUTLS GSLIB
    NETCDF MPFR PUMI HIOP POSIXCLOCKS MFEMBacktrace ZLIB OCCA CEED RAJA UMPIRE
    ADIOS2 CUBLAS CUSPARSE MKL_CPARDISO AMGX CALIPER CODIPACK BENCHMARK PARELAG
    MPI_CXX HIP HIPSPARSE MOONOLITH BLITZ ALGOIM ENZYME)
@@ -513,12 +511,18 @@ foreach(TPL IN LISTS MFEM_TPLS)
    list(APPEND TPL_INCLUDE_DIRS ${${TPL}_INCLUDE_DIRS})
  endif()
 endforeach(TPL)
-list(REVERSE TPL_LIBRARIES)
 list(REMOVE_DUPLICATES TPL_LIBRARIES)
-list(REVERSE TPL_LIBRARIES)
 list(REMOVE_DUPLICATES TPL_INCLUDE_DIRS)
 # message(STATUS "TPL_INCLUDE_DIRS = ${TPL_INCLUDE_DIRS}")

+if (OPENMP_FOUND)
+  message(STATUS "MFEM: using package OpenMP")
+  set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
+  if (MFEM_USE_CUDA)
+    set(CMAKE_CUDA_FLAGS "${CMAKE_CUDA_FLAGS} -Xcompiler=${OpenMP_CXX_FLAGS}")
+  endif()
+endif()
+
 message(STATUS "MFEM build type: CMAKE_BUILD_TYPE = ${CMAKE_BUILD_TYPE}")
 message(STATUS "MFEM version: v${MFEM_VERSION_STRING}")
 message(STATUS "MFEM git string: ${MFEM_GIT_STRING}")
@@ -337,6 +337,10 @@ MFEM_USE_SUNDIALS = YES/NO
   library. When enabled, this option uses the SUNDIALS_* library options,
   see below.

+MFEM_USE_MESQUITE = YES/NO
+   Enable MFEM functionality based on the Mesquite library. When enabled, this
+   option uses the MESQUITE_* library options, see below.
+
 MFEM_USE_SUITESPARSE = YES/NO
   Enable MFEM functionality based on the SuiteSparse library. Currently, this
   option adds the classes UMFPackSolver and KLUSolver (both sparse serial
@@ -632,6 +636,11 @@ The specific libraries and their options are:
  Options: SUNDIALS_OPT, SUNDIALS_LIB.
  Versions: SUNDIALS >= 5.0.0, SUNDIALS >= 5.4.0 for CUDA support.

+- Mesquite (optional), used when MFEM_USE_MESQUITE = YES.
+  URL: http://trilinos.org/oldsite/packages/mesquite
+  Options: MESQUITE_OPT, MESQUITE_LIB.
+  The Mesquite support is deprecated and will be removed in the future.
+
 - SuiteSparse (optional), used when MFEM_USE_SUITESPARSE = YES.
  URL: http://faculty.cse.tamu.edu/davis/suitesparse.html
  Options: SUITESPARSE_OPT, SUITESPARSE_LIB.
@@ -960,6 +969,7 @@ MFEM_USE_LEGACY_OPENMP
 MFEM_USE_OPENMP
 MFEM_USE_MEMALLOC
 MFEM_TIMER_TYPE - Set automatically, can be overwritten.
+MFEM_USE_MESQUITE
 MFEM_USE_SUITESPARSE
 MFEM_USE_SUPERLU
 MFEM_USE_MUMPS
@@ -1024,6 +1034,7 @@ The CMake build system adds auto-detection for the following packages/libraries:
 - HYPRE
 - METIS - The option MFEM_USE_METIS_5 is auto-detected.
 - ParMETIS
+ - MESQUITE
 - SuiteSparse
 - SuperLUDist, STRUMPACK
 - Ginkgo
@@ -212,6 +212,10 @@ IF (DEFINED TPL_ENABLE_SUNDIALS)
  SET(MFEM_USE_SUNDIALS ${TPL_ENABLE_SUNDIALS} CACHE BOOL "Enable SUNDIALS usage" FORCE)
 ENDIF()

+IF (DEFINED TPL_ENABLE_MESQUITE)
+  SET(MFEM_USE_MESQUITE ${TPL_ENABLE_MESQUITE} CACHE BOOL "Enable MESQUITE usage" FORCE)
+ENDIF()
+
 IF (DEFINED TPL_ENABLE_SUITESPARSE)
  SET(MFEM_USE_SUITESPARSE ${TPL_ENABLE_SUITESPARSE} CACHE BOOL "Enable SuiteSparse usage" FORCE)
 ENDIF()
@@ -29,6 +29,7 @@ set(MFEM_USE_LEGACY_OPENMP @MFEM_USE_LEGACY_OPENMP@)
 set(MFEM_USE_MEMALLOC @MFEM_USE_MEMALLOC@)
 set(MFEM_TIMER_TYPE @MFEM_TIMER_TYPE@)
 set(MFEM_USE_SUNDIALS @MFEM_USE_SUNDIALS@)
+set(MFEM_USE_MESQUITE @MFEM_USE_MESQUITE@)
 set(MFEM_USE_SUITESPARSE @MFEM_USE_SUITESPARSE@)
 set(MFEM_USE_SUPERLU @MFEM_USE_SUPERLU@)
 set(MFEM_USE_MUMPS @MFEM_USE_MUMPS@)
@@ -77,6 +77,9 @@
 // Internal MFEM option: enable group/batch allocation for some small objects.
 #cmakedefine MFEM_USE_MEMALLOC

+// Enable MFEM functionality based on the Mesquite library.
+#cmakedefine MFEM_USE_MESQUITE
+
 // Enable MFEM functionality based on the SuiteSparse library.
 #cmakedefine MFEM_USE_SUITESPARSE

@@ -16,22 +16,7 @@

 include(MfemCmakeUtilities)
 mfem_find_package(Caliper CALIPER CALIPER_DIR
-      "include" "caliper/cali.h"
-      "lib" "caliper"
-      "Paths to headers required by Caliper."
-      "Libraries required by Caliper.")
-
-# Append adiak path/lib if the user provided ADIAK_DIR
-if(ADIAK_DIR AND EXISTS ${ADIAK_DIR})
-    find_package(adiak NO_DEFAULT_PATH REQUIRED PATHS ${ADIAK_DIR}/lib/cmake/adiak ${ADIAK_DIR})
-    list(APPEND CALIPER_INCLUDE_DIRS ${adiak_INCLUDE_DIRS})
-    list(APPEND CALIPER_LIBRARIES ${adiak_LIBRARIES})
-endif()
-
-# Append gotcha path/lib if the user provided GOTCHA_DIR
-if(GOTCHA_DIR AND EXISTS ${GOTCHA_DIR})
-    find_package(gotcha NO_DEFAULT_PATH REQUIRED PATHS ${GOTCHA_DIR}/lib/cmake/gotcha ${GOTCHA_DIR})
-    list(APPEND CALIPER_INCLUDE_DIRS ${gotcha_INCLUDE_DIRS})
-    list(APPEND CALIPER_LIBRARIES ${gotcha_LIBRARIES})
-endif()
-
+       	"include" "caliper/cali.h"
+       	"lib" "caliper"
+  	"Paths to headers required by Caliper."
+       	"Libraries required by Caliper.")
@@ -14,21 +14,9 @@
 #   - HDF5_LIBRARIES - The HDF5 libraries
 #   - HDF5_INCLUDE_DIRS - The HDF5 include directories

-# NOTE: Using this FindHDF5.cmake instead of the CMake provided version may lead
-#       to issues with some TPL libraries that depend (or may depend) on HDF5.
-#       For this reason, we should consider removing this file, or at least
-#       making it use the CMake provided version by default and apply the logic
-#       below only when specifically requested by a user.  -V. Dobrev
-
 # First Check for HDF5_DIR
 if(NOT HDF5_DIR)
-  message(FATAL_ERROR
-    "Could not find HDF5. HDF5 support needs explicit HDF5_DIR")
-endif()
-
-if (NOT HDF5_FIND_QUIETLY)
-  message(STATUS "Looking for HDF5 ...")
-  message(STATUS "   in HDF5_DIR = ${HDF5_DIR}")
+    MESSAGE(FATAL_ERROR "Could not find HDF5. HDF5 support needs explicit HDF5_DIR")
 endif()

 # Find includes
@@ -62,9 +50,8 @@ include(FindPackageHandleStandardArgs)

 # Handle the QUIETLY and REQUIRED arguments and set HDF5_FOUND to TRUE if all
 # listed variables are TRUE
-find_package_handle_standard_args(HDF5
-                                  " *** HDF5 not found. Please set HDF5_DIR."
-                                  HDF5_LIBRARIES
+find_package_handle_standard_args(HDF5  DEFAULT_MSG
                                  HDF5_INCLUDE_DIRS
                                  __HDF5_LIBRARY
-                                  __HDF5_HL_LIBRARY)
+                                  __HDF5_HL_LIBRARY
+                                  HDF5_LIBRARIES )
@@ -0,0 +1,20 @@
+# Copyright (c) 2010-2022, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+# Defines the following variables:
+#   - MESQUITE_FOUND
+#   - MESQUITE_LIBRARIES
+#   - MESQUITE_INCLUDE_DIRS
+
+include(MfemCmakeUtilities)
+mfem_find_package(Mesquite MESQUITE MESQUITE_DIR
+  "include" "Mesquite_all_headers.hpp" "lib" "mesquite"
+  "Paths to headers required by Mesquite." "Libraries required by Mesquite.")
@@ -17,24 +17,18 @@
 include(MfemCmakeUtilities)

 # FindHDF5.cmake uses HDF5_ROOT, so we "translate" from the MFEM convention
-# (MFEM's FindHDF5.cmake does not need HDF5_ROOT)
-# set(HDF5_ROOT ${HDF5_DIR} CACHE PATH "")
-
+set(HDF5_ROOT ${HDF5_DIR} CACHE PATH "")
 # We need to guard against the case where HDF5 was already found but without
 # the HL extensions (in which case mfem_find_package will treat the package
 # as already having been found), so we reset the variable to force FindHDF5.cmake
 # to be called for a second time
 set(HDF5_FOUND OFF)
 enable_language(C) # FindHDF5.cmake uses the C compiler
-
 mfem_find_package(NetCDF NETCDF NETCDF_DIR "include" netcdf.h "lib" netcdf
  "Paths to headers required by NetCDF." "Libraries required by NetCDF.")
-
 # The HL extension libraries are in a separate variable and must precede
 # the "regular" hdf5 library, as hdf5_hl depends on hdf5
 # The netcdf library will always be the first element of NETCDF_LIBRARIES
 # and we need to insert after that library but before the hdf5 library, so
 # position 1 is used
-# (MFEM's FindHDF5.cmake does not set HDF5_C_LIBRARY_hdf5_hl and the HL library
-# is already added to NETCDF_LIBRARIES)
-# list(INSERT NETCDF_LIBRARIES 1 ${HDF5_C_LIBRARY_hdf5_hl})
+list(INSERT NETCDF_LIBRARIES 1 ${HDF5_C_LIBRARY_hdf5_hl})
@@ -46,10 +46,6 @@ endfunction()
 # Wrapper for add_executable
 macro(mfem_add_executable NAME)
  add_executable(${NAME} ${ARGN})
-  if (MFEM_USE_CUDA)
-    set_target_properties(${NAME} PROPERTIES
-      CUDA_RESOLVE_DEVICE_SYMBOLS ON)
-  endif()
 endmacro()

 # Wrapper for add_library
@@ -162,12 +158,27 @@ macro(add_mfem_miniapp MFEM_EXE_NAME)

  # Append the additional libraries and options
  if (LIBRARIES_LIST)
-    target_link_libraries(${MFEM_EXE_NAME} PRIVATE ${LIBRARIES_LIST})
+    if(CMAKE_VERSION VERSION_GREATER 2.8.11)
+      target_link_libraries(${MFEM_EXE_NAME} PRIVATE ${LIBRARIES_LIST})
+    else()
+      target_link_libraries(${MFEM_EXE_NAME} ${LIBRARIES_LIST})
+    endif()
  endif()
  if (EXTRA_OPTIONS_LIST)
    string(REPLACE ";" " " EXTRA_OPTIONS_STRING "${EXTRA_OPTIONS_LIST}")
    message(STATUS "${MFEM_EXE_NAME}: add flags \"${EXTRA_OPTIONS_STRING}\"")
-    target_compile_options(${MFEM_EXE_NAME} PRIVATE ${EXTRA_OPTIONS_LIST})
+    if(CMAKE_VERSION VERSION_GREATER 2.8.11)
+      target_compile_options(${MFEM_EXE_NAME} PRIVATE ${EXTRA_OPTIONS_LIST})
+    else()
+      get_target_property(THIS_COMPILE_FLAGS ${MFEM_EXE_NAME} COMPILE_FLAGS)
+      if (THIS_COMPILE_FLAGS)
+        set(THIS_COMPILE_FLAGS "${THIS_COMPILE_FLAGS} ${EXTRA_OPTIONS_STRING}")
+      else()
+        set(THIS_COMPILE_FLAGS "${EXTRA_OPTIONS_STRING}")
+      endif()
+      set_target_properties(${MFEM_EXE_NAME}
+        PROPERTIES COMPILE_FLAGS ${THIS_COMPILE_FLAGS})
+    endif()
  endif()
  if (EXTRA_DEFINES_LIST)
    target_compile_definitions(${MFEM_EXE_NAME} PRIVATE ${EXTRA_DEFINES_LIST})
@@ -176,15 +187,17 @@ macro(add_mfem_miniapp MFEM_EXE_NAME)
  # Handle the MPI separately
  if (MFEM_USE_MPI)
    # Add MPI_CXX_LIBRARIES, in case this target does not link with mfem.
-    target_link_libraries(${MFEM_EXE_NAME} PRIVATE ${MPI_CXX_LIBRARIES})
+    if(CMAKE_VERSION VERSION_GREATER 2.8.11)
+      target_link_libraries(${MFEM_EXE_NAME} PRIVATE ${MPI_CXX_LIBRARIES})
+    else()
+      target_link_libraries(${MFEM_EXE_NAME} ${MPI_CXX_LIBRARIES})
+    endif()

    if (MPI_CXX_INCLUDE_PATH)
      target_include_directories(${MFEM_EXE_NAME} PRIVATE ${MPI_CXX_INCLUDE_PATH})
    endif()
    if (MPI_CXX_COMPILE_FLAGS)
-      separate_arguments(MPI_CXX_COMPILE_ARGS UNIX_COMMAND
-        "${MPI_CXX_COMPILE_FLAGS}")
-      target_compile_options(${MFEM_EXE_NAME} PRIVATE ${MPI_CXX_COMPILE_ARGS})
+      target_compile_options(${MFEM_EXE_NAME} PRIVATE ${MPI_CXX_COMPILE_FLAGS})
    endif()

    if (MPI_CXX_LINK_FLAGS)
@@ -865,7 +878,7 @@ function(mfem_export_mk_files)
  set(CONFIG_MK_BOOL_VARS MFEM_USE_MPI MFEM_USE_METIS MFEM_USE_METIS_5
      MFEM_DEBUG MFEM_USE_EXCEPTIONS MFEM_USE_ZLIB MFEM_USE_LIBUNWIND
      MFEM_USE_LAPACK MFEM_THREAD_SAFE MFEM_USE_LEGACY_OPENMP MFEM_USE_OPENMP
-      MFEM_USE_MEMALLOC MFEM_USE_SUNDIALS MFEM_USE_SUITESPARSE
+      MFEM_USE_MEMALLOC MFEM_USE_SUNDIALS MFEM_USE_MESQUITE MFEM_USE_SUITESPARSE
      MFEM_USE_SUPERLU MFEM_USE_SUPERLU5 MFEM_USE_MUMPS MFEM_USE_STRUMPACK
      MFEM_USE_GINKGO MFEM_USE_AMGX MFEM_USE_GNUTLS MFEM_USE_NETCDF
      MFEM_USE_PETSC MFEM_USE_SLEPC MFEM_USE_MPFR MFEM_USE_SIDRE MFEM_USE_FMS
@@ -85,6 +85,9 @@
 // Enable MFEM functionality based on the SUNDIALS libraries.
 // #define MFEM_USE_SUNDIALS

+// Enable MFEM functionality based on the Mesquite library.
+// #define MFEM_USE_MESQUITE
+
 // Enable MFEM functionality based on the SuiteSparse library.
 // #define MFEM_USE_SUITESPARSE

@@ -29,6 +29,7 @@ MFEM_USE_OPENMP        = @MFEM_USE_OPENMP@
 MFEM_USE_MEMALLOC      = @MFEM_USE_MEMALLOC@
 MFEM_TIMER_TYPE        = @MFEM_TIMER_TYPE@
 MFEM_USE_SUNDIALS      = @MFEM_USE_SUNDIALS@
+MFEM_USE_MESQUITE      = @MFEM_USE_MESQUITE@
 MFEM_USE_SUITESPARSE   = @MFEM_USE_SUITESPARSE@
 MFEM_USE_SUPERLU       = @MFEM_USE_SUPERLU@
 MFEM_USE_SUPERLU5      = @MFEM_USE_SUPERLU5@
@@ -30,6 +30,7 @@ option(MFEM_USE_OPENMP "Enable the OpenMP backend" OFF)
 option(MFEM_USE_LEGACY_OPENMP "Enable legacy OpenMP usage" OFF)
 option(MFEM_USE_MEMALLOC "Enable the internal MEMALLOC option." ON)
 option(MFEM_USE_SUNDIALS "Enable SUNDIALS usage" OFF)
+option(MFEM_USE_MESQUITE "Enable MESQUITE usage" OFF)
 option(MFEM_USE_SUITESPARSE "Enable SuiteSparse usage" OFF)
 option(MFEM_USE_SUPERLU "Enable SuperLU_DIST usage" OFF)
 option(MFEM_USE_SUPERLU5 "Use the old SuperLU_DIST 5.1 version" OFF)
@@ -123,6 +124,9 @@ set(SUNDIALS_DIR "${MFEM_DIR}/../sundials-5.0.0/instdir" CACHE PATH
 # set(SUNDIALS_REQUIRED_PACKAGES "SuiteSparse/KLU/AMD/BTF/COLAMD/config"
 #     CACHE STRING "Additional packages required by SUNDIALS.")

+set(MESQUITE_DIR "${MFEM_DIR}/../mesquite-2.99" CACHE PATH
+    "Path to the Mesquite library.")
+
 set(SuiteSparse_DIR "${MFEM_DIR}/../SuiteSparse" CACHE PATH
    "Path to the SuiteSparse library.")
 set(SuiteSparse_REQUIRED_PACKAGES "BLAS" "METIS"
@@ -184,7 +188,6 @@ set(GNUTLS_DIR "" CACHE PATH "Path to the GnuTLS library.")

 set(GSLIB_DIR "" CACHE PATH "Path to the GSLIB library.")

-set(HDF5_DIR "/usr" CACHE PATH "Path to the HDF5 library.")
 set(NETCDF_DIR "" CACHE PATH "Path to the NetCDF library.")
 set(NetCDF_REQUIRED_PACKAGES "HDF5/C/HL" CACHE STRING
    "Additional packages required by NetCDF.")
@@ -131,6 +131,7 @@ MFEM_USE_LEGACY_OPENMP = NO
 MFEM_USE_MEMALLOC      = YES
 MFEM_TIMER_TYPE        = $(if $(NOTMAC),2,4)
 MFEM_USE_SUNDIALS      = NO
+MFEM_USE_MESQUITE      = NO
 MFEM_USE_SUITESPARSE   = NO
 MFEM_USE_SUPERLU       = NO
 MFEM_USE_SUPERLU5      = NO
@@ -269,6 +270,11 @@ endif
 # If SUNDIALS was built with KLU:
 # MFEM_USE_SUITESPARSE = YES

+# MESQUITE library configuration
+MESQUITE_DIR = @MFEM_DIR@/../mesquite-2.99
+MESQUITE_OPT = -I$(MESQUITE_DIR)/include
+MESQUITE_LIB = -L$(MESQUITE_DIR)/lib -lmesquite
+
 # SuiteSparse library configuration
 LIB_RT = $(if $(NOTMAC),-lrt,)
 SUITESPARSE_DIR = @MFEM_DIR@/../SuiteSparse
@@ -318,9 +324,6 @@ MUMPS_LIB = $(XLINKER)-rpath,$(MUMPS_DIR)/lib -L$(MUMPS_DIR)/lib -ldmumps\

 # STRUMPACK library configuration
 STRUMPACK_DIR = @MFEM_DIR@/../STRUMPACK-build
-ifeq ($(MFEM_USE_STRUMPACK),YES)
-   BASE_FLAGS = -std=c++14
-endif
 STRUMPACK_OPT = -I$(STRUMPACK_DIR)/include $(SCOTCH_OPT)
 # If STRUMPACK was build with OpenMP support, the following may be need:
 # STRUMPACK_OPT += $(OPENMP_OPT)
@@ -473,17 +476,7 @@ OCCA_LIB = $(XLINKER)-rpath,$(OCCA_DIR)/lib -L$(OCCA_DIR)/lib -locca
 # CALIPER library configuration
 CALIPER_DIR = @MFEM_DIR@/../caliper
 CALIPER_OPT = -I$(CALIPER_DIR)/include
-CALIPER_LIB = $(XLINKER)-rpath,$(CALIPER_DIR)/lib64 $(XLINKER)-rpath,$(CALIPER_DIR)/lib -L$(CALIPER_DIR)/lib64 -L$(CALIPER_DIR)/lib -lcaliper
-
-ifdef ADIAK_DIR
-   CALIPER_OPT += -I$(ADIAK_DIR)/include
-   CALIPER_LIB += $(XLINKER)-rpath,$(ADIAK_DIR)/lib64 $(XLINKER)-rpath,$(ADIAK_DIR)/lib -L$(ADIAK_DIR)/lib64 -L$(ADIAK_DIR)/lib -ladiak
-endif
-ifdef GOTCHA_DIR
-   CALIPER_OPT += -I$(GOTCHA_DIR)/include
-   CALIPER_LIB += $(XLINKER)-rpath,$(GOTCHA_DIR)/lib64 $(XLINKER)-rpath,$(GOTCHA_DIR)/lib -L$(GOTCHA_DIR)/lib64 -L$(GOTCHA_DIR)/lib -lgotcha
-endif
-
+CALIPER_LIB = $(XLINKER)-rpath,$(CALIPER_DIR)/lib64 -L$(CALIPER_DIR)/lib64 -lcaliper

 # BLITZ library configuration
 BLITZ_DIR = @MFEM_DIR@/../blitz
@@ -38,7 +38,7 @@ PROJECT_NAME           = "MFEM"
 # could be handy for archiving the generated documentation or if some version
 # control system is used.

-PROJECT_NUMBER         = v4.5.1
+PROJECT_NUMBER         = v4.5.0

 # Using the PROJECT_BRIEF tag one can provide an optional one line description
 # for a project that appears at the top of each page and should give viewer a
@@ -763,54 +763,36 @@ INPUT                  = @MFEM_SOURCE_DIR@/doc/CodeDocumentation.dox \
                         @MFEM_SOURCE_DIR@/config \
                         @MFEM_SOURCE_DIR@/general \
                         @MFEM_SOURCE_DIR@/linalg \
-                         @MFEM_SOURCE_DIR@/linalg/simd \
                         @MFEM_SOURCE_DIR@/mesh \
-                         @MFEM_SOURCE_DIR@/mesh/submesh \
                         @MFEM_SOURCE_DIR@/fem \
-                         @MFEM_SOURCE_DIR@/fem/ceed \
-                         @MFEM_SOURCE_DIR@/fem/ceed/integrators \
-                         @MFEM_SOURCE_DIR@/fem/ceed/integrators/convection \
-                         @MFEM_SOURCE_DIR@/fem/ceed/integrators/diffusion \
-                         @MFEM_SOURCE_DIR@/fem/ceed/integrators/mass \
-                         @MFEM_SOURCE_DIR@/fem/ceed/integrators/nlconvection \
-                         @MFEM_SOURCE_DIR@/fem/ceed/interface \
-                         @MFEM_SOURCE_DIR@/fem/ceed/solvers \
+                         @MFEM_SOURCE_DIR@/fem/moonolith \
                         @MFEM_SOURCE_DIR@/fem/fe \
                         @MFEM_SOURCE_DIR@/fem/lor \
-                         @MFEM_SOURCE_DIR@/fem/moonolith \
-                         @MFEM_SOURCE_DIR@/fem/qinterp \
-                         @MFEM_SOURCE_DIR@/fem/tmop \
                         @MFEM_SOURCE_DIR@/examples \
-                         @MFEM_SOURCE_DIR@/examples/amgx \
                         @MFEM_SOURCE_DIR@/examples/caliper \
+                         @MFEM_SOURCE_DIR@/examples/amgx \
                         @MFEM_SOURCE_DIR@/examples/ginkgo \
-                         @MFEM_SOURCE_DIR@/examples/hiop \
                         @MFEM_SOURCE_DIR@/examples/moonolith \
+                         @MFEM_SOURCE_DIR@/examples/hiop \
                         @MFEM_SOURCE_DIR@/examples/petsc \
                         @MFEM_SOURCE_DIR@/examples/pumi \
                         @MFEM_SOURCE_DIR@/examples/sundials \
                         @MFEM_SOURCE_DIR@/examples/superlu \
                         @MFEM_SOURCE_DIR@/miniapps/adjoint \
-                         @MFEM_SOURCE_DIR@/miniapps/autodiff \
                         @MFEM_SOURCE_DIR@/miniapps/common \
                         @MFEM_SOURCE_DIR@/miniapps/electromagnetics \
                         @MFEM_SOURCE_DIR@/miniapps/gslib \
-                         @MFEM_SOURCE_DIR@/miniapps/hooke \
-                         @MFEM_SOURCE_DIR@/miniapps/hooke/kernels \
-                         @MFEM_SOURCE_DIR@/miniapps/hooke/materials \
-                         @MFEM_SOURCE_DIR@/miniapps/hooke/operators \
-                         @MFEM_SOURCE_DIR@/miniapps/hooke/preconditioners \
                         @MFEM_SOURCE_DIR@/miniapps/meshing \
                         @MFEM_SOURCE_DIR@/miniapps/mtop \
-                         @MFEM_SOURCE_DIR@/miniapps/multidomain \
+                         @MFEM_SOURCE_DIR@/miniapps/autodiff \
                         @MFEM_SOURCE_DIR@/miniapps/navier \
                         @MFEM_SOURCE_DIR@/miniapps/nurbs \
-                         @MFEM_SOURCE_DIR@/miniapps/parelag \
                         @MFEM_SOURCE_DIR@/miniapps/performance \
                         @MFEM_SOURCE_DIR@/miniapps/shifted \
                         @MFEM_SOURCE_DIR@/miniapps/solvers \
                         @MFEM_SOURCE_DIR@/miniapps/tools \
-                         @MFEM_SOURCE_DIR@/miniapps/toys
+                         @MFEM_SOURCE_DIR@/miniapps/toys \
+                         @MFEM_SOURCE_DIR@/miniapps/parelag

 # This tag can be used to specify the character encoding of the source files
 # that doxygen parses. Internally doxygen uses the UTF-8 encoding. Doxygen uses
@@ -33,13 +33,13 @@ add_mfem_examples(CALIPER_EXE_SRCS ${PREFIX})
 if (MFEM_ENABLE_TESTING)
  foreach(SRC_FILE ${CALIPER_EXE_SRCS})
    get_filename_component(SRC_FILENAME ${SRC_FILE} NAME)
-    string(REPLACE ".cpp" "" TEST_NAME ${PREFIX}${SRC_FILENAME})
+    string(REPLACE ".cpp" "" TEST_NAME ${SRC_FILENAME})

    set(THIS_TEST_OPTIONS "-no-vis")

    if (NOT (${TEST_NAME} MATCHES ".*p$"))
      add_test(NAME ${TEST_NAME}_ser
-        COMMAND $<TARGET_FILE:${TEST_NAME}> ${THIS_TEST_OPTIONS})
+        COMMAND ${TEST_NAME} ${THIS_TEST_OPTIONS})
    else()
      add_test(NAME ${TEST_NAME}_np=${MFEM_MPI_NP}
        COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
@@ -30,7 +30,7 @@
 //
 // Device sample runs:
 //               ex1 -pa -d cuda
-//               ex1 -fa -d cuda
+//             * ex1 -fa -d cuda
 //               ex1 -pa -d raja-cuda
 //             * ex1 -pa -d raja-hip
 //               ex1 -pa -d occa-cuda
@@ -192,14 +192,7 @@ int main(int argc, char *argv[])
   //    domain integrator.
   BilinearForm a(&fespace);
   if (pa) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
-   if (fa)
-   {
-      a.SetAssemblyLevel(AssemblyLevel::FULL);
-      // Sort the matrix column indices when running on GPU or with OpenMP (i.e.
-      // when Device::IsEnabled() returns true). This makes the results
-      // bit-for-bit deterministic at the cost of somewhat longer run time.
-      a.EnableSparseMatrixSorting(Device::IsEnabled());
-   }
+   if (fa) { a.SetAssemblyLevel(AssemblyLevel::FULL); }
   a.AddDomainIntegrator(new DiffusionIntegrator(one));

   // 10. Assemble the bilinear form and the corresponding linear system,
@@ -30,7 +30,7 @@
 //
 // Device sample runs:
 //               mpirun -np 4 ex1p -pa -d cuda
-//               mpirun -np 4 ex1p -fa -d cuda
+//             * mpirun -np 4 ex1p -fa -d cuda
 //               mpirun -np 4 ex1p -pa -d occa-cuda
 //               mpirun -np 4 ex1p -pa -d raja-omp
 //               mpirun -np 4 ex1p -pa -d ceed-cpu
@@ -219,14 +219,7 @@ int main(int argc, char *argv[])
   //     Diffusion domain integrator.
   ParBilinearForm a(&fespace);
   if (pa) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
-   if (fa)
-   {
-      a.SetAssemblyLevel(AssemblyLevel::FULL);
-      // Sort the matrix column indices when running on GPU or with OpenMP (i.e.
-      // when Device::IsEnabled() returns true). This makes the results
-      // bit-for-bit deterministic at the cost of somewhat longer run time.
-      a.EnableSparseMatrixSorting(Device::IsEnabled());
-   }
+   if (fa) { a.SetAssemblyLevel(AssemblyLevel::FULL); }
   a.AddDomainIntegrator(new DiffusionIntegrator(one));

   // 12. Assemble the parallel bilinear form and the corresponding linear
@@ -14,5 +14,4 @@
 -mg_levels_esteig_ksp_type cg
 -mg_levels_esteig_ksp_max_it 10
 -mg_levels_ksp_chebyshev_esteig 0,0.05,0,1.05
-pc_gamg_use_sa_esteig 0
 -mg_levels_pc_type sor
@@ -124,7 +124,6 @@ void BilinearForm::SetAssemblyLevel(AssemblyLevel assembly_level)
      case AssemblyLevel::LEGACY:
         break;
      case AssemblyLevel::FULL:
-         SetDiagonalPolicy( DIAG_ONE ); // Only diagonal policy supported on device
         ext = new FABilinearFormExtension(this);
         break;
      case AssemblyLevel::ELEMENT:
@@ -80,9 +80,6 @@ protected:
   /** @brief Extension for supporting Full Assembly (FA), Element Assembly (EA),
       Partial Assembly (PA), or Matrix Free assembly (MF). */
   BilinearFormExtension *ext;
-   /** Indicates if the sparse matrix is sorted after assembly when using
-       Full Assembly (FA). */
-   bool sort_sparse_matrix = false;

   /** @brief Indicates the Mesh::sequence corresponding to the current state of
       the BilinearForm. */
@@ -184,21 +181,6 @@ public:
       If used, this method must be called before assembly. */
   void SetAssemblyLevel(AssemblyLevel assembly_level);

-   /** @brief Force the sparse matrix column indices to be sorted when using
-       AssemblyLevel::FULL.
-
-       When assembling on device the assembly algorithm uses atomic operations
-       to insert values in the sparse matrix, which can result in different
-       column index orderings across runs. Calling this method with @a enable_it
-       set to @a true forces a sorting algorithm to be called at the end of the
-       assembly procedure to ensure sorted column indices (and therefore
-       deterministic results).
-   */
-   void EnableSparseMatrixSorting(bool enable_it)
-   {
-      sort_sparse_matrix = enable_it;
-   }
-
   /// Returns the assembly level
   AssemblyLevel GetAssemblyLevel() const { return assembly; }

@@ -18,6 +18,8 @@
 #include "pgridfunc.hpp"
 #include "ceed/interface/util.hpp"

+#include "../general/nvtx.hpp"
+
 namespace mfem
 {

@@ -289,6 +291,10 @@ void PABilinearFormExtension::SetupRestrictionOperators(const L2FaceValues m)

 void PABilinearFormExtension::Assemble()
 {
+#undef MFEM_NVTX_COLOR
+#define MFEM_NVTX_COLOR NavyBlue
+   NVTX("HO Assemble");
+
   SetupRestrictionOperators(L2FaceValues::DoubleValued);

   Array<BilinearFormIntegrator*> &integrators = *a->GetDBFI();
@@ -383,6 +389,10 @@ void PABilinearFormExtension::FormLinearSystem(const Array<int> &ess_tdof_list,

 void PABilinearFormExtension::Mult(const Vector &x, Vector &y) const
 {
+#undef MFEM_NVTX_COLOR
+#define MFEM_NVTX_COLOR MediumSpringGreen
+   NVTX("HO Apply");
+
   Array<BilinearFormIntegrator*> &integrators = *a->GetDBFI();

   const int iSz = integrators.Size();
@@ -955,10 +965,6 @@ void FABilinearFormExtension::Assemble()
      }
      a->mat = mat;
   }
-   if ( a->sort_sparse_matrix )
-   {
-      a->mat->SortColumnIndices();
-   }
 }


@@ -2338,7 +2338,7 @@ void MixedCurlIntegrator::AssembleElementMatrix2(
   if (spaceH1)
   {
      dshape.SetSize(trial_dof,dim);
-      curlshape.SetSize(trial_dof,dim);
+      curlshape.SetSize(dim*trial_dof,1);
      dimc = dim;
   }
   else
@@ -2367,7 +2367,7 @@ void MixedCurlIntegrator::AssembleElementMatrix2(
      if (spaceH1)
      {
         trial_fe.CalcPhysDShape(Trans, dshape);
-         dshape.GradToVectorCurl2D(curlshape);
+         dshape.GradToCurl(curlshape);
      }
      else
      {
@@ -3059,8 +3059,8 @@ public:

 /** Integrator for the DG form:

-        - < {(Q grad(u)).n}, [v] > + sigma < [u], {(Q grad(v)).n} >
-        + kappa < {h^{-1} Q} [u], [v] >
+    - < {(Q grad(u)).n}, [v] > + sigma < [u], {(Q grad(v)).n} >
+    + kappa < {h^{-1} Q} [u], [v] >,

    where Q is a scalar or matrix diffusion coefficient and u, v are the trial
    and test spaces, respectively. The parameters sigma and kappa determine the
@@ -1686,7 +1686,7 @@ static void PADiffusionApply(const int dim,
         case 0x77: return SmemPADiffusionApply2D<7,7,4>(NE,symm,B,G,D,X,Y);
         case 0x88: return SmemPADiffusionApply2D<8,8,2>(NE,symm,B,G,D,X,Y);
         case 0x99: return SmemPADiffusionApply2D<9,9,2>(NE,symm,B,G,D,X,Y);
-         default:   return PADiffusionApply2D(NE,symm,B,G,Bt,Gt,D,X,Y,D1D,Q1D);
+            // default:   return PADiffusionApply2D(NE,symm,B,G,Bt,Gt,D,X,Y,D1D,Q1D);
      }
   }

@@ -1704,7 +1704,14 @@ static void PADiffusionApply(const int dim,
         case 0x67: return SmemPADiffusionApply3D<6,7>(NE,symm,B,G,D,X,Y);
         case 0x78: return SmemPADiffusionApply3D<7,8>(NE,symm,B,G,D,X,Y);
         case 0x89: return SmemPADiffusionApply3D<8,9>(NE,symm,B,G,D,X,Y);
-         default:   return PADiffusionApply3D(NE,symm,B,G,Bt,Gt,D,X,Y,D1D,Q1D);
+         case 0x33: return SmemPADiffusionApply3D<3,3>(NE,symm,B,G,D,X,Y);
+         case 0x44: return SmemPADiffusionApply3D<4,4>(NE,symm,B,G,D,X,Y);
+         case 0x55: return SmemPADiffusionApply3D<5,5>(NE,symm,B,G,D,X,Y);
+         case 0x66: return SmemPADiffusionApply3D<6,6>(NE,symm,B,G,D,X,Y);
+         case 0x77: return SmemPADiffusionApply3D<7,7>(NE,symm,B,G,D,X,Y);
+         case 0x88: return SmemPADiffusionApply3D<8,8>(NE,symm,B,G,D,X,Y);
+         case 0x99: return SmemPADiffusionApply3D<9,9>(NE,symm,B,G,D,X,Y);
+            // default:   return PADiffusionApply3D(NE,symm,B,G,Bt,Gt,D,X,Y,D1D,Q1D);
      }
   }
   MFEM_ABORT("Unknown kernel: 0x"<<std::hex << id << std::dec);
@@ -47,7 +47,6 @@ void InitRestrictionWithIndices(const FiniteElementSpace &fes,

 /** @brief Initialize a strided CeedElemRestriction

-    @param[in] fes Input finite element space.
    @param[in] nelem is the number of elements.
    @param[in] nqpts is the total number of quadrature points.
    @param[in] qdatasize is the number of data per quadrature point.
@@ -140,11 +140,7 @@ int CeedOperatorGetActiveField(CeedOperator oper, CeedOperatorField *field)
   CeedOperator *subops;
   if (isComposite)
   {
-#if CEED_VERSION_GE(0, 10, 2)
-      ierr = CeedCompositeOperatorGetSubList(oper, &subops); CeedChk(ierr);
-#else
      ierr = CeedOperatorGetSubList(oper, &subops); CeedChk(ierr);
-#endif
      ierr = CeedOperatorGetQFunction(subops[0], &qf); CeedChk(ierr);
   }
   else
@@ -275,14 +275,9 @@ CeedOperator CoarsenCeedCompositeOperator(
                                      &op_coarse); PCeedChk(ierr);

   int nsub;
-   CeedOperator *subops;
-#if CEED_VERSION_GE(0, 10, 2)
-   ierr = CeedCompositeOperatorGetNumSub(op, &nsub); PCeedChk(ierr);
-   ierr = CeedCompositeOperatorGetSubList(op, &subops); PCeedChk(ierr);
-#else
   ierr = CeedOperatorGetNumSub(op, &nsub); PCeedChk(ierr);
+   CeedOperator *subops;
   ierr = CeedOperatorGetSubList(op, &subops); PCeedChk(ierr);
-#endif
   for (int isub=0; isub<nsub; ++isub)
   {
      CeedOperator subop = subops[isub];
@@ -310,13 +310,8 @@ int CeedOperatorFullAssemble(CeedOperator op, SparseMatrix **mat)
   {
      CeedInt numsub;
      CeedOperator *subops;
-#if CEED_VERSION_GE(0, 10, 2)
-      CeedCompositeOperatorGetNumSub(op, &numsub);
-      ierr = CeedCompositeOperatorGetSubList(op, &subops); CeedChk(ierr);
-#else
      CeedOperatorGetNumSub(op, &numsub);
      ierr = CeedOperatorGetSubList(op, &subops); CeedChk(ierr);
-#endif
      for (int i = 0; i < numsub; ++i)
      {
         ierr = CeedSingleOperatorFullAssemble(subops[i], out); CeedChk(ierr);
@@ -66,7 +66,6 @@ int CeedBasisATPMGCoarsen(CeedBasis basisin, CeedBasis* basisout,
    @param[in] coarse_er         CeedElemRestriction for coarse operator
                                 (see CeedATPMGElemRestriction)
    @param[out] coarse_basis_out CeedBasis for coarser operator
-    @param[out] basis_ctof_out   CeedBasis describing interpolation from coarse to fine
    @param[out] out              coarsened CeedOperator
 */
 int CeedATPMGOperator(CeedOperator oper, int order_reduction,
@@ -724,7 +724,7 @@ void VisItDataCollection::ParseVisItRootString(const std::string& json)

   // Set the DataCollection::name using the mesh path
   std::string path = mesh.get("path").get<std::string>();
-   size_t right_sep = path.rfind('_');
+   size_t right_sep = path.find('_');
   if (right_sep == std::string::npos)
   {
      error = READ_ERROR;
@@ -338,12 +338,11 @@ public:
   /// Set the precision (number of digits) used for the text output of doubles
   void SetPrecision(int prec) { precision = prec; }
   /// Set the number of digits used for both the cycle and the MPI rank
-   virtual void SetPadDigits(int digits)
-   { pad_digits_cycle=pad_digits_rank = digits; }
+   void SetPadDigits(int digits) { pad_digits_cycle=pad_digits_rank = digits; }
   /// Set the number of digits used for the cycle
-   virtual void SetPadDigitsCycle(int digits) { pad_digits_cycle = digits; }
+   void SetPadDigitsCycle(int digits) { pad_digits_cycle = digits; }
   /// Set the number of digits used for the MPI rank in filenames
-   virtual void SetPadDigitsRank(int digits) { pad_digits_rank = digits; }
+   void SetPadDigitsRank(int digits) { pad_digits_rank = digits; }
   /// Set the desired output mesh and data format.
   /** See the enumeration #Format for valid options. Derived classes can define
       their own format enumerations and override this method to perform input
@@ -442,29 +441,21 @@ public:
 #endif

   /// Set/change the mesh associated with the collection
-   virtual void SetMesh(Mesh *new_mesh) override;
+   virtual void SetMesh(Mesh *new_mesh);

 #ifdef MFEM_USE_MPI
   /// Set/change the mesh associated with the collection.
-   virtual void SetMesh(MPI_Comm comm, Mesh *new_mesh) override;
+   virtual void SetMesh(MPI_Comm comm, Mesh *new_mesh);
 #endif

   /// Add a grid function to the collection and update the root file
-   virtual void RegisterField(const std::string& field_name,
-                              GridFunction *gf) override;
+   virtual void RegisterField(const std::string& field_name, GridFunction *gf);

   /// Add a quadrature function to the collection and update the root file.
   /** Visualization of quadrature function is not supported in VisIt(3.12).
       A patch has been sent to VisIt developers in June 2020. */
   virtual void RegisterQField(const std::string& q_field_name,
-                               QuadratureFunction *qf) override;
-
-   /// Set the number of digits used for both the cycle and the MPI rank
-   /// @note VisIt seems to require 6 pad digits for the MPI rank. Therefore,
-   /// this function uses this default value. This behavior can be overridden
-   /// by calling SetPadDigitsCycle() and SetPadDigitsRank() instead.
-   virtual void SetPadDigits(int digits) override
-   { pad_digits_cycle=digits; pad_digits_rank=6; }
+                               QuadratureFunction *qf);

   /// Set VisIt parameter: default levels of detail for the MultiresControl
   void SetLevelsOfDetail(int levels_of_detail);
@@ -477,13 +468,13 @@ public:
   void DeleteAll();

   /// Save the collection and a VisIt root file
-   virtual void Save() override;
+   virtual void Save();

   /// Save a VisIt root file for the collection
   void SaveRootFile();

   /// Load the collection based on its VisIt data (described in its root file)
-   virtual void Load(int cycle_ = 0) override;
+   virtual void Load(int cycle_ = 0);

   /// We will delete the mesh and fields if we own them
   virtual ~VisItDataCollection() {}
@@ -339,11 +339,8 @@ ND_TriDofTransformation::TransformDual(double *v) const
 void
 ND_TriDofTransformation::InvTransformDual(double *v) const
 {
-   // Return immediately when no face DoFs are present
-   if (nfdofs < 2) { return; }
-
-   MFEM_VERIFY(Fo.Size() >= 1,
-               "Face orientations are unset in ND_TriDofTransformation");
+   int nedofs = order; // number of DoFs per edge
+   int nfdofs = order*(order-1); // number of DoFs per face

   double data[2];
   Vector v2(data, 2);
@@ -435,11 +432,8 @@ ND_TetDofTransformation::TransformDual(double *v) const
 void
 ND_TetDofTransformation::InvTransformDual(double *v) const
 {
-   // Return immediately when no face DoFs are present
-   if (nfdofs < 2) { return; }
-
-   MFEM_VERIFY(Fo.Size() >= 4,
-               "Face orientations are unset in ND_TetDofTransformation");
+   int nedofs = order; // number of DoFs per edge
+   int nfdofs = order*(order-1); // number of DoFs per face

   double data[2];
   Vector v2(data, 2);
@@ -618,7 +618,7 @@ void InvertLinearTrans(ElementTransformation &trans,

   double store[3];
   Vector v(store, x.Size());
-   pt.Get(store, x.Size());
+   pt.Get(v, x.Size());
   v -= x;

   trans.InverseJacobian().Mult(v, x);
@@ -954,8 +954,8 @@ void VectorFiniteElement::Project_RT(
   {
      Trans.SetIntPoint(&Nodes.IntPoint(k));
      // dof_k = nk^t adj(J) xk
-      dofs(k) = Trans.AdjugateJacobian().InnerProduct(
-                   &vc[k*sdim], nk + d2n[k]*dim);
+      Vector vk(vc.GetData()+k*sdim, sdim);
+      dofs(k) = Trans.AdjugateJacobian().InnerProduct(vk, nk + d2n[k]*dim);
      if (!square_J) { dofs(k) /= Trans.Weight(); }
   }
 }
@@ -1171,8 +1171,9 @@ void VectorFiniteElement::Project_ND(
   for (int k = 0; k < dof; k++)
   {
      Trans.SetIntPoint(&Nodes.IntPoint(k));
+      Vector vk(vc.GetData()+k*dim, dim);
      // dof_k = xk^t J tk
-      dofs(k) = Trans.Jacobian().InnerProduct(tk + d2t[k]*dim, &vc[k*dim]);
+      dofs(k) = Trans.Jacobian().InnerProduct(tk + d2t[k]*dim, vk);
   }
 }

@@ -1319,7 +1320,7 @@ void VectorFiniteElement::LocalL2Projection_RT(
      double w = ip.weight;
      this->CalcVShape(ip, fine_shape);
      Trans.Transform(ip, v);
-      tr_ip.Set(v.GetData(), dim);
+      tr_ip.Set(v, dim);
      cfe.CalcVShape(tr_ip, coarse_shape);

      AddMult_a_AAt(w, fine_shape, fine_mass);
@@ -1406,7 +1407,7 @@ void VectorFiniteElement::LocalL2Projection_ND(
      const IntegrationPoint &ip = ir.IntPoint(i);
      this->CalcVShape(ip, fine_shape);
      Trans.Transform(ip, v);
-      tr_ip.Set(v.GetData(), dim);
+      tr_ip.Set(v, dim);
      cfe.CalcVShape(tr_ip, coarse_shape);

      AddMult_a_AAt(ip.weight, fine_shape, fine_mass);
@@ -256,7 +256,7 @@ protected:

 public:
   /// Enumeration for range_type and deriv_range_type
-   enum RangeType { UNKNOWN_RANGE_TYPE = -1, SCALAR, VECTOR };
+   enum RangeType { SCALAR, VECTOR };

   /** @brief Enumeration for MapType: defines how reference functions are
       mapped to physical space.
@@ -270,8 +270,6 @@ public:
   */
   enum MapType
   {
-      UNKNOWN_MAP_TYPE = -1, /**< Used to distinguish an unset MapType variable
-                                  from the known values below. */
      VALUE,     /**< For scalar fields; preserves point values
                          \f$ u(x) = \hat u(\hat x) \f$ */
      INTEGRAL,  /**< For scalar fields; preserves volume integrals
@@ -1093,11 +1091,6 @@ public:
   // { CalcLegendre(p, x, u); }
   { CalcChebyshev(p, x, u); }

-   /** @brief Evaluate the values of a hierarchical 1D basis at point x
-       hierarchical = k-th basis function is degree k polynomial */
-   static void CalcBasis(const int p, const double x, Vector &u)
-   { CalcBasis(p, x, u.GetData()); }
-
   /// Evaluate the values and derivatives of a hierarchical 1D basis at point @a x
   static void CalcBasis(const int p, const double x, double *u, double *d)
   // { CalcMono(p, x, u, d); }
@@ -1105,11 +1098,6 @@ public:
   // { CalcLegendre(p, x, u, d); }
   { CalcChebyshev(p, x, u, d); }

-   /** @brief Evaluate the values and derivatives of a hierarchical 1D basis at
-       point @a x. */
-   static void CalcBasis(const int p, const double x, Vector &u, Vector &d)
-   { CalcBasis(p, x, u.GetData(), d.GetData()); }
-
   /// Evaluate the values, derivatives and second derivatives of a hierarchical 1D basis at point x
   static void CalcBasis(const int p, const double x, double *u, double *d,
                         double *dd)
@@ -1118,12 +1106,6 @@ public:
   // { CalcLegendre(p, x, u, d); }
   { CalcChebyshev(p, x, u, d, dd); }

-   /** @brief Evaluate the values, derivatives and second derivatives of a
-       hierarchical 1D basis at point @a x. */
-   static void CalcBasis(const int p, const double x, Vector &u, Vector &d,
-                         Vector &dd)
-   { CalcBasis(p, x, u.GetData(), d.GetData(), dd.GetData()); }
-
   /// Evaluate a representation of a Delta function at point x
   static double CalcDelta(const int p, const double x)
   { return pow(x, (double) p); }
@@ -1153,24 +1135,12 @@ public:
   static void CalcBernstein(const int p, const double x, double *u)
   { CalcBinomTerms(p, x, 1. - x, u); }

-   /** @brief Compute the values of the Bernstein basis functions of order
-       @a p at coordinate @a x and store the results in the already allocated
-       @a u array. */
-   static void CalcBernstein(const int p, const double x, Vector &u)
-   { CalcBernstein(p, x, u.GetData()); }
-
   /** @brief Compute the values and derivatives of the Bernstein basis functions
       of order @a p at coordinate @a x and store the results in the already allocated
       @a u and @a d arrays. */
   static void CalcBernstein(const int p, const double x, double *u, double *d)
   { CalcBinomTerms(p, x, 1. - x, u, d); }

-   /** @brief Compute the values and derivatives of the Bernstein basis
-       functions of order @a p at coordinate @a x and store the results in the
-       already allocated @a u and @a d arrays. */
-   static void CalcBernstein(const int p, const double x, Vector &u, Vector &d)
-   { CalcBernstein(p, x, u.GetData(), d.GetData()); }
-
   static void CalcLegendre(const int p, const double x, double *u);
   static void CalcLegendre(const int p, const double x, double *u, double *d);

@@ -319,9 +319,6 @@ void ND_HexahedronElement::CalcVShape(const IntegrationPoint &ip,

   if (obasis1d.IsIntegratedType())
   {
-#ifdef MFEM_THREAD_SAFE
-      Vector dshape_cx(p + 1), dshape_cy(p + 1), dshape_cz(p + 1);
-#endif
      cbasis1d.Eval(ip.x, shape_cx, dshape_cx);
      cbasis1d.Eval(ip.y, shape_cy, dshape_cy);
      cbasis1d.Eval(ip.z, shape_cz, dshape_cz);
@@ -664,9 +661,6 @@ void ND_QuadrilateralElement::CalcVShape(const IntegrationPoint &ip,

   if (obasis1d.IsIntegratedType())
   {
-#ifdef MFEM_THREAD_SAFE
-      Vector dshape_cx(p + 1), dshape_cy(p + 1);
-#endif
      cbasis1d.Eval(ip.x, shape_cx, dshape_cx);
      cbasis1d.Eval(ip.y, shape_cy, dshape_cy);
      obasis1d.ScaleIntegrated(false);
@@ -212,7 +212,7 @@ void BiQuadPos2DFiniteElement::GetLocalInterpolation(
 void BiQuadPos2DFiniteElement::Project(
   Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
 {
-   double *d = dofs.GetData();
+   double *d = dofs;

   for (int i = 0; i < 9; i++)
   {
@@ -382,8 +382,8 @@ void H1Pos_QuadrilateralElement::CalcShape(const IntegrationPoint &ip,
   Vector shape_x(p+1), shape_y(p+1);
 #endif

-   Poly_1D::CalcBernstein(p, ip.x, shape_x);
-   Poly_1D::CalcBernstein(p, ip.y, shape_y);
+   Poly_1D::CalcBernstein(p, ip.x, shape_x.GetData() );
+   Poly_1D::CalcBernstein(p, ip.y, shape_y.GetData() );

   // Reorder so that vertices are at the beginning of the list
   for (int o = 0, j = 0; j <= p; j++)
@@ -402,8 +402,8 @@ void H1Pos_QuadrilateralElement::CalcDShape(const IntegrationPoint &ip,
   Vector shape_x(p+1), shape_y(p+1), dshape_x(p+1), dshape_y(p+1);
 #endif

-   Poly_1D::CalcBernstein(p, ip.x, shape_x, dshape_x);
-   Poly_1D::CalcBernstein(p, ip.y, shape_y, dshape_y);
+   Poly_1D::CalcBernstein(p, ip.x, shape_x.GetData(), dshape_x.GetData() );
+   Poly_1D::CalcBernstein(p, ip.y, shape_y.GetData(), dshape_y.GetData() );

   // Reorder so that vertices are at the beginning of the list
   for (int o = 0, j = 0; j <= p; j++)
@@ -150,9 +150,6 @@ void RT_QuadrilateralElement::CalcVShape(const IntegrationPoint &ip,

   if (obasis1d.IsIntegratedType())
   {
-#ifdef MFEM_THREAD_SAFE
-      Vector dshape_cx(pp1 + 1), dshape_cy(pp1 + 1);
-#endif
      cbasis1d.Eval(ip.x, shape_cx, dshape_cx);
      cbasis1d.Eval(ip.y, shape_cy, dshape_cy);
      obasis1d.ScaleIntegrated(false);
@@ -481,9 +478,6 @@ void RT_HexahedronElement::CalcVShape(const IntegrationPoint &ip,

   if (obasis1d.IsIntegratedType())
   {
-#ifdef MFEM_THREAD_SAFE
-      Vector dshape_cx(pp1 + 1), dshape_cy(pp1 + 1), dshape_cz(pp1 + 1);
-#endif
      cbasis1d.Eval(ip.x, shape_cx, dshape_cx);
      cbasis1d.Eval(ip.y, shape_cy, dshape_cy);
      cbasis1d.Eval(ip.z, shape_cz, dshape_cz);
@@ -22,70 +22,6 @@ namespace mfem

 using namespace std;

-const FiniteElement *
-FiniteElementCollection::FiniteElementForDim(int dim) const
-{
-   for (int g = Geometry::DimStart[dim]; g < Geometry::DimStart[dim+1]; g++)
-   {
-      const FiniteElement *fe = FiniteElementForGeometry((Geometry::Type)g);
-      if (fe != NULL)
-      {
-         return fe;
-      }
-   }
-   return NULL;
-}
-
-int FiniteElementCollection::GetRangeType(int dim) const
-{
-   const FiniteElement *fe = FiniteElementForDim(dim);
-   if (fe)
-   {
-      return fe->GetRangeType();
-   }
-   return FiniteElement::UNKNOWN_RANGE_TYPE;
-}
-
-int FiniteElementCollection::GetDerivRangeType(int dim) const
-{
-   const FiniteElement *fe = FiniteElementForDim(dim);
-   if (fe)
-   {
-      return fe->GetDerivRangeType();
-   }
-   return FiniteElement::UNKNOWN_RANGE_TYPE;
-}
-
-int FiniteElementCollection::GetMapType(int dim) const
-{
-   const FiniteElement *fe = FiniteElementForDim(dim);
-   if (fe)
-   {
-      return fe->GetMapType();
-   }
-   return FiniteElement::UNKNOWN_MAP_TYPE;
-}
-
-int FiniteElementCollection::GetDerivType(int dim) const
-{
-   const FiniteElement *fe = FiniteElementForDim(dim);
-   if (fe)
-   {
-      return fe->GetDerivType();
-   }
-   return FiniteElement::NONE;
-}
-
-int FiniteElementCollection::GetDerivMapType(int dim) const
-{
-   const FiniteElement *fe = FiniteElementForDim(dim);
-   if (fe)
-   {
-      return fe->GetDerivMapType();
-   }
-   return FiniteElement::UNKNOWN_MAP_TYPE;
-}
-
 int FiniteElementCollection::HasFaceDofs(Geometry::Type geom, int p) const
 {
   switch (geom)
@@ -51,14 +51,6 @@ public:
   virtual const FiniteElement *
   FiniteElementForGeometry(Geometry::Type GeomType) const = 0;

-   /** @brief Returns the first non-NULL FiniteElement for the given dimension
-
-       @note Repeatedly calls FiniteElementForGeometry in the order defined in
-       the Geometry::Type enumeration.
-   */
-   virtual const FiniteElement *
-   FiniteElementForDim(int dim) const;
-
   virtual int DofForGeometry(Geometry::Type GeomType) const = 0;

   /** @brief Returns an array, say p, that maps a local permuted index i to a
@@ -74,17 +66,6 @@ public:

   virtual int GetContType() const = 0;

-   /** @note The following methods provide the same information as the
-       corresponding methods of the FiniteElement base class.
-       @{
-   */
-   virtual int GetRangeType(int dim) const;
-   virtual int GetDerivRangeType(int dim) const;
-   virtual int GetMapType(int dim) const;
-   virtual int GetDerivType(int dim) const;
-   virtual int GetDerivMapType(int dim) const;
-   /** @} */
-
   int HasFaceDofs(Geometry::Type geom, int p) const;

   virtual const FiniteElement *TraceFiniteElementForGeometry(
@@ -1920,7 +1920,7 @@ DataCollectionToFmsDataCollection(DataCollection *mfem_dc,
                                             FMS_NODAL_GAUSS_CLOSED, 1);
      err |= FmsDataCollectionAddField(*dc, "Coords", &fcoords);
      err |= FmsFieldSet(fcoords, fdcoords, mmesh->SpaceDimension(), FMS_BY_NODES,
-                         FMS_DOUBLE, mverts.HostRead());
+                         FMS_DOUBLE, mverts);
      err |= FmsComponentSetCoordinates(volume, fcoords);
   }

@@ -416,24 +416,13 @@ void GridFunction::GetNodalValues(int i, Array<double> &nval, int vdim) const

   if (FElem->GetRangeType() == FiniteElement::SCALAR)
   {
+      MFEM_ASSERT(FElem->GetMapType() == FiniteElement::VALUE,
+                  "invalid FE map type");
      Vector shape(dof);
-      if (FElem->GetMapType() == FiniteElement::VALUE)
+      for (k = 0; k < n; k++)
      {
-         for (k = 0; k < n; k++)
-         {
-            FElem->CalcShape(ElemVert->IntPoint(k), shape);
-            nval[k] = shape * (&loc_data[dof * vdim]);
-         }
-      }
-      else
-      {
-         ElementTransformation *Tr = fes->GetElementTransformation(i);
-         for (k = 0; k < n; k++)
-         {
-            Tr->SetIntPoint(&ElemVert->IntPoint(k));
-            FElem->CalcPhysShape(*Tr, shape);
-            nval[k] = shape * (&loc_data[dof * vdim]);
-         }
+         FElem->CalcShape(ElemVert->IntPoint(k), shape);
+         nval[k] = shape * ((const double *)loc_data + dof * vdim);
      }
   }
   else
@@ -506,7 +495,7 @@ void GridFunction::GetVectorValue(int i, const IntegrationPoint &ip,
      val.SetSize(vdim);
      for (int k = 0; k < vdim; k++)
      {
-         val(k) = shape * (&loc_data[dof * k]);
+         val(k) = shape * ((const double *)loc_data + dof * k);
      }
   }
   else
@@ -1038,7 +1027,7 @@ void GridFunction::GetVectorValue(ElementTransformation &T,
      val.SetSize(vdim);
      for (int k = 0; k < vdim; k++)
      {
-         val(k) = shape * (&loc_data[dof * k]);
+         val(k) = shape * ((const double *)loc_data + dof * k);
      }
   }
   else
@@ -1089,7 +1078,7 @@ void GridFunction::GetVectorValues(ElementTransformation &T,

         for (int k = 0; k < vdim; k++)
         {
-            vals(k,j) = shape * (&loc_data[dof * k]);
+            vals(k,j) = shape * ((const double *)loc_data + dof * k);
         }
      }
   }
@@ -1192,7 +1181,8 @@ void GridFunction::GetValuesFrom(const GridFunction &orig_func)
         orig_fe->CalcShape(ip, shape);
         for (d = 0; d < vdim; d++)
         {
-            loc_values(d*dof+j) = shape * (&orig_loc_values[d * odof]);
+            loc_values(d*dof+j) =
+               shape * ((const double *)orig_loc_values + d * odof) ;
         }
      }
      if (doftrans)
@@ -1234,7 +1224,8 @@ void GridFunction::GetBdrValuesFrom(const GridFunction &orig_func)
         orig_fe->CalcShape(ip, shape);
         for (d = 0; d < vdim; d++)
         {
-            loc_values(d*dof+j) = shape * (&orig_loc_values[d * odof]);
+            loc_values(d*dof+j) =
+               shape * ((const double *)orig_loc_values + d * odof);
         }
      }
      SetSubVector(vdofs, loc_values);
@@ -4108,16 +4099,16 @@ void TensorProductLegendre(int dim,                // input
   // Map x to [0, 1] to use CalcLegendre since it uses shifted Legendre Polynomials.
   double x1 = (x(0) - xmin(0))/(xmax(0)-xmin(0)), x2, x3;
   Vector poly_x(order+1), poly_y(order+1), poly_z(order+1);
-   poly1d.CalcLegendre(order, x1, poly_x.GetData());
+   poly1d.CalcLegendre(order, x1, poly_x);
   if (dim > 1)
   {
      x2 = (x(1)-xmin(1))/(xmax(1)-xmin(1));
-      poly1d.CalcLegendre(order, x2, poly_y.GetData());
+      poly1d.CalcLegendre(order, x2, poly_y);
   }
   if (dim == 3)
   {
      x3 = (x(2)-xmin(2))/(xmax(2)-xmin(2));
-      poly1d.CalcLegendre(order, x3, poly_z.GetData());
+      poly1d.CalcLegendre(order, x3, poly_z);
   }

   int basis_dimension = static_cast<int>(pow(order+1,dim));
@@ -179,27 +179,24 @@ void FindPointsGSLIB::FindPoints(const Vector &point_pos,
   gsl_ref.SetSize(points_cnt * dim);
   gsl_dist.SetSize(points_cnt);

-   auto xvFill = [&](const double *xv_base[], unsigned xv_stride[], int dim)
+   const double *xv_base[dim];
+   unsigned xv_stride[dim];
+   for (int d = 0; d < dim; d++)
   {
-      for (int d = 0; d < dim; d++)
+      if (point_pos_ordering == Ordering::byNODES)
      {
-         if (point_pos_ordering == Ordering::byNODES)
-         {
-            xv_base[d] = point_pos.GetData() + d*points_cnt;
-            xv_stride[d] = sizeof(double);
-         }
-         else
-         {
-            xv_base[d] = point_pos.GetData() + d;
-            xv_stride[d] = dim*sizeof(double);
-         }
+         xv_base[d] = point_pos.GetData() + d*points_cnt;
+         xv_stride[d] = sizeof(double);
      }
-   };
+      else
+      {
+         xv_base[d] = point_pos.GetData() + d;
+         xv_stride[d] = dim*sizeof(double);
+      }
+   }
+
   if (dim == 2)
   {
-      const double *xv_base[2];
-      unsigned xv_stride[2];
-      xvFill(xv_base, xv_stride, dim);
      findpts_2(gsl_code.GetData(), sizeof(unsigned int),
                gsl_proc.GetData(), sizeof(unsigned int),
                gsl_elem.GetData(), sizeof(unsigned int),
@@ -207,11 +204,8 @@ void FindPointsGSLIB::FindPoints(const Vector &point_pos,
                gsl_dist.GetData(), sizeof(double),
                xv_base, xv_stride, points_cnt, fdata2D);
   }
-   else  // dim == 3
+   else
   {
-      const double *xv_base[3];
-      unsigned xv_stride[3];
-      xvFill(xv_base, xv_stride, dim);
      findpts_3(gsl_code.GetData(), sizeof(unsigned int),
                gsl_proc.GetData(), sizeof(unsigned int),
                gsl_elem.GetData(), sizeof(unsigned int),
@@ -1086,7 +1080,7 @@ void FindPointsGSLIB::InterpolateGeneral(const GridFunction &field_in,

            sarray_transfer(struct send_pt, sendpt, proc, 1, cr);
            sdpt = (struct send_pt *)sendpt->ptr;
-            for (int index = 0; index < static_cast<int>(sendpt->n); index++)
+            for (int index = 0; index < sendpt->n; index++)
            {
               int idx = field_in.FESpace()->GetOrdering() == Ordering::byNODES ?
                         sdpt->index + j*nptorig :
@@ -1200,27 +1194,24 @@ void OversetFindPointsGSLIB::FindPoints(const Vector &point_pos,
   gsl_ref.SetSize(points_cnt * dim);
   gsl_dist.SetSize(points_cnt);

-   auto xvFill = [&](const double *xv_base[], unsigned xv_stride[], int dim)
+   const double *xv_base[dim];
+   unsigned xv_stride[dim];
+   for (int d = 0; d < dim; d++)
   {
-      for (int d = 0; d < dim; d++)
+      if (point_pos_ordering == Ordering::byNODES)
      {
-         if (point_pos_ordering == Ordering::byNODES)
-         {
-            xv_base[d] = point_pos.GetData() + d*points_cnt;
-            xv_stride[d] = sizeof(double);
-         }
-         else
-         {
-            xv_base[d] = point_pos.GetData() + d;
-            xv_stride[d] = dim*sizeof(double);
-         }
+         xv_base[d] = point_pos.GetData() + d*points_cnt;
+         xv_stride[d] = sizeof(double);
      }
-   };
+      else
+      {
+         xv_base[d] = point_pos.GetData() + d;
+         xv_stride[d] = dim*sizeof(double);
+      }
+   }
+
   if (dim == 2)
   {
-      const double *xv_base[2];
-      unsigned xv_stride[2];
-      xvFill(xv_base, xv_stride, dim);
      findptsms_2(gsl_code.GetData(), sizeof(unsigned int),
                  gsl_proc.GetData(), sizeof(unsigned int),
                  gsl_elem.GetData(), sizeof(unsigned int),
@@ -1230,11 +1221,8 @@ void OversetFindPointsGSLIB::FindPoints(const Vector &point_pos,
                  point_id.GetData(), sizeof(unsigned int), &match,
                  points_cnt, fdata2D);
   }
-   else  // dim == 3
+   else
   {
-      const double *xv_base[3];
-      unsigned xv_stride[3];
-      xvFill(xv_base, xv_stride, dim);
      findptsms_3(gsl_code.GetData(), sizeof(unsigned int),
                  gsl_proc.GetData(), sizeof(unsigned int),
                  gsl_elem.GetData(), sizeof(unsigned int),
@@ -1032,12 +1032,12 @@ void WhiteGaussianNoiseDomainLFIntegrator::AssembleRHSElementVect
      massinteg.AssembleElementMatrix(el, Tr, *M);
      CholeskyFactors chol(M->Data());
      chol.Factor(M->Height());
-      chol.LMult(n,1,elvect.GetData());
+      chol.LMult(n,1,elvect);
   }
   else
   {
      CholeskyFactors chol(L[iel]->Data());
-      chol.LMult(n,1,elvect.GetData());
+      chol.LMult(n,1,elvect);
   }
 }

@@ -365,8 +365,11 @@ FiniteElementSpace &LORBase::GetFESpace() const

 void LORBase::AssembleSystem(BilinearForm &a_ho, const Array<int> &ess_dofs)
 {
-   A.Clear();
-   delete a;
+   if (a)
+   {
+      A.Clear();
+      delete a;
+   }
   if (BatchedLORAssembly::FormIsSupported(a_ho))
   {
      // Skip forming the space
@@ -13,6 +13,9 @@
 #include "../../general/forall.hpp"
 #include "../../fem/pbilinearform.hpp"

+#define MFEM_NVTX_COLOR DeepSkyBlue
+#include "../../general/nvtx.hpp"
+
 namespace mfem
 {

@@ -78,8 +81,10 @@ void BatchedLOR_ADS::Form3DFaceToEdge(Array<int> &face2edge)
   }
 }

-void BatchedLOR_ADS::FormCurlMatrix()
+void BatchedLOR_ADS::FormCurlMatrixLocal()
 {
+   NVTX("Discrete Curl");
+
   // The curl matrix maps from LOR edges to LOR faces. Given a quadrilateral
   // face (defined by its four edges) f_i = (e_j1, e_j2, e_j3, e_j4), the
   // matrix has nonzeros A(i, jk), so there are always exactly four nonzeros
@@ -87,10 +92,9 @@ void BatchedLOR_ADS::FormCurlMatrix()
   const int nface_dof = face_fes.GetNDofs();
   const int nedge_dof = edge_fes.GetNDofs();

-   SparseMatrix C_local;
   C_local.OverrideSize(nface_dof, nedge_dof);
-
-   C_local.GetMemoryI().New(nedge_dof+1, Device::GetDeviceMemoryType());
+   EnsureCapacity(C_local.GetMemoryI(), nedge_dof+1,
+                  Device::GetDeviceMemoryType());
   // Each row always has four nonzeros
   const int nnz = 4*nedge_dof;
   auto I = C_local.WriteI();
@@ -120,8 +124,8 @@ void BatchedLOR_ADS::FormCurlMatrix()
   const auto f2e = Reshape(face2edge.Read(), 4, nface_per_el);

   // Fill J and data
-   C_local.GetMemoryJ().New(nnz, Device::GetDeviceMemoryType());
-   C_local.GetMemoryData().New(nnz, Device::GetDeviceMemoryType());
+   EnsureCapacity(C_local.GetMemoryJ(), nnz, Device::GetDeviceMemoryType());
+   EnsureCapacity(C_local.GetMemoryData(), nnz, Device::GetDeviceMemoryType());

   auto J = C_local.WriteJ();
   auto V = C_local.WriteData();
@@ -146,6 +150,11 @@ void BatchedLOR_ADS::FormCurlMatrix()
         V[i*4 + k] = sgn*sgn_f*sgn_e;
      }
   });
+}
+
+void BatchedLOR_ADS::FormCurlMatrix()
+{
+   FormCurlMatrixLocal();

   // Create a block diagonal parallel matrix
   OperatorHandle C_diag(Operator::Hypre_ParCSR);
@@ -179,6 +188,8 @@ void BatchedLOR_ADS::FormCurlMatrix()
   }
   C->CopyRowStarts();
   C->CopyColStarts();
+
+   C_local.Clear();
 }

 HypreParMatrix *BatchedLOR_ADS::StealCurlMatrix()
@@ -36,7 +36,9 @@ protected:
   ND_FECollection edge_fec; ///< The associated Nedelec collection.
   ParFiniteElementSpace edge_fes; ///< The associated Nedelec space.
   BatchedLOR_AMS ams; ///< The associated AMS object.
-   HypreParMatrix *C; ///< The discrete curl matrix.
+   HypreParMatrix *C = nullptr; ///< The discrete curl matrix.
+
+   SparseMatrix C_local;

   /// Form the local elementwise discrete curl matrix.
   void Form3DFaceToEdge(Array<int> &face2edge);
@@ -64,6 +66,9 @@ public:

   /// Form the discrete curl matrix (not part of the public API).
   void FormCurlMatrix();
+
+   void FormCurlMatrixLocal();
+
   ~BatchedLOR_ADS();
 };

@@ -13,6 +13,9 @@
 #include "../../general/forall.hpp"
 #include "../../fem/pbilinearform.hpp"

+#define MFEM_NVTX_COLOR DeepSkyBlue
+#include "../../general/nvtx.hpp"
+
 namespace mfem
 {

@@ -136,8 +139,9 @@ void BatchedLOR_AMS::Form3DEdgeToVertex(Array<int> &edge2vert)
   }
 }

-void BatchedLOR_AMS::FormGradientMatrix()
+void BatchedLOR_AMS::FormGradientMatrixLocal()
 {
+   NVTX("Discrete Gradient");
   // The gradient matrix maps from LOR vertices to LOR edges. Given an edge
   // (defined by its two vertices) e_i = (v_j1, v_j2), the matrix has nonzeros
   // A(i, j1) = -1 and A(i, j2) = 1, so there are always exactly two nonzeros
@@ -145,10 +149,10 @@ void BatchedLOR_AMS::FormGradientMatrix()
   const int nedge_dof = edge_fes.GetNDofs();
   const int nvert_dof = vert_fes.GetNDofs();

-   SparseMatrix G_local;
   G_local.OverrideSize(nedge_dof, nvert_dof);
+   EnsureCapacity(G_local.GetMemoryI(), nedge_dof+1,
+                  Device::GetDeviceMemoryType());

-   G_local.GetMemoryI().New(nedge_dof+1, Device::GetDeviceMemoryType());
   // Each row always has two nonzeros
   const int nnz = 2*nedge_dof;
   auto I = G_local.WriteI();
@@ -180,8 +184,8 @@ void BatchedLOR_AMS::FormGradientMatrix()
   const auto e2v = Reshape(edge2vertex.Read(), 2, nedge_per_el);

   // Fill J and data
-   G_local.GetMemoryJ().New(nnz, Device::GetDeviceMemoryType());
-   G_local.GetMemoryData().New(nnz, Device::GetDeviceMemoryType());
+   EnsureCapacity(G_local.GetMemoryJ(), nnz, Device::GetDeviceMemoryType());
+   EnsureCapacity(G_local.GetMemoryData(), nnz, Device::GetDeviceMemoryType());

   auto J = G_local.WriteJ();
   auto V = G_local.WriteData();
@@ -204,6 +208,11 @@ void BatchedLOR_AMS::FormGradientMatrix()
      V[i*2 + 0] = -sgn;
      V[i*2 + 1] = sgn;
   });
+}
+
+void BatchedLOR_AMS::FormGradientMatrix()
+{
+   FormGradientMatrixLocal();

   // Create a block diagonal parallel matrix
   OperatorHandle G_diag(Operator::Hypre_ParCSR);
@@ -237,6 +246,8 @@ void BatchedLOR_AMS::FormGradientMatrix()
   }
   G->CopyRowStarts();
   G->CopyColStarts();
+
+   G_local.Clear();
 }

 template <typename T>
@@ -281,7 +292,7 @@ void BatchedLOR_AMS::FormCoordinateVectors(const Vector &X_vert)
   const MemoryClass mc = GetHypreMemoryClass();
   bool dev = (mc == MemoryClass::DEVICE);

-   xyz_tvec = new Vector(ntdofs*dim);
+   if (xyz_tvec == nullptr) { xyz_tvec = new Vector(ntdofs*dim); }

   auto xyz_tv = Reshape(HypreWrite(xyz_tvec->GetMemory()), ntdofs, dim);
   const auto xyz_e =
@@ -302,21 +313,24 @@ void BatchedLOR_AMS::FormCoordinateVectors(const Vector &X_vert)
   });

   // Make x, y, z HypreParVectors point to T-vector data
-   HYPRE_BigInt glob_size = vert_fes.GlobalTrueVSize();
-   HYPRE_BigInt *cols = vert_fes.GetTrueDofOffsets();
+   if (x == nullptr)
+   {
+      HYPRE_BigInt glob_size = vert_fes.GlobalTrueVSize();
+      HYPRE_BigInt *cols = vert_fes.GetTrueDofOffsets();

-   double *d_x_ptr = xyz_tv + 0*ntdofs;
-   x = new HypreParVector(vert_fes.GetComm(), glob_size, d_x_ptr, cols, dev);
-   double *d_y_ptr = xyz_tv + 1*ntdofs;
-   y = new HypreParVector(vert_fes.GetComm(), glob_size, d_y_ptr, cols, dev);
-   if (dim == 3)
-   {
-      double *d_z_ptr = xyz_tv + 2*ntdofs;
-      z = new HypreParVector(vert_fes.GetComm(), glob_size, d_z_ptr, cols, dev);
-   }
-   else
-   {
-      z = NULL;
+      double *d_x_ptr = xyz_tv + 0*ntdofs;
+      x = new HypreParVector(vert_fes.GetComm(), glob_size, d_x_ptr, cols, dev);
+      double *d_y_ptr = xyz_tv + 1*ntdofs;
+      y = new HypreParVector(vert_fes.GetComm(), glob_size, d_y_ptr, cols, dev);
+      if (dim == 3)
+      {
+         double *d_z_ptr = xyz_tv + 2*ntdofs;
+         z = new HypreParVector(vert_fes.GetComm(), glob_size, d_z_ptr, cols, dev);
+      }
+      else
+      {
+         z = NULL;
+      }
   }
 }

@@ -33,12 +33,14 @@ protected:
   const int order; ///< Polynomial degree.
   H1_FECollection vert_fec; ///< The corresponding H1 collection.
   ParFiniteElementSpace vert_fes; ///< The corresponding H1 space.
-   Vector *xyz_tvec; ///< Mesh vertex coordinates in true-vector format.
-   HypreParMatrix *G; ///< Discrete gradient matrix.
+   Vector *xyz_tvec = nullptr; ///< Mesh vertex coordinates in true-vector format.
+   HypreParMatrix *G = nullptr; ///< Discrete gradient matrix.
+
+   SparseMatrix G_local;

   /// @name Mesh coordinate vectors in HypreParVector format
   ///@{
-   HypreParVector *x, *y, *z;
+   HypreParVector *x = nullptr, *y = nullptr, *z = nullptr;
   ///@}

   /// @name Construct the local (elementwise) discrete gradient
@@ -89,6 +91,9 @@ public:

   /// Construct the discrete gradient matrix (not part of the public API).
   void FormGradientMatrix();
+
+   void FormGradientMatrixLocal();
+
   ~BatchedLOR_AMS();
 };

@@ -15,6 +15,8 @@
 #include <climits>
 #include "../pbilinearform.hpp"

+#include "../../general/nvtx.hpp"
+
 // Specializations
 #include "lor_h1.hpp"
 #include "lor_nd.hpp"
@@ -70,8 +72,13 @@ bool BatchedLORAssembly::FormIsSupported(BilinearForm &a)
 }

 void BatchedLORAssembly::FormLORVertexCoordinates(FiniteElementSpace &fes_ho,
-                                                  Vector &X_vert)
+                                                  Vector &X_vert,
+                                                  Vector *evec)
 {
+#undef MFEM_NVTX_COLOR
+#define MFEM_NVTX_COLOR DeepSkyBlue
+   NVTX("LOR Coordinates");
+
   Mesh &mesh_ho = *fes_ho.GetMesh();
   mesh_ho.EnsureNodes();

@@ -87,9 +94,22 @@ void BatchedLORAssembly::FormLORVertexCoordinates(FiniteElementSpace &fes_ho,
   const Operator *nodal_restriction =
      nodal_fes->GetElementRestriction(ElementDofOrdering::LEXICOGRAPHIC);

+   Vector *tmp_evec = nullptr;
+   Vector *nodal_evec;
+
+   if (evec)
+   {
+      nodal_evec = evec;
+      nodal_evec->SetSize(nodal_restriction->Height());
+   }
+   else
+   {
+      tmp_evec = new Vector(nodal_restriction->Height());
+      nodal_evec = tmp_evec;
+   }
+
   // Map from nodal L-vector to E-vector
-   Vector nodal_evec(nodal_restriction->Height());
-   nodal_restriction->Mult(*nodal_gf, nodal_evec);
+   nodal_restriction->Mult(*nodal_gf, *nodal_evec);

   IntegrationRule ir = GetCollocatedIntRule(fes_ho);

@@ -98,7 +118,9 @@ void BatchedLORAssembly::FormLORVertexCoordinates(FiniteElementSpace &fes_ho,
   const QuadratureInterpolator *quad_interp =
      nodal_fes->GetQuadratureInterpolator(ir);
   quad_interp->SetOutputLayout(QVectorLayout::byVDIM);
-   quad_interp->Values(nodal_evec, X_vert);
+   quad_interp->Values(*nodal_evec, X_vert);
+
+   delete tmp_evec;
 }

 // The following two functions (GetMinElt and GetAndIncrementNnzIndex) are
@@ -372,11 +394,11 @@ void BatchedLORAssembly::SparseIJToCSR(OperatorHandle &A) const

   A_mat->OverrideSize(nvdof, nvdof);

-   A_mat->GetMemoryI().New(nvdof+1, Device::GetDeviceMemoryType());
+   EnsureCapacity(A_mat->GetMemoryI(), nvdof+1, Device::GetDeviceMemoryType());
   int nnz = FillI(*A_mat);

-   A_mat->GetMemoryJ().New(nnz, Device::GetDeviceMemoryType());
-   A_mat->GetMemoryData().New(nnz, Device::GetDeviceMemoryType());
+   EnsureCapacity(A_mat->GetMemoryJ(), nnz, Device::GetDeviceMemoryType());
+   EnsureCapacity(A_mat->GetMemoryData(), nnz, Device::GetDeviceMemoryType());
   FillJAndData(*A_mat);
 }

@@ -455,7 +477,6 @@ void BatchedLORAssembly::ParAssemble(
   BilinearForm &a, const Array<int> &ess_dofs, OperatorHandle &A)
 {
   // Assemble the system matrix local to this partition
-   OperatorHandle A_local;
   AssembleWithoutBC(a, A_local);

   ParBilinearForm *pa =
@@ -471,6 +492,9 @@ void BatchedLORAssembly::ParAssemble(
 void BatchedLORAssembly::Assemble(
   BilinearForm &a, const Array<int> ess_dofs, OperatorHandle &A)
 {
+#undef MFEM_NVTX_COLOR
+#define MFEM_NVTX_COLOR NavyBlue
+   NVTX("LOR Assemble");
 #ifdef MFEM_USE_MPI
   if (dynamic_cast<ParFiniteElementSpace*>(&fes_ho))
   {
@@ -53,6 +53,8 @@ protected:
   /// nonzero).
   Array<int> sparse_mapping;

+   OperatorHandle A_local; // Cache this
+
 public:
   /// Construct the batched assembly object corresponding to @a fes_ho_.
   BatchedLORAssembly(FiniteElementSpace &fes_ho_);
@@ -68,12 +70,12 @@ public:

   /// Compute the vertices of the LOR mesh and place the result in @a X_vert.
   static void FormLORVertexCoordinates(FiniteElementSpace &fes_ho,
-                                        Vector &X_vert);
+                                        Vector &X_vert,
+                                        Vector *evec = nullptr);

   /// Return the vertices of the LOR mesh in E-vector format
   const Vector &GetLORVertexCoordinates() { return X_vert; }

-protected:
   /// After assembling the "sparse IJ" format, convert it to CSR.
   void SparseIJToCSR(OperatorHandle &A) const;

@@ -117,12 +119,12 @@ public:
 /// If the capacity of @a mem is not large enough, delete it and allocate new
 /// memory with size @a capacity.
 template <typename T>
-void EnsureCapacity(Memory<T> &mem, int capacity)
+void EnsureCapacity(Memory<T> &mem, int capacity, MemoryType mt)
 {
   if (mem.Capacity() < capacity)
   {
      mem.Delete();
-      mem.New(capacity, mem.GetMemoryType());
+      mem.New(capacity, mt);
   }
 }

@@ -205,6 +205,49 @@ void BatchedLOR_H1::Assemble2D()
   }
 }

+template<int ORDER>
+static void SparseMapping3D(Array<int> &sparse_mapping)
+{
+   static constexpr int nnz_per_row = 27;
+   static constexpr int nd1d = ORDER + 1;
+   static constexpr int ndof_per_el = nd1d*nd1d*nd1d;
+
+   sparse_mapping.SetSize(nnz_per_row*ndof_per_el);
+   sparse_mapping = -1;
+   auto map = Reshape(sparse_mapping.HostReadWrite(), nnz_per_row, ndof_per_el);
+
+   for (int iz=0; iz<nd1d; ++iz)
+   {
+      const int jz_begin = (iz > 0) ? iz - 1 : 0;
+      const int jz_end = (iz < ORDER) ? iz + 1 : ORDER;
+      for (int iy=0; iy<nd1d; ++iy)
+      {
+         const int jy_begin = (iy > 0) ? iy - 1 : 0;
+         const int jy_end = (iy < ORDER) ? iy + 1 : ORDER;
+         for (int ix=0; ix<nd1d; ++ix)
+         {
+            const int jx_begin = (ix > 0) ? ix - 1 : 0;
+            const int jx_end = (ix < ORDER) ? ix + 1 : ORDER;
+
+            const int ii_el = ix + nd1d*(iy + nd1d*iz);
+
+            for (int jz=jz_begin; jz<=jz_end; ++jz)
+            {
+               for (int jy=jy_begin; jy<=jy_end; ++jy)
+               {
+                  for (int jx=jx_begin; jx<=jx_end; ++jx)
+                  {
+                     const int jj_off = (jx-ix+1) + 3*(jy-iy+1) + 9*(jz-iz+1);
+                     const int jj_el = jx + nd1d*(jy + nd1d*jz);
+                     map(jj_off, ii_el) = jj_el;
+                  }
+               }
+            }
+         }
+      }
+   }
+}
+
 template <int ORDER>
 void BatchedLOR_H1::Assemble3D()
 {
@@ -330,21 +373,16 @@ void BatchedLOR_H1::Assemble3D()
                  }
               }

-               //MFEM_UNROLL(2)
               for (int iqx=0; iqx<2; ++iqx)
               {
-                  //MFEM_UNROLL(2)
                  for (int jz=0; jz<2; ++jz)
                  {
                     // Note loop starts at iz=jz here, taking advantage of
                     // symmetries.
-                     //MFEM_UNROLL(2)
                     for (int iz=jz; iz<2; ++iz)
                     {
-                        //MFEM_UNROLL(2)
                        for (int iqy=0; iqy<2; ++iqy)
                        {
-                           //MFEM_UNROLL(2)
                           for (int iqz=0; iqz<2; ++iqz)
                           {
                              const double mq = const_mq ? MQ(0,0,0,0) : MQ(kx+iqx, ky+iqy, kz+iqz, iel_ho);
@@ -379,10 +417,8 @@ void BatchedLOR_H1::Assemble3D()
                              double wdetJ = Q(6,iqz,iqy,iqx);
                              mass_A(iqy,iz,jz,iqx) += mq*wdetJ*biz*bjz;
                           }
-                           //MFEM_UNROLL(2)
                           for (int jy=0; jy<2; ++jy)
                           {
-                              //MFEM_UNROLL(2)
                              for (int iy=0; iy<2; ++iy)
                              {
                                 const double biy = (iy == iqy) ? 1.0 : 0.0;
@@ -405,16 +441,12 @@ void BatchedLOR_H1::Assemble3D()
                              }
                           }
                        }
-                        //MFEM_UNROLL(2)
                        for (int jy=0; jy<2; ++jy)
                        {
-                           //MFEM_UNROLL(2)
                           for (int jx=0; jx<2; ++jx)
                           {
-                              //MFEM_UNROLL(2)
                              for (int iy=0; iy<2; ++iy)
                              {
-                                 //MFEM_UNROLL(2)
                                 for (int ix=0; ix<2; ++ix)
                                 {
                                    const double bix = (ix == iqx) ? 1.0 : 0.0;
@@ -482,40 +514,176 @@ void BatchedLOR_H1::Assemble3D()
         }
      }
   });
+   SparseMapping3D<ORDER>(sparse_mapping);
+}

-   sparse_mapping.SetSize(nnz_per_row*ndof_per_el);
-   sparse_mapping = -1;
-   auto map = Reshape(sparse_mapping.HostReadWrite(), nnz_per_row, ndof_per_el);
-   for (int iz=0; iz<nd1d; ++iz)
+template <>
+void BatchedLOR_H1::Assemble3D<1>()
+{
+   static constexpr int nv = 8;
+   static constexpr int nd1d = 2;
+   static constexpr int ndof_per_el = 8;
+   static constexpr int nnz_per_row = 27;
+   static constexpr int sz_local_mat = nv*nv;
+
+   const int nel_ho = fes_ho.GetNE();
+
+   const bool const_mq = c1.Size() == 1;
+   const auto MQ = const_mq
+                   ? Reshape(c1.Read(), 1, 1, 1, 1)
+                   : Reshape(c1.Read(), nd1d, nd1d, nd1d, nel_ho);
+   const bool const_dq = c2.Size() == 1;
+   const auto DQ = const_dq
+                   ? Reshape(c2.Read(), 1, 1, 1, 1)
+                   : Reshape(c2.Read(), nd1d, nd1d, nd1d, nel_ho);
+
+   sparse_ij.SetSize(nel_ho*ndof_per_el*nnz_per_row);
+   auto V = Reshape(sparse_ij.Write(), nnz_per_row, nd1d, nd1d, nd1d, nel_ho);
+
+   const auto X = X_vert.Read();
+
+   MFEM_FORALL_3D(iel_ho, nel_ho, 8, 4, 1,
   {
-      const int jz_begin = (iz > 0) ? iz - 1 : 0;
-      const int jz_end = (iz < ORDER) ? iz + 1 : ORDER;
-      for (int iy=0; iy<nd1d; ++iy)
+      static constexpr int e[8] = {0,1,3,2,4,5,7,6};
+      MFEM_SHARED double vx[8], vy[8], vz[8];
+      const int tidz = MFEM_THREAD_ID(z);
+
+      MFEM_SHARED double local_mat_[sz_local_mat];
+      DeviceTensor<4> local_mat(local_mat_, 2,2,2, nv);
+
+      if (tidz == 0)
      {
-         const int jy_begin = (iy > 0) ? iy - 1 : 0;
-         const int jy_end = (iy < ORDER) ? iy + 1 : ORDER;
-         for (int ix=0; ix<nd1d; ++ix)
+         MFEM_FOREACH_THREAD(xyz,x,8)
         {
-            const int jx_begin = (ix > 0) ? ix - 1 : 0;
-            const int jx_end = (ix < ORDER) ? ix + 1 : ORDER;
-
-            const int ii_el = ix + nd1d*(iy + nd1d*iz);
-
-            for (int jz=jz_begin; jz<=jz_end; ++jz)
+            const int z = xyz%2, y = (xyz/2)%2, x = xyz/2/2;
+            MFEM_FOREACH_THREAD(j,y,nnz_per_row)
            {
-               for (int jy=jy_begin; jy<=jy_end; ++jy)
+               if (j < 8) { local_mat(z,y,x,j) = 0.0; }
+               V(j,x,y,z,iel_ho) = 0.0;
+               if (j == 0)
               {
-                  for (int jx=jx_begin; jx<=jx_end; ++jx)
-                  {
-                     const int jj_off = (jx-ix+1) + 3*(jy-iy+1) + 9*(jz-iz+1);
-                     const int jj_el = jx + nd1d*(jy + nd1d*jz);
-                     map(jj_off, ii_el) = jj_el;
-                  }
+                  const int i = x + 2*y + 4*z;
+                  const int ei = 3*(e[i] + 8*iel_ho);
+                  vx[i] = X[ei + 0];
+                  vy[i] = X[ei + 1];
+                  vz[i] = X[ei + 2];
               }
            }
         }
      }
-   }
+      MFEM_SYNC_THREAD;
+
+      MFEM_FOREACH_THREAD(xyz,x,8)
+      {
+         const int qz = xyz%2, qy = (xyz/2)%2, qx = xyz/2/2;
+         static constexpr double w = 1.0/8.0;
+
+         double J_[3*3];
+         DeviceTensor<2> J(J_, 3,3);
+         Jacobian3D(qx,qy,qz, vx,vy,vz, J);
+         const double detJ = Det3D(J);
+         const double w_detJ = w/detJ;
+
+         // adj(J)
+         double A_[3*3];
+         DeviceTensor<2> A(A_, 3, 3);
+         Adjugate3D(J, A);
+
+         const double J11 = w_detJ*(A(0,0)*A(0,0)+A(0,1)*A(0,1)+A(0,2)*A(0,2)); // 1,1
+         const double J21 = w_detJ*(A(0,0)*A(1,0)+A(0,1)*A(1,1)+A(0,2)*A(1,2)); // 2,1
+         const double J31 = w_detJ*(A(0,0)*A(2,0)+A(0,1)*A(2,1)+A(0,2)*A(2,2)); // 3,1
+         const double J12 = J21;
+         const double J22 = w_detJ*(A(1,0)*A(1,0)+A(1,1)*A(1,1)+A(1,2)*A(1,2)); // 2,2
+         const double J32 = w_detJ*(A(1,0)*A(2,0)+A(1,1)*A(2,1)+A(1,2)*A(2,2)); // 3,2
+         const double J13 = J31;
+         const double J23 = J32;
+         const double J33 = w_detJ*(A(2,0)*A(2,0)+A(2,1)*A(2,1)+A(2,2)*A(2,2)); // 3,3
+         const double wdetJ =  w*detJ;
+
+         const double mq = const_mq ? MQ(0,0,0,0) : MQ(qx,qy,qz, iel_ho);
+         const double dq = const_dq ? DQ(0,0,0,0) : DQ(qx,qy,qz, iel_ho);
+
+         MFEM_FOREACH_THREAD(xyz,y,8)
+         {
+            const int jz = xyz%2, jy = (xyz/2)%2, jx = xyz/2/2;
+            const double bjz = (jz == qz) ? 1.0 : 0.0;
+            const double gjz = (jz == 0) ? -1.0 : 1.0;
+            const double bjy = (jy == qy) ? 1.0 : 0.0;
+            const double gjy = (jy == 0) ? -1.0 : 1.0;
+            const double bjx = (jx == qx) ? 1.0 : 0.0;
+            const double gjx = (jx == 0) ? -1.0 : 1.0;
+            const double djx = gjx*bjy*bjz;
+            const double djy = bjx*gjy*bjz;
+            const double djz = bjx*bjy*gjz;
+            const int jj_loc = jx + 2*jy + 4*jz;
+
+            MFEM_FOREACH_THREAD(xyz,z,8)
+            {
+               const int iz = xyz%2, iy = (xyz/2)%2, ix = xyz/2/2;
+               const double biz = (iz == qz) ? 1.0 : 0.0;
+               const double giz = (iz == 0) ? -1.0 : 1.0;
+               const double biy = (iy == qy) ? 1.0 : 0.0;
+               const double giy = (iy == 0) ? -1.0 : 1.0;
+               const double bix = (ix == qx) ? 1.0 : 0.0;
+               const double gix = (ix == 0) ? -1.0 : 1.0;
+               const double dix = gix*biy*biz;
+               const double diy = bix*giy*biz;
+               const double diz = bix*biy*giz;
+               const int ii_loc = ix + 2*iy + 4*iz;
+
+               // Only store the lower-triangular part of
+               // the matrix (by symmetry).
+               if (jj_loc > ii_loc) { continue; }
+
+               double grad_grad = 0.0;
+               grad_grad += dix*djx*J11;
+               grad_grad += diy*djx*J12;
+               grad_grad += diz*djx*J13;
+
+               grad_grad += dix*djy*J21;
+               grad_grad += diy*djy*J22;
+               grad_grad += diz*djy*J23;
+
+               grad_grad += dix*djz*J31;
+               grad_grad += diy*djz*J32;
+               grad_grad += diz*djz*J33;
+
+               const double basis_basis = wdetJ*bix*biy*biz*bjx*bjy*bjz;
+               const double value = dq*grad_grad + mq*basis_basis;
+               AtomicAdd(local_mat(iz,iy,ix, jj_loc), value);
+            } // i
+         } // j
+      } // q
+      MFEM_SYNC_THREAD;
+
+      // Assemble the local matrix into the macro-element sparse matrix
+      // in a format similar to coordinate format.
+      // The (I,J) arrays are implicit (not stored explicitly).
+      if (tidz == 0)
+      {
+         MFEM_FOREACH_THREAD(xyz,x,8)
+         {
+            const int iz = xyz%2, iy = (xyz/2)%2, ix = xyz/2/2;
+            const int ii_loc = ix + 2*iy + 4*iz;
+            MFEM_FOREACH_THREAD(jj_loc,y,8)
+            {
+               const int jx = jj_loc%2, jy = (jj_loc/2)%2, jz = jj_loc/2/2;
+               const int jj_off = (jx-ix+1) + 3*(jy-iy+1) + 9*(jz-iz+1);
+
+               if (jj_loc <= ii_loc)
+               {
+                  AtomicAdd(V(jj_off, ix,iy,iz, iel_ho), local_mat(iz,iy,ix, jj_loc));
+               }
+               else
+               {
+                  AtomicAdd(V(jj_off, ix,iy,iz, iel_ho), local_mat(jz,jy,jx, ii_loc));
+               }
+            }
+         }
+      }
+      MFEM_SYNC_THREAD;
+   });
+   SparseMapping3D<1>(sparse_mapping);
 }

 // Explicit template instantiations
@@ -528,7 +696,7 @@ template void BatchedLOR_H1::Assemble2D<6>();
 template void BatchedLOR_H1::Assemble2D<7>();
 template void BatchedLOR_H1::Assemble2D<8>();

-template void BatchedLOR_H1::Assemble3D<1>();
+//template void BatchedLOR_H1::Assemble3D<1>(); // explicitly specialized
 template void BatchedLOR_H1::Assemble3D<2>();
 template void BatchedLOR_H1::Assemble3D<3>();
 template void BatchedLOR_H1::Assemble3D<4>();
@@ -216,6 +216,90 @@ void BatchedLOR_ND::Assemble2D()
   }
 }

+
+template<int ORDER>
+static void SparseMapping3D(Array<int> &sparse_mapping)
+{
+   static constexpr int nnz_per_row = 33;
+   static constexpr int dim = 3;
+   static constexpr int o = ORDER;
+   static constexpr int op1 = ORDER + 1;
+   static constexpr int ndof_per_el = dim*o*op1*op1;
+
+   sparse_mapping.SetSize(nnz_per_row*ndof_per_el);
+   sparse_mapping = -1;
+   auto map = Reshape(sparse_mapping.HostReadWrite(), nnz_per_row, ndof_per_el);
+
+   for (int ci=0; ci<dim; ++ci)
+   {
+      const int i_off = ci*o*op1*op1;
+      const int id0 = ci;
+      const int id1 = (ci+1)%3;
+      const int id2 = (ci+2)%3;
+
+      const int nxi = (ci == 0) ? o : op1;
+      const int nyi = (ci == 1) ? o : op1;
+
+      for (int i0=0; i0<o; ++i0)
+      {
+         for (int i1=0; i1<op1; ++i1)
+         {
+            for (int i2=0; i2<op1; ++i2)
+            {
+               int ii_lex[3];
+               ii_lex[id0] = i0;
+               ii_lex[id1] = i1;
+               ii_lex[id2] = i2;
+               const int ii_el = i_off + ii_lex[0] + ii_lex[1]*nxi + ii_lex[2]*nxi*nyi;
+
+               for (int cj_rel=0; cj_rel<dim; ++cj_rel)
+               {
+                  const int cj = (ci + cj_rel) % 3;
+                  const int j_off = cj*o*op1*op1;
+
+                  const int nxj = (cj == 0) ? o : op1;
+                  const int nyj = (cj == 1) ? o : op1;
+
+                  const int j0_begin = i0;
+                  const int j0_end = (cj_rel == 0) ? i0 : i0 + 1;
+                  const int j1_begin = (i1 > 0) ? i1-1 : i1;
+                  const int j1_end = (cj_rel == 1)
+                                     ? ((i1 < o) ? i1 : i1-1)
+                                     : ((i1 < o) ? i1+1 : i1);
+                  const int j2_begin = (i2 > 0) ? i2-1 : i2;
+                  const int j2_end = (cj_rel == 2)
+                                     ? ((i2 < o) ? i2 : i2-1)
+                                     : ((i2 < o) ? i2+1 : i2);
+
+                  for (int j0=j0_begin; j0<=j0_end; ++j0)
+                  {
+                     const int d0 = j0 - i0;
+                     for (int j1=j1_begin; j1<=j1_end; ++j1)
+                     {
+                        const int d1 = j1 - i1 + 1;
+                        for (int j2=j2_begin; j2<=j2_end; ++j2)
+                        {
+                           const int d2 = j2 - i2 + 1;
+                           int jj_lex[3];
+                           jj_lex[id0] = j0;
+                           jj_lex[id1] = j1;
+                           jj_lex[id2] = j2;
+                           const int jj_el = j_off + jj_lex[0] + jj_lex[1]*nxj + jj_lex[2]*nxj*nyj;
+                           int jj_off;
+                           if (cj_rel == 0) { jj_off = d1 + 3*d2; }
+                           else if (cj_rel == 1) { jj_off = 9 + d0 + 2*d1 + 4*d2; }
+                           else /* if (cj_rel == 2) */ { jj_off = 21 + d0 + 2*d1 + 6*d2; }
+                           map(jj_off, ii_el) = jj_el;
+                        }
+                     }
+                  }
+               }
+            }
+         }
+      }
+   }
+}
+
 template <int ORDER>
 void BatchedLOR_ND::Assemble3D()
 {
@@ -244,7 +328,7 @@ void BatchedLOR_ND::Assemble3D()
   sparse_ij.SetSize(nnz_per_row*ndof_per_el*nel_ho);
   auto V = Reshape(sparse_ij.Write(), nnz_per_row, o*op1*op1, dim, nel_ho);

-   auto X = X_vert.Read();
+   const auto X = X_vert.Read();

   // Last thread dimension is lowered to avoid "too many resources" error
   MFEM_FORALL_3D(iel_ho, nel_ho, ORDER, ORDER, (ORDER>6)?4:ORDER,
@@ -488,78 +572,231 @@ void BatchedLOR_ND::Assemble3D()
         }
      }
   });
+   SparseMapping3D<ORDER>(sparse_mapping);
+}

-   sparse_mapping.SetSize(nnz_per_row*ndof_per_el);
-   sparse_mapping = -1;
-   auto map = Reshape(sparse_mapping.HostReadWrite(), nnz_per_row, ndof_per_el);
-   for (int ci=0; ci<dim; ++ci)
+template<>
+void BatchedLOR_ND::Assemble3D<1>()
+{
+   static constexpr int ne = 12; // number of edges in hexahedron
+   static constexpr int dim = 3;
+   static constexpr int o = 1;
+   static constexpr int op1 = 2;
+   static constexpr int ndof_per_el = dim*o*op1*op1;
+   static constexpr int nnz_per_row = 33;
+   static constexpr int sz_local_mat = ne*ne;
+
+   const int nel_ho = fes_ho.GetNE();
+   const bool const_mq = c1.Size() == 1;
+   const auto MQ = const_mq
+                   ? Reshape(c1.Read(), 1, 1, 1, 1)
+                   : Reshape(c1.Read(), op1, op1, op1, nel_ho);
+   const bool const_dq = c2.Size() == 1;
+   const auto DQ = const_dq
+                   ? Reshape(c2.Read(), 1, 1, 1, 1)
+                   : Reshape(c2.Read(), op1, op1, op1, nel_ho);
+
+   sparse_ij.SetSize(nnz_per_row*ndof_per_el*nel_ho);
+   auto V = Reshape(sparse_ij.Write(), nnz_per_row, o*op1*op1, dim, nel_ho);
+
+   const auto X = X_vert.Read();
+
+   MFEM_FORALL_3D(iel_ho, nel_ho, 8, 1, 4,
   {
-      const int i_off = ci*o*op1*op1;
-      const int id0 = ci;
-      const int id1 = (ci+1)%3;
-      const int id2 = (ci+2)%3;
-
-      const int nxi = (ci == 0) ? o : op1;
-      const int nyi = (ci == 1) ? o : op1;
-
-      for (int i0=0; i0<o; ++i0)
+      MFEM_FOREACH_THREAD(iz,z,o) // 1
      {
-         for (int i1=0; i1<op1; ++i1)
+         MFEM_FOREACH_THREAD(iy,y,op1) // 2
         {
-            for (int i2=0; i2<op1; ++i2)
+            MFEM_FOREACH_THREAD(ix,x,op1) // 2
            {
-               int ii_lex[3];
-               ii_lex[id0] = i0;
-               ii_lex[id1] = i1;
-               ii_lex[id2] = i2;
-               const int ii_el = i_off + ii_lex[0] + ii_lex[1]*nxi + ii_lex[2]*nxi*nyi;
-
-               for (int cj_rel=0; cj_rel<dim; ++cj_rel)
+               for (int c=0; c<dim; ++c)
               {
-                  const int cj = (ci + cj_rel) % 3;
-                  const int j_off = cj*o*op1*op1;
-
-                  const int nxj = (cj == 0) ? o : op1;
-                  const int nyj = (cj == 1) ? o : op1;
-
-                  const int j0_begin = i0;
-                  const int j0_end = (cj_rel == 0) ? i0 : i0 + 1;
-                  const int j1_begin = (i1 > 0) ? i1-1 : i1;
-                  const int j1_end = (cj_rel == 1)
-                                     ? ((i1 < o) ? i1 : i1-1)
-                                     : ((i1 < o) ? i1+1 : i1);
-                  const int j2_begin = (i2 > 0) ? i2-1 : i2;
-                  const int j2_end = (cj_rel == 2)
-                                     ? ((i2 < o) ? i2 : i2-1)
-                                     : ((i2 < o) ? i2+1 : i2);
-
-                  for (int j0=j0_begin; j0<=j0_end; ++j0)
+                  for (int j=0; j<nnz_per_row; ++j)
                  {
-                     const int d0 = j0 - i0;
-                     for (int j1=j1_begin; j1<=j1_end; ++j1)
-                     {
-                        const int d1 = j1 - i1 + 1;
-                        for (int j2=j2_begin; j2<=j2_end; ++j2)
-                        {
-                           const int d2 = j2 - i2 + 1;
-                           int jj_lex[3];
-                           jj_lex[id0] = j0;
-                           jj_lex[id1] = j1;
-                           jj_lex[id2] = j2;
-                           const int jj_el = j_off + jj_lex[0] + jj_lex[1]*nxj + jj_lex[2]*nxj*nyj;
-                           int jj_off;
-                           if (cj_rel == 0) { jj_off = d1 + 3*d2; }
-                           else if (cj_rel == 1) { jj_off = 9 + d0 + 2*d1 + 4*d2; }
-                           else /* if (cj_rel == 2) */ { jj_off = 21 + d0 + 2*d1 + 6*d2; }
-                           map(jj_off, ii_el) = jj_el;
-                        }
-                     }
+                     V(j,ix+iy*op1+iz*op1*op1,c,iel_ho) = 0.0;
                  }
               }
            }
         }
      }
-   }
+      MFEM_SYNC_THREAD;
+
+      MFEM_SHARED double local_mat_[sz_local_mat];
+      DeviceTensor<4> local_mat(local_mat_, 3,4, 3,4);
+      /// should be optimized
+      for (int i=0; i<sz_local_mat; ++i) { local_mat[i] = 0.0; }
+
+      MFEM_SHARED double vx[8], vy[8], vz[8];
+      /// should be optimized
+      LORVertexCoordinates3D<1>(X, iel_ho, 0,0,0, vx,vy,vz);
+      MFEM_SYNC_THREAD;
+
+      MFEM_FOREACH_THREAD(xyz,x,8)
+      {
+         const int qz = xyz%2, qy = (xyz/2)%2, qx = xyz/2/2;
+         static constexpr double w = 1.0/8.0;
+
+         double J_[3*3];
+         DeviceTensor<2> J(J_, 3,3);
+         Jacobian3D(qx,qy,qz, vx,vy,vz, J);
+         const double detJ = Det3D(J);
+         const double w_detJ = w/detJ;
+
+         // adj(J)
+         double A_[3*3];
+         DeviceTensor<2> A(A_, 3,3);
+         Adjugate3D(J, A);
+
+         const double Q0 = w_detJ*(A(0,0)*A(0,0)+A(0,1)*A(0,1)+A(0,2)*A(0,2)); // 1,1
+         const double Q1 = w_detJ*(A(0,0)*A(1,0)+A(0,1)*A(1,1)+A(0,2)*A(1,2)); // 2,1
+         const double Q2 = w_detJ*(A(0,0)*A(2,0)+A(0,1)*A(2,1)+A(0,2)*A(2,2)); // 3,1
+         const double Q3 = w_detJ*(A(1,0)*A(1,0)+A(1,1)*A(1,1)+A(1,2)*A(1,2)); // 2,2
+         const double Q4 = w_detJ*(A(1,0)*A(2,0)+A(1,1)*A(2,1)+A(1,2)*A(2,2)); // 3,2
+         const double Q5 = w_detJ*(A(2,0)*A(2,0)+A(2,1)*A(2,1)+A(2,2)*A(2,2)); // 3,3
+
+         // w J^T J / det(J)
+         const double Q6  = w_detJ*(J(0,0)*J(0,0)+J(1,0)*J(1,0)+J(2,0)*J(2,0)); // 1,1
+         const double Q7  = w_detJ*(J(0,0)*J(0,1)+J(1,0)*J(1,1)+J(2,0)*J(2,1)); // 2,1
+         const double Q8  = w_detJ*(J(0,0)*J(0,2)+J(1,0)*J(1,2)+J(2,0)*J(2,2)); // 3,1
+         const double Q9  = w_detJ*(J(0,1)*J(0,1)+J(1,1)*J(1,1)+J(2,1)*J(2,1)); // 2,2
+         const double Q10 = w_detJ*(J(0,1)*J(0,2)+J(1,1)*J(1,2)+J(2,1)*J(2,2)); // 3,2
+         const double Q11 = w_detJ*(J(0,2)*J(0,2)+J(1,2)*J(1,2)+J(2,2)*J(2,2)); // 3,3
+
+         const double mq = const_mq ? MQ(0,0,0,0) : MQ(qx,qy,qz, iel_ho);
+         const double dq = const_dq ? DQ(0,0,0,0) : DQ(qx,qy,qz, iel_ho);
+
+         // Loop over x,y,z components. 0 => x, 1 => y, 2 => z
+         MFEM_FOREACH_THREAD(cj,y,dim)
+         {
+            const double jq1 = (cj == 0) ? qy : ((cj == 1) ? qz : qx);
+            const double jq2 = (cj == 0) ? qz : ((cj == 1) ? qx : qy);
+
+            const int jd_0 = cj;
+            const int jd_1 = (cj + 1)%3;
+            const int jd_2 = (cj + 2)%3;
+
+            MFEM_FOREACH_THREAD(bj,z,4) // 4 edges in each dim
+            {
+               const int bj1 = bj%2;
+               const int bj2 = bj/2;
+
+               double curl_j[3];
+               curl_j[jd_0] = 0.0;
+               curl_j[jd_1] = ((bj1 == 0) ? jq1 - 1 : -jq1)*((bj2 == 0) ? 1 : -1);
+               curl_j[jd_2] = ((bj2 == 0) ? 1 - jq2 : jq2)*((bj1 == 0) ? 1 : -1);
+
+               double basis_j[3];
+               basis_j[jd_0] = ((bj1 == 0) ? 1 - jq1 : jq1)*((bj2 == 0) ? 1 - jq2 : jq2);
+               basis_j[jd_1] = 0.0;
+               basis_j[jd_2] = 0.0;
+
+               const int jj_loc = bj + 4*cj;
+
+               for (int ci=0; ci<dim; ++ci)
+               {
+                  const double iq1 = (ci == 0) ? qy : ((ci == 1) ? qz : qx);
+                  const double iq2 = (ci == 0) ? qz : ((ci == 1) ? qx : qy);
+                  const int id_0 = ci, id_1 = (ci + 1)%3, id_2 = (ci + 2)%3;
+
+                  for (int bi=0; bi<4; ++bi)
+                  {
+                     const int bi1 = bi%2, bi2 = bi/2;
+
+                     double curl_i[3];
+                     curl_i[id_0] = 0.0;
+                     curl_i[id_1] = ((bi1 == 0) ? iq1 - 1 : -iq1)*((bi2 == 0) ? 1 : -1);
+                     curl_i[id_2] = ((bi2 == 0) ? 1 - iq2 : iq2)*((bi1 == 0) ? 1 : -1);
+
+                     double basis_i[3];
+                     basis_i[id_0] = ((bi1 == 0) ? 1 - iq1 : iq1)*((bi2 == 0) ? 1 - iq2 : iq2);
+                     basis_i[id_1] = 0.0;
+                     basis_i[id_2] = 0.0;
+
+                     const int ii_loc = bi + 4*ci;
+
+                     // Only store the lower-triangular part of
+                     // the matrix (by symmetry).
+                     if (jj_loc > ii_loc) { continue; }
+
+                     double curl_curl = 0.0;
+                     curl_curl += Q6*curl_i[0]*curl_j[0];
+                     curl_curl += Q7*(curl_i[0]*curl_j[1] + curl_i[1]*curl_j[0]);
+                     curl_curl += Q8*(curl_i[0]*curl_j[2] + curl_i[2]*curl_j[0]);
+                     curl_curl += Q9*curl_i[1]*curl_j[1];
+                     curl_curl += Q10*(curl_i[1]*curl_j[2] + curl_i[2]*curl_j[1]);
+                     curl_curl += Q11*curl_i[2]*curl_j[2];
+
+                     double basis_basis = 0.0;
+                     basis_basis += Q0*basis_i[0]*basis_j[0];
+                     basis_basis += Q1*(basis_i[0]*basis_j[1] + basis_i[1]*basis_j[0]);
+                     basis_basis += Q2*(basis_i[0]*basis_j[2] + basis_i[2]*basis_j[0]);
+                     basis_basis += Q3*basis_i[1]*basis_j[1];
+                     basis_basis += Q4*(basis_i[1]*basis_j[2] + basis_i[2]*basis_j[1]);
+                     basis_basis += Q5*basis_i[2]*basis_j[2];
+
+                     const double val = dq*curl_curl + mq*basis_basis;
+                     AtomicAdd(local_mat(ci,bi, cj,bj), val);
+                  } // bi
+               } // ci
+            } // bj
+         } // cj
+      } // q
+      MFEM_SYNC_THREAD;
+
+      // Assemble the local matrix into the macro-element sparse matrix
+      // The nonzeros of the macro-element sparse matrix are ordered as
+      // follows:
+      //
+      // The axes are ordered relative to the direction of the basis
+      // vector, e.g. for x-vectors, the axes are (x,y,z), for
+      // y-vectors the axes are (y,z,x), and for z-vectors the axes are
+      // (z,x,y).
+      //
+      // The nonzeros are then given in "rotated lexicographic"
+      // ordering, according to these axes.
+      if (MFEM_THREAD_ID(x) == 0 && MFEM_THREAD_ID(y) == 0)
+      {
+         MFEM_FOREACH_THREAD(ii_loc,z,ne)
+         {
+            const int ci = ii_loc/4, bi = ii_loc%4;
+            const int i0 = 0, i1 = bi%2, i2 = bi/2;
+            const int id0 = ci, id1 = (ci+1)%3, id2 = (ci+2)%3;
+
+            int ii_lex[3];
+            ii_lex[id0] = i0, ii_lex[id1] = i1, ii_lex[id2] = i2;
+
+            const int nx = (ci == 0) ? o : op1, ny = (ci == 1) ? o : op1;
+            const int ii = ii_lex[0] + (ii_lex[1])*nx + (ii_lex[2])*nx*ny;
+
+            for (int jj_loc=0; jj_loc<ne; ++jj_loc)
+            {
+               const int cj = jj_loc/4, bj = jj_loc%4;
+               // add 3 to take modulus (rather than remainder) when
+               // (cj - ci) is negative
+               const int cj_rel = (3 + cj - ci)%3;
+               const int jd0 = cj_rel, jd1 = (cj_rel+1)%3, jd2 = (cj_rel+2)%3;
+
+               int jj_rel[3];
+               jj_rel[jd0] = 0, jj_rel[jd1] = bj%2, jj_rel[jd2] = bj/2;
+
+               const int d0 = jj_rel[0] - i0;
+               const int d1 = 1 + jj_rel[1] - i1;
+               const int d2 = 1 + jj_rel[2] - i2;
+
+               const int jj_off = (cj_rel == 0) ? d1 + 3*d2 :
+                                  (cj_rel == 1) ? 9 + d0 + 2*d1 + 4*d2 :
+                                  (cj_rel == 2) ? 21 + d0 + 2*d1 + 6*d2 : -1;
+
+               // Symmetry
+               const double val = (jj_loc <= ii_loc)
+                                  ? local_mat(ci,bi, cj,bj)
+                                  : local_mat(cj,bj, ci,bi);
+               AtomicAdd(V(jj_off, ii, ci, iel_ho), val);
+            }
+         }
+      }
+   });
+   SparseMapping3D<1>(sparse_mapping);
 }

 // Explicit template instantiations
@@ -572,7 +809,7 @@ template void BatchedLOR_ND::Assemble2D<6>();
 template void BatchedLOR_ND::Assemble2D<7>();
 template void BatchedLOR_ND::Assemble2D<8>();

-template void BatchedLOR_ND::Assemble3D<1>();
+//template void BatchedLOR_ND::Assemble3D<1>();  // explicitly specialized
 template void BatchedLOR_ND::Assemble3D<2>();
 template void BatchedLOR_ND::Assemble3D<3>();
 template void BatchedLOR_ND::Assemble3D<4>();
@@ -233,6 +233,87 @@ void BatchedLOR_RT::Assemble2D()
   }
 }

+template<int ORDER>
+static void SparseMapping3D(Array<int> &sparse_mapping)
+{
+   static constexpr int nnz_per_row = 11;
+   static constexpr int dim = 3;
+   static constexpr int o = ORDER;
+   static constexpr int op1 = ORDER + 1;
+   static constexpr int ndof_per_el = dim*o*o*op1;
+
+   sparse_mapping.SetSize(nnz_per_row*ndof_per_el);
+   sparse_mapping = -1;
+   auto map = Reshape(sparse_mapping.HostReadWrite(), nnz_per_row, ndof_per_el);
+
+   for (int ci=0; ci<dim; ++ci)
+   {
+      const int i_off = ci*o*o*op1;
+      const int id0 = ci;
+      const int id1 = (ci+1)%3;
+      const int id2 = (ci+2)%3;
+
+      const int nxi = (ci == 0) ? op1 : o;
+      const int nyi = (ci == 1) ? op1 : o;
+
+      for (int i0=0; i0<op1; ++i0)
+      {
+         for (int i1=0; i1<o; ++i1)
+         {
+            for (int i2=0; i2<o; ++i2)
+            {
+               int ii_lex[3];
+               ii_lex[id0] = i0;
+               ii_lex[id1] = i1;
+               ii_lex[id2] = i2;
+               const int ii_el = i_off + ii_lex[0] + ii_lex[1]*nxi + ii_lex[2]*nxi*nyi;
+
+               for (int cj_rel=0; cj_rel<dim; ++cj_rel)
+               {
+                  const int cj = (ci + cj_rel) % 3;
+                  const int j_off = cj*o*o*op1;
+
+                  const int nxj = (cj == 0) ? op1 : o;
+                  const int nyj = (cj == 1) ? op1 : o;
+
+                  const int j0_begin = (i0 > 0) ? i0-1 : i0;
+                  const int j0_end = (cj_rel == 0)
+                                     ? ((i0 < o) ? i0+1 : i0)
+                                     : ((i0 < o) ? i0 : i0-1);
+                  const int j1_begin = i1;
+                  const int j1_end = (cj_rel == 1) ? i1+1 : i1;
+                  const int j2_begin = i2;
+                  const int j2_end = (cj_rel == 2) ? i2+1 : i2;
+
+                  for (int j0=j0_begin; j0<=j0_end; ++j0)
+                  {
+                     const int d0 = 1 + j0 - i0;
+                     for (int j1=j1_begin; j1<=j1_end; ++j1)
+                     {
+                        const int d1 = j1 - i1;
+                        for (int j2=j2_begin; j2<=j2_end; ++j2)
+                        {
+                           const int d2 = j2 - i2;
+                           int jj_lex[3];
+                           jj_lex[id0] = j0;
+                           jj_lex[id1] = j1;
+                           jj_lex[id2] = j2;
+                           const int jj_el = j_off + jj_lex[0] + jj_lex[1]*nxj + jj_lex[2]*nxj*nyj;
+                           int jj_off;
+                           if (cj_rel == 0) { jj_off = d0; }
+                           else if (cj_rel == 1) { jj_off = 3 + d0 + 2*d1; }
+                           else /* if (cj_rel == 2) */ { jj_off = 7 + d0 + 2*d2; }
+                           map(jj_off, ii_el) = jj_el;
+                        }
+                     }
+                  }
+               }
+            }
+         }
+      }
+   }
+}
+
 template <int ORDER>
 void BatchedLOR_RT::Assemble3D()
 {
@@ -261,7 +342,7 @@ void BatchedLOR_RT::Assemble3D()
   sparse_ij.SetSize(nnz_per_row*ndof_per_el*nel_ho);
   auto V = Reshape(sparse_ij.Write(), nnz_per_row, o*o*op1, dim, nel_ho);

-   auto X = X_vert.Read();
+   const auto X = X_vert.Read();

   // Last thread dimension is lowered to avoid "too many resources" error
   MFEM_FORALL_3D(iel_ho, nel_ho, ORDER, ORDER, (ORDER>6)?4:ORDER,
@@ -474,76 +555,185 @@ void BatchedLOR_RT::Assemble3D()
         }
      }
   });
+   SparseMapping3D<ORDER>(sparse_mapping);
+}

-   sparse_mapping.SetSize(nnz_per_row*ndof_per_el);
-   sparse_mapping = -1;
-   auto map = Reshape(sparse_mapping.HostReadWrite(), nnz_per_row, ndof_per_el);
-   for (int ci=0; ci<dim; ++ci)
+template <>
+void BatchedLOR_RT::Assemble3D<1>()
+{
+
+   static constexpr int o = 1;
+   static constexpr int op1 = 2;
+   static constexpr int nf = 6; // number of faces in hexahedron
+   static constexpr int dim = 3;
+   static constexpr int ndof_per_el = dim*o*o*op1;
+   static constexpr int nnz_per_row = 11;
+   static constexpr int sz_local_mat = nf*nf;
+
+   const int nel_ho = fes_ho.GetNE();
+   const bool const_mq = c1.Size() == 1;
+   const auto MQ = const_mq
+                   ? Reshape(c1.Read(), 1, 1, 1, 1)
+                   : Reshape(c1.Read(), op1, op1, op1, nel_ho);
+   const bool const_dq = c2.Size() == 1;
+   const auto DQ = const_dq
+                   ? Reshape(c2.Read(), 1, 1, 1, 1)
+                   : Reshape(c2.Read(), op1, op1, op1, nel_ho);
+
+   sparse_ij.SetSize(nnz_per_row*ndof_per_el*nel_ho);
+   auto V = Reshape(sparse_ij.Write(), nnz_per_row, o*o*op1, dim, nel_ho);
+
+   const auto X = X_vert.Read();
+
+   MFEM_FORALL_3D(iel_ho, nel_ho, 8, 1, 1,
   {
-      const int i_off = ci*o*o*op1;
-      const int id0 = ci;
-      const int id1 = (ci+1)%3;
-      const int id2 = (ci+2)%3;
-
-      const int nxi = (ci == 0) ? op1 : o;
-      const int nyi = (ci == 1) ? op1 : o;
-
-      for (int i0=0; i0<op1; ++i0)
+      MFEM_FOREACH_THREAD(j,x,nnz_per_row)
      {
-         for (int i1=0; i1<o; ++i1)
+         for (int ix = 0; ix < op1; ++ix)
         {
-            for (int i2=0; i2<o; ++i2)
+            for (int c = 0; c < dim; ++c)
            {
-               int ii_lex[3];
-               ii_lex[id0] = i0;
-               ii_lex[id1] = i1;
-               ii_lex[id2] = i2;
-               const int ii_el = i_off + ii_lex[0] + ii_lex[1]*nxi + ii_lex[2]*nxi*nyi;
-
-               for (int cj_rel=0; cj_rel<dim; ++cj_rel)
-               {
-                  const int cj = (ci + cj_rel) % 3;
-                  const int j_off = cj*o*o*op1;
-
-                  const int nxj = (cj == 0) ? op1 : o;
-                  const int nyj = (cj == 1) ? op1 : o;
-
-                  const int j0_begin = (i0 > 0) ? i0-1 : i0;
-                  const int j0_end = (cj_rel == 0)
-                                     ? ((i0 < o) ? i0+1 : i0)
-                                     : ((i0 < o) ? i0 : i0-1);
-                  const int j1_begin = i1;
-                  const int j1_end = (cj_rel == 1) ? i1+1 : i1;
-                  const int j2_begin = i2;
-                  const int j2_end = (cj_rel == 2) ? i2+1 : i2;
-
-                  for (int j0=j0_begin; j0<=j0_end; ++j0)
-                  {
-                     const int d0 = 1 + j0 - i0;
-                     for (int j1=j1_begin; j1<=j1_end; ++j1)
-                     {
-                        const int d1 = j1 - i1;
-                        for (int j2=j2_begin; j2<=j2_end; ++j2)
-                        {
-                           const int d2 = j2 - i2;
-                           int jj_lex[3];
-                           jj_lex[id0] = j0;
-                           jj_lex[id1] = j1;
-                           jj_lex[id2] = j2;
-                           const int jj_el = j_off + jj_lex[0] + jj_lex[1]*nxj + jj_lex[2]*nxj*nyj;
-                           int jj_off;
-                           if (cj_rel == 0) { jj_off = d0; }
-                           else if (cj_rel == 1) { jj_off = 3 + d0 + 2*d1; }
-                           else /* if (cj_rel == 2) */ { jj_off = 7 + d0 + 2*d2; }
-                           map(jj_off, ii_el) = jj_el;
-                        }
-                     }
-                  }
-               }
+               V(j,ix,c,iel_ho) = 0.0;
            }
         }
      }
-   }
+      MFEM_SYNC_THREAD;
+
+      MFEM_FOREACH_THREAD(xyz,x,8)
+      {
+         const int qz = xyz%2, qy = (xyz/2)%2, qx = xyz/2/2;
+         static constexpr double w = 1.0/8.0;
+
+         double local_mat_[sz_local_mat];
+         DeviceTensor<4> local_mat(local_mat_, 3,2, 3,2);
+         for (int i=0; i<sz_local_mat; ++i) { local_mat[i] = 0.0; }
+
+         double vx[8], vy[8], vz[8];
+         LORVertexCoordinates3D<1>(X, iel_ho, 0,0,0, vx, vy, vz);
+
+         double J_[3*3];
+         DeviceTensor<2> J(J_, 3,3);
+         Jacobian3D(qx,qy,qz, vx,vy,vz, J);
+         const double detJ = Det3D(J);
+         const double w_detJ = w/detJ;
+
+         const double Q0 = w_detJ*(J(0,0)*J(0,0)+J(1,0)*J(1,0)+J(2,0)*J(2,0)); // 1,1
+         const double Q1 = w_detJ*(J(0,1)*J(0,0)+J(1,1)*J(1,0)+J(2,1)*J(2,0)); // 2,1
+         const double Q2 = w_detJ*(J(0,2)*J(0,0)+J(1,2)*J(1,0)+J(2,2)*J(2,0)); // 3,1
+         const double Q3 = w_detJ*(J(0,1)*J(0,1)+J(1,1)*J(1,1)+J(2,1)*J(2,1)); // 2,2
+         const double Q4 = w_detJ*(J(0,2)*J(0,1)+J(1,2)*J(1,1)+J(2,2)*J(2,1)); // 3,2
+         const double Q5 = w_detJ*(J(0,2)*J(0,2)+J(1,2)*J(1,2)+J(2,2)*J(2,2)); // 3,3
+         const double Q6 = w_detJ;
+
+         const double mq = const_mq ? MQ(0,0,0,0) : MQ(qx,qy,qz, iel_ho);
+         const double dq = const_dq ? DQ(0,0,0,0) : DQ(qx,qy,qz, iel_ho);
+
+         // Loop over x,y,z components. 0 => x, 1 => y, 2 => z
+         for (int cj=0; cj<dim; ++cj)//MFEM_FOREACH_THREAD(cj,y,dim)
+         {
+            const int jq0 = (cj == 0) ? qx : (cj == 1) ? qy : qz;
+            const int jd0 = cj, jd1 = (cj + 1)%3, jd2 = (cj + 2)%3;
+
+            for (int bj=0; bj<2; ++bj)//MFEM_FOREACH_THREAD(bj,z,2)  // 2 faces in each dim
+            {
+               const double div_j = (bj == 0) ? -1.0 : 1.0;
+
+               double basis_j[3];
+               basis_j[jd0] = (bj == jq0) ? 1.0 : 0.0;
+               basis_j[jd1] = 0.0;
+               basis_j[jd2] = 0.0;
+
+               const int jj_loc = bj + 2*cj;
+
+               for (int ci=0; ci<dim; ++ci)
+               {
+                  const double iq0 = (ci == 0) ? qx : ((ci == 1) ? qy : qz);
+                  const int id0 = ci, id1 = (ci + 1)%3, id2 = (ci + 2)%3;
+
+                  for (int bi=0; bi<2; ++bi)
+                  {
+                     const double div_i = (bi == 0) ? -1.0 : 1.0;
+                     double basis_i[3];
+                     basis_i[id0] = (bi == iq0) ? 1.0 : 0.0;
+                     basis_i[id1] = 0.0;
+                     basis_i[id2] = 0.0;
+
+                     const int ii_loc = bi + 2*ci;
+
+                     // Only store the lower-triangular part of
+                     // the matrix (by symmetry).
+                     if (jj_loc > ii_loc) { continue; }
+
+                     const double div_div = Q6*div_i*div_j;
+
+                     double basis_basis = 0.0;
+                     basis_basis += Q0*basis_i[0]*basis_j[0];
+                     basis_basis += Q1*(basis_i[0]*basis_j[1] + basis_i[1]*basis_j[0]);
+                     basis_basis += Q2*(basis_i[0]*basis_j[2] + basis_i[2]*basis_j[0]);
+                     basis_basis += Q3*basis_i[1]*basis_j[1];
+                     basis_basis += Q4*(basis_i[1]*basis_j[2] + basis_i[2]*basis_j[1]);
+                     basis_basis += Q5*basis_i[2]*basis_j[2];
+
+                     const double val = dq*div_div + mq*basis_basis;
+                     local_mat(ci,bi, cj,bj) += val;
+                  } // bi
+               } // ci
+            } // bj
+         } // cj
+
+         // Assemble the local matrix into the macro-element sparse matrix
+         // The nonzeros of the macro-element sparse matrix are ordered as
+         // follows:
+         //
+         // The axes are ordered relative to the direction of the basis
+         // vector, e.g. for x-vectors, the axes are (x,y,z), for
+         // y-vectors the axes are (y,z,x), and for z-vectors the axes are
+         // (z,x,y).
+         //
+         // The nonzeros are then given in "rotated lexicographic"
+         // ordering, according to these axes.
+         for (int ii_loc=0; ii_loc<nf; ++ii_loc)
+         {
+            const int ci = ii_loc/2, bi = ii_loc%2;
+            const int id0 = ci, id1 = (ci+1)%3, id2 = (ci+2)%3;
+            const int i0 = bi, i1 = 0, i2 = 0;
+
+            int ii_lex[3];
+            ii_lex[id0] = i0, ii_lex[id1] = i1, ii_lex[id2] = i2;
+
+            const int nx = (ci == 0) ? op1 : o;
+            const int ny = (ci == 1) ? op1 : o;
+            const int ii = ii_lex[0] + ii_lex[1]*nx + ii_lex[2]*nx*ny;
+
+            for (int jj_loc=0; jj_loc<nf; ++jj_loc)
+            {
+               const int cj = jj_loc/2, bj = jj_loc%2;
+               // add 3 to take modulus (rather than remainder) when
+               // (cj - ci) is negative
+               const int cj_rel = (3 + cj - ci)%3;
+               const int jd0 = cj_rel, jd1 = (cj_rel+1)%3, jd2 = (cj_rel+2)%3;
+
+               int jj_rel[3];
+               jj_rel[jd0] = bj, jj_rel[jd1] = 0, jj_rel[jd2] = 0;
+
+               const int d0 = jj_rel[0] - i0 + 1;
+               const int d1 = jj_rel[1] - i1;
+               const int d2 = jj_rel[2] - i2;
+
+               const int jj_off = (cj_rel == 0) ?d0 :
+                                  (cj_rel == 1) ?  3 + d0 + 2*d1 :
+                                  (cj_rel == 2) ?  7 + d0 + 2*d2 : -1;
+
+               // Symmetry
+               const double val = (jj_loc <= ii_loc)
+                                  ? local_mat(ci,bi, cj,bj)
+                                  : local_mat(cj,bj, ci,bi);
+               AtomicAdd(V(jj_off, ii, ci, iel_ho), val);
+            } // jj_loc
+         } // ii_loc
+      } // q
+   });
+   SparseMapping3D<1>(sparse_mapping);
 }

 // Explicit template instantiations
@@ -556,7 +746,7 @@ template void BatchedLOR_RT::Assemble2D<6>();
 template void BatchedLOR_RT::Assemble2D<7>();
 template void BatchedLOR_RT::Assemble2D<8>();

-template void BatchedLOR_RT::Assemble3D<1>();
+//template void BatchedLOR_RT::Assemble3D<1>(); // explicitly specialized
 template void BatchedLOR_RT::Assemble3D<2>();
 template void BatchedLOR_RT::Assemble3D<3>();
 template void BatchedLOR_RT::Assemble3D<4>();
@@ -16,6 +16,9 @@
 #include "fem.hpp"
 #include "../general/sort_pairs.hpp"

+#define MFEM_NVTX_COLOR Crimson
+#include "../general/nvtx.hpp"
+
 namespace mfem
 {

@@ -124,6 +127,7 @@ void ParBilinearForm::pAllocMat()
 void ParBilinearForm::ParallelRAP(SparseMatrix &loc_A, OperatorHandle &A,
                                  bool steal_loc_A)
 {
+   NVTX("RAP");
   ParFiniteElementSpace &pfespace = *ParFESpace();

   // Create a block diagonal parallel matrix
@@ -2606,7 +2606,7 @@ int ParFiniteElementSpace
      if (dump < 10)
      {
         char fname[100];
-         snprintf(fname, 100, "dofs%02d.txt", MyRank);
+         sprintf(fname, "dofs%02d.txt", MyRank);
         std::ofstream f(fname);
         DebugDumpDOFs(f, deps, dof_group, dof_owner, finalized);
         dump++;
@@ -151,14 +151,14 @@ void ParGridFunction::ParallelAverage(Vector &tv) const
 {
   MFEM_VERIFY(pfes->Conforming(), "not implemented for NC meshes");
   pfes->GetProlongationMatrix()->MultTranspose(*this, tv);
-   pfes->DivideByGroupSize(tv.HostReadWrite());
+   pfes->DivideByGroupSize(tv);
 }

 void ParGridFunction::ParallelAverage(HypreParVector &tv) const
 {
   MFEM_VERIFY(pfes->Conforming(), "not implemented for NC meshes");
   pfes->GetProlongationMatrix()->MultTranspose(*this, tv);
-   pfes->DivideByGroupSize(tv.HostReadWrite());
+   pfes->DivideByGroupSize(tv);
 }

 HypreParVector *ParGridFunction::ParallelAverage() const
@@ -353,7 +353,7 @@ void ParGridFunction::GetVectorValue(int i, const IntegrationPoint &ip,
         val.SetSize(vdim);
         for (int k = 0; k < vdim; k++)
         {
-            val(k) = shape * (&loc_data[dof * k]);
+            val(k) = shape * ((const double *)loc_data + dof * k);
         }
      }
      else
@@ -468,7 +468,7 @@ void ParGridFunction::GetVectorValue(ElementTransformation &T,
      val.SetSize(vdim);
      for (int k = 0; k < vdim; k++)
      {
-         val(k) = shape * (&loc_data[dof * k]);
+         val(k) = shape * ((const double *)loc_data + dof * k);
      }
   }
   else
@@ -644,9 +644,9 @@ void ParGridFunction::ProjectBdrCoefficient(

   // Count the values globally.
   GroupCommunicator &gcomm = pfes->GroupComm();
-   gcomm.Reduce<int>(values_counter.HostReadWrite(), GroupCommunicator::Sum);
+   gcomm.Reduce<int>(values_counter, GroupCommunicator::Sum);
   // Accumulate the values globally.
-   gcomm.Reduce<double>(values.HostReadWrite(), GroupCommunicator::Sum);
+   gcomm.Reduce<double>(values, GroupCommunicator::Sum);
   // Only the values in the master are guaranteed to be correct!
   for (int i = 0; i < values.Size(); i++)
   {
@@ -682,9 +682,9 @@ void ParGridFunction::ProjectBdrCoefficientTangent(VectorCoefficient &vcoeff,

   // Count the values globally.
   GroupCommunicator &gcomm = pfes->GroupComm();
-   gcomm.Reduce<int>(values_counter.HostReadWrite(), GroupCommunicator::Sum);
+   gcomm.Reduce<int>(values_counter, GroupCommunicator::Sum);
   // Accumulate the values globally.
-   gcomm.Reduce<double>(values.HostReadWrite(), GroupCommunicator::Sum);
+   gcomm.Reduce<double>(values, GroupCommunicator::Sum);
   // Only the values in the master are guaranteed to be correct!
   for (int i = 0; i < values.Size(); i++)
   {
@@ -1109,11 +1109,11 @@ void ParGridFunction::ComputeFlux(
   SumFluxAndCount(blfi, flux, count, wcoef, subdomain);

   // Accumulate flux and counts in parallel
-   ffes->GroupComm().Reduce<double>(flux.HostReadWrite(), GroupCommunicator::Sum);
-   ffes->GroupComm().Bcast<double>(flux.HostReadWrite());
+   ffes->GroupComm().Reduce<double>(flux, GroupCommunicator::Sum);
+   ffes->GroupComm().Bcast<double>(flux);

-   ffes->GroupComm().Reduce<int>(count.HostReadWrite(), GroupCommunicator::Sum);
-   ffes->GroupComm().Bcast<int>(count.HostReadWrite());
+   ffes->GroupComm().Reduce<int>(count, GroupCommunicator::Sum);
+   ffes->GroupComm().Bcast<int>(count);

   // complete averaging
   for (int i = 0; i < count.Size(); i++)
@@ -44,6 +44,31 @@ QuadratureFunction::QuadratureFunction(Mesh *mesh, std::istream &in)
   Load(in, vdim*qspace->GetSize());
 }

+void QuadratureFunction::SetSpace(QuadratureSpaceBase *qspace_, int vdim_)
+{
+   if (qspace_ != qspace)
+   {
+      if (own_qspace) { delete qspace; }
+      qspace = qspace_;
+      own_qspace = false;
+   }
+   vdim = (vdim_ < 0) ? vdim : vdim_;
+   SetSize(vdim*qspace->GetSize());
+}
+
+void QuadratureFunction::SetSpace(
+   QuadratureSpaceBase *qspace_, double *qf_data, int vdim_)
+{
+   if (qspace_ != qspace)
+   {
+      if (own_qspace) { delete qspace; }
+      qspace = qspace_;
+      own_qspace = false;
+   }
+   vdim = (vdim_ < 0) ? vdim : vdim_;
+   NewDataAndSize(qf_data, vdim*qspace->GetSize());
+}
+
 void QuadratureFunction::Save(std::ostream &os) const
 {
   GetSpace()->Save(os);
@@ -262,34 +262,6 @@ inline void QuadratureFunction::GetValues(
   }
 }

-
-inline void QuadratureFunction::SetSpace(QuadratureSpaceBase *qspace_,
-                                         int vdim_)
-{
-   if (qspace_ != qspace)
-   {
-      if (own_qspace) { delete qspace; }
-      qspace = qspace_;
-      own_qspace = false;
-   }
-   vdim = (vdim_ < 0) ? vdim : vdim_;
-   SetSize(vdim*qspace->GetSize());
-}
-
-inline void QuadratureFunction::SetSpace(
-   QuadratureSpaceBase *qspace_, double *qf_data, int vdim_)
-{
-   if (qspace_ != qspace)
-   {
-      if (own_qspace) { delete qspace; }
-      qspace = qspace_;
-      own_qspace = false;
-   }
-   vdim = (vdim_ < 0) ? vdim : vdim_;
-   NewDataAndSize(qf_data, vdim*qspace->GetSize());
-}
-
-
 } // namespace mfem

 #endif
@@ -128,8 +128,8 @@ void TensorValues<QVectorLayout::byNODES>(const int NE,

         default:
         {
-            constexpr int MD = 8;
-            constexpr int MQ = 8;
+            constexpr int MD = 9;
+            constexpr int MQ = 9;
            MFEM_VERIFY(D1D <= MD, "Orders higher than " << MD-1
                        << " are not supported!");
            MFEM_VERIFY(Q1D <= MQ, "Quadrature rules with more than "
@@ -102,8 +102,8 @@ void TensorValues<QVectorLayout::byVDIM>(const int NE,

         default:
         {
-            constexpr int MD = 8;
-            constexpr int MQ = 8;
+            constexpr int MD = 9;
+            constexpr int MQ = 9;
            MFEM_VERIFY(D1D <= MD, "Orders higher than " << MD-1
                        << " are not supported!");
            MFEM_VERIFY(Q1D <= MQ, "Quadrature rules with more than "
@@ -110,8 +110,8 @@ void TensorDerivatives<QVectorLayout::byNODES>(const int NE,
         case 0x348: return Derivatives3D<L,P,3,4,8>(NE,B,G,J,X,Y);
         default:
         {
-            constexpr int MD = 8;
-            constexpr int MQ = 8;
+            constexpr int MD = 9;
+            constexpr int MQ = 9;
            MFEM_VERIFY(D1D <= MD, "Orders higher than " << MD-1
                        << " are not supported!");
            MFEM_VERIFY(Q1D <= MQ, "Quadrature rules with more than "
@@ -16,7 +16,7 @@ namespace mfem

 QuadratureSpaceBase::QuadratureSpaceBase(Mesh &mesh_, Geometry::Type geom,
                                         const IntegrationRule &ir)
-   : mesh(mesh_), order(ir.GetOrder())
+   : mesh(mesh_)
 {
   for (int g = 0; g < Geometry::NumGeom; g++)
   {
@@ -349,7 +349,7 @@ void StaticCondensation::ReduceRHS(const Vector &b, Vector &sc_b) const

      LUFactors lu(const_cast<double*>((const double*)A_data) + A_offsets[i],
                   const_cast<int*>((const int*)A_ipiv) + A_ipiv_offsets[i]);
-      lu.LSolve(npd, 1, b_p.GetData());
+      lu.LSolve(npd, 1, b_p);

      if (symm)
      {
@@ -527,9 +527,8 @@ void StaticCondensation::ComputeSolution(

      LUFactors lu(const_cast<double*>((const double*)A_data) + A_offsets[i],
                   const_cast<int*>((const int*)A_ipiv) + A_ipiv_offsets[i]);
-      lu.LSolve(npd, 1, b_p.GetData());
-      lu.BlockBackSolve(npd, ned, 1, lu.data + npd*npd, s_e.GetData(),
-                        b_p.GetData());
+      lu.LSolve(npd, 1, b_p);
+      lu.BlockBackSolve(npd, ned, 1, lu.data + npd*npd, s_e, b_p);

      for (int j = 0; j < npd; j++)
      {
@@ -1105,59 +1105,6 @@ public:
   virtual ~TMOP_QuadraticLimiter() { }
 };

-/// Exponential limiter function in TMOP_Integrator.
-class TMOP_ExponentialLimiter : public TMOP_LimiterFunction
-{
-public:
-   virtual double Eval(const Vector &x, const Vector &x0, double dist) const
-   {
-      MFEM_ASSERT(x.Size() == x0.Size(), "Bad input.");
-
-      return  exp(10.0*((x.DistanceSquaredTo(x0) / (dist * dist))-1.0));
-   }
-
-   virtual void Eval_d1(const Vector &x, const Vector &x0, double dist,
-                        Vector &d1) const
-   {
-      MFEM_ASSERT(x.Size() == x0.Size(), "Bad input.");
-
-      d1.SetSize(x.Size());
-      double dist_squared = dist*dist;
-      subtract(20.0*exp(10.0*((x.DistanceSquaredTo(x0) / dist_squared) - 1.0)) /
-               dist_squared, x, x0, d1);
-   }
-
-   virtual void Eval_d2(const Vector &x, const Vector &x0, double dist,
-                        DenseMatrix &d2) const
-   {
-      MFEM_ASSERT(x.Size() == x0.Size(), "Bad input.");
-      Vector tmp;
-      tmp.SetSize(x.Size());
-      double dist_squared = dist*dist;
-      double dist_squared_squared = dist_squared*dist_squared;
-      double f = exp(10.0*((x.DistanceSquaredTo(x0) / dist_squared)-1.0));
-
-      subtract(x,x0,tmp);
-      d2.SetSize(x.Size());
-      d2(0,0) = ((400.0*tmp(0)*tmp(0)*f)/dist_squared_squared)+(20.0*f/dist_squared);
-      d2(1,1) = ((400.0*tmp(1)*tmp(1)*f)/dist_squared_squared)+(20.0*f/dist_squared);
-      d2(0,1) = (400.0*tmp(0)*tmp(1)*f)/dist_squared_squared;
-      d2(1,0) = d2(0,1);
-
-      if (x.Size() == 3)
-      {
-         d2(0,2) = (400.0*tmp(0)*tmp(2)*f)/dist_squared_squared;
-         d2(1,2) = (400.0*tmp(1)*tmp(2)*f)/dist_squared_squared;
-         d2(2,0) = d2(0,2);
-         d2(2,1) = d2(1,2);
-         d2(2,2) = ((400.0*tmp(2)*tmp(2)*f)/dist_squared_squared)+(20.0*f/dist_squared);
-      }
-
-   }
-
-   virtual ~TMOP_ExponentialLimiter() { }
-};
-
 class FiniteElementCollection;
 class FiniteElementSpace;
 class ParFiniteElementSpace;
@@ -58,9 +58,8 @@ void TMOP_Integrator::AssemblePA_Limiting()

   MFEM_VERIFY(PA.enabled, "AssemblePA_Limiting but PA is not enabled!");
   MFEM_VERIFY(lim_func, "No TMOP_LimiterFunction specification!")
-   MFEM_VERIFY(dynamic_cast<TMOP_QuadraticLimiter*>(lim_func) ||
-               dynamic_cast<TMOP_ExponentialLimiter*>(lim_func),
-               "Only TMOP_QuadraticLimiter and TMOP_ExponentialLimiter are supported");
+   MFEM_VERIFY(dynamic_cast<TMOP_QuadraticLimiter*>(lim_func),
+               "Only TMOP_QuadraticLimiter is supported");

   const FiniteElementSpace *fes = PA.fes;
   const int NE = PA.ne;
@@ -25,12 +25,8 @@ MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_C0_2D,
                           const int NE,
                           const DenseTensor &j_,
                           const Array<double> &w_,
-                           const Array<double> &b_,
                           const Array<double> &bld_,
-                           const Vector &x0_,
-                           const Vector &x1_,
                           Vector &h0_,
-                           const bool exp_lim,
                           const int d1d,
                           const int q1d)
 {
@@ -45,11 +41,8 @@ MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_C0_2D,
                   Reshape(c0_.Read(), Q1D, Q1D, NE);
   const auto LD = Reshape(lim_dist.Read(), D1D, D1D, NE);
   const auto J = Reshape(j_.Read(), DIM, DIM, Q1D, Q1D, NE);
-   const auto b = Reshape(b_.Read(), Q1D, D1D);
-   const auto bld = Reshape(bld_.Read(), Q1D, D1D);
   const auto W = Reshape(w_.Read(), Q1D, Q1D);
-   const auto X0 = Reshape(x0_.Read(), D1D, D1D, DIM, NE);
-   const auto X1 = Reshape(x1_.Read(), D1D, D1D, DIM, NE);
+   const auto bld = Reshape(bld_.Read(), Q1D, D1D);

   auto H0 = Reshape(h0_.Write(), DIM, DIM, Q1D, Q1D, NE);

@@ -61,37 +54,19 @@ MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_C0_2D,
      constexpr int MQ1 = T_Q1D ? T_Q1D : T_MAX;
      constexpr int MD1 = T_D1D ? T_D1D : T_MAX;

-      MFEM_SHARED double B[MQ1*MD1];
      MFEM_SHARED double BLD[MQ1*MD1];

      MFEM_SHARED double XY[NBZ][MD1*MD1];
      MFEM_SHARED double DQ[NBZ][MD1*MQ1];
      MFEM_SHARED double QQ[NBZ][MQ1*MQ1];

-      MFEM_SHARED double XY0[2][NBZ][MD1*MD1];
-      MFEM_SHARED double DQ0[2][NBZ][MD1*MQ1];
-      MFEM_SHARED double QQ0[2][NBZ][MQ1*MQ1];
-
-      MFEM_SHARED double XY1[2][NBZ][MD1*MD1];
-      MFEM_SHARED double DQ1[2][NBZ][MD1*MQ1];
-      MFEM_SHARED double QQ1[2][NBZ][MQ1*MQ1];
-
      kernels::internal::LoadX<MD1,NBZ>(e,D1D,LD,XY);
-      kernels::internal::LoadX<MD1,NBZ>(e,D1D,X0,XY0);
-      kernels::internal::LoadX<MD1,NBZ>(e,D1D,X1,XY1);

-      kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,b,B);
      kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,bld,BLD);

      kernels::internal::EvalX<MD1,MQ1,NBZ>(D1D,Q1D,BLD,XY,DQ);
      kernels::internal::EvalY<MD1,MQ1,NBZ>(D1D,Q1D,BLD,DQ,QQ);

-      kernels::internal::EvalX<MD1,MQ1,NBZ>(D1D,Q1D,B,XY0,DQ0);
-      kernels::internal::EvalY<MD1,MQ1,NBZ>(D1D,Q1D,B,DQ0,QQ0);
-
-      kernels::internal::EvalX<MD1,MQ1,NBZ>(D1D,Q1D,B,XY1,DQ1);
-      kernels::internal::EvalY<MD1,MQ1,NBZ>(D1D,Q1D,B,DQ1,QQ1);
-
      MFEM_FOREACH_THREAD(qy,y,Q1D)
      {
         MFEM_FOREACH_THREAD(qx,x,Q1D)
@@ -102,37 +77,15 @@ MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_C0_2D,
            const double coeff0 = const_c0 ? C0(0,0,0) : C0(qx,qy,e);
            const double weight_m = weight * lim_normal * coeff0;

-            double D, p0[2], p1[2];
+            double D;
            kernels::internal::PullEval<MQ1,NBZ>(Q1D,qx,qy,QQ,D);
-            kernels::internal::PullEval<MQ1,NBZ>(Q1D,qx,qy,QQ0,p0);
-            kernels::internal::PullEval<MQ1,NBZ>(Q1D,qx,qy,QQ1,p1);
-
            const double dist = D; // GetValues, default comp set to 0

            // lim_func->Eval_d2(p1, p0, d_vals(q), grad_grad);
+            // d2.Diag(1.0 / (dist * dist), x.Size());
+            const double c = 1.0 / (dist * dist);
            double grad_grad[4];
-
-            if (!exp_lim)
-            {
-               // d2.Diag(1.0 / (dist * dist), x.Size());
-               const double c = 1.0 / (dist * dist);
-               kernels::Diag<2>(c, grad_grad);
-            }
-            else
-            {
-               double tmp[2];
-               kernels::Subtract<2>(1.0, p1, p0, tmp);
-               double dsq = kernels::DistanceSquared<2>(p1,p0);
-               double dist_squared = dist*dist;
-               double dist_squared_squared = dist_squared*dist_squared;
-               double f = exp(10.0*((dsq / dist_squared)-1.0));
-               grad_grad[0] = ((400.0*tmp[0]*tmp[0]*f)/dist_squared_squared)+
-                              (20.0*f/dist_squared);
-               grad_grad[1] = (400.0*tmp[0]*tmp[1]*f)/dist_squared_squared;
-               grad_grad[2] = grad_grad[1];
-               grad_grad[3] = ((400.0*tmp[1]*tmp[1]*f)/dist_squared_squared)+
-                              (20.0*f/dist_squared);
-            }
+            kernels::Diag<2>(c, grad_grad);
            ConstDeviceMatrix gg(grad_grad,DIM,DIM);

            for (int i = 0; i < DIM; i++)
@@ -158,17 +111,11 @@ void TMOP_Integrator::AssembleGradPA_C0_2D(const Vector &X) const
   const Vector &LD = PA.LD;
   const DenseTensor &J = PA.Jtr;
   const Array<double> &W   = PA.ir->GetWeights();
-   const Array<double> &B   = PA.maps->B;
   const Array<double> &BLD = PA.maps_lim->B;
   const Vector &C0 = PA.C0;
-   const Vector &X0 = PA.X0;
   Vector &H0 = PA.H0;

-   auto el = dynamic_cast<TMOP_ExponentialLimiter *>(lim_func);
-   const bool exp_lim = (el) ? true : false;
-
-   MFEM_LAUNCH_TMOP_KERNEL(SetupGradPA_C0_2D,id,ln,LD,C0,N,J,W,B,BLD,X0,X,H0,
-                           exp_lim);
+   MFEM_LAUNCH_TMOP_KERNEL(SetupGradPA_C0_2D,id,ln,LD,C0,N,J,W,BLD,H0);
 }

 } // namespace mfem
@@ -25,12 +25,8 @@ MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_Kernel_C0_3D,
                           const int NE,
                           const DenseTensor &j_,
                           const Array<double> &w_,
-                           const Array<double> &b_,
                           const Array<double> &bld_,
-                           const Vector &x0_,
-                           const Vector &x1_,
                           Vector &h0_,
-                           const bool exp_lim,
                           const int d1d,
                           const int q1d)
 {
@@ -44,11 +40,8 @@ MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_Kernel_C0_3D,
                   Reshape(c0_.Read(), Q1D, Q1D, Q1D, NE);
   const auto LD = Reshape(lim_dist.Read(), D1D, D1D, D1D, NE);
   const auto J = Reshape(j_.Read(), DIM, DIM, Q1D, Q1D, Q1D, NE);
-   const auto b = Reshape(b_.Read(), Q1D, D1D);
-   const auto bld = Reshape(bld_.Read(), Q1D, D1D);
   const auto W = Reshape(w_.Read(), Q1D, Q1D, Q1D);
-   const auto X0 = Reshape(x0_.Read(), D1D, D1D, D1D, DIM, NE);
-   const auto X1 = Reshape(x1_.Read(), D1D, D1D, D1D, DIM, NE);
+   const auto bld = Reshape(bld_.Read(), Q1D, D1D);

   auto H0 = Reshape(h0_.Write(), DIM, DIM, Q1D, Q1D, Q1D, NE);

@@ -61,7 +54,6 @@ MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_Kernel_C0_3D,
      constexpr int MD1 = T_D1D ? T_D1D : T_MAX;
      constexpr int MDQ = (MQ1 > MD1) ? MQ1 : MD1;

-      MFEM_SHARED double B[MQ1*MD1];
      MFEM_SHARED double sBLD[MQ1*MD1];
      kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,bld,sBLD);
      ConstDeviceMatrix BLD(sBLD, D1D, Q1D);
@@ -73,34 +65,12 @@ MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_Kernel_C0_3D,
      DeviceCube DQQ(sm0, MD1,MQ1,MQ1);
      DeviceCube QQQ(sm1, MQ1,MQ1,MQ1);

-      MFEM_SHARED double DDD0[3][MD1*MD1*MD1];
-      MFEM_SHARED double DDQ0[3][MD1*MD1*MQ1];
-      MFEM_SHARED double DQQ0[3][MD1*MQ1*MQ1];
-      MFEM_SHARED double QQQ0[3][MQ1*MQ1*MQ1];
-
-      MFEM_SHARED double DDD1[3][MD1*MD1*MD1];
-      MFEM_SHARED double DDQ1[3][MD1*MD1*MQ1];
-      MFEM_SHARED double DQQ1[3][MD1*MQ1*MQ1];
-      MFEM_SHARED double QQQ1[3][MQ1*MQ1*MQ1];
-
      kernels::internal::LoadX(e,D1D,LD,DDD);
-      kernels::internal::LoadX<MD1>(e,D1D,X0,DDD0);
-      kernels::internal::LoadX<MD1>(e,D1D,X1,DDD1);
-
-      kernels::internal::LoadB<MD1,MQ1>(D1D,Q1D,b,B);

      kernels::internal::EvalX(D1D,Q1D,BLD,DDD,DDQ);
      kernels::internal::EvalY(D1D,Q1D,BLD,DDQ,DQQ);
      kernels::internal::EvalZ(D1D,Q1D,BLD,DQQ,QQQ);

-      kernels::internal::EvalX<MD1,MQ1>(D1D,Q1D,B,DDD0,DDQ0);
-      kernels::internal::EvalY<MD1,MQ1>(D1D,Q1D,B,DDQ0,DQQ0);
-      kernels::internal::EvalZ<MD1,MQ1>(D1D,Q1D,B,DQQ0,QQQ0);
-
-      kernels::internal::EvalX<MD1,MQ1>(D1D,Q1D,B,DDD1,DDQ1);
-      kernels::internal::EvalY<MD1,MQ1>(D1D,Q1D,B,DDQ1,DQQ1);
-      kernels::internal::EvalZ<MD1,MQ1>(D1D,Q1D,B,DQQ1,QQQ1);
-
      MFEM_FOREACH_THREAD(qz,z,Q1D)
      {
         MFEM_FOREACH_THREAD(qy,y,Q1D)
@@ -113,44 +83,15 @@ MFEM_REGISTER_TMOP_KERNELS(void, SetupGradPA_Kernel_C0_3D,
               const double coeff0 = const_c0 ? C0(0,0,0,0) : C0(qx,qy,qz,e);
               const double weight_m = weight * lim_normal * coeff0;

-               double D, p0[3], p1[3];
+               double D;
               kernels::internal::PullEval(qx,qy,qz,QQQ,D);
-               kernels::internal::PullEval<MQ1>(Q1D,qx,qy,qz,QQQ0,p0);
-               kernels::internal::PullEval<MQ1>(Q1D,qx,qy,qz,QQQ1,p1);
-
               const double dist = D; // GetValues, default comp set to 0

               // lim_func->Eval_d2(p1, p0, d_vals(q), grad_grad);
-
+               // d2.Diag(1.0 / (dist * dist), x.Size());
+               const double c = 1.0 / (dist * dist);
               double grad_grad[9];
-
-               if (!exp_lim)
-               {
-                  // d2.Diag(1.0 / (dist * dist), x.Size());
-                  const double c = 1.0 / (dist * dist);
-                  kernels::Diag<3>(c, grad_grad);
-               }
-               else
-               {
-                  double tmp[3];
-                  kernels::Subtract<3>(1.0, p1, p0, tmp);
-                  double dsq = kernels::DistanceSquared<3>(p1,p0);
-                  double dist_squared = dist*dist;
-                  double dist_squared_squared = dist_squared*dist_squared;
-                  double f = exp(10.0*((dsq / dist_squared)-1.0));
-                  grad_grad[0] = ((400.0*tmp[0]*tmp[0]*f)/dist_squared_squared)+
-                                 (20.0*f/dist_squared);
-                  grad_grad[1] = (400.0*tmp[0]*tmp[1]*f)/dist_squared_squared;
-                  grad_grad[2] = (400.0*tmp[0]*tmp[2]*f)/dist_squared_squared;
-                  grad_grad[3] = grad_grad[1];
-                  grad_grad[4] = ((400.0*tmp[1]*tmp[1]*f)/dist_squared_squared)+
-                                 (20.0*f/dist_squared);
-                  grad_grad[5] = (400.0*tmp[1]*tmp[2]*f)/dist_squared_squared;
-                  grad_grad[6] = grad_grad[2];
-                  grad_grad[7] = grad_grad[5];
-                  grad_grad[8] = ((400.0*tmp[2]*tmp[2]*f)/dist_squared_squared)+
-                                 (20.0*f/dist_squared);
-               }
+               kernels::Diag<3>(c, grad_grad);
               ConstDeviceMatrix gg(grad_grad,DIM,DIM);

               for (int i = 0; i < DIM; i++)
@@ -176,17 +117,11 @@ void TMOP_Integrator::AssembleGradPA_C0_3D(const Vector &X) const
   const Vector &LD = PA.LD;
   const DenseTensor &J = PA.Jtr;
   const Array<double> &W = PA.ir->GetWeights();
-   const Array<double> &B   = PA.maps->B;
   const Array<double> &BLD = PA.maps_lim->B;
   const Vector &C0 = PA.C0;
-   const Vector &X0 = PA.X0;
   Vector &H0 = PA.H0;

-   auto el = dynamic_cast<TMOP_ExponentialLimiter *>(lim_func);
-   const bool exp_lim = (el) ? true : false;
-
-   MFEM_LAUNCH_TMOP_KERNEL(SetupGradPA_Kernel_C0_3D,id,ln,LD,C0,N,J,W,B,BLD,X0,X,
-                           H0,exp_lim);
+   MFEM_LAUNCH_TMOP_KERNEL(SetupGradPA_Kernel_C0_3D,id,ln,LD,C0,N,J,W,BLD,H0);
 }

 } // namespace mfem
@@ -30,7 +30,6 @@ MFEM_REGISTER_TMOP_KERNELS(void, AddMultPA_Kernel_C0_2D,
                           const Vector &x0_,
                           const Vector &x1_,
                           Vector &y_,
-                           const bool exp_lim,
                           const int d1d,
                           const int q1d)
 {
@@ -111,35 +110,14 @@ MFEM_REGISTER_TMOP_KERNELS(void, AddMultPA_Kernel_C0_2D,
            const double dist = ld; // GetValues, default comp set to 0

            double d1[2];
-            // Eval_d1 (Quadratic Limiter)
+            // Eval_d1
            // subtract(1.0 / (dist * dist), x, x0, d1);
            // z = a * (x - y)
            // grad = a * (x - x0)
-
-            // Eval_d1 (Exponential Limiter)
-            // double dist_squared = dist*dist;
-            // subtract(20.0*exp(10.0*((x.DistanceSquaredTo(x0) / dist_squared) - 1.0)) /
-            // dist_squared, x, x0, d1);
-            // z = a * (x - y)
-            // grad = a * (x - x0)
-
-            double a = 0.0;
+            const double a = 1.0 / (dist * dist);
            const double w = weight * lim_normal * coeff0;
-            const double dist_squared = dist * dist;
-
-            if (!exp_lim)
-            {
-               a =  1.0 / dist_squared;
-            }
-            else
-            {
-               double dsq = kernels::DistanceSquared<2>(p1,p0) / dist_squared;
-               a = 20.0*exp(10.0*(dsq - 1.0))/dist_squared;
-            }
            kernels::Subtract<2>(w*a, p1, p0, d1);
            kernels::internal::PushEval<MQ1,NBZ>(Q1D,qx,qy,d1,QQ0);
-
-
         }
      }
      MFEM_SYNC_THREAD;
@@ -165,11 +143,8 @@ void TMOP_Integrator::AddMultPA_C0_2D(const Vector &X, Vector &Y) const
   MFEM_VERIFY(PA.maps_lim->nqpt == Q1D, "");
   const Vector &X0 = PA.X0;
   const Vector &C0 = PA.C0;
-   auto el = dynamic_cast<TMOP_ExponentialLimiter *>(lim_func);
-   const bool exp_lim = (el) ? true : false;

-   MFEM_LAUNCH_TMOP_KERNEL(AddMultPA_Kernel_C0_2D,id,ln,LD,C0,N,J,W,B,BLD,X0,X,Y,
-                           exp_lim);
+   MFEM_LAUNCH_TMOP_KERNEL(AddMultPA_Kernel_C0_2D,id,ln,LD,C0,N,J,W,B,BLD,X0,X,Y);
 }

 } // namespace mfem
@@ -30,7 +30,6 @@ MFEM_REGISTER_TMOP_KERNELS(void, AddMultPA_Kernel_C0_3D,
                           const Vector &x0_,
                           const Vector &x1_,
                           Vector &y_,
-                           const bool exp_lim,
                           const int d1d,
                           const int q1d)
 {
@@ -118,32 +117,13 @@ MFEM_REGISTER_TMOP_KERNELS(void, AddMultPA_Kernel_C0_3D,
               kernels::internal::PullEval<MQ1>(Q1D,qx,qy,qz,QQQ1,p1);

               double d1[3];
-               // Eval_d1 (Quadratic Limiter)
+               // Eval_d1
               // subtract(1.0 / (dist * dist), x, x0, d1);
               // z = a * (x - y)
               // grad = a * (x - x0)
-
-               // Eval_d1 (Exponential Limiter)
-               // double dist_squared = dist*dist;
-               // subtract(20.0*exp(10.0*((x.DistanceSquaredTo(x0) / dist_squared) - 1.0)) /
-               // dist_squared, x, x0, d1);
-               // z = a * (x - y)
-               // grad = a * (x - x0)
               const double dist = D; // GetValues, default comp set to 0
-               double a = 0.0;
+               const double a = 1.0 / (dist * dist);
               const double w = weight * lim_normal * coeff0;
-               const double dist_squared = dist * dist;
-
-               if (!exp_lim)
-               {
-                  a =  1.0 / dist_squared;
-               }
-               else
-               {
-                  double dsq = kernels::DistanceSquared<3>(p1,p0) / dist_squared;
-                  a = 20.0*exp(10.0*(dsq - 1.0))/dist_squared;
-               }
-
               kernels::Subtract<3>(w*a, p1, p0, d1);
               kernels::internal::PushEval<MQ1>(Q1D,qx,qy,qz,d1,QQQ0);
            }
@@ -173,11 +153,8 @@ void TMOP_Integrator::AddMultPA_C0_3D(const Vector &X, Vector &Y) const
   MFEM_VERIFY(PA.maps_lim->nqpt == Q1D, "");
   const Vector &X0 = PA.X0;
   const Vector &C0 = PA.C0;
-   auto el = dynamic_cast<TMOP_ExponentialLimiter *>(lim_func);
-   const bool exp_lim = (el) ? true : false;

-   MFEM_LAUNCH_TMOP_KERNEL(AddMultPA_Kernel_C0_3D,id,ln,LD,C0,N,J,W,B,BLD,X0,X,Y,
-                           exp_lim);
+   MFEM_LAUNCH_TMOP_KERNEL(AddMultPA_Kernel_C0_3D,id,ln,LD,C0,N,J,W,B,BLD,X0,X,Y);
 }

 } // namespace mfem
@@ -31,7 +31,6 @@ MFEM_REGISTER_TMOP_KERNELS(double, EnergyPA_C0_2D,
                           const Vector &x1_,
                           const Vector &ones,
                           Vector &energy,
-                           const bool exp_lim,
                           const int d1d,
                           const int q1d)
 {
@@ -108,20 +107,9 @@ MFEM_REGISTER_TMOP_KERNELS(double, EnergyPA_C0_2D,
            kernels::internal::PullEval<MQ1,NBZ>(Q1D,qx,qy,QQ0,p0);
            kernels::internal::PullEval<MQ1,NBZ>(Q1D,qx,qy,QQ1,p1);
            const double dist = ld; // GetValues, default comp set to 0
-            double id2 = 0.0;
-            double dsq = 0.0;
-            if (!exp_lim)
-            {
-               id2 = 0.5 / (dist*dist);
-               dsq = kernels::DistanceSquared<2>(p1,p0) * id2;
-               E(qx,qy,e) = weight * lim_normal * dsq * coeff0;
-            }
-            else
-            {
-               id2 = 1.0 / (dist*dist);
-               dsq = kernels::DistanceSquared<2>(p1,p0) * id2;
-               E(qx,qy,e) = weight * lim_normal * exp(10.0*(dsq-1.0)) * coeff0;
-            }
+            const double id2 = 0.5 / (dist*dist);
+            const double dsq = kernels::DistanceSquared<2>(p1,p0) * id2;
+            E(qx,qy,e) = weight * lim_normal * dsq * coeff0;
         }
      }
   });
@@ -147,11 +135,7 @@ double TMOP_Integrator::GetLocalStateEnergyPA_C0_2D(const Vector &X) const
   const Vector &O = PA.O;
   Vector &E = PA.E;

-   auto el = dynamic_cast<TMOP_ExponentialLimiter *>(lim_func);
-   const bool exp_lim = (el) ? true : false;
-
-   MFEM_LAUNCH_TMOP_KERNEL(EnergyPA_C0_2D,id,ln,LD,C0,N,J,W,B,BLD,X0,X,O,E,
-                           exp_lim);
+   MFEM_LAUNCH_TMOP_KERNEL(EnergyPA_C0_2D,id,ln,LD,C0,N,J,W,B,BLD,X0,X,O,E);
 }

 } // namespace mfem
@@ -31,7 +31,6 @@ MFEM_REGISTER_TMOP_KERNELS(double, EnergyPA_C0_3D,
                           const Vector &x1_,
                           const Vector &ones,
                           Vector &energy,
-                           const bool exp_lim,
                           const int d1d,
                           const int q1d)
 {
@@ -119,20 +118,10 @@ MFEM_REGISTER_TMOP_KERNELS(double, EnergyPA_C0_3D,
               kernels::internal::PullEval<MQ1>(Q1D,qx,qy,qz,QQQ1,p1);

               const double dist = D; // GetValues, default comp set to 0
-               double id2 = 0.0;
-               double dsq = 0.0;
-               if (!exp_lim)
-               {
-                  id2 = 0.5 / (dist*dist);
-                  dsq = kernels::DistanceSquared<3>(p1,p0) * id2;
-                  E(qx,qy,qz,e) = weight * lim_normal * dsq * coeff0;
-               }
-               else
-               {
-                  id2 = 1.0 / (dist*dist);
-                  dsq = kernels::DistanceSquared<3>(p1,p0) * id2;
-                  E(qx,qy,qz,e) = weight * lim_normal * exp(10.0*(dsq-1.0)) * coeff0;
-               }
+               const double id2 = 0.5 / (dist*dist);
+
+               const double dsq = kernels::DistanceSquared<3>(p1,p0) * id2;
+               E(qx,qy,qz,e) = weight * lim_normal * dsq * coeff0;
            }
         }
      }
@@ -159,11 +148,7 @@ double TMOP_Integrator::GetLocalStateEnergyPA_C0_3D(const Vector &X) const
   const Vector &O = PA.O;
   Vector &E = PA.E;

-   auto el = dynamic_cast<TMOP_ExponentialLimiter *>(lim_func);
-   const bool exp_lim = (el) ? true : false;
-
-   MFEM_LAUNCH_TMOP_KERNEL(EnergyPA_C0_3D,id,ln,LD,C0,N,J,W,B,BLD,X0,X,O,E,
-                           exp_lim);
+   MFEM_LAUNCH_TMOP_KERNEL(EnergyPA_C0_3D,id,ln,LD,C0,N,J,W,B,BLD,X0,X,O,E);
 }

 } // namespace mfem
@@ -0,0 +1,178 @@
+// Copyright (c) 2010-2022, 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_DEBUG_HPP
+#define MFEM_DEBUG_HPP
+
+#include <string>
+#include <cstring>
+#include <iomanip>
+#include <iostream>
+
+#include "globals.hpp"
+
+#ifdef MFEM_USE_MPI
+#include <mpi.h>
+#endif
+
+namespace mfem
+{
+
+class Debug
+{
+   const bool debug = false;
+public:
+   inline Debug() {}
+
+   inline Debug(const int mpi_rank,
+                const char *FILE, const int LINE,
+                const char *FUNC, int COLOR): debug(true)
+   {
+      if (!debug) { return; }
+      const char *base = Strrnchr(FILE,'/', 2);
+      const char *file = base ? base + 1 : FILE;
+      const uint8_t color = COLOR ? COLOR : 20 + Checksum8(FILE) % 210;
+      mfem::out << "\033[38;5;" << std::to_string(color) << "m";
+      mfem::out << mpi_rank << std::setw(30) << file << ":";
+      mfem::out << "\033[2m" << std::setw(4) << LINE << "\033[22m: ";
+      if (FUNC) { mfem::out << "[" << FUNC << "] "; }
+      mfem::out << "\033[1m";
+   }
+
+   ~Debug()
+   {
+      if (!debug) { return; }
+      mfem::out << "\033[m";
+      mfem::out << std::endl;
+   }
+
+   template <typename T>
+   inline void operator<<(const T &arg) const noexcept { mfem::out << arg; }
+
+   template<typename T, typename... Args>
+   inline void operator()(const char *fmt, const T &arg,
+                          Args... args) const noexcept
+   {
+      if (!debug) { return; }
+      for (; *fmt != '\0'; fmt++ )
+      {
+         if (*fmt == '%')
+         {
+            fmt++;
+            const char c = *fmt;
+            if (c == 'p') { operator<<(arg); }
+            if (c == 's' || c == 'd' || c == 'f') { operator<<(arg); }
+            if (c == 'x' || c == 'X')
+            {
+               mfem::out << std::hex;
+               if (c == 'X') { mfem::out << std::uppercase; }
+               operator<<(arg);
+               mfem::out << std::nouppercase << std::dec;
+            }
+            if (c == '.')
+            {
+               fmt++;
+               const char c2 = *fmt;
+               char num[8] = { 0 };
+               for (int k = 0; *fmt != '\0'; fmt++, k++)
+               {
+                  if (*fmt == 'e' || *fmt == 'f') { break; }
+                  if (*fmt < 0x30 || *fmt > 0x39) { break; }
+                  num[k] = *fmt;
+               }
+               const int fx = std::atoi(num);
+               if (c2 == 'e') { mfem::out << std::scientific; }
+               if (c2 == 'f') { mfem::out << std::fixed; }
+               mfem::out << std::setprecision(fx);
+               operator<<(arg);
+               mfem::out << std::setprecision(6);
+            }
+            return operator()(fmt + 1, args...);
+         }
+         operator<<(*fmt);
+      }
+   }
+
+   template<typename T>
+   inline void operator()(const T &arg) const noexcept
+   {
+      if (!debug) { return; }
+      operator<<(arg);
+   }
+
+   inline void operator()() const noexcept { }
+
+public:
+   static const Debug Set(const char *FILE, const int LINE, const char *FUNC,
+                          int COLOR = 0)
+   {
+      static int mpi_dbg = 0, mpi_rank = 0;
+      static bool env_mpi = false, env_dbg = false;
+      static bool ini_dbg = false;
+      if (!ini_dbg)
+      {
+         const char *DBG = getenv("MFEM_DEBUG");
+         const char *MPI = getenv("MFEM_DEBUG_MPI");
+         env_dbg = DBG != nullptr;
+         env_mpi = MPI != nullptr;
+#ifdef MFEM_USE_MPI
+         int mpi_ini = false;
+         MPI_Initialized(&mpi_ini);
+         if (mpi_ini) { MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank); }
+         mpi_dbg = atoi(env_mpi ? MPI : "0");
+#endif
+         ini_dbg = true;
+      }
+      const bool debug = (env_dbg && (!env_mpi || mpi_rank == mpi_dbg));
+      return debug ? Debug(mpi_rank, FILE, LINE, FUNC, COLOR) : Debug();
+   }
+
+private:
+   inline uint8_t Checksum8(const char *bfr)
+   {
+      unsigned int chk = 0;
+      size_t len = strlen(bfr);
+      for (; len; len--,bfr++) { chk += static_cast<unsigned int>(*bfr); }
+      return (uint8_t) chk;
+   }
+
+   inline const char *Strrnchr(const char *s, const unsigned char c, int n)
+   {
+      size_t len = strlen(s);
+      char *p = const_cast<char*>(s) + len - 1;
+      for (; n; n--,p--,len--)
+      {
+         for (; len; p--,len--)
+            if (*p == c) { break; }
+         if (!len) { return nullptr; }
+         if (n == 1) { return p; }
+      }
+      return nullptr;
+   }
+
+};
+
+#ifndef MFEM_DEBUG_COLOR
+#define MFEM_DEBUG_COLOR 0
+#endif
+
+#define dbg(...) \
+    mfem::Debug::Set(__FILE__,__LINE__,__FUNCTION__,MFEM_DEBUG_COLOR).\
+    operator()(__VA_ARGS__)
+
+} // mfem namespace
+
+#define DBG(...) { printf("\033[33m");  \
+                   printf(__VA_ARGS__); \
+                   printf(" \n\033[m"); \
+                   fflush(0); }
+
+#endif // MFEM_DEBUG_HPP
@@ -61,10 +61,6 @@ void mfem_error(const char *msg = NULL);
 /// Function called by the macro MFEM_WARNING.
 void mfem_warning(const char *msg = NULL);

-#ifdef MFEM_USE_ENZYME
-static void* __enzyme_inactive_global_err = (void*)mfem_error;
-static void* __enzyme_inactive_global_warn = (void*)mfem_warning;
-#endif
 }

 #ifndef _MFEM_FUNC_NAME
@@ -310,11 +310,10 @@ inline uintptr_t MmuLengthP(const void *ptr, const size_t bytes)
 /// The protected access error, used for the host
 static void MmuError(int, siginfo_t *si, void*)
 {
-   constexpr size_t buf_size = 64;
   fflush(0);
-   char str[buf_size];
+   char str[64];
   const void *ptr = si->si_addr;
-   snprintf(str, buf_size, "Error while accessing address %p!", ptr);
+   sprintf(str, "Error while accessing address %p!", ptr);
   mfem::out << std::endl << "An illegal memory access was made!";
   MFEM_ABORT(str);
 }
@@ -908,7 +907,7 @@ void MemoryManager::SetDeviceMemoryType_(void *h_ptr, unsigned flags,
   }
 }

-void MemoryManager::Delete_(void *h_ptr, MemoryType h_mt, unsigned flags)
+MemoryType MemoryManager::Delete_(void *h_ptr, MemoryType h_mt, unsigned flags)
 {
   const bool alias = flags & Mem::ALIAS;
   const bool registered = flags & Mem::REGISTERED;
@@ -925,7 +924,7 @@ void MemoryManager::Delete_(void *h_ptr, MemoryType h_mt, unsigned flags)
   MFEM_ASSERT(registered || !(owns_host || owns_device || owns_internal) ||
               (!(owns_device || owns_internal) && h_ptr == nullptr),
               "invalid Memory state");
-   if (!mm.exists || !registered) { return; }
+   if (!mm.exists || !registered) { return h_mt; }
   if (alias)
   {
      if (owns_internal)
@@ -946,6 +945,7 @@ void MemoryManager::Delete_(void *h_ptr, MemoryType h_mt, unsigned flags)
         mm.Erase(h_ptr, owns_device);
      }
   }
+   return h_mt;
 }

 void MemoryManager::DeleteDevice_(void *h_ptr, unsigned & flags)
@@ -644,8 +644,9 @@ private: // Static methods used by the Memory<T> class
   static void SetDeviceMemoryType_(void *h_ptr, unsigned flags,
                                    MemoryType d_mt);

-   /// Un-register and free memory identified by its host pointer.
-   static void Delete_(void *h_ptr, MemoryType mt, unsigned flags);
+   /// Un-register and free memory identified by its host pointer. Returns the
+   /// memory type of the host pointer.
+   static MemoryType Delete_(void *h_ptr, MemoryType mt, unsigned flags);

   /// Free device memory identified by its host pointer
   static void DeleteDevice_(void *h_ptr, unsigned & flags);
@@ -827,23 +828,6 @@ public:

   static MemoryType GetHostMemoryType() { return host_mem_type; }
   static MemoryType GetDeviceMemoryType() { return device_mem_type; }
-
-#ifdef MFEM_USE_ENZYME
-   static void myfree(void* mem, MemoryType MT, unsigned &flags)
-   {
-      MemoryManager::Delete_(mem, MT, flags);
-   }
-   __attribute__((used))
-   inline static void* __enzyme_allocation_like1[4] = {(void*)static_cast<void*(*)(void*, size_t, MemoryType, unsigned&)>(MemoryManager::New_),
-                                                       (void*)1, (void*)"-1,2,3", (void*)myfree
-                                                      };
-   __attribute__((used))
-   inline static void* __enzyme_allocation_like2[4] = {(void*)static_cast<void*(*)(void*, size_t, MemoryType, MemoryType, unsigned, unsigned&)>(MemoryManager::New_),
-                                                       (void*)1, (void*)"-1,2,4", (void*)MemoryManager::Delete_
-                                                      };
-   __attribute__((used))
-   inline static void* __enzyme_function_like[2] = {(void*)MemoryManager::Delete_, (void*)"free"};
-#endif
 };


@@ -1023,12 +1007,8 @@ inline void Memory<T>::Delete()
   const bool mt_host = h_mt == MemoryType::HOST;
   const bool std_delete = !registered && mt_host;

-   if (!std_delete)
-   {
-      MemoryManager::Delete_((void*)h_ptr, h_mt, flags);
-   }
-
-   if (mt_host)
+   if (std_delete ||
+       MemoryManager::Delete_((void*)h_ptr, h_mt, flags) == MemoryType::HOST)
   {
      if (flags & OWNS_HOST) { delete [] h_ptr; }
   }
@@ -0,0 +1,312 @@
+// Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at
+// the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights
+// reserved. See files LICENSE and NOTICE for details.
+//
+// This file is part of CEED, a collection of benchmarks, miniapps, software
+// libraries and APIs for efficient high-order finite element and spectral
+// element discretizations for exascale applications. For more information and
+// source code availability see http://github.com/ceed.
+//
+// The CEED research is supported by the Exascale Computing Project 17-SC-20-SC,
+// a collaborative effort of two U.S. Department of Energy organizations (Office
+// of Science and the National Nuclear Security Administration) responsible for
+// the planning and preparation of a capable exascale ecosystem, including
+// software, applications, hardware, advanced system engineering and early
+// testbed platforms, in support of the nation's exascale computing imperative.
+
+#ifndef MFEM_NVVP_HPP
+#define MFEM_NVVP_HPP
+
+#include <string>
+#include <cassert>
+
+#include "globals.hpp"
+#include "backends.hpp"
+
+#if defined(MFEM_USE_CUDA) //||1
+#include <cuda.h>
+#include <nvToolsExt.h>
+#include <cudaProfiler.h>
+#include <cuda_runtime.h>
+#else // MFEM_USE_CUDA
+typedef struct
+{
+   uint16_t version;
+   uint16_t size;
+   uint32_t category;
+   int32_t colorType;
+   uint32_t color;
+   int32_t payloadType;
+   int32_t reserved0;
+   int32_t messageType;
+   struct { const char* ascii; } message;
+} nvtxEventAttributes_t;
+#define NVTX_VERSION 2
+#define NVTX_COLOR_ARGB 1
+#define NVTX_MESSAGE_TYPE_ASCII 1
+#define NVTX_EVENT_ATTRIB_STRUCT_SIZE \
+    ((uint16_t)(sizeof(nvtxEventAttributes_t)))
+inline int nvtxRangePushEx(const nvtxEventAttributes_t*) { return 0; }
+inline int nvtxRangePop(void) { return 0; }
+#endif // MFEM_USE_CUDA
+
+namespace mfem
+{
+
+// https://en.wikipedia.org/wiki/Web_colors#Extended_colors
+typedef enum
+{
+   Black, NavyBlue, DarkBlue, MediumBlue, Blue, DarkGreen, WebGreen, Teal,
+   DarkCyan, DeepSkyBlue, DarkTurquoise, MediumSpringGreen, Green, Lime,
+   SpringGreen, Aqua, Cyan, MidnightBlue, DodgerBlue, LightSeaGreen,
+   ForestGreen, SeaGreen, DarkSlateGray, LimeGreen, MediumSeaGreen,
+   Turquoise, RoyalBlue, SteelBlue, DarkSlateBlue, MediumTurquoise, Indigo,
+   DarkOliveGreen, CadetBlue, Cornflower, RebeccaPurple, MediumAquamarine,
+   DimGray, SlateBlue, OliveDrab, SlateGray, LightSlateGray,
+   MediumSlateBlue, LawnGreen, WebMaroon, WebPurple, Chartreuse,
+   Aquamarine, Olive, WebGray, SkyBlue, LightSkyBlue, BlueViolet, DarkRed,
+   DarkMagenta, SaddleBrown, DarkSeaGreen, LightGreen, MediumPurple,
+   DarkViolet, PaleGreen, DarkOrchid, YellowGreen, Purple, Sienna, Brown,
+   DarkGray, LightBlue, GreenYellow, PaleTurquoise, Maroon,
+   LightSteelBlue, PowderBlue, Firebrick, DarkGoldenrod, MediumOrchid,
+   RosyBrown, DarkKhaki, Gray, Silver, MediumVioletRed, IndianRed, Peru,
+   Chocolate, Tan, LightGray, Thistle, Orchid, Goldenrod, PaleVioletRed,
+   Crimson, Gainsboro, Plum, Burlywood, LightCyan, Lavender, DarkSalmon,
+   Violet, PaleGoldenrod, LightCoral, Khaki, AliceBlue, Honeydew, Azure,
+   SandyBrown, Wheat, Beige, WhiteSmoke, MintCream, GhostWhite, Salmon,
+   AntiqueWhite, Linen, LightGoldenrod, OldLace, Red, Fuchsia, Magenta,
+   DeepPink, OrangeRed, Tomato, HotPink, Coral, DarkOrange, LightSalmon,
+   Orange, LightPink, Pink, Gold, PeachPuff, NavajoWhite, Moccasin,
+   Bisque, MistyRose, BlanchedAlmond, PapayaWhip, LavenderBlush, Seashell,
+   Cornsilk, LemonChiffon, FloralWhite, Snow, Yellow, LightYellow, Ivory,
+   White
+} COLOR_NAMES;
+
+static constexpr uint32_t HEX_COLORS[] =
+{
+   0x000000, 0x000080, 0x00008B, 0x0000CD, 0x0000FF, 0x006400, 0x008000,
+   0x008080, 0x008B8B, 0x00BFFF, 0x00CED1, 0x00FA9A, 0x00FF00, 0x00FF00,
+   0x00FF7F, 0x00FFFF, 0x00FFFF, 0x191970, 0x1E90FF, 0x20B2AA, 0x228B22,
+   0x2E8B57, 0x2F4F4F, 0x32CD32, 0x3CB371, 0x40E0D0, 0x4169E1, 0x4682B4,
+   0x483D8B, 0x48D1CC, 0x4B0082, 0x556B2F, 0x5F9EA0, 0x6495ED, 0x663399,
+   0x66CDAA, 0x696969, 0x6A5ACD, 0x6B8E23, 0x708090, 0x778899, 0x7B68EE,
+   0x7CFC00, 0x7F0000, 0x7F007F, 0x7FFF00, 0x7FFFD4, 0x808000, 0x808080,
+   0x87CEEB, 0x87CEFA, 0x8A2BE2, 0x8B0000, 0x8B008B, 0x8B4513, 0x8FBC8F,
+   0x90EE90, 0x9370DB, 0x9400D3, 0x98FB98, 0x9932CC, 0x9ACD32, 0xA020F0,
+   0xA0522D, 0xA52A2A, 0xA9A9A9, 0xADD8E6, 0xADFF2F, 0xAFEEEE, 0xB03060,
+   0xB0C4DE, 0xB0E0E6, 0xB22222, 0xB8860B, 0xBA55D3, 0xBC8F8F, 0xBDB76B,
+   0xBEBEBE, 0xC0C0C0, 0xC71585, 0xCD5C5C, 0xCD853F, 0xD2691E, 0xD2B48C,
+   0xD3D3D3, 0xD8BFD8, 0xDA70D6, 0xDAA520, 0xDB7093, 0xDC143C, 0xDCDCDC,
+   0xDDA0DD, 0xDEB887, 0xE0FFFF, 0xE6E6FA, 0xE9967A, 0xEE82EE, 0xEEE8AA,
+   0xF08080, 0xF0E68C, 0xF0F8FF, 0xF0FFF0, 0xF0FFFF, 0xF4A460, 0xF5DEB3,
+   0xF5F5DC, 0xF5F5F5, 0xF5FFFA, 0xF8F8FF, 0xFA8072, 0xFAEBD7, 0xFAF0E6,
+   0xFAFAD2, 0xFDF5E6, 0xFF0000, 0xFF00FF, 0xFF00FF, 0xFF1493, 0xFF4500,
+   0xFF6347, 0xFF69B4, 0xFF7F50, 0xFF8C00, 0xFFA07A, 0xFFA500, 0xFFB6C1,
+   0xFFC0CB, 0xFFD700, 0xFFDAB9, 0xFFDEAD, 0xFFE4B5, 0xFFE4C4, 0xFFE4E1,
+   0xFFEBCD, 0xFFEFD5, 0xFFF0F5, 0xFFF5EE, 0xFFF8DC, 0xFFFACD, 0xFFFAF0,
+   0xFFFAFA, 0xFFFF00, 0xFFFFE0, 0xFFFFF0, 0xFFFFFF
+};
+
+static const int NUM_HEX_COLORS = sizeof(HEX_COLORS)/sizeof(uint32_t);
+
+class Nvtx
+{
+   const bool nvtx = false;
+   const bool enforce_kernel_sync = false;
+   const char *base, *file;
+   const uint32_t color = Black;
+   mutable std::string ascii;
+   mutable nvtxEventAttributes_t event;
+
+public:
+   Nvtx() { }
+
+   Nvtx(bool enforce_kernel_sync,
+        const char *FILE, const int LINE, const char *FUNC, uint32_t COLOR):
+      nvtx(true),
+      enforce_kernel_sync(enforce_kernel_sync),
+      base(Strrnchr(FILE,'/', 2)),
+      file(base ? base + 1 : FILE),
+      color(COLOR),
+      ascii(file),
+      event({})
+   {
+      event.version = NVTX_VERSION;
+      event.size = NVTX_EVENT_ATTRIB_STRUCT_SIZE;
+      event.colorType = NVTX_COLOR_ARGB;
+      event.color = HEX_COLORS[color % NUM_HEX_COLORS];
+      event.messageType = NVTX_MESSAGE_TYPE_ASCII;
+
+      ascii += ":";
+      ascii += std::to_string(LINE);
+      ascii += ":[";
+      ascii += FUNC;
+      ascii += "] ";
+   }
+
+   Nvtx(const char *title, int color = Wheat, bool enforce_kernel_sync = true):
+      nvtx(true),
+      enforce_kernel_sync(enforce_kernel_sync),
+      color(color),
+      ascii(title),
+      event({})
+   {
+      event.version = NVTX_VERSION;
+      event.size = NVTX_EVENT_ATTRIB_STRUCT_SIZE;
+      event.colorType = NVTX_COLOR_ARGB;
+      event.color = HEX_COLORS[color % NUM_HEX_COLORS];
+      event.messageType = NVTX_MESSAGE_TYPE_ASCII;
+      event.message.ascii = ascii.c_str();
+      nvtxRangePushEx(&event); // push
+   }
+
+   ~Nvtx()
+   {
+      if (!nvtx) { return; }
+      if (enforce_kernel_sync)
+      {
+         nvtxEventAttributes_t eks = {};
+         eks.version = NVTX_VERSION;
+         eks.size = NVTX_EVENT_ATTRIB_STRUCT_SIZE;
+         eks.category = 0; // user value
+         eks.colorType = NVTX_COLOR_ARGB;
+         eks.messageType = NVTX_MESSAGE_TYPE_ASCII;
+         eks.message.ascii = "Sync";
+         eks.color = HEX_COLORS[Yellow];
+         nvtxRangePushEx(&eks);
+         MFEM_STREAM_SYNC;
+         nvtxRangePop();
+      }
+      nvtxRangePop(); // pop
+   }
+
+   // used through MFEM_NVTX
+   inline void operator()() const noexcept
+   {
+      event.message.ascii = ascii.c_str();
+      nvtxRangePushEx(&event); // push
+   }
+
+   template<typename T>
+   inline void operator()(const T &arg) const noexcept
+   {
+      if (!nvtx) { return; }
+      operator<<(arg);
+      event.message.ascii = ascii.c_str();
+      nvtxRangePushEx(&event); // push
+   }
+
+   template<typename T, typename... Args>
+   inline void operator()(const char *fmt, const T &arg,
+                          Args... args) const noexcept
+   {
+      if (!nvtx) { return; }
+
+      for (; *fmt != '\0'; fmt++ )
+      {
+         if (*fmt == '%')
+         {
+            fmt++;
+            const char c = *fmt;
+            if (c == 'p') { operator<<(arg); }
+            if (c == 's' || c == 'd' || c == 'f') { operator<<(arg); }
+            if (c == 'x' || c == 'X')
+            {
+               mfem::out << std::hex;
+               if (c == 'X') { mfem::out << std::uppercase; }
+               operator<<(arg);
+               mfem::out << std::nouppercase << std::dec;
+            }
+            if (c == '.')
+            {
+               fmt++;
+               const char c2 = *fmt;
+               char num[8] = { 0 };
+               for (int k = 0; *fmt != '\0'; fmt++, k++)
+               {
+                  if (*fmt == 'e' || *fmt == 'f') { break; }
+                  if (*fmt < 0x30 || *fmt > 0x39) { break; }
+                  num[k] = *fmt;
+               }
+               const int fx = std::atoi(num);
+               if (c2 == 'e') { mfem::out << std::scientific; }
+               if (c2 == 'f') { mfem::out << std::fixed; }
+               mfem::out << std::setprecision(fx);
+               operator<<(arg);
+               mfem::out << std::setprecision(6);
+            }
+            return operator()(fmt + 1, args...);
+         }
+         operator<<(*fmt);
+      }
+      // should never be here
+      assert(false);
+   }
+
+   template <typename T>
+   inline void operator<<(const T &arg) const noexcept
+   {
+      if (!nvtx) { return; }
+      ascii += arg;
+   }
+
+   inline void operator<<(const int &arg) const noexcept
+   {
+      if (!nvtx) { return; }
+      ascii += std::to_string(arg);
+   }
+
+
+public:
+   static const Nvtx Set(const char *FILE, const int LINE, const char *FUNC,
+                         uint32_t COLOR)
+   {
+      static bool env_nvtx = false;
+      static bool env_eks = false;
+      static bool ini_nvtx = false;
+      if (!ini_nvtx)
+      {
+         env_nvtx = getenv("MFEM_NVTX") != nullptr;
+         env_eks = getenv("MFEM_EKS") != nullptr;
+         ini_nvtx = true;
+      }
+      return env_nvtx ? Nvtx(env_eks, FILE, LINE, FUNC, COLOR) : Nvtx();
+   }
+
+private:
+   inline const char *Strrnchr(const char *s, const unsigned char c, int n)
+   {
+      size_t len = strlen(s);
+      char *p = const_cast<char*>(s) + len - 1;
+      for (; n; n--,p--,len--)
+      {
+         for (; len; p--,len--)
+            if (*p == c) { break; }
+         if (!len) { return nullptr; }
+         if (n == 1) { return p; }
+      }
+      return nullptr;
+   }
+};
+
+#ifndef MFEM_NVTX_COLOR
+#define MFEM_NVTX_COLOR SeaGreen
+#endif
+
+// Helpers for generating unique variable names
+#define NVTX_PRIVATE_NAME(name) NVTX_PRIVATE_CONCAT(name, __LINE__)
+#define NVTX_PRIVATE_CONCAT(a, b) NVTX_PRIVATE_CONCAT2(a, b)
+#define NVTX_PRIVATE_CONCAT2(a, b) a##b
+
+// temporary object which is only alive in the expression
+// __PRETTY_FUNCTION__
+#define NVTX(...) \
+    mfem::Nvtx NVTX_PRIVATE_NAME(nvtx) = \
+    mfem::Nvtx::Set(__FILE__,__LINE__,__FUNCTION__,MFEM_NVTX_COLOR);\
+    NVTX_PRIVATE_NAME(nvtx).operator()(__VA_ARGS__)
+
+#define MFEM_NVTX NVTX()
+
+} // namespace mfem
+
+#endif // MFEM_NVVP_HPP
@@ -112,6 +112,9 @@ const char *GetConfigStr()
 #ifdef MFEM_USE_MEMALLOC
      "MFEM_USE_MEMALLOC\n"
 #endif
+#ifdef MFEM_USE_MESQUITE
+      "MFEM_USE_MESQUITE\n"
+#endif
 #ifdef MFEM_USE_METIS
      "MFEM_USE_METIS\n"
 #endif
@@ -60,9 +60,9 @@ namespace strict_fstream
 {

 // Overloaded error checks to handle POSIX and GNU strerror_r
-inline char* check_strerror_r(int r, char* buff, size_t buff_size)
+inline char* check_strerror_r(int r, char* buff, int err)
 {
-   if (r) { snprintf(buff, buff_size, "unknown error: %d", r); }
+   if (r) { sprintf(buff, "unknown error: %d", err); }
   return buff;
 }

@@ -51,14 +51,14 @@ Eliminator::Eliminator(const SparseMatrix& B, const Array<int>& lagrange_tdofs_,
 void Eliminator::Eliminate(const Vector& vin, Vector& vout) const
 {
   Bp.Mult(vin, vout);
-   Bsinverse.Solve(Bs.Height(), 1, vout.GetData());
+   Bsinverse.Solve(Bs.Height(), 1, vout);
   vout *= -1.0;
 }

 void Eliminator::EliminateTranspose(const Vector& vin, Vector& vout) const
 {
   Vector work(vin);
-   BsTinverse.Solve(Bs.Height(), 1, work.GetData());
+   BsTinverse.Solve(Bs.Height(), 1, work);
   Bp.MultTranspose(work, vout);
   vout *= -1.0;
 }
@@ -66,14 +66,14 @@ void Eliminator::EliminateTranspose(const Vector& vin, Vector& vout) const
 void Eliminator::LagrangeSecondary(const Vector& vin, Vector& vout) const
 {
   vout = vin;
-   Bsinverse.Solve(Bs.Height(), 1, vout.GetData());
+   Bsinverse.Solve(Bs.Height(), 1, vout);
 }

 void Eliminator::LagrangeSecondaryTranspose(const Vector& vin,
                                            Vector& vout) const
 {
   vout = vin;
-   BsTinverse.Solve(Bs.Height(), 1, vout.GetData());
+   BsTinverse.Solve(Bs.Height(), 1, vout);
 }

 void Eliminator::ExplicitAssembly(DenseMatrix& mat) const
@@ -55,9 +55,7 @@ void CPardisoSolver::SetOperator(const Operator &op)
   auto parcsr_op = static_cast<hypre_ParCSRMatrix *>(
                       const_cast<HypreParMatrix &>(hypreParMat));

-   hypreParMat.HostRead();
   hypre_CSRMatrix *csr_op = hypre_MergeDiagAndOffd(parcsr_op);
-   hypreParMat.HypreRead();
 #if MFEM_HYPRE_VERSION >= 21600
   hypre_CSRMatrixBigJtoJ(csr_op);
 #endif
@@ -175,26 +175,12 @@ void DenseMatrix::Mult(const double *x, double *y) const
   kernels::Mult(height, width, Data(), x, y);
 }

-void DenseMatrix::Mult(const double *x, Vector &y) const
-{
-   MFEM_ASSERT(height == y.Size(), "incompatible dimensions");
-
-   Mult(x, y.GetData());
-}
-
-void DenseMatrix::Mult(const Vector &x, double *y) const
-{
-   MFEM_ASSERT(width == x.Size(), "incompatible dimensions");
-
-   Mult(x.GetData(), y);
-}
-
 void DenseMatrix::Mult(const Vector &x, Vector &y) const
 {
   MFEM_ASSERT(height == y.Size() && width == x.Size(),
               "incompatible dimensions");

-   Mult(x.GetData(), y.GetData());
+   Mult((const double *)x, (double *)y);
 }

 double DenseMatrix::operator *(const DenseMatrix &m) const
@@ -227,26 +213,12 @@ void DenseMatrix::MultTranspose(const double *x, double *y) const
   }
 }

-void DenseMatrix::MultTranspose(const double *x, Vector &y) const
-{
-   MFEM_ASSERT(width == y.Size(), "incompatible dimensions");
-
-   MultTranspose(x, y.GetData());
-}
-
-void DenseMatrix::MultTranspose(const Vector &x, double *y) const
-{
-   MFEM_ASSERT(height == x.Size(), "incompatible dimensions");
-
-   MultTranspose(x.GetData(), y);
-}
-
 void DenseMatrix::MultTranspose(const Vector &x, Vector &y) const
 {
   MFEM_ASSERT(height == x.Size() && width == y.Size(),
               "incompatible dimensions");

-   MultTranspose(x.GetData(), y.GetData());
+   MultTranspose((const double *)x, (double *)y);
 }

 void DenseMatrix::AddMult(const Vector &x, Vector &y) const
@@ -254,8 +226,8 @@ void DenseMatrix::AddMult(const Vector &x, Vector &y) const
   MFEM_ASSERT(height == y.Size() && width == x.Size(),
               "incompatible dimensions");

-   const double *xp = x.GetData(), *d_col = data;
-   double *yp = y.GetData();
+   const double *xp = x, *d_col = data;
+   double *yp = y;
   for (int col = 0; col < width; col++)
   {
      double x_col = xp[col];
@@ -290,8 +262,8 @@ void DenseMatrix::AddMult_a(double a, const Vector &x, Vector &y) const
   MFEM_ASSERT(height == y.Size() && width == x.Size(),
               "incompatible dimensions");

-   const double *xp = x.GetData(), *d_col = data;
-   double *yp = y.GetData();
+   const double *xp = x, *d_col = data;
+   double *yp = y;
   for (int col = 0; col < width; col++)
   {
      const double x_col = a*xp[col];
@@ -1209,7 +1181,7 @@ void DenseMatrix::SingularValues(Vector &sv) const
   int         n            = Width();
   double      *a           = copy_of_this.data;
   sv.SetSize(min(m, n));
-   double      *s           = sv.GetData();
+   double      *s           = sv;
   double      *u           = NULL;
   double      *vt          = NULL;
   double      *work        = NULL;
@@ -1487,10 +1459,10 @@ void DenseMatrix::GradToCurl(DenseMatrix &curl)
         int j = i+n;

         // curl of (Ui,0)
-         curl(i,0) = -y;
+         curl(i,0) = y;

         // curl of (0,Ui)
-         curl(j,0) = x;
+         curl(j,0) = -x;
      }
   }
   else
@@ -1523,20 +1495,6 @@ void DenseMatrix::GradToCurl(DenseMatrix &curl)
   }
 }

-void DenseMatrix::GradToVectorCurl2D(DenseMatrix &curl)
-{
-   MFEM_VERIFY(Width() == 2,
-               "DenseMatrix::GradToVectorCurl2D(...): dimension must be 2")
-
-   int n = Height();
-   // rotate gradient
-   for (int i = 0; i < n; i++)
-   {
-      curl(i,0) = (*this)(i,1);
-      curl(i,1) = -(*this)(i,0);
-   }
-}
-
 void DenseMatrix::GradToDiv(Vector &div)
 {
   MFEM_ASSERT(Width()*Height() == div.Size(), "incompatible Vector 'div'!");
@@ -4218,10 +4176,10 @@ const
 {
   int n = SizeI(), ne = SizeK();
   const int *I = elem_dof.GetI(), *J = elem_dof.GetJ(), *dofs;
-   const double *d_col = mfem::HostRead(tdata, n*SizeJ()*ne);
+   const double *d_col = tdata;
   double *yp = y.HostReadWrite();
   double x_col;
-   const double *xp = x.HostRead();
+   const double *xp = x;
   // the '4' here can be tuned for given platform and compiler
   if (n <= 4)
   {
@@ -141,24 +141,12 @@ public:
   /// Matrix vector multiplication.
   void Mult(const double *x, double *y) const;

-   /// Matrix vector multiplication.
-   void Mult(const double *x, Vector &y) const;
-
-   /// Matrix vector multiplication.
-   void Mult(const Vector &x, double *y) const;
-
   /// Matrix vector multiplication.
   virtual void Mult(const Vector &x, Vector &y) const;

   /// Multiply a vector with the transpose matrix.
   void MultTranspose(const double *x, double *y) const;

-   /// Multiply a vector with the transpose matrix.
-   void MultTranspose(const double *x, Vector &y) const;
-
-   /// Multiply a vector with the transpose matrix.
-   void MultTranspose(const Vector &x, double *y) const;
-
   /// Multiply a vector with the transpose matrix.
   virtual void MultTranspose(const Vector &x, Vector &y) const;

@@ -192,7 +180,7 @@ public:

   /// Compute y^t A x
   double InnerProduct(const Vector &x, const Vector &y) const
-   { return InnerProduct(x.GetData(), y.GetData()); }
+   { return InnerProduct((const double *)x, (const double *)y); }

   /// Returns a pointer to the inverse matrix.
   virtual MatrixInverse *Inverse() const;
@@ -248,13 +236,6 @@ public:
   /// Take the 2-norm of the columns of A and store in v
   void Norm2(double *v) const;

-   /// Take the 2-norm of the columns of A and store in v
-   void Norm2(Vector &v) const
-   {
-      MFEM_ASSERT(v.Size() == Width(), "incompatible Vector size!");
-      Norm2(v.GetData());
-   }
-
   /// Compute the norm ||A|| = max_{ij} |A_{ij}|
   double MaxMaxNorm() const;

@@ -344,13 +325,8 @@ public:
   /** Given a DShape matrix (from a scalar FE), stored in *this, returns the
       CurlShape matrix. If *this is a N by D matrix, then curl is a D*N by
       D*(D-1)/2 matrix. The size of curl must be set outside. The dimension D
-       can be either 2 or 3. In 2D this computes the scalar-valued curl of a
-       2D vector field */
+       can be either 2 or 3. */
   void GradToCurl(DenseMatrix &curl);
-   /** Given a DShape matrix (from a scalar FE), stored in *this, returns the
-       CurlShape matrix. This computes the vector-valued curl of a scalar field.
-       *this is N by 2 matrix and curl is N by 2 matrix as well. */
-   void GradToVectorCurl2D(DenseMatrix &curl);
   /** Given a DShape matrix (from a scalar FE), stored in *this,
       returns the DivShape vector. If *this is a N by dim matrix,
       then div is a dim*N vector. The size of div must be set
@@ -92,8 +92,6 @@ public:
    *  internally in HiOp. */
   virtual bool get_starting_point(const hiop::size_type &n, double *x0);

-   using hiop::hiopInterfaceBase::get_starting_point;
-
   virtual bool get_vars_info(const hiop::size_type &n, double *xlow, double* xupp,
                              NonlinearityType* type);

@@ -138,8 +136,6 @@ public:
                          const hiop::index_type *idx_cons,
                          const double *x, bool new_x, double *cons);

-   using hiop::hiopInterfaceBase::eval_cons;
-
   /** Evaluates the Jacobian of the subset of constraints indicated by
    *  idx_cons. The idx_cons is assumed to be of size num_cons.
    *  Example: if cons[c] = C(x)[idx_cons[c]] where c = 0 .. num_cons-1, then
@@ -156,8 +152,6 @@ public:
                              const hiop::index_type *idx_cons,
                              const double *x, bool new_x, double *Jac);

-   using hiop::hiopInterfaceDenseConstraints::eval_Jac_cons;
-
   /** Specifies column partitioning for distributed memory vectors.
    *  Process p owns vector entries with indices cols[p] to cols[p+1]-1,
    *  where p = 0 .. nranks-1. The cols array is of size nranks + 1.
@@ -22,6 +22,8 @@
 #include <cmath>
 #include <cstdlib>

+#include "../general/nvtx.hpp"
+
 using namespace std;

 namespace mfem
@@ -1732,6 +1734,10 @@ void HypreParMatrix::ResetTranspose() const
 HYPRE_Int HypreParMatrix::Mult(HypreParVector &x, HypreParVector &y,
                               double a, double b) const
 {
+#undef MFEM_NVTX_COLOR
+#define MFEM_NVTX_COLOR Crimson
+   NVTX("SpMV");
+
   x.HypreRead();
   (b == 0.0) ? y.HypreWrite() : y.HypreReadWrite();
   return hypre_ParCSRMatrixMatvec(a, A, x, b, y);
@@ -1739,6 +1745,10 @@ HYPRE_Int HypreParMatrix::Mult(HypreParVector &x, HypreParVector &y,

 void HypreParMatrix::Mult(double a, const Vector &x, double b, Vector &y) const
 {
+#undef MFEM_NVTX_COLOR
+#define MFEM_NVTX_COLOR Crimson
+   NVTX("SpMV");
+
   MFEM_ASSERT(x.Size() == Width(), "invalid x.Size() = " << x.Size()
               << ", expected size = " << Width());
   MFEM_ASSERT(y.Size() == Height(), "invalid y.Size() = " << y.Size()
@@ -1797,6 +1807,10 @@ void HypreParMatrix::Mult(double a, const Vector &x, double b, Vector &y) const
 void HypreParMatrix::MultTranspose(double a, const Vector &x,
                                   double b, Vector &y) const
 {
+#undef MFEM_NVTX_COLOR
+#define MFEM_NVTX_COLOR Crimson
+   NVTX("SpMV_T");
+
   MFEM_ASSERT(x.Size() == Height(), "invalid x.Size() = " << x.Size()
               << ", expected size = " << Height());
   MFEM_ASSERT(y.Size() == Width(), "invalid y.Size() = " << y.Size()
@@ -3880,6 +3894,10 @@ void HypreSolver::Setup(const HypreParVector &b, HypreParVector &x) const
 {
   if (setup_called) { return; }

+#undef MFEM_NVTX_COLOR
+#define MFEM_NVTX_COLOR Indigo
+   NVTX("AMG Setup");
+
   MFEM_VERIFY(A != NULL, "HypreParMatrix A is missing");

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

 void HypreSolver::Mult(const HypreParVector &b, HypreParVector &x) const
 {
+#undef MFEM_NVTX_COLOR
+#define MFEM_NVTX_COLOR MediumSpringGreen
+   NVTX("AMG V-cycle");
+
   HYPRE_Int err_flag;
   if (A == NULL)
   {
@@ -42,9 +42,7 @@ void MUMPSSolver::SetOperator(const Operator &op)

   auto parcsr_op = (hypre_ParCSRMatrix *) const_cast<HypreParMatrix &>(*APtr);

-   APtr->HostRead();
   hypre_CSRMatrix *csr_op = hypre_MergeDiagAndOffd(parcsr_op);
-   APtr->HypreRead();
 #if MFEM_HYPRE_VERSION >= 21600
   hypre_CSRMatrixBigJtoJ(csr_op);
 #endif
@@ -1335,12 +1335,7 @@ void PetscParMatrix::ConvertOperator(MPI_Comm comm, const Operator &op, Mat* A,
      PetscBool ismatis;
 #endif

-#if PETSC_VERSION_LT(3,18,0)
      ierr = PetscObjectTypeCompare((PetscObject)(pA->A),MATTRANSPOSEMAT,&istrans);
-#else
-      ierr = PetscObjectTypeCompare((PetscObject)(pA->A),MATTRANSPOSEVIRTUAL,
-                                    &istrans);
-#endif
      CCHKERRQ(pA->GetComm(),ierr);
      if (!istrans)
      {
@@ -2855,43 +2850,37 @@ void PetscBCHandler::ZeroBC(const Vector &x, Vector &y)
 // PetscLinearSolver methods

 PetscLinearSolver::PetscLinearSolver(MPI_Comm comm, const std::string &prefix,
-                                     bool wrapin, bool iter_mode)
-   : PetscSolver(), Solver(0,iter_mode), wrap(wrapin)
+                                     bool wrapin)
+   : PetscSolver(), Solver(), wrap(wrapin)
 {
   KSP ksp;
   ierr = KSPCreate(comm,&ksp); CCHKERRQ(comm,ierr);
   obj  = (PetscObject)ksp;
   ierr = PetscObjectGetClassId(obj,&cid); PCHKERRQ(obj,ierr);
   ierr = KSPSetOptionsPrefix(ksp, prefix.c_str()); PCHKERRQ(ksp, ierr);
-   ierr = KSPSetInitialGuessNonzero(ksp, (PetscBool)iterative_mode);
-   PCHKERRQ(ksp, ierr);
 }

 PetscLinearSolver::PetscLinearSolver(const PetscParMatrix &A,
-                                     const std::string &prefix, bool iter_mode)
-   : PetscSolver(), Solver(0,iter_mode), wrap(false)
+                                     const std::string &prefix)
+   : PetscSolver(), Solver(), wrap(false)
 {
   KSP ksp;
   ierr = KSPCreate(A.GetComm(),&ksp); CCHKERRQ(A.GetComm(),ierr);
   obj  = (PetscObject)ksp;
   ierr = PetscObjectGetClassId(obj,&cid); PCHKERRQ(obj,ierr);
   ierr = KSPSetOptionsPrefix(ksp, prefix.c_str()); PCHKERRQ(ksp, ierr);
-   ierr = KSPSetInitialGuessNonzero(ksp, (PetscBool)iterative_mode);
-   PCHKERRQ(ksp, ierr);
   SetOperator(A);
 }

 PetscLinearSolver::PetscLinearSolver(const HypreParMatrix &A, bool wrapin,
-                                     const std::string &prefix, bool iter_mode)
-   : PetscSolver(), Solver(0,iter_mode), wrap(wrapin)
+                                     const std::string &prefix)
+   : PetscSolver(), Solver(), wrap(wrapin)
 {
   KSP ksp;
   ierr = KSPCreate(A.GetComm(),&ksp); CCHKERRQ(A.GetComm(),ierr);
   obj  = (PetscObject)ksp;
   ierr = PetscObjectGetClassId(obj, &cid); PCHKERRQ(obj, ierr);
   ierr = KSPSetOptionsPrefix(ksp, prefix.c_str()); PCHKERRQ(ksp, ierr);
-   ierr = KSPSetInitialGuessNonzero(ksp, (PetscBool)iterative_mode);
-   PCHKERRQ(ksp, ierr);
   SetOperator(A);
 }

@@ -3100,12 +3089,12 @@ void PetscLinearSolver::MultKernel(const Vector &b, Vector &x, bool trans) const
      }
   }
   B->PlaceMemory(b.GetMemory());
+   X->PlaceMemory(x.GetMemory(),iterative_mode);

   Customize();

-   PetscBool flg;
-   ierr = KSPGetInitialGuessNonzero(ksp, &flg);
-   X->PlaceMemory(x.GetMemory(),flg);
+   ierr = KSPSetInitialGuessNonzero(ksp, (PetscBool)iterative_mode);
+   PCHKERRQ(ksp, ierr);

   // Solve the system.
   if (trans)
@@ -3140,9 +3129,8 @@ PetscLinearSolver::~PetscLinearSolver()

 // PetscPCGSolver methods

-PetscPCGSolver::PetscPCGSolver(MPI_Comm comm, const std::string &prefix,
-                               bool iter_mode)
-   : PetscLinearSolver(comm,prefix,iter_mode)
+PetscPCGSolver::PetscPCGSolver(MPI_Comm comm, const std::string &prefix)
+   : PetscLinearSolver(comm,prefix)
 {
   KSP ksp = (KSP)obj;
   ierr = KSPSetType(ksp,KSPCG); PCHKERRQ(ksp,ierr);
@@ -3150,9 +3138,8 @@ PetscPCGSolver::PetscPCGSolver(MPI_Comm comm, const std::string &prefix,
   ierr = KSPSetNormType(ksp,KSP_NORM_NATURAL); PCHKERRQ(ksp,ierr);
 }

-PetscPCGSolver::PetscPCGSolver(PetscParMatrix& A, const std::string &prefix,
-                               bool iter_mode)
-   : PetscLinearSolver(A,prefix,iter_mode)
+PetscPCGSolver::PetscPCGSolver(PetscParMatrix& A, const std::string &prefix)
+   : PetscLinearSolver(A,prefix)
 {
   KSP ksp = (KSP)obj;
   ierr = KSPSetType(ksp,KSPCG); PCHKERRQ(ksp,ierr);
@@ -3161,8 +3148,8 @@ PetscPCGSolver::PetscPCGSolver(PetscParMatrix& A, const std::string &prefix,
 }

 PetscPCGSolver::PetscPCGSolver(HypreParMatrix& A, bool wrap,
-                               const std::string &prefix, bool iter_mode)
-   : PetscLinearSolver(A,wrap,prefix,iter_mode)
+                               const std::string &prefix)
+   : PetscLinearSolver(A,wrap,prefix)
 {
   KSP ksp = (KSP)obj;
   ierr = KSPSetType(ksp,KSPCG); PCHKERRQ(ksp,ierr);
@@ -744,16 +744,16 @@ private:

 public:
   PetscLinearSolver(MPI_Comm comm, const std::string &prefix = std::string(),
-                     bool wrap = true, bool iter_mode = false);
+                     bool wrap = true);
   PetscLinearSolver(const PetscParMatrix &A,
-                     const std::string &prefix = std::string(), bool iter_mode = false);
+                     const std::string &prefix = std::string());
   /// Constructs a solver using a HypreParMatrix.
   /** If @a wrap is true, then the MatMult ops of HypreParMatrix are wrapped.
       No preconditioner can be automatically constructed from PETSc. If
       @a wrap is false, the HypreParMatrix is converted into a the AIJ
       PETSc format, which is suitable for most preconditioning methods. */
   PetscLinearSolver(const HypreParMatrix &A, bool wrap = true,
-                     const std::string &prefix = std::string(), bool iter_mode = false);
+                     const std::string &prefix = std::string());
   virtual ~PetscLinearSolver();

   /// Sets the operator to be used for mat-vec operations and
@@ -780,12 +780,10 @@ public:
 class PetscPCGSolver : public PetscLinearSolver
 {
 public:
-   PetscPCGSolver(MPI_Comm comm, const std::string &prefix = std::string(),
-                  bool iter_mode = false);
-   PetscPCGSolver(PetscParMatrix &A, const std::string &prefix = std::string(),
-                  bool iter_mode = false);
-   PetscPCGSolver(HypreParMatrix &A, bool wrap = true,
-                  const std::string &prefix = std::string(), bool iter_mode = false);
+   PetscPCGSolver(MPI_Comm comm, const std::string &prefix = std::string());
+   PetscPCGSolver(PetscParMatrix &A, const std::string &prefix = std::string());
+   PetscPCGSolver(HypreParMatrix &A,bool wrap=true,
+                  const std::string &prefix = std::string());
 };


@@ -20,6 +20,9 @@
 #include <cmath>
 #include <set>

+#define MFEM_NVTX_COLOR DarkGreen
+#include "../general/nvtx.hpp"
+
 namespace mfem
 {

@@ -710,6 +713,8 @@ void CGSolver::UpdateVectors()

 void CGSolver::Mult(const Vector &b, Vector &x) const
 {
+   NVTX("CG");
+
   int i;
   double r0, den, nom, nom0, betanom, alpha, beta;

@@ -3034,7 +3039,7 @@ void BlockILU::Mult(const Vector &b, Vector &x) const
      }
      LUFactors A_ii_inv(&DB(0,0,i), &ipiv[i*block_size]);
      // x_i = D_ii^{-1} x_i
-      A_ii_inv.Solve(block_size, 1, xi.GetData());
+      A_ii_inv.Solve(block_size, 1, xi);
   }
 }

@@ -3197,8 +3202,7 @@ void UMFPackSolver::Mult(const Vector &b, Vector &x) const
   {
      int status =
         umfpack_di_solve(UMFPACK_At, mat->HostReadI(), mat->HostReadJ(),
-                          mat->HostReadData(), x.HostWrite(), b.HostRead(),
-                          Numeric, Control, Info);
+                          mat->HostReadData(), x, b, Numeric, Control, Info);
      umfpack_di_report_info(Control, Info);
      if (status < 0)
      {
@@ -3209,9 +3213,8 @@ void UMFPackSolver::Mult(const Vector &b, Vector &x) const
   else
   {
      SuiteSparse_long status =
-         umfpack_dl_solve(UMFPACK_At, AI, AJ, mat->HostReadData(),
-                          x.HostWrite(), b.HostRead(), Numeric, Control,
-                          Info);
+         umfpack_dl_solve(UMFPACK_At, AI, AJ, mat->HostReadData(), x, b,
+                          Numeric, Control, Info);
      umfpack_dl_report_info(Control, Info);
      if (status < 0)
      {
@@ -3232,8 +3235,7 @@ void UMFPackSolver::MultTranspose(const Vector &b, Vector &x) const
   {
      int status =
         umfpack_di_solve(UMFPACK_A, mat->HostReadI(), mat->HostReadJ(),
-                          mat->HostReadData(), x.HostWrite(), b.HostRead(),
-                          Numeric, Control, Info);
+                          mat->HostReadData(), x, b, Numeric, Control, Info);
      umfpack_di_report_info(Control, Info);
      if (status < 0)
      {
@@ -3245,9 +3247,8 @@ void UMFPackSolver::MultTranspose(const Vector &b, Vector &x) const
   else
   {
      SuiteSparse_long status =
-         umfpack_dl_solve(UMFPACK_A, AI, AJ, mat->HostReadData(),
-                          x.HostWrite(), b.HostRead(), Numeric, Control,
-                          Info);
+         umfpack_dl_solve(UMFPACK_A, AI, AJ, mat->HostReadData(), x, b,
+                          Numeric, Control, Info);
      umfpack_dl_report_info(Control, Info);
      if (status < 0)
      {
@@ -899,346 +899,6 @@ public:
   virtual void Mult(const Vector &xt, Vector &x) const;
 };

-/** Defines operators and constraints for the following nonlinear contact
- *  problem:
- *
- *    Find x in R^n and lambda in R^m such that
- *    K(x) - \nabla g(x)^T lambda = f
- *    0 <= g(x) \perp lambda >= 0
- *    
- *    Above, the perpendicularity ("\perp") equation is understood pointwise,
- *    i.e., g_i(x) lambda_i = 0, for all i=1,2,...,m.
- *
- *  This is class is to be derived by concrete contact problems and will be
- * called by the contact solver to evaluate the contact problem.
- */
-class ContactProblem
-{
-protected:
-   // Problem sizes
-   const int m, n;
-   // Right-hand side of the contact problem
-   mutable Vector f;
-  
-public:
-   /// Constructor taking the problem sizes as parameters.
-   ContactProblem(int m, int n);
-
-   inline int GetNumConstraints() const { return m; }
-   inline int GetNumDOFs () { return n; }
-
-   // [TODO:] an alternative design would be via a pure virtual method `Eval()`
-   // a members+getters for the problem's rhs, K, gradK, g, and gradG. The
-   // downside is that for nonlinear IPMs, sometimes g and K are evaluated without
-   // needing to evaluate their derivatives, gradG and gradK.
-  
-   /// Returns the rhs of the (first set of equations of the ) contact problem.
-   virtual void GetRhs(Vector &f_out) { f_out = f; }
-  
-   /** This method evaluates the nonlinear function K from the first set of contact
-    * equations. The contact solver calls this function. Concrete contact 
-    * problems need to implement this callback. The method should return false
-    * when the evaluation failed, otherwise should return true. */
-   virtual bool EvalK(const Vector &x, Vector &Katx) = 0;
-
-   /** This method evaluates the gradient of the function K at x. The contact 
-    * solver calls this function. Concrete contact problems need to implement this
-    * callback. The method should return false when the evaluation  failed, otherwise 
-    * should return true. */
-   virtual bool EvalGradK(const Vector &x, Operator &gradKatx) = 0;
-  
-   /** This method evaluates the contact gap function g at x. 
-    * The contact solver calls this function and concrete contact problems will
-    * implement this callback. The method should return false when the 
-    * evaluation failed, otherwise should return true. */
-   virtual bool Evalg(const Vector &x, Vector &g) = 0;
-
-   /** This method evaluates the gradient (nonlinear operator) of the gap 
-    * function g at x. The contact solver calls this function. Concrete contact 
-    * problems need to implement this callback. The method should return false
-    * when the evaluation failed, otherwise should return true. */
-   virtual bool EvalGradg(const Vector &x, Operator &gradG) = 0;
-
-   //TODO: Hessian of the Lagrangian
-};
-
-//forward declaration of the "contact" linear solver
-class ContactLinearSolver;
-
-/// Abstract solver for ContactProblems.
-class ContactSolver : public IterativeSolver
-{
-protected:
-   ContactProblem *problem;
-   ContactLinearSolver* linear_solver;
-public:
-   ContactSolver()
-     : IterativeSolver(), problem(nullptr), linear_solver(nullptr) { }
-#ifdef MFEM_USE_MPI
-   ContactSolver(MPI_Comm comm_)
-     : IterativeSolver(comm_), problem(nullptr), linear_solver(nullptr) { }
-#endif
-   virtual ~ContactSolver() { }
-  
-   /** This method is virtual as solvers might need to perform some initial
-    *  actions (e.g., validation) with the ContactProblem. */
-   virtual void SetContactProblem(ContactProblem &prob)
-   { problem = &prob; }
-
-   virtual void SetLinearSolver(ContactLinearSolver &lsolver)
-   { linear_solver = &lsolver; }
-  
-   /** This method performs the numerical solve of the complementarity 
-    * problem. Potential implementations of this purely virtual method 
-    * will be available via derived classes, e.g. implementations of 
-    * Uzawa (Augmented Lagrangian), interior-point method, etc. 
-    * 
-    * [TODO: work around the two input parameters to pass the initial 
-    * point and return the solution. Are there multivector vectors in 
-    * MFEM?] */
-   ///[TODO: clarify the const `Mult`;
-   virtual void Mult(const Vector &xt, Vector &x) const = 0;
-
-   virtual void SetPreconditioner(Solver &pr)
-   { MFEM_ABORT("Not meaningful for this solver."); }
-   virtual void SetOperator(const Operator &op)
-   { MFEM_ABORT("Not meaningful for this solver."); }
-};
-
-#if 0
-// illustration of the contact nonlinear (outer) loop and interplay with
-// linear solver and contact problem
-void ContactSolverIPM::Mult(const Vector &xin, Vector &xout) const
-{
-  x_iter.SetToZero();
-  s_iter.SetToZero();
-  l_iter.SetToZero();
-  R.SetToZero();
-  mu = 1.0;
-  int num_iter = 0;
-  
-  //solve loop
-  while(true) {
-    //evaluate contact problem at x_iter
-    problem->GetRhs(f);
-    problem->EvalK(x_iter, K);
-    problem->EvalGradK(x_iter, GradK);
-    problem->Evalg(x_iter, g);
-    problem->EvalGradG(x_iter, gradG);
-
-    //evaluate residuals of the contact problem (illustration for IPM)
-    // rx = f + gradG^T l_iter - K
-    // rs = s - g
-    // rc = mu*ones(m,1) - l_iter .* s_iter
-    this->ComputeResiduals(rx, rs, rc);
-
-    //check residual norms, max # of iterations, etc.
-    if(this->StoppingCriteriaMet(num_iter, rx, rs, rc)) {
-      break;
-    }
-    num_iter++;
-    
-    
-    //
-    //set up linear system and solve to obtain search
-    //direction/step  dx, ds, dl
-    //
-
-    //right hand side for linear system for the 
-    // predictor step rc_pred = - l_iter .* s_iter
-    this->ComputePredictorResidual(rc_pred)
-    
-    linear_solver->SetHessian(K);
-    linear_solver->SetGradG(gradG)
-    linear_solver->SetRhs(rx, rs, rc_pred);
-    linear_solver->SetLambda(l_iter);
-    linear_solver->SetS(s_iter);
-    linear_solver->SetPrimalReg(R);
-
-    linear_solver->Mult(dx, ds, dl);
-    //probably more like
-    //linear_solver->Mult(rx, rs, rc_pred, dx, ds, dl);
-    
-    //was the linear solve successful?
-    while(!lin_sol_ok) {
-      this->ComputeRegularization(R);
-      linear_solver->SetPrimalReg(R);
-      linear_solver->Mult(dx, ds, dl);
-      //check linear solve ok
-    }
-
-    //uses x_iter, s_iter, l_iter, dx, ds, dl
-    double alpha = this->LineSearchBlocking();
-    //uses x_iter, s_iter, l_iter, dx, ds, dl, alpha;
-    double sigma = this->ComputeCenteringParameter(alpha); 
-    
-    //right hand side for linear system for the 
-    // corrector step rc_pred = mu*sigma*ones(m,1) + ds.*dl - l_iter .* s_iter
-    //also uses s_iter, l_iter, ds, dl,
-    this->ComputeCorrectorResidual(mu, sigma, rc_corr);
-
-    linear_solver->SetRhs(rx, rs, rc_corr);
-    linear_solver->Mult(dx, ds, dl);
-    
-    //was the linear solve successful?
-    while(!lin_sol_ok) {
-      this->ComputeRegularization(R);
-      linear_solver->SetPrimalReg(R);
-      linear_solver->Mult(dx, ds, dl);
-      //check linear solve ok
-    }
-
-    //uses x_iter, s_iter, l_iter, dx, ds, dl
-    double alpha_corr = this->LineSearchCorrector();
-
-    //x = x+alpha*dx  s=s+alpha*ds  l=l+alpha*dl
-    this->UpdateIter(alpha_corr);
-  } //end of contact loop
-
-  //return x_iter, s_iter, l_iter as the solution
-  //[TODO:]
-}
-
-#endif
-
-#if 0
-//illustration of the "main" driver
-int main()
-{
-  //dream big
-  //
-  mfem::ContactProblem* problem = new mfem::apps::ContactProblemMortar(/*mesh_files*/);
-  //this is curently doable
-  //mfem::ContactProblem* problem = new mfem::apps::ContactLinearProblemFromFiles(/*matrix_files*/);
-
-  //linear solver based on options
-  mfem::ContactLinearSolver* lin_solver = nullptr;
-  if(that) {
-    lin_solver = new mfem::ContactLinearSolverCondensed();
-  } else {
-    //testing
-    lin_solver = new mfem::ContactLinearSolverCondensendDirect();
-  }
-  mfem::ContactSolver* solver = new mfem::ContactSolverIPM();
-
-  solver->SetContactProblem(*problem);
-  solver->SetLinearSolver(*lin_solver);
-
-  solver->Mult();
-  //error control = ?
-
-  //[TODO: what does the interface(s) need to accomdate solving a sequence of (related) problems
-  
-  //clean up
-  delete lin_solver;
-  delete problem;
-}
-#endif
-
-/** Abstract linear solver for internal use in ContactSolver. Implementations
- * of this class will solve linearizations of the contact problem in the form
- * shown below. The unknowns are vectors dx, ds, and dlambda
- * 
- *    (H+R)dx - gradG^T dlambda = r1
- *    gradG dx - ds = r2
- *    Lambda ds + S dlambda = r3
- * 
- * Above, H is an nxn symm. psd matrix, for example can be
- *    H = gradK(x)-\sum \lambda_i \nabla^2 g_i(x), 
- * or
- *    H = gradK(x),
- * depending on the contact solver. gradG is an mxn matrix and Lambda and
- * S are diagonal matrices with nonnegative entries. The regularization R 
- * is also a diagonal matrix with nonnegative entries.
- * 
- */
-class ContactLinearSolver : public Solver
-{
-public:
-  // [TODO: better names?]
-  void SetHessian(const Operator &gradK);
-  void SetRhs(const Vector &r1, const Vector &r2, const Vector &r3);
-  void SetGradG(const Operator& gradG);
-  void SetLambda(const Vector& Lambda);
-  void SetS(const Vector& S);
-  void SetPrimalReg(const Vector& R);
-  //to be continued...
-   
-   virtual void Mult(const Vector &x_in, Vector &x_out) const = 0;
-protected:
-};
-
-/** Implementation of the contact linear solve in the form of
- * 
- *    (H + R + gradG^T S^{-1} Lambda gradG) dx = 
- *                r1 + gradG^T S^{-1} (Lambda r2 + r3)
- *    ds = gradG dx - r2
- *    dlambda = -S^{-1} Lambda ds + S^{-1}r3
- * 
- * Essentially, the solve is a Gauss elimination of the linearization
- * described in ContactLinearSolver. The matrix from the first set
- * of equations is generally symm. positive definite. The solve
- * done by this class can/is allowed to fail, in which case the
- * subsequent linear solve will be done with a "more" p.d. R.
- * [TODO: ask about error control]
- *
- * This applies to both interior-point and Uzawa linearizations. To
- * use this class with Uzawa, one needs to set S to the identity 
- * matrix and Lambda to a diagonal matrix having entry (i,i) equal to
- * one when the corresponding gap is active (i.e., g_i(x)=0) or to 
- * zero  otherwise, for all i=1,2,...,m. [TODO: what about r3?]
- * 
- * The solver uses AMG [TODO: update me].
- */
-class ContactLinearSolverCondensed : public ContactLinearSolver
-{
-public:
-   virtual void Mult(const Vector &dx_in, Vector &dx_out) const
-   {
-     MFEM_ASSERT(false, "to be implemented");
-   }
-protected:
-};
-
-/** Implementation of the contact linear solve in the form of the
- * general symmetric indefinite linear system shown above in the
- * ContactLinearSolver class using direct linear solvers. This is
- * mostly for testing and debugging purposes.
- * 
- * [TODO: update me] Using a symmetric direct linear solver such
- * as MA57, Pardiso, WSMP, etc. would do it. LU solvers OK.
- */
-class ContactLinearSolverDirect : public ContactLinearSolver
-{
-public:
-   virtual void Mult(const Vector &x_in, Vector &x_out) const
-   {
-     MFEM_ASSERT(false, "to be implemented");
-   }
-protected:
-};
-
-/** Implementation of the contact linear solve in the form of the
- * symmetric p.d. linear system shown above in the
- * ContactLinearSolverCondensed class using direct linear solvers. 
- * This is mostly for testing and debugging purposes.
- * 
- * [TODO: update me] Using a Cholesky direct linear solver such
- * as Cholmod would do it. LU or symmetric indefinite direct 
- * solvers OK.
- * 
- */
-class ContactLinearSolverCondensendDirect : public ContactLinearSolverCondensed
-{
-public:
-   virtual void Mult(const Vector &x_in, Vector &x_out) const
-   {
-     MFEM_ASSERT(false, "to be implemented");
-   }
-protected:
-};
-
-
 /** Block ILU solver:
 *  Performs a block ILU(k) approximate factorization with specified block
 *  size. Currently only k=0 is supported. This is useful as a preconditioner
@@ -15,7 +15,6 @@
 #include "../general/forall.hpp"
 #include "../general/table.hpp"
 #include "../general/sort_pairs.hpp"
-#include "../general/backends.hpp"

 #include <iostream>
 #include <iomanip>
@@ -463,108 +462,25 @@ void SparseMatrix::SortColumnIndices()
      return;
   }

-#ifdef MFEM_USE_CUDA_OR_HIP
-   if ( Device::Allows( Backend::CUDA_MASK ))
+   const int * Ip=HostReadI();
+   HostReadWriteJ();
+   HostReadWriteData();
+
+   Array<Pair<int,double> > row;
+   for (int j = 0, i = 0; i < height; i++)
   {
-#if defined(MFEM_USE_CUDA)
-      size_t pBufferSizeInBytes = 0;
-      void *pBuffer = NULL;
-
-      const int n = Height();
-      const int m = Width();
-      const int nnzA = J.Capacity();
-      double * d_a_sorted = ReadWriteData();
-      const int * d_ia = ReadI();
-      int * d_ja_sorted = ReadWriteJ();
-      csru2csrInfo_t sortInfoA;
-
-      cusparseMatDescr_t matA_descr;
-      cusparseCreateMatDescr( &matA_descr );
-      cusparseSetMatIndexBase( matA_descr, CUSPARSE_INDEX_BASE_ZERO );
-      cusparseSetMatType( matA_descr, CUSPARSE_MATRIX_TYPE_GENERAL );
-
-      cusparseCreateCsru2csrInfo( &sortInfoA );
-
-      cusparseDcsru2csr_bufferSizeExt( handle, n, m, nnzA, d_a_sorted, d_ia,
-                                       d_ja_sorted, sortInfoA,
-                                       &pBufferSizeInBytes);
-
-      CuMemAlloc( &pBuffer, pBufferSizeInBytes );
-
-      cusparseDcsru2csr( handle, n, m, nnzA, matA_descr, d_a_sorted, d_ia,
-                         d_ja_sorted, sortInfoA, pBuffer);
-
-      cusparseDestroyCsru2csrInfo( sortInfoA );
-      cusparseDestroyMatDescr( matA_descr );
-
-      CuMemFree( pBuffer );
-#endif
-   }
-   else if ( Device::Allows( Backend::HIP_MASK ))
-   {
-#if defined(MFEM_USE_HIP)
-      size_t pBufferSizeInBytes = 0;
-      void *pBuffer = NULL;
-      int *P = NULL;
-
-      const int n = Height();
-      const int m = Width();
-      const int nnzA = J.Capacity();
-      double * d_a_sorted = ReadWriteData();
-      const int * d_ia = ReadI();
-      int * d_ja_sorted = ReadWriteJ();
-
-      hipsparseMatDescr_t descrA;
-      hipsparseCreateMatDescr( &descrA );
-      // FIXME: There is not in-place version of csr sort in hipSPARSE currently, so we make
-      //        a temporary copy of the data for gthr, sort that, and then copy the sorted values
-      //        back to the array being returned. Where there is an in-place version available,
-      //        we should use it.
-      Array< double > a_tmp( nnzA );
-      double *d_a_tmp = a_tmp.Write();
-
-      hipsparseXcsrsort_bufferSizeExt(handle, n, m, nnzA, d_ia, d_ja_sorted,
-                                      &pBufferSizeInBytes);
-
-      HipMemAlloc( &pBuffer, pBufferSizeInBytes );
-      HipMemAlloc( (void**)&P, nnzA * sizeof(int) );
-
-      hipsparseCreateIdentityPermutation(handle, nnzA, P);
-      hipsparseXcsrsort(handle, n, m, nnzA, descrA, d_ia, d_ja_sorted, P, pBuffer);
-
-      hipsparseDgthr(handle, nnzA, d_a_sorted, d_a_tmp, P,
-                     HIPSPARSE_INDEX_BASE_ZERO);
-
-      A.CopyFrom( a_tmp.GetMemory(), nnzA );
-      hipsparseDestroyMatDescr( descrA );
-
-      HipMemFree( pBuffer );
-      HipMemFree( P );
-#endif
-   }
-   else
-#endif // MFEM_USE_CUDA_OR_HIP
-   {
-      const int * Ip=HostReadI();
-      HostReadWriteJ();
-      HostReadWriteData();
-
-      Array<Pair<int,double> > row;
-      for (int j = 0, i = 0; i < height; i++)
+      int end = Ip[i+1];
+      row.SetSize(end - j);
+      for (int k = 0; k < row.Size(); k++)
      {
-         int end = Ip[i+1];
-         row.SetSize(end - j);
-         for (int k = 0; k < row.Size(); k++)
-         {
-            row[k].one = J[j+k];
-            row[k].two = A[j+k];
-         }
-         row.Sort();
-         for (int k = 0; k < row.Size(); k++, j++)
-         {
-            J[j] = row[k].one;
-            A[j] = row[k].two;
-         }
+         row[k].one = J[j+k];
+         row[k].two = A[j+k];
+      }
+      row.Sort();
+      for (int k = 0; k < row.Size(); k++, j++)
+      {
+         J[j] = row[k].one;
+         A[j] = row[k].two;
      }
   }
   isSorted = true;
@@ -23,11 +23,7 @@
 #include "densemat.hpp"

 #if defined(MFEM_USE_HIP)
-#if (HIP_VERSION_MAJOR * 100 + HIP_VERSION_MINOR) < 502
 #include <hipsparse.h>
-#else
-#include <hipsparse/hipsparse.h>
-#endif
 #endif


@@ -56,9 +56,7 @@ STRUMPACKRowLocMatrix::STRUMPACKRowLocMatrix(const HypreParMatrix & hypParMat)

   // Create the CSRMatrixMPI A_ by borrowing the internal data from a
   // hypre_CSRMatrix.
-   hypParMat.HostRead();
   hypre_CSRMatrix * csr_op = hypre_MergeDiagAndOffd(parcsr_op);
-   hypParMat.HypreRead();
   hypre_CSRMatrixSetDataOwner(csr_op,0);
 #if MFEM_HYPRE_VERSION >= 21600
   // For now, this method assumes that HYPRE_Int is int. Also, csr_op->num_cols
@@ -203,8 +201,8 @@ void STRUMPACKSolver::Mult( const Vector & x, Vector & y ) const
   MFEM_ASSERT(y.Size() == Height(), "invalid y.Size() = " << y.Size()
               << ", expected size = " << Height());

-   double*  yPtr = y.HostWrite();
-   const double*  xPtr = x.HostRead();
+   double*  yPtr = (double*)y;
+   double*  xPtr = (double*)(const_cast<Vector&>(x));

   solver_->options().set_verbose( factor_verbose_ );
   ReturnCode ret = solver_->factor();
@@ -221,10 +219,6 @@ void STRUMPACKSolver::Mult( const Vector & x, Vector & y ) const
         MFEM_ABORT("STRUMPACK:  Matrix reordering failed!");
      }
      break;
-      default:
-      {
-         MFEM_ABORT("STRUMPACK: 'factor()' error code = " << ret);
-      }
   }
   solver_->options().set_verbose( solve_verbose_ );
   solver_->solve(xPtr, yPtr);
@@ -149,9 +149,7 @@ SuperLURowLocMatrix::SuperLURowLocMatrix( const HypreParMatrix & hypParMat )

   // Create the SuperMatrix A by borrowing the internal data from a
   // hypre_CSRMatrix.
-   hypParMat.HostRead();
   hypre_CSRMatrix * csr_op = hypre_MergeDiagAndOffd(parcsr_op);
-   hypParMat.HypreRead();
   hypre_CSRMatrixSetDataOwner(csr_op,0);
 #if MFEM_HYPRE_VERSION >= 21600
   // For now, this method assumes that HYPRE_BigInt is int. Also, csr_op->num_cols
--- a/Show More
+++ b/Show More
Author	SHA1	Message	Date
camierjs	58a2c7ac15	Merge master in lor-gpu	2022-10-23 09:54:23 -07:00
camierjs	2a04311099	LOR H1/ND/RT explicit specializations	2022-09-15 07:28:55 -07:00
camierjs	026c2ae2dd	Add lor_h1 batched kernel for order 1	2022-09-13 17:26:27 -07:00
Will Pazner	b9de6d2b37	Fix issue with memory reuse	2022-09-13 12:56:41 -07:00
Will Pazner	1cca575225	Merge remote-tracking branch 'origin/lor-gpu' into lor-gpu	2022-09-08 16:42:02 -07:00
Will Pazner	a44abee1ba	Fix some issues with CEED	2022-09-08 16:33:48 -07:00
Will Pazner	55497c158f	Add script to run AMR case	2022-09-08 16:33:14 -07:00
camierjs	39185e89bc	Bring simplified lor_h1 kernel	2022-09-06 10:10:59 -07:00
Will Pazner	cc6b3c025d	Merge remote-tracking branch 'origin/quadrature-coeff-lor' into lor-gpu	2022-09-02 16:02:17 -07:00
Will Pazner	e2be4b8d9b	Merge remote-tracking branch 'origin/face-quadrature-space-coefficient' into lor-gpu	2022-09-02 16:01:54 -07:00
Will Pazner	120963609b	Fix comment	2022-09-02 08:12:56 -07:00
Will PaznerandYohann Dudouit	ea310df39d	Replace some MFEM_ASSERT with MFEM_VERIFY Co-authored-by: Yohann Dudouit <dudouit1@llnl.gov>	2022-09-02 08:06:49 -07:00
Will Pazner	bae885021f	make style	2022-09-01 17:01:11 -07:00
Will Pazner	798f529929	NVTX	2022-09-01 16:57:28 -07:00
Will Pazner	72e5222a98	make style	2022-08-27 17:13:45 -07:00
Will Pazner	db9772acd4	Bench LOR AMS and ADS	2022-08-27 17:13:05 -07:00
Will Pazner	7a09a3d772	Reuse memory in discrete gradient and curl	2022-08-27 17:12:52 -07:00
Will Pazner	45402504ea	Allow memory reuse in FormLORVertexCoordinates	2022-08-27 17:12:34 -07:00
Will Pazner	24fd0b1c1e	Reuse memory in LOR benchmarks	2022-08-25 13:32:31 -07:00
Will Pazner	59fa852dc6	ND, RT, AMS, ADS LOR benchmarking	2022-08-25 11:40:52 -07:00
Will Pazner	eee4d9c25c	Specializations	2022-08-24 14:18:43 -07:00
Will Pazner	218c56945d	Add AMG to benchmarks	2022-08-24 13:56:30 -07:00
Will Pazner	7d31215c91	LOR benchmark	2022-08-24 13:37:17 -07:00
Will Pazner	f378aed68d	gitignore	2022-08-24 13:37:08 -07:00
Will Pazner	e4318fe0fb	Increase MD and MQ limits	2022-08-24 13:36:10 -07:00
Will Pazner	c504463321	Add warm up in LOR benchmark	2022-08-24 12:00:13 -07:00
Will Pazner	cba110b1bc	Reuse memory in LOR assembly	2022-08-24 12:00:02 -07:00
Will Pazner	474aa14907	Work on LOR benchmark	2022-08-24 11:40:17 -07:00
Will Pazner	86265d0848	Output total number of elements in plor_solvers	2022-08-24 09:48:21 -07:00
Will Pazner	152cae0c0e	Don't add NVTX for early-return hypre setup	2022-08-24 09:47:49 -07:00
Will Pazner	f98c2b9df5	Ignore more nsys files	2022-08-24 09:47:20 -07:00
Will Pazner	049ba615b3	Don't compute error unless also visualizing	2022-08-23 15:00:17 -07:00
Will Pazner	3f8c7aa5e1	Ignore nsys-rep files	2022-08-23 14:55:30 -07:00
Will Pazner	dd6164a9b4	Separate AMG setup in plor_solvers	2022-08-23 14:54:19 -07:00
Will Pazner	6974ebd250	Turn off visualization by default in plor_solvers	2022-08-23 14:39:36 -07:00
Will Pazner	947c51ff75	NVTX instrumentation	2022-08-23 14:39:24 -07:00
Will Pazner	d754cf4914	Add LOR benchmark	2022-08-23 07:49:44 -07:00
Will Pazner	5aa5cb551e	Add debug.hpp and nvtx.hpp	2022-08-23 07:48:54 -07:00
Will Pazner	7ce22d38c5	Support general coefficient types in batched LOR	2022-08-05 16:30:21 -07:00