4 #ifndef DUNE_PDELAB_ADAPTIVITY_ADAPTIVITY_HH 5 #define DUNE_PDELAB_ADAPTIVITY_ADAPTIVITY_HH 7 #include<dune/common/exceptions.hh> 13 #include<unordered_map> 14 #include<dune/common/dynmatrix.hh> 15 #include<dune/geometry/quadraturerules.hh> 26 #include<dune/grid/io/file/vtk/subsamplingvtkwriter.hh> 32 template<
typename GFS>
36 typedef typename GFS::Traits::GridView::template Codim<0>::Entity
Cell;
41 typedef std::array<std::size_t,TypeTree::TreeInfo<GFS>::leafCount + 1>
LeafOffsets;
47 assert(leaf_offsets.back() > 0);
54 if (leaf_offsets.back() == 0)
59 std::partial_sum(leaf_offsets.begin(),leaf_offsets.end(),leaf_offsets.begin());
61 assert(leaf_offsets.back() ==
_lfs.
size());
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;
88 typedef typename Dune::QuadratureRule<DF,dim>::const_iterator QRIterator;
90 template<
typename GFS,
typename TreePath>
91 void leaf(
const GFS& gfs, TreePath treePath)
93 auto& fem = gfs.finiteElementMap();
95 size_type local_size = fe.localBasis().size();
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));
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)
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)
164 , _inverse_mass_matrices(GlobalGeometryTypeIndex::size(Element::dimension))
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,
186 TypeTree::applyToTree(_gfs,calculate_mass_matrices);
188 return inverse_mass_matrices;
195 std::vector<MassMatrices> _inverse_mass_matrices;
199 template<
typename GFS,
typename DOFVector,
typename TransferMap>
201 :
public TypeTree::TreeVisitor
202 ,
public TypeTree::DynamicTraversal
210 using IDSet =
typename EntitySet::Traits::GridView::Grid::LocalIdSet;
213 static const int dim = Geometry::mydimension;
214 typedef typename DOFVector::ElementType
RF;
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();
233 using Range =
typename LFSLeaf::Traits::GridFunctionSpace::Traits::FiniteElementMap::
234 Traits::FiniteElement::Traits::LocalBasisType::Traits::RangeType;
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();
245 for (
const auto& ip : QuadratureRules<DF,dim>::rule(_current.type(),_int_order))
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 for (size_type i = 0; i < coarse_phi.size(); ++i)
265 for (size_type j = 0; j < inverse_mass_matrix.M(); ++j)
266 x.axpy(inverse_mass_matrix[i][j],coarse_phi[j]);
267 (*_u_coarse)[coarse_offset + i] += factor * (x * val);
280 _u_view.bind(_lfs_cache);
281 _u_coarse = &_transfer_map[_id_set.id(
_element)];
283 _u_view.read(*_u_coarse);
291 while (_ancestor.mightVanish())
294 if (!_ancestor.hasFather())
297 _ancestor = _ancestor.father();
299 _u_coarse = &_transfer_map[_id_set.id(_ancestor)];
301 if (_u_coarse->size() > 0)
304 std::fill(_u_coarse->begin(),_u_coarse->end(),RF(0));
306 for (
const auto& child : descendantElements(_ancestor,max_level))
318 _u_view.bind(_lfs_cache);
319 _u_fine.resize(_lfs_cache.size());
320 _u_view.read(_u_fine);
323 TypeTree::applyToTree(
_lfs,*
this);
330 Projection& projection,
332 LeafOffsetCache& leaf_offset_cache,
333 TransferMap& transfer_map,
334 std::size_t int_order = 2)
337 , _id_set(gfs.gridView().grid().localIdSet())
338 , _projection(projection)
340 , _transfer_map(transfer_map)
343 , _int_order(int_order)
354 typename DOFVector::template ConstLocalView<LFSCache>
_u_view;
366 template<
typename GFS,
typename DOFVector,
typename CountVector>
368 :
public TypeTree::TreeVisitor
369 ,
public TypeTree::DynamicTraversal
377 using IDSet =
typename EntitySet::Traits::GridView::Grid::LocalIdSet;
380 typedef typename DOFVector::ElementType
RF;
385 using DF =
typename EntitySet::Traits::CoordinateField;
387 template<
typename FiniteElement>
390 using Range =
typename FiniteElement::Traits::LocalBasisType::Traits::RangeType;
392 template<
typename X,
typename Y>
395 _phi.resize(_finite_element.localBasis().size());
396 _finite_element.localBasis().evaluateFunction(_coarse_geometry.local(_fine_geometry.global(x)),_phi);
398 for (size_type i = 0; i < _phi.size(); ++i)
399 y.axpy(_dofs[_offset + i],_phi[i]);
403 : _finite_element(finite_element)
404 , _coarse_geometry(coarse_geometry)
405 , _fine_geometry(fine_geometry)
414 mutable std::vector<Range>
_phi;
420 template<
typename LeafLFS,
typename TreePath>
421 void leaf(
const LeafLFS& leaf_lfs, TreePath treePath)
423 using FiniteElement =
typename LeafLFS::Traits::FiniteElementType;
425 auto& fem = leaf_lfs.gridFunctionSpace().finiteElementMap();
432 _u_tmp.resize(fe.localBasis().size());
433 std::fill(_u_tmp.begin(),_u_tmp.end(),RF(0.0));
434 fe.localInterpolation().interpolate(f,_u_tmp);
435 std::copy(_u_tmp.begin(),_u_tmp.end(),_u_fine.begin() + element_offset);
443 _ancestor = ancestor;
444 _u_coarse = &u_coarse;
448 _u_view.bind(_lfs_cache);
451 if (_id_set.id(element) == _id_set.id(ancestor))
454 _u_view.add(*_u_coarse);
458 _u_fine.resize(_lfs_cache.size());
459 std::fill(_u_fine.begin(),_u_fine.end(),RF(0));
461 TypeTree::applyToTree(
_lfs,*
this);
462 _u_view.add(_u_fine);
466 _uc_view.bind(_lfs_cache);
467 _counts.resize(_lfs_cache.size(),1);
468 _uc_view.add(_counts);
472 replay_visitor(
const GFS& gfs, DOFVector& u, CountVector& uc, LeafOffsetCache& leaf_offset_cache)
475 , _id_set(gfs.entitySet().gridView().grid().localIdSet())
487 typename DOFVector::template LocalView<LFSCache>
_u_view;
488 typename CountVector::template LocalView<LFSCache>
_uc_view;
512 template<
class Gr
id,
class GFSU,
class U,
class Projection>
515 typedef typename Grid::LeafGridView LeafGridView;
516 typedef typename LeafGridView::template Codim<0>
517 ::template Partition<Dune::Interior_Partition>::Iterator LeafIterator;
518 typedef typename Grid::template Codim<0>::Entity Element;
519 typedef typename Grid::LocalIdSet IDSet;
520 typedef typename IDSet::IdType ID;
523 typedef std::unordered_map<ID,std::vector<typename U::ElementType> >
MapType;
539 void backupData(Grid& grid, GFSU& gfsu, Projection& projection, U& u, MapType& transfer_map)
546 for(
const auto& cell : elements(gfsu.entitySet(),Partitions::interior))
555 void replayData(Grid& grid, GFSU& gfsu, Projection& projection, U& u,
const MapType& transfer_map)
557 const IDSet& id_set = grid.localIdSet();
560 CountVector uc(gfsu,0);
566 for (
const auto& cell : elements(gfsu.entitySet(),Partitions::interior))
568 Element ancestor = cell;
570 typename MapType::const_iterator map_it;
571 while ((map_it = transfer_map.find(id_set.id(ancestor))) == transfer_map.end())
573 if (!ancestor.hasFather())
575 "transferMap of GridAdaptor didn't contain ancestor of element with id " << id_set.id(ancestor));
576 ancestor = ancestor.father();
579 visitor(cell,ancestor,map_it->second);
583 DOFHandle addHandle1(gfsu,u);
584 gfsu.entitySet().gridView().communicate(addHandle1,
585 Dune::InteriorBorder_InteriorBorder_Interface,Dune::ForwardCommunication);
587 CountHandle addHandle2(gfsu,uc);
588 gfsu.entitySet().gridView().communicate(addHandle2,
589 Dune::InteriorBorder_InteriorBorder_Interface,Dune::ForwardCommunication);
592 typename CountVector::iterator ucit = uc.begin();
593 for (
typename U::iterator uit = u.begin(), uend = u.end(); uit != uend; ++uit, ++ucit)
594 (*uit) /= ((*ucit) > 0 ? (*ucit) : 1.0);
621 template<
class Gr
id,
class GFS,
class X>
625 Projection projection(gfs,int_order);
634 grid_adaptor.
backupData(grid,gfs,projection,x1,transferMap1);
644 grid_adaptor.
replayData(grid,gfs,projection,x1,transferMap1);
661 template<
class Gr
id,
class GFS,
class X>
662 void adapt_grid (Grid& grid, GFS& gfs, X& x1, X& x2,
int int_order)
665 Projection projection(gfs,int_order);
674 grid_adaptor.
backupData(grid,gfs,projection,x1,transferMap1);
676 grid_adaptor.
backupData(grid,gfs,projection,x2,transferMap2);
686 grid_adaptor.
replayData(grid,gfs,projection,x1,transferMap1);
688 grid_adaptor.
replayData(grid,gfs,projection,x2,transferMap2);
698 template <
typename G,
typename... X>
699 struct GFSWithVectors
703 using Tuple = std::tuple<X&...>;
705 GFSWithVectors (GFS& gfs,
int integrationOrder, X&... x) :
707 _integrationOrder(integrationOrder),
712 int _integrationOrder;
717 template <
typename Gr
id>
718 void iteratePacks(Grid& grid);
719 template <
typename Grid,
typename X,
typename... XS>
720 void iteratePacks(Grid& grid, X& x, XS&... xs);
725 template<std::size_t I = 0,
typename Grid,
typename X,
typename... XS>
727 iterateTuple(Grid& grid, X& x, XS&... xs)
730 iteratePacks(grid,xs...);
745 template<std::size_t I = 0,
typename Grid,
typename X,
typename... XS>
747 iterateTuple(Grid& grid, X& x, XS&... xs)
750 using GFS =
typename X::GFS;
751 using Tuple =
typename X::Tuple;
752 using V =
typename std::decay<typename std::tuple_element<I,Tuple>::type>::type;
759 Projection projection(x._gfs,x._integrationOrder);
764 gridAdaptor.
backupData(grid,x._gfs,projection,std::get<I>(x._tuple),transferMap);
768 iterateTuple<I + 1, Grid, X, XS...>(grid,x,xs...);
772 std::get<I>(x._tuple) = V(x._gfs,0.0);
773 gridAdaptor.
replayData(grid,x._gfs,projection,std::get<I>(x._tuple),transferMap);
779 template <
typename Gr
id>
780 void iteratePacks(Grid& grid)
801 template <
typename Grid,
typename X,
typename... XS>
802 void iteratePacks(Grid& grid, X& x, XS&... xs)
804 iterateTuple(grid,x,xs...);
821 template <
typename GFS,
typename... X>
824 impl::GFSWithVectors<GFS,X...> gfsWithVectors(gfs, integrationOrder, x...);
825 return gfsWithVectors;
838 template <
typename Grid,
typename... X>
845 impl::iteratePacks(grid,x...);
853 void error_fraction(
const T& x,
typename T::ElementType alpha,
typename T::ElementType beta,
854 typename T::ElementType& eta_alpha,
typename T::ElementType& eta_beta,
int verbose=0)
857 std::cout <<
"+++ error fraction: alpha=" << alpha <<
" beta=" << beta << std::endl;
859 typedef typename T::ElementType NumberType;
860 NumberType total_error = x.one_norm();
861 NumberType max_error = x.infinity_norm();
862 NumberType eta_alpha_left = 0.0;
863 NumberType eta_alpha_right = max_error;
864 NumberType eta_beta_left = 0.0;
865 NumberType eta_beta_right = max_error;
866 for (
int j=1; j<=steps; j++)
868 eta_alpha = 0.5*(eta_alpha_left+eta_alpha_right);
869 eta_beta = 0.5*(eta_beta_left+eta_beta_right);
870 NumberType sum_alpha=0.0;
871 NumberType sum_beta=0.0;
872 unsigned int alpha_count = 0;
873 unsigned int beta_count = 0;
874 for (
const auto& error : x)
876 if (error >=eta_alpha) { sum_alpha += error; alpha_count++;}
877 if (error < eta_beta) { sum_beta += error; beta_count++;}
881 std::cout <<
"+++ " << j <<
" eta_alpha=" << eta_alpha <<
" alpha_fraction=" << sum_alpha/total_error
882 <<
" elements: " << alpha_count <<
" of " << x.N() << std::endl;
883 std::cout <<
"+++ " << j <<
" eta_beta=" << eta_beta <<
" beta_fraction=" << sum_beta/total_error
884 <<
" elements: " << beta_count <<
" of " << x.N() << std::endl;
886 if (std::abs(alpha-sum_alpha/total_error) <= 0.01 && std::abs(beta-sum_beta/total_error) <= 0.01)
break;
887 if (sum_alpha>alpha*total_error)
888 eta_alpha_left = eta_alpha;
890 eta_alpha_right = eta_alpha;
891 if (sum_beta>beta*total_error)
892 eta_beta_right = eta_beta;
894 eta_beta_left = eta_beta;
898 std::cout <<
"+++ refine_threshold=" << eta_alpha
899 <<
" coarsen_threshold=" << eta_beta << std::endl;
905 void element_fraction(
const T& x,
typename T::ElementType alpha,
typename T::ElementType beta,
906 typename T::ElementType& eta_alpha,
typename T::ElementType& eta_beta,
int verbose=0)
909 typedef typename T::ElementType NumberType;
910 NumberType total_error =x.N();
911 NumberType max_error = x.infinity_norm();
912 NumberType eta_alpha_left = 0.0;
913 NumberType eta_alpha_right = max_error;
914 NumberType eta_beta_left = 0.0;
915 NumberType eta_beta_right = max_error;
916 for (
int j=1; j<=steps; j++)
918 eta_alpha = 0.5*(eta_alpha_left+eta_alpha_right);
919 eta_beta = 0.5*(eta_beta_left+eta_beta_right);
920 NumberType sum_alpha=0.0;
921 NumberType sum_beta=0.0;
922 unsigned int alpha_count = 0;
923 unsigned int beta_count = 0;
925 for (
const auto& error : x)
927 if (error>=eta_alpha) { sum_alpha += 1.0; alpha_count++;}
928 if (error< eta_beta) { sum_beta +=1.0; beta_count++;}
932 std::cout << j <<
" eta_alpha=" << eta_alpha <<
" alpha_fraction=" << sum_alpha/total_error
933 <<
" elements: " << alpha_count <<
" of " << x.N() << std::endl;
934 std::cout << j <<
" eta_beta=" << eta_beta <<
" beta_fraction=" << sum_beta/total_error
935 <<
" elements: " << beta_count <<
" of " << x.N() << std::endl;
937 if (std::abs(alpha-sum_alpha/total_error) <= 0.01 && std::abs(beta-sum_beta/total_error) <= 0.01)
break;
938 if (sum_alpha>alpha*total_error)
939 eta_alpha_left = eta_alpha;
941 eta_alpha_right = eta_alpha;
942 if (sum_beta>beta*total_error)
943 eta_beta_right = eta_beta;
945 eta_beta_left = eta_beta;
949 std::cout <<
"+++ refine_threshold=" << eta_alpha
950 <<
" coarsen_threshold=" << eta_beta << std::endl;
960 typedef typename T::ElementType NumberType;
961 NumberType total_error = x.one_norm();
962 NumberType total_elements = x.N();
963 NumberType max_error = x.infinity_norm();
964 std::vector<NumberType> left(bins,0.0);
965 std::vector<NumberType> right(bins,max_error*(1.0+1
e-8));
966 std::vector<NumberType> eta(bins);
967 std::vector<NumberType> target(bins);
968 for (
unsigned int k=0; k<bins; k++)
969 target[k]= (k+1)/((NumberType)bins);
970 for (
int j=1; j<=steps; j++)
972 for (
unsigned int k=0; k<bins; k++)
973 eta[k]= 0.5*(left[k]+right[k]);
974 std::vector<NumberType> sum(bins,0.0);
975 std::vector<int> count(bins,0);
977 for (
typename T::const_iterator it = x.begin(),
982 for (
unsigned int k=0; k<bins; k++)
994 for (
unsigned int k=0; k<bins; k++)
995 if (sum[k]<=target[k]*total_error)
1000 std::vector<NumberType> sum(bins,0.0);
1001 std::vector<int> count(bins,0);
1002 for (
unsigned int i=0; i<x.N(); i++)
1003 for (
unsigned int k=0; k<bins; k++)
1009 std::cout <<
"+++ error distribution" << std::endl;
1010 std::cout <<
"+++ number of elements: " << x.N() << std::endl;
1011 std::cout <<
"+++ max element error: " << max_error << std::endl;
1012 std::cout <<
"+++ total error: " << total_error << std::endl;
1013 std::cout <<
"+++ bin #elements eta sum/total " << std::endl;
1014 for (
unsigned int k=0; k<bins; k++)
1015 std::cout <<
"+++ " << k+1 <<
" " << count[k] <<
" " << eta[k] <<
" " << sum[k]/total_error << std::endl;
1018 template<
typename Gr
id,
typename X>
1019 void mark_grid (Grid &grid,
const X& x,
typename X::ElementType refine_threshold,
1020 typename X::ElementType coarsen_threshold,
int min_level = 0,
int max_level = std::numeric_limits<int>::max(),
int verbose=0)
1022 typedef typename Grid::LeafGridView GV;
1024 GV gv = grid.leafGridView();
1026 unsigned int refine_cnt=0;
1027 unsigned int coarsen_cnt=0;
1029 typedef typename X::GridFunctionSpace GFS;
1032 typedef typename X::template ConstLocalView<LFSCache> XView;
1034 LFS lfs(x.gridFunctionSpace());
1035 LFSCache lfs_cache(lfs);
1038 for(
const auto& cell : elements(gv))
1042 x_view.bind(lfs_cache);
1044 if (x_view[0]>=refine_threshold && cell.level() < max_level)
1049 if (x_view[0]<=coarsen_threshold && cell.level() > min_level)
1057 std::cout <<
"+++ mark_grid: " << refine_cnt <<
" marked for refinement, " 1058 << coarsen_cnt <<
" marked for coarsening" << std::endl;
1062 template<
typename Gr
id,
typename X>
1064 typename X::ElementType coarsen_threshold,
int verbose=0)
1066 typedef typename Grid::LeafGridView GV;
1068 GV gv = grid.leafGridView();
1070 unsigned int coarsen_cnt=0;
1072 typedef typename X::GridFunctionSpace GFS;
1075 typedef typename X::template ConstLocalView<LFSCache> XView;
1077 LFS lfs(x.gridFunctionSpace());
1078 LFSCache lfs_cache(lfs);
1081 for(
const auto& cell : elements(gv))
1085 x_view.bind(lfs_cache);
1087 if (x_view[0]>=refine_threshold)
1092 if (x_view[0]<=coarsen_threshold)
1099 std::cout <<
"+++ mark_grid_for_coarsening: " 1100 << coarsen_cnt <<
" marked for coarsening" << std::endl;
1108 double optimistic_factor;
1109 double coarsen_limit;
1110 double balance_limit;
1115 double refine_fraction_while_refinement;
1116 double coarsen_fraction_while_refinement;
1117 double coarsen_fraction_while_coarsening;
1118 double timestep_decrease_factor;
1119 double timestep_increase_factor;
1129 bool have_decreased_time_step;
1130 bool have_refined_grid;
1133 double accumulated_estimated_error_squared;
1134 double minenergy_rate;
1138 : scaling(16.0), optimistic_factor(1.0), coarsen_limit(0.5), balance_limit(0.33333),
1139 tol(tol_), T(T_), verbose(verbose_), no_adapt(false),
1140 refine_fraction_while_refinement(0.7),
1141 coarsen_fraction_while_refinement(0.2),
1142 coarsen_fraction_while_coarsening(0.2),
1143 timestep_decrease_factor(0.5), timestep_increase_factor(1.5),
1144 accept(false), adapt_dt(false), adapt_grid(false), newdt(1.0),
1145 have_decreased_time_step(false), have_refined_grid(false),
1146 accumulated_estimated_error_squared(0.0),
1153 timestep_decrease_factor=
s;
1158 timestep_increase_factor=
s;
1163 refine_fraction_while_refinement=
s;
1168 coarsen_fraction_while_refinement=
s;
1173 coarsen_fraction_while_coarsening=
s;
1198 optimistic_factor=
s;
1248 return accumulated_estimated_error_squared;
1254 have_decreased_time_step =
false;
1255 have_refined_grid =
false;
1258 template<
typename GM,
typename X>
1259 void evaluate_estimators (GM& grid,
double time,
double dt,
const X& eta_space,
const X& eta_time,
double energy_timeslab)
1266 double spatial_error = eta_space.one_norm();
1267 double temporal_error = scaling*eta_time.one_norm();
1268 double sum = spatial_error + temporal_error;
1270 double allowed = tol*tol*(energy_timeslab+minenergy_rate*dt);
1271 q_s = spatial_error/sum;
1272 q_t = temporal_error/sum;
1279 <<
" sum/allowed=" << sum/allowed
1281 <<
" estimated error=" << sqrt(accumulated_estimated_error_squared+sum)
1282 <<
" energy_rate=" << energy_timeslab/dt
1289 accumulated_estimated_error_squared += sum;
1290 if (verbose>1) std::cout <<
"+++ no adapt mode" << std::endl;
1299 if (verbose>1) std::cout <<
"+++ accepting time step" << std::endl;
1300 accumulated_estimated_error_squared += sum;
1303 if (sum<coarsen_limit*allowed)
1306 if (q_t<balance_limit)
1309 newdt = timestep_increase_factor*dt;
1311 if (verbose>1) std::cout <<
"+++ spatial error dominates: increase time step" << std::endl;
1315 if (q_s>balance_limit)
1318 newdt = timestep_increase_factor*dt;
1320 if (verbose>1) std::cout <<
"+++ increasing time step" << std::endl;
1323 double eta_refine, eta_coarsen;
1324 if (verbose>1) std::cout <<
"+++ mark grid for coarsening" << std::endl;
1327 coarsen_fraction_while_coarsening,eta_refine,eta_coarsen);
1335 if (q_t<balance_limit)
1338 newdt = timestep_increase_factor*dt;
1340 if (verbose>1) std::cout <<
"+++ spatial error dominates: increase time step" << std::endl;
1347 if (verbose>1) std::cout <<
"+++ will redo time step" << std::endl;
1348 if (q_t>1-balance_limit)
1351 newdt = timestep_decrease_factor*dt;
1353 have_decreased_time_step =
true;
1354 if (verbose>1) std::cout <<
"+++ decreasing time step only" << std::endl;
1358 if (q_t<balance_limit)
1360 if (!have_decreased_time_step)
1363 newdt = timestep_increase_factor*dt;
1365 if (verbose>1) std::cout <<
"+++ increasing time step" << std::endl;
1371 newdt = timestep_decrease_factor*dt;
1373 have_decreased_time_step =
true;
1374 if (verbose>1) std::cout <<
"+++ decreasing time step" << std::endl;
1377 double eta_refine, eta_coarsen;
1378 if (verbose>1) std::cout <<
"+++ BINGO mark grid for refinement and coarsening" << std::endl;
1381 coarsen_fraction_while_refinement,eta_refine,eta_coarsen,0);
1394 #endif // DUNE_PDELAB_ADAPTIVITY_ADAPTIVITY_HH std::array< std::size_t, TypeTree::TreeInfo< GFS >::leafCount+1 > LeafOffsets
Definition: adaptivity.hh:41
const LeafOffsets & operator[](GeometryType gt) const
Definition: adaptivity.hh:43
Element _ancestor
Definition: adaptivity.hh:351
GFS::Traits::GridView::template Codim< 0 >::Entity Cell
Definition: adaptivity.hh:36
DOFVector::ElementType RF
Definition: adaptivity.hh:380
const QuadratureRule< DF, dim > & _quadrature_rule
Definition: adaptivity.hh:129
L2Projection(const GFS &gfs, int intorder=2)
The constructor.
Definition: adaptivity.hh:161
static const int dim
Definition: adaptivity.hh:84
bool adaptDT() const
Definition: adaptivity.hh:1216
void backupData(Grid &grid, GFSU &gfsu, Projection &projection, U &u, MapType &transfer_map)
Definition: adaptivity.hh:539
const Entity & e
Definition: localfunctionspace.hh:120
void leaf(const LeafLFS &leaf_lfs, TreePath treePath)
Definition: adaptivity.hh:421
typename EntitySet::Element Element
Definition: adaptivity.hh:211
backup_visitor(const GFS &gfs, Projection &projection, const DOFVector &u, LeafOffsetCache &leaf_offset_cache, TransferMap &transfer_map, std::size_t int_order=2)
Definition: adaptivity.hh:329
LocalDOFVector * _u_coarse
Definition: adaptivity.hh:356
typename FiniteElement::Traits::LocalBasisType::Traits::RangeType Range
Definition: adaptivity.hh:390
Class for automatic adaptation of the grid.
Definition: adaptivity.hh:513
void setCoarsenFractionWhileRefinement(double s)
Definition: adaptivity.hh:1166
LocalFunctionSpace< GFS > LFS
Definition: adaptivity.hh:37
void adapt_grid(Grid &grid, GFS &gfs, X &x1, int int_order)
adapt a grid, corresponding function space and solution vectors
Definition: adaptivity.hh:622
Definition: adaptivity.hh:200
Geometry _fine_geometry
Definition: adaptivity.hh:412
typename EntitySet::Traits::CoordinateField DF
Definition: adaptivity.hh:223
typename impl::BackendVectorSelector< GridFunctionSpace, FieldType >::Type Vector
alias of the return type of BackendVectorSelector
Definition: backend/interface.hh:106
const LocalDOFVector * _u_coarse
Definition: adaptivity.hh:489
Definition: adaptivity.hh:33
LFSIndexCache< LFS > LFSCache
Definition: adaptivity.hh:373
void setAdaptationOff()
Definition: adaptivity.hh:1206
CountVector::template LocalView< LFSCache > _uc_view
Definition: adaptivity.hh:488
void startTimeStep()
Definition: adaptivity.hh:1252
LocalDOFVector _u_fine
Definition: adaptivity.hh:492
void setBalanceLimit(double s)
Definition: adaptivity.hh:1186
Definition: adaptivity.hh:388
double endT() const
Definition: adaptivity.hh:1241
replay_visitor(const GFS &gfs, DOFVector &u, CountVector &uc, LeafOffsetCache &leaf_offset_cache)
Definition: adaptivity.hh:472
void leaf(const LFSLeaf &leaf_lfs, TreePath treePath)
Definition: adaptivity.hh:226
LeafOffsetCache & _leaf_offset_cache
Definition: adaptivity.hh:490
LFSIndexCache< LFS > LFSCache
Definition: adaptivity.hh:206
const Cell & _element
Definition: adaptivity.hh:127
coarse_function(const FiniteElement &finite_element, Geometry coarse_geometry, Geometry fine_geometry, const LocalDOFVector &dofs, size_type offset)
Definition: adaptivity.hh:402
double accumulatedErrorSquared() const
Definition: adaptivity.hh:1246
double newDT() const
Definition: adaptivity.hh:1226
void setCoarsenFractionWhileCoarsening(double s)
Definition: adaptivity.hh:1171
void setTimeStepIncreaseFactor(double s)
Definition: adaptivity.hh:1156
Dune::PDELab::LeafOffsetCache< GFS > LeafOffsetCache
Definition: adaptivity.hh:374
size_type _leaf_index
Definition: adaptivity.hh:491
static const unsigned int value
Definition: gridfunctionspace/tags.hh:139
DOFVector::template ConstLocalView< LFSCache > _u_view
Definition: adaptivity.hh:354
LeafOffsetCache & _leaf_offset_cache
Definition: adaptivity.hh:357
LocalCountVector _counts
Definition: adaptivity.hh:494
TimeAdaptationStrategy(double tol_, double T_, int verbose_=0)
Definition: adaptivity.hh:1137
LocalDOFVector _u_fine
Definition: adaptivity.hh:360
L2Projection< typename LFS::Traits::GridFunctionSpace, DOFVector > Projection
Definition: adaptivity.hh:218
Element _element
Definition: adaptivity.hh:350
std::size_t size_type
Definition: adaptivity.hh:222
const std::size_t offset
Definition: localfunctionspace.hh:74
GridAdaptor(const GFSU &gfs)
The constructor.
Definition: adaptivity.hh:530
DOFVector::template LocalView< LFSCache > _u_view
Definition: adaptivity.hh:487
MassMatrices & _mass_matrices
Definition: adaptivity.hh:128
size_type _leaf_index
Definition: adaptivity.hh:130
bool adaptGrid() const
Definition: adaptivity.hh:1221
LFS _lfs
Definition: adaptivity.hh:70
void evaluate(const X &x, Y &y) const
Definition: adaptivity.hh:393
Element _ancestor
Definition: adaptivity.hh:486
const IDSet & _id_set
Definition: adaptivity.hh:484
TransferMap & _transfer_map
Definition: adaptivity.hh:355
Geometry _coarse_geometry
Definition: adaptivity.hh:411
std::vector< RF > LocalDOFVector
Definition: adaptivity.hh:381
std::vector< typename CountVector::ElementType > LocalCountVector
Definition: adaptivity.hh:382
const LocalDOFVector & _dofs
Definition: adaptivity.hh:413
std::unordered_map< ID, std::vector< typename U::ElementType > > MapType
Definition: adaptivity.hh:523
std::array< MassMatrix, TypeTree::TreeInfo< GFS >::leafCount > MassMatrices
Definition: adaptivity.hh:155
typename Element::Geometry Geometry
Definition: adaptivity.hh:379
void element_fraction(const T &x, typename T::ElementType alpha, typename T::ElementType beta, typename T::ElementType &eta_alpha, typename T::ElementType &eta_beta, int verbose=0)
Definition: adaptivity.hh:905
For backward compatibility – Do not use this!
Definition: adaptivity.hh:28
const IDSet & _id_set
Definition: adaptivity.hh:349
typename GFS::Traits::EntitySet EntitySet
Definition: adaptivity.hh:209
const MassMatrices & inverseMassMatrices(const Element &e)
Calculate the inverse local mass matrix, used in the local L2 projection.
Definition: adaptivity.hh:172
typename EntitySet::Traits::CoordinateField DF
Definition: adaptivity.hh:385
Dune::PDELab::LeafOffsetCache< GFS > LeafOffsetCache
Definition: adaptivity.hh:207
typename EntitySet::Element Element
Definition: adaptivity.hh:378
LocalDOFVector _u_tmp
Definition: adaptivity.hh:493
size_type _leaf_index
Definition: adaptivity.hh:359
void operator()(const Element &element, const Element &ancestor, const LocalDOFVector &u_coarse)
Definition: adaptivity.hh:440
Projection::MassMatrices MassMatrices
Definition: adaptivity.hh:219
Element _current
Definition: adaptivity.hh:352
TransferMap::mapped_type LocalDOFVector
Definition: adaptivity.hh:215
Element _element
Definition: adaptivity.hh:485
double qt() const
Definition: adaptivity.hh:1236
void replayData(Grid &grid, GFSU &gfsu, Projection &projection, U &u, const MapType &transfer_map)
Definition: adaptivity.hh:555
size_type _int_order
Definition: adaptivity.hh:358
void setMinEnergyRate(double s)
Definition: adaptivity.hh:1176
void setCoarsenLimit(double s)
Definition: adaptivity.hh:1181
typename EntitySet::Traits::GridView::Grid::LocalIdSet IDSet
Definition: adaptivity.hh:377
void update(const Cell &e)
Definition: adaptivity.hh:51
double qs() const
Definition: adaptivity.hh:1231
Element::Geometry Geometry
Definition: adaptivity.hh:212
size_type _offset
Definition: adaptivity.hh:415
Definition: adaptivity.hh:145
LocalFunctionSpace< GFS > LFS
Definition: adaptivity.hh:372
LocalFunctionSpace< GFS > LFS
Definition: adaptivity.hh:205
void adaptGrid(Grid &grid, X &... x)
Adapt grid and multiple function spaces with corresponding vectors.
Definition: adaptivity.hh:839
LFS _lfs
Definition: adaptivity.hh:482
void setRefineFractionWhileRefinement(double s)
Definition: adaptivity.hh:1161
LFS _lfs
Definition: adaptivity.hh:347
typename EntitySet::Traits::GridView::Grid::LocalIdSet IDSet
Definition: adaptivity.hh:210
Projection & _projection
Definition: adaptivity.hh:353
typename GFS::Traits::EntitySet EntitySet
Definition: adaptivity.hh:376
Iterator extract_lfs_leaf_sizes(const LFS &lfs, Iterator it)
Definition: lfsindexcache.hh:165
LeafOffsetCache(const GFS &gfs)
Definition: adaptivity.hh:65
Definition: adaptivity.hh:1104
std::vector< Range > _phi
Definition: adaptivity.hh:414
Projection::MassMatrix MassMatrix
Definition: adaptivity.hh:220
void operator()(const Element &element)
Definition: adaptivity.hh:274
const GFS & gridFunctionSpace() const
Returns the GridFunctionSpace underlying this LocalFunctionSpace.
Definition: localfunctionspace.hh:278
Definition: adaptivity.hh:367
std::size_t size_type
Definition: adaptivity.hh:384
Base class for all PDELab exceptions.
Definition: exceptions.hh:17
const FiniteElement & _finite_element
Definition: adaptivity.hh:410
bool acceptTimeStep() const
Definition: adaptivity.hh:1211
Definition: genericdatahandle.hh:665
void mark_grid_for_coarsening(Grid &grid, const X &x, typename X::ElementType refine_threshold, typename X::ElementType coarsen_threshold, int verbose=0)
Definition: adaptivity.hh:1063
LFSCache _lfs_cache
Definition: adaptivity.hh:348
void setTimeStepDecreaseFactor(double s)
Definition: adaptivity.hh:1151
std::size_t index
Definition: interpolate.hh:118
void setTemporalScaling(double s)
Definition: adaptivity.hh:1191
void error_fraction(const T &x, typename T::ElementType alpha, typename T::ElementType beta, typename T::ElementType &eta_alpha, typename T::ElementType &eta_beta, int verbose=0)
Definition: adaptivity.hh:853
void mark_grid(Grid &grid, const X &x, typename X::ElementType refine_threshold, typename X::ElementType coarsen_threshold, int min_level=0, int max_level=std::numeric_limits< int >::max(), int verbose=0)
Definition: adaptivity.hh:1019
LFSCache _lfs_cache
Definition: adaptivity.hh:483
DOFVector::ElementType RF
Definition: adaptivity.hh:214
impl::GFSWithVectors< GFS, X... > transferSolutions(GFS &gfs, int integrationOrder, X &... x)
Pack function space and vectors for grid adaption.
Definition: adaptivity.hh:822
void setAdaptationOn()
Definition: adaptivity.hh:1201
const std::string s
Definition: function.hh:830
void setOptimisticFactor(double s)
Definition: adaptivity.hh:1196
void error_distribution(const T &x, unsigned int bins)
Definition: adaptivity.hh:957
std::vector< LeafOffsets > _leaf_offset_cache
Definition: adaptivity.hh:71
DynamicMatrix< typename U::ElementType > MassMatrix
Definition: adaptivity.hh:154
void evaluate_estimators(GM &grid, double time, double dt, const X &eta_space, const X &eta_time, double energy_timeslab)
Definition: adaptivity.hh:1259
Traits::IndexContainer::size_type size() const
get current size
Definition: localfunctionspace.hh:220