 |
MechSys
1.0
Computing library for simulations in continuum and discrete mechanics
|
|
MechSys
1.0
Computing library for simulations in continuum and discrete mechanics
|
00001 /************************************************************************ 00002 * MechSys - Open Library for Mechanical Systems * 00003 * Copyright (C) 2005 Dorival M. Pedroso, Raul Durand * 00004 * Copyright (C) 2009 Sergio Galindo * 00005 * * 00006 * This program is free software: you can redistribute it and/or modify * 00007 * it under the terms of the GNU General Public License as published by * 00008 * the Free Software Foundation, either version 3 of the License, or * 00009 * any later version. * 00010 * * 00011 * This program is distributed in the hope that it will be useful, * 00012 * but WITHOUT ANY WARRANTY; without even the implied warranty of * 00013 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * 00014 * GNU General Public License for more details. * 00015 * * 00016 * You should have received a copy of the GNU General Public License * 00017 * along with this program. If not, see <http://www.gnu.org/licenses/> * 00018 ************************************************************************/ 00019 00020 #ifndef MECHSYS_FEM_HEX8_H 00021 #define MECHSYS_FEM_HEX8_H 00022 00023 // MechSys 00024 #include <mechsys/fem/node.h> 00025 #include <mechsys/fem/geomelem.h> 00026 00027 namespace FEM 00028 { 00029 00030 class Hex8 : public GeomElem 00031 { 00032 public: 00033 // Auxiliar structure to map local face IDs to local node IDs 00034 static const int Face2Node[6][4]; // 6 faces, 4 nodes/face 00035 00036 // Constructor 00037 Hex8 (int NDim); 00038 00039 // Derived methods 00040 void SetIPs (int TotNIP); 00041 size_t FNode (size_t IdxFace, size_t IdxFNode) const { return Face2Node[IdxFace][IdxFNode]; } 00042 void Shape (double r, double s, double t) const; 00043 void Derivs (double r, double s, double t) const; 00044 void FaceShape (double r, double s) const; 00045 void FaceDerivs (double r, double s) const; 00046 void NatCoords (Mat_t & C) const; 00047 }; 00048 00049 00050 /* Local IDs 00051 Nodes Faces 00052 z 00053 | 4 7 00054 ,+--y @________________@ +________________+ 00055 x' ,'| ,'| ,'| ,'| 00056 ,' | ,' | ,' | ___ ,' | 00057 ,' | ,' | ,' |,'5,' [0],' | 00058 5 ,' | 6 ,' | ,' |~~~ ,' | 00059 @'===============@' | +'===============+' ,'| | 00060 | | | | | ,'| | | |3| | 00061 | | | | | |2| | | |,' | 00062 | 0 @______|_________@ | |,' +______|_________+ 00063 | ,' | ,' 3 | ,' | ,' 00064 | ,' | ,' | ,' [1] ___| ,' 00065 | ,' | ,' | ,' ,'4,'| ,' 00066 | ,' | ,' | ,' ~~~ | ,' 00067 @________________@' +________________+' 00068 1 2 00069 */ 00070 const int Hex8::Face2Node[6][4] = {{ 0, 4, 7, 3 }, 00071 { 1, 2, 6, 5 }, 00072 { 0, 1, 5, 4 }, 00073 { 2, 3, 7, 6 }, 00074 { 0, 3, 2, 1 }, 00075 { 4, 5, 6, 7 }}; 00076 00077 00079 00080 00081 inline Hex8::Hex8 (int NDim) 00082 : GeomElem(NDim, /*NN*/8, /*NFN*/4, "Hex8") 00083 { 00084 SetIPs (8); 00085 } 00086 00087 inline void Hex8::SetIPs (int TotNIP) 00088 { 00089 if (TotNIP== 8) IPs = HEX_IP2; 00090 else if (TotNIP== 27) IPs = HEX_IP3; 00091 else if (TotNIP== 64) IPs = HEX_IP4; 00092 else if (TotNIP==125) IPs = HEX_IP5; 00093 else throw new Fatal("Hex8::SetIPs: Total number of integration points = %d is invalid",TotNIP); 00094 00095 NIP = TotNIP; 00096 FIPs = QUAD_IP2; 00097 NFIP = 4; 00098 } 00099 00100 inline void Hex8::Shape (double r, double s, double t) const 00101 { 00102 /* t 00103 * ^ 00104 * | 00105 * 4 7 00106 * @________________@ 00107 * ,'| ,'| 00108 * ,' | ,' | 00109 * ,' | ,' | 00110 * 5 ,' | 6 ,' | 00111 * @'===============@' | 00112 * | | | | 00113 * | | | | 00114 * | 0 @_____ | ________@ 3 --> s 00115 * | ,' | ,' 00116 * | ,' | ,' 00117 * | ,' | ,' 00118 * | ,' | ,' 00119 * @________________@' 00120 * 1 2 00121 * ,' 00122 * |_ 00123 * r 00124 */ 00125 N(0) = 0.125*(1.0-r-s+r*s-t+s*t+r*t-r*s*t); 00126 N(1) = 0.125*(1.0+r-s-r*s-t+s*t-r*t+r*s*t); 00127 N(2) = 0.125*(1.0+r+s+r*s-t-s*t-r*t-r*s*t); 00128 N(3) = 0.125*(1.0-r+s-r*s-t-s*t+r*t+r*s*t); 00129 N(4) = 0.125*(1.0-r-s+r*s+t-s*t-r*t+r*s*t); 00130 N(5) = 0.125*(1.0+r-s-r*s+t-s*t+r*t-r*s*t); 00131 N(6) = 0.125*(1.0+r+s+r*s+t+s*t+r*t+r*s*t); 00132 N(7) = 0.125*(1.0-r+s-r*s+t+s*t-r*t-r*s*t); 00133 } 00134 00135 inline void Hex8::Derivs (double r, double s, double t) const 00136 { 00137 /* _ _ T 00138 * | dNi | 00139 * dN = | --- | , where cj = r, s 00140 * |_ dcj _| 00141 * 00142 * dN(j,i), j=>local coordinate and i=>shape function 00143 */ 00144 dNdR(0,0) = 0.125*(-1.0+s+t-s*t); dNdR(1,0)=0.125*(-1.0+r+t-r*t); dNdR(2,0)=0.125*(-1.0+r+s-r*s); 00145 dNdR(0,1) = 0.125*(+1.0-s-t+s*t); dNdR(1,1)=0.125*(-1.0-r+t+r*t); dNdR(2,1)=0.125*(-1.0-r+s+r*s); 00146 dNdR(0,2) = 0.125*(+1.0+s-t-s*t); dNdR(1,2)=0.125*(+1.0+r-t-r*t); dNdR(2,2)=0.125*(-1.0-r-s-r*s); 00147 dNdR(0,3) = 0.125*(-1.0-s+t+s*t); dNdR(1,3)=0.125*(+1.0-r-t+r*t); dNdR(2,3)=0.125*(-1.0+r-s+r*s); 00148 dNdR(0,4) = 0.125*(-1.0+s-t+s*t); dNdR(1,4)=0.125*(-1.0+r-t+r*t); dNdR(2,4)=0.125*(+1.0-r-s+r*s); 00149 dNdR(0,5) = 0.125*(+1.0-s+t-s*t); dNdR(1,5)=0.125*(-1.0-r-t-r*t); dNdR(2,5)=0.125*(+1.0+r-s-r*s); 00150 dNdR(0,6) = 0.125*(+1.0+s+t+s*t); dNdR(1,6)=0.125*(+1.0+r+t+r*t); dNdR(2,6)=0.125*(+1.0+r+s+r*s); 00151 dNdR(0,7) = 0.125*(-1.0-s-t-s*t); dNdR(1,7)=0.125*(+1.0-r+t-r*t); dNdR(2,7)=0.125*(+1.0-r+s-r*s); 00152 } 00153 00154 inline void Hex8::FaceShape (double r, double s) const 00155 { 00156 /* 3 2 00157 * @--------------------@ 00158 * | (1,1)| 00159 * | s ^ | 00160 * | | | 00161 * | | | 00162 * | +----> r | 00163 * | (0,0) | 00164 * | | 00165 * | | 00166 * |(-1,-1) | 00167 * @--------------------@ 00168 * 0 1 00169 */ 00170 FN(0) = 0.25*(1.0-r-s+r*s); 00171 FN(1) = 0.25*(1.0+r-s-r*s); 00172 FN(2) = 0.25*(1.0+r+s+r*s); 00173 FN(3) = 0.25*(1.0-r+s-r*s); 00174 } 00175 00176 inline void Hex8::FaceDerivs (double r, double s) const 00177 { 00178 /* _ _ T 00179 * | dNi | 00180 * FdNdR = | --- | , where Rj = r, s 00181 * |_ dRj _| 00182 * 00183 * FdNdR(j,i), j=>local coordinate and i=>shape function 00184 */ 00185 00186 FdNdR(0,0) = 0.25*(-1.0+s); FdNdR(1,0) = 0.25*(-1.0+r); 00187 FdNdR(0,1) = 0.25*(+1.0-s); FdNdR(1,1) = 0.25*(-1.0-r); 00188 FdNdR(0,2) = 0.25*(+1.0+s); FdNdR(1,2) = 0.25*(+1.0+r); 00189 FdNdR(0,3) = 0.25*(-1.0-s); FdNdR(1,3) = 0.25*(+1.0-r); 00190 } 00191 00192 inline void Hex8::NatCoords (Mat_t & C) const 00193 { 00194 C.change_dim(8,4); 00195 C = -1.0, -1.0, -1.0, 1.0, 00196 +1.0, -1.0, -1.0, 1.0, 00197 +1.0, +1.0, -1.0, 1.0, 00198 -1.0, +1.0, -1.0, 1.0, 00199 -1.0, -1.0, +1.0, 1.0, 00200 +1.0, -1.0, +1.0, 1.0, 00201 +1.0, +1.0, +1.0, 1.0, 00202 -1.0, +1.0, +1.0, 1.0; 00203 } 00204 00205 00207 00208 00209 // Allocate a new element 00210 GeomElem * Hex8Maker (int NDim) { return new Hex8(NDim); } 00211 00212 // Register element 00213 int Hex8Register () 00214 { 00215 GeomElemFactory["Hex8"] = Hex8Maker; 00216 GEOM.Set ("Hex8", (double)GEOM.Keys.Size()); 00217 return 0; 00218 } 00219 00220 // Call register 00221 int __Hex8_dummy_int = Hex8Register(); 00222 00223 }; // namespace FEM 00224 00225 #endif // MECHSYS_FEM_HEX8_H
1.7.6.1
