[Dune] Maximal level difference
Christian Engwer
christi at uni-hd.de
Mon Nov 29 15:28:45 CET 2010
Hi Carsten,
I think it is good to add this code to the utilities. I didn't yet
look at the actual implementation, though.
Regarding the implementation in UG, I always like this feature and
used it. I really want to be able to have arbitrarily hanging nodes,
well actually I want to have control on this in my program and not let
the grid decide. So I hope that the current behaviour of UGGrid will
still be available.
Ciao
Christian
On Mon, Nov 29, 2010 at 02:44:10PM +0100, Carsten Gräser wrote:
> Sorry, I forgot another issue regarding this:
>
> Am 29.11.2010 14:35, schrieb Carsten Gräser:
> > Dear dune,
> > in dune-subgrid we implemented a simple method to automatically
> > add refinement marks and remove coarsening marks such that
> > the level of neighboring elements differs at most by a user
> > defined number. (ALUGrid e.g. does this internally).
> >
> > As someone asked for this functionality in UGGrid I extracted
> > this using the grid interface only. This might also be helpfull
> > for other third party grid implementations.
> >
> > Are there any objections to add this to dune/grid/utility?
>
> We could also add a method
>
> setMaxLevelDifference(unsigned int level)
>
> to UGGrid and call the utility-function from within
> UGGrid::preAdapt() internally. On the one hand this
> would make it a little easier to write portable code,
> on the other hand it's essentially syntactic candy
> that blows up the UGGrid interface (by one method).
>
> Furthermore I'm not sure if the name is good. Currently
> the function is called
>
> template<class GridType>
> static void enforceNeighborLevelDifference(GridType& grid, int maxLevelDiff);
>
> Any opinions?
>
> Regards,
> Carsten
>
> PS: You can find the full implementation attached.
>
> --
> ----------------------------------------------------------------------
> Carsten Gräser | phone: +49-30 / 838-75349
> Freie Universität Berlin | fax : +49-30 / 838-54977
> Institut für Mathematik | email: graeser at math.fu-berlin.de
> Arnimallee 6 |
> 14195 Berlin, Germany | URL : http://page.mi.fu-berlin.de/graeser
> ----------------------------------------------------------------------
> #ifndef ENFORCE_NEIGHBOR_LEVEL_DIFFERENCE_HH
> #define ENFORCE_NEIGHBOR_LEVEL_DIFFERENCE_HH
>
>
>
> template<class GridType>
> class VertexNeighborLevelProvider
> {
> static const int dim=GridType::dimension;
> typedef typename GridType::LeafGridView GridView;
> typedef typename GridType::LeafIndexSet IndexSet;
> typedef typename GridType::template Codim<0>::Entity Element;
>
> public:
>
> VertexNeighborLevelProvider(const GridType& grid) :
> grid_(grid),
> gridView_(grid.leafView()),
> indexSet_(gridView_.indexSet()),
> maxElementLevel_(grid.size(dim), 0)
> {
> typedef typename GridType::template Codim<0>::LeafIterator ElementLeafIterator;
> ElementLeafIterator it = gridView_.template begin<0>();
> ElementLeafIterator end = gridView_.template end<0>();
> for (; it!=end; ++it)
> setElementLevel(*it, it->level() + grid_.getMark(*it));
> }
>
> char getMaxNeighborLevel(const Element& e) const
> {
> char maxNeighborLevel = 0;
> for (int i = 0; i<e.template count<dim>(); ++i)
> {
> int nodeIndex = indexSet_.subIndex(e, i, dim);
> if (maxNeighborLevel < maxElementLevel_[nodeIndex])
> maxNeighborLevel = maxElementLevel_[nodeIndex];
> }
> return maxNeighborLevel;
> }
>
> void setElementLevel(const Element& e, char newLevel)
> {
> for (int i = 0; i<e.template count<dim>(); ++i)
> {
> int nodeIndex = indexSet_.subIndex(e, i, dim);
> if (maxElementLevel_[nodeIndex] < newLevel)
> maxElementLevel_[nodeIndex] = newLevel;
> }
> }
>
> private:
> const GridType& grid_;
> const GridView gridView_;
> const IndexSet& indexSet_;
> std::vector<char> maxElementLevel_;
> };
>
>
>
> template<class GridType>
> class IntersectionNeighborLevelProvider
> {
> static const int dim=GridType::dimension;
> typedef typename GridType::LeafGridView GridView;
> typedef typename GridType::LeafIndexSet IndexSet;
> typedef typename GridType::template Codim<0>::Entity Element;
> typedef typename GridType::template Codim<0>::EntityPointer ElementPointer;
>
> public:
>
> IntersectionNeighborLevelProvider(const GridType& grid) :
> grid_(grid)
> {}
>
> char getMaxNeighborLevel(const Element& e) const
> {
> typename GridType::LeafIntersectionIterator it = e.ileafbegin();
> typename GridType::LeafIntersectionIterator end = e.ileafend();
> char maxNeighborLevel = 0;
> for(; it!=end; ++it)
> {
> if (not(it->boundary()))
> {
> ElementPointer neighbor = it->outside();
>
> char level = neighbor->level() + grid_.getMark(*neighbor);
> maxNeighborLevel = (maxNeighborLevel < level) ? level : maxNeighborLevel;
> }
> }
> return maxNeighborLevel;
> }
>
> void setElementLevel(const Element& e, char newLevel)
> {}
>
> private:
> const GridType& grid_;
> };
>
>
>
> enum NeighborLevelDifferenceMode {VertexMode, IntersectionMode};
>
>
>
> /** \brief Enforce maximal level difference of neighboring elements
> *
> * This method adds additional refinement marks or removes coarsening marks
> * to ensure that the maximal level difference is respected after the next
> * adaptation.
> *
> * \tparam GridType Type of the grid to consider
> * \param grid Set and remove marks with respect to this grid
> * \param maxLevelDiff Maximal allowed level difference of neighboring elements
> * \param mode VertexMode/IntersectionMode ensure maximal level difference
> * of elements sharing a vertex/intersection.
> */
> template<class GridType>
> static void enforceNeighborLevelDifference(GridType& grid, int maxLevelDiff, enum NeighborLevelDifferenceMode mode=IntersectionMode)
> {
> if (mode==ByVertex)
> {
> VertexNeighborLevelProvider<GridType> neighbors(grid);
> enforceNeighborLevelDifference(grid, maxLevelDiff, neighbors);
> }
> else if (mode==ByIntersection)
> {
> IntersectionNeighborLevelProvider<GridType> neighbors(grid);
> enforceNeighborLevelDifference(grid, maxLevelDiff, neighbors);
> }
> }
>
>
>
> template<class GridType, class NeighborLevelProvider>
> static void enforceNeighborLevelDifference(GridType& grid, int maxLevelDiff, NeighborLevelProvider& neighborLevelProvider)
> {
> typedef typename GridType::template Codim<0>::LevelIterator ElementLevelIterator;
> typedef typename GridType::template Codim<0>::LeafIterator ElementLeafIterator;
> typedef typename GridType::template Codim<0>::Entity::HierarchicIterator HierarchicIterator;
> typedef typename GridType::LevelIndexSet LevelIndexSet;
> typedef typename GridType::template Codim<0>::Entity::EntityPointer ElementPointer;
>
> for (int level=grid.maxLevel(); level>=0; --level)
> {
> bool hasImplicitRefinementMarks = false;
> bool hasCoarseningMarks = false;
>
> ElementLevelIterator it = grid.template lbegin<0>(level);
> ElementLevelIterator end = grid.template lend<0>(level);
> for (; it!=end; ++it)
> {
> if(it->isLeaf())
> {
> // compute minimal allowed level for this element
> int minAllowedLevel = neighborLevelProvider.getMaxNeighborLevel(*it) - maxLevelDiff;
>
> int oldMark = grid.getMark(*it);
> int mark = oldMark;
>
> // set refinement mark and remove coarsening mark if neccessary
> if (level+mark < minAllowedLevel)
> {
> mark = minAllowedLevel - level;
> assert((-1 <= mark) and (mark <= 1));
> grid.mark(mark, *it);
> neighborLevelProvider.setElementLevel(*it, level+mark);
> }
>
> if (mark < 0)
> hasCoarseningMarks = true;
>
> if (mark > oldMark)
> hasImplicitRefinementMarks = true;
> }
> }
>
> // if there are coarsening marks left and we have set
> // implicit refinement marks check again if coarsening is possible
> if (hasCoarseningMarks and hasImplicitRefinementMarks)
> {
> it = grid.template lbegin<0>(level);
> end = grid.template lend<0>(level);
> for (; it!=end; ++it)
> {
> if(it->isLeaf())
> {
> // compute minimal allowed level for this element
> int minAllowedLevel = neighborLevelProvider.getMaxNeighborLevel(*it) - maxLevelDiff;
>
> int mark = grid.getMark(*it);
>
> if (level+mark < minAllowedLevel)
> {
> mark = minAllowedLevel - level;
> assert((-1 <= mark) and (mark <= 1));
> grid.mark(mark, *it);
> neighborLevelProvider.setElementLevel(*it, level+mark);
> }
> }
> }
> }
> }
> }
>
>
>
> #endif
>
> _______________________________________________
> Dune mailing list
> Dune at dune-project.org
> http://lists.dune-project.org/mailman/listinfo/dune
More information about the Dune
mailing list