Merge branch 'master' of https://github.com/Loop3D/LoopStructural

Author: lachlangrose

.github/workflows/release-please.yml

@@ -6,7 +6,7 @@ jobs:
     name: Linting
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v4
       - name: Install dependencies
         run: |
           python -m pip install --upgrade pip
@@ -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"
 
@@ -30,8 +30,9 @@ jobs:
         os: ${{ fromJSON(vars.BUILD_OS)}}
         python-version: ${{ fromJSON(vars.PYTHON_VERSIONS)}}
     steps:
-      - uses: actions/checkout@v2
-      - uses: conda-incubator/setup-miniconda@v2
+      - uses: actions/checkout@v4
+      - uses: conda-incubator/setup-miniconda@v3
+
         with:
           python-version: ${{ matrix.python }}
       - name: Installing dependencies
@@ -50,7 +51,7 @@ jobs:
     runs-on: ubuntu-latest
     #needs: continuous-integration
     steps:
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v4
       - run: |
           cp CHANGELOG.md docs/source/getting_started/CHANGELOG.md
           docker build . -t=lsdocs -f docs/Dockerfile 
@@ -65,12 +66,12 @@ jobs:
     needs: documentation-test
     if: github.ref == 'refs/heads/master'
     steps:
-      - uses: GoogleCloudPlatform/release-please-action@v3.5.0
+      - uses: GoogleCloudPlatform/release-please-action@v4
         id: release
         with:
           release-type: python
           package-name: LoopStructural
-          #version-file: LoopStructural/version.py
+          version-file: LoopStructural/version.py
     outputs:
       release_created: ${{ steps.release.outputs.release_created }}
   # if a release is created then run the deploy scripts for github.io, conda, pypi and docker
@@ -79,7 +80,7 @@ jobs:
     needs: release-please
     if: ${{ needs.release-please.outputs.release_created }}
     steps:
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v4
       - run: |
           cp CHANGELOG.md docs/source/getting_started/CHANGELOG.md
           docker build . -t=lsdocs -f docs/Dockerfile 
@@ -100,12 +101,12 @@ 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 }}
 
-      - uses: actions/checkout@v2
+      - uses: actions/checkout@v4
       - name: update submodules
         #       shell: bash -l {0}
         run: |
@@ -120,29 +121,29 @@ 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
+          name: conda-${{ matrix.os }}-${{ matrix.python-version }}
           path: conda
 
   make_sdist:
     needs: ["documentation-test", "continuous-integration"]
     name: Make SDist
     runs-on: ubuntu-latest
     steps:
-      - uses: actions/checkout@v3
+      - uses: actions/checkout@v4
 
       - name: Build SDist
         run: |
           pip install numpy
           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,21 +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:
-          name: conda
           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:

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":

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`
 

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):