Skip to content

Renamed nx to dof when nx represents the degree of freedom of an interpolator, to avoid confusion with nx the 1st coordinate of a normal #250

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
merged 6 commits into from
Jul 1, 2025
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
28 changes: 14 additions & 14 deletions LoopStructural/interpolators/_discrete_interpolator.py
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand All @@ -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:
Expand All @@ -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
Expand Down Expand Up @@ -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
)
)

Expand Down Expand Up @@ -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"):
Expand Down Expand Up @@ -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,
Expand Down Expand Up @@ -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
Expand All @@ -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,
}
Expand Down Expand Up @@ -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"):
Expand All @@ -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)

Expand Down Expand Up @@ -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
Expand All @@ -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()

Expand All @@ -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
Expand Down
20 changes: 10 additions & 10 deletions LoopStructural/interpolators/_finite_difference_interpolator.py
Original file line number Diff line number Diff line change
Expand Up @@ -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, :]]
Expand Down Expand Up @@ -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, :]]
Expand All @@ -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(
Expand Down Expand Up @@ -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, :]]
Expand Down Expand Up @@ -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, :]]
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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()

Expand All @@ -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]))

Expand All @@ -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
Expand Down
2 changes: 1 addition & 1 deletion LoopStructural/interpolators/_geological_interpolator.py
Original file line number Diff line number Diff line change
Expand Up @@ -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

Expand Down
2 changes: 1 addition & 1 deletion LoopStructural/interpolators/_surfe_wrapper.py
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down
6 changes: 3 additions & 3 deletions LoopStructural/modelling/core/geological_model.py
Original file line number Diff line number Diff line change
Expand Up @@ -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(
Expand Down
Original file line number Diff line number Diff line change
Expand Up @@ -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):
Expand Down
2 changes: 1 addition & 1 deletion tests/unit/interpolator/test_discrete_interpolator.py
Original file line number Diff line number Diff line change
@@ -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):
Expand Down