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