MarmotTensorExponential.h
Go to the documentation of this file.
1 /* ---------------------------------------------------------------------
2  * _
3  * _ __ ___ __ _ _ __ _ __ ___ ___ | |_
4  * | '_ ` _ \ / _` | '__| '_ ` _ \ / _ \| __|
5  * | | | | | | (_| | | | | | | | | (_) | |_
6  * |_| |_| |_|\__,_|_| |_| |_| |_|\___/ \__|
7  *
8  * Unit of Strength of Materials and Structural Analysis
9  * University of Innsbruck,
10  * 2020 - today
11  *
12  * festigkeitslehre@uibk.ac.at
13  *
14  * Matthias Neuner matthias.neuner@uibk.ac.at
15  *
16  * This file is part of the MAteRialMOdellingToolbox (marmot).
17  *
18  * This library is free software; you can redistribute it and/or
19  * modify it under the terms of the GNU Lesser General Public
20  * License as published by the Free Software Foundation; either
21  * version 2.1 of the License, or (at your option) any later version.
22  *
23  * The full text of the license can be found in the file LICENSE.md at
24  * the top level directory of marmot.
25  * ---------------------------------------------------------------------
26  */
27 #pragma once
28 #include "Fastor/Fastor.h"
29 #include "Marmot/MarmotMath.h"
30 #include <iostream>
31 
32 namespace Marmot {
33  namespace ContinuumMechanics::TensorUtility {
34  namespace TensorExponential {
35  struct ExponentialMapFailed : std::exception {};
36 
37  template < typename T, size_t tensorSize >
38  Fastor::Tensor< T, tensorSize, tensorSize > computeTensorExponential(
39  const Fastor::Tensor< T, tensorSize, tensorSize >& theTensor,
40  int maxIterations,
41  double tolerance,
42  double alternativeTolerance )
43  {
44  using TensorNN = Fastor::Tensor< T, tensorSize, tensorSize >;
45 
46  TensorNN theExponential;
47 
48  int n;
49  int facN = 1;
50 
51  // term 1 ( n = 0 )
52  theExponential.eye();
53 
54  // series 2 ( n = 1 )
55  theExponential += theTensor;
56 
57  TensorNN tensorPower = theTensor;
58  // loop starts at n = 2
59  //
60  double norm = 0.0;
61  for ( n = 2; n < maxIterations; ++n ) {
62  facN *= n;
63  tensorPower = tensorPower % theTensor;
64  const TensorNN series = 1. / facN * tensorPower;
65  theExponential += series;
66 
67  // workaround for autodiff::dual
68  Fastor::Tensor< double, tensorSize, tensorSize > series_real;
69  for ( size_t i = 0; i < tensorSize; i++ ) {
70  for ( size_t j = 0; j < tensorSize; j++ ) {
71  series_real( i, j ) = Math::makeReal( series( i, j ) );
72  }
73  }
74  norm = Fastor::norm( series_real );
75 
76  if ( norm <= tolerance )
77  break;
78  }
79 
80  if ( n == maxIterations && norm > alternativeTolerance ) {
81  throw ExponentialMapFailed();
82  }
83 
84  return theExponential;
85  }
86 
87  template < typename T, std::size_t tensorSize >
89  Fastor::Tensor< T, tensorSize, tensorSize > theExponential;
90  Fastor::Tensor< T, tensorSize, tensorSize, tensorSize, tensorSize > theExponentialDerivative;
91  };
92 
93  namespace FirstOrderDerived {
94  template < typename T, size_t tensorSize >
96  const Fastor::Tensor< T, tensorSize, tensorSize >& theTensor,
97  int maxIterations,
98  double tolerance )
99  {
100  // tensor exponential;
101  using Tensor = Fastor::Tensor< T, tensorSize, tensorSize >;
103 
104  std::vector< Tensor > tensorPowers;
105 
106  int n;
107  int facN = 1;
108 
109  // term 1 ( n = 0 )
110  theResult.theExponential.eye();
111  tensorPowers.push_back( theResult.theExponential ); // I
112 
113  // series 2 ( n = 1 )
114  theResult.theExponential += theTensor;
115  tensorPowers.push_back( theTensor ); // T
116 
117  // loop starts at n = 2
118  for ( n = 2; n < maxIterations; ++n ) {
119  facN *= n;
120  const Tensor tensorPower = tensorPowers.back() % theTensor;
121  const Tensor series = 1. / facN * tensorPower;
122  theResult.theExponential += series;
123 
124  tensorPowers.push_back( tensorPower );
125 
126  if ( Fastor::norm( series ) <= tolerance )
127  break;
128  }
129 
130  // derivative
131  const int N = n;
132 
133  facN = 1;
134  Fastor::Tensor< T, tensorSize, tensorSize, tensorSize, tensorSize > innerSeries;
135  Fastor::Tensor< T, tensorSize, tensorSize, tensorSize, tensorSize > outerSeries;
136  outerSeries.zeros();
137  for ( n = 1; n <= N; ++n ) {
138  innerSeries.zeros();
139  for ( int m = 1; m <= n; ++m ) {
140  innerSeries += Fastor::outer( tensorPowers[m - 1], tensorPowers[n - m] ); // ik_lj
141  }
142  facN *= n;
143  outerSeries += 1. / facN * innerSeries;
144  }
145  theResult.theExponentialDerivative = Fastor::permute< Fastor::Index< 0, 3, 1, 2 > >(
146  outerSeries ); // ik_lj -> ij_kl
147 
148  return theResult;
149  }
150  } // namespace FirstOrderDerived
151  } // namespace TensorExponential
152  } // namespace ContinuumMechanics::TensorUtility
153 } // namespace Marmot
Marmot::FastorIndices::i
Fastor::Index< i_ > i
Definition: MarmotFastorTensorBasics.h:137
Marmot::ContinuumMechanics::TensorUtility::TensorExponential::TensorExponentialResult::theExponentialDerivative
Fastor::Tensor< T, tensorSize, tensorSize, tensorSize, tensorSize > theExponentialDerivative
Definition: MarmotTensorExponential.h:90
MarmotMath.h
Marmot::ContinuumMechanics::TensorUtility::TensorExponential::computeTensorExponential
Fastor::Tensor< T, tensorSize, tensorSize > computeTensorExponential(const Fastor::Tensor< T, tensorSize, tensorSize > &theTensor, int maxIterations, double tolerance, double alternativeTolerance)
Definition: MarmotTensorExponential.h:38
Marmot::Math::makeReal
double makeReal(const double &value)
Definition: MarmotMath.cpp:37
Marmot::FastorIndices::j
Fastor::Index< j_ > j
Definition: MarmotFastorTensorBasics.h:182
Marmot
This file includes functions needed for calculations with stress and strain tensors written in voigt ...
Definition: MarmotTesting.h:37
Marmot::ContinuumMechanics::TensorUtility::TensorExponential::ExponentialMapFailed
Definition: MarmotTensorExponential.h:35
Marmot::ContinuumMechanics::TensorUtility::TensorExponential::FirstOrderDerived::computeTensorExponential
TensorExponentialResult< T, tensorSize > computeTensorExponential(const Fastor::Tensor< T, tensorSize, tensorSize > &theTensor, int maxIterations, double tolerance)
Definition: MarmotTensorExponential.h:95
Marmot::FastorIndices::m
Fastor::Index< m_ > m
Definition: MarmotFastorTensorBasics.h:210
Marmot::ContinuumMechanics::TensorUtility::TensorExponential::TensorExponentialResult
Definition: MarmotTensorExponential.h:88
Marmot::FiniteElement::Spatial1D::Bar2::N
NSized N(double xi)
Definition: MarmotFiniteElement1D.cpp:15
Marmot::ContinuumMechanics::TensorUtility::TensorExponential::TensorExponentialResult::theExponential
Fastor::Tensor< T, tensorSize, tensorSize > theExponential
Definition: MarmotTensorExponential.h:89