diff --git a/dune/fufem/Makefile.am b/dune/fufem/Makefile.am index b4134e9483878ae4189ab5e2c2192960fd5109f5..4f3f14d28bacd53ecc842070d97385a095b7a4bc 100644 --- a/dune/fufem/Makefile.am +++ b/dune/fufem/Makefile.am @@ -12,6 +12,7 @@ ag_common_HEADERS = arcofcircle.hh discretizationerror.hh \ dgindexset.hh optms.hh surfmassmatrix.hh \ boundarywriter.hh differencenormsquared.hh improvegrid.hh \ orientedsubface.hh symmetrictensor.hh ulis_tools.hh \ - facehierarchy.hh any.hh + facehierarchy.hh any.hh \ + globalintersectioniterator.hh boundaryiterator.hh include $(top_srcdir)/am/global-rules diff --git a/dune/fufem/boundaryiterator.hh b/dune/fufem/boundaryiterator.hh new file mode 100644 index 0000000000000000000000000000000000000000..85da51f42675a5208a53791f4a1c44c63bac8895 --- /dev/null +++ b/dune/fufem/boundaryiterator.hh @@ -0,0 +1,66 @@ +#ifndef BOUNDARY_ITERATOR_HH +#define BOUNDARY_ITERATOR_HH + +#include "globalintersectioniterator.hh" + + + +/** \brief Iterator of boundary intersections of grid view + * + * \tparam GridView The grid view on which this boundary patch lives +*/ +template <class GridView> +class BoundaryIterator + : public GlobalIntersectionIterator<GridView, BoundaryIterator<GridView> > +{ + typedef GlobalIntersectionIterator<GridView, BoundaryIterator<GridView> > Base; + enum {dim=GridView::dimension}; + + typedef typename GridView::IntersectionIterator IntersectionIterator; + typedef typename GridView::template Codim<0>::Iterator ElementIterator; + typedef typename GridView::template Codim<0>::Entity Element; + +public: + typedef typename Base::PositionFlag PositionFlag; + + /** \brief The type of objects we are iterating over */ + typedef typename Base::Intersection Intersection; + + /** \brief Create begin or end iterator for given grid view + * + * \param gridView Iterate over the intersections of this grid view + * \param flag Create begin or end iterator for PositionFlag = begin or end, respectively. + */ + BoundaryIterator(const GridView& gridView, PositionFlag flag) : + Base(gridView, flag) + { + Base::initialIncrement(); + } + + /** \brief Create iterator for given grid view + * + * \param gridView Iterate over the intersections of this grid view + * \param eIt Element iterator to the first element to consider + * \param endEIt Element iterator after the last element to consider + */ + BoundaryIterator(const GridView& gridView, const ElementIterator& eIt, const ElementIterator& endEIt) : + Base(gridView, eIt, endEIt) + { + Base::initialIncrement(); + } + + bool skipElement(const Element& e) + { + return false; + } + + bool skipIntersection(const Intersection& i) + { + return not(i.boundary()); + } +}; + + + +#endif + diff --git a/dune/fufem/boundarypatch.hh b/dune/fufem/boundarypatch.hh index 11c466f8df12980141531d024434d6c4d493b953..99c95a05c53a6e46ea027e2c07528a45ee12bbb1 100644 --- a/dune/fufem/boundarypatch.hh +++ b/dune/fufem/boundarypatch.hh @@ -5,7 +5,6 @@ #include <set> #include <memory> -#include <dune/common/iteratorfacades.hh> #include <dune/common/bitsetvector.hh> #include <dune/grid/common/genericreferenceelements.hh> @@ -13,269 +12,13 @@ #include <dune/grid/common/gridview.hh> #include <dune/grid/common/gridfactory.hh> -// Forward declaration -template <class GridView> class BoundaryPatchIterator; - - - - - -/** \brief Base class for iterators on intersections of a grid view - * - * \tparam GridView The grid view on which this boundary patch lives - * \tparam Impl The Implementation derived from this class (CRTP). - * - * The derived class Impl has to provide to methods - * - * bool skipElement(const Element& e) - * bool skipIntersection(const Intersection& e) - * - * that specify which elements and intersections should be skipped. - * It must call intitialIncrement() in its constructor in order - * to find the first valid intersection. - */ -template <class GridView, class Impl> -class GlobalIntersectionIterator - : public Dune::ForwardIteratorFacade<GlobalIntersectionIterator<GridView, Impl>, const typename GridView::Intersection> -{ - enum {dim=GridView::dimension}; - - typedef typename GridView::IntersectionIterator IntersectionIterator; - typedef typename GridView::template Codim<0>::Iterator ElementIterator; - -public: - - enum PositionFlag {begin, end}; +#include "globalintersectioniterator.hh" +#include "boundaryiterator.hh" - /** \brief The type of objects we are iterating over - */ - typedef typename IntersectionIterator::Intersection Intersection; - - /** \brief Create begin or end iterator for given grid view - * - * The iterator does not necessarily point to a valid intersection - * after creation! In order to guarantee this you must - * call initialIncrement() from the constructor in the - * derived class. This can not be done here since it - * requires the skipXYZ() methods of the derived class - * that is not yet initialized. - * - * \param gridView Iterate over the intersections of this grid view - * \param flag Create begin or end iterator for PositionFlag = begin or end, respectively. - */ - GlobalIntersectionIterator(const GridView& gridView, PositionFlag flag) : - gridView_(gridView), - eIt_(gridView.template end<0>()), - endEIt_(gridView.template end<0>()), - nIt_(0) - { - if (flag==begin) - { - eIt_ = gridView.template begin<0>(); - nIt_ = new IntersectionIterator(gridView.ibegin(*eIt_)); - } - } - - /** \brief Create iterator for given grid view - * - * The iterator does not necessarily point to a valid intersection - * after creation! In order to guarantee this you must - * call initialIncrement() from the constructor in the - * derived class. This can not be done here since it - * requires the skipXYZ() methods of the derived class - * that is not yet initialized. - * - * \param gridView Iterate over the intersections of this grid view - * \param eIt Element iterator to the first element to consider - * \param endEIt Element iterator after the last element to consider - */ - GlobalIntersectionIterator(const GridView& gridView, const ElementIterator& eIt, const ElementIterator& endEIt) : - gridView_(gridView), - eIt_(eIt), - endEIt_(endEIt), - nIt_(0) - { - if (eIt_!=endEIt_) - nIt_ = new IntersectionIterator(gridView.ibegin(*eIt_)); - } - - /** \brief Copy constructor - * - * \param other Copy from this iterator - */ - GlobalIntersectionIterator(const GlobalIntersectionIterator& other) : - gridView_(other.gridView_), - eIt_(other.eIt_), - endEIt_(other.endEIt_), - nIt_(0) - - { - if (other.nIt_) - nIt_ = new IntersectionIterator(*other.nIt_); - } - - ~GlobalIntersectionIterator() { - delete(nIt_); - } - - /** \brief Increment the pointer to the next intersection - */ - void increment() - { - // Increment until a face in the patch is found - do { - ++(*nIt_); - - // if end intersection reached find next intersection - if ((*nIt_) == gridView_.iend(*eIt_)) - { - do { - ++eIt_; - } while (eIt_ != endEIt_ and asImpl().skipElement(*eIt_)); - if (eIt_ != endEIt_) - (*nIt_) = gridView_.ibegin(*eIt_); - } - - } while (eIt_ != endEIt_ and asImpl().skipIntersection(**nIt_)); - } - - bool equals(const GlobalIntersectionIterator& other) const - { - return (eIt_ == other.eIt_) and ((eIt_ == endEIt_) or ((*nIt_) == (*other.nIt_))); - } - - const Intersection& dereference() const - { - return **nIt_; - } - - const GlobalIntersectionIterator& operator=(const GlobalIntersectionIterator& other) - { - eIt_ = other.eIt_; - endEIt_ = other.endEIt_; - - if (nIt_) - delete nIt_; - - if (other.nIt_) - nIt_ = new IntersectionIterator(*other.nIt_); - else - nIt_ = 0; - - return *this; - } - - bool containsInsideSubentity(int subEntity, int codim) const - { - if (codim==0) // disregard element dofs - return false; - int faceIdx = (*nIt_)->indexInInside(); - - if (codim==1) - return faceIdx == subEntity; - - const Dune::GenericReferenceElement<double,dim>& refElement - = Dune::GenericReferenceElements<double, dim>::general(eIt_->type()); - - for (int j = 0; j<refElement.size(faceIdx, 1, codim); j++) - { - if (refElement.subEntity(faceIdx, 1, j, codim) == subEntity) - return true; - } - return false; - } - -protected: - - /** \brief Find first valid intersection - * - * Checks if current intersection is valid and iterates - * to the first valid intersection if not. - * You must call this method in the constructor of the - * derived class! - */ - void initialIncrement() - { - if (eIt_ != endEIt_ and asImpl().skipIntersection(**nIt_)) - increment(); - } - - const Impl& asImpl() const - { - return static_cast<const Impl&>(*this); - } - - Impl& asImpl() - { - return static_cast<Impl&>(*this); - } - - ElementIterator eIt_; - ElementIterator endEIt_; - - IntersectionIterator* nIt_; - - GridView gridView_; -}; - - - -/** \brief Iterator of boundary intersections of grid view - * - * \tparam GridView The grid view on which this boundary patch lives -*/ -template <class GridView> -class BoundaryIterator - : public GlobalIntersectionIterator<GridView, BoundaryIterator<GridView> > -{ - typedef GlobalIntersectionIterator<GridView, BoundaryIterator<GridView> > Base; - enum {dim=GridView::dimension}; - - typedef typename GridView::IntersectionIterator IntersectionIterator; - typedef typename GridView::template Codim<0>::Iterator ElementIterator; - typedef typename GridView::template Codim<0>::Entity Element; - -public: - typedef typename Base::PositionFlag PositionFlag; - - /** \brief The type of objects we are iterating over */ - typedef typename Base::Intersection Intersection; - - /** \brief Create begin or end iterator for given grid view - * - * \param gridView Iterate over the intersections of this grid view - * \param flag Create begin or end iterator for PositionFlag = begin or end, respectively. - */ - BoundaryIterator(const GridView& gridView, PositionFlag flag) : - Base(gridView, flag) - { - Base::initialIncrement(); - } - - /** \brief Create iterator for given grid view - * - * \param gridView Iterate over the intersections of this grid view - * \param eIt Element iterator to the first element to consider - * \param endEIt Element iterator after the last element to consider - */ - BoundaryIterator(const GridView& gridView, const ElementIterator& eIt, const ElementIterator& endEIt) : - Base(gridView, eIt, endEIt) - { - Base::initialIncrement(); - } - - bool skipElement(const Element& e) - { - return false; - } - - bool skipIntersection(const Intersection& i) - { - return not(i.boundary()); - } -}; +// Forward declaration +template <class GridView> class BoundaryPatchIterator; diff --git a/dune/fufem/globalintersectioniterator.hh b/dune/fufem/globalintersectioniterator.hh new file mode 100644 index 0000000000000000000000000000000000000000..a1d10dfd3dc6336ff7a0c426e2ec9c4fd9faa0e7 --- /dev/null +++ b/dune/fufem/globalintersectioniterator.hh @@ -0,0 +1,212 @@ +#ifndef GLOBAL_INTERSECTION_ITERATOR_HH +#define GLOBAL_INTERSECTION_ITERATOR_HH + +#include <dune/common/iteratorfacades.hh> + +#include <dune/grid/common/genericreferenceelements.hh> + + + +/** \brief Base class for iterators on intersections of a grid view + * + * \tparam GridView The grid view on which this boundary patch lives + * \tparam Impl The Implementation derived from this class (CRTP). + * + * The derived class Impl has to provide to methods + * + * bool skipElement(const Element& e) + * bool skipIntersection(const Intersection& e) + * + * that specify which elements and intersections should be skipped. + * It must call intitialIncrement() in its constructor in order + * to find the first valid intersection. + */ +template <class GridView, class Impl> +class GlobalIntersectionIterator + : public Dune::ForwardIteratorFacade<GlobalIntersectionIterator<GridView, Impl>, const typename GridView::Intersection> +{ + enum {dim=GridView::dimension}; + + typedef typename GridView::IntersectionIterator IntersectionIterator; + typedef typename GridView::template Codim<0>::Iterator ElementIterator; + +public: + + enum PositionFlag {begin, end}; + + /** \brief The type of objects we are iterating over + */ + typedef typename IntersectionIterator::Intersection Intersection; + + /** \brief Create begin or end iterator for given grid view + * + * The iterator does not necessarily point to a valid intersection + * after creation! In order to guarantee this you must + * call initialIncrement() from the constructor in the + * derived class. This can not be done here since it + * requires the skipXYZ() methods of the derived class + * that is not yet initialized. + * + * \param gridView Iterate over the intersections of this grid view + * \param flag Create begin or end iterator for PositionFlag = begin or end, respectively. + */ + GlobalIntersectionIterator(const GridView& gridView, PositionFlag flag) : + gridView_(gridView), + eIt_(gridView.template end<0>()), + endEIt_(gridView.template end<0>()), + nIt_(0) + { + if (flag==begin) + { + eIt_ = gridView.template begin<0>(); + nIt_ = new IntersectionIterator(gridView.ibegin(*eIt_)); + } + } + + /** \brief Create iterator for given grid view + * + * The iterator does not necessarily point to a valid intersection + * after creation! In order to guarantee this you must + * call initialIncrement() from the constructor in the + * derived class. This can not be done here since it + * requires the skipXYZ() methods of the derived class + * that is not yet initialized. + * + * \param gridView Iterate over the intersections of this grid view + * \param eIt Element iterator to the first element to consider + * \param endEIt Element iterator after the last element to consider + */ + GlobalIntersectionIterator(const GridView& gridView, const ElementIterator& eIt, const ElementIterator& endEIt) : + gridView_(gridView), + eIt_(eIt), + endEIt_(endEIt), + nIt_(0) + { + if (eIt_!=endEIt_) + nIt_ = new IntersectionIterator(gridView.ibegin(*eIt_)); + } + + /** \brief Copy constructor + * + * \param other Copy from this iterator + */ + GlobalIntersectionIterator(const GlobalIntersectionIterator& other) : + gridView_(other.gridView_), + eIt_(other.eIt_), + endEIt_(other.endEIt_), + nIt_(0) + + { + if (other.nIt_) + nIt_ = new IntersectionIterator(*other.nIt_); + } + + ~GlobalIntersectionIterator() { + delete(nIt_); + } + + /** \brief Increment the pointer to the next intersection + */ + void increment() + { + // Increment until a face in the patch is found + do { + ++(*nIt_); + + // if end intersection reached find next intersection + if ((*nIt_) == gridView_.iend(*eIt_)) + { + do { + ++eIt_; + } while (eIt_ != endEIt_ and asImpl().skipElement(*eIt_)); + if (eIt_ != endEIt_) + (*nIt_) = gridView_.ibegin(*eIt_); + } + + } while (eIt_ != endEIt_ and asImpl().skipIntersection(**nIt_)); + } + + bool equals(const GlobalIntersectionIterator& other) const + { + return (eIt_ == other.eIt_) and ((eIt_ == endEIt_) or ((*nIt_) == (*other.nIt_))); + } + + const Intersection& dereference() const + { + return **nIt_; + } + + const GlobalIntersectionIterator& operator=(const GlobalIntersectionIterator& other) + { + eIt_ = other.eIt_; + endEIt_ = other.endEIt_; + + if (nIt_) + delete nIt_; + + if (other.nIt_) + nIt_ = new IntersectionIterator(*other.nIt_); + else + nIt_ = 0; + + return *this; + } + + bool containsInsideSubentity(int subEntity, int codim) const + { + if (codim==0) // disregard element dofs + return false; + + int faceIdx = (*nIt_)->indexInInside(); + + if (codim==1) + return faceIdx == subEntity; + + const Dune::GenericReferenceElement<double,dim>& refElement + = Dune::GenericReferenceElements<double, dim>::general(eIt_->type()); + + for (int j = 0; j<refElement.size(faceIdx, 1, codim); j++) + { + if (refElement.subEntity(faceIdx, 1, j, codim) == subEntity) + return true; + } + return false; + } + +protected: + + /** \brief Find first valid intersection + * + * Checks if current intersection is valid and iterates + * to the first valid intersection if not. + * You must call this method in the constructor of the + * derived class! + */ + void initialIncrement() + { + if (eIt_ != endEIt_ and asImpl().skipIntersection(**nIt_)) + increment(); + } + + const Impl& asImpl() const + { + return static_cast<const Impl&>(*this); + } + + Impl& asImpl() + { + return static_cast<Impl&>(*this); + } + + ElementIterator eIt_; + ElementIterator endEIt_; + + IntersectionIterator* nIt_; + + GridView gridView_; +}; + + + +#endif +