Merge pull request #120 from nschloe/degen-cells-2d

Fix degenerate cells in 2D

Author: nschloe

setup.cfg

@@ -1,6 +1,6 @@
 [metadata]
 name = pygalmesh
-version = 0.8.1
+version = 0.8.2
 url = https://github.com/nschloe/pygalmesh
 author = Nico Schlömer
 author_email = nico.schloemer@gmail.com

src/generate_2d.cpp

@@ -68,16 +68,22 @@ generate_2d(
   std::map<Vertex_handle, int> vertex_index;
   std::vector<std::array<double, 2>> out_points(cdt.number_of_vertices());
   k = 0;
-  for (auto vit = cdt.vertices_begin(); vit!= cdt.vertices_end(); ++vit) {
+  for (auto vit: cdt.finite_vertex_handles()) {
     out_points[k][0] = vit->point()[0];
     out_points[k][1] = vit->point()[1];
     vertex_index[vit] = k;
     k++;
   }
 
-  std::vector<std::array<int, 3>> out_cells(cdt.number_of_faces());
+  // https://github.com/CGAL/cgal/issues/5068#issuecomment-706213606
+  auto nb_faces = 0u;
+  for (auto fit: cdt.finite_face_handles()) {
+    if(fit->is_in_domain()) ++nb_faces;
+  }
+  std::vector<std::array<int, 3>> out_cells(nb_faces);
   k = 0;
-  for (auto fit = cdt.faces_begin(); fit!= cdt.faces_end(); ++fit) {
+  for (auto fit: cdt.finite_face_handles()) {
+    if(!fit->is_in_domain()) continue;
     out_cells[k][0] = vertex_index[fit->vertex(0)];
     out_cells[k][1] = vertex_index[fit->vertex(1)];
     out_cells[k][2] = vertex_index[fit->vertex(2)];

test/helpers.py

@@ -24,7 +24,13 @@ def compute_triangle_areas(vertices, triangles):
     e0 = vertices[triangles[:, 1]] - vertices[triangles[:, 0]]
     e1 = vertices[triangles[:, 2]] - vertices[triangles[:, 1]]
 
-    e0_cross_e1 = numpy.cross(e0, e1)
-    areas = 0.5 * numpy.sqrt(_row_dot(e0_cross_e1, e0_cross_e1))
-
-    return areas
+    assert e0.shape == e1.shape
+    if e0.shape[1] == 2:
+        z_component_of_e0_cross_e1 = numpy.cross(e0, e1)
+        cross_magnitude = z_component_of_e0_cross_e1
+    else:
+        assert e0.shape[1] == 3
+        e0_cross_e1 = numpy.cross(e0, e1)
+        cross_magnitude = numpy.sqrt(_row_dot(e0_cross_e1, e0_cross_e1))
+
+    return 0.5 * cross_magnitude

test/test_2d.py

@@ -2,8 +2,10 @@
 
 import pygalmesh
 
+from helpers import compute_triangle_areas
 
-def test_2d():
+
+def test_rectangle():
     points = numpy.array([[0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0]])
     constraints = [[0, 1], [1, 2], [2, 3], [3, 0]]
 
@@ -30,5 +32,20 @@ def test_2d():
     # mesh.write("rect.svg")
 
 
+def test_disk():
+    h = 0.1
+    n = int(2 * numpy.pi / h)
+    points = numpy.array(
+        [
+            [numpy.cos(alpha), numpy.sin(alpha)]
+            for alpha in numpy.linspace(0.0, 2 * numpy.pi, n + 1, endpoint=False)
+        ]
+    )
+    constraints = [[k, k + 1] for k in range(n)] + [[n, 0]]
+    mesh = pygalmesh.generate_2d(points, constraints, edge_size=h, num_lloyd_steps=0)
+    areas = compute_triangle_areas(mesh.points, mesh.get_cells_type("triangle"))
+    assert numpy.all(areas > 1.0e-5)
+
+
 if __name__ == "__main__":
-    test_2d()
+    test_disk()

test/test_remesh_surface.py

@@ -52,6 +52,7 @@ def test_inr():
     assert abs(max(mesh.points[:, 2]) - ref[4]) < tol * abs(ref[4])
     assert abs(min(mesh.points[:, 2]) - ref[5]) < tol * abs(ref[5])
 
+    print(helpers.compute_triangle_areas(mesh.points, mesh.get_cells_type("triangle")))
     vol = sum(
         helpers.compute_triangle_areas(mesh.points, mesh.get_cells_type("triangle"))
     )