From 56bf9df2abb70796beb9772b481d09e683caf9c1 Mon Sep 17 00:00:00 2001 From: Huite Bootsma Date: Fri, 30 Aug 2024 13:49:12 +0200 Subject: [PATCH] Add option to shift origin temporarily. Fixes #20 --- docs/api/changelog.rst | 4 ++++ pandamesh/common.py | 14 +++++++++++ pandamesh/gmsh_mesher.py | 17 +++++++++++--- pandamesh/triangle_mesher.py | 20 ++++++++++++++-- tests/test_common.py | 13 +++++++++++ tests/test_meshers.py | 45 ++++++++++++++++++++++++------------ 6 files changed, 93 insertions(+), 20 deletions(-) diff --git a/docs/api/changelog.rst b/docs/api/changelog.rst index 5a1667d..81f2756 100644 --- a/docs/api/changelog.rst +++ b/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 ------------------ diff --git a/pandamesh/common.py b/pandamesh/common.py index a413535..fa66a3e 100644 --- a/pandamesh/common.py +++ b/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 diff --git a/pandamesh/gmsh_mesher.py b/pandamesh/gmsh_mesher.py index add116c..3e88e4b 100644 --- a/pandamesh/gmsh_mesher.py +++ b/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() diff --git a/pandamesh/triangle_mesher.py b/pandamesh/triangle_mesher.py index bc1d0e6..2d7bc6b 100644 --- a/pandamesh/triangle_mesher.py +++ b/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"] diff --git a/tests/test_common.py b/tests/test_common.py index 3dd281d..b030540 100644 --- a/tests/test_common.py +++ b/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]) diff --git a/tests/test_meshers.py b/tests/test_meshers.py index 57d2200..5c8df79 100644 --- a/tests/test_meshers.py +++ b/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():