Evaluate fem.Function at point

dolfinx
1

Hello, I would like to evaluate a fem.Function at a point. My code is the following:

import dolfinx
import mpi4py
import dolfinx.fem as fem
import numpy as np

L = 1. # total length
d = L/10. # thickness
h = d/10. # size of a cell
my_domain = dolfinx.mesh.create_rectangle(comm=mpi4py.MPI.COMM_WORLD,points=((0.0, -0.5d), (L, 0.5d)), n=(int(L/h), int(d/h)),cell_type=dolfinx.mesh.CellType.triangle)

V = fem.functionspace(my_domain, (“Lagrange”, 1, (my_domain.geometry.dim, )))
u_test=fem.Function(V)
evaluate_at_points(np.array([[1.0L], [0.5d], [0]]), u_test)

import mpi4py
import numpy as np
import dolfinx # FEM in python
from dolfinx import fem
from typing import Optional
import ufl
import typing

def evaluate_at_points(points:np.ndarray|typing.Sequence, function: fem.Function) -> Optional[np.ndarray]:
    mesh = function.function_space.mesh

    if isinstance(points, list):
        points = np.array(points)
    if len(points.shape) == 1:
        if points.size < 3:
            # Pad point with zeros
            _points = np.zeros((3,1), dtype=mesh.geometry.x.dtype)
            _points[:len(points),0] = points
        else:
            if len(points) % 3 != 0 and mesh.geometry.dim == 2:
                # Pad points with extra 0
                _points = np.zeros((3, len(points)//2), dtype=mesh.geometry.x.dtype)
                _points[:2, :] = np.array(points).reshape(2, -1)
            elif len(points) % 3 == 0:
                _points = np.zeros((3, len(points)//3), dtype=mesh.geometry.x.dtype)
                _points[:, :] = np.array(points).reshape(3, -1)
            else:
                raise RuntimeError("Received list of points that cannot be formatted as a (n, 3) array")
    else:
        _points = np.array(points, dtype=mesh.geometry.x.dtype)        

    comm = mesh.comm
    if comm.rank == 0:
        input_points = _points.T
    else:
        input_points = np.empty((0, 3), dtype=points.dtype)

    src_owner, dest_owner, dest_points, dest_cells = dolfinx.cpp.geometry.determine_point_ownership(mesh._cpp_object, input_points, 1e-6)
    values = function.eval(np.array(dest_points).reshape(-1, 3), dest_cells).reshape(-1, function.function_space.dofmap.bs)
    if comm.rank != 0:
        assert np.allclose(dest_owner, 0)
    gathered_values = comm.gather(values, root=0)
    src_counter = np.zeros(len(np.unique(src_owner)), dtype=np.int32)
    bs = function.function_space.dofmap.bs
    values = np.zeros((input_points.shape[0], bs),dtype=function.x.array.dtype)
    if comm.rank ==0:
        for i, owner in enumerate(src_owner):
            if owner == -1:
                print(f"Could not find point in mesh for {input_points[i]}")
                continue
            values[i] = gathered_values[owner][src_counter[owner]]
            src_counter[owner] += 1
    else:
        return None    
    return values

I get the following error:

src_owner, dest_owner, dest_points, dest_cells = dolfinx.cpp.geometry.determine_point_ownership(mesh._cpp_object, input_points, 1e-6)^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

TypeError: cannot unpack non-iterable dolfinx.cpp.geometry.PointOwnershipData_float64 object

Any help would be very welcome!

2

We have just improved the Python API for compute point ownership. It is now a structure:

mimicking the structure of the C++ object.

One should now use the Python-wrapped documented function: dolfinx/python/dolfinx/geometry.py at 1ef0f824f34bcead9ee851b1a333309a1a34dd8d · FEniCS/dolfinx · GitHub

1 Like