1#ifndef DUNE_FEM_MARKING_MAXIMUM_HH
2#define DUNE_FEM_MARKING_MAXIMUM_HH
8#include <dune/grid/common/rangegenerators.hh>
21 namespace GridAdaptation
46 template <
class Gr
id,
class Indicator>
50 const double theta,
int maxLevel = -1 )
54 typedef Dune::ReferenceElements< typename Grid::ctype, Grid::dimension > ReferenceElements;
56 typename Indicator::RangeType value;
58 double totalError( 0 ), maxError( 0 );
59 for (
const auto &e : indicator.space())
61 if (!e.isLeaf())
continue;
62 localIndicator.bind(e);
63 const auto ¢er = ReferenceElements::general( e.type() ).position(0,0);
64 localIndicator.evaluate(center,value);
67 maxError =
max( maxError, eta );
69 maxError = grid.comm().max( maxError );
70 totalError = grid.comm().sum( totalError );
77 std::pair< double, double > a( 0, totalError ), b( maxError, 0 );
78 double m = (a.first + b.first) /
double( 2 );
79 while( (m > a.first) && (a.second > theta*totalError) )
83 std::pair< double, double > c( 0, 0 );
84 for (
const auto &e : indicator.space())
86 if (!e.isLeaf())
continue;
87 localIndicator.bind(e);
88 const auto ¢er = ReferenceElements::general( e.type() ).position(0,0);
89 localIndicator.evaluate(center,value);
90 double eta = value[0];
94 c.first =
max( c.first, eta );
96 c.first = grid.comm().max( c.first );
97 c.second = grid.comm().sum( m );
99 if( a.second > c.second )
102 if( c.second < theta*totalError )
105 a = std::make_pair( m, c.second );
106 m = (a.first + b.first) /
double( 2 );
116 for (
const auto &e : indicator.space())
118 if (!e.isLeaf())
continue;
119 localIndicator.bind(e);
120 const auto ¢er = ReferenceElements::general( e.type() ).position(0,0);
121 localIndicator.evaluate(center,value);
122 double eta = value[0];
126 m = grid.comm().max( m );
132 std::size_t marked = 0;
133 for (
const auto &e : indicator.space())
135 if (!e.isLeaf())
continue;
137 localIndicator.bind(e);
138 const auto ¢er = ReferenceElements::general( e.type() ).position(0,0);
139 localIndicator.evaluate(center,value);
140 double eta = value[0];
144 if (e.level()<maxLevel || maxLevel==-1)
153 return std::make_pair(marked,0);
166 template <
class Gr
id,
class Indicator>
167 static inline std::pair< double, double >
168 computeSumAndMaxGridWalk( Grid& grid,
const Indicator& indicator,
169 const double nu, std::vector< double >& buckets )
171 typedef Dune::ReferenceElements< typename Grid::ctype, Grid::dimension > ReferenceElements;
173 typename Indicator::RangeType value;
174 double maxIndicator = 0;
175 double sumIndicator = 0;
176 for (
const auto &e : indicator.space())
178 if (!e.isLeaf())
continue;
179 localIndicator.bind(e);
180 const auto ¢er = ReferenceElements::general( e.type() ).position(0,0);
181 localIndicator.evaluate(center,value);
182 double eta = value[0];
183 maxIndicator =
std::max(maxIndicator,eta);
188 maxIndicator = indicator.space().gridPart().comm().max( maxIndicator );
189 sumIndicator = indicator.space().gridPart().comm().sum( sumIndicator );
191 for (
const auto &e : indicator.space())
193 if (!e.isLeaf())
continue;
194 localIndicator.bind(e);
195 const auto ¢er = ReferenceElements::general( e.type() ).position(0,0);
196 localIndicator.evaluate(center,value);
197 double eta = value[0];
198 int index = int(eta/maxIndicator*1./nu);
200 assert( index < buckets.size() );
201 buckets[index] += eta;
205 indicator.space().gridPart().comm().sum( buckets.data(), buckets.size() );
207 return std::make_pair( sumIndicator, maxIndicator );
210 template <
class Gr
id,
class Indicator>
211 static inline std::pair< double, double >
212 computeSumAndMax( Grid& grid,
const Indicator& indicator,
213 const double nu, std::vector< double >& buckets )
215 return computeSumAndMaxGridWalk( grid, indicator, nu, buckets );
218 template <
class Gr
id,
class Imp>
219 static inline std::pair< double, double >
221 const double nu, std::vector< double >& buckets )
224 if( indicator.
space().order() > 0 )
225 return computeSumAndMaxGridWalk( grid, indicator, nu, buckets );
227 double maxIndicator = 0;
228 double sumIndicator = 0;
233 maxIndicator =
std::max(maxIndicator,d);
238 maxIndicator = indicator.
space().gridPart().comm().max( maxIndicator );
239 sumIndicator = indicator.
space().gridPart().comm().sum( sumIndicator );
246 int index = int(d/maxIndicator*1./nu);
247 assert( index < buckets.size() );
252 indicator.
space().gridPart().comm().sum( buckets.data(), buckets.size() );
254 return std::make_pair( sumIndicator, maxIndicator );
258 template <
class Gr
id,
class Indicator>
262 const double tolerance,
int maxLevel = -1,
266 maxLevel = ( maxLevel < 0 ) ? std::numeric_limits< int >::max() : maxLevel;
270 typedef Dune::ReferenceElements< typename Grid::ctype, Grid::dimension > ReferenceElements;
272 typename Indicator::RangeType value;
277 std::vector<double> buckets(std::ceil(1./nu)+1);
278 auto sumMax = detail::computeSumAndMax( grid, indicator, nu, buckets );
279 double sumIndicator = sumMax.first;
280 double maxIndicator = sumMax.second;
286 for (
int index=buckets.size()-1;
287 index >= 0 && sum < (1-tolerance)*(1-tolerance)*sumIndicator;
290 sum += buckets[index];
298 const double gammaMaxIndicator = gamma*maxIndicator ;
299 for (
const auto &e : indicator.space())
301 if (!e.isLeaf())
continue;
303 localIndicator.bind(e);
304 const auto ¢er = ReferenceElements::general( e.type() ).position(0,0);
305 localIndicator.evaluate(center,value);
306 double eta = value[0];
307 if (eta > gammaMaxIndicator )
309 if (e.level()<maxLevel)
320 assert( sum >= (1-tolerance)*(1-tolerance)*sumIndicator);
321 return std::make_pair( grid.comm().sum(refMarked), grid.comm().sum(crsMarked) );
Dune::Fem::Double abs(const Dune::Fem::Double &a)
Definition: double.hh:942
double max(const Dune::Fem::Double &v, const double p)
Definition: double.hh:965
Definition: bindguard.hh:11
const GridEntityAccess< Entity >::GridEntityType & gridEntity(const Entity &entity)
Definition: gridpart.hh:412
IteratorRange< typename DF::DofIteratorType > dofs(DF &df)
Iterates over all DOFs.
Definition: rangegenerators.hh:76
typename Impl::ConstLocalFunction< GridFunction >::Type ConstLocalFunction
Definition: const.hh:604
static double max(const Double &v, const double p)
Definition: double.hh:398
static std::pair< int, int > layeredDoerflerMarking(Grid &grid, const Indicator &indicator, const double tolerance, int maxLevel=-1, double nu=0.05)
Definition: doerfler.hh:261
static std::pair< int, int > doerflerMarking(Grid &grid, const Indicator &indicator, const double theta, int maxLevel=-1)
doerflerMarking
Definition: doerfler.hh:49
Definition: common/discretefunction.hh:86
const DiscreteFunctionSpaceType & space() const
obtain a reference to the corresponding DiscreteFunctionSpace
Definition: common/discretefunction.hh:214
static const int keep
Definition: marking/default.hh:19
static const int refine
Definition: marking/default.hh:18