DUNE PDELab (git)

adaptivity.hh
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
8
9#include<array>
10#include<limits>
11#include<vector>
12#include<map>
13#include<unordered_map>
16#include<dune/pdelab/gridfunctionspace/genericdatahandle.hh>
17#include<dune/pdelab/gridfunctionspace/localfunctionspace.hh>
18
19#include<dune/pdelab/common/function.hh>
20// for InterpolateBackendStandard
21#include<dune/pdelab/gridfunctionspace/interpolate.hh>
22// for intersectionoperator
23#include<dune/pdelab/localoperator/defaultimp.hh>
24#include<dune/pdelab/localoperator/flags.hh>
25
27
28namespace Dune {
29 namespace PDELab {
30
31
32 template<typename GFS>
33 struct LeafOffsetCache
34 {
35
36 typedef typename GFS::Traits::GridView::template Codim<0>::Entity Cell;
37 typedef LocalFunctionSpace<GFS> LFS;
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
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>
200 struct backup_visitor
201 : public TypeTree::TreeVisitor
203 {
204
205 typedef LocalFunctionSpace<GFS> LFS;
206 typedef LFSIndexCache<LFS> LFSCache;
207 typedef Dune::PDELab::LeafOffsetCache<GFS> LeafOffsetCache;
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
218 typedef L2Projection<typename LFS::Traits::GridFunctionSpace,DOFVector> Projection;
219 typedef typename Projection::MassMatrices MassMatrices;
220 typedef typename Projection::MassMatrix MassMatrix;
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);
282 _u_coarse = &_transfer_map[_id_set.id(_element)];
283 _u_coarse->resize(_lfs.size());
284 _u_view.read(*_u_coarse);
285 _u_view.unbind();
286
287 _leaf_offset_cache.update(_element);
288
289 size_type max_level = _lfs.gridFunctionSpace().gridView().grid().maxLevel();
290
291 _ancestor = _element;
292 while (_ancestor.mightVanish())
293 {
294 // work around UG bug!
295 if (!_ancestor.hasFather())
296 break;
297
298 _ancestor = _ancestor.father();
299
300 _u_coarse = &_transfer_map[_id_set.id(_ancestor)];
301 // don't project more than once
302 if (_u_coarse->size() > 0)
303 continue;
304
305 _leaf_offset_cache.update(_ancestor);
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);
319 _leaf_offset_cache.update(_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
350 LFS _lfs;
351 LFSCache _lfs_cache;
352 const IDSet& _id_set;
353 Element _element;
354 Element _ancestor;
355 Element _current;
356 Projection& _projection;
357 typename DOFVector::template ConstLocalView<LFSCache> _u_view;
358 TransferMap& _transfer_map;
359 LocalDOFVector* _u_coarse;
360 LeafOffsetCache& _leaf_offset_cache;
361 size_type _int_order;
362 size_type _leaf_index;
363 LocalDOFVector _u_fine;
364
365 };
366
367
368
369 template<typename GFS, typename DOFVector, typename CountVector>
370 struct replay_visitor
371 : public TypeTree::TreeVisitor
372 , public TypeTree::DynamicTraversal
373 {
374
375 typedef LocalFunctionSpace<GFS> LFS;
376 typedef LFSIndexCache<LFS> LFSCache;
377 typedef Dune::PDELab::LeafOffsetCache<GFS> LeafOffsetCache;
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>
391 struct coarse_function
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;
421 Geometry _coarse_geometry;
422 Geometry _fine_geometry;
423 const LocalDOFVector& _dofs;
424 mutable std::vector<Range> _phi;
425 size_type _offset;
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);
456 _leaf_offset_cache.update(_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
492 LFS _lfs;
493 LFSCache _lfs_cache;
494 const IDSet& _id_set;
495 Element _element;
496 Element _ancestor;
497 typename DOFVector::template LocalView<LFSCache> _u_view;
498 typename CountVector::template LocalView<LFSCache> _uc_view;
499 const LocalDOFVector* _u_coarse;
500 LeafOffsetCache& _leaf_offset_cache;
501 size_type _leaf_index;
502 LocalDOFVector _u_fine;
503 LocalDOFVector _u_tmp;
504 LocalCountVector _counts;
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())
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,
596 typedef Dune::PDELab::AddDataHandle<GFSU,CountVector> CountHandle;
597 CountHandle addHandle2(gfsu,uc);
598 gfsu.entitySet().gridView().communicate(addHandle2,
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>
736 inline typename std::enable_if<I == std::tuple_size<typename X::Tuple>::value, void>::type
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>
756 inline typename std::enable_if<I < std::tuple_size<typename X::Tuple>::value, void>::type
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
1114 class TimeAdaptationStrategy
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
1161 void setTimeStepDecreaseFactor (double s)
1162 {
1163 timestep_decrease_factor=s;
1164 }
1165
1166 void setTimeStepIncreaseFactor (double s)
1167 {
1168 timestep_increase_factor=s;
1169 }
1170
1171 void setRefineFractionWhileRefinement (double s)
1172 {
1173 refine_fraction_while_refinement=s;
1174 }
1175
1176 void setCoarsenFractionWhileRefinement (double s)
1177 {
1178 coarsen_fraction_while_refinement=s;
1179 }
1180
1181 void setCoarsenFractionWhileCoarsening (double s)
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
1211 void setAdaptationOn ()
1212 {
1213 no_adapt = false;
1214 }
1215
1216 void setAdaptationOff ()
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
1262 void startTimeStep ()
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 using std::sqrt;
1272 accept=false;
1273 adapt_dt=false;
1274 adapt_grid=false;
1275 newdt=dt;
1276
1277 double spatial_error = eta_space.one_norm();
1278 double temporal_error = scaling*eta_time.one_norm();
1279 double sum = spatial_error + temporal_error;
1280 //double allowed = optimistic_factor*(tol*tol-accumulated_estimated_error_squared)*dt/(T-time);
1281 double allowed = tol*tol*(energy_timeslab+minenergy_rate*dt);
1282 q_s = spatial_error/sum;
1283 q_t = temporal_error/sum;
1284
1285 // print some statistics
1286 if (verbose>0)
1287 std::cout << "+++"
1288 << " q_s=" << q_s
1289 << " q_t=" << q_t
1290 << " sum/allowed=" << sum/allowed
1291 // << " allowed=" << allowed
1292 << " estimated error=" << sqrt(accumulated_estimated_error_squared+sum)
1293 << " energy_rate=" << energy_timeslab/dt
1294 << std::endl;
1295
1296 // for simplicity: a mode that does no adaptation at all
1297 if (no_adapt)
1298 {
1299 accept = true;
1300 accumulated_estimated_error_squared += sum;
1301 if (verbose>1) std::cout << "+++ no adapt mode" << std::endl;
1302 return;
1303 }
1304
1305 // the adaptation strategy
1306 if (sum<=allowed)
1307 {
1308 // we will accept this time step
1309 accept = true;
1310 if (verbose>1) std::cout << "+++ accepting time step" << std::endl;
1311 accumulated_estimated_error_squared += sum;
1312
1313 // check if grid size or time step needs to be adapted
1314 if (sum<coarsen_limit*allowed)
1315 {
1316 // the error is too small, i.e. the computation is inefficient
1317 if (q_t<balance_limit)
1318 {
1319 // spatial error is dominating => increase time step
1320 newdt = timestep_increase_factor*dt;
1321 adapt_dt = true;
1322 if (verbose>1) std::cout << "+++ spatial error dominates: increase time step" << std::endl;
1323 }
1324 else
1325 {
1326 if (q_s>balance_limit)
1327 {
1328 // step sizes balanced: coarsen in time
1329 newdt = timestep_increase_factor*dt;
1330 adapt_dt = true;
1331 if (verbose>1) std::cout << "+++ increasing time step" << std::endl;
1332 }
1333 // coarsen grid in space
1334 double eta_refine, eta_coarsen;
1335 if (verbose>1) std::cout << "+++ mark grid for coarsening" << std::endl;
1336 //error_distribution(eta_space,20);
1337 Dune::PDELab::error_fraction(eta_space,coarsen_fraction_while_coarsening,
1338 coarsen_fraction_while_coarsening,eta_refine,eta_coarsen);
1339 Dune::PDELab::mark_grid_for_coarsening(grid,eta_space,eta_refine,eta_coarsen,verbose);
1340 adapt_grid = true;
1341 }
1342 }
1343 else
1344 {
1345 // modify at least the time step
1346 if (q_t<balance_limit)
1347 {
1348 // spatial error is dominating => increase time step
1349 newdt = timestep_increase_factor*dt;
1350 adapt_dt = true;
1351 if (verbose>1) std::cout << "+++ spatial error dominates: increase time step" << std::endl;
1352 }
1353 }
1354 }
1355 else
1356 {
1357 // error is too large, we need to do something
1358 if (verbose>1) std::cout << "+++ will redo time step" << std::endl;
1359 if (q_t>1-balance_limit)
1360 {
1361 // temporal error is dominating => decrease time step only
1362 newdt = timestep_decrease_factor*dt;
1363 adapt_dt = true;
1364 have_decreased_time_step = true;
1365 if (verbose>1) std::cout << "+++ decreasing time step only" << std::endl;
1366 }
1367 else
1368 {
1369 if (q_t<balance_limit)
1370 {
1371 if (!have_decreased_time_step)
1372 {
1373 // time steps size not balanced (too small)
1374 newdt = timestep_increase_factor*dt;
1375 adapt_dt = true;
1376 if (verbose>1) std::cout << "+++ increasing time step" << std::endl;
1377 }
1378 }
1379 else
1380 {
1381 // step sizes balanced: refine in time as well
1382 newdt = timestep_decrease_factor*dt;
1383 adapt_dt = true;
1384 have_decreased_time_step = true;
1385 if (verbose>1) std::cout << "+++ decreasing time step" << std::endl;
1386 }
1387 // refine grid in space
1388 double eta_refine, eta_coarsen;
1389 if (verbose>1) std::cout << "+++ BINGO mark grid for refinement and coarsening" << std::endl;
1390 //error_distribution(eta_space,20);
1391 Dune::PDELab::error_fraction(eta_space,refine_fraction_while_refinement,
1392 coarsen_fraction_while_refinement,eta_refine,eta_coarsen,0);
1393 Dune::PDELab::mark_grid(grid,eta_space,eta_refine,eta_coarsen,verbose);
1394 adapt_grid = true;
1395 }
1396 }
1397 }
1398 };
1399
1400
1401
1402 } // namespace PDELab
1403} // namespace Dune
1404
1405#endif // DUNE_PDELAB_ADAPTIVITY_ADAPTIVITY_HH
Construct a matrix with a dynamic size.
Definition: dynmatrix.hh:61
Base class for Dune-Exceptions.
Definition: exceptions.hh:96
Wrapper class for geometries.
Definition: geometry.hh:71
static constexpr int mydimension
geometry dimension
Definition: geometry.hh:94
Compute indices for geometry types, taking the dimension into account.
Definition: typeindex.hh:90
static constexpr std::size_t index(const GeometryType &gt)
Compute the index for the given geometry type over all dimensions.
Definition: typeindex.hh:138
Grid abstract base class.
Definition: grid.hh:375
GridFamily::Traits::LocalIdSet LocalIdSet
A type that is a model of Dune::IdSet which provides a unique and persistent numbering for all entiti...
Definition: grid.hh:509
const LocalIdSet & localIdSet() const
return const reference to the grids local id set
Definition: grid.hh:612
GridFamily::Traits::LeafGridView LeafGridView
type of view for leaf grid
Definition: grid.hh:400
Class for automatic adaptation of the grid.
Definition: adaptivity.hh:524
GridAdaptor(const GFSU &gfs)
The constructor.
Definition: adaptivity.hh:540
Definition: adaptivity.hh:146
const MassMatrices & inverseMassMatrices(const Element &e)
Calculate the inverse local mass matrix, used in the local L2 projection.
Definition: adaptivity.hh:172
L2Projection(const GFS &gfs, int intorder=2)
The constructor.
Definition: adaptivity.hh:161
Abstract base class for quadrature rules.
Definition: quadraturerules.hh:214
static const QuadratureRule & rule(const GeometryType &t, int p, QuadratureType::Enum qt=QuadratureType::GaussLegendre)
select the appropriate QuadratureRule for GeometryType t and order p
Definition: quadraturerules.hh:326
A few common exception classes.
This file implements a dense matrix with dynamic numbers of rows and columns.
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
bool gt(const T &first, const T &second, typename EpsilonType< T >::Type epsilon)
test if first greater than second
Definition: float_cmp.cc:158
@ ForwardCommunication
communicate as given in InterfaceType
Definition: gridenums.hh:171
@ InteriorBorder_InteriorBorder_Interface
send/receive interior and border entities
Definition: gridenums.hh:87
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:484
constexpr auto treePath(const T &... t)
Constructs a new HybridTreePath from the given indices.
Definition: treepath.hh:326
void applyToTree(Tree &&tree, Visitor &&visitor)
Apply visitor to TypeTree.
Definition: traversal.hh:239
ImplementationDefined child(Node &&node, Indices... indices)
Extracts the child of a node given by a sequence of compile-time and run-time indices.
Definition: childextraction.hh:128
constexpr Interior interior
PartitionSet for the interior partition.
Definition: partitionset.hh:271
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
Traits class containing decayed types.
Definition: adaptivity.hh:398
Mixin base class for visitors that only need a dynamic TreePath during traversal.
Definition: visitor.hh:426
Convenience base class for visiting the entire tree.
Definition: visitor.hh:435
Provides subsampled file i/o for the visualization toolkit.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 12, 23:30, 2024)