First working draft of Subsystem

- Subsystem class now works. It can be used independently like a
  System or act as a Line within a System as part of its lineList.
- System.addSubsystem can integrate an already created Subsystem
  in the System's lineList.
- More to do in terms of easy initialization of Subsystems.
- Added a few additional functions as a start for tracking mooring
  loads and costs, as well as eventually anchor types.

moorpy/Catenary.py

@@ -86,7 +86,7 @@ def catenary(XF, ZF, L, EA, W, CB=0, alpha=0, HF0=0, VF0=0, Tol=0.000001, nNodes
     if W < 0:
         W = -W
         ZF = -ZF
-        CB = -10000.   # <<< TODO: set this to the distance to sea surface <<<
+        CB = 0 #-10000.   # <<< TODO: could set hA, hB to distances to sea surface <<<
         flipFlag = True
     else:
         flipFlag = False
@@ -167,6 +167,8 @@ def dV_dZ_s(z0, H):   # height off seabed to evaluate at (infinite if 0), horizo
             HF = (XF/L -1.0)*EA + 0.5*CB*W*L
             HA = np.max([0.0, HF - CB*W*L])
         else:                                   # case 3: seabed friction and zero anchor tension
+            if EA*CB*W*(XF-L) < 0: # something went wrong
+                breakpoint() 
             HF = np.sqrt(2*EA*CB*W*(XF-L))
             HA = 0.0
 
@@ -213,7 +215,7 @@ def dV_dZ_s(z0, H):   # height off seabed to evaluate at (infinite if 0), horizo
                     else:                           # the tension is nonzero
                         Xs[I] = s[I] + CB*W/EA*(s[I] - xB)**2
                         Te[I] = HF - CB*W*(L-s[I])
-
+        
 
     # ProfileType 4 case - fully slack
     elif (W > 0.0) and (L >= XF/np.cos(np.radians(alpha)) + LHanging):
@@ -1383,7 +1385,9 @@ def step_func_cat(X, args, Y, info, Ytarget, err, tols, iter, maxIter):
     print(f" F_lateral / ((strain+tol)EA) is {F_lateral/(((d+Tol)/L-1)*EA):6.2e} !!!!!!")
     print(f" F_lateral / ((strain-tol)EA) is {F_lateral/(((d-Tol)/L-1)*EA):6.2e} !!!!!!")
 
-    (fAH1, fAV1, fBH1, fBV1, info1) = catenary(XF, ZF, L, EA, W, CB=-20, Tol=Tol, MaxIter=40, plots=2)
+    #(fAH1, fAV1, fBH1, fBV1, info1) = catenary(XF, ZF, L, EA, W, CB=-20, Tol=Tol, MaxIter=40, plots=2)
+
+    (fAH1, fAV1, fBH1, fBV1, info1) = catenary(400.00297786197154, 0.0, 400.0, 700000000.0, -221.18695627569764, CB=0.0, alpha=-0.0, HF0=5211.258450277256, VF0=0.0, Tol=2e-05, MaxIter=100, plots=1)
 
     print((fAH1, fAV1, fBH1, fBV1))
 

moorpy/line.py

@@ -47,6 +47,7 @@ def __init__(self, mooringSys, num, L, lineType, nSegs=100, cb=0, isRod=0, attac
         self.L = L              # line unstretched length (may be modified if using nonlinear elasticity) [m]
         self.L0 = L             # line reference unstretched length [m]
         self.type = lineType    # dictionary of a System.lineTypes entry
+        self.cost = {}          # empty dictionary to contain cost information
 
         self.EA = self.type['EA']  # use the default stiffness value for now (may be modified if using nonlinear elasticity) [N]
 
@@ -953,6 +954,12 @@ def getPosition(self, s):
         return np.vstack([ Xs, Ys, Zs])
 
 
+    def getCost(self):
+        '''Fill in and returns a cost dictionary for this Line object.'''
+        self.cost = {}  # clear any old cost numbers
+        self.cost['material'] = self.type['cost']*self.L0
+        return cost
+
     def attachLine(self, lineID, endB):
         pass
 
@@ -1004,7 +1011,7 @@ def from2Dto3Drotated(K2D, F, L, R):
     L : float
         Line horizontal distance end-to-end [m]
     R : 3x3 matrix
-        Rotation matrix from global frame to plane to the local
+        Rotation matrix from global frame to the local
         X-Z plane of the line
         
     Returns

moorpy/point.py

@@ -2,8 +2,6 @@
 
 import numpy as np
 
-
-
 class Point():
     '''A class for any object in the mooring system that can be described by three translational coorindates'''
 
@@ -50,7 +48,10 @@ def __init__(self, mooringSys, num, type, r, m=0, v=0, fExt=np.zeros(3), DOFs=[0
         self.number = num
         self.type = type                # 1: fixed/attached to something, 0 free to move, or -1 coupled externally
         self.r = np.array(r, dtype=np.float_)
-
+        self.entity = {type:''}         # dict for entity (e.g. anchor) info
+        self.cost = {}                  # empty dictionary to contain cost info
+        self.loads = {}                 # empty dictionary to contain load info
+
         self.m  = float(m)
         self.v  = float(v)
         self.CdA= float(CdA)
@@ -332,13 +333,38 @@ def getStiffnessA(self, lines_only=False, xyz=False):
             # if on seabed, apply a large stiffness to help out system equilibrium solve (if it's transitioning off, keep it a small step to start with)    
             if self.r[2] == -self.sys.depth:
                 K[2,2] += 1.0e12
-
+        if sum(np.isnan(K).ravel()) > 0: breakpoint()
         if xyz:                     # if asked to output all DOFs, do it
             return K
         else:                       # otherwise only return rows/columns of active DOFs
             return K[:,self.DOFs][self.DOFs,:]
 
+
+    def getCost(self):
+        '''Fill in and returns a cost dictionary for this Point object.
+        So far it only applies for if the point is an anchor.'''
+
+        from moorpy.MoorProps import getAnchorCost
+
+        self.cost = {'material':0}  # clear any old cost numbers and start with 0
+
+        # figure out if it should be an anchor if it isn't already defined
+        if self.entity['type'] == '':
+            depth, _ = self.sys.getDepthFromBathymetry(self.r[0], self.r[1]) 
+            if self.r[3] == depth and self.type==1:  # if it's fixed on the seabed
+                self.entity['type'] = 'anchor'       # assume it's an anchor
+                if self.FA[2] == 0:
+                    self.entity['anchor_type'] = 'drag-embedment'
+                else:
+                    self.entity['anchor_type'] = 'suction'
+
+        # calculate costs if it's an anchor (using simple model)
+        if self.entity['type'] == 'anchor':
+            self.cost['material'] = getAnchorCost(self.loads['fx_max'], 
+                                                  self.loads['fz_max'],
+                                             type=self.entity['anchor_type'])
 
+        return cost