Merge pull request #146 from Loop3D/intrusions

Intrusions

Author: lachlangrose

.github/workflows/release-please.yml

@@ -18,7 +18,7 @@ jobs:
       - name: Installing dependencies
         shell: bash -l {0}
         run: |
-          conda install -c conda-forge cython numpy scipy scikit-image scikit-learn pyamg flake8 pytest networkx osqp numba statsmodels scikit-fmm matplotlib -y
+          conda install -c conda-forge cython numpy scipy scikit-image scikit-learn pyamg flake8 pytest networkx osqp matplotlib -y
       - name: Checking formatting of code
         shell: bash -l {0}
         run: |

Dockerfile

@@ -13,17 +13,15 @@ RUN apt-get update -qq && \
     libopenmpi-dev \ 
     make 
 # RUN conda install -c conda-forge python=3.9 -y
-RUN conda install -c conda-forge -c loop3d pip \
+RUN conda install -c conda-forge -c loop3d\
+    pip \
     map2model\
     hjson\
     owslib\
     beartype\
     gdal=3.5.2\
     rasterio=1.2.10 \
     meshio\
-    scikit-fmm\
-    statsmodels\
-    numba\
     scikit-learn \
     cython \
     numpy \
@@ -41,11 +39,16 @@ RUN conda install -c conda-forge -c loop3d pip \
     folium \
     jupyterlab \
     nodejs \
+    rasterio\
+    geopandas\
     -y
 
 RUN pip install ipyfilechooser
 RUN jupyter nbextension enable --py --sys-prefix ipyleaflet
-RUN pip install lavavu-osmesa geostatspy mplstereonet
+RUN pip install lavavu-osmesa mplstereonet
+
+ENV LD_LIBRARY_PATH=/opt/conda/lib/python3.10/site-packages/lavavu/
+
 
 ENV NB_USER jovyan
 ENV NB_UID 1000
@@ -80,6 +83,7 @@ RUN pip install -e LoopProjectFile
 # WORKDIR Tomofast-x
 # RUN make
 WORKDIR ${HOME}/notebooks
+
 # RUN pip install -e LoopStructural
 CMD ["jupyter", "lab", "--ip='0.0.0.0'", "--NotebookApp.token=''", "--no-browser" ]
 

DockerfileDev

@@ -13,7 +13,8 @@ RUN apt-get update -qq && \
     libopenmpi-dev \ 
     make 
 # RUN conda install -c conda-forge python=3.9 -y
-RUN conda install -c conda-forge pip \
+RUN conda install -c conda-forge "python<=3.9" \
+    pip \
     scikit-learn \
     cython \
     numpy \
@@ -31,11 +32,12 @@ RUN conda install -c conda-forge pip \
     folium \
     jupyterlab \
     nodejs \
+    rasterio\
     -y
 
 RUN pip install ipyfilechooser
 RUN jupyter nbextension enable --py --sys-prefix ipyleaflet
-RUN pip install lavavu-osmesa
+RUN pip install lavavu-osmesa==1.8.32 pyevtk
 
 ENV NB_USER jovyan
 ENV NB_UID 1000

LoopStructural/interpolators/_discrete_interpolator.py

@@ -7,8 +7,6 @@
 import numpy as np
 from scipy.sparse import coo_matrix, bmat, eye
 from scipy.sparse import linalg as sla
-from scipy.sparse.linalg import norm
-from sklearn.preprocessing import normalize
 from ..interpolators import InterpolatorType
 
 from ..interpolators import GeologicalInterpolator
@@ -88,6 +86,7 @@ def region(self) -> np.ndarray:
         np.ndarray
 
         """
+
         return self.region_function(self.support.nodes).astype(bool)
 
     @property
@@ -126,7 +125,8 @@ def set_region(self, region=None):
         """
         # evaluate the region function on the support to determine
         # which nodes are inside update region map and degrees of freedom
-        self.region_function = region
+        print("Cannot use region")
+        # self.region_function = region
         logger.info(
             "Interpolation now uses region and has {} degrees of freedom".format(
                 self.nx
@@ -800,7 +800,9 @@ def _solve(self, solver="cg", **kwargs):
             )
             return
         self.valid = True
-        logging.info(f"Solving interpolation: {self.propertyname} took: {time()-starttime}")
+        logging.info(
+            f"Solving interpolation: {self.propertyname} took: {time()-starttime}"
+        )
 
     def update(self):
         """

LoopStructural/interpolators/_geological_interpolator.py

@@ -24,15 +24,9 @@ def __init__(self):
         this should be callable by using any of these functions. This will
         enable interpolators to be interchanged.
         """
+        self.data = {}  # None
+        self.clean()  # init data structure
 
-        self.data = {
-            "gradient": np.zeros((0, 7)),
-            "value": np.zeros((0, 5)),
-            "normal": np.zeros((0, 7)),
-            "tangent": np.zeros((0, 7)),
-            "interface": np.zeros((0, 5)),
-            "inequality": np.zeros((0, 6)),
-        }
         self.n_g = 0
         self.n_i = 0
         self.n_n = 0
@@ -200,7 +194,7 @@ def get_data_locations(self):
         numpy array
             Nx3 - X,Y,Z location of all data points
         """
-        return np.vstack([d for d in self.data.values()[:, :3]])
+        return np.vstack([d[:, :3] for d in self.data.values()])
 
     def get_interface_constraints(self):
         """Get the location of interface constraints
@@ -231,14 +225,23 @@ def solve_system(self, **kwargs):
     def update(self):
         return False
 
-    def reset(self):
+    def clean(self):
         """
         Removes all of the data from an interpolator
 
         Returns
         -------
 
         """
+        self.data = {
+            "gradient": np.zeros((0, 7)),
+            "value": np.zeros((0, 5)),
+            "normal": np.zeros((0, 7)),
+            "tangent": np.zeros((0, 7)),
+            "interface": np.zeros((0, 5)),
+            "inequality": np.zeros((0, 6)),
+        }
+        self.up_to_date = False
         self.n_g = 0
         self.n_i = 0
         self.n_n = 0

LoopStructural/interpolators/piecewiselinear_interpolator.py

@@ -186,7 +186,9 @@ def add_constant_gradient(
             gi[:] = -1
             gi[self.region] = np.arange(0, self.nx)
             idc = gi[idc]
+
             outside = ~np.any(idc == -1, axis=1)
+
             # w/=A.shape[0]
             self.add_constraints_to_least_squares(
                 A[outside, :], B[outside], idc[outside, :], w=w, name="regularisation"
@@ -239,12 +241,20 @@ def add_gradient_constraints(self, w=1.0):
             outside = ~np.any(idc == -1, axis=1)
             w *= points[:, 6]
             self.add_constraints_to_least_squares(
-                A[outside, :], B[outside], idc[outside, :], w=w, name="gradient strike"
+                A[outside, :],
+                B[outside],
+                idc[outside, :],
+                w=w[outside],
+                name="gradient strike",
             )
             A = np.einsum("ji,ijk->ik", dip_vector, element_gradients)
             # A *= vol[:, None]
             self.add_constraints_to_least_squares(
-                A[outside, :], B[outside], idc[outside, :], w=w, name="gradient dip"
+                A[outside, :],
+                B[outside],
+                idc[outside, :],
+                w=w[outside],
+                name="gradient dip",
             )
 
     def add_norm_constraints(self, w=1.0):
@@ -299,7 +309,11 @@ def add_norm_constraints(self, w=1.0):
             # w /= 3
 
             self.add_constraints_to_least_squares(
-                d_t[outside, :, :], points[outside, 3:6], idc[outside], w=w, name="norm"
+                d_t[outside, :, :],
+                points[outside, 3:6],
+                idc[outside],
+                w=w[outside],
+                name="norm",
             )
 
     def add_value_constraints(self, w=1.0):  # for now weight all value points the same
@@ -339,9 +353,9 @@ def add_value_constraints(self, w=1.0):  # for now weight all value points the s
             w *= points[inside, 4]
             self.add_constraints_to_least_squares(
                 A[outside, :],
-                points[inside, :][outside, 3] * vol[:],
+                points[inside, :][outside, 3] * vol[outside],
                 idc[outside, :],
-                w=w,
+                w=w[outside],
                 name="value",
             )
 

LoopStructural/interpolators/supports/_3d_base_structured.py

@@ -218,7 +218,7 @@ def check_position(self, pos):
             return False
         return pos
 
-    def _global_indicies(self, indexes, nsteps):
+    def _global_indicies(self, indexes: np.ndarray, nsteps: np.ndarray) -> np.ndarray:
         """
         Convert from cell indexes to global cell index
 
@@ -285,7 +285,9 @@ def position_to_cell_corners(self, pos):
         inside = self.inside(pos)
         ix, iy, iz = self.position_to_cell_index(pos)
         cornersx, cornersy, cornersz = self.cell_corner_indexes(ix, iy, iz)
-        globalidx = self.global_indicies(np.dstack([cornersx, cornersy, cornersz]).T)
+        globalidx = self.global_node_indicies(
+            np.dstack([cornersx, cornersy, cornersz]).T
+        )
         # if global index is not inside the support set to -1
         globalidx[~inside] = -1
         return globalidx, inside
@@ -360,7 +362,7 @@ def global_index_to_node_index(self, global_index):
         z_index = global_index // self.nsteps[0, None] // self.nsteps[1, None]
         return x_index, y_index, z_index
 
-    def global_node_indicies(self, indexes):
+    def global_node_indicies(self, indexes) -> np.ndarray:
         """
         Convert from node indexes to global node index
 
@@ -374,7 +376,7 @@ def global_node_indicies(self, indexes):
         """
         return self._global_indicies(indexes, self.nsteps)
 
-    def global_cell_indicies(self, indexes):
+    def global_cell_indicies(self, indexes) -> np.ndarray:
         """
         Convert from cell indexes to global cell index
 

LoopStructural/interpolators/supports/_3d_structured_grid.py

@@ -141,25 +141,6 @@ def position_to_dof_coefs(self, pos):
         weights = self.trilinear(x_local, y_local, local_z)
         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]
-            + self.nsteps[None, None, 0] * self.nsteps[None, None, 1] * indexes[:, :, 2]
-        )
-
     def neighbour_global_indexes(self, mask=None, **kwargs):
         """
         Get neighbour indexes
@@ -475,6 +456,11 @@ def get_element_for_location(self, pos):
             [description]
         """
         vertices, inside = self.position_to_cell_vertices(pos)
+        vertices = np.array(vertices)
+        vertices = vertices.reshape((vertices.shape[1], 8, 3))
         elements, inside = self.position_to_cell_corners(pos)
         a = self.position_to_dof_coefs(pos)
         return vertices, a.T, elements, inside
+
+    def get_elements(self):
+        return