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 <dune/solvers/solvers/loopsolver.hh>
#include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
#include <dune/solvers/iterationsteps/multigridstep.hh>
#include "../spatial-solving/contact/nbodycontacttransfer.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")),
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 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& contactCouplings = nBodyAssembler_.getContactCouplings();
for (size_t i=0; i<contactCouplings.size(); i++) {
print(*contactCouplings[i]->nonmortarBoundary().getVertices(), "nonmortarBoundaries:");
// assemble full global contact problem
Matrix bilinearForm;
nBodyAssembler_.assembleJacobian(matrices_ptr, bilinearForm);
nBodyAssembler_.assembleRightHandSide(velocityRHSs, totalRhs);
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(lower, "lower obstacles:");
//print(upper, "upper obstacles:");
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, upper, totalv0, vUpper); */
//print(vLower, "vLower obstacles:");
//print(vUpper, "vUpper obstacles:");
//std::cout << "- Problem assembled: success" << std::endl;
//using ZeroSolverFactory = SolverFactory<Functional, IgnoreVector>;
//Functional J(bilinearForm, totalRhs, ZeroNonlinearity(), vLower, vUpper);
//ZeroSolverFactory solverFactory(parset_.sub("solver.tnnmg"), J, mgSolver, ignoreNodes_);
Functional J(bilinearForm, totalRhs, globalFriction_, lower, upper);
Factory solverFactory(parset_.sub("solver.tnnmg"), J, ignoreNodes_);
solverFactory.build(linearSolver);
//std::cout << "- Functional and TNNMG step setup: success" << std::endl;
LoopSolver<Vector> velocityProblemSolver(*step.get(), velocityMaxIterations_,
velocityTolerance_, energyNorm,
Vector totalVelocityIterate;
nBodyAssembler_.nodalToTransformed(velocityIterates, totalVelocityIterate);
// project in onto admissible set
const size_t blocksize = Vector::block_type::dimension;
for (size_t i=0; i<totalVelocityIterate.size(); i++) {
for (size_t j=0; j<blocksize; j++) {
if (totalVelocityIterate[i][j] < lower[i][j]) {
totalVelocityIterate[i][j] = lower[i][j];
}
if (totalVelocityIterate[i][j] > upper[i][j]) {
totalVelocityIterate[i][j] = upper[i][j];
}
}
}
Vector old_v = totalVelocityIterate;
for (fixedPointIteration = 0; fixedPointIteration < fixedPointMaxIterations_;
//print(velocityIterates, "velocityIterates:");
//print(totalVelocityIterate, "totalVelocityIterate:");
//std::cout << "- FixedPointIteration iterate" << std::endl;
//std::cout << "- Velocity problem set" << std::endl;
velocityProblemSolver.preprocess();
//std::cout << "FixPointIterator: Energy of TNNMG solution: " << J(tnnmgSol) << std::endl;
nBodyAssembler_.postprocess(tnnmgSol, velocityIterates);
//nBodyAssembler_.postprocess(totalVelocityIterate, velocityIterates);
//print(totalVelocityIterate, "totalVelocityIterate:");
//print(velocityIterates, "velocityIterates:");
multigridIterations += velocityProblemSolver.getResult().iterations;
Vector v_m = old_v;
v_m *= 1.0 - lambda_;
Dune::MatrixVector::addProduct(v_m, lambda_, tnnmgSol);
// extract relative velocities in mortar basis
std::vector<Vector> v_rel;
relativeVelocities(v_m, v_rel);
//print(v_m, "v_m: ");
//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;
//TODO: recently added, might be wrong or superfluous
globalFriction_.updateAlpha(alpha);
//print(alpha, "alpha: ");
//std::cout << "-FixedPointIteration finished! " << 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>::relativeVelocities(
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"