LocalADOLCStiffness: use "hessian2" instead of "hess_vec"

With this commit, we do two things:

-  replace "hess_vec" with a simpler call of  "hessian" which
   calls internally "hess_vec".

- use "hessian2" by default instead of "hessian". "hessian2" calls
  internally "hess_mat" which is faster than "hessian/hess_vec".
  If the user encounters probems with the new "hessian2" driver, a
  boolean "useHessian2" can be set to false in the constructor.

dune/elasticity/assemblers/localadolcstiffness.hh

@@ -27,8 +27,12 @@ class LocalADOLCStiffness
 public:
 
     // accepts ADOL_C's "adouble" only
-    LocalADOLCStiffness(std::shared_ptr<Dune::Elasticity::LocalEnergy<LocalView, adouble>> energy)
+    // useHessian2: there are different basic drivers for computing the hessian by ADOL-C, namely "hessian" and "hessian2".
+    //              They yield the same results, but "hessian2" seems to be faster. The user can switch to "hessian" with
+    //              useHessian2=false. See ADOL-C spec for more info about the drivers
+    LocalADOLCStiffness(std::shared_ptr<Dune::Elasticity::LocalEnergy<LocalView, adouble>> energy, bool useHessian2 = true)
     : localEnergy_(energy)
+    , useHessian2_(useHessian2)
     {}
 
     /** \brief Compute the energy at the current configuration */
@@ -45,6 +49,7 @@ public:
                                             Dune::Matrix<double>& localHessian);
 
     std::shared_ptr<Dune::Elasticity::LocalEnergy<LocalView, adouble>> localEnergy_;
+    const bool useHessian2_;
 
 };
 
@@ -110,29 +115,41 @@ assembleGradientAndHessian(const LocalView& localView,
 
     localHessian.setSize(nDoubles,nDoubles);
 
-    std::vector<double> v(nDoubles);
-    std::vector<double> w(nDoubles);
-
-    std::fill(v.begin(), v.end(), 0.0);
-
+    // allocate raw doubles: ADOL-C requires "**double" arguments
+    double* rawHessian[nDoubles];
     for(size_t i=0; i<nDoubles; i++)
     {
-      // Evaluate Hessian multiplied with the i-th canonical basis vector (i.e.: the i-th column of the Hessian)
-      v[i] = 1.;
+      // we need only the lower half of the hessian, thus allocate i+1 doubles only
+      rawHessian[i] = (double*)malloc((i+1)*sizeof(double));
+    }
+
+    // ADOL-C error code
+    int rc;
+
+    // ADOL-C won't accept pointers with const values, therefore the const_cast
+    if ( useHessian2_ )
+      rc = hessian2(rank,nDoubles,const_cast<double*>(localConfiguration.data()),rawHessian);
+    else
+      rc = hessian(rank,nDoubles,const_cast<double*>(localConfiguration.data()),rawHessian);
 
-      int rc= 3;
-      MINDEC(rc, hess_vec(rank, nDoubles, const_cast<double*>(localConfiguration.data()), v.data(), w.data())); // ADOL-C won't accept pointers with const values
     if (rc < 0)
       DUNE_THROW(Dune::Exception, "ADOL-C has returned with error code " << rc << "!");
 
-      for (size_t j=0; j<nDoubles; j++)
-        localHessian[i][j] = w[j];
-
-      // Make v the null vector again
-      v[i] = 0.;
+    // copy the raw values to the localHessian
+    for (size_t i=0; i<nDoubles; i++)
+    {
+      for (size_t j=0; j<i; j++)
+      {
+        // we can only access the lower half of rawHessian!
+        localHessian[i][j] = rawHessian[i][j];
+        localHessian[j][i] = rawHessian[i][j];
+      }
+      localHessian[i][i] = rawHessian[i][i];
     }
 
-    // TODO: implement ADOLC_VECTOR_MODE variant
+    // don't forget the clean everything up
+    for(size_t i=0; i<nDoubles; i++)
+      free(rawHessian[i]);
 }
 
 } // namespace Dune::Elasticity