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 sympy import Symbol,simplify,Rational,Matrix 00026 00027 ## The purpose of this code is to automatically generate 00028 ## the moment of inertia equations. only the locations of 00029 ## The nodes of the elements are required, as well as the 00030 ## triangular sections that are generated. 00031 ## For example, [1,a,b] means a positive triangular section (+1) 00032 ## with the 3 nodes at a, b, and <0,0>. 00033 ## This is easily understood if the triangle creates an area. 00034 ## Some are positive, some are negative, but if you sum them all 00035 ## up you get the total area of the section. This allows you 00036 ## to generate the moments of inertia for complex sections. 00037 00038 def bar(): 00039 """ 00040 y 00041 ^ 00042 | 00043 1----|----2 00044 | | | 00045 | | | h 00046 | o----------> x 00047 | | 00048 | | 00049 4---------3 00050 b 00051 00052 b = w1 00053 h = h1 00054 """ 00055 print "---Bar---" 00056 #A = b*h 00057 h = Symbol('h') 00058 b = Symbol('b') 00059 half = Rational(1,2) 00060 A = b*h*half 00061 00062 p1 = [-b*half, h*half] 00063 #p2 = [ b*half, h*half] 00064 #p3 = [ b*half,-h*half] 00065 #p4 = [-b*half,-h*half] 00066 00067 p2 = addx(p1,b) 00068 p3 = suby(p2,h) 00069 p4 = subx(p3,b) 00070 print "p1 = ",p1 00071 print "p2 = ",p2 00072 print "p3 = ",p3 00073 print "p4 = ",p4 00074 Atri = half*p1[0]*p1[0] 00075 print "Atri = ",Atri 00076 00077 sections = [ [1,p1,p2],[1,p2,p3],[1,p3,p4],[1,p4,p1] ] 00078 Ixx(sections) 00079 00080 def cross(): 00081 """ 00082 y 00083 b1 ^ b3 00084 | 00085 b2 00086 a---p1--b 00087 | | | h1 00088 k--l | c-----d 00089 | | | 00090 | o---------p2-> x h2 00091 | | 00092 j--i f-----e 00093 | | h3 00094 | | 00095 h-------g 00096 """ 00097 #given 00098 #k-l = 0.5 D1 00099 #k-j = D4 00100 #h-g = D2 00101 #b-g = D3 00102 print "---Cross---" 00103 00104 D1 = Symbol('D1') 00105 D2 = Symbol('D2') 00106 D3 = Symbol('D3') 00107 D4 = Symbol('D4') 00108 o2 = Rational(1,2) 00109 # bases 00110 b2 = D2 00111 b1 = b3 = D1*o2 00112 00113 # heights 00114 h2 = D4 00115 h1 = h3 = (D3-D4)*o2 00116 00117 p1 = [0, D3*o2] 00118 p2 = [(D2+D1)*o2, 0] 00119 00120 a = subx(p1,D2*o2) 00121 b = addx(p1,D2*o2) 00122 c = suby(b,h1) 00123 d = addx(c,b3) 00124 e = suby(d,h2) 00125 f = subx(e,b3) 00126 g = suby(f,h3) 00127 h = subx(g,b2) 00128 i = addy(h,h3) 00129 j = subx(i,b1) 00130 k = addy(j,h2) 00131 l = addx(k,b1) 00132 A = addy(l,h1) 00133 #print "a = ",a 00134 #print "b = ",b 00135 00136 #print "d = ",d 00137 #print "e = ",e 00138 00139 #print "g = ",g 00140 #print "h = ",h 00141 00142 #print "k = ",k 00143 #print "j = ",j 00144 00145 #print "A = ",A 00146 sections = [ [1,a,b],[1,b,c],[1,c,d],[1,d,e], 00147 [1,e,f],[1,f,g],[1,g,h],[1,h,i], 00148 [1,i,j],[1,j,k],[1,k,l],[1,l,a] ] 00149 Ixx(sections) 00150 00151 00152 def box(): 00153 """ 00154 y 00155 ^ 00156 | 00157 w1 | w2 w3 00158 | 00159 a-----p1------b 00160 | | | h1 00161 | e--|----f | 00162 | | | | | 00163 | | | | | h2 00164 | | o-------|--> x 00165 | | | | 00166 | h-------g | 00167 | | h3 00168 d-------------c 00169 00170 # given 00171 #a-d = D2 00172 #(f-g)_y = D3 00173 #d-c = D1 00174 #(g-c)_x = D4 00175 """ 00176 print "---Box---" 00177 00178 D1 = Symbol('D1') 00179 D2 = Symbol('D2') 00180 D3 = Symbol('D3') 00181 D4 = Symbol('D4') 00182 00183 # height 00184 h1 = h3 = D3 00185 h2 = D2-h1-h3 00186 00187 # bases 00188 w1 = w3 = D4 00189 w2 = D1-w1-w3 00190 00191 o2 = Rational(1,2) 00192 00193 # outer 00194 p1 = [0,D2*o2] 00195 a = subx(p1,D1*o2) 00196 b = addx(p1,D1*o2) 00197 c = suby(b,D2) 00198 d = subx(c,D1) 00199 00200 # inner 00201 e = suby(p1,D3) 00202 f = addx(e,w2) 00203 g = suby(f,h2) 00204 h = subx(g,w2) 00205 00206 sections = [ [ 1,a,b],[ 1,b,c],[ 1,c,d],[ 1,d,a], 00207 [-1,e,f],[-1,f,g],[-1,g,h],[-1,h,e] ] 00208 Ixx(sections) 00209 00210 def Box1(): 00211 """ 00212 y 00213 ^ 00214 | 00215 w1 | w2 w3 00216 | 00217 a-----p1------b 00218 | | | h1 00219 | e--|----f | 00220 | | | | | 00221 | | | | | h2 00222 | | o-------|--> x 00223 | | | | 00224 | h-------g | 00225 | | h3 00226 d-------------c 00227 00228 # given 00229 #a-b = D1 00230 #a-d = D2 00231 00232 #(b-f)_y = D3 00233 #(g-c)_y = D4 00234 #(g-c)_x = D5 00235 #(d-h)_x = D6 00236 00237 @warning inertia is wrong b/c 0 is assumed to be <0,0> when it should 00238 be the output of the first inertia calculation 00239 """ 00240 print "---Box1---" 00241 00242 D1 = Symbol('D1') 00243 D2 = Symbol('D2') 00244 D3 = Symbol('D3') 00245 D4 = Symbol('D4') 00246 D5 = Symbol('D5') 00247 D6 = Symbol('D6') 00248 00249 # height 00250 h1 = D3 00251 h2 = D2-D3-D4 00252 h3 = D4 00253 00254 # bases 00255 w1 = D6 00256 w2 = D1-D5-D6 00257 w3 = D5 00258 00259 o2 = Rational(1,2) 00260 00261 # outer 00262 p1 = [0,D2*o2] 00263 a = subx(p1,D1*o2) 00264 b = addx(p1,D1*o2) 00265 c = suby(b,D2) 00266 d = subx(c,D1) 00267 00268 # inner 00269 e = suby(p1,D3) 00270 f = addx(e,w2) 00271 g = suby(f,h2) 00272 h = subx(g,w2) 00273 00274 sections = [ [ 1,a,b],[ 1,b,c],[ 1,c,d],[ 1,d,a], 00275 [-1,e,f],[-1,f,g],[-1,g,h],[-1,h,e] ] 00276 Ixx(sections) 00277 00278 00279 def T(): 00280 """ 00281 y 00282 ^ 00283 | 00284 a---k-p1-l---b 00285 | | | h1 00286 | | | 00287 h---g | d---c 00288 | | | 00289 | o----------> x 00290 | | 00291 | | h2 00292 j f----e i 00293 w1 w2 w3 00294 00295 # given 00296 a-b = D1 00297 (b-e)_y = D2 00298 a-h = D3 00299 f-e = D4 00300 @todo Area seems off 00301 @todo Inertias should be relative to the CG 00302 @warning wrong... 00303 """ 00304 print "---T---" 00305 00306 D1 = Symbol('D1') 00307 D2 = Symbol('D2') 00308 D3 = Symbol('D3') 00309 D4 = Symbol('D4') 00310 00311 o2 = Rational(1,2) 00312 w2 = D4 00313 h1 = D3 00314 h2 = D2-h1 00315 w1 = w3 = (D1-D4)*o2 00316 00317 p1 = [0,0] 00318 a = subx(p1,D1*o2) 00319 b = addx(p1,D1*o2) 00320 c = suby(b,h1) 00321 d = subx(c,w3) 00322 e = suby(d,h2) 00323 f = subx(e,w2) 00324 g = addy(f,h2) 00325 h = subx(g,w1) 00326 A = addy(h,h1) 00327 00328 i = suby(c,h2) 00329 j = subx(i,D1) 00330 k = subx(p1,D4*o2) 00331 l = addx(p1,D4*o2) 00332 #print "a = ",a 00333 #print "A = ",A 00334 00335 sections = [ [ 1,a,b],[1,b,c],[1,d,c],[1,d,e],[1,f,e], 00336 [ 1,f,g],[1,g,h],[1,h,a], 00337 [-1,l,d],[-1,d,g],[-1,g,k] ] 00338 Ixx(sections) 00339 00340 00341 def Ixx(sections): 00342 """ 00343 from http://en.wikipedia.org/wiki/Second_moment_of_area 00344 J_{xx} = \frac{1}{12} \sum_{i = 1}^{n-1} ( y_i^2 + y_i y_{i+1} + y_{i+1}^2 ) a_i 00345 J_{yy} = \frac{1}{12} \sum_{i = 1}^{n-1} ( x_i^2 + x_i x_{i+1} + x_{i+1}^2 ) a_i 00346 J_{xy} = \frac{1}{24} \sum_{i = 1}^{n-1} ( x_i y_{i+1} + 2 x_i y_i + 2 x_{i+1} y_{i+1} + x_{i+1} y_i ) a_i 00347 """ 00348 h = Symbol('h') 00349 b = Symbol('b') 00350 Ixx = 0 00351 Iyy = 0 00352 Ixy = 0 00353 half = Rational(1,2) 00354 CG = [0,0] 00355 Area = 0 00356 o3 = Rational(1,3) 00357 for i,section in enumerate(sections): 00358 (sign,p1,p2) = section 00359 #print "p%s = %s" %(i,p1) 00360 #print "p%s = %s" %(i+1,p2) 00361 #Atri = half * abs(p1[0]*p2[1]) 00362 AA = Matrix([p1,p2]) 00363 #print "AA = \n",AA 00364 detAA = AA.det() 00365 msgAA = '%s'%(detAA) 00366 #print "msgAA = ",msgAA 00367 if '-' in msgAA: 00368 detAA *= Rational(-1,1) 00369 Atri = detAA 00370 #Atri = half*b*h 00371 #Atri = (p1[0]*p2[1] - p2[0]*p1[1])*half 00372 ixxi = p1[1]*p1[1] + p1[1]*p2[1] + p2[1]*p2[1] 00373 iyyi = p1[0]*p1[0] + p1[0]*p2[0] + p2[0]*p2[0] 00374 ixyi = p1[0]*p1[1] + 2*p1[0]*p1[1] + 2*p2[0]*p2[1] + p2[0]*p1[1] 00375 Ixxi = ixxi*Atri 00376 Iyyi = iyyi*Atri 00377 Ixyi = ixyi*Atri 00378 #print "Atri = ",Atri 00379 #print "Ixxi = ",Ixxi 00380 Ixx += Ixxi 00381 Iyy += Iyyi 00382 Ixy += Ixyi 00383 #print "Ixx = ",simplify(Ixx),'\n' 00384 CG = [CG[0]+p1[0]*o3+p2[0]*o3,CG[1]+p1[1]*o3+p2[1]*o3] 00385 Area += Atri 00386 CG = [CG[0]/Area,CG[1]/Area] 00387 #print "Ixx = ",simplify(Ixx) 00388 o12 = Rational(1,12) 00389 o24 = Rational(1,24) 00390 Ixx = Ixx*o12 00391 Iyy = Iyy*o12 00392 Ixy = Ixy*o24 00393 J = Ixx+Iyy 00394 print "Area = ",simplify(Area*half) 00395 print "CG = ",simplify(CG) 00396 print "Ixx = ",simplify(Ixx) 00397 print "Iyy = ",simplify(Iyy) 00398 print "Ixy = ",simplify(Ixy) 00399 print "J = ",simplify(J),'\n' 00400 00401 def add(p1,p2,s1=1,s2=1): 00402 x = p1[0]*s1+p2[0]*s2 00403 y = p1[1]*s1+p2[1]*s2 00404 return [simplify(x),simplify(y)] 00405 00406 def sub(p1,p2): 00407 return add(p1,p2,s2=-1) 00408 00409 def subx(p1,p2x): 00410 return sub(p1,[p2x,0]) 00411 00412 def suby(p1,p2y): 00413 return sub(p1,[0,p2y]) 00414 00415 def addx(p1,p2x): 00416 return add(p1,[p2x,0]) 00417 00418 def addy(p1,p2y): 00419 return add(p1,[0,p2y]) 00420 00421 #bar() 00422 #cross() 00423 #box() 00424 T() 00425