Skip to content
Snippets Groups Projects
fixedpointiterator.cc 11.5 KiB
Newer Older
#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include <dune/common/exceptions.hh>

podlesny's avatar
.  
podlesny committed
#include <dune/matrix-vector/axpy.hh>

#include <dune/solvers/norms/energynorm.hh>
#include <dune/solvers/solvers/loopsolver.hh>

podlesny's avatar
.  
podlesny committed
#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>

podlesny's avatar
podlesny committed
#include "../data-structures/enums.hh"
#include "../data-structures/enumparser.hh"
podlesny's avatar
.  
podlesny committed
#include "../utils/tobool.hh"
podlesny's avatar
podlesny committed
#include "../utils/debugutils.hh"

podlesny's avatar
.  
podlesny committed
#include <dune/solvers/solvers/loopsolver.hh>
#include <dune/solvers/iterationsteps/truncatedblockgsstep.hh>
#include <dune/solvers/iterationsteps/multigridstep.hh>

podlesny's avatar
podlesny committed
#include "../spatial-solving/contact/nbodycontacttransfer.hh"
podlesny's avatar
.  
podlesny committed

podlesny's avatar
.  
podlesny committed
#include "tnnmg/functional.hh"
#include "tnnmg/zerononlinearity.hh"
podlesny's avatar
podlesny committed

#include "fixedpointiterator.hh"
podlesny's avatar
.  
podlesny committed

void FixedPointIterationCounter::operator+=(
    FixedPointIterationCounter const &other) {
  iterations += other.iterations;
  multigridIterations += other.multigridIterations;
}

podlesny's avatar
podlesny committed
template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::FixedPointIterator(
podlesny's avatar
podlesny committed
    Dune::ParameterTree const &parset,
podlesny's avatar
podlesny committed
    const NBodyAssembler& nBodyAssembler,
podlesny's avatar
.  
podlesny committed
    const IgnoreVector& ignoreNodes,
podlesny's avatar
.  
podlesny committed
    GlobalFriction& globalFriction,
    const std::vector<const BitVector*>& bodywiseNonmortarBoundaries,
podlesny's avatar
.  
podlesny committed
    const ErrorNorms& errorNorms)
podlesny's avatar
podlesny committed
    : parset_(parset),
podlesny's avatar
podlesny committed
      nBodyAssembler_(nBodyAssembler),
podlesny's avatar
.  
podlesny committed
      ignoreNodes_(ignoreNodes),
      globalFriction_(globalFriction),
podlesny's avatar
.  
podlesny committed
      bodywiseNonmortarBoundaries_(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")),
podlesny's avatar
.  
podlesny committed
      errorNorms_(errorNorms) {}
podlesny's avatar
podlesny committed
template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
template <class LinearSolver>
Elias Pipping's avatar
Elias Pipping committed
FixedPointIterationCounter
podlesny's avatar
podlesny committed
FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::run(
    Updaters updaters, std::shared_ptr<LinearSolver>& linearSolver,
podlesny's avatar
.  
podlesny committed
    const std::vector<Matrix>& velocityMatrices, const std::vector<Vector>& velocityRHSs,
podlesny's avatar
.  
podlesny committed
    std::vector<Vector>& velocityIterates) {

podlesny's avatar
podlesny committed
  //std::cout << "FixedPointIterator::run()" << std::endl;
podlesny's avatar
.  
podlesny committed

podlesny's avatar
.  
podlesny committed
  // debugging
podlesny's avatar
podlesny committed
  /*const auto& contactCouplings = nBodyAssembler_.getContactCouplings();
podlesny's avatar
.  
podlesny committed
  for (size_t i=0; i<contactCouplings.size(); i++) {
    print(*contactCouplings[i]->nonmortarBoundary().getVertices(), "nonmortarBoundaries:");
podlesny's avatar
podlesny committed
  }*/
podlesny's avatar
.  
podlesny committed

podlesny's avatar
podlesny committed
  const auto nBodies = nBodyAssembler_.nGrids();
podlesny's avatar
.  
podlesny committed

  std::vector<const Matrix*> matrices_ptr(nBodies);
podlesny's avatar
.  
podlesny committed
  for (int i=0; i<nBodies; i++) {
podlesny's avatar
.  
podlesny committed
      matrices_ptr[i] = &velocityMatrices[i];
  }
podlesny's avatar
.  
podlesny committed

podlesny's avatar
podlesny committed
  // assemble full global contact problem
  Matrix bilinearForm;
podlesny's avatar
podlesny committed
  nBodyAssembler_.assembleJacobian(matrices_ptr, bilinearForm);
podlesny's avatar
podlesny committed

podlesny's avatar
podlesny committed
  //print(bilinearForm, "bilinearForm:");
podlesny's avatar
.  
podlesny committed

podlesny's avatar
podlesny committed
  Vector totalRhs;
podlesny's avatar
podlesny committed
  nBodyAssembler_.assembleRightHandSide(velocityRHSs, totalRhs);
podlesny's avatar
podlesny committed

podlesny's avatar
podlesny committed
  //print(totalRhs, "totalRhs:");
podlesny's avatar
.  
podlesny committed

podlesny's avatar
podlesny committed
  // get lower and upper obstacles
podlesny's avatar
podlesny committed
  const auto& totalObstacles = nBodyAssembler_.totalObstacles_;
podlesny's avatar
podlesny committed
  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];
    }
  }

podlesny's avatar
podlesny committed
  //print(totalObstacles, "totalObstacles:");
podlesny's avatar
.  
podlesny committed

podlesny's avatar
podlesny committed
  //print(lower, "lower obstacles:");
  //print(upper, "upper obstacles:");
podlesny's avatar
.  
podlesny committed

podlesny's avatar
.  
podlesny committed
  // compute velocity obstacles
podlesny's avatar
podlesny committed
  /*Vector vLower, vUpper;
podlesny's avatar
.  
podlesny committed
  std::vector<Vector> u0, v0;
  updaters.rate_->extractOldVelocity(v0);
  updaters.rate_->extractOldDisplacement(u0);

  Vector totalu0, totalv0;
podlesny's avatar
podlesny committed
  nBodyAssembler_.concatenateVectors(u0, totalu0);
  nBodyAssembler_.concatenateVectors(v0, totalv0);
podlesny's avatar
.  
podlesny committed

  updaters.rate_->velocityObstacles(totalu0, lower, totalv0, vLower);
podlesny's avatar
podlesny committed
  updaters.rate_->velocityObstacles(totalu0, upper, totalv0, vUpper); */
podlesny's avatar
.  
podlesny committed

podlesny's avatar
podlesny committed
  //print(vLower, "vLower obstacles:");
  //print(vUpper, "vUpper obstacles:");
podlesny's avatar
.  
podlesny committed

podlesny's avatar
podlesny committed
  //std::cout << "- Problem assembled: success" << std::endl;
podlesny's avatar
.  
podlesny committed

podlesny's avatar
podlesny committed
  //print(ignoreNodes_, "ignoreNodes:");
podlesny's avatar
.  
podlesny committed

podlesny's avatar
podlesny committed
  // set up functional and TNMMG solver
podlesny's avatar
podlesny committed
  //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);

podlesny's avatar
podlesny committed
  auto step = solverFactory.step();

podlesny's avatar
podlesny committed
  //std::cout << "- Functional and TNNMG step setup: success" << std::endl;
podlesny's avatar
.  
podlesny committed

podlesny's avatar
podlesny committed
  EnergyNorm<Matrix, Vector> energyNorm(bilinearForm);
podlesny's avatar
.  
podlesny committed
  LoopSolver<Vector> velocityProblemSolver(*step.get(), velocityMaxIterations_,
                                           velocityTolerance_, energyNorm,
podlesny's avatar
podlesny committed
                                           verbosity_);
podlesny's avatar
.  
podlesny committed

  size_t fixedPointIteration;
  size_t multigridIterations = 0;
podlesny's avatar
.  
podlesny committed
  std::vector<ScalarVector> alpha(nBodies);
  updaters.state_->extractAlpha(alpha);
podlesny's avatar
podlesny committed

  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_;
       ++fixedPointIteration) {
podlesny's avatar
.  
podlesny committed

podlesny's avatar
podlesny committed
    //print(alpha, "alpha:");
podlesny's avatar
.  
podlesny committed

podlesny's avatar
.  
podlesny committed
    // contribution from nonlinearity
podlesny's avatar
.  
podlesny committed
    globalFriction_.updateAlpha(alpha);

podlesny's avatar
.  
podlesny committed
    //print(velocityIterates, "velocityIterates:");
    //print(totalVelocityIterate, "totalVelocityIterate:");
podlesny's avatar
podlesny committed
    //std::cout << "- FixedPointIteration iterate" << std::endl;
podlesny's avatar
.  
podlesny committed

podlesny's avatar
.  
podlesny committed
    // solve a velocity problem
podlesny's avatar
.  
podlesny committed
    solverFactory.setProblem(totalVelocityIterate);
podlesny's avatar
.  
podlesny committed

podlesny's avatar
podlesny committed
    //std::cout << "- Velocity problem set" << std::endl;
podlesny's avatar
.  
podlesny committed

    velocityProblemSolver.preprocess();
podlesny's avatar
podlesny committed
    //std::cout << "-- Preprocessed" << std::endl;
    velocityProblemSolver.solve();
podlesny's avatar
podlesny committed
    //std::cout << "-- Solved" << std::endl;
podlesny's avatar
.  
podlesny committed

    const auto& tnnmgSol = step->getSol();
podlesny's avatar
podlesny committed

podlesny's avatar
podlesny committed
    //std::cout << "FixPointIterator: Energy of TNNMG solution: " << J(tnnmgSol) << std::endl;
podlesny's avatar
podlesny committed

podlesny's avatar
podlesny committed
    nBodyAssembler_.postprocess(tnnmgSol, velocityIterates);
    //nBodyAssembler_.postprocess(totalVelocityIterate, velocityIterates);
podlesny's avatar
podlesny committed
    //print(totalVelocityIterate, "totalVelocityIterate:");
    //print(velocityIterates, "velocityIterates:");
podlesny's avatar
.  
podlesny committed

    //DUNE_THROW(Dune::Exception, "Just need to stop here!");
podlesny's avatar
.  
podlesny committed

    multigridIterations += velocityProblemSolver.getResult().iterations;

podlesny's avatar
podlesny committed
    Vector v_m = old_v;
    v_m *= 1.0 - lambda_;
    Dune::MatrixVector::addProduct(v_m, lambda_, tnnmgSol);
podlesny's avatar
.  
podlesny committed

podlesny's avatar
.  
podlesny committed
    // extract relative velocities in mortar basis
    std::vector<Vector> v_rel;
podlesny's avatar
podlesny committed
    relativeVelocities(v_m, v_rel);

    //print(v_m, "v_m: ");
podlesny's avatar
.  
podlesny committed
    //print(v_rel, "v_rel");

podlesny's avatar
podlesny committed
    //std::cout << "- State problem set" << std::endl;
podlesny's avatar
.  
podlesny committed

    // solve a state problem
podlesny's avatar
.  
podlesny committed
    updaters.state_->solve(v_rel);
podlesny's avatar
.  
podlesny committed

podlesny's avatar
podlesny committed
    //std::cout << "-- Solved" << std::endl;
podlesny's avatar
.  
podlesny committed

podlesny's avatar
.  
podlesny committed
    std::vector<ScalarVector> newAlpha(nBodies);
podlesny's avatar
.  
podlesny committed
    updaters.state_->extractAlpha(newAlpha);

podlesny's avatar
podlesny committed
    //print(newAlpha, "new alpha:");

podlesny's avatar
.  
podlesny committed
    bool breakCriterion = true;
podlesny's avatar
.  
podlesny committed
    for (int i=0; i<nBodies; i++) {
podlesny's avatar
.  
podlesny committed
        if (alpha[i].size()==0 || newAlpha[i].size()==0)
            continue;

podlesny's avatar
podlesny committed
        //print(alpha[i], "alpha i:");
        //print(newAlpha[i], "new alpha i:");
podlesny's avatar
.  
podlesny committed
        if (errorNorms_[i]->diff(alpha[i], newAlpha[i]) >= fixedPointTolerance_) {
podlesny's avatar
.  
podlesny committed
            breakCriterion = false;
podlesny's avatar
podlesny committed
            //std::cout << "fixedPoint error: " << errorNorms_[i]->diff(alpha[i], newAlpha[i]) << std::endl;
podlesny's avatar
.  
podlesny committed
            break;
        }
    }
podlesny's avatar
.  
podlesny committed
    if (lambda_ < 1e-12 or breakCriterion) {
podlesny's avatar
podlesny committed
      //std::cout << "-FixedPointIteration finished! " << (lambda_ < 1e-12 ? "lambda" : "breakCriterion") << std::endl;
      fixedPointIteration++;
podlesny's avatar
.  
podlesny committed
    alpha = newAlpha;
podlesny's avatar
.  
podlesny committed

podlesny's avatar
podlesny committed
  //TODO: recently added, might be wrong or superfluous
  globalFriction_.updateAlpha(alpha);

  //print(alpha, "alpha: ");

  //std::cout << "-FixedPointIteration finished! " << std::endl;
podlesny's avatar
.  
podlesny committed

  if (fixedPointIteration == fixedPointMaxIterations_)
    DUNE_THROW(Dune::Exception, "FPI failed to converge");

podlesny's avatar
.  
podlesny committed
  updaters.rate_->postProcess(velocityIterates);
  // 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;
podlesny's avatar
podlesny committed
/*std::ostream &operator<<(std::ostream &stream,
                         FixedPointIterationCounter const &fpic) {
  return stream << "(" << fpic.iterations << "," << fpic.multigridIterations
                << ")";
podlesny's avatar
podlesny committed
}*/
podlesny's avatar
podlesny committed
template <class Factory, class NBodyAssembler, class Updaters, class ErrorNorms>
void FixedPointIterator<Factory, NBodyAssembler, Updaters, ErrorNorms>::relativeVelocities(
podlesny's avatar
.  
podlesny committed
    const Vector& v,
    std::vector<Vector>& v_rel) const {
podlesny's avatar
.  
podlesny committed

podlesny's avatar
podlesny committed
    const size_t nBodies = nBodyAssembler_.nGrids();
podlesny's avatar
.  
podlesny committed
   // const auto& contactCouplings = nBodyAssembler.getContactCouplings();
podlesny's avatar
.  
podlesny committed

podlesny's avatar
.  
podlesny committed
    size_t globalIdx = 0;
podlesny's avatar
.  
podlesny committed
    v_rel.resize(nBodies);
podlesny's avatar
.  
podlesny committed
    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++;
        }
    }

podlesny's avatar
.  
podlesny committed

podlesny's avatar
.  
podlesny committed
   /*
podlesny's avatar
.  
podlesny committed
        boundaryNodes =
podlesny's avatar
.  
podlesny committed

podlesny's avatar
.  
podlesny committed
        const auto gridView = contactCouplings[couplingIndices[0]]->nonmortarBoundary().gridView();
podlesny's avatar
.  
podlesny committed

podlesny's avatar
.  
podlesny committed
        Dune::MultipleCodimMultipleGeomTypeMapper<
            decltype(gridView), Dune::MCMGVertexLayout> const vertexMapper(gridView, Dune::mcmgVertexLayout());
podlesny's avatar
.  
podlesny committed

podlesny's avatar
.  
podlesny committed
        for (auto it = gridView.template begin<block_size>(); it != gridView.template end<block_size>(); ++it) {
            const auto i = vertexMapper.index(*it);
podlesny's avatar
.  
podlesny committed

podlesny's avatar
.  
podlesny committed
            for (size_t j=0; j<couplingIndices.size(); j++) {
                const auto couplingIdx = couplingIndices[j];
podlesny's avatar
.  
podlesny committed

podlesny's avatar
.  
podlesny committed
                if (not contactCouplings[couplingIdx]->nonmortarBoundary().containsVertex(i))
                  continue;
podlesny's avatar
.  
podlesny committed

podlesny's avatar
.  
podlesny committed
                localToGlobal_.emplace_back(i);
                restrictions_.emplace_back(weights[bodyIdx][i], weightedNormalStress[bodyIdx][i],
                                          couplings[i]->frictionData()(geoToPoint(it->geometry())));
                break;
podlesny's avatar
.  
podlesny committed
            }

podlesny's avatar
.  
podlesny committed
            globalIdx++;
podlesny's avatar
.  
podlesny committed
        }
podlesny's avatar
.  
podlesny committed
        maxIndex_[bodyIdx] = globalIdx;
    }*/
podlesny's avatar
.  
podlesny committed
}

#include "fixedpointiterator_tmpl.cc"