Fix some CuDSS lifetime issues

refactor tuple for generic size

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,20 @@ 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 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 +123,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 +140,12 @@ Miscellaneous
   using the new method ApplyDofSigns() in class ParFiniteElementSpace: the
   method will return immediately if no sign flips are needed.
 
+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.

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=
@@ -382,7 +395,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),\
@@ -665,3 +678,15 @@ 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)

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

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/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/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/hybridization_ext.cpp

@@ -368,6 +368,8 @@ void HybridizationExtension::ConstructH()
 
    CAhatInvCt = 0.0;
 
+   // Fill the face-to-face adjacency array. Two faces are adjacent if they are
+   // incident to a common element.
    mfem::forall(nf, [=] MFEM_HOST_DEVICE (int fi)
    {
       const int begin_f = d_face_face_offsets[fi];
@@ -403,6 +405,12 @@ void HybridizationExtension::ConstructH()
             }
          }
       }
+      // Fill unused entries with -1 to indicate invalid
+      const int end_f = d_face_face_offsets[fi + 1];
+      for (int i = begin_f + idx; i < end_f; ++i)
+      {
+         d_face_to_face[i] = -1;
+      }
    });
 
    mfem::forall(nf, [=] MFEM_HOST_DEVICE (int fi)
@@ -412,6 +420,7 @@ void HybridizationExtension::ConstructH()
       for (int idx_j = begin; idx_j < end; ++idx_j)
       {
          const int fj = d_face_to_face[idx_j];
+         if (fj < 0) { break; }
          for (int ei = 0; ei < 2; ++ei)
          {
             const int e = d_face_to_el(0, ei, fi);

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_mixedvecgrad_pa.cpp

@@ -672,6 +672,8 @@ void MixedVectorGradientIntegrator::AssemblePA(const FiniteElementSpace
    const NodalTensorFiniteElement *trial_el =
       dynamic_cast<const NodalTensorFiniteElement*>(trial_fel);
    MFEM_VERIFY(trial_el != NULL, "Only NodalTensorFiniteElement is supported!");
+   MFEM_VERIFY(trial_el->GetMapType() == FiniteElement::VALUE,
+               "Only value map type is supported!");
 
    const VectorTensorFiniteElement *test_el =
       dynamic_cast<const VectorTensorFiniteElement*>(test_fel);

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

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/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

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/cudss.cpp

@@ -87,6 +87,9 @@ CuDSSSolver::CuDSSSolver(MPI_Comm comm_) : mpi_comm(comm_)
 
 CuDSSSolver::~CuDSSSolver()
 {
+   // Sync the stream to make sure any pending asynchronous operations have
+   // completed.
+   MFEM_STREAM_SYNC;
    // Destroy the system Matrix, RHS vector and solution vector
    if (Ac)
    {
@@ -99,7 +102,6 @@ CuDSSSolver::~CuDSSSolver()
    MFEM_CUDSS_CHECK(cudssDataDestroy(handle, solverData));
    MFEM_CUDSS_CHECK(cudssConfigDestroy(solverConfig));
 
-
    MFEM_CUDSS_CHECK(cudssDestroy(handle));
    handle = nullptr;
 
@@ -125,6 +127,9 @@ void CuDSSSolver::InitCuDSS()
    // Create the cuDSS handle
    MFEM_CUDSS_CHECK(cudssCreate(&handle));
 
+   // Set CuDSS to use MFEM's default stream of 0.
+   MFEM_CUDSS_CHECK(cudssSetStream(handle, 0));
+
 #ifdef MFEM_USE_OPENMP
    // NOTE: Set the threading layer library name to NULL so that cuDSS picks
    // it from the environment variable "CUDSS_THREADING_LIB"
@@ -251,27 +256,42 @@ void CuDSSSolver::SetMatrixCuDSS(int *csr_offsets, int *csr_columns,
       Ac = std::make_unique<cudssMatrix_t>();
       // Create empty RHS and solution vectors
       SetNumRHS(1);
-      // Allocate device memory for csr values
-      CuMemAlloc(&csr_values_d, nnz * sizeof(real_t));
    }
 
+   if (cuDSSObjectInitialized && !reorder_reuse)
+   {
+      MFEM_STREAM_SYNC;
+      MFEM_CUDSS_CHECK(cudssMatrixDestroy(*Ac));
+   }
+
+   // Allocate device memory for csr values. Unless reuse is specified, the
+   // nnz may be different, so we will free and reallocate.
+   if (csr_values_d == NULL || !reorder_reuse)
+   {
+      if (csr_values_d != NULL) { CuMemFree(csr_values_d); }
+      CuMemAlloc(&csr_values_d, nnz * sizeof(real_t));
+   }
    CuMemcpyDtoD(csr_values_d, csr_values, nnz * sizeof(real_t));
 
+   // We copy and store the I and J arrays, since the CuDSS matrix object
+   // technically needs these to be valid, so we protect against the caller
+   // destroying the original matrix.
+   if (!cuDSSObjectInitialized || !reorder_reuse)
+   {
+      if (csr_offsets_d != NULL) { CuMemFree(csr_offsets_d); }
+      CuMemAlloc(&csr_offsets_d, (n_loc + 1) * sizeof(int));
+      if (csr_columns_d != NULL) { CuMemFree(csr_columns_d); }
+      CuMemAlloc(&csr_columns_d, nnz * sizeof(int));
+      CuMemcpyDtoD(csr_offsets_d, csr_offsets, (n_loc + 1) * sizeof(int));
+      CuMemcpyDtoD(csr_columns_d, csr_columns, nnz * sizeof(int));
+   }
+
    // New cuDSS CSR matrix object and analysis or reuse the one from a previous
    // matrix
    if (!cuDSSObjectInitialized || !reorder_reuse)
    {
       if (reorder_reuse)  // !cuDSSObjectInitialized && reorder_reuse
       {
-         // NOTE: For CuDSS solver to reuse the reordering (skipping analysis
-         // phase), it needs to access the I and J arrays of the **initial**
-         // matrix. Therefore, we need to copy and keep I and J in device memory.
-         CuMemAlloc(&csr_offsets_d, (n_loc + 1) * sizeof(int));
-         CuMemAlloc(&csr_columns_d, nnz * sizeof(int));
-
-         CuMemcpyDtoD(csr_offsets_d, csr_offsets, (n_loc + 1) * sizeof(int));
-         CuMemcpyDtoD(csr_columns_d, csr_columns, nnz * sizeof(int));
-
 #if CUDSS_VERSION >= 800
          MFEM_CUDSS_CHECK(
             cudssMatrixCreateCsr(
@@ -288,21 +308,17 @@ void CuDSSSolver::SetMatrixCuDSS(int *csr_offsets, int *csr_columns,
       }
       else    // !reorder_reuse
       {
-         if (cuDSSObjectInitialized)
-         {
-            MFEM_CUDSS_CHECK(cudssMatrixDestroy(*Ac));
-         }
 #if CUDSS_VERSION >= 800
          MFEM_CUDSS_CHECK(
             cudssMatrixCreateCsr(
-               Ac.get(), n_global, n_global, nnz, csr_offsets, NULL,
-               csr_columns, csr_values_d, CUDSS_INT_T, CUDSS_INT_T, CUDSS_REAL_T,
+               Ac.get(), n_global, n_global, nnz, csr_offsets_d, NULL,
+               csr_columns_d, csr_values_d, CUDSS_INT_T, CUDSS_INT_T, CUDSS_REAL_T,
                mat_type, mview, CUDSS_BASE_ZERO));
 #else
          MFEM_CUDSS_CHECK(
             cudssMatrixCreateCsr(
-               Ac.get(), n_global, n_global, nnz, csr_offsets, NULL,
-               csr_columns, csr_values_d, CUDSS_INT_T, CUDSS_REAL_T,
+               Ac.get(), n_global, n_global, nnz, csr_offsets_d, NULL,
+               csr_columns_d, csr_values_d, CUDSS_INT_T, CUDSS_REAL_T,
                mat_type, mview, CUDSS_BASE_ZERO));
 #endif
       }
@@ -326,6 +342,9 @@ void CuDSSSolver::SetMatrixCuDSS(int *csr_offsets, int *csr_columns,
    // Factorization
    MFEM_CUDSS_CHECK(cudssExecute(handle, CUDSS_PHASE_FACTORIZATION, solverConfig,
                                  solverData, *Ac, yc, xc));
+
+   // In serial, the factorization can execute asynchronously.
+   MFEM_STREAM_SYNC;
 }
 
 void CuDSSSolver::SetOperator(const Operator &op)
@@ -360,6 +379,7 @@ void CuDSSSolver::SetNumRHS(int nrhs_) const
       if (nrhs > 0)
       {
          // Destroy the previous RHS vector and solution vector
+         MFEM_STREAM_SYNC;
          MFEM_CUDSS_CHECK(cudssMatrixDestroy(xc));
          MFEM_CUDSS_CHECK(cudssMatrixDestroy(yc));
       }

linalg/cudss.hpp

@@ -157,8 +157,7 @@ private:
    mutable int nrhs = 0;  // the number of the RHSs
    int nnz = 0;        // the number of non zeros
 
-   // copy and keep the I and J arrays in device memory when skipping analysis
-   // phase
+   // copy and keep the I and J arrays in device memory
    void *csr_offsets_d = NULL; // copy and keep I in device
    void *csr_columns_d = NULL; // copy and keep J in device
    void *csr_values_d = NULL;  // copy and keep csr data in device

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. */

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;
    };
 

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();

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

@@ -1059,4 +1059,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)