GC_SAFETY_ASSERT FAILURE for manifold mesh

[jo-mueller]

Hi @nmwsharp ,

first of all: Great library.

No to the actual problem: I am trying to use the EdgeFlipGeodesicSolver to find a path on a relatively simple mesh (vertices and faces below) and am receiving an GC_SAFETY_ASSERT FAILURE error like this:

Traceback (most recent call last):
  File "<string>", line 1, in <module>
  File "c:\Users\my_user\AppData\Local\mambaforge\envs\stress\lib\site-packages\potpourri3d\mesh.py", line 64, in __init__
    self.bound_solver = pp3db.EdgeFlipGeodesicsManager(V, F)
RuntimeError: GC_SAFETY_ASSERT FAILURE from D:\a\potpourri3d\potpourri3d\deps\geometry-central\src\surface\manifold_surface_mesh.cpp:110 - duplicate edge in list 1 -- 2

Now, I saw in #18 that this is most likely due to a non-manifold edge. So I did some cleaning (I am using vedo for most mesh operations) with mesh.clean() and also sorted and searched all faces for duplicate triangles. When I check the mesh for whether it is a manifold-mesh using vedo like this, I get:

>>> mesh.is_manifold()
True

Here are the vertices and faces, if it helps to trace down the issue. Here is the mesh in question, in case it helps to debug the issue:

test.zip

Best & thanks for any help,
Johannes

[ElisabethBrockhausQC]

Hi @jo-mueller,

I've had the same problem in the MeshVectorHeatSolver. I ended up manually reordering the mesh. It's just a quick a dirty solution, so not very efficient or anything, but maybe it helps you as well. reorder_edges_of_faces() does the trick for me:

Python code

from sys import setrecursionlimit
from typing import Annotated, Literal, TypeVar

import numpy as np
import numpy.typing as npt

DType = TypeVar("DType", bound=np.generic)
ArrayNx3 = Annotated[npt.NDArray[DType], Literal["N", "3"]]


def _get_duplicate_edges(F: ArrayNx3) -> list[list[int]]:
  """Get list of directed edges that belong to more than one face.

  Parameters
  ----------
  F : ArrayNx3
      Array of vertex indices that form the faces of the tessellation.

  Returns
  -------
  list
      List of duplicate edges.
  """
  edges = np.array([[[a, b], [b, c], [c, a]] for a, b, c in F]).reshape(
      F.shape[0] * 3, 2
  )
  e, c = np.unique(edges, axis=0, return_counts=True)
  return e[c > 1].tolist()


def _reorder_neighbours(
  F: ArrayNx3,
  current_face: int,
  faces_to_neighbours_map: dict,
  reordered_faces: list[int],
) -> ArrayNx3:
  """Recursively reorder neighbouring faces until there are no duplicate edges.

  Parameters
  ----------
  F : ArrayNx3
      Array of vertex indices that form the faces of the tessellation.
  current_face : int
      Index of the correctly ordered reference face whose neighbors are reordered.
  faces_to_neighbours_map : dict
      Dictonary mapping face indices to the indices of all not yet ordered neighbours.
  reordered_faces : list[int]
      List of faces that are already in correct order.

  Returns
  -------
  ArrayNx3
      Array of faces where the neighbours of `current_face` are ordered.
  """
  neighbours = faces_to_neighbours_map[current_face]
  for neighbour in neighbours:
      violation = _get_duplicate_edges(F[[current_face] + [neighbour]])
      if violation:
          F[neighbour] = [F[neighbour][2], F[neighbour][1], F[neighbour][0]]
      reordered_faces.append(neighbour)
      faces_to_neighbours_map[current_face].remove(neighbour)
      faces_to_neighbours_map[neighbour].remove(current_face)

  for i in reordered_faces:
      if faces_to_neighbours_map[i]:
          F = _reorder_neighbours(F, i, faces_to_neighbours_map, reordered_faces)
          break

  if not any(faces_to_neighbours_map.values()):
      return F
  else:
      raise RecursionError("Could not reach all vertices.")


def reorder_edges_of_faces(F: ArrayNx3) -> ArrayNx3:
  """Reorder edges of faces for `potpourri3d.MeshVectorHeatSolver`.

  Parameters
  ----------
  F : ArrayNx3
      Array of unordered faces

  Returns
  -------
  ArrayNx3
      Array of ordered faces
  """
  setrecursionlimit(5 * F.shape[0])

  empty_list: list[int] = []

  edges = np.array([[[a, b], [b, c], [c, a]] for a, b, c in F]).reshape(
      F.shape[0] * 3, 2
  )
  edges_to_faces_map = {(min(e), max(e)): empty_list.copy() for e in edges}
  for i, f in enumerate(F):
      edges_to_faces_map[(min(f[0], f[1]), max(f[0], f[1]))].append(i)
      edges_to_faces_map[(min(f[1], f[2]), max(f[1], f[2]))].append(i)
      edges_to_faces_map[(min(f[2], f[0]), max(f[2], f[0]))].append(i)

  faces_to_neighbours_map = {i: empty_list.copy() for i in range(F.shape[0])}
  for neighbours in edges_to_faces_map.values():
      if len(neighbours) > 1:
          faces_to_neighbours_map[neighbours[0]].append(neighbours[1])
          faces_to_neighbours_map[neighbours[1]].append(neighbours[0])
      if len(neighbours) > 2:
          faces_to_neighbours_map[neighbours[0]].append(neighbours[2])
          faces_to_neighbours_map[neighbours[2]].append(neighbours[0])
          faces_to_neighbours_map[neighbours[1]].append(neighbours[2])
          faces_to_neighbours_map[neighbours[2]].append(neighbours[1])
      if len(neighbours) > 3:
          raise ValueError("Something weird is happening.")

  return _reorder_neighbours(F, 0, faces_to_neighbours_map, [0])