4#ifndef DUNE_PDELAB_ADAPTIVITY_ADAPTIVITY_HH
5#define DUNE_PDELAB_ADAPTIVITY_ADAPTIVITY_HH
13#include<unordered_map>
16#include<dune/pdelab/gridfunctionspace/genericdatahandle.hh>
17#include<dune/pdelab/gridfunctionspace/localfunctionspace.hh>
19#include<dune/pdelab/common/function.hh>
21#include<dune/pdelab/gridfunctionspace/interpolate.hh>
23#include<dune/pdelab/localoperator/defaultimp.hh>
24#include<dune/pdelab/localoperator/flags.hh>
32 template<
typename GFS>
33 struct LeafOffsetCache
36 typedef typename GFS::Traits::GridView::template Codim<0>::Entity Cell;
37 typedef LocalFunctionSpace<GFS> LFS;
41 typedef std::array<std::size_t,TypeTree::TreeInfo<GFS>::leafCount + 1> LeafOffsets;
43 const LeafOffsets& operator[](GeometryType
gt)
const
47 assert(leaf_offsets.back() > 0);
51 void update(
const Cell& e)
54 if (leaf_offsets.back() == 0)
57 extract_lfs_leaf_sizes(_lfs,leaf_offsets.begin()+1);
59 std::partial_sum(leaf_offsets.begin(),leaf_offsets.end(),leaf_offsets.begin());
61 assert(leaf_offsets.back() == _lfs.size());
65 explicit LeafOffsetCache(
const GFS& gfs)
67 , _leaf_offset_cache(GlobalGeometryTypeIndex::size(Cell::dimension))
71 std::vector<LeafOffsets> _leaf_offset_cache;
78 template<
typename MassMatrices,
typename Cell>
79 struct inverse_mass_matrix_calculator
80 :
public TypeTree::TreeVisitor
81 ,
public TypeTree::DynamicTraversal
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;
90 template<
typename GFS,
typename TreePath>
91 void leaf(
const GFS& gfs, TreePath
treePath)
93 auto& fem = gfs.finiteElementMap();
94 auto& fe = fem.find(_element);
95 size_type local_size = fe.localBasis().size();
97 MassMatrix& mass_matrix = _mass_matrices[_leaf_index];
98 mass_matrix.resize(local_size,local_size);
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));
105 for (
const auto& ip : _quadrature_rule)
107 fe.localBasis().evaluateFunction(ip.position(),phi);
108 const DF factor = ip.weight();
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;
115 mass_matrix.invert();
120 inverse_mass_matrix_calculator(MassMatrices& mass_matrices,
const Cell& element, size_type intorder)
122 , _mass_matrices(mass_matrices)
123 , _quadrature_rule(QuadratureRules<DF,dim>::rule(element.type(),intorder))
127 const Cell& _element;
128 MassMatrices& _mass_matrices;
129 const QuadratureRule<DF,dim>& _quadrature_rule;
130 size_type _leaf_index;
144 template<
class GFS,
class U>
147 using EntitySet =
typename GFS::Traits::EntitySet;
148 using Element =
typename EntitySet::Element;
150 typedef typename U::ElementType DF;
155 typedef std::array<MassMatrix,TypeTree::TreeInfo<GFS>::leafCount> MassMatrices;
163 , _intorder(intorder)
174 auto gt = e.geometry().type();
177 if (inverse_mass_matrices[0].N() > 0)
178 return inverse_mass_matrices;
180 inverse_mass_matrix_calculator<MassMatrices,Element> calculate_mass_matrices(
181 inverse_mass_matrices,
188 return inverse_mass_matrices;
195 std::vector<MassMatrices> _inverse_mass_matrices;
199 template<
typename GFS,
typename DOFVector,
typename TransferMap>
200 struct backup_visitor
205 typedef LocalFunctionSpace<GFS> LFS;
206 typedef LFSIndexCache<LFS> LFSCache;
207 typedef Dune::PDELab::LeafOffsetCache<GFS> LeafOffsetCache;
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;
214 typedef typename DOFVector::ElementType RF;
215 typedef typename TransferMap::mapped_type LocalDOFVector;
218 typedef L2Projection<typename LFS::Traits::GridFunctionSpace,DOFVector> Projection;
219 typedef typename Projection::MassMatrices MassMatrices;
220 typedef typename Projection::MassMatrix MassMatrix;
222 typedef std::size_t size_type;
223 using DF =
typename EntitySet::Traits::CoordinateField;
225 template<
typename LFSLeaf,
typename TreePath>
226 void leaf(
const LFSLeaf& leaf_lfs, TreePath
treePath)
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];
233 using Range =
typename LFSLeaf::Traits::GridFunctionSpace::Traits::FiniteElementMap::
234 Traits::FiniteElement::Traits::LocalBasisType::Traits::RangeType;
236 auto& inverse_mass_matrix = _projection.inverseMassMatrices(_ancestor)[_leaf_index];
238 auto coarse_phi = std::vector<Range>{};
239 auto fine_phi = std::vector<Range>{};
241 auto fine_geometry = _current.geometry();
242 auto coarse_geometry = _ancestor.geometry();
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);
256 auto val = Range{0.0};
257 for (size_type i = 0; i < fine_phi.size(); ++i)
259 val.axpy(_u_fine[fine_offset + i],fine_phi[i]);
262 assert(inverse_mass_matrix.M()==coarse_phi.size());
263 for (size_type i = 0; i < coarse_phi.size(); ++i)
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);
275 void operator()(
const Element& element)
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);
287 _leaf_offset_cache.update(_element);
289 size_type max_level = _lfs.gridFunctionSpace().gridView().grid().maxLevel();
291 _ancestor = _element;
292 while (_ancestor.mightVanish())
295 if (!_ancestor.hasFather())
298 _ancestor = _ancestor.father();
300 _u_coarse = &_transfer_map[_id_set.id(_ancestor)];
302 if (_u_coarse->size() > 0)
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));
309 for (
const auto&
child : descendantElements(_ancestor,max_level))
319 _leaf_offset_cache.update(_current);
321 _u_view.bind(_lfs_cache);
322 _u_fine.resize(_lfs_cache.size());
323 _u_view.read(_u_fine);
332 backup_visitor(
const GFS& gfs,
333 Projection& projection,
335 LeafOffsetCache& leaf_offset_cache,
336 TransferMap& transfer_map,
337 std::size_t int_order = 2)
340 , _id_set(gfs.gridView().grid().localIdSet())
341 , _projection(projection)
343 , _transfer_map(transfer_map)
345 , _leaf_offset_cache(leaf_offset_cache)
346 , _int_order(int_order)
352 const IDSet& _id_set;
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;
369 template<
typename GFS,
typename DOFVector,
typename CountVector>
370 struct replay_visitor
371 :
public TypeTree::TreeVisitor
372 ,
public TypeTree::DynamicTraversal
375 typedef LocalFunctionSpace<GFS> LFS;
376 typedef LFSIndexCache<LFS> LFSCache;
377 typedef Dune::PDELab::LeafOffsetCache<GFS> LeafOffsetCache;
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;
387 typedef std::size_t size_type;
388 using DF =
typename EntitySet::Traits::CoordinateField;
390 template<
typename FiniteElement>
391 struct coarse_function
393 using Range =
typename FiniteElement::Traits::LocalBasisType::Traits::RangeType;
394 using Domain =
typename FiniteElement::Traits::LocalBasisType::Traits::DomainType;
399 using RangeType = std::decay_t<Range>;
402 Range operator()(
const Domain& x)
const
404 _phi.resize(_finite_element.localBasis().size());
405 _finite_element.localBasis().evaluateFunction(_coarse_geometry.local(_fine_geometry.global(x)),_phi);
407 for (size_type i = 0; i < _phi.size(); ++i)
408 y.axpy(_dofs[_offset + i],_phi[i]);
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)
420 const FiniteElement& _finite_element;
423 const LocalDOFVector& _dofs;
424 mutable std::vector<Range> _phi;
430 template<
typename LeafLFS,
typename TreePath>
431 void leaf(
const LeafLFS& leaf_lfs, TreePath
treePath)
433 using FiniteElement =
typename LeafLFS::Traits::FiniteElementType;
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];
439 coarse_function<FiniteElement> f(fem.find(_ancestor),_ancestor.geometry(),_element.geometry(),*_u_coarse,ancestor_offset);
440 auto& fe = fem.find(_element);
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);
450 void operator()(
const Element& element,
const Element& ancestor,
const LocalDOFVector& u_coarse)
453 _ancestor = ancestor;
454 _u_coarse = &u_coarse;
456 _leaf_offset_cache.update(_element);
458 _u_view.bind(_lfs_cache);
461 if (_id_set.id(element) == _id_set.id(ancestor))
464 _u_view.add(*_u_coarse);
468 _u_fine.resize(_lfs_cache.size());
469 std::fill(_u_fine.begin(),_u_fine.end(),RF(0));
472 _u_view.add(_u_fine);
476 _uc_view.bind(_lfs_cache);
477 _counts.resize(_lfs_cache.size(),1);
478 _uc_view.add(_counts);
482 replay_visitor(
const GFS& gfs, DOFVector& u, CountVector& uc, LeafOffsetCache& leaf_offset_cache)
485 , _id_set(gfs.entitySet().gridView().grid().localIdSet())
488 , _leaf_offset_cache(leaf_offset_cache)
494 const IDSet& _id_set;
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;
522 template<
class Gr
id,
class GFSU,
class U,
class Projection>
526 typedef typename LeafGridView::template
Codim<0>
527 ::template Partition<Dune::Interior_Partition>::Iterator LeafIterator;
530 typedef typename IDSet::IdType ID;
533 typedef std::unordered_map<ID,std::vector<typename U::ElementType> > MapType;
541 : _leaf_offset_cache(gfs)
549 void backupData(
Grid& grid, GFSU& gfsu, Projection& projection, U& u, MapType& transfer_map)
551 typedef backup_visitor<GFSU,U,MapType> Visitor;
553 Visitor visitor(gfsu,projection,u,_leaf_offset_cache,transfer_map);
565 void replayData(
Grid& grid, GFSU& gfsu, Projection& projection, U& u,
const MapType& transfer_map)
569 using CountVector = Backend::Vector<GFSU,int>;
570 CountVector uc(gfsu,0);
572 typedef replay_visitor<GFSU,U,CountVector> Visitor;
573 Visitor visitor(gfsu,u,uc,_leaf_offset_cache);
578 Element ancestor = cell;
580 typename MapType::const_iterator map_it;
581 while ((map_it = transfer_map.find(id_set.id(ancestor))) == transfer_map.end())
583 if (!ancestor.hasFather())
585 "transferMap of GridAdaptor didn't contain ancestor of element with id " << id_set.id(ancestor));
586 ancestor = ancestor.father();
589 visitor(cell,ancestor,map_it->second);
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,
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);
609 LeafOffsetCache<GFSU> _leaf_offset_cache;
631 template<
class Gr
id,
class GFS,
class X>
632 void adapt_grid (Grid& grid, GFS& gfs, X& x1,
int int_order)
634 typedef L2Projection<GFS,X> Projection;
635 Projection projection(gfs,int_order);
637 GridAdaptor<Grid,GFS,X,Projection> grid_adaptor(gfs);
643 typename GridAdaptor<Grid,GFS,X,Projection>::MapType transferMap1;
644 grid_adaptor.backupData(grid,gfs,projection,x1,transferMap1);
654 grid_adaptor.replayData(grid,gfs,projection,x1,transferMap1);
671 template<
class Gr
id,
class GFS,
class X>
672 void adapt_grid (Grid& grid, GFS& gfs, X& x1, X& x2,
int int_order)
674 typedef L2Projection<GFS,X> Projection;
675 Projection projection(gfs,int_order);
677 GridAdaptor<Grid,GFS,X,Projection> grid_adaptor(gfs);
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);
696 grid_adaptor.replayData(grid,gfs,projection,x1,transferMap1);
698 grid_adaptor.replayData(grid,gfs,projection,x2,transferMap2);
708 template <
typename G,
typename... X>
709 struct GFSWithVectors
713 using Tuple = std::tuple<X&...>;
715 GFSWithVectors (GFS& gfs,
int integrationOrder, X&... x) :
717 _integrationOrder(integrationOrder),
722 int _integrationOrder;
727 template <
typename Gr
id>
728 void iteratePacks(Grid& grid);
729 template <
typename Grid,
typename X,
typename... XS>
730 void iteratePacks(Grid& grid, X& x, XS&... xs);
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)
740 iteratePacks(grid,xs...);
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)
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;
769 Projection projection(x._gfs,x._integrationOrder);
770 GridAdaptor<Grid,GFS,V,Projection> gridAdaptor(x._gfs);
773 typename GridAdaptor<Grid,GFS,V,Projection>::MapType transferMap;
774 gridAdaptor.backupData(grid,x._gfs,projection,std::get<I>(x._tuple),transferMap);
778 iterateTuple<I + 1, Grid, X, XS...>(grid,x,xs...);
782 std::get<I>(x._tuple) = V(x._gfs,0.0);
783 gridAdaptor.replayData(grid,x._gfs,projection,std::get<I>(x._tuple),transferMap);
789 template <
typename Gr
id>
790 void iteratePacks(Grid& grid)
811 template <
typename Grid,
typename X,
typename... XS>
812 void iteratePacks(Grid& grid, X& x, XS&... xs)
814 iterateTuple(grid,x,xs...);
831 template <
typename GFS,
typename... X>
832 impl::GFSWithVectors<GFS,X...> transferSolutions(GFS& gfs,
int integrationOrder, X&... x)
834 impl::GFSWithVectors<GFS,X...> gfsWithVectors(gfs, integrationOrder, x...);
835 return gfsWithVectors;
848 template <
typename Grid,
typename... X>
849 void adaptGrid(Grid& grid, X&... x)
855 impl::iteratePacks(grid,x...);
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)
867 std::cout <<
"+++ error fraction: alpha=" << alpha <<
" beta=" << beta << std::endl;
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++)
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)
886 if (error >=eta_alpha) { sum_alpha += error; alpha_count++;}
887 if (error < eta_beta) { sum_beta += error; beta_count++;}
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;
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;
900 eta_alpha_right = eta_alpha;
901 if (sum_beta>beta*total_error)
902 eta_beta_right = eta_beta;
904 eta_beta_left = eta_beta;
908 std::cout <<
"+++ refine_threshold=" << eta_alpha
909 <<
" coarsen_threshold=" << eta_beta << std::endl;
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)
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++)
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;
935 for (
const auto& error : x)
937 if (error>=eta_alpha) { sum_alpha += 1.0; alpha_count++;}
938 if (error< eta_beta) { sum_beta +=1.0; beta_count++;}
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;
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;
951 eta_alpha_right = eta_alpha;
952 if (sum_beta>beta*total_error)
953 eta_beta_right = eta_beta;
955 eta_beta_left = eta_beta;
959 std::cout <<
"+++ refine_threshold=" << eta_alpha
960 <<
" coarsen_threshold=" << eta_beta << std::endl;
967 void error_distribution(
const T& x,
unsigned int bins)
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++)
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);
987 for (
typename T::const_iterator it = x.begin(),
992 for (
unsigned int k=0; k<bins; k++)
1004 for (
unsigned int k=0; k<bins; k++)
1005 if (sum[k]<=target[k]*total_error)
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++)
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;
1028 template<
typename Gr
id,
typename X>
1029 void mark_grid (Grid &grid,
const X& x,
typename X::ElementType refine_threshold,
1034 GV gv = grid.leafGridView();
1036 unsigned int refine_cnt=0;
1037 unsigned int coarsen_cnt=0;
1039 typedef typename X::GridFunctionSpace GFS;
1040 typedef LocalFunctionSpace<GFS> LFS;
1041 typedef LFSIndexCache<LFS> LFSCache;
1042 typedef typename X::template ConstLocalView<LFSCache> XView;
1044 LFS lfs(x.gridFunctionSpace());
1045 LFSCache lfs_cache(lfs);
1048 for(
const auto& cell : elements(gv))
1052 x_view.bind(lfs_cache);
1054 if (x_view[0]>=refine_threshold && cell.level() < max_level)
1059 if (x_view[0]<=coarsen_threshold && cell.level() > min_level)
1067 std::cout <<
"+++ mark_grid: " << refine_cnt <<
" marked for refinement, "
1068 << coarsen_cnt <<
" marked for coarsening" << std::endl;
1072 template<
typename Gr
id,
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)
1078 GV gv = grid.leafGridView();
1080 unsigned int coarsen_cnt=0;
1082 typedef typename X::GridFunctionSpace GFS;
1083 typedef LocalFunctionSpace<GFS> LFS;
1084 typedef LFSIndexCache<LFS> LFSCache;
1085 typedef typename X::template ConstLocalView<LFSCache> XView;
1087 LFS lfs(x.gridFunctionSpace());
1088 LFSCache lfs_cache(lfs);
1091 for(
const auto& cell : elements(gv))
1095 x_view.bind(lfs_cache);
1097 if (x_view[0]>=refine_threshold)
1102 if (x_view[0]<=coarsen_threshold)
1109 std::cout <<
"+++ mark_grid_for_coarsening: "
1110 << coarsen_cnt <<
" marked for coarsening" << std::endl;
1114 class TimeAdaptationStrategy
1118 double optimistic_factor;
1119 double coarsen_limit;
1120 double balance_limit;
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;
1139 bool have_decreased_time_step;
1140 bool have_refined_grid;
1143 double accumulated_estimated_error_squared;
1144 double minenergy_rate;
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),
1161 void setTimeStepDecreaseFactor (
double s)
1163 timestep_decrease_factor=s;
1166 void setTimeStepIncreaseFactor (
double s)
1168 timestep_increase_factor=s;
1171 void setRefineFractionWhileRefinement (
double s)
1173 refine_fraction_while_refinement=s;
1176 void setCoarsenFractionWhileRefinement (
double s)
1178 coarsen_fraction_while_refinement=s;
1181 void setCoarsenFractionWhileCoarsening (
double s)
1183 coarsen_fraction_while_coarsening=s;
1186 void setMinEnergyRate (
double s)
1191 void setCoarsenLimit (
double s)
1196 void setBalanceLimit (
double s)
1201 void setTemporalScaling (
double s)
1206 void setOptimisticFactor (
double s)
1208 optimistic_factor=s;
1211 void setAdaptationOn ()
1216 void setAdaptationOff ()
1221 bool acceptTimeStep ()
const
1226 bool adaptDT ()
const
1231 bool adaptGrid ()
const
1236 double newDT ()
const
1256 double accumulatedErrorSquared ()
const
1258 return accumulated_estimated_error_squared;
1262 void startTimeStep ()
1264 have_decreased_time_step =
false;
1265 have_refined_grid =
false;
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)
1277 double spatial_error = eta_space.one_norm();
1278 double temporal_error = scaling*eta_time.one_norm();
1279 double sum = spatial_error + temporal_error;
1281 double allowed = tol*tol*(energy_timeslab+minenergy_rate*dt);
1282 q_s = spatial_error/sum;
1283 q_t = temporal_error/sum;
1290 <<
" sum/allowed=" << sum/allowed
1292 <<
" estimated error=" << sqrt(accumulated_estimated_error_squared+sum)
1293 <<
" energy_rate=" << energy_timeslab/dt
1300 accumulated_estimated_error_squared += sum;
1301 if (verbose>1) std::cout <<
"+++ no adapt mode" << std::endl;
1310 if (verbose>1) std::cout <<
"+++ accepting time step" << std::endl;
1311 accumulated_estimated_error_squared += sum;
1314 if (sum<coarsen_limit*allowed)
1317 if (q_t<balance_limit)
1320 newdt = timestep_increase_factor*dt;
1322 if (verbose>1) std::cout <<
"+++ spatial error dominates: increase time step" << std::endl;
1326 if (q_s>balance_limit)
1329 newdt = timestep_increase_factor*dt;
1331 if (verbose>1) std::cout <<
"+++ increasing time step" << std::endl;
1334 double eta_refine, eta_coarsen;
1335 if (verbose>1) std::cout <<
"+++ mark grid for coarsening" << std::endl;
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);
1346 if (q_t<balance_limit)
1349 newdt = timestep_increase_factor*dt;
1351 if (verbose>1) std::cout <<
"+++ spatial error dominates: increase time step" << std::endl;
1358 if (verbose>1) std::cout <<
"+++ will redo time step" << std::endl;
1359 if (q_t>1-balance_limit)
1362 newdt = timestep_decrease_factor*dt;
1364 have_decreased_time_step =
true;
1365 if (verbose>1) std::cout <<
"+++ decreasing time step only" << std::endl;
1369 if (q_t<balance_limit)
1371 if (!have_decreased_time_step)
1374 newdt = timestep_increase_factor*dt;
1376 if (verbose>1) std::cout <<
"+++ increasing time step" << std::endl;
1382 newdt = timestep_decrease_factor*dt;
1384 have_decreased_time_step =
true;
1385 if (verbose>1) std::cout <<
"+++ decreasing time step" << std::endl;
1388 double eta_refine, eta_coarsen;
1389 if (verbose>1) std::cout <<
"+++ BINGO mark grid for refinement and coarsening" << std::endl;
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);
Construct a matrix with a dynamic size.
Definition: dynmatrix.hh:59
Base class for Dune-Exceptions.
Definition: exceptions.hh:94
Wrapper class for geometries.
Definition: geometry.hh:67
@ mydimension
Definition: geometry.hh:90
Compute indices for geometry types, taking the dimension into account.
Definition: typeindex.hh:88
static constexpr std::size_t index(const GeometryType >)
Compute the index for the given geometry type over all dimensions.
Definition: typeindex.hh:136
Grid abstract base class.
Definition: grid.hh:372
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:512
const LocalIdSet & localIdSet() const
return const reference to the grids local id set
Definition: grid.hh:615
GridFamily::Traits::LeafGridView LeafGridView
type of view for leaf grid
Definition: grid.hh:403
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:152
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:280
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:216
bool gt(const T &first, const T &second, typename EpsilonType< T >::Type epsilon)
test if first greater than second
Definition: float_cmp.cc:156
auto max(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::max()
Definition: defaults.hh:79
constexpr HybridTreePath< T... > treePath(const T &... t)
Constructs a new HybridTreePath from the given indices.
Definition: treepath.hh:188
void applyToTree(Tree &&tree, Visitor &&visitor)
Apply visitor to TypeTree.
Definition: traversal.hh:237
constexpr Interior interior
PartitionSet for the interior partition.
Definition: partitionset.hh:270
Dune namespace.
Definition: alignedallocator.hh:11
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:424
Convenience base class for visiting the entire tree.
Definition: visitor.hh:433
Provides subsampled file i/o for the visualization toolkit.