-
Oliver Sander authored
Otherwise the code won't compile when the two dimensions differ.
Oliver Sander authoredOtherwise the code won't compile when the two dimensions differ.
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