WIP: use SymmetricMatrix from fufem for hessian

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

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];
   }