pyNastran  0.5.0
pyNastran BDF Reader/Writer, OP2 Parser, and GUI
panairGridPatch.py
Go to the documentation of this file.
00001 ## GNU Lesser General Public License
00002 ## 
00003 ## Program pyNastran - a python interface to NASTRAN files
00004 ## Copyright (C) 2011-2012  Steven Doyle, Al Danial
00005 ## 
00006 ## Authors and copyright holders of pyNastran
00007 ## Steven Doyle <mesheb82@gmail.com>
00008 ## Al Danial    <al.danial@gmail.com>
00009 ## 
00010 ## This file is part of pyNastran.
00011 ## 
00012 ## pyNastran is free software: you can redistribute it and/or modify
00013 ## it under the terms of the GNU Lesser General Public License as published by
00014 ## the Free Software Foundation, either version 3 of the License, or
00015 ## (at your option) any later version.
00016 ## 
00017 ## pyNastran is distributed in the hope that it will be useful,
00018 ## but WITHOUT ANY WARRANTY; without even the implied warranty of
00019 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
00020 ## GNU General Public License for more details.
00021 ## 
00022 ## You should have received a copy of the GNU Lesser General Public License
00023 ## along with pyNastran.  If not, see <http://www.gnu.org/licenses/>.
00024 ## 
00025 import sys
00026 from numpy import array, zeros, cross, transpose
00027 from numpy.linalg import norm
00028 
00029 def sInt(value):
00030     """
00031     int represented as a short float
00032     """
00033     value = "%f" %(value)
00034     return value.rstrip('0')
00035 
00036 
00037 class PanairPatch(object):
00038     def __init__(self,iNetwork,netName,kt,cpNorm,x,y,z,log):
00039         self.log = log
00040 
00041         self.iNetwork = iNetwork # lets it print out in order, does it need a deepcopy???
00042         self.netName = netName.strip()
00043         self.log.debug('netName=|%s|' %(netName))
00044         self.log.debug("****patch.netName=%s" %(self.netName))
00045         self.kt      = kt
00046         self.cpNorm  = cpNorm
00047         self.matchw=0
00048         self.x = x
00049         self.y = y
00050         self.z = z
00051         shape = x.shape
00052         self.nRows = shape[0]
00053         self.nCols = shape[1]
00054         
00055         self.log.debug("shape = %s" %(str(shape)))
00056 
00057     def process(self):
00058         msg = '     network # being processed %3i\n\n' %(self.iNetwork+1)
00059         return msg
00060 
00061     def quickSummary(self,cumPts,cumPn):
00062         msg = ''
00063         if self.kt==1:
00064             src=1; dblt=12; nlopt1=5; nropt1=3; nlopt2=7; nropt2=-2; ipot=2
00065         elif self.kt==5:
00066             src=1; dblt=12; nlopt1=6; nropt1=9; nlopt2=7; nropt2=-2; ipot=2
00067         elif self.kt==18:
00068             src=0; dblt=18; nlopt1=0; nropt1=9; nlopt2=15; nropt2= 2; ipot=2
00069             if self.matchw==1.:  nlopt2=6
00070             self.log.debug("18...matchw = %s" %(self.matchw))
00071         elif self.kt==20:
00072             src=0; dblt=20; nlopt1=0; nropt1=9; nlopt2=6; nropt2= 2; ipot=2
00073         else:
00074             raise NotImplementedError('new kt...kt=%s' %(self.kt))
00075         pts =self.nPoints()
00076         pans=self.nPanels()
00077         #cumPts=33
00078         #cumPn =50
00079         msg +=' %-10s %4s %7s %7s %3s '   % (self.netName, self.iNetwork+1, self.nRows, self.nCols, self.kt,)
00080         msg += '%4s %5s %7s %7s %7s %7s ' % (src, dblt, nlopt1, nropt1, nlopt2, nropt2)
00081         msg += '%7s %7s %7s %7s '         % (ipot, pts, pans, self.cpNorm)
00082         msg += '%7s %7s\n'                % (cumPts, cumPn)
00083         return msg
00084 
00085     def nPanels(self):
00086         return (self.nRows-1)*(self.nCols-1)
00087 
00088     def nPoints(self):
00089         return (self.nRows)*(self.nCols)
00090 
00091     def getPanelPoints(self,iPanel):
00092         r = iPanel%(self.nRows-1)
00093         c = iPanel/(self.nRows-1)
00094 
00095         #print "r=%s c=%s" %(r,c)
00096         p1 = self.getPoint(r,  c  )
00097         p2 = self.getPoint(r,  c+1)
00098         p3 = self.getPoint(r+1,c+1)
00099         p4 = self.getPoint(r+1,c  )
00100         return (p1,p2,p3,p4)
00101 
00102     def getPanelPointIDs(self,iPanel):
00103         r = iPanel%(self.nRows-1)
00104         c = iPanel/(self.nRows-1)
00105 
00106         #print "r=%s c=%s" %(r,c)
00107         p1 = self.getPointID(r,  c  )
00108         p2 = self.getPointID(r,  c+1)
00109         p3 = self.getPointID(r+1,c+1)
00110         p4 = self.getPointID(r+1,c  )
00111         return (p1,p2,p3,p4)
00112 
00113     def getSubpanelProperties(self, iPanel):
00114         (p1,p2,p3,p4) = self.getPanelPoints(iPanel)
00115         p5 = 0.5*(p1+p2)
00116         p6 = 0.5*(p2+p3)
00117         p7 = 0.5*(p3+p4)
00118         p8 = 0.5*(p4+p1)
00119         p9 = 0.25*(p1+p2+p3+p4)  # centroid
00120         
00121         p10 = 0.5*(p5+p6)
00122         p11 = 0.5*(p6+p7)
00123         p12 = 0.5*(p7+p8)
00124         p13 = 0.5*(p8+p5)
00125         N1 = cross(p10-p12,p11-p13)
00126         n1 = N1/norm(N1)
00127         
00128         N2 = cross(p5-p7,p6-p8)
00129         n2 = N2/norm(N2)
00130 
00131     def getPanelProperties(self,iPanel):
00132         (p1,p2,p3,p4) = self.getPanelPoints(iPanel)
00133         a = p1-p3
00134         b = p2-p4
00135         centroid = 0.25*(p1+p2+p3+p4)
00136 
00137         N = cross(a,b)
00138         normN = norm(N)
00139         n = N/normN  #normal vector
00140 
00141         S = 0.5*normN  # area
00142 
00143         u = ( p1+p2-p3-p4)/2.  # longitudinal
00144         p = (-p1+p2+p3-p4)/2.  # transverse
00145 
00146         u = 0.5*( a+b)  # longitudinal vector in local coordinates
00147         p = 0.5*(-a+b)  # transverse vector in local coordinates
00148         o = cross(n,u)  # normal to both vectors in local coordinates
00149 
00150         diameter = norm(a-b)
00151 
00152         return (S,n,centroid,diameter,u,p,o)
00153 
00154     def getPanelAreaNormal(self,iPanel):
00155         (p1,p2,p3,p4) = self.getPanelPoints(iPanel)
00156         a = p1-p3
00157         b = p2-p4
00158 
00159         N = cross(a,b)
00160         normN = norm(N)
00161         n = N/normN  #normal vector
00162 
00163         S = 0.5*normN
00164         return (S,n)
00165 
00166     def getPanelArea(self,iPanel):
00167         # iPanel=200
00168         (p1,p2,p3,p4) = self.getPanelPoints(iPanel)
00169         
00170         a = p1-p3
00171         b = p2-p4
00172         S = 0.5*norm(cross(a,b))
00173         return S
00174 
00175     def getPoint(self,row,col):
00176         return array([self.x[row][col],
00177                       self.y[row][col],
00178                       self.z[row][col]])
00179 
00180     def getPointID(self,row,col):
00181         return col*self.nRows+row
00182 
00183     def getIPoint(self,iPoint):
00184         iRow = iPoint/(self.nCols)
00185         iCol = iPoint%(self.nCols)
00186         #self.log.debug("iPoint=%s iRow=%s iCol=%s" %(iPoint,iRow,iCol))
00187         return self.getPoint(iRow,iCol)
00188 
00189     def getEdge(self,edgeNumber):
00190         """
00191         gets all the points associated with a given edge
00192         
00193                 edge1
00194               0  1  2   -> i (row)
00195         edge4 3  4  5
00196               6  7  8  edge2
00197               9  10 11
00198             |   edge3
00199             j
00200         """
00201         edgeNumber=2
00202         #edgeNumber = 4
00203         #print "edgeNumber=%s" %(edgeNumber)
00204         if edgeNumber == 1:
00205             x = self.x[0][:]
00206             y = self.y[0][:]
00207             z = self.z[0][:] # pretty sure edge 1 is the 0th row
00208             p = [iCol for iCol in xrange(self.nCols)] # good
00209         elif edgeNumber == 2:
00210             self.log.debug("x.shape = %s" % (str(self.x.shape)))
00211             x = self.x[:][self.nCols-1]
00212             y = self.y[:][self.nCols-1]
00213             z = self.z[:][self.nCols-1] # pretty sure edge 2 is the 0th row
00214             p = [iCol*(self.nRows)+(self.nRows-1) for iCol in xrange(self.nCols)]  #
00215             #p = [iRow*(self.nCols)+(self.nCols-1) for iRow in xrange(self.nRows)]  #
00216         elif edgeNumber == 3:
00217             x = self.x[self.nRows-1][:]
00218             y = self.y[self.nRows-1][:]
00219             z = self.z[self.nRows-1][:] # pretty sure edge3 is the last row
00220             p = [iCol+self.nRows*iCol for iCol in xrange(self.nCols)] # good
00221             #p = [(self.nCols-1)*(self.nRows)+iRow for iRow in xrange(self.nRows)]
00222         elif edgeNumber == 4:
00223             x = self.x[:][0]
00224             y = self.y[:][0]
00225             z = self.z[:][0] # pretty sure edge 2 is the 0th row
00226             p = [self.nRows*iCol for iCol in xrange(self.nCols)] # good
00227         else:
00228             raise ValueError('invalid edge; edgeNumber=%s' % (edgeNumber))
00229         self.log.debug("nRows=%s nCols=%s edgeNumber=%s" % (self.nRows, self.nCols, edgeNumber))
00230         #print "nx = ",len(x)
00231         #print "p = ",p
00232         #print "x = ",x
00233         #print "y = ",y
00234         #print "z = ",z
00235         p = [iPoint for iPoint in xrange(self.nPoints())]
00236         for pointID in p:
00237             #pointID = 2
00238             p2 = self.getIPoint(pointID)
00239             #print "point[%s]=%s" %(pointID,p2)
00240 
00241         return (p,x,y,z)
00242 
00243     def getEdges(self):
00244         nx = 2*(self.nRows+self.nCols)-2
00245         p = zeros(nx)
00246         x = zeros(nx)
00247         y = zeros(nx)
00248         z = zeros(nx)
00249         
00250         i = 0
00251         for edgeID in xrange(1,4+1):
00252             (p1,x1,y1,z1) = self.getEdge(edgeID)
00253             nx1 = len(x1)
00254             p[i:i+nx1] = p1[0:nx1]
00255             x[i:i+nx1] = x1[0:nx1]
00256             y[i:i+nx1] = y1[0:nx1]
00257             z[i:i+nx1] = z1[0:nx1]
00258             self.log.debug("-----")
00259         return (p,x,y,z)
00260             
00261     def getElements(self,pointI):
00262         panels = []
00263         #print "pointI=%s" %(pointI)
00264         for iPanel in xrange(self.nPanels()):
00265             panel = self.getPanelPointIDs(iPanel)
00266             panel2 = []
00267             
00268             for p in panel:
00269                 panel2.append(p+pointI)
00270             #print "panel=%s panel2=%s" %(str(panel),str(panel2))
00271             panels.append(panel2)
00272         return panels
00273     
00274     def getPoints(self):
00275         points = []
00276         #self.log.debug("size(X) = %s" %( str( self.x.shape ) ))
00277         #print "size(X) = %s" %( str(X.size())
00278         
00279         for j in xrange(self.nCols):
00280             for i in xrange(self.nRows):
00281                 point = [self.x[i][j],self.y[i][j],self.z[i][j]]
00282                 points.append(point)
00283 
00284         return points,len(points)
00285 
00286     def writeAsPlot3D(self):
00287         out = ''
00288         x = self.x.ravel() # unravel
00289         y = self.y.ravel() # unravel
00290         z = self.z.ravel() # unravel
00291         
00292         for xi in x:
00293             out += "%s " %(xi)
00294         out += "\n"
00295 
00296         for yi in y:
00297             out += "%s " %(yi)
00298         out += "\n"
00299 
00300         for zi in z:
00301             out += "%s " %(zi)
00302         out += "\n"
00303         self.log.debug(out)
00304         #print x
00305         #for c in xrange(self.nCols):
00306         #    nPointsLeft = nFullLines*2+nPartialLines
00307         #    for r in xrange(0,self.nRows,2)
00308         return out
00309         
00310     def rotate(self):
00311         """
00312         not complete...
00313         """
00314         self.x = transpose(self.x)
00315         self.y = transpose(self.y)
00316         self.z = transpose(self.z)
00317         #self.x[0:n][:] = self.x[-n:-1][:] # something like this...
00318 
00319     def __repr__(self):
00320         """
00321         $points - body to wing wakes
00322         =kn                                               cpnorm
00323         1.
00324         =kt
00325         20.
00326         =nm       nn                                                          netname
00327         4.        2.                                                          awbw
00328         """
00329         #x = self.writeAsPlot3D()
00330 
00331         self.log.debug("*******")
00332         header = '$points - surface panels\n'
00333         points = ''
00334 
00335         header += '%-10s%-10s\n'  %('1.', self.cpNorm) #nNetworks is 1
00336         header += '%-10s\n' %(sInt(self.kt))
00337         header += '%-10s%-10s%50s%-10s\n' %(sInt(self.nRows),sInt(self.nCols),'',self.netName)
00338 
00339         #nFullLines = nm/2
00340         #nPartialLines = nm%2
00341         #nLines = nFullLines+nPartialLines
00342 
00343         nFullLines = self.nRows/2
00344         nPartialLines = self.nRows%2
00345         nLines = nFullLines+nPartialLines
00346 
00347         for c in xrange(self.nCols):
00348             nPointsLeft = nFullLines*2+nPartialLines
00349             for r in xrange(0,self.nRows,2):
00350                 if nPointsLeft>1:
00351                     x1 = self.x[r][c]
00352                     y1 = self.y[r][c]
00353                     z1 = self.z[r][c]
00354 
00355                     x2 = self.x[r+1][c]
00356                     y2 = self.y[r+1][c]
00357                     z2 = self.z[r+1][c]
00358                     points += self.writePoints([x1,y1,z1],[x2,y2,z2])
00359                 else:
00360                     x1 = self.x[r][c]
00361                     y1 = self.y[r][c]
00362                     z1 = self.z[r][c]
00363                     points += self.writePoint([x1,y1,z1])
00364                 ###
00365                 nPointsLeft -=2
00366             ###
00367         ###
00368         return header+points
00369 
00370     def writePoints(self,point1,point2):
00371         point1 = self.fixPoint(point1)
00372         point2 = self.fixPoint(point2)
00373         
00374         out = "%-10s"*6 %(point1[0],point1[1],point1[2],  point2[0],point2[1],point2[2])
00375         return out+'\n'
00376 
00377     def writePoint(self,point1):
00378         point1 = self.fixPoint(point1)
00379         out = "%-10s"*3 %(point1[0],point1[1],point1[2])
00380         return out+'\n'
00381     
00382     def fixPoint(self,pointIn):
00383         pointOut = []
00384         for value in pointIn:
00385             sValue = '%s' %(value)
00386             if len(sValue)>10:
00387                 sValue = sValue[0:9]
00388             pointOut.append(sValue.rstrip('0'))
00389             #print "sValue=%s len=%s" %(sValue,len(sValue))
00390         #print "pointOut = ",pointOut
00391         return pointOut
00392 
00393 class PanairWakePatch(PanairPatch):
00394     def __init__(self,iNetwork,netName,options,x,y,z,log):
00395         (kt,cpNorm,matchw,trailedPanel,edgeNumber,xWake,tWake) = options
00396         PanairPatch.__init__(self,iNetwork,netName,kt,cpNorm,x,y,z,log)
00397 
00398         self.log = log
00399         self.matchw = matchw
00400         self.trailedPanel = trailedPanel
00401         self.edgeNumber = edgeNumber
00402         self.xWake = xWake
00403         self.tWake = tWake
00404         self.log.debug("matchw = %s" %(self.matchw))
00405         self.log.debug("wake patch")
00406 
00407     def __repr__(self):
00408         header = '$trailing wakes\n'
00409         points = ''
00410 
00411         header += '%-10s%-10s\n'  %(1., self.cpNorm) #nNetworks is 1
00412         header += '%-10s%-10s\n' %(sInt(self.kt),self.matchw)
00413         #header += '%-10s%-10s%40s%-10s\n' %(self.nRows,self.nCols,'',self.netName)
00414         header += '%-10s%-10s%-10s%-10s%-30s%-10s\n'%(self.trailedPanel,
00415                                                       sInt(self.edgeNumber),
00416                                                       self.xWake,
00417                                                       sInt(self.tWake),
00418                                                       ' ',
00419                                                       self.netName)
00420         return header
00421 
00422 class PanairGridHelper(object):
00423 
00424     def getCases(self,section):
00425         """
00426         $cases - no. of solutions
00427         =nacase
00428         1.
00429         """
00430         self.ncases = int(float(section[1][0:10]))
00431         self.log.debug("ncases = %s" %(self.ncases))
00432         return True
00433 
00434     def getMach(self,section):
00435         """
00436         $mach number
00437         =amach
00438         .6
00439         """
00440         self.mach = float(section[1][0:10])
00441         self.log.debug("mach = %s" %(self.mach))
00442         return True
00443 
00444     def setMach(self,mach):
00445         self.mach = mach
00446 
00447     def writeMach(self):
00448         out  = '$mach number\n'
00449         out += '%-10s' %(self.mach) +'\n'
00450         return out
00451 
00452     def writeCases(self):
00453         out  = '$cases - number of solutions\n'
00454         out += '%-10s' %(sInt(self.ncases)) +'\n'
00455         return out
00456 
00457     def getAlphas(self,section):
00458         """
00459         $angles-of-attack
00460         =alpc
00461         4.
00462         =alpha(1) alpha(2)  alpha(3)
00463         4.        10.       0.
00464         """
00465         self.alphas = []
00466         self.alphaC = float(section[1][0:10]) # alphaCompressibility
00467         sline = section[2].split()
00468         self.alphas = [float(slot) for slot in sline]
00469         self.log.debug("alphaC=%s alphas=%s" %(self.alphaC,self.alphas))
00470         return True
00471 
00472     def setAlphas(self,alphas,alphaC):
00473         self.alphaC = alphaC
00474         self.alphas = alphas
00475         self.ncases = len(alphas)
00476 
00477     def writeAlphas(self):
00478         out  = '$angles-of-attack\n'
00479         out += '%-10s' %(self.alphaC) +'\n'
00480         out += '%-10s'*len(self.alphas) %(tuple(self.alphas)) +'\n'
00481         return out
00482 
00483     def getBetas(self,section):
00484         """
00485         $angles-of-attack
00486         =alpc
00487         4.
00488         =alpha(1) alpha(2)  alpha(3)
00489         4.        10.       0.
00490         """
00491         self.betas = []
00492         self.betaC = float(section[1][0:10]) # betaCompressibility
00493         sline = section[2].split()
00494         self.betas = [float(slot) for slot in sline]
00495         self.log.debug("betaC=%s betas=%s" %(self.betaC,self.betas))
00496         return True
00497 
00498     def setBetas(self,betas,betaC):
00499         self.betaC = betaC
00500         self.betas = betas
00501         self.ncases = len(betas)
00502 
00503     def writeBetas(self):
00504         out  = '$yaw\n'
00505         out += '%-10s' %(self.betaC) +'\n'
00506         out += '%-10s'*len(self.betas) %(tuple(self.betas)) +'\n'
00507         return out
00508 
00509     def getReferenceQuantities(self,section):
00510         """
00511         $references for accumulated forces and moments
00512         =xref     yref      zref      nref
00513         46.       0.        0.
00514         =sref     bref      cref      dref
00515         2400.     60.       40.       90.
00516         """
00517         self.xref = float(section[1][0 :10])  # 0
00518         self.yref = float(section[1][10:20])
00519         self.zref = float(section[1][20:30])
00520        #self.nref = float(section[1][30:40])
00521         self.sref = float(section[2][0 :10])  # 0
00522         self.bref = float(section[2][10:20])
00523         self.cref = float(section[2][20:30])
00524         self.dref = float(section[2][30:40])
00525         self.log.debug("xref=%s yref=%s zref=%s" %(self.xref,self.yref,self.zref))
00526         self.log.debug("sref=%s bref=%s cref=%s dref=%s " %(self.sref,self.bref,self.cref,self.dref))
00527         return True
00528 
00529     def writeReferenceQuantities(self):
00530         out  = '$references for accumulated forces and moments\n'
00531         out += '%-10s'*3 %(self.xref,self.yref,self.zref) +'\n'
00532         out += '%-10s'*4 %(self.sref,self.bref,self.cref,self.dref)+'\n'
00533         return out
00534 
00535     def getEnd(self,section):
00536         self.isEnd = True
00537         self.log.debug("end...")
00538         return True
00539 
00540     def writeEnd(self):
00541         if self.isEnd:
00542             return '$end of panair inputs\n '
 All Classes Namespaces Files Functions Variables