GJK error in interpolation between non matching second ordered 3D meshes

Lingyue_Shen · October 21, 2024, 9:57am

Hi,
I am trying to do non matching mesh interpolation between 2 different 3D second ordered meshes. I add a slight random shift on the points of one mesh and do nothing on the other mesh. The GJK error always appear. I am using dolfinx version 0.8.0 installed with conda forge. Here is my MWE:

import gmsh
import dolfinx
from dolfinx import cpp
from dolfinx.mesh import GhostMode
from basix.ufl import element
from dolfinx.fem import Function, functionspace, create_nonmatching_meshes_interpolation_data
from dolfinx.io import XDMFFile
from mpi4py import MPI
import numpy as np

def create_mesh(L, W, H, gdim=3, mesh_order=2, mesh_size=4):
    model_rank = 0
    if MPI.COMM_WORLD.rank == model_rank:
        gmsh.model.add('model')
        factory = gmsh.model.occ

        factory.add_rectangle(0, 0, 0, L, W)
        factory.synchronize()
        ov = gmsh.model.get_entities(2)
        factory.extrude(ov, 0, 0, H)
        factory.synchronize()
        ov = gmsh.model.get_entities(gdim)

        ps = gmsh.model.addPhysicalGroup(gdim, [item[1] for item in ov])
        gmsh.model.setPhysicalName(gdim, ps, "fluid")

        gmsh.option.setNumber('Mesh.Algorithm3D', 10)
        gmsh.option.setNumber("General.NumThreads", 16)
        gmsh.option.setNumber("Mesh.MeshSizeMin", mesh_size)
        gmsh.option.setNumber("Mesh.MeshSizeMax", mesh_size)
        gmsh.option.setNumber("General.Verbosity", 2)

        gmsh.model.mesh.generate(gdim)
        gmsh.model.mesh.setOrder(mesh_order)
    
    mesh, _, _ = dolfinx.io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, model_rank, gdim=gdim, partitioner=cpp.mesh.create_cell_partitioner(GhostMode.shared_facet))
    return mesh

def nmm_interpolate(f_out, f_in):
        interpolation_data=create_nonmatching_meshes_interpolation_data(
                                                                        f_out.function_space.mesh,
                                                                        f_out.function_space.element,
                                                                        f_in.function_space.mesh,
                                                                        padding=1e-11)
        f_out.interpolate(f_in,nmm_interpolation_data=interpolation_data)       

def create_function_space(mesh, mesh_order):
    Pv = element("Lagrange", mesh.basix_cell(), mesh_order, shape=(mesh.geometry.dim,))
    V  = functionspace(mesh, Pv)
    return V

def deform_the_mesh(msh, V):
    movement = Function(V)
    movement.x.array[:] += np.random.rand(movement.x.array.size) * 0.1
    msh.geometry.x[:, :msh.geometry.dim] += movement.x.array.reshape((-1, msh.geometry.dim))

if MPI.COMM_WORLD.rank == 0:
    gmsh.initialize()

mesh_order = 2
L, W, H    = 50.0, 50.0, 50.0

for i in range(0, 10):
    # define function on mesh_0
    mesh_0     = create_mesh(L, W, H, gdim=3, mesh_order=mesh_order, mesh_size=4)
    V0         = create_function_space(mesh_0, mesh_order)
    f0         = Function(V0)

    # define function on mesh_new
    mesh_new     = create_mesh(L, W, H, gdim=3, mesh_order=mesh_order, mesh_size=4)
    V_new        = create_function_space(mesh_new, mesh_order)
    f_new        = Function(V_new)
    
    # deform mesh_new with a very slight random shift
    deform_the_mesh(mesh_new, V_new)
    
    # interpolation from mesh_0 to mesh_new and inverse
    nmm_interpolate(f_new, f0)
    if MPI.COMM_WORLD.rank == 0:
        print("nmm_interpolate ok", flush = True)
        
    nmm_interpolate(f0, f_new)    
    if MPI.COMM_WORLD.rank == 0:
        print("inverse nmm_interpolate ok", flush = True)

The error message goes like:

nmm_interpolate ok
Traceback (most recent call last):
  File "/data/nas_rxtd1/Lingyue/iso_parametric/MWE_3.py", line 82, in <module>
    nmm_interpolate(f0, f_new)    
    ^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/data/nas_rxtd1/Lingyue/iso_parametric/MWE_3.py", line 40, in nmm_interpolate
    interpolation_data=create_nonmatching_meshes_interpolation_data(
                       ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/home/shenlingyue/miniconda3/envs/fenicsx080/lib/python3.12/site-packages/dolfinx/fem/__init__.py", line 93, in create_nonmatching_meshes_interpolation_data
    *_create_nonmatching_meshes_interpolation_data(
     ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
RuntimeError: GJK error - max iteration limit reached

What I am confused is that nmm_interpolate(f_new, f0) works just fine. but I can’t do it inversely： nmm_interpolate(f0, f_new)
Could you please explain why could this happen and how to fix the problem?

francesco-ballarin · October 22, 2024, 5:28am

Two questions:

please make sure to set a seed, otherwise each one of us running your code will get a different random mesh
are you sure that after the deformation each cell still has positive volume? Can you add an assert for that?

Lingyue_Shen · October 22, 2024, 8:08am

import gmsh
import dolfinx
from dolfinx import cpp
from dolfinx.mesh import GhostMode
from basix.ufl import element
from dolfinx.fem import form, assemble_vector, Function, functionspace, create_nonmatching_meshes_interpolation_data
from dolfinx.io import XDMFFile
from mpi4py import MPI
import numpy as np
from ufl import TestFunction, dx
def create_mesh(L, W, H, gdim=3, mesh_order=2, mesh_size=4):
    model_rank = 0
    if MPI.COMM_WORLD.rank == model_rank:
        gmsh.model.add('model')
        factory = gmsh.model.occ

        factory.add_rectangle(0, 0, 0, L, W)
        factory.synchronize()
        ov = gmsh.model.get_entities(2)
        factory.extrude(ov, 0, 0, H)
        factory.synchronize()
        ov = gmsh.model.get_entities(gdim)

        ps = gmsh.model.addPhysicalGroup(gdim, [item[1] for item in ov])
        gmsh.model.setPhysicalName(gdim, ps, "fluid")

        gmsh.option.setNumber('Mesh.Algorithm3D', 10)
        gmsh.option.setNumber("General.NumThreads", 16)
        gmsh.option.setNumber("Mesh.MeshSizeMin", mesh_size)
        gmsh.option.setNumber("Mesh.MeshSizeMax", mesh_size)
        gmsh.option.setNumber("General.Verbosity", 2)

        gmsh.model.mesh.generate(gdim)
        gmsh.model.mesh.setOrder(mesh_order)
    
    mesh, _, _ = dolfinx.io.gmshio.model_to_mesh(gmsh.model, MPI.COMM_WORLD, model_rank, gdim=gdim, partitioner=cpp.mesh.create_cell_partitioner(GhostMode.shared_facet))
    return mesh

def nmm_interpolate(f_out, f_in):
        interpolation_data=create_nonmatching_meshes_interpolation_data(
                                                                        f_out.function_space.mesh,
                                                                        f_out.function_space.element,
                                                                        f_in.function_space.mesh,
                                                                        padding=1e-11)
        f_out.interpolate(f_in,nmm_interpolation_data=interpolation_data)       

def create_function_space(mesh, mesh_order):
    Pv = element("Lagrange", mesh.basix_cell(), mesh_order, shape=(mesh.geometry.dim,))
    V  = functionspace(mesh, Pv)
    return V

def deform_the_mesh(msh, V):
    movement = Function(V)
    
    # set seed to np.random
    np.random.seed(0)
    movement.x.array[:] += np.random.rand(movement.x.array.size) * 0.1
    msh.geometry.x[:, :msh.geometry.dim] += movement.x.array.reshape((-1, msh.geometry.dim))
    
    # assure the positve volume
    DG0 = functionspace(msh, ("DG", 0))
    v = TestFunction(DG0)
    cell_area_form = form(v*dx)
    cell_volume = assemble_vector(cell_area_form)
    
    assert np.all(cell_volume.array > 0)

if MPI.COMM_WORLD.rank == 0:
    gmsh.initialize()

mesh_order = 2
L, W, H    = 50.0, 50.0, 50.0

for i in range(0, 10):
    # define function on mesh_0
    mesh_0     = create_mesh(L, W, H, gdim=3, mesh_order=mesh_order, mesh_size=4)
    V0         = create_function_space(mesh_0, mesh_order)
    f0         = Function(V0)

    # define function on mesh_new
    mesh_new     = create_mesh(L, W, H, gdim=3, mesh_order=mesh_order, mesh_size=4)
    V_new        = create_function_space(mesh_new, mesh_order)
    f_new        = Function(V_new)
    
    # deform mesh_new with a very slight random shift
    deform_the_mesh(mesh_new, V_new)
    
    # interpolation from mesh_0 to mesh_new and inverse
    nmm_interpolate(f_new, f0)
    if MPI.COMM_WORLD.rank == 0:
        print("nmm_interpolate ok", flush = True)
        
    nmm_interpolate(f0, f_new)    
    if MPI.COMM_WORLD.rank == 0:
        print("inverse nmm_interpolate ok", flush = True)

I add some code to assure the positive volume and set the seed to the random in deform_the_mesh(msh, V). However the error still shows.

chenyongxin · February 26, 2025, 5:41am

I met the same GJK error when do non-matching interpolation between two 2nd order 3D meshes. See 3D ALE-FSI on a local refined 2nd order mesh: seek advices on remeshing for the my use case. With DOLFINx v0.9.0, only to replace

Lingyue_Shen:

def nmm_interpolate(f_out, f_in):
        interpolation_data=create_nonmatching_meshes_interpolation_data(
                                                                        f_out.function_space.mesh,
                                                                        f_out.function_space.element,
                                                                        f_in.function_space.mesh,
                                                                        padding=1e-11)
        f_out.interpolate(f_in,nmm_interpolation_data=interpolation_data)

by

def nmm_interpolate(f_out, f_in):
        index_map = f_out.function_space.mesh.topology.index_map(3)    
        ncells = index_map.size_local + index_map.num_ghosts
        cells = np.arange(ncells, dtype=np.int32)
        interpolation_data = fem.create_interpolation_data(f_out.function_space, f_in.function_space, 
                                                           cells, padding=1e-11)
        f_out.interpolate_nonmatching(f_in, cells, interpolation_data=interpolation_data)

the error can be reproduced. Any suggestion on this? Or maybe at least how can we know which mesh point cause this error?

dokken · February 26, 2025, 5:48am

Without the actual meshes, it is hard to give you any guidance. What happens if you increase the padding tolerance in interpolation data?

chenyongxin · February 26, 2025, 6:01am

I use code of @Lingyue_Shen to reproduce the same error, which is self-contained.

It is very similar to my use case, except that he shakes every point in the demo as quoted, whereas I keep the boundary points fixed in application. That does not affect the results.

Lingyue_Shen:

# set seed to np.random
    np.random.seed(0)
    movement.x.array[:] += np.random.rand(movement.x.array.size) * 0.1
    msh.geometry.x[:, :msh.geometry.dim] += movement.x.array.reshape((-1, msh.geometry.dim))

I gradually increased padding from 1e-11 to 1, and the error persists.

chenyongxin · February 27, 2025, 12:01pm

Here is meshes.zip - Google Drive which contains two 2nd order 3D mesh files: mesh1 and mesh2 (.xdmf + .h5). The MWE is attached in the end.

mesh1 is directly generated from gmsh and its edges are straight. mesh2 occupies the same space as mesh1, and its interior nodal coordinate is appended by a small disturbance to make its edges “curved”.

f1 and f2 are functions built on mesh1 and mesh2 respectively. The non-matching meshes interpolation f2 <- f1 works (interpolate from a straight-edged mesh to a curved-edged mesh), but the other side interpolation f1 <- f2 fails with RuntimeError: GJK error.

If slightly change the MWE by putting mesh2.xdmf to variable mesh1 (same meshes interpolation), the GJK error still persists. No clue why this happens.

MWE:

# mwe.py with DOLFINx v0.9.0

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

mesh_order = 2

def nonmatching_interpolate(f_out:fem.Function, f_in:fem.Function):
    index_map = f_out.function_space.mesh.topology.index_map(3)    
    ncells = index_map.size_local + index_map.num_ghosts
    cells = np.arange(ncells, dtype=np.int32)
    interpolation_data = fem.create_interpolation_data(f_out.function_space, f_in.function_space, 
                                                       cells, padding=1e-6)
    f_out.interpolate_nonmatching(f_in, cells, interpolation_data=interpolation_data)
    f_out.x.scatter_forward()

with io.XDMFFile(MPI.COMM_WORLD, "mesh1.xdmf", "r") as file:
    mesh1 = file.read_mesh()
    
with io.XDMFFile(MPI.COMM_WORLD, "mesh2.xdmf", "r") as file:
    mesh2 = file.read_mesh()
    
V1 = fem.functionspace(mesh1, ("Lagrange", mesh_order, (3,)))
V2 = fem.functionspace(mesh2, ("Lagrange", mesh_order, (3,)))
f1 = fem.Function(V1)
f2 = fem.Function(V2)

nonmatching_interpolate(f2, f1)
print("f2 <- f1 done", flush=True)

nonmatching_interpolate(f1, f2)
print("f1 <- f2 done", flush=True)

dokken · February 27, 2025, 6:19pm

With the mesh above, I can reduce the problem to:

import numpy as np


p0 = np.array([44.2919, 39.7508, 18.8208, ]).reshape(-1, 3)
p1 = np.array([46.0659, 38.6251, 16.3637, 50.0992, 42.2408, 19.9752, 44.2764, 42.3867, 18.239, 50.0397, 38.4649, 19.9843, 47.1963, 40.4563, 19.1057, 50.0017, 40.354, 19.9947, 47.1886, 42.3193, 19.1362, 48.0266, 38.5593, 18.2475, 45.1735, 40.495, 17.3336, 48.0415, 40.4456, 18.1963, ]).reshape(-1, 3)
import dolfinx
from dolfinx.geometry import compute_distance_gjk
compute_distance_gjk(p0, p1)

which yields this plot:

import matplotlib.pyplot as plt
import matplotlib
fig = plt.figure(figsize=(8, 8))
ax = fig.add_subplot(projection='3d')
ax.scatter(p0[:,0], p0[:,1], p0[:,2], c='r', marker='s')
ax.scatter(p1[:,0], p1[:,1], p1[:,2], c='b', marker='o')
# 
plt.show()

This makes it very weird that the algorithm fail. I’ll debug further.

chenyongxin · March 13, 2025, 9:51am

How is it going so far?

I made up some unmature ideas about this issue and workarounds.

Causes:
The GJK distance algorithim in non-matching meshes interpolation is suitable for convex hulls. In my applications, it involves two 2nd order non-matching meshes interpolation for FSI problems, where the middle points on cell edges is shifted to maintain the curvature of solid boundary. So, a pair of cells crossing the solid boundary must have one convex and one concave shapes, which may cause GJK algorithim fails. I tried to increase the maxk in gjk.h, but the error persists.

Workarounds (two independent ways):

If we can replace the compute_distance_gjk by another suitable method for a more general shape (convex + concave)? It is used in cpp.geometry.determine_point_ownership to determine point-cell relations. For this point, I guess it has a heavy workload as many other functions depend on compute_distance_gjk.
If we can put Functions/DOFs on those 2nd order meshes (curved-edge mesh) to their corresponding 1st order meshes (straight-edge mesh) and then do interpolation via these two 1st order meshes? This could resolve GJK issue as cells with only convex shape are involved. To do so, I planed to do a whole mesh refinement for one time, and then assign DOFs value to the new mesh. But 2nd order mesh does not support mesh refinement.

dokken · March 13, 2025, 9:55am

I have forwarded the bug to @Chris_Richardson .

As you can see from the minimal working example, it is rather strange that it fails for this convex hull, as the point is not even close to the hull.
I have not had time to further investigate after reducing the problem to a “simple” example.

chenyongxin · March 19, 2025, 4:36pm

github.com/FEniCS/dolfinx

[BUG]: GJK fails on simple convex hull

opened 07:01PM - 27 Feb 25 UTC

jorgensd

bug

### Summarize the issue Running the following minimal example makes GJK not con…verge. First reported at: https://fenicsproject.discourse.group/t/gjk-error-in-interpolation-between-non-matching-second-ordered-3d-meshes/16086/8?u=dokken ### How to reproduce the bug See below ### Minimal Example (Python) ```Python import numpy as np p0 = np.array([[44.2919, 39.7508, 18.8208]]) p1 = np.array([[46.0659, 38.6251, 16.3637], [44.2764, 42.3867, 18.239], [ 50.0397, 38.4649, 19.9843], [ 47.1963, 40.4563, 19.1057] , [ 48.0266, 38.5593, 18.2475] , [45.1735, 40.495, 17.3336] , [48.0415, 40.4456, 18.1963 ]]) print(p1.shape) import dolfinx from dolfinx.geometry import compute_distance_gjk compute_distance_gjk(p0, p1) ``` ### Output (Python) ```bash Traceback (most recent call last): File "/home/dokken/Documents/src/debug/mwe.py", line 16, in <module> compute_distance_gjk(p0, p1) ~~~~~~~~~~~~~~~~~~~~^^^^^^^^ File "/home/dokken/src/mambaforge/envs/test_mpl/lib/python3.13/site-packages/dolfinx/geometry.py", line 272, in compute_distance_gjk return _cpp.geometry.compute_distance_gjk(p, q) ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~^^^^^^ RuntimeError: GJK error - max iteration limit reached ``` ### Version main branch ### DOLFINx git commit _No response_ ### Installation Tested both docker and conda, stable and nightly ### Additional information _No response_

From Jorgen’s GJK algorithm failure case (1-to-5 points) as quoted, I’ve got a wild guess. @Chris_Richardson Could you please also comment on this?

The failure is not because of a bug, but the feature of your GJK implementation. Specifically, you only choose the real vertex from convex hulls as your support point, rather than the support point to be interpolated.

I try to visualize the 1-to-5 points failure case with points, corresponding current iteration’s simplex connection (pink dashed line) and sub-simplex connection for the next iteration (blue line) + search direction (red arrow). As you can see that the GJK implementation fails as it falls to a infinite loop since the 3rd iteration.

To compare with the failure case, I slightly move the z value of the 3rd point, and result gets converged at the 3rd iteration with the breaking condition if any existing points are the same as w.

Chris_Richardson · March 20, 2025, 9:43am

If you have a solution, I’m happy to review a Pull Request to dolfinx.

chenyongxin · March 21, 2025, 10:02am

Chris, I have no idea how to improve your implementation so far. But if I completely replace your gjk.h by my implementation with the same interface, it worked on both v0.8.0 and v0.9.0 for the GJK failure case.

I coded up >700 lines of C++, and it definitely lacks of consistency and elegancy as you do in FEniCSx. I’d like to recommend you to have a go with the following link, which I followed in my implementation:

github.com/AlexanderFabisch/distance3d

distance3d/gjk/_gjk_jolt.py

master

"""GJK implementation of Jolt Physics.

Based on: A Fast and Robust GJK Implementation for Collision Detection of
Convex Objects - Gino van den Bergen,
http://www.dtecta.com/papers/jgt98convex.pdf

Copyright 2021 Jorrit Rouwe, MIT license
Source: https://github.com/jrouwe/JoltPhysics/blob/master/Jolt/Geometry/GJKClosestPoint.h
"""
import math
from enum import Enum

import numba
import numpy as np
from ..utils import EPSILON, MAX_FLOAT, scalar_triple_product


EPSILON_SQR = EPSILON * EPSILON
ALL_TRUE = np.array([True, True, True, True], dtype=np.dtype("bool"))

This file has been truncated. show original

chenyongxin · March 21, 2025, 10:03am

Jorgen, thanks for reducing the problem to a minimal 1-to-5 point example. Run the MWE:

import numpy as np

p0 = np.array([[44.2919, 39.7508, 18.8208]])
p1 = np.array(
    [
        [46.0659, 38.6251, 16.3637],
        [44.2764, 42.3867, 18.239],
        [50.0397, 38.4649, 19.9843],
        [45.1735, 40.495, 17.3336],
        [48.0266, 38.5593, 18.2475],
    ]
)

from dolfinx.geometry import compute_distance_gjk
print(compute_distance_gjk(p0, p1))

yielding

[-1.03648825 -1.01316152  1.08986529]

And I can also confirm that the aforementioned two 2nd order meshes interpolation problem is just because of the GJK error. With gjk.h fixed, and now it is resolved.

I will be testing it in real applications to further verify the correctness and robustness.

Lingyue_Shen · March 27, 2025, 5:38am

Thanks for the contribution on this issue greatly!. Now I understand what’s wrong with this function. Hope it could be fixed in future release.

Topic		Replies	Views
Interpolation on non-matching meshes dolfinx	13	1641	February 16, 2023
Dolfinx and pygmsh : error when interpolating dolfinx	1	796	August 21, 2020
Non-matching coordinate ordering for CG1 functions defined on submeshes dolfinx mesh , dolfinx	2	42	November 22, 2024
Create nonmatching meshes interpolation data gives cryptic error. Whats wrong? dolfinx interpolation	3	297	April 16, 2024
Error with 1d mesh immersed in 2d with dolfinx dolfinx	3	355	June 22, 2023

GJK error in interpolation between non matching second ordered 3D meshes

Related topics