Skip to content
Snippets Groups Projects
Commit 42e3d743 authored by Oliver Sander's avatar Oliver Sander Committed by sander@PCPOOL.MI.FU-BERLIN.DE
Browse files

remove dead implementation of the Hessian 

[[Imported from SVN: r5947]]
parent 13c58712
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
dune/gfe/localgeodesicfestiffness.hh
+ 0
− 148
View file @ 42e3d743
......@@ -434,8 +434,6 @@ assemble(const Entity& element,
this->A[i][i][j][j] = 1;
#else
#if 1
// ///////////////////////////////////////////////////////////
// Compute gradient by finite-difference approximation
// ///////////////////////////////////////////////////////////
......@@ -477,152 +475,6 @@ assemble(const Entity& element,
for (int k=0; k<blocksize; k++)
this->A[i][j][k] = localSolution[j].projectOntoTangentSpace(this->A[i][j][k]);
#else
double eps = 1e-4;
typedef typename Dune::Matrix<Dune::FieldMatrix<double,blocksize,blocksize> >::row_type::iterator ColumnIterator;
// ///////////////////////////////////////////////////////////
// Compute gradient by finite-difference approximation
// ///////////////////////////////////////////////////////////
std::vector<TargetSpace> forwardSolution = localSolution;
std::vector<TargetSpace> backwardSolution = localSolution;
std::vector<TargetSpace> forwardForwardSolution = localSolution;
std::vector<TargetSpace> forwardBackwardSolution = localSolution;
std::vector<TargetSpace> backwardForwardSolution = localSolution;
std::vector<TargetSpace> backwardBackwardSolution = localSolution;
// ///////////////////////////////////////////////////////////////
// Loop over all blocks of the element matrix
// ///////////////////////////////////////////////////////////////
for (int i=0; i<this->A.N(); i++) {
ColumnIterator cIt = this->A[i].begin();
ColumnIterator cEndIt = this->A[i].end();
for (; cIt!=cEndIt; ++cIt) {
// compute only the upper right triangular matrix
if (cIt.index() < i)
continue;
// ////////////////////////////////////////////////////////////////////////////
// Compute a finite-difference approximation of this hessian matrix block
// ////////////////////////////////////////////////////////////////////////////
for (int j=0; j<blocksize; j++) {
for (int k=0; k<blocksize; k++) {
// compute only the upper right triangular matrix
if (i==cIt.index() && k<j)
continue;
// Diagonal entries
if (i==cIt.index() && j==k) {
infinitesimalVariation(forwardSolution[i], eps, j);
infinitesimalVariation(backwardSolution[i], -eps, j);
double forwardEnergy = energy(element, forwardSolution);
double solutionEnergy = energy(element, localSolution);
double backwardEnergy = energy(element, backwardSolution);
// Second derivative
(*cIt)[j][k] = (forwardEnergy - 2*solutionEnergy + backwardEnergy) / (eps*eps);
forwardSolution[i] = localSolution[i];
backwardSolution[i] = localSolution[i];
} else {
// Off-diagonal entries
infinitesimalVariation(forwardForwardSolution[i], eps, j);
infinitesimalVariation(forwardForwardSolution[cIt.index()], eps, k);
infinitesimalVariation(forwardBackwardSolution[i], eps, j);
infinitesimalVariation(forwardBackwardSolution[cIt.index()], -eps, k);
infinitesimalVariation(backwardForwardSolution[i], -eps, j);
infinitesimalVariation(backwardForwardSolution[cIt.index()], eps, k);
infinitesimalVariation(backwardBackwardSolution[i], -eps, j);
infinitesimalVariation(backwardBackwardSolution[cIt.index()],-eps, k);
double forwardForwardEnergy = energy(element, forwardForwardSolution);
double forwardBackwardEnergy = energy(element, forwardBackwardSolution);
double backwardForwardEnergy = energy(element, backwardForwardSolution);
double backwardBackwardEnergy = energy(element, backwardBackwardSolution);
(*cIt)[j][k] = (forwardForwardEnergy + backwardBackwardEnergy
- forwardBackwardEnergy - backwardForwardEnergy) / (4*eps*eps);
forwardForwardSolution[i] = localSolution[i];
forwardForwardSolution[cIt.index()] = localSolution[cIt.index()];
forwardBackwardSolution[i] = localSolution[i];
forwardBackwardSolution[cIt.index()] = localSolution[cIt.index()];
backwardForwardSolution[i] = localSolution[i];
backwardForwardSolution[cIt.index()] = localSolution[cIt.index()];
backwardBackwardSolution[i] = localSolution[i];
backwardBackwardSolution[cIt.index()] = localSolution[cIt.index()];
}
}
}
}
}
// ///////////////////////////////////////////////////////////////
// Symmetrize the matrix
// This is possible expensive, but I want to be absolute sure
// that the matrix is symmetric.
// ///////////////////////////////////////////////////////////////
for (int i=0; i<this->A.N(); i++) {
ColumnIterator cIt = this->A[i].begin();
ColumnIterator cEndIt = this->A[i].end();
for (; cIt!=cEndIt; ++cIt) {
if (cIt.index()>i)
continue;
if (cIt.index()==i) {
for (int j=1; j<blocksize; j++)
for (int k=0; k<j; k++)
(*cIt)[j][k] = (*cIt)[k][j];
} else {
const Dune::FieldMatrix<double,blocksize,blocksize>& other = this->A[cIt.index()][i];
for (int j=0; j<blocksize; j++)
for (int k=0; k<blocksize; k++)
(*cIt)[j][k] = other[k][j];
}
}
}
#endif
#endif
}
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment