stackedblocksfactory.cc

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

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

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

#include "../problem-data/bc.hh"
#include "../problem-data/grid/cuboidgeometry.hh"
#include "../problem-data/myglobalfrictiondata.hh"

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

#include "stackedblocksfactory.hh"


template <class HostGridType, class VectorType>
void StackedBlocksFactory<HostGridType, VectorType>::setBodies() {

    // set up cuboid geometries
    double const length = 1.00;
    double const height = 0.27;
    double const weakLength = 0.60;

    std::vector<GlobalCoords> origins(this->bodyCount_);

    const auto& subParset = this->parset_.sub("boundary.friction.weakening");

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

        auto weakBound = [&] (const GlobalCoords& origin) {
            GlobalCoords h = {origin[0] + weakLength, origin[1], origin[2]};
            return h;
        };

        origins[0] = {0,0,0};
        GlobalCoords lowerWeakOrigin = {0.2, 0, 0};
        GlobalCoords upperWeakOrigin = {0.2, height, 0};

        cuboidGeometries_[0] = std::make_shared<CuboidGeometry>(origins[0], length, height, depth);
        cuboidGeometries_[0]->addWeakeningPatch(subParset, upperWeakOrigin, weakBound(upperWeakOrigin));

        for (size_t i=1; i<this->bodyCount_-1; i++) {
            origins[i] = cuboidGeometries_[i-1]->upperLeft();
            GlobalCoords lowerWeakOrigin_i = lowerWeakOrigin + origins[i];
            GlobalCoords upperWeakOrigin_i = upperWeakOrigin + origins[i];

            cuboidGeometries_[i] = std::make_shared<CuboidGeometry>(origins[i], length, height, depth);
            cuboidGeometries_[i]->addWeakeningPatch(subParset, upperWeakOrigin_i, weakBound(upperWeakOrigin_i));
            cuboidGeometries_[i]->addWeakeningPatch(subParset, lowerWeakOrigin_i, weakBound(lowerWeakOrigin_i));
        }

        const size_t idx = this->bodyCount_-1;
        origins[idx] = cuboidGeometries_[idx-1]->upperLeft();
        lowerWeakOrigin += origins[idx];

        cuboidGeometries_[idx] = std::make_shared<CuboidGeometry>(origins[idx], length, height, depth);
        cuboidGeometries_[idx]->addWeakeningPatch(subParset, lowerWeakOrigin, weakBound(lowerWeakOrigin));
#elif MY_DIM == 2
        auto weakBound = [&] (const GlobalCoords& origin) {
            GlobalCoords h = {origin[0] + weakLength, origin[1]};
            return h;
        };

        origins[0] = {0,0};
        GlobalCoords lowerWeakOrigin = {0.2, 0};
        GlobalCoords upperWeakOrigin = {0.2, height};

        cuboidGeometries_[0] = std::make_shared<CuboidGeometry>(origins[0], length, height);
        cuboidGeometries_[0]->addWeakeningPatch(subParset, upperWeakOrigin, weakBound(upperWeakOrigin));

        for (size_t i=1; i<this->bodyCount_-1; i++) {
            origins[i] = cuboidGeometries_[i-1]->upperLeft();
            auto height_i = height/3.0;

            GlobalCoords lowerWeakOrigin_i = lowerWeakOrigin + origins[i];
            GlobalCoords upperWeakOrigin_i = {upperWeakOrigin[0], height_i};
            upperWeakOrigin_i += origins[i];

            cuboidGeometries_[i] = std::make_shared<CuboidGeometry>(origins[i], length, height_i);


            cuboidGeometries_[i]->addWeakeningPatch(subParset, upperWeakOrigin_i, weakBound(upperWeakOrigin_i));
            cuboidGeometries_[i]->addWeakeningPatch(subParset, lowerWeakOrigin_i, weakBound(lowerWeakOrigin_i));
        }

        const size_t idx = this->bodyCount_-1;
        origins[idx] = cuboidGeometries_[idx-1]->upperLeft();
        lowerWeakOrigin += origins[idx];

        cuboidGeometries_[idx] = std::make_shared<CuboidGeometry>(origins[idx], length, height);
        cuboidGeometries_[idx]->addWeakeningPatch(subParset, lowerWeakOrigin, weakBound(lowerWeakOrigin));
#else
#error CuboidGeometry only supports 2D and 3D!"
#endif

    // set up reference grids
    gridConstructor_ = std::make_unique<GridsConstructor<HostGridType>>(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
        const auto& weakeningRegions = cuboidGeometry.weakeningRegions();
        for (size_t j=0; j<weakeningRegions.size(); j++) {
            refine(*grids[i], weakeningRegions[j], this->parset_.template get<double>("boundary.friction.smallestDiameter"), CuboidGeometry::lengthScale());
        }

        // determine minDiameter and maxDiameter
        double minDiameter = std::numeric_limits<double>::infinity();
        double maxDiameter = 0.0;
        for (auto &&e : elements(grids[i]->leafGridView())) {
          auto const geometry = e.geometry();
          auto const diam = diameter(geometry);
          minDiameter = std::min(minDiameter, diam);
          maxDiameter = std::max(maxDiameter, diam);
        }
        std::cout << "Grid" << i << " min diameter: " << minDiameter << std::endl;
        std::cout << "Grid" << i << " max diameter: " << maxDiameter << std::endl;
    }

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


template <class HostGridType, class VectorType>
void StackedBlocksFactory<HostGridType, VectorType>::setCouplings() {
    for (size_t i=0; i<this->bodyCount_; i++) {
        const auto& cuboidGeometry = *cuboidGeometries_[i];

        leafFaces_[i] = std::make_shared<LeafFaces>(this->bodies_[i]->gridView(), cuboidGeometry);
        levelFaces_[i] = std::make_shared<LevelFaces>(this->bodies_[i]->grid()->levelGridView(0), cuboidGeometry);
    }

    auto contactProjection = std::make_shared<Dune::Contact::NormalProjection<LeafBoundaryPatch>>();
    std::shared_ptr<typename Base::FrictionCouplingPair::BackEndType> backend = nullptr;


    for (size_t i=0; i<this->couplings_.size(); i++) {
      auto& coupling = this->couplings_[i];
      coupling = std::make_shared<typename Base::FrictionCouplingPair>();

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

      coupling->set(i, i+1, nonmortarPatch_[i], mortarPatch_[i], 0.1, Base::FrictionCouplingPair::CouplingType::STICK_SLIP, contactProjection, backend);
      coupling->setWeakeningPatch(weakPatches_[i]);
      coupling->setFrictionData(std::make_shared<MyGlobalFrictionData<GlobalCoords>>(this->parset_.sub("boundary.friction"), *weakPatches_[i], CuboidGeometry::lengthScale()));
    }
}

template <class HostGridType, class VectorType>
void StackedBlocksFactory<HostGridType, VectorType>::setLevelBodies() {
    const size_t nBodies = this->bodyCount_;

    size_t maxLevel = 0;
    std::vector<size_t> maxLevels(nBodies);
    for (size_t i=0; i<nBodies; i++) {
        const size_t bodyMaxLevel = this->bodies_[i]->grid()->maxLevel();
        maxLevels[i] = bodyMaxLevel;
        maxLevel = std::max(maxLevel, bodyMaxLevel);
    }

    for (int l=maxLevel; l>=0; l--) {
        this->contactNetwork_.addLevel(maxLevels, l);

        for (size_t i=0; i<nBodies; i++) {
            maxLevels[i]--;
        }
    }
}


template <class HostGridType, class VectorType>
void StackedBlocksFactory<HostGridType, VectorType>::setBoundaryConditions() {
    using LeafBoundaryCondition = BoundaryCondition<LeafGridView, 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->nVertices();

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

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

            std::shared_ptr<LeafBoundaryCondition> velocityDirichletBoundary = std::make_shared<LeafBoundaryCondition>("dirichlet");

            velocityDirichletBoundary->setBoundaryPatch(body->gridView(), velocityDirichletNodes);
            velocityDirichletBoundary->setBoundaryFunction(velocityDirichletFunction);
            body->addBoundaryCondition(velocityDirichletBoundary);
        } else {

        }
    }

    // lower Dirichlet Boundary
    const auto& firstBody = this->contactNetwork_.body(0);
    const auto& firstLeafVertexCount = firstBody->nVertices();
    std::shared_ptr<Dune::BitSetVector<dim>> zeroDirichletNodes = std::make_shared<Dune::BitSetVector<dim>>(firstLeafVertexCount);
    for (int j=0; j<firstLeafVertexCount; j++) {
        if (leafFaces_[0]->lower.containsVertex(j)) {
            for (size_t d=0; d<dim; d++) {
              (*zeroDirichletNodes)[j][d] = true;
            }
        }

        #if MY_DIM == 3 //TODO: wrong, needs revision
        if (leafFaces_[0]->front.containsVertex(j) || leafFaces_[0]->back.containsVertex(j))
            (*zeroDirichletNodes)[j][2] = true;
        #endif
    }

    std::shared_ptr<LeafBoundaryCondition> zeroDirichletBoundary = std::make_shared<LeafBoundaryCondition>("dirichlet");
    zeroDirichletBoundary->setBoundaryPatch(firstBody->gridView(), zeroDirichletNodes);

    std::shared_ptr<Function> zeroFunction = std::make_shared<NeumannCondition>();
    zeroDirichletBoundary->setBoundaryFunction(zeroFunction);
    firstBody->addBoundaryCondition(zeroDirichletBoundary);
}

#include "stackedblocksfactory_tmpl.cc"