Understanding Boundary integration for L2 Spaces

[AlejandroMunozManterola]

Good day to everyone,

I have an inquiry on the inner workings of surface integration for L2 spaces. While running some examples, I have found myself finding that the assembly of LinearForms didn't do what I expected in terms of surface integration.

To try to explain the issue, I have gone with a simple example:

Lets say we have a basic 2D mesh, where a 1.0 x 1.0 square is split in four triangles, we ignore three of the boundaries and mark the right one with a specific boundary tag, which is the one we'll be working with.

Now, in the most basic L2 Space, order 1 and dimension 2, for the designated element, the dof distribution would theoretically look like this (please ignore if the numbering of the dof does not equal how MFEM would tag them, it is here for referencing purposes). The red face and dof are the ones that would be assigned to the marked boundary, while the blue ones are the rest of the sub-elements of the element.

So, my goal is to perform a surface integral on the red marked face. As it can be seen in #2320, we cannot use AddBoundaryIntegrator, as this method does not support L2 spaces due to the nature of the spaces themselves where the dof are not on the faces but infinitesimally inside the element and for a single face, the two elements that share it have their own 'individual' faces, that leaves AddBdrFaceIntegrator as our only method to perform these surface integrals (ignoring that we can just modify the code to force-implement our will, lets try to work within MFEM's support if possible).

So I design an integrator that does what I want to do, being explicit, I am performing the following calculation:

https://latex.codecogs.com/svg.image?\iint_{S}n_{x}E_{y}\;e^{jkr'cos{\psi}}

Where nx is the inverted normal on the face (for meshing purposes, I merely need to invert the normal), Ey is a GridFunction defined on the same L2 FES we're using, and the exponential term is merely a FunctionCoefficient we pass to our Integrator as an argument.

Thus, I'd expect my BdrFace integrator to look something like this:

https://github.com/OpenSEMBA/dgtd/blob/d67851d696d742f83b147a3e0820a6712e230e4a/src/mfemExtension/LinearIntegrators.cpp#L153-L186

I understand, the lines about the normal can look awkward, but I am doubling this integrator for both 2D and 3D use, so please ignore it unless its implemetation is incorrect. Likewise, dir_ is an argument passed to the LinearForm when adding this integrator to specify the direction we want in this calculation, and c_ is the previously defined FunctionCoefficient. (If needed, I can post the code here, but I doubt the implementation of the function is wrong)

I am in the understanding that after performing this calculation, I can merely perform the inner product with my GridFunction Ey and obtain the desired outcome.

Now, the issue comes when I go back to the LinearForm::Assemble() method, specifically the if statement in which we delve into the BdrFace integrators.

mfem/fem/linearform.cpp

Lines 271 to 319 in 37e03ab

	if (boundary_face_integs.Size())
	{
	FaceElementTransformations *tr;
	Mesh *mesh = fes->GetMesh();
	// Which 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 < boundary_face_integs.Size(); k++)
	{
	if (boundary_face_integs_marker[k] == NULL)
	{
	bdr_attr_marker = 1;
	break;
	}
	Array<int> &bdr_marker = *boundary_face_integs_marker[k];
	MFEM_ASSERT(bdr_marker.Size() == bdr_attr_marker.Size(),
	"invalid boundary marker for boundary face integrator #"
	<< k << ", counting from zero");
	for (int i = 0; i < bdr_attr_marker.Size(); i++)
	{
	bdr_attr_marker[i] |= bdr_marker[i];
	}
	}
	for (int i = 0; i < mesh->GetNBE(); i++)
	{
	const int bdr_attr = mesh->GetBdrAttribute(i);
	if (bdr_attr_marker[bdr_attr-1] == 0) { continue; }
	tr = mesh->GetBdrFaceTransformations(i);
	if (tr != NULL)
	{
	fes -> GetElementVDofs (tr -> Elem1No, vdofs);
	for (int k = 0; k < boundary_face_integs.Size(); k++)
	{
	if (boundary_face_integs_marker[k] &&
	(*boundary_face_integs_marker[k])[bdr_attr-1] == 0)
	{ continue; }
	boundary_face_integs[k]->
	AssembleRHSElementVect(*fes->GetFE(tr->Elem1No),
	*tr, elemvect);
	AddElementVector (vdofs, elemvect);
	}
	}
	}
	}

Indeed, the ElementTransformation we get is a FaceElementTransformation which has information about the face itself, and both elements on each side. The core of my doubts, and this Issue, is the vdof on which we will add our calculated integral are ALL of the dof that belong to the element, and not just those that belong to the face (like it is in the Boundary Integrator lines).

If we'd used AddBoundaryIntegrator:

mfem/fem/linearform.cpp

Line 253 in 37e03ab

doftrans = fes -> GetBdrElementVDofs (i, vdofs);

If we use AddBdrFaceIntegrator:

mfem/fem/linearform.cpp

Line 305 in 37e03ab

fes -> GetElementVDofs (tr -> Elem1No, vdofs);

Furthermore, when we define the sizes for the shape_ vector and we load these values through CalcShape, we're doing so on the IntegrationPoint defined by 'const IntegrationPoint& eip = Tr.GetElement1IntPoint();', so I find myself at a loss. In this surface integral calculation method, we are taking into account all of the weights (and later on, values) that exist inside the element, but I only require those that are part of the face itself and not the rest of the degrees of freedom that are part of the element.

Am I overthinking this and this is how it simply works for L2 spaces? Do I need to tweak some of this code to make it work, and if so, can you please give me some hints on how to proceed?

As another consideration, the elements I am trying to integrate here do not share faces with each other, only points as it can be seen in #4114, as if they formed a 'spiked crown'. Considering I am not going to temporally evolve these results, but just perform the surface integral, could I project my solutions on a H1 space and then use AddBoundaryIntegrator instead?

I am sorry if the explanation is too confusing, I don't want to make the issue excessively long. If it'd help though, I'll gladly try to mellow down the problem in a more step-by-step manner.

Thank you very much for your help.

Sincerely,
Alejandro M.

[najlkin]

Hi, what you observe is a totally correct behavior. The key is to understand how L2 elements work. Unlike continuous elements, L2 elements are not precisely defined at boundaries of elements and the notion of trace is used instead, which is not defined uniquely and different common choices exist. Practically, this means "open" bases are used for them typically, like the default Gauss-Legendre, which do not have DOFs at the boundaries. This implies that all DOFs essentially contribute to interpolation of the values (even though we may expect the closest contribute the most). This is the reason why the integration is performed on the whole element.

Projection to H1 might be an option, but keep in mind it is a global procedure, so it is a kind of integral fitting. This might cause artificial dislocation of some values, which might or might not be welcomed depending on the application (i.e., do you expect continuous or discontinuous profiles?). I would try to think differently maybe, I do not know your application exactly, but cannot you reformulate the problem in the weak sense? Like instead of calculation of strong divergence/curl, calculate them weakly? That is how L2 elements are typically used in EM codes 😉