1 #ifndef VOIGT_OPERATIONS_H 2 #define VOIGT_OPERATIONS_H 3 #include "MatrixDef.h" 4 #include "MatrixLibrary.h" 5 6 7 8 // Voigt indexing puts a symmetric 3x3 matrix into a 9 // vector form: [0 1 2 3 4 5] 10 // 11 // matrix form: [[ 0 5 4 ] 12 // [ 5 1 3 ] 13 // [ 4 3 2 ]] 14 // 15 // unsymmetric version 16 // vector form: [0 1 2 3 4 5 6 7 8] 17 // 18 // matrix form: [[ 0 5 4 ] 19 // [ 8 1 3 ] 20 // [ 7 6 2 ]] 21 22 namespace voigt3 { 23 24 25 //const int voigt_idx1_symm[] = {0,1,2,1,0,0}; // first packed voigt index 26 //const int voigt_idx2_symm[] = {0,1,2,2,2,1}; // second packed voigt index 27 const int voigt_idx1[] = {0,1,2,1,0,0,2,2,1}; // first packed voigt index 28 const int voigt_idx2[] = {0,1,2,2,2,1,1,0,0}; // second packed voigt index 29 // const int voigt_idx1[] = {0,1,2,0,0,1,1,2,2}; // first packed voigt index 30 // const int voigt_idx2[] = {0,1,2,1,2,2,0,0,1}; // second packed voigt index 31 32 33 //* Computes a symmetric matrix-matrix product 34 //* Inputs 6-length vectors A, B dsymm(const DENS_VEC & A,const DENS_VEC & B)35 inline DENS_VEC dsymm(const DENS_VEC &A, const DENS_VEC &B) 36 { 37 DENS_VEC C(6,false); 38 C(0) = A(0)*B(0)+A(5)*B(5)+A(4)*B(4); 39 C(1) = A(5)*B(5)+A(1)*B(1)+A(3)*B(3); 40 C(2) = A(4)*B(4)+A(3)*B(3)+A(2)*B(2); 41 C(3) = A(5)*B(4)+A(1)*B(3)+A(3)*B(2); 42 C(4) = A(0)*B(4)+A(5)*B(3)+A(4)*B(2); 43 C(5) = A(0)*B(5)+A(5)*B(1)+A(4)*B(3); 44 return C; 45 } 46 47 //* Returns the trace of a 3x3 matrix in symmetric voigt form. tr(const DENS_VEC & A)48 inline double tr(const DENS_VEC &A) 49 { 50 return A(0) + A(1) + A(2); 51 } 52 53 //* Computes the determinant of a 3x3 matrix in symmetric voigt form. det(const DENS_VEC & A)54 inline double det(const DENS_VEC &A) 55 { 56 return A(0) * (A(1)*A(2)-A(3)*A(3)) 57 -A(5) * (A(5)*A(2)-A(3)*A(4)) 58 +A(4) * (A(5)*A(3)-A(1)*A(4)); 59 } 60 61 //* Returns the derivative of C*C in voigt notation. derivative_of_square(const DENS_VEC & C)62 inline DENS_MAT derivative_of_square(const DENS_VEC &C) 63 { 64 DENS_MAT D(6,6); 65 D(0,0)=2.0*C(0); D(0,1)=0.0; D(0,2)=0.0; 66 D(1,0)=0.0; D(1,1)=2.0*C(1); D(1,2)=0.0; 67 D(2,0)=0.0; D(2,1)=0.0; D(2,2)=2.0*C(2); 68 69 D(0,3)=0.0; D(0,4)=2.0*C(4); D(0,5)=2.0*C(5); 70 D(1,3)=2.0*C(3); D(1,4)=0.0; D(1,5)=2.0*C(5); 71 D(2,3)=2.0*C(3); D(2,4)=2.0*C(4); D(2,5)=0.0; 72 73 D(3,0)=0.0; D(3,1)=C(3); D(3,2)=C(3); 74 D(4,0)=C(4); D(4,1)=0.0; D(4,2)=C(4); 75 D(5,0)=C(5); D(5,1)=C(5); D(5,2)=0.0; 76 77 D(3,3)=C(1)+C(2); D(3,4)=C(5); D(3,5)=C(4); 78 D(4,3)=C(5); D(4,4)=C(0)+C(2); D(4,5)=C(3); 79 D(5,3)=C(4); D(5,4)=C(3); D(5,5)=C(0)+C(1); 80 return D; 81 } 82 83 //* Computes the inverse of a 3x3 matrix in symmetric voigt form. inv(const DENS_VEC & A)84 inline DENS_VEC inv(const DENS_VEC &A) 85 { 86 DENS_VEC B(6,false); 87 const double inv_det = 1.0/det(A); 88 B(0) = (A(1)*A(2)-A(3)*A(3))*inv_det; 89 B(1) = (A(0)*A(2)-A(4)*A(4))*inv_det; 90 B(2) = (A(0)*A(1)-A(5)*A(5))*inv_det; 91 B(3) = (A(4)*A(5)-A(0)*A(3))*inv_det; 92 B(4) = (A(5)*A(3)-A(4)*A(1))*inv_det; 93 B(5) = (A(4)*A(3)-A(5)*A(2))*inv_det; 94 return B; 95 } 96 97 //* Returns the identify matrix in voigt form, optionally scaled by a factor. 98 inline DENS_VEC eye(INDEX N=3, double scale=1.0) 99 { 100 const double dij[] = {0.0, scale}; 101 const INDEX voigt_size = N*N-((N*N-N)>>1); // total - symmetric elements 102 DENS_VEC I(voigt_size,false); 103 for (INDEX i=0; i<voigt_size; i++) I(i) = dij[i<N]; 104 return I; 105 } 106 107 108 //* Returns the voigt form of a symmetric matrix. 109 // consistent with voigt_idx1,2 to_voigt(const DENS_MAT & C)110 inline DENS_VEC to_voigt(const DENS_MAT &C) 111 { 112 DENS_VEC B(6,false); 113 B(0)=C(0,0); 114 B(1)=C(1,1); 115 B(2)=C(2,2); 116 B(3)=C(1,2); // take upper triangle entries 117 B(4)=C(0,2); 118 B(5)=C(0,1); 119 return B; 120 } 121 122 //* Returns a symmetric matrix form a voigt form. 123 // consistent with voigt_idx1,2 from_voigt(const DENS_VEC & B)124 inline DENS_MAT from_voigt(const DENS_VEC &B) 125 { 126 DENS_MAT C(3,3,false); 127 C(0,0)=B(0); C(0,1)=B(5); C(0,2)=B(4); 128 C(1,0)=B(5); C(1,1)=B(1); C(1,2)=B(3); 129 C(2,0)=B(4); C(2,1)=B(3); C(2,2)=B(2); 130 return C; 131 } 132 133 //* Returns the voigt form of an unsymmetric matrix. 134 // consistent with voigt_idx1,2 to_voigt_unsymmetric(const DENS_MAT & C)135 inline DENS_VEC to_voigt_unsymmetric(const DENS_MAT &C) 136 { 137 DENS_VEC B(9,false); 138 B(0)=C(0,0); 139 B(1)=C(1,1); 140 B(2)=C(2,2); 141 B(3)=C(1,2); // upper triangle entries 142 B(4)=C(0,2); 143 B(5)=C(0,1); 144 B(6)=C(2,1); // lower triangle entries 145 B(7)=C(2,0); 146 B(8)=C(1,0); 147 return B; 148 } 149 150 //* Returns a symmetric matrix form a voigt form. 151 // consistent with voigt_idx1,2 from_voigt_unsymmetric(const DENS_VEC & B)152 inline DENS_MAT from_voigt_unsymmetric(const DENS_VEC &B) 153 { 154 DENS_MAT C(3,3,false); 155 C(0,0)=B(0); C(0,1)=B(5); C(0,2)=B(4); 156 C(1,0)=B(8); C(1,1)=B(1); C(1,2)=B(3); 157 C(2,0)=B(7); C(2,1)=B(6); C(2,2)=B(2); 158 return C; 159 } 160 161 //* adds the identity to an unsymmetric matrix form add_identity_voigt_unsymmetric(DENS_VEC & B)162 inline void add_identity_voigt_unsymmetric(DENS_VEC &B) 163 { 164 B(0) +=1; 165 B(1) +=1; 166 B(2) +=1; 167 } 168 169 //* Converts voigt vector form to 3x3 matrix for non-symmetric tensor at specified node. vector_to_matrix(const int i,const DENS_MAT & IN,DENS_MAT & OUT)170 inline void vector_to_matrix(const int i, const DENS_MAT & IN, DENS_MAT & OUT) 171 { 172 OUT.reset(3,3); 173 OUT(0,0)=IN(i,0); OUT(0,1)=IN(i,1); OUT(0,2)=IN(i,2); 174 OUT(1,0)=IN(i,3); OUT(1,1)=IN(i,4); OUT(1,2)=IN(i,5); 175 OUT(2,0)=IN(i,6); OUT(2,1)=IN(i,7); OUT(2,2)=IN(i,8); 176 return; 177 } 178 179 //* Converts 3x3 matrix to voigt vector form for non-symmetric tensor at specified node. matrix_to_vector(const int i,const DENS_MAT & IN,DENS_MAT & OUT)180 inline void matrix_to_vector(const int i, const DENS_MAT & IN, DENS_MAT & OUT) 181 { 182 OUT(i,0) = IN(0,0); 183 OUT(i,1) = IN(0,1); 184 OUT(i,2) = IN(0,2); 185 OUT(i,3) = IN(1,0); 186 OUT(i,4) = IN(1,1); 187 OUT(i,5) = IN(1,2); 188 OUT(i,6) = IN(2,0); 189 OUT(i,7) = IN(2,1); 190 OUT(i,8) = IN(2,2); 191 return; 192 } 193 194 //* Converts voigt vector form to 3x3 matrix for symmetric tensor at specified node. vector_to_symm_matrix(const int i,const DENS_MAT & IN,DENS_MAT & OUT)195 inline void vector_to_symm_matrix(const int i, const DENS_MAT & IN, DENS_MAT & OUT) 196 { 197 OUT.reset(3,3); 198 OUT(0,0)=IN(i,0); OUT(0,1)=IN(i,5); OUT(0,2)=IN(i,4); 199 OUT(1,0)=IN(i,5); OUT(1,1)=IN(i,1); OUT(1,2)=IN(i,3); 200 OUT(2,0)=IN(i,4); OUT(2,1)=IN(i,3); OUT(2,2)=IN(i,2); 201 return; 202 } 203 204 //* Converts 3x3 matrix to voigt vector form for symmetric tensor at specified node. symm_matrix_to_vector(const int i,const DENS_MAT & IN,DENS_MAT & OUT)205 inline void symm_matrix_to_vector(const int i, const DENS_MAT & IN, DENS_MAT & OUT) 206 { 207 OUT(i,0) = IN(0,0); 208 OUT(i,1) = IN(1,1); 209 OUT(i,2) = IN(1,2); 210 OUT(i,3) = IN(1,2); 211 OUT(i,4) = IN(0,2); 212 OUT(i,5) = IN(0,1); 213 return; 214 } 215 216 //* Converts voigt vector form to vector at specified node. global_vector_to_vector(const int i,const DENS_MAT & IN)217 inline DENS_VEC global_vector_to_vector(const int i, const DENS_MAT & IN) 218 { 219 DENS_VEC OUT(9); 220 OUT(0)=IN(i,0); OUT(5)=IN(i,1); OUT(4)=IN(i,2); 221 OUT(8)=IN(i,3); OUT(1)=IN(i,4); OUT(3)=IN(i,5); 222 OUT(7)=IN(i,6); OUT(6)=IN(i,7); OUT(2)=IN(i,8); 223 return OUT; 224 } vector_to_global_vector(const int i,const DENS_VEC & IN,DENS_MAT & OUT)225 inline void vector_to_global_vector(const int i, const DENS_VEC & IN, DENS_MAT & OUT) 226 { 227 OUT(i,0) = IN(0); 228 OUT(i,1) = IN(5); 229 OUT(i,2) = IN(4); 230 OUT(i,3) = IN(8); 231 OUT(i,4) = IN(1); 232 OUT(i,5) = IN(3); 233 OUT(i,6) = IN(7); 234 OUT(i,7) = IN(6); 235 OUT(i,8) = IN(2); 236 return; 237 } 238 239 //* Converts vector to DENS_MAT_VEC vector_to_dens_mat_vec(const DENS_MAT & IN,DENS_MAT_VEC & OUT)240 inline void vector_to_dens_mat_vec(const DENS_MAT & IN, DENS_MAT_VEC & OUT) 241 { 242 for (int i=0; i<IN.nRows(); i++) { 243 for (int j=0; j<3; j++) { 244 for (DENS_MAT_VEC::size_type k=0; k<3; k++) { 245 OUT[k](i,j) = IN(i,3*j+k); 246 } 247 } 248 } 249 return; 250 } 251 252 //* Converts DENS_MAT_VEC to vector symm_dens_mat_vec_to_vector(const DENS_MAT_VEC & IN,DENS_MAT & OUT)253 inline void symm_dens_mat_vec_to_vector(const DENS_MAT_VEC & IN, DENS_MAT & OUT) 254 { 255 for (int i=0; i<IN.front().nRows(); i++) { 256 for (int v=0; v<6; v++) { 257 OUT(i,v) = IN[voigt_idx1[v]](i,voigt_idx2[v]); 258 } 259 } 260 return; 261 } 262 263 } 264 #endif 265 266