From 52412c990ceaf940d6be40f66c9452da6a713483 Mon Sep 17 00:00:00 2001
From: Christian Pavese <cpav@dtu.dk>
Date: Tue, 3 May 2016 09:48:15 +0200
Subject: [PATCH] Adding function to compute envelope of envelopes
---
wetb/prepost/Simulations.py | 253 ++++++++++++++++++++++++++++++------
1 file changed, 216 insertions(+), 37 deletions(-)
diff --git a/wetb/prepost/Simulations.py b/wetb/prepost/Simulations.py
index 75610151..894fb7d1 100755
--- a/wetb/prepost/Simulations.py
+++ b/wetb/prepost/Simulations.py
@@ -3294,6 +3294,181 @@ class WeibullParameters(object):
self.Vstep = 2.
self.shape_k = 2.
+def compute_env_of_env(envelope, dlc_list, Nx=300, Nsectors=12,Ntheta=181):
+ """
+ The function computes load envelopes for given channels and a groups of
+ load cases starting from the envelopes computed for single simulations.
+ The output is the envelope of the envelopes of the single simulations.
+ This total envelope is projected on defined polar directions.
+
+ Parameters
+ ----------
+
+ envelope : dict, dictionaries of interpolated envelopes of a given
+ channel (it's important that each entry of the dictonary
+ contains a matrix of the same dimensions). The dictonary
+ is organized by load case
+
+ dlc_list : list, list of load cases
+
+ Nx : int, default=300
+ Number of points for the envelope interpolation
+
+ Nsectors: int, default=12
+ Number of sectors in which the total envelope will be divided. The
+ default is every 30deg
+
+ Ntheta; int, default=181
+ Number of angles in which the envelope is interpolated in polar
+ coordinates.
+
+ Returns
+ -------
+
+ envelope : array (Nsectors x 6),
+ Total envelope projected on the number of angles defined in Nsectors.
+ The envelope is projected in Mx and My and the other cross-sectional
+ moments are fetched accordingly
+
+ """
+
+ # Group all the single DLCs
+ cloud = np.zeros(((Nx+1)*len(envelope),6))
+ for i in range(len(envelope)):
+ cloud[(Nx+1)*i:(Nx+1)*(i+1),:] = envelope[dlc_list[i]]
+ # Compute tital Hull of all the envelopes
+ hull = scipy.spatial.ConvexHull(cloud[:,:2])
+ cc = np.append(cloud[hull.vertices,:2],\
+ cloud[hull.vertices[0],:2].reshape(1,2),axis=0)
+ # Interpolate full envelope
+ cc_x,cc_up,cc_low,cc_int= int_envelope(cc[:,0],cc[:,1],Nx=Nx)
+ # Porject full envelope on given direction
+ cc_proj = proj_envelope(cc_x,cc_up,cc_low,cc_int,Nx,Nsectors,Ntheta)
+
+ env_proj = np.zeros([len(cc_proj),6])
+ env_proj[:,:2] = cc_proj
+
+ # Based on Mx and My, gather the remaining cross-sectional forces and
+ # moments
+ for ich in range(2, 6):
+ s0 = np.array(cloud[hull.vertices, ich]).reshape(-1, 1)
+ s1 = np.array(cloud[hull.vertices[0], ich]).reshape(-1, 1)
+ s0 = np.append(s0, s1, axis=0)
+ cc = np.append(cc, s0, axis=1)
+
+ _,_,_,extra_sensor = int_envelope(cc[:,0],cc[:,ich],Nx)
+ es = np.atleast_2d(np.array(extra_sensor[:,1])).T
+ cc_int = np.append(cc_int,es,axis=1)
+
+ for isec in range(Nsectors):
+ ids = (np.abs(cc_int[:,0]-cc_proj[isec,0])).argmin()
+ env_proj[isec,ich] = (cc_int[ids-1,ich]+cc_int[ids,ich]+\
+ cc_int[ids+1,ich])/3
+
+ return env_proj
+
+def int_envelope(ch1,ch2,Nx):
+ # Function to interpolate envelopes and output arrays of same length
+
+ # Number of points is defined by Nx + 1, where the + 1 is needed to
+ # close the curve
+
+ upper = []
+ lower = []
+
+ indmax = np.argmax(ch1)
+ indmin = np.argmin(ch1)
+ if indmax > indmin:
+ lower = np.array([ch1[indmin:indmax+1],ch2[indmin:indmax+1]]).T
+ upper = np.concatenate((np.array([ch1[indmax:],ch2[indmax:]]).T,\
+ np.array([ch1[:indmin+1],ch2[:indmin+1]]).T),axis=0)
+ else:
+ upper = np.array([ch1[indmax:indmin+1],ch2[indmax:indmin+1]]).T
+ lower = np.concatenate((np.array([ch1[indmin:],ch2[indmin:]]).T,\
+ np.array([ch1[:indmax+1],ch2[:indmax+1]]).T),axis=0)
+
+
+ int_1 = np.linspace(min(min(upper[:,0]),min(lower[:,0])),\
+ max(max(upper[:,0]),max(upper[:,0])),Nx/2+1)
+ upper = np.flipud(upper)
+ int_2_up = np.interp(int_1,np.array(upper[:,0]),np.array(upper[:,1]))
+ int_2_low = np.interp(int_1,np.array(lower[:,0]),np.array(lower[:,1]))
+
+ int_env = np.concatenate((np.array([int_1[:-1],int_2_up[:-1]]).T,\
+ np.array([int_1[::-1],int_2_low[::-1]]).T),axis=0)
+
+ return int_1, int_2_up, int_2_low, int_env
+
+def proj_envelope(env_x,env_up,env_low,env,Nx,Nsectors,Ntheta):
+ # Function to project envelope on given angles
+
+ # Angles of projection is defined by Nsectors
+ # Porjections are obtained in polar coordinates and outputted in
+ # cartesian
+
+ theta_int = np.linspace(-np.pi,np.pi,Ntheta)
+ sectors = np.linspace(-np.pi,np.pi,Nsectors+1)
+ proj = np.zeros([Nsectors,2])
+
+ R_up = np.sqrt(env_x**2+env_up**2)
+ theta_up = np.arctan2(env_up,env_x)
+
+ R_low = np.sqrt(env_x**2+env_low**2)
+ theta_low = np.arctan2(env_low,env_x)
+
+ R = np.concatenate((R_up,R_low))
+ theta = np.concatenate((theta_up,theta_low))
+ R = R[np.argsort(theta)]
+ theta = np.sort(theta)
+
+ R_int = np.interp(theta_int,theta,R,period=2*np.pi)
+
+ for i in range(Nsectors):
+ if sectors[i]>=-np.pi and sectors[i+1]<-np.pi/2:
+ indices = np.where(np.logical_and(theta_int >= sectors[i],\
+ theta_int <= sectors[i+1]))
+ maxR = max(R_int[indices])
+ proj[i+1,0] = maxR*np.cos(sectors[i+1])
+ proj[i+1,1] = maxR*np.sin(sectors[i+1])
+ elif sectors[i]==-np.pi/2:
+ continue
+ elif sectors[i]>-np.pi/2 and sectors[i+1]<=0:
+ indices = np.where(np.logical_and(theta_int >= sectors[i],\
+ theta_int <= sectors[i+1]))
+ maxR = max(R_int[indices])
+ proj[i,0] = maxR*np.cos(sectors[i])
+ proj[i,1] = maxR*np.sin(sectors[i])
+ elif sectors[i]>=0 and sectors[i+1]<np.pi/2:
+ indices = np.where(np.logical_and(theta_int >= sectors[i],\
+ theta_int <= sectors[i+1]))
+ maxR = max(R_int[indices])
+ proj[i+1,0] = maxR*np.cos(sectors[i+1])
+ proj[i+1,1] = maxR*np.sin(sectors[i+1])
+ elif sectors[i]==np.pi/2:
+ continue
+ elif sectors[i]>np.pi/2 and sectors[i+1]<=np.pi:
+ indices = np.where(np.logical_and(theta_int >= sectors[i],\
+ theta_int <= sectors[i+1]))
+ maxR = max(R_int[indices])
+ proj[i,0] = maxR*np.cos(sectors[i])
+ proj[i,1] = maxR*np.sin(sectors[i])
+
+ ind = np.where(sectors==0)
+ proj[ind,0] = max(env[:,0])
+
+ ind = np.where(sectors==np.pi/2)
+ proj[ind,1] = max(env[:,1])
+
+ ind = np.where(sectors==-np.pi)
+ proj[ind,0] = min(env[:,0])
+
+ ind = np.where(sectors==-np.pi/2)
+ proj[ind,1] = min(env[:,1])
+
+
+ proj[-1,:] = proj[0,:]
+
+ return proj
# FIXME: Cases has a memory leek somewhere, this whole thing needs to be
# reconsidered and rely on a DataFrame instead of a dict!
@@ -4891,6 +5066,34 @@ class Cases(object):
return result
def compute_envelope(self, sig, ch_list, int_env=False, Nx=300):
+ """
+ The function computes load envelopes for given signals and a single
+ load case. Starting from Mx and My moments, the other cross-sectional
+ forces are identified.
+
+ Parameters
+ ----------
+
+ sig : list, time-series signal
+
+ ch_list : list, list of channels for enevelope computation
+
+ int_env : boolean, default=False
+ If the logic parameter is True, the function will interpolate the
+ envelope on a given number of points
+
+ Nx : int, default=300
+ Number of points for the envelope interpolation
+
+ Returns
+ -------
+
+ envelope : dictionary,
+ The dictionary has entries refered to the channels selected.
+ Inside the dictonary under each entry there is a matrix with 6
+ columns, each for the sectional forces and moments
+
+ """
envelope= {}
for ch in ch_list:
@@ -4898,16 +5101,23 @@ class Cases(object):
chi1 = self.res.ch_dict[ch[1]]['chi']
s0 = np.array(sig[:, chi0]).reshape(-1, 1)
s1 = np.array(sig[:, chi1]).reshape(-1, 1)
+ # Compute a Convex Hull, the vertices number varies according to
+ # the shape of the poligon
cloud = np.append(s0, s1, axis=1)
hull = scipy.spatial.ConvexHull(cloud)
closed_contour = np.append(cloud[hull.vertices,:],
cloud[hull.vertices[0],:].reshape(1,2),
axis=0)
+
+ # Interpolate envelope for a given number of points
if int_env:
- closed_contour_int = self.int_envelope(closed_contour[:,0],\
- closed_contour[:,1],Nx=Nx)
-
+ _,_,_,closed_contour_int = int_envelope(closed_contour[:,0],
+ closed_contour[:,1],Nx)
+
+ # Based on Mx and My envelope, the other cross-sectional moments
+ # and forces components are identified and appended to the initial
+ # envelope
for ich in range(2, len(ch)):
chix = self.res.ch_dict[ch[ich]]['chi']
s0 = np.array(sig[hull.vertices, chix]).reshape(-1, 1)
@@ -4915,8 +5125,8 @@ class Cases(object):
s0 = np.append(s0, s1, axis=0)
closed_contour = np.append(closed_contour, s0, axis=1)
if int_env:
- extra_sensor = self.int_envelope(closed_contour[:,0],\
- closed_contour[:,ich],Nx=Nx)
+ _,_,_,extra_sensor = int_envelope(closed_contour[:,0],
+ closed_contour[:,ich],Nx)
es = np.atleast_2d(np.array(extra_sensor[:,1])).T
closed_contour_int = np.append(closed_contour_int,es,axis=1)
@@ -4924,40 +5134,9 @@ class Cases(object):
envelope[ch[0]] = closed_contour_int
else:
envelope[ch[0]] = closed_contour
+
return envelope
-
- def int_envelope(ch1,ch2,Nx):
- # Function to interpolate envelopes and output arrays of same length
-
- # Number of points is defined by Nx + 1, where the + 1 is needed to
- # close the curve
-
- upper = []
- lower = []
-
- indmax = np.argmax(ch1)
- indmin = np.argmin(ch1)
- if indmax > indmin:
- lower = np.array([ch1[indmin:indmax+1],ch2[indmin:indmax+1]]).T
- upper = np.concatenate((np.array([ch1[indmax:],ch2[indmax:]]).T,\
- np.array([ch1[:indmin+1],ch2[:indmin+1]]).T),axis=0)
- else:
- upper = np.array([ch1[indmax:indmin+1],ch2[indmax:indmin+1]]).T
- lower = np.concatenate((np.array([ch1[indmin:],ch2[indmin:]]).T,\
- np.array([ch1[:indmax+1],ch2[:indmax+1]]).T),axis=0)
-
-
- int_1 = np.linspace(min(min(upper[:,0]),min(lower[:,0])),\
- max(max(upper[:,0]),max(upper[:,0])),Nx/2+1)
- upper = np.flipud(upper)
- int_2_up = np.interp(int_1,np.array(upper[:,0]),np.array(upper[:,1]))
- int_2_low = np.interp(int_1,np.array(lower[:,0]),np.array(lower[:,1]))
- int_env = np.concatenate((np.array([int_1[:-1],int_2_up[:-1]]).T,\
- np.array([int_1[::-1],int_2_low[::-1]]).T),axis=0)
-
- return int_env
-
def envelope(self, silent=False, ch_list=[], append=''):
"""
Calculate envelopes and save them in a table.
--
GitLab