dune-pdelab  2.4-dev
adaptivity.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 
4 #ifndef DUNE_PDELAB_ADAPTIVITY_HH
5 #define DUNE_PDELAB_ADAPTIVITY_HH
6 
7 #include<dune/common/exceptions.hh>
8 
9 #include<limits>
10 #include<vector>
11 #include<map>
12 #include<unordered_map>
13 #include<dune/common/dynmatrix.hh>
14 #include<dune/geometry/quadraturerules.hh>
17 
19 // for InterpolateBackendStandard
21 // for intersectionoperator
24 
25 #include<dune/grid/io/file/vtk/subsamplingvtkwriter.hh>
26 
27 namespace Dune {
28  namespace PDELab {
29 
30 
31  template<typename GFS>
33  {
34 
35  typedef typename GFS::Traits::GridView::template Codim<0>::Entity Cell;
37 
38  // we need an additional entry because we store offsets and we also want the
39  // offset after the last leaf for size calculations
40  typedef array<std::size_t,TypeTree::TreeInfo<GFS>::leafCount + 1> LeafOffsets;
41 
42  const LeafOffsets& operator[](GeometryType gt) const
43  {
44  const LeafOffsets& leaf_offsets = _leaf_offset_cache[GlobalGeometryTypeIndex::index(gt)];
45  // make sure we have data for this geometry type
46  assert(leaf_offsets.back() > 0);
47  return leaf_offsets;
48  }
49 
50  void update(const Cell& e)
51  {
52  LeafOffsets& leaf_offsets = _leaf_offset_cache[GlobalGeometryTypeIndex::index(e.type())];
53  if (leaf_offsets.back() == 0)
54  {
55  _lfs.bind(e);
56  extract_lfs_leaf_sizes(_lfs,leaf_offsets.begin()+1);
57  // convert to offsets
58  std::partial_sum(leaf_offsets.begin(),leaf_offsets.end(),leaf_offsets.begin());
59  // sanity check
60  assert(leaf_offsets.back() == _lfs.size());
61  }
62  }
63 
64  explicit LeafOffsetCache(const GFS& gfs)
65  : _lfs(gfs)
66  , _leaf_offset_cache(GlobalGeometryTypeIndex::size(Cell::dimension))
67  {}
68 
69  LFS _lfs;
70  std::vector<LeafOffsets> _leaf_offset_cache;
71 
72  };
73 
74 
75  namespace {
76 
77  template<typename MassMatrices,typename Cell>
78  struct inverse_mass_matrix_calculator
79  : public TypeTree::TreeVisitor
80  , public TypeTree::DynamicTraversal
81  {
82 
83  static const int dim = Cell::Geometry::mydimension;
84  typedef std::size_t size_type;
85  typedef typename MassMatrices::value_type MassMatrix;
86  typedef typename MassMatrix::field_type DF;
87  typedef typename Dune::QuadratureRule<DF,dim>::const_iterator QRIterator;
88 
89  template<typename GFS, typename TreePath>
90  void leaf(const GFS& gfs, TreePath treePath)
91  {
92  typedef typename GFS::Traits::FiniteElementMap FEM;
93  const FEM& fem = gfs.finiteElementMap();
94  typedef typename FEM::Traits::FiniteElement FiniteElement;
95  const FiniteElement& fe = fem.find(_element);
96  size_type local_size = fe.localBasis().size();
97 
98  MassMatrix& mass_matrix = _mass_matrices[_leaf_index];
99  mass_matrix.resize(local_size,local_size);
100 
101 
102  std::vector<typename FiniteElement::Traits::LocalBasisType::Traits::RangeType> phi;
103  phi.resize(std::max(phi.size(),local_size));
104 
105  for (const auto& ip : _quadrature_rule)
106  {
107  fe.localBasis().evaluateFunction(ip.position(),phi);
108  const DF factor = ip.weight();
109 
110  for (size_type i = 0; i < local_size; ++i)
111  for (size_type j = 0; j < local_size; ++j)
112  mass_matrix[i][j] += phi[i] * phi[j] * factor;
113  }
114 
115  mass_matrix.invert();
116  ++_leaf_index;
117 
118  }
119 
120  inverse_mass_matrix_calculator(MassMatrices& mass_matrices, const Cell& element, size_type intorder)
121  : _element(element)
122  , _mass_matrices(mass_matrices)
123  , _quadrature_rule(QuadratureRules<DF,dim>::rule(element.type(),intorder))
124  , _leaf_index(0)
125  {}
126 
127  const Cell& _element;
128  MassMatrices& _mass_matrices;
129  const QuadratureRule<DF,dim>& _quadrature_rule;
130  size_type _leaf_index;
131 
132  };
133 
134  } // anonymous namespace
135 
136 
144  template<class GFS, class U>
146  {
147  typedef typename GFS::Traits::GridViewType::Grid Grid;
148  typedef typename Grid::template Codim<0>::Entity Element;
150  typedef typename U::ElementType DF;
151 
152  public:
153 
154  typedef DynamicMatrix<typename U::ElementType> MassMatrix;
155  typedef array<MassMatrix,TypeTree::TreeInfo<GFS>::leafCount> MassMatrices;
156 
161  explicit L2Projection(const GFS& gfs, int intorder = 2)
162  : _gfs(gfs)
163  , _intorder(intorder)
164  , _inverse_mass_matrices(GlobalGeometryTypeIndex::size(Element::dimension))
165  {}
166 
172  const MassMatrices& inverseMassMatrices(const Element& e)
173  {
174  const GeometryType gt = e.geometry().type();
175  MassMatrices& inverse_mass_matrices = _inverse_mass_matrices[GlobalGeometryTypeIndex::index(gt)];
176  // if the matrix isn't empty, it has already been cached
177  if (inverse_mass_matrices[0].N() > 0)
178  return inverse_mass_matrices;
179 
180  inverse_mass_matrix_calculator<MassMatrices,Element> calculate_mass_matrices(
181  inverse_mass_matrices,
182  e,
183  _intorder
184  );
185 
186  TypeTree::applyToTree(_gfs,calculate_mass_matrices);
187 
188  return inverse_mass_matrices;
189  }
190 
191  private:
192 
193  GFS _gfs;
194  int _intorder;
195  std::vector<MassMatrices> _inverse_mass_matrices;
196  };
197 
198 
199  template<typename GFS, typename DOFVector, typename TransferMap>
201  : public TypeTree::TreeVisitor
202  , public TypeTree::DynamicTraversal
203  {
204 
208 
209  typedef typename GFS::Traits::GridView::Grid::LocalIdSet IDSet;
210  typedef typename GFS::Traits::GridView::template Codim<0>::Entity Cell;
211  typedef typename Cell::Geometry Geometry;
212  static const int dim = Geometry::mydimension;
213  typedef typename Cell::HierarchicIterator HierarchicIterator;
214  typedef typename DOFVector::ElementType RF;
215  typedef typename TransferMap::mapped_type LocalDOFVector;
216 
217 
221 
222  typedef std::size_t size_type;
223  typedef typename GFS::Traits::GridView::ctype DF;
224 
225  template<typename LFSLeaf, typename TreePath>
226  void leaf(const LFSLeaf& leaf_lfs, TreePath treePath)
227  {
228 
229  typedef typename LFSLeaf::Traits::GridFunctionSpace::Traits::FiniteElementMap FEM;
230  typedef typename FEM::Traits::FiniteElement FE;
231  const FEM& fem = leaf_lfs.gridFunctionSpace().finiteElementMap();
232  size_type fine_offset = _leaf_offset_cache[_current.type()][_leaf_index];
233  size_type coarse_offset = _leaf_offset_cache[_ancestor.type()][_leaf_index];
234 
235  typedef typename FE::Traits::LocalBasisType::Traits::RangeType Range;
236 
237  const MassMatrix& inverse_mass_matrix = _projection.inverseMassMatrices(_element)[_leaf_index];
238 
239  std::vector<Range> coarse_phi;
240  std::vector<Range> fine_phi;
241 
242  Geometry fine_geometry = _current.geometry();
243  Geometry coarse_geometry = _ancestor.geometry();
244 
245  // iterate over quadrature points
246  for (const auto& ip : QuadratureRules<DF,dim>::rule(_current.type(),_int_order))
247  {
248  typename Geometry::LocalCoordinate coarse_local = coarse_geometry.local(fine_geometry.global(ip.position()));
249  const FE* fe = &fem.find(_current);
250  fe->localBasis().evaluateFunction(ip.position(),fine_phi);
251  fe = &fem.find(_ancestor);
252  fe->localBasis().evaluateFunction(coarse_local,coarse_phi);
253  const DF factor = ip.weight()
254  * fine_geometry.integrationElement(ip.position())
255  / coarse_geometry.integrationElement(coarse_local);
256 
257  Range val(0.0);
258  for (size_type i = 0; i < fine_phi.size(); ++i)
259  {
260  val.axpy(_u_fine[fine_offset + i],fine_phi[i]);
261  }
262 
263  for (size_type i = 0; i < coarse_phi.size(); ++i)
264  {
265  Range x(0.0);
266  for (size_type j = 0; j < inverse_mass_matrix.M(); ++j)
267  x.axpy(inverse_mass_matrix[i][j],coarse_phi[j]);
268  (*_u_coarse)[coarse_offset + i] += factor * (x * val);
269  }
270  }
271 
272  ++_leaf_index;
273  }
274 
275  void operator()(const Cell& element)
276  {
277  _element = element;
278 
279  _lfs.bind(_element);
280  _lfs_cache.update();
281  _u_view.bind(_lfs_cache);
283  _u_coarse->resize(_lfs.size());
284  _u_view.read(*_u_coarse);
285  _u_view.unbind();
286 
288 
289  size_type max_level = _lfs.gridFunctionSpace().gridView().grid().maxLevel();
290 
292  while (_ancestor.mightVanish())
293  {
294  // work around UG bug!
295  if (!_ancestor.hasFather())
296  break;
297 
298  _ancestor = _ancestor.father();
299 
301  // don't project more than once
302  if (_u_coarse->size() > 0)
303  continue;
304  _u_coarse->resize(_leaf_offset_cache[_ancestor.type()].back());
305  std::fill(_u_coarse->begin(),_u_coarse->end(),RF(0));
306 
307  for (const auto& child : descendantElements(_ancestor,max_level))
308  {
309  // only evaluate on entities with data
310  if (child.isLeaf())
311  {
312  _current = child;
313  // reset leaf_index for next run over tree
314  _leaf_index = 0;
315  // load data
316  _lfs.bind(_current);
318  _lfs_cache.update();
319  _u_view.bind(_lfs_cache);
320  _u_fine.resize(_lfs_cache.size());
321  _u_view.read(_u_fine);
322  _u_view.unbind();
323  // do projection on all leafs
324  TypeTree::applyToTree(_lfs,*this);
325  }
326  }
327  }
328  }
329 
330  backup_visitor(const GFS& gfs,
331  Projection& projection,
332  const DOFVector& u,
333  LeafOffsetCache& leaf_offset_cache,
334  TransferMap& transfer_map,
335  std::size_t int_order = 2)
336  : _lfs(gfs)
337  , _lfs_cache(_lfs)
338  , _id_set(gfs.gridView().grid().localIdSet())
339  , _projection(projection)
340  , _u_view(u)
341  , _transfer_map(transfer_map)
342  , _u_coarse(nullptr)
343  , _leaf_offset_cache(leaf_offset_cache)
344  , _int_order(int_order)
345  , _leaf_index(0)
346  {}
347 
348  LFS _lfs;
349  LFSCache _lfs_cache;
350  const IDSet& _id_set;
351  Cell _element;
352  Cell _ancestor;
353  Cell _current;
354  Projection& _projection;
355  typename DOFVector::template ConstLocalView<LFSCache> _u_view;
356  TransferMap& _transfer_map;
357  LocalDOFVector* _u_coarse;
358  LeafOffsetCache& _leaf_offset_cache;
359  size_type _int_order;
360  size_type _leaf_index;
361  LocalDOFVector _u_fine;
362 
363  };
364 
365 
366 
367  template<typename GFS, typename DOFVector, typename CountVector>
369  : public TypeTree::TreeVisitor
370  , public TypeTree::DynamicTraversal
371  {
372 
376 
377  typedef typename LFS::Traits::GridFunctionSpace::Traits::GridView::template Codim<0>::Entity Cell;
378  typedef typename Cell::Geometry Geometry;
379  typedef typename DOFVector::ElementType RF;
380  typedef std::vector<RF> LocalDOFVector;
381  typedef std::vector<typename CountVector::ElementType> LocalCountVector;
382 
383  typedef std::size_t size_type;
384 
385  template<typename FiniteElement>
387  {
388 
389  template<typename X, typename Y>
390  void evaluate(const X& x, Y& y) const
391  {
392  _phi.resize(_finite_element.localBasis().size());
393  _finite_element.localBasis().evaluateFunction(_coarse_geometry.local(_fine_geometry.global(x)),_phi);
394  y = 0;
395  for (size_type i = 0; i < _phi.size(); ++i)
396  y.axpy(_dofs[_offset + i],_phi[i]);
397  }
398 
399  coarse_function(const FiniteElement& finite_element, Geometry coarse_geometry, Geometry fine_geometry, const LocalDOFVector& dofs, size_type offset)
400  : _finite_element(finite_element)
401  , _coarse_geometry(coarse_geometry)
402  , _fine_geometry(fine_geometry)
403  , _dofs(dofs)
404  , _offset(offset)
405  {}
406 
407  const FiniteElement& _finite_element;
409  Geometry _fine_geometry;
410  const LocalDOFVector& _dofs;
411  mutable std::vector<typename FiniteElement::Traits::LocalBasisType::Traits::RangeType> _phi;
412  size_type _offset;
413 
414  };
415 
416 
417  template<typename LeafLFS, typename TreePath>
418  void leaf(const LeafLFS& leaf_lfs, TreePath treePath)
419  {
420 
421  typedef typename LeafLFS::Traits::GridFunctionSpace::Traits::FiniteElementMap FEM;
422  const FEM& fem = leaf_lfs.gridFunctionSpace().finiteElementMap();
423  size_type element_offset = _leaf_offset_cache[_element.type()][_leaf_index];
424  size_type ancestor_offset = _leaf_offset_cache[_ancestor.type()][_leaf_index];
425 
426  coarse_function<typename FEM::Traits::FiniteElement> f(fem.find(_ancestor),_ancestor.geometry(),_element.geometry(),*_u_coarse,ancestor_offset);
427  const typename FEM::Traits::FiniteElement& fe = fem.find(_element);
428 
429  _u_tmp.resize(fe.localBasis().size());
430  std::fill(_u_tmp.begin(),_u_tmp.end(),RF(0.0));
431  fe.localInterpolation().interpolate(f,_u_tmp);
432  std::copy(_u_tmp.begin(),_u_tmp.end(),_u_fine.begin() + element_offset);
433 
434  ++_leaf_index;
435  }
436 
437  void operator()(const Cell& element, const Cell& ancestor, const LocalDOFVector& u_coarse)
438  {
439  _element = element;
440  _ancestor = ancestor;
441  _u_coarse = &u_coarse;
442  _lfs.bind(_element);
444  _lfs_cache.update();
445  _u_view.bind(_lfs_cache);
446 
447  // test identity using ids
448  if (_lfs.gridFunctionSpace().gridView().grid().localIdSet().id(element) ==
449  _lfs.gridFunctionSpace().gridView().grid().localIdSet().id(ancestor))
450  {
451  // no interpolation necessary, just copy the saved data
452  _u_view.add(*_u_coarse);
453  }
454  else
455  {
456  _u_fine.resize(_lfs_cache.size());
457  std::fill(_u_fine.begin(),_u_fine.end(),RF(0));
458  _leaf_index = 0;
459  TypeTree::applyToTree(_lfs,*this);
460  _u_view.add(_u_fine);
461  }
462  _u_view.commit();
463 
464  _uc_view.bind(_lfs_cache);
465  _counts.resize(_lfs_cache.size(),1);
466  _uc_view.add(_counts);
467  _uc_view.commit();
468  }
469 
470  replay_visitor(const GFS& gfs, DOFVector& u, CountVector& uc, LeafOffsetCache& leaf_offset_cache)
471  : _lfs(gfs)
472  , _lfs_cache(_lfs)
473  , _u_view(u)
474  , _uc_view(uc)
475  , _leaf_offset_cache(leaf_offset_cache)
476  , _leaf_index(0)
477  {}
478 
479  LFS _lfs;
480  LFSCache _lfs_cache;
481  Cell _element;
482  Cell _ancestor;
483  typename DOFVector::template LocalView<LFSCache> _u_view;
484  typename CountVector::template LocalView<LFSCache> _uc_view;
485  const LocalDOFVector* _u_coarse;
486  LeafOffsetCache& _leaf_offset_cache;
487  size_type _leaf_index;
488  LocalDOFVector _u_fine;
489  LocalDOFVector _u_tmp;
490  LocalCountVector _counts;
491 
492  };
493 
494 
508  template<class Grid, class GFSU, class U, class Projection>
510  {
511  typedef typename Grid::LeafGridView LeafGridView;
512  typedef typename LeafGridView::template Codim<0>
513  ::template Partition<Dune::Interior_Partition>::Iterator LeafIterator;
514  typedef typename Grid::template Codim<0>::Entity Element;
515  typedef typename Grid::LocalIdSet IDSet;
516  typedef typename IDSet::IdType ID;
517 
518  public:
519  typedef std::unordered_map<ID,std::vector<typename U::ElementType> > MapType;
520 
521 
526  explicit GridAdaptor(const GFSU& gfs)
527  : _leaf_offset_cache(gfs)
528  {}
529 
530  /* @brief @todo
531  *
532  * @param[in] u The solution that will be saved
533  * @param[out] transferMap The map containing the solution during adaptation
534  */
535  void backupData(Grid& grid, GFSU& gfsu, Projection& projection, U& u, MapType& transfer_map)
536  {
537  typedef backup_visitor<GFSU,U,MapType> Visitor;
538 
539  Visitor visitor(gfsu,projection,u,_leaf_offset_cache,transfer_map);
540 
541  // iterate over all elems
542  for(const auto& cell : elements(grid.leafGridView(),Partitions::interior))
543  visitor(cell);
544  }
545 
546  /* @brief @todo
547  *
548  * @param[out] u The solution after adaptation
549  * @param[in] transferMap The map that contains the information for the rebuild of u
550  */
551  void replayData(Grid& grid, GFSU& gfsu, Projection& projection, U& u, const MapType& transfer_map)
552  {
553  const IDSet& id_set = grid.globalIdSet();
554 
555  using CountVector = Backend::Vector<GFSU,int>;
556  CountVector uc(gfsu,0);
557 
558  typedef replay_visitor<GFSU,U,CountVector> Visitor;
559  Visitor visitor(gfsu,u,uc,_leaf_offset_cache);
560 
561  // iterate over all elems
562  LeafGridView leafView = grid.leafGridView();
563  for (const auto& cell : elements(leafView,Partitions::interior))
564  {
565  Element ancestor = cell;
566 
567  typename MapType::const_iterator map_it;
568  while ((map_it = transfer_map.find(id_set.id(ancestor))) == transfer_map.end())
569  {
570  if (!ancestor.hasFather())
571  DUNE_THROW(Exception,
572  "transferMap of GridAdaptor didn't contain ancestor of element with id " << id_set.id(ancestor));
573  ancestor = ancestor.father();
574  }
575 
576  visitor(cell,ancestor,map_it->second);
577  }
578 
579  typedef Dune::PDELab::AddDataHandle<GFSU,U> DOFHandle;
580  DOFHandle addHandle1(gfsu,u);
581  leafView.communicate (addHandle1,
582  Dune::InteriorBorder_InteriorBorder_Interface,Dune::ForwardCommunication);
584  CountHandle addHandle2(gfsu,uc);
585  leafView.communicate (addHandle2,
586  Dune::InteriorBorder_InteriorBorder_Interface,Dune::ForwardCommunication);
587 
588  // normalize multiple-interpolated DOFs by taking the arithmetic average
589  typename CountVector::iterator ucit = uc.begin();
590  for (typename U::iterator uit = u.begin(), uend = u.end(); uit != uend; ++uit, ++ucit)
591  (*uit) /= ((*ucit) > 0 ? (*ucit) : 1.0);
592  }
593 
594  private:
595 
596  LeafOffsetCache<GFSU> _leaf_offset_cache;
597 
598  };
599 
610  template<class Grid, class GFS, class X>
611  void adapt_grid (Grid& grid, GFS& gfs, X& x1, int int_order)
612  {
613  typedef L2Projection<GFS,X> Projection;
614  Projection projection(gfs,int_order);
615 
616  GridAdaptor<Grid,GFS,X,Projection> grid_adaptor(gfs);
617 
618  // prepare the grid for refinement
619  grid.preAdapt();
620 
621  // save u
622  typename GridAdaptor<Grid,GFS,X,Projection>::MapType transferMap1;
623  grid_adaptor.backupData(grid,gfs,projection,x1,transferMap1);
624 
625  // adapt the grid
626  grid.adapt();
627 
628  // update the function spaces
629  gfs.update();
630 
631  // reset u
632  x1 = X(gfs,0.0);
633  grid_adaptor.replayData(grid,gfs,projection,x1,transferMap1);
634 
635  // clean up
636  grid.postAdapt();
637  }
638 
650  template<class Grid, class GFS, class X>
651  void adapt_grid (Grid& grid, GFS& gfs, X& x1, X& x2, int int_order)
652  {
653  typedef L2Projection<GFS,X> Projection;
654  Projection projection(gfs,int_order);
655 
656  GridAdaptor<Grid,GFS,X,Projection> grid_adaptor(gfs);
657 
658  // prepare the grid for refinement
659  grid.preAdapt();
660 
661  // save solution
662  typename GridAdaptor<Grid,GFS,X,Projection>::MapType transferMap1;
663  grid_adaptor.backupData(grid,gfs,projection,x1,transferMap1);
664  typename GridAdaptor<Grid,GFS,X,Projection>::MapType transferMap2;
665  grid_adaptor.backupData(grid,gfs,projection,x2,transferMap2);
666 
667  // adapt the grid
668  grid.adapt();
669 
670  // update the function spaces
671  gfs.update();
672 
673  // interpolate solution
674  x1 = X(gfs,0.0);
675  grid_adaptor.replayData(grid,gfs,projection,x1,transferMap1);
676  x2 = X(gfs,0.0);
677  grid_adaptor.replayData(grid,gfs,projection,x2,transferMap2);
678 
679  // clean up
680  grid.postAdapt();
681  }
682 
683 
684  template<typename T>
685  void error_fraction(const T& x, typename T::ElementType alpha, typename T::ElementType beta,
686  typename T::ElementType& eta_alpha, typename T::ElementType& eta_beta, int verbose=0)
687  {
688  if (verbose>0)
689  std::cout << "+++ error fraction: alpha=" << alpha << " beta=" << beta << std::endl;
690  const int steps=20; // max number of bisection steps
691  typedef typename T::ElementType NumberType;
692  NumberType total_error = x.one_norm();
693  NumberType max_error = x.infinity_norm();
694  NumberType eta_alpha_left = 0.0;
695  NumberType eta_alpha_right = max_error;
696  NumberType eta_beta_left = 0.0;
697  NumberType eta_beta_right = max_error;
698  for (int j=1; j<=steps; j++)
699  {
700  eta_alpha = 0.5*(eta_alpha_left+eta_alpha_right);
701  eta_beta = 0.5*(eta_beta_left+eta_beta_right);
702  NumberType sum_alpha=0.0;
703  NumberType sum_beta=0.0;
704  unsigned int alpha_count = 0;
705  unsigned int beta_count = 0;
706  for (typename T::const_iterator it = x.begin(),
707  end = x.end();
708  it != end;
709  ++it)
710  {
711  if (*it >=eta_alpha) { sum_alpha += *it; alpha_count++;}
712  if (*it < eta_beta) { sum_beta += *it; beta_count++;}
713  }
714  if (verbose>1)
715  {
716  std::cout << "+++ " << j << " eta_alpha=" << eta_alpha << " alpha_fraction=" << sum_alpha/total_error
717  << " elements: " << alpha_count << " of " << x.N() << std::endl;
718  std::cout << "+++ " << j << " eta_beta=" << eta_beta << " beta_fraction=" << sum_beta/total_error
719  << " elements: " << beta_count << " of " << x.N() << std::endl;
720  }
721  if (std::abs(alpha-sum_alpha/total_error) <= 0.01 && std::abs(beta-sum_beta/total_error) <= 0.01) break;
722  if (sum_alpha>alpha*total_error)
723  eta_alpha_left = eta_alpha;
724  else
725  eta_alpha_right = eta_alpha;
726  if (sum_beta>beta*total_error)
727  eta_beta_right = eta_beta;
728  else
729  eta_beta_left = eta_beta;
730  }
731  if (verbose>0)
732  {
733  std::cout << "+++ refine_threshold=" << eta_alpha
734  << " coarsen_threshold=" << eta_beta << std::endl;
735  }
736  }
737 
738 
739  template<typename T>
740  void element_fraction(const T& x, typename T::ElementType alpha, typename T::ElementType beta,
741  typename T::ElementType& eta_alpha, typename T::ElementType& eta_beta, int verbose=0)
742  {
743  const int steps=20; // max number of bisection steps
744  typedef typename T::ElementType NumberType;
745  NumberType total_error =x.N();
746  NumberType max_error = x.infinity_norm();
747  NumberType eta_alpha_left = 0.0;
748  NumberType eta_alpha_right = max_error;
749  NumberType eta_beta_left = 0.0;
750  NumberType eta_beta_right = max_error;
751  for (int j=1; j<=steps; j++)
752  {
753  eta_alpha = 0.5*(eta_alpha_left+eta_alpha_right);
754  eta_beta = 0.5*(eta_beta_left+eta_beta_right);
755  NumberType sum_alpha=0.0;
756  NumberType sum_beta=0.0;
757  unsigned int alpha_count = 0;
758  unsigned int beta_count = 0;
759 
760  for (typename T::const_iterator it = x.begin(),
761  end = x.end();
762  it != end;
763  ++it)
764  {
765  if (*it>=eta_alpha) { sum_alpha += 1.0; alpha_count++;}
766  if (*it< eta_beta) { sum_beta +=1.0; beta_count++;}
767  }
768  if (verbose>1)
769  {
770  std::cout << j << " eta_alpha=" << eta_alpha << " alpha_fraction=" << sum_alpha/total_error
771  << " elements: " << alpha_count << " of " << x.N() << std::endl;
772  std::cout << j << " eta_beta=" << eta_beta << " beta_fraction=" << sum_beta/total_error
773  << " elements: " << beta_count << " of " << x.N() << std::endl;
774  }
775  if (std::abs(alpha-sum_alpha/total_error) <= 0.01 && std::abs(beta-sum_beta/total_error) <= 0.01) break;
776  if (sum_alpha>alpha*total_error)
777  eta_alpha_left = eta_alpha;
778  else
779  eta_alpha_right = eta_alpha;
780  if (sum_beta>beta*total_error)
781  eta_beta_right = eta_beta;
782  else
783  eta_beta_left = eta_beta;
784  }
785  if (verbose>0)
786  {
787  std::cout << "+++ refine_threshold=" << eta_alpha
788  << " coarsen_threshold=" << eta_beta << std::endl;
789  }
790  }
791 
794  template<typename T>
795  void error_distribution(const T& x, int bins)
796  {
797  const int steps=30; // max number of bisection steps
798  typedef typename T::ElementType NumberType;
799  NumberType total_error = x.one_norm();
800  NumberType total_elements = x.N();
801  NumberType max_error = x.infinity_norm();
802  std::vector<NumberType> left(bins,0.0);
803  std::vector<NumberType> right(bins,max_error*(1.0+1e-8));
804  std::vector<NumberType> eta(bins);
805  std::vector<NumberType> target(bins);
806  for (unsigned int k=0; k<bins; k++)
807  target[k]= (k+1)/((NumberType)bins);
808  for (int j=1; j<=steps; j++)
809  {
810  for (unsigned int k=0; k<bins; k++)
811  eta[k]= 0.5*(left[k]+right[k]);
812  std::vector<NumberType> sum(bins,0.0);
813  std::vector<int> count(bins,0);
814 
815  for (typename T::const_iterator it = x.begin(),
816  end = x.end();
817  it != end;
818  ++it)
819  {
820  for (unsigned int k=0; k<bins; k++)
821  if (*it<=eta[k])
822  {
823  sum[k] += *it;
824  count[k] += 1;
825  }
826  }
827  // std::cout << std::endl;
828  // std::cout << "// step " << j << std::endl;
829  // for (unsigned int k=0; k<bins; k++)
830  // std::cout << k+1 << " " << count[k] << " " << eta[k] << " " << right[k]-left[k]
831  // << " " << sum[k]/total_error << " " << target[k] << std::endl;
832  for (unsigned int k=0; k<bins; k++)
833  if (sum[k]<=target[k]*total_error)
834  left[k] = eta[k];
835  else
836  right[k] = eta[k];
837  }
838  std::vector<NumberType> sum(bins,0.0);
839  std::vector<int> count(bins,0);
840  for (unsigned int i=0; i<x.N(); i++)
841  for (unsigned int k=0; k<bins; k++)
842  if (x[i]<=eta[k])
843  {
844  sum[k] += x[i];
845  count[k] += 1;
846  }
847  std::cout << "+++ error distribution" << std::endl;
848  std::cout << "+++ number of elements: " << x.N() << std::endl;
849  std::cout << "+++ max element error: " << max_error << std::endl;
850  std::cout << "+++ total error: " << total_error << std::endl;
851  std::cout << "+++ bin #elements eta sum/total " << std::endl;
852  for (unsigned int k=0; k<bins; k++)
853  std::cout << "+++ " << k+1 << " " << count[k] << " " << eta[k] << " " << sum[k]/total_error << std::endl;
854  }
855 
856  template<typename Grid, typename X>
857  void mark_grid (Grid &grid, const X& x, typename X::ElementType refine_threshold,
858  typename X::ElementType coarsen_threshold, int min_level = 0, int max_level = std::numeric_limits<int>::max(), int verbose=0)
859  {
860  typedef typename Grid::LeafGridView GV;
861 
862  GV gv = grid.leafGridView();
863 
864  unsigned int refine_cnt=0;
865  unsigned int coarsen_cnt=0;
866 
867  typedef typename X::GridFunctionSpace GFS;
868  typedef LocalFunctionSpace<GFS> LFS;
869  typedef LFSIndexCache<LFS> LFSCache;
870  typedef typename X::template ConstLocalView<LFSCache> XView;
871 
872  LFS lfs(x.gridFunctionSpace());
873  LFSCache lfs_cache(lfs);
874  XView x_view(x);
875 
876  for(const auto& cell : elements(gv))
877  {
878  lfs.bind(cell);
879  lfs_cache.update();
880  x_view.bind(lfs_cache);
881 
882  if (x_view[0]>=refine_threshold && cell.level() < max_level)
883  {
884  grid.mark(1,cell);
885  refine_cnt++;
886  }
887  if (x_view[0]<=coarsen_threshold && cell.level() > min_level)
888  {
889  grid.mark(-1,cell);
890  coarsen_cnt++;
891  }
892  x_view.unbind();
893  }
894  if (verbose>0)
895  std::cout << "+++ mark_grid: " << refine_cnt << " marked for refinement, "
896  << coarsen_cnt << " marked for coarsening" << std::endl;
897  }
898 
899 
900  template<typename Grid, typename X>
901  void mark_grid_for_coarsening (Grid &grid, const X& x, typename X::ElementType refine_threshold,
902  typename X::ElementType coarsen_threshold, int verbose=0)
903  {
904  typedef typename Grid::LeafGridView GV;
905 
906  GV gv = grid.leafGridView();
907 
908  unsigned int coarsen_cnt=0;
909 
910  typedef typename X::GridFunctionSpace GFS;
911  typedef LocalFunctionSpace<GFS> LFS;
912  typedef LFSIndexCache<LFS> LFSCache;
913  typedef typename X::template ConstLocalView<LFSCache> XView;
914 
915  LFS lfs(x.gridFunctionSpace());
916  LFSCache lfs_cache(lfs);
917  XView x_view(x);
918 
919  for(const auto& cell : elements(gv))
920  {
921  lfs.bind(cell);
922  lfs_cache.update();
923  x_view.bind(lfs_cache);
924 
925  if (x_view[0]>=refine_threshold)
926  {
927  grid.mark(-1,cell);
928  coarsen_cnt++;
929  }
930  if (x_view[0]<=coarsen_threshold)
931  {
932  grid.mark(-1,cell);
933  coarsen_cnt++;
934  }
935  }
936  if (verbose>0)
937  std::cout << "+++ mark_grid_for_coarsening: "
938  << coarsen_cnt << " marked for coarsening" << std::endl;
939  }
940 
941 
943  {
944  // strategy parameters
945  double scaling;
946  double optimistic_factor;
947  double coarsen_limit;
948  double balance_limit;
949  double tol;
950  double T;
951  int verbose;
952  bool no_adapt;
953  double refine_fraction_while_refinement;
954  double coarsen_fraction_while_refinement;
955  double coarsen_fraction_while_coarsening;
956  double timestep_decrease_factor;
957  double timestep_increase_factor;
958 
959  // results to be reported to the user after evaluating the error
960  bool accept;
961  bool adapt_dt;
962  bool adapt_grid;
963  double newdt;
964  double q_s, q_t;
965 
966  // state variables
967  bool have_decreased_time_step;
968  bool have_refined_grid;
969 
970  // the only state variable: accumulated error
971  double accumulated_estimated_error_squared;
972  double minenergy_rate;
973 
974  public:
975  TimeAdaptationStrategy (double tol_, double T_, int verbose_=0)
976  : scaling(16.0), optimistic_factor(1.0), coarsen_limit(0.5), balance_limit(0.33333),
977  tol(tol_), T(T_), verbose(verbose_), no_adapt(false),
978  refine_fraction_while_refinement(0.7),
979  coarsen_fraction_while_refinement(0.2),
980  coarsen_fraction_while_coarsening(0.2),
981  timestep_decrease_factor(0.5), timestep_increase_factor(1.5),
982  accept(false), adapt_dt(false), adapt_grid(false), newdt(1.0),
983  have_decreased_time_step(false), have_refined_grid(false),
984  accumulated_estimated_error_squared(0.0),
985  minenergy_rate(0.0)
986  {
987  }
988 
990  {
991  timestep_decrease_factor=s;
992  }
993 
995  {
996  timestep_increase_factor=s;
997  }
998 
1000  {
1001  refine_fraction_while_refinement=s;
1002  }
1003 
1005  {
1006  coarsen_fraction_while_refinement=s;
1007  }
1008 
1010  {
1011  coarsen_fraction_while_coarsening=s;
1012  }
1013 
1014  void setMinEnergyRate (double s)
1015  {
1016  minenergy_rate=s;
1017  }
1018 
1019  void setCoarsenLimit (double s)
1020  {
1021  coarsen_limit=s;
1022  }
1023 
1024  void setBalanceLimit (double s)
1025  {
1026  balance_limit=s;
1027  }
1028 
1029  void setTemporalScaling (double s)
1030  {
1031  scaling=s;
1032  }
1033 
1034  void setOptimisticFactor (double s)
1035  {
1036  optimistic_factor=s;
1037  }
1038 
1040  {
1041  no_adapt = false;
1042  }
1043 
1045  {
1046  no_adapt = true;
1047  }
1048 
1049  bool acceptTimeStep () const
1050  {
1051  return accept;
1052  }
1053 
1054  bool adaptDT () const
1055  {
1056  return adapt_dt;
1057  }
1058 
1059  bool adaptGrid () const
1060  {
1061  return adapt_grid;
1062  }
1063 
1064  double newDT () const
1065  {
1066  return newdt;
1067  }
1068 
1069  double qs () const
1070  {
1071  return q_s;
1072  }
1073 
1074  double qt () const
1075  {
1076  return q_t;
1077  }
1078 
1079  double endT() const
1080  {
1081  return T;
1082  }
1083 
1084  double accumulatedErrorSquared () const
1085  {
1086  return accumulated_estimated_error_squared;
1087  }
1088 
1089  // to be called when new time step is done
1091  {
1092  have_decreased_time_step = false;
1093  have_refined_grid = false;
1094  }
1095 
1096  template<typename GM, typename X>
1097  void evaluate_estimators (GM& grid, double time, double dt, const X& eta_space, const X& eta_time, double energy_timeslab)
1098  {
1099  accept=false;
1100  adapt_dt=false;
1101  adapt_grid=false;
1102  newdt=dt;
1103 
1104  double spatial_error = eta_space.one_norm();
1105  double temporal_error = scaling*eta_time.one_norm();
1106  double sum = spatial_error + temporal_error;
1107  //double allowed = optimistic_factor*(tol*tol-accumulated_estimated_error_squared)*dt/(T-time);
1108  double allowed = tol*tol*(energy_timeslab+minenergy_rate*dt);
1109  q_s = spatial_error/sum;
1110  q_t = temporal_error/sum;
1111 
1112  // print some statistics
1113  if (verbose>0)
1114  std::cout << "+++"
1115  << " q_s=" << q_s
1116  << " q_t=" << q_t
1117  << " sum/allowed=" << sum/allowed
1118  // << " allowed=" << allowed
1119  << " estimated error=" << sqrt(accumulated_estimated_error_squared+sum)
1120  << " energy_rate=" << energy_timeslab/dt
1121  << std::endl;
1122 
1123  // for simplicity: a mode that does no adaptation at all
1124  if (no_adapt)
1125  {
1126  accept = true;
1127  accumulated_estimated_error_squared += sum;
1128  if (verbose>1) std::cout << "+++ no adapt mode" << std::endl;
1129  return;
1130  }
1131 
1132  // the adaptation strategy
1133  if (sum<=allowed)
1134  {
1135  // we will accept this time step
1136  accept = true;
1137  if (verbose>1) std::cout << "+++ accepting time step" << std::endl;
1138  accumulated_estimated_error_squared += sum;
1139 
1140  // check if grid size or time step needs to be adapted
1141  if (sum<coarsen_limit*allowed)
1142  {
1143  // the error is too small, i.e. the computation is inefficient
1144  if (q_t<balance_limit)
1145  {
1146  // spatial error is dominating => increase time step
1147  newdt = timestep_increase_factor*dt;
1148  adapt_dt = true;
1149  if (verbose>1) std::cout << "+++ spatial error dominates: increase time step" << std::endl;
1150  }
1151  else
1152  {
1153  if (q_s>balance_limit)
1154  {
1155  // step sizes balanced: coarsen in time
1156  newdt = timestep_increase_factor*dt;
1157  adapt_dt = true;
1158  if (verbose>1) std::cout << "+++ increasing time step" << std::endl;
1159  }
1160  // coarsen grid in space
1161  double eta_refine, eta_coarsen;
1162  if (verbose>1) std::cout << "+++ mark grid for coarsening" << std::endl;
1163  //error_distribution(eta_space,20);
1164  Dune::PDELab::error_fraction(eta_space,coarsen_fraction_while_coarsening,
1165  coarsen_fraction_while_coarsening,eta_refine,eta_coarsen);
1166  Dune::PDELab::mark_grid_for_coarsening(grid,eta_space,eta_refine,eta_coarsen,verbose);
1167  adapt_grid = true;
1168  }
1169  }
1170  else
1171  {
1172  // modify at least the time step
1173  if (q_t<balance_limit)
1174  {
1175  // spatial error is dominating => increase time step
1176  newdt = timestep_increase_factor*dt;
1177  adapt_dt = true;
1178  if (verbose>1) std::cout << "+++ spatial error dominates: increase time step" << std::endl;
1179  }
1180  }
1181  }
1182  else
1183  {
1184  // error is too large, we need to do something
1185  if (verbose>1) std::cout << "+++ will redo time step" << std::endl;
1186  if (q_t>1-balance_limit)
1187  {
1188  // temporal error is dominating => decrease time step only
1189  newdt = timestep_decrease_factor*dt;
1190  adapt_dt = true;
1191  have_decreased_time_step = true;
1192  if (verbose>1) std::cout << "+++ decreasing time step only" << std::endl;
1193  }
1194  else
1195  {
1196  if (q_t<balance_limit)
1197  {
1198  if (!have_decreased_time_step)
1199  {
1200  // time steps size not balanced (too small)
1201  newdt = timestep_increase_factor*dt;
1202  adapt_dt = true;
1203  if (verbose>1) std::cout << "+++ increasing time step" << std::endl;
1204  }
1205  }
1206  else
1207  {
1208  // step sizes balanced: refine in time as well
1209  newdt = timestep_decrease_factor*dt;
1210  adapt_dt = true;
1211  have_decreased_time_step = true;
1212  if (verbose>1) std::cout << "+++ decreasing time step" << std::endl;
1213  }
1214  // refine grid in space
1215  double eta_refine, eta_coarsen;
1216  if (verbose>1) std::cout << "+++ BINGO mark grid for refinement and coarsening" << std::endl;
1217  //error_distribution(eta_space,20);
1218  Dune::PDELab::error_fraction(eta_space,refine_fraction_while_refinement,
1219  coarsen_fraction_while_refinement,eta_refine,eta_coarsen,0);
1220  Dune::PDELab::mark_grid(grid,eta_space,eta_refine,eta_coarsen,verbose);
1221  adapt_grid = true;
1222  }
1223  }
1224  }
1225  };
1226 
1227 
1228 
1229  } // namespace PDELab
1230 } // namespace Dune
1231 
1232 #endif
void adapt_grid(Grid &grid, GFS &gfs, X &x1, int int_order)
adapt a grid, corresponding function space and solution vectors
Definition: adaptivity.hh:611
const std::size_t offset
Definition: localfunctionspace.hh:74
LocalFunctionSpace< GFS > LFS
Definition: adaptivity.hh:373
GFS::Traits::GridView::Grid::LocalIdSet IDSet
Definition: adaptivity.hh:209
void leaf(const LFSLeaf &leaf_lfs, TreePath treePath)
Definition: adaptivity.hh:226
const MassMatrices & inverseMassMatrices(const Element &e)
Calculate the inverse local mass matrix, used in the local L2 projection.
Definition: adaptivity.hh:172
void mark_grid(Grid &grid, const X &x, typename X::ElementType refine_threshold, typename X::ElementType coarsen_threshold, int min_level=0, int max_level=std::numeric_limits< int >::max(), int verbose=0)
Definition: adaptivity.hh:857
void setAdaptationOff()
Definition: adaptivity.hh:1044
array< MassMatrix, TypeTree::TreeInfo< GFS >::leafCount > MassMatrices
Definition: adaptivity.hh:155
LeafOffsetCache & _leaf_offset_cache
Definition: adaptivity.hh:358
Definition: adaptivity.hh:200
LFSIndexCache< LFS > LFSCache
Definition: adaptivity.hh:374
DOFVector::template LocalView< LFSCache > _u_view
Definition: adaptivity.hh:483
const LeafOffsets & operator[](GeometryType gt) const
Definition: adaptivity.hh:42
Definition: adaptivity.hh:368
Definition: adaptivity.hh:942
void evaluate_estimators(GM &grid, double time, double dt, const X &eta_space, const X &eta_time, double energy_timeslab)
Definition: adaptivity.hh:1097
std::size_t size_type
Definition: adaptivity.hh:383
Cell::Geometry Geometry
Definition: adaptivity.hh:378
const LocalDOFVector & _dofs
Definition: adaptivity.hh:410
const Cell & _element
Definition: adaptivity.hh:127
LFSCache _lfs_cache
Definition: adaptivity.hh:480
Definition: adaptivity.hh:32
size_type _offset
Definition: adaptivity.hh:412
void evaluate(const X &x, Y &y) const
Definition: adaptivity.hh:390
bool acceptTimeStep() const
Definition: adaptivity.hh:1049
Projection & _projection
Definition: adaptivity.hh:354
typename impl::BackendVectorSelector< GridFunctionSpace, FieldType >::Type Vector
alias of the return type of BackendVectorSelector
Definition: backend/interface.hh:112
const E & e
Definition: interpolate.hh:172
size_type _int_order
Definition: adaptivity.hh:359
LocalCountVector _counts
Definition: adaptivity.hh:490
void backupData(Grid &grid, GFSU &gfsu, Projection &projection, U &u, MapType &transfer_map)
Definition: adaptivity.hh:535
Projection::MassMatrices MassMatrices
Definition: adaptivity.hh:219
std::size_t size_type
Definition: adaptivity.hh:222
void setOptimisticFactor(double s)
Definition: adaptivity.hh:1034
LocalDOFVector _u_fine
Definition: adaptivity.hh:488
size_type size() const
Definition: lfsindexcache.hh:423
Dune::PDELab::LeafOffsetCache< GFS > LeafOffsetCache
Definition: adaptivity.hh:207
static const int dim
Definition: adaptivity.hh:212
Cell::HierarchicIterator HierarchicIterator
Definition: adaptivity.hh:213
void setCoarsenLimit(double s)
Definition: adaptivity.hh:1019
Cell _ancestor
Definition: adaptivity.hh:482
void setRefineFractionWhileRefinement(double s)
Definition: adaptivity.hh:999
Geometry _fine_geometry
Definition: adaptivity.hh:409
GridAdaptor(const GFSU &gfs)
The constructor.
Definition: adaptivity.hh:526
Iterator extract_lfs_leaf_sizes(const LFS &lfs, Iterator it)
Definition: lfsindexcache.hh:165
LocalDOFVector _u_tmp
Definition: adaptivity.hh:489
LocalFunctionSpace< GFS > LFS
Definition: adaptivity.hh:36
std::vector< LeafOffsets > _leaf_offset_cache
Definition: adaptivity.hh:70
Definition: adaptivity.hh:145
std::vector< typename FiniteElement::Traits::LocalBasisType::Traits::RangeType > _phi
Definition: adaptivity.hh:411
Dune::PDELab::LeafOffsetCache< GFS > LeafOffsetCache
Definition: adaptivity.hh:375
Traits::IndexContainer::size_type size() const
get current size
Definition: localfunctionspace.hh:206
void mark_grid_for_coarsening(Grid &grid, const X &x, typename X::ElementType refine_threshold, typename X::ElementType coarsen_threshold, int verbose=0)
Definition: adaptivity.hh:901
L2Projection(const GFS &gfs, int intorder=2)
The constructor.
Definition: adaptivity.hh:161
CountVector::template LocalView< LFSCache > _uc_view
Definition: adaptivity.hh:484
MassMatrices & _mass_matrices
Definition: adaptivity.hh:128
size_type _leaf_index
Definition: adaptivity.hh:130
DOFVector::template ConstLocalView< LFSCache > _u_view
Definition: adaptivity.hh:355
const FiniteElement & _finite_element
Definition: adaptivity.hh:407
LFS _lfs
Definition: adaptivity.hh:69
void setCoarsenFractionWhileRefinement(double s)
Definition: adaptivity.hh:1004
LFS::Traits::GridFunctionSpace::Traits::GridView::template Codim< 0 >::Entity Cell
Definition: adaptivity.hh:377
const LocalDOFVector * _u_coarse
Definition: adaptivity.hh:485
L2Projection< typename LFS::Traits::GridFunctionSpace, DOFVector > Projection
Definition: adaptivity.hh:218
TransferMap::mapped_type LocalDOFVector
Definition: adaptivity.hh:215
LFSIndexCache< LFS > LFSCache
Definition: adaptivity.hh:206
replay_visitor(const GFS &gfs, DOFVector &u, CountVector &uc, LeafOffsetCache &leaf_offset_cache)
Definition: adaptivity.hh:470
std::vector< RF > LocalDOFVector
Definition: adaptivity.hh:380
void update()
Definition: lfsindexcache.hh:300
LFS _lfs
Definition: adaptivity.hh:479
Base class for all PDELab exceptions.
Definition: exceptions.hh:17
void setTimeStepIncreaseFactor(double s)
Definition: adaptivity.hh:994
std::vector< typename CountVector::ElementType > LocalCountVector
Definition: adaptivity.hh:381
Definition: genericdatahandle.hh:622
LFS _lfs
Definition: adaptivity.hh:348
Projection::MassMatrix MassMatrix
Definition: adaptivity.hh:220
void setAdaptationOn()
Definition: adaptivity.hh:1039
std::unordered_map< ID, std::vector< typename U::ElementType > > MapType
Definition: adaptivity.hh:519
Definition: adaptivity.hh:27
void update(const Cell &e)
Definition: adaptivity.hh:50
void replayData(Grid &grid, GFSU &gfsu, Projection &projection, U &u, const MapType &transfer_map)
Definition: adaptivity.hh:551
GFS::Traits::GridView::template Codim< 0 >::Entity Cell
Definition: adaptivity.hh:35
LocalFunctionSpace< GFS > LFS
Definition: adaptivity.hh:205
Cell _element
Definition: adaptivity.hh:481
void setCoarsenFractionWhileCoarsening(double s)
Definition: adaptivity.hh:1009
LocalDOFVector _u_fine
Definition: adaptivity.hh:361
void element_fraction(const T &x, typename T::ElementType alpha, typename T::ElementType beta, typename T::ElementType &eta_alpha, typename T::ElementType &eta_beta, int verbose=0)
Definition: adaptivity.hh:740
const GFS & gridFunctionSpace() const
Returns the GridFunctionSpace underlying this LocalFunctionSpace.
Definition: localfunctionspace.hh:264
void error_fraction(const T &x, typename T::ElementType alpha, typename T::ElementType beta, typename T::ElementType &eta_alpha, typename T::ElementType &eta_beta, int verbose=0)
Definition: adaptivity.hh:685
void setBalanceLimit(double s)
Definition: adaptivity.hh:1024
static const int dim
Definition: adaptivity.hh:83
void leaf(const LeafLFS &leaf_lfs, TreePath treePath)
Definition: adaptivity.hh:418
double endT() const
Definition: adaptivity.hh:1079
backup_visitor(const GFS &gfs, Projection &projection, const DOFVector &u, LeafOffsetCache &leaf_offset_cache, TransferMap &transfer_map, std::size_t int_order=2)
Definition: adaptivity.hh:330
Cell _current
Definition: adaptivity.hh:353
double newDT() const
Definition: adaptivity.hh:1064
array< std::size_t, TypeTree::TreeInfo< GFS >::leafCount+1 > LeafOffsets
Definition: adaptivity.hh:40
void startTimeStep()
Definition: adaptivity.hh:1090
TransferMap & _transfer_map
Definition: adaptivity.hh:356
Class for automatic adaptation of the grid.
Definition: adaptivity.hh:509
LeafOffsetCache & _leaf_offset_cache
Definition: adaptivity.hh:486
bool adaptGrid() const
Definition: adaptivity.hh:1059
Cell::Geometry Geometry
Definition: adaptivity.hh:211
TimeAdaptationStrategy(double tol_, double T_, int verbose_=0)
Definition: adaptivity.hh:975
double accumulatedErrorSquared() const
Definition: adaptivity.hh:1084
Cell _ancestor
Definition: adaptivity.hh:352
const std::string s
Definition: function.hh:1102
bool adaptDT() const
Definition: adaptivity.hh:1054
void setTemporalScaling(double s)
Definition: adaptivity.hh:1029
void setMinEnergyRate(double s)
Definition: adaptivity.hh:1014
void error_distribution(const T &x, int bins)
Definition: adaptivity.hh:795
LeafOffsetCache(const GFS &gfs)
Definition: adaptivity.hh:64
LFSCache _lfs_cache
Definition: adaptivity.hh:349
const IDSet & _id_set
Definition: adaptivity.hh:350
void operator()(const Cell &element, const Cell &ancestor, const LocalDOFVector &u_coarse)
Definition: adaptivity.hh:437
LocalDOFVector * _u_coarse
Definition: adaptivity.hh:357
Cell _element
Definition: adaptivity.hh:351
size_type _leaf_index
Definition: adaptivity.hh:487
double qs() const
Definition: adaptivity.hh:1069
double qt() const
Definition: adaptivity.hh:1074
DynamicMatrix< typename U::ElementType > MassMatrix
Definition: adaptivity.hh:154
void operator()(const Cell &element)
Definition: adaptivity.hh:275
GFS::Traits::GridView::template Codim< 0 >::Entity Cell
Definition: adaptivity.hh:210
void setTimeStepDecreaseFactor(double s)
Definition: adaptivity.hh:989
coarse_function(const FiniteElement &finite_element, Geometry coarse_geometry, Geometry fine_geometry, const LocalDOFVector &dofs, size_type offset)
Definition: adaptivity.hh:399
size_type _leaf_index
Definition: adaptivity.hh:360
GFS::Traits::GridView::ctype DF
Definition: adaptivity.hh:223
DOFVector::ElementType RF
Definition: adaptivity.hh:214
const QuadratureRule< DF, dim > & _quadrature_rule
Definition: adaptivity.hh:129
Geometry _coarse_geometry
Definition: adaptivity.hh:408
DOFVector::ElementType RF
Definition: adaptivity.hh:379