Fix exterior facet integrals from meshtags

[jorgensd]

#3900 introduced a bug, in the sense that ds(tag) now doesn't discriminate between what is an exterior facet and what is an interior facet.

This PR corrects this behavior by adding a set intersection with the exterior facets while computing the integration domains.

First reported at: https://fenicsproject.discourse.group/t/inconsistent-behaviors-of-the-ds-measure-between-dolfinx-v0-9-0-and-v0-10-0/19352

[schnellerhase]

Is this really a bug? If we assemble an external integral over an interior interface this is exactly the behaviour one would want. Whether this a right thing to do is another question. The integration domain should be dictated by the entity selection, so relying on dolfinx.mesh.locate_entities_boundary for the identification of those and not the silent dropping in the assembler?

[jorgensd]

Is this really a bug? If we assemble an external integral over an interior interface this is exactly the behaviour one would want. Whether this a right thing to do is another question. The integration domain should be dictated by the entity selection, so relying on dolfinx.mesh.locate_entities_boundary for the identification of those and not the silent dropping in the assembler?

The standard behavior for FEniCS for as long as I can tell is that ds means exterior facet.
changing that silently will cause many unintended bugs. It is also what we call it internally.

as we have discussed before, having an arbitrary one sided restriction as default behavior is not a good idea. (As the kernels are not well defined for many inputs, such as gradients, normals or discontinuous fields).

It should be supported through manual packing of integration data, for expert users, but not be the default for our high level interface.

[schnellerhase]

Hmm, but the "simple" interfaces here work fine:

value = dolfinx.fem.assemble_scalar(dolfinx.fem.form(dolfinx.fem.Constant(domain, 1.0)*ufl.ds))
print(value) # 4

value = dolfinx.fem.assemble_scalar(dolfinx.fem.form(dolfinx.fem.Constant(domain, 1.0)*ufl.dx + dolfinx.fem.Constant(domain, 1.0)*ufl.ds))
print(value) # 5

it only breaks down when one passes explicitly marked internal facets.

[jorgensd]

Hmm, but the "simple" interfaces here work fine:

value = dolfinx.fem.assemble_scalar(dolfinx.fem.form(dolfinx.fem.Constant(domain, 1.0)*ufl.ds))
print(value) # 4

value = dolfinx.fem.assemble_scalar(dolfinx.fem.form(dolfinx.fem.Constant(domain, 1.0)*ufl.dx + dolfinx.fem.Constant(domain, 1.0)*ufl.ds))
print(value) # 5

it only breaks down when one passes explicitly marked internal facets.

As I stated, it fails whenever you have a meshtags object as part of ds that tags interior (and exterior) facets. This is because compute integration data at the moment picks the first cell connected to each facet, irregardless of it being interior or exterior.
this can happen in many cases when one reads in meshes, where an interior surface has the same tag as an exterior surface.

Your example also illustrates the inconsistency of ds, as one would expect ds to integrate more than ds(tag), but if the subdomain data has all facets in the grid tagged with the value tag it would integrate over all interior and exterior facets.

The new interface additionally includes ghost entries from the mesh tags on all processes as well, which would lead to inconsistent parallel behavior.