Function.compute_point_values() deprecated, how to replace when moving mesh?

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:

4 Likes