From 0f2983bb3ccc2fc57dc3bb5a712f83276914118e Mon Sep 17 00:00:00 2001 From: Lachlan Grose Date: Mon, 30 Jun 2025 13:17:45 +1000 Subject: [PATCH 1/8] fix: adding constant norm interpolators --- LoopStructural/interpolators/__init__.py | 64 +++--- .../interpolators/_constant_norm.py | 196 ++++++++++++++++++ .../interpolators/_interpolatortype.py | 22 ++ 3 files changed, 252 insertions(+), 30 deletions(-) create mode 100644 LoopStructural/interpolators/_constant_norm.py create mode 100644 LoopStructural/interpolators/_interpolatortype.py diff --git a/LoopStructural/interpolators/__init__.py b/LoopStructural/interpolators/__init__.py index 56e68480..29aa3d9c 100644 --- a/LoopStructural/interpolators/__init__.py +++ b/LoopStructural/interpolators/__init__.py @@ -3,6 +3,7 @@ """ + __all__ = [ "InterpolatorType", "GeologicalInterpolator", @@ -21,39 +22,13 @@ "StructuredGrid2D", "P2UnstructuredTetMesh", ] -from enum import IntEnum +from ._interpolatortype import InterpolatorType, add_interpolator_type + 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 +54,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 +68,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 +95,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 +114,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..9d262a9c --- /dev/null +++ b/LoopStructural/interpolators/_constant_norm.py @@ -0,0 +1,196 @@ +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): + """Initialise the constant norm inteprolator + with a discrete interpolator. + + Parameters + ---------- + interpolator : DiscreteInterpolator + The discrete interpolator to add constant norm to. + """ + self.interpolator = interpolator + self.support = interpolator.support + self.random_subset = False + self.norm_length = 1.0 + 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] + A1 = np.einsum("ij,ijk->ik", v_t, t_g) + + b = np.zeros(A1.shape[0]) + self.norm_length + 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: + """ + + Parameters + ---------- + solver : Optional[Union[Callable[[sparse.csr_matrix, np.ndarray], np.ndarray], str]], optional + _description_, by default None + tol : Optional[float], optional + _description_, by default None + solver_kwargs : dict, optional + _description_, by default {} + + Returns + ------- + bool + _description_ + """ + for i in range(20): + if i > 0: + self.add_constant_norm(w=(0.1 * i) ** 2 + 0.01) + success = self.interpolator.solve_system(solver=solver, tol=tol, solver_kwargs=solver_kwargs) + return True + +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) + + 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 + _description_, by default None + tol : Optional[float], optional + _description_, by default None + solver_kwargs : dict, optional + _description_, by default {} + + Returns + ------- + bool + _description_ + """ + success = True + for i in range(20): + + if i > 0: + self.add_constant_norm(w=(0.1 * i) ** 2 + 0.01) + success = P1Interpolator.solve_system(self, solver, tol, solver_kwargs) + return success +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) + 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 + _description_, by default None + tol : Optional[float], optional + _description_, by default None + solver_kwargs : dict, optional + _description_, by default {} + + Returns + ------- + bool + _description_ + """ + success=True + for i in range(20): + if i > 0: + self.add_constant_norm(w=(0.1 * i) ** 2 + 0.01) + success = FiniteDifferenceInterpolator.solve_system(self, solver, tol, solver_kwargs) + return success \ No newline at end of file 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" + + + From ee51c680125d0064e8edd803108c4271a24f9c06 Mon Sep 17 00:00:00 2001 From: Lachlan Grose Date: Mon, 30 Jun 2025 13:18:05 +1000 Subject: [PATCH 2/8] fix: add visualisation to plot gradient norm --- .../features/_base_geological_feature.py | 36 +++++++++++++++++++ 1 file changed, 36 insertions(+) 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 From 3026e3cde780f2b222752c6d76312ad5d5508159 Mon Sep 17 00:00:00 2001 From: Lachlan Grose Date: Mon, 30 Jun 2025 13:45:22 +1000 Subject: [PATCH 3/8] fix: update vector scaling --- LoopStructural/datatypes/_point.py | 1 - 1 file changed, 1 deletion(-) 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') From b8162186c74a852eb7b5cd61808be7f41c989741 Mon Sep 17 00:00:00 2001 From: Lachlan Grose Date: Mon, 30 Jun 2025 13:45:34 +1000 Subject: [PATCH 4/8] remove add_interpolator_type --- LoopStructural/interpolators/__init__.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/LoopStructural/interpolators/__init__.py b/LoopStructural/interpolators/__init__.py index 29aa3d9c..401a00b8 100644 --- a/LoopStructural/interpolators/__init__.py +++ b/LoopStructural/interpolators/__init__.py @@ -22,8 +22,7 @@ "StructuredGrid2D", "P2UnstructuredTetMesh", ] -from ._interpolatortype import InterpolatorType, add_interpolator_type - +from ._interpolatortype import InterpolatorType from ..utils import getLogger From 4fafcafb7a1ea98d2da0a19b9047df80de931b0d Mon Sep 17 00:00:00 2001 From: Lachlan Grose Date: Mon, 30 Jun 2025 13:48:32 +1000 Subject: [PATCH 5/8] linting --- LoopStructural/interpolators/_constant_norm.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/LoopStructural/interpolators/_constant_norm.py b/LoopStructural/interpolators/_constant_norm.py index 9d262a9c..c7a2bfce 100644 --- a/LoopStructural/interpolators/_constant_norm.py +++ b/LoopStructural/interpolators/_constant_norm.py @@ -90,11 +90,12 @@ def solve_system( bool _description_ """ + success = True for i in range(20): if i > 0: self.add_constant_norm(w=(0.1 * i) ** 2 + 0.01) success = self.interpolator.solve_system(solver=solver, tol=tol, solver_kwargs=solver_kwargs) - return True + return success class ConstantNormP1Interpolator(P1Interpolator, ConstantNormInterpolator): """Constant norm interpolator using P1 base interpolator From 84870b6bf7b729908f9af82d9c2786bd240374fd Mon Sep 17 00:00:00 2001 From: Lachlan Grose Date: Mon, 30 Jun 2025 14:10:20 +1000 Subject: [PATCH 6/8] reduce code duplication --- .../interpolators/_constant_norm.py | 58 +++++++++---------- 1 file changed, 28 insertions(+), 30 deletions(-) diff --git a/LoopStructural/interpolators/_constant_norm.py b/LoopStructural/interpolators/_constant_norm.py index c7a2bfce..d7da45ec 100644 --- a/LoopStructural/interpolators/_constant_norm.py +++ b/LoopStructural/interpolators/_constant_norm.py @@ -16,7 +16,7 @@ class ConstantNormInterpolator: _type_ _description_ """ - def __init__(self, interpolator: DiscreteInterpolator): + def __init__(self, interpolator: DiscreteInterpolator,basetype): """Initialise the constant norm inteprolator with a discrete interpolator. @@ -25,10 +25,12 @@ def __init__(self, interpolator: DiscreteInterpolator): 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 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 @@ -74,27 +76,33 @@ def solve_system( 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 - _description_, by default None + Solver function or name, by default None tol : Optional[float], optional - _description_, by default None + Tolerance for the solver, by default None solver_kwargs : dict, optional - _description_, by default {} + Additional arguments for the solver, by default {} Returns ------- bool - _description_ + Success status of the solver """ success = True - for i in range(20): + for i in range(self.n_iterations): if i > 0: self.add_constant_norm(w=(0.1 * i) ** 2 + 0.01) - success = self.interpolator.solve_system(solver=solver, tol=tol, solver_kwargs=solver_kwargs) + # 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) + else: + raise TypeError("self.interpolator is not an instance of P1Interpolator") + if not success: + break return success class ConstantNormP1Interpolator(P1Interpolator, ConstantNormInterpolator): @@ -116,7 +124,7 @@ def __init__(self, support): _description_ """ P1Interpolator.__init__(self, support) - ConstantNormInterpolator.__init__(self, self) + ConstantNormInterpolator.__init__(self, self, P1Interpolator) def solve_system( self, @@ -129,24 +137,19 @@ def solve_system( Parameters ---------- solver : Optional[Union[Callable[[sparse.csr_matrix, np.ndarray], np.ndarray], str]], optional - _description_, by default None + Solver function or name, by default None tol : Optional[float], optional - _description_, by default None + Tolerance for the solver, by default None solver_kwargs : dict, optional - _description_, by default {} + Additional arguments for the solver, by default {} Returns ------- bool - _description_ + Success status of the solver """ - success = True - for i in range(20): + return ConstantNormInterpolator.solve_system(self, solver=solver, tol=tol, solver_kwargs=solver_kwargs) - if i > 0: - self.add_constant_norm(w=(0.1 * i) ** 2 + 0.01) - success = P1Interpolator.solve_system(self, solver, tol, solver_kwargs) - return success class ConstantNormFDIInterpolator(FiniteDifferenceInterpolator, ConstantNormInterpolator): """Constant norm interpolator using finite difference base interpolator @@ -166,7 +169,7 @@ def __init__(self, support): _description_ """ FiniteDifferenceInterpolator.__init__(self, support) - ConstantNormInterpolator.__init__(self, self) + ConstantNormInterpolator.__init__(self, self, FiniteDifferenceInterpolator) def solve_system( self, solver: Optional[Union[Callable[[sparse.csr_matrix, np.ndarray], np.ndarray], str]] = None, @@ -178,20 +181,15 @@ def solve_system( Parameters ---------- solver : Optional[Union[Callable[[sparse.csr_matrix, np.ndarray], np.ndarray], str]], optional - _description_, by default None + Solver function or name, by default None tol : Optional[float], optional - _description_, by default None + Tolerance for the solver, by default None solver_kwargs : dict, optional - _description_, by default {} + Additional arguments for the solver, by default {} Returns ------- bool - _description_ + Success status of the solver """ - success=True - for i in range(20): - if i > 0: - self.add_constant_norm(w=(0.1 * i) ** 2 + 0.01) - success = FiniteDifferenceInterpolator.solve_system(self, solver, tol, solver_kwargs) - return success \ No newline at end of file + return ConstantNormInterpolator.solve_system(self, solver=solver, tol=tol, solver_kwargs=solver_kwargs) \ No newline at end of file From 39d273b9917cc48d01ecb026b172a75bd817a173 Mon Sep 17 00:00:00 2001 From: Lachlan Grose Date: Fri, 4 Jul 2025 12:26:13 +1000 Subject: [PATCH 7/8] adding option to store solutions for each iteration and constraints --- LoopStructural/interpolators/_constant_norm.py | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/LoopStructural/interpolators/_constant_norm.py b/LoopStructural/interpolators/_constant_norm.py index d7da45ec..456c8a76 100644 --- a/LoopStructural/interpolators/_constant_norm.py +++ b/LoopStructural/interpolators/_constant_norm.py @@ -31,6 +31,9 @@ def __init__(self, interpolator: DiscreteInterpolator,basetype): 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 @@ -60,9 +63,13 @@ def add_constant_norm(self, w:float): ) 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], @@ -99,6 +106,9 @@ def solve_system( # 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: From b45c56a1096db991633c84894c616d3873e4fc4f Mon Sep 17 00:00:00 2001 From: Lachlan Grose Date: Fri, 4 Jul 2025 12:26:33 +1000 Subject: [PATCH 8/8] scale by element size --- LoopStructural/interpolators/_discrete_interpolator.py | 1 + LoopStructural/interpolators/_p1interpolator.py | 8 ++++++-- .../interpolators/supports/_3d_structured_tetra.py | 10 +++++++--- 3 files changed, 14 insertions(+), 5 deletions(-) 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/_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)