dune-fem  2.8-git
doerfler.hh
Go to the documentation of this file.
1 #ifndef DUNE_FEM_MARKING_MAXIMUM_HH
2 #define DUNE_FEM_MARKING_MAXIMUM_HH
3 
4 #include <cstddef>
5 
6 #include <algorithm>
7 
8 #include <dune/grid/common/rangegenerators.hh>
9 
12 
14 
15 namespace Dune
16 {
17 
18  namespace Fem
19  {
20 
21  namespace GridAdaptation
22  {
23 
24  // doerflerMarking
25  // ---------------
26 
46  template <class Grid, class Indicator>
47  static inline
48  std::pair<int, int>
49  doerflerMarking( Grid& grid, const Indicator& indicator,
50  const double theta, int maxLevel = -1 )
51  {
52  using std::max;
53 
54  typedef Dune::ReferenceElements< typename Grid::ctype, Grid::dimension > ReferenceElements;
55  Dune::Fem::ConstLocalFunction<Indicator> localIndicator(indicator);
56  typename Indicator::RangeType value;
57 
58  double totalError( 0 ), maxError( 0 );
59  for (const auto &e : indicator.space())
60  {
61  if (!e.isLeaf()) continue;
62  localIndicator.bind(e);
63  const auto &center = ReferenceElements::general( e.type() ).position(0,0);
64  localIndicator.evaluate(center,value);
65  double eta = std::abs(value[0]);
66  totalError += eta;
67  maxError = max( maxError, eta );
68  }
69  maxError = grid.comm().max( maxError );
70  totalError = grid.comm().sum( totalError );
71 
72  // Let a.first be a real number and denote by a.second the sum of all
73  // local errors greater than a.first.
74  // Now consider two such pairs, a and b, such that
75  // a.second >= theta*totalError > b.second.
76  // We seek to minimize this interval using bisection.
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) )
80  {
81  // c.first: maximum local error less or equal to m
82  // c.second: sum of local errors greater than m
83  std::pair< double, double > c( 0, 0 );
84  for (const auto &e : indicator.space())
85  {
86  if (!e.isLeaf()) continue;
87  localIndicator.bind(e);
88  const auto &center = ReferenceElements::general( e.type() ).position(0,0);
89  localIndicator.evaluate(center,value);
90  double eta = value[0];
91  if( eta > m )
92  c.second += eta;
93  else
94  c.first = max( c.first, eta );
95  }
96  c.first = grid.comm().max( c.first );
97  c.second = grid.comm().sum( m );
98 
99  if( a.second > c.second )
100  {
101  // note: There is no local error in (m.first, c]
102  if( c.second < theta*totalError )
103  b = c;
104  else
105  a = std::make_pair( m, c.second );
106  m = (a.first + b.first) / double( 2 );
107  }
108  else
109  {
110  // The total error could not be reduced. Is it still possible?
111  // Try to find a local error m in (a.first, b.first) and use it as
112  // next interval splitting point.
113  // Note: If such an m exists, it must be greater than the mid point,
114  // because we already tried that one.
115  m = 0;
116  for (const auto &e : indicator.space())
117  {
118  if (!e.isLeaf()) continue;
119  localIndicator.bind(e);
120  const auto &center = ReferenceElements::general( e.type() ).position(0,0);
121  localIndicator.evaluate(center,value);
122  double eta = value[0];
123  if( eta < b.first )
124  m = max( m, eta );
125  }
126  m = grid.comm().max( m );
127  }
128  }
129 
130  // marking all elements with local error > a.first now yields the desired
131  // property.
132  std::size_t marked = 0;
133  for (const auto &e : indicator.space())
134  {
135  if (!e.isLeaf()) continue;
136 
137  localIndicator.bind(e);
138  const auto &center = ReferenceElements::general( e.type() ).position(0,0);
139  localIndicator.evaluate(center,value);
140  double eta = value[0];
141  if( eta <= a.first )
142  continue;
143  const auto &gridEntity = Dune::Fem::gridEntity(e);
144  if (e.level()<maxLevel || maxLevel==-1)
145  {
146  grid.mark(Marker::refine, gridEntity);
147  ++marked;
148  }
149  else
150  grid.mark(Marker::keep, gridEntity);
151  }
152 
153  return std::make_pair(marked,0);
154  }
155 
156 
164  namespace detail
165  {
166  template <class Grid, class Indicator>
167  static inline std::pair< double, double >
168  computeSumAndMaxGridWalk( Grid& grid, const Indicator& indicator,
169  const double nu, std::vector< double >& buckets )
170  {
171  typedef Dune::ReferenceElements< typename Grid::ctype, Grid::dimension > ReferenceElements;
172  Dune::Fem::ConstLocalFunction<Indicator> localIndicator(indicator);
173  typename Indicator::RangeType value;
174  double maxIndicator = 0;
175  double sumIndicator = 0;
176  for (const auto &e : indicator.space())
177  {
178  if (!e.isLeaf()) continue;
179  localIndicator.bind(e);
180  const auto &center = ReferenceElements::general( e.type() ).position(0,0);
181  localIndicator.evaluate(center,value);
182  double eta = value[0];
183  maxIndicator = std::max(maxIndicator,eta);
184  sumIndicator += eta;
185  }
186 
187  // compute global values
188  maxIndicator = indicator.space().gridPart().comm().max( maxIndicator );
189  sumIndicator = indicator.space().gridPart().comm().sum( sumIndicator );
190 
191  for (const auto &e : indicator.space())
192  {
193  if (!e.isLeaf()) continue;
194  localIndicator.bind(e);
195  const auto &center = 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);
199  // std::cout << " " << eta << " " << maxIndicator << " " << index << std::endl;
200  assert( index < buckets.size() );
201  buckets[index] += eta;
202  }
203 
204  // compute global sum of all buckets in parallel
205  indicator.space().gridPart().comm().sum( buckets.data(), buckets.size() );
206 
207  return std::make_pair( sumIndicator, maxIndicator );
208  }
209 
210  template <class Grid, class Indicator>
211  static inline std::pair< double, double >
212  computeSumAndMax( Grid& grid, const Indicator& indicator,
213  const double nu, std::vector< double >& buckets )
214  {
215  return computeSumAndMaxGridWalk( grid, indicator, nu, buckets );
216  }
217 
218  template <class Grid, class Imp>
219  static inline std::pair< double, double >
220  computeSumAndMax( Grid& grid, const Dune::Fem::DiscreteFunctionInterface< Imp >& indicator,
221  const double nu, std::vector< double >& buckets )
222  {
223  // if space is for some reason not constant then default to general method
224  if( indicator.space().order() > 0 )
225  return computeSumAndMaxGridWalk( grid, indicator, nu, buckets );
226 
227  double maxIndicator = 0;
228  double sumIndicator = 0;
229  // but this way we can avoid grid traversal etc.
230  // at the moment we get a VirtualizedGF so this wouldn't work
231  for (const auto &d : Dune::Fem::dofs(indicator) )
232  {
233  maxIndicator = std::max(maxIndicator,d);
234  sumIndicator += d;
235  }
236 
237  // compute global values
238  maxIndicator = indicator.space().gridPart().comm().max( maxIndicator );
239  sumIndicator = indicator.space().gridPart().comm().sum( sumIndicator );
240 
241  // let's assume that indicator is a FV function (would be good to be able to check this)
242  // but this way we can avoid grid traversal etc.
243  // at the moment we get a VirtualizedGF so this wouldn't work
244  for (const auto &d : Dune::Fem::dofs(indicator) )
245  {
246  int index = int(d/maxIndicator*1./nu);
247  assert( index < buckets.size() );
248  buckets[index] += d;
249  }
250 
251  // compute global sum of all buckets in parallel
252  indicator.space().gridPart().comm().sum( buckets.data(), buckets.size() );
253 
254  return std::make_pair( sumIndicator, maxIndicator );
255  }
256  }
257 
258  template <class Grid, class Indicator>
259  static inline
260  std::pair<int, int>
261  layeredDoerflerMarking( Grid& grid, const Indicator& indicator,
262  const double tolerance, int maxLevel = -1,
263  double nu = 0.05)
264  {
265  // if maxLevel < 0 set to maximum value
266  maxLevel = ( maxLevel < 0 ) ? std::numeric_limits< int >::max() : maxLevel;
267 
268  int refMarked = 0;
269  int crsMarked = 0;
270  typedef Dune::ReferenceElements< typename Grid::ctype, Grid::dimension > ReferenceElements;
271  Dune::Fem::ConstLocalFunction<Indicator> localIndicator(indicator);
272  typename Indicator::RangeType value;
273 
274  // Part 1: compute sum and maximum of indicators
275  // Part 2: subdivide [0,maxEta] into nu sized intervals and compute how
276  // much of the overal error is in each interval
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;
281 
282  // Part 3: compute how many buckets we have to use to get
283  // above (1-tolerance)*(1-tolerance)*sumIndicator:
284  double gamma = 1;
285  double sum = 0;
286  for (int index=buckets.size()-1;
287  index >= 0 && sum < (1-tolerance)*(1-tolerance)*sumIndicator;
288  --index)
289  {
290  sum += buckets[index];
291  gamma -= nu;
292  }
293 
294  // no communication should be required and we can simply now
295  // go ahead with the refinement:
296  //
297  sum = 0; // just for checkint that it worked...
298  const double gammaMaxIndicator = gamma*maxIndicator ;
299  for (const auto &e : indicator.space())
300  {
301  if (!e.isLeaf()) continue;
302  const auto &gridEntity = Dune::Fem::gridEntity(e);
303  localIndicator.bind(e);
304  const auto &center = ReferenceElements::general( e.type() ).position(0,0);
305  localIndicator.evaluate(center,value);
306  double eta = value[0];
307  if (eta > gammaMaxIndicator )
308  {
309  if (e.level()<maxLevel)
310  refMarked += grid.mark(Marker::refine, gridEntity);
311  else
312  grid.mark(Marker::keep, gridEntity);
313  // although we might not have marked for refinement due to
314  // level restriction we count this as part of the indicator
315  // taken care of
316  sum += eta;
317  }
318  }
319  // just checking that the algorihtm did the right thing...
320  assert( sum >= (1-tolerance)*(1-tolerance)*sumIndicator);
321  return std::make_pair( grid.comm().sum(refMarked), grid.comm().sum(crsMarked) );
322  }
323 
324  } // end namespace GridAdaptation
325 
326  } // namespace Fem
327 
328 } // namespace Dune
329 
330 #endif // #ifndef DUNE_FEM_MARKING_MAXIMUM_HH
Definition: bindguard.hh:11
Double abs(const Double &a)
Definition: double.hh:871
typename Impl::ConstLocalFunction< GridFunction >::Type ConstLocalFunction
Definition: const.hh:517
IteratorRange< typename DF::DofIteratorType > dofs(DF &df)
Iterates over all DOFs.
Definition: rangegenerators.hh:76
static double max(const Double &v, const double p)
Definition: double.hh:398
const GridEntityAccess< Entity >::GridEntityType & gridEntity(const Entity &entity)
Definition: gridpart.hh:412
static constexpr T max(T a)
Definition: utility.hh:77
static constexpr std::decay_t< T > sum(T a)
Definition: utility.hh:33
static std::pair< int, int > doerflerMarking(Grid &grid, const Indicator &indicator, const double theta, int maxLevel=-1)
doerflerMarking
Definition: doerfler.hh:49
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
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