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

Author: lachlangrose

LoopStructural/interpolators/cython/dsi_helper.pyx

@@ -93,6 +93,91 @@ def cg(double [:,:,:] EG, long long [:,:] neighbours, long long [:,:] elements,d
                     c[ncons,position_to_write] -= norm[i]*e2[i][itr_right]*area
             ncons+=1
     return idc, c, ncons
+def constant_norm(double [:,:,:] EG, long long [:,:] neighbours, long long [:,:] elements,double [:,:] nodes, long long [:] region):
+    cdef int Nc, Na, i,Ns, j, ne, ncons, e, n, neigh
+    Nc = 5 #numer of constraints shared nodes + independent
+    Na = 4 #number of nodes
+    Ns = Na -1
+    ne = len(neighbours)
+    ncons = 0
+    cdef int [:] flag = np.zeros(ne,dtype=np.int32)
+    cdef double [:,:] c = np.zeros((len(neighbours)*4,Nc))
+    cdef long long [:,:] idc = np.zeros((ne*4,5),dtype=np.int64)
+    cdef long long [3] common
+    cdef double [:] norm = np.zeros((3))
+    cdef double [:,:] shared_pts = np.zeros((3,3))
+    cdef double [:] v1 = np.zeros(3)
+    cdef double [:] v2 = np.zeros(3)
+    cdef double [:,:] e1
+    cdef double [:,:] e2
+    cdef double area = 0
+    cdef long long [:] idl  = np.zeros(4,dtype=np.int64)
+    cdef long long [:] idr = np.zeros(4,dtype=np.int64)
+    for e in range(ne):
+        idl = elements[e,:]
+        e1 = EG[e,:,:]
+        flag[e] = 1
+        # if not in region then skip this tetra
+        if region[idl[0]] == 0 or region[idl[1]] == 0 or region[idl[2]] == 0 or region[idl[3]] == 0:
+            continue
+        for n in range(4):
+            neigh = neighbours[e,n]
+            idr = elements[neigh,:]
+            if neigh < 0:
+                continue
+            if flag[neigh]== 1:
+                continue
+
+            # if not in region then skip this tetra
+            if region[idr[0]] == 0 or region[idr[1]] == 0 or region[idr[2]] == 0 or region[idr[3]] == 0:
+                continue
+
+            e2 = EG[neigh,:,:]
+
+            for i in range(Nc):
+                idc[ncons,i] = -1
+
+            i = 0
+            for itr_right in range(Na):
+                for itr_left in range(Na):
+                    if idl[itr_left] == idr[itr_right]:
+                        common[i] = idl[itr_left]
+                        i+=1
+            for j in range(3):
+                for k in range(3):
+                    shared_pts[j][k] = nodes[common[j]][k]#common
+            for i in range(3):
+                v1[i] = shared_pts[0,i] - shared_pts[1,i]
+                v2[i] = shared_pts[2,i]-shared_pts[1,i]
+            norm[0] = v2[2]*v1[1] - v1[2]*v2[1]
+            norm[1] = v1[2]*v2[0] - v1[0]*v2[2]
+            norm[2] = v1[0]*v2[1] - v1[1]*v2[0]
+
+            # we want to weight the cg by the area of the shared face
+            # area of triangle is half area of parallelogram
+            # https://math.stackexchange.com/questions/128991/how-to-calculate-the-area-of-a-3d-triangle
+            area = 0.5*np.linalg.norm(norm)
+            for itr_left in range(Na):
+                idc[ncons,itr_left] = idl[itr_left]
+                c[ncons,itr_left] = np.sqrt(e1[0][itr_left]*e1[0][itr_left]+e1[1][itr_left]*e1[1][itr_left]+e1[2][itr_left]*e1[2][itr_left])
+              
+            next_available_position = Na
+            for itr_right in range(Na):
+                common_index = -1
+                for itr_left in range(Na):
+                    if idc[ncons,itr_left] == idr[itr_right]:
+                        common_index = itr_left
+
+                position_to_write = 0
+                if common_index != -1:
+                    position_to_write = common_index
+                else:
+                    position_to_write = 4#next_available_position
+                    next_available_position+=1
+                idc[ncons,position_to_write] = idr[itr_right]
+                c[ncons,position_to_write] -= np.sqrt(e2[0][itr_right]*e2[0][itr_right]+e2[1][itr_right]*e2[1][itr_right]+e2[2][itr_right]*e2[2][itr_right]) 
+            ncons+=1
+    return idc, c, ncons
 def fold_cg(double [:,:,:] EG, double [:,:] X, long long [:,:] neighbours, long long [:,:] elements,double [:,:] nodes):
     cdef int Nc, Na, i,Ns, j, ne, ncons, e, n, neigh
     Nc = 5 #numer of constraints shared nodes + independent

LoopStructural/interpolators/discrete_fold_interpolator.py

@@ -82,9 +82,9 @@ def update_fold(self, fold):
         logger.error('updating fold, this should be done by accessing the fold attribute')
         self.fold = fold
 
-    def add_fold_constraints(self, fold_orientation=10., fold_axis_w=10., fold_regularisation=.1,
+    def add_fold_constraints(self, fold_orientation=10., fold_axis_w=10., fold_regularisation=[.1,0.01,0.01],
                              fold_normalisation=1.,
-                             fold_norm=1.):
+                             fold_norm=1., step=2):
         """
 
         Parameters
@@ -93,12 +93,15 @@ def add_fold_constraints(self, fold_orientation=10., fold_axis_w=10., fold_regul
             weight for the fold direction/orientation in the least squares system
         fold_axis_w : double
             weight for the fold axis in the least squares system
-        fold_regularisation : double
+        fold_regularisation : list
             weight for the fold regularisation in the least squares system
         fold_normalisation : double
             weight for the fold norm constraint in the least squares system
         fold_norm
             length of the interpolation norm in the least squares system
+        step: int
+            array step for adding constraints
+
 
         Returns
         -------
@@ -115,57 +118,84 @@ def add_fold_constraints(self, fold_orientation=10., fold_axis_w=10., fold_regul
         # calculate the fold geometry for the elements barycentre
         deformed_orientation, fold_axis, dgz = \
             self.fold.get_deformed_orientation(self.support.barycentre())
-
+        element_idx = np.arange(self.support.n_elements)
+        np.random.shuffle(element_idx)
         # calculate element volume for weighting
         vecs = nodes[:, 1:, :] - nodes[:, 0, None, :]
         vol = np.abs(np.linalg.det(vecs)) / 6
         if fold_orientation is not None:
             """
             dot product between vector in deformed ori plane = 0
             """
+            np.random.shuffle(element_idx)
+
             logger.info("Adding fold orientation constraint to %s w = %f"%(self.propertyname, fold_orientation))
-            A = np.einsum('ij,ijk->ik', deformed_orientation, eg)
-            A *= vol[:, None]
+            A = np.einsum('ij,ijk->ik', deformed_orientation[element_idx[::step],:], eg[element_idx[::step],:,:])
+            A *= vol[element_idx[::step], None]
             A *= fold_orientation
-            B = np.zeros(self.support.n_elements)
-            idc = self.support.get_elements()
+            B = np.zeros(A.shape[0])
+            idc = self.support.get_elements()[element_idx[::step],:]
             self.add_constraints_to_least_squares(A, B, idc)
 
         if fold_axis_w is not None:
             """
             dot product between axis and gradient should be 0
             """
+            np.random.shuffle(element_idx)
+
             logger.info("Adding fold axis constraint to %s w = %f"%(self.propertyname,fold_axis_w))
-            A = np.einsum('ij,ijk->ik', fold_axis, eg)
-            A *= vol[:, None]
+            A = np.einsum('ij,ijk->ik', fold_axis[element_idx[::step],:], eg[element_idx[::step],:,:])
+            A *= vol[element_idx[::step], None]
             A *= fold_axis_w
-            B = np.zeros(self.support.n_elements).tolist()
-            self.add_constraints_to_least_squares(A, B, self.support.get_elements())
+            B = np.zeros(A.shape[0]).tolist()
+            idc = self.support.get_elements()[element_idx[::step],:]
+
+            self.add_constraints_to_least_squares(A, B, idc)
 
         if fold_normalisation is not None:
             """
             specify scalar norm in X direction
             """
+            np.random.shuffle(element_idx)
+
             logger.info("Adding fold normalisation constraint to %s w = %f"%(self.propertyname,fold_normalisation))
-            A = np.einsum('ij,ijk->ik', dgz, eg)
-            A *= vol[:, None]
+            A = np.einsum('ij,ijk->ik', dgz[element_idx[::step],:], eg[element_idx[::step],:,:])
+            A *= vol[element_idx[::step], None]
             A *= fold_normalisation
-            B = np.ones(self.support.n_elements)
+            B = np.ones(A.shape[0])
 
             if fold_norm is not None:
                 B[:] = fold_norm
             B *= fold_normalisation
-            B *= vol
-            self.add_constraints_to_least_squares(A, B, self.support.get_elements())
+            B *= vol[element_idx[::step]]
+            idc = self.support.get_elements()[element_idx[::step],:]
+
+            self.add_constraints_to_least_squares(A, B, idc)
 
         if fold_regularisation is not None:
             """
             fold constant gradient  
             """
-            logger.info("Adding fold regularisation constraint to %s w = %f"%(self.propertyname,fold_regularisation))
+            logger.info("Adding fold regularisation constraint to {} w = {} {} {}".format(self.propertyname,
+            fold_regularisation[0],fold_regularisation[1],fold_regularisation[1]))
+
             idc, c, ncons = fold_cg(eg, dgz, self.support.get_neighbours(), self.support.get_elements(), self.support.nodes)
             A = np.array(c[:ncons, :])
-            A *= fold_regularisation
+            A *= fold_regularisation[0]
+            B = np.zeros(A.shape[0])
+            idc = np.array(idc[:ncons, :])
+            self.add_constraints_to_least_squares(A, B, idc)
+
+            idc, c, ncons = fold_cg(eg, deformed_orientation, self.support.get_neighbours(), self.support.get_elements(), self.support.nodes)
+            A = np.array(c[:ncons, :])
+            A *= fold_regularisation[1]
+            B = np.zeros(A.shape[0])
+            idc = np.array(idc[:ncons, :])
+            self.add_constraints_to_least_squares(A, B, idc)
+
+            idc, c, ncons = fold_cg(eg, fold_axis, self.support.get_neighbours(), self.support.get_elements(), self.support.nodes)
+            A = np.array(c[:ncons, :])
+            A *= fold_regularisation[2]
             B = np.zeros(A.shape[0])
             idc = np.array(idc[:ncons, :])
             self.add_constraints_to_least_squares(A, B, idc)

LoopStructural/interpolators/discrete_interpolator.py

@@ -147,12 +147,22 @@ def add_constraints_to_least_squares(self, A, B, idc, name='undefined'):
         A = np.array(A)
         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: {}".format(name))
-            return
+
         nr = A.shape[0]
+        if A.shape != idc.shape:
+            return
+        
         if len(A.shape) > 2:
             nr = A.shape[0] * A.shape[1]
+            A = A.reshape((A.shape[0]*A.shape[1],A.shape[2]))
+            idc = idc.reshape((idc.shape[0]*idc.shape[1],idc.shape[2]))
+        # going to assume if any are nan they are all nan
+        mask = np.any(np.isnan(A),axis=1)
+        A[mask,:] = 0
+        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: {}".format(name))
+            # return
+        
         rows = np.arange(0, nr).astype(int)
         rows += self.c_
         constraint_ids = rows.copy()

LoopStructural/interpolators/piecewiselinear_interpolator.py

@@ -72,7 +72,10 @@ def _setup_interpolator(self, **kwargs):
             self.interpolation_weights[key] = kwargs[key]
         if self.interpolation_weights['cgw'] > 0.:
             self.up_to_date = False
-            self.add_constant_gradient(self.interpolation_weights['cgw'])
+            self.add_constant_gradient(self.interpolation_weights['cgw'],
+            direction_feature=kwargs.get('direction_feature',None),
+            direction_vector=kwargs.get('direction_vector',None)
+            )
             logger.info("Using constant gradient regularisation w = %f"
                         %self.interpolation_weights['cgw'])
         logger.info("Added %i gradient constraints, %i normal constraints,"
@@ -84,7 +87,10 @@ def _setup_interpolator(self, **kwargs):
         self.add_ctr_pts(self.interpolation_weights['cpw'])
         self.add_tangent_ctr_pts(self.interpolation_weights['tpw'])
         self.add_interface_ctr_pts(self.interpolation_weights['ipw'])
-    def add_constant_gradient(self, w=0.1):
+        if 'constant_norm' in kwargs:
+            self.add_constant_norm(w=kwargs['constant_norm'])
+    
+    def add_constant_gradient(self, w= 0.1, direction_vector=None, direction_feature=None):
         """
         Add the constant gradient regularisation to the system
 
@@ -96,8 +102,59 @@ def add_constant_gradient(self, w=0.1):
         -------
 
         """
+        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))
+
+            
+        # 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)
+        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)
@@ -114,6 +171,37 @@ def add_constant_gradient(self, w=0.1):
                                               name='regularisation')
         return
 
+
+    def add_constant_norm(self, w=0.1):
+        """
+        Add the constant gradient regularisation to the system
+
+        Parameters
+        ----------
+        w (double) - weighting of the cg parameter
+
+        Returns
+        -------
+
+        """
+        # iterate over all elements
+        A, idc, B = self.support.get_constant_norm(region=self.region)
+        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='norm_regularisation')
+        return
+
     def add_gradient_ctr_pts(self, w=1.0):
         """
         Adds gradient constraints to the least squares system with a weight
@@ -325,6 +413,7 @@ def add_gradient_orthogonal_constraint(self, points, vector, w=1.0,
             vertices, element_gradients, tetras, inside = self.support.get_tetra_gradient_for_location(points[:,:3])
             #e, inside = self.support.elements_for_array(points[:, :3])
             #nodes = self.support.nodes[self.support.elements[e]]
+            vector /= np.linalg.norm(vector,axis=1)[:,None]
             vecs = vertices[:, 1:, :] - vertices[:, 0, None, :]
             vol = np.abs(np.linalg.det(vecs))  # / 6
             # d_t = self.support.get_elements_gradients(e)
@@ -340,7 +429,7 @@ def add_gradient_orthogonal_constraint(self, points, vector, w=1.0,
             gi[self.region] = np.arange(0, self.nx).astype(int)
             w /= 3
             idc = gi[tetras]
-            B = np.zeros(idc.shape[0])
+            B = np.zeros(idc.shape[0])+B
             outside = ~np.any(idc == -1, axis=1)
             self.add_constraints_to_least_squares(A[outside, :] * w,
                                                   B[outside], idc[outside, :])