supportpatchfactory.hh

#ifndef SUPPORT_PATCH_FACTORY_HH
#define SUPPORT_PATCH_FACTORY_HH

#include <set>
#include <queue>
#include <memory>

#include <dune/common/bitsetvector.hh>
#include <dune/common/fvector.hh>
#include <dune/common/hybridutilities.hh>

#include <dune/fufem/referenceelementhelper.hh>
#include <dune/grid/common/mcmgmapper.hh>

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

//#include "../../data-structures/levelcontactnetwork.hh"
#include "hierarchicleveliterator.hh"

#include "../../utils/debugutils.hh"

template <class LevelContactNetwork>
class NetworkIndexMapper {
private:
    static const int dim = LevelContactNetwork::dim; //TODO
    using ctype = typename LevelContactNetwork::GridView::ctype;

    using FiniteElementCache = Dune::PQkLocalFiniteElementCache<ctype, ctype, dim, 1>; // P1 finite element cache
    using Element = typename LevelContactNetwork::GridView::template Codim<0>::Entity;

    const FiniteElementCache cache_;

    const LevelContactNetwork& levelContactNetwork_;
    const size_t nBodies_;

    std::vector<size_t> vertexOffSet_;
    std::vector<size_t> faceOffSet_;
    std::vector<size_t> elementOffSet_;

public:
    NetworkIndexMapper(const LevelContactNetwork& levelContactNetwork) :
        levelContactNetwork_(levelContactNetwork),
        nBodies_(levelContactNetwork_.nBodies()),
        vertexOffSet_(nBodies_),
        faceOffSet_(nBodies_),
        elementOffSet_(nBodies_) {

        size_t vertexOffSet = 0;
        size_t faceOffSet = 0;
        size_t elementOffSet = 0;

        for (size_t i=0; i<nBodies_; i++) {
            const auto& gridView = levelContactNetwork_.body(i)->gridView();

            vertexOffSet_[i] = vertexOffSet;
            vertexOffSet = vertexOffSet_[i] + gridView.size(dim);

            faceOffSet_[i] = faceOffSet;
            faceOffSet = faceOffSet_[i] + gridView.size(1);

            elementOffSet_[i] = elementOffSet;
            elementOffSet = elementOffSet_[i] + gridView.size(0);
        }
    }

    size_t vertexIndex(const size_t bodyID, const Element& elem, const size_t localVertex) const {
        const auto& gridView = levelContactNetwork_.body(bodyID)->gridView();

        auto localIdx = cache().get(elem.type()).localCoefficients().localKey(localVertex).subEntity();
        return vertexOffSet_[bodyID] + gridView.indexSet().subIndex(elem, localIdx, dim);
    }

    template <class Intersection>
    size_t faceIndex(const size_t bodyID, const Intersection& is) const {
        const auto& gridView = levelContactNetwork_.body(bodyID)->gridView();
        return faceOffSet_[bodyID] + gridView.indexSet().subIndex(is.inside(), is.indexInInside(), 1);
    }

    size_t faceIndex(const size_t bodyID, const Element& elem, const size_t indexInInside) const {
        const auto& gridView = levelContactNetwork_.body(bodyID)->gridView();
        return faceOffSet_[bodyID] + gridView.indexSet().subIndex(elem, indexInInside, 1);
    }

    size_t elementIndex(const size_t bodyID, const Element& elem) const {
        const auto& gridView = levelContactNetwork_.body(bodyID)->gridView();
        return elementOffSet_[bodyID] + gridView.indexSet().index(elem);
    }

    const auto& cache() const {
        return cache_;
    }
};



/**
 *  constructs BitSetVector whose entries are
 *      true,   if vertex is outside of patch or part of patch boundary
 *      false,  if vertex is in interior of patch
 */
template <class LevelContactNetwork>
class SupportPatchFactory
{
public:
    using Patch = Dune::BitSetVector<1>;

private:
    static const int dim = LevelContactNetwork::dim; //TODO
    using ctype = typename LevelContactNetwork::ctype;

    using Element = typename LevelContactNetwork::GridView::template Codim<0>::Entity;

    struct RemoteIntersection {
        size_t bodyID; // used to store bodyID of outside elem
        size_t contactCouplingID;
        size_t intersectionID;

        bool flip = false;

        void set(const size_t _bodyID, const size_t _contactCouplingID, const size_t _intersectionID, const bool _flip = false) {
            bodyID = _bodyID;
            contactCouplingID = _contactCouplingID;
            intersectionID = _intersectionID;
            flip = _flip;
        }

        auto get(const LevelContactNetwork& levelContactNetwork) const {
            return std::move(levelContactNetwork.glue(contactCouplingID)->getIntersection(intersectionID));
        }
    };

    struct BodyElement {
        size_t bodyID;
        Element element;

        void set(const size_t _bodyID, const Element& _elem) {
            bodyID = _bodyID;
            element = _elem;
        }

        void set(const LevelContactNetwork& levelContactNetwork, const RemoteIntersection& rIs) {
            const auto& is = rIs.get(levelContactNetwork);

            if (rIs.flip) {
                set(rIs.bodyID, is.inside());
            } else {
                set(rIs.bodyID, is.outside());
            }
        }
    };

    using NetworkIndexMapper = NetworkIndexMapper<LevelContactNetwork>;

    const size_t nBodies_;

    const LevelContactNetwork& coarseContactNetwork_;
    const LevelContactNetwork& fineContactNetwork_;

    size_t nVertices_ = 0;

    std::map<size_t, std::vector<RemoteIntersection>> coarseIntersectionMapper_; // map faceID to neighbor elements
    std::map<size_t, std::vector<RemoteIntersection>> fineIntersectionMapper_; // map faceID to neighbor elements
    std::map<size_t, std::vector<size_t>> neighborFaceDofs_; // map faceID to neighbor element dofs

    std::set<size_t> interiorContactDofs_;

    const NetworkIndexMapper coarseIndices_;
    const NetworkIndexMapper fineIndices_;

public:
    SupportPatchFactory(const LevelContactNetwork& coarseContactNetwork, const LevelContactNetwork& fineContactNetwork) :
        nBodies_(coarseContactNetwork.nBodies()),
        coarseContactNetwork_(coarseContactNetwork),
        fineContactNetwork_(fineContactNetwork),
        coarseIndices_(coarseContactNetwork_),
        fineIndices_(fineContactNetwork_) {

        assert(nBodies_ == fineContactNetwork_.nBodies());

        for (size_t i=0; i<nBodies_; i++) {
            nVertices_ += fineContactNetwork_.body(i)->nVertices();
        }

        setup(coarseContactNetwork_, coarseIndices_, coarseIntersectionMapper_);
        setup(fineContactNetwork_, fineIndices_, fineIntersectionMapper_);

        // setup neighborFaceDofs_
        for (size_t i=0; i<coarseContactNetwork_.nCouplings(); i++) {
                const auto& coupling = *coarseContactNetwork_.coupling(i);
                const auto& glue = *coarseContactNetwork_.glue(i);

                const size_t nmBody = coupling.gridIdx_[0];
                const size_t mBody = coupling.gridIdx_[1];

                for (const auto& rIs : intersections(glue)) {
                    size_t nmFaceIdx = coarseIndices_.faceIndex(nmBody, rIs.inside(), rIs.indexInInside());
                    size_t mFaceIdx = coarseIndices_.faceIndex(mBody, rIs.outside(), rIs.indexInOutside());

                    std::vector<std::set<size_t>> dofs;
                    remoteIntersectionDofs(coarseIndices_, rIs, nmBody, mBody, dofs);

                    if (dofs[1].size()>0) {
                        neighborFaceDofs_[nmFaceIdx].insert(neighborFaceDofs_[nmFaceIdx].end(), dofs[1].begin(), dofs[1].end());
                    }

                    if (dofs[0].size()>0) {
                        neighborFaceDofs_[mFaceIdx].insert(neighborFaceDofs_[mFaceIdx].end(), dofs[0].begin(), dofs[0].end());
                    }
                }
        }

        for (size_t i=0; i<nBodies_; i++) {
            std::cout << "Coarse body" << i << ":" << std::endl;

            for (const auto& e : elements(coarseContactNetwork_.body(i)->gridView())) {
                std::cout << "elemID: " << coarseIndices_.elementIndex(i, e) << std::endl;
                std::cout << "vertexIDs: ";
                const int dimElement = Element::dimension;
                const auto& refElement = Dune::ReferenceElements<double, dimElement>::general(e.type());

                for (int j=0; j<refElement.size(dim); j++) {
                    std::cout << coarseIndices_.vertexIndex(i, e, j) << " ";
                }
                std::cout << std::endl;

                std::cout << "intersectionIDs: ";
                for (const auto& is : intersections(coarseContactNetwork_.body(i)->gridView(), e)) {
                    std::cout << coarseIndices_.faceIndex(i, e, is.indexInInside()) << " (";

                    std::set<size_t> dofs;
                    intersectionDofs(coarseIndices_, is, i, dofs);
                    for (auto& d : dofs) {
                        std::cout << d << " ";
                    }
                    std::cout << "); ";

                }
                std::cout << std::endl  << "--------------" << std::endl;
            }
        }

        std::cout << std::endl;
        for (auto& kv : neighborFaceDofs_) {
            std::cout << "faceID: " << kv.first << " dofs: ";

            const auto& dofs = kv.second;
            for (auto& d : dofs) {
                std::cout << d << " ";
            }
            std::cout << std::endl;
        }
        /*
            // neighborElementDofs_ contains dofs of connected intersections that are neighbors to a face given by faceID,
            // boundary dofs are contained once, interior dofs multiple times
            // task: remove boundary dofs and redundant multiples of interior dofs
            // TODO: only works in 2D
            for (auto& it : neighborElementDofs) {
                auto& dofs = it.second;

                std::sort(dofs.begin(), dofs.end());
                auto& dof = dofs.begin();
                while (dof != dofs.end()) {
                    auto next = dof;
                    if (*dof != *(++next))
                        dof = dofs.erase(dof);
                }
                dofs.erase(std::unique(dofs.begin(), dofs.end()), dofs.end());
            }
        */
    }

    void build(const size_t bodyID, const Element& coarseElement, const size_t localVertex, Patch& patchDofs, const size_t patchDepth = 0) {
        std::cout << "Building SupportPatch... ";

        std::cout << "elemID: " << coarseIndices_.elementIndex(bodyID, coarseElement) << " vertexID: " << coarseIndices_.vertexIndex(bodyID, coarseElement, localVertex) << std::endl;

        patchDofs.resize(nVertices_);
        patchDofs.setAll();

        BodyElement seedElement;
        seedElement.set(bodyID, coarseElement);

        auto isPatchIntersection = [this](auto& is, const size_t bodyID, const std::set<size_t>& patchVertices, std::set<size_t>& newPatchVertices) {
            std::set<size_t> dofs;
            intersectionDofs(coarseIndices_, is, bodyID, dofs);

            for (auto& dof : dofs) {
                if (patchVertices.count(dof)) {
                    newPatchVertices.insert(dofs.begin(), dofs.end());
                    return true;
                }
            }
            return false;
        };

        auto isRPatchIntersection = [this](std::vector<std::set<size_t>>& dofs, const std::set<size_t>& patchVertices) {
            for (auto& dof : dofs[0]) {
                if (patchVertices.count(dof)) {
                    return true;
                }
            }
            return false;
        };

        // construct coarse patch
        std::vector<BodyElement> patchElements;
        patchElements.emplace_back(seedElement);

        std::set<size_t> visitedElements;
        //visitedElements.insert(coarseIndices_.elementIndex(seedElement.bodyID, seedElement.element));

        std::set<size_t> patchVertices;
        patchVertices.insert(coarseIndices_.vertexIndex(bodyID, coarseElement, localVertex));

        std::queue<BodyElement> patchSeeds;
        patchSeeds.push(seedElement);

        std::set<size_t> visitedBodies;

        for (size_t depth=0; depth<=patchDepth; depth++) {
            std::vector<BodyElement> nextSeeds;
            std::set<size_t> newPatchVertices;

            while (!patchSeeds.empty()) {
                const auto& patchSeed = patchSeeds.front();

                size_t elemIdx = coarseIndices_.elementIndex(patchSeed.bodyID, patchSeed.element);

                if (visitedElements.count(elemIdx)) {
                    patchSeeds.pop();
                    continue;
                }

                visitedBodies.insert(patchSeed.bodyID);

                visitedElements.insert(elemIdx);

                size_t bodyIdx = patchSeed.bodyID;
                const auto& elem = patchSeed.element;
                const auto& gridView = coarseContactNetwork_.body(bodyIdx)->gridView();

                for (const auto& it : intersections(gridView, elem)) {
                    if (isPatchIntersection(it, bodyIdx, patchVertices, newPatchVertices)) {
                        if (it.neighbor()) {
                            BodyElement neighbor;
                            neighbor.set(bodyIdx, it.outside());

                            size_t neighborIdx = coarseIndices_.elementIndex(neighbor.bodyID, neighbor.element);

                            std::cout << "elementID: " << neighborIdx << std::endl;

                            if (visitedElements.count(neighborIdx))
                                continue;

                            patchElements.emplace_back(neighbor);
                            patchSeeds.push(neighbor);
                        } else {
                            size_t faceIdx = coarseIndices_.faceIndex(bodyIdx, it);

                            if (coarseIntersectionMapper_.count(faceIdx)) {
                                const auto& intersections = coarseIntersectionMapper_[faceIdx];
                                for (size_t i=0; i<intersections.size(); i++) {
                                    const auto& intersection = intersections[i];
                                    const auto& rIs = intersection.get(coarseContactNetwork_);

                                    BodyElement neighbor;
                                    neighbor.set(coarseContactNetwork_, intersection);

                                    size_t neighborIdx = coarseIndices_.elementIndex(neighbor.bodyID, neighbor.element);

                                    if (visitedElements.count(neighborIdx) || visitedBodies.count(neighbor.bodyID))
                                        continue;

                                    std::vector<std::set<size_t>> dofs;
                                    remoteIntersectionDofs(coarseIndices_, rIs, bodyIdx, neighbor.bodyID, dofs, intersection.flip);

                                    if (isRPatchIntersection(dofs, patchVertices)) {
                                        patchVertices.insert(dofs[1].begin(), dofs[1].end());
                                        patchElements.emplace_back(neighbor);
                                        patchSeeds.push(neighbor);
                                    } else {
                                        newPatchVertices.insert(dofs[1].begin(), dofs[1].end());
                                        nextSeeds.emplace_back(neighbor);
                                    }

                                    print(patchVertices, "patchVertices: ");
                                }
                            }
                        }
                    } else {
                        if (it.neighbor()) {
                            BodyElement neighbor;
                            neighbor.set(bodyIdx, it.outside());

                            size_t neighborIdx = coarseIndices_.elementIndex(neighbor.bodyID, neighbor.element);

                            if (visitedElements.count(neighborIdx))
                                continue;

                            nextSeeds.emplace_back(neighbor);
                        } else {
                            size_t faceIdx = coarseIndices_.faceIndex(bodyIdx, it);

                            if (coarseIntersectionMapper_.count(faceIdx)) {
                                const auto& intersections = coarseIntersectionMapper_[faceIdx];
                                for (size_t i=0; i<intersections.size(); i++) {
                                    BodyElement neighbor;
                                    neighbor.set(coarseContactNetwork_, intersections[i]);

                                    size_t neighborIdx = coarseIndices_.elementIndex(neighbor.bodyID, neighbor.element);

                                    if (visitedElements.count(neighborIdx) || visitedBodies.count(neighbor.bodyID))
                                        continue;

                                    newPatchVertices.insert(neighborFaceDofs_[faceIdx].begin(), neighborFaceDofs_[faceIdx].end());
                                    nextSeeds.emplace_back(neighbor);
                                }
                            }
                        }
                    }
                }

                patchSeeds.pop();
            }

            for (size_t i=0; i<nextSeeds.size(); i++) {
                patchSeeds.push(nextSeeds[i]);
            }
            patchVertices.insert(newPatchVertices.begin(), newPatchVertices.end());
            print(patchVertices, "patchVertices: ");
        }

        std::cout << "constructing fine patch... " << std::endl;

        // construct fine patch
        using HierarchicLevelIteratorType = HierarchicLevelIterator<typename LevelContactNetwork::GridType>;
        std::set<size_t> visited;
        std::set<size_t> patchBoundary;
        for (size_t i=0; i<patchElements.size(); ++i) {
            const auto& coarseElem = patchElements[i];

            size_t elemIdx = coarseIndices_.elementIndex(coarseElem.bodyID, coarseElem.element);
            if (visited.count(elemIdx))
                continue;

            std::cout << "coarse element: " << elemIdx << std::endl;

            visited.insert(elemIdx);

            const auto& grid = coarseContactNetwork_.body(coarseElem.bodyID)->gridView().grid();
            const auto fineLevel = fineContactNetwork_.body(coarseElem.bodyID)->level();

            // add fine patch elements
            HierarchicLevelIteratorType endHierIt(grid, coarseElem.element, HierarchicLevelIteratorType::PositionFlag::end, fineLevel);
            HierarchicLevelIteratorType hierIt(grid, coarseElem.element, HierarchicLevelIteratorType::PositionFlag::begin, fineLevel);
            for (; hierIt!=endHierIt; ++hierIt) {
                const Element& fineElem = *hierIt;
                addLocalDofs(coarseElem, fineElem, visitedElements, patchBoundary, patchDofs);
            }
        }

        std::cout << "Success!" << std::endl;
    }

    const auto& coarseContactNetwork() {
        return coarseContactNetwork_;
    }

private:
    void setup(const LevelContactNetwork& levelContactNetwork, const NetworkIndexMapper& indices, std::map<size_t, std::vector<RemoteIntersection>>& faceToRIntersections) {
        for (size_t i=0; i<levelContactNetwork.nCouplings(); i++) {
            const auto& coupling = *levelContactNetwork.coupling(i);
            const auto& glue = *levelContactNetwork.glue(i);

            const size_t nmBody = coupling.gridIdx_[0];
            const size_t mBody = coupling.gridIdx_[1];

            size_t rIsID = 0;
            for (const auto& rIs : intersections(glue)) {
                size_t nmFaceIdx = indices.faceIndex(nmBody, rIs.inside(), rIs.indexInInside());

                RemoteIntersection nmIntersection;
                nmIntersection.set(mBody, i, rIsID, false);
                faceToRIntersections[nmFaceIdx].emplace_back(nmIntersection);

                if (rIs.neighbor()) {
                    size_t mFaceIdx = indices.faceIndex(mBody, rIs.outside(), rIs.indexInOutside());

                    RemoteIntersection mIntersection;
                    mIntersection.set(nmBody, i, rIsID, true);
                    faceToRIntersections[mFaceIdx].emplace_back(mIntersection);
                }

                rIsID++;
            }
        }
    }

    template <class Intersection>
    bool containsInsideSubentity(const Element& elem, const Intersection& intersection, int subEntity, int codim) {
	    return ReferenceElementHelper<double, dim>::subEntityContainsSubEntity(elem.type(), intersection.indexInInside(), 1, subEntity, codim);
	}
      
    void addLocalDofs(const BodyElement& coarseElem, const Element& fineElem, const std::set<size_t>& coarsePatch, std::set<size_t>& patchBoundary, Patch& patchDofs) {
	    // insert local dofs
        const auto& gridView = fineContactNetwork_.body(coarseElem.bodyID)->gridView();
        const size_t bodyID = coarseElem.bodyID;

        const auto coarseLevel = coarseContactNetwork_.body(coarseElem.bodyID)->level();

	    // search for parts of the global and inner boundary 
        for (const auto& is : intersections(gridView, fineElem)) {
            if (is.neighbor()) {
                std::set<size_t> isDofs;
                intersectionDofs(fineIndices_, is, bodyID, isDofs);

                auto outsideFather = coarseFather(is.outside(), coarseLevel);
		    
                if (!coarsePatch.count(coarseIndices_.elementIndex(coarseElem.bodyID, outsideFather))) {
                    addToPatchBoundary(isDofs, patchBoundary, patchDofs);
                } else {
                    addToPatch(isDofs, patchBoundary, patchDofs);
                }
            } else {
                size_t faceIdx = fineIndices_.faceIndex(bodyID, is);
                if (fineIntersectionMapper_.count(faceIdx)) {
                    const auto& intersections = fineIntersectionMapper_[faceIdx];

                    for (size_t i=0; i<intersections.size(); i++) {
                        const auto& intersection = intersections[i];
                        const auto& rIs = intersection.get(fineContactNetwork_);
                        const size_t outsideID = intersection.bodyID;

                        std::vector<std::set<size_t>> rIsDofs;
                        remoteIntersectionDofs(fineIndices_, rIs, bodyID, outsideID, rIsDofs, intersection.flip);

                        if (rIs.neighbor()) {
                            Element outsideFather;
                            const size_t outsideLevel = coarseContactNetwork_.body(outsideID)->level();

                            if (intersection.flip)
                                outsideFather = coarseFather(rIs.inside(), outsideLevel);
                            else
                                outsideFather = coarseFather(rIs.outside(), outsideLevel);

                            std::set<size_t> totalRIsDofs(rIsDofs[0]);
                            totalRIsDofs.insert(rIsDofs[1].begin(), rIsDofs[1].end());

                            if (!coarsePatch.count(coarseIndices_.elementIndex(outsideID, outsideFather))) {
                                addToPatchBoundary(totalRIsDofs, patchBoundary, patchDofs);
                            } else {
                                addToPatch(totalRIsDofs, patchBoundary, patchDofs);
                            }
                        } else {
                            addToPatchBoundary(rIsDofs[0], patchBoundary, patchDofs);
                        }
                    }
                } else {
                    std::set<size_t> isDofs;
                    intersectionDofs(fineIndices_, is, bodyID, isDofs);

                    addToPatchBoundary(isDofs, patchBoundary, patchDofs);
                }
            }
        }
	}

    auto coarseFather(const Element& fineElem, const size_t coarseLevel) const {
        Element coarseElem = fineElem;
        while (coarseElem.level() > coarseLevel) {
            coarseElem = coarseElem.father();
        }
        return coarseElem;
    }

    void addToPatchBoundary(const std::set<size_t>& dofs, std::set<size_t>& patchBoundary, Patch& patchDofs) {
        for (auto& dof : dofs) {
            patchBoundary.insert(dof);
            patchDofs[dof][0] = true;
        }
    }

    void addToPatch(const std::set<size_t>& dofs, const std::set<size_t>& patchBoundary, Patch& patchDofs) {
        for (auto& dof : dofs) {
            if (!patchBoundary.count(dof)) {
                patchDofs[dof][0] = false;
            }
        }
    }

    template <class RIntersection>
    void remoteIntersectionDofs(const NetworkIndexMapper& indices, const RIntersection& is, const size_t insideID, const size_t outsideID, std::vector<std::set<size_t>>& dofs, bool flip = false) {
        assert(!flip || (flip && is.neighbor()));

        dofs.resize(2);
        dofs[0].clear();
        dofs[1].clear();

        const auto& isGeo = is.geometry();
        const auto& isRefElement = Dune::ReferenceElements<ctype, dim-1>::general(is.type());

        if (!flip) {
             const auto& inside = is.inside();
             const auto& insideGeo = inside.geometry();
             const auto& insideRefElement = Dune::ReferenceElements<ctype, dim>::general(inside.type());
             const int dimElement = Dune::ReferenceElements<ctype, dim>::dimension;

             for (int i=0; i<insideRefElement.size(is.indexInInside(), 1, dimElement); i++) {
                 size_t idxInElement = insideRefElement.subEntity(is.indexInInside(), 1, i, dimElement);

                 const auto& localCorner = isGeo.local(insideGeo.corner(idxInElement));

                 if (isRefElement.checkInside(localCorner)) {
                     dofs[0].insert(indices.vertexIndex(insideID, inside, idxInElement));
                 }
             }

             if (is.neighbor()) {
                 const auto& outside = is.outside();
                 const auto& outsideGeo = outside.geometry();
                 const auto& outsideRefElement = Dune::ReferenceElements<ctype, dim>::general(outside.type());

                 for (int i=0; i<outsideRefElement.size(is.indexInOutside(), 1, dimElement); i++) {
                     size_t idxInElement = outsideRefElement.subEntity(is.indexInOutside(), 1, i, dimElement);

                     const auto& localCorner = isGeo.local(outsideGeo.corner(idxInElement));

                     if (isRefElement.checkInside(localCorner)) {
                         dofs[1].insert(indices.vertexIndex(outsideID, outside, idxInElement));
                     }
                 }
             }
        } else {
            const auto& inside = is.outside();
            const auto& insideGeo = inside.geometry();
            const auto& insideRefElement = Dune::ReferenceElements<ctype, dim>::general(inside.type());
            const int dimElement = Dune::ReferenceElements<ctype, dim>::dimension;

            for (int i=0; i<insideRefElement.size(is.indexInOutside(), 1, dimElement); i++) {
                size_t idxInElement = insideRefElement.subEntity(is.indexInOutside(), 1, i, dimElement);

                const auto& localCorner = isGeo.local(insideGeo.corner(idxInElement));

                if (isRefElement.checkInside(localCorner)) {
                    dofs[0].insert(indices.vertexIndex(insideID, inside, idxInElement));
                }
            }

            const auto& outside = is.inside();
            const auto& outsideGeo = outside.geometry();
            const auto& outsideRefElement = Dune::ReferenceElements<ctype, dim>::general(outside.type());

            for (int i=0; i<outsideRefElement.size(is.indexInInside(), 1, dimElement); i++) {
                size_t idxInElement = outsideRefElement.subEntity(is.indexInInside(), 1, i, dimElement);

                const auto& localCorner = isGeo.local(outsideGeo.corner(idxInElement));

                if (isRefElement.checkInside(localCorner)) {
                    dofs[1].insert(indices.vertexIndex(outsideID, outside, idxInElement));
                }
            }
        }

    }

    template <class Intersection>
    void intersectionDofs(const NetworkIndexMapper& indices, const Intersection& is, const size_t bodyID, std::set<size_t>& dofs) {
        dofs.clear();

        const Element& elem = is.inside();
        size_t intersectionIndex = is.indexInInside();

        const int dimElement = Element::dimension;
        const auto& refElement = Dune::ReferenceElements<double, dimElement>::general(elem.type());

        for (int i=0; i<refElement.size(intersectionIndex, 1, dimElement); i++) {
            size_t idxInElement = refElement.subEntity(intersectionIndex, 1, i, dimElement);
            dofs.insert(indices.vertexIndex(bodyID, elem, idxInElement));
        }
    }
};
#endif