WIP: Implementation of LocalHyperDualStiffness

CHANGELOG.md

@@ -3,6 +3,7 @@
 - Introduce class `LocalDensity` for material-specific implementations
 - Introduce class `LocalIntegralEnergy` to work with the densities
 - Local energies and FE assemblers use now `dune-functions` power bases instead of scalar `dune-fufem` bases; the key element is the `LocalView` which contains the information for each element
+- Introduce class `LocalHyperDualStiffness` and `hyperdual` to calculate gradient and hessian using hyper-dual numbers
 
 ## Deprecations and removals
 

dune/elasticity/assemblers/CMakeLists.txt

 install(FILES
     feassembler.hh
     localadolcstiffness.hh
+    localhyperdualstiffness.hh
     localenergy.hh
     localfestiffness.hh
     neohookefunctionalassembler.hh

dune/elasticity/assemblers/localhyperdualstiffness.hh

+#pragma once
+
+#include <dune/common/fmatrix.hh>
+#include <dune/istl/matrix.hh>
+
+#include <dune/elasticity/assemblers/localfestiffness.hh>
+#include <dune/elasticity/common/hyperdual.hh>
+
+namespace Dune::Elasticity {
+
+/** \brief Assembles energy gradient and Hessian using Hyperdual Numbers 
+ */
+template<class LocalView>
+class LocalHyperDualStiffness
+    : public LocalFEStiffness<LocalView,double>
+{
+    enum {gridDim=LocalView::GridView::dimension};
+
+public:
+
+    // accepts hyperdual only
+    LocalHyperDualStiffness(std::shared_ptr<Dune::Elasticity::LocalEnergy<LocalView, hyperdual>> 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, hyperdual>> localEnergy_; 
+
+};
+
+
+template<class LocalView>
+double LocalHyperDualStiffness<LocalView>::
+energy(const LocalView& localView,
+       const std::vector<double>& localConfiguration) const
+{
+    double pureEnergy;
+
+    std::vector<hyperdual> localHyperDualConfiguration(localConfiguration.size());
+    hyperdual energy;
+    
+    for (size_t i=0; i<localConfiguration.size(); i++)
+        localHyperDualConfiguration[i] = hyperdual(localConfiguration[i]);        
+        
+    energy = localEnergy_->energy(localView,localHyperDualConfiguration);
+
+    pureEnergy = energy.real();
+
+    return pureEnergy;
+}
+
+
+template<class LocalView>
+void LocalHyperDualStiffness<LocalView>::
+assembleGradientAndHessian(const LocalView& localView,
+                           const std::vector<double>& localConfiguration,
+                           std::vector<double>& localGradient,
+                           Dune::Matrix<double>& localHessian
+                          )
+{
+    size_t nDoubles = localConfiguration.size(); 
+
+    localGradient.resize(nDoubles);
+
+    localHessian.setSize(nDoubles,nDoubles);
+    // Create hyperdual vector localHyperDualConfiguration = double vector localConfiguration
+    std::vector<hyperdual> localHyperDualConfiguration(nDoubles);
+    for (size_t i=0; i<nDoubles; i++)
+       localHyperDualConfiguration[i] = hyperdual(localConfiguration[i]);
+    
+    // Compute gradient and hessian (symmetry is used) using hyperdual function evaluation
+    // localHyperDualConfiguration puts Ones in the eps1 and eps2 component to evaluate the different partial derivatives
+    for (size_t i=0; i<nDoubles; i++) 
+    {
+      localHyperDualConfiguration[i].seteps1(1.0);
+      localHyperDualConfiguration[i].seteps2(1.0);
+    
+      hyperdual temp = localEnergy_->energy(localView, localHyperDualConfiguration);
+      localGradient[i] = temp.eps1();                   
+      localHessian[i][i] = temp.eps1eps2();
+      
+      localHyperDualConfiguration[i].seteps2(0.0);
+
+      for (size_t j=i+1; j<nDoubles; j++) 
+      { 
+        localHyperDualConfiguration[j].seteps2(1.0);
+        localHessian[i][j] = localEnergy_->energy(localView, localHyperDualConfiguration).eps1eps2();
+        localHessian[j][i] = localHessian[i][j];        // Use symmetry of hessian 
+        localHyperDualConfiguration[j].seteps2(0.0);
+      }
+        
+      localHyperDualConfiguration[i].seteps1(0.0);
+    }
+
+}
+
+} // namespace Dune::Elasticity

dune/elasticity/common/CMakeLists.txt

 install(FILES
     directionalspringrobinfunction.hh
     elasticityhelpers.hh
+    hyperdual.hh
     maxnormtrustregion.hh
     nonlinearelasticityproblem.hh
     nonlinearelasticityproblem.cc

dune/elasticity/common/hyperdual.hh

+/*
+ * Class: hyperdual
+ * 
+ * Implementation of hyper-dual numbers
+ *
+ * Written by: Jeffrey A. Fike
+ * Stanford University, Department of Aeronautics and Astronautics
+ * 
+ * Copyright (c) 2006 Jeffrey A. Fike
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a copy
+ * of this software and associated documentation files (the "Software"), to deal
+ * in the Software without restriction, including without limitation the rights
+ * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+ * copies of the Software, and to permit persons to whom the Software is
+ * furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included in
+ * all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+ * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+ * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+ * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
+ * THE SOFTWARE.
+ *
+ */
+
+
+#ifndef DUNE_ELASTICITY_COMMON_HYPERDUAL_HH
+#define DUNE_ELASTICITY_COMMON_HYPERDUAL_HH
+
+#include <iostream>
+#include <math.h>
+
+class hyperdual{
+	double f0,f1,f2,f12;
+public:
+	//creation operators and function to manually set values
+	hyperdual();
+	hyperdual(double x1,double x2,double x3,double x4);
+	hyperdual(double x1);
+	void setvalues(double x1,double x2,double x3,double x4);
+	void setreal(double x1);
+	void seteps1(double x2);
+	void seteps2(double x3);
+	void seteps12(double x4);
+
+	//examine values
+	void view(void);
+	double real(void) const;
+	double eps1(void) const;
+	double eps2(void) const;
+	double eps1eps2(void) const;
+	friend std::ostream& operator<<(std::ostream& output, const hyperdual& rhs);
+
+	//basic manipulation
+	hyperdual operator+ () const;
+	hyperdual operator+ (const hyperdual rhs) const;
+	friend hyperdual operator+ (const double lhs, const hyperdual rhs);
+	hyperdual operator- () const;
+	hyperdual operator- (const hyperdual rhs) const;
+	friend hyperdual operator- (const double lhs, const hyperdual rhs);
+	hyperdual operator* (const hyperdual rhs)const;
+	friend hyperdual operator* (const double lhs, const hyperdual rhs);
+	friend hyperdual operator/ (const hyperdual lhs, const hyperdual rhs);
+	friend hyperdual operator/ (const double lhs, const hyperdual rhs);
+	friend hyperdual operator/ (const hyperdual lhs, const double rhs);
+	hyperdual& operator+= (hyperdual rhs);
+	hyperdual& operator-= (hyperdual rhs);
+	hyperdual& operator*= (hyperdual rhs);
+	hyperdual& operator*= (double rhs);
+	hyperdual& operator/= (hyperdual rhs);
+	hyperdual& operator/= (double rhs);
+
+	//math.h functions
+	friend hyperdual pow (hyperdual x, double a);
+	friend hyperdual pow (hyperdual x, hyperdual a);
+	friend hyperdual exp(hyperdual x);
+	friend hyperdual log(hyperdual x);
+	friend hyperdual sin(hyperdual x);
+	friend hyperdual cos(hyperdual x);
+	friend hyperdual tan(hyperdual x);
+	friend hyperdual asin(hyperdual x);
+	friend hyperdual acos(hyperdual x);
+	friend hyperdual atan(hyperdual x);
+	friend hyperdual sqrt(hyperdual x);
+	friend hyperdual fabs(hyperdual x);
+	friend hyperdual abs(hyperdual x);  
+	friend hyperdual max(hyperdual x1, hyperdual x2);
+	friend hyperdual max(hyperdual x1, double x2);
+	friend hyperdual max(double x1, hyperdual x2);
+	friend hyperdual min(hyperdual x1, hyperdual x2);
+	friend hyperdual min(hyperdual x1, double x2);
+	friend hyperdual min(double x1, hyperdual x2);
+
+	//comparisons
+	friend bool operator> (hyperdual lhs, hyperdual rhs);
+	friend bool operator> (double lhs, hyperdual rhs);
+	friend bool operator> (hyperdual lhs, double rhs);
+	friend bool operator>= (hyperdual lhs, hyperdual rhs);
+	friend bool operator>= (double lhs, hyperdual rhs);
+	friend bool operator>= (hyperdual lhs, double rhs);
+	friend bool operator< (hyperdual lhs, hyperdual rhs);
+	friend bool operator< (double lhs, hyperdual rhs);
+	friend bool operator< (hyperdual lhs, double rhs);
+	friend bool operator<= (hyperdual lhs, hyperdual rhs);
+	friend bool operator<= (double lhs, hyperdual rhs);
+	friend bool operator<= (hyperdual lhs, double rhs);
+	friend bool operator== (hyperdual lhs, hyperdual rhs);
+	friend bool operator== (double lhs, hyperdual rhs);
+	friend bool operator== (hyperdual lhs, double rhs);
+	friend bool operator!= (hyperdual lhs, hyperdual rhs);
+	friend bool operator!= (double lhs, hyperdual rhs);
+	friend bool operator!= (hyperdual lhs, double rhs);
+};
+
+
+hyperdual::hyperdual()
+{
+	f0 = 0.0;
+	f1 = 0.0;
+	f2 = 0.0;
+	f12 = 0.0;
+}
+hyperdual::hyperdual(double x1,double x2,double x3,double x4)
+{
+	f0 = x1;
+	f1 = x2;
+	f2 = x3;
+	f12 = x4;
+}
+hyperdual::hyperdual(double x1)
+{
+	f0 = x1;
+	f1 = 0.0;
+	f2 = 0.0;
+	f12 = 0.0;
+}
+void hyperdual::setvalues(double x1,double x2,double x3,double x4)
+{
+	f0 = x1;
+	f1 = x2;
+	f2 = x3;
+	f12 = x4;
+}
+void hyperdual::setreal(double x1)
+{
+    f0 = x1;
+}
+void hyperdual::seteps1(double x2)
+{
+    f1 = x2;
+}
+void hyperdual::seteps2(double x3)
+{
+    f2 = x3;
+}
+void hyperdual::seteps12(double x4)
+{
+    f12 = x4;
+}
+void hyperdual::view(void)
+{
+	printf("%g  +  %g epsilon1  +  %g epsilon2  +  %g epsilon1 epsilon2\n",f0,f1,f2,f12);
+}
+double hyperdual::real(void) const
+{
+	return f0;
+}
+double hyperdual::eps1(void) const
+{
+	return f1;
+}
+double hyperdual::eps2(void) const
+{
+	return f2;
+}
+double hyperdual::eps1eps2(void) const
+{
+	return f12;
+}
+std::ostream& operator<<(std::ostream& output, const hyperdual& rhs)
+{
+        output << "(" << rhs.f0 << ","<< rhs.f1 << ","<< rhs.f2 << ","<< rhs.f12 << ")";
+        return output;
+}
+
+hyperdual hyperdual::operator+ () const
+{
+	return *this;
+}
+hyperdual hyperdual::operator+ (const hyperdual rhs) const
+{
+	hyperdual temp;
+	temp.f0 = f0 + rhs.f0;
+	temp.f1 = f1 + rhs.f1;
+	temp.f2 = f2 + rhs.f2;
+	temp.f12 = f12 + rhs.f12;
+	return temp;
+}
+hyperdual operator+ (const double lhs, const hyperdual rhs)
+{
+	hyperdual temp;
+	temp.f0 = lhs + rhs.f0;
+	temp.f1 = rhs.f1;
+	temp.f2 = rhs.f2;
+	temp.f12 = rhs.f12;
+	return temp;
+}
+hyperdual hyperdual::operator- () const
+{
+	hyperdual temp;
+	temp.f0 = -f0;
+	temp.f1 = -f1;
+	temp.f2 = -f2;
+	temp.f12 = -f12;
+	return temp;
+}
+hyperdual hyperdual::operator- (const hyperdual rhs) const
+{
+	hyperdual temp;
+	temp.f0 = f0 - rhs.f0;
+	temp.f1 = f1 - rhs.f1;
+	temp.f2 = f2 - rhs.f2;
+	temp.f12 = f12 - rhs.f12;
+	return temp;
+}
+hyperdual operator- (const double lhs, const hyperdual rhs)
+{
+	hyperdual temp;
+	temp.f0 = lhs - rhs.f0;
+	temp.f1 = -rhs.f1;
+	temp.f2 = -rhs.f2;
+	temp.f12 = -rhs.f12;
+	return temp;
+}
+hyperdual hyperdual::operator* (const hyperdual rhs) const
+{
+	hyperdual temp;
+	temp.f0 = f0*rhs.f0;
+	temp.f1 = f0*rhs.f1 + f1*rhs.f0;
+	temp.f2 = f0*rhs.f2 + f2*rhs.f0;
+	temp.f12 = f0*rhs.f12 + f1*rhs.f2 + f2*rhs.f1 + f12*rhs.f0;
+	return temp;
+}
+hyperdual operator* (const double lhs, const hyperdual rhs)
+{
+	hyperdual temp;
+	temp.f0 = lhs*rhs.f0;
+	temp.f1 = lhs*rhs.f1;
+	temp.f2 = lhs*rhs.f2;
+	temp.f12 = lhs*rhs.f12;
+	return temp;
+}
+hyperdual operator/ (const hyperdual lhs, const hyperdual rhs)
+{
+	hyperdual inv;
+	double prod;
+
+	prod = rhs.f0*rhs.f0;
+	inv.f0 = 1.0/rhs.f0;
+	inv.f1 = -rhs.f1/prod;
+	inv.f2 = -rhs.f2/prod;
+	inv.f12 = 2.0*rhs.f1*rhs.f2/(prod*rhs.f0) - rhs.f12/prod;
+	return lhs*inv;
+}
+hyperdual operator/ (const double lhs, const hyperdual rhs)
+{
+	hyperdual temp,inv;
+	temp = hyperdual(lhs)/rhs;
+	return temp;
+}
+hyperdual operator/ (const hyperdual lhs, const double rhs)
+{
+	hyperdual temp;
+	double inv;
+	inv = 1.0/rhs;
+	temp.f0 = inv*lhs.f0;
+	temp.f1 = inv*lhs.f1;
+	temp.f2 = inv*lhs.f2;
+	temp.f12 = inv*lhs.f12;
+	return temp;
+}
+hyperdual& hyperdual::operator+= (hyperdual rhs)
+{
+	f0 += rhs.f0;
+	f1 += rhs.f1;
+	f2 += rhs.f2;
+	f12 += rhs.f12;
+	return *this;
+}
+hyperdual& hyperdual::operator-= (hyperdual rhs)
+{
+	f0 -= rhs.f0;
+	f1 -= rhs.f1;
+	f2 -= rhs.f2;
+	f12 -= rhs.f12;
+	return *this;
+}
+hyperdual& hyperdual::operator*= (hyperdual rhs)
+{
+	double tf0,tf1,tf2,tf12;
+	tf0=f0;
+	tf1=f1;
+	tf2=f2;
+	tf12=f12;
+	f0 = tf0*rhs.f0;
+	f1 = tf0*rhs.f1 + tf1*rhs.f0;
+	f2 = tf0*rhs.f2 + tf2*rhs.f0;
+	f12 = tf0*rhs.f12 + tf1*rhs.f2 + tf2*rhs.f1 + tf12*rhs.f0;
+	return *this;
+}
+hyperdual& hyperdual::operator*= (double rhs)
+{
+	f0 *= rhs;
+	f1 *= rhs;
+	f2 *= rhs;
+	f12 *= rhs;
+	return *this;
+}
+hyperdual& hyperdual::operator/= (hyperdual rhs) 
+{
+    *this = *this / rhs;
+    return *this;
+}
+hyperdual& hyperdual::operator/= (double rhs)
+{
+	f0 /= rhs;
+	f1 /= rhs;
+	f2 /= rhs;
+	f12 /= rhs;
+	return *this;
+}
+hyperdual pow (hyperdual x, double a)
+{
+	hyperdual temp;
+	double deriv,xval,tol;
+	xval = x.f0;
+	tol = 1e-15;
+	if (fabs(xval) < tol)
+	{
+		if (xval >= 0)
+			xval = tol;
+		if (xval < 0)
+			xval = -tol;
+	}
+	deriv = a*std::pow(xval,(a-1));
+	//temp.f0 = pow(xval,a);
+	temp.f0 = std::pow(x.f0,a);  //Use actual x value, only use tol for derivs
+	temp.f1 = x.f1*deriv;
+	temp.f2 = x.f2*deriv;
+	temp.f12 = x.f12*deriv + a*(a-1)*x.f1*x.f2*std::pow(xval,(a-2));
+
+	return temp;
+}
+hyperdual pow (hyperdual x, hyperdual a)
+{
+	return exp(a*log(x));
+}
+hyperdual exp(hyperdual x)
+{
+	hyperdual temp;
+	double deriv;
+	deriv = std::exp(x.f0);
+	temp.f0 = deriv;
+	temp.f1 = deriv*x.f1;
+	temp.f2 = deriv*x.f2;
+	temp.f12 = deriv*(x.f12 + x.f1*x.f2);
+	return temp;
+}
+hyperdual log(hyperdual x)
+{
+	hyperdual temp;
+	double deriv1,deriv2;
+	deriv1 = x.f1/x.f0;
+	deriv2 = x.f2/x.f0;
+	temp.f0 = std::log(x.f0);
+	temp.f1 = deriv1;
+	temp.f2 = deriv2;
+	temp.f12 = x.f12/x.f0 - (deriv1*deriv2);
+	return temp;
+}
+hyperdual sin(hyperdual x)
+{
+	hyperdual temp;
+	double funval,deriv;
+	funval = std::sin(x.f0);
+	deriv = std::cos(x.f0);
+	temp.f0 = funval;
+	temp.f1 = deriv*x.f1;
+	temp.f2 = deriv*x.f2;
+	temp.f12 = deriv*x.f12 - funval*x.f1*x.f2;
+	return temp;
+}
+hyperdual cos(hyperdual x)
+{
+	hyperdual temp;
+	double funval,deriv;
+	funval = std::cos(x.f0);
+	deriv = -std::sin(x.f0);
+	temp.f0 = funval;
+	temp.f1 = deriv*x.f1;
+	temp.f2 = deriv*x.f2;
+	temp.f12 = deriv*x.f12 - funval*x.f1*x.f2;
+	return temp;
+}
+hyperdual tan(hyperdual x)
+{
+	hyperdual temp;
+	double funval,deriv;
+	funval = std::tan(x.f0);
+	deriv  = funval*funval + 1.0;
+	temp.f0 = funval;
+	temp.f1 = deriv*x.f1;
+	temp.f2 = deriv*x.f2;
+	temp.f12 = deriv*x.f12 + x.f1*x.f2*(2*funval*deriv);
+	return temp;
+}
+hyperdual asin(hyperdual x)
+{
+	hyperdual temp;
+	double funval,deriv1,deriv;
+	funval = std::asin(x.f0);
+	deriv1 = 1.0-x.f0*x.f0;
+	deriv = 1.0/std::sqrt(deriv1);
+	temp.f0 = funval;
+	temp.f1 = deriv*x.f1;
+	temp.f2 = deriv*x.f2;
+	temp.f12 = deriv*x.f12 + x.f1*x.f2*(x.f0*std::pow(deriv1,-1.5));
+	return temp;
+}
+hyperdual acos(hyperdual x)
+{
+	hyperdual temp;
+	double funval,deriv1,deriv;
+	funval = std::acos(x.f0);
+	deriv1 = 1.0-x.f0*x.f0;
+	deriv = -1.0/std::sqrt(deriv1);
+	temp.f0 = funval;
+	temp.f1 = deriv*x.f1;
+	temp.f2 = deriv*x.f2;
+	temp.f12 = deriv*x.f12 + x.f1*x.f2*(-x.f0*std::pow(deriv1,-1.5));
+	return temp;
+}
+hyperdual atan(hyperdual x)
+{
+	hyperdual temp;
+	double funval,deriv1,deriv;
+	funval = std::atan(x.f0);
+	deriv1 = 1.0+x.f0*x.f0;
+	deriv = 1.0/deriv1;