Newer
Older
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <dune/common/exceptions.hh>
#include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/solvers/loopsolver.hh>
#include <dune/contact/assemblers/nbodyassembler.hh>
#include <dune/contact/common/dualbasisadapter.hh>
#include <dune/localfunctions/lagrange/pqkfactory.hh>
#include <dune/functions/gridfunctions/gridfunction.hh>
#include <dune/geometry/quadraturerules.hh>
#include <dune/geometry/type.hh>
#include <dune/geometry/referenceelements.hh>
#include <dune/fufem/functions/basisgridfunction.hh>
#include "../data-structures/enums.hh"
#include "../data-structures/enumparser.hh"
#include "fixedpointiterator.hh"
void FixedPointIterationCounter::operator+=(
FixedPointIterationCounter const &other) {
iterations += other.iterations;
multigridIterations += other.multigridIterations;
}
template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorm>
FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorm>::FixedPointIterator(
Dune::ParameterTree const &parset,
NBodyAssembler& nBodyAssembler,
GlobalFriction& globalFriction,
const std::vector<const BitVector*>& bodywiseNonmortarBoundaries,
const std::vector<const ErrorNorm* >& errorNorms)
fixedPointMaxIterations_(parset.get<size_t>("v.fpi.maximumIterations")),
fixedPointTolerance_(parset.get<double>("v.fpi.tolerance")),
lambda_(parset.get<double>("v.fpi.lambda")),
velocityMaxIterations_(parset.get<size_t>("v.solver.maximumIterations")),
velocityTolerance_(parset.get<double>("v.solver.tolerance")),
verbosity_(parset.get<Solver::VerbosityMode>("v.solver.verbosity")),
template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorm>
FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorm>::run(
Updaters updaters, const std::vector<Matrix>& velocityMatrices, const std::vector<Vector>& velocityRHSs,
const auto nBodies = nBodyAssembler_.nGrids();
std::vector<const Matrix*> matrices_ptr(nBodies);
// assemble full global contact problem
Matrix bilinearForm;
nBodyAssembler_.assembleJacobian(matrices_ptr, bilinearForm);
Vector totalRhs;
nBodyAssembler_.assembleRightHandSide(velocityRHSs, totalRhs);
// get lower and upper obstacles
const auto& totalObstacles = nBodyAssembler_.totalObstacles_;
Vector lower(totalObstacles.size());
Vector upper(totalObstacles.size());
for (size_t j=0; j<totalObstacles.size(); ++j) {
const auto& totalObstaclesj = totalObstacles[j];
auto& lowerj = lower[j];
auto& upperj = upper[j];
for (size_t d=0; d<dims; ++d) {
lowerj[d] = totalObstaclesj[d][0];
upperj[d] = totalObstaclesj[d][1];
}
}
print(totalObstacles, "totalObstacles:");
print(lower, "lower obstacles:");
print(upper, "upper obstacles:");
// comput velocity obstacles
Vector vLower, vUpper;
std::vector<Vector> u0, v0;
updaters.rate_->extractOldVelocity(v0);
updaters.rate_->extractOldDisplacement(u0);
Vector totalu0, totalv0;
nBodyAssembler_.concatenateVectors(u0, totalu0);
nBodyAssembler_.concatenateVectors(v0, totalv0);
updaters.rate_->velocityObstacles(totalu0, lower, totalv0, vLower);
updaters.rate_->velocityObstacles(totalu0, upper, totalv0, vUpper);
print(vLower, "vLower obstacles:");
print(vUpper, "vUpper obstacles:");
std::cout << "- Functional and TNNMG step setup: success" << std::endl;
LoopSolver<Vector> velocityProblemSolver(*step.get(), velocityMaxIterations_,
velocityTolerance_, energyNorm,
for (fixedPointIteration = 0; fixedPointIteration < fixedPointMaxIterations_;
nBodyAssembler_.nodalToTransformed(velocityIterates, totalVelocityIterate);
//print(velocityIterates, "velocityIterates:");
//print(totalVelocityIterate, "totalVelocityIterate:");
velocityProblemSolver.preprocess();
const auto& tnnmgSol = step->getSol();
nBodyAssembler_.postprocess(tnnmgSol, velocityIterates);
//nBodyAssembler_.postprocess(totalVelocityIterate, velocityIterates);
//print(totalVelocityIterate, "totalVelocityIterate:");
//print(velocityIterates, "velocityIterates:");
multigridIterations += velocityProblemSolver.getResult().iterations;
// extract relative velocities in mortar basis
std::vector<Vector> v_rel;
//print(v_rel, "v_rel");
std::cout << "- State problem set" << std::endl;
updaters.state_->extractAlpha(newAlpha);
bool breakCriterion = true;
if (alpha[i].size()==0 || newAlpha[i].size()==0)
continue;
//print(alpha[i], "alpha i:");
//print(newAlpha[i], "new alpha i:");
std::cout << "fixedPoint error: " << errorNorms_[i]->diff(alpha[i], newAlpha[i]) << std::endl;
std::cout << "-FixedPointIteration finished! " << (lambda_ < 1e-12 ? "lambda" : "breakCriterion") << std::endl;
if (fixedPointIteration == fixedPointMaxIterations_)
DUNE_THROW(Dune::Exception, "FPI failed to converge");
// Cannot use return { fixedPointIteration, multigridIterations };
// with gcc 4.9.2, see also http://stackoverflow.com/a/37777814/179927
FixedPointIterationCounter ret;
ret.iterations = fixedPointIteration;
ret.multigridIterations = multigridIterations;
return ret;
}
std::ostream &operator<<(std::ostream &stream,
FixedPointIterationCounter const &fpic) {
return stream << "(" << fpic.iterations << "," << fpic.multigridIterations
<< ")";
template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorm>
void FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorm>::relativeVelocities(
const Vector& v,
std::vector<Vector>& v_rel) const {
const size_t nBodies = nBodyAssembler_.nGrids();
// const auto& contactCouplings = nBodyAssembler_.getContactCouplings();
for (size_t bodyIdx=0; bodyIdx<nBodies; bodyIdx++) {
const auto& nonmortarBoundary = *bodywiseNonmortarBoundaries_[bodyIdx];
auto& v_rel_ = v_rel[bodyIdx];
v_rel_.resize(nonmortarBoundary.size());
for (size_t i=0; i<v_rel_.size(); i++) {
if (toBool(nonmortarBoundary[i])) {
v_rel_[i] = v[globalIdx];
}
globalIdx++;
}
}
const auto gridView = contactCouplings[couplingIndices[0]]->nonmortarBoundary().gridView();
Dune::MultipleCodimMultipleGeomTypeMapper<
decltype(gridView), Dune::MCMGVertexLayout> const vertexMapper(gridView, Dune::mcmgVertexLayout());
for (auto it = gridView.template begin<block_size>(); it != gridView.template end<block_size>(); ++it) {
const auto i = vertexMapper.index(*it);
for (size_t j=0; j<couplingIndices.size(); j++) {
const auto couplingIdx = couplingIndices[j];
if (not contactCouplings[couplingIdx]->nonmortarBoundary().containsVertex(i))
continue;
localToGlobal_.emplace_back(i);
restrictions_.emplace_back(weights[bodyIdx][i], weightedNormalStress[bodyIdx][i],
couplings[i]->frictionData()(geoToPoint(it->geometry())));
break;
#include "fixedpointiterator_tmpl.cc"