diff --git a/rod3d.cc b/rod3d.cc
index 225dea78a0a017eff24d2086794a01e8109114c4..0d907b4fccf53dae57e73b54844d1c84cf695aba 100644
--- a/rod3d.cc
+++ b/rod3d.cc
@@ -1,24 +1,24 @@
#include <config.h>
-//#define DUNE_EXPRESSIONTEMPLATES
+#include <dune/common/bitfield.hh>
+#include <dune/common/configparser.hh>
+
#include <dune/grid/onedgrid.hh>
#include <dune/istl/io.hh>
-#include <dune/common/bitfield.hh>
-#include "src/quaternion.hh"
-
-#include "src/rodassembler.hh"
-#include "src/rodsolver.hh"
#include "../solver/iterativesolver.hh"
#include "../common/geomestimator.hh"
#include "../common/energynorm.hh"
-#include <dune/common/configparser.hh>
#include "src/configuration.hh"
+#include "src/roddifference.hh"
#include "src/rodwriter.hh"
+#include "src/quaternion.hh"
+#include "src/rodassembler.hh"
+#include "src/rodsolver.hh"
// Number of degrees of freedom:
// 7 (x, y, z, q_1, q_2, q_3, q_4) for a spatial rod
@@ -27,6 +27,15 @@ const int blocksize = 6;
using namespace Dune;
using std::string;
+double computeEnergyNormSquared(const BlockVector<FieldVector<double,6> >& x,
+ const BCRSMatrix<FieldMatrix<double, 6, 6> >& matrix)
+{
+ BlockVector<FieldVector<double, 6> > tmp(x.size());
+ tmp = 0;
+ matrix.umv(x,tmp);
+ return x*tmp;
+}
+
int main (int argc, char *argv[]) try
{
typedef std::vector<Configuration> SolutionType;
@@ -47,6 +56,7 @@ int main (int argc, char *argv[]) try
const double baseTolerance = parameterSet.get("baseTolerance", double(0));
const double initialTrustRegionRadius = parameterSet.get("initialTrustRegionRadius", double(1));
const int numRodBaseElements = parameterSet.get("numRodBaseElements", int(0));
+ const bool instrumented = parameterSet.get("instrumented", int(0));
// ///////////////////////////////////////
// Create the grid
@@ -63,6 +73,8 @@ int main (int argc, char *argv[]) try
// //////////////////////////
// Initial solution
// //////////////////////////
+ FieldVector<double,3> zAxis(0);
+ zAxis[2] = 1;
for (int i=0; i<x.size(); i++) {
x[i].r[0] = 0; // x
@@ -70,15 +82,18 @@ int main (int argc, char *argv[]) try
x[i].r[2] = double(i)/(x.size()-1); // z
//x[i].r[2] = i+5;
x[i].q = Quaternion<double>::identity();
+ //x[i].q = Quaternion<double>(zAxis, M_PI/2 * double(i)/(x.size()-1));
}
-// x[x.size()-1].r[0] = 0;
-// x[x.size()-1].r[1] = 0;
-// x[x.size()-1].r[2] = 0;
+
#if 1
- FieldVector<double,3> zAxis(0);
- zAxis[2] = 1;
- x[x.size()-1].q = Quaternion<double>(zAxis, M_PI);
+ FieldVector<double,3> xAxis(0);
+ xAxis[0] = 1;
+ x[1].r[2] = 0.25;
+ x.back().r[2] = 0.5;
+ x[0].q = Quaternion<double>(xAxis, -M_PI/2);
+ x.back().q = Quaternion<double>(xAxis, M_PI/2);
+
#endif
std::cout << "Left boundary orientation:" << std::endl;
@@ -110,7 +125,8 @@ int main (int argc, char *argv[]) try
// ///////////////////////////////////////////
RodAssembler<GridType> rodAssembler(grid);
rodAssembler.setShapeAndMaterial(0.01, 0.0001, 0.0001, 2.5e5, 0.3);
-
+ //rodAssembler.setParameters(0,0,0,1,1,1);
+
RodSolver<GridType> rodSolver;
rodSolver.setup(grid,
&rodAssembler,
@@ -121,7 +137,8 @@ int main (int argc, char *argv[]) try
mgTolerance,
mu, nu1, nu2,
baseIterations,
- baseTolerance);
+ baseTolerance,
+ instrumented);
// /////////////////////////////////////////////////////
// Solve!
@@ -142,9 +159,92 @@ int main (int argc, char *argv[]) try
rodAssembler.getStrain(x,strain);
//std::cout << strain << std::endl;
//exit(0);
+
+ //writeRod(x, "rod3d.strain");
+
+ // //////////////////////////////////////////////////////////
+ // Recompute and compare against exact solution
+ // //////////////////////////////////////////////////////////
- writeRod(x, strain, "rod3d.strain");
+ 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);
+ rodAssembler.assembleMatrix(exactSolution, hessian);
+
+ double error = std::numeric_limits<double>::max();
+ double oldError = 0;
+ double totalConvRate = 1;
+
+ SolutionType intermediateSolution(x.size());
+
+ // Compute error of the initial 3d solution
+
+ // This should really be exactSol-initialSol, but we're starting
+ // from zero anyways
+ //oldError += computeEnergyNormSquared(exactSol3d, *hessian3d);
+
+ /** \todo Rod error still missing */
+
+ oldError = std::sqrt(oldError);
+
+
+
+ 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 (int 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
+ // /////////////////////////////////////////////////////
+ typedef BlockVector<FieldVector<double,6> > RodDifferenceType;
+ RodDifferenceType rodDifference = computeRodDifference(exactSolution, intermediateSolution);
+
+ error = std::sqrt(computeEnergyNormSquared(rodDifference, hessian));
+
+
+ double convRate = error / oldError;
+ totalConvRate *= convRate;
+
+ // Output
+ std::cout << "Trust-region iteration: " << i << " error : " << error << ", "
+ << "convrate " << convRate
+ << " total conv rate " << std::pow(totalConvRate, 1/((double)i+1)) << std::endl;
+
+ if (error < 1e-12)
+ break;
+
+ oldError = error;
+
+ }
+
+
+ // //////////////////////////////
} catch (Exception e) {
std::cout << e << std::endl;