Select Git revision
permutationmanager.hh
Code owners
Assign users and groups as approvers for specific file changes. Learn more.
permutationmanager.hh 14.07 KiB
// -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=8 sw=4 sts=4:
#ifndef DUNE_PERMUTATIONMANAGER_HH
#define DUNE_PERMUTATIONMANAGER_HH
#include <algorithm>
#include <vector>
#include <algorithm>
#include <dune/common/fvector.hh>
#include <dune/grid/common/exceptions.hh>
//! This file contains the permutation manager class.
//! 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.
//! everything is handeled by the PermutationManager
// type of nodes:
template< class CoordinateImp, int dim_domain >
class Node{
public:
static const int dim = dim_domain;
typedef CoordinateImp CoordinateType;
CoordinateType global_coordinate_;
int node_index_; //original number of node (no permutation)
int permuted_node_index_; //permuted number of node
public:
// Constructor
inline explicit Node ( const int &node_index, const CoordinateType &global_coordinate )
: global_coordinate_( global_coordinate ),
node_index_( node_index ),
permuted_node_index_( -1 )
{
}
// Constructor
inline explicit Node ( )
: node_index_( -1 ),
permuted_node_index_( -1 )
{
}
inline void set_global_coordinate ( const CoordinateType &global_coordinate )
{
global_coordinate_ = global_coordinate;
}
inline CoordinateType get_global_coordinate () const
{
return global_coordinate_;
}
inline void set_node_index ( const int &node_index )
{
node_index_ = node_index;
}
inline int get_node_index () const
{ return node_index_;
}
inline void set_permuted_node_index ( const int &permuted_node_index )
{
permuted_node_index_ = permuted_node_index;
}
inline int get_permuted_node_index () const
{
return permuted_node_index_;
}
};
template <class CoordinateImp, int dim_domain>
class AxisComparison
{
int primaryAxis_;
public:
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
bool operator() ( const Node<CoordinateImp,dim_domain> &node1, const Node<CoordinateImp,dim_domain> &node2)
{
CoordinateImp coordinate1 = node1.get_global_coordinate();
CoordinateImp coordinate2 = node2.get_global_coordinate();
// 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;
}
};
// the information about the axis of anisotropy is in NodeList, since this list contains a methode "less_equal", which describes the ordering of the nodes.
// this ordering depends on the axis of anisotropy
template< class GridViewImp >
class PermutationManager{
public:
#if 1
typedef GridViewImp GridViewType;
typedef typename GridViewType::Grid GridType;
typedef typename GridViewType::IndexSet IndexSetType;
#endif
static const int dim = GridType::dimension;
static const int dimworld = GridType::dimensionworld;
typedef typename GridViewType :: template Codim<dim> :: Iterator NodeIterator;
typedef typename GridType :: template Codim<dim> :: Geometry NodeGeometry;
typedef Dune::FieldVector<double, dimworld> DomainPointType; // type of point in the domain
typedef Dune::FieldVector<typename GridType::ctype, dimworld> CoordinateType;
typedef Node< CoordinateType, dim > NodeType; // type of nodes
typedef std::vector< NodeType > NodeListType; // array type list of nodes
const GridViewType gridView_;
// Let 'p' denote the permutation of the node 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_;
int number_of_blocks_;
int *block_sizes;
bool blocks_set_up_;
bool permutation_set_up_; // is permutation set up?
public:
// Constructor for setting number of the axis along which we have the anisotropy (0=x, 1=y, 2=z)
inline explicit PermutationManager ( GridViewType gridView, const int axis_of_anisotropy=dim-1)
: 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 ),
permutation_set_up_( false )
{
set_up_node_lists();
}
private:
//! ----------------------------------- the "set-up-node-list"-part -----------------------------------------------------------------
void set_up_node_lists()
{
NodeIterator starting_node = gridView_.template begin< dim >();
NodeIterator end_node = gridView_.template end< dim >();
const IndexSetType &idSet = gridView_.indexSet();
// create the node list for our grid
// ----------- Iterator over the nodes of the grid --------------------
int number_of_node = 0;
for (NodeIterator node = starting_node; node!=end_node; ++node)
{
const NodeGeometry& node_geometry = node->geometry();
CoordinateType global_coordinate_node = node_geometry.corner(0);
node_list_[idSet.index(*node)].set_node_index( idSet.index(*node) );
node_list_[idSet.index(*node)].set_permuted_node_index( idSet.index(*node) );
node_list_[idSet.index(*node)].set_global_coordinate( global_coordinate_node );
original_node_list_[idSet.index(*node)].set_node_index( idSet.index(*node) );
original_node_list_[idSet.index(*node)].set_permuted_node_index( idSet.index(*node) );
original_node_list_[idSet.index(*node)].set_global_coordinate( global_coordinate_node );
number_of_node += 1;
}
// ------------ end: Iterator over the nodes of the grid --------------
}
//! ----------------------------------- end the "set-up-node-list"-part -------------------------------------------------------------
public:
//! ----------------------------------- the public member functions ------------------------------------------------------------------
/** \brief Yield the computed block structure */
void getBlockStructure(std::vector<std::vector<unsigned int> >& blockStructure) const
{
// The total number of blocks
blockStructure.resize(get_number_of_blocks());
// Fill each block
int idx = 0;
for (size_t i=0; i<get_number_of_blocks(); i++) {
blockStructure[i].resize(get_block_size(i));
for (size_t j=0; j<get_block_size(i); j++)
blockStructure[i][j] = get_inverse_permuted_index(idx++);
}
}
//set up the permutation
void set_up_permutation()
{
AxisComparison<CoordinateType, dim> axisComparison(axis_of_anisotropy_);
std::sort(node_list_.begin(), node_list_.end(), axisComparison);
permutation_set_up_ = true;
int maxSize = node_list_.size();
for( int permuted_index = 0; permuted_index < maxSize; ++permuted_index )
{
int original_index = node_list_[ permuted_index ].get_node_index();
original_node_list_[ original_index ].set_permuted_node_index( permuted_index );
}
}
int get_inverse_permuted_index( int permuted_index ) const
{
return node_list_[ permuted_index ].get_node_index();
}
int get_permuted_index( int index ) //inverse_
{
return original_node_list_[ index ].get_permuted_node_index();
}
// set up the finding of the sizes of the diagonal blocks
void set_up_blocks()
{
if ( permutation_set_up_ == true )
{
int max_size = node_list_.size();
int block[max_size];
int block_counter = -1;
// Initialize to a coordinate that doesn't occur in the grid
CoordinateType buffer_coordinate(std::numeric_limits<double>::infinity());
for( int i = 0; i < max_size ; ++i )
{
bool coordinates_partially_equal = true;
for ( int d = 0; d < dim ; ++d )
{
if (d!=axis_of_anisotropy_ and buffer_coordinate[d] != node_list_[i].get_global_coordinate()[d] )
coordinates_partially_equal = false;
}
if ( coordinates_partially_equal == true )
{
block[block_counter] += 1;
}
else
{
for ( int d = 0; d < dim ; ++d )
buffer_coordinate[d] = node_list_[i].get_global_coordinate()[d];
block_counter += 1;
block[block_counter] = 1;
}
}
number_of_blocks_ = block_counter + 1;
block_sizes = new int[ number_of_blocks_ ];
for (int b=0; b<number_of_blocks_; b++)
block_sizes[b] = block[b];
blocks_set_up_ = true;
}
else {
std :: cout << "The permutation has not yet been set up! Use set_up() before set_up_blocks()." << std :: endl;
abort();
}
}
int get_number_of_blocks() const
{
return number_of_blocks_;
}
// get the size of the block_number'th diagonal block
int get_block_size( int block_number) const
{
if ( blocks_set_up_ == true )
{
if ( block_number < number_of_blocks_ )
{ return block_sizes[block_number]; }
else
{ std :: cout << "Request for block number " << block_number << ", but there are only " << number_of_blocks_ << " blocks." << std :: endl;
std :: cout << "Note that the numeration of the blocks starts with 0." << std :: endl;
abort();
}
}
else
{ std :: cout << "The blocks have not yet been set up! Use set_up_blocks() before get_block_size(int block_number )." << std :: endl; }
return 0;
}
void set_up_all()
{
set_up_permutation();
set_up_blocks();
}
void test_permutation( )
{
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_.size(); ++i )
{
std :: cout << "----------------------------------" << std :: endl;
std :: cout << "i = " << i << std :: endl;
std :: cout << "get_inverse_permuted_index(i) = " << get_inverse_permuted_index(i) << std :: endl;
std :: cout << "get_permuted_index(i) = " << get_permuted_index(i) << std :: endl;
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 identity." << std :: endl; }
if (!(i == get_inverse_permuted_index(get_permuted_index(i)) ))
{ std :: cout << "Error! Inverse permutation of permutation is not equal to identity." << std :: endl; }
}
}
void test_permutation_silently( )
{
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_.size(); ++i )
{
if (!(i == get_permuted_index(get_inverse_permuted_index(i)) ))
{ 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 identity." << std :: endl; } }
}
//set up the permutation
void print_information()
{
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_.size(); ++i )
{
std::cout << "Original node[" << get_inverse_permuted_index( i ) << "] = (" << node_list_[i].get_global_coordinate() << ") recieved permuted index " << i << std::endl;
// alternative:
//std::cout << "Original node[" << get_inverse_permuted_index( i ) << "] = (" << original_node_list_[get_inverse_permuted_index( i )].get_global_coordinate() << ") recieved permuted index " << i << std::endl;
}
#if 0
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;
}
#endif
if (permutation_set_up_ == false)
{ std :: cout << "Note that the permutation has not yet been set up. Therefore the output above is still default and the permutation equals the identity." << std :: endl;}
}
};
#endif