fixedpointiterator.cc

#ifdef HAVE_CONFIG_H
#include "config.h"
#endif

#include <dune/common/exceptions.hh>

#include <dune/solvers/common/arithmetic.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 "../enums.hh"
#include "../enumparser.hh"

#include "fixedpointiterator.hh"

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

template <class Factory, class Updaters, class ErrorNorm>
FixedPointIterator<Factory, Updaters, ErrorNorm>::FixedPointIterator(
    Factory &factory, Dune::ParameterTree const &parset,
    std::shared_ptr<Nonlinearity> globalFriction, ErrorNorm const &errorNorm)
    : step_(factory.getStep()),
      parset_(parset),
      globalFriction_(globalFriction),
      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")),
      errorNorm_(errorNorm) {}

template <class Factory, class Updaters, class ErrorNorm>
FixedPointIterationCounter
FixedPointIterator<Factory, Updaters, ErrorNorm>::run(
    Updaters updaters, Matrix const &velocityMatrix, Vector const &velocityRHS,
    Vector &velocityIterate) {
  EnergyNorm<Matrix, Vector> energyNorm(velocityMatrix);
  LoopSolver<Vector> velocityProblemSolver(step_.get(), velocityMaxIterations_,
                                           velocityTolerance_, &energyNorm,
                                           verbosity_, false); // absolute error

  size_t fixedPointIteration;
  size_t multigridIterations = 0;
  std::vector<ScalarVector> alpha;
  updaters.state_->extractAlpha(alpha);
  for (fixedPointIteration = 0; fixedPointIteration < fixedPointMaxIterations_;
       ++fixedPointIteration) {
    // solve a velocity problem
   // globalFriction_->updateAlpha(alpha);
    ConvexProblem convexProblem(1.0, velocityMatrix, *globalFriction_,
                                velocityRHS, velocityIterate);
    BlockProblem velocityProblem(parset_, convexProblem);
    step_->setProblem(velocityIterate, velocityProblem);
    //step_->setProblem(velocityIterate);
    velocityProblemSolver.preprocess();
    velocityProblemSolver.solve();

    multigridIterations += velocityProblemSolver.getResult().iterations;

    std::vector<Vector> v_m;
    updaters.rate_->extractOldVelocity(v_m);

    for (size_t i=0; i<v_m.size(); i++) {
      v_m[i] *= 1.0 - lambda_;
      //Arithmetic::addProduct(v_m[i], lambda_, velocityIterate[i]);
    }

    // compute relative velocities on contact boundaries
    relativeVelocities(v_m);

    // solve a state problem
    updaters.state_->solve(v_m);
    ScalarVector newAlpha;
   /* updaters.state_->extractAlpha(newAlpha);

    if (lambda_ < 1e-12 or
        errorNorm_.diff(alpha, newAlpha) < fixedPointTolerance_) {
      fixedPointIteration++;
      break;
    }
    alpha = newAlpha;*/
  }
  if (fixedPointIteration == fixedPointMaxIterations_)
    DUNE_THROW(Dune::Exception, "FPI failed to converge");

  updaters.rate_->postProcess(velocityIterate);

  // 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 Updaters, class ErrorNorm>
void FixedPointIterator<Factory, Updaters, ErrorNorm>::relativeVelocities(std::vector<Vector>& v_m) const {
  using field_type = typename Factory::Matrix::field_type;

  // adaptation of DualMortarCoupling::setup()

  const size_t dim = DeformedGrid::dimension;
  typedef typename DeformedGrid::LeafGridView GridView;

  //cache of local bases
  typedef Dune::PQkLocalFiniteElementCache<typename DeformedGrid::ctype, field_type, dim,1> FiniteElementCache1;
  FiniteElementCache1 cache1;

  // cache for the dual functions on the boundary
  using DualCache = Dune::Contact::DualBasisAdapter<GridView, field_type>;
  std::unique_ptr<DualCache> dualCache;
  dualCache = std::make_unique< Dune::Contact::DualBasisAdapterGlobal<GridView, field_type> >();
/*
  // define FE grid functions
  std::vector<BasisGridFunction<> > gridFunctions(nBodyAssembler_.nGrids());
  for (size_t i=0; i<gridFunctions.size(); i++) {

  }

  const auto& contactCouplings = nBodyAssembler_.getContactCouplings();
  for (size_t i=0; i<contactCouplings.size(); i++) {
    auto contactCoupling = contactCouplings[i];
    auto glue = contactCoupling->getGlue();

    // loop over all intersections
    for (const auto& rIs : intersections(glue)) {
        const auto& inside = rIs.inside();

        if (!nonmortarBoundary_.contains(rIs.inside(),rIs.indexInInside()))
            continue;

        const auto& outside = rIs.outside();

        // types of the elements supporting the boundary segments in question
        Dune::GeometryType nonmortarEType = inside.type();
        Dune::GeometryType mortarEType    = outside.type();

        const auto& domainRefElement = Dune::ReferenceElements<ctype, dim>::general(nonmortarEType);
        const auto& targetRefElement = Dune::ReferenceElements<ctype, dim>::general(mortarEType);

        int noOfMortarVec = targetRefElement.size(dim);

        Dune::GeometryType nmFaceType = domainRefElement.type(rIs.indexInInside(),1);
        Dune::GeometryType mFaceType  = targetRefElement.type(rIs.indexInOutside(),1);

        // Select a quadrature rule
        // 2 in 2d and for integration over triangles in 3d.  If one (or both) of the two faces involved
        // are quadrilaterals, then the quad order has to be risen to 3 (4).
        int quadOrder = 2 + (!nmFaceType.isSimplex()) + (!mFaceType.isSimplex());
        const auto& quadRule = Dune::QuadratureRules<ctype, dim-1>::rule(rIs.type(), quadOrder);

        const auto& mortarFiniteElement = cache1.get(mortarEType);
        dualCache->bind(inside, rIs.indexInInside());

        std::vector<Dune::FieldVector<field_type,1> > mortarQuadValues, dualQuadValues;

        const auto& rGeom = rIs.geometry();
        const auto& rGeomOutside = rIs.geometryOutside();
        const auto& rGeomInInside = rIs.geometryInInside();
        const auto& rGeomInOutside = rIs.geometryInOutside();

        int nNonmortarFaceNodes = domainRefElement.size(rIs.indexInInside(),1,dim);
        std::vector<int> nonmortarFaceNodes;
        for (int i=0; i<nNonmortarFaceNodes; i++) {
          int faceIdxi = domainRefElement.subEntity(rIs.indexInInside(), 1, i, dim);
          nonmortarFaceNodes.push_back(faceIdxi);
        }

        for (const auto& quadPt : quadRule) {

            // compute integration element of overlap
            ctype integrationElement = rGeom.integrationElement(quadPt.position());

            // quadrature point positions on the reference element
            Dune::FieldVector<ctype,dim> nonmortarQuadPos = rGeomInInside.global(quadPt.position());
            Dune::FieldVector<ctype,dim> mortarQuadPos    = rGeomInOutside.global(quadPt.position());

            // The current quadrature point in world coordinates
            Dune::FieldVector<field_type,dim> nonmortarQpWorld = rGeom.global(quadPt.position());
            Dune::FieldVector<field_type,dim> mortarQpWorld    = rGeomOutside.global(quadPt.position());;

            // the gap direction (normal * gapValue)
            Dune::FieldVector<field_type,dim> gapVector = mortarQpWorld  - nonmortarQpWorld;

            //evaluate all shapefunctions at the quadrature point
            //nonmortarFiniteElement.localBasis().evaluateFunction(nonmortarQuadPos,nonmortarQuadValues);
            mortarFiniteElement.localBasis().evaluateFunction(mortarQuadPos,mortarQuadValues);
            dualCache->evaluateFunction(nonmortarQuadPos,dualQuadValues);

            // loop over all Lagrange multiplier shape functions
            for (int j=0; j<nNonmortarFaceNodes; j++) {

                int globalDomainIdx = indexSet0.subIndex(inside,nonmortarFaceNodes[j],dim);
                int rowIdx = globalToLocal[globalDomainIdx];

                weakObstacle_[rowIdx][0] += integrationElement * quadPt.weight()
                    * dualQuadValues[nonmortarFaceNodes[j]] * (gapVector*avNormals[globalDomainIdx]);

                // loop over all mortar shape functions
                for (int k=0; k<noOfMortarVec; k++) {

                    int colIdx  = indexSet1.subIndex(outside, k, dim);
                    if (!mortarBoundary_.containsVertex(colIdx))
                        continue;

                    // Integrate over the product of two shape functions
                    field_type mortarEntry =  integrationElement* quadPt.weight()* dualQuadValues[nonmortarFaceNodes[j]]* mortarQuadValues[k];

                    Dune::MatrixVector::addToDiagonal(mortarLagrangeMatrix_[rowIdx][colIdx], mortarEntry);

                }

            }

        }






  }

      ///////////////////////////////////
      //  reducing nonmortar boundary
      /////////////////////////////////

      // Get all fine grid boundary segments that are totally covered by the grid-glue segments
      typedef std::pair<int,int> Pair;
      std::map<Pair,ctype> coveredArea, fullArea;

      // initialize with area of boundary faces
      for (const auto& bIt : nonmortarBoundary_) {
          const Pair p(indexSet0.index(bIt.inside()),bIt.indexInInside());
          fullArea[p] = bIt.geometry().volume();
          coveredArea[p] = 0;
      }

      // sum up the remote intersection areas to find out which are totally covered
      for (const auto& rIs : intersections(glue))
          coveredArea[Pair(indexSet0.index(rIs.inside()),rIs.indexInInside())] += rIs.geometry().volume();

      // add all fine grid faces that are totally covered by the contact mapping
      for (const auto& bIt : nonmortarBoundary_) {
          const auto& inside = bIt.inside();
          if(coveredArea[Pair(indexSet0.index(inside),bIt.indexInInside())]/
              fullArea[Pair(indexSet0.index(inside),bIt.indexInInside())] >= coveredArea_)
              boundaryWithMapping.addFace(inside, bIt.indexInInside());
      }

      //writeBoundary(boundaryWithMapping,debugPath_ + "relevantNonmortar");


      // \todo replace by all fine grid segments which are totally covered by the RemoteIntersections.
      //for (const auto& rIs : intersections(glue))
      //    boundaryWithMapping.addFace(rIs.inside(),rIs.indexInInside());

      printf("contact mapping could be built for %d of %d boundary segments.\n",
             boundaryWithMapping.numFaces(), nonmortarBoundary_.numFaces());

      nonmortarBoundary_ = boundaryWithMapping;
      mortarBoundary_.setup(gridView1);
      for (const auto& rIs : intersections(glue))
          if (nonmortarBoundary_.contains(rIs.inside(),rIs.indexInInside()))
              mortarBoundary_.addFace(rIs.outside(),rIs.indexInOutside());


      // Assemble the diagonal matrix coupling the nonmortar side with the lagrange multiplyers there
      assembleNonmortarLagrangeMatrix();

      // The weak obstacle vector
      weakObstacle_.resize(nonmortarBoundary_.numVertices());
      weakObstacle_ = 0;

      // ///////////////////////////////////////////////////////////
      //   Get the occupation structure for the mortar matrix
      // ///////////////////////////////////////////////////////////

      // todo Also restrict mortar indices and don't use the whole grid level.
      Dune::MatrixIndexSet mortarIndices(nonmortarBoundary_.numVertices(), grid1_->size(dim));

      // Create mapping from the global set of block dofs to the ones on the contact boundary
      std::vector<int> globalToLocal;
      nonmortarBoundary_.makeGlobalToLocal(globalToLocal);

      // loop over all intersections
      for (const auto& rIs : intersections(glue)) {

          if (!nonmortarBoundary_.contains(rIs.inside(),rIs.indexInInside()))
              continue;

          const auto& inside = rIs.inside();
          const auto& outside = rIs.outside();

          const auto& domainRefElement = Dune::ReferenceElements<ctype, dim>::general(inside.type());
          const auto& targetRefElement = Dune::ReferenceElements<ctype, dim>::general(outside.type());

          int nDomainVertices = domainRefElement.size(dim);
          int nTargetVertices = targetRefElement.size(dim);

          for (int j=0; j<nDomainVertices; j++) {

              int localDomainIdx = globalToLocal[indexSet0.subIndex(inside,j,dim)];

              // if the vertex is not contained in the restricted contact boundary then dismiss it
              if (localDomainIdx == -1)
                  continue;

              for (int k=0; k<nTargetVertices; k++) {
                  int globalTargetIdx = indexSet1.subIndex(outside,k,dim);
                  if (!mortarBoundary_.containsVertex(globalTargetIdx))
                      continue;

                  mortarIndices.add(localDomainIdx, globalTargetIdx);
              }
          }
      }

      mortarIndices.exportIdx(mortarLagrangeMatrix_);

      // Clear it
      mortarLagrangeMatrix_ = 0;


      //cache of local bases
      FiniteElementCache1 cache1;

      std::unique_ptr<DualCache> dualCache;
      dualCache = std::make_unique< Dune::Contact::DualBasisAdapterGlobal<GridView0, field_type> >();

      std::vector<Dune::FieldVector<ctype,dim> > avNormals;
      avNormals = nonmortarBoundary_.getNormals();


      }

      // ///////////////////////////////////////
      //    Compute M = D^{-1} \hat{M}
      // ///////////////////////////////////////

      Dune::BCRSMatrix<MatrixBlock>& M  = mortarLagrangeMatrix_;
      Dune::BDMatrix<MatrixBlock>& D    = nonmortarLagrangeMatrix_;

      // First compute D^{-1}
      D.invert();

      // Then the matrix product D^{-1} \hat{M}
      for (auto rowIt = M.begin(); rowIt != M.end(); ++rowIt) {
          const auto rowIndex = rowIt.index();
          for (auto& entry : *rowIt)
              entry.leftmultiply(D[rowIndex][rowIndex]);
      }

      // weakObstacles in transformed basis = D^{-1}*weakObstacle_
      for(size_t rowIdx=0; rowIdx<weakObstacle_.size(); rowIdx++)
          weakObstacle_[rowIdx] *= D[rowIdx][rowIdx][0][0];

      gridGlueBackend_->clear(); */
}

#include "fixedpointiterator_tmpl.cc"