diff --git a/utils/bisection-compatibility/estimate-closure.cc b/utils/bisection-compatibility/estimate-closure.cc index b53ffa5e156afdb05418eff3265acf7528065505..1f1b4ad4c39156a2c18af4267433a0408bd19b20 100644 --- a/utils/bisection-compatibility/estimate-closure.cc +++ b/utils/bisection-compatibility/estimate-closure.cc @@ -15,18 +15,37 @@ template< class GridType > -void estimateClosure ( GridType &grid ) +void estimateClosure ( GridType &grid , int level = 0 ) { - const std::size_t macroSize =grid.levelGridView(0).size(0); + //refine up to level + for(int i = 0; i < level; ++i) + { + for( const auto & entity : Dune::elements( grid.leafGridView( ) ) ) + { + if( entity.level() < level ) + { + grid.mark( 1, entity ); + } + else if( entity.level() > level ) + { + grid.mark( -1, entity ); + } + } + grid.preAdapt(); + grid.adapt(); + grid.postAdapt(); + } + const size_t macroSize = grid.leafGridView().size( 0 ); + const size_t nElements = grid.levelGridView( level ).size( 0 ); - typedef typename GridType::LevelIndexSet LevelIndexSetType; - typedef typename LevelIndexSetType::IndexType IdType; + typedef typename GridType::GlobalIdSet IndexSetType; + typedef typename IndexSetType::IdType IdType; - const LevelIndexSetType & macroIdSet = grid.levelIndexSet(0); + const IndexSetType & macroIdSet = grid.globalIdSet(); std::map< IdType, size_t > elementClosure; typedef typename GridType::LeafGridView LeafGridViewType; - Dune::VTKSequenceWriter< LeafGridViewType > vtkout( grid.leafGridView(), "solution", "./", ".", Dune::VTK::nonconforming ); + Dune::VTKSequenceWriter< LeafGridViewType > vtkout( grid.leafGridView(), "solution" + level, "./", ".", Dune::VTK::nonconforming ); std::shared_ptr< VolumeData< LeafGridViewType > > ptr( new VolumeData< LeafGridViewType > ( ) ); vtkout.addCellData( ptr ); @@ -35,10 +54,14 @@ void estimateClosure ( GridType &grid ) size_t maxClosure = 0; //loop over all macro elements. - for( const auto & macroEntity : Dune::elements( grid.levelGridView(0) ) ) + for( const auto & macroEntity : Dune::elements( grid.levelGridView( level ) ) ) { //if elementClosure calculated - continue - const IdType macroId = macroIdSet.index(macroEntity); + const IdType macroId = macroIdSet.id( macroEntity ); + if( elementClosure.find( macroId ) != elementClosure.end() ) + { + continue; + } //mark for refinement grid.mark( 1, macroEntity ); //adapt @@ -47,25 +70,49 @@ void estimateClosure ( GridType &grid ) grid.postAdapt(); //loop over macro elements size_t closure = grid.leafGridView().size(0) - macroSize; + for( const auto & entity : Dune::elements( grid.levelGridView( level ) ) ) + { + const auto id = macroIdSet.id( entity ); + if( ! entity.isLeaf() ) + { + elementClosure[ id ] = elementClosure[ id ] == 0 ? closure : std::min( closure, elementClosure[ id ] ); + } + } elementClosure[ macroId ] = closure; if(closure > maxClosure) { maxClosure = closure; + std::cout << "New maximum: " << maxClosure << std::endl; vtkout.write( double( closure ) ); } - while( grid.leafGridView().size(0) != macroSize ) + std::cout << "Before: "<< grid.leafGridView().size(0) << ", after: "; + //refine up to level + for(int i = 0; i < closure * (level + 1); ++i) + { + if( grid.leafGridView().size(0) == macroSize ) break; + for( const auto & entity : Dune::elements( grid.leafGridView( ) ) ) + { + grid.mark( -1, entity ); + } + grid.preAdapt(); + grid.adapt(); + grid.postAdapt(); + } + //refine up to level + for(int i = 0; i < level ; ++i) { - for( const auto & entity : Dune::elements( grid.leafGridView() ) ) + for( const auto & entity : Dune::elements( grid.leafGridView( ) ) ) { - if(entity.level() > 0 ) + if( entity.level() < level ) { - grid.mark( -1, entity ); + grid.mark( 1, entity ); } } grid.preAdapt(); grid.adapt(); grid.postAdapt(); } + std::cout << grid.leafGridView().size(0) << std::endl; } size_t minClosure = 1e10; @@ -77,7 +124,7 @@ void estimateClosure ( GridType &grid ) minClosure = std::min( closure.second , minClosure ); avgClosure += closure.second ; } - avgClosure /= macroSize; + avgClosure /= nElements; std::cout << "Closure (min, max, avg): " << minClosure << " " << maxClosure << " " << avgClosure << std::endl << std::endl; } @@ -96,6 +143,7 @@ int main (int argc , char **argv) { GridType & grid = *gridPtr; grid.loadBalance(); estimateClosure( grid ); + estimateClosure( grid, 3); } } catch( Dune::Exception &e )