From da2e2a7937e4922dbf3c9b5bc82fc774ffa6c50e Mon Sep 17 00:00:00 2001
From: Oliver Sander <sander@igpm.rwth-aachen.de>
Date: Wed, 9 Jul 2014 11:51:25 +0000
Subject: [PATCH] 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 | 69 ++++++-----------------
1 file changed, 16 insertions(+), 53 deletions(-)
diff --git a/dune/gfe/localgeodesicfeadolcstiffness.hh b/dune/gfe/localgeodesicfeadolcstiffness.hh
index 3bbbcdd8..819ac2f2 100644
--- a/dune/gfe/localgeodesicfeadolcstiffness.hh
+++ b/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
--
GitLab