Commit 7056169f authored by Sebastian Eichelbaum's avatar Sebastian Eichelbaum
Browse files

[MERGE]

parents b5d26e0a 14e50187
## This is a sample configuration file for OpenWalnut.
## Uncomment the options you are interested in.
[general]
allowOnlyOneFiberDataSet = yes # This will prevent you from accidently loading multiple fiber data sets.
[modules]
## use this to specify the default module to add during load.
## It is a comma seperated list. If this is not specified the default empty is assumed.
......
//---------------------------------------------------------------------------
//
// Project: OpenWalnut ( http://www.openwalnut.org )
//
// Copyright 2009 OpenWalnut Community, BSV@Uni-Leipzig and CNCF@MPI-CBS
// For more information see http://www.openwalnut.org/copying
//
// This file is part of OpenWalnut.
//
// OpenWalnut is free software: you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// OpenWalnut is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>.
//
//---------------------------------------------------------------------------
#include "WBinom.h"
//---------------------------------------------------------------------------
//
// Project: OpenWalnut ( http://www.openwalnut.org )
//
// Copyright 2009 OpenWalnut Community, BSV@Uni-Leipzig and CNCF@MPI-CBS
// For more information see http://www.openwalnut.org/copying
//
// This file is part of OpenWalnut.
//
// OpenWalnut is free software: you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// OpenWalnut is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>.
//
//---------------------------------------------------------------------------
#ifndef WBINOM_H
#define WBINOM_H
#include <string>
namespace wmath
{
/**
* Implements compile-time calculation of binomial coefficients.
*
*
* WBinom< n, k >::value = n! / ( k!(n-k)! ).
*
* \note For k > n, compilation fails.
*/
template< std::size_t n, std::size_t k >
struct WBinom
{
enum
{
value = WBinom< n - 1, k - 1 >::value + WBinom< n - 1, k >::value
};
};
/**
* Specialization for n = k.
*/
template< std::size_t n >
struct WBinom< n, n >
{
enum
{
value = 1
};
};
/**
* Specialization for k = 0.
*/
template< std::size_t n >
struct WBinom< n, 0 >
{
enum
{
value = 1
};
};
/**
* This specialization of the WBinom struct is needed to avoid
* infinite recursion in case of k > n. The compiler should abort
* compilation with an error message.
*/
template< std::size_t k >
struct WBinom< 0, k >
{
};
} // namespace wmath
#endif // WBINOM_H
//---------------------------------------------------------------------------
//
// Project: OpenWalnut ( http://www.openwalnut.org )
//
// Copyright 2009 OpenWalnut Community, BSV@Uni-Leipzig and CNCF@MPI-CBS
// For more information see http://www.openwalnut.org/copying
//
// This file is part of OpenWalnut.
//
// OpenWalnut is free software: you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// OpenWalnut is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>.
//
//---------------------------------------------------------------------------
#include "WTensor.h"
//---------------------------------------------------------------------------
//
// Project: OpenWalnut ( http://www.openwalnut.org )
//
// Copyright 2009 OpenWalnut Community, BSV@Uni-Leipzig and CNCF@MPI-CBS
// For more information see http://www.openwalnut.org/copying
//
// This file is part of OpenWalnut.
//
// OpenWalnut is free software: you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// OpenWalnut is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>.
//
//---------------------------------------------------------------------------
#ifndef WTENSOR_H
#define WTENSOR_H
#include <vector>
#include <iostream>
#include "WTensorBase.h"
namespace wmath
{
// ############################# class WTensor<> #############################################
/**
* Implements a tensor that has the same number of components in every
* direction.
*
* The first template parameter is the order of the tensor.
* The second template parameter is the dimension of the tensor, i.e. the number of components
* in each direction.
* The third template parameter is the datatype of the components, which is double by default.
*
* \note The type Data_T may not throw exceptions on construction, destruction or
* during operator =.
*
* Access to specific elements of the tensor can be achieved in 2 ways:
*
* - operator (), whose number of parameters matches the template parameter order.
* - operator [], whose parameter is either an array or a std::vector of appropriate size.
*
* \note The datatype of the array or std::vector can be any type castable to std::size_t.
* \note There is no bounds checking for the array version of operator [].
* \note Operator () is not supported for orders larger than 6.
*
* Examples:
*
* - Construct a tensor of order 2 and dimension 3 (i.e. a 3x3-Matrix):
*
* wmath::WTensor< 2, 3 > w;
*
* - Change Element (2,0) to 4.0:
*
* w( 2, 0 ) = 4.0;
*
* - Construct a 4D-vector:
*
* wmath::WTensor< 1, 4 > v;
*
* - Access v at position 2:
*
* v( 2 ) = ...;
*
* std::vector< int > i( 1, 2 );
* v[ i ] = ...;
*/
template< std::size_t order, std::size_t dim, typename Data_T = double >
class WTensor : public WTensorFunc< WTensorBase, order, dim, Data_T >
{
};
/**
* Specialization for dim = 0. This essentially prohibits
* instantiation of a 0-dimensional tensor.
*/
template< std::size_t order, typename Data_T >
class WTensor< order, 0, Data_T >
{
};
// ###################### special case WTensor<> of order 0 ##############################
/**
* The order 0 version of WTensor. This essentialy encapsulates a scalar.
*/
template< std::size_t dim, typename Data_T >
class WTensor< 0, dim, Data_T >
{
public:
/**
* Standard constructor.
*/
WTensor();
/**
* Get the dimension of this tensor.
*
* \return The dimension of this tensor.
*/
std::size_t getDimension() const;
/**
* Get the order of this tensor.
*
* \return The order of this tensor.
*/
std::size_t getOrder() const;
/**
* Access via std::vector of indices. Only for compatibility.
*
* \param indices A std::vector of indices.
*
* \return A reference to the value of this tensor.
*/
template< typename Index_T >
Data_T& operator[] ( std::vector< Index_T > const& indices );
/**
* Access via std::vector of indices. Only for compatibility.
*
* \param indices A std::vector of indices.
*
* \return A reference to the value of this tensor.
*/
template< typename Index_T >
Data_T const& operator[] ( std::vector< Index_T > const& indices ) const;
/**
* Access operator.
*
* \return A reference to the value of this tensor.
*/
Data_T& operator() ();
/**
* Access operator.
*
* \return A reference to the value of this tensor.
*/
Data_T const& operator() () const;
private:
/**
* The data element.
*/
Data_T m_data;
};
template< std::size_t dim, typename Data_T >
WTensor< 0, dim, Data_T >::WTensor()
: m_data( Data_T() )
{
}
template< std::size_t dim, typename Data_T >
std::size_t WTensor< 0, dim, Data_T >::getDimension() const
{
return 0;
}
template< std::size_t dim, typename Data_T >
std::size_t WTensor< 0, dim, Data_T >::getOrder() const
{
return 0;
}
template< std::size_t dim, typename Data_T > template< typename Index_T >
Data_T& WTensor< 0, dim, Data_T >::operator[] ( std::vector< Index_T > const& indices )
{
return m_data;
}
template< std::size_t dim, typename Data_T > template< typename Index_T >
Data_T const& WTensor< 0, dim, Data_T >::operator[] ( std::vector< Index_T > const& indices ) const
{
return m_data;
}
template< std::size_t dim, typename Data_T >
Data_T& WTensor< 0, dim, Data_T >::operator() ()
{
return m_data;
}
template< std::size_t dim, typename Data_T >
Data_T const& WTensor< 0, dim, Data_T >::operator() () const
{
return m_data;
}
// ######################## stream output operators #################################
/**
* Write a tensor of order 0 to an output stream.
*
* \param o An output stream.
* \param t A WTensor.
*
* \return The output stream.
*/
template< std::size_t dim, typename Data_T >
std::ostream& operator << ( std::ostream& o, WTensor< 0, dim, Data_T > const& t )
{
o << t() << std::endl;
return o;
}
/**
* Write a tensor of order 1 to an output stream.
*
* \param o An output stream.
* \param t A WTensor.
*
* \return The output stream.
*/
template< std::size_t dim, typename Data_T >
std::ostream& operator << ( std::ostream& o, WTensor< 1, dim, Data_T > const& t )
{
for( std::size_t k = 0; k < dim; ++k )
{
o << t( k ) << " ";
}
o << std::endl;
return o;
}
/**
* Write a tensor of order 2 to an output stream.
*
* \param o An output stream.
* \param t A WTensor.
*
* \return The output stream.
*/
template< std::size_t dim, typename Data_T >
std::ostream& operator << ( std::ostream& o, WTensor< 2, dim, Data_T > const& t )
{
for( std::size_t k = 0; k < dim; ++k )
{
for( std::size_t l = 0; l < dim; ++l )
{
o << t( k, l ) << " ";
}
o << std::endl;
}
return o;
}
} // namespace wmath
#endif // WTENSOR_H
//---------------------------------------------------------------------------
//
// Project: OpenWalnut ( http://www.openwalnut.org )
//
// Copyright 2009 OpenWalnut Community, BSV@Uni-Leipzig and CNCF@MPI-CBS
// For more information see http://www.openwalnut.org/copying
//
// This file is part of OpenWalnut.
//
// OpenWalnut is free software: you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// OpenWalnut is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>.
//
//---------------------------------------------------------------------------
#include "WTensorBase.h"
This diff is collapsed.
//---------------------------------------------------------------------------
//
// Project: OpenWalnut ( http://www.openwalnut.org )
//
// Copyright 2009 OpenWalnut Community, BSV@Uni-Leipzig and CNCF@MPI-CBS
// For more information see http://www.openwalnut.org/copying
//
// This file is part of OpenWalnut.
//
// OpenWalnut is free software: you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// OpenWalnut is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>.
//
//---------------------------------------------------------------------------
#include <vector>
#include <cmath>
#include "../WAssert.h"
#include "../WLimits.h"
#include "WTensorFunctions.h"
namespace wmath
{
template< typename Data_T >
void jacobiEigenvector3D( WTensorSym< 2, 3, Data_T > const& mat,
std::vector< Data_T >* eigenValues,
std::vector< WTensor< 1, 3, Data_T > >* eigenVectors )
{
WTensorSym< 2, 3, Data_T > in = WTensorSym< 2, 3, Data_T >( mat );
WTensor< 2, 3, Data_T > ev;
int iter = 50;
Data_T evp[ 3 ];
Data_T evq[ 3 ];
while( iter >= 0 )
{
int p = 1;
int q = 0;
for( int i = 0; i < 2; ++i )
{
if( fabs( in( 2, i ) ) > fabs( in( p, q ) ) )
{
p = 2;
q = i;
}
}
if( fabs( in( p, q ) ) == 0.0 )
{
for( int i = 0; i < 3; ++i )
{
eigenValues->at( i ) = in( i, i );
for( int j = 0; j < 3; ++j )
{
eigenVectors->at( i )( j ) = ev( i, j );
}
}
return;
}
Data_T r = in( q, q ) - in( p, p );
Data_T o = r / ( 2.0 * in( p, q ) );
Data_T t;
Data_T signofo = ( o < 0.0 ? -1.0 : 1.0 );
if( o * o > wlimits::MAX_DOUBLE )
{
t = signofo / ( 2.0 * fabs( o ) );
}
else
{
t = signofo / ( fabs( o ) + sqrt( o * o + 1.0 ) );
}
Data_T c;
if( t * t < wlimits::DBL_EPS )
{
c = 1.0;
}
else
{
c = 1.0 / sqrt( t * t + 1.0 );
}
Data_T s = c * t;
int k = 0;
while( k == q || k == p )
{
++k;
}
WAssert( k < 3, "" );
Data_T u = ( 1.0 - c ) / s;
Data_T x = in( k, p );
Data_T y = in( k, q );
in( p, k ) = in( k, p ) = x - s * ( y + u * x );
in( q, k ) = in( k, q ) = y + s * ( x - u * y );
x = in( p, p );
y = in( q, q );
in( p, p ) = x - t * in( p, q );
in( q, q ) = y + t * in( p, q );
in( q, p ) = in( p, q ) = 0.0;
for( int l = 0; l < 3; ++l )
{
evp[ l ] = ev( l, p );
evq[ l ] = ev( l, q );
}
for( int l = 0; l < 3; ++l )
{
ev( l, p ) = c * evp[ l ] - s * evq[ l ];
ev( l, q ) = s * evp[ l ] + c * evq[ l ];
}
--iter;
}
WAssert( iter >= 0, "Jacobi eigenvector iteration did not converge." );
}
} // namespace wmath
//---------------------------------------------------------------------------
//
// Project: OpenWalnut ( http://www.openwalnut.org )
//
// Copyright 2009 OpenWalnut Community, BSV@Uni-Leipzig and CNCF@MPI-CBS
// For more information see http://www.openwalnut.org/copying
//
// This file is part of OpenWalnut.
//
// OpenWalnut is free software: you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// OpenWalnut is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public License
// along with OpenWalnut. If not, see <http://www.gnu.org/licenses/>.
//
//---------------------------------------------------------------------------
#ifndef WTENSORFUNCTIONS_H
#define WTENSORFUNCTIONS_H
#include <vector>
#include "WTensorSym.h"
namespace wmath
{
/**
* Compute all eigenvalues as well as the corresponding eigenvectors of a
* symmetric real Matrix.
*
* \pre Data_T must be castable to double.
*
* \param[in] mat A real symmetric matrix.