Commit 077ea068 authored by Steffen Schotthöfer's avatar Steffen Schotthöfer
Browse files

quadrature overhaul. Fixed coordinate trafo, fixed reduced accuracy issue for...

quadrature overhaul. Fixed coordinate trafo, fixed reduced accuracy issue for lookup tables, fixed integration on sphere


Former-commit-id: 2a9d6816
parent 1400f222
/*!
* @file config.h
* @brief Class to compute 1D Gauss Legendre Quadrature
* @author J. Kusch
*
* Disclaimer: This class structure was copied and modifed with open source permission from SU2 v7.0.3 https://su2code.github.io/
*/
#ifndef QGAUSSLEGENDRE1D_H
#define QGAUSSLEGENDRE1D_H
......@@ -15,10 +23,26 @@ class QGaussLegendre1D : public QuadratureBase
QGaussLegendre1D( unsigned quadOrder );
virtual ~QGaussLegendre1D() {}
inline void SetName() override { _name = "Tensorized Gauss-Legendre quadrature."; }
inline void SetName() override { _name = "Tensorized Gauss-Legendre quadrature 1D."; }
inline void SetNq() override { _nq = _order; }
void SetPointsAndWeights() override;
void SetConnectivity() override;
/*! @brief Integrates f(x,y,z) with the quadrature.
* @param double(f)( double x0, double x1, double x2 ) : density function that depends on a three spatial dimensions.
* @returns double result: result of the quadrature rule */
double Integrate( double( f )( double x0, double x1, double x2 ) ) override;
/*! @brief Integrates f(x,y,z) with the quadrature.
* @param double(f)( double my, double phi ) : density function that depends on a spherical coordinates.
* @returns double result: result of the quadrature rule */
double IntegrateSpherical( double( f )( double my, double phi ) ) override;
/*! @brief Integrates vector valued f(x,y,z) with the quadrature. Each dimension is integrated by itself.
* @param : double(f)( double x0, double x1, double x2 ) : density function that depends on a three spatial dimensions.
* @param : len : lenght of vector
* @returns double result: result of the quadrature rule (vector valued) */
std::vector<double> Integrate( std::vector<double>( f )( double x0, double x1, double x2 ), unsigned /* len */ ) override;
};
#endif // QGAUSSLEGENDRE1D_H
......@@ -20,8 +20,6 @@ class QGaussLegendreTensorized : public QuadratureBase
void SetNq() override;
void SetPointsAndWeights() override;
void SetConnectivity() override;
inline VectorVector GetPointsSphere() const override { return _pointsSphere; }
};
#endif // QGAUSSLEGENDRETENSORIZED_H
......@@ -14,7 +14,6 @@ class QMonteCarlo : public QuadratureBase
inline void SetNq() override { _nq = GetOrder(); }
void SetPointsAndWeights() override;
void SetConnectivity() override;
VectorVector GetPointsSphere() const override;
};
#endif // QMONTECARLO_H
......@@ -25,8 +25,6 @@ class QProduct : public QuadratureBase
inline void SetNq() override { _nq = 4 * pow( GetOrder(), 2 ); }
void SetPointsAndWeights() override;
void SetConnectivity() override;
inline VectorVector GetPointsSphere() const override { return _pointsSphere; }
};
#endif // PRODUCTQUADRATURE_H
......@@ -35,13 +35,18 @@ class QuadratureBase
/*! @brief Integrates f(x,y,z) with the quadrature.
* @param double(f)( double x0, double x1, double x2 ) : density function that depends on a three spatial dimensions.
* @returns double result: result of the quadrature rule */
double Integrate( double( f )( double x0, double x1, double x2 ) );
virtual double Integrate( double( f )( double x0, double x1, double x2 ) );
/*! @brief Integrates f(x,y,z) with the quadrature.
* @param double(f)( double my, double phi ) : density function that depends on a spherical coordinates.
* @returns double result: result of the quadrature rule */
virtual double IntegrateSpherical( double( f )( double my, double phi ) );
/*! @brief Integrates vector valued f(x,y,z) with the quadrature. Each dimension is integrated by itself.
* @param : double(f)( double x0, double x1, double x2 ) : density function that depends on a three spatial dimensions.
* @param : len : lenght of vector
* @returns double result: result of the quadrature rule (vector valued) */
std::vector<double> Integrate( std::vector<double>( f )( double x0, double x1, double x2 ), unsigned len );
virtual std::vector<double> Integrate( std::vector<double>( f )( double x0, double x1, double x2 ), unsigned len );
// Quadrature Hub
/*! @brief Creates a quadrature rule with a given name and a given order.
......@@ -60,8 +65,10 @@ class QuadratureBase
inline unsigned GetOrder() const { return _order; } /*! @returns unsigned _order: order of the quadrature */
inline unsigned GetNq() const { return _nq; } /*! @returns unsigned _nq: number of gridpoints of the quadrature */
inline VectorVector GetPoints() const { return _points; } /*! @returns VectorVector _points: coordinates of gridpoints of the quadrature */
virtual VectorVector GetPointsSphere() const; /*! @returns VectorVector _pointsSphere: "---- " in spherical coordinates (my, phi)*/
inline Vector GetWeights() const { return _weights; } /*! @returns Vector _weights: weights of gridpoints of the quadrature */
inline VectorVector GetPointsSphere() const {
return _pointsSphere;
} /*! @returns VectorVector _pointsSphere: "---- " in spherical coordinates (my, phi)*/
inline Vector GetWeights() const { return _weights; } /*! @returns Vector _weights: weights of gridpoints of the quadrature */
inline VectorVectorU GetConnectivity() const {
return _connectivity;
} /*! @returns VectorVectorU _connectivity: connectivity of gridpoints of the quadrature */
......@@ -84,10 +91,9 @@ class QuadratureBase
unsigned _order; /*! @brief order of the quadrature */
unsigned _nq; /*! @brief number of gridpoints of the quadrature */
VectorVector _points; /*! @brief gridpoints of the quadrature */
VectorVector _pointsSphere; /*! @brief (my,phi)gridpoints of the quadrature in spherical cordinates */
Vector _weights; /*! @brief weights of the gridpoints of the quadrature */
VectorVectorU _connectivity; /*! @brief connectivity of the gripoints of the quadrature */
VectorVector _pointsSphere; /*! @brief (my,phi)gridpoints of the quadrature in spherical cordinates */
};
#endif // QUADRATURE_H
......@@ -51,12 +51,13 @@ void QGaussLegendre1D::SetPointsAndWeights() {
// resize points and weights
_points.resize( _nq );
_weights.resize( _nq );
unsigned dim = 3;
unsigned dim = 1;
for( unsigned k = 0; k < _nq; ++k ) {
_points[k].resize( dim );
_points[k][0] = nodes1D[k];
_weights[k] = weights1D[k];
}
_pointsSphere = _points;
}
void QGaussLegendre1D::SetConnectivity() { // TODO
......@@ -139,4 +140,30 @@ double QGaussLegendre1D::Pythag( const double a, const double b ) {
: ( absb == 0.0 ? 0.0 : absb * std::sqrt( 1.0 + ( absa / absb ) * ( absa / absb ) ) ) );
}
bool QGaussLegendre1D::CheckOrder() { return true; }
bool QGaussLegendre1D::CheckOrder() {
return true; // All orders viable
}
double QGaussLegendre1D::Integrate( double( f )( double x0, double x1, double x2 ) ) { // Not Safe!
double result = 0.0;
double x = 0.0;
double y = 0.0;
double z = 0.0;
double w = 0.0;
for( unsigned i = 0; i < _nq; i++ ) {
x = _points[i][0];
w = _weights[i];
result += w * f( x, y, z );
}
return result;
}
double QGaussLegendre1D::IntegrateSpherical( double( f )( double my, double phi ) ) {
ErrorMessages::Error( "This method is not applicable for 1D quadratures.\n", CURRENT_FUNCTION );
return f( 0, 0 );
}
std::vector<double> QGaussLegendre1D::Integrate( std::vector<double>( f )( double x0, double x1, double x2 ), unsigned /* len */ ) {
ErrorMessages::Error( "This method is not applicable for 1D quadratures.\n", CURRENT_FUNCTION );
return { f( 0, 0, 0 ) };
}
#include "quadratures/qlookupquadrature.h"
#include "toolboxes/errormessages.h"
#include <fstream>
#include <math.h>
#include <sstream>
QLookupQuadrature::QLookupQuadrature( Config* settings ) : QuadratureBase( settings ) {}
......@@ -39,6 +40,7 @@ void QLookupQuadrature::SetPointsAndWeights() {
_weights.resize( _nq );
_points.resize( _nq );
_pointsSphere.resize( _nq );
double sumWeights = 0;
......@@ -72,12 +74,31 @@ void QLookupQuadrature::SetPointsAndWeights() {
if( ss.peek() == ',' ) ss.ignore();
}
}
// normalize the points
double norm = 0;
for( unsigned idx_point = 0; idx_point < _nq; idx_point++ ) {
norm = sqrt( _points[idx_point][0] * _points[idx_point][0] + _points[idx_point][1] * _points[idx_point][1] +
_points[idx_point][2] * _points[idx_point][2] );
_points[idx_point][0] /= norm;
_points[idx_point][1] /= norm;
_points[idx_point][2] /= norm;
}
// Correct the scaling of the weights
for( unsigned idx = 0; idx < _nq; idx++ ) {
_weights[idx] /= sumWeights;
_weights[idx] *= 4.0 * M_PI;
}
// Transform _points to _pointsSphere
// Transform _points to _pointsSphere ==>transform (x,y,z) into (my,phi)
for( unsigned idx = 0; idx < _nq; idx++ ) {
_pointsSphere[idx].resize( 2 ); // (my,phi)
_pointsSphere[idx][0] = _points[idx][2]; // my = z
_pointsSphere[idx][1] = atan2( _points[idx][1], _points[idx][0] ); // phi in [-pi,pi]
// adapt intervall s.t. phi in [0,2pi]
if( _pointsSphere[idx][1] < 0 ) {
_pointsSphere[idx][1] = 2 * M_PI + _pointsSphere[idx][1];
}
}
}
#include "quadratures/qmontecarlo.h"
#include <math.h>
QMonteCarlo::QMonteCarlo( Config* settings ) : QuadratureBase( settings ) {
SetName();
SetNq();
......@@ -15,38 +17,49 @@ QMonteCarlo::QMonteCarlo( unsigned quadOrder ) : QuadratureBase( quadOrder ) {
void QMonteCarlo::SetPointsAndWeights() {
// Nq random points on the sphere.
_points.resize( GetNq() );
_pointsSphere.resize( GetNq() );
_points.resize( GetNq(), Vector( 3, 0.0 ) );
_pointsSphere.resize( GetNq(), Vector( 2, 0.0 ) );
std::default_random_engine generator;
std::normal_distribution<double> distribution( 0.0, 1.0 );
std::uniform_real_distribution<double> distributionUni( 0.0, 1.0 );
// Tranlation in spherical coordinates
double my, phi, x, y, z, temp, norm;
double x, y, z, norm;
for( unsigned i = 0; i < GetNq(); i++ ) {
// https://mathworld.wolfram.com/SpherePointPicking.html describes how to generate random points on the sphere.
x = distribution( generator );
y = distribution( generator );
z = distribution( generator );
norm = sqrt( x * x + y * y + z * z );
Vector point( { x / norm, y / norm, z / norm } );
_points[i] = point;
x = distribution( generator );
y = distribution( generator );
z = distribution( generator );
norm = sqrt( x * x + y * y + z * z );
_points[i][0] = x / norm;
_points[i][1] = y / norm;
_points[i][2] = z / norm;
// Version for spherical coordinates (not 1:1 correspondence to karthesian in terms of distribution!)
temp = distributionUni( generator );
my = 2 * temp - 1;
phi = 2 * M_PI * temp;
Vector pointSphere( { my, phi } );
_pointsSphere[i] = pointSphere;
// temp = distributionUni( generator );
// my = 2 * temp - 1;
// phi = 2 * M_PI * temp;
//
// Vector pointSphere( { my, phi } );
// _pointsSphere[i] = pointSphere;
}
// Transform _points to _pointsSphere ==>transform (x,y,z) into (my,phi)
for( unsigned idx = 0; idx < _nq; idx++ ) {
_pointsSphere[idx].resize( 2 ); // (my,phi)
_pointsSphere[idx][0] = _points[idx][2]; // my = z
_pointsSphere[idx][1] = atan2( _points[idx][1], _points[idx][0] ); // phi in [-pi,pi]
// adapt intervall s.t. phi in [0,2pi]
if( _pointsSphere[idx][1] < 0 ) {
_pointsSphere[idx][1] = 2 * M_PI + _pointsSphere[idx][1];
}
}
// Equal weights
_weights = Vector( GetNq(), 4.0 * M_PI / GetNq() );
_weights = Vector( _nq, 4.0 * M_PI / ( (double)_nq ) );
}
VectorVector QMonteCarlo::GetPointsSphere() const { return _pointsSphere; }
void QMonteCarlo::SetConnectivity() { // TODO
// Not initialized for this quadrature.
VectorVectorU connectivity;
......
......@@ -51,19 +51,32 @@ QuadratureBase* QuadratureBase::Create( QUAD_NAME name, unsigned quadOrder ) {
double QuadratureBase::Integrate( double( f )( double x0, double x1, double x2 ) ) {
double result = 0.0;
double x = 0.0;
double y = 0.0;
double z = 0.0;
double w = 0.0;
for( unsigned i = 0; i < _nq; i++ ) {
double x = _points[i][0];
double y = _points[i][1];
double z = _points[i][2];
double w = _weights[i];
x = _points[i][0];
y = _points[i][1];
z = _points[i][2];
w = _weights[i];
result += w * f( x, y, z );
}
return result;
}
VectorVector QuadratureBase::GetPointsSphere() const {
ErrorMessages::Error( "Quadrature points in spherical coordinates are not supported\nfor this quadrature. Exiting", CURRENT_FUNCTION );
return _pointsSphere;
double QuadratureBase::IntegrateSpherical( double( f )( double my, double phi ) ) {
double result = 0.0;
double my = 0.0;
double phi = 0.0;
double w = 0.0;
for( unsigned i = 0; i < _nq; i++ ) {
my = _pointsSphere[i][0];
phi = _pointsSphere[i][1];
w = _weights[i];
result += w * f( my, phi );
}
return result;
}
std::vector<double> QuadratureBase::Integrate( std::vector<double>( f )( double x0, double x1, double x2 ), unsigned len ) {
......
This diff is collapsed.
......@@ -238,20 +238,20 @@ TEST_CASE( "test spherical harmonics basis ", "[spherical_harmonics]" ) {
}
}
double Omega_x( double my, double phi ) { return sqrt( 1 - my * my ) * sin( phi ); }
double Omega_y( double my, double phi ) { return sqrt( 1 - my * my ) * cos( phi ); }
double Omega_z( double my ) { return my; }
double Omega_xBase( double my, double phi ) { return sqrt( 1 - my * my ) * sin( phi ); }
double Omega_yBase( double my, double phi ) { return sqrt( 1 - my * my ) * cos( phi ); }
double Omega_zBase( double my ) { return my; }
double SphericalMonomial_0( double /* my */, double /* phi */ ) { return 1; }
double SphericalMonomial_1( double my, double /*phi*/ ) { return Omega_z( my ); } // omega_z
double SphericalMonomial_2( double my, double phi ) { return Omega_y( my, phi ); } // omega_y
double SphericalMonomial_3( double my, double phi ) { return Omega_x( my, phi ); } // omega_x
double SphericalMonomial_4( double my, double /*phi*/ ) { return Omega_z( my ) * Omega_z( my ); } // omega_z^2
double SphericalMonomial_5( double my, double phi ) { return Omega_y( my, phi ) * Omega_z( my ); } // omega_y*omega_z
double SphericalMonomial_6( double my, double phi ) { return Omega_y( my, phi ) * Omega_y( my, phi ); } // omega_y^2
double SphericalMonomial_7( double my, double phi ) { return Omega_x( my, phi ) * Omega_z( my ); } // omega_x*omega_z
double SphericalMonomial_8( double my, double phi ) { return Omega_x( my, phi ) * Omega_y( my, phi ); } // omega_x*omega_y
double SphericalMonomial_9( double my, double phi ) { return Omega_x( my, phi ) * Omega_x( my, phi ); } // omega_x^2
double SphericalMonomial_1( double my, double /*phi*/ ) { return Omega_zBase( my ); } // omega_z
double SphericalMonomial_2( double my, double phi ) { return Omega_yBase( my, phi ); } // omega_y
double SphericalMonomial_3( double my, double phi ) { return Omega_xBase( my, phi ); } // omega_x
double SphericalMonomial_4( double my, double /*phi*/ ) { return Omega_zBase( my ) * Omega_zBase( my ); } // omega_z^2
double SphericalMonomial_5( double my, double phi ) { return Omega_yBase( my, phi ) * Omega_zBase( my ); } // omega_y*omega_z
double SphericalMonomial_6( double my, double phi ) { return Omega_yBase( my, phi ) * Omega_yBase( my, phi ); } // omega_y^2
double SphericalMonomial_7( double my, double phi ) { return Omega_xBase( my, phi ) * Omega_zBase( my ); } // omega_x*omega_z
double SphericalMonomial_8( double my, double phi ) { return Omega_xBase( my, phi ) * Omega_yBase( my, phi ); } // omega_x*omega_y
double SphericalMonomial_9( double my, double phi ) { return Omega_xBase( my, phi ) * Omega_xBase( my, phi ); } // omega_x^2
TEST_CASE( "test spherical monomial basis", "[spherical_monomials]" ) {
unsigned maxMomentDegree = 2; //==> 6+3+1 basis functions
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment