diff --git a/Makefile.am b/Makefile.am
index cb1fd9ff4ff6d3da55673752406282d8ae430d82..ea8fba908986b09fb6ab846d92619a5f4ff5f383 100644
--- a/Makefile.am
+++ b/Makefile.am
@@ -9,7 +9,7 @@ SUBDIRS = m4 src test
LDADD = $(IPOPT_LDFLAGS) $(IPOPT_LIBS)
AM_CPPFLAGS += $(IPOPT_CPPFLAGS)
-noinst_PROGRAMS = rodobstacle rod3d harmonicmaps dirneucoupling neudircoupling rod-eoc
+noinst_PROGRAMS = rodobstacle rod3d harmonicmaps dirneucoupling rod-eoc
rodobstacle_SOURCES = rodobstacle.cc
rod3d_SOURCES = rod3d.cc
@@ -26,11 +26,6 @@ dirneucoupling_CXXFLAGS = $(UG_CPPFLAGS) $(AMIRAMESH_CPPFLAGS) $(IPOPT_CPPFLAGS)
dirneucoupling_LDADD = $(UG_LDFLAGS) $(AMIRAMESH_LDFLAGS) $(UG_LIBS) $(AMIRAMESH_LIBS) \
$(IPOPT_LDFLAGS) $(IPOPT_LIBS) $(PSURFACE_LDFLAGS) $(PSURFACE_LIBS)
-neudircoupling_SOURCES = neudircoupling.cc
-neudircoupling_CXXFLAGS = $(UG_CPPFLAGS) $(AMIRAMESH_CPPFLAGS) $(IPOPT_CPPFLAGS) $(PSURFACE_CPPFLAGS)
-neudircoupling_LDADD = $(UG_LDFLAGS) $(AMIRAMESH_LDFLAGS) $(UG_LIBS) $(AMIRAMESH_LIBS) \
- $(IPOPT_LDFLAGS) $(IPOPT_LIBS) $(PSURFACE_LDFLAGS) $(PSURFACE_LIBS)
-
# don't follow the full GNU-standard
# we need automake 1.5
AUTOMAKE_OPTIONS = foreign 1.5
diff --git a/neudircoupling.cc b/neudircoupling.cc
deleted file mode 100644
index f36db04986b3af12b54ba2daae25552f96066cc2..0000000000000000000000000000000000000000
--- a/neudircoupling.cc
+++ /dev/null
@@ -1,743 +0,0 @@
-#include <config.h>
-
-#include <dune/common/bitsetvector.hh>
-#include <dune/common/configparser.hh>
-
-#include <dune/grid/onedgrid.hh>
-#include <dune/grid/uggrid.hh>
-#include <dune/grid/io/file/amirameshreader.hh>
-#include <dune/grid/io/file/amirameshwriter.hh>
-
-#include <dune/istl/solvers.hh>
-
-#include <dune/solvers/iterationsteps/multigridstep.hh>
-#include <dune/solvers/solvers/loopsolver.hh>
-#include <dune/solvers/iterationsteps/projectedblockgsstep.hh>
-#ifdef HAVE_IPOPT
-#include <dune/solvers/solvers/quadraticipopt.hh>
-#endif
-#include <dune/ag-common/readbitfield.hh>
-#include <dune/solvers/norms/energynorm.hh>
-#include <dune/ag-common/boundarypatch.hh>
-#include <dune/ag-common/prolongboundarypatch.hh>
-#include <dune/ag-common/sampleonbitfield.hh>
-#include <dune/ag-common/neumannassembler.hh>
-#include <dune/ag-common/computestress.hh>
-
-#include <dune/ag-common/functionspacebases/q1nodalbasis.hh>
-#include <dune/ag-common/assemblers/operatorassembler.hh>
-#include <dune/ag-common/assemblers/localassemblers/stvenantkirchhoffassembler.hh>
-
-#include "src/quaternion.hh"
-#include "src/rodassembler.hh"
-#include "src/rigidbodymotion.hh"
-#include "src/averageinterface.hh"
-#include "src/riemanniantrsolver.hh"
-#include "src/geodesicdifference.hh"
-#include "src/rodwriter.hh"
-#include "src/makestraightrod.hh"
-
-// Space dimension
-const int dim = 3;
-
-using namespace Dune;
-using std::string;
-using std::vector;
-
-// Some types that I need
-typedef vector<RigidBodyMotion<dim> > RodSolutionType;
-typedef BlockVector<FieldVector<double, 6> > RodDifferenceType;
-
-
-int main (int argc, char *argv[]) try
-{
- // Some types that I need
- typedef BCRSMatrix<FieldMatrix<double, dim, dim> > MatrixType;
- typedef BlockVector<FieldVector<double, dim> > VectorType;
-
- // parse data file
- ConfigParser parameterSet;
- if (argc==2)
- parameterSet.parseFile(argv[1]);
- else
- parameterSet.parseFile("neudircoupling.parset");
-
- // read solver settings
- const int numLevels = parameterSet.get<int>("numLevels");
- const double ddTolerance = parameterSet.get<double>("ddTolerance");
- const int maxDirichletNeumannSteps = parameterSet.get<int>("maxDirichletNeumannSteps");
- const double trTolerance = parameterSet.get<double>("trTolerance");
- const int maxTrustRegionSteps = parameterSet.get<int>("maxTrustRegionSteps");
- const int trVerbosity = parameterSet.get<int>("trVerbosity");
- 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 double initialTrustRegionRadius = parameterSet.get<double>("initialTrustRegionRadius");
- const double damping = parameterSet.get<double>("damping");
- string resultPath = parameterSet.get("resultPath", "");
-
- // Problem settings
- string path = parameterSet.get<string>("path");
- string objectName = parameterSet.get<string>("gridFile");
- string dirichletNodesFile = parameterSet.get<string>("dirichletNodes");
- string dirichletValuesFile = parameterSet.get<string>("dirichletValues");
- string interfaceNodesFile = parameterSet.get<string>("interfaceNodes");
- const int numRodBaseElements = parameterSet.get<int>("numRodBaseElements");
-
- double E = parameterSet.get<double>("E");
- double nu = parameterSet.get<double>("nu");
-
- // rod material parameters
- double rodA = parameterSet.get<double>("rodA");
- double rodJ1 = parameterSet.get<double>("rodJ1");
- double rodJ2 = parameterSet.get<double>("rodJ2");
- double rodE = parameterSet.get<double>("rodE");
- double rodNu = parameterSet.get<double>("rodNu");
-
- Dune::array<FieldVector<double,3>,2> rodRestEndPoint;
- rodRestEndPoint[0][0] = parameterSet.get<double>("rodRestEndPoint0X");
- rodRestEndPoint[0][1] = parameterSet.get<double>("rodRestEndPoint0Y");
- rodRestEndPoint[0][2] = parameterSet.get<double>("rodRestEndPoint0Z");
- rodRestEndPoint[1][0] = parameterSet.get<double>("rodRestEndPoint1X");
- rodRestEndPoint[1][1] = parameterSet.get<double>("rodRestEndPoint1Y");
- rodRestEndPoint[1][2] = parameterSet.get<double>("rodRestEndPoint1Z");
-
- // ///////////////////////////////////////
- // Create the rod grid
- // ///////////////////////////////////////
- typedef OneDGrid RodGridType;
- RodGridType rodGrid(numRodBaseElements, 0, (rodRestEndPoint[1]-rodRestEndPoint[0]).two_norm());
-
- // ///////////////////////////////////////
- // Create the grid for the 3d object
- // ///////////////////////////////////////
- typedef UGGrid<dim> GridType;
- GridType grid;
- grid.setRefinementType(GridType::COPY);
-
- AmiraMeshReader<GridType>::read(grid, path + objectName);
-
- // //////////////////////////////////////
- // Globally refine grids
- // //////////////////////////////////////
-
- rodGrid.globalRefine(numLevels-1);
- grid.globalRefine(numLevels-1);
-
- RodSolutionType rodX(rodGrid.size(1));
-
- // //////////////////////////
- // Initial solution
- // //////////////////////////
-
- makeStraightRod(rodX, rodGrid.size(1), rodRestEndPoint[0], rodRestEndPoint[1]);
-
- // /////////////////////////////////////////
- // Read Dirichlet values
- // /////////////////////////////////////////
- rodX.back().r[0] = parameterSet.get("dirichletValueX", rodRestEndPoint[1][0]);
- rodX.back().r[1] = parameterSet.get("dirichletValueY", rodRestEndPoint[1][1]);
- rodX.back().r[2] = parameterSet.get("dirichletValueZ", rodRestEndPoint[1][2]);
-
- FieldVector<double,3> axis;
- axis[0] = parameterSet.get("dirichletAxisX", double(0));
- axis[1] = parameterSet.get("dirichletAxisY", double(0));
- axis[2] = parameterSet.get("dirichletAxisZ", double(0));
- double angle = parameterSet.get("dirichletAngle", double(0));
-
- rodX.back().q = Rotation<3,double>(axis, M_PI*angle/180);
-
- // Backup initial rod iterate for later reference
- RodSolutionType initialIterateRod = rodX;
-
- int toplevel = rodGrid.maxLevel();
-
- // /////////////////////////////////////////////////////
- // Determine the Dirichlet nodes
- // /////////////////////////////////////////////////////
- VectorType coarseDirichletValues(grid.size(0, dim));
- AmiraMeshReader<int>::readFunction(coarseDirichletValues, path + dirichletValuesFile);
-
- LevelBoundaryPatch<GridType> coarseDirichletBoundary(grid,0);
- readBoundaryPatch(coarseDirichletBoundary, path + dirichletNodesFile);
-
- LeafBoundaryPatch<GridType> dirichletBoundary;
- PatchProlongator<GridType>::prolong(coarseDirichletBoundary, dirichletBoundary);
-
- BitSetVector<dim> dirichletNodes(grid.size(dim));
- for (int i=0; i<dirichletNodes.size(); i++)
- dirichletNodes[i] = dirichletBoundary.containsVertex(i);
-
- VectorType dirichletValues;
- sampleOnBitField(grid, coarseDirichletValues, dirichletValues, dirichletNodes);
-
- // /////////////////////////////////////////////////////
- // Determine the interface boundary
- // /////////////////////////////////////////////////////
- std::vector<LevelBoundaryPatch<GridType> > interfaceBoundary;
- interfaceBoundary.resize(numLevels);
- interfaceBoundary[0].setup(grid, 0);
- readBoundaryPatch(interfaceBoundary[0], path + interfaceNodesFile);
- PatchProlongator<GridType>::prolong(interfaceBoundary);
-
- // //////////////////////////////////////////
- // Assemble 3d linear elasticity problem
- // //////////////////////////////////////////
-
- typedef Q1NodalBasis<GridType::LeafGridView,double> FEBasis;
- FEBasis basis(grid.leafView());
- OperatorAssembler<FEBasis,FEBasis> assembler(basis, basis);
-
- StVenantKirchhoffAssembler<GridType, FEBasis::LocalFiniteElement, FEBasis::LocalFiniteElement> localAssembler(E, nu);
- MatrixType stiffnessMatrix3d;
-
- assembler.assemble(localAssembler, stiffnessMatrix3d);
-
- // ///////////////////////////////////////////////////////////////////////
- // Assemble the mass matrix of the interface boundary.
- // It is needed to compute the strong normal stresses resulting from
- // the Dirichlet boundary conditions.
- // ///////////////////////////////////////////////////////////////////////
-
- MatrixType surfaceMassMatrix;
- assembleSurfaceMassMatrix<GridType::LevelGridView,dim>(interfaceBoundary.back(), surfaceMassMatrix);
-
- std::vector<int> globalToLocal;
- interfaceBoundary.back().makeGlobalToLocal(globalToLocal);
-
- // ////////////////////////////////////////////////////////////
- // Create solution and rhs vectors
- // ////////////////////////////////////////////////////////////
-
- VectorType x3d(grid.size(toplevel,dim));
- VectorType rhs3d(grid.size(toplevel,dim));
-
- // No external forces
- rhs3d = 0;
-
- // Set initial solution
- x3d = 0;
- for (int i=0; i<x3d.size(); i++)
- for (int j=0; j<dim; j++)
- if (dirichletNodes[i][j])
- x3d[i][j] = dirichletValues[i][j];
-
- // ///////////////////////////////////////////////////////////////////
- // Add the interface boundary nodes to the set of Dirichlet nodes
- // ///////////////////////////////////////////////////////////////////
-
- for (int i=0; i<dirichletNodes.size(); i++)
- for (int j=0; j<dim; j++)
- dirichletNodes[i][j] = dirichletNodes[i][j] || interfaceBoundary.back().containsVertex(i);
-
- // ///////////////////////////////////////////
- // Dirichlet nodes for the rod problem
- // ///////////////////////////////////////////
-
- BitSetVector<6> rodDirichletNodes(rodGrid.size(1));
- rodDirichletNodes.unsetAll();
-
- //rodDirichletNodes[0] = true;
- rodDirichletNodes.back() = true;
-
- // ///////////////////////////////////////////
- // Create a solver for the rod problem
- // ///////////////////////////////////////////
-
- RodLocalStiffness<RodGridType::LeafGridView,double> rodLocalStiffness(rodGrid.leafView(),
- rodA, rodJ1, rodJ2, rodE, rodNu);
-
- RodAssembler<RodGridType::LeafGridView> rodAssembler(rodGrid.leafView(), &rodLocalStiffness);
-
- RiemannianTrustRegionSolver<RodGridType,RigidBodyMotion<3> > rodSolver;
- rodSolver.setup(rodGrid,
- &rodAssembler,
- rodX,
- rodDirichletNodes,
- trTolerance,
- maxTrustRegionSteps,
- initialTrustRegionRadius,
- multigridIterations,
- mgTolerance,
- mu, nu1, nu2,
- baseIterations,
- baseTolerance,
- false);
-
- switch (trVerbosity) {
- case 0:
- rodSolver.verbosity_ = Solver::QUIET; break;
- case 1:
- rodSolver.verbosity_ = Solver::REDUCED; break;
- default:
- rodSolver.verbosity_ = Solver::FULL; break;
- }
-
-
- // ////////////////////////////////
- // Create a multigrid solver
- // ////////////////////////////////
-
- // First create a gauss-seidel base solver
-#ifdef HAVE_IPOPT
- QuadraticIPOptSolver<MatrixType,VectorType> baseSolver;
-#endif
- baseSolver.verbosity_ = NumProc::QUIET;
- baseSolver.tolerance_ = baseTolerance;
-
- // Make pre and postsmoothers
- BlockGSStep<MatrixType, VectorType> presmoother, postsmoother;
-
- MultigridStep<MatrixType, VectorType> multigridStep(stiffnessMatrix3d, x3d, rhs3d, 1);
-
- multigridStep.setMGType(mu, nu1, nu2);
- multigridStep.ignoreNodes_ = &dirichletNodes;
- multigridStep.basesolver_ = &baseSolver;
- multigridStep.setSmoother(&presmoother, &postsmoother);
- multigridStep.verbosity_ = Solver::QUIET;
-
-
-
- EnergyNorm<MatrixType, VectorType> energyNorm(multigridStep);
-
- ::LoopSolver<VectorType> solver(&multigridStep,
- // IPOpt doesn't like to be started in the solution
- (numLevels!=1) ? multigridIterations : 1,
- mgTolerance,
- &energyNorm,
- Solver::QUIET);
-
- // ////////////////////////////////////
- // Create the transfer operators
- // ////////////////////////////////////
- for (int k=0; k<multigridStep.mgTransfer_.size(); k++)
- delete(multigridStep.mgTransfer_[k]);
-
- multigridStep.mgTransfer_.resize(toplevel);
-
- for (int i=0; i<multigridStep.mgTransfer_.size(); i++){
- CompressedMultigridTransfer<VectorType>* newTransferOp = new CompressedMultigridTransfer<VectorType>;
- newTransferOp->setup(grid,i,i+1);
- multigridStep.mgTransfer_[i] = newTransferOp;
- }
-
- // /////////////////////////////////////////////////////
- // Dirichlet-Neumann Solver
- // /////////////////////////////////////////////////////
-
- // Init interface value
- RigidBodyMotion<3> referenceInterface = rodX[0];
- //RigidBodyMotion<3> lambda = referenceInterface;
- FieldVector<double,3> lambdaForce(0);
- FieldVector<double,3> lambdaTorque(0);
-
- lambdaForce[2] = -5000;
-
- //
- double normOfOldCorrection = 0;
- int dnStepsActuallyTaken = 0;
- for (int i=0; i<maxDirichletNeumannSteps; i++) {
-
- std::cout << "----------------------------------------------------" << std::endl;
- std::cout << " Dirichlet-Neumann Step Number: " << i << std::endl;
- std::cout << "----------------------------------------------------" << std::endl;
-
- // Backup of the current solution for the error computation later on
- VectorType oldSolution3d = x3d;
- RodSolutionType oldSolutionRod = rodX;
-
- // //////////////////////////////////////////////////
- // Neumann step for the rod
- // //////////////////////////////////////////////////
-
- rodSolver.setInitialSolution(rodX);
- rodAssembler.setNeumannData(lambdaForce, lambdaTorque, FieldVector<double,3>(0), FieldVector<double,3>(0));
- rodSolver.solve();
-
- rodX = rodSolver.getSol();
-
-// for (int j=0; j<rodX.size(); j++)
-// std::cout << rodX[j] << std::endl;
-
- // Get resultant force, just for checking
- BitSetVector<1> couplingBitfield(rodX.size(),false);
- couplingBitfield[0] = true;
- LeafBoundaryPatch<RodGridType> couplingBoundary(rodGrid, couplingBitfield);
-
- FieldVector<double,dim> resultantForceDebug, resultantTorqueDebug;
- resultantForceDebug = rodAssembler.getResultantForce(couplingBoundary, rodX, resultantTorqueDebug);
-
- // Flip orientation
- resultantForceDebug *= -1;
- resultantTorqueDebug *= -1;
-
- std::cout << "debugging: resultant force: " << resultantForceDebug
- << " norm: " << resultantForceDebug.two_norm() << std::endl;
- std::cout << "debugging: resultant torque: " << resultantTorqueDebug
- << " norm: " << resultantTorqueDebug.two_norm() << std::endl;
-
-
- // ///////////////////////////////////////////////////////////
- // Extract Dirichlet values and transfer it to the 3d object
- // ///////////////////////////////////////////////////////////
-
- // Using that index 0 is always the left boundary for a uniformly refined OneDGrid
- RigidBodyMotion<3> resultantConfiguration = rodX[0];
-
- std::cout << "Resultant configuration: " << resultantConfiguration << std::endl;
-
- // Compute difference to the reference interface
- /** \todo This is a group operation --> put it into the RigidBodyMotion class */
- RigidBodyMotion<3> differenceToReferenceInterface = referenceInterface;
- differenceToReferenceInterface.q.invert();
- differenceToReferenceInterface.r *= -1;
- differenceToReferenceInterface.q.mult(resultantConfiguration.q);
- differenceToReferenceInterface.r += resultantConfiguration.r;
-
-
-
- GridType::Codim<dim>::LeafIterator vIt = grid.leafbegin<dim>();
- GridType::Codim<dim>::LeafIterator vEndIt = grid.leafend<dim>();
-
- for (; vIt!=vEndIt; ++vIt) {
-
- unsigned int idx = grid.leafIndexSet().index(*vIt);
-
- // Consider only vertices on the interface boundary
- if (interfaceBoundary.back().containsVertex(idx)) {
-
- // apply the rigid body motion to the vertex position and subtract the old position
- FieldMatrix<double,3,3> rotationMatrix;
- differenceToReferenceInterface.q.matrix(rotationMatrix);
- rotationMatrix.mv(vIt->geometry().corner(0), x3d[idx]);
- x3d[idx] += differenceToReferenceInterface.r;
- x3d[idx] -= vIt->geometry().corner(0);
-
- }
-
- }
-
- // ///////////////////////////////////////////////////////////
- // Solve the Dirichlet problem for the 3d body
- // ///////////////////////////////////////////////////////////
- multigridStep.setProblem(stiffnessMatrix3d, x3d, rhs3d, grid.maxLevel()+1);
-
- solver.preprocess();
- multigridStep.preprocess();
-
- solver.solve();
-
- x3d = multigridStep.getSol();
-
- // ///////////////////////////////////////////////////////////
- // Extract new interface resultant force and torque
- // ///////////////////////////////////////////////////////////
-
- FieldVector<double,3> resultantForce(0);
- FieldVector<double,3> resultantTorque(0);
-
- // the weak normal stress, or, in other words, the residual
- VectorType weakNormalStress = rhs3d;
- stiffnessMatrix3d.mmv(x3d, weakNormalStress);
-
- // consider only the coefficients on the interface boundary
- VectorType localWeakNormalStress(interfaceBoundary.back().numVertices());
- for (int j=0; j<globalToLocal.size(); j++)
- if (globalToLocal[j] != -1)
- localWeakNormalStress[globalToLocal[j]] = weakNormalStress[j];
-
- // Compute the strong normal stress, which is the weak stress divided by the surface mass matrix
- VectorType localStrongNormalStress = localWeakNormalStress; // initial value
-
- // Make small cg solver
- MatrixAdapter<MatrixType,VectorType,VectorType> op(surfaceMassMatrix);
- SeqILU0<MatrixType,VectorType,VectorType> ilu0(surfaceMassMatrix,1.0);
- CGSolver<VectorType> cgsolver(op,ilu0,1E-4,100,0);
- Dune::InverseOperatorResult statistics;
- cgsolver.apply(localStrongNormalStress, localWeakNormalStress, statistics);
-
- VectorType strongNormalStress(weakNormalStress.size());
- strongNormalStress = 0;
- for (int j=0; j<globalToLocal.size(); j++)
- if (globalToLocal[j] != -1)
- strongNormalStress[j] = localStrongNormalStress[globalToLocal[j]];
-
-
- computeTotalForceAndTorque(interfaceBoundary.back(), strongNormalStress, resultantConfiguration.r,
- resultantForce, resultantTorque);
-
- std::cout << "average force: " << resultantForce << std::endl;
- std::cout << "average torque: " << resultantTorque << std::endl;
-
- // ///////////////////////////////////////////////////////////
- // Compute new damped interface value
- // ///////////////////////////////////////////////////////////
- for (int j=0; j<dim; j++) {
- lambdaForce[j] = (1-damping) * lambdaForce[j] + damping * resultantForce[j];
- lambdaTorque[j] = (1-damping) * lambdaTorque[j] + damping * resultantTorque[j];
- }
-
- std::cout << "Lambda force: " << lambdaForce << std::endl;
- std::cout << "Lambda torque: " << lambdaTorque << std::endl;
-
- // ////////////////////////////////////////////////////////////////////////
- // Write the two iterates to disk for later convergence rate measurement
- // ////////////////////////////////////////////////////////////////////////
-
- // First the 3d body
- std::stringstream iAsAscii;
- iAsAscii << i;
- std::string iSolFilename = resultPath + "tmp/intermediate3dSolution_" + iAsAscii.str();
-
- LeafAmiraMeshWriter<GridType> amiraMeshWriter;
- amiraMeshWriter.addVertexData(x3d, grid.leafView());
- amiraMeshWriter.write(iSolFilename);
-
- // Then the rod
- iSolFilename = resultPath + "tmp/intermediateRodSolution_" + iAsAscii.str();
-
- FILE* fpRod = fopen(iSolFilename.c_str(), "wb");
- if (!fpRod)
- DUNE_THROW(SolverError, "Couldn't open file " << iSolFilename << " for writing");
-
- for (int j=0; j<rodX.size(); j++) {
-
- for (int k=0; k<dim; k++)
- fwrite(&rodX[j].r[k], sizeof(double), 1, fpRod);
-
- for (int k=0; k<4; k++) // 3d hardwired here!
- fwrite(&rodX[j].q[k], sizeof(double), 1, fpRod);
-
- }
-
- fclose(fpRod);
-
- // ////////////////////////////////////////////
- // Compute error in the energy norm
- // ////////////////////////////////////////////
-
- // the 3d body
- double oldNorm = EnergyNorm<MatrixType,VectorType>::normSquared(oldSolution3d, stiffnessMatrix3d);
- oldSolution3d -= x3d;
- double normOfCorrection = EnergyNorm<MatrixType,VectorType>::normSquared(oldSolution3d, stiffnessMatrix3d);
-
- double max3dRelCorrection = 0;
- for (size_t j=0; j<x3d.size(); j++)
- for (int k=0; k<dim; k++)
- max3dRelCorrection = std::max(max3dRelCorrection,
- std::fabs(oldSolution3d[j][k])/ std::fabs(x3d[j][k]));
-
- // the rod
- RodDifferenceType rodDiff = computeGeodesicDifference(oldSolutionRod, rodX);
- double maxRodRelCorrection = 0;
- for (size_t j=0; j<rodX.size(); j++)
- for (int k=0; k<dim; k++)
- maxRodRelCorrection = std::max(maxRodRelCorrection,
- std::fabs(rodDiff[j][k])/ std::fabs(rodX[j].r[k]));
-
- // Absolute corrections
- double maxRodCorrection = computeGeodesicDifference(oldSolutionRod, rodX).infinity_norm();
- double max3dCorrection = oldSolution3d.infinity_norm();
-
-
- std::cout << "rod correction: " << maxRodCorrection
- << " rod rel correction: " << maxRodRelCorrection
- << " 3d correction: " << max3dCorrection
- << " 3d rel correction: " << max3dRelCorrection << std::endl;
-
-
- oldNorm = std::sqrt(oldNorm);
-
- normOfCorrection = std::sqrt(normOfCorrection);
-
- double relativeError = normOfCorrection / oldNorm;
-
- double convRate = normOfCorrection / normOfOldCorrection;
-
- normOfOldCorrection = normOfCorrection;
-
- // Output
- std::cout << "DD iteration: " << i << " -- ||u^{n+1} - u^n|| / ||u^n||: " << relativeError << ", "
- << "convrate " << convRate << "\n";
-
- dnStepsActuallyTaken = i;
-
- //if (relativeError < ddTolerance)
- if (std::max(max3dRelCorrection,maxRodRelCorrection) < ddTolerance)
- break;
-
- }
-
-
-
- // //////////////////////////////////////////////////////////
- // Recompute and compare against exact solution
- // //////////////////////////////////////////////////////////
- VectorType exactSol3d = x3d;
- RodSolutionType exactSolRod = rodX;
-
- // //////////////////////////////////////////////////////////
- // Compute hessian of the rod functional at the exact solution
- // for use of the energy norm it creates.
- // //////////////////////////////////////////////////////////
-
- BCRSMatrix<FieldMatrix<double, 6, 6> > hessianRod;
- MatrixIndexSet indices(exactSolRod.size(), exactSolRod.size());
- rodAssembler.getNeighborsPerVertex(indices);
- indices.exportIdx(hessianRod);
- rodAssembler.assembleMatrix(exactSolRod, hessianRod);
-
-
- double error = std::numeric_limits<double>::max();
- double oldError = 0;
-
- VectorType intermediateSol3d(x3d.size());
- RodSolutionType intermediateSolRod(rodX.size());
-
- // Compute error of the initial 3d solution
-
- // This should really be exactSol-initialSol, but we're starting
- // from zero anyways
- oldError += EnergyNorm<MatrixType,VectorType>::normSquared(exactSol3d, stiffnessMatrix3d);
-
- // Error of the initial rod iterate
- RodDifferenceType rodDifference = computeGeodesicDifference(initialIterateRod, exactSolRod);
- oldError += EnergyNorm<BCRSMatrix<FieldMatrix<double,6,6> >,RodDifferenceType>::normSquared(rodDifference, hessianRod);
-
- oldError = std::sqrt(oldError);
-
- // Store the history of total conv rates so we can filter out numerical
- // dirt in the end.
- std::vector<double> totalConvRate(maxDirichletNeumannSteps+1);
- totalConvRate[0] = 1;
-
-
- double oldConvRate = 100;
- bool firstTime = true;
- std::stringstream levelAsAscii, dampingAsAscii;
- levelAsAscii << numLevels;
- dampingAsAscii << damping;
- std::string filename = resultPath + "convrate_" + levelAsAscii.str() + "_" + dampingAsAscii.str();
-
- int i;
- for (i=0; i<dnStepsActuallyTaken; i++) {
-
- // /////////////////////////////////////////////////////
- // Read iteration from file
- // /////////////////////////////////////////////////////
-
- // Read 3d solution from file
- std::stringstream iAsAscii;
- iAsAscii << i;
- std::string iSolFilename = resultPath + "tmp/intermediate3dSolution_" + iAsAscii.str();
-
- AmiraMeshReader<int>::readFunction(intermediateSol3d, iSolFilename);
-
- // Read rod solution from file
- iSolFilename = resultPath + "tmp/intermediateRodSolution_" + iAsAscii.str();
-
- FILE* fpInt = fopen(iSolFilename.c_str(), "rb");
- if (!fpInt)
- DUNE_THROW(IOError, "Couldn't open intermediate solution '" << iSolFilename << "'");
- for (int j=0; j<intermediateSolRod.size(); j++) {
- fread(&intermediateSolRod[j].r, sizeof(double), dim, fpInt);
- fread(&intermediateSolRod[j].q, sizeof(double), 4, fpInt);
- }
-
- fclose(fpInt);
-
-
-
- // /////////////////////////////////////////////////////
- // Compute error
- // /////////////////////////////////////////////////////
-
- VectorType solBackup0 = intermediateSol3d;
- solBackup0 -= exactSol3d;
-
- RodDifferenceType rodDifference = computeGeodesicDifference(exactSolRod, intermediateSolRod);
-
- error = std::sqrt(EnergyNorm<MatrixType,VectorType>::normSquared(solBackup0, stiffnessMatrix3d)
- +
- EnergyNorm<BCRSMatrix<FieldMatrix<double,6,6> >,RodDifferenceType>::normSquared(rodDifference, hessianRod));
-
-
- double convRate = error / oldError;
- totalConvRate[i+1] = totalConvRate[i] * convRate;
-
- // Output
- std::cout << "DD iteration: " << i << " error : " << error << ", "
- << "convrate " << convRate
- << " total conv rate " << std::pow(totalConvRate[i+1], 1/((double)i+1)) << std::endl;
-
- // Convergence rates tend to stay fairly constant after a few initial iterates.
- // Only when we hit numerical dirt are they starting to wiggle around wildly.
- // We use this to detect 'the' convergence rate as the last averaged rate before
- // we hit the dirt.
- if (convRate > oldConvRate + 0.1 && i > 5 && firstTime) {
-
- std::cout << "Damping: " << damping
- << " convRate: " << std::pow(totalConvRate[i], 1/((double)i))
- << std::endl;
-
- std::ofstream convRateFile(filename.c_str());
- convRateFile << damping << " "
- << std::pow(totalConvRate[i], 1/((double)i))
- << std::endl;
-
- firstTime = false;
- }
-
- if (error < 1e-12)
- break;
-
- oldError = error;
- oldConvRate = convRate;
-
- }
-
- // Convergence without dirt: write the overall convergence rate now
- if (firstTime) {
- int backTrace = std::min(size_t(4),totalConvRate.size());
- std::cout << "Damping: " << damping
- << " convRate: " << std::pow(totalConvRate[i+1-backTrace], 1/((double)i+1-backTrace))
- << std::endl;
-
- std::ofstream convRateFile(filename.c_str());
- convRateFile << damping << " "
- << std::pow(totalConvRate[i+1-backTrace], 1/((double)i+1-backTrace))
- << std::endl;
- }
-
-
- // //////////////////////////////
- // Delete temporary memory
- // //////////////////////////////
- std::string removeTmpCommand = "rm -rf " + resultPath + "tmp/intermediate*";
- system(removeTmpCommand.c_str());
-
- // //////////////////////////////
- // Output result
- // //////////////////////////////
- LeafAmiraMeshWriter<GridType> amiraMeshWriter(grid);
- amiraMeshWriter.addVertexData(x3d, grid.leafView());
-
- BlockVector<FieldVector<double,1> > stress;
- Stress<GridType>::getStress(grid, x3d, stress, E, nu);
- amiraMeshWriter.addCellData(stress, grid.leafView());
-
- amiraMeshWriter.write(resultPath + "grid.result");
-
-
-
- writeRod(rodX, resultPath + "rod3d.result");
-
- } catch (Exception e) {
-
- std::cout << e << std::endl;
-
- }