Merge pull request #71 from Loop3D/dataprocessor

Dataprocessor

Author: lachlangrose

DockerfileDev

@@ -46,6 +46,7 @@ USER ${NB_USER}
 RUN mkdir notebooks
 RUN mkdir LoopStructural
 RUN mkdir map2loop
+RUN mkdir LoopProjectFile
 
 # RUN pip install -e LoopStructural
 CMD ["jupyter", "notebook", "--ip='0.0.0.0'", "--NotebookApp.token=''", "--no-browser" ]

LoopStructural/__init__.py

@@ -15,7 +15,7 @@
 ch.setFormatter(formatter)
 ch.setLevel(logging.WARNING)
 loggers = {}
-__version__ = '1.1.0'
+__version__ = '1.2.0'
 from .modelling.core.geological_model import GeologicalModel
 from .utils import log_to_console, log_to_file, getLogger
 logger = getLogger(__name__)

LoopStructural/datasets/data/refolded_fold.csv

@@ -79,58 +79,62 @@ def _setup_interpolator(self, **kwargs):
         for key in kwargs:
             self.up_to_date = False
             if 'regularisation' in kwargs:
-                self.interpolation_weights['dxy'] = kwargs[
-                                                        'regularisation'] * 0.7
-                self.interpolation_weights['dyz'] = kwargs[
-                                                        'regularisation'] * 0.7
-                self.interpolation_weights['dxz'] = kwargs[
-                                                        'regularisation'] * 0.7
-                self.interpolation_weights['dxx'] = kwargs[
-                                                        'regularisation'] * 1.
-                self.interpolation_weights['dyy'] = kwargs[
-                                                        'regularisation'] * 1.
-                self.interpolation_weights['dzz'] = kwargs[
-                                                        'regularisation'] * 1.
+                self.interpolation_weights['dxy'] = 0.1*kwargs[
+                                                        'regularisation'] 
+                self.interpolation_weights['dyz'] = 0.1*kwargs[
+                                                        'regularisation'] 
+                self.interpolation_weights['dxz'] = 0.1*kwargs[
+                                                        'regularisation'] 
+                self.interpolation_weights['dxx'] = 0.1*kwargs[
+                                                        'regularisation'] 
+                self.interpolation_weights['dyy'] = 0.1*kwargs[
+                                                        'regularisation'] 
+                self.interpolation_weights['dzz'] = 0.1*kwargs[
+                                                        'regularisation'] 
             self.interpolation_weights[key] = kwargs[key]
         # if we want to define the operators manually
         if 'operators' in kwargs:
             for n, o in kwargs['operators'].items():
                 self.assemble_inner(o[0], o[1])
         # otherwise just use defaults
         if 'operators' not in kwargs:
+            
             operator = Operator.Dxy_mask
-
-            self.assemble_inner(operator, np.sqrt(2 * self.vol) *
-                                self.interpolation_weights['dxy'])
+            weight =  self.interpolation_weights['dxy'] / \
+                             1#(4*self.support.step_vector[0]*self.support.step_vector[1])
+            self.assemble_inner(operator, weight )
             operator = Operator.Dyz_mask
-            self.assemble_inner(operator, np.sqrt(2 * self.vol) *
-                                self.interpolation_weights['dyz'])
+            weight = self.interpolation_weights['dyz'] / \
+                             1#(4*self.support.step_vector[1]*self.support.step_vector[2])
+            self.assemble_inner(operator, weight)
             operator = Operator.Dxz_mask
-            self.assemble_inner(operator, np.sqrt(2 * self.vol) *
-                                self.interpolation_weights['dxz'])
+            weight =  self.interpolation_weights['dxz'] / \
+                             1#(4*self.support.step_vector[0]*self.support.step_vector[2])
+            self.assemble_inner(operator, weight)
             operator = Operator.Dxx_mask
+            weight = self.interpolation_weights['dxx'] \
+                             / 1#self.support.step_vector[0]**2
             self.assemble_inner(operator,
-                                np.sqrt(self.vol) * self.interpolation_weights[
-                                    'dxx'])
+                                weight)
             operator = Operator.Dyy_mask
-            self.assemble_inner(operator,
-                                np.sqrt(self.vol) * self.interpolation_weights[
-                                    'dyy'])
+            weight =  self.interpolation_weights['dyy'] / \
+                             1#self.support.step_vector[1]**2
+            self.assemble_inner(operator,weight)
             operator = Operator.Dzz_mask
-            self.assemble_inner(operator,
-                                np.sqrt(self.vol) * self.interpolation_weights[
-                                    'dzz'])
+            weight = self.interpolation_weights['dzz'] / \
+                             1#self.support.step_vector[2]**2
+            self.assemble_inner(operator,weight)
         self.add_norm_constraint(
-            np.sqrt(self.vol) * self.interpolation_weights['npw'])
+            self.interpolation_weights['npw'])
         self.add_gradient_constraint(
-            np.sqrt(self.vol) * self.interpolation_weights['gpw'])
+             self.interpolation_weights['gpw'])
         self.add_vaue_constraint(
-            np.sqrt(self.vol) * self.interpolation_weights['cpw'])
+             self.interpolation_weights['cpw'])
         self.add_tangent_ctr_pts(
-            np.sqrt(self.vol) * self.interpolation_weights['tpw']
+            self.interpolation_weights['tpw']
         )
         self.add_interface_ctr_pts(
-            np.sqrt(self.vol)*self.interpolation_weights['ipw']
+            self.interpolation_weights['ipw']
         )
 
     def copy(self):
@@ -172,7 +176,7 @@ def add_vaue_constraint(self, w=1.):
             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)
             self.add_constraints_to_least_squares(a.T * w,
                                                   points[inside, 3] * w,
                                                   idc[inside, :],
@@ -259,6 +263,8 @@ def add_gradient_constraint(self, w=1.):
             inside = np.logical_and(~np.any(idc == -1, axis=1), inside)
 
             T = self.support.calcul_T(points[inside, :3])
+            norm = np.linalg.norm(T,axis=2)
+            T/=norm[:,:,None]
             strike_vector, dip_vector = get_vectors(points[inside, 3:6])
             A = np.einsum('ij,ijk->ik', strike_vector.T, T)
 
@@ -284,7 +290,6 @@ def add_norm_constraint(self, w=1.):
         if points.shape[0] > 0:
             # calculate unit vector for orientation data
             # points[:,3:]/=np.linalg.norm(points[:,3:],axis=1)[:,None]
-
             node_idx, inside = self.support.position_to_cell_corners(
                 points[:, :3])
             gi = np.zeros(self.support.n_nodes)
@@ -299,16 +304,17 @@ def add_norm_constraint(self, w=1.):
             # this means we are only constraining direction of grad not the
             # magnitude
             T = self.support.calcul_T(points[inside, :3])
-
+            # T*=np.product(self.support.step_vector)
+            # T/=self.support.step_vector[0]
             w /= 3
             self.add_constraints_to_least_squares(T[:, 0, :] * w,
-                                                  points[inside, 3] * w,
+                                                  points[inside, 3] * w ,
                                                   idc[inside, :], name='norm')
             self.add_constraints_to_least_squares(T[:, 1, :] * w,
-                                                  points[inside, 4] * w,
+                                                  points[inside, 4] * w ,
                                                   idc[inside, :], name='norm')
             self.add_constraints_to_least_squares(T[:, 2, :] * w,
-                                                  points[inside, 5] * w,
+                                                  points[inside, 5] * w ,
                                                   idc[inside, :], name='norm')
 
     def add_gradient_orthogonal_constraint(self, points, vector, w=1.0,
@@ -345,9 +351,15 @@ def add_gradient_orthogonal_constraint(self, points, vector, w=1.0,
             idc[inside, :] = gi[node_idx[inside, :]]
             inside = np.logical_and(~np.any(idc == -1, axis=1), inside)
 
+            #normalise element vector to unit vector for dot product
             T = self.support.calcul_T(points[inside, :3])
+            norm = np.linalg.norm(T,axis=1)
+            T/=norm[:,None,:]
+            # normalise vector to unit vector for dot product
+            vector[inside,:3] /= np.linalg.norm(vector[inside,:3],axis=1)[:,None]
+            # dot product of vector and element gradient 
             A = np.einsum('ij,ijk->ik', vector[inside, :3], T)
-
+            
             B = np.zeros(points[inside, :].shape[0])
             self.add_constraints_to_least_squares(A * w, B, idc[inside, :], name='gradient orthogonal')
 
@@ -390,9 +402,11 @@ def assemble_inner(self, operator, w):
         gi[:] = -1
         gi[self.region] = np.arange(0, self.nx)
         idc = gi[idc]
-        inside = ~np.any(idc == -1, axis=1)
+        inside = ~np.any(idc == -1, axis=1)#np.ones(a.shape[0],dtype=bool)#
+        # a[idc==-1] = 0
+        # idc[idc==-1] = 0
         B = np.zeros(global_indexes.shape[1])
-        self.add_constraints_to_least_squares(a[inside, :] * w,
+        self.add_constraints_to_least_squares(a[inside, :] * w ,
                                               B[inside],
                                               idc[inside, :],
                                               name='regularisation'

LoopStructural/interpolators/finite_difference_interpolator.py

@@ -33,7 +33,7 @@ class Operator(object):
         [-0.25, 0, 0.25],
         [0, 0, 0],
         [0.25, 0, -0.25]
-    ], z])
+    ], z]) / np.sqrt(2)
     Dxz_mask = Dxy_mask.swapaxes(0, 1)
     Dyz_mask = Dxy_mask.swapaxes(0, 2)
 

LoopStructural/interpolators/operator.py

@@ -290,7 +290,7 @@ def add_norm_ctr_pts(self, w=1.0):
             # d_t = self.support.get_elements_gradients(e)
             norm = np.zeros((element_gradients.shape[0],element_gradients.shape[1]))
             norm[inside,:] = np.linalg.norm(element_gradients[inside,:,:], axis=2)
-            element_gradients /= norm[:, :, None]
+            # element_gradients /= norm[:, :, None]
 
             d_t = element_gradients
             d_t[inside,:,:] *= vol[inside, None, None]

LoopStructural/interpolators/piecewiselinear_interpolator.py

@@ -40,19 +40,18 @@ def barycentre(self):
 
 
     def cell_centres(self, global_index):
-        """[summary]
-
-        [extended_summary]
+        """get the centre of specified cells
 
+        
         Parameters
         ----------
-        global_index : [type]
-            [description]
+        global_index : array/list
+            container of integer global indexes to cells
 
         Returns
         -------
-        [type]
-            [description]
+        numpy array
+            Nx3 array of cell centres
         """
         ix, iy, iz = self.global_index_to_cell_index(global_index)
         x = self.origin[None, 0] + self.step_vector[None, 0] * .5 + \
@@ -103,15 +102,14 @@ def position_to_local_coordinates(self, pos):
 
         """
         # TODO check if inside mesh
-        pos = self.rotate(pos)
+        # pos = self.rotate(pos)
         # calculate local coordinates for positions
         local_x = ((pos[:, 0] - self.origin[None, 0]) % self.step_vector[
             None, 0]) / self.step_vector[None, 0]
         local_y = ((pos[:, 1] - self.origin[None, 1]) % self.step_vector[
             None, 1]) / self.step_vector[None, 1]
         local_z = ((pos[:, 2] - self.origin[None, 2]) % self.step_vector[
             None, 2]) / self.step_vector[None, 2]
-
         return local_x, local_y, local_z
 
     def position_to_dof_coefs(self, pos):
@@ -325,9 +323,8 @@ def calcul_T(self, pos):
         # x, y, z = self.node_indexes_to_position(cellx, celly, cellz)
         T = np.zeros((pos.shape[0], 3, 8))
         x, y, z = self.position_to_local_coordinates(pos)
-        # div = self.step_vector[0] * self.step_vector[1] * self.step_vector[2]
-
-        T[:, 0, 0] = -(1 - y) * (1 - z)  # v000
+    
+        T[:, 0, 0] = (1 - z) * (y- 1)  # v000
         T[:, 0, 1] = (1 - y) * (1 - z)  # (y[:, 3] - pos[:, 1]) / div
         T[:, 0, 2] = -y * (1 - z)  # (pos[:, 1] - y[:, 0]) / div
         T[:, 0, 3] = -(1 - y) * z  # (pos[:, 1] - y[:, 1]) / div
@@ -336,7 +333,7 @@ def calcul_T(self, pos):
         T[:, 0, 6] = y * (1 - z)
         T[:, 0, 7] = y * z
 
-        T[:, 1, 0] = - (1 - x) * (1 - z)
+        T[:, 1, 0] =  (x - 1) * (1 - z)
         T[:, 1, 1] = - x * (1 - z)
         T[:, 1, 2] = (1 - x) * (1 - z)
         T[:, 1, 3] = -(1 - x) * z
@@ -353,14 +350,7 @@ def calcul_T(self, pos):
         T[:, 2, 5] = (1 - x) * y
         T[:, 2, 6] = - x * y
         T[:, 2, 7] = x * y
-        return T
+        T/=self.step_vector[0]
+
+        return T 
 
-    def slice(self, propertyname, isovalue, region):
-        logger.error("function has been removed, please use the modelviewer class")
-        return
-        #
-        # verts, faces, normals, values = marching_cubes(
-        #     self.properties[propertyname].reshape(self.nsteps, order='F'),
-        #     isovalue,
-        #     spacing=self.step_vector)
-        # return faces, verts + self.origin[None, :]

LoopStructural/interpolators/structured_grid.py

@@ -115,4 +115,7 @@ def evaluate_gradient(self, evaluation_points):
         if evaluation_points[~mask, :].shape[0] > 0:
             evaluated[~mask,:] = self.surfe.EvaluateVectorInterpolantAtPoints(
                 evaluation_points[~mask])
-        return evaluated
+        return 
+    @property
+    def nx(self):
+        return self.get_data_locations().shape[0]

LoopStructural/interpolators/surfe_wrapper.py

@@ -3,3 +3,4 @@
 
 """
 
+from .input.process_data import ProcessInputData