Fix averaging of face penalty factor in LOR DG

Returning Array instead of pointer allows use on device

fem/CMakeLists.txt

@@ -229,8 +229,13 @@ set(HDRS
   lor/lor_ams.hpp
   lor/lor_batched.hpp
   lor/lor_h1.hpp
+  lor/lor_dg.hpp
   lor/lor_nd.hpp
   lor/lor_rt.hpp
+  lor/lor_h1_impl.hpp
+  lor/lor_dg_impl.hpp
+  lor/lor_nd_impl.hpp
+  lor/lor_rt_impl.hpp
   lor/lor_util.hpp
   multigrid.hpp
   nonlinearform.hpp

fem/bilininteg.hpp

@@ -3441,6 +3441,8 @@ public:
 
    const IntegrationRule &GetRule(int order, Geometry::Type geom);
 
+   real_t GetPenaltyParameter() const { return kappa; }
+
 private:
    void SetupPA(const FiniteElementSpace &fes, FaceType type);
 };

fem/fe/face_map_utils.hpp

@@ -259,6 +259,30 @@ inline void FaceIdxToVolIdx3D(const int index, const int size1d,
    i = yz_plane ? level : _i;
 }
 
+MFEM_HOST_DEVICE
+inline int FaceIdxToVolIdx(int dim, int i, int size1d, int face0, int face1,
+                           int side, int orientation)
+{
+   if (dim == 2)
+   {
+      int ix, iy;
+      internal::FaceIdxToVolIdx2D(i, size1d, face0, face1, side, ix, iy);
+      return ix + iy*size1d;
+   }
+   else if (dim == 3)
+   {
+      int ix, iy, iz;
+      internal::FaceIdxToVolIdx3D(i, size1d, face0, face1, side, orientation,
+                                  ix, iy, iz);
+      return ix + size1d*iy + size1d*size1d*iz;
+   }
+   else
+   {
+      MFEM_ABORT_KERNEL("Invalid dimension");
+      return -1;
+   }
+};
+
 } // namespace internal
 
 } // namespace mfem

fem/lininteg.hpp

@@ -544,7 +544,7 @@ public:
     Specifically, given the Dirichlet data $u_D$, the linear form assembles the
     following integrals on the boundary:
    $$
-    \sigma \langle u_D, (Q \nabla v)) \cdot n \rangle + \kappa \langle {h^{-1} Q} u_D, v \rangle,
+    \sigma \langle u_D, (Q \nabla v) \cdot n \rangle + \kappa \langle {h^{-1} Q} u_D, v \rangle,
    $$
     where Q is a scalar or matrix diffusion coefficient and v is the test
     function. The parameters $\sigma$ and $\kappa$ should be the same as the ones

fem/lor/lor_batched.cpp

@@ -14,9 +14,11 @@
 #include "../../general/forall.hpp"
 #include <climits>
 #include "../pbilinearform.hpp"
+#include "../../fem/fe/face_map_utils.hpp"
 
 // Specializations
 #include "lor_h1.hpp"
+#include "lor_dg.hpp"
 #include "lor_nd.hpp"
 #include "lor_rt.hpp"
 
@@ -54,17 +56,18 @@ bool BatchedLORAssembly::FormIsSupported(BilinearForm &a)
    // Batched LOR requires all tensor elements
    if (!UsesTensorBasis(*a.FESpace())) { return false; }
 
-   if (dynamic_cast<const H1_FECollection*>(fec))
+   if (dynamic_cast<const H1_FECollection*>(fec) ||
+       dynamic_cast<const DG_FECollection*>(fec))
    {
-      if (HasIntegrators<DiffusionIntegrator, MassIntegrator>(a)) { return true; }
+      return HasIntegrators<DiffusionIntegrator, MassIntegrator>(a);
    }
    else if (dynamic_cast<const ND_FECollection*>(fec))
    {
-      if (HasIntegrators<CurlCurlIntegrator, VectorFEMassIntegrator>(a)) { return true; }
+      return HasIntegrators<CurlCurlIntegrator, VectorFEMassIntegrator>(a);
    }
    else if (dynamic_cast<const RT_FECollection*>(fec))
    {
-      if (HasIntegrators<DivDivIntegrator, VectorFEMassIntegrator>(a)) { return true; }
+      return HasIntegrators<DivDivIntegrator, VectorFEMassIntegrator>(a);
    }
    return false;
 }
@@ -75,12 +78,14 @@ void BatchedLORAssembly::FormLORVertexCoordinates(FiniteElementSpace &fes_ho,
    Mesh &mesh_ho = *fes_ho.GetMesh();
    mesh_ho.EnsureNodes();
 
+   const bool dg = fes_ho.IsDGSpace();
+
    // Get nodal points at the LOR vertices
    const int dim = mesh_ho.Dimension();
    const int sdim = mesh_ho.SpaceDimension();
    const int nel_ho = mesh_ho.GetNE();
    const int order = fes_ho.GetMaxElementOrder();
-   const int nd1d = order + 1;
+   const int nd1d = dg ? order + 2 : order + 1;
    const int ndof_per_el = static_cast<int>(pow(nd1d, dim));
 
    const GridFunction *nodal_gf = mesh_ho.GetNodes();
@@ -92,7 +97,8 @@ void BatchedLORAssembly::FormLORVertexCoordinates(FiniteElementSpace &fes_ho,
    Vector nodal_evec(nodal_restriction->Height());
    nodal_restriction->Mult(*nodal_gf, nodal_evec);
 
-   IntegrationRule ir = GetCollocatedIntRule(fes_ho);
+   const IntegrationRule ir = GetLobattoIntRule(
+                                 mesh_ho.GetTypicalElementGeometry(), nd1d);
 
    // Map from nodal E-vector to Q-vector at the LOR vertex points
    X_vert.SetSize(sdim*ndof_per_el*nel_ho);
@@ -159,6 +165,7 @@ int BatchedLORAssembly::FillI(SparseMatrix &A) const
    const auto K = dof_glob2loc_offsets_.Read();
    const auto map = Reshape(sparse_mapping.Read(), nnz_per_row, ndof_per_el);
 
+
    auto I = A.WriteI();
 
    mfem::forall(nvdof + 1, [=] MFEM_HOST_DEVICE (int ii) { I[ii] = 0; });
@@ -358,6 +365,174 @@ void BatchedLORAssembly::FillJAndData(SparseMatrix &A) const
    });
 }
 
+void BatchedLORAssembly::SparseIJToCSR_DG(OperatorHandle &A) const
+{
+   const int ndof_per_el = fes_ho.GetFE(0)->GetDof();
+   const int nel_ho = fes_ho.GetNE();
+   const int nnz_per_row = sparse_ij.Size()/ndof_per_el/nel_ho;
+   const int dim = fes_ho.GetMesh()->Dimension();
+   const int nrows = nel_ho*ndof_per_el;
+   const int p = fes_ho.GetMaxElementOrder();
+   const int pp1 = p + 1;
+   const int nnz = nrows*nnz_per_row;
+
+   const int face_nbr_vsize = [this]()
+   {
+#ifdef MFEM_USE_MPI
+      if (auto *par_fes = dynamic_cast<ParFiniteElementSpace*>(&fes_ho))
+      {
+         return par_fes->GetFaceNbrVSize();
+      }
+#endif
+      return 0;
+   }();
+
+   // If A contains an existing SparseMatrix, reuse it (and try to reuse its
+   // I, J, A arrays if they are big enough)
+   SparseMatrix *A_mat = A.Is<SparseMatrix>();
+   if (!A_mat)
+   {
+      A_mat = new SparseMatrix;
+      A.Reset(A_mat);
+   }
+
+   A_mat->OverrideSize(nrows, nrows + face_nbr_vsize);
+
+   EnsureCapacity(A_mat->GetMemoryI(), nrows + 1);
+   EnsureCapacity(A_mat->GetMemoryJ(), nnz);
+   EnsureCapacity(A_mat->GetMemoryData(), nnz);
+
+   Array<int> nbr_info(nel_ho*3*2*dim);
+   auto h_nbr_info = Reshape(nbr_info.HostWrite(), nel_ho, 2*dim, 3);
+   const int num_faces = fes_ho.GetMesh()->GetNumFaces();
+   for (int f = 0; f < num_faces; f++)
+   {
+      Mesh::FaceInformation finfo = fes_ho.GetMesh()->GetFaceInformation(f);
+      int e0 = finfo.element[0].index;
+      int f0 = finfo.element[0].local_face_id;
+      if (finfo.IsBoundary())
+      {
+         h_nbr_info(e0,f0,0) = -1;
+         h_nbr_info(e0,f0,1)= -1;
+         h_nbr_info(e0,f0,2)= -1;
+      }
+      else if (finfo.IsShared())
+      {
+         // Face neighbors elements are indexed after the last local element
+         h_nbr_info(e0,f0,0) = nel_ho + finfo.element[1].index;
+         h_nbr_info(e0,f0,1)= finfo.element[1].orientation;
+         h_nbr_info(e0,f0,2)= finfo.element[1].local_face_id;
+      }
+      else if (finfo.IsInterior())
+      {
+         int e1 = finfo.element[1].index;
+         int f1 = finfo.element[1].local_face_id;
+         h_nbr_info(e0,f0,0) = e1;
+         h_nbr_info(e0,f0,1)= finfo.element[1].orientation;
+         h_nbr_info(e0,f0,2)= f1;
+         h_nbr_info(e1,f1,0) = e0;
+         h_nbr_info(e1,f1,1) = finfo.element[1].orientation;
+         h_nbr_info(e1,f1,2) = f0;
+      }
+   };
+
+   auto h_I = A_mat->HostWriteI();
+   h_I[0] = 0;
+   for (int i = 0; i < nrows; ++i)
+   {
+      const int iel_ho = i / ndof_per_el;
+      const int iloc = i % ndof_per_el;
+      static const int lex_map_2[4] = {3, 1, 0, 2};
+      static const int lex_map_3[6] = {4, 2, 1, 3, 0, 5};
+      const int local_i[3] = {iloc % pp1, (iloc/pp1)%pp1, iloc/pp1/pp1};
+      int bdr_count = 0;
+      for (int n_idx = 0; n_idx < dim; ++n_idx)
+      {
+         for (int e_i = 0; e_i < 2; ++e_i)
+         {
+            const int j_lex = e_i + n_idx*2;
+            const int f = (dim == 3) ? lex_map_3[j_lex]:lex_map_2[j_lex];
+            const bool boundary = (local_i[n_idx] == e_i * p);
+            if (boundary)
+            {
+               int neighbor_idx = h_nbr_info(iel_ho, f, 0);
+               if (neighbor_idx == -1)
+               {
+                  ++bdr_count;
+               }
+            }
+         }
+      }
+      h_I[i+1] = h_I[i] + (nnz_per_row - bdr_count);
+   }
+
+   const auto V = Reshape(sparse_ij.Read(), nnz_per_row, ndof_per_el, nel_ho);
+   auto J = A_mat->WriteJ();
+   auto AV = A_mat->WriteData();
+   auto I = A_mat->ReadI();
+
+   auto d_nbr_info = Reshape(nbr_info.Read(), nel_ho, 2*dim, 3);
+   mfem::forall(nrows, [=] MFEM_HOST_DEVICE (int i)
+   {
+      const int e = i / ndof_per_el;
+      const int iloc = i % ndof_per_el;
+      const int local_x = iloc % pp1;
+      const int local_y = (iloc/pp1)%pp1;
+      const int local_z = iloc/pp1/pp1;
+      const int local_i[3] = {local_x, local_y, local_z};
+      int offset = I[i];
+      static const int lex_map_2[4] = {3, 1, 0, 2};
+      static const int lex_map_3[6] = {4,2,1,3,0,5};
+      const int *lex_map = (dim == 2) ? lex_map_2 : lex_map_3;
+      AV[offset] = V(0, iloc, e);
+      J[offset] = i;
+      ++offset;
+      for (int n_idx = 0; n_idx < dim; ++n_idx)
+      {
+         // qi is the face lexicographic index, obtained by taking the
+         // lexicographic index of the coordinates ommiting n_idx.
+         int qi = 0;
+         int stride = 1;
+         for (int d = 0; d < dim; ++d)
+         {
+            if (d != n_idx)
+            {
+               qi += local_i[d]*stride;
+               stride *= pp1;
+            }
+         }
+         for (int e_i = 0; e_i < 2; ++e_i)
+         {
+            const int j_lex = e_i + n_idx*2;
+            const int f = lex_map[j_lex];
+            const bool bdr = (local_i[n_idx] == e_i * p);
+            if (bdr)
+            {
+               const int nbr_e = d_nbr_info(e, f, 0);
+               const int nbr_ori = d_nbr_info(e, f, 1);
+               const int nbr_f = d_nbr_info(e, f, 2);
+               if (nbr_e != -1)
+               {
+                  const int nbr_loc_idx = internal::FaceIdxToVolIdx(
+                                             dim, qi, pp1, f, nbr_f, 1, nbr_ori);
+                  J[offset] = nbr_e*ndof_per_el + nbr_loc_idx;
+                  AV[offset] = V(f+1, iloc, e);
+                  ++offset;
+               }
+            }
+            else
+            {
+               int shift = (e_i == 0) ? -1 : 1;
+               for (int n = 0; n < n_idx; ++n) { shift *= pp1; }
+               J[offset] = i + shift;
+               AV[offset] = V(f+1, iloc, e);
+               ++offset;
+            }
+         }
+      }
+   });
+}
+
 void BatchedLORAssembly::SparseIJToCSR(OperatorHandle &A) const
 {
    const int nvdof = fes_ho.GetVSize();
@@ -372,12 +547,11 @@ void BatchedLORAssembly::SparseIJToCSR(OperatorHandle &A) const
    }
 
    A_mat->OverrideSize(nvdof, nvdof);
+   EnsureCapacity(A_mat->GetMemoryI(), nvdof + 1);
 
-   A_mat->GetMemoryI().New(nvdof+1, Device::GetDeviceMemoryType());
-   int nnz = FillI(*A_mat);
-
-   A_mat->GetMemoryJ().New(nnz, Device::GetDeviceMemoryType());
-   A_mat->GetMemoryData().New(nnz, Device::GetDeviceMemoryType());
+   const int nnz = FillI(*A_mat);
+   EnsureCapacity(A_mat->GetMemoryJ(), nnz);
+   EnsureCapacity(A_mat->GetMemoryData(), nnz);
    FillJAndData(*A_mat);
 }
 
@@ -431,6 +605,19 @@ void BatchedLORAssembly::AssembleWithoutBC(BilinearForm &a, OperatorHandle &A)
    // Assemble the matrix, depending on what the form is.
    // This fills in the arrays sparse_ij and sparse_mapping.
    const FiniteElementCollection *fec = fes_ho.FEColl();
+
+   // Handle DG case separately, because assembly of CSR matrix requires
+   // handling face terms.
+   if (dynamic_cast<const DG_FECollection*>(fec))
+   {
+      if (HasIntegrators<DiffusionIntegrator, MassIntegrator>(a))
+      {
+         AssemblyKernel<BatchedLOR_DG>(a);
+      }
+      SparseIJToCSR_DG(A);
+      return;
+   }
+
    if (dynamic_cast<const H1_FECollection*>(fec))
    {
       if (HasIntegrators<DiffusionIntegrator, MassIntegrator>(a))
@@ -453,10 +640,47 @@ void BatchedLORAssembly::AssembleWithoutBC(BilinearForm &a, OperatorHandle &A)
       }
    }
 
-   return SparseIJToCSR(A);
+   SparseIJToCSR(A);
 }
 
 #ifdef MFEM_USE_MPI
+void BatchedLORAssembly::ParAssemble_DG(SparseMatrix &A_local,
+                                        OperatorHandle &A)
+{
+   auto &par_fes = static_cast<ParFiniteElementSpace&>(fes_ho);
+
+   // handle the case when 'a' contains off-diagonal
+   const int lvsize = par_fes.GetVSize();
+   const Array<HYPRE_BigInt> &face_nbr_glob_ldof =
+      par_fes.GetFaceNbrGlobalDofMapArray();
+   const HYPRE_BigInt ldof_offset = par_fes.GetMyDofOffset();
+
+   const int nnz_local = A_local.NumNonZeroElems();
+   Array<HYPRE_BigInt> glob_J(nnz_local);
+
+   const HYPRE_BigInt *d_face_nbr_glob_ldof = face_nbr_glob_ldof.Read();
+   const int *d_J = A_local.ReadJ();
+   HYPRE_BigInt *d_glob_J = glob_J.Write();
+
+   mfem::forall(nnz_local, [=] MFEM_HOST_DEVICE (int i)
+   {
+      if (d_J[i] < lvsize)
+      {
+         d_glob_J[i] = d_J[i] + ldof_offset;
+      }
+      else
+      {
+         d_glob_J[i] = d_face_nbr_glob_ldof[d_J[i] - lvsize];
+      }
+   });
+
+   A.Reset(new HypreParMatrix(
+              par_fes.GetComm(), lvsize, par_fes.GlobalVSize(),
+              par_fes.GlobalVSize(), A_local.HostReadWriteI(),
+              glob_J.HostReadWrite(), A_local.HostReadWriteData(),
+              par_fes.GetDofOffsets(), par_fes.GetDofOffsets()));
+}
+
 void BatchedLORAssembly::ParAssemble(
    BilinearForm &a, const Array<int> &ess_dofs, OperatorHandle &A)
 {
@@ -464,13 +688,18 @@ void BatchedLORAssembly::ParAssemble(
    OperatorHandle A_local;
    AssembleWithoutBC(a, A_local);
 
-   ParBilinearForm *pa =
-      dynamic_cast<ParBilinearForm*>(&a);
-
-   pa->ParallelRAP(*A_local.As<SparseMatrix>(), A, true);
-
-   A.As<HypreParMatrix>()->EliminateBC(ess_dofs,
-                                       Operator::DiagonalPolicy::DIAG_ONE);
+   if (dynamic_cast<const DG_FECollection*>(fes_ho.FEColl()))
+   {
+      ParAssemble_DG(*A_local.As<SparseMatrix>(), A);
+   }
+   else
+   {
+      ParBilinearForm *pa =
+         dynamic_cast<ParBilinearForm*>(&a);
+      pa->ParallelRAP(*A_local.As<SparseMatrix>(), A, true);
+      A.As<HypreParMatrix>()->EliminateBC(ess_dofs,
+                                          Operator::DiagonalPolicy::DIAG_ONE);
+   }
 }
 #endif
 
@@ -497,12 +726,22 @@ BatchedLORAssembly::BatchedLORAssembly(FiniteElementSpace &fes_ho_)
    FormLORVertexCoordinates(fes_ho, X_vert);
 }
 
-IntegrationRule GetCollocatedIntRule(FiniteElementSpace &fes)
+IntegrationRule GetLobattoIntRule(Geometry::Type geom, int nd1d)
 {
    IntegrationRules irs(0, Quadrature1D::GaussLobatto);
-   const Geometry::Type geom = fes.GetMesh()->GetTypicalElementGeometry();
-   const int nd1d = fes.GetMaxElementOrder() + 1;
    return irs.Get(geom, 2*nd1d - 3);
 }
 
+IntegrationRule GetCollocatedIntRule(FiniteElementSpace &fes)
+{
+   const Geometry::Type geom = fes.GetMesh()->GetTypicalElementGeometry();
+   return GetLobattoIntRule(geom, fes.GetMaxElementOrder() + 1);
+}
+
+IntegrationRule GetCollocatedFaceIntRule(FiniteElementSpace &fes)
+{
+   const Geometry::Type geom = fes.GetMesh()->GetTypicalFaceGeometry();
+   return GetLobattoIntRule(geom, fes.GetMaxElementOrder() + 1);
+}
+
 } // namespace mfem

fem/lor/lor_batched.hpp

@@ -25,6 +25,7 @@ namespace mfem
 /// supported, currently:
 ///
 ///  - H1 diffusion + mass
+///  - DG diffusion + mass (in progress)
 ///  - ND curl-curl + mass
 ///  - RT div-div + mass
 ///
@@ -73,6 +74,9 @@ public:
    /// Return the vertices of the LOR mesh in E-vector format
    const Vector &GetLORVertexCoordinates() { return X_vert; }
 
+   /// Specialized implementation of SparseIJToCSR for DG spaces.
+   void SparseIJToCSR_DG(OperatorHandle &A) const;
+
 protected:
    /// After assembling the "sparse IJ" format, convert it to CSR.
    void SparseIJToCSR(OperatorHandle &A) const;
@@ -105,6 +109,9 @@ public:
    void FillJAndData(SparseMatrix &A) const;
 
 #ifdef MFEM_USE_MPI
+   /// Assemble the parallel DG matrix (with shared faces).
+   void ParAssemble_DG(SparseMatrix &A_local, OperatorHandle &A);
+
    /// Assemble the system in parallel and place the result in @a A.
    void ParAssemble(BilinearForm &a, const Array<int> &ess_dofs,
                     OperatorHandle &A);
@@ -128,9 +135,8 @@ void EnsureCapacity(Memory<T> &mem, int capacity)
 
 /// Return the first domain integrator in the form @a i of type @a T.
 template <typename T>
-static T *GetIntegrator(BilinearForm &a)
+static T *GetIntegrator(Array<BilinearFormIntegrator*> *integs)
 {
-   Array<BilinearFormIntegrator*> *integs = a.GetDBFI();
    if (integs != NULL)
    {
       for (auto *i : *integs)
@@ -144,8 +150,32 @@ static T *GetIntegrator(BilinearForm &a)
    return nullptr;
 }
 
+template <typename T>
+static T *GetIntegrator(BilinearForm &a)
+{
+   return GetIntegrator<T>(a.GetDBFI());
+}
+
+template <typename T>
+static T *GetInteriorFaceIntegrator(BilinearForm &a)
+{
+   return GetIntegrator<T>(a.GetFBFI());
+}
+
+/// @brief Return the Gauss-Lobatto rule for geometry @a geom with @a nd1d
+/// points per dimension.
+IntegrationRule GetLobattoIntRule(Geometry::Type geom, int nd1d);
+
+/// @brief Return the Gauss-Lobatto rule collocated with the element nodes.
+///
+/// Assumes @a fes uses Gauss-Lobatto basis.
 IntegrationRule GetCollocatedIntRule(FiniteElementSpace &fes);
 
+/// @brief Return the Gauss-Lobatto rule collocated with face nodes.
+///
+/// Assumes @a fes uses Gauss-Lobatto basis.
+IntegrationRule GetCollocatedFaceIntRule(FiniteElementSpace &fes);
+
 template <typename INTEGRATOR>
 void ProjectLORCoefficient(BilinearForm &a, CoefficientVector &coeff_vector)
 {

fem/lor/lor_dg.hpp

@@ -0,0 +1,79 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_LOR_DG
+#define MFEM_LOR_DG
+
+#include "lor_batched.hpp"
+
+namespace mfem
+{
+
+// BatchedLORKernel specialization for DG spaces. Not user facing. See the
+// classes BatchedLORAssembly and BatchedLORKernel .
+class BatchedLOR_DG : BatchedLORKernel
+{
+   IntegrationRule ir_face; ///< Collocated Gauss-Lobatto face quadrature rule.
+   real_t kappa; ///< DG penalty parameter.
+public:
+   template <int ORDER, int SDIM> void Assemble2D();
+   template <int ORDER> void Assemble3D();
+   BatchedLOR_DG(BilinearForm &a,
+                 FiniteElementSpace &fes_ho_,
+                 Vector &X_vert_,
+                 Vector &sparse_ij_,
+                 Array<int> &sparse_mapping_)
+      : BatchedLORKernel(fes_ho_, X_vert_, sparse_ij_, sparse_mapping_),
+        ir_face(GetLobattoIntRule(fes_ho_.GetMesh()->GetTypicalFaceGeometry(),
+                                  fes_ho_.GetMaxElementOrder() + 1))
+   {
+      ProjectLORCoefficient<MassIntegrator>(a, c1);
+      ProjectLORCoefficient<DiffusionIntegrator>(a, c2);
+
+      auto *integ = GetInteriorFaceIntegrator<DGDiffusionIntegrator>(a);
+      if (integ)
+      {
+         kappa = integ->GetPenaltyParameter();
+      }
+      else
+      {
+         kappa = 0.0;
+      }
+   }
+
+   /// @brief Compute and return the face info array.
+   ///
+   /// The face info array has shape (6, nf), where @a nf is the number of
+   /// faces. For each face @a i, the column (:,i) has entries (e0, f0, o0, e1,
+   /// f1, o1), where @a e is adjacent element, @a f is the local face index,
+   /// and @a o is the orientation. For boundary and shared faces, (e1, f1, o1)
+   /// are all set to -1.
+   Array<int> GetFaceInfo() const;
+
+   /// @brief Compute and return the boundary penalty factor.
+   ///
+   /// The returned vector has shape (nq, nf), where @a nq is the number of
+   /// nodes per face, and @a nf is the number of faces.
+   ///
+   /// The boundary penalty factor is $J_f / h = J_f^2 / J_e$ (since $h = J_e /
+   /// J_f$), where $J_f$ is the face Jacobian determinant, and $J_e$ is the
+   /// element Jacobian determinant.
+   Vector GetBdrPenaltyFactor() const;
+
+   /// Assemble the face penalty terms in the matrix @a sparse_ij.
+   void AssembleFaceTerms();
+};
+
+}
+
+#include "lor_dg_impl.hpp"
+
+#endif

fem/lor/lor_dg_impl.hpp

@@ -0,0 +1,391 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "lor_util.hpp"
+#include "../../linalg/dtensor.hpp"
+#include "../../general/forall.hpp"
+#include "../../mesh/face_nbr_geom.hpp"
+#include "lor_dg.hpp"
+
+namespace mfem
+{
+
+Array<int> BatchedLOR_DG::GetFaceInfo() const
+{
+   Mesh &mesh = *fes_ho.GetMesh();
+   const int nf = mesh.GetNumFaces();
+   Array<int> face_info(nf * 6); // (e0, f0, o0, e1, f1, o1)
+   auto h_face_info = Reshape(face_info.HostWrite(), 6, nf);
+   for (int f = 0; f < nf; ++f)
+   {
+      auto finfo = mesh.GetFaceInformation(f);
+      h_face_info(0, f) = finfo.element[0].index;
+      h_face_info(1, f) = finfo.element[0].local_face_id;
+      h_face_info(2, f) = finfo.element[0].orientation;
+      if (finfo.IsLocal()) // Interior, non-shared face
+      {
+         h_face_info(3, f) = finfo.element[1].index;
+         h_face_info(4, f) = finfo.element[1].local_face_id;
+         h_face_info(5, f) = finfo.element[1].orientation;
+      }
+      else
+      {
+         h_face_info(3, f) = -1;
+         h_face_info(4, f) = -1;
+         h_face_info(5, f) = -1;
+      }
+   }
+   return face_info;
+}
+
+Vector BatchedLOR_DG::GetBdrPenaltyFactor() const
+{
+   Mesh &mesh = *fes_ho.GetMesh();
+
+   const int nf = mesh.GetNumFaces();
+   Array<int> f_int(mesh.GetNFbyType(FaceType::Interior));
+   Array<int> f_bdr(mesh.GetNFbyType(FaceType::Boundary));
+   {
+      int i_int = 0;
+      int i_bdr = 0;
+      for (int i = 0; i < nf; ++i)
+      {
+         const auto f = mesh.GetFaceInformation(i);
+         if (f.IsBoundary())
+         {
+            f_bdr[i_bdr] = i;
+            ++i_bdr;
+         }
+         else if (f.IsInterior())
+         {
+            f_int[i_int] = i;
+            ++i_int;
+         }
+      }
+   }
+
+   const auto geom = fes_ho.GetMesh()->GetGeometricFactors(
+                        ir, GeometricFactors::DETERMINANTS);
+
+   const int nq = ir_face.Size();
+   Vector face_Jh(nq * nf);
+   for (const FaceType ft : {FaceType::Interior, FaceType::Boundary})
+   {
+      const int nft = mesh.GetNFbyType(ft);
+      auto *geom_face = mesh.GetFaceGeometricFactors(
+                           ir_face, FaceGeometricFactors::DETERMINANTS, ft);
+
+      const L2FaceValues fv = (ft == FaceType::Interior)
+                              ? L2FaceValues::DoubleValued
+                              : L2FaceValues::SingleValued;
+      const int m = (fv == L2FaceValues::DoubleValued) ? 2 : 1;
+
+      auto *r = fes_ho.GetFaceRestriction(ElementDofOrdering::LEXICOGRAPHIC, ft, fv);
+      Vector detJ_r(nq * m * nft);
+      r->Mult(geom->detJ, detJ_r);
+
+      const auto *d_i = (ft == FaceType::Interior) ? f_int.Read() : f_bdr.Read();
+      const auto d_detJ_face = Reshape(geom_face->detJ.Read(), nq, nft);
+      const auto d_detJ_r = Reshape(detJ_r.Read(), nq, m, nft);
+      auto d_face_Jh = Reshape(face_Jh.Write(), nq, nf);
+
+      mfem::forall(nft * nq, [=] MFEM_HOST_DEVICE (int ii)
+      {
+         const int i = ii % nq;
+         const int f = ii / nq;
+         const real_t J_el = 0.5*(d_detJ_r(i, 0, f) + d_detJ_r(i, m==2?1:0, f));
+         const real_t J_f = d_detJ_face(i, f);
+         d_face_Jh(i, d_i[f]) = J_f * J_f / J_el;
+      });
+   }
+   return face_Jh;
+}
+
+void BatchedLOR_DG::AssembleFaceTerms()
+{
+   Mesh &mesh = *fes_ho.GetMesh();
+
+   const int nnz_per_row = 1 + mesh.Dimension()*2;
+   const int pp1 = fes_ho.GetMaxElementOrder() + 1;
+   const int nel_ho = mesh.GetNE();
+   const int nf = mesh.GetNumFaces();
+   const int nd_face = ir_face.Size();
+   const int nd = ir.Size();
+   const int dim = mesh.Dimension();
+
+   Array<int> face_info = GetFaceInfo();
+   const auto d_face_info = Reshape(face_info.Read(), 6, nf);
+
+   Vector face_Jh = GetBdrPenaltyFactor();
+   const auto d_face_Jh = Reshape(face_Jh.Read(), nd_face, nf);
+
+   const auto *w_face = ir_face.GetWeights().Read();
+
+   // Penalty parameter (avoid capturing *this in lambda)
+   const real_t d_kappa = kappa;
+
+   // Get diffusion coefficient
+   const bool const_dq = c2.Size() == 1;
+   const auto DQ = const_dq?Reshape(c2.Read(),1,1):Reshape(c2.Read(),nd,nel_ho);
+
+   // Sparse matrix entries
+   auto V = Reshape(sparse_ij.ReadWrite(), nnz_per_row, nd, nel_ho);
+
+   mfem::forall(nf, [=] MFEM_HOST_DEVICE (int f)
+   {
+      const int f_0 = d_face_info(1, f);
+      const int f_1 = d_face_info(4, f);
+      const int nsides = (f_1 >= 0) ? 2 : 1;
+      for (int el_i = 0; el_i < nsides; ++el_i)
+      {
+         const int e = d_face_info(3*el_i, f);
+         const int o = d_face_info(3*el_i + 2, f);
+         const int v_idx = 1 + ((el_i == 0) ? f_0 : f_1);
+         for (int i = 0; i < nd_face; ++i)
+         {
+            const int ii = internal::FaceIdxToVolIdx(dim, i, pp1, f_0, f_1, el_i, o);
+            const real_t Jh = d_face_Jh(i, f);
+            const real_t dq = const_dq ? DQ(0,0) : DQ(ii, e);
+            V(v_idx, ii, e) = -dq*d_kappa*Jh*w_face[i];
+         }
+      }
+   });
+}
+
+template <int ORDER, int SDIM>
+void BatchedLOR_DG::Assemble2D()
+{
+   MFEM_VERIFY(SDIM == 2, "Surface meshes not currently supported for LOR-DG.")
+
+   static constexpr int pp1 = ORDER + 1;
+   static constexpr int ndof_per_el = pp1*pp1;
+   static constexpr int nnz_per_row = 5;
+   const int nel_ho = fes_ho.GetNE();
+
+   // Get element geometric factors; calling before AssembleFaceTerms, since
+   // in AssembleFaceTerms, element Jacobian determinants are used, potentially
+   // saving recomputation.
+   const auto factors = GeometricFactors::DETERMINANTS |
+                        GeometricFactors::JACOBIANS;
+   const auto *geom = fes_ho.GetMesh()->GetGeometricFactors(ir, factors);
+
+   // Sparse matrix entries
+   sparse_ij.SetSize(nnz_per_row*ndof_per_el*nel_ho);
+   sparse_ij.UseDevice(true);
+   sparse_ij = 0.0;
+   auto V = Reshape(sparse_ij.ReadWrite(), nnz_per_row, pp1, pp1, nel_ho);
+
+   AssembleFaceTerms();
+
+   // Populate Gauss-Lobatto quadrature rule of size (p+1)
+   IntegrationRule ir_pp1;
+   QuadratureFunctions1D::GaussLobatto(pp1, &ir_pp1);
+   Vector glx_pp1(pp1), glw_pp1(pp1);
+   for (int i = 0; i < pp1; ++i)
+   {
+      glx_pp1[i] = ir_pp1[i].x;
+      glw_pp1[i] = ir_pp1[i].weight;
+   }
+   const auto *x_pp1 = glx_pp1.Read();
+   const auto *w_1d = glw_pp1.Read();
+
+   // Get coefficients for mass and diffusion
+   const bool const_mq = c1.Size() == 1;
+   const auto MQ = const_mq
+                   ? Reshape(c1.Read(), 1, 1, 1)
+                   : Reshape(c1.Read(), pp1, pp1, nel_ho);
+   const bool const_dq = c2.Size() == 1;
+   const auto DQ = const_dq
+                   ? Reshape(c2.Read(), 1, 1, 1)
+                   : Reshape(c2.Read(), pp1, pp1, nel_ho);
+
+   const auto detJ = Reshape(geom->detJ.Read(), pp1, pp1, nel_ho);
+   const auto J = Reshape(geom->J.Read(), pp1, pp1, 2, 2, nel_ho);
+   const auto W = Reshape(ir.GetWeights().Read(), pp1, pp1);
+
+   mfem::forall(nel_ho, [=] MFEM_HOST_DEVICE (int iel_ho)
+   {
+      for (int iy = 0; iy < pp1; ++iy)
+      {
+         for (int ix = 0; ix < pp1; ++ix)
+         {
+            const real_t mq = const_mq ? MQ(0,0,0) : MQ(ix, iy, iel_ho);
+            const real_t dq = const_dq ? DQ(0,0,0) : DQ(ix, iy, iel_ho);
+
+            for (int n_idx = 0; n_idx < 2; ++n_idx)
+            {
+               for (int e_i = 0; e_i < 2; ++e_i)
+               {
+                  const int i_0 = (n_idx == 0) ? ix + e_i : ix;
+                  const int j_0 = (n_idx == 1) ? iy + e_i : iy;
+
+                  const bool bdr = (n_idx == 0 && (i_0 == 0 || i_0 == pp1)) ||
+                                   (n_idx == 1 && (j_0 == 0 || j_0 == pp1));
+
+                  if (bdr) { continue; }
+
+                  static constexpr int lex_map[] = {4, 2, 1, 3};
+                  const int v_idx_lex = e_i + n_idx*2;
+                  const int v_idx = lex_map[v_idx_lex];
+
+                  const int w_idx = (n_idx == 0) ? iy : ix;
+                  const int x_idx = (n_idx == 0) ? i_0 : j_0;
+
+                  const real_t J1 = J(ix, iy, n_idx, !n_idx, iel_ho);
+                  const real_t J2 = J(ix, iy, !n_idx, !n_idx, iel_ho);
+                  const real_t Jh = (J1*J1 + J2*J2) / detJ(ix, iy, iel_ho);
+
+                  V(v_idx, ix, iy, iel_ho) =
+                     -dq * Jh * w_1d[w_idx] / (x_pp1[x_idx] - x_pp1[x_idx -1]);
+               }
+            }
+            V(0, ix, iy, iel_ho) = mq * detJ(ix, iy, iel_ho) * W(ix, iy);
+            for (int i = 1; i < nnz_per_row; ++i)
+            {
+               V(0, ix, iy, iel_ho) -= V(i, ix, iy, iel_ho);
+            }
+         }
+      }
+   });
+}
+
+template <int ORDER>
+void BatchedLOR_DG::Assemble3D()
+{
+   static constexpr int pp1 = ORDER + 1;
+   static constexpr int ndof_per_el = pp1*pp1*pp1;
+   static constexpr int nnz_per_row = 7;
+   const int nel_ho = fes_ho.GetNE();
+
+   // Get element geometric factors; calling before AssembleFaceTerms, since
+   // in AssembleFaceTerms, element Jacobian determinants are used, potentially
+   // saving recomputation.
+   const auto factors = GeometricFactors::DETERMINANTS |
+                        GeometricFactors::JACOBIANS;
+   const auto geom = fes_ho.GetMesh()->GetGeometricFactors(ir, factors);
+
+   sparse_ij.SetSize(nnz_per_row*ndof_per_el*nel_ho);
+   sparse_ij.UseDevice(true);
+   sparse_ij = 0.0;
+   auto V = Reshape(sparse_ij.Write(), nnz_per_row, pp1, pp1, pp1, nel_ho);
+
+   AssembleFaceTerms();
+
+   // Populate Gauss-Lobatto quadrature rule of size (p+1)
+   IntegrationRule ir_pp1;
+   QuadratureFunctions1D::GaussLobatto(pp1, &ir_pp1);
+   Vector glx_pp1(pp1), glw_pp1(pp1);
+   for (int i = 0; i < pp1; ++i)
+   {
+      glx_pp1[i] = ir_pp1[i].x;
+      glw_pp1[i] = ir_pp1[i].weight;
+   }
+   const auto *x_pp1 = glx_pp1.Read();
+   const auto *w_1d = glw_pp1.Read();
+
+   const bool const_mq = c1.Size() == 1;
+   const auto MQ = const_mq
+                   ? Reshape(c1.Read(), 1, 1, 1, 1)
+                   : Reshape(c1.Read(), pp1, pp1, pp1, nel_ho);
+   const bool const_dq = c2.Size() == 1;
+   const auto DQ = const_dq
+                   ? Reshape(c2.Read(), 1, 1, 1, 1)
+                   : Reshape(c2.Read(), pp1, pp1, pp1, nel_ho);
+   const auto W = Reshape(ir.GetWeights().Read(), pp1, pp1, pp1);
+
+   const auto detJ = Reshape(geom->detJ.Read(), pp1, pp1, pp1, nel_ho);
+   const auto J = Reshape(geom->J.Read(), pp1, pp1, pp1, 3, 3, nel_ho);
+
+   mfem::forall(nel_ho, [=] MFEM_HOST_DEVICE (int iel_ho)
+   {
+      for (int iz = 0; iz < pp1; ++iz)
+      {
+         for (int iy = 0; iy < pp1; ++iy)
+         {
+            for (int ix = 0; ix < pp1; ++ix)
+            {
+               const real_t mq = const_mq ? MQ(0,0,0,0) : MQ(ix, iy, iz, iel_ho);
+               const real_t dq = const_dq ? DQ(0,0,0,0) : DQ(ix, iy, iz, iel_ho);
+
+               const real_t DETJ = detJ(ix, iy, iz, iel_ho);
+
+               for (int n_idx = 0; n_idx < 3; ++n_idx)
+               {
+                  for (int e_i = 0; e_i < 2; ++e_i)
+                  {
+                     static constexpr int lex_map[] = {5,3,2,4,1,6};
+                     const int v_idx_lex = e_i + n_idx*2;
+                     const int v_idx = lex_map[v_idx_lex];
+
+                     const int i_0 = (n_idx == 0) ? ix + e_i : ix;
+                     const int j_0 = (n_idx == 1) ? iy + e_i : iy;
+                     const int k_0 = (n_idx == 2) ? iz + e_i : iz;
+
+                     const bool bdr =
+                        (n_idx == 0 && (i_0 == 0 || i_0 == pp1)) ||
+                        (n_idx == 1 && (j_0 == 0 || j_0 == pp1)) ||
+                        (n_idx == 2 && (k_0 == 0 || k_0 == pp1));
+
+                     if (bdr) { continue; }
+
+                     int x_idx = (n_idx == 0) ? i_0 : (n_idx == 1) ? j_0 : k_0;
+                     int w_idx_1 = (n_idx == 0) ? iy : (n_idx == 1) ? iz : ix;
+                     int w_idx_2 = (n_idx == 0) ? iz : (n_idx == 1) ? ix : iy;
+
+                     const real_t J00 = J(ix, iy, iz, 0, 0, iel_ho);
+                     const real_t J01 = J(ix, iy, iz, 0, 1, iel_ho);
+                     const real_t J02 = J(ix, iy, iz, 0, 2, iel_ho);
+                     const real_t J10 = J(ix, iy, iz, 1, 0, iel_ho);
+                     const real_t J11 = J(ix, iy, iz, 1, 1, iel_ho);
+                     const real_t J12 = J(ix, iy, iz, 1, 2, iel_ho);
+                     const real_t J20 = J(ix, iy, iz, 2, 0, iel_ho);
+                     const real_t J21 = J(ix, iy, iz, 2, 1, iel_ho);
+                     const real_t J22 = J(ix, iy, iz, 2, 2, iel_ho);
+
+                     real_t JinvJinvT_diag = 0.0;
+                     if (n_idx == 0)
+                     {
+                        JinvJinvT_diag = J02*J02*(J11*J11 + J21*J21) + (J12*J21 - J11*J22)*
+                                         (J12*J21 - J11*J22) - 2*J01*J02*(J11*J12 + J21*J22) + J01*J01*
+                                         (J12*J12 + J22*J22);
+                     }
+                     else if (n_idx == 1)
+                     {
+                        JinvJinvT_diag = J02*J02*(J10*J10 + J20*J20) + (J12*J20 - J10*J22)*
+                                         (J12*J20 - J10*J22) - 2*J00*J02*(J10*J12 + J20*J22) + J00*J00*
+                                         (J12*J12 + J22*J22);
+                     }
+                     else if (n_idx == 2)
+                     {
+                        JinvJinvT_diag = J01*J01*(J10*J10 + J20*J20) + (J11*J20 - J10*J21)*
+                                         (J11*J20 - J10*J21) - 2*J00*J01*(J10*J11 + J20*J21) + J00*J00*
+                                         (J11*J11 + J21*J21);
+                     }
+
+                     const real_t Jh = JinvJinvT_diag / DETJ;
+
+                     V(v_idx, ix, iy, iz, iel_ho) = -dq * Jh * w_1d[w_idx_1] * w_1d[w_idx_2] /
+                                                    (x_pp1[x_idx] - x_pp1[x_idx -1]);
+                  }
+               }
+               V(0, ix, iy, iz, iel_ho) = mq * DETJ * W(ix, iy, iz);
+               for (int i = 1; i < 7; ++i)
+               {
+                  V(0, ix, iy, iz, iel_ho) -= V(i, ix, iy, iz, iel_ho);
+               }
+            }
+         }
+      }
+
+   });
+}
+
+} // namespace mfem

fem/pfespace.hpp

@@ -478,6 +478,7 @@ public:
        that the number of DOFs is @a ndofs. */
    const FiniteElement *GetFaceNbrFE(int i, int ndofs = 0) const;
    const FiniteElement *GetFaceNbrFaceFE(int i) const;
+   const Array<HYPRE_BigInt> &GetFaceNbrGlobalDofMapArray() { return face_nbr_glob_dof_map; }
    const HYPRE_BigInt *GetFaceNbrGlobalDofMap() { return face_nbr_glob_dof_map; }
    ElementTransformation *GetFaceNbrElementTransformation(int i) const
    { return pmesh->GetFaceNbrElementTransformation(i); }

miniapps/solvers/lor_solvers.cpp

@@ -110,7 +110,7 @@ int main(int argc, char *argv[])
    else if (string(fe) == "l") { L2 = true; }
    else { MFEM_ABORT("Bad FE type. Must be 'h', 'n', 'r', or 'l'."); }
 
-   real_t kappa = (order+1)*(order+1); // Penalty used for DG discretizations
+   real_t kappa = 10*(order+1)*(order+1); // Penalty used for DG discretizations
 
    Mesh mesh(mesh_file, 1, 1);
    const int dim = mesh.Dimension();
@@ -157,7 +157,7 @@ int main(int argc, char *argv[])
    }
    // Partial assembly not currently supported for DG or for surface meshes with
    // vector finite elements (ND or RT).
-   if (!L2 && (H1 || sdim == dim)) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
+   if (H1 || sdim == dim) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
    a.Assemble();
 
    LinearForm b(&fes);
@@ -168,6 +168,7 @@ int main(int argc, char *argv[])
       // DG boundary conditions are enforced weakly with this integrator.
       b.AddBdrFaceIntegrator(new DGDirichletLFIntegrator(u_coeff, -1.0, kappa));
    }
+   b.UseFastAssembly(true);
    b.Assemble();
 
    GridFunction x(&fes);

miniapps/solvers/plor_solvers.cpp

@@ -113,7 +113,7 @@ int main(int argc, char *argv[])
    else if (string(fe) == "l") { L2 = true; }
    else { MFEM_ABORT("Bad FE type. Must be 'h', 'n', 'r', or 'l'."); }
 
-   real_t kappa = (order+1)*(order+1); // Penalty used for DG discretizations
+   real_t kappa = 10*(order+1)*(order+1); // Penalty used for DG discretizations
 
    Mesh serial_mesh(mesh_file, 1, 1);
    const int dim = serial_mesh.Dimension();
@@ -167,7 +167,7 @@ int main(int argc, char *argv[])
    }
    // Partial assembly not currently supported for DG or for surface meshes with
    // vector finite elements (ND or RT).
-   if (!L2 && (H1 || sdim == dim)) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
+   if (H1 || sdim == dim) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
    a.Assemble();
 
    ParLinearForm b(&fes);
@@ -178,6 +178,7 @@ int main(int argc, char *argv[])
       // DG boundary conditions are enforced weakly with this integrator.
       b.AddBdrFaceIntegrator(new DGDirichletLFIntegrator(u_coeff, -1.0, kappa));
    }
+   b.UseFastAssembly(true);
    b.Assemble();
 
    ParGridFunction x(&fes);

tests/unit/CMakeLists.txt

@@ -103,6 +103,7 @@ set(UNIT_TESTS_SRCS
   fem/test_linearform_ext.cpp
   fem/test_lor.cpp
   fem/test_lor_batched.cpp
+  fem/test_lor_dg.cpp
   fem/test_nonlinearform.cpp
   fem/test_operatorjacobismoother.cpp
   fem/test_oscillation.cpp

tests/unit/fem/make_permuted_mesh.cpp

@@ -0,0 +1,162 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "make_permuted_mesh.hpp"
+
+namespace mfem
+{
+
+Mesh Mesh2D_Orientation(int face_perm_1, int face_perm_2)
+{
+   static const int dim = 2;
+   static const int nv = 6;
+   static const int nel = 2;
+   Mesh mesh(dim, nv, nel);
+   real_t x[dim];
+   x[0] = 0.0;   x[1] = 0.0;
+   mesh.AddVertex(x);
+   x[0] = 1.0;   x[1] = 0.0;
+   mesh.AddVertex(x);
+   x[0] = 2.0;   x[1] = 0.0;
+   mesh.AddVertex(x);
+   x[0] = 0.0;   x[1] = 1.0;
+   mesh.AddVertex(x);
+   x[0] = 1.0;   x[1] = 1.0;
+   mesh.AddVertex(x);
+   x[0] = 2.0;   x[1] = 1.0;
+   mesh.AddVertex(x);
+   int el[4];
+   el[0] = 0;
+   el[1] = 1;
+   el[2] = 4;
+   el[3] = 3;
+   std::rotate(&el[0], &el[face_perm_1], &el[3] + 1);
+
+   mesh.AddQuad(el);
+
+   el[0] = 1;
+   el[1] = 2;
+   el[2] = 5;
+   el[3] = 4;
+   std::rotate(&el[0], &el[face_perm_2], &el[3] + 1);
+   mesh.AddQuad(el);
+
+   mesh.FinalizeQuadMesh(true);
+   mesh.GenerateBoundaryElements();
+   mesh.Finalize();
+   return mesh;
+}
+
+void Rotation3DVertices(int *v, int ref_face, int rot)
+{
+   std::vector<int> face_1, face_2;
+
+   switch (ref_face/2)
+   {
+      case 0:
+         face_1 = {v[0], v[1], v[2], v[3]};
+         face_2 = {v[4], v[5], v[6], v[7]};
+         break;
+      case 1:
+         face_1 = {v[1], v[5], v[6], v[2]};
+         face_2 = {v[0], v[4], v[7], v[3]};
+         break;
+      case 2:
+         face_1 = {v[4], v[5], v[1], v[0]};
+         face_2 = {v[7], v[6], v[2], v[3]};
+         break;
+   }
+   if (ref_face % 2 == 0)
+   {
+      std::reverse(face_1.begin(), face_1.end());
+      std::reverse(face_2.begin(), face_2.end());
+      std::swap(face_1, face_2);
+   }
+
+   std::rotate(face_1.begin(), face_1.begin() + rot, face_1.end());
+   std::rotate(face_2.begin(), face_2.begin() + rot, face_2.end());
+
+   for (int i=0; i<4; ++i)
+   {
+      v[i] = face_1[i];
+      v[i+4] = face_2[i];
+   }
+}
+
+Mesh Mesh3D_Orientation(int face_perm_1, int face_perm_2)
+{
+   static const int dim = 3;
+   static const int nv = 12;
+   static const int nel = 2;
+   Mesh mesh(dim, nv, nel);
+   real_t x[dim];
+   x[0] = 0.0;   x[1] = 0.0;   x[2] = 0.0;
+   mesh.AddVertex(x);
+   x[0] = 1.0;   x[1] = 0.0;   x[2] = 0.0;
+   mesh.AddVertex(x);
+   x[0] = 2.0;   x[1] = 0.0;   x[2] = 0.0;
+   mesh.AddVertex(x);
+   x[0] = 0.0;   x[1] = 1.0;   x[2] = 0.0;
+   mesh.AddVertex(x);
+   x[0] = 1.0;   x[1] = 1.0;   x[2] = 0.0;
+   mesh.AddVertex(x);
+   x[0] = 2.0;   x[1] = 1.0;   x[2] = 0.0;
+   mesh.AddVertex(x);
+   x[0] = 0.0;   x[1] = 0.0;   x[2] = 1.0;
+   mesh.AddVertex(x);
+   x[0] = 1.0;   x[1] = 0.0;   x[2] = 1.0;
+   mesh.AddVertex(x);
+   x[0] = 2.0;   x[1] = 0.0;   x[2] = 1.0;
+   mesh.AddVertex(x);
+   x[0] = 0.0;   x[1] = 1.0;   x[2] = 1.0;
+   mesh.AddVertex(x);
+   x[0] = 1.0;   x[1] = 1.0;   x[2] = 1.0;
+   mesh.AddVertex(x);
+   x[0] = 2.0;   x[1] = 1.0;   x[2] = 1.0;
+   mesh.AddVertex(x);
+
+   int el[8];
+
+   el[0] = 0;
+   el[1] = 1;
+   el[2] = 4;
+   el[3] = 3;
+   el[4] = 6;
+   el[5] = 7;
+   el[6] = 10;
+   el[7] = 9;
+   Rotation3DVertices(el, face_perm_1/4, face_perm_1%4);
+   mesh.AddHex(el);
+
+   el[0] = 1;
+   el[1] = 2;
+   el[2] = 5;
+   el[3] = 4;
+   el[4] = 7;
+   el[5] = 8;
+   el[6] = 11;
+   el[7] = 10;
+   Rotation3DVertices(el, face_perm_2/4, face_perm_2%4);
+   mesh.AddHex(el);
+
+   mesh.FinalizeHexMesh(true);
+   mesh.Finalize();
+   return mesh;
+}
+
+Mesh MeshOrientation(int dim, int o1, int o2)
+{
+   if (dim == 2) { return Mesh2D_Orientation(o1, o2); }
+   else if (dim == 3) { return Mesh3D_Orientation(o1, o2); }
+   else { MFEM_ABORT("Unsupported dimension."); }
+}
+
+}

tests/unit/fem/make_permuted_mesh.hpp

@@ -0,0 +1,25 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_MAKE_PERMUTED_MESH_HPP
+#define MFEM_MAKE_PERMUTED_MESH_HPP
+
+#include "mfem.hpp"
+#include "unit_tests.hpp"
+
+namespace mfem
+{
+
+Mesh MeshOrientation(int dim, int o1, int o2);
+
+}
+
+#endif

tests/unit/fem/test_face_permutation.cpp

@@ -11,148 +11,10 @@
 
 #include "mfem.hpp"
 #include "unit_tests.hpp"
+#include "make_permuted_mesh.hpp"
 
 using namespace mfem;
 
-Mesh *mesh_2d_orientation(int face_perm_1, int face_perm_2)
-{
-   static const int dim = 2;
-   static const int nv = 6;
-   static const int nel = 2;
-   Mesh *mesh = new Mesh(dim, nv, nel);
-   real_t x[dim];
-   x[0] = 0.0;   x[1] = 0.0;
-   mesh->AddVertex(x);
-   x[0] = 1.0;   x[1] = 0.0;
-   mesh->AddVertex(x);
-   x[0] = 2.0;   x[1] = 0.0;
-   mesh->AddVertex(x);
-   x[0] = 0.0;   x[1] = 1.0;
-   mesh->AddVertex(x);
-   x[0] = 1.0;   x[1] = 1.0;
-   mesh->AddVertex(x);
-   x[0] = 2.0;   x[1] = 1.0;
-   mesh->AddVertex(x);
-   int el[4];
-   el[0] = 0;
-   el[1] = 1;
-   el[2] = 4;
-   el[3] = 3;
-   std::rotate(&el[0], &el[face_perm_1], &el[3] + 1);
-
-   mesh->AddQuad(el);
-
-   el[0] = 1;
-   el[1] = 2;
-   el[2] = 5;
-   el[3] = 4;
-   std::rotate(&el[0], &el[face_perm_2], &el[3] + 1);
-   mesh->AddQuad(el);
-
-   mesh->FinalizeQuadMesh(true);
-   mesh->GenerateBoundaryElements();
-   mesh->Finalize();
-   return mesh;
-}
-
-void rotate_3d_vertices(int *v, int ref_face, int rot)
-{
-   std::vector<int> face_1, face_2;
-
-   switch (ref_face/2)
-   {
-      case 0:
-         face_1 = {v[0], v[1], v[2], v[3]};
-         face_2 = {v[4], v[5], v[6], v[7]};
-         break;
-      case 1:
-         face_1 = {v[1], v[5], v[6], v[2]};
-         face_2 = {v[0], v[4], v[7], v[3]};
-         break;
-      case 2:
-         face_1 = {v[4], v[5], v[1], v[0]};
-         face_2 = {v[7], v[6], v[2], v[3]};
-         break;
-   }
-   if (ref_face % 2 == 0)
-   {
-      std::reverse(face_1.begin(), face_1.end());
-      std::reverse(face_2.begin(), face_2.end());
-      std::swap(face_1, face_2);
-   }
-
-   std::rotate(face_1.begin(), face_1.begin() + rot, face_1.end());
-   std::rotate(face_2.begin(), face_2.begin() + rot, face_2.end());
-
-   for (int i=0; i<4; ++i)
-   {
-      v[i] = face_1[i];
-      v[i+4] = face_2[i];
-   }
-}
-
-Mesh *mesh_3d_orientation(int face_perm_1, int face_perm_2)
-{
-   static const int dim = 3;
-   static const int nv = 12;
-   static const int nel = 2;
-   Mesh *mesh = new Mesh(dim, nv, nel);
-   real_t x[dim];
-   x[0] = 0.0;   x[1] = 0.0;   x[2] = 0.0;
-   mesh->AddVertex(x);
-   x[0] = 1.0;   x[1] = 0.0;   x[2] = 0.0;
-   mesh->AddVertex(x);
-   x[0] = 2.0;   x[1] = 0.0;   x[2] = 0.0;
-   mesh->AddVertex(x);
-   x[0] = 0.0;   x[1] = 1.0;   x[2] = 0.0;
-   mesh->AddVertex(x);
-   x[0] = 1.0;   x[1] = 1.0;   x[2] = 0.0;
-   mesh->AddVertex(x);
-   x[0] = 2.0;   x[1] = 1.0;   x[2] = 0.0;
-   mesh->AddVertex(x);
-   x[0] = 0.0;   x[1] = 0.0;   x[2] = 1.0;
-   mesh->AddVertex(x);
-   x[0] = 1.0;   x[1] = 0.0;   x[2] = 1.0;
-   mesh->AddVertex(x);
-   x[0] = 2.0;   x[1] = 0.0;   x[2] = 1.0;
-   mesh->AddVertex(x);
-   x[0] = 0.0;   x[1] = 1.0;   x[2] = 1.0;
-   mesh->AddVertex(x);
-   x[0] = 1.0;   x[1] = 1.0;   x[2] = 1.0;
-   mesh->AddVertex(x);
-   x[0] = 3.0;   x[1] = 1.0;   x[2] = 1.0;
-   mesh->AddVertex(x);
-
-   int el[8];
-
-   el[0] = 0;
-   el[1] = 1;
-   el[2] = 4;
-   el[3] = 3;
-   el[4] = 6;
-   el[5] = 7;
-   el[6] = 10;
-   el[7] = 9;
-   rotate_3d_vertices(el, face_perm_1/4, face_perm_1%4);
-   mesh->AddHex(el);
-
-   el[0] = 1;
-   el[1] = 2;
-   el[2] = 5;
-   el[3] = 4;
-   el[4] = 7;
-   el[5] = 8;
-   el[6] = 11;
-   el[7] = 10;
-   rotate_3d_vertices(el, face_perm_2/4, face_perm_2%4);
-   mesh->AddHex(el);
-
-   mesh->FinalizeHexMesh(true);
-   mesh->GenerateBoundaryElements();
-   mesh->Finalize();
-   return mesh;
-}
-
 real_t x_fn(const Vector &xvec) { return xvec[0]; }
 real_t y_fn(const Vector &xvec) { return xvec[1]; }
 real_t z_fn(const Vector &xvec) { return xvec[2]; }
@@ -214,10 +76,9 @@ TEST_CASE("2D Face Permutation", "[Face Permutation]")
    {
       for (int fp1=0; fp1<4; ++fp1)
       {
-         Mesh *mesh = mesh_2d_orientation(fp1, fp2);
-         real_t error = TestFaceRestriction(*mesh, order);
+         Mesh mesh = MeshOrientation(2, fp1, fp2);
+         real_t error = TestFaceRestriction(mesh, order);
          max_err = std::max(max_err, error);
-         delete mesh;
       }
    }
    REQUIRE(max_err < 1e-15);
@@ -231,10 +92,9 @@ TEST_CASE("3D Face Permutation", "[Face Permutation]")
    {
       for (int fp1=0; fp1<24; ++fp1)
       {
-         Mesh *mesh = mesh_3d_orientation(fp1, fp2);
-         real_t error = TestFaceRestriction(*mesh, order);
+         Mesh mesh = MeshOrientation(3, fp1, fp2);
+         real_t error = TestFaceRestriction(mesh, order);
          max_err = std::max(max_err, error);
-         delete mesh;
       }
    }
    REQUIRE(max_err < 1e-15);

tests/unit/fem/test_lor_batched.cpp

@@ -11,6 +11,7 @@
 
 #include "mfem.hpp"
 #include "unit_tests.hpp"
+#include "../linalg/test_same_matrices.hpp"
 #include "../../fem/lor/lor_ads.hpp"
 #include "../../fem/lor/lor_ams.hpp"
 #include <memory>
@@ -25,47 +26,6 @@ using namespace mfem;
 namespace lor_batched
 {
 
-void TestSameMatrices(SparseMatrix &A1, const SparseMatrix &A2,
-                      HYPRE_BigInt *cmap1=nullptr,
-                      std::unordered_map<HYPRE_BigInt,int> *cmap2inv=nullptr)
-{
-   REQUIRE(A1.Height() == A2.Height());
-   int n = A1.Height();
-
-   const int *I1 = A1.HostReadI();
-   const int *J1 = A1.HostReadJ();
-   const real_t *V1 = A1.HostReadData();
-
-   A2.HostReadI();
-   A2.HostReadJ();
-   A2.HostReadData();
-
-   real_t error = 0.0;
-
-   for (int i=0; i<n; ++i)
-   {
-      for (int jj=I1[i]; jj<I1[i+1]; ++jj)
-      {
-         int j = J1[jj];
-         if (cmap1)
-         {
-            if (cmap2inv->count(cmap1[j]) > 0)
-            {
-               j = (*cmap2inv)[cmap1[j]];
-            }
-            else
-            {
-               error = std::max(error, std::fabs(V1[jj]));
-               continue;
-            }
-         }
-         error = std::max(error, std::fabs(V1[jj] - A2(i,j)));
-      }
-   }
-
-   REQUIRE(error == MFEM_Approx(0.0, 1e-10));
-}
-
 template <typename FE_COLL>
 FE_COLL *NewLOR_FE_Collection(int order, int dim)
 {
@@ -157,30 +117,6 @@ TEST_CASE("LOR Batched RT", "[LOR][BatchedLOR][CUDA]")
 
 #ifdef MFEM_USE_MPI
 
-void TestSameMatrices(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);
-
-   TestSameMatrices(diag1, diag2);
-
-   if (cmap1)
-   {
-      std::unordered_map<HYPRE_BigInt,int> cmap2inv;
-      for (int i=0; i<offd2.Width(); ++i) { cmap2inv[cmap2[i]] = i; }
-      TestSameMatrices(offd1, offd2, cmap1, &cmap2inv);
-   }
-   else
-   {
-      TestSameMatrices(offd1, offd2);
-   }
-}
-
 template <typename FE_COLL, typename INTEG_1, typename INTEG_2>
 void ParTestBatchedLOR()
 {

tests/unit/fem/test_lor_dg.cpp

@@ -0,0 +1,384 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "mfem.hpp"
+#include "unit_tests.hpp"
+#include "make_permuted_mesh.hpp"
+#include "../linalg/test_same_matrices.hpp"
+
+using namespace mfem;
+
+class DG_LOR_DiffusionPreconditioner : public BilinearFormIntegrator
+{
+   Mesh &mesh;
+   double kappa;
+   int p;
+   IntegrationRule gl_p, gl_pp1;
+   Vector shape1, shape2, nor, ni;
+
+public:
+   DG_LOR_DiffusionPreconditioner(Mesh &mesh_, int p_, double kappa_)
+      : mesh(mesh_), kappa(kappa_), p(p_)
+   {
+      QuadratureFunctions1D::GaussLobatto(p+1, &gl_p);
+      QuadratureFunctions1D::GaussLobatto(p+2, &gl_pp1);
+   }
+
+   double PenaltyFactor(int idx1, int idx2)
+   {
+      int pp1 = p + 1;
+
+      int x1 = idx1 % pp1;
+      int y1 = (idx1 / pp1) % pp1;
+      int z1 = (idx1 / pp1) / pp1;
+
+      int x2 = idx2 % pp1;
+      int y2 = (idx2 / pp1) % pp1;
+      int z2 = (idx2 / pp1) / pp1;
+
+      int dim = mesh.Dimension();
+
+      auto compute_factor = [&](int i1, int i2)
+      {
+         int j = std::min(i1, i2);
+         if (i1 == i2)
+         {
+            double w = gl_p[j].weight;
+            double k = gl_pp1[i1+1].x - gl_pp1[i1].x;
+            return w/k;
+         }
+         else
+         {
+            double h = gl_p[j+1].x - gl_p[j].x;
+            double k1 = gl_pp1[i1+1].x - gl_pp1[i1].x;
+            double k2 = gl_pp1[i2+1].x - gl_pp1[i2].x;
+            double avg = 0.5*k1 + 0.5*k2;
+            return avg/h;
+         }
+      };
+
+      double factor = compute_factor(x1, x2);
+      if (dim >= 2) { factor *= compute_factor(y1, y2); }
+      if (dim == 3) { factor *= compute_factor(z1, z2); }
+
+      return factor;
+   }
+
+   double BdrPenaltyFactor(int idx, int f)
+   {
+      int pp1 = p+1;
+
+      int x = idx % pp1;
+      int y = (idx / pp1) % pp1;
+      int z = (idx / pp1) / pp1;
+
+      int dim = mesh.Dimension();
+
+      auto subcell_size = [&](int i)
+      {
+         return gl_pp1[i+1].x - gl_pp1[i].x;
+      };
+
+      double factor = (p+1)*(p+1);
+      if (dim == 1)
+      {
+         factor *= subcell_size(x);
+      }
+      else if (dim == 2)
+      {
+         int ni, nj;
+         ni = (f == 1 || f == 3) ? x : y;
+         nj = (f == 1 || f == 3) ? y : x;
+         factor *= subcell_size(ni)/subcell_size(nj)*gl_p[nj].weight;
+      }
+      else if (dim == 3)
+      {
+         int ni, nj, nk;
+         if (f == 2 || f == 4) { ni = x; nj = y; nk = z; }
+         else if (f == 1 || f == 3) { ni = y; nj = x; nk = z; }
+         else { ni = z; nj = x; nk = y; }
+         factor *= subcell_size(ni)/subcell_size(nj)/subcell_size(nk);
+         factor *= gl_p[nj].weight*gl_p[nk].weight;
+      }
+
+      return factor;
+   }
+
+   using BilinearFormIntegrator::AssembleFaceMatrix;
+   virtual void AssembleFaceMatrix(const FiniteElement &el1,
+                                   const FiniteElement &el2,
+                                   FaceElementTransformations &Trans,
+                                   DenseMatrix &elmat)
+   {
+      int dim, ndof1, ndof2, ndofs;
+      double w, wq = 0.0;
+
+      dim = el1.GetDim();
+      ndof1 = el1.GetDof();
+
+      nor.SetSize(dim);
+      ni.SetSize(dim);
+
+      shape1.SetSize(ndof1);
+      if (Trans.Elem2No >= 0)
+      {
+         ndof2 = el2.GetDof();
+         shape2.SetSize(ndof2);
+      }
+      else
+      {
+         ndof2 = 0;
+      }
+
+      int face_no;
+      if (ndof2) { face_no = Trans.ElementNo; }
+      else { face_no = mesh.GetBdrElementFaceIndex(Trans.ElementNo); }
+
+      int info1, info2;
+      mesh.GetFaceInfos(face_no, &info1, &info2);
+      int local_face = info1/64;
+
+      const CoarseFineTransformations &cftr = mesh.GetRefinementTransforms();
+
+      double factor;
+      bool interior = false;
+      if (Trans.Elem2No >= 0 && Trans.Elem2No < mesh.GetNE())
+      {
+         int parent_el1 = cftr.embeddings[Trans.Elem1No].parent;
+         int parent_el2 = cftr.embeddings[Trans.Elem2No].parent;
+         if (parent_el1 == parent_el2)
+         {
+            interior = true;
+            factor = PenaltyFactor(cftr.embeddings[Trans.Elem1No].matrix,
+                                   cftr.embeddings[Trans.Elem2No].matrix);
+         }
+      }
+      if (!interior)
+      {
+         factor = kappa*BdrPenaltyFactor(cftr.embeddings[Trans.Elem1No].matrix,
+                                         local_face);
+      }
+
+      ndofs = ndof1 + ndof2;
+      elmat.SetSize(ndofs);
+      elmat = 0.0;
+
+      const IntegrationRule *ir = IntRule;
+      if (ir == NULL) { ir = &IntRules.Get(Trans.GetGeometryType(), 1); }
+
+      for (int p = 0; p < ir->GetNPoints(); p++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(p);
+         Trans.SetAllIntPoints(&ip);
+         const IntegrationPoint &eip1 = Trans.GetElement1IntPoint();
+         const IntegrationPoint &eip2 = Trans.GetElement2IntPoint();
+
+         if (dim == 1) { nor(0) = 2*eip1.x - 1.0; }
+         else { CalcOrtho(Trans.Jacobian(), nor); }
+
+         el1.CalcShape(eip1, shape1);
+         w = ip.weight;
+         ni.Set(w, nor);
+
+         double h_face = nor.Norml2();
+         double h_el = Trans.Elem1->Weight();
+         double h = h_el/h_face; // perpendicular element size
+
+         if (ndof2)
+         {
+            el2.CalcShape(eip2, shape2);
+            double h_el_2 = Trans.Elem2->Weight();
+            h = 0.5*h + 0.5*h_el_2/h_face; // average both element sizes
+         }
+
+         if (interior)
+         {
+            wq = w*factor*h_face/h;
+         }
+         else
+         {
+            wq = w*factor*h_face/h;
+         }
+         for (int i = 0; i < ndof1; i++)
+         {
+            const double wsi = wq*shape1(i);
+            for (int j = 0; j < ndof1; j++)
+            {
+               elmat(i, j) += wsi * shape1(j);
+            }
+         }
+         if (ndof2)
+         {
+            for (int i = 0; i < ndof2; i++)
+            {
+               const double wsi = wq*shape2(i);
+               for (int j = 0; j < ndof1; j++)
+               {
+                  elmat(ndof1 + i, j) -= wsi * shape1(j);
+                  elmat(j, ndof1 + i) -= wsi * shape1(j);
+               }
+               for (int j = 0; j < ndof2; j++)
+               {
+                  elmat(ndof1 + i, ndof1 + j) += wsi * shape2(j);
+               }
+            }
+         }
+      }
+   }
+};
+
+class DG_LOR_MassPreconditioner : public BilinearFormIntegrator
+{
+   Mesh &mesh_ho, &mesh_lor;
+   const int p;
+   IntegrationRule gll;
+   Coefficient *Q;
+
+public:
+   DG_LOR_MassPreconditioner(Mesh &mesh_ho_,
+                             Mesh &mesh_lor_,
+                             int p_,
+                             Coefficient *Q_)
+      : mesh_ho(mesh_ho_),
+        mesh_lor(mesh_lor_),
+        p(p_),
+        Q(Q_)
+   {
+      QuadratureFunctions1D::GaussLobatto(p+1, &gll);
+   }
+
+   void AssembleElementMatrix(const FiniteElement &el,
+                              ElementTransformation &Tr,
+                              DenseMatrix &elmat)
+   {
+      const CoarseFineTransformations &cftr = mesh_lor.GetRefinementTransforms();
+      const int parent_el = cftr.embeddings[Tr.ElementNo].parent;
+      // We use the point matrix index to identify the local LOR element index
+      // within the high-order coarse element.
+      const int lor_index = cftr.embeddings[Tr.ElementNo].matrix;
+
+      // Assuming piecewise constant
+      elmat.SetSize(1);
+
+      const int dim = mesh_ho.Dimension();
+      IntegrationPoint ip;
+      if (dim == 2)
+      {
+         const int iy = lor_index / (p + 1);
+         const int ix = lor_index  % (p + 1);
+         ip.x = gll[ix].x;
+         ip.y = gll[iy].x;
+
+         elmat(0,0) = gll[ix].weight * gll[iy].weight;
+      }
+      else if (dim == 3)
+      {
+         const int iz = lor_index / (p + 1) / (p + 1);
+         const int iy = (lor_index / (p + 1)) % (p + 1);
+         const int ix = lor_index  % (p + 1);
+
+         ip.x = gll[ix].x;
+         ip.y = gll[iy].x;
+         ip.z = gll[iz].x;
+
+         elmat(0,0) = gll[ix].weight * gll[iy].weight * gll[iz].weight;
+      }
+
+      ElementTransformation &Tr_ho = *mesh_ho.GetElementTransformation(parent_el);
+      Tr_ho.SetIntPoint(&ip);
+      const real_t detJ = Tr_ho.Weight();
+      elmat(0,0) *= detJ;
+
+      if (Q)
+      {
+         elmat(0,0) *= Q->Eval(Tr_ho, ip);
+      }
+   }
+};
+
+static void TestBatchedLOR_DG(Mesh &mesh, int order)
+{
+   DG_FECollection fec(order, mesh.Dimension(), BasisType::GaussLobatto);
+   FiniteElementSpace fespace(&mesh, &fec);
+
+   // Set up some coefficients using H1 grid functions
+   H1_FECollection h1fec(2, mesh.Dimension());
+   FiniteElementSpace h1fes(&mesh, &h1fec);
+   GridFunction gf1(&h1fes), gf2(&h1fes);
+   gf1.Randomize(1);
+   gf2.Randomize(2);
+   GridFunctionCoefficient mass_coeff(&gf1);
+   GridFunctionCoefficient diff_coeff(&gf2);
+
+   ConstantCoefficient one(1.0);
+   constexpr int sigma = -1.0;
+   const int eta = 2;
+   const int kappa = eta * (order + 1) * (order + 1);
+   BilinearForm a(&fespace);
+   a.AddDomainIntegrator(new DiffusionIntegrator);
+   a.AddInteriorFaceIntegrator(new DGDiffusionIntegrator(one, sigma, kappa));
+   a.AddBdrFaceIntegrator(new DGDiffusionIntegrator(one, sigma, kappa));
+
+   Array<int> ess_dofs; // Empty
+   LORDiscretization lor(fespace);
+   lor.AssembleSystem(a, ess_dofs);
+   SparseMatrix &A1 = lor.GetAssembledMatrix();
+
+   FiniteElementSpace &fes_lor = lor.GetFESpace();
+   Mesh &mesh_lor = *fes_lor.GetMesh();
+   BilinearForm a_lor(&fes_lor);
+   a_lor.AddBdrFaceIntegrator(new DG_LOR_DiffusionPreconditioner(
+                                 mesh_lor, order, eta));
+   a_lor.AddInteriorFaceIntegrator(new DG_LOR_DiffusionPreconditioner(
+                                      mesh_lor, order, eta));
+
+   a_lor.Assemble();
+   a_lor.Finalize();
+   SparseMatrix &A2 = a_lor.SpMat();
+
+   TestSameMatrices(A1, A2);
+}
+
+TEST_CASE("LOR Batched DG Orientation", "[LOR][BatchedLOR][CUDA]")
+{
+   const int order = 3;
+   const int dim = launch_all_non_regression_tests ? GENERATE(2, 3) : 2;
+   const int orientation1 = GENERATE_COPY(range(0, dim == 2 ? 4 : 24));
+   const int orientation2 = GENERATE_COPY(range(0, dim == 2 ? 4 : 24));
+
+   CAPTURE(order, dim, orientation1, orientation2);
+
+   Mesh mesh = MeshOrientation(dim, orientation1, orientation2);
+   TestBatchedLOR_DG(mesh, order);
+}
+
+TEST_CASE("LOR Batched DG", "[LOR][BatchedLOR][CUDA]")
+{
+   const int order = 3;
+   const auto mesh_fname = GENERATE(
+                              "../../data/beam-quad.mesh",
+                              "../../data/l-shape.mesh",
+                              "../../data/beam-hex.mesh",
+                              "../../data/fichera.mesh"
+                           );
+   CAPTURE(mesh_fname);
+   Mesh mesh = Mesh::LoadFromFile(mesh_fname);
+
+   mesh.Transform([](const Vector &xin, Vector &xout)
+   {
+      for (int d = 0; d < xin.Size(); ++d)
+      {
+         xout[d] = xin[d] * (1.0 + d / 3.0);
+      }
+   });
+
+   TestBatchedLOR_DG(mesh, order);
+}

tests/unit/linalg/test_same_matrices.hpp

@@ -0,0 +1,96 @@
+// Copyright (c) 2010-2025, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_TEST_SAME_MATRICES_HPP
+#define MFEM_TEST_SAME_MATRICES_HPP
+
+#include "mfem.hpp"
+
+namespace mfem
+{
+
+#ifndef MFEM_USE_MPI
+#define HYPRE_BigInt int
+#endif // MFEM_USE_MPI
+
+inline void TestSameMatrices(
+   SparseMatrix &A1, const SparseMatrix &A2,
+   HYPRE_BigInt *cmap1=nullptr,
+   std::unordered_map<HYPRE_BigInt,int> *cmap2inv=nullptr)
+{
+   REQUIRE(A1.Height() == A2.Height());
+   int n = A1.Height();
+
+   const int *I1 = A1.HostReadI();
+   const int *J1 = A1.HostReadJ();
+   const real_t *V1 = A1.HostReadData();
+
+   A2.HostReadI();
+   A2.HostReadJ();
+   A2.HostReadData();
+
+   real_t error = 0.0;
+
+   for (int i=0; i<n; ++i)
+   {
+      for (int jj=I1[i]; jj<I1[i+1]; ++jj)
+      {
+         int j = J1[jj];
+         if (cmap1)
+         {
+            if (cmap2inv->count(cmap1[j]) > 0)
+            {
+               j = (*cmap2inv)[cmap1[j]];
+            }
+            else
+            {
+               error = std::max(error, std::fabs(V1[jj]));
+               continue;
+            }
+         }
+         error = std::max(error, std::fabs(V1[jj] - A2(i,j)));
+      }
+   }
+
+   REQUIRE(error == MFEM_Approx(0.0, 1e-10));
+}
+
+#ifdef MFEM_USE_MPI
+
+inline void TestSameMatrices(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);
+
+   TestSameMatrices(diag1, diag2);
+
+   if (cmap1)
+   {
+      std::unordered_map<HYPRE_BigInt,int> cmap2inv;
+      for (int i=0; i<offd2.Width(); ++i) { cmap2inv[cmap2[i]] = i; }
+      TestSameMatrices(offd1, offd2, cmap1, &cmap2inv);
+   }
+   else
+   {
+      TestSameMatrices(offd1, offd2);
+   }
+}
+
+#endif // MFEM_USE_MPI
+
+} // namespace mfem
+
+#endif