Use the vectors of the orthonormal frame as directions in hess_mat

Previously we called hess_mat (vector-mode second derivatives) with the
canonical basis vectors, and then rotated the resulting matrix.
This is wasteful, as more directions are used than necessary.
This patch uses the minimal number of directions.

[[Imported from SVN: r9818]]

dune/gfe/localgeodesicfeadolcstiffness.hh

@@ -244,9 +244,9 @@ assembleGradientAndHessian(const Entity& element,
     for (size_t i=0; i<nDofs; i++)
         orthonormalFrame[i] = localSolution[i].orthonormalFrame();
 
-#ifndef ADOLC_VECTOR_MODE
     Dune::Matrix<Dune::FieldMatrix<double,blocksize, embeddedBlocksize> > embeddedHessian(nDofs,nDofs);
 
+#ifndef ADOLC_VECTOR_MODE
     std::vector<double> v(nDoubles);
     std::vector<double> w(nDoubles);
 
@@ -270,42 +270,9 @@ assembleGradientAndHessian(const Entity& element,
         // Make v the null vector again
         std::fill(&v[i*embeddedBlocksize], &v[(i+1)*embeddedBlocksize], 0.0);
       }
-
-    // From this, compute the Hessian with respect to the manifold (which we assume here is embedded
-    // isometrically in a Euclidean space.
-    // For the detailed explanation of the following see: Absil, Mahoney, Trumpf, "An extrinsic look
-    // at the Riemannian Hessian".
-
-    typedef typename TargetSpace::EmbeddedTangentVector EmbeddedTangentVector;
-    typedef typename TargetSpace::TangentVector         TangentVector;
-
-    this->A_.setSize(nDofs,nDofs);
-
-    for (size_t col=0; col<nDofs; col++) {
-
-        for (size_t subCol=0; subCol<blocksize; subCol++) {
-
-            EmbeddedTangentVector z = orthonormalFrame[col][subCol];
-
-            // P_x \partial^2 f z
-            for (size_t row=0; row<nDofs; row++) {
-                TangentVector semiEmbeddedProduct;
-                embeddedHessian[row][col].mv(z,semiEmbeddedProduct);
-
-                for (int subRow=0; subRow<blocksize; subRow++)
-                    this->A_[row][col][subRow][subCol] = semiEmbeddedProduct[subRow];
-            }
-
-        }
-
-    }
-
 #else
-    Dune::Matrix<Dune::FieldMatrix<double,embeddedBlocksize, embeddedBlocksize> > embeddedHessian(nDofs,nDofs);
-
     int n = nDoubles;
-    //std::cout << "n: " << n << std::endl;
-    int nDirections = nDoubles;
+    int nDirections = nDofs * blocksize;
     double* tangent[nDoubles];
     for(size_t i=0; i<nDoubles; i++)
         tangent[i] = (double*)malloc(nDirections*sizeof(double));
@@ -314,26 +281,27 @@ assembleGradientAndHessian(const Entity& element,
     for(size_t i=0; i<nDoubles; i++)
         rawHessian[i] = (double*)malloc(nDirections*sizeof(double));
 
-    //std::cout << "nDirections: " << nDirections << std::endl;
-    for (int i=0; i<n; i++)
-      for (int j=0; j<nDirections; j++)
-        tangent[i][j] = (i==j);
+    for (int j=0; j<nDirections; j++)
+    {
+      for (int i=0; i<n; i++)
+        tangent[i][j] = 0.0;
 
-    hess_mat(1,nDoubles,nDirections,xp.data(),tangent,rawHessian);
+      for (int i=0; i<embeddedBlocksize; i++)
+        tangent[(j/blocksize)*embeddedBlocksize+i][j] = orthonormalFrame[j/blocksize][j%blocksize][i];
+    }
 
-    //std::cout << "Hallo Welt!" << std::endl;
+    hess_mat(1,nDoubles,nDirections,xp.data(),tangent,rawHessian);
 
     // Copy Hessian into Dune data type
     for(size_t i=0; i<nDoubles; i++)
-      for (size_t j=0; j<nDoubles; j++)
-        embeddedHessian[i/embeddedBlocksize][j/embeddedBlocksize][i%embeddedBlocksize][j%embeddedBlocksize] = rawHessian[i][j];
-
+      for (size_t j=0; j<nDirections; j++)
+        embeddedHessian[j/blocksize][i/embeddedBlocksize][j%blocksize][i%embeddedBlocksize] = rawHessian[i][j];
 
     for(size_t i=0; i<nDoubles; i++) {
         free(rawHessian[i]);
         free(tangent[i]);
     }
-    //printmatrix(std::cout, embeddedHessian, "hessian", "--");
+#endif
 
     // From this, compute the Hessian with respect to the manifold (which we assume here is embedded
     // isometrically in a Euclidean space.
@@ -352,22 +320,17 @@ assembleGradientAndHessian(const Entity& element,
             EmbeddedTangentVector z = orthonormalFrame[col][subCol];
 
             // P_x \partial^2 f z
-            std::vector<EmbeddedTangentVector> semiEmbeddedProduct(nDofs);
-
             for (size_t row=0; row<nDofs; row++) {
-                embeddedHessian[row][col].mv(z,semiEmbeddedProduct[row]);
-                //tmp1 = localSolution[row].projectOntoTangentSpace(tmp1);
-                TangentVector tmp2;
-                orthonormalFrame[row].mv(semiEmbeddedProduct[row],tmp2);
+                TangentVector semiEmbeddedProduct;
+                embeddedHessian[row][col].mv(z,semiEmbeddedProduct);
 
                 for (int subRow=0; subRow<blocksize; subRow++)
-                    this->A_[row][col][subRow][subCol] = tmp2[subRow];
+                    this->A_[row][col][subRow][subCol] = semiEmbeddedProduct[subRow];
             }
 
         }
 
     }
-#endif
 
     //////////////////////////////////////////////////////////////////////////////////////
     //  Further correction due to non-planar configuration space