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 
612  template<class Grid, class GFS, class X>
613  void adapt_grid (Grid& grid, GFS& gfs, X& x1, int int_order)
614  {
615  typedef L2Projection<GFS,X> Projection;
616  Projection projection(gfs,int_order);
617 
618  GridAdaptor<Grid,GFS,X,Projection> grid_adaptor(gfs);
619 
620  // prepare the grid for refinement
621  grid.preAdapt();
622 
623  // save u
624  typename GridAdaptor<Grid,GFS,X,Projection>::MapType transferMap1;
625  grid_adaptor.backupData(grid,gfs,projection,x1,transferMap1);
626 
627  // adapt the grid
628  grid.adapt();
629 
630  // update the function spaces
631  gfs.update(true);
632 
633  // reset u
634  x1 = X(gfs,0.0);
635  grid_adaptor.replayData(grid,gfs,projection,x1,transferMap1);
636 
637  // clean up
638  grid.postAdapt();
639  }
640 
652  template<class Grid, class GFS, class X>
653  void adapt_grid (Grid& grid, GFS& gfs, X& x1, X& x2, int int_order)
654  {
655  typedef L2Projection<GFS,X> Projection;
656  Projection projection(gfs,int_order);
657 
658  GridAdaptor<Grid,GFS,X,Projection> grid_adaptor(gfs);
659 
660  // prepare the grid for refinement
661  grid.preAdapt();
662 
663  // save solution
664  typename GridAdaptor<Grid,GFS,X,Projection>::MapType transferMap1;
665  grid_adaptor.backupData(grid,gfs,projection,x1,transferMap1);
666  typename GridAdaptor<Grid,GFS,X,Projection>::MapType transferMap2;
667  grid_adaptor.backupData(grid,gfs,projection,x2,transferMap2);
668 
669  // adapt the grid
670  grid.adapt();
671 
672  // update the function spaces
673  gfs.update(true);
674 
675  // interpolate solution
676  x1 = X(gfs,0.0);
677  grid_adaptor.replayData(grid,gfs,projection,x1,transferMap1);
678  x2 = X(gfs,0.0);
679  grid_adaptor.replayData(grid,gfs,projection,x2,transferMap2);
680 
681  // clean up
682  grid.postAdapt();
683  }
684 
685 #ifndef DOXYGEN
686  namespace impl{
687 
688  // Struct for storing a GridFunctionSpace, corrosponding vectors and integration order
689  template <typename G, typename... X>
690  struct GFSWithVectors
691  {
692  // Export types
693  using GFS = G;
694  using Tuple = std::tuple<X&...>;
695 
696  GFSWithVectors (GFS& gfs, int integrationOrder, X&... x) :
697  _gfs(gfs),
698  _integrationOrder(integrationOrder),
699  _tuple(x...)
700  {}
701 
702  GFS& _gfs;
703  int _integrationOrder;
704  Tuple _tuple;
705  };
706 
707  // Forward declarations needed for the recursion
708  template <typename Grid>
709  void iteratePacks(Grid& grid);
710  template <typename Grid, typename X, typename... XS>
711  void iteratePacks(Grid& grid, X& x, XS&... xs);
712 
713  // This function is called after the last vector of the tuple. Here
714  // the next pack is called. On the way back we update the current
715  // function space.
716  template<std::size_t I = 0, typename Grid, typename X, typename... XS>
718  iterateTuple(Grid& grid, X& x, XS&... xs)
719  {
720  // Iterate next pack
721  iteratePacks(grid,xs...);
722 
723  // On our way back we need to update the current function space
724  x._gfs.update(true);
725  }
726 
727  /* In this function we store the data of the current vector (indicated
728  * by template parameter I) of the current pack. After recursively
729  * iterating through the other packs and vectors we replay the data.
730  *
731  * @tparam I std:size_t used for tmp
732  * @tparam Grid Grid type
733  * @tparam X Current pack
734  * @tparam ...XS Remaining packs
735  */
736  template<std::size_t I = 0, typename Grid, typename X, typename... XS>
738  iterateTuple(Grid& grid, X& x, XS&... xs)
739  {
740  // Get some basic types
741  using GFS = typename X::GFS;
742  using Tuple = typename X::Tuple;
743  using V = typename std::decay<typename std::tuple_element<I,Tuple>::type>::type;
744  // // alternative:
745  // auto v = std::get<I>(x._tuple);
746  // using V = decltype(v);
747 
748  // Setup classes for data restoring
749  typedef Dune::PDELab::L2Projection <GFS,V> Projection;
750  Projection projection(x._gfs,x._integrationOrder);
751  GridAdaptor<Grid,GFS,V,Projection> gridAdaptor(x._gfs);
752 
753  // Store vector data
754  typename GridAdaptor<Grid,GFS,V,Projection>::MapType transferMap;
755  gridAdaptor.backupData(grid,x._gfs,projection,std::get<I>(x._tuple),transferMap);
756 
757  // Recursively iterate through remaining vectors (and packs). Grid
758  // adaption will be done at the end of recursion.
759  iterateTuple<I + 1, Grid, X, XS...>(grid,x,xs...);
760 
761  // Play back data. Note: At this point the function space was
762  // already updatet.
763  std::get<I>(x._tuple) = V(x._gfs,0.0);
764  gridAdaptor.replayData(grid,x._gfs,projection,std::get<I>(x._tuple),transferMap);
765  }
766 
767  // This gets called after the last pack. After this function call we
768  // have visited every vector of every pack and we will go back through
769  // the recursive function calls.
770  template <typename Grid>
771  void iteratePacks(Grid& grid)
772  {
773  // Adapt the grid
774  grid.adapt();
775  }
776 
777  /* Use template meta programming to iterate over packs at compile time
778  *
779  * In order to adapt our grid and all vectors of all packs we need to
780  * do the following:
781  * - Iterate over all vectors of all packs.
782  * - Store the data from the vectors where things could change.
783  * - Adapt our grid.
784  * - Update function spaces and restore data.
785  *
786  * The key point is that we need the object that stores the data to
787  * replay it. Because of that we can not just iterate over the packs
788  * and within each pack iterate over the vectors but we have to make
789  * one big recursion. Therefore we iterate over the vectors of the
790  * current pack.
791  */
792  template <typename Grid, typename X, typename... XS>
793  void iteratePacks(Grid& grid, X& x, XS&... xs)
794  {
795  iterateTuple(grid,x,xs...);
796  }
797 
798  } // namespace impl
799 #endif // DOXYGEN
800 
812  template <typename GFS, typename... X>
813  impl::GFSWithVectors<GFS,X...> transferSolutions(GFS& gfs, int integrationOrder, X&... x)
814  {
815  impl::GFSWithVectors<GFS,X...> gfsWithVectors(gfs, integrationOrder, x...);
816  return gfsWithVectors;
817  }
818 
829  template <typename Grid, typename... X>
830  void adaptGrid(Grid& grid, X&... x)
831  {
832  // Prepare the grid for refinement
833  grid.preAdapt();
834 
835  // Iterate over packs
836  impl::iteratePacks(grid,x...);
837 
838  // Clean up
839  grid.postAdapt();
840  }
841 
842 
843  template<typename T>
844  void error_fraction(const T& x, typename T::ElementType alpha, typename T::ElementType beta,
845  typename T::ElementType& eta_alpha, typename T::ElementType& eta_beta, int verbose=0)
846  {
847  if (verbose>0)
848  std::cout << "+++ error fraction: alpha=" << alpha << " beta=" << beta << std::endl;
849  const int steps=20; // max number of bisection steps
850  typedef typename T::ElementType NumberType;
851  NumberType total_error = x.one_norm();
852  NumberType max_error = x.infinity_norm();
853  NumberType eta_alpha_left = 0.0;
854  NumberType eta_alpha_right = max_error;
855  NumberType eta_beta_left = 0.0;
856  NumberType eta_beta_right = max_error;
857  for (int j=1; j<=steps; j++)
858  {
859  eta_alpha = 0.5*(eta_alpha_left+eta_alpha_right);
860  eta_beta = 0.5*(eta_beta_left+eta_beta_right);
861  NumberType sum_alpha=0.0;
862  NumberType sum_beta=0.0;
863  unsigned int alpha_count = 0;
864  unsigned int beta_count = 0;
865  for (typename T::const_iterator it = x.begin(),
866  end = x.end();
867  it != end;
868  ++it)
869  {
870  if (*it >=eta_alpha) { sum_alpha += *it; alpha_count++;}
871  if (*it < eta_beta) { sum_beta += *it; beta_count++;}
872  }
873  if (verbose>1)
874  {
875  std::cout << "+++ " << j << " eta_alpha=" << eta_alpha << " alpha_fraction=" << sum_alpha/total_error
876  << " elements: " << alpha_count << " of " << x.N() << std::endl;
877  std::cout << "+++ " << j << " eta_beta=" << eta_beta << " beta_fraction=" << sum_beta/total_error
878  << " elements: " << beta_count << " of " << x.N() << std::endl;
879  }
880  if (std::abs(alpha-sum_alpha/total_error) <= 0.01 && std::abs(beta-sum_beta/total_error) <= 0.01) break;
881  if (sum_alpha>alpha*total_error)
882  eta_alpha_left = eta_alpha;
883  else
884  eta_alpha_right = eta_alpha;
885  if (sum_beta>beta*total_error)
886  eta_beta_right = eta_beta;
887  else
888  eta_beta_left = eta_beta;
889  }
890  if (verbose>0)
891  {
892  std::cout << "+++ refine_threshold=" << eta_alpha
893  << " coarsen_threshold=" << eta_beta << std::endl;
894  }
895  }
896 
897 
898  template<typename T>
899  void element_fraction(const T& x, typename T::ElementType alpha, typename T::ElementType beta,
900  typename T::ElementType& eta_alpha, typename T::ElementType& eta_beta, int verbose=0)
901  {
902  const int steps=20; // max number of bisection steps
903  typedef typename T::ElementType NumberType;
904  NumberType total_error =x.N();
905  NumberType max_error = x.infinity_norm();
906  NumberType eta_alpha_left = 0.0;
907  NumberType eta_alpha_right = max_error;
908  NumberType eta_beta_left = 0.0;
909  NumberType eta_beta_right = max_error;
910  for (int j=1; j<=steps; j++)
911  {
912  eta_alpha = 0.5*(eta_alpha_left+eta_alpha_right);
913  eta_beta = 0.5*(eta_beta_left+eta_beta_right);
914  NumberType sum_alpha=0.0;
915  NumberType sum_beta=0.0;
916  unsigned int alpha_count = 0;
917  unsigned int beta_count = 0;
918 
919  for (typename T::const_iterator it = x.begin(),
920  end = x.end();
921  it != end;
922  ++it)
923  {
924  if (*it>=eta_alpha) { sum_alpha += 1.0; alpha_count++;}
925  if (*it< eta_beta) { sum_beta +=1.0; beta_count++;}
926  }
927  if (verbose>1)
928  {
929  std::cout << j << " eta_alpha=" << eta_alpha << " alpha_fraction=" << sum_alpha/total_error
930  << " elements: " << alpha_count << " of " << x.N() << std::endl;
931  std::cout << j << " eta_beta=" << eta_beta << " beta_fraction=" << sum_beta/total_error
932  << " elements: " << beta_count << " of " << x.N() << std::endl;
933  }
934  if (std::abs(alpha-sum_alpha/total_error) <= 0.01 && std::abs(beta-sum_beta/total_error) <= 0.01) break;
935  if (sum_alpha>alpha*total_error)
936  eta_alpha_left = eta_alpha;
937  else
938  eta_alpha_right = eta_alpha;
939  if (sum_beta>beta*total_error)
940  eta_beta_right = eta_beta;
941  else
942  eta_beta_left = eta_beta;
943  }
944  if (verbose>0)
945  {
946  std::cout << "+++ refine_threshold=" << eta_alpha
947  << " coarsen_threshold=" << eta_beta << std::endl;
948  }
949  }
950 
953  template<typename T>
954  void error_distribution(const T& x, unsigned int bins)
955  {
956  const int steps=30; // max number of bisection steps
957  typedef typename T::ElementType NumberType;
958  NumberType total_error = x.one_norm();
959  NumberType total_elements = x.N();
960  NumberType max_error = x.infinity_norm();
961  std::vector<NumberType> left(bins,0.0);
962  std::vector<NumberType> right(bins,max_error*(1.0+1e-8));
963  std::vector<NumberType> eta(bins);
964  std::vector<NumberType> target(bins);
965  for (unsigned int k=0; k<bins; k++)
966  target[k]= (k+1)/((NumberType)bins);
967  for (int j=1; j<=steps; j++)
968  {
969  for (unsigned int k=0; k<bins; k++)
970  eta[k]= 0.5*(left[k]+right[k]);
971  std::vector<NumberType> sum(bins,0.0);
972  std::vector<int> count(bins,0);
973 
974  for (typename T::const_iterator it = x.begin(),
975  end = x.end();
976  it != end;
977  ++it)
978  {
979  for (unsigned int k=0; k<bins; k++)
980  if (*it<=eta[k])
981  {
982  sum[k] += *it;
983  count[k] += 1;
984  }
985  }
986  // std::cout << std::endl;
987  // std::cout << "// step " << j << std::endl;
988  // for (unsigned int k=0; k<bins; k++)
989  // std::cout << k+1 << " " << count[k] << " " << eta[k] << " " << right[k]-left[k]
990  // << " " << sum[k]/total_error << " " << target[k] << std::endl;
991  for (unsigned int k=0; k<bins; k++)
992  if (sum[k]<=target[k]*total_error)
993  left[k] = eta[k];
994  else
995  right[k] = eta[k];
996  }
997  std::vector<NumberType> sum(bins,0.0);
998  std::vector<int> count(bins,0);
999  for (unsigned int i=0; i<x.N(); i++)
1000  for (unsigned int k=0; k<bins; k++)
1001  if (x[i]<=eta[k])
1002  {
1003  sum[k] += x[i];
1004  count[k] += 1;
1005  }
1006  std::cout << "+++ error distribution" << std::endl;
1007  std::cout << "+++ number of elements: " << x.N() << std::endl;
1008  std::cout << "+++ max element error: " << max_error << std::endl;
1009  std::cout << "+++ total error: " << total_error << std::endl;
1010  std::cout << "+++ bin #elements eta sum/total " << std::endl;
1011  for (unsigned int k=0; k<bins; k++)
1012  std::cout << "+++ " << k+1 << " " << count[k] << " " << eta[k] << " " << sum[k]/total_error << std::endl;
1013  }
1014 
1015  template<typename Grid, typename X>
1016  void mark_grid (Grid &grid, const X& x, typename X::ElementType refine_threshold,
1017  typename X::ElementType coarsen_threshold, int min_level = 0, int max_level = std::numeric_limits<int>::max(), int verbose=0)
1018  {
1019  typedef typename Grid::LeafGridView GV;
1020 
1021  GV gv = grid.leafGridView();
1022 
1023  unsigned int refine_cnt=0;
1024  unsigned int coarsen_cnt=0;
1025 
1026  typedef typename X::GridFunctionSpace GFS;
1027  typedef LocalFunctionSpace<GFS> LFS;
1028  typedef LFSIndexCache<LFS> LFSCache;
1029  typedef typename X::template ConstLocalView<LFSCache> XView;
1030 
1031  LFS lfs(x.gridFunctionSpace());
1032  LFSCache lfs_cache(lfs);
1033  XView x_view(x);
1034 
1035  for(const auto& cell : elements(gv))
1036  {
1037  lfs.bind(cell);
1038  lfs_cache.update();
1039  x_view.bind(lfs_cache);
1040 
1041  if (x_view[0]>=refine_threshold && cell.level() < max_level)
1042  {
1043  grid.mark(1,cell);
1044  refine_cnt++;
1045  }
1046  if (x_view[0]<=coarsen_threshold && cell.level() > min_level)
1047  {
1048  grid.mark(-1,cell);
1049  coarsen_cnt++;
1050  }
1051  x_view.unbind();
1052  }
1053  if (verbose>0)
1054  std::cout << "+++ mark_grid: " << refine_cnt << " marked for refinement, "
1055  << coarsen_cnt << " marked for coarsening" << std::endl;
1056  }
1057 
1058 
1059  template<typename Grid, typename X>
1060  void mark_grid_for_coarsening (Grid &grid, const X& x, typename X::ElementType refine_threshold,
1061  typename X::ElementType coarsen_threshold, int verbose=0)
1062  {
1063  typedef typename Grid::LeafGridView GV;
1064 
1065  GV gv = grid.leafGridView();
1066 
1067  unsigned int coarsen_cnt=0;
1068 
1069  typedef typename X::GridFunctionSpace GFS;
1070  typedef LocalFunctionSpace<GFS> LFS;
1071  typedef LFSIndexCache<LFS> LFSCache;
1072  typedef typename X::template ConstLocalView<LFSCache> XView;
1073 
1074  LFS lfs(x.gridFunctionSpace());
1075  LFSCache lfs_cache(lfs);
1076  XView x_view(x);
1077 
1078  for(const auto& cell : elements(gv))
1079  {
1080  lfs.bind(cell);
1081  lfs_cache.update();
1082  x_view.bind(lfs_cache);
1083 
1084  if (x_view[0]>=refine_threshold)
1085  {
1086  grid.mark(-1,cell);
1087  coarsen_cnt++;
1088  }
1089  if (x_view[0]<=coarsen_threshold)
1090  {
1091  grid.mark(-1,cell);
1092  coarsen_cnt++;
1093  }
1094  }
1095  if (verbose>0)
1096  std::cout << "+++ mark_grid_for_coarsening: "
1097  << coarsen_cnt << " marked for coarsening" << std::endl;
1098  }
1099 
1100 
1102  {
1103  // strategy parameters
1104  double scaling;
1105  double optimistic_factor;
1106  double coarsen_limit;
1107  double balance_limit;
1108  double tol;
1109  double T;
1110  int verbose;
1111  bool no_adapt;
1112  double refine_fraction_while_refinement;
1113  double coarsen_fraction_while_refinement;
1114  double coarsen_fraction_while_coarsening;
1115  double timestep_decrease_factor;
1116  double timestep_increase_factor;
1117 
1118  // results to be reported to the user after evaluating the error
1119  bool accept;
1120  bool adapt_dt;
1121  bool adapt_grid;
1122  double newdt;
1123  double q_s, q_t;
1124 
1125  // state variables
1126  bool have_decreased_time_step;
1127  bool have_refined_grid;
1128 
1129  // the only state variable: accumulated error
1130  double accumulated_estimated_error_squared;
1131  double minenergy_rate;
1132 
1133  public:
1134  TimeAdaptationStrategy (double tol_, double T_, int verbose_=0)
1135  : scaling(16.0), optimistic_factor(1.0), coarsen_limit(0.5), balance_limit(0.33333),
1136  tol(tol_), T(T_), verbose(verbose_), no_adapt(false),
1137  refine_fraction_while_refinement(0.7),
1138  coarsen_fraction_while_refinement(0.2),
1139  coarsen_fraction_while_coarsening(0.2),
1140  timestep_decrease_factor(0.5), timestep_increase_factor(1.5),
1141  accept(false), adapt_dt(false), adapt_grid(false), newdt(1.0),
1142  have_decreased_time_step(false), have_refined_grid(false),
1143  accumulated_estimated_error_squared(0.0),
1144  minenergy_rate(0.0)
1145  {
1146  }
1147 
1149  {
1150  timestep_decrease_factor=s;
1151  }
1152 
1154  {
1155  timestep_increase_factor=s;
1156  }
1157 
1159  {
1160  refine_fraction_while_refinement=s;
1161  }
1162 
1164  {
1165  coarsen_fraction_while_refinement=s;
1166  }
1167 
1169  {
1170  coarsen_fraction_while_coarsening=s;
1171  }
1172 
1173  void setMinEnergyRate (double s)
1174  {
1175  minenergy_rate=s;
1176  }
1177 
1178  void setCoarsenLimit (double s)
1179  {
1180  coarsen_limit=s;
1181  }
1182 
1183  void setBalanceLimit (double s)
1184  {
1185  balance_limit=s;
1186  }
1187 
1188  void setTemporalScaling (double s)
1189  {
1190  scaling=s;
1191  }
1192 
1193  void setOptimisticFactor (double s)
1194  {
1195  optimistic_factor=s;
1196  }
1197 
1199  {
1200  no_adapt = false;
1201  }
1202 
1204  {
1205  no_adapt = true;
1206  }
1207 
1208  bool acceptTimeStep () const
1209  {
1210  return accept;
1211  }
1212 
1213  bool adaptDT () const
1214  {
1215  return adapt_dt;
1216  }
1217 
1218  bool adaptGrid () const
1219  {
1220  return adapt_grid;
1221  }
1222 
1223  double newDT () const
1224  {
1225  return newdt;
1226  }
1227 
1228  double qs () const
1229  {
1230  return q_s;
1231  }
1232 
1233  double qt () const
1234  {
1235  return q_t;
1236  }
1237 
1238  double endT() const
1239  {
1240  return T;
1241  }
1242 
1243  double accumulatedErrorSquared () const
1244  {
1245  return accumulated_estimated_error_squared;
1246  }
1247 
1248  // to be called when new time step is done
1250  {
1251  have_decreased_time_step = false;
1252  have_refined_grid = false;
1253  }
1254 
1255  template<typename GM, typename X>
1256  void evaluate_estimators (GM& grid, double time, double dt, const X& eta_space, const X& eta_time, double energy_timeslab)
1257  {
1258  accept=false;
1259  adapt_dt=false;
1260  adapt_grid=false;
1261  newdt=dt;
1262 
1263  double spatial_error = eta_space.one_norm();
1264  double temporal_error = scaling*eta_time.one_norm();
1265  double sum = spatial_error + temporal_error;
1266  //double allowed = optimistic_factor*(tol*tol-accumulated_estimated_error_squared)*dt/(T-time);
1267  double allowed = tol*tol*(energy_timeslab+minenergy_rate*dt);
1268  q_s = spatial_error/sum;
1269  q_t = temporal_error/sum;
1270 
1271  // print some statistics
1272  if (verbose>0)
1273  std::cout << "+++"
1274  << " q_s=" << q_s
1275  << " q_t=" << q_t
1276  << " sum/allowed=" << sum/allowed
1277  // << " allowed=" << allowed
1278  << " estimated error=" << sqrt(accumulated_estimated_error_squared+sum)
1279  << " energy_rate=" << energy_timeslab/dt
1280  << std::endl;
1281 
1282  // for simplicity: a mode that does no adaptation at all
1283  if (no_adapt)
1284  {
1285  accept = true;
1286  accumulated_estimated_error_squared += sum;
1287  if (verbose>1) std::cout << "+++ no adapt mode" << std::endl;
1288  return;
1289  }
1290 
1291  // the adaptation strategy
1292  if (sum<=allowed)
1293  {
1294  // we will accept this time step
1295  accept = true;
1296  if (verbose>1) std::cout << "+++ accepting time step" << std::endl;
1297  accumulated_estimated_error_squared += sum;
1298 
1299  // check if grid size or time step needs to be adapted
1300  if (sum<coarsen_limit*allowed)
1301  {
1302  // the error is too small, i.e. the computation is inefficient
1303  if (q_t<balance_limit)
1304  {
1305  // spatial error is dominating => increase time step
1306  newdt = timestep_increase_factor*dt;
1307  adapt_dt = true;
1308  if (verbose>1) std::cout << "+++ spatial error dominates: increase time step" << std::endl;
1309  }
1310  else
1311  {
1312  if (q_s>balance_limit)
1313  {
1314  // step sizes balanced: coarsen in time
1315  newdt = timestep_increase_factor*dt;
1316  adapt_dt = true;
1317  if (verbose>1) std::cout << "+++ increasing time step" << std::endl;
1318  }
1319  // coarsen grid in space
1320  double eta_refine, eta_coarsen;
1321  if (verbose>1) std::cout << "+++ mark grid for coarsening" << std::endl;
1322  //error_distribution(eta_space,20);
1323  Dune::PDELab::error_fraction(eta_space,coarsen_fraction_while_coarsening,
1324  coarsen_fraction_while_coarsening,eta_refine,eta_coarsen);
1325  Dune::PDELab::mark_grid_for_coarsening(grid,eta_space,eta_refine,eta_coarsen,verbose);
1326  adapt_grid = true;
1327  }
1328  }
1329  else
1330  {
1331  // modify at least the time step
1332  if (q_t<balance_limit)
1333  {
1334  // spatial error is dominating => increase time step
1335  newdt = timestep_increase_factor*dt;
1336  adapt_dt = true;
1337  if (verbose>1) std::cout << "+++ spatial error dominates: increase time step" << std::endl;
1338  }
1339  }
1340  }
1341  else
1342  {
1343  // error is too large, we need to do something
1344  if (verbose>1) std::cout << "+++ will redo time step" << std::endl;
1345  if (q_t>1-balance_limit)
1346  {
1347  // temporal error is dominating => decrease time step only
1348  newdt = timestep_decrease_factor*dt;
1349  adapt_dt = true;
1350  have_decreased_time_step = true;
1351  if (verbose>1) std::cout << "+++ decreasing time step only" << std::endl;
1352  }
1353  else
1354  {
1355  if (q_t<balance_limit)
1356  {
1357  if (!have_decreased_time_step)
1358  {
1359  // time steps size not balanced (too small)
1360  newdt = timestep_increase_factor*dt;
1361  adapt_dt = true;
1362  if (verbose>1) std::cout << "+++ increasing time step" << std::endl;
1363  }
1364  }
1365  else
1366  {
1367  // step sizes balanced: refine in time as well
1368  newdt = timestep_decrease_factor*dt;
1369  adapt_dt = true;
1370  have_decreased_time_step = true;
1371  if (verbose>1) std::cout << "+++ decreasing time step" << std::endl;
1372  }
1373  // refine grid in space
1374  double eta_refine, eta_coarsen;
1375  if (verbose>1) std::cout << "+++ BINGO mark grid for refinement and coarsening" << std::endl;
1376  //error_distribution(eta_space,20);
1377  Dune::PDELab::error_fraction(eta_space,refine_fraction_while_refinement,
1378  coarsen_fraction_while_refinement,eta_refine,eta_coarsen,0);
1379  Dune::PDELab::mark_grid(grid,eta_space,eta_refine,eta_coarsen,verbose);
1380  adapt_grid = true;
1381  }
1382  }
1383  }
1384  };
1385 
1386 
1387 
1388  } // namespace PDELab
1389 } // namespace Dune
1390 
1391 #endif // DUNE_PDELAB_ADAPTIVITY_ADAPTIVITY_HH
Class for automatic adaptation of the grid.
Definition: adaptivity.hh:512
void setCoarsenFractionWhileRefinement(double s)
Definition: adaptivity.hh:1163
LFSCache _lfs_cache
Definition: adaptivity.hh:482
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:1060
std::array< MassMatrix, TypeTree::TreeInfo< GFS >::leafCount > MassMatrices
Definition: adaptivity.hh:154
const QuadratureRule< DF, dim > & _quadrature_rule
Definition: adaptivity.hh:128
Definition: adaptivity.hh:199
const MassMatrices & inverseMassMatrices(const Element &e)
Calculate the inverse local mass matrix, used in the local L2 projection.
Definition: adaptivity.hh:171
static const int dim
Definition: adaptivity.hh:83
const LocalDOFVector * _u_coarse
Definition: adaptivity.hh:488
Definition: adaptivity.hh:32
void setOptimisticFactor(double s)
Definition: adaptivity.hh:1193
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:1016
Base class for all PDELab exceptions.
Definition: exceptions.hh:17
impl::GFSWithVectors< GFS, X... > transferSolutions(GFS &gfs, int integrationOrder, X &... x)
Pack function space and vectors for grid adaption.
Definition: adaptivity.hh:813
typename EntitySet::Traits::GridView::Grid::LocalIdSet IDSet
Definition: adaptivity.hh:209
void setAdaptationOff()
Definition: adaptivity.hh:1203
void evaluate_estimators(GM &grid, double time, double dt, const X &eta_space, const X &eta_time, double energy_timeslab)
Definition: adaptivity.hh:1256
void adapt_grid(Grid &grid, GFS &gfs, X &x1, int int_order)
adapt a grid, corresponding function space and solution vectors
Definition: adaptivity.hh:613
Iterator extract_lfs_leaf_sizes(const LFS &lfs, Iterator it)
Definition: lfsindexcache.hh:165
static const unsigned int value
Definition: gridfunctionspace/tags.hh:139
Element _element
Definition: adaptivity.hh:349
Element _current
Definition: adaptivity.hh:351
void error_distribution(const T &x, unsigned int bins)
Definition: adaptivity.hh:954
DOFVector::ElementType RF
Definition: adaptivity.hh:379
void operator()(const Element &element)
Definition: adaptivity.hh:273
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:844
const std::string s
Definition: function.hh:1101
Definition: genericdatahandle.hh:622
void leaf(const LeafLFS &leaf_lfs, TreePath treePath)
Definition: adaptivity.hh:420
LeafOffsetCache & _leaf_offset_cache
Definition: adaptivity.hh:489
double accumulatedErrorSquared() const
Definition: adaptivity.hh:1243
double newDT() const
Definition: adaptivity.hh:1223
void setCoarsenFractionWhileCoarsening(double s)
Definition: adaptivity.hh:1168
typename FiniteElement::Traits::LocalBasisType::Traits::RangeType Range
Definition: adaptivity.hh:389
void setTimeStepIncreaseFactor(double s)
Definition: adaptivity.hh:1153
Element::Geometry Geometry
Definition: adaptivity.hh:211
size_type _leaf_index
Definition: adaptivity.hh:490
Geometry _fine_geometry
Definition: adaptivity.hh:411
LocalFunctionSpace< GFS > LFS
Definition: adaptivity.hh:204
void adaptGrid(Grid &grid, X &... x)
Adapt grid and multiple function spaces with corresponding vectors.
Definition: adaptivity.hh:830
LocalCountVector _counts
Definition: adaptivity.hh:493
Projection & _projection
Definition: adaptivity.hh:352
LFSIndexCache< LFS > LFSCache
Definition: adaptivity.hh:372
const Cell & _element
Definition: adaptivity.hh:126
Projection::MassMatrices MassMatrices
Definition: adaptivity.hh:218
TransferMap::mapped_type LocalDOFVector
Definition: adaptivity.hh:214
std::vector< LeafOffsets > _leaf_offset_cache
Definition: adaptivity.hh:70
CountVector::template LocalView< LFSCache > _uc_view
Definition: adaptivity.hh:487
void startTimeStep()
Definition: adaptivity.hh:1249
LocalDOFVector _u_fine
Definition: adaptivity.hh:491
void setBalanceLimit(double s)
Definition: adaptivity.hh:1183
GridAdaptor(const GFSU &gfs)
The constructor.
Definition: adaptivity.hh:529
Projection::MassMatrix MassMatrix
Definition: adaptivity.hh:219
double endT() const
Definition: adaptivity.hh:1238
replay_visitor(const GFS &gfs, DOFVector &u, CountVector &uc, LeafOffsetCache &leaf_offset_cache)
Definition: adaptivity.hh:471
void evaluate(const X &x, Y &y) const
Definition: adaptivity.hh:392
Element _ancestor
Definition: adaptivity.hh:485
coarse_function(const FiniteElement &finite_element, Geometry coarse_geometry, Geometry fine_geometry, const LocalDOFVector &dofs, size_type offset)
Definition: adaptivity.hh:401
std::vector< RF > LocalDOFVector
Definition: adaptivity.hh:380
typename EntitySet::Element Element
Definition: adaptivity.hh:210
std::vector< typename CountVector::ElementType > LocalCountVector
Definition: adaptivity.hh:381
const LocalDOFVector & _dofs
Definition: adaptivity.hh:412
std::unordered_map< ID, std::vector< typename U::ElementType > > MapType
Definition: adaptivity.hh:522
LFSCache _lfs_cache
Definition: adaptivity.hh:347
LocalFunctionSpace< GFS > LFS
Definition: adaptivity.hh:36
Dune::PDELab::LeafOffsetCache< GFS > LeafOffsetCache
Definition: adaptivity.hh:373
typename Element::Geometry Geometry
Definition: adaptivity.hh:378
DOFVector::template ConstLocalView< LFSCache > _u_view
Definition: adaptivity.hh:353
MassMatrices & _mass_matrices
Definition: adaptivity.hh:127
size_type _leaf_index
Definition: adaptivity.hh:129
const std::size_t offset
Definition: localfunctionspace.hh:74
typename EntitySet::Traits::CoordinateField DF
Definition: adaptivity.hh:384
typename EntitySet::Element Element
Definition: adaptivity.hh:377
TimeAdaptationStrategy(double tol_, double T_, int verbose_=0)
Definition: adaptivity.hh:1134
LocalDOFVector _u_tmp
Definition: adaptivity.hh:492
size_type _leaf_index
Definition: adaptivity.hh:358
const GFS & gridFunctionSpace() const
Returns the GridFunctionSpace underlying this LocalFunctionSpace.
Definition: localfunctionspace.hh:264
LeafOffsetCache(const GFS &gfs)
Definition: adaptivity.hh:64
DOFVector::template LocalView< LFSCache > _u_view
Definition: adaptivity.hh:486
bool adaptGrid() const
Definition: adaptivity.hh:1218
double qt() const
Definition: adaptivity.hh:1233
L2Projection(const GFS &gfs, int intorder=2)
The constructor.
Definition: adaptivity.hh:160
const IDSet & _id_set
Definition: adaptivity.hh:483
Geometry _coarse_geometry
Definition: adaptivity.hh:410
typename EntitySet::Traits::GridView::Grid::LocalIdSet IDSet
Definition: adaptivity.hh:376
double qs() const
Definition: adaptivity.hh:1228
typename impl::BackendVectorSelector< GridFunctionSpace, FieldType >::Type Vector
alias of the return type of BackendVectorSelector
Definition: backend/interface.hh:106
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:899
Definition: adaptivity.hh:144
LocalFunctionSpace< GFS > LFS
Definition: adaptivity.hh:371
DOFVector::ElementType RF
Definition: adaptivity.hh:213
typename EntitySet::Traits::CoordinateField DF
Definition: adaptivity.hh:222
LeafOffsetCache & _leaf_offset_cache
Definition: adaptivity.hh:356
typename GFS::Traits::EntitySet EntitySet
Definition: adaptivity.hh:208
Dune::PDELab::LeafOffsetCache< GFS > LeafOffsetCache
Definition: adaptivity.hh:206
LFS _lfs
Definition: adaptivity.hh:481
void operator()(const Element &element, const Element &ancestor, const LocalDOFVector &u_coarse)
Definition: adaptivity.hh:439
typename GFS::Traits::EntitySet EntitySet
Definition: adaptivity.hh:375
LocalDOFVector _u_fine
Definition: adaptivity.hh:359
DynamicMatrix< typename U::ElementType > MassMatrix
Definition: adaptivity.hh:153
Element _element
Definition: adaptivity.hh:484
array< std::size_t, TypeTree::TreeInfo< GFS >::leafCount+1 > LeafOffsets
Definition: adaptivity.hh:40
For backward compatibility – Do not use this!
Definition: adaptivity.hh:27
Traits::IndexContainer::size_type size() const
get current size
Definition: localfunctionspace.hh:206
Definition: adaptivity.hh:1101
std::size_t size_type
Definition: adaptivity.hh:221
Element _ancestor
Definition: adaptivity.hh:350
const Entity & e
Definition: localfunctionspace.hh:111
Definition: adaptivity.hh:366
void replayData(Grid &grid, GFSU &gfsu, Projection &projection, U &u, const MapType &transfer_map)
Definition: adaptivity.hh:554
void setMinEnergyRate(double s)
Definition: adaptivity.hh:1173
void setCoarsenLimit(double s)
Definition: adaptivity.hh:1178
TransferMap & _transfer_map
Definition: adaptivity.hh:354
LocalDOFVector * _u_coarse
Definition: adaptivity.hh:355
size_type _offset
Definition: adaptivity.hh:414
void setAdaptationOn()
Definition: adaptivity.hh:1198
const IDSet & _id_set
Definition: adaptivity.hh:348
void setRefineFractionWhileRefinement(double s)
Definition: adaptivity.hh:1158
LFS _lfs
Definition: adaptivity.hh:346
L2Projection< typename LFS::Traits::GridFunctionSpace, DOFVector > Projection
Definition: adaptivity.hh:217
std::vector< Range > _phi
Definition: adaptivity.hh:413
const LeafOffsets & operator[](GeometryType gt) const
Definition: adaptivity.hh:42
GFS::Traits::GridView::template Codim< 0 >::Entity Cell
Definition: adaptivity.hh:35
LFS _lfs
Definition: adaptivity.hh:69
std::size_t size_type
Definition: adaptivity.hh:383
void leaf(const LFSLeaf &leaf_lfs, TreePath treePath)
Definition: adaptivity.hh:225
bool adaptDT() const
Definition: adaptivity.hh:1213
void backupData(Grid &grid, GFSU &gfsu, Projection &projection, U &u, MapType &transfer_map)
Definition: adaptivity.hh:538
size_type _int_order
Definition: adaptivity.hh:357
const FiniteElement & _finite_element
Definition: adaptivity.hh:409
bool acceptTimeStep() const
Definition: adaptivity.hh:1208
LFSIndexCache< LFS > LFSCache
Definition: adaptivity.hh:205
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 update(const Cell &e)
Definition: adaptivity.hh:50
void setTimeStepDecreaseFactor(double s)
Definition: adaptivity.hh:1148
void setTemporalScaling(double s)
Definition: adaptivity.hh:1188