Merge pull request #89 from Loop3D/unstructured_tetra.py

adding an unstructured tetrahedron class

Author: lachlangrose

LoopStructural/interpolators/__init__.py

@@ -8,6 +8,7 @@
 from LoopStructural.interpolators.piecewiselinear_interpolator import \
     PiecewiseLinearInterpolator
 from LoopStructural.interpolators.structured_tetra import TetMesh
+from LoopStructural.interpolators.unstructured_tetra import UnStructuredTetMesh
 from LoopStructural.interpolators.structured_grid import StructuredGrid
 from LoopStructural.interpolators.geological_interpolator import GeologicalInterpolator
 from LoopStructural.interpolators.discrete_interpolator import DiscreteInterpolator

LoopStructural/interpolators/base_linear_tetra.py

@@ -0,0 +1,24 @@
+# class LinearTetrahedron:
+#     def __init__(self):
+#         pass
+    
+#     def element_gradient(self, elements):
+#         ps = self.
+#         m = np.array(
+#             [[(ps[:, 1, 0] - ps[:, 0, 0]), (ps[:, 1, 1] - ps[:, 0, 1]),
+#               (ps[:, 1, 2] - ps[:, 0, 2])],
+#              [(ps[:, 2, 0] - ps[:, 0, 0]), (ps[:, 2, 1] - ps[:, 0, 1]),
+#               (ps[:, 2, 2] - ps[:, 0, 2])],
+#              [(ps[:, 3, 0] - ps[:, 0, 0]), (ps[:, 3, 1] - ps[:, 0, 1]),
+#               (ps[:, 3, 2] - ps[:, 0, 2])]])
+#         I = np.array(
+#             [[-1., 1., 0., 0.],
+#              [-1., 0., 1., 0.],
+#              [-1., 0., 0., 1.]])
+#         m = np.swapaxes(m, 0, 2)
+#         element_gradients = np.linalg.inv(m)
+
+#         element_gradients = element_gradients.swapaxes(1, 2)
+#         element_gradients = element_gradients @ I
+
+#         return element_gradients[elements,:,:]

LoopStructural/interpolators/base_structured_3d_support.py

@@ -22,7 +22,7 @@ def __init__(self,
         # we use property decorators to update these when different parts of
         # the geometry need to change
         # inisialise the private attributes
-        self._nsteps = np.array(nsteps)
+        self._nsteps = np.array(nsteps,dtype=int)
         self._step_vector = np.array(step_vector)
         self._origin = np.array(origin)  
         self.supporttype='Base'
@@ -181,4 +181,125 @@ def check_position(self, pos):
         if len(pos.shape) != 2:
             print("Position array needs to be a list of points or a point")
             return False
-        return pos
+        return pos
+
+    def global_cell_indicies(self, indexes):
+        """
+        Convert from cell indexes to global cell index
+
+        Parameters
+        ----------
+        indexes
+
+        Returns
+        -------
+
+        """
+        indexes = np.array(indexes).swapaxes(0, 2)
+        return indexes[:, :, 0] + self.nsteps_cells[None, None, 0] \
+               * indexes[:, :, 1] + self.nsteps_cells[None, None, 0] * \
+               self.nsteps_cells[None, None, 1] * indexes[:, :, 2]
+
+    def cell_corner_indexes(self, x_cell_index, y_cell_index, z_cell_index):
+        """
+        Returns the indexes of the corners of a cell given its location xi,
+        yi, zi
+
+        Parameters
+        ----------
+        x_cell_index
+        y_cell_index
+        z_cell_index
+
+        Returns
+        -------
+
+        """
+        xcorner = np.array([0, 1, 0, 0, 1, 0, 1, 1])
+        ycorner = np.array([0, 0, 1, 0, 0, 1, 1, 1])
+        zcorner = np.array([0, 0, 0, 1, 1, 1, 0, 1])
+        xcorners = x_cell_index[:, None] + xcorner[None, :]
+        ycorners = y_cell_index[:, None] + ycorner[None, :]
+        zcorners = z_cell_index[:, None] + zcorner[None, :]
+        return xcorners, ycorners, zcorners
+
+    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)
+        # if global index is not inside the support set to -1
+        globalidx[~inside] = -1
+        return globalidx, inside
+    
+    def node_indexes_to_position(self, xindex, yindex, zindex):
+
+        x = self.origin[0] + self.step_vector[0] * xindex
+        y = self.origin[1] + self.step_vector[1] * yindex
+        z = self.origin[2] + self.step_vector[2] * zindex
+
+        return x, y, z
+        
+    def global_index_to_cell_index(self, global_index):
+        """
+        Convert from global indexes to xi,yi,zi
+
+        Parameters
+        ----------
+        global_index
+
+        Returns
+        -------
+
+        """
+        # 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]
+        z_index = global_index // self.nsteps_cells[0, None] // \
+                  self.nsteps_cells[1, None]
+        return x_index, y_index, z_index
+
+
+    def global_index_to_node_index(self, global_index):
+        """
+        Convert from global indexes to xi,yi,zi
+
+        Parameters
+        ----------
+        global_index
+
+        Returns
+        -------
+
+        """
+        # 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[0, None]
+        y_index = global_index // self.nsteps[0, None] % \
+                  self.nsteps[1, None]
+        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):
+        """
+        Convert from node indexes to global node 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]

LoopStructural/interpolators/finite_difference_interpolator.py

@@ -266,7 +266,7 @@ def add_gradient_constraints(self, w=1.):
             # normalise constraint vector and scale element matrix by this
             norm = np.linalg.norm(points[:,3:6],axis=1)
             points[:,3:6]/=norm[:,None]
-            T/=norm[:,None,None]
+            T/=norm[inside,None,None]
             # calculate two orthogonal vectors to constraint (strike and dip vector)
             strike_vector, dip_vector = get_vectors(points[inside, 3:6])
             A = np.einsum('ij,ijk->ik', strike_vector.T, T)

LoopStructural/interpolators/piecewiselinear_interpolator.py

@@ -4,6 +4,7 @@
 import logging
 
 import numpy as np
+from LoopStructural.interpolators.cython.dsi_helper import cg, constant_norm, fold_cg
 
 from LoopStructural.interpolators.discrete_interpolator import \
     DiscreteInterpolator
@@ -17,23 +18,23 @@ class PiecewiseLinearInterpolator(DiscreteInterpolator):
     """
 
     """
-    def __init__(self, mesh):
+    def __init__(self, support):
         """
         Piecewise Linear Interpolator
         Approximates scalar field by finding coefficients to a piecewise linear
         equation on a tetrahedral mesh. Uses constant gradient regularisation.
 
         Parameters
         ----------
-        mesh - TetMesh
+        support - TetMesh
             interpolation support
         """
 
         self.shape = 'rectangular'
-        DiscreteInterpolator.__init__(self, mesh)
+        DiscreteInterpolator.__init__(self, support)
         # whether to assemble a rectangular matrix or a square matrix
         self.interpolator_type = 'PLI'
-        self.support = mesh
+        self.support = support
 
         self.interpolation_weights = {'cgw': 0.1, 'cpw': 1., 'npw': 1.,
                                       'gpw': 1., 'tpw': 1., 'ipw': 1.}
@@ -103,73 +104,58 @@ def add_constant_gradient(self, w= 0.1, direction_vector=None, direction_feature
 
         """
         if direction_feature is not None:
-            print('dir fe')
             direction_vector = direction_feature.evaluate_gradient(self.support.barycentre())
         if direction_vector is not None:
             if direction_vector.shape[0] == 1:
                 # if using a constant direction, tile array so it works for cg calc
                 direction_vector = np.tile(direction_vector,(self.support.barycentre().shape[0],1))
+        if direction_vector is not None:
+            logger.info("Running constant gradient")
+            elements_gradients = self.support.get_element_gradients(np.arange(self.ntetra))
+            if elements_gradients.shape[0] != direction_vector.shape[0]:
+                logger.error('Cannot add directional CG, vector field is not the correct length')
+                return
+            region = self.region.astype('int64')
+
+            neighbours = self.support.get_neighbours()
+            elements = self.support.get_elements()
+            idc, c, ncons = fold_cg(elements_gradients, direction_vector, neighbours.astype('int64'), elements.astype('int64'), self.nodes)
+
+            idc = np.array(idc[:ncons, :])
+            A = np.array(c[:ncons, :])
+            B = np.zeros(c.shape[0])
+            gi = np.zeros(self.support.n_nodes)
+            gi[:] = -1
+            gi[self.region] = np.arange(0, self.nx)
+            idc = gi[idc]
+            outside = ~np.any(idc == -1, axis=1)
 
-            
-        # iterate over all elements
-        A, idc, B = self.support.get_constant_gradient(region=self.region,direction=direction_vector)
-        A = np.array(A)
-        B = np.array(B)
-        idc = np.array(idc)
-
-        gi = np.zeros(self.support.n_nodes)
-        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, :] * w,
-                                              B[outside] * w, idc[outside, :],
-                                              name='regularisation')
-        return
-
-    def add_direction_constant_gradient(self, w= 0.1, direction_vector=None, direction_feature=None):
-        """
-        Add the constant gradient regularisation to the system where regularisation is projected
-        on a vector
-
-        Parameters
-        ----------
-        w (double) - weighting of the cg parameter
-        direction_vector
-        direction_feature
-
-        Returns
-        -------
-
-        """
-        if direction_feature:
-            print('dir fe')
-            direction_vector = direction_feature.evaluate_gradient(self.support.barycentre())
-        if direction_vector:
-            if direction_vector.shape[0] == 1:
-                # if using a constant direction, tile array so it works for cg calc
-                direction_vector = np.tile(direction_vector,(self.support.barycentre().shape[0],1))
-
-            
-        # iterate over all elements
-        A, idc, B = self.support.get_constant_gradient(region=self.region,direction=direction_vector)
-        A = np.array(A)
-        B = np.array(B)
-        idc = np.array(idc)
-
-        gi = np.zeros(self.support.n_nodes)
-        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, :] * w,
-                                              B[outside] * w, idc[outside, :],
-                                              name='directional regularisation')
-        return
+            # w/=A.shape[0]
+            self.add_constraints_to_least_squares(A[outside, :] * w,
+                                                B[outside] * w, idc[outside, :],
+                                                name='direction_regularisation')
+        else:
+            logger.info("Running constant gradient")
+            elements_gradients = self.support.get_element_gradients(np.arange(self.support.ntetra))
+            region = self.region.astype('int64')
+
+            neighbours = self.support.get_neighbours()
+            elements = self.support.get_elements()
+            idc, c, ncons = cg(elements_gradients, neighbours.astype('int64'), elements.astype('int64'), self.support.nodes,
+                               region.astype('int64'))
+
+            idc = np.array(idc[:ncons, :])
+            A = np.array(c[:ncons, :])
+            B = np.zeros(A.shape[0])
+            gi = np.zeros(self.support.n_nodes)
+            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, :] * w,
+                                                B[outside] * w, idc[outside, :],
+                                                name='regularisation')
 
 
     def add_gradient_constraints(self, w=1.0):
@@ -212,7 +198,7 @@ def add_gradient_constraints(self, w=1.0):
             gi = np.zeros(self.support.n_nodes).astype(int)
             gi[:] = -1
             gi[self.region] = np.arange(0, self.nx).astype(int)
-            w /= 3
+            # w /= 3
             idc = gi[tetras]
             B = np.zeros(idc.shape[0])
             outside = ~np.any(idc == -1, axis=1)
@@ -329,7 +315,6 @@ def add_interface_constraints(self, w=1.0):  # for now weight all value points t
         points = self.get_interface_constraints()
         if points.shape[0] > 1:
             vertices, c, tetras, inside = self.support.get_tetra_for_location(points[:,:3])
-            # print(points[inside,:].shape)
 
             gi = np.zeros(self.support.n_nodes)
             gi[:] = -1