pyNastran  0.5.0
pyNastran BDF Reader/Writer, OP2 Parser, and GUI
dmig.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 # pylint: disable=C0103,R0902,R0904,R0914
00026 
00027 from __future__ import division, print_function
00028 from itertools import izip
00029 from math import log
00030 #from math import (sin,sinh,cos,cosh,tan,tanh,sqrt,atan,atan2,acosh,acos,asin,
00031 #                  asinh,atanh) #,atanh2   # going to be used by DEQATN
00032 
00033 from numpy import zeros #average
00034 from scipy.sparse import coo_matrix
00035 
00036 from pyNastran.bdf.cards.baseCard import BaseCard
00037 
00038 def ssq(*listA):
00039     """
00040     sum of squares
00041     @note used for DEQATN
00042     """
00043     out = 0.
00044     for x in listA:
00045         out += x*x
00046     return out
00047 
00048 def sum2(*listA):
00049     """
00050     sum of listA
00051     @note used for DEQATN
00052     """
00053     return sum(listA)
00054 
00055 def mod(x,y):
00056     """
00057     x%y
00058     @note used for DEQATN
00059     """
00060     return x%y
00061 
00062 def logx(x,y):
00063     """
00064     log base x of y
00065     @note used for DEQATN
00066     """
00067     log(y,x)
00068 
00069 def dim(x,y):
00070     """
00071     @note used for DEQATN
00072     """
00073     return x-min(x,y)
00074 
00075 def db(p,pref):
00076     """
00077     sound pressure in decibels
00078     would capitalize it, but you wouldnt be able to call the function...
00079     """
00080     return 20.*log(p/pref)
00081 
00082 class DEQATN(BaseCard):# needs work...
00083     type = 'DEQATN'
00084     def __init__(self,card=None,data=None):
00085         newCard = ''
00086         foundNone = False
00087         for field in card.card:
00088             if foundNone is False and field is not None:
00089                 newCard += field+','
00090                 foundNone = True
00091             elif foundNone is True and field is not None:
00092                 newCard += field
00093 
00094         #if len(card.card)>1:
00095             #print "card.card = ",card.card
00096             #line0 = ','.join(card.card)
00097         #else:
00098             #line0 = ''.join(card.card)
00099         line0 = newCard
00100         self.eqID = line0[8:16]
00101         
00102         assert len(self.eqID)==8,'len(eqID)==%s' %(len(self.eqID))
00103         eq = line0[16:]
00104         eq = eq.replace(' ','').lower()
00105         (self.name,self.eq) = eq.split('=')
00106         #print("EQ = %s" %(self.eq))
00107         
00108     def evaluate(self,args):
00109         #eqLow = self.eq.lower()
00110         #eval(self.eq)
00111         pass
00112 
00113     def __repr__(self):
00114         eq = self.name+'='+self.eq
00115         eqLine = eq[0:56]
00116         eq = eq[56:]
00117         fields = ['DEQATN  ','%8s'%(self.eqID), eqLine]
00118 
00119         if len(eq):
00120             eqLine = eq[0:72]
00121             eq = eq[72:]
00122             fields += ['        '+eqLine]
00123         return ''.join(fields)
00124 
00125 class NastranMatrix(BaseCard):
00126     """
00127     Base class for the DMIG, DMIJ, DMIJI, DMIK matrices
00128     """
00129     def __init__(self, card=None, data=None):
00130         self.name = card.field(1)
00131         #zero
00132 
00133         ## 4-Lower Triangular; 5=Upper Triangular; 6=Symmetric; 8=Identity (m=nRows, n=m)
00134         self.ifo   = int(card.field(3))
00135         ## 1-Real, Single Precision; 2=Real,Double Precision; 3=Complex, Single; 4=Complex, Double
00136         self.tin   = int(card.field(4))
00137         ## 0-Set by cell precision
00138         self.tout  = int(card.field(5,0))
00139         
00140         self.polar = card.field(6)
00141         self.ncol  = card.field(8)
00142 
00143         self.GCj = []
00144         self.GCi = []
00145         self.Real = []
00146         if self.isComplex():
00147             self.Complex = []
00148     
00149     def writeCodeAster(self):
00150         """
00151         assume set 1 = MAAX1,MAAX2, etc. and 100/n % on each
00152         """
00153         
00154         # for real combination
00155         comm =  'K_Mtx_AB=COMB_MATR_ASSE(COMB_R=(\n'
00156         comm += '    _F(MATR_ASSE = K_Mtx_A,COEF_R = 1.),\n'
00157         comm += '    _F(MATR_ASSE = K_Mtx_B,COEF_R = 1.)));\n'
00158 
00159         # for complex combination
00160 
00161         comm += "K_Mtx_AB=COMB_MATR_ASSE(COMB_C=(\n"
00162         comm += "_F(MATR_ASSE=K_Mtx_A,COEF_C=('RI',0.7,0.3,),)\n"
00163         comm += "_F(MATR_ASSE=K_Mtx_B,COEF_C=('RI',0.7,0.3,),),),);\n"
00164         comm = 'K_Mtx=ASSE_MATRICE(MATR_ELEM=ElMtx_K,NUME_DDL=%s,);'
00165         return comm
00166 
00167     def addColumn(self,card=None,data=None):
00168         #print "hi column"
00169         Gj = card.field(2)
00170         Cj = card.field(3)
00171         
00172         nFields = card.nFields()
00173         nLoops = (nFields-5)//4
00174         minLoops = nLoops-1
00175         if minLoops<=0:
00176             minLoops = 1
00177         #assert nFields <= 8,'nFields=%s' %(nFields)
00178 
00179         #print "minLoops = ",minLoops
00180         #print "nLoops   = ",nLoops
00181         for i in xrange(minLoops):
00182             self.GCj.append((Gj,Cj))
00183 
00184         if self.isComplex():
00185             for i in xrange(minLoops):
00186                 n = 5+4*i
00187                 Gi = card.field(n)
00188                 Ci = card.field(n+1)
00189                 self.GCi.append((Gi,Ci))
00190                 self.Real.append(card.field(n+2))
00191                 self.Complex.append(card.field(n+3))
00192         else:
00193             for i in xrange(minLoops):
00194                 n = 5+4*i
00195                 Gi = card.field(n)
00196                 Ci = card.field(n+1)
00197                 self.GCi.append((Gi,Ci))
00198                 self.Real.append(card.field(n+2))
00199         
00200         assert len(self.GCj)==len(self.GCi),'(len(GCj)=%s len(GCi)=%s' %(len(self.GCj),len(self.GCi))
00201         #if self.isComplex():
00202             #self.Complex(card.field(v)
00203 
00204     def getMatrix(self,isSparse=False):
00205         """
00206         builds the Matrix
00207         @param self the object pointer
00208         @param isSparse should the matrix be returned as a sparse matrix (default=True).  Slower for dense matrices.
00209         @retval M the matrix
00210         @retval rows dictionary of keys=rowID,    values=(Grid,Component) for the matrix
00211         @retval cols dictionary of keys=columnID, values=(Grid,Component) for the matrix
00212         @warning isSparse WILL fail
00213         """
00214         i=0
00215         rows = {}
00216         rowsReversed = {}
00217         for GCi in self.GCi:
00218             if GCi not in rows:
00219                 rows[GCi] = i
00220                 rowsReversed[i] = GCi
00221                 i+=1
00222         #nRows = len(rows2)
00223 
00224         j=0
00225         cols = {}
00226         colsReversed = {}
00227         for GCj in self.GCj:
00228             if GCj not in cols:
00229                 cols[GCj] = j
00230                 colsReversed[j] = GCj
00231                 j+=1
00232         #nCols = len(cols2)
00233         
00234         #A = ss.lil_matrix((3,3), dtype='d') # double precision
00235 
00236         #rows=[]; cols=[]; data=[]
00237         #for i in xrange(3):
00238         #    for j in xrange(3):
00239         #        k = float((i+1)*(j+1))
00240         #        rows.append(i)
00241         #        cols.append(j)
00242         #        data.append(k)
00243         #        A[i,j] = k
00244         #    ###
00245         ###
00246 
00247         #isSparse = False
00248         if isSparse:
00249             data = []
00250             rows2 = []
00251             cols2 = []
00252 
00253             nrows = 0
00254             ncols = 0
00255             if self.isComplex():
00256                 for (GCj, GCi, reali, complexi) in izip(self.GCj, self.GCi, self.Real, self.Complex):
00257                     i = rows[GCi]
00258                     j = cols[GCj]
00259                     nrows = max(i, nrows)
00260                     ncols = max(j, ncols)
00261                     rows2.append(i)
00262                     cols2.append(j)
00263                     data.append(complex(reali, complexi))
00264             else:
00265                 for (GCj, GCi) in izip(self.GCj, self.GCi):
00266                     i = rows[GCi]
00267                     j = cols[GCj]
00268                     nrows = max(i,nrows)
00269                     ncols = max(j,ncols)
00270                     rows2.append(i)
00271                     cols2.append(j)
00272                 data = self.Real
00273             #print("i=%s j=%s len(rows2)=%s len(cols2)=%s len(data)=%s"
00274             #    % (i, j, len(self.GCi), len(self.GCj), len(data)))
00275             # ,dtype=Format
00276             #print(rows2)
00277 
00278             #print("nrows=%s ncols=%s" %(nrows, ncols))
00279             if self.ifo in [1,6]:
00280                 nrows = max(nrows,ncols)
00281                 ncols = nrows
00282             #print("nrows=%s ncols=%s" %(nrows, ncols))
00283             
00284             dType = self.getDType(self.tin)
00285             #A = coo_matrix( (entries,(rows,cols)),shape=(nrows,ncols),dtype=dType) # test
00286             M = coo_matrix( (data, (self.GCi, self.GCj)), shape=(nrows,ncols), dtype=dType)
00287             #M = coo_matrix( (data,(self.GCi,self.GCj)),shape=(i,j)) # old
00288             #M = coo_matrix( (data,(self.GCi,self.GCj)),shape=(nrows,ncols))
00289             #print M.todense()
00290             #print M
00291         else:
00292             if self.isComplex():
00293                 M = zeros((i,j), dtype='complex128')
00294                 for (GCj, GCi, reali, complexi) in izip(self.GCj, self.GCi, self.Real, self.Complex):
00295                     i = rows[GCi]
00296                     j = cols[GCj]
00297                     M[i,j] = complex(reali,complexi)
00298             else:
00299                 M = zeros((i,j), dtype='float64')
00300                 for (GCj, GCi, reali) in izip(self.GCj, self.GCi, self.Real):
00301                     i = rows[GCi]
00302                     j = cols[GCj]
00303                     M[i,j] = reali
00304                 ###
00305             ###
00306         #print(M)
00307         return (M,rowsReversed,colsReversed)
00308 
00309     def rename(self,newName):
00310         self.name = newName
00311 
00312     def isReal(self):
00313         return not self.isComplex()
00314 
00315     def isComplex(self):
00316         if self.tin in [3,4]:
00317             return True
00318         return False
00319 
00320     def getDType(self,Type):
00321         if Type==1:
00322             dType = 'float32'
00323         elif Type==2:
00324             dType = 'float64'
00325         elif Type==3:
00326             dType = 'complex64'
00327         elif Type==4:
00328             dType = 'complex128'
00329         elif Type==0:
00330             if self.isComplex():
00331                 dType = 'complex128'
00332             else:
00333                 dType = 'float64'
00334         else:
00335             raise RuntimeError("invalid option for matrix format")
00336         return dType
00337 
00338     def __repr__(self):
00339         """
00340         @todo support double precision
00341         """
00342         msg = '\n$'+'-'*80
00343         msg += '\n$ %s Matrix %s\n' %(self.type, self.name)
00344         fields = [self.type, self.name, 0, self.ifo, self.tin, self.tout, self.polar, None, self.ncol]
00345         msg += self.printCard(fields)
00346 
00347         if self.isComplex():
00348             for (GCi, GCj, reali, imagi) in izip(self.GCi, self.GCj, self.Real, self.Complex):
00349                 fields = [self.type, self.name, GCj[0], GCj[1], None, GCi[0], GCi[1], reali, imagi]
00350                 msg += self.printCard(fields)
00351         else:
00352             for (GCi, GCj, reali) in izip(self.GCi, self.GCj, self.Real):
00353                 fields = [self.type, self.name, GCj[0], GCj[1], None, GCi[0], GCi[1], reali, None]
00354                 msg += self.printCard(fields)
00355         return msg
00356 
00357 
00358 class DMIG(NastranMatrix):
00359     """
00360     Defines direct input matrices related to grid, extra, and/or scalar points.
00361     The matrix is defined by a single header entry and one or more column
00362     entries. A column entry is required for each column with nonzero elements.
00363     """
00364     type = 'DMIG'
00365     def __init__(self, card=None, data=None):
00366         NastranMatrix.__init__(self, card, data)
00367 
00368 class DMIJ(NastranMatrix):
00369     """
00370     Direct Matrix Input at js-Set of the Aerodynamic Mesh
00371     Defines direct input matrices related to collation degrees-of-freedom
00372     (js-set) of aerodynamic mesh points for CAERO1, CAERO3, CAERO4 and CAERO5
00373     and for the slender body elements of CAERO2. These include W2GJ, FA2J and
00374     input pressures and downwashes associated with AEPRESS and AEDW entries.
00375     The matrix is described by a single header entry and one or more column
00376     entries. A column entry is required for each column with nonzero elements.
00377     For entering data for the interference elements of a CAERO2, use DMIJI
00378     or DMI.
00379     """
00380     type = 'DMIJ'
00381     def __init__(self,card=None,data=None):
00382         NastranMatrix.__init__(self,card,data)
00383 
00384 class DMIJI(NastranMatrix):
00385     """
00386     Direct Matrix Input at js-Set of the Interference Body
00387     Defines direct input matrices related to collation degrees-of-freedom
00388     (js-set) of aerodynamic mesh points for the interference elements of CAERO2.
00389     These include W2GJ, FA2J and input pressures and downwashes associated with
00390     AEPRESS and AEDW entries. The matrix is described by a single header entry
00391     and one or more column entries. A column entry is required for each column
00392     with nonzero elements.  For entering data for the slender elements of a
00393     CAERO2, or a CAERO1, 3, 4 or 5 use DMIJ or DMI.
00394     """
00395     type = 'DMIJI'
00396     def __init__(self, card=None, data=None):
00397         NastranMatrix.__init__(self, card, data)
00398 
00399 class DMIK(NastranMatrix):
00400     """
00401     Direct Matrix Input at ks-Set of the Aerodynamic Mesh
00402     Defines direct input matrices related to physical (displacement)
00403     degrees-of-freedom (ks-set) of aerodynamic grid points. These include WKK,
00404     WTFACT and input forces associated with AEFORCE entries. The matrix is
00405     described by a single header entry and one or more column entries. A column
00406     entry is required for each column with nonzero elements.
00407     """
00408     type = 'DMIK'
00409     def __init__(self,card=None,data=None):
00410         NastranMatrix.__init__(self,card,data)
00411 
00412 class DMI(BaseCard):
00413     type = 'DMI'
00414     def __init__(self,card=None,data=None):
00415         self.name = card.field(1)
00416         #zero
00417         
00418         ## Form of the matrix:  1=Square (not symmetric); 2=Rectangular;
00419         ## 3=Diagonal (m=nRows,n=1);  4=Lower Triangular; 5=Upper Triangular;
00420         ## 6=Symmetric; 8=Identity (m=nRows, n=m)
00421         self.form  = int(card.field(3))
00422         
00423         ## 1-Real, Single Precision; 2=Real,Double Precision;
00424         ## 3=Complex, Single; 4=Complex, Double
00425         self.tin   = int(card.field(4))
00426         
00427         ## 0-Set by cell precision
00428         self.tout  = int(card.field(5, 0))
00429         
00430         self.nRows = int(card.field(7))
00431         self.nCols = int(card.field(8))
00432 
00433         self.GCj = []
00434         self.GCi = []
00435         self.Real = []
00436         
00437         if self.isComplex():
00438             self.Complex = []
00439 
00440     def addColumn(self, card=None, data=None):
00441 
00442         if not self.isComplex(): # real
00443             self.readReal(card)
00444 
00445     def readReal(self, card):
00446         ## column number
00447         j = card.field(2)
00448 
00449         # counter
00450         i = 0
00451         fields = card.fields(3)
00452 
00453         # Real, starts at A(i1,j), goes to A(i2,j) in a column
00454         while i<len(fields):
00455             i1 = fields[i]
00456             #print "i1 = ",i1
00457             if isinstance(i1,int):
00458                 i+=1
00459                 isDoneReadingFloats = False
00460                 while not isDoneReadingFloats and i<len(fields):
00461                     #print "i=%s len(fields)=%s" %(i,len(fields))
00462                     realValue = fields[i]
00463                     if isinstance(realValue, int):
00464                         isDoneReadingFloats = True
00465                     elif isinstance(realValue, float):
00466                         self.GCj.append(j)
00467                         self.GCi.append(i1)
00468                         self.Real.append(realValue)
00469                         #print "i=%s j=%s value=%s" %(i1,j,realValue)
00470                         i+=1
00471                     else:
00472                         #print "*i=%s j=%s value=%s type=%s" %(i1,j,realValue,type(realValue))
00473                         realValue = self.Real[-1]
00474                         #print "*i=%s j=%s value=%s" %(i1,j,realValue)
00475                         endI = fields[i+1]
00476                         #print "*i=%s endI=%s j=%s value=%s" %(i1,endI,j,realValue)
00477                         for ii in xrange(i1, endI+1):
00478                             self.GCj.append(j)
00479                             self.GCi.append(ii)
00480                             self.Real.append(realValue)
00481                         ###
00482                         #print "i = ",3+i
00483                         #print 'field i=',fields[i]
00484                         i+=1
00485                         isDoneReadingFloats = True
00486                     ###
00487             ###
00488 
00489     def readComplex(self,card):
00490         msg = 'complex matrices not supported in the DMI reader...'
00491         raise NotImplementedError(msg)
00492         ## column number
00493         j = card.field(2)
00494 
00495         # counter
00496         i = 0
00497         fields = card.fields(3)
00498 
00499         # Complex, starts at A(i1,j)+imag*A(i1,j), goes to A(i2,j) in a column
00500         while i<len(fields):
00501             i1 = fields[i]
00502             i+=1
00503             isDoneReadingFloats = False
00504             asdf
00505             while not isDoneReadingFloats and i<len(fields):
00506                 #print("i=%s len(fields)=%s" %(i, len(fields)))
00507                 realValue = fields[i]
00508                 if isinstance(floatValue, int):
00509                     isDoneReadingFloats = True
00510                 elif isinstance(realValue, float):
00511                     complexValue = fields[i+1]
00512                     self.GCj.append(j)
00513                     self.GCi.append(i1)
00514                     self.Real.append(realValue)
00515                     self.Complex.append(complexValue)
00516                     i+=2
00517                 else:
00518                     asdf
00519                 ###
00520             ### while floats
00521         ### while fields
00522 
00523     def rename(self,newName):
00524         self.name = newName
00525 
00526     def isReal(self):
00527         return not self.isComplex()
00528 
00529     def isComplex(self):
00530         if self.tin in [3,4]:
00531             return True
00532         return False
00533 
00534     def __repr__(self):
00535         """
00536         @todo support shortened output format.  There's a stupidly low 1000
00537         DMI cap, I assume this is entries and not matrices.
00538         @todo support double precision
00539         """
00540         msg = '\n$'+'-'*80
00541         msg += '\n$ %s Matrix %s\n' %(self.type,self.name)
00542         fields = [self.type, self.name, 0, self.form, self.tin, self.tout, None, self.nRows, self.nCols]
00543         msg += self.printCard(fields)
00544         #msg += self.printCard(fields,size=16,isD=False)
00545 
00546         if self.isComplex():
00547             for (GCi, GCj, reali, imagi) in izip(self.GCi, self.GCj, self.Real, self.Complex):
00548                 fields = [self.type, self.name, GCj, GCi, reali, imagi]
00549                 msg += self.printCard(fields)
00550         else:
00551             for (GCi, GCj, reali) in izip(self.GCi, self.GCj, self.Real):
00552                 fields = [self.type, self.name, GCj, GCi, reali]
00553                 msg += self.printCard(fields)
00554         return msg
00555     
 All Classes Namespaces Files Functions Variables