Newer
Older
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif
#include <dune/common/exceptions.hh>
#include <dune/tectonic/utils/reductionfactors.hh>
#include "functionalfactory.hh"
#include "nonlinearsolver.hh"
void FixedPointIterationCounter::operator+=(
FixedPointIterationCounter const &other) {
iterations += other.iterations;
multigridIterations += other.multigridIterations;
}
template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::FixedPointIterator(
GlobalFriction& globalFriction,
const std::vector<const BitVector*>& bodywiseNonmortarBoundaries,
fixedPointMaxIterations_(parset.get<size_t>("v.fpi.maximumIterations")),
fixedPointTolerance_(parset.get<double>("v.fpi.tolerance")),
lambda_(parset.get<double>("v.fpi.lambda")),
solverParset_(parset.sub("v.solver")),
template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
template <class LinearSolver>
FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::run(
Updaters updaters, std::shared_ptr<LinearSolver>& linearSolver,
const std::vector<Matrix>& velocityMatrices, const std::vector<Vector>& velocityRHSs,
//std::cout << "FixedPointIterator::run()" << std::endl;
const auto nBodies = nBodyAssembler_.nGrids();
/*const auto& contactCouplings = nBodyAssembler_.getContactCouplings();
for (size_t i=0; i<contactCouplings.size(); i++) {
print(*contactCouplings[i]->nonmortarBoundary().getVertices(), "nonmortarBoundaries:");
FunctionalFactory<std::decay_t<decltype(nBodyAssembler_)>, decltype(globalFriction_), Matrix, Vector> functionalFactory(nBodyAssembler_, globalFriction_, velocityMatrices, velocityRHSs);
functionalFactory.build();
auto functional = functionalFactory.functional();
NonlinearSolver<std::remove_reference_t<decltype(*functional)>, IgnoreVector> solver(parset_.sub("solver.tnnmg"), linearSolver, functional, ignoreNodes_);
nBodyAssembler_.nodalToTransformed(velocityIterates, totalVelocityIterate);
for (fixedPointIteration = 0; fixedPointIteration < fixedPointMaxIterations_;
auto res = solver.solve(solverParset_, totalVelocityIterate);
nBodyAssembler_.postprocess(totalVelocityIterate, velocityIterates);
//nBodyAssembler_.postprocess(totalVelocityIterate, velocityIterates);
//print(totalVelocityIterate, "totalVelocityIterate:");
//print(velocityIterates, "velocityIterates:");
multigridIterations += res.iterations;
Dune::MatrixVector::addProduct(v_m, lambda_, totalVelocityIterate);
// extract relative velocities in mortar basis
std::vector<Vector> v_rel;
//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 << std::endl;
//std::cout << "-FixedPointIteration finished! " << (lambda_ < 1e-12 ? "lambda" : "breakCriterion") << std::endl;
//TODO: recently added, might be wrong or superfluous
globalFriction_.updateAlpha(alpha);
//print(alpha, "alpha: ");
//std::cout << "-FixedPointIteration finished with " << fixedPointIteration << " iterations, lambda " << lambda_ << "! " << 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;
FixedPointIterationCounter const &fpic) {
return stream << "(" << fpic.iterations << "," << fpic.multigridIterations
<< ")";
template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
void FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::split(
for (size_t i=0; i<splitV_.size(); i++) {
splitV_[i] = v[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"