From 4a493f6e1d3cfd9cc762d2a9b461179997534090 Mon Sep 17 00:00:00 2001 From: AxMeNi <159522803+AxMeNi@users.noreply.github.com> Date: Mon, 16 Jun 2025 09:50:52 +1000 Subject: [PATCH 1/4] Modified str method of _bounding_box for clarity --- LoopStructural/datatypes/_bounding_box.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/LoopStructural/datatypes/_bounding_box.py b/LoopStructural/datatypes/_bounding_box.py index a2e1eb25..c4ea65eb 100644 --- a/LoopStructural/datatypes/_bounding_box.py +++ b/LoopStructural/datatypes/_bounding_box.py @@ -547,10 +547,10 @@ def reproject(self, xyz, inplace=False): return xyz + self.global_origin def __repr__(self): - return f"BoundingBox({self.origin}, {self.maximum}, {self.nsteps})" + return f"BoundingBox(origin:{self.origin}, maximum:{self.maximum}, nsteps:{self.nsteps})" def __str__(self): - return f"BoundingBox({self.origin}, {self.maximum}, {self.nsteps})" + return f"BoundingBox(origin:{self.origin}, maximum:{self.maximum}, nsteps:{self.nsteps})" def __eq__(self, other): if not isinstance(other, BoundingBox): From 0c1c16cce88580d1168daf141f7230ee5ceee23b Mon Sep 17 00:00:00 2001 From: AxMeNi <159522803+AxMeNi@users.noreply.github.com> Date: Wed, 18 Jun 2025 14:23:34 +1000 Subject: [PATCH 2/4] Update _3d_base_structured.py --- .../supports/_3d_base_structured.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/LoopStructural/interpolators/supports/_3d_base_structured.py b/LoopStructural/interpolators/supports/_3d_base_structured.py index 4c91322c..621da7be 100644 --- a/LoopStructural/interpolators/supports/_3d_base_structured.py +++ b/LoopStructural/interpolators/supports/_3d_base_structured.py @@ -373,6 +373,22 @@ def cell_corner_indexes(self, cell_indexes: np.ndarray) -> np.ndarray: return corner_indexes def position_to_cell_corners(self, pos): + """Get the global indices of the vertices (corners) of the cell containing each point. + + Parameters + ---------- + pos : np.array + (N, 3) array of xyz coordinates representing the positions of N points. + + Returns + ------- + globalidx : np.array + (N, 8) array of global indices corresponding to the 8 corner nodes of the cell + each point lies in. If a point lies outside the support, its corresponding entry + will be set to -1. + inside : np.array + (N,) boolean array indicating whether each point is inside the support domain. + """ cell_indexes, inside = self.position_to_cell_index(pos) corner_indexes = self.cell_corner_indexes(cell_indexes) globalidx = self.global_node_indices(corner_indexes) From 22f73fffa97017a126ee3e2c9a91c6534c5267a6 Mon Sep 17 00:00:00 2001 From: AxMeNi <159522803+AxMeNi@users.noreply.github.com> Date: Wed, 18 Jun 2025 14:25:19 +1000 Subject: [PATCH 3/4] Update _2d_structured_grid.py --- .../supports/_2d_structured_grid.py | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/LoopStructural/interpolators/supports/_2d_structured_grid.py b/LoopStructural/interpolators/supports/_2d_structured_grid.py index 61e8d825..3731c306 100644 --- a/LoopStructural/interpolators/supports/_2d_structured_grid.py +++ b/LoopStructural/interpolators/supports/_2d_structured_grid.py @@ -372,6 +372,22 @@ def node_indexes_to_position(self, node_indexes: np.ndarray) -> np.ndarray: return xy def position_to_cell_corners(self, pos): + """Get the global indices of the vertices (corner) nodes of the cell containing each point. + + Parameters + ---------- + pos : np.array + (N, 2) array of xy coordinates representing the positions of N points. + + Returns + ------- + globalidx : np.array + (N, 4) array of global indices corresponding to the 4 corner nodes of the cell + each point lies in. If a point lies outside the support, its corresponding entry + will be set to -1. + inside : np.array + (N,) boolean array indicating whether each point is inside the support domain. + """ corner_index, inside = self.position_to_cell_index(pos) corners = self.cell_corner_indexes(corner_index) globalidx = self.global_node_indices(corners) From ab989cb9b7a22e339c945483f6316493c2f1fbb3 Mon Sep 17 00:00:00 2001 From: Axel Mengelle-Nicole Date: Thu, 19 Jun 2025 14:01:22 +1000 Subject: [PATCH 4/4] Renamed nx to dof when it represents the degree of freedom of an interpolator --- .../interpolators/_discrete_interpolator.py | 28 +++++++++---------- .../_finite_difference_interpolator.py | 20 ++++++------- .../interpolators/_geological_interpolator.py | 2 +- .../interpolators/_surfe_wrapper.py | 2 +- .../supports/_2d_base_unstructured.py | 2 +- .../modelling/core/geological_model.py | 6 ++-- .../builders/_geological_feature_builder.py | 2 +- .../test_discrete_interpolator.py | 2 +- 8 files changed, 32 insertions(+), 32 deletions(-) diff --git a/LoopStructural/interpolators/_discrete_interpolator.py b/LoopStructural/interpolators/_discrete_interpolator.py index 60f19932..45ac63ab 100644 --- a/LoopStructural/interpolators/_discrete_interpolator.py +++ b/LoopStructural/interpolators/_discrete_interpolator.py @@ -42,7 +42,7 @@ def __init__(self, support, data={}, c=None, up_to_date=False): self.shape = "rectangular" if self.shape == "square": - self.B = np.zeros(self.nx) + self.B = np.zeros(self.dof) self.c_ = 0 self.solver = "cg" @@ -60,7 +60,7 @@ def __init__(self, support, data={}, c=None, up_to_date=False): self.non_linear_constraints = [] self.constraints = {} self.interpolation_weights = {} - logger.info("Creating discrete interpolator with {} degrees of freedom".format(self.nx)) + logger.info("Creating discrete interpolator with {} degrees of freedom".format(self.dof)) self.type = InterpolatorType.BASE_DISCRETE def set_nelements(self, nelements: int) -> int: @@ -78,7 +78,7 @@ def n_elements(self) -> int: return self.support.n_elements @property - def nx(self) -> int: + def dof(self) -> int: """Number of degrees of freedom for the interpolator Returns @@ -125,7 +125,7 @@ def set_region(self, region=None): # self.region_function = region logger.info( "Cannot use region at the moment. Interpolation now uses region and has {} degrees of freedom".format( - self.nx + self.dof ) ) @@ -175,7 +175,7 @@ def reset(self): """ self.constraints = {} self.c_ = 0 - self.regularisation_scale = np.ones(self.nx) + self.regularisation_scale = np.ones(self.dof) logger.info("Resetting interpolation constraints") def add_constraints_to_least_squares(self, A, B, idc, w=1.0, name="undefined"): @@ -245,7 +245,7 @@ def add_constraints_to_least_squares(self, A, B, idc, w=1.0, name="undefined"): rows = np.tile(rows, (A.shape[-1], 1)).T self.constraints[name] = { 'matrix': sparse.coo_matrix( - (A.flatten(), (rows.flatten(), idc.flatten())), shape=(n_rows, self.nx) + (A.flatten(), (rows.flatten(), idc.flatten())), shape=(n_rows, self.dof) ).tocsc(), 'b': B.flatten(), 'w': w, @@ -286,7 +286,7 @@ def add_inequality_constraints_to_matrix( A : numpy array matrix of coefficients bounds : numpy array - nx3 lower, upper, 1 + n*3 lower, upper, 1 idc : numpy array index of constraints in the matrix Returns @@ -296,14 +296,14 @@ def add_inequality_constraints_to_matrix( # map from mesh node index to region node index gi = np.zeros(self.support.n_nodes, dtype=int) gi[:] = -1 - gi[self.region] = np.arange(0, self.nx, dtype=int) + gi[self.region] = np.arange(0, self.dof, dtype=int) idc = gi[idc] rows = np.arange(0, idc.shape[0]) rows = np.tile(rows, (A.shape[-1], 1)).T self.ineq_constraints[name] = { 'matrix': sparse.coo_matrix( - (A.flatten(), (rows.flatten(), idc.flatten())), shape=(rows.shape[0], self.nx) + (A.flatten(), (rows.flatten(), idc.flatten())), shape=(rows.shape[0], self.dof) ).tocsc(), "bounds": bounds, } @@ -431,7 +431,7 @@ def add_inequality_feature( np.ones((value.shape[0], 1)), l, u, - np.arange(0, self.nx, dtype=int), + np.arange(0, self.dof, dtype=int), ) def add_equality_constraints(self, node_idx, values, name="undefined"): @@ -454,7 +454,7 @@ def add_equality_constraints(self, node_idx, values, name="undefined"): # map from mesh node index to region node index gi = np.zeros(self.support.n_nodes) gi[:] = -1 - gi[self.region] = np.arange(0, self.nx) + gi[self.region] = np.arange(0, self.dof) idc = gi[node_idx] outside = ~(idc == -1) @@ -515,7 +515,7 @@ def add_equality_block(self, A, B): if len(self.equal_constraints) > 0: ATA = A.T.dot(A) ATB = A.T.dot(B) - logger.info(f"Equality block is {self.eq_const_c} x {self.nx}") + logger.info(f"Equality block is {self.eq_const_c} x {self.dof}") # solving constrained least squares using # | ATA CT | |c| = b # | C 0 | |y| d @@ -540,7 +540,7 @@ def add_equality_block(self, A, B): C = sparse.coo_matrix( (np.array(a), (np.array(rows), cols)), - shape=(self.eq_const_c, self.nx), + shape=(self.eq_const_c, self.dof), dtype=float, ).tocsr() @@ -557,7 +557,7 @@ def build_inequality_matrix(self): mats.append(c['matrix']) bounds.append(c['bounds']) if len(mats) == 0: - return sparse.csr_matrix((0, self.nx), dtype=float), np.zeros((0, 3)) + return sparse.csr_matrix((0, self.dof), dtype=float), np.zeros((0, 3)) Q = sparse.vstack(mats) bounds = np.vstack(bounds) return Q, bounds diff --git a/LoopStructural/interpolators/_finite_difference_interpolator.py b/LoopStructural/interpolators/_finite_difference_interpolator.py index e51faa2d..84e1f59f 100644 --- a/LoopStructural/interpolators/_finite_difference_interpolator.py +++ b/LoopStructural/interpolators/_finite_difference_interpolator.py @@ -160,7 +160,7 @@ def add_value_constraints(self, w=1.0): # print(points[inside,:].shape) gi = np.zeros(self.support.n_nodes, dtype=int) gi[:] = -1 - gi[self.region] = np.arange(0, self.nx, dtype=int) + gi[self.region] = np.arange(0, self.dof, dtype=int) idc = np.zeros(node_idx.shape) idc[:] = -1 idc[inside, :] = gi[node_idx[inside, :]] @@ -204,7 +204,7 @@ def add_interface_constraints(self, w=1.0): ) gi = np.zeros(self.support.n_nodes, dtype=int) gi[:] = -1 - gi[self.region] = np.arange(0, self.nx, dtype=int) + gi[self.region] = np.arange(0, self.dof, dtype=int) idc = np.zeros(node_idx.shape).astype(int) idc[:] = -1 idc[inside, :] = gi[node_idx[inside, :]] @@ -227,11 +227,11 @@ def add_interface_constraints(self, w=1.0): interface_A = np.hstack([A[mask, :][ij[:, 0], :], -A[mask, :][ij[:, 1], :]]) interface_idc = np.hstack([idc[mask, :][ij[:, 0], :], idc[mask, :][ij[:, 1], :]]) # now map the index from global to region create array size of mesh - # initialise as np.nan, then map points inside region to 0->nx + # initialise as np.nan, then map points inside region to 0->dof gi = np.zeros(self.support.n_nodes).astype(int) gi[:] = -1 - gi[self.region] = np.arange(0, self.nx) + gi[self.region] = np.arange(0, self.dof) interface_idc = gi[interface_idc] outside = ~np.any(interface_idc == -1, axis=1) self.add_constraints_to_least_squares( @@ -266,7 +266,7 @@ def add_gradient_constraints(self, w=1.0): # magnitude gi = np.zeros(self.support.n_nodes) gi[:] = -1 - gi[self.region] = np.arange(0, self.nx) + gi[self.region] = np.arange(0, self.dof) idc = np.zeros(node_idx.shape) idc[:] = -1 idc[inside, :] = gi[node_idx[inside, :]] @@ -331,7 +331,7 @@ def add_norm_constraints(self, w=1.0): ) gi = np.zeros(self.support.n_nodes) gi[:] = -1 - gi[self.region] = np.arange(0, self.nx) + gi[self.region] = np.arange(0, self.dof) idc = np.zeros(node_idx.shape) idc[:] = -1 idc[inside, :] = gi[node_idx[inside, :]] @@ -421,7 +421,7 @@ def add_gradient_orthogonal_constraints( # magnitude gi = np.zeros(self.support.n_nodes) gi[:] = -1 - gi[self.region] = np.arange(0, self.nx) + gi[self.region] = np.arange(0, self.dof) idc = np.zeros(node_idx.shape) idc[:] = -1 @@ -472,7 +472,7 @@ def add_gradient_orthogonal_constraints( # """ # # First get the global indicies of the pairs of neighbours this should be an - # # Nx27 array for 3d and an Nx9 array for 2d + # # N*27 array for 3d and an N*9 array for 2d # global_indexes = self.support.neighbour_global_indexes() @@ -490,7 +490,7 @@ def assemble_inner(self, operator, w, name='regularisation'): """ # First get the global indicies of the pairs of neighbours this should be an - # Nx27 array for 3d and an Nx9 array for 2d + # N*27 array for 3d and an N*9 array for 2d global_indexes = self.support.neighbour_global_indexes() # np.array([ii,jj])) @@ -499,7 +499,7 @@ def assemble_inner(self, operator, w, name='regularisation'): gi = np.zeros(self.support.n_nodes) gi[:] = -1 - gi[self.region] = np.arange(0, self.nx) + gi[self.region] = np.arange(0, self.dof) idc = gi[idc] inside = ~np.any(idc == -1, axis=1) # np.ones(a.shape[0],dtype=bool)# # a[idc==-1] = 0 diff --git a/LoopStructural/interpolators/_geological_interpolator.py b/LoopStructural/interpolators/_geological_interpolator.py index 840faf6a..0368698c 100644 --- a/LoopStructural/interpolators/_geological_interpolator.py +++ b/LoopStructural/interpolators/_geological_interpolator.py @@ -96,7 +96,7 @@ def to_json(self): json["constraints"] = self.constraints json["data"] = self.data json["type"] = self.type - # json["dof"] = self.nx + # json["dof"] = self.dof json["up_to_date"] = self.up_to_date return json diff --git a/LoopStructural/interpolators/_surfe_wrapper.py b/LoopStructural/interpolators/_surfe_wrapper.py index 4784434a..6b7bef4b 100644 --- a/LoopStructural/interpolators/_surfe_wrapper.py +++ b/LoopStructural/interpolators/_surfe_wrapper.py @@ -203,7 +203,7 @@ def evaluate_gradient(self, evaluation_points): return @property - def nx(self): + def dof(self): return self.get_data_locations().shape[0] @property def n_elements(self)->int: diff --git a/LoopStructural/interpolators/supports/_2d_base_unstructured.py b/LoopStructural/interpolators/supports/_2d_base_unstructured.py index 82868d0b..d0d53e0e 100644 --- a/LoopStructural/interpolators/supports/_2d_base_unstructured.py +++ b/LoopStructural/interpolators/supports/_2d_base_unstructured.py @@ -30,7 +30,7 @@ def __init__(self, elements, vertices, neighbours, aabb_nsteps=None): self.order = 1 elif self.elements.shape[1] == 6: self.order = 2 - self.nx = self.vertices.shape[0] + self.dof = self.vertices.shape[0] self.neighbours = neighbours self.minimum = np.min(self.nodes, axis=0) self.maximum = np.max(self.nodes, axis=0) diff --git a/LoopStructural/modelling/core/geological_model.py b/LoopStructural/modelling/core/geological_model.py index e8a10a8f..e2b210b8 100644 --- a/LoopStructural/modelling/core/geological_model.py +++ b/LoopStructural/modelling/core/geological_model.py @@ -1670,15 +1670,15 @@ def update(self, verbose=False, progressbar=True): for f in self.features: if f.type == FeatureType.FAULT: nfeatures += 3 - total_dof += f[0].interpolator.nx * 3 + total_dof += f[0].interpolator.dof * 3 continue if isinstance(f, StructuralFrame): nfeatures += 3 - total_dof += f[0].interpolator.nx * 3 + total_dof += f[0].interpolator.dof * 3 continue if f.type == FeatureType.INTERPOLATED: nfeatures += 1 - total_dof += f.interpolator.nx + total_dof += f.interpolator.dof continue if verbose: print( diff --git a/LoopStructural/modelling/features/builders/_geological_feature_builder.py b/LoopStructural/modelling/features/builders/_geological_feature_builder.py index 937068c9..9cd1b3c5 100644 --- a/LoopStructural/modelling/features/builders/_geological_feature_builder.py +++ b/LoopStructural/modelling/features/builders/_geological_feature_builder.py @@ -452,7 +452,7 @@ def set_interpolation_geometry(self, origin, maximum, rotation=None): self.interpolator.support.rotation_xy = rotation self._up_to_date = False - while self.interpolator.nx < 100: + while self.interpolator.dof < 100: self.interpolator.support.step_vector = self.interpolator.support.step_vector * 0.9 def check_interpolation_geometry(self, data, buffer=0.3): diff --git a/tests/unit/interpolator/test_discrete_interpolator.py b/tests/unit/interpolator/test_discrete_interpolator.py index 981b1bbe..8aa84f62 100644 --- a/tests/unit/interpolator/test_discrete_interpolator.py +++ b/tests/unit/interpolator/test_discrete_interpolator.py @@ -1,5 +1,5 @@ def test_nx(interpolator, data): - assert interpolator.nx == 21 * 21 * 21 + assert interpolator.dof == 21 * 21 * 21 def test_region(interpolator, data, region_func):