Commit 71cbbcea authored by Michael Kuron's avatar Michael Kuron
Browse files

First attempt at allowing PE to run at double precision while LB runs at float

parent 111b212c
......@@ -61,6 +61,7 @@ include( CTest )
# Build options
option ( WALBERLA_DOUBLE_ACCURACY "Floating point accuracy, defaults to double" ON )
option ( WALBERLA_PE_DOUBLE_ACCURACY "PE Floating point accuracy, defaults to double" ON )
option ( WALBERLA_ENABLE_GUI "Compile with GUI" )
option ( WALBERLA_BUILD_TESTS "Build Testcases" )
......
......@@ -22,6 +22,7 @@
#pragma once
#include "waLBerlaDefinitions.h"
#include "core/math/MathTrait.h"
#include <boost/cstdint.hpp>
#include <boost/make_shared.hpp>
......@@ -171,6 +172,11 @@ typedef double real_t;
#else
typedef float real_t;
#endif
#ifdef WALBERLA_PE_DOUBLE_ACCURACY
typedef double pe_real_t;
#else
typedef float pe_real_t;
#endif
template< typename T > inline real_t real_c ( T t ) { return numeric_cast< real_t >(t); } ///< cast to type real_t using "real_c(x)"
template< typename T > inline double double_c( T t ) { return numeric_cast< double >(t); } ///< cast to type double
......@@ -214,20 +220,11 @@ inline bool realIsEqual( const real_t a, const real_t b, const real_t eps = real
return std::fabs( a - b ) < eps;
}
inline bool floatIsEqual( long double lhs, long double rhs, const long double epsilon = real_comparison::Epsilon<long double>::value )
{
return std::fabs( lhs - rhs ) < epsilon;
}
inline bool floatIsEqual( double lhs, double rhs, const double epsilon = real_comparison::Epsilon<double>::value )
{
return std::fabs( lhs - rhs ) < epsilon;
}
inline bool floatIsEqual( float lhs, float rhs, const float epsilon = real_comparison::Epsilon<float>::value )
template< typename T, typename U >
inline bool floatIsEqual( T lhs, U rhs, const typename math::MathTrait<T,U>::LowType epsilon = real_comparison::Epsilon<typename math::MathTrait<T,U>::LowType>::value )
{
return std::fabs( lhs - rhs ) < epsilon;
typedef typename math::MathTrait<T,U>::LowType LT;
return std::fabs( LT(lhs) - LT(rhs) ) < epsilon;
}
......
......@@ -1745,16 +1745,16 @@ namespace walberla {
namespace debug {
namespace check_functions_detail {
template< >
inline bool check_float_equal( const math::Matrix3<real_t> & lhs, const math::Matrix3<real_t> & rhs )
template< typename T, typename U >
inline bool check_float_equal( const math::Matrix3<T> & lhs, const math::Matrix3<U> & rhs )
{
return floatIsEqual( lhs[0], rhs[0] ) && floatIsEqual( lhs[1], rhs[1] ) && floatIsEqual( lhs[2], rhs[2] )
&& floatIsEqual( lhs[3], rhs[3] ) && floatIsEqual( lhs[4], rhs[4] ) && floatIsEqual( lhs[5], rhs[5] )
&& floatIsEqual( lhs[6], rhs[6] ) && floatIsEqual( lhs[7], rhs[7] ) && floatIsEqual( lhs[8], rhs[8] );
}
template< >
inline bool check_float_equal_eps( const math::Matrix3<real_t> & lhs, const math::Matrix3<real_t> & rhs, const real_t epsilon )
template< typename T, typename U >
inline bool check_float_equal_eps( const math::Matrix3<T> & lhs, const math::Matrix3<U> & rhs, const typename math::MathTrait<T,U>::LowType epsilon )
{
return floatIsEqual( lhs[0], rhs[0], epsilon ) && floatIsEqual( lhs[1], rhs[1], epsilon ) && floatIsEqual( lhs[2], rhs[2], epsilon )
&& floatIsEqual( lhs[3], rhs[3], epsilon ) && floatIsEqual( lhs[4], rhs[4], epsilon ) && floatIsEqual( lhs[5], rhs[5], epsilon )
......
......@@ -115,7 +115,7 @@ void FunctionParser::parse( const std::string & eq )
if (isConstant_)
{
const double value = evaluate(std::map<std::string, double>());
isZero_ = floatIsEqual(value, 0);
isZero_ = floatIsEqual(value, 0.0);
}
else
{
......
......@@ -271,7 +271,7 @@ inline Quaternion<Type>::Quaternion( Vector3<Axis> axis, Type angle )
static_assert(boost::is_floating_point<Axis>::value, "Axis has to be floating point!" );
auto axisLength = axis.length();
if( (floatIsEqual(axisLength, 0)) || (math::equal(std::fabs(angle), real_c(0))) ) {
if( (floatIsEqual(axisLength, 0.0)) || (math::equal(std::fabs(angle), real_c(0))) ) {
reset();
return;
}
......@@ -1192,14 +1192,14 @@ namespace walberla {
namespace debug {
namespace check_functions_detail {
template< >
inline bool check_float_equal( const math::Quaternion<real_t> & lhs, const math::Quaternion<real_t> & rhs )
template< typename T, typename U >
inline bool check_float_equal( const math::Quaternion<T> & lhs, const math::Quaternion<U> & rhs )
{
return floatIsEqual( lhs[0], rhs[0] ) && floatIsEqual( lhs[1], rhs[1] ) && floatIsEqual( lhs[2], rhs[2] ) && floatIsEqual( lhs[3], rhs[3] );
}
template< >
inline bool check_float_equal_eps( const math::Quaternion<real_t> & lhs, const math::Quaternion<real_t> & rhs, const real_t epsilon )
template< typename T, typename U >
inline bool check_float_equal_eps( const math::Quaternion<T> & lhs, const math::Quaternion<U> & rhs, const typename math::MathTrait<T,U>::LowType epsilon )
{
return floatIsEqual( lhs[0], rhs[0], epsilon ) && floatIsEqual( lhs[1], rhs[1], epsilon ) && floatIsEqual( lhs[2], rhs[2], epsilon ) && floatIsEqual( lhs[3], rhs[3], epsilon );
}
......
......@@ -1944,14 +1944,14 @@ namespace walberla {
namespace debug {
namespace check_functions_detail {
template< >
inline bool check_float_equal( const math::Vector3<real_t> & lhs, const math::Vector3<real_t> & rhs )
template< typename T, typename U >
inline bool check_float_equal( const math::Vector3<T> & lhs, const math::Vector3<U> & rhs )
{
return floatIsEqual( lhs[0], rhs[0] ) && floatIsEqual( lhs[1], rhs[1] ) && floatIsEqual( lhs[2], rhs[2] );
}
template< >
inline bool check_float_equal_eps( const math::Vector3<real_t> & lhs, const math::Vector3<real_t> & rhs, const real_t epsilon )
template< typename T, typename U >
inline bool check_float_equal_eps( const math::Vector3<T> & lhs, const math::Vector3<U> & rhs, const typename math::MathTrait<T,U>::LowType epsilon )
{
return floatIsEqual( lhs[0], rhs[0], epsilon ) && floatIsEqual( lhs[1], rhs[1], epsilon ) && floatIsEqual( lhs[2], rhs[2], epsilon );
}
......
......@@ -61,32 +61,33 @@ namespace geometry {
* \image html lineBox.png
* \image latex lineBox.eps "Calculation of the closest points between a line and a box" width=455pt
*/
void getClosestLineBoxPoints( const Vector3<real_t>& p1, const Vector3<real_t>& p2, const Vector3<real_t>& c, const Matrix3<real_t>& R,
const Vector3<real_t>& side, Vector3<real_t>& lret, Vector3<real_t>& bret )
template <typename T>
void getClosestLineBoxPoints( const Vector3<T>& p1, const Vector3<T>& p2, const Vector3<T>& c, const Matrix3<T>& R,
const Vector3<T>& side, Vector3<T>& lret, Vector3<T>& bret )
{
using namespace walberla::math;
//----- Note: All computations will be performed in the box-relative coordinate-system -----
// Calculating the global and relative direction of the line p1--p2
const Vector3<real_t> l( p2 - p1 );
Vector3<real_t> tmp( R * l );
const Vector3<T> l( p2 - p1 );
Vector3<T> tmp( R * l );
// Saving the sign of the direction p1--p2
const real_t sign[] = { math::sign(tmp[0]), math::sign(tmp[1]), math::sign(tmp[2]) };
const T sign[] = { math::sign(tmp[0]), math::sign(tmp[1]), math::sign(tmp[2]) };
// Calculating the absolute values of the direction direction
const real_t v [] = { sign[0]*tmp[0], sign[1]*tmp[1], sign[2]*tmp[2] };
const real_t v2[] = { sq( v[0] ) , sq( v[1] ) , sq( v[2] ) };
const T v [] = { sign[0]*tmp[0], sign[1]*tmp[1], sign[2]*tmp[2] };
const T v2[] = { sq( v[0] ) , sq( v[1] ) , sq( v[2] ) };
// Calculating the starting point of the line p1--p2 in box-relative coordinates
tmp = p1 - c;
const real_t s[] = { sign[0]*( R[0]*tmp[0] + R[3]*tmp[1] + R[6]*tmp[2] ),
const T s[] = { sign[0]*( R[0]*tmp[0] + R[3]*tmp[1] + R[6]*tmp[2] ),
sign[1]*( R[1]*tmp[0] + R[4]*tmp[1] + R[7]*tmp[2] ),
sign[2]*( R[2]*tmp[0] + R[5]*tmp[1] + R[8]*tmp[2] ) };
// Calculating the half lengths of the box
const real_t h[] = { real_t(0.5)*side[0], real_t(0.5)*side[1], real_t(0.5)*side[2] };
const T h[] = { T(0.5)*side[0], T(0.5)*side[1], T(0.5)*side[2] };
// Estimating the region of the starting point depending on which side of the
......@@ -95,10 +96,10 @@ void getClosestLineBoxPoints( const Vector3<real_t>& p1, const Vector3<real_t>&
// are no more. Additionally, d|d|^2/dt is computed for t=0. If it is >= 0
// then p1 is the closest point since the line points away from the box.
int region [] = { 0, 0, 0 };
real_t tanchor[] = { 2, 2, 2 }; // Invalid t values; t cannot be greater than 1
real_t dd2dt( 0 );
T tanchor[] = { 2, 2, 2 }; // Invalid t values; t cannot be greater than 1
T dd2dt( 0 );
if( v[0] > Limits<real_t>::epsilon() )
if( v[0] > Limits<T>::epsilon() )
{
if( s[0] < -h[0] ) {
region[0] = -1;
......@@ -115,7 +116,7 @@ void getClosestLineBoxPoints( const Vector3<real_t>& p1, const Vector3<real_t>&
}
}
if( v[1] > Limits<real_t>::epsilon() )
if( v[1] > Limits<T>::epsilon() )
{
if( s[1] < -h[1] ) {
region[1] = -1;
......@@ -132,7 +133,7 @@ void getClosestLineBoxPoints( const Vector3<real_t>& p1, const Vector3<real_t>&
}
}
if( v[2] > Limits<real_t>::epsilon() )
if( v[2] > Limits<T>::epsilon() )
{
if( s[2] < -h[2] ) {
region[2] = -1;
......@@ -151,26 +152,26 @@ void getClosestLineBoxPoints( const Vector3<real_t>& p1, const Vector3<real_t>&
// Calculating the value t for the closest point on the line
real_t t( 0 );
T t( 0 );
if( dd2dt < -Limits<real_t>::accuracy() )
if( dd2dt < -Limits<T>::accuracy() )
{
do {
// Finding the t value for the next clip plane/line intersection
real_t next_t( 1 );
T next_t( 1 );
if( ( tanchor[0] > t ) && ( tanchor[0] < real_t(1) ) && ( tanchor[0] < next_t ) ) {
if( ( tanchor[0] > t ) && ( tanchor[0] < T(1) ) && ( tanchor[0] < next_t ) ) {
next_t = tanchor[0];
}
if( ( tanchor[1] > t ) && ( tanchor[1] < real_t(1) ) && ( tanchor[1] < next_t ) ) {
if( ( tanchor[1] > t ) && ( tanchor[1] < T(1) ) && ( tanchor[1] < next_t ) ) {
next_t = tanchor[1];
}
if( ( tanchor[2] > t ) && ( tanchor[2] < real_t(1) ) && ( tanchor[2] < next_t ) ) {
if( ( tanchor[2] > t ) && ( tanchor[2] < T(1) ) && ( tanchor[2] < next_t ) ) {
next_t = tanchor[2];
}
// Computing d|d|^2/dt for the next t
real_t next_dd2dt( 0 );
T next_dd2dt( 0 );
if( region[0] != 0 ) {
next_dd2dt += v2[0] * ( next_t - tanchor[0] );
......@@ -205,11 +206,11 @@ void getClosestLineBoxPoints( const Vector3<real_t>& p1, const Vector3<real_t>&
t = next_t;
dd2dt = next_dd2dt;
}
while( t < real_t(1) );
while( t < T(1) );
}
WALBERLA_ASSERT_GREATER_EQUAL( t, real_t(0), "Invalid line point" );
WALBERLA_ASSERT_LESS_EQUAL( t, real_t(1), "Invalid line point" );
WALBERLA_ASSERT_GREATER_EQUAL( t, T(0), "Invalid line point" );
WALBERLA_ASSERT_LESS_EQUAL( t, T(1), "Invalid line point" );
// Computing the closest point on the line
lret = p1 + t * l;
......@@ -229,6 +230,12 @@ void getClosestLineBoxPoints( const Vector3<real_t>& p1, const Vector3<real_t>&
bret = ( R * tmp ) + c;
}
template
void getClosestLineBoxPoints<float>( const Vector3<float>& p1, const Vector3<float>& p2, const Vector3<float>& c, const Matrix3<float>& R,
const Vector3<float>& side, Vector3<float>& lret, Vector3<float>& bret );
template
void getClosestLineBoxPoints<double>( const Vector3<double>& p1, const Vector3<double>& p2, const Vector3<double>& c, const Matrix3<double>& R,
const Vector3<double>& side, Vector3<double>& lret, Vector3<double>& bret );
//*************************************************************************************************
......@@ -248,45 +255,46 @@ void getClosestLineBoxPoints( const Vector3<real_t>& p1, const Vector3<real_t>&
* involving the endpoint of at least one line will be returned. This will also work correctly
* for zero length lines, e.g. if \a a1 = \a a2 and/or \a b1 = \a b2.
*/
void getClosestLineSegmentPoints( const Vector3<real_t>& a1, const Vector3<real_t>& a2, const Vector3<real_t>& b1, const Vector3<real_t>& b2,
Vector3<real_t>& cp1, Vector3<real_t>& cp2 )
template <typename T>
void getClosestLineSegmentPoints( const Vector3<T>& a1, const Vector3<T>& a2, const Vector3<T>& b1, const Vector3<T>& b2,
Vector3<T>& cp1, Vector3<T>& cp2 )
{
using namespace walberla::math;
//----- Checking the vertex-vertex features -----
const Vector3<real_t> a1a2( a2 - a1 );
const Vector3<real_t> b1b2( b2 - b1 );
const Vector3<real_t> a1b1( b1 - a1 );
const real_t da1 ( a1a2 * a1b1 );
const real_t db1 ( b1b2 * a1b1 );
if( da1 <= +Limits<real_t>::fpuAccuracy() && db1 >= -Limits<real_t>::fpuAccuracy() ) {
const Vector3<T> a1a2( a2 - a1 );
const Vector3<T> b1b2( b2 - b1 );
const Vector3<T> a1b1( b1 - a1 );
const T da1 ( a1a2 * a1b1 );
const T db1 ( b1b2 * a1b1 );
if( da1 <= +Limits<T>::fpuAccuracy() && db1 >= -Limits<T>::fpuAccuracy() ) {
cp1 = a1;
cp2 = b1;
return;
}
const Vector3<real_t> a1b2( b2 - a1 );
const real_t da2 ( a1a2 * a1b2 );
const real_t db2 ( b1b2 * a1b2 );
if( da2 <= +Limits<real_t>::fpuAccuracy() && db2 <= +Limits<real_t>::fpuAccuracy() ) {
const Vector3<T> a1b2( b2 - a1 );
const T da2 ( a1a2 * a1b2 );
const T db2 ( b1b2 * a1b2 );
if( da2 <= +Limits<T>::fpuAccuracy() && db2 <= +Limits<T>::fpuAccuracy() ) {
cp1 = a1;
cp2 = b2;
return;
}
const Vector3<real_t> a2b1( b1 - a2 );
const real_t da3 ( a1a2 * a2b1 );
const real_t db3 ( b1b2 * a2b1 );
if( da3 >= -Limits<real_t>::fpuAccuracy() && db3 >= -Limits<real_t>::fpuAccuracy() ) {
const Vector3<T> a2b1( b1 - a2 );
const T da3 ( a1a2 * a2b1 );
const T db3 ( b1b2 * a2b1 );
if( da3 >= -Limits<T>::fpuAccuracy() && db3 >= -Limits<T>::fpuAccuracy() ) {
cp1 = a2;
cp2 = b1;
return;
}
const Vector3<real_t> a2b2( b2 - a2 );
const real_t da4 ( a1a2 * a2b2 );
const real_t db4 ( b1b2 * a2b2 );
if( da4 >= -Limits<real_t>::fpuAccuracy() && db4 <= +Limits<real_t>::fpuAccuracy() ) {
const Vector3<T> a2b2( b2 - a2 );
const T da4 ( a1a2 * a2b2 );
const T db4 ( b1b2 * a2b2 );
if( da4 >= -Limits<T>::fpuAccuracy() && db4 <= +Limits<T>::fpuAccuracy() ) {
cp1 = a2;
cp2 = b2;
return;
......@@ -297,44 +305,44 @@ void getClosestLineSegmentPoints( const Vector3<real_t>& a1, const Vector3<real_
// If one or both of the line segments have zero length, we will never get here. Therefore
// we don't have to worry about possible divisions by zero in the following calculations.
Vector3<real_t> n, k;
Vector3<T> n, k;
const real_t la( a1a2 * a1a2 );
if( da1 >= -Limits<real_t>::fpuAccuracy() && da3 <= +Limits<real_t>::fpuAccuracy() ) {
const T la( a1a2 * a1a2 );
if( da1 >= -Limits<T>::fpuAccuracy() && da3 <= +Limits<T>::fpuAccuracy() ) {
k = (da1 / la) * a1a2;
n = a1b1 - k;
if( b1b2 * n >= -Limits<real_t>::fpuAccuracy() ) {
if( b1b2 * n >= -Limits<T>::fpuAccuracy() ) {
cp1 = a1 + k;
cp2 = b1;
return;
}
}
if( da2 >= -Limits<real_t>::fpuAccuracy() && da4 <= +Limits<real_t>::fpuAccuracy() ) {
if( da2 >= -Limits<T>::fpuAccuracy() && da4 <= +Limits<T>::fpuAccuracy() ) {
k = (da2 / la) * a1a2;
n = a1b2 - k;
if( b1b2 * n <= +Limits<real_t>::fpuAccuracy() ) {
if( b1b2 * n <= +Limits<T>::fpuAccuracy() ) {
cp1 = a1 + k;
cp2 = b2;
return;
}
}
const real_t lb( b1b2 * b1b2 );
if( db1 <= +Limits<real_t>::fpuAccuracy() && db2 >= -Limits<real_t>::fpuAccuracy() ) {
const T lb( b1b2 * b1b2 );
if( db1 <= +Limits<T>::fpuAccuracy() && db2 >= -Limits<T>::fpuAccuracy() ) {
k = (-db1 / lb) * b1b2;
n = -a1b1 - k;
if( a1a2 * n >= -Limits<real_t>::fpuAccuracy() ) {
if( a1a2 * n >= -Limits<T>::fpuAccuracy() ) {
cp1 = a1;
cp2= b1 + k;
return;
}
}
if( db3 <= +Limits<real_t>::fpuAccuracy() && db4 >= -Limits<real_t>::fpuAccuracy() ) {
if( db3 <= +Limits<T>::fpuAccuracy() && db4 >= -Limits<T>::fpuAccuracy() ) {
k = (-db3 / lb) * b1b2;
n = -a2b1 - k;
if( a1a2 * n <= +Limits<real_t>::fpuAccuracy() ) {
if( a1a2 * n <= +Limits<T>::fpuAccuracy() ) {
cp1 = a2;
cp2 = b1 + k;
return;
......@@ -344,18 +352,24 @@ void getClosestLineSegmentPoints( const Vector3<real_t>& a1, const Vector3<real_
//----- Checking the edge-edge features -----
const real_t scale( a1a2 * b1b2 );
real_t div( la * lb - math::sq(scale) );
const T scale( a1a2 * b1b2 );
T div( la * lb - math::sq(scale) );
WALBERLA_ASSERT_GREATER( div, real_t(0), std::setprecision(16) << "Invalid division\n" << a1 << "\n" << a2 << "\n" << b1 << "\n" << b2 );
div = real_t(1) / div;
WALBERLA_ASSERT_GREATER( div, T(0), std::setprecision(16) << "Invalid division\n" << a1 << "\n" << a2 << "\n" << b1 << "\n" << b2 );
div = T(1) / div;
const real_t s( ( lb * da1 - scale * db1 ) * div );
const real_t t( ( scale * da1 - la * db1 ) * div );
const T s( ( lb * da1 - scale * db1 ) * div );
const T t( ( scale * da1 - la * db1 ) * div );
cp1 = a1 + s * a1a2;
cp2 = b1 + t * b1b2;
}
template
void getClosestLineSegmentPoints<float>( const Vector3<float>& a1, const Vector3<float>& a2, const Vector3<float>& b1, const Vector3<float>& b2,
Vector3<float>& cp1, Vector3<float>& cp2 );
template
void getClosestLineSegmentPoints<double>( const Vector3<double>& a1, const Vector3<double>& a2, const Vector3<double>& b1, const Vector3<double>& b2,
Vector3<double>& cp1, Vector3<double>& cp2 );
//*************************************************************************************************
......@@ -377,30 +391,37 @@ void getClosestLineSegmentPoints( const Vector3<real_t>& a1, const Vector3<real_
\f[ s = \frac{\det(o_2-o_1,d_2,d_1 \times d_2)}{\left \Vert d_1 \times d_2 \right \| ^2 } \f]
\f[ t = \frac{\det(o_2-o_1,d_1,d_1 \times d_2)}{\left \Vert d_1 \times d_2 \right \| ^2 } \f]
*/
void intersectLines( const Vector3<real_t>& o1, const Vector3<real_t>& d1, const Vector3<real_t>& o2, const Vector3<real_t>& d2,
real_t& s, real_t& t )
template <typename T>
void intersectLines( const Vector3<T>& o1, const Vector3<T>& d1, const Vector3<T>& o2, const Vector3<T>& d2,
T& s, T& t )
{
using namespace walberla::math;
const real_t sqrlen( ( d1 % d2 ).sqrLength() );
const T sqrlen( ( d1 % d2 ).sqrLength() );
if( floatIsEqual(sqrlen, 0) )
if( floatIsEqual(sqrlen, 0.0) )
{
s = real_t(0),
t = real_t(0);
s = T(0),
t = T(0);
}
else
{
const real_t isqrlen( real_t(1) / sqrlen );
const Vector3<real_t> p( o2 - o1 );
const real_t dot( d1 * d2 );
const real_t a ( d1 * p );
const real_t b ( -d2 * p );
const T isqrlen( T(1) / sqrlen );
const Vector3<T> p( o2 - o1 );
const T dot( d1 * d2 );
const T a ( d1 * p );
const T b ( -d2 * p );
s = ( a * ( d2 * d2 ) + dot*b ) * isqrlen;
t = ( b * ( d1 * d1 ) + dot*a ) * isqrlen;
}
}
template
void intersectLines<float>( const Vector3<float>& o1, const Vector3<float>& d1, const Vector3<float>& o2, const Vector3<float>& d2,
float& s, float& t );
template
void intersectLines<double>( const Vector3<double>& o1, const Vector3<double>& d1, const Vector3<double>& o2, const Vector3<double>& d2,
double& s, double& t );
//*************************************************************************************************
} // namespace math
......
......@@ -47,14 +47,17 @@ namespace geometry {
//*************************************************************************************************
/*!\name Geometry functions */
//@{
void getClosestLineBoxPoints( const Vector3<real_t>& p1, const Vector3<real_t>& p2, const Vector3<real_t>& c, const Matrix3<real_t>& R,
const Vector3<real_t>& side, Vector3<real_t>& lret, Vector3<real_t>& bret);
template <typename T>
void getClosestLineBoxPoints( const Vector3<T>& p1, const Vector3<T>& p2, const Vector3<T>& c, const Matrix3<T>& R,
const Vector3<T>& side, Vector3<T>& lret, Vector3<T>& bret);
void getClosestLineSegmentPoints( const Vector3<real_t>& a1, const Vector3<real_t>& a2, const Vector3<real_t>& b1, const Vector3<real_t>& b2,
Vector3<real_t>& cp1, Vector3<real_t>& cp2 );
template <typename T>
void getClosestLineSegmentPoints( const Vector3<T>& a1, const Vector3<T>& a2, const Vector3<T>& b1, const Vector3<T>& b2,
Vector3<T>& cp1, Vector3<T>& cp2 );
void intersectLines( const Vector3<real_t>& o1, const Vector3<real_t>& d1, const Vector3<real_t>& o2, const Vector3<real_t>& d2,
real_t& s, real_t& t );
template <typename T>
void intersectLines( const Vector3<T>& o1, const Vector3<T>& d1, const Vector3<T>& o2, const Vector3<T>& d2,
T& s, T& t );
inline bool originInTetrahedron( const Vector3<real_t>& A, const Vector3<real_t>& B, const Vector3<real_t>& C, const Vector3<real_t>& D );
inline bool pointInTetrahedron ( const Vector3<real_t>& A, const Vector3<real_t>& B, const Vector3<real_t>& C, const Vector3<real_t>& D,
......
......@@ -47,10 +47,11 @@ bool hasNeighborOwner(const BlockT& block, const Owner& owner)
* Returns -1 if no valid block is found otherwise the process rank of the containing block is returned.
*/
template <class BlockT>
Owner findContainingProcess(const BlockT& block, math::Vector3<real_t> pt)
Owner findContainingProcess(const BlockT& block, math::Vector3<pe_real_t> pt)
{
if (block.getAABB().contains(pt)) return Owner(int_c(block.getProcess()), block.getId().getID());
if (!block.getBlockStorage().getDomain().contains(pt)) block.getBlockStorage().mapToPeriodicDomain(pt);
Vector3<real_t> pr(real_c(pt[0]), real_c(pt[1]), real_c(pt[2]));
if (!block.getBlockStorage().getDomain().contains(pr)) block.getBlockStorage().mapToPeriodicDomain(pr);
for( uint_t i = uint_t(0); i != block.getNeighborhoodSize(); ++i )
{
if (block.getNeighborAABB(i).contains(pt)) return Owner(int_c(block.getNeighborProcess(i)), block.getNeighborId(i).getID());
......
......@@ -5,7 +5,11 @@
#
###################################################################################################
if(WALBERLA_DOUBLE_ACCURACY)
if (WALBERLA_DOUBLE_ACCURACY AND NOT WALBERLA_PE_DOUBLE_ACCURACY)
message(FATAL_ERROR "PE cannot use double precision if waLBerla uses single precision")
endif()
if(WALBERLA_PE_DOUBLE_ACCURACY)
set(CCD_DOUBLE ON)
set(CCD_SINGLE OFF)
else()
......
......@@ -44,7 +44,7 @@ namespace pe {
* For more details about the calculation of the motion of a rigid body, see the description of
* the RigidBody::calcMotion() function.
*/
const real_t sleepThreshold = real_c( 0 );
const pe_real_t sleepThreshold = real_c( 0 );
//*************************************************************************************************
......@@ -57,7 +57,7 @@ const real_t sleepThreshold = real_c( 0 );
* recency-weighted average equal to the current motion of the rigid body, a bias of one
* ignores the current motion altogether.
*/
const real_t sleepBias = real_c( 0.5 );
const pe_real_t sleepBias = real_c( 0.5 );
//*************************************************************************************************
} // namespace pe
......
......@@ -72,34 +72,34 @@ bool Material::activateMaterials()
// Initializing the coefficients of restitution
// | Iron | Copper | Granite | Oak | Fir |
const real_t cor[5][5] = {
{ static_cast<real_t>(0.25), static_cast<real_t>(0.25), static_cast<real_t>(0.25), static_cast<real_t>(0.25), static_cast<real_t>(0.25) }, // Iron
{ static_cast<real_t>(0.25), static_cast<real_t>(0.25), static_cast<real_t>(0.25), static_cast<real_t>(0.25), static_cast<real_t>(0.25) }, // Copper
{ static_cast<real_t>(0.25), static_cast<real_t>(0.25), static_cast<real_t>(0.25), static_cast<real_t>(0.25), static_cast<real_t>(0.25) }, // Granite
{ static_cast<real_t>(0.25), static_cast<real_t>(0.25), static_cast<real_t>(0.25), static_cast<real_t>(0.25), static_cast<real_t>(0.25) }, // Oak
{ static_cast<real_t>(0.25), static_cast<real_t>(0.25), static_cast<real_t>(0.25), static_cast<real_t>(0.25), static_cast<real_t>(0.25) } // Fir
const pe_real_t cor[5][5] = {
{ static_cast<pe_real_t>(0.25), static_cast<pe_real_t>(0.25), static_cast<pe_real_t>(0.25), static_cast<pe_real_t>(0.25), static_cast<pe_real_t>(0.25) }, // Iron
{ static_cast<pe_real_t>(0.25), static_cast<pe_real_t>(0.25), static_cast<pe_real_t>(0.25), static_cast<pe_real_t>(0.25), static_cast<pe_real_t>(0.25) }, // Copper
{ static_cast<pe_real_t>(0.25), static_cast<pe_real_t>(0.25), static_cast<pe_real_t>(0.25), static_cast<pe_real_t>(0.25), static_cast<pe_real_t>(0.25) }, // Granite
{ static_cast<pe_real_t>(0.25), static_cast<pe_real_t>(0.25), static_cast<pe_real_t>(0.25), static_cast<