I run into a bit surprising implementation issue when trying to implement a material parameter function with distinct expressions on different subdomains, using version `0.7.3`

. Namely, I want to work with a discretization of

lets say in 2D, with k for simplicity being constant on each \Omega_i \subset \Omega, \bigcup_i \Omega_i = \Omega. My motivation is testing of numerical methods and inverse analysis so I would like to be able to change \Omega_i's and values of k in a comfortable way.

I found out that his issue had already been discussed on the forum, for example here or here and is covered in the tutorial.

Unfortunately, the solutions that I had found are somewhat lacking. The linked forum posts propose a solution based on the numbering of cells being same numbering of dofs for P(0) function. This works for now but might break with some future release.

There is an option of writing the function k as an expression that is interpolated onto a suitable `dolfinx.fem.Function`

but writing the expression gets a bit cumbersome as the subdomain partitioning must be afaik encoded into Heaviside functions since the interpolation is not accepting `if else`

contructions.

The tutorial proposes a usage of a marker function (effectively a characteristic fuction for each \Omega_i) and `dolfinx.mesh.locate_entities`

to define the subdomains. This approach seemed the most systematic but from my experiments and if I understand the code correctly, `dolfinx.mesh.locate_entities`

considers a mesh element to satisfy the marker condition only if **all** of its vertices satisfy the condition. This creates a situation where even the tutorial division of square two halves using marker functions

```
def Omega_0(x):
return x[1] <= 0.5
def Omega_1(x):
return x[1] >= 0.5
```

to the corresponding cells by

```
cells_0 = locate_entities(mesh, mesh.topology.dim, Omega_0)
cells_1 = locate_entities(mesh, mesh.topology.dim, Omega_1)
```

runs into troubles for mesh generated by

```
mesh = create_unit_square(MPI.COMM_WORLD, 11, 11)
```

as there is now a whole â€ślineâ€ť of cells that belong to neither cell list/subdomain. I wanted to fix this by writing a custom `locate_entities`

function that would only check midpoints of elements but I am not capable of that since `locate_entities`

is implemented in C++.

The last option I see would be to make a distinct mesh that respects the subdomain boundaries for each distinct subdomain configuration. This again is a bit cumbersome.

My question is if there some better way? If there is not, which approach is going to be supported in future development of dolfinx?