New dymamic functions from Serag now working with Subsystem:

- Edited new Line dynamic methods to streamline passing of arguments
  and to only load the new functions on demand.
- Added Subsystem.getDynamicMatrices method - copied closely from the
  CompositeLine method. This allows the new methods to work with
  Subsystem (all other methods are already inherited from Line).
- Added Subsystem.makeFromSystem method - adapted from CompositeLine
  initialization to allow similar functionality for Subsystem. Not tested.
- Made Examples/dynamic.py script to test some things out simply.

examples/dynamic.py

@@ -0,0 +1,129 @@
+# Hasty example with a Subsystem and the new dynamic features
+
+import numpy as np
+import matplotlib.pyplot as plt
+import moorpy as mp
+
+
+# ===== Create a MoorPy system to test with =====
+
+# ----- choose some system geometry parameters -----
+
+depth     = 200                             # water depth [m]
+angles    = np.radians([60, 300])      # line headings list [rad]
+rAnchor   = 600                            # anchor radius/spacing [m]
+zFair     = -21                             # fairlead z elevation [m]
+rFair     = 20                              # fairlead radius [m]
+lineLength= 650                            # line unstretched length [m]
+typeName  = "chain1"                        # identifier string for the line type
+
+
+# ----- Now a Subsystem in a larger System with a floating body -----
+
+# Create new MoorPy System and set its depth
+ms = mp.System(depth=depth)
+
+# add a line type
+ms.setLineType(dnommm=120, material='chain', name=typeName, Cd=1.4, Ca=1, CdAx=0.4, CaAx=0)  # this would be 120 mm chain
+
+# Add a free, body at [0,0,0] to the system (including some properties to make it hydrostatically stiff)
+ms.addBody(0, np.zeros(6), m=1e6, v=1e3, rM=100, AWP=1e3)
+
+# For each line heading, set the anchor point, the fairlead point, and the line itself
+for i, angle in enumerate(angles):
+
+    # create end Points for the line
+    ms.addPoint(1, [rAnchor*np.cos(angle), rAnchor*np.sin(angle), -depth])   # create anchor point (type 0, fixed)
+    ms.addPoint(1, [  rFair*np.cos(angle),   rFair*np.sin(angle),  zFair])   # create fairlead point (type 0, fixed)
+
+    # attach the fairlead Point to the Body (so it's fixed to the Body rather than the ground)
+    ms.bodyList[0].attachPoint(2*i+2, [rFair*np.cos(angle), rFair*np.sin(angle), zFair]) 
+
+    # add a Line going between the anchor and fairlead Points
+    ms.addLine(lineLength, typeName, pointA=2*i+1, pointB=2*i+2)
+
+# ----- Now add a SubSystem line! -----
+ss = mp.Subsystem(mooringSys=ms, depth=depth, spacing=rAnchor, rBfair=[10,0,-20])
+
+# set up the line types
+ms.setLineType(180, 'chain', name='one', Cd=1.4, Ca=1, CdAx=0.4, CaAx=0)
+ms.setLineType( 50, 'chain', name='two', Cd=1.4, Ca=1, CdAx=0.4, CaAx=0)
+
+
+# set up the lines and points and stuff
+ls = [350, 300]
+ts = ['one', 'two']
+ss.makeGeneric(lengths=ls, types=ts)
+ss.lineList[1].nNodes = 10  # reduce nodes on rope line for easier viewing
+ss.initialize()  # update things after changing node number
+
+# add points that the subSystem will attach to...
+ms.addPoint(1, [-rAnchor, 100, -depth])   # Point 5 - create anchor point (type 0, fixed)
+ms.addPoint(1, [ -rFair ,   0,  zFair])   # Point 6 - create fairlead point (type 0, fixed)
+ms.bodyList[0].attachPoint(6, [-rFair, 0, zFair])  # attach the fairlead Point to the Body 
+
+# string the Subsystem between the two points!
+ms.lineList.append(ss)  # add the SubSystem to the System's lineList
+ss.number = 3
+ms.pointList[4].attachLine(3, 0)  # attach it to the respective points
+ms.pointList[5].attachLine(3, 1)  # attach it to the respective points
+
+
+# ----- run the model to check that the Subsystem is working -----
+
+ms.initialize()                                             # make sure everything's connected
+
+ms.solveEquilibrium()                                       # equilibrate
+fig, ax = ms.plot()                                         # plot the system in original configuration
+
+ms.bodyList[0].f6Ext = np.array([3e6, 0, 0, 0, 0, 0])       # apply an external force on the body 
+ms.solveEquilibrium()                                      # equilibrate
+fig, ax = ms.plot(ax=ax, color='red')                       # plot the system in displaced configuration (on the same plot, in red)
+
+print(f"Body offset position is {ms.bodyList[0].r6}")
+
+
+
+# ===== Now try out Serag's dynamic tension functions with it =====
+
+# dynamic solve of some lines
+kbot = 3E+06
+cbot = 3E+05
+
+moorpy_freqs = []
+moorpy_fairten_psds = []
+moorpy_ten_stds = []
+
+omegas = np.array([ 0.02,  0.04,  0.06,  0.08 ])
+S_zeta = np.array([ 10.0 ,  10.0 ,  10.0 ,  10.0  ])
+RAO_fl = np.array([[ 2.0 ,  0.0 ,  0.0 ],
+                   [ 2.0 ,  0.0 ,  0.0 ],
+                   [ 2.0 ,  0.0 ,  0.0 ],
+                   [ 2.0 ,  0.0 ,  0.0 ]])
+
+T_nodes_psd_fd,T_nodes_std_fd,s,r_static,r_dynamic,r_total,X = ms.lineList[1].dynamicSolve(
+                                omegas,S_zeta,RAO_A=0,RAO_B=RAO_fl,depth=-ms.depth,
+                                kbot=kbot,cbot=cbot, seabed_tol=1e-4, tol=0.01, iters=100, w=0.8)
+
+fig2,ax2 = plt.subplots(1,1)
+
+plt.plot(T_nodes_std_fd,'-r',label = 'MoorPy (FD)')
+plt.xlabel('Node number (anchor to fairlead)')
+plt.ylabel('Fairlead tension std. dev [$N$]')
+
+# now try a Subsystem
+T_nodes_psd_fd,T_nodes_std_fd,s,r_static,r_dynamic,r_total,X = ms.lineList[2].dynamicSolve(
+                                omegas,S_zeta,RAO_A=0,RAO_B=RAO_fl,depth=-ms.depth,
+                                kbot=kbot,cbot=cbot, seabed_tol=1e-4, tol=0.01, iters=100, w=0.8)
+
+fig2,ax2 = plt.subplots(1,1)
+plt.plot(T_nodes_std_fd,'-r',label = 'MoorPy (FD)')
+plt.xlabel('Node number (anchor to fairlead)')
+plt.ylabel('Fairlead tension std. dev [$N$]')
+
+
+# Some mode shape plots
+ms.lineList[1].getModes(plot_modes=7, amp_factor=1000)  # with a Line
+ms.lineList[2].getModes(plot_modes=7, amp_factor=1000)  # with a Subsystem
+
+plt.show()

moorpy/dynamic_tension_functions.py

@@ -5,6 +5,7 @@
 import matplotlib.pyplot as plt
 from collections.abc import Iterable
 
+
 def get_horizontal_oop_vec(p1,p2):
     """
     Evaluates the horizontal out of plane vector given the coordinates of two points.
@@ -22,6 +23,7 @@ def get_horizontal_oop_vec(p1,p2):
         n_op = oop_vec/np.linalg.norm(oop_vec)
     return n_op
 
+
 def get_dynamic_matrices(Line, omegas, S_zeta,r_dynamic,depth,kbot,cbot,seabed_tol=1e-4):
     """
     Evaluates dynamic matrices for a Line object.
@@ -238,6 +240,7 @@ def get_dynamic_matrices(Line, omegas, S_zeta,r_dynamic,depth,kbot,cbot,seabed_t
 
     return M,A,B,K,r_mean,EA_segs
 
+
 def get_dynamic_tension(Line,omegas,S_zeta,RAO_A,RAO_B,depth,kbot,cbot,seabed_tol=1e-4,tol = 0.01,iters=100, w = 0.8, conv_time=True):
     """
     Evaluates dynamic tension at all the nodes for an instance of MoorPy's Line or CompositeLine classes.
@@ -365,6 +368,7 @@ def get_dynamic_tension(Line,omegas,S_zeta,RAO_A,RAO_B,depth,kbot,cbot,seabed_to
 
     return T_nodes_psd,T_nodes_std,s,r_static,r_dynamic,r_total,X
 
+
 def get_modes(line,fix_A=True,fix_B=True,plot_modes=False,amp_factor=1,adj_view = False,kbot=3E+06,cbot=3E+05,seabed_tol=1E-04):
     M,A,_,K,r_nodes,_ = line.getDynamicMatrices(np.ones(1), np.ones(1),0.,line.sys.depth,kbot,cbot,seabed_tol=seabed_tol)
 
@@ -377,7 +381,7 @@ def get_modes(line,fix_A=True,fix_B=True,plot_modes=False,amp_factor=1,adj_view
         n2 = -1
     else:
         n2 = r_nodes.shape[0]
-
+    
     eigvals,eigvecs = la.eig(K[3*n1:3*n2,3*n1:3*n2],M[3*n1:3*n2,3*n1:3*n2]+A[3*n1:3*n2,3*n1:3*n2])
     stable_eigvals = eigvals[np.real(eigvals)>0]
     stable_eigvecs = eigvecs[:,np.real(eigvals)>0]

moorpy/line.py

@@ -7,7 +7,7 @@
 from moorpy.helpers import (unitVector, LineError, CatenaryError, 
                      rotationMatrix, makeTower, read_mooring_file, 
                      quiver_data_to_segments, printVec, printMat)
-from moorpy.dynamic_tension_functions import get_dynamic_matrices,get_dynamic_tension,get_modes
+
 from os import path
 
 
@@ -1004,27 +1004,29 @@ def revertToStaticStiffness(self):
         # revert to original line length
         self.L = self.L0
 
-    def getDynamicMatrices(self, omegas, S_zeta,r_dynamic,depth,kbot,cbot,seabed_tol=1e-4):
-        M,A,B,K,r_mean,EA_segs = get_dynamic_matrices(self, omegas, S_zeta,r_dynamic,depth,kbot,cbot,seabed_tol=seabed_tol)
-        return M,A,B,K,r_mean,EA_segs
-
-    def dynamicSolve(self,omegas,S_zeta,RAO_A,RAO_B,depth,kbot,cbot,seabed_tol=1e-4,tol = 0.01,iters=100, w = 0.8):
-        T_nodes_psd,T_nodes_std,s,r_static,r_dynamic,r_total,X = get_dynamic_tension(self,omegas,S_zeta,RAO_A,RAO_B,depth,kbot,cbot,
-                                                                                     seabed_tol=seabed_tol,tol = tol,iters=iters, w = w)
-        return T_nodes_psd,T_nodes_std,s,r_static,r_dynamic,r_total,X
 
-    def getModes(self,fix_A=True,fix_B=True,plot_modes=False,amp_factor=1,adj_view = False,kbot=3E+06,cbot=3E+05,seabed_tol=1E-04):
+    def getDynamicMatrices(self, *args, **kwargs):
+        '''Compute M,A,B,K matrices for Line. See get_dynamic_matrices().'''
+        from moorpy.dynamic_tension_functions import get_dynamic_matrices
 
-        if plot_modes:
-            freqs,mode_shapes,r_nodes,M,A,K,fig,ax = get_modes(self,fix_A=fix_A,fix_B=fix_B,plot_modes=plot_modes,amp_factor=amp_factor,adj_view = adj_view,
-                                                               kbot=kbot,cbot=cbot,seabed_tol=seabed_tol)
-            return freqs,mode_shapes,r_nodes,M,A,K,fig,ax
-        else:
-            freqs,mode_shapes,r_nodes,M,A,K = get_modes(self,fix_A=fix_A,fix_B=fix_B,plot_modes=plot_modes,amp_factor=amp_factor,adj_view = adj_view,
-                                                        kbot=kbot,cbot=cbot,seabed_tol=seabed_tol)
-            return freqs,mode_shapes,r_nodes,M,A,K
+        return get_dynamic_matrices(self, *args, **kwargs)
 
 
+    def dynamicSolve(self, *args, **kwargs):
+        '''Compute complex amplitudes of line nodes. See get_dynamic_tension().'''
+        from moorpy.dynamic_tension_functions import get_dynamic_tension
+
+        return get_dynamic_tension(self, *args, **kwargs)
+
+
+    def getModes(self,*args, **kwargs):
+        '''Compute (and optionally plot) the line's mode shapes.
+        See get_modes().'''
+        from moorpy.dynamic_tension_functions import get_modes
+
+        return get_modes(self, *args, **kwargs)
+
+
 def from2Dto3Drotated(K2D, F, L, R): 
     '''Initialize a line end's analytic stiffness matrix in the 
     plane of the catenary then rotate the matrix to be about the 

moorpy/subsystem.py

@@ -160,7 +160,7 @@ def makeGeneric(self, lengths, types, points=[], shared=False):
         if shared:
             self.addPoint(-1, rA, DOFs=[2]) # add shared line point, free only to move in z
         else:
-            self.addPoint(-1, rA, DOFs=[0,2])                # add anchor point
+            self.addPoint(-1, rA, DOFs=[0,2])  # add anchor point
 
         # Go through each line segment and add its upper point, add the line, and connect the line to the points
         for i in range(self.nLines):
@@ -194,6 +194,88 @@ def makeGeneric(self, lengths, types, points=[], shared=False):
 
         # if a points list is provided, apply any mass or other properties it contains?
 
+        self.nLines = len(self.lineList)
+        self.nNodes = np.sum([line.nNodes for line in self.lineList]) - self.nLines + 1
+
+
+    def makeFromSystem(self, sys, point_id):
+        '''Populate a Subsystem based on a series of mooring lines in an
+        existing system (sys), starting at an end point (point_id).
+        Taken from Serag's CompositeLine class.
+        In this version, the subsystem still needs to be added to a system
+        afterward.
+        '''
+
+        if len(self.pointList)+len(self.lineList) > 0:
+            raise Exception('makeFromSystem can only be called for an empty Subsystem.')
+
+        point = sys.pointList[point_id-1] # Starting point id
+
+        # check starting point to make sure it is either coupled or fixed and that it is not connected to more than 1 line.
+        if point.type == 0:
+            raise Exception(f'Starting point ({point.number}) must not be free (change point type to -1 or 1).')
+        elif len(point.attached)>1:
+            raise Exception(f'Starting point ({point.number}) cannot be connected to more than 1 line')
+
+        # Save point A info
+        self.rA = np.array(point.r)
+        self.addPoint(1, self.rA)  # add anchor point
+
+        # Move through the points along the composite
+        while True:
+            # make sure that the point is free
+            if len(point.attached) > 2:
+                raise Exception(f'Point {point.number} is attached to more than two lines.')
+
+            # get the line id and line object
+            line_id = point.attached[-1] # get next line's id
+            line = sys.lineList[line_id - 1] # get line object
+            self.lineTypes[line.type['name']] = dict(line.type)  # copy the lineType dict
+            self.addLine(line.L0, self.lineTypes[line.type['name']]) # add the line
+
+
+            # get the next point
+            attached_points = line.attached.copy() # get the IDs of the points attached to the line
+            pointA_id = point.number # get first point ID
+            attached_points.remove(pointA_id) # remove first point from attached point list
+            pointB_id = attached_points[0] # get second point id
+            point = sys.pointList[pointB_id-1] # get second point object
+
+            if line(point.attached) == 1:  # must be the endpoint
+                self.addPoint(-1, point.r)
+            else:  # intermediate point along line
+                self.addPoint(0, point.r, DOFs=[0,2]) # may need to ensure there's no y component
+
+            # Attach the line ends to the points
+            self.pointList[-2].attachLine(len(self.lineList), 0)
+            self.pointList[-1].attachLine(len(self.lineList), 1)
+
+            # break from the loop when a point with a single attachment is reached
+            if len(point.attached) == 1:
+                break
+
+        # make sure that the last point is not a free point
+        if point.type == 0:
+            raise Exception(f'Last point ({point.number}) is a free point.')
+
+        self.rB = np.array(point.r)
+        self.nLines = len(self.lineList)
+        self.nNodes = np.sum([line.nNodes for line in self.lineList]) - self.nLines + 1
+
+
+    def initialize(self):
+        '''Initializes the subsystem objects to their initial positions, and
+        counts number of nodes.'''
+
+        self.nDOF, self.nCpldDOF, _ = self.getDOFs()
+
+        self.nNodes = np.sum([line.nNodes for line in self.lineList]) - self.nLines + 1
+
+        for point in self.pointList:
+            point.setPosition(point.r)
+
+        self.staticSolve()
+
 
     def setEndPosition(self, r, endB, sink=False):
         '''Sets either end position of the subsystem in the global/system
@@ -395,6 +477,59 @@ def revertToStaticStiffness(self):
         System.revertToStaticStiffness(self)
 
 
+    # ----- Function for dynamic frequency-domain tensions -----
+
+    def getDynamicMatrices(self, omegas, S_zeta, r_dynamic, depth, kbot, cbot,
+                           seabed_tol=1e-4):
+        '''Compute M,A,B,K matrices for the Subsystem. This calls 
+        get_dynamic_matrices() for each Line in the Subsystem then combines
+        the results. Note that this method overrides the Line method. Other
+        Line methods used for dynamics can be used directly in Subsystem.
+        '''
+        self.nNodes = np.sum([line.nNodes for line in self.lineList]) - self.nLines + 1
+
+        EA_segs = np.zeros(self.nNodes-1) # extensional stiffness of the segments
+        n_dofs = 3*self.nNodes # number of dofs
+        M = np.zeros([n_dofs,n_dofs], dtype='float')
+        A = np.zeros([n_dofs,n_dofs], dtype='float')
+        B = np.zeros([n_dofs,n_dofs], dtype='float')
+        K = np.zeros([n_dofs,n_dofs], dtype='float')
+        r_mean = np.zeros([self.nNodes,3], dtype='float')
+        r_dynamic = np.ones((len(omegas),self.nNodes,3),dtype='float')*r_dynamic
+        v_dynamic = 1j*omegas[:,None,None]*r_dynamic
+
+        n = 0  # starting index of the next line's entries in the matrices
+
+        for line in self.lineList:
+            n1 = int(n/3) 
+            n2 = n1 + line.nNodes
+
+            # Filling matrices for line (dof n to dof 3xline_nodes+n)
+            M_n,A_n,B_n,K_n,r_n,EA_segs_n = line.getDynamicMatrices(omegas, S_zeta,r_dynamic[:,n1:n2,:],depth,kbot,cbot,seabed_tol=seabed_tol)
+            M[n:3*line.nNodes+n,n:3*line.nNodes+n] += M_n
+            A[n:3*line.nNodes+n,n:3*line.nNodes+n] += A_n
+            B[n:3*line.nNodes+n,n:3*line.nNodes+n] += B_n
+            K[n:3*line.nNodes+n,n:3*line.nNodes+n] += K_n
+
+            # Attachment point properties
+            attachment = self.pointList[line.attached[-1]-1] # attachment point
+            attachment_idx = n2 - 1 # last node index
+            sigma_vp = np.sqrt(np.trapz(np.abs(v_dynamic[:,attachment_idx,:])**2*S_zeta[:,None],omegas,axis=0)) # standard deviations of the global components of the attachment point's velocity
+
+            M[3*line.nNodes - 3:3*line.nNodes, 3*line.nNodes - 3:3*line.nNodes] += attachment.m*np.eye(3) 
+            A[3*line.nNodes - 3:3*line.nNodes, 3*line.nNodes - 3:3*line.nNodes] += attachment.Ca* self.rho * attachment.v * np.eye(3)
+            B[3*line.nNodes - 3:3*line.nNodes, 3*line.nNodes - 3:3*line.nNodes] += 0.5* self.rho * attachment.CdA * np.pi*(3*attachment.v/4/np.pi)**(2/3) \
+                                                                                   * np.eye(3) * np.sqrt(8/np.pi) * np.diag(sigma_vp)
+
+            # Static line properties
+            r_mean[n1:n2,:] = r_n
+            EA_segs[n1:n2-1] = EA_segs_n 
+
+            n += 3*line.nNodes - 3 # next line starting node add the number of dofs of the current line minus 3 to get the last shared node
+
+        return M,A,B,K,r_mean,EA_segs
+
+
     # ---- Extra convenience functions (subsystem should be in equilibrium) -----