Confusion regarding L2_FECollection, Boundary Integrators

[samuelpmishLLNL]

I'm trying to use mfem to build some tools for doing surface integrals in serac/smith, and running into several confusing issues related to using an L2_FECollection with different boundary integrators:

1. A finite element space made from an L2_FECollection doesn't appear to define boundary or face elements:

   
   int p = 2;
   int dim = 2;
   {
   auto fec = mfem::H1_FECollection(p, dim);
   mfem::FiniteElementSpace fespace(some_2D_mesh, &fec);
   auto foo = fespace.GetBE(0); // works
   auto bar = fespace.GetFaceElement(0); // works
   }

{ auto fec = mfem::L2_FECollection(p, dim); mfem::FiniteElementSpace fespace(some_2D_mesh, &fec); auto foo = fespace.GetBE(0); // null pointer (?) auto bar = fespace.GetFaceElement(0); // null pointer (?) }

Is this expected? Conceptually, whether or not an element exists seems like a property of the mesh (which is the same in both cases).

2. It is not clear what the difference is between `AddBoundaryIntegrator` and `AddBdrFaceIntegrator`. The doxygen documentation has only tautological descriptions of these functions:
  - [`LinearForm::AddBdrFaceIntegrator`](https://mfem.github.io/doxygen/html/classmfem_1_1LinearForm.html#af084ba37899ab6b938aae58a614ee36e): "Adds new Boundary Face Integrator"
  - [`LinearForm::AddBoundaryIntegrator`](https://mfem.github.io/doxygen/html/classmfem_1_1LinearForm.html#a287bc86983b37c829a3312aab0dcbd53): "Adds new Boundary Integrator"

3. It seems like boundary integrators are expected to be used with exactly one of `LinearForm::AddBdrFaceIntegrator` or `LinearForm::AddBoundaryIntegrator`, but there is nothing in place to convey or enforce this, resulting in mysterious segfaults later on. e.g.
```cpp
auto fec = mfem::L2_FECollection(p, dim);
mfem::ParFiniteElementSpace fespace(&mesh, &fec);

// this example below is well-formed and compiles without warning or error
{
  mfem::LinearForm f(&fespace);
  f.AddBdrFaceIntegrator(new mfem::BoundaryLFIntegrator(some_coef));
  f.Assemble(); // works!
}

// the example below is apparently ill-formed but compiles without warning or error
{
  mfem::LinearForm f(&fespace);
  f.AddBoundaryIntegrator(new mfem::BoundaryLFIntegrator(some_coef)); // this is not expected to work (?) but does not warn or error
  f.Assemble(); // segfault as a result of issues (1) and (2) above
}

It is very confusing to have a function named AddBoundaryIntegrator that does not support something called Boundary***Integrator (and yet, still compiles without issue).

Similar problems exist with BilinearForm::AddBoundaryIntegrator vs. BilinearForm::AddBdrFaceIntegrator. If BoundaryIntegrators and BdrFaceIntegrators are fundamentally different things, making them different types would completely prevent this issue, by failing to compile. Alternatively, Add***Integrator could verify that its argument is supported, and emit an error otherwise.

4. The documentation on the supported kinds of integrators for LinearForm appears to be incorrect. It says that BoundaryLFNormalIntegrator is supported on L2 spaces, but it doesn't appear to implement AssembleRHSElementVect (FiniteElement &, FaceElementTransformations &, Vector &), making it unusable by both LinearForm::AddBdrFaceIntegrator and LinearForm::AddBoundaryIntegrator:

auto fec = mfem::L2_FECollection(p, dim);
mfem::ParFiniteElementSpace fespace(&mesh, &fec);

{
  mfem::LinearForm f(&fespace);
  f.AddBdrFaceIntegrator(new mfem::BoundaryNormalLFIntegrator(some_vector_coef));
  f.Assemble(); // runtime error: LinearFormIntegrator::AssembleRHSElementVect(...) unimplemented
}

{
  mfem::LinearForm f(&fespace);
  f.AddBoundaryIntegrator(new mfem::BoundaryNormalLFIntegrator(some_vector_coef));
  f.Assemble(); // segfault, as in (3) above
}

[pazner]

I think you bring up some valid points. I don't have all the answers to your questions, but I can try to clarify a couple of the points.

1. A finite element space made from an L2_FECollection doesn't appear to define boundary or face elements: […] Is this expected? Conceptually, whether or not an element exists seems like a property of the mesh (which is the same in both cases).

I believe this is an intentional design decision. An interior face in an L2 space is not a well-defined finite element, because the trace of the DG field on that face is double-valued.

Now, on the domain boundary, the trace actually is single-valued, but I don't believe MFEM supports working with the boundary as finite elements, partly because (depending on the DG basis), there are not necessarily any "boundary DOFs" in the element (all the DOFs may be interior, as in the case of Gauss-Legendre basis).

2. It is not clear what the difference is between AddBoundaryIntegrator and AddBdrFaceIntegrator. The doxygen documentation has only tautological descriptions of these functions:

My understanding is that AddBoundaryIntegrator is intended for working with finite element spaces that support boundary elements (i.e. that have well-defined DOFs on the boundary). AddBdrFaceIntegrator is supposed to work with so-called "face integrators", just like AddInteriorFaceIntegrator. The main purpose of these face integrators is to implement DG integrators that support double-valued traces. As you mention further down in your comment, the "face integrators" (both Bdr and Interior) are the ones that work with FaceElementTransformations, while the "boundary integrators" work with an ElementTransformation that represents the boundary element.

I completely agree that the documentation could be improved.

3. It seems like boundary integrators are expected to be used with exactly one of LinearForm::AddBdrFaceIntegrator or LinearForm::AddBoundaryIntegrator, but there is nothing in place to convey or enforce this, resulting in mysterious segfaults later on. […]

   It is very confusing to have a function named AddBoundaryIntegrator that does not support something called Boundary***Integrator (and yet, still compiles without issue).

   Similar problems exist with BilinearForm::AddBoundaryIntegrator vs. BilinearForm::AddBdrFaceIntegrator. If BoundaryIntegrators and BdrFaceIntegrators are fundamentally different things, making them different types would completely prevent this issue, by failing to compile. Alternatively, Add***Integrator could verify that its argument is supported, and emit an error otherwise.

I largely agree. Making them different types would be nice, so we could get compiler errors vs. runtime errors. It is hard to anticipate how feasible this refactor would be, but it would be good to think about.

Other (less comprehensive) solutions could include better documentation or runtime verification (with a better error message).

Note that in your example, the problem is not that BoundaryLFIntegrator does not work with AddBoundaryIntegrator, but rather that AddBoundaryIntegrator does not work with L2 spaces. It would definitely be confusing if AddBoundaryIntegrator did not work with an integrator named Boundary***Integrator, and we should avoid creating integrators named like that, but that is not what is happening in this case.

4. The documentation on the supported kinds of integrators for LinearForm appears to be incorrect. […]

This should be fixed, I agree.

[v-dobrev]

I'll add a "todo" label with the intent to add better documentation for the various Add*Integrator methods in classes like LinearForm, BilinearForm, MixedBilinearForm, etc.

[stale[bot]]

:warning: This issue or PR has been automatically marked as stale because it has not had any activity in the last month. If no activity occurs in the next week, it will be automatically closed. Thank you for your contributions.

[stale[bot]]

:warning: This issue or PR has been automatically marked as stale because it has not had any activity in the last month. If no activity occurs in the next week, it will be automatically closed. Thank you for your contributions.