 |
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_TRI15_H 00021 #define MECHSYS_FEM_TRI15_H 00022 00023 // MechSys 00024 #include <mechsys/fem/node.h> 00025 #include <mechsys/fem/geomelem.h> 00026 00027 namespace FEM 00028 { 00029 00030 class Tri15: public GeomElem 00031 { 00032 public: 00033 // Auxiliar structure to map local face IDs to local node IDs 00034 static const int Face2Node[3][5]; // 3 edges, 5 nodes/edge 00035 00036 // Constructor 00037 Tri15 (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 00053 y 2 00054 | @ @ 00055 | / \ / \ 00056 +--x / \ / \ 00057 10 @ @ 9 / \ 00058 / \ / \ 00059 / 14 \ / \ 00060 5 @ @ @ 4 2 / \ 1 00061 / \ / \ 00062 / \ / \ 00063 11 @ @ @ @ 8 / \ 00064 / 12 13 \ / \ 00065 / \ / \ 00066 @-----@-----@-----@-----@ @-----------------------@ 00067 0 6 3 7 1 0 00068 */ 00069 const int Tri15::Face2Node[3][5]= {{ 0, 1, 3, 6, 7 }, 00070 { 1, 2, 4, 8, 9 }, 00071 { 2, 0, 5, 10, 11 }}; // order of nodes is important 00072 00073 00075 00076 00077 inline Tri15::Tri15 (int NDim) 00078 : GeomElem(NDim, /*NN*/15, /*NFN*/5, "Tri15") 00079 { 00080 SetIPs (16); 00081 } 00082 00083 inline void Tri15::SetIPs (int TotNIP) 00084 { 00085 if (TotNIP==12) IPs = TRI_IP12; 00086 else if (TotNIP==13) IPs = TRI_IP12; 00087 else if (TotNIP==16) IPs = TRI_IP16; 00088 else if (TotNIP==25) IPs = TRI_IP25; 00089 else throw new Fatal("tri6::SetIntPoints: Total number of integration points = %d is invalid",TotNIP); 00090 00091 NIP = TotNIP; 00092 FIPs = LIN_IP5; 00093 NFIP = 5; 00094 } 00095 00096 inline void Tri15::Shape (double r, double s, double t) const 00097 { 00098 /* s 00099 * ^ 00100 * | 00101 * 2 00102 * @,(0,1) 00103 * | ', 00104 * | ', 9 00105 * 10 @ @, 00106 * | 14 ', 4 00107 * 5 @ @ @ 00108 * | ', 8 00109 * 11 @ 12@ @ '@ 00110 * | 13 ', 00111 * |(0,0) ', (1,0) 00112 * @----@----@----@----@ --> r 00113 * 0 6 3 7 1 00114 */ 00115 double pt1 = 42.666666666666667; 00116 double pt2 = 10.666666666666667; 00117 double cc = 1.0 - r - s; 00118 double t1 = r - 0.25; 00119 double t2 = r - 0.5; 00120 double t3 = r - 0.75; 00121 double t4 = s - 0.25; 00122 double t5 = s - 0.5; 00123 double t6 = s - 0.75; 00124 double t7 = cc - 0.25; 00125 double t8 = cc - 0.5; 00126 double t9 = cc - 0.75; 00127 N(0) = pt2*cc*t7*t8*t9; 00128 N(1) = pt2*r*t1*t2*t3; 00129 N(2) = pt2*s*t4*t5*t6; 00130 N(3) = 64.0*cc*r*t1*t7; 00131 N(4) = 64.0*r*s*t1*t4; 00132 N(5) = 64.0*s*cc*t4*t7; 00133 N(6) = pt1*cc*r*t7*t8; 00134 N(7) = pt1*cc*r*t1*t2; 00135 N(8) = pt1*r*s*t1*t2; 00136 N(9) = pt1*r*s*t4*t5; 00137 N(10) = pt1*s*cc*t4*t5; 00138 N(11) = pt1*s*cc*t7*t8; 00139 N(12) = 128.0*r*s*cc*t7; 00140 N(13) = 128.0*r*s*t1*cc; 00141 N(14) = 128.0*r*s*cc*t4; 00142 } 00143 00144 inline void Tri15::Derivs (double r, double s, double t) const 00145 { 00146 /* _ _ T 00147 * | dNi | 00148 * dNdR = | --- | , where Rj = r, s 00149 * |_ dRj _| 00150 * 00151 * dNdR(j,i), j=>local coordinate and i=>shape function 00152 */ 00153 double pt1 = 42.666666666666667; 00154 double pt2 = 10.666666666666667; 00155 double cc = 1.0 - r - s; 00156 double t1 = r - 0.25; 00157 double t2 = r - 0.5; 00158 double t3 = r - 0.75; 00159 double t4 = s - 0.25; 00160 double t5 = s - 0.5; 00161 double t6 = s - 0.75; 00162 double t7 = cc - 0.25; 00163 double t8 = cc - 0.5; 00164 double t9 = cc - 0.75; 00165 00166 dNdR(0, 0) = - pt2 * (t8 * t9 * (t7 + cc) + cc * t7 * (t8 + t9) ); 00167 dNdR(0, 1) = pt2 * (t2 * t3 * (t1 + r) + r * t1 * (t3 + t2) ); 00168 dNdR(0, 2) = 0.0; 00169 dNdR(0, 3) = 64.0 * (cc * t7 * (t1 + r) - r * t1 * (t7 + cc) ); 00170 dNdR(0, 4) = 64.0 * s * t4 * (t1 + r); 00171 dNdR(0, 5) = - 64.0 * s * t4 * (t7 + cc); 00172 dNdR(0, 6) = pt1 * (cc * t7 * t8 - r * (t8 * (t7 + cc) + cc * t7) ); 00173 dNdR(0, 7) = pt1 * (cc * (t2 * (t1 + r) + r * t1) - r * t1 * t2); 00174 dNdR(0, 8) = pt1 * s * (t2 * (t1 + r) + r * t1); 00175 dNdR(0, 9) = pt1 * s * t4 * t5; 00176 dNdR(0, 10) = - pt1 * s * t4 * t5; 00177 dNdR(0, 11) = - pt1 * s * (t8 * (t7 + cc) + cc * t7); 00178 dNdR(0, 12) = 128.0 * s * (cc * t7 - r * (t7 + cc) ); 00179 dNdR(0, 13) = 128.0 * s * (cc * (t1 + r) - r * t1); 00180 dNdR(0, 14) = 128.0 * s * t4 * (cc - r); 00181 00182 dNdR(1, 0) = - pt2 * (t8 * t9 * (t7 + cc) + cc * t7 * (t8 + t9) ); 00183 dNdR(1, 1) = 0.0; 00184 dNdR(1, 2) = pt2 * (t5 * t6 * (t4 + s) + s * t4 * (t6 + t5) ); 00185 dNdR(1, 3) = - 64.0 * r * t1 * (t7 + cc); 00186 dNdR(1, 4) = 64.0 * r * t1 * (t4 + s); 00187 dNdR(1, 5) = 64.0 * (cc * t7 * (t4 + s) - s * t4 * (t7 + cc) ); 00188 dNdR(1, 6) = - pt1 * r * (t8 * (t7 + cc) + cc * t7); 00189 dNdR(1, 7) = - pt1 * r * t1 * t2; 00190 dNdR(1, 8) = pt1 * r * t1 * t2; 00191 dNdR(1, 9) = pt1 * r * (t5 * (t4 + s) + s * t4); 00192 dNdR(1, 10) = pt1 * ((cc * (t5 * (t4 + s) + s * t4)) - s * t4 * t5); 00193 dNdR(1, 11) = pt1 * (cc * t7 * t8 - s * (t8 * (t7 + cc) + cc * t7)); 00194 dNdR(1, 12) = 128.0 * r * (cc * t7 - s * (cc + t7) ); 00195 dNdR(1, 13) = 128.0 * r * t1 * (cc - s); 00196 dNdR(1, 14) = 128.0 * r * (cc * (t4 + s) - s * t4); 00197 } 00198 00199 inline void Tri15::FaceShape (double r, double s) const 00200 { 00201 /* 00202 * @-----@-----@-----@-----@-> r 00203 * 0 3 2 4 1 00204 * | | | 00205 * r=-1 -1/2 r=0 1/2 r=+1 00206 */ 00207 FN(0) = (r-1.)*(1.-2.*r)*r*(-1.-2.*r)/6.; 00208 FN(1) = (1.-2.*r)*r*(-1.-2.*r)*(1.+r)/6.; 00209 FN(2) = (1.-r)*(1.-2.*r)*(-1.-2.*r)*(-1.-r); 00210 FN(3) = 4.*(1.-r)*(1.-2.*r)*r*(-1.-r)/3.; 00211 FN(4) = 4.*(1.-r)*r*(-1.-2.*r)*(-1.-r)/3.; 00212 } 00213 00214 inline void Tri15::FaceDerivs (double r, double s) const 00215 { 00216 /* _ _ T 00217 * | dNi | 00218 * FdNdR = | --- | , where Rj = r, s 00219 * |_ dRj _| 00220 * 00221 * FdNdR(j,i), j=>local coordinate and i=>shape function 00222 */ 00223 FdNdR(0,0) = -((1.-2.*r)*(r-1.)*r)/3.-((-2.*r-1.)*(r-1.)*r)/3.+((-2.*r-1.)*(1.-2.*r)*r)/6.+((-2.*r-1.)*(1.-2.*r)*(r-1.))/6.; 00224 FdNdR(0,1) = -((1.-2.*r)*r*(r+1.))/3.-((-2.*r-1.)*r*(r+1.))/3.+((-2.*r-1.)*(1.-2.*r)*(r+1.))/6.+((-2.*r-1.)*(1.-2.*r)*r)/6.; 00225 FdNdR(0,2) = -2.*(1.-2.*r)*(-r-1.)*(1.-r)-2.*(-2.*r-1.)*(-r-1.)*(1.-r)-(-2.*r-1.)*(1.-2.*r)*(1.-r)-(-2.*r-1.)*(1.-2.*r)*(-r-1.); 00226 FdNdR(0,3) = -(8.*(-r-1.)*(1.-r)*r)/3.-(4.*(1.-2.*r)*(1.-r)*r)/3.-(4.*(1.-2.*r)*(-r-1.)*r)/3.+(4.*(1.-2.*r)*(-r-1.)*(1.-r))/3.; 00227 FdNdR(0,4) = -(8.*(-r-1.)*(1.-r)*r)/3.-(4.*(-2.*r-1.)*(1.-r)*r)/3.-(4.*(-2.*r-1.)*(-r-1.)*r)/3.+(4.*(-2.*r-1.)*(-r-1.)*(1.-r))/3.; 00228 } 00229 00230 inline void Tri15::NatCoords (Mat_t & C) const 00231 { 00232 throw new Fatal("Tri15::NatCoords: method is not available yet"); 00233 } 00234 00235 00237 00238 00239 // Allocate a new element 00240 GeomElem * Tri15Maker (int NDim) { return new Tri15(NDim); } 00241 00242 // Register element 00243 int Tri15Register () 00244 { 00245 GeomElemFactory["Tri15"] = Tri15Maker; 00246 GEOM.Set ("Tri15", (double)GEOM.Keys.Size()); 00247 return 0; 00248 } 00249 00250 // Call register 00251 int __Tri15_dummy_int = Tri15Register(); 00252 00253 }; // namespace FEM 00254 00255 #endif // MECHSYS_FEM_TRI15_H
1.7.6.1
