Second order meshes interpolation error

General
, ,
5

Good to know this improved way to fastly assign the whole function with the use of mesh.entities_to_geometry and contrain interpolation region with cells0 in xyz.interpolate.

Thanks for confirming this. I believe this is the reason to cause the negative interpolated node in my application above. The following validation code shows a simple test, which maybe useful for the prospective people: a mesh with only one single triangle cell is made and a P2 function with all greater than 0 values is built on it. Several points inside the mesh are evaluated to be negative values as interpolated points on the new mesh:

# P2_eval.py
# $ python P2_eval.py
# Test on DOLFINx v0.9.0

import numpy as np
from mpi4py import MPI
from dolfinx import fem, mesh, io, geometry

import basix.ufl
import dolfinx
import ufl
import matplotlib.pyplot as plt

def evaluate_at_points(_mesh, f, points):

    if points.shape[1] == 2:
        z_coords = np.zeros((points.shape[0], 1))
        points = np.hstack([points, z_coords])

    bb_tree = geometry.bb_tree(_mesh, _mesh.topology.dim)
    cell_candidates = geometry.compute_collisions_points(bb_tree, points)
    colliding_cells = geometry.compute_colliding_cells(_mesh, cell_candidates, points)

    cells = []
    points_on_proc = []
    for i, point in enumerate(points):
        if len(colliding_cells.links(i)) > 0:
            cells.append(colliding_cells.links(i)[0])
            points_on_proc.append(point)

    values = f.eval(np.array(points_on_proc), np.array(cells))
    
    return values, np.array(points_on_proc)

comm = MPI.COMM_WORLD
nodes = np.array(
    [[0,   0  ],
     [1,   1  ],
     [0,   2  ],
     [0.5, 1.5],
     [0,   1  ],
     [0.5, 0.5],
    ],
    dtype=np.float64,
)

connectivity = np.array([[0, 1, 2, 3, 4, 5]], dtype=np.int64)

c_el = ufl.Mesh(basix.ufl.element("Lagrange", "triangle", 2, shape=(nodes.shape[1],)))
_mesh = dolfinx.mesh.create_mesh(MPI.COMM_SELF, cells=connectivity, x=nodes, e=c_el)

Q = fem.functionspace(_mesh, ("Lagrange", 2))
f = fem.Function(Q)

def x_tip(x):
    return x[0] > 0.9

corner_point = mesh.locate_entities_boundary(_mesh, 0, x_tip)
dof = fem.locate_dofs_topological(Q, 0, corner_point)
f.x.array[dof] = 1

n = 100
points = np.ones((n, 2))
points[:, 0] = np.linspace(0, 1, n)
vals, _ = evaluate_at_points(_mesh, f, points)

fig = plt.figure()
plt.plot(points[:, 0], vals, "o-")
plt.xlabel("x")
plt.ylabel("value")
plt.legend("y = 1 probe")
plt.grid("on")
plt.savefig("evals.png", dpi = 300)

with io.XDMFFile(comm, "f.xdmf", "w") as xdmf:
    xdmf.write_mesh(_mesh)
    xdmf.write_function(f)

Yielding

image
image1652×746 92.3 KB