From 71eb2f88b41760826edbc02e81d97b627a53e3dc Mon Sep 17 00:00:00 2001 From: Joseph Capriotti Date: Tue, 22 Nov 2022 16:07:56 -0800 Subject: [PATCH 1/9] Add new refine surface, bounding box, and point refine methods. --- discretize/_extensions/tree_ext.pyx | 25 ++- discretize/tree_mesh.py | 305 ++++++++++++++++++++++++++++ discretize/utils/mesh_utils.py | 52 +++++ 3 files changed, 374 insertions(+), 8 deletions(-) diff --git a/discretize/_extensions/tree_ext.pyx b/discretize/_extensions/tree_ext.pyx index d7d617264..ac9bd2bfb 100644 --- a/discretize/_extensions/tree_ext.pyx +++ b/discretize/_extensions/tree_ext.pyx @@ -521,9 +521,9 @@ cdef class _TreeMesh: ---------- points : (N, dim) array_like The centers of the refinement balls - radii : (N) array_like + radii : float or (N) array_like of float A 1D array defining the radius for each ball - levels : (N) array_like of int + levels : int or (N) array_like of int A 1D array defining the maximum refinement level for each ball finalize : bool, optional Whether to finalize after refining @@ -566,12 +566,19 @@ cdef class _TreeMesh: if points.shape[1] != self.dim: raise ValueError(f"points array must be (N, {self.dim})") cdef double[:, :] cs = points - cdef double[:] rs = np.require(np.atleast_1d(radii), dtype=np.float64, + radii = np.require(np.atleast_1d(radii), dtype=np.float64, requirements='C') + if radii.shape[0] == 1: + radii = np.full(points.shape[0], radii[0], dtype=np.float64) + cdef double[:] rs = radii if points.shape[0] != rs.shape[0]: raise ValueError("radii length must match the points array's first dimension") - cdef int[:] ls = np.require(np.atleast_1d(levels), dtype=np.int32, + + levels = np.require(np.atleast_1d(levels), dtype=np.int32, requirements='C') + if levels.shape[0] == 1: + levels = np.full(points.shape[0], levels[0], dtype=np.int32) + cdef int[:] ls = levels if points.shape[0] != ls.shape[0]: raise ValueError("level length must match the points array's first dimension") @@ -603,11 +610,11 @@ cdef class _TreeMesh: The minimum location of the boxes x1s : (N, dim) array_like The maximum location of the boxes - levels : (N) array_like of int + levels : int or (N) array_like of int The level to refine intersecting cells to finalize : bool, optional Whether to finalize after refining - diagonal_balance : bool or None, optional + 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`. @@ -640,7 +647,6 @@ cdef class _TreeMesh: >>> rect = patches.Rectangle([0.8, 0.8], 0.1, 0.2, facecolor='none', edgecolor='k', linewidth=3) >>> ax.add_patch(rect) >>> plt.show() - """ x0s = np.require(np.atleast_2d(x0s), dtype=np.float64, requirements='C') @@ -654,8 +660,11 @@ cdef class _TreeMesh: raise ValueError(f"x0s and x1s must have the same length") cdef double[:, :] x0 = x0s cdef double[:, :] x1 = x1s - cdef int[:] ls = np.require(np.atleast_1d(levels), dtype=np.int32, + levels = np.require(np.atleast_1d(levels), dtype=np.int32, requirements='C') + if levels.shape[0] == 1: + levels = np.full(x0.shape[0], levels[0], dtype=np.int32) + cdef int[:] ls = levels if x0.shape[0] != ls.shape[0]: raise ValueError("level length must match the points array's first dimension") diff --git a/discretize/tree_mesh.py b/discretize/tree_mesh.py index c3d63f8ce..385bda146 100644 --- a/discretize/tree_mesh.py +++ b/discretize/tree_mesh.py @@ -94,6 +94,7 @@ import scipy.sparse as sp import warnings from discretize.utils.code_utils import deprecate_property +from scipy.spatial import Delaunay class TreeMesh( @@ -359,6 +360,310 @@ def origin(self, value): # then update the TreeMesh with the hidden value self._set_origin(self._origin) + def refine_bounding_box(self, points, level, padding_cells_by_level=None, finalize=True, diagonal_balance=None): + """Refine within a bounding box based on the maximum and minimum extent of scattered points. + + This function refines the tree mesh based on the bounding box defined by the + maximum and minimum extent of the scattered input `points`. It will refine + the tree mesh for each level given. + + It also optionally pads the bounding box at each level based on the number of + padding cells at each dimension. + + Parameters + ---------- + points : (N, dim) array_like + The bounding box will be the maximum and minimum extent of these points + level : (n_level) array_like of int + The level of the treemesh to refine the bounding box to. If an array it will + refine the box to each level (mostly useful in combination with `padding_cells_by_level`) + Negative values index backwards, (e.g. `-1` is `max_level`). + padding_cells_by_level : None or (n_level) iterable, optional + The number of cells to pad the bounding box at each level of refinement. + Each entry may be either a single number, which means to pad equally in each + dimension, or it may be an (`dim`) `array_like` for variable padding along + each dimension. None implies no padding. + finalize : bool, optional + Whether to finalize the mesh after the call. + 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`. + + Examples + -------- + Given a set of points, we want to refine the tree mesh with the bounding box + that surrounds those points. The arbitrary points we use for this example are + uniformly scattered between [3/8, 5/8] in the first and second dimension. + + >>> import discretize + >>> import matplotlib.pyplot as plt + >>> import matplotlib.patches as patches + >>> mesh = discretize.TreeMesh([32, 32]) + >>> points = np.random.rand(20, 2) * 0.25 + 3/8 + + Now we want to refine to the maximum level, with a no padding the in `x` + direction and `2` cells in `y`, and the second highest level we want 2 padding + cells in each direction equally beyond that. + + >>> levels = [-1, -2] + >>> padding = [[0, 2], 2] + >>> mesh.refine_bounding_box(points, levels, padding) + + For demonstration we, overlay the bounding box to show the effects of padding. + + >>> ax = mesh.plot_grid() + >>> rect = patches.Rectangle([3/8, 3/8], 1/4, 1/4, facecolor='none', edgecolor='r', linewidth=3) + >>> ax.add_patch(rect) + >>> plt.show() + """ + bsw = np.min(np.atleast_2d(xyz), axis=0) + tnw = np.max(np.atleast_2d(xyz), axis=0) + level = np.atleast_1d(level) + + # pad based on the number of cells at each level + if padding_cells_by_level is None: + padding_cells_by_level = np.zeros_like(level) + padding_cells_by_level = np.asarray(padding_cells_by_level, dtype=object) + padding_cells_by_level = np.atleast_1d(padding_cells_by_level) + if len(level) != len(padding_cells_by_level): + raise ValueError("level and padding_cells_by_level must be the same length") + if np.any(level > self.max_level): + raise IndexError(f"Level beyond max octree level, {self.max_level}") + sorted = np.argsort(level)[::-1] + level = level[sorted] + padding_cells_by_level = padding_cells_by_level[sorted] + + h_min = np.r_[[h.min() for h in self.h]] + x0 = [] + xF = [] + for lv, n_pad in zip(level, padding_cells_by_level): + if lv < 0: + lv = (self.max_level + 1) - (abs(lv) % (self.max_level + 1)) + padding_at_level = n_pad * h_min * 2 ** (self.max_level - lv) + bsw = bsw - padding_at_level + tnw = tnw + padding_at_level + x0.append(bsw) + xF.append(tnw) + + self.refine_box(x0, xF, level, finalize=finalize, diagonal_balance=diagonal_balance) + + def refine_points(self, points, level, padding_cells_by_level=None, finalize=True, diagonal_balance=None): + """Refine the mesh at given points to the prescribed level. + + This function refines the tree mesh around the `points`. It will refine + the tree mesh for each level given. It also optionally radially pads around each + point at each level with the number of padding cells given. + + Parameters + ---------- + points : (N, dim) array_like + The points to be refined around. + level : int or (n_level) array_like of int + The level of the tree mesh to refine to. If an array, it will + refine at each point to each level (mostly useful in combination with + `padding_cells_by_level`). Negative values index backwards, (e.g. `-1` is + `max_level`). + padding_cells_by_level : None or (n_level) iterable, optional + The number of cells each level to pad around the point. The base cell size + is the largest value amongst the smallest cell widths for each dimension. + None implies no padding. + finalize : bool, optional + Whether to finalize the mesh after the call. + 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`. + + Examples + -------- + Given a set of points, we want to refine the tree mesh around these points to + a prescribed level with a certain amount of padding. + + >>> import discretize + >>> import matplotlib.pyplot as plt + >>> import matplotlib.patches as patches + >>> mesh = discretize.TreeMesh([32, 32]) + + Now we want to refine to the maximum level with 1 cell padding around the point, + and then refine at the second highest level with at least 3 cell beyond that. + + >>> points = np.array([[0.1, 0.3], [0.6, 0.8]]) + >>> levels = [-1, -2] + >>> padding = [1, 3] + >>> mesh.refine_points(points, levels, padding) + + >>> ax = mesh.plot_grid() + >>> ax.scatter(*points.T, color='C1') + >>> plt.show() + """ + level = np.atleast_1d(level) + + # pad based on the number of cells at each level + if padding_cells_by_level is None: + padding_cells_by_level = np.zeros_like(level) + padding_cells_by_level = np.atleast_1d(padding_cells_by_level) + if len(level) != len(padding_cells_by_level): + raise ValueError("level and padding_cells_by_level must be the same length") + if np.any(level > self.max_level): + raise IndexError(f"Level beyond max octree level, {self.max_level}") + sorted = np.argsort(level)[::-1] + level = level[sorted] + padding_cells_by_level = padding_cells_by_level[sorted] + + h_min = np.max([h.min() for h in self.h]) + radius_at_level = 0.0 + for lv, n_pad in zip(level, padding_cells_by_level): + if lv < 0: + lv = (self.max_level + 1) - (abs(lv) % (self.max_level + 1)) + radius_at_level += n_pad * h_min * 2 ** (self.max_level - lv) + self.refine_ball( + points, radius_at_level, lv, finalize=False, diagonal_balance=diagonal_balance + ) + + if finalize: + self.finalize() + + def refine_surface(self, xyz, level, padding_cells_by_level=None, pad_up=False, pad_down=True, finalize=True, diagonal_balance=None): + """Refine along a surface triangulated from xyz to the prescribed level. + + This function refines the mesh based on a triangulated surface from the input + points. Every cell that intersects the surface will be refined to the given + level. + + It also optionally pads the bounding box at each level based on the number of + padding cells at each dimension. + + Parameters + ---------- + xyz : (N, dim) array_like + The points defining the surface. Will be triangulated along the horizontal + dimensions. + level : int or (n_level) array_like of int + The level of the tree mesh to refine to. If an array, it will + refine at the surface to each level (mostly useful in combination with + `padding_cells_by_level`). Negative values index backwards, (e.g. `-1` is + `max_level`). + padding_cells_by_level : None or (n_level) iterable, optional + The number of cells to pad the surface with at each level of refinement. + Each entry may be either a single number, which means to pad equally in each + dimension, or it may be an (`dim`) `array_like` for variable padding along + each dimension. None implies no padding. + pad_up : bool, optional + Whether to pad above the surface. + pad_down : bool, optional + Whether to pad below the surface. + finalize : bool, optional + Whether to finalize the mesh after the call. + 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`. + + Examples + -------- + In 2D we define the surface as a line segment, which we would like to refine + along. + + >>> import discretize + >>> import matplotlib.pyplot as plt + >>> import matplotlib.patches as patches + >>> mesh = discretize.TreeMesh([32, 32]) + + This surface is a simple sine curve. Which we would like to pad with at least + 2 cells vertically below at the maximum level, with 3 cells below that at the + second highest level. + + >>> x = np.linspace(0.2, 0.8, 51) + >>> z = 0.25*np.sin(2*np.pi*x)+0.5 + >>> xz = np.c_[x, z] + >>> mesh.refine_surface(xz, [-1, -2], [[0, 2], [0, 3]]) + + >>> ax = mesh.plot_grid() + >>> ax.plot(x, z, color='C1') + >>> plt.show() + + In 3D we define a grid of surface locations with there corresponding elevations. + In this example we pad 2 cells at the finest level below the surface, and 3 + cells down at the next level. + + >>> mesh = discretize.TreeMesh([32, 32, 32]) + >>> x, y = np.mgrid[0.2:0.8:21j, 0.2:0.8:21j] + >>> z = 0.125*np.sin(2*np.pi*x) + 0.5 + 0.125 * np.cos(2 * np.pi * y) + >>> points = np.stack([x, y, z], axis=-1).reshape(-1, 3) + >>> mesh.refine_surface(points, [-1, -2], [[0, 0, 2], [0, 0, 3]]) + + >>> v = mesh.cell_levels_by_index(np.arange(mesh.n_cells)) + >>> fig, axs = plt.subplots(1, 3, figsize=(12,4)) + >>> mesh.plot_slice(v, ax=axs[0], normal='x', grid=True, clim=[2, 5]) + >>> mesh.plot_slice(v, ax=axs[1], normal='y', grid=True, clim=[2, 5]) + >>> mesh.plot_slice(v, ax=axs[2], normal='z', grid=True, clim=[2, 5]) + >>> plt.show() + + """ + level = np.atleast_1d(level) + # pad based on the number of cells at each level + if padding_cells_by_level is None: + padding_cells_by_level = np.zeros_like(level) + padding_cells_by_level = np.asarray(padding_cells_by_level, dtype=object) + padding_cells_by_level = np.atleast_1d(padding_cells_by_level) + if len(level) != len(padding_cells_by_level): + raise ValueError("level and padding_cells_by_level must be the same length") + if np.any(level > self.max_level): + raise IndexError(f"Level beyond max octree level, {self.max_level}") + sorted = np.argsort(level)[::-1] + level = level[sorted] + padding_cells_by_level = padding_cells_by_level[sorted] + + if self.dim == 2: + sorted = np.argsort(xyz[:, 0]) + xyz = xyz[sorted] + 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] + else: + if isinstance(xyz, tuple): + xyz, simps = xyz + else: + triang = Delaunay(xyz[:, :2]) + simps = triang.simplices + n_ps = len(xyz) + simps1 = np.c_[ + simps[:, 0], simps[:, 1], simps[:, 2], simps[:, 0] + ] + [0, 0, 0, n_ps] + simps2 = np.c_[ + simps[:, 0], simps[:, 1], simps[:, 2], simps[:, 1] + ] + [n_ps, n_ps, n_ps, 0] + simps3 = np.c_[ + simps[:, 1], simps[:, 2], simps[:, 0], simps[:, 2] + ] + [0, 0, n_ps, n_ps] + simps = np.r_[simps1, simps2, simps3] + + # 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((2, n_ps, self.dim)) + points[:, :, -1] = xyz[:, -1] + + h_min = np.r_[[h.min() for h in self.h]] + pad = 0.0 + for lv, n_pad in zip(level, padding_cells_by_level): + if lv < 0: + lv = (self.max_level + 1) - (abs(lv) % (self.max_level + 1)) + pad += n_pad * h_min * 2 ** (self.max_level - lv) + if pad_up: + points[0, :, -1] = xyz[:, -1] + pad[-1] + if pad_down: + points[1, :, -1] = xyz[:, -1] - pad[-1] + expansion = (half_width + pad[:-1]) / half_width + points[:, :, :-1] = expansion * (xyz[:, :-1] - center) + center + tetrahedrons = points.reshape((-1, self.dim))[simps] + self.refine_tetrahedron(tetrahedrons, lv, finalize=False, diagonal_balance=diagonal_balance) + + if finalize: + self.finalize() + @property def vntF(self): """ diff --git a/discretize/utils/mesh_utils.py b/discretize/utils/mesh_utils.py index 3a0f0601b..49b73f342 100644 --- a/discretize/utils/mesh_utils.py +++ b/discretize/utils/mesh_utils.py @@ -591,6 +591,14 @@ def refine_tree_xyz( hull, the cells will have a width of *2h* . Within a distance of *nc3 x (4h)* , the cells will have a width of *4h* . Etc... + .. deprecated:: 0.9.0 + `refine_tree_xyz` will be removed in a future version of discretize. It is + replaced by `discretize.TreeMesh.refine_surface`, `discretize.TreeMesh.refine_bounding_box`, + and `discretize.TreeMesh.refine_points`, to separate the calling convetions, + and improve the individual documentation. Those methods are more explicit + about which levels of the TreeMesh you are refining, and provide more + flexibility for padding cells in each dimension. + Parameters ---------- mesh : discretize.TreeMesh @@ -626,6 +634,15 @@ def refine_tree_xyz( discretize.TreeMesh The refined tree mesh + See Also + -------- + discretize.TreeMesh.refine_surface + Recommended to use instead of this function for the `surface` option. + discretize.TreeMesh.refine_bounding_box + Recommended to use instead of this function for the `box` option. + discretize.TreeMesh.refine_points + Recommended to use instead of this function for the `radial` option. + Examples -------- Here we use the **refine_tree_xyz** function refine a tree mesh @@ -693,8 +710,20 @@ def refine_tree_xyz( octree_levels = np.asarray(octree_levels) octree_levels_padding = np.asarray(octree_levels_padding) + # levels = mesh.max_level - np.arange(len(octree_levels)) + # non_zeros = octree_levels != 0 + # levels = levels[non_zeros] + # Trigger different refine methods if method.lower() == "radial": + # padding = octree_levels[non_zeros] + # mesh.refine_points(xyz, levels, padding, finalize=finalize,) + warnings.warn( + DeprecationWarning, + "The radial option is deprecated as of `0.9.0` please update your code to " + "use the `TreeMesh.refine_points` functionality. It will be removed in a " + "future version of discretize." + ) # Compute the outer limits of each octree level rMax = np.cumsum( @@ -712,6 +741,17 @@ def refine_tree_xyz( mesh.finalize() elif method.lower() == "surface": + warnings.warn( + DeprecationWarning, + "The surface option is deprecated as of `0.9.0` please update your code to " + "use the `TreeMesh.refine_surface` functionality. It will be removed in a " + "future version of discretize." + ) + # padding = np.zeros((len(octree_levels), mesh.dim)) + # padding[:, -1] = np.maximum(octree_levels - 1, 0) + # padding[:, :-1] = octree_levels_padding[:, None] + # padding = padding[non_zeros] + # mesh.refine_surface(xyz, levels, padding, finalize=finalize, pad_down=True, pad_up=False) # Compute centroid centroid = np.mean(xyz, axis=0) @@ -823,6 +863,18 @@ def refine_tree_xyz( mesh.finalize() elif method.lower() == "box": + warnings.warn( + DeprecationWarning, + "The box option is deprecated as of `0.9.0` please update your code to " + "use the `TreeMesh.refine_bounding_box` functionality. It will be removed in a " + "future version of discretize." + ) + # padding = np.zeros((len(octree_levels), mesh.dim)) + # padding[:, -1] = np.maximum(octree_levels - 1, 0) + # padding[:, :-1] = octree_levels_padding[:, None] + # padding = padding[non_zeros] + # mesh.refine_bounding_box(xyz, levels, padding, finalize=finalize) + # Define the data extent [bottom SW, top NE] bsw = np.min(xyz, axis=0) tne = np.max(xyz, axis=0) From d2a4fd9078f5de00bded2e559222de44eea52f77 Mon Sep 17 00:00:00 2001 From: Joseph Capriotti Date: Tue, 22 Nov 2022 16:35:41 -0800 Subject: [PATCH 2/9] remove refine_tree_xyz in the examples. --- discretize/_extensions/tree_ext.pyx | 2 -- discretize/tree_mesh.py | 4 +-- tutorials/mesh_generation/4_tree_mesh.py | 39 ++++++++++++------------ 3 files changed, 21 insertions(+), 24 deletions(-) diff --git a/discretize/_extensions/tree_ext.pyx b/discretize/_extensions/tree_ext.pyx index ac9bd2bfb..114eecc31 100644 --- a/discretize/_extensions/tree_ext.pyx +++ b/discretize/_extensions/tree_ext.pyx @@ -2242,7 +2242,6 @@ cdef class _TreeMesh: of the x and y-boundary cells. >>> from discretize import TreeMesh - >>> from discretize.utils import refine_tree_xyz >>> import numpy as np >>> import matplotlib.pyplot as plt @@ -2341,7 +2340,6 @@ cdef class _TreeMesh: of the x and y-boundary faces. >>> from discretize import TreeMesh - >>> from discretize.utils import refine_tree_xyz >>> import numpy as np >>> import matplotlib.pyplot as plt diff --git a/discretize/tree_mesh.py b/discretize/tree_mesh.py index 385bda146..7fb0739e5 100644 --- a/discretize/tree_mesh.py +++ b/discretize/tree_mesh.py @@ -416,8 +416,8 @@ def refine_bounding_box(self, points, level, padding_cells_by_level=None, finali >>> ax.add_patch(rect) >>> plt.show() """ - bsw = np.min(np.atleast_2d(xyz), axis=0) - tnw = np.max(np.atleast_2d(xyz), axis=0) + bsw = np.min(np.atleast_2d(points), axis=0) + tnw = np.max(np.atleast_2d(points), axis=0) level = np.atleast_1d(level) # pad based on the number of cells at each level diff --git a/tutorials/mesh_generation/4_tree_mesh.py b/tutorials/mesh_generation/4_tree_mesh.py index a857a9de9..bf08195ad 100644 --- a/tutorials/mesh_generation/4_tree_mesh.py +++ b/tutorials/mesh_generation/4_tree_mesh.py @@ -20,7 +20,7 @@ - The number of base mesh cells in x, y and z must all be powers of 2 - We cannot refine the mesh to create cells smaller than those defining the base mesh - The range of cell sizes in the tree mesh depends on the number of base mesh cells in x, y and z - + """ @@ -33,7 +33,7 @@ # from discretize import TreeMesh -from discretize.utils import mkvc, refine_tree_xyz +from discretize.utils import mkvc import matplotlib.pyplot as plt import numpy as np @@ -60,10 +60,11 @@ xp, yp = np.meshgrid([120.0, 240.0], [80.0, 160.0]) xy = np.c_[mkvc(xp), mkvc(yp)] # mkvc creates vectors -# Discretize to finest cell size within rectangular box -mesh = refine_tree_xyz(mesh, xy, octree_levels=[2, 2], method="box", finalize=False) - -mesh.finalize() # Must finalize tree mesh before use +# Discretize to finest cell size within rectangular box, with padding in the z direction +# at the finest and second finest levels. +levels = [-1, -2] +padding = [[0, 2], [0, 2]] +mesh.refine_bounding_box(xy, levels, padding) mesh.plotGrid(show_it=True) @@ -101,15 +102,15 @@ xx = mesh.vectorNx yy = -3 * np.exp((xx ** 2) / 100 ** 2) + 50.0 pts = np.c_[mkvc(xx), mkvc(yy)] -mesh = refine_tree_xyz( - mesh, pts, octree_levels=[2, 2], method="surface", finalize=False -) +levels = [-1, -2] +padding = [[0, 2], [0, 2]] +mesh.refine_surface(pts, levels, padding, finalize=False) # Refine mesh near points xx = np.array([0.0, 10.0, 0.0, -10.0]) yy = np.array([-20.0, -10.0, 0.0, -10]) pts = np.c_[mkvc(xx), mkvc(yy)] -mesh = refine_tree_xyz(mesh, pts, octree_levels=[2, 2], method="radial", finalize=False) +mesh.refine_points(pts, [-1, -2], [2, 2], finalize=False) mesh.finalize() @@ -148,15 +149,15 @@ xx = mesh.vectorNx yy = -3 * np.exp((xx ** 2) / 100 ** 2) + 50.0 pts = np.c_[mkvc(xx), mkvc(yy)] -mesh = refine_tree_xyz( - mesh, pts, octree_levels=[2, 2], method="surface", finalize=False -) +levels = [-1 , -2] +padding = [[0, 2], [0, 2]] +mesh.refine_surface(pts, levels, padding, finalize=False) # Refine near points xx = np.array([0.0, 10.0, 0.0, -10.0]) yy = np.array([-20.0, -10.0, 0.0, -10]) pts = np.c_[mkvc(xx), mkvc(yy)] -mesh = refine_tree_xyz(mesh, pts, octree_levels=[2, 2], method="radial", finalize=False) +mesh.refine_points(pts, [-1, -2], [2, 2], finalize=False) mesh.finalize() @@ -213,16 +214,14 @@ [xx, yy] = np.meshgrid(mesh.vectorNx, mesh.vectorNy) zz = -3 * np.exp((xx ** 2 + yy ** 2) / 100 ** 2) + 50.0 pts = np.c_[mkvc(xx), mkvc(yy), mkvc(zz)] -mesh = refine_tree_xyz( - mesh, pts, octree_levels=[2, 2], method="surface", finalize=False -) +levels = [-1, -2] +padding = [[0, 0, 2], [0, 0, 2]] +mesh.refine_surface(pts, levels, padding, finalize=False) # Refine box xp, yp, zp = np.meshgrid([-40.0, 40.0], [-40.0, 40.0], [-60.0, 0.0]) xyz = np.c_[mkvc(xp), mkvc(yp), mkvc(zp)] - -mesh = refine_tree_xyz(mesh, xyz, octree_levels=[2, 2], method="box", finalize=False) - +mesh.refine_bounding_box(xyz, levels, padding, finalize=False) mesh.finalize() # The bottom west corner From 8f8a8e56c5200b15c635a3015bbadf6aa6d24e08 Mon Sep 17 00:00:00 2001 From: Joseph Capriotti Date: Tue, 22 Nov 2022 16:59:16 -0800 Subject: [PATCH 3/9] fix deprecation warning --- discretize/utils/mesh_utils.py | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/discretize/utils/mesh_utils.py b/discretize/utils/mesh_utils.py index 49b73f342..16fffe59b 100644 --- a/discretize/utils/mesh_utils.py +++ b/discretize/utils/mesh_utils.py @@ -719,10 +719,10 @@ def refine_tree_xyz( # padding = octree_levels[non_zeros] # mesh.refine_points(xyz, levels, padding, finalize=finalize,) warnings.warn( - DeprecationWarning, "The radial option is deprecated as of `0.9.0` please update your code to " "use the `TreeMesh.refine_points` functionality. It will be removed in a " "future version of discretize." + DeprecationWarning, ) # Compute the outer limits of each octree level @@ -742,10 +742,10 @@ def refine_tree_xyz( elif method.lower() == "surface": warnings.warn( - DeprecationWarning, "The surface option is deprecated as of `0.9.0` please update your code to " "use the `TreeMesh.refine_surface` functionality. It will be removed in a " - "future version of discretize." + "future version of discretize.", + DeprecationWarning, ) # padding = np.zeros((len(octree_levels), mesh.dim)) # padding[:, -1] = np.maximum(octree_levels - 1, 0) @@ -864,10 +864,10 @@ def refine_tree_xyz( elif method.lower() == "box": warnings.warn( - DeprecationWarning, "The box option is deprecated as of `0.9.0` please update your code to " "use the `TreeMesh.refine_bounding_box` functionality. It will be removed in a " - "future version of discretize." + "future version of discretize.", + DeprecationWarning, ) # padding = np.zeros((len(octree_levels), mesh.dim)) # padding[:, -1] = np.maximum(octree_levels - 1, 0) From 64a699c45f07dbc60d446c87780562fa82fe8aa2 Mon Sep 17 00:00:00 2001 From: Joseph Capriotti Date: Wed, 23 Nov 2022 00:00:29 -0800 Subject: [PATCH 4/9] fix syntax --- discretize/utils/mesh_utils.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/discretize/utils/mesh_utils.py b/discretize/utils/mesh_utils.py index 16fffe59b..ef4dd11ee 100644 --- a/discretize/utils/mesh_utils.py +++ b/discretize/utils/mesh_utils.py @@ -721,7 +721,7 @@ def refine_tree_xyz( warnings.warn( "The radial option is deprecated as of `0.9.0` please update your code to " "use the `TreeMesh.refine_points` functionality. It will be removed in a " - "future version of discretize." + "future version of discretize.", DeprecationWarning, ) From c7e5bfab5d32b07dceb94556b83f805dbb2af6ca Mon Sep 17 00:00:00 2001 From: Joseph Capriotti Date: Wed, 23 Nov 2022 09:43:13 -0800 Subject: [PATCH 5/9] fix working error test --- tests/tree/test_refine.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/tree/test_refine.py b/tests/tree/test_refine.py index 5f5718b65..ed600972d 100644 --- a/tests/tree/test_refine.py +++ b/tests/tree/test_refine.py @@ -367,7 +367,7 @@ def test_box_errors(): # incorrect number of levels with pytest.raises(ValueError): - mesh.refine_box(x0s2d, x1s2d, [mesh.max_level], finalize=False) + mesh.refine_box(x0s2d, x1s2d, [-1, -1, -1], finalize=False) def test_ball_errors(): From 7f3a95612bd188ff59bdf70800fd98e6cb1e9e75 Mon Sep 17 00:00:00 2001 From: Joseph Capriotti Date: Wed, 23 Nov 2022 12:28:35 -0800 Subject: [PATCH 6/9] add tests for new functions --- discretize/tree_mesh.py | 18 ++++--- tests/tree/test_refine.py | 103 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 115 insertions(+), 6 deletions(-) diff --git a/discretize/tree_mesh.py b/discretize/tree_mesh.py index 7fb0739e5..77bbf0ae0 100644 --- a/discretize/tree_mesh.py +++ b/discretize/tree_mesh.py @@ -403,7 +403,7 @@ def refine_bounding_box(self, points, level, padding_cells_by_level=None, finali Now we want to refine to the maximum level, with a no padding the in `x` direction and `2` cells in `y`, and the second highest level we want 2 padding - cells in each direction equally beyond that. + cells in each direction beyond that. >>> levels = [-1, -2] >>> padding = [[0, 2], 2] @@ -464,7 +464,7 @@ def refine_points(self, points, level, padding_cells_by_level=None, finalize=Tru `padding_cells_by_level`). Negative values index backwards, (e.g. `-1` is `max_level`). padding_cells_by_level : None or (n_level) iterable, optional - The number of cells each level to pad around the point. The base cell size + The number of cells at each level to pad around the point. The base cell size is the largest value amongst the smallest cell widths for each dimension. None implies no padding. finalize : bool, optional @@ -529,14 +529,16 @@ def refine_surface(self, xyz, level, padding_cells_by_level=None, pad_up=False, points. Every cell that intersects the surface will be refined to the given level. - It also optionally pads the bounding box at each level based on the number of - padding cells at each dimension. + 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. Parameters ---------- - xyz : (N, dim) array_like + xyz : (N, dim) array_like or tuple The points defining the surface. Will be triangulated along the horizontal - dimensions. + dimensions. You are able to supply your own triangulation in 3D by inputing + a tuple of (`xyz`, `triangle_indices`). level : int or (n_level) array_like of int The level of the tree mesh to refine to. If an array, it will refine at the surface to each level (mostly useful in combination with @@ -613,6 +615,7 @@ def refine_surface(self, xyz, level, padding_cells_by_level=None, pad_up=False, padding_cells_by_level = padding_cells_by_level[sorted] if self.dim == 2: + xyz = np.asarray(xyz) sorted = np.argsort(xyz[:, 0]) xyz = xyz[sorted] n_ps = len(xyz) @@ -623,7 +626,10 @@ def refine_surface(self, xyz, level, padding_cells_by_level=None, pad_up=False, 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) diff --git a/tests/tree/test_refine.py b/tests/tree/test_refine.py index ed600972d..96631a02b 100644 --- a/tests/tree/test_refine.py +++ b/tests/tree/test_refine.py @@ -400,3 +400,106 @@ def test_insert_errors(): # incorrect number of levels with pytest.raises(ValueError): mesh.insert_cells(x0s2d, [1, 1, 3], finalize=False) + + +def test_bounding_box(): + mesh1 = discretize.TreeMesh([32, 32]) + + xyz = np.random.rand(20, 2) + mesh1.refine_bounding_box(xyz, -1, 2) + + padding = (2/32) + x0 = xyz.min(axis=0) - padding + xF = xyz.max(axis=0) + padding + + mesh2 = discretize.TreeMesh([32, 32]) + mesh2.refine_box(x0, xF, -1) + + assert mesh1.equals(mesh2) + + +def test_bounding_box_errors(): + mesh1 = discretize.TreeMesh([32, 32]) + + xyz = np.random.rand(20, 3) + with pytest.raises(ValueError): + mesh1.refine_bounding_box(xyz, [-1, -2], [2, 3, 4]) + + with pytest.raises(IndexError): + mesh1.refine_bounding_box(xyz, 20) + + +def test_refine_points(): + mesh1 = discretize.TreeMesh([32, 32]) + point = [0.5, 0.5] + level = -1 + n_cell_pad = 3 + mesh1.refine_points(point, level, n_cell_pad) + + ball_rad = 1/32 * n_cell_pad + mesh2 = discretize.TreeMesh([32, 32]) + mesh2.refine_ball(point, ball_rad, level) + + assert mesh1.equals(mesh2) + + +def test_refine_points_errors(): + mesh1 = discretize.TreeMesh([32, 32]) + + point = [0.5, 0.5] + + with pytest.raises(ValueError): + mesh1.refine_points(point, [-1, -2], [2, 2, 2]) + + with pytest.raises(IndexError): + mesh1.refine_points(point, [20, -2], [2, 2]) + + +def test_refine_surface2D(): + + mesh1 = discretize.TreeMesh([32, 32]) + points = [[0.3, 0.3], [0.7, 0.3]] + mesh1.refine_surface(points, -1, [[0, 2]], pad_up=True, pad_down=True) + + mesh2 = discretize.TreeMesh([32, 32]) + x0 = [0.3, 0.3 - 2/32] + xF = [0.7, 0.3 + 2/32] + mesh2.refine_box(x0, xF, -1) + + assert mesh1.equals(mesh2) + + +def test_refine_surface3D(): + mesh1 = discretize.TreeMesh([32, 32, 32]) + points = [ + [0.3, 0.3, 0.5], + [0.7, 0.3, 0.5], + [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) + + mesh2 = discretize.TreeMesh([32, 32, 32]) + pad = np.array([1, 2, 3])/32 + x0 = [0.3, 0.3, 0.5] - pad + xF = [0.7, 0.7, 0.5] + pad + mesh2.refine_box(x0, xF, -1) + + assert mesh1.equals(mesh2) + + 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) + + assert mesh1.equals(mesh3) + + +def test_refine_surface_errors(): + mesh = discretize.TreeMesh([32, 32]) + points = [[0.3, 0.3], [0.7, 0.3]] + + with pytest.raises(ValueError): + mesh.refine_surface(points, [-1, -2], [0, 1, 2, 3]) + + with pytest.raises(IndexError): + mesh.refine_surface(points, 20) From 597a37fbc6c86dd5ceac27273c7a8befe44bb113 Mon Sep 17 00:00:00 2001 From: Joseph Capriotti Date: Wed, 23 Nov 2022 14:18:39 -0800 Subject: [PATCH 7/9] modify test to hit last line of code --- tests/tree/test_refine.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/tree/test_refine.py b/tests/tree/test_refine.py index 96631a02b..0fd42a92b 100644 --- a/tests/tree/test_refine.py +++ b/tests/tree/test_refine.py @@ -452,7 +452,7 @@ def test_refine_points_errors(): mesh1.refine_points(point, [-1, -2], [2, 2, 2]) with pytest.raises(IndexError): - mesh1.refine_points(point, [20, -2], [2, 2]) + mesh1.refine_points(point, [20, -2]) def test_refine_surface2D(): From 441419277d482b8b887c30b2ef5806961815acbd Mon Sep 17 00:00:00 2001 From: Joseph Capriotti Date: Wed, 30 Nov 2022 20:21:53 -0800 Subject: [PATCH 8/9] refine updates --- discretize/_extensions/tree.cpp | 177 ++++++++++++++++++++++ discretize/_extensions/tree.h | 8 + discretize/_extensions/tree.pxd | 1 + discretize/_extensions/tree_ext.pyx | 90 +++++++++++ discretize/tree_mesh.py | 185 +++++++++++------------ tests/tree/test_refine.py | 145 +++++++++++++++--- tutorials/mesh_generation/4_tree_mesh.py | 18 +-- 7 files changed, 500 insertions(+), 124 deletions(-) diff --git a/discretize/_extensions/tree.cpp b/discretize/_extensions/tree.cpp index a05de1d3a..c68c0fafb 100644 --- a/discretize/_extensions/tree.cpp +++ b/discretize/_extensions/tree.cpp @@ -685,6 +685,163 @@ void Cell::refine_triangle( } } +void Cell::refine_vert_triang_prism( + node_map_t& nodes, + double* x0, double* x1, double* x2, double h, + double* e0, double* e1, double* e2, double* t_norm, + int_t p_level, double *xs, double *ys, double* zs, bool diag_balance +){ + // Return If I'm at max_level or p_level + if (level >= p_level || level == max_level){ + return; + } + // check all the AABB faces + double v0[3], v1[3], v2[3], half[3]; + double vmin, vmax; + double p0, p1, p2, p3, pmin, pmax, rad; + for(int_t i=0; i < n_dim; ++i){ + v0[i] = x0[i] - location[i]; + v1[i] = x1[i] - location[i]; + vmin = std::min(v0[i], v1[i]); + vmax = std::max(v0[i], v1[i]); + v2[i] = x2[i] - location[i]; + vmin = std::min(vmin, v2[i]); + vmax = std::max(vmax, v2[i]); + if(i == 2){ + vmax += h; + } + half[i] = location[i] - points[0]->location[i]; + + // Bounding box check + if (vmin > half[i] || vmax < -half[i]){ + return; + } + } + // first do the 3 edge cross tests that apply in 2D and 3D + + // edge 0 cross z_hat + //p0 = e0[1] * v0[0] - e0[0] * v0[1]; + p1 = e0[1] * v1[0] - e0[0] * v1[1]; + p2 = e0[1] * v2[0] - e0[0] * v2[1]; + pmin = std::min(p1, p2); + pmax = std::max(p1, p2); + rad = std::abs(e0[1]) * half[0] + std::abs(e0[0]) * half[1]; + if (pmin > rad || pmax < -rad){ + return; + } + + // edge 1 cross z_hat + p0 = e1[1] * v0[0] - e1[0] * v0[1]; + p1 = e1[1] * v1[0] - e1[0] * v1[1]; + //p2 = e1[1] * v2[0] - e1[0] * v2[1]; + pmin = std::min(p0, p1); + pmax = std::max(p0, p1); + rad = std::abs(e1[1]) * half[0] + std::abs(e1[0]) * half[1]; + if (pmin > rad || pmax < -rad){ + return; + } + + // edge 2 cross z_hat + //p0 = e2[1] * v0[0] - e2[0] * v0[1]; + p1 = e2[1] * v1[0] - e2[0] * v1[1]; + p2 = e2[1] * v2[0] - e2[0] * v2[1]; + pmin = std::min(p1, p2); + pmax = std::max(p1, p2); + rad = std::abs(e2[1]) * half[0] + std::abs(e2[0]) * half[1]; + if (pmin > rad || pmax < -rad){ + return; + } + + // edge 0 cross x_hat + p0 = e0[2] * v0[1] - e0[1] * v0[2]; + p1 = e0[2] * v0[1] - e0[1] * (v0[2] + h); + p2 = e0[2] * v2[1] - e0[1] * v2[2]; + p3 = e0[2] * v2[1] - e0[1] * (v2[2] + h); + pmin = std::min(std::min(std::min(p0, p1), p2), p3); + pmax = std::max(std::max(std::max(p0, p1), p2), p3); + rad = std::abs(e0[2]) * half[1] + std::abs(e0[1]) * half[2]; + if (pmin > rad || pmax < -rad){ + return; + } + // edge 0 cross y_hat + p0 = -e0[2] * v0[0] + e0[0] * v0[2]; + p1 = -e0[2] * v0[0] + e0[0] * (v0[2] + h); + p2 = -e0[2] * v2[0] + e0[0] * v2[2]; + p3 = -e0[2] * v2[0] + e0[0] * (v2[2] + h); + pmin = std::min(std::min(std::min(p0, p1), p2), p3); + pmax = std::max(std::max(std::max(p0, p1), p2), p3); + rad = std::abs(e0[2]) * half[0] + std::abs(e0[0]) * half[2]; + if (pmin > rad || pmax < -rad){ + return; + } + // edge 1 cross x_hat + p0 = e1[2] * v0[1] - e1[1] * v0[2]; + p1 = e1[2] * v0[1] - e1[1] * (v0[2] + h); + p2 = e1[2] * v2[1] - e1[1] * v2[2]; + p3 = e1[2] * v2[1] - e1[1] * (v2[2] + h); + pmin = std::min(std::min(std::min(p0, p1), p2), p3); + pmax = std::max(std::max(std::max(p0, p1), p2), p3); + rad = std::abs(e1[2]) * half[1] + std::abs(e1[1]) * half[2]; + if (pmin > rad || pmax < -rad){ + return; + } + // edge 1 cross y_hat + p0 = -e1[2] * v0[0] + e1[0] * v0[2]; + p1 = -e1[2] * v0[0] + e1[0] * (v0[2] + h); + p2 = -e1[2] * v2[0] + e1[0] * v2[2]; + p3 = -e1[2] * v2[0] + e1[0] * (v2[2] + h); + pmin = std::min(std::min(std::min(p0, p1), p2), p3); + pmax = std::max(std::max(std::max(p0, p1), p2), p3); + rad = std::abs(e1[2]) * half[0] + std::abs(e1[0]) * half[2]; + if (pmin > rad || pmax < -rad){ + return; + } + // edge 2 cross x_hat + p0 = e2[2] * v0[1] - e2[1] * v0[2]; + p1 = e2[2] * v0[1] - e2[1] * (v0[2] + h); + p2 = e2[2] * v1[1] - e2[1] * v1[2]; + p3 = e2[2] * v1[1] - e2[1] * (v1[2] + h); + pmin = std::min(std::min(std::min(p0, p1), p2), p3); + pmax = std::max(std::max(std::max(p0, p1), p2), p3); + rad = std::abs(e2[2]) * half[1] + std::abs(e2[1]) * half[2]; + if (pmin > rad || pmax < -rad){ + return; + } + // edge 2 cross y_hat + p0 = -e2[2] * v0[0] + e2[0] * v0[2]; + p1 = -e2[2] * v0[0] + e2[0] * (v0[2] + h); + p2 = -e2[2] * v1[0] + e2[0] * v1[2]; + p3 = -e2[2] * v1[0] + e2[0] * (v1[2] + h); + pmin = std::min(std::min(std::min(p0, p1), p2), p3); + pmax = std::max(std::max(std::max(p0, p1), p2), p3); + rad = std::abs(e2[2]) * half[0] + std::abs(e2[0]) * half[2]; + if (pmin > rad || pmax < -rad){ + return; + } + + // triangle normal axis + p0 = t_norm[0] * v0[0] + t_norm[1] * v0[1] + t_norm[2] * v0[2]; + p1 = t_norm[0] * v0[0] + t_norm[1] * v0[1] + t_norm[2] * (v0[2] + h); + pmin = std::min(p0, p1); + pmax = std::max(p0, p1); + rad = std::abs(t_norm[0]) * half[0] + std::abs(t_norm[1]) * half[1] + std::abs(t_norm[2]) * half[2]; + if (pmin > rad || pmax < -rad){ + return; + } + // the axes defined by the three vertical prism faces + // should already be tested by the e0, e1, e2 cross z_hat tests + + // If here, then I intersect the triangle! + if(is_leaf()){ + divide(nodes, xs, ys, zs, true, diag_balance); + } + for(int_t i = 0; i < (1<refine_vert_triang_prism( + nodes, x0, x1, x2, h, e0, e1, e2, t_norm, p_level, xs, ys, zs, diag_balance + ); + } +} + void Cell::refine_tetra( node_map_t& nodes, double* x0, double* x1, double* x2, double* x3, @@ -1347,6 +1504,26 @@ void Tree::refine_triangle(double* x0, double* x1, double* x2, int_t p_level, bo ); }; +void Tree::refine_vert_triang_prism(double* x0, double* x1, double* x2, double h, int_t p_level, bool diagonal_balance){ + double e0[3], e1[3], e2[3], t_norm[3]; + for(int_t i=0; i 2){ + t_norm[0] = e0[1] * e1[2] - e0[2] * e1[1]; + t_norm[1] = e0[2] * e1[0] - e0[0] * e1[2]; + t_norm[2] = e0[0] * e1[1] - e0[1] * e1[0]; + } + for(int_t iz=0; izrefine_vert_triang_prism( + nodes, x0, x1, x2, h, e0, e1, e2, t_norm, p_level, xs, ys, zs, diagonal_balance + ); +}; + void Tree::refine_tetra(double* x0, double* x1, double* x2, double* x3, int_t p_level, bool diagonal_balance){ double t_edges[6][3]; double face_normals[4][3]; diff --git a/discretize/_extensions/tree.h b/discretize/_extensions/tree.h index 8f6fbe99b..f3b0bef93 100644 --- a/discretize/_extensions/tree.h +++ b/discretize/_extensions/tree.h @@ -132,6 +132,11 @@ class Cell{ double* x0, double* x1, double* x2, double* e0, double* e1, double* e2, double* t_norm, int_t p_level, double *xs, double *ys, double* zs, bool diag_balance=false ); + void refine_vert_triang_prism(node_map_t& nodes, + double* x0, double* x1, double* x2, double h, + double* e0, double* e1, double* e2, double* t_norm, + int_t p_level, double *xs, double *ys, double* zs, bool diag_balance=false + ); void refine_tetra( node_map_t& nodes, double* x0, double* x1, double* x2, double* x3, @@ -179,6 +184,9 @@ class Tree{ void refine_tetra( double* x0, double* x1, double* x2, double* x3, int_t p_level, bool diag_balance=false ); + void refine_vert_triang_prism( + double* x0, double* x1, double* x2, double h, int_t p_level, bool diagonal_balance=false + ); void number(); void finalize_lists(); diff --git a/discretize/_extensions/tree.pxd b/discretize/_extensions/tree.pxd index c7586a841..bb2fc35ce 100644 --- a/discretize/_extensions/tree.pxd +++ b/discretize/_extensions/tree.pxd @@ -89,6 +89,7 @@ cdef extern from "tree.h": void refine_box(double*, double*, int_t, bool) void refine_line(double*, double*, int_t, bool) void refine_triangle(double*, double*, double*, int_t, bool) + void refine_vert_triang_prism(double*, double*, double*, double, int_t, bool) void refine_tetra(double*, double*, double*, double*, int_t, bool) void number() void initialize_roots() diff --git a/discretize/_extensions/tree_ext.pyx b/discretize/_extensions/tree_ext.pyx index 114eecc31..da9696e9a 100644 --- a/discretize/_extensions/tree_ext.pyx +++ b/discretize/_extensions/tree_ext.pyx @@ -836,6 +836,96 @@ cdef class _TreeMesh: if finalize: self.finalize() + @cython.cdivision(True) + def refine_vertical_trianglular_prism(self, triangle, h, levels, finalize=True, diagonal_balance=None): + """Refines the :class:`~discretize.TreeMesh` along the trianglular prism to the desired level + + Refines the TreeMesh by determining if a cell intersects the given trianglular prism(s) + to the prescribed level(s). + + Parameters + ---------- + triangle : (N, 3, dim) array_like + The nodes of the bottom triangle(s). + h : (N) array_like + The height of the prism(s). + levels : int or (N) array_like of int + The level to refine intersecting cells to. + finalize : bool, optional + Whether to finalize after refining + diagonal_balance : bool or None, optional + Whether to balance cells diagonally in the refinement, `None` implies using + the same setting used to instantiate the TreeMesh`. + + Examples + -------- + We create a simple mesh and refine the TreeMesh such that all cells that + intersect the line segment path are at the given levels. + + >>> import discretize + >>> import matplotlib.pyplot as plt + >>> import matplotlib.patches as patches + >>> tree_mesh = discretize.TreeMesh([32, 32]) + >>> tree_mesh.max_level + 5 + + Next we define the points along the line and the level we want to refine to, + and refine the mesh. + + >>> triangle = [[0.14, 0.31], [0.32, 0.96], [0.23, 0.87]] + >>> levels = 5 + >>> tree_mesh.refine_triangle(triangle, levels) + + Now lets look at the mesh, and overlay the line on it to ensure it refined + where we wanted it to. + + >>> ax = tree_mesh.plot_grid() + >>> tri = patches.Polygon(triangle, fill=False) + >>> ax.add_patch(tri) + >>> plt.show() + + """ + if self.dim == 2: + raise NotImplementedError("refine_vertical_trianglular_prism only implemented in 3D.") + triangle = np.require(np.atleast_2d(triangle), dtype=np.float64, requirements="C") + if triangle.ndim == 2: + triangle = triangle[None, ...] + if triangle.shape[-1] != self.dim or triangle.shape[-2] != 3: + raise ValueError(f"triangle array must be (N, 3, {self.dim})") + cdef double[:, :, :] tris = triangle + + h = np.require(np.atleast_1d(h), dtype=np.float64, requirements="C") + levels = np.require(np.atleast_1d(levels), dtype=np.int32, + requirements='C') + cdef int n_triangles = triangle.shape[0]; + if levels.shape[0] == 1: + levels = np.full(n_triangles, levels[0], dtype=np.int32) + if h.shape[0] == 1: + h = np.full(n_triangles, h[0], dtype=np.float64) + if n_triangles != levels.shape[0]: + raise ValueError(f"inconsistent number of triangles {n_triangles} and levels {levels.shape[0]}") + if n_triangles != h.shape[0]: + raise ValueError(f"inconsistent number of triangles {n_triangles} and heights {h.shape[0]}") + if np.any(h < 0): + raise ValueError("All heights must be positive.") + + cdef int[:] ls = levels + cdef double[:] hs = h + + if diagonal_balance is None: + diagonal_balance = self._diagonal_balance + cdef bool diag_balance = diagonal_balance + + cdef int l + cdef int max_level = self.max_level + for i in range(n_triangles): + l = ls[i] + if l < 0: + l = (max_level + 1) - (abs(l) % (max_level + 1)) + self.tree.refine_vert_triang_prism(&tris[i, 0, 0], &tris[i, 1, 0], &tris[i, 2, 0], hs[i], l, diag_balance) + if finalize: + self.finalize() + @cython.cdivision(True) def refine_tetrahedron(self, tetra, levels, finalize=True, diagonal_balance=None): """Refines the :class:`~discretize.TreeMesh` along the tetrahedron to the desired level diff --git a/discretize/tree_mesh.py b/discretize/tree_mesh.py index 77bbf0ae0..9209bb8a4 100644 --- a/discretize/tree_mesh.py +++ b/discretize/tree_mesh.py @@ -360,7 +360,7 @@ def origin(self, value): # then update the TreeMesh with the hidden value self._set_origin(self._origin) - def refine_bounding_box(self, points, level, padding_cells_by_level=None, finalize=True, diagonal_balance=None): + def refine_bounding_box(self, points, level=-1, padding_cells_by_level=None, finalize=True, diagonal_balance=None): """Refine within a bounding box based on the maximum and minimum extent of scattered points. This function refines the tree mesh based on the bounding box defined by the @@ -374,15 +374,16 @@ def refine_bounding_box(self, points, level, padding_cells_by_level=None, finali ---------- points : (N, dim) array_like The bounding box will be the maximum and minimum extent of these points - level : (n_level) array_like of int - The level of the treemesh to refine the bounding box to. If an array it will - refine the box to each level (mostly useful in combination with `padding_cells_by_level`) - Negative values index backwards, (e.g. `-1` is `max_level`). - padding_cells_by_level : None or (n_level) iterable, optional + level : int, optional + The level of the treemesh to refine the bounding box to. + Negative values index tree levels backwards, (e.g. `-1` is `max_level`). + padding_cells_by_level : None, int, (n_level) array_like or (n_level, dim) array_like, optional The number of cells to pad the bounding box at each level of refinement. - Each entry may be either a single number, which means to pad equally in each - dimension, or it may be an (`dim`) `array_like` for variable padding along - each dimension. None implies no padding. + If a single number, each level below ``level`` will be padded with those + number of cells. If array_like, `n_level`s below ``level`` will be padded, + where if a 1D array, each dimension will be padded with the same number of + cells, or 2D array supports variable padding along each dimension. None + implies no padding. finalize : bool, optional Whether to finalize the mesh after the call. diagonal_balance : None or bool, optional @@ -405,9 +406,8 @@ def refine_bounding_box(self, points, level, padding_cells_by_level=None, finali direction and `2` cells in `y`, and the second highest level we want 2 padding cells in each direction beyond that. - >>> levels = [-1, -2] - >>> padding = [[0, 2], 2] - >>> mesh.refine_bounding_box(points, levels, padding) + >>> padding = [[0, 2], [2, 2]] + >>> mesh.refine_bounding_box(points, -1, padding) For demonstration we, overlay the bounding box to show the effects of padding. @@ -418,36 +418,37 @@ def refine_bounding_box(self, points, level, padding_cells_by_level=None, finali """ bsw = np.min(np.atleast_2d(points), axis=0) tnw = np.max(np.atleast_2d(points), axis=0) - level = np.atleast_1d(level) + if level < 0: + level = (self.max_level + 1) - (abs(level) % (self.max_level + 1)) + if level > self.max_level: + raise IndexError(f"Level beyond max octree level, {self.max_level}") # pad based on the number of cells at each level if padding_cells_by_level is None: - padding_cells_by_level = np.zeros_like(level) - padding_cells_by_level = np.asarray(padding_cells_by_level, dtype=object) + padding_cells_by_level = np.array([0]) + padding_cells_by_level = np.asarray(padding_cells_by_level) + if padding_cells_by_level.ndim == 0 and level > 1: + padding_cells_by_level = np.full(level - 1, padding_cells_by_level) padding_cells_by_level = np.atleast_1d(padding_cells_by_level) - if len(level) != len(padding_cells_by_level): - raise ValueError("level and padding_cells_by_level must be the same length") - if np.any(level > self.max_level): - raise IndexError(f"Level beyond max octree level, {self.max_level}") - sorted = np.argsort(level)[::-1] - level = level[sorted] - padding_cells_by_level = padding_cells_by_level[sorted] + if padding_cells_by_level.ndim == 2 and padding_cells_by_level.shape[1] != self.dim: + raise ValueError("incorrect dimension for padding_cells_by_level.") h_min = np.r_[[h.min() for h in self.h]] x0 = [] xF = [] - for lv, n_pad in zip(level, padding_cells_by_level): - if lv < 0: - lv = (self.max_level + 1) - (abs(lv) % (self.max_level + 1)) + ls = [] + # The zip will short circuit on the shorter array. + for lv, n_pad in zip(np.arange(level, 1, -1), padding_cells_by_level): padding_at_level = n_pad * h_min * 2 ** (self.max_level - lv) bsw = bsw - padding_at_level tnw = tnw + padding_at_level x0.append(bsw) xF.append(tnw) + ls.append(lv) - self.refine_box(x0, xF, level, finalize=finalize, diagonal_balance=diagonal_balance) + self.refine_box(x0, xF, ls, finalize=finalize, diagonal_balance=diagonal_balance) - def refine_points(self, points, level, padding_cells_by_level=None, finalize=True, diagonal_balance=None): + def refine_points(self, points, level=-1, padding_cells_by_level=None, finalize=True, diagonal_balance=None): """Refine the mesh at given points to the prescribed level. This function refines the tree mesh around the `points`. It will refine @@ -458,14 +459,13 @@ def refine_points(self, points, level, padding_cells_by_level=None, finalize=Tru ---------- points : (N, dim) array_like The points to be refined around. - level : int or (n_level) array_like of int - The level of the tree mesh to refine to. If an array, it will - refine at each point to each level (mostly useful in combination with - `padding_cells_by_level`). Negative values index backwards, (e.g. `-1` is - `max_level`). - padding_cells_by_level : None or (n_level) iterable, optional - The number of cells at each level to pad around the point. The base cell size - is the largest value amongst the smallest cell widths for each dimension. + level : int, optional + The level of the tree mesh to refine to. Negative values index tree levels + backwards, (e.g. `-1` is `max_level`). + padding_cells_by_level : None, int, (n_level) array_like, optional + The number of cells to pad the bounding box at each level of refinement. + If a single number, each level below ``level`` will be padded with those + number of cells. If array_like, `n_level`s below ``level`` will be padded. None implies no padding. finalize : bool, optional Whether to finalize the mesh after the call. @@ -487,42 +487,41 @@ def refine_points(self, points, level, padding_cells_by_level=None, finalize=Tru and then refine at the second highest level with at least 3 cell beyond that. >>> points = np.array([[0.1, 0.3], [0.6, 0.8]]) - >>> levels = [-1, -2] >>> padding = [1, 3] - >>> mesh.refine_points(points, levels, padding) + >>> mesh.refine_points(points, -1, padding) >>> ax = mesh.plot_grid() >>> ax.scatter(*points.T, color='C1') >>> plt.show() """ - level = np.atleast_1d(level) + if level < 0: + level = (self.max_level + 1) - (abs(level) % (self.max_level + 1)) + if level > self.max_level: + raise IndexError(f"Level beyond max octree level, {self.max_level}") # pad based on the number of cells at each level if padding_cells_by_level is None: - padding_cells_by_level = np.zeros_like(level) + padding_cells_by_level = np.array([0]) + padding_cells_by_level = np.asarray(padding_cells_by_level) + if padding_cells_by_level.ndim == 0 and level > 1: + padding_cells_by_level = np.full(level - 1, padding_cells_by_level) padding_cells_by_level = np.atleast_1d(padding_cells_by_level) - if len(level) != len(padding_cells_by_level): - raise ValueError("level and padding_cells_by_level must be the same length") - if np.any(level > self.max_level): - raise IndexError(f"Level beyond max octree level, {self.max_level}") - sorted = np.argsort(level)[::-1] - level = level[sorted] - padding_cells_by_level = padding_cells_by_level[sorted] h_min = np.max([h.min() for h in self.h]) radius_at_level = 0.0 - for lv, n_pad in zip(level, padding_cells_by_level): - if lv < 0: - lv = (self.max_level + 1) - (abs(lv) % (self.max_level + 1)) + for lv, n_pad in zip(np.arange(level, 1, -1), padding_cells_by_level): radius_at_level += n_pad * h_min * 2 ** (self.max_level - lv) - self.refine_ball( - points, radius_at_level, lv, finalize=False, diagonal_balance=diagonal_balance - ) + if radius_at_level == 0: + self.insert_cells(points, lv, finalize=False, diagonal_balance=diagonal_balance) + else: + self.refine_ball( + points, radius_at_level, lv, finalize=False, diagonal_balance=diagonal_balance + ) if finalize: self.finalize() - def refine_surface(self, xyz, level, padding_cells_by_level=None, pad_up=False, pad_down=True, finalize=True, diagonal_balance=None): + def refine_surface(self, xyz, level=-1, padding_cells_by_level=None, pad_up=False, pad_down=True, finalize=True, diagonal_balance=None): """Refine along a surface triangulated from xyz to the prescribed level. This function refines the mesh based on a triangulated surface from the input @@ -540,15 +539,15 @@ def refine_surface(self, xyz, level, padding_cells_by_level=None, pad_up=False, dimensions. You are able to supply your own triangulation in 3D by inputing a tuple of (`xyz`, `triangle_indices`). level : int or (n_level) array_like of int - The level of the tree mesh to refine to. If an array, it will - refine at the surface to each level (mostly useful in combination with - `padding_cells_by_level`). Negative values index backwards, (e.g. `-1` is - `max_level`). - padding_cells_by_level : None or (n_level) iterable, optional - The number of cells to pad the surface with at each level of refinement. - Each entry may be either a single number, which means to pad equally in each - dimension, or it may be an (`dim`) `array_like` for variable padding along - each dimension. None implies no padding. + The level of the tree mesh to refine to. Negative values index tree levels + backwards, (e.g. `-1` is `max_level`). + padding_cells_by_level : None, int, (n_level) array_like or (n_level, dim) array_like, optional + The number of cells to pad the bounding box at each level of refinement. + If a single number, each level below ``level`` will be padded with those + number of cells. If array_like, `n_level`s below ``level`` will be padded, + where if a 1D array, each dimension will be padded with the same number of + cells, or 2D array supports variable padding along each dimension. None + implies no padding. pad_up : bool, optional Whether to pad above the surface. pad_down : bool, optional @@ -576,7 +575,7 @@ def refine_surface(self, xyz, level, padding_cells_by_level=None, pad_up=False, >>> x = np.linspace(0.2, 0.8, 51) >>> z = 0.25*np.sin(2*np.pi*x)+0.5 >>> xz = np.c_[x, z] - >>> mesh.refine_surface(xz, [-1, -2], [[0, 2], [0, 3]]) + >>> mesh.refine_surface(xz, -1, [[0, 2], [0, 3]]) >>> ax = mesh.plot_grid() >>> ax.plot(x, z, color='C1') @@ -590,7 +589,7 @@ def refine_surface(self, xyz, level, padding_cells_by_level=None, pad_up=False, >>> x, y = np.mgrid[0.2:0.8:21j, 0.2:0.8:21j] >>> z = 0.125*np.sin(2*np.pi*x) + 0.5 + 0.125 * np.cos(2 * np.pi * y) >>> points = np.stack([x, y, z], axis=-1).reshape(-1, 3) - >>> mesh.refine_surface(points, [-1, -2], [[0, 0, 2], [0, 0, 3]]) + >>> mesh.refine_surface(points, -1, [[0, 0, 2], [0, 0, 3]]) >>> v = mesh.cell_levels_by_index(np.arange(mesh.n_cells)) >>> fig, axs = plt.subplots(1, 3, figsize=(12,4)) @@ -598,21 +597,21 @@ def refine_surface(self, xyz, level, padding_cells_by_level=None, pad_up=False, >>> mesh.plot_slice(v, ax=axs[1], normal='y', grid=True, clim=[2, 5]) >>> mesh.plot_slice(v, ax=axs[2], normal='z', grid=True, clim=[2, 5]) >>> plt.show() - """ - level = np.atleast_1d(level) + if level < 0: + level = (self.max_level + 1) - (abs(level) % (self.max_level + 1)) + if level > self.max_level: + raise IndexError(f"Level beyond max octree level, {self.max_level}") + # pad based on the number of cells at each level if padding_cells_by_level is None: - padding_cells_by_level = np.zeros_like(level) - padding_cells_by_level = np.asarray(padding_cells_by_level, dtype=object) + padding_cells_by_level = np.array([0]) + padding_cells_by_level = np.asarray(padding_cells_by_level) + if padding_cells_by_level.ndim == 0 and level > 1: + padding_cells_by_level = np.full(level - 1, padding_cells_by_level) padding_cells_by_level = np.atleast_1d(padding_cells_by_level) - if len(level) != len(padding_cells_by_level): - raise ValueError("level and padding_cells_by_level must be the same length") - if np.any(level > self.max_level): - raise IndexError(f"Level beyond max octree level, {self.max_level}") - sorted = np.argsort(level)[::-1] - level = level[sorted] - padding_cells_by_level = padding_cells_by_level[sorted] + if padding_cells_by_level.ndim == 2 and padding_cells_by_level.shape[1] != self.dim: + raise ValueError("incorrect dimension for padding_cells_by_level.") if self.dim == 2: xyz = np.asarray(xyz) @@ -633,39 +632,33 @@ def refine_surface(self, xyz, level, padding_cells_by_level=None, pad_up=False, triang = Delaunay(xyz[:, :2]) simps = triang.simplices n_ps = len(xyz) - simps1 = np.c_[ - simps[:, 0], simps[:, 1], simps[:, 2], simps[:, 0] - ] + [0, 0, 0, n_ps] - simps2 = np.c_[ - simps[:, 0], simps[:, 1], simps[:, 2], simps[:, 1] - ] + [n_ps, n_ps, n_ps, 0] - simps3 = np.c_[ - simps[:, 1], simps[:, 2], simps[:, 0], simps[:, 2] - ] + [0, 0, n_ps, n_ps] - simps = np.r_[simps1, simps2, simps3] # 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((2, n_ps, self.dim)) - points[:, :, -1] = xyz[:, -1] + points = np.empty((n_ps, self.dim)) + points[:, -1] = xyz[:, -1] h_min = np.r_[[h.min() for h in self.h]] pad = 0.0 - for lv, n_pad in zip(level, padding_cells_by_level): - if lv < 0: - lv = (self.max_level + 1) - (abs(lv) % (self.max_level + 1)) + 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 pad_up: - points[0, :, -1] = xyz[:, -1] + pad[-1] + h += pad[-1] if pad_down: - points[1, :, -1] = xyz[:, -1] - pad[-1] - expansion = (half_width + pad[:-1]) / half_width - points[:, :, :-1] = expansion * (xyz[:, :-1] - center) + center - tetrahedrons = points.reshape((-1, self.dim))[simps] - self.refine_tetrahedron(tetrahedrons, lv, finalize=False, diagonal_balance=diagonal_balance) + 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(triangles, lv, finalize=False, diagonal_balance=diagonal_balance) + else: + triangles = points[simps] + self.refine_vertical_trianglular_prism(triangles, h, lv, finalize=False, diagonal_balance=diagonal_balance) if finalize: self.finalize() diff --git a/tests/tree/test_refine.py b/tests/tree/test_refine.py index 0fd42a92b..8e4b5068f 100644 --- a/tests/tree/test_refine.py +++ b/tests/tree/test_refine.py @@ -402,29 +402,115 @@ def test_insert_errors(): mesh.insert_cells(x0s2d, [1, 1, 3], finalize=False) +def test_refine_triang_prism(): + xyz = np.array([ + [0.41, 0.21, 0.11], + [0.21, 0.61, 0.22], + [0.71, 0.71, 0.31], + ]) + h = 0.48 + + simps = np.array([[0, 1, 2]]) + + n_ps = len(xyz) + simps1 = np.c_[ + simps[:, 0], simps[:, 1], simps[:, 2], simps[:, 0] + ] + [0, 0, 0, n_ps] + simps2 = np.c_[ + simps[:, 0], simps[:, 1], simps[:, 2], simps[:, 1] + ] + [n_ps, n_ps, n_ps, 0] + simps3 = np.c_[ + simps[:, 1], simps[:, 2], simps[:, 0], simps[:, 2] + ] + [0, 0, n_ps, n_ps] + simps = np.r_[simps1, simps2, simps3] + + points = np.r_[xyz, xyz + [0, 0, h]] + mesh1 = discretize.TreeMesh([32, 32, 32]) + mesh1.refine_tetrahedron(points[simps], -1) + + mesh2 = discretize.TreeMesh([32, 32, 32]) + mesh2.refine_vertical_trianglular_prism(xyz, h, -1) + + assert mesh1.equals(mesh2) + + +def test_refine_triang_prism_errors(): + xyz = np.array([ + [0.41, 0.21, 0.11], + [0.21, 0.61, 0.22], + [0.71, 0.71, 0.31], + ]) + h = 0.48 + + mesh = discretize.TreeMesh([32, 32]) + # Not implemented for 2D + with pytest.raises(NotImplementedError): + mesh.refine_vertical_trianglular_prism(xyz, h, -1) + + mesh = discretize.TreeMesh([32, 32, 32]) + # incorrect triangle dimensions + with pytest.raises(ValueError): + mesh.refine_vertical_trianglular_prism(xyz[:, :-1], h, -1) + + # incorrect levels and triangles + ps = np.random.rand(10, 3, 3) + with pytest.raises(ValueError): + mesh.refine_vertical_trianglular_prism(ps, h, [-1, -2]) + + # incorrect heights and triangles + ps = np.random.rand(10, 3, 3) + with pytest.raises(ValueError): + mesh.refine_vertical_trianglular_prism(ps, [h, h], -1) + + # negative heights + ps = np.random.rand(10, 3, 3) + with pytest.raises(ValueError): + mesh.refine_vertical_trianglular_prism(ps, -h, -1) + + def test_bounding_box(): - mesh1 = discretize.TreeMesh([32, 32]) - xyz = np.random.rand(20, 2) - mesh1.refine_bounding_box(xyz, -1, 2) + # No padding + xyz = np.random.rand(20, 2) * 0.25 + 3/8 + mesh1 = discretize.TreeMesh([32, 32]) + mesh1.refine_bounding_box(xyz, -1, None) - padding = (2/32) - x0 = xyz.min(axis=0) - padding - xF = xyz.max(axis=0) + padding + x0 = xyz.min(axis=0) + xF = xyz.max(axis=0) mesh2 = discretize.TreeMesh([32, 32]) mesh2.refine_box(x0, xF, -1) assert mesh1.equals(mesh2) + # Padding at all levels + n_cell_pad = 2 + x0 = xyz.min(axis=0) + xF = xyz.max(axis=0) + + mesh1 = discretize.TreeMesh([32, 32]) + mesh1.refine_bounding_box(xyz, -1, n_cell_pad) + + mesh2 = discretize.TreeMesh([32, 32]) + for lv in range(mesh2.max_level, 1, -1): + padding = n_cell_pad * (2**(mesh2.max_level - lv)/32) + x0 -= padding + xF += padding + mesh2.refine_box(x0, xF, lv, finalize=False) + mesh2.finalize() + + assert mesh1.equals(mesh2) + def test_bounding_box_errors(): mesh1 = discretize.TreeMesh([32, 32]) xyz = np.random.rand(20, 3) + # incorrect padding shape with pytest.raises(ValueError): - mesh1.refine_bounding_box(xyz, [-1, -2], [2, 3, 4]) + mesh1.refine_bounding_box(xyz, -1, [[2, 3, 4]]) + # bad level with pytest.raises(IndexError): mesh1.refine_bounding_box(xyz, 20) @@ -434,11 +520,22 @@ def test_refine_points(): point = [0.5, 0.5] level = -1 n_cell_pad = 3 + mesh1.refine_points(point, level, None) + + mesh2 = discretize.TreeMesh([32, 32]) + mesh2.insert_cells(point, -1) + + assert mesh1.equals(mesh2) + + mesh1 = discretize.TreeMesh([32, 32]) mesh1.refine_points(point, level, n_cell_pad) - ball_rad = 1/32 * n_cell_pad mesh2 = discretize.TreeMesh([32, 32]) - mesh2.refine_ball(point, ball_rad, level) + ball_rad = 0.0 + for lv in range(mesh2.max_level, 1, -1): + ball_rad += 2**(mesh2.max_level - lv)/32 * n_cell_pad + mesh2.refine_ball(point, ball_rad, lv, finalize=False) + mesh2.finalize() assert mesh1.equals(mesh2) @@ -448,26 +545,40 @@ def test_refine_points_errors(): point = [0.5, 0.5] - with pytest.raises(ValueError): - mesh1.refine_points(point, [-1, -2], [2, 2, 2]) - with pytest.raises(IndexError): - mesh1.refine_points(point, [20, -2]) + mesh1.refine_points(point, 20) def test_refine_surface2D(): mesh1 = discretize.TreeMesh([32, 32]) points = [[0.3, 0.3], [0.7, 0.3]] - mesh1.refine_surface(points, -1, [[0, 2]], pad_up=True, pad_down=True) + mesh1.refine_surface(points, -1, None, pad_up=True, pad_down=True) mesh2 = discretize.TreeMesh([32, 32]) - x0 = [0.3, 0.3 - 2/32] - xF = [0.7, 0.3 + 2/32] + x0 = [0.3, 0.3] + xF = [0.7, 0.3] mesh2.refine_box(x0, xF, -1) assert mesh1.equals(mesh2) + 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) + + mesh2 = discretize.TreeMesh([32, 32]) + x0 = np.r_[0.3, 0.3] + xF = np.r_[0.7, 0.3] + for lv in range(mesh2.max_level, 1, -1): + pad = 2**(mesh2.max_level - lv)/32 * n_cell_pad + x0 -= pad + xF += pad + mesh2.refine_box(x0, xF, lv, finalize=False) + mesh2.finalize() + + assert mesh1.equals(mesh2) + def test_refine_surface3D(): mesh1 = discretize.TreeMesh([32, 32, 32]) @@ -499,7 +610,7 @@ def test_refine_surface_errors(): points = [[0.3, 0.3], [0.7, 0.3]] with pytest.raises(ValueError): - mesh.refine_surface(points, [-1, -2], [0, 1, 2, 3]) + mesh.refine_surface(points, -1, [[0, 1, 2, 3]]) with pytest.raises(IndexError): mesh.refine_surface(points, 20) diff --git a/tutorials/mesh_generation/4_tree_mesh.py b/tutorials/mesh_generation/4_tree_mesh.py index bf08195ad..5792861eb 100644 --- a/tutorials/mesh_generation/4_tree_mesh.py +++ b/tutorials/mesh_generation/4_tree_mesh.py @@ -62,9 +62,8 @@ # Discretize to finest cell size within rectangular box, with padding in the z direction # at the finest and second finest levels. -levels = [-1, -2] padding = [[0, 2], [0, 2]] -mesh.refine_bounding_box(xy, levels, padding) +mesh.refine_bounding_box(xy, level=-1, padding_cells_by_level=padding) mesh.plotGrid(show_it=True) @@ -102,15 +101,14 @@ xx = mesh.vectorNx yy = -3 * np.exp((xx ** 2) / 100 ** 2) + 50.0 pts = np.c_[mkvc(xx), mkvc(yy)] -levels = [-1, -2] padding = [[0, 2], [0, 2]] -mesh.refine_surface(pts, levels, padding, finalize=False) +mesh.refine_surface(pts, padding_cells_by_level=padding, finalize=False) # Refine mesh near points xx = np.array([0.0, 10.0, 0.0, -10.0]) yy = np.array([-20.0, -10.0, 0.0, -10]) pts = np.c_[mkvc(xx), mkvc(yy)] -mesh.refine_points(pts, [-1, -2], [2, 2], finalize=False) +mesh.refine_points(pts, padding_cells_by_level=[2, 2], finalize=False) mesh.finalize() @@ -149,15 +147,14 @@ xx = mesh.vectorNx yy = -3 * np.exp((xx ** 2) / 100 ** 2) + 50.0 pts = np.c_[mkvc(xx), mkvc(yy)] -levels = [-1 , -2] padding = [[0, 2], [0, 2]] -mesh.refine_surface(pts, levels, padding, finalize=False) +mesh.refine_surface(pts, padding_cells_by_level=padding, finalize=False) # Refine near points xx = np.array([0.0, 10.0, 0.0, -10.0]) yy = np.array([-20.0, -10.0, 0.0, -10]) pts = np.c_[mkvc(xx), mkvc(yy)] -mesh.refine_points(pts, [-1, -2], [2, 2], finalize=False) +mesh.refine_points(pts, padding_cells_by_level=[2, 2], finalize=False) mesh.finalize() @@ -214,14 +211,13 @@ [xx, yy] = np.meshgrid(mesh.vectorNx, mesh.vectorNy) zz = -3 * np.exp((xx ** 2 + yy ** 2) / 100 ** 2) + 50.0 pts = np.c_[mkvc(xx), mkvc(yy), mkvc(zz)] -levels = [-1, -2] padding = [[0, 0, 2], [0, 0, 2]] -mesh.refine_surface(pts, levels, padding, finalize=False) +mesh.refine_surface(pts, padding_cells_by_level=padding, finalize=False) # Refine box xp, yp, zp = np.meshgrid([-40.0, 40.0], [-40.0, 40.0], [-60.0, 0.0]) xyz = np.c_[mkvc(xp), mkvc(yp), mkvc(zp)] -mesh.refine_bounding_box(xyz, levels, padding, finalize=False) +mesh.refine_bounding_box(xyz, padding_cells_by_level=padding, finalize=False) mesh.finalize() # The bottom west corner From b3d059cd6158b2da6316d1fb137ff5c9f92bfbff Mon Sep 17 00:00:00 2001 From: Joseph Capriotti Date: Thu, 1 Dec 2022 11:34:31 -0800 Subject: [PATCH 9/9] documentation updates --- discretize/_extensions/tree_ext.pyx | 26 ++++++++------ discretize/tree_mesh.py | 54 ++++++++++++++++++++++++----- 2 files changed, 61 insertions(+), 19 deletions(-) diff --git a/discretize/_extensions/tree_ext.pyx b/discretize/_extensions/tree_ext.pyx index da9696e9a..e6f94afdb 100644 --- a/discretize/_extensions/tree_ext.pyx +++ b/discretize/_extensions/tree_ext.pyx @@ -857,6 +857,10 @@ cdef class _TreeMesh: Whether to balance cells diagonally in the refinement, `None` implies using the same setting used to instantiate the TreeMesh`. + See Also + -------- + refine_surface + Examples -------- We create a simple mesh and refine the TreeMesh such that all cells that @@ -865,23 +869,25 @@ cdef class _TreeMesh: >>> import discretize >>> import matplotlib.pyplot as plt >>> import matplotlib.patches as patches - >>> tree_mesh = discretize.TreeMesh([32, 32]) + >>> tree_mesh = discretize.TreeMesh([32, 32, 32]) >>> tree_mesh.max_level 5 - Next we define the points along the line and the level we want to refine to, - and refine the mesh. + Next we define the bottom points of the prism, its heights, and the level we + want to refine to, then refine the mesh. - >>> triangle = [[0.14, 0.31], [0.32, 0.96], [0.23, 0.87]] + >>> triangle = [[0.14, 0.31, 0.21], [0.32, 0.96, 0.34], [0.87, 0.23, 0.12]] + >>> height = 0.35 >>> levels = 5 - >>> tree_mesh.refine_triangle(triangle, levels) + >>> mesh.refine_vertical_trianglular_prism(triangle, height, levels) - Now lets look at the mesh, and overlay the line on it to ensure it refined - where we wanted it to. + Now lets look at the mesh. - >>> ax = tree_mesh.plot_grid() - >>> tri = patches.Polygon(triangle, fill=False) - >>> ax.add_patch(tri) + >>> v = mesh.cell_levels_by_index(np.arange(mesh.n_cells)) + >>> fig, axs = plt.subplots(1, 3, figsize=(12,4)) + >>> mesh.plot_slice(v, ax=axs[0], normal='x', grid=True, clim=[2, 5]) + >>> mesh.plot_slice(v, ax=axs[1], normal='y', grid=True, clim=[2, 5]) + >>> mesh.plot_slice(v, ax=axs[2], normal='z', grid=True, clim=[2, 5]) >>> plt.show() """ diff --git a/discretize/tree_mesh.py b/discretize/tree_mesh.py index 9209bb8a4..ce8668613 100644 --- a/discretize/tree_mesh.py +++ b/discretize/tree_mesh.py @@ -111,8 +111,32 @@ class TreeMesh( larger than some base cell dimension. Unlike the :class:`~discretize.TensorMesh` class, gridded locations and numerical operators for instances of ``TreeMesh`` cannot be simply constructed using tensor products. Furthermore, each cell - is an instance of ``TreeMesh`` is an instance of the - :class:`~discretize.tree_mesh.TreeCell` . + is an instance of :class:`~discretize.tree_mesh.TreeCell` . + + Each cell of a `TreeMesh` has a certain level associated with it which describes its + height in the structure of the tree. The tree mesh is refined by specifying what + level you want in a given location. They start at level 0, defining the largest + possible cell on the mesh, and go to a level of `max_level` describing the + finest possible cell on the mesh. TreeMesh` contains several refinement functions + used to design the mesh, starting with a general purpose refinement function, along + with several faster specialized refinement functions based on fundamental geometric + entities: + + - `refine` -> The general purpose refinement function supporting arbitrary defined functions + - `insert_cells` + - `refine_ball` + - `refine_line` + - `refine_triangle` + - `refine_box` + - `refine_tetrahedron` + - `refine_vertical_trianglular_prism` + - `refine_bounding_box` + - `refine_points` + - `refine_surface` + + Like array indexing in python, you can also supply negative indices as a level + arguments to these functions to index levels in a reveresed order (i.e. -1 is + equivalent to `max_level`). Parameters ---------- @@ -390,6 +414,10 @@ def refine_bounding_box(self, points, level=-1, padding_cells_by_level=None, fin Whether to balance cells diagonally in the refinement, `None` implies using the same setting used to instantiate the TreeMesh`. + See Also + -------- + refine_box + Examples -------- Given a set of points, we want to refine the tree mesh with the bounding box @@ -402,8 +430,8 @@ def refine_bounding_box(self, points, level=-1, padding_cells_by_level=None, fin >>> mesh = discretize.TreeMesh([32, 32]) >>> points = np.random.rand(20, 2) * 0.25 + 3/8 - Now we want to refine to the maximum level, with a no padding the in `x` - direction and `2` cells in `y`, and the second highest level we want 2 padding + Now we want to refine to the maximum level, with no padding the in `x` + direction and `2` cells in `y`. At the second highest level we want 2 padding cells in each direction beyond that. >>> padding = [[0, 2], [2, 2]] @@ -473,6 +501,10 @@ def refine_points(self, points, level=-1, padding_cells_by_level=None, finalize= Whether to balance cells diagonally in the refinement, `None` implies using the same setting used to instantiate the TreeMesh`. + See Also + -------- + refine_ball, insert_cells + Examples -------- Given a set of points, we want to refine the tree mesh around these points to @@ -484,7 +516,7 @@ def refine_points(self, points, level=-1, padding_cells_by_level=None, finalize= >>> mesh = discretize.TreeMesh([32, 32]) Now we want to refine to the maximum level with 1 cell padding around the point, - and then refine at the second highest level with at least 3 cell beyond that. + and then refine at the second highest level with at least 3 cells beyond that. >>> points = np.array([[0.1, 0.3], [0.6, 0.8]]) >>> padding = [1, 3] @@ -536,9 +568,9 @@ def refine_surface(self, xyz, level=-1, padding_cells_by_level=None, pad_up=Fals ---------- xyz : (N, dim) array_like or tuple The points defining the surface. Will be triangulated along the horizontal - dimensions. You are able to supply your own triangulation in 3D by inputing + dimensions. You are able to supply your own triangulation in 3D by passing a tuple of (`xyz`, `triangle_indices`). - level : int or (n_level) array_like of int + level : int, optional The level of the tree mesh to refine to. Negative values index tree levels backwards, (e.g. `-1` is `max_level`). padding_cells_by_level : None, int, (n_level) array_like or (n_level, dim) array_like, optional @@ -558,6 +590,10 @@ def refine_surface(self, xyz, level=-1, padding_cells_by_level=None, pad_up=Fals Whether to balance cells diagonally in the refinement, `None` implies using the same setting used to instantiate the TreeMesh`. + See Also + -------- + refine_triangle, refine_vertical_trianglular_prism + Examples -------- In 2D we define the surface as a line segment, which we would like to refine @@ -568,8 +604,8 @@ def refine_surface(self, xyz, level=-1, padding_cells_by_level=None, pad_up=Fals >>> import matplotlib.patches as patches >>> mesh = discretize.TreeMesh([32, 32]) - This surface is a simple sine curve. Which we would like to pad with at least - 2 cells vertically below at the maximum level, with 3 cells below that at the + This surface is a simple sine curve, which we would like to pad with at least + 2 cells vertically below at the maximum level, and 3 cells below that at the second highest level. >>> x = np.linspace(0.2, 0.8, 51)