Repository: meshpro/dmsh (language: python)

Issue:
Duplicate points
I have found that `dmsh` sometimes generates duplicate points, which makes it difficult to use with meshless methods as this can result in singular matrices. An example is given below which produces three duplicate points. 
```
import dmsh
import optimesh

geo = dmsh.Union([dmsh.Rectangle(-1.0, +1.0,-0.2, +0.2),
                  dmsh.Rectangle(+0.5, +1.0,-0.6, +1.0),
                  dmsh.Circle([-0.5, 0.0], 0.4)])
geo = dmsh.Difference(geo, dmsh.Circle([+0.5, 0.2], 0.2))
r = 0.1   
pts, cells = dmsh.generate(geo, edge_size=r, tol=1.0e-10)
# optionally optimize the mesh
pts, cells = optimesh.optimize_points_cells(pts, cells, "CVT (full)",
                                            1.0e-10, 100)

dmsh.helpers.show(pts, cells, geo)
```
I note there is a helper function `unique_rows()` that may help, but is there a function call that can delete duplicate points and also the associated cells?



Discussion / hints:
Can you try and simplify this more? E.g., is optimesh required for this? How about the tolerance? A simpler mesh?

In any case, duplicate points are bad and should indeed now happen.

_Edit:_ Also, which points are duplicated?
I have tried various values for `r`, setting `tol=1.0e-14` and eliminating `optimesh`, but still get duplicate points.

The following function removes duplicates from the pts array, but that still leaves the associated cells to be removed.
```
def testForDups(pts) :
  # ref: https://stackoverflow.com/questions/31097247/remove-duplicate-rows-of-a-numpy-array
  new_array = [tuple(row) for row in pts]
  pts_unique = np.unique(new_array,axis=0)
  
  return pts_unique
```
And the following function returns the indices of duplicate rows,
```
def findDups(pts) :
  # Ref: https://stackoverflow.com/questions/48099479/how-do-you-find-and-save-duplicated-rows-in-a-numpy-array
  unq, count = np.unique(pts, axis=0, return_counts=True)
  repeated_groups = unq[count > 1]
  
  pts_dupIndices = []
  for repeated_group in repeated_groups:
      repeated_idx = np.argwhere(np.all(pts == repeated_group, axis=1))
      pts_dupIndices.append(repeated_idx.ravel())
  return(np.array(pts_dupIndices))
```
for the above example, it returns,
```
[[ 1 13]
 [ 2 16]]
```
Which indicates that rows `1` and `13` of `pts` are duplicated and rows `2` and `16` are also duplicated.
Further experimentation shows that in the above example, the duplicate points are generated due to the union of two objects which have part of their boundaries exactly overlapping.  By moving the overlapping boundary part of one object so that it lies within the other object boundary, duplicate points are avoided. Therefore, changing the above example to,
```
import dmsh
import optimesh

geo = dmsh.Union([dmsh.Rectangle(-1.0, +0.9,-0.2, +0.2),
                  dmsh.Rectangle(+0.5, +1.0,-0.6, +1.0),
                  dmsh.Circle([-0.5, 0.0], 0.4)])
geo = dmsh.Difference(geo, dmsh.Circle([+0.5, 0.2], 0.2))
r = 0.1   
pts, cells = dmsh.generate(geo, edge_size=r, tol=1.0e-10)
# optionally optimize the mesh
pts, cells = optimesh.optimize_points_cells(pts, cells, "CVT (full)",
                                            1.0e-10, 100)

dmsh.helpers.show(pts, cells, geo)
```
avoids this particular problem. However, this approach is unlikely to solve all problems.
Yeah, I have another example which doesn't have overlapping boundaries but still has duplicate points.  The geometry is that of [DFG benchmark 2D-2](http://www.featflow.de/en/benchmarks/cfdbenchmarking/flow/dfg_benchmark2_re100.html), as reused in the [FEniCS tutorial](https://fenicsproject.org/pub/tutorial/html/._ftut1009.html#___sec78).  Looking to redo that in pure Python, without mshr (gdmcbain/fenics-tuto-in-skfem#5), I tried dmsh but with current main this generates 1670 points of which only 1662 are unique.

```python
import dmsh
import numpy as np

geo = dmsh.Difference(
    dmsh.Rectangle(0.0, 2.2, 0.0, 0.41),
    dmsh.Circle([.2, .2], 0.05)
)

points, triangles = dmsh.generate(geo, 0.025, tol=1e-9)

tmp = np.ascontiguousarray(points)  # test from skfem.Mesh._validate
if points.shape[0] != np.unique(tmp.view([('', tmp.dtype)] * tmp.shape[1])).shape[0]:
    print(points.shape[0], np.unique(tmp.view([('', tmp.dtype)] * tmp.shape[1])).shape[0])
```
I should add that this example was working (no duplicate points) last year; I only noticed it broken a couple of weeks ago.  It still works with
```shell
pip install -U dmsh==0.2.2 meshplex==0.13.6
```
(Not pinning meshplex raises AttributeError: 'MeshTri' object has no attribute 'node_coords'.)

Fix patch:
diff --git a/dmsh/helpers.py b/dmsh/helpers.py
index 2d192f6..c9fbb57 100644
--- a/dmsh/helpers.py
+++ b/dmsh/helpers.py
@@ -32,9 +32,19 @@ def show(pts, cells, geo, title=None, full_screen=True):
     import matplotlib.pyplot as plt
 
     eps = 1.0e-10
+    # hightlight outside points in C3, and points which aren't part of any cell in C4
+    is_part_of_cell = np.zeros(len(pts), dtype=bool)
+    is_part_of_cell[cells.flat] = True
+
     is_inside = geo.dist(pts.T) < eps
-    plt.plot(pts[is_inside, 0], pts[is_inside, 1], ".")
-    plt.plot(pts[~is_inside, 0], pts[~is_inside, 1], ".", color="r")
+
+    sp = pts[is_inside & is_part_of_cell]
+    plt.plot(sp[:, 0], sp[:, 1], ".", color="C0")
+    sp = pts[~is_inside]
+    plt.plot(sp[:, 0], sp[:, 1], ".", color="C3")
+    sp = pts[~is_part_of_cell]
+    plt.plot(sp[:, 0], sp[:, 1], ".", color="k")
+    # plt.plot(pts[~is_inside, 0], pts[~is_part_of_cell, 1], ".", color="k")
     plt.triplot(pts[:, 0], pts[:, 1], cells)
     plt.axis("square")
 
diff --git a/dmsh/main.py b/dmsh/main.py
index f42aff3..b64652a 100644
--- a/dmsh/main.py
+++ b/dmsh/main.py
@@ -36,21 +36,32 @@ def _create_cells(pts, geo):
     return cells
 
 
-def _recell(mesh, geo, flip_tol):
+def _recell_and_boundary_step(mesh, geo, flip_tol):
     # We could do a _create_cells() here, but inverted boundary cell removal plus Lawson
     # flips produce the same result and are much cheaper. This is because, most of the
-    # time, there are no cells to be removed and no edges to be flipped.  Because the
-    # rest of the algorithm is so cheap. The flip is still a fairly expensive operation.
-    # First kick out all boundary cells whose barycenters are not in the geometry.
-    mesh.remove_boundary_cells(
-        lambda is_boundary_cell: geo.dist(mesh.compute_centroids(is_boundary_cell).T)
-        > 0.0
-    )
+    # time, there are no cells to be removed and no edges to be flipped. (The flip is
+    # still a fairly expensive operation.)
+    while True:
+        idx = mesh.is_boundary_point
+        points_new = mesh.points.copy()
+        points_new[idx] = geo.boundary_step(points_new[idx].T).T
+        mesh.points = points_new
+        #
+        num_removed_cells = mesh.remove_boundary_cells(
+            lambda is_bdry_cell: mesh.compute_signed_cell_areas(is_bdry_cell) < 1.0e-10
+        )
+        #
+        # The flip has to come right after the boundary cell removal to prevent
+        # "degenerate cell" errors.
+        mesh.flip_until_delaunay(tol=flip_tol)
+        #
+        if num_removed_cells == 0:
+            break
+
+    # Last kick out all boundary cells whose barycenters are not in the geometry.
     mesh.remove_boundary_cells(
-        lambda is_boundary_cell: mesh.compute_signed_cell_areas(is_boundary_cell)
-        < 1.0e-10
+        lambda is_bdry_cell: geo.dist(mesh.compute_centroids(is_bdry_cell).T) > 0.0
     )
-    mesh.flip_until_delaunay(tol=flip_tol)
 
 
 def create_staggered_grid(h, bounding_box):
@@ -112,7 +123,7 @@ def create_staggered_grid(h, bounding_box):
 
 def generate(
     geo,
-    edge_size: float,
+    edge_size,  # float or function that returns float
     # smoothing_method="distmesh",
     tol: float = 1.0e-5,
     random_seed: int = 0,
@@ -285,10 +296,13 @@ def distmesh_smoothing(
         # # alpha = np.min(max_step / step_lengths)
         # # update *= alpha
 
+        points_old = mesh.points.copy()
+
         # update coordinates
         points_new = mesh.points + update
         # leave feature points untouched
         points_new[:num_feature_points] = mesh.points[:num_feature_points]
+        mesh.points = points_new
         # Some mesh boundary points may have been moved off of the domain boundary,
         # either because they were pushed outside or because they just became boundary
         # points by way of cell removal. Move them all (back) onto the domain boundary.
@@ -297,14 +311,9 @@ def distmesh_smoothing(
         # Alternative: Push all boundary points (the ones _inside_ the geometry as well)
         # back to the boundary.
         # idx = is_outside | is_boundary_point
-        idx = mesh.is_boundary_point
-        points_new[idx] = geo.boundary_step(points_new[idx].T).T
-
-        diff = points_new - mesh.points
-        mesh.points = points_new
-
-        _recell(mesh, geo, flip_tol)
+        _recell_and_boundary_step(mesh, geo, flip_tol)
 
+        diff = points_new - points_old
         move2 = np.einsum("ij,ij->i", diff, diff)
         if verbose:
             print("max_move: {:.6e}".format(np.sqrt(np.max(move2))))
@@ -312,4 +321,8 @@ def distmesh_smoothing(
             break
 
     # print("num steps:  ", k)
+
+    # The cell removal steps in _recell_and_boundary_step() might create points which
+    # aren't part of any cell (dangling points). Remove them now.
+    mesh.remove_dangling_points()
     return mesh
diff --git a/setup.cfg b/setup.cfg
index f1e24d9..a70e7f0 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -1,6 +1,6 @@
 [metadata]
 name = dmsh
-version = 0.2.10
+version = 0.2.11
 author = Nico Schlömer
 author_email = nico.schloemer@gmail.com
 description = High-quality 2D mesh generator based on distmesh
@@ -30,7 +30,7 @@ classifiers =
 packages = find:
 install_requires =
     importlib_metadata;python_version<"3.8"
-    meshplex >= 0.15.5, < 0.16.0
+    meshplex >= 0.15.8, < 0.16.0
     numpy
     scipy
 python_requires = >=3.6


Test patch:
diff --git a/test/test_circle.py b/test/test_circle.py
index 9e77bc9..152dde9 100644
--- a/test/test_circle.py
+++ b/test/test_circle.py
@@ -8,11 +8,14 @@
 
 @pytest.mark.parametrize(
     "radius,ref_norms",
-    [(0.1, [327.95194, 14.263721, 1.0]), (0.4, [18.899253166, 3.70111746, 1.0])],
+    [
+        (0.1, [3.2592107070061820e02, 1.4190745248684369e01, 1.0000000000000000e00]),
+        (0.4, [18.899253166, 3.70111746, 1.0]),
+    ],
 )
 def test_circle(radius, ref_norms, show=False):
     geo = dmsh.Circle([0.0, 0.0], 1.0)
-    X, cells = dmsh.generate(geo, radius, show=show)
+    X, cells = dmsh.generate(geo, radius, show=show, max_steps=100)
     meshplex.MeshTri(X, cells).show()
 
     # make sure the origin is part of the mesh
@@ -28,7 +31,7 @@ def test_circle(radius, ref_norms, show=False):
 )
 def test_degenerate_circle(target_edge_length):
     geo = dmsh.Circle([0.0, 0.0], 1.0)
-    X, cells = dmsh.generate(geo, target_edge_length, show=False)
+    X, cells = dmsh.generate(geo, target_edge_length, show=False, max_steps=200)
 
     mesh = meshplex.MeshTri(X, cells)
     min_q = np.min(mesh.q_radius_ratio)
diff --git a/test/test_difference.py b/test/test_difference.py
index 51a06fe..7f54d18 100644
--- a/test/test_difference.py
+++ b/test/test_difference.py
@@ -6,11 +6,11 @@
 
 def test_difference(show=False):
     geo = dmsh.Difference(dmsh.Circle([-0.5, 0.0], 1.0), dmsh.Circle([+0.5, 0.0], 1.0))
-    X, cells = dmsh.generate(geo, 0.1, show=show)
+    X, cells = dmsh.generate(geo, 0.1, show=show, max_steps=100)
 
     geo.plot()
 
-    ref_norms = [2.9445582949135735e02, 1.5856370081862632e01, 1.4999999025187867e00]
+    ref_norms = [2.9409044729708609e02, 1.5855488859739937e01, 1.5000000000000000e00]
     assert_norm_equality(X.flatten(), ref_norms, 1.0e-6)
     return X, cells
 
diff --git a/test/test_halfspace.py b/test/test_halfspace.py
index 0034fe0..af4bd73 100644
--- a/test/test_halfspace.py
+++ b/test/test_halfspace.py
@@ -11,9 +11,9 @@ def test_halfspace(show=False):
             dmsh.Circle([0.0, 0.0], 1.0),
         ]
     )
-    X, cells = dmsh.generate(geo, 0.1, show=show)
+    X, cells = dmsh.generate(geo, 0.1, show=show, max_steps=100)
 
-    ref_norms = [1.6445956676826719e02, 1.0032819728269299e01, 9.9962353342404042e-01]
+    ref_norms = [1.6399670188761661e02, 1.0011048291798387e01, 9.9959986881486440e-01]
     assert_norm_equality(X.flatten(), ref_norms, 1.0e-6)
     return X, cells
 
diff --git a/test/test_intersection.py b/test/test_intersection.py
index 5b8cf26..e5e4e03 100644
--- a/test/test_intersection.py
+++ b/test/test_intersection.py
@@ -8,11 +8,11 @@ def test_intersection(show=False):
     geo = dmsh.Intersection(
         [dmsh.Circle([0.0, -0.5], 1.0), dmsh.Circle([0.0, +0.5], 1.0)]
     )
-    X, cells = dmsh.generate(geo, 0.1, show=show, tol=1.0e-10)
+    X, cells = dmsh.generate(geo, 0.1, show=show, tol=1.0e-10, max_steps=100)
 
     geo.plot()
 
-    ref_norms = [8.6619344595913475e01, 6.1599895121114274e00, 8.6602540378466342e-01]
+    ref_norms = [8.6491736892894920e01, 6.1568624411912278e00, 8.6602540378466342e-01]
     assert_norm_equality(X.flatten(), ref_norms, 1.0e-10)
     return X, cells
 
@@ -26,9 +26,9 @@ def test_intersection_circles(show=False):
             dmsh.Circle([np.cos(angles[2]), np.sin(angles[2])], 1.5),
         ]
     )
-    X, cells = dmsh.generate(geo, 0.1, show=show, tol=1.0e-10)
+    X, cells = dmsh.generate(geo, 0.1, show=show, tol=1.0e-10, max_steps=100)
 
-    ref_norms = [6.7650924332929776e01, 5.0559320825183374e00, 7.2474487138537913e-01]
+    ref_norms = [6.7661318585210836e01, 5.0568863746561723e00, 7.2474487138537913e-01]
     assert_norm_equality(X.flatten(), ref_norms, 1.0e-10)
     return X, cells
 
diff --git a/test/test_large.py b/test/test_large.py
index ea7109f..6510be9 100644
--- a/test/test_large.py
+++ b/test/test_large.py
@@ -9,9 +9,9 @@ def test_large(show=False):
     c = dmsh.Circle([0.0, 0.0], 3)
     geo = dmsh.Difference(r, c)
 
-    X, cells = dmsh.generate(geo, 2.0, tol=1.0e-5, max_steps=10000, show=show)
+    X, cells = dmsh.generate(geo, 2.0, tol=1.0e-5, max_steps=100, show=show)
 
-    ref_norms = [4.6181926470943245e03, 2.4141659372697231e02, 2.0000000000000000e01]
+    ref_norms = [4.6292581642363657e03, 2.4187329297982635e02, 2.0000000000000000e01]
     assert_norm_equality(X.flatten(), ref_norms, 1.0e-4)
 
 
diff --git a/test/test_pacman.py b/test/test_pacman.py
index 96bd3ad..8e59fd1 100644
--- a/test/test_pacman.py
+++ b/test/test_pacman.py
@@ -8,9 +8,9 @@ def test_pacman(show=False):
         dmsh.Circle([0.0, 0.0], 1.0),
         dmsh.Polygon([[0.0, 0.0], [1.5, 0.4], [1.5, -0.4]]),
     )
-    X, cells = dmsh.generate(geo, 0.1, show=show, tol=1.0e-5)
+    X, cells = dmsh.generate(geo, 0.1, show=show, tol=1.0e-5, max_steps=100)
 
-    ref_norms = [3.0307531341596325e02, 1.3612316633640821e01, 9.9999999953623819e-01]
+    ref_norms = [3.0173012692535394e02, 1.3565685453257570e01, 9.9999999999884770e-01]
     assert_norm_equality(X.flatten(), ref_norms, 1.0e-10)
     return X, cells
 
diff --git a/test/test_polygon.py b/test/test_polygon.py
index aece7a0..f538c7b 100644
--- a/test/test_polygon.py
+++ b/test/test_polygon.py
@@ -17,9 +17,9 @@ def test(show=False):
         ]
     )
     # geo.show()
-    X, cells = dmsh.generate(geo, 0.1, show=show)
+    X, cells = dmsh.generate(geo, 0.1, show=show, max_steps=100)
 
-    ref_norms = [4.1453626885304385e02, 2.1858069585354144e01, 2.0000000000000000e00]
+    ref_norms = [4.1426056822140765e02, 2.1830112296142847e01, 2.0000000000000000e00]
     assert_norm_equality(X.flatten(), ref_norms, 1.0e-5)
     return X, cells
 
diff --git a/test/test_quarter_annulus.py b/test/test_quarter_annulus.py
index adf28f0..40dbbd2 100644
--- a/test/test_quarter_annulus.py
+++ b/test/test_quarter_annulus.py
@@ -14,10 +14,13 @@ def test_quarter_annulus():
     quarter = dmsh.Intersection([diff0, rect])
 
     points, cells = dmsh.generate(
-        quarter, edge_size=lambda x: h + 0.1 * np.abs(disk0.dist(x)), tol=1.0e-10
+        quarter,
+        edge_size=lambda x: h + 0.1 * np.abs(disk0.dist(x)),
+        tol=1.0e-10,
+        max_steps=100,
     )
 
-    ref_norms = [8.0420341247274109e01, 6.7183669347644024e00, 1.0000000000000000e00]
+    ref_norms = [8.0232179592990462e01, 6.6832464479565372e00, 1.0000000000000000e00]
     assert_norm_equality(points.flatten(), ref_norms, 1.0e-10)
     return points, cells
 
diff --git a/test/test_rectangle.py b/test/test_rectangle.py
index cde3e94..0c52586 100644
--- a/test/test_rectangle.py
+++ b/test/test_rectangle.py
@@ -38,13 +38,31 @@ def test_boundary_step():
 
 def test_rectangle(show=False):
     geo = dmsh.Rectangle(-1.0, +2.0, -1.0, +1.0)
-    X, cells = dmsh.generate(geo, 0.1, show=show)
+    X, cells = dmsh.generate(geo, 0.1, show=show, max_steps=100)
 
-    ref_norms = [9.7542855197092831e02, 3.1710489987948261e01, 2.0000000000000000e00]
+    ref_norms = [9.7172325705673779e02, 3.1615286239175994e01, 2.0000000000000000e00]
     assert_norm_equality(X.flatten(), ref_norms, 1.0e-10)
     return X, cells
 
 
+def test_duplicate_points(show=False):
+    # https://github.com/nschloe/dmsh/issues/66
+    # geo = dmsh.Rectangle(0.0, 1.8, 0.0, 0.41)
+
+    # points, cells = dmsh.generate(geo, 0.2, tol=2e-2, show=show)
+
+    # is_part_of_cell = np.zeros(len(points), dtype=bool)
+    # is_part_of_cell[cells.flat] = True
+    # assert np.all(is_part_of_cell)
+
+    geo = dmsh.Rectangle(0.0, 1.4, 0.0, 0.41)
+    points, cells = dmsh.generate(geo, 0.025, tol=1e-5, show=show, max_steps=1)
+    is_part_of_cell = np.zeros(len(points), dtype=bool)
+    is_part_of_cell[cells.flat] = True
+    assert np.all(is_part_of_cell)
+
+
 if __name__ == "__main__":
+    # test_duplicate_points(show=True)
     X, cells = test_rectangle(show=False)
     save("rectangle.png", X, cells)
diff --git a/test/test_rectangle_hole.py b/test/test_rectangle_hole.py
index bd4fc56..b42ca1f 100644
--- a/test/test_rectangle_hole.py
+++ b/test/test_rectangle_hole.py
@@ -8,9 +8,11 @@ def test_rectangle_hole(show=False):
         dmsh.Rectangle(60, 330, 380, 650), dmsh.Rectangle(143, 245, 440, 543)
     )
 
-    X, cells = dmsh.generate(geo, 20, tol=1.0e-5, show=show, flip_tol=1.0e-10)
+    X, cells = dmsh.generate(
+        geo, 20, tol=1.0e-5, show=show, flip_tol=1.0e-10, max_steps=100
+    )
 
-    ref_norms = [1.2895976881893826e05, 7.6183307832433229e03, 6.5000000000000000e02]
+    ref_norms = [1.2931633675576400e05, 7.6377328985582844e03, 6.5000000000000000e02]
     assert_norm_equality(X.flatten(), ref_norms, 1.0e-10)
 
 
@@ -19,9 +21,9 @@ def test_rectangle_hole2(show=False):
         dmsh.Rectangle(0.0, 5.0, 0.0, 5.0),
         dmsh.Polygon([[1, 1], [4, 1], [4, 4], [1, 4]]),
     )
-    X, cells = dmsh.generate(geo, 1.0, show=show, tol=1.0e-3)
+    X, cells = dmsh.generate(geo, 1.0, show=show, tol=1.0e-3, max_steps=100)
 
-    ref_norms = [1.4346826773382367e02, 2.3308127585352914e01, 5.0000000000000000e00]
+    ref_norms = [1.3990406144096474e02, 2.2917592510234346e01, 5.0000000000000000e00]
     assert_norm_equality(X.flatten(), ref_norms, 1.0e-2)
 
 
diff --git a/test/test_refinement_point_line.py b/test/test_refinement_point_line.py
index d089c99..8b821af 100644
--- a/test/test_refinement_point_line.py
+++ b/test/test_refinement_point_line.py
@@ -12,9 +12,9 @@ def test(show=False):
     def edge_size(x):
         return 0.03 + 0.1 * p1.dist(x)
 
-    X, cells = dmsh.generate(geo, edge_size, show=show, tol=1.0e-10)
+    X, cells = dmsh.generate(geo, edge_size, show=show, tol=1.0e-10, max_steps=100)
 
-    ref_norms = [3.8768005962507056e02, 1.5684756178925376e01, 1.0000000000000000e00]
+    ref_norms = [3.8828407463996700e02, 1.5681640918183025e01, 1.0000000000000000e00]
     assert_norm_equality(X.flatten(), ref_norms, 1.0e-3)
     return X, cells
 
diff --git a/test/test_rotation.py b/test/test_rotation.py
index 0a0da07..517de20 100644
--- a/test/test_rotation.py
+++ b/test/test_rotation.py
@@ -6,9 +6,9 @@
 
 def test(show=False):
     geo = dmsh.Rotation(dmsh.Rectangle(-1.0, +2.0, -1.0, +1.0), 0.1 * np.pi)
-    X, cells = dmsh.generate(geo, 0.1, show=show, tol=1.0e-10)
+    X, cells = dmsh.generate(geo, 0.1, show=show, tol=1.0e-10, max_steps=100)
 
-    ref_norms = [9.5352192762855123e02, 3.1344318120249987e01, 2.2111300269652543e00]
+    ref_norms = [9.4730152857365385e02, 3.1160562530932285e01, 2.2111300269652543e00]
     assert_norm_equality(X.flatten(), ref_norms, 1.0e-10)
     return X, cells
 
diff --git a/test/test_scaling.py b/test/test_scaling.py
index 198a21b..3677346 100644
--- a/test/test_scaling.py
+++ b/test/test_scaling.py
@@ -5,9 +5,9 @@
 
 def test(show=False):
     geo = dmsh.Scaling(dmsh.Rectangle(-1.0, +2.0, -1.0, +1.0), 2.0)
-    X, cells = dmsh.generate(geo, 0.1, show=show, tol=1.0e-5)
+    X, cells = dmsh.generate(geo, 0.1, show=show, tol=1.0e-5, max_steps=100)
 
-    ref_norms = [7.7120646314665200e03, 1.2509238764771823e02, 4.0000000000000000e00]
+    ref_norms = [7.6829959173892494e03, 1.2466061090733828e02, 4.0000000000000000e00]
     assert_norm_equality(X.flatten(), ref_norms, 1.0e-7)
     return X, cells
 
diff --git a/test/test_square_hole_refined.py b/test/test_square_hole_refined.py
index 5afc687..8ec7707 100644
--- a/test/test_square_hole_refined.py
+++ b/test/test_square_hole_refined.py
@@ -10,7 +10,11 @@ def test(show=False):
     geo = dmsh.Difference(r, c)
 
     X, cells = dmsh.generate(
-        geo, lambda pts: np.abs(c.dist(pts)) / 5 + 0.05, show=show, tol=1.0e-10
+        geo,
+        lambda pts: np.abs(c.dist(pts)) / 5 + 0.05,
+        show=show,
+        tol=1.0e-10,
+        max_steps=100,
     )
 
     ref_norms = [2.48e02, 1.200e01, 1.0]
diff --git a/test/test_stretch.py b/test/test_stretch.py
index f480f55..39fd549 100644
--- a/test/test_stretch.py
+++ b/test/test_stretch.py
@@ -3,11 +3,11 @@
 import dmsh
 
 
-def test(show=True):
+def test(show=False):
     geo = dmsh.Stretch(dmsh.Rectangle(-1.0, +2.0, -1.0, +1.0), [1.0, 1.0])
-    X, cells = dmsh.generate(geo, 0.3, show=show, tol=1.0e-3)
+    X, cells = dmsh.generate(geo, 0.3, show=show, tol=1.0e-3, max_steps=100)
 
-    ref_norms = [1.9391178579025609e02, 1.5890693098212086e01, 2.6213203435596428e00]
+    ref_norms = [1.9006907971528796e02, 1.5666202908904914e01, 2.6213203435596428e00]
     assert_norm_equality(X.flatten(), ref_norms, 1.0e-2)
     return X, cells
 
diff --git a/test/test_translation.py b/test/test_translation.py
index 1d5f2b5..14974e1 100644
--- a/test/test_translation.py
+++ b/test/test_translation.py
@@ -5,9 +5,9 @@
 
 def test(show=False):
     geo = dmsh.Translation(dmsh.Rectangle(-1.0, +2.0, -1.0, +1.0), [1.0, 1.0])
-    X, cells = dmsh.generate(geo, 0.1, show=show)
+    X, cells = dmsh.generate(geo, 0.1, show=show, max_steps=100)
 
-    ref_norms = [1.7525000000000000e03, 5.5677284113304673e01, 3.0000000000000000e00]
+    ref_norms = [1.7524999999999998e03, 5.5612899955332637e01, 3.0000000000000000e00]
     assert_norm_equality(X.flatten(), ref_norms, 1.0e-7)
 
 
diff --git a/test/test_union.py b/test/test_union.py
index 6145095..7015ad9 100644
--- a/test/test_union.py
+++ b/test/test_union.py
@@ -6,11 +6,11 @@
 
 def test_union_circles(show=False):
     geo = dmsh.Union([dmsh.Circle([-0.5, 0.0], 1.0), dmsh.Circle([+0.5, 0.0], 1.0)])
-    X, cells = dmsh.generate(geo, 0.15, show=show, tol=1.0e-5)
+    X, cells = dmsh.generate(geo, 0.15, show=show, tol=1.0e-5, max_steps=100)
 
     geo.plot()
 
-    ref_norms = [3.0087320808275558e02, 1.5784709390397147e01, 1.5000000000000000e00]
+    ref_norms = [3.0080546580519666e02, 1.5775854476745508e01, 1.5000000000000000e00]
     assert_norm_equality(X.flatten(), ref_norms, 1.0e-10)
     return X, cells
 
@@ -19,9 +19,9 @@ def test_union_rectangles(show=False):
     geo = dmsh.Union(
         [dmsh.Rectangle(-1.0, +0.5, -1.0, +0.5), dmsh.Rectangle(-0.5, +1.0, -0.5, +1.0)]
     )
-    X, cells = dmsh.generate(geo, 0.15, show=show, tol=1.0e-5)
+    X, cells = dmsh.generate(geo, 0.15, show=show, tol=1.0e-5, max_steps=100)
 
-    ref_norms = [1.8420561245539048e02, 1.1269215671323987e01, 1.0000000000000000e00]
+    ref_norms = [1.8417796811774514e02, 1.1277323166424049e01, 1.0000000000000000e00]
     assert_norm_equality(X.flatten(), ref_norms, 1.0e-10)
     return X, cells
 
@@ -35,9 +35,9 @@ def test_union_three_circles(show=False):
             dmsh.Circle([np.cos(angles[2]), np.sin(angles[2])], 1.0),
         ]
     )
-    X, cells = dmsh.generate(geo, 0.2, show=show, tol=1.0e-5)
+    X, cells = dmsh.generate(geo, 0.2, show=show, tol=1.0e-5, max_steps=100)
 
-    ref_norms = [4.0370458355600903e02, 2.1155724692606825e01, 2.0000000000000000e00]
+    ref_norms = [4.0359760255235619e02, 2.1162741423521961e01, 2.0000000000000000e00]
     assert_norm_equality(X.flatten(), ref_norms, 1.0e-10)
     return X, cells