Fix double counting of an effect + deal with matrices w/ null rows/cols

- In system.getSystemDynamicMatrices, the effect of other attached points (the loop on point2) is already lumped in point.getDynamicMatrices
-  Replaced np.linalg.inv with np.linalg.pinv (which computes the pseudo-inverse of the matrix) because we can have some cases where rows and columns of the matrix are zero (e.g., when adopting Ca=0 in a given direction)
- The procedure to lump the dynamic matrices in a point involves inverting the matrices twice. If the input matrix is not exactly symmetrical, this can greatly amplify the spurious asymmetry. For this reason, we make the matrices computed with the lumped mass symmetrical.
- For some reason that I don't know yet, the procedure does not work if we use the whole matrix. We need to remove the rows and columns associated with the dofs of the node opposite to the node where we want to lump the matrices.

moorpy/point.py

@@ -341,7 +341,7 @@ def getStiffnessA(self, lines_only=False, xyz=False):
         else:                       # otherwise only return rows/columns of active DOFs
             return K[:,self.DOFs][self.DOFs,:]
 
-    def getDynamicMatrices(self, omegas, S_zeta, r_dynamic_init, depth, kbot, cbot, xyz=False):
+    def getDynamicMatrices(self, omegas, S_zeta, r_dynamic_init, depth, kbot, cbot):
         '''Gets inertia, added mass, damping, and stiffness matrices of Point due only to mooring lines (with other objects fixed)
         using a lumped mass model.
 
@@ -355,15 +355,19 @@ def getDynamicMatrices(self, omegas, S_zeta, r_dynamic_init, depth, kbot, cbot,
         #print("Getting Point "+str(self.number)+" analytic stiffness matrix...")
 
         def lump_matrix(matrix, endB):
-            if endB == 1:
-                matrix_woAnc = matrix[3:, 3:]
-                matrix_inv = np.linalg.inv(matrix_woAnc)
+            # The matrices should be symmetrical, but they can be slightly off due to numerical errors.
+            # Because we are going to invert them twice, we force them to be symmetrical to avoid amplifying the errors. 
+            matrix = (matrix + matrix.T)/2
+            if endB == 1:            
+                matrix_woAnc = matrix[3:, 3:] # For some reason that I don't know yet, we need to remove the elements in the opposite side of the matrix. Otherwise it doesn't work
+                matrix_inv = np.linalg.pinv(matrix_woAnc)
                 matrix_inv_coupled = matrix_inv[-3:, -3:]
             else:
                 matrix_woAnc = matrix[:-3, :-3]
-                matrix_inv = np.linalg.inv(matrix_woAnc)
+                matrix_inv = np.linalg.pinv(matrix_woAnc)
                 matrix_inv_coupled = matrix_inv[0:3, 0:3]
-            return np.linalg.inv(matrix_inv_coupled)
+            matrix_out = np.linalg.inv(matrix_inv_coupled)
+            return matrix_out 
 
         M = np.zeros([3,3])
         A = np.zeros([3,3])
@@ -382,11 +386,7 @@ def lump_matrix(matrix, endB):
             B += lump_matrix(B_all, endB)
             K += lump_matrix(K_all, endB)
 
-        if sum(np.isnan(K).ravel()) > 0: breakpoint()
-        if xyz:                     # if asked to output all DOFs, do it
-            return M, A, B, K
-        else:                       # otherwise only return rows/columns of active DOFs            
-            return M[:,self.DOFs][self.DOFs,:], A[:,self.DOFs][self.DOFs,:], B[:,self.DOFs][self.DOFs,:], K[:,self.DOFs][self.DOFs,:]
+        return M, A, B, K
 
 
     def getCost(self):

moorpy/system.py

@@ -3065,40 +3065,7 @@ def getSystemDynamicMatrices(self, omegas, S_zeta, r_dynamic_init, depth, kbot,
                 A[i:i+n,i:i+n] += A1
                 B[i:i+n,i:i+n] += B1 
                 K[i:i+n,i:i+n] += K1
-
-                # # go through attached lines and add cross-coupling terms
-                for lineID in point.attached:
-
-                    j = i + n
-
-                    # go through movable points to see if one is attached
-                    for point2 in self.pointList[self.pointList.index(point)+1: ]:
-                        if point2.type in d:
-                            if lineID in point2.attached:  # if this point is at the other end of the line
-
-                                # get cross-coupling stiffness of line: force on end attached to point1 due to motion of other end
-                                MB, AB, BB, KB = point2.getDynamicMatrices(omegas, S_zeta, r_dynamic_init, depth, kbot, cbot)
-                                if point.attachedEndB != 1:
-                                    MB = MB.T
-                                    AB = AB.T
-                                    BB = BB.T
-                                    KB = KB.T
-
-                                # Trim stiffness matrix to only use the enabled DOFs of each point
-                                MB = MB[point.DOFs,:][:,point2.DOFs]
-                                AB = AB[point.DOFs,:][:,point2.DOFs]
-                                BB = BB[point.DOFs,:][:,point2.DOFs]
-                                KB = KB[point.DOFs,:][:,point2.DOFs]                                
-
-                                M[i:i+n          , j:j+point2.nDOF] += MB    # force on P1 due to movement of P2
-                                M[j:j+point2.nDOF, i:i+n          ] += MB.T  # mirror (f on P2 due to x of P1)
-                                A[i:i+n          , j:j+point2.nDOF] += AB  
-                                A[j:j+point2.nDOF, i:i+n          ] += AB.T
-                                B[i:i+n          , j:j+point2.nDOF] += BB  
-                                B[j:j+point2.nDOF, i:i+n          ] += BB.T
-                                K[i:i+n          , j:j+point2.nDOF] += KB  
-                                K[j:j+point2.nDOF, i:i+n          ] += KB.T                                
-                            j += point2.nDOF                  # if this point has DOFs we're considering, then count them                                
+
                 i += n
 
         return M, A, B, K