MarmotAutomaticDifferentiationForFastor.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  * Alexander Dummer alexander.dummer@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 
28 #pragma once
29 #include "Fastor/Fastor.h"
30 #include "Marmot/MarmotAutomaticDifferentiation.h"
31 #include "Marmot/MarmotFastorTensorBasics.h"
32 #include "autodiff/forward/dual.hpp"
33 #include "autodiff/forward/dual/eigen.hpp"
34 #include <autodiff/forward/dual/dual.hpp>
35 #include <functional>
36 
37 namespace Marmot {
38 
39  namespace AutomaticDifferentiation {
40 
41  using namespace autodiff;
42 
43  template < size_t order, size_t... Rest >
44  Fastor::Tensor< HigherOrderDual< order + 1, double >, Rest... > increaseDualOrderWithShift(
45  const Fastor::Tensor< HigherOrderDual< order, double >, Rest... >& in )
46  {
47  using in_scalar_type = HigherOrderDual< order, double >;
48  using out_scalar_type = HigherOrderDual< order + 1, double >;
49 
50  Fastor::Tensor< out_scalar_type, Rest... > out( 0.0 );
51  out_scalar_type* out_data = out.data();
52  in_scalar_type* in_data = in.data();
53 
54  for ( Fastor::FASTOR_INDEX i = 0; i < in.size(); ++i ) {
55  out_data[out.get_mem_index( i )] = increaseDualOrderWithShift< order >( in_data[in.get_mem_index( i )] );
56  }
57  return out;
58  }
59 
60  template < size_t... Rest >
61  using tensor_to_scalar_function_type = std::function< dual( const Fastor::Tensor< dual, Rest... >& T ) >;
62 
63  template < size_t... Rest >
64  Fastor::Tensor< double, Rest... > df_dT( const tensor_to_scalar_function_type< Rest... >& f,
65  const Fastor::Tensor< double, Rest... >& T )
66  {
67  Fastor::Tensor< double, Rest... > df_dT( 0.0 );
68  Fastor::Tensor< dual, Rest... > T_right = makeDual( T );
69 
70  double* df_dT_data = df_dT.data();
71  dual* T_right_data = T_right.data();
72 
73  for ( Fastor::FASTOR_INDEX i = 0; i < T.size(); ++i ) {
74  const int T_right_mem_idx = T_right.get_mem_index( i );
75  seed< 1 >( T_right_data[T_right_mem_idx], 1.0 );
76  df_dT_data[df_dT.get_mem_index( i )] = derivative< 1 >( f( T_right ) );
77  seed< 1 >( T_right_data[T_right_mem_idx], 0.0 );
78  }
79 
80  return df_dT;
81  }
82 
83  template < size_t order, size_t... Rest >
84  using tensor_to_scalar_function_type_arbitrary_dual_order = std::function< HigherOrderDual< order, double >(
85  const Fastor::Tensor< HigherOrderDual< order, double >, Rest... >& T ) >;
86 
87  template < size_t order, size_t... Rest >
88  std::pair< HigherOrderDual< order, double >, Fastor::Tensor< HigherOrderDual< order, double >, Rest... > > df_dT(
89  const tensor_to_scalar_function_type_arbitrary_dual_order< order + 1, Rest... >& f,
90  const Fastor::Tensor< HigherOrderDual< order, double >, Rest... >& T )
91  {
92 
93  using scalartype = HigherOrderDual< order, double >;
94  using higherOrderScalartype = HigherOrderDual< order + 1, double >;
95 
96  higherOrderScalartype f_;
97  Fastor::Tensor< scalartype, Rest... > df_dT_( 0.0 );
98  Fastor::Tensor< higherOrderScalartype, Rest... > T_right = increaseDualOrderWithShift< order >( T );
99 
100  higherOrderScalartype* T_right_data = T_right.data();
101  scalartype* df_dT_data = df_dT_.data();
102 
103  for ( Fastor::FASTOR_INDEX i = 0; i < T.size(); ++i ) {
104  seed< 1 >( T_right_data[T_right.get_mem_index( i )], 1.0 );
105  f_ = f( T_right );
106  df_dT_data[df_dT_.get_mem_index( i )] = decreaseDualOrderWithShift< order + 1 >( f_ );
107  seed< 1 >( T_right_data[T_right.get_mem_index( i )], 0.0 );
108  }
109  f_ = f( T_right );
110  return { decreaseDualOrder< order + 1 >( f_ ), df_dT_ };
111  }
112 
113  template < size_t... RestF, size_t... RestT >
114  std::pair< Fastor::Tensor< double, RestF... >, Fastor::Tensor< double, RestF..., RestT... > > dF_dT(
115  std::function< Fastor::Tensor< dual, RestF... >( const Fastor::Tensor< dual, RestT... >& ) >& F,
116  const Fastor::Tensor< double, RestT... >& T )
117  {
118 
119  Fastor::Tensor< double, RestF... > F_( 0.0 );
120  Fastor::Tensor< double, RestF..., RestT... > dF_dT_( 0.0 );
121  Fastor::Tensor< dual, RestT... > T_right = makeDual( T );
122 
123  Fastor::Tensor< dual, RestF... > F_at_T_right( 0.0 );
124 
125  double* dF_dT_data = dF_dT_.data();
126  dual* T_right_data = T_right.data();
127  dual* F_at_T_right_data = F_at_T_right.data();
128 
129  for ( Fastor::FASTOR_INDEX i = 0; i < T.size(); ++i ) {
130  const int T_right_mem_idx = T_right.get_mem_index( i );
131  T_right_data[T_right_mem_idx].grad += 1.0;
132  F_at_T_right = F( T_right );
133 
134  for ( Fastor::FASTOR_INDEX j = 0; j < F_at_T_right.size(); ++j ) {
135  dF_dT_data[dF_dT_.get_mem_index( j * T.size() + i )] = F_at_T_right_data[F_at_T_right.get_mem_index( j )]
136  .grad;
137  }
138  T_right_data[T_right_mem_idx].grad -= 1.0;
139  }
140 
141  F_ = makeReal( F_at_T_right );
142 
143  return { F_, dF_dT_ };
144  }
145 
146  namespace SecondOrder {
147 
148  template < size_t dim >
149  using tensor_to_scalar_function_type = std::function< dual2nd( const Fastor::Tensor< dual2nd, dim, dim >& T ) >;
150 
151  template < size_t dim >
152  std::tuple< double, Fastor::Tensor< double, dim, dim >, Fastor::Tensor< double, dim, dim, dim, dim > > d2f_dT2(
153  const tensor_to_scalar_function_type< dim >& F,
154  const Fastor::Tensor< double, dim, dim >& T )
155  {
156  double F_;
157  dual2nd F_right;
158  Fastor::Tensor< double, dim, dim > dF_dT_;
159  Fastor::Tensor< double, dim, dim, dim, dim > d2F_dT2;
160  Fastor::Tensor< dual2nd, dim, dim > T_right = makeHigherOrderDual< 2 >( T );
161 
162  for ( size_t i = 0; i < dim; i++ ) {
163  for ( size_t j = 0; j < dim; j++ ) {
164  seed< 1 >( T_right( i, j ), 1.0 );
165 
166  for ( size_t k = 0; k < dim; k++ ) {
167  for ( size_t l = 0; l < dim; l++ ) {
168 
169  seed< 2 >( T_right( k, l ), 1.0 );
170  F_right = F( T_right );
171  d2F_dT2( i, j, k, l ) = derivative< 2 >( F_right );
172 
173  seed< 2 >( T_right( k, l ), 0.0 );
174  }
175  }
176  dF_dT_( i, j ) = derivative< 1 >( F_right );
177  F_ = double( F_right );
178  seed< 1 >( T_right( i, j ), 0.0 );
179  }
180  }
181 
182  return { F_, dF_dT_, d2F_dT2 };
183  }
184 
185  template < size_t dim >
186  using tensor_and_scalar_to_scalar_function_type = std::function<
187  dual2nd( const Fastor::Tensor< dual2nd, dim, dim >& T, const dual2nd scalar ) >;
188 
189  template < size_t dim >
190  Fastor::Tensor< double, dim, dim > d2f_dTensor_dScalar( const tensor_and_scalar_to_scalar_function_type< dim >& F,
191  const Fastor::Tensor< double, dim, dim >& T,
192  const double scalar )
193  {
194  Fastor::Tensor< double, dim, dim > d2F_dTdScalar;
195  Fastor::Tensor< dual2nd, dim, dim > T_right = makeHigherOrderDual< 2 >( T );
196 
197  dual2nd scalar_right( scalar );
198  seed< 2 >( scalar_right, 1.0 );
199 
200  for ( size_t i = 0; i < dim; i++ ) {
201  for ( size_t j = 0; j < dim; j++ ) {
202 
203  seed< 1 >( T_right( i, j ), 1.0 );
204 
205  d2F_dTdScalar( i, j ) = derivative< 2 >( F( T_right, scalar_right ) );
206 
207  seed< 1 >( T_right( i, j ), 0.0 );
208  }
209  }
210 
211  return d2F_dTdScalar;
212  }
213  } // namespace SecondOrder
214 
215  } // namespace AutomaticDifferentiation
216 
217 } // namespace Marmot
Marmot::FastorIndices::i
Fastor::Index< i_ > i
Definition: MarmotFastorTensorBasics.h:137
Marmot::AutomaticDifferentiation::SecondOrder::d2f_dTensor_dScalar
Fastor::Tensor< double, dim, dim > d2f_dTensor_dScalar(const tensor_and_scalar_to_scalar_function_type< dim > &F, const Fastor::Tensor< double, dim, dim > &T, const double scalar)
Definition: MarmotAutomaticDifferentiationForFastor.h:190
Marmot::AutomaticDifferentiation::increaseDualOrderWithShift
autodiff::HigherOrderDual< order+1, double > increaseDualOrderWithShift(const autodiff::HigherOrderDual< order, double > &in)
Definition: MarmotAutomaticDifferentiation.h:69
MarmotFastorTensorBasics.h
Marmot::AutomaticDifferentiation::df_dT
Fastor::Tensor< double, Rest... > df_dT(const tensor_to_scalar_function_type< Rest... > &f, const Fastor::Tensor< double, Rest... > &T)
Definition: MarmotAutomaticDifferentiationForFastor.h:64
Marmot::FastorIndices::l
Fastor::Index< l_ > l
Definition: MarmotFastorTensorBasics.h:207
Marmot::AutomaticDifferentiation::dF_dT
std::pair< Fastor::Tensor< double, RestF... >, Fastor::Tensor< double, RestF..., RestT... > > dF_dT(std::function< Fastor::Tensor< dual, RestF... >(const Fastor::Tensor< dual, RestT... > &) > &F, const Fastor::Tensor< double, RestT... > &T)
Definition: MarmotAutomaticDifferentiationForFastor.h:114
Marmot::makeReal
Fastor::Tensor< double, Rest... > makeReal(const Fastor::Tensor< T, Rest... > &in)
Definition: MarmotFastorTensorBasics.h:399
Marmot::FastorIndices::in
Fastor::Index< i_, n_ > in
Definition: MarmotFastorTensorBasics.h:179
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::AutomaticDifferentiation::SecondOrder::tensor_and_scalar_to_scalar_function_type
std::function< dual2nd(const Fastor::Tensor< dual2nd, dim, dim > &T, const dual2nd scalar) > tensor_and_scalar_to_scalar_function_type
Definition: MarmotAutomaticDifferentiationForFastor.h:187
Marmot::makeDual
Fastor::Tensor< autodiff::dual, Rest... > makeDual(const Fastor::Tensor< T, Rest... > &in)
Definition: MarmotFastorTensorBasics.h:413
Marmot::AutomaticDifferentiation::SecondOrder::tensor_to_scalar_function_type
std::function< dual2nd(const Fastor::Tensor< dual2nd, dim, dim > &T) > tensor_to_scalar_function_type
Definition: MarmotAutomaticDifferentiationForFastor.h:149
Marmot::FastorIndices::k
Fastor::Index< k_ > k
Definition: MarmotFastorTensorBasics.h:195
Marmot::FiniteElement::EAS::F
Eigen::MatrixXd F(const Eigen::MatrixXd &J)
Marmot::AutomaticDifferentiation::tensor_to_scalar_function_type_arbitrary_dual_order
std::function< HigherOrderDual< order, double >(const Fastor::Tensor< HigherOrderDual< order, double >, Rest... > &T) > tensor_to_scalar_function_type_arbitrary_dual_order
Definition: MarmotAutomaticDifferentiationForFastor.h:85
Marmot::AutomaticDifferentiation::SecondOrder::d2f_dT2
std::tuple< double, Fastor::Tensor< double, dim, dim >, Fastor::Tensor< double, dim, dim, dim, dim > > d2f_dT2(const tensor_to_scalar_function_type< dim > &F, const Fastor::Tensor< double, dim, dim > &T)
Definition: MarmotAutomaticDifferentiationForFastor.h:152
MarmotAutomaticDifferentiation.h
Marmot::AutomaticDifferentiation::tensor_to_scalar_function_type
std::function< dual(const Fastor::Tensor< dual, Rest... > &T) > tensor_to_scalar_function_type
Definition: MarmotAutomaticDifferentiationForFastor.h:61