Merge pull request #31 from Loop3D/lg_dev

Lg dev

LoopStructural/__init__.py

@@ -30,4 +30,4 @@
 #     # temp_file = tempfile.tempdir+Path('/default-loop-structural-logfile.log')
 #     log_to_file(temp_file)
 log_to_console()
-__version__ = '1.0.1'
+__version__ = '1.0.3'

LoopStructural/interpolators/discrete_interpolator.py

@@ -148,7 +148,7 @@ def add_constraints_to_least_squares(self, A, B, idc, name='undefined'):
         B = np.array(B)
         idc = np.array(idc)
         if np.any(np.isnan(idc)) or np.any(np.isnan(A)) or np.any(np.isnan(B)):
-            logger.warning("Constraints contain nan not adding constraints")
+            logger.warning("Constraints contain nan not adding constraints: {}".format(name))
             return
         nr = A.shape[0]
         if len(A.shape) > 2:

LoopStructural/interpolators/finite_difference_interpolator.py

@@ -203,8 +203,8 @@ def add_gradient_constraint(self, w=1.):
             A = np.einsum('ij,ijk->ik', strike_vector.T, T)
 
             B = np.zeros(points[inside, :].shape[0])
-            # self.add_constraints_to_least_squares(A * w, B, idc[inside, :])
-            A += np.einsum('ij,ijk->ik', dip_vector.T, T)
+            self.add_constraints_to_least_squares(A * w, B, idc[inside, :])
+            A = np.einsum('ij,ijk->ik', dip_vector.T, T)
             self.add_constraints_to_least_squares(A * w, B, idc[inside, :])
 
     def add_norm_constraint(self, w=1.):

LoopStructural/interpolators/piecewiselinear_interpolator.py

@@ -154,12 +154,12 @@ def add_gradient_ctr_pts(self, w=1.0):
             idc = gi[tetras]
             B = np.zeros(idc.shape[0])
             outside = ~np.any(idc == -1, axis=1)
-            # self.add_constraints_to_least_squares(A[outside, :] * w,
-            #                                       B[outside], idc[outside, :],
-            #                                       name = 'gradient')
-            A2 = np.einsum('ji,ijk->ik', dip_vector, element_gradients)
-            A2 *= vol[:, None]
-            A+=A2
+            self.add_constraints_to_least_squares(A[outside, :] * w,
+                                                  B[outside], idc[outside, :],
+                                                  name = 'gradient')
+            A = np.einsum('ji,ijk->ik', dip_vector, element_gradients)
+            A *= vol[:, None]
+            # A+=A2
             self.add_constraints_to_least_squares(A[outside, :] * w,
                                           B[outside], idc[outside, :],
                                                   name='gradient')
@@ -192,19 +192,17 @@ def add_norm_ctr_pts(self, w=1.0):
             # e, inside = self.support.elements_for_array(points[:, :3])
             # nodes = self.support.nodes[self.support.elements[e]]
             vol = np.zeros(element_gradients.shape[0])
-            vecs = vertices[inside, 1:, :] - vertices[inside, 0, None, :]
-            vol[inside] = np.abs(np.linalg.det(vecs))  # / 6
+            vecs = vertices[:, 1:, :] - vertices[:, 0, None, :]
+            vol = np.abs(np.linalg.det(vecs))  # / 6
             # 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]
 
             d_t = element_gradients
             d_t[inside,:,:] *= vol[inside, None, None]
-            # w*=10^11
-
             # add in the element gradient matrix into the inte
-            idc = np.tile(tetras[inside,:], (3, 1, 1))
+            idc = np.tile(tetras[:,:], (3, 1, 1))
             idc = idc.swapaxes(0,1)
             # idc = self.support.elements[e]
             gi = np.zeros(self.support.n_nodes).astype(int)
@@ -216,8 +214,8 @@ def add_norm_ctr_pts(self, w=1.0):
             outside = outside[:, 0]
             w /= 3
 
-            self.add_constraints_to_least_squares(d_t[inside,:,:][outside, :, :] * w,
-                                                  points[inside,:][outside, 3:6] * w *
+            self.add_constraints_to_least_squares(d_t[outside, :, :] * w,
+                                                  points[outside, 3:6] * w *
                                                   vol[outside, None],
                                                   idc[outside],
                                                   name='norm')

LoopStructural/modelling/core/geological_model.py

@@ -230,7 +230,10 @@ def set_model_data(self, data):
                 self.data.loc[mask, 'strike'], self.data.loc[mask, 'dip'])
             self.data.drop(['strike', 'dip'], axis=1, inplace=True)
         #     self.data.loc
-
+        # if 'nx' in self.data and 'ny' in self.data and 'nz' in self.data:
+        #     mask = np.all(~np.isnan(self.data.loc[:, ['nx', 'ny','nz']]),
+        #                   axis=1)
+        #     self.data.loc[mask,['nx', 'ny','nz']] /= self.scale_factor
     def extend_model_data(self, newdata):
         """
         Extends the data frame
@@ -314,16 +317,20 @@ def get_interpolator(self, interpolatortype='PLI', nelements=1e5,
         # add a buffer to the interpolation domain, this is necessary for
         # faults but also generally a good
         # idea to avoid boundary problems
+        # buffer = bb[1, :]
+        buffer = (bb[1,:]-bb[0,:])*buffer
         bb[0, :] -= buffer  # *(bb[1,:]-bb[0,:])
         bb[1, :] += buffer  # *(bb[1,:]-bb[0,:])
+        box_vol = (bb[1, 0]-bb[0, 0]) * (bb[1, 1]-bb[0, 1]) * (bb[1, 2]-bb[0, 2])
         if interpolatortype == "PLI":
             nelements /= 5
-            ele_vol = bb[1, 0] * bb[1, 1] * bb[1, 2] / nelements
+            ele_vol = box_vol / nelements
             # calculate the step vector of a regular cube
             step_vector = np.zeros(3)
             step_vector[:] = ele_vol ** (1. / 3.)
+            # step_vector /= np.array([1,1,2])
             # number of steps is the length of the box / step vector
-            nsteps = ((bb[1, :] - bb[0, :]) / step_vector).astype(int)
+            nsteps = np.ceil((bb[1, :] - bb[0, :]) / step_vector).astype(int)
             # create a structured grid using the origin and number of steps
             mesh_id = 'mesh_{}'.format(nelements)
             mesh = self.support.get(mesh_id,
@@ -338,12 +345,12 @@ def get_interpolator(self, interpolatortype='PLI', nelements=1e5,
 
         if interpolatortype == 'FDI':
             # find the volume of one element
-            ele_vol = bb[1, 0] * bb[1, 1] * bb[1, 2] / nelements
+            ele_vol = box_vol / nelements
             # calculate the step vector of a regular cube
             step_vector = np.zeros(3)
             step_vector[:] = ele_vol ** (1. / 3.)
             # number of steps is the length of the box / step vector
-            nsteps = ((bb[1, :] - bb[0, :]) / step_vector).astype(int)
+            nsteps = np.ceil((bb[1, :] - bb[0, :]) / step_vector).astype(int)
             # create a structured grid using the origin and number of steps
             grid_id = 'grid_{}'.format(nelements)
             grid = self.support.get(grid_id, StructuredGrid(origin=bb[0, :],
@@ -357,12 +364,12 @@ def get_interpolator(self, interpolatortype='PLI', nelements=1e5,
 
         if interpolatortype == "DFI":  # "fold" in kwargs:
             nelements /= 5
-            ele_vol = bb[1, 0] * bb[1, 1] * bb[1, 2] / nelements
+            ele_vol = box_vol / nelements
             # calculate the step vector of a regular cube
             step_vector = np.zeros(3)
             step_vector[:] = ele_vol ** (1. / 3.)
             # number of steps is the length of the box / step vector
-            nsteps = ((bb[1, :] - bb[0, :]) / step_vector).astype(int)
+            nsteps = np.ceil((bb[1, :] - bb[0, :]) / step_vector).astype(int)
             # create a structured grid using the origin and number of steps
             mesh = kwargs.get('mesh', TetMesh(origin=bb[0, :], nsteps=nsteps,
                                               step_vector=step_vector))
@@ -1020,7 +1027,7 @@ def voxet(self, nsteps=(50, 50, 25)):
         """
         return {'bounding_box': self.bounding_box, 'nsteps': nsteps}
 
-    def regular_grid(self, nsteps=(50, 50, 25), shuffle = True, rescale=True):
+    def regular_grid(self, nsteps=(50, 50, 25), shuffle = True, rescale=False):
         """
         Return a regular grid within the model bounding box
 

LoopStructural/modelling/features/geological_feature.py

@@ -186,12 +186,15 @@ def evaluate_gradient_misfit(self):
         """
         grad = self.interpolator.get_gradient_constraints()
         norm = self.interpolator.get_norm_constraints()
+
         dot = []
         if grad.shape[0] > 0:
+            grad /=np.linalg.norm(grad,axis=1)[:,None]
             model_grad = self.evaluate_gradient(grad[:,:3])
             dot.append(np.einsum('ij,ij->i',model_grad,grad[:,:3:6]).tolist())
 
         if norm.shape[0] > 0:
+            norm /=np.linalg.norm(norm,axis=1)[:,None]
             model_norm = self.evaluate_gradient(norm[:, :3])
             dot.append(np.einsum('ij,ij->i', model_norm, norm[:,:3:6]))
 

LoopStructural/utils/helper.py

@@ -251,13 +251,14 @@ def create_box(bounding_box, nsteps):
     return points, tri
 
 def get_vectors(normal):
-    normal /= np.linalg.norm(normal,axis=1)[:,None]
+    length =  np.linalg.norm(normal,axis=1)[:,None]
+    normal /= length#np.linalg.norm(normal,axis=1)[:,None]
     strikedip = normal_vector_to_strike_and_dip(normal)
     strike_vec = get_strike_vector(strikedip[:, 0])
     strike_vec /= np.linalg.norm(strike_vec,axis=0)[None,:]
     dip_vec = np.cross(strike_vec, normal, axisa=0, axisb=1).T  # (strikedip[:, 0], strikedip[:, 1])
     dip_vec /= np.linalg.norm(dip_vec,axis=0)[None,:]
-    return strike_vec, dip_vec
+    return strike_vec*length.T, dip_vec*length.T
 
 
 def get_strike_vector(strike):

LoopStructural/utils/map2loop.py

@@ -20,7 +20,7 @@ def process_map2loop(m2l_directory, flags={}):
     tangents = pd.read_csv(m2l_directory + '/tmp/raw_contacts.csv')
     groups = pd.read_csv(m2l_directory + '/tmp/all_sorts.csv', index_col=0)
     contact_orientations = pd.read_csv(m2l_directory + '/output/orientations.csv')
-    formation_thickness = pd.read_csv(m2l_directory + '/output/formation_thicknesses.csv')
+    # formation_thickness = pd.read_csv)
     contacts = pd.read_csv(m2l_directory + '/output/contacts4.csv')
     displacements = pd.read_csv(m2l_directory + '/output/fault_displacements3.csv')
     fault_orientations = pd.read_csv(m2l_directory + '/output/fault_orientations.csv')
@@ -43,8 +43,10 @@ def process_map2loop(m2l_directory, flags={}):
     except:
         for g in groups['group'].unique():
             supergroups[g] = g
-    supergroups.pop('\n')
-
+    try:
+        supergroups.pop('\n')
+    except KeyError:
+        pass
 
 
     bb = pd.read_csv(m2l_directory+'/tmp/bbox.csv')
@@ -56,19 +58,38 @@ def process_map2loop(m2l_directory, flags={}):
     tangents.drop(['angle', 'lsx', 'lsy'], inplace=True, axis=1)
 
     # convert azimuth and dip to gx, gy, gz
+
+
+    # calculate scalar field values
+    i = 0
+    thickness = {}
+    max_thickness = 0
+    with open(m2l_directory + '/output/formation_summary_thicknesses.csv') as file:
+        for l in file:
+            if i>1:
+                linesplit = l.split(',')
+                thickness[linesplit[0]] = float(linesplit[1])
+                # normalise the thicknesses
+                if float(linesplit[1]) > max_thickness:
+                    max_thickness=float(linesplit[1])
+    #             print(l.split(',')[1])
+            i+=1
+    # for k in thickness.keys():
+    #     thickness[k] /= max_thickness
+
     from LoopStructural.utils.helper import strike_dip_vector
     contact_orientations['strike'] = contact_orientations['azimuth'] - 90
-    contact_orientations['gx'] = np.nan
-    contact_orientations['gy'] = np.nan
-    contact_orientations['gz'] = np.nan
-    contact_orientations[['gx', 'gy', 'gz']] = strike_dip_vector(contact_orientations['strike'],
-                                                                 contact_orientations['dip'])
+    contact_orientations['nx'] = np.nan
+    contact_orientations['ny'] = np.nan
+    contact_orientations['nz'] = np.nan
+    contact_orientations[['nx', 'ny', 'nz']] = strike_dip_vector(contact_orientations['strike'],
+                                                                 contact_orientations['dip'])*max_thickness 
     contact_orientations.drop(['strike', 'dip', 'azimuth'], inplace=True, axis=1)
+    # with open(m2l_directory + '/output/formation_summary_thicknesses.csv') as file:
 
-    # calculate scalar field values
-    thickness = {}
-    for f in formation_thickness['formation'].unique():
-        thickness[f] = np.mean(formation_thickness[formation_thickness['formation'] == f]['thickness'])
+    # thickness = {}
+    # for f in formation_thickness['formation'].unique():
+    #     thickness[f] = formation_thickness[formation_thickness['formation'] == f]['thickness'])
 
     strat_val = {}
     stratigraphic_column = {}
@@ -83,10 +104,10 @@ def process_map2loop(m2l_directory, flags={}):
         for c in groups.loc[groups['group number'] == i, 'code']:
             strat_val[c] = np.nan
             if c in thickness:
-                stratigraphic_column[g][c] = {'min': val[g], 'max': val[g] + thickness[c], 'id': unit_id}
+                stratigraphic_column[g][c] = {'max': val[g], 'min': val[g] - thickness[c], 'id': unit_id}
                 unit_id += 1
                 strat_val[c] = val[g]
-                val[g] += thickness[c]
+                val[g] -= thickness[c]
     group_name = None
     for g, i in stratigraphic_column.items():
         if len(i) ==0:
@@ -159,7 +180,7 @@ def process_map2loop(m2l_directory, flags={}):
             'bounding_box':bb,
             'strat_va':strat_val}
 
-def build_model(m2l_data, skip_faults = False, unconformities=False, fault_params = None, foliation_params=None):
+def build_model(m2l_data, skip_faults = False, unconformities=False, fault_params = None, foliation_params=None, rescale = True):
     """[summary]
 
     [extended_summary]
@@ -191,7 +212,8 @@ def build_model(m2l_data, skip_faults = False, unconformities=False, fault_param
     boundary_points[1, 1] = m2l_data['bounding_box']['maxy']
     boundary_points[1, 2] = m2l_data['bounding_box']['upper']
 
-    model = GeologicalModel(boundary_points[0, :], boundary_points[1, :])
+    model = GeologicalModel(boundary_points[0, :], boundary_points[1, :], rescale=rescale)
+    # m2l_data['data']['val'] /= model.scale_factor
     model.set_model_data(m2l_data['data'])
     if not skip_faults:
         faults = []

LoopStructural/visualisation/model_visualisation.py

@@ -74,6 +74,8 @@ def deep_clean(self):
 
         [extended_summary]
         """
+        self.lv.clear()
+        self.lv.cleardata()
         pass
 
     def add_section(self, geological_feature=None, axis='x', value=None, **kwargs):
@@ -191,7 +193,7 @@ def add_isosurface(self, geological_feature, value = None, isovalue=None,
         # do isosurfacing of support using marching tetras/cubes
         x = np.linspace(self.bounding_box[0, 0], self.bounding_box[1, 0], self.nsteps[0])
         y = np.linspace(self.bounding_box[0, 1], self.bounding_box[1, 1], self.nsteps[1])
-        z = np.linspace(self.bounding_box[0, 2], self.bounding_box[1, 2], self.nsteps[2])
+        z = np.linspace(self.bounding_box[1, 2], self.bounding_box[0, 2], self.nsteps[2])
         xx, yy, zz = np.meshgrid(x, y, z, indexing='ij')
         points = np.array([xx.flatten(), yy.flatten(), zz.flatten()]).T
         val = geological_feature.evaluate_value(points)
@@ -521,7 +523,7 @@ def add_value_data(self, position, value, name, **kwargs):
         if "pointsize" not in kwargs:
             kwargs["pointsize"] = 4
         # set the colour map to diverge unless user decides otherwise
-        cmap = kwargs.get('cmap', "diverge")
+        cmap = kwargs.get('cmap', "spot")
         p = self.lv.points(name, **kwargs)
         p.vertices(position)
         p.values(value, "v")
@@ -576,6 +578,7 @@ def interactive(self, popout=False):
             self.lv.control.ObjectList()
             self.lv.interactive()
 
+
     def set_zscale(self,zscale):
         """ Set the vertical scale for lavavu
 
@@ -617,14 +620,17 @@ def save(self, fname, **kwargs):
         """
         self.lv.image(fname, **kwargs)
 
-    def display(self):
+    def display(self, fname=None, **kwargs):
         """
         Calls the lv object display function. Shows a static image of the viewer inline.
 
         Returns
         -------
 
         """
+        if fname:
+            self.lv.image(fname, **kwargs)
+
         self.lv.display()
 
     def image(self, name, **kwargs):
@@ -702,3 +708,25 @@ def rotate(self, r):
 
         """
         self.lv.rotate(r)
+
+    @property
+    def rotation(self):
+        """Accessor for the viewer rotation
+        Returns
+        -------
+        list
+            x,y,z rotations
+        """
+        return self.lv['xyzrotate']
+
+    @rotation.setter
+    def rotation(self,xyz):
+        """Set the rotation of the viewer
+
+        Parameters
+        ----------
+        xyz : list like
+            x y z rotations
+        """
+        self.lv.rotation(xyz)
+