Merge branch 'main' into 1.0.0_alpha

Author: jcapriot

.ci/azure/test.yml

@@ -12,6 +12,9 @@ jobs:
       linux-Python312:
         image: ubuntu-latest
         python.version: '3.12'
+      linux-Python313:
+        image: ubuntu-latest
+        python.version: '3.13'
       osx-Python310:
         image: macOS-latest
         python.version: '3.10'
@@ -21,6 +24,9 @@ jobs:
       osx-Python312:
         image: macOS-latest
         python.version: '3.12'
+      osx-Python313:
+        image: macOS-latest
+        python.version: '3.13'
       win-Python310:
         image: windows-latest
         python.version: '3.10'
@@ -30,6 +36,9 @@ jobs:
       win-Python312:
         image: windows-latest
         python.version: '3.12'
+      win-Python313:
+        image: windows-latest
+        python.version: '3.13'
   displayName: "${{ variables.image }} ${{ variables.python.version }}"
   pool:
     vmImage: $(image)

LICENSE

@@ -1,6 +1,6 @@
 The MIT License (MIT)
 
-Copyright (c) 2013-2024 SimPEG Developers
+Copyright (c) 2013-2025 SimPEG Developers
 
 Permission is hereby granted, free of charge, to any person obtaining a copy of
 this software and associated documentation files (the "Software"), to deal in

discretize/_extensions/tree_ext.pyx

@@ -6,7 +6,7 @@ cimport numpy as np
 from libc.stdlib cimport malloc, free
 from libcpp.vector cimport vector
 from libcpp cimport bool
-from numpy.math cimport INFINITY
+from libc.math cimport INFINITY
 
 from .tree cimport int_t, Tree as c_Tree, PyWrapper, Node, Edge, Face, Cell as c_Cell
 from . cimport geom
@@ -15,6 +15,9 @@ import scipy.sparse as sp
 import numpy as np
 from .interputils_cython cimport _bisect_left, _bisect_right
 
+class TreeMeshNotFinalizedError(RuntimeError):
+    """Raise when a TreeMesh is not finalized."""
+
 
 cdef class TreeCell:
     """A class for defining cells within instances of :class:`~discretize.TreeMesh`.
@@ -1959,6 +1962,7 @@ cdef class _TreeMesh:
         (n_cells, dim) numpy.ndarray of float
             Gridded cell center locations
         """
+        self._error_if_not_finalized("cell_centers")
         cdef np.float64_t[:, :] gridCC
         cdef np.int64_t ii, ind, dim
         if self._cell_centers is None:
@@ -1984,6 +1988,7 @@ cdef class _TreeMesh:
         (n_nodes, dim) numpy.ndarray of float
             Gridded non-hanging node locations
         """
+        self._error_if_not_finalized("nodes")
         cdef np.float64_t[:, :] gridN
         cdef Node *node
         cdef np.int64_t ii, ind, dim
@@ -2012,6 +2017,7 @@ cdef class _TreeMesh:
         (n_hanging_nodes, dim) numpy.ndarray of float
             Gridded hanging node locations
         """
+        self._error_if_not_finalized("hanging_nodes")
         cdef np.float64_t[:, :] gridN
         cdef Node *node
         cdef np.int64_t ii, ind, dim
@@ -2038,6 +2044,7 @@ cdef class _TreeMesh:
         (n_boundary_nodes, dim) numpy.ndarray of float
             Gridded boundary node locations
         """
+        self._error_if_not_finalized("boundary_nodes")
         nodes = self.nodes
         x0, xF = self._xs[0], self._xs[-1]
         y0, yF = self._ys[0], self._ys[-1]
@@ -2069,6 +2076,7 @@ cdef class _TreeMesh:
         (n_cells, dim) numpy.ndarray of float
             Gridded cell dimensions
         """
+        self._error_if_not_finalized("h_gridded")
         if self._h_gridded is not None:
             return self._h_gridded
         cdef np.float64_t[:, :] gridCH
@@ -2097,6 +2105,7 @@ cdef class _TreeMesh:
         (n_edges_x, dim) numpy.ndarray of float
             Gridded locations of all non-hanging x-edges
         """
+        self._error_if_not_finalized("edges_x")
         cdef np.float64_t[:, :] gridEx
         cdef Edge *edge
         cdef np.int64_t ii, ind, dim
@@ -2124,6 +2133,7 @@ cdef class _TreeMesh:
         (n_hanging_edges_x, dim) numpy.ndarray of float
             Gridded locations of all hanging x-edges
         """
+        self._error_if_not_finalized("hanging_edges_x")
         cdef np.float64_t[:, :] gridhEx
         cdef Edge *edge
         cdef np.int64_t ii, ind, dim
@@ -2149,6 +2159,7 @@ cdef class _TreeMesh:
         (n_edges_y, dim) numpy.ndarray of float
             Gridded locations of all non-hanging y-edges
         """
+        self._error_if_not_finalized("edges_y")
         cdef np.float64_t[:, :] gridEy
         cdef Edge *edge
         cdef np.int64_t ii, ind, dim
@@ -2176,6 +2187,7 @@ cdef class _TreeMesh:
         (n_haning_edges_y, dim) numpy.ndarray of float
             Gridded locations of all hanging y-edges
         """
+        self._error_if_not_finalized("hanging_edges_y")
         cdef np.float64_t[:, :] gridhEy
         cdef Edge *edge
         cdef np.int64_t ii, ind, dim
@@ -2201,6 +2213,7 @@ cdef class _TreeMesh:
         (n_edges_z, dim) numpy.ndarray of float
             Gridded locations of all non-hanging z-edges
         """
+        self._error_if_not_finalized("edges_z")
         cdef np.float64_t[:, :] gridEz
         cdef Edge *edge
         cdef np.int64_t ii, ind, dim
@@ -2228,6 +2241,7 @@ cdef class _TreeMesh:
         (n_hanging_edges_z, dim) numpy.ndarray of float
             Gridded locations of all hanging z-edges
         """
+        self._error_if_not_finalized("hanging_edges_z")
         cdef np.float64_t[:, :] gridhEz
         cdef Edge *edge
         cdef np.int64_t ii, ind, dim
@@ -2255,6 +2269,7 @@ cdef class _TreeMesh:
         (n_boundary_edges, dim) numpy.ndarray of float
             Gridded boundary edge locations
         """
+        self._error_if_not_finalized("boundary_edges")
         edges_x = self.edges_x
         edges_y = self.edges_y
         x0, xF = self._xs[0], self._xs[-1]
@@ -2294,6 +2309,7 @@ cdef class _TreeMesh:
         (n_faces_x, dim) numpy.ndarray of float
             Gridded locations of all non-hanging x-faces
         """
+        self._error_if_not_finalized("faces_x")
         if(self._dim == 2): return self.edges_y
 
         cdef np.float64_t[:, :] gridFx
@@ -2323,6 +2339,7 @@ cdef class _TreeMesh:
         (n_faces_y, dim) numpy.ndarray of float
             Gridded locations of all non-hanging y-faces
         """
+        self._error_if_not_finalized("faces_y")
         if(self._dim == 2): return self.edges_x
         cdef np.float64_t[:, :] gridFy
         cdef Face *face
@@ -2351,6 +2368,7 @@ cdef class _TreeMesh:
         (n_faces_z, dim) numpy.ndarray of float
             Gridded locations of all non-hanging z-faces
         """
+        self._error_if_not_finalized("faces_z")
         if(self._dim == 2): return self.cell_centers
 
         cdef np.float64_t[:, :] gridFz
@@ -2380,6 +2398,7 @@ cdef class _TreeMesh:
         (n_hanging_faces_x, dim) numpy.ndarray of float
             Gridded locations of all hanging x-faces
         """
+        self._error_if_not_finalized("hanging_faces_x")
         if(self._dim == 2): return self.hanging_edges_y
 
         cdef np.float64_t[:, :] gridFx
@@ -2407,6 +2426,7 @@ cdef class _TreeMesh:
         (n_hanging_faces_y, dim) numpy.ndarray of float
             Gridded locations of all hanging y-faces
         """
+        self._error_if_not_finalized("hanging_faces_y")
         if(self._dim == 2): return self.hanging_edges_x
 
         cdef np.float64_t[:, :] gridhFy
@@ -2434,6 +2454,7 @@ cdef class _TreeMesh:
         (n_hanging_faces_z, dim) numpy.ndarray of float
             Gridded locations of all hanging z-faces
         """
+        self._error_if_not_finalized("hanging_faces_z")
         if(self._dim == 2): return np.array([])
 
         cdef np.float64_t[:, :] gridhFz
@@ -2463,6 +2484,7 @@ cdef class _TreeMesh:
         (n_boundary_faces, dim) numpy.ndarray of float
             Gridded boundary face locations
         """
+        self._error_if_not_finalized("boundary_faces")
         faces_x = self.faces_x
         faces_y = self.faces_y
         x0, xF = self._xs[0], self._xs[-1]
@@ -2492,6 +2514,7 @@ cdef class _TreeMesh:
         (n_boundary_faces, dim) numpy.ndarray of float
             Outward normals of boundary faces
         """
+        self._error_if_not_finalized("boundary_face_outward_normals")
         faces_x = self.faces_x
         faces_y = self.faces_y
         x0, xF = self._xs[0], self._xs[-1]
@@ -2534,6 +2557,7 @@ cdef class _TreeMesh:
               - *3D:* Returns the cell volumes
 
         """
+        self._error_if_not_finalized("cell_volumes")
         cdef np.float64_t[:] vol
         if self._cell_volumes is None:
             self._cell_volumes = np.empty(self.n_cells, dtype=np.float64)
@@ -2558,6 +2582,7 @@ cdef class _TreeMesh:
               respectively
             - *3D:* returns the x, y and z-face areas in order
         """
+        self._error_if_not_finalized("face_areas")
         if self._dim == 2 and self._face_areas is None:
             self._face_areas = np.r_[self.edge_lengths[self.n_edges_x:], self.edge_lengths[:self.n_edges_x]]
         cdef np.float64_t[:] area
@@ -2600,6 +2625,7 @@ cdef class _TreeMesh:
             - *2D:* returns the x-edge and y-edge lengths in order
             - *3D:* returns the x, y and z-edge lengths in order
         """
+        self._error_if_not_finalized("edge_lengths")
         cdef np.float64_t[:] edge_l
         cdef Edge *edge
         cdef int_t ind, offset
@@ -3564,6 +3590,7 @@ cdef class _TreeMesh:
         (n_total_faces_x, n_cells) scipy.sparse.csr_matrix
             The stencil for the x-component of the cell gradient
         """
+        self._error_if_not_finalized("stencil_cell_gradient_x")
         if getattr(self, '_stencil_cell_gradient_x', None) is not None:
             return self._stencil_cell_gradient_x
         cdef np.int64_t[:] I = np.zeros(2*self.n_total_faces_x, dtype=np.int64)
@@ -3641,6 +3668,7 @@ cdef class _TreeMesh:
         (n_total_faces_y, n_cells) scipy.sparse.csr_matrix
             The stencil for the y-component of the cell gradient
         """
+        self._error_if_not_finalized("stencil_cell_gradient_y")
         if getattr(self, '_stencil_cell_gradient_y', None) is not None:
             return self._stencil_cell_gradient_y
 
@@ -3719,6 +3747,7 @@ cdef class _TreeMesh:
         (n_total_faces_z, n_cells) scipy.sparse.csr_matrix
             The stencil for the z-component of the cell gradient
         """
+        self._error_if_not_finalized("stencil_cell_gradient_z")
         if getattr(self, '_stencil_cell_gradient_z', None) is not None:
             return self._stencil_cell_gradient_z
 
@@ -5019,6 +5048,7 @@ cdef class _TreeMesh:
 
             phi_f = Acf @ phi_c
         """
+        self._error_if_not_finalized("average_cell_to_face")
         if self._average_cell_to_face is not None:
             return self._average_cell_to_face
         stacks = [self.average_cell_to_face_x, self.average_cell_to_face_y]
@@ -5079,6 +5109,7 @@ cdef class _TreeMesh:
         in the x-direction. For boundary faces, nearest neighbor is used to extrapolate
         the values.
         """
+        self._error_if_not_finalized("average_cell_vector_to_face")
         if self._average_cell_vector_to_face is not None:
             return self._average_cell_vector_to_face
         stacks = [self.average_cell_to_face_x, self.average_cell_to_face_y]
@@ -5102,6 +5133,7 @@ cdef class _TreeMesh:
         (n_faces_x, n_cells) scipy.sparse.csr_matrix
             The scalar averaging operator from cell centers to x faces
         """
+        self._error_if_not_finalized("average_cell_to_face_x")
         if self._average_cell_to_face_x is not None:
             return self._average_cell_to_face_x
         cdef np.int64_t[:] I = np.zeros(2*self.n_total_faces_x, dtype=np.int64)
@@ -5218,6 +5250,7 @@ cdef class _TreeMesh:
         (n_faces_y, n_cells) scipy.sparse.csr_matrix
             The scalar averaging operator from cell centers to y faces
         """
+        self._error_if_not_finalized("average_cell_to_face_y")
         if self._average_cell_to_face_y is not None:
             return self._average_cell_to_face_y
         cdef np.int64_t[:] I = np.zeros(2*self.n_total_faces_y, dtype=np.int64)
@@ -5334,6 +5367,7 @@ cdef class _TreeMesh:
         (n_faces_z, n_cells) scipy.sparse.csr_matrix
             The scalar averaging operator from cell centers to z faces
         """
+        self._error_if_not_finalized("average_cell_to_face_z")
         if self.dim == 2:
             raise Exception('TreeMesh has no z-faces in 2D')
         if self._average_cell_to_face_z is not None:
@@ -5424,6 +5458,7 @@ cdef class _TreeMesh:
         scipy.sparse.csr_matrix
             (n_boundary_faces, n_faces) Projection matrix with shape
         """
+        self._error_if_not_finalized("project_face_to_boundary_face")
         faces_x = self.faces_x
         faces_y = self.faces_y
 
@@ -5459,6 +5494,7 @@ cdef class _TreeMesh:
         (n_boundary_edges, n_edges) scipy.sparse.csr_matrix
             Projection matrix with shape
         """
+        self._error_if_not_finalized("project_edge_to_boundary_edge")
         edges_x = self.edges_x
         edges_y = self.edges_y
 
@@ -5501,6 +5537,7 @@ cdef class _TreeMesh:
         (n_boundary_nodes, n_nodes) scipy.sparse.csr_matrix
             Projection matrix with shape
         """
+        self._error_if_not_finalized("project_node_to_boundary_node")
         nodes = self.nodes
         x0, xF = self._xs[0], self._xs[-1]
         y0, yF = self._ys[0], self._ys[-1]
@@ -6964,6 +7001,17 @@ cdef class _TreeMesh:
         """
         return self.get_cells_in_aabb(*rectangle.reshape(self.dim, 2).T)
 
+    def _error_if_not_finalized(self, method: str):
+        """
+        Raise error if mesh is not finalized.
+        """
+        if not self.finalized:
+            msg = (
+                f"`{type(self).__name__}.{method}` requires a finalized mesh. "
+                "Use the `finalize()` method to finalize it."
+            )
+            raise TreeMeshNotFinalizedError(msg)
+
     def _require_ndarray_with_dim(self, name, arr, ndim=1, dtype=None, requirements=None):
         """Returns an ndarray that has dim along it's last dimension, with ndim dims,
 

discretize/mixins/mesh_io.py

@@ -466,7 +466,7 @@ def read_UBC(TreeMesh, file_name, directory=None):
         else:
             max_level = min(ls) + 1
 
-        mesh = TreeMesh(hs, origin=origin)
+        mesh = TreeMesh(hs, origin=origin, diagonal_balance=False)
         levels = indArr[:, -1]
         indArr = indArr[:, :-1]
 

discretize/mixins/mpl_mod.py

@@ -2053,7 +2053,7 @@ def __plot_slice_tree(
         normal[normalInd] = 1
 
         # create a temporary TreeMesh with the slice through
-        temp_mesh = discretize.TreeMesh(h2d, x2d)
+        temp_mesh = discretize.TreeMesh(h2d, x2d, diagonal_balance=False)
         level_diff = self.max_level - temp_mesh.max_level
 
         # get list of cells which intersect the slicing plane