From 7e6da5b417702d29ccee8435b0ce734d62046b0a Mon Sep 17 00:00:00 2001 From: Lachlan Grose Date: Mon, 2 Sep 2024 18:39:44 +1000 Subject: [PATCH] fix: gradient of faulted feature will not be aliased by the interpolation grid --- .../modelling/features/_geological_feature.py | 56 +++++++++++++++---- .../features/fault/_fault_segment.py | 52 +++-------------- LoopStructural/utils/maths.py | 51 +++++++++++++++++ 3 files changed, 104 insertions(+), 55 deletions(-) diff --git a/LoopStructural/modelling/features/_geological_feature.py b/LoopStructural/modelling/features/_geological_feature.py index bf1b3d2a..557134d3 100644 --- a/LoopStructural/modelling/features/_geological_feature.py +++ b/LoopStructural/modelling/features/_geological_feature.py @@ -2,6 +2,7 @@ Geological features """ +from LoopStructural.utils.maths import regular_tetraherdron_for_points, gradient_from_tetrahedron from ...modelling.features import BaseFeature from ...utils import getLogger from ...modelling.features import FeatureType @@ -131,7 +132,9 @@ def evaluate_value(self, pos: np.ndarray, ignore_regions=False, fillnan=None) -> v[nanmask] = v[~nanmask][i] return v - def evaluate_gradient(self, pos: np.ndarray, ignore_regions=False) -> np.ndarray: + def evaluate_gradient( + self, pos: np.ndarray, ignore_regions=False, element_scale_parameter=None + ) -> np.ndarray: """ Parameters @@ -147,6 +150,17 @@ def evaluate_gradient(self, pos: np.ndarray, ignore_regions=False) -> np.ndarray if pos.shape[1] != 3: raise LoopValueError("Need Nx3 array of xyz points to evaluate gradient") logger.info(f'Calculating gradient for {self.name}') + if element_scale_parameter is None: + if self.model is not None: + element_scale_parameter = np.min(self.model.bounding_box.step_vector) / 10 + else: + element_scale_parameter = 1 + else: + try: + element_scale_parameter = float(element_scale_parameter) + except ValueError: + logger.error("element_scale_parameter must be a float") + element_scale_parameter = 1 self.builder.up_to_date() @@ -156,16 +170,38 @@ def evaluate_gradient(self, pos: np.ndarray, ignore_regions=False) -> np.ndarray # evaluate the faults on the nodes of the faulted feature support # then evaluate the gradient at these points if len(self.faults) > 0: + # generate a regular tetrahedron for each point + # we will then move these points by the fault and then recalculate the gradient. + # this should work... + resolved = False + tetrahedron = regular_tetraherdron_for_points(pos, element_scale_parameter) + + while resolved: + for f in self.faults: + v = ( + f[0] + .evaluate_value(tetrahedron.reshape(-1, 3), fillnan='nearest') + .reshape(tetrahedron.shape[0], 4) + ) + flag = np.logical_or(np.all(v > 0, axis=1), np.all(v < 0, axis=1)) + if np.any(~flag): + logger.warning( + f"Points are too close to fault {f[0].name}. Refining the tetrahedron" + ) + element_scale_parameter *= 0.5 + tetrahedron = regular_tetraherdron_for_points(pos, element_scale_parameter) + + resolved = True + + tetrahedron_faulted = self._apply_faults(np.array(tetrahedron.reshape(-1, 3))).reshape( + tetrahedron.shape + ) + + values = self.interpolator.evaluate_value(tetrahedron_faulted.reshape(-1, 3)).reshape( + (-1, 4) + ) + v[mask, :] = gradient_from_tetrahedron(tetrahedron_faulted[mask, :, :], values[mask]) - if issubclass(type(self.interpolator), DiscreteInterpolator): - points = self.interpolator.support.nodes - else: - raise NotImplementedError( - "Faulted feature gradients are only supported by DiscreteInterpolator at the moment." - ) - points_faulted = self._apply_faults(points) - values = self.interpolator.evaluate_value(points_faulted) - v[mask, :] = self.interpolator.support.evaluate_gradient(pos[mask, :], values) return v pos = self._apply_faults(pos) if mask.dtype not in [int, bool]: diff --git a/LoopStructural/modelling/features/fault/_fault_segment.py b/LoopStructural/modelling/features/fault/_fault_segment.py index eb770bcb..544a58ae 100644 --- a/LoopStructural/modelling/features/fault/_fault_segment.py +++ b/LoopStructural/modelling/features/fault/_fault_segment.py @@ -1,3 +1,4 @@ +from LoopStructural.utils.maths import regular_tetraherdron_for_points, gradient_from_tetrahedron from ....modelling.features.fault._fault_function_feature import ( FaultDisplacementFeature, ) @@ -390,51 +391,12 @@ def apply_to_vectors(self, vector: np.ndarray, scale_parameter: float = 1.0) -> # the nodes of the tetrahedron are then restored by the fault and then gradient is # recalculated using the updated node positions but the original corner values - regular_tetrahedron = np.array( - [ - [np.sqrt(8 / 9), 0, -1 / 3], - [-np.sqrt(2 / 9), np.sqrt(2 / 3), -1 / 3], - [-np.sqrt(2 / 9), -np.sqrt(2 / 3), -1 / 3], - [0, 0, 1], - ] - ) - regular_tetrahedron *= scale_parameter - xyz = vector[:, :3] - tetrahedron = np.zeros((xyz.shape[0], 4, 3)) - tetrahedron[:] = xyz[:, None, :] - tetrahedron[:, :, :] += regular_tetrahedron[None, :, :] - - vectors = vector[:, 3:] - corners = np.einsum('ikj,ij->ik', tetrahedron - xyz[:, None, :], vectors) - tetrahedron = tetrahedron.reshape(-1, 3) - tetrahedron = self.apply_to_points(tetrahedron) - tetrahedron = tetrahedron.reshape(-1, 4, 3) - m = np.array( - [ - [ - (tetrahedron[:, 1, 0] - tetrahedron[:, 0, 0]), - (tetrahedron[:, 1, 1] - tetrahedron[:, 0, 1]), - (tetrahedron[:, 1, 2] - tetrahedron[:, 0, 2]), - ], - [ - (tetrahedron[:, 2, 0] - tetrahedron[:, 0, 0]), - (tetrahedron[:, 2, 1] - tetrahedron[:, 0, 1]), - (tetrahedron[:, 2, 2] - tetrahedron[:, 0, 2]), - ], - [ - (tetrahedron[:, 3, 0] - tetrahedron[:, 0, 0]), - (tetrahedron[:, 3, 1] - tetrahedron[:, 0, 1]), - (tetrahedron[:, 3, 2] - tetrahedron[:, 0, 2]), - ], - ] - ) - I = np.array([[-1.0, 1.0, 0.0, 0.0], [-1.0, 0.0, 1.0, 0.0], [-1.0, 0.0, 0.0, 1.0]]) - m = np.swapaxes(m, 0, 2) - element_gradients = np.linalg.inv(m) - - element_gradients = element_gradients.swapaxes(1, 2) - element_gradients = element_gradients @ I - v = np.sum(element_gradients * corners[:, None, :], axis=2) + tetrahedron = regular_tetraherdron_for_points(vector[:, :3], scale_parameter) + corners = np.einsum('ikj,ij->ik', tetrahedron - vector[:, None, :3], vector[:, 3:]) + + tetrahedron = self.apply_to_points(tetrahedron.reshape(-1, 3)).reshape(-1, 4, 3) + v = gradient_from_tetrahedron(tetrahedron, corners) + return v def add_abutting_fault(self, abutting_fault_feature, positive=None): diff --git a/LoopStructural/utils/maths.py b/LoopStructural/utils/maths.py index 1d4d8f98..caf8247d 100644 --- a/LoopStructural/utils/maths.py +++ b/LoopStructural/utils/maths.py @@ -244,3 +244,54 @@ def get_dip_vector(strike, dip): ] ) return v + + +def regular_tetraherdron_for_points(xyz, scale_parameter): + regular_tetrahedron = np.array( + [ + [np.sqrt(8 / 9), 0, -1 / 3], + [-np.sqrt(2 / 9), np.sqrt(2 / 3), -1 / 3], + [-np.sqrt(2 / 9), -np.sqrt(2 / 3), -1 / 3], + [0, 0, 1], + ] + ) + regular_tetrahedron *= scale_parameter + tetrahedron = np.zeros((xyz.shape[0], 4, 3)) + tetrahedron[:] = xyz[:, None, :] + tetrahedron[:, :, :] += regular_tetrahedron[None, :, :] + + return tetrahedron + + +def gradient_from_tetrahedron(tetrahedron, value): + """ + Calculate the gradient from a tetrahedron + """ + tetrahedron = tetrahedron.reshape(-1, 4, 3) + m = np.array( + [ + [ + (tetrahedron[:, 1, 0] - tetrahedron[:, 0, 0]), + (tetrahedron[:, 1, 1] - tetrahedron[:, 0, 1]), + (tetrahedron[:, 1, 2] - tetrahedron[:, 0, 2]), + ], + [ + (tetrahedron[:, 2, 0] - tetrahedron[:, 0, 0]), + (tetrahedron[:, 2, 1] - tetrahedron[:, 0, 1]), + (tetrahedron[:, 2, 2] - tetrahedron[:, 0, 2]), + ], + [ + (tetrahedron[:, 3, 0] - tetrahedron[:, 0, 0]), + (tetrahedron[:, 3, 1] - tetrahedron[:, 0, 1]), + (tetrahedron[:, 3, 2] - tetrahedron[:, 0, 2]), + ], + ] + ) + I = np.array([[-1.0, 1.0, 0.0, 0.0], [-1.0, 0.0, 1.0, 0.0], [-1.0, 0.0, 0.0, 1.0]]) + m = np.swapaxes(m, 0, 2) + element_gradients = np.linalg.inv(m) + + element_gradients = element_gradients.swapaxes(1, 2) + element_gradients = element_gradients @ I + v = np.sum(element_gradients * value[:, None, :], axis=2) + return v