 |
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_HEX20_H 00021 #define MECHSYS_FEM_HEX20_H 00022 00023 // MechSys 00024 #include <mechsys/fem/node.h> 00025 #include <mechsys/fem/geomelem.h> 00026 00027 namespace FEM 00028 { 00029 00030 class Hex20 : public GeomElem 00031 { 00032 public: 00033 // Auxiliar structure to map local face IDs to local node IDs 00034 static const int Face2Node[6][8]; // 6 faces, 8 nodes/face 00035 00036 // Constructor 00037 Hex20 (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 Vertices Faces 00052 t 00053 | 4 15 7 00054 ,+--s @-------@--------@ +----------------+ 00055 r' ,'| ,'| ,'| ,'| 00056 12 @' | 14 ,' | ,' | ___ ,' | 00057 ,' |16 ,@ |19 ,' |,'5,' [0],' | 00058 5 ,' @ 6 ,' @ ,' |~~~ ,' | 00059 @'=======@=======@' | +'===============+' ,'| | 00060 | 13 | | | | ,'| | | |3| | 00061 | | | 11 | | |2| | | |,' | 00062 17 | 0 @- - - | @- - - -@ | |,' +- - - | +- - - -+ 00063 @ ,' @ ,' 3 | ,' | ,' 00064 | 8 @' 18 | ,' | ,' [1] ___| ,' 00065 | ,' | ,@ 10 | ,' ,'4,'| ,' 00066 | ,' | ,' | ,' ~~~ | ,' 00067 @-------@--------@' +----------------+' 00068 1 9 2 00069 */ 00070 const int Hex20::Face2Node[6][8] = {{ 0, 4, 7, 3, 16, 15, 19, 11 }, 00071 { 1, 2, 6, 5, 9, 18, 13, 17 }, 00072 { 0, 1, 5, 4, 8, 17, 12, 16 }, 00073 { 2, 3, 7, 6, 10, 19, 14, 18 }, 00074 { 0, 3, 2, 1, 11, 10, 9, 8 }, 00075 { 4, 5, 6, 7, 12, 13, 14, 15 }}; 00076 00077 00079 00080 00081 inline Hex20::Hex20 (int NDim) 00082 : GeomElem(NDim, /*NN*/20, /*NFN*/8, "Hex20") 00083 { 00084 SetIPs (27); 00085 } 00086 00087 inline void Hex20::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("Hex20::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 Hex20::Shape (double r, double s, double t) const 00101 { 00102 double rp1=1.0+r; double rm1=1.0-r; 00103 double sp1=1.0+s; double sm1=1.0-s; 00104 double tp1=1.0+t; double tm1=1.0-t; 00105 00106 N( 0) = 0.125*rm1*sm1*tm1*(-r-s-t-2); 00107 N( 1) = 0.125*rp1*sm1*tm1*( r-s-t-2); 00108 N( 2) = 0.125*rp1*sp1*tm1*( r+s-t-2); 00109 N( 3) = 0.125*rm1*sp1*tm1*(-r+s-t-2); 00110 N( 4) = 0.125*rm1*sm1*tp1*(-r-s+t-2); 00111 N( 5) = 0.125*rp1*sm1*tp1*( r-s+t-2); 00112 N( 6) = 0.125*rp1*sp1*tp1*( r+s+t-2); 00113 N( 7) = 0.125*rm1*sp1*tp1*(-r+s+t-2); 00114 N( 8) = 0.25*(1-r*r)*sm1*tm1; 00115 N( 9) = 0.25*rp1*(1-s*s)*tm1; 00116 N(10) = 0.25*(1-r*r)*sp1*tm1; 00117 N(11) = 0.25*rm1*(1-s*s)*tm1; 00118 N(12) = 0.25*(1-r*r)*sm1*tp1; 00119 N(13) = 0.25*rp1*(1-s*s)*tp1; 00120 N(14) = 0.25*(1-r*r)*sp1*tp1; 00121 N(15) = 0.25*rm1*(1-s*s)*tp1; 00122 N(16) = 0.25*rm1*sm1*(1-t*t); 00123 N(17) = 0.25*rp1*sm1*(1-t*t); 00124 N(18) = 0.25*rp1*sp1*(1-t*t); 00125 N(19) = 0.25*rm1*sp1*(1-t*t); 00126 } 00127 00128 inline void Hex20::Derivs(double r, double s, double t) const 00129 { 00130 /* _ _ T 00131 * | dNi | 00132 * dN = | --- | , where cj = r, s 00133 * |_ dcj _| 00134 * 00135 * dN(j,i), j=>local coordinate and i=>shape function 00136 */ 00137 double rp1=1.0+r; double rm1=1.0-r; 00138 double sp1=1.0+s; double sm1=1.0-s; 00139 double tp1=1.0+t; double tm1=1.0-t; 00140 00141 //Derivatives with respect to r 00142 dNdR(0, 0)= -0.125*sm1*tm1*(-r-s-t-2)-0.125*rm1*sm1*tm1; 00143 dNdR(0, 1)= 0.125*sm1*tm1*( r-s-t-2)+0.125*rp1*sm1*tm1; 00144 dNdR(0, 2)= 0.125*sp1*tm1*( r+s-t-2)+0.125*rp1*sp1*tm1; 00145 dNdR(0, 3)= -0.125*sp1*tm1*(-r+s-t-2)-0.125*rm1*sp1*tm1; 00146 dNdR(0, 4)= -0.125*sm1*tp1*(-r-s+t-2)-0.125*rm1*sm1*tp1; 00147 dNdR(0, 5)= 0.125*sm1*tp1*( r-s+t-2)+0.125*rp1*sm1*tp1; 00148 dNdR(0, 6)= 0.125*sp1*tp1*( r+s+t-2)+0.125*rp1*sp1*tp1; 00149 dNdR(0, 7)= -0.125*sp1*tp1*(-r+s+t-2)-0.125*rm1*sp1*tp1; 00150 dNdR(0, 8)= -0.5*r*sm1*tm1; 00151 dNdR(0, 9)= 0.25*(1-s*s)*tm1; 00152 dNdR(0,10)= -0.5*r*sp1*tm1; 00153 dNdR(0,11)= -0.25*(1-s*s)*tm1; 00154 dNdR(0,12)= -0.5*r*sm1*tp1; 00155 dNdR(0,13)= 0.25*(1-s*s)*tp1; 00156 dNdR(0,14)= -0.5*r*sp1 *tp1; 00157 dNdR(0,15)= -0.25*(1-s*s)*tp1; 00158 dNdR(0,16)= -0.25*sm1*(1-t*t); 00159 dNdR(0,17)= 0.25*sm1*(1-t*t); 00160 dNdR(0,18)= 0.25*sp1*(1-t*t); 00161 dNdR(0,19)= -0.25*sp1*(1-t*t); 00162 00163 //Derivatives with respect to s 00164 dNdR(1, 0)= -0.125*rm1*tm1*(-r-s-t-2)-0.125*rm1*sm1*tm1; 00165 dNdR(1, 1)= -0.125*rp1*tm1*( r-s-t-2)-0.125*rp1*sm1*tm1; 00166 dNdR(1, 2)= 0.125*rp1*tm1*( r+s-t-2)+0.125*rp1*sp1*tm1; 00167 dNdR(1, 3)= 0.125*rm1*tm1*(-r+s-t-2)+0.125*rm1*sp1*tm1; 00168 dNdR(1, 4)= -0.125*rm1*tp1*(-r-s+t-2)-0.125*rm1*sm1*tp1; 00169 dNdR(1, 5)= -0.125*rp1*tp1*( r-s+t-2)-0.125*rp1*sm1*tp1; 00170 dNdR(1, 6)= 0.125*rp1*tp1*( r+s+t-2)+0.125*rp1*sp1*tp1; 00171 dNdR(1, 7)= 0.125*rm1*tp1*(-r+s+t-2)+0.125*rm1*sp1*tp1; 00172 dNdR(1, 8)= -0.25*(1-r*r)*tm1; 00173 dNdR(1, 9)= -0.5*s*rp1*tm1; 00174 dNdR(1,10)= 0.25*(1-r*r)*tm1; 00175 dNdR(1,11)= -0.5*s*rm1*tm1; 00176 dNdR(1,12)= -0.25*(1-r*r)*tp1; 00177 dNdR(1,13)= -0.5*s*rp1*tp1; 00178 dNdR(1,14)= 0.25*(1-r*r)*tp1; 00179 dNdR(1,15)= -0.5*s*rm1*tp1; 00180 dNdR(1,16)= -0.25*rm1*(1-t*t); 00181 dNdR(1,17)= -0.25*rp1*(1-t*t); 00182 dNdR(1,18)= 0.25*rp1*(1-t*t); 00183 dNdR(1,19)= 0.25*rm1*(1-t*t); 00184 00185 //Derivatives with respect to t 00186 dNdR(2, 0)= -0.125*rm1*sm1*(-r-s-t-2)-0.125*rm1*sm1*tm1; 00187 dNdR(2, 1)= -0.125*rp1*sm1*( r-s-t-2)-0.125*rp1*sm1*tm1; 00188 dNdR(2, 2)= -0.125*rp1*sp1*( r+s-t-2)-0.125*rp1*sp1*tm1; 00189 dNdR(2, 3)= -0.125*rm1*sp1*(-r+s-t-2)-0.125*rm1*sp1*tm1; 00190 dNdR(2, 4)= 0.125*rm1*sm1*(-r-s+t-2)+0.125*rm1*sm1*tp1; 00191 dNdR(2, 5)= 0.125*rp1*sm1*( r-s+t-2)+0.125*rp1*sm1*tp1; 00192 dNdR(2, 6)= 0.125*rp1*sp1*( r+s+t-2)+0.125*rp1*sp1*tp1; 00193 dNdR(2, 7)= 0.125*rm1*sp1*(-r+s+t-2)+0.125*rm1*sp1*tp1; 00194 dNdR(2, 8)= -0.25*(1-r*r)*sm1; 00195 dNdR(2, 9)= -0.25*rp1*(1-s*s); 00196 dNdR(2,10)= -0.25*(1-r*r)*sp1; 00197 dNdR(2,11)= -0.25*rm1*(1-s*s); 00198 dNdR(2,12)= 0.25*(1-r*r)*sm1; 00199 dNdR(2,13)= 0.25*rp1*(1-s*s); 00200 dNdR(2,14)= 0.25*(1-r*r)*sp1; 00201 dNdR(2,15)= 0.25*rm1*(1-s*s); 00202 dNdR(2,16)= -0.5*t*rm1*sm1; 00203 dNdR(2,17)= -0.5*t*rp1*sm1; 00204 dNdR(2,18)= -0.5*t*rp1*sp1; 00205 dNdR(2,19)= -0.5*t*rm1*sp1; 00206 } 00207 00208 inline void Hex20::FaceShape(double r, double s) const 00209 { 00210 /* 3 6 2 00211 * @---------@----------@ 00212 * | (1,1)| 00213 * | s ^ | 00214 * | | | 00215 * | | | 00216 * 7 @ +----> r @ 5 00217 * | (0,0) | 00218 * | | 00219 * | | 00220 * |(-1,-1) | 00221 * @---------@----------@ 00222 * 0 4 1 00223 */ 00224 double rp1=1.0+r; double rm1=1.0-r; 00225 double sp1=1.0+s; double sm1=1.0-s; 00226 00227 FN(0) = 0.25*rm1*sm1*(rm1+sm1-3.0); 00228 FN(1) = 0.25*rp1*sm1*(rp1+sm1-3.0); 00229 FN(2) = 0.25*rp1*sp1*(rp1+sp1-3.0); 00230 FN(3) = 0.25*rm1*sp1*(rm1+sp1-3.0); 00231 FN(4) = 0.50*sm1*(1.0-r*r); 00232 FN(5) = 0.50*rp1*(1.0-s*s); 00233 FN(6) = 0.50*sp1*(1.0-r*r); 00234 FN(7) = 0.50*rm1*(1.0-s*s); 00235 } 00236 00237 inline void Hex20::FaceDerivs(double r, double s) const 00238 { 00239 /* _ _ T 00240 * | dNi | 00241 * FdN = | --- | , where cj = r, s 00242 * |_ dcj _| 00243 * 00244 * FdNdR(j,i), j=>local coordinate and i=>shape function 00245 */ 00246 double rp1=1.0+r; double rm1=1.0-r; 00247 double sp1=1.0+s; double sm1=1.0-s; 00248 00249 FdNdR(0,0) = - 0.25 * sm1 * (rm1 + rm1 + sm1 - 3.0); 00250 FdNdR(0,1) = 0.25 * sm1 * (rp1 + rp1 + sm1 - 3.0); 00251 FdNdR(0,2) = 0.25 * sp1 * (rp1 + rp1 + sp1 - 3.0); 00252 FdNdR(0,3) = - 0.25 * sp1 * (rm1 + rm1 + sp1 - 3.0); 00253 FdNdR(0,4) = - r * sm1; 00254 FdNdR(0,5) = 0.50 * (1.0 - s * s); 00255 FdNdR(0,6) = - r * sp1; 00256 FdNdR(0,7) = - 0.5 * (1.0 - s * s); 00257 00258 FdNdR(1,0) = - 0.25 * rm1 * (sm1 + rm1 + sm1 - 3.0); 00259 FdNdR(1,1) = - 0.25 * rp1 * (sm1 + rp1 + sm1 - 3.0); 00260 FdNdR(1,2) = 0.25 * rp1 * (sp1 + rp1 + sp1 - 3.0); 00261 FdNdR(1,3) = 0.25 * rm1 * (sp1 + rm1 + sp1 - 3.0); 00262 FdNdR(1,4) = - 0.50 * (1.0 - r * r); 00263 FdNdR(1,5) = - s * rp1; 00264 FdNdR(1,6) = 0.50 * (1.0 - r * r); 00265 FdNdR(1,7) = - s * rm1; 00266 } 00267 00268 inline void Hex20::NatCoords (Mat_t & C) const 00269 { 00270 C.change_dim(20,4); 00271 C = -1.0, -1.0, -1.0, 1.0, // nodes 0 to 7 00272 1.0, -1.0, -1.0, 1.0, 00273 1.0, 1.0, -1.0, 1.0, 00274 -1.0, 1.0, -1.0, 1.0, 00275 -1.0, -1.0, 1.0, 1.0, 00276 1.0, -1.0, 1.0, 1.0, 00277 1.0, 1.0, 1.0, 1.0, 00278 -1.0, 1.0, 1.0, 1.0, 00279 00280 0.0, -1.0, -1.0, 1.0, // nodes 8 to 11 00281 1.0, 0.0, -1.0, 1.0, 00282 0.0, 1.0, -1.0, 1.0, 00283 -1.0, 0.0, -1.0, 1.0, 00284 00285 0.0, -1.0, 1.0, 1.0, // nodes 12 to 15 00286 1.0, 0.0, 1.0, 1.0, 00287 0.0, 1.0, 1.0, 1.0, 00288 -1.0, 0.0, 1.0, 1.0, 00289 00290 -1.0, -1.0, 0.0, 1.0, // nodes 16 to 19 00291 1.0, -1.0, 0.0, 1.0, 00292 1.0, 1.0, 0.0, 1.0, 00293 -1.0, 1.0, 0.0, 1.0; 00294 } 00295 00296 00298 00299 00300 // Allocate a new element 00301 GeomElem * Hex20Maker (int NDim) { return new Hex20(NDim); } 00302 00303 // Register element 00304 int Hex20Register () 00305 { 00306 GeomElemFactory["Hex20"] = Hex20Maker; 00307 GEOM.Set ("Hex20", (double)GEOM.Keys.Size()); 00308 return 0; 00309 } 00310 00311 // Call register 00312 int __Hex20_dummy_int = Hex20Register(); 00313 00314 }; // namespace FEM 00315 00316 #endif // MECHSYS_FEM_HEX20_H
1.7.6.1
