fork download
  1. #include<iostream>
  2. #include<vector>
  3. #include<cmath>
  4. #include<iomanip>
  5.  
  6.  
  7. std::vector<std::vector<int>> mat_voigt_ind(const int N)
  8. {
  9. /** this function will create pointer vector of matrix mat
  10.   * in voigt notation
  11.   */
  12.  
  13. /** for every element A(i,j) of given
  14.   * it will stores it index in the voigt vetor
  15.   * given that matrix is square matrix of size N x N*/
  16. std::vector<std::vector<int>> mat = std::vector<std::vector<int>>(N, std::vector<int>(N, 0)) ;
  17.  
  18.  
  19. // starting and end index of the column to be filled
  20. std::size_t row_s = 0 ;
  21. std::size_t row_e = N - 1 ;
  22.  
  23. std::size_t col_s = 0 ;
  24. std::size_t col_e = N - 1 ;
  25.  
  26. /**we will fill diagonal first
  27.   * then the last column
  28.   * then first row
  29.   *
  30.   * now next iteration
  31.   */
  32.  
  33. // starting index for matrix
  34. int i = row_s ;
  35. int j = col_s ;
  36.  
  37. // k upto (N*N + N)/2 means we will run loop only for upper triangulat matrix
  38. for(std::size_t k = 0; k < (N*N + N)/2; k++)
  39. {
  40. mat[i][j] = k ;
  41.  
  42. /** if non diagonal el, then add the el of lower triangle also */
  43. if(i != j)
  44. mat[j][i] = k + (N*N - N)/2 ;
  45.  
  46. /** now update the index for matrix */
  47.  
  48. // last column
  49. if(j == col_e and i > row_s)
  50. i-- ;
  51.  
  52. // first row
  53. else if(i == row_s and j > col_s + 1)
  54. j-- ;
  55.  
  56. // condition for next cycle
  57. else if(i == row_s and j == col_s + 1)
  58. {
  59. row_s += 1 ;
  60. row_e -= 2 ;
  61.  
  62. col_s += 2 ;
  63. col_e -= 1 ;
  64.  
  65. i = row_s ;
  66. j = col_s ;
  67. }
  68.  
  69. // along the diagonal
  70. else
  71. {
  72. i++ ;
  73. j++ ;
  74. }
  75. }
  76.  
  77. return mat ;
  78. }
  79.  
  80.  
  81. template <class T>
  82. std::vector<std::vector<T>> dyadic(const std::vector<std::vector<T>>& A, const std::vector<std::vector<T>>& B)
  83. {
  84. /** it will retun C as A dyadic B
  85.   * given dimensions of A and B are equal
  86.   */
  87.  
  88. int N = A.size() ;
  89.  
  90. std::vector<std::vector<T>> C(N*N, std::vector<float>(N*N,0)) ; // N is the size of matrix A, i.e. which is calling the function
  91.  
  92. /** if size of A and B are not equal */
  93. if(A.size() != B.size())
  94. return C ;
  95.  
  96. /** generate a matrix of same size as A and B
  97.   * which will store the index of voigt's matrix for that index of matrix
  98.   */
  99. std::vector<std::vector<int>> voigt_mat = mat_voigt_ind(N) ;
  100.  
  101. /** Cijkl = Aij * Bkl */
  102. for(int i= 0; i<N; i++)
  103. for(int j= 0; j<N; j++)
  104. for(int k= 0; k<N; k++)
  105. for(int l= 0; l<N; l++)
  106. C[voigt_mat[i][j]][voigt_mat[k][l]] = A[i][j] * B[k][l] ;
  107. // C ij kl = A ij * B kl
  108.  
  109. return C ;
  110. }
  111.  
  112. template <class T>
  113. std::vector<std::vector<T>> dyadic_up(const std::vector<std::vector<T>>& A, const std::vector<std::vector<T>>& B)
  114. {
  115. /** it will retun C as A dyadic up B
  116.   * given dimensions of A and B are equal
  117.   */
  118.  
  119. int N = A.size() ;
  120.  
  121. std::vector<std::vector<T>> C(N*N, std::vector<float>(N*N,0)) ; // N is the size of matrix A, i.e. which is calling the function
  122.  
  123. /** if size of A and B are not equal */
  124. if(A.size() != B.size())
  125. return C ;
  126.  
  127. /** generate a matrix of same size as A and B
  128.   * which will store the index of voigt's matrix for that index of matrix
  129.   */
  130. std::vector<std::vector<int>> voigt_mat = mat_voigt_ind(N) ;
  131.  
  132. /** Cijkl = Aik * Bjl */
  133. for(int i= 0; i<N; i++)
  134. for(int j= 0; j<N; j++)
  135. for(int k= 0; k<N; k++)
  136. for(int l= 0; l<N; l++)
  137. C[voigt_mat[i][j]][voigt_mat[k][l]] = A[i][k] * B[j][l] ;
  138. // C ij kl = A ik * B jl
  139.  
  140. return C ;
  141. }
  142.  
  143. template <class T>
  144. std::vector<std::vector<T>> dyadic_down(const std::vector<std::vector<T>>& A, const std::vector<std::vector<T>>& B)
  145. {
  146. /** it will retun C as A dyadic down B
  147.   * given dimensions of A and B are equal
  148.   */
  149.  
  150. int N = A.size() ;
  151.  
  152. std::vector<std::vector<T>> C(N*N, std::vector<float>(N*N,0)) ; // N is the size of matrix A, i.e. which is calling the function
  153.  
  154. /** if size of A and B are not equal */
  155. if(A.size() != B.size())
  156. return C ;
  157.  
  158. /** generate a matrix of same size as A and B
  159.   * which will store the index of voigt's matrix for that index of matrix
  160.   */
  161. std::vector<std::vector<int>> voigt_mat = mat_voigt_ind(N) ;
  162.  
  163. /** Cijkl = Aik * Bjl */
  164. for(int i= 0; i<N; i++)
  165. for(int j= 0; j<N; j++)
  166. for(int k= 0; k<N; k++)
  167. for(int l= 0; l<N; l++)
  168. C[voigt_mat[i][j]][voigt_mat[k][l]] = A[i][l] * B[j][k] ;
  169. // C ij kl = A il * B jk
  170.  
  171. return C ;
  172. }
  173.  
  174. template <class T>
  175. void printMatrix(const std::vector<std::vector<T>>& C, int N = 3)
  176. {
  177. for(int i=0; i<C.size(); i++)
  178. {
  179. for(int j= 0; j<C[i].size(); j++)
  180. {
  181. std::cout<<std::setw(10) << std::fixed << std::setprecision(4)<<C[i][j] ;
  182.  
  183. if((j+1)%N == 0)
  184. std::cout<<"\t" ;
  185. }
  186.  
  187.  
  188. if((i+1)%N == 0)
  189. std::cout<<std::endl ;
  190.  
  191. std::cout<<std::endl ;
  192. }
  193. }
  194.  
  195.  
  196. int main()
  197. {
  198. constexpr int N = 3 ;
  199.  
  200. std::vector<std::vector<float>> A(N, std::vector<float>(N,0)) ;
  201. std::vector<std::vector<float>> B(N, std::vector<float>(N,0)) ;
  202.  
  203.  
  204. /** Problem no 5.3
  205.   * considering alpha as 0
  206.   */
  207. float alp = M_PI/2 ;
  208.  
  209. // Rotation Matrix about V1
  210. A = { {1 , 0 , 0}, {0 , cos(alp), sin(alp)}, {0, -sin(alp), cos(alp)}} ;
  211.  
  212. std::cout<<"\nTransformation Matrix for second order tensor for rotation about e1:"<<std::endl ;
  213. printMatrix(dyadic_up(A, A)) ;
  214.  
  215. // Rotation Matrix about V2
  216. B = { {cos(alp), 0, -sin(alp)}, {0 , 1, 0}, {sin(alp), 0, cos(alp)}} ;
  217.  
  218. std::cout<<"\nTransformation Matrix for second order tensor for rotation about e2:"<<std::endl ;
  219. printMatrix(dyadic_up(B, B)) ;
  220.  
  221.  
  222. /** for problem 5.4 */
  223. A = { {5, 3, 0}, {2, 4, 6}, {0, 6, 2}} ;
  224. B = { {4, 5, 2}, {0, 3, 4}, {0, 2, 5}} ;
  225.  
  226. std::cout<<"\nA dyadic B:"<<std::endl ;
  227. printMatrix(dyadic(A, B)) ;
  228.  
  229. std::cout<<"\nA dyadic up B:"<<std::endl ;
  230. printMatrix(dyadic_up(A, B)) ;
  231.  
  232. return 0 ;
  233. }
  234.  
  235.  
Success #stdin #stdout 0s 5316KB
stdin
Standard input is empty
stdout
Transformation Matrix for second order tensor for rotation about e1:
    1.0000    0.0000    0.0000	    0.0000    0.0000    0.0000	    0.0000    0.0000    0.0000	
    0.0000    0.0000    1.0000	   -0.0000    0.0000   -0.0000	   -0.0000    0.0000   -0.0000	
    0.0000    1.0000    0.0000	    0.0000   -0.0000   -0.0000	    0.0000   -0.0000   -0.0000	

    0.0000    0.0000   -0.0000	    0.0000   -0.0000   -0.0000	   -1.0000    0.0000   -0.0000	
    0.0000   -0.0000   -0.0000	   -0.0000   -0.0000   -1.0000	   -0.0000    0.0000    0.0000	
    0.0000   -0.0000    0.0000	    0.0000    1.0000   -0.0000	   -0.0000    0.0000    0.0000	

    0.0000    0.0000   -0.0000	   -1.0000    0.0000   -0.0000	    0.0000   -0.0000   -0.0000	
    0.0000   -0.0000   -0.0000	   -0.0000    0.0000    0.0000	   -0.0000   -0.0000   -1.0000	
    0.0000   -0.0000    0.0000	   -0.0000    0.0000    0.0000	    0.0000    1.0000   -0.0000	


Transformation Matrix for second order tensor for rotation about e2:
    0.0000    0.0000    1.0000	   -0.0000    0.0000   -0.0000	   -0.0000    0.0000   -0.0000	
    0.0000    1.0000    0.0000	    0.0000    0.0000    0.0000	    0.0000    0.0000    0.0000	
    1.0000    0.0000    0.0000	   -0.0000   -0.0000    0.0000	   -0.0000   -0.0000    0.0000	

    0.0000    0.0000   -0.0000	   -0.0000   -0.0000    0.0000	    0.0000    0.0000    1.0000	
   -0.0000    0.0000    0.0000	   -0.0000    0.0000   -0.0000	   -0.0000   -1.0000    0.0000	
   -0.0000    0.0000   -0.0000	    0.0000   -0.0000   -0.0000	   -1.0000   -0.0000    0.0000	

    0.0000    0.0000   -0.0000	    0.0000    0.0000    1.0000	   -0.0000   -0.0000    0.0000	
   -0.0000    0.0000    0.0000	   -0.0000   -1.0000    0.0000	   -0.0000    0.0000   -0.0000	
   -0.0000    0.0000   -0.0000	   -1.0000   -0.0000    0.0000	    0.0000   -0.0000   -0.0000	


A dyadic B:
   20.0000   15.0000   25.0000	   20.0000   10.0000   25.0000	   10.0000    0.0000    0.0000	
   16.0000   12.0000   20.0000	   16.0000    8.0000   20.0000	    8.0000    0.0000    0.0000	
    8.0000    6.0000   10.0000	    8.0000    4.0000   10.0000	    4.0000    0.0000    0.0000	

   24.0000   18.0000   30.0000	   24.0000   12.0000   30.0000	   12.0000    0.0000    0.0000	
    0.0000    0.0000    0.0000	    0.0000    0.0000    0.0000	    0.0000    0.0000    0.0000	
   12.0000    9.0000   15.0000	   12.0000    6.0000   15.0000	    6.0000    0.0000    0.0000	

   24.0000   18.0000   30.0000	   24.0000   12.0000   30.0000	   12.0000    0.0000    0.0000	
    0.0000    0.0000    0.0000	    0.0000    0.0000    0.0000	    0.0000    0.0000    0.0000	
    8.0000    6.0000   10.0000	    8.0000    4.0000   10.0000	    4.0000    0.0000    0.0000	


A dyadic up B:
   20.0000   15.0000    0.0000	    6.0000   10.0000   25.0000	    0.0000    0.0000   12.0000	
    0.0000   12.0000   24.0000	   16.0000    8.0000    6.0000	   18.0000    0.0000    0.0000	
    0.0000   12.0000   10.0000	   30.0000    0.0000    0.0000	    4.0000    0.0000    0.0000	

    0.0000    8.0000   30.0000	   20.0000   10.0000    4.0000	   12.0000    0.0000    0.0000	
    0.0000    6.0000    0.0000	   15.0000   25.0000   10.0000	    0.0000    0.0000    0.0000	
    0.0000    9.0000    0.0000	   12.0000   20.0000   15.0000	    0.0000    0.0000    0.0000	

    0.0000   18.0000    8.0000	   24.0000    0.0000    0.0000	    6.0000    0.0000    0.0000	
    0.0000   30.0000    4.0000	   12.0000    0.0000    0.0000	   10.0000    8.0000   24.0000	
    8.0000   20.0000   12.0000	    8.0000    4.0000   10.0000	   30.0000   24.0000   16.0000