 |
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_ANISOINVS_H 00021 #define MECHSYS_ANISOINVS_H 00022 00023 // Std Lib 00024 #include <iostream> 00025 #include <fstream> 00026 00027 // MechSys 00028 #include <mechsys/util/maps.h> 00029 #include <mechsys/util/fatal.h> 00030 #include <mechsys/linalg/matvec.h> 00031 #include <mechsys/linalg/jacobirot.h> 00032 00033 class AnisoInvs 00034 { 00035 public: 00036 // Constructor & Destructor 00037 AnisoInvs (double b, double Alpha, Vec3_t const & a, bool Obliq=true, bool Check=true, bool UseJacobi=false); 00038 ~AnisoInvs () {} 00039 00040 // Methods 00041 void Calc (Vec_t const & Sig, bool WithDerivs=false); 00042 void CheckSym (char const * Name, Ten2_t const & T) const; 00043 void CheckSym (char const * Name, Ten3_t const & T) const; 00044 00045 // Constants 00046 double b, Alpha; 00047 Vec3_t a; 00048 bool Obliq; 00049 bool Check; 00050 bool UseJacobi; 00051 double Tol; 00052 double Zero; 00053 Vec3_t au; 00054 00055 // Data 00056 static double sp, sq; // invariants 00057 static Vec3_t N123, N, Nu; // SMP 00058 static Vec3_t n, nu; // AMP 00059 static Mat3_t Q; // rotation matrix 00060 static Vec3_t L; // eigenvalues 00061 static Vec3_t v0,v1,v2; // eigenvectors 00062 static Vec3_t t, p, q; // traction and projections 00063 static Mat3_t P; // projector 00064 00065 // Derivs 00066 #ifdef HAS_TENSORS 00067 static Ten1_t V0, V1, V2; // eigenvectors 00068 static Ten2_t E0, E1, E2; // eigenprojectors 00069 static Ten3_t dv0dSig, dv1dSig, dv2dSig; // eigenvectors 00070 static Ten2_t dN123dL; // SMP(123) 00071 static Ten3_t dNdSig, dNudSig; // SMP 00072 static Ten2_t dNudN; // SMP 00073 static Ten3_t dndSig, dnudSig; // AMP 00074 static Ten2_t dnudn; // AMP 00075 static Ten2_t tSig, I; // tensor Sig and Identity 00076 static Ten1_t tN, tNu, tn, tnu; // tensors: normal vectors 00077 static Ten3_t M, dtdSig; // auxiliary tensor and traction 00078 static Ten1_t tt, tp, tq; // tensor: traction and projections 00079 static Ten2_t tP; // tensor: projector 00080 static Ten2_t dsdSig; // auxiliary tensor for dpdSig 00081 static Ten3_t dpdSig, dqdSig; // projections 00082 static Ten2_t tdspdSig, tdsqdSig; // invariants (tensor) 00083 static Vec_t dspdSig, dsqdSig; // invariants (Mandel's basis) 00084 #endif 00085 00086 private: 00087 bool _allocated; 00088 }; 00089 00090 // Data 00091 double AnisoInvs::sp; double AnisoInvs::sq; 00092 Vec3_t AnisoInvs::N123; Vec3_t AnisoInvs::N; Vec3_t AnisoInvs::Nu; 00093 Vec3_t AnisoInvs::n; Vec3_t AnisoInvs::nu; 00094 Mat3_t AnisoInvs::Q; 00095 Vec3_t AnisoInvs::L; 00096 Vec3_t AnisoInvs::v0; Vec3_t AnisoInvs::v1; Vec3_t AnisoInvs::v2; 00097 Vec3_t AnisoInvs::t; Vec3_t AnisoInvs::p; Vec3_t AnisoInvs::q; 00098 Mat3_t AnisoInvs::P; 00099 00100 // Derivs 00101 #ifdef HAS_TENSORS 00102 Ten1_t AnisoInvs::V0; Ten1_t AnisoInvs::V1; Ten1_t AnisoInvs::V2; 00103 Ten3_t AnisoInvs::dv0dSig; Ten3_t AnisoInvs::dv1dSig; Ten3_t AnisoInvs::dv2dSig; 00104 Ten2_t AnisoInvs::E0; Ten2_t AnisoInvs::E1; Ten2_t AnisoInvs::E2; 00105 Ten2_t AnisoInvs::dN123dL; 00106 Ten3_t AnisoInvs::dNdSig; Ten3_t AnisoInvs::dNudSig; 00107 Ten2_t AnisoInvs::dNudN; 00108 Ten3_t AnisoInvs::dndSig; Ten3_t AnisoInvs::dnudSig; 00109 Ten2_t AnisoInvs::dnudn; 00110 Ten2_t AnisoInvs::tSig; Ten2_t AnisoInvs::I; 00111 Ten1_t AnisoInvs::tN; Ten1_t AnisoInvs::tNu; Ten1_t AnisoInvs::tn; Ten1_t AnisoInvs::tnu; 00112 Ten3_t AnisoInvs::M; Ten3_t AnisoInvs::dtdSig; 00113 Ten1_t AnisoInvs::tt; Ten1_t AnisoInvs::tp; Ten1_t AnisoInvs::tq; 00114 Ten2_t AnisoInvs::tP; 00115 Ten2_t AnisoInvs::dsdSig; 00116 Ten3_t AnisoInvs::dpdSig; Ten3_t AnisoInvs::dqdSig; 00117 Ten2_t AnisoInvs::tdspdSig; Ten2_t AnisoInvs::tdsqdSig; 00118 Vec_t AnisoInvs::dspdSig; Vec_t AnisoInvs::dsqdSig; 00119 #endif 00120 00121 00123 00124 00125 inline AnisoInvs::AnisoInvs (double Theb, double TheAlpha, Vec3_t const & Thea, bool TheObliq, bool TheCheck, bool TheUseJacobi) 00126 : b(Theb), Alpha(TheAlpha), a(Thea), Obliq(TheObliq), Check(TheCheck), UseJacobi(TheUseJacobi), Tol(1.0e-8), Zero(1.0e-14), _allocated(false) 00127 { 00128 if (!_allocated) 00129 { 00130 I.SetDiagonal (1.0); 00131 _allocated = true; 00132 } 00133 00134 // check 00135 double norm_a = Norm(a); 00136 if (norm_a<Tol) throw new Fatal("AnisoInvs::AnisoInvs 'a' vector must be non zero. a=[%g,%g,%g]. norm(a)=%g",a(0),a(1),a(2),norm_a); 00137 au = a / norm_a; 00138 } 00139 00140 inline void AnisoInvs::Calc (Vec_t const & Sig, bool WithDerivs) 00141 { 00142 // principal values and eigenvectors 00143 if (UseJacobi) 00144 { 00145 // rotation matrix 00146 JacobiRot (Sig, Q, L); 00147 00148 // eigenvectors 00149 v0(0) = Q(0,0); v1(0) = Q(0,1); v2(0) = Q(0,2); 00150 v0(1) = Q(1,0); v1(1) = Q(1,1); v2(1) = Q(1,2); 00151 v0(2) = Q(2,0); v1(2) = Q(2,1); v2(2) = Q(2,2); 00152 } 00153 else 00154 { 00155 // eigenvectors 00156 Eig (Sig, L, v0, v1, v2); 00157 00158 // rotation matrix 00159 Q(0,0) = v0(0); Q(0,1) = v1(0); Q(0,2) = v2(0); 00160 Q(1,0) = v0(1); Q(1,1) = v1(1); Q(1,2) = v2(1); 00161 Q(2,0) = v0(2); Q(2,1) = v1(2); Q(2,2) = v2(2); 00162 } 00163 //if (L(0)>0.0 || L(1)>0.0 || L(2)>0.0) throw new Fatal("AnisoInvs::Calc: This method only works when all principal values are negative (compression octant). L=[%g,%g,%g]",L(0),L(1),L(2)); 00164 00165 // normal vectors 00166 bool zero = (fabs(L(0))<Tol || fabs(L(1))<Tol || fabs(L(2))<Tol); 00167 if (zero) N123 = 1.0, 1.0, 1.0; // SMP(123) 00168 else N123 = 1.0/pow(fabs(L(0)),b), 1.0/pow(fabs(L(1)),b), 1.0/pow(fabs(L(2)),b); // SMP(123) 00169 00170 Mat3_t Qt; 00171 Qt(0,0) = Q(0,0); Qt(0,1) = Q(1,0); Qt(0,2) = Q(2,0); 00172 Qt(1,0) = Q(0,1); Qt(1,1) = Q(1,1); Qt(1,2) = Q(2,1); 00173 Qt(2,0) = Q(0,2); Qt(2,1) = Q(1,2); Qt(2,2) = Q(2,2); 00174 00175 N = product (Q, N123); // SMP(xyz) 00176 //N = product (Qt, N123); // SMP(xyz) 00177 Nu = N / Norm(N); // SMP 00178 n = Alpha*au + Nu; // AMP 00179 nu = n / Norm(n); // AMP 00180 00181 // traction 00182 Mult (Sig, nu, t); // t = Sig * nu 00183 00184 // P projector 00185 double s = 0.0; 00186 if (Obliq) 00187 { 00188 s = 1.0 / dot(N, n); 00189 Dyad (s, N, n, P); // P = s*(N dy n) 00190 } 00191 else Dyad (nu, nu, P); 00192 00193 // invariants 00194 p = product (P, t); 00195 q = t - p; 00196 sp = Norm(p); 00197 sq = Norm(q); 00198 00199 // derivatives 00200 if (WithDerivs) 00201 { 00202 #ifdef HAS_TENSORS 00203 // normal to SMP in principal system 00204 if (zero) dN123dL.SetDiagonal (1.0); 00205 else dN123dL = -b/(L(0)*pow(fabs(L(0)),b)), 0.0, 0.0, 00206 0.0, -b/(L(1)*pow(fabs(L(1)),b)), 0.0, 00207 0.0, 0.0, -b/(L(2)*pow(fabs(L(2)),b)); 00208 00209 // eigenvectors and eigenprojectors 00210 Vec2Tensor (v0,V0); 00211 Vec2Tensor (v1,V1); 00212 Vec2Tensor (v2,V2); 00213 EigenVecDerivs (L, V0,V1,V2, dv0dSig,dv1dSig,dv2dSig); 00214 E0 = (V0 & V0); 00215 E1 = (V1 & V1); 00216 E2 = (V2 & V2); 00217 00218 // normal to SMP in xyz system 00219 dNdSig = ((V0&E0)*dN123dL[0][0]) + ((V1&E1)*dN123dL[1][1]) + ((V2&E2)*dN123dL[2][2]) + (dv0dSig*N123(0)) + (dv1dSig*N123(1)) + (dv2dSig*N123(2)); 00220 00221 // unit normal to SMP in xyz system 00222 UnitVecDeriv (N, Nu, dNudN, Tol); 00223 dNudSig = dNudN * dNdSig; 00224 00225 // normal to AMP in xyz system 00226 dndSig = dNudSig; 00227 00228 // unit normal to AMP in xyz system 00229 UnitVecDeriv (n, nu, dnudn, Tol); 00230 dnudSig = dnudn * dndSig; 00231 00232 // convert to tensors 00233 Mat2Tensor (P, tP); 00234 Ten2Tensor (Sig, tSig); 00235 Vec2Tensor (t, tt); 00236 Vec2Tensor (N, tN); 00237 Vec2Tensor (Nu, tNu); 00238 Vec2Tensor (n, tn); 00239 Vec2Tensor (nu, tnu); 00240 Vec2Tensor (p, tp); 00241 Vec2Tensor (q, tq); 00242 00243 // traction 00244 for (size_t i=0; i<3; ++i) 00245 for (size_t j=0; j<3; ++j) 00246 for (size_t k=0; k<3; ++k) 00247 M[i][j][k] = 0.5*(I[i][j]*nu[k] + I[i][k]*nu[j]); 00248 dtdSig = M + (tSig * dnudSig); 00249 00250 // normal projection 00251 if (Obliq) 00252 { 00253 double m = tn * tt; 00254 dsdSig = (tn*dNdSig + tN*dndSig)*(-s*s); 00255 dpdSig = (tN&dsdSig)*m + dNdSig*(s*m) + (tN&(tt*dndSig))*s + tP*dtdSig; 00256 } 00257 else 00258 { 00259 double m = tnu * tt; 00260 dpdSig = dnudSig*m + (tnu&(tt*dnudSig)) + tP*dtdSig; 00261 } 00262 00263 // on-plane projection 00264 dqdSig = dtdSig - dpdSig; 00265 00266 // check symmetry 00267 if (Check) 00268 { 00269 CheckSym ("dNdSig", dNdSig ); 00270 CheckSym ("dNudSig", dNudSig); 00271 CheckSym ("dndSig", dndSig ); 00272 CheckSym ("dnudSig", dnudSig); 00273 CheckSym ("dsdSig", dsdSig ); 00274 CheckSym ("dpdSig", dpdSig ); 00275 CheckSym ("M", M ); 00276 CheckSym ("dtdSig", dtdSig ); 00277 CheckSym ("dqdSig", dqdSig ); 00278 } 00279 00280 // invariants 00281 if (sp>Tol) tdspdSig = (tp*(1.0/sp)) * dpdSig; else tdspdSig.SetDiagonal (1.0); 00282 if (sq>Tol) tdsqdSig = (tq*(1.0/sq)) * dqdSig; else tdsqdSig.SetDiagonal (1.0); 00283 size_t ncp = size(Sig); 00284 Tensor2Ten (tdspdSig, dspdSig, ncp); 00285 Tensor2Ten (tdsqdSig, dsqdSig, ncp); 00286 #else 00287 throw new Fatal("AnisoInvs::Calc: Tensors library is required in order to calculate derivatives"); 00288 #endif 00289 } 00290 } 00291 00292 inline void AnisoInvs::CheckSym (char const * Name, Ten2_t const & T) const 00293 { 00294 for (size_t i=0; i<3; ++i) 00295 for (size_t j=0; j<3; ++j) 00296 { 00297 double error = fabs(T[i][j] - T[i][j]); 00298 if (error>Zero) throw new Fatal("CheckSym: Second order tensor %s is not symmetric. error=%g",Name,error); 00299 } 00300 } 00301 00302 inline void AnisoInvs::CheckSym (char const * Name, Ten3_t const & T) const 00303 { 00304 for (size_t i=0; i<3; ++i) 00305 for (size_t j=0; j<3; ++j) 00306 for (size_t k=0; k<3; ++k) 00307 { 00308 double error = fabs(T[i][j][k] - T[i][k][j]); 00309 if (error>Zero) throw new Fatal("CheckSym: Third order tensor %s is not minor symmetric. error=%g",Name,error); 00310 } 00311 } 00312 00313 #endif // MECHSYS_ANISOINVS_H
1.7.6.1
