hierarchicleafiterator.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 HIERARCHIC_LEAF_ITERATOR_HH
#define HIERARCHIC_LEAF_ITERATOR_HH
 
#include <dune/common/fvector.hh>
#include <dune/common/iteratorfacades.hh>
 
#include <dune/geometry/multilineargeometry.hh>
#include <dune/geometry/referenceelements.hh>
 
template <class GridImp>
class HierarchicLeafIterator :
    public Dune::ForwardIteratorFacade<HierarchicLeafIterator<GridImp>, const typename GridImp::template Codim<0>::Entity>
{
    typedef typename GridImp::template Codim<0>::Entity Element;
    typedef typename GridImp::HierarchicIterator HierarchicIterator;
    enum {dim = GridImp::dimension};
    enum {dimworld = GridImp::dimensionworld};
public:
    typedef Dune::CachedMultiLinearGeometry<typename GridImp::ctype, dim, dimworld> LocalGeometry;
 
    enum PositionFlag {begin, end};
 
    HierarchicLeafIterator(const GridImp& grid, const Element& element,
                            PositionFlag flag, bool nested = true)
        : element_(element), maxlevel_(grid.maxLevel()), hIt_(element.hend(maxlevel_)),
            flag_(flag), nested_(nested)
    {
 
        // if the element itself is leaf, then we don't have to iterate over the children
        if (flag==begin && !element_.isLeaf()) {
            hIt_ = element_.hbegin(maxlevel_);
 
            //NOTE This class by now assumes that possible non-nestedness of the grid levels only arises
            // due to boundary parametrisation
            if (!nested_)  {
                // Check if the element is a boundary element, and set the nested flag correspondingly
                // If the element is not at the boundary, then we can neglect the nested flag
                typename GridImp::LevelGridView levelView = grid.levelGridView(element_.level());
                nested_ = true;
 
                for (const auto lIt : intersections(levelView,element_))
                    if (lIt.boundary()) {
                        nested_ = false;
                        break;
                    }
            }
 
            if (!hIt_->isLeaf())
                increment();
            else
              leafChild_ = *hIt_;
        }
    }
 
    bool equals(const HierarchicLeafIterator<GridImp>& other) const {
        return  (element_ == other.element_)  &&
                ((flag_==other.flag_ && flag_==end) || (hIt_ == other.hIt_ && flag_==begin && other.flag_==flag_));
    }
 
    void increment() {
 
        HierarchicIterator hEndIt = element_.hend(maxlevel_);
 
        if (hIt_ == hEndIt) {
            flag_ = end;
            return;
        }
 
        // Increment until we hit a leaf child element
        do {
            ++hIt_;
        } while(hIt_ != hEndIt && (!hIt_->isLeaf()));
 
        if (hIt_ == hEndIt)
            flag_ = end;
        else
            leafChild_ = *hIt_;
    }
 
    const Element& dereference() const {
        if (flag_ == end)
            DUNE_THROW(Dune::Exception,"HierarchicLeafIterator: Cannot dereference end iterator!");
        return (element_.isLeaf()) ? element_ : leafChild_;
    }
 
    LocalGeometry geometryInAncestor() {
 
        //NOTE: We assume here that the type of the child and the ancestors are the same!
        auto ref = Dune::ReferenceElements<typename GridImp::ctype,dim>::general(element_.type());
 
        // store local coordinates of the leaf child element within the coarse starting element
        std::vector<Dune::FieldVector<typename GridImp::ctype,dim> > corners(ref.size(dim));
        for (int i=0; i<corners.size(); i++)
            corners[i] = ref.position(i,dim);
 
        // If the element is leaf, then return an Identity mapping
        if (element_.isLeaf())
            return LocalGeometry(ref,corners);
 
        if (nested_) {
            const typename Element::Geometry leafGeom = leafChild_->geometry();
            const typename Element::Geometry coarseGeom = element_.geometry();
 
            for (int i=0; i<corners.size();i++)
                corners[i] = coarseGeom.local(leafGeom.corner(i));
            return LocalGeometry(ref,corners);
        }
 
        Element father(leafChild_);
        while (father != element_) {
 
            const typename Element::LocalGeometry fatherGeom = father.geometryInFather();
            father = father->father();
 
            for (int i=0; i<corners.size(); i++)
                corners[i] = fatherGeom.global(corners[i]);
        }
 
        return LocalGeometry(ref,corners);
    }
 
private:
 
    const Element& element_;
 
    int maxlevel_;
 
    HierarchicIterator hIt_;
 
    Element leafChild_;
 
    PositionFlag flag_;
 
    bool nested_;
 
};
#endif