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

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

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

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

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

#include "strikeslipfactory.hh"

template <class HostGridType, class VectorTEMPLATE>
void StrikeSlipFactory<HostGridType, VectorTEMPLATE>::setBodies() {
    // set up cuboid geometries

    std::array<double, 2> lengths = {this->parset_.template get<double>("body0.length"), this->parset_.template get<double>("body1.length")};
    std::array<double, 2> heights = {this->parset_.template get<double>("body0.height"), this->parset_.template get<double>("body1.height")};

    GlobalCoords origin(0);

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

#if MY_DIM == 3 // TODO: not implemented
        //std::array<double, 2> depths = {this->parset_.template get<double>("body0.depth"), this->parset_.template get<double>("body1.depth")};

        //triangleGeometries_[0] = std::make_shared<TriangleGeometry>(origins[0], lengths[0], heights[0], depths[0]);
        //triangleGeometries_[0]->addWeakeningPatch(frictionParset, {origins[0][0], origins[0][1]+ heights[0], 0}, {origins[0][0] + lengths[0], origins[0][1]+ heights[0], 0});

        //triangleGeometries_[1] = std::make_shared<TriangleGeometry>(origins[1], lengths[1], heights[1], depths[1]);
        //triangleGeometries_[1]->addWeakeningPatch(frictionParset, origins[1], {origins[1][0] + lengths[1], origins[1][1], 0});
#elif MY_DIM == 2
        triangleGeometries_[0] = std::make_shared<TriangleGeometry>(origin, lengths[0], heights[0]);
        triangleGeometries_[0]->addWeakeningPatch(frictionParset, triangleGeometries_[0]->A(), triangleGeometries_[0]->C());

        triangleGeometries_[1] = std::make_shared<TriangleGeometry>(triangleGeometries_[0]->C(), -lengths[1], -heights[1]);
        triangleGeometries_[1]->addWeakeningPatch(frictionParset, triangleGeometries_[1]->A(), triangleGeometries_[1]->C());
#else
#error CuboidGeometry only supports 2D and 3D!"
#endif

    // set up reference grids
    gridConstructor_ = std::make_unique<GridsConstructor<HostGridType>>(triangleGeometries_);
    auto& grids = gridConstructor_->getGrids();

    std::array<double, 2> smallestDiameter = {this->parset_.template get<double>("body0.smallestDiameter"), this->parset_.template get<double>("body1.smallestDiameter")};

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

        // define weak patch and refine grid
        const auto& weakeningRegions = triangleGeometry.weakeningRegions();
        for (size_t j=0; j<weakeningRegions.size(); j++) {
            refine(*grids[i], weakeningRegions[j], smallestDiameter[i], TriangleGeometry::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;
    }

#if MY_DIM == 2
    bodyData_[0] = std::make_shared<MyBodyData<dim>>(this->parset_.sub("body0"), 0.0, zenith_());
    this->bodies_[0] = std::make_shared<typename Base::LeafBody>(bodyData_[0], grids[0]);

    bodyData_[1] = std::make_shared<MyBodyData<dim>>(this->parset_.sub("body1"), 0.0, zenith_());
    this->bodies_[1] = std::make_shared<typename Base::LeafBody>(bodyData_[1], grids[1]);
#else
    bodyData_[0] = std::make_shared<MyBodyData<dim>>(this->parset_.sub("body0"), this->parset_.template get<double>("general.gravity"), zenith_());
    this->bodies_[0] = std::make_shared<typename Base::LeafBody>(bodyData_[0], grids[0]);

    bodyData_[1] = std::make_shared<MyBodyData<dim>>(this->parset_.sub("body1"), this->parset_.template get<double>("general.gravity"), zenith_());
    this->bodies_[1] = std::make_shared<typename Base::LeafBody>(bodyData_[1], grids[1]);
#endif
}

template <class HostGridType, class VectorTEMPLATE>
void StrikeSlipFactory<HostGridType, VectorTEMPLATE>::setLevelBodies() {
    const size_t maxLevel = std::max({this->bodies_[0]->grid()->maxLevel(), this->bodies_[1]->grid()->maxLevel()});

    for (size_t l=0; l<=maxLevel; l++) {
        std::vector<size_t> bodyLevels(2, l);
        this->contactNetwork_.addLevel(bodyLevels, l);
    }
}

template <class HostGridType, class VectorTEMPLATE>
void StrikeSlipFactory<HostGridType, VectorTEMPLATE>::setCouplings() {
    for (size_t i=0; i<this->bodyCount_; i++) {
        const auto& triangleGeometry = *triangleGeometries_[i];
        leafFaces_[i] = std::make_shared<LeafFaces>(this->bodies_[i]->gridView(), triangleGeometry);
        levelFaces_[i] = std::make_shared<LevelFaces>(this->bodies_[i]->grid()->levelGridView(0), triangleGeometry);
    }

    auto contactProjection = std::make_shared<Dune::Contact::NormalProjection<LeafBoundaryPatch>>();

    nonmortarPatch_[0] = std::make_shared<LevelBoundaryPatch>(levelFaces_[1]->b);
    mortarPatch_[0] = std::make_shared<LevelBoundaryPatch>(levelFaces_[0]->b);
    weakPatches_[0] = std::make_shared<WeakeningRegion>(triangleGeometries_[1]->weakeningRegions()[0]);

    auto& coupling = this->couplings_[0];
    coupling = std::make_shared<typename Base::FrictionCouplingPair>();

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

template <class HostGridType, class VectorTEMPLATE>
void StrikeSlipFactory<HostGridType, VectorTEMPLATE>::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>(-1.0*this->parset_.template get<double>("boundary.dirichlet.finalVelocity"), 0.25);

    const double lengthScale = TriangleGeometry::lengthScale();

    // body0
    const auto& body0 = this->contactNetwork_.body(0);
    const auto& leafVertexCount0 = body0->nVertices();
    std::shared_ptr<Dune::BitSetVector<dim>> zeroDirichletNodes = std::make_shared<Dune::BitSetVector<dim>>(leafVertexCount0);
    std::shared_ptr<Dune::BitSetVector<dim>> velocityDirichletNodes = std::make_shared<Dune::BitSetVector<dim>>(leafVertexCount0);

    const auto& gridView0 = body0->gridView();
    const auto& indexSet0 = gridView0.indexSet();

    for (const auto& e : elements(gridView0)) {
        for (const auto& is : intersections(gridView0, e)) {
            if (!is.boundary())
                continue;

            auto isCenter = is.geometry().center();

            // Dirichlet boundary (c)
            if (leafFaces_[0]->c.contains(is)) {
                isCenter -=  triangleGeometries_[0]->A();

                if (isCenter.two_norm()/std::abs(triangleGeometries_[0]->length()) > 0.1) {
                    const auto inside = is.inside();

                    auto refElement = Dune::ReferenceElements<double, dim>::general(inside.type());
                    int n = refElement.size(is.indexInInside(),1,dim);

                    for (int i=0; i<n; i++) {
                        int faceIdxi = refElement.subEntity(is.indexInInside(), 1, i, dim);
                        int globalIdx = indexSet0.template subIndex(inside,faceIdxi,dim);
                        (*zeroDirichletNodes)[globalIdx][1] = true;
                    }
                }
            }

            // velocity Dirichlet boundary (a)
            if (leafFaces_[0]->a.contains(is)) {
                isCenter -=  triangleGeometries_[0]->C();

                if (isCenter.two_norm()/std::abs(triangleGeometries_[0]->height()) > 0.1) {
                    const auto inside = is.inside();

                    auto refElement = Dune::ReferenceElements<double, dim>::general(inside.type());
                    int n = refElement.size(is.indexInInside(),1,dim);

                    for (int i=0; i<n; i++) {
                        int faceIdxi = refElement.subEntity(is.indexInInside(), 1, i, dim);
                        int globalIdx = indexSet0.template subIndex(inside,faceIdxi,dim);
                        (*velocityDirichletNodes)[globalIdx][0] = true;
                    }
                }
            }
        }
    }

    std::shared_ptr<Function> zeroFunction = std::make_shared<NeumannCondition>();

    std::shared_ptr<LeafBoundaryCondition> zeroDirichletBoundary = std::make_shared<LeafBoundaryCondition>("dirichlet");
    zeroDirichletBoundary->setBoundaryPatch(body0->gridView(), zeroDirichletNodes);
    zeroDirichletBoundary->setBoundaryFunction(zeroFunction);
    body0->addBoundaryCondition(zeroDirichletBoundary);

    std::shared_ptr<LeafBoundaryCondition> velocityDirichletBoundary = std::make_shared<LeafBoundaryCondition>("dirichlet");
    velocityDirichletBoundary->setBoundaryPatch(body0->gridView(), velocityDirichletNodes);
    velocityDirichletBoundary->setBoundaryFunction(velocityDirichletFunction);
    body0->addBoundaryCondition(velocityDirichletBoundary);

    // body1
    const auto& body1 = this->contactNetwork_.body(1);
    const auto& leafVertexCount1 = body1->nVertices();
    std::shared_ptr<Dune::BitSetVector<dim>> zeroDirichletNodes1 = std::make_shared<Dune::BitSetVector<dim>>(leafVertexCount1);
    std::shared_ptr<Dune::BitSetVector<dim>> loadNeumannNodes = std::make_shared<Dune::BitSetVector<dim>>(leafVertexCount1);

    const auto& gridView1 = body1->gridView();
    const auto& indexSet1 = gridView1.indexSet();

    for (const auto& e : elements(gridView1)) {
        for (const auto& is : intersections(gridView1, e)) {
            if (!is.boundary())
                continue;

            auto isCenter = is.geometry().center();

            // Neumann boundary, normal load (c)
            if (leafFaces_[1]->c.contains(is)) {
                isCenter -=  triangleGeometries_[1]->A();

                if (isCenter.two_norm()/std::abs(triangleGeometries_[1]->length()) > 0.1) {
                    const auto inside = is.inside();

                    auto refElement = Dune::ReferenceElements<double, dim>::general(inside.type());
                    int n = refElement.size(is.indexInInside(),1,dim);

                    for (int i=0; i<n; i++) {
                        int faceIdxi = refElement.subEntity(is.indexInInside(), 1, i, dim);
                        int globalIdx = indexSet1.template subIndex(inside,faceIdxi,dim);
                        (*loadNeumannNodes)[globalIdx][1] = true;
                    }
                }
            }

            // Dirichlet boundary (a)
            if (leafFaces_[1]->a.contains(is)) {
                isCenter -=  triangleGeometries_[1]->C();

                if (isCenter.two_norm()/std::abs(triangleGeometries_[1]->height()) > 0.1) {
                    const auto inside = is.inside();

                    auto refElement = Dune::ReferenceElements<double, dim>::general(inside.type());
                    int n = refElement.size(is.indexInInside(),1,dim);

                    for (int i=0; i<n; i++) {
                        int faceIdxi = refElement.subEntity(is.indexInInside(), 1, i, dim);
                        int globalIdx = indexSet1.template subIndex(inside,faceIdxi,dim);
                        (*zeroDirichletNodes1)[globalIdx][0] = true;
                    }
                }
            }
        }
    }

    std::shared_ptr<Function> constantFunction = std::make_shared<NeumannCondition>(this->parset_.template get<double>("boundary.neumann.sigmaN"));

    auto loadNeumannBoundary = std::make_shared<LeafBoundaryCondition>("neumann");
    loadNeumannBoundary->setBoundaryPatch(body1->gridView(), loadNeumannNodes);
    loadNeumannBoundary->setBoundaryFunction(constantFunction);
    body1->addBoundaryCondition(loadNeumannBoundary);

    std::shared_ptr<LeafBoundaryCondition> zeroDirichletBoundary1 = std::make_shared<LeafBoundaryCondition>("dirichlet");
    zeroDirichletBoundary1->setBoundaryPatch(body1->gridView(), zeroDirichletNodes1);
    zeroDirichletBoundary1->setBoundaryFunction(zeroFunction);
    body1->addBoundaryCondition(zeroDirichletBoundary1);

    // body0, body1: natural boundary conditions
    for (size_t i=0; i<this->bodyCount_; i++) {
        const auto& body = this->contactNetwork_.body(i);

        // Neumann boundary
        auto neumannBoundary = std::make_shared<LeafBoundaryCondition>("neumann");
        auto neumannPatch = std::make_shared<LeafBoundaryPatch>(body->gridView(), true);
        neumannBoundary->setBoundaryPatch(neumannPatch);
        neumannBoundary->setBoundaryFunction(neumannFunction);
        body->addBoundaryCondition(neumannBoundary);
    }
}

#include "strikeslipfactory_tmpl.cc"