4 #ifndef DUNE_PDELAB_ADAPTIVITY_HH 5 #define DUNE_PDELAB_ADAPTIVITY_HH 7 #include<dune/common/exceptions.hh> 12 #include<unordered_map> 13 #include<dune/common/dynmatrix.hh> 14 #include<dune/geometry/quadraturerules.hh> 25 #include<dune/grid/io/file/vtk/subsamplingvtkwriter.hh> 31 template<
typename GFS>
35 typedef typename GFS::Traits::GridView::template Codim<0>::Entity
Cell;
40 typedef array<std::size_t,TypeTree::TreeInfo<GFS>::leafCount + 1>
LeafOffsets;
44 const LeafOffsets& leaf_offsets =
_leaf_offset_cache[GlobalGeometryTypeIndex::index(gt)];
46 assert(leaf_offsets.back() > 0);
52 LeafOffsets& leaf_offsets =
_leaf_offset_cache[GlobalGeometryTypeIndex::index(e.type())];
53 if (leaf_offsets.back() == 0)
58 std::partial_sum(leaf_offsets.begin(),leaf_offsets.end(),leaf_offsets.begin());
60 assert(leaf_offsets.back() ==
_lfs.
size());
77 template<
typename MassMatrices,
typename Cell>
78 struct inverse_mass_matrix_calculator
79 :
public TypeTree::TreeVisitor
80 ,
public TypeTree::DynamicTraversal
83 static const int dim = Cell::Geometry::mydimension;
84 typedef std::size_t size_type;
85 typedef typename MassMatrices::value_type MassMatrix;
86 typedef typename MassMatrix::field_type DF;
87 typedef typename Dune::QuadratureRule<DF,dim>::const_iterator QRIterator;
89 template<
typename GFS,
typename TreePath>
90 void leaf(
const GFS& gfs, TreePath treePath)
92 auto& fem = gfs.finiteElementMap();
94 size_type local_size = fe.localBasis().size();
97 mass_matrix.resize(local_size,local_size);
99 using Range =
typename GFS::Traits::FiniteElementMap::Traits::
100 FiniteElement::Traits::LocalBasisType::Traits::RangeType;
101 std::vector<Range> phi;
102 phi.resize(std::max(phi.size(),local_size));
106 fe.localBasis().evaluateFunction(ip.position(),phi);
107 const DF factor = ip.weight();
109 for (size_type i = 0; i < local_size; ++i)
110 for (size_type j = 0; j < local_size; ++j)
111 mass_matrix[i][j] += phi[i] * phi[j] * factor;
114 mass_matrix.invert();
119 inverse_mass_matrix_calculator(MassMatrices& mass_matrices,
const Cell& element, size_type intorder)
143 template<
class GFS,
class U>
146 using EntitySet =
typename GFS::Traits::EntitySet;
147 using Element =
typename EntitySet::Element;
149 typedef typename U::ElementType DF;
154 typedef std::array<MassMatrix,TypeTree::TreeInfo<GFS>::leafCount>
MassMatrices;
162 , _intorder(intorder)
163 , _inverse_mass_matrices(GlobalGeometryTypeIndex::size(Element::dimension))
173 auto gt = e.geometry().type();
174 auto& inverse_mass_matrices = _inverse_mass_matrices[GlobalGeometryTypeIndex::index(gt)];
176 if (inverse_mass_matrices[0].N() > 0)
177 return inverse_mass_matrices;
179 inverse_mass_matrix_calculator<MassMatrices,Element> calculate_mass_matrices(
180 inverse_mass_matrices,
185 TypeTree::applyToTree(_gfs,calculate_mass_matrices);
187 return inverse_mass_matrices;
194 std::vector<MassMatrices> _inverse_mass_matrices;
198 template<
typename GFS,
typename DOFVector,
typename TransferMap>
200 :
public TypeTree::TreeVisitor
201 ,
public TypeTree::DynamicTraversal
209 using IDSet =
typename EntitySet::Traits::GridView::Grid::LocalIdSet;
212 static const int dim = Geometry::mydimension;
213 typedef typename DOFVector::ElementType
RF;
222 using DF =
typename EntitySet::Traits::CoordinateField;
224 template<
typename LFSLeaf,
typename TreePath>
225 void leaf(
const LFSLeaf& leaf_lfs, TreePath treePath)
228 auto& fem = leaf_lfs.gridFunctionSpace().finiteElementMap();
232 using Range =
typename LFSLeaf::Traits::GridFunctionSpace::Traits::FiniteElementMap::
233 Traits::FiniteElement::Traits::LocalBasisType::Traits::RangeType;
237 auto coarse_phi = std::vector<Range>{};
238 auto fine_phi = std::vector<Range>{};
240 auto fine_geometry = _current.geometry();
241 auto coarse_geometry = _ancestor.geometry();
244 for (
const auto& ip : QuadratureRules<DF,dim>::rule(_current.type(),_int_order))
246 auto coarse_local = coarse_geometry.local(fine_geometry.global(ip.position()));
247 auto fe = &fem.find(_current);
248 fe->localBasis().evaluateFunction(ip.position(),fine_phi);
249 fe = &fem.find(_ancestor);
250 fe->localBasis().evaluateFunction(coarse_local,coarse_phi);
251 const DF factor = ip.weight()
252 * fine_geometry.integrationElement(ip.position())
253 / coarse_geometry.integrationElement(coarse_local);
255 auto val = Range{0.0};
256 for (size_type i = 0; i < fine_phi.size(); ++i)
258 val.axpy(_u_fine[fine_offset + i],fine_phi[i]);
261 for (size_type i = 0; i < coarse_phi.size(); ++i)
264 for (size_type j = 0; j < inverse_mass_matrix.M(); ++j)
265 x.axpy(inverse_mass_matrix[i][j],coarse_phi[j]);
266 (*_u_coarse)[coarse_offset + i] += factor * (x * val);
279 _u_view.bind(_lfs_cache);
280 _u_coarse = &_transfer_map[_id_set.id(
_element)];
282 _u_view.read(*_u_coarse);
290 while (_ancestor.mightVanish())
293 if (!_ancestor.hasFather())
296 _ancestor = _ancestor.father();
298 _u_coarse = &_transfer_map[_id_set.id(_ancestor)];
300 if (_u_coarse->size() > 0)
303 std::fill(_u_coarse->begin(),_u_coarse->end(),RF(0));
305 for (
const auto& child : descendantElements(_ancestor,max_level))
317 _u_view.bind(_lfs_cache);
318 _u_fine.resize(_lfs_cache.size());
319 _u_view.read(_u_fine);
322 TypeTree::applyToTree(
_lfs,*
this);
329 Projection& projection,
331 LeafOffsetCache& leaf_offset_cache,
332 TransferMap& transfer_map,
333 std::size_t int_order = 2)
336 , _id_set(gfs.gridView().grid().localIdSet())
337 , _projection(projection)
339 , _transfer_map(transfer_map)
342 , _int_order(int_order)
353 typename DOFVector::template ConstLocalView<LFSCache>
_u_view;
365 template<
typename GFS,
typename DOFVector,
typename CountVector>
367 :
public TypeTree::TreeVisitor
368 ,
public TypeTree::DynamicTraversal
376 using IDSet =
typename EntitySet::Traits::GridView::Grid::LocalIdSet;
379 typedef typename DOFVector::ElementType
RF;
384 using DF =
typename EntitySet::Traits::CoordinateField;
386 template<
typename FiniteElement>
389 using Range =
typename FiniteElement::Traits::LocalBasisType::Traits::RangeType;
391 template<
typename X,
typename Y>
394 _phi.resize(_finite_element.localBasis().size());
395 _finite_element.localBasis().evaluateFunction(_coarse_geometry.local(_fine_geometry.global(x)),_phi);
397 for (size_type i = 0; i < _phi.size(); ++i)
398 y.axpy(_dofs[_offset + i],_phi[i]);
402 : _finite_element(finite_element)
403 , _coarse_geometry(coarse_geometry)
404 , _fine_geometry(fine_geometry)
413 mutable std::vector<Range>
_phi;
419 template<
typename LeafLFS,
typename TreePath>
420 void leaf(
const LeafLFS& leaf_lfs, TreePath treePath)
422 using FiniteElement =
typename LeafLFS::Traits::FiniteElementType;
424 auto& fem = leaf_lfs.gridFunctionSpace().finiteElementMap();
431 _u_tmp.resize(fe.localBasis().size());
432 std::fill(_u_tmp.begin(),_u_tmp.end(),RF(0.0));
433 fe.localInterpolation().interpolate(f,_u_tmp);
434 std::copy(_u_tmp.begin(),_u_tmp.end(),_u_fine.begin() + element_offset);
442 _ancestor = ancestor;
443 _u_coarse = &u_coarse;
447 _u_view.bind(_lfs_cache);
450 if (_id_set.id(element) == _id_set.id(ancestor))
453 _u_view.add(*_u_coarse);
457 _u_fine.resize(_lfs_cache.size());
458 std::fill(_u_fine.begin(),_u_fine.end(),RF(0));
460 TypeTree::applyToTree(
_lfs,*
this);
461 _u_view.add(_u_fine);
465 _uc_view.bind(_lfs_cache);
466 _counts.resize(_lfs_cache.size(),1);
467 _uc_view.add(_counts);
471 replay_visitor(
const GFS& gfs, DOFVector& u, CountVector& uc, LeafOffsetCache& leaf_offset_cache)
474 , _id_set(gfs.entitySet().gridView().grid().localIdSet())
486 typename DOFVector::template LocalView<LFSCache>
_u_view;
487 typename CountVector::template LocalView<LFSCache>
_uc_view;
511 template<
class Gr
id,
class GFSU,
class U,
class Projection>
514 typedef typename Grid::LeafGridView LeafGridView;
515 typedef typename LeafGridView::template Codim<0>
516 ::template Partition<Dune::Interior_Partition>::Iterator LeafIterator;
517 typedef typename Grid::template Codim<0>::Entity Element;
518 typedef typename Grid::LocalIdSet IDSet;
519 typedef typename IDSet::IdType ID;
522 typedef std::unordered_map<ID,std::vector<typename U::ElementType> >
MapType;
538 void backupData(Grid& grid, GFSU& gfsu, Projection& projection, U& u, MapType& transfer_map)
545 for(
const auto& cell : elements(gfsu.entitySet(),Partitions::interior))
554 void replayData(Grid& grid, GFSU& gfsu, Projection& projection, U& u,
const MapType& transfer_map)
556 const IDSet& id_set = grid.localIdSet();
559 CountVector uc(gfsu,0);
565 for (
const auto& cell : elements(gfsu.entitySet(),Partitions::interior))
567 Element ancestor = cell;
569 typename MapType::const_iterator map_it;
570 while ((map_it = transfer_map.find(id_set.id(ancestor))) == transfer_map.end())
572 if (!ancestor.hasFather())
574 "transferMap of GridAdaptor didn't contain ancestor of element with id " << id_set.id(ancestor));
575 ancestor = ancestor.father();
578 visitor(cell,ancestor,map_it->second);
582 DOFHandle addHandle1(gfsu,u);
583 gfsu.entitySet().gridView().communicate(addHandle1,
584 Dune::InteriorBorder_InteriorBorder_Interface,Dune::ForwardCommunication);
586 CountHandle addHandle2(gfsu,uc);
587 gfsu.entitySet().gridView().communicate(addHandle2,
588 Dune::InteriorBorder_InteriorBorder_Interface,Dune::ForwardCommunication);
591 typename CountVector::iterator ucit = uc.begin();
592 for (
typename U::iterator uit = u.begin(), uend = u.end(); uit != uend; ++uit, ++ucit)
593 (*uit) /= ((*ucit) > 0 ? (*ucit) : 1.0);
612 template<
class Gr
id,
class GFS,
class X>
616 Projection projection(gfs,int_order);
625 grid_adaptor.
backupData(grid,gfs,projection,x1,transferMap1);
635 grid_adaptor.
replayData(grid,gfs,projection,x1,transferMap1);
652 template<
class Gr
id,
class GFS,
class X>
653 void adapt_grid (Grid& grid, GFS& gfs, X& x1, X& x2,
int int_order)
656 Projection projection(gfs,int_order);
665 grid_adaptor.
backupData(grid,gfs,projection,x1,transferMap1);
667 grid_adaptor.
backupData(grid,gfs,projection,x2,transferMap2);
677 grid_adaptor.
replayData(grid,gfs,projection,x1,transferMap1);
679 grid_adaptor.
replayData(grid,gfs,projection,x2,transferMap2);
689 template <
typename G,
typename... X>
690 struct GFSWithVectors
694 using Tuple = std::tuple<X&...>;
696 GFSWithVectors (GFS& gfs,
int integrationOrder, X&... x) :
698 _integrationOrder(integrationOrder),
703 int _integrationOrder;
708 template <
typename Gr
id>
709 void iteratePacks(Grid& grid);
710 template <
typename Grid,
typename X,
typename... XS>
711 void iteratePacks(Grid& grid, X& x, XS&... xs);
716 template<std::size_t I = 0,
typename Grid,
typename X,
typename... XS>
718 iterateTuple(Grid& grid, X& x, XS&... xs)
721 iteratePacks(grid,xs...);
736 template<std::size_t I = 0,
typename Grid,
typename X,
typename... XS>
738 iterateTuple(Grid& grid, X& x, XS&... xs)
741 using GFS =
typename X::GFS;
742 using Tuple =
typename X::Tuple;
743 using V =
typename std::decay<typename std::tuple_element<I,Tuple>::type>::type;
750 Projection projection(x._gfs,x._integrationOrder);
755 gridAdaptor.
backupData(grid,x._gfs,projection,std::get<I>(x._tuple),transferMap);
759 iterateTuple<I + 1, Grid, X, XS...>(grid,x,xs...);
763 std::get<I>(x._tuple) = V(x._gfs,0.0);
764 gridAdaptor.
replayData(grid,x._gfs,projection,std::get<I>(x._tuple),transferMap);
770 template <
typename Gr
id>
771 void iteratePacks(Grid& grid)
792 template <
typename Grid,
typename X,
typename... XS>
793 void iteratePacks(Grid& grid, X& x, XS&... xs)
795 iterateTuple(grid,x,xs...);
812 template <
typename GFS,
typename... X>
815 impl::GFSWithVectors<GFS,X...> gfsWithVectors(gfs, integrationOrder, x...);
816 return gfsWithVectors;
829 template <
typename Grid,
typename... X>
836 impl::iteratePacks(grid,x...);
844 void error_fraction(
const T& x,
typename T::ElementType alpha,
typename T::ElementType beta,
845 typename T::ElementType& eta_alpha,
typename T::ElementType& eta_beta,
int verbose=0)
848 std::cout <<
"+++ error fraction: alpha=" << alpha <<
" beta=" << beta << std::endl;
850 typedef typename T::ElementType NumberType;
851 NumberType total_error = x.one_norm();
852 NumberType max_error = x.infinity_norm();
853 NumberType eta_alpha_left = 0.0;
854 NumberType eta_alpha_right = max_error;
855 NumberType eta_beta_left = 0.0;
856 NumberType eta_beta_right = max_error;
857 for (
int j=1; j<=steps; j++)
859 eta_alpha = 0.5*(eta_alpha_left+eta_alpha_right);
860 eta_beta = 0.5*(eta_beta_left+eta_beta_right);
861 NumberType sum_alpha=0.0;
862 NumberType sum_beta=0.0;
863 unsigned int alpha_count = 0;
864 unsigned int beta_count = 0;
865 for (
typename T::const_iterator it = x.begin(),
870 if (*it >=eta_alpha) { sum_alpha += *it; alpha_count++;}
871 if (*it < eta_beta) { sum_beta += *it; beta_count++;}
875 std::cout <<
"+++ " << j <<
" eta_alpha=" << eta_alpha <<
" alpha_fraction=" << sum_alpha/total_error
876 <<
" elements: " << alpha_count <<
" of " << x.N() << std::endl;
877 std::cout <<
"+++ " << j <<
" eta_beta=" << eta_beta <<
" beta_fraction=" << sum_beta/total_error
878 <<
" elements: " << beta_count <<
" of " << x.N() << std::endl;
880 if (std::abs(alpha-sum_alpha/total_error) <= 0.01 && std::abs(beta-sum_beta/total_error) <= 0.01)
break;
881 if (sum_alpha>alpha*total_error)
882 eta_alpha_left = eta_alpha;
884 eta_alpha_right = eta_alpha;
885 if (sum_beta>beta*total_error)
886 eta_beta_right = eta_beta;
888 eta_beta_left = eta_beta;
892 std::cout <<
"+++ refine_threshold=" << eta_alpha
893 <<
" coarsen_threshold=" << eta_beta << std::endl;
899 void element_fraction(
const T& x,
typename T::ElementType alpha,
typename T::ElementType beta,
900 typename T::ElementType& eta_alpha,
typename T::ElementType& eta_beta,
int verbose=0)
903 typedef typename T::ElementType NumberType;
904 NumberType total_error =x.N();
905 NumberType max_error = x.infinity_norm();
906 NumberType eta_alpha_left = 0.0;
907 NumberType eta_alpha_right = max_error;
908 NumberType eta_beta_left = 0.0;
909 NumberType eta_beta_right = max_error;
910 for (
int j=1; j<=steps; j++)
912 eta_alpha = 0.5*(eta_alpha_left+eta_alpha_right);
913 eta_beta = 0.5*(eta_beta_left+eta_beta_right);
914 NumberType sum_alpha=0.0;
915 NumberType sum_beta=0.0;
916 unsigned int alpha_count = 0;
917 unsigned int beta_count = 0;
919 for (
typename T::const_iterator it = x.begin(),
924 if (*it>=eta_alpha) { sum_alpha += 1.0; alpha_count++;}
925 if (*it< eta_beta) { sum_beta +=1.0; beta_count++;}
929 std::cout << j <<
" eta_alpha=" << eta_alpha <<
" alpha_fraction=" << sum_alpha/total_error
930 <<
" elements: " << alpha_count <<
" of " << x.N() << std::endl;
931 std::cout << j <<
" eta_beta=" << eta_beta <<
" beta_fraction=" << sum_beta/total_error
932 <<
" elements: " << beta_count <<
" of " << x.N() << std::endl;
934 if (std::abs(alpha-sum_alpha/total_error) <= 0.01 && std::abs(beta-sum_beta/total_error) <= 0.01)
break;
935 if (sum_alpha>alpha*total_error)
936 eta_alpha_left = eta_alpha;
938 eta_alpha_right = eta_alpha;
939 if (sum_beta>beta*total_error)
940 eta_beta_right = eta_beta;
942 eta_beta_left = eta_beta;
946 std::cout <<
"+++ refine_threshold=" << eta_alpha
947 <<
" coarsen_threshold=" << eta_beta << std::endl;
957 typedef typename T::ElementType NumberType;
958 NumberType total_error = x.one_norm();
959 NumberType total_elements = x.N();
960 NumberType max_error = x.infinity_norm();
961 std::vector<NumberType> left(bins,0.0);
962 std::vector<NumberType> right(bins,max_error*(1.0+1
e-8));
963 std::vector<NumberType> eta(bins);
964 std::vector<NumberType> target(bins);
965 for (
unsigned int k=0; k<bins; k++)
966 target[k]= (k+1)/((NumberType)bins);
967 for (
int j=1; j<=steps; j++)
969 for (
unsigned int k=0; k<bins; k++)
970 eta[k]= 0.5*(left[k]+right[k]);
971 std::vector<NumberType> sum(bins,0.0);
972 std::vector<int> count(bins,0);
974 for (
typename T::const_iterator it = x.begin(),
979 for (
unsigned int k=0; k<bins; k++)
991 for (
unsigned int k=0; k<bins; k++)
992 if (sum[k]<=target[k]*total_error)
997 std::vector<NumberType> sum(bins,0.0);
998 std::vector<int> count(bins,0);
999 for (
unsigned int i=0; i<x.N(); i++)
1000 for (
unsigned int k=0; k<bins; k++)
1006 std::cout <<
"+++ error distribution" << std::endl;
1007 std::cout <<
"+++ number of elements: " << x.N() << std::endl;
1008 std::cout <<
"+++ max element error: " << max_error << std::endl;
1009 std::cout <<
"+++ total error: " << total_error << std::endl;
1010 std::cout <<
"+++ bin #elements eta sum/total " << std::endl;
1011 for (
unsigned int k=0; k<bins; k++)
1012 std::cout <<
"+++ " << k+1 <<
" " << count[k] <<
" " << eta[k] <<
" " << sum[k]/total_error << std::endl;
1015 template<
typename Gr
id,
typename X>
1016 void mark_grid (Grid &grid,
const X& x,
typename X::ElementType refine_threshold,
1017 typename X::ElementType coarsen_threshold,
int min_level = 0,
int max_level = std::numeric_limits<int>::max(),
int verbose=0)
1019 typedef typename Grid::LeafGridView GV;
1021 GV gv = grid.leafGridView();
1023 unsigned int refine_cnt=0;
1024 unsigned int coarsen_cnt=0;
1026 typedef typename X::GridFunctionSpace GFS;
1029 typedef typename X::template ConstLocalView<LFSCache> XView;
1031 LFS lfs(x.gridFunctionSpace());
1032 LFSCache lfs_cache(lfs);
1035 for(
const auto& cell : elements(gv))
1039 x_view.bind(lfs_cache);
1041 if (x_view[0]>=refine_threshold && cell.level() < max_level)
1046 if (x_view[0]<=coarsen_threshold && cell.level() > min_level)
1054 std::cout <<
"+++ mark_grid: " << refine_cnt <<
" marked for refinement, " 1055 << coarsen_cnt <<
" marked for coarsening" << std::endl;
1059 template<
typename Gr
id,
typename X>
1061 typename X::ElementType coarsen_threshold,
int verbose=0)
1063 typedef typename Grid::LeafGridView GV;
1065 GV gv = grid.leafGridView();
1067 unsigned int coarsen_cnt=0;
1069 typedef typename X::GridFunctionSpace GFS;
1072 typedef typename X::template ConstLocalView<LFSCache> XView;
1074 LFS lfs(x.gridFunctionSpace());
1075 LFSCache lfs_cache(lfs);
1078 for(
const auto& cell : elements(gv))
1082 x_view.bind(lfs_cache);
1084 if (x_view[0]>=refine_threshold)
1089 if (x_view[0]<=coarsen_threshold)
1096 std::cout <<
"+++ mark_grid_for_coarsening: " 1097 << coarsen_cnt <<
" marked for coarsening" << std::endl;
1105 double optimistic_factor;
1106 double coarsen_limit;
1107 double balance_limit;
1112 double refine_fraction_while_refinement;
1113 double coarsen_fraction_while_refinement;
1114 double coarsen_fraction_while_coarsening;
1115 double timestep_decrease_factor;
1116 double timestep_increase_factor;
1126 bool have_decreased_time_step;
1127 bool have_refined_grid;
1130 double accumulated_estimated_error_squared;
1131 double minenergy_rate;
1135 : scaling(16.0), optimistic_factor(1.0), coarsen_limit(0.5), balance_limit(0.33333),
1136 tol(tol_), T(T_), verbose(verbose_), no_adapt(false),
1137 refine_fraction_while_refinement(0.7),
1138 coarsen_fraction_while_refinement(0.2),
1139 coarsen_fraction_while_coarsening(0.2),
1140 timestep_decrease_factor(0.5), timestep_increase_factor(1.5),
1141 accept(false), adapt_dt(false), adapt_grid(false), newdt(1.0),
1142 have_decreased_time_step(false), have_refined_grid(false),
1143 accumulated_estimated_error_squared(0.0),
1150 timestep_decrease_factor=
s;
1155 timestep_increase_factor=
s;
1160 refine_fraction_while_refinement=
s;
1165 coarsen_fraction_while_refinement=
s;
1170 coarsen_fraction_while_coarsening=
s;
1195 optimistic_factor=
s;
1245 return accumulated_estimated_error_squared;
1251 have_decreased_time_step =
false;
1252 have_refined_grid =
false;
1255 template<
typename GM,
typename X>
1256 void evaluate_estimators (GM& grid,
double time,
double dt,
const X& eta_space,
const X& eta_time,
double energy_timeslab)
1263 double spatial_error = eta_space.one_norm();
1264 double temporal_error = scaling*eta_time.one_norm();
1265 double sum = spatial_error + temporal_error;
1267 double allowed = tol*tol*(energy_timeslab+minenergy_rate*dt);
1268 q_s = spatial_error/sum;
1269 q_t = temporal_error/sum;
1276 <<
" sum/allowed=" << sum/allowed
1278 <<
" estimated error=" << sqrt(accumulated_estimated_error_squared+sum)
1279 <<
" energy_rate=" << energy_timeslab/dt
1286 accumulated_estimated_error_squared += sum;
1287 if (verbose>1) std::cout <<
"+++ no adapt mode" << std::endl;
1296 if (verbose>1) std::cout <<
"+++ accepting time step" << std::endl;
1297 accumulated_estimated_error_squared += sum;
1300 if (sum<coarsen_limit*allowed)
1303 if (q_t<balance_limit)
1306 newdt = timestep_increase_factor*dt;
1308 if (verbose>1) std::cout <<
"+++ spatial error dominates: increase time step" << std::endl;
1312 if (q_s>balance_limit)
1315 newdt = timestep_increase_factor*dt;
1317 if (verbose>1) std::cout <<
"+++ increasing time step" << std::endl;
1320 double eta_refine, eta_coarsen;
1321 if (verbose>1) std::cout <<
"+++ mark grid for coarsening" << std::endl;
1324 coarsen_fraction_while_coarsening,eta_refine,eta_coarsen);
1332 if (q_t<balance_limit)
1335 newdt = timestep_increase_factor*dt;
1337 if (verbose>1) std::cout <<
"+++ spatial error dominates: increase time step" << std::endl;
1344 if (verbose>1) std::cout <<
"+++ will redo time step" << std::endl;
1345 if (q_t>1-balance_limit)
1348 newdt = timestep_decrease_factor*dt;
1350 have_decreased_time_step =
true;
1351 if (verbose>1) std::cout <<
"+++ decreasing time step only" << std::endl;
1355 if (q_t<balance_limit)
1357 if (!have_decreased_time_step)
1360 newdt = timestep_increase_factor*dt;
1362 if (verbose>1) std::cout <<
"+++ increasing time step" << std::endl;
1368 newdt = timestep_decrease_factor*dt;
1370 have_decreased_time_step =
true;
1371 if (verbose>1) std::cout <<
"+++ decreasing time step" << std::endl;
1374 double eta_refine, eta_coarsen;
1375 if (verbose>1) std::cout <<
"+++ BINGO mark grid for refinement and coarsening" << std::endl;
1378 coarsen_fraction_while_refinement,eta_refine,eta_coarsen,0);
const Cell & _element
Definition: adaptivity.hh:126
std::vector< RF > LocalDOFVector
Definition: adaptivity.hh:380
std::vector< typename CountVector::ElementType > LocalCountVector
Definition: adaptivity.hh:381
typename EntitySet::Traits::GridView::Grid::LocalIdSet IDSet
Definition: adaptivity.hh:376
void error_distribution(const T &x, unsigned int bins)
Definition: adaptivity.hh:954
std::size_t size_type
Definition: adaptivity.hh:221
bool adaptGrid() const
Definition: adaptivity.hh:1218
LocalDOFVector * _u_coarse
Definition: adaptivity.hh:355
Iterator extract_lfs_leaf_sizes(const LFS &lfs, Iterator it)
Definition: lfsindexcache.hh:165
replay_visitor(const GFS &gfs, DOFVector &u, CountVector &uc, LeafOffsetCache &leaf_offset_cache)
Definition: adaptivity.hh:471
typename EntitySet::Traits::CoordinateField DF
Definition: adaptivity.hh:222
void setCoarsenFractionWhileRefinement(double s)
Definition: adaptivity.hh:1163
static const unsigned int value
Definition: gridfunctionspace/tags.hh:177
typename GFS::Traits::EntitySet EntitySet
Definition: adaptivity.hh:375
LocalCountVector _counts
Definition: adaptivity.hh:493
Class for automatic adaptation of the grid.
Definition: adaptivity.hh:512
LFSIndexCache< LFS > LFSCache
Definition: adaptivity.hh:372
Element::Geometry Geometry
Definition: adaptivity.hh:211
Base class for all PDELab exceptions.
Definition: exceptions.hh:17
const LocalDOFVector & _dofs
Definition: adaptivity.hh:412
Element _ancestor
Definition: adaptivity.hh:485
Projection & _projection
Definition: adaptivity.hh:352
Element _ancestor
Definition: adaptivity.hh:350
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:844
size_type _int_order
Definition: adaptivity.hh:357
DynamicMatrix< typename U::ElementType > MassMatrix
Definition: adaptivity.hh:153
static const int dim
Definition: adaptivity.hh:83
Definition: genericdatahandle.hh:622
void evaluate(const X &x, Y &y) const
Definition: adaptivity.hh:392
DOFVector::template ConstLocalView< LFSCache > _u_view
Definition: adaptivity.hh:353
double qs() const
Definition: adaptivity.hh:1228
LFS _lfs
Definition: adaptivity.hh:69
void startTimeStep()
Definition: adaptivity.hh:1249
const FiniteElement & _finite_element
Definition: adaptivity.hh:409
void adapt_grid(Grid &grid, GFS &gfs, X &x1, int int_order)
adapt a grid, corresponding function space and solution vectors
Definition: adaptivity.hh:613
Definition: adaptivity.hh:27
L2Projection< typename LFS::Traits::GridFunctionSpace, DOFVector > Projection
Definition: adaptivity.hh:217
TransferMap::mapped_type LocalDOFVector
Definition: adaptivity.hh:214
Projection::MassMatrices MassMatrices
Definition: adaptivity.hh:218
GridAdaptor(const GFSU &gfs)
The constructor.
Definition: adaptivity.hh:529
LocalFunctionSpace< GFS > LFS
Definition: adaptivity.hh:36
void operator()(const Element &element, const Element &ancestor, const LocalDOFVector &u_coarse)
Definition: adaptivity.hh:439
LeafOffsetCache & _leaf_offset_cache
Definition: adaptivity.hh:356
typename Element::Geometry Geometry
Definition: adaptivity.hh:378
std::unordered_map< ID, std::vector< typename U::ElementType > > MapType
Definition: adaptivity.hh:522
CountVector::template LocalView< LFSCache > _uc_view
Definition: adaptivity.hh:487
void replayData(Grid &grid, GFSU &gfsu, Projection &projection, U &u, const MapType &transfer_map)
Definition: adaptivity.hh:554
size_type _leaf_index
Definition: adaptivity.hh:490
Definition: adaptivity.hh:366
Element _element
Definition: adaptivity.hh:349
double qt() const
Definition: adaptivity.hh:1233
Definition: adaptivity.hh:144
const std::size_t offset
Definition: localfunctionspace.hh:74
GFS::Traits::GridView::template Codim< 0 >::Entity Cell
Definition: adaptivity.hh:35
const Entity & e
Definition: localfunctionspace.hh:111
void leaf(const LeafLFS &leaf_lfs, TreePath treePath)
Definition: adaptivity.hh:420
LeafOffsetCache & _leaf_offset_cache
Definition: adaptivity.hh:489
Definition: adaptivity.hh:1101
MassMatrices & _mass_matrices
Definition: adaptivity.hh:127
Element _current
Definition: adaptivity.hh:351
TimeAdaptationStrategy(double tol_, double T_, int verbose_=0)
Definition: adaptivity.hh:1134
LocalFunctionSpace< GFS > LFS
Definition: adaptivity.hh:371
Traits::IndexContainer::size_type size() const
get current size
Definition: localfunctionspace.hh:206
typename EntitySet::Traits::GridView::Grid::LocalIdSet IDSet
Definition: adaptivity.hh:209
LocalDOFVector _u_tmp
Definition: adaptivity.hh:492
array< std::size_t, TypeTree::TreeInfo< GFS >::leafCount+1 > LeafOffsets
Definition: adaptivity.hh:40
void setRefineFractionWhileRefinement(double s)
Definition: adaptivity.hh:1158
Geometry _coarse_geometry
Definition: adaptivity.hh:410
Element _element
Definition: adaptivity.hh:484
Definition: adaptivity.hh:387
void evaluate_estimators(GM &grid, double time, double dt, const X &eta_space, const X &eta_time, double energy_timeslab)
Definition: adaptivity.hh:1256
void setTimeStepIncreaseFactor(double s)
Definition: adaptivity.hh:1153
const MassMatrices & inverseMassMatrices(const Element &e)
Calculate the inverse local mass matrix, used in the local L2 projection.
Definition: adaptivity.hh:171
Dune::PDELab::LeafOffsetCache< GFS > LeafOffsetCache
Definition: adaptivity.hh:206
void backupData(Grid &grid, GFSU &gfsu, Projection &projection, U &u, MapType &transfer_map)
Definition: adaptivity.hh:538
void adaptGrid(Grid &grid, X &...x)
Adapt grid and multiple function spaces with corresponding vectors.
Definition: adaptivity.hh:830
const LeafOffsets & operator[](GeometryType gt) const
Definition: adaptivity.hh:42
std::size_t size_type
Definition: adaptivity.hh:383
Projection::MassMatrix MassMatrix
Definition: adaptivity.hh:219
void setTimeStepDecreaseFactor(double s)
Definition: adaptivity.hh:1148
void setTemporalScaling(double s)
Definition: adaptivity.hh:1188
LFSCache _lfs_cache
Definition: adaptivity.hh:482
void update(const Cell &e)
Definition: adaptivity.hh:50
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:1016
const IDSet & _id_set
Definition: adaptivity.hh:348
void setCoarsenLimit(double s)
Definition: adaptivity.hh:1178
bool adaptDT() const
Definition: adaptivity.hh:1213
DOFVector::template LocalView< LFSCache > _u_view
Definition: adaptivity.hh:486
Geometry _fine_geometry
Definition: adaptivity.hh:411
LFS _lfs
Definition: adaptivity.hh:346
LeafOffsetCache(const GFS &gfs)
Definition: adaptivity.hh:64
typename EntitySet::Traits::CoordinateField DF
Definition: adaptivity.hh:384
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:328
const IDSet & _id_set
Definition: adaptivity.hh:483
typename EntitySet::Element Element
Definition: adaptivity.hh:377
bool acceptTimeStep() const
Definition: adaptivity.hh:1208
double accumulatedErrorSquared() const
Definition: adaptivity.hh:1243
Definition: adaptivity.hh:32
size_type _leaf_index
Definition: adaptivity.hh:358
void setMinEnergyRate(double s)
Definition: adaptivity.hh:1173
LocalFunctionSpace< GFS > LFS
Definition: adaptivity.hh:204
LocalDOFVector _u_fine
Definition: adaptivity.hh:491
DOFVector::ElementType RF
Definition: adaptivity.hh:213
double newDT() const
Definition: adaptivity.hh:1223
void setOptimisticFactor(double s)
Definition: adaptivity.hh:1193
LFSCache _lfs_cache
Definition: adaptivity.hh:347
DOFVector::ElementType RF
Definition: adaptivity.hh:379
typename impl::BackendVectorSelector< GridFunctionSpace, FieldType >::Type Vector
alias of the return type of BackendVectorSelector
Definition: backend/interface.hh:113
LFSIndexCache< LFS > LFSCache
Definition: adaptivity.hh:205
TransferMap & _transfer_map
Definition: adaptivity.hh:354
size_type _leaf_index
Definition: adaptivity.hh:129
typename GFS::Traits::EntitySet EntitySet
Definition: adaptivity.hh:208
const std::string s
Definition: function.hh:1102
void setAdaptationOn()
Definition: adaptivity.hh:1198
Definition: adaptivity.hh:199
std::vector< LeafOffsets > _leaf_offset_cache
Definition: adaptivity.hh:70
LFS _lfs
Definition: adaptivity.hh:481
void leaf(const LFSLeaf &leaf_lfs, TreePath treePath)
Definition: adaptivity.hh:225
typename FiniteElement::Traits::LocalBasisType::Traits::RangeType Range
Definition: adaptivity.hh:389
coarse_function(const FiniteElement &finite_element, Geometry coarse_geometry, Geometry fine_geometry, const LocalDOFVector &dofs, size_type offset)
Definition: adaptivity.hh:401
typename EntitySet::Element Element
Definition: adaptivity.hh:210
const GFS & gridFunctionSpace() const
Returns the GridFunctionSpace underlying this LocalFunctionSpace.
Definition: localfunctionspace.hh:264
std::vector< Range > _phi
Definition: adaptivity.hh:413
size_type _offset
Definition: adaptivity.hh:414
impl::GFSWithVectors< GFS, X... > transferSolutions(GFS &gfs, int integrationOrder, X &...x)
Pack function space and vectors for grid adaption.
Definition: adaptivity.hh:813
std::array< MassMatrix, TypeTree::TreeInfo< GFS >::leafCount > MassMatrices
Definition: adaptivity.hh:154
void setAdaptationOff()
Definition: adaptivity.hh:1203
const QuadratureRule< DF, dim > & _quadrature_rule
Definition: adaptivity.hh:128
LocalDOFVector _u_fine
Definition: adaptivity.hh:359
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:1060
L2Projection(const GFS &gfs, int intorder=2)
The constructor.
Definition: adaptivity.hh:160
double endT() const
Definition: adaptivity.hh:1238
const LocalDOFVector * _u_coarse
Definition: adaptivity.hh:488
void setBalanceLimit(double s)
Definition: adaptivity.hh:1183
Dune::PDELab::LeafOffsetCache< GFS > LeafOffsetCache
Definition: adaptivity.hh:373
void setCoarsenFractionWhileCoarsening(double s)
Definition: adaptivity.hh:1168
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:899
void operator()(const Element &element)
Definition: adaptivity.hh:273