From 57011003fa30fdc7738976a7b829c14d42bd25ff Mon Sep 17 00:00:00 2001 From: Lachlan Grose Date: Mon, 4 Mar 2024 13:50:09 +1100 Subject: [PATCH] 2d interpolation (#176) * fix: adding json exports * fix: moving and renaming conversion functions also adding some type hints and a generic numeric type * tests: adding tests for conversions need to add some better tests not just 0 and 1s * fix: remove process_map2loop not required, can be done using processor * fix: adding synthetic horizontal layers model * fix: adding corners property to bounding box * fix: substituting pli for p1 same interpolator but this implementation uses numpy instead of cython * tests: test that horizontal layers are working * fix: making structured tetra compatible with p1 interp * fix: making p1 interpolator work * fix: cast input to np array * fix: adding type enums * fix: adding sparse matrix neighbour calculator this means no more cython code * fix: adding some type hints * refactor: :heavy_minus_sign: removing model builder there is no real advantage to using a builder design pattern for loopstructural * fix: :sparkles: adding ability to calculate normalised value of a feature this just makes it easier to scale scalar field between 0 and 1 * fix: adding map for interpolator to support type * fix: enabling support builder with interpolator factory * fix: updating imports for maths functions * fix: adding support creator from bbox * fix: when reset clear the constraints dictionary * fix: reset when setting up interpolator * fix: removing old pli interpolator * fix: fixing tetra indexing and masks * fix: removing unused imports * fix: speeding up element getter * fix: adding onGeometryChange function for supports * fix: reset change number of constraints * fix: adding gradient orthogonal to p1interpolator * fix: changing instrusion example to rescale points * fix: setting model tolerance using new bounding box * fix: updating fold to use p1 interpolator * fix: renaming variables * fix: adding name and vector args to min edge jumps * fix: removing old imports * fix: fixing orthogonal constraints for p1 * fix: moving fault setup into fault builder frm model * fix: flake8 error * fix: add bb origin to surfaces * fix: fold constraints added in setup interpolator this allows for interpolator to be reset and doesn't remove the fold constraints when resetting * fix: temp fix for issue with np types * Replace reqs with correct theme Theme being loaded by requirements.txt was wrong (an old one?) So have updated with the pydata_sphinx_theme so now builds correctly * Adding PyPI custom link to navigation * Replacing square icon with (nicer) round icon * Adding navbar options and metadata Image is a placeholder for now, but this is an (onvs opinionated) view on a nicer format. Adds a limit to the number of nav items, adds a logo, changes the title to look cleaner (and not say "documentation" at the end), adds metadata * Fixing svg img * Updating the contributors docs * removing accidental tab * docs: fixing missing automodule * fix: adding length to bbox * fix: removing dof and name from interpolator json * docs: adding explicit imports for submodules Without explicit import of submodules and using __ALL__ the majority of the code isn't imported so sphinx doesn't document it. * docs: adding citations to documentation * fix: removing reference to dtm creator and using get interpolator correctly * fix: face table resize when geometry reset * fix: elements property for base structured support * fix: adding imports to api * fix: changing from interpolator type string to enum * fix: removing cython * fix: only axial foliation is constrained by fold frame * Update release-please.yml * Update release-please.yml * Update release-please.yml * ci: updating to use pip wheel . * Update release-please.yml * fix: disable progress bar for model update * ci: remove numpy version for conda * fix: changing to interpolator factory required changing how fold interpolator is setup. Rather than initialising with a fold, the fold is added to the interpolator added map between stringinterpolator names and the interpolatortypes enum * fix: identify points that fall on face of shared tetras. Automatically chooses the second one. * fix: rename surface vtk property to vtk from pyvista * fix: add docstring and cast isinside to an array * docs: adding docs to surface datatype * fix: removing cython code * fix: removing model api for time being * build: adding conda builds to git ignore * typo * ci: using conda convert to create osx * ci: removing cython * ci: typo in yml * ci: again * ci: again * fix: removing analysis module, will create separate ptyhon library * style: style fixes by ruff and autoformatting by black * fix: changing import of bounding box * fix: migrating np.random.shuffle to use generator object * fix: removing axis/angle from apply faults * fix: calculate gradient of faulted feature returns true gradient * fix: bb import update and some formatting * style: style fixes by ruff and autoformatting by black * fix: adding a method to rotate vectors that are faulted * fix: store constraints as matrix rather than dict the matrices can then be stacked together. no real advantage to the previous approach. It is useful if you want to use the matrix for a constraint outside of the typical LS workflow. * fix: removing unused variable * fix: adding --user flag and removing cython from test suite * fix: changing discrete support to return element index previously there was a mix of returning the index or the cells in the index. making supports consistent * fix: global indexes for elements * fix: linting * style: style fixes by ruff and autoformatting by black * fix: changing bounding box import * ci: updating action version numbers * fix: moving face table and aabb to separate file * fix: adding abc for base support * fix: updating 2d grid to follow same pattern as 3d * fix: updating 2d tri mesh to same pattern as 3d tetra * fix: typo vaue to value * fix: adding abstract methods for interpolation functions * fix: changing add methods to abstract method names * fix: use abstract base class for supports * fix: don't swap axes for dn, not consistent with 3d * docs: adding type hints * fix: adding get element gradient for location * removing prints * fix: updating p2 2d support * fix: linting * docs: type hints * fix: adding implementation of required abstract methods * tests: updating 2d local coords to use array * fix: adding gradient orthogonal abstract method * fix: adding methods for setting inequality constraints * fix: adding element size abstract method * fix: removing unused operator code * fix: adding faults setter/getter * fix: making intrusion feature compatible with base feature * fix: linting * fix: count number of constraints before adding previously not counting constraints before tying to add * style: style fixes by ruff and autoformatting by black * tests: don't test base feature anymore * style: style fixes by ruff and autoformatting by black --------- Co-authored-by: Sam Co-authored-by: lachlangrose --- .github/workflows/release-please.yml | 21 +- LoopStructural/interpolators/_builders.py | 10 +- .../interpolators/_discrete_interpolator.py | 82 +---- .../_finite_difference_interpolator.py | 30 +- .../interpolators/_geological_interpolator.py | 36 ++- LoopStructural/interpolators/_operator.py | 11 - .../interpolators/_p1interpolator.py | 15 +- .../interpolators/_p2interpolator.py | 28 +- .../supports/_2d_base_unstructured.py | 252 ++++++++++++---- .../supports/_2d_p1_unstructured.py | 14 +- .../supports/_2d_p2_unstructured.py | 88 +++--- .../supports/_2d_structured_grid.py | 282 +++++++++++------- .../supports/_3d_base_structured.py | 12 +- .../supports/_3d_unstructured_tetra.py | 34 ++- .../interpolators/supports/_aabb.py | 77 +++++ .../interpolators/supports/_base_support.py | 114 +++++++ .../interpolators/supports/_face_table.py | 70 +++++ .../features/_base_geological_feature.py | 31 +- .../builders/_geological_feature_builder.py | 17 +- .../modelling/intrusions/intrusion_feature.py | 16 +- tests/integration/test_interpolator.py | 12 + .../interpolator/test_2d_discrete_support.py | 6 +- .../unit/modelling/test_geological_feature.py | 10 +- 23 files changed, 896 insertions(+), 372 deletions(-) create mode 100644 LoopStructural/interpolators/supports/_aabb.py create mode 100644 LoopStructural/interpolators/supports/_base_support.py create mode 100644 LoopStructural/interpolators/supports/_face_table.py diff --git a/.github/workflows/release-please.yml b/.github/workflows/release-please.yml index ddff4af02..1332c827c 100644 --- a/.github/workflows/release-please.yml +++ b/.github/workflows/release-please.yml @@ -17,7 +17,7 @@ jobs: - name: Lint with ruff run: | ruff check LoopStructural --fix - - uses: stefanzweifel/git-auto-commit-action@v4 + - uses: stefanzweifel/git-auto-commit-action@v5 with: commit_message: "style: style fixes by ruff and autoformatting by black" @@ -31,7 +31,8 @@ jobs: python-version: ${{ fromJSON(vars.PYTHON_VERSIONS)}} steps: - uses: actions/checkout@v4 - - uses: conda-incubator/setup-miniconda@v2 + - uses: conda-incubator/setup-miniconda@v3 + with: python-version: ${{ matrix.python }} - name: Installing dependencies @@ -100,7 +101,7 @@ jobs: os: ["ubuntu-latest"] python-version: ${{ fromJSON(vars.PYTHON_VERSIONS)}} steps: - - uses: conda-incubator/setup-miniconda@v2 + - uses: conda-incubator/setup-miniconda@v3 with: auto-update-conda: true python-version: ${{ matrix.python-version }} @@ -120,7 +121,7 @@ jobs: conda convert -p all conda/linux-64/*.tar.bz2 -f -o conda conda install anaconda-client -y - name: upload artifacts - uses: actions/upload-artifact@v3 + uses: actions/upload-artifact@v4 with: name: conda-${{ matrix.os }}-${{ matrix.python-version }} path: conda @@ -138,11 +139,11 @@ jobs: python setup.py sdist pip wheel -w wheelhouse . - - uses: actions/upload-artifact@v3 + - uses: actions/upload-artifact@v4 with: name: dist path: dist/*.tar.gz - - uses: actions/upload-artifact@v3 + - uses: actions/upload-artifact@v4 with: name: dist path: wheelhouse/*.whl @@ -152,20 +153,20 @@ jobs: runs-on: ubuntu-latest if: ${{ needs.release-please.outputs.release_created }} steps: - - uses: actions/download-artifact@v3 + - uses: actions/download-artifact@v4 with: name: dist path: dist - - uses: actions/download-artifact@v3 + - uses: actions/download-artifact@v4 with: path: conda - - uses: pypa/gh-action-pypi-publish@v1.6.4 + - uses: pypa/gh-action-pypi-publish@v1 with: skip_existing: true verbose: true user: ${{ secrets.PYPI_USERNAME }} password: ${{ secrets.PYPI_PASSWORD }} - - uses: conda-incubator/setup-miniconda@v2 + - uses: conda-incubator/setup-miniconda@v3 - name: upload all files to conda-forge shell: bash -l {0} env: diff --git a/LoopStructural/interpolators/_builders.py b/LoopStructural/interpolators/_builders.py index bae46c703..414e04abd 100644 --- a/LoopStructural/interpolators/_builders.py +++ b/LoopStructural/interpolators/_builders.py @@ -10,7 +10,7 @@ StructuredGrid, TetMesh, ) -from LoopStructural.utils import BoundingBox +from LoopStructural.datatypes import BoundingBox from LoopStructural.utils.logging import getLogger logger = getLogger(__name__) @@ -61,6 +61,7 @@ def get_interpolator( ) return P1Interpolator(support) + return P1Interpolator(support) if interpolatortype == "P2": if support is not None: logger.info( @@ -69,9 +70,7 @@ def get_interpolator( ) return P2Interpolator(support) else: - raise ValueError( - "Cannot create P2 interpolator without support, try using PLI" - ) + raise ValueError("Cannot create P2 interpolator without support, try using PLI") if interpolatortype == "FDI": # find the volume of one element @@ -96,8 +95,7 @@ def get_interpolator( grid = StructuredGrid(origin=origin, nsteps=nsteps, step_vector=step_vector) logger.info( - f"Creating regular grid with {grid.n_elements} elements \n" - "for modelling using FDI" + f"Creating regular grid with {grid.n_elements} elements \n" "for modelling using FDI" ) return FiniteDifferenceInterpolator(grid) if interpolatortype == "DFI": diff --git a/LoopStructural/interpolators/_discrete_interpolator.py b/LoopStructural/interpolators/_discrete_interpolator.py index f68efd7b9..caff12cda 100644 --- a/LoopStructural/interpolators/_discrete_interpolator.py +++ b/LoopStructural/interpolators/_discrete_interpolator.py @@ -2,6 +2,7 @@ Discrete interpolator base for least squares """ +from abc import abstractmethod import logging from time import time @@ -39,16 +40,12 @@ def __init__(self, support, data={}, c=None, up_to_date=False): else np.zeros(self.support.n_nodes) ) self.region_function = lambda xyz: np.ones(xyz.shape[0], dtype=bool) - # self.region_map[self.region] = np.array(range(0, - # len(self.region_map[self.region]))) + self.shape = "rectangular" if self.shape == "square": self.B = np.zeros(self.nx) self.c_ = 0 - # self.A = [] # sparse matrix storage coo format - # self.col = [] - # self.row = [] # sparse matrix storage - # self.w = [] + self.solver = "cg" self.eq_const_C = [] @@ -221,6 +218,12 @@ def add_constraints_to_least_squares(self, A, B, idc, w=1.0, name="undefined"): 'w': w, } + @abstractmethod + def add_gradient_orthogonal_constraints( + self, points: np.ndarray, vectors: np.ndarray, w: float = 1.0 + ): + pass + def calculate_residual_for_constraints(self): """Calculates Ax-B for all constraints added to the interpolator This could be a proxy to identify which constraints are controlling the model @@ -252,7 +255,6 @@ def remove_constraints_from_least_squares(self, name="undefined", constraint_ids ------- """ - pass def add_equality_constraints(self, node_idx, values, name="undefined"): @@ -288,72 +290,6 @@ def add_equality_constraints(self, node_idx, values, name="undefined"): } self.eq_const_c += idc[outside].shape[0] - def add_inequality_constraints_to_matrix(self, A, l, u, idc, name="undefined"): - """Adds constraints for a matrix where the linear function - l < Ax > u constrains the objective function - - - Parameters - ---------- - A : numpy array - matrix of coefficients - l : numpy array - lower bounds - u : numpy array - upper bounds - idc : numpy array - index of constraints - Returns - ------- - - """ - # 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) - idc = gi[idc] - rows = np.arange(self.ineq_const_c, self.ineq_const_c + idc.shape[0]) - rows = np.tile(rows, (A.shape[-1], 1)).T - self.ineq_constraints[name] = {"A": A, "l": l, "col": idc, "u": u, "row": rows} - self.ineq_const_c += idc.shape[0] - - def add_inequality_feature(self, feature, lower=True, mask=None): - """Add an inequality constraint to the interpolator using an existing feature. - This will make the interpolator greater than or less than the exising feature. - Evaluate the feature at the interpolation nodes. - Can provide a boolean mask to restrict to only some parts - - Parameters - ---------- - feature : BaseFeature - the feature that will be used to constraint the interpolator - lower : bool, optional - lower or upper constraint, by default True - mask : np.ndarray, optional - restrict the nodes to evaluate on, by default None - """ - # add inequality value for the nodes of the mesh - # flag lower determines whether the feature is a lower bound or upper bound - # mask is just a boolean array determining which nodes to apply it to - - value = feature(self.support.nodes) - if mask is None: - mask = np.ones(value.shape[0], dtype=bool) - l = np.zeros(value.shape[0]) - np.inf - u = np.zeros(value.shape[0]) + np.inf - mask = np.logical_and(mask, ~np.isnan(value)) - if lower: - l[mask] = value[mask] - if not lower: - u[mask] = value[mask] - - self.add_inequality_constraints_to_matrix( - np.ones((value.shape[0], 1)), - l, - u, - np.arange(0, self.nx, dtype=int), - ) - def add_tangent_constraints(self, w=1.0): """Adds the constraints :math:`f(X)\cdotT=0` diff --git a/LoopStructural/interpolators/_finite_difference_interpolator.py b/LoopStructural/interpolators/_finite_difference_interpolator.py index 2f9919548..de733a802 100644 --- a/LoopStructural/interpolators/_finite_difference_interpolator.py +++ b/LoopStructural/interpolators/_finite_difference_interpolator.py @@ -121,7 +121,7 @@ def setup_interpolator(self, **kwargs): self.assemble_inner(operator, weight) self.add_norm_constraints(self.interpolation_weights["npw"]) self.add_gradient_constraints(self.interpolation_weights["gpw"]) - self.add_vaue_constraints(self.interpolation_weights["cpw"]) + self.add_value_constraints(self.interpolation_weights["cpw"]) self.add_tangent_constraints(self.interpolation_weights["tpw"]) self.add_interface_constraints(self.interpolation_weights["ipw"]) self.add_inequality_constraints() @@ -136,7 +136,7 @@ def copy(self): """ return FiniteDifferenceInterpolator(self.support) - def add_vaue_constraints(self, w=1.0): + def add_value_constraints(self, w=1.0): """ Parameters @@ -151,7 +151,9 @@ def add_vaue_constraints(self, w=1.0): points = self.get_value_constraints() # check that we have added some points if points.shape[0] > 0: - node_idx, inside = self.support.position_to_cell_corners(points[:, :3]) + node_idx, inside = self.support.position_to_cell_corners( + points[:, : self.support.dimension] + ) # print(points[inside,:].shape) gi = np.zeros(self.support.n_nodes, dtype=int) gi[:] = -1 @@ -160,14 +162,14 @@ def add_vaue_constraints(self, w=1.0): idc[:] = -1 idc[inside, :] = gi[node_idx[inside, :]] inside = np.logical_and(~np.any(idc == -1, axis=1), inside) - a = self.support.position_to_dof_coefs(points[inside, :3]) - # a*=w - # a/=np.product(self.support.step_vector) + a = self.support.position_to_dof_coefs(points[inside, : self.support.dimension]) + # a *= w + # a/=self.support.enp.product(self.support.step_vector) self.add_constraints_to_least_squares( a, points[inside, 3], idc[inside, :], - w=w * points[inside, 4], + w=w * points[inside, 4] * self.support.element_size, name="value", ) if np.sum(inside) <= 0: @@ -233,7 +235,9 @@ def add_interface_constraints(self, w=1.0): # get elements for points points = self.get_interface_constraints() if points.shape[0] > 1: - vertices, c, tetras, inside = self.support.get_element_for_location(points[:, :3]) + vertices, c, tetras, inside = self.support.get_element_for_location( + points[:, : self.support.dimension] + ) # calculate volume of tetras # vecs = vertices[inside, 1:, :] - vertices[inside, 0, None, :] # vol = np.abs(np.linalg.det(vecs)) / 6 @@ -391,7 +395,7 @@ def add_norm_constraints(self, w=1.0): self.up_to_date = False def add_gradient_orthogonal_constraints( - self, points: np.ndarray, vector: np.ndarray, w: float = 1.0, B: float = 0 + self, points: np.ndarray, vectors: np.ndarray, w: float = 1.0, B: float = 0 ): """ constraints scalar field to be orthogonal to a given vector @@ -424,8 +428,8 @@ def add_gradient_orthogonal_constraints( idc[inside, :] = gi[node_idx[inside, :]] inside = np.logical_and(~np.any(idc == -1, axis=1), inside) # normalise vector and scale element gradient matrix by norm as well - norm = np.linalg.norm(vector, axis=1) - vector[norm > 0, :] /= norm[norm > 0, None] + norm = np.linalg.norm(vectors, axis=1) + vectors[norm > 0, :] /= norm[norm > 0, None] # normalise element vector to unit vector for dot product ( @@ -437,7 +441,7 @@ def add_gradient_orthogonal_constraints( T[norm > 0, :, :] /= norm[norm > 0, None, None] # dot product of vector and element gradient = 0 - A = np.einsum("ij,ijk->ik", vector[inside, :3], T) + A = np.einsum("ij,ijk->ik", vectors[inside, :3], T) B = np.zeros(points[inside, :].shape[0]) + B self.add_constraints_to_least_squares( A, B, idc[inside, :], w=w * self.vol, name="gradient orthogonal" @@ -461,7 +465,7 @@ def add_regularisation(self, operator, w=0.1): ------- """ - self.assemble_inner(operator) + self.assemble_inner(operator, w) # self.assemble_borders() # def assemble_borders(self, operator, w): diff --git a/LoopStructural/interpolators/_geological_interpolator.py b/LoopStructural/interpolators/_geological_interpolator.py index c8ef841b4..e03180a1a 100644 --- a/LoopStructural/interpolators/_geological_interpolator.py +++ b/LoopStructural/interpolators/_geological_interpolator.py @@ -106,8 +106,8 @@ def set_value_constraints(self, points: np.ndarray): """ points = self.check_array(points) - if points.shape[1] < 4: - raise ValueError("Value points must at least have X,Y,Z,val") + if points.shape[1] < 5: + raise ValueError("Value points must at least have X,Y,Z,val,w") self.data["value"] = points self.n_i = points.shape[0] self.up_to_date = False @@ -172,10 +172,20 @@ def set_interface_constraints(self, points: np.ndarray): self.data["interface"] = points self.up_to_date = False - def set_inequality_constraints(self, points: np.ndarray): + def set_value_inequality_constraints(self, points: np.ndarray): + if points.shape[1] < 5: + raise ValueError("Inequality constraints must at least have X,Y,Z,lower,upper") self.data["inequality"] = points self.up_to_date = False + def set_inequality_pairs_constraints(self, pointsa: np.ndarray, pointsb: np.ndarray): + if pointsa.shape[1] < 3: + raise ValueError("Inequality pairs constraints must at least have X,Y,Z") + if pointsb.shape[1] < 3: + raise ValueError("Inequality pairs constraints must at least have X,Y,Z") + self.data["inequality_pairs"] = {"a": pointsa, "b": pointsb} + self.up_to_date = False + def get_value_constraints(self): """ @@ -274,6 +284,26 @@ def evaluate_gradient(self, locations: np.ndarray): def reset(self): pass + @abstractmethod + def add_value_constraints(self, w: float = 1.0): + pass + + @abstractmethod + def add_gradient_constraints(self, w: float = 1.0): + pass + + @abstractmethod + def add_norm_constraints(self, w: float = 1.0): + pass + + @abstractmethod + def add_tangent_constraints(self, w: float = 1.0): + pass + + @abstractmethod + def add_interface_constraints(self, w: float = 1.0): + pass + def to_dict(self): return { "type": self.type, diff --git a/LoopStructural/interpolators/_operator.py b/LoopStructural/interpolators/_operator.py index ce1489bb6..ed50d61d1 100644 --- a/LoopStructural/interpolators/_operator.py +++ b/LoopStructural/interpolators/_operator.py @@ -36,14 +36,3 @@ class Operator(object): [[0, 0, 0], [0, 1, 0], [0, 0, 0]], # third plane ] ) - - # Hessian is - # Dxx Dxy Dxz - # Dxy Dyy Dzy - # Dzx Dzy Dzz - - # det of hessian - # Det_Hessian_operator = - # Dxx_mask.dot(Dyy_mask.dot(Dzz_mask)) - Dxx_mask.dot(Dyz_mask.dot(Dyz_mask) - - # Dxy_mask.dot(Dxy_mask.dot(Dzz_mask)) - Dxy_mask.dot(Dxz_mask.dot(Dyz_mask)) + - # Dxz_mask.dot(Dxy_mask.dot(Dyz_mask)) - Dxz_mask.dot(Dyy_mask.dot(Dxz_mask)) diff --git a/LoopStructural/interpolators/_p1interpolator.py b/LoopStructural/interpolators/_p1interpolator.py index 7ad7f1d25..931e5cb96 100644 --- a/LoopStructural/interpolators/_p1interpolator.py +++ b/LoopStructural/interpolators/_p1interpolator.py @@ -38,10 +38,10 @@ def __init__(self, mesh): "ipw": 1.0, } - def add_gradient_ctr_pts(self, w=1.0): + def add_gradient_constraints(self, w=1.0): pass - def add_norm_ctr_pts(self, w=1.0): + 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]) @@ -64,7 +64,7 @@ def add_norm_ctr_pts(self, w=1.0): self.up_to_date = False pass - def add_ctr_pts(self, w=1.0): + def add_value_constraints(self, w=1.0): points = self.get_value_constraints() if points.shape[0] > 1: N, elements, inside = self.support.evaluate_shape(points[:, :3]) @@ -165,9 +165,9 @@ def setup_interpolator(self, **kwargs): "%i tangent constraints and %i value constraints" % (self.n_g, self.n_n, self.n_t, self.n_i) ) - self.add_gradient_ctr_pts(self.interpolation_weights["gpw"]) - self.add_norm_ctr_pts(self.interpolation_weights["npw"]) - self.add_ctr_pts(self.interpolation_weights["cpw"]) + self.add_gradient_constraints(self.interpolation_weights["gpw"]) + self.add_norm_constraints(self.interpolation_weights["npw"]) + self.add_value_constraints(self.interpolation_weights["cpw"]) self.add_tangent_constraints(self.interpolation_weights["tpw"]) # self.add_interface_constraints(self.interpolation_weights["ipw"]) @@ -213,3 +213,6 @@ def add_gradient_orthogonal_constraints( gradient constraints not added: outside of model bounding box" ) self.up_to_date = False + + def add_interface_constraints(self, w: float = 1): + raise NotImplementedError diff --git a/LoopStructural/interpolators/_p2interpolator.py b/LoopStructural/interpolators/_p2interpolator.py index aef2259e8..a5bf6e970 100644 --- a/LoopStructural/interpolators/_p2interpolator.py +++ b/LoopStructural/interpolators/_p2interpolator.py @@ -3,7 +3,7 @@ """ import logging -from typing import Optional +from typing import Optional, Callable import numpy as np @@ -92,7 +92,7 @@ def setup_interpolator(self, **kwargs): def copy(self): return P2Interpolator(self.support) - def add_gradient_constraints(self, w=1.0): + def add_gradient_constraints(self, w: float = 1.0): points = self.get_gradient_constraints() if points.shape[0] > 0: grad, elements = self.support.evaluate_shape_derivatives(points[:, :3]) @@ -105,7 +105,9 @@ def add_gradient_constraints(self, w=1.0): elements = self.support[elements[inside]] self.add_constraints_to_least_squares(A * wt[:, None], B, elements, name="gradient") - def add_gradient_orthogonal_constraints(self, points, vector, w=1.0, B=0): + def add_gradient_orthogonal_constraints( + self, points: np.ndarray, vector: np.ndarray, w=1.0, B=0 + ): """ constraints scalar field to be orthogonal to a given vector @@ -133,7 +135,7 @@ def add_gradient_orthogonal_constraints(self, points, vector, w=1.0, B=0): A * wt[:, None], B, elements, name="gradient orthogonal" ) - def add_norm_constraints(self, w=1.0): + def add_norm_constraints(self, w: float = 1.0): points = self.get_norm_constraints() if points.shape[0] > 0: grad, elements = self.support.evaluate_shape_derivatives(points[:, :3]) @@ -150,7 +152,7 @@ def add_norm_constraints(self, w=1.0): name="norm", ) - def add_value_constraints(self, w=1.0): + def add_value_constraints(self, w: float = 1.0): points = self.get_value_constraints() if points.shape[0] > 1: N, elements, mask = self.support.evaluate_shape(points[:, :3]) @@ -167,7 +169,10 @@ def add_value_constraints(self, w=1.0): ) def minimise_grad_steepness( - self, w: float = 0.1, maskall: bool = False, wtfunc: Optional[callable] = None + self, + w: float = 0.1, + maskall: bool = False, + wtfunc: Optional[Callable[[np.ndarray], np.ndarray]] = None, ): """This constraint minimises the second derivative of the gradient mimimising the 2nd derivative should prevent high curvature solutions @@ -206,7 +211,11 @@ def minimise_grad_steepness( ) def minimise_edge_jumps( - self, w: float = 0.1, wtfunc: callable = None, vector_func: callable = None + self, + w: float = 0.1, + wtfunc: Optional[Callable[[np.ndarray], np.ndarray]] = None, + vector_func: Optional[Callable[[np.ndarray], np.ndarray]] = None, + quadrature_points: Optional[int] = None, ): """_summary_ @@ -222,7 +231,7 @@ def minimise_edge_jumps( # NOTE: imposes \phi_T1(xi)-\phi_T2(xi) dot n =0 # iterate over all triangles - cp, weight = self.support.get_quadrature_points(3) + cp, weight = self.support.get_quadrature_points() norm = self.support.shared_element_norm @@ -263,3 +272,6 @@ def evaluate_d2(self, evaluation_points: np.ndarray) -> np.ndarray: if evaluation_points[~mask, :].shape[0] > 0: evaluated[~mask] = self.support.evaluate_d2(evaluation_points[~mask], self.c) return evaluated + + def add_interface_constraints(self, w: float = 1): + raise NotImplementedError diff --git a/LoopStructural/interpolators/supports/_2d_base_unstructured.py b/LoopStructural/interpolators/supports/_2d_base_unstructured.py index 538503666..1fe6cc672 100644 --- a/LoopStructural/interpolators/supports/_2d_base_unstructured.py +++ b/LoopStructural/interpolators/supports/_2d_base_unstructured.py @@ -2,54 +2,94 @@ Tetmesh based on cartesian grid for piecewise linear interpolation """ +from abc import abstractmethod import logging - +from typing import Tuple import numpy as np +from scipy import sparse from . import SupportType +from ._2d_structured_grid import StructuredGrid2D +from ._base_support import BaseSupport +from ._aabb import _initialise_aabb +from ._face_table import _init_face_table logger = logging.getLogger(__name__) -class BaseUnstructured2d: +class BaseUnstructured2d(BaseSupport): """ """ - def __init__(self, elements, vertices, neighbours): + dimension = 2 + + def __init__(self, elements, vertices, neighbours, aabb_nsteps=None): self.type = SupportType.BaseUnstructured2d - self.elements = elements + self._elements = elements self.vertices = vertices if self.elements.shape[1] == 3: self.order = 1 elif self.elements.shape[1] == 6: self.order = 2 - self.nelements = self.elements.shape[0] self.nx = self.vertices.shape[0] - self.n_nodes = self.nx self.neighbours = neighbours + self.minimum = np.min(self.nodes, axis=0) + self.maximum = np.max(self.nodes, axis=0) + length = self.maximum - self.minimum + self.minimum -= length * 0.1 + self.maximum += length * 0.1 + if aabb_nsteps is None: + box_vol = np.prod(self.maximum - self.minimum) + element_volume = box_vol / (len(self.elements) / 20) + # calculate the step vector of a regular cube + step_vector = np.zeros(2) + step_vector[:] = element_volume ** (1.0 / 2.0) + # number of steps is the length of the box / step vector + aabb_nsteps = np.ceil((self.maximum - self.minimum) / step_vector).astype(int) + # make sure there is at least one cell in every dimension + aabb_nsteps[aabb_nsteps < 2] = 2 + step_vector = (self.maximum - self.minimum) / (aabb_nsteps - 1) + self.aabb_grid = StructuredGrid2D(self.minimum, nsteps=aabb_nsteps, step_vector=step_vector) + # make a big table to store which tetra are in which element. + # if this takes up too much memory it could be simplified by using sparse matrices or dict but + # at the expense of speed + self.aabb_table = sparse.csr_matrix( + (self.aabb_grid.n_elements, len(self.elements)), dtype=bool + ) + self.shared_element_relationships = np.zeros( + (self.neighbours[self.neighbours >= 0].flatten().shape[0], 2), dtype=int + ) + self.shared_elements = np.zeros( + (self.neighbours[self.neighbours >= 0].flatten().shape[0], self.dimension), dtype=int + ) + + _init_face_table(self) + _initialise_aabb(self) + + @property + def elements(self): + return self._elements + + def onGeometryChange(self): + pass - # build an array of edges and edge relationships - self.edges = np.zeros((self.nelements * 3, 2), dtype=int) - self.edge_relationships = np.zeros((self.nelements * 3, 2), dtype=int) - edge_index = 0 - flag = np.zeros(self.elements.shape[0], dtype=bool) - for i, t in enumerate(self.elements): - flag[i] = True - for n in self.neighbours[i]: - if n < 0: - continue - if flag[n]: - continue - edge_node_index = 0 - self.edge_relationships[edge_index, 0] = i - self.edge_relationships[edge_index, 1] = n - for v in t: - if v in self.elements[n, :3]: - self.edges[edge_index, edge_node_index] = v - edge_node_index += 1 - - edge_index += 1 - self.edges = self.edges[:edge_index, :] - self.edge_relationships = self.edge_relationships[:edge_index, :] + @property + def n_elements(self): + return self.elements.shape[0] + + @property + def n_nodes(self): + return self.vertices.shape[0] + + def inside(self, pos): + if pos.shape[1] > self.dimension: + logger.warning(f"Converting {pos.shape[1]} to 3d using first {self.dimension} columns") + pos = pos[:, : self.dimension] + + inside = np.ones(pos.shape[0]).astype(bool) + for i in range(self.dimension): + inside *= pos[:, i] > self.origin[None, i] + inside *= pos[:, i] < self.maximum[None, i] + return inside @property def ncps(self): @@ -85,11 +125,55 @@ def barycentre(self): ------- """ - element_idx = np.arange(0, self.nelements) + element_idx = np.arange(0, self.n_elements) elements = self.elements[element_idx] barycentre = np.sum(self.nodes[elements][:, :3, :], axis=1) / 3.0 return barycentre + @property + def shared_element_norm(self): + """ + Get the normal to all of the shared elements + """ + elements = self.shared_elements + v1 = self.nodes[elements[:, 1], :] - self.nodes[elements[:, 0], :] + norm = np.zeros_like(v1) + norm[:, 0] = v1[:, 1] + norm[:, 1] = -v1[:, 0] + return norm + + @property + def shared_element_size(self): + """ + Get the size of the shared elements + """ + elements = self.shared_elements + v1 = self.nodes[elements[:, 1], :] - self.nodes[elements[:, 0], :] + return np.linalg.norm(v1, axis=1) + + @property + def element_size(self): + v1 = self.nodes[self.elements[:, 1], :] - self.nodes[self.elements[:, 0], :] + v2 = self.nodes[self.elements[:, 2], :] - self.nodes[self.elements[:, 0], :] + # cross product isn't defined in 2d, numpy returns the magnitude of the orthogonal vector. + return 0.5 * np.cross(v1, v2, axisa=1, axisb=1) + + @abstractmethod + def evaluate_shape(self, locations) -> Tuple[np.ndarray, np.ndarray, np.ndarray]: + """ + Evaluate the shape functions at the locations + + Parameters + ---------- + locations - numpy array + locations to evaluate + + Returns + ------- + + """ + pass + def element_area(self, elements): tri_points = self.nodes[self.elements[elements, :], :] M_t = np.ones((tri_points.shape[0], 3, 3)) @@ -97,7 +181,7 @@ def element_area(self, elements): area = np.abs(np.linalg.det(M_t)) * 0.5 return area - def evaluate_value(self, pos, values): + def evaluate_value(self, evaluation_points: np.ndarray, property_array: np.ndarray): """ Evaluate value of interpolant @@ -112,15 +196,18 @@ def evaluate_value(self, pos, values): ------- """ - values = np.zeros(pos.shape[0]) - values[:] = np.nan - c, tri = self.evaluate_shape(pos[:, :2]) + pos = np.asarray(evaluation_points) + return_values = np.zeros(pos.shape[0]) + return_values[:] = np.nan + _verts, c, tri, inside = self.get_element_for_location(pos[:, :2]) inside = tri >= 0 # vertices, c, elements, inside = self.get_elements_for_location(pos) - values[inside] = np.sum(c[inside, :] * values[self.elements[tri[inside], :]], axis=1) - return values + return_values[inside] = np.sum( + c[inside, :] * property_array[self.elements[tri[inside], :]], axis=1 + ) + return return_values - def evaluate_gradient(self, pos, prop): + def evaluate_gradient(self, evaluation_points, property_array): """ Evaluate the gradient of an interpolant at the locations @@ -136,18 +223,19 @@ def evaluate_gradient(self, pos, prop): ------- """ - values = np.zeros(pos.shape) + values = np.zeros(evaluation_points.shape) values[:] = np.nan - element_gradients, tri = self.evaluate_shape_derivatives(pos[:, :2]) + element_gradients, tri, inside = self.evaluate_shape_derivatives(evaluation_points[:, :2]) inside = tri >= 0 values[inside, :] = ( - element_gradients[inside, :, :] - * self.properties[prop][self.elements[tri[inside], :, None]] + element_gradients[inside, :, :] * property_array[self.elements[tri[inside], :, None]] ).sum(1) return values - def get_element_for_location(self, pos): + def get_element_for_location( + self, points: np.ndarray + ) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: """ Determine the elements from a numpy array of points @@ -161,18 +249,68 @@ def get_element_for_location(self, pos): ------- """ - M = np.ones((self.elements.shape[0], 3, 3)) - M[:, :, 1:] = self.vertices[self.elements, :][:, :3, :] - points_ = np.ones((pos.shape[0], 3)) - points_[:, 1:] = pos - minv = np.linalg.inv(M) - c = np.einsum("kij,li->lkj", minv, points_) - isinside = np.all(c >= 0, axis=2) - ix, iy = np.where(isinside) - element_idx = np.zeros(pos.shape[0], dtype=int) - element_idx[:] = -1 - element_idx[ix] = iy - c_return = np.zeros((pos.shape[0], 3)) - c_return[:] = np.nan - c_return[ix, :] = c[isinside, :] - return c_return, element_idx + verts = np.zeros((points.shape[0], self.dimension + 1, self.dimension)) + bc = np.zeros((points.shape[0], self.dimension + 1)) + tetras = np.zeros(points.shape[0], dtype="int64") + inside = np.zeros(points.shape[0], dtype=bool) + npts = 0 + npts_step = int(1e4) + # break into blocks of 10k points + while npts < points.shape[0]: + cell_index, inside = self.aabb_grid.position_to_cell_index( + points[: npts + npts_step, :] + ) + global_index = self.aabb_grid.global_cell_indices(cell_index) + tetra_indices = self.aabb_table[global_index[inside], :].tocoo() + # tetra_indices[:] = -1 + row = tetra_indices.row + col = tetra_indices.col + # using returned indexes calculate barycentric coords to determine which tetra the points are in + vertices = self.nodes[self.elements[col, : self.dimension + 1]] + pos = points[row, : self.dimension] + row = tetra_indices.row + col = tetra_indices.col + # using returned indexes calculate barycentric coords to determine which tetra the points are in + vpa = pos[:, :] - vertices[:, 0, :] + vba = vertices[:, 1, :] - vertices[:, 0, :] + vca = vertices[:, 2, :] - vertices[:, 0, :] + d00 = np.einsum('ij,ij->i', vba, vba) + d01 = np.einsum('ij,ij->i', vba, vca) + d11 = np.einsum('ij,ij->i', vca, vca) + d20 = np.einsum('ij,ij->i', vpa, vba) + d21 = np.einsum('ij,ij->i', vpa, vca) + denom = d00 * d11 - d01 * d01 + c = np.zeros((denom.shape[0], 3)) + c[:, 0] = (d11 * d20 - d01 * d21) / denom + c[:, 1] = (d00 * d21 - d01 * d20) / denom + c[:, 2] = 1.0 - c[:, 0] - c[:, 1] + + mask = np.all(c >= 0, axis=1) + + verts[: npts + npts_step, :, :][row[mask], :, :] = vertices[mask, :, :] + bc[: npts + npts_step, :][row[mask], :] = c[mask, :] + tetras[: npts + npts_step][row[mask]] = col[mask] + inside[: npts + npts_step][row[mask]] = True + npts += npts_step + tetra_return = np.zeros((points.shape[0])).astype(int) + tetra_return[:] = -1 + tetra_return[inside] = tetras[inside] + return verts, bc, tetra_return, inside + + def get_element_gradient_for_location( + self, pos: np.ndarray + ) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: + """ + Get the element gradients for a location + + Parameters + ---------- + pos : np.array + location to evaluate + + Returns + ------- + + """ + verts, c, tri, inside = self.get_element_for_location(pos) + return self.evaluate_shape_derivatives(pos, tri) diff --git a/LoopStructural/interpolators/supports/_2d_p1_unstructured.py b/LoopStructural/interpolators/supports/_2d_p1_unstructured.py index 3748587ae..f54172131 100644 --- a/LoopStructural/interpolators/supports/_2d_p1_unstructured.py +++ b/LoopStructural/interpolators/supports/_2d_p1_unstructured.py @@ -22,9 +22,13 @@ def evaluate_shape_derivatives(self, locations, elements=None): """ compute dN/ds (1st row), dN/dt(2nd row) """ + inside = None + if elements is not None: + inside = np.zeros(self.n_elements, dtype=bool) + inside[elements] = True locations = np.array(locations) if elements is None: - c, tri = self.get_element_for_location(locations) + vertices, c, tri, inside = self.get_element_for_location(locations) else: tri = elements M = np.ones((elements.shape[0], 3, 3)) @@ -51,14 +55,14 @@ def evaluate_shape_derivatives(self, locations, elements=None): d_n = np.linalg.inv(jac) # d_n = d_n.swapaxes(1,2) d_n = d_n @ dN - d_n = d_n.swapaxes(2, 1) + # d_n = d_n.swapaxes(2, 1) # d_n = np.dot(np.linalg.inv(jac),dN) - return d_n, tri + return d_n, tri, inside def evaluate_shape(self, locations): locations = np.array(locations) - c, tri = self.get_element_for_location(locations) + vertices, c, tri, inside = self.get_element_for_location(locations) # c = np.dot(np.array([1,x,y]),np.linalg.inv(M)) # convert to barycentric coordinates # order of bary coord is (1-s-t,s,t) N = c # np.zeros((c.shape[0],3)) #evaluate shape functions at barycentric coordinates - return N, tri + return N, tri, inside diff --git a/LoopStructural/interpolators/supports/_2d_p2_unstructured.py b/LoopStructural/interpolators/supports/_2d_p2_unstructured.py index e1909f51d..99688be13 100644 --- a/LoopStructural/interpolators/supports/_2d_p2_unstructured.py +++ b/LoopStructural/interpolators/supports/_2d_p2_unstructured.py @@ -18,7 +18,7 @@ def __init__(self, elements, vertices, neighbours): BaseUnstructured2d.__init__(self, elements, vertices, neighbours) self.type = SupportType.P2Unstructured2d # hessian of shape functions - self.H = np.array( + self.hessian = np.array( [ [[4, 4, 0, 0, 0, -8], [4, 0, 0, 4, -4, -4]], [[4, 0, 0, 4, -4, -4], [4, 0, 4, 0, -8, 0]], @@ -130,43 +130,43 @@ def evaluate_d2_shape(self, indexes): # + self.hN[None, 1, :] * (jac[:, 1, 0] * jac[:, 1, 1])[:, None] # ) - # def evaluate_shape_d2(self, indexes): - # """evaluate second derivatives of shape functions in s and t + def evaluate_shape_d2(self, indexes): + """evaluate second derivatives of shape functions in s and t - # Parameters - # ---------- - # M : [type] - # [description] + Parameters + ---------- + M : [type] + [description] - # Returns - # ------- - # [type] - # [description] - # """ + Returns + ------- + [type] + [description] + """ - # vertices = self.nodes[self.elements[indexes], :] + vertices = self.nodes[self.elements[indexes], :] - # jac = np.array( - # [ - # [ - # (vertices[:, 1, 0] - vertices[:, 0, 0]), - # (vertices[:, 1, 1] - vertices[:, 0, 1]), - # ], - # [ - # vertices[:, 2, 0] - vertices[:, 0, 0], - # vertices[:, 2, 1] - vertices[:, 0, 1], - # ], - # ] - # ).T - # jac = np.linalg.inv(jac) - # jac = jac * jac + jac = np.array( + [ + [ + (vertices[:, 1, 0] - vertices[:, 0, 0]), + (vertices[:, 1, 1] - vertices[:, 0, 1]), + ], + [ + vertices[:, 2, 0] - vertices[:, 0, 0], + vertices[:, 2, 1] - vertices[:, 0, 1], + ], + ] + ).T + jac = np.linalg.inv(jac) + jac = jac * jac - # d2_prod = np.einsum("lij,ik->lik", jac, self.hN) - # d2Const = d2_prod[:, 0, :] + d2_prod[:, 1, :] - # xxConst = d2_prod[:, 0, :] - # yyConst = d2_prod[:, 1, :] + d2_prod = np.einsum("lij,ik->lik", jac, self.hN) + # d2Const = d2_prod[:, 0, :] + d2_prod[:, 1, :] + xxConst = d2_prod[:, 0, :] + yyConst = d2_prod[:, 1, :] - # return xxConst, yyConst + return xxConst, yyConst def evaluate_shape_derivatives(self, locations, elements=None): """ @@ -174,7 +174,7 @@ def evaluate_shape_derivatives(self, locations, elements=None): """ locations = np.array(locations) if elements is None: - c, tri = self.get_element_for_location(locations) + verts, c, tri, inside = self.get_element_for_location(locations) else: tri = elements M = np.ones((elements.shape[0], 3, 3)) @@ -214,13 +214,13 @@ def evaluate_shape_derivatives(self, locations, elements=None): d_n = np.linalg.inv(jac) # d_n = d_n.swapaxes(1,2) d_n = d_n @ dN - d_n = d_n.swapaxes(2, 1) + # d_n = d_n.swapaxes(2, 1) # d_n = np.dot(np.linalg.inv(jac),dN) return d_n, tri def evaluate_shape(self, locations): locations = np.array(locations) - c, tri = self.get_element_for_location(locations) + verts, c, tri, inside = self.get_element_for_location(locations) # c = np.dot(np.array([1,x,y]),np.linalg.inv(M)) # convert to barycentric coordinates # order of bary coord is (1-s-t,s,t) N = np.zeros((c.shape[0], 6)) # evaluate shape functions at barycentric coordinates @@ -231,9 +231,9 @@ def evaluate_shape(self, locations): N[:, 4] = 4 * c[:, 2] * c[:, 0] # 4t(1-s-t) N[:, 5] = 4 * c[:, 1] * c[:, 0] # 4s(1-s-t) - return N, tri + return N, tri, inside - def evaluate_d2(self, pos, prop): + def evaluate_d2(self, pos, property_array): """ Evaluate value of interpolant @@ -256,15 +256,15 @@ def evaluate_d2(self, pos, prop): inside = tri > 0 # vertices, c, elements, inside = self.get_elements_for_location(pos) values[inside] = np.sum( - xxConst[inside, :] * self.properties[prop][self.elements[tri[inside], :]], + xxConst[inside, :] * property_array[self.elements[tri[inside], :]], axis=1, ) values[inside] += np.sum( - yyConst[inside, :] * self.properties[prop][self.elements[tri[inside], :]], + yyConst[inside, :] * property_array[self.elements[tri[inside], :]], axis=1, ) values[inside] += np.sum( - xyConst[inside, :] * self.properties[prop][self.elements[tri[inside], :]], + xyConst[inside, :] * property_array[self.elements[tri[inside], :]], axis=1, ) @@ -272,16 +272,16 @@ def evaluate_d2(self, pos, prop): def get_quadrature_points(self, npts=2): if npts == 2: - v1 = self.nodes[self.edges][:, 0, :] - v2 = self.nodes[self.edges][:, 1, :] + v1 = self.nodes[self.shared_elements][:, 0, :] + v2 = self.nodes[self.shared_elements][:, 1, :] cp = np.zeros((v1.shape[0], self.ncps, 2)) cp[:, 0] = 0.25 * v1 + 0.75 * v2 cp[:, 1] = 0.75 * v1 + 0.25 * v2 - return cp + return cp, np.ones(cp.shape) raise NotImplementedError("Only 2 point quadrature is implemented") def get_edge_normal(self, e): - v = self.nodes[self.edges][:, 0, :] - self.nodes[self.edges][:, 1, :] + v = self.nodes[self.shared_elements][:, 0, :] - self.nodes[self.shared_elements][:, 1, :] # e_len = np.linalg.norm(v, axis=1) normal = np.array([v[:, 1], -v[:, 0]]).T normal /= np.linalg.norm(normal, axis=1)[:, None] diff --git a/LoopStructural/interpolators/supports/_2d_structured_grid.py b/LoopStructural/interpolators/supports/_2d_structured_grid.py index 9c24d4566..4766a93f1 100644 --- a/LoopStructural/interpolators/supports/_2d_structured_grid.py +++ b/LoopStructural/interpolators/supports/_2d_structured_grid.py @@ -5,15 +5,19 @@ import logging +from typing import Tuple import numpy as np from . import SupportType +from ._base_support import BaseSupport logger = logging.getLogger(__name__) -class StructuredGrid2D: +class StructuredGrid2D(BaseSupport): """ """ + dimension = 2 + def __init__( self, origin=np.zeros(2), @@ -34,13 +38,11 @@ def __init__( self.origin = np.array(origin) self.maximum = origin + self.nsteps * self.step_vector - self.n_nodes = self.nsteps[0] * self.nsteps[1] self.dim = 2 self.nsteps_cells = self.nsteps - 1 self.n_cell_x = self.nsteps[0] - 1 self.n_cell_y = self.nsteps[1] - 1 self.properties = {} - self.n_elements = self.n_cell_x * self.n_cell_y # calculate the node positions using numpy (this should probably not # be stored as it defeats @@ -59,6 +61,18 @@ def nodes(self): xx, yy = np.meshgrid(x, y, indexing="ij") return np.array([xx.flatten(order="F"), yy.flatten(order="F")]).T + @property + def n_nodes(self): + return self.nsteps[0] * self.nsteps[1] + + @property + def n_elements(self): + return self.nsteps_cells[0] * self.nsteps_cells[1] + + @property + def element_size(self): + return np.prod(self.step_vector) + @property def barycentre(self): return self.cell_centres(np.arange(self.n_elements)) @@ -66,6 +80,12 @@ def barycentre(self): # @property # def barycentre(self): # return self.cell_centres(np.arange(self.n_elements)) + @property + def elements(self) -> np.ndarray: + global_index = np.arange(self.n_elements) + cell_indexes = self.global_index_to_cell_index(global_index) + + return self.global_node_indices(self.cell_corner_indexes(cell_indexes)) def print_geometry(self): print("Origin: %f %f %f" % (self.origin[0], self.origin[1], self.origin[2])) @@ -75,7 +95,7 @@ def print_geometry(self): max = self.origin + self.nsteps_cells * self.step_vector print("Max extent: %f %f %f" % (max[0], max[1], max[2])) - def cell_centres(self, global_index): + def cell_centres(self, global_index: np.ndarray) -> np.ndarray: """[summary] [extended_summary] @@ -90,12 +110,22 @@ def cell_centres(self, global_index): [type] [description] """ - ix, iy = self.global_index_to_cell_index(global_index) - x = self.origin[None, 0] + self.step_vector[None, 0] * 0.5 + self.step_vector[None, 0] * ix - y = self.origin[None, 1] + self.step_vector[None, 1] * 0.5 + self.step_vector[None, 1] * iy - return np.array([x, y]).T + cell_indexes = self.global_index_to_cell_index(global_index) + cell_centres = np.zeros((cell_indexes.shape[0], 2)) + + cell_centres[:, 0] = ( + self.origin[None, 0] + + self.step_vector[None, 0] * 0.5 + + self.step_vector[None, 0] * cell_indexes[:, 0] + ) + cell_centres[:, 1] = ( + self.origin[None, 1] + + self.step_vector[None, 1] * 0.5 + + self.step_vector[None, 1] * cell_indexes[:, 1] + ) + return cell_centres - def position_to_cell_index(self, pos): + def position_to_cell_index(self, pos: np.ndarray) -> Tuple[np.ndarray, np.ndarray]: """[summary] [extended_summary] @@ -110,16 +140,14 @@ def position_to_cell_index(self, pos): [type] [description] """ + inside = self.inside(pos) + cell_indexes = np.zeros((pos.shape[0], 2)) + cell_indexes[:, 0] = pos[:, 0] - self.origin[None, 0] + cell_indexes[:, 1] = pos[:, 1] - self.origin[None, 1] + cell_indexes /= self.step_vector[None, :] + return cell_indexes.astype(int), inside - pos = self.check_position(pos) - - ix = pos[:, 0] - self.origin[None, 0] - iy = pos[:, 1] - self.origin[None, 1] - ix = ix // self.step_vector[None, 0] - iy = iy // self.step_vector[None, 1] - return ix.astype(int), iy.astype(int) - - def inside(self, pos): + def inside(self, pos: np.ndarray) -> np.ndarray: # check whether point is inside box inside = np.ones(pos.shape[0]).astype(bool) for i in range(self.dim): @@ -130,7 +158,7 @@ def inside(self, pos): ) return inside - def check_position(self, pos): + def check_position(self, pos: np.ndarray) -> np.ndarray: """[summary] [extended_summary] @@ -149,13 +177,13 @@ def check_position(self, pos): if len(pos.shape) == 1: pos = np.array([pos]) if len(pos.shape) != 2: - print("Position array needs to be a list of points or a point") - return False + raise ValueError("Position array needs to be a list of points or a point") + return pos - def bilinear(self, x, y): + def bilinear(self, local_coords: np.ndarray) -> np.ndarray: """ - returns the trilinear interpolation for the local coordinates + returns the bilinear interpolation for the local coordinates Parameters ---------- x - double, array of doubles @@ -167,9 +195,17 @@ def bilinear(self, x, y): array of interpolation coefficients """ - return np.array([(1 - x) * (1 - y), x * (1 - y), (1 - x) * y, x * y]) - def position_to_local_coordinates(self, pos): + return np.array( + [ + (1 - local_coords[:, 0]) * (1 - local_coords[:, 1]), + local_coords[:, 0] * (1 - local_coords[:, 1]), + (1 - local_coords[:, 0]) * local_coords[:, 1], + local_coords[:, 0] * local_coords[:, 1], + ] + ).T + + def position_to_local_coordinates(self, pos: np.ndarray) -> np.ndarray: """ Convert from global to local coordinates within a cel Parameters @@ -184,16 +220,17 @@ def position_to_local_coordinates(self, pos): # TODO check if inside mesh # calculate local coordinates for positions - local_x = ( + local_coords = np.zeros(pos.shape) + local_coords[:, 0] = ( (pos[:, 0] - self.origin[None, 0]) % self.step_vector[None, 0] ) / self.step_vector[None, 0] - local_y = ( + local_coords[:, 1] = ( (pos[:, 1] - self.origin[None, 1]) % self.step_vector[None, 1] ) / self.step_vector[None, 1] - return local_x, local_y + return local_coords - def position_to_dof_coefs(self, pos): + def position_to_dof_coefs(self, pos: np.ndarray): """ global posotion to interpolation coefficients Parameters @@ -204,25 +241,10 @@ def position_to_dof_coefs(self, pos): ------- """ - x_local, y_local = self.position_to_local_coordinates(pos) - weights = self.bilinear(x_local, y_local) + local_coords = self.position_to_local_coordinates(pos) + weights = self.bilinear(local_coords) return weights - def global_indicies(self, indexes): - """ - xi, yi, zi to global index - - Parameters - ---------- - indexes - - Returns - ------- - - """ - indexes = np.array(indexes).swapaxes(0, 2) - return indexes[:, :, 0] + self.nsteps[None, None, 0] * indexes[:, :, 1] - def neighbour_global_indexes(self, mask=None, **kwargs): """ Get neighbour indexes @@ -258,7 +280,7 @@ def neighbour_global_indexes(self, mask=None, **kwargs): np.int64 ) - def cell_corner_indexes(self, x_cell_index, y_cell_index): + def cell_corner_indexes(self, cell_indexes: np.ndarray) -> np.ndarray: """ Returns the indexes of the corners of a cell given its location xi, yi, zi @@ -273,11 +295,16 @@ def cell_corner_indexes(self, x_cell_index, y_cell_index): ------- """ + corner_indexes = np.zeros((cell_indexes.shape[0], 4, 2), dtype=np.int64) xcorner = np.array([0, 1, 0, 1]) ycorner = np.array([0, 0, 1, 1]) - xcorners = x_cell_index[:, None] + xcorner[None, :] - ycorners = y_cell_index[:, None] + ycorner[None, :] - return xcorners, ycorners + corner_indexes[:, :, 0] = ( + cell_indexes[:, None, 0] + corner_indexes[:, :, 0] + xcorner[None, :] + ) + corner_indexes[:, :, 1] = ( + cell_indexes[:, None, 1] + corner_indexes[:, :, 1] + ycorner[None, :] + ) + return corner_indexes def global_index_to_cell_index(self, global_index): """ @@ -294,27 +321,52 @@ def global_index_to_cell_index(self, global_index): # determine the ijk indices for the global index. # remainder when dividing by nx = i # remained when dividing modulus of nx by ny is j - - x_index = global_index % self.nsteps_cells[0, None] - y_index = global_index // self.nsteps_cells[0, None] % self.nsteps_cells[1, None] - return x_index, y_index - - def node_indexes_to_position(self, xindex, yindex): - x = self.origin[0] + self.step_vector[0] * xindex - y = self.origin[1] + self.step_vector[1] * yindex - - return x, y + cell_indexes = np.zeros((global_index.shape[0], 2), dtype=np.int64) + cell_indexes[:, 0] = global_index % self.nsteps_cells[0, None] + cell_indexes[:, 1] = global_index // self.nsteps_cells[0, None] % self.nsteps_cells[1, None] + return cell_indexes + + def global_index_to_node_index(self, global_index): + cell_indexes = np.zeros((global_index.shape[0], 2), dtype=np.int64) + cell_indexes[:, 0] = global_index % self.nsteps[0, None] + cell_indexes[:, 1] = global_index // self.nsteps[0, None] % self.nsteps[1, None] + return cell_indexes + + def _global_indices(self, indexes: np.ndarray, nsteps: np.ndarray) -> np.ndarray: + if len(indexes.shape) == 1: + raise ValueError("Indexes must be a 2D array") + if indexes.shape[-1] != 2: + raise ValueError("Last dimension of cell indexing needs to be ijk indexing") + original_shape = indexes.shape + indexes = indexes.reshape(-1, 2) + gi = indexes[:, 0] + nsteps[0] * indexes[:, 1] + return gi.reshape(original_shape[:-1]) + + def global_cell_indices(self, indexes: np.ndarray) -> np.ndarray: + return self._global_indices(indexes, self.nsteps_cells) + + def global_node_indices(self, indexes: np.ndarray) -> np.ndarray: + return self._global_indices(indexes, self.nsteps) + + def node_indexes_to_position(self, node_indexes: np.ndarray) -> np.ndarray: + + original_shape = node_indexes.shape + node_indexes = node_indexes.reshape((-1, 2)) + xy = np.zeros((node_indexes.shape[0], 2), dtype=float) + xy[:, 0] = self.origin[0] + self.step_vector[0] * node_indexes[:, 0] + xy[:, 1] = self.origin[1] + self.step_vector[1] * node_indexes[:, 1] + xy = xy.reshape(original_shape) + return xy def position_to_cell_corners(self, pos): - inside = self.inside(pos) - ix, iy = self.position_to_cell_index(pos) - cornersx, cornersy = self.cell_corner_indexes(ix, iy) - globalidx = self.global_indicies(np.dstack([cornersx, cornersy]).T) + corner_index, inside = self.position_to_cell_index(pos) + corners = self.cell_corner_indexes(corner_index) + globalidx = self.global_node_indices(corners) # if global index is not inside the support set to -1 globalidx[~inside] = -1 return globalidx, inside - def evaluate_value(self, evaluation_points, property): + def evaluate_value(self, evaluation_points: np.ndarray, property_array: np.ndarray): """ Evaluate the value of of the property at the locations. Trilinear interpolation dot corner values @@ -332,55 +384,79 @@ def evaluate_value(self, evaluation_points, property): v = np.zeros(idc.shape) v[:, :] = np.nan - v[inside, :] = self.position_to_dof_coefs(evaluation_points[inside, :]).T - v[inside, :] *= property[idc[inside, :]] + v[inside, :] = self.position_to_dof_coefs(evaluation_points[inside, :]) + v[inside, :] *= property_array[idc[inside, :]] return np.sum(v, axis=1) - def evaluate_gradient(self, evaluation_points, property): - idc, inside = self.position_to_cell_corners(evaluation_points) - T = np.zeros((idc.shape[0], 2, 4)) - T[inside, :, :] = self.calcul_T(evaluation_points[inside, :]) + def evaluate_gradient(self, evaluation_points, property_array): + T = np.zeros((evaluation_points.shape[0], 2, 4)) + _vertices, T, elements, inside = self.get_element_gradient_for_location(evaluation_points) # indices = np.array([self.position_to_cell_index(evaluation_points)]) # idc = self.global_indicies(indices.swapaxes(0,1)) # print(idc) - T[inside, 0, :] *= property[idc[inside, :]] - T[inside, 1, :] *= property[idc[inside, :]] + T[inside, 0, :] *= property_array[self.elements[elements[inside]]] + T[inside, 1, :] *= property_array[self.elements[elements[inside]]] # T[inside, 2, :] *= self.properties[property_name][idc[inside, :]] return np.array([np.sum(T[:, 0, :], axis=1), np.sum(T[:, 1, :], axis=1)]).T - def calcul_T(self, pos): + def get_element_gradient_for_location( + self, pos + ) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: """ Calculates the gradient matrix at location pos :param pos: numpy array of location Nx3 :return: Nx3x4 matrix """ - # 6_ _ _ _ 8 - # /| /| - # 4 /_| 5/ | - # | 2|_ _|_| 7 - # | / | / - # |/_ _ _|/ - # 0 1 - # - # xindex, yindex, zindex = self.position_to_cell_index(pos) - # cellx, celly, cellz = self.cell_corner_indexes(xindex, yindex,zindex) - # x, y, z = self.node_indexes_to_position(cellx, celly, cellz) + pos = np.asarray(pos) T = np.zeros((pos.shape[0], 2, 4)) - x, y = self.position_to_local_coordinates(pos) - - xindex, yindex = self.position_to_cell_index(pos) - cellx, celly = self.cell_corner_indexes(xindex, yindex) - x, y = self.node_indexes_to_position(cellx, celly) - # div = self.step_vector[0] * self.step_vector[1] * self.step_vector[2] - div = self.step_vector[0] * self.step_vector[1] - - T[:, 0, 0] = -(y[:, 2] - pos[:, 1]) / div - T[:, 0, 1] = (y[:, 3] - pos[:, 1]) / div - T[:, 0, 2] = -(pos[:, 1] - y[:, 0]) / div - T[:, 0, 3] = (pos[:, 1] - y[:, 1]) / div - - T[:, 1, 0] = -(x[:, 1] - pos[:, 0]) / div - T[:, 1, 1] = -(pos[:, 0] - x[:, 0]) / div - T[:, 1, 2] = (x[:, 3] - pos[:, 0]) / div - T[:, 1, 3] = (pos[:, 0] - x[:, 2]) / div - return T + local_coords = self.position_to_local_coordinates(pos) + vertices, inside = self.position_to_cell_corners(pos) + elements, inside = self.position_to_cell_index(pos) + elements = self.global_cell_indices(elements) + + T[:, 0, 0] = -(1 - local_coords[:, 1]) + T[:, 0, 1] = 1 - local_coords[:, 1] + T[:, 0, 2] = -local_coords[:, 1] + T[:, 0, 3] = local_coords[:, 1] + + T[:, 1, 0] = -(1 - local_coords[:, 0]) + T[:, 1, 1] = -local_coords[:, 0] + T[:, 1, 2] = 1 - local_coords[:, 0] + T[:, 1, 3] = local_coords[:, 0] + + return vertices, T, elements, inside + + def get_element_for_location( + self, pos: np.ndarray + ) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: + + vertices, inside = self.position_to_cell_vertices(pos) + vertices = np.array(vertices) + # print("ver", vertices.shape) + # vertices = vertices.reshape((vertices.shape[1], 8, 3)) + elements, inside = self.position_to_cell_corners(pos) + elements, inside = self.position_to_cell_index(pos) + elements = self.global_cell_indices(elements) + a = self.position_to_dof_coefs(pos) + return vertices, a, elements, inside + + def position_to_cell_vertices(self, pos): + """Get the vertices of the cell a point is in + + Parameters + ---------- + pos : np.array + Nx3 array of xyz locations + + Returns + ------- + np.array((N,3),dtype=float), np.array(N,dtype=int) + vertices, inside + """ + gi, inside = self.position_to_cell_corners(pos) + + node_indexes = self.global_index_to_node_index(gi.flatten()) + return self.node_indexes_to_position(node_indexes), inside + + def onGeometryChange(self): + pass diff --git a/LoopStructural/interpolators/supports/_3d_base_structured.py b/LoopStructural/interpolators/supports/_3d_base_structured.py index 84ebac94d..4772963eb 100644 --- a/LoopStructural/interpolators/supports/_3d_base_structured.py +++ b/LoopStructural/interpolators/supports/_3d_base_structured.py @@ -1,15 +1,19 @@ from LoopStructural.utils.exceptions import LoopException -from abc import ABCMeta, abstractmethod +from abc import abstractmethod import numpy as np from LoopStructural.utils import getLogger from . import SupportType logger = getLogger(__name__) +from ._base_support import BaseSupport -class BaseStructuredSupport(metaclass=ABCMeta): + +class BaseStructuredSupport(BaseSupport): """ """ + dimension = 3 + def __init__( self, origin=np.zeros(3), @@ -438,6 +442,10 @@ def global_cell_indices(self, indexes) -> np.ndarray: """ return self._global_indicies(indexes, self.nsteps_cells) + @property + def element_size(self): + return np.prod(self.step_vector) + @property def vtk(self): try: diff --git a/LoopStructural/interpolators/supports/_3d_unstructured_tetra.py b/LoopStructural/interpolators/supports/_3d_unstructured_tetra.py index 8f338ef2b..05e0d931e 100644 --- a/LoopStructural/interpolators/supports/_3d_unstructured_tetra.py +++ b/LoopStructural/interpolators/supports/_3d_unstructured_tetra.py @@ -4,19 +4,23 @@ from ast import Tuple + import numpy as np from scipy.sparse import csr_matrix, coo_matrix, tril from . import StructuredGrid from LoopStructural.utils import getLogger from . import SupportType +from ._base_support import BaseSupport logger = getLogger(__name__) -class UnStructuredTetMesh: +class UnStructuredTetMesh(BaseSupport): """ """ + dimension = 3 + def __init__( self, nodes: np.ndarray, @@ -42,17 +46,16 @@ def __init__( force nsteps for aabb, by default None """ self.type = SupportType.UnStructuredTetMesh - self.nodes = np.array(nodes) - if self.nodes.shape[1] != 3: + self._nodes = np.array(nodes) + if self._nodes.shape[1] != 3: raise ValueError("Nodes must be 3D") - self.n_nodes = self.nodes.shape[0] self.neighbours = np.array(neighbours, dtype=np.int64) if self.neighbours.shape[1] != 4: raise ValueError("Neighbours array is too big") - self.elements = np.array(elements, dtype=np.int64) + self._elements = np.array(elements, dtype=np.int64) if self.elements.shape[0] != self.neighbours.shape[0]: raise ValueError("Number of elements and neighbours do not match") - self.barycentre = np.sum(self.nodes[self.elements[:, :4]][:, :, :], axis=1) / 4.0 + self._barycentre = np.sum(self.nodes[self.elements[:, :4]][:, :, :], axis=1) / 4.0 self.minimum = np.min(self.nodes, axis=0) self.maximum = np.max(self.nodes, axis=0) length = self.maximum - self.minimum @@ -84,6 +87,25 @@ def __init__( self._init_face_table() self._initialise_aabb() + @property + def nodes(self): + return self._nodes + + @property + def elements(self): + return self._elements + + @property + def barycentre(self): + return self._barycentre + + @property + def n_nodes(self): + return self.nodes.shape[0] + + def onGeometryChange(self): + pass + def _init_face_table(self): """ Fill table containing elements that share a face, and another diff --git a/LoopStructural/interpolators/supports/_aabb.py b/LoopStructural/interpolators/supports/_aabb.py new file mode 100644 index 000000000..5574cf297 --- /dev/null +++ b/LoopStructural/interpolators/supports/_aabb.py @@ -0,0 +1,77 @@ +import numpy as np +from scipy import sparse + + +def _initialise_aabb(grid): + """assigns the tetras to the grid cells where the bounding box + of the tetra element overlaps the grid cell. + It could be changed to use the separating axis theorem, however this would require + significantly more calculations. (12 more I think).. #TODO test timing + """ + # calculate the bounding box for all tetraherdon in the mesh + # find the min/max extents for xyz + # tetra_bb = np.zeros((grid.elements.shape[0], 19, 3)) + minx = np.min(grid.nodes[grid.elements[:, :4], 0], axis=1) + maxx = np.max(grid.nodes[grid.elements[:, :4], 0], axis=1) + miny = np.min(grid.nodes[grid.elements[:, :4], 1], axis=1) + maxy = np.max(grid.nodes[grid.elements[:, :4], 1], axis=1) + + cell_indexes = grid.aabb_grid.global_index_to_cell_index(np.arange(grid.aabb_grid.n_elements)) + corners = grid.aabb_grid.cell_corner_indexes(cell_indexes) + positions = grid.aabb_grid.node_indexes_to_position(corners) + ## Because we known the node orders just select min/max from each + # coordinate. Use these to check whether the tetra is in the cell + x_boundary = positions[:, [0, 1], 0] + y_boundary = positions[:, [0, 2], 1] + a = np.logical_and( + minx[None, :] > x_boundary[:, None, 0], + minx[None, :] < x_boundary[:, None, 1], + ) # min point between cell + b = np.logical_and( + maxx[None, :] < x_boundary[:, None, 1], + maxx[None, :] > x_boundary[:, None, 0], + ) # max point between cell + c = np.logical_and( + minx[None, :] < x_boundary[:, None, 0], + maxx[None, :] > x_boundary[:, None, 0], + ) # min point < than cell & max point > cell + + x_logic = np.logical_or(np.logical_or(a, b), c) + + a = np.logical_and( + miny[None, :] > y_boundary[:, None, 0], + miny[None, :] < y_boundary[:, None, 1], + ) # min point between cell + b = np.logical_and( + maxy[None, :] < y_boundary[:, None, 1], + maxy[None, :] > y_boundary[:, None, 0], + ) # max point between cell + c = np.logical_and( + miny[None, :] < y_boundary[:, None, 0], + maxy[None, :] > y_boundary[:, None, 0], + ) # min point < than cell & max point > cell + + y_logic = np.logical_or(np.logical_or(a, b), c) + logic = np.logical_and(x_logic, y_logic) + + if grid.dimension == 3: + z_boundary = positions[:, [0, 6], 2] + minz = np.min(grid.nodes[grid.elements[:, :4], 2], axis=1) + maxz = np.max(grid.nodes[grid.elements[:, :4], 2], axis=1) + a = np.logical_and( + minz[None, :] > z_boundary[:, None, 0], + minz[None, :] < z_boundary[:, None, 1], + ) # min point between cell + b = np.logical_and( + maxz[None, :] < z_boundary[:, None, 1], + maxz[None, :] > z_boundary[:, None, 0], + ) # max point between cell + c = np.logical_and( + minz[None, :] < z_boundary[:, None, 0], + maxz[None, :] > z_boundary[:, None, 0], + ) # min point < than cell & max point > cell + + z_logic = np.logical_or(np.logical_or(a, b), c) + logic = np.logical_and(logic, z_logic) + + grid.aabb_table = sparse.csr_matrix(logic) diff --git a/LoopStructural/interpolators/supports/_base_support.py b/LoopStructural/interpolators/supports/_base_support.py new file mode 100644 index 000000000..b465a75f6 --- /dev/null +++ b/LoopStructural/interpolators/supports/_base_support.py @@ -0,0 +1,114 @@ +from abc import ABCMeta, abstractmethod +import numpy as np +from typing import Tuple + + +class BaseSupport(metaclass=ABCMeta): + """ + Base support class + """ + + @abstractmethod + def __init__(self): + """ + This class is the base + """ + + @abstractmethod + def evaluate_value(self, evaluation_points: np.ndarray, property_array: np.ndarray): + """ + Evaluate the value of the support at the evaluation points + """ + pass + + @abstractmethod + def evaluate_gradient(self, evaluation_points: np.ndarray, property_array: np.ndarray): + """ + Evaluate the gradient of the support at the evaluation points + """ + pass + + @abstractmethod + def inside(self, pos): + """ + Check if a position is inside the support + """ + pass + + @abstractmethod + def onGeometryChange(self): + """ + Called when the geometry changes + """ + pass + + @abstractmethod + def get_element_for_location( + self, pos: np.ndarray + ) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: + """ + Get the element for a location + """ + pass + + @abstractmethod + def get_element_gradient_for_location( + self, pos: np.ndarray + ) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]: + pass + + @property + @abstractmethod + def elements(self): + """ + Return the elements + """ + pass + + @property + @abstractmethod + def n_elements(self): + """ + Return the number of elements + """ + pass + + @property + @abstractmethod + def n_nodes(self): + """ + Return the number of points + """ + pass + + @property + @abstractmethod + def nodes(self): + """ + Return the nodes + """ + pass + + @property + @abstractmethod + def barycentre(self): + """ + Return the number of dimensions + """ + pass + + @property + @abstractmethod + def dimension(self): + """ + Return the number of dimensions + """ + pass + + @property + @abstractmethod + def element_size(self): + """ + Return the element size + """ + pass diff --git a/LoopStructural/interpolators/supports/_face_table.py b/LoopStructural/interpolators/supports/_face_table.py new file mode 100644 index 000000000..8ea406775 --- /dev/null +++ b/LoopStructural/interpolators/supports/_face_table.py @@ -0,0 +1,70 @@ +import numpy as np +from scipy import sparse + + +def _init_face_table(grid): + """ + Fill table containing elements that share a face, and another + table that contains the nodes for a face. + """ + # need to identify the shared nodes for pairs of elements + # we do this by creating a sparse matrix that has N rows (number of elements) + # and M columns (number of nodes). + # We then fill the location where a node is in an element with true + # Then we create a table for the pairs of elements in the mesh + # we have the neighbour relationships, which are the 4 neighbours for each element + # create a new table that shows the element index repeated four times + # flatten both of these arrays so we effectively have a table with pairs of neighbours + # disgard the negative neighbours because these are border neighbours + rows = np.tile(np.arange(grid.n_elements)[:, None], (1, grid.dimension + 1)) + elements = grid.elements + neighbours = grid.neighbours + # add array of bool to the location where there are elements for each node + + # use this to determine shared faces + + element_nodes = sparse.coo_matrix( + ( + np.ones(elements.shape[0] * (grid.dimension + 1)), + (rows.ravel(), elements[:, : grid.dimension + 1].ravel()), + ), + shape=(grid.n_elements, grid.n_nodes), + dtype=bool, + ).tocsr() + n1 = np.tile(np.arange(neighbours.shape[0], dtype=int)[:, None], (1, grid.dimension + 1)) + n1 = n1.flatten() + n2 = neighbours.flatten() + n1 = n1[n2 >= 0] + n2 = n2[n2 >= 0] + el_rel = np.zeros((grid.neighbours.flatten().shape[0], 2), dtype=int) + el_rel[:] = -1 + el_rel[np.arange(n1.shape[0]), 0] = n1 + el_rel[np.arange(n1.shape[0]), 1] = n2 + el_rel = el_rel[el_rel[:, 0] >= 0, :] + + # el_rel2 = np.zeros((grid.neighbours.flatten().shape[0], 2), dtype=int) + grid.shared_element_relationships[:] = -1 + el_pairs = sparse.coo_matrix((np.ones(el_rel.shape[0]), (el_rel[:, 0], el_rel[:, 1]))).tocsr() + i, j = sparse.tril(el_pairs).nonzero() + grid.shared_element_relationships[: len(i), 0] = i + grid.shared_element_relationships[: len(i), 1] = j + + grid.shared_element_relationships = grid.shared_element_relationships[ + grid.shared_element_relationships[:, 0] >= 0, : + ] + + faces = element_nodes[grid.shared_element_relationships[:, 0], :].multiply( + element_nodes[grid.shared_element_relationships[:, 1], :] + ) + shared_faces = faces[np.array(np.sum(faces, axis=1) == grid.dimension).flatten(), :] + row, col = shared_faces.nonzero() + row = row[row.argsort()] + col = col[row.argsort()] + shared_face_index = np.zeros((shared_faces.shape[0], grid.dimension), dtype=int) + shared_face_index[:] = -1 + shared_face_index[row.reshape(-1, grid.dimension)[:, 0], :] = col.reshape(-1, grid.dimension) + grid.shared_elements[np.arange(grid.shared_element_relationships.shape[0]), :] = ( + shared_face_index + ) + # resize + grid.shared_elements = grid.shared_elements[: len(grid.shared_element_relationships), :] diff --git a/LoopStructural/modelling/features/_base_geological_feature.py b/LoopStructural/modelling/features/_base_geological_feature.py index e5a03ac8e..bb46b0a02 100644 --- a/LoopStructural/modelling/features/_base_geological_feature.py +++ b/LoopStructural/modelling/features/_base_geological_feature.py @@ -1,14 +1,15 @@ +from abc import ABCMeta, abstractmethod from LoopStructural.modelling.features import FeatureType from LoopStructural.utils import getLogger from LoopStructural.utils.typing import NumericInput -# from LoopStructural import GeologicalModel + import numpy as np logger = getLogger(__name__) -class BaseFeature: +class BaseFeature(metaclass=ABCMeta): """ Base class for geological features. """ @@ -34,13 +35,35 @@ def __init__(self, name: str, model=None, faults: list = [], regions: list = [], self.name = name self.type = FeatureType.BASE self.regions = regions - self.faults = faults + self._faults = [] + if faults: + self.faults = faults self._model = model self.builder = builder self.faults_enabled = True self._min = None self._max = None + @property + def faults(self): + return self._faults + + @faults.setter + def faults(self, faults: list): + _faults = [] + try: + for f in faults: + if not issubclass(type(f), BaseFeature): + raise TypeError("Faults must be a list of BaseFeature") + _faults.append(f) + except TypeError: + logger.error( + f'Faults must be a list of BaseFeature \n Trying to set using {type(faults)}' + ) + raise TypeError("Faults must be a list of BaseFeature") + + self._faults = _faults + def to_json(self): """ Returns a json representation of the geological feature @@ -136,6 +159,7 @@ def __call__(self, xyz): """ return self.evaluate_value(xyz) + @abstractmethod def evaluate_value(self, pos): """ Evaluate the feature at a given position. @@ -195,6 +219,7 @@ def _apply_faults(self, evaluation_points: np.ndarray, reverse: bool = False) -> evaluation_points = f.apply_to_points(evaluation_points, reverse=reverse) return evaluation_points + @abstractmethod def evaluate_gradient(self, pos): """ Evaluate the gradient of the feature at a given position. diff --git a/LoopStructural/modelling/features/builders/_geological_feature_builder.py b/LoopStructural/modelling/features/builders/_geological_feature_builder.py index a8d42bc23..2a90b0274 100644 --- a/LoopStructural/modelling/features/builders/_geological_feature_builder.py +++ b/LoopStructural/modelling/features/builders/_geological_feature_builder.py @@ -212,7 +212,7 @@ def add_data_to_interpolator(self, constrained=False, force_constrained=False, * if constrained: # Change normals to gradients mask = np.all(~np.isnan(data.loc[:, normal_vec_names()]), axis=1) - if mask.shape[0] > 0: + if mask.sum() > 0: data.loc[mask, gradient_vec_names()] = data.loc[mask, normal_vec_names()].to_numpy( float ) @@ -226,7 +226,7 @@ def add_data_to_interpolator(self, constrained=False, force_constrained=False, * if not constrained or force_constrained: # change gradient constraints to normal vector constraints mask = np.all(~np.isnan(data.loc[:, gradient_vec_names()]), axis=1) - if mask.shape[0] > 0: + if mask.sum() > 0: data.loc[mask, normal_vec_names()] = data.loc[mask, gradient_vec_names()].to_numpy( float ) @@ -240,7 +240,7 @@ def add_data_to_interpolator(self, constrained=False, force_constrained=False, * # self.interpolator.reset() mask = ~np.isnan(data.loc[:, val_name()].to_numpy(float)) # add value constraints - if mask.shape[0] > 0: + if mask.sum() > 0: value_data = data.loc[mask[:, 0], xyz_names() + val_name() + weight_name()].to_numpy( float ) @@ -248,7 +248,7 @@ def add_data_to_interpolator(self, constrained=False, force_constrained=False, * # add gradient constraints mask = np.all(~np.isnan(data.loc[:, gradient_vec_names()].to_numpy(float)), axis=1) - if mask.shape[0] > 0: + if mask.sum() > 0: gradient_data = data.loc[ mask, xyz_names() + gradient_vec_names() + weight_name() ].to_numpy(float) @@ -256,7 +256,7 @@ def add_data_to_interpolator(self, constrained=False, force_constrained=False, * # add normal vector data mask = np.all(~np.isnan(data.loc[:, normal_vec_names()].to_numpy(float)), axis=1) - if mask.shape[0] > 0: + if mask.sum() > 0: normal_data = data.loc[mask, xyz_names() + normal_vec_names() + weight_name()].to_numpy( float ) @@ -264,7 +264,7 @@ def add_data_to_interpolator(self, constrained=False, force_constrained=False, * # add tangent data mask = np.all(~np.isnan(data.loc[:, tangent_vec_names()].to_numpy(float)), axis=1) - if mask.shape[0] > 0: + if mask.sum() > 0: tangent_data = data.loc[ mask, xyz_names() + tangent_vec_names() + weight_name() ].to_numpy(float) @@ -272,15 +272,16 @@ def add_data_to_interpolator(self, constrained=False, force_constrained=False, * # add interface constraints mask = np.all(~np.isnan(data.loc[:, interface_name()].to_numpy(float)), axis=1) - if mask.shape[0] > 0: + if mask.sum() > 0: interface_data = data.loc[ mask, xyz_names() + interface_name() + weight_name() ].to_numpy(float) self.interpolator.set_interface_constraints(interface_data) # add inequality constraints mask = np.all(~np.isnan(data.loc[:, inequality_name()].to_numpy(float)), axis=1) - if mask.shape[0] > 0: + if mask.sum() > 0: inequality_data = data.loc[mask, xyz_names() + inequality_name()].to_numpy(float) + print(inequality_data) self.interpolator.set_inequality_constraints(inequality_data) self.data_added = True diff --git a/LoopStructural/modelling/intrusions/intrusion_feature.py b/LoopStructural/modelling/intrusions/intrusion_feature.py index d539cc4b3..4b76b1c89 100644 --- a/LoopStructural/modelling/intrusions/intrusion_feature.py +++ b/LoopStructural/modelling/intrusions/intrusion_feature.py @@ -228,7 +228,11 @@ def interpolate_vertical_thresholds(self, points_coord1, points_coord2): return threshold_values, residual_values, conceptual_values - def evaluate_value(self, points): + def evaluate_gradient(self, pos): + ## LG TODO check whether it can be implemented + raise NotImplementedError("Cannot calculate gradient of Intrusion") + + def evaluate_value(self, pos): """ Computes a distance scalar field to the intrusion contact (isovalue = 0). @@ -244,12 +248,12 @@ def evaluate_value(self, points): self.builder.up_to_date() # compute coordinates values for each evaluated point - intrusion_coord0_pts = self.intrusion_frame[0].evaluate_value(points) - intrusion_coord1_pts = self.intrusion_frame[1].evaluate_value(points) - intrusion_coord2_pts = self.intrusion_frame[2].evaluate_value(points) + intrusion_coord0_pts = self.intrusion_frame[0].evaluate_value(pos) + intrusion_coord1_pts = self.intrusion_frame[1].evaluate_value(pos) + intrusion_coord2_pts = self.intrusion_frame[2].evaluate_value(pos) self.evaluated_points = [ - points, + pos, intrusion_coord0_pts, intrusion_coord1_pts, intrusion_coord2_pts, @@ -276,7 +280,7 @@ def evaluate_value(self, points): if len(self.assisting_faults) > 0: fault = self.assisting_faults.get("structure") weight = self.assisting_faults.get("asymmetry_weight", 1) - evaluation_points_in_fault = fault[0].evaluate_value(points) + evaluation_points_in_fault = fault[0].evaluate_value(pos) c0_maxside_threshold[evaluation_points_in_fault >= 0] = ( c0_maxside_threshold[evaluation_points_in_fault >= 0] * weight ) diff --git a/tests/integration/test_interpolator.py b/tests/integration/test_interpolator.py index a9cde3399..7c1835f69 100644 --- a/tests/integration/test_interpolator.py +++ b/tests/integration/test_interpolator.py @@ -113,6 +113,16 @@ def test_horizontal_layers(interpolatortype, nelements): assert np.all(np.isclose(model["strati"].evaluate_value(data[["X", "Y", "Z"]]), data["val"])) +def test_horizontal_layers(interpolatortype, nelements): + data, bb = load_horizontal() + model = GeologicalModel(bb[0, :], bb[1, :]) + + model.data = data + model.create_and_add_foliation("strati", interpolatortype=interpolatortype, nelements=1e4) + + assert np.all(np.isclose(model["strati"].evaluate_value(data[["X", "Y", "Z"]]), data["val"])) + + if __name__ == "__main__": test_create_model() test_add_data() @@ -126,4 +136,6 @@ def test_horizontal_layers(interpolatortype, nelements): test_model_with_data_outside_of_bounding_box() test_horizontal_layers("FDI", 1000) test_horizontal_layers("PLI", 1000) + test_horizontal_layers("FDI", 1000) + test_horizontal_layers("PLI", 1000) print("ok") diff --git a/tests/unit/interpolator/test_2d_discrete_support.py b/tests/unit/interpolator/test_2d_discrete_support.py index a718c32d6..b9eebde25 100644 --- a/tests/unit/interpolator/test_2d_discrete_support.py +++ b/tests/unit/interpolator/test_2d_discrete_support.py @@ -48,9 +48,9 @@ def test_get_element_2d(): def test_global_to_local_coordinates2d(): grid = StructuredGrid2D() point = np.array([[1.2, 1.5, 1.7]]) - lx, ly = grid.position_to_local_coordinates(point) - assert np.isclose(lx[0], 0.2) - assert np.isclose(ly[0], 0.5) + local_coords = grid.position_to_local_coordinates(point) + assert np.isclose(local_coords[0, 0], 0.2) + assert np.isclose(local_coords[0, 1], 0.5) def test_get_element_outside2d(): diff --git a/tests/unit/modelling/test_geological_feature.py b/tests/unit/modelling/test_geological_feature.py index 1c904fad2..21d8ff3dd 100644 --- a/tests/unit/modelling/test_geological_feature.py +++ b/tests/unit/modelling/test_geological_feature.py @@ -9,9 +9,9 @@ def test_constructors(): # test constructors work and that the types are set correctly - base_feature = BaseFeature("test", None, [], [], None) - assert base_feature.type == FeatureType.BASE - assert base_feature.name == "test" + # base_feature = GeologicalFeature("test", None, [], [], None) + # assert base_feature.type == FeatureType.BASE + # assert base_feature.name == "test" feature = GeologicalFeature("test", None, [], [], None) assert feature.type == FeatureType.INTERPOLATED assert feature.name == "test" @@ -29,7 +29,7 @@ def test_constructors(): def test_toggle_faults(): - base_feature = BaseFeature("test", None, [], [], None) + base_feature = GeologicalFeature("test", None, [], [], None) assert base_feature.faults_enabled is True base_feature.toggle_faults() assert base_feature.faults_enabled is False @@ -38,7 +38,7 @@ def test_toggle_faults(): def test_tojson(): - base_feature = BaseFeature("test", None, [], [], None) + base_feature = GeologicalFeature("test", None, [], [], None) import json from LoopStructural.utils import LoopJSONEncoder