facehierarchy.hh

Go to the documentation of this file.

// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set ts=8 sw=4 et sts=4:
 
// SPDX-FileCopyrightText: Copyright © DUNE-FUFEM Project contributors, see file AUTHORS.md
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception OR LGPL-3.0-or-later
 
#ifndef DUNE_FUFEM_FACE_HIERARCHY_HH
#define DUNE_FUFEM_FACE_HIERARCHY_HH
 
#include <dune/geometry/referenceelements.hh>
 
 
 
// forward declarations
template<class GridType> class Face;
template<class GridType> struct FaceHierarchy;
template<class GridType> class HierarchicFaceIterator;
 
 
 
template<class GridType>
struct FaceHierarchy
{
 
    public:
        typedef typename GridType::template Codim<0>::Entity Element;
        typedef typename GridType::LevelGridView::Intersection LevelIntersection;
        typedef typename GridType::LevelGridView::IntersectionIterator LevelIntersectionIterator;
 
        enum {dim = Element::dimension};
 
        static bool isSubFace(const Element& father, const int fatherFace, const Element& son, const int sonFace)
        {
            typedef typename Element::LocalGeometry::GlobalCoordinate LocalFatherCoordinate;
 
            auto fatherRefElement = Dune::ReferenceElements<double, dim>::general(father.type());
            auto sonRefElement = Dune::ReferenceElements<double, dim>::general(son.type());
 
            LocalFatherCoordinate centerInFather = son.geometryInFather().global(sonRefElement.position(sonFace, 1));
 
            const auto& mapping = fatherRefElement.template geometry<1>(fatherFace);
 
            centerInFather -= mapping.global(mapping.local(centerInFather));
 
            return centerInFather.two_norm2() < 1e-4;
        }
 
        static int findSubFace(const Element& father, const int fatherFace, const Element& child)
        {
 
            for(unsigned int i=0; i < child.subEntities(1); ++i)
                if (isSubFace(father, fatherFace, child, i))
                    return i;
            return -1;
        }
 
        static int findFatherFace(const GridType& grid, const Element& element, const int face)
        {
            const Element& father = element.father();
            for(unsigned int i=0; i < father->subEntities(1); ++i)
                if (isSubFace(father, i, element, face))
                    return i;
            return -1;
        }
 
        static LevelIntersectionIterator findFatherFaceIntersection(const GridType& grid, const Element& element, const int face)
        {
            const Element& father = element.father();
 
            LevelIntersectionIterator nIt = grid.levelGridView(father.level()).ibegin(father);
            LevelIntersectionIterator nEnd = grid.levelGridView(father.level()).iend(father);
            for(; nIt != nEnd; ++nIt)
            {
                if (isSubFace(father, nIt->indexInInside(), element, face))
                    return nIt;
            }
            DUNE_THROW(Dune::Exception, "No father face found!");
        }
 
        static Element finestCoarseNeighbor(const GridType& grid, const Element& element, int face)
        {
            LevelIntersectionIterator nIt = grid.levelGridView(element.level()).ibegin(element);
            LevelIntersectionIterator nEnd = grid.levelGridView(element.level()).iend(element);
            for(; nIt != nEnd; ++nIt)
            {
                if ((nIt->neighbor()) and (nIt->indexInInside()==face))
                    return nIt->outside();
            }
 
            Element checkElement(element);
            while(checkElement.level()>0)
            {
                LevelIntersectionIterator fatherIntersection = findFatherFaceIntersection(grid, checkElement, face);
                if (fatherIntersection->neighbor())
                    return fatherIntersection->outside();
                checkElement = checkElement.father();
                face = fatherIntersection->indexInInside();
            }
            return element;
        }
 
};
 
 
 
template<class GridType>
class Face
{
        typedef typename GridType::LevelGridView::IntersectionIterator LevelIntersectionIterator;
        enum {dim = GridType::dimension};
 
    public:
        typedef typename GridType::template Codim<0>::Entity Element;
        typedef typename GridType::LevelGridView::Intersection Intersection;
        typedef typename GridType::LevelGridView::IntersectionIterator IntersectionIterator;
        typedef HierarchicFaceIterator<GridType> HierarchicIterator;
        typedef Face<GridType> FaceType;
 
        Face(const GridType& grid, const Element &e, int index=-1):
            grid_(&grid),
            element_(e),
            index_(index)
        {}
 
        template<class OtherFace>
        Face(const OtherFace& other):
            grid_(&other.grid()),
            element_(other.element()),
            index_(other.index())
        {}
 
        const Element& element() const
        {
            return element_;
        }
 
        void setElement(const Element& e)
        {
            element_ = e;
        }
 
        int index() const
        {
            return index_;
        }
 
        void setIndex(int i)
        {
            index_ = i;
        }
 
        int level() const
        {
            return element().level();
        }
 
        const GridType& grid() const
        {
            return *grid_;
        }
 
        bool isInitialized() const
        {
            return index_>=0;
        }
 
        HierarchicIterator hbegin(int maxLevel) const
        {
            return HierarchicIterator(*grid_, *this, maxLevel, HierarchicIterator::begin);
        }
 
        HierarchicIterator hend(int maxLevel) const
        {
            return HierarchicIterator(*grid_, *this, maxLevel, HierarchicIterator::end);
        }
 
        FaceType father() const
        {
            FaceType fatherFace(*grid_, element().father());
            for(int i=0; i < fatherFace.element().subEntities(1); ++i)
            {
                if (isSubFace(fatherFace.element(), i, element(), index_))
                {
                    fatherFace.setIndex(i);
                    return fatherFace;
                }
            }
            return fatherFace;
        }
 
        template<class OtherFace>
        bool isFather(const OtherFace& fatherFace) const
        {
            if (not(isInitialized()) or not(fatherFace.isInitialized()))
                return false;
            if (fatherFace.element() != *element().father())
                return false;
            return FaceHierarchy<GridType>::isSubFace(fatherFace.element(), fatherFace.index(), *element_, index_);
        }
 
        template<class OtherFace>
        bool isChild(const OtherFace& childFace) const
        {
            return FaceHierarchy<GridType>::isSubFace(element(), index_, childFace.element(), childFace.index());
        }
 
        IntersectionIterator intersection() const
        {
            LevelIntersectionIterator nIt = grid_->levelGridView(level()).ibegin(element());
            LevelIntersectionIterator nEnd = grid_->levelGridView(level()).iend(element());
            for(; nIt != nEnd; ++nIt)
            {
                if (nIt->indexInInside()==index_)
                    return nIt;
            }
        }
 
        IntersectionIterator fatherIntersection() const
        {
            return FaceHierarchy<GridType>::findFatherFaceIntersection(*grid_, element(), index_);
        }
 
#if 0
        Face<GridType> outside() const
        {
            Intersection i = intersection();
            if (i.neighbor())
                return Face<GridType>(*grid_, i.outside(), i.indexInOutside());
            return Face<GridType>(*grid_);
        }
 
        Face<GridType> coarseNeighbor() const
        {
            Face<GridType> result = outside();
 
            if (result.isInitialized())
                return result;
 
            auto checkElement = element_;
            int checkFace = index_;
            while(checkElement.level()>0)
            {
                Intersection fatherInt = fatherIntersection(*grid_, checkElement, checkFace);
                if (fatherInt.neighbor())
                    return Face<GridType>(*grid_, fatherInt.outside(), fatherInt.indexInOutside());
                checkElement = checkElement->father();
                checkFace = fatherInt.indexInInside();
            }
            return Face<GridType>(*grid_);
        }
#endif
 
    protected:
        const GridType* grid_;
        Element element_;
        int index_;
};
 
template<class GridType>
class HierarchicFaceIterator :
    public Dune::ForwardIteratorFacade< HierarchicFaceIterator<GridType>, const Face<GridType> >
{
        enum {dim=GridType::dimension};
 
        typedef typename GridType::HierarchicIterator HierarchicElementIterator;
        typedef typename GridType::template Codim<0>::Entity Element;
 
 
    public:
 
        enum PositionFlag {begin, end};
 
        using FaceType = Face<GridType>;
 
        HierarchicFaceIterator(const GridType& grid, const FaceType& fatherFace, const int maxLevel, PositionFlag flag) :
            grid_(grid),
            maxLevel_(maxLevel),
            fatherFace_(fatherFace),
            face_(fatherFace)
        {
            if (flag==begin)
            {
                stack_.reserve(grid.maxLevel()-fatherFace_.element().level());
                incrementFrom(fatherFace_.element());
            }
        }
 
        void increment()
        {
            incrementFrom(*stack_.back().it_);
        }
 
        bool equals(const HierarchicFaceIterator& other) const
        {
            return ((stack_.size() == 0) and (other.stack_.size() == 0));
        }
 
        const FaceType& dereference() const
        {
            return face_;
        }
 
        const Element element() const
        {
            return *stack_.back().it_;
        }
 
        int index() const
        {
            return stack_.back().index_;
        }
 
 
    protected:
 
        // The following methods implement a depth first of
        // sub faces in all descendant elements as depicted
        // below:
        //
        // *-3-*-4-*-6-*-7-*
        // *---2---*---5---*
        // *-------1-------*
        //
        // All ancestors are stored in the stack in order
        // to visit their siblings later on.
        //
        // To this end we need to store the triple ChildFaceIterator
        // containing the element iterator, the face index, and the
        // end element iterator.
 
        // increment starting from given element
        // this is needed since the constructor increments
        // from the entity pointer in father face and increment()
        // increments from the current back of the stack
        void incrementFrom(const Element& e)
        {
            // If there are children to consider go up in the tree
            // and put child iterator on the stack.
            //
            // Notice that we only consider children of elements
            // that contain a subface since the current back of
            // the stack pointed to such an element.
            int level = e.level();
            if ((level < maxLevel_) and not(e.isLeaf()))
            {
                stack_.push_back(ChildFaceIterator(e.hbegin(level+1), e.hend(level+1)));
                if (findChildInBack())
                    return;
            }
 
            if (stack_.size() == 0)
                return;
 
            // Iterate on current level and go down in the tree
            // until find a child face or reach the end (=stack is empty)
            do
            {
                ++(stack_.back().it_);
                while (stack_.back().it_ == stack_.back().end_)
                {
                    stack_.pop_back();
                    if (stack_.size() == 0)
                        return;
                    ++(stack_.back().it_);
                }
            } while(not(findChildInBack()));
        }
 
        // find child face in back of stack
        // and set the corresponding face index
        // return whether a child face was found
        bool findChildInBack()
        {
            bool found = false;
            if (stack_.size() == 1)
                found = stack_.back().selectChildFaceOf(fatherFace_.element(), fatherFace_.index());
            else if (stack_.size() > 1)
                found = stack_.back().selectChildFaceOf(*stack_[stack_.size()-2].it_, stack_[stack_.size()-2].index_);
            if (found)
                face_ = stack_.back().toFace(grid_);
            return found;
        }
 
        // encapsulate hierarchic iterator, end iterator and face index
        struct ChildFaceIterator
        {
            ChildFaceIterator(
                const HierarchicElementIterator& it,
                const HierarchicElementIterator& end,
                int face = -1) :
                it_(it),
                end_(end),
                index_(face)
            {}
 
            FaceType toFace(const GridType& grid) const
            {
                return FaceType(grid, *it_, index_);
            }
 
            bool selectChildFaceOf(const Element& father, int fatherFace)
            {
                index_ = FaceHierarchy<GridType>::findSubFace(father, fatherFace, *it_);
                return (index_ >= 0);
            }
 
            HierarchicElementIterator it_;
            HierarchicElementIterator end_;
            int index_;
        };
 
 
        // guess what
        const GridType& grid_;
 
        // maximal level to consider
        int maxLevel_;
 
        // the father face we start from
        const FaceType& fatherFace_;
 
        // the face the iterator points to
        FaceType face_;
 
        // a stack of hierarchic iterators with face indices
        std::vector<ChildFaceIterator> stack_;
};
 
#endif