pyNastran  0.5.0
pyNastran BDF Reader/Writer, OP2 Parser, and GUI
generalMath.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 from __future__ import print_function
00026 import sys
00027 
00028 from numpy import (float32, float64, complex64, complex128,
00029                    array, cross, allclose, zeros, matrix, insert, diag)
00030 from numpy.linalg import norm #, solve
00031 
00032 from scipy.linalg      import solve_banded
00033 from scipy.interpolate import splrep,splev
00034 from scipy.integrate   import quad
00035 
00036 if sys.version_info < (3,0):
00037     # "fixes" bug where scipy screws up return code handling
00038     import scipy.weave
00039 
00040 def integrateLine(x,y):
00041     """
00042     Integrates a line of length 1.0
00043     @param x the independent variable
00044     @param y the   dependent variable
00045     """
00046     if len(set(y))==1:
00047         return y[0] # (x1-x0 = 1., so yBar*1 = yBar)
00048     try:
00049         assert len(x)==len(y), 'x=%s y=%s' %(x,y)
00050         spline = buildSpline(x,y)
00051         #y = splev(xi,spline)
00052         out = quad(scipy.splev,0.,1.,args=(spline))  # now integrate the area
00053         A = out[0]
00054     except:
00055         print('spline Error x=%s y=%s' %(x,y))
00056         raise
00057     return A
00058 
00059 def evaluatePositiveSpline(x,*args):
00060     spline,minValue = args
00061     y = splev([x],spline)
00062     return max(y,minValue)
00063 
00064 def buildSpline(x,y):
00065     """
00066     builds a cubic spline or 1st order spline if there are less than 3 terms
00067     @param x the independent variable
00068     @param y the   dependent variable
00069     @note a 1st order spline is the same as linear interpolation
00070     """
00071     if len(x)<3:
00072         spline = splrep(x,y,k=1)  # build a linearly interpolated representation
00073     else:
00074         spline = splrep(x,y)  # build a cubic spline representation
00075     return spline
00076 
00077 def integratePositiveLine(x,y,minValue=0.):
00078     """
00079     Integrates a line of length 1.0
00080     @param x the independent variable
00081     @param y the   dependent variable
00082     """
00083     if len(set(y))==1:
00084         return y[0] # (x1-x0 = 1., so yBar*1 = yBar)
00085     try:
00086         assert len(x)==len(y), 'x=%s y=%s' %(x,y)
00087         spline = buildSpline(x,y)
00088         out = quad(evaluatePositiveSpline,0.,1.,args=(spline,minValue))  # now integrate the area
00089         A = out[0]
00090     except:
00091         raise RuntimeError('spline Error x=%s y=%s' %(x,y))
00092     return A
00093 
00094 def reduceMatrix(matA,nids):
00095     """
00096     takes a list of ids and removes those rows and cols
00097     """
00098     nRows = len(nids)
00099     matB = matrix(zeros((nRows,nRows),'d') )
00100 
00101     for i,irow in enumerate(nids):
00102         for j,jcol in enumerate(nids):  
00103             matB[i,j] = matA[irow,jcol]
00104     return matB
00105 
00106 def isListRanged(a,List,b):
00107     """
00108     Returns true if a<= x <= b
00109     or a-x < 0 < b-x
00110     """
00111     for x in List:
00112         if not isFloatRanged(a,x,b):
00113             return False
00114         ###
00115     ###
00116     return True
00117 
00118 def isFloatRanged(a,x,b):
00119     """
00120     Returns true if a<= x <= b
00121     or a-x < 0 < b-x
00122     """
00123     if not a<x:
00124         if not allclose(x,a):
00125            return False 
00126         ###
00127     ###
00128     if not x<b:
00129         if not allclose(x,b):
00130             return False
00131         ###
00132     ###
00133     return True
00134 
00135 def printAnnotatedMatrix(A,rowNames=None,tol=1e-8):
00136     """
00137     takes a list/dictionary and annotates the row number with that value
00138     indicies go from 0 to N
00139     """
00140     if rowNames==None: rowNames=[i for i in xrange(A.shape[0])]
00141     B = array(A)
00142     msg = ''
00143     (nr,nc) = B.shape
00144     for i in xrange(nr):
00145         rowName = str(rowNames[i])
00146         msg += '%-2s' %(rowName) +' '+ListPrint(B[i,:],tol)+'\n'
00147     return msg
00148 
00149 def printMatrix(A,tol=1e-8):
00150     B = array(A)
00151     msg = ''
00152     (nr,nc) = B.shape
00153     for i in xrange(nr):
00154         msg += ListPrint(B[i,:],tol)+'\n'
00155     #for row in A:
00156     #    msg += ListPrint(row)+'\n'
00157     ###
00158     return msg
00159 ###
00160 
00161 def ListPrint(listA,tol=1e-8):
00162     if len(listA)==0:
00163         return '[]'
00164     ###
00165 
00166     msg = '['
00167     for a in listA:
00168         if isinstance(a,str):
00169             msg += ' %s,' %(a)
00170         elif isinstance(a,float) or isinstance(a,float32) or isinstance(a,float64):
00171             if abs(a)<tol:
00172                 a = 0.
00173             msg += ' %-3.2g,' %(a)
00174         elif isinstance(a,int):
00175             if abs(a)<tol:
00176                 a = 0.
00177             msg += ' %3i,' %(a)
00178         elif isinstance(a,complex64) or isinstance(a,complex128):
00179             r= '%.4g'%(a.real)
00180             i= '%+.4gj'%(a.imag)
00181 
00182             if abs(a.real)<tol:
00183                 r = '0'
00184             if abs(a.imag)<tol:
00185                 i = ''
00186             temp = '%4s%4s' %(r,i)
00187             msg += ' %8s,' %(temp.strip())
00188         else:
00189             try:
00190                 print("type(a) is not supported...%s" %(type(a)))
00191                 msg += ' %g,' %(a)
00192             except TypeError:
00193                 print("a = |%s|" %(a))
00194                 raise
00195             ###
00196         ###
00197     ###
00198     msg = msg[:-1]
00199     msg += ' ]'
00200     return msg
00201 
00202 def augmentedIdentity(A):
00203     """
00204     Creates an Identity Matrix augmented with zeros.
00205     The location of the extra zeros depends on A.
00206     """
00207     (nx,ny) = A.shape
00208     #(ny,nx) = A.shape
00209     I = zeros([nx,ny],'d')
00210     
00211     for i in xrange(nx):
00212         if i==nx or i==ny:
00213             break
00214         I[i][i] = 1.
00215     return I
00216 
00217 def solveTridag(A, D):
00218      # Find the diagonals
00219      ud = insert(diag(A,1), 0, 0) # upper diagonal
00220      d  = diag(A) # main diagonal
00221      ld = insert(diag(A,-1), len(d)-1, 0) # lower diagonal
00222    
00223      # simplified matrix
00224      ab = matrix([ud,d,ld])
00225      #print "ab = ",ab
00226      return solve_banded((1, 1), ab, D,overwrite_ab=True,overwrite_b=True)
00227 
00228 def Area(a,b):
00229     return 0.5*norm(cross(a,b))
00230 
00231 def AreaNormal(nodes):
00232     """
00233     Returns area,unitNormal
00234     n = Normal = a x b
00235     Area   = 1/2 * |a x b|
00236     V = <v1,v2,v3>
00237     |V| = sqrt(v1^0.5+v2^0.5+v3^0.5) = norm(V)
00238     
00239     Area = 0.5 * |n|
00240     unitNormal = n/|n|
00241     """
00242     (n0,n1,n2) = nodes
00243     a = n0-n1
00244     b = n0-n2
00245     vector = cross(a,b)
00246     length = norm(vector)
00247     normal = vector/length
00248     area = 0.5*length
00249     if allclose(norm(normal),1.)==False:
00250         print("a = ",a)
00251         print("b = ",b)
00252         print("normal = ",normal)
00253         print("length = ",length)
00254         raise RuntimeError('check...')
00255     return (area,normal)
00256 
00257 def Triangle_AreaCentroidNormal(nodes):
00258     """Returns area,centroid,unitNormal"""
00259     (area,normal) = AreaNormal(nodes)
00260     centroid = Centroid(*nodes)
00261     return (area,centroid,normal)
00262 
00263 def Normal(a,b):
00264     """finds the unit normal vector of 2 vectors"""
00265     vector = cross(a,b)
00266     length = norm(vector)
00267     normal = vector/length
00268     assert allclose(norm(normal),1.)
00269     return normal
00270 
00271 def Centroid(A,B,C):
00272     """returns the centroid of a triangle"""
00273     #print "type(A,B,C) = ",type(A),type(C),type(B)
00274     centroid = (A+B+C)/3.
00275     return centroid
 All Classes Namespaces Files Functions Variables