dualmortarcoupling.cc

// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set ts=8 sw=4 et sts=4:
#include <dune/common/exceptions.hh>

#include <dune/istl/matrixindexset.hh>
#include <dune/istl/io.hh>

#include <dune/geometry/quadraturerules.hh>
#include <dune/geometry/type.hh>
#include <dune/geometry/referenceelements.hh>

#include <dune/localfunctions/lagrange/pqkfactory.hh>


#include <dune/grid-glue/merging/contactmerge.hh>
//#include <dune/grid-glue/adapter/gridgluevtkwriter.hh>
//#include <dune/grid-glue/adapter/gridglueamirawriter.hh>

#include <dune/contact/common/dualbasisadapter.hh>
#include <dune/matrix-vector/addtodiagonal.hh>

template <class field_type, class GridType0, class GridType1>
void DualMortarCoupling<field_type, GridType0, GridType1>::reducePatches() {
    typedef typename GridType0::LeafGridView GridView0;
    typedef typename GridType1::LeafGridView GridView1;

    using Element0 = typename GridView0::template Codim<0>::Entity;
    using Element1 = typename GridView1::template Codim<0>::Entity;

    auto desc0 = [&] (const Element0& e, unsigned int face) {
        return nonmortarBoundary_.contains(e,face);
    };

    auto desc1 = [&] (const Element1& e, unsigned int face) {
        return mortarBoundary_.contains(e,face);
    };

    GridView0 gridView0 = grid0_->leafGridView();
    GridView1 gridView1 = grid1_->leafGridView();

    auto extract0 = std::make_shared<Extractor0>(gridView0,desc0);
    auto extract1 = std::make_shared<Extractor1>(gridView1,desc1);

    if (!gridGlueBackend_)
        gridGlueBackend_ = std::make_shared< Dune::GridGlue::ContactMerge<dimworld, ctype> >(overlap_);
    glue_ = std::make_shared<Glue>(extract0, extract1, gridGlueBackend_);
    auto& glue = *glue_;
    glue.build();

    //std::cout << glue.size() << " remote intersections found." << std::endl;
    //GridGlueAmiraWriter::write<GlueType>(glue,debugPath_);

    // Restrict the hasObstacle fields to the part of the nonmortar boundary
    // where we could actually create a contact mapping
    LeafBoundaryPatch0 boundaryWithMapping(gridView0);

    const auto& indexSet0 = gridView0.indexSet();
    const auto& indexSet1 = gridView1.indexSet();

    ///////////////////////////////////
    //  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());
}

template <class field_type, class GridType0, class GridType1>
void DualMortarCoupling<field_type, GridType0, GridType1>::setCrosspoints(const std::set<size_t>& allCrosspoints) {
    crosspoints_.clear();

    if (!allCrosspoints.empty()) {
        Dune::BitSetVector<1> nonmortarPatchBoundary;
        nonmortarBoundary_.getPatchBoundaryVertices(nonmortarPatchBoundary);

        for (auto p : allCrosspoints) {
            if (nonmortarPatchBoundary[p][0]) {
                crosspoints_.emplace(p);
            }
        }
    }
}

template <class field_type, class GridType0, class GridType1>
void DualMortarCoupling<field_type, GridType0, GridType1>::assembleNonmortarLagrangeMatrix()
{  
    // clear matrix
    const auto nNonmortarVertices = nonmortarBoundary_.numVertices(); // - crosspoints_.size();
    nonmortarLagrangeMatrix_ = Dune::BDMatrix<MatrixBlock>(nNonmortarVertices);
    nonmortarLagrangeMatrix_ = 0;

    const auto& indexSet = grid0_->leafIndexSet();

    // loop over all faces of the boundary patch
    for (const auto& nIt : nonmortarBoundary_) {

        const auto& geometry = nIt.geometry();
        const field_type numCorners = geometry.corners();

        ctype intElem = geometry.integrationElement(Dune::FieldVector<ctype,dim-1>(0));

        field_type sfI = (numCorners==3) ? intElem/6.0 : intElem/numCorners;

        // turn scalar element mass matrix into vector-valued one
        // and add element mass matrix to the global matrix

        // Get global node ids
        const auto& inside = nIt.inside();
        const auto& refElement = Dune::ReferenceElements<ctype, dim>::general(inside.type());

        // Add localMatrix to nonmortarLagrangeMat
        for (int i=0; i<refElement.size(nIt.indexInInside(),1,dim); i++) {
            // we can use subEntity here because we add all indices anyway
            int v = globalToLocal_[indexSet.subIndex(inside,refElement.subEntity(nIt.indexInInside(),1,
                                                                          i,dim),dim)];
            if (v == -1)
               continue;

            Dune::MatrixVector::addToDiagonal(nonmortarLagrangeMatrix_[v][v],sfI);
        }

    }

}

template <class field_type, class GridType0, class GridType1>
void DualMortarCoupling<field_type, GridType0,GridType1>::setup()
{
    typedef typename GridType0::LeafGridView GridView0;
    typedef typename GridType1::LeafGridView GridView1;

    using FiniteElementCache1 = Dune::PQkLocalFiniteElementCache<typename GridType1::ctype, field_type, GridType1::dimension,1>;

    // cache for the dual functions on the boundary
    using DualCache = Dune::Contact::DualBasisAdapter<GridView0, field_type>;
    using ConstantDualCache = Dune::PQkLocalFiniteElementCache<typename GridType0::ctype, field_type, GridType0::dimension, 0>;

    GridView0 gridView0 = grid0_->leafGridView();
    GridView1 gridView1 = grid1_->leafGridView();

    const auto& indexSet0 = gridView0.indexSet();
    const auto& indexSet1 = gridView1.indexSet();

    // Create mapping from the global set of block dofs to the ones on the contact boundary
    const auto& vertices = *nonmortarBoundary_.getVertices();

    globalToLocal_.resize(vertices.size());
    hasObstacle_.resize(vertices.size());
    int idx = 0;
    for (size_t i=0; i<globalToLocal_.size(); i++) {
        bool nonmortarVertex = vertices[i][0]; // and crosspoints_.count(i)==0;
        globalToLocal_[i] = nonmortarVertex ? idx++ : -1;
        hasObstacle_[i][0] = nonmortarVertex ? true : false;
    }

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

    // The weak obstacle vector
    const auto nNonmortarVertices = nonmortarBoundary_.numVertices();// - crosspoints_.size();
    weakObstacle_.resize(nNonmortarVertices);
    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(nNonmortarVertices, grid1_->size(dim));

    // 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 or a crosspoint, 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> >();

    // adapt dual basis close to crosspoints using ConstantDualCache
    ConstantDualCache constantDualCache;

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

    // loop over all intersections and compute the matrix entriescrosspoints_
    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 typename FiniteElementCache0::FiniteElementType& nonmortarFiniteElement = cache0.get(nonmortarEType);
        const auto& mortarFiniteElement = cache1.get(mortarEType);

        const auto& constantDualFE = constantDualCache.get(nonmortarEType);
        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();

        bool isCrosspointFace = false;
        std::vector<int> nonmortarFaceNodes;
        for (int i=0; i<domainRefElement.size(rIs.indexInInside(),1,dim); i++) {
          int faceIdxi = domainRefElement.subEntity(rIs.indexInInside(), 1, i, dim);
          int globalIdx = indexSet0.subIndex(inside,faceIdxi,dim);

          /*if (crosspoints_.count(globalIdx)>0) {
            isCrosspointFace = true;
          } else {*/
            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);

            /*if (isCrosspointFace) {
              constantDualFE.localBasis().evaluateFunction(nonmortarQuadPos,dualQuadValues);

              // loop over all Lagrange multiplier shape functions
              for (size_t j=0; j<nonmortarFaceNodes.size(); j++) {

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

                  weakObstacle_[rowIdx][0] += integrationElement * quadPt.weight()
                      * dualQuadValues[0] * (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[0]* mortarQuadValues[k];

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

                  }

               }
            } else {*/
              dualCache->evaluateFunction(nonmortarQuadPos,dualQuadValues);

              // loop over all Lagrange multiplier shape functions
              for (size_t j=0; j<nonmortarFaceNodes.size(); 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);

                  }

               }
            //}
        }

    }

    // ///////////////////////////////////////
    //    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();
}