pyNastran  0.5.0
pyNastran BDF Reader/Writer, OP2 Parser, and GUI
oes_solids.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 
00026 from __future__ import division, print_function
00027 import sys
00028 from numpy import array
00029 from numpy.linalg import eigh
00030 
00031 from .oes_objects import stressObject,strainObject
00032 
00033 
00034 class SolidStressObject(stressObject):
00035     """
00036     # sCode=0
00037             N O N L I N E A R   S T R E S S E S   I N   T E T R A H E D R O N   S O L I D   E L E M E N T S   ( T E T R A )
00038     ELEMENT GRID/   POINT                         STRESSES/ TOTAL STRAINS                          EQUIVALENT EFF. STRAIN  EFF. CREEP
00039        ID   GAUSS     ID       X           Y           Z           XY          YZ          ZX        STRESS   PLAS/NLELAS   STRAIN
00040         3    GRID   CENTER  6.6667E+02  2.6667E+03  2.6667E+03 -1.3333E+03  2.6667E+03 -1.3333E+03  6.0000E+03  1.5000E-04   0.0
00041                             1.6667E-05  6.6667E-05  6.6667E-05 -6.6667E-05  1.3333E-04 -6.6667E-05
00042 
00043     # sCode=1
00044                           S T R E S S E S   I N   H E X A H E D R O N   S O L I D   E L E M E N T S   ( H E X A )
00045                    CORNER        ------CENTER AND CORNER POINT STRESSES---------       DIR.  COSINES       MEAN                   
00046     ELEMENT-ID    GRID-ID        NORMAL              SHEAR             PRINCIPAL       -A-  -B-  -C-     PRESSURE       VON MISES 
00047         1               0GRID CS  8 GP
00048                    CENTER  X   4.499200E+02  XY  -5.544791E+02   A   1.000000E+04  LX 0.00 0.69-0.72  -3.619779E+03    9.618462E+03
00049                            Y   4.094179E+02  YZ   5.456968E-12   B  -1.251798E+02  LY 0.00 0.72 0.69
00050                            Z   1.000000E+04  ZX  -4.547474E-13   C   9.845177E+02  LZ 1.00 0.00 0.00
00051     """
00052     def __init__(self,dataCode,isSort1,iSubcase,dt=None):
00053         stressObject.__init__(self,dataCode,iSubcase)
00054 
00055         self.eType = {}
00056         self.code = [self.formatCode,self.sortCode,self.sCode]
00057 
00058         self.cid = {} # gridGauss
00059         self.oxx = {}
00060         self.oyy = {}
00061         self.ozz = {}
00062         self.txy = {}
00063         self.tyz = {}
00064         self.txz = {}
00065         self.o1 = {}
00066         self.o2 = {}
00067         self.o3 = {}
00068         self.ovmShear = {}
00069 
00070         self.dt = dt
00071         if isSort1:
00072             if dt is not None:
00073                 self.add = self.addSort1
00074                 self.addNewEid = self.addNewEidSort1
00075             ###
00076         else:
00077             assert dt is not None
00078             self.add = self.addSort2
00079             self.addNewEid = self.addNewEidSort2
00080         ###
00081 
00082     def addF06Data(self,data,transient):
00083         if transient is None:
00084             if not hasattr(self,'data'):
00085                 self.data = []
00086             self.data += data
00087         else:
00088             if not hasattr(self,'data'):
00089                 self.data = {}
00090             if dt not in self.data:
00091                 self.data[dt] = []
00092             for line in data:
00093                 self.data[dt] += data
00094             ###
00095         ###
00096 
00097     def processF06Data(self):
00098         """
00099                           S T R E S S E S   I N   P E N T A H E D R O N   S O L I D   E L E M E N T S   ( P E N T A )
00100                        CORNER        ------CENTER AND CORNER POINT STRESSES---------       DIR.  COSINES       MEAN                   
00101         ELEMENT-ID    GRID-ID        NORMAL              SHEAR             PRINCIPAL       -A-  -B-  -C-     PRESSURE       VON MISES 
00102                 2           0GRID CS  6 GP
00103                        CENTER  X  -1.829319E+03  XY   7.883865E+01   A   1.033182E+04  LX-0.12 0.71 0.70  -2.115135E+03    1.232595E+04
00104                                Y  -1.825509E+03  YZ  -1.415218E+03   B  -2.080181E+03  LY-0.12 0.69-0.71
00105                                Z   1.000023E+04  ZX  -1.415218E+03   C  -1.906232E+03  LZ 0.99 0.16 0.00
00106         """
00107         eMap = {'CTETRA':5,'CPENTA':7,'CHEXA':9,'HEXA':9,'PENTA':7,'TETRA':5,}   # +1 for the centroid
00108         if self.nonlinearFactor is None:
00109             ipack = []
00110             i = 0
00111             n = 0
00112             while n<len(self.data):
00113                 line = self.data[n]
00114                 #print n,line
00115 
00116                 eType = line[0]
00117                 eid   = int(line[1])
00118                 #print "eType = ",eType
00119                 nNodes = eMap[eType]
00120                 #nodeID = None
00121                 self.eType[eid] = eType
00122                 self.oxx[eid] = {}
00123                 self.oyy[eid] = {}
00124                 self.ozz[eid] = {}
00125                 self.txy[eid] = {}
00126                 self.tyz[eid] = {}
00127                 self.txz[eid] = {}
00128                 self.o1[eid]  = {}
00129                 self.o2[eid]  = {}
00130                 self.o3[eid]  = {}
00131                 self.ovmShear[eid] = {}
00132                 n += 1
00133                 for j in xrange(nNodes):
00134                     #print self.data[n]
00135                     (blank,nodeID,x,oxx,xy,txy,a,o1,lx,d1,d2,d3,pressure,ovmShear) = self.data[n]
00136                     (blank,blank, y,oyy,yz,tyz,b,o2,ly,d1,d2,d3,blank,blank) = self.data[n+1]
00137                     (blank,blank, z,ozz,zx,txz,c,o3,lz,d1,d2,d3,blank,blank) = self.data[n+2]
00138                     if    nodeID.strip()=='CENTER': nodeID='C'
00139                     else: nodeID = int(nodeID)
00140 
00141                     self.oxx[eid][nodeID] = float(oxx)
00142                     self.oyy[eid][nodeID] = float(oyy)
00143                     self.ozz[eid][nodeID] = float(ozz)
00144                     self.txy[eid][nodeID] = float(txy)
00145                     self.tyz[eid][nodeID] = float(tyz)
00146                     self.txz[eid][nodeID] = float(txz)
00147                     self.o1[eid][nodeID] = float(o1)
00148                     self.o2[eid][nodeID] = float(o2)
00149                     self.o3[eid][nodeID] = float(o3)
00150                     self.ovmShear[eid][nodeID] = float(ovmShear)
00151                     n+=3
00152                 ###
00153             ###
00154             return
00155         raise NotImplementedError()
00156 
00157         (dtName,dt) = transient
00158         self.dataCode['name'] = dtName
00159         if dt not in self.gridTypes:
00160             self.addNewTransient()
00161 
00162         for line in data:
00163             (nodeID,gridType,t1,t2,t3,r1,r2,r3) = line
00164             self.gridTypes[dt][nodeID]    = array([t1,t2,t3])
00165             self.translations[dt][nodeID] = array([t1,t2,t3])
00166             self.rotations[dt][nodeID]    = array([r1,r2,r3])
00167         ###
00168         del self.data
00169 
00170     def deleteTransient(self,dt):
00171         del self.oxx[dt]
00172         del self.oyy[dt]
00173         del self.ozz[dt]
00174         del self.txy[dt]
00175         del self.tyz[dt]
00176         del self.txz[dt]
00177         del self.o1[dt]
00178         del self.o2[dt]
00179         del self.o3[dt]
00180 
00181     def getTransients(self):
00182         k = self.oxx.keys()
00183         k.sort()
00184         return k
00185 
00186     def addNewTransient(self,dt):
00187         """
00188         initializes the transient variables
00189         """
00190         self.oxx[dt] = {}
00191         self.oyy[dt] = {}
00192         self.ozz[dt] = {}
00193         self.txy[dt] = {}
00194         self.tyz[dt] = {}
00195         self.txz[dt] = {}
00196         self.o1[dt] = {}
00197         self.o2[dt] = {}
00198         self.o3[dt] = {}
00199 
00200         #self.aCos[dt] = {}
00201         #self.bCos[dt] = {}
00202         #self.cCos[dt] = {}
00203         #self.pressure[dt] = {}
00204         self.ovmShear[dt]  = {}
00205 
00206     def addNewEid(self,eType,cid,dt,eid,nodeID,oxx,oyy,ozz,txy,tyz,txz,o1,o2,o3,aCos,bCos,cCos,pressure,ovm):
00207         #print "Solid Stress add..."
00208         #assert eid not in self.oxx
00209         assert cid >= 0
00210         assert eid >= 0
00211         self.eType[eid] = eType
00212         self.cid[eid]  = cid
00213         self.oxx[eid]  = {nodeID: oxx}
00214         self.oyy[eid]  = {nodeID: oyy}
00215         self.ozz[eid]  = {nodeID: ozz}
00216         self.txy[eid]  = {nodeID: txy}
00217         self.tyz[eid]  = {nodeID: tyz}
00218         self.txz[eid]  = {nodeID: txz}
00219         self.o1[eid]   = {nodeID: o1}
00220         self.o2[eid]   = {nodeID: o2}
00221         self.o3[eid]   = {nodeID: o3}
00222         #self.aCos[eid] = {nodeID: aCos}
00223         #self.bCos[eid] = {nodeID: bCos}
00224         #self.cCos[eid] = {nodeID: cCos}
00225         #self.pressure[eid] = {nodeID: pressure}
00226         self.ovmShear[eid]  = {nodeID: ovm}
00227         msg = "*eid=%s nodeID=%s vm=%g" %(eid,nodeID,ovm)
00228         #print msg
00229         if nodeID==0: raise Exception(msg)
00230 
00231     def addNewEidSort1(self,eType,cid,dt,eid,nodeID,oxx,oyy,ozz,txy,tyz,txz,o1,o2,o3,aCos,bCos,cCos,pressure,ovm):
00232         #print "Solid Stress add transient..."
00233         assert cid >= 0
00234         assert eid >= 0
00235 
00236         #print "dt=%s eid=%s eType=%s" %(dt,eid,eType)
00237         if dt not in self.oxx:
00238             self.addNewTransient(dt)
00239         assert eid not in self.oxx[dt],self.oxx[dt]
00240 
00241         self.eType[eid] = eType
00242         self.cid[eid]   = cid
00243         self.oxx[dt][eid] = {nodeID: oxx}
00244         self.oyy[dt][eid] = {nodeID: oyy}
00245         self.ozz[dt][eid] = {nodeID: ozz}
00246         self.txy[dt][eid] = {nodeID: txy}
00247         self.tyz[dt][eid] = {nodeID: tyz}
00248         self.txz[dt][eid] = {nodeID: txz}
00249         
00250         self.o1[dt][eid]  = {nodeID: o1}
00251         self.o2[dt][eid]  = {nodeID: o2}
00252         self.o3[dt][eid]  = {nodeID: o3}
00253 
00254         #self.aCos[dt][eid] = {nodeID: aCos}
00255         #self.bCos[dt][eid] = {nodeID: bCos}
00256         #self.cCos[dt][eid] = {nodeID: cCos}
00257         #self.pressure[dt][eid] = {nodeID: pressure}
00258         self.ovmShear[dt][eid]  = {nodeID: ovm}
00259         msg = "*eid=%s nodeID=%s vm=%g" %(eid,nodeID,ovm)
00260         #print msg
00261         if nodeID==0: raise Exception(msg)
00262 
00263     def add(self,dt,eid,nodeID,oxx,oyy,ozz,txy,tyz,txz,o1,o2,o3,aCos,bCos,cCos,pressure,ovm):
00264         #print "***add"
00265         msg = "eid=%s nodeID=%s vm=%g" %(eid,nodeID,ovm)
00266         #print msg
00267         #print self.oxx
00268         #print self.fiberDistance
00269         self.oxx[eid][nodeID] = oxx
00270         self.oyy[eid][nodeID] = oyy
00271         self.ozz[eid][nodeID] = ozz
00272 
00273         self.txy[eid][nodeID] = txy
00274         self.tyz[eid][nodeID] = tyz
00275         self.txz[eid][nodeID] = txz
00276 
00277         self.o1[eid][nodeID] = o1
00278         self.o2[eid][nodeID] = o2
00279         self.o3[eid][nodeID] = o3
00280 
00281         #self.aCos[eid][nodeID] = aCos
00282         #self.bCos[eid][nodeID] = bCos
00283         #self.cCos[eid][nodeID] = cCos
00284         #self.pressure[eid][nodeID] = pressure
00285         self.ovmShear[eid][nodeID] = ovm
00286         if nodeID==0: raise Exception(msg)
00287 
00288     def addSort1(self,dt,eid,nodeID,oxx,oyy,ozz,txy,tyz,txz,o1,o2,o3,aCos,bCos,cCos,pressure,ovm):
00289         #print "***add"
00290         msg = "eid=%s nodeID=%s vm=%g" %(eid,nodeID,ovm)
00291         #print msg
00292         #print self.oxx
00293         #print self.fiberDistance
00294 
00295         #self.eType[eid] = eType
00296         #print "eid=%s nid=%s oxx=%s" %(eid,nodeID,oxx)
00297         self.oxx[dt][eid][nodeID] = oxx
00298         self.oyy[dt][eid][nodeID] = oyy
00299         self.ozz[dt][eid][nodeID] = ozz
00300         #print self.oxx
00301         self.txy[dt][eid][nodeID] = txy
00302         self.tyz[dt][eid][nodeID] = tyz
00303         self.txz[dt][eid][nodeID] = txz
00304 
00305         self.o1[dt][eid][nodeID] = o1
00306         self.o2[dt][eid][nodeID] = o2
00307         self.o3[dt][eid][nodeID] = o3
00308 
00309         #self.aCos[dt][eid][nodeID] = aCos
00310         #self.bCos[dt][eid][nodeID] = bCos
00311         #self.cCos[dt][eid][nodeID] = cCos
00312         #self.pressure[dt][eid][nodeID] = pressure
00313         self.ovmShear[dt][eid][nodeID] = ovm
00314         if nodeID==0: raise Exception(msg)
00315 
00316     def getHeaders(self):
00317         headers = ['oxx','oyy','ozz','txy','tyz','txz']
00318         if self.isVonMises():
00319             headers.append('oVonMises')
00320         else:
00321             headers.append('oMaxShear')
00322         return headers
00323 
00324     def pressure(self,o1,o2,o3):
00325         """
00326         returns the hydrostatic pressure
00327         (o1+o2+o3)/-3.
00328         http://en.wikipedia.org/wiki/Stress_%28mechanics%29
00329         """
00330         return (o1+o2+o3)/-3.
00331 
00332     def directionalVectors(self,oxx,oyy,ozz,txy,tyz,txz):
00333         A = [[oxx,txy,txz],
00334              [txy,oyy,tyz],
00335              [txz,tyz,ozz]]
00336         (Lambda,v) = eigh(A) # a hermitian matrix is a symmetric-real matrix
00337         return v
00338     
00339     def getF06_Header(self):
00340         if self.isVonMises():
00341             vonMises = 'VON MISES'
00342         else:
00343             vonMises = 'MAX SHEAR'
00344 
00345         tetraMsg = ['                   S T R E S S E S   I N    T E T R A H E D R O N   S O L I D   E L E M E N T S   ( C T E T R A )\n',
00346                     '0                CORNER        ------CENTER AND CORNER POINT STRESSES---------       DIR.  COSINES       MEAN                   \n',
00347                     '  ELEMENT-ID    GRID-ID        NORMAL              SHEAR             PRINCIPAL       -A-  -B-  -C-     PRESSURE       %s \n' %(vonMises)]
00348 
00349         pentaMsg = ['                    S T R E S S E S   I N   P E N T A H E D R O N   S O L I D   E L E M E N T S   ( P E N T A )\n',
00350                     '0                CORNER        ------CENTER AND CORNER POINT STRESSES---------       DIR.  COSINES       MEAN                   \n',
00351                     '  ELEMENT-ID    GRID-ID        NORMAL              SHEAR             PRINCIPAL       -A-  -B-  -C-     PRESSURE       %s \n' %(vonMises)]
00352 
00353         hexaMsg = ['                      S T R E S S E S   I N   H E X A H E D R O N   S O L I D   E L E M E N T S   ( H E X A )\n',
00354                    '0                CORNER        ------CENTER AND CORNER POINT STRESSES---------       DIR.  COSINES       MEAN                   \n',
00355                    '  ELEMENT-ID    GRID-ID        NORMAL              SHEAR             PRINCIPAL       -A-  -B-  -C-     PRESSURE       %s \n' %(vonMises)]
00356 
00357         #nNodes = {'CTETRA':4,'CPENTA':6,'CHEXA':8,'HEXA':8,'PENTA':6,'TETRA':4,}
00358         tetraEids = []
00359         hexaEids  = []
00360         pentaEids = []
00361         for eid,eType in sorted(self.eType.iteritems()):
00362             if eType=='CTETRA' or eType=='TETRA':
00363                 tetraEids.append(eid)
00364             elif eType=='CPENTA' or eType=='PENTA':
00365                 pentaEids.append(eid)
00366             elif eType=='CHEXA' or eType=='HEXA':
00367                 hexaEids.append(eid)
00368             else:
00369                 raise NotImplementedError('eType=|%r|' %(eType))
00370             ###
00371         ###
00372         return (tetraMsg,pentaMsg,hexaMsg,tetraEids,hexaEids,pentaEids)
00373         
00374     def writeF06(self,header,pageStamp,pageNum=1,f=None,isMagPhase=False):
00375         if self.nonlinearFactor is not None:
00376             return self.writeF06Transient(header,pageStamp,pageNum,f)
00377         msg = []
00378 
00379         (tetraMsg,pentaMsg,hexaMsg,tetraEids,hexaEids,pentaEids) = self.getF06_Header()
00380         #nNodes = {'CTETRA':4,'CPENTA':6,'CHEXA':8,'HEXA':8,'PENTA':6,'TETRA':4,}
00381         if tetraEids:
00382             msg += self.writeElement('CTETRA',4,tetraEids,header,tetraMsg,f)
00383             msg.append(pageStamp+str(pageNum)+'\n')
00384             pageNum+=1
00385         if hexaEids:
00386             msg += self.writeElement('CHEXA',8,hexaEids, header,hexaMsg,f)
00387             msg.append(pageStamp+str(pageNum)+'\n')
00388             pageNum+=1
00389         if pentaEids:
00390             msg += self.writeElement('CPENTA',6,pentaEids,header,pentaMsg,f)
00391             msg.append(pageStamp+str(pageNum)+'\n')
00392             pageNum+=1
00393         ###
00394         return (''.join(msg),pageNum-1)
00395 
00396     def writeF06Transient(self,header,pageStamp,pageNum=1,f=None,isMagPhase=False):
00397         msg = []
00398         (tetraMsg,pentaMsg,hexaMsg,tetraEids,hexaEids,pentaEids) = self.getF06_Header()
00399         dts = self.oxx.keys()
00400         for dt in dts:
00401             if tetraEids:
00402                 msg += self.writeElementTransient('CTETRA',4,tetraEids,dt,header,tetraMsg,f)
00403                 msg.append(pageStamp+str(pageNum)+'\n')
00404                 pageNum+=1
00405             if hexaEids:
00406                 msg += self.writeElementTransient('CHEXA',8,hexaEids, dt,header,hexaMsg,f)
00407                 msg.append(pageStamp+str(pageNum)+'\n')
00408                 pageNum+=1
00409             if pentaEids:
00410                 msg += self.writeElementTransient('CPENTA',6,pentaEids,dt,header,pentaMsg,f)
00411                 msg.append(pageStamp+str(pageNum)+'\n')
00412                 pageNum+=1
00413             ###
00414         return (''.join(msg),pageNum-1)
00415 
00416     def writeElement(self,eType,nNodes,eids,header,tetraMsg,f=None):
00417         msg = header+tetraMsg
00418         for eid in eids:
00419             #eType = self.eType[eid]
00420 
00421             k = self.oxx[eid].keys()
00422             k.remove('C')
00423             k.sort()
00424             msg.append('0  %8s           0GRID CS  %i GP\n' %(eid,nNodes) )
00425             for nid in ['C']+k:
00426                 oxx = self.oxx[eid][nid]
00427                 oyy = self.oyy[eid][nid]
00428                 ozz = self.ozz[eid][nid]
00429                 txy = self.txy[eid][nid]
00430                 tyz = self.tyz[eid][nid]
00431                 txz = self.txz[eid][nid]
00432 
00433                 o1 = self.o1[eid][nid]
00434                 o2 = self.o2[eid][nid]
00435                 o3 = self.o3[eid][nid]
00436                 ovm = self.ovmShear[eid][nid]
00437                 p = (o1+o2+o3)/-3.
00438 
00439                 if nid=='C': nid='CENTER'
00440                 #print "nid = |%r|" %(nid)
00441                 A = [[oxx,txy,txz],
00442                      [txy,oyy,tyz],
00443                      [txz,tyz,ozz]]
00444                 (Lambda,v) = eigh(A) # a hermitian matrix is a symmetric-real matrix
00445 
00446                 ([oxx,oyy,ozz,txy,tyz,txz,o1,o2,o3,p,ovm],isAllZeros) = self.writeFloats13E([oxx,oyy,ozz,txy,tyz,txz,o1,o2,o3,p,ovm])
00447                 msg.append('0              %8s  X  %13s  XY  %13s   A  %13s  LX%5.2f%5.2f%5.2f  %13s   %-s\n' %(nid,oxx,txy,o1,v[0,1],v[0,2],v[0,0],p,ovm.strip()))
00448                 msg.append('               %8s  Y  %13s  YZ  %13s   B  %13s  LY%5.2f%5.2f%5.2f\n'             %('', oyy,tyz,o2,v[1,1],v[1,2],v[1,0]))
00449                 msg.append('               %8s  Z  %13s  ZX  %13s   C  %13s  LZ%5.2f%5.2f%5.2f\n'             %('', ozz,txz,o3,v[2,1],v[2,2],v[2,0]))
00450             ###
00451         ###
00452         if f is not None:
00453             f.write(''.join(msg))
00454         ###
00455         return ''.join(msg)
00456 
00457     def writeElementTransient(self,eType,nNodes,eids,dt,header,tetraMsg,f=None):
00458         dtLine = '%14s = %12.5E\n'%(self.dataCode['name'],dt)
00459         header[1] = dtLine
00460         msg = header+tetraMsg
00461         for eid in eids:
00462             #eType = self.eType[eid]
00463 
00464             k = self.oxx[dt][eid].keys()
00465             k.remove('C')
00466             k.sort()
00467             msg.append('0  %8s           0GRID CS  %i GP\n' %(eid,nNodes))
00468             for nid in ['C']+k:
00469                 oxx = self.oxx[dt][eid][nid]
00470                 oyy = self.oyy[dt][eid][nid]
00471                 ozz = self.ozz[dt][eid][nid]
00472                 txy = self.txy[dt][eid][nid]
00473                 tyz = self.tyz[dt][eid][nid]
00474                 txz = self.txz[dt][eid][nid]
00475 
00476                 o1 = self.o1[dt][eid][nid]
00477                 o2 = self.o2[dt][eid][nid]
00478                 o3 = self.o3[dt][eid][nid]
00479                 ovm = self.ovmShear[dt][eid][nid]
00480                 p = (o1+o2+o3)/-3.
00481 
00482                 if nid=='C': nid='CENTER'
00483                 #print "nid = |%r|" %(nid)
00484                 A = [[oxx,txy,txz],
00485                      [txy,oyy,tyz],
00486                      [txz,tyz,ozz]]
00487                 (Lambda,v) = eigh(A) # a hermitian matrix is a symmetric-real matrix
00488 
00489                 ([oxx,oyy,ozz,txy,tyz,txz,o1,o2,o3,p,ovm],isAllZeros) = self.writeFloats13E([oxx,oyy,ozz,txy,tyz,txz,o1,o2,o3,p,ovm])
00490                 msg.append('0              %8s  X  %13s  XY  %13s   A  %13s  LX%5.2f%5.2f%5.2f  %13s   %-s\n' %(nid,oxx,txy,o1,v[0,1],v[0,2],v[0,0],p,ovm.strip()))
00491                 msg.append('               %8s  Y  %13s  YZ  %13s   B  %13s  LY%5.2f%5.2f%5.2f\n'             %('', oyy,tyz,o2,v[1,1],v[1,2],v[1,0]))
00492                 msg.append('               %8s  Z  %13s  ZX  %13s   C  %13s  LZ%5.2f%5.2f%5.2f\n'             %('', ozz,txz,o3,v[2,1],v[2,2],v[2,0]))
00493             ###
00494         ###
00495         if f is not None:
00496             f.write(''.join(msg))
00497             msg = ['']
00498         ###
00499         return ''.join(msg)
00500 
00501     def __repr__(self):
00502         #print "self.dt = ",self.dt
00503         #print "self.nf = ",self.nonlinearFactor
00504         if self.nonlinearFactor is not None:
00505             return self.__reprTransient__()
00506 
00507         msg = '---SOLID STRESS---\n'
00508         headers = self.getHeaders()
00509         msg += '%-6s %6s %8s ' %('EID','eType','nodeID')
00510         for header in headers:
00511             msg += '%9s ' %(header)
00512         msg += '\n'
00513         for eid,oxxNodes in sorted(self.oxx.iteritems()):
00514             eType = self.eType[eid]
00515             for nid in sorted(oxxNodes):
00516                 oxx = self.oxx[eid][nid]
00517                 oyy = self.oyy[eid][nid]
00518                 ozz = self.ozz[eid][nid]
00519                 txy = self.txy[eid][nid]
00520                 tyz = self.tyz[eid][nid]
00521                 txz = self.txz[eid][nid]
00522 
00523                 #o1 = self.o1[eid][nid]
00524                 #o2 = self.o2[eid][nid]
00525                 #o3 = self.o3[eid][nid]
00526                 ovm = self.ovmShear[eid][nid]
00527                 msg += '%-6i %6s %8s ' %(eid,eType,nid)
00528                 vals = [oxx,oyy,ozz,txy,tyz,txz,ovm]
00529                 for val in vals:
00530                     if abs(val)<1e-6:
00531                         msg += '%9s ' %('0')
00532                     else:
00533                         msg += '%9i ' %(val)
00534                     ###
00535                 msg += '\n'
00536                 #msg += "eid=%-4s eType=%-6s nid=%-2i oxx=%-5i oyy=%-5i ozz=%-5i txy=%-5i tyz=%-5i txz=%-5i ovm=%-5i\n" %(eid,eType,nid,oxx,oyy,ozz,txy,tyz,txz,ovm)
00537             ###
00538         ###
00539         return msg
00540 
00541     def __reprTransient__(self):
00542         msg = '---SOLID STRESS---\n'
00543         msg += '%-6s %6s %8s ' %('EID','eType','nodeID')
00544         headers = self.getHeaders()
00545         for header in headers:
00546             msg += '%9s ' %(header)
00547         msg += '\n'
00548 
00549         for dt,oxxs in sorted(self.oxx.iteritems()):
00550             msg += '%s = %g\n' %(self.dataCode['name'],dt)
00551             for eid,oxxNodes in sorted(oxxs.iteritems()):
00552                 eType = self.eType[eid]
00553                 for nid in sorted(oxxNodes):
00554                     oxx = self.oxx[dt][eid][nid]
00555                     oyy = self.oyy[dt][eid][nid]
00556                     ozz = self.ozz[dt][eid][nid]
00557                     txy = self.txy[dt][eid][nid]
00558                     tyz = self.tyz[dt][eid][nid]
00559                     txz = self.txz[dt][eid][nid]
00560 
00561                     #o1 = self.o1[dt][eid][nid]
00562                     #o2 = self.o2[dt][eid][nid]
00563                     #o3 = self.o3[dt][eid][nid]
00564 
00565                     ovm = self.ovmShear[dt][eid][nid]
00566                     msg += '%-6i %6s %8s ' %(eid,eType,nid)
00567                     vals = [oxx,oyy,ozz,txy,tyz,txz,ovm]
00568                     for val in vals:
00569                         if abs(val)<1e-6:
00570                             msg += '%9s ' %('0')
00571                         else:
00572                             msg += '%9i ' %(val)
00573                         ###
00574                     msg += '\n'
00575                     #msg += "eid=%-4s eType=%-6s nid=%-2i oxx=%-5i oyy=%-5i ozz=%-5i txy=%-5i tyz=%-5i txz=%-5i ovm=%-5i\n" %(eid,eType,nid,oxx,oyy,ozz,txy,tyz,txz,ovm)
00576                 ###
00577             ###
00578         ###
00579         return msg
00580 
00581 class SolidStrainObject(strainObject):
00582     """
00583     # code=[1,0,11]
00584                           S T R A I N S   I N   H E X A H E D R O N   S O L I D   E L E M E N T S   ( H E X A )
00585                    CORNER        ------CENTER AND CORNER POINT STRESSES---------       DIR.  COSINES       MEAN                   
00586     ELEMENT-ID    GRID-ID        NORMAL              SHEAR             PRINCIPAL       -A-  -B-  -C-     PRESSURE       VON MISES 
00587             1           0GRID CS  8 GP
00588                    CENTER  X   4.499200E+02  XY  -5.544791E+02   A   1.000000E+04  LX 0.00 0.69-0.72  -3.619779E+03    9.618462E+03
00589                            Y   4.094179E+02  YZ   5.456968E-12   B  -1.251798E+02  LY 0.00 0.72 0.69
00590                            Z   1.000000E+04  ZX  -4.547474E-13   C   9.845177E+02  LZ 1.00 0.00 0.00
00591 
00592     # code=[1,0,10]
00593                        S T R A I N S   I N    T E T R A H E D R O N   S O L I D   E L E M E N T S   ( C T E T R A )
00594                    CORNER        ------CENTER AND CORNER POINT  STRAINS---------       DIR.  COSINES       MEAN         OCTAHEDRAL
00595     ELEMENT-ID    GRID-ID        NORMAL              SHEAR             PRINCIPAL       -A-  -B-  -C-     PRESSURE         SHEAR   
00596             4           0GRID CS  4 GP
00597                    CENTER  X  -2.288232E-04  XY   1.240506E-04   A   9.631978E-04  LX-0.10-0.71-0.70  -1.601805E-04    5.692614E-04
00598                            Y  -2.289814E-04  YZ  -2.369997E-04   B  -2.909276E-04  LY-0.10 0.71-0.70
00599                            Z   9.383460E-04  ZX  -2.369997E-04   C  -1.917288E-04  LZ 0.99 0.00-0.15
00600     """
00601     def __init__(self,dataCode,isSort1,iSubcase,dt=None):
00602         strainObject.__init__(self,dataCode,iSubcase)
00603         self.eType = {}
00604 
00605         self.code = [self.formatCode,self.sortCode,self.sCode]
00606 
00607         self.cid = {}
00608         self.exx = {}
00609         self.eyy = {}
00610         self.ezz = {}
00611         self.exy = {}
00612         self.eyz = {}
00613         self.exz = {}
00614         self.e1 = {}
00615         self.e2 = {}
00616         self.e3 = {}
00617         #self.aCos = {}
00618         #self.bCos = {}
00619         #self.cCos = {}
00620         #self.pressure = {}
00621         self.evmShear = {}
00622         
00623 
00624         self.nonlinearFactor = dt
00625         if isSort1:
00626             if dt is not None:
00627                 self.add = self.addSort1
00628                 self.addNewEid = self.addNewEidSort1
00629             ###
00630         else:
00631             assert dt is not None
00632             self.add = self.addSort2
00633             self.addNewEid = self.addNewEidSort2
00634         ###
00635 
00636     def addF06Data(self,data,transient):
00637         if transient is None:
00638             if not hasattr(self,'data'):
00639                 self.data = []
00640             self.data += data
00641         else:
00642             if not hasattr(self,'data'):
00643                 self.data = {}
00644             if dt not in self.data:
00645                 self.data[dt] = []
00646             for line in data:
00647                 self.data[dt] += data
00648             ###
00649         ###
00650 
00651     def processF06Data(self):
00652         """
00653                           S T R E S S E S   I N   P E N T A H E D R O N   S O L I D   E L E M E N T S   ( P E N T A )
00654                        CORNER        ------CENTER AND CORNER POINT STRESSES---------       DIR.  COSINES       MEAN                   
00655         ELEMENT-ID    GRID-ID        NORMAL              SHEAR             PRINCIPAL       -A-  -B-  -C-     PRESSURE       VON MISES 
00656                 2           0GRID CS  6 GP
00657                        CENTER  X  -1.829319E+03  XY   7.883865E+01   A   1.033182E+04  LX-0.12 0.71 0.70  -2.115135E+03    1.232595E+04
00658                                Y  -1.825509E+03  YZ  -1.415218E+03   B  -2.080181E+03  LY-0.12 0.69-0.71
00659                                Z   1.000023E+04  ZX  -1.415218E+03   C  -1.906232E+03  LZ 0.99 0.16 0.00
00660         """
00661         eMap = {'CTETRA':5,'CPENTA':7,'CHEXA':9,'HEXA':9,'PENTA':7,'TETRA':5,}   # +1 for the centroid
00662         if self.nonlinearFactor is None:
00663             pack = []
00664             i = 0
00665             n = 0
00666             while n<len(self.data):
00667                 line = self.data[n]
00668                 #print n,line
00669 
00670                 eType = line[0]
00671                 eid   = int(line[1])
00672                 #print "eType = ",eType
00673                 nNodes = eMap[eType]
00674                 #nodeID = None
00675                 self.eType[eid] = eType
00676                 self.exx[eid] = {}
00677                 self.eyy[eid] = {}
00678                 self.ezz[eid] = {}
00679                 self.exy[eid] = {}
00680                 self.eyz[eid] = {}
00681                 self.exz[eid] = {}
00682                 self.e1[eid]  = {}
00683                 self.e2[eid]  = {}
00684                 self.e3[eid]  = {}
00685                 self.evmShear[eid] = {}
00686                 n += 1
00687                 for j in xrange(nNodes):
00688                     #print self.data[n]
00689                     (blank,nodeID,x,exx,xy,exy,a,e1,lx,d1,d2,d3,pressure,evmShear) = self.data[n]
00690                     (blank,blank, y,eyy,yz,eyz,b,e2,ly,d1,d2,d3,blank,blank) = self.data[n+1]
00691                     (blank,blank, z,ezz,zx,exz,c,e3,lz,d1,d2,d3,blank,blank) = self.data[n+2]
00692                     if    nodeID.strip()=='CENTER': nodeID='C'
00693                     else: nodeID = int(nodeID)
00694                     self.exx[eid][nodeID] = float(exx)
00695                     self.eyy[eid][nodeID] = float(eyy)
00696                     self.ezz[eid][nodeID] = float(ezz)
00697                     self.exy[eid][nodeID] = float(exy)
00698                     self.eyz[eid][nodeID] = float(eyz)
00699                     self.exz[eid][nodeID] = float(exz)
00700                     self.e1[eid][nodeID] = float(e1)
00701                     self.e2[eid][nodeID] = float(e2)
00702                     self.e3[eid][nodeID] = float(e3)
00703                     self.evmShear[eid][nodeID] = float(evmShear)
00704                     n+=3
00705                 ###
00706             ###
00707             return
00708         raise NotImplementedError()
00709         (dtName,dt) = transient
00710         self.dataCode['name'] = dtName
00711         if dt not in self.gridTypes:
00712             self.addNewTransient()
00713 
00714         for line in data:
00715             (nodeID,gridType,t1,t2,t3,r1,r2,r3) = line
00716             self.gridTypes[dt][nodeID]    = array([t1,t2,t3])
00717             self.translations[dt][nodeID] = array([t1,t2,t3])
00718             self.rotations[dt][nodeID]    = array([r1,r2,r3])
00719         ###
00720         del self.data
00721 
00722     def deleteTransient(self,dt):
00723         del self.exx[dt]
00724         del self.eyy[dt]
00725         del self.ezz[dt]
00726         del self.exy[dt]
00727         del self.eyz[dt]
00728         del self.exz[dt]
00729         del self.e1[dt]
00730         del self.e2[dt]
00731         del self.e3[dt]
00732 
00733     def getTransients(self):
00734         k = self.exx.keys()
00735         k.sort()
00736         return k
00737 
00738     def addNewTransient(self,dt):
00739         """
00740         initializes the transient variables
00741         """
00742         self.nonlinearFactor = dt
00743         self.exx[dt] = {}
00744         self.eyy[dt] = {}
00745         self.ezz[dt] = {}
00746         self.exy[dt] = {}
00747         self.eyz[dt] = {}
00748         self.exz[dt] = {}
00749         self.e1[dt] = {}
00750         self.e2[dt] = {}
00751         self.e3[dt] = {}
00752 
00753         #self.aCos[dt] = {}
00754         #self.bCos[dt] = {}
00755         #self.cCos[dt] = {}
00756         #self.pressure[dt] = {}
00757         self.evmShear[dt]  = {}
00758         ###
00759 
00760     def addNewEid(self,eType,cid,dt,eid,nodeID,exx,eyy,ezz,exy,eyz,exz,e1,e2,e3,aCos,bCos,cCos,pressure,evm):
00761         #print "Solid Strain add..."
00762         assert eid not in self.exx
00763         assert cid >= 0
00764         assert eid >= 0
00765         self.eType[eid] = eType
00766         self.cid[eid]  = cid
00767         self.exx[eid]  = {nodeID: exx}
00768         self.eyy[eid]  = {nodeID: eyy}
00769         self.ezz[eid]  = {nodeID: ezz}
00770         self.exy[eid]  = {nodeID: exy}
00771         self.eyz[eid]  = {nodeID: eyz}
00772         self.exz[eid]  = {nodeID: exz}
00773         self.e1[eid]   = {nodeID: e1}
00774         self.e2[eid]   = {nodeID: e2}
00775         self.e3[eid]   = {nodeID: e3}
00776         #self.aCos[eid] = {nodeID: aCos}
00777         #self.bCos[eid] = {nodeID: bCos}
00778         #self.cCos[eid] = {nodeID: cCos}
00779         #self.pressure[eid] = {nodeID: pressure}
00780         self.evmShear[eid]      = {nodeID: evm}
00781         msg = "*eid=%s nodeID=%s evmShear=%g" %(eid,nodeID,evm)
00782         #print msg
00783         if nodeID==0: raise Exception(msg)
00784 
00785     def addNewEidSort1(self,eType,cid,dt,eid,nodeID,exx,eyy,ezz,exy,eyz,exz,e1,e2,e3,aCos,bCos,cCos,pressure,evm):
00786         #print "Solid Strain add..."
00787         assert cid >= 0
00788         assert eid >= 0
00789         
00790         if dt not in self.exx:
00791             self.addNewTransient(dt)
00792         
00793         assert eid not in self.exx[dt],self.exx[dt]
00794         self.eType[eid] = eType
00795         self.cid[eid]  = cid
00796         self.exx[dt][eid]  = {nodeID: exx}
00797         self.eyy[dt][eid]  = {nodeID: eyy}
00798         self.ezz[dt][eid]  = {nodeID: ezz}
00799         self.exy[dt][eid]  = {nodeID: exy}
00800         self.eyz[dt][eid]  = {nodeID: eyz}
00801         self.exz[dt][eid]  = {nodeID: exz}
00802         self.e1[dt][eid]   = {nodeID: e1}
00803         self.e2[dt][eid]   = {nodeID: e2}
00804         self.e3[dt][eid]   = {nodeID: e3}
00805         #self.aCos[dt][eid] = {nodeID: aCos}
00806         #self.bCos[dt][eid] = {nodeID: bCos}
00807         #self.cCos[dt][eid] = {nodeID: cCos}
00808         #self.pressure[dt][eid] = {nodeID: pressure}
00809         self.evmShear[dt][eid]      = {nodeID: evm}
00810         msg = "*eid=%s nodeID=%s evmShear=%g" %(eid,nodeID,evm)
00811         #print msg
00812         if nodeID==0: raise Exception(msg)
00813 
00814     def add(self,dt,eid,nodeID,exx,eyy,ezz,exy,eyz,exz,e1,e2,e3,aCos,bCos,cCos,pressure,evm):
00815         #print "***add"
00816         msg = "eid=%s nodeID=%s vm=%g" %(eid,nodeID,evm)
00817         #print msg
00818         #print self.exx
00819         #print self.fiberDistance
00820         self.exx[eid][nodeID] = exx
00821         self.eyy[eid][nodeID] = eyy
00822         self.ezz[eid][nodeID] = ezz
00823 
00824         self.exy[eid][nodeID] = exy
00825         self.eyz[eid][nodeID] = eyz
00826         self.exz[eid][nodeID] = exz
00827         self.e1[eid][nodeID] = e1
00828         self.e2[eid][nodeID] = e2
00829         self.e3[eid][nodeID] = e3
00830 
00831         #self.aCos[eid][nodeID] = aCos
00832         #self.bCos[eid][nodeID] = bCos
00833         #self.cCos[eid][nodeID] = cCos
00834         #self.pressure[eid][nodeID] = pressure
00835         self.evmShear[eid][nodeID] = evm
00836 
00837         if nodeID==0: raise Exception(msg)
00838 
00839     def addSort1(self,dt,eid,nodeID,exx,eyy,ezz,exy,eyz,exz,e1,e2,e3,aCos,bCos,cCos,pressure,evm):
00840         #print "***add"
00841         msg = "eid=%s nodeID=%s vm=%g" %(eid,nodeID,evm)
00842         #print msg
00843         #print self.exx
00844         #print self.fiberDistance
00845         
00846         self.exx[dt][eid][nodeID] = exx
00847         self.eyy[dt][eid][nodeID] = eyy
00848         self.ezz[dt][eid][nodeID] = ezz
00849 
00850         self.exy[dt][eid][nodeID] = exy
00851         self.eyz[dt][eid][nodeID] = eyz
00852         self.exz[dt][eid][nodeID] = exz
00853 
00854         self.e1[dt][eid][nodeID] = e1
00855         self.e2[dt][eid][nodeID] = e2
00856         self.e3[dt][eid][nodeID] = e3
00857 
00858         #self.aCos[dt][eid][nodeID] = aCos
00859         #self.bCos[dt][eid][nodeID] = bCos
00860         #self.cCos[dt][eid][nodeID] = cCos
00861         #self.pressure[dt][eid][nodeID] = pressure
00862         self.evmShear[dt][eid][nodeID] = evm
00863 
00864         if nodeID==0: raise Exception(msg)
00865 
00866     def getHeaders(self):
00867         headers = ['exx','eyy','ezz','exy','eyz','exz']
00868         if self.isVonMises():
00869             headers.append('evm')
00870         else:
00871             headers.append('maxShear')
00872         return headers
00873 
00874     def ovm(self,o11,o22,o33,o12,o13,o23):
00875         """http://en.wikipedia.org/wiki/Von_Mises_yield_criterion"""
00876         ovm = 0.5 * ( (o11-o22)**2+(o22-o33)**2+(o11-o33)**2+6*(o23**2+o13**2+o12**2))
00877         return ovm
00878 
00879     #def ovmPlane(self,o11,o22,o12):
00880         #"""http://en.wikipedia.org/wiki/Von_Mises_yield_criterion"""
00881         #ovm = (o11**2+o22**2-o1*o2+3*o12**2)**0.5
00882         #return ovm
00883 
00884     def octahedral(self,o11,o22,o33,o12,o13,o23):
00885         """http://en.wikipedia.org/wiki/Von_Mises_yield_criterion"""
00886         ovm = self.ovm(o11,o22,o33,o12,o13,o23)
00887         return ovm/3*2**0.5
00888 
00889     def pressure(self,e1,e2,e3):
00890         """
00891         returns the hydrostatic pressure
00892         (e1+e2+e3)/-3.
00893         http://en.wikipedia.org/wiki/Stress_%28mechanics%29
00894         """
00895         return (e1+e2+e3)/-3.
00896 
00897     def directionalVectors(self,exx,eyy,ezz,exy,eyz,exz):
00898         A = [[exx,exy,exz],
00899              [exy,eyy,eyz],
00900              [exz,eyz,ezz]]
00901         (Lambda,v) = eigh(A) # a hermitian matrix is a symmetric-real matrix
00902         return v
00903 
00904     def getF06_Header(self):
00905         if self.isVonMises():
00906             vonMises = 'VON MISES'
00907         else:
00908             vonMises = 'MAX SHEAR'
00909 
00910         tetraMsg = ['                     S T R A I N S   I N    T E T R A H E D R O N   S O L I D   E L E M E N T S   ( C T E T R A )\n',
00911                    '0                CORNER        ------CENTER AND CORNER POINT STRAINS---------       DIR.  COSINES       MEAN\n',
00912                     '  ELEMENT-ID    GRID-ID        NORMAL              SHEAR             PRINCIPAL       -A-  -B-  -C-     PRESSURE       %s \n' %(vonMises)]
00913 
00914         pentaMsg = ['                      S T R A I N S   I N   P E N T A H E D R O N   S O L I D   E L E M E N T S   ( P E N T A )\n',
00915                    '0                CORNER        ------CENTER AND CORNER POINT STRAINS---------       DIR.  COSINES       MEAN\n',
00916                     '  ELEMENT-ID    GRID-ID        NORMAL              SHEAR             PRINCIPAL       -A-  -B-  -C-     PRESSURE       %s \n' %(vonMises)]
00917 
00918         hexaMsg  = ['                      S T R A I N S   I N   H E X A H E D R O N   S O L I D   E L E M E N T S   ( H E X A )\n',
00919                    '0                CORNER        ------CENTER AND CORNER POINT STRAINS---------       DIR.  COSINES       MEAN\n',
00920                    '  ELEMENT-ID    GRID-ID        NORMAL              SHEAR             PRINCIPAL       -A-  -B-  -C-     PRESSURE       %s \n' %(vonMises)]
00921 
00922 
00923         #nNodes = {'CTETRA':4,'CPENTA':6,'CHEXA':8,'HEXA':8,'PENTA':6,'TETRA':4,}
00924         tetraEids = []
00925         hexaEids  = []
00926         pentaEids = []
00927         for eid,eType in sorted(self.eType.iteritems()):
00928             if eType=='CTETRA' or eType=='TETRA':
00929                 tetraEids.append(eid)
00930             elif eType=='CPENTA' or eType=='PENTA':
00931                 pentaEids.append(eid)
00932             elif eType=='CHEXA' or eType=='HEXA':
00933                 hexaEids.append(eid)
00934             else:
00935                 raise NotImplementedError('eType=|%r|' %(eType))
00936             ###
00937         ###
00938         return (tetraMsg,pentaMsg,hexaMsg,tetraEids,hexaEids,pentaEids)
00939 
00940     def writeF06(self,header,pageStamp,pageNum=1,f=None,isMagPhase=False):
00941         if self.nonlinearFactor is not None:
00942             return self.writeF06Transient(header,pageStamp,pageNum,f)
00943         msg = []
00944 
00945         (tetraMsg,pentaMsg,hexaMsg,tetraEids,hexaEids,pentaEids) = self.getF06_Header()
00946         #nNodes = {'CTETRA':4,'CPENTA':6,'CHEXA':8,'HEXA':8,'PENTA':6,'TETRA':4,}
00947         if tetraEids:
00948             msg += self.writeElement('CTETRA',4,tetraEids,header,tetraMsg,f)
00949             msg.append(pageStamp+str(pageNum)+'\n')
00950             pageNum+=1
00951         if hexaEids:
00952             msg += self.writeElement('CHEXA',8,hexaEids,header,hexaMsg,f)
00953             msg.append(pageStamp+str(pageNum)+'\n')
00954             pageNum+=1
00955         if pentaEids:
00956             msg += self.writeElement('CPENTA',6,pentaEids,header,pentaMsg,f)
00957             msg.append(pageStamp+str(pageNum)+'\n')
00958             pageNum+=1
00959         ###
00960         return (''.join(msg),pageNum-1)
00961 
00962     def writeF06Transient(self,header,pageStamp,pageNum=1,f=None,isMagPhase=False):
00963         msg = []
00964         (tetraMsg,pentaMsg,hexaMsg,tetraEids,hexaEids,pentaEids) = self.getF06_Header()
00965         dts = self.exx.keys()
00966         for dt in dts:
00967             if tetraEids:
00968                 msg += self.writeElementTransient('CTETRA',4,tetraEids,dt,header,tetraMsg,f)
00969                 msg.append(pageStamp+str(pageNum)+'\n')
00970                 pageNum+=1
00971             if hexaEids:
00972                 msg += self.writeElementTransient('CHEXA',8,hexaEids, dt,header,hexaMsg,f)
00973                 msg.append(pageStamp+str(pageNum)+'\n')
00974                 pageNum+=1
00975             if pentaEids:
00976                 msg += self.writeElementTransient('CPENTA',6,pentaEids,dt,header,pentaMsg,f)
00977                 msg.append(pageStamp+str(pageNum)+'\n')
00978                 pageNum+=1
00979             ###
00980         return (''.join(msg),pageNum-1)
00981 
00982     def writeElement(self,eType,nNodes,eids,header,tetraMsg,f=None):
00983         msg = header+tetraMsg
00984         for eid in eids:
00985             #eType = self.eType[eid]
00986 
00987             k = self.exx[eid].keys()
00988             k.remove('C')
00989             k.sort()
00990             msg.append('0  %8s           0GRID CS  %i GP\n' %(eid,nNodes) )
00991             for nid in ['C']+k:
00992                 exx = self.exx[eid][nid]
00993                 eyy = self.eyy[eid][nid]
00994                 ezz = self.ezz[eid][nid]
00995                 exy = self.exy[eid][nid]
00996                 eyz = self.eyz[eid][nid]
00997                 exz = self.exz[eid][nid]
00998 
00999                 e1 = self.e1[eid][nid]
01000                 e2 = self.e2[eid][nid]
01001                 e3 = self.e3[eid][nid]
01002                 evm = self.evmShear[eid][nid]
01003                 p = (e1+e2+e3)/-3.
01004 
01005                 if nid=='C': nid='CENTER'
01006                 #print "nid = |%r|" %(nid)
01007                 A = [[exx,exy,exz],
01008                      [exy,eyy,eyz],
01009                      [exz,eyz,ezz]]
01010                 (Lambda,v) = eigh(A) # a hermitian matrix is a symmetric-real matrix
01011 
01012                 ([exx,eyy,ezz,exy,eyz,exz,e1,e2,e3,p,evm],isAllZeros) = self.writeFloats13E([exx,eyy,ezz,exy,eyz,exz,e1,e2,e3,p,evm])
01013                 msg.append('0              %8s  X  %13s  XY  %13s   A  %13s  LX%5.2f%5.2f%5.2f  %13s   %-s\n' %(nid,exx,exy,e1,v[0,1],v[0,2],v[0,0],p,evm.strip()))
01014                 msg.append('               %8s  Y  %13s  YZ  %13s   B  %13s  LY%5.2f%5.2f%5.2f\n'             %('', eyy,eyz,e2,v[1,1],v[1,2],v[1,0]))
01015                 msg.append('               %8s  Z  %13s  ZX  %13s   C  %13s  LZ%5.2f%5.2f%5.2f\n'             %('', ezz,exz,e3,v[2,1],v[2,2],v[2,0]))
01016             ###
01017         ###
01018         if f is not None:
01019             f.write(''.join(msg))
01020             msg = ['']
01021         ###
01022         return ''.join(msg)
01023 
01024     def writeElementTransient(self,eType,nNodes,eids,dt,header,tetraMsg,pageStamp=1,f=None):
01025         dtLine = '%14s = %12.5E\n'%(self.dataCode['name'],dt)
01026         header[1] = dtLine
01027         msg = header+tetraMsg
01028         for eid in eids:
01029             #eType = self.eType[eid]
01030 
01031             k = self.exx[dt][eid].keys()
01032             k.remove('C')
01033             k.sort()
01034             msgA  = '0  %8s           0GRID CS  %i GP\n' %(eid,nNodes)
01035             for nid in ['C']+k:
01036                 exx = self.exx[dt][eid][nid]
01037                 eyy = self.eyy[dt][eid][nid]
01038                 ezz = self.ezz[dt][eid][nid]
01039                 exy = self.exy[dt][eid][nid]
01040                 eyz = self.eyz[dt][eid][nid]
01041                 exz = self.exz[dt][eid][nid]
01042 
01043                 e1 = self.e1[dt][eid][nid]
01044                 e2 = self.e2[dt][eid][nid]
01045                 e3 = self.e3[dt][eid][nid]
01046                 evm = self.evmShear[dt][eid][nid]
01047                 p = (e1+e2+e3)/-3.
01048 
01049                 if nid=='C': nid='CENTER'
01050                 #print "nid = |%r|" %(nid)
01051                 A = [[exx,exy,exz],
01052                      [exy,eyy,eyz],
01053                      [exz,eyz,ezz]]
01054                 (Lambda,v) = eigh(A) # a hermitian matrix is a symmetric-real matrix
01055 
01056                 ([exx,eyy,ezz,exy,eyz,exz,e1,e2,e3,p,evm],isAllZeros) = self.writeFloats13E([exx,eyy,ezz,exy,eyz,exz,e1,e2,e3,p,evm])
01057                 msgA += '0              %8s  X  %13s  XY  %13s   A  %13s  LX%5.2f%5.2f%5.2f  %13s   %-s\n' %(nid,exx,exy,e1,v[0,1],v[0,2],v[0,0],p,evm.strip())
01058                 msgA += '               %8s  Y  %13s  YZ  %13s   B  %13s  LY%5.2f%5.2f%5.2f\n'             %('', eyy,eyz,e2,v[1,1],v[1,2],v[1,0])
01059                 msgA += '               %8s  Z  %13s  ZX  %13s   C  %13s  LZ%5.2f%5.2f%5.2f\n'             %('', ezz,exz,e3,v[2,1],v[2,2],v[2,0])
01060             ###
01061         ###
01062         if f is not None:
01063             f.write(''.join(msg))
01064             msg = ['']
01065         ###
01066         return ''.join(msg)
01067 
01068     def __repr__(self):
01069         if self.nonlinearFactor is not None:
01070             return self.__reprTransient__()
01071 
01072         msg = '---SOLID STRESS---\n'
01073         headers = self.getHeaders()
01074         msg += '%-6s %6s %8s ' %('EID','eType','nodeID')
01075         for header in headers:
01076             msg += '%9s ' %(header)
01077         msg += '\n'
01078         for eid,oxxNodes in sorted(self.oxx.iteritems()):
01079             eType = self.eType[eid]
01080             for nid in sorted(oxxNodes):
01081                 oxx = self.oxx[eid][nid]
01082                 oyy = self.oyy[eid][nid]
01083                 ozz = self.ozz[eid][nid]
01084                 txy = self.txy[eid][nid]
01085                 tyz = self.tyz[eid][nid]
01086                 txz = self.txz[eid][nid]
01087 
01088                 #o1 = self.o1[eid][nid]
01089                 #o2 = self.o2[eid][nid]
01090                 #o3 = self.o3[eid][nid]
01091                 ovm = self.ovmShear[eid][nid]
01092                 msg += '%-6i %6s %8s ' %(eid,eType,nid)
01093                 vals = [oxx,oyy,ozz,txy,tyz,txz,ovm]
01094                 for val in vals:
01095                     if abs(val)<1e-6:
01096                         msg += '%9s ' %('0')
01097                     else:
01098                         msg += '%9i ' %(val)
01099                     ###
01100                 msg += '\n'
01101                 #msg += "eid=%-4s eType=%-6s nid=%-2i oxx=%-5i oyy=%-5i ozz=%-5i txy=%-5i tyz=%-5i txz=%-5i ovm=%-5i\n" %(eid,eType,nid,oxx,oyy,ozz,txy,tyz,txz,ovm)
01102             ###
01103         ###
01104         return msg
01105 
01106     def __repr__(self):
01107         if self.nonlinearFactor is not None:
01108             return self.__reprTransient__()
01109 
01110         msg = '---SOLID STRAIN---\n'
01111         headers = self.getHeaders()
01112         msg += '%-6s %6s %8s ' %('EID','eType','nodeID')
01113         for header in headers:
01114             msg += '%10s ' %(header)
01115         msg += '\n'
01116         for eid,exxNodes in sorted(self.exx.iteritems()):
01117             eType = self.eType[eid]
01118             for nid in sorted(exxNodes):
01119                 exx = self.exx[eid][nid]
01120                 eyy = self.eyy[eid][nid]
01121                 ezz = self.ezz[eid][nid]
01122                 exy = self.exy[eid][nid]
01123                 eyz = self.eyz[eid][nid]
01124                 exz = self.exz[eid][nid]
01125                 evm = self.evmShear[eid][nid]
01126                 msg += '%-6i %6s %8s ' %(eid,eType,nid)
01127                 vals = [exx,eyy,ezz,exy,eyz,exz,evm]
01128                 for val in vals:
01129                     if abs(val)<1e-6:
01130                         msg += '%10s ' %('0')
01131                     else:
01132                         msg += '%10.3e ' %(val)
01133                     ###
01134                 msg += '\n'
01135                 #msg += "eid=%-4s eType=%-6s nid=%-2i exx=%-5i eyy=%-5i ezz=%-5i exy=%-5i eyz=%-5i exz=%-5i evm=%-5i\n" %(eid,eType,nid,exx,eyy,ezz,exy,eyz,exz,evm)
01136             ###
01137         ###
01138         return msg
01139 
01140     def __reprTransient__(self):
01141         #return ''
01142         msg = ['---SOLID STRAIN---\n']
01143         headers = self.getHeaders()
01144         msg2 = '%-6s %6s %8s ' %('EID','eType','nodeID')
01145         for header in headers:
01146             msg2 += '%10s ' %(header)
01147         msg2 += '\n'
01148         msg.append(msg2)
01149         for dt,exxs in sorted(self.exx.iteritems()):
01150             msg.append('%s = %g\n' %(self.dataCode['name'],dt))
01151             for eid,exxNodes in sorted(exxs.iteritems()):
01152                 eType = self.eType[eid]
01153                 msg2 = ''
01154                 for nid in sorted(exxNodes):
01155                     exx = self.exx[dt][eid][nid]
01156                     eyy = self.eyy[dt][eid][nid]
01157                     ezz = self.ezz[dt][eid][nid]
01158                     exy = self.exy[dt][eid][nid]
01159                     eyz = self.eyz[dt][eid][nid]
01160                     exz = self.exz[dt][eid][nid]
01161                     evm = self.evmShear[dt][eid][nid]
01162                     msg2 += '%-6i %6s %8s ' %(eid,eType,nid)
01163                     vals = [exx,eyy,ezz,exy,eyz,exz,evm]
01164                     for val in vals:
01165                         if abs(val)<1e-6:
01166                             msg2 += '%10s ' %('0')
01167                         else:
01168                             msg2 += '%10.3e ' %(val)
01169                         ###
01170                     msg2 += '\n'
01171                     msg.append(msg2)
01172                     #msg += "eid=%-4s eType=%-6s nid=%-2i exx=%-5i eyy=%-5i ezz=%-5i exy=%-5i eyz=%-5i exz=%-5i evm=%-5i\n" %(eid,eType,nid,exx,eyy,ezz,exy,eyz,exz,evm)
01173                 ###
01174             ###
01175         ###
01176         return ''.join(msg)
 All Classes Namespaces Files Functions Variables