pyNastran  0.5.0
pyNastran BDF Reader/Writer, OP2 Parser, and GUI
solid.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 # 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
 All Classes Namespaces Files Functions Variables