Hi everyone,
I have met a problem in interpolating a function from a mesh to the new one. The version I use is dolfinx 0.7.2. The code is like:
domain_a = file_input1.read_mesh(name=“Grid”)
W_a = fem.functionspace(domain_a, (“CG”, 2))
u_a = fem.Function(W_a)
u_b.interpolate(u_a, nmm_interpolation_data=fem.create_nonmatching_meshes_interpolation_data(domain_b._cpp_object, W_b.element, domain_a._cpp_object, 0))
Some spots in the center of the interpolation picture can be seen, which represent the tolerance between the origin values and the interpolated values. I would like to have your ideas about how to eliminate this tolerance, or in another word, is there something I can do or adjust to make the interpolation perfect?
Thanks in advance!
Best regards,
dokken
November 27, 2023, 3:20pm
2
jeremych:
u_b.interpolate(u_a, nmm_interpolation_data=fem.create_nonmatching_meshes_interpolation_data(domain_b._cpp_object, W_b.element, domain_a._cpp_object, 0))
Increase the last argument of fem.create_nonmatching_meshes_interpolation_data from 0 to 1e-6
Hello Dr. Dokken,
Thank you so much for your answer.
Best regards,
dokken
November 28, 2023, 9:17am
4
When you have non matching meshes, you can have a point in one cell that is on the border of a set of cells in the other mesh.
The GJK collision algorithm that we use for finding colliding cells, does not work at machine precision, so this point may be 1e-7-1e-8 outside all cells in the new mesh.
opened 12:25PM - 19 Sep 23 UTC
closed 03:59PM - 03 Oct 23 UTC
bug
### How to reproduce the bug
I stumpled across this issue when using the interp… olate function for non-matching meshes. When I have only a 'small' shift of coordinates of the two meshes, the interpolate function returns a zero. This behavior is clear for boundary values, but it also happens for values inside both domains. I could track this issue back to the geometry.compute_colliding_cells function.
In the minimal example I presented the issue:
The 'evaluate' function I define is based on this fenicsx tutorial: https://jsdokken.com/dolfinx-tutorial/chapter1/membrane_code.html#making-curve-plots-throughout-the-domain
I am creating two meshes and shift one of those meshes slightly by a random distribution (I added the seed for reproducability). I am then defining some function on mesh 1 and interpolate it to mesh 2. At one position in the function defined on mesh 2 a zero is interpolated, although it is inside both meshes. You can see this clearly in paraview, I also print the problematic DOF to the console. I also write the two functions to paraview where you can have a closed look. Making the eps parameters smaller results in more of these candidates.
Then I am using the 'evaluate' function to further track down the problem (I guess the interpolate function is using something similar). While the cell_candidates are still correct (the bounding box tree finds both elements which might be colliding with the problematic point), the compute_colliding_cells function returns an empty array of cells. I guess because it cannot decide which cell to choose, so it chooses none.
### Minimal Example (Python)
```Python
from dolfinx import io, fem, mesh, geometry
from mpi4py import MPI
import numpy as np
# Evaluate function values as implemented in
# https://jsdokken.com/dolfinx-tutorial/chapter1/membrane_code.html#making-curve-plots-throughout-the-domain
def evaluate_function(points, u):
domain = u.function_space.mesh
bb_tree = geometry.BoundingBoxTree(domain, domain.topology.dim)
cells = []
points_on_proc = []
# Find cells whose bounding-box collide with the points
cell_candidates = geometry.compute_collisions(bb_tree, points)
print("Cell candidates (still correct):", cell_candidates.links(6166))
# Choose one of the cells that contains the point
colliding_cells = geometry.compute_colliding_cells(domain, cell_candidates, points)
print("Colliding cells (wrong):", colliding_cells.links(6166))
for i, point in enumerate(points):
if len(colliding_cells.links(i))>0:
points_on_proc.append(point)
cells.append(colliding_cells.links(i)[0])
points_on_proc = np.array(points_on_proc, dtype=np.float64)
values = u.eval(points_on_proc, cells)
return values
file1 = io.XDMFFile(MPI.COMM_WORLD, 'testout/output1.xdmf', 'w')
file2 = io.XDMFFile(MPI.COMM_WORLD, 'testout/output2.xdmf', 'w')
# Create two initially equal meshes
msh1 = mesh.create_rectangle(comm=MPI.COMM_WORLD,
points=((0.0, 0.0), (1.0, 1.0)), n=(100, 100),
cell_type=mesh.CellType.triangle)
msh2 = mesh.create_rectangle(comm=MPI.COMM_WORLD,
points=((0.0, 0.0), (1.0, 1.0)), n=(100, 100),
cell_type=mesh.CellType.triangle)
# Shift the geometry of the second mesh slighly in any direction
eps = 5.0e-7
np.random.seed(42)
rand_array = eps * np.random.rand(*msh2.geometry.x.shape)
rand_array[:, -1] = 0.0
msh2.geometry.x[:] += rand_array
V1 = fem.FunctionSpace(msh1, ('CG', 1))
V2 = fem.FunctionSpace(msh2, ('CG', 1))
u1 = fem.Function(V1)
u1.interpolate(lambda x: x[0] ** 2 + x[1]) # Arbitrary function value
u2 = fem.Function(V2)
u2.interpolate(u1)
np.set_printoptions(precision=12)
print("Problematic point of u2 (at index 6166)", msh2.geometry.x[6166])
print("Original point (u1 is defined at) ", msh1.geometry.x[6166])
print("Magnitude of shift:", np.linalg.norm(msh1.geometry.x[6166]-msh2.geometry.x[6166]))
print("Function value at u2[6166]:", u2.x.array[6166])
print("Function value at u1[6166]:", u1.x.array[6166])
# Test the evaluate function
vals = evaluate_function(msh2.geometry.x, u1)
# Write and close files
file1.write_mesh(msh1)
file1.write_function(u1)
file2.write_mesh(msh2)
file2.write_function(u2)
file1.close()
file2.close()
```
### Output (Python)
_No response_
### Version
0.6.0
### DOLFINx git commit
24f86a9ce57df6978070dbee22b3eae8bb77235f
### Installation
I used docker on windows 11
### Additional information
_No response_
