From b71f252330a680d9b8864ce59ad8d3a6179715b5 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Fri, 30 Aug 2024 11:10:32 -0700 Subject: [PATCH 1/6] Add `TensorMesh.cell_bounds` property Add a new `cell_bounds` property to `TensorMesh` that returns the bounds of each cell in the mesh. This implements a faster way to obtain the bounds of each cell in the mesh than trying to access them through `TensorCell.bounds`. --- discretize/tensor_mesh.py | 32 ++++++++++++++++++++++++++++++++ tests/base/test_tensor.py | 14 +++++++++++--- 2 files changed, 43 insertions(+), 3 deletions(-) diff --git a/discretize/tensor_mesh.py b/discretize/tensor_mesh.py index 842dd0106..5862cf93d 100644 --- a/discretize/tensor_mesh.py +++ b/discretize/tensor_mesh.py @@ -590,6 +590,38 @@ def face_boundary_indices(self): indzu = self.gridFz[:, 2] == max(self.gridFz[:, 2]) return indxd, indxu, indyd, indyu, indzd, indzu + @property + def cell_bounds(self): + """The bounds of each cell. + + Return a 2D array with the coordinates that define the bounds of each + cell in the mesh. Each row of the array contains the bounds for + a particular cell in the following order: ``x1``, ``x2``, ``y1``, + ``y2``, ``z1``, ``z2``, where ``x1 < x2``, ``y1 < y2`` and ``z1 < z2``. + """ + # Define indexing for meshgrid and order for ravel + indexing = "ij" + order = "F" + + if self.dim == 1: + x = self.nodes_x + bounds = (x[:-1], x[1:]) + elif self.dim == 2: + x, y = self.nodes_x, self.nodes_y + x1, y1 = np.meshgrid(x[:-1], y[:-1], indexing=indexing) + x2, y2 = np.meshgrid(x[1:], y[1:], indexing=indexing) + bounds = (x1, x2, y1, y2) + else: + x, y, z = self.nodes_x, self.nodes_y, self.nodes_z + x1, y1, z1 = np.meshgrid(x[:-1], y[:-1], z[:-1], indexing=indexing) + x2, y2, z2 = np.meshgrid(x[1:], y[1:], z[1:], indexing=indexing) + bounds = (x1, x2, y1, y2, z1, z2) + + cell_bounds = np.hstack( + tuple(v.ravel(order=order)[:, np.newaxis] for v in bounds) + ) + return cell_bounds + @property def cell_nodes(self): """The index of all nodes for each cell. diff --git a/tests/base/test_tensor.py b/tests/base/test_tensor.py index 9349cafa7..3d5e0039d 100644 --- a/tests/base/test_tensor.py +++ b/tests/base/test_tensor.py @@ -251,9 +251,9 @@ def test_serialization(self): self.assertTrue(np.all(self.mesh2.gridCC == mesh.gridCC)) -class TestTensorMeshCellNodes: +class TestTensorMeshProperties: """ - Test TensorMesh.cell_nodes + Test some of the properties in TensorMesh """ @pytest.fixture(params=[1, 2, 3], ids=["dims-1", "dims-2", "dims-3"]) @@ -261,17 +261,25 @@ def mesh(self, request): """Sample TensorMesh.""" if request.param == 1: h = [10] + origin = (-35.5,) elif request.param == 2: h = [10, 15] + origin = (-35.5, 105.3) else: h = [10, 15, 20] - return discretize.TensorMesh(h) + origin = (-35.5, 105.3, -27.3) + return discretize.TensorMesh(h, origin=origin) def test_cell_nodes(self, mesh): """Test TensorMesh.cell_nodes.""" expected_cell_nodes = np.array([cell.nodes for cell in mesh]) np.testing.assert_allclose(mesh.cell_nodes, expected_cell_nodes) + def test_cell_bounds(self, mesh): + """Test TensorMesh.cell_bounds.""" + expected_cell_bounds = np.array([cell.bounds for cell in mesh]) + np.testing.assert_allclose(mesh.cell_bounds, expected_cell_bounds) + class TestPoissonEqn(discretize.tests.OrderTest): name = "Poisson Equation" From 4c45e59e194699af7911da6611f71f4e5fae2e05 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Thu, 5 Sep 2024 10:22:53 -0700 Subject: [PATCH 2/6] Simplify the cell_nodes function Instead of building the nodes from scratch with meshgrid, use the existing `cell_centers` and `h_gridded` properties of the `TensorMesh`. --- discretize/tensor_mesh.py | 28 +++++++++++++--------------- 1 file changed, 13 insertions(+), 15 deletions(-) diff --git a/discretize/tensor_mesh.py b/discretize/tensor_mesh.py index 00899e7c8..8aef0e93f 100644 --- a/discretize/tensor_mesh.py +++ b/discretize/tensor_mesh.py @@ -602,27 +602,25 @@ def cell_bounds(self): a particular cell in the following order: ``x1``, ``x2``, ``y1``, ``y2``, ``z1``, ``z2``, where ``x1 < x2``, ``y1 < y2`` and ``z1 < z2``. """ - # Define indexing for meshgrid and order for ravel - indexing = "ij" - order = "F" - + centers, widths = self.cell_centers, self.h_gridded if self.dim == 1: - x = self.nodes_x - bounds = (x[:-1], x[1:]) + bounds = (centers - widths.ravel() / 2, centers + widths.ravel() /2) elif self.dim == 2: - x, y = self.nodes_x, self.nodes_y - x1, y1 = np.meshgrid(x[:-1], y[:-1], indexing=indexing) - x2, y2 = np.meshgrid(x[1:], y[1:], indexing=indexing) + x1 = centers[:, 0] - widths[:, 0] / 2 + x2 = centers[:, 0] + widths[:, 0] / 2 + y1 = centers[:, 1] - widths[:, 1] / 2 + y2 = centers[:, 1] + widths[:, 1] / 2 bounds = (x1, x2, y1, y2) else: - x, y, z = self.nodes_x, self.nodes_y, self.nodes_z - x1, y1, z1 = np.meshgrid(x[:-1], y[:-1], z[:-1], indexing=indexing) - x2, y2, z2 = np.meshgrid(x[1:], y[1:], z[1:], indexing=indexing) + x1 = centers[:, 0] - widths[:, 0] / 2 + x2 = centers[:, 0] + widths[:, 0] / 2 + y1 = centers[:, 1] - widths[:, 1] / 2 + y2 = centers[:, 1] + widths[:, 1] / 2 + z1 = centers[:, 2] - widths[:, 2] / 2 + z2 = centers[:, 2] + widths[:, 2] / 2 bounds = (x1, x2, y1, y2, z1, z2) - cell_bounds = np.hstack( - tuple(v.ravel(order=order)[:, np.newaxis] for v in bounds) - ) + cell_bounds = np.hstack(tuple(v[:, np.newaxis] for v in bounds)) return cell_bounds @property From e1c92b006c74ed2a7e0385c20139542326c5fb75 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Thu, 5 Sep 2024 10:33:56 -0700 Subject: [PATCH 3/6] Run black --- discretize/tensor_mesh.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/discretize/tensor_mesh.py b/discretize/tensor_mesh.py index 8aef0e93f..6adb82e52 100644 --- a/discretize/tensor_mesh.py +++ b/discretize/tensor_mesh.py @@ -604,7 +604,7 @@ def cell_bounds(self): """ centers, widths = self.cell_centers, self.h_gridded if self.dim == 1: - bounds = (centers - widths.ravel() / 2, centers + widths.ravel() /2) + bounds = (centers - widths.ravel() / 2, centers + widths.ravel() / 2) elif self.dim == 2: x1 = centers[:, 0] - widths[:, 0] / 2 x2 = centers[:, 0] + widths[:, 0] / 2 From 49f794b3c64fe802f44448533fc1d2f3d15e7509 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 15 Oct 2024 15:07:33 -0700 Subject: [PATCH 4/6] Cleaner implementation of cell_bounds Co-authored-by: Joseph Capriotti --- discretize/tensor_mesh.py | 28 +++++++++------------------- 1 file changed, 9 insertions(+), 19 deletions(-) diff --git a/discretize/tensor_mesh.py b/discretize/tensor_mesh.py index 6adb82e52..54343f13f 100644 --- a/discretize/tensor_mesh.py +++ b/discretize/tensor_mesh.py @@ -602,25 +602,15 @@ def cell_bounds(self): a particular cell in the following order: ``x1``, ``x2``, ``y1``, ``y2``, ``z1``, ``z2``, where ``x1 < x2``, ``y1 < y2`` and ``z1 < z2``. """ - centers, widths = self.cell_centers, self.h_gridded - if self.dim == 1: - bounds = (centers - widths.ravel() / 2, centers + widths.ravel() / 2) - elif self.dim == 2: - x1 = centers[:, 0] - widths[:, 0] / 2 - x2 = centers[:, 0] + widths[:, 0] / 2 - y1 = centers[:, 1] - widths[:, 1] / 2 - y2 = centers[:, 1] + widths[:, 1] / 2 - bounds = (x1, x2, y1, y2) - else: - x1 = centers[:, 0] - widths[:, 0] / 2 - x2 = centers[:, 0] + widths[:, 0] / 2 - y1 = centers[:, 1] - widths[:, 1] / 2 - y2 = centers[:, 1] + widths[:, 1] / 2 - z1 = centers[:, 2] - widths[:, 2] / 2 - z2 = centers[:, 2] + widths[:, 2] / 2 - bounds = (x1, x2, y1, y2, z1, z2) - - cell_bounds = np.hstack(tuple(v[:, np.newaxis] for v in bounds)) + nodes = self.nodes.reshape((*self.shape_nodes, -1), order="F") + + min_nodes = nodes[(slice(-1), )*self.dim] + min_nodes = min_nodes.reshape((self.n_cells, -1), order="F") + max_nodes = nodes[(slice(1, None), )*self.dim] + max_nodes = max_nodes.reshape((self.n_cells, -1), order="F") + + cell_bounds = np.stack((min_nodes, max_nodes), axis=-1) + cell_bounds = cell_bounds.reshape((self.n_cells, -1)) return cell_bounds @property From dfe5a54d7e0deac8aabbe7e06d9b5428594f61f9 Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 15 Oct 2024 15:10:05 -0700 Subject: [PATCH 5/6] Replace assert_allclose for assert_equal --- tests/base/test_tensor.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/tests/base/test_tensor.py b/tests/base/test_tensor.py index b9ff8ac2c..b0c30dec4 100644 --- a/tests/base/test_tensor.py +++ b/tests/base/test_tensor.py @@ -275,12 +275,12 @@ def mesh(self, request): def test_cell_nodes(self, mesh): """Test TensorMesh.cell_nodes.""" expected_cell_nodes = np.array([cell.nodes for cell in mesh]) - np.testing.assert_allclose(mesh.cell_nodes, expected_cell_nodes) + np.testing.assert_equal(mesh.cell_nodes, expected_cell_nodes) def test_cell_bounds(self, mesh): """Test TensorMesh.cell_bounds.""" expected_cell_bounds = np.array([cell.bounds for cell in mesh]) - np.testing.assert_allclose(mesh.cell_bounds, expected_cell_bounds) + np.testing.assert_equal(mesh.cell_bounds, expected_cell_bounds) class TestPoissonEqn(discretize.tests.OrderTest): From 7bd7d1b2701904c51e57ef3e0dd0f29364b806fa Mon Sep 17 00:00:00 2001 From: Santiago Soler Date: Tue, 15 Oct 2024 15:14:41 -0700 Subject: [PATCH 6/6] Run black --- discretize/tensor_mesh.py | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/discretize/tensor_mesh.py b/discretize/tensor_mesh.py index 54343f13f..c8ddb3278 100644 --- a/discretize/tensor_mesh.py +++ b/discretize/tensor_mesh.py @@ -603,12 +603,12 @@ def cell_bounds(self): ``y2``, ``z1``, ``z2``, where ``x1 < x2``, ``y1 < y2`` and ``z1 < z2``. """ nodes = self.nodes.reshape((*self.shape_nodes, -1), order="F") - - min_nodes = nodes[(slice(-1), )*self.dim] + + min_nodes = nodes[(slice(-1),) * self.dim] min_nodes = min_nodes.reshape((self.n_cells, -1), order="F") - max_nodes = nodes[(slice(1, None), )*self.dim] + max_nodes = nodes[(slice(1, None),) * self.dim] max_nodes = max_nodes.reshape((self.n_cells, -1), order="F") - + cell_bounds = np.stack((min_nodes, max_nodes), axis=-1) cell_bounds = cell_bounds.reshape((self.n_cells, -1)) return cell_bounds