Skip to content
Snippets Groups Projects
Commit d4352d83 authored by Sander, Oliver's avatar Sander, Oliver
Browse files

Rod3d: Remove the code that measures the solver convergence

It is very unlikely that somebody will use that ever again.
parent db0e3113
No related branches found
No related tags found
Loading
This commit is part of merge request !67. Comments created here will be created in the context of that merge request.
......@@ -14,7 +14,6 @@
#include <dune/gfe/cosseratvtkwriter.hh>
#include <dune/gfe/rigidbodymotion.hh>
#include <dune/gfe/geodesicdifference.hh>
#include <dune/gfe/rotation.hh>
#include <dune/gfe/rodassembler.hh>
#include <dune/gfe/riemanniantrsolver.hh>
......@@ -170,83 +169,6 @@ int main (int argc, char *argv[]) try
BlockVector<FieldVector<double, 6> > strain(x.size()-1);
rodAssembler.getStrain(x,strain);
// If convergence measurement is not desired stop here
if (!instrumented)
exit(0);
// //////////////////////////////////////////////////////////
// Recompute and compare against exact solution
// //////////////////////////////////////////////////////////
SolutionType exactSolution = x;
// //////////////////////////////////////////////////////////
// Compute hessian of the rod functional at the exact solution
// for use of the energy norm it creates.
// //////////////////////////////////////////////////////////
BCRSMatrix<FieldMatrix<double, 6, 6> > hessian;
MatrixIndexSet indices(exactSolution.size(), exactSolution.size());
rodAssembler.getNeighborsPerVertex(indices);
indices.exportIdx(hessian);
BlockVector<FieldVector<double,6> > dummyRhs(x.size());
rodAssembler.assembleGradientAndHessian(exactSolution, dummyRhs, hessian);
double error = std::numeric_limits<double>::max();
SolutionType intermediateSolution(x.size());
// Create statistics file
std::ofstream statisticsFile((resultPath + "trStatistics").c_str());
// Compute error of the initial iterate
typedef BlockVector<FieldVector<double,6> > RodDifferenceType;
RodDifferenceType rodDifference = computeGeodesicDifference(exactSolution, initialIterate);
double oldError = std::sqrt(EnergyNorm<BCRSMatrix<FieldMatrix<double, blocksize, blocksize> >, BlockVector<FieldVector<double,blocksize> > >::normSquared(rodDifference, hessian));
int i;
for (i=0; i<maxTrustRegionSteps; i++) {
// /////////////////////////////////////////////////////
// Read iteration from file
// /////////////////////////////////////////////////////
char iSolFilename[100];
sprintf(iSolFilename, "tmp/intermediateSolution_%04d", i);
FILE* fp = fopen(iSolFilename, "rb");
if (!fp)
DUNE_THROW(IOError, "Couldn't open intermediate solution '" << iSolFilename << "'");
for (size_t j=0; j<intermediateSolution.size(); j++) {
fread(&intermediateSolution[j].r, sizeof(double), 3, fp);
fread(&intermediateSolution[j].q, sizeof(double), 4, fp);
}
fclose(fp);
// /////////////////////////////////////////////////////
// Compute error
// /////////////////////////////////////////////////////
rodDifference = computeGeodesicDifference(exactSolution, intermediateSolution);
error = std::sqrt(EnergyNorm<BCRSMatrix<FieldMatrix<double, blocksize, blocksize> >, BlockVector<FieldVector<double,blocksize> > >::normSquared(rodDifference, hessian));
double convRate = error / oldError;
// Output
std::cout << "Trust-region iteration: " << i << " error : " << error << ", "
<< "convrate " << convRate << std::endl;
statisticsFile << i << " " << error << " " << convRate << std::endl;
if (error < 1e-12)
break;
oldError = error;
}
} catch (Exception& e)
{
std::cout << e.what() << std::endl;
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment