 |
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_TET10_H 00021 #define MECHSYS_FEM_TET10_H 00022 00023 // MechSys 00024 #include <mechsys/fem/node.h> 00025 #include <mechsys/fem/geomelem.h> 00026 00027 namespace FEM 00028 { 00029 00030 class Tet10 : public GeomElem 00031 { 00032 public: 00033 // Auxiliar structure to map local face IDs to local node IDs 00034 static const int Face2Node[4][6]; // 4 faces, 6 nodes/face 00035 00036 // Constructor 00037 Tet10 (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 /* TODO: check this 00050 const int Tet10::Face2Node[4][6] = {{ }, 00051 { }, 00052 { }, 00053 { }}; 00054 */ 00055 00057 00058 00059 inline Tet10::Tet10 (int NDim) 00060 : GeomElem(NDim, /*NN*/10, /*NFN*/6, "Tet10") 00061 { 00062 SetIPs (4); 00063 } 00064 00065 inline void Tet10::SetIPs (int TotNIP) 00066 { 00067 if (TotNIP== 4) IPs = TET_IP4; 00068 else throw new Fatal("Tet10::SetIPs: Total number of integration points = %d is invalid",TotNIP); 00069 00070 NIP = TotNIP; 00071 FIPs = TRI_IP6; 00072 NFIP = 6; 00073 } 00074 00075 inline void Tet10::Shape(double r, double s, double t) const 00076 { 00077 /* t 00078 * | 00079 * | 00080 * | 3 00081 * @, 00082 * /|` 00083 * || `, 00084 * / | ', 00085 * | | \ 00086 * / | `. 00087 * | | `, 9 00088 * / @ 7 `@ 00089 * | | \ 00090 * / | `. 00091 * | | ', 00092 * 8 @ | \ 00093 * | @.,,_ 6 `. 00094 * | / 0 ``'-.,,@_ `. 00095 * | / ``''-.,,_ ', 2 00096 * | / ``'@.,,, 00097 * | ' ,.-`` ``''- s 00098 * | ,@ 4 _,-'` 00099 * ' / ,.'` 00100 * | / _.@`` 00101 * '/ ,-'` 5 00102 * |/ ,.-`` 00103 * / _,-`` 00104 * .@ '` 00105 * / 1 00106 * / 00107 * / 00108 * r 00109 */ 00110 double u = 1.0 - r - s - t; 00111 00112 // corners 00113 N(0) = u*(2.0*u - 1.0); 00114 N(1) = r*(2.0*r - 1.0); 00115 N(2) = s*(2.0*s - 1.0); 00116 N(3) = t*(2.0*t - 1.0); 00117 00118 // midedge 00119 N(4) = 4.0 * u * r; 00120 N(5) = 4.0 * r * s; 00121 N(6) = 4.0 * s * u; 00122 N(7) = 4.0 * u * t; 00123 N(8) = 4.0 * r * t; 00124 N(9) = 4.0 * s * t; 00125 } 00126 00127 inline void Tet10::Derivs(double r, double s, double t) const 00128 { 00129 /* _ _ T 00130 * | dNi | 00131 * Derivs = | --- | , where cj = r, s, t 00132 * |_ dcj _| 00133 * 00134 * Derivs(j,i), j=>local coordinate and i=>shape function 00135 */ 00136 // r-derivatives: dN0/dr to dN9/dr 00137 dNdR(0,0) = 4.0*(r + s + t) - 3.0; 00138 dNdR(0,1) = 4.0*r - 1.0; 00139 dNdR(0,2) = 0.0; 00140 dNdR(0,3) = 0.0; 00141 dNdR(0,4) = 4.0 - 8.0*r - 4.0*s - 4.0*t; 00142 dNdR(0,5) = 4.0*s; 00143 dNdR(0,6) = -4.0*s; 00144 dNdR(0,7) = -4.0*t; 00145 dNdR(0,8) = 4.0*t; 00146 dNdR(0,9) = 0.0; 00147 00148 // s-derivatives: dN0/ds to dN9/ds 00149 dNdR(1,0) = 4.0*(r + s + t) - 3.0; 00150 dNdR(1,1) = 0.0; 00151 dNdR(1,2) = 4.0*s - 1.0; 00152 dNdR(1,3) = 0.0; 00153 dNdR(1,4) = -4.0*r; 00154 dNdR(1,5) = 4.0*r; 00155 dNdR(1,6) = 4.0 - 4.0*r - 8.0*s - 4.0*t; 00156 dNdR(1,7) = -4.0*t; 00157 dNdR(1,8) = 0.0; 00158 dNdR(1,9) = 4.0*t; 00159 00160 // t-derivatives: dN0/dt to dN9/dt 00161 dNdR(2,0) = 4.0*(r + s + t) - 3.0; 00162 dNdR(2,1) = 0.0; 00163 dNdR(2,2) = 0.0; 00164 dNdR(2,3) = 4.0*t - 1.0; 00165 dNdR(2,4) = -4.0*r; 00166 dNdR(2,5) = 0.0; 00167 dNdR(2,6) = -4.0*s; 00168 dNdR(2,7) = 4.0 - 4.0*r - 4.0*s - 8.0*t; 00169 dNdR(2,8) = 4.0*r; 00170 dNdR(2,9) = 4.0*s; 00171 } 00172 00173 inline void Tet10::FaceShape(double r, double s) const 00174 { 00175 /* s 00176 * ^ 00177 * | 00178 * 2 @. 00179 * | '. 00180 * | '. 4 00181 * 5 @ @. 00182 * | '. 00183 * | '. 00184 * @-----@-----@-> r 00185 * 0 3 1 00186 */ 00187 FN(0) = 1.0-(r+s)*(3.0-2.0*(r+s)); 00188 FN(1) = r*(2.0*r-1.0); 00189 FN(2) = s*(2.0*s-1.0); 00190 FN(3) = 4.0*r*(1.0-(r+s)); 00191 FN(4) = 4.0*r*s; 00192 FN(5) = 4.0*s*(1.0-(r+s)); 00193 } 00194 00195 inline void Tet10::FaceDerivs(double r, double s) const 00196 { 00197 /* _ _ T 00198 * | dNi | 00199 * Derivs = | --- | , where cj = r, s 00200 * |_ dcj _| 00201 * 00202 * Derivs(j,i), j=>local coordinate and i=>shape function 00203 */ 00204 FdNdR(0,0) = -3.0 + 4.0 * (r + s); FdNdR(1,0) = -3.0 + 4.0*(r + s); 00205 FdNdR(0,1) = 4.0 * r - 1.; FdNdR(1,1) = 0.0 ; 00206 FdNdR(0,2) = 0.0; FdNdR(1,2) = 4.0 * s - 1.0; 00207 FdNdR(0,3) = 4.0 - 8.0 * r - 4.0 * s; FdNdR(1,3) = -4.0 * r; 00208 FdNdR(0,4) = 4.0 * s; FdNdR(1,4) = 4.0 * r; 00209 FdNdR(0,5) = -4.0 * s; FdNdR(1,5) = 4.0 - 4.0 * r - 8.0*s; 00210 } 00211 00212 inline void Tet10::NatCoords (Mat_t & C) const 00213 { 00214 throw new Fatal("Tet10::NatCoords: method is not available yet"); 00215 } 00216 00217 00219 00220 00221 // Allocate a new element 00222 GeomElem * Tet10Maker (int NDim) { return new Tet10(NDim); } 00223 00224 // Register element 00225 int Tet10Register () 00226 { 00227 GeomElemFactory["Tet10"] = Tet10Maker; 00228 GEOM.Set ("Tet10", (double)GEOM.Keys.Size()); 00229 return 0; 00230 } 00231 00232 // Call register 00233 int __Tet10_dummy_int = Tet10Register(); 00234 00235 }; // namespace FEM 00236 00237 #endif // MECHSYS_FEM_TET10_H
1.7.6.1
