ug4
Loading...
Searching...
No Matches
math_tensor_functions_common_impl.hpp
Go to the documentation of this file.
1/*
2 * Copyright (c) 2014-2015: G-CSC, Goethe University Frankfurt
3 * Author: Raphael Prohl
4 *
5 * This file is part of UG4.
6 *
7 * UG4 is free software: you can redistribute it and/or modify it under the
8 * terms of the GNU Lesser General Public License version 3 (as published by the
9 * Free Software Foundation) with the following additional attribution
10 * requirements (according to LGPL/GPL v3 §7):
11 *
12 * (1) The following notice must be displayed in the Appropriate Legal Notices
13 * of covered and combined works: "Based on UG4 (www.ug4.org/license)".
14 *
15 * (2) The following notice must be displayed at a prominent place in the
16 * terminal output of covered works: "Based on UG4 (www.ug4.org/license)".
17 *
18 * (3) The following bibliography is recommended for citation and must be
19 * preserved in all covered files:
20 * "Reiter, S., Vogel, A., Heppner, I., Rupp, M., and Wittum, G. A massively
21 * parallel geometric multigrid solver on hierarchically distributed grids.
22 * Computing and visualization in science 16, 4 (2013), 151-164"
23 * "Vogel, A., Reiter, S., Rupp, M., Nägel, A., and Wittum, G. UG4 -- a novel
24 * flexible software system for simulating pde based models on high performance
25 * computers. Computing and visualization in science 16, 4 (2013), 165-179"
26 *
27 * This program is distributed in the hope that it will be useful,
28 * but WITHOUT ANY WARRANTY; without even the implied warranty of
29 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
30 * GNU Lesser General Public License for more details.
31 */
32
33#ifndef __H__UG__COMMON__MATH_TENSOR_FUNCTIONS_COMMON_IMPL__
34#define __H__UG__COMMON__MATH_TENSOR_FUNCTIONS_COMMON_IMPL__
35
36#include <cmath>
37#include <iostream>
38#include <iomanip>
39#include <cassert>
40#include "math_tensor.h"
41#include "common/assert.h"
43
44namespace ug{
45
47// transformation of a tensors
49
50// transformation of a tensor of rank 2 (R^dim x R^dim) into a vector (R^(dim^2))
51template <std::size_t TDim, std::size_t TDimSQ>
52void
54 const MathMatrix<TDim, TDim>& tens2)
55{
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");
60
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];
64 }
65}
66
67// transformation of a vector (R^(dim^2)) into a tensor of rank 2 (R^dim x R^dim)
68template <std::size_t TDim, std::size_t TDimSQ>
69void
72{
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");
77
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];
81 }
82}
83
84// transformation of a tensor of rank 4 (R^dim x R^dim x R^dim x R^dim)
85// into a matrix (R^(dim^2) x R^(dim^2))
86template <std::size_t TDim, std::size_t TDimSQ>
87void
90{
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");
95
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];
101 }
102}
103
104
105// transformation of a matrix (R^(dim^2) x R^(dim^2))
106// into a tensor of rank 4 (R^dim x R^dim x R^dim x R^dim)
107template <std::size_t TDim, std::size_t TDimSQ>
108void
111{
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");
116
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);
122 }
123}
124
125
127// Addition of tensors
129
130// adds a fourth order tensor:
131// tens4_out = tens4a + tens4b
132template <std::size_t TDim>
133void
137{
138 static const size_t dim = TDim;
139
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];
146 }
147}
148
149
151// Subtraction of tensors
153
154// subtracts a fourth order tensor:
155// tens4_out = tens4a - tens4b
156template <std::size_t TDim>
157void
161{
162 static const size_t dim = TDim;
163
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];
170 }
171}
172
173
175// Transpose of a tensor
177
178
179// transposes a fourth order tensor:
180// tens4_out = (tens4)^T
181template <std::size_t TDim>
182void
185{
186 static const size_t dim = TDim;
187
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];
193 }
194}
195
196
198// Invert a rank-4-tensor
200
201// This function inverts a tensor of order 4 by inverting the associated matrix
202// tens4_out = (tens4)^{-1}
203template <std::size_t TDim>
204void
207{
208 static const size_t dim = TDim;
209 static const size_t dimSQ = dim * dim;
210
212
213// get associated matrix from tens4
214 Tens4ToMat(mat, tens4);
215
216// we now create a matrix, where we store the inverse matrix
218
219// get the inverse
220 if(!GetInverse(inv, mat))
221 UG_THROW("Could not compute inverse.");
222
223// get associated tensor of rank 4 from the inverse matrix
224 MatToTens4(tens4_out, inv);
225}
226
227
228/*
229 * This function solves
230 *
231 * A : X = RHS
232 *
233 * X, RHS: tensors of second order
234 * A: tensor of fourth order
235 *
236 * The inversion of the 4th order tensors (R^dim x R^dim x R^dim x R^dim) therein is done by associating
237 * a R^(dim^2) x R^(dim^2) Matrix to the 4th order tensor and by associating a R^(dim^2)
238 * vector to the second order tensors.
239 * Hence, we can solve this 'associated' matrix-vector system (A_mat : x_vec = rhs_vec)
240 */
241template <std::size_t TDim>
242void
245 const MathMatrix<TDim, TDim>& rhs)
246{
247 static const size_t dim = TDim;
248 static const size_t dimSQ = dim * dim;
249
252
253// transform MathTensors of order 4 resp. 2 to DenseMatrices resp. DenseVectors:
254 Tens4ToMat(A_mat,A);
255 Tens2ToVec(rhs_vec,rhs);
256
257 typename block_traits<DenseMatrix< FixedArray2<number, dimSQ, dimSQ> > >::inverse_type A_matI;
258
259// get the inverse
260 if(!GetInverse(A_matI, A_mat))
261 UG_THROW("Could not compute inverse.");
262
263// x_vec = A_matI * rhs_vec
264 MatMult(x_vec, 1.0, A_matI, rhs_vec);
265// get MathMatrix X from DenseVector x:
266 VecToTens2(X,x_vec);
267}
268
269
271// Contractions of rank-4-tensors
273
274// this function computes the contraction of a 4th order tensor by a second order tensor
275// tens2_out = tens4 : tens2
276template <std::size_t TDim>
277void
280 const MathMatrix<TDim, TDim>& tens2)
281{
282 static const size_t dim = TDim;
283
284 for(size_t i = 0; i < dim; ++i)
285 for(size_t j = 0; j < dim; ++j)
286 {
287 tens2_out[i][j] = 0.0;
288
289 for(size_t k = 0; k < dim; ++k)
290 for(size_t l = 0; l < dim; ++l)
291 {
292 tens2_out[i][j] += tens4[i][j][k][l] * tens2[k][l];
293 }
294 }
295}
296
297// this function computes the contraction of a 4th order tensor by a fourth order tensor
298// tens4_out = tens4a : tens4b
299template <std::size_t TDim>
300void
304{
305 static const size_t dim = TDim;
306
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)
311 {
312 tens4_out[i][j][k][l] = 0.0;
313
314 for(size_t m = 0; m < dim; ++m)
315 for(size_t n = 0; n < dim; ++n)
316 {
317 tens4_out[i][j][k][l] += tens4a[i][j][m][n] * tens4b[m][n][k][l];
318 }
319 }
320}
321
322// this function computes
323// tens4_out = tens4a : tens4b : tens4c
324template <std::size_t TDim>
325void
330{
331 static const size_t dim = TDim;
332
334
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)
339 {
340 tens4_out[i][j][k][l] = 0.0;
341
342 for(size_t m = 0; m < dim; ++m)
343 for(size_t n = 0; n < dim; ++n)
344 {
345 help[m][n][k][l] = 0.0;
346
347 for(size_t r = 0; r < dim; ++r)
348 for(size_t s = 0; s < dim; ++s)
349 {
350 help[m][n][k][l] += tens4b[m][n][r][s] * tens4c[r][s][k][l];
351 }
352
353 tens4_out[i][j][k][l] += tens4a[i][j][m][n] * help[m][n][k][l];
354 }
355 }
356}
357
359// Some tensor func
361
362template <std::size_t TDim>
363void
365{
366 static const size_t dim = TDim;
367
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;
373
374}
375
376// this function computes the 4th order identity tensor
377template <std::size_t TDim>
378void
380{
381 static const size_t dim = TDim;
382
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)
387 {
388 Ident[i][j][k][l] = 0.0;
389 if ((i==k) && (j==l)) Ident[i][j][k][l] = 1.0;
390 }
391}
392
393// this function computes the symmetric 4th order identity tensor
394template <std::size_t TDim>
395void
397{
398 static const size_t dim = TDim;
399
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)
404 {
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;
408 }
409}
410
411
412} // end of namespace
413
414#endif /* __H__UG__COMMON__MATH_TENSOR_FUNCTIONS_COMMON_IMPL__ */
parameterString s
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
the ug namespace
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