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