[Dune] Maximal level difference
Carsten Gräser
graeser at math.fu-berlin.de
Mon Nov 29 15:46:05 CET 2010
Hi Christian,
Am 29.11.2010 15:28, schrieb Christian Engwer:
> 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.
it's not intended to restrict the functionality. Currently you can only
have maxLevelDifference=\infty in UGGrid. The function would allow
to be more restrictive if you like.
Best,
Carsten
>
> 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);
>> }
>> }
>> }
>> }
>> }
>> }
More information about the Dune
mailing list