Merge branch 'master' into codespell-action

Bump RAJA required C++ version in CMake

.github/workflows/repo-check.yml

@@ -176,3 +176,14 @@ jobs:
         run: |
           echo "::warning::branch-history check failed, but the" \
             "'branch-history-override' label is set."
+
+  codespell:
+    runs-on: ubuntu-latest
+    steps:
+    - name: Codespell with annotations
+      uses: codespell-project/actions-codespell@master
+      with:
+        check_filenames: true
+        check_hidden: true
+        ignore_words_list: allright,ba,equil,esy,fo,hda,lod,nd,ned,numer,ot,pres,ro,seh,shat,solfes,strat,tbe,te,warmup
+

.gitignore

@@ -260,6 +260,7 @@ miniapps/meshing/polar-nc
 miniapps/meshing/mesh-quality
 miniapps/meshing/hpref
 miniapps/meshing/phpref
+miniapps/meshing/pref321
 miniapps/meshing/mobius-strip.mesh
 miniapps/meshing/klein-bottle.mesh
 miniapps/meshing/toroid-*.mesh

CHANGELOG

@@ -46,8 +46,19 @@ Discretization improvements
 
 - Extend FindPointsGSLIB to support surface meshes.
 
+- Added support for complex-valued mixed bilinear forms via the new classes
+  MixedSesquilinearForm and ParMixedSesquilinearForm, mirroring the existing
+  SesquilinearForm classes. Rectangular complex operators are now also
+  handled correctly by ComplexSparseMatrix::GetSystemMatrix and
+  ComplexHypreParMatrix::GetSystemMatrix, which previously assumed equal
+  trial and test spaces.
+
 Meshing improvements
 --------------------
+- Added support for nonuniform anisotropic mesh refinement on parallel quad/hex
+  meshes with arbitrary spacing in each direction. This enables in particular
+  3:1 refinement in parallel, as demonstrated in the new meshing miniapp pref321.
+
 - Added option to guarantee mesh validity during TMOP-based r-adaptivity, using
   bounds on the determinant of the mesh transformation Jacobian.
 
@@ -87,6 +98,9 @@ GPU computing
 
 - Allow specifying GPU kernel launch bounds for native and RAJA GPU backends.
 
+- Changed VectorFEMassIntegrator to use kernel specialization dispatch for
+  partial assembly.
+
 - Added support for FiniteElement::MapType::INTEGRAL spaces to
   QuadratureInterpolator.
 

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_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.

doc/CodeDocumentation.dox

@@ -201,6 +201,7 @@ namespace mfem {
  * - <a class="el" href="nurbs__naca__cmesh_8cpp_source.html">NURBS NACA Mesher</a>: generate NURBS based mesh around a NACA foil
  * - <a class="el" href="nurbs__printfunc_8cpp_source.html">NURBS Printer</a>: print the NURBS-basis
  * - <a class="el" href="nurbs__mesh_info_8cpp_source.html">NURBS Mesh info</a>: print the info of a NURBS mesh
+ * - <a class="el" href="nurbs__surface_8cpp_source.html">NURBS Surface</a>: interpolate a 3D Surface in a NURBS Patch
  *
  * <H3>Miniapps</H3>
  * - <a class="el" href="volta_8cpp_source.html">Volta</a>: simple electrostatics simulation code
@@ -245,6 +246,9 @@ namespace mfem {
  * - <a class="el" href="pdiffusion_8cpp_source.html">DPG Diffusion example</a>: DPG formulation for the diffusion problem
  * - <a class="el" href="pmaxwell_8cpp_source.html">DPG Maxwell example</a>: DPG formulation for the indefinite Maxwell problem
  * - <a class="el" href="lor__elast_8cpp_source.html">LOR Elasticity</a>: solve linear elasticity with LOR preconditioning on GPUs
+ * - <a class="el" href="reflector_8cpp_source.html">Reflector Miniapp</a>: reflect a mesh about a plane
+ * - <a class="el" href="ref321_8cpp_source.html">3:1 Refinement Miniapp</a>: perform 3:1 anisotropic mesh refinements
+ * - <a class="el" href="pref321_8cpp_source.html">3:1 Refinement Miniapp</a>: parallel 3:1 anisotropic mesh refinements
  *
  * See also the <a class="el" href="https://mfem.org/examples/">examples documentation</a> online.
  */

examples/ex27.cpp

@@ -177,7 +177,7 @@ int main(int argc, char *argv[])
    Array<int> ess_tdof_list(0);
    if (h1 && mesh->bdr_attributes.Size())
    {
-      // For a continuous basis the linear system must be modified to enforce an
+      // For a continuous basis the linear system must be modifed to enforce an
       // essential (Dirichlet) boundary condition. In the DG case this is not
       // necessary as the boundary condition will only be enforced weakly.
       fespace.GetEssentialTrueDofs(dbc_bdr, ess_tdof_list);

fem/bilininteg.hpp

@@ -3003,11 +3003,10 @@ public:
     vector (diagonal matrix), or matrix), trial function $u$ is in $H(curl$ or
     $H(div)$, and test function $v$ is in $H(curl$, $H(div)$, or $v=(v_1,\dots,v_n)$, where
     $v_i$ are in $H^1$. */
-class VectorFEMassIntegrator: public BilinearFormIntegrator
+class VectorFEMassIntegrator : public BilinearFormIntegrator
 {
 private:
-   void Init(Coefficient *q, DiagonalMatrixCoefficient *dq, MatrixCoefficient *mq)
-   { Q = q; DQ = dq; MQ = mq; }
+   void Init(Coefficient *q, DiagonalMatrixCoefficient *dq, MatrixCoefficient *mq);
 
 #ifndef MFEM_THREAD_SAFE
    Vector shape;
@@ -3030,7 +3029,8 @@ protected:
    const DofToQuad *mapsOtest;     ///< Not owned. DOF-to-quad map, open.
    const DofToQuad *mapsCtest;     ///< Not owned. DOF-to-quad map, closed.
    const GeometricFactors *geom;   ///< Not owned
-   int dim, ne, nq, dofs1D, dofs1Dtest, quad1D, trial_fetype, test_fetype;
+   int dim, ne, nq, dofs1D, dofs1Dtest, quad1D;
+   FiniteElement::DerivType trial_fetype, test_fetype;
    bool symmetric = true; ///< False if using a nonsymmetric matrix coefficient
 
 public:
@@ -3061,6 +3061,29 @@ public:
                    const bool add) override;
 
    const Coefficient *GetCoefficient() const { return Q; }
+
+   using ApplyKernelType =
+      void (*)(const int NE, bool symmetric, const bool scalar_coeff,
+               const Array<real_t> &trialBO, const Array<real_t> &trialBC,
+               const Array<real_t> &testBOt, const Array<real_t> &testBCt,
+               const Vector &pa_data, const Vector &x, Vector &y,
+               const int triald1d, const int testd1d, const int q1d);
+
+   /// parameters: trial_fetype, test_fetype, ndims, trial_d1d, test_d1d, q1d
+   MFEM_REGISTER_KERNELS(ApplyPAKernels, ApplyKernelType,
+                         (FiniteElement::DerivType, FiniteElement::DerivType,
+                          int, int, int, int));
+
+   struct Kernels { Kernels(); };
+
+   template <FiniteElement::DerivType TrialType,
+             FiniteElement::DerivType TestType, int DIM, int TRIAL_D1D,
+             int TEST_D1D, int Q1D>
+   static void AddSpecialization()
+   {
+      ApplyPAKernels::Specialization<TrialType, TestType, DIM, TRIAL_D1D,
+                     TEST_D1D, Q1D>::Add();
+   }
 };
 
 /** Integrator for $(Q \nabla \cdot u, v)$ where $u=(u_1,\cdots,u_n)$ and all $u_i$ are in the same

fem/complex_fem.cpp

@@ -718,8 +718,8 @@ SesquilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list,
 
    // A = A_r + i A_i
    A.Clear();
-   if ( A_r.Type() == Operator::MFEM_SPARSEMAT ||
-        A_i.Type() == Operator::MFEM_SPARSEMAT )
+   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>(),
@@ -779,8 +779,8 @@ SesquilinearForm::FormSystemMatrix(const Array<int> &ess_tdof_list,
 
    // A = A_r + i A_i
    A.Clear();
-   if ( A_r.Type() == Operator::MFEM_SPARSEMAT ||
-        A_i.Type() == Operator::MFEM_SPARSEMAT )
+   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>(),
@@ -843,6 +843,426 @@ SesquilinearForm::Update(FiniteElementSpace *nfes)
    if ( blfi ) { blfi->Update(nfes); }
 }
 
+bool
+MixedSesquilinearForm::RealInteg()
+{
+   int nint = mblfr->GetDBFI()->Size() + mblfr->GetBBFI()->Size() +
+              mblfr->GetFBFI()->Size() + mblfr->GetBFBFI()->Size() +
+              mblfr->GetTFBFI()->Size() + mblfr->GetBTFBFI()->Size();
+   return (nint != 0);
+}
+
+bool
+MixedSesquilinearForm::ImagInteg()
+{
+   int nint = mblfi->GetDBFI()->Size() + mblfi->GetBBFI()->Size() +
+              mblfi->GetFBFI()->Size() + mblfi->GetBFBFI()->Size() +
+              mblfi->GetTFBFI()->Size() + mblfi->GetBTFBFI()->Size();
+   return (nint != 0);
+}
+
+MixedSesquilinearForm::MixedSesquilinearForm(FiniteElementSpace * trial_fes,
+                                             FiniteElementSpace * test_fes,
+                                             ComplexOperator::Convention convention)
+   : conv(convention),
+     mblfr(new mfem::MixedBilinearForm(trial_fes, test_fes)),
+     mblfi(new mfem::MixedBilinearForm(trial_fes, test_fes))
+{
+}
+
+MixedSesquilinearForm::MixedSesquilinearForm(FiniteElementSpace * trial_fes,
+                                             FiniteElementSpace * test_fes,
+                                             MixedBilinearForm * bfr,
+                                             MixedBilinearForm * bfi,
+                                             ComplexOperator::Convention convention)
+   : conv(convention),
+     mblfr(new MixedBilinearForm(trial_fes, test_fes, bfr)),
+     mblfi(new MixedBilinearForm(trial_fes, test_fes, bfi))
+{
+}
+
+MixedSesquilinearForm::~MixedSesquilinearForm()
+{
+   delete mblfr;
+   delete mblfi;
+}
+
+void
+MixedSesquilinearForm::AddDomainIntegrator(BilinearFormIntegrator * bfi_real,
+                                           BilinearFormIntegrator * bfi_imag)
+{
+   if (bfi_real)
+   {
+      mblfr->AddDomainIntegrator(bfi_real);
+   }
+   if (bfi_imag)
+   {
+      mblfi->AddDomainIntegrator(bfi_imag);
+   }
+}
+
+void
+MixedSesquilinearForm::AddDomainIntegrator(BilinearFormIntegrator * bfi_real,
+                                           BilinearFormIntegrator * bfi_imag,
+                                           Array<int> & elem_marker)
+{
+   if (bfi_real)
+   {
+      mblfr->AddDomainIntegrator(bfi_real, elem_marker);
+   }
+   if (bfi_imag)
+   {
+      mblfi->AddDomainIntegrator(bfi_imag, elem_marker);
+   }
+}
+
+void
+MixedSesquilinearForm::AddBoundaryIntegrator(BilinearFormIntegrator * bfi_real,
+                                             BilinearFormIntegrator * bfi_imag)
+{
+   if (bfi_real)
+   {
+      mblfr->AddBoundaryIntegrator(bfi_real);
+   }
+   if (bfi_imag)
+   {
+      mblfi->AddBoundaryIntegrator(bfi_imag);
+   }
+}
+
+void
+MixedSesquilinearForm::AddBoundaryIntegrator(BilinearFormIntegrator * bfi_real,
+                                             BilinearFormIntegrator * bfi_imag,
+                                             Array<int> & bdr_marker)
+{
+   if (bfi_real)
+   {
+      mblfr->AddBoundaryIntegrator(bfi_real, bdr_marker);
+   }
+   if (bfi_imag)
+   {
+      mblfi->AddBoundaryIntegrator(bfi_imag, bdr_marker);
+   }
+}
+
+void
+MixedSesquilinearForm::AddInteriorFaceIntegrator(BilinearFormIntegrator *
+                                                 bfi_real,
+                                                 BilinearFormIntegrator * bfi_imag)
+{
+   if (bfi_real)
+   {
+      mblfr->AddInteriorFaceIntegrator(bfi_real);
+   }
+   if (bfi_imag)
+   {
+      mblfi->AddInteriorFaceIntegrator(bfi_imag);
+   }
+}
+
+void
+MixedSesquilinearForm::AddBdrFaceIntegrator(BilinearFormIntegrator * bfi_real,
+                                            BilinearFormIntegrator * bfi_imag)
+{
+   if (bfi_real)
+   {
+      mblfr->AddBdrFaceIntegrator(bfi_real);
+   }
+   if (bfi_imag)
+   {
+      mblfi->AddBdrFaceIntegrator(bfi_imag);
+   }
+}
+
+void
+MixedSesquilinearForm::AddBdrFaceIntegrator(BilinearFormIntegrator * bfi_real,
+                                            BilinearFormIntegrator * bfi_imag,
+                                            Array<int> & bdr_marker)
+{
+   if (bfi_real)
+   {
+      mblfr->AddBdrFaceIntegrator(bfi_real, bdr_marker);
+   }
+   if (bfi_imag)
+   {
+      mblfi->AddBdrFaceIntegrator(bfi_imag, bdr_marker);
+   }
+}
+
+void MixedSesquilinearForm::AddTraceFaceIntegrator(BilinearFormIntegrator *
+                                                   bfi_real,
+                                                   BilinearFormIntegrator * bfi_imag)
+{
+   if (bfi_real)
+   {
+      mblfr->AddTraceFaceIntegrator(bfi_real);
+   }
+   if (bfi_imag)
+   {
+      mblfi->AddTraceFaceIntegrator(bfi_imag);
+   }
+}
+
+void MixedSesquilinearForm::AddBdrTraceFaceIntegrator(BilinearFormIntegrator
+                                                      *bfi_real,
+                                                      BilinearFormIntegrator *bfi_imag)
+{
+   if (bfi_real)
+   {
+      mblfr->AddBdrTraceFaceIntegrator(bfi_real);
+   }
+   if (bfi_imag)
+   {
+      mblfi->AddBdrTraceFaceIntegrator(bfi_imag);
+   }
+}
+
+void MixedSesquilinearForm::AddBdrTraceFaceIntegrator(BilinearFormIntegrator
+                                                      *bfi_real,
+                                                      BilinearFormIntegrator *bfi_imag,
+                                                      Array<int> &bdr_marker)
+{
+   if (bfi_real)
+   {
+      mblfr->AddBdrTraceFaceIntegrator(bfi_real, bdr_marker);
+   }
+   if (bfi_imag)
+   {
+      mblfi->AddBdrTraceFaceIntegrator(bfi_imag, bdr_marker);
+   }
+}
+
+void
+MixedSesquilinearForm::Assemble(int skip_zeros)
+{
+   mblfr->Assemble(skip_zeros);
+   mblfi->Assemble(skip_zeros);
+}
+
+void
+MixedSesquilinearForm::Finalize(int skip_zeros)
+{
+   mblfr->Finalize(skip_zeros);
+   mblfi->Finalize(skip_zeros);
+}
+
+ComplexSparseMatrix *
+MixedSesquilinearForm::AssembleComplexSparseMatrix()
+{
+   return new mfem::ComplexSparseMatrix(
+             &mblfr->SpMat(), &mblfi->SpMat(), false, false, conv);
+}
+
+void
+MixedSesquilinearForm::FormRectangularLinearSystem(const Array<int> &
+                                                   ess_trial_tdof_list,
+                                                   const Array<int> & ess_test_tdof_list,
+                                                   Vector & x,
+                                                   Vector & b,
+                                                   OperatorHandle & A,
+                                                   Vector & X,
+                                                   Vector & B)
+{
+   FiniteElementSpace * fes_trial = mblfr->TrialFESpace();
+   FiniteElementSpace * fes_test = mblfr->TestFESpace();
+   const int vsize_trial = fes_trial->GetVSize();
+   const int vsize_test = fes_test->GetVSize();
+
+   // Allocate temporary Vector
+   Vector b_0;
+   b_0.UseDevice(true);
+   b_0.SetSize(vsize_test);
+   b_0 = 0.0;
+
+   // Extract the real and imaginary parts of the input Vectors
+   MFEM_ASSERT(x.Size() == 2 * vsize_trial,
+               "Input GridFunction of incorrect size!");
+   x.Read();
+   Vector x_r;
+   x_r.MakeRef(x, 0, vsize_trial);
+   Vector x_i;
+   x_i.MakeRef(x, vsize_trial, vsize_trial);
+
+   MFEM_ASSERT(b.Size() == 2 * vsize_test, "Input LinearForm of incorrect size!");
+   b.Read();
+   Vector b_r;
+   b_r.MakeRef(b, 0, vsize_test);
+   Vector b_i;
+   b_i.MakeRef(b, vsize_test, vsize_test);
+
+   if (conv == ComplexOperator::BLOCK_SYMMETRIC)
+   {
+      b_i *= -1.0;
+   }
+
+   const int tvsize_trial = fes_trial->GetTrueVSize();
+   const int tvsize_test = fes_test->GetTrueVSize();
+   OperatorHandle A_r, A_i;
+
+   X.UseDevice(true);
+   X.SetSize(2 * tvsize_trial);
+   X = 0.0;
+
+   B.UseDevice(true);
+   B.SetSize(2 * tvsize_test);
+   B = 0.0;
+
+   Vector X_r;
+   X_r.MakeRef(X, 0, tvsize_trial);
+   Vector X_i;
+   X_i.MakeRef(X, tvsize_trial, tvsize_trial);
+   Vector B_r;
+   B_r.MakeRef(B, 0, tvsize_test);
+   Vector B_i;
+   B_i.MakeRef(B, tvsize_test, tvsize_test);
+
+   Vector X_0, B_0;
+
+   if (RealInteg())
+   {
+      b_0 = b_r;
+      mblfr->FormRectangularLinearSystem(
+         ess_trial_tdof_list, ess_test_tdof_list, x_r, b_0, A_r, X_0, B_0);
+      X_r = X_0;
+      B_r = B_0;
+
+      b_0 = b_i;
+      mblfr->FormRectangularLinearSystem(
+         ess_trial_tdof_list, ess_test_tdof_list, x_i, b_0, A_r, X_0, B_0);
+      X_i = X_0;
+      B_i = B_0;
+
+      if (ImagInteg())
+      {
+         b_0 = 0.0;
+         mblfi->FormRectangularLinearSystem(
+            ess_trial_tdof_list, ess_test_tdof_list, x_i, b_0, A_i, X_0, B_0);
+         B_r -= B_0;
+
+         b_0 = 0.0;
+         mblfi->FormRectangularLinearSystem(
+            ess_trial_tdof_list, ess_test_tdof_list, x_r, b_0, A_i, X_0, B_0);
+         B_i += B_0;
+      }
+   }
+   else if (ImagInteg())
+   {
+      b_0 = b_i;
+      mblfi->FormRectangularLinearSystem(
+         ess_trial_tdof_list, ess_test_tdof_list, x_r, b_0, A_i, X_0, B_0);
+      X_r = X_0;
+      B_i = B_0;
+
+      b_0 = b_r;
+      b_0 *= -1.0;
+      mblfi->FormRectangularLinearSystem(
+         ess_trial_tdof_list, ess_test_tdof_list, x_i, b_0, A_i, X_0, B_0);
+      X_i = X_0;
+      B_r = B_0;
+      B_r *= -1.0;
+   }
+   else
+   {
+      MFEM_ABORT("Real and Imaginary part of the Mixed Sesquilinear form are empty");
+   }
+
+   if (conv == ComplexOperator::BLOCK_SYMMETRIC)
+   {
+      B_i *= -1.0;
+      b_i *= -1.0;
+   }
+
+   x_r.SyncAliasMemory(x);
+   x_i.SyncAliasMemory(x);
+   b_r.SyncAliasMemory(b);
+   b_i.SyncAliasMemory(b);
+
+   X_r.SyncAliasMemory(X);
+   X_i.SyncAliasMemory(X);
+   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_hyp =
+         new ComplexSparseMatrix(A_r.As<SparseMatrix>(),
+                                 A_i.As<SparseMatrix>(),
+                                 A_r.OwnsOperator(),
+                                 A_i.OwnsOperator(),
+                                 conv);
+      A.Reset<ComplexSparseMatrix>(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);
+}
+
+void
+MixedSesquilinearForm::FormRectangularSystemMatrix(const mfem::Array<int> &
+                                                   ess_trial_tdof_list,
+                                                   const mfem::Array<int> & ess_test_tdof_list,
+                                                   mfem::OperatorHandle & A)
+{
+   OperatorHandle A_r, A_i;
+   if (RealInteg())
+   {
+      mblfr->FormRectangularSystemMatrix(ess_trial_tdof_list, ess_test_tdof_list,
+                                         A_r);
+   }
+   if (ImagInteg())
+   {
+      mblfi->FormRectangularSystemMatrix(ess_trial_tdof_list, ess_test_tdof_list,
+                                         A_i);
+   }
+   if (!RealInteg() && !ImagInteg())
+   {
+      MFEM_ABORT("Both Real and Imaginary part of the Mixed Sesquilinear form are empty");
+   }
+
+   // 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_hyp =
+         new ComplexSparseMatrix(A_r.As<SparseMatrix>(),
+                                 A_i.As<SparseMatrix>(),
+                                 A_r.OwnsOperator(),
+                                 A_i.OwnsOperator(),
+                                 conv);
+      A.Reset<ComplexSparseMatrix>(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);
+}
+
+void
+MixedSesquilinearForm::Update()
+{
+   mblfr->Update();
+   mblfi->Update();
+}
+
 
 #ifdef MFEM_USE_MPI
 
@@ -1614,8 +2034,8 @@ ParSesquilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list,
 
    // A = A_r + i A_i
    A.Clear();
-   if ( A_r.Type() == Operator::Hypre_ParCSR ||
-        A_i.Type() == Operator::Hypre_ParCSR )
+   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>(),
@@ -1682,8 +2102,8 @@ ParSesquilinearForm::FormSystemMatrix(const Array<int> &ess_tdof_list,
 
    // A = A_r + i A_i
    A.Clear();
-   if ( A_r.Type() == Operator::Hypre_ParCSR ||
-        A_i.Type() == Operator::Hypre_ParCSR )
+   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>(),
@@ -1741,6 +2161,434 @@ ParSesquilinearForm::Update(FiniteElementSpace *nfes)
    if ( pblfi ) { pblfi->Update(nfes); }
 }
 
+bool
+ParMixedSesquilinearForm::RealInteg()
+{
+   int nint = pmblfr->GetDBFI()->Size() + pmblfr->GetBBFI()->Size() +
+              pmblfr->GetFBFI()->Size() + pmblfr->GetBFBFI()->Size() +
+              pmblfr->GetTFBFI()->Size() + pmblfr->GetBTFBFI()->Size();
+   return (nint != 0);
+}
+
+bool
+ParMixedSesquilinearForm::ImagInteg()
+{
+   int nint = pmblfi->GetDBFI()->Size() + pmblfi->GetBBFI()->Size() +
+              pmblfi->GetFBFI()->Size() + pmblfi->GetBFBFI()->Size() +
+              pmblfi->GetTFBFI()->Size() + pmblfi->GetBTFBFI()->Size();
+   return (nint != 0);
+}
+
+ParMixedSesquilinearForm::ParMixedSesquilinearForm(ParFiniteElementSpace *
+                                                   trial_fes,
+                                                   ParFiniteElementSpace * test_fes,
+                                                   ComplexOperator::Convention convention)
+   : conv(convention),
+     pmblfr(new ParMixedBilinearForm(trial_fes, test_fes)),
+     pmblfi(new ParMixedBilinearForm(trial_fes, test_fes))
+{
+}
+
+ParMixedSesquilinearForm::ParMixedSesquilinearForm(ParFiniteElementSpace *
+                                                   trial_fes,
+                                                   ParFiniteElementSpace * test_fes,
+                                                   ParMixedBilinearForm * pbfr,
+                                                   ParMixedBilinearForm * pbfi,
+                                                   ComplexOperator::Convention convention)
+   : conv(convention),
+     pmblfr(new ParMixedBilinearForm(trial_fes, test_fes, pbfr)),
+     pmblfi(new ParMixedBilinearForm(trial_fes, test_fes, pbfi))
+{
+}
+
+ParMixedSesquilinearForm::~ParMixedSesquilinearForm()
+{
+   delete pmblfr;
+   delete pmblfi;
+}
+
+void
+ParMixedSesquilinearForm::AddDomainIntegrator(BilinearFormIntegrator * bfi_real,
+                                              BilinearFormIntegrator * bfi_imag)
+{
+   if (bfi_real)
+   {
+      pmblfr->AddDomainIntegrator(bfi_real);
+   }
+   if (bfi_imag)
+   {
+      pmblfi->AddDomainIntegrator(bfi_imag);
+   }
+}
+
+void
+ParMixedSesquilinearForm::AddDomainIntegrator(BilinearFormIntegrator * bfi_real,
+                                              BilinearFormIntegrator * bfi_imag,
+                                              Array<int> & elem_marker)
+{
+   if (bfi_real)
+   {
+      pmblfr->AddDomainIntegrator(bfi_real, elem_marker);
+   }
+   if (bfi_imag)
+   {
+      pmblfi->AddDomainIntegrator(bfi_imag, elem_marker);
+   }
+}
+
+void
+ParMixedSesquilinearForm::AddBoundaryIntegrator(BilinearFormIntegrator *
+                                                bfi_real,
+                                                BilinearFormIntegrator * bfi_imag)
+{
+   if (bfi_real)
+   {
+      pmblfr->AddBoundaryIntegrator(bfi_real);
+   }
+   if (bfi_imag)
+   {
+      pmblfi->AddBoundaryIntegrator(bfi_imag);
+   }
+}
+
+void
+ParMixedSesquilinearForm::AddBoundaryIntegrator(BilinearFormIntegrator *
+                                                bfi_real,
+                                                BilinearFormIntegrator * bfi_imag,
+                                                Array<int> & bdr_marker)
+{
+   if (bfi_real)
+   {
+      pmblfr->AddBoundaryIntegrator(bfi_real, bdr_marker);
+   }
+   if (bfi_imag)
+   {
+      pmblfi->AddBoundaryIntegrator(bfi_imag, bdr_marker);
+   }
+}
+
+void
+ParMixedSesquilinearForm::AddInteriorFaceIntegrator(BilinearFormIntegrator *
+                                                    bfi_real,
+                                                    BilinearFormIntegrator * bfi_imag)
+{
+   if (bfi_real)
+   {
+      pmblfr->AddInteriorFaceIntegrator(bfi_real);
+   }
+   if (bfi_imag)
+   {
+      pmblfi->AddInteriorFaceIntegrator(bfi_imag);
+   }
+}
+
+void
+ParMixedSesquilinearForm::AddBdrFaceIntegrator(BilinearFormIntegrator *
+                                               bfi_real,
+                                               BilinearFormIntegrator * bfi_imag)
+{
+   if (bfi_real)
+   {
+      pmblfr->AddBdrFaceIntegrator(bfi_real);
+   }
+   if (bfi_imag)
+   {
+      pmblfi->AddBdrFaceIntegrator(bfi_imag);
+   }
+}
+
+void
+ParMixedSesquilinearForm::AddBdrFaceIntegrator(BilinearFormIntegrator *
+                                               bfi_real,
+                                               BilinearFormIntegrator * bfi_imag,
+                                               Array<int> & bdr_marker)
+{
+   if (bfi_real)
+   {
+      pmblfr->AddBdrFaceIntegrator(bfi_real, bdr_marker);
+   }
+   if (bfi_imag)
+   {
+      pmblfi->AddBdrFaceIntegrator(bfi_imag, bdr_marker);
+   }
+}
+
+void ParMixedSesquilinearForm::AddTraceFaceIntegrator(BilinearFormIntegrator *
+                                                      bfi_real,
+                                                      BilinearFormIntegrator * bfi_imag)
+{
+   if (bfi_real)
+   {
+      pmblfr->AddTraceFaceIntegrator(bfi_real);
+   }
+   if (bfi_imag)
+   {
+      pmblfi->AddTraceFaceIntegrator(bfi_imag);
+   }
+}
+
+void ParMixedSesquilinearForm::AddBdrTraceFaceIntegrator(
+   BilinearFormIntegrator *bfi_real,
+   BilinearFormIntegrator *bfi_imag)
+{
+   if (bfi_real)
+   {
+      pmblfr->AddBdrTraceFaceIntegrator(bfi_real);
+   }
+   if (bfi_imag)
+   {
+      pmblfi->AddBdrTraceFaceIntegrator(bfi_imag);
+   }
+}
+
+void ParMixedSesquilinearForm::AddBdrTraceFaceIntegrator(
+   BilinearFormIntegrator *bfi_real,
+   BilinearFormIntegrator *bfi_imag,
+   Array<int> &bdr_marker)
+{
+   if (bfi_real)
+   {
+      pmblfr->AddBdrTraceFaceIntegrator(bfi_real, bdr_marker);
+   }
+   if (bfi_imag)
+   {
+      pmblfi->AddBdrTraceFaceIntegrator(bfi_imag, bdr_marker);
+   }
+}
+
+void
+ParMixedSesquilinearForm::Assemble(int skip_zeros)
+{
+   pmblfr->Assemble(skip_zeros);
+   pmblfi->Assemble(skip_zeros);
+}
+
+void
+ParMixedSesquilinearForm::Finalize(int skip_zeros)
+{
+   pmblfr->Finalize(skip_zeros);
+   pmblfi->Finalize(skip_zeros);
+}
+
+ComplexHypreParMatrix *
+ParMixedSesquilinearForm::ParallelAssemble()
+{
+   return new ComplexHypreParMatrix(
+             pmblfr->ParallelAssemble(), pmblfi->ParallelAssemble(), true, true, conv);
+}
+
+void
+ParMixedSesquilinearForm::FormRectangularLinearSystem(const Array<int> &
+                                                      ess_trial_tdof_list,
+                                                      const Array<int> & ess_test_tdof_list,
+                                                      Vector & x,
+                                                      Vector & b,
+                                                      OperatorHandle & A,
+                                                      Vector & X,
+                                                      Vector & B)
+{
+   FiniteElementSpace * pfes_trial = pmblfr->TrialFESpace();
+   FiniteElementSpace * pfes_test = pmblfr->TestFESpace();
+   const int vsize_trial = pfes_trial->GetVSize();
+   const int vsize_test = pfes_test->GetVSize();
+
+   // Allocate temporary Vector
+   Vector b_0;
+   b_0.UseDevice(true);
+   b_0.SetSize(vsize_test);
+   b_0 = 0.0;
+
+   // Extract the real and imaginary parts of the input Vectors
+   MFEM_ASSERT(x.Size() == 2 * vsize_trial,
+               "Input GridFunction of incorrect size!");
+   x.Read();
+   Vector x_r;
+   x_r.MakeRef(x, 0, vsize_trial);
+   Vector x_i;
+   x_i.MakeRef(x, vsize_trial, vsize_trial);
+
+   MFEM_ASSERT(b.Size() == 2 * vsize_test, "Input LinearForm of incorrect size!");
+   b.Read();
+   Vector b_r;
+   b_r.MakeRef(b, 0, vsize_test);
+   Vector b_i;
+   b_i.MakeRef(b, vsize_test, vsize_test);
+
+   if (conv == ComplexOperator::BLOCK_SYMMETRIC)
+   {
+      b_i *= -1.0;
+   }
+
+   const int tvsize_trial = pfes_trial->GetTrueVSize();
+   const int tvsize_test = pfes_test->GetTrueVSize();
+   OperatorHandle A_r, A_i;
+
+   X.UseDevice(true);
+   X.SetSize(2 * tvsize_trial);
+   X = 0.0;
+
+   B.UseDevice(true);
+   B.SetSize(2 * tvsize_test);
+   B = 0.0;
+
+   Vector X_r;
+   X_r.MakeRef(X, 0, tvsize_trial);
+   Vector X_i;
+   X_i.MakeRef(X, tvsize_trial, tvsize_trial);
+   Vector B_r;
+   B_r.MakeRef(B, 0, tvsize_test);
+   Vector B_i;
+   B_i.MakeRef(B, tvsize_test, tvsize_test);
+
+   Vector X_0, B_0;
+
+   if (RealInteg())
+   {
+      b_0 = b_r;
+      pmblfr->FormRectangularLinearSystem(
+         ess_trial_tdof_list, ess_test_tdof_list, x_r, b_0, A_r, X_0, B_0);
+      X_r = X_0;
+      B_r = B_0;
+
+      b_0 = b_i;
+      pmblfr->FormRectangularLinearSystem(
+         ess_trial_tdof_list, ess_test_tdof_list, x_i, b_0, A_r, X_0, B_0);
+      X_i = X_0;
+      B_i = B_0;
+
+      if (ImagInteg())
+      {
+         b_0 = 0.0;
+         pmblfi->FormRectangularLinearSystem(
+            ess_trial_tdof_list, ess_test_tdof_list, x_i, b_0, A_i, X_0, B_0);
+         B_r -= B_0;
+
+         b_0 = 0.0;
+         pmblfi->FormRectangularLinearSystem(
+            ess_trial_tdof_list, ess_test_tdof_list, x_r, b_0, A_i, X_0, B_0);
+         B_i += B_0;
+      }
+   }
+   else if (ImagInteg())
+   {
+      b_0 = b_i;
+      pmblfi->FormRectangularLinearSystem(
+         ess_trial_tdof_list, ess_test_tdof_list, x_r, b_0, A_i, X_0, B_0);
+      X_r = X_0;
+      B_i = B_0;
+
+      b_0 = b_r;
+      b_0 *= -1.0;
+      pmblfi->FormRectangularLinearSystem(
+         ess_trial_tdof_list, ess_test_tdof_list, x_i, b_0, A_i, X_0, B_0);
+      X_i = X_0;
+      B_r = B_0;
+      B_r *= -1.0;
+   }
+   else
+   {
+      MFEM_ABORT("Real and Imaginary part of the Mixed Sesquilinear form are empty");
+   }
+
+   if (conv == ComplexOperator::BLOCK_SYMMETRIC)
+   {
+      B_i *= -1.0;
+      b_i *= -1.0;
+   }
+
+   x_r.SyncAliasMemory(x);
+   x_i.SyncAliasMemory(x);
+   b_r.SyncAliasMemory(b);
+   b_i.SyncAliasMemory(b);
+
+   X_r.SyncAliasMemory(X);
+   X_i.SyncAliasMemory(X);
+   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);
+}
+
+void
+ParMixedSesquilinearForm::FormRectangularSystemMatrix(const Array<int> &
+                                                      ess_trial_tdof_list,
+                                                      const Array<int> & ess_test_tdof_list,
+                                                      OperatorHandle & A)
+{
+   OperatorHandle A_r, A_i;
+   if (RealInteg())
+   {
+      pmblfr->FormRectangularSystemMatrix(ess_trial_tdof_list, ess_test_tdof_list,
+                                          A_r);
+   }
+   if (ImagInteg())
+   {
+      pmblfi->FormRectangularSystemMatrix(ess_trial_tdof_list, ess_test_tdof_list,
+                                          A_i);
+   }
+   if (!RealInteg() && !ImagInteg())
+   {
+      MFEM_ABORT("Both Real and Imaginary part of the Mixed Sesquilinear form are empty");
+   }
+
+   // 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);
+}
+
+void
+ParMixedSesquilinearForm::Update()
+{
+   pmblfr->Update();
+   pmblfi->Update();
+}
+
+
+
 #endif // MFEM_USE_MPI
 
 }

fem/complex_fem.hpp

@@ -505,6 +505,186 @@ public:
    virtual ~SesquilinearForm();
 };
 
+/** Class for a mixed sesquilinear form
+
+    A mixed sesquilinear form is a generalization of a mixed bilinear form to
+    complex-valued fields. Mixed sesquilinear forms are linear in the second
+    argument but the first argument involves a complex conjugate in the sense
+    that:
+
+                a(alpha u, beta v) = conj(alpha) beta a(u, v)
+
+    The @a convention argument in the class's constructor is documented in the
+    mfem::ComplexOperator class found in linalg/complex_operator.hpp.
+
+    When supplying integrators to the MixedSesquilinearForm either the real or
+    imaginary integrator can be NULL. This indicates that the corresponding
+    portion of the complex-valued material coefficient is equal to zero.
+*/
+class MixedSesquilinearForm
+{
+private:
+   ComplexOperator::Convention conv;
+
+   MixedBilinearForm * mblfr;
+   MixedBilinearForm * mblfi;
+
+   /* These methods check if the real/imag parts of the sesqulinear form are not
+      empty */
+   bool RealInteg();
+   bool ImagInteg();
+
+public:
+   MixedSesquilinearForm(
+      FiniteElementSpace * trial_fes,
+      FiniteElementSpace * test_fes,
+      ComplexOperator::Convention convention = ComplexOperator::HERMITIAN);
+
+   /** @brief Create a MixedSesquilinearForm on the given trial and test
+       FiniteElementSpaces, using the same integrators as the
+       MixedBilinearForms @a bfr and @a bfi.
+
+       The FiniteElementSpace pointers are not owned by the newly constructed
+       object.
+
+       The integrators are copied as pointers and they are not owned by the
+       newly constructed MixedSesquilinearForm. */
+   MixedSesquilinearForm(
+      FiniteElementSpace * trial_fes,
+      FiniteElementSpace * test_fes,
+      MixedBilinearForm * bfr,
+      MixedBilinearForm * bfi,
+      ComplexOperator::Convention convention = ComplexOperator::HERMITIAN);
+
+   ComplexOperator::Convention GetConvention() const { return conv; }
+   void SetConvention(const ComplexOperator::Convention & convention) { conv = convention; }
+
+   /// Set the desired assembly level.
+   /** Valid choices are:
+
+       - AssemblyLevel::LEGACY (default)
+       - AssemblyLevel::FULL
+       - AssemblyLevel::PARTIAL
+       - AssemblyLevel::ELEMENT
+       - AssemblyLevel::NONE
+
+       This method must be called before assembly. */
+   void SetAssemblyLevel(AssemblyLevel assembly_level)
+   {
+      mblfr->SetAssemblyLevel(assembly_level);
+      mblfi->SetAssemblyLevel(assembly_level);
+   }
+
+   MixedBilinearForm & real() { return *mblfr; }
+   MixedBilinearForm & imag() { return *mblfi; }
+   const MixedBilinearForm & real() const { return *mblfr; }
+   const MixedBilinearForm & imag() const { return *mblfi; }
+
+   /// Adds new Domain Integrator.
+   void AddDomainIntegrator(BilinearFormIntegrator * bfi_real,
+                            BilinearFormIntegrator * bfi_imag);
+
+   /// Adds new Domain Integrator, restricted to specific attributes.
+   void AddDomainIntegrator(BilinearFormIntegrator * bfi_real,
+                            BilinearFormIntegrator * bfi_imag,
+                            Array<int> & elem_marker);
+
+   /// Adds new Boundary Integrator.
+   void AddBoundaryIntegrator(BilinearFormIntegrator * bfi_real,
+                              BilinearFormIntegrator * bfi_imag);
+
+   /** @brief Adds new boundary Integrator, restricted to specific boundary
+       attributes.
+
+       Assumes ownership of @a bfi.
+
+       The mfem::array @a bdr_marker is stored internally as a pointer to the given
+       mfem::Array<int> object. */
+   void AddBoundaryIntegrator(BilinearFormIntegrator * bfi_real,
+                              BilinearFormIntegrator * bfi_imag,
+                              Array<int> & bdr_marker);
+
+   /// Adds new interior Face Integrator. Assumes ownership of @a bfi.
+   void AddInteriorFaceIntegrator(BilinearFormIntegrator * bfi_real,
+                                  BilinearFormIntegrator * bfi_imag);
+
+   /// Adds new boundary Face Integrator. Assumes ownership of @a bfi.
+   void AddBdrFaceIntegrator(BilinearFormIntegrator * bfi_real,
+                             BilinearFormIntegrator * bfi_imag);
+
+   /** @brief Adds new boundary Face Integrator, restricted to specific boundary
+       attributes.
+
+       Assumes ownership of @a bfi.
+
+       The mfem::array @a bdr_marker is stored internally as a pointer to the given
+       mfem::Array<int> object. */
+   void AddBdrFaceIntegrator(BilinearFormIntegrator * bfi_real,
+                             BilinearFormIntegrator * bfi_imag,
+                             Array<int> & bdr_marker);
+
+   /** @brief Add a trace face integrator. Assumes ownership of @a bfi.
+
+       This type of integrator assembles terms over all faces of the mesh using
+       the face FE from the trial space and the two adjacent volume FEs from
+       the test space. */
+   void AddTraceFaceIntegrator(BilinearFormIntegrator * bfi_real,
+                               BilinearFormIntegrator * bfi_imag);
+
+   /// Adds a boundary trace face integrator. Assumes ownership of @a bfi.
+   void AddBdrTraceFaceIntegrator(BilinearFormIntegrator * bfi_real,
+                                  BilinearFormIntegrator * bfi_imag);
+
+   /// Adds a boundary trace face integrator. Assumes ownership of @a bfi.
+   void AddBdrTraceFaceIntegrator(BilinearFormIntegrator * bfi_real,
+                                  BilinearFormIntegrator * bfi_imag,
+                                  Array<int> &bdr_marker);
+
+   /// Assemble the local matrix
+   void Assemble(int skip_zeros = 1);
+
+   /// Finalizes the matrix initialization.
+   void Finalize(int skip_zeros = 1);
+
+   /// Updates the internal mixed forms with the new finite element space.
+   virtual void Update();
+
+   /** @brief Return a ComplexSparseMatrix wrapping the local (L-dof) real
+       and imaginary matrices of the form.
+
+       The returned wrapper has to be deleted by the caller, but it does not
+       own the wrapped real and imaginary matrices, which remain owned by
+       this form. */
+   ComplexSparseMatrix *AssembleComplexSparseMatrix();
+
+   /// Return the trial FE space associated with the MixedSesquilinearForm.
+   FiniteElementSpace *TrialFESpace() { return mblfr->TrialFESpace(); }
+
+   /// Read-only access to the associated trial FiniteElementSpace.
+   const FiniteElementSpace *TrialFESpace() const { return mblfr->TrialFESpace(); }
+
+   /// Return the test FE space associated with the MixedSesquilinearForm.
+   FiniteElementSpace *TestFESpace() { return mblfr->TestFESpace(); }
+
+   /// Read-only access to the associated test FiniteElementSpace.
+   const FiniteElementSpace *TestFESpace() const { return mblfr->TestFESpace(); }
+
+
+   void FormRectangularLinearSystem(const Array<int> & ess_trial_tdof_list,
+                                    const Array<int> & ess_test_tdof_list,
+                                    Vector & x,
+                                    Vector & b,
+                                    OperatorHandle & A,
+                                    Vector & X,
+                                    Vector & B);
+
+   void FormRectangularSystemMatrix(const Array<int> & ess_trial_tdof_list,
+                                    const Array<int> & ess_test_tdof_list,
+                                    OperatorHandle & A);
+
+   virtual ~MixedSesquilinearForm();
+};
+
 #ifdef MFEM_USE_MPI
 
 /// Class for parallel complex-valued grid function - real + imaginary part
@@ -921,6 +1101,169 @@ public:
    virtual ~ParSesquilinearForm();
 };
 
+/** Class for a parallel mixed sesquilinear form
+
+    A mixed sesquilinear form is a generalization of a mixed bilinear form to
+    complex-valued fields. Mixed sesquilinear forms are linear in the second
+    argument but the first argument involves a complex conjugate in the sense
+    that:
+
+                a(alpha u, beta v) = conj(alpha) beta a(u, v)
+
+    The @a convention argument in the class's constructor is documented in the
+    mfem::ComplexOperator class found in linalg/complex_operator.hpp.
+
+    When supplying integrators to the ParMixedSesquilinearForm either the real
+    or imaginary integrator can be NULL. This indicates that the corresponding
+    portion of the complex-valued material coefficient is equal to zero.
+*/
+class ParMixedSesquilinearForm
+{
+private:
+   ComplexOperator::Convention conv;
+
+   ParMixedBilinearForm * pmblfr;
+   ParMixedBilinearForm * pmblfi;
+
+   /* These methods check if the real/imag parts of the sesqulinear form are
+    not empty */
+   bool RealInteg();
+   bool ImagInteg();
+
+public:
+   ParMixedSesquilinearForm(
+      ParFiniteElementSpace * trial_fes,
+      ParFiniteElementSpace * test_fes,
+      ComplexOperator::Convention convention = ComplexOperator::HERMITIAN);
+
+   /** @brief Create a ParMixedSesquilinearForm on the given trial and test
+       ParFiniteElementSpaces, using the same integrators as the
+       ParMixedBilinearForms @a pbfr and @a pbfi.
+
+       The ParFiniteElementSpace pointers are not owned by the newly
+       constructed object.
+
+       The integrators are copied as pointers and they are not owned by the
+       newly constructed ParMixedSesquilinearForm. */
+   ParMixedSesquilinearForm(
+      ParFiniteElementSpace * trial_fes,
+      ParFiniteElementSpace * test_fes,
+      ParMixedBilinearForm * pbfr,
+      ParMixedBilinearForm * pbfi,
+      ComplexOperator::Convention convention = ComplexOperator::HERMITIAN);
+
+   ComplexOperator::Convention GetConvention() const { return conv; }
+   void SetConvention(const ComplexOperator::Convention & convention) { conv = convention; }
+
+   /// Set the desired assembly level.
+   /** Valid choices are:
+
+       - AssemblyLevel::LEGACY (default)
+       - AssemblyLevel::FULL
+       - AssemblyLevel::PARTIAL
+       - AssemblyLevel::ELEMENT
+       - AssemblyLevel::NONE
+
+       This method must be called before assembly. */
+   void SetAssemblyLevel(AssemblyLevel assembly_level)
+   {
+      pmblfr->SetAssemblyLevel(assembly_level);
+      pmblfi->SetAssemblyLevel(assembly_level);
+   }
+
+   ParMixedBilinearForm & real() { return *pmblfr; }
+   ParMixedBilinearForm & imag() { return *pmblfi; }
+   const ParMixedBilinearForm & real() const { return *pmblfr; }
+   const ParMixedBilinearForm & imag() const { return *pmblfi; }
+
+   /// Adds new Domain Integrator.
+   void AddDomainIntegrator(BilinearFormIntegrator * bfi_real,
+                            BilinearFormIntegrator * bfi_imag);
+
+   /// Adds new Domain Integrator, restricted to specific attributes.
+   void AddDomainIntegrator(BilinearFormIntegrator * bfi_real,
+                            BilinearFormIntegrator * bfi_imag,
+                            Array<int> & elem_marker);
+
+   /// Adds new Boundary Integrator.
+   void AddBoundaryIntegrator(BilinearFormIntegrator * bfi_real,
+                              BilinearFormIntegrator * bfi_imag);
+
+   /** @brief Adds new boundary Integrator, restricted to specific boundary
+       attributes.
+
+       Assumes ownership of @a bfi.
+
+       The mfem::array @a bdr_marker is stored internally as a pointer to the given
+       mfem::Array<int> object. */
+   void AddBoundaryIntegrator(BilinearFormIntegrator * bfi_real,
+                              BilinearFormIntegrator * bfi_imag,
+                              Array<int> & bdr_marker);
+
+   /// Adds new interior Face Integrator. Assumes ownership of @a bfi.
+   void AddInteriorFaceIntegrator(BilinearFormIntegrator * bfi_real,
+                                  BilinearFormIntegrator * bfi_imag);
+
+   /// Adds new boundary Face Integrator. Assumes ownership of @a bfi.
+   void AddBdrFaceIntegrator(BilinearFormIntegrator * bfi_real,
+                             BilinearFormIntegrator * bfi_imag);
+
+   /** @brief Adds new boundary Face Integrator, restricted to specific boundary
+       attributes.
+
+       Assumes ownership of @a bfi.
+
+       The mfem::array @a bdr_marker is stored internally as a pointer to the given
+       mfem::Array<int> object. */
+   void AddBdrFaceIntegrator(BilinearFormIntegrator * bfi_real,
+                             BilinearFormIntegrator * bfi_imag,
+                             Array<int> & bdr_marker);
+
+   /** @brief Add a trace face integrator. Assumes ownership of @a bfi.
+
+      This type of integrator assembles terms over all faces of the mesh using
+      the face FE from the trial space and the two adjacent volume FEs from
+      the test space. */
+   void AddTraceFaceIntegrator(BilinearFormIntegrator * bfi_real,
+                               BilinearFormIntegrator * bfi_imag);
+
+   /// Adds a boundary trace face integrator. Assumes ownership of @a bfi.
+   void AddBdrTraceFaceIntegrator(BilinearFormIntegrator * bfi_real,
+                                  BilinearFormIntegrator * bfi_imag);
+
+   /// Adds a boundary trace face integrator. Assumes ownership of @a bfi.
+   void AddBdrTraceFaceIntegrator(BilinearFormIntegrator * bfi_real,
+                                  BilinearFormIntegrator * bfi_imag,
+                                  Array<int> &bdr_marker);
+
+   /// Assemble the local matrix
+   void Assemble(int skip_zeros = 1);
+
+   /// Finalizes the matrix initialization.
+   void Finalize(int skip_zeros = 1);
+
+   /// Updates the internal mixed forms with the new finite element space.
+   virtual void Update();
+
+   /// Returns the matrix assembled on the true dofs, i.e. P^t A P.
+   /** The returned matrix has to be deleted by the caller. */
+   ComplexHypreParMatrix * ParallelAssemble();
+
+   void FormRectangularLinearSystem(const Array<int> & ess_trial_tdof_list,
+                                    const Array<int> & ess_test_tdof_list,
+                                    Vector & x,
+                                    Vector & b,
+                                    OperatorHandle & A,
+                                    Vector & X,
+                                    Vector & B);
+
+   void FormRectangularSystemMatrix(const Array<int> & ess_trial_tdof_list,
+                                    const Array<int> & ess_test_tdof_list,
+                                    OperatorHandle & A);
+
+   virtual ~ParMixedSesquilinearForm();
+};
+
 #endif // MFEM_USE_MPI
 
 }

fem/fe/fe_nd.cpp

@@ -1282,12 +1282,49 @@ ND_SegmentElement::ND_SegmentElement(const int p, const int ob_type)
    }
 }
 
+void ND_SegmentElement::CalcShape(const IntegrationPoint &ip,
+                                  Vector &shape) const
+{
+   if (obasis1d.IsIntegratedType()) { obasis1d.ScaleIntegrated(false); }
+   obasis1d.Eval(ip.x, shape);
+}
+
 void ND_SegmentElement::CalcVShape(const IntegrationPoint &ip,
                                    DenseMatrix &shape) const
 {
    Vector vshape(shape.Data(), dof);
 
-   obasis1d.Eval(ip.x, vshape);
+   CalcShape(ip, vshape);
+}
+
+void ND_SegmentElement::ProjectIntegrated(VectorCoefficient &vc,
+                                          ElementTransformation &Trans,
+                                          Vector &dofs) const
+{
+   MFEM_ASSERT(obasis1d.IsIntegratedType(), "Not integrated type");
+   real_t vk[Geometry::MaxDim];
+   Vector xk(vk, vc.GetVDim());
+
+   const real_t *cp = poly1d.ClosedPoints(dof, BasisType::GaussLobatto);
+   const IntegrationRule &ir = IntRules.Get(Geometry::SEGMENT, dof);
+   IntegrationPoint ip;
+
+   for (int i = 0; i < dof; i++)
+   {
+      const real_t h = cp[i+1] - cp[i];
+      real_t val = 0.0;
+
+      for (int q = 0; q < ir.GetNPoints(); q++)
+      {
+         const IntegrationPoint &ip1d = ir.IntPoint(q);
+         ip.x = cp[i] + h*ip1d.x;
+         Trans.SetIntPoint(&ip);
+         vc.Eval(xk, Trans, ip);
+         val += ip1d.weight*Trans.Jacobian().InnerProduct(tk, vk);
+      }
+
+      dofs(i) = val*h;
+   }
 }
 
 const real_t ND_WedgeElement::tk[15] =

fem/fe/fe_nd.hpp

@@ -303,8 +303,7 @@ public:
    /** @brief Construct the ND_SegmentElement of order @a p and open
        BasisType @a ob_type */
    ND_SegmentElement(const int p, const int ob_type = BasisType::GaussLegendre);
-   void CalcShape(const IntegrationPoint &ip, Vector &shape) const override
-   { obasis1d.Eval(ip.x, shape); }
+   void CalcShape(const IntegrationPoint &ip, Vector &shape) const override;
    void CalcVShape(const IntegrationPoint &ip,
                    DenseMatrix &shape) const override;
    void CalcVShape(ElementTransformation &Trans,
@@ -325,7 +324,10 @@ public:
    using FiniteElement::Project;
    void Project(VectorCoefficient &vc,
                 ElementTransformation &Trans, Vector &dofs) const override
-   { Project_ND(tk, dof2tk, vc, Trans, dofs); }
+   {
+      if (obasis1d.IsIntegratedType()) { ProjectIntegrated(vc, Trans, dofs); }
+      else { Project_ND(tk, dof2tk, vc, Trans, dofs); }
+   }
    void ProjectMatrixCoefficient(MatrixCoefficient &mc,
                                  ElementTransformation &T,
                                  Vector &dofs) const override
@@ -338,6 +340,11 @@ public:
                     ElementTransformation &Trans,
                     DenseMatrix &grad) const override
    { ProjectGrad_ND(tk, dof2tk, fe, Trans, grad); }
+
+protected:
+   void ProjectIntegrated(VectorCoefficient &vc,
+                          ElementTransformation &Trans,
+                          Vector &dofs) const;
 };
 
 class ND_WedgeElement : public VectorFiniteElement

fem/fe_coll.hpp

@@ -100,6 +100,10 @@ public:
       return FiniteElementForGeometry(GeomType);
    }
 
+   /** @brief Returns a collection of the trace elements.
+
+       @note The collection is owned by the caller and is NOT deleted in the
+       destructor. */
    virtual FiniteElementCollection *GetTraceCollection() const;
 
    virtual ~FiniteElementCollection();

fem/integ/bilininteg_hcurl_kernels.cpp

@@ -147,18 +147,16 @@ void PAHcurlMassAssembleDiagonal3D(const int D1D,
    }); // end of element loop
 }
 
-void PAHcurlMassApply2D(const int D1D,
-                        const int Q1D,
-                        const int NE,
-                        const bool symmetric,
-                        const Array<real_t> &bo,
-                        const Array<real_t> &bc,
-                        const Array<real_t> &bot,
-                        const Array<real_t> &bct,
-                        const Vector &pa_data,
-                        const Vector &x,
-                        Vector &y)
+void PAHcurlMassApply2D(const int NE, const bool symmetric,
+                        [[maybe_unused]] const bool scalar_coeff,
+                        const Array<real_t> &bo, const Array<real_t> &bc,
+                        const Array<real_t> &bot, const Array<real_t> &bct,
+                        const Vector &pa_data, const Vector &x, Vector &y,
+                        const int D1D, [[maybe_unused]] const int TestD1D,
+                        const int Q1D)
 {
+   MFEM_ASSERT(D1D == TestD1D,
+               "Trial and Test space must have the same number of dofs");
    auto Bo = Reshape(bo.Read(), Q1D, D1D-1);
    auto Bc = Reshape(bc.Read(), Q1D, D1D);
    auto Bot = Reshape(bot.Read(), D1D-1, Q1D);
@@ -277,18 +275,16 @@ void PAHcurlMassApply2D(const int D1D,
    }); // end of element loop
 }
 
-void PAHcurlMassApply3D(const int D1D,
-                        const int Q1D,
-                        const int NE,
-                        const bool symmetric,
-                        const Array<real_t> &bo,
-                        const Array<real_t> &bc,
-                        const Array<real_t> &bot,
-                        const Array<real_t> &bct,
-                        const Vector &pa_data,
-                        const Vector &x,
-                        Vector &y)
+void PAHcurlMassApply3D(const int NE, const bool symmetric,
+                        [[maybe_unused]] const bool scalar_coeff,
+                        const Array<real_t> &bo, const Array<real_t> &bc,
+                        const Array<real_t> &bot, const Array<real_t> &bct,
+                        const Vector &pa_data, const Vector &x, Vector &y,
+                        const int D1D, [[maybe_unused]] const int TestD1D,
+                        const int Q1D)
 {
+   MFEM_VERIFY(D1D == TestD1D,
+               "Trial and test spaces must have same number of dofs");
    MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().HCURL_MAX_D1D,
                "Error: D1D > MAX_D1D");
    MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().HCURL_MAX_Q1D,

fem/integ/bilininteg_hcurl_kernels.hpp

@@ -181,228 +181,312 @@ inline void SmemPAHcurlMassAssembleDiagonal3D(const int d1d,
 }
 
 // PA H(curl) Mass Apply 2D kernel
-void PAHcurlMassApply2D(const int D1D,
-                        const int Q1D,
-                        const int NE,
-                        const bool symmetric,
-                        const Array<real_t> &bo,
-                        const Array<real_t> &bc,
-                        const Array<real_t> &bot,
-                        const Array<real_t> &bct,
-                        const Vector &pa_data,
-                        const Vector &x,
-                        Vector &y);
+void PAHcurlMassApply2D(const int NE, const bool symmetric,
+                        const bool scalar_coeff, const Array<real_t> &bo,
+                        const Array<real_t> &bc, const Array<real_t> &bot,
+                        const Array<real_t> &bct, const Vector &pa_data,
+                        const Vector &x, Vector &y, const int TrialD1D,
+                        const int TestD1D, const int Q1D);
 
 // PA H(curl) Mass Apply 3D kernel
-void PAHcurlMassApply3D(const int D1D,
-                        const int Q1D,
-                        const int NE,
-                        const bool symmetric,
-                        const Array<real_t> &bo,
-                        const Array<real_t> &bc,
-                        const Array<real_t> &bot,
-                        const Array<real_t> &bct,
-                        const Vector &pa_data,
-                        const Vector &x,
-                        Vector &y);
+void PAHcurlMassApply3D(const int NE, const bool symmetric,
+                        [[maybe_unused]] const bool scalar_coeff,
+                        const Array<real_t> &bo, const Array<real_t> &bc,
+                        const Array<real_t> &bot, const Array<real_t> &bct,
+                        const Vector &pa_data, const Vector &x, Vector &y,
+                        const int TrialD1D, [[maybe_unused]] const int TestD1D,
+                        const int Q1D);
 
 // Shared memory PA H(curl) Mass Apply 3D kernel
-template<int T_D1D = 0, int T_Q1D = 0>
-inline void SmemPAHcurlMassApply3D(const int d1d,
-                                   const int q1d,
-                                   const int NE,
-                                   const bool symmetric,
-                                   const Array<real_t> &bo,
-                                   const Array<real_t> &bc,
-                                   const Array<real_t> &bot,
-                                   const Array<real_t> &bct,
-                                   const Vector &pa_data,
-                                   const Vector &x,
-                                   Vector &y)
+template <int T_D1D = 0, int T_Q1D = 0, int TBATCH = 0, bool ACCUMULATE = true>
+inline void SmemPAHcurlMassApply3D(
+   const int NE, const bool symmetric, [[maybe_unused]] const bool scalar_coeff,
+   const Array<real_t> &bo, const Array<real_t> &bc,
+   [[maybe_unused]] const Array<real_t> &bot,
+   [[maybe_unused]] const Array<real_t> &bct, const Vector &pa_data,
+   const Vector &x, Vector &y, const int d1d = 0,
+   [[maybe_unused]] const int test_d1d = 0, const int q1d = 0)
 {
+   const int D1D = T_D1D ? T_D1D : d1d;
+   const int Q1D = T_Q1D ? T_Q1D : q1d;
+
    MFEM_VERIFY(T_D1D || d1d <= DeviceDofQuadLimits::Get().HCURL_MAX_D1D,
                "Error: d1d > HCURL_MAX_D1D");
    MFEM_VERIFY(T_Q1D || q1d <= DeviceDofQuadLimits::Get().HCURL_MAX_Q1D,
                "Error: q1d > HCURL_MAX_Q1D");
-   const int D1D = T_D1D ? T_D1D : d1d;
-   const int Q1D = T_Q1D ? T_Q1D : q1d;
-
+   MFEM_ASSERT(Q1D >= D1D, "Expected Q1D >= D1D");
    const int dataSize = symmetric ? 6 : 9;
 
-   auto Bo = Reshape(bo.Read(), Q1D, D1D-1);
-   auto Bc = Reshape(bc.Read(), Q1D, D1D);
-   auto op = Reshape(pa_data.Read(), Q1D, Q1D, Q1D, dataSize, NE);
-   auto X = Reshape(x.Read(), 3*(D1D-1)*D1D*D1D, NE);
-   auto Y = Reshape(y.ReadWrite(), 3*(D1D-1)*D1D*D1D, NE);
+   // assume trial space == test space
+   auto Bo = bo.Read();
+   auto Bc = bc.Read();
+   auto op =
+      Reshape(pa_data.Read(), Q1D, Q1D, Q1D, dataSize, NE);
+   auto X_ = Reshape(x.Read(), 3 * (D1D - 1) * D1D * D1D, NE);
+   auto y_ = y.ReadWrite();
 
-   mfem::forall_3D(NE, Q1D, Q1D, Q1D, [=] MFEM_HOST_DEVICE (int e)
+   constexpr int MD_ = T_D1D ? T_D1D : DofQuadLimits::HCURL_MAX_D1D;
+   constexpr int MQ_ = T_Q1D ? T_Q1D : DofQuadLimits::HCURL_MAX_Q1D;
+   constexpr int MDQ_ = std::max(MD_, MQ_);
+   constexpr int MB_ = TBATCH ? TBATCH : 1;
+
+   mfem::forall_2D_batch<MDQ_ * MDQ_ * MDQ_ * MB_>(
+      NE, MDQ_ * MDQ_ * MDQ_, 1, MB_, [=] MFEM_HOST_DEVICE(int e)
    {
+#if defined(__CUDA_ARCH__) || defined(__HIP_DEVICE_COMPILE__)
+      constexpr int nbz = TBATCH ? TBATCH : 1;
+      int tidz = MFEM_THREAD_ID(z);
+#else
+      constexpr int nbz = 1;
+      constexpr int tidz = 0;
+#endif
+
       constexpr int VDIM = 3;
       constexpr int MD1D = T_D1D ? T_D1D : DofQuadLimits::HCURL_MAX_D1D;
       constexpr int MQ1D = T_Q1D ? T_Q1D : DofQuadLimits::HCURL_MAX_Q1D;
-      const int D1D = T_D1D ? T_D1D : d1d;
-      const int Q1D = T_Q1D ? T_Q1D : q1d;
+      constexpr int MDQ = std::max(MD1D, MQ1D);
 
-      MFEM_SHARED real_t sBo[MQ1D][MD1D];
-      MFEM_SHARED real_t sBc[MQ1D][MD1D];
+      // nvcc limit work-around: can't have Y_ be captured first in
+      // if constexpr, so capture y_ and construct Y_ locally
+      // only works on GPU
+      auto Y = Reshape(y_, VDIM * (D1D - 1) * D1D * D1D, NE);
 
-      real_t op9[9];
-      MFEM_SHARED real_t sop[9*MQ1D*MQ1D];
-      MFEM_SHARED real_t mass[MQ1D][MQ1D][3];
+      MFEM_SHARED real_t sBo[MDQ * (MD1D - 1)];
+      MFEM_SHARED real_t sBc[MDQ * MD1D];
+      auto BO = Reshape(sBo, Q1D, D1D - 1);
+      auto BC = Reshape(sBc, Q1D, D1D);
 
-      MFEM_SHARED real_t sX[MD1D][MD1D][MD1D];
+      MFEM_SHARED real_t sX[nbz * VDIM * (MD1D - 1) * MD1D * MD1D];
+      MFEM_SHARED real_t sm0[nbz * VDIM * MDQ * MDQ * MDQ];
+      MFEM_SHARED real_t sm1[nbz * VDIM * MDQ * MDQ * MDQ];
 
-      MFEM_FOREACH_THREAD(qx,x,Q1D)
+      real_t(*X)[nbz][(MD1D - 1) * MD1D * MD1D] =
+         (real_t(*)[nbz][(MD1D - 1) * MD1D * MD1D])(sX);
+      // shapes of buffers always use MQ1D to mitigate shared memory bank
+      // conflicts
+      real_t(*DDQ)[nbz][MQ1D][MQ1D][MQ1D] =
+         (real_t(*)[nbz][MQ1D][MQ1D][MQ1D])(sm0);
+      real_t(*DQQ)[nbz][MQ1D][MQ1D][MQ1D] =
+         (real_t(*)[nbz][MQ1D][MQ1D][MQ1D])(sm1);
+      real_t(*QQQ)[nbz][MQ1D][MQ1D][MQ1D] =
+         (real_t(*)[nbz][MQ1D][MQ1D][MQ1D])(sm0);
+      real_t(*QQD)[nbz][MQ1D][MQ1D][MQ1D] =
+         (real_t(*)[nbz][MQ1D][MQ1D][MQ1D])(sm1);
+      real_t(*QDD)[nbz][MQ1D][MQ1D][MQ1D] =
+         (real_t(*)[nbz][MQ1D][MQ1D][MQ1D])(sm0);
+
+      // load dofs into smem
+      const int offset = (D1D - 1) * D1D * D1D;
+      MFEM_FOREACH_THREAD_DIRECT(ix, x, offset)
       {
-         MFEM_FOREACH_THREAD(qy,y,Q1D)
+         for (int dim = 0; dim < VDIM; ++dim)
          {
-            MFEM_FOREACH_THREAD(qz,z,Q1D)
-            {
-               for (int i=0; i<dataSize; ++i)
-               {
-                  op9[i] = op(qx,qy,qz,i,e);
-               }
-            }
+            X[dim][tidz][ix] = X_(ix + dim * offset, e);
          }
       }
-
-      const int tidx = MFEM_THREAD_ID(x);
-      const int tidy = MFEM_THREAD_ID(y);
-      const int tidz = MFEM_THREAD_ID(z);
-
+      // load basis functions data
       if (tidz == 0)
       {
-         MFEM_FOREACH_THREAD(d,y,D1D)
+         MFEM_FOREACH_THREAD_DIRECT(ix, x, D1D * Q1D) { sBc[ix] = Bc[ix]; }
+         MFEM_FOREACH_THREAD_DIRECT(ix, x, (D1D - 1) * Q1D)
          {
-            MFEM_FOREACH_THREAD(q,x,Q1D)
+            sBo[ix] = Bo[ix];
+         }
+      }
+
+      for (int dim0 = 0; dim0 < VDIM; ++dim0)
+      {
+         MFEM_SYNC_THREAD;
+         // sum factor to QQQ = Q_{dim0,dim1} B X_{dim1}
+         for (int dim1 = 0; dim1 < VDIM; ++dim1)
+         {
+            const int D1Dz = (dim1 == 2) ? D1D - 1 : D1D;
+            const int D1Dy = (dim1 == 1) ? D1D - 1 : D1D;
+            const int D1Dx = (dim1 == 0) ? D1D - 1 : D1D;
+
+            // threads assigned to mitigate bank conflicts
+            MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(qx, dy, dz, x, Q1D, D1Dy, D1Dz,
+                                                 Q1D, Q1D, Q1D)
             {
-               sBc[q][d] = Bc(q,d);
-               if (d < D1D-1)
+               real_t u = 0;
+               for (int dx = 0; dx < D1Dx; ++dx)
                {
-                  sBo[q][d] = Bo(q,d);
+                  real_t b;
+                  if (dim1 == 0)
+                  {
+                     b = BO(qx, dx);
+                  }
+                  else
+                  {
+                     b = BC(qx, dx);
+                  }
+                  u += X[dim1][tidz][dx + (dy + dz * D1Dy) * D1Dx] * b;
+               }
+               DDQ[dim1][tidz][dz][dy][qx] = u;
+            }
+         }
+         MFEM_SYNC_THREAD;
+         for (int dim1 = 0; dim1 < VDIM; ++dim1)
+         {
+            const int D1Dz = (dim1 == 2) ? D1D - 1 : D1D;
+            const int D1Dy = (dim1 == 1) ? D1D - 1 : D1D;
+            // const int D1Dx = (dim1 == 0) ? D1D - 1 : D1D;
+            // threads assigned to mitigate bank conflicts
+            MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(qx, qy, dz, x, Q1D, Q1D, D1Dz,
+                                                 Q1D, Q1D, Q1D)
+            {
+               real_t u = 0;
+               for (int dy = 0; dy < D1Dy; ++dy)
+               {
+                  real_t b;
+                  if (dim1 == 1)
+                  {
+                     b = BO(qy, dy);
+                  }
+                  else
+                  {
+                     b = BC(qy, dy);
+                  }
+                  u += DDQ[dim1][tidz][dz][dy][qx] * b;
+               }
+               DQQ[dim1][tidz][dz][qy][qx] = u;
+            }
+         }
+         MFEM_SYNC_THREAD;
+         for (int dim1 = 0; dim1 < VDIM; ++dim1)
+         {
+            const int D1Dz = (dim1 == 2) ? D1D - 1 : D1D;
+            // const int D1Dy = (dim1 == 1) ? D1D - 1 : D1D;
+            // const int D1Dx = (dim1 == 0) ? D1D - 1 : D1D;
+            MFEM_FOREACH_THREAD_DIRECT_3D(qx, qy, qz, x, Q1D, Q1D, Q1D)
+            {
+               real_t u = 0;
+               for (int dz = 0; dz < D1Dz; ++dz)
+               {
+                  real_t b;
+                  if (dim1 == 2)
+                  {
+                     b = BO(qz, dz);
+                  }
+                  else
+                  {
+                     b = BC(qz, dz);
+                  }
+                  u += DQQ[dim1][tidz][dz][qy][qx] * b;
+               }
+               // pa_data is row major
+               int idx;
+               if (symmetric)
+               {
+                  int row;
+                  int col;
+                  if (dim0 > dim1)
+                  {
+                     row = dim1;
+                     col = dim0;
+                  }
+                  else
+                  {
+                     row = dim0;
+                     col = dim1;
+                  }
+                  idx = col + VDIM * row - row * (row + 1) / 2;
+               }
+               else
+               {
+                  idx = dim0 * VDIM + dim1;
+               }
+               QQQ[dim1][tidz][qz][qy][qx] = op(qx, qy, qz, idx, e) * u;
+            }
+         }
+         MFEM_SYNC_THREAD;
+         // sum factor back to Y
+         // Assume bot and bct == bo^t and bc^t respectively (i.e. test ==
+         // trial functions), skip loading them again.
+         {
+            const int D1Dz = (dim0 == 2) ? D1D - 1 : D1D;
+            const int D1Dy = (dim0 == 1) ? D1D - 1 : D1D;
+            const int D1Dx = (dim0 == 0) ? D1D - 1 : D1D;
+            // threads assigned to mitigate bank conflicts
+            MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(dz, qx, qy, x, D1Dz, Q1D, Q1D,
+                                                 Q1D, Q1D, Q1D)
+            {
+               for (int dim1 = 0; dim1 < VDIM; ++dim1)
+               {
+                  real_t u = 0;
+                  for (int qz = 0; qz < Q1D; ++qz)
+                  {
+                     real_t b = 0;
+                     if (dim0 == 2)
+                     {
+                        b = BO(qz, dz);
+                     }
+                     else
+                     {
+                        b = BC(qz, dz);
+                     }
+                     u += QQQ[dim1][tidz][qz][qy][qx] * b;
+                  }
+                  QQD[dim1][tidz][qy][qx][dz] = u;
+               }
+            }
+            MFEM_SYNC_THREAD;
+            // threads assigned to mitigate bank conflicts
+            MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(dy, dz, qx, x, D1Dy, D1Dz, Q1D,
+                                                 Q1D, Q1D, Q1D)
+            {
+               for (int dim1 = 0; dim1 < VDIM; ++dim1)
+               {
+                  real_t u = 0;
+                  for (int qy = 0; qy < Q1D; ++qy)
+                  {
+                     real_t b;
+                     if (dim0 == 1)
+                     {
+                        b = BO(qy, dy);
+                     }
+                     else
+                     {
+                        b = BC(qy, dy);
+                     }
+                     u += QQD[dim1][tidz][qy][qx][dz] * b;
+                  }
+                  QDD[dim1][tidz][qx][dz][dy] = u;
+               }
+            }
+            MFEM_SYNC_THREAD;
+
+            MFEM_FOREACH_THREAD_DIRECT_3D(dx, dy, dz, x, D1Dx, D1Dy, D1Dz)
+            {
+               int ix = dx + D1Dx * (dy + D1Dy * dz);
+               real_t u = 0;
+               for (int qx = 0; qx < Q1D; ++qx)
+               {
+                  real_t b;
+                  if (dim0 == 0)
+                  {
+                     b = BO(qx, dx);
+                  }
+                  else
+                  {
+                     b = BC(qx, dx);
+                  }
+                  for (int dim1 = 0; dim1 < VDIM; ++dim1)
+                  {
+                     u += QDD[dim1][tidz][qx][dz][dy] * b;
+                  }
+               }
+               if constexpr (ACCUMULATE)
+               {
+                  Y(ix + dim0 * offset, e) += u;
+               }
+               else
+               {
+                  Y(ix + dim0 * offset, e) = u;
                }
             }
          }
       }
-      MFEM_SYNC_THREAD;
-
-      for (int qz=0; qz < Q1D; ++qz)
-      {
-         int osc = 0;
-         for (int c = 0; c < VDIM; ++c)  // loop over x, y, z components
-         {
-            const int D1Dz = (c == 2) ? D1D - 1 : D1D;
-            const int D1Dy = (c == 1) ? D1D - 1 : D1D;
-            const int D1Dx = (c == 0) ? D1D - 1 : D1D;
-
-            MFEM_FOREACH_THREAD(dz,z,D1Dz)
-            {
-               MFEM_FOREACH_THREAD(dy,y,D1Dy)
-               {
-                  MFEM_FOREACH_THREAD(dx,x,D1Dx)
-                  {
-                     sX[dz][dy][dx] = X(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e);
-                  }
-               }
-            }
-            MFEM_SYNC_THREAD;
-
-            if (tidz == qz)
-            {
-               for (int i=0; i<dataSize; ++i)
-               {
-                  sop[i + (dataSize*tidx) + (dataSize*Q1D*tidy)] = op9[i];
-               }
-
-               MFEM_FOREACH_THREAD(qy,y,Q1D)
-               {
-                  MFEM_FOREACH_THREAD(qx,x,Q1D)
-                  {
-                     real_t u = 0.0;
-
-                     for (int dz = 0; dz < D1Dz; ++dz)
-                     {
-                        const real_t wz = (c == 2) ? sBo[qz][dz] : sBc[qz][dz];
-                        for (int dy = 0; dy < D1Dy; ++dy)
-                        {
-                           const real_t wy = (c == 1) ? sBo[qy][dy] : sBc[qy][dy];
-                           for (int dx = 0; dx < D1Dx; ++dx)
-                           {
-                              const real_t t = sX[dz][dy][dx];
-                              const real_t wx = (c == 0) ? sBo[qx][dx] : sBc[qx][dx];
-                              u += t * wx * wy * wz;
-                           }
-                        }
-                     }
-
-                     mass[qy][qx][c] = u;
-                  } // qx
-               } // qy
-            } // tidz == qz
-
-            osc += D1Dx * D1Dy * D1Dz;
-            MFEM_SYNC_THREAD;
-         } // c
-
-         MFEM_SYNC_THREAD;  // Sync mass[qy][qx][d] and sop
-
-         osc = 0;
-         for (int c = 0; c < VDIM; ++c)  // loop over x, y, z components
-         {
-            const int D1Dz = (c == 2) ? D1D - 1 : D1D;
-            const int D1Dy = (c == 1) ? D1D - 1 : D1D;
-            const int D1Dx = (c == 0) ? D1D - 1 : D1D;
-
-            real_t dxyz = 0.0;
-
-            MFEM_FOREACH_THREAD(dz,z,D1Dz)
-            {
-               const real_t wz = (c == 2) ? sBo[qz][dz] : sBc[qz][dz];
-
-               MFEM_FOREACH_THREAD(dy,y,D1Dy)
-               {
-                  MFEM_FOREACH_THREAD(dx,x,D1Dx)
-                  {
-                     for (int qy = 0; qy < Q1D; ++qy)
-                     {
-                        const real_t wy = (c == 1) ? sBo[qy][dy] : sBc[qy][dy];
-                        for (int qx = 0; qx < Q1D; ++qx)
-                        {
-                           const int os = (dataSize*qx) + (dataSize*Q1D*qy);
-                           const int id1 = os + ((c == 0) ? 0 : ((c == 1) ? (symmetric ? 1 : 3) :
-                                                                 (symmetric ? 2 : 6))); // O11, O21, O31
-                           const int id2 = os + ((c == 0) ? 1 : ((c == 1) ? (symmetric ? 3 : 4) :
-                                                                 (symmetric ? 4 : 7))); // O12, O22, O32
-                           const int id3 = os + ((c == 0) ? 2 : ((c == 1) ? (symmetric ? 4 : 5) :
-                                                                 (symmetric ? 5 : 8))); // O13, O23, O33
-
-                           const real_t m_c = (sop[id1] * mass[qy][qx][0]) + (sop[id2] * mass[qy][qx][1]) +
-                                              (sop[id3] * mass[qy][qx][2]);
-
-                           const real_t wx = (c == 0) ? sBo[qx][dx] : sBc[qx][dx];
-                           dxyz += m_c * wx * wy * wz;
-                        }
-                     }
-                  }
-               }
-            }
-
-            MFEM_SYNC_THREAD;
-
-            MFEM_FOREACH_THREAD(dz,z,D1Dz)
-            {
-               MFEM_FOREACH_THREAD(dy,y,D1Dy)
-               {
-                  MFEM_FOREACH_THREAD(dx,x,D1Dx)
-                  {
-                     Y(dx + ((dy + (dz * D1Dy)) * D1Dx) + osc, e) += dxyz;
-                  }
-               }
-            }
-
-            osc += D1Dx * D1Dy * D1Dz;
-         } // c loop
-      } // qz
    }); // end of element loop
 }
 

fem/integ/bilininteg_hcurlhdiv_kernels.hpp

@@ -62,6 +62,30 @@ void PAHcurlHdivMassApply2D(const int D1D,
                             const Vector &x_,
                             Vector &y_);
 
+/// H(curl) test, H(div) trial
+inline void
+PAHcurlHdivMassApply2D(const int NE, const bool, const bool scalarCoeff,
+                       const Array<real_t> &Bo_, const Array<real_t> &Bc_,
+                       const Array<real_t> &Bot_, const Array<real_t> &Bct_,
+                       const Vector &op_, const Vector &x_, Vector &y_,
+                       const int D1D, const int D1Dtest, const int Q1D)
+{
+   return PAHcurlHdivMassApply2D(D1D, D1Dtest, Q1D, NE, scalarCoeff, false,
+                                 false, Bo_, Bc_, Bot_, Bct_, op_, x_, y_);
+}
+
+/// H(div) test, H(curl) trial
+inline void
+PAHdivHcurlMassApply2D(const int NE, const bool, const bool scalarCoeff,
+                       const Array<real_t> &Bo_, const Array<real_t> &Bc_,
+                       const Array<real_t> &Bot_, const Array<real_t> &Bct_,
+                       const Vector &op_, const Vector &x_, Vector &y_,
+                       const int D1D, const int D1Dtest, const int Q1D)
+{
+   return PAHcurlHdivMassApply2D(D1D, D1Dtest, Q1D, NE, scalarCoeff, true,
+                                 false, Bo_, Bc_, Bot_, Bct_, op_, x_, y_);
+}
+
 // PA H(curl)-H(div) Mass Apply 3D kernel
 void PAHcurlHdivMassApply3D(const int D1D,
                             const int D1Dtest,
@@ -78,6 +102,30 @@ void PAHcurlHdivMassApply3D(const int D1D,
                             const Vector &x_,
                             Vector &y_);
 
+/// H(curl) test, H(div) trial
+inline void
+PAHcurlHdivMassApply3D(const int NE, const bool, const bool scalarCoeff,
+                       const Array<real_t> &Bo_, const Array<real_t> &Bc_,
+                       const Array<real_t> &Bot_, const Array<real_t> &Bct_,
+                       const Vector &op_, const Vector &x_, Vector &y_,
+                       const int D1D, const int D1Dtest, const int Q1D)
+{
+   PAHcurlHdivMassApply3D(D1D, D1Dtest, Q1D, NE, scalarCoeff, false, false, Bo_,
+                          Bc_, Bot_, Bct_, op_, x_, y_);
+}
+
+/// H(div) test, H(curl) trial
+inline void
+PAHdivHcurlMassApply3D(const int NE, const bool, const bool scalarCoeff,
+                       const Array<real_t> &Bo_, const Array<real_t> &Bc_,
+                       const Array<real_t> &Bot_, const Array<real_t> &Bct_,
+                       const Vector &op_, const Vector &x_, Vector &y_,
+                       const int D1D, const int D1Dtest, const int Q1D)
+{
+   PAHcurlHdivMassApply3D(D1D, D1Dtest, Q1D, NE, scalarCoeff, true, false, Bo_,
+                          Bc_, Bot_, Bct_, op_, x_, y_);
+}
+
 // PA H(curl)-H(div) Curl Apply 3D kernel
 template<int T_D1D = 0, int T_D1D_TEST = 0, int T_Q1D = 0>
 inline void PAHcurlHdivApply3D(const int d1d,

fem/integ/bilininteg_hdiv_kernels.cpp

@@ -294,61 +294,14 @@ void PAHdivMassAssembleDiagonal3D(const int D1D,
    }); // end of element loop
 }
 
-void PAHdivMassApply(const int dim,
-                     const int D1D,
-                     const int Q1D,
-                     const int NE,
-                     const bool symmetric,
-                     const Array<real_t> &Bo,
-                     const Array<real_t> &Bc,
-                     const Array<real_t> &Bot,
-                     const Array<real_t> &Bct,
-                     const Vector &op,
-                     const Vector &x,
-                     Vector &y)
-{
-   const int id = (D1D << 4) | Q1D;
-
-   if (dim == 2)
-   {
-      switch (id)
-      {
-         case 0x22: return SmemPAHdivMassApply2D<2,2>(NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
-         case 0x33: return SmemPAHdivMassApply2D<3,3>(NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
-         case 0x44: return SmemPAHdivMassApply2D<4,4>(NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
-         case 0x55: return SmemPAHdivMassApply2D<5,5>(NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
-         default: // fallback
-            return PAHdivMassApply2D(D1D,Q1D,NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
-      }
-   }
-   else if (dim == 3)
-   {
-      switch (id)
-      {
-         case 0x23: return SmemPAHdivMassApply3D<2,3>(NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
-         case 0x34: return SmemPAHdivMassApply3D<3,4>(NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
-         case 0x45: return SmemPAHdivMassApply3D<4,5>(NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
-         case 0x56: return SmemPAHdivMassApply3D<5,6>(NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
-         case 0x67: return SmemPAHdivMassApply3D<6,7>(NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
-         case 0x78: return SmemPAHdivMassApply3D<7,8>(NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
-         default: // fallback
-            return PAHdivMassApply3D(D1D,Q1D,NE,symmetric,Bo,Bc,Bot,Bct,op,x,y);
-      }
-   }
-}
-
-void PAHdivMassApply2D(const int D1D,
-                       const int Q1D,
-                       const int NE,
-                       const bool symmetric,
-                       const Array<real_t> &Bo_,
-                       const Array<real_t> &Bc_,
-                       const Array<real_t> &Bot_,
-                       const Array<real_t> &Bct_,
-                       const Vector &op_,
-                       const Vector &x_,
-                       Vector &y_)
+void PAHdivMassApply2D(const int NE, const bool symmetric, const bool,
+                       const Array<real_t> &Bo_, const Array<real_t> &Bc_,
+                       const Array<real_t> &Bot_, const Array<real_t> &Bct_,
+                       const Vector &op_, const Vector &x_, Vector &y_,
+                       const int D1D, const int TestD1D, const int Q1D)
 {
+   MFEM_VERIFY(D1D == TestD1D,
+               "Trial and test spaces must have same number of dofs");
    auto Bo = Reshape(Bo_.Read(), Q1D, D1D-1);
    auto Bc = Reshape(Bc_.Read(), Q1D, D1D);
    auto Bot = Reshape(Bot_.Read(), D1D-1, Q1D);
@@ -468,18 +421,14 @@ void PAHdivMassApply2D(const int D1D,
    }); // end of element loop
 }
 
-void PAHdivMassApply3D(const int D1D,
-                       const int Q1D,
-                       const int NE,
-                       const bool symmetric,
-                       const Array<real_t> &Bo_,
-                       const Array<real_t> &Bc_,
-                       const Array<real_t> &Bot_,
-                       const Array<real_t> &Bct_,
-                       const Vector &op_,
-                       const Vector &x_,
-                       Vector &y_)
+void PAHdivMassApply3D(const int NE, const bool symmetric, const bool,
+                       const Array<real_t> &Bo_, const Array<real_t> &Bc_,
+                       const Array<real_t> &Bot_, const Array<real_t> &Bct_,
+                       const Vector &op_, const Vector &x_, Vector &y_,
+                       const int D1D, const int TestD1D, const int Q1D)
 {
+   MFEM_VERIFY(D1D == TestD1D,
+               "Trial and test spaces must have same number of dofs");
    MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().HDIV_MAX_D1D,
                "Error: D1D > HDIV_MAX_D1D");
    MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().HDIV_MAX_Q1D,

fem/integ/bilininteg_hdiv_kernels.hpp

@@ -66,58 +66,29 @@ void PAHdivMassAssembleDiagonal3D(const int D1D,
                                   const Vector &op_,
                                   Vector &diag_);
 
-void PAHdivMassApply(const int dim,
-                     const int D1D,
-                     const int Q1D,
-                     const int NE,
-                     const bool symmetric,
-                     const Array<real_t> &Bo,
-                     const Array<real_t> &Bc,
-                     const Array<real_t> &Bot,
-                     const Array<real_t> &Bct,
-                     const Vector &op,
-                     const Vector &x,
-                     Vector &y);
-
 // PA H(div) Mass Apply 2D kernel
-void PAHdivMassApply2D(const int D1D,
-                       const int Q1D,
-                       const int NE,
-                       const bool symmetric,
-                       const Array<real_t> &Bo_,
-                       const Array<real_t> &Bc_,
-                       const Array<real_t> &Bot_,
-                       const Array<real_t> &Bct_,
-                       const Vector &op_,
-                       const Vector &x_,
-                       Vector &y_);
+void PAHdivMassApply2D(const int NE, const bool symmetric,
+                       const bool scalar_coeff, const Array<real_t> &Bo_,
+                       const Array<real_t> &Bc_, const Array<real_t> &Bot_,
+                       const Array<real_t> &Bct_, const Vector &op_,
+                       const Vector &x_, Vector &y_, const int D1D,
+                       const int TestD1D, const int Q1D);
 
 // PA H(div) Mass Apply 3D kernel
-void PAHdivMassApply3D(const int D1D,
-                       const int Q1D,
-                       const int NE,
-                       const bool symmetric,
-                       const Array<real_t> &Bo_,
-                       const Array<real_t> &Bc_,
-                       const Array<real_t> &Bot_,
-                       const Array<real_t> &Bct_,
-                       const Vector &op_,
-                       const Vector &x_,
-                       Vector &y_);
+void PAHdivMassApply3D(const int NE, const bool symmetric,
+                       const bool scalar_coeff, const Array<real_t> &Bo_,
+                       const Array<real_t> &Bc_, const Array<real_t> &Bot_,
+                       const Array<real_t> &Bct_, const Vector &op_,
+                       const Vector &x_, Vector &y_, const int D1D,
+                       const int TestD1D, const int Q1D);
 
 // Shared memory PA H(div) Mass Apply 2D kernel
-template<int T_D1D = 0, int T_Q1D = 0>
-inline void SmemPAHdivMassApply2D(const int NE,
-                                  const bool symmetric,
-                                  const Array<real_t> &Bo_,
-                                  const Array<real_t> &Bc_,
-                                  const Array<real_t> &Bot_,
-                                  const Array<real_t> &Bct_,
-                                  const Vector &op_,
-                                  const Vector &x_,
-                                  Vector &y_,
-                                  const int d1d = 0,
-                                  const int q1d = 0)
+template <int T_D1D = 0, int T_Q1D = 0>
+inline void SmemPAHdivMassApply2D(
+   const int NE, const bool symmetric, const bool, const Array<real_t> &Bo_,
+   const Array<real_t> &Bc_, const Array<real_t> &Bot_,
+   const Array<real_t> &Bct_, const Vector &op_, const Vector &x_, Vector &y_,
+   const int d1d = 0, const int = 0, const int q1d = 0)
 {
    MFEM_CONTRACT_VAR(Bot_);
    MFEM_CONTRACT_VAR(Bct_);
@@ -280,18 +251,13 @@ inline void SmemPAHdivMassApply2D(const int NE,
 }
 
 // Shared memory PA H(div) Mass Apply 3D kernel
-template<int T_D1D = 0, int T_Q1D = 0>
-inline void SmemPAHdivMassApply3D(const int NE,
-                                  const bool symmetric,
-                                  const Array<real_t> &Bo_,
-                                  const Array<real_t> &Bc_,
-                                  const Array<real_t> &Bot_,
-                                  const Array<real_t> &Bct_,
-                                  const Vector &op_,
-                                  const Vector &x_,
-                                  Vector &y_,
-                                  const int d1d = 0,
-                                  const int q1d = 0)
+template <int T_D1D = 0, int T_Q1D = 0>
+inline void
+SmemPAHdivMassApply3D(const int NE, const bool symmetric, const bool,
+                      const Array<real_t> &Bo_, const Array<real_t> &Bc_,
+                      const Array<real_t> &Bot_, const Array<real_t> &Bct_,
+                      const Vector &op_, const Vector &x_, Vector &y_,
+                      const int d1d = 0, const int = 0, const int q1d = 0)
 {
    MFEM_CONTRACT_VAR(Bot_);
    MFEM_CONTRACT_VAR(Bct_);

fem/integ/bilininteg_vectorfemass_kernels.hpp

@@ -0,0 +1,113 @@
+// 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_BILININTEG_VECTORFEMASS_KERNELS_HPP
+#define MFEM_BILININTEG_VECTORFEMASS_KERNELS_HPP
+
+#include "../../config/config.hpp"
+#include "../bilininteg.hpp"
+
+#include "bilininteg_diffusion_kernels.hpp"
+#include "bilininteg_hcurl_kernels.hpp"
+#include "bilininteg_hdiv_kernels.hpp"
+#include "bilininteg_hcurlhdiv_kernels.hpp"
+
+namespace mfem
+{
+/// \cond DO_NOT_DOCUMENT
+namespace internal
+{
+namespace hcurlmass
+{
+constexpr int NBZ3D(int d1d, int q1d)
+{
+   if (d1d <= 1 || q1d <= 0)
+   {
+      return 1;
+   }
+   // assume q1d >= d1d
+   // z dimension is capped at 64 on nvidia and amd gpus
+   int tmp = std::min((128 + q1d * q1d * q1d - 1) / (q1d * q1d * q1d), 64);
+   int smem_req =
+      sizeof(mfem::real_t) *
+      (3 * ((d1d - 1) * d1d * d1d + 2 * q1d * q1d * q1d) * tmp +
+       q1d * (d1d - 1) + q1d * d1d);
+   // assume GPU has at least 48k shared memory
+   return std::max(std::min(tmp, (48 * 1024 + smem_req - 1) / smem_req), 1);
+}
+} // namespace hcurlmass
+} // namespace internal
+
+template <FiniteElement::DerivType TrialType, FiniteElement::DerivType TestType,
+          int DIM, int TrialD1D, int TestD1D, int Q1D>
+VectorFEMassIntegrator::ApplyKernelType
+VectorFEMassIntegrator::ApplyPAKernels::Kernel()
+{
+   constexpr bool trial_curl = (TrialType == mfem::FiniteElement::CURL);
+   constexpr bool trial_div = (TrialType == mfem::FiniteElement::DIV);
+   constexpr bool test_curl = (TestType == mfem::FiniteElement::CURL);
+   constexpr bool test_div = (TestType == mfem::FiniteElement::DIV);
+
+   if constexpr (DIM == 3)
+   {
+      if constexpr (trial_curl && test_curl)
+      {
+         if (Device::Allows(Backend::DEVICE_MASK))
+         {
+            // assume TrialD1D == TestD1D
+            return internal::SmemPAHcurlMassApply3D<
+                   TrialD1D, Q1D, internal::hcurlmass::NBZ3D(TrialD1D, Q1D)>;
+         }
+         else
+         {
+            return internal::PAHcurlMassApply3D;
+         }
+      }
+      else if constexpr (trial_div && test_div)
+      {
+         // assumes TrialD1D == TestD1D
+         return internal::SmemPAHdivMassApply3D<TrialD1D, Q1D>;
+      }
+      else if constexpr (trial_curl && test_div)
+      {
+         return internal::PAHdivHcurlMassApply3D;
+      }
+      else if constexpr (trial_div && test_curl)
+      {
+         return internal::PAHcurlHdivMassApply3D;
+      }
+   }
+   else if constexpr (DIM == 2) // 2D
+   {
+      if constexpr (trial_curl && test_curl)
+      {
+         return internal::PAHcurlMassApply2D;
+      }
+      else if constexpr (trial_div && test_div)
+      {
+         // assumes TrialD1D == TestD1D
+         return internal::SmemPAHdivMassApply2D<TrialD1D, Q1D>;
+      }
+      else if constexpr (trial_curl && test_div)
+      {
+         return internal::PAHdivHcurlMassApply2D;
+      }
+      else if constexpr (trial_div && test_curl)
+      {
+         return internal::PAHcurlHdivMassApply2D;
+      }
+   }
+   MFEM_ABORT("Unknown kernel.");
+}
+/// \endcond DO_NOT_DOCUMENT
+}
+
+#endif

fem/integ/bilininteg_vectorfemass_pa.cpp

@@ -10,15 +10,123 @@
 // CONTRIBUTING.md for details.
 
 #include "../bilininteg.hpp"
-#include "../gridfunc.hpp"
-#include "../qfunction.hpp"
-#include "bilininteg_diffusion_kernels.hpp"
-#include "bilininteg_hcurl_kernels.hpp"
-#include "bilininteg_hdiv_kernels.hpp"
-#include "bilininteg_hcurlhdiv_kernels.hpp"
+#include "bilininteg_vectorfemass_kernels.hpp"
 
 namespace mfem
 {
+/// \cond DO_NOT_DOCUMENT
+VectorFEMassIntegrator::ApplyKernelType
+VectorFEMassIntegrator::ApplyPAKernels::Fallback(
+   FiniteElement::DerivType TrialType, FiniteElement::DerivType TestType,
+   int dim, int, int, int)
+{
+   const bool trial_curl = (TrialType == mfem::FiniteElement::CURL);
+   const bool trial_div = (TrialType == mfem::FiniteElement::DIV);
+   const bool test_curl = (TestType == mfem::FiniteElement::CURL);
+   const bool test_div = (TestType == mfem::FiniteElement::DIV);
+
+   if (dim == 3)
+   {
+      if (trial_curl && test_curl)
+      {
+         return internal::PAHcurlMassApply3D;
+      }
+      else if (trial_div && test_div)
+      {
+         return internal::PAHdivMassApply3D;
+      }
+      else if (trial_curl && test_div)
+      {
+         return internal::PAHdivHcurlMassApply3D;
+      }
+      else if (trial_div && test_curl)
+      {
+         return internal::PAHcurlHdivMassApply3D;
+      }
+   }
+   else if (dim == 2) // 2D
+   {
+      if (trial_curl && test_curl)
+      {
+         return internal::PAHcurlMassApply2D;
+      }
+      else if (trial_div && test_div)
+      {
+         return internal::PAHdivMassApply2D;
+      }
+      else if (trial_curl && test_div)
+      {
+         return internal::PAHdivHcurlMassApply2D;
+      }
+      else if (trial_div && test_curl)
+      {
+         return internal::PAHcurlHdivMassApply2D;
+      }
+   }
+   MFEM_ABORT("Unknown kernel.");
+}
+/// \endcond DO_NOT_DOCUMENT
+
+VectorFEMassIntegrator::Kernels::Kernels()
+{
+   // h(curl), h(curl)
+   // Q = P + 1 (3D)
+   VectorFEMassIntegrator::AddSpecialization<FiniteElement::CURL,
+                          FiniteElement::CURL, 3, 2, 2, 3>();
+   VectorFEMassIntegrator::AddSpecialization<FiniteElement::CURL,
+                          FiniteElement::CURL, 3, 3, 3, 4>();
+   VectorFEMassIntegrator::AddSpecialization<FiniteElement::CURL,
+                          FiniteElement::CURL, 3, 4, 4, 5>();
+   VectorFEMassIntegrator::AddSpecialization<FiniteElement::CURL,
+                          FiniteElement::CURL, 3, 5, 5, 6>();
+   // Q = P + 2 (3D)
+   VectorFEMassIntegrator::AddSpecialization<FiniteElement::CURL,
+                          FiniteElement::CURL, 3, 2, 2, 4>();
+   VectorFEMassIntegrator::AddSpecialization<FiniteElement::CURL,
+                          FiniteElement::CURL, 3, 3, 3, 5>();
+   VectorFEMassIntegrator::AddSpecialization<FiniteElement::CURL,
+                          FiniteElement::CURL, 3, 4, 4, 6>();
+   VectorFEMassIntegrator::AddSpecialization<FiniteElement::CURL,
+                          FiniteElement::CURL, 3, 5, 5, 7>();
+   // Q = P + 4 (3D)
+   VectorFEMassIntegrator::AddSpecialization<FiniteElement::CURL,
+                          FiniteElement::CURL, 3, 2, 2, 6>();
+   VectorFEMassIntegrator::AddSpecialization<FiniteElement::CURL,
+                          FiniteElement::CURL, 3, 3, 3, 7>();
+   VectorFEMassIntegrator::AddSpecialization<FiniteElement::CURL,
+                          FiniteElement::CURL, 3, 4, 4, 8>();
+   VectorFEMassIntegrator::AddSpecialization<FiniteElement::CURL,
+                          FiniteElement::CURL, 3, 5, 5, 9>();
+   // h(div), h(div)
+   // Q = P (2D)
+   VectorFEMassIntegrator::AddSpecialization<FiniteElement::DIV,
+                          FiniteElement::DIV, 2, 2, 2, 2>();
+   VectorFEMassIntegrator::AddSpecialization<FiniteElement::DIV,
+                          FiniteElement::DIV, 2, 3, 3, 3>();
+   VectorFEMassIntegrator::AddSpecialization<FiniteElement::DIV,
+                          FiniteElement::DIV, 2, 4, 4, 4>();
+   VectorFEMassIntegrator::AddSpecialization<FiniteElement::DIV,
+                          FiniteElement::DIV, 2, 5, 5, 5>();
+
+   // Q = P + 1 (3D)
+   VectorFEMassIntegrator::AddSpecialization<FiniteElement::DIV,
+                          FiniteElement::DIV, 3, 2, 2, 3>();
+   VectorFEMassIntegrator::AddSpecialization<FiniteElement::DIV,
+                          FiniteElement::DIV, 3, 3, 3, 4>();
+   VectorFEMassIntegrator::AddSpecialization<FiniteElement::DIV,
+                          FiniteElement::DIV, 3, 4, 4, 5>();
+   VectorFEMassIntegrator::AddSpecialization<FiniteElement::DIV,
+                          FiniteElement::DIV, 3, 5, 5, 6>();
+}
+
+void VectorFEMassIntegrator::Init(Coefficient *q, DiagonalMatrixCoefficient *dq,
+                                  MatrixCoefficient *mq)
+{
+   static Kernels kernels{};
+   Q = q;
+   DQ = dq;
+   MQ = mq;
+}
 
 void VectorFEMassIntegrator::AssemblePA(const FiniteElementSpace &fes)
 {
@@ -67,8 +175,8 @@ void VectorFEMassIntegrator::AssemblePA(const FiniteElementSpace &trial_fes,
 
    MFEM_VERIFY(dofs1D == mapsO->ndof + 1 && quad1D == mapsO->nqpt, "");
 
-   trial_fetype = trial_el->GetDerivType();
-   test_fetype = test_el->GetDerivType();
+   trial_fetype = static_cast<FiniteElement::DerivType>(trial_el->GetDerivType());
+   test_fetype = static_cast<FiniteElement::DerivType>(test_el->GetDerivType());
 
    const bool trial_curl = (trial_fetype == mfem::FiniteElement::CURL);
    const bool trial_div = (trial_fetype == mfem::FiniteElement::DIV);
@@ -215,225 +323,34 @@ void VectorFEMassIntegrator::AssembleDiagonalPA(Vector& diag)
 
 void VectorFEMassIntegrator::AddMultPA(const Vector &x, Vector &y) const
 {
-   const bool trial_curl = (trial_fetype == mfem::FiniteElement::CURL);
-   const bool trial_div = (trial_fetype == mfem::FiniteElement::DIV);
-   const bool test_curl = (test_fetype == mfem::FiniteElement::CURL);
-   const bool test_div = (test_fetype == mfem::FiniteElement::DIV);
-
-   if (dim == 3)
-   {
-      if (trial_curl && test_curl)
-      {
-         if (Device::Allows(Backend::DEVICE_MASK))
-         {
-            const int ID = (dofs1D << 4) | quad1D;
-            switch (ID)
-            {
-               case 0x23:
-                  return internal::SmemPAHcurlMassApply3D<2,3>(
-                            dofs1D, quad1D, ne, symmetric,
-                            mapsO->B, mapsC->B, mapsO->Bt,
-                            mapsC->Bt, pa_data, x, y);
-               case 0x34:
-                  return internal::SmemPAHcurlMassApply3D<3,4>(
-                            dofs1D, quad1D, ne, symmetric,
-                            mapsO->B, mapsC->B, mapsO->Bt,
-                            mapsC->Bt, pa_data, x, y);
-               case 0x45:
-                  return internal::SmemPAHcurlMassApply3D<4,5>(
-                            dofs1D, quad1D, ne, symmetric,
-                            mapsO->B, mapsC->B, mapsO->Bt,
-                            mapsC->Bt, pa_data, x, y);
-               case 0x56:
-                  return internal::SmemPAHcurlMassApply3D<5,6>(
-                            dofs1D, quad1D, ne, symmetric,
-                            mapsO->B, mapsC->B, mapsO->Bt,
-                            mapsC->Bt, pa_data, x, y);
-               default:
-                  return internal::SmemPAHcurlMassApply3D(
-                            dofs1D, quad1D, ne, symmetric,
-                            mapsO->B, mapsC->B, mapsO->Bt,
-                            mapsC->Bt, pa_data, x, y);
-            }
-         }
-         else
-         {
-            internal::PAHcurlMassApply3D(dofs1D, quad1D, ne, symmetric, mapsO->B, mapsC->B,
-                                         mapsO->Bt, mapsC->Bt, pa_data, x, y);
-         }
-      }
-      else if (trial_div && test_div)
-      {
-         internal::PAHdivMassApply(3, dofs1D, quad1D, ne, symmetric, mapsO->B, mapsC->B,
-                                   mapsO->Bt, mapsC->Bt, pa_data, x, y);
-      }
-      else if (trial_curl && test_div)
-      {
-         const bool scalarCoeff = !(DQ || MQ);
-         internal::PAHcurlHdivMassApply3D(dofs1D, dofs1Dtest, quad1D, ne, scalarCoeff,
-                                          true, false, mapsO->B, mapsC->B, mapsOtest->Bt,
-                                          mapsCtest->Bt, pa_data, x, y);
-      }
-      else if (trial_div && test_curl)
-      {
-         const bool scalarCoeff = !(DQ || MQ);
-         internal::PAHcurlHdivMassApply3D(dofs1D, dofs1Dtest, quad1D, ne, scalarCoeff,
-                                          false, false, mapsO->B, mapsC->B, mapsOtest->Bt,
-                                          mapsCtest->Bt, pa_data, x, y);
-      }
-      else
-      {
-         MFEM_ABORT("Unknown kernel.");
-      }
-   }
-   else // 2D
-   {
-      if (trial_curl && test_curl)
-      {
-         internal::PAHcurlMassApply2D(dofs1D, quad1D, ne, symmetric, mapsO->B, mapsC->B,
-                                      mapsO->Bt, mapsC->Bt, pa_data, x, y);
-      }
-      else if (trial_div && test_div)
-      {
-         internal::PAHdivMassApply(2, dofs1D, quad1D, ne, symmetric, mapsO->B, mapsC->B,
-                                   mapsO->Bt,
-                                   mapsC->Bt, pa_data, x, y);
-      }
-      else if ((trial_curl && test_div) || (trial_div && test_curl))
-      {
-         const bool scalarCoeff = !(DQ || MQ);
-         internal::PAHcurlHdivMassApply2D(dofs1D, dofs1Dtest, quad1D, ne, scalarCoeff,
-                                          trial_curl, false, mapsO->B, mapsC->B,
-                                          mapsOtest->Bt, mapsCtest->Bt, pa_data, x, y);
-      }
-      else
-      {
-         MFEM_ABORT("Unknown kernel.");
-      }
-   }
+   const bool scalar_coeff = !(DQ || MQ);
+   ApplyPAKernels::Run(trial_fetype, test_fetype, dim, dofs1D, dofs1Dtest,
+                       quad1D, ne, symmetric, scalar_coeff, mapsO->B, mapsC->B,
+                       mapsOtest->Bt, mapsCtest->Bt, pa_data, x, y, dofs1D,
+                       dofs1Dtest, quad1D);
 }
 
 void VectorFEMassIntegrator::AddAbsMultPA(const Vector &x, Vector &y) const
 {
-   const bool trial_curl = (trial_fetype == mfem::FiniteElement::CURL);
-   const bool trial_div = (trial_fetype == mfem::FiniteElement::DIV);
-   const bool test_curl = (test_fetype == mfem::FiniteElement::CURL);
-   const bool test_div = (test_fetype == mfem::FiniteElement::DIV);
+   const bool scalar_coeff = !(DQ || MQ);
 
    Vector abs_pa_data(pa_data);
    abs_pa_data.Abs();
 
    Array<real_t> absBo(mapsO->B);
    Array<real_t> absBc(mapsC->B);
-   Array<real_t> absBto(mapsO->Bt);
-   Array<real_t> absBtc(mapsC->Bt);
    Array<real_t> absBto_t(mapsOtest->Bt);
    Array<real_t> absBtc_t(mapsCtest->Bt);
 
    absBo.Abs();
    absBc.Abs();
-   absBto.Abs();
-   absBtc.Abs();
    absBto_t.Abs();
    absBtc_t.Abs();
 
-   if (dim == 3)
-   {
-      if (trial_curl && test_curl)
-      {
-         if (Device::Allows(Backend::DEVICE_MASK))
-         {
-            const int ID = (dofs1D << 4) | quad1D;
-            switch (ID)
-            {
-               case 0x23:
-                  return internal::SmemPAHcurlMassApply3D<2,3>(
-                            dofs1D, quad1D, ne, symmetric,
-                            absBo, absBc, absBto, absBtc,
-                            abs_pa_data, x, y);
-               case 0x34:
-                  return internal::SmemPAHcurlMassApply3D<3,4>(
-                            dofs1D, quad1D, ne, symmetric,
-                            absBo, absBc, absBto, absBtc,
-                            abs_pa_data, x, y);
-               case 0x45:
-                  return internal::SmemPAHcurlMassApply3D<4,5>(
-                            dofs1D, quad1D, ne, symmetric,
-                            absBo, absBc, absBto, absBtc,
-                            abs_pa_data, x, y);
-               case 0x56:
-                  return internal::SmemPAHcurlMassApply3D<5,6>(
-                            dofs1D, quad1D, ne, symmetric,
-                            absBo, absBc, absBto, absBtc,
-                            abs_pa_data, x, y);
-               default:
-                  return internal::SmemPAHcurlMassApply3D(
-                            dofs1D, quad1D, ne, symmetric,
-                            absBo, absBc, absBto, absBtc,
-                            abs_pa_data, x, y);
-            }
-         }
-         else
-         {
-            internal::PAHcurlMassApply3D(dofs1D, quad1D, ne, symmetric,
-                                         absBo, absBc, absBto, absBtc,
-                                         abs_pa_data, x, y);
-         }
-      }
-      else if (trial_div && test_div)
-      {
-         internal::PAHdivMassApply(3, dofs1D, quad1D, ne, symmetric,
-                                   absBo, absBc, absBto, absBtc,
-                                   abs_pa_data, x, y);
-      }
-      else if (trial_curl && test_div)
-      {
-         const bool scalarCoeff = !(DQ || MQ);
-         internal::PAHcurlHdivMassApply3D(dofs1D, dofs1Dtest, quad1D, ne,
-                                          scalarCoeff, true, false,
-                                          absBo, absBc, absBto_t, absBtc_t,
-                                          abs_pa_data, x, y);
-      }
-      else if (trial_div && test_curl)
-      {
-         const bool scalarCoeff = !(DQ || MQ);
-         internal::PAHcurlHdivMassApply3D(dofs1D, dofs1Dtest, quad1D, ne,
-                                          scalarCoeff, false, false,
-                                          absBo, absBc, absBto_t, absBtc_t,
-                                          abs_pa_data, x, y);
-      }
-      else
-      {
-         MFEM_ABORT("Unknown kernel.");
-      }
-   }
-   else // 2D
-   {
-      if (trial_curl && test_curl)
-      {
-         internal::PAHcurlMassApply2D(dofs1D, quad1D, ne, symmetric,
-                                      absBo, absBc, absBto, absBtc,
-                                      abs_pa_data, x, y);
-      }
-      else if (trial_div && test_div)
-      {
-         internal::PAHdivMassApply(2, dofs1D, quad1D, ne, symmetric,
-                                   absBo, absBc, absBto, absBtc,
-                                   abs_pa_data, x, y);
-      }
-      else if ((trial_curl && test_div) || (trial_div && test_curl))
-      {
-         const bool scalarCoeff = !(DQ || MQ);
-         internal::PAHcurlHdivMassApply2D(dofs1D, dofs1Dtest, quad1D, ne,
-                                          scalarCoeff, trial_curl, false,
-                                          absBo, absBc, absBto_t, absBtc_t,
-                                          abs_pa_data, x, y);
-      }
-      else
-      {
-         MFEM_ABORT("Unknown kernel.");
-      }
-   }
+   ApplyPAKernels::Run(trial_fetype, test_fetype, dim, dofs1D, dofs1Dtest,
+                       quad1D, ne, symmetric, scalar_coeff, absBo, absBc,
+                       absBto_t, absBtc_t, abs_pa_data, x, y, dofs1D,
+                       dofs1Dtest, quad1D);
 }
 
 void VectorFEMassIntegrator::AddMultTransposePA(const Vector &x,

fem/intrules.cpp

@@ -542,7 +542,10 @@ void QuadratureFunctions1D::GaussJacobi(const int np, const real_t alpha,
          return;
    }
 
-#ifndef MFEM_USE_MPFR
+#ifdef MFEM_USE_MPFR
+   MFEM_WARNING("MPFR implementation of Gauss-Jacobi quadrature not implemented yet. Falling "
+                "back to double precision implementation...");
+#endif
 
    const int n = np;
    // common constants for Jacobi polynomials
@@ -611,13 +614,6 @@ void QuadratureFunctions1D::GaussJacobi(const int np, const real_t alpha,
                                                 ab + 1) / ((1.0 - xi*xi)*pp*pp) / pow(2, ab);
       // map nodes and weights to the interval [0,1]
    }
-
-#else // MFEM_USE_MPFR is defined
-
-   MFEM_ABORT("MPFR implementation of Gauss-Jacobi quadrature not defined yet");
-
-#endif // MFEM_USE_MPFR
-
 }
 
 

fem/qinterp/grad.hpp

@@ -268,8 +268,9 @@ static void Derivatives3D(const int NE,
       DeviceMatrix B(BG[0], D1D, Q1D);
       DeviceMatrix G(BG[1], D1D, Q1D);
 
-      MFEM_SHARED real_t sm0[3][MQ1*MQ1*MQ1];
-      MFEM_SHARED real_t sm1[3][MQ1*MQ1*MQ1];
+      constexpr int MDQ = MD1 > MQ1 ? MD1 : MQ1;
+      MFEM_SHARED real_t sm0[3][MD1*MD1*MDQ];
+      MFEM_SHARED real_t sm1[3][MD1*MQ1*MQ1];
       DeviceTensor<3> X(sm0[2], D1D, D1D, D1D);
       DeviceTensor<3> DDQ0(sm0[0], D1D, D1D, Q1D);
       DeviceTensor<3> DDQ1(sm0[1], D1D, D1D, Q1D);

general/backends.hpp

@@ -66,6 +66,23 @@ constexpr bool mfem_use_gpu = false;
 #define MFEM_THREAD_SIZE(k) 1
 #define MFEM_FOREACH_THREAD(i,k,N) for(int i=0; i<N; i++)
 #define MFEM_FOREACH_THREAD_DIRECT(i,k,N) MFEM_FOREACH_THREAD(i,k,N)
+// Assigns a thread block shaped (SX,SY,SZ) contiguous in x.
+// Example (3,2,1) block:
+// 0 (0,0), 1 (1,0), 2 (2,0)
+// 3 (1,0), 4 (1,1), 5 (2,1)
+#define MFEM_FOREACH_THREAD_DIRECT_3D(ix, iy, iz, k, SX, SY, SZ)               \
+   for (int iz = 0; iz < SZ; ++iz)                                             \
+      for (int iy = 0; iy < SY; ++iy)                                          \
+         for (int ix = 0; ix < SX; ++ix)
+// Assigns a thread block shaped (OX,OY,OZ) to work on items (SX,SY,SZ),
+// contiguous in x. This intentionally offsets threads within the block to avoid
+// shared memory bank conflicts.
+// Example (3,2,1) block assigned to work on (2,2,1) items:
+// 0 (0,0), 1 (1,0), 2 (N/A)
+// 3 (1,0), 4 (1,1), 5 (N/A)
+#define MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(ix, iy, iz, k, SX, SY, SZ, OX,    \
+                                             OY, OZ)                           \
+   MFEM_FOREACH_THREAD_DIRECT_3D(ix, iy, iz, k, SX, SY, SZ)
 #endif
 
 // 'double' and 'float' atomicAdd implementation for previous versions of CUDA

general/cuda.hpp

@@ -49,6 +49,23 @@ constexpr bool mfem_use_gpu = true;
 #define MFEM_THREAD_SIZE(k) blockDim.k
 #define MFEM_FOREACH_THREAD(i,k,N) for(int i=threadIdx.k; i<N; i+=blockDim.k)
 #define MFEM_FOREACH_THREAD_DIRECT(i,k,N) if(const int i=threadIdx.k; i<N)
+// Assigns a thread block shaped (SX,SY,SZ) contiguous in x.
+// Example (3,2,1) block:
+// 0 (0,0), 1 (1,0), 2 (2,0)
+// 3 (1,0), 4 (1,1), 5 (2,1)
+#define MFEM_FOREACH_THREAD_DIRECT_3D(ix, iy, iz, k, SX, SY, SZ)               \
+   if (int ix = threadIdx.k % (SX), iy = threadIdx.k / (SX), iz = iy / (SY);   \
+       (iy %= (SY)), (threadIdx.k < (SX) * (SY) * (SZ)))
+// Assigns a thread block shaped (OX,OY,OZ) to work on items (SX,SY,SZ),
+// contiguous in x. This intentionally offsets threads within the block to avoid
+// shared memory bank conflicts.
+// Example (3,2,1) block assigned to work on (2,2,1) items:
+// 0 (0,0), 1 (1,0), 2 (N/A)
+// 3 (1,0), 4 (1,1), 5 (N/A)
+#define MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(ix, iy, iz, k, SX, SY, SZ, OX,    \
+                                             OY, OZ)                           \
+   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__)
 

general/hip.hpp

@@ -51,6 +51,25 @@ constexpr bool mfem_use_gpu = true;
    for(int i=hipThreadIdx_ ##k; i<N; i+=hipBlockDim_ ##k)
 #define MFEM_FOREACH_THREAD_DIRECT(i,k,N) \
    if(const int i=hipThreadIdx_ ##k; i<N)
+// Assigns a thread block shaped (SX,SY,SZ) contiguous in x.
+// Example (3,2,1) block:
+// 0 (0,0), 1 (1,0), 2 (2,0)
+// 3 (1,0), 4 (1,1), 5 (2,1)
+#define MFEM_FOREACH_THREAD_DIRECT_3D(ix, iy, iz, k, SX, SY, SZ)               \
+   if (int ix = hipThreadIdx_##k % (SX), iy = hipThreadIdx_##k / (SX),         \
+       iz = iy / (SY);                                                         \
+       (iy %= (SY)), (hipThreadIdx_##k < (SX) * (SY) * (SZ)))
+// Assigns a thread block shaped (OX,OY,OZ) to work on items (SX,SY,SZ),
+// contiguous in x. This intentionally offsets threads within the block to avoid
+// shared memory bank conflicts.
+// Example (3,2,1) block assigned to work on (2,2,1) items:
+// 0 (0,0), 1 (1,0), 2 (N/A)
+// 3 (1,0), 4 (1,1), 5 (N/A)
+#define MFEM_FOREACH_THREAD_DIRECT_3D_OFFSET(ix, iy, iz, k, SX, SY, SZ, OX,    \
+                                             OY, OZ)                           \
+   if (int ix = hipThreadIdx_##k % (OX), iy = hipThreadIdx_##k / (OX),         \
+       iz = iy / (OY);                                                         \
+       (ix < (SX)) && ((iy %= (OY)) < (SY)) && (iz < (SZ)))
 #endif // defined(__HIP_DEVICE_COMPILE__)
 #endif // defined(MFEM_USE_HIP) && defined(__HIP__)
 

linalg/batched/batched.cpp

@@ -13,6 +13,7 @@
 #include "native.hpp"
 #include "gpu_blas.hpp"
 #include "magma.hpp"
+#include "../../general/reducers.hpp"
 
 namespace mfem
 {
@@ -119,4 +120,16 @@ void BatchedLinAlgBase::MultTranspose(const DenseTensor &A, const Vector &x,
    AddMult(A, x, y, 1.0, 0.0, Op::T);
 }
 
+void VerifyBatchedLUInfo(const Array<int> &info_array, const char *message)
+{
+   static Array<int> workspace;
+   int status = 0;
+   const int *d_info = info_array.Read();
+   mfem::reduce(
+      info_array.Size(), status,
+   [=] MFEM_HOST_DEVICE (int i, int &r) { r |= d_info[i]; },
+   BOrReducer<int> {}, true, workspace);
+   MFEM_VERIFY(status == 0, message);
+}
+
 }

linalg/batched/batched.hpp

@@ -141,6 +141,9 @@ public:
    virtual ~BatchedLinAlgBase() { }
 };
 
+/// Check that all batched LU info values are zero.
+void VerifyBatchedLUInfo(const Array<int> &info_array, const char *message);
+
 } // namespace mfem
 
 #endif

linalg/batched/gpu_blas.cpp

@@ -126,7 +126,8 @@ void GPUBlasBatchedLinAlg::LUFactor(DenseTensor &A, Array<int> &P) const
    const blasStatus_t status = MFEM_GPUBLAS_PREFIX(getrfBatched)(
                                   GPUBlas::Handle(), n, d_A_ptrs, n, P.Write(),
                                   info_array.Write(), n_mat);
-   MFEM_VERIFY(status == MFEM_BLAS_SUCCESS, "");
+   MFEM_VERIFY(status == MFEM_BLAS_SUCCESS, "GPU BLAS error.");
+   VerifyBatchedLUInfo(info_array, "Batch LU factorization failed");
 }
 
 void GPUBlasBatchedLinAlg::LUSolve(
@@ -189,12 +190,14 @@ void GPUBlasBatchedLinAlg::Invert(DenseTensor &A) const
    status = MFEM_GPUBLAS_PREFIX(getrfBatched)(
                GPUBlas::Handle(), n, d_LU_ptrs, n, P.Write(),
                info_array.Write(), n_mat);
-   MFEM_VERIFY(status == MFEM_BLAS_SUCCESS, "");
+   MFEM_VERIFY(status == MFEM_BLAS_SUCCESS, "GPU BLAS error.");
+   VerifyBatchedLUInfo(info_array, "Batch LU factorization failed");
 
    status = MFEM_GPUBLAS_PREFIX(getriBatched)(
                GPUBlas::Handle(), n, d_LU_ptrs, n, P.ReadWrite(), d_A_ptrs, n,
                info_array.Write(), n_mat);
-   MFEM_VERIFY(status == MFEM_BLAS_SUCCESS, "");
+   MFEM_VERIFY(status == MFEM_BLAS_SUCCESS, "GPU BLAS error.");
+   VerifyBatchedLUInfo(info_array, "Batch matrix inversion failed");
 }
 
 #endif

linalg/batched/magma.cpp

@@ -99,7 +99,8 @@ void MagmaBatchedLinAlg::LUFactor(DenseTensor &A, Array<int> &P) const
    const magma_int_t status = MFEM_MAGMA_PREFIX(getrf_batched)(
                                  n, n, d_A_ptrs, n, d_P_ptrs,
                                  info_array.Write(), n_mat, Magma::Queue());
-   MFEM_VERIFY(status == MAGMA_SUCCESS, "");
+   MFEM_VERIFY(status == MAGMA_SUCCESS, "MAGMA error.");
+   VerifyBatchedLUInfo(info_array, "Batch LU factorization failed");
 }
 
 void MagmaBatchedLinAlg::LUSolve(
@@ -169,12 +170,14 @@ void MagmaBatchedLinAlg::Invert(DenseTensor &A) const
    status = MFEM_MAGMA_PREFIX(getrf_batched)(
                n, n, d_LU_ptrs, n, d_P_ptrs, info_array.Write(), n_mat,
                Magma::Queue());
-   MFEM_VERIFY(status == MAGMA_SUCCESS, "");
+   MFEM_VERIFY(status == MAGMA_SUCCESS, "MAGMA error.");
+   VerifyBatchedLUInfo(info_array, "Batch LU factorization failed");
 
    status = MFEM_MAGMA_PREFIX(getri_outofplace_batched)(
                n, d_LU_ptrs, n, d_P_ptrs, d_A_ptrs, n, info_array.Write(),
                n_mat, Magma::Queue());
-   MFEM_VERIFY(status == MAGMA_SUCCESS, "");
+   MFEM_VERIFY(status == MAGMA_SUCCESS, "MAGMA error.");
+   VerifyBatchedLUInfo(info_array, "Batch matrix inversion failed");
 }
 
 } // namespace mfem

linalg/complex_operator.cpp

@@ -246,6 +246,10 @@ SparseMatrix * ComplexSparseMatrix::GetSystemMatrix() const
    const int  nrows_i = (A_i)?A_i->Height():0;
    const int    nrows = std::max(nrows_r, nrows_i);
 
+   const int  ncols_r = (A_r)?A_r->Width():0;
+   const int  ncols_i = (A_i)?A_i->Width():0;
+   const int    ncols = std::max(ncols_r, ncols_i);
+
    const int     *I_r = (A_r)?A_r->GetI():NULL;
    const int     *I_i = (A_i)?A_i->GetI():NULL;
 
@@ -280,7 +284,7 @@ SparseMatrix * ComplexSparseMatrix::GetSystemMatrix() const
             J[I[i] + j] = J_r[I_r[i] + j];
             D[I[i] + j] = D_r[I_r[i] + j];
 
-            J[I[i+nrows] + off_i + j] = J_r[I_r[i] + j] + nrows;
+            J[I[i+nrows] + off_i + j] = J_r[I_r[i] + j] + ncols;
             D[I[i+nrows] + off_i + j] = factor*D_r[I_r[i] + j];
          }
       }
@@ -289,7 +293,7 @@ SparseMatrix * ComplexSparseMatrix::GetSystemMatrix() const
          const int off_r = (I_r)?(I_r[i+1] - I_r[i]):0;
          for (int j=0; j<I_i[i+1] - I_i[i]; j++)
          {
-            J[I[i] + off_r + j] =  J_i[I_i[i] + j] + nrows;
+            J[I[i] + off_r + j] =  J_i[I_i[i] + j] + ncols;
             D[I[i] + off_r + j] = -D_i[I_i[i] + j];
 
             J[I[i+nrows] + j] = J_i[I_i[i] + j];
@@ -892,12 +896,12 @@ ComplexHypreParMatrix::getColStartStop(const HypreParMatrix * A_r,
    HYPRE_BigInt loc_start_stop[2];
    offd_col_start_stop = new HYPRE_BigInt[2 * num_recv_procs];
 
-   const HYPRE_BigInt * row_part = (A_r) ? A_r->RowPart() :
-                                   ((A_i) ? A_i->RowPart() : NULL);
+   const HYPRE_BigInt * col_part = (A_r) ? A_r->ColPart() :
+                                   ((A_i) ? A_i->ColPart() : NULL);
 
-   int row_part_ind = (HYPRE_AssumedPartitionCheck()) ? 0 : myid_;
-   loc_start_stop[0] = row_part[row_part_ind];
-   loc_start_stop[1] = row_part[row_part_ind+1];
+   int col_part_ind = (HYPRE_AssumedPartitionCheck()) ? 0 : myid_;
+   loc_start_stop[0] = col_part[col_part_ind];
+   loc_start_stop[1] = col_part[col_part_ind+1];
 
    MPI_Request * req = new MPI_Request[send_procs.size()+recv_procs.size()];
    MPI_Status * stat = new MPI_Status[send_procs.size()+recv_procs.size()];

makefile

@@ -810,6 +810,7 @@ MINIAPPS_SUBDIRS = dpg/util hooke/operators hooke/preconditioners \
    hooke/materials hooke/kernels
 FORMAT_FILES += $(foreach dir,$(TESTS_SUBDIRS),tests/$(dir)/*.?pp)
 FORMAT_FILES += $(foreach dir,$(UNIT_TESTS_SUBDIRS),tests/unit/$(dir)/*.?pp)
+FORMAT_FILES += tests/unit/fem/specializations/*.?pp
 FORMAT_FILES += $(foreach dir,$(MINIAPPS_SUBDIRS),miniapps/$(dir)/*.?pp)
 FORMAT_FILES += config/cmake/config.hpp.in config/config.hpp.in mfem*.hpp
 FORMAT_EXCLUDE = general/tinyxml2.cpp tests/unit/catch.hpp

mesh/mesh.cpp

@@ -667,9 +667,84 @@ void Mesh::GetEdgeTransformation(int EdgeNo,
          }
          EdTr->SetFE(edge_el);
       }
-      else
+      else // L2 Nodes (e.g., periodic mesh), go through the face containing the edge
       {
-         MFEM_ABORT("Not implemented.");
+         // Search for a face that contains this edge
+         GetEdgeFaceTable();
+
+         Array<int> faces_e;
+         edge_face->GetRow(EdgeNo, faces_e);
+
+         MFEM_VERIFY(faces_e.Size() > 0, "Edge not found in any face!");
+         const int face_no = faces_e[0];
+
+         // Get edge local index and orientation
+         Array<int> edges_f, oris_f;
+         GetFaceEdges(face_no, edges_f, oris_f);
+         const int local_idx = edges_f.Find(EdgeNo);
+         MFEM_ASSERT(local_idx >= 0, "Edge not found on the face!");
+         const int edge_ori = oris_f[local_idx] > 0 ? 0 : 1;
+
+         // Get face information
+         const FaceInfo &face_info = faces_info[face_no];
+
+         // Get transformation from face to edge
+         IntegrationPointTransformation LocEdge;
+         int edge_info = EncodeFaceInfo(local_idx, edge_ori);
+         Element::Type face_type = GetFaceElementType(face_no);
+
+         switch (face_type)
+         {
+            case Element::TRIANGLE:
+               GetLocalSegToTriTransformation(LocEdge.Transf, edge_info);
+               break;
+            case Element::QUADRILATERAL:
+               GetLocalSegToQuadTransformation(LocEdge.Transf, edge_info);
+               break;
+            default:
+               MFEM_ABORT("Unsupported face type for edge transformation!");
+         }
+
+         // Get edge element
+         const int order = Nodes->FESpace()->GetElementOrder(face_info.Elem1No);
+         const L2_FECollection *l2_fec = dynamic_cast<const L2_FECollection*>
+                                         (Nodes->FESpace()->FEColl());
+         if (l2_fec)
+         {
+            // L2 elements do not have a defined trace space
+            if (!EdgeTransfElement || EdgeTransfElement->GetOrder() != order
+                || EdgeTransfElement->GetBasisType() != l2_fec->GetBasisType())
+            {
+               EdgeTransfElement = make_unique<L2_SegmentElement>(
+                                      order, l2_fec->GetBasisType());
+            }
+            edge_el = EdgeTransfElement.get();
+         }
+         else
+         {
+            MFEM_ABORT("Unsupported finite element collection.");
+         }
+
+         // Map edge nodes to face reference space
+         IntegrationRule face_ir(edge_el->GetDof());
+         LocEdge.Transform(edge_el->GetNodes(), face_ir);
+
+         // Then, map from face to element
+         IntegrationPointTransformation Loc1;
+         GetLocalFaceTransformation(face_type,
+                                    GetElementType(face_info.Elem1No),
+                                    Loc1.Transf, face_info.Elem1Inf);
+
+         IntegrationRule elem_ir(edge_el->GetDof());
+         Loc1.Transf.ElementNo = face_info.Elem1No;
+         Loc1.Transf.ElementType = ElementTransformation::ELEMENT;
+         Loc1.Transf.mesh = this;
+         Loc1.Transform(face_ir, elem_ir);
+
+         // Finally, get the physical coordinates
+         Nodes->GetVectorValues(Loc1.Transf, elem_ir, pm);
+
+         EdTr->SetFE(edge_el);
       }
    }
 }
@@ -1824,8 +1899,8 @@ void Mesh::Init()
 
 void Mesh::InitTables()
 {
-   el_to_edge =
-      el_to_face = el_to_el = bel_to_edge = face_edge = edge_vertex = NULL;
+   el_to_edge = el_to_face = el_to_el = bel_to_edge = NULL;
+   face_edge = edge_face = edge_vertex = NULL;
    face_to_elem = NULL;
 }
 
@@ -1848,6 +1923,7 @@ void Mesh::DestroyTables()
    }
 
    delete face_edge;
+   delete edge_face;
    delete edge_vertex;
 
    delete face_to_elem;
@@ -1921,6 +1997,7 @@ void Mesh::ResetLazyData()
 {
    delete el_to_el;     el_to_el = NULL;
    delete face_edge;    face_edge = NULL;
+   delete edge_face;    edge_face = NULL;
    delete face_to_elem;    face_to_elem = NULL;
    delete edge_vertex;  edge_vertex = NULL;
    DeleteGeometricFactors();
@@ -2845,6 +2922,7 @@ void Mesh::ReorderElements(const Array<int> &ordering, bool reorder_vertices)
    //                 boundary element ordering
    // - el_to_el    - no need to rebuild
    // - face_edge   - no need to rebuild
+   // - edge_face   - no need to rebuild
    // - edge_vertex - no need to rebuild
    // - geom_factors - no need to rebuild
 
@@ -4598,8 +4676,9 @@ Mesh::Mesh(const Mesh &mesh, bool copy_nodes)
    // Do NOT copy the element-to-element Table, el_to_el
    el_to_el = NULL;
 
-   // Do NOT copy the face-to-edge Table, face_edge
+   // Do NOT copy the face-to-edge Table, face_edge and edge_face
    face_edge = NULL;
+   edge_face = NULL;
    face_to_elem = NULL;
 
    // Copy the edge-to-vertex Table, edge_vertex
@@ -8094,6 +8173,22 @@ Table *Mesh::GetFaceEdgeTable() const
    return (face_edge);
 }
 
+Table *Mesh::GetEdgeFaceTable() const
+{
+   if (edge_face)
+   {
+      return edge_face;
+   }
+
+   if (Dim != 3)
+   {
+      return NULL;
+   }
+
+   edge_face = Transpose(*GetFaceEdgeTable());
+   return edge_face;
+}
+
 Table *Mesh::GetEdgeVertexTable() const
 {
    if (edge_vertex)
@@ -11452,6 +11547,7 @@ void Mesh::Swap(Mesh& other, bool non_geometry)
    mfem::Swap(bel_to_edge, other.bel_to_edge);
    mfem::Swap(be_to_face, other.be_to_face);
    mfem::Swap(face_edge, other.face_edge);
+   mfem::Swap(edge_face, other.edge_face);
    mfem::Swap(face_to_elem, other.face_to_elem);
    mfem::Swap(edge_vertex, other.edge_vertex);
 

mesh/mesh.hpp

@@ -250,16 +250,18 @@ protected:
    Table *bel_to_edge;    // for 3D only
 
    // Note that the following tables are owned by this class and should not be
-   // deleted by the caller. Of these three tables, only face_edge and
+   // deleted by the caller. Of these four tables, only face_edge, edge_face and
    // edge_vertex are returned by access functions.
    mutable Table *face_to_elem;  // Used by FindFaceNeighbors, not returned.
    mutable Table *face_edge;     // Returned by GetFaceEdgeTable().
+   mutable Table *edge_face;     // Returned by GetEdgeFaceTable().
    mutable Table *edge_vertex;   // Returned by GetEdgeVertexTable().
 
    IsoparametricTransformation Transformation, Transformation2;
    IsoparametricTransformation BdrTransformation;
    IsoparametricTransformation FaceTransformation, EdgeTransformation;
    FaceElementTransformations FaceElemTr;
+   mutable std::unique_ptr<L2_SegmentElement> EdgeTransfElement;
 
    // refinement embeddings for forward compatibility with NCMesh
    mutable CoarseFineTransformations CoarseFineTr;
@@ -1731,6 +1733,11 @@ public:
    /// @note The returned object should NOT be deleted by the caller.
    Table *GetFaceEdgeTable() const;
 
+   /// Returns the edge-to-face Table (3D)
+   ///
+   /// @note The returned object should NOT be deleted by the caller.
+   Table *GetEdgeFaceTable() const;
+
    /// Returns the edge-to-vertex Table (3D)
    ///
    /// @note The returned object should NOT be deleted by the caller.

mesh/pmesh.cpp

@@ -4866,6 +4866,13 @@ void ParMesh::Print(std::ostream &os, const std::string &comments) const
       return;
    }
 
+   if (pncmesh && pncmesh->using_scaling)
+   {
+      // For nodes scaling, we write the file in the format MFEM NC mesh v1.1.
+      Printer(os, "", comments);
+      return;
+   }
+
    const Array<int>* s2l_face;
    if (!pncmesh)
    {

mesh/pncmesh.cpp

@@ -28,6 +28,48 @@ namespace mfem
 
 using namespace bin_io;
 
+static int GetHexEdgeSplit(const int* nodes, int v1, int v2);
+
+static bool SameSplitScale(real_t a, real_t b)
+{
+#ifdef MFEM_USE_DOUBLE
+   constexpr real_t rel_tol = 1.0e-8;
+#else
+   constexpr real_t rel_tol = 1.0e-5;
+#endif
+   return std::abs(a - b) <= rel_tol *
+          std::max(real_t(1.0), std::max(std::abs(a), std::abs(b)));
+}
+
+static real_t DirectedHexEdgeScale(const int* nodes, const Refinement &ref,
+                                   int v0, int v1)
+{
+   const int dir = GetHexEdgeSplit(nodes, v0, v1);
+   static const int split_edges[3][4][2] =
+   {
+      {{0, 1}, {3, 2}, {4, 5}, {7, 6}},
+      {{1, 2}, {0, 3}, {5, 6}, {4, 7}},
+      {{0, 4}, {1, 5}, {2, 6}, {3, 7}}
+   };
+
+   for (int i = 0; i < 4; i++)
+   {
+      const int a = nodes[split_edges[dir][i][0]];
+      const int b = nodes[split_edges[dir][i][1]];
+      if (a == v0 && b == v1)
+      {
+         return ref.s[dir];
+      }
+      if (a == v1 && b == v0)
+      {
+         return 1.0 - ref.s[dir];
+      }
+   }
+
+   MFEM_ABORT("Shared face edge does not match the refinement direction.");
+   return 0.0;
+}
+
 ParNCMesh::ParNCMesh(MPI_Comm comm, const NCMesh &ncmesh,
                      const int *partitioning)
    : NCMesh(ncmesh)
@@ -1555,7 +1597,7 @@ bool ParNCMesh::AnisotropicConflict(const Array<Refinement> &refinements,
       ElementNeighborProcessors(elem, ranks);
       for (int j = 0; j < ranks.Size(); j++)
       {
-         send_ref[ranks[j]].AddRefinement(elem, ref.GetType());
+         send_ref[ranks[j]].AddRefinement(elem, ref);
       }
    }
 
@@ -1576,8 +1618,8 @@ bool ParNCMesh::AnisotropicConflict(const Array<Refinement> &refinements,
    for (int i = 0; i < refinements.Size(); i++)
    {
       const Refinement &ref = refinements[i];
-      CheckRefinement(leaf_elements[ref.index], ref.GetType(), refinements,
-                      elemToRef, conflicts);
+      CheckRefinement(leaf_elements[ref.index], ref, refinements, elemToRef,
+                      conflicts);
    }
 
    // Receive (ghost layer) refinements from all neighbors
@@ -1593,7 +1635,9 @@ bool ParNCMesh::AnisotropicConflict(const Array<Refinement> &refinements,
       // check the ghost refinements
       for (int i = 0; i < msg.Size(); i++)
       {
-         CheckRefinement(msg.elements[i], msg.values[i], refinements, elemToRef,
+         Refinement ghost_ref(msg.elements[i], msg.values[i].ref_type);
+         ghost_ref.SetScaleForType(msg.values[i].scale);
+         CheckRefinement(msg.elements[i], ghost_ref, refinements, elemToRef,
                          conflicts);
       }
    }
@@ -1749,7 +1793,7 @@ int FindHexFace(const int* no, int vn1, int vn2, int vn3, int vn4)
 
 // Assumption: v1 and v2 are indices of hex vertices connected by an edge.
 // The return value is {0,1,2} denoting split {X,Y,Z}.
-int GetHexEdgeSplit(const int* nodes, int v1, int v2)
+static int GetHexEdgeSplit(const int* nodes, int v1, int v2)
 {
    Array<int> v(2);
    v[0] = v1;
@@ -1780,7 +1824,8 @@ int GetHexEdgeSplit(const int* nodes, int v1, int v2)
    return edgeDir[edge];
 }
 
-void ParNCMesh::CheckRefAnisoFace(int elem, int vn1, int vn2, int vn3, int vn4,
+void ParNCMesh::CheckRefAnisoFace(const Refinement &ref, int elem,
+                                  int vn1, int vn2, int vn3, int vn4,
                                   const Array<Refinement> &refinements,
                                   const std::map<int, int> &elemToRef,
                                   std::set<int> &conflicts)
@@ -1798,11 +1843,11 @@ void ParNCMesh::CheckRefAnisoFace(int elem, int vn1, int vn2, int vn3, int vn4,
    if (elemToRef.count(nghbIndex) > 0)
    {
       const int refIndex = elemToRef.at(nghbIndex);
-      const Refinement& ref = refinements[refIndex];
+      const Refinement& nghb_ref = refinements[refIndex];
 
       bool refDir[3];
       for (int i=0; i<3; ++i)
-         refDir[i] = ref.s[i] > real_t{0};
+         refDir[i] = nghb_ref.s[i] > real_t{0};
 
       const int localFace = FindHexFace(nghb.node, vn1, vn2, vn3, vn4);
       const int faceDir = GetHexFaceDir(localFace);
@@ -1834,30 +1879,50 @@ void ParNCMesh::CheckRefAnisoFace(int elem, int vn1, int vn2, int vn3, int vn4,
          MFEM_ASSERT(cnt == 2 && hexSplitOnFace >= 0, "");
 
          const int edgeSplit = GetHexEdgeSplit(nghb.node, vn1, vn2);
-         if (edgeSplit != hexSplitOnFace) { conflicts.insert(refIndex); }
+         if (edgeSplit != hexSplitOnFace)
+         {
+            conflicts.insert(refIndex);
+         }
+         else
+         {
+            const real_t elem_scale =
+               DirectedHexEdgeScale(elements[elem].node, ref, vn1, vn2);
+            const real_t nghb_scale =
+               DirectedHexEdgeScale(nghb.node, nghb_ref, vn1, vn2);
+            if (!SameSplitScale(elem_scale, nghb_scale))
+            {
+               conflicts.insert(refIndex);
+            }
+         }
       }
    }
    // The else case is that the neighbor is not refined, so there is no need to
    // check for conflicts.
 }
 
-void ParNCMesh::CheckRefIsoFace(int elem, int vn1, int vn2, int vn3, int vn4,
+void ParNCMesh::CheckRefIsoFace(const Refinement &ref, int elem,
+                                int vn1, int vn2, int vn3, int vn4,
                                 int en1, int en2, int en3, int en4,
                                 const Array<Refinement> &refinements,
                                 const std::map<int, int> &elemToRef,
                                 std::set<int> &conflicts)
 {
-   CheckRefAnisoFace(elem, vn1, vn2, en2, en4, refinements, elemToRef, conflicts);
-   CheckRefAnisoFace(elem, en4, en2, vn3, vn4, refinements, elemToRef, conflicts);
-   CheckRefAnisoFace(elem, vn4, vn1, en1, en3, refinements, elemToRef, conflicts);
-   CheckRefAnisoFace(elem, en3, en1, vn2, vn3, refinements, elemToRef, conflicts);
+   CheckRefAnisoFace(ref, elem, vn1, vn2, en2, en4, refinements, elemToRef,
+                     conflicts);
+   CheckRefAnisoFace(ref, elem, en4, en2, vn3, vn4, refinements, elemToRef,
+                     conflicts);
+   CheckRefAnisoFace(ref, elem, vn4, vn1, en1, en3, refinements, elemToRef,
+                     conflicts);
+   CheckRefAnisoFace(ref, elem, en3, en1, vn2, vn3, refinements, elemToRef,
+                     conflicts);
 }
 
-void ParNCMesh::CheckRefinement(int elem, char ref_type,
+void ParNCMesh::CheckRefinement(int elem, const Refinement &ref,
                                 const Array<Refinement> &refinements,
                                 const std::map<int, int> &elemToRef,
                                 std::set<int> &conflicts)
 {
+   const char ref_type = ref.GetType();
    const Element &el = elements[elem];
    MFEM_ASSERT(el.geom == Geometry::CUBE && el.ref_type == 0,
                "Element must be an unrefined hexahedron");
@@ -1868,46 +1933,46 @@ void ParNCMesh::CheckRefinement(int elem, char ref_type,
    // This follows the logic of NCMesh::RefineElement().
    if (ref_type == Refinement::X) // split along X axis
    {
-      CheckRefAnisoFace(elem, no[0], no[1], no[5], no[4], refinements,
+      CheckRefAnisoFace(ref, elem, no[0], no[1], no[5], no[4], refinements,
                         elemToRef, conflicts);
-      CheckRefAnisoFace(elem, no[2], no[3], no[7], no[6], refinements,
+      CheckRefAnisoFace(ref, elem, no[2], no[3], no[7], no[6], refinements,
                         elemToRef, conflicts);
-      CheckRefAnisoFace(elem, no[4], no[5], no[6], no[7], refinements,
+      CheckRefAnisoFace(ref, elem, no[4], no[5], no[6], no[7], refinements,
                         elemToRef, conflicts);
-      CheckRefAnisoFace(elem, no[3], no[2], no[1], no[0], refinements,
+      CheckRefAnisoFace(ref, elem, no[3], no[2], no[1], no[0], refinements,
                         elemToRef, conflicts);
    }
    else if (ref_type == Refinement::Y) // split along Y axis
    {
-      CheckRefAnisoFace(elem, no[1], no[2], no[6], no[5], refinements,
+      CheckRefAnisoFace(ref, elem, no[1], no[2], no[6], no[5], refinements,
                         elemToRef, conflicts);
-      CheckRefAnisoFace(elem, no[3], no[0], no[4], no[7], refinements,
+      CheckRefAnisoFace(ref, elem, no[3], no[0], no[4], no[7], refinements,
                         elemToRef, conflicts);
-      CheckRefAnisoFace(elem, no[5], no[6], no[7], no[4], refinements,
+      CheckRefAnisoFace(ref, elem, no[5], no[6], no[7], no[4], refinements,
                         elemToRef, conflicts);
-      CheckRefAnisoFace(elem, no[0], no[3], no[2], no[1], refinements,
+      CheckRefAnisoFace(ref, elem, no[0], no[3], no[2], no[1], refinements,
                         elemToRef, conflicts);
    }
    else if (ref_type == Refinement::Z) // split along Z axis
    {
-      CheckRefAnisoFace(elem, no[4], no[0], no[1], no[5], refinements,
+      CheckRefAnisoFace(ref, elem, no[4], no[0], no[1], no[5], refinements,
                         elemToRef, conflicts);
-      CheckRefAnisoFace(elem, no[5], no[1], no[2], no[6], refinements,
+      CheckRefAnisoFace(ref, elem, no[5], no[1], no[2], no[6], refinements,
                         elemToRef, conflicts);
-      CheckRefAnisoFace(elem, no[6], no[2], no[3], no[7], refinements,
+      CheckRefAnisoFace(ref, elem, no[6], no[2], no[3], no[7], refinements,
                         elemToRef, conflicts);
-      CheckRefAnisoFace(elem, no[7], no[3], no[0], no[4], refinements,
+      CheckRefAnisoFace(ref, elem, no[7], no[3], no[0], no[4], refinements,
                         elemToRef, conflicts);
    }
    else if (ref_type == Refinement::XY) // XY split
    {
-      CheckRefAnisoFace(elem, no[0], no[1], no[5], no[4], refinements,
+      CheckRefAnisoFace(ref, elem, no[0], no[1], no[5], no[4], refinements,
                         elemToRef, conflicts);
-      CheckRefAnisoFace(elem, no[1], no[2], no[6], no[5], refinements,
+      CheckRefAnisoFace(ref, elem, no[1], no[2], no[6], no[5], refinements,
                         elemToRef, conflicts);
-      CheckRefAnisoFace(elem, no[2], no[3], no[7], no[6], refinements,
+      CheckRefAnisoFace(ref, elem, no[2], no[3], no[7], no[6], refinements,
                         elemToRef, conflicts);
-      CheckRefAnisoFace(elem, no[3], no[0], no[4], no[7], refinements,
+      CheckRefAnisoFace(ref, elem, no[3], no[0], no[4], no[7], refinements,
                         elemToRef, conflicts);
 
       const int mid01 = GetMidEdgeNode(no[0], no[1]);
@@ -1920,20 +1985,20 @@ void ParNCMesh::CheckRefinement(int elem, char ref_type,
       const int mid67 = GetMidEdgeNode(no[6], no[7]);
       const int mid74 = GetMidEdgeNode(no[7], no[4]);
 
-      CheckRefIsoFace(elem, no[3], no[2], no[1], no[0], mid23, mid12, mid01,
+      CheckRefIsoFace(ref, elem, no[3], no[2], no[1], no[0], mid23, mid12, mid01,
                       mid30, refinements, elemToRef, conflicts);
-      CheckRefIsoFace(elem, no[4], no[5], no[6], no[7], mid45, mid56, mid67,
+      CheckRefIsoFace(ref, elem, no[4], no[5], no[6], no[7], mid45, mid56, mid67,
                       mid74, refinements, elemToRef, conflicts);
    }
    else if (ref_type == Refinement::XZ) // XZ split
    {
-      CheckRefAnisoFace(elem, no[3], no[2], no[1], no[0], refinements,
+      CheckRefAnisoFace(ref, elem, no[3], no[2], no[1], no[0], refinements,
                         elemToRef, conflicts);
-      CheckRefAnisoFace(elem, no[2], no[6], no[5], no[1], refinements,
+      CheckRefAnisoFace(ref, elem, no[2], no[6], no[5], no[1], refinements,
                         elemToRef, conflicts);
-      CheckRefAnisoFace(elem, no[6], no[7], no[4], no[5], refinements,
+      CheckRefAnisoFace(ref, elem, no[6], no[7], no[4], no[5], refinements,
                         elemToRef, conflicts);
-      CheckRefAnisoFace(elem, no[7], no[3], no[0], no[4], refinements,
+      CheckRefAnisoFace(ref, elem, no[7], no[3], no[0], no[4], refinements,
                         elemToRef, conflicts);
 
       const int mid01 = GetMidEdgeNode(no[0], no[1]);
@@ -1946,9 +2011,9 @@ void ParNCMesh::CheckRefinement(int elem, char ref_type,
       const int mid26 = GetMidEdgeNode(no[2], no[6]);
       const int mid37 = GetMidEdgeNode(no[3], no[7]);
 
-      CheckRefIsoFace(elem, no[0], no[1], no[5], no[4], mid01, mid15, mid45,
+      CheckRefIsoFace(ref, elem, no[0], no[1], no[5], no[4], mid01, mid15, mid45,
                       mid04, refinements, elemToRef, conflicts);
-      CheckRefIsoFace(elem, no[2], no[3], no[7], no[6], mid23, mid37, mid67,
+      CheckRefIsoFace(ref, elem, no[2], no[3], no[7], no[6], mid23, mid37, mid67,
                       mid26, refinements, elemToRef, conflicts);
    }
    else if (ref_type == Refinement::YZ) // YZ split
@@ -1963,18 +2028,18 @@ void ParNCMesh::CheckRefinement(int elem, char ref_type,
       const int mid26 = GetMidEdgeNode(no[2], no[6]);
       const int mid37 = GetMidEdgeNode(no[3], no[7]);
 
-      CheckRefAnisoFace(elem, no[4], no[0], no[1], no[5], refinements,
+      CheckRefAnisoFace(ref, elem, no[4], no[0], no[1], no[5], refinements,
                         elemToRef, conflicts);
-      CheckRefAnisoFace(elem, no[0], no[3], no[2], no[1], refinements,
+      CheckRefAnisoFace(ref, elem, no[0], no[3], no[2], no[1], refinements,
                         elemToRef, conflicts);
-      CheckRefAnisoFace(elem, no[3], no[7], no[6], no[2], refinements,
+      CheckRefAnisoFace(ref, elem, no[3], no[7], no[6], no[2], refinements,
                         elemToRef, conflicts);
-      CheckRefAnisoFace(elem, no[7], no[4], no[5], no[6], refinements,
+      CheckRefAnisoFace(ref, elem, no[7], no[4], no[5], no[6], refinements,
                         elemToRef, conflicts);
 
-      CheckRefIsoFace(elem, no[1], no[2], no[6], no[5], mid12, mid26, mid56,
+      CheckRefIsoFace(ref, elem, no[1], no[2], no[6], no[5], mid12, mid26, mid56,
                       mid15, refinements, elemToRef, conflicts);
-      CheckRefIsoFace(elem, no[3], no[0], no[4], no[7], mid30, mid04, mid74,
+      CheckRefIsoFace(ref, elem, no[3], no[0], no[4], no[7], mid30, mid04, mid74,
                       mid37, refinements, elemToRef, conflicts);
    }
    else if (ref_type == Refinement::XYZ) // XYZ split
@@ -1994,17 +2059,17 @@ void ParNCMesh::CheckRefinement(int elem, char ref_type,
       const int mid26 = GetMidEdgeNode(no[2], no[6]);
       const int mid37 = GetMidEdgeNode(no[3], no[7]);
 
-      CheckRefIsoFace(elem, no[3], no[2], no[1], no[0], mid23, mid12, mid01,
+      CheckRefIsoFace(ref, elem, no[3], no[2], no[1], no[0], mid23, mid12, mid01,
                       mid30, refinements, elemToRef, conflicts);
-      CheckRefIsoFace(elem, no[0], no[1], no[5], no[4], mid01, mid15, mid45,
+      CheckRefIsoFace(ref, elem, no[0], no[1], no[5], no[4], mid01, mid15, mid45,
                       mid04, refinements, elemToRef, conflicts);
-      CheckRefIsoFace(elem, no[1], no[2], no[6], no[5], mid12, mid26, mid56,
+      CheckRefIsoFace(ref, elem, no[1], no[2], no[6], no[5], mid12, mid26, mid56,
                       mid15, refinements, elemToRef, conflicts);
-      CheckRefIsoFace(elem, no[2], no[3], no[7], no[6], mid23, mid37, mid67,
+      CheckRefIsoFace(ref, elem, no[2], no[3], no[7], no[6], mid23, mid37, mid67,
                       mid26, refinements, elemToRef, conflicts);
-      CheckRefIsoFace(elem, no[3], no[0], no[4], no[7], mid30, mid04, mid74,
+      CheckRefIsoFace(ref, elem, no[3], no[0], no[4], no[7], mid30, mid04, mid74,
                       mid37, refinements, elemToRef, conflicts);
-      CheckRefIsoFace(elem, no[4], no[5], no[6], no[7], mid45, mid56, mid67,
+      CheckRefIsoFace(ref, elem, no[4], no[5], no[6], no[7], mid45, mid56, mid67,
                       mid74, refinements, elemToRef, conflicts);
    }
    else
@@ -2053,7 +2118,7 @@ void ParNCMesh::Refine(const Array<Refinement> &refinements)
       ElementNeighborProcessors(elem, ranks);
       for (int j = 0; j < ranks.Size(); j++)
       {
-         send_ref[ranks[j]].AddRefinement(elem, ref.GetType());
+         send_ref[ranks[j]].AddRefinement(elem, ref);
       }
    }
 
@@ -2063,8 +2128,9 @@ void ParNCMesh::Refine(const Array<Refinement> &refinements)
    // do local refinements
    for (int i = 0; i < refinements.Size(); i++)
    {
-      const Refinement &ref = refinements[i];
-      NCMesh::RefineElement(leaf_elements[ref.index], ref.GetType());
+      Refinement ref_i = refinements[i];
+      ref_i.index = leaf_elements[refinements[i].index];
+      NCMesh::RefineElement(ref_i);
    }
 
    // receive (ghost layer) refinements from all neighbors
@@ -2080,7 +2146,9 @@ void ParNCMesh::Refine(const Array<Refinement> &refinements)
       // do the ghost refinements
       for (int i = 0; i < msg.Size(); i++)
       {
-         NCMesh::RefineElement(msg.elements[i], msg.values[i]);
+         Refinement ghost_ref(msg.elements[i], msg.values[i].ref_type);
+         ghost_ref.SetScaleForType(msg.values[i].scale);
+         NCMesh::RefineElement(ghost_ref);
       }
    }
 

mesh/pncmesh.hpp

@@ -497,11 +497,27 @@ protected: // implementation
    /** Used by ParNCMesh::Refine() to inform neighbors about refinements at
     *  the processor boundary. This keeps their ghost layers synchronized.
     */
-   class NeighborRefinementMessage : public ElementValueMessage<char, false,
-      VarMessageTag::NEIGHBOR_REFINEMENT_VM>
+   struct NeighborRefinement
+   {
+      char ref_type;
+      real_t scale[3];
+   };
+
+   class NeighborRefinementMessage
+      : public ElementValueMessage<NeighborRefinement, false,
+        VarMessageTag::NEIGHBOR_REFINEMENT_VM>
    {
    public:
-      void AddRefinement(int elem, char ref_type) { Add(elem, ref_type); }
+      void AddRefinement(int elem, const Refinement &ref)
+      {
+         NeighborRefinement data{};
+         data.ref_type = ref.GetType();
+         for (int i = 0; i < 3; i++)
+         {
+            data.scale[i] = ref.s[i];
+         }
+         Add(elem, data);
+      }
       typedef std::map<int, NeighborRefinementMessage> Map;
    };
 
@@ -602,7 +618,8 @@ protected: // implementation
    /** For the face with ordered vertices vn* and neighboring element @a elem,
        check whether the other neighboring element (if it exists) is marked for
        a horizontal refinement conflicting with a vertical split. */
-   void CheckRefAnisoFace(int elem, int vn1, int vn2, int vn3, int vn4,
+   void CheckRefAnisoFace(const Refinement &ref, int elem,
+                          int vn1, int vn2, int vn3, int vn4,
                           const Array<Refinement> &refinements,
                           const std::map<int, int> &elemToRef,
                           std::set<int> &conflicts);
@@ -611,7 +628,8 @@ protected: // implementation
        neighboring element @a elem, check whether the other neighboring element
        (if it exists) is marked for a refinement conflicting with an isotropic
        refinement of the face. */
-   void CheckRefIsoFace(int elem, int vn1, int vn2, int vn3, int vn4,
+   void CheckRefIsoFace(const Refinement &ref, int elem,
+                        int vn1, int vn2, int vn3, int vn4,
                         int en1, int en2, int en3, int en4,
                         const Array<Refinement> &refinements,
                         const std::map<int, int> &elemToRef,
@@ -622,9 +640,8 @@ protected: // implementation
                               const std::map<int, int> &elemToRef,
                               std::set<int> &conflicts);
 
-   /** Check whether the refinement of the element with index @a elem and type
-       @a ref_type would cause a conflict. */
-   void CheckRefinement(int elem, char ref_type,
+   /// Check whether the input refinement would cause a conflict.
+   void CheckRefinement(int elem, const Refinement &ref,
                         const Array<Refinement> &refinements,
                         const std::map<int, int> &elemToRef,
                         std::set<int> &conflicts);

miniapps/meshing/CMakeLists.txt

@@ -151,6 +151,10 @@ if (MFEM_USE_MPI)
     MAIN phpref.cpp
     LIBRARIES mfem)
 
+  add_mfem_miniapp(pref321
+    MAIN pref321.cpp
+    LIBRARIES mfem)
+
   # Add parallel tests.
   if (MFEM_ENABLE_TESTING)
     set(PARALLEL_TESTS
@@ -160,6 +164,7 @@ if (MFEM_USE_MPI)
       fit-node-position
       pminimal-surface
       phpref
+      pref321
     )
     # Meshing miniapps that return MFEM_SKIP_RETURN_VALUE in some cases:
     set(SKIP_TESTS)

miniapps/meshing/makefile

@@ -24,7 +24,7 @@ SEQ_MINIAPPS = mobius-strip klein-bottle toroid trimmer twist mesh-explorer\
 	shaper extruder mesh-optimizer minimal-surface polar-nc reflector\
 	ref321 mesh-quality hpref
 PAR_MINIAPPS = pmesh-optimizer pminimal-surface pmesh-fitting fit-node-position\
-	phpref mesh-bounding-boxes
+	phpref pref321 mesh-bounding-boxes
 ifeq ($(MFEM_USE_MPI),NO)
    MINIAPPS = $(SEQ_MINIAPPS)
 else
@@ -99,6 +99,8 @@ hpref-test-seq: hpref
 	@$(call mfem-test,$<,, Serial hp-refinement)
 phpref-test-par: phpref
 	@$(call mfem-test,$<, $(RUN_MPI), Parallel hp-refinement)
+pref321-test-par: pref321
+	@$(call mfem-test,$<, $(RUN_MPI), Parallel 3:1 refinement)
 mesh-bounding-boxes-test-par: mesh-bounding-boxes
 	@$(call mfem-test,$<, $(RUN_MPI), Parallel bounding boxes)
 ref321-test-seq: ref321

miniapps/meshing/pref321.cpp

@@ -0,0 +1,336 @@
+// 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.
+//
+//      -----------------------------------------------------------------
+//      3:1 Refinement Miniapp: Parallel 3:1 anisotropic mesh refinements
+//      -----------------------------------------------------------------
+//
+// This miniapp performs random 3:1 refinements of a quadrilateral or hexahedral
+// mesh. A diffusion equation is solved in an H1 finite element space defined on
+// the refined mesh, and its continuity is verified across local and shared
+// faces.
+//
+// Compile with: make pref321
+//
+// Sample runs:  mpirun -np 4 pref321 -mm -dim 2 -o 2 -r 100
+//               mpirun -np 4 pref321 -mm -dim 3 -o 2 -r 100
+//               mpirun -np 4 pref321 -m ../../data/star.mesh -o 2 -r 100
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+real_t CheckH1Continuity(ParGridFunction &x);
+
+// Find the two children of parent element `elem` after its refinement in one
+// direction.
+void FindChildren(const Mesh &mesh, int elem, Array<int> &children)
+{
+   const CoarseFineTransformations &cf = mesh.ncmesh->GetRefinementTransforms();
+   MFEM_ASSERT(mesh.GetNE() == cf.embeddings.Size(), "");
+
+   // Note that row `elem` of the table constructed by cf.MakeCoarseToFineTable
+   // is an alternative to this global loop, but constructing the table is also
+   // a global operation with global storage.
+   for (int i = 0; i < mesh.GetNE(); i++)
+   {
+      const int p = cf.embeddings[i].parent;
+      if (p == elem)
+      {
+         children.Append(i);
+      }
+   }
+}
+
+// Refine 3:1 via 2 refinements with scalings 2/3 and 1/2.
+void Refine31(Mesh &mesh, int elem, char type)
+{
+   Array<Refinement> refs; // Refinement is defined in ncmesh.hpp
+   refs.Append(Refinement(elem, type, 2.0 / 3.0));
+   mesh.GeneralRefinement(refs);
+
+   // Find the elements with parent `elem`
+   Array<int> children;
+   FindChildren(mesh, elem, children);
+   MFEM_ASSERT(children.Size() == 2, "");
+
+   const int elem1 = children[0];
+
+   refs.SetSize(0);
+   refs.Append(Refinement(elem1, type)); // Default scaling of 0.5
+   mesh.GeneralRefinement(refs);
+}
+
+// Randomly select elements for 3:1 refinements in random directions.
+void TestAnisoRefRandom(int num_refs, int dim, ParMesh &mesh, int myid,
+                        int seed = 0)
+{
+   std::mt19937 gen(seed);
+   for (int i = 0; i < num_refs; i++)
+   {
+      const int elem = gen() % mesh.GetNE();
+      const int t = gen() % dim;
+      auto type = t == 0 ? Refinement::X :
+                  (t == 1 ? Refinement::Y : Refinement::Z);
+
+      // In 3D, check for conflicts in the parallel refinements.
+      if (dim == 3)
+      {
+         std::set<int> conflicts; // Indices in refs of conflicting elements
+         Array<Refinement> refs;
+         refs.Append(Refinement(elem, type));
+         const bool conflict = mesh.AnisotropicConflict(refs, conflicts);
+         if (conflict)
+         {
+            if (myid == 0)
+               cout << "Anisotropic conflict on iteration " << i
+                    << ", retrying\n";
+            i--;
+            continue;
+         }
+      }
+
+      Refine31(mesh, elem, type);
+   }
+
+   mesh.EnsureNodes();
+   mesh.SetScaledNCMesh();
+}
+
+int main(int argc, char *argv[])
+{
+   Mpi::Init(argc, argv);
+   Hypre::Init();
+
+   const int num_procs = Mpi::WorldSize();
+   const int myid = Mpi::WorldRank();
+
+   // 1. Parse command-line options.
+   const char *mesh_file = "../../data/star.mesh";
+   int order = 1;
+   bool visualization = true;
+   bool makeMesh = false;
+   int num_refs = 1;
+   int tdim = 2;  // Mesh dimension for Cartesian meshes.
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree).");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&makeMesh, "-mm", "--make-mesh", "-no-mm",
+                  "--no-make-mesh", "Create Cartesian mesh");
+   args.AddOption(&tdim, "-dim", "--dimension", "Dimension for Cartesian mesh");
+   args.AddOption(&num_refs, "-r", "--refs", "Number of 3:1 refinements");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (myid == 0)
+      {
+         args.PrintUsage(cout);
+      }
+      return 1;
+   }
+   if (myid == 0)
+   {
+      args.PrintOptions(cout);
+   }
+
+   // 2. Create or read the serial mesh on all ranks, then apply the same
+   //    deterministic 3:1 refinement sequence before partitioning it.
+   Mesh mesh;
+   if (makeMesh)
+   {
+      mesh = tdim == 3 ? Mesh::MakeCartesian3D(2, 2, 2, Element::HEXAHEDRON) :
+             Mesh::MakeCartesian2D(2, 2, Element::QUADRILATERAL);
+   }
+   else
+   {
+      mesh = Mesh::LoadFromFile(mesh_file, 1, 1);
+   }
+
+   const int dim = mesh.Dimension();
+
+   mesh.EnsureNCMesh();
+   mesh.SetScaledNCMesh();
+
+   // 3. Partition the refined serial mesh.
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
+
+   TestAnisoRefRandom(num_refs, dim, pmesh, myid, myid);
+
+   // 4. Define a parallel H1 finite element space and report its global size.
+   H1_FECollection fec(order, dim);
+   ParFiniteElementSpace fespace(&pmesh, &fec);
+   if (myid == 0)
+   {
+      cout << "Number of finite element unknowns: "
+           << fespace.GlobalTrueVSize() << endl;
+   }
+
+   // 5. Assemble and solve the Poisson problem, following ex1p.
+   ParGridFunction x(&fespace);
+   x = 0.0;
+
+   ParLinearForm b(&fespace);
+   ConstantCoefficient one(1.0);
+   b.AddDomainIntegrator(new DomainLFIntegrator(one));
+   b.Assemble();
+
+   ParBilinearForm a(&fespace);
+   a.AddDomainIntegrator(new DiffusionIntegrator());
+   a.Assemble();
+
+   OperatorPtr A;
+   Vector B, X;
+   Array<int> ess_tdof_list;
+   if (pmesh.bdr_attributes.Size())
+   {
+      Array<int> ess_bdr(pmesh.bdr_attributes.Max());
+      ess_bdr = 0;
+      pmesh.MarkExternalBoundaries(ess_bdr);
+      fespace.GetEssentialTrueDofs(ess_bdr, ess_tdof_list);
+   }
+
+   a.FormLinearSystem(ess_tdof_list, x, b, A, X, B);
+
+   HypreBoomerAMG M;
+   CGSolver cg(MPI_COMM_WORLD);
+   cg.SetPreconditioner(M);
+   cg.SetOperator(*A);
+   cg.SetRelTol(1e-12);
+   cg.SetMaxIter(2000);
+   cg.SetPrintLevel(1);
+   cg.Mult(B, X);
+
+   a.RecoverFEMSolution(X, b, x);
+
+   // 6. Verify the continuity of the solution in H1 over local and shared
+   //    faces and compute the global maximum jump.
+   const real_t h1err = CheckH1Continuity(x);
+   if (myid == 0)
+   {
+      cout << "Error of H1 continuity: " << h1err << endl;
+   }
+   MFEM_VERIFY(h1err < 1.0e-7, "H1 discontinuity found");
+
+   // 7. Save the refined mesh and the solution in parallel. This output can
+   //    be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
+   {
+      ostringstream mesh_name, sol_name;
+      mesh_name << "mesh." << setfill('0') << setw(6) << myid;
+      sol_name << "sol." << setfill('0') << setw(6) << myid;
+
+      ofstream mesh_ofs(mesh_name.str().c_str());
+      mesh_ofs.precision(8);
+      pmesh.Print(mesh_ofs);
+
+      ofstream sol_ofs(sol_name.str().c_str());
+      sol_ofs.precision(8);
+      x.Save(sol_ofs);
+   }
+
+   // 8. Send the parallel solution to GLVis.
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int visport = 19916;
+      socketstream sol_sock(vishost, visport);
+      sol_sock << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock.precision(8);
+      sol_sock << "solution\n" << pmesh << x << flush;
+   }
+
+   return 0;
+}
+
+real_t CheckH1Continuity(ParGridFunction &x)
+{
+   const ParFiniteElementSpace *pfes = x.ParFESpace();
+   ParMesh *pmesh = pfes->GetParMesh();
+   const int dim = pmesh->Dimension();
+
+   real_t errorMax = 0.0;
+
+   // Shared-face values require face-neighbor data.
+   x.ExchangeFaceNbrData();
+
+   // First handle faces for which both elements are local to this rank.
+   for (int f = 0; f < pmesh->GetNumFaces(); f++)
+   {
+      const auto info = pmesh->GetFaceInformation(f);
+      if (!info.IsLocal())
+      {
+         continue;
+      }
+
+      FaceElementTransformations *FT = pmesh->GetFaceElementTransformations(f);
+      const int faceOrder = dim == 3 ? pfes->GetFaceOrder(f) :
+                            pfes->GetEdgeOrder(f);
+      const IntegrationRule &ir = IntRules.Get(FT->FaceGeom, 2 * faceOrder);
+
+      for (int i = 0; i < ir.GetNPoints(); i++)
+      {
+         const IntegrationPoint &fip = ir.IntPoint(i);
+         IntegrationPoint ip1, ip2;
+
+         FT->Loc1.Transform(fip, ip1);
+         FT->Loc2.Transform(fip, ip2);
+
+         const real_t v1 = x.GetValue(*FT->Elem1, ip1);
+         const real_t v2 = x.GetValue(*FT->Elem2, ip2);
+         errorMax = std::max(errorMax, std::abs(v1 - v2));
+      }
+   }
+
+   // Then check partition interfaces. Conforming shared faces are handled on
+   // the lower-rank side, while shared slave nonconforming faces are handled
+   // only on the slave side and therefore do not need additional filtering.
+   for (int sf = 0; sf < pmesh->GetNSharedFaces(); sf++)
+   {
+      const int f = pmesh->GetSharedFace(sf);
+      const auto info = pmesh->GetFaceInformation(f);
+      if (!info.IsShared())
+      {
+         continue;
+      }
+
+      FaceElementTransformations *FT = pmesh->GetSharedFaceTransformations(sf);
+      const int faceOrder = dim == 3 ? pfes->GetFaceOrder(f) :
+                            pfes->GetEdgeOrder(f);
+      const IntegrationRule &ir = IntRules.Get(FT->FaceGeom, 2 * faceOrder);
+
+      for (int i = 0; i < ir.GetNPoints(); i++)
+      {
+         const IntegrationPoint &fip = ir.IntPoint(i);
+         IntegrationPoint ip1, ip2;
+
+         FT->Loc1.Transform(fip, ip1);
+         FT->Loc2.Transform(fip, ip2);
+
+         const real_t v1 = x.GetValue(*FT->Elem1, ip1);
+         const real_t v2 = x.GetValue(*FT->Elem2, ip2);
+         errorMax = std::max(errorMax, std::abs(v1 - v2));
+      }
+   }
+
+   MPI_Allreduce(MPI_IN_PLACE, &errorMax, 1, MFEM_MPI_REAL_T, MPI_MAX,
+                 pmesh->GetComm());
+
+   return errorMax;
+}

miniapps/meshing/ref321.cpp

@@ -71,22 +71,14 @@ void Refine31(Mesh & mesh, int elem, char type)
    mesh.GeneralRefinement(refs);
 }
 
-// Deterministic, somewhat random integer generator
-int MyRand(int & s)
-{
-   s++;
-   const double a = 1000 * sin(s * 1.1234 * M_PI);
-   return int(std::abs(a));
-}
-
 // Randomly select elements for 3:1 refinements in random directions.
-void TestAnisoRefRandom(int iter, int dim, Mesh & mesh)
+void TestAnisoRefRandom(int num_refs, int dim, Mesh & mesh)
 {
-   int seed = 0;
-   for (int i = 0; i < iter; i++)
+   std::mt19937 gen(1);
+   for (int i = 0; i < num_refs; i++)
    {
-      const int elem = MyRand(seed) % mesh.GetNE();
-      const int t = MyRand(seed) % dim;
+      const auto elem = gen() % mesh.GetNE();
+      const auto t = gen() % dim;
       auto type = t == 0 ? Refinement::X :
                   (t == 1 ? Refinement::Y : Refinement::Z);
       Refine31(mesh, elem, type);

tests/unit/CMakeLists.txt

@@ -140,6 +140,7 @@ set(UNIT_TESTS_SRCS
   fem/test_lor_batched.cpp
   fem/test_lor_dg.cpp
   fem/test_lor.cpp
+  fem/test_mixedsesqform.cpp
   fem/test_nonlinearform.cpp
   fem/test_operatorjacobismoother.cpp
   fem/test_oscillation.cpp
@@ -168,6 +169,20 @@ set(UNIT_TESTS_SRCS
   fem/test_transfer.cpp
   fem/test_var_order.cpp
   fem/test_white_noise.cpp
+  fem/specializations/test_diffusion_integ.cpp
+  fem/specializations/test_mass_integ.cpp
+  fem/specializations/test_convection_integ.cpp
+  fem/specializations/test_vecmass_integ.cpp
+  fem/specializations/test_curlcurl_integ.cpp
+  fem/specializations/test_vecdiffusion_integ.cpp
+  fem/specializations/test_dgtrace_integ.cpp
+  fem/specializations/test_dgdiffusion_integ.cpp
+  fem/specializations/test_dgmassinv.cpp
+  fem/specializations/test_qinterp_det.cpp
+  fem/specializations/test_qinterp_eval.cpp
+  fem/specializations/test_qinterp_grad.cpp
+  fem/specializations/test_qinterp_tensoreval.cpp
+  fem/specializations/test_qinterp_eval_hdiv.cpp
   enzyme/compatibility.cpp
   # The following are tested separately (keep the comment as a reminder).
   # This list can be updated using (in bash):

tests/unit/data/hex-nc-cross.mesh

@@ -0,0 +1,379 @@
+MFEM NC mesh v1.0
+
+# NCMesh supported geometry types:
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+# PRISM       = 6
+# PYRAMID     = 7
+
+dimension
+3
+
+rank
+0
+
+# rank attr geom ref_type nodes/children
+elements
+75
+0 1 5 0 0 1 5 4 16 17 21 20
+0 1 5 0 16 17 21 20 32 33 37 36
+-1 1 5 7 59 60 61 62 63 64 65 66
+0 1 5 0 1 2 6 5 17 18 22 21
+-1 1 5 7 27 28 29 30 31 32 33 34
+0 1 5 0 21 22 26 25 37 38 42 41
+-1 1 5 7 43 44 45 46 47 48 49 50
+0 1 5 0 4 5 9 8 20 21 25 24
+0 1 5 0 8 9 13 12 24 25 29 28
+0 1 5 0 24 25 29 28 40 41 45 44
+0 1 5 0 9 10 14 13 25 26 30 29
+-1 1 5 7 67 68 69 70 71 72 73 74
+0 1 5 0 41 42 46 45 57 58 62 61
+0 1 5 0 40 41 45 44 56 57 61 60
+0 1 5 0 36 37 41 40 52 53 57 56
+-1 1 5 7 35 36 37 38 39 40 41 42
+0 1 5 0 32 33 37 36 48 49 53 52
+0 1 5 0 33 34 38 37 49 50 54 53
+0 1 5 0 34 35 39 38 50 51 55 54
+0 1 5 0 38 39 43 42 54 55 59 58
+0 1 5 0 42 43 47 46 58 59 63 62
+0 1 5 0 26 27 31 30 42 43 47 46
+0 1 5 0 10 11 15 14 26 27 31 30
+0 1 5 0 6 7 11 10 22 23 27 26
+-1 1 5 7 51 52 53 54 55 56 57 58
+0 1 5 0 18 19 23 22 34 35 39 38
+0 1 5 0 2 3 7 6 18 19 23 22
+0 1 5 0 5 94 208 99 74 209 214 212
+0 1 5 0 94 6 97 208 209 96 210 214
+0 1 5 0 208 97 10 98 214 210 103 211
+0 1 5 0 99 208 98 9 212 214 211 104
+0 1 5 0 74 209 214 212 21 86 213 102
+0 1 5 0 209 96 210 214 86 22 100 213
+0 1 5 0 214 210 103 211 213 100 26 101
+0 1 5 0 212 214 211 104 102 213 101 25
+0 1 5 0 37 89 269 107 156 270 275 273
+0 1 5 0 89 38 105 269 270 159 271 275
+0 1 5 0 269 105 42 106 275 271 144 272
+0 1 5 0 107 269 106 41 273 275 272 143
+0 1 5 0 156 270 275 273 53 157 274 153
+0 1 5 0 270 159 271 275 157 54 158 274
+0 1 5 0 275 271 144 272 274 158 58 139
+0 1 5 0 273 275 272 143 153 274 139 57
+0 1 5 0 20 70 330 111 83 331 336 334
+0 1 5 0 70 21 102 330 331 82 332 336
+0 1 5 0 330 102 25 110 336 332 109 333
+0 1 5 0 111 330 110 24 334 336 333 114
+0 1 5 0 83 331 336 334 36 78 335 113
+0 1 5 0 331 82 332 336 78 37 107 335
+0 1 5 0 336 332 109 333 335 107 41 112
+0 1 5 0 334 336 333 114 113 335 112 40
+0 1 5 0 22 198 387 100 91 388 393 391
+0 1 5 0 198 23 199 387 388 201 389 393
+0 1 5 0 387 199 27 186 393 389 189 390
+0 1 5 0 100 387 186 26 391 393 390 108
+0 1 5 0 91 388 393 391 38 170 392 105
+0 1 5 0 388 201 389 393 170 39 176 392
+0 1 5 0 393 389 189 390 392 176 43 177
+0 1 5 0 391 393 390 108 105 392 177 42
+0 1 5 0 17 84 444 69 81 445 450 448
+0 1 5 0 84 18 85 444 445 90 446 450
+0 1 5 0 444 85 22 86 450 446 91 447
+0 1 5 0 69 444 86 21 448 450 447 82
+0 1 5 0 81 445 450 448 33 87 449 77
+0 1 5 0 445 90 446 450 87 34 88 449
+0 1 5 0 450 446 91 447 449 88 38 89
+0 1 5 0 448 450 447 82 77 449 89 37
+0 1 5 0 25 101 497 121 109 498 503 501
+0 1 5 0 101 26 133 497 498 108 499 503
+0 1 5 0 497 133 30 134 503 499 138 500
+0 1 5 0 121 497 134 29 501 503 500 129
+0 1 5 0 109 498 503 501 41 106 502 126
+0 1 5 0 498 108 499 503 106 42 136 502
+0 1 5 0 503 499 138 500 502 136 46 137
+0 1 5 0 501 503 500 129 126 502 137 45
+
+# attr geom nodes
+boundary
+72
+1 3 4 5 1 0
+1 3 0 1 17 16
+1 3 4 0 16 20
+1 3 16 17 33 32
+1 3 20 16 32 36
+1 3 5 6 2 1
+1 3 1 2 18 17
+1 3 8 9 5 4
+1 3 8 4 20 24
+1 3 12 13 9 8
+1 3 13 12 28 29
+1 3 12 8 24 28
+1 3 29 28 44 45
+1 3 28 24 40 44
+1 3 13 14 10 9
+1 3 14 13 29 30
+1 3 46 45 61 62
+1 3 57 58 62 61
+1 3 45 44 60 61
+1 3 44 40 56 60
+1 3 56 57 61 60
+1 3 40 36 52 56
+1 3 52 53 57 56
+1 3 32 33 49 48
+1 3 36 32 48 52
+1 3 48 49 53 52
+1 3 33 34 50 49
+1 3 49 50 54 53
+1 3 34 35 51 50
+1 3 35 39 55 51
+1 3 50 51 55 54
+1 3 39 43 59 55
+1 3 54 55 59 58
+1 3 43 47 63 59
+1 3 47 46 62 63
+1 3 58 59 63 62
+1 3 27 31 47 43
+1 3 31 30 46 47
+1 3 14 15 11 10
+1 3 11 15 31 27
+1 3 15 14 30 31
+1 3 10 11 7 6
+1 3 7 11 27 23
+1 3 18 19 35 34
+1 3 19 23 39 35
+1 3 6 7 3 2
+1 3 2 3 19 18
+1 3 3 7 23 19
+2 3 99 208 94 5
+2 3 208 97 6 94
+2 3 98 10 97 208
+2 3 9 98 208 99
+2 3 53 157 274 153
+2 3 157 54 158 274
+2 3 274 158 58 139
+2 3 153 274 139 57
+2 3 111 20 83 334
+2 3 24 111 334 114
+2 3 334 83 36 113
+2 3 114 334 113 40
+2 3 23 199 389 201
+2 3 199 27 189 389
+2 3 201 389 176 39
+2 3 389 189 43 176
+2 3 17 84 445 81
+2 3 84 18 90 445
+2 3 81 445 87 33
+2 3 445 90 34 87
+2 3 30 134 500 138
+2 3 134 29 129 500
+2 3 138 500 137 46
+2 3 500 129 45 137
+
+# vert_id p1 p2
+vertex_parents
+102
+69 17 21
+70 20 21
+74 5 21
+77 33 37
+78 36 37
+81 17 33
+82 21 37
+83 20 36
+84 17 18
+85 18 22
+86 21 22
+87 33 34
+88 34 38
+89 37 38
+90 18 34
+91 22 38
+94 5 6
+96 6 22
+97 6 10
+98 9 10
+99 5 9
+100 22 26
+101 25 26
+102 21 25
+103 10 26
+104 9 25
+105 38 42
+106 41 42
+107 37 41
+108 26 42
+109 25 41
+110 24 25
+111 20 24
+112 40 41
+113 36 40
+114 24 40
+121 25 29
+126 41 45
+129 29 45
+133 26 30
+134 29 30
+136 42 46
+137 45 46
+138 30 46
+139 57 58
+143 41 57
+144 42 58
+153 53 57
+156 37 53
+157 53 54
+158 54 58
+159 38 54
+170 38 39
+176 39 43
+177 42 43
+186 26 27
+189 27 43
+198 22 23
+199 23 27
+201 23 39
+208 97 99
+209 74 96
+210 96 103
+211 103 104
+212 74 104
+213 100 102
+214 209 211
+269 105 107
+270 156 159
+271 144 159
+272 143 144
+273 143 156
+274 153 158
+275 270 272
+330 102 111
+331 82 83
+332 82 109
+333 109 114
+334 83 114
+335 107 113
+336 331 333
+387 100 199
+388 91 201
+389 189 201
+390 108 189
+391 91 108
+392 105 176
+393 388 390
+444 69 85
+445 81 90
+446 90 91
+447 82 91
+448 81 82
+449 77 88
+450 445 447
+497 121 133
+498 108 109
+499 108 138
+500 129 138
+501 109 129
+502 126 136
+503 498 500
+
+# root element orientation
+root_state
+27
+0
+1
+1
+15
+15
+6
+6
+22
+15
+8
+12
+10
+10
+18
+18
+13
+7
+22
+22
+15
+16
+16
+16
+7
+8
+6
+21
+
+# top-level node coordinates
+coordinates
+64
+3
+0 0 0
+0.33333333 0 0
+0.66666667 0 0
+1 0 0
+0 0.33333333 0
+0.33333333 0.33333333 0
+0.66666667 0.33333333 0
+1 0.33333333 0
+0 0.66666667 0
+0.33333333 0.66666667 0
+0.66666667 0.66666667 0
+1 0.66666667 0
+0 1 0
+0.33333333 1 0
+0.66666667 1 0
+1 1 0
+0 0 0.33333333
+0.33333333 0 0.33333333
+0.66666667 0 0.33333333
+1 0 0.33333333
+0 0.33333333 0.33333333
+0.33333333 0.33333333 0.33333333
+0.66666667 0.33333333 0.33333333
+1 0.33333333 0.33333333
+0 0.66666667 0.33333333
+0.33333333 0.66666667 0.33333333
+0.66666667 0.66666667 0.33333333
+1 0.66666667 0.33333333
+0 1 0.33333333
+0.33333333 1 0.33333333
+0.66666667 1 0.33333333
+1 1 0.33333333
+0 0 0.66666667
+0.33333333 0 0.66666667
+0.66666667 0 0.66666667
+1 0 0.66666667
+0 0.33333333 0.66666667
+0.33333333 0.33333333 0.66666667
+0.66666667 0.33333333 0.66666667
+1 0.33333333 0.66666667
+0 0.66666667 0.66666667
+0.33333333 0.66666667 0.66666667
+0.66666667 0.66666667 0.66666667
+1 0.66666667 0.66666667
+0 1 0.66666667
+0.33333333 1 0.66666667
+0.66666667 1 0.66666667
+1 1 0.66666667
+0 0 1
+0.33333333 0 1
+0.66666667 0 1
+1 0 1
+0 0.33333333 1
+0.33333333 0.33333333 1
+0.66666667 0.33333333 1
+1 0.33333333 1
+0 0.66666667 1
+0.33333333 0.66666667 1
+0.66666667 0.66666667 1
+1 0.66666667 1
+0 1 1
+0.33333333 1 1
+0.66666667 1 1
+1 1 1
+
+mfem_mesh_end

tests/unit/data/tet-nc-cross.mesh

@@ -0,0 +1,555 @@
+MFEM NC mesh v1.0
+
+# NCMesh supported geometry types:
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+# PRISM       = 6
+# PYRAMID     = 7
+
+dimension
+3
+
+rank
+0
+
+# rank attr geom ref_type nodes/children
+elements
+258
+0 1 4 0 21 0 5 1
+0 1 4 0 21 0 1 17
+0 1 4 0 21 0 17 16
+0 1 4 0 21 0 4 5
+0 1 4 0 21 0 20 4
+0 1 4 0 21 0 16 20
+0 1 4 0 22 1 6 2
+0 1 4 0 22 1 2 18
+0 1 4 0 22 1 18 17
+0 1 4 0 22 1 5 6
+0 1 4 0 22 1 21 5
+0 1 4 0 22 1 17 21
+0 1 4 0 23 2 7 3
+0 1 4 0 23 2 3 19
+0 1 4 0 23 2 19 18
+0 1 4 0 23 2 6 7
+0 1 4 0 23 2 22 6
+0 1 4 0 23 2 18 22
+0 1 4 0 25 4 9 5
+0 1 4 0 25 4 5 21
+0 1 4 0 25 4 21 20
+0 1 4 0 25 4 8 9
+0 1 4 0 25 4 24 8
+0 1 4 0 25 4 20 24
+-1 1 4 7 170 171 172 173 174 175 176 177
+0 1 4 0 26 5 6 22
+0 1 4 0 26 5 22 21
+-1 1 4 7 162 163 164 165 166 167 168 169
+0 1 4 0 26 5 25 9
+0 1 4 0 26 5 21 25
+0 1 4 0 27 6 11 7
+0 1 4 0 27 6 7 23
+0 1 4 0 27 6 23 22
+0 1 4 0 27 6 10 11
+0 1 4 0 27 6 26 10
+0 1 4 0 27 6 22 26
+0 1 4 0 29 8 13 9
+0 1 4 0 29 8 9 25
+0 1 4 0 29 8 25 24
+0 1 4 0 29 8 12 13
+0 1 4 0 29 8 28 12
+0 1 4 0 29 8 24 28
+0 1 4 0 30 9 14 10
+0 1 4 0 30 9 10 26
+0 1 4 0 30 9 26 25
+0 1 4 0 30 9 13 14
+0 1 4 0 30 9 29 13
+0 1 4 0 30 9 25 29
+0 1 4 0 31 10 15 11
+0 1 4 0 31 10 11 27
+0 1 4 0 31 10 27 26
+0 1 4 0 31 10 14 15
+0 1 4 0 31 10 30 14
+0 1 4 0 31 10 26 30
+0 1 4 0 37 16 21 17
+0 1 4 0 37 16 17 33
+0 1 4 0 37 16 33 32
+0 1 4 0 37 16 20 21
+0 1 4 0 37 16 36 20
+0 1 4 0 37 16 32 36
+0 1 4 0 38 17 22 18
+-1 1 4 7 226 227 228 229 230 231 232 233
+-1 1 4 7 234 235 236 237 238 239 240 241
+0 1 4 0 38 17 21 22
+0 1 4 0 38 17 37 21
+0 1 4 0 38 17 33 37
+0 1 4 0 39 18 23 19
+0 1 4 0 39 18 19 35
+0 1 4 0 39 18 35 34
+0 1 4 0 39 18 22 23
+0 1 4 0 39 18 38 22
+0 1 4 0 39 18 34 38
+0 1 4 0 41 20 25 21
+0 1 4 0 41 20 21 37
+0 1 4 0 41 20 37 36
+0 1 4 0 41 20 24 25
+-1 1 4 7 202 203 204 205 206 207 208 209
+-1 1 4 7 194 195 196 197 198 199 200 201
+0 1 4 0 42 21 26 22
+0 1 4 0 42 21 22 38
+0 1 4 0 42 21 38 37
+0 1 4 0 42 21 25 26
+0 1 4 0 42 21 41 25
+0 1 4 0 42 21 37 41
+-1 1 4 7 210 211 212 213 214 215 216 217
+-1 1 4 7 218 219 220 221 222 223 224 225
+0 1 4 0 43 22 39 38
+0 1 4 0 43 22 26 27
+0 1 4 0 43 22 42 26
+0 1 4 0 43 22 38 42
+0 1 4 0 45 24 29 25
+0 1 4 0 45 24 25 41
+0 1 4 0 45 24 41 40
+0 1 4 0 45 24 28 29
+0 1 4 0 45 24 44 28
+0 1 4 0 45 24 40 44
+0 1 4 0 46 25 30 26
+0 1 4 0 46 25 26 42
+0 1 4 0 46 25 42 41
+-1 1 4 7 250 251 252 253 254 255 256 257
+-1 1 4 7 242 243 244 245 246 247 248 249
+0 1 4 0 46 25 41 45
+0 1 4 0 47 26 31 27
+0 1 4 0 47 26 27 43
+0 1 4 0 47 26 43 42
+0 1 4 0 47 26 30 31
+0 1 4 0 47 26 46 30
+0 1 4 0 47 26 42 46
+0 1 4 0 53 32 37 33
+0 1 4 0 53 32 33 49
+0 1 4 0 53 32 49 48
+0 1 4 0 53 32 36 37
+0 1 4 0 53 32 52 36
+0 1 4 0 53 32 48 52
+0 1 4 0 54 33 38 34
+0 1 4 0 54 33 34 50
+0 1 4 0 54 33 50 49
+0 1 4 0 54 33 37 38
+0 1 4 0 54 33 53 37
+0 1 4 0 54 33 49 53
+0 1 4 0 55 34 39 35
+0 1 4 0 55 34 35 51
+0 1 4 0 55 34 51 50
+0 1 4 0 55 34 38 39
+0 1 4 0 55 34 54 38
+0 1 4 0 55 34 50 54
+0 1 4 0 57 36 41 37
+0 1 4 0 57 36 37 53
+0 1 4 0 57 36 53 52
+0 1 4 0 57 36 40 41
+0 1 4 0 57 36 56 40
+0 1 4 0 57 36 52 56
+0 1 4 0 58 37 42 38
+0 1 4 0 58 37 38 54
+-1 1 4 7 178 179 180 181 182 183 184 185
+0 1 4 0 58 37 41 42
+0 1 4 0 58 37 57 41
+-1 1 4 7 186 187 188 189 190 191 192 193
+0 1 4 0 59 38 43 39
+0 1 4 0 59 38 39 55
+0 1 4 0 59 38 55 54
+0 1 4 0 59 38 42 43
+0 1 4 0 59 38 58 42
+0 1 4 0 59 38 54 58
+0 1 4 0 61 40 45 41
+0 1 4 0 61 40 41 57
+0 1 4 0 61 40 57 56
+0 1 4 0 61 40 44 45
+0 1 4 0 61 40 60 44
+0 1 4 0 61 40 56 60
+0 1 4 0 62 41 46 42
+0 1 4 0 62 41 42 58
+0 1 4 0 62 41 58 57
+0 1 4 0 62 41 45 46
+0 1 4 0 62 41 61 45
+0 1 4 0 62 41 57 61
+0 1 4 0 63 42 47 43
+0 1 4 0 63 42 43 59
+0 1 4 0 63 42 59 58
+0 1 4 0 63 42 46 47
+0 1 4 0 63 42 62 46
+0 1 4 0 63 42 58 62
+0 1 4 0 26 125 132 126
+0 1 4 0 125 5 115 128
+0 1 4 0 132 115 9 133
+0 1 4 0 126 128 133 10
+0 1 4 0 125 133 132 126
+0 1 4 0 125 133 126 128
+0 1 4 0 125 133 128 115
+0 1 4 0 125 133 115 132
+0 1 4 0 26 125 126 127
+0 1 4 0 125 5 128 95
+0 1 4 0 126 128 10 129
+0 1 4 0 127 95 129 6
+0 1 4 0 125 129 126 127
+0 1 4 0 125 129 127 95
+0 1 4 0 125 129 95 128
+0 1 4 0 125 129 128 126
+0 1 4 0 58 305 308 309
+0 1 4 0 305 37 283 262
+0 1 4 0 308 283 54 284
+0 1 4 0 309 262 284 53
+0 1 4 0 305 284 308 309
+0 1 4 0 305 284 309 262
+0 1 4 0 305 284 262 283
+0 1 4 0 305 284 283 308
+0 1 4 0 58 305 309 311
+0 1 4 0 305 37 262 297
+0 1 4 0 309 262 53 298
+0 1 4 0 311 297 298 57
+0 1 4 0 305 298 309 311
+0 1 4 0 305 298 311 297
+0 1 4 0 305 298 297 262
+0 1 4 0 305 298 262 309
+0 1 4 0 41 213 217 219
+0 1 4 0 213 20 191 220
+0 1 4 0 217 191 36 222
+0 1 4 0 219 220 222 40
+0 1 4 0 213 222 217 219
+0 1 4 0 213 222 219 220
+0 1 4 0 213 222 220 191
+0 1 4 0 213 222 191 217
+0 1 4 0 41 213 219 218
+0 1 4 0 213 20 220 124
+0 1 4 0 219 220 40 221
+0 1 4 0 218 124 221 24
+0 1 4 0 213 221 219 218
+0 1 4 0 213 221 218 124
+0 1 4 0 213 221 124 220
+0 1 4 0 213 221 220 219
+0 1 4 0 43 230 231 232
+0 1 4 0 230 22 141 110
+0 1 4 0 231 141 27 140
+0 1 4 0 232 110 140 23
+0 1 4 0 230 140 231 232
+0 1 4 0 230 140 232 110
+0 1 4 0 230 140 110 141
+0 1 4 0 230 140 141 231
+0 1 4 0 43 230 232 233
+0 1 4 0 230 22 110 211
+0 1 4 0 232 110 23 204
+0 1 4 0 233 211 204 39
+0 1 4 0 230 204 232 233
+0 1 4 0 230 204 233 211
+0 1 4 0 230 204 211 110
+0 1 4 0 230 204 110 232
+0 1 4 0 38 193 195 196
+0 1 4 0 193 17 93 197
+0 1 4 0 195 93 18 198
+0 1 4 0 196 197 198 34
+0 1 4 0 193 198 195 196
+0 1 4 0 193 198 196 197
+0 1 4 0 193 198 197 93
+0 1 4 0 193 198 93 195
+0 1 4 0 38 193 196 199
+0 1 4 0 193 17 197 184
+0 1 4 0 196 197 34 200
+0 1 4 0 199 184 200 33
+0 1 4 0 193 200 196 199
+0 1 4 0 193 200 199 184
+0 1 4 0 193 200 184 197
+0 1 4 0 193 200 197 196
+0 1 4 0 46 247 253 252
+0 1 4 0 247 25 239 150
+0 1 4 0 253 239 45 238
+0 1 4 0 252 150 238 29
+0 1 4 0 247 238 253 252
+0 1 4 0 247 238 252 150
+0 1 4 0 247 238 150 239
+0 1 4 0 247 238 239 253
+0 1 4 0 46 247 252 248
+0 1 4 0 247 25 150 165
+0 1 4 0 252 150 29 168
+0 1 4 0 248 165 168 30
+0 1 4 0 247 168 252 248
+0 1 4 0 247 168 248 165
+0 1 4 0 247 168 165 150
+0 1 4 0 247 168 150 252
+
+# attr geom nodes
+boundary
+144
+1 2 0 5 1
+1 2 0 1 17
+1 2 0 17 16
+1 2 0 4 5
+1 2 0 20 4
+1 2 0 16 20
+1 2 1 6 2
+1 2 1 2 18
+1 2 1 18 17
+1 2 1 5 6
+1 2 2 7 3
+1 2 23 3 7
+1 2 2 3 19
+1 2 23 19 3
+1 2 2 19 18
+1 2 2 6 7
+1 2 4 9 5
+1 2 4 8 9
+1 2 4 24 8
+1 2 4 20 24
+1 2 6 11 7
+1 2 27 7 11
+1 2 27 23 7
+1 2 6 10 11
+1 2 8 13 9
+1 2 8 12 13
+1 2 29 13 12
+1 2 8 28 12
+1 2 29 12 28
+1 2 8 24 28
+1 2 9 14 10
+1 2 9 13 14
+1 2 30 14 13
+1 2 30 13 29
+1 2 10 15 11
+1 2 31 11 15
+1 2 31 27 11
+1 2 10 14 15
+1 2 31 15 14
+1 2 31 14 30
+1 2 16 17 33
+1 2 16 33 32
+1 2 16 36 20
+1 2 16 32 36
+1 2 39 19 23
+1 2 18 19 35
+1 2 39 35 19
+1 2 18 35 34
+1 2 45 29 28
+1 2 24 44 28
+1 2 45 28 44
+1 2 24 40 44
+1 2 47 27 31
+1 2 47 43 27
+1 2 47 31 30
+1 2 47 30 46
+1 2 32 33 49
+1 2 32 49 48
+1 2 53 48 49
+1 2 32 52 36
+1 2 32 48 52
+1 2 53 52 48
+1 2 33 34 50
+1 2 33 50 49
+1 2 54 49 50
+1 2 54 53 49
+1 2 55 35 39
+1 2 34 35 51
+1 2 55 51 35
+1 2 34 51 50
+1 2 55 50 51
+1 2 55 54 50
+1 2 57 52 53
+1 2 36 56 40
+1 2 36 52 56
+1 2 57 56 52
+1 2 59 39 43
+1 2 59 55 39
+1 2 59 54 55
+1 2 59 58 54
+1 2 61 56 57
+1 2 61 45 44
+1 2 40 60 44
+1 2 61 44 60
+1 2 40 56 60
+1 2 61 60 56
+1 2 62 57 58
+1 2 62 46 45
+1 2 62 45 61
+1 2 62 61 57
+1 2 63 43 47
+1 2 63 59 43
+1 2 63 58 59
+1 2 63 47 46
+1 2 63 46 62
+1 2 63 62 58
+2 2 5 115 128
+2 2 115 9 133
+2 2 128 133 10
+2 2 133 128 115
+2 2 5 128 95
+2 2 128 10 129
+2 2 95 129 6
+2 2 129 95 128
+2 2 58 309 308
+2 2 308 284 54
+2 2 309 53 284
+2 2 284 308 309
+2 2 58 311 309
+2 2 309 298 53
+2 2 311 57 298
+2 2 298 309 311
+2 2 20 191 220
+2 2 191 36 222
+2 2 220 222 40
+2 2 222 220 191
+2 2 20 220 124
+2 2 220 40 221
+2 2 124 221 24
+2 2 221 124 220
+2 2 43 232 231
+2 2 231 140 27
+2 2 232 23 140
+2 2 140 231 232
+2 2 43 233 232
+2 2 232 204 23
+2 2 233 39 204
+2 2 204 232 233
+2 2 17 93 197
+2 2 93 18 198
+2 2 197 198 34
+2 2 198 197 93
+2 2 17 197 184
+2 2 197 34 200
+2 2 184 200 33
+2 2 200 184 197
+2 2 46 252 253
+2 2 253 238 45
+2 2 252 29 238
+2 2 238 253 252
+2 2 46 248 252
+2 2 252 168 29
+2 2 248 30 168
+2 2 168 252 248
+
+# vert_id p1 p2
+vertex_parents
+54
+93 17 18
+95 5 6
+110 22 23
+115 5 9
+124 20 24
+125 5 26
+126 10 26
+127 6 26
+128 5 10
+129 6 10
+132 9 26
+133 9 10
+140 23 27
+141 22 27
+150 25 29
+165 25 30
+168 29 30
+184 17 33
+191 20 36
+193 17 38
+195 18 38
+196 34 38
+197 17 34
+198 18 34
+199 33 38
+200 33 34
+204 23 39
+211 22 39
+213 20 41
+217 36 41
+218 24 41
+219 40 41
+220 20 40
+221 24 40
+222 36 40
+230 22 43
+231 27 43
+232 23 43
+233 39 43
+238 29 45
+239 25 45
+247 25 46
+248 30 46
+252 29 46
+253 45 46
+262 37 53
+283 37 54
+284 53 54
+297 37 57
+298 53 57
+305 37 58
+308 54 58
+309 53 58
+311 57 58
+
+# top-level node coordinates
+coordinates
+64
+3
+0 0 0
+0.33333333 0 0
+0.66666667 0 0
+1 0 0
+0 0.33333333 0
+0.33333333 0.33333333 0
+0.66666667 0.33333333 0
+1 0.33333333 0
+0 0.66666667 0
+0.33333333 0.66666667 0
+0.66666667 0.66666667 0
+1 0.66666667 0
+0 1 0
+0.33333333 1 0
+0.66666667 1 0
+1 1 0
+0 0 0.33333333
+0.33333333 0 0.33333333
+0.66666667 0 0.33333333
+1 0 0.33333333
+0 0.33333333 0.33333333
+0.33333333 0.33333333 0.33333333
+0.66666667 0.33333333 0.33333333
+1 0.33333333 0.33333333
+0 0.66666667 0.33333333
+0.33333333 0.66666667 0.33333333
+0.66666667 0.66666667 0.33333333
+1 0.66666667 0.33333333
+0 1 0.33333333
+0.33333333 1 0.33333333
+0.66666667 1 0.33333333
+1 1 0.33333333
+0 0 0.66666667
+0.33333333 0 0.66666667
+0.66666667 0 0.66666667
+1 0 0.66666667
+0 0.33333333 0.66666667
+0.33333333 0.33333333 0.66666667
+0.66666667 0.33333333 0.66666667
+1 0.33333333 0.66666667
+0 0.66666667 0.66666667
+0.33333333 0.66666667 0.66666667
+0.66666667 0.66666667 0.66666667
+1 0.66666667 0.66666667
+0 1 0.66666667
+0.33333333 1 0.66666667
+0.66666667 1 0.66666667
+1 1 0.66666667
+0 0 1
+0.33333333 0 1
+0.66666667 0 1
+1 0 1
+0 0.33333333 1
+0.33333333 0.33333333 1
+0.66666667 0.33333333 1
+1 0.33333333 1
+0 0.66666667 1
+0.33333333 0.66666667 1
+0.66666667 0.66666667 1
+1 0.66666667 1
+0 1 1
+0.33333333 1 1
+0.66666667 1 1
+1 1 1
+
+mfem_mesh_end

tests/unit/fem/specializations/test_convection_integ.cpp

@@ -0,0 +1,14 @@
+/// Tests which make sure adding user-defined kernel specializations work.
+/// These tests are compile/link-only tests
+
+#include "mfem.hpp"
+#include "unit_tests.hpp"
+
+#include "fem/integ/bilininteg_convection_kernels.hpp"
+
+TEST_CASE("Convection Kernel Specializations", "[Specializations]")
+{
+   using namespace mfem;
+
+   ConvectionIntegrator::AddSpecialization<2, 2, 4>();
+}

tests/unit/fem/specializations/test_curlcurl_integ.cpp

@@ -0,0 +1,14 @@
+/// Tests which make sure adding user-defined kernel specializations work.
+/// These tests are compile/link-only tests
+
+#include "mfem.hpp"
+#include "unit_tests.hpp"
+
+#include "fem/integ/bilininteg_hcurl_kernels.hpp"
+
+TEST_CASE("CurlCurl Kernel Specializations", "[Specializations]")
+{
+   using namespace mfem;
+
+   CurlCurlIntegrator::AddSpecialization<3, 2, 4>();
+}

tests/unit/fem/specializations/test_dgdiffusion_integ.cpp

@@ -0,0 +1,14 @@
+/// Tests which make sure adding user-defined kernel specializations work.
+/// These tests are compile/link-only tests
+
+#include "mfem.hpp"
+#include "unit_tests.hpp"
+
+#include "fem/integ/bilininteg_dgdiffusion_kernels.hpp"
+
+TEST_CASE("DGDiffusion Kernel Specializations", "[Specializations]")
+{
+   using namespace mfem;
+
+   DGDiffusionIntegrator::AddSpecialization<2, 2, 4>();
+}

tests/unit/fem/specializations/test_dgmassinv.cpp

@@ -0,0 +1,14 @@
+/// Tests which make sure adding user-defined kernel specializations work.
+/// These tests are compile/link-only tests
+
+#include "mfem.hpp"
+#include "unit_tests.hpp"
+
+#include "fem/dgmassinv_kernels.hpp"
+
+TEST_CASE("DGMassInverse Kernel Specializations", "[Specializations]")
+{
+   using namespace mfem;
+
+   DGMassInverse::CGKernels::Specialization<2, 1, 2>::Add();
+}

tests/unit/fem/specializations/test_dgtrace_integ.cpp

@@ -0,0 +1,14 @@
+/// Tests which make sure adding user-defined kernel specializations work.
+/// These tests are compile/link-only tests
+
+#include "mfem.hpp"
+#include "unit_tests.hpp"
+
+#include "fem/integ/bilininteg_dgtrace_kernels.hpp"
+
+TEST_CASE("DGTrace Kernel Specializations", "[Specializations]")
+{
+   using namespace mfem;
+
+   DGTraceIntegrator::AddSpecialization<2, 2, 3>();
+}

tests/unit/fem/specializations/test_diffusion_integ.cpp

@@ -0,0 +1,15 @@
+/// Tests which make sure adding user-defined kernel specializations work.
+/// These tests are compile/link-only tests
+
+#include "mfem.hpp"
+#include "unit_tests.hpp"
+
+#include "fem/integ/bilininteg_diffusion_kernels.hpp"
+
+TEST_CASE("Diffusion Kernel Specializations", "[Specializations]")
+{
+   using namespace mfem;
+
+   DiffusionIntegrator::AddSpecialization<2, 1, 5>();
+   DiffusionIntegrator::AddSimplexSpecialization<2, 2, 3>();
+}

tests/unit/fem/specializations/test_mass_integ.cpp

@@ -0,0 +1,15 @@
+/// Tests which make sure adding user-defined kernel specializations work.
+/// These tests are compile/link-only tests
+
+#include "mfem.hpp"
+#include "unit_tests.hpp"
+
+#include "fem/integ/bilininteg_mass_kernels.hpp"
+
+TEST_CASE("Mass Kernel Specializations", "[Specializations]")
+{
+   using namespace mfem;
+
+   MassIntegrator::AddSpecialization<2, 1, 3>();
+   MassIntegrator::AddSimplexSpecialization<2, 2, 4>();
+}

tests/unit/fem/specializations/test_qinterp_det.cpp

@@ -0,0 +1,14 @@
+/// Tests which make sure adding user-defined kernel specializations work.
+/// These tests are compile/link-only tests
+
+#include "mfem.hpp"
+#include "unit_tests.hpp"
+
+#include "fem/qinterp/det.hpp"
+
+TEST_CASE("QInterp Det Kernel Specializations", "[Specializations]")
+{
+   using namespace mfem;
+
+   QuadratureInterpolator::AddDetSpecializations<2, 3, 2, 2>();
+}

tests/unit/fem/specializations/test_qinterp_eval.cpp

@@ -0,0 +1,14 @@
+/// Tests which make sure adding user-defined kernel specializations work.
+/// These tests are compile/link-only tests
+
+#include "mfem.hpp"
+#include "unit_tests.hpp"
+
+#include "fem/qinterp/eval.hpp"
+
+TEST_CASE("QInterp Eval Kernel Specializations", "[Specializations]")
+{
+   using namespace mfem;
+
+   QuadratureInterpolator::AddEvalSpecializations<2, 1, 1, 2>();
+}

tests/unit/fem/specializations/test_qinterp_eval_hdiv.cpp

@@ -0,0 +1,19 @@
+/// Tests which make sure adding user-defined kernel specializations work.
+/// These tests are compile/link-only tests
+
+#include "mfem.hpp"
+#include "unit_tests.hpp"
+
+#include "fem/qinterp/eval_hdiv.hpp"
+
+TEST_CASE("QInterp Eval HDiv Kernel Specializations", "[Specializations]")
+{
+   using namespace mfem;
+
+   QuadratureInterpolator::TensorEvalHDivKernels::Specialization<
+   2, QVectorLayout::byNODES, QuadratureInterpolator::PHYSICAL_VALUES, 2,
+   4>::Add();
+   QuadratureInterpolator::TensorEvalHDivKernels::Specialization<
+   2, QVectorLayout::byNODES, QuadratureInterpolator::PHYSICAL_MAGNITUDES, 2,
+   4>::Add();
+}

tests/unit/fem/specializations/test_qinterp_grad.cpp

@@ -0,0 +1,21 @@
+/// Tests which make sure adding user-defined kernel specializations work.
+/// These tests are compile/link-only tests
+
+#include "mfem.hpp"
+#include "unit_tests.hpp"
+
+#include "fem/qinterp/grad.hpp"
+
+TEST_CASE("QInterp Grad Kernel Specializations", "[Specializations]")
+{
+   using namespace mfem;
+
+   QuadratureInterpolator::AddGradSpecializations<
+   2, QVectorLayout::byNODES, false, 1, 3, 3, 1>();
+   QuadratureInterpolator::AddGradSpecializations<
+   2, QVectorLayout::byNODES, true, 1, 3, 3, 1>();
+   QuadratureInterpolator::AddCollocatedGradSpecializations<
+   2, QVectorLayout::byNODES, false, 1, 2, 1>();
+   QuadratureInterpolator::AddCollocatedGradSpecializations<
+   2, QVectorLayout::byNODES, true, 1, 2, 1>();
+}

tests/unit/fem/specializations/test_qinterp_tensoreval.cpp

@@ -0,0 +1,15 @@
+/// Tests which make sure adding user-defined kernel specializations work.
+/// These tests are compile/link-only tests
+
+#include "mfem.hpp"
+#include "unit_tests.hpp"
+
+#include "fem/qinterp/eval.hpp"
+
+TEST_CASE("QInterp TensorEval Kernel Specializations", "[Specializations]")
+{
+   using namespace mfem;
+
+   QuadratureInterpolator::AddTensorEvalSpecializations<
+   2, QVectorLayout::byNODES, 1, 3, 3, 2>();
+}

tests/unit/fem/specializations/test_vecdiffusion_integ.cpp

@@ -0,0 +1,14 @@
+/// Tests which make sure adding user-defined kernel specializations work.
+/// These tests are compile/link-only tests
+
+#include "mfem.hpp"
+#include "unit_tests.hpp"
+
+#include "fem/integ/bilininteg_vecdiffusion_pa.hpp"
+
+TEST_CASE("VectorDiffusion Kernel Specializations", "[Specializations]")
+{
+   using namespace mfem;
+
+   VectorDiffusionIntegrator::AddSpecialization<2, 2, 2, 3>();
+}

tests/unit/fem/specializations/test_vecmass_integ.cpp

@@ -0,0 +1,14 @@
+/// Tests which make sure adding user-defined kernel specializations work.
+/// These tests are compile/link-only tests
+
+#include "mfem.hpp"
+#include "unit_tests.hpp"
+
+#include "fem/integ/bilininteg_vecmass_pa.hpp"
+
+TEST_CASE("Vector Mass Kernel Specializations", "[Specializations]")
+{
+   using namespace mfem;
+
+   VectorMassIntegrator::VectorMassAddMultPA::Specialization<2, 2, 4>::Add();
+}

tests/unit/fem/test_mixedsesqform.cpp

@@ -0,0 +1,323 @@
+// 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 "mfem.hpp"
+#include "unit_tests.hpp"
+
+using namespace mfem;
+
+namespace
+{
+constexpr real_t a_coef = 1.0;
+constexpr real_t b_coef = 2.0;
+constexpr real_t c_coef = 3.0;
+constexpr real_t omega_val = 10.0;
+
+real_t V_exact_fn(const Vector &x)
+{
+   return a_coef*x[0] + b_coef*x[1] + c_coef*x[2];
+}
+} // namespace
+
+TEST_CASE("Mixed Sesquilinear Form", "[MixedSesquilinearForm]")
+{
+   const bool cross = GENERATE(false, true);
+   const auto conv = GENERATE(ComplexOperator::HERMITIAN,
+                              ComplexOperator::BLOCK_SYMMETRIC);
+   CAPTURE(cross, int(conv));
+
+   Mesh mesh = Mesh::MakeCartesian3D(10, 10, 1, Element::HEXAHEDRON);
+   H1_FECollection fec_h1(1, mesh.Dimension());
+   FiniteElementSpace fespace_h1(&mesh, &fec_h1);
+   ND_FECollection fec_nd(1, mesh.Dimension());
+   FiniteElementSpace fespace_nd(&mesh, &fec_nd);
+
+   Array<int> dbc_bdr(mesh.bdr_attributes.Max());
+   dbc_bdr = 1;
+   Array<int> ess_tdof_list_h1;
+   Array<int> ess_tdof_list_nd;
+   fespace_h1.GetEssentialTrueDofs(dbc_bdr, ess_tdof_list_h1);
+   fespace_nd.GetEssentialTrueDofs(dbc_bdr, ess_tdof_list_nd);
+
+   ConstantCoefficient omega(omega_val);
+   ConstantCoefficient neg_omega(-omega_val);
+   ConstantCoefficient half(0.5);
+   FunctionCoefficient V_exact_real(V_exact_fn);
+   const real_t den = cross ? 2.0*omega_val : omega_val;
+   Vector A_vec({a_coef/den, b_coef/den, c_coef/den});
+   VectorConstantCoefficient A_exact_imag(A_vec);
+   A_vec *= (cross ? -1.0 : 0.0);
+   VectorConstantCoefficient A_exact_real(A_vec);
+
+   ComplexGridFunction V(&fespace_h1);
+   ComplexGridFunction A(&fespace_nd);
+
+   V = 0.0;
+   A = 0.0;
+   V.real().ProjectBdrCoefficient(V_exact_real, dbc_bdr);
+   A.real().ProjectBdrCoefficientTangent(A_exact_real, dbc_bdr);
+   A.imag().ProjectBdrCoefficientTangent(A_exact_imag, dbc_bdr);
+
+   ComplexLinearForm b_h1(&fespace_h1);
+   b_h1 = 0.0;
+   b_h1.Assemble();
+   ComplexLinearForm b_nd(&fespace_nd);
+   b_nd = 0.0;
+   b_nd.Assemble();
+
+   // Add integrators to the blocks
+   MixedSesquilinearForm a_h1_nd(&fespace_h1, &fespace_nd, conv);
+   a_h1_nd.AddDomainIntegrator(cross ? new MixedVectorGradientIntegrator(half)
+                               : new MixedVectorGradientIntegrator,
+                               cross ? new MixedVectorGradientIntegrator(half)
+                               : nullptr);
+   a_h1_nd.Assemble();
+
+   MixedSesquilinearForm a_nd_h1(&fespace_nd, &fespace_h1, conv);
+   a_nd_h1.AddDomainIntegrator(
+      cross ? new MixedVectorWeakDivergenceIntegrator(neg_omega) : nullptr,
+      new MixedVectorWeakDivergenceIntegrator(neg_omega));
+   a_nd_h1.Assemble();
+
+   SesquilinearForm a_h1(&fespace_h1, conv);
+   a_h1.AddDomainIntegrator(new DiffusionIntegrator, nullptr);
+   a_h1.Assemble();
+
+   SesquilinearForm a_nd(&fespace_nd, conv);
+   a_nd.AddDomainIntegrator(new CurlCurlIntegrator, nullptr);
+   a_nd.AddDomainIntegrator(nullptr, new VectorFEMassIntegrator(omega));
+   a_nd.Assemble();
+
+   // Set block offsets (doubled for real+imag)
+   mfem::Array<int> bOffsets(3);
+   bOffsets[0] = 0;
+   bOffsets[1] = 2 * fespace_h1.GetTrueVSize();
+   bOffsets[2] = 2 * fespace_nd.GetTrueVSize();
+   bOffsets.PartialSum();
+
+   OperatorPtr A_h1, A_nd, A_h1_nd, A_nd_h1;
+   BlockVector trueX(bOffsets), trueRHS(bOffsets);
+   Vector B_h1, X_h1, B_nd, X_nd;
+
+   trueX = 0.0;
+   trueRHS = 0.0;
+
+   // Form the diagonal entries
+   a_h1.FormLinearSystem(ess_tdof_list_h1, V, b_h1, A_h1, X_h1, B_h1);
+   a_nd.FormLinearSystem(ess_tdof_list_nd, A, b_nd, A_nd, X_nd, B_nd);
+
+   trueX.GetBlock(0) = X_h1;
+   trueX.GetBlock(1) = X_nd;
+   trueRHS.GetBlock(0) += B_h1;
+   trueRHS.GetBlock(1) += B_nd;
+
+   // Form the off-diagonal entries
+   a_h1_nd.FormRectangularLinearSystem(ess_tdof_list_h1, ess_tdof_list_nd, V, b_nd,
+                                       A_h1_nd, X_h1, B_nd);
+   a_nd_h1.FormRectangularLinearSystem(ess_tdof_list_nd, ess_tdof_list_h1, A, b_h1,
+                                       A_nd_h1, X_nd, B_h1);
+
+   trueRHS.GetBlock(0) += B_h1;
+   trueRHS.GetBlock(1) += B_nd;
+
+   auto *Ah1 = A_h1.As<ComplexSparseMatrix>();
+   auto *And = A_nd.As<ComplexSparseMatrix>();
+   auto *Ah1nd = A_h1_nd.As<ComplexSparseMatrix>();
+   auto *Andh1 = A_nd_h1.As<ComplexSparseMatrix>();
+
+   BlockOperator blockOp(bOffsets);
+   blockOp.SetBlock(0, 0, Ah1);
+   blockOp.SetBlock(1, 1, And);
+   blockOp.SetBlock(0, 1, Andh1);
+   blockOp.SetBlock(1, 0, Ah1nd);
+
+   SparseMatrix *Sh1 = Ah1->GetSystemMatrix();
+   SparseMatrix *Snd = And->GetSystemMatrix();
+   GSSmoother smoothSh1(*Sh1), smoothSnd(*Snd);
+
+   BlockDiagonalPreconditioner P(bOffsets);
+   P.SetDiagonalBlock(0, &smoothSh1);
+   P.SetDiagonalBlock(1, &smoothSnd);
+
+   GMRESSolver gmres;
+   gmres.SetOperator(blockOp);
+   gmres.SetPreconditioner(P);
+   gmres.SetAbsTol(1e-10);
+   gmres.SetMaxIter(2000);
+   gmres.SetKDim(200);
+   gmres.SetPrintLevel(2);
+   gmres.Mult(trueRHS, trueX);
+
+   delete Sh1;
+   delete Snd;
+
+   V = trueX.GetBlock(0);
+   A = trueX.GetBlock(1);
+
+   // Check solution
+   ConstantCoefficient zero(0.0);
+
+   real_t err_V = V.ComputeL2Error(V_exact_real, zero);
+   real_t err_A = A.ComputeL2Error(A_exact_real, A_exact_imag);
+
+   REQUIRE(err_V == MFEM_Approx(0.0, 1e-5));
+   REQUIRE(err_A == MFEM_Approx(0.0, 1e-5));
+}
+
+#ifdef MFEM_USE_MPI
+#ifdef MFEM_USE_SUPERLU
+
+TEST_CASE("Parallel Mixed Sesquilinear Form",
+          "[MixedSesquilinearForm][Parallel]")
+{
+   // See the serial test above for the manufactured solution
+   const bool cross = GENERATE(false, true);
+   const auto conv = GENERATE(ComplexOperator::HERMITIAN,
+                              ComplexOperator::BLOCK_SYMMETRIC);
+   CAPTURE(cross, int(conv));
+
+   Mesh mesh = Mesh::MakeCartesian3D(10, 10, 1, Element::HEXAHEDRON);
+   ParMesh par_mesh(MPI_COMM_WORLD, mesh);
+   H1_FECollection fec_h1(1, mesh.Dimension());
+   ParFiniteElementSpace fespace_h1(&par_mesh, &fec_h1);
+   ND_FECollection fec_nd(1, mesh.Dimension());
+   ParFiniteElementSpace fespace_nd(&par_mesh, &fec_nd);
+
+   Array<int> dbc_bdr(par_mesh.bdr_attributes.Max());
+   dbc_bdr = 1;
+   Array<int> ess_tdof_list_h1;
+   Array<int> ess_tdof_list_nd;
+   fespace_h1.GetEssentialTrueDofs(dbc_bdr, ess_tdof_list_h1);
+   fespace_nd.GetEssentialTrueDofs(dbc_bdr, ess_tdof_list_nd);
+
+   ConstantCoefficient omega(omega_val);
+   ConstantCoefficient neg_omega(-omega_val);
+   ConstantCoefficient half(0.5);
+   FunctionCoefficient V_exact_real(V_exact_fn);
+   const real_t den = cross ? 2.0*omega_val : omega_val;
+   Vector A_vec({a_coef/den, b_coef/den, c_coef/den});
+   VectorConstantCoefficient A_exact_imag(A_vec);
+   A_vec *= (cross ? -1.0 : 0.0);
+   VectorConstantCoefficient A_exact_real(A_vec);
+
+   ParComplexGridFunction V(&fespace_h1);
+   ParComplexGridFunction A(&fespace_nd);
+
+   V = 0.0;
+   A = 0.0;
+   V.real().ProjectBdrCoefficient(V_exact_real, dbc_bdr);
+   A.real().ProjectBdrCoefficientTangent(A_exact_real, dbc_bdr);
+   A.imag().ProjectBdrCoefficientTangent(A_exact_imag, dbc_bdr);
+
+   ParComplexLinearForm b_h1(&fespace_h1);
+   b_h1 = 0.0;
+   b_h1.Assemble();
+   ParComplexLinearForm b_nd(&fespace_nd);
+   b_nd = 0.0;
+   b_nd.Assemble();
+
+   // Add integrators to the blocks
+   ParMixedSesquilinearForm a_h1_nd(&fespace_h1, &fespace_nd, conv);
+   a_h1_nd.AddDomainIntegrator(cross ? new MixedVectorGradientIntegrator(half)
+                               : new MixedVectorGradientIntegrator,
+                               cross ? new MixedVectorGradientIntegrator(half)
+                               : nullptr);
+   a_h1_nd.Assemble();
+
+   ParMixedSesquilinearForm a_nd_h1(&fespace_nd, &fespace_h1, conv);
+   a_nd_h1.AddDomainIntegrator(
+      cross ? new MixedVectorWeakDivergenceIntegrator(neg_omega) : nullptr,
+      new MixedVectorWeakDivergenceIntegrator(neg_omega));
+   a_nd_h1.Assemble();
+
+   ParSesquilinearForm a_h1(&fespace_h1, conv);
+   a_h1.AddDomainIntegrator(new DiffusionIntegrator, nullptr);
+   a_h1.Assemble();
+
+   ParSesquilinearForm a_nd(&fespace_nd, conv);
+   a_nd.AddDomainIntegrator(new CurlCurlIntegrator, nullptr);
+   a_nd.AddDomainIntegrator(nullptr, new VectorFEMassIntegrator(omega));
+   a_nd.Assemble();
+
+   mfem::Array2D<const mfem::HypreParMatrix *> h_blocks;
+   h_blocks.SetSize(2, 2);
+   h_blocks = nullptr;
+
+   // Set block offsets
+   mfem::Array<int> bOffsets(3);
+   bOffsets[0] = 0;
+   bOffsets[1] = 2 * fespace_h1.TrueVSize();
+   bOffsets[2] = 2 * fespace_nd.TrueVSize();
+   bOffsets.PartialSum();
+
+   OperatorPtr A_h1, A_nd, A_h1_nd, A_nd_h1;
+   BlockVector trueX(bOffsets), trueRHS(bOffsets);
+   Vector B_h1, X_h1, B_nd, X_nd;
+
+   trueX = 0.0;
+   trueRHS = 0.0;
+
+   // Form the diagonal entries
+   a_h1.FormLinearSystem(ess_tdof_list_h1, V, b_h1, A_h1, X_h1, B_h1);
+   a_nd.FormLinearSystem(ess_tdof_list_nd, A, b_nd, A_nd, X_nd, B_nd);
+
+   trueX.GetBlock(0) = X_h1;
+   trueX.GetBlock(1) = X_nd;
+   trueRHS.GetBlock(0) += B_h1;
+   trueRHS.GetBlock(1) += B_nd;
+
+   // Form the off-diagonal entries
+   a_h1_nd.FormRectangularLinearSystem(ess_tdof_list_h1, ess_tdof_list_nd, V, b_nd,
+                                       A_h1_nd, X_h1, B_nd);
+   a_nd_h1.FormRectangularLinearSystem(ess_tdof_list_nd, ess_tdof_list_h1, A, b_h1,
+                                       A_nd_h1, X_nd, B_h1);
+
+   trueRHS.GetBlock(0) += B_h1;
+   trueRHS.GetBlock(1) += B_nd;
+
+   h_blocks(0,0) = A_h1.As<ComplexHypreParMatrix>()->GetSystemMatrix();
+   h_blocks(1,1) = A_nd.As<ComplexHypreParMatrix>()->GetSystemMatrix();
+   h_blocks(0,1) = A_nd_h1.As<ComplexHypreParMatrix>()->GetSystemMatrix();
+   h_blocks(1,0) = A_h1_nd.As<ComplexHypreParMatrix>()->GetSystemMatrix();
+
+   OperatorHandle op(HypreParMatrixFromBlocks(h_blocks));
+
+   SuperLURowLocMatrix S_op(*op);
+   SuperLUSolver superlu(MPI_COMM_WORLD);
+   superlu.SetPrintStatistics(false);
+   superlu.SetSymmetricPattern(false);
+   superlu.SetOperator(S_op);
+   superlu.Mult(trueRHS, trueX);
+
+   trueX.GetBlock(0).SyncAliasMemory(trueX);
+   trueX.GetBlock(1).SyncAliasMemory(trueX);
+
+   V.Distribute(trueX.GetBlock(0));
+   A.Distribute(trueX.GetBlock(1));
+
+   // Check solution
+   ConstantCoefficient zero(0.0);
+
+   real_t err_Vr = V.real().ComputeL2Error(V_exact_real);
+   real_t err_Vi = V.imag().ComputeL2Error(zero);
+   real_t err_Ar = A.real().ComputeL2Error(A_exact_real);
+   real_t err_Ai = A.imag().ComputeL2Error(A_exact_imag);
+
+   REQUIRE(err_Vr == MFEM_Approx(0.0, 1e-5));
+   REQUIRE(err_Vi == MFEM_Approx(0.0, 1e-5));
+   REQUIRE(err_Ar == MFEM_Approx(0.0, 1e-5));
+   REQUIRE(err_Ai == MFEM_Approx(0.0, 1e-5));
+}
+
+#endif
+#endif
+
+

tests/unit/fem/test_pa_coeff.cpp

@@ -750,4 +750,313 @@ TEST_CASE("Hcurl/Hdiv Mixed PA Coefficient",
    }
 }
 
+TEST_CASE("3D Bilinear VectorFE Integrators PartialAssembly",
+          "[BilinearFormIntegrator]"
+          "[PartialAssembly]"
+          "[GPU]")
+{
+   auto order = GENERATE(1, 2);
+   CAPTURE(order);
+   dimension = 3;
+
+   FunctionCoefficient q3_coeff(coeffFunction);
+   VectorFunctionCoefficient F3_coeff(dimension, vectorCoeffFunction);
+   MatrixFunctionCoefficient M3_coeff(dimension, asymmetricMatrixCoeffFunction);
+
+   auto mesh_fname =
+      GENERATE("../../data/fichera-amr.mesh", "../../data/fichera-q2.mesh");
+   CAPTURE(mesh_fname);
+   Mesh mesh(mesh_fname);
+   REQUIRE(mesh.Dimension() == dimension);
+   REQUIRE(mesh.SpaceDimension() == dimension);
+
+   SECTION("RT to RT Scalar Coeff")
+   {
+      RT_FECollection fec_rt(order - 1, dimension);
+      FiniteElementSpace fespace_rt(&mesh, &fec_rt);
+
+      BilinearForm bfa(&fespace_rt);
+      bfa.AddDomainIntegrator(new VectorFEMassIntegrator(q3_coeff));
+      bfa.Assemble();
+      bfa.Finalize();
+
+      BilinearForm bpa(&fespace_rt);
+      bpa.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+      bpa.AddDomainIntegrator(new VectorFEMassIntegrator(q3_coeff));
+      bpa.Assemble();
+
+      GridFunction x(&fespace_rt), y_fa(&fespace_rt), y_pa(&fespace_rt);
+      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 RT Diagonal Matrix Coeff")
+   {
+      RT_FECollection fec_rt(order - 1, dimension);
+      FiniteElementSpace fespace_rt(&mesh, &fec_rt);
+
+      BilinearForm bfa(&fespace_rt);
+      bfa.AddDomainIntegrator(new VectorFEMassIntegrator(F3_coeff));
+      bfa.Assemble();
+      bfa.Finalize();
+
+      BilinearForm bpa(&fespace_rt);
+      bpa.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+      bpa.AddDomainIntegrator(new VectorFEMassIntegrator(F3_coeff));
+      bpa.Assemble();
+
+      GridFunction x(&fespace_rt), y_fa(&fespace_rt), y_pa(&fespace_rt);
+      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 RT Matrix Coeff")
+   {
+      RT_FECollection fec_rt(order - 1, dimension);
+      FiniteElementSpace fespace_rt(&mesh, &fec_rt);
+
+      BilinearForm bfa(&fespace_rt);
+      bfa.AddDomainIntegrator(new VectorFEMassIntegrator(M3_coeff));
+      bfa.Assemble();
+      bfa.Finalize();
+
+      BilinearForm bpa(&fespace_rt);
+      bpa.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+      bpa.AddDomainIntegrator(new VectorFEMassIntegrator(M3_coeff));
+      bpa.Assemble();
+
+      GridFunction x(&fespace_rt), y_fa(&fespace_rt), y_pa(&fespace_rt);
+      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 Scalar Coeff")
+   {
+      RT_FECollection fec_rt(order - 1, dimension);
+      FiniteElementSpace fespace_rt(&mesh, &fec_rt);
+      ND_FECollection fec_nd(order, dimension);
+      FiniteElementSpace fespace_nd(&mesh, &fec_nd);
+
+      MixedBilinearForm bfa(&fespace_rt, &fespace_nd);
+      bfa.AddDomainIntegrator(new VectorFEMassIntegrator(q3_coeff));
+      bfa.Assemble();
+      bfa.Finalize();
+
+      MixedBilinearForm bpa(&fespace_rt, &fespace_nd);
+      bpa.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+      bpa.AddDomainIntegrator(new VectorFEMassIntegrator(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")
+   {
+      RT_FECollection fec_rt(order - 1, dimension);
+      FiniteElementSpace fespace_rt(&mesh, &fec_rt);
+      ND_FECollection fec_nd(order, dimension);
+      FiniteElementSpace fespace_nd(&mesh, &fec_nd);
+
+      MixedBilinearForm bfa(&fespace_rt, &fespace_nd);
+      bfa.AddDomainIntegrator(new VectorFEMassIntegrator(F3_coeff));
+      bfa.Assemble();
+      bfa.Finalize();
+
+      MixedBilinearForm bpa(&fespace_rt, &fespace_nd);
+      bpa.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+      bpa.AddDomainIntegrator(new VectorFEMassIntegrator(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) );
+   }
+
+   SECTION("RT to ND Matrix Coeff")
+   {
+      RT_FECollection fec_rt(order - 1, dimension);
+      FiniteElementSpace fespace_rt(&mesh, &fec_rt);
+      ND_FECollection fec_nd(order, dimension);
+      FiniteElementSpace fespace_nd(&mesh, &fec_nd);
+
+      MixedBilinearForm bfa(&fespace_rt, &fespace_nd);
+      bfa.AddDomainIntegrator(new VectorFEMassIntegrator(M3_coeff));
+      bfa.Assemble();
+      bfa.Finalize();
+
+      MixedBilinearForm bpa(&fespace_rt, &fespace_nd);
+      bpa.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+      bpa.AddDomainIntegrator(new VectorFEMassIntegrator(M3_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("ND to RT Scalar Coeff")
+   {
+      RT_FECollection fec_rt(order - 1, dimension);
+      FiniteElementSpace fespace_rt(&mesh, &fec_rt);
+      ND_FECollection fec_nd(order, dimension);
+      FiniteElementSpace fespace_nd(&mesh, &fec_nd);
+
+      MixedBilinearForm bfa(&fespace_nd, &fespace_rt);
+      bfa.AddDomainIntegrator(new VectorFEMassIntegrator(q3_coeff));
+      bfa.Assemble();
+      bfa.Finalize();
+
+      MixedBilinearForm bpa(&fespace_nd, &fespace_rt);
+      bpa.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+      bpa.AddDomainIntegrator(new VectorFEMassIntegrator(q3_coeff));
+      bpa.Assemble();
+
+      GridFunction x(&fespace_nd), y_fa(&fespace_rt), y_pa(&fespace_rt);
+      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("ND to RT Diagonal Matrix Coeff")
+   {
+      RT_FECollection fec_rt(order - 1, dimension);
+      FiniteElementSpace fespace_rt(&mesh, &fec_rt);
+      ND_FECollection fec_nd(order, dimension);
+      FiniteElementSpace fespace_nd(&mesh, &fec_nd);
+
+      MixedBilinearForm bfa(&fespace_nd, &fespace_rt);
+      bfa.AddDomainIntegrator(new VectorFEMassIntegrator(F3_coeff));
+      bfa.Assemble();
+      bfa.Finalize();
+
+      MixedBilinearForm bpa(&fespace_nd, &fespace_rt);
+      bpa.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+      bpa.AddDomainIntegrator(new VectorFEMassIntegrator(F3_coeff));
+      bpa.Assemble();
+
+      GridFunction x(&fespace_nd), y_fa(&fespace_rt), y_pa(&fespace_rt);
+      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("ND to RT Matrix Coeff")
+   {
+      RT_FECollection fec_rt(order - 1, dimension);
+      FiniteElementSpace fespace_rt(&mesh, &fec_rt);
+      ND_FECollection fec_nd(order, dimension);
+      FiniteElementSpace fespace_nd(&mesh, &fec_nd);
+
+      MixedBilinearForm bfa(&fespace_nd, &fespace_rt);
+      bfa.AddDomainIntegrator(new VectorFEMassIntegrator(M3_coeff));
+      bfa.Assemble();
+      bfa.Finalize();
+
+      MixedBilinearForm bpa(&fespace_nd, &fespace_rt);
+      bpa.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+      bpa.AddDomainIntegrator(new VectorFEMassIntegrator(M3_coeff));
+      bpa.Assemble();
+
+      GridFunction x(&fespace_nd), y_fa(&fespace_rt), y_pa(&fespace_rt);
+      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("ND to ND Scalar Coeff")
+   {
+      ND_FECollection fec_nd(order, dimension);
+      FiniteElementSpace fespace_nd(&mesh, &fec_nd);
+
+      BilinearForm bfa(&fespace_nd);
+      bfa.AddDomainIntegrator(new VectorFEMassIntegrator(q3_coeff));
+      bfa.Assemble();
+      bfa.Finalize();
+
+      BilinearForm bpa(&fespace_nd);
+      bpa.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+      bpa.AddDomainIntegrator(new VectorFEMassIntegrator(q3_coeff));
+      bpa.Assemble();
+
+      GridFunction x(&fespace_nd), 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("ND to ND Diagonal Matrix Coeff")
+   {
+      ND_FECollection fec_nd(order, dimension);
+      FiniteElementSpace fespace_nd(&mesh, &fec_nd);
+
+      BilinearForm bfa(&fespace_nd);
+      bfa.AddDomainIntegrator(new VectorFEMassIntegrator(F3_coeff));
+      bfa.Assemble();
+      bfa.Finalize();
+
+      BilinearForm bpa(&fespace_nd);
+      bpa.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+      bpa.AddDomainIntegrator(new VectorFEMassIntegrator(F3_coeff));
+      bpa.Assemble();
+
+      GridFunction x(&fespace_nd), 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("ND to ND Matrix Coeff")
+   {
+      ND_FECollection fec_nd(order, dimension);
+      FiniteElementSpace fespace_nd(&mesh, &fec_nd);
+
+      BilinearForm bfa(&fespace_nd);
+      bfa.AddDomainIntegrator(new VectorFEMassIntegrator(M3_coeff));
+      bfa.Assemble();
+      bfa.Finalize();
+
+      BilinearForm bpa(&fespace_nd);
+      bpa.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+      bpa.AddDomainIntegrator(new VectorFEMassIntegrator(M3_coeff));
+      bpa.Assemble();
+
+      GridFunction x(&fespace_nd), 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/fem/test_project_bdr.cpp

@@ -17,12 +17,63 @@ using namespace mfem;
 namespace project_bdr
 {
 
-void Func_3D_lin(const Vector &x, Vector &v)
+TEST_CASE("3D ProjectBdrCoefficient",
+          "[GridFunction]"
+          "[NCMesh]")
 {
-   v.SetSize(3);
-   v[0] =  1.234 * x[0] - 2.357 * x[1] + 3.572 * x[2];
-   v[1] =  2.537 * x[0] + 4.321 * x[1] - 1.234 * x[2];
-   v[2] = -2.572 * x[0] + 1.321 * x[1] + 3.234 * x[2];
+   const char *mesh_file = GENERATE("data/hex-nc-cross.mesh",
+                                    "data/tet-nc-cross.mesh");
+   CAPTURE(mesh_file);
+   constexpr int order = 3;
+   constexpr real_t freq = 5.0;
+
+   // Attributes
+   Array<int> bdr_attr(2);
+   bdr_attr = 0;
+   bdr_attr[1] = 1;
+
+   // Coefficient
+   FunctionCoefficient coeff([&](const Vector &x)
+   {
+      return cos(freq * M_PI * x[0])
+             * cos(freq * M_PI * x[1])
+             * cos(freq * M_PI * x[2]);
+   });
+
+   // Vertex-based mesh
+   Mesh mesh_v(mesh_file, 1, 1);
+   H1_FECollection fec_v(order, mesh_v.Dimension());
+   FiniteElementSpace fes_v(&mesh_v, &fec_v);
+   GridFunction gf_v(&fes_v);
+   gf_v = 0.0;
+   gf_v.ProjectBdrCoefficient(coeff, bdr_attr);
+
+   // Nodal mesh
+   Mesh mesh_n(mesh_file, 1, 1);
+   mesh_n.SetCurvature(order, true);
+   H1_FECollection fec_n(order, mesh_n.Dimension());
+   FiniteElementSpace fes_n(&mesh_n, &fec_n);
+   GridFunction gf_n(&fes_n);
+   gf_n = 0.0;
+   gf_n.ProjectBdrCoefficient(coeff, bdr_attr);
+
+   gf_n -= gf_v;
+
+   REQUIRE(gf_n.Norml2() == MFEM_Approx(0.0));
+}
+
+void Func_lin(const Vector &x, Vector &v)
+{
+   const int dim = x.Size();
+   v.SetSize(dim);
+   v[0] = 1.234 * x[0] - 2.357 * x[1];
+   v[1] = 2.537 * x[0] + 4.321 * x[1];
+   if (dim == 3)
+   {
+      v[0] += 3.572 * x[2];
+      v[1] -= 1.234 * x[2];
+      v[2] = -2.572 * x[0] + 1.321 * x[1] + 3.234 * x[2];
+   }
 }
 
 TEST_CASE("3D ProjectBdrCoefficientNormal Vector",
@@ -41,7 +92,7 @@ TEST_CASE("3D ProjectBdrCoefficientNormal Vector",
       Mesh mesh = Mesh::MakeCartesian3D(
                      n, n, n, (Element::Type)type, 2.0, 3.0, 5.0);
 
-      VectorFunctionCoefficient funcCoef(dim, Func_3D_lin);
+      VectorFunctionCoefficient funcCoef(dim, Func_lin);
 
       SECTION("3D GetVectorValue tests for element type " +
               std::to_string(type))
@@ -133,7 +184,7 @@ TEST_CASE("3D ProjectBdrCoefficientNormal Scalar",
       Mesh mesh = Mesh::MakeCartesian3D(
                      n, n, n, (Element::Type)type, 2.0, 3.0, 5.0);
 
-      VectorFunctionCoefficient funcCoef(dim, Func_3D_lin);
+      VectorFunctionCoefficient funcCoef(dim, Func_lin);
 
       SECTION("3D GetVectorValue tests for element type " +
               std::to_string(type))
@@ -227,7 +278,7 @@ TEST_CASE("3D ProjectBdrCoefficientTangent",
       Mesh mesh = Mesh::MakeCartesian3D(
                      n, n, n, (Element::Type)type, 2.0, 3.0, 5.0);
 
-      VectorFunctionCoefficient funcCoef(dim, Func_3D_lin);
+      VectorFunctionCoefficient funcCoef(dim, Func_lin);
 
       SECTION("3D GetVectorValue tests for element type " +
               std::to_string(type))
@@ -305,4 +356,50 @@ TEST_CASE("3D ProjectBdrCoefficientTangent",
    }
 }
 
+TEST_CASE("ProjectBdrCoefficientTangent with IntegratedGLL",
+          "[GridFunction]"
+          "[VectorGridFunctionCoefficient]")
+{
+   const int dim = GENERATE(2, 3);
+   CAPTURE(dim);
+   Mesh mesh = (dim == 2) ?
+               Mesh::MakeCartesian2D(1, 2, Element::QUADRILATERAL,
+                                     true, 2.0, 5.0) :
+               Mesh::MakeCartesian3D(1, 1, 2, Element::HEXAHEDRON,
+                                     2.0, 3.0, 5.0);
+   mesh.EnsureNodes();
+   mesh.EnsureNCMesh(false);
+   VectorFunctionCoefficient func_coef(dim, Func_lin);
+   Array<int> all_bdr(mesh.bdr_attributes.Max());
+   all_bdr = 1;
+
+   for (int order = 1; order <= 4; order++)
+   {
+      CAPTURE(order);
+      ND_FECollection nd_fec(order, dim, BasisType::GaussLobatto,
+                             BasisType::IntegratedGLL);
+      FiniteElementSpace nd_fespace(&mesh, &nd_fec);
+      GridFunction volume_projection(&nd_fespace);
+      GridFunction boundary_projection(&nd_fespace);
+
+      volume_projection.ProjectCoefficient(func_coef);
+      boundary_projection = 0.0;
+      boundary_projection.ProjectBdrCoefficientTangent(func_coef, all_bdr);
+
+      Array<int> ess_vdofs;
+      nd_fespace.GetEssentialVDofs(all_bdr, ess_vdofs);
+      real_t max_error = 0.0;
+      for (int i = 0; i < ess_vdofs.Size(); i++)
+      {
+         if (ess_vdofs[i])
+         {
+            max_error = std::max(max_error, std::abs(
+                                    boundary_projection[i] -
+                                    volume_projection[i]));
+         }
+      }
+      REQUIRE(max_error == MFEM_Approx(0.0));
+   }
+}
+
 } // namespace project_bdr

tests/unit/linalg/test_matrix_dense.cpp

@@ -494,6 +494,64 @@ TEST_CASE("Batched Linear Algebra",
    }
 }
 
+#ifdef MFEM_USE_EXCEPTIONS
+namespace
+{
+
+DenseTensor MakeSingularBatchedMatrices()
+{
+   const int n = 3;
+   const int n_mat = 3;
+
+   DenseTensor A_batch(n, n, n_mat);
+   for (int i = 0; i < n_mat; ++i)
+   {
+      DenseMatrix &A = A_batch(i);
+      A = 0.0;
+      for (int j = 0; j < n; ++j)
+      {
+         A(j, j) = 2.0 + i + j;
+      }
+   }
+
+   DenseMatrix &singular = A_batch(1);
+   singular = 0.0;
+   singular(0, 0) = 1.0;
+   singular(2, 2) = 1.0;
+
+   return A_batch;
+}
+
+}
+
+TEST_CASE("Batched LU factorization failure handling",
+          "[DenseMatrix][GPU]")
+{
+   auto backend = GENERATE(BatchedLinAlg::NATIVE,
+                           BatchedLinAlg::GPU_BLAS,
+                           BatchedLinAlg::MAGMA);
+   if (!BatchedLinAlg::IsAvailable(backend)) { return; }
+   CAPTURE(backend);
+
+   SECTION("LUFactor")
+   {
+      DenseTensor A_batch = MakeSingularBatchedMatrices();
+      Array<int> P;
+      REQUIRE_THROWS_WITH(BatchedLinAlg::Get(backend).LUFactor(A_batch, P),
+                          Catch::Matchers::Contains(
+                             "Batch LU factorization failed"));
+   }
+
+   SECTION("Invert")
+   {
+      DenseTensor A_batch = MakeSingularBatchedMatrices();
+      REQUIRE_THROWS_WITH(BatchedLinAlg::Get(backend).Invert(A_batch),
+                          Catch::Matchers::Contains(
+                             "Batch LU factorization failed"));
+   }
+}
+#endif
+
 TEST_CASE("DenseTensor copy", "[DenseMatrix][DenseTensor]")
 {
    DenseTensor t1(2,3,4);