Distance Between Points and Boundary for Tagged Meshes

Hi! I haven’t quite found a final answer about this even tho I got pretty close with tips from other posts.

I have a mesh that is two “continuous potatoes”, one inside the other. Let’s say the mesh is a dolfinx.mesh.Mesh object called mesh. I also have a dolfinx.mesh.MeshTagsMetaClass object called tags. I can find the nodes of the inner potato with tags.find(1) and those of the outer potato with tags.find(2) (for example). For every point inside the outer potato, I want to find the distance between the point and the surface separating the inner potato from the outer potato.

What would be the best way to achieve this? I have attached a simple example mesh with spheres instead of potatoes to make my problem clearer. In the example I would like to compute the distance from the points inside the outer shell to the surface between the outer shell and the inner sphere.

For context, I need this distance to define a dolfinx.fem.Function. Basically, the value of the function at the nodes in the outer potato depends on that distance.

1 Like

You can do this in steps, i.e. start by using:

mesh_vertices_2 = dolfinx.cpp.mesh.entities_to_geometry(mesh._cpp_object, mesh.topology.dim-1,tags.find(2), False)
mesh_nodes_2 = mesh.geometry.x[mesh_vertices_2].reshape(-1, 3)
squared_distance = dolfinx.geometry.squared_distance(mesh, mesh.topology.dim-1, tags.find(1), mesh_nodes_2)

Note that if you want this as part of a function, you either need to assume a 1-1 correspondence between the mesh vertices and the function space(i.e. use a first order Lagrange space), or you need to use the dof-coordinates instead of the mesh vertices to compute these distances.

Amazing, thank you for your help @dokken ! I finally been able to get back to this.

I reckon I will need to do the dof-coordinates thing. With the help of the documentation and some friendly AI I managed to put something together. I cannot currently test it because I get my dolfinx from arch4edu and they haven’t update it yet for the latest python, which I have…

Anyway, I thought I would share an example of what I put together in case you or others have comments about it. I just put a couple of questions in the comments explaining what I am doing.

import numpy as np
from petsc4py.PETSc import ScalarType
import dolfinx
import ufl
from scipy.spatial import KDTree


# I have a mesh here:
mesh: dolfinx.mesh.Mesh # This can be generated in any way. In my case from a MED file.

# I also have tags:
tags: dolfinx.mesh.MeshTagsMetaClass # This can also be generated in any way

# To find the vertices of the inner potato:
tags.find(1)

# To find the vertices of the outer potato:
tags.find(2)

# So, I will go ahead and define my elements and function space.
# First, I will make sure that I am using petsc-complex (my PDE is complex):
assert np.issubdtype(ScalarType, np.complexfloating)

# I want to solve for complex pressure and particle velocity in acoustics. So I will make a mixed space:
pressure_element = ufl.FiniteElement(
    familiy='P',
    cell=mesh.ufl_cell(),
    degree=2
)

pressure_element = ufl.VectorElement(
    familiy='P',
    cell=mesh.ufl_cell(),
    degree=2,
    dim=mesh.topology.dim
)

element = ufl.MixedElement([pressure_element, velocity_element])

f_space = dolfinx.fem.FunctionSpace(mesh=mesh, element=element)

# Cool. Now, I would like to define an adiabatic absorber in the outer potato. That is, the wavenumber should become complex in the outer potato to attenuate sound in there.
ordinary_wave_number: float  # This is the ordinary wave number in the medium

shape_function: Callable[[np.ndarray], np,ndarray]  # This function gives the shape of the imaginary part of the wavenumber in the outer potato.

# The wave number should be given by a function. In the inner potato it is the ordinary one. In the outer it is ordinary + 1j shape_function

wave_number = dolfinx.fem.Function(
    V=f_space,
    name='wave_number',
    dtype=np.complexfloating
)

# I will define the values for the array of this function. Then we can assign them. I assume the function array is as big as the number of DOFs. Is this right?

values = np.full(wave_number.array.size, ordinary_wave_number, dtype=np.complexfloating)

# Let's add the contributions of the shape function to the values over the outer potato. The shape function depends on the distance between the point in the outer potato and the surface between the inner and outer potato. So, we must compute that.

# The coordinates of all the DOFs
dof_coords = f_space.tabulate_dof_coordinates()

# The indexes of the DOFs in the outer potato
outer_dof_idx = dolfinx.fem.locate_dofs_topological(
    V=f_space,
    entity_dim=mesh.topology.dim,
    entities=tags.find(2),
    remote=False
)

# The indexes of the DOFs in the inner potato
inner_dof_idx = dolfinx.fem.locate_dofs_topological(
    V=f_space,
    entity_dim=mesh.topology.dim,
    entities=tags.find(1),
    remote=False
)

# We now search for the nearest inner potato DOF to each DOF in the outer potato, and get the distance
tree = KDTree(data=dof_coordinates[inner_dof_idx].transpose())
distance, idx = tree.query(
    x=dof_coords[outer_dof_idx].transpose(), 
    k=1
)

# As far as I understand, I now have, for each DOF in the outer potato, the distance from that point and the nearest DOF in the inner potato. I guess things could be made more efficient by only searching for the nearest in the surface:

#  inner_dof_surf_idx = dolfinx.fem.locate_dofs_topological(
#      V=f_space,
#      entity_dim=mesh.topology.dim-1,
#      entities=tags.find(4),  # If 4 was the tag for the surface. Or, I could use 1, but that would assume there is only 1 2D surface around the inner potato. In the future I might do more complex things.
#      remote=False
#  )

# I couldn't quite figure out what to do with this. I feel like the arguments I chosen do not make sense:
# squared_distance = dolfinx.geometry.squared_distance(
#     mesh=mesh, 
#     dim=mesh.topology.dim, 
#     entities=tags.find(1), 
#     points=dof_coords[idx]
# )

# distance = np.sqrt(squared_distance)

# And we update the relevant values of the function array (this will only work if the values of the function array are the same order and size as the DOFs. Is that correct?)
values[outer_dof_idx] += shape_function(distance)
wave_number.x.array[:] = values

I guess I will only be able to figure this out when I get dolfinx to run again. In the meanwhile, does that make any sense?

Also, I have a feeling I should prefer dolfinx.geometry.squared_distance over scipy.spatial.KDTree, but I am slightly confused by the arguments it accepts. The same goes for locate_dofs_topological.

For example, if I did:

idx = dolfinx.fem.locate_dofs_topological(
    V=f_space,
    entity_dim=mesh.topology.dim-1,
    entities=tags.find(2),
    remote=False
)

Would this find all the DOFs on the 2D surfaces containing the outer potato? (assuming 3D mesh). I would have 2 of those. One in contact with the inner potato and one in contact with the rest of the universe.

And this:

dolfinx.geometry.squared_distance(
    mesh=mesh, 
    dim=mesh.topology.dim-1, 
    entities=tags.find(1), 
    points=dof_coords[idx]
)

Would this give me the distance between all the DOFs in said surfaces and the inner potato points? What is dim for? I couldn’t quite wrap my head around it with the docs.

This is not needed, it doesn’t do anything.

Pressure element defined twice, no velocity element defined.

As so much of the code clearly hasnt been gone through carefully, i will not spend more time looking through it.