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=C0103,R0902,R0904,R0914 00026 00027 from __future__ import division, print_function 00028 from numpy import dot,cross,array,matrix,zeros 00029 from numpy.linalg import solve 00030 00031 from pyNastran.bdf.cards.elements.elements import Element 00032 from pyNastran.general.generalMath import Area 00033 00034 def Volume4(n1, n2, n3, n4): 00035 """ 00036 V = (a-d) * ((b-d) x (c-d))/6 where x is cross and * is dot 00037 \f[ \large V = {(a-d) \dot \left( (b-d) \times (c-d) \right) }{6} \f] 00038 """ 00039 V = dot((n1-n4), cross(n2-n4, n3-n4))/6. 00040 return V 00041 00042 class SolidElement(Element): 00043 def __init__(self, card, data): 00044 Element.__init__(self, card, data) 00045 00046 def crossReference(self, model): 00047 self.nodes = model.Nodes(self.nodes) 00048 self.pid = model.Property(self.pid) 00049 00050 def Mass(self): 00051 return self.Rho()*self.Volume() 00052 00053 def Mid(self): 00054 return self.pid.Mid() 00055 00056 def Rho(self): 00057 try: 00058 return self.pid.mid.rho 00059 except AttributeError: 00060 print("self.pid = %s" %(self.pid)) 00061 print("self.pid.mid = %s" %(str(self.pid.mid))) 00062 #print("self.pid = %s" %(self.pid)) 00063 #print("self.pid = %s" %(self.pid)) 00064 raise 00065 00066 def isSameCard(self, elem, debug=False): 00067 if self.type!=elem.type: return False 00068 fields1 = [self.eid, self.Pid()]+self.nodes 00069 fields2 = [elem.eid, elem.Pid()]+elem.nodes 00070 if debug: 00071 print("fields1=%s fields2=%s" %(fields1, fields2)) 00072 return self.isSameFields(fields1, fields2) 00073 00074 def rawFields(self): 00075 fields = [self.type, self.eid, self.Pid()]+self.nodeIDs() 00076 return fields 00077 00078 class CHEXA8(SolidElement): 00079 """ 00080 CHEXA EID PID G1 G2 G3 G4 G5 G6 00081 G7 G8 00082 """ 00083 type = 'CHEXA' 00084 asterType = 'HEXA8' 00085 calculixType = 'C3D8' 00086 def __init__(self, card=None, data=None): 00087 SolidElement.__init__(self, card, data) 00088 if card: 00089 #print "hexa = ",card 00090 self.eid = card.field(1) 00091 self.pid = card.field(2) 00092 nids = card.fields(3, 11) 00093 else: 00094 self.eid = data[0] 00095 self.pid = data[1] 00096 nids = data[2:] 00097 assert len(data)==10, 'len(data)=%s data=%s' %(len(data),data) 00098 #print "nids = ",nids 00099 self.prepareNodeIDs(nids) 00100 assert len(self.nodes)==8 00101 00102 def areaCentroid(self, n1, n2, n3, n4): 00103 a = n1-n2 00104 b = n2-n4 00105 area1 = Area(a, b) 00106 c1 = (n1+n2+n4)/3. 00107 00108 a = n2-n4 00109 b = n2-n3 00110 area2 = Area(a, b) 00111 c2 = (n2+n3+n4)/3. 00112 00113 area = area1+area2 00114 centroid = (c1*area1+c2*area2)/area 00115 return(area, centroid) 00116 00117 def Centroid(self): 00118 (n1, n2, n3, n4, n5, n6, n7, n8) = self.nodePositions() 00119 A1,c1 = self.areaCentroid(n1, n2, n3, n4) 00120 A2,c2 = self.areaCentroid(n5, n6, n7, n8) 00121 c = (c1+c2)/2. 00122 return c 00123 00124 def Volume(self): 00125 (n1, n2, n3, n4, n5, n6, n7, n8) = self.nodePositions() 00126 (A1, c1) = self.areaCentroid(n1, n2, n3, n4) 00127 (A2, c2) = self.areaCentroid(n5, n6, n7, n8) 00128 V = (A1+A2)/2.*(c1-c2) 00129 return abs(V) 00130 00131 class CHEXA20(CHEXA8): 00132 """ 00133 CHEXA EID PID G1 G2 G3 G4 G5 G6 00134 G7 G8 G9 G10 G11 G12 G13 G14 00135 G15 G16 G17 G18 G19 G20 00136 """ 00137 type = 'CHEXA' 00138 asterType = 'HEXA20' 00139 calculixType = 'C3D20' 00140 def __init__(self, card=None, data=None): 00141 SolidElement.__init__(self, card, data) 00142 00143 if card: 00144 self.eid = card.field(1) 00145 self.pid = card.field(2) 00146 nids = card.fields(3, 23) 00147 else: 00148 self.eid = data[0] 00149 self.pid = data[1] 00150 nids = data[2:] 00151 self.prepareNodeIDs(nids, allowEmptyNodes=True) 00152 msg = 'len(nids)=%s nids=%s' %(len(nids), nids) 00153 assert len(self.nodes)<=20,msg 00154 00155 def Centroid(self): 00156 (n1, n2, n3, n4, n5, 00157 n6, n7, n8, n9, n10, 00158 n11, n12, n13, n14, n15, 00159 n16, n17, n18, n19, n20) = self.nodePositions() 00160 (A1, c1) = self.areaCentroid(n1, n2, n3, n4) 00161 (A2, c2) = self.areaCentroid(n5, n6, n7, n8) 00162 c = (c1+c2)/2. 00163 return c 00164 00165 def Volume(self): 00166 (n1, n2, n3, n4, n5, 00167 n6, n7, n8, n9, n10, 00168 n11, n12, n13, n14, n15, 00169 n16, n17, n18, n19, n20) = self.nodePositions() 00170 (A1, c1) = self.areaCentroid(n1, n2, n3, n4) 00171 (A2, c2) = self.areaCentroid(n5, n6, n7, n8) 00172 V = (A1+A2)/2.*(c1-c2) 00173 return abs(V) 00174 00175 00176 class CPENTA6(SolidElement): 00177 """ 00178 CPENTA EID PID G1 G2 G3 G4 G5 G6 00179 *----------* 00180 / \ / \ 00181 / A \ / c \ 00182 *---*-----*-----* 00183 V = (A1+A2)/2 * (c1-c2) 00184 C = (c1-c2)/2 00185 """ 00186 type = 'CPENTA' 00187 asterType = 'PENTA6' 00188 calculixType = 'C3D6' 00189 def __init__(self, card=None, data=None): 00190 SolidElement.__init__(self, card, data) 00191 00192 if card: 00193 self.eid = card.field(1) 00194 self.pid = card.field(2) 00195 nids = card.fields(3, 9) 00196 else: 00197 self.eid = data[0] 00198 self.pid = data[1] 00199 nids = data[2:] 00200 assert len(data)==8, 'len(data)=%s data=%s' %(len(data),data) 00201 self.prepareNodeIDs(nids) 00202 assert len(self.nodes)==6 00203 00204 def getFaceNodesAndArea(self, nid, nidOpposite): 00205 nids = self.nodeIDs()[:6] 00206 indx1 = nids.index(nid) 00207 indx2 = nids.index(nidOpposite) 00208 ## offset so it's easier to map the nodes with the QRG 00209 pack = [indx1+1, indx2+1] 00210 pack.sort() 00211 mapper = { # reverse points away from the element 00212 [1,2]: [1,2,3], # close 00213 [2,3]: [1,2,3], 00214 [1,3]: [1,2,3], 00215 00216 [4,5]: [4,5,6], # far-reverse 00217 [5,6]: [4,5,6], 00218 [4,6]: [4,5,6], 00219 00220 [1,5]: [1,2,5,4], # bottom 00221 [2,4]: [1,2,5,4], 00222 00223 [1,6]: [1,3,6,4], # left-reverse 00224 [3,4]: [1,3,6,4], 00225 00226 [2,6]: [2,5,6,3], # right 00227 [3,5]: [2,3,6,5], 00228 } 00229 pack2 = mapper[pack] 00230 if len(pack2)==3: 00231 (n1, n2, n3) = pack2 00232 faceNodeIDs = [n1, n2, n3] 00233 n1i = nids.index(n1-1) 00234 n2i = nids.index(n2-1) 00235 n3i = nids.index(n3-1) 00236 p1 = self.nodes[n1i].Position() 00237 p2 = self.nodes[n2i].Position() 00238 p3 = self.nodes[n3i].Position() 00239 a = p1-p2 00240 b = p1-p3 00241 A = Area(a, b) 00242 else: 00243 (n1, n2, n3, n4) = pack2 00244 n1i = nids.index(n1-1) 00245 n2i = nids.index(n2-1) 00246 n3i = nids.index(n3-1) 00247 n4i = nids.index(n4-1) 00248 faceNodeIDs = [n1, n2, n3, n4] 00249 p1 = self.nodes[n1i].Position() 00250 p2 = self.nodes[n2i].Position() 00251 p3 = self.nodes[n3i].Position() 00252 p4 = self.nodes[n4i].Position() 00253 a = p1-p3 00254 b = p2-p4 00255 A = Area(a, b) 00256 return [faceNodeIDs,A] 00257 00258 def Centroid(self): 00259 (n1, n2, n3, n4, n5, n6) = self.nodePositions() 00260 c1 = (n1+n2+n3)/3. 00261 c2 = (n4+n5+n6)/3. 00262 c = (c1+c2)/2. 00263 return c 00264 00265 def Volume(self): 00266 (n1, n2, n3, n4, n5, n6) = self.nodePositions() 00267 A1 = Area(n3-n1, n2-n1) 00268 A2 = Area(n6-n4, n5-n4) 00269 c1 = (n1+n2+n3)/3. 00270 c2 = (n4+n5+n6)/3. 00271 00272 V = (A1+A2)/2. * (c1-c2) 00273 return abs(V) 00274 00275 class CPENTA15(CPENTA6): 00276 """ 00277 CPENTA EID PID G1 G2 G3 G4 G5 G6 00278 G7 G8 G9 G10 G11 G12 G13 G14 00279 G15 00280 """ 00281 type = 'CPENTA' 00282 asterType = 'PENTA15' 00283 calculixType = 'C3D15' 00284 def __init__(self, card=None, data=None): 00285 SolidElement.__init__(self, card, data) 00286 00287 if card: 00288 self.eid = card.field(1) 00289 self.pid = card.field(2) 00290 nids = card.fields(3, 18) 00291 else: 00292 self.eid = data[0] 00293 self.pid = data[1] 00294 nids = data[2:] 00295 assert len(data)==17, 'len(data)=%s data=%s' %(len(data),data) 00296 self.prepareNodeIDs(nids, allowEmptyNodes=True) 00297 assert len(self.nodes)<=15 00298 00299 def Centroid(self): 00300 (n1, n2, n3, n4, n5, 00301 n6, n7, n8, n9, n10, 00302 n11, n12, n13, n14, n15) = self.nodePositions() 00303 c1 = (n1+n2+n3)/3. 00304 c2 = (n4+n5+n6)/3. 00305 c = (c1-c2)/2. 00306 return c 00307 00308 def Volume(self): 00309 (n1, n2, n3, n4, n5, 00310 n6, n7, n8, n9, n10, 00311 n11, n12, n13, n14, n15) = self.nodePositions() 00312 A1 = Area(n3-n1, n2-n1) 00313 A2 = Area(n6-n4, n5-n4) 00314 c1 = (n1+n2+n3)/3. 00315 c2 = (n4+n5+n6)/3. 00316 00317 V = (A1+A2)/2. * (c1-c2) 00318 return abs(V) 00319 00320 class CTETRA4(SolidElement): 00321 """ 00322 CTETRA EID PID G1 G2 G3 G4 00323 """ 00324 type = 'CTETRA' 00325 asterType = 'TETRA4' 00326 calculixType = 'C3D4' 00327 def __init__(self, card=None, data=None): 00328 SolidElement.__init__(self, card, data) 00329 if card: 00330 self.eid = card.field(1) 00331 self.pid = card.field(2) 00332 nids = card.fields(3, 7) 00333 else: 00334 self.eid = data[0] 00335 self.pid = data[1] 00336 nids = data[2:] 00337 assert len(data)==6, 'len(data)=%s data=%s' %(len(data),data) 00338 self.prepareNodeIDs(nids) 00339 assert len(self.nodes)==4 00340 00341 def Volume(self): 00342 (n1, n2, n3, n4) = self.nodePositions() 00343 return Volume4(n1, n2, n3, n4) 00344 00345 def Centroid(self): 00346 (n1, n2, n3, n4) = self.nodePositions() 00347 return (n1+n2+n3+n4)/4. 00348 00349 def getFaceNodes(self,nid,nidOpposite): 00350 nids = self.nodeIDs()[:4] 00351 indx = nids.index(nidOpposite) 00352 nids.pop(indx) 00353 return nids 00354 00355 def Stiffness(self): 00356 ng = 1 00357 (P,W) = gauss(1) 00358 00359 K = zeros((6,6)) 00360 for i in xrange(ng): 00361 for j in xrange(ng): 00362 for k in xrange(ng): 00363 p = [ P[i],P[j],P[k] ] 00364 w = W[i]*W[j]*W[k] 00365 pz = self.zeta(p) 00366 J = self.Jacobian(p) 00367 #N = self.N() 00368 #B = self.B() 00369 K += w*J*self.BtEB(pz) 00370 #K += w*J*B.T*E*B 00371 ### 00372 ### 00373 ### 00374 return K 00375 00376 def BtEB(self,pz): 00377 #E = self.Dsolid() 00378 00379 #o = E*e 00380 #e = [exx,eyy,ezz,2exy,2eyz,2ezx] 00381 #C = array([[Bxi*E11+Byi*E14+Bzi*E16, Byi*E12+Bxi*E14+Bzi*E15, Bzi*E13+Byi*E15+Bxi*E16] 00382 # [Bxi*E12+Byi*E24+Bzi*E26, Byi*E22+Bxi*E24+Bzi*E25, Bzi*E23+Byi*E25+Bxi*E26] 00383 # [Bxi*E13+Byi*E34+Bzi*E36, Byi*E23+Bxi*E34+Bzi*E35, Bzi*E33+Byi*E35+Bxi*E36] 00384 # [Bxi*E14+Byi*E44+Bzi*E46, Byi*E24+Bxi*E44+Bzi*E45, Bzi*E34+Byi*E45+Bxi*E46] 00385 # [Bxi*E15+Byi*E45+Bzi*E56, Byi*E25+Bxi*E45+Bzi*E55, Bzi*E35+Byi*E55+Bxi*E56] 00386 # [Bxi*E16+Byi*E46+Bzi*E16, Byi*E26+Bxi*E46+Bzi*E56, Bzi*E36+Byi*E56+Bxi*E66]]) 00387 00388 #Qij = array([[Bxj*C11+ Byj*C41+Bzj+Bzj*C61, Bxj*C12+Byj*C42+Bzj*C62, Bxj*C13+Byj*C43+Bzj*C63], 00389 # [Byj*C21+ Bxj*C41+Bzj+Bzj*C51, Byj*C22+Bxj*C42+Bzj*C52, Byj*C23+Bxj*C43+Bzj*C53], 00390 # [Bzj*C31+ Byj*C51+Bzj+Bxj*C61, Bzj*C32+Byj*C52+Bxj*C62, Bzj*C33+Byj*C53+Bxj*C63]]) 00391 Qij = None 00392 return Qij 00393 00394 def zeta(self,p): 00395 p2 = array([1,p[0],p[1],p[2]]) 00396 J = self.Jacobian() 00397 zeta = solve(J,p2) 00398 return zeta 00399 00400 def Jacobian(self): 00401 """ 00402 \f[ \large [J] = 00403 \left[ 00404 \begin{array}{ccc} 00405 1 & 1 & 1 \\ 00406 x_1 & y_1 & z_1 \\ 00407 x_2 & y_2 & z_2 \\ 00408 x_3 & y_3 & z_3 \\ 00409 x_4 & y_4 & z_4 \\ 00410 \end{array} \right] 00411 \f] 00412 @todo 00413 this has got to be wrong 00414 @warning 00415 this has got to be wrong 00416 """ 00417 m = matrix((6,6),'d') 00418 (n1,n2,n3,n4) = self.nodePositions() 00419 m[0][0] = m[0][1] = m[0][2] = m[0][2] = 1. 00420 m[1][0]=n1[0]; m[2][0]=n1[1]; m[3][0]=n1[2]; 00421 m[1][1]=n2[0]; m[2][1]=n2[1]; m[3][1]=n2[2]; 00422 m[1][2]=n3[0]; m[2][2]=n3[1]; m[3][2]=n3[2]; 00423 m[1][3]=n4[0]; m[2][3]=n4[1]; m[3][3]=n4[2]; 00424 return m 00425 00426 class CTETRA10(CTETRA4): 00427 """ 00428 CTETRA EID PID G1 G2 G3 G4 G5 G6 00429 G7 G8 G9 G10 00430 CTETRA 1 1 239 229 516 99 335 103 00431 265 334 101 102 00432 """ 00433 type = 'CTETRA' 00434 asterType = 'TETRA10' 00435 calculixType = 'C3D10' 00436 def __init__(self, card=None, data=None): 00437 SolidElement.__init__(self, card, data) 00438 if card: 00439 self.eid = card.field(1) 00440 self.pid = card.field(2) 00441 nids = card.fields(3, 13) 00442 else: 00443 self.eid = data[0] 00444 self.pid = data[1] 00445 nids = data[2:] 00446 assert len(data)==12, 'len(data)=%s data=%s' %(len(data), data) 00447 self.prepareNodeIDs(nids, allowEmptyNodes=True) 00448 assert len(self.nodes)<=10 00449 00450 def N_10(self, g1, g2, g3, g4): 00451 N1 = g1*(2*g1-1) 00452 N2 = g2*(2*g2-1) 00453 N3 = g3*(2*g3-1) 00454 N4 = g4*(2*g4-1) 00455 N5 = 4*g1*g2 00456 N6 = 4*g2*g3 00457 N7 = 4*g3*g1 00458 N8 = 4*g1*g4 00459 N9 = 4*g2*g4 00460 N10 = 4*g3*g4 00461 return (N1, N2, N3, N4, N5, N6, N7, N8, N9, N10) 00462 00463 def Volume(self): 00464 (n1, n2, n3, n4, n5, n6, n7, n8, n9, n10) = self.nodePositions() 00465 return Volume4(n1, n2, n3, n4) 00466 00467 def Centroid(self): 00468 (n1, n2, n3, n4, n5, n6, n7, n8, n9, n10) = self.nodePositions() 00469 return (n1+n2+n3+n4)/4. 00470 00471 def getFaceNodes(self, nid, nidOpposite): 00472 nids = self.nodeIDs()[:4] 00473 indx = nids.index(nidOpposite) 00474 nids.pop(indx) 00475 return nids