Second order meshes interpolation error

General
, ,
1

I simulate ALE FSI problem in 3D with 2nd order meshes. In my framework, a remeshing process is needed to interpolate data from the old mesh to the new mesh in every several steps. Here a strange interpolation result is presented and seek for help on solving this issue.

To demonstrate this issue, two test sets of meshes (mesh1_order[1|2].bp and mesh2_order[1|2].bp) are provided – https://drive.google.com/file/d/1YjC6tCtiJaxQ7Sg7FVakZ8or23QNNMAc/view?usp=sharing. The MWE below interpolates function f1 on mesh1 to f2 on mesh2. The problem occurs on the interpolated result f2.

In the mesh.bp files, mesh and cell tags are stored. The mesh is a cubic domain ranging from [-1, -1, -1] to [1, 1, 1]. A sphere with radius 0.5 is centered at [0, 0, 0]. The cell tags for the sphere and the rest part are 0 and 1 to represent solid and fluid tags, respectively.

The result can be found in the attached figure. I tested it on DOLFINx v0.9.0. The original function f1 is all 0 around and with [0, 1] value in the spherical solid part. To compare the interpolation performance, 1st order and 2nd order meshes are used. I can sense that the wiggling result in 1st order (figure g,h) is due to numerical error with lower resolution on the sphere. But it seems a mistake is in 2nd order result (figure c,d) as [-0.12, 0] values on f2 is computed in the halo region of the sphere. I have no clue why there is a miscalculation on the second order meshes interpolation.

demo
demo1280×720 185 KB

MWE:

# interpolation.py
# Run with `mpirun -n 3 python interpolation.py`
# Test on DOLFINx v0.9.0

import os
os.environ['OMP_NUM_THREADS'] = '1'

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

def assign_xyz(submesh, xyz) -> None:
    xlims = []
    xmins = []
    for d in range(ndims):
        geom = submesh.geometry.x[:, d]
        xmin = geom.min() if len(geom) > 0 else 1e6
        xmax = geom.max() if len(geom) > 0 else -1
        xmin = comm.allreduce(xmin, op = MPI.MIN)
        xlim = comm.allreduce(xmax, op = MPI.MAX) - xmin
        xmins.append(xmin)
        xlims.append(xlim)

    def exp_xyz(x):
        r = []
        for d in range(ndims):
            r.append((x[d]-xmins[d])/xlims[d])
        return r

    xyz.interpolate(exp_xyz)

def nonmatching_meshes_interpolation(dest: fem.Function, src: fem.Function):
    V_dest = dest.function_space
    V_src = src.function_space
    tdim = V_src.mesh.topology.dim
    dest_mesh_cell_map = V_dest.mesh.topology.index_map(tdim)
    num_cells_on_proc = dest_mesh_cell_map.size_local + dest_mesh_cell_map.num_ghosts
    cells = np.arange(num_cells_on_proc, dtype = np.int32)
    interpolation_data = fem.create_interpolation_data(V_dest, V_src, cells, padding=1e-8)
    dest.interpolate_nonmatching(src, cells, interpolation_data=interpolation_data)
    dest.x.scatter_forward()


comm = MPI.COMM_WORLD
ndims = 3
order = 2 # 1 or 2
file1 = f"mesh1_order{order}.bp"
file2 = f"mesh2_order{order}.bp"

mesh1 = adios4dolfinx.read_mesh(file1, comm, engine="BP4")
mesh2 = adios4dolfinx.read_mesh(file2, comm, engine="BP4")

celltags1 = adios4dolfinx.read_meshtags(file1, mesh1, meshtag_name = "celltags", engine = "BP4")
celltags2 = adios4dolfinx.read_meshtags(file2, mesh2, meshtag_name = "celltags", engine = "BP4")

V1 = fem.functionspace(mesh1, ("Lagrange", order, (ndims,)))
V2 = fem.functionspace(mesh2, ("Lagrange", order, (ndims,)))

f1 = fem.Function(V1)
f2 = fem.Function(V2)

submesh1, cellmap1, _, _ = mesh.create_submesh(mesh1, ndims, celltags1.find(0))
submesh2, cellmap2, _, _ = mesh.create_submesh(mesh2, ndims, celltags2.find(0))

xyz1 = fem.Function(fem.functionspace(submesh1, ("Lagrange", order, (ndims,))))
xyz2 = fem.Function(fem.functionspace(submesh2, ("Lagrange", order, (ndims,))))

# make functions on submesh1 and mesh1
assign_xyz(submesh1, xyz1)
n = len(cellmap1)
f1.interpolate(xyz1, cells0=np.arange(n), cells1=cellmap1)
f1.x.scatter_forward()

# interpolate into f2 with nonmatching meshes interpolation
nonmatching_meshes_interpolation(f2, f1)

# make function on submesh2 with identical mesh interpolation
n = len(cellmap2)
xyz2.interpolate(f2, cells0=cellmap2, cells1=np.arange(n))
xyz2.x.scatter_forward()

# visualization
with io.XDMFFile(comm, f"f1_order{order}.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh1)
    xdmf.write_function(f1)

with io.XDMFFile(comm, f"f2_order{order}.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh2)
    xdmf.write_function(f2)

# assistant visualization
with io.XDMFFile(comm, f"xyz1_order{order}.xdmf", "w") as xdmf:
    xdmf.write_mesh(submesh1)
    xdmf.write_function(xyz1)

with io.XDMFFile(comm, f"xyz2_order{order}.xdmf", "w") as xdmf:
    xdmf.write_mesh(submesh2)
    xdmf.write_function(xyz2)
2

This is a problem with your mesh.
If you inspect the files f1_order2.xdmf and f2_order2.xdmf carefully, you will observe that the geometry of f1_order2.xdmf has a weird kink in the offending cell

image
image1920×1085 405 KB

This causes interpolation for that node to happen in the cell that is on the outside of the sphere, yielding a negative value.
Here is the script Im running (in serial) with DOLFINx main branch:

# interpolation.py
# Run with `mpirun -n 3 python interpolation.py`
# Test on DOLFINx v0.9.0

import os

os.environ["OMP_NUM_THREADS"] = "1"

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


def assign_xyz(submesh, xyz) -> None:
    xlims = []
    xmins = []
    for d in range(ndims):
        geom = submesh.geometry.x[:, d]
        xmin = geom.min() if len(geom) > 0 else 1e6
        xmax = geom.max() if len(geom) > 0 else -1
        xmin = comm.allreduce(xmin, op=MPI.MIN)
        xlim = comm.allreduce(xmax, op=MPI.MAX) - xmin
        xmins.append(xmin)
        xlims.append(xlim)

    def exp_xyz(x):
        r = []
        for d in range(ndims):
            r.append((x[d] - xmins[d]) / xlims[d])
        return r

    xyz.interpolate(exp_xyz)


def nonmatching_meshes_interpolation(dest: fem.Function, src: fem.Function):
    V_dest = dest.function_space
    V_src = src.function_space
    tdim = V_src.mesh.topology.dim
    dest_mesh_cell_map = V_dest.mesh.topology.index_map(tdim)
    num_cells_on_proc = dest_mesh_cell_map.size_local + dest_mesh_cell_map.num_ghosts
    cells = np.arange(num_cells_on_proc, dtype=np.int32)
    interpolation_data = fem.create_interpolation_data(
        V_dest, V_src, cells, padding=1e-8
    )
    dest.interpolate_nonmatching(src, cells, interpolation_data=interpolation_data)
    dest.x.scatter_forward()


comm = MPI.COMM_WORLD
ndims = 3
order = 2  # 1 or 2
file1 = f"mesh1_order{order}.bp"
file2 = f"mesh2_order{order}.bp"

mesh1 = adios4dolfinx.read_mesh(file1, comm, engine="BP4")
mesh2 = adios4dolfinx.read_mesh(file2, comm, engine="BP4")

celltags1 = adios4dolfinx.read_meshtags(
    file1, mesh1, meshtag_name="celltags", engine="BP4"
)
celltags2 = adios4dolfinx.read_meshtags(
    file2, mesh2, meshtag_name="celltags", engine="BP4"
)

V1 = fem.functionspace(mesh1, ("Lagrange", order, (ndims,)))
V2 = fem.functionspace(mesh2, ("Lagrange", order, (ndims,)))

f1 = fem.Function(V1)
f2 = fem.Function(V2)

submesh1, cellmap1, _, _ = mesh.create_submesh(mesh1, ndims, celltags1.find(0))
submesh2, cellmap2, _, _ = mesh.create_submesh(mesh2, ndims, celltags2.find(0))

xyz1 = fem.Function(fem.functionspace(submesh1, ("Lagrange", order, (ndims,))))
xyz2 = fem.Function(fem.functionspace(submesh2, ("Lagrange", order, (ndims,))))

# make functions on submesh1 and mesh1
assign_xyz(submesh1, xyz1)
sub_to_parent1 = cellmap1.sub_topology_to_topology(
    np.arange(
        submesh1.topology.index_map(submesh1.topology.dim).size_local
        + submesh1.topology.index_map(submesh1.topology.dim).num_ghosts,
        dtype=np.int32,
    ),
    inverse=False,
)
n = len(sub_to_parent1)
f1.interpolate(xyz1, cells0=np.arange(n, dtype=np.int32), cells1=sub_to_parent1)
f1.x.scatter_forward()

# interpolate into f2 with nonmatching meshes interpolation
nonmatching_meshes_interpolation(f2, f1)

# make function on submesh2 with identical mesh interpolation
sub_to_parent2 = cellmap2.sub_topology_to_topology(
    np.arange(
        submesh2.topology.index_map(submesh2.topology.dim).size_local
        + submesh2.topology.index_map(submesh2.topology.dim).num_ghosts,
        dtype=np.int32,
    ),
    inverse=False,
)
n = len(sub_to_parent2)
xyz2.interpolate(f2, cells0=sub_to_parent2, cells1=np.arange(n, dtype=np.int32))
xyz2.x.scatter_forward()

# visualization
with io.XDMFFile(comm, f"f1_order{order}.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh1)
    xdmf.write_function(f1)

with io.XDMFFile(comm, f"f2_order{order}.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh2)
    xdmf.write_function(f2)

# assistant visualization
with io.XDMFFile(comm, f"xyz1_order{order}.xdmf", "w") as xdmf:
    xdmf.write_mesh(submesh1)
    xdmf.write_function(xyz1)

with io.XDMFFile(comm, f"xyz2_order{order}.xdmf", "w") as xdmf:
    xdmf.write_mesh(submesh2)
    xdmf.write_function(xyz2)
3

Could you please point out which point causes the problem? For me, gmsh is used in a normal way to generate this 2nd order mesh. Still lost in how to improve it.

My question is why a cell with all greater than 0 nodes could result in a negative interpolated node in 2nd order meshes interpolation. This situation does not happen in 1st order interpolation. Is that because a quadratic basis function of P2 element has negative values, which could cause the point just outside the solid to be evaluated as a negative one?

4

Yes, the values of basis functions for a second order basis function is not necessarily bounded by the values at the degrees of freedom.
For instance, if we just inspect f1 after interpolating from xyz1 to f1 with scifem.compute_extrema, we get:

import scifem

extrema, node = scifem.compute_extrema(f1.sub(0), min)

(Pdb++) extrema
np.float64(-0.14891819580613833)
(Pdb++) node
array([ 0.10721618, -0.09852694, -0.54036315])

which is at:

image
image1920×1085 272 KB

We can even do this without having a submesh to mesh interpolation by calling:

import os

os.environ["OMP_NUM_THREADS"] = "1"

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


def assign_xyz(domain, xyz, cells0=None) -> None:
    xlims = []
    xmins = []
    for d in range(ndims):
        if cells0 is not None:
            nodes = np.unique(
                mesh.entities_to_geometry(domain, domain.topology.dim, cells0).flatten()
            )
            geom = domain.geometry.x[nodes, d]
        else:
            geom = domain.geometry.x[:, d]
        xmin = geom.min() if len(geom) > 0 else 1e6
        xmax = geom.max() if len(geom) > 0 else -1
        xmin = comm.allreduce(xmin, op=MPI.MIN)
        xlim = comm.allreduce(xmax, op=MPI.MAX) - xmin
        xmins.append(xmin)
        xlims.append(xlim)

    def exp_xyz(x):
        r = []
        for d in range(ndims):
            r.append((x[d] - xmins[d]) / xlims[d])
        return r

    xyz.interpolate(exp_xyz, cells0=cells0)
    xyz.x.scatter_forward()


comm = MPI.COMM_WORLD
ndims = 3
order = 2
file1 = f"mesh1_order{order}.bp"

mesh1 = adios4dolfinx.read_mesh(file1, comm, engine="BP4")

celltags1 = adios4dolfinx.read_meshtags(
    file1, mesh1, meshtag_name="celltags", engine="BP4"
)

V1 = fem.functionspace(mesh1, ("Lagrange", order, (ndims,)))

f1 = fem.Function(V1)

assign_xyz(mesh1, f1, cells0=celltags1.find(0))

extrema, node = scifem.compute_extrema(f1.sub(0), min, num_candidates=10, p=1)
print(extrema, node)

# visualization
with io.XDMFFile(comm, f"f1_order{order}.xdmf", "w") as xdmf:
    xdmf.write_mesh(mesh1)
    xdmf.write_function(f1)
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