/* 2008 (c) Dorival M. Pedroso */#ifndef MPM_TENSORS_H#define MPM_TENSORS_H// STL#include <cmath> // for sqrt// Local#include "defs.h"/* Second order tensors: *      y             Sy *      |             | *      |__x      ____|____________ *    ,'        ,'    |          ,'| *   z        ,'     \|/       ,'  | *          ,'  <-----'---   ,'    |  *        ,'   Sxy         ,'      | *      ,'_______________,'   |    | *      |                |    |    | *      |                |    | <------- Sx *      |        __,     |   \|    | *      |        ,'|     |    '   ,' *      |      ,'        |  Sxy ,' *      |    ,'          |   ,' *      |  Sz            | ,' *      |________________|' * * (Symmetric) Tensor components (Mandel's basis) *   / Sx   Sxy  Sxz \ *   | Sxy  Sy   Syz | => { Sx, Sy, Sz, SQ2*Sxy, SQ2*Syz, SQ2*Sxz } *   \ Sxz  Syz  Sz  / * * (Assymmetric) Tensor components *   / Sx   Sxy  Sxz \ *   | Syx  Sy   Syz | => { Sx, Sy, Sz, Sxy, Syz, Sxz, Syx, Szy, Szx } *   \ Szx  Szy  Sz  / *//////////////////////////////////////////////////////////////////////////////////////////// Entities //////** Second order identity tensor (symmetric/Mandel's basis). */STensor2 SymI;/** Forth order identity tensor (symmetric/Mandel's basis). */STensor4 SymII;/** Forth order tensor (symmetric/Mandel's basis). */STensor4 SymIdyI;/** Forth order symmetric-deviatoric tensor (symmetric/Mandel's basis). */STensor4 SymPsd;/** Forth order isotropic tensor (symmetric/Mandel's basis). */STensor4 SymPiso;///////////////////////////////////////////////////////////////////////////////////////// Operations //////** Euclidian norm */inline double Norm (Vector3D const & v){	return sqrt(v(0)+v(1)+v(2));}/** Add scaled tensors Z = a*X + b*Y. */inline void AddScaled (double const & a, STensor4 const & X,                       double const & b, STensor4 const & Y,  STensor4 & Z){	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);	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);	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);	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);	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);	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);}/** Symmetric of a tensor S = a*sym(T) = a * 0.5*(T + trn(T)). */inline void Sym (double a, ATensor2 const & T,  STensor2 & S){	S(0) = a*T(0);	S(1) = a*T(1);	S(2) = a*T(2);	S(3) = a*0.5*(T(3)+T(6))*SQ2;	S(4) = a*0.5*(T(4)+T(7))*SQ2;	S(5) = a*0.5*(T(5)+T(8))*SQ2;}/////////////////////////////////////////////////////////////////////// 3D///////////////////////////////////////// Asymmetric//////////// Only assignment/** 1) Dot product: v=A*u  =>  v(i)=A(i,k)*u(k) */inline void Dot (ATensor2 const & A, Vector3D const & u,  Vector3D & v){	v(0) = A(0)*u(0)+A(3)*u(1)+A(5)*u(2);	v(1) = A(6)*u(0)+A(1)*u(1)+A(4)*u(2);	v(2) = A(8)*u(0)+A(7)*u(1)+A(2)*u(2);} /** 2) Dot product: v=u*A  =>  v(i)=u(k)*A(k,i) */inline void Dot (Vector3D const & u, ATensor2 const & A,  Vector3D & v){	v(0) = u(0)*A(0)+u(1)*A(6)+u(2)*A(8);	v(1) = u(0)*A(3)+u(1)*A(1)+u(2)*A(7);	v(2) = u(0)*A(5)+u(1)*A(4)+u(2)*A(2);}/** 3) Dot product: y = T*x  =>  y = T:x  =>  y(i,j)=T(i,j,k,l)*x(k,l) */inline void Dot (ATensor4 const & T, ATensor2 const & x,  ATensor2 & y){	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);	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);	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);	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);	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);	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);	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);	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);	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);}/** a) Dyadic product: A = u dyad v  =>  A(i,j)=u(i)*v(j) */inline void Dyad (Vector3D const & u, Vector3D const & v,  ATensor2 & A){	A(0) = u(0)*v(0);	A(3) = u(0)*v(1);	A(5) = u(0)*v(2);	A(6) = u(1)*v(0);	A(1) = u(1)*v(1);	A(4) = u(1)*v(2);	A(8) = u(2)*v(0);	A(7) = u(2)*v(1);	A(2) = u(2)*v(2);}/** b) Dyadic with Dot: R=v dyad (F*u)  =>  R(i,j)=v(i)*F(j,k)*u(k) */inline void DyadDot (Vector3D const & v, ATensor2 const & F, Vector3D const & u,  ATensor2 & R){	R(0) = v(0)*F(0)*u(0)+v(0)*F(3)*u(1)+v(0)*F(5)*u(2);	R(3) = v(0)*F(6)*u(0)+v(0)*F(1)*u(1)+v(0)*F(4)*u(2);	R(5) = v(0)*F(8)*u(0)+v(0)*F(7)*u(1)+v(0)*F(2)*u(2);	R(6) = v(1)*F(0)*u(0)+v(1)*F(3)*u(1)+v(1)*F(5)*u(2);	R(1) = v(1)*F(6)*u(0)+v(1)*F(1)*u(1)+v(1)*F(4)*u(2);	R(4) = v(1)*F(8)*u(0)+v(1)*F(7)*u(1)+v(1)*F(2)*u(2);	R(8) = v(2)*F(0)*u(0)+v(2)*F(3)*u(1)+v(2)*F(5)*u(2);	R(7) = v(2)*F(6)*u(0)+v(2)*F(1)*u(1)+v(2)*F(4)*u(2);	R(2) = v(2)*F(8)*u(0)+v(2)*F(7)*u(1)+v(2)*F(2)*u(2);}/** c) Dyadic with Dot: R=(v dyad u)*F  =>  R(i,j)=v(i)*u(k)*F(k,j) */inline void DyadDot (Vector3D const & v, Vector3D const & u, ATensor2 const & F,  ATensor2 & R){	R(0) = v(0)*u(0)*F(0)+v(0)*u(1)*F(6)+v(0)*u(2)*F(8);	R(3) = v(0)*u(0)*F(3)+v(0)*u(1)*F(1)+v(0)*u(2)*F(7);	R(5) = v(0)*u(0)*F(5)+v(0)*u(1)*F(4)+v(0)*u(2)*F(2);	R(6) = v(1)*u(0)*F(0)+v(1)*u(1)*F(6)+v(1)*u(2)*F(8);	R(1) = v(1)*u(0)*F(3)+v(1)*u(1)*F(1)+v(1)*u(2)*F(7);	R(4) = v(1)*u(0)*F(5)+v(1)*u(1)*F(4)+v(1)*u(2)*F(2);	R(8) = v(2)*u(0)*F(0)+v(2)*u(1)*F(6)+v(2)*u(2)*F(8);	R(7) = v(2)*u(0)*F(3)+v(2)*u(1)*F(1)+v(2)*u(2)*F(7);	R(2) = v(2)*u(0)*F(5)+v(2)*u(1)*F(4)+v(2)*u(2)*F(2);}//////////////// With update/** 4) Dot product (with update): v += A*u  =>  v(i)+=A(i,k)*u(k) */inline void DotUp (ATensor2 const & A, Vector3D const & u,  Vector3D & v){	v(0) += A(0)*u(0)+A(3)*u(1)+A(5)*u(2);	v(1) += A(6)*u(0)+A(1)*u(1)+A(4)*u(2);	v(2) += A(8)*u(0)+A(7)*u(1)+A(2)*u(2);} /** 5) Dot product (with update): v += u*A  =>  v(i)+=u(k)*A(k,i) */inline void DotUp (Vector3D const & u, ATensor2 const & A,  Vector3D & v){	v(0) += u(0)*A(0)+u(1)*A(6)+u(2)*A(8);	v(1) += u(0)*A(3)+u(1)*A(1)+u(2)*A(7);	v(2) += u(0)*A(5)+u(1)*A(4)+u(2)*A(2);}/** 6) Dot product (with update): y += T*x  =>  y = T:x  =>  y(i,j)+=T(i,j,k,l)*x(k,l) */inline void DotUp (ATensor4 const & T, ATensor2 const & x,  ATensor2 & y){	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);	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);	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);	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);	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);	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);	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);	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);	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);}/** d) Dyadic product (with update): A += u dyad v  =>  A(i,j)+=u(i)*v(j) */inline void DyadUp (Vector3D const & u, Vector3D const & v,  ATensor2 & A){	A(0) += u(0)*v(0);	A(3) += u(0)*v(1);	A(5) += u(0)*v(2);	A(6) += u(1)*v(0);	A(1) += u(1)*v(1);	A(4) += u(1)*v(2);	A(8) += u(2)*v(0);	A(7) += u(2)*v(1);	A(2) += u(2)*v(2);}/** e) Dyadic with Dot (with update): R += v dyad (F*u)  =>  R(i,j)+=v(i)*F(j,k)*u(k) */inline void DyadDotUp (Vector3D const & v, ATensor2 const & F, Vector3D const & u,  ATensor2 & R){	R(0) += v(0)*F(0)*u(0)+v(0)*F(3)*u(1)+v(0)*F(5)*u(2);	R(3) += v(0)*F(6)*u(0)+v(0)*F(1)*u(1)+v(0)*F(4)*u(2);	R(5) += v(0)*F(8)*u(0)+v(0)*F(7)*u(1)+v(0)*F(2)*u(2);	R(6) += v(1)*F(0)*u(0)+v(1)*F(3)*u(1)+v(1)*F(5)*u(2);	R(1) += v(1)*F(6)*u(0)+v(1)*F(1)*u(1)+v(1)*F(4)*u(2);	R(4) += v(1)*F(8)*u(0)+v(1)*F(7)*u(1)+v(1)*F(2)*u(2);	R(8) += v(2)*F(0)*u(0)+v(2)*F(3)*u(1)+v(2)*F(5)*u(2);	R(7) += v(2)*F(6)*u(0)+v(2)*F(1)*u(1)+v(2)*F(4)*u(2);	R(2) += v(2)*F(8)*u(0)+v(2)*F(7)*u(1)+v(2)*F(2)*u(2);}/** f) Dyadic with Dot (with update): R += (v dyad u)*F  =>  R(i,j)+=v(i)*u(k)*F(k,j) */inline void DyadDotUp (Vector3D const & v, Vector3D const & u, ATensor2 const & F,  ATensor2 & R){	R(0) += v(0)*u(0)*F(0)+v(0)*u(1)*F(6)+v(0)*u(2)*F(8);	R(3) += v(0)*u(0)*F(3)+v(0)*u(1)*F(1)+v(0)*u(2)*F(7);	R(5) += v(0)*u(0)*F(5)+v(0)*u(1)*F(4)+v(0)*u(2)*F(2);	R(6) += v(1)*u(0)*F(0)+v(1)*u(1)*F(6)+v(1)*u(2)*F(8);	R(1) += v(1)*u(0)*F(3)+v(1)*u(1)*F(1)+v(1)*u(2)*F(7);	R(4) += v(1)*u(0)*F(5)+v(1)*u(1)*F(4)+v(1)*u(2)*F(2);	R(8) += v(2)*u(0)*F(0)+v(2)*u(1)*F(6)+v(2)*u(2)*F(8);	R(7) += v(2)*u(0)*F(3)+v(2)*u(1)*F(1)+v(2)*u(2)*F(7);	R(2) += v(2)*u(0)*F(5)+v(2)*u(1)*F(4)+v(2)*u(2)*F(2);}////////////////////////////////////////// Symmetric//////////// Only assignment/** 7) Dot product: v = A*u  =>  v(i)=A(i,k)*u(k) */inline void Dot (STensor2 const & A, Vector3D const & u,  Vector3D & v){	v(0) = A(0)*u(0)     + A(3)*u(1)/SQ2 + A(5)*u(2)/SQ2;	v(1) = A(3)*u(0)/SQ2 + A(1)*u(1)     + A(4)*u(2)/SQ2;	v(2) = A(5)*u(0)/SQ2 + A(4)*u(1)/SQ2 + A(2)*u(2);}/** 8) Dot product: v = u*A  =>  v(i)=u(k)*A(k,i) */inline void Dot (Vector3D const & u, STensor2 const & A,  Vector3D & v){	v(0) = u(0)*A(0)     + u(1)*A(3)/SQ2 + u(2)*A(5)/SQ2;	v(1) = u(0)*A(3)/SQ2 + u(1)*A(1)     + u(2)*A(4)/SQ2;	v(2) = u(0)*A(5)/SQ2 + u(1)*A(4)/SQ2 + u(2)*A(2);}/** 9) Dot product: y = T*x  =>  y = T:x  =>  y(i,j)=T(i,j,k,l)*x(k,l) */inline void Dot (STensor4 const & T, STensor2 const & x,  STensor2 & y){	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);	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);	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);	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);	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);	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);}//////////////// With update/** 10) Dot product (with update): v += A*u  =>  v(i)+=A(i,k)*u(k) */inline void DotUp (STensor2 const & A, Vector3D const & u,  Vector3D & v){	v(0) += A(0)*u(0)     + A(3)*u(1)/SQ2 + A(5)*u(2)/SQ2;	v(1) += A(3)*u(0)/SQ2 + A(1)*u(1)     + A(4)*u(2)/SQ2;	v(2) += A(5)*u(0)/SQ2 + A(4)*u(1)/SQ2 + A(2)*u(2);}/** 11) Dot product (with update): v += u*A  =>  v(i)+=u(k)*A(k,i) */inline void DotUp (Vector3D const & u, STensor2 const & A,  Vector3D & v){	v(0) += u(0)*A(0)     + u(1)*A(3)/SQ2 + u(2)*A(5)/SQ2;	v(1) += u(0)*A(3)/SQ2 + u(1)*A(1)     + u(2)*A(4)/SQ2;	v(2) += u(0)*A(5)/SQ2 + u(1)*A(4)/SQ2 + u(2)*A(2);}/** 11x) Scaled Dot product (with update): v += s*u*A  =>  v(i)+=u(k)*A(k,i) */inline void ScDotUp (double s, Vector3D const & u, STensor2 const & A,  Vector3D & v){	v(0) += s*( u(0)*A(0)     + u(1)*A(3)/SQ2 + u(2)*A(5)/SQ2 );	v(1) += s*( u(0)*A(3)/SQ2 + u(1)*A(1)     + u(2)*A(4)/SQ2 );	v(2) += s*( u(0)*A(5)/SQ2 + u(1)*A(4)/SQ2 + u(2)*A(2)     );}/** 12) Dot product: y += T*x  =>  y += T:x  =>  y(i,j)+=T(i,j,k,l)*x(k,l) */inline void DotUp (STensor4 const & T, STensor2 const & x,  STensor2 & y){	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);	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);	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);	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);	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);	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);}///////////////////////////////////////////////////////////////////////////////////// Initialization //////** Initialize (global) tensors. */int __initialize_tensors (){	SymI    =  1.0, 1.0, 1.0, 0.0, 0.0, 0.0;	SymII   =  1.0, 0.0, 0.0, 0.0, 0.0, 0.0,	           0.0, 1.0, 0.0, 0.0, 0.0, 0.0,	           0.0, 0.0, 1.0, 0.0, 0.0, 0.0,	           0.0, 0.0, 0.0, 1.0, 0.0, 0.0,	           0.0, 0.0, 0.0, 0.0, 1.0, 0.0,	           0.0, 0.0, 0.0, 0.0, 0.0, 1.0;	SymIdyI =  1.0, 1.0, 1.0, 0.0, 0.0, 0.0,	           1.0, 1.0, 1.0, 0.0, 0.0, 0.0,	           1.0, 1.0, 1.0, 0.0, 0.0, 0.0,	           0.0, 0.0, 0.0, 0.0, 0.0, 0.0,	           0.0, 0.0, 0.0, 0.0, 0.0, 0.0,	           0.0, 0.0, 0.0, 0.0, 0.0, 0.0;	SymPsd  =  2.0/3.0, -1.0/3.0, -1.0/3.0, 0.0, 0.0, 0.0,	          -1.0/3.0,  2.0/3.0, -1.0/3.0, 0.0, 0.0, 0.0,	          -1.0/3.0, -1.0/3.0,  2.0/3.0, 0.0, 0.0, 0.0,	               0.0,      0.0,      0.0, 1.0, 0.0, 0.0,	               0.0,      0.0,      0.0, 0.0, 1.0, 0.0,	               0.0,      0.0,      0.0, 0.0, 0.0, 1.0;	SymPiso =  1.0/3.0, 1.0/3.0, 1.0/3.0, 0.0, 0.0, 0.0,	           1.0/3.0, 1.0/3.0, 1.0/3.0, 0.0, 0.0, 0.0,	           1.0/3.0, 1.0/3.0, 1.0/3.0, 0.0, 0.0, 0.0,	               0.0,     0.0,     0.0, 0.0, 0.0, 0.0,	               0.0,     0.0,     0.0, 0.0, 0.0, 0.0,	               0.0,     0.0,     0.0, 0.0, 0.0, 0.0;	return 0;}/** Dummy variable used only for tensors initialization. */int __dummy_initialize_tensors = __initialize_tensors ();#endif // MPM_TENSORS_H/* 2008 (c) Dorival M. Pedroso */