#include <config.h>
#include <dune/grid/onedgrid.hh>
#include <dune/grid/uggrid.hh>
#include <dune/istl/io.hh>
#include <dune/istl/solvers.hh>
#include <dune/grid/io/file/amirameshreader.hh>
#include <dune/grid/io/file/amirameshwriter.hh>
#include <dune/common/bitsetvector.hh>
#include <dune/common/configparser.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/fufem/readbitfield.hh>
#include <dune/solvers/norms/energynorm.hh>
#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/prolongboundarypatch.hh>
#include <dune/fufem/sampleonbitfield.hh>
#include <dune/fufem/neumannassembler.hh>
#include <dune/fufem/computestress.hh>
#include <dune/fufem/functionspacebases/q1nodalbasis.hh>
#include <dune/fufem/assemblers/operatorassembler.hh>
#include <dune/fufem/assemblers/localassemblers/stvenantkirchhoffassembler.hh>
#include <dune/gfe/quaternion.hh>
#include <dune/gfe/rodassembler.hh>
#include <dune/gfe/rigidbodymotion.hh>
#include <dune/gfe/averageinterface.hh>
#include <dune/gfe/riemanniantrsolver.hh>
#include <dune/gfe/geodesicdifference.hh>
#include <dune/gfe/rodwriter.hh>
#include <dune/gfe/makestraightrod.hh>
#include <dune/gfe/coupling/rodcontinuumcomplex.hh>
// Space dimension
const int dim = 3;
using namespace Dune;
using std::string;
using std::vector;
// Some types that I need
//typedef BCRSMatrix<FieldMatrix<double, dim, dim> > OperatorType;
//typedef BlockVector<FieldVector<double, dim> > VectorType;
typedef vector<RigidBodyMotion<dim> > RodSolutionType;
typedef BlockVector<FieldVector<double, 6> > RodDifferenceType;
template <class GridView, class MatrixType, class VectorType>
class LinearizedContinuumNeumannToDirichletMap
{
public:
/** \brief Constructor
*
*/
LinearizedContinuumNeumannToDirichletMap(const BoundaryPatchBase<GridView>& interface,
const VectorType& weakVolumeAndNeumannTerm,
const VectorType& dirichletValues,
const LinearLocalAssembler<typename GridView::Grid,
typename P1NodalBasis<GridView,double>::LocalFiniteElement,
typename P1NodalBasis<GridView,double>::LocalFiniteElement,
Dune::FieldMatrix<double,3,3> >* localAssembler,
const shared_ptr< ::LoopSolver<VectorType> >& solver
)
: interface_(interface),
weakVolumeAndNeumannTerm_(weakVolumeAndNeumannTerm),
dirichletValues_(dirichletValues),
solver_(solver),
localAssembler_(localAssembler)
{}
/** \brief Map a Neumann value to a Dirichlet value
*
* \param currentIterate The continuum configuration that we are linearizing about
*
* \return The infinitesimal movement of the interface
*/
FieldVector<double,6> apply(const VectorType& currentIterate,
const Dune::FieldVector<double,3>& force,
const Dune::FieldVector<double,3>& torque,
const Dune::FieldVector<double,3>& centerOfTorque
)
{
////////////////////////////////////////////////////
// Assemble the linearized problem
////////////////////////////////////////////////////
/** \todo We are actually assembling standard linear elasticity,
* i.e. the linearization at zero
*/
typedef P1NodalBasis<GridView,double> P1Basis;
P1Basis basis(interface_.gridView());
OperatorAssembler<P1Basis,P1Basis> assembler(basis, basis);
MatrixType stiffnessMatrix;
assembler.assemble(*localAssembler_, stiffnessMatrix);
/////////////////////////////////////////////////////////////////////
// Turn the input force and torque into a Neumann boundary field
/////////////////////////////////////////////////////////////////////
VectorType neumannValues(stiffnessMatrix.N());
neumannValues = 0;
//
computeAveragePressure<GridView>(force, torque,
interface_,
centerOfTorque,
neumannValues);
// The weak form of the Neumann data
VectorType rhs = weakVolumeAndNeumannTerm_;
assembleAndAddNeumannTerm<GridView, VectorType>(interface_,
neumannValues,
rhs);
// Solve the Neumann problem for the continuum
VectorType x = dirichletValues_;
assert( (dynamic_cast<LinearIterationStep<MatrixType,VectorType>* >(solver_->iterationStep_)) );
dynamic_cast<LinearIterationStep<MatrixType,VectorType>* >(solver_->iterationStep_)->setProblem(stiffnessMatrix, x, rhs);
//solver.preprocess();
solver_->iterationStep_->preprocess();
solver_->solve();
x = solver_->iterationStep_->getSol();
std::cout << "x:\n" << x << std::endl;
// Average the continuum displacement on the coupling boundary
RigidBodyMotion<3> averageInterface;
computeAverageInterface(interface_, x, averageInterface);
std::cout << "Average interface: " << averageInterface << std::endl;
// Compute the linearization
FieldVector<double,6> interfaceCorrection;
interfaceCorrection[0] = averageInterface.r[0];
interfaceCorrection[1] = averageInterface.r[1];
interfaceCorrection[2] = averageInterface.r[2];
FieldVector<double,3> infinitesimalRotation = Rotation<3,double>::expInv(averageInterface.q);
interfaceCorrection[3] = infinitesimalRotation[0];
interfaceCorrection[4] = infinitesimalRotation[1];
interfaceCorrection[5] = infinitesimalRotation[2];
return interfaceCorrection;
}
private:
const VectorType& weakVolumeAndNeumannTerm_;
const VectorType& dirichletValues_;
const shared_ptr< ::LoopSolver<VectorType> > solver_;
const BoundaryPatchBase<GridView>& interface_;
const LinearLocalAssembler<typename GridView::Grid,
typename P1NodalBasis<GridView,double>::LocalFiniteElement,
typename P1NodalBasis<GridView,double>::LocalFiniteElement,
Dune::FieldMatrix<double,3,3> >* localAssembler_;
};
template <class GridView, class VectorType>
class LinearizedRodNeumannToDirichletMap
{
public:
/** \brief Constructor
*
*/
LinearizedRodNeumannToDirichletMap(const BoundaryPatchBase<GridView>& interface,
LocalGeodesicFEStiffness<GridView, RigidBodyMotion<3> >* localAssembler)
: interface_(interface),
localAssembler_(localAssembler)
{}
/** \brief Map a Neumann value to a Dirichlet value
*
* \param currentIterate The rod configuration that we are linearizing about
*
* \return The Dirichlet value
*/
Dune::FieldVector<double,6> apply(const RodSolutionType& currentIterate,
const Dune::FieldVector<double,3>& force,
const Dune::FieldVector<double,3>& torque,
const Dune::FieldVector<double,3>& centerOfTorque
)
{
////////////////////////////////////////////////////
// Assemble the linearized rod problem
////////////////////////////////////////////////////
typedef BCRSMatrix<FieldMatrix<double,6,6> > MatrixType;
GeodesicFEAssembler<GridView, RigidBodyMotion<3> > assembler(interface_.gridView(),
localAssembler_);
MatrixType stiffnessMatrix;
assembler.assembleMatrix(currentIterate,
stiffnessMatrix,
true // assemble occupation pattern
);
VectorType rhs(currentIterate.size());
rhs = 0;
assembler.assembleGradient(currentIterate, rhs);
// The right hand side is the _negative_ gradient
rhs *= -1;
/////////////////////////////////////////////////////////////////////
// Turn the input force and torque into a Neumann boundary field
/////////////////////////////////////////////////////////////////////
// The weak form of the Neumann data
rhs[0][0] += force[0];
rhs[0][1] += force[1];
rhs[0][2] += force[2];
rhs[0][3] += torque[0];
rhs[0][4] += torque[1];
rhs[0][5] += torque[2];
///////////////////////////////////////////////////////////
// Modify matrix and rhs to contain one Dirichlet node
///////////////////////////////////////////////////////////
int dIdx = rhs.size()-1; // hardwired index of the single Dirichlet node
rhs[dIdx] = 0;
stiffnessMatrix[dIdx] = 0;
stiffnessMatrix[dIdx][dIdx] = Dune::ScaledIdentityMatrix<double,6>(1);
//////////////////////////////////////////////////
// Solve the Neumann problem
//////////////////////////////////////////////////
VectorType x(rhs.size());
x = 0; // initial iterate
// Technicality: turn the matrix into a linear operator
Dune::MatrixAdapter<MatrixType,VectorType,VectorType> op(stiffnessMatrix);
// A preconditioner
Dune::SeqILU0<MatrixType,VectorType,VectorType> ilu0(stiffnessMatrix,1.0);
// A preconditioned conjugate-gradient solver
Dune::CGSolver<VectorType> cg(op,ilu0,1E-4,100,2);
// Object storing some statistics about the solving process
Dune::InverseOperatorResult statistics;
// Solve!
cg.apply(x, rhs, statistics);
std::cout << "x:\n" << x << std::endl;
std::cout << "Linear rod interface correction: " << x[0] << std::endl;
return x[0];
}
private:
const BoundaryPatchBase<GridView>& interface_;
LocalGeodesicFEStiffness<GridView, RigidBodyMotion<3> >* localAssembler_;
};
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("dirneucoupling.parset");
// read solver settings
const int numLevels = parameterSet.get<int>("numLevels");
string ddType = parameterSet.get<string>("ddType");
string preconditioner = parameterSet.get<string>("preconditioner");
const double ddTolerance = parameterSet.get<double>("ddTolerance");
const int maxDDIterations = parameterSet.get<int>("maxDirichletNeumannSteps");
const double damping = parameterSet.get<double>("damping");
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");
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] = parameterSet.get<FieldVector<double,3> >("rodRestEndPoint0");
rodRestEndPoint[1] = parameterSet.get<FieldVector<double,3> >("rodRestEndPoint1");
//////////////////////////////////////////////////////////////////
// Print the algorithm type so we have it in the log files
//////////////////////////////////////////////////////////////////
std::cout << "Algorithm: " << ddType << std::endl;
if (ddType=="RichardsonIteration")
std::cout << "Preconditioner: " << preconditioner << std::endl;
// ///////////////////////////////////////////////
// Create the rod grid and continuum grid
// ///////////////////////////////////////////////
typedef OneDGrid RodGridType;
typedef UGGrid<dim> GridType;
RodContinuumComplex<RodGridType,GridType> complex;
complex.rodGrids_["rod"] = shared_ptr<RodGridType>
(new RodGridType(numRodBaseElements, 0, (rodRestEndPoint[1]-rodRestEndPoint[0]).two_norm()));
complex.continuumGrids_["continuum"] = shared_ptr<GridType>(AmiraMeshReader<GridType>::read(path + objectName));
complex.continuumGrids_["continuum"]->setRefinementType(GridType::COPY);
// //////////////////////////////////////
// Globally refine grids
// //////////////////////////////////////
complex.rodGrids_["rod"]->globalRefine(numLevels-1);
complex.continuumGrids_["continuum"]->globalRefine(numLevels-1);
std::vector<BitSetVector<dim> > dirichletNodes(1);
RodSolutionType rodX(complex.rodGrids_["rod"]->size(1));
// //////////////////////////
// Initial solution
// //////////////////////////
makeStraightRod(rodX, complex.rodGrids_["rod"]->size(1), rodRestEndPoint[0], rodRestEndPoint[1]);
// /////////////////////////////////////////
// Read Dirichlet values
// /////////////////////////////////////////
rodX.back().r = parameterSet.get("dirichletValue", rodRestEndPoint[1]);
FieldVector<double,3> axis = parameterSet.get("dirichletAxis", FieldVector<double,3>(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 = complex.rodGrids_["rod"]->maxLevel();
// /////////////////////////////////////////////////////
// Determine the Dirichlet nodes
// /////////////////////////////////////////////////////
std::vector<VectorType> dirichletValues;
dirichletValues.resize(toplevel+1);
dirichletValues[0].resize(complex.continuumGrids_["continuum"]->size(0, dim));
AmiraMeshReader<int>::readFunction(dirichletValues[0], path + dirichletValuesFile);
std::vector<LevelBoundaryPatch<GridType> > dirichletBoundary;
dirichletBoundary.resize(numLevels);
dirichletBoundary[0].setup(*complex.continuumGrids_["continuum"], 0);
readBoundaryPatch(dirichletBoundary[0], path + dirichletNodesFile);
PatchProlongator<GridType>::prolong(dirichletBoundary);
dirichletNodes.resize(toplevel+1);
for (int i=0; i<=toplevel; i++) {
dirichletNodes[i].resize( complex.continuumGrids_["continuum"]->size(i,dim));
for (int j=0; j<complex.continuumGrids_["continuum"]->size(i,dim); j++)
dirichletNodes[i][j] = dirichletBoundary[i].containsVertex(j);
}
sampleOnBitField(*complex.continuumGrids_["continuum"], dirichletValues, dirichletNodes);
// /////////////////////////////////////////////////////
// Determine the interface boundary
// /////////////////////////////////////////////////////
std::vector<LevelBoundaryPatch<GridType> > interfaceBoundary;
interfaceBoundary.resize(numLevels);
interfaceBoundary[0].setup(*complex.continuumGrids_["continuum"], 0);
readBoundaryPatch(interfaceBoundary[0], path + interfaceNodesFile);
PatchProlongator<GridType>::prolong(interfaceBoundary);
// //////////////////////////////////////////
// Assemble 3d linear elasticity problem
// //////////////////////////////////////////
typedef P1NodalBasis<GridType::LeafGridView,double> FEBasis;
FEBasis basis(complex.continuumGrids_["continuum"]->leafView());
OperatorAssembler<FEBasis,FEBasis> assembler(basis, basis);
StVenantKirchhoffAssembler<GridType, FEBasis::LocalFiniteElement, FEBasis::LocalFiniteElement> localAssembler(E, nu);
MatrixType stiffnessMatrix3d;
assembler.assemble(localAssembler, stiffnessMatrix3d);
// ////////////////////////////////////////////////////////////
// Create solution and rhs vectors
// ////////////////////////////////////////////////////////////
VectorType x3d(complex.continuumGrids_["continuum"]->size(toplevel,dim));
VectorType rhs3d(complex.continuumGrids_["continuum"]->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[toplevel][i][j])
x3d[i][j] = dirichletValues[toplevel][i][j];
// ///////////////////////////////////////////
// Dirichlet nodes for the rod problem
// ///////////////////////////////////////////
BitSetVector<6> rodDirichletNodes(complex.rodGrids_["rod"]->size(1));
rodDirichletNodes.unsetAll();
rodDirichletNodes[0] = true;
rodDirichletNodes.back() = true;
// ///////////////////////////////////////////
// Create a solver for the rod problem
// ///////////////////////////////////////////
RodLocalStiffness<RodGridType::LeafGridView,double> rodLocalStiffness(complex.rodGrids_["rod"]->leafView(),
rodA, rodJ1, rodJ2, rodE, rodNu);
RodAssembler<RodGridType::LeafGridView,3> rodAssembler(complex.rodGrids_["rod"]->leafView(), &rodLocalStiffness);
RiemannianTrustRegionSolver<RodGridType,RigidBodyMotion<3> > rodSolver;
rodSolver.setup(*complex.rodGrids_["rod"],
&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.back();
multigridStep.basesolver_ = &baseSolver;
multigridStep.setSmoother(&presmoother, &postsmoother);
multigridStep.verbosity_ = Solver::QUIET;
EnergyNorm<MatrixType, VectorType> energyNorm(multigridStep);
shared_ptr< ::LoopSolver<VectorType> > solver = shared_ptr< ::LoopSolver<VectorType> >( new ::LoopSolver<VectorType>(&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(*complex.continuumGrids_["continuum"],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);
//
double normOfOldCorrection = 0;
int dnStepsActuallyTaken = 0;
for (int i=0; i<maxDDIterations; i++) {
std::cout << "----------------------------------------------------" << std::endl;
std::cout << " Domain-Decomposition- 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;
RigidBodyMotion<3> averageInterface;
if (ddType=="FixedPointIteration") {
// //////////////////////////////////////////////////
// Dirichlet step for the rod
// //////////////////////////////////////////////////
rodX[0] = lambda;
rodSolver.setInitialSolution(rodX);
rodSolver.solve();
rodX = rodSolver.getSol();
// for (int j=0; j<rodX.size(); j++)
// std::cout << rodX[j] << std::endl;
// ///////////////////////////////////////////////////////////
// Extract Neumann values and transfer it to the 3d object
// ///////////////////////////////////////////////////////////
BitSetVector<1> couplingBitfield(rodX.size(),false);
// Using that index 0 is always the left boundary for a uniformly refined OneDGrid
couplingBitfield[0] = true;
LeafBoundaryPatch<RodGridType> couplingBoundary(*complex.rodGrids_["rod"], couplingBitfield);
FieldVector<double,dim> resultantForce, resultantTorque;
resultantForce = rodAssembler.getResultantForce(couplingBoundary, rodX, resultantTorque);
// Flip orientation
resultantForce *= -1;
resultantTorque *= -1;
std::cout << "resultant force: " << resultantForce << std::endl;
std::cout << "resultant torque: " << resultantTorque << std::endl;
VectorType neumannValues(rhs3d.size());
// Using that index 0 is always the left boundary for a uniformly refined OneDGrid
computeAveragePressure<GridType::LevelGridView>(resultantForce, resultantTorque,
interfaceBoundary[complex.continuumGrids_["continuum"]->maxLevel()],
rodX[0].r,
neumannValues);
rhs3d = 0;
assembleAndAddNeumannTerm<GridType::LevelGridView, VectorType>(interfaceBoundary[complex.continuumGrids_["continuum"]->maxLevel()],
neumannValues,
rhs3d);
// ///////////////////////////////////////////////////////////
// Solve the Neumann problem for the continuum
// ///////////////////////////////////////////////////////////
multigridStep.setProblem(stiffnessMatrix3d, x3d, rhs3d, complex.continuumGrids_["continuum"]->maxLevel()+1);
solver->preprocess();
multigridStep.preprocess();
solver->solve();
x3d = multigridStep.getSol();
// ///////////////////////////////////////////////////////////
// Extract new interface position and orientation
// ///////////////////////////////////////////////////////////
computeAverageInterface(interfaceBoundary[toplevel], x3d, averageInterface);
//averageInterface.r[0] = averageInterface.r[1] = 0;
//averageInterface.q = Quaternion<double>::identity();
std::cout << "average interface: " << averageInterface << std::endl;
std::cout << "director 0: " << averageInterface.q.director(0) << std::endl;
std::cout << "director 1: " << averageInterface.q.director(1) << std::endl;
std::cout << "director 2: " << averageInterface.q.director(2) << std::endl;
std::cout << std::endl;
//////////////////////////////////////////////////////////////
// Compute new damped interface value
//////////////////////////////////////////////////////////////
for (int j=0; j<dim; j++)
lambda.r[j] = (1-damping) * lambda.r[j]
+ damping * (referenceInterface.r[j] + averageInterface.r[j]);
lambda.q = Rotation<3,double>::interpolate(lambda.q,
referenceInterface.q.mult(averageInterface.q),
damping);
} else if (ddType=="RichardsonIteration") {
///////////////////////////////////////////////////////////////////
// One preconditioned Richardson step
///////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////
// Evaluate the Dirichlet-to-Neumann map for the rod
///////////////////////////////////////////////////////////////////
// solve a Dirichlet problem for the rod
rodX[0] = lambda;
rodSolver.setInitialSolution(rodX);
rodSolver.solve();
rodX = rodSolver.getSol();
// Extract Neumann values
BitSetVector<1> couplingBitfield(rodX.size(),false);
// Using that index 0 is always the left boundary for a uniformly refined OneDGrid
couplingBitfield[0] = true;
LeafBoundaryPatch<RodGridType> couplingBoundary(*complex.rodGrids_["rod"], couplingBitfield);
FieldVector<double,dim> resultantForce, resultantTorque;
resultantForce = rodAssembler.getResultantForce(couplingBoundary, rodX, resultantTorque);
// Flip orientation
resultantForce *= -1;
resultantTorque *= -1;
std::cout << "resultant force: " << resultantForce << std::endl;
std::cout << "resultant torque: " << resultantTorque << std::endl;
///////////////////////////////////////////////////////////////////
// Evaluate the Dirichlet-to-Neumann map for the continuum
///////////////////////////////////////////////////////////////////
// Turn \lambda \in TSE(3) into a Dirichlet value for the continuum
RigidBodyMotion<3> relativeMovement;
relativeMovement.r = lambda.r - referenceInterface.r;
relativeMovement.q = referenceInterface.q;
relativeMovement.q.invert();
relativeMovement.q = lambda.q.mult(relativeMovement.q);
setRotation(interfaceBoundary.back(), x3d, relativeMovement);
// Solve the Dirichlet problem
multigridStep.setProblem(stiffnessMatrix3d, x3d, rhs3d, complex.continuumGrids_["continuum"]->maxLevel()+1);
solver->preprocess();
multigridStep.preprocess();
solver->solve();
x3d = multigridStep.getSol();
// Integrate over the residual on the coupling boundary to obtain
// an element of T^*SE.
FieldVector<double,3> continuumForce, continuumTorque;
VectorType residual = rhs3d;
stiffnessMatrix3d.mmv(x3d,residual);
/** \todo Is referenceInterface.r the correct center of rotation? */
computeTotalForceAndTorque(interfaceBoundary.back(), residual, referenceInterface.r,
continuumForce, continuumTorque);
///////////////////////////////////////////////////////////////
// Compute the overall Steklov-Poincare residual
///////////////////////////////////////////////////////////////
FieldVector<double,3> residualForce = resultantForce + continuumForce;
FieldVector<double,3> residualTorque = resultantTorque + continuumTorque;
///////////////////////////////////////////////////////////////
// Apply the preconditioner
///////////////////////////////////////////////////////////////
FieldVector<double,6> interfaceCorrection;
if (preconditioner=="DirichletNeumann") {
////////////////////////////////////////////////////////////////////
// Preconditioner is the linearized Neumann-to-Dirichlet map
// of the continuum.
////////////////////////////////////////////////////////////////////
LinearizedContinuumNeumannToDirichletMap<GridType::LevelGridView,MatrixType,VectorType>
linContNtDMap(interfaceBoundary.back(),
rhs3d,
dirichletValues.back(),
&localAssembler,
solver);
interfaceCorrection = linContNtDMap.apply(x3d, residualForce, residualTorque, rodX[0].r);
} else if (preconditioner=="NeumannDirichlet") {
////////////////////////////////////////////////////////////////////
// Preconditioner is the linearized Neumann-to-Dirichlet map
// of the rod.
////////////////////////////////////////////////////////////////////
typedef BlockVector<FieldVector<double,6> > RodVectorType;
LinearizedRodNeumannToDirichletMap<RodGridType::LeafGridView,RodVectorType> linRodNtDMap(couplingBoundary,
&rodLocalStiffness);
interfaceCorrection = linRodNtDMap.apply(rodX, residualForce, residualTorque, rodX[0].r);
} else if (preconditioner=="NeumannNeumann") {
////////////////////////////////////////////////////////////////////
// Preconditioner is a convex combination of the linearized
// Neumann-to-Dirichlet map of the continuum and the linearized
// Neumann-to-Dirichlet map of the rod.
////////////////////////////////////////////////////////////////////
LinearizedContinuumNeumannToDirichletMap<GridType::LevelGridView,MatrixType,VectorType>
linContNtDMap(interfaceBoundary.back(),
rhs3d,
dirichletValues.back(),
&localAssembler,
solver);
typedef BlockVector<FieldVector<double,6> > RodVectorType;
LinearizedRodNeumannToDirichletMap<RodGridType::LeafGridView,RodVectorType> linRodNtDMap(couplingBoundary,
&rodLocalStiffness);
array<double, 2> alpha = parameterSet.get<array<double,2> >("neumannNeumannWeights");
FieldVector<double,6> continuumCorrection = linContNtDMap.apply(x3d, residualForce, residualTorque, rodX[0].r);
FieldVector<double,6> rodCorrection = linRodNtDMap.apply(rodX, residualForce, residualTorque, rodX[0].r);
for (int j=0; j<6; j++)
interfaceCorrection[j] = (alpha[0] * continuumCorrection[j] + alpha[1] * rodCorrection[j])
/ alpha[0] + alpha[1];
} else if (preconditioner=="RobinRobin") {
DUNE_THROW(NotImplemented, "Robin-Robin preconditioner not implemented yet");
} else
DUNE_THROW(NotImplemented, preconditioner << " is not a known preconditioner");
///////////////////////////////////////////////////////////////////////////////
// Apply the damped correction to the current interface value
///////////////////////////////////////////////////////////////////////////////
interfaceCorrection *= damping;
lambda = RigidBodyMotion<3>::exp(lambda, interfaceCorrection);
} else
DUNE_THROW(NotImplemented, ddType << " is not a known domain decomposition algorithm");
std::cout << "Lambda: " << lambda << 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, complex.continuumGrids_["continuum"]->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(maxDDIterations+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(*complex.continuumGrids_["continuum"]);
amiraMeshWriter.addVertexData(x3d, complex.continuumGrids_["continuum"]->leafView());
BlockVector<FieldVector<double,1> > stress;
Stress<GridType>::getStress(*complex.continuumGrids_["continuum"], x3d, stress, E, nu);
amiraMeshWriter.addCellData(stress, complex.continuumGrids_["continuum"]->leafView());
amiraMeshWriter.write(resultPath + "grid.result");
writeRod(rodX, resultPath + "rod3d.result");
} catch (Exception e) {
std::cout << e << std::endl;
}