MarmotMath.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  * Magdalena Schreter magdalena.schreter@uibk.ac.at
16  * Alexander Dummer alexander.dummer@uibk.ac.at
17  *
18  * This file is part of the MAteRialMOdellingToolbox (marmot).
19  *
20  * This library is free software; you can redistribute it and/or
21  * modify it under the terms of the GNU Lesser General Public
22  * License as published by the Free Software Foundation; either
23  * version 2.1 of the License, or (at your option) any later version.
24  *
25  * The full text of the license can be found in the file LICENSE.md at
26  * the top level directory of marmot.
27  * ---------------------------------------------------------------------
28  */
29 
30 #pragma once
31 #include "Marmot/MarmotConstants.h"
32 #include "Marmot/MarmotTypedefs.h"
33 #include "autodiff/forward/dual.hpp"
34 #include "autodiff/forward/real.hpp"
35 #include <algorithm>
36 #include <autodiff/forward/dual/dual.hpp>
37 #include <complex>
38 
39 namespace Marmot {
40  namespace Math {
43  template < typename T >
44  bool isNaN( T x )
45  {
46  return x != x;
47  }
48 
51  double linearInterpolation( double x, double x0, double x1, double y0, double y1 );
52 
55  double exp( double x );
56 
59  int getExponentPowerTen( const double x );
60 
63  constexpr double radToDeg( const double alpha )
64  {
65  return alpha * 180 / Marmot::Constants::Pi;
66  }
67 
70  constexpr double degToRad( const double alpha )
71  {
72  return alpha / 180 * Marmot::Constants::Pi;
73  }
74 
77  constexpr double macauly( double scalar )
78  {
79  return scalar >= 0 ? scalar : 0.0;
80  }
81 
84  constexpr int heaviside( double scalar )
85  {
86  return scalar >= 0 ? 1 : 0;
87  }
88 
91  template < typename T >
92  constexpr T sgn( T val )
93  {
94  return val / std::abs( val );
95  }
96 
97  double makeReal( const double& value );
98  double makeReal( const std::complex< double >& value );
99  double makeReal( const autodiff::real& value );
100 
101  template < typename T, typename G >
102  double makeReal( const autodiff::detail::Dual< T, G >& number )
103  {
104  return double( number );
105  }
106 
107  template < typename T, int... Rest >
108  Eigen::Matrix< double, Rest... > makeReal( const Eigen::Matrix< T, Rest... > mat )
109  {
110  Eigen::Matrix< double, Rest... > out;
111  const size_t m = static_cast< size_t >( mat.rows() );
112  const size_t n = static_cast< size_t >( mat.cols() );
113 
114  for ( size_t i = 0; i < m; i++ ) {
115  for ( size_t j = 0; j < n; j++ ) {
116  out( i, j ) = makeReal( mat( i, j ) );
117  }
118  }
119  return out;
120  }
121  template < typename T >
122  Eigen::VectorXd makeReal( Eigen::Vector< T, Eigen::Dynamic > in )
123  {
124 
125  int inSize = in.size();
126  Eigen::VectorXd out( inSize );
127  for ( int i = 0; i < inSize; i++ ) {
128  out( i ) = double( in( i ) );
129  }
130  return out;
131  }
132 
136  template < int nRows, int nCols >
137  Eigen::Matrix< double, nRows, nCols > macaulyMatrix( const Eigen::Matrix< double, nRows, nCols >& mat )
138  {
139  Eigen::Matrix< double, nRows, nCols > positivePart = Eigen::Matrix< double, nRows, nCols >::Zero();
140  for ( int i = 0; i < nRows; i++ ) {
141  for ( int j = 0; j < nCols; j++ ) {
142  positivePart( i, j ) = macauly( mat( i, j ) );
143  }
144  }
145  return positivePart;
146  }
147 
150  template < typename functionType, typename yType, typename... Args >
151  yType explicitEuler( yType yN, const double dt, functionType fRate, Args&&... fRateArgs )
152  {
153  return yN + fRate( yN, fRateArgs... ) * dt;
154  }
155 
161  template < int ySize, typename functionType, typename... Args >
162  Eigen::Matrix< double, ySize, 1 > semiImplicitEuler( Eigen::Matrix< double, ySize, 1 > yN,
163  const double dt,
164  functionType fRate,
165  Args&&... fRateArgs )
166  {
167  Eigen::Matrix< double, ySize, ySize > fS = Eigen::Matrix< double, ySize, ySize >::Zero();
168  Eigen::Matrix< double, ySize, ySize > Iy = Eigen::Matrix< double, ySize, ySize >::Identity();
169 
170  Eigen::VectorXd leftX( ySize );
171  Eigen::VectorXd rightX( ySize );
172 
173  for ( size_t i = 0; i < ySize; i++ ) {
174  double volatile h = std::max( 1.0, std::abs( yN( i ) ) ) * Constants::cubicRootEps();
175  leftX = yN;
176  leftX( i ) -= h;
177  rightX = yN;
178  rightX( i ) += h;
179  fS.col( i ) = 1. / ( 2. * h ) * ( fRate( rightX, fRateArgs... ) - fRate( leftX, fRateArgs... ) );
180  }
181 
182  return yN + ( Iy - dt * fS ).inverse() * dt * fRate( yN, fRateArgs... );
183  }
184 
188  template < int nRows, int nCols, typename functionType >
189  Eigen::Matrix< double, nRows, nCols > centralDiff( functionType f, const Eigen::Matrix< double, nCols, 1 >& X )
190  {
191  Eigen::Matrix< double, nRows, nCols > dXdY = Eigen::Matrix< double, nRows, nCols >::Zero();
192 
193  Eigen::VectorXd leftX( nCols );
194  Eigen::VectorXd rightX( nCols );
195 
196  for ( size_t i = 0; i < nCols; i++ ) {
197  double volatile h = std::max( 1.0, std::abs( X( i ) ) ) * Constants::cubicRootEps();
198  leftX = X;
199  leftX( i ) -= h;
200  rightX = X;
201  rightX( i ) += h;
202  dXdY.col( i ) = 1. / ( 2. * h ) * ( f( rightX ) - f( leftX ) );
203  }
204 
205  return dXdY;
206  }
207 
212  template < typename functionType, typename yType, typename... Args >
213  yType explicitEulerRichardson( yType yN, const double dt, functionType fRate, Args&&... fRateArgs )
214  {
215  yType fN = fRate( yN, fRateArgs... );
216  yType u = yN + fN * dt;
217  yType v = yN + fN * dt / 2.;
218  yType w = v + fRate( v, fRateArgs... ) * dt / 2.;
219 
220  return 2. * w - u;
221  }
222 
227  template < int ySize, typename functionType, typename... Args >
228  std::tuple< Eigen::Matrix< double, ySize, 1 >, double > explicitEulerRichardsonWithErrorEstimator(
229  Eigen::Matrix< double, ySize, 1 > yN,
230  const double dt,
231  const double TOL,
232  functionType fRate,
233  Args&&... fRateArgs )
234  {
235 
236  typedef Eigen::Matrix< double, ySize, 1 > ySized;
237  ySized fN = fRate( yN, fRateArgs... );
238  ySized u = yN + fN * dt;
239  ySized v = yN + fN * dt / 2.;
240  ySized w = v + fRate( v, fRateArgs... ) * dt / 2.;
241  ySized yNew = 2. * w - u;
242 
243  // error estimator
244  const double AERR = 1.0;
245  const double aI = AERR / TOL;
246  const double rI = 1.0;
247  double scaling = 0;
248  ySized ESTVec = ySized::Zero();
249 
250  for ( int i = 0; i < ySize; i++ ) {
251  scaling = aI + rI * std::max( abs( yNew( i ) ), abs( yN( i ) ) );
252  ESTVec( i ) = abs( w( i ) - u( i ) ) / abs( scaling );
253  }
254 
255  const double EST = ESTVec.maxCoeff();
256  const double tauNew = dt * std::min( 2., std::max( 0.2, 0.9 * std::sqrt( TOL / EST ) ) );
257 
258  return std::make_tuple( yNew, tauNew );
259  }
260 
264  Matrix3d directionCosines( const Matrix3d& transformedCoordinateSystem );
265 
270 
274  Matrix3d orthonormalCoordinateSystem( const Vector3d& n1, const Vector3d& n2 );
275  } // namespace Math
276 } // namespace Marmot
Marmot::Math::radToDeg
constexpr double radToDeg(const double alpha)
Definition: MarmotMath.h:63
Marmot::Math::explicitEulerRichardsonWithErrorEstimator
std::tuple< Eigen::Matrix< double, ySize, 1 >, double > explicitEulerRichardsonWithErrorEstimator(Eigen::Matrix< double, ySize, 1 > yN, const double dt, const double TOL, functionType fRate, Args &&... fRateArgs)
Definition: MarmotMath.h:228
Marmot::Math::macaulyMatrix
Eigen::Matrix< double, nRows, nCols > macaulyMatrix(const Eigen::Matrix< double, nRows, nCols > &mat)
Definition: MarmotMath.h:137
MarmotConstants.h
Marmot::Math::exp
double exp(double x)
Definition: MarmotMath.cpp:13
Marmot::Math::degToRad
constexpr double degToRad(const double alpha)
Definition: MarmotMath.h:70
Marmot::Math::sgn
constexpr T sgn(T val)
Definition: MarmotMath.h:92
MarmotTypedefs.h
Marmot::Math::isNaN
bool isNaN(T x)
Definition: MarmotMath.h:44
Marmot::Vector3d
Eigen::Matrix< double, 3, 1 > Vector3d
Definition: MarmotTypedefs.h:42
Marmot::Math::makeReal
double makeReal(const double &value)
Definition: MarmotMath.cpp:37
Marmot::Math::linearInterpolation
double linearInterpolation(double x, double x0, double x1, double y0, double y1)
Definition: MarmotMath.cpp:7
Marmot::Matrix3d
Eigen::Matrix< double, 3, 3 > Matrix3d
Definition: MarmotTypedefs.h:40
Marmot
This file includes functions needed for calculations with stress and strain tensors written in voigt ...
Definition: MarmotTesting.h:32
Marmot::Math::explicitEulerRichardson
yType explicitEulerRichardson(yType yN, const double dt, functionType fRate, Args &&... fRateArgs)
Definition: MarmotMath.h:213
Marmot::Math::explicitEuler
yType explicitEuler(yType yN, const double dt, functionType fRate, Args &&... fRateArgs)
Definition: MarmotMath.h:151
Marmot::Math::heaviside
constexpr int heaviside(double scalar)
Definition: MarmotMath.h:84
Marmot::Math::centralDiff
Eigen::Matrix< double, nRows, nCols > centralDiff(functionType f, const Eigen::Matrix< double, nCols, 1 > &X)
Definition: MarmotMath.h:189
Marmot::Math::orthonormalCoordinateSystem
Matrix3d orthonormalCoordinateSystem(Vector3d &normalVector)
Definition: MarmotMath.cpp:53
Marmot::Math::macauly
constexpr double macauly(double scalar)
Definition: MarmotMath.h:77
Marmot::Math::directionCosines
Matrix3d directionCosines(const Matrix3d &transformedCoordinateSystem)
Definition: MarmotMath.cpp:96
Marmot::Constants::cubicRootEps
double cubicRootEps()
Definition: MarmotConstants.h:36
Marmot::Constants::Pi
constexpr double Pi
Definition: MarmotConstants.h:34
Marmot::Math::semiImplicitEuler
Eigen::Matrix< double, ySize, 1 > semiImplicitEuler(Eigen::Matrix< double, ySize, 1 > yN, const double dt, functionType fRate, Args &&... fRateArgs)
Definition: MarmotMath.h:162
Marmot::Math::getExponentPowerTen
int getExponentPowerTen(const double x)
Definition: MarmotMath.cpp:43