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 os 00026 import sys 00027 from math import ceil, sqrt 00028 from itertools import izip 00029 from numpy import array 00030 00031 #from pyNastran.general.general import ListPrint 00032 from struct import unpack 00033 from pyNastran.op2.fortranFile import FortranFile 00034 from pyNastran.general.general import is_binary 00035 00036 def convertToFloat(svalues): 00037 values = [] 00038 for value in svalues: 00039 values.append(float(value)) 00040 return values 00041 00042 #------------------------------------------------------------------ 00043 00044 class Cart3DAsciiReader(object): 00045 modelType = 'cart3d' 00046 isStructured = False 00047 isOutwardNormals = True 00048 00049 def __init__(self,log=None,debug=False): 00050 self.isHalfModel = True 00051 self.cartType = None # grid, result 00052 self.nPoints = None 00053 self.nElements = None 00054 self.infile = None 00055 self.infilename = None 00056 self.readHalf = False 00057 00058 if log is None: 00059 from pyNastran.general.logger import dummyLogger 00060 if debug: 00061 word = 'debug' 00062 else: 00063 word = 'info' 00064 loggerObj = dummyLogger() 00065 log = loggerObj.startLog(word) # or info 00066 self.log = log 00067 00068 def readCart3d(self,infilename): 00069 """extracts the points, elements, and Cp""" 00070 self.infilename = infilename 00071 self.log.info("---starting reading cart3d file...|%s|---" %(self.infilename)) 00072 self.infile = open(self.infilename,'r') 00073 self.readHeader() 00074 points = self.getPoints() 00075 elements = self.getElements(bypass=False) 00076 regions = self.getRegions(bypass=False) 00077 loads = self.readResults(0,self.infile) 00078 #if self.cartType=='results': 00079 00080 #loads = self.getLoads(['Cp']) 00081 self.log.info("nPoints=%s nElements=%s" %(self.nPoints,self.nElements)) 00082 self.log.info("---finished reading cart3d file...|%s|---" %(self.infilename)) 00083 return (points,elements,regions,loads) 00084 00085 def getYZeroNodes(self,nodes): 00086 yZeroNodes = {} 00087 outerNodes = {} 00088 isInnerNode = {} 00089 00090 #nodeTypes = {} 00091 for (iNode,node) in sorted(nodes.items()): 00092 if node[1]==0: 00093 yZeroNodes[iNode] = node 00094 isInnerNode[iNode] = True 00095 else: 00096 outerNodes[iNode] = [node[0],-node[1],node[2]] 00097 isInnerNode[iNode] = False 00098 ### 00099 ### 00100 return (yZeroNodes,outerNodes,isInnerNode) 00101 00102 def makeFullModel(self,nodes,elements,regions,loads): 00103 """assumes a y=0 line of nodes""" 00104 self.log.info('---starting makeFullModel---') 00105 00106 maxNode = len(nodes.keys()) 00107 maxElements = len(elements.keys()) 00108 (yZeroNodes,outerNodes,isInnerNode) = self.getYZeroNodes(nodes) 00109 00110 #nodes2={} 00111 nodes2Map = {} 00112 iNode2=1 00113 for (iNode,node) in sorted(nodes.items()): 00114 if isInnerNode[iNode]: 00115 #pass 00116 #nodes2[iNode] = node # dont need to set this... 00117 nodes2Map[iNode] = iNode 00118 else: 00119 #testElements2[maxNode+iNode2] = outerNodes[iNode] 00120 nodes[maxNode+iNode2] = outerNodes[iNode] 00121 nodes2Map[iNode] = maxNode+iNode2 00122 iNode2+=1 00123 ### 00124 ### 00125 00126 # based on the conversion (nodes2map) object, renumber each element 00127 # stick the result in the location for the new element 00128 #elementsB={} 00129 #regionsB={} 00130 00131 nThree = 147120+10 #len(elements)*2 ## todo why is there +10??? 00132 allSet = set([i+1 for i in xrange(nThree)]) 00133 elementKeysSet = set(elements.keys()) 00134 missingKeys = list(allSet.difference(elementKeysSet)) 00135 missingKeys.sort() 00136 #print "missingKeys = ",missingKeys 00137 #print "nMissing = ",len(missingKeys) 00138 iMissing = 0 00139 for (iElement,element) in sorted(elements.items()): 00140 #(n1,n2,n3) = element 00141 region = regions[iElement] 00142 element2 = [] 00143 00144 for nID in element: 00145 nID2 = nodes2Map[nID] 00146 element2.append(nID2) 00147 element2.reverse() # the normal will be flipped otherwise 00148 00149 eMissing = missingKeys[iMissing] 00150 #print "eMissing = ",eMissing 00151 elements[eMissing] = element2 00152 regions[eMissing] = region 00153 iMissing += 1 00154 ### 00155 00156 for i in xrange(1,len(elements)): 00157 assert i in elements,'i=%s is not in elements...' %(i) 00158 #assert i in regions ,'i=%s is not in regions...' %(i) 00159 self.log.info('---finished makeFullModel---') 00160 sys.stdout.flush() 00161 00162 return (nodes,elements,regions,loads) 00163 00164 def makeHalfModel(self,nodes,elements,regions,Cp): 00165 nNodesStart = len(nodes) 00166 00167 self.log.info('---starting makeHalfModel---') 00168 for (iNode,node) in sorted(nodes.items()): 00169 if node[1]<0: 00170 #print iNode,'...',node 00171 del nodes[iNode] 00172 try: 00173 del Cp[iNode] # assume loads=0 00174 except: 00175 pass 00176 ### 00177 ### 00178 sys.stdout.flush() 00179 00180 for (iElement,element) in elements.items(): 00181 #print "iElement = |%s|" %(iElement),type(iElement) 00182 inNodeList = [True,True,True] 00183 for (i,node) in enumerate(element): 00184 if node not in nodes: 00185 inNodeList[i] = False 00186 ### 00187 if False in inNodeList: 00188 #print iElement,element,inNodeList,type(iElement) 00189 del elements[iElement] 00190 del regions[iElement] 00191 ### 00192 ### 00193 assert len(nodes)!=nNodesStart 00194 self.log.info('---finished makeHalfModel---') 00195 return (nodes,elements,regions,Cp) 00196 00197 def makeMirrorModel(self,nodes,elements,regions,loads): 00198 #nNodesStart = len(nodes) 00199 00200 self.log.info('---starting makeMirrorModel---') 00201 for (iNode,node) in sorted(nodes.items()): 00202 if node[1]<0: 00203 #print iNode,'...',node 00204 del nodes[iNode] 00205 #del Cp[iNode] # assume loads=0 00206 sys.stdout.flush() 00207 00208 for (iElement,element) in elements.items(): 00209 #print "iElement = |%s|" %(iElement),type(iElement) 00210 inNodeList = [True,True,True] 00211 for (i,node) in enumerate(element): 00212 if node not in nodes: 00213 inNodeList[i] = False 00214 ### 00215 if False in inNodeList: 00216 #print iElement,element,inNodeList,type(iElement) 00217 del elements[iElement] 00218 del regions[iElement] 00219 ### 00220 ### 00221 #assert len(nodes)!=nNodesStart 00222 self.log.info('---finished makeMirrorModel---') 00223 00224 return self.makeFullModel(nodes,elements,regions,loads) 00225 #return self.renumberMesh(nodes,elements,regions,loads) 00226 00227 def renumberMesh(self,nodes,elements,regions,loads): 00228 iNodeCounter = 1 00229 iElementCounter = 1 00230 00231 NodeOldToNew = {} 00232 #ElementOldToNew = {} 00233 nodes2 = {} 00234 elements2 = {} 00235 regions2 = {} 00236 for (iNode,node) in sorted(nodes.items()): 00237 NodeOldToNew[iNode] = iNodeCounter 00238 nodes2[iNodeCounter] = node 00239 iNodeCounter += 1 00240 00241 00242 for (iElement,element) in sorted(elements.items()): 00243 element2 = [] 00244 for nID in element: 00245 nID2 = NodeOldToNew[nID] 00246 element2.append(nID2) 00247 elements2[iElementCounter] = element2 00248 regions2[iElementCounter] = regions[iElement] 00249 iElementCounter += 1 00250 00251 return (nodes,elements,regions,loads) 00252 00253 def writeOutfile(self,outfilename,points,elements,regions): 00254 assert len(points)>0 00255 self.log.info("---writing cart3d file...|%s|---" %(outfilename)) 00256 f = open(outfilename,'wb') 00257 00258 self.writeHeader(f,points,elements) 00259 self.writePoints(f,points) 00260 self.writeElements(f,points,elements) 00261 self.writeRegions(f,regions) 00262 f.close() 00263 00264 def writeRegions(self,f,regions): 00265 msg = '' 00266 for (iElement,region) in sorted(regions.items()): 00267 msg += "%s\n" %(region) 00268 f.write(msg) 00269 00270 def writeElements(self,f,nodes,elements): 00271 msg = '' 00272 strElements = len(str(len(elements))) 00273 Format = " %%%ss " %(strElements) 00274 #Format = '%-10s '+Format*3 + '\n' 00275 #Format = Format*3 + '\n' 00276 #print "Format = ",Format 00277 00278 maxNodes = [] 00279 for (iElement,element) in sorted(elements.items()): 00280 #self.log.info("element = %s" %(element)) 00281 #msg += Format %(iElement,element[0],element[1],element[2]) 00282 msg += "%s %s %s\n" %(element[0], element[1], element[2]) 00283 for nid in element: 00284 assert nid in nodes, 'nid=%s is missing' % (nid) 00285 00286 maxNode = max(element) 00287 maxNodes.append(maxNode) 00288 self.log.info("maxNodeID = %s" %(max(maxNodes))) 00289 f.write(msg) 00290 00291 def writePoints(self,f,points): 00292 msg = '' 00293 for (iPoint,point) in sorted(points.items()): 00294 #msg += "%-10s %-15s %-15s %-15s\n" %(iPoint,point[0],point[1],point[2]) 00295 msg += "%6.6f %6.6f %6.6f\n" %(point[0],point[1],point[2]) 00296 if iPoint%1000: 00297 f.write(msg) 00298 msg = '' 00299 f.write(msg) 00300 00301 def writeHeader(self,f,points,elements): 00302 nPoints = len(points) 00303 nElements = len(elements) 00304 msg = "%s %s\n" %(nPoints,nElements) 00305 f.write(msg) 00306 00307 def readHeader(self): 00308 line = self.infile.readline() 00309 sline = line.strip().split() 00310 if len(sline)==2: 00311 self.nPoints,self.nElements = int(sline[0]),int(sline[1]) 00312 self.cartType = 'grid' 00313 elif len(sline)==3: 00314 self.nPoints,self.nElements,self.nResults = int(sline[0]),int(sline[1]),int(sline[2]) 00315 self.cartType = 'result' 00316 else: 00317 raise ValueError('invalid result type') 00318 00319 if self.readHalf: 00320 self.nElementsRead = self.nElements//2 00321 self.nElementsSkip = self.nElements//2 00322 else: 00323 self.nElementsRead = self.nElements 00324 self.nElementsSkip = 0 00325 00326 def getPoints(self): 00327 """ 00328 A point is defined by x,y,z and the ID is the location in points. 00329 """ 00330 points = {} 00331 p=0 00332 data = [] 00333 while p<self.nPoints: 00334 data += self.infile.readline().strip().split() 00335 while len(data)>2: 00336 x = data.pop(0) 00337 y = data.pop(0) 00338 z = data.pop(0) 00339 points[p+1] = array([float(x),float(y),float(z)]) 00340 p+=1 00341 ### 00342 00343 maxX = self.getMax(points,0) 00344 maxY = self.getMax(points,1) 00345 maxZ = self.getMax(points,2) 00346 00347 minX = self.getMin(points,0) 00348 minY = self.getMin(points,1) 00349 minZ = self.getMin(points,2) 00350 00351 self.log.info("X max=%g min=%g" %(maxX,minX)) 00352 self.log.info("Y max=%g min=%g" %(maxY,minY)) 00353 self.log.info("Z max=%g min=%g" %(maxZ,minZ)) 00354 return points 00355 00356 def getMin(self,points,i): 00357 minX = points[1][i] 00358 for key,point in points.items(): 00359 minX = min(minX,point[i]) 00360 return minX 00361 00362 def getMax(self,points,i): 00363 maxX = points[1][i] 00364 for key,point in points.items(): 00365 maxX = max(maxX,point[i]) 00366 return maxX 00367 00368 def getElements(self,bypass=False): 00369 """ 00370 An element is defined by n1,n2,n3 and the ID is the location in elements. 00371 """ 00372 assert bypass==False 00373 elements = {} 00374 00375 print "nElementsRead=%s nElementsSkip=%s" %(self.nElementsRead,self.nElementsSkip) 00376 00377 e=0 00378 data = [] 00379 while e<self.nElementsRead: 00380 data += self.infile.readline().strip().split() 00381 while len(data)>2: 00382 n1 = int(data.pop(0)) 00383 n2 = int(data.pop(0)) 00384 n3 = int(data.pop(0)) 00385 elements[e+1] = [n1,n2,n3] 00386 e+=1 00387 00388 e=0 00389 while e<self.nElementsSkip: 00390 data += self.infile.readline().strip().split() 00391 while len(data)>2: 00392 data.pop() # popping from the end is faster 00393 data.pop() 00394 data.pop() 00395 e+=1 00396 00397 return elements 00398 00399 def getRegions(self,bypass=True): 00400 regions = {} 00401 if bypass: 00402 for i in xrange(self.nElements): 00403 self.infile.readline() 00404 else: 00405 r=0 00406 data = [] 00407 while r<self.nElementsRead: 00408 data += self.infile.readline().strip().split() 00409 while len(data)>0: 00410 regions[r+1] = int(data.pop(0)) 00411 r+=1 00412 00413 r=0 00414 while r<self.nElementsSkip: 00415 data += self.infile.readline().strip().split() 00416 while len(data)>0: 00417 data.pop() 00418 r+=1 00419 ### 00420 return regions 00421 00422 def readResults(self,i,infile): 00423 """ 00424 Reads the Cp results. 00425 Cp = (p - 1/gamma) / (0.5*M_inf*M_inf) 00426 (pV)^2 = (pu)^2+(pv)^2+(pw)^2 00427 M^2 = (pV)^2/rho^2 00428 00429 # ??? 00430 p = (gamma-1)*(e- (rhoU**2+rhoV**2+rhoW**2)/(2.*rho)) 00431 00432 # ??? 00433 rho,rhoU,rhoV,rhoW,rhoE 00434 """ 00435 Cp = {} 00436 Mach = {} 00437 U = {} 00438 V = {} 00439 W = {} 00440 if(self.cartType=='result'): 00441 nResultLines = int(ceil(self.nResults/5.))-1 00442 #print "nResultLines = ",nResultLines 00443 for pointNum in xrange(self.nPoints): 00444 #print "pointNum = ", pointNum 00445 sline = infile.readline().strip().split() # rho rhoU,rhoV,rhoW,pressure/rhoE/E 00446 i+=1 00447 for n in xrange(nResultLines): 00448 #print "nResultLines = ",n 00449 sline += infile.readline().strip().split() # Cp 00450 #infile.readline() 00451 i+=1 00452 #print "sline = ",sline 00453 #gamma = 1.4 00454 #else: 00455 # p=0. 00456 sline = self.getList(sline) 00457 00458 #p=0 00459 cp = sline[0] 00460 rho = float(sline[1]) 00461 if(rho>abs(0.000001)): 00462 rhoU = float(sline[2]) 00463 rhoV = float(sline[3]) 00464 rhoW = float(sline[4]) 00465 #rhoE = float(sline[5]) 00466 VV = (rhoU)**2+(rhoV)**2+(rhoW)**2/rho**2 00467 mach = sqrt(VV) 00468 if mach>10: 00469 print "nid=%s Cp=%s mach=%s rho=%s rhoU=%s rhoV=%s rhoW=%s" %(pointNum,cp,mach,rho,rhoU,rhoV,rhoW) 00470 00471 Cp[pointNum+1] = cp 00472 Mach[pointNum+1] = mach 00473 U[pointNum+1] = rhoU/rho 00474 V[pointNum+1] = rhoV/rho 00475 W[pointNum+1] = rhoW/rho 00476 #print "pt=%s i=%s Cp=%s p=%s" %(pointNum,i,sline[0],p) 00477 ### 00478 ### 00479 loads = {} 00480 if Cp: 00481 loads['Cp'] = Cp 00482 if Mach: 00483 loads['Mach'] = Mach 00484 if U: 00485 loads['U'] = U 00486 loads['V'] = V 00487 loads['W'] = W 00488 return loads 00489 ### 00490 00491 def getList(self,sline): 00492 """Takes a list of strings and converts them to floats.""" 00493 try: 00494 sline2 = convertToFloat(sline) 00495 except ValueError: 00496 print("sline = %s" %(sline)) 00497 raise SyntaxError('cannot parse %s' %(sline)) 00498 return sline2 00499 00500 def getValue(self,sline): 00501 """Converts a string to a float.""" 00502 value = float(sline[1]) 00503 return value 00504 ### 00505 00506 def exportToNastran(self,fname,points,elements,regions): 00507 from pyNastran.bdf.fieldWriter import printCard 00508 00509 f = open(fname, 'wb') 00510 for nid, grid in enumerate(1, points): 00511 (x, y, z) = grid 00512 f.write(printCard(['GRID', nid, '', x, y, z])) 00513 ### 00514 00515 e = 1e7 00516 nu = 0.3 00517 g = '' 00518 thickness = 0.25 00519 setRegions = list(set(regions)) 00520 for pidMid in setRegions: 00521 f.write(printCard(['MAT1',pidMid,e,g,nu])) 00522 f.write(printCard(['PSHELL',pidMid,pidMid,thickness])) 00523 ### 00524 00525 for eid,(nodes,region) in enumerate(1,izip(elements,regions)): 00526 (n1,n2,n3) = nodes 00527 f.write(printCard(['CTRIA3',eid,region,n1,n2,n3])) 00528 ### 00529 f.close() 00530 00531 #------------------------------------------------------------------ 00532 def getSegments(self,nodes,elements): 00533 segments = {} # key=eid, 00534 lengths = {} 00535 for eid,e in elements.iteritems(): 00536 a = tuple(sorted([ e[0],e[1] ])) # segments of e 00537 b = tuple(sorted([ e[1],e[2] ])) 00538 c = tuple(sorted([ e[2],e[0] ])) 00539 print eid,a,b,c 00540 00541 #a.sort() 00542 #b.sort() 00543 #c.sort() 00544 if a in segments: 00545 segments[a].append(eid) 00546 else: 00547 segments[a] = [eid] 00548 #print(a) 00549 #print nodes[1] 00550 lengths[a] = nodes[a[1]]-nodes[a[0]] 00551 00552 if b in segments: 00553 segments[b].append(eid) 00554 else: 00555 segments[b] = [eid] 00556 lengths[b] = nodes[b[1]]-nodes[b[0]] 00557 00558 if c in segments: 00559 segments[c].append(eid) 00560 else: 00561 segments[c] = [eid] 00562 lengths[c] = nodes[c[1]]-nodes[c[0]] 00563 00564 return segments,lengths 00565 00566 def makeQuads(self,nodes,elements): 00567 segments,lengths = self.getSegments(nodes,elements) 00568 for eid,e in elements.iteritems(): 00569 a = tuple(sorted([ e[0],e[1] ])) # segments of e 00570 b = tuple(sorted([ e[1],e[2] ])) 00571 c = tuple(sorted([ e[2],e[0] ])) 00572 #a.sort() 00573 #b.sort() 00574 #c.sort() 00575 print eid,e 00576 print segments[a] 00577 print lengths[a] 00578 print len(segments[a]) 00579 00580 eidA = self.getSegment(a,eid,segments) 00581 eidB = self.getSegment(b,eid,segments) 00582 eidC = self.getSegment(c,eid,segments) 00583 print "eidA=%s eidB=%s eidC=%s" %(eidA,eidB,eidC) 00584 if eidA: 00585 i=0 00586 e2 = elements[eidA] 00587 self.checkQuad(nodes,eid,eidA,e,e2,a,b,c,i) 00588 del segments[a] 00589 if eidB: 00590 i=1 00591 e2 = elements[eidB] 00592 self.checkQuad(nodes,eid,eidB,e,e2,a,b,c,i) 00593 del segments[b] 00594 if eidC: 00595 i=2 00596 e2 = elements[eidC] 00597 self.checkQuad(nodes,eid,eidC,e,e2,a,b,c,i) 00598 del segments[c] 00599 00600 print "------" 00601 #break 00602 #for segment in segments: 00603 asdf 00604 00605 def checkQuad(self,nodes,eid,eidA,e,e2,a,b,c,i): 00606 """ 00607 A----B 00608 | \ e| 00609 |e2 \| 00610 C----D 00611 00612 two tests 00613 1. folding angle A-B x A-C 00614 2a. abs(A-C) - abs(B-D) = 0 (abs to prevent 2L) 00615 2b. abs(A-B) - abs(C-D) = 0 00616 """ 00617 00618 iplus1 = i+1 00619 iplus2 = i+2 00620 00621 if iplus1 > 2: 00622 iplus1 -= 3 00623 if iplus2 > 2: 00624 iplus2 -= 3 00625 print (i,iplus1) 00626 print (iplus1,iplus2) 00627 print (iplus2,i) 00628 AD = nodes[e[i]] -nodes[e[iplus1]] 00629 AB = nodes[e[iplus1]]-nodes[e[iplus2]] 00630 BD = nodes[e[iplus2]]-nodes[e[i]] 00631 00632 print AD 00633 print AB 00634 print BD 00635 print e2 00636 j = e2.index(e[i]) 00637 00638 jplus1 = j+1 00639 jplus2 = j+2 00640 if jplus1 > 2: 00641 jplus1 -= 3 00642 if jplus2 > 2: 00643 jplus2 -= 3 00644 00645 print "DA = ",e[j],e[jplus1] 00646 DA = nodes[e[j]] -nodes[e[jplus1]] 00647 print DA 00648 00649 asdf 00650 00651 def getSegment(self,a,eid,segments): 00652 if a in segments: 00653 aElems = segments[a] 00654 print aElems 00655 i = aElems.index(eid) 00656 #print i 00657 aElems.pop(i) 00658 #print aElems 00659 eidA = aElems[0] 00660 #eidA = elements[a] 00661 print "eidA = ",eidA 00662 return eidA 00663 return None 00664 00665 00666 #------------------------------------------------------------------ 00667 class Cart3DBinaryReader(FortranFile,Cart3DAsciiReader): 00668 modelType = 'cart3d' 00669 isStructured = False 00670 isOutwardNormals = True 00671 00672 def __init__(self,log=None,debug=False): 00673 Cart3DAsciiReader.__init__(self) 00674 FortranFile.__init__(self) 00675 00676 self.isHalfModel = True 00677 self.cartType = None # grid, result 00678 self.nPoints = None 00679 self.nElements = None 00680 self.infile = None 00681 self.infilename = None 00682 self.readHalf = False 00683 self.isNodeArray = True 00684 00685 self.makeOp2Debug = False 00686 self.n = 0 00687 00688 if log is None: 00689 from pyNastran.general.logger import dummyLogger 00690 loggerObj = dummyLogger() 00691 if debug: 00692 log = loggerObj.startLog('debug') # or info 00693 else: 00694 log = loggerObj.startLog('info') # or info 00695 self.log = log 00696 00697 def readCart3d(self,infileName): 00698 self.infilename = infileName 00699 self.op2 = open(infileName,'rb') 00700 (self.nPoints,self.nElements) = self.readHeader() 00701 points = self.readNodes(self.nPoints) 00702 elements = self.readElements(self.nElements) 00703 regions = self.readRegions(self.nElements) 00704 loads = {} 00705 #print "size = ",size 00706 #print self.printSection2(200,'>') 00707 00708 00709 #print self.printSection2(200,'>') 00710 return (points,elements,regions,loads) 00711 00712 def readHeader(self): 00713 data = self.op2.read(4) 00714 size, = unpack('>i',data) 00715 #print "size = ",size 00716 00717 data = self.op2.read(size) 00718 so4 = size//4 # size over 4 00719 if so4==3: 00720 (nPoints,nElements,nResults) = unpack('>iii',data) 00721 print "nPoints=%s nElements=%s nResults=%s" %(nPoints,nElements,nResults) 00722 elif so4==2: 00723 (nPoints,nElements) = unpack('>ii',data) 00724 print "nPoints=%s nElements=%s" %(nPoints,nElements) 00725 else: 00726 raise RuntimeError('in the wrong spot...endian...') 00727 self.op2.read(8) # end of first block, start of second block 00728 00729 return (nPoints,nElements) 00730 00731 def readNodes(self,nPoints): 00732 #print "starting Nodes" 00733 isBuffered = True 00734 size = nPoints*12 # 12=3*4 all the points 00735 Format = '>'+'f'*3000 # 3000 floats; 1000 points 00736 00737 nodes = {} 00738 00739 np=1 00740 while size>12000: # 12k = 4 bytes/float*3 floats/point*1000 points 00741 #print "size = ",size 00742 n = np+999 00743 #print "nStart=%s np=%s" %(n,np) 00744 data = self.op2.read(4*3000) 00745 nodeXYZs = list(unpack(Format,data)) 00746 while nodeXYZs: 00747 z = nodeXYZs.pop() 00748 y = nodeXYZs.pop() 00749 x = nodeXYZs.pop() 00750 node = array([x,y,z]) 00751 assert n not in nodes,'nid=%s in nodes' %(n) 00752 nodes[n] = node 00753 n-=1 00754 np+=1 00755 size -= 4*3000 00756 00757 assert size>=0 00758 #print "size = ",size 00759 #while size>0: # 4k is 1000 points 00760 n = nPoints 00761 if size>0: 00762 #leftover = size-(nPoints-np)*12 00763 data = self.op2.read(size) 00764 #print "leftover=%s size//4=%s" %(leftover,size//4) 00765 Format = '>'+'f'*(size//4) 00766 00767 #print "len(data) = ",len(data) 00768 nodeXYZs = list(unpack(Format, data)) 00769 while nodeXYZs: 00770 z = nodeXYZs.pop() 00771 y = nodeXYZs.pop() 00772 x = nodeXYZs.pop() 00773 node = array([x,y,z]) 00774 assert n not in nodes,'nid=%s in nodes' %(n) 00775 nodes[n] = node 00776 n-=1 00777 np+=1 00778 ### 00779 size=0 00780 00781 #for p,point in sorted(nodes.iteritems()): 00782 #print "%s %s %s" %(tuple(point)) 00783 00784 if isBuffered: 00785 pass 00786 else: 00787 raise RuntimeError('unBuffered') 00788 00789 for nid in xrange(1,nPoints+1): 00790 assert nid in nodes,'nid=%s not in nodes' %(nid) 00791 self.op2.read(8) # end of second block, start of third block 00792 return nodes 00793 00794 def readElements(self,nElements): 00795 self.nElementsRead = nElements 00796 self.nElementsSkip = 0 00797 #print "starting Elements" 00798 isBuffered = True 00799 size = nElements*12 # 12=3*4 all the elements 00800 Format = '>'+'i'*3000 00801 00802 elements = {} 00803 00804 ne=0 00805 while size>12000: # 4k is 1000 elements 00806 #print "size = ",size 00807 e=ne+1000 00808 data = self.op2.read(4*3000) 00809 nodes = list(unpack(Format,data)) 00810 while nodes: 00811 n3 = nodes.pop() 00812 n2 = nodes.pop() 00813 n1 = nodes.pop() 00814 element = [n1,n2,n3] 00815 elements[e] = element 00816 e-=1 00817 ne+=1 00818 size -= 4*3000 00819 00820 assert size>=0 00821 #print "size = ",size 00822 #while size>0: # 4k is 1000 elements 00823 if size>0: 00824 #leftover = size-(nElements-ne)*12 00825 data = self.op2.read(size) 00826 #print "leftover=%s size//4=%s" %(leftover,size//4) 00827 Format = '>'+'i'*(size//4) 00828 00829 #print "len(data) = ",len(data) 00830 nodes = list(unpack(Format,data)) 00831 e=nElements 00832 while nodes: 00833 n3 = nodes.pop() 00834 n2 = nodes.pop() 00835 n1 = nodes.pop() 00836 element = [n1,n2,n3] 00837 elements[e] = element 00838 e-=1 00839 ne+=1 00840 ### 00841 size=0 00842 00843 #for p,point in sorted(nodes.iteritems()): 00844 #print "%s %s %s" %(tuple(point)) 00845 00846 if isBuffered: 00847 pass 00848 else: 00849 raise RuntimeError('unBuffered') 00850 00851 self.op2.read(8) # end of third (element) block, start of regions (fourth) block 00852 #print "finished Elements" 00853 return elements 00854 00855 def readRegions(self,nElements): 00856 #print "starting Regions" 00857 #isBuffered = True 00858 size = nElements*4 # 12=3*4 all the elements 00859 Format = '>'+'i'*3000 00860 00861 regions = {} 00862 00863 nr=0 00864 while size>12000: # 12k is 3000 elements 00865 #print "size = ",size 00866 data = self.op2.read(4*3000) 00867 regionData = list(unpack(Format,data)) 00868 00869 r=nr+3000 00870 #print "len(regionData) = ",len(regionData) 00871 while regionData: 00872 regions[r] = regionData.pop() 00873 r-=1 00874 nr+=1 00875 size -= 4*3000 00876 #print "size = ",size 00877 00878 assert size>=0 00879 #print "size = ",size 00880 00881 if size>0: 00882 #leftover = size-(nElements-nr)*4 00883 data = self.op2.read(size) 00884 #print "leftover=%s size//4=%s" %(leftover,size//4) 00885 Format = '>'+'i'*(size//4) 00886 00887 #print "len(data) = ",len(data) 00888 regionData = list(unpack(Format,data)) 00889 r=nElements 00890 while regionData: 00891 regions[r] = regionData.pop() 00892 r-=1 00893 nr+=1 00894 ### 00895 size=0 00896 00897 00898 self.op2.read(4) # end of regions (fourth) block 00899 #print "finished Regions" 00900 return regions 00901 00902 #------------------------------------------------------------------ 00903 00904 def genericCart3DReader(infileName,log=None,debug=False): 00905 print "infileName = ",infileName 00906 f = open(infileName,'rb') 00907 data = f.read(4) 00908 f.close() 00909 00910 if is_binary(infileName): 00911 00912 #eight, = unpack('>i',data) 00913 #if eight==8: # is binary 00914 #print "Binary eight = ",eight 00915 obj = Cart3DBinaryReader(log,debug) 00916 else: 00917 #print "Ascii eight = ",eight 00918 obj = Cart3DAsciiReader(log,debug) 00919 00920 return obj 00921 00922 00923 00924 00925 if __name__=='__main__': 00926 # binary 00927 cart3dGeom = os.path.join('models','spike.a.tri') 00928 #outfilename = os.path.join('models','threePlugs2.tri') 00929 cart = genericCart3DReader(cart3dGeom) 00930 (points,elements,regions,loads) = cart.readCart3d(cart3dGeom) 00931 #cart.makeQuads(points,elements) 00932 00933 cart.writeOutfile(outfilename,points,elements,regions) 00934 00935 # ascii 00936 cart3dGeom = os.path.join('models','threePlugs.a.tri') 00937 outfilename = os.path.join('models','threePlugs2.a.tri') 00938 cart2 = genericCart3DReader(cart3dGeom) 00939 (points,elements,regions,loads) = cart2.readCart3d(cart3dGeom) 00940 cart2.writeOutfile(outfilename,points,elements,regions) 00941 #print points 00942 00943 if 0: 00944 basepath = os.getcwd() 00945 configpath = os.path.join(basepath,'inputs') 00946 workpath = os.path.join(basepath,'outputsFinal') 00947 #bdfModel = os.path.join(configpath,'aeroModel.bdf') 00948 #assert os.path.exists(bdfModel),'|%s| doesnt exist' %(bdfModel) 00949 os.chdir(workpath) 00950 print "basepath",basepath 00951 00952 00953 cart3dGeom = os.path.join(configpath,'Cart3d_bwb2.i.tri') 00954 cart3dGeom2 = os.path.join(workpath,'Cart3d_half.i.tri') 00955 cart3dGeom3 = os.path.join(workpath,'Cart3d_full.i.tri') 00956 #cart3dGeom4 = os.path.join(workpath,'Cart3d_full2.i.tri') 00957 00958 cart = Cart3DAsciiReader() 00959 (points,elements,regions,loads) = cart.readCart3d(cart3dGeom) 00960 (points,elements,regions,loads) = cart.makeHalfModel(points,elements,regions,loads) 00961 cart.writeOutfile(cart3dGeom2,points,elements,regions) 00962 00963 #cart = Cart3DAsciiReader(cart3dGeom) # bJet.a.tri 00964 #(cartPoints,elements,regions,loads) = cart.read() 00965 #points2 = {} 00966 #for (iPoint,point) in sorted(points.items()): 00967 #(x,y,z) = point 00968 #print "p=%s x=%s y=%s z=%s z2=%s" %(iPoint,x,y,z,z+x/10.) 00969 #points2[iPoint] = [x,y,z+x/10.] 00970 #(points,elements,regions,loads) = cart.makeMirrorModel(points2,elements,regions,loads) 00971 00972 cart2 = Cart3DAsciiReader(cart3dGeom2) 00973 (points,elements,regions,loads) = cart2.read() 00974 00975 #makeFullModel 00976 (points,elements,regions,loads) = cart2.makeFullModel(points,elements,regions,loads) # half model values going in 00977 cart2.writeOutfile(cart3dGeom3,points,elements,regions) 00978 00979 #cart3 = Cart3DAsciiReader(cart3dGeom2) 00980 #(points,elements,regions,loads) = cart3.read() 00981 #(points,elements,regions,loads) = cart3.makeMirrorModel(points,elements,regions,loads) 00982 00983 00984 #print "loads = ",ListPrint(loads),len(loads) 00985 #cartOutfile = os.path.join(workpath,'bJet.a.tri_test') 00986 #cart.writeInfile(cartOutfile,cartPoints,elements,regions) 00987