 |
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 * 00004 * * 00005 * This program is free software: you can redistribute it and/or modify * 00006 * it under the terms of the GNU General Public License as published by * 00007 * the Free Software Foundation, either version 3 of the License, or * 00008 * any later version. * 00009 * * 00010 * This program is distributed in the hope that it will be useful, * 00011 * but WITHOUT ANY WARRANTY; without even the implied warranty of * 00012 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * 00013 * GNU General Public License for more details. * 00014 * * 00015 * You should have received a copy of the GNU General Public License * 00016 * along with this program. If not, see <http://www.gnu.org/licenses/> * 00017 ************************************************************************/ 00018 00019 /* Tensors - Copyright (C) 2007 Dorival de Moraes Pedroso */ 00020 00021 #ifndef MPM_TENSORS_H 00022 #define MPM_TENSORS_H 00023 00024 // STL 00025 #include <cmath> // for sqrt ... 00026 #include <cfloat> // for DBL_EPSILON 00027 #include <sstream> 00028 00029 // MechSys 00030 #include <mechsys/linalg/matvec.h> 00031 00032 // Local 00033 #include <mechsys/mpm/defs.h> 00034 #include <mechsys/mpm/jacobirot.h> 00035 #include <mechsys/mpm/fmtnum.h> 00036 00037 namespace MPM { 00038 00039 /* Second order tensors: 00040 * y Sy 00041 * | | 00042 * |__x ____|____________ 00043 * ,' ,' | ,'| 00044 * z ,' \|/ ,' | 00045 * ,' <-----'--- ,' | 00046 * ,' Sxy ,' | 00047 * ,'_______________,' | | 00048 * | | | | 00049 * | | | <------- Sx 00050 * | __, | \| | 00051 * | ,'| | ' ,' 00052 * | ,' | Sxy ,' 00053 * | ,' | ,' 00054 * | Sz | ,' 00055 * |________________|' 00056 * 00057 * (Symmetric) Tensor components (Mandel's basis) 00058 * / Sx Sxy Sxz \ 00059 * | Sxy Sy Syz | => { Sx, Sy, Sz, SQ2*Sxy, SQ2*Syz, SQ2*Sxz } 00060 * \ Sxz Syz Sz / 00061 * 00062 * (Assymmetric) Tensor components 00063 * / Sx Sxy Sxz \ 00064 * | Syx Sy Syz | => { Sx, Sy, Sz, Sxy, Syz, Sxz, Syx, Szy, Szx } 00065 * \ Szx Szy Sz / 00066 */ 00067 00068 00070 00072 STensor2 SymI; 00073 00075 STensor4 SymII; 00076 00078 STensor4 SymIdyI; 00079 00081 STensor4 SymPsd; 00082 00084 STensor4 SymPiso; 00085 00086 00088 00091 inline double NormV (Vector3D const & v) 00092 { 00093 return sqrt(v(0)*v(0) + v(1)*v(1) + v(2)*v(2)); 00094 } 00095 00096 inline double NormT (STensor2 const & T) 00097 { 00098 return sqrt(T(0)*T(0) + T(1)*T(1) + T(2)*T(2) + T(3)*T(3) + T(4)*T(4) + T(5)*T(5)); 00099 } 00100 00102 inline void AddScaled (double const & a, STensor4 const & X, 00103 double const & b, STensor4 const & Y, STensor4 & Z) 00104 { 00105 Z(0,0)=a*X(0,0)+b*Y(0,0); Z(0,1)=a*X(0,1)+b*Y(0,1); Z(0,2)=a*X(0,2)+b*Y(0,2); Z(0,3)=a*X(0,3)+b*Y(0,3); Z(0,4)=a*X(0,4)+b*Y(0,4); Z(0,5)=a*X(0,5)+b*Y(0,5); 00106 Z(1,0)=a*X(1,0)+b*Y(1,0); Z(1,1)=a*X(1,1)+b*Y(1,1); Z(1,2)=a*X(1,2)+b*Y(1,2); Z(1,3)=a*X(1,3)+b*Y(1,3); Z(1,4)=a*X(1,4)+b*Y(1,4); Z(1,5)=a*X(1,5)+b*Y(1,5); 00107 Z(2,0)=a*X(2,0)+b*Y(2,0); Z(2,1)=a*X(2,1)+b*Y(2,1); Z(2,2)=a*X(2,2)+b*Y(2,2); Z(2,3)=a*X(2,3)+b*Y(2,3); Z(2,4)=a*X(2,4)+b*Y(2,4); Z(2,5)=a*X(2,5)+b*Y(2,5); 00108 Z(3,0)=a*X(3,0)+b*Y(3,0); Z(3,1)=a*X(3,1)+b*Y(3,1); Z(3,2)=a*X(3,2)+b*Y(3,2); Z(3,3)=a*X(3,3)+b*Y(3,3); Z(3,4)=a*X(3,4)+b*Y(3,4); Z(3,5)=a*X(3,5)+b*Y(3,5); 00109 Z(4,0)=a*X(4,0)+b*Y(4,0); Z(4,1)=a*X(4,1)+b*Y(4,1); Z(4,2)=a*X(4,2)+b*Y(4,2); Z(4,3)=a*X(4,3)+b*Y(4,3); Z(4,4)=a*X(4,4)+b*Y(4,4); Z(4,5)=a*X(4,5)+b*Y(4,5); 00110 Z(5,0)=a*X(5,0)+b*Y(5,0); Z(5,1)=a*X(5,1)+b*Y(5,1); Z(5,2)=a*X(5,2)+b*Y(5,2); Z(5,3)=a*X(5,3)+b*Y(5,3); Z(5,4)=a*X(5,4)+b*Y(5,4); Z(5,5)=a*X(5,5)+b*Y(5,5); 00111 } 00112 00113 00115 00117 inline void Sym (double a, ATensor2 const & T, STensor2 & S) 00118 { 00119 S(0) = a*T(0); 00120 S(1) = a*T(1); 00121 S(2) = a*T(2); 00122 S(3) = a*0.5*(T(3)+T(6))*SQ2; 00123 S(4) = a*0.5*(T(4)+T(7))*SQ2; 00124 S(5) = a*0.5*(T(5)+T(8))*SQ2; 00125 } 00126 00128 inline double Det (ATensor2 const & B) 00129 { 00130 return B(0)*(B(1)*B(2)-B(4)*B(7))-B(3)*(B(6)*B(2)-B(4)*B(8))+B(5)*(B(6)*B(7)-B(1)*B(8)); 00131 } 00132 00134 inline double Det (STensor2 const & T) 00135 { 00136 return T(0)*T(1)*T(2) + T(3)*T(4)*T(5)/SQ2 - T(0)*T(4)*T(4)/2.0 - T(1)*T(5)*T(5)/2.0 - T(2)*T(3)*T(3)/2.0; 00137 } 00138 00140 inline double Cam_p (STensor2 const & S) 00141 { 00142 return (S(0)+S(1)+S(2))/3.0; 00143 } 00144 00146 inline double Cam_q (STensor2 const & S) 00147 { 00148 return sqrt(( (S(0)-S(1))*(S(0)-S(1)) 00149 + (S(1)-S(2))*(S(1)-S(2)) 00150 + (S(2)-S(0))*(S(2)-S(0)) + 3.0*(S(3)*S(3)+S(4)*S(4)+S(5)*S(5)) )/2.0); 00151 } 00152 00153 inline void Stress_p_q (STensor2 const & Sig, double & p, double & q) 00154 { 00155 p = (Sig(0)+Sig(1)+Sig(2))/3.0; 00156 q = sqrt(( (Sig(0)-Sig(1))*(Sig(0)-Sig(1)) 00157 + (Sig(1)-Sig(2))*(Sig(1)-Sig(2)) 00158 + (Sig(2)-Sig(0))*(Sig(2)-Sig(0)) + 3.0*(Sig(3)*Sig(3) + Sig(4)*Sig(4) + Sig(5)*Sig(5)) )/2.0); 00159 } 00160 00161 inline void Eigenvals (STensor2 const & T, double L[3]) 00162 { 00163 // Calculate eigenvalues and eigenvectors 00164 if (JacobiRot(T,L)==-1) 00165 throw new Fatal("Eigenvals: Jacobi rotation did not converge. T=<%.6f %.6f %.6f %.6f %.6f %.6f>",T(0),T(1),T(2),T(3)/SQ2,T(4)/SQ2,T(5)/SQ2); 00166 } 00167 00168 inline void Eigenvp (STensor2 const & T, double L[3], STensor2 P[3]) 00169 { 00170 double V0[3]; // eigenvector 00171 double V1[3]; // eigenvector 00172 double V2[3]; // eigenvector 00173 00174 // Calculate eigenvalues and eigenvectors 00175 if (JacobiRot(T,V0,V1,V2,L)==-1) 00176 throw new Fatal("Eigenvp: Jacobi rotation did not converge. T=<%.6f %.6f %.6f %.6f %.6f %.6f>",T(0),T(1),T(2),T(3)/SQ2,T(4)/SQ2,T(5)/SQ2); 00177 00178 // Calculate eigenprojectors (Pi = Vi dyad Vi) 00179 // Eigenprojector 0 00180 P[0](0) = V0[0] * V0[0]; 00181 P[0](1) = V0[1] * V0[1]; 00182 P[0](2) = V0[2] * V0[2]; 00183 P[0](3) = V0[0] * V0[1] * SQ2; 00184 P[0](4) = V0[1] * V0[2] * SQ2; 00185 P[0](5) = V0[0] * V0[2] * SQ2; 00186 // Eigenprojector 1 00187 P[1](0) = V1[0] * V1[0]; 00188 P[1](1) = V1[1] * V1[1]; 00189 P[1](2) = V1[2] * V1[2]; 00190 P[1](3) = V1[0] * V1[1] * SQ2; 00191 P[1](4) = V1[1] * V1[2] * SQ2; 00192 P[1](5) = V1[0] * V1[2] * SQ2; 00193 // Eigenprojector 2 00194 P[2](0) = V2[0] * V2[0]; 00195 P[2](1) = V2[1] * V2[1]; 00196 P[2](2) = V2[2] * V2[2]; 00197 P[2](3) = V2[0] * V2[1] * SQ2; 00198 P[2](4) = V2[1] * V2[2] * SQ2; 00199 P[2](5) = V2[0] * V2[2] * SQ2; 00200 } 00201 00202 inline void Stress_psmp_qsmp (STensor2 const & Sig, double & psmp, double & qsmp) 00203 { 00204 // Principal values 00205 double L[3]; 00206 Eigenvals (Sig, L); 00207 00208 // Invariants 00209 if (L[0]>0.0 && L[1]>0.0 && L[2]>0.0) 00210 { 00211 // Characteristic invariants 00212 double Is[4]; 00213 Is[0] = L[0]+L[1]+L[2]; 00214 Is[1] = L[0]*L[1]+L[1]*L[2]+L[2]*L[0]; 00215 Is[2] = L[0]*L[1]*L[2]; 00216 Is[3] = sqrt(L[0])+sqrt(L[1])+sqrt(L[2]); 00217 00218 // psmp and qsmp invariants 00219 if (Is[1]>0.0) 00220 { 00221 psmp = Is[3]*sqrt(Is[2]/Is[1]); 00222 double del = L[0]*L[0]+L[1]*L[1]+L[2]*L[2]-psmp*psmp; 00223 if (fabs(del)<DBL_EPSILON*100.0) del = 0.0; 00224 if (del>=0.0) qsmp = sqrt(del); 00225 else qsmp=0.0; //throw new Fatal(_("Stress_psmp_qsmp: del=s1^2+s2^2+s3^2-psmp^2=%f must be positive. s1=%f, s2=%f, s3=%f, psmp=%f"),del,L[0],L[1],L[2],psmp); 00226 } 00227 else throw new Fatal("Stress_psmp_qsmp: I2=%f must be positive",Is[1]); 00228 } 00229 else throw new Fatal("Stress_psmp_qsmp: The principal stresses must be positive. s1=%f, s2=%f, s3=%f",L[0],L[1],L[2]); 00230 } 00231 00232 inline void Derivs_psmp_qsmp (STensor2 const & Sig, double & psmp, double & qsmp, double DpsmpDs[3], double DqsmpDs[3], STensor2 Proj[3]) 00233 { 00234 // Eigenvalues and Eigenprojectors 00235 double L[3]; 00236 Eigenvp (Sig, L, Proj); 00237 00238 // Invariants and derivatives 00239 if (L[0]>0.0 && L[1]>0.0 && L[2]>0.0) 00240 { 00241 // Characteristic invariants 00242 double Is[4]; 00243 Is[0] = L[0]+L[1]+L[2]; 00244 Is[1] = L[0]*L[1]+L[1]*L[2]+L[2]*L[0]; 00245 Is[2] = L[0]*L[1]*L[2]; 00246 Is[3] = sqrt(L[0])+sqrt(L[1])+sqrt(L[2]); 00247 00248 // Derivatives and invariants 00249 if (Is[1]>0.0) 00250 { 00251 // Invariants 00252 psmp = Is[3]*sqrt(Is[2]/Is[1]); 00253 double del = L[0]*L[0]+L[1]*L[1]+L[2]*L[2]-psmp*psmp; 00254 if (fabs(del)<DBL_EPSILON*100.0) del = 0.0; 00255 if (del>=0.0) qsmp = sqrt(del); 00256 else throw new Fatal("Derivs_psmp_qsmp: del=s1^2+s2^2+s3^2-psmp^2=%f must be positive. s1=%f, s2=%f, s3=%f, psmp=%f",del,L[0],L[1],L[2],psmp); 00257 00258 // Derivatives 00259 double dI1ds[3], dI2ds[3], dI3ds[3]; 00260 for (int k=0; k<3; ++k) 00261 { 00262 // Derivatives of the characteristic invariants 00263 dI1ds[k] = Is[0]-L[k]; 00264 dI2ds[k] = Is[2]/L[k]; 00265 dI3ds[k] = 0.5/sqrt(L[k]); 00266 00267 // Derivatives of the psmp and qsmp invariants 00268 DpsmpDs[k] = (-0.5*Is[3]*sqrt(Is[2])/pow(Is[1],1.5))*dI1ds[k] + (0.5*Is[3]/sqrt(Is[1]*Is[2]))*dI2ds[k] + sqrt(Is[2]/Is[1])*dI3ds[k]; 00269 if (qsmp>DBL_EPSILON*100.0) DqsmpDs[k] = (1.0/qsmp)*(L[k]-psmp*DpsmpDs[k]); 00270 else DqsmpDs[k] = 0.0; 00271 } 00272 } 00273 else throw new Fatal("Derivs_psmp_qsmp: I2=%f must be positive",Is[1]); 00274 } 00275 else throw new Fatal("Derivs_psmp_qsmp: The principal stresses must be positive. s1=%f, s2=%f, s3=%f",L[0],L[1],L[2]); 00276 } 00277 00278 inline void Stress_ssmp_tsmp (STensor2 const & Sig, double & ssmp, double & tsmp) 00279 { 00280 // Principal values 00281 double L[3]; 00282 Eigenvals (Sig, L); 00283 00284 // Invariants 00285 if (L[0]>0.0 && L[1]>0.0 && L[2]>0.0) 00286 { 00287 // Characteristic invariants 00288 double Is[4]; 00289 Is[0] = L[0]+L[1]+L[2]; 00290 Is[1] = L[0]*L[1]+L[1]*L[2]+L[2]*L[0]; 00291 Is[2] = L[0]*L[1]*L[2]; 00292 Is[3] = sqrt(L[0])+sqrt(L[1])+sqrt(L[2]); 00293 00294 // ssmp and tsmp invariants 00295 if (Is[1]>0.0) 00296 { 00297 ssmp = 3.0*Is[2]/Is[1]; 00298 double del = Is[0]*Is[2]/Is[1]-ssmp*ssmp; 00299 if (fabs(del)<DBL_EPSILON*100.0) del = 0.0; 00300 if (del>=0.0) tsmp = sqrt(del); 00301 else throw new Fatal("Stress_ssmp_tsmp: del=I1*I3/I2-ssmp*ssmp=%f must be positive. I1=%f, I2=%f, I3=%f, ssmp=%f",del,Is[0],Is[1],Is[2],ssmp); 00302 } 00303 else throw new Fatal("Stress_ssmp_tsmp: I2=%f must be positive",Is[1]); 00304 } 00305 else throw new Fatal("Stress_ssmp_tsmp: The principal stresses must be positive. s1=%f, s2=%f, s3=%f",L[0],L[1],L[2]); 00306 } 00307 00308 inline void Derivs_ssmp_tsmp (STensor2 const & Sig, double & ssmp, double & tsmp, double DssmpDs[3], double DtsmpDs[3], STensor2 Proj[3]) 00309 { 00310 // Eigenvalues and Eigenprojectors 00311 double L[3]; 00312 Eigenvp (Sig, L, Proj); 00313 00314 // Invariants and derivatives 00315 if (L[0]>0.0 && L[1]>0.0 && L[2]>0.0) 00316 { 00317 // Characteristic invariants 00318 double Is[4]; 00319 Is[0] = L[0]+L[1]+L[2]; 00320 Is[1] = L[0]*L[1]+L[1]*L[2]+L[2]*L[0]; 00321 Is[2] = L[0]*L[1]*L[2]; 00322 Is[3] = sqrt(L[0])+sqrt(L[1])+sqrt(L[2]); 00323 00324 // Derivatives and invariants 00325 if (Is[1]>0.0) 00326 { 00327 // Invariants 00328 ssmp = 3.0*Is[2]/Is[1]; 00329 double del = Is[0]*Is[2]/Is[1]-ssmp*ssmp; 00330 if (fabs(del)<DBL_EPSILON*100.0) del = 0.0; 00331 if (del>=0.0) tsmp = sqrt(del); 00332 else throw new Fatal("Derivs_ssmp_tsmp: del=I1*I3/I2-ssmp*ssmp=%f must be positive. I1=%f, I2=%f, I3=%f, ssmp=%f",del,Is[0],Is[1],Is[2],ssmp); 00333 00334 // Derivatives 00335 double dI0ds[3], dI1ds[3], dI2ds[3]; 00336 for (int k=0; k<3; ++k) 00337 { 00338 // Derivatives of the characteristic invariants 00339 dI0ds[k] = 1.0; 00340 dI1ds[k] = Is[0]-L[k]; 00341 dI2ds[k] = Is[2]/L[k]; 00342 00343 // Derivatives of the ssmp and tsmp invariants 00344 DssmpDs[k] = (-3.0*Is[2]/(Is[1]*Is[1]))*dI1ds[k] + (3.0/Is[1])*dI2ds[k]; 00345 if (tsmp>DBL_EPSILON*100.0) DtsmpDs[k] = (0.5/tsmp)*((Is[2]/Is[1])*dI0ds[k]-(Is[0]*Is[2]/(Is[1]*Is[1]))*dI1ds[k]+(Is[0]/Is[1])*dI2ds[k]-2.0*ssmp*DssmpDs[k]); 00346 else DtsmpDs[k] = 0.0; 00347 } 00348 } 00349 else throw new Fatal("Derivs_ssmp_tsmp: I2=%f must be positive",Is[1]); 00350 } 00351 else throw new Fatal("Derivs_ssmp_tsmp: The principal stresses must be positive. s1=%f, s2=%f, s3=%f",L[0],L[1],L[2]); 00352 } 00353 00354 inline void Derivs_p_q (STensor2 const & Sig, double & p, double & q, double DpDs[3], double DqDs[3], STensor2 Proj[3]) 00355 { 00356 // Eigenvalues and Eigenprojectors 00357 double L[3]; 00358 Eigenvp (Sig, L, Proj); 00359 00360 // Invariants 00361 p = (L[0]+L[1]+L[2])/3.0; 00362 q = sqrt((L[0]-L[1])*(L[0]-L[1])+(L[1]-L[2])*(L[1]-L[2])+(L[2]-L[0])*(L[2]-L[0]))/SQ2; 00363 00364 // Derivatives 00365 for (int k=0; k<3; ++k) 00366 { 00367 DpDs[k] = 1.0/3.0; 00368 if (q>DBL_EPSILON*100) DqDs[k] = 3.0*(L[k]-p)/(2.0*q); 00369 else DqDs[k] = 0.0; 00370 } 00371 } 00372 00373 inline void Stress_p_q_t (STensor2 const & Sig, double & p, double & q, double & t) 00374 { 00375 Stress_p_q (Sig,p,q); 00376 STensor2 S; S = Sig - SymI * p; 00377 if (q<=DBL_EPSILON*100.0) t = -1.0; // 3th==-90 => th=30 00378 else t = 27.0*Det(S)/(2.0*q*q*q); 00379 if (t>= 1.0) t = 1.0; 00380 if (t<=-1.0) t = -1.0; 00381 } 00382 00383 inline void Derivs_p_q_t (STensor2 const & Sig, double & p, double & q, double & t, double DpDs[3], double DqDs[3], double DtDs[3], STensor2 Proj[3]) 00384 { 00385 // Eigenvalues and Eigenprojectors 00386 double L[3]; 00387 Eigenvp (Sig, L, Proj); 00388 00389 // Invariants 00390 p = (L[0]+L[1]+L[2])/3.0; 00391 q = sqrt((L[0]-L[1])*(L[0]-L[1])+(L[1]-L[2])*(L[1]-L[2])+(L[2]-L[0])*(L[2]-L[0]))/SQ2; 00392 double S[3]; for (int k=0; k<3; ++k) S[k] = L[k]-p; 00393 if (q>DBL_EPSILON*100.0) t = 27.0*S[0]*S[1]*S[2]/(2.0*q*q*q); 00394 else t = 0.0; 00395 00396 // Derivatives 00397 double B[3]; 00398 B[0] = L[2]-L[1]; 00399 B[1] = L[0]-L[2]; 00400 B[2] = L[1]-L[0]; 00401 double l = (L[0]-L[1])*(L[1]-L[2])*(L[2]-L[0]); 00402 for (int k=0; k<3; ++k) 00403 { 00404 DpDs[k] = 1.0/3.0; 00405 if (q>DBL_EPSILON*100) 00406 { 00407 DqDs[k] = 3.0*S[k]/(2.0*q); 00408 DtDs[k] = 27.0*l*B[k]/(4.0*pow(q,5.0)); 00409 } 00410 else 00411 { 00412 DqDs[k] = 0.0; 00413 DtDs[k] = 0.0; 00414 } 00415 } 00416 } 00417 00418 inline void Strain_Ev_Ed (STensor2 const & Eps, double & Ev, double & Ed) 00419 { 00420 Ev = Eps(0)+Eps(1)+Eps(2); 00421 Ed = sqrt(2.0*( (Eps(0)-Eps(1))*(Eps(0)-Eps(1)) 00422 + (Eps(1)-Eps(2))*(Eps(1)-Eps(2)) 00423 + (Eps(2)-Eps(0))*(Eps(2)-Eps(0)) + 3.0*(Eps(3)*Eps(3) + Eps(4)*Eps(4) + Eps(5)*Eps(5)) ))/3.0; 00424 } 00425 00427 inline void OutputHeader (std::ostream & Os, bool NewLine=true) 00428 { 00429 Os << _8s<< "Sx" << _8s<< "Sy" << _8s<< "Sz" << _8s<< "Sxy" << _8s<< "Syz" << _8s<< "Sxz"; 00430 Os << _8s<< "Ex" << _8s<< "Ey" << _8s<< "Ez" << _8s<< "Exy" << _8s<< "Eyz" << _8s<< "Exz"; 00431 Os << _8s<< "p" << _8s<< "q" << _8s<< "Ev" << _8s<< "Ed"; 00432 if (NewLine) Os << std::endl; 00433 } 00434 00436 inline void Output (STensor2 const & Sig, STensor2 const & Eps, std::ostream & Os, bool NewLine=true) 00437 { 00438 double p,q,ev,ed; 00439 Stress_p_q (Sig,p,q); 00440 Strain_Ev_Ed (Eps,ev,ed); 00441 Os << _8s<< Sig(0) << _8s<< Sig(1) << _8s<< Sig(2) << _8s<< Sig(3)/SQ2 << _8s<< Sig(4)/SQ2 << _8s<< Sig(5)/SQ2; 00442 Os << _8s<< Eps(0)*100.0 << _8s<< Eps(1)*100.0 << _8s<< Eps(2)*100.0 << _8s<< Eps(3)*100.0/SQ2 << _8s<< Eps(4)*100.0/SQ2 << _8s<< Eps(5)*100.0/SQ2; 00443 Os << _8s<< p << _8s<< q << _8s<< ev*100.0 << _8s<< ed*100.0; 00444 if (NewLine) Os << std::endl; 00445 } 00446 00448 inline void Hid2Sig (Vector3D const & RAZ, Vector3D & XYZ) 00449 { 00450 double SA = -RAZ(0)*sin(RAZ(1)); 00451 double SB = RAZ(0)*cos(RAZ(1)); 00452 XYZ = RAZ(2) + 2.0*SB/SQ6 , 00453 RAZ(2) - SB/SQ6 - SA/SQ2, 00454 RAZ(2) - SB/SQ6 + SA/SQ2; 00455 } 00456 00457 00459 00461 00463 00465 inline void Dot (ATensor2 const & A, Vector3D const & u, Vector3D & v) 00466 { 00467 v(0) = A(0)*u(0)+A(3)*u(1)+A(5)*u(2); 00468 v(1) = A(6)*u(0)+A(1)*u(1)+A(4)*u(2); 00469 v(2) = A(8)*u(0)+A(7)*u(1)+A(2)*u(2); 00470 } 00471 00473 inline void Dot (Vector3D const & u, ATensor2 const & A, Vector3D & v) 00474 { 00475 v(0) = u(0)*A(0)+u(1)*A(6)+u(2)*A(8); 00476 v(1) = u(0)*A(3)+u(1)*A(1)+u(2)*A(7); 00477 v(2) = u(0)*A(5)+u(1)*A(4)+u(2)*A(2); 00478 } 00479 00481 inline void Dot (ATensor4 const & T, ATensor2 const & x, ATensor2 & y) 00482 { 00483 y(0) = T(0,0)*x(0)+T(0,3)*x(3)+T(0,5)*x(5)+T(0,6)*x(6)+T(0,1)*x(1)+T(0,4)*x(4)+T(0,8)*x(8)+T(0,7)*x(7)+T(0,2)*x(2); 00484 y(3) = T(3,0)*x(0)+T(3,3)*x(3)+T(3,5)*x(5)+T(3,6)*x(6)+T(3,1)*x(1)+T(3,4)*x(4)+T(3,8)*x(8)+T(3,7)*x(7)+T(3,2)*x(2); 00485 y(5) = T(5,0)*x(0)+T(5,3)*x(3)+T(5,5)*x(5)+T(5,6)*x(6)+T(5,1)*x(1)+T(5,4)*x(4)+T(5,8)*x(8)+T(5,7)*x(7)+T(5,2)*x(2); 00486 y(6) = T(6,0)*x(0)+T(6,3)*x(3)+T(6,5)*x(5)+T(6,6)*x(6)+T(6,1)*x(1)+T(6,4)*x(4)+T(6,8)*x(8)+T(6,7)*x(7)+T(6,2)*x(2); 00487 y(1) = T(1,0)*x(0)+T(1,3)*x(3)+T(1,5)*x(5)+T(1,6)*x(6)+T(1,1)*x(1)+T(1,4)*x(4)+T(1,8)*x(8)+T(1,7)*x(7)+T(1,2)*x(2); 00488 y(4) = T(4,0)*x(0)+T(4,3)*x(3)+T(4,5)*x(5)+T(4,6)*x(6)+T(4,1)*x(1)+T(4,4)*x(4)+T(4,8)*x(8)+T(4,7)*x(7)+T(4,2)*x(2); 00489 y(8) = T(8,0)*x(0)+T(8,3)*x(3)+T(8,5)*x(5)+T(8,6)*x(6)+T(8,1)*x(1)+T(8,4)*x(4)+T(8,8)*x(8)+T(8,7)*x(7)+T(8,2)*x(2); 00490 y(7) = T(7,0)*x(0)+T(7,3)*x(3)+T(7,5)*x(5)+T(7,6)*x(6)+T(7,1)*x(1)+T(7,4)*x(4)+T(7,8)*x(8)+T(7,7)*x(7)+T(7,2)*x(2); 00491 y(2) = T(2,0)*x(0)+T(2,3)*x(3)+T(2,5)*x(5)+T(2,6)*x(6)+T(2,1)*x(1)+T(2,4)*x(4)+T(2,8)*x(8)+T(2,7)*x(7)+T(2,2)*x(2); 00492 } 00493 00495 inline void Dot (ATensor2 const & A, ATensor2 const & B, ATensor2 & C) 00496 { 00497 C(0) = A(0)*B(0)+A(3)*B(6)+A(5)*B(8); 00498 C(3) = A(0)*B(3)+A(3)*B(1)+A(5)*B(7); 00499 C(5) = A(0)*B(5)+A(3)*B(4)+A(5)*B(2); 00500 C(6) = A(6)*B(0)+A(1)*B(6)+A(4)*B(8); 00501 C(1) = A(6)*B(3)+A(1)*B(1)+A(4)*B(7); 00502 C(4) = A(6)*B(5)+A(1)*B(4)+A(4)*B(2); 00503 C(8) = A(8)*B(0)+A(7)*B(6)+A(2)*B(8); 00504 C(7) = A(8)*B(3)+A(7)*B(1)+A(2)*B(7); 00505 C(2) = A(8)*B(5)+A(7)*B(4)+A(2)*B(2); 00506 } 00507 00509 inline void Dyad (Vector3D const & u, Vector3D const & v, ATensor2 & A) 00510 { 00511 A(0) = u(0)*v(0); 00512 A(3) = u(0)*v(1); 00513 A(5) = u(0)*v(2); 00514 A(6) = u(1)*v(0); 00515 A(1) = u(1)*v(1); 00516 A(4) = u(1)*v(2); 00517 A(8) = u(2)*v(0); 00518 A(7) = u(2)*v(1); 00519 A(2) = u(2)*v(2); 00520 } 00521 00523 inline void DyadDot (Vector3D const & v, ATensor2 const & F, Vector3D const & u, ATensor2 & R) 00524 { 00525 R(0) = v(0)*F(0)*u(0)+v(0)*F(3)*u(1)+v(0)*F(5)*u(2); 00526 R(3) = v(0)*F(6)*u(0)+v(0)*F(1)*u(1)+v(0)*F(4)*u(2); 00527 R(5) = v(0)*F(8)*u(0)+v(0)*F(7)*u(1)+v(0)*F(2)*u(2); 00528 R(6) = v(1)*F(0)*u(0)+v(1)*F(3)*u(1)+v(1)*F(5)*u(2); 00529 R(1) = v(1)*F(6)*u(0)+v(1)*F(1)*u(1)+v(1)*F(4)*u(2); 00530 R(4) = v(1)*F(8)*u(0)+v(1)*F(7)*u(1)+v(1)*F(2)*u(2); 00531 R(8) = v(2)*F(0)*u(0)+v(2)*F(3)*u(1)+v(2)*F(5)*u(2); 00532 R(7) = v(2)*F(6)*u(0)+v(2)*F(1)*u(1)+v(2)*F(4)*u(2); 00533 R(2) = v(2)*F(8)*u(0)+v(2)*F(7)*u(1)+v(2)*F(2)*u(2); 00534 } 00535 00537 inline void DyadDot (Vector3D const & v, Vector3D const & u, ATensor2 const & F, ATensor2 & R) 00538 { 00539 R(0) = v(0)*u(0)*F(0)+v(0)*u(1)*F(6)+v(0)*u(2)*F(8); 00540 R(3) = v(0)*u(0)*F(3)+v(0)*u(1)*F(1)+v(0)*u(2)*F(7); 00541 R(5) = v(0)*u(0)*F(5)+v(0)*u(1)*F(4)+v(0)*u(2)*F(2); 00542 R(6) = v(1)*u(0)*F(0)+v(1)*u(1)*F(6)+v(1)*u(2)*F(8); 00543 R(1) = v(1)*u(0)*F(3)+v(1)*u(1)*F(1)+v(1)*u(2)*F(7); 00544 R(4) = v(1)*u(0)*F(5)+v(1)*u(1)*F(4)+v(1)*u(2)*F(2); 00545 R(8) = v(2)*u(0)*F(0)+v(2)*u(1)*F(6)+v(2)*u(2)*F(8); 00546 R(7) = v(2)*u(0)*F(3)+v(2)*u(1)*F(1)+v(2)*u(2)*F(7); 00547 R(2) = v(2)*u(0)*F(5)+v(2)*u(1)*F(4)+v(2)*u(2)*F(2); 00548 } 00549 00551 00553 inline void DotUp (ATensor2 const & A, Vector3D const & u, Vector3D & v) 00554 { 00555 v(0) += A(0)*u(0)+A(3)*u(1)+A(5)*u(2); 00556 v(1) += A(6)*u(0)+A(1)*u(1)+A(4)*u(2); 00557 v(2) += A(8)*u(0)+A(7)*u(1)+A(2)*u(2); 00558 } 00559 00561 inline void DotUp (Vector3D const & u, ATensor2 const & A, Vector3D & v) 00562 { 00563 v(0) += u(0)*A(0)+u(1)*A(6)+u(2)*A(8); 00564 v(1) += u(0)*A(3)+u(1)*A(1)+u(2)*A(7); 00565 v(2) += u(0)*A(5)+u(1)*A(4)+u(2)*A(2); 00566 } 00567 00569 inline void DotUp (ATensor4 const & T, ATensor2 const & x, ATensor2 & y) 00570 { 00571 y(0) += T(0,0)*x(0)+T(0,3)*x(3)+T(0,5)*x(5)+T(0,6)*x(6)+T(0,1)*x(1)+T(0,4)*x(4)+T(0,8)*x(8)+T(0,7)*x(7)+T(0,2)*x(2); 00572 y(3) += T(3,0)*x(0)+T(3,3)*x(3)+T(3,5)*x(5)+T(3,6)*x(6)+T(3,1)*x(1)+T(3,4)*x(4)+T(3,8)*x(8)+T(3,7)*x(7)+T(3,2)*x(2); 00573 y(5) += T(5,0)*x(0)+T(5,3)*x(3)+T(5,5)*x(5)+T(5,6)*x(6)+T(5,1)*x(1)+T(5,4)*x(4)+T(5,8)*x(8)+T(5,7)*x(7)+T(5,2)*x(2); 00574 y(6) += T(6,0)*x(0)+T(6,3)*x(3)+T(6,5)*x(5)+T(6,6)*x(6)+T(6,1)*x(1)+T(6,4)*x(4)+T(6,8)*x(8)+T(6,7)*x(7)+T(6,2)*x(2); 00575 y(1) += T(1,0)*x(0)+T(1,3)*x(3)+T(1,5)*x(5)+T(1,6)*x(6)+T(1,1)*x(1)+T(1,4)*x(4)+T(1,8)*x(8)+T(1,7)*x(7)+T(1,2)*x(2); 00576 y(4) += T(4,0)*x(0)+T(4,3)*x(3)+T(4,5)*x(5)+T(4,6)*x(6)+T(4,1)*x(1)+T(4,4)*x(4)+T(4,8)*x(8)+T(4,7)*x(7)+T(4,2)*x(2); 00577 y(8) += T(8,0)*x(0)+T(8,3)*x(3)+T(8,5)*x(5)+T(8,6)*x(6)+T(8,1)*x(1)+T(8,4)*x(4)+T(8,8)*x(8)+T(8,7)*x(7)+T(8,2)*x(2); 00578 y(7) += T(7,0)*x(0)+T(7,3)*x(3)+T(7,5)*x(5)+T(7,6)*x(6)+T(7,1)*x(1)+T(7,4)*x(4)+T(7,8)*x(8)+T(7,7)*x(7)+T(7,2)*x(2); 00579 y(2) += T(2,0)*x(0)+T(2,3)*x(3)+T(2,5)*x(5)+T(2,6)*x(6)+T(2,1)*x(1)+T(2,4)*x(4)+T(2,8)*x(8)+T(2,7)*x(7)+T(2,2)*x(2); 00580 } 00581 00583 inline void DyadUp (Vector3D const & u, Vector3D const & v, ATensor2 & A) 00584 { 00585 A(0) += u(0)*v(0); 00586 A(3) += u(0)*v(1); 00587 A(5) += u(0)*v(2); 00588 A(6) += u(1)*v(0); 00589 A(1) += u(1)*v(1); 00590 A(4) += u(1)*v(2); 00591 A(8) += u(2)*v(0); 00592 A(7) += u(2)*v(1); 00593 A(2) += u(2)*v(2); 00594 } 00595 00597 inline void ScDyadUp (double s, Vector3D const & u, Vector3D const & v, ATensor2 & A) 00598 { 00599 A(0) += s*u(0)*v(0); 00600 A(3) += s*u(0)*v(1); 00601 A(5) += s*u(0)*v(2); 00602 A(6) += s*u(1)*v(0); 00603 A(1) += s*u(1)*v(1); 00604 A(4) += s*u(1)*v(2); 00605 A(8) += s*u(2)*v(0); 00606 A(7) += s*u(2)*v(1); 00607 A(2) += s*u(2)*v(2); 00608 } 00609 00611 inline void DyadDotUp (Vector3D const & v, ATensor2 const & F, Vector3D const & u, ATensor2 & R) 00612 { 00613 R(0) += v(0)*F(0)*u(0)+v(0)*F(3)*u(1)+v(0)*F(5)*u(2); 00614 R(3) += v(0)*F(6)*u(0)+v(0)*F(1)*u(1)+v(0)*F(4)*u(2); 00615 R(5) += v(0)*F(8)*u(0)+v(0)*F(7)*u(1)+v(0)*F(2)*u(2); 00616 R(6) += v(1)*F(0)*u(0)+v(1)*F(3)*u(1)+v(1)*F(5)*u(2); 00617 R(1) += v(1)*F(6)*u(0)+v(1)*F(1)*u(1)+v(1)*F(4)*u(2); 00618 R(4) += v(1)*F(8)*u(0)+v(1)*F(7)*u(1)+v(1)*F(2)*u(2); 00619 R(8) += v(2)*F(0)*u(0)+v(2)*F(3)*u(1)+v(2)*F(5)*u(2); 00620 R(7) += v(2)*F(6)*u(0)+v(2)*F(1)*u(1)+v(2)*F(4)*u(2); 00621 R(2) += v(2)*F(8)*u(0)+v(2)*F(7)*u(1)+v(2)*F(2)*u(2); 00622 } 00623 00625 inline void DyadDotUp (Vector3D const & v, Vector3D const & u, ATensor2 const & F, ATensor2 & R) 00626 { 00627 R(0) += v(0)*u(0)*F(0)+v(0)*u(1)*F(6)+v(0)*u(2)*F(8); 00628 R(3) += v(0)*u(0)*F(3)+v(0)*u(1)*F(1)+v(0)*u(2)*F(7); 00629 R(5) += v(0)*u(0)*F(5)+v(0)*u(1)*F(4)+v(0)*u(2)*F(2); 00630 R(6) += v(1)*u(0)*F(0)+v(1)*u(1)*F(6)+v(1)*u(2)*F(8); 00631 R(1) += v(1)*u(0)*F(3)+v(1)*u(1)*F(1)+v(1)*u(2)*F(7); 00632 R(4) += v(1)*u(0)*F(5)+v(1)*u(1)*F(4)+v(1)*u(2)*F(2); 00633 R(8) += v(2)*u(0)*F(0)+v(2)*u(1)*F(6)+v(2)*u(2)*F(8); 00634 R(7) += v(2)*u(0)*F(3)+v(2)*u(1)*F(1)+v(2)*u(2)*F(7); 00635 R(2) += v(2)*u(0)*F(5)+v(2)*u(1)*F(4)+v(2)*u(2)*F(2); 00636 } 00637 00639 00641 00643 inline void Dot (STensor2 const & A, Vector3D const & u, Vector3D & v) 00644 { 00645 v(0) = A(0)*u(0) + A(3)*u(1)/SQ2 + A(5)*u(2)/SQ2; 00646 v(1) = A(3)*u(0)/SQ2 + A(1)*u(1) + A(4)*u(2)/SQ2; 00647 v(2) = A(5)*u(0)/SQ2 + A(4)*u(1)/SQ2 + A(2)*u(2); 00648 } 00649 00651 inline void Dot (Vector3D const & u, STensor2 const & A, Vector3D & v) 00652 { 00653 v(0) = u(0)*A(0) + u(1)*A(3)/SQ2 + u(2)*A(5)/SQ2; 00654 v(1) = u(0)*A(3)/SQ2 + u(1)*A(1) + u(2)*A(4)/SQ2; 00655 v(2) = u(0)*A(5)/SQ2 + u(1)*A(4)/SQ2 + u(2)*A(2); 00656 } 00657 00659 inline void Dot (STensor4 const & T, STensor2 const & x, STensor2 & y) 00660 { 00661 y(0) = T(0,0)*x(0)+T(0,1)*x(1)+T(0,2)*x(2)+T(0,3)*x(3)+T(0,4)*x(4)+T(0,5)*x(5); 00662 y(1) = T(1,0)*x(0)+T(1,1)*x(1)+T(1,2)*x(2)+T(1,3)*x(3)+T(1,4)*x(4)+T(1,5)*x(5); 00663 y(2) = T(2,0)*x(0)+T(2,1)*x(1)+T(2,2)*x(2)+T(2,3)*x(3)+T(2,4)*x(4)+T(2,5)*x(5); 00664 y(3) = T(3,0)*x(0)+T(3,1)*x(1)+T(3,2)*x(2)+T(3,3)*x(3)+T(3,4)*x(4)+T(3,5)*x(5); 00665 y(4) = T(4,0)*x(0)+T(4,1)*x(1)+T(4,2)*x(2)+T(4,3)*x(3)+T(4,4)*x(4)+T(4,5)*x(5); 00666 y(5) = T(5,0)*x(0)+T(5,1)*x(1)+T(5,2)*x(2)+T(5,3)*x(3)+T(5,4)*x(4)+T(5,5)*x(5); 00667 } 00668 00670 inline void ScDot (double s, STensor2 const & x, STensor4 const & T, STensor2 & y) 00671 { 00672 y(0) = s*(x(0)*T(0,0)+x(1)*T(1,0)+x(2)*T(2,0)+x(3)*T(3,0)+x(4)*T(4,0)+x(5)*T(5,0)); 00673 y(1) = s*(x(0)*T(0,1)+x(1)*T(1,1)+x(2)*T(2,1)+x(3)*T(3,1)+x(4)*T(4,1)+x(5)*T(5,1)); 00674 y(2) = s*(x(0)*T(0,2)+x(1)*T(1,2)+x(2)*T(2,2)+x(3)*T(3,2)+x(4)*T(4,2)+x(5)*T(5,2)); 00675 y(3) = s*(x(0)*T(0,3)+x(1)*T(1,3)+x(2)*T(2,3)+x(3)*T(3,3)+x(4)*T(4,3)+x(5)*T(5,3)); 00676 y(4) = s*(x(0)*T(0,4)+x(1)*T(1,4)+x(2)*T(2,4)+x(3)*T(3,4)+x(4)*T(4,4)+x(5)*T(5,4)); 00677 y(5) = s*(x(0)*T(0,5)+x(1)*T(1,5)+x(2)*T(2,5)+x(3)*T(3,5)+x(4)*T(4,5)+x(5)*T(5,5)); 00678 } 00679 00681 00683 inline void DotUp (STensor2 const & A, Vector3D const & u, Vector3D & v) 00684 { 00685 v(0) += A(0)*u(0) + A(3)*u(1)/SQ2 + A(5)*u(2)/SQ2; 00686 v(1) += A(3)*u(0)/SQ2 + A(1)*u(1) + A(4)*u(2)/SQ2; 00687 v(2) += A(5)*u(0)/SQ2 + A(4)*u(1)/SQ2 + A(2)*u(2); 00688 } 00689 00691 inline void ScDotUp (double s, STensor2 const & A, Vector3D const & u, Vector3D & v) 00692 { 00693 v(0) += s*( A(0)*u(0) + A(3)*u(1)/SQ2 + A(5)*u(2)/SQ2 ); 00694 v(1) += s*( A(3)*u(0)/SQ2 + A(1)*u(1) + A(4)*u(2)/SQ2 ); 00695 v(2) += s*( A(5)*u(0)/SQ2 + A(4)*u(1)/SQ2 + A(2)*u(2) ); 00696 } 00697 00699 inline void ScDotUp (double s, Vec_t const & A, Vector3D const & u, Vector3D & v) 00700 { 00701 if (size(A)==4) 00702 { 00703 v(0) += s*( A(0)*u(0) + A(3)*u(1)/SQ2 ); 00704 v(1) += s*( A(3)*u(0)/SQ2 + A(1)*u(1) ); 00705 v(2) += s*( + A(2)*u(2) ); 00706 } 00707 else 00708 { 00709 v(0) += s*( A(0)*u(0) + A(3)*u(1)/SQ2 + A(5)*u(2)/SQ2 ); 00710 v(1) += s*( A(3)*u(0)/SQ2 + A(1)*u(1) + A(4)*u(2)/SQ2 ); 00711 v(2) += s*( A(5)*u(0)/SQ2 + A(4)*u(1)/SQ2 + A(2)*u(2) ); 00712 } 00713 } 00714 00716 inline void DotUp (Vector3D const & u, STensor2 const & A, Vector3D & v) 00717 { 00718 v(0) += u(0)*A(0) + u(1)*A(3)/SQ2 + u(2)*A(5)/SQ2; 00719 v(1) += u(0)*A(3)/SQ2 + u(1)*A(1) + u(2)*A(4)/SQ2; 00720 v(2) += u(0)*A(5)/SQ2 + u(1)*A(4)/SQ2 + u(2)*A(2); 00721 } 00722 00724 inline void ScDotUp (double s, Vector3D const & u, STensor2 const & A, Vector3D & v) 00725 { 00726 v(0) += s*( u(0)*A(0) + u(1)*A(3)/SQ2 + u(2)*A(5)/SQ2 ); 00727 v(1) += s*( u(0)*A(3)/SQ2 + u(1)*A(1) + u(2)*A(4)/SQ2 ); 00728 v(2) += s*( u(0)*A(5)/SQ2 + u(1)*A(4)/SQ2 + u(2)*A(2) ); 00729 } 00730 00732 inline void DotUp (STensor4 const & T, STensor2 const & x, STensor2 & y) 00733 { 00734 y(0) += T(0,0)*x(0)+T(0,1)*x(1)+T(0,2)*x(2)+T(0,3)*x(3)+T(0,4)*x(4)+T(0,5)*x(5); 00735 y(1) += T(1,0)*x(0)+T(1,1)*x(1)+T(1,2)*x(2)+T(1,3)*x(3)+T(1,4)*x(4)+T(1,5)*x(5); 00736 y(2) += T(2,0)*x(0)+T(2,1)*x(1)+T(2,2)*x(2)+T(2,3)*x(3)+T(2,4)*x(4)+T(2,5)*x(5); 00737 y(3) += T(3,0)*x(0)+T(3,1)*x(1)+T(3,2)*x(2)+T(3,3)*x(3)+T(3,4)*x(4)+T(3,5)*x(5); 00738 y(4) += T(4,0)*x(0)+T(4,1)*x(1)+T(4,2)*x(2)+T(4,3)*x(3)+T(4,4)*x(4)+T(4,5)*x(5); 00739 y(5) += T(5,0)*x(0)+T(5,1)*x(1)+T(5,2)*x(2)+T(5,3)*x(3)+T(5,4)*x(4)+T(5,5)*x(5); 00740 } 00741 00743 inline void DotAdd (STensor2 const & r, STensor4 const & T, STensor2 const & x, STensor2 & y) 00744 { 00745 y(0) = r(0) + T(0,0)*x(0)+T(0,1)*x(1)+T(0,2)*x(2)+T(0,3)*x(3)+T(0,4)*x(4)+T(0,5)*x(5); 00746 y(1) = r(1) + T(1,0)*x(0)+T(1,1)*x(1)+T(1,2)*x(2)+T(1,3)*x(3)+T(1,4)*x(4)+T(1,5)*x(5); 00747 y(2) = r(2) + T(2,0)*x(0)+T(2,1)*x(1)+T(2,2)*x(2)+T(2,3)*x(3)+T(2,4)*x(4)+T(2,5)*x(5); 00748 y(3) = r(3) + T(3,0)*x(0)+T(3,1)*x(1)+T(3,2)*x(2)+T(3,3)*x(3)+T(3,4)*x(4)+T(3,5)*x(5); 00749 y(4) = r(4) + T(4,0)*x(0)+T(4,1)*x(1)+T(4,2)*x(2)+T(4,3)*x(3)+T(4,4)*x(4)+T(4,5)*x(5); 00750 y(5) = r(5) + T(5,0)*x(0)+T(5,1)*x(1)+T(5,2)*x(2)+T(5,3)*x(3)+T(5,4)*x(4)+T(5,5)*x(5); 00751 } 00752 00753 00754 // 14) s = xt * A * y 00755 inline double Reduce (STensor2 const & x, STensor4 const & A, STensor2 const & y) 00756 { 00757 return x(0)*(A(0,0)*y(0) + A(0,1)*y(1) + A(0,2)*y(2) + A(0,3)*y(3) + A(0,4)*y(4) + A(0,5)*y(5)) + 00758 x(1)*(A(1,0)*y(0) + A(1,1)*y(1) + A(1,2)*y(2) + A(1,3)*y(3) + A(1,4)*y(4) + A(1,5)*y(5)) + 00759 x(2)*(A(2,0)*y(0) + A(2,1)*y(1) + A(2,2)*y(2) + A(2,3)*y(3) + A(2,4)*y(4) + A(2,5)*y(5)) + 00760 x(3)*(A(3,0)*y(0) + A(3,1)*y(1) + A(3,2)*y(2) + A(3,3)*y(3) + A(3,4)*y(4) + A(3,5)*y(5)) + 00761 x(4)*(A(4,0)*y(0) + A(4,1)*y(1) + A(4,2)*y(2) + A(4,3)*y(3) + A(4,4)*y(4) + A(4,5)*y(5)) + 00762 x(5)*(A(5,0)*y(0) + A(5,1)*y(1) + A(5,2)*y(2) + A(5,3)*y(3) + A(5,4)*y(4) + A(5,5)*y(5)); 00763 } 00764 00765 // 15) D = a * (A*x) dyadic (y*B) + C 00766 inline void Generate (double const & a, STensor4 const & A, STensor2 const & x, STensor2 const & y, STensor4 const & B, STensor4 const & C, STensor4 & D) 00767 { 00768 //C(i,j)=a* (A(i,k)*x(k)) * (y(l)*B(l,j)) +D(i,j); 00769 00770 D(0,0)=a*(A(0,0)*x(0)+A(0,1)*x(1)+A(0,2)*x(2)+A(0,3)*x(3)+A(0,4)*x(4)+A(0,5)*x(5)) * (y(0)*B(0,0)+y(1)*B(1,0)+y(2)*B(2,0)+y(3)*B(3,0)+y(4)*B(4,0)+y(5)*B(5,0)) + C(0,0); 00771 D(0,1)=a*(A(0,0)*x(0)+A(0,1)*x(1)+A(0,2)*x(2)+A(0,3)*x(3)+A(0,4)*x(4)+A(0,5)*x(5)) * (y(0)*B(0,1)+y(1)*B(1,1)+y(2)*B(2,1)+y(3)*B(3,1)+y(4)*B(4,1)+y(5)*B(5,1)) + C(0,1); 00772 D(0,2)=a*(A(0,0)*x(0)+A(0,1)*x(1)+A(0,2)*x(2)+A(0,3)*x(3)+A(0,4)*x(4)+A(0,5)*x(5)) * (y(0)*B(0,2)+y(1)*B(1,2)+y(2)*B(2,2)+y(3)*B(3,2)+y(4)*B(4,2)+y(5)*B(5,2)) + C(0,2); 00773 D(0,3)=a*(A(0,0)*x(0)+A(0,1)*x(1)+A(0,2)*x(2)+A(0,3)*x(3)+A(0,4)*x(4)+A(0,5)*x(5)) * (y(0)*B(0,3)+y(1)*B(1,3)+y(2)*B(2,3)+y(3)*B(3,3)+y(4)*B(4,3)+y(5)*B(5,3)) + C(0,3); 00774 D(0,4)=a*(A(0,0)*x(0)+A(0,1)*x(1)+A(0,2)*x(2)+A(0,3)*x(3)+A(0,4)*x(4)+A(0,5)*x(5)) * (y(0)*B(0,4)+y(1)*B(1,4)+y(2)*B(2,4)+y(3)*B(3,4)+y(4)*B(4,4)+y(5)*B(5,4)) + C(0,4); 00775 D(0,5)=a*(A(0,0)*x(0)+A(0,1)*x(1)+A(0,2)*x(2)+A(0,3)*x(3)+A(0,4)*x(4)+A(0,5)*x(5)) * (y(0)*B(0,5)+y(1)*B(1,5)+y(2)*B(2,5)+y(3)*B(3,5)+y(4)*B(4,5)+y(5)*B(5,5)) + C(0,5); 00776 00777 D(1,0)=a*(A(1,0)*x(0)+A(1,1)*x(1)+A(1,2)*x(2)+A(1,3)*x(3)+A(1,4)*x(4)+A(1,5)*x(5)) * (y(0)*B(0,0)+y(1)*B(1,0)+y(2)*B(2,0)+y(3)*B(3,0)+y(4)*B(4,0)+y(5)*B(5,0)) + C(1,0); 00778 D(1,1)=a*(A(1,0)*x(0)+A(1,1)*x(1)+A(1,2)*x(2)+A(1,3)*x(3)+A(1,4)*x(4)+A(1,5)*x(5)) * (y(0)*B(0,1)+y(1)*B(1,1)+y(2)*B(2,1)+y(3)*B(3,1)+y(4)*B(4,1)+y(5)*B(5,1)) + C(1,1); 00779 D(1,2)=a*(A(1,0)*x(0)+A(1,1)*x(1)+A(1,2)*x(2)+A(1,3)*x(3)+A(1,4)*x(4)+A(1,5)*x(5)) * (y(0)*B(0,2)+y(1)*B(1,2)+y(2)*B(2,2)+y(3)*B(3,2)+y(4)*B(4,2)+y(5)*B(5,2)) + C(1,2); 00780 D(1,3)=a*(A(1,0)*x(0)+A(1,1)*x(1)+A(1,2)*x(2)+A(1,3)*x(3)+A(1,4)*x(4)+A(1,5)*x(5)) * (y(0)*B(0,3)+y(1)*B(1,3)+y(2)*B(2,3)+y(3)*B(3,3)+y(4)*B(4,3)+y(5)*B(5,3)) + C(1,3); 00781 D(1,4)=a*(A(1,0)*x(0)+A(1,1)*x(1)+A(1,2)*x(2)+A(1,3)*x(3)+A(1,4)*x(4)+A(1,5)*x(5)) * (y(0)*B(0,4)+y(1)*B(1,4)+y(2)*B(2,4)+y(3)*B(3,4)+y(4)*B(4,4)+y(5)*B(5,4)) + C(1,4); 00782 D(1,5)=a*(A(1,0)*x(0)+A(1,1)*x(1)+A(1,2)*x(2)+A(1,3)*x(3)+A(1,4)*x(4)+A(1,5)*x(5)) * (y(0)*B(0,5)+y(1)*B(1,5)+y(2)*B(2,5)+y(3)*B(3,5)+y(4)*B(4,5)+y(5)*B(5,5)) + C(1,5); 00783 00784 D(2,0)=a*(A(2,0)*x(0)+A(2,1)*x(1)+A(2,2)*x(2)+A(2,3)*x(3)+A(2,4)*x(4)+A(2,5)*x(5)) * (y(0)*B(0,0)+y(1)*B(1,0)+y(2)*B(2,0)+y(3)*B(3,0)+y(4)*B(4,0)+y(5)*B(5,0)) + C(2,0); 00785 D(2,1)=a*(A(2,0)*x(0)+A(2,1)*x(1)+A(2,2)*x(2)+A(2,3)*x(3)+A(2,4)*x(4)+A(2,5)*x(5)) * (y(0)*B(0,1)+y(1)*B(1,1)+y(2)*B(2,1)+y(3)*B(3,1)+y(4)*B(4,1)+y(5)*B(5,1)) + C(2,1); 00786 D(2,2)=a*(A(2,0)*x(0)+A(2,1)*x(1)+A(2,2)*x(2)+A(2,3)*x(3)+A(2,4)*x(4)+A(2,5)*x(5)) * (y(0)*B(0,2)+y(1)*B(1,2)+y(2)*B(2,2)+y(3)*B(3,2)+y(4)*B(4,2)+y(5)*B(5,2)) + C(2,2); 00787 D(2,3)=a*(A(2,0)*x(0)+A(2,1)*x(1)+A(2,2)*x(2)+A(2,3)*x(3)+A(2,4)*x(4)+A(2,5)*x(5)) * (y(0)*B(0,3)+y(1)*B(1,3)+y(2)*B(2,3)+y(3)*B(3,3)+y(4)*B(4,3)+y(5)*B(5,3)) + C(2,3); 00788 D(2,4)=a*(A(2,0)*x(0)+A(2,1)*x(1)+A(2,2)*x(2)+A(2,3)*x(3)+A(2,4)*x(4)+A(2,5)*x(5)) * (y(0)*B(0,4)+y(1)*B(1,4)+y(2)*B(2,4)+y(3)*B(3,4)+y(4)*B(4,4)+y(5)*B(5,4)) + C(2,4); 00789 D(2,5)=a*(A(2,0)*x(0)+A(2,1)*x(1)+A(2,2)*x(2)+A(2,3)*x(3)+A(2,4)*x(4)+A(2,5)*x(5)) * (y(0)*B(0,5)+y(1)*B(1,5)+y(2)*B(2,5)+y(3)*B(3,5)+y(4)*B(4,5)+y(5)*B(5,5)) + C(2,5); 00790 00791 D(3,0)=a*(A(3,0)*x(0)+A(3,1)*x(1)+A(3,2)*x(2)+A(3,3)*x(3)+A(3,4)*x(4)+A(3,5)*x(5)) * (y(0)*B(0,0)+y(1)*B(1,0)+y(2)*B(2,0)+y(3)*B(3,0)+y(4)*B(4,0)+y(5)*B(5,0)) + C(3,0); 00792 D(3,1)=a*(A(3,0)*x(0)+A(3,1)*x(1)+A(3,2)*x(2)+A(3,3)*x(3)+A(3,4)*x(4)+A(3,5)*x(5)) * (y(0)*B(0,1)+y(1)*B(1,1)+y(2)*B(2,1)+y(3)*B(3,1)+y(4)*B(4,1)+y(5)*B(5,1)) + C(3,1); 00793 D(3,2)=a*(A(3,0)*x(0)+A(3,1)*x(1)+A(3,2)*x(2)+A(3,3)*x(3)+A(3,4)*x(4)+A(3,5)*x(5)) * (y(0)*B(0,2)+y(1)*B(1,2)+y(2)*B(2,2)+y(3)*B(3,2)+y(4)*B(4,2)+y(5)*B(5,2)) + C(3,2); 00794 D(3,3)=a*(A(3,0)*x(0)+A(3,1)*x(1)+A(3,2)*x(2)+A(3,3)*x(3)+A(3,4)*x(4)+A(3,5)*x(5)) * (y(0)*B(0,3)+y(1)*B(1,3)+y(2)*B(2,3)+y(3)*B(3,3)+y(4)*B(4,3)+y(5)*B(5,3)) + C(3,3); 00795 D(3,4)=a*(A(3,0)*x(0)+A(3,1)*x(1)+A(3,2)*x(2)+A(3,3)*x(3)+A(3,4)*x(4)+A(3,5)*x(5)) * (y(0)*B(0,4)+y(1)*B(1,4)+y(2)*B(2,4)+y(3)*B(3,4)+y(4)*B(4,4)+y(5)*B(5,4)) + C(3,4); 00796 D(3,5)=a*(A(3,0)*x(0)+A(3,1)*x(1)+A(3,2)*x(2)+A(3,3)*x(3)+A(3,4)*x(4)+A(3,5)*x(5)) * (y(0)*B(0,5)+y(1)*B(1,5)+y(2)*B(2,5)+y(3)*B(3,5)+y(4)*B(4,5)+y(5)*B(5,5)) + C(3,5); 00797 00798 D(4,0)=a*(A(4,0)*x(0)+A(4,1)*x(1)+A(4,2)*x(2)+A(4,3)*x(3)+A(4,4)*x(4)+A(4,5)*x(5)) * (y(0)*B(0,0)+y(1)*B(1,0)+y(2)*B(2,0)+y(3)*B(3,0)+y(4)*B(4,0)+y(5)*B(5,0)) + C(4,0); 00799 D(4,1)=a*(A(4,0)*x(0)+A(4,1)*x(1)+A(4,2)*x(2)+A(4,3)*x(3)+A(4,4)*x(4)+A(4,5)*x(5)) * (y(0)*B(0,1)+y(1)*B(1,1)+y(2)*B(2,1)+y(3)*B(3,1)+y(4)*B(4,1)+y(5)*B(5,1)) + C(4,1); 00800 D(4,2)=a*(A(4,0)*x(0)+A(4,1)*x(1)+A(4,2)*x(2)+A(4,3)*x(3)+A(4,4)*x(4)+A(4,5)*x(5)) * (y(0)*B(0,2)+y(1)*B(1,2)+y(2)*B(2,2)+y(3)*B(3,2)+y(4)*B(4,2)+y(5)*B(5,2)) + C(4,2); 00801 D(4,3)=a*(A(4,0)*x(0)+A(4,1)*x(1)+A(4,2)*x(2)+A(4,3)*x(3)+A(4,4)*x(4)+A(4,5)*x(5)) * (y(0)*B(0,3)+y(1)*B(1,3)+y(2)*B(2,3)+y(3)*B(3,3)+y(4)*B(4,3)+y(5)*B(5,3)) + C(4,3); 00802 D(4,4)=a*(A(4,0)*x(0)+A(4,1)*x(1)+A(4,2)*x(2)+A(4,3)*x(3)+A(4,4)*x(4)+A(4,5)*x(5)) * (y(0)*B(0,4)+y(1)*B(1,4)+y(2)*B(2,4)+y(3)*B(3,4)+y(4)*B(4,4)+y(5)*B(5,4)) + C(4,4); 00803 D(4,5)=a*(A(4,0)*x(0)+A(4,1)*x(1)+A(4,2)*x(2)+A(4,3)*x(3)+A(4,4)*x(4)+A(4,5)*x(5)) * (y(0)*B(0,5)+y(1)*B(1,5)+y(2)*B(2,5)+y(3)*B(3,5)+y(4)*B(4,5)+y(5)*B(5,5)) + C(4,5); 00804 00805 D(5,0)=a*(A(5,0)*x(0)+A(5,1)*x(1)+A(5,2)*x(2)+A(5,3)*x(3)+A(5,4)*x(4)+A(5,5)*x(5)) * (y(0)*B(0,0)+y(1)*B(1,0)+y(2)*B(2,0)+y(3)*B(3,0)+y(4)*B(4,0)+y(5)*B(5,0)) + C(5,0); 00806 D(5,1)=a*(A(5,0)*x(0)+A(5,1)*x(1)+A(5,2)*x(2)+A(5,3)*x(3)+A(5,4)*x(4)+A(5,5)*x(5)) * (y(0)*B(0,1)+y(1)*B(1,1)+y(2)*B(2,1)+y(3)*B(3,1)+y(4)*B(4,1)+y(5)*B(5,1)) + C(5,1); 00807 D(5,2)=a*(A(5,0)*x(0)+A(5,1)*x(1)+A(5,2)*x(2)+A(5,3)*x(3)+A(5,4)*x(4)+A(5,5)*x(5)) * (y(0)*B(0,2)+y(1)*B(1,2)+y(2)*B(2,2)+y(3)*B(3,2)+y(4)*B(4,2)+y(5)*B(5,2)) + C(5,2); 00808 D(5,3)=a*(A(5,0)*x(0)+A(5,1)*x(1)+A(5,2)*x(2)+A(5,3)*x(3)+A(5,4)*x(4)+A(5,5)*x(5)) * (y(0)*B(0,3)+y(1)*B(1,3)+y(2)*B(2,3)+y(3)*B(3,3)+y(4)*B(4,3)+y(5)*B(5,3)) + C(5,3); 00809 D(5,4)=a*(A(5,0)*x(0)+A(5,1)*x(1)+A(5,2)*x(2)+A(5,3)*x(3)+A(5,4)*x(4)+A(5,5)*x(5)) * (y(0)*B(0,4)+y(1)*B(1,4)+y(2)*B(2,4)+y(3)*B(3,4)+y(4)*B(4,4)+y(5)*B(5,4)) + C(5,4); 00810 D(5,5)=a*(A(5,0)*x(0)+A(5,1)*x(1)+A(5,2)*x(2)+A(5,3)*x(3)+A(5,4)*x(4)+A(5,5)*x(5)) * (y(0)*B(0,5)+y(1)*B(1,5)+y(2)*B(2,5)+y(3)*B(3,5)+y(4)*B(4,5)+y(5)*B(5,5)) + C(5,5); 00811 } 00812 00813 00815 00817 int __initialize_tensors () 00818 { 00819 SymI = 1.0, 1.0, 1.0, 0.0, 0.0, 0.0; 00820 SymII = 1.0, 0.0, 0.0, 0.0, 0.0, 0.0, 00821 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 00822 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 00823 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 00824 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 00825 0.0, 0.0, 0.0, 0.0, 0.0, 1.0; 00826 SymIdyI = 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 00827 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 00828 1.0, 1.0, 1.0, 0.0, 0.0, 0.0, 00829 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 00830 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 00831 0.0, 0.0, 0.0, 0.0, 0.0, 0.0; 00832 SymPsd = 2.0/3.0, -1.0/3.0, -1.0/3.0, 0.0, 0.0, 0.0, 00833 -1.0/3.0, 2.0/3.0, -1.0/3.0, 0.0, 0.0, 0.0, 00834 -1.0/3.0, -1.0/3.0, 2.0/3.0, 0.0, 0.0, 0.0, 00835 0.0, 0.0, 0.0, 1.0, 0.0, 0.0, 00836 0.0, 0.0, 0.0, 0.0, 1.0, 0.0, 00837 0.0, 0.0, 0.0, 0.0, 0.0, 1.0; 00838 SymPiso = 1.0/3.0, 1.0/3.0, 1.0/3.0, 0.0, 0.0, 0.0, 00839 1.0/3.0, 1.0/3.0, 1.0/3.0, 0.0, 0.0, 0.0, 00840 1.0/3.0, 1.0/3.0, 1.0/3.0, 0.0, 0.0, 0.0, 00841 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 00842 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 00843 0.0, 0.0, 0.0, 0.0, 0.0, 0.0; 00844 return 0; 00845 } 00846 00848 int __dummy_initialize_tensors = __initialize_tensors (); 00849 00850 }; // namespace MPM 00851 00852 #endif // MPM_TENSORS_H
1.7.6.1
