Add bool ownership flag to AddBdrConstraintIntegrator

Prototype of "maybe owning" smart pointer semantics.

examples/ex4.cpp

@@ -65,6 +65,7 @@ int main(int argc, char *argv[])
    bool static_cond = false;
    bool hybridization = false;
    bool pa = false;
+   bool ea = false;
    const char *device_config = "cpu";
    bool visualization = 1;
 
@@ -83,18 +84,14 @@ int main(int argc, char *argv[])
                   "--no-hybridization", "Enable hybridization.");
    args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
                   "--no-partial-assembly", "Enable Partial Assembly.");
+   args.AddOption(&ea, "-ea", "--element-assembly", "-no-ea",
+                  "--no-element-assembly", "Enable Element Assembly.");
    args.AddOption(&device_config, "-d", "--device",
                   "Device configuration string, see Device::Configure().");
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
-   args.Parse();
-   if (!args.Good())
-   {
-      args.PrintUsage(cout);
-      return 1;
-   }
-   args.PrintOptions(cout);
+   args.ParseCheck();
    kappa = freq * M_PI;
 
    // 2. Enable hardware devices such as GPUs, and programming models such as
@@ -166,6 +163,7 @@ int main(int argc, char *argv[])
    Coefficient *beta  = new ConstantCoefficient(1.0);
    BilinearForm *a = new BilinearForm(fespace);
    if (pa) { a->SetAssemblyLevel(AssemblyLevel::PARTIAL); }
+   if (ea) { a->SetAssemblyLevel(AssemblyLevel::ELEMENT); }
    a->AddDomainIntegrator(new DivDivIntegrator(*alpha));
    a->AddDomainIntegrator(new VectorFEMassIntegrator(*beta));
 

examples/ex4p.cpp

@@ -71,6 +71,7 @@ int main(int argc, char *argv[])
    bool static_cond = false;
    bool hybridization = false;
    bool pa = false;
+   bool ea = false;
    const char *device_config = "cpu";
    bool visualization = 1;
 
@@ -89,24 +90,14 @@ int main(int argc, char *argv[])
                   "--no-hybridization", "Enable hybridization.");
    args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
                   "--no-partial-assembly", "Enable Partial Assembly.");
+   args.AddOption(&ea, "-ea", "--element-assembly", "-no-ea",
+                  "--no-element-assembly", "Enable Element Assembly.");
    args.AddOption(&device_config, "-d", "--device",
                   "Device configuration string, see Device::Configure().");
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
-   args.Parse();
-   if (!args.Good())
-   {
-      if (myid == 0)
-      {
-         args.PrintUsage(cout);
-      }
-      return 1;
-   }
-   if (myid == 0)
-   {
-      args.PrintOptions(cout);
-   }
+   args.ParseCheck();
    kappa = freq * M_PI;
 
    // 3. Enable hardware devices such as GPUs, and programming models such as
@@ -194,6 +185,7 @@ int main(int argc, char *argv[])
    Coefficient *beta  = new ConstantCoefficient(1.0);
    ParBilinearForm *a = new ParBilinearForm(fespace);
    if (pa) { a->SetAssemblyLevel(AssemblyLevel::PARTIAL); }
+   if (ea) { a->SetAssemblyLevel(AssemblyLevel::ELEMENT); }
    a->AddDomainIntegrator(new DivDivIntegrator(*alpha));
    a->AddDomainIntegrator(new VectorFEMassIntegrator(*beta));
 

fem/CMakeLists.txt

@@ -35,6 +35,7 @@ set(SRCS
   integ/bilininteg_mass_ea.cpp
   integ/bilininteg_mixedcurl_pa.cpp
   integ/bilininteg_mixedvecgrad_pa.cpp
+  integ/bilininteg_trace_jump_ea.cpp
   integ/bilininteg_transpose_ea.cpp
   integ/bilininteg_vecdiffusion_mf.cpp
   integ/bilininteg_vecdiffusion_pa.cpp
@@ -46,6 +47,7 @@ set(SRCS
   integ/bilininteg_diffusion_kernels.cpp
   integ/bilininteg_elasticity_kernels.cpp
   integ/bilininteg_hcurl_kernels.cpp
+  integ/bilininteg_hdiv_ea.cpp
   integ/bilininteg_hdiv_kernels.cpp
   integ/bilininteg_hcurlhdiv_kernels.cpp
   integ/bilininteg_mass_kernels.cpp
@@ -80,6 +82,7 @@ set(SRCS
   geom.cpp
   gridfunc.cpp
   hybridization.cpp
+  hybridization_ext.cpp
   intrules.cpp
   intrules_cut.cpp
   ceed/interface/basis.cpp
@@ -188,6 +191,7 @@ set(HDRS
   geom.hpp
   gridfunc.hpp
   hybridization.hpp
+  hybridization_ext.hpp
   intrules.hpp
   intrules_cut.hpp
   kernel_dispatch.hpp

fem/bilinearform.cpp

@@ -71,15 +71,11 @@ BilinearForm::BilinearForm(FiniteElementSpace * f)
    sequence = f->GetSequence();
    mat = mat_e = NULL;
    extern_bfs = 0;
-   element_matrices = NULL;
-   static_cond = NULL;
-   hybridization = NULL;
    precompute_sparsity = 0;
    diag_policy = DIAG_KEEP;
 
    assembly = AssemblyLevel::LEGACY;
    batch = 1;
-   ext = NULL;
 }
 
 BilinearForm::BilinearForm (FiniteElementSpace * f, BilinearForm * bf, int ps)
@@ -89,15 +85,11 @@ BilinearForm::BilinearForm (FiniteElementSpace * f, BilinearForm * bf, int ps)
    sequence = f->GetSequence();
    mat_e = NULL;
    extern_bfs = 1;
-   element_matrices = NULL;
-   static_cond = NULL;
-   hybridization = NULL;
    precompute_sparsity = ps;
    diag_policy = DIAG_KEEP;
 
    assembly = AssemblyLevel::LEGACY;
    batch = 1;
-   ext = NULL;
 
    // Copy the pointers to the integrators
    domain_integs = bf->domain_integs;
@@ -127,16 +119,16 @@ void BilinearForm::SetAssemblyLevel(AssemblyLevel assembly_level)
          break;
       case AssemblyLevel::FULL:
          SetDiagonalPolicy( DIAG_ONE ); // Only diagonal policy supported on device
-         ext = new FABilinearFormExtension(this);
+         ext.reset(new FABilinearFormExtension(this));
          break;
       case AssemblyLevel::ELEMENT:
-         ext = new EABilinearFormExtension(this);
+         ext.reset(new EABilinearFormExtension(this));
          break;
       case AssemblyLevel::PARTIAL:
-         ext = new PABilinearFormExtension(this);
+         ext.reset(new PABilinearFormExtension(this));
          break;
       case AssemblyLevel::NONE:
-         ext = new MFBilinearFormExtension(this);
+         ext.reset(new MFBilinearFormExtension(this));
          break;
       default:
          MFEM_ABORT("BilinearForm: unknown assembly level");
@@ -145,14 +137,13 @@ void BilinearForm::SetAssemblyLevel(AssemblyLevel assembly_level)
 
 void BilinearForm::EnableStaticCondensation()
 {
-   delete static_cond;
    if (assembly != AssemblyLevel::LEGACY)
    {
-      static_cond = NULL;
+      static_cond.reset();
       MFEM_WARNING("Static condensation not supported for this assembly level");
       return;
    }
-   static_cond = new StaticCondensation(fes);
+   static_cond.reset(new StaticCondensation(fes));
    if (static_cond->ReducesTrueVSize())
    {
       bool symmetric = false;      // TODO
@@ -161,8 +152,7 @@ void BilinearForm::EnableStaticCondensation()
    }
    else
    {
-      delete static_cond;
-      static_cond = NULL;
+      static_cond.reset();
    }
 }
 
@@ -170,15 +160,18 @@ void BilinearForm::EnableHybridization(FiniteElementSpace *constr_space,
                                        BilinearFormIntegrator *constr_integ,
                                        const Array<int> &ess_tdof_list)
 {
-   delete hybridization;
-   if (assembly != AssemblyLevel::LEGACY)
+   if (assembly != AssemblyLevel::LEGACY && assembly != AssemblyLevel::ELEMENT)
    {
       delete constr_integ;
-      hybridization = NULL;
+      hybridization.reset();
       MFEM_WARNING("Hybridization not supported for this assembly level");
       return;
    }
-   hybridization = new Hybridization(fes, constr_space);
+   hybridization.reset(new Hybridization(fes, constr_space));
+   if (assembly == AssemblyLevel::ELEMENT)
+   {
+      hybridization->EnableDeviceExecution();
+   }
    hybridization->SetConstraintIntegrator(constr_integ);
    hybridization->Init(ess_tdof_list);
 }
@@ -231,8 +224,8 @@ void BilinearForm::Finalize (int skip_zeros)
       if (!static_cond) { mat->Finalize(skip_zeros); }
       if (mat_e) { mat_e->Finalize(skip_zeros); }
       if (static_cond) { static_cond->Finalize(); }
-      if (hybridization) { hybridization->Finalize(); }
    }
+   if (hybridization) { hybridization->Finalize(); }
 }
 
 void BilinearForm::AddDomainIntegrator(BilinearFormIntegrator *bfi)
@@ -465,6 +458,10 @@ void BilinearForm::Assemble(int skip_zeros)
    if (ext)
    {
       ext->Assemble();
+      if (hybridization)
+      {
+         hybridization->AssembleElementMatrices(GetElementMatrices());
+      }
       return;
    }
 
@@ -834,7 +831,19 @@ void BilinearForm::FormLinearSystem(const Array<int> &ess_tdof_list, Vector &x,
 {
    if (ext)
    {
-      ext->FormLinearSystem(ess_tdof_list, x, b, A, X, B, copy_interior);
+      if (hybridization)
+      {
+         FormSystemMatrix(ess_tdof_list, A);
+         ConstrainedOperator A_constrained(this, ess_tdof_list);
+         A_constrained.EliminateRHS(x, b);
+         hybridization->ReduceRHS(b, B);
+         X.SetSize(B.Size());
+         X = 0.0;
+      }
+      else
+      {
+         ext->FormLinearSystem(ess_tdof_list, x, b, A, X, B, copy_interior);
+      }
       return;
    }
    const SparseMatrix *P = fes->GetConformingProlongation();
@@ -902,7 +911,16 @@ void BilinearForm::FormSystemMatrix(const Array<int> &ess_tdof_list,
 {
    if (ext)
    {
-      ext->FormSystemMatrix(ess_tdof_list, A);
+      if (hybridization)
+      {
+         const int remove_zeros = 0;
+         Finalize(remove_zeros);
+         A.Reset(&hybridization->GetMatrix(), false);
+      }
+      else
+      {
+         ext->FormSystemMatrix(ess_tdof_list, A);
+      }
       return;
    }
 
@@ -943,7 +961,7 @@ void BilinearForm::FormSystemMatrix(const Array<int> &ess_tdof_list,
 void BilinearForm::RecoverFEMSolution(const Vector &X,
                                       const Vector &b, Vector &x)
 {
-   if (ext)
+   if (ext && !hybridization)
    {
       ext->RecoverFEMSolution(X, b, x);
       return;
@@ -1000,16 +1018,26 @@ void BilinearForm::RecoverFEMSolution(const Vector &X,
 
 void BilinearForm::ComputeElementMatrices()
 {
-   if (element_matrices || domain_integs.Size() == 0 || fes->GetNE() == 0)
+   if (element_matrices) { return; }
+
+   if (auto *ea_ext = dynamic_cast<EABilinearFormExtension*>(ext.get()))
    {
+      element_matrices.reset(new DenseTensor);
+      ea_ext->GetElementMatrices(*element_matrices, ElementDofOrdering::NATIVE, true);
+      return;
+   }
+
+   if (domain_integs.Size() == 0 || fes->GetNE() == 0)
+   {
+      element_matrices.reset(new DenseTensor);
       return;
    }
 
    int num_elements = fes->GetNE();
    int num_dofs_per_el = fes->GetFE(0)->GetDof() * fes->GetVDim();
 
-   element_matrices = new DenseTensor(num_dofs_per_el, num_dofs_per_el,
-                                      num_elements);
+   element_matrices.reset(new DenseTensor(num_dofs_per_el, num_dofs_per_el,
+                                          num_elements));
 
    DenseMatrix tmp;
    IsoparametricTransformation eltrans;
@@ -1040,6 +1068,12 @@ void BilinearForm::ComputeElementMatrices()
    }
 }
 
+const DenseTensor &BilinearForm::GetElementMatrices()
+{
+   ComputeElementMatrices(); // Won't recompute if element_matrices exists
+   return *element_matrices;
+}
+
 void BilinearForm::EliminateEssentialBC(const Array<int> &bdr_attr_is_ess,
                                         const Vector &sol, Vector &rhs,
                                         DiagonalPolicy dpolicy)
@@ -1227,15 +1261,13 @@ void BilinearForm::Update(FiniteElementSpace *nfes)
    delete mat_e;
    mat_e = NULL;
    FreeElementMatrices();
-   delete static_cond;
-   static_cond = NULL;
+   static_cond.reset();
 
    if (full_update)
    {
       delete mat;
       mat = NULL;
-      delete hybridization;
-      hybridization = NULL;
+      hybridization.reset();
       sequence = fes->GetSequence();
    }
    else
@@ -1258,9 +1290,6 @@ BilinearForm::~BilinearForm()
 {
    delete mat_e;
    delete mat;
-   delete element_matrices;
-   delete static_cond;
-   delete hybridization;
 
    if (!extern_bfs)
    {
@@ -1272,8 +1301,6 @@ BilinearForm::~BilinearForm()
       for (k=0; k < boundary_face_integs.Size(); k++)
       { delete boundary_face_integs[k]; }
    }
-
-   delete ext;
 }
 
 
@@ -1300,7 +1327,6 @@ MixedBilinearForm::MixedBilinearForm (FiniteElementSpace *tr_fes,
    mat = NULL;
    mat_e = NULL;
    extern_bfs = 1;
-   ext = NULL;
 
    // Copy the pointers to the integrators
    domain_integs = mbf->domain_integs;
@@ -1330,22 +1356,22 @@ void MixedBilinearForm::SetAssemblyLevel(AssemblyLevel assembly_level)
       case AssemblyLevel::LEGACY:
          break;
       case AssemblyLevel::FULL:
-         // ext = new FAMixedBilinearFormExtension(this);
+         // ext.reset(new FAMixedBilinearFormExtension(this));
          // Use the original BilinearForm implementation for now
          break;
       case AssemblyLevel::ELEMENT:
-         mfem_error("Element assembly not supported yet... stay tuned!");
-         // ext = new EAMixedBilinearFormExtension(this);
+         MFEM_ABORT("Element assembly not supported yet... stay tuned!");
+         // ext.reset(new EAMixedBilinearFormExtension(this));
          break;
       case AssemblyLevel::PARTIAL:
-         ext = new PAMixedBilinearFormExtension(this);
+         ext.reset(new PAMixedBilinearFormExtension(this));
          break;
       case AssemblyLevel::NONE:
-         mfem_error("Matrix-free action not supported yet... stay tuned!");
-         // ext = new MFMixedBilinearFormExtension(this);
+         MFEM_ABORT("Matrix-free action not supported yet... stay tuned!");
+         // ext.reset(new MFMixedBilinearFormExtension(this));
          break;
       default:
-         mfem_error("Unknown assembly level");
+         MFEM_ABORT("Unknown assembly level");
    }
 }
 
@@ -2342,7 +2368,6 @@ MixedBilinearForm::~MixedBilinearForm()
       for (i = 0; i < boundary_trace_face_integs.Size(); i++)
       { delete boundary_trace_face_integs[i]; }
    }
-   delete ext;
 }
 
 void DiscreteLinearOperator::SetAssemblyLevel(AssemblyLevel assembly_level)
@@ -2359,16 +2384,16 @@ void DiscreteLinearOperator::SetAssemblyLevel(AssemblyLevel assembly_level)
          // Use the original implementation for now
          break;
       case AssemblyLevel::ELEMENT:
-         mfem_error("Element assembly not supported yet... stay tuned!");
+         MFEM_ABORT("Element assembly not supported yet... stay tuned!");
          break;
       case AssemblyLevel::PARTIAL:
-         ext = new PADiscreteLinearOperatorExtension(this);
+         ext.reset(new PADiscreteLinearOperatorExtension(this));
          break;
       case AssemblyLevel::NONE:
-         mfem_error("Matrix-free action not supported yet... stay tuned!");
+         MFEM_ABORT("Matrix-free action not supported yet... stay tuned!");
          break;
       default:
-         mfem_error("Unknown assembly level");
+         MFEM_ABORT("Unknown assembly level");
    }
 }
 

fem/bilinearform.hpp

@@ -83,7 +83,7 @@ protected:
    /** @brief Extension for supporting Full Assembly (FA),
        Element Assembly (EA),Partial Assembly (PA),
        or Matrix Free assembly (MF). */
-   BilinearFormExtension *ext;
+   std::unique_ptr<BilinearFormExtension> ext;
 
    /** Indicates if the sparse matrix is sorted after assembly when using
        Full Assembly (FA). */
@@ -122,10 +122,10 @@ protected:
    mutable DenseMatrix elemmat;
    mutable Array<int>  vdofs;
 
-   DenseTensor *element_matrices; ///< Owned.
+   std::unique_ptr<DenseTensor> element_matrices;
 
-   StaticCondensation *static_cond; ///< Owned.
-   Hybridization *hybridization; ///< Owned.
+   std::unique_ptr<StaticCondensation> static_cond;
+   std::unique_ptr<Hybridization> hybridization;
 
    /** @brief This data member allows one to specify what should be done to the
        diagonal matrix entries and corresponding RHS values upon elimination of
@@ -148,13 +148,11 @@ protected:
    BilinearForm() : Matrix (0)
    {
       fes = NULL; sequence = -1;
-      mat = mat_e = NULL; extern_bfs = 0; element_matrices = NULL;
-      static_cond = NULL; hybridization = NULL;
+      mat = mat_e = NULL; extern_bfs = 0;
       precompute_sparsity = 0;
       diag_policy = DIAG_KEEP;
       assembly = AssemblyLevel::LEGACY;
       batch = 1;
-      ext = NULL;
    }
 
 private:
@@ -214,7 +212,7 @@ public:
    /// Returns the assembly level
    AssemblyLevel GetAssemblyLevel() const { return assembly; }
 
-   Hybridization *GetHybridization() const { return hybridization; }
+   Hybridization *GetHybridization() const { return hybridization.get(); }
 
    /** @brief Enable the use of static condensation. For details see the
        description for class StaticCondensation in fem/staticcond.hpp This
@@ -224,7 +222,7 @@ public:
 
    /** @brief Check if static condensation was actually enabled by a previous
        call to EnableStaticCondensation(). */
-   bool StaticCondensationIsEnabled() const { return static_cond; }
+   bool StaticCondensationIsEnabled() const { return static_cond.get() != nullptr; }
 
    /// Return the trace FE space associated with static condensation.
    FiniteElementSpace *SCFESpace() const
@@ -569,12 +567,20 @@ public:
    void RecoverFEMSolution(const Vector &X, const Vector &b,
                            Vector &x) override;
 
-   /// Compute and store internally all element matrices.
+   /// @brief Compute and store internally all element matrices.
+   ///
+   /// If AssemblyLevel::ELEMENT is selected with SetAssemblyLeve(), this will
+   /// use effecient (device-accelerated) assembly of the element matrices.
    void ComputeElementMatrices();
 
    /// Free the memory used by the element matrices.
-   void FreeElementMatrices()
-   { delete element_matrices; element_matrices = NULL; }
+   void FreeElementMatrices() { element_matrices.reset(); }
+
+   /// @brief Return a DenseTensor containing the assembled element matrices.
+   ///
+   /// If AssemblyLevel::ELEMENT is selected with SetAssemblyLeve(), this will
+   /// use effecient (device-accelerated) assembly of the element matrices.
+   const DenseTensor &GetElementMatrices();
 
    /// Compute the element matrix of the given element
    /** The element matrix is computed by calling the domain integrators
@@ -753,7 +759,7 @@ protected:
 
    /** Extension for supporting Full Assembly (FA), Element Assembly (EA),
        Partial Assembly (PA), or Matrix Free assembly (MF). */
-   MixedBilinearFormExtension *ext;
+   std::unique_ptr<MixedBilinearFormExtension> ext;
 
    /** @brief Indicates the BilinearFormIntegrator%s stored in
        MixedBilinearForm#domain_integs, MixedBilinearForm#boundary_integs,

fem/bilinearform_ext.cpp

@@ -16,6 +16,7 @@
 #include "bilinearform.hpp"
 #include "pbilinearform.hpp"
 #include "pgridfunc.hpp"
+#include "fe/face_map_utils.hpp"
 #include "ceed/interface/util.hpp"
 
 namespace mfem
@@ -864,54 +865,139 @@ void EABilinearFormExtension::Assemble()
    ne = trial_fes->GetMesh()->GetNE();
    elemDofs = trial_fes->GetFE(0)->GetDof();
 
-   ea_data.SetSize(ne*elemDofs*elemDofs, Device::GetMemoryType());
-   ea_data.UseDevice(true);
+   Vector ea_data_tmp;
 
-   Array<BilinearFormIntegrator*> &integrators = *a->GetDBFI();
-   const int integratorCount = integrators.Size();
-   if ( integratorCount == 0 )
+   auto add_with_markers = [&](const Vector &ea_1, Vector &ea_2, const int ne_,
+                               const Array<int> &markers, const Array<int> &attrs,
+                               const bool add)
    {
-      ea_data = 0.0;
-   }
-   for (int i = 0; i < integratorCount; ++i)
+      if (ne_ == 0) { return; }
+      const int sz = ea_1.Size() / ne_;
+      const int *d_m = markers.Read();
+      const int *d_a = attrs.Read();
+      const auto d_ea_1 = Reshape(ea_1.Read(), sz, ne_);
+      auto d_ea_2 = Reshape(add ? ea_2.ReadWrite() : ea_2.Write(), sz, ne_);
+
+      mfem::forall(sz*ne_, [=] MFEM_HOST_DEVICE (int idx)
+      {
+         const int i = idx % sz;
+         const int e = idx / sz;
+         const real_t val = d_m[d_a[e] - 1] ? d_ea_1(i, e) : 0.0;
+         if (add)
+         {
+            d_ea_2(i, e) += val;
+         }
+         else
+         {
+            d_ea_2(i, e) = val;
+         }
+      });
+   };
+
    {
-      integrators[i]->AssembleEA(*a->FESpace(), ea_data, i);
+      ea_data.SetSize(ne*elemDofs*elemDofs);
+      ea_data.UseDevice(true);
+      Array<BilinearFormIntegrator*> &integrators = *a->GetDBFI();
+      Array<Array<int>*> &markers_array = *a->GetDBFI_Marker();
+
+      if (integrators.Size() == 0) { ea_data = 0.0; }
+
+      for (int i = 0; i < integrators.Size(); ++i)
+      {
+         const bool add = (i > 0);
+         const Array<int> *markers = markers_array[i];
+         if (markers == nullptr)
+         {
+            integrators[i]->AssembleEA(*a->FESpace(), ea_data, add);
+         }
+         else
+         {
+            ea_data_tmp.SetSize(ea_data.Size());
+            integrators[i]->AssembleEA(*a->FESpace(), ea_data_tmp, false);
+            add_with_markers(ea_data_tmp, ea_data, ne, *markers,
+                             elem_attributes, add);
+         }
+      }
    }
 
    faceDofs = trial_fes ->
               GetTraceElement(0, trial_fes->GetMesh()->GetFaceGeometry(0)) ->
               GetDof();
 
-   MFEM_VERIFY(a->GetBBFI()->Size() == 0,
-               "Element assembly does not support AddBoundaryIntegrator yet.");
-
-   Array<BilinearFormIntegrator*> &intFaceIntegrators = *a->GetFBFI();
-   const int intFaceIntegratorCount = intFaceIntegrators.Size();
-   if (intFaceIntegratorCount>0)
    {
-      nf_int = trial_fes->GetNFbyType(FaceType::Interior);
-      ea_data_int.SetSize(2*nf_int*faceDofs*faceDofs, Device::GetMemoryType());
-      ea_data_ext.SetSize(2*nf_int*faceDofs*faceDofs, Device::GetMemoryType());
-   }
-   for (int i = 0; i < intFaceIntegratorCount; ++i)
-   {
-      intFaceIntegrators[i]->AssembleEAInteriorFaces(*a->FESpace(),
-                                                     ea_data_int,
-                                                     ea_data_ext,
-                                                     i);
+      Array<BilinearFormIntegrator*> &bdr_integs = *a->GetBBFI();
+      Array<Array<int>*> &markers_array = *a->GetBBFI_Marker();
+      const int n_bdr_integs = bdr_integs.Size();
+      if (n_bdr_integs > 0)
+      {
+         nf_bdr = trial_fes->GetNFbyType(FaceType::Boundary);
+         ea_data_bdr.SetSize(nf_bdr*faceDofs*faceDofs);
+      }
+      for (int i = 0; i < n_bdr_integs; ++i)
+      {
+         const bool add = (i > 0);
+         const Array<int> *markers = markers_array[i];
+         if (markers == nullptr)
+         {
+            bdr_integs[i]->AssembleEABoundary(*a->FESpace(), ea_data_bdr, add);
+         }
+         else
+         {
+            ea_data_tmp.SetSize(ea_data_bdr.Size());
+            bdr_integs[i]->AssembleEABoundary(*a->FESpace(), ea_data_tmp, add);
+            add_with_markers(ea_data_tmp, ea_data_bdr, nf_bdr, *markers,
+                             bdr_attributes, add);
+         }
+      }
    }
 
-   Array<BilinearFormIntegrator*> &bdrFaceIntegrators = *a->GetBFBFI();
-   const int boundFaceIntegratorCount = bdrFaceIntegrators.Size();
-   if (boundFaceIntegratorCount>0)
    {
-      nf_bdr = trial_fes->GetNFbyType(FaceType::Boundary);
-      ea_data_bdr.SetSize(nf_bdr*faceDofs*faceDofs, Device::GetMemoryType());
-      ea_data_bdr = 0.0;
+      Array<BilinearFormIntegrator*> &intFaceIntegrators = *a->GetFBFI();
+      const int intFaceIntegratorCount = intFaceIntegrators.Size();
+      if (intFaceIntegratorCount>0)
+      {
+         nf_int = trial_fes->GetNFbyType(FaceType::Interior);
+         ea_data_int.SetSize(2*nf_int*faceDofs*faceDofs);
+         ea_data_ext.SetSize(2*nf_int*faceDofs*faceDofs);
+      }
+      for (int i = 0; i < intFaceIntegratorCount; ++i)
+      {
+         const bool add = (i > 0);
+         intFaceIntegrators[i]->AssembleEAInteriorFaces(*a->FESpace(),
+                                                        ea_data_int,
+                                                        ea_data_ext,
+                                                        add);
+      }
    }
-   for (int i = 0; i < boundFaceIntegratorCount; ++i)
+
    {
-      bdrFaceIntegrators[i]->AssembleEABoundaryFaces(*a->FESpace(),ea_data_bdr,i);
+      Array<BilinearFormIntegrator*> &bdr_face_integs = *a->GetBFBFI();
+      Array<Array<int>*> &markers_array = *a->GetBFBFI_Marker();
+      const int n_bdr_face_integs = bdr_face_integs.Size();
+      if (n_bdr_face_integs > 0)
+      {
+         nf_bdr = trial_fes->GetNFbyType(FaceType::Boundary);
+         ea_data_bdr.SetSize(nf_bdr*faceDofs*faceDofs);
+      }
+      for (int i = 0; i < n_bdr_face_integs; ++i)
+      {
+         const bool add = (i > 0);
+         const Array<int> *markers = markers_array[i];
+         if (markers == nullptr)
+         {
+            bdr_face_integs[i]->AssembleEABoundaryFaces(
+               *a->FESpace(), ea_data_bdr, add);
+         }
+         else
+         {
+            ea_data_tmp.SetSize(ea_data_bdr.Size());
+            bdr_face_integs[i]->AssembleEABoundaryFaces(*a->FESpace(),
+                                                        ea_data_tmp,
+                                                        add);
+            add_with_markers(ea_data_tmp, ea_data_bdr, nf_bdr, *markers,
+                             bdr_attributes, add);
+         }
+      }
    }
 
    if (factorize_face_terms && int_face_restrict_lex)
@@ -1023,34 +1109,29 @@ void EABilinearFormExtension::Mult(const Vector &x, Vector &y) const
    }
 
    // Treatment of boundary faces
-   Array<BilinearFormIntegrator*> &bdrFaceIntegrators = *a->GetBFBFI();
-   const int bFISz = bdrFaceIntegrators.Size();
-   if (!factorize_face_terms && bdr_face_restrict_lex && bFISz>0)
+   if (!factorize_face_terms && bdr_face_restrict_lex && ea_data_bdr.Size() > 0)
    {
       // Apply the Boundary Face Restriction
       bdr_face_restrict_lex->Mult(x, bdr_face_X);
-      if (bdr_face_X.Size()>0)
+      bdr_face_Y = 0.0;
+      // Apply the boundary face matrices
+      const int NDOFS = faceDofs;
+      auto X = Reshape(bdr_face_X.Read(), NDOFS, nf_bdr);
+      auto Y = Reshape(bdr_face_Y.ReadWrite(), NDOFS, nf_bdr);
+      auto A = Reshape(ea_data_bdr.Read(), NDOFS, NDOFS, nf_bdr);
+      mfem::forall(nf_bdr*NDOFS, [=] MFEM_HOST_DEVICE (int glob_j)
       {
-         bdr_face_Y = 0.0;
-         // Apply the boundary face matrices
-         const int NDOFS = faceDofs;
-         auto X = Reshape(bdr_face_X.Read(), NDOFS, nf_bdr);
-         auto Y = Reshape(bdr_face_Y.ReadWrite(), NDOFS, nf_bdr);
-         auto A = Reshape(ea_data_bdr.Read(), NDOFS, NDOFS, nf_bdr);
-         mfem::forall(nf_bdr*NDOFS, [=] MFEM_HOST_DEVICE (int glob_j)
+         const int f = glob_j/NDOFS;
+         const int j = glob_j%NDOFS;
+         real_t res = 0.0;
+         for (int i = 0; i < NDOFS; i++)
          {
-            const int f = glob_j/NDOFS;
-            const int j = glob_j%NDOFS;
-            real_t res = 0.0;
-            for (int i = 0; i < NDOFS; i++)
-            {
-               res += A(i, j, f)*X(i, f);
-            }
-            Y(j, f) += res;
-         });
-         // Apply the Boundary Face Restriction transposed
-         bdr_face_restrict_lex->AddMultTransposeInPlace(bdr_face_Y, y);
-      }
+            res += A(i, j, f)*X(i, f);
+         }
+         Y(j, f) += res;
+      });
+      // Apply the Boundary Face Restriction transposed
+      bdr_face_restrict_lex->AddMultTransposeInPlace(bdr_face_Y, y);
    }
 }
 
@@ -1151,35 +1232,164 @@ void EABilinearFormExtension::MultTranspose(const Vector &x, Vector &y) const
    }
 
    // Treatment of boundary faces
-   Array<BilinearFormIntegrator*> &bdrFaceIntegrators = *a->GetBFBFI();
-   const int bFISz = bdrFaceIntegrators.Size();
-   if (!factorize_face_terms && bdr_face_restrict_lex && bFISz>0)
+   if (!factorize_face_terms && bdr_face_restrict_lex && ea_data_bdr.Size() > 0)
    {
       // Apply the Boundary Face Restriction
       bdr_face_restrict_lex->Mult(x, bdr_face_X);
-      if (bdr_face_X.Size()>0)
+      bdr_face_Y = 0.0;
+      // Apply the boundary face matrices transposed
+      const int NDOFS = faceDofs;
+      auto X = Reshape(bdr_face_X.Read(), NDOFS, nf_bdr);
+      auto Y = Reshape(bdr_face_Y.ReadWrite(), NDOFS, nf_bdr);
+      auto A = Reshape(ea_data_bdr.Read(), NDOFS, NDOFS, nf_bdr);
+      mfem::forall(nf_bdr*NDOFS, [=] MFEM_HOST_DEVICE (int glob_j)
       {
-         bdr_face_Y = 0.0;
-         // Apply the boundary face matrices transposed
-         const int NDOFS = faceDofs;
-         auto X = Reshape(bdr_face_X.Read(), NDOFS, nf_bdr);
-         auto Y = Reshape(bdr_face_Y.ReadWrite(), NDOFS, nf_bdr);
-         auto A = Reshape(ea_data_bdr.Read(), NDOFS, NDOFS, nf_bdr);
-         mfem::forall(nf_bdr*NDOFS, [=] MFEM_HOST_DEVICE (int glob_j)
+         const int f = glob_j/NDOFS;
+         const int j = glob_j%NDOFS;
+         real_t res = 0.0;
+         for (int i = 0; i < NDOFS; i++)
          {
-            const int f = glob_j/NDOFS;
-            const int j = glob_j%NDOFS;
-            real_t res = 0.0;
-            for (int i = 0; i < NDOFS; i++)
-            {
-               res += A(j, i, f)*X(i, f);
-            }
-            Y(j, f) += res;
-         });
-         // Apply the Boundary Face Restriction transposed
-         bdr_face_restrict_lex->AddMultTransposeInPlace(bdr_face_Y, y);
+            res += A(j, i, f)*X(i, f);
+         }
+         Y(j, f) += res;
+      });
+      // Apply the Boundary Face Restriction transposed
+      bdr_face_restrict_lex->AddMultTransposeInPlace(bdr_face_Y, y);
+   }
+}
+
+void EABilinearFormExtension::GetElementMatrices(
+   DenseTensor &element_matrices, ElementDofOrdering ordering, bool add_bdr)
+{
+   // Ensure the EA data is assembled
+   if (ea_data.Size() == 0) { Assemble(); }
+
+   const int ndofs = elemDofs;
+   element_matrices.SetSize(ndofs, ndofs, ne);
+   const int N = element_matrices.TotalSize();
+
+   const auto d_ea_data = Reshape(ea_data.Read(), ndofs, ndofs, ne);
+   auto d_element_matrices = Reshape(element_matrices.Write(),
+                                     ndofs, ndofs,
+                                     ne);
+
+   const int *d_dof_map = nullptr;
+   Array<int> dof_map;
+   if (ordering == ElementDofOrdering::NATIVE)
+   {
+      const TensorBasisElement* tbe =
+         dynamic_cast<const TensorBasisElement*>(trial_fes->GetFE(0));
+      if (tbe)
+      {
+         // Deep copy to avoid issues with host device (see similar comment in
+         // HybridizationExtension::ConstructC).
+         dof_map = tbe->GetDofMap();
+         d_dof_map = dof_map.Read();
       }
    }
+
+   if (d_dof_map)
+   {
+      // Reordering required
+      mfem::forall(N, [=] MFEM_HOST_DEVICE (int idx)
+      {
+         const int e = idx / ndofs / ndofs;
+         const int i = idx % ndofs;
+         const int j = (idx / ndofs) % ndofs;
+         const int ii_s = d_dof_map[i];
+         const int ii = (ii_s >= 0) ? ii_s : -1 - ii_s;
+         const int s_i = (ii_s >= 0) ? 1 : -1;
+         const int jj_s = d_dof_map[j];
+         const int jj = (jj_s >= 0) ? jj_s : -1 - jj_s;
+         const int s_j = (jj_s >= 0) ? 1 : -1;
+         d_element_matrices(ii, jj, e) = s_i*s_j*d_ea_data(j, i, e);
+      });
+   }
+   else
+   {
+      // No reordering required
+      mfem::forall(N, [=] MFEM_HOST_DEVICE (int idx)
+      {
+         const int e = idx / ndofs / ndofs;
+         const int i = idx % ndofs;
+         const int j = (idx / ndofs) % ndofs;
+         d_element_matrices(i, j, e) = d_ea_data(j, i, e);
+      });
+   }
+
+   if (add_bdr && ea_data_bdr.Size() > 0)
+   {
+      const int ndof_face = faceDofs;
+      const auto d_ea_bdr = Reshape(ea_data_bdr.Read(),
+                                    ndof_face, ndof_face, nf_bdr);
+
+      // Get all the local face maps (mapping from lexicographic face index to
+      // lexicographic volume index, depending on the local face index).
+      const Mesh &mesh = *trial_fes->GetMesh();
+      const int dim = mesh.Dimension();
+      const int n_faces_per_el = 2*dim; // assuming tensor product
+      Array<int> face_maps(ndof_face * n_faces_per_el);
+      for (int lf_i = 0; lf_i < n_faces_per_el; ++lf_i)
+      {
+         Array<int> face_map(ndof_face);
+         trial_fes->GetFE(0)->GetFaceMap(lf_i, face_map);
+         for (int i = 0; i < ndof_face; ++i)
+         {
+            face_maps[i + lf_i*ndof_face] = face_map[i];
+         }
+      }
+
+      Array<int> face_info(nf_bdr * 2);
+      {
+         int fidx = 0;
+         for (int f = 0; f < mesh.GetNumFaces(); ++f)
+         {
+            Mesh::FaceInformation finfo = mesh.GetFaceInformation(f);
+            if (!finfo.IsBoundary()) { continue; }
+            face_info[0 + fidx*2] = finfo.element[0].local_face_id;
+            face_info[1 + fidx*2] = finfo.element[0].index;
+            fidx++;
+         }
+      }
+
+      const auto d_face_maps = Reshape(face_maps.Read(), ndof_face, n_faces_per_el);
+      const auto d_face_info = Reshape(face_info.Read(), 2, nf_bdr);
+
+      const bool reorder = (ordering == ElementDofOrdering::NATIVE);
+
+      mfem::forall_2D(nf_bdr, ndof_face, ndof_face, [=] MFEM_HOST_DEVICE (int f)
+      {
+         const int lf_i = d_face_info(0, f);
+         const int e = d_face_info(1, f);
+         // Loop over face indices in "native ordering"
+         MFEM_FOREACH_THREAD(i_lex_face, x, ndof_face)
+         {
+            // Convert from lexicographic face DOF to volume DOF
+            const int i_lex = d_face_maps(i_lex_face, lf_i);
+
+            const int ii_s = d_dof_map[i_lex];
+            const int ii = (ii_s >= 0) ? ii_s : -1 - ii_s;
+
+            const int i = reorder ? ii : i_lex;
+            const int s_i = (ii_s < 0 && reorder) ? -1 : 1;
+
+            MFEM_FOREACH_THREAD(j_lex_face, y, ndof_face)
+            {
+               // Convert from lexicographic face DOF to volume DOF
+               const int j_lex = d_face_maps(j_lex_face, lf_i);
+
+               const int jj_s = d_dof_map[j_lex];
+               const int jj = (jj_s >= 0) ? jj_s : -1 - jj_s;
+
+               const int j = reorder ? jj : j_lex;
+               const int s_j = (jj_s < 0 && reorder) ? -1 : 1;
+
+               AtomicAdd(d_element_matrices(i, j, e),
+                         s_i*s_j*d_ea_bdr(i_lex_face, j_lex_face, f));
+            }
+         }
+      });
+   }
 }
 
 // Data and methods for fully-assembled bilinear forms

fem/bilinearform_ext.hpp

@@ -154,6 +154,17 @@ public:
    void Assemble() override;
    void Mult(const Vector &x, Vector &y) const override;
    void MultTranspose(const Vector &x, Vector &y) const override;
+
+   /// @brief Populates @a element_matrices with the element matrices.
+   ///
+   /// The element matrices are converted from row-major (how they are stored in
+   /// @a ea_data) to column-major format.
+   ///
+   /// If @a ordering is ElementDofOrdering::NATIVE, then the matrices are
+   /// reordered from the lexicographic ordering used internally.
+   void GetElementMatrices(DenseTensor &element_matrices,
+                           ElementDofOrdering ordering,
+                           bool add_bdr);
 };
 
 /// Data and methods for fully-assembled bilinear forms

fem/bilininteg.cpp

@@ -72,6 +72,14 @@ void BilinearFormIntegrator::AssembleEA(const FiniteElementSpace &fes,
               "   is not implemented for this class.");
 }
 
+void BilinearFormIntegrator::AssembleEABoundary(const FiniteElementSpace &fes,
+                                                Vector &emat,
+                                                const bool add)
+{
+   MFEM_ABORT("BilinearFormIntegrator::AssembleEABoundary(...)\n"
+              "   is not implemented for this class.");
+}
+
 void BilinearFormIntegrator::AssembleEAInteriorFaces(const FiniteElementSpace
                                                      &fes,
                                                      Vector &ea_data_int,
@@ -82,6 +90,16 @@ void BilinearFormIntegrator::AssembleEAInteriorFaces(const FiniteElementSpace
               "   is not implemented for this class.");
 }
 
+void BilinearFormIntegrator::AssembleEAInteriorFaces(
+   const FiniteElementSpace &trial_fes,
+   const FiniteElementSpace &test_fes,
+   Vector &emat,
+   const bool add)
+{
+   MFEM_ABORT("BilinearFormIntegrator::AssembleEAInteriorFaces(...)\n"
+              "   is not implemented for this class.");
+}
+
 void BilinearFormIntegrator::AssembleEABoundaryFaces(const FiniteElementSpace
                                                      &fes,
                                                      Vector &ea_data_bdr,

fem/bilininteg.hpp

@@ -124,11 +124,25 @@ public:
    /// Assemble diagonal and add it to Vector @a diag.
    virtual void AssembleDiagonalMF(Vector &diag);
 
+   virtual void AssembleEABoundary(const FiniteElementSpace &fes,
+                                   Vector &ea_data_bdr,
+                                   const bool add = true);
+
    virtual void AssembleEAInteriorFaces(const FiniteElementSpace &fes,
                                         Vector &ea_data_int,
                                         Vector &ea_data_ext,
                                         const bool add = true);
 
+   /// @brief Method defining element assembly for mixed trace integrators.
+   ///
+   /// This is the element assembly analogue of AssembleFaceMatrix(const
+   /// FiniteElement&, const FiniteElement&, const FiniteElement&,
+   /// FaceElementTransformations&, DenseMatrix&).
+   virtual void AssembleEAInteriorFaces(const FiniteElementSpace &trial_fes,
+                                        const FiniteElementSpace &test_fes,
+                                        Vector &emat,
+                                        const bool add = true);
+
    virtual void AssembleEABoundaryFaces(const FiniteElementSpace &fes,
                                         Vector &ea_data_bdr,
                                         const bool add = true);
@@ -383,6 +397,7 @@ public:
    void AssembleEA(const FiniteElementSpace &fes, Vector &emat,
                    const bool add) override;
 
+   using BilinearFormIntegrator::AssembleEAInteriorFaces;
    void AssembleEAInteriorFaces(const FiniteElementSpace &fes,
                                 Vector &ea_data_int,
                                 Vector &ea_data_ext,
@@ -494,6 +509,7 @@ public:
    void AssembleEA(const FiniteElementSpace &fes, Vector &emat,
                    const bool add) override;
 
+   using BilinearFormIntegrator::AssembleEAInteriorFaces;
    void AssembleEAInteriorFaces(const FiniteElementSpace &fes,
                                 Vector &ea_data_int,
                                 Vector &ea_data_ext,
@@ -2334,6 +2350,8 @@ protected:
    const FaceGeometricFactors *face_geom; ///< Not owned
    int dim, ne, nq, dofs1D, quad1D;
 
+   void AssembleEA_(Vector &ea, const bool add);
+
 public:
 
    using ApplyKernelType = void(*)(const int, const Array<real_t>&,
@@ -2374,7 +2392,10 @@ public:
    void AssembleEA(const FiniteElementSpace &fes, Vector &emat,
                    const bool add) override;
 
-   void AssembleDiagonalPA(Vector &diag) override;
+   virtual void AssembleEABoundary(const FiniteElementSpace &fes, Vector &emat,
+                                   const bool add) override;
+
+   virtual void AssembleDiagonalPA(Vector &diag) override;
 
    void AssembleDiagonalMF(Vector &diag) override;
 
@@ -2874,12 +2895,14 @@ public:
                                ElementTransformation &Trans,
                                DenseMatrix &elmat) override;
 
-   void AssemblePA(const FiniteElementSpace &fes) override;
-   void AssemblePA(const FiniteElementSpace &trial_fes,
-                   const FiniteElementSpace &test_fes) override;
-   void AddMultPA(const Vector &x, Vector &y) const override;
-   void AddMultTransposePA(const Vector &x, Vector &y) const override;
-   void AssembleDiagonalPA(Vector& diag) override;
+   virtual void AssemblePA(const FiniteElementSpace &fes) override;
+   virtual void AssemblePA(const FiniteElementSpace &trial_fes,
+                           const FiniteElementSpace &test_fes) override;
+   virtual void AddMultPA(const Vector &x, Vector &y) const override;
+   virtual void AddMultTransposePA(const Vector &x, Vector &y) const override;
+   virtual void AssembleDiagonalPA(Vector& diag) override;
+   virtual void AssembleEA(const FiniteElementSpace &fes, Vector &emat,
+                           const bool add) override;
 
    const Coefficient *GetCoefficient() const { return Q; }
 };
@@ -2938,11 +2961,6 @@ class DivDivIntegrator: public BilinearFormIntegrator
 protected:
    Coefficient *Q;
 
-   using BilinearFormIntegrator::AssemblePA;
-   void AssemblePA(const FiniteElementSpace &fes) override;
-   void AddMultPA(const Vector &x, Vector &y) const override;
-   void AssembleDiagonalPA(Vector& diag) override;
-
 private:
 #ifndef MFEM_THREAD_SAFE
    Vector divshape, te_divshape;
@@ -2969,6 +2987,13 @@ public:
                                ElementTransformation &Trans,
                                DenseMatrix &elmat) override;
 
+   using BilinearFormIntegrator::AssemblePA;
+   void AssemblePA(const FiniteElementSpace &fes) override;
+   void AddMultPA(const Vector &x, Vector &y) const override;
+   void AssembleDiagonalPA(Vector& diag) override;
+   void AssembleEA(const FiniteElementSpace &fes, Vector &emat,
+                   const bool add) override;
+
    const Coefficient *GetCoefficient() const { return Q; }
 };
 
@@ -3270,6 +3295,7 @@ public:
 
    void AddMultPA(const Vector&, Vector&) const override;
 
+   using BilinearFormIntegrator::AssembleEAInteriorFaces;
    void AssembleEAInteriorFaces(const FiniteElementSpace& fes,
                                 Vector &ea_data_int,
                                 Vector &ea_data_ext,
@@ -3578,6 +3604,12 @@ public:
                            const FiniteElement &test_fe2,
                            FaceElementTransformations &Trans,
                            DenseMatrix &elmat) override;
+
+   using BilinearFormIntegrator::AssembleEAInteriorFaces;
+   void AssembleEAInteriorFaces(const FiniteElementSpace &trial_fes,
+                                const FiniteElementSpace &test_fes,
+                                Vector &emat,
+                                const bool add = true) override;
 };
 
 /** Integrator for the DPG form:$ \langle v, w \rangle $ over a face (the interface) where

fem/fe/face_map_utils.hpp

@@ -69,16 +69,16 @@ inline int ToLexOrdering2D(const int face_id, const int size1d, const int i)
 }
 
 /// @brief Given a face DOF index on a shared face, ordered lexicographically
-/// relative to element 1, return the corresponding face DOF index ordered
-/// lexicographically relative to element 2.
+/// relative to element the element (where the local face is face_id), and
+/// return the corresponding face DOF index ordered lexicographically relative
+/// to the face itself.
 MFEM_HOST_DEVICE
-inline int PermuteFace2D(const int face_id1, const int face_id2,
-                         const int orientation, const int size1d,
-                         const int index)
+inline int PermuteFace2D(const int face_id, const int orientation,
+                         const int size1d, const int index)
 {
    int new_index;
    // Convert from element 1 lex ordering to native ordering
-   if (face_id1 == 2 || face_id1 == 3)
+   if (face_id == 2 || face_id == 3)
    {
       new_index = size1d-1-index;
    }
@@ -91,7 +91,18 @@ inline int PermuteFace2D(const int face_id1, const int face_id2,
    {
       new_index = size1d-1-new_index;
    }
-   // Covert to element 2 lex ordering
+   return new_index;
+}
+
+/// @brief Given a face DOF index on a shared face, ordered lexicographically
+/// relative to element 1, return the corresponding face DOF index ordered
+/// lexicographically relative to element 2.
+MFEM_HOST_DEVICE
+inline int PermuteFace2D(const int face_id1, const int face_id2,
+                         const int orientation, const int size1d,
+                         const int index)
+{
+   const int new_index = PermuteFace2D(face_id1, orientation, size1d, index);
    return ToLexOrdering2D(face_id2, size1d, new_index);
 }
 
@@ -116,26 +127,22 @@ inline int ToLexOrdering3D(const int face_id, const int size1d, const int i,
    }
 }
 
-/// @brief Given the index of a face DOF in lexicographic ordering relative
-/// element 1, permute the index so that it is lexicographically ordered
-/// relative to element 2.
-///
-/// The given face corresponds to local face index @a face_id1 relative to
-/// element 1, and @a face_id2 (with @a orientation) relative to element 2.
+/// @brief Given the index of a face DOF in lexicographic ordering relative the
+/// element (where the local face id is @a face_id), permute the index so that
+/// it is lexicographically ordered relative to the face itself.
 MFEM_HOST_DEVICE
-inline int PermuteFace3D(const int face_id1, const int face_id2,
-                         const int orientation,
+inline int PermuteFace3D(const int face_id, const int orientation,
                          const int size1d, const int index)
 {
    int i=0, j=0, new_i=0, new_j=0;
    i = index%size1d;
    j = index/size1d;
    // Convert from lex ordering
-   if (face_id1==3 || face_id1==4)
+   if (face_id==3 || face_id==4)
    {
       i = size1d-1-i;
    }
-   else if (face_id1==0)
+   else if (face_id==0)
    {
       j = size1d-1-j;
    }
@@ -175,6 +182,23 @@ inline int PermuteFace3D(const int face_id1, const int face_id2,
          new_j = (size1d-1-j);
          break;
    }
+   return new_i + new_j*size1d;
+}
+
+/// @brief Given the index of a face DOF in lexicographic ordering relative
+/// element 1, permute the index so that it is lexicographically ordered
+/// relative to element 2.
+///
+/// The given face corresponds to local face index @a face_id1 relative to
+/// element 1, and @a face_id2 (with @a orientation) relative to element 2.
+MFEM_HOST_DEVICE
+inline int PermuteFace3D(const int face_id1, const int face_id2,
+                         const int orientation,
+                         const int size1d, const int index)
+{
+   const int new_index = PermuteFace3D(face_id1, orientation, size1d, index);
+   const int new_i = new_index%size1d;
+   const int new_j = new_index/size1d;
    return ToLexOrdering3D(face_id2, size1d, new_i, new_j);
 }
 

fem/fespace.cpp

@@ -1417,6 +1417,10 @@ const FaceRestriction *FiniteElementSpace::GetFaceRestriction(
             res = new NCL2FaceRestriction(*this, f_ordering, type, m);
          }
       }
+      else if (dynamic_cast<const DG_Interface_FECollection*>(fec))
+      {
+         res = new L2InterfaceFaceRestriction(*this, f_ordering, type);
+      }
       else
       {
          res = new ConformingFaceRestriction(*this, f_ordering, type);

fem/hybridization.cpp

@@ -10,6 +10,7 @@
 // CONTRIBUTING.md for details.
 
 #include "hybridization.hpp"
+#include "hybridization_ext.hpp"
 #include "gridfunc.hpp"
 
 #ifdef MFEM_USE_MPI
@@ -29,36 +30,21 @@ namespace mfem
 
 Hybridization::Hybridization(FiniteElementSpace *fespace,
                              FiniteElementSpace *c_fespace)
-   : fes(fespace), c_fes(c_fespace), c_bfi(NULL), extern_bdr_constr_integs(0),
-     Ct(NULL), H(NULL), Af_data(NULL), Af_ipiv(NULL)
+   : fes(*fespace), c_fes(*c_fespace)
 {
 #ifdef MFEM_USE_MPI
-   pC = P_pc = NULL;
    pH.SetType(Operator::Hypre_ParCSR);
 #endif
 }
 
-Hybridization::~Hybridization()
+void Hybridization::EnableDeviceExecution()
 {
-#ifdef MFEM_USE_MPI
-   delete P_pc;
-   delete pC;
-#endif
-   delete [] Af_ipiv;
-   delete [] Af_data;
-   delete H;
-   delete Ct;
-   delete c_bfi;
-   if (!extern_bdr_constr_integs)
-   {
-      for (int k=0; k < boundary_constraint_integs.Size(); k++)
-      { delete boundary_constraint_integs[k]; }
-   }
+   ext.reset(new HybridizationExtension(*this));
 }
 
 void Hybridization::ConstructC()
 {
-   const int NE = fes->GetNE();
+   const int NE = fes.GetNE();
    int num_hat_dofs = hat_offsets[NE];
    Array<int> vdofs, c_vdofs;
 
@@ -73,7 +59,7 @@ void Hybridization::ConstructC()
 
    int c_num_face_nbr_dofs = 0;
 #ifdef MFEM_USE_MPI
-   ParFiniteElementSpace *c_pfes = dynamic_cast<ParFiniteElementSpace*>(c_fes);
+   ParFiniteElementSpace *c_pfes = dynamic_cast<ParFiniteElementSpace*>(&c_fes);
    ParMesh *pmesh = c_pfes ? c_pfes->GetParMesh() : NULL;
    HYPRE_BigInt num_shared_slave_faces = 0, glob_num_shared_slave_faces = 0;
    if (c_pfes)
@@ -105,15 +91,15 @@ void Hybridization::ConstructC()
    }
 #endif
 
-   const int c_vsize = c_fes->GetVSize();
-   Ct = new SparseMatrix(num_hat_dofs, c_vsize + c_num_face_nbr_dofs);
+   const int c_vsize = c_fes.GetVSize();
+   Ct.reset(new SparseMatrix(num_hat_dofs, c_vsize + c_num_face_nbr_dofs));
 
    if (c_bfi)
    {
       const int skip_zeros = 1;
       DenseMatrix elmat;
       FaceElementTransformations *FTr;
-      Mesh *mesh = fes->GetMesh();
+      Mesh *mesh = fes.GetMesh();
       int num_faces = mesh->GetNumFaces();
       for (int i = 0; i < num_faces; i++)
       {
@@ -133,17 +119,17 @@ void Hybridization::ConstructC()
          {
             vdofs[s1+j] = o2 + j;
          }
-         c_fes->GetFaceVDofs(i, c_vdofs);
-         c_bfi->AssembleFaceMatrix(*c_fes->GetFaceElement(i),
-                                   *fes->GetFE(FTr->Elem1No),
-                                   *fes->GetFE(FTr->Elem2No),
+         c_fes.GetFaceVDofs(i, c_vdofs);
+         c_bfi->AssembleFaceMatrix(*c_fes.GetFaceElement(i),
+                                   *fes.GetFE(FTr->Elem1No),
+                                   *fes.GetFE(FTr->Elem2No),
                                    *FTr, elmat);
          // zero-out small elements in elmat
          elmat.Threshold(mtol * elmat.MaxMaxNorm());
          Ct->AddSubMatrix(vdofs, c_vdofs, elmat, skip_zeros);
       }
 
-      if (boundary_constraint_integs.Size())
+      if (!boundary_constraint_integs.empty())
       {
          const FiniteElement *fe1, *fe2;
          const FiniteElement *face_el;
@@ -152,14 +138,14 @@ void Hybridization::ConstructC()
          Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
                                     mesh->bdr_attributes.Max() : 0);
          bdr_attr_marker = 0;
-         for (int k = 0; k < boundary_constraint_integs.Size(); k++)
+         for (size_t k = 0; k < boundary_constraint_integs.size(); k++)
          {
-            if (boundary_constraint_integs_marker[k] == NULL)
+            if (boundary_constraint_integs_marker[k].IsEmpty())
             {
                bdr_attr_marker = 1;
                break;
             }
-            Array<int> &bdr_marker = *boundary_constraint_integs_marker[k];
+            Array<int> &bdr_marker = boundary_constraint_integs_marker[k];
             MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
                         "invalid boundary marker for boundary face integrator #"
                         << k << ", counting from zero");
@@ -169,7 +155,7 @@ void Hybridization::ConstructC()
             }
          }
 
-         for (int i = 0; i < fes->GetNBE(); i++)
+         for (int i = 0; i < fes.GetNBE(); i++)
          {
             const int bdr_attr = mesh->GetBdrAttribute(i);
             if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
@@ -186,17 +172,17 @@ void Hybridization::ConstructC()
                vdofs[j] = o1 + j;
             }
             int iface = mesh->GetBdrElementFaceIndex(i);
-            c_fes->GetFaceVDofs(iface, c_vdofs);
-            face_el = c_fes->GetFaceElement(iface);
-            fe1 = fes -> GetFE (FTr -> Elem1No);
+            c_fes.GetFaceVDofs(iface, c_vdofs);
+            face_el = c_fes.GetFaceElement(iface);
+            fe1 = fes.GetFE(FTr -> Elem1No);
             // The fe2 object is really a dummy and not used on the boundaries,
             // but we can't dereference a NULL pointer, and we don't want to
             // actually make a fake element.
             fe2 = fe1;
-            for (int k = 0; k < boundary_constraint_integs.Size(); k++)
+            for (size_t k = 0; k < boundary_constraint_integs.size(); k++)
             {
-               if (boundary_constraint_integs_marker[k] &&
-                   (*boundary_constraint_integs_marker[k])[bdr_attr-1] == 0) { continue; }
+               if (boundary_constraint_integs_marker[k].Size() &&
+                   boundary_constraint_integs_marker[k][bdr_attr-1] == 0) { continue; }
 
                boundary_constraint_integs[k]->AssembleFaceMatrix(*face_el, *fe1, *fe2, *FTr,
                                                                  elmat);
@@ -220,8 +206,8 @@ void Hybridization::ConstructC()
             {
                FTr = pmesh->GetFaceElementTransformations(face_no);
                MFEM_ASSERT(FTr->Elem2No < 0, "");
-               face_fe = c_fes->GetFaceElement(face_no);
-               c_fes->GetFaceVDofs(face_no, c_vdofs);
+               face_fe = c_fes.GetFaceElement(face_no);
+               c_fes.GetFaceVDofs(face_no, c_vdofs);
             }
             else
             {
@@ -242,7 +228,7 @@ void Hybridization::ConstructC()
             {
                vdofs[j] = o1 + j;
             }
-            fe = fes->GetFE(FTr->Elem1No);
+            fe = fes.GetFE(FTr->Elem1No);
             c_bfi->AssembleFaceMatrix(*face_fe, *fe, *fe, *FTr, elmat);
             // zero-out small elements in elmat
             elmat.Threshold(mtol * elmat.MaxMaxNorm());
@@ -272,12 +258,12 @@ void Hybridization::ConstructC()
                                Ct->GetI(), Ct_J.GetData(), Ct->GetData(),
                                Ct_rows, c_pfes->GetDofOffsets());
             Ct_J.DeleteAll();
-            pC = pCt.Transpose();
+            pC.reset(pCt.Transpose());
          }
          if (pmesh->Nonconforming())
          {
             // TODO - Construct P_pc directly in the pH format
-            P_pc = c_pfes->GetPartialConformingInterpolation();
+            P_pc.reset(c_pfes->GetPartialConformingInterpolation());
          }
       }
 #endif
@@ -294,15 +280,21 @@ void Hybridization::Init(const Array<int> &ess_tdof_list)
 {
    if (Ct) { return; }
 
+   if (ext)
+   {
+      ext->Init(ess_tdof_list);
+      return;
+   }
+
    // count the number of dofs in the discontinuous version of fes:
-   const int NE = fes->GetNE();
+   const int NE = fes.GetNE();
    Array<int> vdofs;
    int num_hat_dofs = 0;
    hat_offsets.SetSize(NE+1);
    hat_offsets[0] = 0;
    for (int i = 0; i < NE; i++)
    {
-      fes->GetElementVDofs(i, vdofs);
+      fes.GetElementVDofs(i, vdofs);
       num_hat_dofs += vdofs.Size();
       hat_offsets[i+1] = num_hat_dofs;
    }
@@ -318,7 +310,7 @@ void Hybridization::Init(const Array<int> &ess_tdof_list)
       C->PrintMatlab(C_file);
       delete C;
 
-      const SparseMatrix *P = fes->GetConformingProlongation();
+      const SparseMatrix *P = fes.GetConformingProlongation();
       if (P)
       {
          std::ofstream P_file("P_matrix.txt");
@@ -333,11 +325,11 @@ void Hybridization::Init(const Array<int> &ess_tdof_list)
    hat_dofs_marker.SetSize(num_hat_dofs);
    Array<int> free_tdof_marker;
 #ifdef MFEM_USE_MPI
-   ParFiniteElementSpace *pfes = dynamic_cast<ParFiniteElementSpace*>(fes);
+   ParFiniteElementSpace *pfes = dynamic_cast<ParFiniteElementSpace*>(&fes);
    free_tdof_marker.SetSize(pfes ? pfes->TrueVSize() :
-                            fes->GetConformingVSize());
+                            fes.GetConformingVSize());
 #else
-   free_tdof_marker.SetSize(fes->GetConformingVSize());
+   free_tdof_marker.SetSize(fes.GetConformingVSize());
 #endif
    free_tdof_marker = 1;
    for (int i = 0; i < ess_tdof_list.Size(); i++)
@@ -348,38 +340,38 @@ void Hybridization::Init(const Array<int> &ess_tdof_list)
 #ifdef MFEM_USE_MPI
    if (!pfes)
    {
-      const SparseMatrix *cP = fes->GetConformingProlongation();
+      const SparseMatrix *cP = fes.GetConformingProlongation();
       if (!cP)
       {
          free_vdofs_marker.MakeRef(free_tdof_marker);
       }
       else
       {
-         free_vdofs_marker.SetSize(fes->GetVSize());
+         free_vdofs_marker.SetSize(fes.GetVSize());
          cP->BooleanMult(free_tdof_marker, free_vdofs_marker);
       }
    }
    else
    {
       HypreParMatrix *P = pfes->Dof_TrueDof_Matrix();
-      free_vdofs_marker.SetSize(fes->GetVSize());
+      free_vdofs_marker.SetSize(fes.GetVSize());
       P->BooleanMult(1, free_tdof_marker, 0, free_vdofs_marker);
    }
 #else
-   const SparseMatrix *cP = fes->GetConformingProlongation();
+   const SparseMatrix *cP = fes.GetConformingProlongation();
    if (!cP)
    {
       free_vdofs_marker.MakeRef(free_tdof_marker);
    }
    else
    {
-      free_vdofs_marker.SetSize(fes->GetVSize());
+      free_vdofs_marker.SetSize(fes.GetVSize());
       cP->BooleanMult(free_tdof_marker, free_vdofs_marker);
    }
 #endif
    for (int i = 0; i < NE; i++)
    {
-      fes->GetElementVDofs(i, vdofs);
+      fes.GetElementVDofs(i, vdofs);
       FiniteElementSpace::AdjustVDofs(vdofs);
       for (int j = 0; j < vdofs.Size(); j++)
       {
@@ -442,18 +434,18 @@ void Hybridization::Init(const Array<int> &ess_tdof_list)
 #undef MFEM_DEBUG_HERE
 #endif
 
-   Af_data = new real_t[Af_offsets[NE]];
-   Af_ipiv = new int[Af_f_offsets[NE]];
+   Af_data.SetSize(Af_offsets[NE]);
+   Af_ipiv.SetSize(Af_f_offsets[NE]);
 
 #ifdef MFEM_DEBUG
    // check that Ref = 0
-   const SparseMatrix *R = fes->GetRestrictionMatrix();
+   const SparseMatrix *R = fes.GetRestrictionMatrix();
    if (!R) { return; }
-   Array<int> vdof_marker(fes->GetVSize()); // 0 - f, 1 - e
+   Array<int> vdof_marker(fes.GetVSize()); // 0 - f, 1 - e
    vdof_marker = 0;
    for (int i = 0; i < NE; i++)
    {
-      fes->GetElementVDofs(i, vdofs);
+      fes.GetElementVDofs(i, vdofs);
       FiniteElementSpace::AdjustVDofs(vdofs);
       for (int j = 0; j < vdofs.Size(); j++)
       {
@@ -519,11 +511,17 @@ void Hybridization::GetBDofs(int el, int &num_idofs, Array<int> &b_dofs) const
 
 void Hybridization::AssembleMatrix(int el, const DenseMatrix &A)
 {
+   if (ext)
+   {
+      ext->AssembleMatrix(el, A);
+      return;
+   }
+
    Array<int> i_dofs, b_dofs;
 
    GetIBDofs(el, i_dofs, b_dofs);
 
-   DenseMatrix A_ii(Af_data + Af_offsets[el], i_dofs.Size(), i_dofs.Size());
+   DenseMatrix A_ii(&Af_data[Af_offsets[el]], i_dofs.Size(), i_dofs.Size());
    DenseMatrix A_ib(A_ii.Data() + i_dofs.Size()*i_dofs.Size(),
                     i_dofs.Size(), b_dofs.Size());
    DenseMatrix A_bi(A_ib.Data() + i_dofs.Size()*b_dofs.Size(),
@@ -557,12 +555,32 @@ void Hybridization::AssembleMatrix(int el, const DenseMatrix &A)
    }
 }
 
+void Hybridization::AssembleElementMatrices(const class DenseTensor &el_mats)
+{
+   if (ext)
+   {
+      ext->AssembleElementMatrices(el_mats);
+      return;
+   }
+
+   for (int e = 0; e < el_mats.SizeK(); ++e)
+   {
+      AssembleMatrix(e, el_mats(e));
+   }
+}
+
 void Hybridization::AssembleBdrMatrix(int bdr_el, const DenseMatrix &A)
 {
+   if (ext)
+   {
+      ext->AssembleBdrMatrix(bdr_el, A);
+      return;
+   }
+
    // Not tested.
 #ifdef MFEM_DEBUG
    Array<int> vdofs, bvdofs;
-   fes->GetBdrElementVDofs(bdr_el, bvdofs);
+   fes.GetBdrElementVDofs(bdr_el, bvdofs);
 #endif
 
    int el;
@@ -570,19 +588,19 @@ void Hybridization::AssembleBdrMatrix(int bdr_el, const DenseMatrix &A)
    Array<int> i_dofs, b_dofs, e2f;
 
    {
-      int info, vdim = fes->GetVDim();
+      int info, vdim = fes.GetVDim();
       Array<int> lvdofs;
-      Mesh *mesh = fes->GetMesh();
+      Mesh *mesh = fes.GetMesh();
       mesh->GetBdrElementAdjacentElement(bdr_el, el, info);
       e2f.SetSize(hat_offsets[el+1]-hat_offsets[el], -1);
       lvdofs.Reserve(A.Height());
-      fes->FEColl()->SubDofOrder(mesh->GetElementBaseGeometry(el),
-                                 mesh->Dimension()-1, info, lvdofs);
+      fes.FEColl()->SubDofOrder(mesh->GetElementBaseGeometry(el),
+                                mesh->Dimension()-1, info, lvdofs);
       // Convert local element dofs to local element vdofs.
       Ordering::DofsToVDofs<Ordering::byNODES>(e2f.Size()/vdim, vdim, lvdofs);
       MFEM_ASSERT(lvdofs.Size() == A.Height(), "internal error");
 #ifdef MFEM_DEBUG
-      fes->GetElementVDofs(el, vdofs);
+      fes.GetElementVDofs(el, vdofs);
       for (int i = 0; i < lvdofs.Size(); i++)
       {
          int bd = lvdofs[i];
@@ -601,7 +619,7 @@ void Hybridization::AssembleBdrMatrix(int bdr_el, const DenseMatrix &A)
 
    GetIBDofs(el, i_dofs, b_dofs);
 
-   DenseMatrix A_ii(Af_data + Af_offsets[el], i_dofs.Size(), i_dofs.Size());
+   DenseMatrix A_ii(&Af_data[Af_offsets[el]], i_dofs.Size(), i_dofs.Size());
    DenseMatrix A_ib(A_ii.Data() + i_dofs.Size()*i_dofs.Size(),
                     i_dofs.Size(), b_dofs.Size());
    DenseMatrix A_bi(A_ib.Data() + i_dofs.Size()*b_dofs.Size(),
@@ -647,15 +665,24 @@ void Hybridization::AssembleBdrMatrix(int bdr_el, const DenseMatrix &A)
 
 void Hybridization::ComputeH()
 {
+   if (ext)
+   {
+      ext->ConstructH();
+      return;
+   }
+
    const int skip_zeros = 1;
    Array<int> c_dof_marker(Ct->Width());
    Array<int> b_dofs, c_dofs;
-   const int NE = fes->GetNE();
+   const int NE = fes.GetNE();
    DenseMatrix Cb_t, Sb_inv_Cb_t, Hb;
 #ifndef MFEM_USE_MPI
-   H = new SparseMatrix(Ct->Width());
+   H.reset(new SparseMatrix(Ct->Width()));
 #else
-   H = pC ? NULL : new SparseMatrix(Ct->Width());
+   if (!pC)
+   {
+      H.reset(new SparseMatrix(Ct->Width()));
+   }
    // V = Sb^{-1} Cb^T, for parallel non-conforming meshes
    SparseMatrix *V = pC ? new SparseMatrix(Ct->Height(), Ct->Width()) : NULL;
 #endif
@@ -667,7 +694,7 @@ void Hybridization::ComputeH()
       int i_dofs_size;
       GetBDofs(el, i_dofs_size, b_dofs);
 
-      LUFactors LU_ii(Af_data + Af_offsets[el], Af_ipiv + Af_f_offsets[el]);
+      LUFactors LU_ii(&Af_data[Af_offsets[el]], Af_ipiv + Af_f_offsets[el]);
       real_t *A_ib_data = LU_ii.data + i_dofs_size*i_dofs_size;
       real_t *A_bi_data = A_ib_data + i_dofs_size*b_dofs.Size();
       LUFactors LU_bb(A_bi_data + i_dofs_size*b_dofs.Size(),
@@ -737,7 +764,7 @@ void Hybridization::ComputeH()
 #ifndef MFEM_USE_MPI
    H->Finalize(skip_zeros, fix_empty_rows);
 #else
-   ParFiniteElementSpace *c_pfes = dynamic_cast<ParFiniteElementSpace*>(c_fes);
+   ParFiniteElementSpace *c_pfes = dynamic_cast<ParFiniteElementSpace*>(&c_fes);
    if (!pC)
    {
       H->Finalize(skip_zeros, fix_empty_rows);
@@ -745,12 +772,11 @@ void Hybridization::ComputeH()
 
       OperatorHandle pP(pH.Type()), dH(pH.Type());
       // TODO - construct P_pc / Dof_TrueDof_Matrix directly in the pH format
-      pP.ConvertFrom(P_pc ? P_pc : c_pfes->Dof_TrueDof_Matrix());
+      pP.ConvertFrom(P_pc ? P_pc.get() : c_pfes->Dof_TrueDof_Matrix());
       dH.MakeSquareBlockDiag(c_pfes->GetComm(),c_pfes->GlobalVSize(),
-                             c_pfes->GetDofOffsets(), H);
+                             c_pfes->GetDofOffsets(), H.get());
       pH.MakePtAP(dH, pP);
-      delete H;
-      H = NULL;
+      H.reset();
    }
    else
    {
@@ -758,7 +784,7 @@ void Hybridization::ComputeH()
       V->Finalize();
       Array<HYPRE_BigInt> V_J(V->NumNonZeroElems());
       MFEM_ASSERT(c_pfes, "");
-      const int c_vsize = c_fes->GetVSize();
+      const int c_vsize = c_fes.GetVSize();
       HYPRE_BigInt c_ldof_offset = c_pfes->GetMyDofOffset();
       const HYPRE_BigInt *c_face_nbr_glob_ldof = c_pfes->GetFaceNbrGlobalDofMap();
       int *J = V->GetJ();
@@ -779,11 +805,11 @@ void Hybridization::ComputeH()
          // safely delete V_J and V:
          V_J.DeleteAll();
          delete V;
-         lpH = ParMult(pC, &pV);
+         lpH = ParMult(pC.get(), &pV);
       }
       OperatorHandle pP(pH.Type()), plpH(pH.Type());
       // TODO - construct P_pc directly in the pH format
-      pP.ConvertFrom(P_pc);
+      pP.ConvertFrom(P_pc.get());
       plpH.ConvertFrom(lpH);
       MFEM_VERIFY(pH.Type() != Operator::PETSC_MATIS, "To be implemented");
       pH.MakePtAP(plpH, pP);
@@ -806,18 +832,18 @@ void Hybridization::MultAfInv(const Vector &b, const Vector &lambda, Vector &bf,
 {
    // b1 = Rf^t b (assuming that Ref = 0)
    Vector b1;
-   const SparseMatrix *R = fes->GetRestrictionMatrix();
+   const SparseMatrix *R = fes.GetRestrictionMatrix();
    if (!R)
    {
       b1.SetDataAndSize(b.GetData(), b.Size());
    }
    else
    {
-      b1.SetSize(fes->GetVSize());
+      b1.SetSize(fes.GetVSize());
       R->MultTranspose(b, b1);
    }
 
-   const int NE = fes->GetMesh()->GetNE();
+   const int NE = fes.GetMesh()->GetNE();
    Array<int> vdofs, i_dofs, b_dofs;
    Vector el_vals, bf_i, i_vals, b_vals;
    bf.SetSize(hat_offsets[NE]);
@@ -825,7 +851,7 @@ void Hybridization::MultAfInv(const Vector &b, const Vector &lambda, Vector &bf,
    {
 #ifdef MFEM_USE_MPI
       ParFiniteElementSpace *c_pfes =
-         dynamic_cast<ParFiniteElementSpace*>(c_fes);
+         dynamic_cast<ParFiniteElementSpace*>(&c_fes);
       if (!c_pfes)
       {
          Ct->Mult(lambda, bf);
@@ -833,7 +859,7 @@ void Hybridization::MultAfInv(const Vector &b, const Vector &lambda, Vector &bf,
       else
       {
          Vector L(c_pfes->GetVSize());
-         (P_pc ? P_pc : c_pfes->GetProlongationMatrix())->Mult(lambda, L);
+         (P_pc ? P_pc.get() : c_pfes->GetProlongationMatrix())->Mult(lambda, L);
          pC ? pC->MultTranspose(L, bf) : Ct->Mult(L, bf);
       }
 #else
@@ -845,7 +871,7 @@ void Hybridization::MultAfInv(const Vector &b, const Vector &lambda, Vector &bf,
    vdof_marker = false;
    for (int i = 0; i < NE; i++)
    {
-      fes->GetElementVDofs(i, vdofs);
+      fes.GetElementVDofs(i, vdofs);
       b1.GetSubVector(vdofs, el_vals);
       for (int j = 0; j < vdofs.Size(); j++)
       {
@@ -863,7 +889,9 @@ void Hybridization::MultAfInv(const Vector &b, const Vector &lambda, Vector &bf,
       el_vals.GetSubVector(i_dofs, i_vals);
       el_vals.GetSubVector(b_dofs, b_vals);
 
-      LUFactors LU_ii(Af_data + Af_offsets[i], Af_ipiv + Af_f_offsets[i]);
+      real_t *Af_data_ptr = const_cast<real_t*>(&Af_data[Af_offsets[i]]);
+      int *Af_ipiv_ptr = const_cast<int*>(&Af_ipiv[Af_f_offsets[i]]);
+      LUFactors LU_ii(Af_data_ptr, Af_ipiv_ptr);
       real_t *U_ib = LU_ii.data + i_dofs.Size()*i_dofs.Size();
       real_t *L_bi = U_ib + i_dofs.Size()*b_dofs.Size();
       LUFactors LU_bb(L_bi + b_dofs.Size()*i_dofs.Size(),
@@ -884,13 +912,19 @@ void Hybridization::MultAfInv(const Vector &b, const Vector &lambda, Vector &bf,
 
 void Hybridization::ReduceRHS(const Vector &b, Vector &b_r) const
 {
+   if (ext)
+   {
+      ext->ReduceRHS(b, b_r);
+      return;
+   }
+
    // bf = Af^{-1} Rf^t b
    Vector bf;
    MultAfInv(b, b, bf, 0);
 
    // b_r = Cf bf
 #ifdef MFEM_USE_MPI
-   ParFiniteElementSpace *c_pfes = dynamic_cast<ParFiniteElementSpace*>(c_fes);
+   ParFiniteElementSpace *c_pfes = dynamic_cast<ParFiniteElementSpace*>(&c_fes);
    if (!c_pfes)
    {
       b_r.SetSize(Ct->Width());
@@ -908,7 +942,7 @@ void Hybridization::ReduceRHS(const Vector &b, Vector &b_r) const
          Ct->MultTranspose(bf, bl);
       }
       b_r.SetSize(pH.Ptr()->Height());
-      (P_pc ? P_pc : c_pfes->GetProlongationMatrix())->MultTranspose(bl, b_r);
+      (P_pc ? P_pc.get() : c_pfes->GetProlongationMatrix())->MultTranspose(bl, b_r);
    }
 #else
    b_r.SetSize(Ct->Width());
@@ -919,28 +953,34 @@ void Hybridization::ReduceRHS(const Vector &b, Vector &b_r) const
 void Hybridization::ComputeSolution(const Vector &b, const Vector &sol_r,
                                     Vector &sol) const
 {
-   // bf = Af^{-1} ( Rf^t - Cf^t sol_r )
+   if (ext)
+   {
+      ext->ComputeSolution(b, sol_r, sol);
+      return;
+   }
+
+   // bf = Af^{-1} ( Rf^t b - Cf^t sol_r )
    Vector bf;
    MultAfInv(b, sol_r, bf, 1);
 
    // sol = Rf bf
    GridFunction s;
-   const SparseMatrix *R = fes->GetRestrictionMatrix();
+   const SparseMatrix *R = fes.GetRestrictionMatrix();
    if (!R)
    {
-      MFEM_ASSERT(sol.Size() == fes->GetVSize(), "");
-      s.MakeRef(fes, sol, 0);
+      MFEM_ASSERT(sol.Size() == fes.GetVSize(), "");
+      s.MakeRef(&fes, sol, 0);
    }
    else
    {
-      s.SetSpace(fes);
+      s.SetSpace(&fes);
       R->MultTranspose(sol, s);
    }
-   const int NE = fes->GetMesh()->GetNE();
+   const int NE = fes.GetMesh()->GetNE();
    Array<int> vdofs;
    for (int i = 0; i < NE; i++)
    {
-      fes->GetElementVDofs(i, vdofs);
+      fes.GetElementVDofs(i, vdofs);
       for (int j = hat_offsets[i]; j < hat_offsets[i+1]; j++)
       {
          if (hat_dofs_marker[j] == 1) { continue; } // skip essential b.c.
@@ -957,11 +997,15 @@ void Hybridization::ComputeSolution(const Vector &b, const Vector &sol_r,
 
 void Hybridization::Reset()
 {
-   delete H;
-   H = NULL;
+   H.reset();
 #ifdef MFEM_USE_MPI
    pH.Clear();
 #endif
+   if (ext) { ext->Reset(); }
 }
 
+// Set to default in cpp file because of use of incomplete type
+// (HybridizationExtension) in unique_ptr member data.
+Hybridization::~Hybridization() = default;
+
 }

fem/hybridization.hpp

@@ -15,10 +15,26 @@
 #include "../config/config.hpp"
 #include "fespace.hpp"
 #include "bilininteg.hpp"
+#include <memory>
 
 namespace mfem
 {
 
+template <typename T>
+std::shared_ptr<T> Owning(T *t) { return std::shared_ptr<T>(t); }
+
+template <typename T>
+std::shared_ptr<T> NonOwning(T *t)
+{
+   return std::shared_ptr<T>(t, [](T*) {});
+}
+
+template <typename T>
+std::shared_ptr<T> OptionallyOwning(T *t, bool own)
+{
+   return own ? Owning(t) : NonOwning(t);
+}
+
 /** @brief Auxiliary class Hybridization, used to implement BilinearForm
     hybridization.
 
@@ -60,34 +76,43 @@ namespace mfem
     (weak) continuity constraints between neighboring elements. */
 class Hybridization
 {
+   friend class HybridizationExtension;
 protected:
-   FiniteElementSpace *fes, *c_fes;
-   BilinearFormIntegrator *c_bfi;
-
-   /// Set of constraint boundary face integrators to be applied.
-   Array<BilinearFormIntegrator*> boundary_constraint_integs;
-   Array<Array<int>*>             boundary_constraint_integs_marker;
-   /// Indicates if the boundary_constraint_integs integrators are owned externally
-   bool extern_bdr_constr_integs;
-
-   SparseMatrix *Ct, *H;
+   FiniteElementSpace &fes; ///< The finite element space.
+   FiniteElementSpace &c_fes; ///< The constraint finite element space.
+   /// Extension for device execution.
+   std::unique_ptr<class HybridizationExtension> ext;
+   /// The constraint integrator.
+   std::unique_ptr<BilinearFormIntegrator> c_bfi;
+   /// The constraint boundary face integrators.
+   std::vector<std::shared_ptr<BilinearFormIntegrator>> boundary_constraint_integs;
+   /// Boundary markers for constraint face integrators.
+   std::vector<Array<int>> boundary_constraint_integs_marker;
+   /// The constraint matrix.
+   std::unique_ptr<SparseMatrix> Ct;
+   /// The Schur complement system for the Lagrange multiplier.
+   std::unique_ptr<SparseMatrix> H;
 
    Array<int> hat_offsets, hat_dofs_marker;
    Array<int> Af_offsets, Af_f_offsets;
-   real_t *Af_data;
-   int *Af_ipiv;
+   Array<real_t> Af_data;
+   Array<int> Af_ipiv;
 
 #ifdef MFEM_USE_MPI
-   HypreParMatrix *pC, *P_pc; // for parallel non-conforming meshes
+   std::unique_ptr<HypreParMatrix> pC, P_pc; // for parallel non-conforming meshes
    OperatorHandle pH;
 #endif
 
+   /// Construct the constraint matrix.
    void ConstructC();
 
+   /// Returns the local indices of the i-dofs and b-dofs of element @a el.
    void GetIBDofs(int el, Array<int> &i_dofs, Array<int> &b_dofs) const;
 
+   /// Returns global indices of the b-dofs of element @a el.
    void GetBDofs(int el, int &num_idofs, Array<int> &b_dofs) const;
 
+   /// Construct the Schur complement system.
    void ComputeH();
 
    // Compute depending on mode:
@@ -101,43 +126,58 @@ protected:
                   int mode) const;
 
 public:
-   /// Constructor
+   /// Constructor.
    Hybridization(FiniteElementSpace *fespace, FiniteElementSpace *c_fespace);
-   /// Destructor
+
+   /// Destructor.
    ~Hybridization();
 
-   /** Set the integrator that will be used to construct the constraint matrix
-       C. The Hybridization object assumes ownership of the integrator, i.e. it
-       will delete the integrator when destroyed. */
+   /// Turns on device execution.
+   void EnableDeviceExecution();
+
+   /// @brief Set the integrator that will be used to construct the constraint
+   /// matrix C.
+   ///
+   /// The Hybridization object assumes ownership of the integrator, i.e. it
+   /// will delete the integrator when destroyed.
    void SetConstraintIntegrator(BilinearFormIntegrator *c_integ)
-   { delete c_bfi; c_bfi = c_integ; }
+   { c_bfi.reset(c_integ); }
 
-   /** Add the boundary face integrator that will be used to construct the
-       constraint matrix C. The Hybridization object assumes ownership of the
-       integrator, i.e. it will delete the integrator when destroyed. */
-   void AddBdrConstraintIntegrator(BilinearFormIntegrator *c_integ)
+   /// @brief Add a boundary face integrator that will be used to construct the
+   /// constraint matrix C.
+   ///
+   /// The integrator will apply to the boundaries specified using the marker
+   /// array @a bdr_marker. If @a bdr_marker is empty (its default value) then
+   /// the integrator will be applied on all boundaries.
+   void AddBdrConstraintIntegrator(
+      const std::shared_ptr<BilinearFormIntegrator> &c_integ,
+      const Array<int> &bdr_marker = Array<int>())
    {
-      boundary_constraint_integs.Append(c_integ);
-      boundary_constraint_integs_marker.Append(
-         NULL); // NULL marker means apply everywhere
+      boundary_constraint_integs.push_back(c_integ);
+      boundary_constraint_integs_marker.push_back(bdr_marker);
    }
+
+   /// @brief Add the boundary face integrator that will be used to construct
+   /// the constraint matrix C.
+   ///
+   /// If @a own is true (its default value), then the Hybridization object
+   /// assumes ownership of the integrator.
+   ///
+   /// @sa AddBdrConstraintIntegrator().
    void AddBdrConstraintIntegrator(BilinearFormIntegrator *c_integ,
-                                   Array<int> &bdr_marker)
+                                   const Array<int> &bdr_marker = Array<int>(),
+                                   bool own = true)
    {
-      boundary_constraint_integs.Append(c_integ);
-      boundary_constraint_integs_marker.Append(&bdr_marker);
+      AddBdrConstraintIntegrator(OptionallyOwning(c_integ, own), bdr_marker);
    }
 
-   /// Access all integrators added with AddBdrConstraintIntegrator().
-   Array<BilinearFormIntegrator*> *GetBCBFI() { return &boundary_constraint_integs; }
+   /// Access the integrators added with AddBdrConstraintIntegrator().
+   const std::vector<std::shared_ptr<BilinearFormIntegrator>>
+   &GetBdrConstraintIntegrators() const { return boundary_constraint_integs; }
 
-   /// Access all boundary markers added with AddBdrConstraintIntegrator().
-   /** If no marker was specified when the integrator was added, the
-       corresponding pointer (to Array<int>) will be NULL. */
-   Array<Array<int>*> *GetBCBFI_Marker() { return &boundary_constraint_integs_marker; }
-
-   /// Indicate that boundary constraint integrators are not owned
-   void UseExternalBdrConstraintIntegrators() { extern_bdr_constr_integs = true; }
+   /// Access the boundary markers added with AddBdrConstraintIntegrator().
+   const std::vector<Array<int>> &GetBdrConstraintIntegratorMarker() const
+   { return boundary_constraint_integs_marker; }
 
    /// Prepare the Hybridization object for assembly.
    void Init(const Array<int> &ess_tdof_list);
@@ -145,6 +185,9 @@ public:
    /// Assemble the element matrix A into the hybridized system matrix.
    void AssembleMatrix(int el, const DenseMatrix &A);
 
+   /// Assemble all of the element matrices given in the form of a DenseTensor.
+   void AssembleElementMatrices(const class DenseTensor &el_mats);
+
    /// Assemble the boundary element matrix A into the hybridized system matrix.
    void AssembleBdrMatrix(int bdr_el, const DenseMatrix &A);
 
@@ -158,30 +201,34 @@ public:
    /// Return the parallel hybridized matrix.
    HypreParMatrix &GetParallelMatrix() { return *pH.Is<HypreParMatrix>(); }
 
-   /** @brief Return the parallel hybridized matrix in the format specified by
-       SetOperatorType(). */
+   /// @brief Return the parallel hybridized matrix in the format specified by
+   /// SetOperatorType().
    void GetParallelMatrix(OperatorHandle &H_h) const { H_h = pH; }
 
    /// Set the operator type id for the parallel hybridized matrix/operator.
    void SetOperatorType(Operator::Type tid) { pH.SetType(tid); }
 #endif
 
-   /** Perform the reduction of the given r.h.s. vector, b, to a r.h.s vector,
-       b_r, for the hybridized system. */
+   /// @brief Perform the reduction of the given right-hand side @a b to a
+   /// right-hand side vector @a b_r for the hybridized system.
    void ReduceRHS(const Vector &b, Vector &b_r) const;
 
-   /** Reconstruct the solution of the original system, sol, from solution of
-       the hybridized system, sol_r, and the original r.h.s. vector, b.
-       It is assumed that the vector sol has the right essential b.c. */
+   /// @brief Reconstruct the solution of the original system @a sol from
+   /// solution of the hybridized system @a sol_r and the original right-hand
+   /// side @a b.
+   ///
+   /// It is assumed that the vector sol has the correct essential boundary
+   /// conditions.
    void ComputeSolution(const Vector &b, const Vector &sol_r,
                         Vector &sol) const;
 
-   /** @brief Destroy the current hybridization matrix while preserving the
-       computed constraint matrix and the set of essential true dofs. After
-       Reset(), a new hybridized matrix can be assembled via AssembleMatrix()
-       and Finalize(). The Mesh and FiniteElementSpace objects are assumed to be
-       un-modified. If that is not the case, a new Hybridization object must be
-       created. */
+   /// @brief Destroy the current hybridization matrix while preserving the
+   /// computed constraint matrix and the set of essential true dofs.
+   ///
+   /// After Reset(), a new hybridized matrix can be assembled via
+   /// AssembleMatrix() and Finalize(). The Mesh and FiniteElementSpace objects
+   /// are assumed to be unmodified. If that is not the case, a new
+   /// Hybridization object must be created.
    void Reset();
 };
 

fem/hybridization_ext.hpp

@@ -0,0 +1,119 @@
+// Copyright (c) 2010-2024, 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_HYBRIDIZATION_EXT
+#define MFEM_HYBRIDIZATION_EXT
+
+#include "../config/config.hpp"
+#include "../general/array.hpp"
+#include "../linalg/vector.hpp"
+
+namespace mfem
+{
+
+/// @brief Extension class supporting Hybridization on device (GPU).
+///
+/// Similar to BilinearFormExtension and LinearFormExtension, this extension
+/// class provides device execution capabilities for the Hybridization class.
+///
+/// As with the other extension classes, a limitation of this class is that it
+/// requires meshes consisting only of tensor-product elements, and finite
+/// element spaces without variable polynomial degrees.
+class HybridizationExtension
+{
+   friend class Hybridization;
+public:
+   enum DofType : char
+   {
+      ESSENTIAL,
+      BOUNDARY,
+      INTERIOR
+   };
+protected:
+   class Hybridization &h; ///< The associated Hybridization object.=
+   int num_hat_dofs; ///< Number of Lagrange multipliers.
+   mutable Vector tmp1, tmp2; ///< Temporary vectors.
+
+   Array<int> hat_dof_gather_map;
+   Array<DofType> hat_dof_marker;
+
+   Array<int> el_to_face;
+   Array<int> face_to_el;
+   Vector Ct_mat; ///< Constraint matrix (transposed) stored element-wise.
+
+   Array<int> idofs, bdofs;
+
+   Vector Ahat, Ahat_ii, Ahat_ib, Ahat_bi, Ahat_bb;
+   Array<int> Ahat_ii_piv, Ahat_bb_piv;
+
+public:
+   /// Construct the constraint matrix.
+   void ConstructC();
+
+   template <int MID, int MBD>
+   void FactorElementMatrices(Vector &AhatInvCt_mat);
+
+   /// Form the Schur complement matix $H$.
+   void ConstructH();
+
+   /// Compute the action of C^t x.
+   void MultCt(const Vector &x, Vector &y) const;
+
+   /// Compute the action of C x.
+   void MultC(const Vector &x, Vector &y) const;
+
+   /// @brief Assemble the element matrix A into the hybridized system matrix.
+   ///
+   /// @warning Using the interface will be very slow. AssembleElementMatrices()
+   /// should be used instead.
+   void AssembleMatrix(int el, const class DenseMatrix &elmat);
+
+   /// @brief Assemble the boundary element matrix A into the hybridized system
+   /// matrix.
+   ///
+   /// @warning Using the interface will be very slow. AssembleElementMatrices()
+   /// should be used instead.
+   void AssembleBdrMatrix(int bdr_el, const class DenseMatrix &elmat);
+
+   /// Invert and store the element matrices Ahat.
+   void AssembleElementMatrices(const class DenseTensor &el_mats);
+
+   /// Apply the action of R mapping from "hat DOFs" to T-vector
+   void MultR(const Vector &b, Vector &b_hat) const;
+
+   /// Apply the action of R^t mapping into the "hat DOF" space.
+   void MultRt(const Vector &b, Vector &b_hat) const;
+
+   /// Apply the elementwise A_hat^{-1}.
+   void MultAhatInv(Vector &x) const;
+
+   /// Constructor.
+   HybridizationExtension(class Hybridization &hybridization_);
+
+   /// Prepare for assembly; form the constraint matrix.
+   void Init(const Array<int> &ess_tdof_list);
+
+   /// @brief Given a right-hand side on the original space, compute the
+   /// corresponding right-hand side for the Lagrange multipliers.
+   void ReduceRHS(const Vector &b, Vector &b_r) const;
+
+   /// @brief Given Lagrange multipliers @a sol_r and the original right-hand
+   /// side @a b, recover the solution @a sol on the original finite element
+   /// space.
+   void ComputeSolution(const Vector &b, const Vector &sol_r, Vector &sol) const;
+
+   /// Destroys the stored element matrices.
+   void Reset() { Ahat = 0.0; }
+};
+
+}
+
+#endif

fem/integ/bilininteg_hdiv_ea.cpp

@@ -0,0 +1,336 @@
+// Copyright (c) 2010-2024, 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 "../../general/forall.hpp"
+#include "../bilininteg.hpp"
+#include "../gridfunc.hpp"
+
+namespace mfem
+{
+
+// For H(div) mass, Bo and Bc are the basis evaluation operators, and the
+// pa_data corresponds to a (potentially symmetric) matrix coefficient.
+// coeff_dim must be 3 or 4 depending on symmetry.
+//
+// For div-div, Bc is the derivative evaluation operator, and pa_data
+// corresponds to a scalar coefficient. coeff_dim must be 1.
+//
+// These two integrators are distinguished using coeff_dim.
+template<int T_D1D = 0, int T_Q1D = 0>
+static void EAHdivAssemble2D(const int NE,
+                             const Array<real_t> &Bo_,
+                             const Array<real_t> &Bc_,
+                             const int coeff_dim,
+                             const Vector &pa_data,
+                             Vector &ea_data,
+                             const bool add,
+                             const int d1d = 0,
+                             const int q1d = 0)
+{
+   const int D1D = T_D1D ? T_D1D : d1d;
+   const int Q1D = T_Q1D ? T_Q1D : q1d;
+   MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().HDIV_MAX_D1D, "");
+   MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().HDIV_MAX_Q1D, "");
+   const int NDOF = 2*(D1D-1)*D1D;
+   const auto Bo = Reshape(Bo_.Read(), Q1D, D1D-1);
+   const auto Bc = Reshape(Bc_.Read(), Q1D, D1D);
+   const auto D = Reshape(pa_data.Read(), Q1D, Q1D, coeff_dim, NE);
+   const bool symmetric = (coeff_dim == 3);
+   auto M = Reshape(add ? ea_data.ReadWrite() : ea_data.Write(), NDOF, NDOF, NE);
+   mfem::forall_2D(NE, NDOF, 1, [=] MFEM_HOST_DEVICE (int e)
+   {
+      constexpr int MD1 = T_D1D ? T_D1D : DofQuadLimits::HDIV_MAX_D1D;
+      constexpr int MQ1 = T_Q1D ? T_Q1D : DofQuadLimits::HDIV_MAX_Q1D;
+      // Load Bo and Bc matrices into registers
+      real_t r_Bo[MQ1][MD1];
+      real_t r_Bc[MQ1][MD1];
+      for (int d = 0; d < D1D; d++)
+      {
+         for (int q = 0; q < Q1D; q++)
+         {
+            if (d < D1D - 1) { r_Bo[q][d] = Bo(q,d); }
+            r_Bc[q][d] = Bc(q,d);
+         }
+      }
+      // Store PA data in shared memory
+      MFEM_SHARED real_t s_D[4][MQ1][MQ1];
+      MFEM_FOREACH_THREAD(idx_q, x, Q1D*Q1D)
+      {
+         const int qx = idx_q % Q1D;
+         const int qy = idx_q / Q1D;
+         if (coeff_dim == 1)
+         {
+            const real_t val = D(qx, qy, 0, e);
+            for (int i = 0; i < 4; ++i) { s_D[i][qx][qy] = val; }
+         }
+         else
+         {
+            s_D[0][qx][qy] = D(qx, qy, 0, e);
+            s_D[1][qx][qy] = D(qx, qy, 1, e);
+            s_D[2][qx][qy] = (symmetric) ? s_D[1][qx][qy] : D(qx, qy, 2, e);
+            s_D[3][qx][qy] = (symmetric) ? D(qx, qy, 2, e) : D(qx, qy, 3, e);
+         }
+      }
+      MFEM_SYNC_THREAD;
+      // Assemble (one row per thread)
+      MFEM_FOREACH_THREAD(idx_i, x, NDOF)
+      {
+         const int ic = idx_i / D1D / (D1D-1);
+         const int idx_ii = idx_i % (D1D * (D1D-1));
+         const int ix = (ic == 0) ? idx_ii%D1D : idx_ii%(D1D-1);
+         const int iy = (ic == 0) ? idx_ii/D1D : idx_ii/(D1D-1);
+
+         const real_t (&Bi1)[MQ1][MD1] = (ic == 0) ? r_Bc : r_Bo;
+         const real_t (&Bi2)[MQ1][MD1] = (ic == 0) ? r_Bo : r_Bc;
+
+         for (int idx_j = 0; idx_j < NDOF; ++idx_j)
+         {
+            const int jc = idx_j / (D1D*(D1D-1));
+            const int idx_jj = idx_j % (D1D * (D1D-1));
+            const int jx = (jc == 0) ? idx_jj%D1D : idx_jj%(D1D-1);
+            const int jy = (jc == 0) ? idx_jj/D1D : idx_jj/(D1D-1);
+
+            const real_t (&Bj1)[MQ1][MD1] = (jc == 0) ? r_Bc : r_Bo;
+            const real_t (&Bj2)[MQ1][MD1] = (jc == 0) ? r_Bo : r_Bc;
+
+            real_t val = 0.0;
+            for (int qx = 0; qx < Q1D; ++qx)
+            {
+               for (int qy = 0; qy < Q1D; ++qy)
+               {
+                  const double coeff = s_D[ic + jc*2][qx][qy];
+                  val += coeff*Bi1[qx][ix]*Bi2[qy][iy]*Bj1[qx][jx]*Bj2[qy][jy];
+               }
+            }
+            if (add)
+            {
+               M(idx_i, idx_j, e) += val;
+            }
+            else
+            {
+               M(idx_i, idx_j, e) = val;
+            }
+         }
+      }
+   });
+}
+
+// For H(div) mass, Bo and Bc are the basis evaluation operators, and the
+// pa_data corresponds to a (potentially symmetric) matrix coefficient.
+// coeff_dim must be 6 or 9 depending on symmetry.
+//
+// For div-div, Bc is the derivative evaluation operator, and pa_data
+// corresponds to a scalar coefficient. coeff_dim must be 1.
+//
+// These two integrators are distinguished using coeff_dim.
+template<int T_D1D = 0, int T_Q1D = 0>
+static void EAHdivAssemble3D(const int NE,
+                             const Array<real_t> &Bo_,
+                             const Array<real_t> &Bc_,
+                             const int coeff_dim,
+                             const Vector &pa_data,
+                             Vector &ea_data,
+                             const bool add,
+                             const int d1d = 0,
+                             const int q1d = 0)
+{
+   const int D1D = T_D1D ? T_D1D : d1d;
+   const int Q1D = T_Q1D ? T_Q1D : q1d;
+   MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().HDIV_MAX_D1D, "");
+   MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().HDIV_MAX_Q1D, "");
+   const int NDOF_C = (D1D-1)*(D1D-1)*D1D;
+   const int NDOF = 3*NDOF_C;
+   const auto Bo = Reshape(Bo_.Read(), Q1D, D1D-1);
+   const auto Bc = Reshape(Bc_.Read(), Q1D, D1D);
+   const auto D = Reshape(pa_data.Read(), Q1D, Q1D, Q1D, coeff_dim, NE);
+   const bool symmetric = (coeff_dim == 6);
+   auto M = Reshape(add ? ea_data.ReadWrite() : ea_data.Write(), NDOF, NDOF, NE);
+   mfem::forall_2D(NE, NDOF, 1, [=] MFEM_HOST_DEVICE (int e)
+   {
+      constexpr int MD1 = T_D1D ? T_D1D : DofQuadLimits::HDIV_MAX_D1D;
+      constexpr int MQ1 = T_Q1D ? T_Q1D : DofQuadLimits::HDIV_MAX_Q1D;
+      // Load Bo and Bc matrices into registers
+      real_t r_Bo[MQ1][MD1];
+      real_t r_Bc[MQ1][MD1];
+      for (int d = 0; d < D1D; d++)
+      {
+         for (int q = 0; q < Q1D; q++)
+         {
+            if (d < D1D - 1) { r_Bo[q][d] = Bo(q,d); }
+            r_Bc[q][d] = Bc(q,d);
+         }
+      }
+      // Store PA data in shared memory
+      MFEM_SHARED real_t s_D[9][MQ1][MQ1][MQ1];
+      MFEM_FOREACH_THREAD(idx_q, x, Q1D*Q1D*Q1D)
+      {
+         const int qx = idx_q % Q1D;
+         const int qy = (idx_q / Q1D) % Q1D;
+         const int qz = (idx_q / Q1D) / Q1D;
+         if (coeff_dim == 1)
+         {
+            const real_t val = D(qx,qy,qz,0,e);
+            for (int i = 0; i < 9; ++i) { s_D[i][qx][qy][qz] = val; }
+         }
+         else
+         {
+            s_D[0][qx][qy][qz] = D(qx,qy,qz,0,e);
+            s_D[1][qx][qy][qz] = D(qx,qy,qz,1,e);
+            s_D[2][qx][qy][qz] = D(qx,qy,qz,2,e);
+            s_D[3][qx][qy][qz] = symmetric ? s_D[1][qx][qy][qz] : D(qx,qy,qz,3,e);
+            s_D[4][qx][qy][qz] = symmetric ? D(qx,qy,qz,3,e) : D(qx,qy,qz,4,e);
+            s_D[5][qx][qy][qz] = symmetric ? D(qx,qy,qz,4,e) : D(qx,qy,qz,5,e);
+            s_D[6][qx][qy][qz] = symmetric ? s_D[2][qx][qy][qz] : D(qx,qy,qz,6,e);
+            s_D[7][qx][qy][qz] = symmetric ? s_D[5][qx][qy][qz] : D(qx,qy,qz,7,e);
+            s_D[8][qx][qy][qz] = symmetric ? D(qx,qy,qz,5,e) : D(qx,qy,qz,8,e);
+         }
+      }
+      MFEM_SYNC_THREAD;
+      // Assemble (one row per thread)
+      MFEM_FOREACH_THREAD(idx_i, x, NDOF)
+      {
+         const int ic = idx_i / NDOF_C;
+         const int idx_ii = idx_i % NDOF_C;
+
+         const int nx_i = (ic == 0) ? D1D : D1D-1;
+         const int ny_i = (ic == 1) ? D1D : D1D-1;
+
+         const int ix = idx_ii % nx_i;
+         const int iy = (idx_ii / nx_i) % ny_i;
+         const int iz = (idx_ii / nx_i) / ny_i;
+
+         const real_t (&Bi1)[MQ1][MD1] = (ic == 0) ? r_Bc : r_Bo;
+         const real_t (&Bi2)[MQ1][MD1] = (ic == 1) ? r_Bc : r_Bo;
+         const real_t (&Bi3)[MQ1][MD1] = (ic == 2) ? r_Bc : r_Bo;
+
+         for (int idx_j = 0; idx_j < NDOF; ++idx_j)
+         {
+            const int jc = idx_j / NDOF_C;
+            const int idx_jj = idx_j % NDOF_C;
+
+            const int nx_j = (jc == 0) ? D1D : D1D-1;
+            const int ny_j = (jc == 1) ? D1D : D1D-1;
+
+            const int jx = idx_jj % nx_j;
+            const int jy = (idx_jj / nx_j) % ny_j;
+            const int jz = (idx_jj / nx_j) / ny_j;
+
+            const real_t (&Bj1)[MQ1][MD1] = (jc == 0) ? r_Bc : r_Bo;
+            const real_t (&Bj2)[MQ1][MD1] = (jc == 1) ? r_Bc : r_Bo;
+            const real_t (&Bj3)[MQ1][MD1] = (jc == 2) ? r_Bc : r_Bo;
+
+            real_t val = 0.0;
+            for (int qx = 0; qx < Q1D; ++qx)
+            {
+               for (int qy = 0; qy < Q1D; ++qy)
+               {
+                  for (int qz = 0; qz < Q1D; ++qz)
+                  {
+                     const double coeff = s_D[ic + jc*3][qx][qy][qz];
+                     val += coeff*Bi1[qx][ix]*Bi2[qy][iy]*Bi3[qz][iz]*
+                            Bj1[qx][jx]*Bj2[qy][jy]*Bj3[qz][jz];
+                  }
+               }
+            }
+            if (add)
+            {
+               M(idx_i, idx_j, e) += val;
+            }
+            else
+            {
+               M(idx_i, idx_j, e) = val;
+            }
+         }
+      }
+   });
+}
+
+void VectorFEMassIntegrator::AssembleEA(const FiniteElementSpace &fes,
+                                        Vector &ea_data,
+                                        const bool add)
+{
+   AssemblePA(fes);
+
+   if (trial_fetype != mfem::FiniteElement::DIV ||
+       test_fetype != mfem::FiniteElement::DIV)
+   {
+      MFEM_ABORT("Unsupported kernel.");
+   }
+
+   const Array<real_t> &Bo = mapsO->B;
+   const Array<real_t> &Bc = mapsC->B;
+
+   if (dim == 2)
+   {
+      const int coeff_dim = symmetric ? 3 : 4;
+      auto kernel = EAHdivAssemble2D<0,0>;
+      switch ((dofs1D << 4 ) | quad1D)
+      {
+         case 0x22: kernel = EAHdivAssemble2D<2,2>; break;
+         case 0x33: kernel = EAHdivAssemble2D<3,3>; break;
+         case 0x44: kernel = EAHdivAssemble2D<4,4>; break;
+         case 0x55: kernel = EAHdivAssemble2D<5,5>; break;
+      }
+      return kernel(ne,Bo,Bc,coeff_dim,pa_data,ea_data,add,dofs1D,quad1D);
+   }
+   else if (dim == 3)
+   {
+      const int coeff_dim = symmetric ? 6 : 9;
+      auto kernel = EAHdivAssemble3D<0,0>;
+      switch ((dofs1D << 4 ) | quad1D)
+      {
+         case 0x23: kernel = EAHdivAssemble3D<2,3>; break;
+         case 0x34: kernel = EAHdivAssemble3D<3,4>; break;
+         case 0x45: kernel = EAHdivAssemble3D<4,5>; break;
+         case 0x56: kernel = EAHdivAssemble3D<5,6>; break;
+      }
+      return kernel(ne,Bo,Bc,coeff_dim,pa_data,ea_data,add,dofs1D,quad1D);
+   }
+   MFEM_ABORT("Unknown kernel.");
+}
+
+void DivDivIntegrator::AssembleEA(const FiniteElementSpace &fes,
+                                  Vector &ea_data,
+                                  const bool add)
+{
+   AssemblePA(fes);
+
+   const Array<real_t> &Bo = mapsO->B;
+   const Array<real_t> &Gc = mapsC->G;
+
+   if (dim == 2)
+   {
+      auto kernel = EAHdivAssemble2D<0,0>;
+      switch ((dofs1D << 4 ) | quad1D)
+      {
+         case 0x22: kernel = EAHdivAssemble2D<2,2>; break;
+         case 0x33: kernel = EAHdivAssemble2D<3,3>; break;
+         case 0x44: kernel = EAHdivAssemble2D<4,4>; break;
+         case 0x55: kernel = EAHdivAssemble2D<5,5>; break;
+      }
+      return kernel(ne,Bo,Gc,1,pa_data,ea_data,add,dofs1D,quad1D);
+   }
+   else if (dim == 3)
+   {
+      auto kernel = EAHdivAssemble3D<0,0>;
+      switch ((dofs1D << 4 ) | quad1D)
+      {
+         case 0x23: kernel = EAHdivAssemble3D<2,3>; break;
+         case 0x34: kernel = EAHdivAssemble3D<3,4>; break;
+         case 0x45: kernel = EAHdivAssemble3D<4,5>; break;
+         case 0x56: kernel = EAHdivAssemble3D<5,6>; break;
+      }
+      return kernel(ne,Bo,Gc,1,pa_data,ea_data,add,dofs1D,quad1D);
+   }
+   MFEM_ABORT("Unknown kernel.");
+}
+
+}

fem/integ/bilininteg_mass_ea.cpp

@@ -12,248 +12,67 @@
 #include "../../general/forall.hpp"
 #include "../bilininteg.hpp"
 #include "../gridfunc.hpp"
+#include "bilininteg_mass_kernels.hpp"
 
 namespace mfem
 {
 
-template<int T_D1D = 0, int T_Q1D = 0>
-static void EAMassAssemble1D(const int NE,
-                             const Array<real_t> &basis,
-                             const Vector &padata,
-                             Vector &eadata,
-                             const bool add,
-                             const int d1d = 0,
-                             const int q1d = 0)
+void MassIntegrator::AssembleEA_(Vector &ea_data,
+                                 const bool add)
 {
-   const int D1D = T_D1D ? T_D1D : d1d;
-   const int Q1D = T_Q1D ? T_Q1D : q1d;
-   MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().MAX_D1D, "");
-   MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().MAX_Q1D, "");
-   auto B = Reshape(basis.Read(), Q1D, D1D);
-   auto D = Reshape(padata.Read(), Q1D, NE);
-   auto M = Reshape(add ? eadata.ReadWrite() : eadata.Write(), D1D, D1D, NE);
-   mfem::forall_2D(NE, D1D, D1D, [=] MFEM_HOST_DEVICE (int e)
+   using internal::EAMassAssemble1D;
+   using internal::EAMassAssemble2D;
+   using internal::EAMassAssemble3D;
+
+   const Array<real_t> &B = maps->B;
+   if (dim == 1)
    {
-      const int D1D = T_D1D ? T_D1D : d1d;
-      const int Q1D = T_Q1D ? T_Q1D : q1d;
-      constexpr int MQ1 = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
-      real_t r_Bi[MQ1];
-      real_t r_Bj[MQ1];
-      for (int q = 0; q < Q1D; q++)
+      auto kernel = EAMassAssemble1D<0,0>;
+      switch ((dofs1D << 4 ) | quad1D)
       {
-         r_Bi[q] = B(q,MFEM_THREAD_ID(x));
-         r_Bj[q] = B(q,MFEM_THREAD_ID(y));
+         case 0x22: kernel = EAMassAssemble1D<2,2>; break;
+         case 0x33: kernel = EAMassAssemble1D<3,3>; break;
+         case 0x44: kernel = EAMassAssemble1D<4,4>; break;
+         case 0x55: kernel = EAMassAssemble1D<5,5>; break;
+         case 0x66: kernel = EAMassAssemble1D<6,6>; break;
+         case 0x77: kernel = EAMassAssemble1D<7,7>; break;
+         case 0x88: kernel = EAMassAssemble1D<8,8>; break;
+         case 0x99: kernel = EAMassAssemble1D<9,9>; break;
       }
-      MFEM_FOREACH_THREAD(i1,x,D1D)
-      {
-         MFEM_FOREACH_THREAD(j1,y,D1D)
-         {
-            real_t val = 0.0;
-            for (int k1 = 0; k1 < Q1D; ++k1)
-            {
-               val += r_Bi[k1] * r_Bj[k1] * D(k1, e);
-            }
-            if (add)
-            {
-               M(i1, j1, e) += val;
-            }
-            else
-            {
-               M(i1, j1, e) = val;
-            }
-         }
-      }
-   });
-}
-
-template<int T_D1D = 0, int T_Q1D = 0>
-static void EAMassAssemble2D(const int NE,
-                             const Array<real_t> &basis,
-                             const Vector &padata,
-                             Vector &eadata,
-                             const bool add,
-                             const int d1d = 0,
-                             const int q1d = 0)
-{
-   const int D1D = T_D1D ? T_D1D : d1d;
-   const int Q1D = T_Q1D ? T_Q1D : q1d;
-   MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().MAX_D1D, "");
-   MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().MAX_Q1D, "");
-   auto B = Reshape(basis.Read(), Q1D, D1D);
-   auto D = Reshape(padata.Read(), Q1D, Q1D, NE);
-   auto M = Reshape(add ? eadata.ReadWrite() : eadata.Write(), D1D, D1D, D1D, D1D,
-                    NE);
-   mfem::forall_2D(NE, D1D, D1D, [=] MFEM_HOST_DEVICE (int e)
+      return kernel(ne,B,pa_data,ea_data,add,dofs1D,quad1D);
+   }
+   else if (dim == 2)
    {
-      const int D1D = T_D1D ? T_D1D : d1d;
-      const int Q1D = T_Q1D ? T_Q1D : q1d;
-      constexpr int MD1 = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
-      constexpr int MQ1 = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
-      real_t r_B[MQ1][MD1];
-      for (int d = 0; d < D1D; d++)
+      auto kernel = EAMassAssemble2D<0,0>;
+      switch ((dofs1D << 4 ) | quad1D)
       {
-         for (int q = 0; q < Q1D; q++)
-         {
-            r_B[q][d] = B(q,d);
-         }
+         case 0x22: kernel = EAMassAssemble2D<2,2>; break;
+         case 0x33: kernel = EAMassAssemble2D<3,3>; break;
+         case 0x44: kernel = EAMassAssemble2D<4,4>; break;
+         case 0x55: kernel = EAMassAssemble2D<5,5>; break;
+         case 0x66: kernel = EAMassAssemble2D<6,6>; break;
+         case 0x77: kernel = EAMassAssemble2D<7,7>; break;
+         case 0x88: kernel = EAMassAssemble2D<8,8>; break;
+         case 0x99: kernel = EAMassAssemble2D<9,9>; break;
       }
-      MFEM_SHARED real_t s_D[MQ1][MQ1];
-      MFEM_FOREACH_THREAD(k1,x,Q1D)
-      {
-         MFEM_FOREACH_THREAD(k2,y,Q1D)
-         {
-            s_D[k1][k2] = D(k1,k2,e);
-         }
-      }
-      MFEM_SYNC_THREAD;
-      MFEM_FOREACH_THREAD(i1,x,D1D)
-      {
-         MFEM_FOREACH_THREAD(i2,y,D1D)
-         {
-            for (int j1 = 0; j1 < D1D; ++j1)
-            {
-               for (int j2 = 0; j2 < D1D; ++j2)
-               {
-                  real_t val = 0.0;
-                  for (int k1 = 0; k1 < Q1D; ++k1)
-                  {
-                     for (int k2 = 0; k2 < Q1D; ++k2)
-                     {
-                        val += r_B[k1][i1] * r_B[k1][j1]
-                               * r_B[k2][i2] * r_B[k2][j2]
-                               * s_D[k1][k2];
-                     }
-                  }
-                  if (add)
-                  {
-                     M(i1, i2, j1, j2, e) += val;
-                  }
-                  else
-                  {
-                     M(i1, i2, j1, j2, e) = val;
-                  }
-               }
-            }
-         }
-      }
-   });
-}
-
-template<int T_D1D = 0, int T_Q1D = 0>
-static void EAMassAssemble3D(const int NE,
-                             const Array<real_t> &basis,
-                             const Vector &padata,
-                             Vector &eadata,
-                             const bool add,
-                             const int d1d = 0,
-                             const int q1d = 0)
-{
-   const int D1D = T_D1D ? T_D1D : d1d;
-   const int Q1D = T_Q1D ? T_Q1D : q1d;
-   MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().MAX_D1D, "");
-   MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().MAX_Q1D, "");
-   auto B = Reshape(basis.Read(), Q1D, D1D);
-   auto D = Reshape(padata.Read(), Q1D, Q1D, Q1D, NE);
-   auto M = Reshape(add ? eadata.ReadWrite() : eadata.Write(), D1D, D1D, D1D, D1D,
-                    D1D, D1D, NE);
-   mfem::forall_3D(NE, D1D, D1D, D1D, [=] MFEM_HOST_DEVICE (int e)
+      return kernel(ne,B,pa_data,ea_data,add,dofs1D,quad1D);
+   }
+   else if (dim == 3)
    {
-      const int D1D = T_D1D ? T_D1D : d1d;
-      const int Q1D = T_Q1D ? T_Q1D : q1d;
-      constexpr int MD1 = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
-      constexpr int MQ1 = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
-      constexpr int DQ = T_D1D * T_Q1D;
-
-      // For quadratic and lower it's better to use registers but for higher-order you start to
-      // spill and it's better to use shared memory
-      constexpr bool USE_REG = DQ != 0 && DQ <= 12;
-      constexpr int MD1r = USE_REG ? MD1 : 1;
-      constexpr int MQ1r = USE_REG ? MQ1 : 1;
-      constexpr int MD1s = USE_REG ? 1 : MD1;
-      constexpr int MQ1s = USE_REG ? 1 : MQ1;
-
-      MFEM_SHARED real_t s_B[MQ1s][MD1s];
-      real_t r_B[MQ1r][MD1r];
-      real_t (*l_B)[MD1] = nullptr;
-      if (USE_REG)
+      auto kernel = EAMassAssemble3D<0,0>;
+      switch ((dofs1D << 4 ) | quad1D)
       {
-         for (int d = 0; d < D1D; d++)
-         {
-            for (int q = 0; q < Q1D; q++)
-            {
-               r_B[q][d] = B(q,d);
-            }
-         }
-         l_B = (real_t (*)[MD1])r_B;
+         case 0x23: kernel = EAMassAssemble3D<2,3>; break;
+         case 0x34: kernel = EAMassAssemble3D<3,4>; break;
+         case 0x45: kernel = EAMassAssemble3D<4,5>; break;
+         case 0x56: kernel = EAMassAssemble3D<5,6>; break;
+         case 0x67: kernel = EAMassAssemble3D<6,7>; break;
+         case 0x78: kernel = EAMassAssemble3D<7,8>; break;
+         case 0x89: kernel = EAMassAssemble3D<8,9>; break;
       }
-      else
-      {
-         if (MFEM_THREAD_ID(z) == 0)
-         {
-            MFEM_FOREACH_THREAD(d,x,D1D)
-            {
-               MFEM_FOREACH_THREAD(q,y,Q1D)
-               {
-                  s_B[q][d] = B(q,d);
-               }
-            }
-         }
-         l_B = (real_t (*)[MD1])s_B;
-      }
-
-      MFEM_SHARED real_t s_D[MQ1][MQ1][MQ1];
-      MFEM_FOREACH_THREAD(k1,x,Q1D)
-      {
-         MFEM_FOREACH_THREAD(k2,y,Q1D)
-         {
-            MFEM_FOREACH_THREAD(k3,z,Q1D)
-            {
-               s_D[k1][k2][k3] = D(k1,k2,k3,e);
-            }
-         }
-      }
-      MFEM_SYNC_THREAD;
-      MFEM_FOREACH_THREAD(i1,x,D1D)
-      {
-         MFEM_FOREACH_THREAD(i2,y,D1D)
-         {
-            MFEM_FOREACH_THREAD(i3,z,D1D)
-            {
-               for (int j1 = 0; j1 < D1D; ++j1)
-               {
-                  for (int j2 = 0; j2 < D1D; ++j2)
-                  {
-                     for (int j3 = 0; j3 < D1D; ++j3)
-                     {
-                        real_t val = 0.0;
-                        for (int k1 = 0; k1 < Q1D; ++k1)
-                        {
-                           for (int k2 = 0; k2 < Q1D; ++k2)
-                           {
-                              for (int k3 = 0; k3 < Q1D; ++k3)
-                              {
-                                 val += l_B[k1][i1] * l_B[k1][j1]
-                                        * l_B[k2][i2] * l_B[k2][j2]
-                                        * l_B[k3][i3] * l_B[k3][j3]
-                                        * s_D[k1][k2][k3];
-                              }
-                           }
-                        }
-                        if (add)
-                        {
-                           M(i1, i2, i3, j1, j2, j3, e) += val;
-                        }
-                        else
-                        {
-                           M(i1, i2, i3, j1, j2, j3, e) = val;
-                        }
-                     }
-                  }
-               }
-            }
-         }
-      }
-   });
+      return kernel(ne,B,pa_data,ea_data,add,dofs1D,quad1D);
+   }
+   MFEM_ABORT("Unknown kernel.");
 }
 
 void MassIntegrator::AssembleEA(const FiniteElementSpace &fes,
@@ -261,56 +80,15 @@ void MassIntegrator::AssembleEA(const FiniteElementSpace &fes,
                                 const bool add)
 {
    AssemblePA(fes);
-   ne = fes.GetMesh()->GetNE();
-   const Array<real_t> &B = maps->B;
-   if (dim == 1)
-   {
-      switch ((dofs1D << 4 ) | quad1D)
-      {
-         case 0x22: return EAMassAssemble1D<2,2>(ne,B,pa_data,ea_data,add);
-         case 0x33: return EAMassAssemble1D<3,3>(ne,B,pa_data,ea_data,add);
-         case 0x44: return EAMassAssemble1D<4,4>(ne,B,pa_data,ea_data,add);
-         case 0x55: return EAMassAssemble1D<5,5>(ne,B,pa_data,ea_data,add);
-         case 0x66: return EAMassAssemble1D<6,6>(ne,B,pa_data,ea_data,add);
-         case 0x77: return EAMassAssemble1D<7,7>(ne,B,pa_data,ea_data,add);
-         case 0x88: return EAMassAssemble1D<8,8>(ne,B,pa_data,ea_data,add);
-         case 0x99: return EAMassAssemble1D<9,9>(ne,B,pa_data,ea_data,add);
-         default:   return EAMassAssemble1D(ne,B,pa_data,ea_data,add,
-                                               dofs1D,quad1D);
-      }
-   }
-   else if (dim == 2)
-   {
-      switch ((dofs1D << 4 ) | quad1D)
-      {
-         case 0x22: return EAMassAssemble2D<2,2>(ne,B,pa_data,ea_data,add);
-         case 0x33: return EAMassAssemble2D<3,3>(ne,B,pa_data,ea_data,add);
-         case 0x44: return EAMassAssemble2D<4,4>(ne,B,pa_data,ea_data,add);
-         case 0x55: return EAMassAssemble2D<5,5>(ne,B,pa_data,ea_data,add);
-         case 0x66: return EAMassAssemble2D<6,6>(ne,B,pa_data,ea_data,add);
-         case 0x77: return EAMassAssemble2D<7,7>(ne,B,pa_data,ea_data,add);
-         case 0x88: return EAMassAssemble2D<8,8>(ne,B,pa_data,ea_data,add);
-         case 0x99: return EAMassAssemble2D<9,9>(ne,B,pa_data,ea_data,add);
-         default:   return EAMassAssemble2D(ne,B,pa_data,ea_data,add,
-                                               dofs1D,quad1D);
-      }
-   }
-   else if (dim == 3)
-   {
-      switch ((dofs1D << 4 ) | quad1D)
-      {
-         case 0x23: return EAMassAssemble3D<2,3>(ne,B,pa_data,ea_data,add);
-         case 0x34: return EAMassAssemble3D<3,4>(ne,B,pa_data,ea_data,add);
-         case 0x45: return EAMassAssemble3D<4,5>(ne,B,pa_data,ea_data,add);
-         case 0x56: return EAMassAssemble3D<5,6>(ne,B,pa_data,ea_data,add);
-         case 0x67: return EAMassAssemble3D<6,7>(ne,B,pa_data,ea_data,add);
-         case 0x78: return EAMassAssemble3D<7,8>(ne,B,pa_data,ea_data,add);
-         case 0x89: return EAMassAssemble3D<8,9>(ne,B,pa_data,ea_data,add);
-         default:   return EAMassAssemble3D(ne,B,pa_data,ea_data,add,
-                                               dofs1D,quad1D);
-      }
-   }
-   MFEM_ABORT("Unknown kernel.");
+   if (ne > 0) { AssembleEA_(ea_data, add); }
+}
+
+void MassIntegrator::AssembleEABoundary(const FiniteElementSpace &fes,
+                                        Vector &ea_data,
+                                        const bool add)
+{
+   AssemblePABoundary(fes);
+   if (ne > 0) { AssembleEA_(ea_data, add); }
 }
 
 }

fem/integ/bilininteg_mass_kernels.hpp

@@ -1132,6 +1132,244 @@ inline void SmemPAMassApply3D(const int NE,
    });
 }
 
+template<int T_D1D = 0, int T_Q1D = 0>
+inline void EAMassAssemble1D(const int NE,
+                             const Array<real_t> &basis,
+                             const Vector &padata,
+                             Vector &eadata,
+                             const bool add,
+                             const int d1d = 0,
+                             const int q1d = 0)
+{
+   const int D1D = T_D1D ? T_D1D : d1d;
+   const int Q1D = T_Q1D ? T_Q1D : q1d;
+   MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().MAX_D1D, "");
+   MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().MAX_Q1D, "");
+   auto B = Reshape(basis.Read(), Q1D, D1D);
+   auto D = Reshape(padata.Read(), Q1D, NE);
+   auto M = Reshape(add ? eadata.ReadWrite() : eadata.Write(), D1D, D1D, NE);
+   mfem::forall_2D(NE, D1D, D1D, [=] MFEM_HOST_DEVICE (int e)
+   {
+      const int D1D = T_D1D ? T_D1D : d1d;
+      const int Q1D = T_Q1D ? T_Q1D : q1d;
+      constexpr int MQ1 = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
+      MFEM_FOREACH_THREAD(i1,x,D1D)
+      {
+         real_t r_Bi[MQ1];
+         for (int q = 0; q < Q1D; q++) { r_Bi[q] = B(q,i1); }
+         MFEM_FOREACH_THREAD(j1,y,D1D)
+         {
+            real_t r_Bj[MQ1];
+            for (int q = 0; q < Q1D; q++) { r_Bj[q] = B(q,j1); }
+
+            real_t val = 0.0;
+            for (int k1 = 0; k1 < Q1D; ++k1)
+            {
+               val += r_Bi[k1] * r_Bj[k1] * D(k1, e);
+            }
+            if (add)
+            {
+               M(i1, j1, e) += val;
+            }
+            else
+            {
+               M(i1, j1, e) = val;
+            }
+         }
+      }
+   });
+}
+
+template<int T_D1D = 0, int T_Q1D = 0>
+inline void EAMassAssemble2D(const int NE,
+                             const Array<real_t> &basis,
+                             const Vector &padata,
+                             Vector &eadata,
+                             const bool add,
+                             const int d1d = 0,
+                             const int q1d = 0)
+{
+   const int D1D = T_D1D ? T_D1D : d1d;
+   const int Q1D = T_Q1D ? T_Q1D : q1d;
+   MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().MAX_D1D, "");
+   MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().MAX_Q1D, "");
+   auto B = Reshape(basis.Read(), Q1D, D1D);
+   auto D = Reshape(padata.Read(), Q1D, Q1D, NE);
+   auto M = Reshape(add ? eadata.ReadWrite() : eadata.Write(), D1D, D1D, D1D, D1D,
+                    NE);
+   mfem::forall_2D(NE, D1D, D1D, [=] MFEM_HOST_DEVICE (int e)
+   {
+      const int D1D = T_D1D ? T_D1D : d1d;
+      const int Q1D = T_Q1D ? T_Q1D : q1d;
+      constexpr int MD1 = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
+      constexpr int MQ1 = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
+      real_t r_B[MQ1][MD1];
+      for (int d = 0; d < D1D; d++)
+      {
+         for (int q = 0; q < Q1D; q++)
+         {
+            r_B[q][d] = B(q,d);
+         }
+      }
+      MFEM_SHARED real_t s_D[MQ1][MQ1];
+      MFEM_FOREACH_THREAD(k1,x,Q1D)
+      {
+         MFEM_FOREACH_THREAD(k2,y,Q1D)
+         {
+            s_D[k1][k2] = D(k1,k2,e);
+         }
+      }
+      MFEM_SYNC_THREAD;
+      MFEM_FOREACH_THREAD(i1,x,D1D)
+      {
+         MFEM_FOREACH_THREAD(i2,y,D1D)
+         {
+            for (int j1 = 0; j1 < D1D; ++j1)
+            {
+               for (int j2 = 0; j2 < D1D; ++j2)
+               {
+                  real_t val = 0.0;
+                  for (int k1 = 0; k1 < Q1D; ++k1)
+                  {
+                     for (int k2 = 0; k2 < Q1D; ++k2)
+                     {
+                        val += r_B[k1][i1] * r_B[k1][j1]
+                               * r_B[k2][i2] * r_B[k2][j2]
+                               * s_D[k1][k2];
+                     }
+                  }
+                  if (add)
+                  {
+                     M(i1, i2, j1, j2, e) += val;
+                  }
+                  else
+                  {
+                     M(i1, i2, j1, j2, e) = val;
+                  }
+               }
+            }
+         }
+      }
+   });
+}
+
+template<int T_D1D = 0, int T_Q1D = 0>
+inline void EAMassAssemble3D(const int NE,
+                             const Array<real_t> &basis,
+                             const Vector &padata,
+                             Vector &eadata,
+                             const bool add,
+                             const int d1d = 0,
+                             const int q1d = 0)
+{
+   const int D1D = T_D1D ? T_D1D : d1d;
+   const int Q1D = T_Q1D ? T_Q1D : q1d;
+   MFEM_VERIFY(D1D <= DeviceDofQuadLimits::Get().MAX_D1D, "");
+   MFEM_VERIFY(Q1D <= DeviceDofQuadLimits::Get().MAX_Q1D, "");
+   auto B = Reshape(basis.Read(), Q1D, D1D);
+   auto D = Reshape(padata.Read(), Q1D, Q1D, Q1D, NE);
+   auto M = Reshape(add ? eadata.ReadWrite() : eadata.Write(), D1D, D1D, D1D, D1D,
+                    D1D, D1D, NE);
+   mfem::forall_3D(NE, D1D, D1D, D1D, [=] MFEM_HOST_DEVICE (int e)
+   {
+      const int D1D = T_D1D ? T_D1D : d1d;
+      const int Q1D = T_Q1D ? T_Q1D : q1d;
+      constexpr int MD1 = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
+      constexpr int MQ1 = T_Q1D ? T_Q1D : DofQuadLimits::MAX_Q1D;
+      constexpr int DQ = T_D1D * T_Q1D;
+
+      // For quadratic and lower it's better to use registers but for higher-order you start to
+      // spill and it's better to use shared memory
+      constexpr bool USE_REG = DQ != 0 && DQ <= 12;
+      constexpr int MD1r = USE_REG ? MD1 : 1;
+      constexpr int MQ1r = USE_REG ? MQ1 : 1;
+      constexpr int MD1s = USE_REG ? 1 : MD1;
+      constexpr int MQ1s = USE_REG ? 1 : MQ1;
+
+      MFEM_SHARED real_t s_B[MQ1s][MD1s];
+      real_t r_B[MQ1r][MD1r];
+      real_t (*l_B)[MD1] = nullptr;
+      if (USE_REG)
+      {
+         for (int d = 0; d < D1D; d++)
+         {
+            for (int q = 0; q < Q1D; q++)
+            {
+               r_B[q][d] = B(q,d);
+            }
+         }
+         l_B = (real_t (*)[MD1])r_B;
+      }
+      else
+      {
+         if (MFEM_THREAD_ID(z) == 0)
+         {
+            MFEM_FOREACH_THREAD(d,x,D1D)
+            {
+               MFEM_FOREACH_THREAD(q,y,Q1D)
+               {
+                  s_B[q][d] = B(q,d);
+               }
+            }
+         }
+         l_B = (real_t (*)[MD1])s_B;
+      }
+
+      MFEM_SHARED real_t s_D[MQ1][MQ1][MQ1];
+      MFEM_FOREACH_THREAD(k1,x,Q1D)
+      {
+         MFEM_FOREACH_THREAD(k2,y,Q1D)
+         {
+            MFEM_FOREACH_THREAD(k3,z,Q1D)
+            {
+               s_D[k1][k2][k3] = D(k1,k2,k3,e);
+            }
+         }
+      }
+      MFEM_SYNC_THREAD;
+      MFEM_FOREACH_THREAD(i1,x,D1D)
+      {
+         MFEM_FOREACH_THREAD(i2,y,D1D)
+         {
+            MFEM_FOREACH_THREAD(i3,z,D1D)
+            {
+               for (int j1 = 0; j1 < D1D; ++j1)
+               {
+                  for (int j2 = 0; j2 < D1D; ++j2)
+                  {
+                     for (int j3 = 0; j3 < D1D; ++j3)
+                     {
+                        real_t val = 0.0;
+                        for (int k1 = 0; k1 < Q1D; ++k1)
+                        {
+                           for (int k2 = 0; k2 < Q1D; ++k2)
+                           {
+                              for (int k3 = 0; k3 < Q1D; ++k3)
+                              {
+                                 val += l_B[k1][i1] * l_B[k1][j1]
+                                        * l_B[k2][i2] * l_B[k2][j2]
+                                        * l_B[k3][i3] * l_B[k3][j3]
+                                        * s_D[k1][k2][k3];
+                              }
+                           }
+                        }
+                        if (add)
+                        {
+                           M(i1, i2, i3, j1, j2, j3, e) += val;
+                        }
+                        else
+                        {
+                           M(i1, i2, i3, j1, j2, j3, e) = val;
+                        }
+                     }
+                  }
+               }
+            }
+         }
+      }
+   });
+}
+
 } // namespace internal
 
 namespace

fem/integ/bilininteg_mass_pa.cpp

@@ -29,7 +29,8 @@ void MassIntegrator::AssemblePA(const FiniteElementSpace &fes)
    // Assuming the same element type
    fespace = &fes;
    Mesh *mesh = fes.GetMesh();
-   if (mesh->GetNE() == 0) { return; }
+   ne = fes.GetMesh()->GetNE();
+   if (ne == 0) { return; }
    const FiniteElement &el = *fes.GetFE(0);
    ElementTransformation *T0 = mesh->GetElementTransformation(0);
    const IntegrationRule *ir = IntRule ? IntRule : &GetRule(el, el, *T0);
@@ -50,7 +51,6 @@ void MassIntegrator::AssemblePA(const FiniteElementSpace &fes)
    }
    int map_type = el.GetMapType();
    dim = mesh->Dimension();
-   ne = fes.GetMesh()->GetNE();
    nq = ir->GetNPoints();
    geom = mesh->GetGeometricFactors(*ir, GeometricFactors::DETERMINANTS, mt);
    maps = &el.GetDofToQuad(*ir, DofToQuad::TENSOR);
@@ -123,14 +123,14 @@ void MassIntegrator::AssemblePABoundary(const FiniteElementSpace &fes)
    // Assuming the same element type
    fespace = &fes;
    Mesh *mesh = fes.GetMesh();
-   if (mesh->GetNBE() == 0) { return; }
+   ne = mesh->GetNFbyType(FaceType::Boundary);
+   if (ne == 0) { return; }
    const FiniteElement &el = *fes.GetBE(0);
    ElementTransformation *T0 = mesh->GetBdrElementTransformation(0);
    const IntegrationRule *ir = IntRule ? IntRule : &GetRule(el, el, *T0);
 
    int map_type = el.GetMapType();
    dim = el.GetDim(); // Dimension of the boundary element, *not* the mesh
-   ne = fes.GetMesh()->GetNFbyType(FaceType::Boundary);
    nq = ir->GetNPoints();
    face_geom = mesh->GetFaceGeometricFactors(*ir, GeometricFactors::DETERMINANTS,
                                              FaceType::Boundary, mt);

fem/integ/bilininteg_trace_jump_ea.cpp

@@ -0,0 +1,158 @@
+// Copyright (c) 2010-2024, 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 "../../general/forall.hpp"
+#include "../fe/face_map_utils.hpp"
+#include "../bilininteg.hpp"
+#include "bilininteg_mass_kernels.hpp"
+
+namespace mfem
+{
+
+void NormalTraceJumpIntegrator::AssembleEAInteriorFaces(
+   const FiniteElementSpace &trial_fes,
+   const FiniteElementSpace &test_fes,
+   Vector &emat,
+   const bool add)
+{
+   Mesh &mesh = *trial_fes.GetMesh();
+   const int dim = mesh.Dimension();
+   const FaceType ftype = FaceType::Interior;
+   const int nf = mesh.GetNFbyType(ftype);
+
+   const Geometry::Type geom = mesh.GetFaceGeometry(0);
+   const int trial_order = trial_fes.GetMaxElementOrder();
+   const int test_order = test_fes.GetMaxElementOrder();
+   const int qorder = test_order + trial_order - 1;
+   const IntegrationRule &ir = IntRule ? *IntRule : IntRules.Get(geom, qorder);
+   const int nquad = ir.Size();
+
+   Vector pa_data(nquad * nf);
+   {
+      const auto d_w = ir.GetWeights().Read();
+      auto d_pa_data = Reshape(pa_data.Write(), nquad, nf);
+      mfem::forall(nquad * nf, [=] MFEM_HOST_DEVICE (int idx)
+      {
+         const int q = idx % nquad;
+         const int f = idx / nquad;
+         d_pa_data(q, f) = d_w[q];
+      });
+   }
+
+   const FiniteElement &trial_face_el = *trial_fes.GetFaceElement(0);
+   const auto maps = &trial_face_el.GetDofToQuad(ir, DofToQuad::TENSOR);
+   const int ndof_face = trial_face_el.GetDof();
+
+   const Array<real_t> &B = maps->B;
+   const int d1d = maps->ndof;
+   const int q1d = maps->nqpt;
+
+   Vector mass_emat(ndof_face*ndof_face*nf);
+
+   // Note: dim is the element dimension, and we integrate over the faces (one
+   // dimension less)
+   if (dim == 2)
+   {
+      internal::EAMassAssemble1D(nf, B, pa_data, mass_emat, false, d1d, q1d);
+   }
+   else if (dim == 3)
+   {
+      internal::EAMassAssemble2D(nf, B, pa_data, mass_emat, false, d1d, q1d);
+   }
+   else
+   {
+      MFEM_ABORT("Unknown kernel.");
+   }
+
+   const FiniteElement &test_el = *test_fes.GetFE(0);
+   const int n_faces_per_el = 2*dim; // assuming tensor product
+   // Get all the local face maps (mapping from lexicographic face index to
+   // lexicographic volume index, depending on the local face index).
+   Array<int> face_maps(ndof_face * n_faces_per_el);
+   for (int lf_i = 0; lf_i < n_faces_per_el; ++lf_i)
+   {
+      Array<int> face_map(ndof_face);
+      test_el.GetFaceMap(lf_i, face_map);
+      for (int i = 0; i < ndof_face; ++i)
+      {
+         face_maps[i + lf_i*ndof_face] = face_map[i];
+      }
+   }
+
+   Array<int> face_info(nf * 4);
+   {
+      int fidx = 0;
+      for (int f = 0; f < mesh.GetNumFaces(); ++f)
+      {
+         Mesh::FaceInformation finfo = mesh.GetFaceInformation(f);
+         if (!finfo.IsInterior()) { continue; }
+         face_info[0 + fidx*4] = finfo.element[0].local_face_id;
+         face_info[1 + fidx*4] = finfo.element[0].orientation;
+         face_info[2 + fidx*4] = finfo.element[1].local_face_id;
+         face_info[3 + fidx*4] = finfo.element[1].orientation;
+         fidx++;
+      }
+   }
+
+   const int ndof_vol = test_el.GetDof();
+   const auto d_face_maps = Reshape(face_maps.Read(), ndof_face, n_faces_per_el);
+   const auto d_face_info = Reshape(face_info.Read(), 2, 2, nf);
+
+   real_t *d_emat;
+   if (add)
+   {
+      d_emat = emat.ReadWrite();
+   }
+   if (!add)
+   {
+      d_emat = emat.Write();
+      mfem::forall(emat.Size(), [=] MFEM_HOST_DEVICE (int i) { d_emat[i] = 0.0; });
+   }
+
+   const auto face_mats = Reshape(mass_emat.Read(), ndof_face, ndof_face, nf);
+   auto el_mats = Reshape(d_emat, ndof_vol, ndof_face, 2, nf);
+
+   auto permute_face = [=] MFEM_HOST_DEVICE(int local_face_id, int orient,
+                                            int size1d, int index)
+   {
+      if (dim == 2)
+      {
+         return internal::PermuteFace2D(local_face_id, orient, size1d, index);
+      }
+      else // dim == 3
+      {
+         return internal::PermuteFace3D(local_face_id, orient, size1d, index);
+      }
+   };
+
+   mfem::forall_3D(nf, ndof_face, ndof_face, 2, [=] MFEM_HOST_DEVICE (int f)
+   {
+      MFEM_FOREACH_THREAD(el_i, z, 2)
+      {
+         const int lf_i = d_face_info(0, el_i, f);
+         const int orient = d_face_info(1, el_i, f);
+         // Loop over face indices in "native ordering"
+         MFEM_FOREACH_THREAD(i_lex, x, ndof_face)
+         {
+            // Convert to lexicographic relative to the face itself
+            const int i_face = permute_face(lf_i, orient, d1d, i_lex);
+            // Convert from lexicographic face DOF to volume DOF
+            const int i = d_face_maps(i_lex, lf_i);
+            MFEM_FOREACH_THREAD(j, y, ndof_face)
+            {
+               el_mats(i, j, el_i, f) += face_mats(i_face, j, f);
+            }
+         }
+      }
+   });
+}
+
+}

fem/pbilinearform.cpp

@@ -438,9 +438,31 @@ void ParBilinearForm::FormLinearSystem(
    const Array<int> &ess_tdof_list, Vector &x, Vector &b,
    OperatorHandle &A, Vector &X, Vector &B, int copy_interior)
 {
+   const Operator &P = *pfes->GetProlongationMatrix();
+   const SparseMatrix &R = *pfes->GetRestrictionMatrix();
+
    if (ext)
    {
-      ext->FormLinearSystem(ess_tdof_list, x, b, A, X, B, copy_interior);
+      if (hybridization)
+      {
+         HypreParVector true_X(pfes), true_B(pfes);
+         P.MultTranspose(b, true_B);
+         R.Mult(x, true_X);
+
+         FormSystemMatrix(ess_tdof_list, A);
+         ConstrainedOperator *A_constrained;
+         Operator::FormConstrainedSystemOperator(ess_tdof_list, A_constrained);
+         A_constrained->EliminateRHS(true_X, true_B);
+         delete A_constrained;
+         R.MultTranspose(true_B, b);
+         hybridization->ReduceRHS(true_B, B);
+         X.SetSize(B.Size());
+         X = 0.0;
+      }
+      else
+      {
+         ext->FormLinearSystem(ess_tdof_list, x, b, A, X, B, copy_interior);
+      }
       return;
    }
 
@@ -448,9 +470,6 @@ void ParBilinearForm::FormLinearSystem(
    // eliminated part of the matrix.
    FormSystemMatrix(ess_tdof_list, A);
 
-   const Operator &P = *pfes->GetProlongationMatrix();
-   const SparseMatrix &R = *pfes->GetRestrictionMatrix();
-
    // Transform the system and perform the elimination in B, based on the
    // essential BC values from x. Restrict the BC part of x in X, and set the
    // non-BC part to zero. Since there is no good initial guess for the Lagrange
@@ -495,7 +514,16 @@ void ParBilinearForm::FormSystemMatrix(const Array<int> &ess_tdof_list,
 {
    if (ext)
    {
-      ext->FormSystemMatrix(ess_tdof_list, A);
+      if (hybridization)
+      {
+         const int remove_zeros = 0;
+         Finalize(remove_zeros);
+         hybridization->GetParallelMatrix(A);
+      }
+      else
+      {
+         ext->FormSystemMatrix(ess_tdof_list, A);
+      }
       return;
    }
 
@@ -541,7 +569,7 @@ void ParBilinearForm::FormSystemMatrix(const Array<int> &ess_tdof_list,
 void ParBilinearForm::RecoverFEMSolution(
    const Vector &X, const Vector &b, Vector &x)
 {
-   if (ext)
+   if (ext && !hybridization)
    {
       ext->RecoverFEMSolution(X, b, x);
       return;

fem/pfespace.cpp

@@ -564,6 +564,10 @@ const FaceRestriction *ParFiniteElementSpace::GetFaceRestriction(
             res = new ParNCL2FaceRestriction(*this, f_ordering, type, m);
          }
       }
+      else if (dynamic_cast<const DG_Interface_FECollection*>(fec))
+      {
+         res = new L2InterfaceFaceRestriction(*this, f_ordering, type);
+      }
       else
       {
          if (Conforming())

fem/restriction.cpp

@@ -56,7 +56,7 @@ ElementRestriction::ElementRestriction(const FiniteElementSpace &f,
          dynamic_cast<const TensorBasisElement*>(fe);
       const Array<int> &fe_dof_map = el->GetDofMap();
       MFEM_VERIFY(fe_dof_map.Size() > 0, "invalid dof map");
-      dof_map = fe_dof_map.GetData();
+      dof_map = fe_dof_map.HostRead();
    }
    const Table& e2dTable = fes.GetElementToDofTable();
    const int* element_map = e2dTable.GetJ();
@@ -2284,6 +2284,89 @@ void NCL2FaceRestriction::ComputeGatherIndices()
    gather_offsets[0] = 0;
 }
 
+L2InterfaceFaceRestriction::L2InterfaceFaceRestriction(
+   const FiniteElementSpace& fes_,
+   const ElementDofOrdering ordering_,
+   const FaceType type_)
+   : fes(fes_),
+     ordering(ordering_),
+     type(type_),
+     nfaces(fes.GetNFbyType(type)),
+     vdim(fes.GetVDim()),
+     byvdim(fes.GetOrdering() == Ordering::byVDIM),
+     face_dofs(nfaces > 0 ? fes.GetFaceElement(0)->GetDof() : 0),
+     nfdofs(face_dofs*nfaces),
+     ndofs(fes.GetNDofs())
+{
+   height = nfdofs;
+   width = ndofs;
+
+   const Table &face2dof = fes.GetFaceToDofTable();
+
+   const Mesh &mesh = *fes.GetMesh();
+   int face_idx = 0;
+   gather_map.SetSize(nfdofs);
+   for (int f = 0; f < mesh.GetNumFaces(); ++f)
+   {
+      Mesh::FaceInformation face = mesh.GetFaceInformation(f);
+      if (!face.IsOfFaceType(type)) { continue; }
+      for (int i = 0; i < face_dofs; ++i)
+      {
+         gather_map[i + face_idx*face_dofs] = face2dof.GetJ()[i + f*face_dofs];
+      }
+      ++face_idx;
+   }
+}
+
+void L2InterfaceFaceRestriction::Mult(const Vector &x, Vector &y) const
+{
+   const int nd = face_dofs;
+   const int nf = nfaces;
+   const int vd = vdim;
+   const bool t = byvdim;
+   const int *map = gather_map.Read();
+
+   const auto d_x = Reshape(x.Read(), t?vd:ndofs, t?ndofs:vd);
+   auto d_y = Reshape(y.Write(), nd, vd, nf);
+
+   mfem::forall(nd*nf, [=] MFEM_HOST_DEVICE (int i)
+   {
+      const int j = map[i];
+      for (int c = 0; c < vd; ++c)
+      {
+         d_y(i % nd, c, i / nd) = d_x(t?c:j, t?j:c);
+      }
+   });
+}
+
+void L2InterfaceFaceRestriction::AddMultTranspose(
+   const Vector &x, Vector &y, const real_t a) const
+{
+   const int nd = face_dofs;
+   const int nf = nfaces;
+   const int vd = vdim;
+   const bool t = byvdim;
+   const int *map = gather_map.Read();
+
+   const auto d_x = Reshape(x.Read(), nd, vd, nf);
+   auto d_y = Reshape(y.Write(), t?vd:ndofs, t?ndofs:vd);
+
+   mfem::forall(ndofs, [=] MFEM_HOST_DEVICE (int i) { d_y[i] = 0.0; });
+   mfem::forall(nd*nf, [=] MFEM_HOST_DEVICE (int i)
+   {
+      const int j = map[i];
+      for (int c = 0; c < vd; ++c)
+      {
+         d_y(t?c:j, t?j:c) = d_x(i % nd, c, i / nd);
+      }
+   });
+}
+
+const Array<int> &L2InterfaceFaceRestriction::GatherMap() const
+{
+   return gather_map;
+}
+
 Vector GetLVectorFaceNbrData(
    const FiniteElementSpace &fes, const Vector &x, FaceType ftype)
 {

fem/restriction.hpp

@@ -262,6 +262,12 @@ public:
    {
       MFEM_ABORT("Not implemented for this restriction operator.");
    }
+
+   /// @brief Low-level access to the underlying gather map.
+   virtual const Array<int> &GatherMap() const
+   {
+      MFEM_ABORT("Not implemented for this restriction operator.");
+   }
 };
 
 /// @brief Operator that extracts face degrees of freedom for H1, ND, or RT
@@ -1066,6 +1072,62 @@ public:
    void DoubleValuedNonconformingTransposeInterpolationInPlace(Vector& x) const;
 };
 
+/// Operator that extracts face degrees of freedom for L2 interface spaces.
+/** Objects of this type are typically created and owned by FiniteElementSpace
+    objects, see FiniteElementSpace::GetFaceRestriction(). */
+class L2InterfaceFaceRestriction : public FaceRestriction
+{
+protected:
+   const FiniteElementSpace &fes; ///< The finite element space
+   const ElementDofOrdering ordering; ///< Requested ordering
+   const FaceType type; ///< Face type (interior or boundary)
+   const int nfaces; ///< Number of faces of the requested type
+   const int vdim; ///< vdim of the space
+   const bool byvdim; ///< DOF ordering (by nodes or by vdim)
+   const int face_dofs; ///< Number of dofs on each face
+   const int nfdofs; ///< Total number of dofs on the faces (E-vector size)
+   const int ndofs; ///< Number of dofs in the space (L-vector size)
+   Array<int> gather_map; ///< Gather map
+
+public:
+   /** @brief Constructs an L2InterfaceFaceRestriction.
+
+       @param[in] fes_       The FiniteElementSpace on which this operates
+       @param[in] ordering_  Request a specific face dof ordering
+       @param[in] type_      Request internal or boundary faces dofs */
+   L2InterfaceFaceRestriction(const FiniteElementSpace& fes_,
+                              const ElementDofOrdering ordering_,
+                              const FaceType type_);
+
+   /** @brief Scatter the degrees of freedom, i.e. goes from L-Vector to
+       face E-Vector.
+
+       @param[in]  x The L-vector degrees of freedom.
+       @param[out] y The face E-Vector degrees of freedom with size (face_dofs,
+                     vdim, nf), where nf is the number of interior or boundary
+                     faces requested by @a type in the constructor. The
+                     face_dofs are ordered according to the given
+                     ElementDofOrdering. */
+   void Mult(const Vector &x, Vector &y) const override;
+
+   using FaceRestriction::AddMultTranspose;
+
+   /** @brief Gather the degrees of freedom, i.e. goes from face E-Vector to
+       L-Vector.
+
+       @param[in]     x The face E-Vector degrees of freedom with size
+                        (face_dofs, vdim, nf), where nf is the number of
+                        interior or boundary faces requested by @a type in the
+                        constructor. The face_dofs should be ordered according
+                        to the given ElementDofOrdering
+       @param[in,out] y The L-vector degrees of freedom.
+       @param[in]     a Scalar coefficient for addition. */
+   void AddMultTranspose(const Vector &x, Vector &y,
+                         const real_t a = 1.0) const override;
+
+   const Array<int> &GatherMap() const override;
+};
+
 /** @brief Convert a dof face index from Native ordering to lexicographic
     ordering for quads and hexes.
 

linalg/batched/native.cpp

@@ -146,7 +146,6 @@ void NativeBatchedLinAlg::Invert(DenseTensor &A) const
 
 void NativeBatchedLinAlg::LUFactor(DenseTensor &A, Array<int> &P) const
 {
-   constexpr real_t tol = 0.0; // Make this user-adjustable?
    const int m = A.SizeI();
    const int NE = A.SizeK();
    P.SetSize(m*NE);
@@ -159,52 +158,8 @@ void NativeBatchedLinAlg::LUFactor(DenseTensor &A, Array<int> &P) const
 
    mfem::forall(NE, [=] MFEM_HOST_DEVICE (int e)
    {
-      for (int i = 0; i < m; i++)
-      {
-         // pivoting
-         {
-            int piv = i;
-            real_t a = fabs(data_all(piv,i,e));
-            for (int j = i+1; j < m; j++)
-            {
-               const real_t b = fabs(data_all(j,i,e));
-               if (b > a)
-               {
-                  a = b;
-                  piv = j;
-               }
-            }
-            ipiv_all(i,e) = piv;
-            if (piv != i)
-            {
-               // swap rows i and piv in both L and U parts
-               for (int j = 0; j < m; j++)
-               {
-                  kernels::internal::Swap<real_t>(data_all(i,j,e), data_all(piv,j,e));
-               }
-            }
-         } // pivot end
-
-         if (std::abs(data_all(i,i,e)) <= tol)
-         {
-            d_pivot_flag[0] = false;
-         }
-
-         const real_t a_ii_inv = 1.0 / data_all(i,i,e);
-         for (int j = i+1; j < m; j++)
-         {
-            data_all(j,i,e) *= a_ii_inv;
-         }
-
-         for (int k = i+1; k < m; k++)
-         {
-            const real_t a_ik = data_all(i,k,e);
-            for (int j = i+1; j < m; j++)
-            {
-               data_all(j,k,e) -= a_ik * data_all(j,i,e);
-            }
-         }
-      }
+      const bool flag = kernels::LUFactor(&data_all(0,0,e), m, &ipiv_all(0,e));
+      if (!flag) { d_pivot_flag[0] = false; }
    });
 
    MFEM_VERIFY(pivot_flag.HostRead()[0], "Batch LU factorization failed");

linalg/densemat.cpp

@@ -3487,20 +3487,7 @@ void LUFactors::LSolve(int m, int n, real_t *X) const
    real_t *x = X;
    for (int k = 0; k < n; k++)
    {
-      // X <- P X
-      for (int i = 0; i < m; i++)
-      {
-         mfem::Swap<real_t>(x[i], x[ipiv[i]-ipiv_base]);
-      }
-      // X <- L^{-1} X
-      for (int j = 0; j < m; j++)
-      {
-         const real_t x_j = x[j];
-         for (int i = j+1; i < m; i++)
-         {
-            x[i] -= data[i+j*m] * x_j;
-         }
-      }
+      kernels::LSolve(data, m, ipiv, x);
       x += m;
    }
 }
@@ -3508,17 +3495,9 @@ void LUFactors::LSolve(int m, int n, real_t *X) const
 void LUFactors::USolve(int m, int n, real_t *X) const
 {
    real_t *x = X;
-   // X <- U^{-1} X
    for (int k = 0; k < n; k++)
    {
-      for (int j = m-1; j >= 0; j--)
-      {
-         const real_t x_j = ( x[j] /= data[j+j*m] );
-         for (int i = 0; i < j; i++)
-         {
-            x[i] -= data[i+j*m] * x_j;
-         }
-      }
+      kernels::USolve(data, m, x);
       x += m;
    }
 }
@@ -3661,44 +3640,13 @@ void LUFactors::GetInverseMatrix(int m, real_t *X) const
 void LUFactors::SubMult(int m, int n, int r, const real_t *A21,
                         const real_t *X1, real_t *X2)
 {
-   // X2 <- X2 - A21 X1
-   for (int k = 0; k < r; k++)
-   {
-      for (int j = 0; j < m; j++)
-      {
-         const real_t x1_jk = X1[j+k*m];
-         for (int i = 0; i < n; i++)
-         {
-            X2[i+k*n] -= A21[i+j*n] * x1_jk;
-         }
-      }
-   }
+   kernels::SubMult(m, n, r, A21, X1, X2);
 }
 
 void LUFactors::BlockFactor(
    int m, int n, real_t *A12, real_t *A21, real_t *A22) const
 {
-   // A12 <- L^{-1} P A12
-   LSolve(m, n, A12);
-   // A21 <- A21 U^{-1}
-   for (int j = 0; j < m; j++)
-   {
-      const real_t u_jj_inv = 1.0/data[j+j*m];
-      for (int i = 0; i < n; i++)
-      {
-         A21[i+j*n] *= u_jj_inv;
-      }
-      for (int k = j+1; k < m; k++)
-      {
-         const real_t u_jk = data[j+k*m];
-         for (int i = 0; i < n; i++)
-         {
-            A21[i+k*n] -= A21[i+j*n] * u_jk;
-         }
-      }
-   }
-   // A22 <- A22 - A21 A12
-   SubMult(m, n, n, A21, A12, A22);
+   kernels::BlockFactor(data, m, ipiv, n, A12, A21, A22);
 }
 
 void LUFactors::BlockForwSolve(int m, int n, int r, const real_t *L21,

linalg/densemat.hpp

@@ -1320,7 +1320,6 @@ void BatchLUFactor(DenseTensor &Mlu, Array<int> &P, const real_t TOL = 0.0);
     dimension m x n. */
 void BatchLUSolve(const DenseTensor &Mlu, const Array<int> &P, Vector &X);
 
-
 // Inline methods
 
 inline real_t &DenseMatrix::operator()(int i, int j)

linalg/kernels.hpp

@@ -450,6 +450,30 @@ void MultAtB(const int Aheight, const int Awidth, const int Bwidth,
    AddMultAtB(Aheight, Awidth, Bwidth, Adata, Bdata, AtBdata, TB(1.0), TA(0.0));
 }
 
+/** @brief Multiply the transpose of a matrix of size @a Aheight x @a Awidth
+    and data @a Adata with a matrix of size @a Aheight x @a Bwidth and data @a
+    Bdata: At * B. Add the result to the matrix with data @a AtBdata. */
+template<typename TA, typename TB, typename TC>
+MFEM_HOST_DEVICE inline
+void AddMultAtB(const int Aheight, const int Awidth, const int Bwidth,
+                const TA *Adata, const TB *Bdata, TC *AtBdata)
+{
+   TC *c = AtBdata;
+   for (int i = 0; i < Bwidth; ++i)
+   {
+      for (int j = 0; j < Awidth; ++j)
+      {
+         TC val = 0.0;
+         for (int k = 0; k < Aheight; ++k)
+         {
+            val += Adata[j * Aheight + k] * Bdata[i * Aheight + k];
+         }
+         *c += val;
+         c++;
+      }
+   }
+}
+
 /// Given a matrix of size 2x1, 3x1, or 3x2, compute the left inverse.
 template<int HEIGHT, int WIDTH> MFEM_HOST_DEVICE
 void CalcLeftInverse(const real_t *data, real_t *left_inv);
@@ -1665,24 +1689,21 @@ have_aa:
    return sqrt(fabs(aa))*mult; // take abs before we sort?
 }
 
-
-/// Assuming L.U = P.A for a factored matrix (m x m),
-//  compute x <- A x
+/// @brief Assuming L.U = P.A factored matrix of size (m x m), compute
+/// X <- L^{-1} P X, for a vector X of length m.
 //
 // @param [in] data LU factorization of A
 // @param [in] m square matrix height
-// @param [in] ipiv array storing pivot information
+// @param [in] ipiv array storing pivots
 // @param [in, out] x vector storing right-hand side and then solution
 MFEM_HOST_DEVICE
-inline void LUSolve(const real_t *data, const int m, const int *ipiv,
-                    real_t *x)
+inline void LSolve(const real_t *data, const int m, const int *ipiv, real_t *x)
 {
    // X <- P X
    for (int i = 0; i < m; i++)
    {
       internal::Swap<real_t>(x[i], x[ipiv[i]]);
    }
-
    // X <- L^{-1} X
    for (int j = 0; j < m; j++)
    {
@@ -1692,8 +1713,17 @@ inline void LUSolve(const real_t *data, const int m, const int *ipiv,
          x[i] -= data[i + j * m] * x_j;
       }
    }
+}
 
-   // X <- U^{-1} X
+/// @brief Assuming L.U = P.A factored matrix of size (m x m), compute
+/// X <- U^{-1} X, for a vector X of length m.
+//
+// @param [in] data LU factorization of A
+// @param [in] m square matrix height
+// @param [in, out] x vector storing right-hand side and then solution
+MFEM_HOST_DEVICE
+inline void USolve(const real_t *data, const int m, real_t *x)
+{
    for (int j = m - 1; j >= 0; j--)
    {
       const real_t x_j = (x[j] /= data[j + j * m]);
@@ -1704,6 +1734,148 @@ inline void LUSolve(const real_t *data, const int m, const int *ipiv,
    }
 }
 
+/// @brief Assuming L.U = P.A for a factored matrix (m x m),
+//  compute x <- A x
+//
+// @param [in] data LU factorization of A
+// @param [in] m square matrix height
+// @param [in] ipiv array storing pivot information
+// @param [in, out] x vector storing right-hand side and then solution
+MFEM_HOST_DEVICE
+inline void LUSolve(const real_t *data, const int m, const int *ipiv, real_t *x)
+{
+   LSolve(data, m, ipiv, x);
+   USolve(data, m, x);
+}
+
+
+/// @brief Given an (n x m) matrix A21, compute X2 <- X2 - A21 X1, for matrices
+/// X1, and X2 of size (m x r) and (n x r), respectively.
+MFEM_HOST_DEVICE
+inline void SubMult(const int m, const int n, const int r, const real_t *A21,
+                    const real_t *X1, real_t *X2)
+{
+   // X2 <- X2 - A21 X1
+   for (int k = 0; k < r; k++)
+   {
+      for (int j = 0; j < m; j++)
+      {
+         const real_t x1_jk = X1[j+k*m];
+         for (int i = 0; i < n; i++)
+         {
+            X2[i+k*n] -= A21[i+j*n] * x1_jk;
+         }
+      }
+   }
+}
+
+/// Assuming P.A = L.U factored data of size (m x m), compute the 2x2 block
+/// decomposition:
+///      | P 0 | |  A  A12 | = |  L  0 | | U U12 |
+///      | 0 I | | A21 A22 |   | L21 I | | 0 S22 |
+///   where A12, A21, and A22 are matrices of size (m x n), (n x m), and
+///   (n x n), respectively. The blocks are overwritten as follows:
+///      A12 <- U12 = L^{-1} P A12
+///      A21 <- L21 = A21 U^{-1}
+///      A22 <- S22 = A22 - L21 U12.
+///   The block S22 is the Schur complement.
+MFEM_HOST_DEVICE
+inline void BlockFactor(const real_t *data, int m, const int *ipiv,
+                        int n, real_t *A12, real_t *A21, real_t *A22)
+{
+   // A12 <- L^{-1} P A12
+   for (int i = 0; i < n; ++i)
+   {
+      LSolve(data, m, ipiv, A12 + i*m);
+   }
+   // A21 <- A21 U^{-1}
+   for (int j = 0; j < m; j++)
+   {
+      const real_t u_jj_inv = 1.0/data[j+j*m];
+      for (int i = 0; i < n; i++)
+      {
+         A21[i+j*n] *= u_jj_inv;
+      }
+      for (int k = j+1; k < m; k++)
+      {
+         const real_t u_jk = data[j+k*m];
+         for (int i = 0; i < n; i++)
+         {
+            A21[i+k*n] -= A21[i+j*n] * u_jk;
+         }
+      }
+   }
+   // A22 <- A22 - A21 A12
+   SubMult(m, n, n, A21, A12, A22);
+}
+
+/// @brief Compute the LU factorization of the m x m matrix @a A.
+///
+/// Factorize the matrix of size (m x m) overwriting it with the LU factors. The
+/// factorization is such that L.U = P.A, where A is the original matrix and P
+/// is a permutation matrix represented by ipiv.
+///
+/// @param [in, out] A matrix
+/// @param [in] m size of the square matrix
+/// @param [out] ipiv array of pivots (length m)
+/// @param [in] tol optional fuzzy comparison tolerance. Defaults to 0.0.
+///
+/// @return true if the factorization succeeds, false otherwise (zero pivot).
+MFEM_HOST_DEVICE
+inline bool LUFactor(real_t *A, const int m, int *ipiv, const real_t tol=0.0)
+{
+   bool pivot_flag = true;
+
+   for (int i = 0; i < m; i++)
+   {
+      // pivoting
+      {
+         int piv = i;
+         real_t a = fabs(A[piv + m*i]);
+         for (int j = i+1; j < m; j++)
+         {
+            const real_t b = fabs(A[j + m*i]);
+            if (b > a)
+            {
+               a = b;
+               piv = j;
+            }
+         }
+         ipiv[i] = piv;
+         if (piv != i)
+         {
+            // swap rows i and piv in both L and U parts
+            for (int j = 0; j < m; j++)
+            {
+               internal::Swap<real_t>(A[i + m*j], A[piv + m*j]);
+            }
+         }
+      } // pivot end
+
+      if (abs(A[i + m*i]) <= tol)
+      {
+         pivot_flag = false;
+      }
+
+      const real_t a_ii_inv = 1.0 / A[i + m*i];
+      for (int j = i+1; j < m; j++)
+      {
+         A[j + m*i] *= a_ii_inv;
+      }
+
+      for (int k = i+1; k < m; k++)
+      {
+         const real_t a_ik = A[i + m*k];
+         for (int j = i+1; j < m; j++)
+         {
+            A[j + m*k] -= a_ik * A[j + m*i];
+         }
+      }
+   }
+
+   return pivot_flag;
+}
+
 } // namespace kernels
 
 } // namespace mfem

miniapps/hdiv-linear-solver/grad_div.cpp

@@ -158,6 +158,8 @@ int main(int argc, char *argv[])
 
    if (use_ams)
    {
+      x.ProjectCoefficient(u_vec_coeff);
+
       if (Mpi::Root()) { cout << "\nAMS solver... " << flush; }
       tic_toc.Clear(); tic_toc.Start();
 
@@ -168,8 +170,6 @@ int main(int argc, char *argv[])
 
       OperatorHandle A;
       Vector B, X;
-      b.Assemble();
-      x.ProjectCoefficient(u_vec_coeff);
       a.FormLinearSystem(ess_rt_dofs, x, b, A, X, B);
       HypreParMatrix &Ah = *A.As<HypreParMatrix>();
 
@@ -189,6 +189,9 @@ int main(int argc, char *argv[])
       RT_FECollection fec_rt_lor(order-1, dim, b1, b2_lor);
       ParFiniteElementSpace fes_rt_lor(&mesh, &fec_rt_lor);
 
+      ParGridFunction x_lor(&fes_rt_lor);
+      x_lor.ProjectCoefficient(u_vec_coeff);
+
       ParLinearForm b_lor(&fes_rt_lor);
       b_lor.AddDomainIntegrator(new VectorFEDomainLFIntegrator(f_vec_coeff));
       b_lor.UseFastAssembly(true);
@@ -203,9 +206,6 @@ int main(int argc, char *argv[])
       a.AddDomainIntegrator(new VectorFEMassIntegrator(beta_coeff));
       a.Assemble();
 
-      ParGridFunction x_lor(&fes_rt_lor);
-      x_lor.ProjectCoefficient(u_vec_coeff);
-
       OperatorHandle A;
       Vector B, X;
       a.FormLinearSystem(ess_rt_dofs, x_lor, b_lor, A, X, B);
@@ -222,6 +222,11 @@ int main(int argc, char *argv[])
 
    if (use_hybridization)
    {
+      // Don't include cuBLAS setup time in the hybridization timings
+      GPUBlas::Handle();
+
+      x.ProjectCoefficient(u_vec_coeff);
+
       if (Mpi::Root()) { cout << "\nHybridization solver... " << flush; }
       tic_toc.Clear(); tic_toc.Start();
 
@@ -231,13 +236,12 @@ int main(int argc, char *argv[])
       ParBilinearForm a(&fes_rt);
       a.AddDomainIntegrator(new DivDivIntegrator(alpha_coeff));
       a.AddDomainIntegrator(new VectorFEMassIntegrator(beta_coeff));
+      a.SetAssemblyLevel(AssemblyLevel::ELEMENT);
       a.EnableHybridization(&fes_hb, new NormalTraceJumpIntegrator, ess_rt_dofs);
       a.Assemble();
 
       OperatorHandle A;
       Vector B, X;
-      b.Assemble();
-      x.ProjectCoefficient(u_vec_coeff);
       a.FormLinearSystem(ess_rt_dofs, x, b, A, X, B);
 
       HypreBoomerAMG amg_hb(*A.As<HypreParMatrix>());

tests/unit/fem/test_assembly_levels.cpp

@@ -11,6 +11,7 @@
 
 #include "unit_tests.hpp"
 #include "mfem.hpp"
+#include "linalg/dtensor.hpp"
 #include <fstream>
 #include <iostream>
 
@@ -101,20 +102,35 @@ void test_assembly_level(const char *meshname,
         << ", order=" << order << ", q_order=" << q_order << ", DG=" << dg
         << ", pb=" << getString(pb) << ", assembly=" << getString(assembly));
    Mesh mesh(meshname, 1, 1);
+   mesh.RemoveInternalBoundaries();
    mesh.EnsureNodes();
-   int dim = mesh.Dimension();
+   const int dim = mesh.Dimension();
 
-   FiniteElementCollection *fec;
+   for (int e = 0; e < mesh.GetNE(); ++e)
+   {
+      mesh.SetAttribute(e, 1 + (e % 2));
+   }
+   for (int be = 0; be < mesh.GetNBE(); ++be)
+   {
+      mesh.SetBdrAttribute(be, 1 + (be % 2));
+   }
+   mesh.SetAttributes();
+
+   Array<int> elem_marker({1, 0}), bdr_marker({1, 0});
+   // Periodic meshes = no boundary attributes, don't use markers
+   if (mesh.bdr_attributes.Size() == 0) { bdr_marker.DeleteAll(); }
+
+   std::unique_ptr<FiniteElementCollection> fec;
    if (dg)
    {
-      fec = new L2_FECollection(order, dim, BasisType::GaussLobatto);
+      fec.reset(new L2_FECollection(order, dim, BasisType::GaussLobatto));
    }
    else
    {
-      fec = new H1_FECollection(order, dim);
+      fec.reset(new H1_FECollection(order, dim));
    }
 
-   FiniteElementSpace fespace(&mesh, fec);
+   FiniteElementSpace fespace(&mesh, fec.get());
 
    BilinearForm k_test(&fespace);
    BilinearForm k_ref(&fespace);
@@ -126,12 +142,19 @@ void test_assembly_level(const char *meshname,
    const bool use_ir = q_order_inc > 0;
    const IntegrationRule *ir =
       use_ir ? &IntRules.Get(mesh.GetElementGeometry(0), q_order) : nullptr;
+   const IntegrationRule &ir_face =
+      IntRules.Get(mesh.GetFaceGeometry(0), q_order);
 
    switch (pb)
    {
       case Problem::Mass:
-         k_ref.AddDomainIntegrator(new MassIntegrator(one,ir));
-         k_test.AddDomainIntegrator(new MassIntegrator(one,ir));
+         k_ref.AddDomainIntegrator(new MassIntegrator(one,ir), elem_marker);
+         k_test.AddDomainIntegrator(new MassIntegrator(one,ir), elem_marker);
+         if (!dg && mesh.Conforming() && assembly != AssemblyLevel::FULL)
+         {
+            k_ref.AddBoundaryIntegrator(new MassIntegrator(one, &ir_face), bdr_marker);
+            k_test.AddBoundaryIntegrator(new MassIntegrator(one, &ir_face), bdr_marker);
+         }
          break;
       case Problem::Convection:
          AddConvectionIntegrators(k_ref, vel_coeff, dg);
@@ -168,8 +191,6 @@ void test_assembly_level(const char *meshname,
    y_test -= y_ref;
 
    REQUIRE(y_test.Norml2() < 1.e-12);
-
-   delete fec;
 }
 
 TEST_CASE("H1 Assembly Levels", "[AssemblyLevel], [PartialAssembly], [CUDA]")
@@ -228,6 +249,142 @@ TEST_CASE("H1 Assembly Levels", "[AssemblyLevel], [PartialAssembly], [CUDA]")
    }
 } // H1 Assembly Levels test case
 
+TEST_CASE("H(div) Element Assembly", "[AssemblyLevel][CUDA]")
+{
+   const auto fname = GENERATE(
+                         "../../data/inline-quad.mesh",
+                         "../../data/star-q3.mesh",
+                         "../../data/inline-hex.mesh",
+                         "../../data/fichera-q3.mesh"
+                      );
+   const auto order = GENERATE(1, 2);
+   const auto problem = GENERATE(Problem::Mass, Problem::Diffusion);
+
+   CAPTURE(fname, order, getString(problem));
+
+   Mesh mesh(fname);
+   const int dim = mesh.Dimension();
+   const int ne = mesh.GetNE();
+
+   RT_FECollection fec(order - 1, dim);
+   FiniteElementSpace fes(&mesh, &fec);
+
+   std::unique_ptr<BilinearFormIntegrator> integ;
+   if (problem == Problem::Mass) { integ.reset(new VectorFEMassIntegrator); }
+   else if (problem == Problem::Diffusion) { integ.reset(new DivDivIntegrator); }
+
+   const FiniteElement &fe = *fes.GetFE(0);
+   {
+      ElementTransformation &T = *mesh.GetElementTransformation(0);
+      integ->SetIntegrationRule(MassIntegrator::GetRule(fe, fe, T));
+   }
+
+   const TensorBasisElement *tbe =
+      dynamic_cast<const TensorBasisElement*>(&fe);
+   MFEM_VERIFY(tbe, "");
+   const int ndof = fes.GetFE(0)->GetDof();
+   const Array<int> &dof_map = tbe->GetDofMap();
+
+   Vector ea_data(ne*ndof*ndof);
+   integ->AssembleEA(fes, ea_data, false);
+   const auto ea_mats = Reshape(ea_data.HostRead(), ndof, ndof, ne);
+
+   DenseMatrix elmat;
+   for (int e = 0; e < ne; ++e)
+   {
+      const FiniteElement &el = *fes.GetFE(e);
+      ElementTransformation &T = *mesh.GetElementTransformation(e);
+      integ->AssembleElementMatrix(el, T, elmat);
+
+      for (int i = 0; i < ndof; ++i)
+      {
+         const int ii_s = dof_map[i];
+         const int ii = ii_s >= 0 ? ii_s : -1 - ii_s;
+         const int s_i = ii_s >= 0 ? 1 : -1;
+         for (int j = 0; j < ndof; ++j)
+         {
+            const int jj_s = dof_map[j];
+            const int jj = jj_s >= 0 ? jj_s : -1 - jj_s;
+            const int s_j = jj_s >= 0 ? 1 : -1;
+            elmat(ii, jj) -= s_i*s_j*ea_mats(i, j, e);
+         }
+      }
+
+      REQUIRE(elmat.MaxMaxNorm() == MFEM_Approx(0.0, 1e-10));
+   }
+}
+
+TEST_CASE("NormalTraceJumpIntegrator Element Assembly", "[AssemblyLevel][CUDA]")
+{
+   const auto fname = GENERATE(
+                         "../../data/inline-quad.mesh",
+                         "../../data/star-q3.mesh",
+                         "../../data/inline-hex.mesh",
+                         "../../data/fichera-q3.mesh"
+                      );
+   const int order = GENERATE(1, 2, 3);
+
+   CAPTURE(fname, order);
+
+   Mesh mesh(fname);
+   const int dim = mesh.Dimension();
+
+   RT_FECollection fec(order - 1, dim);
+   FiniteElementSpace fes(&mesh, &fec);
+
+   DG_Interface_FECollection hfec(order - 1, dim);
+   FiniteElementSpace hfes(&mesh, &hfec);
+
+   NormalTraceJumpIntegrator integ;
+
+   const int nf = mesh.GetNFbyType(FaceType::Interior);
+   const int ndof_trial = hfes.GetFaceElement(0)->GetDof();
+   const int ndof_test = fes.GetFE(0)->GetDof();
+   Vector emat(ndof_trial*ndof_test*2*nf);
+   integ.AssembleEAInteriorFaces(hfes, fes, emat, false);
+
+   const TensorBasisElement *tbe =
+      dynamic_cast<const TensorBasisElement*>(fes.GetFE(0));
+   MFEM_VERIFY(tbe, "");
+   const Array<int> &dof_map = tbe->GetDofMap();
+
+   const auto e_mat = Reshape(emat.HostRead(), ndof_test, ndof_trial, 2, nf);
+
+   int fidx = 0;
+   for (int f = 0; f < mesh.GetNumFaces(); ++f)
+   {
+      const Mesh::FaceInformation info = mesh.GetFaceInformation(f);
+      if (!info.IsInterior()) { continue; }
+
+      const int el1 = info.element[0].index;
+      const int el2 = info.element[1].index;
+
+      FaceElementTransformations *FTr = mesh.GetInteriorFaceTransformations(f);
+
+      DenseMatrix elmat;
+      integ.AssembleFaceMatrix(*hfes.GetFaceElement(f),
+                               *fes.GetFE(el1),
+                               *fes.GetFE(el2),
+                               *FTr, elmat);
+      elmat.Threshold(1e-12 * elmat.MaxMaxNorm());
+      for (int ie = 0; ie < 2; ++ie)
+      {
+         for (int i_lex = 0; i_lex < ndof_test; ++i_lex)
+         {
+            const int i_s = dof_map[i_lex];
+            const int i = (i_s >= 0) ? i_s : -1 - i_s;
+            for (int j = 0; j < ndof_trial; ++j)
+            {
+               elmat(i + ie*ndof_test, j) -= e_mat(i_lex, j, ie, fidx);
+            }
+         }
+      }
+      REQUIRE(elmat.MaxMaxNorm() == MFEM_Approx(0.0));
+
+      fidx++;
+   }
+}
+
 TEST_CASE("L2 Assembly Levels", "[AssemblyLevel], [PartialAssembly], [CUDA]")
 {
    const bool dg = true;
@@ -332,44 +489,6 @@ void CompareMatricesNonZeros(SparseMatrix &A1, const SparseMatrix &A2,
    REQUIRE(error == MFEM_Approx(0.0, 1e-10));
 }
 
-#ifdef MFEM_USE_MPI
-
-void CompareMatricesNonZeros(HypreParMatrix &A1, const HypreParMatrix &A2)
-{
-   HYPRE_BigInt *cmap1, *cmap2;
-   SparseMatrix diag1, offd1, diag2, offd2;
-
-   A1.GetDiag(diag1);
-   A2.GetDiag(diag2);
-   A1.GetOffd(offd1, cmap1);
-   A2.GetOffd(offd2, cmap2);
-
-   CompareMatricesNonZeros(diag1, diag2);
-
-   if (cmap1)
-   {
-      std::unordered_map<HYPRE_BigInt,int> cmap2inv;
-      for (int i=0; i<offd2.Width(); ++i) { cmap2inv[cmap2[i]] = i; }
-      CompareMatricesNonZeros(offd1, offd2, cmap1, &cmap2inv);
-   }
-   else
-   {
-      CompareMatricesNonZeros(offd1, offd2);
-   }
-}
-
-void TestSameHypreMatrices(OperatorHandle &A1, OperatorHandle &A2)
-{
-   HypreParMatrix *M1 = A1.Is<HypreParMatrix>();
-   HypreParMatrix *M2 = A2.Is<HypreParMatrix>();
-
-   REQUIRE(M1 != NULL);
-   REQUIRE(M2 != NULL);
-
-   CompareMatricesNonZeros(*M1, *M2);
-   CompareMatricesNonZeros(*M2, *M1);
-}
-
 void TestSameSparseMatrices(OperatorHandle &A1, OperatorHandle &A2)
 {
    SparseMatrix *M1 = A1.Is<SparseMatrix>();
@@ -445,6 +564,44 @@ TEST_CASE("Serial H1 Full Assembly", "[AssemblyLevel], [CUDA]")
    REQUIRE(B1.Normlinf() == MFEM_Approx(0.0));
 }
 
+#ifdef MFEM_USE_MPI
+
+void CompareMatricesNonZeros(HypreParMatrix &A1, const HypreParMatrix &A2)
+{
+   HYPRE_BigInt *cmap1, *cmap2;
+   SparseMatrix diag1, offd1, diag2, offd2;
+
+   A1.GetDiag(diag1);
+   A2.GetDiag(diag2);
+   A1.GetOffd(offd1, cmap1);
+   A2.GetOffd(offd2, cmap2);
+
+   CompareMatricesNonZeros(diag1, diag2);
+
+   if (cmap1)
+   {
+      std::unordered_map<HYPRE_BigInt,int> cmap2inv;
+      for (int i=0; i<offd2.Width(); ++i) { cmap2inv[cmap2[i]] = i; }
+      CompareMatricesNonZeros(offd1, offd2, cmap1, &cmap2inv);
+   }
+   else
+   {
+      CompareMatricesNonZeros(offd1, offd2);
+   }
+}
+
+void TestSameHypreMatrices(OperatorHandle &A1, OperatorHandle &A2)
+{
+   HypreParMatrix *M1 = A1.Is<HypreParMatrix>();
+   HypreParMatrix *M2 = A2.Is<HypreParMatrix>();
+
+   REQUIRE(M1 != NULL);
+   REQUIRE(M2 != NULL);
+
+   CompareMatricesNonZeros(*M1, *M2);
+   CompareMatricesNonZeros(*M2, *M1);
+}
+
 TEST_CASE("Parallel H1 Full Assembly", "[AssemblyLevel], [Parallel], [CUDA]")
 {
    auto order = GENERATE(1, 2, 3);

tests/unit/fem/test_bilinearform.cpp

@@ -142,3 +142,28 @@ TEST_CASE("FormLinearSystem/SolutionScope",
       REQUIRE(AsConst(sol)(bdr_dof) == 0.0);
    }
 }
+
+TEST_CASE("GetElementMatrices", "[BilinearForm]")
+{
+   const int order = 3;
+   Mesh mesh = Mesh::MakeCartesian2D(3, 3, Element::QUADRILATERAL);
+   H1_FECollection fec(order, mesh.Dimension());
+   FiniteElementSpace fes(&mesh, &fec);
+
+   BilinearForm a(&fes);
+   a.AddDomainIntegrator(new MassIntegrator);
+   const DenseTensor &el_mat = a.GetElementMatrices();
+
+   BilinearForm a_ea(&fes);
+   a_ea.AddDomainIntegrator(new MassIntegrator);
+   a_ea.SetAssemblyLevel(AssemblyLevel::ELEMENT);
+   const DenseTensor &el_mat_ea = a_ea.GetElementMatrices();
+
+   for (int e = 0; e < mesh.GetNE(); ++e)
+   {
+      DenseMatrix m = el_mat(e);
+      const DenseMatrix &m_ea = el_mat_ea(e);
+      m -= m_ea;
+      REQUIRE(m.MaxMaxNorm() == MFEM_Approx(0.0));
+   }
+}