dualmortarcoupling.hh

// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set ts=8 sw=4 et sts=4:
#ifndef DUNE_TECTONIC_SPATIAL_SOLVING_CONTACT_DUAL_MORTAR_COUPLING_HH
#define DUNE_TECTONIC_SPATIAL_SOLVING_CONTACT_DUAL_MORTAR_COUPLING_HH

#include <memory>

#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/bdmatrix.hh>
#include <dune/istl/bvector.hh>

#include <dune/fufem/boundarypatch.hh>
#include <dune/fufem/boundarypatchprolongator.hh>

#include <dune/grid-glue/gridglue.hh>
#include <dune/grid-glue/extractors/codim1extractor.hh>
#include <dune/grid-glue/merging/merger.hh>

/** \brief Assembles the transfer operator for two-body contact.
 * Extension of Dune::Contact::DualMortarCoupling by crosspoints.
 * Note: so war only works in 2D!
 */
template<class field_type, class GridType0, class GridType1=GridType0>
class DualMortarCoupling {

    // /////////////////////
    // private types
    // /////////////////////

    enum {dim = GridType0::dimension};

    enum {dimworld = GridType0::dimensionworld};

  protected:
    typedef Dune::FieldMatrix<field_type, dim, dim> MatrixBlock;
    typedef Dune::BCRSMatrix<MatrixBlock> MatrixType;
    using ObstacleVector = Dune::BlockVector<Dune::FieldVector<field_type, 1> >;

    typedef typename GridType0::ctype ctype;

    typedef BoundaryPatch<typename GridType0::LeafGridView> LeafBoundaryPatch0;
    typedef BoundaryPatch<typename GridType1::LeafGridView> LeafBoundaryPatch1;

    typedef BoundaryPatch<typename GridType0::LevelGridView> LevelBoundaryPatch0;
    typedef BoundaryPatch<typename GridType1::LevelGridView> LevelBoundaryPatch1;

    using Extractor0 = Dune::GridGlue::Codim1Extractor<typename GridType0::LeafGridView>;
    using Extractor1 = Dune::GridGlue::Codim1Extractor<typename GridType1::LeafGridView>;

public:

    using Glue = Dune::GridGlue::GridGlue<Extractor0, Extractor1>;

    DualMortarCoupling(field_type overlap = 1e-2, field_type coveredArea = 1.0 - 1e-2)
        : overlap_(overlap), gridGlueBackend_(nullptr), coveredArea_(coveredArea)
    {}

    DualMortarCoupling(const GridType0& grid0, const GridType1& grid1,
                       field_type overlap = 1e-2, field_type coveredArea = 1.0 - 1e-2)
        : grid0_(&grid0), grid1_(&grid1), overlap_(overlap),
          gridGlueBackend_(nullptr), coveredArea_(coveredArea)
    {}

    virtual ~DualMortarCoupling() {
        /* Nothing. */
    }

    void reducePatches();

    void setCrosspoints(const std::set<size_t>& allCrosspoints);

    /** \brief Sets up the contact coupling operator on the grid leaf level */
    virtual void setup();

    /** \brief Assemble nonmortar matrix on the leaf level. */
    void assembleNonmortarLagrangeMatrix();

    /** \brief Get the dual mortar coupling matrix of the leaf level. */
    virtual const MatrixType& mortarLagrangeMatrix() const {
        return mortarLagrangeMatrix_;
    }

    /** \brief Get the dual mortar coupling matrix of the leaf level. */
    virtual MatrixType& mortarLagrangeMatrix() {
        return mortarLagrangeMatrix_;
    }

    /** \brief Get the non-mortar matrix. */
    const Dune::BDMatrix<MatrixBlock>& nonmortarLagrangeMatrix() const {
        return nonmortarLagrangeMatrix_;
    }

    /** \brief Setup leaf nonmortar and mortar patches. */
    virtual void setupContactPatch(const LevelBoundaryPatch0& coarseNonmortar,
                           const LevelBoundaryPatch1& coarseMortar) {

        BoundaryPatchProlongator<GridType0>::prolong(coarseNonmortar,nonmortarBoundary_);
        BoundaryPatchProlongator<GridType1>::prolong(coarseMortar,mortarBoundary_);
    }

    /** \brief Return the leafnonmortar boundary. */
    const LeafBoundaryPatch0& nonmortarBoundary() const {
        return nonmortarBoundary_;
    }

    /** \brief Return the leafmortar boundary. */
    const LeafBoundaryPatch1& mortarBoundary() const {
        return mortarBoundary_;
    }

    /** \brief Set the non-mortar and mortar grid. */
    void setGrids(const GridType0& grid0, const GridType1& grid1) {
        grid0_ = &grid0;
        grid1_ = &grid1;
    }

    /** \brief Get non-mortar grid. */
    const GridType0* getGrid0() const {
        return grid0_;
    }

    /** \brief Get mortar grid. */
    const GridType1* getGrid1() const {
        return grid1_;
    }

    /** \brief Set the percentage of covered area each face must at least have
     *         by the grid-glue projection.
     */
    void setCoveredArea(field_type coveredArea) {
        coveredArea_ = coveredArea;
    }

    void setOverlap(field_type overlap) {
        overlap_ = overlap;
    }

    /** \brief Get the GridGlue. */
    std::shared_ptr<Glue> getGlue() {
        return glue_;
    }

    const ObstacleVector& getWeakObstacle() const {
        return weakObstacle_;
    }

    const std::vector<int>& globalToLocal() const {
        return globalToLocal_;
    }

   /* const auto& crosspoints() const {
        return crosspoints_;
    }*/

    // /////////////////////////////////////////////
    //   Data members
    // /////////////////////////////////////////////

    /** \brief The two grids involved in the two-body contact problem. */
    const GridType0* grid0_;
    const GridType1* grid1_;

    /** \brief For each dof a bit specifying whether the dof carries an obstacle or not. */
    LeafBoundaryPatch0 nonmortarBoundary_;

    /** \brief The mortar boundary. */
    LeafBoundaryPatch1 mortarBoundary_;

    /** \brief The crosspoints of the nonmortar body. */
    std::set<size_t> crosspoints_;

    /** \brief Nonmortar boundary patch global to local map. */
    std::vector<int> globalToLocal_;

    /** \brief The matrix coupling mortar side and Lagrange multipliers. */
    MatrixType mortarLagrangeMatrix_;

    /** \brief The diagonal matrix coupling nonmortar side and Lagrange multipliers */
    Dune::BDMatrix<MatrixBlock> nonmortarLagrangeMatrix_;

    /** \brief The weak obstacles. */
    ObstacleVector weakObstacle_;

    /** \brief Allow some small penetration of the bodies. */
    field_type overlap_;

    /** \brief A grid-glue merger. If this is not set ContactMerge is used. */
    std::shared_ptr< Dune::GridGlue::Merger<field_type,dim-1,dim-1,dim> > gridGlueBackend_;

    // Path where to write the merged grids and the reduced contact boundaries
    //std::string debugPath_;

    std::shared_ptr<Glue> glue_;

    /** \brief Dismiss all faces that are not at least covered by the grid-glue projection for this
     *         much percentage ranging between one - for total coverage and zero for taking all faces.
     */
    field_type coveredArea_;
};

#include "dualmortarcoupling.cc"

#endif