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

code/include/quadratures/qgausslegendre1D.h

+/*!
+ * @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

code/include/quadratures/qgausslegendretensorized.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

code/include/quadratures/qmontecarlo.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

code/include/quadratures/qproduct.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

code/include/quadratures/quadraturebase.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

code/src/quadratures/qgausslegendre1D.cpp

@@ -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 ) };
+}

code/src/quadratures/qlookupquadrature.cpp

 #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];
+        }
+    }
 }

code/src/quadratures/qmontecarlo.cpp

 #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;

code/src/quadratures/quadraturebase.cpp

@@ -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 ) {

code/tests/test_sphericalBasis.cpp

@@ -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