diff --git a/discretize/tree_mesh.py b/discretize/tree_mesh.py index 0dc7fe92d..d4c3c415b 100644 --- a/discretize/tree_mesh.py +++ b/discretize/tree_mesh.py @@ -101,6 +101,7 @@ import scipy.sparse as sp from discretize.utils.code_utils import deprecate_property from scipy.spatial import Delaunay +from discretize.utils.mesh_utils import extend_surface_boundary class TreeMesh( @@ -644,6 +645,7 @@ def refine_surface( pad_down=True, finalize=True, diagonal_balance=None, + extrapolate_surface=None, ): """Refine along a surface triangulated from xyz to the prescribed level. @@ -652,8 +654,14 @@ def refine_surface( level. It also optionally pads the surface at each level based on the number of - padding cells at each dimension. It does so by stretching the bounding box of - the surface in each dimension. + padding cells at each dimension. Vertically, the surface is shifted down + (and/or up) by the padding amount. Horizontally, the padding is done in one + of two ways, set by ``extrapolate_surface``. If ``False`` (the current + default), the bounding box of the surface is stretched in each dimension. If + ``True``, points are added outward from every node along the boundary of the + triangulation, along the outward bisector of the adjacent boundary edges + (equivalent to a nearest neighbor extrapolation of the surface), see + :func:`~discretize.utils.extend_surface_boundary`. Parameters ---------- @@ -680,10 +688,26 @@ def refine_surface( diagonal_balance : None or bool, optional Whether to balance cells diagonally in the refinement, `None` implies using the same setting used to instantiate the TreeMesh`. + extrapolate_surface : None or bool, optional + Whether to pad horizontally by extrapolating the surface outward from the + nodes on its boundary, instead of stretching its bounding box. `None` + implies `False`, but issues a warning if horizontal padding is requested, + as the default will change to `True` in discretize v1.0. This has no + effect if there is no horizontal padding. See Also -------- - refine_triangle, refine_vertical_trianglular_prism + refine_triangle, refine_vertical_trianglular_prism, discretize.utils.extend_surface_boundary + + Notes + ----- + When ``extrapolate_surface=True`` and the surface has a non-convex boundary, + the outward bisectors of different boundary nodes can intersect if the + horizontal padding is large compared to + the size of local features of the boundary (e.g. notches or concavities). In + that case the padded region can fold over itself and the result may be + unreliable. Keep the padding small relative to such features, or use a convex + boundary. Examples -------- @@ -743,46 +767,77 @@ def refine_surface( ): raise ValueError("incorrect dimension for padding_cells_by_level.") + has_h_pad = ( + np.any(padding_cells_by_level[:, :-1]) + if padding_cells_by_level.ndim == 2 + else np.any(padding_cells_by_level) + ) if self.dim == 2: xyz = np.asarray(xyz) sorter = np.argsort(xyz[:, 0]) xyz = xyz[sorter] - n_ps = len(xyz) - inds = np.arange(n_ps) - simps1 = np.c_[inds[:-1], inds[1:], inds[:-1]] + [0, 0, n_ps] - simps2 = np.c_[inds[:-1], inds[1:], inds[1:]] + [n_ps, n_ps, 0] - simps = np.r_[simps1, simps2] + simps = np.c_[np.arange(len(xyz) - 1), np.arange(1, len(xyz))] + elif isinstance(xyz, tuple): + xyz, simps = xyz + xyz = np.asarray(xyz) + simps = np.asarray(simps) else: - if isinstance(xyz, tuple): - xyz, simps = xyz - xyz = np.asarray(xyz) - simps = np.asarray(simps) - else: - xyz = np.asarray(xyz) - triang = Delaunay(xyz[:, :2]) - simps = triang.simplices - n_ps = len(xyz) - - # calculate bounding box for padding - bb_min = np.min(xyz, axis=0)[:-1] - bb_max = np.max(xyz, axis=0)[:-1] - half_width = (bb_max - bb_min) / 2 - center = (bb_max + bb_min) / 2 - points = np.empty((n_ps, self.dim)) - points[:, -1] = xyz[:, -1] + xyz = np.asarray(xyz) + simps = Delaunay(xyz[:, :2]).simplices + n_ps = len(xyz) + + if extrapolate_surface is None: + extrapolate_surface = False + if has_h_pad: + warnings.warn( + "In discretize v1.0 refine_surface will change the default " + "behavior of horizontal padding to extrapolate the surface " + "outwards from the nodes on its boundary, instead of stretching " + "its bounding box. To use the new behavior now, set " + "extrapolate_surface=True. To keep the current behavior, " + "explicitly set extrapolate_surface=False.", + FutureWarning, + stacklevel=2, + ) + + # Add points outward from the boundary of the triangulation, so that padding + # is a scaling of their offsets, without needing to re-triangulate. + if has_h_pad and extrapolate_surface: + xyz, simps, pad_dirs = extend_surface_boundary((xyz, simps)) + else: + pad_dirs = np.zeros((n_ps, self.dim - 1)) + stretch = has_h_pad and not extrapolate_surface + if stretch: + # bounding box of the surface, to stretch for padding + bb_min = np.min(xyz, axis=0)[:-1] + bb_max = np.max(xyz, axis=0)[:-1] + half_width = (bb_max - bb_min) / 2 + center = (bb_max + bb_min) / 2 + n_pts = len(xyz) + + if self.dim == 2: + simps1 = np.c_[simps, simps[:, 0]] + [0, 0, n_pts] + simps2 = np.c_[simps, simps[:, 1]] + [n_pts, n_pts, 0] + simps = np.r_[simps1, simps2] + + points = np.empty((n_pts, self.dim)) h_min = np.r_[[h.min() for h in self.h]] pad = 0.0 for lv, n_pad in zip(np.arange(level, 1, -1), padding_cells_by_level): pad += n_pad * h_min * 2 ** (self.max_level - lv) h = 0 + if stretch: + horizontal_expansion = (half_width + pad[:-1]) / half_width + points[:, :-1] = horizontal_expansion * (xyz[:, :-1] - center) + center + else: + points[:, :-1] = xyz[:, :-1] + pad[:-1] * pad_dirs + points[:, -1] = xyz[:, -1] if pad_up: h += pad[-1] if pad_down: h += pad[-1] points[:, -1] = xyz[:, -1] - pad[-1] - horizontal_expansion = (half_width + pad[:-1]) / half_width - points[:, :-1] = horizontal_expansion * (xyz[:, :-1] - center) + center if self.dim == 2: triangles = np.r_[points, points + [0, h]][simps] self.refine_triangle( diff --git a/discretize/utils/__init__.py b/discretize/utils/__init__.py index 310ed515c..63d2ba30d 100644 --- a/discretize/utils/__init__.py +++ b/discretize/utils/__init__.py @@ -92,6 +92,7 @@ active_from_xyz mesh_builder_xyz slice_mesh + extend_surface_boundary Utilities for Curvilinear Meshes -------------------------------- @@ -139,6 +140,7 @@ mesh_builder_xyz, example_simplex_mesh, slice_mesh, + extend_surface_boundary, ) from discretize.utils.curvilinear_utils import ( example_curvilinear_grid, diff --git a/discretize/utils/mesh_utils.py b/discretize/utils/mesh_utils.py index fa2c7b2fa..2351b5409 100644 --- a/discretize/utils/mesh_utils.py +++ b/discretize/utils/mesh_utils.py @@ -1362,3 +1362,176 @@ def example_simplex_mesh(rect_shape): closest_points = deprecate_function( closest_points_index, "closest_points", removal_version="1.0.0", error=True ) + + +def _outward_boundary_directions(xy, simps, max_miter=5.0): + """Outward offset directions for every boundary node of a 2D triangulation. + + Parameters + ---------- + xy : (n_points, 2) numpy.ndarray + The horizontal locations of the points. + simps : (n_simplices, 3) numpy.ndarray + The triangles of the triangulation. + max_miter : float + Maximum allowed length of the offset direction, (for very sharp corners). + + Returns + ------- + nodes : (n_boundary,) numpy.ndarray + Indices of the nodes on the boundary. + directions : (n_boundary, 2) numpy.ndarray + Outward directions such that offsetting a node by ``d * directions`` moves + the adjacent boundary edges out by a perpendicular distance ``d`` (a miter). + edges : (n_edges, 2) numpy.ndarray + The boundary edges (as node indices). + """ + edges = np.r_[simps[:, [0, 1]], simps[:, [1, 2]], simps[:, [2, 0]]] + third = np.r_[simps[:, 2], simps[:, 0], simps[:, 1]] + _, inverse, counts = np.unique( + np.sort(edges, axis=1), axis=0, return_inverse=True, return_counts=True + ) + on_boundary = counts[inverse.reshape(-1)] == 1 + edges = edges[on_boundary] + third = third[on_boundary] + + tangent = xy[edges[:, 1]] - xy[edges[:, 0]] + normals = np.c_[tangent[:, 1], -tangent[:, 0]] + length = np.linalg.norm(normals, axis=1) + length[length == 0] = 1.0 + normals /= length[:, None] + # point away from the triangle the edge belongs to + inward = np.sum(normals * (xy[third] - xy[edges[:, 0]]), axis=1) > 0 + normals[inward] *= -1 + + nodes = np.unique(edges) + local = np.searchsorted(nodes, edges) + n_sum = np.zeros((len(nodes), 2)) + n_count = np.zeros(len(nodes)) + for k in range(2): + np.add.at(n_sum, local[:, k], normals) + np.add.at(n_count, local[:, k], 1) + + bisector = n_sum / n_count[:, None] + b_len = np.linalg.norm(bisector, axis=1) + b_len[b_len == 0] = 1.0 + bisector /= b_len[:, None] + + # mean of cos(half angle) between the bisector and adjacent edge normals + cos_half = np.zeros(len(nodes)) + for k in range(2): + np.add.at(cos_half, local[:, k], np.sum(bisector[local[:, k]] * normals, 1)) + cos_half /= n_count + scale = 1.0 / np.clip(cos_half, 1.0 / max_miter, 1.0) + return nodes, bisector * scale[:, None], edges + + +def extend_surface_boundary(xyz, max_miter=5.0): + """Extend a triangulated surface outward from every node on its boundary. + + A new point is added for every node along the boundary of the triangulation + (not just the convex hull). Each new point starts at the same location as the + node it comes from, and is offset horizontally along the outward bisector of the + boundary edges adjacent to that node. The new points are triangulated with the + original ones once, so extending the surface by a distance ``d`` (a nearest + neighbor extrapolation of the surface) is just a scaling of the returned + directions. + + Parameters + ---------- + xyz : (n_points, dim) array_like or tuple + Surface points, where ``dim`` is 2 or 3. The last column is the vertical + coordinate. A point cloud is triangulated (Delaunay) along the horizontal + dimensions (in 2D it is sorted along x and connected as a line). Alternatively + supply your own connectivity as a tuple of ``(xyz, simplices)``, where + ``simplices`` are line segments (2D) or triangles (3D) indexing ``xyz``. + max_miter : float, optional + Maximum length of the offset direction, which limits how far the extension + of very sharp boundary corners can reach. + + Returns + ------- + xyz_ext : (n_points + n_new, dim) numpy.ndarray + The original (sorted, in 2D) points followed by the new points, which have + the vertical coordinate of the node they extend. + simplices_ext : (n_simplices + n_new_simplices, dim) numpy.ndarray + The simplices connecting the original and new points. + directions : (n_points + n_new, dim - 1) numpy.ndarray + Horizontal offset directions, which are zero for the original points. The + extended surface for a distance ``d`` (a scalar, or one value per horizontal + dimension) is ``xyz_ext[:, :-1] + d * directions``. The directions are + scaled such that boundary edges move out by a perpendicular distance ``d``. + + Raises + ------ + ValueError + If ``xyz`` has 2 columns and all of its points have the same x location. + + Notes + ----- + For surfaces with a non-convex boundary, the outward bisectors of different + boundary nodes can intersect if the distance is large compared to the size of + local features of the boundary (e.g. notches or concavities). The extended + surface can then fold over itself and be unreliable. + + Examples + -------- + >>> import numpy as np + >>> from discretize.utils import extend_surface_boundary + >>> x, y = np.mgrid[0:1:3j, 0:1:3j] + >>> xyz = np.stack([x, y, x + y], axis=-1).reshape(-1, 3) + >>> xyz_ext, simps, dirs = extend_surface_boundary(xyz) + >>> padded = xyz_ext[:, :-1] + 0.1 * dirs + >>> padded.shape + (17, 2) + """ + if isinstance(xyz, tuple): + xyz, simplices = xyz + simplices = np.asarray(simplices) + else: + simplices = None + xyz = np.asarray(xyz, dtype=float) + dim = xyz.shape[1] + if dim not in (2, 3): + raise ValueError("xyz must have 2 or 3 columns.") + n_ps = len(xyz) + if dim == 2: + if np.ptp(xyz[:, 0]) == 0: + raise ValueError( + "The surface has no horizontal extent, so it has no outward direction." + ) + if simplices is None: + xyz = xyz[np.argsort(xyz[:, 0])] + simplices = np.c_[np.arange(n_ps - 1), np.arange(1, n_ps)] + # boundary nodes are those on a single segment + ends, counts = np.unique(simplices, return_counts=True) + nodes = ends[counts == 1] + # neighbor of each end point, to point away from + neighbor = np.empty(n_ps, dtype=int) + neighbor[simplices[:, 0]] = simplices[:, 1] + neighbor[simplices[:, 1]] = simplices[:, 0] + dirs = np.sign(xyz[nodes, :1] - xyz[neighbor[nodes], :1]) + # vertical segment at an end: fall back to which side of the middle it is on + tie = dirs == 0 + dirs[tie] = np.sign(xyz[nodes, :1] - xyz[:, :1].mean())[tie] + dirs[dirs == 0] = 1.0 + new_inds = n_ps + np.arange(len(nodes)) + simplices = np.r_[simplices, np.c_[nodes, new_inds]] + else: + if simplices is None: + simplices = Delaunay(xyz[:, :2]).simplices + nodes, dirs, edges = _outward_boundary_directions( + xyz[:, :2], simplices, max_miter + ) + new_of_node = np.full(n_ps, -1) + new_of_node[nodes] = n_ps + np.arange(len(nodes)) + i, j = edges[:, 0], edges[:, 1] + simplices = np.r_[ + simplices, + np.c_[i, j, new_of_node[j]], + np.c_[i, new_of_node[j], new_of_node[i]], + ] + xyz_ext = np.r_[xyz, xyz[nodes]] + directions = np.zeros((len(xyz_ext), dim - 1)) + directions[n_ps:] = dirs + return xyz_ext, simplices, directions diff --git a/tests/tree/test_refine.py b/tests/tree/test_refine.py index f8dfc2259..4c110cea4 100644 --- a/tests/tree/test_refine.py +++ b/tests/tree/test_refine.py @@ -1,3 +1,4 @@ +import warnings import re import discretize import numpy as np @@ -439,10 +440,13 @@ def test_refine_points_errors(): mesh1.refine_points(point, 20) -def test_refine_surface2D(): +@pytest.mark.parametrize("extrapolate", [False, True]) +def test_refine_surface2D(extrapolate): mesh1 = discretize.TreeMesh([32, 32]) points = [[0.3, 0.3], [0.7, 0.3]] - mesh1.refine_surface(points, -1, None, pad_up=True, pad_down=True) + mesh1.refine_surface( + points, -1, None, pad_up=True, pad_down=True, extrapolate_surface=extrapolate + ) mesh2 = discretize.TreeMesh([32, 32]) x0 = [0.3, 0.3] @@ -454,7 +458,14 @@ def test_refine_surface2D(): mesh1 = discretize.TreeMesh([32, 32]) points = [[0.3, 0.3], [0.7, 0.3]] n_cell_pad = 2 - mesh1.refine_surface(points, -1, n_cell_pad, pad_up=True, pad_down=True) + mesh1.refine_surface( + points, + -1, + n_cell_pad, + pad_up=True, + pad_down=True, + extrapolate_surface=extrapolate, + ) mesh2 = discretize.TreeMesh([32, 32]) x0 = np.r_[0.3, 0.3] @@ -469,7 +480,8 @@ def test_refine_surface2D(): assert mesh1.equals(mesh2) -def test_refine_surface3D(): +@pytest.mark.parametrize("extrapolate", [False, True]) +def test_refine_surface3D(extrapolate): mesh1 = discretize.TreeMesh([32, 32, 32]) points = [ [0.3, 0.3, 0.5], @@ -477,7 +489,14 @@ def test_refine_surface3D(): [0.3, 0.7, 0.5], [0.7, 0.7, 0.5], ] - mesh1.refine_surface(points, -1, [[1, 2, 3]], pad_up=True, pad_down=True) + mesh1.refine_surface( + points, + -1, + [[1, 2, 3]], + pad_up=True, + pad_down=True, + extrapolate_surface=extrapolate, + ) mesh2 = discretize.TreeMesh([32, 32, 32]) pad = np.array([1, 2, 3]) / 32 @@ -489,11 +508,115 @@ def test_refine_surface3D(): mesh3 = discretize.TreeMesh([32, 32, 32]) simps = [[0, 1, 2], [1, 2, 3]] - mesh3.refine_surface((points, simps), -1, [[1, 2, 3]], pad_up=True, pad_down=True) + mesh3.refine_surface( + (points, simps), + -1, + [[1, 2, 3]], + pad_up=True, + pad_down=True, + extrapolate_surface=extrapolate, + ) assert mesh1.equals(mesh3) +def test_refine_surface3D_boundary_nodes(): + # 3x3 grid, includes mid-edge boundary nodes, should match the padded box + x, y = np.mgrid[0.3:0.7:3j, 0.3:0.7:3j] + points = np.stack([x, y, np.full_like(x, 0.5)], axis=-1).reshape(-1, 3) + pad_cells = [[1, 2, 3]] + + mesh1 = discretize.TreeMesh([32, 32, 32]) + mesh1.refine_surface( + points, -1, pad_cells, pad_up=True, pad_down=True, extrapolate_surface=True + ) + + mesh2 = discretize.TreeMesh([32, 32, 32]) + pad = np.array([1, 2, 3]) / 32 + mesh2.refine_box([0.3, 0.3, 0.5] - pad, [0.7, 0.7, 0.5] + pad, -1) + assert mesh1.equals(mesh2) + + +def test_refine_surface3D_nonconvex(): + # L-shaped surface, padding follows the boundary rather than the bounding box + x, y = np.mgrid[0.2:0.8:13j, 0.2:0.8:13j] + keep = ~((x > 0.5) & (y > 0.5)) + points = np.stack([x[keep], y[keep], np.full(keep.sum(), 0.5)], axis=-1) + # triangulate only the L (drop triangles in the removed corner) + from scipy.spatial import Delaunay + + simps = Delaunay(points[:, :2]).simplices + cent = points[simps, :2].mean(axis=1) + simps = simps[~((cent[:, 0] > 0.5) & (cent[:, 1] > 0.5))] + + mesh1 = discretize.TreeMesh([64, 64, 64]) + mesh1.refine_surface((points, simps), -1, [[2, 2, 2]], extrapolate_surface=True) + + # the removed corner, away from the L, is not refined to the finest level + inds = np.atleast_1d(mesh1.get_containing_cells([[0.75, 0.75, 0.5]])) + assert mesh1.cell_levels_by_index(inds[0]) < mesh1.max_level + # points just outside the L edges are + inds = np.atleast_1d( + mesh1.get_containing_cells([[0.19, 0.3, 0.5], [0.3, 0.19, 0.5]]) + ) + lvl = [mesh1.cell_levels_by_index(i) for i in inds] + assert np.all(np.asarray(lvl) == mesh1.max_level) + + +def test_extend_surface_boundary(): + from discretize.utils import extend_surface_boundary + from scipy.spatial import Delaunay + + x, y = np.mgrid[0.0:1.0:3j, 0.0:1.0:3j] + xyz = np.stack([x, y, x + y], axis=-1).reshape(-1, 3) + xyz_ext, simps, dirs = extend_surface_boundary(xyz) + # every one of the 8 boundary nodes gets a copy + assert len(xyz_ext) == 9 + 8 + assert np.all(dirs[:9] == 0) + # copies keep their parent's elevation + assert np.all(np.isin(xyz_ext[9:, 2], xyz[:, 2])) + padded = xyz_ext[:, :-1] + 0.1 * dirs + assert np.isclose(padded[:, 0].min(), -0.1) + assert np.isclose(padded[:, 0].max(), 1.1) + assert simps.max() == len(xyz_ext) - 1 + + # user supplied triangulation gives the same result + simps_in = Delaunay(xyz[:, :2]).simplices + out2 = extend_surface_boundary((xyz, simps_in)) + out1 = extend_surface_boundary(xyz) + for a, b in zip(out1, out2): + np.testing.assert_allclose(a, b) + + xz = np.array([[0.7, 1.0], [0.3, 2.0]]) + xz_seg = extend_surface_boundary((xz, [[1, 0]])) + np.testing.assert_allclose(xz_seg[0][:2], xz) + np.testing.assert_allclose(xz_seg[2][2:, 0], [1.0, -1.0]) + xz_ext, simps, dirs = extend_surface_boundary(xz) + padded = xz_ext[:, :-1] + 0.1 * dirs + np.testing.assert_allclose(padded[:, 0], [0.3, 0.7, 0.2, 0.8]) + np.testing.assert_allclose(xz_ext[2:, 1], [2.0, 1.0]) + + +def test_refine_surface_extrapolate_warning(): + points = [[0.3, 0.3], [0.7, 0.3]] + + # default is the old behavior, and warns when there is horizontal padding + mesh1 = discretize.TreeMesh([32, 32], diagonal_balance=False) + with pytest.warns(FutureWarning, match="extrapolate_surface"): + mesh1.refine_surface(points, -1, 2) + mesh2 = discretize.TreeMesh([32, 32], diagonal_balance=False) + mesh2.refine_surface(points, -1, 2, extrapolate_surface=False) + assert mesh1.equals(mesh2) + + # no warning without horizontal padding, or when explicitly set + mesh3 = discretize.TreeMesh([32, 32], diagonal_balance=False) + with warnings.catch_warnings(): + warnings.simplefilter("error") + mesh3.refine_surface(points, -1, None) + mesh3.refine_surface(points, -1, [[0, 2]]) + mesh3.refine_surface(points, -1, 2, extrapolate_surface=True) + + def test_refine_surface_errors(): mesh = discretize.TreeMesh([32, 32]) points = [[0.3, 0.3], [0.7, 0.3]]