diff --git a/dune/solvers/common/permutationmanager.hh b/dune/solvers/common/permutationmanager.hh index b9db3b865188ba3d725b2a2e7b12f434d28578bd..1adc7597bf37b965dd5a87c93249dfcb1bde7d0d 100644 --- a/dune/solvers/common/permutationmanager.hh +++ b/dune/solvers/common/permutationmanager.hh @@ -24,9 +24,6 @@ // int get_number_of_blocks() // number of diagonal blocks // int get_block_size( int block_number) // get the size of the block_number'th diagonal block - -// PermutationManager only for dim = 2 or dim = 3. - //! Idea of implementation: define an ordering for the nodes (with respect to the axis of anisotropy) //! and then use a merge sort algorithm to renumerate the nodes along the lines which are parallel to this axis. //! to do so, we introduce a 'NodeList' class. @@ -104,13 +101,16 @@ template <class CoordinateImp, int dim_domain> class AxisComparison { - int axisOfAnisotropy_; + int primaryAxis_; public: - AxisComparison(int axisOfAnisotropy) - : axisOfAnisotropy_(axisOfAnisotropy) - {} + AxisComparison(int primaryAxis) + : primaryAxis_(primaryAxis) + { + if ( primaryAxis_ >= dim_domain) + DUNE_THROW(Dune::GridError, "Number of axis of anisotropy > space dim."); + } // lexicographic order '<' on the set of nodes (depending on the dimension and the axis of anisotropy) // this order guarantees that the nodes on the anisotropic lines follow each other @@ -119,56 +119,44 @@ public: CoordinateImp coordinate1 = node1.get_global_coordinate(); CoordinateImp coordinate2 = node2.get_global_coordinate(); - - int primary_axis = axisOfAnisotropy_; - int secondary_axis = -1; - int third_axis = -1; - - if (dim_domain == 2) - { - if ( primary_axis == 1 ) - {secondary_axis = 0;} - else - {secondary_axis = 1;} - - if ( primary_axis > 1) - {std :: cout << "Error. Number of axis of anisotropy > space dim." << std :: endl; - abort();} - - if ( (coordinate1[secondary_axis] < coordinate2[secondary_axis]) || - ((coordinate1[secondary_axis] == coordinate2[secondary_axis])&&(coordinate1[primary_axis] < coordinate2[primary_axis])) ) - {return true;} - else - {return false;} - } - - if (dim_domain == 3) - { - - if ( primary_axis == 2 ) { - secondary_axis = 1; - third_axis = 0;} - - if ( primary_axis == 1 ) { - secondary_axis = 0; - third_axis = 2;} - - if ( primary_axis == 0 ) { - secondary_axis = 1; - third_axis = 2;} - - if ( (coordinate1[third_axis] < coordinate2[third_axis]) || - ((coordinate1[third_axis] == coordinate2[third_axis])&&(coordinate1[secondary_axis] < coordinate2[secondary_axis])) || - ((coordinate1[third_axis] == coordinate2[third_axis])&&(coordinate1[secondary_axis] == coordinate2[secondary_axis])&&(coordinate1[primary_axis] < coordinate2[primary_axis])) ) - {return true;} - else - {return false;} - } - - std :: cout << "Class only implemented for dim=2 or dim=3" << std :: cout; + + // axis "rank": primaryAxis = axis_of_anisotropy + // axis_of_anisotropy, 0, 1, ..., (axis_of_anisotropy-1), (axis_of_anisotropy+1), ..., dim-1 + + int current_axis = dim_domain-1; + + while ( current_axis >= 0 ) + { + if ( current_axis == primaryAxis_ ) + current_axis -= 1; + + if ( current_axis < 0 ) + break; + + if ( coordinate1[ current_axis ] < coordinate2[ current_axis ] ) + return true; + else if ( coordinate1[ current_axis ] > coordinate2[ current_axis ] ) + return false; + + current_axis -= 1; + + } + + // the case current_axis = primaryAxis_: + + if ( coordinate1[ primaryAxis_ ] < coordinate2[ primaryAxis_ ] ) + return true; + +#if 0 + else if ( coordinate1[ primaryAxis_ ] > coordinate2[ primaryAxis_ ] ) + { return false; } + + std :: cout << "Error in PermutationManager. Coordinate 1 equals Coordinate 2. Less-equal relation can not be found." << std :: endl; abort(); +#endif + return false; - + } }; @@ -198,17 +186,15 @@ public: typedef Dune::FieldVector<double, dim> DomainPointType; // type of point in the domain typedef Node< DomainPointType , dim > NodeType; // type of nodes - typedef std::vector< NodeType > NodeListType; // array type list of nodes + typedef std::vector< NodeType > NodeListType; // array type list of nodes - typedef Dune::FieldVector<typename GridType::ctype, GridType::dimensionworld> CoordinateType; + typedef Dune::FieldVector<typename GridType::ctype, GridType::dimensionworld> CoordinateType; - - GridType *grid_; - const int grid_level_; + const GridViewType gridView_; // Let 'p' denote the permutation of the node numeration - NodeListType &node_list_; - NodeListType &original_node_list_; //original node list numeration + NodeListType node_list_; + NodeListType original_node_list_; //original node list numeration // it holds: node_list_[p(i)] = original_node_list_[i] int axis_of_anisotropy_; @@ -228,11 +214,10 @@ public: } // Constructor for default axis of anisotropy = dim-1 - inline explicit PermutationManager ( GridType &grid, const int grid_level ) - : grid_( &grid ), - grid_level_( grid_level ), - node_list_( *(new NodeListType( grid.size(grid_level,dim) )) ), // grid_level = level, dim = codim --> the set of nodes at a certain level - original_node_list_( *(new NodeListType( grid.size(grid_level,dim) )) ), + inline explicit PermutationManager ( GridViewType gridView) + : gridView_( gridView ), + node_list_( gridView.size(dim) ), + original_node_list_( gridView.size(dim) ), axis_of_anisotropy_( dim-1 ), number_of_blocks_( 0 ), blocks_set_up_( false ), @@ -242,11 +227,10 @@ public: } // Constructor for setting number of the axis along which we have the anisotropy (0=x, 1=y, 2=z) - inline explicit PermutationManager ( GridType &grid, const int grid_level, const int axis_of_anisotropy) - : grid_( &grid ), - grid_level_( grid_level ), - node_list_( *(new NodeListType( grid.size(grid_level,dim))) ), // grid_level = level, dim = codim --> the set of nodes at a certain level - original_node_list_( *(new NodeListType( grid.size(grid_level,dim))) ), + inline explicit PermutationManager ( GridViewType gridView, const int axis_of_anisotropy) + : gridView_( gridView ), + node_list_( gridView.size(dim) ), + original_node_list_( gridView.size(dim) ), axis_of_anisotropy_( axis_of_anisotropy ), number_of_blocks_( 0 ), blocks_set_up_( false ), @@ -261,12 +245,10 @@ private: void set_up_node_lists() { - GridViewType level_grid_view = grid_->levelView( grid_level_ ); - - NodeIterator starting_node = level_grid_view.template begin< dim >(); - NodeIterator end_node = level_grid_view.template end< dim >(); + NodeIterator starting_node = gridView_.template begin< dim >(); + NodeIterator end_node = gridView_.template end< dim >(); - const IndexSetType &idSet = level_grid_view.indexSet(); + const IndexSetType &idSet = gridView_.indexSet(); // create the node list for our grid // ----------- Iterator over the nodes of the grid -------------------- @@ -423,7 +405,7 @@ public: if (permutation_set_up_ == false) { std :: cout << "Note that the permutation has not yet been set up. Therefore the following output is still default and the permutation equals the identity." << std :: endl;} - for( int i = 0; i < node_list_.get_size(); ++i ) + for( int i = 0; i < node_list_.size(); ++i ) { std :: cout << "----------------------------------" << std :: endl; std :: cout << "i = " << i << std :: endl; @@ -432,9 +414,9 @@ public: std :: cout << "----------------------------------" << std :: endl; if (!(i == get_permuted_index(get_inverse_permuted_index(i)) )) - { std :: cout << "Error! Permutation of inverse permutation is not equal to identiy. Check Permutation Manager for grid level " << grid_level_ << "." << std :: endl; } + { std :: cout << "Error! Permutation of inverse permutation is not equal to identity." << std :: endl; } if (!(i == get_inverse_permuted_index(get_permuted_index(i)) )) - { std :: cout << "Error! Inverse permutation of permutation is not equal to identiy. Check Permutation Manager for grid level " << grid_level_ << "." << std :: endl; } + { std :: cout << "Error! Inverse permutation of permutation is not equal to identity." << std :: endl; } } } @@ -443,12 +425,12 @@ public: if (permutation_set_up_ == false) { std :: cout << "Note that the permutation has not yet been set up. Therefore the following output is still default and the permutation equals the identity." << std :: endl;} - for( int i = 0; i < node_list_.get_size(); ++i ) + for( int i = 0; i < node_list_.size(); ++i ) { if (!(i == get_permuted_index(get_inverse_permuted_index(i)) )) - { std :: cout << "Error! Permutation of inverse permutation is not equal to identiy. Check Permutation Manager for grid level " << grid_level_ << "." << std :: endl; } + { std :: cout << "Error! Permutation of inverse permutation is not equal to identity." << std :: endl; } if (!(i == get_inverse_permuted_index(get_permuted_index(i)) )) - { std :: cout << "Error! Inverse permutation of permutation is not equal to identiy. Check Permutation Manager for grid level " << grid_level_ << "." << std :: endl; } + { std :: cout << "Error! Inverse permutation of permutation is not equal to identity." << std :: endl; } } } @@ -459,7 +441,7 @@ public: if (permutation_set_up_ == false) { std :: cout << "Note that the permutation has not yet been set up. Therefore the following output is still default and the permutation equals the identity." << std :: endl;} - for( int i = 0; i < node_list_.get_size(); ++i ) + for( int i = 0; i < node_list_.size(); ++i ) { std::cout << "Original node[" << get_inverse_permuted_index( i ) << "] = (" << node_list_[i].get_global_coordinate() << ") recieved permuted index " << i << std::endl; // alternative: @@ -467,7 +449,7 @@ public: } #if 0 - for( int i = 0; i < node_list_.get_size(); ++i ) + for( int i = 0; i < node_list_.size(); ++i ) { std::cout << "Permuted node[" << i << "] = (" << original_node_list_[get_inverse_permuted_index( i )].get_global_coordinate() << ") has original index " << get_inverse_permuted_index( i ) << std::endl; }