some cleanup of secondDerivativeOfDistanceSquaredWRTSecondArgument

[[Imported from SVN: r7050]]

dune/gfe/unitvector.hh

@@ -160,22 +160,16 @@ public:
      */
     static Dune::FieldMatrix<double,N,N> secondDerivativeOfDistanceSquaredWRTSecondArgument(const UnitVector& p, const UnitVector& q) {
 
-        Dune::FieldMatrix<double,N,N> result;
-
         double sp = p.data_ * q.data_;
-
-        // Compute vector A (see notes)
-        Dune::FieldMatrix<double,1,N> row;
-        row[0] = q.projectOntoTangentSpace(p.globalCoordinates());
-        row *= secondDerivativeOfArcCosSquared(sp);
-
-        Dune::FieldMatrix<double,N,1> column;
-        for (int i=0; i<N; i++)
-            column[i] = p.globalCoordinates()[i] - q.globalCoordinates()[i]*sp;
+        
+        Dune::FieldVector<double,N> pProjected = q.projectOntoTangentSpace(p.globalCoordinates());
 
         Dune::FieldMatrix<double,N,N> A;
-        // A = row * column
-        Dune::FMatrixHelp::multMatrix(column,row,A);
+        for (int i=0; i<N; i++)
+            for (int j=0; j<N; j++)
+                A[i][j] = pProjected[i]*pProjected[j];
+        
+        A *= secondDerivativeOfArcCosSquared(sp);
 
         // Compute matrix B (see notes)
         Dune::FieldMatrix<double,N,N> B;
@@ -184,7 +178,7 @@ public:
                 B[i][j] = (i==j)*sp + p.data_[i]*q.data_[j];
 
         // Bring it all together
-        result = A;
+        Dune::FieldMatrix<double,N,N> result = A;
         result.axpy(-1*derivativeOfArcCosSquared(sp), B);
 
         for (int i=0; i<N; i++)