Hi,
I’m updating some old code and found that the Function.compute_point_values() method has been deprecated (dolfinx image 7387132cd03b, 22Jul22). I was trying to move a mesh as described by Nate here, how can I do this now?
Thanks,
Alex
Hi,
I’m updating some old code and found that the Function.compute_point_values() method has been deprecated (dolfinx image 7387132cd03b, 22Jul22). I was trying to move a mesh as described by Nate here, how can I do this now?
Thanks,
Alex
Managed to work it out myself. Here’s the solution for future reference:
def move_mesh(mesh,u,dt):
'''Move 2D mesh using finite element velocity vector field.
Cf the tutorial at https://jorgensd.github.io/dolfinx-tutorial/chapter1/membrane_code.html'''
x = mesh.geometry.x
gdim = mesh.geometry.dim
points = x.T
# Initialise cell search
bb_tree = BoundingBoxTree(mesh, mesh.topology.dim)
cells = []
points_on_proc = []
# Find cells whose bounding-box collide with the the points
cell_candidates = compute_collisions(bb_tree, points.T)
# Choose one of the cells that contains the point
colliding_cells = compute_colliding_cells(mesh, cell_candidates, points.T)
for i, point in enumerate(points.T):
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)
# Evaluate u
u_values = u.eval(points_on_proc, cells)
#Add zero values for 3rd dimension
if gdim==2:
z_values = np.zeros(u_values.shape[0])
u_values = np.vstack((u_values.T,z_values)).T
#Move mesh coordinates
x += dt*u_values
return mesh
It’s probably a lot easier to interpolate the displacement field or an absolute coordinate change to the mesh’s coordinate element:
import dolfinx
from mpi4py import MPI
import numpy as np
mesh = dolfinx.mesh.create_unit_square(MPI.COMM_WORLD, 4, 4)
Vex = mesh.ufl_domain().ufl_coordinate_element()
V = dolfinx.fem.FunctionSpace(mesh, Vex)
u = dolfinx.fem.Function(V)
u.interpolate(lambda x: np.stack((x[1]**2, x[0]**2)))
mesh.geometry.x[:,:mesh.geometry.dim] = u.x.array.reshape((-1, mesh.geometry.dim))
Serial:
Parallel: