From d4352d83da2a72755f45ea3d34fd490ece70ea1d Mon Sep 17 00:00:00 2001
From: Oliver Sander <oliver.sander@tu-dresden.de>
Date: Thu, 5 Nov 2020 17:14:19 +0100
Subject: [PATCH] Rod3d: Remove the code that measures the solver convergence

It is very unlikely that somebody will use that ever again.
---
 src/rod3d.cc | 78 ----------------------------------------------------
 1 file changed, 78 deletions(-)

diff --git a/src/rod3d.cc b/src/rod3d.cc
index 4f2f86f5..f3e77d74 100644
--- a/src/rod3d.cc
+++ b/src/rod3d.cc
@@ -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;
-- 
GitLab