Add one-sided integration example as example of custom integration

[jorgensd]

# Computing a onesided integral over an interior facet with consistent orientation
# Enabled by: https://github.com/FEniCS/dolfinx/pull/2269
# Copyright 2023 Jørgen S. Dokken
# SPDX: MIT

import dolfinx
import numpy as np
import ufl
from mpi4py import MPI

mesh = dolfinx.mesh.create_unit_square(
    MPI.COMM_WORLD, 10, 10, ghost_mode=dolfinx.mesh.GhostMode.shared_facet)

# Create connectivties required for defining integration entities
tdim = mesh.topology.dim
fdim = tdim - 1
mesh.topology.create_entities(fdim)
mesh.topology.create_connectivity(fdim, tdim)
mesh.topology.create_connectivity(tdim, fdim)


# Get number of cells on process
cell_map = mesh.topology.index_map(tdim)
num_cells = cell_map.size_local + cell_map.num_ghosts

# Create markers for each size of the interface
cell_values = np.ones(num_cells, dtype=np.int32)
cell_values[dolfinx.mesh.locate_entities(
    mesh, tdim, lambda x: x[0] <= 0.5+1e-13)] = 2
ct = dolfinx.mesh.meshtags(mesh, tdim, np.arange(
    num_cells, dtype=np.int32), cell_values)

facet_map = mesh.topology.index_map(fdim)
num_facets = facet_map.size_local + facet_map.num_ghosts

# Create facet markers
facet_values = np.ones(num_facets, dtype=np.int32)
facet_values[dolfinx.mesh.locate_entities(
    mesh, fdim, lambda x: np.isclose(x[0], 0.5))] = 2
ft = dolfinx.mesh.meshtags(mesh, fdim, np.arange(
    num_facets, dtype=np.int32), facet_values)

# Give a set of facets marked with a value (in this case 2), get a consistent orientation for an interior integral
facets_to_integrate = ft.find(2)

f_to_c = mesh.topology.connectivity(fdim, tdim)
c_to_f = mesh.topology.connectivity(tdim, fdim)
# Compute integration entities for a single facet of a cell.
# Each facet is represented as a tuple (cell_index, local_facet_index), where cell_index is local to process
# local_facet_index is the local indexing of a facet for a given cell
integration_entities = []
for i, facet in enumerate(facets_to_integrate):
    # Only loop over facets owned by the process to avoid duplicate integration
    if facet >= facet_map.size_local:
        continue
    # Find cells connected to facet
    cells = f_to_c.links(facet)
    # Get value of cells
    marked_cells = ct.values[cells]
    # Get the cell marked with 2
    correct_cell = np.flatnonzero(marked_cells == 2)

    assert len(correct_cell) == 1
    # Get local index of facet
    local_facets = c_to_f.links(cells[correct_cell[0]])
    local_index = np.flatnonzero(local_facets == facet)
    assert len(local_index) == 1

    # Append integration entities
    integration_entities.append(cells[correct_cell[0]])
    integration_entities.append(local_index[0])

# Create custom integration measure for one-sided integrals
breakpoint()
ds = ufl.Measure("ds", domain=mesh, subdomain_data=[
                 (8, np.asarray(integration_entities, dtype=np.int32))])
n = ufl.FacetNormal(mesh)
x = ufl.SpatialCoordinate(mesh)
# Exact integral is [y/2**2]_0^1= 1/2
L = ufl.dot(ufl.as_vector((x[1], 0)), n)*ds(8)
L_compiled = dolfinx.fem.form(L)

print(
    f"Correct integral: {mesh.comm.allreduce(dolfinx.fem.assemble_scalar(L_compiled), op=MPI.SUM)}")


# Create reference implementation where we use a restricted two-sided integral with no notion of orientation
dS = ufl.Measure("dS", domain=mesh, subdomain_data=ft, subdomain_id=2)
n = ufl.FacetNormal(mesh)
x = ufl.SpatialCoordinate(mesh)
L2 = ufl.dot(ufl.as_vector((x[1], 0)), n("+"))*dS
L2_compiled = dolfinx.fem.form(L2)
print(
    f"Wrong integral: {mesh.comm.allreduce(dolfinx.fem.assemble_scalar(L2_compiled), op=MPI.SUM)}")

Source: https://fenicsproject.discourse.group/t/wrong-facetnormal-vector-on-internal-boundaries/12887/2?u=dokken

[molel-gt]

I am trying this out on dolfinx=0.8.0 and running into the error
(8, np.asarray(integration_entities, dtype=np.int32))]) ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ ValueError: setting an array element with a sequence. The requested array has an inhomogeneous shape after 1 dimensions. The detected shape was (20,) + inhomogeneous part.

I believe that the line integration_entities.append(local_index) should instead be integration_entities.append(local_index[0])

[jorgensd]

@molel-gt I've updated the code to take this into account. Thanks for the feedback.

[pantolin]

I was wondering if there is not a neater way of doing this.

Beyond all the code needed above to find pairs of cells and local facets from the given facets, how this special way of defining tags (without a MeshTags object) is managed here and here using try/except seems a little bit hacky to me.

In my opinion, it would be just simpler to use ufl.ds with a provided list of (exterior and/or exterior facets). To the best of my knowledge, the only thing preventing this to work is the way in which fem::compute_integration_domains is implemented. In particular this part of the code that filters out non exterior facets.

However, removing this filter has a (non-minor) downside: if ufl.ds gets invoked without providing mesh tags, the integral gets evaluated for every facet in the domain, and not only on the exterior facets like it is the case right now. So, it changes what is somehow its default behaviour. This could be prevented through a flag (in the metadata of ufl.ds?) that defaults to the current behaviour.

[jorgensd]

I was wondering if there is not a neater way of doing this.

Beyond all the code needed above to find pairs of cells and local facets from the given facets, how this special way of defining tags (without a MeshTags object) is managed here and here using try/except seems a little bit hacky to me.

I am not sure what you are getting at here. The last of the try excepts: https://github.com/FEniCS/dolfinx/blob/4627a17522e73020dc7982e1a030da5297605f09/python/dolfinx/fem/forms.py#L359C1-L361C56 is due to the fact that ufl accepts both integer subdomain markers and tuples of integers as subdomain markers (which is very convenient when wanting to integrate over many subdomains, i.e. ds((2,5,7,8,3,1)) instead of having to write ds(2), ds(5), ds(7).....

For the first try except: This is how it has been decided to handle multiple type inputs in DOLFINx to functions (to avoid having duplicate signatures all around). One could have used:

if isinstance(subdomain, MeshTags):
   # do something
elif isinstance(subdomain, list):
   # do other thing
else:
   raise RuntimeError("Unknown input")

but some of the developers do not like the "isinstance" handling, as one doesn't want to strongly type Python interface,
and be closer to a Protocol interface: https://typing.python.org/en/latest/spec/protocol.html.

In my opinion, it would be just simpler to use ufl.ds with a provided list of (exterior and/or exterior facets). To the best of my knowledge, the only thing preventing this to work is the way in which fem::compute_integration_domains is implemented. In particular this part of the code that filters out non exterior facets.

I am not sure how you would do that for an interior integral? In the example above I do a one-sided facet integration on an interior boundary. This will always be connected to two cells, and one has to decide what cell one would like the integral to be over. Sending in a facet index by itself cannot determine this, it has to be a (cell, local_index) tuple.

However, removing this filter has a (non-minor) downside: if ufl.ds gets invoked without providing mesh tags, the integral gets evaluated for every facet in the domain, and not only on the exterior facets like it is the case right now. So, it changes what is somehow its default behaviour. This could be prevented through a flag (in the metadata of ufl.ds?) that defaults to the current behaviour.

Adding it to the metadata seems to me as a more hacky way of doing it.

Please note that this is a discussion that could be elevated into DOLFINx itself, and not the tutorial.

[pantolin]

I was wondering if there is not a neater way of doing this.
Beyond all the code needed above to find pairs of cells and local facets from the given facets, how this special way of defining tags (without a MeshTags object) is managed here and here using try/except seems a little bit hacky to me.

I am not sure what you are getting at here. The last of the try excepts: https://github.com/FEniCS/dolfinx/blob/4627a17522e73020dc7982e1a030da5297605f09/python/dolfinx/fem/forms.py#L359C1-L361C56 is due to the fact that ufl accepts both integer subdomain markers and tuples of integers as subdomain markers (which is very convenient when wanting to integrate over many subdomains, i.e. ds((2,5,7,8,3,1)) instead of having to write ds(2), ds(5), ds(7).....

Running your example above with the debugger brings me to that line. My guess is that integral.subdomain_id() evaluates to8, that is different from everywhere, so you enter that try/except. Once inside the try, for sid in integral.subdomain_id() fails as 8 is not iterable, so you end up at the except.

For the first try except: This is how it has been decided to handle multiple type inputs in DOLFINx to functions (to avoid having duplicate signatures all around). One could have used:

if isinstance(subdomain, MeshTags):

do something

elif isinstance(subdomain, list):

do other thing

else:
raise RuntimeError("Unknown input")
but some of the developers do not like the "isinstance" handling, as one doesn't want to strongly type Python interface, and be closer to a Protocol interface: https://typing.python.org/en/latest/spec/protocol.html.

Wouldn't be my choice, but I understand this is an agreed decision.

In my opinion, it would be just simpler to use ufl.ds with a provided list of (exterior and/or exterior facets). To the best of my knowledge, the only thing preventing this to work is the way in which fem::compute_integration_domains is implemented. In particular this part of the code that filters out non exterior facets.

I am not sure how you would do that for an interior integral? In the example above I do a one-sided facet integration on an interior boundary. This will always be connected to two cells, and one has to decide what cell one would like the integral to be over. Sending in a facet index by itself cannot determine this, it has to be a (cell, local_index) tuple.

Yes, I agree that finding the right side of interior facets may not be straightforward, but once found that side is uniquely identified by an index that can be fed directly to ds/dS.

However, removing this filter has a (non-minor) downside: if ufl.ds gets invoked without providing mesh tags, the integral gets evaluated for every facet in the domain, and not only on the exterior facets like it is the case right now. So, it changes what is somehow its default behaviour. This could be prevented through a flag (in the metadata of ufl.ds?) that defaults to the current behaviour.

Adding it to the metadata seems to me as a more hacky way of doing it.

I agree this solution may be far from ideal.

Please note that this is a discussion that could be elevated into DOLFINx itself, and not the tutorial.

You're right, it came natural to write here, just below the example.