From 2f5906596111999fad29b047eba70beea0f083a5 Mon Sep 17 00:00:00 2001 From: Lucas-Carmo Date: Fri, 24 May 2024 17:42:23 -0600 Subject: [PATCH] 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 | 22 +++++++++++----------- moorpy/system.py | 35 +---------------------------------- 2 files changed, 12 insertions(+), 45 deletions(-) diff --git a/moorpy/point.py b/moorpy/point.py index 637ca10..086be00 100644 --- a/moorpy/point.py +++ b/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) + # + matrix = (matrix + matrix.T)/2 # To make sure it is symmetric (it is not due to numerical errors) + 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) + # matrix_out = (matrix_out + matrix_out.T)/2 # To make sure it is symmetric (it is not due to numerical errors) + 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): diff --git a/moorpy/system.py b/moorpy/system.py index 90570f5..e831c02 100644 --- a/moorpy/system.py +++ b/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