Merge 56bf9df into f8f93db

docs/api/changelog.rst

@@ -41,6 +41,10 @@ Changed
   :class:`pandamesh.GmshMesher`, :class:`pandamesh.MeshAlgorithm`, or
   :class:`pandamesh.SubdivisionAlgorithm` will now accept one of these enums,
   or the enum member name as a string.
+- :class:`pandamesh.TriangleMesher` and :class:`pandamesh.GmshMesher` now take
+  a ``shift_origin`` argument to temporarily shift the coordinate system to the
+  centroid of the geometries' bounding box to mitigate floating point precision
+  problems. This is enabled by default.
 
 [0.1.6] 2024-07-17
 ------------------

pandamesh/common.py

@@ -212,6 +212,20 @@ def separate(
     return polygons, linestrings, points
 
 
+def central_origin(
+    gdf: gpd.GeoDataFrame, shift_origin: bool
+) -> Tuple[gpd.GeoDataFrame, float, float]:
+    if shift_origin:
+        xmin, ymin, xmax, ymax = gdf.total_bounds
+        xoff = 0.5 * (xmin + xmax)
+        yoff = 0.5 * (ymin + ymax)
+        moved = gdf.copy()
+        moved["geometry"] = gdf["geometry"].translate(xoff=-xoff, yoff=-yoff)
+        return moved, xoff, yoff
+    else:
+        return gdf, 0.0, 0.0
+
+
 def to_ugrid(vertices: FloatArray, faces: IntArray) -> "xugrid.Ugrid2d":  # type: ignore # noqa pragma: no cover
     try:
         import xugrid

pandamesh/gmsh_mesher.py

@@ -11,6 +11,7 @@
 from pandamesh.common import (
     FloatArray,
     IntArray,
+    central_origin,
     check_geodataframe,
     gmsh,
     repr,
@@ -93,7 +94,13 @@ class GmshMesher(MesherBase):
     Parameters
     ----------
     gdf: gpd.GeoDataFrame
-        GeoDataFrame containing the vector geometry.
+        GeoDataFrame containing the vector geometry. Must contain a "cellsize"
+        column.
+    shift_origin: bool, optional, default is True.
+        If True, temporarily shifts the coordinate system origin to the centroid
+        of the geometry's bounding box during mesh generation. This helps mitigate
+        floating-point precision issues. The resulting mesh vertices are
+        automatically translated back to the original coordinate system.
     read_config_files: bool
         Gmsh initialization option: Read system Gmsh configuration files
         (gmshrc and gmsh-options).
@@ -104,14 +111,15 @@ class GmshMesher(MesherBase):
     def __init__(
         self,
         gdf: gpd.GeoDataFrame,
+        shift_origin: bool = True,
         read_config_files: bool = True,
-        run: bool = False,
         interruptible: bool = True,
     ) -> None:
         self._initialize_gmsh(
             read_config_files=read_config_files, interruptible=interruptible
         )
         check_geodataframe(gdf)
+        gdf, self._xoff, self._yoff = central_origin(gdf, shift_origin)
         polygons, linestrings, points = separate(gdf)
 
         # Include geometry into gmsh
@@ -443,7 +451,10 @@ def _vertices(self):
         # getNodes returns: node_tags, coord, parametric_coord
         _, vertices, _ = gmsh.model.mesh.getNodes()
         # Return x and y
-        return vertices.reshape((-1, 3))[:, :2]
+        vertices = vertices.reshape((-1, 3))[:, :2]
+        vertices[:, 0] += self._xoff
+        vertices[:, 1] += self._yoff
+        return np.ascontiguousarray(vertices)
 
     def _faces(self):
         element_types, _, node_tags = gmsh.model.mesh.getElements()

pandamesh/triangle_mesher.py

@@ -6,6 +6,7 @@
 from pandamesh.common import (
     FloatArray,
     IntArray,
+    central_origin,
     check_geodataframe,
     repr,
     separate,
@@ -48,10 +49,22 @@ class TriangleMesher(MesherBase):
 
     For more details on Triangle, see:
     https://www.cs.cmu.edu/~quake/triangle.defs.html
+
+    Parameters
+    ----------
+    gdf: gpd.GeoDataFrame
+        GeoDataFrame containing the vector geometry. Must contain a "cellsize"
+        column.
+    shift_origin: bool, optional, default is True.
+        If True, temporarily shifts the coordinate system origin to the centroid
+        of the geometry's bounding box during mesh generation. This helps mitigate
+        floating-point precision issues. The resulting mesh vertices are
+        automatically translated back to the original coordinate system.
     """
 
-    def __init__(self, gdf: gpd.GeoDataFrame) -> None:
+    def __init__(self, gdf: gpd.GeoDataFrame, shift_origin: bool = True) -> None:
         check_geodataframe(gdf)
+        gdf, self._xoff, self._yoff = central_origin(gdf, shift_origin)
         polygons, linestrings, points = separate(gdf)
         self.vertices, self.segments, self.regions = collect_geometry(
             polygons, linestrings, points
@@ -196,4 +209,7 @@ def generate(self) -> Tuple[FloatArray, IntArray]:
             tri["regions"] = self.regions
 
         result = triangle.triangulate(tri=tri, opts=options)
-        return result["vertices"], result["triangles"]
+        vertices = result["vertices"]
+        vertices[:, 0] += self._xoff
+        vertices[:, 1] += self._yoff
+        return vertices, result["triangles"]

tests/test_common.py

@@ -240,6 +240,19 @@ def test_separate():
         common.separate(gdf)
 
 
+def test_central_origin():
+    gdf = gpd.GeoDataFrame(geometry=[d])
+    back, x, y = common.central_origin(gdf, False)
+    assert gdf is back
+    assert x == 0
+    assert y == 0
+
+    back, x, y = common.central_origin(gdf, True)
+    assert np.allclose(x, 2.5)
+    assert np.allclose(y, 0.5)
+    assert np.array_equal(back.total_bounds, [-0.5, -0.5, 0.5, 0.5])
+
+
 def test_grouper():
     a = np.array([1, 2, 3, 4])
     a_index = np.array([0, 0, 1, 1])

tests/test_meshers.py

@@ -31,6 +31,11 @@
 other_donut = sg.Polygon(outer, holes=[other_hole])
 
 
+def bounds(vertices):
+    x, y = vertices.T
+    return x.min(), y.min(), x.max(), y.max()
+
+
 def area(vertices, triangles):
     """
     Compute the area of every triangle in the mesh.
@@ -42,53 +47,61 @@ def area(vertices, triangles):
     return 0.5 * np.abs(u[:, 0] * v[:, 1] - u[:, 1] * v[:, 0])
 
 
-def triangle_generate(gdf: gpd.GeoDataFrame):
-    mesher = pm.TriangleMesher(gdf)
+def triangle_generate(gdf: gpd.GeoDataFrame, shift: bool):
+    mesher = pm.TriangleMesher(gdf, shift_origin=shift)
     return mesher.generate()
 
 
-def gmsh_generate(gdf: gpd.GeoDataFrame):
-    mesher = pm.GmshMesher(gdf)
+def gmsh_generate(gdf: gpd.GeoDataFrame, shift: bool):
+    mesher = pm.GmshMesher(gdf, shift_origin=shift)
     return mesher.generate()
 
 
 @pytest.mark.parametrize("generate", [triangle_generate, gmsh_generate])
-def test_empty(generate):
+@pytest.mark.parametrize("shift", [False, True])
+def test_empty(generate, shift):
     gdf = gpd.GeoDataFrame(geometry=[line], data={"cellsize": [1.0]})
     with pytest.raises(ValueError, match="No polygons provided"):
-        generate(gdf)
+        generate(gdf, shift)
 
 
 @pytest.mark.parametrize("generate", [triangle_generate, gmsh_generate])
-def test_basic(generate):
+@pytest.mark.parametrize("shift", [False, True])
+def test_basic(generate, shift):
     polygon = sg.Polygon(outer)
     gdf = gpd.GeoDataFrame(geometry=[polygon])
     gdf["cellsize"] = 1.0
-    vertices, triangles = generate(gdf)
+    vertices, triangles = generate(gdf, shift)
     mesh_area = area(vertices, triangles).sum()
     assert np.allclose(mesh_area, polygon.area)
+    assert np.allclose(bounds(vertices), gdf.total_bounds)
 
 
 @pytest.mark.parametrize("generate", [triangle_generate, gmsh_generate])
-def test_hole(generate):
+@pytest.mark.parametrize("shift", [False, True])
+def test_hole(generate, shift):
     gdf = gpd.GeoDataFrame(geometry=[donut])
     gdf["cellsize"] = 1.0
-    vertices, triangles = generate(gdf)
+    vertices, triangles = generate(gdf, shift)
     mesh_area = area(vertices, triangles).sum()
     assert np.allclose(mesh_area, donut.area)
+    assert np.allclose(bounds(vertices), gdf.total_bounds)
 
 
 @pytest.mark.parametrize("generate", [triangle_generate, gmsh_generate])
-def test_partial_hole(generate):
+@pytest.mark.parametrize("shift", [False, True])
+def test_partial_hole(generate, shift):
     gdf = gpd.GeoDataFrame(geometry=[other_donut, other_inner])
     gdf["cellsize"] = 1.0
-    vertices, triangles = generate(gdf)
+    vertices, triangles = generate(gdf, shift)
     mesh_area = area(vertices, triangles).sum()
     assert np.allclose(mesh_area, other_donut.area + other_inner.area)
+    assert np.allclose(bounds(vertices), gdf.total_bounds)
 
 
 @pytest.mark.parametrize("generate", [triangle_generate, gmsh_generate])
-def test_adjacent_donut(generate):
+@pytest.mark.parametrize("shift", [False, True])
+def test_adjacent_donut(generate, shift):
     inner_coords2 = inner_coords.copy()
     outer_coords2 = outer_coords.copy()
     inner_coords2[:, 0] += 10.0
@@ -99,9 +112,10 @@ def test_adjacent_donut(generate):
 
     gdf = gpd.GeoDataFrame(geometry=[donut, donut2])
     gdf["cellsize"] = [1.0, 0.5]
-    vertices, triangles = generate(gdf)
+    vertices, triangles = generate(gdf, shift)
     mesh_area = area(vertices, triangles).sum()
     assert np.allclose(mesh_area, 2 * donut.area)
+    assert np.allclose(bounds(vertices), gdf.total_bounds)
 
     # With a line at y=8.0 and points in the left polygon, at y=2.0
     line1 = sg.LineString([(0.25, 8.0), (9.75, 8.0)])
@@ -112,9 +126,10 @@ def test_adjacent_donut(generate):
     gdf = gpd.GeoDataFrame(geometry=[donut, donut2, line1, line2, *points])
     gdf["cellsize"] = 1.0
 
-    vertices, triangles = generate(gdf)
+    vertices, triangles = generate(gdf, shift)
     mesh_area = area(vertices, triangles).sum()
     assert np.allclose(mesh_area, 2 * donut.area)
+    assert np.allclose(bounds(vertices), gdf.total_bounds)
 
 
 def test_triangle_properties():