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 # pylint: disable=R0904,R0902,E1101,E1103 00026 00027 from __future__ import division, print_function 00028 #import sys 00029 00030 from numpy import matrix, zeros, ones, array, transpose, dot 00031 from numpy.linalg import norm 00032 #from pyNastran.general.generalMath import printMatrix 00033 00034 from pyNastran.bdf.errors import CardInstantiationError 00035 from pyNastran.bdf.fieldWriter import set_blank_if_default 00036 from pyNastran.bdf.cards.baseCard import Element, Mid 00037 00038 00039 class RodElement(Element): # CROD, CONROD, CTUBE 00040 def __init__(self, card, data): 00041 Element.__init__(self, card, data) 00042 00043 def crossReference(self, model): 00044 self.nodes = model.Nodes(self.nodes) 00045 self.pid = model.Property(self.pid) 00046 00047 def Rho(self): 00048 """returns the material density \f$ \rho \f$""" 00049 #print str(self.pid),type(self.pid) 00050 #raise NotImplementedMethodError('implement self.Rho() for %s' %(self.type)) 00051 return self.pid.mid.rho 00052 00053 def Length(self): 00054 """ 00055 Returns the length of the element 00056 \f[ \large \sqrt{ (n_{x2}-n_{x1})^2+(n_{y2}-n_{y1})^2+(n_{z2}-n_{z1})^2 } \f] 00057 @param self the object pointer 00058 """ 00059 L = norm(self.nodes[1].Position()-self.nodes[0].Position()) 00060 return L 00061 00062 def Mass(self): 00063 """ 00064 returns the mass of the element 00065 \f[ \large mass = \left( \rho A + nsm \right) L \f] 00066 """ 00067 L = self.Length() 00068 mass = (self.Rho()*self.Area() + self.Nsm())*L 00069 return mass 00070 00071 00072 class LineElement(Element): # CBAR, CBEAM, CBEAM3, CBEND 00073 def __init__(self, card, data): 00074 Element.__init__(self, card, data) 00075 00076 def C(self): 00077 """torsional constant""" 00078 return self.pid.C() 00079 00080 def Area(self): 00081 """returns the area of the element face""" 00082 raise NotImplementedError('implement self.Area() for %s' %(self.type)) 00083 00084 def E(self): 00085 """returns the Young's Modulus \f$ E \f$""" 00086 return self.pid.mid.E() 00087 00088 def G(self): 00089 """returns the Shear Modulus \f$ G \f$""" 00090 return self.pid.mid.G() 00091 00092 def J(self): 00093 """returns the Polar Moment of Inertia. \f$ J \f$""" 00094 return self.pid.J() 00095 00096 def I11(self): 00097 """returns the Moment of Inertia. \f$ I_{11} \f$""" 00098 return self.pid.I11() 00099 00100 def I22(self): 00101 """returns the Moment of Inertia. \f$ I_{22} \f$""" 00102 return self.pid.I22() 00103 00104 def I12(self): 00105 """returns the Moment of Inertia. \f$ I_{12} \f$""" 00106 return self.pid.I12() 00107 00108 def Nu(self): 00109 """returns Poisson's Ratio \f$ \nu \f$""" 00110 return self.pid.mid.nu 00111 00112 def Rho(self): 00113 """returns the material density \f$ \rho \f$""" 00114 #print str(self.pid),type(self.pid) 00115 #raise NotImplementedMethodError('implement self.Rho() for %s' %(self.type)) 00116 return self.pid.mid.rho 00117 00118 def Nsm(self): 00119 """Placeholder method for the non-structural mass""" 00120 raise NotImplementedError('implement self.Area() for %s' %(self.type)) 00121 00122 def MassPerLength(self): 00123 """Returns the mass per unit length""" 00124 return self.pid.MassPerLength() 00125 00126 def Mass(self): 00127 """ 00128 returns the mass of the element 00129 00130 \f[ \large mass = \left( \rho A + nsm \right) L \f] 00131 """ 00132 L = self.Length() 00133 mass = (self.Rho()*self.Area() + self.Nsm())*L 00134 return mass 00135 00136 def crossReference(self, model): 00137 self.nodes = model.Nodes(self.nodes) 00138 self.pid = model.Property(self.pid) 00139 00140 def Length(self): 00141 """ 00142 Returns the length of the element 00143 \f[ \large \sqrt{ (n_{x2}-n_{x1})^2+(n_{y2}-n_{y1})^2+(n_{z2}-n_{z1})^2 } \f] 00144 @param self the object pointer 00145 """ 00146 L = norm(self.nodes[1].Position()-self.nodes[0].Position()) 00147 return L 00148 00149 def k_Axial(self): 00150 """ 00151 Returns the axial stiffness matrix. 00152 00153 \f[ \large k_{Axial} = \frac{AE}{2L} 00154 \left[ 00155 \begin{array}{cc} 00156 1 & -1 \\ 00157 -1 & 1 00158 \end{array} \right] 00159 \f] 00160 """ 00161 raise NotImplementedError() 00162 L = self.Length() 00163 E = self.E() 00164 A = self.Area() 00165 kMag = A*E/(2*L) 00166 K = ones(1,1) 00167 K[0,1] = K[1,0] = -1 00168 return kMag*K 00169 00170 def k_Torsion(self): # not done 00171 """ 00172 Returns the torsional stiffness matrix. 00173 00174 \f[ \large k_{Axial} = \frac{L}{GJ} 00175 \left[ 00176 \begin{array}{cc} 00177 1 & -1 \\ 00178 -1 & 1 00179 \end{array} \right] 00180 \f] 00181 @warning formula not verified 00182 """ 00183 raise NotImplementedError() 00184 L = self.Length() 00185 G = self.G() 00186 J = self.J() 00187 #A = self.Area() 00188 #kMag = A*E/(2*L) 00189 kMag = L/(G*J) 00190 K = ones(1, 1) 00191 K[0, 1] = K[1, 0] = -1 00192 return kMag*K 00193 00194 def k_Bending(self): 00195 """ 00196 Returns the bending stiffness matrix. 00197 00198 \f[ \large k_{Bending} = \frac{EI}{L^3} 00199 \left[ 00200 \begin{array}{cccc} 00201 12 & 6L & -12 & 6L \\ 00202 6L & 4L^2 & -6L & 2L^2 \\ 00203 -12 & -6L & 12 & -6L \\ 00204 6L & 2L^2 & -6L & 4L^2 00205 \end{array} \right] 00206 \f] 00207 """ 00208 raise NotImplementedError() 00209 L = self.Length() 00210 E = self.E() 00211 I = self.I() 00212 LL = L*L 00213 LLL = L*LL 00214 sL = 6*L 00215 tLL = 2*LL 00216 fLL = 4*LL 00217 kMag = E*I/LLL 00218 00219 #K = Matrix(zeros(4,4)) 00220 K = matrix([[12., sL, -12, sL], 00221 [sL, fLL, -sL, tLL], 00222 [-12, -sL, 12., -sL], 00223 [sL, tLL, -sL, fLL]]) 00224 #M[1,0] = sL 00225 #M[2,0] = -12. 00226 #M[3,0] = sL 00227 00228 #M[2,4] = -sL 00229 #M[1,1] = M[3,3] = fLL 00230 00231 return kMag*K 00232 00233 class CROD(RodElement): 00234 type = 'CROD' 00235 def __init__(self, card=None, data=None): 00236 RodElement.__init__(self, card, data) 00237 if card: 00238 self.eid = int(card.field(1)) 00239 self.pid = int(card.field(2, self.eid)) 00240 nids = card.fields(3, 5) 00241 else: 00242 self.eid = data[0] 00243 self.pid = data[1] 00244 nids = data[2:4] 00245 ### 00246 self.prepareNodeIDs(nids) 00247 assert len(self.nodes) == 2 00248 00249 def Centroid(self): 00250 return (self.nodes[0].Position()+self.nodes[1].Position())/2. 00251 00252 def Mid(self): 00253 return self.pid.Mid() 00254 00255 def Area(self): 00256 return self.pid.A 00257 00258 def Nsm(self): 00259 return self.pid.nsm 00260 00261 def MassPerLength(self): 00262 massPerLength = self.pid.mid.rho*self.pid.A + self.pid.nsm 00263 return massPerLength 00264 00265 def Rmatrix(self, model, is3D): 00266 """ 00267 where \f$ [R]_{ij} \f$ is the tranformation matrix 00268 \f[ \large [R]_{ij} \left[ 00269 \begin{array}{ccc} 00270 g_x \dot e_x & g_x \dot e_y & g_x \dot e_z \\ 00271 g_y \dot e_x & g_y \dot e_y & g_y \dot e_z \\ 00272 g_z \dot e_x & g_z \dot e_y & g_z \dot e_z 00273 \end{array} \right] 00274 \f] 00275 """ 00276 (n1, n2) = self.nodeIDs() 00277 p1 = model.Node(n1).xyz 00278 p2 = model.Node(n2).xyz 00279 v1 = p2-p1 00280 v1 = v1/norm(v1) 00281 (l,m,n) = v1 00282 00283 v1x = array([v1[0], 0., 0.]) 00284 v1y = array([ 0., v1[1], 0.]) 00285 v1z = array([ 0., 0., v1[2]]) 00286 00287 g1x = array([1., 0., 0.]) 00288 g1y = array([0., 1., 0.]) 00289 g1z = array([0., 0., 1.]) 00290 00291 if is3D: 00292 R = matrix([ #global rod 00293 [dot(v1x,g1x), dot(v1y,g1x), dot(v1z,g1x)], 00294 [dot(v1x,g1y), dot(v1y,g1y), dot(v1z,g1y)], 00295 [dot(v1x,g1z), dot(v1y,g1z), dot(v1z,g1z)], 00296 ]) # rod 00297 #R = matrix([ 00298 # [], 00299 # [], 00300 # [], 00301 # ]) 00302 00303 else: 00304 R = matrix([ # there can be no z component 00305 [dot(v1x,g1x), dot(v1y,g1x)], 00306 [dot(v1x,g1y), dot(v1y,g1y)], 00307 ]) # rod 00308 #R = matrix([ # there can be no z component 00309 # [dot(v1x,g1x),dot(v1y,g1x)], 00310 # [dot(v1x,g1y),dot(v1y,g1y)], 00311 # ]) # rod 00312 return R 00313 00314 def Lambda(self, model, debug=True): 00315 """ 00316 2d [l,m,0,0] 00317 [0,0,l,m] 00318 00319 3d [l,m,n,0,0,0] 00320 [0,0,0,l,m,n] 00321 """ 00322 is3D = False 00323 #R = self.Rmatrix(model,is3D) 00324 00325 (n1, n2) = self.nodeIDs() 00326 p1 = model.Node(n1).Position() 00327 p2 = model.Node(n2).Position() 00328 v1 = p2-p1 00329 if debug: 00330 print("v1=%s" %(v1)) 00331 v1 = v1/norm(v1) 00332 (l, m, n) = v1 00333 if is3D: 00334 Lambda = matrix(zeros((2,6),'d')) # 3D 00335 else: 00336 Lambda = matrix(zeros((2,4),'d')) 00337 00338 #print("R = \n",R) 00339 Lambda[0,0] = Lambda[1,2] = l 00340 Lambda[0,1] = Lambda[1,3] = m 00341 00342 if is3D: 00343 Lambda[0,2] = Lambda[1,5] = n # 3D 00344 if debug: 00345 print("Lambda = \n"+str(Lambda)) 00346 return Lambda 00347 00348 #def Stiffness1(self,model): 00349 #nIJV = [(nodes[0],1),(nodes[0],2),(nodes[1],1),] 00350 00351 def Stiffness(self, model, grav, is3D): # CROD/CONROD 00352 print("----------------") 00353 Lambda = self.Lambda(model) 00354 #print("Lambda = \n"+str(Lambda)) 00355 00356 k = self.Stiffness1D(model) #/250000. 00357 #print R 00358 #print(k) 00359 #print("Lambda.shape = ",Lambda.shape) 00360 K = dot(dot(transpose(Lambda), k), Lambda) 00361 #Fg = dot(dot(transpose(Lambda),grav),Lambda) 00362 00363 #print('size(grav) =',grav.shape) 00364 mass = self.Mass() 00365 mg = -grav*mass 00366 if is3D: 00367 Fg = [mg[0],mg[1],mg[2],mg[0],mg[1],mg[2]] 00368 else: 00369 Fg = [mg[0],mg[1],mg[0],mg[1]] 00370 #print("Fg = ",Fg) 00371 print("mass = ",mass) 00372 print("Fg = ",Fg) 00373 #print(K) 00374 print("K[%s] = \n%s\n" %(self.eid,K)) 00375 #print("Fg[%s] = %s\n" %(self.eid,Fg)) 00376 00377 nodes = self.nodeIDs() 00378 nIJV = [(nodes[0],1),(nodes[0],2),(nodes[1],1),(nodes[1],2),] 00379 nGrav = [(nodes[0],1),(nodes[0],2),(nodes[0],3),(nodes[1],1),(nodes[1],2),(nodes[1],3)] 00380 00381 return(K, nIJV, Fg, nGrav) 00382 00383 def displacementStressF06(self, model, q, dofs): 00384 pass 00385 00386 def displacementStress(self, model, q, dofs): 00387 (n1, n2) = self.nodeIDs() 00388 Lambda = self.Lambda(model, debug=False) 00389 n11 = dofs[(n1,1)] 00390 n12 = dofs[(n1,2)] 00391 n21 = dofs[(n2,1)] 00392 n22 = dofs[(n2,2)] 00393 print("type=%s n1=%s n2=%s" %(self.type, n1, n2)) 00394 #print("n11=%s n12=%s n21=%s n22=%s" %(n11,n12,n21,n22)) 00395 00396 q2 = array( [q[n11], q[n12], q[n21], q[n22]] ) 00397 print("q[%s] = %s" %(self.eid,q2)) 00398 #print("Lambda = \n"+str(Lambda)) 00399 00400 #print "Lsize = ",Lambda.shape 00401 #print "qsize = ",q.shape 00402 u = dot(array(Lambda), q2) 00403 #L = self.Length() 00404 EL = self.E()/self.Length() 00405 00406 #stressX = -EL*u[0]+EL*u[1] 00407 stressX = EL*(-u[0]+u[1]) 00408 print("stressX = %s [psi]\n" %(stressX)) 00409 return stressX 00410 00411 def Stiffness1D(self, model): # CROD/CONROD 00412 """ 00413 @todo remove this method after making sure it still works 00414 """ 00415 #L = norm(r) 00416 (n1,n2) = self.nodeIDs() 00417 #node1 = model.Node(n1) 00418 #node2 = model.Node(n2) 00419 00420 p1 = model.Node(n1).xyz 00421 p2 = model.Node(n2).xyz 00422 #print "node1 = ",node1 00423 #print "node2 = ",node2 00424 L = norm(p1-p2) 00425 00426 if L==0.0: 00427 msg = 'invalid CROD length=0.0\n%s' %(self.__repr__()) 00428 raise ZeroDivisionError(msg) 00429 00430 A = self.Area() 00431 #mat = self.mid 00432 E = self.E() 00433 print("A = ",A) 00434 print("E = ",E) 00435 print("L = ",L) 00436 #ki = 1. 00437 ki = A*E/L 00438 #ki = 250000. 00439 #knorm = 250000. 00440 K = ki*matrix([[1.,-1.],[-1.,1.]]) # rod 00441 00442 print("A=%g E=%g L=%g AE/L=%g" %(A,E,L,A*E/L)) 00443 #print "K = \n",K 00444 return K 00445 00446 def rawFields(self): 00447 fields = ['CROD',self.eid,self.Pid()]+self.nodeIDs() 00448 return fields 00449 00450 def reprFields(self): 00451 return self.rawFields() 00452 00453 class CTUBE(RodElement): 00454 type = 'CTUBE' 00455 def __init__(self, card=None, data=None): 00456 RodElement.__init__(self, card, data) 00457 if card: 00458 self.eid = int(card.field(1)) 00459 self.pid = int(card.field(2, self.eid)) 00460 nids = card.fields(3, 5) 00461 else: 00462 self.eid = data[0] 00463 self.pid = data[1] 00464 nids = data[2:4] 00465 ### 00466 self.prepareNodeIDs(nids) 00467 assert len(self.nodes) == 2 00468 00469 def Area(self): 00470 return self.pid.Area() 00471 00472 def Centroid(self): 00473 """@todo improve the formuala for CTUBE centroid""" 00474 return (self.nodes[0].Position()+self.nodes[1].Position())/2. 00475 00476 def rawFields(self): 00477 fields = ['CTUBE', self.eid, self.Pid()]+self.nodeIDs() 00478 return fields 00479 ### 00480 00481 class CONROD(RodElement): 00482 type = 'CONROD' 00483 def __init__(self, card=None, data=None): 00484 RodElement.__init__(self, card, data) 00485 if card: 00486 self.eid = int(card.field(1)) 00487 #print "self.eid = ",self.eid 00488 nids = card.fields(2,4) 00489 00490 self.mid = int(card.field(4)) 00491 self.A = float(card.field(5)) 00492 self.j = float(card.field(6, 0.0)) 00493 self.c = float(card.field(7, 0.0)) 00494 self.nsm = float(card.field(8, 0.0)) 00495 else: 00496 self.eid = data[0] 00497 nids = data[1:3] 00498 self.mid = data[3] 00499 self.A = data[4] 00500 self.j = data[5] 00501 self.c = data[6] 00502 self.nsm = data[7] 00503 ### 00504 self.prepareNodeIDs(nids) 00505 assert len(self.nodes) == 2 00506 #print self.nodes 00507 00508 def crossReference(self, model): 00509 self.nodes = model.Nodes(self.nodes) 00510 self.mid = model.Material(self.mid) 00511 00512 def Centroid(self): 00513 return (self.nodes[0].Position()+self.nodes[1].Position())/2. 00514 00515 def Mid(self): 00516 return Mid(self) 00517 00518 def Pid(self): 00519 return None 00520 00521 def MassPerLength(self): 00522 massPerLength = self.mid.rho*self.A + self.nsm 00523 return massPerLength 00524 00525 def C(self): 00526 """torsional constant""" 00527 return self.c 00528 00529 def Area(self): 00530 return self.A 00531 00532 def J(self): 00533 """returns the Polar Moment of Inertia. \f$ J \f$""" 00534 return self.j 00535 00536 def Nsm(self): 00537 """Placeholder method for the non-structural mass""" 00538 return self.nsm 00539 00540 def E(self): 00541 """returns the Young's Modulus \f$ E \f$""" 00542 return self.mid.E() 00543 00544 def G(self): 00545 """returns the Shear Modulus \f$ G \f$""" 00546 return self.mid.G() 00547 00548 def Nu(self): 00549 """returns Poisson's Ratio \f$ \nu \f$""" 00550 return self.mid.nu 00551 00552 def Rho(self): 00553 """returns the material density \f$ \rho \f$""" 00554 return self.mid.rho 00555 00556 def writeCodeAster(self): 00557 msg = '' 00558 msg += " POUTRE=_F(GROUP_MA='CONROD_%s',\n" %(self.eid) 00559 msg += " SECTION='CERCLE', # circular section\n" 00560 if self.Thickness(): 00561 msg += " CARA=('R','EP'), # radius, thickness\n" 00562 msg += " VALE=(%g,%g),\n" %(self.Radius(),self.Thickness()) 00563 else: 00564 msg += " CARA=('R') # radius\n" 00565 msg += " VALE=(%g),\n" %(self.Radius()) 00566 ### 00567 return msg 00568 00569 def rawFields(self): 00570 fields = ['CONROD', self.eid]+self.nodeIDs()+[self.Mid(), self.A, self.j, self.c, self.nsm] 00571 return fields 00572 00573 def reprFields(self): 00574 j = set_blank_if_default(self.j, 0.0) 00575 c = set_blank_if_default(self.c, 0.0) 00576 nsm = set_blank_if_default(self.nsm, 0.0) 00577 fields = ['CONROD', self.eid]+self.nodeIDs()+[self.Mid(), self.A, j, c, nsm] 00578 return fields 00579 00580 class CBAR(LineElement): 00581 """ 00582 CBAR EID PID GA GB X1 X2 X3 OFFT 00583 PA PB W1A W2A W3A W1B W2B W3B 00584 or 00585 CBAR EID PID GA GB G0 OFFT 00586 PA PB W1A W2A W3A W1B W2B W3B 00587 00588 CBAR 22062 4 21798 21799 0.0+0 0.0+0 -1. 00589 0.0+0 0.0+0 -9. 0.0+0 0.0+0 -9. 00590 """ 00591 type = 'CBAR' 00592 asterType = 'CBAR' 00593 def __init__(self, card=None, data=None): 00594 LineElement.__init__(self, card, data) 00595 if card: 00596 self.eid = int(card.field(1)) 00597 self.pid = int(card.field(2, self.eid)) 00598 self.ga = int(card.field(3)) 00599 self.gb = int(card.field(4)) 00600 self.initX_G0(card) 00601 00602 self.offt = card.field(8, 'GGG') 00603 #print 'self.offt = |%s|' %(self.offt) 00604 00605 self.pa = card.field(9, 0) 00606 self.pb = card.field(10, 0) 00607 00608 self.w1a = float(card.field(11, 0.0)) 00609 self.w2a = float(card.field(12, 0.0)) 00610 self.w3a = float(card.field(13, 0.0)) 00611 00612 self.w1b = float(card.field(14, 0.0)) 00613 self.w2b = float(card.field(15, 0.0)) 00614 self.w3b = float(card.field(16, 0.0)) 00615 else: ## @todo verify 00616 #data = [[eid,pid,ga,gb,pa,pb,w1a,w2a,w3a,w1b,w2b,w3b],[f,g0]] 00617 #data = [[eid,pid,ga,gb,pa,pb,w1a,w2a,w3a,w1b,w2b,w3b],[f,x1,x2,x3]] 00618 00619 main = data[0] 00620 00621 flag = data[1][0] 00622 if flag in [0, 1]: 00623 self.g0 = None 00624 self.x1 = data[1][1] 00625 self.x2 = data[1][2] 00626 self.x3 = data[1][3] 00627 else: 00628 self.g0 = data[1][1] 00629 self.x1 = None 00630 self.x2 = None 00631 self.x3 = None 00632 00633 self.eid = main[0] 00634 self.pid = main[1] 00635 self.ga = main[2] 00636 self.gb = main[3] 00637 #self.offt = str(data[4]) # GGG 00638 self.offt = 'GGG'## @todo offt can be an integer; translate to char 00639 self.pa = main[4] 00640 self.pb = main[5] 00641 00642 self.w1a = main[6] 00643 self.w2a = main[7] 00644 self.w3a = main[8] 00645 00646 self.w1b = main[9] 00647 self.w2b = main[10] 00648 self.w3b = main[11] 00649 ### 00650 #print("offt = %s" %(self.offt)) 00651 if not isinstance(self.offt, str): 00652 raise SyntaxError('invalid offt expected a string of length 3 ' 00653 'offt=|%s|' %(self.offt)) 00654 00655 msg = 'invalid offt parameter of %s...offt=%s' %(self.type, self.offt) 00656 # B,G,O 00657 assert self.offt[0] in ['G', 'B'], msg 00658 assert self.offt[1] in ['G', 'O'], msg 00659 assert self.offt[2] in ['G', 'O'], msg 00660 00661 def Mid(self): 00662 return self.pid.Mid() 00663 00664 def Area(self): 00665 A = self.pid.Area() 00666 assert isinstance(A, float) 00667 return A 00668 00669 def Length(self): 00670 L = norm(self.gb.Position()-self.ga.Position()) 00671 assert isinstance(L, float) 00672 return L 00673 00674 def Nsm(self): 00675 nsm = self.pid.Nsm() 00676 assert isinstance(nsm, float) 00677 return nsm 00678 00679 def I1(self): 00680 return self.pid.I1() 00681 00682 def I2(self): 00683 return self.pid.I2() 00684 00685 def Centroid(self): 00686 return (self.ga.Position()+self.gb.Position())/2. 00687 00688 def initX_G0(self, card): 00689 field5 = card.field(5) 00690 if isinstance(field5, int): 00691 self.g0 = field5 00692 self.x1 = None 00693 self.x2 = None 00694 self.x3 = None 00695 elif isinstance(field5, float): 00696 self.g0 = None 00697 self.x1 = float(card.field(5, 0.0)) 00698 self.x2 = float(card.field(6, 0.0)) 00699 self.x3 = float(card.field(7, 0.0)) 00700 else: 00701 #msg = 'field5 on %s is the wrong type...id=%s field5=%s type=%s' %(self.type,self.eid,field5,type(field5)) 00702 #raise InvalidFieldError(msg) 00703 self.g0 = None 00704 self.x1 = 0. 00705 self.x2 = 0. 00706 self.x3 = 0. 00707 #if self.eid==14100238: 00708 #print "g0=%s x1=%s x2=%s x3=%s" %(self.g0,self.x1,self.x2,self.x3) 00709 00710 def crossReference(self,mesh): 00711 """ 00712 set g0-ga to x1,x2,x3 00713 """ 00714 #if self.g0: 00715 # v = nodes[self.g0].Position()-nodes[self.ga].Position() 00716 # self.x1 = v[0] 00717 # self.x2 = v[1] 00718 # self.x3 = v[2] 00719 ### 00720 self.ga = mesh.Node(self.ga) 00721 self.gb = mesh.Node(self.gb) 00722 self.pid = mesh.Property(self.pid) 00723 00724 #def updateNodes(self,nodes): 00725 # """@todo maybe improve""" 00726 # self.crossReference(self,nodes) 00727 00728 def Ga(self): 00729 if isinstance(self.ga,int): 00730 return self.ga 00731 else: 00732 return self.ga.nid 00733 ### 00734 00735 def Gb(self): 00736 if isinstance(self.gb,int): 00737 return self.gb 00738 else: 00739 return self.gb.nid 00740 ### 00741 00742 def getX_G0_defaults(self): 00743 if self.g0: 00744 return (self.g0, None, None) 00745 else: 00746 #x1 = set_blank_if_default(self.x1, 0.0) 00747 #x2 = set_blank_if_default(self.x2, 0.0) 00748 #x3 = set_blank_if_default(self.x3, 0.0) 00749 return (self.x1, self.x2, self.x3) 00750 ### 00751 00752 def nodeIDs(self): 00753 return [self.Ga(),self.Gb()] 00754 00755 def Stiffness(self, model, grav, is3D): # CBAR 00756 print("----------------") 00757 Lambda = self.Lambda(model) 00758 #print("Lambda = \n"+str(Lambda)) 00759 00760 k = self.Stiffness1D(model) #/250000. 00761 #print R 00762 #print(k) 00763 print("Lambda.shape = ",Lambda.shape) 00764 K = dot(dot(transpose(Lambda), k), Lambda) 00765 #Fg = dot(dot(transpose(Lambda),grav),Lambda) 00766 00767 #print('size(grav) =',grav.shape) 00768 mass = self.Mass() 00769 mg = -grav*mass 00770 if is3D: 00771 Fg = [mg[0], mg[1], mg[2], mg[0], mg[1], mg[2]] 00772 else: 00773 Fg = [mg[0], mg[1], mg[0], mg[1]] 00774 #print("Fg = ",Fg) 00775 print("mass = ",mass) 00776 print("Fg = ",Fg) 00777 #print(K) 00778 print("K[%s] = \n%s\n" %(self.eid,K)) 00779 #print("Fg[%s] = %s\n" %(self.eid,Fg)) 00780 00781 nodes = self.nodeIDs() 00782 # u1 v1 theta1 u2,v2,w2 00783 nIJV = [ 00784 (nodes[0],1),(nodes[0],2),(nodes[0],5), (nodes[1],1),(nodes[1],2),(nodes[1],5), # X1 00785 (nodes[0],1),(nodes[0],2),(nodes[0],5), (nodes[1],1),(nodes[1],2),(nodes[1],5), # Y1 00786 (nodes[0],1),(nodes[0],2),(nodes[0],5), (nodes[1],1),(nodes[1],2),(nodes[1],5), # M1 00787 00788 (nodes[0],1),(nodes[0],2),(nodes[0],5), (nodes[1],1),(nodes[1],2),(nodes[1],5), # X2 00789 (nodes[0],1),(nodes[0],2),(nodes[0],5), (nodes[1],1),(nodes[1],2),(nodes[1],5), # Y2 00790 (nodes[0],1),(nodes[0],2),(nodes[0],5), (nodes[1],1),(nodes[1],2),(nodes[1],5), # M2 00791 ] 00792 #print("nIJV = ",nIJV) 00793 nGrav = [(nodes[0],1),(nodes[0],2),(nodes[0],3),(nodes[1],1),(nodes[1],2),(nodes[1],3)] 00794 00795 return(K, nIJV, Fg, nGrav) 00796 00797 def Stiffness1D(self, model): # CBAR 00798 #L = norm(r) 00799 (n1, n2) = self.nodeIDs() 00800 #node1 = model.Node(n1) 00801 #node2 = model.Node(n2) 00802 00803 p1 = model.Node(n1).xyz 00804 p2 = model.Node(n2).xyz 00805 #print "node1 = ",node1 00806 #print "node2 = ",node2 00807 L = norm(p1-p2) 00808 00809 if L==0.0: 00810 msg = 'invalid CBAR length=0.0\n%s' %(self.__repr__()) 00811 raise ZeroDivisionError(msg) 00812 00813 A = self.Area() 00814 #mat = self.mid 00815 E = self.E() 00816 I1 = self.I1() 00817 I2 = self.I2() 00818 print("A = ", A) 00819 print("E = ", E) 00820 print("L = ", L) 00821 print("I1 = ", I1) 00822 print("I2 = ", I2) 00823 #ki = 1. 00824 #ki = A*E/L 00825 #ki = 250000. 00826 #knorm = 250000. 00827 #K = ki*matrix([[1.,-1.],[-1.,1.]]) # rod 00828 #P = 10.416666667 00829 P = 1. 00830 k1 = A*E/L 00831 k2 = P*L**3/(E*I1) 00832 print("A=%g E=%g L=%g I1=%g I2=%G AE/L=%g L^3/E*Iz=%g" %(A,E,L,I1,I2,k1,k2)) 00833 #print "K = \n",K 00834 #return K 00835 00836 00837 L = self.Length() 00838 E = self.E() 00839 #I = self.I1() 00840 00841 LL = L*L 00842 #LLL = L*LL 00843 #kMag = E*I/LLL 00844 sL = 6*L 00845 tLL = 2*LL 00846 fLL = 4*LL 00847 00848 #K = zeros(4,4)) 00849 K = array([[0.,12., sL, 0.,-12, sL], 00850 [0.,sL, fLL, 0.,-sL, tLL], 00851 [0., 0, 0, 0., 0, 0], 00852 [0.,-12, -sL, 0.,12., -sL], 00853 [0.,sL, tLL, 0.,-sL, fLL], 00854 [0., 0, 0, 0., 0, 0], 00855 ]) 00856 print("k =\n"+str(K)) 00857 return K 00858 00859 def Lambda(self, model, debug=True): # CBAR from CROD/CONROD 00860 """ 00861 2d [l,m,0,0, l,m,0,0] 00862 [0,0,l,m, 0,0,l,m] L*k = 2x4*W4x4 00863 00864 3d [l,m,n,0,0,0, l,m,n,0,0,0] 00865 [0,0,0,l,m,n, 0,0,0,l,m,n] 00866 """ 00867 is3D = False 00868 #R = self.Rmatrix(model,is3D) 00869 00870 (n1,n2) = self.nodeIDs() 00871 p1 = model.Node(n1).Position() 00872 p2 = model.Node(n2).Position() 00873 v1 = p2-p1 00874 if debug: 00875 print("v1=%s" % (v1)) 00876 v1 = v1/norm(v1) 00877 #v2 = v2/norm(v2) 00878 #v3 = v3/norm(v3) 00879 (l, m, n) = v1 00880 #(o,p,q) = v2 00881 #(r,s,t) = v3 00882 print("l=%s m=%s n=%s" % (l, m, n)) 00883 if is3D: 00884 Lambda = matrix([[l,m,n,0,0,0,], 00885 [m,l,n,0,0,0,], 00886 [0,0,1,0,0,0,], 00887 [0,0,0,l,m,n,], 00888 [0,0,0,m,l,n,], 00889 [0,0,0,0,0,1,], 00890 ]) 00891 else: 00892 Lambda = matrix([[l,m,n,0,0,0,], 00893 [m,l,n,0,0,0,], 00894 [0,0,1,0,0,0,], 00895 [0,0,0,l,m,n,], 00896 [0,0,0,m,l,n,], 00897 [0,0,0,0,0,1,], 00898 ]) 00899 if debug: 00900 print("Lambda* = \n"+str(Lambda)) 00901 return Lambda 00902 00903 def rawFields(self): 00904 """@todo not perfectly accurate""" 00905 (x1, x2, x3) = self.getX_G0_defaults() 00906 offt = set_blank_if_default(self.offt, 'GGG') 00907 fields = ['CBAR',self.eid,self.Pid(),self.Ga(),self.Gb(),x1,x2,x3,offt, 00908 self.pa,self.pb,self.w1a,self.w2a,self.w3a,self.w1b,self.w2b,self.w3b] 00909 00910 return fields 00911 00912 def reprFields(self): 00913 pa = set_blank_if_default(self.pa, 0) 00914 pb = set_blank_if_default(self.pb, 0) 00915 00916 w1a = set_blank_if_default(self.w1a, 0.0) 00917 w2a = set_blank_if_default(self.w2a, 0.0) 00918 w3a = set_blank_if_default(self.w3a, 0.0) 00919 w1b = set_blank_if_default(self.w1b, 0.0) 00920 w2b = set_blank_if_default(self.w2b, 0.0) 00921 w3b = set_blank_if_default(self.w3b, 0.0) 00922 (x1, x2, x3) = self.getX_G0_defaults() 00923 offt = set_blank_if_default(self.offt,'GGG') 00924 fields = ['CBAR',self.eid,self.Pid(),self.Ga(),self.Gb(),x1,x2,x3,offt, 00925 pa,pb,w1a,w2a,w3a,w1b,w2b,w3b] 00926 00927 return fields 00928 00929 class CBEAM3(CBAR): 00930 """ 00931 Defines a three-node beam element 00932 """ 00933 type = 'CBEAM3' 00934 def __init__(self, card=None, data=None): 00935 LineElement.__init__(self, card, data) 00936 if card: 00937 self.eid = int(card.field(1)) 00938 self.eid = int(card.field(1)) 00939 self.pid = int(card.field(2, self.eid)) 00940 self.ga = int(card.field(3)) 00941 self.gb = int(card.field(4)) 00942 self.gc = int(card.field(5)) 00943 00944 self.initX_G0(card) 00945 00946 self.w1a = float(card.field(9, 0.0)) 00947 self.w2a = float(card.field(10, 0.0)) 00948 self.w3a = float(card.field(11, 0.0)) 00949 00950 self.w1b = float(card.field(12, 0.0)) 00951 self.w2b = float(card.field(13, 0.0)) 00952 self.w3b = float(card.field(14, 0.0)) 00953 00954 self.w1c = float(card.field(15, 0.0)) 00955 self.w2c = float(card.field(16, 0.0)) 00956 self.w3c = float(card.field(17, 0.0)) 00957 00958 self.twa = card.field(18, 0.) 00959 self.twb = card.field(19, 0.) 00960 self.twc = card.field(20, 0.) 00961 00962 self.sa = card.field(21) 00963 self.sb = card.field(22) 00964 self.sc = card.field(23) 00965 else: 00966 raise NotImplementedError(data) 00967 ### 00968 00969 def crossReference(self,model): 00970 self.ga = model.Node(self.ga) 00971 self.gb = model.Node(self.gb) 00972 self.gc = model.Node(self.gc) 00973 self.pid = model.Property(self.pid) 00974 00975 def Length(self): 00976 """ 00977 L = gb-ga 00978 @todo improve formula 00979 """ 00980 L = norm(self.gb.Position()-self.ga.Position()) 00981 assert isinstance(L, float) 00982 return L 00983 00984 def rawFields(self): 00985 (x1,x2,x3) = self.getX_G0_defaults() 00986 (ga,gb,gc) = self.nodeIDs() 00987 fields = ['CBEAM3',self.eid,self.Pid(),ga,gb,gc,x1,x2,x3, 00988 self.w1a,self.w2a,self.w3a, self.w1b,self.w2b,self.w3b, self.w1c,self.w2c,self.w3c, 00989 self.twa,self.twb,self.twc,self.sa,self.sb,self.sc] 00990 return fields 00991 00992 def reprFields(self): 00993 w1a = set_blank_if_default(self.w1a, 0.0) 00994 w2a = set_blank_if_default(self.w2a, 0.0) 00995 w3a = set_blank_if_default(self.w3a, 0.0) 00996 w1b = set_blank_if_default(self.w1b, 0.0) 00997 w2b = set_blank_if_default(self.w2b, 0.0) 00998 w3b = set_blank_if_default(self.w3b, 0.0) 00999 w1c = set_blank_if_default(self.w1c, 0.0) 01000 w2c = set_blank_if_default(self.w2c, 0.0) 01001 w3c = set_blank_if_default(self.w3c, 0.0) 01002 01003 twa = set_blank_if_default(self.twa, 0.0) 01004 twb = set_blank_if_default(self.twb, 0.0) 01005 twc = set_blank_if_default(self.twc, 0.0) 01006 01007 (x1, x2, x3) = self.getX_G0_defaults() 01008 (ga, gb, gc) = self.nodeIDs() 01009 fields = ['CBEAM3', self.eid, self.Pid(), ga, gb, x1, x2, x3, 01010 w1a, w2a, w3a, w1b, w2b, w3b, w1c, w2c, w3c, 01011 twa, twb, twc, self.sa, self.sb, self.sc] 01012 return fields 01013 01014 01015 class CBEAM(CBAR): 01016 """ 01017 CBEAM EID PID GA GB X1 X2 X3 OFFT/BIT 01018 PA PB W1A W2A W3A W1B W2B W3B 01019 SA SB 01020 or 01021 CBEAM EID PID GA GB G0 - - OFFT/BIT 01022 PA PB W1A W2A W3A W1B W2B W3B 01023 SA SB 01024 """ 01025 type = 'CBEAM' 01026 def __init__(self, card=None, data=None): 01027 LineElement.__init__(self, card, data) 01028 if card: 01029 self.eid = int(card.field(1)) 01030 self.pid = int(card.field(2, self.eid)) 01031 self.ga = int(card.field(3)) 01032 self.gb = int(card.field(4)) 01033 01034 self.initX_G0(card) 01035 self.initOfftBit(card) 01036 self.pa = card.field(9) 01037 self.pb = card.field(10) 01038 01039 self.w1a = float(card.field(11, 0.0)) 01040 self.w2a = float(card.field(12, 0.0)) 01041 self.w3a = float(card.field(13, 0.0)) 01042 01043 self.w1b = float(card.field(14, 0.0)) 01044 self.w2b = float(card.field(15, 0.0)) 01045 self.w3b = float(card.field(16, 0.0)) 01046 01047 self.sa = card.field(17) 01048 self.sb = card.field(18) 01049 01050 else: ## @todo verify 01051 #data = [[eid,pid,ga,gb,sa,sb, pa,pb,w1a,w2a,w3a,w1b,w2b,w3b],[f,g0]] 01052 #data = [[eid,pid,ga,gb,sa,sb, pa,pb,w1a,w2a,w3a,w1b,w2b,w3b],[f,x1,x2,x3]] 01053 01054 main = data[0] 01055 01056 flag = data[1][0] 01057 if flag in [0,1]: 01058 self.g0 = None 01059 self.x1 = data[1][1] 01060 self.x2 = data[1][2] 01061 self.x3 = data[1][3] 01062 else: 01063 self.g0 = data[1][1] 01064 self.x1 = None 01065 self.x2 = None 01066 self.x3 = None 01067 01068 self.eid = main[0] 01069 self.pid = main[1] 01070 self.ga = main[2] 01071 self.gb = main[3] 01072 self.sa = main[4] 01073 self.sb = main[5] 01074 01075 self.isOfft = True ## @todo is this correct??? 01076 #self.offt = str(data[6]) # GGG 01077 self.offt = 'GGG' ## @todo is this correct??? 01078 01079 self.pa = main[6] 01080 self.pb = main[7] 01081 01082 self.w1a = main[8] 01083 self.w2a = main[9] 01084 self.w3a = main[10] 01085 01086 self.w1b = main[11] 01087 self.w2b = main[12] 01088 self.w3b = main[13] 01089 ### 01090 01091 def initOfftBit(self,card): 01092 field8 = card.field(8) 01093 if isinstance(field8, float): 01094 self.isOfft = False 01095 self.offt = None 01096 self.bit = field8 01097 elif field8 is None: 01098 self.isOfft = True 01099 self.offt = 'GGG' # default 01100 self.bit = None 01101 elif isinstance(field8, str): 01102 self.isOfft = True 01103 self.bit = None 01104 self.offt = field8 01105 #print "self.offt = ",self.offt 01106 assert self.offt[0] in ['G', 'B', 'O'], 'invalid offt parameter of CBEAM...offt=%s' %(self.offt) 01107 assert self.offt[1] in ['G', 'B', 'O'], 'invalid offt parameter of CBEAM...offt=%s' %(self.offt) 01108 assert self.offt[2] in ['G', 'B', 'O'], 'invalid offt parameter of CBEAM...offt=%s' %(self.offt) 01109 else: 01110 msg = 'field8 on %s card is not a string(offt) or bit (float)...field8=%s\n' %(self.type, field8) 01111 raise CardInstantiationError(msg) 01112 ### 01113 01114 def Mid(self): 01115 return self.pid.Mid() 01116 01117 def Area(self): 01118 return self.pid.Area() 01119 01120 def Nsm(self): 01121 #print "CBEAM pid = ",str(self.pid) 01122 return self.pid.Nsm() 01123 01124 def getOfft_Bit_defaults(self): 01125 if self.isOfft: 01126 field8 = self.offt 01127 else: 01128 field8 = set_blank_if_default(self.bit, 0.0) 01129 return field8 01130 01131 def crossReference(self,model): 01132 self.ga = model.Node(self.ga) 01133 self.gb = model.Node(self.gb) 01134 self.pid = model.Property(self.pid) 01135 01136 def Stiffness(self,model, r, A, E, I): # CBEAM 01137 """ 01138 from makeTruss??? 01139 """ 01140 Ke = zeros((6,6),'d') 01141 L = r 01142 AE = A*E 01143 EI = E*I 01144 01145 if 1: 01146 Ke[0, 0] = AE/L 01147 Ke[3, 0] = -AE/L 01148 Ke[0, 3] = -AE/L 01149 Ke[3, 3] = AE/L 01150 01151 Ke[1, 1] = 12*EI/L**3 01152 Ke[1, 2] = 6*EI/L**2 01153 Ke[2, 1] = Ke[1, 2] # 6*EI/L**2 01154 Ke[2, 2] = 4*EI/L 01155 01156 Ke[1, 4] = -Ke[1, 1] # -12*EI/L**3 01157 Ke[1, 5] = Ke[1, 2] # 6*EI/L**2 01158 Ke[2, 4] = -Ke[1, 2] # -6*EI/L**2 01159 Ke[2, 5] = 2*EI/L 01160 01161 Ke[4, 1] = -Ke[1, 4] # -12*EI/L**3 01162 Ke[4, 2] = Ke[2, 4] # -6*EI/L**2 01163 Ke[5, 1] = Ke[1, 2] # 6*EI/L**2 01164 Ke[5, 2] = Ke[2, 5] # 2*EI/L 01165 01166 Ke[4, 4] = Ke[1, 1] # 12*EI/L**3 01167 Ke[4, 5] = Ke[2, 4] # -6*EI/L**2 01168 Ke[5, 4] = Ke[2, 4] # -6*EI/L**2 01169 Ke[5, 5] = Ke[2, 2] # 4*EI/L 01170 else: 01171 Ke[0, 0] = AE 01172 Ke[3, 0] = -AE 01173 Ke[0, 3] = -AE 01174 Ke[3, 3] = AE 01175 01176 Ke[1, 1] = 12*EI/L**2 01177 Ke[1, 2] = 6*EI/L 01178 Ke[2, 1] = Ke[1, 2] # 6*EI/L**2 01179 Ke[2, 2] = 4*EI 01180 01181 Ke[1, 4] = -Ke[1, 1] # -12*EI/L**3 01182 Ke[1, 5] = Ke[1, 2] # 6*EI/L**2 01183 Ke[2, 4] = -Ke[1, 2] # -6*EI/L**2 01184 Ke[2, 5] = 2*EI 01185 01186 Ke[4, 1] = -Ke[1, 4] # -12*EI/L**3 01187 Ke[4, 2] = Ke[2, 4] # -6*EI/L**2 01188 Ke[5, 1] = Ke[1, 2] # 6*EI/L**2 01189 Ke[5, 2] = Ke[2, 5] # 2*EI/L 01190 01191 Ke[4, 4] = Ke[1, 1] # 12*EI/L**3 01192 Ke[4, 5] = Ke[2, 4] # -6*EI/L**2 01193 Ke[5, 4] = Ke[2, 4] # -6*EI/L**2 01194 Ke[5, 5] = Ke[2, 2] # 4*EI/L 01195 01196 Ke = Ke/L 01197 ### 01198 return Ke 01199 01200 def rawFields(self): 01201 (x1,x2,x3) = self.getX_G0_defaults() 01202 offt = self.getOfft_Bit_defaults() 01203 ga,gb = self.nodeIDs() 01204 fields = ['CBEAM', self.eid, self.Pid(), ga, gb, x1, x2, x3, offt, 01205 self.pa, self.pb, self.w1a, self.w2a, self.w3a, 01206 self.w1b, self.w2b, self.w3b, self.sa, self.sb] 01207 return fields 01208 01209 def reprFields(self): 01210 w1a = set_blank_if_default(self.w1a, 0.0) 01211 w2a = set_blank_if_default(self.w2a, 0.0) 01212 w3a = set_blank_if_default(self.w3a, 0.0) 01213 w1b = set_blank_if_default(self.w1b, 0.0) 01214 w2b = set_blank_if_default(self.w2b, 0.0) 01215 w3b = set_blank_if_default(self.w3b, 0.0) 01216 (x1, x2, x3) = self.getX_G0_defaults() 01217 offt = self.getOfft_Bit_defaults() 01218 ga,gb = self.nodeIDs() 01219 fields = ['CBEAM', self.eid, self.Pid(), ga, gb, x1, x2, x3, offt, 01220 self.pa, self.pb, w1a, w2a, w3a, 01221 w1b, w2b, w3b, self.sa, self.sb] 01222 return fields 01223 01224 class CBEND(LineElement): 01225 type = 'CBEND' 01226 def __init__(self, card=None, data=None): 01227 LineElement.__init__(self, card, data) 01228 self.eid = int(card.field(1)) 01229 self.pid = card.field(2,self.eid) 01230 self.ga = card.field(3) 01231 self.gb = card.field(4) 01232 x1Go = card.field(5) 01233 #self.initX_G0(card) 01234 if isinstance(x1Go,int): 01235 self.g0 = x1Go 01236 self.x1 = None 01237 self.x2 = None 01238 self.x3 = None 01239 elif isinstance(x1Go,float): 01240 self.g0 = None 01241 self.x1 = x1Go 01242 self.x2 = card.field(6) 01243 self.x3 = card.field(7) 01244 else: 01245 raise ValueError('invalid x1Go=|%s| on CBEND' % (x1Go)) 01246 self.geom = card.field(8) 01247 assert self.geom in [1, 2, 3, 4], 'geom is invalid geom=|%s|' % (self.geom) 01248 01249 def Area(self): 01250 return self.pid.Area() 01251 01252 def rawFields(self): 01253 (x1, x2, x3) = self.getX_G0_defaults() 01254 fields = ['CBEND', self.eid, self.Pid(), self.Ga(), self.Gb(), 01255 x1, x2, x3, self.geom] 01256 return fields 01257 01258 def reprFields(self): 01259 return self.rawFields()