tensors.h

crack modeling with xfem

/* 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 */