dune-pdelab  2.5-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_ADAPTIVITY_HH
5 #define DUNE_PDELAB_ADAPTIVITY_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  auto& fem = gfs.finiteElementMap();
93  auto& fe = fem.find(_element);
94  size_type local_size = fe.localBasis().size();
95 
96  MassMatrix& mass_matrix = _mass_matrices[_leaf_index];
97  mass_matrix.resize(local_size,local_size);
98 
99  using Range = typename GFS::Traits::FiniteElementMap::Traits::
100  FiniteElement::Traits::LocalBasisType::Traits::RangeType;
101  std::vector<Range> phi;
102  phi.resize(std::max(phi.size(),local_size));
103 
104  for (const auto& ip : _quadrature_rule)
105  {
106  fe.localBasis().evaluateFunction(ip.position(),phi);
107  const DF factor = ip.weight();
108 
109  for (size_type i = 0; i < local_size; ++i)
110  for (size_type j = 0; j < local_size; ++j)
111  mass_matrix[i][j] += phi[i] * phi[j] * factor;
112  }
113 
114  mass_matrix.invert();
115  ++_leaf_index;
116 
117  }
118 
119  inverse_mass_matrix_calculator(MassMatrices& mass_matrices, const Cell& element, size_type intorder)
120  : _element(element)
121  , _mass_matrices(mass_matrices)
122  , _quadrature_rule(QuadratureRules<DF,dim>::rule(element.type(),intorder))
123  , _leaf_index(0)
124  {}
125 
126  const Cell& _element;
127  MassMatrices& _mass_matrices;
128  const QuadratureRule<DF,dim>& _quadrature_rule;
129  size_type _leaf_index;
130 
131  };
132 
133  } // anonymous namespace
134 
135 
143  template<class GFS, class U>
145  {
146  using EntitySet = typename GFS::Traits::EntitySet;
147  using Element = typename EntitySet::Element;
149  typedef typename U::ElementType DF;
150 
151  public:
152 
153  typedef DynamicMatrix<typename U::ElementType> MassMatrix;
154  typedef std::array<MassMatrix,TypeTree::TreeInfo<GFS>::leafCount> MassMatrices;
155 
160  explicit L2Projection(const GFS& gfs, int intorder = 2)
161  : _gfs(gfs)
162  , _intorder(intorder)
163  , _inverse_mass_matrices(GlobalGeometryTypeIndex::size(Element::dimension))
164  {}
165 
171  const MassMatrices& inverseMassMatrices(const Element& e)
172  {
173  auto gt = e.geometry().type();
174  auto& inverse_mass_matrices = _inverse_mass_matrices[GlobalGeometryTypeIndex::index(gt)];
175  // if the matrix isn't empty, it has already been cached
176  if (inverse_mass_matrices[0].N() > 0)
177  return inverse_mass_matrices;
178 
179  inverse_mass_matrix_calculator<MassMatrices,Element> calculate_mass_matrices(
180  inverse_mass_matrices,
181  e,
182  _intorder
183  );
184 
185  TypeTree::applyToTree(_gfs,calculate_mass_matrices);
186 
187  return inverse_mass_matrices;
188  }
189 
190  private:
191 
192  GFS _gfs;
193  int _intorder;
194  std::vector<MassMatrices> _inverse_mass_matrices;
195  };
196 
197 
198  template<typename GFS, typename DOFVector, typename TransferMap>
200  : public TypeTree::TreeVisitor
201  , public TypeTree::DynamicTraversal
202  {
203 
207 
208  using EntitySet = typename GFS::Traits::EntitySet;
209  using IDSet = typename EntitySet::Traits::GridView::Grid::LocalIdSet;
210  using Element = typename EntitySet::Element;
211  typedef typename Element::Geometry Geometry;
212  static const int dim = Geometry::mydimension;
213  typedef typename DOFVector::ElementType RF;
214  typedef typename TransferMap::mapped_type LocalDOFVector;
215 
216 
220 
221  typedef std::size_t size_type;
222  using DF = typename EntitySet::Traits::CoordinateField;
223 
224  template<typename LFSLeaf, typename TreePath>
225  void leaf(const LFSLeaf& leaf_lfs, TreePath treePath)
226  {
227 
228  auto& fem = leaf_lfs.gridFunctionSpace().finiteElementMap();
229  auto fine_offset = _leaf_offset_cache[_current.type()][_leaf_index];
230  auto coarse_offset = _leaf_offset_cache[_ancestor.type()][_leaf_index];
231 
232  using Range = typename LFSLeaf::Traits::GridFunctionSpace::Traits::FiniteElementMap::
233  Traits::FiniteElement::Traits::LocalBasisType::Traits::RangeType;
234 
235  auto& inverse_mass_matrix = _projection.inverseMassMatrices(_element)[_leaf_index];
236 
237  auto coarse_phi = std::vector<Range>{};
238  auto fine_phi = std::vector<Range>{};
239 
240  auto fine_geometry = _current.geometry();
241  auto coarse_geometry = _ancestor.geometry();
242 
243  // iterate over quadrature points
244  for (const auto& ip : QuadratureRules<DF,dim>::rule(_current.type(),_int_order))
245  {
246  auto coarse_local = coarse_geometry.local(fine_geometry.global(ip.position()));
247  auto fe = &fem.find(_current);
248  fe->localBasis().evaluateFunction(ip.position(),fine_phi);
249  fe = &fem.find(_ancestor);
250  fe->localBasis().evaluateFunction(coarse_local,coarse_phi);
251  const DF factor = ip.weight()
252  * fine_geometry.integrationElement(ip.position())
253  / coarse_geometry.integrationElement(coarse_local);
254 
255  auto val = Range{0.0};
256  for (size_type i = 0; i < fine_phi.size(); ++i)
257  {
258  val.axpy(_u_fine[fine_offset + i],fine_phi[i]);
259  }
260 
261  for (size_type i = 0; i < coarse_phi.size(); ++i)
262  {
263  auto x = Range{0.0};
264  for (size_type j = 0; j < inverse_mass_matrix.M(); ++j)
265  x.axpy(inverse_mass_matrix[i][j],coarse_phi[j]);
266  (*_u_coarse)[coarse_offset + i] += factor * (x * val);
267  }
268  }
269 
270  ++_leaf_index;
271  }
272 
273  void operator()(const Element& element)
274  {
275  _element = element;
276 
277  _lfs.bind(_element);
278  _lfs_cache.update();
279  _u_view.bind(_lfs_cache);
280  _u_coarse = &_transfer_map[_id_set.id(_element)];
281  _u_coarse->resize(_lfs.size());
282  _u_view.read(*_u_coarse);
283  _u_view.unbind();
284 
286 
287  size_type max_level = _lfs.gridFunctionSpace().gridView().grid().maxLevel();
288 
289  _ancestor = _element;
290  while (_ancestor.mightVanish())
291  {
292  // work around UG bug!
293  if (!_ancestor.hasFather())
294  break;
295 
296  _ancestor = _ancestor.father();
297 
298  _u_coarse = &_transfer_map[_id_set.id(_ancestor)];
299  // don't project more than once
300  if (_u_coarse->size() > 0)
301  continue;
302  _u_coarse->resize(_leaf_offset_cache[_ancestor.type()].back());
303  std::fill(_u_coarse->begin(),_u_coarse->end(),RF(0));
304 
305  for (const auto& child : descendantElements(_ancestor,max_level))
306  {
307  // only evaluate on entities with data
308  if (child.isLeaf())
309  {
310  _current = child;
311  // reset leaf_index for next run over tree
312  _leaf_index = 0;
313  // load data
314  _lfs.bind(_current);
315  _leaf_offset_cache.update(_current);
316  _lfs_cache.update();
317  _u_view.bind(_lfs_cache);
318  _u_fine.resize(_lfs_cache.size());
319  _u_view.read(_u_fine);
320  _u_view.unbind();
321  // do projection on all leafs
322  TypeTree::applyToTree(_lfs,*this);
323  }
324  }
325  }
326  }
327 
328  backup_visitor(const GFS& gfs,
329  Projection& projection,
330  const DOFVector& u,
331  LeafOffsetCache& leaf_offset_cache,
332  TransferMap& transfer_map,
333  std::size_t int_order = 2)
334  : _lfs(gfs)
335  , _lfs_cache(_lfs)
336  , _id_set(gfs.gridView().grid().localIdSet())
337  , _projection(projection)
338  , _u_view(u)
339  , _transfer_map(transfer_map)
340  , _u_coarse(nullptr)
341  , _leaf_offset_cache(leaf_offset_cache)
342  , _int_order(int_order)
343  , _leaf_index(0)
344  {}
345 
346  LFS _lfs;
347  LFSCache _lfs_cache;
348  const IDSet& _id_set;
352  Projection& _projection;
353  typename DOFVector::template ConstLocalView<LFSCache> _u_view;
354  TransferMap& _transfer_map;
355  LocalDOFVector* _u_coarse;
356  LeafOffsetCache& _leaf_offset_cache;
357  size_type _int_order;
358  size_type _leaf_index;
359  LocalDOFVector _u_fine;
360 
361  };
362 
363 
364 
365  template<typename GFS, typename DOFVector, typename CountVector>
367  : public TypeTree::TreeVisitor
368  , public TypeTree::DynamicTraversal
369  {
370 
374 
375  using EntitySet = typename GFS::Traits::EntitySet;
376  using IDSet = typename EntitySet::Traits::GridView::Grid::LocalIdSet;
377  using Element = typename EntitySet::Element;
378  using Geometry = typename Element::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  using DF = typename EntitySet::Traits::CoordinateField;
385 
386  template<typename FiniteElement>
388  {
389  using Range = typename FiniteElement::Traits::LocalBasisType::Traits::RangeType;
390 
391  template<typename X, typename Y>
392  void evaluate(const X& x, Y& y) const
393  {
394  _phi.resize(_finite_element.localBasis().size());
395  _finite_element.localBasis().evaluateFunction(_coarse_geometry.local(_fine_geometry.global(x)),_phi);
396  y = 0;
397  for (size_type i = 0; i < _phi.size(); ++i)
398  y.axpy(_dofs[_offset + i],_phi[i]);
399  }
400 
401  coarse_function(const FiniteElement& finite_element, Geometry coarse_geometry, Geometry fine_geometry, const LocalDOFVector& dofs, size_type offset)
402  : _finite_element(finite_element)
403  , _coarse_geometry(coarse_geometry)
404  , _fine_geometry(fine_geometry)
405  , _dofs(dofs)
406  , _offset(offset)
407  {}
408 
409  const FiniteElement& _finite_element;
412  const LocalDOFVector& _dofs;
413  mutable std::vector<Range> _phi;
414  size_type _offset;
415 
416  };
417 
418 
419  template<typename LeafLFS, typename TreePath>
420  void leaf(const LeafLFS& leaf_lfs, TreePath treePath)
421  {
422  using FiniteElement = typename LeafLFS::Traits::FiniteElementType;
423 
424  auto& fem = leaf_lfs.gridFunctionSpace().finiteElementMap();
425  auto element_offset = _leaf_offset_cache[_element.type()][_leaf_index];
426  auto ancestor_offset = _leaf_offset_cache[_ancestor.type()][_leaf_index];
427 
428  coarse_function<FiniteElement> f(fem.find(_ancestor),_ancestor.geometry(),_element.geometry(),*_u_coarse,ancestor_offset);
429  auto& fe = fem.find(_element);
430 
431  _u_tmp.resize(fe.localBasis().size());
432  std::fill(_u_tmp.begin(),_u_tmp.end(),RF(0.0));
433  fe.localInterpolation().interpolate(f,_u_tmp);
434  std::copy(_u_tmp.begin(),_u_tmp.end(),_u_fine.begin() + element_offset);
435 
436  ++_leaf_index;
437  }
438 
439  void operator()(const Element& element, const Element& ancestor, const LocalDOFVector& u_coarse)
440  {
441  _element = element;
442  _ancestor = ancestor;
443  _u_coarse = &u_coarse;
444  _lfs.bind(_element);
446  _lfs_cache.update();
447  _u_view.bind(_lfs_cache);
448 
449  // test identity using ids
450  if (_id_set.id(element) == _id_set.id(ancestor))
451  {
452  // no interpolation necessary, just copy the saved data
453  _u_view.add(*_u_coarse);
454  }
455  else
456  {
457  _u_fine.resize(_lfs_cache.size());
458  std::fill(_u_fine.begin(),_u_fine.end(),RF(0));
459  _leaf_index = 0;
460  TypeTree::applyToTree(_lfs,*this);
461  _u_view.add(_u_fine);
462  }
463  _u_view.commit();
464 
465  _uc_view.bind(_lfs_cache);
466  _counts.resize(_lfs_cache.size(),1);
467  _uc_view.add(_counts);
468  _uc_view.commit();
469  }
470 
471  replay_visitor(const GFS& gfs, DOFVector& u, CountVector& uc, LeafOffsetCache& leaf_offset_cache)
472  : _lfs(gfs)
473  , _lfs_cache(_lfs)
474  , _id_set(gfs.entitySet().gridView().grid().localIdSet())
475  , _u_view(u)
476  , _uc_view(uc)
477  , _leaf_offset_cache(leaf_offset_cache)
478  , _leaf_index(0)
479  {}
480 
481  LFS _lfs;
482  LFSCache _lfs_cache;
483  const IDSet& _id_set;
486  typename DOFVector::template LocalView<LFSCache> _u_view;
487  typename CountVector::template LocalView<LFSCache> _uc_view;
488  const LocalDOFVector* _u_coarse;
489  LeafOffsetCache& _leaf_offset_cache;
490  size_type _leaf_index;
491  LocalDOFVector _u_fine;
492  LocalDOFVector _u_tmp;
493  LocalCountVector _counts;
494 
495  };
496 
497 
511  template<class Grid, class GFSU, class U, class Projection>
513  {
514  typedef typename Grid::LeafGridView LeafGridView;
515  typedef typename LeafGridView::template Codim<0>
516  ::template Partition<Dune::Interior_Partition>::Iterator LeafIterator;
517  typedef typename Grid::template Codim<0>::Entity Element;
518  typedef typename Grid::LocalIdSet IDSet;
519  typedef typename IDSet::IdType ID;
520 
521  public:
522  typedef std::unordered_map<ID,std::vector<typename U::ElementType> > MapType;
523 
524 
529  explicit GridAdaptor(const GFSU& gfs)
530  : _leaf_offset_cache(gfs)
531  {}
532 
533  /* @brief @todo
534  *
535  * @param[in] u The solution that will be saved
536  * @param[out] transferMap The map containing the solution during adaptation
537  */
538  void backupData(Grid& grid, GFSU& gfsu, Projection& projection, U& u, MapType& transfer_map)
539  {
540  typedef backup_visitor<GFSU,U,MapType> Visitor;
541 
542  Visitor visitor(gfsu,projection,u,_leaf_offset_cache,transfer_map);
543 
544  // iterate over all elems
545  for(const auto& cell : elements(gfsu.entitySet(),Partitions::interior))
546  visitor(cell);
547  }
548 
549  /* @brief @todo
550  *
551  * @param[out] u The solution after adaptation
552  * @param[in] transferMap The map that contains the information for the rebuild of u
553  */
554  void replayData(Grid& grid, GFSU& gfsu, Projection& projection, U& u, const MapType& transfer_map)
555  {
556  const IDSet& id_set = grid.localIdSet();
557 
558  using CountVector = Backend::Vector<GFSU,int>;
559  CountVector uc(gfsu,0);
560 
561  typedef replay_visitor<GFSU,U,CountVector> Visitor;
562  Visitor visitor(gfsu,u,uc,_leaf_offset_cache);
563 
564  // iterate over all elems
565  for (const auto& cell : elements(gfsu.entitySet(),Partitions::interior))
566  {
567  Element ancestor = cell;
568 
569  typename MapType::const_iterator map_it;
570  while ((map_it = transfer_map.find(id_set.id(ancestor))) == transfer_map.end())
571  {
572  if (!ancestor.hasFather())
573  DUNE_THROW(Exception,
574  "transferMap of GridAdaptor didn't contain ancestor of element with id " << id_set.id(ancestor));
575  ancestor = ancestor.father();
576  }
577 
578  visitor(cell,ancestor,map_it->second);
579  }
580 
581  typedef Dune::PDELab::AddDataHandle<GFSU,U> DOFHandle;
582  DOFHandle addHandle1(gfsu,u);
583  gfsu.entitySet().gridView().communicate(addHandle1,
584  Dune::InteriorBorder_InteriorBorder_Interface,Dune::ForwardCommunication);
586  CountHandle addHandle2(gfsu,uc);
587  gfsu.entitySet().gridView().communicate(addHandle2,
588  Dune::InteriorBorder_InteriorBorder_Interface,Dune::ForwardCommunication);
589 
590  // normalize multiple-interpolated DOFs by taking the arithmetic average
591  typename CountVector::iterator ucit = uc.begin();
592  for (typename U::iterator uit = u.begin(), uend = u.end(); uit != uend; ++uit, ++ucit)
593  (*uit) /= ((*ucit) > 0 ? (*ucit) : 1.0);
594  }
595 
596  private:
597 
599 
600  };
601 
620  template<class Grid, class GFS, class X>
621  void adapt_grid (Grid& grid, GFS& gfs, X& x1, int int_order)
622  {
623  typedef L2Projection<GFS,X> Projection;
624  Projection projection(gfs,int_order);
625 
626  GridAdaptor<Grid,GFS,X,Projection> grid_adaptor(gfs);
627 
628  // prepare the grid for refinement
629  grid.preAdapt();
630 
631  // save u
632  typename GridAdaptor<Grid,GFS,X,Projection>::MapType transferMap1;
633  grid_adaptor.backupData(grid,gfs,projection,x1,transferMap1);
634 
635  // adapt the grid
636  grid.adapt();
637 
638  // update the function spaces
639  gfs.update(true);
640 
641  // reset u
642  x1 = X(gfs,0.0);
643  grid_adaptor.replayData(grid,gfs,projection,x1,transferMap1);
644 
645  // clean up
646  grid.postAdapt();
647  }
648 
660  template<class Grid, class GFS, class X>
661  void adapt_grid (Grid& grid, GFS& gfs, X& x1, X& x2, int int_order)
662  {
663  typedef L2Projection<GFS,X> Projection;
664  Projection projection(gfs,int_order);
665 
666  GridAdaptor<Grid,GFS,X,Projection> grid_adaptor(gfs);
667 
668  // prepare the grid for refinement
669  grid.preAdapt();
670 
671  // save solution
672  typename GridAdaptor<Grid,GFS,X,Projection>::MapType transferMap1;
673  grid_adaptor.backupData(grid,gfs,projection,x1,transferMap1);
674  typename GridAdaptor<Grid,GFS,X,Projection>::MapType transferMap2;
675  grid_adaptor.backupData(grid,gfs,projection,x2,transferMap2);
676 
677  // adapt the grid
678  grid.adapt();
679 
680  // update the function spaces
681  gfs.update(true);
682 
683  // interpolate solution
684  x1 = X(gfs,0.0);
685  grid_adaptor.replayData(grid,gfs,projection,x1,transferMap1);
686  x2 = X(gfs,0.0);
687  grid_adaptor.replayData(grid,gfs,projection,x2,transferMap2);
688 
689  // clean up
690  grid.postAdapt();
691  }
692 
693 #ifndef DOXYGEN
694  namespace impl{
695 
696  // Struct for storing a GridFunctionSpace, corrosponding vectors and integration order
697  template <typename G, typename... X>
698  struct GFSWithVectors
699  {
700  // Export types
701  using GFS = G;
702  using Tuple = std::tuple<X&...>;
703 
704  GFSWithVectors (GFS& gfs, int integrationOrder, X&... x) :
705  _gfs(gfs),
706  _integrationOrder(integrationOrder),
707  _tuple(x...)
708  {}
709 
710  GFS& _gfs;
711  int _integrationOrder;
712  Tuple _tuple;
713  };
714 
715  // Forward declarations needed for the recursion
716  template <typename Grid>
717  void iteratePacks(Grid& grid);
718  template <typename Grid, typename X, typename... XS>
719  void iteratePacks(Grid& grid, X& x, XS&... xs);
720 
721  // This function is called after the last vector of the tuple. Here
722  // the next pack is called. On the way back we update the current
723  // function space.
724  template<std::size_t I = 0, typename Grid, typename X, typename... XS>
726  iterateTuple(Grid& grid, X& x, XS&... xs)
727  {
728  // Iterate next pack
729  iteratePacks(grid,xs...);
730 
731  // On our way back we need to update the current function space
732  x._gfs.update(true);
733  }
734 
735  /* In this function we store the data of the current vector (indicated
736  * by template parameter I) of the current pack. After recursively
737  * iterating through the other packs and vectors we replay the data.
738  *
739  * @tparam I std:size_t used for tmp
740  * @tparam Grid Grid type
741  * @tparam X Current pack
742  * @tparam ...XS Remaining packs
743  */
744  template<std::size_t I = 0, typename Grid, typename X, typename... XS>
746  iterateTuple(Grid& grid, X& x, XS&... xs)
747  {
748  // Get some basic types
749  using GFS = typename X::GFS;
750  using Tuple = typename X::Tuple;
751  using V = typename std::decay<typename std::tuple_element<I,Tuple>::type>::type;
752  // // alternative:
753  // auto v = std::get<I>(x._tuple);
754  // using V = decltype(v);
755 
756  // Setup classes for data restoring
757  typedef Dune::PDELab::L2Projection <GFS,V> Projection;
758  Projection projection(x._gfs,x._integrationOrder);
759  GridAdaptor<Grid,GFS,V,Projection> gridAdaptor(x._gfs);
760 
761  // Store vector data
762  typename GridAdaptor<Grid,GFS,V,Projection>::MapType transferMap;
763  gridAdaptor.backupData(grid,x._gfs,projection,std::get<I>(x._tuple),transferMap);
764 
765  // Recursively iterate through remaining vectors (and packs). Grid
766  // adaption will be done at the end of recursion.
767  iterateTuple<I + 1, Grid, X, XS...>(grid,x,xs...);
768 
769  // Play back data. Note: At this point the function space was
770  // already updatet.
771  std::get<I>(x._tuple) = V(x._gfs,0.0);
772  gridAdaptor.replayData(grid,x._gfs,projection,std::get<I>(x._tuple),transferMap);
773  }
774 
775  // This gets called after the last pack. After this function call we
776  // have visited every vector of every pack and we will go back through
777  // the recursive function calls.
778  template <typename Grid>
779  void iteratePacks(Grid& grid)
780  {
781  // Adapt the grid
782  grid.adapt();
783  }
784 
785  /* Use template meta programming to iterate over packs at compile time
786  *
787  * In order to adapt our grid and all vectors of all packs we need to
788  * do the following:
789  * - Iterate over all vectors of all packs.
790  * - Store the data from the vectors where things could change.
791  * - Adapt our grid.
792  * - Update function spaces and restore data.
793  *
794  * The key point is that we need the object that stores the data to
795  * replay it. Because of that we can not just iterate over the packs
796  * and within each pack iterate over the vectors but we have to make
797  * one big recursion. Therefore we iterate over the vectors of the
798  * current pack.
799  */
800  template <typename Grid, typename X, typename... XS>
801  void iteratePacks(Grid& grid, X& x, XS&... xs)
802  {
803  iterateTuple(grid,x,xs...);
804  }
805 
806  } // namespace impl
807 #endif // DOXYGEN
808 
820  template <typename GFS, typename... X>
821  impl::GFSWithVectors<GFS,X...> transferSolutions(GFS& gfs, int integrationOrder, X&... x)
822  {
823  impl::GFSWithVectors<GFS,X...> gfsWithVectors(gfs, integrationOrder, x...);
824  return gfsWithVectors;
825  }
826 
837  template <typename Grid, typename... X>
838  void adaptGrid(Grid& grid, X&... x)
839  {
840  // Prepare the grid for refinement
841  grid.preAdapt();
842 
843  // Iterate over packs
844  impl::iteratePacks(grid,x...);
845 
846  // Clean up
847  grid.postAdapt();
848  }
849 
850 
851  template<typename T>
852  void error_fraction(const T& x, typename T::ElementType alpha, typename T::ElementType beta,
853  typename T::ElementType& eta_alpha, typename T::ElementType& eta_beta, int verbose=0)
854  {
855  if (verbose>0)
856  std::cout << "+++ error fraction: alpha=" << alpha << " beta=" << beta << std::endl;
857  const int steps=20; // max number of bisection steps
858  typedef typename T::ElementType NumberType;
859  NumberType total_error = x.one_norm();
860  NumberType max_error = x.infinity_norm();
861  NumberType eta_alpha_left = 0.0;
862  NumberType eta_alpha_right = max_error;
863  NumberType eta_beta_left = 0.0;
864  NumberType eta_beta_right = max_error;
865  for (int j=1; j<=steps; j++)
866  {
867  eta_alpha = 0.5*(eta_alpha_left+eta_alpha_right);
868  eta_beta = 0.5*(eta_beta_left+eta_beta_right);
869  NumberType sum_alpha=0.0;
870  NumberType sum_beta=0.0;
871  unsigned int alpha_count = 0;
872  unsigned int beta_count = 0;
873  for (typename T::const_iterator it = x.begin(),
874  end = x.end();
875  it != end;
876  ++it)
877  {
878  if (*it >=eta_alpha) { sum_alpha += *it; alpha_count++;}
879  if (*it < eta_beta) { sum_beta += *it; beta_count++;}
880  }
881  if (verbose>1)
882  {
883  std::cout << "+++ " << j << " eta_alpha=" << eta_alpha << " alpha_fraction=" << sum_alpha/total_error
884  << " elements: " << alpha_count << " of " << x.N() << std::endl;
885  std::cout << "+++ " << j << " eta_beta=" << eta_beta << " beta_fraction=" << sum_beta/total_error
886  << " elements: " << beta_count << " of " << x.N() << std::endl;
887  }
888  if (std::abs(alpha-sum_alpha/total_error) <= 0.01 && std::abs(beta-sum_beta/total_error) <= 0.01) break;
889  if (sum_alpha>alpha*total_error)
890  eta_alpha_left = eta_alpha;
891  else
892  eta_alpha_right = eta_alpha;
893  if (sum_beta>beta*total_error)
894  eta_beta_right = eta_beta;
895  else
896  eta_beta_left = eta_beta;
897  }
898  if (verbose>0)
899  {
900  std::cout << "+++ refine_threshold=" << eta_alpha
901  << " coarsen_threshold=" << eta_beta << std::endl;
902  }
903  }
904 
905 
906  template<typename T>
907  void element_fraction(const T& x, typename T::ElementType alpha, typename T::ElementType beta,
908  typename T::ElementType& eta_alpha, typename T::ElementType& eta_beta, int verbose=0)
909  {
910  const int steps=20; // max number of bisection steps
911  typedef typename T::ElementType NumberType;
912  NumberType total_error =x.N();
913  NumberType max_error = x.infinity_norm();
914  NumberType eta_alpha_left = 0.0;
915  NumberType eta_alpha_right = max_error;
916  NumberType eta_beta_left = 0.0;
917  NumberType eta_beta_right = max_error;
918  for (int j=1; j<=steps; j++)
919  {
920  eta_alpha = 0.5*(eta_alpha_left+eta_alpha_right);
921  eta_beta = 0.5*(eta_beta_left+eta_beta_right);
922  NumberType sum_alpha=0.0;
923  NumberType sum_beta=0.0;
924  unsigned int alpha_count = 0;
925  unsigned int beta_count = 0;
926 
927  for (typename T::const_iterator it = x.begin(),
928  end = x.end();
929  it != end;
930  ++it)
931  {
932  if (*it>=eta_alpha) { sum_alpha += 1.0; alpha_count++;}
933  if (*it< eta_beta) { sum_beta +=1.0; beta_count++;}
934  }
935  if (verbose>1)
936  {
937  std::cout << j << " eta_alpha=" << eta_alpha << " alpha_fraction=" << sum_alpha/total_error
938  << " elements: " << alpha_count << " of " << x.N() << std::endl;
939  std::cout << j << " eta_beta=" << eta_beta << " beta_fraction=" << sum_beta/total_error
940  << " elements: " << beta_count << " of " << x.N() << std::endl;
941  }
942  if (std::abs(alpha-sum_alpha/total_error) <= 0.01 && std::abs(beta-sum_beta/total_error) <= 0.01) break;
943  if (sum_alpha>alpha*total_error)
944  eta_alpha_left = eta_alpha;
945  else
946  eta_alpha_right = eta_alpha;
947  if (sum_beta>beta*total_error)
948  eta_beta_right = eta_beta;
949  else
950  eta_beta_left = eta_beta;
951  }
952  if (verbose>0)
953  {
954  std::cout << "+++ refine_threshold=" << eta_alpha
955  << " coarsen_threshold=" << eta_beta << std::endl;
956  }
957  }
958 
961  template<typename T>
962  void error_distribution(const T& x, unsigned int bins)
963  {
964  const int steps=30; // max number of bisection steps
965  typedef typename T::ElementType NumberType;
966  NumberType total_error = x.one_norm();
967  NumberType total_elements = x.N();
968  NumberType max_error = x.infinity_norm();
969  std::vector<NumberType> left(bins,0.0);
970  std::vector<NumberType> right(bins,max_error*(1.0+1e-8));
971  std::vector<NumberType> eta(bins);
972  std::vector<NumberType> target(bins);
973  for (unsigned int k=0; k<bins; k++)
974  target[k]= (k+1)/((NumberType)bins);
975  for (int j=1; j<=steps; j++)
976  {
977  for (unsigned int k=0; k<bins; k++)
978  eta[k]= 0.5*(left[k]+right[k]);
979  std::vector<NumberType> sum(bins,0.0);
980  std::vector<int> count(bins,0);
981 
982  for (typename T::const_iterator it = x.begin(),
983  end = x.end();
984  it != end;
985  ++it)
986  {
987  for (unsigned int k=0; k<bins; k++)
988  if (*it<=eta[k])
989  {
990  sum[k] += *it;
991  count[k] += 1;
992  }
993  }
994  // std::cout << std::endl;
995  // std::cout << "// step " << j << std::endl;
996  // for (unsigned int k=0; k<bins; k++)
997  // std::cout << k+1 << " " << count[k] << " " << eta[k] << " " << right[k]-left[k]
998  // << " " << sum[k]/total_error << " " << target[k] << std::endl;
999  for (unsigned int k=0; k<bins; k++)
1000  if (sum[k]<=target[k]*total_error)
1001  left[k] = eta[k];
1002  else
1003  right[k] = eta[k];
1004  }
1005  std::vector<NumberType> sum(bins,0.0);
1006  std::vector<int> count(bins,0);
1007  for (unsigned int i=0; i<x.N(); i++)
1008  for (unsigned int k=0; k<bins; k++)
1009  if (x[i]<=eta[k])
1010  {
1011  sum[k] += x[i];
1012  count[k] += 1;
1013  }
1014  std::cout << "+++ error distribution" << std::endl;
1015  std::cout << "+++ number of elements: " << x.N() << std::endl;
1016  std::cout << "+++ max element error: " << max_error << std::endl;
1017  std::cout << "+++ total error: " << total_error << std::endl;
1018  std::cout << "+++ bin #elements eta sum/total " << std::endl;
1019  for (unsigned int k=0; k<bins; k++)
1020  std::cout << "+++ " << k+1 << " " << count[k] << " " << eta[k] << " " << sum[k]/total_error << std::endl;
1021  }
1022 
1023  template<typename Grid, typename X>
1024  void mark_grid (Grid &grid, const X& x, typename X::ElementType refine_threshold,
1025  typename X::ElementType coarsen_threshold, int min_level = 0, int max_level = std::numeric_limits<int>::max(), int verbose=0)
1026  {
1027  typedef typename Grid::LeafGridView GV;
1028 
1029  GV gv = grid.leafGridView();
1030 
1031  unsigned int refine_cnt=0;
1032  unsigned int coarsen_cnt=0;
1033 
1034  typedef typename X::GridFunctionSpace GFS;
1035  typedef LocalFunctionSpace<GFS> LFS;
1036  typedef LFSIndexCache<LFS> LFSCache;
1037  typedef typename X::template ConstLocalView<LFSCache> XView;
1038 
1039  LFS lfs(x.gridFunctionSpace());
1040  LFSCache lfs_cache(lfs);
1041  XView x_view(x);
1042 
1043  for(const auto& cell : elements(gv))
1044  {
1045  lfs.bind(cell);
1046  lfs_cache.update();
1047  x_view.bind(lfs_cache);
1048 
1049  if (x_view[0]>=refine_threshold && cell.level() < max_level)
1050  {
1051  grid.mark(1,cell);
1052  refine_cnt++;
1053  }
1054  if (x_view[0]<=coarsen_threshold && cell.level() > min_level)
1055  {
1056  grid.mark(-1,cell);
1057  coarsen_cnt++;
1058  }
1059  x_view.unbind();
1060  }
1061  if (verbose>0)
1062  std::cout << "+++ mark_grid: " << refine_cnt << " marked for refinement, "
1063  << coarsen_cnt << " marked for coarsening" << std::endl;
1064  }
1065 
1066 
1067  template<typename Grid, typename X>
1068  void mark_grid_for_coarsening (Grid &grid, const X& x, typename X::ElementType refine_threshold,
1069  typename X::ElementType coarsen_threshold, int verbose=0)
1070  {
1071  typedef typename Grid::LeafGridView GV;
1072 
1073  GV gv = grid.leafGridView();
1074 
1075  unsigned int coarsen_cnt=0;
1076 
1077  typedef typename X::GridFunctionSpace GFS;
1078  typedef LocalFunctionSpace<GFS> LFS;
1079  typedef LFSIndexCache<LFS> LFSCache;
1080  typedef typename X::template ConstLocalView<LFSCache> XView;
1081 
1082  LFS lfs(x.gridFunctionSpace());
1083  LFSCache lfs_cache(lfs);
1084  XView x_view(x);
1085 
1086  for(const auto& cell : elements(gv))
1087  {
1088  lfs.bind(cell);
1089  lfs_cache.update();
1090  x_view.bind(lfs_cache);
1091 
1092  if (x_view[0]>=refine_threshold)
1093  {
1094  grid.mark(-1,cell);
1095  coarsen_cnt++;
1096  }
1097  if (x_view[0]<=coarsen_threshold)
1098  {
1099  grid.mark(-1,cell);
1100  coarsen_cnt++;
1101  }
1102  }
1103  if (verbose>0)
1104  std::cout << "+++ mark_grid_for_coarsening: "
1105  << coarsen_cnt << " marked for coarsening" << std::endl;
1106  }
1107 
1108 
1110  {
1111  // strategy parameters
1112  double scaling;
1113  double optimistic_factor;
1114  double coarsen_limit;
1115  double balance_limit;
1116  double tol;
1117  double T;
1118  int verbose;
1119  bool no_adapt;
1120  double refine_fraction_while_refinement;
1121  double coarsen_fraction_while_refinement;
1122  double coarsen_fraction_while_coarsening;
1123  double timestep_decrease_factor;
1124  double timestep_increase_factor;
1125 
1126  // results to be reported to the user after evaluating the error
1127  bool accept;
1128  bool adapt_dt;
1129  bool adapt_grid;
1130  double newdt;
1131  double q_s, q_t;
1132 
1133  // state variables
1134  bool have_decreased_time_step;
1135  bool have_refined_grid;
1136 
1137  // the only state variable: accumulated error
1138  double accumulated_estimated_error_squared;
1139  double minenergy_rate;
1140 
1141  public:
1142  TimeAdaptationStrategy (double tol_, double T_, int verbose_=0)
1143  : scaling(16.0), optimistic_factor(1.0), coarsen_limit(0.5), balance_limit(0.33333),
1144  tol(tol_), T(T_), verbose(verbose_), no_adapt(false),
1145  refine_fraction_while_refinement(0.7),
1146  coarsen_fraction_while_refinement(0.2),
1147  coarsen_fraction_while_coarsening(0.2),
1148  timestep_decrease_factor(0.5), timestep_increase_factor(1.5),
1149  accept(false), adapt_dt(false), adapt_grid(false), newdt(1.0),
1150  have_decreased_time_step(false), have_refined_grid(false),
1151  accumulated_estimated_error_squared(0.0),
1152  minenergy_rate(0.0)
1153  {
1154  }
1155 
1157  {
1158  timestep_decrease_factor=s;
1159  }
1160 
1162  {
1163  timestep_increase_factor=s;
1164  }
1165 
1167  {
1168  refine_fraction_while_refinement=s;
1169  }
1170 
1172  {
1173  coarsen_fraction_while_refinement=s;
1174  }
1175 
1177  {
1178  coarsen_fraction_while_coarsening=s;
1179  }
1180 
1181  void setMinEnergyRate (double s)
1182  {
1183  minenergy_rate=s;
1184  }
1185 
1186  void setCoarsenLimit (double s)
1187  {
1188  coarsen_limit=s;
1189  }
1190 
1191  void setBalanceLimit (double s)
1192  {
1193  balance_limit=s;
1194  }
1195 
1196  void setTemporalScaling (double s)
1197  {
1198  scaling=s;
1199  }
1200 
1201  void setOptimisticFactor (double s)
1202  {
1203  optimistic_factor=s;
1204  }
1205 
1207  {
1208  no_adapt = false;
1209  }
1210 
1212  {
1213  no_adapt = true;
1214  }
1215 
1216  bool acceptTimeStep () const
1217  {
1218  return accept;
1219  }
1220 
1221  bool adaptDT () const
1222  {
1223  return adapt_dt;
1224  }
1225 
1226  bool adaptGrid () const
1227  {
1228  return adapt_grid;
1229  }
1230 
1231  double newDT () const
1232  {
1233  return newdt;
1234  }
1235 
1236  double qs () const
1237  {
1238  return q_s;
1239  }
1240 
1241  double qt () const
1242  {
1243  return q_t;
1244  }
1245 
1246  double endT() const
1247  {
1248  return T;
1249  }
1250 
1251  double accumulatedErrorSquared () const
1252  {
1253  return accumulated_estimated_error_squared;
1254  }
1255 
1256  // to be called when new time step is done
1258  {
1259  have_decreased_time_step = false;
1260  have_refined_grid = false;
1261  }
1262 
1263  template<typename GM, typename X>
1264  void evaluate_estimators (GM& grid, double time, double dt, const X& eta_space, const X& eta_time, double energy_timeslab)
1265  {
1266  accept=false;
1267  adapt_dt=false;
1268  adapt_grid=false;
1269  newdt=dt;
1270 
1271  double spatial_error = eta_space.one_norm();
1272  double temporal_error = scaling*eta_time.one_norm();
1273  double sum = spatial_error + temporal_error;
1274  //double allowed = optimistic_factor*(tol*tol-accumulated_estimated_error_squared)*dt/(T-time);
1275  double allowed = tol*tol*(energy_timeslab+minenergy_rate*dt);
1276  q_s = spatial_error/sum;
1277  q_t = temporal_error/sum;
1278 
1279  // print some statistics
1280  if (verbose>0)
1281  std::cout << "+++"
1282  << " q_s=" << q_s
1283  << " q_t=" << q_t
1284  << " sum/allowed=" << sum/allowed
1285  // << " allowed=" << allowed
1286  << " estimated error=" << sqrt(accumulated_estimated_error_squared+sum)
1287  << " energy_rate=" << energy_timeslab/dt
1288  << std::endl;
1289 
1290  // for simplicity: a mode that does no adaptation at all
1291  if (no_adapt)
1292  {
1293  accept = true;
1294  accumulated_estimated_error_squared += sum;
1295  if (verbose>1) std::cout << "+++ no adapt mode" << std::endl;
1296  return;
1297  }
1298 
1299  // the adaptation strategy
1300  if (sum<=allowed)
1301  {
1302  // we will accept this time step
1303  accept = true;
1304  if (verbose>1) std::cout << "+++ accepting time step" << std::endl;
1305  accumulated_estimated_error_squared += sum;
1306 
1307  // check if grid size or time step needs to be adapted
1308  if (sum<coarsen_limit*allowed)
1309  {
1310  // the error is too small, i.e. the computation is inefficient
1311  if (q_t<balance_limit)
1312  {
1313  // spatial error is dominating => increase time step
1314  newdt = timestep_increase_factor*dt;
1315  adapt_dt = true;
1316  if (verbose>1) std::cout << "+++ spatial error dominates: increase time step" << std::endl;
1317  }
1318  else
1319  {
1320  if (q_s>balance_limit)
1321  {
1322  // step sizes balanced: coarsen in time
1323  newdt = timestep_increase_factor*dt;
1324  adapt_dt = true;
1325  if (verbose>1) std::cout << "+++ increasing time step" << std::endl;
1326  }
1327  // coarsen grid in space
1328  double eta_refine, eta_coarsen;
1329  if (verbose>1) std::cout << "+++ mark grid for coarsening" << std::endl;
1330  //error_distribution(eta_space,20);
1331  Dune::PDELab::error_fraction(eta_space,coarsen_fraction_while_coarsening,
1332  coarsen_fraction_while_coarsening,eta_refine,eta_coarsen);
1333  Dune::PDELab::mark_grid_for_coarsening(grid,eta_space,eta_refine,eta_coarsen,verbose);
1334  adapt_grid = true;
1335  }
1336  }
1337  else
1338  {
1339  // modify at least the time step
1340  if (q_t<balance_limit)
1341  {
1342  // spatial error is dominating => increase time step
1343  newdt = timestep_increase_factor*dt;
1344  adapt_dt = true;
1345  if (verbose>1) std::cout << "+++ spatial error dominates: increase time step" << std::endl;
1346  }
1347  }
1348  }
1349  else
1350  {
1351  // error is too large, we need to do something
1352  if (verbose>1) std::cout << "+++ will redo time step" << std::endl;
1353  if (q_t>1-balance_limit)
1354  {
1355  // temporal error is dominating => decrease time step only
1356  newdt = timestep_decrease_factor*dt;
1357  adapt_dt = true;
1358  have_decreased_time_step = true;
1359  if (verbose>1) std::cout << "+++ decreasing time step only" << std::endl;
1360  }
1361  else
1362  {
1363  if (q_t<balance_limit)
1364  {
1365  if (!have_decreased_time_step)
1366  {
1367  // time steps size not balanced (too small)
1368  newdt = timestep_increase_factor*dt;
1369  adapt_dt = true;
1370  if (verbose>1) std::cout << "+++ increasing time step" << std::endl;
1371  }
1372  }
1373  else
1374  {
1375  // step sizes balanced: refine in time as well
1376  newdt = timestep_decrease_factor*dt;
1377  adapt_dt = true;
1378  have_decreased_time_step = true;
1379  if (verbose>1) std::cout << "+++ decreasing time step" << std::endl;
1380  }
1381  // refine grid in space
1382  double eta_refine, eta_coarsen;
1383  if (verbose>1) std::cout << "+++ BINGO mark grid for refinement and coarsening" << std::endl;
1384  //error_distribution(eta_space,20);
1385  Dune::PDELab::error_fraction(eta_space,refine_fraction_while_refinement,
1386  coarsen_fraction_while_refinement,eta_refine,eta_coarsen,0);
1387  Dune::PDELab::mark_grid(grid,eta_space,eta_refine,eta_coarsen,verbose);
1388  adapt_grid = true;
1389  }
1390  }
1391  }
1392  };
1393 
1394 
1395 
1396  } // namespace PDELab
1397 } // namespace Dune
1398 
1399 #endif // DUNE_PDELAB_ADAPTIVITY_ADAPTIVITY_HH
Element _current
Definition: adaptivity.hh:351
std::vector< typename CountVector::ElementType > LocalCountVector
Definition: adaptivity.hh:381
void setMinEnergyRate(double s)
Definition: adaptivity.hh:1181
TransferMap::mapped_type LocalDOFVector
Definition: adaptivity.hh:214
CountVector::template LocalView< LFSCache > _uc_view
Definition: adaptivity.hh:487
Element _ancestor
Definition: adaptivity.hh:485
void setAdaptationOn()
Definition: adaptivity.hh:1206
const QuadratureRule< DF, dim > & _quadrature_rule
Definition: adaptivity.hh:128
const LocalDOFVector * _u_coarse
Definition: adaptivity.hh:488
static const int dim
Definition: adaptivity.hh:83
void replayData(Grid &grid, GFSU &gfsu, Projection &projection, U &u, const MapType &transfer_map)
Definition: adaptivity.hh:554
const Entity & e
Definition: localfunctionspace.hh:111
double accumulatedErrorSquared() const
Definition: adaptivity.hh:1251
LFS _lfs
Definition: adaptivity.hh:481
std::unordered_map< ID, std::vector< typename U::ElementType > > MapType
Definition: adaptivity.hh:522
LFSIndexCache< LFS > LFSCache
Definition: adaptivity.hh:372
std::size_t size_type
Definition: adaptivity.hh:383
Geometry _coarse_geometry
Definition: adaptivity.hh:410
typename Element::Geometry Geometry
Definition: adaptivity.hh:378
LocalFunctionSpace< GFS > LFS
Definition: adaptivity.hh:371
LocalDOFVector _u_fine
Definition: adaptivity.hh:491
void setCoarsenLimit(double s)
Definition: adaptivity.hh:1186
void setBalanceLimit(double s)
Definition: adaptivity.hh:1191
const GFS & gridFunctionSpace() const
Returns the GridFunctionSpace underlying this LocalFunctionSpace.
Definition: localfunctionspace.hh:264
Dune::PDELab::LeafOffsetCache< GFS > LeafOffsetCache
Definition: adaptivity.hh:206
double qt() const
Definition: adaptivity.hh:1241
const std::string s
Definition: function.hh:830
LocalFunctionSpace< GFS > LFS
Definition: adaptivity.hh:204
LFS _lfs
Definition: adaptivity.hh:69
LFS _lfs
Definition: adaptivity.hh:346
typename EntitySet::Traits::CoordinateField DF
Definition: adaptivity.hh:384
Traits::IndexContainer::size_type size() const
get current size
Definition: localfunctionspace.hh:206
typename GFS::Traits::EntitySet EntitySet
Definition: adaptivity.hh:375
void adapt_grid(Grid &grid, GFS &gfs, X &x1, int int_order)
adapt a grid, corresponding function space and solution vectors
Definition: adaptivity.hh:621
void operator()(const Element &element)
Definition: adaptivity.hh:273
size_type _int_order
Definition: adaptivity.hh:357
const IDSet & _id_set
Definition: adaptivity.hh:483
const Cell & _element
Definition: adaptivity.hh:126
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:907
typename impl::BackendVectorSelector< GridFunctionSpace, FieldType >::Type Vector
alias of the return type of BackendVectorSelector
Definition: backend/interface.hh:106
Element::Geometry Geometry
Definition: adaptivity.hh:211
bool acceptTimeStep() const
Definition: adaptivity.hh:1216
DOFVector::ElementType RF
Definition: adaptivity.hh:379
L2Projection(const GFS &gfs, int intorder=2)
The constructor.
Definition: adaptivity.hh:160
void startTimeStep()
Definition: adaptivity.hh:1257
const LocalDOFVector & _dofs
Definition: adaptivity.hh:412
DOFVector::ElementType RF
Definition: adaptivity.hh:213
TimeAdaptationStrategy(double tol_, double T_, int verbose_=0)
Definition: adaptivity.hh:1142
size_type _leaf_index
Definition: adaptivity.hh:490
replay_visitor(const GFS &gfs, DOFVector &u, CountVector &uc, LeafOffsetCache &leaf_offset_cache)
Definition: adaptivity.hh:471
LocalFunctionSpace< GFS > LFS
Definition: adaptivity.hh:36
typename EntitySet::Traits::GridView::Grid::LocalIdSet IDSet
Definition: adaptivity.hh:209
Projection & _projection
Definition: adaptivity.hh:352
Definition: adaptivity.hh:32
void setCoarsenFractionWhileRefinement(double s)
Definition: adaptivity.hh:1171
Projection::MassMatrix MassMatrix
Definition: adaptivity.hh:219
LFSCache _lfs_cache
Definition: adaptivity.hh:482
Base class for all PDELab exceptions.
Definition: exceptions.hh:17
const std::size_t offset
Definition: localfunctionspace.hh:74
Dune::PDELab::LeafOffsetCache< GFS > LeafOffsetCache
Definition: adaptivity.hh:373
std::vector< LeafOffsets > _leaf_offset_cache
Definition: adaptivity.hh:70
Element _ancestor
Definition: adaptivity.hh:350
bool adaptDT() const
Definition: adaptivity.hh:1221
MassMatrices & _mass_matrices
Definition: adaptivity.hh:127
size_type _leaf_index
Definition: adaptivity.hh:129
size_type _offset
Definition: adaptivity.hh:414
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:852
void backupData(Grid &grid, GFSU &gfsu, Projection &projection, U &u, MapType &transfer_map)
Definition: adaptivity.hh:538
For backward compatibility – Do not use this!
Definition: adaptivity.hh:27
std::vector< RF > LocalDOFVector
Definition: adaptivity.hh:380
Geometry _fine_geometry
Definition: adaptivity.hh:411
void setTemporalScaling(double s)
Definition: adaptivity.hh:1196
LFSCache _lfs_cache
Definition: adaptivity.hh:347
void setTimeStepIncreaseFactor(double s)
Definition: adaptivity.hh:1161
void adaptGrid(Grid &grid, X &... x)
Adapt grid and multiple function spaces with corresponding vectors.
Definition: adaptivity.hh:838
const FiniteElement & _finite_element
Definition: adaptivity.hh:409
LocalDOFVector _u_tmp
Definition: adaptivity.hh:492
void evaluate_estimators(GM &grid, double time, double dt, const X &eta_space, const X &eta_time, double energy_timeslab)
Definition: adaptivity.hh:1264
const MassMatrices & inverseMassMatrices(const Element &e)
Calculate the inverse local mass matrix, used in the local L2 projection.
Definition: adaptivity.hh:171
coarse_function(const FiniteElement &finite_element, Geometry coarse_geometry, Geometry fine_geometry, const LocalDOFVector &dofs, size_type offset)
Definition: adaptivity.hh:401
void setAdaptationOff()
Definition: adaptivity.hh:1211
typename FiniteElement::Traits::LocalBasisType::Traits::RangeType Range
Definition: adaptivity.hh:389
typename EntitySet::Traits::CoordinateField DF
Definition: adaptivity.hh:222
impl::GFSWithVectors< GFS, X... > transferSolutions(GFS &gfs, int integrationOrder, X &... x)
Pack function space and vectors for grid adaption.
Definition: adaptivity.hh:821
void error_distribution(const T &x, unsigned int bins)
Definition: adaptivity.hh:962
void operator()(const Element &element, const Element &ancestor, const LocalDOFVector &u_coarse)
Definition: adaptivity.hh:439
void leaf(const LeafLFS &leaf_lfs, TreePath treePath)
Definition: adaptivity.hh:420
typename EntitySet::Element Element
Definition: adaptivity.hh:210
bool adaptGrid() const
Definition: adaptivity.hh:1226
Definition: adaptivity.hh:144
DynamicMatrix< typename U::ElementType > MassMatrix
Definition: adaptivity.hh:153
void leaf(const LFSLeaf &leaf_lfs, TreePath treePath)
Definition: adaptivity.hh:225
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:328
void setTimeStepDecreaseFactor(double s)
Definition: adaptivity.hh:1156
void setCoarsenFractionWhileCoarsening(double s)
Definition: adaptivity.hh:1176
LocalDOFVector * _u_coarse
Definition: adaptivity.hh:355
LocalCountVector _counts
Definition: adaptivity.hh:493
void setRefineFractionWhileRefinement(double s)
Definition: adaptivity.hh:1166
DOFVector::template ConstLocalView< LFSCache > _u_view
Definition: adaptivity.hh:353
LeafOffsetCache(const GFS &gfs)
Definition: adaptivity.hh:64
Definition: adaptivity.hh:1109
std::vector< Range > _phi
Definition: adaptivity.hh:413
typename EntitySet::Element Element
Definition: adaptivity.hh:377
Definition: adaptivity.hh:366
std::array< MassMatrix, TypeTree::TreeInfo< GFS >::leafCount > MassMatrices
Definition: adaptivity.hh:154
GFS::Traits::GridView::template Codim< 0 >::Entity Cell
Definition: adaptivity.hh:35
L2Projection< typename LFS::Traits::GridFunctionSpace, DOFVector > Projection
Definition: adaptivity.hh:217
double qs() const
Definition: adaptivity.hh:1236
double newDT() const
Definition: adaptivity.hh:1231
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:1068
Element _element
Definition: adaptivity.hh:484
Element _element
Definition: adaptivity.hh:349
std::size_t size_type
Definition: adaptivity.hh:221
GridAdaptor(const GFSU &gfs)
The constructor.
Definition: adaptivity.hh:529
array< std::size_t, TypeTree::TreeInfo< GFS >::leafCount+1 > LeafOffsets
Definition: adaptivity.hh:40
double endT() const
Definition: adaptivity.hh:1246
LFSIndexCache< LFS > LFSCache
Definition: adaptivity.hh:205
TransferMap & _transfer_map
Definition: adaptivity.hh:354
Iterator extract_lfs_leaf_sizes(const LFS &lfs, Iterator it)
Definition: lfsindexcache.hh:165
Class for automatic adaptation of the grid.
Definition: adaptivity.hh:512
std::size_t index
Definition: interpolate.hh:118
void setOptimisticFactor(double s)
Definition: adaptivity.hh:1201
Definition: adaptivity.hh:199
Definition: genericdatahandle.hh:622
static const unsigned int value
Definition: gridfunctionspace/tags.hh:139
LeafOffsetCache & _leaf_offset_cache
Definition: adaptivity.hh:489
typename EntitySet::Traits::GridView::Grid::LocalIdSet IDSet
Definition: adaptivity.hh:376
LeafOffsetCache & _leaf_offset_cache
Definition: adaptivity.hh:356
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:1024
void update(const Cell &e)
Definition: adaptivity.hh:50
const IDSet & _id_set
Definition: adaptivity.hh:348
typename GFS::Traits::EntitySet EntitySet
Definition: adaptivity.hh:208
const LeafOffsets & operator[](GeometryType gt) const
Definition: adaptivity.hh:42
LocalDOFVector _u_fine
Definition: adaptivity.hh:359
size_type _leaf_index
Definition: adaptivity.hh:358
Projection::MassMatrices MassMatrices
Definition: adaptivity.hh:218
DOFVector::template LocalView< LFSCache > _u_view
Definition: adaptivity.hh:486
void evaluate(const X &x, Y &y) const
Definition: adaptivity.hh:392