From 695efe8938a123e9520b3024a0d1a6c9976d6709 Mon Sep 17 00:00:00 2001
From: Jonathan Drechsel <jonathan.drechsel@gmx.de>
Date: Wed, 24 Jun 2020 11:35:00 +0200
Subject: [PATCH] Implementation of LocalHyperDualStiffness
- implementation of hyper-dual numbers is added
- LocalHyperDualStiffness is a new alternative to LocalADOLCStiffness that uses hyper-dual numbers instead of AD by ADOL-C
- localhyperdualstiffnesstest verifies the results of LocalHyperDualStiffness
---
CHANGELOG.md | 1 +
dune/elasticity/assemblers/CMakeLists.txt | 1 +
.../assemblers/localhyperdualstiffness.hh | 109 +++
dune/elasticity/common/CMakeLists.txt | 1 +
dune/elasticity/common/hyperdual.hh | 634 ++++++++++++++++++
test/CMakeLists.txt | 3 +
test/localhyperdualstiffnesstest.cc | 114 ++++
7 files changed, 863 insertions(+)
create mode 100644 dune/elasticity/assemblers/localhyperdualstiffness.hh
create mode 100644 dune/elasticity/common/hyperdual.hh
create mode 100644 test/localhyperdualstiffnesstest.cc
diff --git a/CHANGELOG.md b/CHANGELOG.md
index 1c3b976..23d7579 100644
--- a/CHANGELOG.md
+++ b/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
diff --git a/dune/elasticity/assemblers/CMakeLists.txt b/dune/elasticity/assemblers/CMakeLists.txt
index fe27d4d..7425e0e 100644
--- a/dune/elasticity/assemblers/CMakeLists.txt
+++ b/dune/elasticity/assemblers/CMakeLists.txt
@@ -1,6 +1,7 @@
install(FILES
feassembler.hh
localadolcstiffness.hh
+ localhyperdualstiffness.hh
localenergy.hh
localfestiffness.hh
neohookefunctionalassembler.hh
diff --git a/dune/elasticity/assemblers/localhyperdualstiffness.hh b/dune/elasticity/assemblers/localhyperdualstiffness.hh
new file mode 100644
index 0000000..a91b68e
--- /dev/null
+++ b/dune/elasticity/assemblers/localhyperdualstiffness.hh
@@ -0,0 +1,109 @@
+#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
diff --git a/dune/elasticity/common/CMakeLists.txt b/dune/elasticity/common/CMakeLists.txt
index a46e7dc..df13ce7 100644
--- a/dune/elasticity/common/CMakeLists.txt
+++ b/dune/elasticity/common/CMakeLists.txt
@@ -1,6 +1,7 @@
install(FILES
directionalspringrobinfunction.hh
elasticityhelpers.hh
+ hyperdual.hh
maxnormtrustregion.hh
nonlinearelasticityproblem.hh
nonlinearelasticityproblem.cc
diff --git a/dune/elasticity/common/hyperdual.hh b/dune/elasticity/common/hyperdual.hh
new file mode 100644
index 0000000..524b785
--- /dev/null
+++ b/dune/elasticity/common/hyperdual.hh
@@ -0,0 +1,634 @@
+/*
+ * 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;
+ temp.f0 = funval;
+ temp.f1 = deriv*x.f1;
+ temp.f2 = deriv*x.f2;
+ temp.f12 = deriv*x.f12 + x.f1*x.f2*(-2*x.f0/(deriv1*deriv1));
+ return temp;
+}
+hyperdual sqrt(hyperdual x)
+{
+ return pow(x,0.5);
+}
+hyperdual fabs(hyperdual x)
+{
+ hyperdual temp;
+ if (x < 0.0)
+ temp = -x;
+ else
+ temp = x;
+ return temp;
+}
+
+hyperdual abs(hyperdual x)
+{
+ return fabs(x);
+}
+hyperdual max(hyperdual x1, hyperdual x2)
+{
+ hyperdual temp;
+ if (x1>x2)
+ temp = x1;
+ else
+ temp = x2;
+ return temp;
+}
+hyperdual max(hyperdual x1, double x2)
+{
+ hyperdual temp;
+ if (x1>x2)
+ temp = x1;
+ else
+ temp = x2;
+ return temp;
+}
+hyperdual max(double x1, hyperdual x2)
+{
+ hyperdual temp;
+ if (x1>x2)
+ temp = x1;
+ else
+ temp = x2;
+ return temp;
+}
+hyperdual min(hyperdual x1, hyperdual x2)
+{
+ hyperdual temp;
+ if (x1<x2)
+ temp = x1;
+ else
+ temp = x2;
+ return temp;
+}
+hyperdual min(hyperdual x1, double x2)
+{
+ hyperdual temp;
+ if (x1<x2)
+ temp = x1;
+ else
+ temp = x2;
+ return temp;
+}
+hyperdual min(double x1, hyperdual x2)
+{
+ hyperdual temp;
+ if (x1<x2)
+ temp = x1;
+ else
+ temp = x2;
+ return temp;
+}
+
+bool operator> (hyperdual lhs, hyperdual rhs)
+{
+ return (lhs.f0 > rhs.f0);
+}
+bool operator> (double lhs, hyperdual rhs)
+{
+ return (lhs > rhs.f0);
+}
+bool operator> (hyperdual lhs, double rhs)
+{
+ return (lhs.f0 > rhs);
+}
+bool operator>= (hyperdual lhs, hyperdual rhs)
+{
+ return (lhs.f0 >= rhs.f0);
+}
+bool operator>= (double lhs, hyperdual rhs)
+{
+ return (lhs >= rhs.f0);
+}
+bool operator>= (hyperdual lhs, double rhs)
+{
+ return (lhs.f0 >= rhs);
+}
+bool operator< (hyperdual lhs, hyperdual rhs)
+{
+ return (lhs.f0 < rhs.f0);
+}
+bool operator< (double lhs, hyperdual rhs)
+{
+ return (lhs < rhs.f0);
+}
+bool operator< (hyperdual lhs, double rhs)
+{
+ return (lhs.f0 < rhs);
+}
+bool operator<= (hyperdual lhs, hyperdual rhs)
+{
+ return (lhs.f0 <= rhs.f0);
+}
+bool operator<= (double lhs, hyperdual rhs)
+{
+ return (lhs <= rhs.f0);
+}
+bool operator<= (hyperdual lhs, double rhs)
+{
+ return (lhs.f0 <= rhs);
+}
+bool operator== (hyperdual lhs, hyperdual rhs)
+{
+ return (lhs.f0 == rhs.f0);
+}
+bool operator== (double lhs, hyperdual rhs)
+{
+ return (lhs == rhs.f0);
+}
+bool operator== (hyperdual lhs, double rhs)
+{
+ return (lhs.f0 == rhs);
+}
+bool operator!= (hyperdual lhs, hyperdual rhs)
+{
+ return (lhs.f0 != rhs.f0);
+}
+bool operator!= (double lhs, hyperdual rhs)
+{
+ return (lhs != rhs.f0);
+}
+bool operator!= (hyperdual lhs, double rhs)
+{
+ return (lhs.f0 != rhs);
+}
+
+
+// Dune related definitions with hyperdual numbers
+namespace Dune
+{
+ template<>
+ struct IsNumber<hyperdual>
+ : std::true_type
+ {};
+
+ /* Specializations of the PromotionTraits class from dune-common for the hyperdual type.
+ * This is needed, e.g., to be able to use the FieldVector implementation of dot products
+ * between a FieldVector<double> and a FieldVector<hyperdual>.
+ */
+ template <>
+ struct PromotionTraits<double,hyperdual>
+ {
+ using PromotedType = hyperdual;
+ };
+
+ template <>
+ struct PromotionTraits<hyperdual,double>
+ {
+ using PromotedType = hyperdual;
+ };
+}
+
+
+#endif
diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt
index f39aa17..6503516 100644
--- a/test/CMakeLists.txt
+++ b/test/CMakeLists.txt
@@ -7,3 +7,6 @@ dune_add_test(SOURCES materialtest
CMAKE_GUARD ADOLC_FOUND)
add_dune_ug_flags(materialtest)
add_dune_adolc_flags(materialtest)
+
+dune_add_test(SOURCES localhyperdualstiffnesstest
+ CMAKE_GUARD ADOLC_FOUND)
diff --git a/test/localhyperdualstiffnesstest.cc b/test/localhyperdualstiffnesstest.cc
new file mode 100644
index 0000000..f0c8887
--- /dev/null
+++ b/test/localhyperdualstiffnesstest.cc
@@ -0,0 +1,114 @@
+
+#include <config.h>
+
+#include <iostream>
+
+#include <dune/common/test/testsuite.hh>
+#include <dune/common/parametertree.hh>
+#include <dune/common/parametertreeparser.hh>
+#include <dune/common/timer.hh>
+
+#include <dune/grid/uggrid.hh>
+#include <dune/grid/utility/structuredgridfactory.hh>
+
+#include <dune/istl/bvector.hh>
+#include <dune/istl/bcrsmatrix.hh>
+
+#include <dune/functions/functionspacebases/lagrangebasis.hh>
+#include <dune/functions/functionspacebases/powerbasis.hh>
+
+#include <dune/elasticity/assemblers/feassembler.hh>
+#include <dune/elasticity/assemblers/localadolcstiffness.hh>
+#include <dune/elasticity/assemblers/localhyperdualstiffness.hh>
+#include <dune/elasticity/materials/stvenantkirchhoffdensity.hh>
+#include <dune/elasticity/materials/localintegralenergy.hh>
+
+using namespace Dune;
+
+template<int dim>
+TestSuite hyperdualEqualsADOL_C()
+{
+ TestSuite t;
+
+ using Vector = BlockVector<FieldVector<double, dim>>;
+ using Matrix = BCRSMatrix<FieldMatrix<double,dim,dim>>;
+
+ using GridType = UGGrid<dim>;
+
+ std::shared_ptr<GridType> grid;
+
+ FieldVector<double,dim> lower(0), upper(1);
+
+ std::array<unsigned int,dim> elements;
+ elements.fill(10);
+
+ grid = StructuredGridFactory<GridType>::createCubeGrid(lower, upper, elements);
+
+ using GridView = typename GridType::LeafGridView;
+ GridView gridView = grid->leafGridView();
+
+ using namespace Dune::Functions::BasisFactory;
+ auto basis = makeBasis(
+ gridView,
+ power<dim>(
+ lagrange<1>()
+ ));
+
+ auto localView = basis.localView();
+ using LocalView = decltype(localView);
+ using Basis = decltype(basis);
+
+ std::shared_ptr<Elasticity::LocalFEStiffness<LocalView,double>> localStiffness;
+
+ std::string parameterString = "mu = 2.7191e+6\nlambda = 4.4364e+6";
+ std::istringstream stream(parameterString);
+ ParameterTree parameterSet;
+ ParameterTreeParser::readINITree(stream, parameterSet);
+
+ auto hyperdualElasticDensity = std::make_shared<Elasticity::StVenantKirchhoffDensity<dim,hyperdual>>(parameterSet);
+ auto adolcElasticDensity = std::make_shared<Elasticity::StVenantKirchhoffDensity<dim,adouble>>(parameterSet);
+
+ auto hyperdualElasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<LocalView, hyperdual>>(hyperdualElasticDensity);
+ auto adolcElasticEnergy = std::make_shared<Elasticity::LocalIntegralEnergy<LocalView, adouble>>(adolcElasticDensity);
+
+ auto hyperdualLocalStiffness = std::make_shared<Elasticity::LocalHyperDualStiffness<LocalView>>(hyperdualElasticEnergy);
+ auto adolcLocalStiffness = std::make_shared<Elasticity::LocalADOLCStiffness<LocalView>>(adolcElasticEnergy);
+
+ Elasticity::FEAssembler<Basis,Vector> hyperdualAssembler(basis, hyperdualLocalStiffness);
+ Elasticity::FEAssembler<Basis,Vector> adolcAssembler(basis, adolcLocalStiffness);
+
+ Matrix hyperdualHessian, adolcHessian;
+ Vector hyperdualGradient, adolcGradient;
+
+ Vector x(basis.size());
+ x = 42.;
+
+ Timer timer;
+ hyperdualAssembler.assembleGradientAndHessian(x,hyperdualGradient,hyperdualHessian);
+ std::cout << dim <<"D: hyperdual assembler took " << timer.elapsed() << " sec." << std::endl;
+ timer.reset();
+ adolcAssembler.assembleGradientAndHessian(x,adolcGradient,adolcHessian);
+ std::cout << dim << "D: ADOL-C assembler took " << timer.elapsed() << " sec." << std::endl;
+
+
+ hyperdualHessian -= adolcHessian;
+ hyperdualGradient -= adolcGradient;
+
+ // check the relative error
+ t.check(hyperdualHessian.frobenius_norm()/adolcHessian.frobenius_norm() < 1e-12 && hyperdualGradient.two_norm()/adolcGradient.two_norm() < 1e-12);
+
+ return t;
+
+}
+
+int main(int argc, char **argv)
+{
+ Dune::MPIHelper::instance(argc, argv);
+
+ TestSuite t;
+
+ t.subTest(hyperdualEqualsADOL_C<2>());
+ t.subTest(hyperdualEqualsADOL_C<3>());
+
+ return t.exit();
+}
--
GitLab