Interpolate_nonmatching gives wrong results

General
1

Dear all,

I am currently working on a geometry consisting on a set of connected beams arranged in a particular way. I need to define a normalized vector field which is parallel to the axis of each beam everywhere in the geometry. To do so I created two meshes: one for the full 3D body and a 1D mesh embedded in 3D space in which each beam is replaced by a straight line (sort of the skeleton of the structure). My idea is to use Jacobian to get the tangent in the 1D mesh and then interpolate it to the 3D mesh. However, when interpolating between the meshes, the interpolated vector field is wrong: while the original is always parallel to the 1D beams, the interpolated final field simply points at a single direction (see image):

MWE_both
MWE_both2500×1406 255 KB

The left figure is the correct vector field in the 1D mesh and the right figure is the incorrectly interpolated field in the 3D body mesh. Why is the vector field wrongly interpolated? Is there a better way to approach this? The dolfinx code (run in dolfinx_nightly through the docker image) is:

import dolfinx
import numpy as np
import ufl
from mpi4py import MPI
from basix.ufl import quadrature_element, element


mesh_name = 'pMM2x'

#Import Meshes from xdmf file 
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, mesh_name+".xdmf", 'r') as xdmf_infile:
    body = xdmf_infile.read_mesh(name = 'Grid')
    
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, mesh_name+"_skeleton.xdmf", 'r') as xdmf_infile:
    skeleton = xdmf_infile.read_mesh(name = 'Grid')

#Define relevant vector elements and spaces
el0 = element("DG", skeleton.basix_cell(), 0, shape=(3,))
V0 = dolfinx.fem.functionspace(skeleton, el0)
el1 = element("DG", body.basix_cell(), 0, shape=(3,))
V1 = dolfinx.fem.functionspace(body, el1)
padding = 1e-14

#Compute tangent vector 
dx_dX = ufl.Jacobian(skeleton)[:, 0]
t_ufl = dx_dX / ufl.sqrt(ufl.inner(dx_dX, dx_dX))
t = dolfinx.fem.Function(V0, name="Tangent_vector")
t.interpolate(dolfinx.fem.Expression(t_ufl, V0.element.interpolation_points))
t.x.scatter_forward()

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "skeleton_tangent.xdmf", "w") as xdmf_outfile:
    xdmf_outfile.write_mesh(skeleton)
    t.name = "Tangent"
    xdmf_outfile.write_function(t)

fine_mesh_cell_map = body.topology.index_map(body.topology.dim)
num_cells_on_proc = fine_mesh_cell_map.size_local + fine_mesh_cell_map.num_ghosts
cells = np.arange(num_cells_on_proc, dtype=np.int32)
interpolation_data = dolfinx.fem.create_interpolation_data(V1, V0, cells, padding=padding)

n = dolfinx.fem.Function(V1, name="Director")
n.interpolate_nonmatching(t, cells, interpolation_data)
n.x.scatter_forward()

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "body_tangent.xdmf", "w") as xdmf_outfile:
    xdmf_outfile.write_mesh(body)
    n.name = "Tangent"
    xdmf_outfile.write_function(n)

Thanks in advance for any help/insights that you can provide!

2

Without the meshes, there is no way to debug this.

Could you start by interpolating a scalar value from one grid to the other, to ensure that works?

3

I actually wanted to upload the h5 and xdmf files as well but apparently I can only upload image files in the forum (its my first time asking here). Is there a way to do so? In the meantime here is the gmsh python code that produces the two meshes:


import numpy as np
import gmsh
import meshio

# Initialize Gmsh
gmsh.initialize()

# Create a new model
gmsh.model.add("model3D")
gmsh.option.setNumber("Geometry.OCCBooleanSimplify", 0)
mesh_name = 'pMM2x'
#geom.characteristic_length_max = mesh_size
# Define parameters
r = 0.05  # Small radius
R = 0.25  # Big radius
H = 1.0   # Height
a = 4 * H / np.sqrt(3)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 2*r)
def add_link(ϕ, rotate_axis):
    # Define the three points
    points = [
        [0.0, 0.0, 0.0],
        [r, 0.0, 0.0],      # Point 1
        [R, 0.0, H / 2.0], # Point 2
        [r, 0.0, H],        # Point 3
        [0.0, 0.0, H]       # Point 4
    ]

    # Define the axis of rotation (e.g., the z-axis)
    axis = [0.0, 0.0, 1.0]
    point_on_axis = [0.0, 0.0, 0.0]

    # Add points to the model
    geom_points = [gmsh.model.occ.addPoint(*p) for p in points]
    sym_axis = gmsh.model.occ.addLine(geom_points[-1], geom_points[0])
    #print(sym_axis)

    # Create lines connecting the points
    lines = []
    for i in range(len(geom_points) - 1):
        lines.append(gmsh.model.occ.addLine(geom_points[i], geom_points[i + 1]))
    lines.append(gmsh.model.occ.addLine(geom_points[-1], geom_points[0]))

    
    # Create a wire (curve loop) from the lines
    wire = gmsh.model.occ.addCurveLoop(lines)

    # Create a surface from the wire
    surface = gmsh.model.occ.addPlaneSurface([wire])

    # Rotate the surface around the specified axis to create a 3D body
    θ = 2 * np.pi / 3
    s = surface
    bodies = []
    for i in range(3):
        # Revolve the surface
        result = gmsh.model.occ.revolve(
            [(2, s)],  # Surface to revolve
            point_on_axis[0], point_on_axis[1], point_on_axis[2],  # Point on axis
            axis[0], axis[1], axis[2],  # Axis of rotation
            θ  # Angle of rotation
        )
        # Extract the resulting entities
        top = result[0][1]  # The new surface created by the revolution
        body = result[1][1]  # The 3D volume created by the revolution
        bodies.append(body)
        s = top  # Update the surface for the next revolution

    # Add spheres at the ends
    s1 = gmsh.model.occ.addSphere(points[0][0], points[0][1], points[0][2], r)
    s2 = gmsh.model.occ.addSphere(points[-1][0], points[-1][1], points[-1][2], r)

    # Combine all bodies into one
    b = gmsh.model.occ.fuse([(3, s1), (3, s2)], [(3, body) for body in bodies])

    # Rotate the combined body by the specified angle and axis
    gmsh.model.occ.rotate(
        b[0],  # The combined body
        point_on_axis[0], point_on_axis[1], point_on_axis[2],  # Point on axis
        rotate_axis[0], rotate_axis[1], rotate_axis[2],  # Axis of rotation
        ϕ  # Angle of rotation
    )

    gmsh.model.occ.rotate(
        [(1,sym_axis)],  # The combined body
        point_on_axis[0], point_on_axis[1], point_on_axis[2],  # Point on axis
        rotate_axis[0], rotate_axis[1], rotate_axis[2],  # Axis of rotation
        ϕ  # Angle of rotation
    )

    return b[0], sym_axis  # Return the tag of the resulting body

def create_metaatom(translate):
    normal_planes = [[1.0, -1.0, 0.0], [1.0, -1.0, 0.0], [1.0, 1.0, 0.0], [1.0, 1.0, 0.0]]
    ϕ = [np.arccos(1 / np.sqrt(3)), 2.0 * np.pi - np.arccos(1 / np.sqrt(3)), np.arccos(-1 / np.sqrt(3)), 2.0 * np.pi - np.arccos(-1 / np.sqrt(3))]
    links = []
    sk = []

    for i, n in enumerate(normal_planes):
        link, ax = add_link(ϕ[i], n)
        #print(ax)
        gmsh.model.occ.translate([(1,ax)], translate[0], translate[1], translate[2])
        links.append(link[0][:])
        sk.append(ax)

    # Combine all links into one meta-atom
    meta_atom = gmsh.model.occ.fuse([links[0]], links[1:])
    # Translate the meta-atom
    gmsh.model.occ.translate(meta_atom[0], translate[0], translate[1], translate[2])

    return meta_atom[0][:], sk  # Return the tag of the resulting meta-atom

def create_SC_site(center):
    translations = [[0.0, 0.0, 0.0], [a / 2.0, a / 2.0, 0.0], [0.0, a / 2.0, a / 2.0], [a / 2.0, 0.0, a / 2.0]]
    unit_cell = []
    skeleton = []

    for t in translations:
        l, sk = create_metaatom(t)
        for link in sk:
            gmsh.model.occ.translate([(1,link)], center[0], center[1], center[2])
        #gmsh.model.occ.translate(sk, center[0], center[1], center[2])
        unit_cell.append(l[0][:])
        skeleton += sk
        #skeleton.append(sk[0][:])

    # Combine all meta-atoms into one unit cell
    u_cell = gmsh.model.occ.fuse([unit_cell[0]], unit_cell[1:])
    #cell_skeleton = gmsh.model.occ.fuse([skeleton[0]], skeleton[1:])

    # Translate the unit cell
    gmsh.model.occ.translate(u_cell[0], center[0], center[1], center[2])
    #gmsh.model.occ.translate(cell_skeleton[0], center[0], center[1], center[2])

    return u_cell[0][:], skeleton  # Return the tag of the resulting unit cell

# # Create two unit cells
N = 2 
MM = []
MM_skeleton = []
for i in range(N):
    #for j in range(N):
        #for k in range(N):
    l, bone = create_SC_site([i*a, 0.0, 0.0])
    MM.append(l[0][:])
    MM_skeleton += bone
# #l1 = create_SC_site([a, 0.0, 0.0])
# print(MM)
union = gmsh.model.occ.fuse([MM[0]], MM[1:])
#skeleton = gmsh.model.occ.fuse([MM_skeleton[0]], MM_skeleton[1:])
print(MM_skeleton)

# Synchronize the model
gmsh.model.occ.synchronize()
marker = 1
gmsh.model.addPhysicalGroup(3,[union[0][:][0][1]], marker)
gmsh.model.addPhysicalGroup(1,MM_skeleton, marker)

# Save the geometry to a file
gmsh.write(mesh_name+".geo_unrolled")

# Generate the mesh (optional)
gmsh.model.mesh.generate()

# Save the mesh to a file (optional)
gmsh.write(mesh_name+".msh")

# Finalize Gmsh
gmsh.finalize()

msh = meshio.read(mesh_name+".msh")
#Transform mesh to xdmf file. No boundary conditions are necessary so no surface meshes are imported
body_element_tag = "tetra"
line_data = msh.cell_data_dict["gmsh:physical"][body_element_tag]
meshio.write(mesh_name+".xdmf",
    meshio.Mesh(points=msh.points,
        cells={body_element_tag: msh.cells_dict[body_element_tag]},
        cell_data={"bnd_marker": [line_data]}
    )
)

surface_element_tag = "line"
line_data = msh.cell_data_dict["gmsh:physical"][surface_element_tag]
meshio.write(mesh_name+"_skeleton.xdmf",
    meshio.Mesh(points=msh.points,
        cells={surface_element_tag: msh.cells_dict[surface_element_tag]},
        cell_data={"bnd_marker": [line_data]}
    )
)

I also tried to interpolate the scalar function x0^2+x1^2 and it is also wrongly interpolated:

MWE_joint
MWE_joint2500×1406 100 KB

4

Looking more closely at this post, one thing strikes me:
Does the 1D structure align with the edges of the 3D space?

If not, I don’t see how one can expect an accurate interpolation from 1D to 3D, as all fields would have to be extrapolated to match the interpolation points of the 3D structure.
In your code, it seems like they are matching, and they give me a result I would expect with

image
image1467×1247 114 KB

Code:

import numpy as np
import gmsh
import meshio

# Initialize Gmsh
gmsh.initialize()

# Create a new model
gmsh.model.add("model3D")
gmsh.option.setNumber("Geometry.OCCBooleanSimplify", 0)
mesh_name = "pMM2x"
# geom.characteristic_length_max = mesh_size
# Define parameters
r = 0.05  # Small radius
R = 0.25  # Big radius
H = 1.0  # Height
a = 4 * H / np.sqrt(3)
gmsh.option.setNumber("Mesh.CharacteristicLengthMax", 2 * r)


def add_link(ϕ, rotate_axis):
    # Define the three points
    points = [
        [0.0, 0.0, 0.0],
        [r, 0.0, 0.0],  # Point 1
        [R, 0.0, H / 2.0],  # Point 2
        [r, 0.0, H],  # Point 3
        [0.0, 0.0, H],  # Point 4
    ]

    # Define the axis of rotation (e.g., the z-axis)
    axis = [0.0, 0.0, 1.0]
    point_on_axis = [0.0, 0.0, 0.0]

    # Add points to the model
    geom_points = [gmsh.model.occ.addPoint(*p) for p in points]
    sym_axis = gmsh.model.occ.addLine(geom_points[-1], geom_points[0])
    # print(sym_axis)

    # Create lines connecting the points
    lines = []
    for i in range(len(geom_points) - 1):
        lines.append(gmsh.model.occ.addLine(geom_points[i], geom_points[i + 1]))
    lines.append(gmsh.model.occ.addLine(geom_points[-1], geom_points[0]))

    # Create a wire (curve loop) from the lines
    wire = gmsh.model.occ.addCurveLoop(lines)

    # Create a surface from the wire
    surface = gmsh.model.occ.addPlaneSurface([wire])

    # Rotate the surface around the specified axis to create a 3D body
    θ = 2 * np.pi / 3
    s = surface
    bodies = []
    for i in range(3):
        # Revolve the surface
        result = gmsh.model.occ.revolve(
            [(2, s)],  # Surface to revolve
            point_on_axis[0],
            point_on_axis[1],
            point_on_axis[2],  # Point on axis
            axis[0],
            axis[1],
            axis[2],  # Axis of rotation
            θ,  # Angle of rotation
        )
        # Extract the resulting entities
        top = result[0][1]  # The new surface created by the revolution
        body = result[1][1]  # The 3D volume created by the revolution
        bodies.append(body)
        s = top  # Update the surface for the next revolution

    # Add spheres at the ends
    s1 = gmsh.model.occ.addSphere(points[0][0], points[0][1], points[0][2], r)
    s2 = gmsh.model.occ.addSphere(points[-1][0], points[-1][1], points[-1][2], r)

    # Combine all bodies into one
    b = gmsh.model.occ.fuse([(3, s1), (3, s2)], [(3, body) for body in bodies])

    # Rotate the combined body by the specified angle and axis
    gmsh.model.occ.rotate(
        b[0],  # The combined body
        point_on_axis[0],
        point_on_axis[1],
        point_on_axis[2],  # Point on axis
        rotate_axis[0],
        rotate_axis[1],
        rotate_axis[2],  # Axis of rotation
        ϕ,  # Angle of rotation
    )

    gmsh.model.occ.rotate(
        [(1, sym_axis)],  # The combined body
        point_on_axis[0],
        point_on_axis[1],
        point_on_axis[2],  # Point on axis
        rotate_axis[0],
        rotate_axis[1],
        rotate_axis[2],  # Axis of rotation
        ϕ,  # Angle of rotation
    )

    return b[0], sym_axis  # Return the tag of the resulting body


def create_metaatom(translate):
    normal_planes = [
        [1.0, -1.0, 0.0],
        [1.0, -1.0, 0.0],
        [1.0, 1.0, 0.0],
        [1.0, 1.0, 0.0],
    ]
    ϕ = [
        np.arccos(1 / np.sqrt(3)),
        2.0 * np.pi - np.arccos(1 / np.sqrt(3)),
        np.arccos(-1 / np.sqrt(3)),
        2.0 * np.pi - np.arccos(-1 / np.sqrt(3)),
    ]
    links = []
    sk = []

    for i, n in enumerate(normal_planes):
        link, ax = add_link(ϕ[i], n)
        # print(ax)
        gmsh.model.occ.translate([(1, ax)], translate[0], translate[1], translate[2])
        links.append(link[0][:])
        sk.append(ax)

    # Combine all links into one meta-atom
    meta_atom = gmsh.model.occ.fuse([links[0]], links[1:])
    # Translate the meta-atom
    gmsh.model.occ.translate(meta_atom[0], translate[0], translate[1], translate[2])

    return meta_atom[0][:], sk  # Return the tag of the resulting meta-atom


def create_SC_site(center):
    translations = [
        [0.0, 0.0, 0.0],
        [a / 2.0, a / 2.0, 0.0],
        [0.0, a / 2.0, a / 2.0],
        [a / 2.0, 0.0, a / 2.0],
    ]
    unit_cell = []
    skeleton = []

    for t in translations:
        l, sk = create_metaatom(t)
        for link in sk:
            gmsh.model.occ.translate([(1, link)], center[0], center[1], center[2])
        # gmsh.model.occ.translate(sk, center[0], center[1], center[2])
        unit_cell.append(l[0][:])
        skeleton += sk
        # skeleton.append(sk[0][:])

    # Combine all meta-atoms into one unit cell
    u_cell = gmsh.model.occ.fuse([unit_cell[0]], unit_cell[1:])
    # cell_skeleton = gmsh.model.occ.fuse([skeleton[0]], skeleton[1:])

    # Translate the unit cell
    gmsh.model.occ.translate(u_cell[0], center[0], center[1], center[2])
    # gmsh.model.occ.translate(cell_skeleton[0], center[0], center[1], center[2])

    return u_cell[0][:], skeleton  # Return the tag of the resulting unit cell


# # Create two unit cells
N = 2
MM = []
MM_skeleton = []
for i in range(N):
    # for j in range(N):
    # for k in range(N):
    l, bone = create_SC_site([i * a, 0.0, 0.0])
    MM.append(l[0][:])
    MM_skeleton += bone
# #l1 = create_SC_site([a, 0.0, 0.0])
# print(MM)
union = gmsh.model.occ.fuse([MM[0]], MM[1:])
# skeleton = gmsh.model.occ.fuse([MM_skeleton[0]], MM_skeleton[1:])
print(MM_skeleton)

# Synchronize the model
gmsh.model.occ.synchronize()
marker = 1
gmsh.model.addPhysicalGroup(3, [union[0][:][0][1]], marker)
gmsh.model.addPhysicalGroup(1, MM_skeleton, marker)

# Save the geometry to a file
gmsh.write(mesh_name + ".geo_unrolled")

# Generate the mesh (optional)
gmsh.model.mesh.generate()

# Save the mesh to a file (optional)
gmsh.write(mesh_name + ".msh")

# Finalize Gmsh
gmsh.finalize()

msh = meshio.read(mesh_name + ".msh")
# Transform mesh to xdmf file. No boundary conditions are necessary so no surface meshes are imported
body_element_tag = "tetra"
line_data = msh.cell_data_dict["gmsh:physical"][body_element_tag]
meshio.write(
    mesh_name + ".xdmf",
    meshio.Mesh(
        points=msh.points,
        cells={body_element_tag: msh.cells_dict[body_element_tag]},
        cell_data={"bnd_marker": [line_data]},
    ),
)

surface_element_tag = "line"
line_data = msh.cell_data_dict["gmsh:physical"][surface_element_tag]
meshio.write(
    mesh_name + "_skeleton.xdmf",
    meshio.Mesh(
        points=msh.points,
        cells={surface_element_tag: msh.cells_dict[surface_element_tag]},
        cell_data={"bnd_marker": [line_data]},
    ),
)
import dolfinx
import numpy as np
from mpi4py import MPI
from basix.ufl import element


mesh_name = "pMM2x"

# Import Meshes from xdmf file
with dolfinx.io.XDMFFile(MPI.COMM_WORLD, mesh_name + ".xdmf", "r") as xdmf_infile:
    body = xdmf_infile.read_mesh(name="Grid")

with dolfinx.io.XDMFFile(
    MPI.COMM_WORLD, mesh_name + "_skeleton.xdmf", "r"
) as xdmf_infile:
    skeleton = xdmf_infile.read_mesh(name="Grid")


# Define relevant vector elements and spaces
el0 = element("Lagrange", skeleton.basix_cell(), 1)
V0 = dolfinx.fem.functionspace(skeleton, el0)
el1 = element("Lagrange", body.basix_cell(), 1)
V1 = dolfinx.fem.functionspace(body, el1)
padding = 1e-14

t = dolfinx.fem.Function(V0, name="Tangent")
t.interpolate(lambda x: x[0] + 2 * x[1] - x[2])

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "skeleton_tangent.xdmf", "w") as xdmf_outfile:
    xdmf_outfile.write_mesh(skeleton)
    t.name = "Tangent"
    xdmf_outfile.write_function(t)

fine_mesh_cell_map = body.topology.index_map(body.topology.dim)
num_cells_on_proc = fine_mesh_cell_map.size_local + fine_mesh_cell_map.num_ghosts
cells = np.arange(num_cells_on_proc, dtype=np.int32)
interpolation_data = dolfinx.fem.create_interpolation_data(
    V1, V0, cells, padding=padding
)

n = dolfinx.fem.Function(V1, name="Director")
n.interpolate_nonmatching(t, cells, interpolation_data)
n.x.scatter_forward()

with dolfinx.io.XDMFFile(MPI.COMM_WORLD, "body_tangent.xdmf", "w") as xdmf_outfile:
    xdmf_outfile.write_mesh(body)
    n.name = "Tangent"
    xdmf_outfile.write_function(n)

If you turn up the padding, you get something different (padding = 1)

image
image2009×1347 152 KB

as one then extrapolates data from the closest part of the skeleton for each cell of the mesh.

5

I wrote the code so that the links and the beams are created and translated simultaneously, so they should match. Here is a picture of the two superimposed meshes showing that they do in fact match :

match
match944×698 128 KB

6

As stated above, I am not sure what you expect to happen wherever there is no line close to the cell in question. If you allow extrapolation, it gives results across the whole domain

7

Ah sorry, I had misunderstood the question. Yes I want the field to be extrapolated throughout each beam (i.e. “beamwise” constant so to say), so the above result is the right one. This worked beautifully, also for vector interpolation (see image). Thanks @dokken!

MWE_final
MWE_final1013×698 201 KB

As a final question, is there any particular intuition as to how high I should set this padding parameter or is it a trial and error basis?

8

The padding relates to how far one will search around each cell of the 3D mesh for a 1D structure. Thus, if you know the width of the body, that should probably be the max search width.