pyNastran
0.5.0
pyNastran BDF Reader/Writer, OP2 Parser, and GUI
|
00001 ## GNU Lesser General Public License 00002 ## 00003 ## Program pyNastran - a python interface to NASTRAN files 00004 ## Copyright (C) 2011-2012 Steven Doyle, Al Danial 00005 ## 00006 ## Authors and copyright holders of pyNastran 00007 ## Steven Doyle <mesheb82@gmail.com> 00008 ## Al Danial <al.danial@gmail.com> 00009 ## 00010 ## This file is part of pyNastran. 00011 ## 00012 ## pyNastran is free software: you can redistribute it and/or modify 00013 ## it under the terms of the GNU Lesser General Public License as published by 00014 ## the Free Software Foundation, either version 3 of the License, or 00015 ## (at your option) any later version. 00016 ## 00017 ## pyNastran is distributed in the hope that it will be useful, 00018 ## but WITHOUT ANY WARRANTY; without even the implied warranty of 00019 ## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 00020 ## GNU General Public License for more details. 00021 ## 00022 ## You should have received a copy of the GNU Lesser General Public License 00023 ## along with pyNastran. If not, see <http://www.gnu.org/licenses/>. 00024 ## 00025 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