address pedantic self-assign warning in macos-latest-make-opt-seq-int32 CI test

examples/CMakeLists.txt

@@ -40,6 +40,8 @@ list(APPEND ALL_EXE_SRCS
   ex30.cpp
   ex31.cpp
   ex33.cpp
+  ex34.cpp
+  ex35.cpp
   )
 
 if (MFEM_USE_MPI)

examples/ex34.cpp

@@ -0,0 +1,412 @@
+//                                MFEM Example 34
+//
+// Compile with: make ex34
+//
+// Sample runs:  ex34
+//
+// Description:  This example code demonstrates the use of MFEM to define a
+//               discontinuous Galerkin (DG) finite element discretization of
+//               the Laplace problem -Delta u = f with Dirichlet boundary
+//               conditions. Finite element spaces of any order, including zero
+//               on regular grids, are supported. The example highlights the
+//               use of coupling solution domains though custom physics defined
+//               on internal boundaries.
+//
+//               We recommend viewing examples 1 and 14 before viewing this
+//               example.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+class InteriorLFIntegrator : public LinearFormIntegrator
+{
+public:
+   InteriorLFIntegrator(Coefficient &Q)
+      : Q(Q)
+   {}
+
+   void AssembleRHSElementVect(const FiniteElement &el1,
+                               const FiniteElement &el2,
+                               FaceElementTransformations &trans,
+                               Vector &mesh_coords_bar) override;
+
+   void AssembleRHSElementVect(const FiniteElement &el,
+                               ElementTransformation &trans,
+                               Vector &elvect) override
+   {
+      mfem_error("AssembleRHSElementVect(...)");
+   }
+
+private:
+   Coefficient &Q;
+
+#ifndef MFEM_THREAD_SAFE
+   Vector shape1;
+   Vector shape2;
+#endif
+};
+
+Mesh generate_mesh(int ref, int internal_bdr_attr = 5);
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   int ref_levels = 0;
+   int order = 1;
+   int sol_order = 3;
+   double jump = -2;
+   double sigma = -1.0;
+   double kappa = -1.0;
+   double eta = 0.0;
+   bool visualization = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&ref_levels, "-r", "--refine",
+                  "Number of times to refine the mesh uniformly, -1 for auto.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) >= 0.");
+   args.AddOption(&sigma, "-s", "--sigma",
+                  "One of the three DG penalty parameters, typically +1/-1."
+                  " See the documentation of class DGDiffusionIntegrator.");
+   args.AddOption(&kappa, "-k", "--kappa",
+                  "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(&sol_order, "-so", "--solution_order",
+                  "Polynomial order of the exact solution >= 0.");
+   args.AddOption(&jump, "-j", "--jump",
+                  "Value of the discontinuity between the material regions.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   if (kappa < 0)
+   {
+      kappa = (order+1)*(order+1);
+   }
+   if (sol_order < 0)
+   {
+      sol_order = 1;
+   }
+   args.PrintOptions(cout);
+
+   // 2. Construct the (serial) mesh and refine it if requested.
+   auto mesh = generate_mesh(ref_levels);
+
+   int dim = mesh.Dimension();
+
+   if (mesh.NURBSext)
+   {
+      mesh.SetCurvature(max(order, 1));
+   }
+
+   // 3. Define a finite element space on the mesh. Here we use discontinuous
+   //    finite elements of the specified order >= 0.
+   DG_FECollection fec(order, dim);
+   FiniteElementSpace fespace(&mesh, &fec);
+   cout << "Number of unknowns: " << fespace.GetVSize() << endl;
+
+   // 4. Set up the linear form b(.) which corresponds to the right-hand side of
+   //    the FEM linear system.
+   LinearForm b(&fespace);
+
+   Array<int> p1_attr_marker(mesh.attributes.Max());
+   p1_attr_marker = 0;
+   p1_attr_marker[0] = 1;
+
+   FunctionCoefficient p1_source([sol_order](const Vector &p)
+   {
+      const double x = p(0);
+      const double val = -(sol_order - 1)*sol_order*pow(x, sol_order-2);
+      return val;
+   });
+   b.AddDomainIntegrator(new DomainLFIntegrator(p1_source), p1_attr_marker);
+
+   Array<int> p2_attr_marker(mesh.attributes.Max());
+   p2_attr_marker = 0;
+   p2_attr_marker[1] = 1;
+
+   FunctionCoefficient p2_source([sol_order](const Vector &p)
+   {
+      const double x = p(0);
+      double val = -(sol_order - 1)*sol_order*pow(x - 2, sol_order-2);
+      if (sol_order % 2 == 0)
+      {
+         val *= -1.0;
+      }
+      return val;
+   });
+   b.AddDomainIntegrator(new DomainLFIntegrator(p2_source), p2_attr_marker);
+
+   ConstantCoefficient one(1.0);
+
+   Array<int> p1_bdr_attr_marker(mesh.bdr_attributes.Max());
+   p1_bdr_attr_marker = 0;
+   p1_bdr_attr_marker[0] = 1;
+
+   ConstantCoefficient left_bc_val(0.0);
+   b.AddBdrFaceIntegrator(
+      new DGDirichletLFIntegrator(left_bc_val, one, sigma, kappa),
+      p1_bdr_attr_marker);
+
+   Array<int> p2_bdr_attr_marker(mesh.bdr_attributes.Max());
+   p2_bdr_attr_marker = 0;
+   p2_bdr_attr_marker[1] = 1;
+
+   ConstantCoefficient right_bc_val(2.0 + jump);
+   b.AddBdrFaceIntegrator(
+      new DGDirichletLFIntegrator(right_bc_val, one, sigma, kappa),
+      p2_bdr_attr_marker);
+
+   Array<int> internal_bdr_attr_marker(mesh.bdr_attributes.Max());
+   internal_bdr_attr_marker = 0;
+   internal_bdr_attr_marker[4] = 1;
+
+   ConstantCoefficient interface_flux(sol_order);
+   b.AddInternalBoundaryFaceIntegrator(
+      new InteriorLFIntegrator(interface_flux),
+      internal_bdr_attr_marker);
+
+   b.Assemble();
+
+   // 5. Define the solution vector x as a finite element grid function
+   //    corresponding to fespace. Initialize x with initial guess of zero.
+   GridFunction x(&fespace);
+   x = 0.0;
+
+   // 6. 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(&fespace);
+   a.AddDomainIntegrator(new DiffusionIntegrator(one));
+
+   a.AddInteriorFaceIntegrator(new DGDiffusionIntegrator(one, sigma, kappa));
+   a.AddBdrFaceIntegrator(new DGDiffusionIntegrator(one, sigma, kappa),
+                          p1_bdr_attr_marker);
+   a.AddBdrFaceIntegrator(new DGDiffusionIntegrator(one, sigma, kappa),
+                          p2_bdr_attr_marker);
+   if (eta > 0)
+   {
+      a.AddInteriorFaceIntegrator(new DGDiffusionBR2Integrator(fespace, eta));
+      a.AddBdrFaceIntegrator(new DGDiffusionBR2Integrator(fespace, eta));
+   }
+
+   // 7. Negate the DG interface terms along the internal boundary so that the
+   //    only coupling between domains is from the chosen model (constant flux
+   //    in this case).
+   ProductCoefficient neg_one(-1.0, one);
+   a.AddInternalBoundaryFaceIntegrator(new DGDiffusionIntegrator(neg_one, sigma,
+                                                                 kappa),
+                                       internal_bdr_attr_marker);
+   if (eta > 0)
+   {
+      a.AddInternalBoundaryFaceIntegrator(new DGDiffusionBR2Integrator(fespace,
+                                                                       neg_one, eta),
+                                          internal_bdr_attr_marker);
+   }
+
+   a.Assemble();
+   a.Finalize();
+   const SparseMatrix &A = a.SpMat();
+
+#ifndef MFEM_USE_SUITESPARSE
+   // 8. 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)
+   {
+      PCG(A, M, b, x, 1, 500, 1e-12, 0.0);
+   }
+   else
+   {
+      GMRES(A, M, b, x, 1, 500, 500, 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);
+#endif
+
+   // 9. 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);
+   ofstream sol_ofs("sol.gf");
+   sol_ofs.precision(8);
+   x.Save(sol_ofs);
+
+   // 10. 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;
+   }
+
+   return 0;
+}
+
+void InteriorLFIntegrator::AssembleRHSElementVect(
+   const FiniteElement &el1,
+   const FiniteElement &el2,
+   FaceElementTransformations &trans,
+   Vector &elvect)
+{
+   int ndof1 = el1.GetDof();
+   int ndof2 = el2.GetDof();
+   int ndof = ndof1 + ndof2;
+
+#ifdef MFEM_THREAD_SAFE
+   Vector shape1;
+   Vector shape2;
+#endif
+   shape1.SetSize(ndof1);
+   shape2.SetSize(ndof2);
+
+   const auto *ir = IntRule;
+   if (ir == NULL)
+   {
+      int order = 2 * max(el1.GetOrder(), el2.GetOrder());
+      ir = &IntRules.Get(trans.GetGeometryType(), order);
+   }
+
+   elvect.SetSize(ndof);
+   Vector elvect1(elvect.GetData(), ndof1);
+   Vector elvect2(elvect.GetData() + ndof1, ndof2);
+   elvect = 0.0;
+   for (int i = 0; i < ir->GetNPoints(); i++)
+   {
+      const auto &ip = ir->IntPoint(i);
+
+      // Set the integration point in the face and the neighboring element
+      trans.SetAllIntPoints(&ip);
+
+      const double w = ip.weight * trans.Weight();
+
+      // Access the neighboring element's integration point
+      const auto &eip1 = trans.GetElement1IntPoint();
+      const auto &eip2 = trans.GetElement2IntPoint();
+
+      double Q_val = Q.Eval(trans, ip);
+
+      el1.CalcShape(eip1, shape1);
+      el2.CalcShape(eip2, shape2);
+
+      elvect1.Add(Q_val * w, shape1);
+      elvect2.Add(-Q_val * w, shape2);
+   }
+}
+
+Mesh generate_mesh(int ref, int internal_bdr_attr)
+{
+   int nxy = 4 * (ref+1);
+   auto mesh = Mesh::MakeCartesian2D(nxy, nxy, Element::TRIANGLE, true, 2.0, 1.0);
+   // auto mesh = Mesh::MakeCartesian2D(nxy, nxy, Element::QUADRILATERAL, true, 2.0, 1.0);
+
+   // assign element attributes to left and right sides
+   for (int i = 0; i < mesh.GetNE(); ++i)
+   {
+      auto *elem = mesh.GetElement(i);
+
+      Array<int> verts;
+      elem->GetVertices(verts);
+
+      bool left = true;
+      for (int j = 0; j < verts.Size(); ++j)
+      {
+         auto *vtx = mesh.GetVertex(verts[j]);
+         if (vtx[0] <= 1.0)
+         {
+            continue;
+         }
+         else
+         {
+            left = false;
+         }
+      }
+      if (left)
+      {
+         elem->SetAttribute(1);
+      }
+      else
+      {
+         elem->SetAttribute(2);
+      }
+   }
+
+   // assign boundary element attributes to left and right sides
+   for (int i = 0; i < mesh.GetNBE(); ++i)
+   {
+      auto *elem = mesh.GetBdrElement(i);
+
+      Array<int> verts;
+      elem->GetVertices(verts);
+
+      bool left = true;
+      bool right = true;
+      bool top = true;
+      bool bottom = true;
+      for (int j = 0; j < verts.Size(); ++j)
+      {
+         auto *vtx = mesh.GetVertex(verts[j]);
+         left = left && abs(vtx[0] - 0.0) < 1e-12;
+         right = right && abs(vtx[0] - 2.0) < 1e-12;
+         top = top && abs(vtx[1] - 1.0) < 1e-12;
+         bottom = bottom && abs(vtx[1] - 0.0) < 1e-12;
+      }
+      if (left)
+      {
+         elem->SetAttribute(1);
+      }
+      else if (right)
+      {
+         elem->SetAttribute(2);
+      }
+      else if (top)
+      {
+         elem->SetAttribute(3);
+      }
+      else if (bottom)
+      {
+         elem->SetAttribute(4);
+      }
+   }
+
+   // add internal boundary elements
+   for (int i = 0; i < mesh.GetNumFaces(); ++i)
+   {
+      int e1, e2;
+      mesh.GetFaceElements(i, &e1, &e2);
+      if (e1 >= 0 && e2 >= 0 && mesh.GetAttribute(e1) != mesh.GetAttribute(e2))
+      {
+         // This is the internal face between attributes.
+         auto *new_elem = mesh.GetFace(i)->Duplicate(&mesh);
+         new_elem->SetAttribute(internal_bdr_attr);
+         mesh.AddBdrElement(new_elem);
+      }
+   }
+
+   mesh.FinalizeTopology(); // Finalize to build relevant tables
+   mesh.Finalize();
+   mesh.SetAttributes();
+
+   return mesh;
+}

examples/ex35.cpp

@@ -0,0 +1,451 @@
+//                                MFEM Example 36
+//
+// Compile with: make ex36
+//
+// Sample runs:  ex36
+//
+// Description:  This example code demonstrates the use of MFEM to define a
+//               discontinuous Galerkin (DG) finite element discretization of
+//               the Laplace problem -Delta u = f with Dirichlet boundary
+//               conditions. Finite element spaces of any order, including zero
+//               on regular grids, are supported. The example highlights the
+//               use of coupling solution domains though custom physics defined
+//               on internal boundaries.
+//
+//               We recommend viewing examples 1, 14, and 34 before viewing this
+//               example.
+
+#include "mfem.hpp"
+#include <fstream>
+#include <iostream>
+
+using namespace std;
+using namespace mfem;
+
+class InteriorMassIntegrator : public BilinearFormIntegrator
+{
+public:
+   InteriorMassIntegrator(Coefficient &Q)
+      : Q(Q)
+   {}
+
+   void AssembleFaceMatrix(const FiniteElement &el1,
+                           const FiniteElement &el2,
+                           FaceElementTransformations &trans,
+                           DenseMatrix &elmat) override;
+
+   using BilinearFormIntegrator::AssembleFaceMatrix;
+
+private:
+   Coefficient &Q;
+
+#ifndef MFEM_THREAD_SAFE
+   Vector shape1;
+   Vector shape2;
+   DenseMatrix elmat11;
+   DenseMatrix elmat12;
+   DenseMatrix elmat21;
+   DenseMatrix elmat22;
+#endif
+};
+
+Mesh generate_mesh(int ref, int internal_bdr_attr = 5);
+
+int main(int argc, char *argv[])
+{
+   // 1. Parse command-line options.
+   int ref_levels = 0;
+   int order = 1;
+   int sol_order = 3;
+   double jump = -2;
+   double sigma = -1.0;
+   double kappa = -1.0;
+   double eta = 0.0;
+   bool visualization = 1;
+
+   OptionsParser args(argc, argv);
+   args.AddOption(&ref_levels, "-r", "--refine",
+                  "Number of times to refine the mesh uniformly, -1 for auto.");
+   args.AddOption(&order, "-o", "--order",
+                  "Finite element order (polynomial degree) >= 0.");
+   args.AddOption(&sigma, "-s", "--sigma",
+                  "One of the three DG penalty parameters, typically +1/-1."
+                  " See the documentation of class DGDiffusionIntegrator.");
+   args.AddOption(&kappa, "-k", "--kappa",
+                  "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(&sol_order, "-so", "--solution_order",
+                  "Polynomial order of the exact solution >= 0.");
+   args.AddOption(&jump, "-j", "--jump",
+                  "Value of the discontinuity between the material regions.");
+   args.AddOption(&visualization, "-vis", "--visualization", "-no-vis",
+                  "--no-visualization",
+                  "Enable or disable GLVis visualization.");
+   args.Parse();
+   if (!args.Good())
+   {
+      args.PrintUsage(cout);
+      return 1;
+   }
+   if (kappa < 0)
+   {
+      kappa = (order+1)*(order+1);
+   }
+   if (sol_order < 0)
+   {
+      sol_order = 1;
+   }
+   args.PrintOptions(cout);
+
+   // 2. Construct the (serial) mesh and refine it if requested.
+   auto mesh = generate_mesh(ref_levels);
+
+   int dim = mesh.Dimension();
+
+   if (mesh.NURBSext)
+   {
+      mesh.SetCurvature(max(order, 1));
+   }
+
+   // 3. Define a finite element space on the mesh. Here we use discontinuous
+   //    finite elements of the specified order >= 0.
+   DG_FECollection fec(order, dim);
+   FiniteElementSpace fespace(&mesh, &fec);
+   cout << "Number of unknowns: " << fespace.GetVSize() << endl;
+
+   // 4. Set up the linear form b(.) which corresponds to the right-hand side of
+   //    the FEM linear system.
+   LinearForm b(&fespace);
+
+   Array<int> p1_attr_marker(mesh.attributes.Max());
+   p1_attr_marker = 0;
+   p1_attr_marker[0] = 1;
+
+   FunctionCoefficient p1_source([sol_order](const Vector &p)
+   {
+      const double x = p(0);
+      const double val = -(sol_order - 1)*sol_order*pow(x, sol_order-2);
+      return val;
+   });
+   b.AddDomainIntegrator(new DomainLFIntegrator(p1_source), p1_attr_marker);
+
+   Array<int> p2_attr_marker(mesh.attributes.Max());
+   p2_attr_marker = 0;
+   p2_attr_marker[1] = 1;
+
+   FunctionCoefficient p2_source([sol_order](const Vector &p)
+   {
+      const double x = p(0);
+      double val = -(sol_order - 1)*sol_order*pow(x - 2, sol_order-2);
+      if (sol_order % 2 == 0)
+      {
+         val *= -1.0;
+      }
+      return val;
+   });
+   b.AddDomainIntegrator(new DomainLFIntegrator(p2_source), p2_attr_marker);
+
+   ConstantCoefficient one(1.0);
+
+   Array<int> p1_bdr_attr_marker(mesh.bdr_attributes.Max());
+   p1_bdr_attr_marker = 0;
+   p1_bdr_attr_marker[0] = 1;
+
+   ConstantCoefficient left_bc_val(0.0);
+   b.AddBdrFaceIntegrator(
+      new DGDirichletLFIntegrator(left_bc_val, one, sigma, kappa),
+      p1_bdr_attr_marker);
+
+   Array<int> p2_bdr_attr_marker(mesh.bdr_attributes.Max());
+   p2_bdr_attr_marker = 0;
+   p2_bdr_attr_marker[1] = 1;
+
+   ConstantCoefficient right_bc_val(2.0 + jump);
+   b.AddBdrFaceIntegrator(
+      new DGDirichletLFIntegrator(right_bc_val, one, sigma, kappa),
+      p2_bdr_attr_marker);
+
+   b.Assemble();
+
+   // 5. Define the solution vector x as a finite element grid function
+   //    corresponding to fespace. Initialize x with initial guess of zero.
+   GridFunction x(&fespace);
+   x = 0.0;
+
+   // 6. 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(&fespace);
+   a.AddDomainIntegrator(new DiffusionIntegrator(one));
+
+   a.AddInteriorFaceIntegrator(new DGDiffusionIntegrator(one, sigma, kappa));
+   a.AddBdrFaceIntegrator(new DGDiffusionIntegrator(one, sigma, kappa),
+                          p1_bdr_attr_marker);
+   a.AddBdrFaceIntegrator(new DGDiffusionIntegrator(one, sigma, kappa),
+                          p2_bdr_attr_marker);
+   if (eta > 0)
+   {
+      a.AddInteriorFaceIntegrator(new DGDiffusionBR2Integrator(fespace, eta));
+      a.AddBdrFaceIntegrator(new DGDiffusionBR2Integrator(fespace, eta));
+   }
+
+   // 7. Negate the DG interface terms along the internal boundary so that the
+   //    only coupling between domains is from the chosen model (constant flux
+   //    in this case).
+   Array<int> internal_bdr_attr_marker(mesh.bdr_attributes.Max());
+   internal_bdr_attr_marker = 0;
+   internal_bdr_attr_marker[4] = 1;
+
+   ProductCoefficient neg_one(-1.0, one);
+   a.AddInternalBoundaryFaceIntegrator(new DGDiffusionIntegrator(neg_one, sigma,
+                                                                 kappa),
+                                       internal_bdr_attr_marker);
+   if (eta > 0)
+   {
+      a.AddInternalBoundaryFaceIntegrator(new DGDiffusionBR2Integrator(fespace,
+                                                                       neg_one, eta),
+                                          internal_bdr_attr_marker);
+   }
+
+   ConstantCoefficient mass_coeff(sol_order / jump);
+   a.AddInternalBoundaryFaceIntegrator(new InteriorMassIntegrator(mass_coeff),
+                                       internal_bdr_attr_marker);
+
+   a.Assemble();
+   a.Finalize();
+   const SparseMatrix &A = a.SpMat();
+
+#ifndef MFEM_USE_SUITESPARSE
+   // 8. 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 && !(jump < 0))
+   {
+      PCG(A, M, b, x, 1, 500, 1e-12, 0.0);
+   }
+   else
+   {
+      GMRES(A, M, b, x, 1, 500, 500, 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);
+#endif
+
+   // 9. 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);
+   ofstream sol_ofs("sol.gf");
+   sol_ofs.precision(8);
+   x.Save(sol_ofs);
+
+   // 10. 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;
+   }
+
+   return 0;
+}
+
+void InteriorMassIntegrator::AssembleFaceMatrix(
+   const FiniteElement &el1,
+   const FiniteElement &el2,
+   FaceElementTransformations &trans,
+   DenseMatrix &elmat)
+{
+   int ndof1 = el1.GetDof();
+   int ndof2 = el2.GetDof();
+   int ndof = ndof1 + ndof2;
+
+#ifdef MFEM_THREAD_SAFE
+   Vector shape1;
+   Vector shape2;
+   DenseMatrix elmat11;
+   DenseMatrix elmat12;
+   DenseMatrix elmat21;
+   DenseMatrix elmat22;
+#endif
+   shape1.SetSize(ndof1);
+   shape2.SetSize(ndof2);
+
+   elmat11.SetSize(ndof1);
+   elmat12.SetSize(ndof1, ndof2);
+   elmat21.SetSize(ndof2, ndof1);
+   elmat22.SetSize(ndof2);
+
+   const auto *ir = IntRule;
+   if (ir == NULL)
+   {
+      int order = 2 * max(el1.GetOrder(), el2.GetOrder());
+      ir = &IntRules.Get(trans.GetGeometryType(), order);
+   }
+
+   elmat.SetSize(ndof);
+   elmat = 0.0;
+   for (int i = 0; i < ir->GetNPoints(); i++)
+   {
+      const auto &ip = ir->IntPoint(i);
+
+      // Set the integration point in the face and the neighboring element
+      trans.SetAllIntPoints(&ip);
+
+      const double w = ip.weight * trans.Weight();
+
+      // Access the neighboring element's integration point
+      const auto &eip1 = trans.GetElement1IntPoint();
+      const auto &eip2 = trans.GetElement2IntPoint();
+
+      el1.CalcShape(eip1, shape1);
+      el2.CalcShape(eip2, shape2);
+
+      const double Q_val = Q.Eval(trans, ip);
+
+      elmat11 = 0.0;
+      AddMult_a_VVt(Q_val * w, shape1, elmat11);
+
+      elmat12 = 0.0;
+      AddMult_a_VWt(-Q_val * w, shape2, shape1, elmat12);
+
+      elmat21 = 0.0;
+      AddMult_a_VWt(-Q_val * w, shape1, shape2, elmat21);
+
+      elmat22 = 0.0;
+      AddMult_a_VVt(Q_val * w, shape2, elmat22);
+
+      for (int j = 0; j < ndof1; ++j)
+      {
+         for (int k = 0; k < ndof1; ++k)
+         {
+            elmat(j, k) += elmat11(j, k);
+         }
+      }
+
+      for (int j = 0; j < ndof1; ++j)
+      {
+         for (int k = 0; k < ndof2; ++k)
+         {
+            elmat(j, k + ndof1) += elmat12(j, k);
+            elmat(k + ndof1, j) += elmat21(k, j);
+         }
+      }
+
+      for (int j = 0; j < ndof2; ++j)
+      {
+         for (int k = 0; k < ndof2; ++k)
+         {
+            elmat(j + ndof1, k + ndof1) += elmat22(j, k);
+         }
+      }
+   }
+}
+
+Mesh generate_mesh(int ref, int internal_bdr_attr)
+{
+   int nxy = 4 * (ref+1);
+   auto mesh = Mesh::MakeCartesian2D(nxy, nxy, Element::TRIANGLE, true, 2.0, 1.0);
+   // auto mesh = Mesh::MakeCartesian2D(nxy, nxy, Element::QUADRILATERAL, true, 2.0, 1.0);
+
+   // assign element attributes to left and right sides
+   for (int i = 0; i < mesh.GetNE(); ++i)
+   {
+      auto *elem = mesh.GetElement(i);
+
+      Array<int> verts;
+      elem->GetVertices(verts);
+
+      bool left = true;
+      for (int j = 0; j < verts.Size(); ++j)
+      {
+         auto *vtx = mesh.GetVertex(verts[j]);
+         if (vtx[0] <= 1.0)
+         {
+            continue;
+         }
+         else
+         {
+            left = false;
+         }
+      }
+      if (left)
+      {
+         elem->SetAttribute(1);
+      }
+      else
+      {
+         elem->SetAttribute(2);
+      }
+   }
+
+   // assign boundary element attributes to left and right sides
+   for (int i = 0; i < mesh.GetNBE(); ++i)
+   {
+      auto *elem = mesh.GetBdrElement(i);
+
+      Array<int> verts;
+      elem->GetVertices(verts);
+
+      bool left = true;
+      bool right = true;
+      bool top = true;
+      bool bottom = true;
+      for (int j = 0; j < verts.Size(); ++j)
+      {
+         auto *vtx = mesh.GetVertex(verts[j]);
+         left = left && abs(vtx[0] - 0.0) < 1e-12;
+         right = right && abs(vtx[0] - 2.0) < 1e-12;
+         top = top && abs(vtx[1] - 1.0) < 1e-12;
+         bottom = bottom && abs(vtx[1] - 0.0) < 1e-12;
+      }
+      if (left)
+      {
+         elem->SetAttribute(1);
+      }
+      else if (right)
+      {
+         elem->SetAttribute(2);
+      }
+      else if (top)
+      {
+         elem->SetAttribute(3);
+      }
+      else if (bottom)
+      {
+         elem->SetAttribute(4);
+      }
+   }
+
+   // add internal boundary elements
+   for (int i = 0; i < mesh.GetNumFaces(); ++i)
+   {
+      int e1, e2;
+      mesh.GetFaceElements(i, &e1, &e2);
+      if (e1 >= 0 && e2 >= 0 && mesh.GetAttribute(e1) != mesh.GetAttribute(e2))
+      {
+         // This is the internal face between attributes.
+         auto *new_elem = mesh.GetFace(i)->Duplicate(&mesh);
+         new_elem->SetAttribute(internal_bdr_attr);
+         mesh.AddBdrElement(new_elem);
+      }
+   }
+
+   mesh.FinalizeTopology(); // Finalize to build relevant tables
+   mesh.Finalize();
+   mesh.SetAttributes();
+
+   return mesh;
+}

examples/makefile

@@ -23,7 +23,7 @@ MFEM_LIB_FILE = mfem_is_not_built
 
 SEQ_EXAMPLES = ex0 ex1 ex2 ex3 ex4 ex5 ex6 ex7 ex8 ex9 ex10 ex14 ex15 ex16 \
   ex17 ex18 ex19 ex20 ex21 ex22 ex23 ex24 ex25 ex26 ex27 ex28 ex29 ex30 \
-  ex31 ex33
+  ex31 ex33 ex34 ex35
 PAR_EXAMPLES = ex0p ex1p ex2p ex3p ex4p ex5p ex6p ex7p ex8p ex9p ex10p ex11p \
   ex12p ex13p ex14p ex15p ex16p ex17p ex18p ex19p ex20p ex21p ex22p ex24p \
   ex25p ex26p ex27p ex28p ex29p ex30p ex31p ex32p ex33p

fem/bilinearform.cpp

@@ -14,6 +14,7 @@
 #include "fem.hpp"
 #include "../general/device.hpp"
 #include <cmath>
+#include <cstddef>
 
 namespace mfem
 {
@@ -109,6 +110,9 @@ BilinearForm::BilinearForm (FiniteElementSpace * f, BilinearForm * bf, int ps)
    boundary_face_integs = bf->boundary_face_integs;
    boundary_face_integs_marker = bf->boundary_face_integs_marker;
 
+   internal_boundary_face_integs = bf->internal_boundary_face_integs;
+   internal_boundary_face_integs_marker = bf->internal_boundary_face_integs_marker;
+
    AllocMat();
 }
 
@@ -278,6 +282,22 @@ void BilinearForm::AddBdrFaceIntegrator(BilinearFormIntegrator *bfi,
    boundary_face_integs_marker.Append(&bdr_marker);
 }
 
+void BilinearForm::AddInternalBoundaryFaceIntegrator(BilinearFormIntegrator
+                                                     *bfi)
+{
+   internal_boundary_face_integs.Append(bfi);
+   // nullptr -> all attributes are active
+   internal_boundary_face_integs_marker.Append(nullptr);
+}
+
+void BilinearForm::AddInternalBoundaryFaceIntegrator(BilinearFormIntegrator
+                                                     *bfi,
+                                                     Array<int> &internal_bdr_attr_marker)
+{
+   internal_boundary_face_integs.Append(bfi);
+   internal_boundary_face_integs_marker.Append(&internal_bdr_attr_marker);
+}
+
 void BilinearForm::ComputeElementMatrix(int i, DenseMatrix &elmat)
 {
    if (element_matrices)
@@ -630,6 +650,59 @@ void BilinearForm::Assemble(int skip_zeros)
       }
    }
 
+   if (internal_boundary_face_integs.Size())
+   {
+      // Which internal boundary attributes need to be processed?
+      Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
+                                 mesh->bdr_attributes.Max() : 0);
+      bdr_attr_marker = 0;
+      for (int k = 0; k < internal_boundary_face_integs.Size(); k++)
+      {
+         if (internal_boundary_face_integs_marker[k] == NULL)
+         {
+            bdr_attr_marker = 1;
+            break;
+         }
+         auto &bdr_marker = *internal_boundary_face_integs_marker[k];
+         MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
+                     "invalid boundary marker for internal boundary face "
+                     "integrator #" << k << ", counting from zero");
+         for (int i = 0; i < bdr_attr_marker.Size(); i++)
+         {
+            bdr_attr_marker[i] |= bdr_marker[i];
+         }
+      }
+
+      Array<int> vdofs2;
+      for (int i = 0; i < mesh->GetNBE(); i++)
+      {
+         const int bdr_attr = mesh->GetBdrAttribute(i);
+         if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
+
+         auto *tr = mesh->GetInternalBdrFaceTransformations(i);
+         if (tr != nullptr)
+         {
+            fes->GetElementVDofs(tr->Elem1No, vdofs);
+            fes->GetElementVDofs(tr->Elem2No, vdofs2);
+            vdofs.Append(vdofs2);
+            const auto *fe1 = fes->GetFE(tr->Elem1No);
+            const auto *fe2 = fes->GetFE(tr->Elem2No);
+            for (int k = 0; k < internal_boundary_face_integs.Size(); k++)
+            {
+               if (internal_boundary_face_integs_marker[k] &&
+                   (*internal_boundary_face_integs_marker[k])[bdr_attr - 1] == 0)
+               {
+                  continue;
+               }
+
+               internal_boundary_face_integs[k]->AssembleFaceMatrix(
+                  *fe1, *fe2, *tr, elemmat);
+               mat->AddSubMatrix(vdofs, vdofs, elemmat, skip_zeros);
+            }
+         }
+      }
+   }
+
 #ifdef MFEM_USE_LEGACY_OPENMP
    if (free_element_matrices)
    {
@@ -1143,6 +1216,10 @@ BilinearForm::~BilinearForm()
       { delete interior_face_integs[k]; }
       for (k=0; k < boundary_face_integs.Size(); k++)
       { delete boundary_face_integs[k]; }
+      for (int i = 0; i < internal_boundary_face_integs.Size(); i++)
+      {
+         delete internal_boundary_face_integs[i];
+      }
    }
 
    delete ext;

fem/bilinearform.hpp

@@ -113,6 +113,10 @@ protected:
    Array<BilinearFormIntegrator*> boundary_face_integs;
    Array<Array<int>*> boundary_face_integs_marker; ///< Entries are not owned.
 
+   /// Set of internal boundary face integrators to be applied.
+   Array<BilinearFormIntegrator*> internal_boundary_face_integs;
+   Array<Array<int>*> internal_boundary_face_integs_marker; ///< Entries not owned.
+
    DenseMatrix elemmat;
    Array<int>  vdofs;
 
@@ -416,6 +420,18 @@ public:
    void AddBdrFaceIntegrator(BilinearFormIntegrator *bfi,
                              Array<int> &bdr_marker);
 
+   /// @brief Add new internal boundary face integrator. Assumes ownership of
+   /// @a bfi.
+   void AddInternalBoundaryFaceIntegrator(BilinearFormIntegrator *bfi);
+
+   /** @brief Add new internal boundary face integrator, restricted to the given
+       boundary attributes.
+
+       Assumes ownership of @a bfi. The array @a internal_bdr_attr_marker is
+       stored internally as a pointer to the given Array<int> object. */
+   void AddInternalBoundaryFaceIntegrator(BilinearFormIntegrator *bfi,
+                                          Array<int> &internal_bdr_attr_marker);
+
    /// Sets all sparse values of \f$ M \f$ and \f$ M_e \f$ to 'a'.
    void operator=(const double a)
    {

fem/linearform.cpp

@@ -36,6 +36,9 @@ LinearForm::LinearForm(FiniteElementSpace *f, LinearForm *lf)
 
    boundary_face_integs = lf->boundary_face_integs;
    boundary_face_integs_marker = lf->boundary_face_integs_marker;
+
+   internal_boundary_face_integs = lf->internal_boundary_face_integs;
+   internal_boundary_face_integs_marker = lf->internal_boundary_face_integs_marker;
 }
 
 void LinearForm::AddDomainIntegrator(LinearFormIntegrator *lfi)
@@ -101,6 +104,20 @@ void LinearForm::AddInteriorFaceIntegrator(LinearFormIntegrator *lfi)
    interior_face_integs.Append(lfi);
 }
 
+void LinearForm::AddInternalBoundaryFaceIntegrator(LinearFormIntegrator *lfi)
+{
+   internal_boundary_face_integs.Append(lfi);
+   // nullptr -> all attributes are active
+   internal_boundary_face_integs_marker.Append(nullptr);
+}
+
+void LinearForm::AddInternalBoundaryFaceIntegrator(LinearFormIntegrator *lfi,
+                                                   Array<int> &internal_bdr_attr_marker)
+{
+   internal_boundary_face_integs.Append(lfi);
+   internal_boundary_face_integs_marker.Append(&internal_bdr_attr_marker);
+}
+
 bool LinearForm::SupportsDevice() const
 {
    // return false for NURBS meshes, so we don’t convert it to non-NURBS
@@ -121,7 +138,10 @@ bool LinearForm::SupportsDevice() const
    if (!IntegratorsSupportDevice(domain_integs)) { return false; }
    if (!IntegratorsSupportDevice(boundary_integs)) { return false; }
    if (boundary_face_integs.Size() > 0 || interior_face_integs.Size() > 0 ||
-       domain_delta_integs.Size() > 0) { return false; }
+       domain_delta_integs.Size() > 0 || internal_boundary_face_integs.Size() > 0)
+   {
+      return false;
+   }
 
    if (boundary_integs.Size() > 0)
    {
@@ -339,6 +359,61 @@ void LinearForm::Assemble()
          }
       }
    }
+
+   if (internal_boundary_face_integs.Size())
+   {
+      auto *mesh = fes->GetMesh();
+
+      // Which internal boundary attributes need to be processed?
+      Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
+                                 mesh->bdr_attributes.Max() : 0);
+      bdr_attr_marker = 0;
+      for (int k = 0; k < internal_boundary_face_integs.Size(); k++)
+      {
+         if (internal_boundary_face_integs_marker[k] == NULL)
+         {
+            bdr_attr_marker = 1;
+            break;
+         }
+         auto &bdr_marker = *internal_boundary_face_integs_marker[k];
+         MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
+                     "invalid boundary marker for internal boundary face "
+                     "integrator #" << k << ", counting from zero");
+         for (int i = 0; i < bdr_attr_marker.Size(); i++)
+         {
+            bdr_attr_marker[i] |= bdr_marker[i];
+         }
+      }
+
+      Array<int> vdofs2;
+      for (int i = 0; i < mesh->GetNBE(); i++)
+      {
+         const int bdr_attr = mesh->GetBdrAttribute(i);
+         if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
+
+         auto *tr = mesh->GetInternalBdrFaceTransformations(i);
+         if (tr != nullptr)
+         {
+            fes->GetElementVDofs(tr->Elem1No, vdofs);
+            fes->GetElementVDofs(tr->Elem2No, vdofs2);
+            vdofs.Append(vdofs2);
+            const auto *fe1 = fes->GetFE(tr->Elem1No);
+            const auto *fe2 = fes->GetFE(tr->Elem2No);
+            for (int k = 0; k < internal_boundary_face_integs.Size(); k++)
+            {
+               if (internal_boundary_face_integs_marker[k] &&
+                   (*internal_boundary_face_integs_marker[k])[bdr_attr - 1] == 0)
+               {
+                  continue;
+               }
+
+               internal_boundary_face_integs[k]->AssembleRHSElementVect(
+                  *fe1, *fe2, *tr, elemvect);
+               AddElementVector(vdofs, elemvect);
+            }
+         }
+      }
+   }
 }
 
 void LinearForm::Update()
@@ -429,6 +504,10 @@ LinearForm::~LinearForm()
       { delete boundary_face_integs[k]; }
       for (k=0; k < interior_face_integs.Size(); k++)
       { delete interior_face_integs[k]; }
+      for (int i = 0; i < internal_boundary_face_integs.Size(); i++)
+      {
+         delete internal_boundary_face_integs[i];
+      }
    }
 
    delete ext;

fem/linearform.hpp

@@ -65,6 +65,10 @@ protected:
    /// Set of Internal Face Integrators to be applied.
    Array<LinearFormIntegrator*> interior_face_integs;
 
+   /// Set of internal boundary face integrators to be applied.
+   Array<LinearFormIntegrator*> internal_boundary_face_integs;
+   Array<Array<int>*> internal_boundary_face_integs_marker; ///< Entries not owned.
+
    /// The element ids where the centers of the delta functions lie
    Array<int> domain_delta_integs_elem_id;
 
@@ -162,6 +166,18 @@ public:
    /// Adds new Interior Face Integrator. Assumes ownership of @a lfi.
    void AddInteriorFaceIntegrator(LinearFormIntegrator *lfi);
 
+   /// @brief Add new internal boundary face integrator. Assumes ownership of
+   /// @a lfi.
+   void AddInternalBoundaryFaceIntegrator(LinearFormIntegrator *lfi);
+
+   /** @brief Add new internal boundary face integrator, restricted to the given
+       boundary attributes.
+
+       Assumes ownership of @a lfi. The array @a internal_bdr_attr_marker is
+       stored internally as a pointer to the given Array<int> object. */
+   void AddInternalBoundaryFaceIntegrator(LinearFormIntegrator *lfi,
+                                          Array<int> &internal_bdr_attr_marker);
+
    /** @brief Access all integrators added with AddDomainIntegrator() which are
        not DeltaLFIntegrator%s or they are DeltaLFIntegrator%s with non-delta
        coefficients. */

fem/nonlinearform.cpp

@@ -89,6 +89,8 @@ double NonlinearForm::GetGridFunctionEnergy(const Vector &x) const
    {
       MFEM_VERIFY(!fnfi.Size(), "Interior faces terms not yet implemented!");
       MFEM_VERIFY(!bfnfi.Size(), "Boundary face terms not yet implemented!");
+      MFEM_VERIFY(!internal_boundary_face_integs.Size(),
+                  "Internal boundary face terms not yet implemented!");
       return ext->GetGridFunctionEnergy(x);
    }
 
@@ -125,6 +127,11 @@ double NonlinearForm::GetGridFunctionEnergy(const Vector &x) const
       MFEM_ABORT("TODO: add energy contribution from boundary face terms");
    }
 
+   if (internal_boundary_face_integs.Size())
+   {
+      MFEM_ABORT("TODO: add energy contribution from internal boundary face terms");
+   }
+
    return energy;
 }
 
@@ -274,6 +281,62 @@ void NonlinearForm::Mult(const Vector &x, Vector &y) const
       }
    }
 
+   if (internal_boundary_face_integs.Size())
+   {
+      // Which internal boundary attributes need to be processed?
+      Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
+                                 mesh->bdr_attributes.Max() : 0);
+      bdr_attr_marker = 0;
+      for (int k = 0; k < internal_boundary_face_integs.Size(); k++)
+      {
+         if (internal_boundary_face_integs_marker[k] == NULL)
+         {
+            bdr_attr_marker = 1;
+            break;
+         }
+         auto &bdr_marker = *internal_boundary_face_integs_marker[k];
+         MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
+                     "invalid boundary marker for internal boundary face "
+                     "integrator #" << k << ", counting from zero");
+         for (int i = 0; i < bdr_attr_marker.Size(); i++)
+         {
+            bdr_attr_marker[i] |= bdr_marker[i];
+         }
+      }
+
+      Array<int> vdofs2;
+      for (int i = 0; i < mesh->GetNBE(); i++)
+      {
+         const int bdr_attr = mesh->GetBdrAttribute(i);
+         if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
+
+         auto *tr = mesh->GetInternalBdrFaceTransformations(i);
+         if (tr != nullptr)
+         {
+            fes->GetElementVDofs(tr->Elem1No, vdofs);
+            fes->GetElementVDofs(tr->Elem2No, vdofs2);
+            vdofs.Append(vdofs2);
+
+            px.GetSubVector(vdofs, el_x);
+
+            const auto *fe1 = fes->GetFE(tr->Elem1No);
+            const auto *fe2 = fes->GetFE(tr->Elem2No);
+            for (int k = 0; k < internal_boundary_face_integs.Size(); k++)
+            {
+               if (internal_boundary_face_integs_marker[k] &&
+                   (*internal_boundary_face_integs_marker[k])[bdr_attr - 1] == 0)
+               {
+                  continue;
+               }
+
+               internal_boundary_face_integs[k]->AssembleFaceVector(
+                  *fe1, *fe2, *tr, el_x, el_y);
+               py.AddElementVector(vdofs, el_y);
+            }
+         }
+      }
+   }
+
    if (Serial())
    {
       if (cP) { cP->MultTranspose(py, y); }
@@ -422,6 +485,62 @@ Operator &NonlinearForm::GetGradient(const Vector &x) const
       }
    }
 
+   if (internal_boundary_face_integs.Size())
+   {
+      // Which internal boundary attributes need to be processed?
+      Array<int> bdr_attr_marker(mesh->bdr_attributes.Size() ?
+                                 mesh->bdr_attributes.Max() : 0);
+      bdr_attr_marker = 0;
+      for (int k = 0; k < internal_boundary_face_integs.Size(); k++)
+      {
+         if (internal_boundary_face_integs_marker[k] == NULL)
+         {
+            bdr_attr_marker = 1;
+            break;
+         }
+         auto &bdr_marker = *internal_boundary_face_integs_marker[k];
+         MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
+                     "invalid boundary marker for internal boundary face "
+                     "integrator #" << k << ", counting from zero");
+         for (int i = 0; i < bdr_attr_marker.Size(); i++)
+         {
+            bdr_attr_marker[i] |= bdr_marker[i];
+         }
+      }
+
+      Array<int> vdofs2;
+      for (int i = 0; i < mesh->GetNBE(); i++)
+      {
+         const int bdr_attr = mesh->GetBdrAttribute(i);
+         if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
+
+         auto *tr = mesh->GetInternalBdrFaceTransformations(i);
+         if (tr != nullptr)
+         {
+            fes->GetElementVDofs(tr->Elem1No, vdofs);
+            fes->GetElementVDofs(tr->Elem2No, vdofs2);
+            vdofs.Append(vdofs2);
+
+            px.GetSubVector(vdofs, el_x);
+
+            const auto *fe1 = fes->GetFE(tr->Elem1No);
+            const auto *fe2 = fes->GetFE(tr->Elem2No);
+            for (int k = 0; k < internal_boundary_face_integs.Size(); k++)
+            {
+               if (internal_boundary_face_integs_marker[k] &&
+                   (*internal_boundary_face_integs_marker[k])[bdr_attr - 1] == 0)
+               {
+                  continue;
+               }
+
+               internal_boundary_face_integs[k]->AssembleFaceGrad(
+                  *fe1, *fe2, *tr, el_x, elmat);
+               Grad->AddSubMatrix(vdofs, vdofs, elmat, skip_zeros);
+            }
+         }
+      }
+   }
+
    if (!Grad->Finalized())
    {
       Grad->Finalize(skip_zeros);
@@ -474,6 +593,10 @@ NonlinearForm::~NonlinearForm()
    for (int i = 0; i <  dnfi.Size(); i++) { delete  dnfi[i]; }
    for (int i = 0; i <  fnfi.Size(); i++) { delete  fnfi[i]; }
    for (int i = 0; i < bfnfi.Size(); i++) { delete bfnfi[i]; }
+   for (int i = 0; i < internal_boundary_face_integs.Size(); i++)
+   {
+      delete internal_boundary_face_integs[i];
+   }
    delete ext;
 }
 

fem/nonlinearform.hpp

@@ -45,6 +45,10 @@ protected:
    Array<NonlinearFormIntegrator*> bfnfi; // owned
    Array<Array<int>*>              bfnfi_marker; // not owned
 
+   /// Set of internal boundary face integrators to be applied.
+   Array<NonlinearFormIntegrator*> internal_boundary_face_integs;
+   Array<Array<int>*> internal_boundary_face_integs_marker; ///< Entries not owned.
+
    mutable SparseMatrix *Grad, *cGrad; // owned
    /// Gradient Operator when not assembled as a matrix.
    mutable OperatorHandle hGrad; // has internal ownership flag
@@ -138,6 +142,33 @@ public:
    const Array<NonlinearFormIntegrator*> &GetBdrFaceIntegrators() const
    { return bfnfi; }
 
+   /// @brief Add new internal boundary face integrator. Assumes ownership of
+   /// @a nfi.
+   void AddInternalBoundaryFaceIntegrator(NonlinearFormIntegrator *nfi)
+   {
+      internal_boundary_face_integs.Append(nfi);
+      // nullptr -> all attributes are active
+      internal_boundary_face_integs_marker.Append(nullptr);
+   }
+
+   /** @brief Add new internal boundary face integrator, restricted to the given
+       boundary attributes.
+
+       Assumes ownership of @a nfi. The array @a internal_bdr_attr_marker is
+       stored internally as a pointer to the given Array<int> object. */
+   void AddInternalBoundaryFaceIntegrator(NonlinearFormIntegrator *nfi,
+                                          Array<int> &internal_bdr_attr_marker)
+   {
+      internal_boundary_face_integs.Append(nfi);
+      internal_boundary_face_integs_marker.Append(&internal_bdr_attr_marker);
+   }
+
+   /** @brief Access all boundary face integrators added with
+       AddBdrFaceIntegrator(). */
+   const Array<NonlinearFormIntegrator*> &GetInternalBoundaryFaceIntegrators()
+   const
+   { return internal_boundary_face_integs; }
+
    /// Specify essential boundary conditions.
    /** This method calls FiniteElementSpace::GetEssentialTrueDofs() and stores
        the result internally for use by other methods. If the @a rhs pointer is

fem/pbilinearform.cpp

@@ -264,6 +264,11 @@ void ParBilinearForm::AssembleSharedFaces(int skip_zeros)
 
 void ParBilinearForm::Assemble(int skip_zeros)
 {
+   if (internal_boundary_face_integs.Size())
+   {
+      MFEM_ABORT("TODO: assemble contributions from shared internal boundary terms");
+   }
+
    if (interior_face_integs.Size())
    {
       pfes->ExchangeFaceNbrData();

fem/pnonlinearform.cpp

@@ -38,6 +38,11 @@ double ParNonlinearForm::GetParGridFunctionEnergy(const Vector &x) const
       MFEM_ABORT("TODO: add energy contribution from shared faces");
    }
 
+   if (internal_boundary_face_integs.Size())
+   {
+      MFEM_ABORT("TODO: add energy contributions from shared internal boundary terms");
+   }
+
    MPI_Allreduce(&loc_energy, &glob_energy, 1, MPI_DOUBLE, MPI_SUM,
                  ParFESpace()->GetComm());
 
@@ -46,6 +51,11 @@ double ParNonlinearForm::GetParGridFunctionEnergy(const Vector &x) const
 
 void ParNonlinearForm::Mult(const Vector &x, Vector &y) const
 {
+   if (internal_boundary_face_integs.Size() != 0)
+   {
+      MFEM_ABORT("TODO: assemble contributions from shared internal boundary terms");
+   }
+
    NonlinearForm::Mult(x, y); // x --(P)--> aux1 --(A_local)--> aux2
 
    if (fnfi.Size())
@@ -116,6 +126,11 @@ Operator &ParNonlinearForm::GetGradient(const Vector &x) const
 
    OperatorHandle dA(pGrad.Type()), Ph(pGrad.Type());
 
+   if (internal_boundary_face_integs.Size() != 0)
+   {
+      MFEM_ABORT("TODO: assemble contributions from shared internal boundary terms");
+   }
+
    if (fnfi.Size() == 0)
    {
       dA.MakeSquareBlockDiag(pfes->GetComm(), pfes->GlobalVSize(),

mesh/mesh.cpp

@@ -1117,6 +1117,25 @@ FaceElementTransformations *Mesh::GetBdrFaceTransformations(int BdrElemNo)
    return tr;
 }
 
+FaceElementTransformations *Mesh::GetInternalBdrFaceTransformations(
+   int IntBdrElemNo)
+{
+   int fn = GetBdrFace(IntBdrElemNo);
+
+   // Check if the face is not interior
+   if (!FaceIsTrueInterior(fn))
+   {
+      return nullptr;
+   }
+
+   auto *tr = GetFaceElementTransformations(fn, 31);
+   tr->Attribute = boundary[IntBdrElemNo]->GetAttribute();
+   tr->ElementNo = IntBdrElemNo;
+   tr->ElementType = ElementTransformation::BDR_FACE;
+   tr->mesh = this;
+   return tr;
+}
+
 int Mesh::GetBdrFace(int BdrElemNo) const
 {
    int fn;

mesh/mesh.hpp

@@ -1389,6 +1389,10 @@ public:
    /// @note The returned object should NOT be deleted by the caller.
    FaceElementTransformations *GetBdrFaceTransformations (int BdrElemNo);
 
+   /// Builds the transformation defining the given internal face.
+   /// @note The returned object should NOT be deleted by the caller.
+   FaceElementTransformations *GetInternalBdrFaceTransformations(int IntBdrElemNo);
+
    /// Return the local face index for the given boundary face.
    int GetBdrFace(int BdrElemNo) const;