Error during interpolation using lagrange elements

dolfinx
1

Hello,

yesterday I downloaded the latest (nightly) build of the dolfinX docker container and rerun one of my scripts, that worked correctly, using dolfinX v0.2.0a0. Using the new version of dolfinX, my script failed, when I tried to interpolate a function at a set of points. Some testing showed, that my interpolation routine (based on Implementation — FEniCSx tutorial) worked correctly on courser meshes but failed, when I refine them. This failure only happens, when I run the test case within a docker container started from my Ubuntu 20.04 machine. Switching to a Debian machine seams to overcome this issue.

I set up a short test case where I defined a rectangular domain, initialized a function F1 by a constant value and tried to interpolate the Function values at every coordinate of a degree of freedom. Using the latest version of dolfinX, this worked fine for meshes <=90 elements per edge. Thereafter, the function compute_collisions_point is not able to locate some points, which should belong to the considered mesh. These points are mostly located on the edges of the domain and have coordinates like e.g. x[0]=-9.15E-19 (value taken from a mesh with 91 elements).

Could there be a problem with the tolerances, used for comparison of point and element coordinates or are there any suggestions how I could fix this issue?

Best wishes,
Maximilian

import numpy as np
from mpi4py import MPI
from petsc4py import PETSc

import dolfinx
import dolfinx.geometry
import dolfinx.io
import ufl
from dolfinx.cpp.mesh import CellType

# Pointwise evaluation of function
def eval_poitvalues(uh, x):
    """
    Interpolation of dolfinx.Function at specified point.
    Source: https://jorgensd.github.io/dolfinx-tutorial/chapter1/membrane_code.html

    *Arguments:*
            uh: function from which quantities spud be interpolated
            x:  spacial coordinates of points (ndarray, size: nPoints x dimSpace)

    *Output*
            values: Interpolated values at specified points (ndarray, size: nPoints x dimFuncSpace)
    """
    # Get mesh
    mesh = uh.function_space.mesh

    # Create searchable data structure of current mesh
    bb_tree = dolfinx.geometry.BoundingBoxTree(mesh, mesh.topology.dim)

    # Initialize subset arrays
    cells = []
    points_on_proc = []

    # Select points on current processor
    for point in x:
        # Find cells that are close to the point
        cell_candidates = dolfinx.geometry.compute_collisions_point(bb_tree, point)
        # Choose one of the cells that contains the point
        cell = dolfinx.geometry.select_colliding_cells(mesh, cell_candidates, point, 1)
        # Only use evaluate for points on current processor
        if len(cell) == 1:
            points_on_proc.append(point)
            cells.append(cell[0])

    # Interpolate values at points on current processor
    points_on_proc = np.array(points_on_proc, dtype=np.float64)
    values = uh.eval(points_on_proc, cells)

    return values

# Create mesh
n_elmt = 91
mesh   = dolfinx.RectangleMesh(MPI.COMM_WORLD, [[0, 0, 0], [1, 1, 0]], [n_elmt, n_elmt], cell_type=CellType.quadrilateral)

# Create function spaces
order_k = 1
Pk      = ufl.VectorElement('CG', mesh.ufl_cell(), order_k)
V1      = dolfinx.FunctionSpace(mesh, Pk)


# Create interpolation functions
F1_k   = dolfinx.Function(V1)
F2_k   = dolfinx.Function(V1)

# Interpolate function into F1_k
with F1_k.vector.localForm() as F1_loc:
    F1_loc.set(2.5)

# Interpolate F1_k into F2_k
x_target          = V1.tabulate_dof_coordinates()
funcVal           = eval_poitvalues(F1_k, x_target)

if int(len(x_target)*2) == len(funcVal.ravel()):
    print('interpolation correct')
    with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "testInt-F1k.xdmf", "w") as xdmf:
        xdmf.write_mesh(mesh)
        xdmf.write_function(F1_k)
    F2_k.vector.array = funcVal.ravel()
    with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "testInt-F2k.xdmf", "w") as xdmf:
        xdmf.write_mesh(mesh)
        xdmf.write_function(F2_k)
else:
    print(int(len(x_target)*2))
    print(len(funcVal.ravel()))
2

A slight modification to your code yields the correct result:

# Pointwise evaluation of function
def eval_pointvalues(uh, x, distance_tolerance = 1e-15):
    """
    Interpolation of dolfinx.Function at specified point.
    Source: https://jorgensd.github.io/dolfinx-tutorial/chapter1/membrane_code.html

    *Arguments:*
            uh: function from which quantities spud be interpolated
            x:  spacial coordinates of points (ndarray, size: nPoints x dimSpace)
            distance_tolerance: Tolerance for how close an entity has to be to consider it a collision

    *Output*
            values: Interpolated values at specified points (ndarray, size: nPoints x dimFuncSpace)
    """
    # Get mesh
    mesh = uh.function_space.mesh

    # Create searchable data structure of current mesh
    bb_tree = dolfinx.geometry.BoundingBoxTree(mesh, mesh.topology.dim)

    # Initialize subset arrays
    cells = []
    points_on_proc = []

    # Select points on current processor
    for point in x:
        # Find cells that are close to the point
        cell_candidates = dolfinx.geometry.compute_collisions_point(bb_tree, point)
        if len(cell_candidates)==0:
            entity, distance = dolfinx.geometry.compute_closest_entity(bb_tree, point, mesh)
            if distance < distance_tolerance:
                cell = [entity]
        else:
            # Choose one of the cells that contains the point
            cell = dolfinx.geometry.select_colliding_cells(mesh, cell_candidates, point, 1)
        # Only use evaluate for points on current processor
        if len(cell) == 1:
            points_on_proc.append(point)
            cells.append(cell[0])

    # Interpolate values at points on current processor
    points_on_proc = np.array(points_on_proc, dtype=np.float64)
    values = uh.eval(points_on_proc, cells)

    return values
1 Like
3

Hallo dokken,

thanks for your reply. Your modification solved my problem.

Best wishes,
Maximilian