diff --git a/LoopStructural/datatypes/_point.py b/LoopStructural/datatypes/_point.py index 4542262d..c4f11ae3 100644 --- a/LoopStructural/datatypes/_point.py +++ b/LoopStructural/datatypes/_point.py @@ -157,7 +157,6 @@ def vtk( try: locations = bb.project(locations) _projected = True - scale = bb.scale_by_projection_factor(scale) except Exception as e: logger.error(f'Failed to project points to bounding box: {e}') logger.error('Using unprojected points, this may cause issues with the glyphing') diff --git a/LoopStructural/interpolators/__init__.py b/LoopStructural/interpolators/__init__.py index 56e68480..401a00b8 100644 --- a/LoopStructural/interpolators/__init__.py +++ b/LoopStructural/interpolators/__init__.py @@ -3,6 +3,7 @@ """ + __all__ = [ "InterpolatorType", "GeologicalInterpolator", @@ -21,39 +22,12 @@ "StructuredGrid2D", "P2UnstructuredTetMesh", ] -from enum import IntEnum +from ._interpolatortype import InterpolatorType from ..utils import getLogger logger = getLogger(__name__) - -class InterpolatorType(IntEnum): - """ - Enum for the different interpolator types - - 1-9 should cover interpolators with supports - 9+ are data supported - """ - - BASE = 0 - BASE_DISCRETE = 1 - FINITE_DIFFERENCE = 2 - DISCRETE_FOLD = 3 - PIECEWISE_LINEAR = 4 - PIECEWISE_QUADRATIC = 5 - BASE_DATA_SUPPORTED = 10 - SURFE = 11 - - -interpolator_string_map = { - "FDI": InterpolatorType.FINITE_DIFFERENCE, - "PLI": InterpolatorType.PIECEWISE_LINEAR, - "P2": InterpolatorType.PIECEWISE_QUADRATIC, - "P1": InterpolatorType.PIECEWISE_LINEAR, - "DFI": InterpolatorType.DISCRETE_FOLD, - 'surfe': InterpolatorType.SURFE, -} from ..interpolators._geological_interpolator import GeologicalInterpolator from ..interpolators._discrete_interpolator import DiscreteInterpolator from ..interpolators.supports import ( @@ -79,7 +53,7 @@ class InterpolatorType(IntEnum): ) from ..interpolators._p2interpolator import P2Interpolator from ..interpolators._p1interpolator import P1Interpolator - +from ..interpolators._constant_norm import ConstantNormP1Interpolator, ConstantNormFDIInterpolator try: from ..interpolators._surfe_wrapper import SurfeRBFInterpolator except ImportError: @@ -93,6 +67,24 @@ def __init__(self, *args, **kwargs): raise ImportError( "Surfe cannot be imported. Please install Surfe. pip install surfe/ conda install -c loop3d surfe" ) + +# Ensure compatibility between the fallback and imported class +SurfeRBFInterpolator = SurfeRBFInterpolator + + +interpolator_string_map = { + "FDI": InterpolatorType.FINITE_DIFFERENCE, + "PLI": InterpolatorType.PIECEWISE_LINEAR, + "P2": InterpolatorType.PIECEWISE_QUADRATIC, + "P1": InterpolatorType.PIECEWISE_LINEAR, + "DFI": InterpolatorType.DISCRETE_FOLD, + 'surfe': InterpolatorType.SURFE, + "FDI_CN": InterpolatorType.FINITE_DIFFERENCE_CONSTANT_NORM, + "P1_CN": InterpolatorType.PIECEWISE_LINEAR_CONSTANT_NORM, + +} + +# Define the mapping after all imports interpolator_map = { InterpolatorType.BASE: GeologicalInterpolator, InterpolatorType.BASE_DISCRETE: DiscreteInterpolator, @@ -102,6 +94,8 @@ def __init__(self, *args, **kwargs): InterpolatorType.PIECEWISE_QUADRATIC: P2Interpolator, InterpolatorType.BASE_DATA_SUPPORTED: GeologicalInterpolator, InterpolatorType.SURFE: SurfeRBFInterpolator, + InterpolatorType.PIECEWISE_LINEAR_CONSTANT_NORM: ConstantNormP1Interpolator, + InterpolatorType.FINITE_DIFFERENCE_CONSTANT_NORM: ConstantNormFDIInterpolator, } support_interpolator_map = { @@ -119,9 +113,18 @@ def __init__(self, *args, **kwargs): 3: SupportType.DataSupported, 2: SupportType.DataSupported, }, + InterpolatorType.PIECEWISE_LINEAR_CONSTANT_NORM:{ + 3: SupportType.TetMesh, + 2: SupportType.P1Unstructured2d, + }, + InterpolatorType.FINITE_DIFFERENCE_CONSTANT_NORM: { + 3: SupportType.StructuredGrid, + 2: SupportType.StructuredGrid2D, + } } from ._interpolator_factory import InterpolatorFactory from ._interpolator_builder import InterpolatorBuilder -# from ._api import LoopInterpolator + + diff --git a/LoopStructural/interpolators/_constant_norm.py b/LoopStructural/interpolators/_constant_norm.py new file mode 100644 index 00000000..456c8a76 --- /dev/null +++ b/LoopStructural/interpolators/_constant_norm.py @@ -0,0 +1,205 @@ +import numpy as np + +from LoopStructural.interpolators._discrete_interpolator import DiscreteInterpolator +from LoopStructural.interpolators._finite_difference_interpolator import FiniteDifferenceInterpolator +from ._p1interpolator import P1Interpolator +from typing import Optional, Union, Callable +from scipy import sparse +from LoopStructural.utils import rng + +class ConstantNormInterpolator: + """Adds a non linear constraint to an interpolator to constrain + the norm of the gradient to be a set value. + + Returns + ------- + _type_ + _description_ + """ + def __init__(self, interpolator: DiscreteInterpolator,basetype): + """Initialise the constant norm inteprolator + with a discrete interpolator. + + Parameters + ---------- + interpolator : DiscreteInterpolator + The discrete interpolator to add constant norm to. + """ + self.basetype = basetype + self.interpolator = interpolator + self.support = interpolator.support + self.random_subset = False + self.norm_length = 1.0 + self.n_iterations = 20 + self.store_solution_history = False + self.solution_history = []#np.zeros((self.n_iterations, self.support.n_nodes)) + self.gradient_constraint_store = [] + def add_constant_norm(self, w:float): + """Add a constraint to the interpolator to constrain the norm of the gradient + to be a set value + + Parameters + ---------- + w : float + weighting of the constraint + """ + if "constant norm" in self.interpolator.constraints: + _ = self.interpolator.constraints.pop("constant norm") + + element_indices = np.arange(self.support.elements.shape[0]) + if self.random_subset: + rng.shuffle(element_indices) + element_indices = element_indices[: int(0.1 * self.support.elements.shape[0])] + vertices, gradient, elements, inside = self.support.get_element_gradient_for_location( + self.support.barycentre[element_indices] + ) + + t_g = gradient[:, :, :] + # t_n = gradient[self.support.shared_element_relationships[:, 1], :, :] + v_t = np.einsum( + "ijk,ik->ij", + t_g, + self.interpolator.c[self.support.elements[elements]], + ) + + v_t = v_t / np.linalg.norm(v_t, axis=1)[:, np.newaxis] + self.gradient_constraint_store.append(np.hstack([self.support.barycentre[element_indices],v_t])) + A1 = np.einsum("ij,ijk->ik", v_t, t_g) + volume = self.support.element_size[element_indices] + A1 = A1 / volume[:, np.newaxis] # normalise by element size + + b = np.zeros(A1.shape[0]) + self.norm_length + b = b / volume # normalise by element size + idc = np.hstack( + [ + self.support.elements[elements], + ] + ) + self.interpolator.add_constraints_to_least_squares(A1, b, idc, w=w, name="constant norm") + + def solve_system( + self, + solver: Optional[Union[Callable[[sparse.csr_matrix, np.ndarray], np.ndarray], str]] = None, + tol: Optional[float] = None, + solver_kwargs: dict = {}, + ) -> bool: + """Solve the system of equations iteratively for the constant norm interpolator. + + Parameters + ---------- + solver : Optional[Union[Callable[[sparse.csr_matrix, np.ndarray], np.ndarray], str]], optional + Solver function or name, by default None + tol : Optional[float], optional + Tolerance for the solver, by default None + solver_kwargs : dict, optional + Additional arguments for the solver, by default {} + + Returns + ------- + bool + Success status of the solver + """ + success = True + for i in range(self.n_iterations): + if i > 0: + self.add_constant_norm(w=(0.1 * i) ** 2 + 0.01) + # Ensure the interpolator is cast to P1Interpolator before calling solve_system + if isinstance(self.interpolator, self.basetype): + success = self.basetype.solve_system(self.interpolator, solver=solver, tol=tol, solver_kwargs=solver_kwargs) + if self.store_solution_history: + + self.solution_history.append(self.interpolator.c) + else: + raise TypeError("self.interpolator is not an instance of P1Interpolator") + if not success: + break + return success + +class ConstantNormP1Interpolator(P1Interpolator, ConstantNormInterpolator): + """Constant norm interpolator using P1 base interpolator + + Parameters + ---------- + P1Interpolator : class + The P1Interpolator class. + ConstantNormInterpolator : class + The ConstantNormInterpolator class. + """ + def __init__(self, support): + """Initialise the constant norm P1 interpolator. + + Parameters + ---------- + support : _type_ + _description_ + """ + P1Interpolator.__init__(self, support) + ConstantNormInterpolator.__init__(self, self, P1Interpolator) + + def solve_system( + self, + solver: Optional[Union[Callable[[sparse.csr_matrix, np.ndarray], np.ndarray], str]] = None, + tol: Optional[float] = None, + solver_kwargs: dict = {}, + ) -> bool: + """Solve the system of equations for the constant norm P1 interpolator. + + Parameters + ---------- + solver : Optional[Union[Callable[[sparse.csr_matrix, np.ndarray], np.ndarray], str]], optional + Solver function or name, by default None + tol : Optional[float], optional + Tolerance for the solver, by default None + solver_kwargs : dict, optional + Additional arguments for the solver, by default {} + + Returns + ------- + bool + Success status of the solver + """ + return ConstantNormInterpolator.solve_system(self, solver=solver, tol=tol, solver_kwargs=solver_kwargs) + +class ConstantNormFDIInterpolator(FiniteDifferenceInterpolator, ConstantNormInterpolator): + """Constant norm interpolator using finite difference base interpolator + + Parameters + ---------- + FiniteDifferenceInterpolator : class + The FiniteDifferenceInterpolator class. + ConstantNormInterpolator : class + The ConstantNormInterpolator class. + """ + def __init__(self, support): + """Initialise the constant norm finite difference interpolator. + + Parameters + ---------- + support : _type_ + _description_ + """ + FiniteDifferenceInterpolator.__init__(self, support) + ConstantNormInterpolator.__init__(self, self, FiniteDifferenceInterpolator) + def solve_system( + self, + solver: Optional[Union[Callable[[sparse.csr_matrix, np.ndarray], np.ndarray], str]] = None, + tol: Optional[float] = None, + solver_kwargs: dict = {}, + ) -> bool: + """Solve the system of equations for the constant norm finite difference interpolator. + + Parameters + ---------- + solver : Optional[Union[Callable[[sparse.csr_matrix, np.ndarray], np.ndarray], str]], optional + Solver function or name, by default None + tol : Optional[float], optional + Tolerance for the solver, by default None + solver_kwargs : dict, optional + Additional arguments for the solver, by default {} + + Returns + ------- + bool + Success status of the solver + """ + return ConstantNormInterpolator.solve_system(self, solver=solver, tol=tol, solver_kwargs=solver_kwargs) \ No newline at end of file diff --git a/LoopStructural/interpolators/_discrete_interpolator.py b/LoopStructural/interpolators/_discrete_interpolator.py index 60f19932..f861ce47 100644 --- a/LoopStructural/interpolators/_discrete_interpolator.py +++ b/LoopStructural/interpolators/_discrete_interpolator.py @@ -219,6 +219,7 @@ def add_constraints_to_least_squares(self, A, B, idc, w=1.0, name="undefined"): # normalise by rows of A # Should this be done? It should make the solution more stable length = np.linalg.norm(A, axis=1) + # length[length>0] = 1. B[length > 0] /= length[length > 0] # going to assume if any are nan they are all nan mask = np.any(np.isnan(A), axis=1) diff --git a/LoopStructural/interpolators/_interpolatortype.py b/LoopStructural/interpolators/_interpolatortype.py new file mode 100644 index 00000000..a68e0c91 --- /dev/null +++ b/LoopStructural/interpolators/_interpolatortype.py @@ -0,0 +1,22 @@ +from enum import Enum + +class InterpolatorType(Enum): + """ + Enum for the different interpolator types + + Each value is a unique identifier. + """ + + BASE = "BASE" + BASE_DISCRETE = "BASE_DISCRETE" + FINITE_DIFFERENCE = "FINITE_DIFFERENCE" + DISCRETE_FOLD = "DISCRETE_FOLD" + PIECEWISE_LINEAR = "PIECEWISE_LINEAR" + PIECEWISE_QUADRATIC = "PIECEWISE_QUADRATIC" + BASE_DATA_SUPPORTED = "BASE_DATA_SUPPORTED" + SURFE = "SURFE" + PIECEWISE_LINEAR_CONSTANT_NORM = "PIECEWISE_LINEAR_CONSTANT_NORM" + FINITE_DIFFERENCE_CONSTANT_NORM = "FINITE_DIFFERENCE_CONSTANT_NORM" + + + diff --git a/LoopStructural/interpolators/_p1interpolator.py b/LoopStructural/interpolators/_p1interpolator.py index 90da04a2..2d9b7099 100644 --- a/LoopStructural/interpolators/_p1interpolator.py +++ b/LoopStructural/interpolators/_p1interpolator.py @@ -45,7 +45,7 @@ def add_norm_constraints(self, w=1.0): points = self.get_norm_constraints() if points.shape[0] > 0: grad, elements, inside = self.support.evaluate_shape_derivatives(points[:, :3]) - size = self.support.element_size[elements[inside]] + size = self.support.element_scale[elements[inside]] wt = np.ones(size.shape[0]) wt *= w # s* size elements = np.tile(self.support.elements[elements[inside]], (3, 1, 1)) @@ -70,6 +70,7 @@ def add_value_constraints(self, w=1.0): if points.shape[0] > 1: N, elements, inside = self.support.evaluate_shape(points[:, :3]) size = self.support.element_size[elements[inside]] + wt = np.ones(size.shape[0]) wt *= w # * size self.add_constraints_to_least_squares( @@ -110,13 +111,16 @@ def minimise_edge_jumps(self, w=0.1, vector_func=None, vector=None, name="edge j const_n = -np.einsum("ij,ijk->ik", norm, Dn) # const_t_cp2 = np.einsum('ij,ikj->ik',normal,cp2_Dt) # const_n_cp2 = -np.einsum('ij,ikj->ik',normal,cp2_Dn) - + # shared_element_size = self.support.shared_element_size + # const_t /= shared_element_size[:, None] # normalise by element size + # const_n /= shared_element_size[:, None] # normalise by element size const = np.hstack([const_t, const_n]) # get vertex indexes tri_cp1 = np.hstack([self.support.elements[tri1], self.support.elements[tri2]]) # tri_cp2 = np.hstack([self.support.elements[cp2_tri1],self.support.elements[tri2]]) # add cp1 and cp2 to the least squares system + self.add_constraints_to_least_squares( const, np.zeros(const.shape[0]), diff --git a/LoopStructural/interpolators/supports/_3d_structured_tetra.py b/LoopStructural/interpolators/supports/_3d_structured_tetra.py index 70ecba37..baf00b91 100644 --- a/LoopStructural/interpolators/supports/_3d_structured_tetra.py +++ b/LoopStructural/interpolators/supports/_3d_structured_tetra.py @@ -85,10 +85,14 @@ def element_size(self): ) return np.abs(np.linalg.det(vecs)) / 6 - + @property def element_scale(self): - return self.element_size / np.mean(self.element_size) + size = self.element_size + size-= np.min(size) + size/= np.max(size) + size+=1. + return size @property def barycentre(self) -> np.ndarray: @@ -189,7 +193,7 @@ def shared_element_size(self): """ norm = self.shared_element_norm return 0.5 * np.linalg.norm(norm, axis=1) - + @property def shared_element_scale(self): return self.shared_element_size / np.mean(self.shared_element_size) diff --git a/LoopStructural/modelling/features/_base_geological_feature.py b/LoopStructural/modelling/features/_base_geological_feature.py index e49291a1..56171c53 100644 --- a/LoopStructural/modelling/features/_base_geological_feature.py +++ b/LoopStructural/modelling/features/_base_geological_feature.py @@ -347,6 +347,42 @@ def scalar_field(self, bounding_box=None): grid.cell_properties[self.name] = value return grid + def gradient_norm_scalar_field(self, bounding_box=None): + """Create a scalar field for the gradient norm of the feature + + Parameters + ---------- + bounding_box : Optional[BoundingBox], optional + bounding box to evaluate the scalar field in, by default None + + Returns + ------- + np.ndarray + scalar field of the gradient norm + """ + if bounding_box is None: + if self.model is None: + raise ValueError("Must specify bounding box") + bounding_box = self.model.bounding_box + grid = bounding_box.structured_grid(name=self.name) + value = np.linalg.norm( + self.evaluate_gradient(bounding_box.regular_grid(local=False, order='F')), + axis=1, + ) + if self.model is not None: + value = np.linalg.norm( + self.evaluate_gradient( + self.model.scale(bounding_box.regular_grid(local=False, order='F')) + ), + axis=1, + ) + grid.properties[self.name] = value + + value = np.linalg.norm( + self.evaluate_gradient(bounding_box.cell_centres(order='F')), axis=1 + ) + grid.cell_properties[self.name] = value + return grid def vector_field(self, bounding_box=None, tolerance=0.05, scale=1.0): """Create a vector field for the feature