pyNastran  0.5.0
pyNastran BDF Reader/Writer, OP2 Parser, and GUI
cart3d_reader.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 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     
 All Classes Namespaces Files Functions Variables