Some additions and tweaks to get ex14p to work with empty partitions

Fix integration rule for the parallel test

examples/ex14.cpp

@@ -19,6 +19,10 @@
 //               ex14 -m ../data/amr-hex.mesh
 //               ex14 -m ../data/fichera-amr.mesh
 //
+// Device sample runs:
+//               ex14 -pa -d cuda -o 3
+//               ex14 -pa -d cuda -o 3 -m ../data/fichera.mesh
+//
 // Description:  This example code demonstrates the use of MFEM to define a
 //               discontinuous Galerkin (DG) finite element discretization of
 //               the Laplace problem -Delta u = 1 with homogeneous Dirichlet
@@ -46,7 +50,9 @@ int main(int argc, char *argv[])
    real_t sigma = -1.0;
    real_t kappa = -1.0;
    real_t eta = 0.0;
+   bool pa = false;
    bool visualization = 1;
+   const char *device_config = "cpu";
 
    OptionsParser args(argc, argv);
    args.AddOption(&mesh_file, "-m", "--mesh",
@@ -62,9 +68,13 @@ int main(int argc, char *argv[])
                   "One of the three DG penalty parameters, should be positive."
                   " Negative values are replaced with (order+1)^2.");
    args.AddOption(&eta, "-e", "--eta", "BR2 penalty parameter.");
+   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
+                  "--no-partial-assembly", "Enable Partial Assembly.");
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
    args.Parse();
    if (!args.Good())
    {
@@ -77,117 +87,134 @@ int main(int argc, char *argv[])
    }
    args.PrintOptions(cout);
 
-   // 2. Read the mesh from the given mesh file. We can handle triangular,
+   // 2. Enable hardware devices such as GPUs, and programming models such as
+   //    CUDA, OCCA, RAJA and OpenMP based on command line options.
+   Device device(device_config);
+   device.Print();
+
+   // 3. Read the mesh from the given mesh file. We can handle triangular,
    //    quadrilateral, tetrahedral and hexahedral meshes with the same code.
    //    NURBS meshes are projected to second order meshes.
-   Mesh *mesh = new Mesh(mesh_file, 1, 1);
-   int dim = mesh->Dimension();
+   Mesh mesh(mesh_file);
+   const int dim = mesh.Dimension();
 
-   // 3. Refine the mesh to increase the resolution. In this example we do
+   // 4. Refine the mesh to increase the resolution. In this example we do
    //    'ref_levels' of uniform refinement. By default, or if ref_levels < 0,
    //    we choose it to be the largest number that gives a final mesh with no
    //    more than 50,000 elements.
    {
       if (ref_levels < 0)
       {
-         ref_levels = (int)floor(log(50000./mesh->GetNE())/log(2.)/dim);
+         ref_levels = (int)floor(log(50000./mesh.GetNE())/log(2.)/dim);
       }
       for (int l = 0; l < ref_levels; l++)
       {
-         mesh->UniformRefinement();
+         mesh.UniformRefinement();
       }
    }
-   if (mesh->NURBSext)
+   if (mesh.NURBSext)
    {
-      mesh->SetCurvature(max(order, 1));
+      mesh.SetCurvature(max(order, 1));
    }
 
-   // 4. Define a finite element space on the mesh. Here we use discontinuous
+   // 5. Define a finite element space on the mesh. Here we use discontinuous
    //    finite elements of the specified order >= 0.
-   FiniteElementCollection *fec = new DG_FECollection(order, dim);
-   FiniteElementSpace *fespace = new FiniteElementSpace(mesh, fec);
-   cout << "Number of unknowns: " << fespace->GetVSize() << endl;
+   DG_FECollection fec(order, dim, BasisType::GaussLobatto);
+   FiniteElementSpace fespace(&mesh, &fec);
+   cout << "Number of unknowns: " << fespace.GetVSize() << endl;
 
-   // 5. Set up the linear form b(.) which corresponds to the right-hand side of
+   // 6. Set up the linear form b(.) which corresponds to the right-hand side of
    //    the FEM linear system.
-   LinearForm *b = new LinearForm(fespace);
+   LinearForm b(&fespace);
    ConstantCoefficient one(1.0);
    ConstantCoefficient zero(0.0);
-   b->AddDomainIntegrator(new DomainLFIntegrator(one));
-   b->AddBdrFaceIntegrator(
+   b.AddDomainIntegrator(new DomainLFIntegrator(one));
+   b.AddBdrFaceIntegrator(
       new DGDirichletLFIntegrator(zero, one, sigma, kappa));
-   b->Assemble();
+   b.Assemble();
 
-   // 6. Define the solution vector x as a finite element grid function
+   // 7. Define the solution vector x as a finite element grid function
    //    corresponding to fespace. Initialize x with initial guess of zero.
-   GridFunction x(fespace);
+   GridFunction x(&fespace);
    x = 0.0;
 
-   // 7. Set up the bilinear form a(.,.) on the finite element space
+   // 8. Set up the bilinear form a(.,.) on the finite element space
    //    corresponding to the Laplacian operator -Delta, by adding the Diffusion
    //    domain integrator and the interior and boundary DG face integrators.
    //    Note that boundary conditions are imposed weakly in the form, so there
    //    is no need for dof elimination. After assembly and finalizing we
    //    extract the corresponding sparse matrix A.
-   BilinearForm *a = new BilinearForm(fespace);
-   a->AddDomainIntegrator(new DiffusionIntegrator(one));
-   a->AddInteriorFaceIntegrator(new DGDiffusionIntegrator(one, sigma, kappa));
-   a->AddBdrFaceIntegrator(new DGDiffusionIntegrator(one, sigma, kappa));
+   BilinearForm a(&fespace);
+   a.AddDomainIntegrator(new DiffusionIntegrator(one));
+   a.AddInteriorFaceIntegrator(new DGDiffusionIntegrator(one, sigma, kappa));
+   a.AddBdrFaceIntegrator(new DGDiffusionIntegrator(one, sigma, kappa));
    if (eta > 0)
    {
-      a->AddInteriorFaceIntegrator(new DGDiffusionBR2Integrator(*fespace, eta));
-      a->AddBdrFaceIntegrator(new DGDiffusionBR2Integrator(*fespace, eta));
+      MFEM_VERIFY(!pa, "BR2 not yet compatible with partial assembly.");
+      a.AddInteriorFaceIntegrator(new DGDiffusionBR2Integrator(fespace, eta));
+      a.AddBdrFaceIntegrator(new DGDiffusionBR2Integrator(fespace, eta));
    }
-   a->Assemble();
-   a->Finalize();
-   const SparseMatrix &A = a->SpMat();
+   if (pa) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
+   a.Assemble();
+   a.Finalize();
 
-#ifndef MFEM_USE_SUITESPARSE
-   // 8. Define a simple symmetric Gauss-Seidel preconditioner and use it to
+   // 9. Define a simple symmetric Gauss-Seidel preconditioner and use it to
    //    solve the system Ax=b with PCG in the symmetric case, and GMRES in the
-   //    non-symmetric one.
-   GSSmoother M(A);
-   if (sigma == -1.0)
+   //    non-symmetric one. (Note that tolerances are squared: 1e-24 corresponds
+   //    to a relative tolerance of 1e-12).
+   //
+   //    If MFEM was compiled with SuiteSparse, use UMFPACK to solve the system.
+   if (pa)
    {
-      PCG(A, M, *b, x, 1, 500, 1e-12, 0.0);
+      const Operator &A = a;
+      if (sigma == -1.0)
+      {
+         CG(A, b, x, 1, 500, 1e-24, 0.0);
+      }
+      else
+      {
+         MFEM_ABORT("The case of PA with sigma != -1 is not yet supported.");
+      }
    }
    else
    {
-      GMRES(A, M, *b, x, 1, 500, 10, 1e-12, 0.0);
-   }
+      const SparseMatrix &A = a.SpMat();
+#ifndef MFEM_USE_SUITESPARSE
+      GSSmoother M(A);
+      if (sigma == -1.0)
+      {
+         PCG(A, M, b, x, 1, 500, 1e-24, 0.0);
+      }
+      else
+      {
+         GMRES(A, M, b, x, 1, 500, 10, 1e-24, 0.0);
+      }
 #else
-   // 8. If MFEM was compiled with SuiteSparse, use UMFPACK to solve the system.
-   UMFPackSolver umf_solver;
-   umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
-   umf_solver.SetOperator(A);
-   umf_solver.Mult(*b, x);
+      UMFPackSolver umf_solver;
+      umf_solver.Control[UMFPACK_ORDERING] = UMFPACK_ORDERING_METIS;
+      umf_solver.SetOperator(A);
+      umf_solver.Mult(b, x);
 #endif
+   }
 
-   // 9. Save the refined mesh and the solution. This output can be viewed later
-   //    using GLVis: "glvis -m refined.mesh -g sol.gf".
+   // 10. Save the refined mesh and the solution. This output can be viewed
+   //     later using GLVis: "glvis -m refined.mesh -g sol.gf".
    ofstream mesh_ofs("refined.mesh");
    mesh_ofs.precision(8);
-   mesh->Print(mesh_ofs);
+   mesh.Print(mesh_ofs);
    ofstream sol_ofs("sol.gf");
    sol_ofs.precision(8);
    x.Save(sol_ofs);
 
-   // 10. Send the solution by socket to a GLVis server.
+   // 11. Send the solution by socket to a GLVis server.
    if (visualization)
    {
       char vishost[] = "localhost";
       int  visport   = 19916;
       socketstream sol_sock(vishost, visport);
       sol_sock.precision(8);
-      sol_sock << "solution\n" << *mesh << x << flush;
+      sol_sock << "solution\n" << mesh << x << flush;
    }
 
-   // 11. Free the used memory.
-   delete a;
-   delete b;
-   delete fespace;
-   delete fec;
-   delete mesh;
-
    return 0;
 }

examples/ex14p.cpp

@@ -18,6 +18,10 @@
 //               mpirun -np 4 ex14p -m ../data/amr-quad.mesh -rs 3
 //               mpirun -np 4 ex14p -m ../data/amr-hex.mesh
 //
+// Device sample runs:
+//               mpirun -np 4 ex14p -pa -d cuda -o 3
+//               mpirun -np 4 ex14p -pa -d cuda -m ../data/fichera.mesh -o 3
+//
 // Description:  This example code demonstrates the use of MFEM to define a
 //               discontinuous Galerkin (DG) finite element discretization of
 //               the Laplace problem -Delta u = 1 with homogeneous Dirichlet
@@ -38,11 +42,14 @@ using namespace mfem;
 
 class CustomSolverMonitor : public IterativeSolverMonitor
 {
+private:
+   const ParMesh &pmesh;
+   ParGridFunction &pgf;
 public:
-   CustomSolverMonitor(const ParMesh *m,
-                       ParGridFunction *f) :
-      pmesh(m),
-      pgf(f) {}
+   CustomSolverMonitor(const ParMesh &pmesh_,
+                       ParGridFunction &pgf_) :
+      pmesh(pmesh_),
+      pgf(pgf_) {}
 
    void MonitorSolution(int i, real_t norm, const Vector &x, bool final)
    {
@@ -50,30 +57,24 @@ public:
       int  visport   = 19916;
       int  num_procs, myid;
 
-      MPI_Comm_size(pmesh->GetComm(),&num_procs);
-      MPI_Comm_rank(pmesh->GetComm(),&myid);
+      MPI_Comm_size(pmesh.GetComm(), &num_procs);
+      MPI_Comm_rank(pmesh.GetComm(), &myid);
 
-      pgf->SetFromTrueDofs(x);
+      pgf.SetFromTrueDofs(x);
 
       socketstream sol_sock(vishost, visport);
       sol_sock << "parallel " << num_procs << " " << myid << "\n";
       sol_sock.precision(8);
-      sol_sock << "solution\n" << *pmesh << *pgf
+      sol_sock << "solution\n" << pmesh << pgf
                << "window_title 'Iteration no " << i << "'"
                << "keys rRjlc\n" << flush;
    }
-
-private:
-   const ParMesh *pmesh;
-   ParGridFunction *pgf;
 };
 
 int main(int argc, char *argv[])
 {
    // 1. Initialize MPI and HYPRE.
    Mpi::Init(argc, argv);
-   int num_procs = Mpi::WorldSize();
-   int myid = Mpi::WorldRank();
    Hypre::Init();
 
    // 2. Parse command-line options.
@@ -84,7 +85,9 @@ int main(int argc, char *argv[])
    real_t sigma = -1.0;
    real_t kappa = -1.0;
    real_t eta = 0.0;
+   bool pa = false;
    bool visualization = 1;
+   const char *device_config = "cpu";
 
    OptionsParser args(argc, argv);
    args.AddOption(&mesh_file, "-m", "--mesh",
@@ -103,13 +106,17 @@ int main(int argc, char *argv[])
                   "One of the three DG penalty parameters, should be positive."
                   " Negative values are replaced with (order+1)^2.");
    args.AddOption(&eta, "-e", "--eta", "BR2 penalty parameter.");
+   args.AddOption(&pa, "-pa", "--partial-assembly", "-no-pa",
+                  "--no-partial-assembly", "Enable Partial Assembly.");
    args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
                   "--no-visualization",
                   "Enable or disable GLVis visualization.");
+   args.AddOption(&device_config, "-d", "--device",
+                  "Device configuration string, see Device::Configure().");
    args.Parse();
    if (!args.Good())
    {
-      if (myid == 0)
+      if (Mpi::Root())
       {
          args.PrintUsage(cout);
       }
@@ -119,16 +126,19 @@ int main(int argc, char *argv[])
    {
       kappa = (order+1)*(order+1);
    }
-   if (myid == 0)
+   if (Mpi::Root())
    {
       args.PrintOptions(cout);
    }
 
+   Device device(device_config);
+   if (Mpi::Root()) { device.Print(); }
+
    // 3. 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. NURBS meshes are projected to second order meshes.
-   Mesh *mesh = new Mesh(mesh_file, 1, 1);
-   int dim = mesh->Dimension();
+   Mesh mesh(mesh_file);
+   int dim = mesh.Dimension();
 
    // 4. Refine the serial mesh on all processors to increase the resolution. In
    //    this example we do 'ser_ref_levels' of uniform refinement. By default,
@@ -137,53 +147,53 @@ int main(int argc, char *argv[])
    {
       if (ser_ref_levels < 0)
       {
-         ser_ref_levels = (int)floor(log(10000./mesh->GetNE())/log(2.)/dim);
+         ser_ref_levels = (int)floor(log(10000./mesh.GetNE())/log(2.)/dim);
       }
       for (int l = 0; l < ser_ref_levels; l++)
       {
-         mesh->UniformRefinement();
+         mesh.UniformRefinement();
       }
    }
-   if (mesh->NURBSext)
+   if (mesh.NURBSext)
    {
-      mesh->SetCurvature(max(order, 1));
+      mesh.SetCurvature(max(order, 1));
    }
 
    // 5. 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.
-   ParMesh *pmesh = new ParMesh(MPI_COMM_WORLD, *mesh);
-   delete mesh;
+   ParMesh pmesh(MPI_COMM_WORLD, mesh);
+   mesh.Clear();
    {
       for (int l = 0; l < par_ref_levels; l++)
       {
-         pmesh->UniformRefinement();
+         pmesh.UniformRefinement();
       }
    }
 
    // 6. Define a parallel finite element space on the parallel mesh. Here we
    //    use discontinuous finite elements of the specified order >= 0.
-   FiniteElementCollection *fec = new DG_FECollection(order, dim);
-   ParFiniteElementSpace *fespace = new ParFiniteElementSpace(pmesh, fec);
-   HYPRE_BigInt size = fespace->GlobalTrueVSize();
-   if (myid == 0)
+   DG_FECollection fec(order, dim, BasisType::GaussLobatto);
+   ParFiniteElementSpace fespace(&pmesh, &fec);
+   HYPRE_BigInt size = fespace.GlobalTrueVSize();
+   if (Mpi::Root())
    {
       cout << "Number of unknowns: " << size << endl;
    }
 
    // 7. Set up the parallel linear form b(.) which corresponds to the
    //    right-hand side of the FEM linear system.
-   ParLinearForm *b = new ParLinearForm(fespace);
+   ParLinearForm b(&fespace);
    ConstantCoefficient one(1.0);
    ConstantCoefficient zero(0.0);
-   b->AddDomainIntegrator(new DomainLFIntegrator(one));
-   b->AddBdrFaceIntegrator(
+   b.AddDomainIntegrator(new DomainLFIntegrator(one));
+   b.AddBdrFaceIntegrator(
       new DGDirichletLFIntegrator(zero, one, sigma, kappa));
-   b->Assemble();
+   b.Assemble();
 
    // 8. Define the solution vector x as a parallel finite element grid function
    //    corresponding to fespace. Initialize x with initial guess of zero.
-   ParGridFunction x(fespace);
+   ParGridFunction x(&fespace);
    x = 0.0;
 
    // 9. Set up the bilinear form a(.,.) on the finite element space
@@ -192,42 +202,51 @@ int main(int argc, char *argv[])
    //    Note that boundary conditions are imposed weakly in the form, so there
    //    is no need for dof elimination. After serial and parallel assembly we
    //    extract the corresponding parallel matrix A.
-   ParBilinearForm *a = new ParBilinearForm(fespace);
-   a->AddDomainIntegrator(new DiffusionIntegrator(one));
-   a->AddInteriorFaceIntegrator(new DGDiffusionIntegrator(one, sigma, kappa));
-   a->AddBdrFaceIntegrator(new DGDiffusionIntegrator(one, sigma, kappa));
+   ParBilinearForm a(&fespace);
+   a.AddDomainIntegrator(new DiffusionIntegrator(one));
+   a.AddInteriorFaceIntegrator(new DGDiffusionIntegrator(one, sigma, kappa));
+   a.AddBdrFaceIntegrator(new DGDiffusionIntegrator(one, sigma, kappa));
    if (eta > 0)
    {
-      a->AddInteriorFaceIntegrator(new DGDiffusionBR2Integrator(*fespace, eta));
-      a->AddBdrFaceIntegrator(new DGDiffusionBR2Integrator(*fespace, eta));
+      MFEM_VERIFY(!pa, "BR2 not yet compatible with partial assembly.");
+      a.AddInteriorFaceIntegrator(new DGDiffusionBR2Integrator(fespace, eta));
+      a.AddBdrFaceIntegrator(new DGDiffusionBR2Integrator(fespace, eta));
    }
-   a->Assemble();
-   a->Finalize();
+   if (pa) { a.SetAssemblyLevel(AssemblyLevel::PARTIAL); }
+   a.Assemble();
+   a.Finalize();
 
    // 10. Define the parallel (hypre) matrix and vectors representing a(.,.),
    //     b(.) and the finite element approximation.
-   HypreParMatrix *A = a->ParallelAssemble();
-   HypreParVector *B = b->ParallelAssemble();
-   HypreParVector *X = x.ParallelProject();
+   OperatorHandle A;
 
-   delete a;
-   delete b;
-
-   // 11. Depending on the symmetry of A, define and apply a parallel PCG or
-   //     GMRES solver for AX=B using the BoomerAMG preconditioner from hypre.
-   HypreSolver *amg = new HypreBoomerAMG(*A);
-   if (sigma == -1.0)
+   std::unique_ptr<HypreBoomerAMG> amg;
+   if (pa)
    {
-      HyprePCG pcg(*A);
-      pcg.SetTol(1e-12);
-      pcg.SetMaxIter(500);
-      pcg.SetPrintLevel(2);
-      pcg.SetPreconditioner(*amg);
-      pcg.Mult(*B, *X);
+      A.Reset(&a, false);
    }
    else
    {
-      CustomSolverMonitor monitor(pmesh, &x);
+      A.SetType(Operator::Hypre_ParCSR);
+      a.ParallelAssemble(A);
+      amg.reset(new HypreBoomerAMG(*A.As<HypreParMatrix>()));
+   }
+
+   // 11. Depending on the symmetry of A, define and apply a parallel PCG or
+   //     GMRES solver for AX=B using the BoomerAMG preconditioner from hypre.
+   if (sigma == -1.0)
+   {
+      CGSolver cg(MPI_COMM_WORLD);
+      cg.SetRelTol(1e-12);
+      cg.SetMaxIter(500);
+      cg.SetPrintLevel(1);
+      cg.SetOperator(*A);
+      if (amg) { cg.SetPreconditioner(*amg); }
+      cg.Mult(b, x);
+   }
+   else
+   {
+      CustomSolverMonitor monitor(pmesh, x);
       GMRESSolver gmres(MPI_COMM_WORLD);
       gmres.SetAbsTol(0.0);
       gmres.SetRelTol(1e-12);
@@ -235,51 +254,37 @@ int main(int argc, char *argv[])
       gmres.SetKDim(10);
       gmres.SetPrintLevel(1);
       gmres.SetOperator(*A);
-      gmres.SetPreconditioner(*amg);
+      if (amg) { gmres.SetPreconditioner(*amg); }
       gmres.SetMonitor(monitor);
-      gmres.Mult(*B, *X);
+      gmres.Mult(b, x);
    }
-   delete amg;
 
-   // 12. Extract the parallel grid function corresponding to the finite element
-   //     approximation X. This is the local solution on each processor.
-   x = *X;
-
-   // 13. Save the refined mesh and the solution in parallel. This output can
+   // 12. Save the refined mesh and the solution in parallel. This output can
    //     be viewed later using GLVis: "glvis -np <np> -m mesh -g sol".
    {
       ostringstream mesh_name, sol_name;
-      mesh_name << "mesh." << setfill('0') << setw(6) << myid;
-      sol_name << "sol." << setfill('0') << setw(6) << myid;
+      mesh_name << "mesh." << setfill('0') << setw(6) << Mpi::WorldRank();
+      sol_name << "sol." << setfill('0') << setw(6) << Mpi::WorldRank();
 
       ofstream mesh_ofs(mesh_name.str().c_str());
       mesh_ofs.precision(8);
-      pmesh->Print(mesh_ofs);
+      pmesh.Print(mesh_ofs);
 
       ofstream sol_ofs(sol_name.str().c_str());
       sol_ofs.precision(8);
       x.Save(sol_ofs);
    }
 
-   // 14. Send the solution by socket to a GLVis server.
+   // 13. Send the solution by socket to a GLVis server.
    if (visualization)
    {
       char vishost[] = "localhost";
       int  visport   = 19916;
       socketstream sol_sock(vishost, visport);
-      sol_sock << "parallel " << num_procs << " " << myid << "\n";
+      sol_sock << "parallel " << Mpi::WorldSize() << " " << Mpi::WorldRank() << "\n";
       sol_sock.precision(8);
-      sol_sock << "solution\n" << *pmesh << x << flush;
+      sol_sock << "solution\n" << pmesh << x << flush;
    }
 
-   // 15. Free the used memory.
-   delete X;
-   delete B;
-   delete A;
-
-   delete fespace;
-   delete fec;
-   delete pmesh;
-
    return 0;
 }

fem/CMakeLists.txt

@@ -18,6 +18,7 @@ set(SRCS
   integ/bilininteg_convection_pa.cpp
   integ/bilininteg_convection_ea.cpp
   integ/bilininteg_curlcurl_pa.cpp
+  integ/bilininteg_dgdiffusion_pa.cpp
   integ/bilininteg_dgtrace_pa.cpp
   integ/bilininteg_dgtrace_ea.cpp
   integ/bilininteg_diffusion_mf.cpp
@@ -117,6 +118,7 @@ set(SRCS
   quadinterpolator.cpp
   quadinterpolator_face.cpp
   restriction.cpp
+  normal_deriv_restriction.cpp
   staticcond.cpp
   tmop.cpp
   tmop/tmop_pa.cpp
@@ -228,6 +230,7 @@ set(HDRS
   quadinterpolator.hpp
   quadinterpolator_face.hpp
   restriction.hpp
+  normal_deriv_restriction.hpp
   fespacehierarchy.hpp
   staticcond.hpp
   tbilinearform.hpp

fem/bilinearform_ext.cpp

@@ -282,6 +282,22 @@ void PABilinearFormExtension::SetupRestrictionOperators(const L2FaceValues m)
       int_face_X.SetSize(int_face_restrict_lex->Height(), Device::GetMemoryType());
       int_face_Y.SetSize(int_face_restrict_lex->Height(), Device::GetMemoryType());
       int_face_Y.UseDevice(true); // ensure 'int_face_Y = 0.0' is done on device
+
+      bool needs_normal_derivs = false;
+      auto &integs = *a->GetFBFI();
+      for (int i = 0; i < integs.Size(); ++i)
+      {
+         if (integs[i]->RequiresFaceNormalDerivatives())
+         {
+            needs_normal_derivs = true;
+            break;
+         }
+      }
+      if (needs_normal_derivs)
+      {
+         int_face_dXdn.SetSize(int_face_restrict_lex->Height());
+         int_face_dYdn.SetSize(int_face_restrict_lex->Height());
+      }
    }
 
    const bool has_bdr_integs = (a->GetBFBFI()->Size() > 0 ||
@@ -296,6 +312,22 @@ void PABilinearFormExtension::SetupRestrictionOperators(const L2FaceValues m)
       bdr_face_Y.SetSize(bdr_face_restrict_lex->Height(), Device::GetMemoryType());
       bdr_face_Y.UseDevice(true); // ensure 'faceBoundY = 0.0' is done on device
 
+      bool needs_normal_derivs = false;
+      auto &integs = *a->GetBFBFI();
+      for (int i = 0; i < integs.Size(); ++i)
+      {
+         if (integs[i]->RequiresFaceNormalDerivatives())
+         {
+            needs_normal_derivs = true;
+            break;
+         }
+      }
+      if (needs_normal_derivs)
+      {
+         bdr_face_dXdn.SetSize(bdr_face_restrict_lex->Height());
+         bdr_face_dYdn.SetSize(bdr_face_restrict_lex->Height());
+      }
+
       const Mesh &mesh = *trial_fes->GetMesh();
       // See LinearFormExtension::Update for explanation of f_to_be logic.
       std::unordered_map<int,int> f_to_be;
@@ -542,8 +574,8 @@ void PABilinearFormExtension::Mult(const Vector &x, Vector &y) const
          localY = 0.0;
          for (int i = 0; i < iSz; ++i)
          {
-            AddMultWithMarkers(*integrators[i], localX, elem_markers[i], elem_attributes,
-                               false, localY);
+            AddMultWithMarkers(*integrators[i], localX, elem_markers[i],
+                               elem_attributes, false, localY);
          }
          elem_restrict->MultTranspose(localY, y);
       }
@@ -558,14 +590,39 @@ void PABilinearFormExtension::Mult(const Vector &x, Vector &y) const
    if (int_face_restrict_lex && iFISz>0)
    {
       int_face_restrict_lex->Mult(x, int_face_X);
-      if (int_face_X.Size()>0)
+      if (int_face_dXdn.Size() > 0)
+      {
+         int_face_restrict_lex->NormalDerivativeMult(x, int_face_dXdn);
+      }
+      if (int_face_X.Size() > 0)
       {
          int_face_Y = 0.0;
+
+         // if normal derivatives are needed by at least one integrator...
+         if (int_face_dYdn.Size() > 0)
+         {
+            int_face_dYdn = 0.0;
+         }
+
          for (int i = 0; i < iFISz; ++i)
          {
-            intFaceIntegrators[i]->AddMultPA(int_face_X, int_face_Y);
+            if (intFaceIntegrators[i]->RequiresFaceNormalDerivatives())
+            {
+               intFaceIntegrators[i]->AddMultPAFaceNormalDerivatives(
+                  int_face_X, int_face_dXdn,
+                  int_face_Y, int_face_dYdn);
+            }
+            else
+            {
+               intFaceIntegrators[i]->AddMultPA(int_face_X, int_face_Y);
+            }
          }
          int_face_restrict_lex->AddMultTransposeInPlace(int_face_Y, y);
+         if (int_face_dYdn.Size() > 0)
+         {
+            int_face_restrict_lex->NormalDerivativeAddMultTranspose(
+               int_face_dYdn, y);
+         }
       }
    }
 
@@ -579,9 +636,19 @@ void PABilinearFormExtension::Mult(const Vector &x, Vector &y) const
       Array<Array<int>*> &bdr_markers = *a->GetBBFI_Marker();
       Array<Array<int>*> &bdr_face_markers = *a->GetBFBFI_Marker();
       bdr_face_restrict_lex->Mult(x, bdr_face_X);
-      if (bdr_face_X.Size()>0)
+      if (bdr_face_dXdn.Size() > 0)
+      {
+         bdr_face_restrict_lex->NormalDerivativeMult(x, bdr_face_dXdn);
+      }
+      if (bdr_face_X.Size() > 0)
       {
          bdr_face_Y = 0.0;
+
+         // if normal derivatives are needed by at least one integrator...
+         if (bdr_face_dYdn.Size() > 0)
+         {
+            bdr_face_dYdn = 0.0;
+         }
          for (int i = 0; i < n_bdr_integs; ++i)
          {
             AddMultWithMarkers(*bdr_integs[i], bdr_face_X, bdr_markers[i], bdr_attributes,
@@ -589,10 +656,23 @@ void PABilinearFormExtension::Mult(const Vector &x, Vector &y) const
          }
          for (int i = 0; i < n_bdr_face_integs; ++i)
          {
-            AddMultWithMarkers(*bdr_face_integs[i], bdr_face_X, bdr_face_markers[i],
-                               bdr_attributes, false, bdr_face_Y);
+            if (bdr_face_integs[i]->RequiresFaceNormalDerivatives())
+            {
+               AddMultNormalDerivativesWithMarkers(
+                  *bdr_face_integs[i], bdr_face_X, bdr_face_dXdn,
+                  bdr_face_markers[i], bdr_attributes, bdr_face_Y, bdr_face_dYdn);
+            }
+            else
+            {
+               AddMultWithMarkers(*bdr_face_integs[i], bdr_face_X, bdr_face_markers[i],
+                                  bdr_attributes, false, bdr_face_Y);
+            }
          }
          bdr_face_restrict_lex->AddMultTransposeInPlace(bdr_face_Y, y);
+         if (bdr_face_dYdn.Size() > 0)
+         {
+            bdr_face_restrict_lex->NormalDerivativeAddMultTranspose(bdr_face_dYdn, y);
+         }
       }
    }
 }
@@ -693,6 +773,37 @@ static void AddWithMarkers_(
    });
 }
 
+void PABilinearFormExtension::AddMultNormalDerivativesWithMarkers(
+   const BilinearFormIntegrator &integ,
+   const Vector &x,
+   const Vector &dxdn,
+   const Array<int> *markers,
+   const Array<int> &attributes,
+   Vector &y,
+   Vector &dydn) const
+{
+   if (markers)
+   {
+      tmp_evec.SetSize(y.Size() + dydn.Size());
+      tmp_evec = 0.0;
+      Vector tmp_y(tmp_evec, 0, y.Size());
+      Vector tmp_dydn(tmp_evec, y.Size(), dydn.Size());
+
+      integ.AddMultPAFaceNormalDerivatives(x, dxdn, tmp_y, tmp_dydn);
+
+      const int ne = attributes.Size();
+      const int nd_1 = x.Size() / ne;
+      const int nd_2 = dxdn.Size() / ne;
+
+      AddWithMarkers_(ne, nd_1, tmp_y, *markers, attributes, y);
+      AddWithMarkers_(ne, nd_2, tmp_dydn, *markers, attributes, dydn);
+   }
+   else
+   {
+      integ.AddMultPAFaceNormalDerivatives(x, dxdn, y, dydn);
+   }
+}
+
 void PABilinearFormExtension::AddMultWithMarkers(
    const BilinearFormIntegrator &integ,
    const Vector &x,

fem/bilinearform_ext.hpp

@@ -74,6 +74,8 @@ protected:
    mutable Vector localX, localY;
    mutable Vector int_face_X, int_face_Y;
    mutable Vector bdr_face_X, bdr_face_Y;
+   mutable Vector int_face_dXdn, int_face_dYdn;
+   mutable Vector bdr_face_dXdn, bdr_face_dYdn;
    const Operator *elem_restrict; // Not owned
    const FaceRestriction *int_face_restrict_lex; // Not owned
    const FaceRestriction *bdr_face_restrict_lex; // Not owned
@@ -113,6 +115,23 @@ protected:
                            const Array<int> &attributes,
                            const bool transpose,
                            Vector &y) const;
+
+   /// @brief Performs the same function as AddMultWithMarkers, but takes as
+   /// input and output face normal derivatives.
+   ///
+   /// This is required when the integrator requires face normal derivatives,
+   /// for example, DGDiffusionIntegrator.
+   ///
+   /// This is called when the integrator's member function
+   /// BilinearFormIntegrator::RequiresFaceNormalDerivatives() returns true.
+   void AddMultNormalDerivativesWithMarkers(
+      const BilinearFormIntegrator &integ,
+      const Vector &x,
+      const Vector &dxdn,
+      const Array<int> *markers,
+      const Array<int> &attributes,
+      Vector &y,
+      Vector &dydn) const;
 };
 
 /// Data and methods for element-assembled bilinear forms

fem/bilininteg.cpp

@@ -189,6 +189,12 @@ void BilinearFormIntegrator::AssembleTraceFaceMatrix (int elem,
               " Integrator class.");
 }
 
+void BilinearFormIntegrator::AddMultPAFaceNormalDerivatives(
+   const Vector &x, const Vector &dxdn, Vector &y, Vector &dydn) const
+{
+   MFEM_ABORT("Not implemented.");
+}
+
 void BilinearFormIntegrator::AssembleElementVector(
    const FiniteElement &el, ElementTransformation &Tr, const Vector &elfun,
    Vector &elvect)
@@ -3423,7 +3429,7 @@ void DGDiffusionIntegrator::AssembleFaceMatrix(
    const FiniteElement &el1, const FiniteElement &el2,
    FaceElementTransformations &Trans, DenseMatrix &elmat)
 {
-   int dim, ndof1, ndof2, ndofs;
+   int ndof1, ndof2, ndofs;
    bool kappa_is_nonzero = (kappa != 0.);
    real_t w, wq = 0.0;
 
@@ -3466,17 +3472,9 @@ void DGDiffusionIntegrator::AssembleFaceMatrix(
    const IntegrationRule *ir = IntRule;
    if (ir == NULL)
    {
-      // a simple choice for the integration order; is this OK?
-      int order;
-      if (ndof2)
-      {
-         order = 2*max(el1.GetOrder(), el2.GetOrder());
-      }
-      else
-      {
-         order = 2*el1.GetOrder();
-      }
-      ir = &IntRules.Get(Trans.GetGeometryType(), order);
+      const int order = (ndof2) ? max(el1.GetOrder(),
+                                      el2.GetOrder()) : el1.GetOrder();
+      ir = &GetRule(order, Trans);
    }
 
    // assemble: < {(Q \nabla u).n},[v] >      --> elmat
@@ -3654,6 +3652,13 @@ void DGDiffusionIntegrator::AssembleFaceMatrix(
    }
 }
 
+const IntegrationRule &DGDiffusionIntegrator::GetRule(
+   int order, FaceElementTransformations &T)
+{
+   // order is typically the maximum of the order of the left and right elements
+   // neighboring the given face.
+   return IntRules.Get(T.GetGeometryType(), 2*order);
+}
 
 // static method
 void DGElasticityIntegrator::AssembleBlock(

fem/bilininteg.hpp

@@ -266,6 +266,39 @@ public:
                                     Vector &flux, Vector *d_energy = NULL)
    { return 0.0; }
 
+   /** @brief For bilinear forms on element faces, specifies if the normal
+              derivatives are needed on the faces or just the face restriction.
+
+       @details if RequiresFaceNormalDerivatives() == true, then
+                AddMultPAFaceNormalDerivatives(...) should be invoked in place
+                of AddMultPA(...) and L2NormalDerivativeFaceRestriction should
+                be used to compute the normal derivatives. This is used for some
+                DG integrators, for example DGDiffusionIntegrator.
+
+       @returns whether normal derivatives appear in the bilinear form.
+   */
+   virtual bool RequiresFaceNormalDerivatives() const { return false; }
+
+   /// Method for partially assembled action.
+   /** @brief For bilinear forms on element faces that depend on the normal
+              derivative on the faces, computes the action of integrator to the
+              face values @a x and reference-normal derivatives @a dxdn and adds
+              the result to @a y and @a dydn.
+
+      @details This method can be called only after the method AssemblePA() has
+               been called.
+
+      @param[in]     x E-vector of face values (provided by
+                       FaceRestriction::Mult)
+      @param[in]     dxdn E-vector of face reference-normal derivatives
+                          (provided by FaceRestriction::NormalDerivativeMult)
+      @param[in,out] y E-vector of face values to add action to.
+      @param[in,out] dydn E-vector of face reference-normal derivative values to
+                          add action to.
+   */
+   virtual void AddMultPAFaceNormalDerivatives(const Vector &x, const Vector &dxdn,
+                                               Vector &y, Vector &dydn) const;
+
    virtual ~BilinearFormIntegrator() { }
 };
 
@@ -3229,6 +3262,13 @@ protected:
    Vector shape1, shape2, dshape1dn, dshape2dn, nor, nh, ni;
    DenseMatrix jmat, dshape1, dshape2, mq, adjJ;
 
+
+   // PA extension
+   Vector pa_data; // (Q, h, dot(n,J)|el0, dot(n,J)|el1)
+   const DofToQuad *maps; ///< Not owned
+   int dim, nf, nq, dofs1D, quad1D;
+   IntegrationRules irs{0, Quadrature1D::GaussLobatto};
+
 public:
    DGDiffusionIntegrator(const real_t s, const real_t k)
       : Q(NULL), MQ(NULL), sigma(s), kappa(k) { }
@@ -3237,10 +3277,26 @@ public:
    DGDiffusionIntegrator(MatrixCoefficient &q, const real_t s, const real_t k)
       : Q(NULL), MQ(&q), sigma(s), kappa(k) { }
    using BilinearFormIntegrator::AssembleFaceMatrix;
-   virtual void AssembleFaceMatrix(const FiniteElement &el1,
-                                   const FiniteElement &el2,
-                                   FaceElementTransformations &Trans,
-                                   DenseMatrix &elmat);
+   void AssembleFaceMatrix(const FiniteElement &el1,
+                           const FiniteElement &el2,
+                           FaceElementTransformations &Trans,
+                           DenseMatrix &elmat) override;
+
+   bool RequiresFaceNormalDerivatives() const override { return true; }
+
+   using BilinearFormIntegrator::AssemblePA;
+
+   void AssemblePAInteriorFaces(const FiniteElementSpace &fes) override;
+
+   void AssemblePABoundaryFaces(const FiniteElementSpace &fes) override;
+
+   void AddMultPAFaceNormalDerivatives(const Vector &x, const Vector &dxdn,
+                                       Vector &y, Vector &dydn) const override;
+
+   const IntegrationRule &GetRule(int order, FaceElementTransformations &T);
+
+private:
+   void SetupPA(const FiniteElementSpace &fes, FaceType type);
 };
 
 /** Integrator for the "BR2" diffusion stabilization term

fem/fe/face_map_utils.hpp

@@ -13,6 +13,7 @@
 #define MFEM_FACE_MAP_UTILS_HPP
 
 #include "../../general/array.hpp"
+#include "../../general/backends.hpp"
 #include <utility> // std::pair
 #include <vector>
 
@@ -51,6 +52,189 @@ void FillFaceMap(const int n_face_dofs_per_component,
 void GetTensorFaceMap(const int dim, const int order, const int face_id,
                       Array<int> &face_map);
 
+/// @brief Given a face DOF index in native (counter-clockwise) ordering, return
+/// the corresponding DOF index in lexicographic ordering (for a quadrilateral
+/// element).
+MFEM_HOST_DEVICE
+inline int ToLexOrdering2D(const int face_id, const int size1d, const int i)
+{
+   if (face_id==2 || face_id==3)
+   {
+      return size1d-1-i;
+   }
+   else
+   {
+      return i;
+   }
+}
+
+/// @brief Given a face DOF index on a shared face, ordered lexicographically
+/// relative to element 1, return the corresponding face DOF index ordered
+/// lexicographically relative to element 2.
+MFEM_HOST_DEVICE
+inline int PermuteFace2D(const int face_id1, const int face_id2,
+                         const int orientation, const int size1d,
+                         const int index)
+{
+   int new_index;
+   // Convert from element 1 lex ordering to native ordering
+   if (face_id1 == 2 || face_id1 == 3)
+   {
+      new_index = size1d-1-index;
+   }
+   else
+   {
+      new_index = index;
+   }
+   // Permute based on face orientations
+   if (orientation == 1)
+   {
+      new_index = size1d-1-new_index;
+   }
+   // Covert to element 2 lex ordering
+   return ToLexOrdering2D(face_id2, size1d, new_index);
+}
+
+/// @brief Given a face DOF index in native (counter-clockwise) ordering, return
+/// the corresponding DOF index in lexicographic ordering (for a hexahedral
+/// element).
+MFEM_HOST_DEVICE
+inline int ToLexOrdering3D(const int face_id, const int size1d, const int i,
+                           const int j)
+{
+   if (face_id==2 || face_id==1 || face_id==5)
+   {
+      return i + j*size1d;
+   }
+   else if (face_id==3 || face_id==4)
+   {
+      return (size1d-1-i) + j*size1d;
+   }
+   else // face_id==0
+   {
+      return i + (size1d-1-j)*size1d;
+   }
+}
+
+/// @brief Given the index of a face DOF in lexicographic ordering relative
+/// element 1, permute the index so that it is lexicographically ordered
+/// relative to element 2.
+///
+/// The given face corresponds to local face index @a face_id1 relative to
+/// element 1, and @a face_id2 (with @a orientation) relative to element 2.
+MFEM_HOST_DEVICE
+inline int PermuteFace3D(const int face_id1, const int face_id2,
+                         const int orientation,
+                         const int size1d, const int index)
+{
+   int i=0, j=0, new_i=0, new_j=0;
+   i = index%size1d;
+   j = index/size1d;
+   // Convert from lex ordering
+   if (face_id1==3 || face_id1==4)
+   {
+      i = size1d-1-i;
+   }
+   else if (face_id1==0)
+   {
+      j = size1d-1-j;
+   }
+   // Permute based on face orientations
+   switch (orientation)
+   {
+      case 0:
+         new_i = i;
+         new_j = j;
+         break;
+      case 1:
+         new_i = j;
+         new_j = i;
+         break;
+      case 2:
+         new_i = j;
+         new_j = (size1d-1-i);
+         break;
+      case 3:
+         new_i = (size1d-1-i);
+         new_j = j;
+         break;
+      case 4:
+         new_i = (size1d-1-i);
+         new_j = (size1d-1-j);
+         break;
+      case 5:
+         new_i = (size1d-1-j);
+         new_j = (size1d-1-i);
+         break;
+      case 6:
+         new_i = (size1d-1-j);
+         new_j = i;
+         break;
+      case 7:
+         new_i = i;
+         new_j = (size1d-1-j);
+         break;
+   }
+   return ToLexOrdering3D(face_id2, size1d, new_i, new_j);
+}
+
+/// @brief Given a face DOF (or quadrature) index ordered lexicographically
+/// relative to element 1, return the associated (i, j) coordinates.
+///
+/// The returned coordinates will be relative to element 1 or element 2
+/// according to the value of side (side == 0 corresponds element 1).
+MFEM_HOST_DEVICE
+inline void FaceIdxToVolIdx2D(const int qi, const int nq, const int face_id0,
+                              const int face_id1, const int side, int &i, int &j)
+{
+   // Note: in 2D, a consistently ordered mesh will always have the element 2
+   // face reversed relative to element 1, so orientation is determined entirely
+   // by side. (In 3D, separate orientation information is needed).
+   const int orientation = side;
+
+   const int face_id = (side == 0) ? face_id0 : face_id1;
+   const int edge_idx = (side == 0) ? qi : PermuteFace2D(face_id0, face_id1,
+                                                         orientation, nq, qi);
+
+   const int level = (face_id == 0 || face_id == 3) ? 0 : (nq-1);
+   const bool x_axis = (face_id == 0 || face_id == 2);
+
+   i = x_axis ? edge_idx : level;
+   j = x_axis ? level : edge_idx;
+}
+
+/// @brief Given a face DOF (or quadrature) index ordered lexicographically
+/// relative to element 1, return the associated (i, j, k) coordinates.
+///
+/// The returned coordinates will be relative to element 1 or element 2
+/// according to the value of side (side == 0 corresponds element 1).
+MFEM_HOST_DEVICE
+inline void FaceIdxToVolIdx3D(const int index, const int size1d,
+                              const int face_id0, const int face_id1,
+                              const int side, const int orientation,
+                              int& i, int& j, int& k)
+{
+   MFEM_VERIFY_KERNEL(face_id1 >= 0 || side == 0,
+                      "Accessing second side but face_id1 is not valid.");
+
+   const int face_id = (side == 0) ? face_id0 : face_id1;
+   const int fidx = (side == 0) ? index
+                    : PermuteFace3D(face_id0, face_id1, orientation, size1d, index);
+
+   const bool xy_plane = (face_id == 0 || face_id == 5);
+   const bool yz_plane = (face_id == 2 || face_id == 4);
+
+   const int level = (face_id == 0 || face_id == 1 || face_id == 4)
+                     ? 0 : (size1d-1);
+
+   const int _i = fidx % size1d;
+   const int _j = fidx / size1d;
+
+   k = xy_plane ? level : _j;
+   j = yz_plane ? _i : xy_plane ? _j : level;
+   i = yz_plane ? level : _i;
+}
+
 } // namespace internal
 
 } // namespace mfem

fem/fespace.cpp

@@ -3201,6 +3201,16 @@ const FiniteElement *FiniteElementSpace::GetFE(int i) const
    return FE;
 }
 
+const FiniteElement *FiniteElementSpace::GetTypicalFE() const
+{
+   if (mesh->GetNE() > 0) { return GetFE(0); }
+
+   Geometry::Type geom = mesh->GetTypicalElementGeometry();
+   const FiniteElement *fe = fec->FiniteElementForGeometry(geom);
+   MFEM_VERIFY(fe != nullptr, "Could not determine a typical FE!");
+   return fe;
+}
+
 const FiniteElement *FiniteElementSpace::GetBE(int i) const
 {
    int order = fec->GetOrder();

fem/fespace.hpp

@@ -1168,6 +1168,13 @@ public:
         an empty partition. */
    virtual const FiniteElement *GetFE(int i) const;
 
+   /** @brief Return GetFE(0) if the local mesh is not empty; otherwise return a
+       typical FE based on the Geometry types in the global mesh.
+
+       This method can be used as a replacement for GetFE(0) that will be valid
+       even if the local mesh is empty. */
+   const FiniteElement *GetTypicalFE() const;
+
    /** @brief Returns pointer to the FiniteElement in the FiniteElementCollection
         associated with i'th boundary face in the mesh object. */
    const FiniteElement *GetBE(int i) const;
@@ -1341,18 +1348,19 @@ public:
    virtual ~FiniteElementSpace();
 };
 
-/// @brief Return true if the mesh contains only one topology and the elements are tensor elements.
+/// @brief Return true if the mesh contains only one topology and the elements
+/// are tensor elements.
 inline bool UsesTensorBasis(const FiniteElementSpace& fes)
 {
    Mesh & mesh = *fes.GetMesh();
    const bool mixed = mesh.GetNumGeometries(mesh.Dimension()) > 1;
-   // Potential issue: empty local mesh --> no element 0.
    return !mixed &&
-          dynamic_cast<const mfem::TensorBasisElement *>(fes.GetFE(0))!=nullptr;
+          dynamic_cast<const mfem::TensorBasisElement *>(
+             fes.GetTypicalFE()) != nullptr;
 }
 
-/// @brief Return LEXICOGRAPHIC if mesh contains only one topology and the elements are tensor
-/// elements, otherwise, return NATIVE.
+/// @brief Return LEXICOGRAPHIC if mesh contains only one topology and the
+/// elements are tensor elements, otherwise, return NATIVE.
 ElementDofOrdering GetEVectorOrdering(const FiniteElementSpace& fes);
 
 }

fem/integ/bilininteg_diffusion_pa.cpp

@@ -26,8 +26,7 @@ void DiffusionIntegrator::AssemblePA(const FiniteElementSpace &fes)
    // Assuming the same element type
    fespace = &fes;
    Mesh *mesh = fes.GetMesh();
-   if (mesh->GetNE() == 0) { return; }
-   const FiniteElement &el = *fes.GetFE(0);
+   const FiniteElement &el = *fes.GetTypicalFE();
    const IntegrationRule *ir = IntRule ? IntRule : &GetRule(el, el);
    if (DeviceCanUseCeed())
    {

fem/normal_deriv_restriction.cpp

@@ -0,0 +1,858 @@
+// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "normal_deriv_restriction.hpp"
+#include "fespace.hpp"
+#include "pgridfunc.hpp"
+#include "fe/face_map_utils.hpp"
+#include "../general/forall.hpp"
+
+namespace mfem
+{
+
+/// Compute the face index to volume index map "face_to_vol" in 2D
+static void NormalDerivativeSetupFaceIndexMap2D(
+   int nf, int d, const Array<int>& face_to_elem, Array<int>& face_to_vol)
+{
+   const auto f2e = Reshape(face_to_elem.HostRead(), 2, 2, nf);
+   auto f2v = Reshape(face_to_vol.HostWrite(), d, 2, nf);
+
+   for (int f = 0; f < nf; ++f)
+   {
+      const int fid0 = f2e(0, 1, f);
+      const int fid1 = f2e(1, 1, f);
+      for (int side = 0; side < 2; ++side)
+      {
+         const int el = f2e(side, 0, f);
+
+         if (el < 0)
+         {
+            for (int p = 0; p < d; ++p)
+            {
+               f2v(p, side, f) = -1;
+            }
+         }
+         else
+         {
+            for (int p = 0; p < d; ++p)
+            {
+               int i, j;
+               internal::FaceIdxToVolIdx2D(p, d, fid0, fid1, side, i, j);
+
+               f2v(p, side, f) = i + d * j;
+            }
+         }
+      }
+   }
+}
+
+/// Compute the face index to volume index map "face_to_vol" in 3D
+static void NormalDerivativeSetupFaceIndexMap3D(
+   int nf, int d, const Array<int>& face_to_elem, Array<int>& face_to_vol)
+{
+   const auto f2e = Reshape(face_to_elem.HostRead(), 2, 3, nf);
+   auto f2v = Reshape(face_to_vol.HostWrite(), d*d, 2, nf);
+
+   for (int f = 0; f < nf; ++f)
+   {
+      const int fid0 = f2e(0, 1, f);
+      const int fid1 = f2e(1, 1, f);
+      for (int side = 0; side < 2; ++side)
+      {
+         const int el = f2e(side, 0, f);
+         const int orientation = f2e(side, 2, f);
+
+         if (el < 0)
+         {
+            for (int p = 0; p < d*d; ++p)
+            {
+               f2v(p, side, f) = -1;
+            }
+         }
+         else
+         {
+            for (int p = 0; p < d*d; ++p)
+            {
+               int i, j, k; // 3D lexicographic index of quad point p
+               internal::FaceIdxToVolIdx3D(p, d, fid0, fid1, side, orientation, i, j, k);
+
+               f2v(p, side, f) = i + d * (j + d * k);
+            }
+         }
+      }
+   }
+}
+
+L2NormalDerivativeFaceRestriction::L2NormalDerivativeFaceRestriction(
+   const FiniteElementSpace &fes_,
+   const ElementDofOrdering f_ordering,
+   const FaceType face_type_)
+   : fes(fes_),
+     face_type(face_type_),
+     dim(fes.GetMesh()->Dimension()),
+     nf(fes.GetNFbyType(face_type)),
+     ne(fes.GetNE())
+{
+   MFEM_VERIFY(f_ordering == ElementDofOrdering::LEXICOGRAPHIC,
+               "Non-lexicographic ordering not currently supported in "
+               "L2NormalDerivativeFaceRestriction.");
+
+   Mesh &mesh = *fes.GetMesh();
+
+   const FiniteElement &fe = *fes.GetFE(0);
+   const int d = fe.GetDofToQuad(fe.GetNodes(), DofToQuad::TENSOR).ndof;
+
+   if (dim == 2)
+   {
+      // (el0, el1, fid0, fid1)
+      face_to_elem.SetSize(nf * 4);
+      face_to_vol.SetSize(2 * nf * d);
+   }
+   else if (dim == 3)
+   {
+      // (el0, el1, fid0, fid1, or0, or1)
+      face_to_elem.SetSize(nf * 6);
+      face_to_vol.SetSize(2 * nf * d * d);
+   }
+   else
+   {
+      MFEM_ABORT("Unsupported dimension.");
+   }
+   auto f2e = Reshape(face_to_elem.HostWrite(), 2, (dim == 2) ? 2 : 3, nf);
+
+   // Populate the face_to_elem array. The elem_indicator will be used to count
+   // the number of elements that are adjacent to faces of the given type.
+   Array<int> elem_indicator(ne);
+   elem_indicator = 0;
+
+   int f_ind = 0;
+   for (int f = 0; f < fes.GetNF(); ++f)
+   {
+      Mesh::FaceInformation face = mesh.GetFaceInformation(f);
+
+      if (face.IsOfFaceType(face_type))
+      {
+         f2e(0, 0, f_ind) = face.element[0].index;
+         f2e(0, 1, f_ind) = face.element[0].local_face_id;
+         if (dim == 3)
+         {
+            f2e(0, 2, f_ind) = face.element[0].orientation;
+         }
+
+         elem_indicator[face.element[0].index] = 1;
+
+         if (face_type == FaceType::Interior)
+         {
+            const int el_idx_1 = face.element[1].index;
+            if (face.IsShared())
+            {
+               // Indicate shared face by index >= ne
+               f2e(1, 0, f_ind) = ne + el_idx_1;
+            }
+            else
+            {
+               // Face is not shared
+               f2e(1, 0, f_ind) = el_idx_1;
+               elem_indicator[el_idx_1] = 1;
+            }
+            f2e(1, 1, f_ind) = face.element[1].local_face_id;
+
+            if (dim == 3)
+            {
+               f2e(1, 2, f_ind) = face.element[1].orientation;
+            }
+         }
+         else
+         {
+            f2e(1, 0, f_ind) = -1;
+            f2e(1, 1, f_ind) = -1;
+
+            if (dim == 3)
+            {
+               f2e(1, 2, f_ind) = -1;
+            }
+         }
+
+         f_ind++;
+      }
+   }
+
+   // evaluate face to vol map
+   if (dim == 2)
+   {
+      NormalDerivativeSetupFaceIndexMap2D(nf, d, face_to_elem, face_to_vol);
+   }
+   else if (dim == 3)
+   {
+      NormalDerivativeSetupFaceIndexMap3D(nf, d, face_to_elem, face_to_vol);
+   }
+
+   // Number of elements adjacent to faces of face_type
+   ne_type = elem_indicator.Sum();
+
+   // In 2D: (el, f0,f1,f2,f3, s0,s1,s2,s3)
+   // In 3D: (el, f0,f1,f2,f3,f4,f5, s0,s1,s2,s3,s4,s5)
+   const int elem_data_sz = (dim == 2) ? 9 : 13;
+
+   elem_to_face.SetSize(elem_data_sz * ne_type);
+   elem_to_face = -1;
+
+   auto e2f = Reshape(elem_to_face.HostWrite(), elem_data_sz, ne_type);
+   elem_indicator.PartialSum();
+
+   const int nsides = (face_type == FaceType::Interior) ? 2 : 1;
+   const int side_begin = (dim == 2) ? 5 : 7;
+   for (int f = 0; f < nf; ++f)
+   {
+      for (int side = 0; side < nsides; ++side)
+      {
+         const int el = f2e(side, 0, f);
+         // Skip shared faces
+         if (el < ne)
+         {
+            const int face_id = f2e(side, 1, f);
+
+            const int e = elem_indicator[el] - 1;
+            e2f(0, e) = el;
+            e2f(1 + face_id, e) = f;
+            e2f(side_begin + face_id, e) = side;
+         }
+      }
+   }
+}
+
+void L2NormalDerivativeFaceRestriction::Mult(const Vector &x, Vector &y) const
+{
+   if (nf == 0) { return; }
+   switch (dim)
+   {
+      case 2:
+      {
+         const int d1d = fes.GetElementOrder(0) + 1;
+         switch (d1d)
+         {
+            case 1: Mult2D<1>(x, y); break;
+            case 2: Mult2D<2>(x, y); break;
+            case 3: Mult2D<3>(x, y); break;
+            case 4: Mult2D<4>(x, y); break;
+            case 5: Mult2D<5>(x, y); break;
+            case 6: Mult2D<6>(x, y); break;
+            case 7: Mult2D<7>(x, y); break;
+            case 8: Mult2D<8>(x, y); break;
+            default: Mult2D(x, y); break;
+         }
+      }
+      break;
+      case 3:
+      {
+         const int d1d = fes.GetElementOrder(0) + 1;
+         switch (d1d)
+         {
+            case 1: Mult3D<1>(x, y); break;
+            case 2: Mult3D<2>(x, y); break;
+            case 3: Mult3D<3>(x, y); break;
+            case 4: Mult3D<4>(x, y); break;
+            case 5: Mult3D<5>(x, y); break;
+            case 6: Mult3D<6>(x, y); break;
+            case 7: Mult3D<7>(x, y); break;
+            case 8: Mult3D<8>(x, y); break;
+            default: Mult3D(x, y); break; // fallback
+         }
+         break;
+      }
+      default: MFEM_ABORT("Dimension not supported."); break;
+   }
+}
+
+void L2NormalDerivativeFaceRestriction::AddMultTranspose(
+   const Vector &x, Vector &y, const real_t a) const
+{
+   if (nf == 0) { return; }
+   switch (dim)
+   {
+      case 2:
+      {
+         const int d1d = fes.GetElementOrder(0) + 1;
+         switch (d1d)
+         {
+            case 1: AddMultTranspose2D<1>(x, y, a); break;
+            case 2: AddMultTranspose2D<2>(x, y, a); break;
+            case 3: AddMultTranspose2D<3>(x, y, a); break;
+            case 4: AddMultTranspose2D<4>(x, y, a); break;
+            case 5: AddMultTranspose2D<5>(x, y, a); break;
+            case 6: AddMultTranspose2D<6>(x, y, a); break;
+            case 7: AddMultTranspose2D<7>(x, y, a); break;
+            case 8: AddMultTranspose2D<8>(x, y, a); break;
+            default: AddMultTranspose2D(x, y, a); break;
+         }
+      }
+      break;
+      case 3:
+      {
+         const int d1d = fes.GetElementOrder(0) + 1;
+         switch (d1d)
+         {
+            case 1: AddMultTranspose3D<1>(x, y, a); break;
+            case 2: AddMultTranspose3D<2>(x, y, a); break;
+            case 3: AddMultTranspose3D<3>(x, y, a); break;
+            case 4: AddMultTranspose3D<4>(x, y, a); break;
+            case 5: AddMultTranspose3D<5>(x, y, a); break;
+            case 6: AddMultTranspose3D<6>(x, y, a); break;
+            case 7: AddMultTranspose3D<7>(x, y, a); break;
+            case 8: AddMultTranspose3D<8>(x, y, a); break;
+            default: AddMultTranspose3D(x, y, a); break; // fallback
+         }
+         break;
+      }
+      default: MFEM_ABORT("Not yet implemented"); break;
+   }
+}
+
+template <int T_D1D>
+void L2NormalDerivativeFaceRestriction::Mult2D(const Vector &x, Vector &y) const
+{
+   int ne_shared = 0;
+   const real_t *face_nbr_data = nullptr;
+
+#ifdef MFEM_USE_MPI
+   std::unique_ptr<ParGridFunction> x_gf;
+   if (const auto *pfes = dynamic_cast<const ParFiniteElementSpace*>(&fes))
+   {
+      if (face_type == FaceType::Interior)
+      {
+         x_gf.reset(new ParGridFunction(const_cast<ParFiniteElementSpace*>(pfes),
+                                        const_cast<Vector&>(x), 0));
+         x_gf->ExchangeFaceNbrData();
+         face_nbr_data = x_gf->FaceNbrData().Read();
+         ne_shared = pfes->GetParMesh()->GetNFaceNeighborElements();
+      }
+   }
+#endif
+
+   const int vd = fes.GetVDim();
+   const bool t = fes.GetOrdering() == Ordering::byVDIM;
+   const int num_elem = ne;
+
+   const FiniteElement &fe = *fes.GetFE(0);
+   const DofToQuad &maps = fe.GetDofToQuad(fe.GetNodes(), DofToQuad::TENSOR);
+
+   const int q = maps.nqpt;
+   const int d = maps.ndof;
+
+   MFEM_VERIFY(q == d, "");
+   MFEM_VERIFY(T_D1D == d || T_D1D == 0, "");
+
+   // derivative of 1D basis function
+   const auto G_ = Reshape(maps.G.Read(), q, d);
+   // (el0, el1, fid0, fid1)
+   const auto f2e = Reshape(face_to_elem.Read(), 2, 2, nf);
+
+   const auto f2v = Reshape(face_to_vol.Read(), q, 2, nf);
+
+   // if byvdim, d_x has shape (vdim, nddof, nddof, ne)
+   // otherwise, d_x has shape (nddof, nddof, ne, vdim)
+   const auto d_x = Reshape(x.Read(), t?vd:d, d, t?d:ne, t?ne:vd);
+   const auto d_x_shared = Reshape(face_nbr_data,
+                                   t?vd:d, d, t?d:ne_shared, t?ne_shared:vd);
+   auto d_y = Reshape(y.Write(), q, vd, 2, nf);
+
+   mfem::forall_2D(nf, 2, q, [=] MFEM_HOST_DEVICE (int f) -> void
+   {
+      constexpr int MD = (T_D1D) ? T_D1D : DofQuadLimits::MAX_D1D;
+
+      MFEM_SHARED real_t G_s[MD*MD];
+      DeviceMatrix G(G_s, q, d);
+
+      MFEM_SHARED int E[2];
+      MFEM_SHARED int FID[2];
+      MFEM_SHARED int F2V[2][MD];
+
+      if (MFEM_THREAD_ID(x) == 0)
+      {
+         MFEM_FOREACH_THREAD(j, y, d)
+         {
+            for (int i = 0; i < q; ++i)
+            {
+               G(i, j) = G_(i, j);
+            }
+         }
+      }
+
+      MFEM_FOREACH_THREAD(side, x, 2)
+      {
+         if (MFEM_THREAD_ID(y) == 0)
+         {
+            E[side] = f2e(side, 0, f);
+            FID[side] = f2e(side, 1, f);
+         }
+
+         MFEM_FOREACH_THREAD(j, y, d)
+         {
+            F2V[side][j] = f2v(j, side, f);
+         }
+      }
+      MFEM_SYNC_THREAD;
+
+      MFEM_FOREACH_THREAD(side, x, 2)
+      {
+         const int el = E[side];
+         const bool shared = (el >= num_elem);
+         const auto &d_x_e = shared ? d_x_shared : d_x;
+         const int el_idx = shared ? el - num_elem : el;
+
+         const int face_id = FID[side];
+
+         MFEM_FOREACH_THREAD(p, y, q)
+         {
+            if (el < 0)
+            {
+               for (int c = 0; c < vd; ++c)
+               {
+                  d_y(p, c, side, f) = 0.0;
+               }
+            }
+            else
+            {
+               const int ij = F2V[side][p];
+               const int i = ij % q;
+               const int j = ij / q;
+
+               for (int c=0; c < vd; ++c)
+               {
+                  real_t grad_n = 0;
+                  for (int kk=0; kk < d; ++kk)
+                  {
+                     const int k = (face_id == 0 || face_id == 2) ? i : kk;
+                     const int l = (face_id == 0 || face_id == 2) ? kk : j;
+                     const real_t g = (face_id == 0 || face_id == 2) ? G(j,l) : G(i,k);
+                     grad_n += g * d_x_e(t?c:k, t?k:l, t?l:el_idx, t?el_idx:c);
+                  }
+                  d_y(p, c, side, f) = grad_n;
+               }
+            }
+         }
+      }
+   });
+}
+
+template <int T_D1D>
+void L2NormalDerivativeFaceRestriction::Mult3D(const Vector &x, Vector &y) const
+{
+   int ne_shared = 0;
+   const real_t *face_nbr_data = nullptr;
+
+#ifdef MFEM_USE_MPI
+   std::unique_ptr<ParGridFunction> x_gf;
+   if (const auto *pfes = dynamic_cast<const ParFiniteElementSpace*>(&fes))
+   {
+      if (face_type == FaceType::Interior)
+      {
+         x_gf.reset(new ParGridFunction(const_cast<ParFiniteElementSpace*>(pfes),
+                                        const_cast<Vector&>(x), 0));
+         x_gf->ExchangeFaceNbrData();
+         face_nbr_data = x_gf->FaceNbrData().Read();
+         ne_shared = pfes->GetParMesh()->GetNFaceNeighborElements();
+      }
+   }
+#endif
+
+   const int vd = fes.GetVDim();
+   const bool t = fes.GetOrdering() == Ordering::byVDIM;
+   const int num_elem = ne;
+
+   const FiniteElement &fe = *fes.GetFE(0);
+   const DofToQuad &maps = fe.GetDofToQuad(fe.GetNodes(), DofToQuad::TENSOR);
+
+   const int q = maps.nqpt;
+   const int d = maps.ndof;
+   const int q2d = q * q;
+
+   MFEM_VERIFY(q == d, "");
+   MFEM_VERIFY(T_D1D == d || T_D1D == 0, "");
+
+   const auto G_ = Reshape(maps.G.Read(), q, d);
+   // (el0, el1, fid0, fid1, or0, or1)
+   const auto f2e = Reshape(face_to_elem.Read(), 2, 3, nf);
+   const auto f2v = Reshape(face_to_vol.Read(), q2d, 2, nf);
+
+   // t ? (vdim, d, d, d, ne) : (d, d, d, ne, vdim)
+   const auto d_x = Reshape(x.Read(), t?vd:d, d, d, t?d:ne, t?ne:vd);
+   const auto d_x_shared = Reshape(face_nbr_data,
+                                   t?vd:d, d, d, t?d:ne_shared, t?ne_shared:vd);
+   auto d_y = Reshape(y.Write(), q2d, vd, 2, nf);
+
+   mfem::forall_2D(nf, q2d, 2, [=] MFEM_HOST_DEVICE (int f) -> void
+   {
+      static constexpr int MD = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
+
+      MFEM_SHARED real_t G_s[MD*MD];
+      DeviceMatrix G(G_s, d, q);
+
+      MFEM_SHARED int E[2];
+      MFEM_SHARED int FID[2];
+      MFEM_SHARED int F2V[2][MD*MD];
+
+      // Load G matrix into shared memory
+      if (MFEM_THREAD_ID(y) == 0)
+      {
+         MFEM_FOREACH_THREAD(j, x, d*q)
+         {
+            const int p = j % q;
+            const int k = j / q;
+            G(k, p) = G_(p, k);
+         }
+      }
+
+      MFEM_FOREACH_THREAD(side, y, 2)
+      {
+         if (MFEM_THREAD_ID(x) == 0)
+         {
+            E[side] = f2e(side, 0, f);
+            FID[side] = f2e(side, 1, f);
+         }
+         MFEM_FOREACH_THREAD(j, x, q2d)
+         {
+            F2V[side][j] = f2v(j, side, f);
+         }
+      }
+      MFEM_SYNC_THREAD;
+
+      MFEM_FOREACH_THREAD(side, y, 2)
+      {
+         const int el = E[side];
+         const bool shared = (el >= num_elem);
+         const auto &d_x_e = shared ? d_x_shared : d_x;
+         const int el_idx = shared ? el - num_elem : el;
+
+         const int face_id = FID[side];
+
+         // Is this face parallel to the x-y plane in reference coordinates?
+         const bool xy_plane = (face_id == 0 || face_id == 5);
+         const bool xz_plane = (face_id == 1 || face_id == 3);
+         const bool yz_plane = (face_id == 2 || face_id == 4);
+
+         MFEM_FOREACH_THREAD(p, x, q2d)
+         {
+            if (el_idx < 0)
+            {
+               for (int c = 0; c < vd; ++c)
+               {
+                  d_y(p, c, side, f) = 0.0;
+               }
+            }
+            else
+            {
+               const int ijk = F2V[side][p];
+               const int k = ijk / q2d;
+               const int i = ijk % q;
+               const int j = (ijk - q2d*k) / q;
+
+               // the fixed 1D index of the normal component of the face
+               // quadrature point
+               const int g_row = yz_plane ? i : xz_plane ? j : k;
+
+               for (int c = 0; c < vd; ++c)
+               {
+                  real_t grad_n = 0.0;
+
+                  for (int kk = 0; kk < d; ++kk)
+                  {
+                     // (l, m, n) 3D lexicographic index of interior points used
+                     // in evaluating normal derivatives
+                     const int l = yz_plane ? kk : i;
+                     const int m = xz_plane ? kk : j;
+                     const int n = xy_plane ? kk : k;
+
+                     const real_t g = G(kk, g_row);
+
+                     grad_n += g * d_x_e(t?c:l, t?l:m, t?m:n, t?n:el_idx, t?el_idx:c);
+                  }
+                  d_y(p, c, side, f) = grad_n;
+               }
+            }
+         }
+      }
+   });
+}
+
+template <int T_D1D>
+void L2NormalDerivativeFaceRestriction::AddMultTranspose2D(
+   const Vector &y, Vector &x, const real_t a) const
+{
+   const int vd = fes.GetVDim();
+   const bool t = fes.GetOrdering() == Ordering::byVDIM;
+
+   const FiniteElement &fe = *fes.GetFE(0);
+   const DofToQuad &maps = fe.GetDofToQuad(fe.GetNodes(), DofToQuad::TENSOR);
+
+   const int q = maps.nqpt;
+   const int d = maps.ndof;
+
+   // derivative of 1D basis function
+   auto G_ = Reshape(maps.G.Read(), q, d);
+
+   // entries of e2f: (el,f0,f1,f2,f3,s0,s1,s2,s3)
+   auto e2f = Reshape(elem_to_face.Read(), 9, ne_type);
+
+   auto f2v = Reshape(face_to_vol.Read(), d, 2, nf);
+
+   // if byvdim, d_x has shape (vdim, nddof, nddof, ne)
+   // otherwise, d_x has shape (nddof, nddof, ne, vdim)
+   auto d_x = Reshape(x.ReadWrite(), t?vd:d, d, t?d:ne, t?ne:vd);
+   auto d_y = Reshape(y.Read(), q, vd, 2, nf);
+
+   mfem::forall_2D(ne_type, d, d, [=] MFEM_HOST_DEVICE (int e)
+   {
+      constexpr int MD = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
+
+      MFEM_SHARED real_t y_s[MD];
+      MFEM_SHARED int pp[MD];
+      MFEM_SHARED int jj;
+      if (MFEM_THREAD_ID(x) == 0 && MFEM_THREAD_ID(y) == 0) { jj = 0; }
+
+      MFEM_SHARED real_t BG[MD*MD];
+      DeviceMatrix G(BG, q, d);
+
+      MFEM_SHARED real_t x_s[MD*MD];
+      DeviceMatrix xx(x_s, d, d);
+
+      MFEM_SHARED int el; // global element index
+      MFEM_SHARED int faces[4];
+      MFEM_SHARED int sides[4];
+
+      MFEM_FOREACH_THREAD(i,x,d)
+      {
+         MFEM_FOREACH_THREAD(p,y,q)
+         {
+            G(p,i) = a * G_(p,i);
+            xx(p,i) = 0.0;
+         }
+      }
+
+      if (MFEM_THREAD_ID(y) == 0)
+      {
+         if (MFEM_THREAD_ID(x) == 0)
+         {
+            el = e2f(0, e);
+         }
+
+         MFEM_FOREACH_THREAD(i, x, 4)
+         {
+            faces[i] = e2f(1 + i, e);
+            sides[i] = e2f(5 + i, e);
+         }
+      }
+      MFEM_SYNC_THREAD;
+
+      for (int face_id=0; face_id < 4; ++face_id)
+      {
+         const int f = faces[face_id];
+
+         if (f < 0) { continue; }
+
+         const int side = sides[face_id];
+
+         if (MFEM_THREAD_ID(y) == 0)
+         {
+            MFEM_FOREACH_THREAD(p,x,d)
+            {
+               y_s[p] = d_y(p, 0, side, f);
+
+               const int ij = f2v(p, side, f);
+               const int i = ij % q;
+               const int j = ij / q;
+
+               pp[(face_id == 0 || face_id == 2) ? i : j] = p;
+               if (MFEM_THREAD_ID(x) == 0)
+               {
+                  jj = (face_id == 0 || face_id == 2) ? j : i;
+               }
+            }
+         }
+         MFEM_SYNC_THREAD;
+
+         MFEM_FOREACH_THREAD(k,x,d)
+         {
+            MFEM_FOREACH_THREAD(l,y,d)
+            {
+               const int p = (face_id == 0 || face_id == 2) ? pp[k] : pp[l];
+               const int kk = (face_id == 0 || face_id == 2) ? l : k;
+               const real_t g = G(jj, kk);
+               xx(k,l) += g * y_s[p];
+            }
+         }
+      }
+      MFEM_SYNC_THREAD;
+
+      MFEM_FOREACH_THREAD(k,x,d)
+      {
+         MFEM_FOREACH_THREAD(l,y,d)
+         {
+            const int c = 0;
+            d_x(t?c:k, t?k:l, t?l:el, t?el:c) += xx(k,l);
+         }
+      }
+   });
+}
+
+template <int T_D1D>
+void L2NormalDerivativeFaceRestriction::AddMultTranspose3D(
+   const Vector &y, Vector &x, const real_t a) const
+{
+   const int vd = fes.GetVDim();
+   const bool t = fes.GetOrdering() == Ordering::byVDIM;
+
+   MFEM_VERIFY(vd == 1, "vdim > 1 not supported.");
+
+   const FiniteElement &fe = *fes.GetFE(0);
+   const DofToQuad &maps = fe.GetDofToQuad(fe.GetNodes(), DofToQuad::TENSOR);
+
+   const int q = maps.nqpt;
+   const int d = maps.ndof;
+   const int q2d = q * q;
+
+   MFEM_VERIFY(q == d, "");
+   MFEM_VERIFY(T_D1D == d || T_D1D == 0, "");
+
+   auto G_ = Reshape(maps.G.Read(), q, d);
+
+   // (el, f0,f1,f2,f3,f4,f5, s0,s1,s2,s3,s4,s5)
+   auto e2f = Reshape(elem_to_face.Read(), 13, ne_type);
+
+   auto f2v = Reshape(face_to_vol.Read(), q2d, 2, nf);
+
+   auto d_x = Reshape(x.ReadWrite(), t?vd:d, d, d, t?d:ne, t?ne:vd);
+   const auto d_y = Reshape(y.Read(), q2d, vd, 2, nf);
+
+   mfem::forall_2D(ne_type, q, q, [=] MFEM_HOST_DEVICE (int e) -> void
+   {
+      static constexpr int MD = T_D1D ? T_D1D : DofQuadLimits::MAX_D1D;
+
+      MFEM_SHARED int pp[MD][MD];
+      MFEM_SHARED real_t y_s[MD*MD];
+      MFEM_SHARED int jj;
+      if (MFEM_THREAD_ID(x) == 0 && MFEM_THREAD_ID(y) == 0) { jj = 0; }
+
+      MFEM_SHARED real_t xx_s[MD*MD*MD];
+      auto xx = Reshape(xx_s, d, d, d);
+
+      MFEM_SHARED real_t G_s[MD*MD];
+      DeviceMatrix G(G_s, q, d);
+
+      MFEM_SHARED int el;
+      MFEM_SHARED int faces[6];
+      MFEM_SHARED int sides[6];
+
+      // Load G into shared memory
+      MFEM_FOREACH_THREAD(j, x, d)
+      {
+         MFEM_FOREACH_THREAD(i, y, q)
+         {
+            G(i, j) = a * G_(i, j);
+            G(i, j) = a * G_(i, j);
+            G(i, j) = a * G_(i, j);
+         }
+      }
+
+      if (MFEM_THREAD_ID(y) == 0)
+      {
+         if (MFEM_THREAD_ID(x) == 0)
+         {
+            el = e2f(0, e); // global element index
+         }
+
+         MFEM_FOREACH_THREAD(i, x, 6)
+         {
+            faces[i] = e2f(1 + i, e);
+            sides[i] = e2f(7 + i, e);
+         }
+      }
+
+      MFEM_FOREACH_THREAD(k, x, d)
+      {
+         MFEM_FOREACH_THREAD(j, y, d)
+         {
+            for (int i = 0; i < d; ++i)
+            {
+               xx(i, j, k) = 0.0;
+            }
+         }
+      }
+      MFEM_SYNC_THREAD;
+
+      for (int face_id = 0; face_id < 6; ++face_id)
+      {
+         const int f = faces[face_id];
+
+         if (f < 0)
+         {
+            continue;
+         }
+
+         const int side = sides[face_id];
+
+         // is this face parallel to the x-y plane in reference coordinates?
+         const bool xy_plane = (face_id == 0 || face_id == 5);
+         const bool xz_plane = (face_id == 1 || face_id == 3);
+
+         MFEM_FOREACH_THREAD(p1, x, q)
+         {
+            MFEM_FOREACH_THREAD(p2, y, q)
+            {
+               const int p = p1 + q * p2;
+               y_s[p] = d_y(p, 0, side, f);
+
+               const int ijk = f2v(p, side, f);
+               const int k = ijk / q2d;
+               const int i = ijk % q;
+               const int j = (ijk - q2d*k) / q;
+
+               pp[(xy_plane || xz_plane) ? i : j][(xy_plane) ? j : k] = p;
+               if (MFEM_THREAD_ID(x) == 0 && MFEM_THREAD_ID(y) == 0)
+               {
+                  jj = (xy_plane) ? k : (xz_plane) ? j : i;
+               }
+            }
+         }
+         MFEM_SYNC_THREAD;
+
+         MFEM_FOREACH_THREAD(n, x, d)
+         {
+            MFEM_FOREACH_THREAD(m, y, d)
+            {
+               for (int l = 0; l < d; ++l)
+               {
+                  const int p = (xy_plane) ? pp[l][m] : (xz_plane) ? pp[l][n] : pp[m][n];
+                  const int kk = (xy_plane) ? n : (xz_plane) ? m : l;
+                  const real_t g = G(jj, kk);
+                  xx(l, m, n) += g * y_s[p];
+               }
+            }
+         }
+      }
+      MFEM_SYNC_THREAD;
+
+      // map back to global array
+      MFEM_FOREACH_THREAD(n, x, d)
+      {
+         MFEM_FOREACH_THREAD(m, y, d)
+         {
+            for (int l = 0; l < d; ++l)
+            {
+               const int c = 0;
+               d_x(t?c:l, t?l:m, t?m:n, t?n:el, t?el:c) += xx(l, m, n);
+            }
+         }
+      }
+   });
+}
+
+} // namespace mfem

fem/normal_deriv_restriction.hpp

@@ -0,0 +1,85 @@
+// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_NORMAL_DERIV_RESTRICTION
+#define MFEM_NORMAL_DERIV_RESTRICTION
+
+#include "../mesh/mesh.hpp"
+
+namespace mfem
+{
+
+class FiniteElementSpace;
+enum class ElementDofOrdering;
+
+/// @brief Class to compute face normal derivatives (in reference coordinate) of
+/// an L2 grid function (used internally by L2FaceRestriction).
+class L2NormalDerivativeFaceRestriction
+{
+protected:
+   const FiniteElementSpace &fes; ///< The L2 finite element space.
+   const FaceType face_type; ///< Face type: either boundary or interior.
+   const int dim; ///< Dimension of the mesh.
+   const int nf; ///< Number of faces of the given @a face_type.
+   const int ne; ///< Number of elements.
+   int ne_type; ///< Number of elements with faces of type face type
+
+   Array<int> face_to_elem; ///< Face-wise information array.
+   Array<int> elem_to_face; ///< Element-wise information array.
+   Array<int> face_to_vol; ///< maps face index to volume index
+
+public:
+   /// @brief Constructor.
+   /// @param[in] fes_ The associated FiniteElementSpace (should be L2/DG).
+   /// @param[in] f_ordering Request a specific face dof ordering. Currently
+   /// only ElementDofOrdering::LEXICOGRAPHIC is supported.
+   /// @param[in] face_type_ Type of faces to compute restriction (interior or boundary).
+   L2NormalDerivativeFaceRestriction(const FiniteElementSpace &fes_,
+                                     const ElementDofOrdering f_ordering,
+                                     const FaceType face_type_);
+
+   /// @brief Computes the normal derivatives on the @a face_type faces of the mesh.
+   /// @param[in] x The L-vector degrees of freedom.
+   /// @param[out] y The face E(like)-vector degrees of freedom of the format
+   /// (face_dofs x vdim x 2 x nf) where nf is the number of faces of type @a
+   /// face_type. The face_dofs are ordered according to @a f_ordering specified
+   /// in the constructor.
+   void Mult(const Vector &x, Vector &y) const;
+
+   /// @brief Computes the transpose of the action of Mult(), accumulating into
+   ///        @a y with coefficient @a a.
+   /// @param x Face E-vector layout (face_dofs x vdim x 2 x nf).
+   /// @param y L-vector layout.
+   /// @param a Optional coefficient (y = y + a*R^t*x)
+   void AddMultTranspose(const Vector &x, Vector &y,
+                         const real_t a = 1.0) const;
+
+   /// @name Internal compute kernels. Public because of nvcc restriction.
+   ///@{
+
+   template <int T_D1D = 0>
+   void Mult2D(const Vector &x, Vector &y) const;
+
+   template <int T_D1D = 0>
+   void AddMultTranspose2D(const Vector &x, Vector &y, const real_t a) const;
+
+   template <int T_D1D = 0>
+   void Mult3D(const Vector &x, Vector &y) const;
+
+   template <int T_D1D = 0>
+   void AddMultTranspose3D(const Vector &x, Vector &y, const real_t a) const;
+
+   /// @}
+};
+
+}
+
+#endif // MFEM_RESTRICTION

fem/prestriction.cpp

@@ -18,6 +18,7 @@
 #include "pgridfunc.hpp"
 #include "pfespace.hpp"
 #include "fespace.hpp"
+#include "fe/face_map_utils.hpp"
 #include "../general/forall.hpp"
 
 namespace mfem
@@ -277,21 +278,22 @@ void ParNCH1FaceRestriction::ComputeGatherIndices(
    gather_offsets[0] = 0;
 }
 
-ParL2FaceRestriction::ParL2FaceRestriction(const ParFiniteElementSpace &fes,
+ParL2FaceRestriction::ParL2FaceRestriction(const ParFiniteElementSpace &pfes_,
                                            ElementDofOrdering f_ordering,
                                            FaceType type,
                                            L2FaceValues m,
                                            bool build)
-   : L2FaceRestriction(fes, f_ordering, type, m, false)
+   : L2FaceRestriction(pfes_, f_ordering, type, m, false),
+     pfes(pfes_)
 {
    if (!build) { return; }
    if (nf==0) { return; }
 
-   CheckFESpace(f_ordering);
+   CheckFESpace();
 
-   ComputeScatterIndicesAndOffsets(f_ordering, type);
+   ComputeScatterIndicesAndOffsets();
 
-   ComputeGatherIndices(f_ordering, type);
+   ComputeGatherIndices();
 }
 
 ParL2FaceRestriction::ParL2FaceRestriction(const ParFiniteElementSpace &fes,
@@ -307,8 +309,6 @@ void ParL2FaceRestriction::DoubleValuedConformingMult(
    MFEM_ASSERT(
       m == L2FaceValues::DoubleValued,
       "This method should be called when m == L2FaceValues::DoubleValued.");
-   const ParFiniteElementSpace &pfes =
-      static_cast<const ParFiniteElementSpace&>(this->fes);
    ParGridFunction x_gf;
    x_gf.MakeRef(const_cast<ParFiniteElementSpace*>(&pfes),
                 const_cast<Vector&>(x), 0);
@@ -567,13 +567,9 @@ void ParL2FaceRestriction::FillJAndData(const Vector &ea_data,
    });
 }
 
-void ParL2FaceRestriction::ComputeScatterIndicesAndOffsets(
-   const ElementDofOrdering f_ordering,
-   const FaceType type)
+void ParL2FaceRestriction::ComputeScatterIndicesAndOffsets()
 {
    Mesh &mesh = *fes.GetMesh();
-   const ParFiniteElementSpace &pfes =
-      static_cast<const ParFiniteElementSpace&>(this->fes);
 
    // Initialization of the offsets
    for (int i = 0; i <= ndofs; ++i)
@@ -622,9 +618,7 @@ void ParL2FaceRestriction::ComputeScatterIndicesAndOffsets(
 }
 
 
-void ParL2FaceRestriction::ComputeGatherIndices(
-   const ElementDofOrdering f_ordering,
-   const FaceType type)
+void ParL2FaceRestriction::ComputeGatherIndices()
 {
    Mesh &mesh = *fes.GetMesh();
 
@@ -666,11 +660,11 @@ ParNCL2FaceRestriction::ParNCL2FaceRestriction(const ParFiniteElementSpace &fes,
    if (nf==0) { return; }
    x_interp.UseDevice(true);
 
-   CheckFESpace(f_ordering);
+   CheckFESpace();
 
-   ComputeScatterIndicesAndOffsets(f_ordering, type);
+   ComputeScatterIndicesAndOffsets();
 
-   ComputeGatherIndices(f_ordering, type);
+   ComputeGatherIndices();
 }
 
 void ParNCL2FaceRestriction::SingleValuedNonconformingMult(
@@ -979,9 +973,7 @@ void ParNCL2FaceRestriction::FillJAndData(const Vector &ea_data,
    MFEM_ABORT("Not yet implemented.");
 }
 
-void ParNCL2FaceRestriction::ComputeScatterIndicesAndOffsets(
-   const ElementDofOrdering f_ordering,
-   const FaceType type)
+void ParNCL2FaceRestriction::ComputeScatterIndicesAndOffsets()
 {
    Mesh &mesh = *fes.GetMesh();
 
@@ -1064,9 +1056,7 @@ void ParNCL2FaceRestriction::ComputeScatterIndicesAndOffsets(
    interpolations.InitializeNCInterpConfig();
 }
 
-void ParNCL2FaceRestriction::ComputeGatherIndices(
-   const ElementDofOrdering f_ordering,
-   const FaceType type)
+void ParNCL2FaceRestriction::ComputeGatherIndices()
 {
    Mesh &mesh = *fes.GetMesh();
 

fem/prestriction.hpp

@@ -139,9 +139,11 @@ public: // For nvcc
 class ParL2FaceRestriction : virtual public L2FaceRestriction
 {
 protected:
+   const ParFiniteElementSpace &pfes;
+
    /** @brief Constructs an ParL2FaceRestriction.
 
-       @param[in] fes        The ParFiniteElementSpace on which this operates
+       @param[in] pfes_      The ParFiniteElementSpace on which this operates
        @param[in] f_ordering Request a specific face dof ordering
        @param[in] type       Request internal or boundary faces dofs
        @param[in] m          Request the face dofs for elem1, or both elem1 and
@@ -149,7 +151,7 @@ protected:
        @param[in] build      Request the ParL2FaceRestriction to compute the
                              scatter/gather indices. False should only be used
                              when inheriting from ParL2FaceRestriction. */
-   ParL2FaceRestriction(const ParFiniteElementSpace& fes,
+   ParL2FaceRestriction(const ParFiniteElementSpace &pfes_,
                         ElementDofOrdering f_ordering,
                         FaceType type,
                         L2FaceValues m,
@@ -229,22 +231,14 @@ public:
 private:
    /** @brief Compute the scatter indices: L-vector to E-vector, and the offsets
        for the gathering: E-vector to L-vector.
-
-       @param[in] f_ordering Request a specific face dof ordering.
-       @param[in] type       Request internal or boundary faces dofs.
    */
-   void ComputeScatterIndicesAndOffsets(const ElementDofOrdering f_ordering,
-                                        const FaceType type);
+   void ComputeScatterIndicesAndOffsets();
 
    /** @brief Compute the gather indices: E-vector to L-vector.
 
        Note: Requires the gather offsets to be computed.
-
-       @param[in] f_ordering Request a specific face dof ordering.
-       @param[in] type       Request internal or boundary faces dofs.
    */
-   void ComputeGatherIndices(const ElementDofOrdering f_ordering,
-                             const FaceType type);
+   void ComputeGatherIndices();
 
 public:
    /** @brief Scatter the degrees of freedom, i.e. goes from L-Vector to
@@ -381,22 +375,14 @@ private:
    /** @brief Compute the scatter indices: L-vector to E-vector, the offsets
        for the gathering: E-vector to L-vector, and the interpolators from
        coarse to fine face for master non-comforming faces.
-
-       @param[in] f_ordering Request a specific face dof ordering.
-       @param[in] type       Request internal or boundary faces dofs.
    */
-   void ComputeScatterIndicesAndOffsets(const ElementDofOrdering f_ordering,
-                                        const FaceType type);
+   void ComputeScatterIndicesAndOffsets();
 
    /** @brief Compute the gather indices: E-vector to L-vector.
 
        Note: Requires the gather offsets to be computed.
-
-       @param[in] f_ordering Request a specific face dof ordering.
-       @param[in] type       Request internal or boundary faces dofs.
    */
-   void ComputeGatherIndices(const ElementDofOrdering f_ordering,
-                             const FaceType type);
+   void ComputeGatherIndices();
 
 public:
    /** @brief Scatter the degrees of freedom, i.e. goes from L-Vector to

fem/qspace.cpp

@@ -139,9 +139,9 @@ QuadratureSpace::QuadratureSpace(Mesh *mesh_, std::istream &in)
 }
 
 QuadratureSpace::QuadratureSpace(Mesh &mesh_, const IntegrationRule &ir)
-   : QuadratureSpaceBase(mesh_, mesh_.GetElementGeometry(0), ir)
+   : QuadratureSpaceBase(mesh_, mesh_.GetTypicalElementGeometry(), ir)
 {
-   MFEM_VERIFY(mesh.GetNumGeometries(mesh.Dimension()) == 1,
+   MFEM_VERIFY(mesh.GetNumGeometries(mesh.Dimension()) <= 1,
                "Constructor not valid for mixed meshes");
    ConstructOffsets();
 }

fem/qspace.hpp

@@ -215,6 +215,10 @@ public:
    /// quadrature point, oriented relative to "element 1".
    int GetPermutedIndex(int idx, int iq) const override;
 
+   /// @brief Get the face index (in the standard Mesh numbering) associated
+   /// with face @a idx in the FaceQuadratureSpace.
+   int GetMeshFaceIndex(int idx) const { return face_indices[idx]; }
+
    /// @brief Returns the index associated with the face described by @a T.
    ///
    /// The index may differ from the mesh face or boundary element index

fem/restriction.cpp

@@ -10,16 +10,15 @@
 // CONTRIBUTING.md for details.
 
 #include "restriction.hpp"
+#include "normal_deriv_restriction.hpp"
 #include "gridfunc.hpp"
 #include "fespace.hpp"
-#include "../general/forall.hpp"
-#include <climits>
-
-#ifdef MFEM_USE_MPI
-
 #include "pfespace.hpp"
+#include "qspace.hpp"
+#include "fe/face_map_utils.hpp"
+#include "../general/forall.hpp"
 
-#endif
+#include <climits>
 
 namespace mfem
 {
@@ -907,112 +906,6 @@ void ConformingFaceRestriction::SetFaceDofsGatherIndices(
    }
 }
 
-static int ToLexOrdering2D(const int face_id, const int size1d, const int i)
-{
-   if (face_id==2 || face_id==3)
-   {
-      return size1d-1-i;
-   }
-   else
-   {
-      return i;
-   }
-}
-
-static int PermuteFace2D(const int face_id1, const int face_id2,
-                         const int orientation,
-                         const int size1d, const int index)
-{
-   int new_index;
-   // Convert from lex ordering
-   if (face_id1==2 || face_id1==3)
-   {
-      new_index = size1d-1-index;
-   }
-   else
-   {
-      new_index = index;
-   }
-   // Permute based on face orientations
-   if (orientation==1)
-   {
-      new_index = size1d-1-new_index;
-   }
-   return ToLexOrdering2D(face_id2, size1d, new_index);
-}
-
-static int ToLexOrdering3D(const int face_id, const int size1d, const int i,
-                           const int j)
-{
-   if (face_id==2 || face_id==1 || face_id==5)
-   {
-      return i + j*size1d;
-   }
-   else if (face_id==3 || face_id==4)
-   {
-      return (size1d-1-i) + j*size1d;
-   }
-   else // face_id==0
-   {
-      return i + (size1d-1-j)*size1d;
-   }
-}
-
-static int PermuteFace3D(const int face_id1, const int face_id2,
-                         const int orientation,
-                         const int size1d, const int index)
-{
-   int i=0, j=0, new_i=0, new_j=0;
-   i = index%size1d;
-   j = index/size1d;
-   // Convert from lex ordering
-   if (face_id1==3 || face_id1==4)
-   {
-      i = size1d-1-i;
-   }
-   else if (face_id1==0)
-   {
-      j = size1d-1-j;
-   }
-   // Permute based on face orientations
-   switch (orientation)
-   {
-      case 0:
-         new_i = i;
-         new_j = j;
-         break;
-      case 1:
-         new_i = j;
-         new_j = i;
-         break;
-      case 2:
-         new_i = j;
-         new_j = (size1d-1-i);
-         break;
-      case 3:
-         new_i = (size1d-1-i);
-         new_j = j;
-         break;
-      case 4:
-         new_i = (size1d-1-i);
-         new_j = (size1d-1-j);
-         break;
-      case 5:
-         new_i = (size1d-1-j);
-         new_j = (size1d-1-i);
-         break;
-      case 6:
-         new_i = (size1d-1-j);
-         new_j = i;
-         break;
-      case 7:
-         new_i = i;
-         new_j = (size1d-1-j);
-         break;
-   }
-   return ToLexOrdering3D(face_id2, size1d, new_i, new_j);
-}
-
 // Permute dofs or quads on a face for e2 to match with the ordering of e1
 int PermuteFaceL2(const int dim, const int face_id1,
                   const int face_id2, const int orientation,
@@ -1023,9 +916,9 @@ int PermuteFaceL2(const int dim, const int face_id1,
       case 1:
          return 0;
       case 2:
-         return PermuteFace2D(face_id1, face_id2, orientation, size1d, index);
+         return internal::PermuteFace2D(face_id1, face_id2, orientation, size1d, index);
       case 3:
-         return PermuteFace3D(face_id1, face_id2, orientation, size1d, index);
+         return internal::PermuteFace3D(face_id1, face_id2, orientation, size1d, index);
       default:
          MFEM_ABORT("Unsupported dimension.");
          return 0;
@@ -1038,6 +931,7 @@ L2FaceRestriction::L2FaceRestriction(const FiniteElementSpace &fes,
                                      const L2FaceValues m,
                                      bool build)
    : fes(fes),
+     ordering(f_ordering),
      nf(fes.GetNFbyType(type)),
      ne(fes.GetNE()),
      vdim(fes.GetVDim()),
@@ -1045,7 +939,7 @@ L2FaceRestriction::L2FaceRestriction(const FiniteElementSpace &fes,
      face_dofs(nf > 0 ?
                fes.GetTraceElement(0, fes.GetMesh()->GetFaceGeometry(0))->GetDof()
                : 0),
-     elem_dofs(fes.GetFE(0)->GetDof()),
+     elem_dofs(ne > 0 ? fes.GetFE(0)->GetDof() : 0),
      nfdofs(nf*face_dofs),
      ndofs(fes.GetNDofs()),
      type(type),
@@ -1060,11 +954,9 @@ L2FaceRestriction::L2FaceRestriction(const FiniteElementSpace &fes,
    width = fes.GetVSize();
    if (!build) { return; }
 
-   CheckFESpace(f_ordering);
-
-   ComputeScatterIndicesAndOffsets(f_ordering,type);
-
-   ComputeGatherIndices(f_ordering, type);
+   CheckFESpace();
+   ComputeScatterIndicesAndOffsets();
+   ComputeGatherIndices();
 }
 
 L2FaceRestriction::L2FaceRestriction(const FiniteElementSpace &fes,
@@ -1327,7 +1219,7 @@ void L2FaceRestriction::AddFaceMatricesToElementMatrices(const Vector &fea_data,
    }
 }
 
-void L2FaceRestriction::CheckFESpace(const ElementDofOrdering f_ordering)
+void L2FaceRestriction::CheckFESpace()
 {
 #ifdef MFEM_USE_MPI
 
@@ -1343,7 +1235,7 @@ void L2FaceRestriction::CheckFESpace(const ElementDofOrdering f_ordering)
 
 #ifdef MFEM_DEBUG
    // If fespace == L2
-   const FiniteElement *fe0 = fes.GetFE(0);
+   const FiniteElement *fe0 = fes.GetTypicalFE();
    const TensorBasisElement *tfe = dynamic_cast<const TensorBasisElement*>(fe0);
    MFEM_VERIFY(tfe != NULL &&
                (tfe->GetBasisType()==BasisType::GaussLobatto ||
@@ -1351,7 +1243,7 @@ void L2FaceRestriction::CheckFESpace(const ElementDofOrdering f_ordering)
                "Only Gauss-Lobatto and Bernstein basis are supported in "
                "L2FaceRestriction.");
    if (nf==0) { return; }
-   const bool dof_reorder = (f_ordering == ElementDofOrdering::LEXICOGRAPHIC);
+   const bool dof_reorder = (ordering == ElementDofOrdering::LEXICOGRAPHIC);
    if (!dof_reorder)
    {
       MFEM_ABORT("Non-Tensor L2FaceRestriction not yet implemented.");
@@ -1371,9 +1263,7 @@ void L2FaceRestriction::CheckFESpace(const ElementDofOrdering f_ordering)
 #endif
 }
 
-void L2FaceRestriction::ComputeScatterIndicesAndOffsets(
-   const ElementDofOrdering f_ordering,
-   const FaceType face_type)
+void L2FaceRestriction::ComputeScatterIndicesAndOffsets()
 {
    Mesh &mesh = *fes.GetMesh();
    // Initialization of the offsets
@@ -1389,16 +1279,16 @@ void L2FaceRestriction::ComputeScatterIndicesAndOffsets(
       Mesh::FaceInformation face = mesh.GetFaceInformation(f);
       MFEM_ASSERT(!face.IsShared(),
                   "Unexpected shared face in L2FaceRestriction.");
-      if ( face.IsOfFaceType(face_type) )
+      if ( face.IsOfFaceType(type) )
       {
          SetFaceDofsScatterIndices1(face,f_ind);
          if ( m==L2FaceValues::DoubleValued )
          {
-            if ( face_type==FaceType::Interior && face.IsInterior() )
+            if ( type==FaceType::Interior && face.IsInterior() )
             {
                PermuteAndSetFaceDofsScatterIndices2(face,f_ind);
             }
-            else if ( face_type==FaceType::Boundary && face.IsBoundary() )
+            else if ( type==FaceType::Boundary && face.IsBoundary() )
             {
                SetBoundaryDofsScatterIndices2(face,f_ind);
             }
@@ -1415,9 +1305,7 @@ void L2FaceRestriction::ComputeScatterIndicesAndOffsets(
    }
 }
 
-void L2FaceRestriction::ComputeGatherIndices(
-   const ElementDofOrdering f_ordering,
-   const FaceType face_type)
+void L2FaceRestriction::ComputeGatherIndices()
 {
    Mesh &mesh = *fes.GetMesh();
    // Computation of gather_indices
@@ -1427,11 +1315,11 @@ void L2FaceRestriction::ComputeGatherIndices(
       Mesh::FaceInformation face = mesh.GetFaceInformation(f);
       MFEM_ASSERT(!face.IsShared(),
                   "Unexpected shared face in L2FaceRestriction.");
-      if ( face.IsOfFaceType(face_type) )
+      if ( face.IsOfFaceType(type) )
       {
          SetFaceDofsGatherIndices1(face,f_ind);
          if ( m==L2FaceValues::DoubleValued &&
-              face_type==FaceType::Interior &&
+              type==FaceType::Interior &&
               face.IsLocal())
          {
             PermuteAndSetFaceDofsGatherIndices2(face,f_ind);
@@ -1598,6 +1486,28 @@ void L2FaceRestriction::PermuteAndSetFaceDofsGatherIndices2(
    }
 }
 
+void L2FaceRestriction::NormalDerivativeMult(const Vector &x, Vector &y) const
+{
+   EnsureNormalDerivativeRestriction();
+   normal_deriv_restr->Mult(x, y);
+}
+
+void L2FaceRestriction::NormalDerivativeAddMultTranspose(const Vector &x,
+                                                         Vector &y) const
+{
+   EnsureNormalDerivativeRestriction();
+   normal_deriv_restr->AddMultTranspose(x, y);
+}
+
+void L2FaceRestriction::EnsureNormalDerivativeRestriction() const
+{
+   if (!normal_deriv_restr)
+   {
+      normal_deriv_restr.reset(
+         new L2NormalDerivativeFaceRestriction(fes, ordering, type));
+   }
+}
+
 InterpolationManager::InterpolationManager(const FiniteElementSpace &fes,
                                            ElementDofOrdering ordering,
                                            FaceType type)
@@ -1775,11 +1685,11 @@ NCL2FaceRestriction::NCL2FaceRestriction(const FiniteElementSpace &fes,
    if (!build) { return; }
    x_interp.UseDevice(true);
 
-   CheckFESpace(f_ordering);
+   CheckFESpace();
 
-   ComputeScatterIndicesAndOffsets(f_ordering, type);
+   ComputeScatterIndicesAndOffsets();
 
-   ComputeGatherIndices(f_ordering, type);
+   ComputeGatherIndices();
 }
 
 NCL2FaceRestriction::NCL2FaceRestriction(const FiniteElementSpace &fes,
@@ -2259,18 +2169,16 @@ int ToLexOrdering(const int dim, const int face_id, const int size1d,
       case 1:
          return 0;
       case 2:
-         return ToLexOrdering2D(face_id, size1d, index);
+         return internal::ToLexOrdering2D(face_id, size1d, index);
       case 3:
-         return ToLexOrdering3D(face_id, size1d, index%size1d, index/size1d);
+         return internal::ToLexOrdering3D(face_id, size1d, index%size1d, index/size1d);
       default:
          MFEM_ABORT("Unsupported dimension.");
          return 0;
    }
 }
 
-void NCL2FaceRestriction::ComputeScatterIndicesAndOffsets(
-   const ElementDofOrdering f_ordering,
-   const FaceType type)
+void NCL2FaceRestriction::ComputeScatterIndicesAndOffsets()
 {
    Mesh &mesh = *fes.GetMesh();
 
@@ -2334,9 +2242,7 @@ void NCL2FaceRestriction::ComputeScatterIndicesAndOffsets(
    interpolations.InitializeNCInterpConfig();
 }
 
-void NCL2FaceRestriction::ComputeGatherIndices(
-   const ElementDofOrdering f_ordering,
-   const FaceType type)
+void NCL2FaceRestriction::ComputeGatherIndices()
 {
    Mesh &mesh = *fes.GetMesh();
    // Computation of gather_indices

fem/restriction.hpp

@@ -14,6 +14,7 @@
 
 #include "../linalg/operator.hpp"
 #include "../mesh/mesh.hpp"
+#include "normal_deriv_restriction.hpp"
 
 namespace mfem
 {
@@ -21,6 +22,8 @@ namespace mfem
 class FiniteElementSpace;
 enum class ElementDofOrdering;
 
+class FaceQuadratureSpace;
+
 /// Abstract base class that defines an interface for element restrictions.
 class ElementRestrictionOperator : public Operator
 {
@@ -220,6 +223,45 @@ public:
       y = 0.0;
       AddMultTranspose(x, y);
    }
+
+   /** @brief For each face, sets @a y to the partial derivative of @a x with
+              respect to the reference coordinate whose direction is
+              perpendicular to the face on the reference element.
+
+    @details This is not the normal derivative in physical coordinates, but can
+             be mapped to the physical normal derivative using the element
+             Jacobian and the tangential derivatives (in reference coordinates)
+             which can be computed from the face values (provided by Mult).
+
+             Note that due to the polynomial degree of the element mapping, the
+             physical normal derivative may be a higher degree polynomial than
+             the restriction of the values to the face. However, the normal
+             derivative in reference coordinates has degree-1, and therefore can
+             be exactly represented with the degrees of freedom of a face
+             E-vector.
+
+    @param[in]     x The L-vector degrees of freedom.
+    @param[in,out] y The reference normal derivative degrees of freedom. Is
+                     E-vector like.
+    */
+   virtual void NormalDerivativeMult(const Vector &x, Vector &y) const
+   {
+      MFEM_ABORT("Not implemented for this restriction operator.");
+   }
+
+   /** @brief Add the face reference-normal derivative degrees of freedom in @a
+              x to the element degrees of freedom in @a y.
+
+       @details see NormalDerivativeMult.
+
+       @param[in]     x The degrees of freedom of the face reference-normal
+                        derivative. Is E-vector like.
+       @param[in,out] y The L-vector degrees of freedom.
+   */
+   virtual void NormalDerivativeAddMultTranspose(const Vector &x, Vector &y) const
+   {
+      MFEM_ABORT("Not implemented for this restriction operator.");
+   }
 };
 
 /// @brief Operator that extracts face degrees of freedom for H1, ND, or RT
@@ -365,6 +407,7 @@ class L2FaceRestriction : public FaceRestriction
 {
 protected:
    const FiniteElementSpace &fes;
+   const ElementDofOrdering ordering;
    const int nf; // Number of faces of the requested type
    const int ne; // Number of elements
    const int vdim; // vdim
@@ -379,6 +422,7 @@ protected:
    Array<int> scatter_indices2; // Scattering indices for element 2 on each face
    Array<int> gather_offsets; // offsets for the gathering indices of each dof
    Array<int> gather_indices; // gathering indices for each dof
+   mutable std::unique_ptr<L2NormalDerivativeFaceRestriction> normal_deriv_restr;
 
    /** @brief Constructs an L2FaceRestriction.
 
@@ -487,34 +531,49 @@ public:
    virtual void AddFaceMatricesToElementMatrices(const Vector &fea_data,
                                                  Vector &ea_data) const;
 
+   /** @brief Scatter the degrees of freedom, i.e. goes from L-Vector to
+      face E-Vector.
+
+      @param[in]  x The L-vector degrees of freedom.
+      @param[out] y The face E-Vector degrees of freedom with the given format:
+                    (face_dofs x vdim x 2 x nf) where nf is the number of
+                    interior or boundary faces requested by @a type in the
+                    constructor. The face_dofs are ordered according to the
+                    given ElementDofOrdering. */
+   void NormalDerivativeMult(const Vector &x, Vector &y) const override;
+
+   /** @brief Add the face reference-normal derivative degrees of freedom in @a
+              x to the element degrees of freedom in @a y.
+
+       @details see NormalDerivativeMult.
+
+       @param[in]     x The degrees of freedom of the face reference-normal
+                        derivative. Is E-vector like.
+       @param[in,out] y The L-vector degrees of freedom.
+   */
+   void NormalDerivativeAddMultTranspose(const Vector &x,
+                                         Vector &y) const override;
 private:
    /** @brief Compute the scatter indices: L-vector to E-vector, and the offsets
        for the gathering: E-vector to L-vector.
-
-       @param[in] f_ordering Request a specific face dof ordering.
-       @param[in] type       Request internal or boundary faces dofs.
    */
-   void ComputeScatterIndicesAndOffsets(const ElementDofOrdering f_ordering,
-                                        const FaceType type);
+   void ComputeScatterIndicesAndOffsets();
 
    /** @brief Compute the gather indices: E-vector to L-vector.
 
        Note: Requires the gather offsets to be computed.
-
-       @param[in] f_ordering Request a specific face dof ordering.
-       @param[in] type       Request internal or boundary faces dofs.
    */
-   void ComputeGatherIndices(const ElementDofOrdering f_ordering,
-                             const FaceType type);
+   void ComputeGatherIndices();
+
+   /// Create the internal normal derivative restriction operator if needed.
+   void EnsureNormalDerivativeRestriction() const;
 
 protected:
    mutable Array<int> face_map; // Used in the computation of GetFaceDofs
 
    /** @brief Verify that L2FaceRestriction is built from an L2 FESpace.
-
-       @param[in] f_ordering The requested face dof ordering.
    */
-   void CheckFESpace(const ElementDofOrdering f_ordering);
+   void CheckFESpace();
 
    /** @brief Set the scattering indices of elem1, and increment the offsets for
        the face described by the @a face. The ordering of the face dofs of elem1
@@ -938,22 +997,14 @@ private:
    /** @brief Compute the scatter indices: L-vector to E-vector, the offsets
        for the gathering: E-vector to L-vector, and the interpolators from
        coarse to fine face for master non-comforming faces.
-
-       @param[in] f_ordering Request a specific face dof ordering.
-       @param[in] type       Request internal or boundary faces dofs.
    */
-   void ComputeScatterIndicesAndOffsets(const ElementDofOrdering f_ordering,
-                                        const FaceType type);
+   void ComputeScatterIndicesAndOffsets();
 
    /** @brief Compute the gather indices: E-vector to L-vector.
 
        Note: Requires the gather offsets to be computed.
-
-       @param[in] f_ordering Request a specific face dof ordering.
-       @param[in] type       Request internal or boundary faces dofs.
    */
-   void ComputeGatherIndices(const ElementDofOrdering f_ordering,
-                             const FaceType type);
+   void ComputeGatherIndices();
 
 public:
    /** @brief Scatter the degrees of freedom, i.e. goes from L-Vector to
@@ -1015,7 +1066,6 @@ public:
    void DoubleValuedNonconformingTransposeInterpolationInPlace(Vector& x) const;
 };
 
-
 /** @brief Convert a dof face index from Native ordering to lexicographic
     ordering for quads and hexes.
 

general/table.hpp

@@ -95,7 +95,7 @@ public:
        not called, it returns the number of possible connections established
        by the used constructor. Otherwise, it is exactly the number of
        established connections before calling Finalize(). */
-   inline int Size_of_connections() const { return I[size]; }
+   inline int Size_of_connections() const { HostReadI(); return I[size]; }
 
    /** Returns index of the connection between element i of TYPE I and
        element j of TYPE II. If there is no connection between element i

mesh/CMakeLists.txt

@@ -12,6 +12,7 @@
 set(SRCS
   attribute_sets.cpp
   element.cpp
+  face_nbr_geom.cpp
   gmsh.cpp
   hexahedron.cpp
   mesh.cpp
@@ -37,6 +38,7 @@ set(SRCS
 set(HDRS
   attribute_sets.hpp
   element.hpp
+  face_nbr_geom.hpp
   gmsh.hpp
   hexahedron.hpp
   mesh.hpp

mesh/face_nbr_geom.cpp

@@ -0,0 +1,127 @@
+// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#include "face_nbr_geom.hpp"
+#include "../general/forall.hpp"
+
+namespace mfem
+{
+
+FaceNeighborGeometricFactors::FaceNeighborGeometricFactors(
+   const GeometricFactors &geom_) : num_neighbor_elems(0), geom(geom_)
+{
+#ifdef MFEM_USE_MPI
+   if (const ParMesh *par_mesh = dynamic_cast<const ParMesh*>(geom.mesh))
+   {
+      const int flags = geom.computed_factors;
+      const int dim = par_mesh->Dimension();
+      const int sdim = par_mesh->SpaceDimension();
+
+      const_cast<ParMesh*>(par_mesh)->ExchangeFaceNbrData();
+
+      num_neighbor_elems = par_mesh->GetNFaceNeighborElements();
+
+      if (flags & GeometricFactors::COORDINATES)
+      {
+         ExchangeFaceNbrQVectors(geom.X, X, sdim);
+      }
+      if (flags & GeometricFactors::JACOBIANS)
+      {
+         ExchangeFaceNbrQVectors(geom.J, J, dim*sdim);
+      }
+      if (flags & GeometricFactors::DETERMINANTS)
+      {
+         ExchangeFaceNbrQVectors(geom.detJ, detJ, 1);
+      }
+
+      // Free memory of work arrays.
+      send_data.Destroy();
+      send_offsets.DeleteAll();
+      recv_offsets.DeleteAll();
+   }
+#endif
+}
+
+void FaceNeighborGeometricFactors::ExchangeFaceNbrQVectors(
+   const Vector &x_local, Vector &x_shared, const int vdim)
+{
+#ifdef MFEM_USE_MPI
+
+   const int nq = geom.IntRule->Size();
+   const int ndof_per_el = vdim * nq;
+
+   const ParMesh *mesh = static_cast<const ParMesh*>(geom.mesh);
+
+   const int n_face_nbr = mesh->GetNFaceNeighbors();
+   if (n_face_nbr == 0) { return; }
+
+   const int ne_send = mesh->send_face_nbr_elements.Size_of_connections();
+
+   x_shared.SetSize(ndof_per_el * num_neighbor_elems);
+   send_offsets.SetSize(n_face_nbr + 1);
+   send_data.SetSize(ndof_per_el * ne_send);
+   auto h_send_data = Reshape(send_data.HostWrite(), ndof_per_el, ne_send);
+   const auto h_x_local = Reshape(x_local.HostRead(), ndof_per_el, mesh->GetNE());
+   int idx = 0;
+   Array<int> row;
+   for (int i = 0; i < n_face_nbr; ++i)
+   {
+      send_offsets[i] = ndof_per_el*idx;
+      mesh->send_face_nbr_elements.GetRow(i, row);
+      for (const int el : row)
+      {
+         for (int j = 0; j < ndof_per_el; ++j)
+         {
+            h_send_data(j, idx) = h_x_local(j, el);
+         }
+         ++idx;
+      }
+   }
+   send_offsets[n_face_nbr] = ndof_per_el*idx;
+   MFEM_ASSERT(send_offsets[n_face_nbr] == send_data.Size(), "");
+
+   recv_offsets.SetSize(n_face_nbr + 1);
+   for (int i = 0; i < n_face_nbr + 1; ++i)
+   {
+      recv_offsets[i] = ndof_per_el * mesh->face_nbr_elements_offset[i];
+   }
+   MFEM_ASSERT(recv_offsets[n_face_nbr] == x_shared.Size(), "");
+
+   MPI_Comm comm = mesh->GetComm();
+   std::vector<MPI_Request> send_reqs(n_face_nbr);
+   std::vector<MPI_Request> recv_reqs(n_face_nbr);
+   std::vector<MPI_Status> statuses(n_face_nbr);
+
+   bool mpi_gpu_aware = Device::GetGPUAwareMPI();
+   const auto send_data_ptr = mpi_gpu_aware ? send_data.Read() :
+                              send_data.HostRead();
+   auto x_shared_ptr = mpi_gpu_aware ? x_shared.Write() : x_shared.HostWrite();
+
+   for (int i = 0; i < n_face_nbr; ++i)
+   {
+      const int nbr_rank = mesh->GetFaceNbrRank(i);
+      const int tag = 0;
+
+      MPI_Isend(send_data_ptr + send_offsets[i],
+                send_offsets[i+1] - send_offsets[i],
+                MPITypeMap<real_t>::mpi_type, nbr_rank, tag, comm, &send_reqs[i]);
+
+      MPI_Irecv(x_shared_ptr + recv_offsets[i],
+                recv_offsets[i+1] - recv_offsets[i],
+                MPITypeMap<real_t>::mpi_type, nbr_rank, tag, comm, &recv_reqs[i]);
+   }
+
+   MPI_Waitall(n_face_nbr, send_reqs.data(), statuses.data());
+   MPI_Waitall(n_face_nbr, recv_reqs.data(), statuses.data());
+#endif
+}
+
+} // namespace mfem

mesh/face_nbr_geom.hpp

@@ -0,0 +1,65 @@
+// Copyright (c) 2010-2024, Lawrence Livermore National Security, LLC. Produced
+// at the Lawrence Livermore National Laboratory. All Rights reserved. See files
+// LICENSE and NOTICE for details. LLNL-CODE-806117.
+//
+// This file is part of the MFEM library. For more information and source code
+// availability visit https://mfem.org.
+//
+// MFEM is free software; you can redistribute it and/or modify it under the
+// terms of the BSD-3 license. We welcome feedback and contributions, see file
+// CONTRIBUTING.md for details.
+
+#ifndef MFEM_FACE_NBR_GEOM
+#define MFEM_FACE_NBR_GEOM
+
+#include "../config/config.hpp"
+#include "mesh.hpp"
+#include "pmesh.hpp"
+
+namespace mfem
+{
+
+/// @brief Class for accessing the geometric factors of face neighbor elements
+/// (i.e. across boundaries of MPI mesh partitions).
+///
+/// @sa GeometricFactors
+class FaceNeighborGeometricFactors
+{
+public:
+   int num_neighbor_elems; ///< Number of face neighbor elements.
+
+   /// @name Geometric factor data arrays
+   /// These are stored with layout (NQ, VDIM, NE). See the documentation of
+   /// GeometricFactors for more details.
+   ///@{
+
+   Vector X; ///< Physical coordinates of the mesh.
+   Vector J; ///< Jacobian matrices
+   Vector detJ; ///< Jacobian determinants
+
+   ///@}
+
+   /// Communicate (if needed) to gather the face neighbor geometric factors.
+   FaceNeighborGeometricFactors(const GeometricFactors &geom_);
+
+protected:
+   const GeometricFactors &geom; ///< The GeometricFactors of the Mesh.
+
+   /// @name Internal work arrays, used for MPI communication
+   ///@{
+
+   Vector send_data;
+   Array<int> send_offsets, recv_offsets;
+
+   ///@}
+
+   /// @brief Given a Q-vector @a x_local with @a vdim components, fill the
+   /// face-neighbor Q-vector @a x_shared by communicating with neighboring MPI
+   /// partitions.
+   void ExchangeFaceNbrQVectors(const Vector &x_local, Vector &x_shared,
+                                const int vdim);
+};
+
+} // namespace mfem
+
+#endif

mesh/mesh.cpp

@@ -1484,6 +1484,34 @@ Array<int> Mesh::GetFaceToBdrElMap() const
    return face_to_be;
 }
 
+Geometry::Type Mesh::GetTypicalElementGeometry() const
+{
+   if (GetNE() > 0) { return GetElementGeometry(0); }
+
+   const int dim = Dimension();
+   if (dim == 1)
+   {
+      return Geometry::SEGMENT;
+   }
+   const int mesh_gen = MeshGenerator();
+   Geometry::Type geom = Geometry::INVALID;
+   if (dim == 2)
+   {
+      geom = ((mesh_gen & 1) ? Geometry::TRIANGLE :
+              ((mesh_gen & 2) ? Geometry::SQUARE : Geometry::INVALID));
+   }
+   else if (dim == 3)
+   {
+      geom = ((mesh_gen & 1) ? Geometry::TETRAHEDRON :
+              ((mesh_gen & 2) ? Geometry::CUBE :
+               ((mesh_gen & 4) ? Geometry::PRISM :
+                ((mesh_gen & 8) ? Geometry::PYRAMID : Geometry::INVALID))));
+   }
+   MFEM_VERIFY(geom != Geometry::INVALID,
+               "Could not determine a typical element Geometry!");
+   return geom;
+}
+
 void Mesh::Init()
 {
    // in order of declaration:
@@ -13329,7 +13357,7 @@ void GeometricFactors::Compute(const GridFunction &nodes,
 {
 
    const FiniteElementSpace *fespace = nodes.FESpace();
-   const FiniteElement *fe = fespace->GetFE(0);
+   const FiniteElement *fe = fespace->GetTypicalFE();
    const int dim  = fe->GetDim();
    const int vdim = fespace->GetVDim();
    const int NE   = fespace->GetNE();

mesh/mesh.hpp

@@ -1370,6 +1370,13 @@ public:
       return elements[i]->GetGeometryType();
    }
 
+   /** @brief If the local mesh is not empty, return GetElementGeometry(0);
+       otherwise, return a typical Geometry present in the global mesh.
+
+       This method can be used to replace calls like GetElementGeometry(0) in
+       order to handle empty local meshes better. */
+   Geometry::Type GetTypicalElementGeometry() const;
+
    Geometry::Type GetBdrElementGeometry(int i) const
    {
       return boundary[i]->GetGeometryType();

tests/unit/fem/test_pa_kernels.cpp

@@ -713,4 +713,120 @@ TEST_CASE("PA Boundary Mass", "[PartialAssembly], [CUDA]")
    REQUIRE(y_fa.Normlinf() == MFEM_Approx(0.0));
 }
 
+namespace
+{
+template <typename T> struct ParTypeHelper { };
+template <> struct ParTypeHelper<FiniteElementSpace>
+{
+   using GF_t = GridFunction;
+   using BLF_t = BilinearForm;
+};
+#ifdef MFEM_USE_MPI
+template <> struct ParTypeHelper<ParFiniteElementSpace>
+{
+   using GF_t = ParGridFunction;
+   using BLF_t = ParBilinearForm;
+};
+#endif
+}
+
+template <typename FES>
+void test_dg_diffusion(FES &fes)
+{
+   using GF_t = typename ParTypeHelper<FES>::GF_t;
+   using BLF_t = typename ParTypeHelper<FES>::BLF_t;
+
+   GF_t x(&fes), y_fa(&fes), y_pa(&fes);
+   x.Randomize(1);
+
+   ConstantCoefficient pi(3.14159);
+
+   const real_t sigma = -1.0;
+   const real_t kappa = 10.0;
+
+   IntegrationRules irs(0, Quadrature1D::GaussLobatto);
+   const IntegrationRule &ir = irs.Get(fes.GetMesh()->GetFaceGeometry(0),
+                                       2*fes.GetMaxElementOrder());
+
+   BLF_t blf_fa(&fes);
+   blf_fa.AddInteriorFaceIntegrator(new DGDiffusionIntegrator(pi, sigma, kappa));
+   blf_fa.AddBdrFaceIntegrator(new DGDiffusionIntegrator(pi, sigma, kappa));
+   (*blf_fa.GetFBFI())[0]->SetIntegrationRule(ir);
+   (*blf_fa.GetBFBFI())[0]->SetIntegrationRule(ir);
+   blf_fa.Assemble();
+   blf_fa.Finalize();
+   OperatorHandle A_fa;
+   Array<int> empty;
+   blf_fa.FormSystemMatrix(empty, A_fa);
+   A_fa->Mult(x, y_fa);
+
+   BLF_t blf_pa(&fes);
+   blf_pa.SetAssemblyLevel(AssemblyLevel::PARTIAL);
+   blf_pa.AddInteriorFaceIntegrator(new DGDiffusionIntegrator(pi, sigma, kappa));
+   blf_pa.AddBdrFaceIntegrator(new DGDiffusionIntegrator(pi, sigma, kappa));
+   (*blf_pa.GetFBFI())[0]->SetIntegrationRule(ir);
+   (*blf_pa.GetBFBFI())[0]->SetIntegrationRule(ir);
+   blf_pa.Assemble();
+   blf_pa.Mult(x, y_pa);
+
+   y_fa -= y_pa;
+
+   REQUIRE(y_fa.Normlinf() == MFEM_Approx(0.0));
+}
+
+std::vector<std::string> get_dg_test_meshes()
+{
+   std::vector<std::string> mesh_filenames =
+   {
+      "../../data/star.mesh",
+      "../../data/star-q3.mesh",
+      "../../data/fichera.mesh",
+      "../../data/fichera-q3.mesh",
+   };
+   const bool have_data_dir = mfem_data_dir != "";
+   if (have_data_dir)
+   {
+      mesh_filenames.push_back(mfem_data_dir + "/gmsh/v22/unstructured_quad.v22.msh");
+      mesh_filenames.push_back(mfem_data_dir + "/gmsh/v22/unstructured_hex.v22.msh");
+   }
+   return mesh_filenames;
+}
+
+TEST_CASE("PA DG Diffusion", "[PartialAssembly], [CUDA]")
+{
+   const auto mesh_fname = GENERATE_COPY(from_range(get_dg_test_meshes()));
+   const int order = GENERATE(1, 2);
+   CAPTURE(order, mesh_fname);
+
+   Mesh mesh = Mesh::LoadFromFile(mesh_fname.c_str());
+   const int dim = mesh.Dimension();
+
+   DG_FECollection fec(order, dim, BasisType::GaussLobatto);
+   FiniteElementSpace fes(&mesh, &fec);
+
+   test_dg_diffusion(fes);
+}
+
+#ifdef MFEM_USE_MPI
+
+TEST_CASE("Parallel PA DG Diffusion", "[PartialAssembly][Parallel][CUDA]")
+{
+   const auto mesh_fname = GENERATE_COPY(from_range(get_dg_test_meshes()));
+   const int order = GENERATE(1, 2);
+   CAPTURE(order, mesh_fname);
+
+   Mesh serial_mesh = Mesh::LoadFromFile(mesh_fname.c_str());
+   ParMesh mesh(MPI_COMM_WORLD, serial_mesh);
+   serial_mesh.Clear();
+
+   const int dim = mesh.Dimension();
+
+   DG_FECollection fec(order, dim, BasisType::GaussLobatto);
+   ParFiniteElementSpace fes(&mesh, &fec);
+
+   test_dg_diffusion(fes);
+}
+
+#endif
+
 } // namespace pa_kernels