raise MeshplexError("Degenerate cells.")

[rob-rb]

Dear Nico,

I really appreciate your work a lot. I had to adapt the main.py for my workflow. I have one geometry where the workflow now works great (it is for a schwarz-p-surface). now I try to get it work for other minimal surfaces. I have the issue that once in a while the delaunay_flip_code is raising the expection above.

I tried the following fix

mesh.points = new_points
n = len(mesh.cells["points"])
dels = numpy.array(range(n))[mesh.cell_volumes < 1e-2]
mesh.remove_cells(dels)
mesh.flip_until_delaunay()

it removes some cells but the error persists. Could you give me a hint what went wrong?

It's funny actually you write in the readme:

only works for triangular meshes, flat and on a surface

but for me your library seems to work as well for 3d-coords without a surface (Actually the surface it "should" converge to is uknown beforehand)

[rob-rb]

printing the mesh.cell_volumes i get: 0.0004 as a minimum before the breakdown happens. It seems tria_areas in mesh_plex and cell_volumes are not the same. Maybe there are issues with self-interersections i create in the update process

[rob-rb]

yep I tried to reduce max_step when getting this error. The problem copying the mesh by

def copy_mesh(mesh):
    return meshplex.MeshTri(np.array(mesh.points, copy=True), np.array(mesh.cells["points"], copy=True))

leads to the error:

/opt/conda/lib/python3.8/site-packages/meshplex/mesh_tri.py in create_edges(self) 508 a_unique, inv, cts = unique_rows(a) 509 --> 510 assert np.all(cts < 3), "No edge has more than 2 cells. Are cells listed twice?" 511 512 self._is_boundary_edge_local = (cts[inv] == 1).reshape(s[1:])

AssertionError: No edge has more than 2 cells. Are cells listed twice?

sounds like a self-intersection?

[rob-rb]

I tried to spot who is changing the topology of my mesh: it happens calling mesh.flip_until_delaunay(). After the call there is an edge with 4-cells. Funny enough I can make another 4 iterations until the mesh.flip_until_delaunay() generates the title expection. Maybe it would be better to raise an error directly at the topology change? Nico soll ich das versuchen in Meshplex einzufügen?

[nschloe]

It's hard to debug these issues without being able to reproduce them. That's why it's important to always attach an MWE, a minimal working example. This includes the code you used to generate the issue and perhaps other data (a mesh file etc.). In most cases, 50 percent of the work is making the example minimal such that only the essence of the issue remains.

Can you try and produce that?

[rob-rb]

it's gonna be hard work. i'll try

[rob-rb]

The following code fails on the last line. Why?

surface.zip

import meshplex
import numpy as np

from meshplex import *

mesh = meshplex.reader.read("surface.obj")
print("1")
mesh.flip_until_delaunay()
print("done")
points = mesh.points
faces = list(mesh.cells["points"])
mesh = MeshTri(points, faces)
print("done")
mesh.flip_until_delaunay()
print("2")

[nschloe]

For me, this gives

1
done
done
Traceback (most recent call last):
  File "/tmp/mm.py", line 14, in <module>
    mesh.flip_until_delaunay()
  File "/home/nschloe/rcs/meshplex/meshplex/mesh_tri.py", line 1097, in flip_until_delaunay
    if np.all(self.ce_ratios[~self.is_boundary_edge_local] > -0.5 * tol):
  File "/home/nschloe/rcs/meshplex/meshplex/mesh_tri.py", line 429, in is_boundary_edge_local
    self.create_edges()
  File "/home/nschloe/rcs/meshplex/meshplex/mesh_tri.py", line 480, in create_edges
    assert np.all(cts < 3), "No edge has more than 2 cells. Are cells listed twice?"
AssertionError: No edge has more than 2 cells. Are cells listed twice?

[rob-rb]

yep for me too.

and at the same time I expected the following code

points = mesh.points faces = list(mesh.cells["points"]) mesh = MeshTri(points, faces)

to be an identity operator. But it is not (when i introduce a second mesh.flip_until_delaunay() after print("1") it works!). Why?

[nschloe]

yep for me too.

It looks like cells are listed twice then. :) That's usually not what you want.

points = mesh.points
faces = list(mesh.cells["points"])
mesh = MeshTri(points, faces)

to be an identity operator.

What exactly changes?

[rob-rb]

No no the cells are unique:

I introduced this somewhat naive test:

import meshplex
import numpy as np
from meshplex import *

mesh = meshplex.reader.read("surface.obj")
print("1")
mesh.flip_until_delaunay()
print("done")
points = mesh.points
faces0 = list(mesh.cells["points"])

from collections import defaultdict
def control_uniqueness(faces):
    ht = defaultdict(int)
    for f in faces:
        ht[tuple(f)] += 1

    for k, count in ht.items():
        if count>1:
            raise RuntimeError("face is doubled")
    print("faces are unique")
control_uniqueness(faces0)
mesh = MeshTri(points, faces0)
faces1 = list(mesh.cells["points"])

control_uniqueness(faces1)
n = len(faces0)
for j in range(n):
    if tuple(faces1[j])!= tuple(faces0[j]):
        raise RuntimeError("not equal")

mesh.flip_until_delaunay()
print("2")

which has the output

1
done
faces are unique
faces are unique
109027 109027
done
Traceback (most recent call last):
  File "C:/bar/spherene/minimal_optimesh_fail/show_problem.py", line 33, in <module>
    mesh.flip_until_delaunay()
  File "C:\Anaconda\lib\site-packages\meshplex\mesh_tri.py", line 1140, in flip_until_delaunay
    if np.all(self.ce_ratios[~self.is_boundary_edge_local] > -0.5 * tol):
  File "C:\Anaconda\lib\site-packages\meshplex\mesh_tri.py", line 459, in is_boundary_edge_local
    self.create_edges()
  File "C:\Anaconda\lib\site-packages\meshplex\mesh_tri.py", line 510, in create_edges
    assert np.all(cts < 3), "No edge has more than 2 cells. Are cells listed twice?"
AssertionError: No edge has more than 2 cells. Are cells listed twice?

[nschloe]

Note to self: After edge flipping, one cell is listed twice:

import meshplex

mesh = meshplex.read("surface.obj")
mesh.flip_until_delaunay()

print(mesh.cells["points"][208])
print(mesh.cells["points"][222])

[17888 17887 17105]
[17887 17888 17105]

Probably an edge flip bug then.

[rob-rb]

ja stimmt ich hatte einen Fehler in control_uniqueness: jetzt:

def control_uniqueness(faces):
    ht = defaultdict(int)
    for f in faces:
        f = list(f)
        f.sort()
        if f[0]==f[1] or f[1]==f[2]:
            raise RuntimeError("?")
        ht[tuple(f)] += 1
    err = False
    for k, count in ht.items():
        if count > 1:
            err = True
            print("face doubled ",k)
            #raise RuntimeError("face is doubled")
    if not err:
        print("faces are unique")

kriege ich auch face doubled (17105, 17887, 17888)

[nschloe]

Smaller MWE:

import meshplex
import numpy as np

points = np.array([
 [ 0.0, 0.0, 0.0],
 [ 1.0, 0.0, 0.0],
 [ 0.5, 0.3, 0.3],
 [ 0.5, -0.3, 0.3],
])
cells = np.array([
 [1, 2, 0],
 [2, 3, 0],
 [3, 1, 0],
])

mesh = meshplex.MeshTri(points, cells)
print(mesh.ce_ratios_per_interior_facet)
mesh.write("out1.vtk")
mesh.flip_until_delaunay()
mesh.write("out2.vtk")

print(mesh.points)
print(mesh.cells["points"])

mesh.create_facets()

The long edge (bottom) gets flipped and the face on the right-hand side is duplicated.

[nschloe]

The issue is fixed with https://github.com/nschloe/meshplex/pull/119. (Edge flips which would duplicate cells/edges aren't performed.)

[nschloe]

I've now settled for a more consistent fix. The reality is that duplicate cells can happen when performing edge flips on manifold meshes. The data structures are consistent too, but when trying to reinstate the mesh from points and cells only, edge construction will fail. The fix is to call the new remove_duplicate_cells() before constructing edges.