Preffered supported way of implementing a varying material parameters

variational formulation
,
1

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

\int_\Omega k(x) \left(\nabla u, \nabla v\right) \mathrm{d}x = \int_\Omega f v \mathrm{d}x,

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?

2

You can get the midpoints with

dolfinx.mesh.compute_midpoints(mesh: Mesh, dim: int, entities: npt.NDArray[np.int32])

So pseudocode would be (I can’t run this right now):

num_cells = mesh.topology.index_map(mesh.topology.dim).size_local
mps = dolfinx.mesh.compute_midpoints(mesh, mesh.topology.dim, np.arange(num_cells, type=np.int32))
cells_0 = np.where(mps[:,1] < 0.5)[0]
cells_1 = np.where(mps[:,1] >= 0.5)[0]

You can extend this to be compatible with your Omega_* functions by taking the transpose of the midpoints.

2 Likes
3

Thanks, this works nicely to fix my current issues.

However, in best case, it would be nice to know what is the approach that you as developers and maintainers prefer and want to support.

4

The aim of DOLFINx is to make it easy for users to choose whatever method they seem fit, and not be bound to a “one-fits-all” approach.
As you would like to check midpoints, using compute_midpoints seems like a natural fit.
If you want to combine midpoints and vertices in checks, I would call them separately and use numpy set operations (Set routines — NumPy v1.26 Manual) to find intersections, unions etc.

1 Like
5

That is quite an ambitious goal, I wish you all the best in realization. Thanks for the numpy link!