Updates to year on copyright statement on makefiles

Currently just a cleaned-up version of ex9. Will add diffusion and
IMEX time integration.

examples/CMakeLists.txt

@@ -46,6 +46,7 @@ list(APPEND ALL_EXE_SRCS
   ex38.cpp
   ex39.cpp
   ex40.cpp
+  ex41.cpp
   )
 
 if (MFEM_USE_MPI)
@@ -89,6 +90,7 @@ if (MFEM_USE_MPI)
     ex37p.cpp
     ex39p.cpp
     ex40p.cpp
+    ex41p.cpp
   )
 endif()
 
@@ -151,7 +153,7 @@ if (MFEM_ENABLE_TESTING)
       ex1 ex3 ex4 ex5 ex6 ex9 ex14 ex22 ex24 ex25 ex26 ex34
       # parallel examples with device support:
       ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex9p ex13p ex14p ex22p ex24p ex25p
-      ex26p ex34p ex35p)
+      ex26p ex34p ex35p ex41p)
   set(MFEM_TEST_DEVICE)
   if (MFEM_USE_CUDA)
     set(MFEM_TEST_DEVICE "cuda")

examples/ex41.cpp

@@ -0,0 +1,455 @@
+//                                MFEM Example 41
+//
+// Compile with: make ex41
+//
+// Sample runs: ex41
+//              ex41 -p 0 -r 2 -dt 0.01 -tf 10
+//              ex41 -m ../data/periodic-hexagon.mesh -p 0 -r 2 -dt 0.005 -tf 10
+//              ex41 -m ../data/periodic-square.mesh -p 1 -r 2 -dt 0.005 -tf 9
+//              ex41 -m ../data/periodic-hexagon.mesh -p 1 -r 2 -dt 0.005 -tf 9
+//              ex41 -m ../data/amr-quad.mesh -p 1 -r 2 -dt 0.002 -tf 9
+//              ex41 -m ../data/star-q3.mesh -p 1 -r 2 -dt 0.001 -tf 9
+//              ex41 -m ../data/star-mixed.mesh -p 1 -r 2 -dt 0.005 -tf 9
+//              ex41 -m ../data/disc-nurbs.mesh -p 1 -r 3 -dt 0.005 -tf 9
+//              ex41 -m ../data/disc-nurbs.mesh -p 2 -r 3 -dt 0.005 -tf 9
+//              ex41 -m ../data/periodic-square.mesh -p 3 -r 4 -dt 0.0025 -tf 9 -vs 20
+//              ex41 -m ../data/periodic-cube.mesh -p 0 -r 2 -o 2 -dt 0.01 -tf 8
+//
+// Device sample runs:
+//
+// Description:  This example code solves the time-dependent advection-diffusion
+//               equation du/dt + v.grad(u) - a div(grad(u)) = 0, where v is a
+//               given fluid velocity, a is the diffusion coefficient, and
+//               u0(x)=u(0,x) is a given initial condition.
+//
+//               The example demonstrates the use of Discontinuous Galerkin (DG)
+//               bilinear forms in MFEM (face integrators), and the use of IMEX ODE time integrators.
+
+#include "mfem.hpp"
+
+using namespace std;
+using namespace mfem;
+
+// Choice for the problem setup. The fluid velocity, initial condition and
+// inflow boundary condition are chosen based on this parameter.
+int problem;
+
+// Velocity coefficient
+void velocity_function(const Vector &x, Vector &v);
+
+// Initial condition
+real_t u0_function(const Vector &x);
+
+// Inflow boundary condition - 0 for this example
+real_t inflow_function(const Vector &x);
+
+// Mesh bounding box
+Vector bb_min, bb_max;
+
+class DG_Solver : public Solver
+{
+private:
+   SparseMatrix &M, &K, &S, A;
+   CGSolver linear_solver;
+   BlockILU prec;
+   real_t dt;
+public:
+   DG_Solver(SparseMatrix &M_, SparseMatrix &K_, SparseMatrix &S_,
+             const FiniteElementSpace &fes)
+      : M(M_),
+        K(K_),
+        S(S_),
+        prec(fes.GetTypicalFE()->GetDof(),
+             BlockILU::Reordering::MINIMUM_DISCARDED_FILL),
+        dt(1.0)
+   {
+      linear_solver.iterative_mode = false;
+      linear_solver.SetRelTol(1e-9);
+      linear_solver.SetAbsTol(0.0);
+      linear_solver.SetMaxIter(100);
+      linear_solver.SetPrintLevel(0);
+      linear_solver.SetPreconditioner(prec);
+   }
+
+   void SetTimeStep(real_t dt_)
+   {
+      if (dt_ != dt)
+      {
+         dt = dt_;
+         // Form operator A = M + dt*S
+         A = S;
+         A *= dt;
+         A += M;
+
+         // this will also call SetOperator on the preconditioner
+         linear_solver.SetOperator(A);
+      }
+   }
+
+   void SetOperator(const Operator &op) override
+   {
+      linear_solver.SetOperator(op);
+   }
+
+   void Mult(const Vector &x, Vector &y) const override
+   {
+      linear_solver.Mult(x, y);
+   }
+};
+
+/** A time-dependent operator for the right-hand side of the ODE. The DG weak
+    form of the advection-diffusion equation is (M + dt S) du/dt = Su - K u + b, where M and K are the mass
+    and advection matrices, and b describes the flow on the boundary. In the case of IMEX evolution, the diffusion term is treated
+    implicitly, and the advection term is treated explicitly.  */
+class IMEX_Evolution : public SplitTimeDependentOperator
+{
+private:
+   BilinearForm &M, &K, &S;
+   const Vector &b;
+   unique_ptr<Solver> M_prec;
+   CGSolver M_solver;
+   unique_ptr<DG_Solver> dg_solver;
+
+   mutable Vector z;
+
+public:
+   IMEX_Evolution(BilinearForm &M_, BilinearForm &K_, BilinearForm &S_,
+                  const Vector &b_);
+
+   void Mult1(const Vector &x, Vector &y) const;
+   void ImplicitSolve2(const real_t dt, const Vector &x, Vector &k) override;
+};
+
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   problem = 0;
+   const char *mesh_file = "../data/periodic-square.mesh";
+   int ref_levels = 2;
+   int order = 3;
+   const char *device_config = "cpu";
+   int ode_solver_type = 58;
+   real_t t_final = 10.0;
+   real_t dt = 0.001;
+   bool paraview = false;
+   int vis_steps = 50;
+   real_t diffusion_term = 0.01;
+   real_t kappa = (order+1)*(order+1);
+   real_t sigma = -1.0;
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
+   args.AddOption(&problem, "-p", "--problem",
+                  "Problem setup to use. See options in velocity_function().");
+   args.AddOption(&ref_levels, "-r", "--refine",
+                  "Number of times to refine the mesh uniformly.");
+   args.AddOption(&order, "-o", "--order", "Order of the finite elements.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&ode_solver_type, "-s", "--ode-solver",
+                  ODESolver::Types.c_str());
+   args.AddOption(&t_final, "-tf", "--t-final", "Final time; start time is 0.");
+   args.AddOption(&dt, "-dt", "--time-step", "Time step.");
+   args.AddOption(&paraview, "-paraview", "--paraview-datafiles", "-no-paraview",
+                  "--no-paraview-datafiles",
+                  "Save data files for ParaView (paraview.org) visualization.");
+   args.AddOption(&vis_steps, "-vs", "--visualization-steps",
+                  "Visualize every n-th timestep.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   if (kappa < 0)
+   {
+      kappa = (order+1)*(order+1);
+   }
+   args.PrintOptions(cout);
+   Device device(device_config);
+   device.Print();
+
+   // 2. Read the mesh from the given mesh file. We can handle geometrically
+   //    periodic meshes in this code.
+   Mesh mesh(mesh_file);
+   const int dim = mesh.Dimension();
+
+   // 3. Define the Split ODE solver used for time integration. The IMEX solvers currently
+   // available are: 55 - Forward Backward Euler, 56 - IMEXRK2(2,2,2), 57 - IMEXRK2(2,3,2), and
+   // 58 - IMEX_DIRK_RK3.
+   unique_ptr<SplitODESolver> ode_solver = SplitODESolver::Select(ode_solver_type);
+
+   // 4. Refine the mesh to increase the resolution. In this example we do
+   //    'ref_levels' of uniform refinement, where 'ref_levels' is a
+   //    command-line parameter.
+   for (int lev = 0; lev < ref_levels; lev++) {mesh.UniformRefinement();}
+   if (mesh.NURBSext) {mesh.SetCurvature(max(order, 1));}
+   mesh.GetBoundingBox(bb_min, bb_max, max(order, 1));
+
+   // 5. Define the discontinuous DG finite element space of the given
+   //    polynomial order on the refined mesh.
+   DG_FECollection fec(order, dim, BasisType::GaussLobatto);
+   FiniteElementSpace fes(&mesh, &fec);
+
+   cout << "Number of unknowns: " << fes.GetVSize() << endl;
+
+   // 6. Set up and assemble the bilinear and linear forms corresponding to the
+   //    DG discretization. The DGTraceIntegrator involves integrals over mesh
+   //    interior faces.
+   VectorFunctionCoefficient velocity(dim, velocity_function);
+   FunctionCoefficient inflow(inflow_function); //Zero for now
+   ConstantCoefficient diff_coeff(diffusion_term);
+
+   BilinearForm m(&fes);
+   BilinearForm k(&fes);
+   BilinearForm s(&fes);
+
+   m.AddDomainIntegrator(new MassIntegrator);
+
+   constexpr real_t alpha = -1.0;
+   k.AddDomainIntegrator(new ConvectionIntegrator(velocity, alpha));
+   k.AddInteriorFaceIntegrator(new NonconservativeDGTraceIntegrator(velocity,
+                                                                    alpha));
+   k.AddBdrFaceIntegrator(new NonconservativeDGTraceIntegrator(velocity, alpha));
+
+   s.AddDomainIntegrator(new DiffusionIntegrator(diff_coeff));
+   s.AddInteriorFaceIntegrator(new DGDiffusionIntegrator(diff_coeff, sigma,
+                                                         kappa));
+   s.AddBdrFaceIntegrator(new DGDiffusionIntegrator(diff_coeff, sigma, kappa));
+
+   LinearForm b(&fes);
+   b.AddBdrFaceIntegrator(new BoundaryFlowIntegrator(inflow, velocity, alpha));
+
+   int skip_zeros = 0;
+   m.Assemble(skip_zeros);
+   k.Assemble(skip_zeros);
+   s.Assemble(skip_zeros);
+   b.Assemble();
+
+   m.Finalize(skip_zeros);
+   k.Finalize(skip_zeros);
+   s.Finalize(skip_zeros);
+
+   // 7. Define the initial conditions.
+   FunctionCoefficient u0(u0_function);
+   GridFunction u(&fes);
+   u.ProjectCoefficient(u0);
+
+   // 8. Set up paraview visualization, if desired.
+   unique_ptr<ParaViewDataCollection> pv;
+   if (paraview)
+   {
+      pv = make_unique<ParaViewDataCollection>("Example41", &mesh);
+      pv->SetPrefixPath("ParaView");
+      pv->RegisterField("solution", &u);
+      pv->SetLevelsOfDetail(order);
+      pv->SetDataFormat(VTKFormat::BINARY);
+      pv->SetHighOrderOutput(true);
+      pv->SetCycle(0);
+      pv->SetTime(0.0);
+      pv->Save();
+   }
+
+   // 9. Define the time-dependent evolution operator describing the ODE
+   //    right-hand side, and perform time-integration (looping over the time
+   //    iterations, ti, with a time-step dt).
+   IMEX_Evolution adv(m, k, s, b);
+
+   real_t t = 0.0;
+   adv.SetTime(t);
+   ode_solver->Init(adv);
+
+   bool done = false;
+   for (int ti = 0; !done; )
+   {
+      real_t dt_real = min(dt, t_final - t);
+      ode_solver->Step(u, t, dt_real);
+      ti++;
+
+      done = (t >= t_final - 1e-8*dt);
+
+      if (done || ti % vis_steps == 0)
+      {
+         cout << "time step: " << ti << ", time: " << t << endl;
+         if (paraview)
+         {
+            pv->SetCycle(ti);
+            pv->SetTime(t);
+            pv->Save();
+         }
+      }
+   }
+   return 0;
+}
+
+
+// Implementation of class IMEX_Evolution
+IMEX_Evolution::IMEX_Evolution(BilinearForm &M_, BilinearForm &K_,
+                               BilinearForm &S_, const Vector &b_)
+   : SplitTimeDependentOperator(M_.FESpace()->GetTrueVSize()),
+     M(M_), K(K_), S(S_), b(b_), z(height)
+{
+   Array<int> ess_tdof_list;
+   if (M.GetAssemblyLevel() == AssemblyLevel::LEGACY)
+   {
+      M_prec = make_unique<DSmoother>(M.SpMat());
+      M_solver.SetOperator(M.SpMat());
+      dg_solver = make_unique<DG_Solver>(M.SpMat(), K.SpMat(), S.SpMat(),
+                                         *M.FESpace());
+   }
+   else
+   {
+      M_prec = make_unique<OperatorJacobiSmoother>(M, ess_tdof_list);
+      M_solver.SetOperator(M);
+      dg_solver = NULL;
+   }
+   M_solver.SetPreconditioner(*M_prec);
+   M_solver.iterative_mode = false;
+   M_solver.SetRelTol(1e-9);
+   M_solver.SetAbsTol(0.0);
+   M_solver.SetMaxIter(100);
+   M_solver.SetPrintLevel(0);
+}
+
+void IMEX_Evolution::Mult1(const Vector &x, Vector &y) const
+{
+   // Perform the explicit step
+   // y = M^{-1} (K x + b)
+   K.Mult(x, z);
+   z += b;
+   M_solver.Mult(z, y);
+}
+
+void IMEX_Evolution::ImplicitSolve2(const real_t dt, const Vector &x, Vector &k)
+{
+   // Perform the implicit step
+   // solve for k, k = -(M+dt S)^{-1} S x
+   MFEM_VERIFY(dg_solver != NULL,
+               "Implicit time integration is not supported with partial assembly");
+   S.Mult(x, z);
+   z*= -1.0;
+   dg_solver->SetTimeStep(dt);
+   dg_solver->Mult(z, k);
+}
+
+// Velocity coefficient
+void velocity_function(const Vector &x, Vector &v)
+{
+   int dim = x.Size();
+
+   // map to the reference [-1,1] domain
+   Vector X(dim);
+   for (int i = 0; i < dim; i++)
+   {
+      real_t center = (bb_min[i] + bb_max[i]) * 0.5;
+      X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
+   }
+   switch (problem)
+   {
+      case 0:
+      {
+         // Translations in 1D, 2D, and 3D
+         switch (dim)
+         {
+            case 1: v(0) = 1.0; break;
+            case 2: v(0) = sqrt(2./3.); v(1) = sqrt(1./3.); break;
+            case 3: v(0) = sqrt(3./6.); v(1) = sqrt(2./6.); v(2) = sqrt(1./6.);
+               break;
+         }
+         break;
+      }
+      case 1:
+      case 2:
+      {
+         // Clockwise rotation in 2D around the origin
+         const real_t w = M_PI/2;
+         switch (dim)
+         {
+            case 1: v(0) = 1.0; break;
+            case 2: v(0) = w*X(1); v(1) = -w*X(0); break;
+            case 3: v(0) = w*X(1); v(1) = -w*X(0); v(2) = 0.0; break;
+         }
+         break;
+      }
+      case 3:
+      {
+         // Clockwise twisting rotation in 2D around the origin
+         const real_t w = M_PI/2;
+         real_t d = max((X(0)+1.)*(1.-X(0)),0.) * max((X(1)+1.)*(1.-X(1)),0.);
+         d = d*d;
+         switch (dim)
+         {
+            case 1: v(0) = 1.0; break;
+            case 2: v(0) = d*w*X(1); v(1) = -d*w*X(0); break;
+            case 3: v(0) = d*w*X(1); v(1) = -d*w*X(0); v(2) = 0.0; break;
+         }
+         break;
+      }
+   }
+}
+
+// Initial condition
+real_t u0_function(const Vector &x)
+{
+   int dim = x.Size();
+
+   // map to the reference [-1,1] domain
+   Vector X(dim);
+   for (int i = 0; i < dim; i++)
+   {
+      real_t center = (bb_min[i] + bb_max[i]) * 0.5;
+      X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
+   }
+
+   switch (problem)
+   {
+      case 0:
+      case 1:
+      {
+         switch (dim)
+         {
+            case 1:
+               return exp(-40.*pow(X(0)-0.5,2));
+            case 2:
+            case 3:
+            {
+               real_t rx = 0.45, ry = 0.25, cx = 0., cy = -0.2, w = 10.;
+               if (dim == 3)
+               {
+                  const real_t s = (1. + 0.25*cos(2*M_PI*X(2)));
+                  rx *= s;
+                  ry *= s;
+               }
+               return ( std::erfc(w*(X(0)-cx-rx))*std::erfc(-w*(X(0)-cx+rx)) *
+                        std::erfc(w*(X(1)-cy-ry))*std::erfc(-w*(X(1)-cy+ry)) )/16;
+            }
+         }
+      }
+      case 2:
+      {
+         real_t x_ = X(0), y_ = X(1), rho, phi;
+         rho = std::hypot(x_, y_);
+         phi = atan2(y_, x_);
+         return pow(sin(M_PI*rho),2)*sin(3*phi);
+      }
+      case 3:
+      {
+         const real_t f = M_PI;
+         return sin(f*X(0))*sin(f*X(1));
+      }
+   }
+   return 0.0;
+}
+
+// Inflow boundary condition (zero for the problems considered in this example)
+real_t inflow_function(const Vector &x)
+{
+   switch (problem)
+   {
+      case 0:
+      case 1:
+      case 2:
+      case 3: return 0.0;
+   }
+   return 0.0;
+}

examples/ex41p.cpp

@@ -0,0 +1,666 @@
+//                                MFEM Example 41 - Parallel Version
+//
+// Compile with: make ex41p
+//
+// Sample runs: mpirun -np 4 ex41p
+//              mpirun -np 4 ex41p -p 0 -dt 0.01 -tf 10
+//              mpirun -np 4 ex41p -m ../data/periodic-hexagon.mesh -p 0 -dt 0.005 -tf 10
+//              mpirun -np 4 ex41p -m ../data/periodic-square.mesh -p 1 -dt 0.005 -tf 9
+//              mpirun -np 4 ex41p -m ../data/periodic-hexagon.mesh -p 1  -dt 0.005 -tf 9
+//              mpirun -np 4 ex41p -m ../data/star-q3.mesh -p 1 -rp 1 -dt 0.001 -tf 9
+//              mpirun -np 4 ex41p -m ../data/disc-nurbs.mesh -p 1 -rp 1 -dt 0.005 -tf 9
+//              mpirun -np 4 ex41p -m ../data/disc-nurbs.mesh -p 2 -rp 1 -dt 0.005 -tf 9
+//              mpirun -np 4 ex41p -m ../data/periodic-square.mesh -rp 2 -dt 0.0025 -tf 9 -vs 20
+//              mpirun -np 4 ex41p -m ../data/periodic-cube.mesh -p 0 -rs 2 -o 2 -dt 0.01 -tf 8
+//
+// Device sample runs:
+//
+// Description:  This example code solves the time-dependent advection-diffusion
+//               equation du/dt + v.grad(u) - a div(grad(u)) = 0, where v is a
+//               given fluid velocity, a is the diffusion coefficient, and
+//               u0(x)=u(0,x) is a given initial condition.
+//
+//               The example demonstrates the use of Discontinuous Galerkin (DG)
+//               bilinear forms in MFEM (face integrators), DG-LOR Preconditioning and the use of IMEX ODE time integrators.
+
+#include "mfem.hpp"
+
+using namespace std;
+using namespace mfem;
+
+// Choice for the problem setup. The fluid velocity, initial condition and
+// inflow boundary condition are chosen based on this parameter.
+int problem;
+
+// Velocity coefficient
+void velocity_function(const Vector &x, Vector &v);
+
+// Initial condition
+real_t u0_function(const Vector &x);
+
+// Inflow boundary condition
+real_t inflow_function(const Vector &x);
+
+// Mesh bounding box
+Vector bb_min, bb_max;
+
+
+class DG_Solver : public Solver
+{
+private:
+   HypreParMatrix &M, &K, &S;
+   HypreParMatrix *A;
+   CGSolver linear_solver;
+   real_t dt;
+   SparseMatrix M_diag;
+public:
+   DG_Solver(HypreParMatrix &M_, HypreParMatrix &K_, HypreParMatrix &S_,
+             const FiniteElementSpace &fes)
+      : M(M_),
+        K(K_),
+        S(S_),
+        A(nullptr),
+        linear_solver(M.GetComm()),
+        dt(1.0)
+   {
+      linear_solver.iterative_mode = false;
+      linear_solver.SetRelTol(1e-9);
+      linear_solver.SetAbsTol(0.0);
+      linear_solver.SetMaxIter(100);
+      linear_solver.SetPrintLevel(0);
+
+      M.GetDiag(M_diag);
+   }
+
+   void SetTimeStep(real_t dt_)
+   {
+      if (dt_ != dt)
+      {
+         delete A;
+         dt = dt_;
+         // // Form operator A = M + dt*S
+         A = Add(dt, S, 1.0, M);
+         linear_solver.SetOperator(*A);
+      }
+   }
+
+   void SetOperator(const Operator &op) override
+   {
+      linear_solver.SetOperator(op);
+   }
+
+   void Mult(const Vector &x, Vector &y) const override
+   {
+      linear_solver.Mult(x, y);
+   }
+
+   void SetPreconditioner(Solver &precond)
+   {
+      linear_solver.SetPreconditioner(precond);
+   }
+
+   ~DG_Solver() override
+   {
+      delete A;
+   }
+};
+
+/** A time-dependent operator for the right-hand side of the ODE. The DG weak
+    form of the advection-diffusion equation is (M + dt S) du/dt = Su - K u + b, where M and K are the mass
+    and advection matrices, and b describes the flow on the boundary. In the case of IMEX evolution, the diffusion term is treated
+    implicitly, and the advection term is treated explicitly.  */
+class IMEX_Evolution : public SplitTimeDependentOperator
+{
+private:
+   OperatorHandle M, K, S, A;
+   const Vector &b;
+   Solver *M_prec;
+   CGSolver M_solver;
+   DG_Solver *dg_solver;
+   LORSolver<HypreBoomerAMG>* lor_solver;
+
+   mutable Vector z;
+   mutable Vector w;
+
+public:
+   IMEX_Evolution(ParBilinearForm &M_, ParBilinearForm &K_, ParBilinearForm &S_,
+                  const Vector &b_, ParBilinearForm &A_);
+
+   virtual
+   ~IMEX_Evolution()
+   {
+      delete dg_solver;
+      delete lor_solver;
+      delete M_prec;
+   }
+
+   virtual
+   void Mult1(const Vector &x, Vector &y) const override;
+   virtual
+   void ImplicitSolve2(const real_t dt, const Vector &x, Vector &k) override;
+};
+
+
+int main(int argc, char *argv[])
+{
+
+   // 1. Initialize MPI and HYPRE.
+   Mpi::Init();
+   int num_procs = Mpi::WorldSize();
+   int myid = Mpi::WorldRank();
+   Hypre::Init();
+
+   // 2. Parse command-line options.
+   problem = 0;
+   const char *mesh_file = "../data/periodic-square.mesh";
+   int ser_ref_levels = 2;
+   int par_ref_levels = 0;
+   int order = 3;
+   bool pa = false;
+   bool ea = false;
+   bool fa = false;
+   const char *device_config = "cpu";
+   int ode_solver_type = 58;
+   real_t t_final = 10.0;
+   real_t dt = 0.001;
+   bool paraview = false;
+   int vis_steps = 50;
+   bool adios2 = false;
+   bool binary = false;
+   real_t diffusion_term = 0.01;
+   real_t kappa = -1.0;
+   real_t sigma = -1.0;
+   bool visualization = false;
+   //    #if MFEM_HYPRE_VERSION >= 21800
+   //    PrecType prec_type = PrecType::AIR;
+   // #else
+   //    PrecType prec_type = PrecType::ILU;
+   // #endif
+   int precision = 16;
+   cout.precision(precision);
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&problem, "-p", "--problem",
+                  "Problem setup to use. See options in velocity_function().");
+   args.AddOption(&ser_ref_levels, "-rs", "--refine-serial",
+                  "Number of times to refine the mesh uniformly in serial.");
+   args.AddOption(&par_ref_levels, "-rp", "--refine-parallel",
+                  "Number of times to refine the mesh uniformly in parallel.");
+   args.AddOption(&order, "-o", "--order",
+                  "Order (degree) of the finite elements.");
+   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(&fa, "-fa", "--full-assembly", "-no-fa",
+                  "--no-full-assembly", "Enable Full Assembly.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
+   args.AddOption(&ode_solver_type, "-s", "--ode-solver",
+                  SplitODESolver::Types.c_str());
+   args.AddOption(&t_final, "-tf", "--t-final",
+                  "Final time; start time is 0.");
+   args.AddOption(&dt, "-dt", "--time-step",
+                  "Time step.");
+   args.AddOption(&paraview, "-paraview", "--paraview-datafiles", "-no-paraview",
+                  "--no-paraview-datafiles",
+                  "Save data files for ParaView (paraview.org) visualization.");
+   args.AddOption(&adios2, "-adios2", "--adios2-streams", "-no-adios2",
+                  "--no-adios2-streams",
+                  "Save data using adios2 streams.");
+   args.AddOption(&binary, "-binary", "--binary-datafiles", "-ascii",
+                  "--ascii-datafiles",
+                  "Use binary (Sidre) or ascii format for VisIt data files.");
+   args.AddOption(&vis_steps, "-vs", "--visualization-steps",
+                  "Visualize every n-th timestep.");
+   args.Parse();
+   if (!args.Good())
+   {
+      if (Mpi::Root())
+      {
+         args.PrintUsage(cout);
+      }
+      return 1;
+   }
+   if (Mpi::Root())
+   {
+      args.PrintOptions(cout);
+   }
+   if (kappa < 0)
+   {
+      kappa = (order+1)*(order+1);
+   }
+   Device device(device_config);
+   if (Mpi::Root()) { device.Print(); }
+
+   // 3. Read the mesh from the given mesh file. We can handle geometrically
+   //    periodic meshes in this code.
+   Mesh *mesh = new Mesh(mesh_file);
+   const int dim = mesh->Dimension();
+
+   // 4. Define the Split ODE solver used for time integration. The IMEX solvers currently
+   // available are: 55 - Forward Backward Euler, 56 - IMEXRK2(2,2,2), 57 - IMEXRK2(2,3,2), and
+   // 58 - IMEX_DIRK_RK3.
+   unique_ptr<SplitODESolver> ode_solver = SplitODESolver::Select(ode_solver_type);
+
+   // 5. Refine the mesh to increase the resolution. In this example we do
+   //    'ref_levels' of uniform refinement, where 'ref_levels' is a
+   //    command-line parameter.
+   for (int lev = 0; lev < ser_ref_levels; lev++) { mesh->UniformRefinement(); }
+   if (mesh->NURBSext)
+   {
+      mesh->SetCurvature(max(order, 1));
+   }
+   mesh->GetBoundingBox(bb_min, bb_max, max(order, 1));
+
+
+   // 6. Define the parallel mesh by a partitioning of the serial mesh. Refine
+   //    this mesh further in parallel to increase the resolution. Once the
+   //    parallel mesh is defined, the serial mesh can be deleted.
+   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
+   delete mesh;
+   for (int lev = 0; lev < par_ref_levels; lev++)
+   {
+      pmesh->UniformRefinement();
+   }
+
+   // 7. Define the discontinuous DG finite element space of the given
+   //    polynomial order on the refined mesh.
+   DG_FECollection fec(order, dim, BasisType::GaussLobatto);
+   ParFiniteElementSpace *fes = new ParFiniteElementSpace(pmesh, &fec);
+
+   HYPRE_BigInt global_vSize = fes->GlobalTrueVSize();
+   if (Mpi::Root())
+   {
+      cout << "Number of unknowns: " << global_vSize << endl;
+   }
+
+   // 8. Set up and assemble the bilinear and linear forms corresponding to the
+   //    DG discretization. The DGTraceIntegrator involves integrals over mesh
+   //    interior faces.
+   VectorFunctionCoefficient velocity(dim, velocity_function);
+   FunctionCoefficient inflow(inflow_function);
+   ConstantCoefficient diff_coeff(diffusion_term);
+   ConstantCoefficient dt_diff_coeff(dt*diffusion_term);
+
+   ParBilinearForm *m = new ParBilinearForm(fes);
+   ParBilinearForm *k = new ParBilinearForm(fes);
+   ParBilinearForm *s = new ParBilinearForm(fes);
+
+   if (pa)
+   {
+      m->SetAssemblyLevel(AssemblyLevel::PARTIAL);
+      k->SetAssemblyLevel(AssemblyLevel::PARTIAL);
+      s->SetAssemblyLevel(AssemblyLevel::PARTIAL);
+   }
+   else if (ea)
+   {
+      m->SetAssemblyLevel(AssemblyLevel::ELEMENT);
+      k->SetAssemblyLevel(AssemblyLevel::ELEMENT);
+      s->SetAssemblyLevel(AssemblyLevel::ELEMENT);
+   }
+   else if (fa)
+   {
+      m->SetAssemblyLevel(AssemblyLevel::FULL);
+      k->SetAssemblyLevel(AssemblyLevel::FULL);
+      s->SetAssemblyLevel(AssemblyLevel::FULL);
+   }
+
+   m->AddDomainIntegrator(new MassIntegrator);
+
+   constexpr real_t alpha = -1.0;
+   k->AddDomainIntegrator(new ConvectionIntegrator(velocity, alpha));
+   k->AddInteriorFaceIntegrator(new NonconservativeDGTraceIntegrator(velocity,
+                                                                     alpha));
+   k->AddBdrFaceIntegrator(new NonconservativeDGTraceIntegrator(velocity, alpha));
+
+   s->AddDomainIntegrator(new DiffusionIntegrator(diff_coeff));
+   s->AddInteriorFaceIntegrator(new DGDiffusionIntegrator(diff_coeff, sigma,
+                                                          kappa));
+   s->AddBdrFaceIntegrator(new DGDiffusionIntegrator(diff_coeff, sigma, kappa));
+
+   ParLinearForm *b = new ParLinearForm(fes);
+   b->AddBdrFaceIntegrator(new BoundaryFlowIntegrator(inflow, velocity, alpha));
+
+   //For the preconditioner - create billinear form corresponding to operator (M + dt S)
+   ParBilinearForm *a = new ParBilinearForm(fes);
+   a->AddDomainIntegrator(new MassIntegrator);
+   a->AddDomainIntegrator(new DiffusionIntegrator(dt_diff_coeff));
+   a->AddInteriorFaceIntegrator(new DGDiffusionIntegrator(dt_diff_coeff, sigma,
+                                                          kappa));
+   a->AddBdrFaceIntegrator(new DGDiffusionIntegrator(dt_diff_coeff, sigma, kappa));
+
+   int skip_zeros = 0;
+   m->Assemble(skip_zeros);
+   k->Assemble(skip_zeros);
+   s->Assemble(skip_zeros);
+   a->Assemble();
+   b->Assemble();
+
+   m->Finalize(skip_zeros);
+   k->Finalize(skip_zeros);
+   s->Finalize(skip_zeros);
+   a->Finalize(skip_zeros);
+   HypreParVector *B = b->ParallelAssemble();
+
+   // 9. Define the initial conditions. Set up visualization (if desired).
+   FunctionCoefficient u0(u0_function);
+   ParGridFunction *u = new ParGridFunction(fes);
+   u->ProjectCoefficient(u0);
+   HypreParVector *U = u->GetTrueDofs();
+   {
+      ostringstream mesh_name, sol_name;
+      mesh_name << "ex41-mesh." << setfill('0') << setw(6) << myid;
+      sol_name << "ex41-init." << setfill('0') << setw(6) << myid;
+      ofstream omesh(mesh_name.str().c_str());
+      omesh.precision(precision);
+      pmesh->Print(omesh);
+      ofstream osol(sol_name.str().c_str());
+      osol.precision(precision);
+      u->Save(osol);
+   }
+   ParaViewDataCollection *pd = NULL;
+   if (paraview)
+   {
+      pd = new ParaViewDataCollection("Example41P", pmesh);
+      pd->SetPrefixPath("ParaView");
+      pd->RegisterField("solution", u);
+      pd->SetLevelsOfDetail(order);
+      pd->SetDataFormat(VTKFormat::BINARY);
+      pd->SetHighOrderOutput(true);
+      pd->SetCycle(0);
+      pd->SetTime(0.0);
+      pd->Save();
+   }
+   socketstream sout;
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      sout.open(vishost, visport);
+      if (!sout)
+      {
+         if (Mpi::Root())
+         {
+            cout << "Unable to connect to GLVis server at "
+                 << vishost << ':' << visport << endl;
+         }
+         visualization = false;
+         if (Mpi::Root())
+         {
+            cout << "GLVis visualization disabled.\n";
+         }
+      }
+      else
+      {
+         sout << "parallel " << num_procs << " " << myid << "\n";
+         sout.precision(precision);
+         sout << "solution\n" << *pmesh << *u;
+         sout << "pause\n";
+         sout << flush;
+         if (Mpi::Root())
+         {
+            cout << "GLVis visualization paused."
+                 << " Press space (in the GLVis window) to resume it.\n";
+         }
+      }
+   }
+
+   // 10. Define the time-dependent evolution operator describing the ODE
+   //    right-hand side, and perform time-integration (looping over the time
+   //    iterations, ti, with a time-step dt).
+   IMEX_Evolution adv(*m, *k, *s, *B, *a);
+
+   real_t t = 0.0;
+   adv.SetTime(t);
+   ode_solver->Init(adv);
+
+
+   bool done = false;
+   for (int ti = 0; !done; )
+   {
+      real_t dt_real = min(dt, t_final - t);
+      ode_solver->Step(*U, t, dt_real);
+      ti++;
+
+      done = (t >= t_final - 1e-8*dt);
+
+      if (done || ti % vis_steps == 0)
+      {
+         if (Mpi::Root())
+         {
+            cout << "time step: " << ti << ", time: " << t << endl;
+         }
+         *u = *U;
+         if (visualization)
+         {
+            sout << "parallel " << num_procs << " " << myid << "\n";
+            sout << "solution\n" << *pmesh << *u << flush;
+         }
+         if (paraview)
+         {
+            pd->SetCycle(ti);
+            pd->SetTime(t);
+            pd->Save();
+         }
+      }
+   }
+
+   {
+      *u = *U;
+      ostringstream sol_name;
+      sol_name << "ex41-final." << setfill('0') << setw(6) << myid;
+      ofstream osol(sol_name.str().c_str());
+      osol.precision(precision);
+      u->Save(osol);
+   }
+
+   // 11. Free the used memory.
+   delete U;
+   delete u;
+   delete B;
+   delete a;
+   delete b;
+   delete s;
+   delete k;
+   delete m;
+   delete fes;
+   delete pmesh;
+
+   return 0;
+}
+
+// Implementation of class IMEX_Evolution
+IMEX_Evolution::IMEX_Evolution(ParBilinearForm &M_, ParBilinearForm &K_,
+                               ParBilinearForm &S_, const Vector &b_, ParBilinearForm &A_)
+   : SplitTimeDependentOperator(M_.ParFESpace()->GetTrueVSize()), b(b_),
+     M_solver(M_.ParFESpace()->GetComm()), z(height), w(height)
+{
+   if (M_.GetAssemblyLevel()==AssemblyLevel::LEGACY)
+   {
+      M.Reset(M_.ParallelAssemble(), true);
+      K.Reset(K_.ParallelAssemble(), true);
+      S.Reset(S_.ParallelAssemble(), true);
+   }
+   else
+   {
+      M.Reset(&M_, false);
+      K.Reset(&K_, false);
+      S.Reset(&S_, false);
+   }
+
+   M_solver.SetOperator(*M);
+
+   Array<int> ess_tdof_list;
+   if (M_.GetAssemblyLevel() == AssemblyLevel::LEGACY)
+   {
+      A.Reset(A_.ParallelAssemble(), true);
+      HypreParMatrix &M_mat = *M.As<HypreParMatrix>();
+      HypreParMatrix &K_mat = *K.As<HypreParMatrix>();
+      HypreParMatrix &S_mat = *S.As<HypreParMatrix>();
+      HypreSmoother *hypre_prec = new HypreSmoother(M_mat, HypreSmoother::Jacobi);
+      M_prec = hypre_prec;
+
+      dg_solver = new DG_Solver(M_mat, K_mat, S_mat, *M_.FESpace());
+      lor_solver = new LORSolver<HypreBoomerAMG>(A_, ess_tdof_list);
+      lor_solver->GetSolver().SetSystemsOptions(A_.ParFESpace()->GetVDim(), true);
+      dg_solver -> SetPreconditioner(*lor_solver);
+   }
+   else
+   {
+      M_prec = new OperatorJacobiSmoother(M_, ess_tdof_list);
+      dg_solver = NULL;
+   }
+   M_solver.SetPreconditioner(*M_prec);
+   M_solver.iterative_mode = false;
+   M_solver.SetRelTol(1e-9);
+   M_solver.SetAbsTol(0.0);
+   M_solver.SetMaxIter(100);
+   M_solver.SetPrintLevel(0);
+}
+
+void IMEX_Evolution::Mult1(const Vector &x, Vector &y) const
+{
+   // Perform the explicit step
+   // y = M^{-1} (K x + b)
+   K->Mult(x, z);
+   z += b;
+   M_solver.Mult(z, y);
+}
+
+void IMEX_Evolution::ImplicitSolve2(const real_t dt, const Vector &x, Vector &k)
+{
+   // Perform the implicit step
+   // solve for k, k = -(M+dt S)^{-1} S x
+   MFEM_VERIFY(dg_solver != NULL,
+               "Implicit time integration is not supported with partial assembly");
+   S->Mult(x, z);
+   z*= -1.0;
+   dg_solver->SetTimeStep(dt);
+   dg_solver->Mult(z, k);
+}
+
+// Velocity coefficient
+void velocity_function(const Vector &x, Vector &v)
+{
+   int dim = x.Size();
+
+   // map to the reference [-1,1] domain
+   Vector X(dim);
+   for (int i = 0; i < dim; i++)
+   {
+      real_t center = (bb_min[i] + bb_max[i]) * 0.5;
+      X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
+   }
+
+   switch (problem)
+   {
+      case 0:
+      {
+         // Translations in 1D, 2D, and 3D
+         switch (dim)
+         {
+            case 1: v(0) = 1.0; break;
+            case 2: v(0) = sqrt(2./3.); v(1) = sqrt(1./3.); break;
+            case 3: v(0) = sqrt(3./6.); v(1) = sqrt(2./6.); v(2) = sqrt(1./6.);
+               break;
+         }
+         break;
+      }
+      case 1:
+      case 2:
+      {
+         // Clockwise rotation in 2D around the origin
+         const real_t w = M_PI/2;
+         switch (dim)
+         {
+            case 1: v(0) = 1.0; break;
+            case 2: v(0) = w*X(1); v(1) = -w*X(0); break;
+            case 3: v(0) = w*X(1); v(1) = -w*X(0); v(2) = 0.0; break;
+         }
+         break;
+      }
+      case 3:
+      {
+         // Clockwise twisting rotation in 2D around the origin
+         const real_t w = M_PI/2;
+         real_t d = max((X(0)+1.)*(1.-X(0)),0.) * max((X(1)+1.)*(1.-X(1)),0.);
+         d = d*d;
+         switch (dim)
+         {
+            case 1: v(0) = 1.0; break;
+            case 2: v(0) = d*w*X(1); v(1) = -d*w*X(0); break;
+            case 3: v(0) = d*w*X(1); v(1) = -d*w*X(0); v(2) = 0.0; break;
+         }
+         break;
+      }
+   }
+}
+
+// Initial condition
+real_t u0_function(const Vector &x)
+{
+   int dim = x.Size();
+
+   // map to the reference [-1,1] domain
+   Vector X(dim);
+   for (int i = 0; i < dim; i++)
+   {
+      real_t center = (bb_min[i] + bb_max[i]) * 0.5;
+      X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
+   }
+
+   switch (problem)
+   {
+      case 0:
+      case 1:
+      {
+         switch (dim)
+         {
+            case 1:
+               return exp(-40.*pow(X(0)-0.5,2));
+            case 2:
+            case 3:
+            {
+               real_t rx = 0.45, ry = 0.25, cx = 0., cy = -0.2, w = 10.;
+               if (dim == 3)
+               {
+                  const real_t s = (1. + 0.25*cos(2*M_PI*X(2)));
+                  rx *= s;
+                  ry *= s;
+               }
+               return ( std::erfc(w*(X(0)-cx-rx))*std::erfc(-w*(X(0)-cx+rx)) *
+                        std::erfc(w*(X(1)-cy-ry))*std::erfc(-w*(X(1)-cy+ry)) )/16;
+            }
+         }
+      }
+      case 2:
+      {
+         real_t x_ = X(0), y_ = X(1), rho, phi;
+         rho = std::hypot(x_, y_);
+         phi = atan2(y_, x_);
+         return pow(sin(M_PI*rho),2)*sin(3*phi);
+      }
+      case 3:
+      {
+         const real_t f = M_PI;
+         return sin(f*X(0))*sin(f*X(1));
+      }
+   }
+   return 0.0;
+}
+
+// Inflow boundary condition (zero for the problems considered in this example)
+real_t inflow_function(const Vector &x)
+{
+   switch (problem)
+   {
+      case 0:
+      case 1:
+      case 2:
+      case 3: return 0.0;
+   }
+   return 0.0;
+}

examples/makefile

@@ -22,14 +22,14 @@ MFEM_LIB_FILE = mfem_is_not_built
 
 SEQ_EXAMPLES = ex0 ex1 ex2 ex3 ex4 ex5 ex6 ex7 ex8 ex9 ex10 ex14 ex15 ex16 \
    ex17 ex18 ex19 ex20 ex21 ex22 ex23 ex24 ex25 ex26 ex27 ex28 ex29 ex30 \
-   ex31 ex33 ex34 ex36 ex37 ex38 ex39 ex40
+   ex31 ex33 ex34 ex36 ex37 ex38 ex39 ex40 ex41 
 PAR_EXAMPLES = ex0p ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex8p ex9p ex10p ex11p \
    ex12p ex13p ex14p ex15p ex16p ex17p ex18p ex19p ex20p ex21p ex22p ex24p \
    ex25p ex26p ex27p ex28p ex29p ex30p ex31p ex32p ex33p ex34p ex35p ex36p \
-   ex37p ex39p ex40p
+   ex37p ex39p ex40p ex41p
 SEQ_DEVICE_EXAMPLES = ex1 ex3 ex4 ex5 ex6 ex9 ex14 ex22 ex24 ex25 ex26 ex34
 PAR_DEVICE_EXAMPLES = ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex9p ex13p ex14p \
-   ex22p ex24p ex25p ex26p ex34p ex35p
+   ex22p ex24p ex25p ex26p ex34p ex35p ex41p
 
 ifeq ($(MFEM_USE_LAPACK),YES)
    SEQ_EXAMPLES += ex38

fem/CMakeLists.txt

@@ -241,8 +241,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

@@ -3459,6 +3459,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.cpp

@@ -48,7 +48,7 @@ void LORBase::AddIntegratorsAndMarkers(BilinearForm &a_from,
    for (int i=0; i<integrators->Size(); ++i)
    {
       BilinearFormIntegrator *integrator = (*integrators)[i];
-      if (markers[i] != nullptr)
+      if (markers[i])
       {
          (a_to.*add_integrator_marker)(integrator, *markers[i]);
       }

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,177 @@ 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 = [&]()
+   {
+#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);
+   }
+
+   // The second argument (nrows + face_nbr_vsize) accounts for additional
+   // columns contributed by DG face neighbors in parallel finite element
+   // spaces. In serial, face_nbr_vsize is set to 0.
+   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 +550,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 +608,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 +643,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 +691,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
 
@@ -504,12 +736,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

@@ -481,6 +481,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); }

linalg/ode.cpp

@@ -28,6 +28,11 @@ std::string ODESolver::ImplicitTypes  =
    "        GA      : 40 -- 50  - Generalized-alpha,\n\t"
    "        AM      : 51 - AM1, 52 - AM2, 53 - AM3, 54 - AM4\n";
 
+std::string ODESolver::IMEXTypes =
+   "\n\tIMEX solver: \n\t"
+   "           55 - Forward Backward Euler, 56 - IMEXRK2(2,2,2), 57 - IMEXRK2(2,3,2)\n\t"
+   "           58 - IMEX_DIRK_RK3\n";
+
 std::string ODESolver::Types = ODESolver::ExplicitTypes +
                                ODESolver::ImplicitTypes;
 
@@ -41,6 +46,10 @@ std::unique_ptr<ODESolver> ODESolver::Select(int ode_solver_type)
    {
       return SelectImplicit(ode_solver_type);
    }
+   // else
+   // {
+   //    return SelectIMEX(ode_solver_type);
+   // }
 }
 
 std::unique_ptr<ODESolver> ODESolver::SelectExplicit(int ode_solver_type)
@@ -106,6 +115,20 @@ std::unique_ptr<ODESolver> ODESolver::SelectImplicit(int ode_solver_type)
    }
 }
 
+std::unique_ptr<SplitODESolver> SplitODESolver::Select(int ode_solver_type)
+{
+   using ode_ptr = std::unique_ptr<SplitODESolver>;
+   switch (ode_solver_type)
+   {
+      case 55: return ode_ptr(new IMEXExpImplEuler);
+      case 56: return ode_ptr(new IMEXRK2);
+      case 57: return ode_ptr(new IMEXRK2_3StageExplicit);
+      case 58: return ode_ptr(new IMEX_DIRK_RK3);
+      default: MFEM_ABORT("Unknown ODE solver type: " << ode_solver_type );
+   }
+}
+
+
 
 void ODEStateDataVector::SetSize( int vsize, MemoryType m_t)
 {
@@ -1277,4 +1300,195 @@ void GeneralizedAlpha2Solver::Step(Vector &x, Vector &dxdt,
    t += dt;
 }
 
+
+void SplitODESolver::Init(SplitTimeDependentOperator &f_)
+{
+   this->f = &f_;
+   mem_type = GetMemoryType(f_.GetMemoryClass());
 }
+
+void IMEXExpImplEuler::Init(SplitTimeDependentOperator &f_)
+{
+   SplitODESolver::Init(f_);
+   int n = f->Width();
+   k1.SetSize(n, mem_type);
+   k2.SetSize(n, mem_type);
+}
+
+void IMEXExpImplEuler::Step(Vector &x, real_t &t, real_t &dt)
+{
+   f->SetTime(t);
+   f->Mult1(x, k1);
+   f->SetTime(t+dt);
+   f->ImplicitSolve2(dt, x, k2);
+   f->SetTime(t);
+   x.Add(dt, k1);
+   x.Add(dt, k2);
+   t += dt;
+}
+
+void IMEXRK2::Init(SplitTimeDependentOperator &f_)
+{
+   SplitODESolver::Init(f_);
+   int n = f->Width();
+   k1_exp.SetSize(n, mem_type);
+   k2_exp.SetSize(n, mem_type);
+   k2_imp.SetSize(n, mem_type);
+   k3_imp.SetSize(n, mem_type);
+   y.SetSize(n, mem_type);
+   z.SetSize(n, mem_type);
+}
+
+void IMEXRK2::Step(Vector &x, real_t &t, real_t &dt)
+{
+   double gamma = 1 - sqrt(2)/2;
+   double delta = 1 - 1/(2*gamma);
+
+   f->SetTime(t);
+
+   //K1 exp is just f_1(t, x)
+   f->Mult1(x, k1_exp);
+
+   //K2 exp is f_1(t + gamma dt, x + dt gamma K1)
+   f->SetTime(t + gamma*dt);
+   add(x, dt*gamma, k1_exp, y);
+   f->Mult1(y, k2_exp);
+
+   //K2_imp = f_2(t + gamma dt, x + dt gamma K2_imp)
+   f->ImplicitSolve2(dt*gamma, x, k2_imp);
+
+   //K3_imp = f_2(t+dt,x + dt(1-gamma)K2_imp + dt gamma K3_imp)
+   f -> SetTime(t + dt);
+   add(x, dt*(1-gamma), k2_imp, z);
+   f->ImplicitSolve2(dt*gamma, z, k3_imp);
+
+   //add it all up
+   x.Add(dt*delta, k1_exp);
+   x.Add(dt*(1-delta), k2_exp);
+   x.Add(dt*(1-gamma), k2_imp);
+   x.Add(dt*gamma, k3_imp);
+   t += dt;
+}
+
+void IMEXRK2_3StageExplicit::Init(SplitTimeDependentOperator &f_)
+{
+   SplitODESolver::Init(f_);
+   int n = f->Width();
+   k1_exp.SetSize(n, mem_type);
+   k2_exp.SetSize(n, mem_type);
+   k3_exp.SetSize(n, mem_type);
+   k2_imp.SetSize(n, mem_type);
+   k3_imp.SetSize(n, mem_type);
+   y.SetSize(n, mem_type);
+   z.SetSize(n, mem_type);
+   w.SetSize(n, mem_type);
+}
+
+void IMEXRK2_3StageExplicit::Step(Vector &x, real_t &t, real_t &dt)
+{
+   // WIP
+   double gamma = 1 - sqrt(2)/2;
+   double delta = -2*sqrt(2)/3;
+
+   f->SetTime(t);
+
+   //K1 exp is just f_1(t, x)
+   f->Mult1(x, k1_exp);
+
+   //K2 exp is f_1(t + gamma dt, x + dt gamma K1)
+   f->SetTime(t + gamma*dt);
+   add(x, dt*gamma, k1_exp, y);
+   f->Mult1(y, k2_exp);
+
+   //K3 Exp is f_1(t + dt, x + dt gamma K1_exp + dt (1-gamma) K2_exp)
+   f->SetTime(t + dt);
+   add(x, dt*delta, k1_exp, y);
+   add(y, dt*(1-delta), k2_exp, w);
+   f->Mult1(w, k3_exp);
+
+   //K2_imp = f_2(t + gamma dt, x + dt gamma K2_imp)
+   f->SetTime(t + gamma*dt);
+   f->ImplicitSolve2(dt*gamma, x, k2_imp);
+
+   //K3_imp = f_2(t+dt,x + dt(1-gamma)K2_imp + dt gamma K3_imp)
+   f -> SetTime(t + dt);
+   add(x, dt*(1-gamma), k2_imp, z);
+   f->ImplicitSolve2(dt*gamma, z, k3_imp);
+
+   //add it all up
+   x.Add(dt*delta, k2_exp);
+   x.Add(dt*(1-delta), k3_exp);
+   x.Add(dt*(1-gamma), k2_imp);
+   x.Add(dt*gamma, k3_imp);
+   t += dt;
+}
+
+void IMEX_DIRK_RK3::Init(SplitTimeDependentOperator &f_)
+{
+   SplitODESolver::Init(f_);
+   int n = f->Width();
+   k1_exp.SetSize(n, mem_type);
+   k2_exp.SetSize(n, mem_type);
+   k3_exp.SetSize(n, mem_type);
+   k4_exp.SetSize(n, mem_type);
+   k2_imp.SetSize(n, mem_type);
+   k3_imp.SetSize(n, mem_type);
+   k4_imp.SetSize(n, mem_type);
+   y.SetSize(n, mem_type);
+   z.SetSize(n, mem_type);
+   w.SetSize(n, mem_type);
+   v.SetSize(n, mem_type);
+   u.SetSize(n, mem_type);
+}
+
+void IMEX_DIRK_RK3::Step(Vector &x, real_t &t, real_t &dt)
+{
+   double gamma = 0.4358665215;
+   double b1 = 1.208496649;
+   double b2 = -0.644363171;
+   double a_31 = 0.3212788860;
+   double a_32 = 0.3966543747;
+   double a_41 = -0.105858296;
+   double a_42 = 0.5529291479;
+   double a_43 = 0.5529291479;
+
+   //K1_exp
+   f->SetTime(t);
+   f->Mult1(x, k1_exp);
+
+   //K2_imp, K2_exp
+   f->SetTime(t + gamma*dt);
+   add(x, dt*gamma, k1_exp, y);
+   f->Mult1(y, k2_exp);
+   f->ImplicitSolve2(dt*gamma, x, k2_imp);
+
+   //K3_imp, K3_exp
+   f->SetTime(t + (1+gamma)/2*dt);
+   add(x, dt*a_31, k1_exp, y);
+   add(y, dt*a_32, k2_exp, w);
+   f->Mult1(w, k3_exp);
+   add(x, dt*(1-gamma)/2, k2_imp, z);
+   f->ImplicitSolve2(dt*gamma, z, k3_imp);
+
+   //K4_imp, K4_exp
+   f->SetTime(t+dt);
+   add(x, dt*a_41, k1_exp, y);
+   add(y, dt*a_42, k2_exp, w);
+   add(w, dt*a_43, k3_exp, v);
+   f->Mult1(v, k4_exp);
+   add(x, dt*b1, k2_imp, z);
+   add(z, dt*b2, k3_imp, u);
+   f->ImplicitSolve2(dt*gamma, u, k4_imp);
+
+   //add it all together
+   x.Add(dt*b1, k2_exp);
+   x.Add(dt*b2, k3_exp);
+   x.Add(dt*gamma, k4_exp);
+   x.Add(dt*b1, k2_imp);
+   x.Add(dt*b2, k3_imp);
+   x.Add(dt*gamma, k4_imp);
+   t += dt;
+
+}
+
+}

linalg/ode.hpp

@@ -184,6 +184,7 @@ public:
    // Help info for ODESolver options
    static MFEM_EXPORT std::string ExplicitTypes;
    static MFEM_EXPORT std::string ImplicitTypes;
+   static MFEM_EXPORT std::string IMEXTypes;
    static MFEM_EXPORT std::string Types;
 
    /// Function for selecting the desired ODESolver (Explicit and Implicit)
@@ -203,6 +204,12 @@ public:
    static MFEM_EXPORT std::unique_ptr<ODESolver> SelectImplicit(
       const int ode_solver_type);
 
+
+   /// Function for selecting the desired IMEX ODESolver
+   /// Returns an ODESolver pointer based on an type
+   /// Caller gets ownership of the object and is responsible for its deletion
+   //static MFEM_EXPORT std::unique_ptr<SplitODESolver> SelectIMEX(const int ode_solver_type);
+
    virtual ~ODESolver() { }
 };
 
@@ -931,6 +938,91 @@ public:
 
 };
 
+/// Class for solving systems of split ODEs: dx/dt = f_1(x,t) + f_2(x,t)
+class SplitODESolver : public ODESolver
+{
+protected:
+   /// Pointer to the associated SplitTimeDependentOperator.
+   SplitTimeDependentOperator *f;  // f(.,t) : R^n --> R^n
+   MemoryType mem_type;
+
+public:
+   SplitODESolver() : f(NULL) { mem_type = Device::GetHostMemoryType(); }
+
+   /// Associate a SplitTimeDependentOperator with the ODE solver. Overrides Init from ODESolver
+   /** This method has to be called:
+       - Before the first call to Step().
+       - When the dimensions of the associated SplitTimeDependentOperator change.
+       - When a time stepping sequence has to be restarted.
+       - To change the associated SplitTimeDependentOperator. */
+   virtual void Init(SplitTimeDependentOperator &f);
+
+   static MFEM_EXPORT std::unique_ptr<SplitODESolver> Select(
+      const int ode_solver_type);
+
+   virtual ~SplitODESolver() { }
+};
+
+class IMEXExpImplEuler : public SplitODESolver
+{
+private:
+   Vector k1; Vector k2;
+public:
+   void Init(SplitTimeDependentOperator &f_) override;
+
+   void Step(Vector &x, real_t &t, real_t &dt) override;
+};
+
+
+/**IMEX RK2 Method from "On the Stability of IMEX Upwind gSBP Schemes for 1D Linear Advection‑Difusion Equations" by Sigrun Ortleb.
+ * Same as (2,2,2) from "Implicit-explicit Runge-Kutta
+    methods for time-dependent partial differential equations" by Ascher, Ruuth
+    and Spiteri, Applied Numerical Mathematics (1997).**/
+class IMEXRK2 : public SplitODESolver
+{
+private:
+   Vector k1_exp; Vector k2_exp; Vector k2_imp; Vector k3_imp;
+   //helper vectors
+   Vector y; Vector z;
+public:
+   void Init(SplitTimeDependentOperator &f_) override;
+
+   void Step(Vector &x, real_t &t, real_t &dt) override;
+};
+
+
+/**(2,3,2) from "Implicit-explicit Runge-Kutta
+    methods for time-dependent partial differential equations" by Ascher, Ruuth
+    and Spiteri, Applied Numerical Mathematics (1997).**/
+class IMEXRK2_3StageExplicit : public SplitODESolver
+{
+private:
+   Vector k1_exp; Vector k2_exp; Vector k3_exp; Vector k2_imp; Vector k3_imp;
+   //helper vectors
+   Vector y; Vector z; Vector w;
+public:
+   void Init(SplitTimeDependentOperator &f_) override;
+
+   void Step(Vector &x, real_t &t, real_t &dt) override;
+};
+
+/**(3,4,3) from "Implicit-explicit Runge-Kutta
+    methods for time-dependent partial differential equations" by Ascher, Ruuth
+    and Spiteri, Applied Numerical Mathematics (1997).**/
+class IMEX_DIRK_RK3 : public SplitODESolver
+{
+private:
+   Vector k1_exp; Vector k2_exp; Vector k3_exp; Vector k4_exp;
+   Vector k2_imp; Vector k3_imp; Vector k4_imp;
+   //helper vectors
+   Vector y; Vector z; Vector w; Vector u; Vector v;
+public:
+   void Init(SplitTimeDependentOperator &f_) override;
+
+   void Step(Vector &x, real_t &t, real_t &dt) override;
+};
+
+
 }
 
-#endif
+#endif

linalg/operator.cpp

@@ -365,6 +365,29 @@ void SecondOrderTimeDependentOperator::ImplicitSolve(const real_t dt0,
    mfem_error("SecondOrderTimeDependentOperator::ImplicitSolve() is not overridden!");
 }
 
+void SplitTimeDependentOperator::Mult1(const Vector &, Vector &) const
+{
+   mfem_error("SplitTimeDependentOperator::Mult1() is not overridden!");
+}
+
+void SplitTimeDependentOperator::ImplicitSolve1(const real_t, const Vector &,
+                                                Vector &)
+{
+   mfem_error("SplitTimeDependentOperator::ImplicitSolve1() is not overridden!");
+}
+
+void SplitTimeDependentOperator::Mult2(const Vector &, Vector &) const
+{
+   mfem_error("SplitTimeDependentOperator::Mult2() is not overridden!");
+}
+
+void SplitTimeDependentOperator::ImplicitSolve2(const real_t, const Vector &,
+                                                Vector &)
+{
+   mfem_error("SplitTimeDependentOperator::ImplicitSolve2() is not overridden!");
+}
+
+
 SumOperator::SumOperator(const Operator *A, const real_t alpha,
                          const Operator *B, const real_t beta,
                          bool ownA, bool ownB)
@@ -911,4 +934,4 @@ real_t PowerMethod::EstimateLargestEigenvalue(Operator& opr, Vector& v0,
    return eigenvalue;
 }
 
-}
+}

linalg/operator.hpp

@@ -786,6 +786,37 @@ public:
    virtual ~SecondOrderTimeDependentOperator() { }
 };
 
+// F(y,\frac{dy}{dt},t) = G1(y,t) + G2(y,t)
+class SplitTimeDependentOperator : public TimeDependentOperator
+{
+   // virtual void Mult(...) sum of Mult1 and Mult2
+public:
+   /** @brief Construct a "square" SplitTimeDependentOperator
+   y = F(x,dxdt,t), where x, dxdt and y have the same dimension @a n. */
+   explicit SplitTimeDependentOperator(int n = 0, real_t t_ = 0.0,
+                                       Type type_ = EXPLICIT)
+      : TimeDependentOperator(n, t_,type_) { }
+
+   /** @brief Construct a SplitTimeDependentOperator y = F(x,dxdt,t),
+   where x, dxdt and y have the same dimension @a n. */
+   SplitTimeDependentOperator(int h, int w, real_t t_ = 0.0,
+                              Type type_ = EXPLICIT)
+      : TimeDependentOperator(h, w, t_,type_) { }
+   virtual void Mult1(const Vector &u, Vector &k) const; //M^{-1} G1
+   virtual void Mult2(const Vector &u, Vector &k) const; //M^{-1} G2
+   virtual void ImplicitSolve1(const real_t gamma, const Vector &u,
+                               Vector &k);  //F(u + gamma k, k, t) = G1(u + gamma k, t)
+   virtual void ImplicitSolve2(const real_t gamma, const Vector &u,
+                               Vector &k);  //F(u + gamma k, k, t) = G2(u + gamma k, t)
+
+   virtual ~SplitTimeDependentOperator() { }
+
+   // - Fully explicit, need Mult1 and Mult2 (or just plain Mult)
+   // - IMEX, need Mult1 and ImplicitSolve2
+   // - "Split implicit", need ImplicitSolve1 and ImplicitSolve2
+   // - Fully implicit, need ImplicitSolve
+};
+
 
 /// Base class for solvers
 class Solver : public Operator
@@ -1176,4 +1207,4 @@ public:
 
 }
 
-#endif
+#endif

miniapps/CMakeLists.txt

@@ -36,6 +36,4 @@ add_subdirectory(parelag)
 add_subdirectory(tribol)
 add_subdirectory(hooke)
 add_subdirectory(dpg)
-add_subdirectory(hdiv-linear-solver)
-add_subdirectory(dfem)
-add_subdirectory(diag-smoothers)
+add_subdirectory(hdiv-linear-solver)

miniapps/mtop/CMakeLists.txt

@@ -10,22 +10,21 @@
 # CONTRIBUTING.md for details.
 
 list(APPEND SEQMTOP_COMMON_SOURCES
-	paramnonlinearform.cpp
-        mtop_integrators.cpp)
+	darcy_heat_transfer_ex.cpp)
 
-list(APPEND SEQMTOP_COMMON_HEADERS
-	paramnonlinearform.hpp
-        mtop_integrators.hpp)
+# list(APPEND SEQMTOP_COMMON_HEADERS
+# 	paramnonlinearform.hpp
+#         mtop_integrators.hpp)
 
 convert_filenames_to_full_paths(SEQMTOP_COMMON_SOURCES)
-convert_filenames_to_full_paths(SEQMTOP_COMMON_HEADERS)
+//convert_filenames_to_full_paths(SEQMTOP_COMMON_HEADERS)
 
 set(SEQMTOP_COMMON_FILES
 	EXTRA_SOURCES ${SEQMTOP_COMMON_SOURCES}
 	EXTRA_HEADERS ${SEQMTOP_COMMON_HEADERS})
 
 add_mfem_miniapp(seqheat
-	MAIN seqheat.cpp
+	MAIN darcy_heat_transfer_ex.cpp
 	${SEQMTOP_COMMON_FILES}
 	LIBRARIES mfem)
 
@@ -51,4 +50,4 @@ add_mfem_miniapp(parheat
      ${PARMTOP_COMMON_FILES}
      LIBRARIES mfem)
 
-endif ()
+endif ()

miniapps/mtop/darcy_heat_transfer_ex.cpp

@@ -0,0 +1,758 @@
+//                               MFEM Darcy Test Run
+//
+// Compile with: make darcy_heat_transfer_ex
+//
+//
+// Description:  This code performs the forward and backward adjoint solve for advection diffusion, where the velocity field is given by Darcy
+
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+
+using namespace std;
+using namespace mfem;
+
+
+// *****Funtion definitions for the Advection-Diffusion solve******
+
+// Velocity coefficient
+void velocity_function(const Vector &x, Vector &v);
+
+// Initial condition
+double theta0_function(const Vector &x);
+
+// true solution
+real_t theta_exact(const Vector &x, real_t t);
+
+// rhs
+double forcing_function(const Vector &x, real_t t);
+
+// inflow
+double inflow_function(const Vector &x);
+
+real_t f_natural(const Vector & x);
+
+// Mesh bounding box
+Vector bb_min, bb_max;
+
+class DG_Solver : public Solver
+{
+private:
+   SparseMatrix &M, &K, &S, A;
+   CGSolver linear_solver;
+   BlockILU prec;
+   real_t dt;
+public:
+   DG_Solver(SparseMatrix &M_, SparseMatrix &K_, SparseMatrix &S_,
+             const FiniteElementSpace &fes)
+      : M(M_),
+        K(K_),
+        S(S_),
+        prec(fes.GetTypicalFE()->GetDof(),
+             BlockILU::Reordering::MINIMUM_DISCARDED_FILL),
+        dt(1.0)
+   {
+      linear_solver.iterative_mode = false;
+      linear_solver.SetRelTol(1e-9);
+      linear_solver.SetAbsTol(0.0);
+      linear_solver.SetMaxIter(100);
+      linear_solver.SetPrintLevel(0);
+      linear_solver.SetPreconditioner(prec);
+   }
+
+   void SetTimeStep(real_t dt_)
+   {
+      if (dt_ != dt)
+      {
+         dt = dt_;
+         // Form operator A = M + dt*S
+         A = S;
+         A *= dt;
+         A += M;
+
+         // this will also call SetOperator on the preconditioner
+         linear_solver.SetOperator(A);
+      }
+   }
+
+   void SetOperator(const Operator &op) override
+   {
+      linear_solver.SetOperator(op);
+   }
+
+   void Mult(const Vector &x, Vector &y) const override
+   {
+      linear_solver.Mult(x, y);
+   }
+};
+
+
+/** A time-dependent operator for the right-hand side of the ODE. The DG weak
+    form of the advection-diffusion equation is (M + dt S) du/dt = Su - K u + b, where M and K are the mass
+    and advection matrices, and b describes the flow on the boundary. In the case of IMEX evolution, the diffusion term is treated
+    implicitly, and the advection term is treated explicitly.  */
+class IMEX_Evolution : public SplitTimeDependentOperator
+{
+private:
+   BilinearForm &M, &K, &S;
+   const Vector &b;
+   unique_ptr<Solver> M_prec;
+   CGSolver M_solver;
+   unique_ptr<DG_Solver> dg_solver;
+
+   mutable Vector z;
+
+public:
+   IMEX_Evolution(BilinearForm &M_, BilinearForm &K_, BilinearForm &S_,
+                  const Vector &b_);
+
+   void Mult1(const Vector &x, Vector &y) const;
+   void ImplicitSolve2(const real_t dt, const Vector &x, Vector &k) override;
+};
+
+// *****Define the analytical solution and forcing terms / boundary conditions for Darcy*****
+void uFun_ex(const Vector & x, Vector & u);
+real_t pFun_ex(const Vector & x);
+void fFun(const Vector & x, Vector & f);
+real_t gFun(const Vector & x);
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   const char *mesh_file =
+      "square-extended.mesh"; //reference square, but extended to be [-1, 1] x [-1, 1]
+   int order_darcy = 1;
+   int ref_levels = 2;
+   int order_ad = 3;
+   int ode_solver_type = 55;
+   double t_final = 10.0;
+   double d_coef = 0.01;
+   double dt = 0.01;
+   double sigma = -1.0;
+   double kappa = -1.0;
+   bool visualization = true;
+   bool visit = false;
+   bool binary = false;
+   int vis_steps = 5;
+   bool paraview = false;
+   int precision = 16;
+   const char *device_config = "cpu";
+   cout.precision(precision);
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&mesh_file, "-m", "--mesh",
+                  "Mesh file to use.");
+   args.AddOption(&ref_levels, "-r", "--refine",
+                  "Number of times to refine the mesh uniformly.");
+   args.AddOption(&order_darcy, "-od", "--order_darcy",
+                  "Order (degree) of the finite elements for darcy solve.");
+   args.AddOption(&order_ad, "-oad", "--order_ad",
+                  "Order (degree) of the finite elements for advection diffusion.");
+   args.AddOption(&ode_solver_type, "-s", "--ode-solver",
+                  "55 - Forward Backward Euler, 56 - IMEXRK2(2,2,2), 57 - IMEXRK2(2,3,2), 58 - IMEX_DIRK_RK3\n");
+   args.AddOption(&t_final, "-tf", "--t-final",
+                  "Final time; start time is 0.");
+   args.AddOption(&dt, "-dt", "--time-step",
+                  "Time step.");
+   args.AddOption(&d_coef, "-d", "--diff-coef",
+                  "Diffusion coefficient.");
+   args.AddOption(&sigma, "-s", "--sigma",
+                  "One of the two DG penalty parameters, typically +1/-1."
+                  " See the documentation of class DGDiffusionIntegrator.");
+   args.AddOption(&kappa, "-k", "--kappa",
+                  "One of the two DG penalty parameters, should be positive."
+                  " Negative values are replaced with (order+1)^2.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.AddOption(&visit, "-visit", "--visit-datafiles", "-no-visit",
+                  "--no-visit-datafiles",
+                  "Save data files for VisIt (visit.llnl.gov) visualization.");
+   args.AddOption(&binary, "-binary", "--binary-datafiles", "-ascii",
+                  "--ascii-datafiles",
+                  "Use binary (Sidre) or ascii format for VisIt data files.");
+   args.AddOption(&vis_steps, "-vs", "--visualization-steps",
+                  "Visualize every n-th timestep.");
+   args.AddOption(&paraview, "-paraview", "--paraview-datafiles", "-no-paraview",
+                  "--no-paraview-datafiles",
+                  "Save data files for ParaView (paraview.org) visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   if (kappa < 0)
+   {
+      kappa = (order_ad+1)*(order_ad+1);
+   }
+   args.PrintOptions(cout);
+   Device device(device_config);
+   device.Print();
+
+   // 2. Define the ODE solver used for time integration. Several explicit, implicit and IMEX
+   //    Runge-Kutta methods are available.
+   unique_ptr<SplitODESolver> ode_solver = SplitODESolver::Select(ode_solver_type);
+   unique_ptr<SplitODESolver> ode_solver_adj = SplitODESolver::Select(
+                                                  ode_solver_type);
+   // 3. Read the mesh from the given mesh file.
+   Mesh mesh(mesh_file, 1, 1);
+   int dim = mesh.Dimension();
+
+   // 4. Refine the mesh in serial to increase the resolution. In this example
+   //    we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
+   //    a command-line parameter.
+   for (int lev = 0; lev < ref_levels; lev++) {mesh.UniformRefinement();}
+   if (mesh.NURBSext) {mesh.SetCurvature(max(order_ad, 1));}
+   mesh.GetBoundingBox(bb_min, bb_max, max(order_ad, 1));
+
+   // ********DARCY SOLVE
+   // 5. Define a finite element space on the mesh. Here we use the
+   //    Raviart-Thomas finite elements of the specified order.
+   FiniteElementCollection *hdiv_coll(new RT_FECollection(order_darcy, dim));
+   FiniteElementCollection *l2_coll(new L2_FECollection(order_darcy, dim));
+
+   FiniteElementSpace *R_space = new FiniteElementSpace(&mesh, hdiv_coll);
+   FiniteElementSpace *W_space = new FiniteElementSpace(&mesh, l2_coll);
+
+   // 6. Define the BlockStructure of the problem, i.e. define the array of
+   //    offsets for each variable. The last component of the Array is the sum
+   //    of the dimensions of each block.
+   Array<int> block_offsets(3); // number of variables + 1
+   block_offsets[0] = 0;
+   block_offsets[1] = R_space->GetVSize();
+   block_offsets[2] = W_space->GetVSize();
+   block_offsets.PartialSum();
+
+   std::cout << "***********************************************************\n";
+   std::cout << "dim(R) = " << block_offsets[1] - block_offsets[0] << "\n";
+   std::cout << "dim(W) = " << block_offsets[2] - block_offsets[1] << "\n";
+   std::cout << "dim(R+W) = " << block_offsets.Last() << "\n";
+   std::cout << "***********************************************************\n";
+
+   // 7. Define the coefficients, analytical solution, and rhs of the Darcy PDE.
+   ConstantCoefficient one(1.0);
+
+   VectorFunctionCoefficient fcoeff(dim, fFun);
+   FunctionCoefficient fnatcoeff(f_natural);
+   FunctionCoefficient gcoeff(gFun);
+
+   VectorFunctionCoefficient ucoeff(dim, uFun_ex);
+   FunctionCoefficient pcoeff(pFun_ex);
+
+   // 8. Allocate memory for solution and rhs of Darcy
+   MemoryType mt = device.GetMemoryType();
+   BlockVector x(block_offsets, mt), rhs(block_offsets, mt);
+
+   LinearForm *fform(new LinearForm);
+   fform->Update(R_space, rhs.GetBlock(0), 0);
+   fform->AddDomainIntegrator(new VectorFEDomainLFIntegrator(fcoeff));
+   fform->AddBoundaryIntegrator(new VectorFEBoundaryFluxLFIntegrator(fnatcoeff));
+   fform->Assemble();
+   fform->SyncAliasMemory(rhs);
+
+   LinearForm *gform(new LinearForm);
+   gform->Update(W_space, rhs.GetBlock(1), 0);
+   gform->AddDomainIntegrator(new DomainLFIntegrator(gcoeff));
+   gform->Assemble();
+   gform->SyncAliasMemory(rhs);
+
+   // 9. Assemble the finite element matrices for the Darcy operator
+   //
+   //                            D = [ M  B^T ]
+   //                                [ B   0  ]
+   //     where:
+   //
+   //     M = \int_\Omega k u_h \cdot v_h d\Omega   u_h, v_h \in R_h
+   //     B   = -\int_\Omega \div u_h q_h d\Omega   u_h \in R_h, q_h \in W_h
+   BilinearForm *mVarf(new BilinearForm(R_space));
+   mVarf->AddDomainIntegrator(new VectorFEMassIntegrator(one));
+   mVarf->Assemble();
+   MixedBilinearForm *bVarf(new MixedBilinearForm(R_space, W_space));
+   bVarf->AddDomainIntegrator(new VectorFEDivergenceIntegrator);
+   bVarf->Assemble();
+   mVarf->Finalize();
+   bVarf->Finalize();
+   BlockOperator darcyOp(block_offsets);
+   TransposeOperator *Bt = NULL;
+   SparseMatrix &M(mVarf->SpMat());
+   SparseMatrix &B(bVarf->SpMat());
+   B *= -1.;
+   Bt = new TransposeOperator(&B);
+   darcyOp.SetBlock(0,0, &M);
+   darcyOp.SetBlock(0,1, Bt);
+   darcyOp.SetBlock(1,0, &B);
+
+   // 10. Construct the operators for preconditioner
+   //
+   //                 P = [ diag(M)         0         ]
+   //                     [  0       B diag(M)^-1 B^T ]
+   //
+   //     Here we use Symmetric Gauss-Seidel to approximate the inverse of the
+   //     pressure Schur Complement
+   SparseMatrix *MinvBt = NULL;
+   Vector Md(mVarf->Height());
+
+   BlockDiagonalPreconditioner darcyPrec(block_offsets);
+   Solver *invM, *invS;
+   SparseMatrix *S = NULL;
+   // SparseMatrix &M(mVarf->SpMat());
+   M.GetDiag(Md);
+   Md.HostReadWrite();
+   // SparseMatrix &B(bVarf->SpMat());
+   MinvBt = Transpose(B);
+   for (int i = 0; i < Md.Size(); i++)
+   {
+      MinvBt->ScaleRow(i, 1./Md(i));
+   }
+   S = Mult(B, *MinvBt);
+   invM = new DSmoother(M);
+#ifndef MFEM_USE_SUITESPARSE
+   invS = new GSSmoother(*S);
+#else
+   invS = new UMFPackSolver(*S);
+#endif
+   invM->iterative_mode = false;
+   invS->iterative_mode = false;
+
+   darcyPrec.SetDiagonalBlock(0, invM);
+   darcyPrec.SetDiagonalBlock(1, invS);
+
+   // 11. Solve the linear system with MINRES.
+   //     Check the norm of the unpreconditioned residual.
+   int maxIter(1000);
+   real_t rtol(1.e-6);
+   real_t atol(1.e-10);
+
+   MINRESSolver solver;
+   solver.SetAbsTol(atol);
+   solver.SetRelTol(rtol);
+   solver.SetMaxIter(maxIter);
+   solver.SetOperator(darcyOp);
+   solver.SetPreconditioner(darcyPrec);
+   solver.SetPrintLevel(1);
+   x = 0.0;
+   solver.Mult(rhs, x);
+
+   if (solver.GetConverged())
+   {
+      std::cout << "MINRES converged in " << solver.GetNumIterations()
+                << " iterations with a residual norm of "
+                << solver.GetFinalNorm() << ".\n";
+   }
+   else
+   {
+      std::cout << "MINRES did not converge in " << solver.GetNumIterations()
+                << " iterations. Residual norm is " << solver.GetFinalNorm()
+                << ".\n";
+   }
+
+   // 12. Create the grid functions u and p. Compute the L2 error norms.
+   GridFunction u, p;
+   u.MakeRef(R_space, x.GetBlock(0), 0);
+   p.MakeRef(W_space, x.GetBlock(1), 0);
+
+   int order_quad = max(2, 2*order_darcy+1);
+   const IntegrationRule *irs[Geometry::NumGeom];
+   for (int i=0; i < Geometry::NumGeom; ++i)
+   {
+      irs[i] = &(IntRules.Get(i, order_quad));
+   }
+
+   real_t err_u  = u.ComputeL2Error(ucoeff, irs);
+   real_t norm_u = ComputeLpNorm(2., ucoeff, mesh, irs);
+   real_t err_p  = p.ComputeL2Error(pcoeff, irs);
+   real_t norm_p = ComputeLpNorm(2., pcoeff, mesh, irs);
+
+   std::cout << "|| u_h - u_ex || / || u_ex || = " << err_u / norm_u << "\n";
+   std::cout << "|| p_h - p_ex || / || p_ex || = " << err_p / norm_p << "\n";
+
+   if (visualization)
+   {
+      char vishost[] = "localhost";
+      int  visport   = 19916;
+      socketstream u_sock(vishost, visport);
+      u_sock.precision(8);
+      u_sock << "solution\n" << mesh << u << "window_title 'Velocity'" << endl;
+      socketstream p_sock(vishost, visport);
+      p_sock.precision(8);
+      p_sock << "solution\n" << mesh << p << "window_title 'Pressure'" << endl;
+   }
+
+
+   // ******Forward Advection-Diffusion solve
+   // 13. Define the DG finite element space on the
+   //    refined mesh of the given polynomial order.
+   DG_FECollection fec(order_ad, dim, BasisType::GaussLobatto);
+   FiniteElementSpace fes(&mesh, &fec);
+   int num_dofs = fes.GetNDofs();
+
+   cout << "Number of unknowns (advection diffusion problem): " << fes.GetVSize()
+        << endl;
+
+   // 14. Set up and assemble the parallel bilinear and linear forms (and the
+   //    parallel hypre matrices) corresponding to the DG discretization. The
+   //    DGTraceIntegrator involves integrals over mesh interior faces.
+   const GridFunction* u_pointer = &u;
+   VectorGridFunctionCoefficient velocity(u_pointer);
+   FunctionCoefficient inflow(inflow_function);
+   ConstantCoefficient diff_coef(d_coef);
+
+   BilinearForm m(&fes);
+   m.AddDomainIntegrator(new MassIntegrator);
+
+   BilinearForm k(&fes);
+   k.AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0));
+   k.AddInteriorFaceIntegrator(new NonconservativeDGTraceIntegrator(velocity,
+                                                                    -1.0));
+   k.AddBdrFaceIntegrator(new NonconservativeDGTraceIntegrator(velocity, -1.0));
+
+   BilinearForm s(&fes);
+   s.AddDomainIntegrator(new DiffusionIntegrator(diff_coef));
+   s.AddInteriorFaceIntegrator(new DGDiffusionIntegrator(diff_coef, sigma, kappa));
+   s.AddBdrFaceIntegrator(new DGDiffusionIntegrator(diff_coef, sigma, kappa));
+
+   LinearForm b(&fes);
+   b.AddBdrFaceIntegrator(new BoundaryFlowIntegrator(inflow, velocity, -1.0));
+   //b.AddBdrFaceIntegrator(new DGDirichletLFIntegrator(U, diff_coef, sigma, kappa));
+
+   int skip_zeros = 0;
+   m.Assemble(skip_zeros);
+   k.Assemble(skip_zeros);
+   s.Assemble(skip_zeros);
+   b.Assemble();
+
+   m.Finalize(skip_zeros);
+   k.Finalize(skip_zeros);
+   s.Finalize(skip_zeros);
+
+   // 15. Define the initial conditions, save the corresponding grid function to
+   //    a file and (optionally) save data in the VisIt format and initialize
+   //    GLVis visualization.
+   FunctionCoefficient theta0(theta0_function);
+   GridFunction theta(&fes);
+   theta.ProjectCoefficient(theta0);
+
+   // Set up visualization, if desired.
+   ParaViewDataCollection *pd_forward = NULL;
+   if (paraview)
+   {
+      pd_forward = new ParaViewDataCollection("darcy-adv-diff-forward", &mesh);
+      pd_forward->SetPrefixPath("ParaView");
+      pd_forward->RegisterField("solution_forward", &theta);
+      pd_forward->SetLevelsOfDetail(order_ad);
+      pd_forward->SetDataFormat(VTKFormat::BINARY);
+      pd_forward->SetHighOrderOutput(true);
+      pd_forward->SetCycle(0);
+      pd_forward->SetTime(0.0);
+      pd_forward->Save();
+   }
+
+   // 16. Define the time-dependent evolution operator describing the ODE
+   //    right-hand side, and perform time-integration (looping over the time
+   //    iterations, ti, with a time-step dt).
+   IMEX_Evolution adv(m, k, s, b);
+
+   real_t t = 0.0;
+   adv.SetTime(t);
+   ode_solver->Init(adv);
+
+   int n_steps = (int)ceil(t_final / dt);
+   double dt_real = t_final / n_steps;
+   // Vector err_vec(n_steps-1);
+
+   std::vector<GridFunction> theta_gf_vector;
+   theta_gf_vector.push_back(theta);
+
+   for (int ti = 0; ti < n_steps; ti++)
+   {
+      ode_solver->Step(theta, t, dt_real);
+      theta_gf_vector.push_back(theta);
+      if (ti % vis_steps == 0 || ti == n_steps -1)
+      {
+         cout << "time step: " << ti << ", time: " << t << endl;
+         if (paraview)
+         {
+            pd_forward->SetCycle(ti);
+            pd_forward->SetTime(t);
+            pd_forward->Save();
+         }
+      }
+   }
+
+   // ******Backward Advection-Diffusion solve
+   // 17. Define the DG finite element space on the
+   //    refined mesh of the given polynomial order.
+   DG_FECollection fec_adjoint(order_ad, dim);
+   FiniteElementSpace fes_adjoint(&mesh, &fec_adjoint);
+
+   // 18. Set up and assemble the parallel bilinear and linear forms (and the
+   //    parallel hypre matrices) corresponding to the DG discretization. The
+   //    DGTraceIntegrator involves integrals over mesh interior faces.
+   ConstantCoefficient zero(0.0);
+   GridFunctionCoefficient theta_coeff(&(theta_gf_vector[n_steps-1]));
+   FunctionCoefficient inflow_adj(inflow_function); //zero for now
+   ConstantCoefficient diff_coef_adj(-d_coef);
+
+   // FunctionCoefficient theta_exact_coeff(theta_exact);
+   BilinearForm m_adj(&fes_adjoint);
+   m_adj.AddDomainIntegrator(new MassIntegrator);
+
+   BilinearForm k_adj(&fes_adjoint);
+   k_adj.AddDomainIntegrator(new ConvectionIntegrator(velocity, -1.0));
+   k_adj.AddInteriorFaceIntegrator(new NonconservativeDGTraceIntegrator(velocity,
+                                                                        -1.0));
+   k_adj.AddBdrFaceIntegrator(new NonconservativeDGTraceIntegrator(velocity,
+                                                                   -1.0));
+
+   BilinearForm s_adj(&fes_adjoint);
+   s_adj.AddDomainIntegrator(new DiffusionIntegrator(diff_coef_adj));
+   s_adj.AddInteriorFaceIntegrator(new DGDiffusionIntegrator(diff_coef_adj, sigma,
+                                                             kappa));
+   s_adj.AddBdrFaceIntegrator(new DGDiffusionIntegrator(diff_coef_adj, sigma,
+                                                        kappa));
+
+   LinearForm b_adj(&fes_adjoint);
+   b_adj.AddDomainIntegrator(new DomainLFIntegrator(theta_coeff));
+   //b.AddBdrFaceIntegrator(new DGDirichletLFIntegrator(zero, diff_coef, sigma, kappa));
+
+   //int skip_zeros = 0;
+   m_adj.Assemble(skip_zeros);
+   m_adj.Finalize(skip_zeros);
+   k_adj.Assemble(skip_zeros);
+   k_adj.Finalize(skip_zeros);
+   s_adj.Assemble(skip_zeros);
+   s_adj.Finalize(skip_zeros);
+   b_adj.Assemble();
+
+   // 19. Define the initial conditions, save the corresponding grid function to
+   //    a file and (optionally) save data in the VisIt format and initialize
+   //    GLVis visualization.
+   GridFunction lam(&fes_adjoint);
+   lam.ProjectCoefficient(zero);
+   ParaViewDataCollection *pd_backward = NULL;
+   if (paraview)
+   {
+      pd_backward = new ParaViewDataCollection("darcy-adv-diff-backward", &mesh);
+      pd_backward->SetPrefixPath("ParaView");
+      pd_backward->RegisterField("solution-backward", &lam);
+      pd_backward->SetLevelsOfDetail(order_ad);
+      pd_backward->SetDataFormat(VTKFormat::BINARY);
+      pd_backward->SetHighOrderOutput(true);
+      pd_backward->SetCycle(0);
+      pd_backward->SetTime(t_final);
+      pd_backward->Save();
+   }
+
+   // 20. Define the time-dependent evolution operator describing the ODE
+   //    right-hand side, and perform time-integration (looping over the time
+   //    iterations, ti, with a time-step dt).
+   IMEX_Evolution adv_adj(m_adj, k_adj, s_adj, b_adj);
+
+   real_t t_adj = t_final;
+   adv_adj.SetTime(t_adj);
+   ode_solver_adj->Init(adv_adj);
+
+   // int n_steps = (int)ceil(t_final / dt);
+   double dt_real_adj = -dt;
+   std::cout << "dt back = " << dt_real_adj << std::endl;
+   //Vector err_vec(n_steps-1);
+
+   for (int ti = 0; ti < n_steps; ti++)
+   {
+      ode_solver_adj->Step(lam, t_adj, dt_real_adj);
+      Vector lam_vals(num_dofs);
+      Vector theta_values(num_dofs);
+      const GridFunction* theta_gf = theta_coeff.GetGridFunction();
+      theta_gf->GetTrueDofs(theta_values);
+      lam.GetTrueDofs(lam_vals);
+      theta_coeff = *(new GridFunctionCoefficient(&(theta_gf_vector[n_steps - ti -
+                                                                            1])));
+      b_adj = *(new LinearForm(&fes_adjoint));
+      b_adj.AddDomainIntegrator(new DomainLFIntegrator(theta_coeff));
+      b_adj.Assemble();
+      if (ti % vis_steps == 0 || ti == n_steps - 1)
+      {
+         cout << "time step: " << ti << ", time: " << t_adj << endl;
+         if (paraview)
+         {
+            pd_backward->SetCycle(ti);
+            pd_backward->SetTime(t_adj);
+            pd_backward->Save();
+         }
+      }
+   }
+
+   // 21. Free the used memory.
+   // delete &ode_solver;
+   // delete &adv;
+   // delete &adv_adj;
+   delete fform;
+   delete gform;
+   delete invM;
+   delete invS;
+   delete S;
+   delete Bt;
+   delete MinvBt;
+   delete mVarf;
+   delete bVarf;
+   delete W_space;
+   delete R_space;
+   delete l2_coll;
+   delete hdiv_coll;
+   // delete &b_adj;
+   // delete &theta_coeff;
+
+   return 0;
+}
+
+void uFun_ex(const Vector & x, Vector & u)
+{
+   real_t xi(x(0));
+   real_t yi(x(1));
+   real_t zi(0.0);
+   if (x.Size() == 3)
+   {
+      zi = x(2);
+   }
+
+   u(0) = - exp(xi)*sin(yi)*cos(zi);
+   u(1) = - exp(xi)*cos(yi)*cos(zi);
+
+   if (x.Size() == 3)
+   {
+      u(2) = exp(xi)*sin(yi)*sin(zi);
+   }
+}
+
+// Change if needed
+real_t pFun_ex(const Vector & x)
+{
+   real_t xi(x(0));
+   real_t yi(x(1));
+   real_t zi(0.0);
+
+   if (x.Size() == 3)
+   {
+      zi = x(2);
+   }
+
+   return exp(xi)*sin(yi)*cos(zi);
+}
+
+void fFun(const Vector & x, Vector & f)
+{
+   f = 0.0;
+}
+
+real_t gFun(const Vector & x)
+{
+   if (x.Size() == 3)
+   {
+      return -pFun_ex(x);
+   }
+   else
+   {
+      return 0;
+   }
+}
+
+real_t f_natural(const Vector & x)
+{
+   return (-pFun_ex(x));
+}
+
+// Implementation of class IMEX_Evolution
+IMEX_Evolution::IMEX_Evolution(BilinearForm &M_, BilinearForm &K_,
+                               BilinearForm &S_, const Vector &b_)
+   : SplitTimeDependentOperator(M_.FESpace()->GetTrueVSize()),
+     M(M_), K(K_), S(S_), b(b_), z(height)
+{
+   Array<int> ess_tdof_list;
+   if (M.GetAssemblyLevel() == AssemblyLevel::LEGACY)
+   {
+      M_prec = make_unique<DSmoother>(M.SpMat());
+      M_solver.SetOperator(M.SpMat());
+      dg_solver = make_unique<DG_Solver>(M.SpMat(), K.SpMat(), S.SpMat(),
+                                         *M.FESpace());
+   }
+   else
+   {
+      M_prec = make_unique<OperatorJacobiSmoother>(M, ess_tdof_list);
+      M_solver.SetOperator(M);
+      dg_solver = NULL;
+   }
+   M_solver.SetPreconditioner(*M_prec);
+   M_solver.iterative_mode = false;
+   M_solver.SetRelTol(1e-9);
+   M_solver.SetAbsTol(0.0);
+   M_solver.SetMaxIter(100);
+   M_solver.SetPrintLevel(0);
+}
+
+void IMEX_Evolution::Mult1(const Vector &x, Vector &y) const
+{
+   // Perform the explicit step
+   // y = M^{-1} (K x + b)
+   K.Mult(x, z);
+   z += b;
+   M_solver.Mult(z, y);
+}
+
+void IMEX_Evolution::ImplicitSolve2(const real_t dt, const Vector &x, Vector &k)
+{
+   // Perform the implicit step
+   // solve for k, k = -(M+dt S)^{-1} S x
+   MFEM_VERIFY(dg_solver != NULL,
+               "Implicit time integration is not supported with partial assembly");
+   S.Mult(x, z);
+   z*= -1.0;
+   dg_solver->SetTimeStep(dt);
+   dg_solver->Mult(z, k);
+}
+
+
+// Initial condition
+double theta0_function(const Vector &x)
+{
+   int dim = x.Size();
+   // map to the reference [-1,1] domain
+   Vector X(dim);
+   // for (int i = 0; i < dim; i++)
+   // {
+   //    double center = (bb_min[i] + bb_max[i]) * 0.5;
+   //    X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
+   // }
+
+   double rx = 0.45, ry = 0.25, cx = 0., cy = -0.2, w = 10.;
+   if (dim == 3)
+   {
+      const double s = (1. + 0.25*cos(2*M_PI*x(2)));
+      rx *= s;
+      ry *= s;
+   }
+   return ( erfc(w*(x(0)-cx-rx))*erfc(-w*(x(0)-cx+rx))*erfc(w*(x(1)-cy-ry))*erfc(
+               -w*(x(1)-cy+ry)) )/16;
+}
+
+//forcing term
+real_t forcing_function(const Vector &x, real_t t)
+{
+   int dim = x.Size();
+   //map to the reference [-1,1] domain
+   Vector X(dim);
+   for (int i = 0; i < dim; i++)
+   {
+      double center = (bb_min[i] + bb_max[i]) * 0.5;
+      X(i) = 2 * (x(i) - center) / (bb_max[i] - bb_min[i]);
+   }
+   return 0.0;
+}
+
+// Inflow boundary condition (zero for the problems considered in this example)
+double inflow_function(const Vector &x)
+{
+   return 0.0;
+}

miniapps/mtop/makefile

@@ -12,10 +12,11 @@
 # Use the MFEM build directory
 MFEM_DIR ?= ../..
 MFEM_BUILD_DIR ?= ../..
-MFEM_INSTALL_DIR ?= ../../mfem
 SRC = $(if $(MFEM_DIR:../..=),$(MFEM_DIR)/miniapps/mtop/,)
-CONFIG_MK = $(or $(wildcard $(MFEM_BUILD_DIR)/config/config.mk),\
-   $(wildcard $(MFEM_INSTALL_DIR)/share/mfem/config.mk))
+CONFIG_MK = $(MFEM_BUILD_DIR)/config/config.mk
+# Use the MFEM install directory
+# MFEM_INSTALL_DIR = ../../mfem
+# CONFIG_MK = $(MFEM_INSTALL_DIR)/share/mfem/config.mk
 
 # Include defaults.mk to get XLINKER
 DEFAULTS_MK = $(MFEM_DIR)/config/defaults.mk
@@ -24,11 +25,10 @@ include $(DEFAULTS_MK)
 MFEM_LIB_FILE = mfem_is_not_built
 -include $(CONFIG_MK)
 
-MTOP_COMMON_SRC = mtop_integrators.cpp paramnonlinearform.cpp pparamnonlinearform.cpp
-
+MTOP_COMMON_SRC =  darcy_heat_transfer_ex.cpp
 MTOP_COMMON_OBJ = $(MTOP_COMMON_SRC:.cpp=.o)
 
-SEQ_MINIAPPS = seqheat
+SEQ_MINIAPPS = seqheat darcy_heat_transfer_ex
 PAR_MINIAPPS = parheat
 ifeq ($(MFEM_USE_MPI),NO)
    MINIAPPS = $(SEQ_MINIAPPS)
@@ -76,4 +76,4 @@ clean-build:
 	rm -rf *.dSYM *.TVD.*breakpoints
 
 clean-exec:
-	@rm -rf SeqHeat* ParHeat*
+	@rm -rf SeqHeat* ParHeat* for_adv_diff_solve* darcy_heat_transfer_ex*

miniapps/mtop/mtop_integrators.cpp

@@ -1,390 +0,0 @@
-// 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 "mtop_integrators.hpp"
-
-namespace mfem
-{
-
-real_t ParametricLinearDiffusion::GetElementEnergy(const
-                                                   Array<const FiniteElement *> &el,
-                                                   const Array<const FiniteElement *> &pel,
-                                                   ElementTransformation &Tr,
-                                                   const Array<const Vector *> &elfun,
-                                                   const Array<const Vector *> &pelfun)
-{
-   int dof_u0 = el[0]->GetDof();
-   int dof_r0 = pel[0]->GetDof();
-
-   int dim = el[0]->GetDim();
-   int spaceDim = Tr.GetSpaceDim();
-   if (dim != spaceDim)
-   {
-      mfem::mfem_error("ParametricLinearDiffusion::GetElementEnergy"
-                       " is not defined on manifold meshes");
-   }
-
-   // shape functions
-   Vector shu0(dof_u0);
-   Vector shr0(dof_r0);
-   DenseMatrix dsu0(dof_u0,dim);
-   DenseMatrix B(dof_u0, 4);
-   B=0.0;
-
-   real_t w;
-
-   Vector param(1); param=0.0;
-   Vector uu(4); uu=0.0;
-
-   real_t energy =0.0;
-
-   const IntegrationRule *ir;
-   {
-      int order= 2 * el[0]->GetOrder() + Tr.OrderGrad(el[0])
-                 +pel[0]->GetOrder();
-      ir=&IntRules.Get(Tr.GetGeometryType(),order);
-   }
-
-   for (int i = 0; i < ir->GetNPoints(); i++)
-   {
-      const IntegrationPoint &ip = ir->IntPoint(i);
-      Tr.SetIntPoint(&ip);
-      w=Tr.Weight();
-      w = ip.weight * w;
-
-      el[0]->CalcPhysDShape(Tr,dsu0);
-      el[0]->CalcPhysShape(Tr,shu0);
-      pel[0]->CalcPhysShape(Tr,shr0);
-
-      param[0]=shr0*(*pelfun[0]);
-
-      // set the matrix B
-      for (int jj=0; jj<dim; jj++)
-      {
-         B.SetCol(jj,dsu0.GetColumn(jj));
-      }
-      B.SetCol(3,shu0);
-      B.MultTranspose(*elfun[0],uu);
-      energy=energy+w * qfun.QEnergy(Tr,ip,param,uu);
-   }
-   return energy;
-}
-
-
-void ParametricLinearDiffusion::AssembleElementVector(const
-                                                      Array<const FiniteElement *> &el,
-                                                      const Array<const FiniteElement *> &pel,
-                                                      ElementTransformation &Tr,
-                                                      const Array<const Vector *> &elfun,
-                                                      const Array<const Vector *> &pelfun,
-                                                      const Array<Vector *> &elvec)
-{
-   int dof_u0 = el[0]->GetDof();
-   int dof_r0 = pel[0]->GetDof();
-
-   int dim = el[0]->GetDim();
-
-   elvec[0]->SetSize(dof_u0);
-   *elvec[0]=0.0;
-   int spaceDim = Tr.GetSpaceDim();
-   if (dim != spaceDim)
-   {
-      mfem::mfem_error("ParametricLinearDiffusion::AssembleElementVector"
-                       " is not defined on manifold meshes");
-   }
-
-   // shape functions
-   Vector shu0(dof_u0);
-   Vector shr0(dof_r0);
-   DenseMatrix dsu0(dof_u0,dim);
-   DenseMatrix B(dof_u0, 4);
-   B=0.0;
-
-   real_t w;
-
-   Vector param(1); param=0.0;
-   Vector uu(4); uu=0.0;
-   Vector rr(4);
-   Vector lvec; lvec.SetSize(dof_u0);
-
-   const IntegrationRule *ir = nullptr;
-   int order= 2 * el[0]->GetOrder() + Tr.OrderGrad(el[0])
-              +pel[0]->GetOrder();
-   ir=&IntRules.Get(Tr.GetGeometryType(),order);
-
-
-   for (int i = 0; i < ir->GetNPoints(); i++)
-   {
-      const IntegrationPoint &ip = ir->IntPoint(i);
-      Tr.SetIntPoint(&ip);
-      w=Tr.Weight();
-      w = ip.weight * w;
-
-      el[0]->CalcPhysDShape(Tr,dsu0);
-      el[0]->CalcPhysShape(Tr,shu0);
-      pel[0]->CalcPhysShape(Tr,shr0);
-
-      param[0]=shr0*(*pelfun[0]);
-
-      // set the matrix B
-      for (int jj=0; jj<dim; jj++)
-      {
-         B.SetCol(jj,dsu0.GetColumn(jj));
-      }
-      B.SetCol(3,shu0);
-      B.MultTranspose(*elfun[0],uu);
-      qfun.QResidual(Tr,ip,param, uu, rr);
-
-      B.Mult(rr,lvec);
-      elvec[0]->Add(w,lvec);
-   }
-}
-
-void ParametricLinearDiffusion::AssembleElementGrad(const
-                                                    Array<const FiniteElement *> &el,
-                                                    const Array<const FiniteElement *> &pel,
-                                                    ElementTransformation &Tr,
-                                                    const Array<const Vector *> &elfun,
-                                                    const Array<const Vector *> &pelfun,
-                                                    const Array2D<DenseMatrix *> &elmats)
-{
-   int dof_u0 = el[0]->GetDof();
-   int dof_r0 = pel[0]->GetDof();
-
-   int dim = el[0]->GetDim();
-
-   DenseMatrix* K=elmats(0,0);
-   K->SetSize(dof_u0,dof_u0);
-   (*K)=0.0;
-
-   int spaceDim = Tr.GetSpaceDim();
-   if (dim != spaceDim)
-   {
-      mfem::mfem_error("ParametricLinearDiffusion::AssembleElementGrad"
-                       " is not defined on manifold meshes");
-   }
-
-   // shape functions
-   Vector shu0(dof_u0);
-   Vector shr0(dof_r0);
-   DenseMatrix dsu0(dof_u0,dim);
-   DenseMatrix B(dof_u0, 4);
-   DenseMatrix A(dof_u0, 4);
-   B=0.0;
-   real_t w;
-
-   Vector param(1); param=0.0;
-   Vector uu(4); uu=0.0;
-   DenseMatrix hh(4,4);
-   Vector lvec; lvec.SetSize(dof_u0);
-
-   const IntegrationRule *ir = nullptr;
-   int order= 2 * el[0]->GetOrder() + Tr.OrderGrad(el[0])
-              +pel[0]->GetOrder();
-   ir=&IntRules.Get(Tr.GetGeometryType(),order);
-
-   for (int i = 0; i < ir->GetNPoints(); i++)
-   {
-      const IntegrationPoint &ip = ir->IntPoint(i);
-      Tr.SetIntPoint(&ip);
-      w = Tr.Weight();
-      w = ip.weight * w;
-
-      el[0]->CalcPhysDShape(Tr,dsu0);
-      el[0]->CalcPhysShape(Tr,shu0);
-      pel[0]->CalcPhysShape(Tr,shr0);
-
-      param[0]=shr0*(*pelfun[0]);
-
-      // set the matrix B
-      for (int jj=0; jj<dim; jj++)
-      {
-         B.SetCol(jj,dsu0.GetColumn(jj));
-      }
-      B.SetCol(3,shu0);
-      B.MultTranspose(*elfun[0],uu);
-      qfun.QGradResidual(Tr,ip,param,uu,hh);
-      Mult(B,hh,A);
-      AddMult_a_ABt(w,A,B,*K);
-   }
-}
-
-void ParametricLinearDiffusion::AssemblePrmElementVector(
-   const Array<const FiniteElement *> &el,
-   const Array<const FiniteElement *> &pel,
-   ElementTransformation &Tr,
-   const Array<const Vector *> &elfun,
-   const Array<const Vector *> &alfun,
-   const Array<const Vector *> &pelfun,
-   const Array<Vector *> &elvec)
-{
-   int dof_u0 = el[0]->GetDof();
-   int dof_r0 = pel[0]->GetDof();
-
-   int dim = el[0]->GetDim();
-   Vector& e0 = *(elvec[0]);
-
-   e0.SetSize(dof_r0);
-   e0=0.0;
-
-   int spaceDim = Tr.GetSpaceDim();
-   if (dim != spaceDim)
-   {
-      mfem::mfem_error("ParametricLinearDiffusion::AssemblePrmElementVector"
-                       " is not defined on manifold meshes");
-   }
-
-   // shape functions
-   Vector shu0(dof_u0);
-   Vector shr0(dof_r0);
-   DenseMatrix dsu0(dof_u0,dim);
-   DenseMatrix B(dof_u0, 4);
-   B=0.0;
-
-   real_t w;
-
-   Vector param(1); param=0.0;
-   Vector uu(4); uu=0.0;
-   Vector aa(4); aa=0.0;
-   Vector rr(1);
-   Vector lvec0; lvec0.SetSize(dof_r0);
-
-   const IntegrationRule *ir;
-   {
-      int order= 2 * el[0]->GetOrder() + Tr.OrderGrad(el[0])
-                 +pel[0]->GetOrder();
-      ir=&IntRules.Get(Tr.GetGeometryType(),order);
-   }
-
-   for (int i = 0; i < ir->GetNPoints(); i++)
-   {
-      const IntegrationPoint &ip = ir->IntPoint(i);
-      Tr.SetIntPoint(&ip);
-      w=Tr.Weight();
-      w = ip.weight * w;
-
-      el[0]->CalcPhysDShape(Tr,dsu0);
-      el[0]->CalcPhysShape(Tr,shu0);
-      pel[0]->CalcPhysShape(Tr,shr0);
-
-      param[0]=shr0*(*pelfun[0]);
-
-      // set the matrix B
-      for (int jj=0; jj<dim; jj++)
-      {
-         B.SetCol(jj,dsu0.GetColumn(jj));
-      }
-      B.SetCol(3,shu0);
-      B.MultTranspose(*elfun[0],uu);
-      B.MultTranspose(*alfun[0],aa);
-
-      qfun.AQResidual(Tr, ip, param, uu, aa, rr);
-
-      lvec0=shr0;
-      lvec0*=rr[0];
-
-      e0.Add(w,lvec0);
-   }
-}
-
-real_t DiffusionObjIntegrator::GetElementEnergy(const
-                                                Array<const FiniteElement *> &el,
-                                                ElementTransformation &Tr,
-                                                const Array<const Vector *> &elfun)
-{
-   int dof_u0 = el[0]->GetDof();
-   int dim = el[0]->GetDim();
-   int spaceDim = Tr.GetSpaceDim();
-   if (dim != spaceDim)
-   {
-      mfem::mfem_error("DiffusionObjIntegrator::GetElementEnergy"
-                       " is not defined on manifold meshes");
-   }
-
-   // shape functions
-   Vector shu0(dof_u0);
-
-   real_t w;
-   real_t val;
-
-   real_t energy = 0.0;
-
-   const IntegrationRule *ir;
-   {
-      int order= 2 * el[0]->GetOrder() + Tr.OrderGrad(el[0]);
-      ir=&IntRules.Get(Tr.GetGeometryType(),order);
-   }
-
-   for (int i = 0; i < ir->GetNPoints(); i++)
-   {
-      const IntegrationPoint &ip = ir->IntPoint(i);
-      Tr.SetIntPoint(&ip);
-      w=Tr.Weight();
-
-      w = ip.weight * w;
-
-      el[0]->CalcPhysShape(Tr,shu0);
-
-      val=shu0*(*elfun[0]);
-      energy=energy + w * val * val;
-   }
-   return 0.5*energy;
-}
-
-void DiffusionObjIntegrator::AssembleElementVector(const
-                                                   Array<const FiniteElement *> &el,
-                                                   ElementTransformation &Tr,
-                                                   const Array<const Vector *> &elfun,
-                                                   const Array<Vector *> &elvec)
-{
-   int dof_u0 = el[0]->GetDof();
-   int dim = el[0]->GetDim();
-   int spaceDim = Tr.GetSpaceDim();
-
-   elvec[0]->SetSize(dof_u0);
-   *elvec[0]=0.0;
-
-   if (dim != spaceDim)
-   {
-      mfem::mfem_error("DiffusionObjIntegrator::GetElementEnergy"
-                       " is not defined on manifold meshes");
-   }
-
-   // shape functions
-   Vector shu0(dof_u0);
-
-   real_t w;
-   real_t val;
-
-   const IntegrationRule *ir;
-   {
-      int order= 2 * el[0]->GetOrder() + Tr.OrderGrad(el[0]);
-      ir=&IntRules.Get(Tr.GetGeometryType(),order);
-   }
-
-   for (int i = 0; i < ir->GetNPoints(); i++)
-   {
-      const IntegrationPoint &ip = ir->IntPoint(i);
-      Tr.SetIntPoint(&ip);
-      w=Tr.Weight();
-
-      w = ip.weight * w;
-
-      el[0]->CalcPhysShape(Tr,shu0);
-
-      val=shu0*(*elfun[0]);
-
-      elvec[0]->Add(w*val,shu0);
-   }
-}
-
-
-} // end mfem namespace

miniapps/mtop/mtop_integrators.hpp

@@ -1,233 +0,0 @@
-// 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 MTOPINTEGRATORS_HPP
-#define MTOPINTEGRATORS_HPP
-
-#include "mfem.hpp"
-#include "paramnonlinearform.hpp"
-
-#include <map>
-
-namespace mfem
-{
-
-/// Base class for representing function at integration points.
-class BaseQFunction
-{
-public:
-   virtual ~BaseQFunction() {}
-
-   /// Returns a user defined string identifying the function.
-   virtual std::string GetType()=0;
-
-   // Returns the energy at an integration point.
-   virtual
-   real_t QEnergy(ElementTransformation &T, const IntegrationPoint &ip,
-                  mfem::Vector &dd, mfem::Vector &uu)
-   {
-      return 0.0;
-   }
-
-   // Returns the residual at an integration point.
-   virtual
-   void QResidual(ElementTransformation &T, const IntegrationPoint &ip,
-                  mfem::Vector &dd, mfem::Vector &uu, mfem::Vector &rr)=0;
-
-   /// Returns the gradient of the residual at a integration point.
-   virtual
-   void QGradResidual(ElementTransformation &T, const IntegrationPoint &ip,
-                      mfem::Vector &dd, mfem::Vector &uu, mfem::DenseMatrix &hh)=0;
-
-   /// Returns the gradient of the residual with respect to the design
-   /// parameters, multiplied by the adjoint.
-   virtual
-   void AQResidual(ElementTransformation &T, const IntegrationPoint &ip,
-                   mfem::Vector &dd, mfem::Vector &uu,
-                   mfem::Vector &aa, mfem::Vector &rr)=0;
-
-};
-
-/* QLinearDiffusion implements methods for computing the energy, the residual,
- * gradient of the residual and the product of the adjoint fields with the
- * derivative of the residual with respect to the parameters. All computations
- * are performed at a integration point. Therefore the vectors (vv,uu,aa,rr ..)
- * hold the fields' values and the fields' derivatives at the integration
- * point. For example for a single scalar parametric field representing the
- * density in topology optimization the vector dd will have size one and the
- * element will be the density at the integration point. The map between state
- * and parameter is not fixed and depends on the implementation of the QFunction
- * class. */
-class QLinearDiffusion:public BaseQFunction
-{
-public:
-   QLinearDiffusion(mfem::Coefficient& diffco, mfem::Coefficient& hsrco,
-                    real_t pp=1.0, real_t minrho=1e-7, real_t betac=4.0, real_t etac=0.5):
-      diff(diffco),load(hsrco), powerc(pp), rhomin(minrho), beta(betac), eta(etac)
-   {
-
-   }
-
-   std::string GetType() override
-   {
-      return "QLinearDiffusion";
-   }
-
-   real_t QEnergy(ElementTransformation &T, const IntegrationPoint &ip,
-                  Vector &dd, Vector &uu) override
-   {
-      // dd[0] - density
-      // uu[0] - grad_x
-      // uu[1] - grad_y
-      // uu[2] - grad_z
-      // uu[3] - temperature/scalar field
-
-      real_t di=diff.Eval(T,ip);
-      real_t ll=load.Eval(T,ip);
-      // Computes the physical density using projection.
-      real_t rz=0.5+0.5*std::tanh(beta*(dd[0]-eta)); //projection
-      // Computes the diffusion coefficient at the integration point.
-      real_t fd=di*(std::pow(rz,powerc)+rhomin);
-      // Computes the sum of the energy and the product of the temperature and
-      // the external input at the integration point.
-      real_t rez = 0.5*(uu[0]*uu[0]+uu[1]*uu[1]+uu[2]*uu[2])*fd-uu[3]*ll;
-      return rez;
-   }
-
-   /// Returns the derivative of QEnergy with respect to the state vector uu.
-   void QResidual(ElementTransformation &T, const IntegrationPoint &ip,
-                  Vector &dd, Vector &uu, Vector &rr) override
-   {
-      real_t di=diff.Eval(T,ip);
-      real_t ll=load.Eval(T,ip);
-      real_t rz=0.5+0.5*std::tanh(beta*(dd[0]-eta));
-      real_t fd=di*(std::pow(rz,powerc)+rhomin);
-
-      rr[0]=uu[0]*fd;
-      rr[1]=uu[1]*fd;
-      rr[2]=uu[2]*fd;
-      rr[3]=-ll;
-   }
-
-
-   // Returns the derivative, with respect to the density, of the product of
-   // the adjoint field with the residual at the integration point ip.
-   void AQResidual(ElementTransformation &T, const IntegrationPoint &ip,
-                   Vector &dd, Vector &uu, Vector &aa, Vector &rr) override
-   {
-      real_t di=diff.Eval(T,ip);
-      real_t tt=std::tanh(beta*(dd[0]-eta));
-      real_t rz=0.5+0.5*tt;
-      real_t fd=di*powerc*std::pow(rz,powerc-1.0)*0.5*(1.0-tt*tt)*beta;
-
-      rr[0] = -(aa[0]*uu[0]+aa[1]*uu[1]+aa[2]*uu[2])*fd;
-   }
-
-   // Returns the gradient of the residual with respect to the state vector at
-   // the integration point ip.
-   void QGradResidual(ElementTransformation &T, const IntegrationPoint &ip,
-                      Vector &dd, Vector &uu, DenseMatrix &hh) override
-   {
-      real_t di=diff.Eval(T,ip);
-      real_t tt=std::tanh(beta*(dd[0]-eta));
-      real_t rz=0.5+0.5*tt;
-      real_t fd=di*(std::pow(rz,powerc)+rhomin);
-      hh=0.0;
-
-      hh(0,0)=fd;
-      hh(1,1)=fd;
-      hh(2,2)=fd;
-      hh(3,3)=0.0;
-   }
-
-private:
-   mfem::Coefficient& diff; //diffusion coefficient
-   mfem::Coefficient& load; //load coefficient
-   real_t powerc; //penalization coefficient
-   real_t rhomin; //lower bound for the density
-   real_t beta;   //controls the sharpness of the projection
-   real_t eta;    //projection threshold for tanh
-};
-
-/// Provides implementation of an integrator for linear diffusion with
-/// parametrization provided by a density field. The setup is standard for
-/// topology optimization problems.
-class ParametricLinearDiffusion: public ParametricBNLFormIntegrator
-{
-public:
-   ParametricLinearDiffusion(BaseQFunction& qfunm): qfun(qfunm)
-   {
-
-   }
-
-   /// Computes the local energy.
-   real_t GetElementEnergy(const Array<const FiniteElement *> &el,
-                           const Array<const FiniteElement *> &pel,
-                           ElementTransformation &Tr,
-                           const Array<const Vector *> &elfun,
-                           const Array<const Vector *> &pelfun) override;
-
-   /// Computes the element's residual.
-   void AssembleElementVector(const Array<const FiniteElement *> &el,
-                              const Array<const FiniteElement *> &pel,
-                              ElementTransformation &Tr,
-                              const Array<const Vector *> &elfun,
-                              const Array<const Vector *> &pelfun,
-                              const Array<Vector *> &elvec) override;
-
-   /// Computes the stiffness/tangent matrix.
-   void AssembleElementGrad(const Array<const FiniteElement *> &el,
-                            const Array<const FiniteElement *> &pel,
-                            ElementTransformation &Tr,
-                            const Array<const Vector *> &elfun,
-                            const Array<const Vector *> &pelfun,
-                            const Array2D<DenseMatrix *> &elmats) override;
-
-   /// Computes the product of the adjoint solution and the derivative of the
-   /// residual with respect to the parametric fields.
-   void AssemblePrmElementVector(const Array<const FiniteElement *> &el,
-                                 const Array<const FiniteElement *> &pel,
-                                 ElementTransformation &Tr,
-                                 const Array<const Vector *> &elfun,
-                                 const Array<const Vector *> &alfun,
-                                 const Array<const Vector *> &pelfun,
-                                 const Array<Vector *> &elvec) override;
-private:
-   BaseQFunction& qfun;
-};
-
-
-/// Computes an example of nonlinear objective
-/// $\int \rm{field}*\rm{field}*\rm{weight})\rm{d}\Omega_e$.
-class DiffusionObjIntegrator:public BlockNonlinearFormIntegrator
-{
-public:
-
-   DiffusionObjIntegrator()
-   {
-
-   }
-
-   /// Returns the objective contribution at element level.
-   real_t GetElementEnergy(const Array<const FiniteElement *> &el,
-                           ElementTransformation &Tr,
-                           const Array<const Vector *> &elfun) override;
-
-   /// Returns the gradient of the objective contribution at element level.
-   void AssembleElementVector(const Array<const FiniteElement *> &el,
-                              ElementTransformation &Tr,
-                              const Array<const Vector *> &elfun,
-                              const Array<Vector *> &elvec) override;
-};
-
-}
-
-#endif

miniapps/mtop/paramnonlinearform.hpp

@@ -1,300 +0,0 @@
-// 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_PRMNONLINEARFORM
-#define MFEM_PRMNONLINEARFORM
-
-#include "mfem.hpp"
-
-namespace mfem
-{
-
-/** The abstract base class ParametricBNLFormIntegrator is a generalization of
-    the BlockNonlinearFormIntegrator class suitable for block state and
-    parameter vectors. */
-class ParametricBNLFormIntegrator
-{
-public:
-   /// Compute the local energy
-   virtual real_t GetElementEnergy(const Array<const FiniteElement *>&el,
-                                   const Array<const FiniteElement *>&pel,
-                                   ElementTransformation &Tr,
-                                   const Array<const Vector *>&elfun,
-                                   const Array<const Vector *>&pelfun);
-
-   /// Perform the local action of the BlockNonlinearFormIntegrator
-   virtual void AssembleElementVector(const Array<const FiniteElement *> &el,
-                                      const Array<const FiniteElement *>&pel,
-                                      ElementTransformation &Tr,
-                                      const Array<const Vector *> &elfun,
-                                      const Array<const Vector *>&pelfun,
-                                      const Array<Vector *> &elvec);
-
-   /// Perform the local action of the BlockNonlinearFormIntegrator on element
-   /// faces
-   virtual void AssembleFaceVector(const Array<const FiniteElement *> &el1,
-                                   const Array<const FiniteElement *> &el2,
-                                   const Array<const FiniteElement *> &pel1,
-                                   const Array<const FiniteElement *> &pel2,
-                                   FaceElementTransformations &Tr,
-                                   const Array<const Vector *> &elfun,
-                                   const Array<const Vector *>&pelfun,
-                                   const Array<Vector *> &elvect);
-
-   /// Perform the local action on the parameters of the BNLFormIntegrator
-   virtual void AssemblePrmElementVector(const Array<const FiniteElement *> &el,
-                                         const Array<const FiniteElement *>&pel,
-                                         ElementTransformation &Tr,
-                                         const Array<const Vector *> &elfun,
-                                         const Array<const Vector *> &alfun,
-                                         const Array<const Vector *>&pelfun,
-                                         const Array<Vector *> &pelvec);
-
-   /// Perform the local action on the parameters of the BNLFormIntegrator on
-   /// faces
-   virtual void AssemblePrmFaceVector(const Array<const FiniteElement *> &el1,
-                                      const Array<const FiniteElement *> &el2,
-                                      const Array<const FiniteElement *> &pel1,
-                                      const Array<const FiniteElement *> &pel2,
-                                      FaceElementTransformations &Tr,
-                                      const Array<const Vector *> &elfun,
-                                      const Array<const Vector *> &alfun,
-                                      const Array<const Vector *>&pelfun,
-                                      const Array<Vector *> &pelvect);
-
-   /// Assemble the local gradient matrix
-   virtual void AssembleElementGrad(const Array<const FiniteElement*> &el,
-                                    const Array<const FiniteElement *>&pel,
-                                    ElementTransformation &Tr,
-                                    const Array<const Vector *> &elfun,
-                                    const Array<const Vector *>&pelfun,
-                                    const Array2D<DenseMatrix *> &elmats);
-
-   /// Assemble the local gradient matrix on faces of the elements
-   virtual void AssembleFaceGrad(const Array<const FiniteElement *>&el1,
-                                 const Array<const FiniteElement *>&el2,
-                                 const Array<const FiniteElement *> &pel1,
-                                 const Array<const FiniteElement *> &pel2,
-                                 FaceElementTransformations &Tr,
-                                 const Array<const Vector *> &elfun,
-                                 const Array<const Vector *>&pelfun,
-                                 const Array2D<DenseMatrix *> &elmats);
-
-
-   virtual ~ParametricBNLFormIntegrator() { }
-};
-
-
-/** @brief A class representing a general parametric block nonlinear operator
-    defined on the Cartesian product of multiple FiniteElementSpace%s. */
-class ParametricBNLForm : public Operator
-{
-protected:
-   /// FE spaces on which the form lives.
-   Array<FiniteElementSpace*> fes;
-
-   /// FE spaces for the parametric fields
-   Array<FiniteElementSpace*> paramfes;
-
-   int paramheight;
-   int paramwidth;
-
-   /// Set of Domain Integrators to be assembled (added).
-   Array<ParametricBNLFormIntegrator*> dnfi;
-
-   /// Set of interior face Integrators to be assembled (added).
-   Array<ParametricBNLFormIntegrator*> fnfi;
-
-   /// Set of Boundary Face Integrators to be assembled (added).
-   Array<ParametricBNLFormIntegrator*> bfnfi;
-   Array<Array<int>*>           bfnfi_marker;
-
-   /** Auxiliary block-vectors for wrapping input and output vectors or holding
-       GridFunction-like block-vector data (e.g. in parallel). */
-   mutable BlockVector xs, ys;
-   mutable BlockVector prmxs, prmys;
-
-   /** Auxiliary block-vectors for holding GridFunction-like block-vector data
-       (e.g. in parallel). */
-   mutable BlockVector xsv;
-
-   /** Auxiliary block-vectors for holding GridFunction-like block-vector data
-       for the parameter fields (e.g. in parallel). */
-   mutable BlockVector xdv;
-   /** Auxiliary block-vectors for holding GridFunction-like block-vector data
-       for the adjoint fields (e.g. in parallel). */
-   mutable BlockVector adv;
-
-   mutable Array2D<SparseMatrix*> Grads, cGrads;
-   mutable BlockOperator *BlockGrad;
-
-   // A list of the offsets
-   Array<int> block_offsets;
-   Array<int> block_trueOffsets;
-   // A list with the offsets for the parametric fields
-   Array<int> paramblock_offsets;
-   Array<int> paramblock_trueOffsets;
-
-   // Array of Arrays of tdofs for each space in 'fes'
-   Array<Array<int> *> ess_tdofs;
-
-   // Array of Arrays of tdofs for each space in 'paramfes'
-   Array<Array<int> *> paramess_tdofs;
-
-   /// Array of pointers to the prolongation matrix of fes, may be NULL
-   Array<const Operator *> P;
-
-   /// Array of pointers to the prolongation matrix of paramfes, may be NULL
-   Array<const Operator *> Pparam;
-
-   /// Array of results of dynamic-casting P to SparseMatrix pointer
-   Array<const SparseMatrix *> cP;
-
-   /// Array of results of dynamic-casting Pparam to SparseMatrix pointer
-   Array<const SparseMatrix *> cPparam;
-
-   /// Indicator if the Operator is part of a parallel run
-   bool is_serial = true;
-
-   /// Indicator if the Operator needs prolongation on assembly
-   bool needs_prolongation = false;
-
-   /// Indicator if the Operator needs prolongation on assembly
-   bool prmneeds_prolongation = false;
-
-   mutable BlockVector aux1, aux2;
-
-   mutable BlockVector prmaux1, prmaux2;
-
-   const BlockVector &Prolongate(const BlockVector &bx) const;
-
-   const BlockVector &ParamProlongate(const BlockVector &bx) const;
-
-   real_t GetEnergyBlocked(const BlockVector &bx, const BlockVector &dx) const;
-
-
-   /// Specialized version of Mult() for BlockVector%s
-   /// Block L-Vector to Block L-Vector
-   void MultBlocked(const BlockVector &bx, const BlockVector &dx,
-                    BlockVector &by) const;
-
-   /// Specialized version of Mult() for BlockVector%s
-   /// Block L-Vector to Block L-Vector
-   /// bx - state vector, ax - adjoint vector, dx - parametric fields
-   /// dy = ax' d(residual(bx))/d(dx)
-   void MultParamBlocked(const BlockVector &bx, const BlockVector & ax,
-                         const BlockVector &dx, BlockVector &dy) const;
-
-
-   /// Specialized version of GetGradient() for BlockVector
-   void ComputeGradientBlocked(const BlockVector &bx, const BlockVector &dx) const;
-
-public:
-   /// Construct an empty BlockNonlinearForm. Initialize with SetSpaces().
-   ParametricBNLForm();
-
-   /// Construct a BlockNonlinearForm on the given set of FiniteElementSpace%s.
-   ParametricBNLForm(Array<FiniteElementSpace *> &statef,
-                     Array<FiniteElementSpace *> &paramf);
-
-   /// Return the @a k-th FE space of the ParametricBNLForm.
-   FiniteElementSpace *FESpace(int k) { return fes[k]; }
-
-   /// Return the @a k-th parametric FE space of the ParametricBNLForm.
-   FiniteElementSpace *ParamFESpace(int k) { return paramfes[k]; }
-
-
-   /// Return the @a k-th FE space of the BlockNonlinearForm (const version).
-   const FiniteElementSpace *FESpace(int k) const { return fes[k]; }
-
-   /// Return the @a k-th parametric FE space of the BlockNonlinearForm (const
-   /// version).
-   const FiniteElementSpace *ParamFESpace(int k) const { return paramfes[k]; }
-
-   /// Return the integrators
-   Array<ParametricBNLFormIntegrator*>& GetDNFI() { return dnfi;}
-
-
-   /// (Re)initialize the ParametricBNLForm.
-   /** After a call to SetSpaces(), the essential b.c. must be set again. */
-   void SetSpaces(Array<FiniteElementSpace *> &statef,
-                  Array<FiniteElementSpace *> &paramf);
-
-   /// Return the regular dof offsets.
-   const Array<int> &GetBlockOffsets() const { return block_offsets; }
-
-   /// Return the true-dof offsets.
-   const Array<int> &GetBlockTrueOffsets() const { return block_trueOffsets; }
-
-   /// Return the regular dof offsets for the parameters.
-   const Array<int> &ParamGetBlockOffsets() const { return paramblock_offsets; }
-
-   /// Return the true-dof offsets for the parameters.
-   const Array<int> &ParamGetBlockTrueOffsets() const { return paramblock_trueOffsets; }
-
-   /// Adds new Domain Integrator.
-   void AddDomainIntegrator(ParametricBNLFormIntegrator *nlfi)
-   { dnfi.Append(nlfi); }
-
-   /// Adds new Interior Face Integrator.
-   void AddInteriorFaceIntegrator(ParametricBNLFormIntegrator *nlfi)
-   { fnfi.Append(nlfi); }
-
-   /// Adds new Boundary Face Integrator.
-   void AddBdrFaceIntegrator(ParametricBNLFormIntegrator *nlfi)
-   { bfnfi.Append(nlfi); bfnfi_marker.Append(NULL); }
-
-   /** @brief Adds new Boundary Face Integrator, restricted to specific boundary
-       attributes. */
-   void AddBdrFaceIntegrator(ParametricBNLFormIntegrator *nlfi,
-                             Array<int> &bdr_marker);
-
-   /// Set the essential boundary conditions.
-   virtual void SetEssentialBC(const Array<Array<int> *>&bdr_attr_is_ess,
-                               Array<Vector *> &rhs);
-
-   /// Set the essential boundary conditions on the parametric fields.
-   virtual void SetParamEssentialBC(const Array<Array<int> *>&bdr_attr_is_ess,
-                                    Array<Vector *> &rhs);
-
-
-   /// Computes the energy for a state vector x.
-   virtual real_t GetEnergy(const Vector &x) const;
-
-   /// Method is only called in serial, the parallel version calls MultBlocked
-   /// directly.
-   void Mult(const Vector &x, Vector &y) const override;
-
-   /// Method is only called in serial, the parallel version calls MultBlocked
-   /// directly.
-   virtual void ParamMult(const Vector &x, Vector &y) const;
-
-   /// Method is only called in serial, the parallel version calls
-   /// GetGradientBlocked directly.
-   BlockOperator &GetGradient(const Vector &x) const override;
-
-   /// Set the state fields
-   virtual void SetStateFields(const Vector &xv) const;
-
-   /// Set the adjoint fields
-   virtual void SetAdjointFields(const Vector &av) const;
-
-   /// Set the parameters/design fields
-   virtual void SetParamFields(const Vector &dv) const;
-
-   /// Destructor.
-   virtual ~ParametricBNLForm();
-
-};
-
-}
-
-#endif

miniapps/mtop/parheat.cpp

@@ -1,354 +0,0 @@
-// 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.
-//
-//     ----------------------------------------------------------------
-//     ParHeat Miniapp: Gradients of PDE constrained objective function
-//     ----------------------------------------------------------------
-//                           (Parallel Version)
-//
-// The following example computes the gradients of a specified objective
-// function with respect to parametric fields. The objective function is having
-// the following form f(u(\rho)) where u(\rho) is a solution of a specific state
-// problem (in the example that is the diffusion equation), and \rho is a
-// parametric field discretized by finite elements. The parametric field (also
-// called density in topology optimization) controls the coefficients of the
-// state equation. For the considered case, the density controls the diffusion
-// coefficient within the computational domain.
-//
-// For more information, the users are referred to:
-//
-//    Hinze, M.; Pinnau, R.; Ulbrich, M. & Ulbrich, S.
-//    Optimization with PDE Constraints
-//    Springer Netherlands, 2009
-//
-//    Bendsøe, M. P. & Sigmund, O.
-//    Topology Optimization - Theory, Methods and Applications
-//    Springer Verlag, Berlin Heidelberg, 2003
-//
-// Compile with: make parheat
-//
-// Sample runs:
-//
-//    mpirun -np 4 parheat --visualization
-//    mpirun -np 4 parheat --visualization -m ../../data/beam-quad.mesh
-
-#include "mfem.hpp"
-#include <fstream>
-#include <iostream>
-
-#include "pparamnonlinearform.hpp"
-#include "mtop_integrators.hpp"
-
-using namespace mfem;
-
-int main(int argc, char *argv[])
-{
-   // 1. Initialize MPI and HYPRE.
-   mfem::Mpi::Init(argc, argv);
-   int myrank = mfem::Mpi::WorldRank();
-   mfem::Hypre::Init();
-
-   // Parse command-line options.
-   const char *mesh_file = "../../data/star.mesh";
-   int order = 1;
-   bool static_cond = false;
-   int ser_ref_levels = 1;
-   int par_ref_levels = 1;
-   real_t newton_rel_tol = 1e-7;
-   real_t newton_abs_tol = 1e-12;
-   int newton_iter = 10;
-   int print_level = 1;
-   bool visualization = false;
-
-   mfem::OptionsParser args(argc, argv);
-   args.AddOption(&mesh_file, "-m", "--mesh",
-                  "Mesh file to use.");
-   args.AddOption(&ser_ref_levels,
-                  "-rs",
-                  "--refine-serial",
-                  "Number of times to refine the mesh uniformly in serial.");
-   args.AddOption(&par_ref_levels,
-                  "-rp",
-                  "--refine-parallel",
-                  "Number of times to refine the mesh uniformly in parallel.");
-   args.AddOption(&order, "-o", "--order",
-                  "Finite element order (polynomial degree) or -1 for"
-                  " isoparametric space.");
-   args.AddOption(&visualization,
-                  "-vis",
-                  "--visualization",
-                  "-no-vis",
-                  "--no-visualization",
-                  "Enable or disable GLVis visualization.");
-   args.AddOption(&static_cond, "-sc", "--static-condensation", "-no-sc",
-                  "--no-static-condensation", "Enable static condensation.");
-   args.AddOption(&newton_rel_tol,
-                  "-rel",
-                  "--relative-tolerance",
-                  "Relative tolerance for the Newton solve.");
-   args.AddOption(&newton_abs_tol,
-                  "-abs",
-                  "--absolute-tolerance",
-                  "Absolute tolerance for the Newton solve.");
-   args.AddOption(&newton_iter,
-                  "-it",
-                  "--newton-iterations",
-                  "Maximum iterations for the Newton solve.");
-   args.Parse();
-   if (!args.Good())
-   {
-      if (myrank == 0)
-      {
-         args.PrintUsage(std::cout);
-      }
-      return 1;
-   }
-
-   if (myrank == 0)
-   {
-      args.PrintOptions(std::cout);
-   }
-
-   // Read the (serial) mesh from the given mesh file on all processors.  We
-   // can handle triangular, quadrilateral, tetrahedral, hexahedral, surface
-   // and volume meshes with the same code.
-   mfem::Mesh mesh(mesh_file, 1, 1);
-   int dim = mesh.Dimension();
-
-   // Refine the serial mesh on all processors to increase the resolution. In
-   // this example we do 'ref_levels' of uniform refinement. We choose
-   // 'ref_levels' to be the largest number that gives a final mesh with no
-   // more than 10,000 elements.
-   {
-      int ref_levels =
-         (int)floor(log(10000./mesh.GetNE())/log(2.)/dim);
-      for (int l = 0; l < ref_levels; l++)
-      {
-         mesh.UniformRefinement();
-      }
-   }
-
-   // Define a parallel mesh by a partitioning of the serial mesh. Refine
-   // this mesh further in parallel to increase the resolution. Once the
-   // parallel mesh is defined, the serial mesh can be deleted.
-   mfem::ParMesh pmesh(MPI_COMM_WORLD, mesh);
-   mesh.Clear();
-   {
-      for (int l = 0; l < par_ref_levels; l++)
-      {
-         pmesh.UniformRefinement();
-      }
-   }
-
-   // Define the Diffusion coefficient.
-   mfem::ConstantCoefficient* diffco=new mfem::ConstantCoefficient(1.0);
-   // Define the Heat source.
-   mfem::ConstantCoefficient* loadco=new mfem::ConstantCoefficient(1.0);
-   // Define the q-function.
-   mfem::QLinearDiffusion* qfun=new mfem::QLinearDiffusion(*diffco,*loadco,1.0,
-                                                           1e-7,4.0,0.5);
-
-   // Define FE collection and space for the state solution.
-   mfem::H1_FECollection sfec(order, dim);
-   mfem::ParFiniteElementSpace* sfes=new mfem::ParFiniteElementSpace(&pmesh,&sfec,
-                                                                     1);
-   // Define FE collection and space for the density field.
-   mfem::L2_FECollection pfec(order, dim);
-   mfem::ParFiniteElementSpace* pfes=new mfem::ParFiniteElementSpace(&pmesh,&pfec,
-                                                                     1);
-
-   // Define the arrays for the nonlinear form.
-   mfem::Array<mfem::ParFiniteElementSpace*> asfes;
-   mfem::Array<mfem::ParFiniteElementSpace*> apfes;
-
-   asfes.Append(sfes);
-   apfes.Append(pfes);
-
-   // Define parametric block nonlinear form using single scalar H1 field
-   // and L2 scalar density field.
-   mfem::ParParametricBNLForm* nf=new mfem::ParParametricBNLForm(asfes,apfes);
-   // Add a parametric integrator.
-   nf->AddDomainIntegrator(new mfem::ParametricLinearDiffusion(*qfun));
-
-   // Define true block vectors for state, adjoint, resudual.
-   mfem::BlockVector solbv; solbv.Update(nf->GetBlockTrueOffsets());    solbv=0.0;
-   mfem::BlockVector adjbv; adjbv.Update(nf->GetBlockTrueOffsets());    adjbv=0.0;
-   mfem::BlockVector resbv; resbv.Update(nf->GetBlockTrueOffsets());    resbv=0.0;
-   // Define true block vectors for parametric field and gradients.
-   mfem::BlockVector prmbv; prmbv.Update(nf->ParamGetBlockTrueOffsets());
-   prmbv=0.0;
-   mfem::BlockVector grdbv; grdbv.Update(nf->ParamGetBlockTrueOffsets());
-   grdbv=0.0;
-
-   // Set the BCs for the physics.
-   mfem::Array<mfem::Array<int> *> ess_bdr;
-   mfem::Array<mfem::Vector*>      ess_rhs;
-   ess_bdr.Append(new mfem::Array<int>(pmesh.bdr_attributes.Max()));
-   ess_rhs.Append(nullptr);
-   (*ess_bdr[0]) = 1;
-   nf->SetEssentialBC(ess_bdr,ess_rhs);
-   delete ess_bdr[0];
-
-   // Set the density field to 0.5.
-   prmbv=0.5;
-   // Set the density as parametric field in the parametric BNLForm.
-   nf->SetParamFields(prmbv); //set the density
-
-   // Compute the stiffness/tangent matrix for density prmbv=0.5.
-   mfem::BlockOperator *A = &nf->GetGradient(solbv);
-   mfem::HypreBoomerAMG* prec=new mfem::HypreBoomerAMG();
-   prec->SetPrintLevel(print_level);
-   // Use only block (0,0) as in this case we have a single field.
-   prec->SetOperator(A->GetBlock(0,0));
-
-   // Construct block preconditioner for the BNLForm.
-   mfem::BlockDiagonalPreconditioner *blpr = new mfem::BlockDiagonalPreconditioner(
-      nf->GetBlockTrueOffsets());
-   blpr->SetDiagonalBlock(0,prec);
-
-   // Define the solvers.
-   mfem::GMRESSolver *gmres;
-   gmres = new mfem::GMRESSolver(MPI_COMM_WORLD);
-   gmres->SetAbsTol(newton_abs_tol/10);
-   gmres->SetRelTol(newton_rel_tol/10);
-   gmres->SetMaxIter(100);
-   gmres->SetPrintLevel(print_level);
-   gmres->SetPreconditioner(*blpr);
-   gmres->SetOperator(*A);
-
-
-   // Solve the problem.
-   solbv=0.0;
-   nf->Mult(solbv,resbv); resbv.Neg(); //compute RHS
-   gmres->Mult(resbv, solbv);
-
-   // Compute the energy of the state system.
-   real_t energy = nf->GetEnergy(solbv);
-   if (myrank==0)
-   {
-      std::cout << "energy =" << energy << std::endl;
-   }
-
-   // Define the block nonlinear form utilized for representing the objective -
-   // use the state array from the BNLForm.
-   mfem::ParBlockNonlinearForm* ob=new mfem::ParBlockNonlinearForm(asfes);
-   // Add the integrator for the objective.
-   ob->AddDomainIntegrator(new mfem::DiffusionObjIntegrator());
-
-   // Compute the objective.
-   real_t obj=ob->GetEnergy(solbv);
-   if (myrank==0)
-   {
-      std::cout << "Objective =" << obj << std::endl;
-   }
-
-   // Solve the adjoint.
-   {
-      mfem::BlockVector adjrhs; adjrhs.Update(nf->GetBlockTrueOffsets());  adjrhs=0.0;
-      // Compute the RHS for the adjoint, i.e., the gradients with respect to
-      // the parametric fields.
-      ob->Mult(solbv, adjrhs);
-      // Get the tangent matrix from the state problem. We do not need to
-      // transpose the operator for diffusion.  Compute the adjoint solution.
-      gmres->Mult(adjrhs, adjbv);
-   }
-
-   // Compute gradients.
-   // First set the adjoint field.
-   nf->SetAdjointFields(adjbv);
-   // Set the state field.
-   nf->SetStateFields(solbv);
-   // Call the parametric Mult.
-   nf->ParamMult(prmbv, grdbv);
-
-   // Dump out the data.
-   if (visualization)
-   {
-      mfem::ParaViewDataCollection *dacol=new mfem::ParaViewDataCollection("ParHeat",
-                                                                           &pmesh);
-      mfem::ParGridFunction gfgrd(pfes); gfgrd.SetFromTrueDofs(grdbv.GetBlock(0));
-      mfem::ParGridFunction gfdns(pfes); gfdns.SetFromTrueDofs(prmbv.GetBlock(0));
-      // Define state grid function.
-      mfem::ParGridFunction gfsol(sfes); gfsol.SetFromTrueDofs(solbv.GetBlock(0));
-      mfem::ParGridFunction gfadj(sfes); gfadj.SetFromTrueDofs(adjbv.GetBlock(0));
-
-      dacol->SetLevelsOfDetail(order);
-      dacol->RegisterField("sol", &gfsol);
-      dacol->RegisterField("adj", &gfadj);
-      dacol->RegisterField("dns", &gfdns);
-      dacol->RegisterField("grd", &gfgrd);
-
-      dacol->SetTime(1.0);
-      dacol->SetCycle(1);
-      dacol->Save();
-
-      delete dacol;
-   }
-
-   // FD check
-   {
-      mfem::BlockVector prtbv;
-      mfem::BlockVector tmpbv;
-      prtbv.Update(nf->ParamGetBlockTrueOffsets());
-      tmpbv.Update(nf->ParamGetBlockTrueOffsets());
-      prtbv.GetBlock(0).Randomize();
-      prtbv*=1.0;
-      real_t lsc=1.0;
-
-      real_t gQoI=ob->GetEnergy(solbv);
-      real_t lQoI;
-
-      real_t nd=mfem::InnerProduct(MPI_COMM_WORLD,prtbv,prtbv);
-      real_t td=mfem::InnerProduct(MPI_COMM_WORLD,prtbv,grdbv);
-      td=td/nd;
-
-      for (int l = 0; l < 10; l++)
-      {
-         lsc/=10.0;
-         prtbv/=10.0;
-         add(prmbv,prtbv,tmpbv);
-         nf->SetParamFields(tmpbv);
-         // Solve the physics.
-         solbv=0.0;
-         nf->Mult(solbv,resbv); resbv.Neg(); //compute RHS
-         A = &nf->GetGradient(solbv);
-         prec->SetPrintLevel(0);
-         prec->SetOperator(A->GetBlock(0,0));
-         gmres->SetOperator(*A);
-         gmres->SetPrintLevel(0);
-         gmres->Mult(resbv,solbv);
-         // Compute the objective.
-         lQoI=ob->GetEnergy(solbv);
-         real_t ld=(lQoI-gQoI)/lsc;
-         if (myrank==0)
-         {
-            std::cout << "dx=" << lsc <<" FD approximation=" << ld/nd
-                      << " adjoint gradient=" << td
-                      << " err=" << std::fabs(ld/nd-td) << std::endl;
-         }
-      }
-   }
-
-   delete ob;
-   delete gmres;
-   delete blpr;
-   delete prec;
-
-   delete nf;
-   delete pfes;
-   delete sfes;
-
-   delete qfun;
-   delete loadco;
-   delete diffco;
-
-   return 0;
-}

miniapps/mtop/pparamnonlinearform.cpp

@@ -1,362 +0,0 @@
-// 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 "pparamnonlinearform.hpp"
-
-#ifdef MFEM_USE_MPI
-
-namespace mfem
-{
-
-ParParametricBNLForm::ParParametricBNLForm(Array<ParFiniteElementSpace *>
-                                           &statef,
-                                           Array<ParFiniteElementSpace *> &paramf)
-   :ParametricBNLForm()
-{
-   pBlockGrad = nullptr;
-   SetParSpaces(statef,paramf);
-}
-
-void ParParametricBNLForm::SetParSpaces(Array<ParFiniteElementSpace *> &statef,
-                                        Array<ParFiniteElementSpace *> &paramf)
-{
-   delete pBlockGrad;
-   pBlockGrad = nullptr;
-
-   for (int s1=0; s1<fes.Size(); ++s1)
-   {
-      for (int s2=0; s2<fes.Size(); ++s2)
-      {
-         delete phBlockGrad(s1,s2);
-      }
-   }
-
-   Array<FiniteElementSpace *> serialSpaces(statef.Size());
-   Array<FiniteElementSpace *> prmserialSpaces(paramf.Size());
-   for (int s=0; s<statef.Size(); s++)
-   {
-      serialSpaces[s] = (FiniteElementSpace *) statef[s];
-   }
-   for (int s=0; s<paramf.Size(); s++)
-   {
-      prmserialSpaces[s] = (FiniteElementSpace *) paramf[s];
-   }
-
-   SetSpaces(serialSpaces,prmserialSpaces);
-
-   phBlockGrad.SetSize(fes.Size(), fes.Size());
-
-   for (int s1=0; s1<fes.Size(); ++s1)
-   {
-      for (int s2=0; s2<fes.Size(); ++s2)
-      {
-         phBlockGrad(s1,s2) = new OperatorHandle(Operator::Hypre_ParCSR);
-      }
-   }
-}
-
-ParFiniteElementSpace * ParParametricBNLForm::ParFESpace(int k)
-{
-   return (ParFiniteElementSpace *)fes[k];
-}
-
-const ParFiniteElementSpace *ParParametricBNLForm::ParFESpace(int k) const
-{
-   return (const ParFiniteElementSpace *)fes[k];
-}
-
-
-ParFiniteElementSpace * ParParametricBNLForm::ParParamFESpace(int k)
-{
-   return (ParFiniteElementSpace *)paramfes[k];
-}
-
-const ParFiniteElementSpace *ParParametricBNLForm::ParParamFESpace(int k) const
-{
-   return (const ParFiniteElementSpace *)paramfes[k];
-}
-
-// Here, rhs is a true dof vector
-void ParParametricBNLForm::SetEssentialBC(const
-                                          Array<Array<int> *>&bdr_attr_is_ess,
-                                          Array<Vector *> &rhs)
-{
-   Array<Vector *> nullarray(fes.Size());
-   nullarray = NULL;
-
-   ParametricBNLForm::SetEssentialBC(bdr_attr_is_ess, nullarray);
-
-   for (int s = 0; s < fes.Size(); ++s)
-   {
-      if (rhs[s])
-      {
-         rhs[s]->SetSubVector(*ess_tdofs[s], 0.0);
-      }
-   }
-}
-
-void ParParametricBNLForm::SetParamEssentialBC(const
-                                               Array<Array<int> *>&bdr_attr_is_ess,
-                                               Array<Vector *> &rhs)
-{
-   Array<Vector *> nullarray(fes.Size());
-   nullarray = NULL;
-
-   ParametricBNLForm::SetParamEssentialBC(bdr_attr_is_ess, nullarray);
-
-   for (int s = 0; s < paramfes.Size(); ++s)
-   {
-      if (rhs[s])
-      {
-         rhs[s]->SetSubVector(*paramess_tdofs[s], 0.0);
-      }
-   }
-}
-
-real_t ParParametricBNLForm::GetEnergy(const Vector &x) const
-{
-   xs_true.Update(const_cast<Vector&>(x), block_trueOffsets);
-   xs.Update(block_offsets);
-
-   for (int s = 0; s < fes.Size(); ++s)
-   {
-      fes[s]->GetProlongationMatrix()->Mult(xs_true.GetBlock(s), xs.GetBlock(s));
-   }
-
-   real_t enloc = ParametricBNLForm::GetEnergyBlocked(xs,xdv);
-   real_t englo = 0.0;
-
-   MPI_Allreduce(&enloc, &englo, 1, MPITypeMap<real_t>::mpi_type, MPI_SUM,
-                 ParFESpace(0)->GetComm());
-
-   return englo;
-}
-
-void ParParametricBNLForm::Mult(const Vector &x, Vector &y) const
-{
-   xs_true.Update(const_cast<Vector&>(x), block_trueOffsets);
-   ys_true.Update(y, block_trueOffsets);
-   xs.Update(block_offsets);
-   ys.Update(block_offsets);
-
-   for (int s=0; s<fes.Size(); ++s)
-   {
-      fes[s]->GetProlongationMatrix()->Mult(
-         xs_true.GetBlock(s), xs.GetBlock(s));
-   }
-
-   ParametricBNLForm::MultBlocked(xs, xdv, ys);
-
-   if (fnfi.Size() > 0)
-   {
-      MFEM_ABORT("TODO: assemble contributions from shared face terms");
-   }
-
-   for (int s=0; s<fes.Size(); ++s)
-   {
-      fes[s]->GetProlongationMatrix()->MultTranspose(
-         ys.GetBlock(s), ys_true.GetBlock(s));
-
-      ys_true.GetBlock(s).SetSubVector(*ess_tdofs[s], 0.0);
-   }
-}
-
-/// Block T-Vector to Block T-Vector
-void ParParametricBNLForm::ParamMult(const Vector &x, Vector &y) const
-{
-   xs_true.Update(const_cast<Vector&>(x), paramblock_trueOffsets);
-   ys_true.Update(y, paramblock_trueOffsets);
-   prmxs.Update(paramblock_offsets);
-   prmys.Update(paramblock_offsets);
-
-   for (int s=0; s<paramfes.Size(); ++s)
-   {
-      paramfes[s]->GetProlongationMatrix()->Mult(
-         xs_true.GetBlock(s), prmxs.GetBlock(s));
-   }
-
-   ParametricBNLForm::MultParamBlocked(xsv,adv,xdv,prmys);
-
-   if (fnfi.Size() > 0)
-   {
-      MFEM_ABORT("TODO: assemble contributions from shared face terms");
-   }
-
-   for (int s=0; s<paramfes.Size(); ++s)
-   {
-      paramfes[s]->GetProlongationMatrix()->MultTranspose(
-         prmys.GetBlock(s), ys_true.GetBlock(s));
-
-      ys_true.GetBlock(s).SetSubVector(*paramess_tdofs[s], 0.0);
-   }
-
-}
-
-/// Return the local gradient matrix for the given true-dof vector x
-const BlockOperator & ParParametricBNLForm::GetLocalGradient(
-   const Vector &x) const
-{
-   xs_true.Update(const_cast<Vector&>(x), block_trueOffsets);
-   xs.Update(block_offsets);
-
-   for (int s=0; s<fes.Size(); ++s)
-   {
-      fes[s]->GetProlongationMatrix()->Mult(
-         xs_true.GetBlock(s), xs.GetBlock(s));
-   }
-
-   ParametricBNLForm::ComputeGradientBlocked(xs,
-                                             xdv); // (re)assemble Grad with b.c.
-
-   delete BlockGrad;
-   BlockGrad = new BlockOperator(block_offsets);
-
-   for (int i = 0; i < fes.Size(); ++i)
-   {
-      for (int j = 0; j < fes.Size(); ++j)
-      {
-         BlockGrad->SetBlock(i, j, Grads(i, j));
-      }
-   }
-   return *BlockGrad;
-}
-
-// Set the operator type id for the parallel gradient matrix/operator.
-void ParParametricBNLForm::SetGradientType(Operator::Type tid)
-{
-   for (int s1=0; s1<fes.Size(); ++s1)
-   {
-      for (int s2=0; s2<fes.Size(); ++s2)
-      {
-         phBlockGrad(s1,s2)->SetType(tid);
-      }
-   }
-}
-
-BlockOperator & ParParametricBNLForm::GetGradient(const Vector &x) const
-{
-   if (pBlockGrad == NULL)
-   {
-      pBlockGrad = new BlockOperator(block_trueOffsets);
-   }
-
-   Array<const ParFiniteElementSpace *> pfes(fes.Size());
-
-   for (int s1=0; s1<fes.Size(); ++s1)
-   {
-      pfes[s1] = ParFESpace(s1);
-
-      for (int s2=0; s2<fes.Size(); ++s2)
-      {
-         phBlockGrad(s1,s2)->Clear();
-      }
-   }
-
-   GetLocalGradient(x); // gradients are stored in 'Grads'
-
-   if (fnfi.Size() > 0)
-   {
-      MFEM_ABORT("TODO: assemble contributions from shared face terms");
-   }
-
-   for (int s1=0; s1<fes.Size(); ++s1)
-   {
-      for (int s2=0; s2<fes.Size(); ++s2)
-      {
-         OperatorHandle dA(phBlockGrad(s1,s2)->Type()),
-                        Ph(phBlockGrad(s1,s2)->Type()),
-                        Rh(phBlockGrad(s1,s2)->Type());
-
-         if (s1 == s2)
-         {
-            dA.MakeSquareBlockDiag(pfes[s1]->GetComm(), pfes[s1]->GlobalVSize(),
-                                   pfes[s1]->GetDofOffsets(), Grads(s1,s1));
-            Ph.ConvertFrom(pfes[s1]->Dof_TrueDof_Matrix());
-            phBlockGrad(s1,s1)->MakePtAP(dA, Ph);
-
-            OperatorHandle Ae;
-            Ae.EliminateRowsCols(*phBlockGrad(s1,s1), *ess_tdofs[s1]);
-         }
-         else
-         {
-            dA.MakeRectangularBlockDiag(pfes[s1]->GetComm(),
-                                        pfes[s1]->GlobalVSize(),
-                                        pfes[s2]->GlobalVSize(),
-                                        pfes[s1]->GetDofOffsets(),
-                                        pfes[s2]->GetDofOffsets(),
-                                        Grads(s1,s2));
-            Rh.ConvertFrom(pfes[s1]->Dof_TrueDof_Matrix());
-            Ph.ConvertFrom(pfes[s2]->Dof_TrueDof_Matrix());
-
-            phBlockGrad(s1,s2)->MakeRAP(Rh, dA, Ph);
-
-            phBlockGrad(s1,s2)->EliminateRows(*ess_tdofs[s1]);
-            phBlockGrad(s1,s2)->EliminateCols(*ess_tdofs[s2]);
-         }
-
-         pBlockGrad->SetBlock(s1, s2, phBlockGrad(s1,s2)->Ptr());
-      }
-   }
-
-   return *pBlockGrad;
-}
-
-ParParametricBNLForm::~ParParametricBNLForm()
-{
-   delete pBlockGrad;
-   for (int s1=0; s1<fes.Size(); ++s1)
-   {
-      for (int s2=0; s2<fes.Size(); ++s2)
-      {
-         delete phBlockGrad(s1,s2);
-      }
-   }
-}
-
-
-void ParParametricBNLForm::SetStateFields(const Vector &xv) const
-{
-   xs_true.Update(const_cast<Vector&>(xv), block_trueOffsets);
-   xsv.Update(block_offsets);
-   for (int s=0; s<fes.Size(); ++s)
-   {
-      fes[s]->GetProlongationMatrix()->Mult(
-         xs_true.GetBlock(s), xsv.GetBlock(s));
-   }
-}
-
-
-void ParParametricBNLForm::SetAdjointFields(const Vector &av) const
-{
-   xs_true.Update(const_cast<Vector&>(av), block_trueOffsets);
-   adv.Update(block_offsets);
-   for (int s=0; s<fes.Size(); ++s)
-   {
-      fes[s]->GetProlongationMatrix()->Mult(
-         xs_true.GetBlock(s), adv.GetBlock(s));
-   }
-}
-
-void ParParametricBNLForm::SetParamFields(const Vector &dv) const
-{
-   xs_true.Update(const_cast<Vector&>(dv),paramblock_trueOffsets);
-   xdv.Update(paramblock_offsets);
-   for (int s=0; s<paramfes.Size(); ++s)
-   {
-      paramfes[s]->GetProlongationMatrix()->Mult(
-         xs_true.GetBlock(s), xdv.GetBlock(s));
-   }
-}
-
-}
-
-#endif

miniapps/mtop/pparamnonlinearform.hpp

@@ -1,114 +0,0 @@
-// 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_PPRMNONLINEARFORM
-#define MFEM_PPRMNONLINEARFORM
-
-
-#ifdef MFEM_USE_MPI
-
-#include "mfem.hpp"
-#include "paramnonlinearform.hpp"
-
-namespace mfem
-{
-
-/** @brief A class representing a general parametric parallel block nonlinear
-    operator defined on the Cartesian product of multiple
-    ParFiniteElementSpace%s. */
-/** The ParParametricBNLForm takes as input, and returns as output, vectors on
-    the true dofs. */
-class ParParametricBNLForm : public ParametricBNLForm
-{
-protected:
-   mutable BlockVector xs_true, ys_true;
-   mutable Array2D<OperatorHandle *> phBlockGrad;
-   mutable BlockOperator *pBlockGrad;
-
-public:
-   /// Computes the energy of the system
-   real_t GetEnergy(const Vector &x) const override;
-
-   /// Construct an empty ParParametricBNLForm. Initialize with SetParSpaces().
-   ParParametricBNLForm() : pBlockGrad(nullptr) { }
-
-   /** @brief Construct a ParParametricBNLForm on the given set of
-       parametric and state ParFiniteElementSpace%s. */
-   ParParametricBNLForm(Array<ParFiniteElementSpace *> &statef,
-                        Array<ParFiniteElementSpace *> &paramf);
-
-   /// Return the @a k-th parallel FE state space of the ParParametricBNLForm.
-   ParFiniteElementSpace *ParFESpace(int k);
-   /** @brief Return the @a k-th parallel FE state space of the
-       ParParametricBNLForm (const version). */
-   const ParFiniteElementSpace *ParFESpace(int k) const;
-
-   /// Return the @a k-th parallel FE parameters space of the
-   /// ParParametricBNLForm.
-   ParFiniteElementSpace *ParParamFESpace(int k);
-   /** @brief Return the @a k-th parallel FE parameters space of the
-       ParParametricBNLForm (const version). */
-   const ParFiniteElementSpace *ParParamFESpace(int k) const;
-
-   /** @brief Set the parallel FE spaces for the state and the parametric
-    *  fields. After a call to SetParSpaces(), the essential b.c. and the
-    *  gradient-type (if different from the default) must be set again. */
-   void SetParSpaces(Array<ParFiniteElementSpace *> &statef,
-                     Array<ParFiniteElementSpace *> &paramf);
-
-   /// Set the state essential BCs. Here, rhs is a true dof vector!
-   void SetEssentialBC(const Array<Array<int> *>&bdr_attr_is_ess,
-                       Array<Vector *> &rhs) override;
-
-   // Set the essential BCs for the parametric fields. Here, rhs is a true dof
-   // vector!
-   void SetParamEssentialBC(const Array<Array<int> *>&bdr_attr_is_ess,
-                            Array<Vector *> &rhs) override;
-
-
-   /** @brief Calculates the residual for a state input given by block T-Vector.
-    *  The result is Block T-Vector! The parametric fields should be set in
-    *  advance by calling SetParamFields(). */
-   void Mult(const Vector &x, Vector &y) const override;
-
-   /** @brief Calculates the product of the adjoint field and the derivative of
-    *  the state residual with respect to the parametric fields. The adjoint and
-    *  the state fields should be set in advance by calling SetAdjointFields()
-    *  and SetStateFields(). The input and the result are block T-Vectors!*/
-   void ParamMult(const Vector &x, Vector &y) const override;
-
-   /// Return the local block gradient matrix for the given true-dof vector x
-   const BlockOperator &GetLocalGradient(const Vector &x) const;
-
-   /// Return the block gradient matrix for the given true-dof vector x
-   BlockOperator &GetGradient(const Vector &x) const override;
-
-   /** @brief Set the operator type id for the blocks of the parallel gradient
-       matrix/operator. The default type is Operator::Hypre_ParCSR. */
-   void SetGradientType(Operator::Type tid);
-
-   /// Destructor.
-   virtual ~ParParametricBNLForm();
-
-   /// Set the state fields
-   void SetStateFields(const Vector &xv) const override;
-
-   /// Set the adjoint fields
-   void SetAdjointFields(const Vector &av) const override;
-
-   /// Set the parameters/design fields
-   void SetParamFields(const Vector &dv) const override;
-};
-
-}
-
-#endif
-#endif

miniapps/mtop/seqheat.cpp

@@ -1,308 +0,0 @@
-// 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.
-//
-//     ----------------------------------------------------------------
-//     SeqHeat Miniapp: Gradients of PDE constrained objective function
-//     ----------------------------------------------------------------
-//                          (Sequential Version)
-//
-// The following example computes the gradients of a specified objective
-// function with respect to parametric fields. The objective function is having
-// the following form f(u(\rho)) where u(\rho) is a solution of a specific state
-// problem (in the example that is the diffusion equation), and \rho is a
-// parametric field discretized by finite elements. The parametric field (also
-// called density in topology optimization) controls the coefficients of the
-// state equation. For the considered case, the density controls the diffusion
-// coefficient within the computational domain.
-//
-// For more information, the users are referred to:
-//
-//    Hinze, M.; Pinnau, R.; Ulbrich, M. & Ulbrich, S.
-//    Optimization with PDE Constraints
-//    Springer Netherlands, 2009
-//
-//    Bendsøe, M. P. & Sigmund, O.
-//    Topology Optimization - Theory, Methods and Applications
-//    Springer Verlag, Berlin Heidelberg, 2003
-//
-// Compile with: make seqheat
-//
-// Sample runs:
-//
-//    seqheat -m ../../data/star-mixed.mesh
-//    seqheat --visualization
-
-#include "mfem.hpp"
-#include <fstream>
-#include <iostream>
-
-#include "mtop_integrators.hpp"
-
-using namespace mfem;
-
-int main(int argc, char *argv[])
-{
-   const char *mesh_file = "../../data/star.vtk";
-   int ser_ref_levels = 1;
-   int order = 2;
-   bool visualization = false;
-   real_t newton_rel_tol = 1e-4;
-   real_t newton_abs_tol = 1e-6;
-   int newton_iter = 10;
-   int print_level = 0;
-
-   mfem::OptionsParser args(argc, argv);
-   args.AddOption(&mesh_file, "-m", "--mesh", "Mesh file to use.");
-   args.AddOption(&ser_ref_levels,
-                  "-rs",
-                  "--refine-serial",
-                  "Number of times to refine the mesh uniformly in serial.");
-   args.AddOption(&order,
-                  "-o",
-                  "--order",
-                  "Order (degree) of the finite elements.");
-   args.AddOption(&visualization,
-                  "-vis",
-                  "--visualization",
-                  "-no-vis",
-                  "--no-visualization",
-                  "Enable or disable GLVis visualization.");
-   args.AddOption(&newton_rel_tol,
-                  "-rel",
-                  "--relative-tolerance",
-                  "Relative tolerance for the Newton solve.");
-   args.AddOption(&newton_abs_tol,
-                  "-abs",
-                  "--absolute-tolerance",
-                  "Absolute tolerance for the Newton solve.");
-   args.AddOption(&newton_iter,
-                  "-it",
-                  "--newton-iterations",
-                  "Maximum iterations for the Newton solve.");
-   args.Parse();
-   if (!args.Good())
-   {
-      args.PrintUsage(std::cout);
-      return 1;
-   }
-   args.PrintOptions(std::cout);
-
-   // Read the (serial) mesh from the given mesh file on all processors. We
-   // can handle triangular, quadrilateral, tetrahedral and hexahedral meshes
-   // with the same code.
-   mfem::Mesh *mesh = new mfem::Mesh(mesh_file, 1, 1);
-   int dim = mesh->Dimension();
-
-   // Refine the mesh in serial to increase the resolution. In this example
-   // we do 'ser_ref_levels' of uniform refinement, where 'ser_ref_levels' is
-   // a command-line parameter.
-   for (int lev = 0; lev < ser_ref_levels; lev++)
-   {
-      mesh->UniformRefinement();
-   }
-
-   // Diffusion coefficient
-   mfem::ConstantCoefficient* diffco=new mfem::ConstantCoefficient(1.0);
-   // Heat source
-   mfem::ConstantCoefficient* loadco=new mfem::ConstantCoefficient(1.0);
-   // Define the q-function
-   mfem::QLinearDiffusion* qfun=new mfem::QLinearDiffusion(*diffco,*loadco,1.0,
-                                                           1e-7,4.0,0.5);
-
-   // Define FE collection and space for the state solution
-   mfem::H1_FECollection sfec(order, dim);
-   mfem::FiniteElementSpace* sfes=new mfem::FiniteElementSpace(mesh,&sfec,1);
-   // Define FE collection and space for the density field
-   mfem::L2_FECollection pfec(order, dim);
-   mfem::FiniteElementSpace* pfes=new mfem::FiniteElementSpace(mesh,&pfec,1);
-
-   // Define the arrays for the nonlinear form
-   mfem::Array<mfem::FiniteElementSpace*> asfes;
-   mfem::Array<mfem::FiniteElementSpace*> apfes;
-
-   asfes.Append(sfes);
-   apfes.Append(pfes);
-   // Define parametric block nonlinear form using single scalar H1 field
-   // and L2 scalar density field
-   mfem::ParametricBNLForm* nf=new mfem::ParametricBNLForm(asfes,apfes);
-   // Add the parametric integrator
-   nf->AddDomainIntegrator(new mfem::ParametricLinearDiffusion(*qfun));
-
-   // Define true block vectors for state, adjoint, residual
-   mfem::BlockVector solbv; solbv.Update(nf->GetBlockTrueOffsets());    solbv=0.0;
-   mfem::BlockVector adjbv; adjbv.Update(nf->GetBlockTrueOffsets());    adjbv=0.0;
-   mfem::BlockVector resbv; resbv.Update(nf->GetBlockTrueOffsets());    resbv=0.0;
-   // Define true block vectors for parametric field and gradients
-   mfem::BlockVector prmbv; prmbv.Update(nf->ParamGetBlockTrueOffsets());
-   prmbv=0.0;
-   mfem::BlockVector grdbv; grdbv.Update(nf->ParamGetBlockTrueOffsets());
-   grdbv=0.0;
-
-   // Set the BC for the physics
-   mfem::Array<mfem::Array<int> *> ess_bdr;
-   mfem::Array<mfem::Vector*>      ess_rhs;
-   ess_bdr.Append(new mfem::Array<int>(mesh->bdr_attributes.Max()));
-   ess_rhs.Append(nullptr);
-   (*ess_bdr[0]) = 1;
-   nf->SetEssentialBC(ess_bdr,ess_rhs);
-   delete ess_bdr[0];
-
-   // Define the linear solvers
-   mfem::GMRESSolver *gmres;
-   gmres = new mfem::GMRESSolver();
-   gmres->SetAbsTol(newton_abs_tol/10);
-   gmres->SetRelTol(newton_rel_tol/10);
-   gmres->SetMaxIter(300);
-   gmres->SetPrintLevel(print_level);
-
-   // Define the Newton solver
-   mfem::NewtonSolver *ns;
-   ns = new mfem::NewtonSolver();
-   ns->iterative_mode = true;
-   ns->SetSolver(*gmres);
-   ns->SetOperator(*nf);
-   ns->SetPrintLevel(print_level);
-   ns->SetRelTol(newton_rel_tol);
-   ns->SetAbsTol(newton_abs_tol);
-   ns->SetMaxIter(newton_iter);
-
-   // Solve the problem
-   // Set the density to 0.5
-   prmbv=0.5;
-   nf->SetParamFields(prmbv); // Set the density
-   // Define the RHS
-   mfem::Vector b;
-   solbv=0.0;
-   // Newton solve
-   ns->Mult(b, solbv);
-
-   // Compute the residual
-   nf->Mult(solbv,resbv);
-   std::cout<<"Norm residual="<<resbv.Norml2()<<std::endl;
-
-   // Compute the energy of the state system
-   real_t energy = nf->GetEnergy(solbv);
-   std::cout<<"energy ="<< energy<<std::endl;
-
-   // Define the block nonlinear form utilized for representing the
-   // objective. The input is the state array asfes defined earlier.
-   mfem::BlockNonlinearForm* ob=new mfem::BlockNonlinearForm(asfes);
-
-   // Add the integrator for the objective
-   ob->AddDomainIntegrator(new mfem::DiffusionObjIntegrator());
-
-   // Compute the objective
-   real_t obj=ob->GetEnergy(solbv);
-   std::cout<<"Objective ="<<obj<<std::endl;
-
-   // Solve the adjoint
-   {
-      mfem::BlockVector adjrhs; adjrhs.Update(nf->GetBlockTrueOffsets());  adjrhs=0.0;
-      // Compute the RHS for the adjoint
-      ob->Mult(solbv, adjrhs);
-      // Get the tangent matrix from the state problem
-      mfem::BlockOperator& A=nf->GetGradient(solbv);
-      // We do not need to transpose the operator for diffusion
-      gmres->SetOperator(A.GetBlock(0,0));
-      // Compute the adjoint solution
-      gmres->Mult(adjrhs.GetBlock(0), adjbv.GetBlock(0));
-   }
-
-   // Compute gradients
-   nf->SetAdjointFields(adjbv);
-   nf->SetStateFields(solbv);
-   nf->ParamMult(prmbv, grdbv);
-
-   // Dump out the data
-   if (visualization)
-   {
-      mfem::ParaViewDataCollection *dacol=new mfem::ParaViewDataCollection("SeqHeat",
-                                                                           mesh);
-      mfem::GridFunction gfgrd(pfes); gfgrd.SetFromTrueDofs(grdbv.GetBlock(0));
-      mfem::GridFunction gfdns(pfes); gfdns.SetFromTrueDofs(prmbv.GetBlock(0));
-      // Define state grid function
-      mfem::GridFunction gfsol(sfes); gfsol.SetFromTrueDofs(solbv.GetBlock(0));
-      mfem::GridFunction gfadj(sfes); gfadj.SetFromTrueDofs(adjbv.GetBlock(0));
-
-      dacol->SetLevelsOfDetail(order);
-      dacol->RegisterField("sol", &gfsol);
-      dacol->RegisterField("adj", &gfadj);
-      dacol->RegisterField("dns", &gfdns);
-      dacol->RegisterField("grd", &gfgrd);
-
-      dacol->SetTime(1.0);
-      dacol->SetCycle(1);
-      dacol->Save();
-
-      delete dacol;
-   }
-
-   // FD check
-   {
-      // Perturbation vector
-      mfem::BlockVector prtbv;
-      mfem::BlockVector tmpbv;
-      prtbv.Update(nf->ParamGetBlockTrueOffsets());
-      tmpbv.Update(nf->ParamGetBlockTrueOffsets());
-      // Generate the perturbation
-      prtbv.GetBlock(0).Randomize();
-      prtbv*=1.0;
-      // Scaling parameter
-      real_t lsc=1.0;
-
-      // Compute initial objective
-      real_t gQoI=ob->GetEnergy(solbv);
-      real_t lQoI;
-
-      // Norm of the perturbation
-      real_t nd=mfem::InnerProduct(prtbv,prtbv);
-      // Projection of the adjoint gradient on the perturbation
-      real_t td=mfem::InnerProduct(prtbv,grdbv);
-      // Normalize the directional derivative
-      td=td/nd;
-
-      for (int l = 0; l < 10; l++)
-      {
-         lsc/=10.0;
-         // Scale the perturbation
-         prtbv/=10.0;
-         // Add the perturbation to the original density
-         add(prmbv,prtbv,tmpbv);
-         nf->SetParamFields(tmpbv);
-         // Solve the physics
-         ns->Mult(b,solbv);
-         // Compute the objective
-         lQoI=ob->GetEnergy(solbv);
-         // FD approximation
-         real_t ld=(lQoI-gQoI)/lsc;
-         std::cout << "dx=" << lsc << " FD gradient=" << ld/nd
-                   << " adjoint gradient=" << td
-                   << " err=" << std::fabs(ld/nd-td) << std::endl;
-      }
-   }
-
-   delete ob;
-
-   delete ns;
-   delete gmres;
-
-   delete nf;
-   delete pfes;
-   delete sfes;
-
-   delete qfun;
-   delete loadco;
-   delete diffco;
-
-   delete mesh;
-
-   return 0;
-}

miniapps/mtop/square-extended.mesh

@@ -0,0 +1,35 @@
+MFEM mesh v1.0
+
+#
+# MFEM Geometry Types (see fem/geom.hpp):
+#
+# POINT       = 0
+# SEGMENT     = 1
+# TRIANGLE    = 2
+# SQUARE      = 3
+# TETRAHEDRON = 4
+# CUBE        = 5
+# PRISM       = 6
+#
+
+dimension
+2
+
+elements
+1
+1 3 0 1 2 3
+
+boundary
+4
+1 1 0 1
+2 1 1 2
+3 1 2 3
+4 1 3 0
+
+vertices
+4
+2
+-1 -1
+1 -1
+1 1
+-1 1

miniapps/solvers/lor_solvers.cpp

@@ -109,7 +109,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();
@@ -156,7 +156,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);
@@ -167,6 +167,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));
    }
+   if (H1) { b.UseFastAssembly(true); }
    b.Assemble();
 
    GridFunction x(&fes);

miniapps/solvers/plor_solvers.cpp

@@ -112,7 +112,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();
@@ -166,7 +166,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);
@@ -177,6 +177,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));
    }
+   if (H1) { b.UseFastAssembly(true); }
    b.Assemble();
 
    ParGridFunction x(&fes);

tests/unit/CMakeLists.txt

@@ -63,6 +63,7 @@ set(UNIT_TESTS_SRCS
   mesh/test_vtu.cpp
   mesh/test_nurbs.cpp
   mesh/test_exodus_writer.cpp
+  fem/make_permuted_mesh.cpp
   fem/test_1d_bilininteg.cpp
   fem/test_2d_bilininteg.cpp
   fem/test_3d_bilininteg.cpp
@@ -107,6 +108,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][GPU]")
 
 #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,382 @@
+// 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;
+
+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);
+
+      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 q = 0; q < ir->GetNPoints(); q++)
+      {
+         const IntegrationPoint &ip = ir->IntPoint(q);
+         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;
+
+         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 real_t 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