diff --git a/dune.module b/dune.module index 8bb8868b3ed53af6907535634bdc8eebbb9c5e24..0106adffc5280db4915ad44ed3f90e1e2d1dae4f 100644 --- a/dune.module +++ b/dune.module @@ -7,5 +7,5 @@ Module:dune-elasticity Version: 2.8-git Maintainer: oliver.sander@tu-dresden.de, youett@mi.fu-berlin.de #depending on -Depends:dune-common dune-grid dune-istl dune-solvers dune-fufem dune-geometry dune-functions +Depends:dune-common dune-grid dune-istl dune-solvers dune-fufem dune-geometry dune-functions dune-autodiff Suggests: dune-uggrid dune-parmg diff --git a/dune/elasticity/assemblers/localmagicdoublestiffness.hh b/dune/elasticity/assemblers/localmagicdoublestiffness.hh new file mode 100644 index 0000000000000000000000000000000000000000..af72d52e259abc5ec60ff96ea5fe16fd4c715e31 --- /dev/null +++ b/dune/elasticity/assemblers/localmagicdoublestiffness.hh @@ -0,0 +1,113 @@ +#pragma once + +#include <dune/common/fmatrix.hh> +#include <dune/istl/matrix.hh> + +#include <dune/elasticity/assemblers/localfestiffness.hh> + +#include <dune/autodiff/magic.hh> + +namespace Dune::Elasticity { + + +/** \brief Assembles energy gradient and Hessian using magic Numbers + */ +template<class LocalView> +class LocalMagicDoubleStiffness + : public LocalFEStiffness<LocalView,double> +{ + enum {gridDim=LocalView::GridView::dimension}; + +public: + + LocalMagicDoubleStiffness(std::shared_ptr<Dune::Elasticity::LocalEnergy<LocalView, Magic::magicdouble>> energy) + : localEnergy_(energy) + {} + + /** \brief Compute the energy at the current configuration */ + virtual double energy (const LocalView& localView, + const std::vector<double>& localConfiguration) const; + + /** \brief Assemble the local stiffness matrix at the current position + * + * This uses the automatic differentiation using hyperdual numbers + */ + virtual void assembleGradientAndHessian(const LocalView& localView, + const std::vector<double>& localConfiguration, + std::vector<double>& localGradient, + Dune::Matrix<double>& localHessian); + + std::shared_ptr<Dune::Elasticity::LocalEnergy<LocalView, Magic::magicdouble>> localEnergy_; + +}; + + +template<class LocalView> +double LocalMagicDoubleStiffness<LocalView>:: +energy(const LocalView& localView, + const std::vector<double>& localConfiguration) const +{ + std::vector<Magic::magicdouble> localMagicDoubleConfiguration(localConfiguration.size()); + + for (size_t i=0; i<localConfiguration.size(); i++) + localMagicDoubleConfiguration[i] = Magic::magicdouble(localConfiguration[i]); + + auto energy = localEnergy_->energy(localView,localMagicDoubleConfiguration); + + return energy.value(); +} + + +template<class LocalView> +void LocalMagicDoubleStiffness<LocalView>:: +assembleGradientAndHessian(const LocalView& localView, + const std::vector<double>& localConfiguration, + std::vector<double>& localGradient, + Dune::Matrix<double>& localHessian + ) +{ + std::cerr << " assembleGradientAndHessian\n"; + + size_t nDoubles = localConfiguration.size(); + + localGradient.resize(nDoubles); + + std::vector<double> direction1(nDoubles,0.0); + std::vector<double> direction2(nDoubles,0.0); + + localHessian.setSize(nDoubles,nDoubles); + std::vector<Magic::magicdouble> localMagicDoubleConfiguration; + for (size_t i=0; i<nDoubles; i++) + localMagicDoubleConfiguration.push_back( Magic::magicdouble(Magic::inputdouble( &localConfiguration[i], &direction1[i], &direction2[i] ) ) ); + +// std::cerr << " check1\n"; + + // this is done only once + auto temp = localEnergy_->energy(localView, localMagicDoubleConfiguration); + +// std::cerr << " check2\n"; + + for (size_t i=0; i<nDoubles; i++) + { + direction1[i] = 1.0; + direction2[i] = 1.0; + + localGradient[i] = temp.gradient(); + localHessian[i][i] = temp.hessian(); + + direction2[i] = 0.0; + + for (size_t j=i+1; j<nDoubles; j++) + { + direction2[j] = 1.0; + localHessian[i][j] = temp.hessian(); + localHessian[j][i] = localHessian[i][j]; // Use symmetry of hessian + direction2[j] = 0.0; + } + + direction1[i] = 0.0; + } + +} + +} // namespace Dune::Elasticity