Newer
Older
#include <config.h>
#include <dune/common/bitsetvector.hh>
#include <dune/common/configparser.hh>
#include <dune/grid/uggrid.hh>
#include <dune/grid/io/file/amirameshreader.hh>
#include <dune-solvers/solvers/iterativesolver.hh>
#include <dune-solvers/norms/energynorm.hh>
#include "src/geodesicdifference.hh"
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
#include "src/rodwriter.hh"
#include "src/rotation.hh"
#include "src/geodesicfeassembler.hh"
#include "src/rodsolver.hh"
// grid dimension
const int dim = 3;
// Image space of the geodesic fe functions
typedef Rotation<3,double> TargetSpace;
// Tangent vector of the image space
const int blocksize = TargetSpace::TangentVector::size;
using namespace Dune;
int main (int argc, char *argv[]) try
{
typedef std::vector<TargetSpace> SolutionType;
// parse data file
ConfigParser parameterSet;
if (argc==2)
parameterSet.parseFile(argv[1]);
else
parameterSet.parseFile("harmonicmaps.parset");
// read solver settings
const int numLevels = parameterSet.get<int>("numLevels");
const double tolerance = parameterSet.get<double>("tolerance");
const int maxTrustRegionSteps = parameterSet.get<int>("maxTrustRegionSteps");
const double initialTrustRegionRadius = parameterSet.get<double>("initialTrustRegionRadius");
const int multigridIterations = parameterSet.get<int>("numIt");
const int nu1 = parameterSet.get<int>("nu1");
const int nu2 = parameterSet.get<int>("nu2");
const int mu = parameterSet.get<int>("mu");
const int baseIterations = parameterSet.get<int>("baseIt");
const double mgTolerance = parameterSet.get<double>("mgTolerance");
const double baseTolerance = parameterSet.get<double>("baseTolerance");
const bool instrumented = parameterSet.get<bool>("instrumented");
std::string resultPath = parameterSet.get("resultPath", "");
// read problem settings
std::string path = parameterSet.get<std::string>("path");
std::string gridFile = parameterSet.get<std::string>("gridFile");
// ///////////////////////////////////////
// Create the grid
// ///////////////////////////////////////
typedef UGGrid<dim> GridType;
GridType grid;
grid.setRefinementType(GridType::COPY);
AmiraMeshReader<GridType>::read(grid, path + gridFile);
grid.globalRefine(numLevels-1);
SolutionType x(grid.size(dim));
// //////////////////////////
// Initial solution
// //////////////////////////
for (int i=0; i<x.size(); i++) {
x[i] = Rotation<3,double>::identity();
}
// backup for error measurement later
SolutionType initialIterate = x;
// /////////////////////////////////////////
// Read Dirichlet values
// /////////////////////////////////////////
BitSetVector<1> allNodes(grid.size(dim));
allNodes.setAll();
LeafBoundaryPatch<GridType> dirichletBoundary(grid, allNodes);
BitSetVector<blocksize> dirichletNodes(grid.size(1));
for (int i=0; i<dirichletNodes.size(); i++)
dirichletNodes[i] = dirichletBoundary.containsVertex(i);
// ////////////////////////////////////////////////////////////
// Create an assembler for the Harmonic Energy Functional
// ////////////////////////////////////////////////////////////
GeodesicFEAssembler<GridType::LeafGridView,TargetSpace> assembler(grid.leafView());
// /////////////////////////////////////////////////
// Create a Riemannian trust-region solver
// /////////////////////////////////////////////////
#if 0
RodSolver<GridType> rodSolver;
rodSolver.setup(grid,
&assembler,
x,
dirichletNodes,
tolerance,
maxTrustRegionSteps,
initialTrustRegionRadius,
multigridIterations,
mgTolerance,
mu, nu1, nu2,
baseIterations,
baseTolerance,
instrumented);
// /////////////////////////////////////////////////////
// Solve!
// /////////////////////////////////////////////////////
std::cout << "Energy: " << assembler.computeEnergy(x) << std::endl;
rodSolver.setInitialSolution(x);
rodSolver.solve();
x = rodSolver.getSol();
#endif
// //////////////////////////////
// Output result
// //////////////////////////////
#if 0
writeRod(x, resultPath + "rod3d.result");
#endif
// //////////////////////////////////////////////////////////
// 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, blocksize, blocksize> > hessian;
// MatrixIndexSet indices(exactSolution.size(), exactSolution.size());
// assembler.getNeighborsPerVertex(indices);
// indices.exportIdx(hessian);
assembler.assembleMatrix(exactSolution, 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,blocksize> > DifferenceType;
DifferenceType geodesicDifference = computeGeodesicDifference(exactSolution, initialIterate);
double oldError = std::sqrt(EnergyNorm<BCRSMatrix<FieldMatrix<double, blocksize, blocksize> >, BlockVector<FieldVector<double,blocksize> > >::normSquared(geodesicDifference, 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 (int j=0; j<intermediateSolution.size(); j++) {
fread(&intermediateSolution[j], sizeof(double), 4, fp);
}
fclose(fp);
// /////////////////////////////////////////////////////
// Compute error
// /////////////////////////////////////////////////////
geodesicDifference = computeGeodesicDifference(exactSolution, intermediateSolution);
error = std::sqrt(EnergyNorm<BCRSMatrix<FieldMatrix<double, blocksize, blocksize> >, BlockVector<FieldVector<double,blocksize> > >::normSquared(geodesicDifference, 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 << std::endl;
}