Merge branch 'master' into glvis_stream

Add weighted LOR transfer

CHANGELOG

@@ -79,6 +79,11 @@ Linear and nonlinear solvers
   PRefinement multigrid methods for problems posed on trace spaces (see e.g. the
   DPG miniapps).
 
+- Added new class MultiVector: an array of Vectors of different sizes where each
+  Vector can be allocated independently. Also, added associated methods in class
+  Operator: MultMV, MultTransposeMV, and GetGradientMV, that use MultiVector
+  objects for input and/or output parameters. [PR #5249]
+
 GPU computing
 -------------
 - Improved partial assembly for VectorDivergenceIntegrator with shared-memory
@@ -92,6 +97,22 @@ GPU computing
 
 - Added device assembly support for 3D H(curl) VectorFEDomainLFIntegrator.
 
+- Added partial assembly support for MixedScalarWeakGradientIntegrator.
+
+- Added partial assembly support for MixedDotProductIntegrator.
+
+- Added partial assembly support for MixedScalarCrossProductIntegrator.
+
+- Added partial assembly support for MixedScalarWeakCrossProductIntegrator.
+
+- Added partial assembly support for MixedVectorGradientIntegrator for H1->RT.
+
+- Added support for device partial assembly CurlInterpolator.
+  This supports 2D and 3D variants:
+    2D H1 (out-of-plane) to RT (in-plane)
+    2D ND (in-plane) to Integral L2 (out-of-plane)
+    3D ND to RT
+
 - Added NVIDIA cuDSS library interface. Implementation examples have been
   added to ex1 and ex1p. See https://developer.nvidia.com/cudss for more
   details. Supported versions >= 0.6.0.
@@ -104,6 +125,9 @@ GPU computing
 - Added support for FiniteElement::MapType::INTEGRAL spaces to
   QuadratureInterpolator.
 
+- Added support for FiniteElement::MapType::INTEGRAL spaces to
+  MixedScalarCurlIntegrator.
+
 New and updated examples and miniapps
 -------------------------------------
 - The Lorentz miniapp (in miniapps/electromagnetics) has been updated to
@@ -118,6 +142,20 @@ Miscellaneous
   using the new method ApplyDofSigns() in class ParFiniteElementSpace: the
   method will return immediately if no sign flips are needed.
 
+- Added support for coefficient-weighted LOR transfer in
+  L2ProjectionGridTransfer. The transfer conserves the weighted mass, for
+  example when transferring velocity while conserving density-weighted momentum.
+  This is illustrated in the lor-transfer and plor-transfer miniapps.
+
+- Added support for saving DataCollection output on the node-local storage,
+  instead of requiring that the filesystem is shared among all the ranks.
+
+API changes
+-----------
+- Removed ProjectGrad from 2D RT elements. Users should use ProjectCurl instead.
+  This also fixes a bug where ProjectCurl was returning the negative curl,
+  identical to ProjectGrad.
+
 
 Version 4.9, released on Dec 11, 2025
 =====================================

CMakeLists.txt

@@ -88,18 +88,9 @@ if (MFEM_USE_STRUMPACK OR MFEM_USE_MUMPS)
   # Just needed to find the MPI_Fortran libraries to link with
   set(XSDK_ENABLE_Fortran ON)
 endif()
-# Ginkgo requires C++17:
-if ((MFEM_USE_GINKGO) AND ("${CMAKE_CXX_STANDARD}" LESS "17"))
-  set(CMAKE_CXX_STANDARD 17 CACHE STRING "C++ standard to use." FORCE)
-# Google Benchmark, SUNDIALS, STRUMPACK, Tribol, RAJA and Umpire require C++14:
-elseif ((MFEM_USE_BENCHMARK OR
-         MFEM_USE_SUNDIALS OR
-         MFEM_USE_STRUMPACK OR
-         MFEM_USE_TRIBOL OR
-         MFEM_USE_RAJA OR
-         MFEM_USE_UMPIRE) AND
-    ("${CMAKE_CXX_STANDARD}" LESS "14"))
-  set(CMAKE_CXX_STANDARD 14 CACHE STRING "C++ standard to use." FORCE)
+# RAJA requires C++20:
+if ((MFEM_USE_UMPIRE OR MFEM_USE_RAJA) AND ("${CMAKE_CXX_STANDARD}" LESS "20"))
+  set(CMAKE_CXX_STANDARD 20 CACHE STRING "C++ standard to use." FORCE)
 endif()
 
 # Include xSDK default CMake file.
@@ -608,6 +599,11 @@ if (MFEM_USE_ENZYME)
     set(ENZYME_INCLUDE_DIRS ${ENZYME_DIR}/include)
 endif()
 
+# GLVis
+if (MFEM_USE_GLVIS AND NOT MFEM_FETCH_GLVIS)
+  find_package(GLVis REQUIRED)
+endif()
+
 # MFEM_TIMER_TYPE
 if (NOT DEFINED MFEM_TIMER_TYPE)
   if (APPLE)
@@ -649,7 +645,7 @@ set(MFEM_TPLS OPENMP HYPRE LAPACK BLAS SuperLUDist STRUMPACK METIS SuiteSparse
     NETCDF MPFR PUMI HIOP POSIXCLOCKS MFEMBacktrace ZLIB OCCA CEED RAJA UMPIRE
     ADIOS2 MKL_CPARDISO MKL_PARDISO AMGX MAGMA CUSPARSE CUBLAS CUDSS CALIPER CODIPACK
     BENCHMARK PARELAG TRIBOL MPI_CXX HIP HIPBLAS HIPSPARSE MOONOLITH BLITZ
-    ALGOIM ENZYME CUDA::cudart)
+    ALGOIM ENZYME GLVIS CUDA::cudart)
 
 # Add all created targets and *_FOUND libraries in the variables TPL_TARGETS and
 # TPL_LIBRARIES, respectively.
@@ -829,6 +825,11 @@ endif()
 
 set(MFEM_CUSTOM_TARGET_PREFIX CACHE STRING "")
 
+if (MFEM_USE_GLVIS AND MFEM_FETCH_GLVIS)
+    # needs to be after mfem_add_library(mfem) to be able to depend on it
+    find_package(GLVis REQUIRED)
+endif()
+
 #-------------------------------------------------------------------------------
 # Examples, miniapps, benchmarks and testing
 #-------------------------------------------------------------------------------
@@ -929,6 +930,10 @@ target_include_directories(mfem BEFORE
   PUBLIC
     $<INSTALL_INTERFACE:${INSTALL_INCLUDE_DIR}>)
 
+if (MFEM_USE_GLVIS AND MFEM_FETCH_GLVIS)
+    target_link_libraries(mfem PUBLIC GLVIS)
+endif()
+
 # The 'install' target will not depend on 'all'.
 # set(CMAKE_SKIP_INSTALL_ALL_DEPENDENCY TRUE)
 

INSTALL

@@ -626,6 +626,9 @@ MFEM_USE_ENZYME = YES/NO
    config/defaults.mk. For more detailed instructions, see the section "Specific
    options for Enzyme" below.
 
+MFEM_USE_GLVIS = YES/NO
+   Enables using a GLVis stream directly instead of the socket one.
+
 MFEM_BUILD_TAG = (any value)
    An optional tag to characterize the build. Exported to config/config.mk.
    Can be used to identify the MFEM build from other makefiles.
@@ -1092,6 +1095,7 @@ MFEM_USE_BENCHMARK
 MFEM_USE_PARELAG
 MFEM_USE_TRIBOL
 MFEM_USE_ENZYME
+MFEM_USE_GLVIS
 
 The following options are CMake specific:
 
@@ -1102,6 +1106,7 @@ MFEM_FETCH_TPLS      - Enable fetching of all supported third-party libraries.
 MFEM_FETCH_GSLIB     - Enable fetching of gslib.
 MFEM_FETCH_HYPRE     - Enable fetching of hypre.
 MFEM_FETCH_METIS     - Enable fetching of metis.
+MFEM_FETCH_GLVIS     - Enable fetching of GLVis.
 
 External libraries (CMake):
 ---------------------------

config/cmake/config.hpp.in

@@ -222,4 +222,7 @@
 // Enable Enzyme for AD
 #cmakedefine MFEM_USE_ENZYME
 
+// Enable GLVis
+#cmakedefine MFEM_USE_GLVIS
+
 #endif // MFEM_CONFIG_HEADER

config/cmake/modules/FindGLVis.cmake

@@ -0,0 +1,135 @@
+# Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+# Defines the following variables:
+#   - GLVIS_FOUND
+#   - GLVIS_LIBRARIES
+#   - GLVIS_INCLUDE_DIRS
+
+if (MFEM_FETCH_GLVIS OR MFEM_FETCH_TPLS)
+    message(STATUS "GLVis: Fetch/ExternalProject")
+
+    get_directory_property(COMPILE_OPTS COMPILE_OPTIONS)
+    string(REPLACE ";" " " COMPILE_CXX_FLAGS "${COMPILE_OPTS}")
+
+    # get_directory_property(COMPILE_DEFINITIONS COMPILE_DEFINITIONS)
+    # string(REPLACE "\"" "\\\"" COMPILE_DEFS_QUOTED_STR "${COMPILE_DEFINITIONS}")
+    # string(REPLACE ";" " -D" COMPILE_DEFS_STR "-D${COMPILE_DEFS_QUOTED_STR}")
+    # string(JOIN " " COMPILE_CXX_FLAGS ${COMPILE_OPTS_STR} ${COMPILE_DEFS_STR})
+
+    cmake_host_system_information(RESULT NCPU QUERY NUMBER_OF_LOGICAL_CORES)
+
+    add_library(GLVIS STATIC IMPORTED)
+
+    include(ExternalProject)
+
+    set(FETCH_DIR ${CMAKE_CURRENT_BINARY_DIR}/fetch)
+    set(FETCH_GLVIS "${FETCH_DIR}/glvis")
+
+    ExternalProject_Add(glvis
+        GIT_REPOSITORY https://github.com/GLVis/glvis.git
+        GIT_TAG stream_sessions
+        GIT_SHALLOW TRUE
+        UPDATE_DISCONNECTED TRUE
+        CMAKE_GENERATOR ${CMAKE_GENERATOR}
+        PREFIX ${FETCH_GLVIS}
+        SOURCE_DIR ${FETCH_GLVIS}/src
+        STAMP_DIR ${FETCH_GLVIS}/stamp
+        BINARY_DIR ${FETCH_GLVIS}/build
+        DEPENDS mfem
+        CMAKE_ARGS
+            # -DCMAKE_VERBOSE_MAKEFILE=ON
+            -DMFEM_DIR=${CMAKE_CURRENT_BINARY_DIR}
+            -DCMAKE_BUILD_TYPE=${CMAKE_BUILD_TYPE}
+            -DCMAKE_CXX_COMPILER=${CMAKE_CXX_COMPILER}
+            -DCMAKE_CXX_FLAGS:STRING=${COMPILE_CXX_FLAGS}
+            -DGLVIS_BUILD_LIB_ONLY=ON
+        BUILD_COMMAND ${CMAKE_COMMAND}
+            --build ${FETCH_GLVIS}/build
+            --config $<CONFIG>
+            --target glvis glvis_logo
+            --parallel ${NCPU}
+        BUILD_BYPRODUCTS 
+            ${FETCH_GLVIS}/build/lib/libglvis.a
+            ${FETCH_GLVIS}/build/share/libglvis_logo.a
+        INSTALL_COMMAND "")
+        
+    set_target_properties(GLVIS PROPERTIES
+        IMPORTED_LOCATION ${FETCH_GLVIS}/build/lib/libglvis.a)
+
+    find_package(OpenGL REQUIRED)
+    find_package(GLEW   REQUIRED)
+    find_package(SDL2   REQUIRED)
+    find_package(PNG    REQUIRED)
+    find_package(Freetype REQUIRED)
+    find_package(Fontconfig REQUIRED)
+    if(APPLE)
+        find_library(COCOA_LIBRARY Cocoa)
+    endif()
+    target_link_libraries(GLVIS INTERFACE
+        ${FETCH_GLVIS}/build/share/libglvis_logo.a
+        OpenGL::GL
+        GLEW::GLEW
+        SDL2::SDL2
+        PNG::PNG
+        Freetype::Freetype
+        Fontconfig::Fontconfig)
+    if(APPLE)
+        target_link_libraries(GLVIS INTERFACE ${COCOA_LIBRARY})
+    endif()
+
+    set(GLVIS_FOUND TRUE)
+    return()
+endif()
+
+message(STATUS "[🔵 GLVis 🔵] Find pre-installed package")
+
+include(MfemCmakeUtilities)
+mfem_find_package(GLVis GLVIS GLVIS_DIR
+        "include" "lib/glwindow.hpp"
+        "lib" "build/lib/libglvis.a"
+        "Paths to headers required by GLVis"
+        "Libraries required by GLVis")
+
+if (GLVIS_FOUND)
+        set(GLVIS_INCLUDE_DIRS ${GLVIS_INCLUDE_DIRS}/lib)
+
+        find_library(GLVIS_LOGO_LIBRARY
+                NAMES glvis_logo
+                PATHS ${GLVIS_DIR}/build/share
+                NO_DEFAULT_PATH 
+                REQUIRED)
+        list(APPEND GLVIS_LIBRARIES ${GLVIS_LOGO_LIBRARY})
+
+        find_package(OpenGL REQUIRED)
+        list(APPEND GLVIS_LIBRARIES OpenGL::GL)
+
+        find_package(GLEW REQUIRED)
+        list(APPEND GLVIS_LIBRARIES GLEW::GLEW)
+        
+        find_package(SDL2 REQUIRED)
+        list(APPEND GLVIS_LIBRARIES SDL2::SDL2)
+
+        find_package(PNG REQUIRED)
+        list(APPEND GLVIS_LIBRARIES PNG::PNG)
+
+        find_package(Freetype REQUIRED)
+        list(APPEND GLVIS_LIBRARIES Freetype::Freetype)
+
+        find_package(Fontconfig REQUIRED)
+        list(APPEND GLVIS_LIBRARIES Fontconfig::Fontconfig)
+
+        find_library(COCOA_LIBRARY Cocoa)
+        list(APPEND GLVIS_LIBRARIES ${COCOA_LIBRARY})
+endif()
+        
+message(STATUS "GLVIS_INCLUDE_DIRS: ${GLVIS_INCLUDE_DIRS}")
+message(STATUS "GLVIS_LIBRARIES: ${GLVIS_LIBRARIES}")

config/cmake/modules/MfemCmakeUtilities.cmake

@@ -884,7 +884,7 @@ function(mfem_export_mk_files)
       MFEM_USE_ADIOS2 MFEM_USE_MKL_CPARDISO MFEM_USE_MKL_PARDISO
       MFEM_USE_ADFORWARD MFEM_USE_CODIPACK MFEM_USE_BENCHMARK MFEM_USE_PARELAG
       MFEM_USE_TRIBOL MFEM_USE_MOONOLITH MFEM_USE_ALGOIM MFEM_USE_ENZYME
-      MFEM_USE_HDF5)
+      MFEM_USE_HDF5 MFEM_USE_GLVIS)
   foreach(var ${CONFIG_MK_BOOL_VARS})
     if (${var})
       set(${var} YES)

config/config.hpp

@@ -18,7 +18,7 @@
 #define MFEM_CONFIG_HPP
 
 #ifdef MFEM_CONFIG_FILE
-#include MFEM_CONFIG_FILE
+#include MFEM_CONFIG_FILE // IWYU pragma: export
 #else
 #include "_config.hpp"
 #endif

config/config.hpp.in

@@ -222,4 +222,7 @@
 // Enable the Enzyme LLVM plugin
 // #define MFEM_USE_ENZYME
 
+// Enable GLVis.
+// #define MFEM_USE_GLVIS
+
 #endif // MFEM_CONFIG_HEADER

config/config.mk.in

@@ -72,6 +72,7 @@ MFEM_USE_BENCHMARK     = @MFEM_USE_BENCHMARK@
 MFEM_USE_PARELAG       = @MFEM_USE_PARELAG@
 MFEM_USE_TRIBOL        = @MFEM_USE_TRIBOL@
 MFEM_USE_ENZYME        = @MFEM_USE_ENZYME@
+MFEM_USE_GLVIS         = @MFEM_USE_GLVIS@
 
 # Compiler, compile options, and link options
 MFEM_CXX       = @MFEM_CXX@

config/defaults.cmake

@@ -72,6 +72,7 @@ option(MFEM_USE_BENCHMARK "Enable Google Benchmark" OFF)
 option(MFEM_USE_PARELAG "Enable ParELAG" OFF)
 option(MFEM_USE_TRIBOL "Enable Tribol" OFF)
 option(MFEM_USE_ENZYME "Enable Enzyme" OFF)
+option(MFEM_USE_GLVIS "Enable GLVis" OFF)
 
 # Optional overrides for autodetected MPIEXEC and MPIEXEC_NUMPROC_FLAG
 # set(MFEM_MPIEXEC "mpirun" CACHE STRING "Command for running MPI tests")
@@ -96,6 +97,7 @@ option(MFEM_FETCH_TPLS "Enable fetching of all supported third-party libraries"
 option(MFEM_FETCH_GSLIB "Enable fetching of GSLIB" OFF)
 option(MFEM_FETCH_HYPRE "Enable fetching of hypre" OFF)
 option(MFEM_FETCH_METIS "Enable fetching of METIS" OFF)
+option(MFEM_FETCH_GLVIS "Enable fetching of GLVis" OFF)
 
 # Setting CXX/MPICXX on the command line or in user.cmake will overwrite the
 # autodetected C++ compiler.
@@ -278,6 +280,8 @@ set(Tribol_REQUIRED_PACKAGES "Axom/core/mint/slam/slic" CACHE STRING
 
 set(ENZYME_DIR "${MFEM_DIR}/../enzyme" CACHE PATH "Path to Enzyme")
 
+set(GLVIS_DIR "${MFEM_DIR}/../glvis" CACHE PATH "Path to GLVis")
+
 set(BLAS_INCLUDE_DIRS "" CACHE STRING "Path to BLAS headers.")
 set(BLAS_LIBRARIES "" CACHE STRING "The BLAS library.")
 set(LAPACK_INCLUDE_DIRS "" CACHE STRING "Path to LAPACK headers.")

config/defaults.mk

@@ -28,11 +28,8 @@ MPICXX = mpicxx
 BASE_FLAGS  = -std=c++17
 OPTIM_FLAGS = -O3 $(BASE_FLAGS)
 
-# Shadow warnings for clang only; GCC's -Wshadow flags more.
-SHADOW_WARNING_FLAG = $(if $(findstring clang,\
- $(shell $(MFEM_HOST_CXX) --version 2>/dev/null)),-Wshadow,)
-WARNING_FLAGS = -pedantic -Wall $(SHADOW_WARNING_FLAG)
-
+# The variable WARNING_FLAGS depends on which compiler is used, and is defined
+# later in this file.
 DEBUG_FLAGS = $(strip -g $(addprefix $(XCOMPILER),$(WARNING_FLAGS)) $(BASE_FLAGS))
 
 # Prefixes for passing flags to the compiler and linker when using CXX or MPICXX
@@ -52,6 +49,10 @@ SHARED = NO
 #
 # If you set MFEM_USE_ENZYME=YES, must use CUDA_CXX=clang++
 CUDA_CXX = nvcc
+# CUDA compute capability used during compilation, e.g. sm_60. Multiple
+# architectures can be requested as a comma-separated list, e.g. sm_70,sm_80.
+# A single value may also be one of the nvcc special values "all",
+# "all-major", or "native".
 CUDA_ARCH = sm_60
 # Base CUDA install directory, only needed if building with clang+cuda:
 # The default setting is:
@@ -60,11 +61,23 @@ CUDA_ARCH = sm_60
 # 3. Use /usr/local/cuda
 CUDA_DIR = $(or $(CUDA_HOME),$(patsubst %/,%,$(dir \
  $(patsubst %/,%,$(dir $(shell command -v nvcc))))),/usr/local/cuda)
+# Derive nvcc/clang architecture flags from CUDA_ARCH. A comma-separated list
+# expands into one -gencode / --cuda-gpu-arch flag per architecture; otherwise
+# use the -arch / --cuda-gpu-arch shorthand.
+MFEM_COMMA := ,
+CUDA_ARCH_NUMS = $(patsubst sm_%,%,$(subst $(MFEM_COMMA), ,$(CUDA_ARCH)))
+NVCC_ARCH_FLAGS = $(strip $(if $(findstring $(MFEM_COMMA),$(CUDA_ARCH)),\
+ $(foreach arch,$(CUDA_ARCH_NUMS),\
+ -gencode arch=compute_$(arch)$(MFEM_COMMA)code=sm_$(arch)),\
+ -arch=$(CUDA_ARCH)))
+CLANG_ARCH_FLAGS = $(strip $(if $(findstring $(MFEM_COMMA),$(CUDA_ARCH)),\
+ $(foreach arch,$(CUDA_ARCH_NUMS),--cuda-gpu-arch=sm_$(arch)),\
+ --cuda-gpu-arch=$(CUDA_ARCH)))
 # flags for clang+cuda
-CLANG_CUDA_FLAGS = -xcuda --cuda-path=$(CUDA_DIR) --cuda-gpu-arch=$(CUDA_ARCH)
+CLANG_CUDA_FLAGS = -xcuda --cuda-path=$(CUDA_DIR) $(CLANG_ARCH_FLAGS)
 # flags for nvcc
 NVCC_FLAGS = -x=cu --expt-extended-lambda --expt-relaxed-constexpr \
- -arch=$(CUDA_ARCH) -isystem "$(CUDA_DIR)/include"
+ $(NVCC_ARCH_FLAGS) -isystem "$(CUDA_DIR)/include"
 # Prefixes for passing flags to the host compiler and linker when using
 # CUDA_CXX=nvcc
 CUDA_XCOMPILER = -Xcompiler=
@@ -194,6 +207,7 @@ MFEM_USE_BENCHMARK     = NO
 MFEM_USE_PARELAG       = NO
 MFEM_USE_TRIBOL        = NO
 MFEM_USE_ENZYME        = NO
+MFEM_USE_GLVIS         = NO
 
 # Process MFEM_PRECISION -> MFEM_USE_SINGLE, MFEM_USE_DOUBLE
 ifneq ($(filter double Double DOUBLE,$(MFEM_PRECISION)),)
@@ -382,7 +396,7 @@ CUDSS_LIBRARY_DIR = $(CUDSS_DIR)/lib
 CUDSS_OPT         = -I$(CUDSS_INCLUDE_DIR)
 CUDSS_LIB         = \
  $(XLINKER)-rpath,$(CUDSS_LIBRARY_DIR) -L$(CUDSS_LIBRARY_DIR) -lcudss
-# The cuDSS communication and threading libraries. 
+# The cuDSS communication and threading libraries.
 MFEM_CUDSS_COMM_LIB = $(abspath $(wildcard $(or $(CUDSS_COMM_LIB),\
    $(subst @MFEM_DIR@,$(MFEM_DIR), $(CUDSS_LIBRARY_DIR)/libcudss_commlayer_openmpi.so))))
 MFEM_CUDSS_THREADING_LIB = $(abspath $(wildcard $(or $(CUDSS_THREADING_LIB),\
@@ -660,8 +674,23 @@ endif
 ENZYME_OPT = -fplugin=$(ENZYME_PLUGIN)
 ENZYME_LIB =
 
+# GLVis configuration
+include $(GLVIS_MK)
+
 # If YES, enable some informational messages
 VERBOSE = NO
 
 # Optional build tag
 MFEM_BUILD_TAG = $(shell uname -snm)
+
+# Enable -pedantic flag only for gcc or clang. nvcc complains with -pedantic
+# because of line directives.
+PEDANTIC_FLAG = $(if \
+   $(findstring NVIDIA,$(shell $(MFEM_CXX) --version 2>&1)),, \
+   $(if $(or \
+      $(findstring gcc version,$(shell $(MFEM_CXX) -v 2>&1)), \
+      $(findstring clang version,$(shell $(MFEM_CXX) -v 2>&1))),-pedantic,))
+# Enable shadow warnings for clang only; GCC's -Wshadow flags more.
+SHADOW_WARNING_FLAG = $(if $(findstring clang,\
+ $(shell $(MFEM_HOST_CXX) --version 2>/dev/null)),-Wshadow,)
+WARNING_FLAGS = $(PEDANTIC_FLAG) -Wall $(SHADOW_WARNING_FLAG)

config/glvis.mk

@@ -0,0 +1,116 @@
+# GLVis library - Adapted from GLVis' makefile
+
+# Macro that searches for a file in a list of directories returning the first
+# directory that contains the file.
+# $(1) - the file to search for
+# $(2) - list of directories to search
+define find_dir
+$(patsubst %/$(1),%,$(firstword $(wildcard $(foreach d,$(2),$(d)/$(1)))))
+endef
+
+# Macro to find the proper library sub-directory, 'lib64' or 'lib', given a
+# tentative prefix and a library name. Returns empty path if prefix is empty,
+# '/usr', or the library is not found.
+# $(1) - the prefix to search, e.g. $(SDL_DIR)
+# $(2) - library name without 'lib' prefix, e.g. 'SDL2'
+define dir2lib
+$(if $(filter-out /usr,$(1)),$(patsubst %/,%,$(dir $(firstword $(wildcard\
+ $(1)/lib64/lib$(2).* $(1)/lib/lib$(2).*)))))
+endef
+
+BREW_PREFIX := $(if $(NOTMAC),,$(shell brew --prefix 2> /dev/null))
+
+FREETYPE_SEARCH_PATHS = $(BREW_PREFIX) /usr /opt/X11
+FREETYPE_SEARCH_FILE = include/freetype2/ft2build.h
+FREETYPE_DIR = $(call find_dir,$(FREETYPE_SEARCH_FILE),$(FREETYPE_SEARCH_PATHS))
+FREETYPE_LIB_DIR = $(call dir2lib,$(FREETYPE_DIR),freetype)
+FREETYPE_LIBS = -lfreetype -lfontconfig
+
+# If GLEW is in /usr, there's no need to add search paths
+GLEW_SEARCH_PATHS = /usr/local $(BREW_PREFIX) $(abspath ../glew)
+GLEW_SEARCH_FILE = include/GL/glew.h
+GLEW_DIR ?= $(call find_dir,$(GLEW_SEARCH_FILE),$(GLEW_SEARCH_PATHS))
+GLEW_LIB_DIR = $(call dir2lib,$(GLEW_DIR),GLEW)
+GLEW_LIBS = -lGLEW
+
+# If SDL is in /usr, there's no need to add search paths
+SDL_SEARCH_PATHS := /usr/local $(BREW_PREFIX) $(abspath ../SDL2)
+SDL_SEARCH_FILE = include/SDL2/SDL.h
+SDL_DIR ?= $(call find_dir,$(SDL_SEARCH_FILE),$(SDL_SEARCH_PATHS))
+SDL_LIB_DIR = $(call dir2lib,$(SDL_DIR),SDL2)
+SDL_LIBS = -lSDL2
+
+# If GLM is in /usr/include, there's no need to add search paths
+GLM_SEARCH_PATHS = /usr/local/include \
+ $(if $(BREW_PREFIX),$(BREW_PREFIX)/include) $(abspath ../glm)
+GLM_SEARCH_FILE = glm/glm.hpp
+GLM_DIR ?= $(call find_dir,$(GLM_SEARCH_FILE),$(GLM_SEARCH_PATHS))
+
+# If OpenGL is in /usr, there's no need to add search paths
+OPENGL_SEARCH_PATHS = /usr/local /opt/local
+OPENGL_SEARCH_FILE = include/GL/gl.h
+OPENGL_DIR ?= $(call find_dir,$(OPENGL_SEARCH_FILE),$(OPENGL_SEARCH_PATHS))
+OPENGL_LIB_DIR = $(if $(NOTMAC),$(call dir2lib,$(OPENGL_DIR),GL))
+OPENGL_LIBS = $(if $(NOTMAC),-lGL,-framework OpenGL -framework Cocoa)
+
+# Regarding -DGLEW_NO_GLU, see https://github.com/nigels-com/glew/issues/192
+GL_OPTS ?= $(if $(FREETYPE_DIR),-I$(FREETYPE_DIR)/include/freetype2) \
+ $(if $(SDL_DIR),-I$(SDL_DIR)/include) \
+ $(if $(GLEW_DIR),-I$(GLEW_DIR)/include) -DGLEW_NO_GLU \
+ $(if $(GLM_DIR),-I$(GLM_DIR)) \
+ $(if $(OPENGL_DIR),-I$(OPENGL_DIR)/include)
+
+rpath=-Wl,-rpath,
+GL_LIBS ?= $(if $(FREETYPE_LIB_DIR),-L$(FREETYPE_LIB_DIR)) \
+ $(if $(SDL_LIB_DIR),-L$(SDL_LIB_DIR) $(rpath)$(SDL_LIB_DIR)) \
+ $(if $(NOTMAC),$(if $(OPENGL_LIB_DIR),-L$(OPENGL_LIB_DIR) \
+   $(rpath)$(OPENGL_LIB_DIR))) \
+ $(if $(GLEW_LIB_DIR),-L$(GLEW_LIB_DIR) $(rpath)$(GLEW_LIB_DIR)) \
+ $(FREETYPE_LIBS) $(SDL_LIBS) $(GLEW_LIBS)
+
+GLVIS_FLAGS += $(GL_OPTS)
+GLVIS_LIBS  += $(GL_LIBS)
+
+# Take screenshots internally with libtiff, libpng, or sdl2?
+GLVIS_USE_LIBTIFF ?= NO
+GLVIS_USE_LIBPNG  ?= YES
+TIFF_OPTS = -DGLVIS_USE_LIBTIFF -I/sw/include
+TIFF_LIBS = -L/sw/lib -ltiff
+PNG_OPTS = -DGLVIS_USE_LIBPNG
+PNG_LIBS = -lpng
+ifeq ($(GLVIS_USE_LIBTIFF),YES)
+   GLVIS_FLAGS += $(TIFF_OPTS)
+   GLVIS_LIBS  += $(TIFF_LIBS)
+else ifeq ($(GLVIS_USE_LIBPNG),YES)
+   GLVIS_FLAGS += $(PNG_OPTS)
+   GLVIS_LIBS  += $(PNG_LIBS)
+else
+   # no flag --> SDL screenshots
+endif
+
+# EGL headless rendering
+GLVIS_USE_EGL ?= NO
+EGL_OPTS = -DGLVIS_USE_EGL
+EGL_LIBS = -lEGL
+ifeq ($(GLVIS_USE_EGL),YES)
+	GLVIS_FLAGS += $(EGL_OPTS)
+	GLVIS_LIBS  += $(EGL_LIBS)
+endif
+
+# CGL headless rendering
+GLVIS_USE_CGL ?= $(if $(NOTMAC),NO,YES)
+CGL_OPTS = -DGLVIS_USE_CGL
+ifeq ($(GLVIS_USE_CGL),YES)
+	GLVIS_FLAGS += $(CGL_OPTS)
+endif
+
+PTHREAD_LIB = -lpthread
+
+GLVIS_LIBS += $(PTHREAD_LIB)
+GLVIS_LIBS += $(if $(NOTMAC),-lmfem)
+GLVIS_LIBS := $(sort $(GLVIS_LIBS))
+GLVIS_LIBS += $(OPENGL_LIBS)
+
+GLVIS_DIR = @MFEM_DIR@/../glvis
+GLVIS_OPT = 
+GLVIS_LIB = -L$(GLVIS_DIR)/lib -lglvis $(GLVIS_LIBS)

doc/CodeDocumentation.conf.in

@@ -1083,7 +1083,8 @@ EXCLUDE_PATTERNS       =
 # ANamespace::AClass, ANamespace::*Test
 
 EXCLUDE_SYMBOLS        = mfem::internal \
-                         mfem::kernels::internal
+                         mfem::kernels::internal \
+                         mfem::future::detail
 
 # The EXAMPLE_PATH tag can be used to specify one or more files or directories
 # that contain example code fragments that are included (see the \include

examples/CMakeLists.txt

@@ -230,6 +230,11 @@ if (MFEM_USE_GINKGO)
   add_subdirectory(ginkgo)
 endif()
 
+# Include the examples/glvis directory if GLVis is enabled.
+if (MFEM_USE_GLVIS)
+  add_subdirectory(glvis)
+endif()
+
 # Include the examples/hiop directory if HiOp is enabled
 if (MFEM_USE_HIOP)
   add_subdirectory(hiop)

examples/glvis/CMakeLists.txt

@@ -0,0 +1,60 @@
+# Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+set(GLVIS_EXAMPLES_SRCS ex0.cpp)
+if (MFEM_USE_MPI)
+  list(APPEND GLVIS_EXAMPLES_SRCS ex0p.cpp)
+endif()
+
+# Include the source directory where mfem.hpp and mfem-performance.hpp are.
+include_directories(BEFORE ${PROJECT_BINARY_DIR})
+
+# Add "test_glvis" target, see below.
+add_custom_target(test_glvis
+  ${CMAKE_CTEST_COMMAND} -R glvis USES_TERMINAL)
+
+# Add one executable per cpp file, adding "glvis_" as prefix.
+# Sets "test_glvis" as a target that depends on the given examples.
+set(PFX glvis_)
+add_mfem_examples(GLVIS_EXAMPLES_SRCS ${PFX} "" test_glvis)
+
+# Testing.
+# The GLVis tests can be run separately using the target "test_glvis"
+# which builds the examples and runs:
+#   ctest -R glvis
+
+if (MFEM_ENABLE_TESTING)
+  # Command line options for the tests.
+  set(EX0_TEST_OPTS -m ../../data/square-disc.mesh)
+  set(EX0P_TEST_OPTS ${EX0_TEST_OPTS})
+
+  # Add the tests: one test per source file.
+  foreach(SRC_FILE ${GLVIS_EXAMPLES_SRCS})
+    get_filename_component(SRC_FILENAME ${SRC_FILE} NAME)
+    string(REPLACE ".cpp" "" TEST_NAME ${SRC_FILENAME})
+    string(TOUPPER ${TEST_NAME} UP_TEST_NAME)
+    set(TEST_NAME ${PFX}${TEST_NAME})
+
+    set(THIS_TEST_OPTIONS ${${UP_TEST_NAME}_TEST_OPTS})
+    # message(STATUS "Test ${TEST_NAME} options: ${THIS_TEST_OPTIONS}")
+
+    if (NOT (${TEST_NAME} MATCHES ".*p$"))
+      add_test(NAME ${TEST_NAME}_ser
+        COMMAND ${TEST_NAME} ${THIS_TEST_OPTIONS})
+    else()
+      add_test(NAME ${TEST_NAME}_np=${MFEM_MPI_NP}
+        COMMAND ${MPIEXEC} ${MPIEXEC_NUMPROC_FLAG} ${MFEM_MPI_NP}
+        ${MPIEXEC_PREFLAGS}
+        $<TARGET_FILE:${TEST_NAME}> ${THIS_TEST_OPTIONS}
+        ${MPIEXEC_POSTFLAGS})
+    endif()
+  endforeach()
+endif()

examples/glvis/README

@@ -0,0 +1,27 @@
+                    Finite Element Discretization Library
+                                   __
+                       _ __ ___   / _|  ___  _ __ ___
+                      | '_ ` _ \ | |_  / _ \| '_ ` _ \
+                      | | | | | ||  _||  __/| | | | | |
+                      |_| |_| |_||_|   \___||_| |_| |_|
+
+                               https://mfem.org
+
+This directory contains modifications of the example codes that illustrate the
+use of MFEM features based on GLVis in non-server mode.
+
+To build these examples, make sure that MFEM is configured with the option
+"MFEM_USE_GLVIS = YES", see the top-level INSTALL file for details.
+
+Unlike the main examples in examples/, which save files or send data to a
+GLVis server, the codes here use the mfem::glvis_stream class to stream mesh
+and solution data directly into GLVis.
+Currently this directory contains serial and parallel versions of example 0
+(ex0 and ex0p).
+
+We recommend comparing the original example codes with the corresponding files
+in the current directory.
+
+From this directory, the codes can be built with "make" and tested with
+"make test". With CMake, the executables are named glvis_ex0 and glvis_ex0p,
+and the tests can be run with "make test_glvis" or "ctest -R glvis".

examples/glvis/ex0.cpp

@@ -0,0 +1,82 @@
+//                                MFEM Example 0
+//
+// Compile with: make ex0
+//
+// Sample runs:  ex0
+//               ex0 -m ../data/fichera.mesh
+//               ex0 -m ../data/square-disc.mesh -o 2
+//
+// Description: This example code demonstrates the most basic usage of MFEM to
+//              define a simple finite element discretization of the Poisson
+//              problem -Delta u = 1 with zero Dirichlet boundary conditions.
+//              General 2D/3D mesh files and finite element polynomial degrees
+//              can be specified by command line options.
+
+#include "mfem.hpp"
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command line options.
+   string mesh_file = "../data/star.mesh";
+   int order = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order", "Finite element polynomial degree");
+   args.ParseCheck();
+
+   // 2. Read the mesh from the given mesh file, and refine once uniformly.
+   Mesh mesh(mesh_file);
+   mesh.UniformRefinement();
+
+   // 3. Define a finite element space on the mesh. Here we use H1 continuous
+   //    high-order Lagrange finite elements of the given order.
+   H1_FECollection fec(order, mesh.Dimension());
+   FiniteElementSpace fespace(&mesh, &fec);
+   cout << "Number of unknowns: " << fespace.GetTrueVSize() << endl;
+
+   // 4. Extract the list of all the boundary DOFs. These will be marked as
+   //    Dirichlet in order to enforce zero boundary conditions.
+   Array<int> boundary_dofs;
+   fespace.GetBoundaryTrueDofs(boundary_dofs);
+
+   // 5. Define the solution x as a finite element grid function in fespace. Set
+   //    the initial guess to zero, which also sets the boundary conditions.
+   GridFunction x(&fespace);
+   x = 0.0;
+
+   // 6. Set up the linear form b(.) corresponding to the right-hand side.
+   ConstantCoefficient one(1.0);
+   LinearForm b(&fespace);
+   b.AddDomainIntegrator(new DomainLFIntegrator(one));
+   b.Assemble();
+
+   // 7. Set up the bilinear form a(.,.) corresponding to the -Delta operator.
+   BilinearForm a(&fespace);
+   a.AddDomainIntegrator(new DiffusionIntegrator);
+   a.Assemble();
+
+   // 8. Form the linear system A X = B. This includes eliminating boundary
+   //    conditions, applying AMR constraints, and other transformations.
+   SparseMatrix A;
+   Vector B, X;
+   a.FormLinearSystem(boundary_dofs, x, b, A, X, B);
+
+   // 9. Solve the system using PCG with symmetric Gauss-Seidel preconditioner.
+   GSSmoother M(A);
+   PCG(A, M, B, X, 1, 200, 1e-12, 0.0);
+
+   // 10. Recover the solution x as a grid function and save to file.
+   a.RecoverFEMSolution(X, b, x);
+
+   // 11. Send the solution to a non-server mode GLVis.
+   glvis_stream glvis;
+   glvis.precision(8);
+   glvis << "solution\n" << mesh << x << flush;
+
+   return 0;
+}

examples/glvis/ex0p.cpp

@@ -0,0 +1,103 @@
+//                       MFEM Example 0 - Parallel Version
+//
+// Compile with: make ex0p
+//
+// Sample runs:  mpirun -np 4 ex0p
+//               mpirun -np 4 ex0p -m ../data/fichera.mesh
+//               mpirun -np 4 ex0p -m ../data/square-disc.mesh -o 2
+//
+// Description: This example code demonstrates the most basic parallel usage of
+//              MFEM to define a simple finite element discretization of the
+//              Poisson problem -Delta u = 1 with zero Dirichlet boundary
+//              conditions. General 2D/3D serial mesh files and finite element
+//              polynomial degrees can be specified by command line options.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+int main(int argc, char *argv[])
+{
+   // 1. Initialize MPI and HYPRE.
+   Mpi::Init(argc, argv);
+   Hypre::Init();
+
+   // 2. Parse command line options.
+   string mesh_file = "../data/star.mesh";
+   int order = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order", "Finite element polynomial degree");
+   args.ParseCheck();
+
+   // 3. Read the serial mesh from the given mesh file.
+   Mesh serial_mesh(mesh_file);
+
+   // 4. Define a parallel mesh by a partitioning of the serial mesh. Refine
+   //    this mesh once in parallel to increase the resolution.
+   ParMesh mesh(MPI_COMM_WORLD, serial_mesh);
+   serial_mesh.Clear(); // the serial mesh is no longer needed
+   mesh.UniformRefinement();
+
+   // 5. Define a finite element space on the mesh. Here we use H1 continuous
+   //    high-order Lagrange finite elements of the given order.
+   H1_FECollection fec(order, mesh.Dimension());
+   ParFiniteElementSpace fespace(&mesh, &fec);
+   HYPRE_BigInt total_num_dofs = fespace.GlobalTrueVSize();
+   if (Mpi::Root())
+   {
+      cout << "Number of unknowns: " << total_num_dofs << endl;
+   }
+
+   // 6. Extract the list of all the boundary DOFs. These will be marked as
+   //    Dirichlet in order to enforce zero boundary conditions.
+   Array<int> boundary_dofs;
+   fespace.GetBoundaryTrueDofs(boundary_dofs);
+
+   // 7. Define the solution x as a finite element grid function in fespace. Set
+   //    the initial guess to zero, which also sets the boundary conditions.
+   ParGridFunction x(&fespace);
+   x = 0.0;
+
+   // 8. Set up the linear form b(.) corresponding to the right-hand side.
+   ConstantCoefficient one(1.0);
+   ParLinearForm b(&fespace);
+   b.AddDomainIntegrator(new DomainLFIntegrator(one));
+   b.Assemble();
+
+   // 9. Set up the bilinear form a(.,.) corresponding to the -Delta operator.
+   ParBilinearForm a(&fespace);
+   a.AddDomainIntegrator(new DiffusionIntegrator);
+   a.Assemble();
+
+   // 10. Form the linear system A X = B. This includes eliminating boundary
+   //     conditions, applying AMR constraints, parallel assembly, etc.
+   HypreParMatrix A;
+   Vector B, X;
+   a.FormLinearSystem(boundary_dofs, x, b, A, X, B);
+
+   // 11. Solve the system using PCG with hypre's BoomerAMG preconditioner.
+   HypreBoomerAMG M(A);
+   CGSolver cg(MPI_COMM_WORLD);
+   cg.SetRelTol(1e-12);
+   cg.SetMaxIter(2000);
+   cg.SetPrintLevel(1);
+   cg.SetPreconditioner(M);
+   cg.SetOperator(A);
+   cg.Mult(B, X);
+
+   // 12. Recover the solution x as a grid function.
+   a.RecoverFEMSolution(X, b, x);
+
+   // 13. Send the solution to a non-server mode GLVis.
+   glvis_stream glvis;
+   glvis << "parallel " << Mpi::WorldSize() << " " << Mpi::WorldRank() << "\n";
+   glvis.precision(8);
+   glvis << "solution\n" << mesh << x << flush;
+
+   return 0;
+}

examples/glvis/makefile

@@ -0,0 +1,74 @@
+# Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+# at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+# LICENSE and NOTICE for details. LLNL-CODE-806117.
+#
+# This file is part of the MFEM library. For more information and source code
+# availability visit https://mfem.org.
+#
+# MFEM is free software; you can redistribute it and/or modify it under the
+# terms of the BSD-3 license. We welcome feedback and contributions, see file
+# CONTRIBUTING.md for details.
+
+# Use the MFEM build directory
+MFEM_DIR ?= ../..
+MFEM_BUILD_DIR ?= ../..
+MFEM_INSTALL_DIR ?= ../../mfem
+SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/examples/glvis/,)
+CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
+   $(wildcard $(MFEM_INSTALL_DIR)/share/mfem/config.mk))
+
+MFEM_LIB_FILE = mfem_is_not_built
+-include $(CONFIG_MK)
+
+SEQ_EXAMPLES = ex0
+PAR_EXAMPLES = ex0p
+ifeq ($(MFEM_USE_MPI),NO)
+   EXAMPLES = $(SEQ_EXAMPLES)
+else
+   EXAMPLES = $(PAR_EXAMPLES) $(SEQ_EXAMPLES)
+endif
+
+.SUFFIXES:
+.SUFFIXES: .o .cpp .mk
+.PHONY: all clean clean-build clean-exec
+
+# Remove built-in rule
+%: %.cpp
+
+# Replace the default implicit rule for *.cpp files
+%: $(SRC)%.cpp $(MFEM_LIB_FILE) $(CONFIG_MK)
+	$(MFEM_CXX) $(MFEM_FLAGS) $< -o $@ $(MFEM_LIBS)
+
+all: $(EXAMPLES)
+
+ifeq ($(MFEM_USE_GLVIS),NO)
+$(EXAMPLES):
+	$(error MFEM is not configured with GLVIS)
+endif
+
+MFEM_TESTS = EXAMPLES
+include $(MFEM_TEST_MK)
+
+# Testing
+RUN_MPI = $(MFEM_MPIEXEC) $(MFEM_MPIEXEC_NP) $(MFEM_MPI_NP)
+EX0_ARGS := -m ../../data/square-disc.mesh
+
+ex0-test-seq: ex0
+	@$(call mfem-test,$<,, Serial GLVis example,$(EX0_ARGS),SKIP-NO-VIS)
+ex0p-test-par: ex0p
+	@$(call mfem-test,$<, $(RUN_MPI), Parallel GLVis example,$(EX0_ARGS),SKIP-NO-VIS)
+
+# Testing: "test" target and mfem-test* variables are defined in config/test.mk
+
+# Generate an error message if the MFEM library is not built and exit
+$(MFEM_LIB_FILE):
+	$(error The MFEM library is not built)
+
+clean: clean-build clean-exec
+
+clean-build:
+	rm -f *.o *~ $(SEQ_EXAMPLES) $(PAR_EXAMPLES)
+	rm -rf *.dSYM *.TVD.*breakpoints
+
+clean-exec:
+	@:

fem/bilinearform.cpp

@@ -1255,6 +1255,31 @@ void BilinearForm::Mult(const Vector &x, Vector &y) const
    }
 }
 
+void BilinearForm::AddMult(const Vector &x, Vector &y, const real_t a) const
+{
+   if (ext)
+   {
+      ext->AddMult(x, y, a);
+   }
+   else
+   {
+      mat->AddMult(x, y, a);
+   }
+}
+
+void BilinearForm::AddMultTranspose(const Vector &x, Vector &y,
+                                    const real_t a) const
+{
+   if (ext)
+   {
+      ext->AddMultTranspose(x, y, a);
+   }
+   else
+   {
+      mat->AddMultTranspose(x, y, a);
+   }
+}
+
 void BilinearForm::MultTranspose(const Vector & x, Vector & y) const
 {
    if (ext)

fem/bilinearform.hpp

@@ -307,8 +307,8 @@ public:
    { mat->Mult(x, y); mat_e->AddMult(x, y); }
 
    /// Add the matrix vector multiple to a vector:  $ y += a M x $
-   void AddMult(const Vector &x, Vector &y, const real_t a = 1.0) const override
-   { mat -> AddMult (x, y, a); }
+   void AddMult(const Vector &x, Vector &y,
+                const real_t a = 1.0) const override;
 
    /** @brief Add the original uneliminated matrix vector multiple to a vector.
        The original matrix is $ M + Me $ so we have:
@@ -318,8 +318,7 @@ public:
 
    /// Add the matrix transpose vector multiplication:  $ y += a M^T x $
    void AddMultTranspose(const Vector & x, Vector & y,
-                         const real_t a = 1.0) const override
-   { mat->AddMultTranspose(x, y, a); }
+                         const real_t a = 1.0) const override;
 
    /** @brief Add the original uneliminated matrix transpose vector
        multiple to a vector. The original matrix is $ M + M_e $

fem/bilinearform_ext.cpp

@@ -1997,7 +1997,11 @@ void PADiscreteLinearOperatorExtension::Assemble()
    }
    else
    {
-      mfem_error("A real ElementRestriction is required in this setting!");
+      const L2ElementRestriction* l2_elem_restrict =
+         dynamic_cast<const L2ElementRestriction*>(elem_restrict_test);
+      MFEM_VERIFY(l2_elem_restrict,
+                  "A real ElementRestriction is required in this setting!");
+      test_multiplicity = 1.0;
    }
 
    auto tm = test_multiplicity.ReadWrite();
@@ -2036,7 +2040,13 @@ void PADiscreteLinearOperatorExtension::AddMult(
    }
    else
    {
-      mfem_error("In this setting you need a real ElementRestriction!");
+      const L2ElementRestriction* l2_elem_restrict =
+         dynamic_cast<const L2ElementRestriction*>(elem_restrict_test);
+      MFEM_VERIFY(l2_elem_restrict,
+                  "In this setting you need a real ElementRestriction!");
+      tempY.SetSize(y.Size());
+      l2_elem_restrict->MultTranspose(localTest, tempY);
+      y += tempY;
    }
 }
 

fem/coefficient.hpp

@@ -1055,7 +1055,8 @@ public:
 
 typedef VectorCoefficient DiagonalMatrixCoefficient;
 
-/// Base class for Matrix Coefficients that optionally depend on time and space.
+/** Base class for matrix-valued coefficients that optionally depend on time
+    and space. */
 class MatrixCoefficient
 {
 protected:
@@ -1102,6 +1103,9 @@ public:
    /// the quadrature points. The matrix will be transposed or not according to
    /// the boolean argument @a transpose.
    ///
+   /// The stored entries use the same row/column convention as `Eval()`,
+   /// unless `transpose == true`, in which case `K^T` is stored instead.
+   ///
    /// The @a vdim of the QuadratureFunction should be equal to the height times
    /// the width of the matrix.
    virtual void Project(QuadratureFunction &qf, bool transpose=false);

fem/complex_fem.cpp

@@ -588,6 +588,38 @@ SesquilinearForm::AssembleComplexSparseMatrix()
                                   false, false, conv);
 }
 
+void
+SesquilinearForm::BuildComplexOperator(OperatorHandle &A_r,
+                                       OperatorHandle &A_i,
+                                       OperatorHandle &A) const
+{
+   // A = A_r + i A_i
+   A.Clear();
+   if ((!A_r.Ptr() || A_r.Type() == Operator::MFEM_SPARSEMAT) &&
+       (!A_i.Ptr() || A_i.Type() == Operator::MFEM_SPARSEMAT))
+   {
+      ComplexSparseMatrix * A_sp =
+         new ComplexSparseMatrix(A_r.As<SparseMatrix>(),
+                                 A_i.As<SparseMatrix>(),
+                                 A_r.OwnsOperator(),
+                                 A_i.OwnsOperator(),
+                                 conv);
+      A.Reset<ComplexSparseMatrix>(A_sp, true);
+   }
+   else
+   {
+      ComplexOperator * A_op =
+         new ComplexOperator(A_r.Ptr(),
+                             A_i.Ptr(),
+                             A_r.OwnsOperator(),
+                             A_i.OwnsOperator(),
+                             conv);
+      A.Reset<ComplexOperator>(A_op, true);
+   }
+   A_r.SetOperatorOwner(false);
+   A_i.SetOperatorOwner(false);
+}
+
 void
 SesquilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list,
                                    Vector &x, Vector &b,
@@ -716,31 +748,7 @@ SesquilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list,
    B_r.SyncAliasMemory(B);
    B_i.SyncAliasMemory(B);
 
-   // A = A_r + i A_i
-   A.Clear();
-   if ((!A_r.Ptr() || A_r.Type() == Operator::MFEM_SPARSEMAT) &&
-       (!A_i.Ptr() || A_i.Type() == Operator::MFEM_SPARSEMAT))
-   {
-      ComplexSparseMatrix * A_sp =
-         new ComplexSparseMatrix(A_r.As<SparseMatrix>(),
-                                 A_i.As<SparseMatrix>(),
-                                 A_r.OwnsOperator(),
-                                 A_i.OwnsOperator(),
-                                 conv);
-      A.Reset<ComplexSparseMatrix>(A_sp, true);
-   }
-   else
-   {
-      ComplexOperator * A_op =
-         new ComplexOperator(A_r.Ptr(),
-                             A_i.Ptr(),
-                             A_r.OwnsOperator(),
-                             A_i.OwnsOperator(),
-                             conv);
-      A.Reset<ComplexOperator>(A_op, true);
-   }
-   A_r.SetOperatorOwner(false);
-   A_i.SetOperatorOwner(false);
+   BuildComplexOperator(A_r, A_i, A);
 }
 
 void
@@ -777,31 +785,7 @@ SesquilinearForm::FormSystemMatrix(const Array<int> &ess_tdof_list,
       }
    }
 
-   // A = A_r + i A_i
-   A.Clear();
-   if ((!A_r.Ptr() || A_r.Type() == Operator::MFEM_SPARSEMAT) &&
-       (!A_i.Ptr() || A_i.Type() == Operator::MFEM_SPARSEMAT))
-   {
-      ComplexSparseMatrix * A_sp =
-         new ComplexSparseMatrix(A_r.As<SparseMatrix>(),
-                                 A_i.As<SparseMatrix>(),
-                                 A_r.OwnsOperator(),
-                                 A_i.OwnsOperator(),
-                                 conv);
-      A.Reset<ComplexSparseMatrix>(A_sp, true);
-   }
-   else
-   {
-      ComplexOperator * A_op =
-         new ComplexOperator(A_r.Ptr(),
-                             A_i.Ptr(),
-                             A_r.OwnsOperator(),
-                             A_i.OwnsOperator(),
-                             conv);
-      A.Reset<ComplexOperator>(A_op, true);
-   }
-   A_r.SetOperatorOwner(false);
-   A_i.SetOperatorOwner(false);
+   BuildComplexOperator(A_r, A_i, A);
 }
 
 void
@@ -1893,6 +1877,81 @@ ParSesquilinearForm::ParallelAssemble()
                                     true, true, conv);
 }
 
+void
+ParSesquilinearForm::BuildComplexOperator(OperatorHandle &A_r,
+                                          OperatorHandle &A_i,
+                                          OperatorHandle &A) const
+{
+   // A = A_r + i A_i
+   A.Clear();
+   if ((!A_r.Ptr() || A_r.Type() == Operator::Hypre_ParCSR) &&
+       (!A_i.Ptr() || A_i.Type() == Operator::Hypre_ParCSR))
+   {
+      ComplexHypreParMatrix * A_hyp =
+         new ComplexHypreParMatrix(A_r.As<HypreParMatrix>(),
+                                   A_i.As<HypreParMatrix>(),
+                                   A_r.OwnsOperator(),
+                                   A_i.OwnsOperator(),
+                                   conv);
+      A.Reset<ComplexHypreParMatrix>(A_hyp, true);
+   }
+   else
+   {
+      ComplexOperator * A_op =
+         new ComplexOperator(A_r.As<Operator>(),
+                             A_i.As<Operator>(),
+                             A_r.OwnsOperator(),
+                             A_i.OwnsOperator(),
+                             conv);
+      A.Reset<ComplexOperator>(A_op, true);
+   }
+   A_r.SetOperatorOwner(false);
+   A_i.SetOperatorOwner(false);
+}
+
+namespace
+{
+struct ZeroDiagonalHypreKernel
+{
+   const int *ess_tdof_list;
+   const HYPRE_Int *diag_i;
+   real_t *diag_data;
+
+   void MFEM_HOST_DEVICE operator()(int k) const
+   {
+      const int j = ess_tdof_list[k];
+      diag_data[diag_i[j]] = 0.0;
+   }
+};
+}
+
+void
+ParSesquilinearForm::SetImaginaryEssentialDiagonalToZero(
+   const Array<int> &ess_tdof_list, OperatorHandle &A)
+{
+   if (A.Type() == Operator::Hypre_ParCSR)
+   {
+      const int n = ess_tdof_list.Size();
+      HypreParMatrix *Ah;
+      A.Get(Ah);
+      hypre_ParCSRMatrix *Aih = *Ah;
+      Ah->HypreReadWrite();
+      const int *d_ess_tdof_list =
+         ess_tdof_list.GetMemory().Read(GetHypreForallMemoryClass(), n);
+      HYPRE_Int *d_diag_i = Aih->diag->i;
+      real_t *d_diag_data = Aih->diag->data;
+      mfem::hypre_forall(n, ZeroDiagonalHypreKernel
+      {
+         d_ess_tdof_list, d_diag_i, d_diag_data
+      });
+   }
+   else
+   {
+      A.As<ConstrainedOperator>()->SetDiagonalPolicy
+      (mfem::Operator::DiagonalPolicy::DIAG_ZERO);
+   }
+}
+
 void
 ParSesquilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list,
                                       Vector &x, Vector &b,
@@ -1993,27 +2052,7 @@ ParSesquilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list,
       });
       // Modify off-diagonal blocks (imaginary parts of the matrix) to conform
       // with standard essential BC treatment
-      if (A_i.Type() == Operator::Hypre_ParCSR)
-      {
-         HypreParMatrix * Ah;
-         A_i.Get(Ah);
-         hypre_ParCSRMatrix *Aih = *Ah;
-         Ah->HypreReadWrite();
-         const int *d_ess_tdof_list =
-            ess_tdof_list.GetMemory().Read(GetHypreForallMemoryClass(), n);
-         HYPRE_Int *d_diag_i = Aih->diag->i;
-         real_t *d_diag_data = Aih->diag->data;
-         mfem::hypre_forall(n, [=] MFEM_HOST_DEVICE (int k)
-         {
-            const int j = d_ess_tdof_list[k];
-            d_diag_data[d_diag_i[j]] = 0.0;
-         });
-      }
-      else
-      {
-         A_i.As<ConstrainedOperator>()->SetDiagonalPolicy
-         (mfem::Operator::DiagonalPolicy::DIAG_ZERO);
-      }
+      SetImaginaryEssentialDiagonalToZero(ess_tdof_list, A_i);
    }
 
    if (conv == ComplexOperator::BLOCK_SYMMETRIC)
@@ -2032,31 +2071,7 @@ ParSesquilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list,
    B_r.SyncAliasMemory(B);
    B_i.SyncAliasMemory(B);
 
-   // A = A_r + i A_i
-   A.Clear();
-   if ((!A_r.Ptr() || A_r.Type() == Operator::Hypre_ParCSR) &&
-       (!A_i.Ptr() || A_i.Type() == Operator::Hypre_ParCSR))
-   {
-      ComplexHypreParMatrix * A_hyp =
-         new ComplexHypreParMatrix(A_r.As<HypreParMatrix>(),
-                                   A_i.As<HypreParMatrix>(),
-                                   A_r.OwnsOperator(),
-                                   A_i.OwnsOperator(),
-                                   conv);
-      A.Reset<ComplexHypreParMatrix>(A_hyp, true);
-   }
-   else
-   {
-      ComplexOperator * A_op =
-         new ComplexOperator(A_r.As<Operator>(),
-                             A_i.As<Operator>(),
-                             A_r.OwnsOperator(),
-                             A_i.OwnsOperator(),
-                             conv);
-      A.Reset<ComplexOperator>(A_op, true);
-   }
-   A_r.SetOperatorOwner(false);
-   A_i.SetOperatorOwner(false);
+   BuildComplexOperator(A_r, A_i, A);
 }
 
 void
@@ -2081,50 +2096,10 @@ ParSesquilinearForm::FormSystemMatrix(const Array<int> &ess_tdof_list,
    {
       // Modify off-diagonal blocks (imaginary parts of the matrix) to conform
       // with standard essential BC treatment
-      if ( A_i.Type() == Operator::Hypre_ParCSR )
-      {
-         int n = ess_tdof_list.Size();
-         HypreParMatrix * Ah;
-         A_i.Get(Ah);
-         hypre_ParCSRMatrix * Aih = *Ah;
-         for (int k = 0; k < n; k++)
-         {
-            int j = ess_tdof_list[k];
-            Aih->diag->data[Aih->diag->i[j]] = 0.0;
-         }
-      }
-      else
-      {
-         A_i.As<ConstrainedOperator>()->SetDiagonalPolicy
-         (mfem::Operator::DiagonalPolicy::DIAG_ZERO);
-      }
+      SetImaginaryEssentialDiagonalToZero(ess_tdof_list, A_i);
    }
 
-   // A = A_r + i A_i
-   A.Clear();
-   if ((!A_r.Ptr() || A_r.Type() == Operator::Hypre_ParCSR) &&
-       (!A_i.Ptr() || A_i.Type() == Operator::Hypre_ParCSR))
-   {
-      ComplexHypreParMatrix * A_hyp =
-         new ComplexHypreParMatrix(A_r.As<HypreParMatrix>(),
-                                   A_i.As<HypreParMatrix>(),
-                                   A_r.OwnsOperator(),
-                                   A_i.OwnsOperator(),
-                                   conv);
-      A.Reset<ComplexHypreParMatrix>(A_hyp, true);
-   }
-   else
-   {
-      ComplexOperator * A_op =
-         new ComplexOperator(A_r.As<Operator>(),
-                             A_i.As<Operator>(),
-                             A_r.OwnsOperator(),
-                             A_i.OwnsOperator(),
-                             conv);
-      A.Reset<ComplexOperator>(A_op, true);
-   }
-   A_r.SetOperatorOwner(false);
-   A_i.SetOperatorOwner(false);
+   BuildComplexOperator(A_r, A_i, A);
 }
 
 void

fem/complex_fem.hpp

@@ -392,6 +392,9 @@ private:
    bool RealInteg();
    bool ImagInteg();
 
+   void BuildComplexOperator(OperatorHandle &A_r, OperatorHandle &A_i,
+                             OperatorHandle &A) const;
+
 public:
    SesquilinearForm(FiniteElementSpace *fes,
                     ComplexOperator::Convention
@@ -986,6 +989,12 @@ private:
    bool RealInteg();
    bool ImagInteg();
 
+   void SetImaginaryEssentialDiagonalToZero(
+      const Array<int> &ess_tdof_list, OperatorHandle &A);
+
+   void BuildComplexOperator(OperatorHandle &A_r, OperatorHandle &A_i,
+                             OperatorHandle &A) const;
+
 public:
    ParSesquilinearForm(ParFiniteElementSpace *pf,
                        ComplexOperator::Convention

fem/datacollection.cpp

@@ -38,9 +38,24 @@ int DataCollection::create_directory(const std::string &dir_name,
    // create directories recursively
    const char path_delim = '/';
    std::string::size_type pos = 0;
-   int err_flag;
+   int err_flag = 0;
 #ifdef MFEM_USE_MPI
    const ParMesh *pmesh = dynamic_cast<const ParMesh*>(mesh);
+   // In addition to the global root, let the lowest rank on each shared-memory
+   // node create the directory too, so that node-local (non-shared) filesystems
+   // get it on every node rather than only where the global root lives. On a
+   // shared filesystem the extra mkdir() hits EEXIST and is tolerated below.
+   bool node_root = true;
+   if (pmesh)
+   {
+      MPI_Comm node_comm;
+      MPI_Comm_split_type(pmesh->GetComm(), MPI_COMM_TYPE_SHARED, myid,
+                          MPI_INFO_NULL, &node_comm);
+      int node_rank;
+      MPI_Comm_rank(node_comm, &node_rank);
+      node_root = (node_rank == 0);
+      MPI_Comm_free(&node_comm);
+   }
 #endif
 
    do
@@ -52,7 +67,7 @@ int DataCollection::create_directory(const std::string &dir_name,
       err_flag = mkdir(subdir.c_str(), 0777);
       err_flag = (err_flag && (errno != EEXIST)) ? 1 : 0;
 #else
-      if (myid == 0 || pmesh == NULL)
+      if (node_root || pmesh == NULL)
       {
          err_flag = mkdir(subdir.c_str(), 0777);
          err_flag = (err_flag && (errno != EEXIST)) ? 1 : 0;
@@ -64,7 +79,8 @@ int DataCollection::create_directory(const std::string &dir_name,
 #ifdef MFEM_USE_MPI
    if (pmesh)
    {
-      MPI_Bcast(&err_flag, 1, MPI_INT, 0, pmesh->GetComm());
+      MPI_Allreduce(MPI_IN_PLACE, &err_flag, 1, MPI_INT, MPI_MAX,
+                    pmesh->GetComm());
    }
 #endif
 

fem/dfem/doperator.cpp

@@ -51,4 +51,52 @@ DifferentiableOperator::DifferentiableOperator(
    }
 }
 
+
+void FDJacobian::Mult(const Vector &v, Vector &y) const
+{
+   // See [1] for choice of eps.
+   //
+   // [1] Woodward, C.S., Gardner, D.J. and Evans, K.J., 2015. On the use of
+   // finite difference matrix-vector products in Newton-Krylov solvers for
+   // implicit climate dynamics with spectral elements. Procedia Computer
+   // Science, 51, pp.2036-2045.
+   real_t eps;
+   if (fixed_eps > 0.0)
+   {
+      eps = fixed_eps;
+   }
+   else
+   {
+      const real_t vnorm_local = v.Norml2();
+      real_t vnorm;
+      MPI_Allreduce(&vnorm_local, &vnorm, 1, MPITypeMap<real_t>::mpi_type, MPI_SUM,
+                    MPI_COMM_WORLD);
+      eps = lambda * (lambda + xnorm / vnorm);
+   }
+
+   // x + eps * v
+   {
+      const auto d_v = v.Read();
+      const auto d_x = x.Read();
+      auto d_xpev = xpev.Write();
+      mfem::forall(x.Size(), [=] MFEM_HOST_DEVICE (int i)
+      {
+         d_xpev[i] = d_x[i] + eps * d_v[i];
+      });
+   }
+
+   // y = f(x + eps * v)
+   op.Mult(xpev, y);
+
+   // y = (f(x + eps * v) - f(x)) / eps
+   {
+      const auto d_f = f.Read();
+      auto d_y = y.ReadWrite();
+      mfem::forall(f.Size(), [=] MFEM_HOST_DEVICE (int i)
+      {
+         d_y[i] = (d_y[i] - d_f[i]) / eps;
+      });
+   }
+}
+
 #endif // MFEM_USE_MPI

fem/dfem/doperator.hpp

@@ -697,17 +697,18 @@ void DifferentiableOperator::AddIntegrator(
 
    // The explicit captures are necessary to avoid dependency on
    // the specific instance of this class (this pointer).
-   restriction_callback =
-      [=, solutions = this->solutions, parameters = this->parameters]
-      (std::vector<Vector> &sol,
-       const std::vector<Vector> &par,
-       std::vector<Vector> &f)
+   restriction_callback = [element_dof_ordering,
+                           solutions_ = this->solutions,
+                           parameters_ = this->parameters]
+   (std::vector<Vector> &sol,
+    const std::vector<Vector> &par,
+    std::vector<Vector> &f)
    {
-      restriction<entity_t>(solutions, sol, f,
+      restriction<entity_t>(solutions_, sol, f,
                             element_dof_ordering);
-      restriction<entity_t>(parameters, par, f,
+      restriction<entity_t>(parameters_, par, f,
                             element_dof_ordering,
-                            solutions.size());
+                            solutions_.size());
    };
 
    prolongation_transpose = get_prolongation_transpose(
@@ -835,19 +836,19 @@ void DifferentiableOperator::AddIntegrator(
 
          // capture by ref:
          &restriction_cb = this->restriction_callback,
-         &fields_e = this->fields_e,
-         &residual_e = this->residual_e,
-         &output_restriction_transpose = this->output_restriction_transpose
+         &fields_e_ = this->fields_e,
+         &residual_e_ = this->residual_e,
+         &output_restriction_transpose_ = this->output_restriction_transpose
       ]
       (std::vector<Vector> &sol, const std::vector<Vector> &par, Vector &res)
       mutable // mutable: needed to modify 'shmem_cache'
    {
-      restriction_cb(sol, par, fields_e);
+      restriction_cb(sol, par, fields_e_);
 
-      residual_e = 0.0;
-      auto ye = Reshape(residual_e.ReadWrite(), test_vdim, num_test_dof, num_entities);
+      residual_e_ = 0.0;
+      auto ye = Reshape(residual_e_.ReadWrite(), test_vdim, num_test_dof, num_entities);
 
-      auto wrapped_fields_e = wrap_fields(fields_e,
+      auto wrapped_fields_e = wrap_fields(fields_e_,
                                           action_shmem_info.field_sizes,
                                           num_entities);
 
@@ -878,7 +879,7 @@ void DifferentiableOperator::AddIntegrator(
             y, fhat, output_fop, output_dtq_shmem[0],
             scratch_shmem, dimension, use_sum_factorization);
       }, num_entities, thread_blocks, action_shmem_info.total_size, shmem_cache.ReadWrite());
-      output_restriction_transpose(residual_e, res);
+      output_restriction_transpose_(residual_e_, res);
    });
 
    // Without this compile-time check, some valid instantiations of this method
@@ -1193,7 +1194,7 @@ void DifferentiableOperator::AddIntegrator(
 
                // capture by ref:
                &qpdc_mem = derivative_qp_caches_ref,
-               &fields = fields_ref
+               &fields_ = fields_ref
             ](std::vector<Vector> &f_e, SparseMatrix *&A) mutable
          {
             auto wrapped_fields_e = wrap_fields(f_e, shmem_info.field_sizes,
@@ -1241,14 +1242,14 @@ void DifferentiableOperator::AddIntegrator(
             {
                if (input_is_dependent[s])
                {
-                  trial_field = &fields[input_to_field[s]];
+                  trial_field = &fields_[input_to_field[s]];
                }
             }
 
             auto trial_fes = *std::get_if<const ParFiniteElementSpace *>
                              (&trial_field->data);
             auto test_fes = *std::get_if<const ParFiniteElementSpace *>
-                            (&fields[output_to_field[0]].data);
+                            (&fields_[output_to_field[0]].data);
 
             A = new SparseMatrix(test_fes->GetVSize(), trial_fes->GetVSize());
 
@@ -1334,7 +1335,7 @@ void DifferentiableOperator::AddIntegrator(
                input_to_field,
                output_to_field,
                &spmatcb = assemble_derivative_sparsematrix_callbacks_ref,
-               &fields = fields_ref
+               &fields_ = fields_ref
             ](std::vector<Vector> &f_e, HypreParMatrix *&A) mutable
          {
             SparseMatrix *spmat = nullptr;
@@ -1366,14 +1367,14 @@ void DifferentiableOperator::AddIntegrator(
             {
                if (input_is_dependent[s])
                {
-                  trial_field = &fields[input_to_field[s]];
+                  trial_field = &fields_[input_to_field[s]];
                }
             }
 
             auto trial_fes = *std::get_if<const ParFiniteElementSpace *>
                              (&trial_field->data);
             auto test_fes = *std::get_if<const ParFiniteElementSpace *>
-                            (&fields[output_to_field[0]].data);
+                            (&fields_[output_to_field[0]].data);
 
             if (same_test_and_trial)
             {

fem/dfem/util.hpp

@@ -597,7 +597,7 @@ struct ThreadBlocks
    int z = 1;
 };
 
-#if defined(MFEM_USE_CUDA_OR_HIP)
+#if defined(MFEM_USE_CUDA_OR_HIP_LANG)
 template <typename func_t>
 __global__ void forall_kernel_shmem(func_t f, int n)
 {
@@ -617,10 +617,11 @@ void forall(func_t f,
             int num_shmem = 0,
             real_t *shmem = nullptr)
 {
-   if (Device::Allows(Backend::CUDA_MASK) ||
-       Device::Allows(Backend::HIP_MASK))
+   internal::RequireKernelCompilation();
+
+#if defined(MFEM_USE_CUDA_OR_HIP_LANG)
+   if (Device::Allows(Backend::CUDA_MASK | Backend::HIP_MASK))
    {
-#if defined(MFEM_USE_CUDA_OR_HIP)
       // int gridsize = (N + Z - 1) / Z;
       int num_bytes = num_shmem * sizeof(decltype(shmem));
       dim3 block_size(blocks.x, blocks.y, blocks.z);
@@ -631,9 +632,10 @@ void forall(func_t f,
       MFEM_GPU_CHECK(hipGetLastError());
 #endif
       MFEM_DEVICE_SYNC;
-#endif
+      return;
    }
-   else if (Device::Allows(Backend::CPU_MASK))
+#endif
+   if (Device::Allows(Backend::CPU_MASK))
    {
       MFEM_ASSERT(!((bool)num_shmem != (bool)shmem),
                   "Backend::CPU needs a pre-allocated shared memory block");
@@ -671,52 +673,7 @@ public:
                     MPI_COMM_WORLD);
    }
 
-   void Mult(const Vector &v, Vector &y) const override
-   {
-      // See [1] for choice of eps.
-      //
-      // [1] Woodward, C.S., Gardner, D.J. and Evans, K.J., 2015. On the use of
-      // finite difference matrix-vector products in Newton-Krylov solvers for
-      // implicit climate dynamics with spectral elements. Procedia Computer
-      // Science, 51, pp.2036-2045.
-      real_t eps;
-      if (fixed_eps > 0.0)
-      {
-         eps = fixed_eps;
-      }
-      else
-      {
-         const real_t vnorm_local = v.Norml2();
-         real_t vnorm;
-         MPI_Allreduce(&vnorm_local, &vnorm, 1, MPITypeMap<real_t>::mpi_type, MPI_SUM,
-                       MPI_COMM_WORLD);
-         eps = lambda * (lambda + xnorm / vnorm);
-      }
-
-      // x + eps * v
-      {
-         const auto d_v = v.Read();
-         const auto d_x = x.Read();
-         auto d_xpev = xpev.Write();
-         mfem::forall(x.Size(), [=] MFEM_HOST_DEVICE (int i)
-         {
-            d_xpev[i] = d_x[i] + eps * d_v[i];
-         });
-      }
-
-      // y = f(x + eps * v)
-      op.Mult(xpev, y);
-
-      // y = (f(x + eps * v) - f(x)) / eps
-      {
-         const auto d_f = f.Read();
-         auto d_y = y.ReadWrite();
-         mfem::forall(f.Size(), [=] MFEM_HOST_DEVICE (int i)
-         {
-            d_y[i] = (d_y[i] - d_f[i]) / eps;
-         });
-      }
-   }
+   void Mult(const Vector &v, Vector &y) const override;
 
    virtual MemoryClass GetMemoryClass() const override
    {

fem/fe/fe_base.cpp

@@ -1316,13 +1316,14 @@ void VectorFiniteElement::Project_RT(
    }
 }
 
-void VectorFiniteElement::ProjectGrad_RT(
+void VectorFiniteElement::ProjectCurl2D_RT(
    const real_t *nk, const Array<int> &d2n, const FiniteElement &fe,
    ElementTransformation &Trans, DenseMatrix &grad) const
 {
+   // 2D "ProjectCurl_RT"
    if (dim != 2)
    {
-      mfem_error("VectorFiniteElement::ProjectGrad_RT works only in 2D!");
+      mfem_error("VectorFiniteElement::ProjectCurl2D_RT works only in 2D!");
    }
 
    DenseMatrix dshape(fe.GetDof(), fe.GetDim());
@@ -1333,8 +1334,8 @@ void VectorFiniteElement::ProjectGrad_RT(
    for (int k = 0; k < dof; k++)
    {
       fe.CalcDShape(Nodes.IntPoint(k), dshape);
-      tk[0] = nk[d2n[k]*dim+1];
-      tk[1] = -nk[d2n[k]*dim];
+      tk[0] = -nk[d2n[k]*dim+1];
+      tk[1] = nk[d2n[k]*dim];
       dshape.Mult(tk, grad_k);
       for (int j = 0; j < grad_k.Size(); j++)
       {
@@ -1381,7 +1382,7 @@ void VectorFiniteElement::ProjectCurl_ND(
    }
 }
 
-void VectorFiniteElement::ProjectCurl_RT(
+void VectorFiniteElement::ProjectCurl3D_RT(
    const real_t *nk, const Array<int> &d2n, const FiniteElement &fe,
    ElementTransformation &Trans, DenseMatrix &curl) const
 {

fem/fe/fe_base.hpp

@@ -957,10 +957,11 @@ protected:
                    const FiniteElement &fe, ElementTransformation &Trans,
                    DenseMatrix &I) const;
 
-   // rotated gradient in 2D
-   void ProjectGrad_RT(const real_t *nk, const Array<int> &d2n,
-                       const FiniteElement &fe, ElementTransformation &Trans,
-                       DenseMatrix &grad) const;
+   // Input is a scalar representing the Z (out of plane) component, Output is
+   // the X-Y (in-plane) RT curl
+   void ProjectCurl2D_RT(const real_t *nk, const Array<int> &d2n,
+                         const FiniteElement &fe, ElementTransformation &Trans,
+                         DenseMatrix &grad) const;
 
    // Compute the curl as a discrete operator from ND FE (fe) to ND FE (this).
    // The natural FE for the range is RT, so this is an approximation.
@@ -968,9 +969,9 @@ protected:
                        const FiniteElement &fe, ElementTransformation &Trans,
                        DenseMatrix &curl) const;
 
-   void ProjectCurl_RT(const real_t *nk, const Array<int> &d2n,
-                       const FiniteElement &fe, ElementTransformation &Trans,
-                       DenseMatrix &curl) const;
+   void ProjectCurl3D_RT(const real_t *nk, const Array<int> &d2n,
+                         const FiniteElement &fe, ElementTransformation &Trans,
+                         DenseMatrix &curl) const;
 
    /** @brief Project a vector coefficient onto the ND basis functions
        @param tk    Edge tangent vectors for this element type
@@ -1446,6 +1447,8 @@ public:
                                 dof2quad_array_open);
    }
 
+   const Poly_1D::Basis &GetOpenBasis1D() const { return obasis1d; }
+
    virtual ~VectorTensorFiniteElement();
 };
 

fem/fe/fe_rt.hpp

@@ -73,16 +73,11 @@ public:
    void Project(const FiniteElement &fe, ElementTransformation &Trans,
                 DenseMatrix &I) const override
    { Project_RT(nk, dof2nk, fe, Trans, I); }
-   // Gradient + rotation = Curl: H1 -> H(div)
-   void ProjectGrad(const FiniteElement &fe,
-                    ElementTransformation &Trans,
-                    DenseMatrix &grad) const override
-   { ProjectGrad_RT(nk, dof2nk, fe, Trans, grad); }
    // Curl = Gradient + rotation: H1 -> H(div)
    void ProjectCurl(const FiniteElement &fe,
                     ElementTransformation &Trans,
                     DenseMatrix &curl) const override
-   { ProjectGrad_RT(nk, dof2nk, fe, Trans, curl); }
+   { ProjectCurl2D_RT(nk, dof2nk, fe, Trans, curl); }
 
    void GetFaceMap(const int face_id, Array<int> &face_map) const override;
 
@@ -148,7 +143,7 @@ public:
    void ProjectCurl(const FiniteElement &fe,
                     ElementTransformation &Trans,
                     DenseMatrix &curl) const override
-   { ProjectCurl_RT(nk, dof2nk, fe, Trans, curl); }
+   { ProjectCurl3D_RT(nk, dof2nk, fe, Trans, curl); }
 
    /// @brief Return the mapping from lexicographically ordered face DOFs to
    /// lexicographically ordered element DOFs corresponding to local face
@@ -210,16 +205,11 @@ public:
    void Project(const FiniteElement &fe, ElementTransformation &Trans,
                 DenseMatrix &I) const override
    { Project_RT(nk, dof2nk, fe, Trans, I); }
-   // Gradient + rotation = Curl: H1 -> H(div)
-   void ProjectGrad(const FiniteElement &fe,
-                    ElementTransformation &Trans,
-                    DenseMatrix &grad) const override
-   { ProjectGrad_RT(nk, dof2nk, fe, Trans, grad); }
    // Curl = Gradient + rotation: H1 -> H(div)
    void ProjectCurl(const FiniteElement &fe,
                     ElementTransformation &Trans,
                     DenseMatrix &curl) const override
-   { ProjectGrad_RT(nk, dof2nk, fe, Trans, curl); }
+   { ProjectCurl2D_RT(nk, dof2nk, fe, Trans, curl); }
 };
 
 
@@ -274,7 +264,7 @@ public:
    void ProjectCurl(const FiniteElement &fe,
                     ElementTransformation &Trans,
                     DenseMatrix &curl) const override
-   { ProjectCurl_RT(nk, dof2nk, fe, Trans, curl); }
+   { ProjectCurl3D_RT(nk, dof2nk, fe, Trans, curl); }
 };
 
 class RT_WedgeElement : public VectorFiniteElement
@@ -332,7 +322,7 @@ public:
    void ProjectCurl(const FiniteElement &fe,
                     ElementTransformation &Trans,
                     DenseMatrix &curl) const override
-   { ProjectCurl_RT(nk, dof2nk, fe, Trans, curl); }
+   { ProjectCurl3D_RT(nk, dof2nk, fe, Trans, curl); }
 };
 
 /** Arbitrary order H(Div) basis functions defined on pyramid-shaped elements
@@ -428,7 +418,7 @@ public:
    virtual void ProjectCurl(const FiniteElement &fe,
                             ElementTransformation &Trans,
                             DenseMatrix &curl) const
-   { ProjectCurl_RT(nk, dof2nk, fe, Trans, curl); }
+   { ProjectCurl3D_RT(nk, dof2nk, fe, Trans, curl); }
 
    void CalcRawVShape(const IntegrationPoint &ip,
                       DenseMatrix &shape) const;

fem/gslib.cpp

@@ -556,7 +556,7 @@ void obboxsurf_calc_3(Vector &bb,
    gslib::lagrange_fun *const lag = gslib::gll_lag_setup(work, n);
    lag(I0, work, n, 1, 0);
 
-   for (int ie = 0; ie < nel; ie++,x+=n2,y+=n2,z+=n2)
+   for (int ie = 0; (unsigned)ie < nel; ie++,x+=n2,y+=n2,z+=n2)
    {
       struct gslib::dbl_range ab[3];
       struct gslib::dbl_range tb[3];
@@ -780,7 +780,7 @@ void obboxedge_calc_2(Vector &bb,
    gslib::lagrange_fun *const lag = gslib::gll_lag_setup(work, nr);
    lag(I0r, work, nr,1, 0);
 
-   for (int ie = 0; ie < nel; ie++,x+=nr,y+=nr)
+   for (int ie = 0; (unsigned)ie < nel; ie++,x+=nr,y+=nr)
    {
       double x0[2], A[4];
       struct gslib::dbl_range ab[2], tb[2];
@@ -892,7 +892,7 @@ void obboxedge_calc_3(Vector &bb,
    gslib::lagrange_fun *const lag = gslib::gll_lag_setup(work, nr);
    lag(I0r, work, nr, 1, 0);
 
-   for (int ie = 0; ie < nel; ie++,x+=nr,y+=nr,z+=nr)
+   for (int ie = 0; (unsigned)ie < nel; ie++,x+=nr,y+=nr,z+=nr)
    {
       double x0[3], A[9], Ai[9];
       struct gslib::dbl_range ab[3], tb[3];
@@ -4518,7 +4518,7 @@ Mesh* FindPointsGSLIB::GetBoundingBoxMesh(int type)
    int eidx = 0;
    if (myid == save_rank)
    {
-      for (int p = 0; p < gsl_comm->np; p++)
+      for (int p = 0; (unsigned)p < gsl_comm->np; p++)
       {
          if (static_cast<unsigned int>(p) != save_rank)
          {

fem/integ/bilininteg_convection_pa.cpp

@@ -178,6 +178,8 @@ void ConvectionIntegrator::AssemblePA(const FiniteElementSpace &fes)
    // Assumes tensor-product elements
    Mesh *mesh = fes.GetMesh();
    const FiniteElement &el = *fes.GetTypicalFE();
+   MFEM_VERIFY(el.GetMapType() == FiniteElement::VALUE,
+               "Only value map type currently supported");
    ElementTransformation &Trans = *mesh->GetTypicalElementTransformation();
    const IntegrationRule *ir = IntRule ? IntRule : &GetRule(el, Trans);
    if (DeviceCanUseCeed())

fem/integ/bilininteg_hcurl_kernels.cpp

@@ -785,6 +785,23 @@ void PAHcurlL2Setup2D(const int Q1D,
    });
 }
 
+void PAHcurlL2IntSetup2D(const int Q1D, const int NE, const Array<real_t> &w,
+                         Vector &coeff, const Vector &detJ, Vector &op)
+{
+   const int NQ = Q1D*Q1D;
+   auto W = w.Read();
+   auto C = Reshape(coeff.Read(), NQ, NE);
+   auto J = Reshape(detJ.Read(), NQ, NE);
+   auto y = Reshape(op.Write(), NQ, NE);
+   mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
+   {
+      for (int q = 0; q < NQ; ++q)
+      {
+         y(q,e) = W[q] * C(q,e) / J(q,e);
+      }
+   });
+}
+
 void PAHcurlL2Setup3D(const int NQ,
                       const int coeffDim,
                       const int NE,

fem/integ/bilininteg_hcurl_kernels.hpp

@@ -1889,13 +1889,17 @@ inline void SmemPACurlCurlApply3D(const int d1d,
    ForallWrap<3>(true, NE, device_kernel, host_kernel, Q1D, Q1D, Q1D);
 }
 
-// PA H(curl)-L2 Assemble 2D kernel
+// PA H(curl)-L2 value Assemble 2D kernel
 void PAHcurlL2Setup2D(const int Q1D,
                       const int NE,
                       const Array<real_t> &w,
                       Vector &coeff,
                       Vector &op);
 
+// PA H(curl)-L2 integral Assemble 2D kernel
+void PAHcurlL2IntSetup2D(const int Q1D, const int NE, const Array<real_t> &w,
+                         Vector &coeff, const Vector &detJ, Vector &op);
+
 // PA H(curl)-L2 Assemble 3D kernel
 void PAHcurlL2Setup3D(const int NQ,
                       const int coeffDim,

fem/integ/bilininteg_hcurlhdiv_kernels.hpp

@@ -864,8 +864,656 @@ inline void PAHcurlHdivApplyTranspose3D(const int d1d,
    }); // end of element loop
 }
 
+namespace curlinterp
+{
+constexpr int NBZ3D(int ndof_o, int nquad_o, int mdq)
+{
+   if (ndof_o <= 0 || nquad_o <= 0)
+   {
+      return 1;
+   }
+   int ndof_c = ndof_o + 1;
+   int nquad_c = nquad_o + 1;
+   // z dimension is capped at 64 on nvidia and amd gpus
+   int tmp =
+      std::min((128 + mdq * mdq * (mdq - 1) - 1) / (mdq * mdq * (mdq - 1)), 64);
+   int smem_req =
+      sizeof(mfem::real_t) *
+      ((3 * ndof_c * ndof_c * ndof_o + 2 * 2 * mdq * mdq * mdq) * tmp +
+       ndof_c * nquad_o + ndof_c * nquad_c + ndof_o * nquad_o);
+   // assume GPU has at least 48k shared memory
+   return std::max(std::min(tmp, (48 * 1024 + smem_req - 1) / smem_req), 1);
+}
+}
+
+template <int T_NDOF_O, int T_NQUAD_O>
+void CurlInterpolatorApply3DSmem(const int ne, const int ndof_o,
+                                 const int nquad_o, const Vector &pa,
+                                 const Vector &x_, Vector &y_)
+{
+   constexpr int mnd_o = T_NDOF_O ? T_NDOF_O : DofQuadLimits::HCURL_MAX_D1D - 1;
+   constexpr int mnq_o =
+      T_NQUAD_O ? T_NQUAD_O : DofQuadLimits::HDIV_MAX_D1D - 1;
+   constexpr int mndq = std::max(mnd_o + 1, mnq_o + 1);
+   constexpr int tbatch = curlinterp::NBZ3D(T_NDOF_O, T_NQUAD_O, mndq);
+   MFEM_VERIFY(ndof_o <= mnd_o, "Error: H(curl) order larger than supported");
+   MFEM_VERIFY(nquad_o <= mnq_o, "Error: H(div) order larger than supported");
+   int mnq = std::max(ndof_o + 1, nquad_o + 1);
+   auto pa_data = pa.Read();
+   auto x_d = x_.Read();
+   auto y_d = y_.ReadWrite();
+   mfem::forall_2D_batch<mndq * mndq * (mndq - 1) * tbatch>(
+      ne, mnq * mnq * (mnq - 1), 1, tbatch, [=] MFEM_HOST_DEVICE(int e)
+   {
+      constexpr int MND_O =
+         T_NDOF_O ? T_NDOF_O : DofQuadLimits::HCURL_MAX_D1D - 1;
+      constexpr int MNQ_O =
+         T_NQUAD_O ? T_NQUAD_O : DofQuadLimits::HDIV_MAX_D1D - 1;
+      constexpr int MNDQ = std::max(MND_O + 1, MNQ_O + 1);
+#if defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__)
+      constexpr int nbz = curlinterp::NBZ3D(T_NDOF_O, T_NQUAD_O, MNDQ);
+      int tidz = MFEM_THREAD_ID(z);
+      // Make mnq a local variable since capturing would result in different
+      // captures between host/device versions, and spuriously fails
+      int mnq = std::max(ndof_o + 1, nquad_o + 1);
+#else
+      constexpr int nbz = 1;
+      constexpr int tidz = 0;
+#endif
+      const int NDOF_O = T_NDOF_O ? T_NDOF_O : ndof_o;
+      const int NQUAD_O = T_NQUAD_O ? T_NQUAD_O : nquad_o;
+      const int NDOF_C = NDOF_O + 1;
+      const int NQUAD_C = NQUAD_O + 1;
+      MFEM_SHARED real_t
+      sBG[(MND_O + 1) * MNQ_O + (MND_O + 1) * (MNQ_O + 1) + MND_O * MNQ_O];
+      auto X_ = Reshape(x_d, 3 * NDOF_C * NDOF_C * NDOF_O, ne);
+      auto Y = Reshape(y_d, 3 * NQUAD_C * NQUAD_O * NQUAD_O, ne);
+      auto Gco = Reshape(sBG, NQUAD_O, NDOF_C);
+      auto Bcc = Reshape(sBG + NDOF_C * NQUAD_O, NQUAD_C, NDOF_C);
+      auto Boo =
+         Reshape(sBG + NDOF_C * NQUAD_O + NDOF_C * NQUAD_C, NQUAD_O, NDOF_O);
+      MFEM_SHARED real_t X[3][nbz][MND_O * (MND_O + 1) * (MND_O + 1)];
+      MFEM_SHARED real_t sm0[nbz * 2 * MNDQ * MNDQ * MNDQ];
+      MFEM_SHARED real_t sm1[nbz * 2 * MNDQ * MNDQ * MNDQ];
+
+      // shapes of buffers always use MNDQ to mitigate shared memory bank
+      // conflicts
+      real_t(*DDQ)[nbz][MNDQ][MNDQ][MNDQ] =
+         (real_t(*)[nbz][MNDQ][MNDQ][MNDQ])(sm0);
+      real_t(*DQQ)[nbz][MNDQ][MNDQ][MNDQ] =
+         (real_t(*)[nbz][MNDQ][MNDQ][MNDQ])(sm1);
+      real_t(*QQQ)[nbz][MNDQ][MNDQ][MNDQ] =
+         (real_t(*)[nbz][MNDQ][MNDQ][MNDQ])(sm0);
+      const int offset = NDOF_O * NDOF_C * NDOF_C;
+      const int offsetq = NQUAD_C * NQUAD_O * NQUAD_O;
+      MFEM_FOREACH_THREAD_DIRECT(ix, x, offset)
+      {
+         for (int dim = 0; dim < 3; ++dim)
+         {
+            X[dim][tidz][ix] = X_(ix + dim * offset, e);
+         }
+      }
+      // load basis functions data
+      if (tidz == 0)
+      {
+         auto npts = NDOF_C * NQUAD_O + NDOF_C * NQUAD_C + NDOF_O * NQUAD_O;
+         MFEM_FOREACH_THREAD(ix, x, npts) { sBG[ix] = pa_data[ix]; }
+      }
+      MFEM_SYNC_THREAD;
+
+      // x: Vz Bcc Gco Boo - Vy Bcc Boo Gco
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(qx, dy, dz, x, NQUAD_C, NDOF_C,
+                                           NDOF_O, mnq, mnq, mnq - 1)
+      {
+         real_t u = 0;
+         for (int dx = 0; dx < NDOF_C; ++dx)
+         {
+            u += X[2][tidz][dx + (dy + dz * NDOF_C) * NDOF_C] * Bcc(qx, dx);
+         }
+         DDQ[0][tidz][dz][dy][qx] = u;
+      }
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(qx, dy, dz, x, NQUAD_C, NDOF_O,
+                                           NDOF_C, mnq, mnq - 1, mnq)
+      {
+         real_t u = 0;
+         for (int dx = 0; dx < NDOF_C; ++dx)
+         {
+            u += X[1][tidz][dx + (dy + dz * NDOF_O) * NDOF_C] * Bcc(qx, dx);
+         }
+         DDQ[1][tidz][dz][dy][qx] = u;
+      }
+      MFEM_SYNC_THREAD;
+
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(qx, qy, dz, x, NQUAD_C, NQUAD_O,
+                                           NDOF_O, mnq, mnq, mnq - 1)
+      {
+         real_t u = 0;
+         for (int dy = 0; dy < NDOF_C; ++dy)
+         {
+            u += DDQ[0][tidz][dz][dy][qx] * Gco(qy, dy);
+         }
+         DQQ[0][tidz][dz][qy][qx] = u;
+      }
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(qx, qy, dz, x, NQUAD_C, NQUAD_O,
+                                           NDOF_C, mnq, mnq - 1, mnq)
+      {
+         real_t u = 0;
+         for (int dy = 0; dy < NDOF_O; ++dy)
+         {
+            u += DDQ[1][tidz][dz][dy][qx] * Boo(qy, dy);
+         }
+         DQQ[1][tidz][dz][qy][qx] = u;
+      }
+      MFEM_SYNC_THREAD;
+
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(qx, qy, qz, x, NQUAD_C, NQUAD_O,
+                                           NQUAD_O, mnq, mnq, mnq - 1)
+      {
+         real_t u = 0;
+         for (int dz = 0; dz < NDOF_O; ++dz)
+         {
+            u += DQQ[0][tidz][dz][qy][qx] * Boo(qz, dz);
+         }
+         QQQ[0][tidz][qz][qy][qx] = u;
+      }
+      MFEM_SYNC_THREAD;
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(qx, qy, qz, x, NQUAD_C, NQUAD_O,
+                                           NQUAD_O, mnq, mnq, mnq - 1)
+      {
+         real_t u = 0;
+         for (int dz = 0; dz < NDOF_C; ++dz)
+         {
+            u += DQQ[1][tidz][dz][qy][qx] * Gco(qz, dz);
+         }
+         Y(qx + (qy + qz * NQUAD_O) * NQUAD_C, e) =
+            QQQ[0][tidz][qz][qy][qx] - u;
+      }
+      MFEM_SYNC_THREAD;
+
+      // y: Vx Boo Bcc Gco - Vz Gco Bcc Boo
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(qx, dy, dz, x, NQUAD_O, NDOF_C,
+                                           NDOF_C, mnq - 1, mnq, mnq)
+      {
+         real_t u = 0;
+         for (int dx = 0; dx < NDOF_O; ++dx)
+         {
+            u += X[0][tidz][dx + (dy + dz * NDOF_C) * NDOF_O] * Boo(qx, dx);
+         }
+         DDQ[0][tidz][dz][dy][qx] = u;
+      }
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(qx, dy, dz, x, NQUAD_O, NDOF_C,
+                                           NDOF_O, mnq, mnq, mnq - 1)
+      {
+         real_t u = 0;
+         for (int dx = 0; dx < NDOF_C; ++dx)
+         {
+            u += X[2][tidz][dx + (dy + dz * NDOF_C) * NDOF_C] * Gco(qx, dx);
+         }
+         DDQ[1][tidz][dz][dy][qx] = u;
+      }
+      MFEM_SYNC_THREAD;
+
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(qx, qy, dz, x, NQUAD_O, NQUAD_C,
+                                           NDOF_C, mnq - 1, mnq, mnq)
+      {
+         real_t u = 0;
+         for (int dy = 0; dy < NDOF_C; ++dy)
+         {
+            u += DDQ[0][tidz][dz][dy][qx] * Bcc(qy, dy);
+         }
+         DQQ[0][tidz][dz][qy][qx] = u;
+      }
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(qx, qy, dz, x, NQUAD_O, NQUAD_C,
+                                           NDOF_O, mnq - 1, mnq, mnq)
+      {
+         real_t u = 0;
+         for (int dy = 0; dy < NDOF_C; ++dy)
+         {
+            u += DDQ[1][tidz][dz][dy][qx] * Bcc(qy, dy);
+         }
+         DQQ[1][tidz][dz][qy][qx] = u;
+      }
+      MFEM_SYNC_THREAD;
+
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(qx, qy, qz, x, NQUAD_O, NQUAD_C,
+                                           NQUAD_O, mnq, mnq, mnq - 1)
+      {
+         real_t u = 0;
+         for (int dz = 0; dz < NDOF_C; ++dz)
+         {
+            u += DQQ[0][tidz][dz][qy][qx] * Gco(qz, dz);
+         }
+         QQQ[0][tidz][qz][qy][qx] = u;
+      }
+      MFEM_SYNC_THREAD;
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(qx, qy, qz, x, NQUAD_O, NQUAD_C,
+                                           NQUAD_O, mnq, mnq, mnq - 1)
+      {
+         real_t u = 0;
+         for (int dz = 0; dz < NDOF_O; ++dz)
+         {
+            u += DQQ[1][tidz][dz][qy][qx] * Boo(qz, dz);
+         }
+         Y(qx + (qy + qz * NQUAD_C) * NQUAD_O + offsetq, e) =
+            QQQ[0][tidz][qz][qy][qx] - u;
+      }
+      MFEM_SYNC_THREAD;
+
+      // z: Vy Gco Boo Bcc - Vx Boo Gco Bcc
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(qx, dy, dz, x, NQUAD_O, NDOF_O,
+                                           NDOF_C, mnq, mnq - 1, mnq)
+      {
+         real_t u = 0;
+         for (int dx = 0; dx < NDOF_C; ++dx)
+         {
+            u += X[1][tidz][dx + (dy + dz * NDOF_O) * NDOF_C] * Gco(qx, dx);
+         }
+         DDQ[0][tidz][dz][dy][qx] = u;
+      }
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(qx, dy, dz, x, NQUAD_O, NDOF_C,
+                                           NDOF_C, mnq - 1, mnq, mnq)
+      {
+         real_t u = 0;
+         for (int dx = 0; dx < NDOF_O; ++dx)
+         {
+            u += X[0][tidz][dx + (dy + dz * NDOF_C) * NDOF_O] * Boo(qx, dx);
+         }
+         DDQ[1][tidz][dz][dy][qx] = u;
+      }
+      MFEM_SYNC_THREAD;
+
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(qx, qy, dz, x, NQUAD_O, NQUAD_O,
+                                           NDOF_C, mnq, mnq - 1, mnq)
+      {
+         real_t u = 0;
+         for (int dy = 0; dy < NDOF_O; ++dy)
+         {
+            u += DDQ[0][tidz][dz][dy][qx] * Boo(qy, dy);
+         }
+         DQQ[0][tidz][dz][qy][qx] = u;
+      }
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(qx, qy, dz, x, NQUAD_O, NQUAD_O,
+                                           NDOF_C, mnq, mnq - 1, mnq)
+      {
+         real_t u = 0;
+         for (int dy = 0; dy < NDOF_C; ++dy)
+         {
+            u += DDQ[1][tidz][dz][dy][qx] * Gco(qy, dy);
+         }
+         DQQ[1][tidz][dz][qy][qx] = u;
+      }
+      MFEM_SYNC_THREAD;
+
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(qx, qy, qz, x, NQUAD_O, NQUAD_O,
+                                           NQUAD_C, mnq, mnq - 1, mnq)
+      {
+         real_t u = 0;
+         for (int dz = 0; dz < NDOF_C; ++dz)
+         {
+            u += DQQ[0][tidz][dz][qy][qx] * Bcc(qz, dz);
+         }
+         QQQ[0][tidz][qz][qy][qx] = u;
+      }
+      MFEM_SYNC_THREAD;
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(qx, qy, qz, x, NQUAD_O, NQUAD_O,
+                                           NQUAD_C, mnq, mnq - 1, mnq)
+      {
+         real_t u = 0;
+         for (int dz = 0; dz < NDOF_C; ++dz)
+         {
+            u += DQQ[1][tidz][dz][qy][qx] * Bcc(qz, dz);
+         }
+         Y(qx + (qy + qz * NQUAD_O) * NQUAD_O + 2 * offsetq, e) =
+            QQQ[0][tidz][qz][qy][qx] - u;
+      }
+      MFEM_SYNC_THREAD;
+   });
+}
+
+template <int T_NDOF_O, int T_NQUAD_O>
+void CurlInterpolatorTApply3DSmem(const int ne, const int ndof_o,
+                                  const int nquad_o, const Vector &pa,
+                                  const Vector &x_, Vector &y_)
+{
+   constexpr int mnd_o = T_NDOF_O ? T_NDOF_O : DofQuadLimits::HCURL_MAX_D1D - 1;
+   constexpr int mnq_o =
+      T_NQUAD_O ? T_NQUAD_O : DofQuadLimits::HDIV_MAX_D1D - 1;
+   constexpr int mndq = std::max(mnd_o + 1, mnq_o + 1);
+   constexpr int tbatch = curlinterp::NBZ3D(T_NDOF_O, T_NQUAD_O, mndq);
+   MFEM_VERIFY(ndof_o <= mnd_o, "Error: H(curl) order larger than supported");
+   MFEM_VERIFY(nquad_o <= mnq_o, "Error: H(div) order larger than supported");
+   int mnq = std::max(ndof_o + 1, nquad_o + 1);
+   auto pa_data = pa.Read();
+   auto x_d = x_.Read();
+   auto y_d = y_.ReadWrite();
+   mfem::forall_2D_batch<mndq * mndq * (mndq - 1) * tbatch>(
+      ne, mnq * mnq * (mnq - 1), 1, tbatch, [=] MFEM_HOST_DEVICE(int e)
+   {
+      constexpr int MND_O =
+         T_NDOF_O ? T_NDOF_O : DofQuadLimits::HCURL_MAX_D1D - 1;
+      constexpr int MNQ_O =
+         T_NQUAD_O ? T_NQUAD_O : DofQuadLimits::HDIV_MAX_D1D - 1;
+      constexpr int MNDQ = std::max(MND_O + 1, MNQ_O + 1);
+#if defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__)
+      constexpr int nbz = curlinterp::NBZ3D(T_NDOF_O, T_NQUAD_O, MNDQ);
+      int tidz = MFEM_THREAD_ID(z);
+      // Make mnq a local variable since capturing would result in different
+      // captures between host/device versions, and spuriously fails
+      int mnq = std::max(ndof_o + 1, nquad_o + 1);
+#else
+      constexpr int nbz = 1;
+      constexpr int tidz = 0;
+#endif
+      const int NDOF_O = T_NDOF_O ? T_NDOF_O : ndof_o;
+      const int NQUAD_O = T_NQUAD_O ? T_NQUAD_O : nquad_o;
+      const int NDOF_C = NDOF_O + 1;
+      const int NQUAD_C = NQUAD_O + 1;
+      MFEM_SHARED real_t
+      sBG[(MND_O + 1) * MNQ_O + (MND_O + 1) * (MNQ_O + 1) + MND_O * MNQ_O];
+      auto X_ = Reshape(x_d, 3 * NQUAD_C * NQUAD_O * NQUAD_O, ne);
+      auto Y = Reshape(y_d, 3 * NDOF_C * NDOF_C * NDOF_O, ne);
+      auto Gco = Reshape(sBG, NQUAD_O, NDOF_C);
+      auto Bcc = Reshape(sBG + NDOF_C * NQUAD_O, NQUAD_C, NDOF_C);
+      auto Boo =
+         Reshape(sBG + NDOF_C * NQUAD_O + NDOF_C * NQUAD_C, NQUAD_O, NDOF_O);
+      MFEM_SHARED real_t X[3][nbz][MNQ_O * MNQ_O * (MNQ_O + 1)];
+      MFEM_SHARED real_t sm0[nbz * 2 * MNDQ * MNDQ * MNDQ];
+      MFEM_SHARED real_t sm1[nbz * 2 * MNDQ * MNDQ * MNDQ];
+
+      // shapes of buffers always use MNDQ to mitigate shared memory bank
+      // conflicts
+      real_t(*QQD)[nbz][MNDQ][MNDQ][MNDQ] =
+         (real_t(*)[nbz][MNDQ][MNDQ][MNDQ])(sm0);
+      real_t(*QDD)[nbz][MNDQ][MNDQ][MNDQ] =
+         (real_t(*)[nbz][MNDQ][MNDQ][MNDQ])(sm1);
+      real_t(*DDD)[nbz][MNDQ][MNDQ][MNDQ] =
+         (real_t(*)[nbz][MNDQ][MNDQ][MNDQ])(sm0);
+      const int offset = NDOF_O * NDOF_C * NDOF_C;
+      const int offsetq = NQUAD_C * NQUAD_O * NQUAD_O;
+      MFEM_FOREACH_THREAD_DIRECT(ix, x, offsetq)
+      {
+         for (int dim = 0; dim < 3; ++dim)
+         {
+            X[dim][tidz][ix] = X_(ix + dim * offsetq, e);
+         }
+      }
+      // load basis functions data
+      if (tidz == 0)
+      {
+         auto npts = NDOF_C * NQUAD_O + NDOF_C * NQUAD_C + NDOF_O * NQUAD_O;
+         MFEM_FOREACH_THREAD(ix, x, npts) { sBG[ix] = pa_data[ix]; }
+      }
+      MFEM_SYNC_THREAD;
+
+      // x: Vy Boo Bcc Gco - Vz Boo Gco Bcc
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(dz, qx, qy, x, NDOF_C, NQUAD_O,
+                                           NQUAD_C, mnq, mnq - 1, mnq)
+      {
+         real_t u = 0;
+         for (int qz = 0; qz < NQUAD_O; ++qz)
+         {
+            u += X[1][tidz][qx + (qy + qz * NQUAD_C) * NQUAD_O] * Gco(qz, dz);
+         }
+         QQD[0][tidz][qy][qx][dz] = u;
+      }
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(dz, qx, qy, x, NDOF_C, NQUAD_O,
+                                           NQUAD_O, mnq, mnq, mnq - 1)
+      {
+         real_t u = 0;
+         for (int qz = 0; qz < NQUAD_C; ++qz)
+         {
+            u += X[2][tidz][qx + (qy + qz * NQUAD_O) * NQUAD_O] * Bcc(qz, dz);
+         }
+         QQD[1][tidz][qy][qx][dz] = u;
+      }
+      MFEM_SYNC_THREAD;
+
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(dy, dz, qx, x, NDOF_C, NDOF_C,
+                                           NQUAD_O, mnq, mnq, mnq - 1)
+      {
+         real_t u = 0;
+         for (int qy = 0; qy < NQUAD_C; ++qy)
+         {
+            u += QQD[0][tidz][qy][qx][dz] * Bcc(qy, dy);
+         }
+         QDD[0][tidz][qx][dz][dy] = u;
+      }
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(dy, dz, qx, x, NDOF_C, NDOF_C,
+                                           NQUAD_O, mnq, mnq, mnq - 1)
+      {
+         real_t u = 0;
+         for (int qy = 0; qy < NQUAD_O; ++qy)
+         {
+            u += QQD[1][tidz][qy][qx][dz] * Gco(qy, dy);
+         }
+         QDD[1][tidz][qx][dz][dy] = u;
+      }
+      MFEM_SYNC_THREAD;
+
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(dx, dy, dz, x, NDOF_O, NDOF_C,
+                                           NDOF_C, mnq - 1, mnq, mnq)
+      {
+         real_t u = 0;
+         for (int qx = 0; qx < NQUAD_O; ++qx)
+         {
+            u += QDD[0][tidz][qx][dz][dy] * Boo(qx, dx);
+         }
+         DDD[0][tidz][dz][dy][dx] = u;
+      }
+      MFEM_SYNC_THREAD;
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(dx, dy, dz, x, NDOF_O, NDOF_C,
+                                           NDOF_C, mnq - 1, mnq, mnq)
+      {
+         real_t u = 0;
+         for (int qx = 0; qx < NQUAD_O; ++qx)
+         {
+            u += QDD[1][tidz][qx][dz][dy] * Boo(qx, dx);
+         }
+         Y(dx + (dy + dz * NDOF_C) * NDOF_O, e) =
+            DDD[0][tidz][dz][dy][dx] - u;
+      }
+      MFEM_SYNC_THREAD;
+
+      // y: Vz Gco Boo Bcc - Vx Bcc Boo Gco
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(dz, qx, qy, x, NDOF_C, NQUAD_O,
+                                           NQUAD_O, mnq, mnq, mnq - 1)
+      {
+         real_t u = 0;
+         for (int qz = 0; qz < NQUAD_C; ++qz)
+         {
+            u += X[2][tidz][qx + (qy + qz * NQUAD_O) * NQUAD_O] * Bcc(qz, dz);
+         }
+         QQD[0][tidz][qy][qx][dz] = u;
+      }
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(dz, qx, qy, x, NDOF_C, NQUAD_C,
+                                           NQUAD_O, mnq, mnq, mnq - 1)
+      {
+         real_t u = 0;
+         for (int qz = 0; qz < NQUAD_O; ++qz)
+         {
+            u += X[0][tidz][qx + (qy + qz * NQUAD_O) * NQUAD_C] * Gco(qz, dz);
+         }
+         QQD[1][tidz][qy][qx][dz] = u;
+      }
+      MFEM_SYNC_THREAD;
+
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(dy, dz, qx, x, NDOF_O, NDOF_C,
+                                           NQUAD_O, mnq, mnq, mnq - 1)
+      {
+         real_t u = 0;
+         for (int qy = 0; qy < NQUAD_O; ++qy)
+         {
+            u += QQD[0][tidz][qy][qx][dz] * Boo(qy, dy);
+         }
+         QDD[0][tidz][qx][dz][dy] = u;
+      }
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(dy, dz, qx, x, NDOF_O, NDOF_C,
+                                           NQUAD_C, mnq - 1, mnq, mnq)
+      {
+         real_t u = 0;
+         for (int qy = 0; qy < NQUAD_O; ++qy)
+         {
+            u += QQD[1][tidz][qy][qx][dz] * Boo(qy, dy);
+         }
+         QDD[1][tidz][qx][dz][dy] = u;
+      }
+      MFEM_SYNC_THREAD;
+
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(dx, dy, dz, x, NDOF_C, NDOF_O,
+                                           NDOF_C, mnq, mnq - 1, mnq)
+      {
+         real_t u = 0;
+         for (int qx = 0; qx < NQUAD_O; ++qx)
+         {
+            u += QDD[0][tidz][qx][dz][dy] * Gco(qx, dx);
+         }
+         DDD[0][tidz][dz][dy][dx] = u;
+      }
+      MFEM_SYNC_THREAD;
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(dx, dy, dz, x, NDOF_C, NDOF_O,
+                                           NDOF_C, mnq, mnq - 1, mnq)
+      {
+         real_t u = 0;
+         for (int qx = 0; qx < NQUAD_C; ++qx)
+         {
+            u += QDD[1][tidz][qx][dz][dy] * Bcc(qx, dx);
+         }
+         Y(dx + (dy + dz * NDOF_O) * NDOF_C + offset, e) =
+            DDD[0][tidz][dz][dy][dx] - u;
+      }
+      MFEM_SYNC_THREAD;
+
+      // z: Vx Bcc Gco Boo - Vy Gco Bcc Boo
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(dz, qx, qy, x, NDOF_O, NQUAD_C,
+                                           NQUAD_O, mnq, mnq, mnq - 1)
+      {
+         real_t u = 0;
+         for (int qz = 0; qz < NQUAD_O; ++qz)
+         {
+            u += X[0][tidz][qx + (qy + qz * NQUAD_O) * NQUAD_C] * Boo(qz, dz);
+         }
+         QQD[0][tidz][qy][qx][dz] = u;
+      }
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(dz, qx, qy, x, NDOF_O, NQUAD_O,
+                                           NQUAD_C, mnq, mnq - 1, mnq)
+      {
+         real_t u = 0;
+         for (int qz = 0; qz < NQUAD_O; ++qz)
+         {
+            u += X[1][tidz][qx + (qy + qz * NQUAD_C) * NQUAD_O] * Boo(qz, dz);
+         }
+         QQD[1][tidz][qy][qx][dz] = u;
+      }
+      MFEM_SYNC_THREAD;
+
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(dy, dz, qx, x, NDOF_C, NDOF_O,
+                                           NQUAD_C, mnq, mnq - 1, mnq)
+      {
+         real_t u = 0;
+         for (int qy = 0; qy < NQUAD_O; ++qy)
+         {
+            u += QQD[0][tidz][qy][qx][dz] * Gco(qy, dy);
+         }
+         QDD[0][tidz][qx][dz][dy] = u;
+      }
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(dy, dz, qx, x, NDOF_C, NDOF_O,
+                                           NQUAD_O, mnq, mnq, mnq - 1)
+      {
+         real_t u = 0;
+         for (int qy = 0; qy < NQUAD_C; ++qy)
+         {
+            u += QQD[1][tidz][qy][qx][dz] * Bcc(qy, dy);
+         }
+         QDD[1][tidz][qx][dz][dy] = u;
+      }
+      MFEM_SYNC_THREAD;
+
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(dx, dy, dz, x, NDOF_C, NDOF_C,
+                                           NDOF_O, mnq, mnq, mnq - 1)
+      {
+         real_t u = 0;
+         for (int qx = 0; qx < NQUAD_C; ++qx)
+         {
+            u += QDD[0][tidz][qx][dz][dy] * Bcc(qx, dx);
+         }
+         DDD[0][tidz][dz][dy][dx] = u;
+      }
+      MFEM_SYNC_THREAD;
+      // threads assigned to mitigate bank conflicts
+      MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(dx, dy, dz, x, NDOF_C, NDOF_C,
+                                           NDOF_O, mnq, mnq, mnq - 1)
+      {
+         real_t u = 0;
+         for (int qx = 0; qx < NQUAD_O; ++qx)
+         {
+            u += QDD[1][tidz][qx][dz][dy] * Gco(qx, dx);
+         }
+         Y(dx + (dy + dz * NDOF_C) * NDOF_C + 2 * offset, e) =
+            DDD[0][tidz][dz][dy][dx] - u;
+      }
+      MFEM_SYNC_THREAD;
+   });
+}
+
 } // namespace internal
 
+template <int DIM, int NDOF_O, int NQUAD_O>
+CurlInterpolator::ApplyKernelType
+CurlInterpolator::ApplyPAKernels::Kernel()
+{
+   if constexpr (DIM == 3)
+   {
+      return internal::CurlInterpolatorApply3DSmem<NDOF_O, NQUAD_O>;
+   }
+   MFEM_ABORT("Bad dimension!");
+}
+
+template <int DIM, int NDOF_O, int NQUAD_O>
+CurlInterpolator::ApplyKernelType
+CurlInterpolator::ApplyTPAKernels::Kernel()
+{
+   if constexpr (DIM == 3)
+   {
+      return internal::CurlInterpolatorTApply3DSmem<NDOF_O, NQUAD_O>;
+   }
+   MFEM_ABORT("Bad dimension!");
+}
+
 } // namespace mfem
 
 /// \endcond DO_NOT_DOCUMENT

fem/integ/bilininteg_interp_pa.cpp

@@ -14,9 +14,218 @@
 #include "../gridfunc.hpp"
 #include "../qfunction.hpp"
 
+#include "bilininteg_hcurlhdiv_kernels.hpp"
+
 namespace mfem
 {
 
+namespace
+{
+
+void PAHcurlApplyCurl2D(const int c_dofs1D,
+                        const int o_dofs1D,
+                        const int NE,
+                        const Array<real_t> &Bo_,
+                        const Array<real_t> &Gc_,
+                        const Vector &x_,
+                        Vector &y_)
+{
+   auto Bo = Reshape(Bo_.Read(), o_dofs1D, o_dofs1D);
+   auto Gc = Reshape(Gc_.Read(), o_dofs1D, c_dofs1D);
+   auto X = Reshape(x_.Read(), 2 * c_dofs1D * o_dofs1D, NE);
+   auto Y = Reshape(y_.ReadWrite(), o_dofs1D, o_dofs1D, NE);
+
+   mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
+   {
+      for (int iy = 0; iy < c_dofs1D; ++iy)
+      {
+         for (int ix = 0; ix < o_dofs1D; ++ix)
+         {
+            const real_t xv = X(ix + iy * o_dofs1D, e);
+            for (int oy = 0; oy < o_dofs1D; ++oy)
+            {
+               const real_t gy = Gc(oy, iy);
+               for (int ox = 0; ox < o_dofs1D; ++ox)
+               {
+                  Y(ox, oy, e) -= Bo(ox, ix) * gy * xv;
+               }
+            }
+         }
+      }
+
+      const int y_nd = c_dofs1D * o_dofs1D;
+      for (int iy = 0; iy < o_dofs1D; ++iy)
+      {
+         for (int ix = 0; ix < c_dofs1D; ++ix)
+         {
+            const real_t xv = X(y_nd + ix + iy * c_dofs1D, e);
+            for (int oy = 0; oy < o_dofs1D; ++oy)
+            {
+               const real_t by = Bo(oy, iy);
+               for (int ox = 0; ox < o_dofs1D; ++ox)
+               {
+                  Y(ox, oy, e) += Gc(ox, ix) * by * xv;
+               }
+            }
+         }
+      }
+   });
+}
+
+void PAHcurlApplyCurlTranspose2D(const int c_dofs1D,
+                                 const int o_dofs1D,
+                                 const int NE,
+                                 const Array<real_t> &Bo_,
+                                 const Array<real_t> &Gc_,
+                                 const Vector &x_,
+                                 Vector &y_)
+{
+   auto Bo = Reshape(Bo_.Read(), o_dofs1D, o_dofs1D);
+   auto Gc = Reshape(Gc_.Read(), o_dofs1D, c_dofs1D);
+   auto X = Reshape(x_.Read(), o_dofs1D, o_dofs1D, NE);
+   auto Y = Reshape(y_.ReadWrite(), 2 * c_dofs1D * o_dofs1D, NE);
+
+   mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
+   {
+      for (int dy = 0; dy < c_dofs1D; ++dy)
+      {
+         for (int dx = 0; dx < o_dofs1D; ++dx)
+         {
+            real_t sum = 0.0;
+            for (int oy = 0; oy < o_dofs1D; ++oy)
+            {
+               const real_t gy = Gc(oy, dy);
+               for (int ox = 0; ox < o_dofs1D; ++ox)
+               {
+                  sum -= Bo(ox, dx) * gy * X(ox, oy, e);
+               }
+            }
+            Y(dx + dy * o_dofs1D, e) += sum;
+         }
+      }
+
+      const int y_nd = c_dofs1D * o_dofs1D;
+      for (int dy = 0; dy < o_dofs1D; ++dy)
+      {
+         for (int dx = 0; dx < c_dofs1D; ++dx)
+         {
+            real_t sum = 0.0;
+            for (int oy = 0; oy < o_dofs1D; ++oy)
+            {
+               const real_t by = Bo(oy, dy);
+               for (int ox = 0; ox < o_dofs1D; ++ox)
+               {
+                  sum += Gc(ox, dx) * by * X(ox, oy, e);
+               }
+            }
+            Y(y_nd + dx + dy * c_dofs1D, e) += sum;
+         }
+      }
+   });
+}
+
+void PAHdivApplyCurl2D(const int c_dofs1D,
+                       const int o_dofs1D,
+                       const int NE,
+                       const Array<real_t> &Bc_,
+                       const Array<real_t> &Gc_,
+                       const Vector &x_,
+                       Vector &y_)
+{
+   auto Bc = Reshape(Bc_.Read(), c_dofs1D, c_dofs1D);
+   auto Gc = Reshape(Gc_.Read(), o_dofs1D, c_dofs1D);
+   auto X = Reshape(x_.Read(), c_dofs1D, c_dofs1D, NE);
+   auto Y = Reshape(y_.ReadWrite(), 2 * c_dofs1D * o_dofs1D, NE);
+
+   mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
+   {
+      for (int iy = 0; iy < c_dofs1D; ++iy)
+      {
+         for (int ix = 0; ix < c_dofs1D; ++ix)
+         {
+            const real_t xv = X(ix, iy, e);
+            for (int oy = 0; oy < o_dofs1D; ++oy)
+            {
+               const real_t gy = Gc(oy, iy);
+               for (int ox = 0; ox < c_dofs1D; ++ox)
+               {
+                  Y(ox + oy * c_dofs1D, e) += Bc(ox, ix) * gy * xv;
+               }
+            }
+         }
+      }
+
+      const int y_nd = c_dofs1D * o_dofs1D;
+      for (int iy = 0; iy < c_dofs1D; ++iy)
+      {
+         for (int ix = 0; ix < c_dofs1D; ++ix)
+         {
+            const real_t xv = X(ix, iy, e);
+            for (int oy = 0; oy < c_dofs1D; ++oy)
+            {
+               const real_t by = Bc(oy, iy);
+               for (int ox = 0; ox < o_dofs1D; ++ox)
+               {
+                  Y(y_nd + ox + oy * o_dofs1D, e) -= Gc(ox, ix) * by * xv;
+               }
+            }
+         }
+      }
+   });
+}
+
+void PAHdivApplyCurlTranspose2D(const int c_dofs1D,
+                                const int o_dofs1D,
+                                const int NE,
+                                const Array<real_t> &Bc_,
+                                const Array<real_t> &Gc_,
+                                const Vector &x_,
+                                Vector &y_)
+{
+   auto Bc = Reshape(Bc_.Read(), c_dofs1D, c_dofs1D);
+   auto Gc = Reshape(Gc_.Read(), o_dofs1D, c_dofs1D);
+   auto X = Reshape(x_.Read(), 2 * c_dofs1D * o_dofs1D, NE);
+   auto Y = Reshape(y_.ReadWrite(), c_dofs1D, c_dofs1D, NE);
+
+   mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
+   {
+      for (int dy = 0; dy < o_dofs1D; ++dy)
+      {
+         for (int dx = 0; dx < c_dofs1D; ++dx)
+         {
+            const real_t xv = X(dx + dy * c_dofs1D, e);
+            for (int iy = 0; iy < c_dofs1D; ++iy)
+            {
+               const real_t gy = Gc(dy, iy);
+               for (int ix = 0; ix < c_dofs1D; ++ix)
+               {
+                  Y(ix, iy, e) += Bc(dx, ix) * gy * xv;
+               }
+            }
+         }
+      }
+
+      const int y_nd = c_dofs1D * o_dofs1D;
+      for (int dy = 0; dy < c_dofs1D; ++dy)
+      {
+         for (int dx = 0; dx < o_dofs1D; ++dx)
+         {
+            const real_t xv = X(y_nd + dx + dy * o_dofs1D, e);
+            for (int iy = 0; iy < c_dofs1D; ++iy)
+            {
+               const real_t by = Bc(dy, iy);
+               for (int ix = 0; ix < c_dofs1D; ++ix)
+               {
+                  Y(ix, iy, e) -= Gc(dx, ix) * by * xv;
+               }
+            }
+         }
+      }
+   });
+}
+
+}
+
 // Apply to x corresponding to DOFs in H^1 (domain) the (topological) gradient
 // to get a dof in H(curl) (range). You can think of the range as the "test" space
 // and the domain as the "trial" space, but there's no integration.
@@ -1950,4 +2159,266 @@ void IdentityInterpolator::AddMultTransposePA(const Vector &x, Vector &y) const
    }
 }
 
+void CurlInterpolator::AssemblePA(const FiniteElementSpace &dom_fes,
+                                  const FiniteElementSpace &ran_fes)
+{
+   Mesh *mesh = dom_fes.GetMesh();
+   dim = mesh->Dimension();
+   ne = dom_fes.GetNE();
+   pa_mode_2d = 0;
+   MFEM_VERIFY(ne == ran_fes.GetNE(),
+               "Different meshes for domain and range spaces");
+
+   if (dim == 2)
+   {
+      pa_data.SetSize(0);
+      const FiniteElement *dom_fel = dom_fes.GetTypicalFE();
+      const FiniteElement *ran_fel = ran_fes.GetTypicalFE();
+      const bool hcurl_to_scalar =
+         dynamic_cast<const VectorTensorFiniteElement*>(dom_fel) != NULL &&
+         dom_fel->GetDerivType() == FiniteElement::CURL &&
+         dynamic_cast<const TensorBasisElement*>(ran_fel) != NULL &&
+         ran_fel->GetRangeType() == FiniteElement::SCALAR;
+      const bool scalar_to_hdiv =
+         dynamic_cast<const TensorBasisElement*>(dom_fel) != NULL &&
+         dom_fel->GetRangeType() == FiniteElement::SCALAR &&
+         dynamic_cast<const VectorTensorFiniteElement*>(ran_fel) != NULL &&
+         ran_fel->GetDerivType() == FiniteElement::DIV;
+
+      MFEM_VERIFY(hcurl_to_scalar || scalar_to_hdiv,
+                  "2D CurlInterpolator PA supports H(curl)->scalar and scalar->H(div) only.");
+
+      int closed_basis_type = -1;
+      int open_basis_type = -1;
+      if (hcurl_to_scalar)
+      {
+         const auto *trial_fec = dynamic_cast<const ND_FECollection*>(dom_fes.FEColl());
+         const auto *range_fec = dynamic_cast<const L2_FECollection*>(ran_fes.FEColl());
+         MFEM_VERIFY(trial_fec != NULL, "H(curl) domain must use ND_FECollection.");
+         MFEM_VERIFY(range_fec != NULL, "Scalar range must use L2_FECollection.");
+         MFEM_VERIFY(ran_fel->GetMapType() == FiniteElement::INTEGRAL,
+                     "2D H(curl)->scalar CurlInterpolator PA supports integral-map scalar range spaces only.");
+         closed_basis_type = trial_fec->GetClosedBasisType();
+         open_basis_type = trial_fec->GetOpenBasisType();
+         MFEM_VERIFY(range_fec->GetBasisType() == open_basis_type,
+                     "Domain/range open basis types do not match.");
+         pa_mode_2d = 1;
+      }
+      else
+      {
+         const auto *trial_fec = dynamic_cast<const H1_FECollection*>(dom_fes.FEColl());
+         const auto *range_fec = dynamic_cast<const RT_FECollection*>(ran_fes.FEColl());
+         MFEM_VERIFY(trial_fec != NULL, "Scalar domain must use H1_FECollection.");
+         MFEM_VERIFY(range_fec != NULL, "H(div) range must use RT_FECollection.");
+         closed_basis_type = trial_fec->GetBasisType();
+         open_basis_type = range_fec->GetOpenBasisType();
+         MFEM_VERIFY(range_fec->GetClosedBasisType() == closed_basis_type,
+                     "Domain/range closed basis types do not match.");
+         pa_mode_2d = 2;
+      }
+
+      const int order = hcurl_to_scalar
+                        ? dynamic_cast<const VectorTensorFiniteElement*>(dom_fel)->GetOrder()
+                        : dynamic_cast<const NodalTensorFiniteElement*>(dom_fel)->GetOrder();
+      c_dofs1D = order + 1;
+      o_dofs1D = order;
+
+      closed_dofquad_fe.reset(new H1_SegmentElement(order, closed_basis_type));
+      open_dofquad_fe.reset(new L2_SegmentElement(order - 1, open_basis_type));
+
+      mfem::QuadratureFunctions1D qf1d;
+      mfem::IntegrationRule closed_ir;
+      closed_ir.SetSize(c_dofs1D);
+      qf1d.GaussLobatto(c_dofs1D, &closed_ir);
+
+      mfem::IntegrationRule open_ir;
+      open_ir.SetSize(o_dofs1D);
+      qf1d.GaussLegendre(o_dofs1D, &open_ir);
+
+      maps_C_C = &closed_dofquad_fe->GetDofToQuad(closed_ir, DofToQuad::TENSOR);
+      maps_O_C = &closed_dofquad_fe->GetDofToQuad(open_ir, DofToQuad::TENSOR);
+      maps_O_O = &open_dofquad_fe->GetDofToQuad(open_ir, DofToQuad::TENSOR);
+
+      MFEM_VERIFY(maps_C_C->ndof == c_dofs1D && maps_C_C->nqpt == c_dofs1D, "");
+      MFEM_VERIFY(maps_O_C->ndof == c_dofs1D && maps_O_C->nqpt == o_dofs1D, "");
+      MFEM_VERIFY(maps_O_O->ndof == o_dofs1D && maps_O_O->nqpt == o_dofs1D, "");
+      return;
+   }
+
+   closed_dofquad_fe.reset();
+   open_dofquad_fe.reset();
+   maps_C_C = nullptr;
+   maps_O_C = nullptr;
+   maps_O_O = nullptr;
+
+   const VectorTensorFiniteElement *dom_el =
+      dynamic_cast<const VectorTensorFiniteElement *>(dom_fes.GetTypicalFE());
+   const VectorTensorFiniteElement *ran_el =
+      dynamic_cast<const VectorTensorFiniteElement *>(ran_fes.GetTypicalFE());
+   MFEM_VERIFY(dom_el != NULL, "Only VectorTensorFiniteElement is supported!");
+   MFEM_VERIFY(ran_el != NULL, "Only VectorTensorFiniteElement is supported!");
+   MFEM_VERIFY(dom_el->GetDerivType() == FiniteElement::CURL,
+               "Domain space must be H(curl)");
+   MFEM_VERIFY(ran_el->GetDerivType() == FiniteElement::DIV,
+               "Range space must be H(div)");
+
+   const int dims = dom_el->GetDim();
+   MFEM_VERIFY(dims == 3, "");
+
+   ndof_o = dom_el->GetOrder();
+   int ndof_c = ndof_o + 1;
+   nquad_o = ran_el->GetOrder();
+   int nquad_c = nquad_o + 1;
+
+   // extract the tensor product range dof locations
+   std::vector<real_t> qc(nquad_c);
+   std::vector<real_t> qo(nquad_o);
+   {
+      const IntegrationRule &ran_nodes = ran_el->GetNodes();
+      const Array<int> &quad_map = ran_el->GetDofMap();
+      for (int i = 0; i < nquad_c; ++i)
+      {
+         int idx = UnsignIndex(quad_map[i]);
+         qc[i] = ran_nodes.IntPoint(idx).x;
+      }
+      int offset = ndof_c * ndof_o * ndof_o;
+      for (int i = 0; i < nquad_o; ++i)
+      {
+         int idx = UnsignIndex(quad_map[i + offset]);
+         qo[i] = ran_nodes.IntPoint(idx).x;
+      }
+   }
+
+   // evaluate closed/open 1D basis (and their derivatives) at closed and
+   // open quads
+   // storage order: GCO, BCC, BOO
+   pa_data.SetSize(ndof_c * nquad_o + ndof_c * nquad_c + ndof_o * nquad_o);
+   auto ptr = pa_data.HostWrite();
+   auto &cbasis1d = dom_el->GetBasis1D();
+   auto &obasis1d = dom_el->GetOpenBasis1D();
+   Vector b, g;
+   b.SetSize(ndof_c);
+   g.SetSize(ndof_c);
+   for (int j = 0; j < nquad_o; ++j)
+   {
+      cbasis1d.Eval(qo[j], b, g);
+      for (int i = 0; i < ndof_c; ++i)
+      {
+         ptr[j + i * nquad_o] = g[i];
+      }
+   }
+   ptr += nquad_o * ndof_c;
+
+   for (int j = 0; j < nquad_c; ++j)
+   {
+      cbasis1d.Eval(qc[j], b);
+      for (int i = 0; i < ndof_c; ++i)
+      {
+         ptr[j + i * nquad_c] = b[i];
+      }
+   }
+   ptr += ndof_c * nquad_c;
+
+   b.SetSize(ndof_o);
+   for (int j = 0; j < nquad_o; ++j)
+   {
+      obasis1d.Eval(qo[j], b);
+      for (int i = 0; i < ndof_o; ++i)
+      {
+         ptr[j + i * nquad_o] = b[i];
+      }
+   }
+}
+
+CurlInterpolator::Kernels::Kernels()
+{
+   CurlInterpolator::AddSpecialization<3, 1, 1>();
+   CurlInterpolator::AddSpecialization<3, 2, 2>();
+   CurlInterpolator::AddSpecialization<3, 3, 3>();
+   CurlInterpolator::AddSpecialization<3, 4, 4>();
+   CurlInterpolator::AddSpecialization<3, 5, 5>();
+}
+
+CurlInterpolator::CurlInterpolator() { static Kernels kernels{}; }
+
+void CurlInterpolator::AddMultPA(const Vector &x, Vector &y) const
+{
+   if (dim == 2)
+   {
+      MFEM_VERIFY(maps_C_C != nullptr && maps_O_C != nullptr,
+                  "2D CurlInterpolator PA data is not assembled.");
+      if (pa_mode_2d == 1)
+      {
+         MFEM_VERIFY(maps_O_O != nullptr,
+                     "2D CurlInterpolator scalar curl map is not assembled.");
+         PAHcurlApplyCurl2D(c_dofs1D, o_dofs1D, ne, maps_O_O->B, maps_O_C->G,
+                            x, y);
+      }
+      else if (pa_mode_2d == 2)
+      {
+         PAHdivApplyCurl2D(c_dofs1D, o_dofs1D, ne, maps_C_C->B, maps_O_C->G,
+                           x, y);
+      }
+      else
+      {
+         MFEM_ABORT("Unsupported 2D CurlInterpolator mode.");
+      }
+      return;
+   }
+
+   ApplyPAKernels::Run(dim, ndof_o, nquad_o, ne, ndof_o, nquad_o, pa_data, x, y);
+}
+
+void CurlInterpolator::AddMultTransposePA(const Vector &x, Vector &y) const
+{
+   if (dim == 2)
+   {
+      MFEM_VERIFY(maps_C_C != nullptr && maps_O_C != nullptr,
+                  "2D CurlInterpolator PA data is not assembled.");
+      if (pa_mode_2d == 1)
+      {
+         MFEM_VERIFY(maps_O_O != nullptr,
+                     "2D CurlInterpolator scalar curl map is not assembled.");
+         PAHcurlApplyCurlTranspose2D(c_dofs1D, o_dofs1D, ne, maps_O_O->B,
+                                     maps_O_C->G, x, y);
+      }
+      else if (pa_mode_2d == 2)
+      {
+         PAHdivApplyCurlTranspose2D(c_dofs1D, o_dofs1D, ne, maps_C_C->B,
+                                    maps_O_C->G, x, y);
+      }
+      else
+      {
+         MFEM_ABORT("Unsupported 2D CurlInterpolator mode.");
+      }
+      return;
+   }
+
+   ApplyTPAKernels::Run(dim, ndof_o, nquad_o, ne, ndof_o, nquad_o, pa_data, x, y);
+}
+
+/// \cond DO_NOT_DOCUMENT
+
+CurlInterpolator::ApplyKernelType
+CurlInterpolator::ApplyPAKernels::Fallback(int DIM, int, int)
+{
+   if (DIM == 3)
+   {
+      return internal::CurlInterpolatorApply3DSmem<0, 0>;
+   }
+   MFEM_ABORT("Bad dimension!");
+}
+
+CurlInterpolator::ApplyKernelType
+CurlInterpolator::ApplyTPAKernels::Fallback(int DIM, int, int)
+{
+   if (DIM == 3)
+   {
+      return internal::CurlInterpolatorTApply3DSmem<0, 0>;
+   }
+   MFEM_ABORT("Bad dimension!");
+}
+
+/// \endcond DO_NOT_DOCUMENT
+
 } // namespace mfem

fem/integ/bilininteg_vecmass_pa.cpp

@@ -22,6 +22,8 @@ void VectorMassIntegrator::AssemblePA(const FiniteElementSpace &fes)
 {
    Mesh *mesh = fes.GetMesh();
    const FiniteElement &el = *fes.GetTypicalFE();
+   MFEM_VERIFY(el.GetMapType() == FiniteElement::VALUE,
+               "Only value map type supported");
    ElementTransformation &Trans = *mesh->GetTypicalElementTransformation();
    const auto *ir = IntRule ? IntRule : &MassIntegrator::GetRule(el, el, Trans);
 

fem/lor/lor_ams.cpp

@@ -94,10 +94,10 @@ void BatchedLOR_AMS::Form2DEdgeToVertex_RT(Array<int> &edge2vert)
          const int iv0 = ix + iy*op1;
          const int iv1 = ix1 + iy1*op1;
 
-         // Rotated gradient in 2D (-dy, dx), so flip the sign for the first
-         // component (c == 0).
-         e2v(0, iedge) = (c == 1) ? iv0 : iv1;
-         e2v(1, iedge) = (c == 1) ? iv1 : iv0;
+         // 2D curl (dy, -dx), so flip the sign for the second
+         // component (c == 1).
+         e2v(0, iedge) = (c == 0) ? iv0 : iv1;
+         e2v(1, iedge) = (c == 0) ? iv1 : iv0;
       }
    }
 }

fem/lor/lor_batched.cpp

@@ -142,8 +142,6 @@ static MFEM_HOST_DEVICE int GetAndIncrementNnzIndex(const int i_L, int* I)
 
 int BatchedLORAssembly::FillI(SparseMatrix &A) const
 {
-   static constexpr int Max = 16;
-
    const int nvdof = fes_ho.GetVSize();
 
    const int ndof_per_el = fes_ho.GetTypicalFE()->GetDof();
@@ -165,6 +163,8 @@ int BatchedLORAssembly::FillI(SparseMatrix &A) const
    const auto K = dof_glob2loc_offsets_.Read();
    const auto map = Reshape(sparse_mapping.Read(), nnz_per_row, ndof_per_el);
 
+   Array<int> ij_elts(dof_glob2loc_.Size() * 2);
+   auto d_ij_elts = Reshape(ij_elts.Write(), dof_glob2loc_.Size(), 2);
 
    auto I = A.WriteI();
 
@@ -176,10 +176,10 @@ int BatchedLORAssembly::FillI(SparseMatrix &A) const
       const int sii = el_dof_lex(ii_el, iel_ho);
       const int ii = (sii >= 0) ? sii : -1 -sii;
       // Get number and list of elements containing this DOF
-      int i_elts[Max];
       const int i_offset = K[ii];
       const int i_next_offset = K[ii+1];
       const int i_ne = i_next_offset - i_offset;
+      int *i_elts = &d_ij_elts(i_offset, 0);
       for (int e_i = 0; e_i < i_ne; ++e_i)
       {
          const int si_E = dof_glob2loc[i_offset+e_i]; // signed
@@ -202,7 +202,7 @@ int BatchedLORAssembly::FillI(SparseMatrix &A) const
          }
          else // assembly required
          {
-            int j_elts[Max];
+            int *j_elts = &d_ij_elts(j_offset, 1);
             for (int e_j = 0; e_j < j_ne; ++e_j)
             {
                const int sj_E = dof_glob2loc[j_offset+e_j]; // signed
@@ -269,7 +269,8 @@ void BatchedLORAssembly::FillJAndData(SparseMatrix &A) const
       mfem::forall(nvdof + 1, [=] MFEM_HOST_DEVICE (int i) { I[i] = I2[i]; });
    }
 
-   static constexpr int Max = 16;
+   Array<int> ij_B_el(dof_glob2loc_.Size() * 4);
+   auto d_ij_B_el = Reshape(ij_B_el.Write(), dof_glob2loc_.Size(), 4);
 
    mfem::forall(ndof_per_el*nel_ho, [=] MFEM_HOST_DEVICE (int i)
    {
@@ -279,11 +280,13 @@ void BatchedLORAssembly::FillJAndData(SparseMatrix &A) const
       const int sii = el_dof_lex(ii_el, iel_ho); // signed
       const int ii = (sii >= 0) ? sii : -1 - sii;
       // Get number and list of elements containing this DOF
-      int i_elts[Max];
-      int i_B[Max];
       const int i_offset = K[ii];
       const int i_next_offset = K[ii+1];
       const int i_ne = i_next_offset - i_offset;
+
+      int *i_elts = &d_ij_B_el(i_offset, 0);
+      int *i_B = &d_ij_B_el(i_offset, 1);
+
       for (int e_i = 0; e_i < i_ne; ++e_i)
       {
          const int si_E = dof_glob2loc[i_offset+e_i]; // signed
@@ -312,8 +315,8 @@ void BatchedLORAssembly::FillJAndData(SparseMatrix &A) const
          }
          else // assembly required
          {
-            int j_elts[Max];
-            int j_B[Max];
+            int *j_elts = &d_ij_B_el(j_offset, 2);
+            int *j_B = &d_ij_B_el(j_offset, 3);
             for (int e_j = 0; e_j < j_ne; ++e_j)
             {
                const int sj_E = dof_glob2loc[j_offset+e_j]; // signed

fem/transfer.cpp

@@ -231,9 +231,11 @@ const Operator &InterpolationGridTransfer::BackwardOperator()
 
 L2ProjectionGridTransfer::L2Projection::L2Projection(
    const FiniteElementSpace &fes_ho_, const FiniteElementSpace &fes_lor_,
+   CoefficientWithOrder coeff_ho_, CoefficientWithOrder coeff_lor_,
    MemoryType d_mt_)
    : Operator(fes_lor_.GetVSize(), fes_ho_.GetVSize()),
-     fes_ho(fes_ho_), fes_lor(fes_lor_), d_mt(d_mt_)
+     fes_ho(fes_ho_), fes_lor(fes_lor_), coeff_ho(coeff_ho_),
+     coeff_lor(coeff_lor_), d_mt(d_mt_)
 { }
 
 void L2ProjectionGridTransfer::L2Projection::BuildHo2Lor(
@@ -263,12 +265,13 @@ void L2ProjectionGridTransfer::L2Projection::ElemMixedMass(
    IntegrationPointTransformation& ip_tr,
    DenseMatrix& M_mixed_el) const
 {
-   int order = fe_lor.GetOrder() + fe_ho.GetOrder() + tr_lor->OrderW();
-   const IntegrationRule* ir = &IntRules.Get(geom, order);
+   int order = fe_lor.GetOrder() + fe_ho.GetOrder() + tr_lor->OrderW() +
+               coeff_ho.order;
+   const IntegrationRule &ir = IntRules.Get(geom, order);
    M_mixed_el = 0.0;
-   for (int i = 0; i < ir->GetNPoints(); i++)
+   for (int i = 0; i < ir.GetNPoints(); i++)
    {
-      const IntegrationPoint& ip_lor = ir->IntPoint(i);
+      const IntegrationPoint& ip_lor = ir.IntPoint(i);
       IntegrationPoint ip_ho;
       ip_tr.Transform(ip_lor, ip_ho);
       Vector shape_lor(fe_lor.GetDof());
@@ -284,23 +287,23 @@ void L2ProjectionGridTransfer::L2Projection::ElemMixedMass(
       {
          w *= tr_lor->Weight();
       }
+      if (coeff_ho)
+      {
+         w *= coeff_ho.coeff->Eval(*tr_ho, ip_ho);
+      }
       shape_lor *= w;
       AddMultVWt(shape_lor, shape_ho, M_mixed_el);
    }
 }
 
-void L2ProjectionGridTransfer::L2Projection::ElemMixedMass(
-   Geometry::Type geom, const FiniteElement& fe_ho,
-   const FiniteElement& fe_lor, ElementTransformation* el_tr,
-   IntegrationPointTransformation& ip_tr,
+void L2ProjectionGridTransfer::L2Projection::ElemMixedEvaluation(
+   Geometry::Type geom, const FiniteElement& fe_ho, const FiniteElement& fe_lor,
+   IntegrationPointTransformation& ip_tr, const IntegrationRule& ir,
    DenseMatrix& B_L, DenseMatrix& B_H) const
 {
-   int order = fe_lor.GetOrder() + fe_ho.GetOrder() + el_tr->OrderW();
-   const IntegrationRule* ir = &IntRules.Get(geom, order);
-
-   for (int i = 0; i < ir->GetNPoints(); i++)
+   for (int i = 0; i < ir.GetNPoints(); i++)
    {
-      const IntegrationPoint& ip_lor = ir->IntPoint(i);
+      const IntegrationPoint& ip_lor = ir.IntPoint(i);
       IntegrationPoint ip_ho;
 
       // maps integration point ip_lor -> ip_ho
@@ -320,7 +323,6 @@ void L2ProjectionGridTransfer::L2Projection::ElemMixedMass(
          B_H(i, j) = shape_ho(j);
       }
    }
-
 }
 
 void L2ProjectionGridTransfer::L2Projection::MixedMassEA(
@@ -328,10 +330,11 @@ void L2ProjectionGridTransfer::L2Projection::MixedMassEA(
    const FiniteElementSpace& fes_lor_ea,
    Vector &M_LH, MemoryType d_mt_)
 {
-   Mesh* mesh_ho = fes_ho_ea.GetMesh();
-   Mesh* mesh_lor = fes_lor_ea.GetMesh();
-   int nel_ho = mesh_ho->GetNE();
-   int nel_lor = mesh_lor->GetNE();
+   Mesh &mesh_ho = *fes_ho_ea.GetMesh();
+   Mesh &mesh_lor = *fes_lor_ea.GetMesh();
+
+   const int nel_ho = mesh_ho.GetNE();
+   const int nel_lor = mesh_lor.GetNE();
 
    if (nel_ho == 0)
    {
@@ -339,11 +342,11 @@ void L2ProjectionGridTransfer::L2Projection::MixedMassEA(
       return;
    }
 
-   const CoarseFineTransformations& cf_tr = mesh_lor->GetRefinementTransforms();
+   const CoarseFineTransformations& cf_tr = mesh_lor.GetRefinementTransforms();
 
    int nref_max = 0;
    Array<Geometry::Type> geoms;
-   mesh_ho->GetGeometries(mesh_ho->Dimension(), geoms);
+   mesh_ho.GetGeometries(mesh_ho.Dimension(), geoms);
    for (int ig = 0; ig < geoms.Size(); ++ig)
    {
       Geometry::Type geom = geoms[ig];
@@ -360,130 +363,226 @@ void L2ProjectionGridTransfer::L2Projection::MixedMassEA(
    {
       // Assume all HO elements are LOR in the same way
       const int iho = 0;
-      {
-         Array<int> lor_els;
-         ho2lor.GetRow(iho, lor_els);
-         int nref = ho2lor.RowSize(iho);
-
-         Geometry::Type geom = mesh_ho->GetElementBaseGeometry(iho);
-         const FiniteElement &fe_ho = *fes_ho_ea.GetFE(iho);
-         const FiniteElement &fe_lor = *fes_lor_ea.GetFE(lor_els[0]);
-
-         // Allocate space for DenseTensors
-         ElementTransformation *el_tr = fes_lor_ea.GetElementTransformation(0);
-         int order = fe_lor.GetOrder() + fe_ho.GetOrder() + el_tr->OrderW();
-         const IntegrationRule* ir_ea = &IntRules.Get(geom, order);
-         int qPts = ir_ea->GetNPoints();
-
-         // Containers for the basis functions sampled at quadrature points
-         B_L.SetSize(qPts, fe_lor.GetDof(), nref, d_mt);
-         B_H.SetSize(qPts, fe_ho.GetDof(), nref, d_mt);
-         D.SetSize(qPts, nref, nel_ho, d_mt);
-
-         const GeometricFactors *geo_facts =
-            mesh_lor->GetGeometricFactors(*ir_ea, GeometricFactors::DETERMINANTS);
-
-         MFEM_ASSERT(nel_ho*nref == nel_lor, "we expect nel_ho*nref == nel_lor");
-
-         // Setup data at quadrature points
-         // TODO add support for user coefficient
-         const auto W = Reshape(ir_ea->GetWeights().Read(), qPts);
-         const auto J = Reshape(geo_facts->detJ.Read(), qPts, nel_lor);
-         const auto d_D = Reshape(D.Write(), qPts, nref, nel_ho);
-
-         mfem::forall(qPts * nref * nel_ho, [=] MFEM_HOST_DEVICE (int tid)
-         {
-            const int q    = tid % qPts;
-            const int iref = (tid / qPts) % nref;
-            const int iho  = (tid / (qPts * nref)) % nel_ho;
-
-            const int lo_el_id = iref + nref*iho;
-            const real_t detJ = J(q, lo_el_id);
-
-            d_D(q, iref, iho) = W(q) * detJ;
-
-         });
-
-         emb_tr.SetIdentityTransformation(geom);
-         const DenseTensor &pmats = cf_tr.point_matrices[geom];
-
-         // Collect the basis functions
-         for (int iref = 0; iref < nref; ++iref)
-         {
-            int ilor = lor_els[iref];
-            // Now assemble the block-row of the mixed mass matrix associated
-            // with integrating HO functions against LOR functions on the LOR
-            // sub-element.
-
-            // Create the transformation that embeds the fine low-order element
-            // within the coarse high-order element in reference space
-            emb_tr.SetPointMat(pmats(cf_tr.embeddings[ilor].matrix));
-
-            DenseMatrix &b_lo = B_L(ilor);
-            DenseMatrix &b_ho = B_H(ilor);
-
-            ElemMixedMass(geom, fe_ho, fe_lor, el_tr, ip_tr, b_lo, b_ho);
-
-         } // loop over subcells of ho element
-         // end of quadrature point setup
-      }
-
-   } // completed setup of basis function and quadrature point
-
-   // Assemble mixed mass matrix
-   {
-      int iho = 0;
       Array<int> lor_els;
       ho2lor.GetRow(iho, lor_els);
-      int nref = ho2lor.RowSize(iho);
+      const int nref = ho2lor.RowSize(iho);
+      MFEM_VERIFY(nel_ho*nref == nel_lor, "we expect nel_ho*nref == nel_lor");
+
+      Geometry::Type geom = mesh_ho.GetElementBaseGeometry(iho);
+
+      emb_tr.SetIdentityTransformation(geom);
+      const DenseTensor &pmats = cf_tr.point_matrices[geom];
 
       const FiniteElement &fe_ho = *fes_ho_ea.GetFE(iho);
       const FiniteElement &fe_lor = *fes_lor_ea.GetFE(lor_els[0]);
-      const int ndof_ho = fe_ho.GetDof();
-      const int ndof_lor = fe_lor.GetDof();
 
-      const int qPts = D.SizeI();
+      // Allocate space for DenseTensors
+      ElementTransformation &el_tr = *mesh_lor.GetTypicalElementTransformation();
+      const int order = fe_lor.GetOrder() + fe_ho.GetOrder() + el_tr.OrderW()
+                        + coeff_ho.order;
+      const IntegrationRule &ir_ea = IntRules.Get(geom, order);
+      const int qPts = ir_ea.GetNPoints();
 
-      M_LH.SetSize(ndof_lor*ndof_ho*nref*nel_ho, d_mt);
+      // Containers for the basis functions sampled at quadrature points
+      B_L.SetSize(qPts, fe_lor.GetDof(), nref, d_mt);
+      B_H.SetSize(qPts, fe_ho.GetDof(), nref, d_mt);
+      D.SetSize(qPts, nref, nel_ho, d_mt);
 
-      // Rows x columns
-      // Recall MFEM is column major
-      // rows x columns is inverted  - matrix is ndof_lor x ndof_ho
-      auto v_M_LH = Reshape(M_LH.Write(), ndof_lor, ndof_ho, nref,
-                            nel_ho);
+      const GeometricFactors *geo_facts =
+         mesh_lor.GetGeometricFactors(ir_ea, GeometricFactors::DETERMINANTS);
 
-      const int fe_ho_ndof = fe_ho.GetDof();
-      const int fe_lor_ndof = fe_lor.GetDof();
+      Vector coeff_vec(qPts*nel_lor);
+      coeff_vec.UseDevice(true);
 
-      auto d_B_L = Reshape(B_L.Read(), qPts, fe_lor_ndof, nref);
-      auto d_B_H = Reshape(B_H.Read(), qPts, fe_ho_ndof, nref);
-      auto d_D   = Reshape(D.Read(), qPts, nref, nel_ho);
+      const int dim = mesh_ho.Dimension();
+      const int nq1d = (int)floor(pow(ir_ea.Size(), 1.0/dim) + 0.5);
+      const int nref_1d = (int)floor(pow(nref, 1.0/dim) + 0.5);
 
-      mfem::forall(fe_ho_ndof*nref*nel_ho, [=] MFEM_HOST_DEVICE (int idx)
+      if (!coeff_ho)
       {
-         const int bh = idx % fe_ho_ndof;
-         const int iref = (idx / fe_ho_ndof) % nref;
-         const int iho = idx / fe_ho_ndof / nref;
-         // (B_lo_dofs x Q) x (Q x B_ho_dofs)
-         for (int bl = 0; bl < fe_lor_ndof; ++bl)
+         coeff_vec = 1.0;
+      }
+      else if (UsesTensorBasis(fes_ho) &&
+               nq1d*nref_1d <= DeviceDofQuadLimits::Get().MAX_Q1D)
+      {
+         // Fast coefficient evaluation for tensor-product case. We create a
+         // "composite" quadrature rule in the high-order element that is the
+         // union of the quadrature rules within each of the low-order-refined
+         // subelements.
+         //
+         // NOTE: if the integration rule order is high and there are many LOR
+         // subelements, this can create a very big quadrature rule. That is
+         // why we need to check that we do not exceed MAX_Q1D. If we do, then
+         // we fall back on the slower "legacy" evaluation.
+
+         // Construct the composite rule as a tensor-product of the 1D LOR rule.
+         IntegrationRule ir_ho = [&]()
          {
-            real_t dot = 0.0;
-            for (int qi=0; qi<qPts; ++qi)
+            IntegrationRule ir_ho_1d(nq1d * nref_1d);
+            for (int iref = 0; iref < nref_1d; ++iref)
             {
-               dot += d_B_L(qi, bl, iref) *  d_D(qi, iref, iho) * d_B_H(qi, bh, iref);
+               const real_t a = pmats(cf_tr.embeddings[iref].matrix)(0,0);
+               const real_t b = pmats(cf_tr.embeddings[iref].matrix)(0,1);
+               for (int iq = 0; iq < nq1d; ++iq)
+               {
+                  ir_ho_1d[iq + iref*nq1d].x = a + ir_ea[iq].x*(b - a);
+               }
+            }
+            if (dim == 1) { return ir_ho_1d; }
+            else if (dim == 2) { return IntegrationRule(ir_ho_1d, ir_ho_1d); }
+            else { return IntegrationRule(ir_ho_1d, ir_ho_1d, ir_ho_1d); }
+         }();
+
+         // Project the high-order coefficient on the high-order composite rule.
+         QuadratureSpace qs(mesh_ho, ir_ho);
+         CoefficientVector coeff_vec_ho(*coeff_ho.coeff, qs);
+
+         // Permute the coefficient values to the expected LOR ordering.
+         const int nq_ho = ir_ho.Size();
+         const auto d_Q_ho = Reshape(coeff_vec_ho.Read(), nq_ho, nel_ho);
+         const auto d_Q = Reshape(coeff_vec.Write(), qPts, nel_lor);
+
+         mfem::forall(nq_ho * nel_ho, [=] MFEM_HOST_DEVICE (int ii)
+         {
+            const int e_ho = ii / nq_ho;
+            const int iq_ho = ii % nq_ho;
+
+            int iq_tensor = iq_ho;
+            int iq_lor = 0;
+            int iref = 0;
+            int iq_stride = 1;
+            int iref_stride = 1;
+            const int nq_ho_1d = nq1d*nref_1d;
+            for (int d = 0; d < dim; ++d)
+            {
+               const int iq_ho_1d = iq_tensor % nq_ho_1d;
+               iq_tensor /= nq_ho_1d;
+
+               iq_lor += (iq_ho_1d % nq1d)*iq_stride;
+               iref += (iq_ho_1d / nq1d)*iref_stride;
+               iq_stride *= nq1d;
+               iref_stride *= nref_1d;
+            }
+            const int e_lor = iref + e_ho*nref;
+
+            d_Q(iq_lor, e_lor) = d_Q_ho(iq_ho, e_ho);
+         });
+      }
+      else
+      {
+         // Legacy/fallback coefficient evaluation for non-tensor-product cases
+         // or when the number of quadrature points is too large for the device
+         // kernels.
+         IntegrationPoint ip_ho;
+         for (int e_ho = 0; e_ho < nel_ho; ++e_ho)
+         {
+            ElementTransformation &ho_tr = *mesh_ho.GetElementTransformation(e_ho);
+            for (int iref = 0; iref < nref; ++iref)
+            {
+               const int e_lor = iref + e_ho*nref;
+               emb_tr.SetPointMat(pmats(cf_tr.embeddings[e_lor].matrix));
+
+               for (int iq = 0; iq < qPts; ++iq)
+               {
+                  const IntegrationPoint &ip_lor = ir_ea[iq];
+                  ip_tr.Transform(ip_lor, ip_ho);
+                  ho_tr.SetIntPoint(&ip_ho);
+                  coeff_vec[iq + e_lor*qPts] = coeff_ho.coeff->Eval(ho_tr, ip_ho);
+               }
             }
-            // column major storage
-            v_M_LH(bl, bh, iref, iho) = dot;
          }
+      }
+
+      // Setup data at quadrature points
+      const auto W = Reshape(ir_ea.GetWeights().Read(), qPts);
+      const auto J = Reshape(geo_facts->detJ.Read(), qPts, nel_lor);
+      const auto d_D = Reshape(D.Write(), qPts, nref, nel_ho);
+      const auto d_Q = Reshape(coeff_vec.Read(), qPts, nel_lor);
+
+      mfem::forall(qPts * nref * nel_ho, [=] MFEM_HOST_DEVICE (int tid)
+      {
+         const int q    = tid % qPts;
+         const int iref = (tid / qPts) % nref;
+         const int iho  = (tid / (qPts * nref)) % nel_ho;
+
+         const int lo_el_id = iref + nref*iho;
+         const real_t detJ = J(q, lo_el_id);
+
+         d_D(q, iref, iho) = W(q) * d_Q(q, lo_el_id) * detJ;
       });
-   } // end of mixed assembly mass matrix
+
+      // Collect the basis functions
+      for (int iref = 0; iref < nref; ++iref)
+      {
+         int ilor = lor_els[iref];
+         // Now assemble the block-row of the mixed mass matrix associated
+         // with integrating HO functions against LOR functions on the LOR
+         // sub-element.
+
+         // Create the transformation that embeds the fine low-order element
+         // within the coarse high-order element in reference space
+         emb_tr.SetPointMat(pmats(cf_tr.embeddings[ilor].matrix));
+
+         DenseMatrix &b_lo = B_L(ilor);
+         DenseMatrix &b_ho = B_H(ilor);
+
+         ElemMixedEvaluation(geom, fe_ho, fe_lor, ip_tr, ir_ea, b_lo, b_ho);
+      } // loop over subcells of ho element
+      // end of quadrature point setup
+   } // completed setup of basis function and quadrature point
+
+   // Assemble mixed mass matrix
+   int iho = 0;
+   Array<int> lor_els;
+   ho2lor.GetRow(iho, lor_els);
+   int nref = ho2lor.RowSize(iho);
+
+   const FiniteElement &fe_ho = *fes_ho_ea.GetFE(iho);
+   const FiniteElement &fe_lor = *fes_lor_ea.GetFE(lor_els[0]);
+   const int ndof_ho = fe_ho.GetDof();
+   const int ndof_lor = fe_lor.GetDof();
+
+   const int qPts = D.SizeI();
+
+   M_LH.SetSize(ndof_lor*ndof_ho*nref*nel_ho, d_mt);
+
+   // Rows x columns
+   // Recall MFEM is column major
+   // rows x columns is inverted  - matrix is ndof_lor x ndof_ho
+   auto v_M_LH = Reshape(M_LH.Write(), ndof_lor, ndof_ho, nref,
+                         nel_ho);
+
+   const int fe_ho_ndof = fe_ho.GetDof();
+   const int fe_lor_ndof = fe_lor.GetDof();
+
+   auto d_B_L = Reshape(B_L.Read(), qPts, fe_lor_ndof, nref);
+   auto d_B_H = Reshape(B_H.Read(), qPts, fe_ho_ndof, nref);
+   auto d_D   = Reshape(D.Read(), qPts, nref, nel_ho);
+
+   mfem::forall(fe_ho_ndof*nref*nel_ho, [=] MFEM_HOST_DEVICE (int idx)
+   {
+      const int bh = idx % fe_ho_ndof;
+      const int iref = (idx / fe_ho_ndof) % nref;
+      const int iho = idx / fe_ho_ndof / nref;
+      // (B_lo_dofs x Q) x (Q x B_ho_dofs)
+      for (int bl = 0; bl < fe_lor_ndof; ++bl)
+      {
+         real_t dot = 0.0;
+         for (int qi=0; qi<qPts; ++qi)
+         {
+            dot += d_B_L(qi, bl, iref) *  d_D(qi, iref, iho) * d_B_H(qi, bh, iref);
+         }
+         // column major storage
+         v_M_LH(bl, bh, iref, iho) = dot;
+      }
+   });
 }
 
 L2ProjectionGridTransfer::L2ProjectionL2Space::L2ProjectionL2Space
 (const FiniteElementSpace &fes_ho_, const FiniteElementSpace &fes_lor_,
+ CoefficientWithOrder coeff_ho_, CoefficientWithOrder coeff_lor_,
  const bool use_ea_, MemoryType d_mt_)
-   : L2Projection(fes_ho_, fes_lor_, d_mt_),
-     use_ea(use_ea_)
+   : L2Projection(fes_ho_, fes_lor_, coeff_ho_, coeff_lor_, d_mt_), use_ea(use_ea_)
 {
    if (use_ea)
    {
@@ -559,7 +658,11 @@ L2ProjectionGridTransfer::L2ProjectionL2Space::L2ProjectionL2Space
       DenseMatrix Minv_lor(ndof_lor*nref, ndof_lor*nref);
       DenseMatrix M_mixed(ndof_lor*nref, ndof_ho);
 
-      MassIntegrator mi;
+      MassIntegrator mi = [&]()
+      {
+         return coeff_lor ? MassIntegrator(*coeff_lor.coeff) : MassIntegrator();
+      }();
+
       DenseMatrix M_lor_el(ndof_lor, ndof_lor);
       DenseMatrixInverse Minv_lor_el(&M_lor_el);
       DenseMatrix M_lor(ndof_lor*nref, ndof_lor*nref);
@@ -577,6 +680,10 @@ L2ProjectionGridTransfer::L2ProjectionL2Space::L2ProjectionL2Space
          // Assemble the low-order refined mass matrix and invert locally
          int ilor = lor_els[iref];
          ElementTransformation *tr_lor = fes_lor.GetElementTransformation(ilor);
+
+         const int order = 2*fe_lor.GetOrder() + tr_lor->OrderW() + coeff_lor.order;
+         mi.SetIntegrationRule(IntRules.Get(geom, order));
+
          mi.AssembleElementMatrix(fe_lor, *tr_lor, M_lor_el);
          M_lor.CopyMN(M_lor_el, iref*ndof_lor, iref*ndof_lor);
          Minv_lor_el.Factor();
@@ -668,25 +775,22 @@ void L2ProjectionGridTransfer::L2ProjectionL2Space::EAL2ProjectionL2Space()
    // Need to compute M_L
    // Note: Using user-inputted M_LH IntegrationRule ir
    // (higher order than needed) in order to re-use coeff
-   MassIntegrator mi;
+   MassIntegrator mi = [&]()
+   {
+      return coeff_lor ? MassIntegrator(*coeff_lor.coeff) : MassIntegrator();
+   }();
+
+   const int order = 2*fes_lor.GetMaxElementOrder()
+                     + mesh_lor->GetTypicalElementTransformation()->OrderW()
+                     + coeff_lor.order;
+   mi.SetIntegrationRule(
+      IntRules.Get(mesh_lor->GetTypicalElementGeometry(), order));
 
    Vector M_ea_lor;
-   int ndof_lor;
-   int ndof_ho;
-   int nref;
-   {
-      int iho = 0;
-      Array<int> lor_els;
-      ho2lor.GetRow(iho, lor_els);
-      nref = ho2lor.RowSize(iho);
-
-      const FiniteElement &fe_ho = *fes_ho.GetFE(iho);
-      const FiniteElement &fe_lor = *fes_lor.GetFE(lor_els[0]);
-      ndof_ho = fe_ho.GetDof();
-      ndof_lor = fe_lor.GetDof();
-
-      M_ea_lor.SetSize(ndof_lor*ndof_lor*nel_lor, d_mt);
-   }
+   const int ndof_lor = fes_lor.GetTypicalFE()->GetDof();
+   const int ndof_ho = fes_ho.GetTypicalFE()->GetDof();
+   const int nref = ho2lor.RowSize(0);
+   M_ea_lor.SetSize(ndof_lor*ndof_lor*nel_lor, d_mt);
 
    const bool add = false;
    mi.AssembleEA(fes_lor, M_ea_lor, add);
@@ -1032,8 +1136,9 @@ void L2ProjectionGridTransfer::L2ProjectionL2Space::EAProlongateTranspose(
 
 L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space(
    const FiniteElementSpace& fes_ho_, const FiniteElementSpace& fes_lor_,
+   CoefficientWithOrder coeff_ho_, CoefficientWithOrder coeff_lor_,
    const bool use_ea_, MemoryType d_mt_)
-   : L2Projection(fes_ho_, fes_lor_, d_mt_),
+   : L2Projection(fes_ho_, fes_lor_, coeff_ho_, coeff_lor_, d_mt_),
      use_ea(use_ea_)
 {
 
@@ -1092,8 +1197,9 @@ L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space(
 
 L2ProjectionGridTransfer::L2ProjectionH1Space::L2ProjectionH1Space(
    const ParFiniteElementSpace& pfes_ho, const ParFiniteElementSpace& pfes_lor,
+   CoefficientWithOrder coeff_ho_, CoefficientWithOrder coeff_lor_,
    const bool use_ea_, MemoryType d_mt_)
-   : L2Projection(pfes_ho, pfes_lor, d_mt_),
+   : L2Projection(pfes_ho, pfes_lor, coeff_ho_, coeff_lor_, d_mt_),
      use_ea(use_ea_), pcg(pfes_ho.GetComm())
 {
 
@@ -1165,12 +1271,12 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::SetupPCG()
 
 void L2ProjectionGridTransfer::L2ProjectionH1Space::EAL2ProjectionH1Space()
 {
-   Mesh* mesh_ho = fes_ho.GetMesh();
-   Mesh* mesh_lor = fes_lor.GetMesh();
-   int nel_ho = mesh_ho->GetNE();
-   int nel_lor = mesh_lor->GetNE();
-   int ndof_ho = fes_ho.GetNDofs();
-   int ndof_lor = fes_lor.GetNDofs();
+   Mesh &mesh_ho = *fes_ho.GetMesh();
+   Mesh &mesh_lor = *fes_lor.GetMesh();
+   const int nel_ho = mesh_ho.GetNE();
+   const int nel_lor = mesh_lor.GetNE();
+   const int ndof_ho = fes_ho.GetNDofs();
+   const int ndof_lor = fes_lor.GetNDofs();
 
    // If the local mesh is empty, skip all computations
    if (nel_ho == 0)
@@ -1178,11 +1284,11 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::EAL2ProjectionH1Space()
       return;
    }
 
-   const CoarseFineTransformations& cf_tr = mesh_lor->GetRefinementTransforms();
+   const CoarseFineTransformations& cf_tr = mesh_lor.GetRefinementTransforms();
 
    int nref_max = 0;
    Array<Geometry::Type> geoms;
-   mesh_ho->GetGeometries(mesh_ho->Dimension(), geoms);
+   mesh_ho.GetGeometries(mesh_ho.Dimension(), geoms);
    for (int ig = 0; ig < geoms.Size(); ++ig)
    {
       Geometry::Type geom = geoms[ig];
@@ -1205,7 +1311,8 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::EAL2ProjectionH1Space()
 
    BilinearForm Mho(fes_ho_scalar.get());
    Mho.SetAssemblyLevel(AssemblyLevel::PARTIAL);
-   Mho.AddDomainIntegrator(new MassIntegrator);
+   Mho.AddDomainIntegrator(coeff_ho ? new MassIntegrator(*coeff_ho.coeff)
+                           : new MassIntegrator);
    Mho.Assemble();
 
    // Processor local lumped Mass
@@ -1215,7 +1322,16 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::EAL2ProjectionH1Space()
 
    BilinearForm Mlor(fes_lor_scalar.get());
    Mlor.SetAssemblyLevel(AssemblyLevel::PARTIAL);
-   Mlor.AddDomainIntegrator(new MassIntegrator);
+   {
+      MassIntegrator *mi = coeff_lor ? new MassIntegrator(*coeff_lor.coeff)
+                           : new MassIntegrator;
+      const int order = 2*fes_lor.GetMaxElementOrder()
+                        + mesh_lor.GetTypicalElementTransformation()->OrderW()
+                        + coeff_lor.order;
+      mi->SetIntegrationRule(
+         IntRules.Get(mesh_lor.GetTypicalElementGeometry(), order));
+      Mlor.AddDomainIntegrator(mi);
+   }
    Mlor.Assemble();
 
    Vector ones_lor(Mlor.Width()); ones_lor = 1.0;
@@ -1228,15 +1344,14 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::EAL2ProjectionH1Space()
    MixedMassEA(fes_ho, fes_lor, M_LH_ea, d_mt);
 
    // Set ownership
-   M_LH_local_op = new H1SpaceMixedMassOperator(fes_ho_scalar.get(),
-                                                fes_lor_scalar.get(),
-                                                &ho2lor,
-                                                &M_LH_ea);
+   M_LH.reset(new H1SpaceMixedMassOperator(fes_ho_scalar.get(),
+                                           fes_lor_scalar.get(),
+                                           &ho2lor,
+                                           &M_LH_ea));
 
    ML_inv_vea.reset(new H1SpaceLumpedMassOperator(fes_ho_scalar.get(),
                                                   fes_lor_scalar.get(),
                                                   ML_inv_ea));
-   M_LH.reset(M_LH_local_op);
    R.reset(new ProductOperator(ML_inv_vea.get(), M_LH.get(), false,
                                false));
 
@@ -1253,18 +1368,18 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::EAL2ProjectionH1Space()
 void L2ProjectionGridTransfer::L2ProjectionH1Space::EAL2ProjectionH1Space
 (const ParFiniteElementSpace& pfes_ho, const ParFiniteElementSpace& pfes_lor)
 {
-   Mesh* mesh_ho = pfes_ho.GetParMesh();
-   Mesh* mesh_lor = pfes_lor.GetParMesh();
-   int nel_ho = mesh_ho->GetNE();
-   int nel_lor = mesh_lor->GetNE();
+   Mesh &mesh_ho = *pfes_ho.GetParMesh();
+   Mesh &mesh_lor = *pfes_lor.GetParMesh();
+   int nel_ho = mesh_ho.GetNE();
+   int nel_lor = mesh_lor.GetNE();
    int ndof_ho = pfes_ho.GetNDofs();
    int ndof_lor = pfes_lor.GetNDofs();
 
-   const CoarseFineTransformations& cf_tr = mesh_lor->GetRefinementTransforms();
+   const CoarseFineTransformations& cf_tr = mesh_lor.GetRefinementTransforms();
 
    int nref_max = 0;
    Array<Geometry::Type> geoms;
-   mesh_ho->GetGeometries(mesh_ho->Dimension(), geoms);
+   mesh_ho.GetGeometries(mesh_ho.Dimension(), geoms);
    for (int ig = 0; ig < geoms.Size(); ++ig)
    {
       Geometry::Type geom = geoms[ig];
@@ -1287,7 +1402,8 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::EAL2ProjectionH1Space
 
    ParBilinearForm pMho(pfes_ho_scalar.get());
    pMho.SetAssemblyLevel(AssemblyLevel::PARTIAL);
-   pMho.AddDomainIntegrator(new MassIntegrator);
+   pMho.AddDomainIntegrator(coeff_ho ? new MassIntegrator(*coeff_ho.coeff)
+                            : new MassIntegrator);
    pMho.Assemble();
 
    // Processor local lumped Mass
@@ -1297,7 +1413,16 @@ void L2ProjectionGridTransfer::L2ProjectionH1Space::EAL2ProjectionH1Space
 
    ParBilinearForm pMlor(pfes_lor_scalar.get());
    pMlor.SetAssemblyLevel(AssemblyLevel::PARTIAL);
-   pMlor.AddDomainIntegrator(new MassIntegrator);
+   {
+      MassIntegrator *mi = coeff_lor ? new MassIntegrator(*coeff_lor.coeff)
+                           : new MassIntegrator;
+      const int order = 2*fes_lor.GetMaxElementOrder()
+                        + mesh_lor.GetTypicalElementTransformation()->OrderW()
+                        + coeff_lor.order;
+      mi->SetIntegrationRule(
+         IntRules.Get(mesh_lor.GetTypicalElementGeometry(), order));
+      pMlor.AddDomainIntegrator(mi);
+   }
    pMlor.Assemble();
 
    Vector ones_lor(pMlor.Width()); ones_lor = 1.0;
@@ -1570,7 +1695,7 @@ std::unique_ptr<SparseMatrix>>
          int ilor = lor_els[iref];
          ElementTransformation* el_tr = fes_lor.GetElementTransformation(ilor);
 
-         int order = 2 * fe_lor.GetOrder() + el_tr->OrderW();
+         int order = 2 * fe_lor.GetOrder() + el_tr->OrderW() + coeff_lor.order;
          const IntegrationRule* ir = &IntRules.Get(geom, order);
          ML_el = 0.0;
          for (int i = 0; i < ir->GetNPoints(); ++i)
@@ -1578,7 +1703,13 @@ std::unique_ptr<SparseMatrix>>
             const IntegrationPoint& ip_lor = ir->IntPoint(i);
             fe_lor.CalcShape(ip_lor, shape_lor);
             el_tr->SetIntPoint(&ip_lor);
-            ML_el += (shape_lor *= (el_tr->Weight() * ip_lor.weight));
+            real_t w = ip_lor.weight;
+            if (coeff_lor)
+            {
+               w *= coeff_lor.coeff->Eval(*el_tr, ip_lor);
+            }
+            shape_lor *= el_tr->Weight() * w;
+            ML_el += shape_lor;
          }
          fes_lor.GetElementDofs(ilor, dofs_lor);
          ML_inv.AddElementVector(dofs_lor, ML_el);
@@ -2024,8 +2155,8 @@ void L2ProjectionGridTransfer::BuildF()
    {
       if (!Parallel())
       {
-         F = new L2ProjectionH1Space(dom_fes, ran_fes,
-                                     use_ea, d_mt);
+         F = new L2ProjectionH1Space(
+            dom_fes, ran_fes, coeff_ho, coeff_lor, use_ea, d_mt);
       }
       else
       {
@@ -2034,15 +2165,15 @@ void L2ProjectionGridTransfer::BuildF()
             static_cast<mfem::ParFiniteElementSpace&>(dom_fes);
          const mfem::ParFiniteElementSpace& ran_pfes =
             static_cast<mfem::ParFiniteElementSpace&>(ran_fes);
-         F = new L2ProjectionH1Space(dom_pfes, ran_pfes,
-                                     use_ea, d_mt);
+         F = new L2ProjectionH1Space(
+            dom_pfes, ran_pfes, coeff_ho, coeff_lor, use_ea, d_mt);
 #endif
       }
    }
    else
    {
-      F = new L2ProjectionL2Space(dom_fes, ran_fes,
-                                  use_ea, d_mt);
+      F = new L2ProjectionL2Space(
+         dom_fes, ran_fes, coeff_ho, coeff_lor, use_ea, d_mt);
    }
 }
 

fem/transfer.hpp

@@ -19,6 +19,8 @@
 #include "pfespace.hpp"
 #endif
 
+#include <cstddef>
+
 namespace mfem
 {
 
@@ -162,6 +164,18 @@ public:
 };
 
 
+struct CoefficientWithOrder
+{
+   Coefficient *coeff;
+   int order;
+   CoefficientWithOrder() : coeff(nullptr), order(0) { }
+   CoefficientWithOrder(std::nullptr_t) : coeff(nullptr), order(0) { }
+   CoefficientWithOrder(Coefficient &coeff_) : coeff(&coeff_), order(1) { }
+   CoefficientWithOrder(Coefficient &coeff_, int order_)
+      : coeff(&coeff_), order(order_) { }
+   operator bool() const { return coeff != nullptr; }
+};
+
 /** @brief Transfer data in L2 and H1 finite element spaces between a coarse
     mesh and an embedded refined mesh using L2 projection. */
 /** The forward, coarse-to-fine, transfer uses L2 projection. The backward,
@@ -207,6 +221,8 @@ public:
    protected:
       const FiniteElementSpace& fes_ho;
       const FiniteElementSpace& fes_lor;
+      CoefficientWithOrder coeff_ho;
+      CoefficientWithOrder coeff_lor;
 
       MemoryType d_mt;
       Array<int> offsets;
@@ -214,8 +230,15 @@ public:
 
       L2Projection(const FiniteElementSpace& fes_ho_,
                    const FiniteElementSpace& fes_lor_,
+                   CoefficientWithOrder coeff_ho_,
+                   CoefficientWithOrder coeff_lor_,
                    MemoryType d_mt_ = Device::GetHostMemoryType());
 
+      L2Projection(const FiniteElementSpace& fes_ho_,
+                   const FiniteElementSpace& fes_lor_,
+                   MemoryType d_mt_ = Device::GetHostMemoryType())
+         : L2Projection(fes_ho_, fes_lor_, nullptr, nullptr, d_mt_) { }
+
       void BuildHo2Lor(int nel_ho, int nel_lor,
                        const CoarseFineTransformations& cf_tr);
 
@@ -225,11 +248,11 @@ public:
                          IntegrationPointTransformation& ip_tr,
                          DenseMatrix& M_mixed_el) const;
 
-      void ElemMixedMass(Geometry::Type geom, const FiniteElement& fe_ho,
-                         const FiniteElement& fe_lor,
-                         ElementTransformation* el_tr,
-                         IntegrationPointTransformation& ip_tr,
-                         DenseMatrix& B_L, DenseMatrix& B_H) const;
+      void ElemMixedEvaluation(Geometry::Type geom, const FiniteElement& fe_ho,
+                               const FiniteElement& fe_lor,
+                               IntegrationPointTransformation& ip_tr,
+                               const IntegrationRule& ir,
+                               DenseMatrix& B_L, DenseMatrix& B_H) const;
    public:
       /* Returns the Mixed Mass M_LH via device element assembly by building the
       basis functions and data at the quadrature points. */
@@ -287,9 +310,17 @@ public:
    public:
       L2ProjectionL2Space(const FiniteElementSpace& fes_ho_,
                           const FiniteElementSpace& fes_lor_,
+                          CoefficientWithOrder coeff_ho_,
+                          CoefficientWithOrder coeff_lor_,
                           const bool use_ea_,
                           MemoryType d_mt_ = Device::GetHostMemoryType());
 
+      L2ProjectionL2Space(const FiniteElementSpace& fes_ho_,
+                          const FiniteElementSpace& fes_lor_,
+                          const bool use_ea_,
+                          MemoryType d_mt_ = Device::GetHostMemoryType())
+         : L2ProjectionL2Space(fes_ho_, fes_lor_, nullptr, nullptr, use_ea_, d_mt_) { }
+
       /*Same as above but assembles and stores R_ea, P_ea */
       void EAL2ProjectionL2Space();
 
@@ -356,13 +387,30 @@ public:
    public:
       L2ProjectionH1Space(const FiniteElementSpace &fes_ho_,
                           const FiniteElementSpace &fes_lor_,
+                          CoefficientWithOrder coeff_ho_,
+                          CoefficientWithOrder coeff_lor_,
                           const bool use_ea_,
                           MemoryType d_mt_ = Device::GetHostMemoryType());
+
+      L2ProjectionH1Space(const FiniteElementSpace& fes_ho_,
+                          const FiniteElementSpace& fes_lor_,
+                          const bool use_ea_,
+                          MemoryType d_mt_ = Device::GetHostMemoryType())
+         : L2ProjectionH1Space(fes_ho_, fes_lor_, nullptr, nullptr, use_ea_, d_mt_) { }
+
 #ifdef MFEM_USE_MPI
       L2ProjectionH1Space(const ParFiniteElementSpace &pfes_ho_,
                           const ParFiniteElementSpace &pfes_lor_,
+                          CoefficientWithOrder coeff_ho_,
+                          CoefficientWithOrder coeff_lor_,
                           const bool use_ea_,
                           MemoryType d_mt_ = Device::GetHostMemoryType());
+
+      L2ProjectionH1Space(const ParFiniteElementSpace& fes_ho_,
+                          const ParFiniteElementSpace& fes_lor_,
+                          const bool use_ea_,
+                          MemoryType d_mt_ = Device::GetHostMemoryType())
+         : L2ProjectionH1Space(fes_ho_, fes_lor_, nullptr, nullptr, use_ea_, d_mt_) { }
 #endif
       /// Same as above but assembles action of R through 4 parts:
       ///   ( )  inv( lumped(M_L) ), which is a diagonal matrix (essentially a vector)
@@ -508,18 +556,38 @@ public:
       virtual ~L2Prolongation() { }
    };
 
+   /// Coefficient for the mixed L2 inner product.
+   CoefficientWithOrder coeff_ho;
+   /// Coefficient for the low-order L2 inner product.
+   CoefficientWithOrder coeff_lor;
    L2Projection   *F; ///< Forward, coarse-to-fine, operator
    L2Prolongation *B; ///< Backward, fine-to-coarse, operator
    bool force_l2_space;
 
 public:
+   /// Construct the unweighted L2 projection grid transfer.
    L2ProjectionGridTransfer(FiniteElementSpace &coarse_fes_,
                             FiniteElementSpace &fine_fes_,
                             bool force_l2_space_ = false,
                             MemoryType d_mt_ = Device::GetHostMemoryType()) // move to method
       : GridTransfer(coarse_fes_, fine_fes_),
-        F(NULL), B(NULL), force_l2_space(force_l2_space_)
-   { }
+        coeff_ho(nullptr), coeff_lor(nullptr), F(nullptr), B(nullptr),
+        force_l2_space(force_l2_space_) { }
+
+   /// @brief Construct the weighted L2 projection grid transfer.
+   ///
+   /// The low-order inner product is weighted by @a coeff_lor, and the mixed
+   /// inner product is weighted by @a coeff_ho.
+   L2ProjectionGridTransfer(FiniteElementSpace &coarse_fes_,
+                            FiniteElementSpace &fine_fes_,
+                            CoefficientWithOrder coeff_ho_,
+                            CoefficientWithOrder coeff_lor_,
+                            bool force_l2_space_ = false,
+                            MemoryType d_mt_ = Device::GetHostMemoryType()) // move to method
+      : GridTransfer(coarse_fes_, fine_fes_),
+        coeff_ho(coeff_ho_), coeff_lor(coeff_lor_), F(nullptr), B(nullptr),
+        force_l2_space(force_l2_space_) { }
+
    virtual ~L2ProjectionGridTransfer();
 
    const Operator &ForwardOperator() override;
@@ -527,6 +595,7 @@ public:
    const Operator &BackwardOperator() override;
 
    bool SupportsBackwardsOperator() const override;
+
 private:
    void BuildF();
 };

general/CMakeLists.txt

@@ -17,6 +17,7 @@ list(APPEND SRCS
   error.cpp
   gecko.cpp
   globals.cpp
+  glvis_stream.cpp
   hash.cpp
   hash_util.cpp
   isockstream.cpp
@@ -45,6 +46,7 @@ list(APPEND HDRS
   error.hpp
   gecko.hpp
   globals.hpp
+  glvis_stream.hpp
   zstr.hpp
   hash.hpp
   hash_util.hpp

general/backends.hpp

@@ -14,7 +14,7 @@
 
 #include "../config/config.hpp"
 
-#if defined(MFEM_USE_CUDA) && defined(__CUDACC__)
+#if defined(MFEM_USE_CUDA)
 #include <cusparse.h>
 #include <library_types.h>
 #include <cuda_runtime.h>
@@ -22,7 +22,7 @@
 #endif
 #include "cuda.hpp"
 
-#if defined(MFEM_USE_HIP) && defined(__HIP__)
+#if defined(MFEM_USE_HIP)
 #include <hip/hip_runtime.h>
 #endif
 #include "hip.hpp"
@@ -45,15 +45,17 @@
 #endif
 
 #if !defined(MFEM_USE_CUDA_OR_HIP)
-constexpr bool mfem_use_gpu = false;
-#define MFEM_DEVICE
-#define MFEM_HOST
-#define MFEM_LAMBDA
-// #define MFEM_HOST_DEVICE // defined in config/config.hpp
 // MFEM_DEVICE_SYNC is made available for debugging purposes
 #define MFEM_DEVICE_SYNC
 // MFEM_STREAM_SYNC is used for UVM and MPI GPU-Aware kernels
 #define MFEM_STREAM_SYNC
+#endif
+
+#if !defined(MFEM_USE_CUDA_OR_HIP_LANG)
+#define MFEM_DEVICE
+#define MFEM_HOST
+#define MFEM_LAMBDA
+// #define MFEM_HOST_DEVICE // defined in config/config.hpp
 #define MFEM_LAUNCH_BOUNDS(...)
 #endif
 
@@ -126,4 +128,23 @@ MFEM_HOST_DEVICE T AtomicAdd(T &add, const T val)
 #endif
 }
 
+namespace mfem::internal
+{
+
+#if defined(MFEM_USE_CUDA_OR_HIP) && !defined(MFEM_USE_CUDA_OR_HIP_LANG)
+static constexpr bool can_compile_kernels = false;
+#else
+static constexpr bool can_compile_kernels = true;
+#endif
+
+template <bool can_compile_kernels = can_compile_kernels>
+void RequireKernelCompilation()
+{
+   static_assert(
+      can_compile_kernels,
+      "The calling function needs to be compiled with CUDA/HIP language!");
+}
+
+}
+
 #endif // MFEM_BACKENDS_HPP

general/cuda.hpp

@@ -18,14 +18,8 @@
 // CUDA block size used by MFEM.
 #define MFEM_CUDA_BLOCKS 256
 
-#if defined(MFEM_USE_CUDA) && defined(__CUDACC__)
+#if defined(MFEM_USE_CUDA)
 #define MFEM_USE_CUDA_OR_HIP
-constexpr bool mfem_use_gpu = true;
-#define MFEM_DEVICE __device__
-#define MFEM_HOST __host__
-#define MFEM_LAMBDA __host__
-#define MFEM_LAUNCH_BOUNDS __launch_bounds__
-// #define MFEM_HOST_DEVICE __host__ __device__ // defined in config/config.hpp
 #define MFEM_DEVICE_SYNC MFEM_GPU_CHECK(cudaDeviceSynchronize())
 #define MFEM_STREAM_SYNC MFEM_GPU_CHECK(cudaStreamSynchronize(0))
 // Define a CUDA error check macro, MFEM_GPU_CHECK(x), where x returns/is of
@@ -40,6 +34,15 @@ constexpr bool mfem_use_gpu = true;
     }                                                                          \
   } while (0)
 
+// Macros defined only when compiling with CUDA language
+#if defined(__CUDACC__)
+#define MFEM_USE_CUDA_OR_HIP_LANG
+#define MFEM_DEVICE __device__
+#define MFEM_HOST __host__
+#define MFEM_LAMBDA __host__
+#define MFEM_LAUNCH_BOUNDS __launch_bounds__
+// #define MFEM_HOST_DEVICE __host__ __device__ // defined in config/config.hpp
+
 // Define the MFEM inner threading macros
 #if defined(__CUDA_ARCH__)
 #define MFEM_SHARED __shared__
@@ -67,12 +70,13 @@ constexpr bool mfem_use_gpu = true;
    if (int ix = threadIdx.k % (OX), iy = threadIdx.k / (OX), iz = iy / (OY);   \
        (ix < (SX)) && ((iy %= (OY)) < (SY)) && (iz < (SZ)))
 #endif // defined(__CUDA_ARCH__)
-#endif // defined(MFEM_USE_CUDA) && defined(__CUDACC__)
+#endif // defined(__CUDACC__)
+#endif // defined(MFEM_USE_CUDA)
 
 namespace mfem
 {
 
-#if defined(MFEM_USE_CUDA) && defined(__CUDACC__)
+#if defined(MFEM_USE_CUDA)
 // Function used by the macro MFEM_GPU_CHECK.
 void mfem_cuda_error(cudaError_t err, const char *expr, const char *func,
                      const char *file, int line);

general/error.cpp

@@ -171,7 +171,7 @@ void mfem_error(const char *msg)
 #ifdef MFEM_USE_EXCEPTIONS
    if (mfem_error_action == MFEM_ERROR_THROW)
    {
-      throw ErrorException(msg);
+      throw ErrorException(msg ? msg : "");
    }
 #endif
 

general/error.hpp

@@ -15,7 +15,7 @@
 #include "../config/config.hpp"
 #include <iomanip>
 #include <sstream>
-#ifdef MFEM_USE_HIP
+#if defined(MFEM_USE_HIP)
 #include <hip/hip_runtime.h>
 #endif
 
@@ -153,21 +153,13 @@ void mfem_warning(const char *msg = NULL);
 
 
 // Additional abort functions for HIP
-#if defined(MFEM_USE_HIP)
-#ifndef __HIP_DEVICE_COMPILE__
-template<typename T>
-__host__ void abort_msg(T & msg)
-{
-   MFEM_ABORT(msg);
-}
-#else
+#if defined(__HIP_DEVICE_COMPILE__)
 template<typename T>
 __device__ void abort_msg(T & msg)
 {
    abort();
 }
 #endif
-#endif
 
 // Abort inside a device kernel
 #if defined(__CUDA_ARCH__)

general/forall.hpp

@@ -1044,6 +1044,8 @@ inline void ForallWrap(const bool use_dev, const int N,
                        const int X=0, const int Y=0, const int Z=0,
                        const int G=0)
 {
+   internal::RequireKernelCompilation();
+
    MFEM_CONTRACT_VAR(X);
    MFEM_CONTRACT_VAR(Y);
    MFEM_CONTRACT_VAR(Z);
@@ -1276,6 +1278,9 @@ inline void hypre_forall_cpu(int N, lambda &&body)
 template<typename lambda>
 inline void hypre_forall_gpu(int N, lambda &&body)
 {
+   internal::RequireKernelCompilation();
+
+#if defined(MFEM_USE_CUDA_OR_HIP_LANG)
 #if defined(HYPRE_USING_CUDA)
    CuWrap1D(N, body);
 #elif defined(HYPRE_USING_HIP)
@@ -1283,6 +1288,7 @@ inline void hypre_forall_gpu(int N, lambda &&body)
 #else
 #error Unknown HYPRE GPU backend!
 #endif
+#endif
 }
 #endif
 

general/glvis_stream.cpp

@@ -0,0 +1,188 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include <cassert>
+#include <istream>
+
+#include "../config/config.hpp"
+
+#ifdef MFEM_USE_GLVIS
+
+#include "glvis_stream.hpp"
+
+#ifdef MFEM_USE_MPI
+#include <mpi.h>
+#include <limits>
+#include <numeric>
+#endif
+
+#include "../fem/geom.hpp"
+thread_local mfem::GeometryRefiner GLVisGeometryRefiner;
+
+// Use local declaration to avoid circular dependency when fetching GLVis
+extern int GLVisStreamSession(
+   bool fix_elem_orient,
+   bool save_coloring,
+   bool keep_attr,
+   bool headless,
+   const std::string& plot_caption,
+   const std::string& data_type,
+   std::vector<std::unique_ptr<std::istream>>&& streams);
+
+namespace mfem
+{
+
+#ifdef MFEM_USE_MPI
+namespace
+{
+glvis_data MakeGlVisData()
+{
+   int size, rank;
+   MPI_Comm_size(MPI_COMM_WORLD, &size);
+   MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+   return glvis_data(size == 1, size, rank == 0);
+}
+} // namespace
+#endif
+
+glvis_stream::glvis_stream(): std::iostream(nullptr),
+#ifdef MFEM_USE_MPI
+   data(MakeGlVisData())
+#else
+   data(true, 1, true)
+#endif
+{
+   std::iostream::rdbuf(data.stream.rdbuf());
+}
+
+glvis_stream& glvis_stream::operator<<(ostream_manipulator pf)
+{
+   pf(static_cast<std::ostream&>(*this));
+   this->flush();
+   this->operator()();
+   return *this;
+}
+
+void glvis_stream::operator()()
+{
+   if (data.serial)
+   {
+      const auto size = this->size();
+      data.offsets.resize(2);
+      data.offsets[0] = 0, data.offsets[1] = size;
+      data.total_size = size;
+   }
+   else
+   {
+      serialize();
+   }
+
+   this->reset(); // reset the local buffer for reuse
+
+   if (data.mpi_root)
+   {
+      MFEM_VERIFY(data.mpi_size >= 0 &&
+                  (size_t) data.mpi_size == data.offsets.size() - 1,
+                  "Invalid MPI size");
+
+      data.streams.clear();
+      data.type.clear();
+
+      // loop over all input streams
+      for (int k = 0; k < data.mpi_size; ++k)
+      {
+         const size_t offset = data.offsets[k];
+         const size_t size = data.offsets[k+1] - data.offsets[k];
+
+         // add a new stream for this rank's data
+         data.streams.emplace_back(std::make_unique<std::stringstream>());
+         data.streams.back()->write(data.stream.str().data() + offset, size);
+
+         auto stream = data.streams.back().get();
+         if (!(*stream)) { break; }
+
+         *stream >> std::ws >> data.type >> std::ws;
+
+         if (data.type == "parallel") // Handle parallel data
+         {
+            int is_mpi_size, is_mpi_rank;
+            *stream >> is_mpi_size >> is_mpi_rank;
+            assert(is_mpi_size == static_cast<int>(data.mpi_size));
+            assert(is_mpi_rank == static_cast<int>(k));
+         }
+         else if (data.type != "mesh" && data.type != "solution")
+         {
+            MFEM_ABORT("Stream: unknown command: " << data.type);
+         }
+      }
+   }
+
+   if (!data.mpi_root) { return; }
+
+   constexpr bool fix_elem_orien = true;
+   constexpr bool save_coloring = true;
+   constexpr bool keep_attr = false;
+   constexpr bool headless = false;
+   const std::string plot_caption {};
+   std::vector<std::unique_ptr<std::istream>> istreams;
+   istreams.reserve(data.streams.size());
+   for (auto &s : data.streams) { istreams.push_back(std::move(s)); }
+   GLVisStreamSession(fix_elem_orien,
+                      save_coloring,
+                      keep_attr,
+                      headless,
+                      plot_caption,
+                      data.type,
+                      std::move(istreams));
+}
+
+void glvis_stream::serialize()
+{
+#ifdef MFEM_USE_MPI
+   const std::string local = data.stream.str();
+   MFEM_VERIFY(local.size() <= static_cast<size_t>
+               (std::numeric_limits<int>::max()),
+               "GLVis stream is too large for MPI_Gatherv");
+   const int local_size = static_cast<int>(local.size());
+
+   std::vector<int> sizes(data.mpi_size);
+   MFEM_VERIFY(MPI_Allgather(&local_size, 1, MPI_INT,
+                             sizes.data(), 1, MPI_INT,
+                             MPI_COMM_WORLD) == MPI_SUCCESS,
+               "MPI_Allgather failed");
+
+   if (data.mpi_root)
+   {
+      data.offsets.resize(data.mpi_size + 1);
+      data.offsets[0] = 0;
+      std::partial_sum(sizes.begin(), sizes.end(), data.offsets.begin() + 1);
+      data.total_size = data.offsets[data.mpi_size];
+   }
+
+   std::vector<char> recvbuf(data.mpi_root ? data.total_size : 0);
+   MFEM_VERIFY(MPI_Gatherv(local.data(), local_size, MPI_CHAR,
+                           data.mpi_root ? recvbuf.data() : nullptr,
+                           data.mpi_root ? sizes.data() : nullptr,
+                           data.mpi_root ? data.offsets.data() : nullptr,
+                           MPI_CHAR, 0, MPI_COMM_WORLD) == MPI_SUCCESS,
+               "MPI_Gatherv failed");
+
+   if (data.mpi_root)
+   {
+      reset();
+      data.stream.write(recvbuf.data(), data.total_size);
+   }
+#endif // MFEM_USE_MPI
+}
+
+} // namespace mfem
+
+#endif // MFEM_USE_GLVIS

general/glvis_stream.hpp

@@ -0,0 +1,83 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+#pragma once
+
+#include <memory>
+#include <sstream>
+#include <string>
+#include <vector>
+
+namespace mfem
+{
+
+struct glvis_data
+{
+   const bool serial;
+   const int mpi_size;
+   const bool mpi_root;
+   std::stringstream stream;
+   std::vector<std::unique_ptr<std::stringstream>> streams;
+   int total_size;
+   std::vector<int> offsets;
+   std::string type;
+
+   glvis_data(const bool serial, const int size, const bool root):
+      serial(serial), mpi_size(size), mpi_root(root),
+      total_size(0), type({}) {}
+};
+
+class glvis_stream: public std::iostream
+{
+   glvis_data data;
+   void serialize();
+
+public:
+   glvis_stream();
+
+   glvis_stream(glvis_stream &&) = delete;
+   glvis_stream(const glvis_stream &) = delete;
+   glvis_stream &operator=(const glvis_stream &) = delete;
+   glvis_stream &operator=(glvis_stream &&) = delete;
+
+   ~glvis_stream() = default;
+
+   size_t size() { return data.stream.tellp(); }
+   std::streamsize precision() const { return std::iostream::precision(); }
+   std::streamsize precision(std::streamsize new_prec)
+   { return std::iostream::precision(new_prec); }
+
+   using ostream_manipulator = std::ostream& (*)(std::ostream&);
+   glvis_stream& operator<<(ostream_manipulator pf);
+
+   template<typename T>
+   glvis_stream& operator<<(const T& val)
+   {
+      static_cast<std::ostream&>(*this) << val;
+      return *this;
+   }
+
+   int open(const char *, int) { return 0; }
+   bool is_open() const { return true; }
+   int close() { return 0; }
+
+   void flush() { std::iostream::flush(); }
+
+   void reset()
+   {
+      data.stream.clear();
+      data.stream.seekg(0, std::ios::beg);
+      data.stream.seekp(0, std::ios::beg);
+   }
+
+   void operator()();
+};
+
+} // namespace mfem

general/hip.hpp

@@ -18,14 +18,8 @@
 // HIP block size used by MFEM.
 #define MFEM_HIP_BLOCKS 256
 
-#if defined(MFEM_USE_HIP) && defined(__HIP__)
+#if defined(MFEM_USE_HIP)
 #define MFEM_USE_CUDA_OR_HIP
-constexpr bool mfem_use_gpu = true;
-#define MFEM_DEVICE __device__
-#define MFEM_HOST __host__
-#define MFEM_LAMBDA __host__ __device__
-#define MFEM_LAUNCH_BOUNDS __launch_bounds__
-// #define MFEM_HOST_DEVICE __host__ __device__ // defined in config/config.hpp
 #define MFEM_DEVICE_SYNC MFEM_GPU_CHECK(hipDeviceSynchronize())
 #define MFEM_STREAM_SYNC MFEM_GPU_CHECK(hipStreamSynchronize(0))
 // Define a HIP error check macro, MFEM_GPU_CHECK(x), where x returns/is of
@@ -40,6 +34,15 @@ constexpr bool mfem_use_gpu = true;
     }                                                                          \
   } while (0)
 
+// Macros defined only when compiling with HIP language
+#if defined(__HIP__)
+#define MFEM_USE_CUDA_OR_HIP_LANG
+#define MFEM_DEVICE __device__
+#define MFEM_HOST __host__
+#define MFEM_LAMBDA __host__ __device__
+#define MFEM_LAUNCH_BOUNDS __launch_bounds__
+// #define MFEM_HOST_DEVICE __host__ __device__ // defined in config/config.hpp
+
 // Define the MFEM inner threading macros
 #if defined(__HIP_DEVICE_COMPILE__)
 #define MFEM_SHARED __shared__
@@ -71,7 +74,8 @@ constexpr bool mfem_use_gpu = true;
        iz = iy / (OY);                                                         \
        (ix < (SX)) && ((iy %= (OY)) < (SY)) && (iz < (SZ)))
 #endif // defined(__HIP_DEVICE_COMPILE__)
-#endif // defined(MFEM_USE_HIP) && defined(__HIP__)
+#endif // defined(__HIP__)
+#endif // defined(MFEM_USE_HIP)
 
 namespace mfem
 {

general/reducers.hpp

@@ -550,10 +550,10 @@ void reduce(int N, T &res, B &&body, const R &reducer, bool use_dev,
 
       int num_mp = Device::NumMultiprocessors(Device::GetId());
 #if defined(MFEM_USE_CUDA)
-      // good value of mp_sat found experimentally on Lassen
+      // good value of mp_sat found experimentally on Lassen (V100)
       constexpr int mp_sat = 8;
 #elif defined(MFEM_USE_HIP)
-      // good value of mp_sat found experimentally on Tuolumne
+      // good value of mp_sat found experimentally on Tuolumne (MI300A)
       constexpr int mp_sat = 4;
 #else
       num_mp = 1;

general/scan.hpp

@@ -15,6 +15,10 @@
 #include "backends.hpp"
 #include "forall.hpp"
 
+#if defined(MFEM_USE_CUDA_OR_HIP) && !defined(MFEM_USE_CUDA_OR_HIP_LANG)
+#error "This header requires compilation with CUDA/HIP language!"
+#else
+
 #ifdef MFEM_USE_CUDA
 #include <cub/device/device_scan.cuh>
 #include <cub/device/device_select.cuh>
@@ -406,4 +410,6 @@ void CopyUnique(bool use_dev, InputIt d_in, OutputIt d_out,
 
 #undef MFEM_CUB_NAMESPACE
 
-#endif
+#endif // defined(MFEM_USE_CUDA_OR_HIP) && !defined(MFEM_USE_CUDA_OR_HIP_LANG)
+
+#endif // MFEM_SCAN_HPP

general/socketstream.cpp

@@ -1087,3 +1087,4 @@ socketstream::~socketstream()
 }
 
 } // namespace mfem
+

general/socketstream.hpp

@@ -13,6 +13,7 @@
 #define MFEM_SOCKETSTREAM
 
 #include "../config/config.hpp"
+
 #include "error.hpp"
 #include "globals.hpp"
 

general/version.cpp

@@ -100,6 +100,9 @@ const char *GetConfigStr()
 #ifdef MFEM_USE_GSLIB
       "MFEM_USE_GSLIB\n"
 #endif
+#ifdef MFEM_USE_GLVIS
+      "MFEM_USE_GLVIS\n"
+#endif
 #ifdef MFEM_USE_HDF5
       "MFEM_USE_HDF5\n"
 #endif

linalg/CMakeLists.txt

@@ -27,6 +27,7 @@ list(APPEND SRCS
   handle.cpp
   matrix.cpp
   mma.cpp
+  multivector.cpp
   ode.cpp
   operator.cpp
   ordering.cpp
@@ -63,6 +64,7 @@ list(APPEND HDRS
   linalg.hpp
   matrix.hpp
   mma.hpp
+  multivector.hpp
   ode.hpp
   operator.hpp
   ordering.hpp

linalg/densemat.hpp

@@ -1136,6 +1136,17 @@ private:
 public:
    DenseTensor() : ni(0), nj(0), nk(0) { }
 
+   DenseTensor(const DenseTensor &other)
+      : tdata(other.tdata), ni(other.ni), nj(other.nj), nk(other.nk) { }
+
+   DenseTensor(DenseTensor &&other)
+      : tdata(std::move(other.tdata)), ni(other.ni), nj(other.nj), nk(other.nk)
+   {
+      // Reset other; other.tdata is reset in Array<T> move constructror.
+      other.Mk.ClearExternalData();
+      other.ni = other.nj = other.nk = 0;
+   }
+
    DenseTensor(int i, int j, int k) : tdata(i*j*k), ni(i), nj(j), nk(k) { }
 
    DenseTensor(real_t *d, int i, int j, int k)
@@ -1144,6 +1155,33 @@ public:
    DenseTensor(int i, int j, int k, MemoryType mt)
       : tdata(i*j*k, mt), ni(i), nj(j), nk(k) { }
 
+   DenseTensor &operator=(const DenseTensor &other)
+   {
+      if (this == &other) { return *this; }
+      Mk.ClearExternalData();
+      tdata = other.tdata;
+      ni = other.ni;
+      nj = other.nj;
+      nk = other.nk;
+      return *this;
+   }
+
+   DenseTensor &operator=(DenseTensor &&other)
+   {
+      if (this == &other) { return *this; }
+      Mk.ClearExternalData();
+      tdata = std::move(other.tdata);
+      ni = other.ni;
+      nj = other.nj;
+      nk = other.nk;
+
+      // Reset other; other.tdata is reset in Array<T> move assignment.
+      other.Mk.ClearExternalData();
+      other.ni = other.nj = other.nk = 0;
+
+      return *this;
+   }
+
    int SizeI() const { return ni; }
    int SizeJ() const { return nj; }
    int SizeK() const { return nk; }

linalg/hypre.cpp

@@ -5842,6 +5842,10 @@ void HypreAMS::MakeGradientAndInterpolation(
    {
       grad->AddTraceFaceInterpolator(new GradientInterpolator);
    }
+   else if (dynamic_cast<const RT_FECollection *>(edge_fec))
+   {
+      grad->AddDomainInterpolator(new CurlInterpolator);
+   }
    else
    {
       grad->AddDomainInterpolator(new GradientInterpolator);

linalg/linalg.hpp

@@ -15,6 +15,7 @@
 // Linear algebra header file
 
 #include "vector.hpp"
+#include "multivector.hpp"
 #include "operator.hpp"
 #include "matrix.hpp"
 #include "sparsemat.hpp"

linalg/multivector.cpp

@@ -0,0 +1,60 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "multivector.hpp"
+
+namespace mfem
+{
+
+MultiVector::MultiVector(const Array<int> &vector_sizes)
+{
+   SetSizes(vector_sizes);
+}
+
+MultiVector::MultiVector(const Array<int> &vector_sizes, MemoryType mt)
+{
+   SetSizes(vector_sizes, mt);
+}
+
+MultiVector::MultiVector(Vector &base, const Array<int> &vector_sizes)
+{
+   MakeRef(base, vector_sizes);
+}
+
+void MultiVector::SetSizes(const Array<int> &vector_sizes)
+{
+   blocks.resize(vector_sizes.Size());
+   for (int i = 0; i < vector_sizes.Size(); i++)
+   {
+      operator[](i).SetSize(vector_sizes[i]);
+   }
+}
+
+void MultiVector::SetSizes(const Array<int> &vector_sizes, MemoryType mt)
+{
+   blocks.resize(vector_sizes.Size());
+   for (int i = 0; i < vector_sizes.Size(); i++)
+   {
+      operator[](i).SetSize(vector_sizes[i], mt);
+   }
+}
+
+void MultiVector::MakeRef(Vector &base, const Array<int> &vector_sizes)
+{
+   blocks.resize(vector_sizes.Size());
+   for (int offset = 0, i = 0; i < vector_sizes.Size(); i++)
+   {
+      blocks[i].emplace<0>(base, offset, vector_sizes[i]);
+      offset += vector_sizes[i];
+   }
+}
+
+} // namespace mfem

linalg/multivector.hpp

@@ -0,0 +1,251 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_MULTIVECTOR_HPP
+#define MFEM_MULTIVECTOR_HPP
+
+#include "../general/array.hpp"
+#include "vector.hpp"
+#include <vector>
+#include <array>
+#include <variant>
+
+namespace mfem
+{
+
+/// Class representing an array of Vectors with generally different sizes.
+/** This class is similar to BlockVector with the following two main
+    differences:
+    - the data for the individual Vector blocks does not need to be part of one
+      big contiguous memory allocation;
+    - this class does not inherit from class Vector (as a consequence of the
+      first bullet).
+
+    Internally, each Vector block is represented as one of the following
+    three options:
+    - (default) a Vector object constructed and owned by this class; this
+      object, in turn, as any Vector object, can own its Memory allocation or
+      refer to a sub-Memory of another Memory object; or
+    - a pointer to an externally allocated Vector or classes derived from
+      Vector.
+    - a pointer to an externally allocated const Vector or classes derived from
+      Vector. This option is helpful for wrapping const Vector objects as a
+      MultiVector that will be then used as a const MultiVector. */
+class MultiVector
+{
+private:
+   std::vector<std::variant<Vector,Vector*,const Vector*>> blocks;
+
+public:
+   /// Create an empty MultiVector with zero blocks.
+   MultiVector() = default;
+
+   /** @brief Create a MultiVector with @a num_blocks blocks. The individual
+       Vector blocks are default initialized, i.e. they all have size zero. */
+   MultiVector(int num_blocks)
+      : blocks(num_blocks) { }
+
+   /** @brief Construct a MultiVector with number of blocks and individual block
+       Vector sizes given by @a vector_sizes.
+
+       @note The memory of the individual Vector blocks is NOT initialized. */
+   MultiVector(const Array<int> &vector_sizes);
+
+   /** @brief Construct a MultiVector with number of blocks and individual block
+       Vector sizes given by @a vector_sizes. All Vector blocks use the
+       MemoryType @a mt.
+
+       @note The memory of the individual Vector blocks is NOT initialized. */
+   MultiVector(const Array<int> &vector_sizes, MemoryType mt);
+
+   /** @brief Construct a MultiVector referencing data within a given monolithic
+       Vector @a base.
+
+       With this constructor, the Memory flags of @a base and of the individual
+       Vector blocks may need to be explicitly synchronized when data is moved
+       between host and device. */
+   MultiVector(Vector &base, const Array<int> &vector_sizes);
+
+   /** @brief Construct a MultiVector referencing multiple Vectors given as
+       arguments.
+
+       The VectorTypes reference arguments are expected to be static_cast-able
+       to (Vector &) which is the case if the types are derived from Vector,
+       e.g. HypreParVector, GridFunction, etc.
+
+       With this constructor, operations on individual Vector blocks are
+       performed directly on the objects @a vs. In particular, there is no need
+       to synchronize the Memory flags of @a vs and the ones of the individual
+       Vector blocks when data is moved between host and device. */
+   template <typename... VectorTypes,
+             std::enable_if_t<
+                std::conjunction_v<
+                   std::is_convertible<VectorTypes&,Vector&>...>, bool> = true>
+   MultiVector(VectorTypes &...vs) { MakeRef(vs...); }
+
+   /** @brief Construct a MultiVector referencing multiple const Vectors given
+       as arguments. Individual blocks are read-only; non-const operator[]
+       will generate an error. */
+   template <typename... VectorTypes,
+             std::enable_if_t<
+                std::conjunction_v<
+                   std::is_convertible<const VectorTypes&,const Vector&>...>,
+                bool> = true>
+   MultiVector(const VectorTypes &...vs) { MakeRef(vs...); }
+
+   /// Return the number of Vectors in the MultiVector.
+   int NumBlocks() const { return blocks.size(); }
+
+   /** @brief Set the number of Vectors in the MultiVector. Existing Vector
+       blocks will remain unmodified. New Vector blocks will be default
+       initialized, i.e. they all have size zero. */
+   void SetNumBlocks(int num_blocks) { blocks.resize(num_blocks); }
+
+   /** @brief Read-write access to the i-th Vector. Generates an error if the
+       i-th block is read-only, i.e. it is a pointer to a const Vector. */
+   inline Vector &operator[](int i);
+
+   /// Read-only access to the i-th Vector.
+   inline const Vector &operator[](int i) const;
+
+   /** @brief Update the MultiVector according to the given @a vector_sizes.
+
+       This method can be used to add or remove blocks. The individual Vector
+       sizes are updated using the method Vector::SetSize(int). */
+   void SetSizes(const Array<int> &vector_sizes);
+
+   /** @brief Update the MultiVector according to the given @a vector_sizes and
+       MemoryType @a mt.
+
+       This method can be used to add or remove blocks. The individual Vector
+       sizes and MemoryType are updated using the method
+       Vector::SetSize(int, MemoryType). */
+   void SetSizes(const Array<int> &vector_sizes, MemoryType mt);
+
+   /** @brief Update the MultiVector to reference data within a given monolithic
+       Vector @a base.
+
+       After calling this method, the Memory flags of @a base and of the
+       individual Vector blocks may need to be explicitly synchronized when data
+       is moved between host and device.*/
+   void MakeRef(Vector &base, const Array<int> &vector_sizes);
+
+   /** @brief Update the @a i-th MultiVector block to reference data within the
+       given monolithic Vector @a base at the given @a offset and with the given
+       @a size.
+
+       After calling this method, the Memory flags of @a base and of the @a i-th
+       Vector block may need to be explicitly synchronized when data is moved
+       between host and device.*/
+   inline void MakeRef(int i, Vector &base, int offset, int size)
+   {
+      blocks[i].emplace<0>(base, offset, size);
+   }
+
+   /** @brief Update the MultiVector to reference multiple Vectors given as
+       arguments.
+
+       The VectorTypes reference arguments are expected to be static_cast-able
+       to (Vector &) which is the case if the types are derived from Vector,
+       e.g. HypreParVector, GridFunction, etc.
+
+       After calling this method, operations on individual Vector blocks are
+       performed directly on the objects @a vs. In particular, there is no need
+       to synchronize the Memory flags of @a vs and the ones of the individual
+       Vector blocks when data is moved between host and device. */
+   template <typename... VectorTypes,
+             std::enable_if_t<
+                std::conjunction_v<
+                   std::is_convertible<VectorTypes&,Vector&>...>, bool> = true>
+   inline void MakeRef(VectorTypes &...vs);
+
+   /** @brief Update the MultiVector to reference multiple const Vectors given
+       as arguments. Individual blocks are read-only; non-const operator[]
+       will generate an error. */
+   template <typename... VectorTypes,
+             std::enable_if_t<
+                std::conjunction_v<
+                   std::is_convertible<const VectorTypes&,const Vector&>...>,
+                bool> = true>
+   inline void MakeRef(const VectorTypes &...vs);
+
+   /** @brief Update the @a i-th MultiVector block to reference the given
+       Vector @a v.
+
+       After calling this method, operations on the @a i-th Vector block are
+       performed directly on the Vector @a v. In particular, there is no need
+       to synchronize the Memory flags of @a v and the ones of the @a i-th
+       Vector blocks when data is moved between host and device. */
+   inline void MakeRef(int i, Vector &v) { blocks[i] = &v; }
+
+   /** @brief Update the @a i-th MultiVector block to reference the given
+       const Vector @a v. The block becomes read-only. */
+   inline void MakeRef(int i, const Vector &v) { blocks[i] = &v; }
+};
+
+// Inline and template methods
+
+inline Vector &MultiVector::operator[](int i)
+{
+   auto &bi = blocks[i];
+   const auto idx = bi.index();
+   if (idx == 0) { return std::get<0>(bi); }
+   if (idx == 1) { return *std::get<1>(bi); }
+   MFEM_ABORT("Non-const access to a const Vector block!");
+}
+
+inline const Vector &MultiVector::operator[](int i) const
+{
+   auto &bi = blocks[i];
+   const auto idx = bi.index();
+   return (idx == 0) ? std::get<0>(bi) :
+          (idx == 1) ? *std::get<1>(bi) :
+          /**/         *std::get<2>(bi);
+}
+
+template <typename... VectorTypes,
+          std::enable_if_t<
+             std::conjunction_v<
+                std::is_convertible<VectorTypes&,Vector&>...>, bool>>
+inline void MultiVector::MakeRef(VectorTypes &...vs)
+{
+   blocks.resize(sizeof...(vs));
+   if constexpr (sizeof...(vs) > 0)
+   {
+      const std::array vs_p{&static_cast<Vector&>(vs)...};
+      for (std::size_t i = 0; i < sizeof...(vs); i++)
+      {
+         blocks[i] = vs_p[i];
+      }
+   }
+}
+
+template <typename... VectorTypes,
+          std::enable_if_t<
+             std::conjunction_v<
+                std::is_convertible<const VectorTypes&,const Vector&>...>,
+             bool>>
+inline void MultiVector::MakeRef(const VectorTypes &...vs)
+{
+   blocks.resize(sizeof...(vs));
+   if constexpr (sizeof...(vs) > 0)
+   {
+      const std::array vs_p{&static_cast<const Vector&>(vs)...};
+      for (std::size_t i = 0; i < sizeof...(vs); i++)
+      {
+         blocks[i] = vs_p[i];
+      }
+   }
+}
+
+} // namespace mfem
+
+#endif // MFEM_MULTIVECTOR_HPP

linalg/operator.cpp

@@ -111,6 +111,21 @@ void Operator::ArrayAddMultTranspose(const Array<const Vector *> &X,
    }
 }
 
+void Operator::MultMV(const MultiVector &, MultiVector &) const
+{
+   MFEM_ABORT("this method is not overridden for this class!");
+}
+
+void Operator::MultTransposeMV(const MultiVector &x, MultiVector &y) const
+{
+   MFEM_ABORT("this method is not overridden for this class!");
+}
+
+Operator &Operator::GetGradientMV(const MultiVector &) const
+{
+   MFEM_ABORT("this method is not overridden for this class!");
+}
+
 void Operator::FormLinearSystem(const Array<int> &ess_tdof_list,
                                 Vector &x, Vector &b,
                                 Operator* &Aout, Vector &X, Vector &B,

linalg/operator.hpp

@@ -13,6 +13,7 @@
 #define MFEM_OPERATOR
 
 #include "vector.hpp"
+#include "multivector.hpp"
 
 namespace mfem
 {
@@ -129,6 +130,20 @@ public:
    virtual void ArrayAddMultTranspose(const Array<const Vector *> &X,
                                       Array<Vector *> &Y, const real_t a = 1.0) const;
 
+   /** @brief Operator application, y = A(x), where the input @a x and the
+       output @a y are MultiVector objects, i.e. they generally use
+       non-contiguous memory representation.
+
+       The base class implementation for the method is to generate an error. */
+   virtual void MultMV(const MultiVector &x, MultiVector &y) const;
+
+   /** @brief Action of the transpose operator, y = A^t(x), where the input @a x
+       and the output @a y are MultiVector objects, i.e. they generally use
+       non-contiguous memory representation.
+
+       The base class implementation for this method is to generate an error. */
+   virtual void MultTransposeMV(const MultiVector &x, MultiVector &y) const;
+
    /** @brief Evaluate the gradient operator at the point @a x. The default
        behavior in class Operator is to generate an error. */
    virtual Operator &GetGradient(const Vector &x) const
@@ -137,6 +152,13 @@ public:
       return const_cast<Operator &>(*this);
    }
 
+   /** @brief Evaluate the gradient operator at the point @a x. The input @a x
+       is provided as a MultiVector, i.e. it generally uses non-contiguous
+       memory representation.
+
+       The base class implementation for the method is to generate an error. */
+   virtual Operator &GetGradientMV(const MultiVector &x) const;
+
    /** @brief Computes the diagonal entries into @a diag. Typically, this
        operation only makes sense for linear Operator%s. In some cases, only an
        approximation of the diagonal is computed. */

makefile

@@ -119,7 +119,8 @@ $(if $(word 2,$(SRC)),$(error Spaces in SRC = "$(SRC)" are not supported))
 MFEM_GIT_STRING = $(shell [ -d $(MFEM_DIR)/.git ] && git -C $(MFEM_DIR) \
    describe --all --long --abbrev=40 --dirty --always 2> /dev/null)
 
-EXAMPLE_SUBDIRS = amgx caliper ginkgo hiop petsc pumi sundials superlu moonolith
+EXAMPLE_SUBDIRS = amgx caliper ginkgo glvis hiop petsc pumi sundials \
+ superlu moonolith
 EXAMPLE_DIRS := examples $(addprefix examples/,$(EXAMPLE_SUBDIRS))
 EXAMPLE_TEST_DIRS := examples
 
@@ -163,6 +164,8 @@ MFEM_BUILD_DIR := $(BUILD_DIR)
 
 CONFIG_MK = $(BLD)config/config.mk
 
+GLVIS_MK = $(SRC)config/glvis.mk
+
 DEFAULTS_MK = $(SRC)config/defaults.mk
 include $(DEFAULTS_MK)
 
@@ -291,6 +294,11 @@ ifeq ($(MFEM_USE_HIP),YES)
    endif
 endif
 
+# GLVis configuration
+ifeq ($(MFEM_USE_GLVIS),YES)
+	GLVIS_DIR:=$(abspath $(subst @MFEM_DIR@,$(if $(MFEM_DIR),$(MFEM_DIR),..),$(GLVIS_DIR)))
+endif
+
 DEP_CXX ?= $(MFEM_CXX)
 
 # Check legacy OpenMP configuration
@@ -307,7 +315,7 @@ endif
 MFEM_REQ_LIB_DEPS = SUPERLU MUMPS METIS FMS CONDUIT SIDRE LAPACK SUNDIALS\
  SUITESPARSE STRUMPACK GINKGO GNUTLS HDF5 NETCDF SLEPC PETSC MPFR PUMI HIOP\
  GSLIB OCCA CEED RAJA UMPIRE MKL_CPARDISO MKL_PARDISO AMGX MAGMA CALIPER PARELAG\
- TRIBOL BENCHMARK MOONOLITH ALGOIM CUDSS
+ TRIBOL BENCHMARK MOONOLITH ALGOIM CUDSS GLVIS
 
 
 PETSC_ERROR_MSG = $(if $(PETSC_FOUND),,. PETSC config not found: $(PETSC_VARS))
@@ -377,7 +385,8 @@ MFEM_DEFINES = MFEM_VERSION MFEM_VERSION_STRING MFEM_GIT_STRING MFEM_USE_MPI\
  MFEM_USE_MAGMA MFEM_USE_MUMPS MFEM_USE_ADFORWARD MFEM_USE_CODIPACK MFEM_USE_CALIPER\
  MFEM_USE_BENCHMARK MFEM_USE_PARELAG MFEM_USE_TRIBOL MFEM_USE_ALGOIM MFEM_USE_ENZYME\
  MFEM_SOURCE_DIR MFEM_INSTALL_DIR MFEM_SHARED_BUILD MFEM_USE_DOUBLE MFEM_USE_SINGLE\
- MFEM_USE_CUDSS MFEM_CUDSS_COMM_LIB MFEM_CUDSS_THREADING_LIB
+ MFEM_USE_CUDSS MFEM_CUDSS_COMM_LIB MFEM_CUDSS_THREADING_LIB\
+ MFEM_USE_GLVIS
 
 # List of makefile variables that will be written to config.mk:
 MFEM_CONFIG_VARS = MFEM_CXX MFEM_HOST_CXX MFEM_CPPFLAGS MFEM_CXXFLAGS\
@@ -477,7 +486,8 @@ OKL_DIRS = fem
 %:	%.cpp
 
 # Default rule.
-lib: $(if $(static),$(BLD)libmfem.a) $(if $(shared),$(BLD)libmfem.$(SO_EXT))
+lib: $(if $(static),$(BLD)libmfem.a) $(if $(shared),$(BLD)libmfem.$(SO_EXT)) \
+	$(if $(filter YES,$(MFEM_USE_GLVIS)),$(if $(static), $(GLVIS_DIR)/lib/libglvis.a))
 
 # Flags used for compiling all source files.
 MFEM_BUILD_FLAGS = $(MFEM_PICFLAG) $(MFEM_CPPFLAGS) $(MFEM_CXXFLAGS)\
@@ -510,6 +520,14 @@ $(BLD)libmfem.$(SO_EXT): $(BLD)libmfem.$(SO_VER)
 	cd $(@D) && ln -sf $(<F) $(@F)
 	@$(MAKE) deprecation-warnings
 
+ifeq ($(MFEM_USE_GLVIS),YES)
+$(GLVIS_DIR)/lib/libglvis.a: $(BLD)libmfem.a
+	$(if $(wildcard $(GLVIS_DIR)/makefile),,$(error No makefile in GLVIS_DIR: $(GLVIS_DIR)))
+	@$(MAKE) -C $(GLVIS_DIR) -j $(shell getconf _NPROCESSORS_ONLN 2>/dev/null || echo 1) \
+		MFEM_DIR=$(BUILD_REAL_DIR) \
+		GLVIS_USE_LOGO=NO GLVIS_USE_LIBPNG=YES lib/libglvis.a
+endif
+
 # If some of the external libraries are build without -fPIC, linking shared MFEM
 # library may fail. In such cases, one may set EXT_LIBS on the command line.
 EXT_LIBS = $(MFEM_EXT_LIBS)
@@ -604,6 +622,7 @@ clean: $(addsuffix /clean,$(EM_DIRS) $(TEST_DIRS))
 
 distclean: clean config/clean doc/clean
 	rm -rf mfem/
+	$(if $(filter YES,$(MFEM_USE_GLVIS)),-$(MAKE) -C $(GLVIS_DIR) distclean)
 
 # User-definable install permissions.
 # Install permissions for everything except directories and binaries:
@@ -628,10 +647,14 @@ INSTALL_SHARED_LIB = $(MFEM_CXX) $(MFEM_LINK_FLAGS) $(INSTALL_SOFLAGS)\
    cd $(PREFIX_LIB) && chmod $(INSTALL_BIN_PERM) libmfem.$(SO_VER) && \
    ( umask $(INSTALLMASK) && ln -sf libmfem.$(SO_VER) libmfem.$(SO_EXT) )
 
-install: $(if $(static),$(BLD)libmfem.a) $(if $(shared),$(BLD)libmfem.$(SO_EXT))
+install: $(if $(static),$(BLD)libmfem.a) \
+			$(if $(shared),$(BLD)libmfem.$(SO_EXT)) \
+			$(if $(filter YES,$(MFEM_USE_GLVIS)),\
+				$(if $(static),$(GLVIS_DIR)/lib/libglvis.a))
 	$(MKINSTALLDIR) $(PREFIX_LIB)
 # install static and/or shared library
 	$(if $(static),$(INSTALLDEF) $(BLD)libmfem.a $(PREFIX_LIB))
+	$(if $(filter YES,$(MFEM_USE_GLVIS)),$(if $(static),$(INSTALLDEF) $(GLVIS_DIR)/lib/libglvis.a) $(PREFIX_LIB))
 	$(if $(shared),$(INSTALL_SHARED_LIB))
 # install top level includes
 	$(MKINSTALLDIR) $(PREFIX_INC)/mfem
@@ -778,6 +801,7 @@ status info:
 	$(info MFEM_USE_PARELAG       = $(MFEM_USE_PARELAG))
 	$(info MFEM_USE_TRIBOL        = $(MFEM_USE_TRIBOL))
 	$(info MFEM_USE_ENZYME        = $(MFEM_USE_ENZYME))
+	$(info MFEM_USE_GLVIS         = $(MFEM_USE_GLVIS))
 	$(info MFEM_CXX               = $(value MFEM_CXX))
 	$(info MFEM_HOST_CXX          = $(value MFEM_HOST_CXX))
 	$(info MFEM_CPPFLAGS          = $(value MFEM_CPPFLAGS))

mesh/pncmesh.cpp

@@ -2624,7 +2624,8 @@ void ParNCMesh::RedistributeElements(Array<int> &new_ranks, int target_elements,
          for (int i = 0; i < rank_neighbors.Size(); i++)
          {
             int elem = rank_neighbors[i];
-            msg.AddElementRank(elem, new_ranks[elements[elem].index]);
+            const Element &el = elements[elem];
+            msg.AddElement(elem, new_ranks[el.index], el.attribute);
          }
 
          msg.Isend(rank, MyComm);
@@ -2647,7 +2648,9 @@ void ParNCMesh::RedistributeElements(Array<int> &new_ranks, int target_elements,
       {
          int ghost_index = elements[msg.elements[i]].index;
          MFEM_ASSERT(element_type[ghost_index] == 2, "");
-         new_ranks[ghost_index] = msg.values[i];
+         const ElementRankAndAttribute &value = msg.values[i];
+         new_ranks[ghost_index] = value.rank;
+         elements[msg.elements[i]].attribute = value.attribute;
       }
    }
 
@@ -2718,7 +2721,7 @@ void ParNCMesh::RedistributeElements(Array<int> &new_ranks, int target_elements,
 
                if ((element_type[el.index] & 1) || el.rank != rank)
                {
-                  msg.AddElementRank(elem, el.rank);
+                  msg.AddElement(elem, el.rank, el.attribute);
                }
                // NOTE: we skip 'ghosts' that are of the receiver's rank because
                // they are not really ghosts and would get sent multiple times,
@@ -2770,10 +2773,12 @@ void ParNCMesh::RedistributeElements(Array<int> &new_ranks, int target_elements,
 
          for (int i = 0; i < msg.Size(); i++)
          {
-            int elem_rank = msg.values[i];
-            elements[msg.elements[i]].rank = elem_rank;
+            const ElementRankAndAttribute &value = msg.values[i];
+            Element &el = elements[msg.elements[i]];
+            el.rank = value.rank;
+            el.attribute = value.attribute;
 
-            if (elem_rank == MyRank) { received_elements++; }
+            if (value.rank == MyRank) { received_elements++; }
          }
 
          // save the ranks we received from, for later use in RecvRebalanceDofs
@@ -2809,7 +2814,10 @@ void ParNCMesh::RedistributeElements(Array<int> &new_ranks, int target_elements,
 
             for (int i = 0; i < msg.Size(); i++)
             {
-               elements[msg.elements[i]].rank = msg.values[i];
+               const ElementRankAndAttribute &value = msg.values[i];
+               Element &el = elements[msg.elements[i]];
+               el.rank = value.rank;
+               el.attribute = value.attribute;
             }
 
             // save the ranks we received from, for later use in RecvRebalanceDofs

mesh/pncmesh.hpp

@@ -531,26 +531,36 @@ protected: // implementation
       typedef std::map<int, NeighborDerefinementMessage> Map;
    };
 
-   /** Used in Step 2 of Rebalance() to synchronize new rank assignments in
-    *  the ghost layer.
+   struct ElementRankAndAttribute
+   {
+      int rank;
+      int attribute;
+   };
+
+   /** Used in RedistributeElements() to synchronize new rank assignments and
+    *  element attributes in the ghost layer.
     */
-   class NeighborElementRankMessage : public ElementValueMessage<int, false,
+   class NeighborElementRankMessage :
+      public ElementValueMessage<ElementRankAndAttribute, false,
       VarMessageTag::NEIGHBOR_ELEMENT_RANK_VM>
    {
    public:
-      void AddElementRank(int elem, int rank) { Add(elem, rank); }
+      void AddElement(int elem, int rank, int attribute)
+      { Add(elem, {rank, attribute}); }
       typedef std::map<int, NeighborElementRankMessage> Map;
    };
 
-   /** Used by Rebalance() to send elements and their ranks. Note that
+   /** Used by Rebalance() to send elements, ranks, and attributes. Note that
     *  RefTypes == true which means the refinement hierarchy will be recreated
     *  on the receiving side.
     */
-   class RebalanceMessage : public ElementValueMessage<int, true,
+   class RebalanceMessage :
+      public ElementValueMessage<ElementRankAndAttribute, true,
       VarMessageTag::REBALANCE_VM>
    {
    public:
-      void AddElementRank(int elem, int rank) { Add(elem, rank); }
+      void AddElement(int elem, int rank, int attribute)
+      { Add(elem, {rank, attribute}); }
       typedef std::map<int, RebalanceMessage> Map;
    };
 

mfem.hpp

@@ -30,6 +30,9 @@
 #ifdef MFEM_USE_ADIOS2
 #include "general/adios2stream.hpp"
 #endif // MFEM_USE_ADIOS2
+#ifdef MFEM_USE_GLVIS
+#include "general/glvis_stream.hpp"
+#endif // MFEM_USE_GLVIS
 #include "general/isockstream.hpp"
 #include "general/osockstream.hpp"
 #include "general/socketstream.hpp"

miniapps/fluids/navier/makefile

@@ -52,6 +52,8 @@ endif
 .SUFFIXES:
 .SUFFIXES: .o .cpp .mk
 .PHONY: all lib-common clean clean-build clean-exec
+# Keeping the *.o files fixes an issue with the MacOS version of 'make'.
+.PRECIOUS: %.o
 
 # Remove built-in rules
 %: %.cpp

miniapps/multidomain/makefile

@@ -68,7 +68,7 @@ multidomain-test-par: multidomain
 multidomain_nd-test-par: multidomain_nd
 	@$(call mfem-test,$<, $(RUN_MPI), Multidomain ND miniapp,-tf 0.001)
 multidomain_rt-test-par: multidomain_rt
-	@$(call mfem-test,$<, $(RUN_MPI), Multidomain RT iniapp,-tf 0.001)
+	@$(call mfem-test,$<, $(RUN_MPI), Multidomain RT miniapp,-tf 0.001)
 
 # Generate an error message if the MFEM library is not built and exit
 $(MFEM_LIB_FILE):

miniapps/parelag/MultilevelHcurlHdivSolver.cpp

@@ -761,7 +761,7 @@ int main(int argc, char *argv[])
       if (visualize)
       {
          hcurlhdiv_dofTrueDof.Distribute(X, x);
-         MultiVector tmp(x.GetData(), 1, x.Size());
+         parelag::MultiVector tmp(x.GetData(), 1, x.Size());
          sequence[0]->show(jform, tmp);
       }
       post_timer.Stop();

miniapps/tools/lor-transfer.cpp

@@ -33,6 +33,7 @@
 //
 // Sample runs:  lor-transfer
 //               lor-transfer -h1
+//               lor-transfer -ea -w
 //               lor-transfer -t
 //               lor-transfer -m ../../data/star-q2.mesh -lref 5 -p 4
 //               lor-transfer -m ../../data/star-mixed.mesh -lref 3 -p 2
@@ -59,11 +60,12 @@ string direction;
 
 // Exact functions to project
 real_t RHO_exact(const Vector &x);
+real_t W_exact(const Vector &x);
+real_t weight(const Vector &x);
 
 // Helper functions
 void visualize(VisItDataCollection &, string, int, int, int visport = 19916);
-real_t compute_mass(FiniteElementSpace *, real_t, VisItDataCollection &,
-                    string);
+real_t compute_mass(GridFunction &, real_t, string, CoefficientWithOrder);
 
 int main(int argc, char *argv[])
 {
@@ -76,6 +78,7 @@ int main(int argc, char *argv[])
    bool useH1 = false;
    int visport = 19916;
    bool use_pointwise_transfer = false;
+   bool use_weighted_transfer = false;
    const char *device_config = "cpu";
    bool use_ea       = false;
 
@@ -98,6 +101,9 @@ int main(int argc, char *argv[])
    args.AddOption(&use_pointwise_transfer, "-t", "--use-pointwise-transfer",
                   "-no-t", "--dont-use-pointwise-transfer",
                   "Use pointwise transfer operators instead of L2 projection.");
+   args.AddOption(&use_weighted_transfer, "-w", "--use-weighted-transfer",
+                  "-no-w", "--dont-use-weighted-transfer",
+                  "Use coefficient-weighted L2 projection.");
    args.AddOption(&device_config, "-d", "--device",
                   "Device configuration string, see Device::Configure().");
    args.AddOption(&use_ea, "-ea", "--ea-version", "-no-ea",
@@ -107,6 +113,15 @@ int main(int argc, char *argv[])
    // Configure device
    Device device(device_config);
 
+   if (use_weighted_transfer && !use_pointwise_transfer)
+   {
+      if (problem != 5)
+      {
+         cout << "Switching to positive problem = 5 for weighted transfer.\n";
+      }
+      problem = 5;
+   }
+
    // Read the mesh from the given mesh file.
    Mesh mesh(mesh_file, 1, 1);
    int dim = mesh.Dimension();
@@ -138,6 +153,14 @@ int main(int argc, char *argv[])
    FiniteElementSpace fespace(&mesh, fec);
    FiniteElementSpace fespace_lor(&mesh_lor, fec_lor);
 
+   FunctionCoefficient weight_fn_coeff(weight);
+   CoefficientWithOrder weight_coeff;
+   if (use_weighted_transfer)
+   {
+      weight_coeff.coeff = &weight_fn_coeff;
+      weight_coeff.order = 2;
+   }
+
    GridFunction rho(&fespace);
    GridFunction rho_lor(&fespace_lor);
 
@@ -165,7 +188,7 @@ int main(int argc, char *argv[])
    rho.SetTrueVector();
    rho.SetFromTrueVector();
 
-   real_t ho_mass = compute_mass(&fespace, -1.0, HO_dc, "HO       ");
+   real_t ho_mass = compute_mass(rho, -1.0, "HO       ", weight_coeff);
    if (vis) { visualize(HO_dc, "HO", Wx, Wy, visport); Wx += offx; }
 
    GridTransfer *gt;
@@ -175,7 +198,8 @@ int main(int argc, char *argv[])
    }
    else
    {
-      gt = new L2ProjectionGridTransfer(fespace, fespace_lor);
+      gt = new L2ProjectionGridTransfer(fespace, fespace_lor, weight_coeff,
+                                        weight_coeff);
    }
 
    // Configure element assembly for device acceleration
@@ -186,9 +210,44 @@ int main(int argc, char *argv[])
    // HO->LOR restriction
    direction = "HO -> LOR @ LOR";
    R.Mult(rho, rho_lor);
-   compute_mass(&fespace_lor, ho_mass, LOR_dc, "R(HO)    ");
+   compute_mass(rho_lor, ho_mass, "R(HO)    ", weight_coeff);
    if (vis) { visualize(LOR_dc, "R(HO)", Wx, Wy, visport); Wx += offx; }
 
+   if (use_weighted_transfer && !use_pointwise_transfer)
+   {
+      // Transfer velocity while conserving rho-weighted momentum.
+      GridFunctionCoefficient rho_coeff(&rho);
+      GridFunctionCoefficient rho_lor_coeff(&rho_lor);
+      ProductCoefficient prod_coeff(weight_fn_coeff, rho_coeff);
+      ProductCoefficient prod_lor_coeff(weight_fn_coeff, rho_lor_coeff);
+      CoefficientWithOrder prod_weight(prod_coeff, order + 2);
+      CoefficientWithOrder prod_lor_weight(prod_lor_coeff, lorder + 2);
+
+      GridFunction w(&fespace), w_lor(&fespace_lor);
+      FunctionCoefficient W(W_exact);
+      w.ProjectCoefficient(W);
+
+      cout << '\n';
+      const real_t ho_momentum = compute_mass(w, -1.0, "rho w HO ", prod_weight);
+
+      L2ProjectionGridTransfer vel_gt(fespace, fespace_lor, prod_weight,
+                                      prod_lor_weight);
+      vel_gt.UseEA(use_ea);
+      vel_gt.ForwardOperator().Mult(w, w_lor);
+      compute_mass(w_lor, ho_momentum, "rho w LOR", prod_lor_weight);
+
+      if (vel_gt.SupportsBackwardsOperator())
+      {
+         GridFunction w_prev = w;
+         vel_gt.BackwardOperator().Mult(w_lor, w);
+         compute_mass(w, ho_momentum, "P(rho w) ", prod_weight);
+
+         w_prev -= w;
+         cout.precision(12);
+         cout << "|w - P(R(w))|_∞     = " << w_prev.Normlinf() << "\n\n";
+      }
+   }
+
    if (gt->SupportsBackwardsOperator())
    {
       const Operator &P = gt->BackwardOperator();
@@ -196,7 +255,7 @@ int main(int argc, char *argv[])
       direction = "HO -> LOR @ HO";
       GridFunction rho_prev = rho;
       P.Mult(rho_lor, rho);
-      compute_mass(&fespace, ho_mass, HO_dc, "P(R(HO)) ");
+      compute_mass(rho, ho_mass, "P(R(HO)) ", weight_coeff);
       if (vis) { visualize(HO_dc, "P(R(HO))", Wx, Wy, visport); Wx = 0; Wy += offy; }
 
       rho_prev -= rho;
@@ -218,7 +277,7 @@ int main(int argc, char *argv[])
    direction = "LOR -> HO @ LOR";
    rho_lor.ProjectCoefficient(RHO);
    GridFunction rho_lor_prev = rho_lor;
-   real_t lor_mass = compute_mass(&fespace_lor, -1.0, LOR_dc, "LOR      ");
+   real_t lor_mass = compute_mass(rho_lor, -1.0, "LOR      ", weight_coeff);
    if (vis) { visualize(LOR_dc, "LOR", Wx, Wy, visport); Wx += offx; }
 
    if (gt->SupportsBackwardsOperator())
@@ -227,14 +286,14 @@ int main(int argc, char *argv[])
       // Prolongate to HO space
       direction = "LOR -> HO @ HO";
       P.Mult(rho_lor, rho);
-      compute_mass(&fespace, lor_mass, HO_dc, "P(LOR)   ");
+      compute_mass(rho, lor_mass, "P(LOR)   ", weight_coeff);
       if (vis) { visualize(HO_dc, "P(LOR)", Wx, Wy, visport); Wx += offx; }
 
       // Restrict back to LOR space. This won't give the original function because
       // the rho_lor doesn't necessarily live in the range of R.
       direction = "LOR -> HO @ LOR";
       R.Mult(rho, rho_lor);
-      compute_mass(&fespace_lor, lor_mass, LOR_dc, "R(P(LOR))");
+      compute_mass(rho_lor, lor_mass, "R(P(LOR))", weight_coeff);
       if (vis) { visualize(LOR_dc, "R(P(LOR))", Wx, Wy, visport); }
 
       rho_lor_prev -= rho_lor;
@@ -270,12 +329,26 @@ real_t RHO_exact(const Vector &x)
          return M_PI/2-atan(5*(2*x.Norml2()-1));
       case 4: // basis function
          return (x.Norml2() < 0.1) ? 1 : 0;
+      case 5: // positive function
+         return 2.0 + 2*x(0)*x(0) + 3*x(1)*x(1) - x(0)*x(1) + 0.1*sin(x.Norml2());
       default:
          return 1.0;
    }
 }
 
 
+real_t W_exact(const Vector &x)
+{
+   return x(1) + 0.25*cos(2*M_PI*x.Norml2());
+}
+
+
+real_t weight(const Vector &x)
+{
+   return x(0)*x(0) + x(1)*x(1) + 1.0;
+}
+
+
 void visualize(VisItDataCollection &dc, string prefix, int x, int y,
                int visport)
 {
@@ -292,21 +365,32 @@ void visualize(VisItDataCollection &dc, string prefix, int x, int y,
 }
 
 
-real_t compute_mass(FiniteElementSpace *L2, real_t massL2,
-                    VisItDataCollection &dc, string prefix)
+real_t compute_mass(GridFunction &gf, real_t oldmass, string prefix,
+                    CoefficientWithOrder mass_coeff)
 {
+   FiniteElementSpace &fes = *gf.FESpace();
+   Mesh &mesh = *fes.GetMesh();
+
+   // Integration order is a * (element order) + b.
+   const int a = 2;
+   const int b = mesh.GetTypicalElementTransformation()->OrderW() +
+                 mass_coeff.order;
+
    ConstantCoefficient one(1.0);
-   LinearForm lf(L2);
-   lf.AddDomainIntegrator(new DomainLFIntegrator(one));
+   Coefficient &coeff = mass_coeff ? *mass_coeff.coeff : one;
+   DomainLFIntegrator *integ = new DomainLFIntegrator(coeff, a, b);
+
+   LinearForm lf(&fes);
+   lf.AddDomainIntegrator(integ);
    lf.Assemble();
 
-   real_t newmass = lf(*dc.GetField("density"));
+   const real_t newmass = lf(gf);
    cout.precision(18);
    cout << space << " " << prefix << " mass   = " << newmass;
-   if (massL2 >= 0)
+   if (oldmass >= 0)
    {
       cout.precision(4);
-      cout << " ("  << fabs(newmass-massL2)*100/massL2 << "%)";
+      cout << " ("  << fabs(newmass-oldmass)*100/oldmass << "%)";
    }
    cout << endl;
    return newmass;

miniapps/tools/plor-transfer.cpp

@@ -33,6 +33,7 @@
 //
 // Sample runs:  plor-transfer
 //               plor-transfer -h1
+//               plor-transfer -ea -w
 //               plor-transfer -t
 //               plor-transfer -m ../../data/star-q2.mesh -lref 5 -p 4
 //               plor-transfer -m ../../data/star-mixed.mesh -lref 3 -p 2
@@ -59,11 +60,12 @@ string direction;
 
 // Exact functions to project
 real_t RHO_exact(const Vector &x);
+real_t W_exact(const Vector &x);
+real_t weight(const Vector &x);
 
 // Helper functions
 void visualize(VisItDataCollection &, string, int, int, int /* visport */);
-real_t compute_mass(ParFiniteElementSpace *, real_t, VisItDataCollection &,
-                    string);
+real_t compute_mass(ParGridFunction &, real_t, string, CoefficientWithOrder);
 
 int main(int argc, char *argv[])
 {
@@ -80,6 +82,7 @@ int main(int argc, char *argv[])
    bool useH1 = false;
    int visport = 19916;
    bool use_pointwise_transfer = false;
+   bool use_weighted_transfer = false;
    const char *device_config = "cpu";
    bool use_ea       = false;
 
@@ -102,6 +105,9 @@ int main(int argc, char *argv[])
    args.AddOption(&use_pointwise_transfer, "-t", "--use-pointwise-transfer",
                   "-no-t", "--dont-use-pointwise-transfer",
                   "Use pointwise transfer operators instead of L2 projection.");
+   args.AddOption(&use_weighted_transfer, "-w", "--use-weighted-transfer",
+                  "-no-w", "--dont-use-weighted-transfer",
+                  "Use coefficient-weighted L2 projection.");
    args.AddOption(&device_config, "-d", "--device",
                   "Device configuration string, see Device::Configure().");
    args.AddOption(&use_ea, "-ea", "--ea-version", "-no-ea",
@@ -112,6 +118,15 @@ int main(int argc, char *argv[])
    Device device(device_config);
    if (Mpi::Root()) { device.Print(); }
 
+   if (use_weighted_transfer && !use_pointwise_transfer)
+   {
+      if (problem != 5 && Mpi::Root())
+      {
+         cout << "Switching to positive problem = 5 for weighted transfer.\n";
+      }
+      problem = 5;
+   }
+
    // Read the mesh from the given mesh file.
    Mesh serial_mesh(mesh_file, 1, 1);
    ParMesh mesh(MPI_COMM_WORLD, serial_mesh);
@@ -154,6 +169,14 @@ int main(int argc, char *argv[])
    ParFiniteElementSpace fespace(&mesh, fec);
    ParFiniteElementSpace fespace_lor(&mesh_lor, fec_lor);
 
+   FunctionCoefficient weight_fn_coeff(weight);
+   CoefficientWithOrder weight_coeff;
+   if (use_weighted_transfer)
+   {
+      weight_coeff.coeff = &weight_fn_coeff;
+      weight_coeff.order = 2;
+   }
+
    ParGridFunction rho(&fespace);
    ParGridFunction rho_lor(&fespace_lor);
 
@@ -183,7 +206,7 @@ int main(int argc, char *argv[])
    rho.SetTrueVector();
    rho.SetFromTrueVector();
 
-   real_t ho_mass = compute_mass(&fespace, -1.0, HO_dc, "HO       ");
+   real_t ho_mass = compute_mass(rho, -1.0, "HO       ", weight_coeff);
    if (vis) { visualize(HO_dc, "HO", Wx, Wy, visport); Wx += offx; }
 
    GridTransfer *gt;
@@ -193,7 +216,8 @@ int main(int argc, char *argv[])
    }
    else
    {
-      gt = new L2ProjectionGridTransfer(fespace, fespace_lor);
+      gt = new L2ProjectionGridTransfer(fespace, fespace_lor, weight_coeff,
+                                        weight_coeff);
    }
 
    // Configure element assembly for device acceleration
@@ -204,7 +228,7 @@ int main(int argc, char *argv[])
    // HO->LOR restriction
    direction = "HO -> LOR @ LOR";
    R.Mult(rho, rho_lor);
-   compute_mass(&fespace_lor, ho_mass, LOR_dc, "R(HO)    ");
+   compute_mass(rho_lor, ho_mass, "R(HO)    ", weight_coeff);
    if (vis) { visualize(LOR_dc, "R(HO)", Wx, Wy, visport); Wx += offx; }
    auto global_max = [](const Vector& v)
    {
@@ -214,6 +238,47 @@ int main(int argc, char *argv[])
       return max;
    };
 
+   if (use_weighted_transfer && !use_pointwise_transfer)
+   {
+      // Transfer velocity while conserving rho-weighted momentum.
+      GridFunctionCoefficient rho_coeff(&rho);
+      GridFunctionCoefficient rho_lor_coeff(&rho_lor);
+      ProductCoefficient prod_coeff(weight_fn_coeff, rho_coeff);
+      ProductCoefficient prod_lor_coeff(weight_fn_coeff, rho_lor_coeff);
+      CoefficientWithOrder prod_weight(prod_coeff, order + 2);
+      CoefficientWithOrder prod_lor_weight(prod_lor_coeff, lorder + 2);
+
+      ParGridFunction w(&fespace), w_lor(&fespace_lor);
+      FunctionCoefficient W(W_exact);
+      w.ProjectCoefficient(W);
+
+      if (Mpi::Root()) { cout << '\n'; }
+      const real_t ho_momentum = compute_mass(w, -1.0, "rho w HO ", prod_weight);
+
+      L2ProjectionGridTransfer vel_gt(fespace, fespace_lor, prod_weight,
+                                      prod_lor_weight);
+      vel_gt.UseEA(use_ea);
+      vel_gt.ForwardOperator().Mult(w, w_lor);
+      compute_mass(w_lor, ho_momentum, "rho w LOR", prod_lor_weight);
+
+      if (vel_gt.SupportsBackwardsOperator())
+      {
+         ParGridFunction w_prev = w;
+         vel_gt.BackwardOperator().Mult(w_lor, w);
+         compute_mass(w, ho_momentum, "P(rho w) ", prod_weight);
+
+         w_prev -= w;
+         Vector w_prev_true(fespace.GetTrueVSize());
+         w_prev.GetTrueDofs(w_prev_true);
+         const real_t l_inf = global_max(w_prev_true);
+         if (Mpi::Root())
+         {
+            cout.precision(12);
+            cout << "|w - P(R(w))|_∞     = " << l_inf << "\n\n";
+         }
+      }
+   }
+
    if (gt->SupportsBackwardsOperator())
    {
       const Operator &P = gt->BackwardOperator();
@@ -221,7 +286,7 @@ int main(int argc, char *argv[])
       direction = "HO -> LOR @ HO";
       ParGridFunction rho_prev = rho;
       P.Mult(rho_lor, rho);
-      compute_mass(&fespace, ho_mass, HO_dc, "P(R(HO)) ");
+      compute_mass(rho, ho_mass, "P(R(HO)) ", weight_coeff);
       if (vis) { visualize(HO_dc, "P(R(HO))", Wx, Wy, visport); Wx = 0; Wy += offy; }
 
       rho_prev -= rho;
@@ -263,7 +328,7 @@ int main(int argc, char *argv[])
    direction = "LOR -> HO @ LOR";
    rho_lor.ProjectCoefficient(RHO);
    ParGridFunction rho_lor_prev = rho_lor;
-   real_t lor_mass = compute_mass(&fespace_lor, -1.0, LOR_dc, "LOR      ");
+   real_t lor_mass = compute_mass(rho_lor, -1.0, "LOR      ", weight_coeff);
    if (vis) { visualize(LOR_dc, "LOR", Wx, Wy, visport); Wx += offx; }
 
    if (gt->SupportsBackwardsOperator())
@@ -272,14 +337,14 @@ int main(int argc, char *argv[])
       // Prolongate to HO space
       direction = "LOR -> HO @ HO";
       P.Mult(rho_lor, rho);
-      compute_mass(&fespace, lor_mass, HO_dc, "P(LOR)   ");
+      compute_mass(rho, lor_mass, "P(LOR)   ", weight_coeff);
       if (vis) { visualize(HO_dc, "P(LOR)", Wx, Wy, visport); Wx += offx; }
 
       // Restrict back to LOR space. This won't give the original function because
       // the rho_lor doesn't necessarily live in the range of R.
       direction = "LOR -> HO @ LOR";
       R.Mult(rho, rho_lor);
-      compute_mass(&fespace_lor, lor_mass, LOR_dc, "R(P(LOR))");
+      compute_mass(rho_lor, lor_mass, "R(P(LOR))", weight_coeff);
       if (vis) { visualize(LOR_dc, "R(P(LOR))", Wx, Wy, visport); }
 
       rho_lor_prev -= rho_lor;
@@ -334,12 +399,26 @@ real_t RHO_exact(const Vector &x)
          return M_PI/2-atan(5*(2*x.Norml2()-1));
       case 4: // basis function
          return (x.Norml2() < 0.1) ? 1 : 0;
+      case 5: // positive function
+         return 2.0 + 2*x(0)*x(0) + 3*x(1)*x(1) - x(0)*x(1) + 0.1*sin(x.Norml2());
       default:
          return 1.0;
    }
 }
 
 
+real_t W_exact(const Vector &x)
+{
+   return x(1) + 0.25*cos(2*M_PI*x.Norml2());
+}
+
+
+real_t weight(const Vector &x)
+{
+   return x(0)*x(0) + x(1)*x(1) + 1.0;
+}
+
+
 void visualize(VisItDataCollection &dc, string prefix, int x, int y,
                int visport)
 {
@@ -358,23 +437,34 @@ void visualize(VisItDataCollection &dc, string prefix, int x, int y,
 }
 
 
-real_t compute_mass(ParFiniteElementSpace *L2, real_t massL2,
-                    VisItDataCollection &dc, string prefix)
+real_t compute_mass(ParGridFunction &gf, real_t oldmass, string prefix,
+                    CoefficientWithOrder mass_coeff)
 {
+   ParFiniteElementSpace &fes = *gf.ParFESpace();
+   Mesh &mesh = *fes.GetMesh();
+
+   // Integration order is a * (element order) + b.
+   const int a = 2;
+   const int b = mesh.GetTypicalElementTransformation()->OrderW() +
+                 mass_coeff.order;
+
    ConstantCoefficient one(1.0);
-   ParLinearForm lf(L2);
-   lf.AddDomainIntegrator(new DomainLFIntegrator(one));
+   Coefficient &coeff = mass_coeff ? *mass_coeff.coeff : one;
+   DomainLFIntegrator *integ = new DomainLFIntegrator(coeff, a, b);
+
+   ParLinearForm lf(&fes);
+   lf.AddDomainIntegrator(integ);
    lf.Assemble();
 
-   real_t newmass = lf(*dc.GetParField("density"));
+   const real_t newmass = lf(gf);
    if (Mpi::Root())
    {
       cout.precision(18);
       cout << space << " " << prefix << " mass   = " << newmass;
-      if (massL2 >= 0)
+      if (oldmass >= 0)
       {
          cout.precision(4);
-         cout << " ("  << fabs(newmass-massL2)*100/massL2 << "%)";
+         cout << " ("  << fabs(newmass-oldmass)*100/oldmass << "%)";
       }
       cout << endl;
    }

tests/benchmarks/bench_ceed.cpp

@@ -32,7 +32,11 @@
 // Custom benchmark arguments generator
 static void CustomArguments(bm::Benchmark *b) noexcept
 {
-   constexpr int MAX_NDOFS = 16 * 1024 * (mfem_use_gpu ? 1024 : 8);
+#if defined(MFEM_USE_CUDA_OR_HIP_LANG)
+   constexpr int MAX_NDOFS = 16 * 1024 * 1024;
+#else
+   constexpr int MAX_NDOFS = 16 * 1024 * 8;
+#endif
 
    const auto orders = { 7, 6, 5, 4, 3, 2, 1 };
 

tests/unit/CMakeLists.txt

@@ -39,6 +39,7 @@ set(UNIT_TESTS_SRCS
   dfem/test_divergence.cpp
   dfem/test_lvector_interface.cpp
   dfem/test_mass.cpp
+  dfem/test_tuple.cpp
   general/test_array.cpp
   general/test_scan.cpp
   general/test_arrays_by_name.cpp

tests/unit/dfem/test_tuple.cpp

@@ -0,0 +1,274 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "../unit_tests.hpp"
+#include "mfem.hpp"
+#ifndef MFEM_USE_MPI
+#include "../../../fem/dfem/tuple.hpp"
+#endif
+
+using namespace mfem;
+using namespace mfem::future;
+
+namespace tuple_test
+{
+
+// A payload that is not a scalar, mimicking what dFEM kernels actually store.
+using vec3 = tensor<real_t, 3>;
+using tuple3 = tuple<real_t, int, vec3>;
+
+// mfem::future::tuple is no longer an aggregate: it derives from tuple_leaf
+// bases so that it can be defined for an arbitrary number of elements. These
+// checks pin down the properties that the aggregate used to provide for free
+// and that device kernels (which capture tuples by value) depend on.
+static_assert(std::is_trivially_copyable<tuple3>::value,
+              "tuple must be trivially copyable to be captured by value in device kernels");
+static_assert(std::is_trivially_destructible<tuple3>::value,
+              "tuple must be trivially destructible");
+static_assert(std::is_trivially_default_constructible<tuple3>::value,
+              "tuple must be trivially default constructible");
+static_assert(std::is_trivially_copy_assignable<tuple3>::value,
+              "tuple must be trivially copy assignable");
+static_assert(sizeof(tuple3) == sizeof(real_t) + sizeof(int) + sizeof(vec3) +
+              (alignof(real_t) - sizeof(int)),
+              "tuple must not be larger than the sum of its (padded) members");
+
+// Size and element types, both through mfem::future and through the std
+// specializations that drive structured bindings.
+static_assert(tuple_size<tuple3>::value == 3, "");
+static_assert(std::tuple_size<tuple3>::value == 3, "");
+static_assert(std::is_same<tuple_element<0, tuple3>::type, real_t>::value, "");
+static_assert(std::is_same<tuple_element<1, tuple3>::type, int>::value, "");
+static_assert(std::is_same<tuple_element<2, tuple3>::type, vec3>::value, "");
+static_assert(std::is_same<std::tuple_element_t<0, tuple3>, real_t>::value, "");
+static_assert(std::is_same<std::tuple_element_t<2, tuple3>, vec3>::value, "");
+
+// get must preserve the value category and constness of its argument.
+static_assert(std::is_same<decltype(get<1>(std::declval<tuple3&>())),
+              int&>::value, "get on an lvalue must return an lvalue reference");
+static_assert(std::is_same<decltype(get<1>(std::declval<const tuple3&>())),
+              const int&>::value,
+              "get on a const lvalue must return a const lvalue reference");
+static_assert(std::is_same<decltype(get<1>(std::declval<tuple3&&>())),
+              int&&>::value, "get on an rvalue must return an rvalue reference");
+static_assert(std::is_same<decltype(get<1>(std::declval<const tuple3&&>())),
+              const int&&>::value,
+              "get on a const rvalue must return a const rvalue reference");
+
+// += and -= must return a reference, not a copy of the whole tuple.
+using tuple2 = tuple<real_t, vec3>;
+static_assert(std::is_same<decltype(std::declval<tuple2&>() +=
+                                       std::declval<const tuple2&>()), tuple2&>::value,
+              "operator+= must return a reference");
+static_assert(std::is_same<decltype(std::declval<tuple2&>() -=
+                                       std::declval<const tuple2&>()), tuple2&>::value,
+              "operator-= must return a reference");
+
+// The element-wise constructor must stay implicit, so that the
+// copy-list-initialization forms that worked with the aggregate keep working.
+static_assert(std::is_convertible<int, tuple<int>>::value,
+              "tuple's element-wise constructor must not be explicit");
+
+// Constructing from an incompatible type must SFINAE out rather than hard-error,
+// so that the constructor does not poison type traits.
+struct not_a_number { };
+static_assert(!std::is_constructible<tuple<int, int>, int, not_a_number>::value,
+              "");
+static_assert(!std::is_constructible<tuple<int, int>, int>::value,
+              "arity mismatch must not be constructible");
+
+// Usable at compile time.
+constexpr tuple<int, real_t> const_tuple {2, 3.0};
+static_assert(get<0>(const_tuple) == 2, "");
+
+// Copy-list-initialization in a return statement (broken by an explicit ctor).
+tuple<int, real_t> returns_braced_init_list() { return {7, 8.0}; }
+
+} // namespace tuple_test
+
+using namespace tuple_test;
+
+TEST_CASE("dFEM tuple structured bindings", "[dFEM]")
+{
+   tuple3 t {1.0, 2, vec3{{3.0, 4.0, 5.0}}};
+
+   SECTION("binding by reference writes through")
+   {
+      auto &[a, b, c] = t;
+      a = 10.0;
+      b = 20;
+      c(0) = 30.0;
+      REQUIRE(get<0>(t) == 10.0_r);
+      REQUIRE(get<1>(t) == 20);
+      REQUIRE(get<2>(t)(0) == 30.0_r);
+   }
+
+   SECTION("binding by value copies")
+   {
+      auto [a, b, c] = t;
+      a = 10.0;
+      b = 20;
+      c(0) = 30.0;
+      REQUIRE(get<0>(t) == 1.0_r);
+      REQUIRE(get<1>(t) == 2);
+      REQUIRE(get<2>(t)(0) == 3.0_r);
+   }
+
+   SECTION("binding to const")
+   {
+      const auto &[a, b, c] = t;
+      REQUIRE(a == 1.0_r);
+      REQUIRE(b == 2);
+      REQUIRE(c(2) == 5.0_r);
+      static_assert(std::is_same<decltype(a), const real_t>::value, "");
+      static_assert(std::is_same<decltype(c), const vec3>::value, "");
+   }
+
+   SECTION("the bindings alias the tuple storage")
+   {
+      auto &[a, b, c] = t;
+      REQUIRE(&a == &get<0>(t));
+      REQUIRE(&b == &get<1>(t));
+      REQUIRE(&c == &get<2>(t));
+   }
+}
+
+TEST_CASE("dFEM tuple construction", "[dFEM]")
+{
+   SECTION("copy-list-initialization")
+   {
+      tuple<int, real_t> a = {1, 2.0};
+      REQUIRE(get<0>(a) == 1);
+      REQUIRE(get<1>(a) == 2.0_r);
+
+      const auto b = returns_braced_init_list();
+      REQUIRE(get<0>(b) == 7);
+      REQUIRE(get<1>(b) == 8.0_r);
+   }
+
+   SECTION("direct initialization and CTAD")
+   {
+      tuple c {1, 2.0_r, vec3{{1.0, 2.0, 3.0}}};
+      static_assert(std::is_same<decltype(c), tuple<int, real_t, vec3>>::value,
+                    "CTAD must decay the arguments");
+      REQUIRE(get<1>(c) == 2.0_r);
+   }
+
+   SECTION("make_tuple")
+   {
+      const auto d = make_tuple(1, 2.0_r);
+      static_assert(std::is_same<decltype(d), const tuple<int, real_t>>::value, "");
+      REQUIRE(get<0>(d) == 1);
+   }
+
+   SECTION("copy and move construction preserve values")
+   {
+      tuple3 t {1.0, 2, vec3{{3.0, 4.0, 5.0}}};
+      tuple3 copy(t);
+      tuple3 moved(std::move(t));
+      REQUIRE(get<1>(copy) == 2);
+      REQUIRE(get<2>(moved)(1) == 4.0_r);
+   }
+
+   SECTION("value initialization zeroes trivial members")
+   {
+      tuple<int, real_t> z {};
+      REQUIRE(get<0>(z) == 0);
+      REQUIRE(get<1>(z) == 0.0_r);
+   }
+}
+
+TEST_CASE("dFEM tuple arithmetic", "[dFEM]")
+{
+   const tuple2 x {1.0, vec3{{1.0, 2.0, 3.0}}};
+   const tuple2 y {2.0, vec3{{4.0, 5.0, 6.0}}};
+
+   SECTION("element-wise binary operators")
+   {
+      const auto sum = x + y;
+      REQUIRE(get<0>(sum) == 3.0_r);
+      REQUIRE(get<1>(sum)(2) == 9.0_r);
+
+      const auto diff = y - x;
+      REQUIRE(get<0>(diff) == 1.0_r);
+      REQUIRE(get<1>(diff)(0) == 3.0_r);
+   }
+
+   SECTION("compound assignment mutates in place and returns a reference")
+   {
+      tuple2 z = x;
+      auto &ref = (z += y);
+      REQUIRE(&ref == &z);
+      REQUIRE(get<0>(z) == 3.0_r);
+      REQUIRE(get<1>(z)(1) == 7.0_r);
+
+      auto &ref2 = (z -= y);
+      REQUIRE(&ref2 == &z);
+      REQUIRE(get<0>(z) == 1.0_r);
+      REQUIRE(get<1>(z)(1) == 2.0_r);
+   }
+
+   SECTION("scalar operators and unary minus")
+   {
+      const auto scaled = 2.0_r * x;
+      REQUIRE(get<0>(scaled) == 2.0_r);
+      REQUIRE(get<1>(scaled)(2) == 6.0_r);
+
+      const auto halved = x / 2.0_r;
+      REQUIRE(get<0>(halved) == 0.5_r);
+
+      const auto negated = -x;
+      REQUIRE(get<0>(negated) == -1.0_r);
+      REQUIRE(get<1>(negated)(0) == -1.0_r);
+   }
+
+   SECTION("apply")
+   {
+      const auto s = apply([](const real_t &a, const vec3 &b) { return a + b(0); },
+      x);
+      REQUIRE(s == 2.0_r);
+   }
+}
+
+// The tuples are captured by value in device kernels, so exercise a round trip
+// through device memory: construct, mutate through structured bindings and read
+// back on the device.
+TEST_CASE("dFEM tuple on device", "[dFEM][GPU]")
+{
+   Vector res(4);
+   auto d_res = res.Write();
+
+   forall(1, [=] MFEM_HOST_DEVICE (int)
+   {
+      tuple3 t {1.0, 2, vec3{{3.0, 4.0, 5.0}}};
+      auto &[a, b, c] = t;
+      a += static_cast<real_t>(b);
+      c(0) = a;
+
+      tuple2 u {get<0>(t), get<2>(t)};
+      u += tuple2 {1.0, vec3{{1.0, 1.0, 1.0}}};
+
+      d_res[0] = get<0>(u);
+      d_res[1] = get<1>(u)(0);
+      d_res[2] = get<1>(u)(1);
+      d_res[3] = static_cast<real_t>(get<1>(t));
+
+      tuple2 v1{0_r, vec3{0_r, 0_r, 0_r}};
+      tuple2 v2{0_r, vec3{0_r, 0_r, 0_r}};
+      [[maybe_unused]] auto v = v1 + v2;
+   });
+
+   res.HostRead();
+   REQUIRE(std::as_const(res)(0) == 4.0_r);
+   REQUIRE(std::as_const(res)(1) == 4.0_r);
+   REQUIRE(std::as_const(res)(2) == 5.0_r);
+   REQUIRE(std::as_const(res)(3) == 2.0_r);
+}

tests/unit/fem/test_2d_bilininteg.cpp

@@ -3451,4 +3451,81 @@ TEST_CASE("2D Bilinear Scalar Weak Curl Cross Integrators",
    }
 }
 
+TEST_CASE("2D Bilinear Scalar Curl Integrator PartialAssembly",
+          "[MixedScalarCurlIntegrator]"
+          "[BilinearFormIntegrator]"
+          "[NonlinearFormIntegrator]"
+          "[GPU]")
+{
+   int order = 2, n = 1, dim = 2;
+   double tol = 1e-9;
+
+   Mesh mesh = Mesh::MakeCartesian2D(n, n, Element::QUADRILATERAL, 1, 2.0, 3.0);
+
+   VectorFunctionCoefficient   F2_coef(dim, F2);
+   FunctionCoefficient         q2_coef(q2);
+
+   SECTION("Operators on ND")
+   {
+      ND_FECollection    fec_nd(order, dim);
+      FiniteElementSpace fespace_nd(&mesh, &fec_nd);
+
+      GridFunction f_nd(&fespace_nd); f_nd.ProjectCoefficient(F2_coef);
+
+      for (int map_type = (int)FiniteElement::VALUE;
+           map_type <= (int)FiniteElement::INTEGRAL; map_type++)
+      {
+         SECTION("Mapping ND to L2 (" +
+                 MapTypeName((FiniteElement::MapType)map_type) + ")")
+         {
+            L2_FECollection    fec_l2(order - 1, dim,
+                                      BasisType::GaussLegendre,
+                                      (FiniteElement::MapType)map_type);
+            FiniteElementSpace fespace_l2(&mesh, &fec_l2);
+
+            Vector tmp_l2(fespace_l2.GetNDofs());
+            Vector tmp_l2_pa(fespace_l2.GetNDofs());
+
+            SECTION("Without Coefficient")
+            {
+               MixedBilinearForm blf_fa(&fespace_nd, &fespace_l2);
+               blf_fa.AddDomainIntegrator(new MixedScalarCurlIntegrator());
+               blf_fa.Assemble();
+               blf_fa.Finalize();
+
+               blf_fa.Mult(f_nd, tmp_l2);
+
+               MixedBilinearForm blf_pa(&fespace_nd, &fespace_l2);
+               blf_pa.SetAssemblyLevel(mfem::AssemblyLevel::PARTIAL);
+               blf_pa.AddDomainIntegrator(new MixedScalarCurlIntegrator());
+               blf_pa.Assemble();
+
+               blf_pa.Mult(f_nd, tmp_l2_pa);
+               tmp_l2_pa -= tmp_l2;
+               REQUIRE(tmp_l2_pa.Normlinf() < tol);
+            }
+            SECTION("With Scalar Coefficient")
+            {
+               MixedBilinearForm blf_fa(&fespace_nd, &fespace_l2);
+               blf_fa.AddDomainIntegrator(
+                  new MixedScalarCurlIntegrator(q2_coef));
+               blf_fa.Assemble();
+               blf_fa.Finalize();
+
+               blf_fa.Mult(f_nd, tmp_l2);
+
+               MixedBilinearForm blf_pa(&fespace_nd, &fespace_l2);
+               blf_pa.SetAssemblyLevel(mfem::AssemblyLevel::PARTIAL);
+               blf_pa.AddDomainIntegrator(new MixedScalarCurlIntegrator(q2_coef));
+               blf_pa.Assemble();
+
+               blf_pa.Mult(f_nd, tmp_l2_pa);
+               tmp_l2_pa -= tmp_l2;
+               REQUIRE(tmp_l2_pa.Normlinf() < tol);
+            }
+         }
+      }
+   }
+}
+
 } // namespace bilininteg_2d

tests/unit/fem/test_lin_interp.cpp

@@ -1069,4 +1069,238 @@ TEST_CASE("Exact Sequence Properties: d(df)=0",
    }
 }
 
+template <class A, class B>
+static void TestCurl(FiniteElementSpace &dom_fes, FiniteElementSpace &ran_fes,
+                     A coeff, B dcoeff)
+{
+   real_t tol = 1e-10;
+   DiscreteLinearOperator CurlFA(&dom_fes, &ran_fes);
+   CurlFA.AddDomainInterpolator(new CurlInterpolator());
+   CurlFA.Assemble();
+   CurlFA.Finalize();
+
+   SparseMatrix &Curl = CurlFA.SpMat();
+   GridFunction x(&dom_fes), y_fa(&ran_fes), y(&ran_fes);
+   x.ProjectCoefficient(coeff);
+   y.ProjectCoefficient(dcoeff);
+   REQUIRE(x.Size() == Curl.Width());
+   REQUIRE(y_fa.Size() == Curl.Height());
+   Curl.Mult(x, y_fa);
+   y_fa -= y;
+   REQUIRE(y_fa.Normlinf() < tol);
+}
+
+template<class Coeff, class TCoeff>
+static void CompareCurlPA(FiniteElementSpace& dom_fes,
+                          FiniteElementSpace &ran_fes,
+                          Coeff coeff, TCoeff tcoeff)
+{
+   real_t tol = 1e-10;
+   DiscreteLinearOperator CurlFA(&dom_fes, &ran_fes);
+   CurlFA.AddDomainInterpolator(new CurlInterpolator());
+   CurlFA.Assemble();
+   CurlFA.Finalize();
+   DiscreteLinearOperator CurlPA(&dom_fes, &ran_fes);
+   CurlPA.AddDomainInterpolator(new CurlInterpolator());
+   CurlPA.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+   CurlPA.Assemble();
+
+   SparseMatrix &Curl = CurlFA.SpMat();
+   GridFunction x(&dom_fes), y_fa(&ran_fes), y_pa(&ran_fes);
+   x.ProjectCoefficient(coeff);
+   REQUIRE(x.Size() == Curl.Width());
+   REQUIRE(y_fa.Size() == Curl.Height());
+   REQUIRE(x.Size() == CurlPA.Width());
+   REQUIRE(y_pa.Size() == CurlPA.Height());
+   Curl.Mult(x, y_fa);
+   CurlPA.Mult(x, y_pa);
+   y_pa -= y_fa;
+   REQUIRE(y_pa.Normlinf() < tol);
+   // transpose
+   y_fa.ProjectCoefficient(tcoeff);
+   GridFunction x_fa(&dom_fes), x_pa(&dom_fes);
+   Curl.MultTranspose(y_fa, x_fa);
+   CurlPA.MultTranspose(y_fa, x_pa);
+   x_pa -= x_fa;
+   REQUIRE(x_pa.Normlinf() < tol);
+}
+
+TEST_CASE("Partial Assemble Linear Interpolator",
+          "[CurlInterpolator]"
+          "[GPU]")
+{
+   constexpr int maxOrder = 3;
+   auto order = GENERATE_COPY(range(1, maxOrder + 1));
+   CAPTURE(order);
+
+   auto dim = GENERATE(2, 3);
+   CAPTURE(dim);
+
+   int n = 3;
+
+   Mesh mesh;
+
+   switch (dim)
+   {
+      case 2:
+         mesh =
+            Mesh::MakeCartesian2D(n, n, Element::QUADRILATERAL, true, 2.0, 3.0);
+         break;
+      case 3:
+         mesh = Mesh::MakeCartesian3D(n, n, n, Element::HEXAHEDRON, 2.0, 3.0, 5.0);
+         break;
+   }
+
+   // domain spaces
+   H1_FECollection fec_h1(order, dim);
+   FiniteElementSpace fespace_h1(&mesh, &fec_h1);
+   ND_FECollection fec_nd(order, dim);
+   FiniteElementSpace fespace_nd(&mesh, &fec_nd);
+
+   // range spaces
+   RT_FECollection fec_rt(order - 1, dim);
+   FiniteElementSpace fespace_rt(&mesh, &fec_rt);
+   L2_FECollection fec_l2(order - 1, dim, BasisType::GaussLegendre,
+                          FiniteElement::INTEGRAL);
+   FiniteElementSpace fespace_l2(&mesh, &fec_l2);
+
+   switch (dim)
+   {
+      case 2:
+      {
+         FunctionCoefficient coeff([](const Vector &x)
+         { return sin(2 * M_PI * x[1] / 3) - cos(2 * M_PI * x[0] / 2); });
+         VectorFunctionCoefficient vcoeff(2, [](const Vector &x, Vector &y)
+         {
+            y.SetSize(2);
+            y[0] = -cos(2 * M_PI * x[1] / 3);
+            y[1] = sin(2 * M_PI * x[0] / 2);
+         });
+         // out of plane H1 -> in-plane RT
+         SECTION("H1 to RT")
+         {
+            CompareCurlPA(fespace_h1, fespace_rt, coeff, vcoeff);
+         }
+         // in-plane ND -> out of plane L2
+         SECTION("ND to L2")
+         {
+            CompareCurlPA(fespace_nd, fespace_l2, vcoeff, coeff);
+         }
+         break;
+      }
+      case 3:
+      {
+         VectorFunctionCoefficient coeff(3, [](const Vector &x, Vector &y)
+         {
+            y.SetSize(3);
+            y[0] = sin(2 * M_PI * x[2] / 5) - cos(2 * M_PI * x[1] / 3);
+            y[1] = sin(2 * M_PI * x[0] / 2) - cos(2 * M_PI * x[2] / 5);
+            y[2] = sin(2 * M_PI * x[1] / 3) - cos(2 * M_PI * x[0] / 2);
+         });
+         CompareCurlPA(fespace_nd, fespace_rt, coeff, coeff);
+         break;
+      }
+   }
+}
+
+TEST_CASE("Curl Linear Interpolator",
+          "[CurlInterpolator]"
+          "[GPU]")
+{
+   int order = 2;
+
+   auto type = (Element::Type)GENERATE(range((int)Element::TRIANGLE,
+                                             (int)Element::PYRAMID + 1));
+   CAPTURE(type);
+
+   int n = 3;
+
+   Mesh mesh;
+
+   int dim;
+
+   if (type < (int)Element::TETRAHEDRON)
+   {
+      dim = 2;
+      mesh = Mesh::MakeCartesian2D(n, n, (Element::Type)type, 1, 2.0, 3.0);
+   }
+   else
+   {
+      dim = 3;
+      mesh = Mesh::MakeCartesian3D(n, n, n, (Element::Type)type,
+                                   2.0, 3.0, 5.0);
+   }
+
+   // domain spaces
+   H1_FECollection fec_h1(order, dim);
+   FiniteElementSpace fespace_h1(&mesh, &fec_h1);
+   ND_FECollection fec_nd(order, dim);
+   FiniteElementSpace fespace_nd(&mesh, &fec_nd);
+
+   // range spaces
+   RT_FECollection fec_rt(order - 1, dim);
+   FiniteElementSpace fespace_rt(&mesh, &fec_rt);
+   L2_FECollection fec_l2(order - 1, dim, BasisType::GaussLegendre,
+                          FiniteElement::INTEGRAL);
+   FiniteElementSpace fespace_l2(&mesh, &fec_l2);
+
+   switch (dim)
+   {
+      case 2:
+      {
+         // out of plane H1 -> in-plane RT
+         SECTION("H1 to RT")
+         {
+            FunctionCoefficient coeff([](const Vector &x)
+            {
+               return 1 - 2 * x[0] + 3 * x[1];
+            });
+            VectorFunctionCoefficient dcoeff(2, [](const Vector &x, Vector &y)
+            {
+               y.SetSize(2);
+               // d Ez/dy
+               y[0] = 3;
+               // -d Ez/dx
+               y[1] = 2;
+            });
+
+            TestCurl(fespace_h1, fespace_rt, coeff, dcoeff);
+         }
+         // in-plane ND -> out of plane L2
+         SECTION("ND to L2")
+         {
+            VectorFunctionCoefficient coeff(2, [](const Vector &x, Vector &y)
+            {
+               y.SetSize(2);
+               y[0] = 1 - 2 * x[0] + 3 * x[1];
+               y[1] = 2 * (1 - 2 * x[0] + 3 * x[1]);
+            });
+            FunctionCoefficient dcoeff([](const Vector &x)
+            { return 2 * (-2) - 3; });
+            TestCurl(fespace_nd, fespace_l2, coeff, dcoeff);
+         }
+         break;
+      }
+      case 3:
+      {
+         VectorFunctionCoefficient coeff(3, [](const Vector &x, Vector &y)
+         {
+            y.SetSize(3);
+            y[0] = 1 + 2 * x[0] - 3 * x[1] + 4 * x[2];
+            y[1] = 4 + 3 * x[0] - 2 * x[1] + 1 * x[2];
+            y[2] = 2 - 1 * x[0] + 4 * x[1] - 3 * x[2];
+         });
+         VectorFunctionCoefficient dcoeff(3, [](const Vector &x, Vector &y)
+         {
+            y.SetSize(3);
+            y[0] = 4 - 1;
+            y[1] = 4 + 1;
+            y[2] = 3 + 3;
+         });
+         TestCurl(fespace_nd, fespace_rt, coeff, dcoeff);
+         break;
+      }
+   }
+}
+
 } // namespace lin_interp

tests/unit/fem/test_lor_batched.cpp

@@ -214,7 +214,14 @@ TEST_CASE("LOR AMS", "[LOR][BatchedLOR][AMS][Parallel][GPU]")
    ParFiniteElementSpace vert_fespace(edge_fespace.GetParMesh(), &vert_fec);
 
    ParDiscreteLinearOperator grad(&vert_fespace, &edge_fespace);
-   grad.AddDomainInterpolator(new GradientInterpolator);
+   if (space_type == RT)
+   {
+      grad.AddDomainInterpolator(new CurlInterpolator);
+   }
+   else
+   {
+      grad.AddDomainInterpolator(new GradientInterpolator);
+   }
    grad.Assemble();
    grad.Finalize();
    std::unique_ptr<HypreParMatrix> G(grad.ParallelAssemble());

tests/unit/fem/test_pa_coeff.cpp

@@ -750,6 +750,89 @@ TEST_CASE("Hcurl/Hdiv Mixed PA Coefficient",
    }
 }
 
+TEST_CASE("Hcurl/Hdiv MixedVectorGradientPA",
+          "[GPU][PartialAssembly][Coefficient]")
+{
+   constexpr real_t tol = 4e-12;
+   dimension = GENERATE(2, 3);
+   // no coeff, scalar coeff, diagonal matrix coeff, full matrix coeff
+   auto coeffType = GENERATE(0, 1, 2, 3);
+   auto order = GENERATE(1, 2, 3);
+   // RT, ND
+   auto vFEType = GENERATE(0, 1);
+   CAPTURE(dimension, coeffType, order, vFEType);
+
+   const int ne = 3;
+   Mesh mesh = MakeCartesianNonaligned(dimension, ne);
+
+   H1_FECollection scalar_fec(order, dimension);
+   FiniteElementSpace s_fespace(&mesh, &scalar_fec);
+
+   std::unique_ptr<FiniteElementCollection> vector_fec;
+
+   switch (vFEType)
+   {
+      case 0:
+         vector_fec.reset(new RT_FECollection(order - 1, dimension));
+         break;
+      case 1:
+         vector_fec.reset(new ND_FECollection(order, dimension));
+         break;
+   }
+   FiniteElementSpace v_fespace(&mesh, vector_fec.get());
+
+   MixedBilinearForm pa_form(&s_fespace, &v_fespace);
+   pa_form.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+   MixedBilinearForm fa_form(&s_fespace, &v_fespace);
+
+   std::unique_ptr<Coefficient> coeff;
+   std::unique_ptr<DiagonalMatrixCoefficient> dq_coeff;
+   std::unique_ptr<MatrixCoefficient> mq_coeff;
+   switch (coeffType)
+   {
+      case 0:
+         pa_form.AddDomainIntegrator(new MixedVectorGradientIntegrator);
+         fa_form.AddDomainIntegrator(new MixedVectorGradientIntegrator);
+         break;
+      case 1:
+         coeff.reset(new FunctionCoefficient(&coeffFunction));
+         pa_form.AddDomainIntegrator(new MixedVectorGradientIntegrator(*coeff));
+         fa_form.AddDomainIntegrator(new MixedVectorGradientIntegrator(*coeff));
+         break;
+      case 2:
+         dq_coeff.reset(new VectorFunctionCoefficient(dimension, &vectorCoeffFunction));
+         pa_form.AddDomainIntegrator(new MixedVectorGradientIntegrator(*dq_coeff));
+         fa_form.AddDomainIntegrator(new MixedVectorGradientIntegrator(*dq_coeff));
+         break;
+      case 3:
+         mq_coeff.reset(new MatrixFunctionCoefficient(
+                           dimension, &asymmetricMatrixCoeffFunction));
+         pa_form.AddDomainIntegrator(new MixedVectorGradientIntegrator(*mq_coeff));
+         fa_form.AddDomainIntegrator(new MixedVectorGradientIntegrator(*mq_coeff));
+         break;
+   }
+   pa_form.Assemble();
+   fa_form.Assemble();
+
+   GridFunction x_fa(&s_fespace), y_fa(&v_fespace), y_pa(&v_fespace);
+   x_fa.Randomize(1234);
+   REQUIRE(x_fa.Size() == pa_form.Width());
+   REQUIRE(x_fa.Size() == fa_form.Width());
+   REQUIRE(y_fa.Size() == fa_form.Height());
+   REQUIRE(y_pa.Size() == pa_form.Height());
+   pa_form.Mult(x_fa, y_pa);
+   fa_form.Mult(x_fa, y_fa);
+   y_pa -= y_fa;
+   REQUIRE(y_pa.Normlinf() <= tol);
+
+   GridFunction x_pa(&s_fespace);
+   y_fa.Randomize(1234);
+   pa_form.MultTranspose(y_fa, x_pa);
+   fa_form.MultTranspose(y_fa, x_fa);
+   x_pa -= x_fa;
+   REQUIRE(x_pa.Normlinf() <= tol);
+}
+
 TEST_CASE("3D Bilinear VectorFE Integrators PartialAssembly",
           "[BilinearFormIntegrator]"
           "[PartialAssembly]"
@@ -1059,4 +1142,111 @@ TEST_CASE("3D Bilinear VectorFE Integrators PartialAssembly",
    }
 }
 
+TEST_CASE("3D Bilinear Weak Curl Integrators Partial Assembly",
+          "[MixedVectorWeakCurlIntegrator]"
+          "[BilinearFormIntegrator]"
+          "[PartialAssembly]"
+          "[GPU]")
+{
+   auto order = GENERATE(1, 2);
+   CAPTURE(order);
+   int dim = 3;
+
+   FunctionCoefficient q3_coeff(coeffFunction);
+   VectorFunctionCoefficient F3_coeff(dim, vectorCoeffFunction);
+
+   auto mesh_fname =
+      GENERATE("../../data/fichera-amr.mesh", "../../data/ball-nurbs.mesh");
+   CAPTURE(mesh_fname);
+   Mesh mesh(mesh_fname);
+   REQUIRE(mesh.Dimension() == dim);
+   REQUIRE(mesh.SpaceDimension() == dim);
+
+   // convert nurbs into piecewise-quadratic curved mesh
+   if (mesh.NURBSext)
+   {
+      mesh.UniformRefinement();
+      mesh.SetCurvature(2);
+   }
+
+   SECTION("RT to ND No Coeff")
+   {
+      ND_FECollection fec_nd(order, dim);
+      FiniteElementSpace fespace_nd(&mesh, &fec_nd);
+      RT_FECollection fec_rt(order - 1, dim);
+      FiniteElementSpace fespace_rt(&mesh, &fec_rt);
+
+      MixedBilinearForm bfa(&fespace_rt, &fespace_nd);
+      bfa.AddDomainIntegrator(new MixedVectorWeakCurlIntegrator);
+      bfa.Assemble();
+      bfa.Finalize();
+
+      MixedBilinearForm bpa(&fespace_rt, &fespace_nd);
+      bpa.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+      bpa.AddDomainIntegrator(new MixedVectorWeakCurlIntegrator);
+      bpa.Assemble();
+
+      GridFunction x(&fespace_rt), y_fa(&fespace_nd), y_pa(&fespace_nd);
+      x.Randomize(1234);
+      REQUIRE(bfa.Height() == y_fa.Size());
+      REQUIRE(bfa.Width() == x.Size());
+      REQUIRE(bpa.Height() == y_fa.Size());
+      REQUIRE(bpa.Width() == x.Size());
+      bfa.Mult(x, y_fa);
+      bpa.Mult(x, y_pa);
+      y_pa -= y_fa;
+      REQUIRE( y_pa.Normlinf() == MFEM_Approx(0_r) );
+   }
+
+   SECTION("RT to ND Scalar Coeff")
+   {
+      ND_FECollection fec_nd(order, dim);
+      FiniteElementSpace fespace_nd(&mesh, &fec_nd);
+      RT_FECollection fec_rt(order - 1, dim);
+      FiniteElementSpace fespace_rt(&mesh, &fec_rt);
+
+      MixedBilinearForm bfa(&fespace_rt, &fespace_nd);
+      bfa.AddDomainIntegrator(new MixedVectorWeakCurlIntegrator(q3_coeff));
+      bfa.Assemble();
+      bfa.Finalize();
+
+      MixedBilinearForm bpa(&fespace_rt, &fespace_nd);
+      bpa.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+      bpa.AddDomainIntegrator(new MixedVectorWeakCurlIntegrator(q3_coeff));
+      bpa.Assemble();
+
+      GridFunction x(&fespace_rt), y_fa(&fespace_nd), y_pa(&fespace_nd);
+      x.Randomize(1234);
+      bfa.Mult(x, y_fa);
+      bpa.Mult(x, y_pa);
+      y_pa -= y_fa;
+      REQUIRE( y_pa.Normlinf() == MFEM_Approx(0_r) );
+   }
+
+   SECTION("RT to ND Diagonal Matrix Coeff")
+   {
+      ND_FECollection fec_nd(order, dim);
+      FiniteElementSpace fespace_nd(&mesh, &fec_nd);
+      RT_FECollection fec_rt(order - 1, dim);
+      FiniteElementSpace fespace_rt(&mesh, &fec_rt);
+
+      MixedBilinearForm bfa(&fespace_rt, &fespace_nd);
+      bfa.AddDomainIntegrator(new MixedVectorWeakCurlIntegrator(F3_coeff));
+      bfa.Assemble();
+      bfa.Finalize();
+
+      MixedBilinearForm bpa(&fespace_rt, &fespace_nd);
+      bpa.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+      bpa.AddDomainIntegrator(new MixedVectorWeakCurlIntegrator(F3_coeff));
+      bpa.Assemble();
+
+      GridFunction x(&fespace_rt), y_fa(&fespace_nd), y_pa(&fespace_nd);
+      x.Randomize(1234);
+      bfa.Mult(x, y_fa);
+      bpa.Mult(x, y_pa);
+      y_pa -= y_fa;
+      REQUIRE( y_pa.Normlinf() == MFEM_Approx(0_r) );
+   }
+}
+
 } // namespace pa_coeff

tests/unit/linalg/test_complex_operator.cpp

@@ -164,12 +164,29 @@ TEST_CASE("ComplexHypreParMatrix GetSystemMatrix",
    a.AddDomainIntegrator(new VectorFEMassIntegrator(one),
                          new VectorFEMassIntegrator(one));
    a.Assemble();
+
+   // 2. Test ParSesquilinearForm::FormSystemMatrix directly and verify that
+   //    essential entries on the imaginary diagonal are zero.
    OperatorPtr Ah;
+   a.FormSystemMatrix(ess_tdof_list, Ah);
+   ComplexHypreParMatrix *A_complex = Ah.Is<ComplexHypreParMatrix>();
+   REQUIRE(A_complex != nullptr);
+   Vector diag;
+   A_complex->imag().GetDiag(diag);
+   const Array<int> &ess_tdofs = ess_tdof_list;
+   const Vector &diag_h = diag;
+   ess_tdofs.HostRead();
+   diag_h.HostRead();
+   for (const int tdof : ess_tdofs)
+   {
+      REQUIRE(diag_h[tdof] == 0.0);
+   }
+
+   // 3. Test the call to ComplexHypreParMatrix::GetSystemMatrix and destroying
+   //    the returned matrix.
    Vector B, X;
    a.FormLinearSystem(ess_tdof_list, x, b, Ah, X, B);
 
-   // 2. Test the call to ComplexHypreParMatrix::GetSystemMatrix and destroying
-   //    the returned matrix.
    HypreParMatrix *A = Ah.As<ComplexHypreParMatrix>()->GetSystemMatrix();
    delete A;
 }

tests/unit/mesh/test_bb_grid_map.cpp

@@ -152,7 +152,7 @@ TEST_CASE("GlobalBBoxTensorGridMap Parallel",
    std::map<int, std::vector<int>> pt_to_procs;
    map.MapPointsToProcs(centers, 1, pt_to_procs);
 
-   REQUIRE(pt_to_procs.size() == nel + 1);
+   REQUIRE(pt_to_procs.size() == (unsigned)nel + 1);
    for (int i = 0; i < nel; i++)
    {
       std::vector<int> procs = pt_to_procs[i];

tests/unit/mesh/test_ncmesh.cpp

@@ -304,6 +304,66 @@ TEST_CASE("pNCMesh PA diagonal",  "[Parallel], [NCMesh]")
    }
 } // test case
 
+TEST_CASE("ParNCMesh Rebalance preserves element attributes",
+          "[Parallel], [NCMesh]")
+{
+   const int rank = Mpi::WorldRank();
+   const int nranks = Mpi::WorldSize();
+   if (nranks < 2) { return; }
+
+   auto mesh_fname = GENERATE("../../data/star.mesh",
+                              "../../data/fichera.mesh");
+   CAPTURE(mesh_fname);
+
+   auto CheckRebalance = [rank, nranks, mesh_fname](bool refine,
+                                                    bool custom_partition)
+   {
+      Mesh mesh(mesh_fname);
+      mesh.EnsureNCMesh();
+      ParMesh pmesh(MPI_COMM_WORLD, mesh);
+
+      const int attribute = 1234 + (custom_partition ? rank : 0);
+      for (int i = 0; i < pmesh.GetNE(); i++)
+      {
+         pmesh.SetAttribute(i, attribute);
+      }
+      pmesh.SetAttributes();
+
+      if (refine)
+      {
+         Array<int> refinements;
+         if (pmesh.GetNE() && (custom_partition || rank == 0))
+         {
+            refinements.Append(0);
+         }
+         pmesh.GeneralRefinement(refinements);
+      }
+
+      int expected_attribute = attribute;
+      if (custom_partition)
+      {
+         // Move every element to the next rank, as in GitHub issue #4009.
+         Array<int> partition(pmesh.GetNE());
+         partition = (rank + 1) % nranks;
+         pmesh.Rebalance(partition);
+         expected_attribute = 1234 + (rank + nranks - 1) % nranks;
+      }
+      else
+      {
+         pmesh.Rebalance();
+      }
+
+      for (int i = 0; i < pmesh.GetNE(); i++)
+      {
+         CHECK(pmesh.GetAttribute(i) == expected_attribute);
+      }
+   };
+
+   SECTION("Custom partition, unrefined") { CheckRebalance(false, true); }
+   SECTION("Custom partition, refined") { CheckRebalance(true, true); }
+   SECTION("Default partition, refined") { CheckRebalance(true, false); }
+}
+
 TEST_CASE("EdgeFaceConstraint", "[Parallel], [NCMesh]")
 {
    auto exact_soln = [](const Vector& x)