WIP: magic number autodiff

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

dune/elasticity/assemblers/localmagicdoublestiffness.hh

+#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