cantorfactory.cc

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

#include <dune/common/fvector.hh>

#include <dune/fufem/geometry/convexpolyhedron.hh>

#include <dune/contact/projections/normalprojection.hh>

#include "../utils/almostequal.hh"

#include "cantorfactory.hh"

#include "../multi-body-problem-data/bc.hh"
#include "../multi-body-problem-data/cuboidgeometry.hh"
#include "../multi-body-problem-data/myglobalfrictiondata.hh"
#include "../multi-body-problem-data/weakpatch.hh"
#include "../multi-body-problem-data/cubegridconstructor.hh"

#include "../utils/diameter.hh"

template <class GridType, int dim>
void CantorFactory<GridType, dim>::setBodies() {



    for (size_t level=0; level<=maxLevel_; level++) {
        const LevelCubes& levelCubes = cubes_[level];

        for (size_t i=0; i<levelCubes.size(); i++) {
            const std::shared_ptr<Cube>& cube = levelCubes[i];

            CubeGridConstructor
        }
    }

    std::vector<LocalVector> origins(this->bodyCount_);
    std::vector<LocalVector> lowerWeakOrigins(this->bodyCount_);
    std::vector<LocalVector> upperWeakOrigins(this->bodyCount_);

#if MY_DIM == 3
        double const depth = 0.60;

        origins[0] = {0,0,0};
        lowerWeakOrigins[0] = {0.2,0,0};
        upperWeakOrigins[0] = {0.2,width,0};

        cuboidGeometries_[0] = std::make_shared<CuboidGeometry>(origins[0], lowerWeakOrigins[0], upperWeakOrigins[0], 0.0, weakLength, length, width, depth);
        for (size_t i=1; i<this->bodyCount_-1; i++) {
            origins[i] = cuboidGeometries_[i-1]->D;
            lowerWeakOrigins[i] = {0.6,i*width,0};
            upperWeakOrigins[i] = {0.6,(i+1)*width,0};

            cuboidGeometries_[i] = std::make_shared<CuboidGeometry>(origins[i], lowerWeakOrigins[i], upperWeakOrigins[i], weakLength, weakLength, length, width, depth);
        }
        const size_t idx = this->bodyCount_-1;
        origins[idx] = cuboidGeometries_[idx-1]->D;
        lowerWeakOrigins[idx] = {0.6,idx*width,0};
        upperWeakOrigins[idx] = {0.6,(idx+1)*width,0};

        cuboidGeometries_[idx] = std::make_shared<CuboidGeometry>(origins[idx], lowerWeakOrigins[idx], upperWeakOrigins[idx], weakLength, 0.0, length, width, depth);
#elif MY_DIM == 2
        origins[0] = {0,0};
        lowerWeakOrigins[0] = {0.2,0};
        upperWeakOrigins[0] = {0.2,width};

        cuboidGeometries_[0] = std::make_shared<CuboidGeometry>(origins[0], lowerWeakOrigins[0], upperWeakOrigins[0], 0.0, weakLength, length, width);
        for (size_t i=1; i<this->bodyCount_-1; i++) {
            origins[i] = cuboidGeometries_[i-1]->D;
            lowerWeakOrigins[i] = {0.6,i*width};
            upperWeakOrigins[i] = {0.6,(i+1)*width};

            cuboidGeometries_[i] = std::make_shared<CuboidGeometry>(origins[i], lowerWeakOrigins[i], upperWeakOrigins[i], weakLength, weakLength, length, width);
        }
        const size_t idx = this->bodyCount_-1;
        origins[idx] = cuboidGeometries_[idx-1]->D;
        lowerWeakOrigins[idx] = {0.6,idx*width};
        upperWeakOrigins[idx] = {0.6,(idx+1)*width};

        cuboidGeometries_[idx] = std::make_shared<CuboidGeometry>(origins[idx], lowerWeakOrigins[idx], upperWeakOrigins[idx], weakLength, 0.0, length, width);
#else
#error CuboidGeometry only supports 2D and 3D!"
#endif

    // set up reference grids
    gridConstructor_ = new GridsConstructor<GridType>(cuboidGeometries_);
    auto& grids = gridConstructor_->getGrids();

    for (size_t i=0; i<this->bodyCount_; i++) {
        const auto& cuboidGeometry = *cuboidGeometries_[i];

        // define weak patch and refine grid
        auto lowerWeakPatch = lowerWeakPatches_[i] = std::make_shared<ConvexPolyhedron<LocalVector>>();
        auto upperWeakPatch = upperWeakPatches_[i] = std::make_shared<ConvexPolyhedron<LocalVector>>();
        getWeakPatch<LocalVector>(this->parset_.sub("boundary.friction.weakening"), cuboidGeometry, *lowerWeakPatch, *upperWeakPatch);
        refine(*grids[i], *lowerWeakPatch, this->parset_.template get<double>("boundary.friction.smallestDiameter"), cuboidGeometry.lengthScale);
        refine(*grids[i], *upperWeakPatch, this->parset_.template get<double>("boundary.friction.smallestDiameter"), cuboidGeometry.lengthScale);
    }

    for (size_t i=0; i<this->bodyCount_; i++) {
        this->bodies_[i] = std::make_shared<typename Base::Body>(bodyData_, grids[i]);
    }
}

template <class GridType, int dim>
void CantorFactory<GridType, dim>::setCouplings() {
    const auto deformedGrids = this->contactNetwork_.deformedGrids();

    for (size_t i=0; i<this->bodyCount_; i++) {
        const auto& cuboidGeometry = *cuboidGeometries_[i];
        leafFaces_[i] = std::make_shared<LeafFaces>(this->contactNetwork_.leafView(i), cuboidGeometry);
        levelFaces_[i] = std::make_shared<LevelFaces>(this->contactNetwork_.levelView(i), cuboidGeometry);
    }

    auto contactProjection = std::make_shared<Dune::Contact::NormalProjection<LeafBoundaryPatch>>();
    std::shared_ptr<typename Base::CouplingPair::BackEndType> backend = nullptr;
    for (size_t i=0; i<this->couplings_.size(); i++) {
      auto& coupling = this->couplings_[i];
      coupling = std::make_shared<typename Base::CouplingPair>();

      nonmortarPatch_[i] = std::make_shared<LevelBoundaryPatch>(levelFaces_[i]->upper);
      mortarPatch_[i] = std::make_shared<LevelBoundaryPatch>(levelFaces_[i+1]->lower);

      coupling->set(i, i+1, nonmortarPatch_[i], mortarPatch_[i], 0.1, Base::CouplingPair::CouplingType::STICK_SLIP, contactProjection, backend);
    }
}

template <class GridType, int dim>
void CantorFactory<GridType, dim>::setBoundaryConditions() {
    using LeafBoundaryCondition = BoundaryCondition<DeformedLeafGridView, dim>;
    using FrictionBoundaryCondition = FrictionBoundaryCondition<DeformedLeafGridView, LocalVector, dim>;

    using Function = Dune::VirtualFunction<double, double>;
    std::shared_ptr<Function> neumannFunction = std::make_shared<NeumannCondition>();
    std::shared_ptr<Function> velocityDirichletFunction = std::make_shared<VelocityDirichletCondition>();

    const double lengthScale = CuboidGeometry::lengthScale;

    for (size_t i=0; i<this->bodyCount_; i++) {
      const auto& body = this->contactNetwork_.body(i);
      const auto& leafVertexCount = body->leafVertexCount();

      std::cout << "Grid" << i << " Number of DOFs: " << leafVertexCount << std::endl;

      // friction boundary
      if (i<this->bodyCount_-1 && upperWeakPatches_[i]->vertices.size()>0) {
        std::shared_ptr<FrictionBoundaryCondition> frictionBoundary = std::make_shared<FrictionBoundaryCondition>();
        frictionBoundary->setBoundaryPatch(leafFaces_[i]->upper);
        frictionBoundary->setWeakeningPatch(upperWeakPatches_[i]);
        frictionBoundary->setFrictionData(std::make_shared<MyGlobalFrictionData<LocalVector>>(this->parset_.sub("boundary.friction"), *upperWeakPatches_[i], lengthScale));
        body->addBoundaryCondition(frictionBoundary);
      }
      if (i>0 && lowerWeakPatches_[i]->vertices.size()>0) {
        std::shared_ptr<FrictionBoundaryCondition> frictionBoundary = std::make_shared<FrictionBoundaryCondition>();
        frictionBoundary->setBoundaryPatch(leafFaces_[i]->lower);
        frictionBoundary->setWeakeningPatch(lowerWeakPatches_[i]);
        frictionBoundary->setFrictionData(std::make_shared<MyGlobalFrictionData<LocalVector>>(this->parset_.sub("boundary.friction"), *lowerWeakPatches_[i], lengthScale));
        body->addBoundaryCondition(frictionBoundary);
      }

      // Neumann boundary
      std::shared_ptr<LeafBoundaryCondition> neumannBoundary = std::make_shared<LeafBoundaryCondition>(std::make_shared<LeafBoundaryPatch>(body->leafView()), neumannFunction, "neumann");
      body->addBoundaryCondition(neumannBoundary);

      // Dirichlet Boundary
      // identify Dirichlet nodes on leaf level
      std::shared_ptr<Dune::BitSetVector<dim>> dirichletNodes = std::make_shared<Dune::BitSetVector<dim>>(leafVertexCount);
      for (int j=0; j<leafVertexCount; j++) {
        if (leafFaces_[i]->right.containsVertex(j))
          (*dirichletNodes)[j][0] = true;

        if (leafFaces_[i]->lower.containsVertex(j))
          (*dirichletNodes)[j][0] = true;

        #if MY_DIM == 3
        if (leafFaces_[i]->front.containsVertex(j) || leafFaces_[i]->back.containsVertex(j))
            dirichletNodes->at(j)[2] = true;
        #endif
      }

      std::shared_ptr<LeafBoundaryCondition> dirichletBoundary = std::make_shared<LeafBoundaryCondition>("dirichlet");
      dirichletBoundary->setBoundaryPatch(body->leafView(), dirichletNodes);
      dirichletBoundary->setBoundaryFunction(velocityDirichletFunction);
      body->addBoundaryCondition(dirichletBoundary);
    }
}


class CantorIndexHierarchy {
    public:
        typedef std::map<std::pair<int,int>, bool> LevelCantorIndices;

    private:
        const size_t maxLevel_;
        std::vector<size_t> levelN_;
        std::vector<LevelCantorIndices> cantorIndexHierarchy_;

        void prolongIndices(size_t currentLevel, size_t newLevel) {
            const LevelCantorIndices& indices = cantorIndexHierarchy_[currentLevel];
            LevelCantorIndices& newIndices = cantorIndexHierarchy_[newLevel];

            size_t offset = levelN_[currentLevel];

            std::map<std::pair<int,int>, bool>::const_iterator it = indices.begin();
            std::map<std::pair<int,int>, bool>::const_iterator endIt = indices.end();
            for (; it!=endIt; ++it) {
                int xID = it->first.first;
                int yID = it->first.second;
                newIndices[std::make_pair(xID, yID)] = true;
                newIndices[std::make_pair(xID+offset, yID)] = true;
                newIndices[std::make_pair(xID, yID+offset)] = true;
            }

            size_t doubleOffset = 2*offset;
            for (size_t i=offset+1; i<=doubleOffset; i=i+2) {
                newIndices[std::make_pair(doubleOffset, i)] = true;
                newIndices[std::make_pair(i, doubleOffset)] = true;
            }
        }

        void print(const LevelCantorIndices& indices) const {
            std::cout << "LevelCantorIndices: " << std::endl;

            std::map<std::pair<int,int>, bool>::const_iterator it = indices.begin();
            std::map<std::pair<int,int>, bool>::const_iterator endIt = indices.end();
            for (; it!=endIt; ++it) {
                std::cout << "(" << it->first.first << ", " << it->first.second << ")"<< std::endl;
            }
        }

    public:
        CantorIndexHierarchy(size_t maxLevel) :
            maxLevel_(maxLevel)
        {
            levelN_.resize(maxLevel_+1);
            cantorIndexHierarchy_.resize(maxLevel_+1);

            // init levelCantorIndices on level 0
            LevelCantorIndices& initCantorIndices = cantorIndexHierarchy_[0];
            initCantorIndices[std::make_pair(0, 1)] = true;
            initCantorIndices[std::make_pair(1, 0)] = true;
            initCantorIndices[std::make_pair(1, 2)] = true;
            initCantorIndices[std::make_pair(2, 1)] = true;

            for (size_t i=0; i<maxLevel_; i++) {
                levelN_[i] = std::pow(2, i+1);
                prolongIndices(i, i+1);
            }

            levelN_[maxLevel_] = std::pow(2, maxLevel_+1);
        }

        const LevelCantorIndices& levelCantorIndices(size_t i) const {
                return cantorIndexHierarchy_[i];
        }

        size_t levelN(const size_t i) const {
            return levelN_[i];
        }
};


template <class GridType>
class CantorFaultFactory
{
    //! Parameter for mapper class
    template<int dim>
    struct FaceMapperLayout
    {
        bool contains (Dune::GeometryType gt)
        {
            return gt.dim() == dim-1;
        }
    };

    protected:
        typedef OscUnitCube<GridType, 2> GridOb;
        static const int dimworld = GridType::dimensionworld;
        static const int dim = GridType::dimension;
        typedef typename GridType::ctype ctype;
        typedef typename GridType::LevelGridView GV;

        typedef typename Dune::MultipleCodimMultipleGeomTypeMapper<GV, FaceMapperLayout > FaceMapper;

    private:
        const std::map<size_t, size_t>& levelToCantorLevel_;
        const int coarseResolution_;
        const size_t maxLevel_;
        const size_t maxCantorLevel_;
        const CantorIndexHierarchy cantorIndexHierarchy_;
        const int coarseGridN_;

        std::shared_ptr<GridType> grid_;
        std::vector<double> levelResolutions_;

        std::shared_ptr<InterfaceNetwork<GridType>> interfaceNetwork_;

private:


    bool computeID(Dune::FieldVector<ctype, dimworld> vertex, const size_t gridN, std::pair<size_t, size_t>& IDs) const {
        ctype xID = vertex[0]*gridN;
        ctype yID = vertex[1]*gridN;

        ctype x = 0;
        ctype y = 0;

        bool xIsInt = almost_equal(std::modf(xID, &x), 0.0, 2);
        bool yIsInt = almost_equal(std::modf(yID, &y), 0.0, 2);

        if (xIsInt && yIsInt) {
            IDs = std::make_pair((size_t) x, (size_t) y);
            return true;
        }

        if (xIsInt) {
            y = std::ceil(yID);
            if ((size_t) y % 2==0)
                y = y-1;
            IDs = std::make_pair((size_t) x, (size_t) y);
            return true;
        }

        if (yIsInt) {
            x = std::ceil(xID);
            if ((size_t) x % 2==0)
                x = x-1;
            IDs = std::make_pair((size_t) x, (size_t) y);
            return true;
        }

        return false;
    }

public:
        //setup
    CantorFaultFactory(const std::map<size_t, size_t>& levelToCantorLevel, const int coarseResolution, const size_t maxLevel, const size_t maxCantorLevel) :
        levelToCantorLevel_(levelToCantorLevel),
        coarseResolution_(coarseResolution),
        maxLevel_(maxLevel),
        maxCantorLevel_(maxCantorLevel),
        cantorIndexHierarchy_(maxCantorLevel_),
        coarseGridN_(std::pow(2, coarseResolution_)),
        interfaceNetwork_(nullptr)
        {
            Dune::UGGrid<dim>::setDefaultHeapSize(4000);
            GridOb unitCube(coarseGridN_);
            grid_ = unitCube.grid();

            grid_->globalRefine(maxLevel_);

            levelResolutions_.resize(maxLevel_+1);

            // init interface network
            interfaceNetwork_ = std::make_shared<InterfaceNetwork<GridType>>(*grid_);

            for (size_t i=0; i<=maxLevel_; i++) {
                if (i>0) {
                    interfaceNetwork_->prolongLevel(i-1, i);
                }

                if (levelToCantorLevel_.count(i)) {

                    levelResolutions_[i] = std::pow(2, coarseResolution_+i);
                    const size_t levelResolution = levelResolutions_[i];
                    const size_t cantorResolution = cantorIndexHierarchy_.levelN(levelToCantorLevel_.at(i));

                    if (2*levelResolution<cantorResolution)
                        DUNE_THROW(Dune::Exception, "Level " + std::to_string(i) + " does not support cantor set with resolution " + std::to_string(cantorResolution));

                    const typename CantorIndexHierarchy::LevelCantorIndices& levelCantorIndices = cantorIndexHierarchy_.levelCantorIndices(levelToCantorLevel_.at(i));
                    std::shared_ptr<LevelInterfaceNetwork<GV>> levelInterfaceNetwork = interfaceNetwork_->levelInterfaceNetworkPtr(i);

                    const GV& gridView = levelInterfaceNetwork->levelGridView();
                    FaceMapper intersectionMapper(gridView);
                    std::vector<bool> intersectionHandled(intersectionMapper.size(),false);

                    for (const auto& elem:elements(gridView)) {
                        for (const auto& isect:intersections(gridView, elem)) {

                            if (intersectionHandled[intersectionMapper.subIndex(elem, isect.indexInInside(),1)])
                                continue;

                            intersectionHandled[intersectionMapper.subIndex(elem, isect.indexInInside(),1)]=true;

                            if (isect.boundary())
                                continue;

                            std::pair<size_t, size_t> intersectionID;
                            bool isAdmissibleID = computeID(isect.geometry().center(), cantorResolution, intersectionID);

                            if (isAdmissibleID && levelCantorIndices.count(intersectionID)) {
                                levelInterfaceNetwork->addIntersection(isect, i);
                            }
                        }
                    }
                }
            }
        }



#include "cantorfactory_tmpl.cc"