pyNastran
0.5.0
pyNastran BDF Reader/Writer, OP2 Parser, and GUI
|
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 '