diff --git a/dune/elasticity/assemblers/localhyperdualstiffness.hh b/dune/elasticity/assemblers/localhyperdualstiffness.hh
index e3918385e15e2cd8d53cb8df61944011a2e472a4..c6ed85aa8623f1152315a418c445f14a5eced349 100644
--- a/dune/elasticity/assemblers/localhyperdualstiffness.hh
+++ b/dune/elasticity/assemblers/localhyperdualstiffness.hh
@@ -86,15 +86,14 @@ assembleGradientAndHessian(const LocalView& localView,
for (size_t i=0; i<nDoubles; i++)
{
localGradient[i] = temp.gradient()[i];
- localHessian[i][i] = temp.hessian()[i][i];
+ localHessian[i][i] = (temp.hessian())(i,i);
for (size_t j=i+1; j<nDoubles; j++)
{
- localHessian[i][j] = temp.hessian()[i][j];
+ localHessian[i][j] = (temp.hessian())(i,j);
localHessian[j][i] = localHessian[i][j]; // Use symmetry of hessian
}
}
-
}
} // namespace Dune::Elasticity
diff --git a/dune/elasticity/common/hyperdual.hh b/dune/elasticity/common/hyperdual.hh
index 6c6df1140860f6a29acb002095ef1145ffc24c83..3273aaf9bd6cdf401ede5a6f9340b4551c4f69ef 100644
--- a/dune/elasticity/common/hyperdual.hh
+++ b/dune/elasticity/common/hyperdual.hh
@@ -36,7 +36,7 @@
#include <math.h>
#include<dune/common/fmatrix.hh>
-#include <boost/hana/size.hpp>
+#include<dune/fufem/symmetricmatrix.hh>
class hyperdual{
@@ -44,7 +44,8 @@ public:
const static int size = 29;
- using Matrix = Dune::FieldMatrix<double,size,size>;
+// using Matrix = Dune::FieldMatrix<double,size,size>;
+ using Matrix = Dune::Fufem::SymmetricMatrix<double,size>;
using Vector = Dune::FieldVector<double,size>;
//creation operators and function to manually set values
@@ -195,7 +196,7 @@ hyperdual::Matrix hyperdual::hessian() const
std::ostream& operator<<(std::ostream& output, const hyperdual& rhs)
{
- output << "value: " << rhs.value() << "\ngradient: " << rhs.gradient() << "\nhessian: " << rhs.hessian();
+// output << "value: " << rhs.value() << "\ngradient: " << rhs.gradient() << "\nhessian: " << rhs.hessian();
return output;
}
@@ -297,7 +298,7 @@ hyperdual& hyperdual::operator*= (hyperdual rhs)
for(int i=0; i<size; i++)
{
for(int j=i; j<size; j++)
- hess[i][j] = rhs.grad[j]*grad[i] + rhs.val*hess[i][j] + grad[j]*rhs.grad[i] + val*rhs.hess[i][j];
+ hess.setEntry(i,j,rhs.grad[j]*grad[i] + rhs.val*hess(i,j) + grad[j]*rhs.grad[i] + val*rhs.hess(i,j));
grad[i] = val*rhs.grad[i] + rhs.val*grad[i];
}
@@ -322,9 +323,9 @@ hyperdual& hyperdual::operator/= (hyperdual rhs)
for(int i=0; i<size; i++)
{
for(int j=i; j<size; j++)
- hess[i][j] = -(rhs.grad[j]*grad[i] + rhs.grad[i]*grad[j] + val*rhs.hess[i][j])/(rhs.val*rhs.val)
- + hess[i][j]/rhs.val
- +2.*val*rhs.grad[i]*rhs.grad[j]/(rhs.val*rhs.val*rhs.val);
+ hess.setEntry(i,j,-(rhs.grad[j]*grad[i] + rhs.grad[i]*grad[j] + val*rhs.hess(i,j))/(rhs.val*rhs.val)
+ + hess(i,j)/rhs.val
+ +2.*val*rhs.grad[i]*rhs.grad[j]/(rhs.val*rhs.val*rhs.val));
grad[i] = grad[i]/rhs.val - val*rhs.grad[i]/(rhs.val*rhs.val);
}
@@ -337,8 +338,8 @@ hyperdual& hyperdual::operator/= (hyperdual rhs)
hyperdual& hyperdual::operator/= (double rhs)
{
val /= rhs;
- grad /= rhs;
- hess /= rhs;
+ grad *= (1./rhs);
+ hess *= (1./rhs);
return *this;
}
@@ -360,7 +361,7 @@ hyperdual pow (hyperdual x, double a)
for(int i=0; i<hyperdual::size; i++)
{
for(int j=i; j<hyperdual::size; j++)
- temp.hess[i][j] = aAMinusOneXPowMinusTwo*x.grad[i]*x.grad[j] + aXPowMinusOne*x.hess[i][j];
+ temp.hess.setEntry(i,j,aAMinusOneXPowMinusTwo*x.grad[i]*x.grad[j] + aXPowMinusOne*x.hess(i,j));
temp.grad[i] = aXPowMinusOne*x.grad[i];
}
@@ -383,7 +384,7 @@ hyperdual exp(hyperdual x)
for(int i=0; i<hyperdual::size; i++)
{
for(int j=i; j<hyperdual::size; j++)
- temp.hess[i][j] = eX*(x.grad[i]*x.grad[j] + x.hess[i][j]);
+ temp.hess.setEntry(i,j,eX*(x.grad[i]*x.grad[j] + x.hess(i,j)));
temp.grad[i] = eX*x.grad[i];
}
@@ -402,7 +403,7 @@ hyperdual log(hyperdual x)
for(int i=0; i<hyperdual::size; i++)
{
for(int j=i; j<hyperdual::size; j++)
- temp.hess[i][j] = -inv*inv*(x.grad[i]*x.grad[j]) + inv*x.hess[i][j];
+ temp.hess.setEntry(i,j,-inv*inv*(x.grad[i]*x.grad[j]) + inv*x.hess(i,j));
temp.grad[i] = inv*x.grad[i];
}
@@ -420,7 +421,7 @@ hyperdual sin(hyperdual x)
for(int i=0; i<hyperdual::size; i++)
{
for(int j=i; j<hyperdual::size; j++)
- temp.hess[i][j] = -sinX*(x.grad[i]*x.grad[j]) + cosX*x.hess[i][j];
+ temp.hess.setEntry(i,j,-sinX*(x.grad[i]*x.grad[j]) + cosX*x.hess(i,j));
temp.grad[i] = cosX*x.grad[i];
}
@@ -438,7 +439,7 @@ hyperdual cos(hyperdual x)
for(int i=0; i<hyperdual::size; i++)
{
for(int j=i; j<hyperdual::size; j++)
- temp.hess[i][j] = -cosX*(x.grad[i]*x.grad[j]) - sinX*x.hess[i][j];
+ temp.hess.setEntry(i,j,-cosX*(x.grad[i]*x.grad[j]) - sinX*x.hess(i,j));
temp.grad[i] = -sinX*x.grad[i];
}
@@ -455,7 +456,7 @@ hyperdual tan(hyperdual x)
for(int i=0; i<hyperdual::size; i++)
{
for(int j=i; j<hyperdual::size; j++)
- temp.hess[i][j] = 2.*tanX*deriv*(x.grad[i]*x.grad[j]) + deriv*x.hess[i][j];
+ temp.hess.setEntry(i,j,2.*tanX*deriv*(x.grad[i]*x.grad[j]) + deriv*x.hess(i,j));
temp.grad[i] = deriv*x.grad[i];
}
@@ -474,7 +475,7 @@ hyperdual asin(hyperdual x)
for(int i=0; i<hyperdual::size; i++)
{
for(int j=i; j<hyperdual::size; j++)
- temp.hess[i][j] = deriv2*(x.grad[i]*x.grad[j]) + deriv*x.hess[i][j];
+ temp.hess.setEntry(i,j,deriv2*(x.grad[i]*x.grad[j]) + deriv*x.hess(i,j));
temp.grad[i] = deriv*x.grad[i];
}
@@ -493,7 +494,7 @@ hyperdual acos(hyperdual x)
for(int i=0; i<hyperdual::size; i++)
{
for(int j=i; j<hyperdual::size; j++)
- temp.hess[i][j] = deriv2*(x.grad[i]*x.grad[j]) + deriv*x.hess[i][j];
+ temp.hess.setEntry(i,j,deriv2*(x.grad[i]*x.grad[j]) + deriv*x.hess(i,j));
temp.grad[i] = deriv*x.grad[i];
}
@@ -512,7 +513,7 @@ hyperdual atan(hyperdual x)
for(int i=0; i<hyperdual::size; i++)
{
for(int j=i; j<hyperdual::size; j++)
- temp.hess[i][j] = deriv2*(x.grad[i]*x.grad[j]) + deriv*x.hess[i][j];
+ temp.hess.setEntry(i,j,deriv2*(x.grad[i]*x.grad[j]) + deriv*x.hess(i,j));
temp.grad[i] = deriv*x.grad[i];
}
@@ -531,7 +532,7 @@ hyperdual sqrt(hyperdual x)
for(int i=0; i<hyperdual::size; i++)
{
for(int j=i; j<hyperdual::size; j++)
- temp.hess[i][j] = deriv2*(x.grad[i]*x.grad[j]) + deriv*x.hess[i][j];
+ temp.hess.setEntry(i,j,deriv2*(x.grad[i]*x.grad[j]) + deriv*x.hess(i,j));
temp.grad[i] = deriv*x.grad[i];
}