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