diff --git a/moorpy/compositeLine.py b/moorpy/compositeLine.py deleted file mode 100644 index bb45e31..0000000 --- a/moorpy/compositeLine.py +++ /dev/null @@ -1,111 +0,0 @@ -import numpy as np -from moorpy.dynamic_tension_functions import get_dynamic_tension, get_modes - -class CompositeLine(): - def __init__(self, sys, point_id, rho=1025.): - - self.rho = rho - 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 ValueError(f'Starting point ({point.number}) must not be free (change point type to -1 or 1).') - elif len(point.attached)>1: - raise ValueError(f'Starting point ({point.number}) cannot be connected to more than 1 line') - - self.pointList = [] - self.lineList = [] - - self.pointA = point - - # Move through the points along the composite - while True: - # make sure that the point is free - if len(point.attached)>2: - raise ValueError(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 - - # append the starting point and the line object - self.pointList.append(point) - self.lineList.append(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 - - # break from the loop when a point with a single attachment is reached - if len(point.attached) == 1: - self.pointList.append(point) - break - # make sure that the last point is not a free point - if point.type == 0: - raise ValueError(f'Last point ({point.number}) is a free point.') # make sure that the last point is not free - - self.pointB = point - self.nNodes = (np.sum([line.nNodes for line in self.lineList]) - (len(self.lineList)-1)) - self.sys = sys - - def getDynamicMatrices(self, omegas, S_zeta,r_dynamic,depth,kbot,cbot,seabed_tol=1e-4): - n_nodes = self.nNodes # number of nodes - EA_segs = np.zeros(n_nodes-1) # extensional stiffness of the segments - n_dofs = 3*n_nodes # 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') - n = 0 - r_mean = np.zeros([n_nodes,3], dtype='float') - r_dynamic = np.ones((len(omegas),n_nodes,3),dtype='float')*r_dynamic - v_dynamic = 1j*omegas[:,None,None]*r_dynamic - - for line in self.lineList: - line_type = line.type - 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 - - 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): - - 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 \ No newline at end of file diff --git a/moorpy/dynamic_tension_functions.py b/moorpy/dynamic_tension_functions.py deleted file mode 100644 index 720743c..0000000 --- a/moorpy/dynamic_tension_functions.py +++ /dev/null @@ -1,455 +0,0 @@ -import numpy as np -from numpy.linalg import solve,norm -import scipy.linalg as la -from datetime import datetime -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. - """ - hor_vec = p2 - np.array([p1[0],p1[1],p2[2]]) - ver_vec = p1 - np.array([p1[0],p1[1],p2[2]]) - - if np.isclose(np.linalg.norm(hor_vec),0): # vertical line - n_op = np.array([1,0,0]) - elif np.isclose(np.linalg.norm(ver_vec),0): # horizontal line - oop_vec = np.cross(hor_vec,np.array([0,0,1])) - n_op = oop_vec/np.linalg.norm(oop_vec) - else: - oop_vec = np.cross(hor_vec,ver_vec) - 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. - - Parameters - ---------- - Line : Line - An instance of MoorPy's Line class - omegas : ndarray - Array of frequencies in rad/s. - S_zeta : ndarray - Wave spectrum array in m^2/(rad/s), must be of the same length as omegas. - r_dynamic : ndarray - A 3d array of the frequency dependent complex amplitudes of line nodes. The array has a shape of (m,n,3) where m is the number of frequencies, - n is the number of nodes, and the last dimension of 3 correspond to the x,y,z components. - depth : float - Water depth. - kbot : float - Vertical stiffness for points lying on the seabed. - cbot : float - Vertical damping for points lying on the seabed. - seabed_tol : float, optional - Distance from seabed within which a node is considered to be lying on the seabed, by default 1e-4 m. - - Returns - ------- - M: ndarray - Mass matrix. - A: ndarray - Added mass matrix. - B: ndarray - Damping matrix. - K: ndarray - Stiffness matrix. - M: ndarray - Mass matrix. - r_mean: ndarray - Mean static positions of the nodes given as a (m,3) array where m is the number of nodes. - EA_segs: ndarray - Extensional stiffness of segments. - """ - # extract line properties - N = Line.nNodes - mden = Line.type['m'] # line mass density function - deq = Line.type['d_vol'] # line volume equivalent diameter - EA = Line.type['EA'] # extensional stiffness - Ca = Line.type['Ca'] # normal added mass coeff - CaAx = Line.type['CaAx'] # tangential added mass coeff - Cd = Line.type['Cd'] # normal drag coeff - CdAx = Line.type['CdAx'] # tangential drag coeff - - # extract mean node coordinates - X_mean,Y_mean,Z_mean,T_mean = Line.getLineCoords(0.0,n=N) # coordinates of line nodes and tension values - r_mean = np.vstack((X_mean,Y_mean,Z_mean)).T # coordinates of line nodes - - # evaluate node velocities - r_dynamic = np.ones((len(omegas),N,3),dtype='float')*r_dynamic - v_dynamic = 1j*omegas[:,None,None]*r_dynamic - - # define out of plane normal - h_op = get_horizontal_oop_vec(r_mean[0],r_mean[-1]) # horizontal out-of-plane vector - hh_op = np.outer(h_op,h_op) - - # intialize line matrices - M = np.zeros([3*N, 3*N], dtype='float') # mass matrix - A = np.zeros([3*N, 3*N], dtype='float') # added mass matrix - B = np.zeros([3*N, 3*N], dtype='float') # linearized viscous damping matrix - K = np.zeros([3*N, 3*N], dtype='float') # stiffness matrix - - # Node 0 - dr_e1 = r_mean[1] - r_mean[0] - L_e1 = la.norm(dr_e1) # element 1 length - t_e1 = (dr_e1)/L_e1 # tangential unit vector - p_e1 = np.cross(t_e1,h_op) # in plane normal unit vector - - - ut_e1 = np.einsum('ij,j->i',v_dynamic[:,0,:],t_e1) # tangential velocity - uh_e1 = np.einsum('ij,j->i',v_dynamic[:,0,:],h_op) # normal horizontal out of plane velocity - up_e1 = np.einsum('ij,j->i',v_dynamic[:,0,:],p_e1) # normal in plane velocity - - sigma_ut_e1 = np.sqrt(np.trapz(np.abs(ut_e1)**2*S_zeta,omegas)) - sigma_uh_e1 = np.sqrt(np.trapz(np.abs(uh_e1)**2*S_zeta,omegas)) - sigma_up_e1 = np.sqrt(np.trapz(np.abs(up_e1)**2*S_zeta,omegas)) - - tt_e1 = np.outer(t_e1,t_e1) # local tangential to global components transformation matrix - pp_e1 = np.outer(p_e1,p_e1) # local normal inplane to global components transformation matrix - - M[0:3,0:3] += mden*L_e1/2*np.eye(3) # element 1 mass contribution - - A_e1 = 1025*np.pi/4*deq**2*L_e1/2*(Ca*(hh_op+pp_e1) + CaAx*tt_e1) # element 1 added mass contribution - - B_e1 = 0.5*1025*deq*L_e1/2*np.sqrt(8/np.pi)*(Cd*(sigma_uh_e1*hh_op + sigma_up_e1*pp_e1) + - CdAx*sigma_ut_e1*tt_e1) # element 1 damping contribution - - K_e1 = EA/L_e1*tt_e1 + (T_mean[0]/L_e1)*(hh_op+pp_e1) # element 1 stiffness (axial + geometric) - - ## assembling element 1 contributions (rows corresponding to node 0) - A[0:3,0:3] += A_e1 - B[0:3,0:3] += B_e1 - K[0:3,0:3] += K_e1 - K[0:3,3:6] += -K_e1 - - ## add seabed contribution to node 0 - if np.isclose(r_mean[0,2],-depth,seabed_tol): - K[2,2] += kbot - B[2,2] += cbot - - # Internal nodes loop (each internal node has contributions from two elements n-1/2 and n+1/2) - for n in range(1, N-1): - - ## backward element (n-1/2) contributions - dr_bw = r_mean[n-1] - r_mean[n] - L_bw = la.norm(dr_bw) # element 1 length - t_bw = (dr_bw)/L_bw # tangential unit vector - p_bw = np.cross(t_bw,h_op) # in plane normal unit vector - - ut_bw = np.einsum('ij,j->i',v_dynamic[:,n,:],t_bw) # tangential velocity - uh_bw = np.einsum('ij,j->i',v_dynamic[:,n,:],h_op) # normal horizontal out of plane velocity - up_bw = np.einsum('ij,j->i',v_dynamic[:,n,:],p_bw) # normal in plane velocity - - sigma_ut_bw = np.sqrt(np.trapz(np.abs(ut_bw)**2*S_zeta,omegas)) - sigma_uh_bw = np.sqrt(np.trapz(np.abs(uh_bw)**2*S_zeta,omegas)) - sigma_up_bw = np.sqrt(np.trapz(np.abs(up_bw)**2*S_zeta,omegas)) - - tt_bw = np.outer(t_bw,t_bw) # local tangential to global components transformation matrix - pp_bw = np.outer(p_bw,p_bw) # local normal inplane to global components transformation matrix - - M[3*n:3*n+3,3*n:3*n+3] += mden*L_bw/2*np.eye(3) # mass contribution from adjacent elements - - A_bw = 1025*np.pi/4*deq**2*L_bw/2*(Ca*(hh_op+pp_bw) + CaAx*tt_bw) # backward element added mass contribution - - B_bw = 0.5*1025*deq*L_bw/2*np.sqrt(8/np.pi)*(Cd*(sigma_uh_bw*hh_op + sigma_up_bw*pp_bw) + - CdAx*sigma_ut_bw*tt_bw) # backward element damping contribution - - K_bw = EA/L_bw*tt_bw + (T_mean[n]/L_bw)*(hh_op+pp_bw) # backward element stiffness (axial + geometric) - - ## forward element (n+1/2) contributions - dr_fw = r_mean[n+1] - r_mean[n] - L_fw = la.norm(dr_fw) # element 1 length - t_fw = (dr_fw)/L_fw # tangential unit vector - p_fw = np.cross(t_fw,h_op) # in plane normal unit vector - - - ut_fw = np.einsum('ij,j->i',v_dynamic[:,n,:],t_fw) # tangential velocity - uh_fw = np.einsum('ij,j->i',v_dynamic[:,n,:],h_op) # normal horizontal out of plane velocity - up_fw = np.einsum('ij,j->i',v_dynamic[:,n,:],p_fw) # normal in plane velocity - - sigma_ut_fw = np.sqrt(np.trapz(np.abs(ut_fw)**2*S_zeta,omegas)) - sigma_uh_fw = np.sqrt(np.trapz(np.abs(uh_fw)**2*S_zeta,omegas)) - sigma_up_fw = np.sqrt(np.trapz(np.abs(up_fw)**2*S_zeta,omegas)) - - tt_fw = np.outer(t_fw,t_fw) # local tangential to global components transformation matrix - pp_fw = np.outer(p_fw,p_fw) # local normal inplane to global components transformation matrix - - M[3*n:3*n+3,3*n:3*n+3] += mden*L_fw/2*np.eye(3) # mass contribution from adjacent elements - - A_fw = 1025*np.pi/4*deq**2*L_fw/2*(Ca*(hh_op+pp_fw) + CaAx*tt_fw) # forward element added mass contribution - - B_fw = 0.5*1025*deq*L_fw/2*np.sqrt(8/np.pi)*(Cd*(sigma_uh_fw*hh_op + sigma_up_fw*pp_fw) + - CdAx*sigma_ut_fw*tt_fw) # forward element damping contribution - - K_fw = EA/L_fw*tt_fw + (T_mean[n]/L_fw)*(hh_op+pp_fw) # forward element stiffness (axial + geometric) - - ## assembling bwd and fwd elements contributions (rows corresponding to node n) - A[3*n:3*n+3,3*n:3*n+3] += A_bw + A_fw - B[3*n:3*n+3,3*n:3*n+3] += B_bw + B_fw - K[3*n:3*n+3,3*n:3*n+3] += K_bw + K_fw - K[3*n:3*n+3,3*n-3:3*n] += -K_bw - K[3*n:3*n+3,3*n+3:3*n+6] += -K_fw - - ## add seabed contribution to node n - if np.isclose(r_mean[n,2],-depth,seabed_tol): - K[3*n+2,3*n+2] += kbot - B[3*n+2,3*n+2] += cbot - - # Node N - dr_eN = r_mean[N-1] - r_mean[N-2] - L_eN = la.norm(dr_eN) # element N length - t_eN = (dr_eN)/L_eN # tangential unit vector - p_eN = np.cross(t_eN,h_op) # in plane normal unit vector - - ut_eN = np.einsum('ij,j->i',v_dynamic[:,N-1,:],t_eN) # tangential velocity - uh_eN = np.einsum('ij,j->i',v_dynamic[:,N-1,:],h_op) # normal horizontal out of plane velocity - up_eN = np.einsum('ij,j->i',v_dynamic[:,N-1,:],p_eN) # normal in plane velocity - - sigma_ut_eN = np.sqrt(np.trapz(np.abs(ut_eN)**2*S_zeta,omegas)) - sigma_uh_eN = np.sqrt(np.trapz(np.abs(uh_eN)**2*S_zeta,omegas)) - sigma_up_eN = np.sqrt(np.trapz(np.abs(up_eN)**2*S_zeta,omegas)) - - tt_eN = np.outer(t_eN,t_eN) # local tangential to global components transformation matrix - pp_eN = np.outer(p_eN,p_eN) # local normal inplane to global components transformation matrix - - M[3*(N-1):3*(N-1)+3,3*(N-1):3*(N-1)+3] += mden*L_eN/2*np.eye(3) # element N mass contribution - - A_eN = 1025*np.pi/4*deq**2*L_eN/2*(Ca*(hh_op+pp_eN) + CaAx*tt_eN) # element N added mass contribution - - B_eN = 0.5*1025*deq*L_eN/2*np.sqrt(8/np.pi)*(Cd*(sigma_uh_eN*hh_op + sigma_up_eN*pp_eN) + - CdAx*sigma_ut_eN*tt_eN) # element N damping contribution - - K_eN = EA/L_eN*tt_eN + (T_mean[N-1]/L_eN)*(hh_op+pp_eN) # element N stiffness (axial + geometric) - - ## assembling element N contributions (rows corresponding to node N) - A[3*(N-1):3*(N-1)+3,3*(N-1):3*(N-1)+3] += A_eN - B[3*(N-1):3*(N-1)+3,3*(N-1):3*(N-1)+3] += B_eN - K[3*(N-1):3*(N-1)+3,3*(N-1):3*(N-1)+3] += K_eN - K[3*(N-1):3*(N-1)+3,3*(N-1)-3:3*(N-1)] += -K_eN - - ## add seabed contribution to node N - if np.isclose(r_mean[N-1,2],-depth,seabed_tol): - K[3*(N-1)+2,3*(N-1)+2] += kbot - B[3*(N-1)+2,3*(N-1)+2] += cbot - - EA_segs = Line.type['EA']*np.ones(Line.nNodes - 1) - - 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. - - Parameters - ---------- - Line : Line/CompositeLine - An instance of MoorPy's Line or CompositeLine classes. - omegas : ndarray - Array of frequencies in rad/s. - S_zeta : ndarray - Wave spectrum array in m^2/(rad/s), must be of the same length as omegas. - RAO_A : ndarray - Translational RAOs for end node A given as a (m,3) array where m is the number of frequencies (must be equal to the length of omegas) . - RAO_B : ndarray - Translational RAOs for end node B given as a (m,3) array where m is the number of frequencies (must be equal to the length of omegas) . - depth : float - Water depth. - kbot : float - Vertical stiffness for points lying on the seabed. - cbot : float - Vertical damping for points lying on the seabed. - seabed_tol : float, optional - Distance from seabed within which a node is considered to be lying on the seabed, by default 1e-4 m. - tol : float, optional - Relative tolerance for iteration, by default 0.01 - iters : int, optional - Maximum number of iterations, by default 100 - w : float, optional - Succesive relaxation coefficient, by default 0.8 - - Returns - ------- - T_nodes_psd: ndarray - Tension spectra at nodes given as (m,n) array where m is the number of frequencies and n is the number of nodes. - T_nodes_std: ndarray - Tension standard deviation at nodes. - s: ndarray: - Node locations along the line. - r_static: ndarray - Nodes mean static position given as (n,3) array where n the number of nodes. - r_dynamic: ndarray - Nodes complex dynamic amplitudes given as (m,n,3) array where m the number of frequencies, n is the number of nodes - r_total: ndarray - Combined static and dynamic positions. - X: ndarray - Solution of the dynamic problem. - """ - N = Line.nNodes - n_dofs = 3*N - - if np.all(RAO_A == 0): - RAO_A = np.zeros_like(RAO_B) - if np.all(RAO_B == 0): - RAO_B = np.zeros_like(RAO_A) - - # intialize iteration matrices - r_dynamic_init = np.ones((len(omegas),N,3)) - M,A,B,K,r_static,EA_segs = Line.getDynamicMatrices(omegas,S_zeta,r_dynamic_init,depth,kbot,cbot,seabed_tol=seabed_tol) # TODO: return EA_segs - X = np.zeros((len(omegas),n_dofs),dtype = 'complex') - r_dynamic = np.zeros(((len(omegas),int(n_dofs/3),3)),dtype = 'complex') - S_Xd = np.zeros((len(omegas),n_dofs),dtype = 'float') - sigma_Xd = np.zeros(n_dofs,dtype = 'float') - sigma_Xd0 = np.zeros(n_dofs,dtype = 'float') - X[:, :3] = RAO_A - X[:,-3:] = RAO_B - - # solving dynamics - start = datetime.now() - for ni in range(iters): - H = - omegas[:,None,None]**2*(M+A)[None,:,:]\ - + 1j*omegas[:,None,None]*B[None,:,:]\ - + K[None,:,:]\ - - F_A = np.einsum('nij,njk->ni',-H[:,3:-3, :3],RAO_A[:,:,None]) - F_B = np.einsum('nij,njk->ni',-H[:,3:-3,-3:],RAO_B[:,:,None]) - F = F_A + F_B - - X[:,3:-3] = solve(H[:,3:-3,3:-3],F) - S_Xd[:] = np.abs(1j*omegas[:,None]*X)**2*S_zeta[:,None] - sigma_Xd[:] = np.sqrt(np.trapz(S_Xd,omegas,axis=0)) - r_dynamic[:] = X.reshape(X.shape[0],int(X.shape[1]/3),3) - if (np.abs(sigma_Xd-sigma_Xd0) <= tol*np.abs(sigma_Xd0)).all(): - break - else: - sigma_Xd0[:] = w * sigma_Xd + (1.-w) * sigma_Xd0 - _,_,B[:],_,_,_ = Line.getDynamicMatrices(omegas,S_zeta,r_dynamic,depth,kbot,cbot,seabed_tol=seabed_tol) - if conv_time: - print(f'Finished {ni} dynamic tension iterations in {datetime.now()-start} seconds (w = {w}).') - - # evaluate tension - dw = np.diff(omegas, - prepend= omegas[0] - (omegas[1]-omegas[0]), - append= omegas[-1] + (omegas[-1]-omegas[-2])) - dw = (dw[1:]+dw[:-1])/2 - wave_amps = np.sqrt(S_zeta*dw) #evaluate wave amplitudes of harmonic components from wave spectrum - - r_dynamic *= wave_amps[:,None,None] - r_total = r_static[None,:,:] + r_dynamic - dr_static = r_static[:-1] - r_static[1:] - dr_dynamic = r_dynamic[:,:-1,:] - r_dynamic[:,1:,:] - tangents = dr_static/la.norm(r_static[:-1] - r_static[1:], axis=-1)[:,None] - L_static = la.norm(dr_static, axis=-1) - dL_dynamic = np.einsum('mni,ni->mn', dr_dynamic, tangents) - eps_segs = np.abs(dL_dynamic)/L_static - - T_segs = EA_segs * eps_segs - T_nodes = np.zeros((len(omegas),N)) - T_nodes[:,0] = T_segs[:,0] - T_nodes[:,1:-1] = (T_segs[:,:-1] + T_segs[:,1:])/2 - T_nodes[:,-1] = T_segs[:,-1] - - # S_T = np.zeros((len(omegas),N)) - # S_T[:,1:] = T_e**2/dw[:,None] - # S_T[:,0] = S_T[:,1] - - T_nodes_psd = T_nodes**2/dw[:,None] - T_nodes_std = np.sqrt(np.trapz(T_nodes_psd,omegas,axis=0)) - - - dr = np.diff(r_static,axis=0) - ds = la.norm(dr,axis=1) - s = np.zeros_like(T_nodes_std) - s = np.cumsum(ds) - - 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) - - if fix_A: - n1 = 1 - else: - n1 = 0 - - if fix_B: - 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] - - idx = stable_eigvals.argsort()[::-1] - stable_eigvals = stable_eigvals[idx] - stable_eigvecs = stable_eigvecs[:,idx] - - freqs = np.sqrt(np.real(stable_eigvals))/2/np.pi - mode_shapes = np.zeros(stable_eigvecs.shape,dtype='float') - - for i in range(stable_eigvecs.shape[1]): - mode_shapes[:,i] = r_nodes[n1:n2].flatten('C') + stable_eigvecs[:,i] - - freqs = np.flip(freqs) - mode_shapes = np.flip(mode_shapes,axis=1) - - if plot_modes: - cols = 4 - rows = plot_modes//cols + bool(plot_modes%cols) - fig,ax = plt.subplots(rows,cols,subplot_kw={"projection": "3d"},figsize=(5*cols,5*rows)) - - i = 0 - for axes in ax: - if not isinstance(axes,Iterable): - axes = [axes] - - for axis in axes: - if i >= plot_modes: - break - - r = r_nodes.copy() - r[n1:n2] = mode_shapes[:,i].reshape([int(len(mode_shapes[:,i])/3),3]) - r = (r-r_nodes)*amp_factor - x = r[:,0] - y = r[:,1] - z = r[:,2] - - x_0 = r_nodes[:,0] - y_0 = r_nodes[:,1] - z_0 = r_nodes[:,2] - - axis.plot(x_0,y_0,z_0,'-ko',label='initial') - axis.plot(x+x_0,y+y_0,z+z_0,'--ro',label='mode shape') - - # R_0 = np.sqrt(x_0**2 + y_0**2) - if adj_view: - # h_min = np.min((x_0,y_0)) - # h_max = np.max((x_0,y_0)) - # axis.set_xlim(h_min,h_max) - # axis.set_ylim(h_min,h_max) - # axis.set_zlim(z_0.min(),z_0.max()) - sigma_x = x.std() - sigma_y = y.std() - sigma_z = z.std() - azim = np.arctan2(sigma_x,sigma_y)*180/np.pi - elev = np.arctan2(np.hypot(sigma_x,sigma_y),sigma_z)*180/np.pi - axis.view_init(elev=elev,azim=azim) - - # axis.set_box_aspect([np.ptp(coord) for coord in [x, y, z]]) - axis.set_xlabel('X (m)') - axis.set_ylabel('Y (m)') - axis.set_zlabel('Z (m)') - axis.set_title(f'f = {freqs[i]:.3f} Hz, T = {1/freqs[i]:.3f} s') - - i+=1 - - fig.tight_layout() - return freqs,mode_shapes,r_nodes,M,A,K,fig,ax - else: - return freqs,mode_shapes,r_nodes,M,A,K diff --git a/moorpy/helpers.py b/moorpy/helpers.py index b6d6ad9..03b789e 100644 --- a/moorpy/helpers.py +++ b/moorpy/helpers.py @@ -1011,7 +1011,8 @@ def lines2subsystem(lines,ms,span=None,case=0): ss.lineTypes[types[-1]] = ms.lineTypes[types[-1]] # use makeGeneric to build the subsystem line - ss.makeGeneric(lengths,types,suspended=case) + nSegs = [ms.lineList[i].nNodes-1 for i in lines] + ss.makeGeneric(lengths,types,suspended=case, nSegs=nSegs) ss.setEndPosition(ms.lineList[lines[0]].rA,endB=0) ss.setEndPosition(ms.lineList[lines[-1]].rB,endB=1) @@ -1467,4 +1468,456 @@ def read_output_file(dirName,fileName, skiplines=-1, hasunits=1, chanlim=999, di unitDict[channels[i]] = units[i] return dataDict, unitDict else: - return data3, ch, channels, units \ No newline at end of file + return data3, ch, channels, units + + +def get_horizontal_oop_vec(p1,p2): + """ + Evaluates the horizontal out of plane vector given the coordinates of two points. + """ + hor_vec = p2 - np.array([p1[0],p1[1],p2[2]]) + ver_vec = p1 - np.array([p1[0],p1[1],p2[2]]) + + if np.isclose(np.linalg.norm(hor_vec),0): # vertical line + n_op = np.array([1,0,0]) + elif np.isclose(np.linalg.norm(ver_vec),0): # horizontal line + oop_vec = np.cross(hor_vec,np.array([0,0,1])) + n_op = oop_vec/np.linalg.norm(oop_vec) + else: + oop_vec = np.cross(hor_vec,ver_vec) + 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. + + Parameters + ---------- + Line : Line + An instance of MoorPy's Line class + omegas : ndarray + Array of frequencies in rad/s. + S_zeta : ndarray + Wave spectrum array in m^2/(rad/s), must be of the same length as omegas. + r_dynamic : ndarray + A 3d array of the frequency dependent complex amplitudes of line nodes. The array has a shape of (m,n,3) where m is the number of frequencies, + n is the number of nodes, and the last dimension of 3 correspond to the x,y,z components. + depth : float + Water depth. + kbot : float + Vertical stiffness for points lying on the seabed. + cbot : float + Vertical damping for points lying on the seabed. + seabed_tol : float, optional + Distance from seabed within which a node is considered to be lying on the seabed, by default 1e-4 m. + + Returns + ------- + M: ndarray + Mass matrix. + A: ndarray + Added mass matrix. + B: ndarray + Damping matrix. + K: ndarray + Stiffness matrix. + M: ndarray + Mass matrix. + r_mean: ndarray + Mean static positions of the nodes given as a (m,3) array where m is the number of nodes. + EA_segs: ndarray + Extensional stiffness of segments. + """ + # extract line properties + N = Line.nNodes + mden = Line.type['m'] # line mass density function + deq = Line.type['d_vol'] # line volume equivalent diameter + EA = Line.type['EA'] # extensional stiffness + Ca = Line.type['Ca'] # normal added mass coeff + CaAx = Line.type['CaAx'] # tangential added mass coeff + Cd = Line.type['Cd'] # normal drag coeff + CdAx = Line.type['CdAx'] # tangential drag coeff + + # extract mean node coordinates + X_mean,Y_mean,Z_mean,T_mean = Line.getLineCoords(0.0,n=N) # coordinates of line nodes and tension values + r_mean = np.vstack((X_mean,Y_mean,Z_mean)).T # coordinates of line nodes + + # evaluate node velocities + r_dynamic = np.ones((len(omegas),N,3),dtype='float')*r_dynamic + v_dynamic = 1j*omegas[:,None,None]*r_dynamic + + # define out of plane normal + h_op = get_horizontal_oop_vec(r_mean[0],r_mean[-1]) # horizontal out-of-plane vector + hh_op = np.outer(h_op,h_op) + + # intialize line matrices + M = np.zeros([3*N, 3*N], dtype='float') # mass matrix + A = np.zeros([3*N, 3*N], dtype='float') # added mass matrix + B = np.zeros([3*N, 3*N], dtype='float') # linearized viscous damping matrix + K = np.zeros([3*N, 3*N], dtype='float') # stiffness matrix + + # Node 0 + dr_e1 = r_mean[1] - r_mean[0] + L_e1 = np.linalg.norm(dr_e1) # element 1 length + t_e1 = (dr_e1)/L_e1 # tangential unit vector + p_e1 = np.cross(t_e1,h_op) # in plane normal unit vector + + + ut_e1 = np.einsum('ij,j->i',v_dynamic[:,0,:],t_e1) # tangential velocity + uh_e1 = np.einsum('ij,j->i',v_dynamic[:,0,:],h_op) # normal horizontal out of plane velocity + up_e1 = np.einsum('ij,j->i',v_dynamic[:,0,:],p_e1) # normal in plane velocity + + sigma_ut_e1 = np.sqrt(np.trapz(np.abs(ut_e1)**2*S_zeta,omegas)) + sigma_uh_e1 = np.sqrt(np.trapz(np.abs(uh_e1)**2*S_zeta,omegas)) + sigma_up_e1 = np.sqrt(np.trapz(np.abs(up_e1)**2*S_zeta,omegas)) + + tt_e1 = np.outer(t_e1,t_e1) # local tangential to global components transformation matrix + pp_e1 = np.outer(p_e1,p_e1) # local normal inplane to global components transformation matrix + + M[0:3,0:3] += mden*L_e1/2*np.eye(3) # element 1 mass contribution + + A_e1 = 1025*np.pi/4*deq**2*L_e1/2*(Ca*(hh_op+pp_e1) + CaAx*tt_e1) # element 1 added mass contribution + + B_e1 = 0.5*1025*deq*L_e1/2*np.sqrt(8/np.pi)*(Cd*(sigma_uh_e1*hh_op + sigma_up_e1*pp_e1) + + CdAx*sigma_ut_e1*tt_e1) # element 1 damping contribution + + K_e1 = EA/L_e1*tt_e1 + (T_mean[0]/L_e1)*(hh_op+pp_e1) # element 1 stiffness (axial + geometric) + + ## assembling element 1 contributions (rows corresponding to node 0) + A[0:3,0:3] += A_e1 + B[0:3,0:3] += B_e1 + K[0:3,0:3] += K_e1 + K[0:3,3:6] += -K_e1 + + ## add seabed contribution to node 0 + if np.isclose(r_mean[0,2],-depth,seabed_tol): + K[2,2] += kbot + B[2,2] += cbot + + # Internal nodes loop (each internal node has contributions from two elements n-1/2 and n+1/2) + for n in range(1, N-1): + + ## backward element (n-1/2) contributions + dr_bw = r_mean[n-1] - r_mean[n] + L_bw = np.linalg.norm(dr_bw) # element 1 length + t_bw = (dr_bw)/L_bw # tangential unit vector + p_bw = np.cross(t_bw,h_op) # in plane normal unit vector + + ut_bw = np.einsum('ij,j->i',v_dynamic[:,n,:],t_bw) # tangential velocity + uh_bw = np.einsum('ij,j->i',v_dynamic[:,n,:],h_op) # normal horizontal out of plane velocity + up_bw = np.einsum('ij,j->i',v_dynamic[:,n,:],p_bw) # normal in plane velocity + + sigma_ut_bw = np.sqrt(np.trapz(np.abs(ut_bw)**2*S_zeta,omegas)) + sigma_uh_bw = np.sqrt(np.trapz(np.abs(uh_bw)**2*S_zeta,omegas)) + sigma_up_bw = np.sqrt(np.trapz(np.abs(up_bw)**2*S_zeta,omegas)) + + tt_bw = np.outer(t_bw,t_bw) # local tangential to global components transformation matrix + pp_bw = np.outer(p_bw,p_bw) # local normal inplane to global components transformation matrix + + M[3*n:3*n+3,3*n:3*n+3] += mden*L_bw/2*np.eye(3) # mass contribution from adjacent elements + + A_bw = 1025*np.pi/4*deq**2*L_bw/2*(Ca*(hh_op+pp_bw) + CaAx*tt_bw) # backward element added mass contribution + + B_bw = 0.5*1025*deq*L_bw/2*np.sqrt(8/np.pi)*(Cd*(sigma_uh_bw*hh_op + sigma_up_bw*pp_bw) + + CdAx*sigma_ut_bw*tt_bw) # backward element damping contribution + + K_bw = EA/L_bw*tt_bw + (T_mean[n]/L_bw)*(hh_op+pp_bw) # backward element stiffness (axial + geometric) + + ## forward element (n+1/2) contributions + dr_fw = r_mean[n+1] - r_mean[n] + L_fw = np.linalg.norm(dr_fw) # element 1 length + t_fw = (dr_fw)/L_fw # tangential unit vector + p_fw = np.cross(t_fw,h_op) # in plane normal unit vector + + + ut_fw = np.einsum('ij,j->i',v_dynamic[:,n,:],t_fw) # tangential velocity + uh_fw = np.einsum('ij,j->i',v_dynamic[:,n,:],h_op) # normal horizontal out of plane velocity + up_fw = np.einsum('ij,j->i',v_dynamic[:,n,:],p_fw) # normal in plane velocity + + sigma_ut_fw = np.sqrt(np.trapz(np.abs(ut_fw)**2*S_zeta,omegas)) + sigma_uh_fw = np.sqrt(np.trapz(np.abs(uh_fw)**2*S_zeta,omegas)) + sigma_up_fw = np.sqrt(np.trapz(np.abs(up_fw)**2*S_zeta,omegas)) + + tt_fw = np.outer(t_fw,t_fw) # local tangential to global components transformation matrix + pp_fw = np.outer(p_fw,p_fw) # local normal inplane to global components transformation matrix + + M[3*n:3*n+3,3*n:3*n+3] += mden*L_fw/2*np.eye(3) # mass contribution from adjacent elements + + A_fw = 1025*np.pi/4*deq**2*L_fw/2*(Ca*(hh_op+pp_fw) + CaAx*tt_fw) # forward element added mass contribution + + B_fw = 0.5*1025*deq*L_fw/2*np.sqrt(8/np.pi)*(Cd*(sigma_uh_fw*hh_op + sigma_up_fw*pp_fw) + + CdAx*sigma_ut_fw*tt_fw) # forward element damping contribution + + K_fw = EA/L_fw*tt_fw + (T_mean[n]/L_fw)*(hh_op+pp_fw) # forward element stiffness (axial + geometric) + + ## assembling bwd and fwd elements contributions (rows corresponding to node n) + A[3*n:3*n+3,3*n:3*n+3] += A_bw + A_fw + B[3*n:3*n+3,3*n:3*n+3] += B_bw + B_fw + K[3*n:3*n+3,3*n:3*n+3] += K_bw + K_fw + K[3*n:3*n+3,3*n-3:3*n] += -K_bw + K[3*n:3*n+3,3*n+3:3*n+6] += -K_fw + + ## add seabed contribution to node n + if np.isclose(r_mean[n,2],-depth,seabed_tol): + K[3*n+2,3*n+2] += kbot + B[3*n+2,3*n+2] += cbot + + # Node N + dr_eN = r_mean[N-1] - r_mean[N-2] + L_eN = np.linalg.norm(dr_eN) # element N length + t_eN = (dr_eN)/L_eN # tangential unit vector + p_eN = np.cross(t_eN,h_op) # in plane normal unit vector + + ut_eN = np.einsum('ij,j->i',v_dynamic[:,N-1,:],t_eN) # tangential velocity + uh_eN = np.einsum('ij,j->i',v_dynamic[:,N-1,:],h_op) # normal horizontal out of plane velocity + up_eN = np.einsum('ij,j->i',v_dynamic[:,N-1,:],p_eN) # normal in plane velocity + + sigma_ut_eN = np.sqrt(np.trapz(np.abs(ut_eN)**2*S_zeta,omegas)) + sigma_uh_eN = np.sqrt(np.trapz(np.abs(uh_eN)**2*S_zeta,omegas)) + sigma_up_eN = np.sqrt(np.trapz(np.abs(up_eN)**2*S_zeta,omegas)) + + tt_eN = np.outer(t_eN,t_eN) # local tangential to global components transformation matrix + pp_eN = np.outer(p_eN,p_eN) # local normal inplane to global components transformation matrix + + M[3*(N-1):3*(N-1)+3,3*(N-1):3*(N-1)+3] += mden*L_eN/2*np.eye(3) # element N mass contribution + + A_eN = 1025*np.pi/4*deq**2*L_eN/2*(Ca*(hh_op+pp_eN) + CaAx*tt_eN) # element N added mass contribution + + B_eN = 0.5*1025*deq*L_eN/2*np.sqrt(8/np.pi)*(Cd*(sigma_uh_eN*hh_op + sigma_up_eN*pp_eN) + + CdAx*sigma_ut_eN*tt_eN) # element N damping contribution + + K_eN = EA/L_eN*tt_eN + (T_mean[N-1]/L_eN)*(hh_op+pp_eN) # element N stiffness (axial + geometric) + + ## assembling element N contributions (rows corresponding to node N) + A[3*(N-1):3*(N-1)+3,3*(N-1):3*(N-1)+3] += A_eN + B[3*(N-1):3*(N-1)+3,3*(N-1):3*(N-1)+3] += B_eN + K[3*(N-1):3*(N-1)+3,3*(N-1):3*(N-1)+3] += K_eN + K[3*(N-1):3*(N-1)+3,3*(N-1)-3:3*(N-1)] += -K_eN + + ## add seabed contribution to node N + if np.isclose(r_mean[N-1,2],-depth,seabed_tol): + K[3*(N-1)+2,3*(N-1)+2] += kbot + B[3*(N-1)+2,3*(N-1)+2] += cbot + + EA_segs = Line.type['EA']*np.ones(Line.nNodes - 1) + + 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. + + Parameters + ---------- + Line : Line/CompositeLine + An instance of MoorPy's Line or CompositeLine classes. + omegas : ndarray + Array of frequencies in rad/s. + S_zeta : ndarray + Wave spectrum array in m^2/(rad/s), must be of the same length as omegas. + RAO_A : ndarray + Translational RAOs for end node A given as a (m,3) array where m is the number of frequencies (must be equal to the length of omegas) . + RAO_B : ndarray + Translational RAOs for end node B given as a (m,3) array where m is the number of frequencies (must be equal to the length of omegas) . + depth : float + Water depth. + kbot : float + Vertical stiffness for points lying on the seabed. + cbot : float + Vertical damping for points lying on the seabed. + seabed_tol : float, optional + Distance from seabed within which a node is considered to be lying on the seabed, by default 1e-4 m. + tol : float, optional + Relative tolerance for iteration, by default 0.01 + iters : int, optional + Maximum number of iterations, by default 100 + w : float, optional + Succesive relaxation coefficient, by default 0.8 + + Returns + ------- + T_nodes_psd: ndarray + Tension spectra at nodes given as (m,n) array where m is the number of frequencies and n is the number of nodes. + T_nodes_std: ndarray + Tension standard deviation at nodes. + s: ndarray: + Node locations along the line. + r_static: ndarray + Nodes mean static position given as (n,3) array where n the number of nodes. + r_dynamic: ndarray + Nodes complex dynamic amplitudes given as (m,n,3) array where m the number of frequencies, n is the number of nodes + r_total: ndarray + Combined static and dynamic positions. + X: ndarray + Solution of the dynamic problem. + """ + N = Line.nNodes + n_dofs = 3*N + + if np.all(RAO_A == 0): + RAO_A = np.zeros_like(RAO_B) + if np.all(RAO_B == 0): + RAO_B = np.zeros_like(RAO_A) + + # intialize iteration matrices + r_dynamic_init = np.ones((len(omegas),N,3)) + M,A,B,K,r_static,EA_segs = Line.getDynamicMatrices(omegas,S_zeta,r_dynamic_init,depth,kbot,cbot,seabed_tol=seabed_tol) # TODO: return EA_segs + X = np.zeros((len(omegas),n_dofs),dtype = 'complex') + r_dynamic = np.zeros(((len(omegas),int(n_dofs/3),3)),dtype = 'complex') + S_Xd = np.zeros((len(omegas),n_dofs),dtype = 'float') + sigma_Xd = np.zeros(n_dofs,dtype = 'float') + sigma_Xd0 = np.zeros(n_dofs,dtype = 'float') + X[:, :3] = RAO_A + X[:,-3:] = RAO_B + + # solving dynamics + start = time.time() + for ni in range(iters): + H = - omegas[:,None,None]**2*(M+A)[None,:,:]\ + + 1j*omegas[:,None,None]*B[None,:,:]\ + + K[None,:,:]\ + + F_A = np.einsum('nij,njk->ni',-H[:,3:-3, :3],RAO_A[:,:,None]) + F_B = np.einsum('nij,njk->ni',-H[:,3:-3,-3:],RAO_B[:,:,None]) + F = F_A + F_B + + X[:,3:-3] = np.linalg.solve(H[:,3:-3,3:-3],F) + S_Xd[:] = np.abs(1j*omegas[:,None]*X)**2*S_zeta[:,None] + sigma_Xd[:] = np.sqrt(np.trapz(S_Xd,omegas,axis=0)) + r_dynamic[:] = X.reshape(X.shape[0],int(X.shape[1]/3),3) + if (np.abs(sigma_Xd-sigma_Xd0) <= tol*np.abs(sigma_Xd0)).all(): + break + else: + sigma_Xd0[:] = w * sigma_Xd + (1.-w) * sigma_Xd0 + _,_,B[:],_,_,_ = Line.getDynamicMatrices(omegas,S_zeta,r_dynamic,depth,kbot,cbot,seabed_tol=seabed_tol) + if conv_time: + print(f'Finished {ni} dynamic tension iterations in {time.time()-start} seconds (w = {w}).') + + # evaluate tension + dw = np.diff(omegas, + prepend= omegas[0] - (omegas[1]-omegas[0]), + append= omegas[-1] + (omegas[-1]-omegas[-2])) + dw = (dw[1:]+dw[:-1])/2 + wave_amps = np.sqrt(S_zeta*dw) #evaluate wave amplitudes of harmonic components from wave spectrum + + r_dynamic *= wave_amps[:,None,None] + r_total = r_static[None,:,:] + r_dynamic + dr_static = r_static[:-1] - r_static[1:] + dr_dynamic = r_dynamic[:,:-1,:] - r_dynamic[:,1:,:] + tangents = dr_static/np.linalg.norm(r_static[:-1] - r_static[1:], axis=-1)[:,None] + L_static = np.linalg.norm(dr_static, axis=-1) + dL_dynamic = np.einsum('mni,ni->mn', dr_dynamic, tangents) + eps_segs = np.abs(dL_dynamic)/L_static + + T_segs = EA_segs * eps_segs + T_nodes = np.zeros((len(omegas),N)) + T_nodes[:,0] = T_segs[:,0] + T_nodes[:,1:-1] = (T_segs[:,:-1] + T_segs[:,1:])/2 + T_nodes[:,-1] = T_segs[:,-1] + + # S_T = np.zeros((len(omegas),N)) + # S_T[:,1:] = T_e**2/dw[:,None] + # S_T[:,0] = S_T[:,1] + + T_nodes_psd = T_nodes**2/dw[:,None] + T_nodes_std = np.sqrt(np.trapz(T_nodes_psd,omegas,axis=0)) + + + dr = np.diff(r_static,axis=0) + ds = np.linalg.norm(dr,axis=1) + s = np.zeros_like(T_nodes_std) + s = np.cumsum(ds) + + 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): + import matplotlib.pyplot as plt + from collections.abc import Iterable + + M,A,_,K,r_nodes,_ = line.getDynamicMatrices(np.ones(1), np.ones(1),0.,line.sys.depth,kbot,cbot,seabed_tol=seabed_tol) + + if fix_A: + n1 = 1 + else: + n1 = 0 + + if fix_B: + 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] + + idx = stable_eigvals.argsort()[::-1] + stable_eigvals = stable_eigvals[idx] + stable_eigvecs = stable_eigvecs[:,idx] + + freqs = np.sqrt(np.real(stable_eigvals))/2/np.pi + mode_shapes = np.zeros(stable_eigvecs.shape,dtype='float') + + for i in range(stable_eigvecs.shape[1]): + mode_shapes[:,i] = r_nodes[n1:n2].flatten('C') + stable_eigvecs[:,i] + + freqs = np.flip(freqs) + mode_shapes = np.flip(mode_shapes,axis=1) + + if plot_modes: + cols = 4 + rows = plot_modes//cols + bool(plot_modes%cols) + fig,ax = plt.subplots(rows,cols,subplot_kw={"projection": "3d"},figsize=(5*cols,5*rows)) + + i = 0 + for axes in ax: + if not isinstance(axes,Iterable): + axes = [axes] + + for axis in axes: + if i >= plot_modes: + break + + r = r_nodes.copy() + r[n1:n2] = mode_shapes[:,i].reshape([int(len(mode_shapes[:,i])/3),3]) + r = (r-r_nodes)*amp_factor + x = r[:,0] + y = r[:,1] + z = r[:,2] + + x_0 = r_nodes[:,0] + y_0 = r_nodes[:,1] + z_0 = r_nodes[:,2] + + axis.plot(x_0,y_0,z_0,'-ko',label='initial') + axis.plot(x+x_0,y+y_0,z+z_0,'--ro',label='mode shape') + + # R_0 = np.sqrt(x_0**2 + y_0**2) + if adj_view: + # h_min = np.min((x_0,y_0)) + # h_max = np.max((x_0,y_0)) + # axis.set_xlim(h_min,h_max) + # axis.set_ylim(h_min,h_max) + # axis.set_zlim(z_0.min(),z_0.max()) + sigma_x = x.std() + sigma_y = y.std() + sigma_z = z.std() + azim = np.arctan2(sigma_x,sigma_y)*180/np.pi + elev = np.arctan2(np.hypot(sigma_x,sigma_y),sigma_z)*180/np.pi + axis.view_init(elev=elev,azim=azim) + + # axis.set_box_aspect([np.ptp(coord) for coord in [x, y, z]]) + axis.set_xlabel('X (m)') + axis.set_ylabel('Y (m)') + axis.set_zlabel('Z (m)') + axis.set_title(f'f = {freqs[i]:.3f} Hz, T = {1/freqs[i]:.3f} s') + + i+=1 + + fig.tight_layout() + return freqs,mode_shapes,r_nodes,M,A,K,fig,ax + else: + return freqs,mode_shapes,r_nodes,M,A,K diff --git a/moorpy/line.py b/moorpy/line.py index 721b5bf..ffd917c 100644 --- a/moorpy/line.py +++ b/moorpy/line.py @@ -6,7 +6,8 @@ from moorpy.nonlinear import nonlinear from moorpy.helpers import (unitVector, LineError, CatenaryError, rotationMatrix, makeTower, read_mooring_file, - quiver_data_to_segments, printVec, printMat) + quiver_data_to_segments, printVec, printMat, + get_dynamic_matrices, get_dynamic_tension, get_modes) from os import path @@ -1013,23 +1014,17 @@ def revertToStaticStiffness(self): 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 - 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) diff --git a/moorpy/subsystem.py b/moorpy/subsystem.py index f5874b1..0ffe509 100644 --- a/moorpy/subsystem.py +++ b/moorpy/subsystem.py @@ -119,7 +119,7 @@ def __init__(self, mooringSys=None, num=0, depth=0, rho=1025, g=9.81, self.qs = 1 # flag that it's a MoorPy analysis, so System methods don't complain. One day should replace this <<< - def makeGeneric(self, lengths, types, suspended=0): + def makeGeneric(self, lengths, types, suspended=0, nSegs=40): '''Creates a cable of n components going between an anchor point and a floating body (or a bridle point). If shared, it goes between two floating bodies. @@ -184,8 +184,16 @@ def makeGeneric(self, lengths, types, suspended=0): else: raise Exception(f"Can't find lineType '{types[i]}' in the SubSystem.") + # add the line segment using the reference to its lineType dict # add the line segment using the reference to its lineType dict - self.addLine(lengths[i], self.lineTypes[i]) + if nSegs is None: + self.addLine(lengths[i], self.lineTypes[i]) + elif isinstance(nSegs, (int, float)): + self.addLine(lengths[i], self.lineTypes[i], nSegs=nSegs) + elif isinstance(nSegs, list): + self.addLine(lengths[i], self.lineTypes[i], nSegs=nSegs[i]) + else: + raise ValueError("Invalid type for nSegs. Expected None, a number, or a list.") # add the upper end point of the segment if i==self.nLines-1: # if this is the upper-most line @@ -239,7 +247,7 @@ def makeFromSystem(self, sys, point_id): 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 + self.addLine(line.L0, self.lineTypes[line.type['name']], nSegs=line.nNodes-1) # add the line # get the next point @@ -249,7 +257,7 @@ def makeFromSystem(self, sys, point_id): 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 + if len(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 diff --git a/moorpy/system.py b/moorpy/system.py index 31763eb..7563101 100644 --- a/moorpy/system.py +++ b/moorpy/system.py @@ -18,7 +18,6 @@ from moorpy.point import Point from moorpy.line import Line from moorpy.lineType import LineType -from moorpy.compositeLine import CompositeLine import matplotlib as mpl #import moorpy.MoorSolve as msolve from moorpy.helpers import (rotationMatrix, rotatePosition, getH, printVec,