From 6e2c55c6e2c57697bcae2770dc2a7c2d2e0895f1 Mon Sep 17 00:00:00 2001 From: Joseph Capriotti Date: Wed, 30 Sep 2026 11:00:18 -0600 Subject: [PATCH 1/2] Make refine_surface not distort the point cloud to handle horizontal extensions. --- discretize/tree_mesh.py | 71 +++++++++----- discretize/utils/__init__.py | 2 + discretize/utils/mesh_utils.py | 173 +++++++++++++++++++++++++++++++++ tests/tree/test_refine.py | 75 ++++++++++++++ 4 files changed, 295 insertions(+), 26 deletions(-) diff --git a/discretize/tree_mesh.py b/discretize/tree_mesh.py index 0dc7fe92d..bcfa062a7 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( @@ -652,8 +653,11 @@ 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. Horizontally, it does so by adding points + 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). Vertically, the surface is shifted + down (and/or up) by the padding amount. Parameters ---------- @@ -685,6 +689,15 @@ def refine_surface( -------- refine_triangle, refine_vertical_trianglular_prism + Notes + ----- + For surfaces with 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 -------- In 2D we define the surface as a line segment, which we would like to refine @@ -743,46 +756,52 @@ 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) + + # 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: + xyz, simps, pad_dirs = extend_surface_boundary((xyz, simps)) + else: + pad_dirs = np.zeros((n_ps, self.dim - 1)) + 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 + 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..f4fd1552d 100644 --- a/tests/tree/test_refine.py +++ b/tests/tree/test_refine.py @@ -494,6 +494,81 @@ def test_refine_surface3D(): 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) + + 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]]) + + # 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_errors(): mesh = discretize.TreeMesh([32, 32]) points = [[0.3, 0.3], [0.7, 0.3]] From e8c3081b49d94336f97fc04a096d816155ae36cd Mon Sep 17 00:00:00 2001 From: Joseph Capriotti Date: Tue, 6 Oct 2026 22:51:27 -0600 Subject: [PATCH 2/2] restore the old method as a default and make the new way opt in --- discretize/tree_mesh.py | 56 ++++++++++++++++++++++++++++------ tests/tree/test_refine.py | 64 ++++++++++++++++++++++++++++++++++----- 2 files changed, 102 insertions(+), 18 deletions(-) diff --git a/discretize/tree_mesh.py b/discretize/tree_mesh.py index bcfa062a7..d4c3c415b 100644 --- a/discretize/tree_mesh.py +++ b/discretize/tree_mesh.py @@ -645,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. @@ -653,11 +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. Horizontally, it does so by adding points - 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). Vertically, the surface is shifted - down (and/or up) by the padding amount. + 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 ---------- @@ -684,15 +688,22 @@ 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 ----- - For surfaces with a non-convex boundary, the outward bisectors of different - boundary nodes can intersect if the horizontal padding is large compared to + 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 @@ -775,12 +786,33 @@ def refine_surface( 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: + 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: @@ -795,7 +827,11 @@ def refine_surface( 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 - points[:, :-1] = xyz[:, :-1] + pad[:-1] * pad_dirs + 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] diff --git a/tests/tree/test_refine.py b/tests/tree/test_refine.py index f4fd1552d..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,7 +508,14 @@ 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) @@ -501,7 +527,9 @@ def test_refine_surface3D_boundary_nodes(): pad_cells = [[1, 2, 3]] mesh1 = discretize.TreeMesh([32, 32, 32]) - mesh1.refine_surface(points, -1, pad_cells, pad_up=True, pad_down=True) + 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 @@ -522,7 +550,7 @@ def test_refine_surface3D_nonconvex(): 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]]) + 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]])) @@ -569,6 +597,26 @@ def test_extend_surface_boundary(): 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]]