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
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:
# To find the vertices of the outer potato:
# 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(
pressure_element = ufl.VectorElement(
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(
# 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(
# The indexes of the DOFs in the inner potato
inner_dof_idx = dolfinx.fem.locate_dofs_topological(
# 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(
# 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(
# 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.
# 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(
# 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
scipy.spatial.KDTree, but I am slightly confused by the arguments it accepts. The same goes for
For example, if I did:
idx = dolfinx.fem.locate_dofs_topological(
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.
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.