The problem of subdomain partitioning

General
1

Hi,guys! I would like to inquire about the customization of ds and dS
Firstly, we perform rectangular partitioning on a disk:（The schematic diagram is as follows）

1
11254×847 142 KB

Now let’s define ds and dS in dolfinx:

boundaries = [(1, lambda x: np.isclose(x[0]**2+x[1]**2, ro**2)),
              (2, lambda x: np.isclose(x[0]**2+x[1]**2, ri**2))]

interiors = [(22,lambda x: (x[0]**2+x[1]**2<=ro**2+0.0001)&(x[0]**2+x[1]**2>=ri**2-0.0001)),
             (23,lambda x: x[0]**2+x[1]**2<=ri**2-0.0001)]

ds,facet_tag=creat_measure(“ds”,msh,fdim,boundaries)
dS,interiors_tag=creat_measure(“dS”,msh,fdim,interiors)
ds_oc=ds((1,2))
dS_oc=dS(22)

This is my ideal definition result:

2
21473×1056 146 KB

But I found that in the code,ds is a subset of dS.
If we define：

interiors = [(22,lambda x: (x[0]**2+x[1]**2<=ro**2-0.0001)&(x[0]**2+x[1]**2>=ri**2+0.0001)),

In this case ds and dS is disjoint .（The schematic diagram is as follows）
But the result seems to be missing some facets.

3
31333×1070 140 KB

I don’t know how dolfinx solved this problem.

2

ds will always be exterior facets, if you do not use custom integration measures (as done in:
Add one-sided integration example as example of custom integration · Issue #158 · jorgensd/dolfinx-tutorial · GitHub)

dS are reserved for two-sided facet integrals, where a facet is connected to two cells.

For either ds or ds, you can provide custom integration domains as illustrated by the aforementioned link in the post.

3

Thanks dokken！
I use custom integration measures:

def creat_measure(integral_type,msh,dim,entitiesmarkers):
    entities_indices, entities_markers = [], []
    for (marker, locator) in entitiesmarkers:
        entities = mesh.locate_entities(msh, dim, locator)
        entities_indices.append(entities)
        entities_markers.append(np.full_like(entities, marker))
    entities_indices = np.hstack(entities_indices).astype(np.int32)
    entities_markers = np.hstack(entities_markers).astype(np.int32)
    sorted_entities = np.argsort(entities_indices)
    entities_tag = mesh.meshtags(msh, dim, entities_indices[sorted_entities], entities_markers[sorted_entities])
    if dim==msh.topology.dim-1:
        msh.topology.create_connectivity(msh.topology.dim-1, msh.topology.dim)
    return ufl.Measure(integral_type, domain=msh, subdomain_data=entities_tag),entities_tag

ds_oc and dS_oc was created by:

ds,facet_tag=creat_measure(“ds”,msh,fdim,boundaries)
dS,interiors_tag=creat_measure(“dS”,msh,fdim,interiors)
ds_oc=ds((1,2))
dS_oc=dS(22)

Is there any difference between ds and dS when using custom integration measures?

4

Yes. dS instructs the code generation such that the generated kernels takes in two cells, one on either side of the facet in question, while ds is a one-sided integral, only taking in a single cell index and its local facet index.

5

Yes, I understand.
What I am worried about is whether the surface of the ds integral will produce incorrect results if it is the boundary of the sub mesh, not the boundary of the mesh?（like this:)

2
21473×1056 146 KB

I know the facets integral result will be assigned to the two units of the link.
dS on the boundary of the sub mesh(The red line in the picture)，Will integral result be allocated back to the two units of the link?
The distinction between ds and dS seems to rely on the topological information of the mesh.