33#ifndef __H__UG__COMMON__MATH_TENSOR_FUNCTIONS_COMMON_IMPL__
34#define __H__UG__COMMON__MATH_TENSOR_FUNCTIONS_COMMON_IMPL__
51template <std::
size_t TDim, std::
size_t TDimSQ>
56 static const size_t dim = TDim;
57 static const size_t dimSQ = TDimSQ;
58 if (dimSQ != dim * dim)
59 UG_THROW(
"Tens2ToVec::Invalid dimensions: " << dimSQ <<
" is not the square of " << dim <<
"! \n");
61 for (
size_t i = 0; i < dim; ++i)
62 for (
size_t j = 0; j < dim; ++j){
63 vec[i*dim + j] = tens2[i][j];
68template <std::
size_t TDim, std::
size_t TDimSQ>
73 static const size_t dim = TDim;
74 static const size_t dimSQ = TDimSQ;
75 if (dimSQ != dim * dim)
76 UG_THROW(
"VecToTens2::Invalid dimensions: " << dimSQ <<
" is not the square of " << dim <<
"! \n");
78 for (
size_t i = 0; i < dim; ++i)
79 for (
size_t j = 0; j < dim; ++j){
80 tens2[i][j] = vec[i*dim + j];
86template <std::
size_t TDim, std::
size_t TDimSQ>
91 static const size_t dim = TDim;
92 static const size_t dimSQ = TDimSQ;
93 if (dimSQ != dim * dim)
94 UG_THROW(
"Tens4ToMat::Invalid dimensions: " << dimSQ <<
" is not the square of " << dim <<
"! \n");
96 for (
size_t i = 0; i < dim; ++i)
97 for (
size_t j = 0; j < dim; ++j)
98 for (
size_t k = 0; k < dim; ++k)
99 for (
size_t l = 0; l < dim; ++l){
100 mat(i*dim + j, k*dim + l) = tens4[i][j][k][l];
107template <std::
size_t TDim, std::
size_t TDimSQ>
112 static const size_t dim = TDim;
113 static const size_t dimSQ = TDimSQ;
114 if (dimSQ != dim * dim)
115 UG_THROW(
"MatToTens4::Invalid dimensions: " << dimSQ <<
" is not the square of " << dim <<
"! \n");
117 for (
size_t i = 0; i < dim; ++i)
118 for (
size_t j = 0; j < dim; ++j)
119 for (
size_t k = 0; k < dim; ++k)
120 for (
size_t l = 0; l < dim; ++l){
121 tens4[i][j][k][l] = mat(i*dim + j, k*dim + l);
132template <std::
size_t TDim>
138 static const size_t dim = TDim;
140 for (
size_t i = 0; i < dim; ++i)
141 for (
size_t j = 0; j < dim; ++j)
142 for (
size_t k = 0; k < dim; ++k)
143 for (
size_t l = 0; l < dim; ++l){
144 tens4_out[i][j][k][l] = tens4a[i][j][k][l]
145 + tens4b[i][j][k][l];
156template <std::
size_t TDim>
162 static const size_t dim = TDim;
164 for (
size_t i = 0; i < dim; ++i)
165 for (
size_t j = 0; j < dim; ++j)
166 for (
size_t k = 0; k < dim; ++k)
167 for (
size_t l = 0; l < dim; ++l){
168 tens4_out[i][j][k][l] = tens4a[i][j][k][l]
169 - tens4b[i][j][k][l];
181template <std::
size_t TDim>
186 static const size_t dim = TDim;
188 for (
size_t i = 0; i < dim; ++i)
189 for (
size_t j = 0; j < dim; ++j)
190 for (
size_t k = 0; k < dim; ++k)
191 for (
size_t l = 0; l < dim; ++l){
192 tens4_out[k][l][i][j] = tens4[i][j][k][l];
203template <std::
size_t TDim>
208 static const size_t dim = TDim;
209 static const size_t dimSQ = dim * dim;
221 UG_THROW(
"Could not compute inverse.");
241template <std::
size_t TDim>
247 static const size_t dim = TDim;
248 static const size_t dimSQ = dim * dim;
261 UG_THROW(
"Could not compute inverse.");
264 MatMult(x_vec, 1.0, A_matI, rhs_vec);
276template <std::
size_t TDim>
282 static const size_t dim = TDim;
284 for(
size_t i = 0; i < dim; ++i)
285 for(
size_t j = 0; j < dim; ++j)
287 tens2_out[i][j] = 0.0;
289 for(
size_t k = 0; k < dim; ++k)
290 for(
size_t l = 0; l < dim; ++l)
292 tens2_out[i][j] += tens4[i][j][k][l] * tens2[k][l];
299template <std::
size_t TDim>
305 static const size_t dim = TDim;
307 for(
size_t i = 0; i < dim; ++i)
308 for(
size_t j = 0; j < dim; ++j)
309 for(
size_t k = 0; k < dim; ++k)
310 for(
size_t l = 0; l < dim; ++l)
312 tens4_out[i][j][k][l] = 0.0;
314 for(
size_t m = 0; m < dim; ++m)
315 for(
size_t n = 0; n < dim; ++n)
317 tens4_out[i][j][k][l] += tens4a[i][j][m][n] * tens4b[m][n][k][l];
324template <std::
size_t TDim>
331 static const size_t dim = TDim;
335 for(
size_t i = 0; i < dim; ++i)
336 for(
size_t j = 0; j < dim; ++j)
337 for(
size_t k = 0; k < dim; ++k)
338 for(
size_t l = 0; l < dim; ++l)
340 tens4_out[i][j][k][l] = 0.0;
342 for(
size_t m = 0; m < dim; ++m)
343 for(
size_t n = 0; n < dim; ++n)
345 help[m][n][k][l] = 0.0;
347 for(
size_t r = 0; r < dim; ++r)
348 for(
size_t s = 0;
s < dim; ++
s)
350 help[m][n][k][l] += tens4b[m][n][r][
s] * tens4c[r][
s][k][l];
353 tens4_out[i][j][k][l] += tens4a[i][j][m][n] * help[m][n][k][l];
362template <std::
size_t TDim>
366 static const size_t dim = TDim;
368 for(
size_t i = 0; i < dim; ++i)
369 for(
size_t j = 0; j < dim; ++j)
370 for(
size_t k = 0; k < dim; ++k)
371 for(
size_t l = 0; l < dim; ++l)
372 tensOut[i][j][k][l] = 0.0;
377template <std::
size_t TDim>
381 static const size_t dim = TDim;
383 for(
size_t i = 0; i < dim; ++i)
384 for(
size_t j = 0; j < dim; ++j)
385 for(
size_t k = 0; k < dim; ++k)
386 for(
size_t l = 0; l < dim; ++l)
388 Ident[i][j][k][l] = 0.0;
389 if ((i==k) && (j==l)) Ident[i][j][k][l] = 1.0;
394template <std::
size_t TDim>
398 static const size_t dim = TDim;
400 for(
size_t i = 0; i < dim; ++i)
401 for(
size_t j = 0; j < dim; ++j)
402 for(
size_t k = 0; k < dim; ++k)
403 for(
size_t l = 0; l < dim; ++l)
405 Ident[i][j][k][l] = 0.0;
406 if ((i==k) && (j==l)) Ident[i][j][k][l] += 0.5;
407 if ((i==l) && (j==k)) Ident[i][j][k][l] += 0.5;
Definition densematrix.h:57
Definition lapack_densematrix_inverse.h:51
Definition densevector.h:101
Definition fixed_array.h:56
Definition fixed_array.h:135
A class for fixed size, dense matrices.
Definition math_matrix.h:63
Definition math_tensor.h:208
void MatToTens4(MathTensor4< TDim, TDim, TDim, TDim > &tens4, const DenseMatrixInverse< FixedArray2< number, TDimSQ, TDimSQ > > &mat)
Definition math_tensor_functions_common_impl.hpp:109
void Tens4Contract(MathMatrix< TDim, TDim > &tens2_out, const MathTensor4< TDim, TDim, TDim, TDim > &tens4, const MathMatrix< TDim, TDim > &tens2)
Definition math_tensor_functions_common_impl.hpp:278
void VecToTens2(MathMatrix< TDim, TDim > &tens2, const DenseVector< FixedArray1< number, TDimSQ > > &vec)
transformation of a vector (R^(dim^2)) into a tensor of rank 2 (R^dim x R^dim)
Definition math_tensor_functions_common_impl.hpp:70
void Tens4Zero(MathTensor4< TDim, TDim, TDim, TDim > &tensOut)
Definition math_tensor_functions_common_impl.hpp:364
void Tens4Identity(MathTensor4< TDim, TDim, TDim, TDim > &Ident)
this function computes the 4th order identity tensor
Definition math_tensor_functions_common_impl.hpp:379
void TransTens4(MathTensor4< TDim, TDim, TDim, TDim > &tens4_out, const MathTensor4< TDim, TDim, TDim, TDim > &tens4)
Definition math_tensor_functions_common_impl.hpp:183
void InvertTensor4(MathTensor4< TDim, TDim, TDim, TDim > &tens4_out, const MathTensor4< TDim, TDim, TDim, TDim > &tens4)
Definition math_tensor_functions_common_impl.hpp:205
void SolveTensorMatrixEquation(MathMatrix< TDim, TDim > &X, const MathTensor4< TDim, TDim, TDim, TDim > &A, const MathMatrix< TDim, TDim > &rhs)
Definition math_tensor_functions_common_impl.hpp:243
void Tens4Add(MathTensor4< TDim, TDim, TDim, TDim > &tens4_out, const MathTensor4< TDim, TDim, TDim, TDim > &tens4a, const MathTensor4< TDim, TDim, TDim, TDim > &tens4b)
Definition math_tensor_functions_common_impl.hpp:134
void Tens4IdentitySym(MathTensor4< TDim, TDim, TDim, TDim > &Ident)
this function computes the symmetric 4th order identity tensor
Definition math_tensor_functions_common_impl.hpp:396
void Tens4Subtract(MathTensor4< TDim, TDim, TDim, TDim > &tens4_out, const MathTensor4< TDim, TDim, TDim, TDim > &tens4a, const MathTensor4< TDim, TDim, TDim, TDim > &tens4b)
Definition math_tensor_functions_common_impl.hpp:158
void Tens4ToMat(DenseMatrix< FixedArray2< number, TDimSQ, TDimSQ > > &mat, const MathTensor4< TDim, TDim, TDim, TDim > &tens4)
Definition math_tensor_functions_common_impl.hpp:88
void Tens2ToVec(DenseVector< FixedArray1< number, TDimSQ > > &vec, const MathMatrix< TDim, TDim > &tens2)
transformation of a tensor of rank 2 (R^dim x R^dim) into a vector (R^(dim^2))
Definition math_tensor_functions_common_impl.hpp:53
#define UG_THROW(msg)
Definition error.h:57
bool GetInverse(block_traits< T >::inverse_type &inv, const T &m)
bool MatMult(vector_t &dest, const number &beta1, const matrix_t &A1, const vector_t &w1)
calculates dest = beta1 * A1;
Definition operations_mat.h:59
Definition communication_policies.h:58