Correction of cross-coupling stiffness terms between Bodies.

- Errors in the stiffness matrix from a Line going between two
  Bodies or a Body and a Point have been fixed.
- In the process, renamed Line.KAB to KBA for consistency.

moorpy/Catenary.py

@@ -50,7 +50,16 @@ def catenary(XF, ZF, L, EA, W, CB=0, alpha=0, HF0=0, VF0=0, Tol=0.000001, nNodes
     Returns
     -------
     : tuple
-        (end 1 horizontal tension, end 1 vertical tension, end 2 horizontal tension, end 2 vertical tension, info dictionary) [N] (positive up)
+        (end 1 horizontal tension, end 1 vertical tension, end 2 horizontal 
+        tension, end 2 vertical tension, info dictionary) [N] (positive up).
+        Info dictionary contains the following:
+        HF and VF - horizontal and vertical tension components of end B [N].
+        stiffnessA - 2D stiffness matrix for end A [N/m].
+        stiffnessB - 2D stiffness matrix for end B [N/m].
+        stiffnessBA - 2D stiffness matrix for force at B due to movement of A [N/m].
+        LBot - length of line section laying on the seabed [m].
+        ProfileType
+        Zextreme - extreme z coordinate of the line section (in direction of wet weight) [m].
     
     '''
 
@@ -187,7 +196,7 @@ def dV_dZ_s(z0, H):   # height off seabed to evaluate at (infinite if 0), horizo
         info["VF"] = 0.0
         info["stiffnessB"]  = np.array([[ dHF_dXF, 0.0], [0.0, dVF_dZF]])
         info["stiffnessA"]  = np.array([[ dHF_dXF, 0.0], [0.0, dVF_dZF]])
-        info["stiffnessAB"] = np.array([[-dHF_dXF, 0.0], [0.0, 0.0]])
+        info["stiffnessBA"] = np.array([[-dHF_dXF, 0.0], [0.0, 0.0]])
         info["LBot"] = L
         info['ProfileType'] = 0
         info['Zextreme'] = 0
@@ -235,7 +244,7 @@ def dV_dZ_s(z0, H):   # height off seabed to evaluate at (infinite if 0), horizo
             info["VF"] = VF
             info["stiffnessB"]  = np.array([[0.0, 0.0], [0.0, dVF_dZF]])
             info["stiffnessA"]  = np.array([[0.0, 0.0], [0.0, W]]) 
-            info["stiffnessAB"] = np.array([[0.0, 0.0], [0.0, 0.0]])
+            info["stiffnessBA"] = np.array([[0.0, 0.0], [0.0, 0.0]])
             info["LBot"] = L - LHanging
             info['ProfileType'] = 4
             info['Zextreme'] = 0
@@ -277,7 +286,7 @@ def dV_dZ_s(z0, H):   # height off seabed to evaluate at (infinite if 0), horizo
             info["VF"] = VF
             info["stiffnessB"]  = np.array([[0.0, 0.0], [0.0, W / np.sqrt(2.0*hB/EA_W + 1.0)]])
             info["stiffnessA"]  = np.array([[0.0, 0.0], [0.0, W / np.sqrt(2.0*hA/EA_W + 1.0)]])
-            info["stiffnessAB"] = np.array([[0.0, 0.0], [0.0, 0.0]])
+            info["stiffnessBA"] = np.array([[0.0, 0.0], [0.0, 0.0]])
             info["LBot"] = L - LHanging
             info['ProfileType'] = 5
             info['Zextreme'] = -hA
@@ -327,7 +336,7 @@ def dV_dZ_s(z0, H):   # height off seabed to evaluate at (infinite if 0), horizo
             info["VF"] = VF
             info["stiffnessB"]  = np.array([[ VF/ZF, 0.0], [0.0, 0.5*W]])
             info["stiffnessA"]  = np.array([[ VF/ZF, 0.0], [0.0, 0.5*W]]) 
-            info["stiffnessAB"] = np.array([[-VF/ZF, 0.0], [0.0,-0.5*W]])
+            info["stiffnessBA"] = np.array([[-VF/ZF, 0.0], [0.0,-0.5*W]])
             info["LBot"] = 0.0
             info['ProfileType'] = 6
             info['Zextreme'] = -hA
@@ -367,7 +376,7 @@ def dV_dZ_s(z0, H):   # height off seabed to evaluate at (infinite if 0), horizo
             info["VF"] = VF
             info["stiffnessB"]  = np.array([[ VF/ZF, 0.0], [0.0, EA/L]])
             info["stiffnessA"]  = np.array([[ VF/ZF, 0.0], [0.0, EA/L]]) 
-            info["stiffnessAB"] = np.array([[-VF/ZF, 0.0], [0.0,-EA/L]])
+            info["stiffnessBA"] = np.array([[-VF/ZF, 0.0], [0.0,-EA/L]])
             info["LBot"] = 0.0
             info['ProfileType'] = 6
             info['Zextreme'] = 0
@@ -530,7 +539,7 @@ def dV_dZ_s(z0, H):   # height off seabed to evaluate at (infinite if 0), horizo
 
                         info4['oths']["stiffnessA"] = np.array(K)
                         info4['oths']["stiffnessB"] = np.array(K)
-                        info4['oths']["stiffnessAB"] = np.array(-K)
+                        info4['oths']["stiffnessBA"] = np.array(-K)
 
                         info4['oths']['ProfileType'] = -1 # flag for the parabolic solution for the coordinates...
                         info4['oths']['error'] = False
@@ -560,7 +569,7 @@ def dV_dZ_s(z0, H):   # height off seabed to evaluate at (infinite if 0), horizo
 
                         info4['oths']["stiffnessA"] = np.array(K)
                         info4['oths']["stiffnessB"] = np.array(K)
-                        info4['oths']["stiffnessAB"] = np.array(-K)
+                        info4['oths']["stiffnessBA"] = np.array(-K)
 
                         info4['oths']['ProfileType'] = -2 # flag for the bilinear solution for the coordinates...
                         info4['oths']['error'] = False
@@ -734,7 +743,7 @@ def step_func_U(X, args, Y, info, Ytarget, err, tols, iter, maxIter):
             info['stiffnessB'] = np.array([[ dH_dX               ,  K2[0,1] *K1[0,0]*dxdH        ], 
                                            [K2[1,0] *K1[0,0]*dxdH,  K2[1,1] -K2[1,0]*dxdH*K2[0,1]]])
 
-            info['stiffnessAB']= np.array([[-K1[0,0] *K2[0,0]*dxdH,  K1[0,1] *K2[0,0]*dxdH ],    # this is the lower-left submatrix, A motions, B reaction forces (corrected/flipped sign of entry 0,1)
+            info['stiffnessBA']= np.array([[-K1[0,0] *K2[0,0]*dxdH,  K1[0,1] *K2[0,0]*dxdH ],    # this is the lower-left submatrix, A motions, B reaction forces (corrected/flipped sign of entry 0,1)
                                            [-K1[0,0] *dxdH*K2[1,0], -K1[0,1] *dxdH*K2[1,0] ]])
 
 
@@ -945,15 +954,15 @@ def step_func_U(X, args, Y, info, Ytarget, err, tols, iter, maxIter):
     # get A and AB stiffness matrices for catenary profiles here based on fairlead (B) stiffness matrix
     if ProfileType == 1:
         info['stiffnessA'] = np.array(info['stiffnessB'])
-        info['stiffnessAB'] = -info['stiffnessB']
+        info['stiffnessBA'] = -info['stiffnessB']
 
     elif ProfileType in [2,3,7]:
         if CB == 0.0:
             info['stiffnessA'] = np.array([[info['stiffnessB'][0,0], 0], [0, dV_dZ_s(Tol, HF)]])  # vertical term is very approximate 
-            info['stiffnessAB'] = np.array([[-info['stiffnessB'][0,0], 0], [0, 0]])  # note: A and AB stiffnesses for this case only valid if zero friction
+            info['stiffnessBA'] = np.array([[-info['stiffnessB'][0,0], 0], [0, 0]])  # note: A and AB stiffnesses for this case only valid if zero friction
         else:
             info['stiffnessA'] = np.ones([2,2]) * np.nan  # if friction, flag to ensure users don't use this
-            info['stiffnessAB'] = np.ones([2,2]) * np.nan  # if friction, flag to ensure users don't use this
+            info['stiffnessBA'] = np.ones([2,2]) * np.nan  # if friction, flag to ensure users don't use this
 
     # un-swap line ends if they've been previously swapped, and apply global sign convention 
     # (vertical force positive-up, horizontal force positive from A to B)
@@ -972,8 +981,8 @@ def step_func_U(X, args, Y, info, Ytarget, err, tols, iter, maxIter):
         info["stiffnessA"][1,0] = -info["stiffnessA"][1,0]
         info["stiffnessB"][0,1] = -info["stiffnessB"][0,1]
         info["stiffnessB"][1,0] = -info["stiffnessB"][1,0]
-        info["stiffnessAB"][0,1] = -info["stiffnessAB"][0,1]  # for cross coupling matrix could also maybe transpose? but should be symmetrical so no need
-        info["stiffnessAB"][1,0] = -info["stiffnessAB"][1,0]
+        info["stiffnessBA"][0,1] = -info["stiffnessBA"][0,1]  # for cross coupling matrix could also maybe transpose? but should be symmetrical so no need
+        info["stiffnessBA"][1,0] = -info["stiffnessBA"][1,0]
 
     else:
         FxA =  HA
@@ -995,8 +1004,8 @@ def step_func_U(X, args, Y, info, Ytarget, err, tols, iter, maxIter):
         info["stiffnessA"][1,0] = -info["stiffnessA"][1,0]
         info["stiffnessB"][0,1] = -info["stiffnessB"][0,1]
         info["stiffnessB"][1,0] = -info["stiffnessB"][1,0]
-        info["stiffnessAB"][0,1] = -info["stiffnessAB"][0,1]
-        info["stiffnessAB"][1,0] = -info["stiffnessAB"][1,0]
+        info["stiffnessBA"][0,1] = -info["stiffnessBA"][0,1]
+        info["stiffnessBA"][1,0] = -info["stiffnessBA"][1,0]
 
         # TODO <<< should add more info <<<
 

moorpy/MoorProps.py

@@ -503,7 +503,7 @@ def getAnchorCost(fx, fz, type="drag-embedment",soil_type='medium clay', method
     # mooring line sizing:  Tension limit for QS: 50% MBS.  Or FOS = 2
 
 
-    return anchorMatCost, anchorInstCost, anchorDecomCost   # [USD]
+    return anchorMatCost, anchorInstCost, anchorDecomCost, info   # [USD]
 
 
 def getAnchorProps(fx, fz, type="drag-embedment", display=0):

moorpy/line.py

@@ -63,7 +63,7 @@ def __init__(self, mooringSys, num, L, lineType, nSegs=100, cb=0, isRod=0, attac
         self.TB = 0
         self.KA = np.zeros([3,3])  # 3D stiffness matrix of end A [N/m]
         self.KB = np.zeros([3,3])  # 3D stiffness matrix of end B
-        self.KAB= np.zeros([3,3])  # 3D stiffness matrix of cross coupling between ends
+        self.KBA= np.zeros([3,3])  # 3D stiffness matrix of cross coupling between ends
 
         self.HF = 0           # fairlead horizontal force saved for next solve
         self.VF = 0           # fairlead vertical force saved for next solve
@@ -852,7 +852,7 @@ def staticSolve(self, reset=False, tol=0.0001, profiles=0):
         # save 3d stiffness matrix in global orientation for both line ends (3 DOF + 3 DOF)
         self.KA  = from2Dto3Drotated(info['stiffnessA'], -fBH, LH, R.T)  # reaction at A due to motion of A
         self.KB  = from2Dto3Drotated(info['stiffnessB'], -fBH, LH, R.T)  # reaction at B due to motion of B
-        self.KAB = from2Dto3Drotated(info['stiffnessAB'], fBH, LH, R.T)  # reaction at B due to motion of A
+        self.KBA = from2Dto3Drotated(info['stiffnessBA'], fBH, LH, R.T)  # reaction at B due to motion of A
 
         # may want to skip stiffness calcs when just getting profiles for plotting...
 
@@ -958,10 +958,12 @@ def getPosition(self, s):
 
 
     def getCost(self):
-        '''Fill in and returns a cost dictionary for this Line object.'''
+        '''Fill in a cost dictionary and return the total cost for this Line object.'''
         self.cost = {}  # clear any old cost numbers
         self.cost['material'] = self.type['cost']*self.L0
-        return cost
+        total_cost = sum(self.cost.values()) 
+        return total_cost
+
 
     def attachLine(self, lineID, endB):
         pass

moorpy/subsystem.py

@@ -195,7 +195,7 @@ def makeGeneric(self, lengths, types, points=[], shared=False):
         # if a points list is provided, apply any mass or other properties it contains?
 
 
-    def setEndPosition(self, r, endB):
+    def setEndPosition(self, r, endB, sink=False):
         '''Sets either end position of the subsystem in the global/system
         reference frame. This is included mainly to mimic the Line method.
 
@@ -205,13 +205,19 @@ def setEndPosition(self, r, endB):
             x,y,z coorindate position vector of the line end [m].
         endB : boolean
             An indicator of whether the r array is at the end or beginning of the line
+        sink : bool
+            If true, and if there is a subsystem, the z position will be on the seabed.
 
         Raises
         ------
         LineError
             If the given endB value is not a 1 or 0
         '''
 
+        if sink: # set z coordinate on seabed
+            z = getDepthFromBathymetry(r[0], r[1])
+            r = np.array([r[0], r[1], z])  # (making a copy of r to not overwrite it)
+
         # set end coordinates in global frame just like for a Line
         if endB == 1:
             self.rB = np.array(r, dtype=np.float_)
@@ -269,7 +275,7 @@ def staticSolve(self, reset=False, tol=0.01, profiles=0):
         R = np.eye(3)
         self.KA_L  = from2Dto3Drotated(K[:2,:2], -self.fB_L[0], LH, R.T)  # reaction at A due to motion of A
         self.KB_L  = from2Dto3Drotated(K[2:,2:], -self.fB_L[0], LH, R.T)  # reaction at B due to motion of B
-        self.KAB_L = from2Dto3Drotated(K[2:,:2], -self.fB_L[0], LH, R.T)  # reaction at B due to motion of A
+        self.KBA_L = from2Dto3Drotated(K[2:,:2], -self.fB_L[0], LH, R.T)  # reaction at B due to motion of A
 
         self.TA = np.linalg.norm(self.fA_L)  # tensions [N]
         self.TB = np.linalg.norm(self.fB_L)
@@ -280,7 +286,7 @@ def staticSolve(self, reset=False, tol=0.01, profiles=0):
 
         self.KA  = np.matmul(np.matmul(self.R, self.KA_L ), self.R.T)
         self.KB  = np.matmul(np.matmul(self.R, self.KB_L ), self.R.T)
-        self.KAB = np.matmul(np.matmul(self.R, self.KAB_L), self.R.T)
+        self.KBA = np.matmul(np.matmul(self.R, self.KBA_L), self.R.T)
 
         #self.LBot = info["LBot"]  <<< this should be calculated considering all lines