added spherical gridpoints for MC and GL quadrature, prepared test cases....

added spherical gridpoints for MC and GL quadrature, prepared test cases. small changes to mnsolver. still does not work

code/include/quadratures/qgausslegendretensorized.h

@@ -18,6 +18,8 @@ class QGaussLegendreTensorized : public QuadratureBase
     inline void SetNq() override { _nq = pow( GetOrder(), 2 ); }
     void SetPointsAndWeights() override;
     void SetConnectivity() override;
+
+    inline VectorVector GetPointsSphere() const override { return _pointsSphere; }
 };
 
 #endif    // QGAUSSLEGENDRETENSORIZED_H

code/include/quadratures/qmontecarlo.h

@@ -13,6 +13,7 @@ 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/quadraturebase.h

@@ -43,6 +43,7 @@ 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 VectorVectorU GetConnectivity() const {
         return _connectivity;
@@ -61,12 +62,15 @@ class QuadratureBase
     virtual void SetPointsAndWeights() = 0;
 
     // Member variables
+    // TODO Config* _settings;           /*! @brief pointer to settings class that manages the solver */
     std::string _name;           /*! @brief name of the quadrature */
     unsigned _order;             /*! @brief order of the quadrature */
     unsigned _nq;                /*! @brief number of gridpoints of the quadrature */
     VectorVector _points;        /*! @brief gridpoints of the quadrature */
     Vector _weights;             /*! @brief weights of the gridpoints of the quadrature */
     VectorVectorU _connectivity; /*! @brief connectivity of the gripoints of the quadrature */
+
+    VectorVector _pointsSphere; /*! @brief gridpoints of the quadrature in spherical cordinates */
 };
 
 #endif    // QUADRATURE_H

code/src/main.cpp

@@ -15,78 +15,94 @@
 #include "quadratures/qmontecarlo.h"
 #include "solvers/sphericalharmonics.h"
 
+double testFunc( double my, double phi ) { return my * my + phi; }
+
+double testFunc2( double x, double y, double z ) { return x + y + z; }
+
 int main( int argc, char** argv ) {
     MPI_Init( &argc, &argv );
 
-    // QGaussLegendreTensorized quad( 4 );
-    // int maxMomentDegree = 2;
-    // SphericalHarmonics testBase( maxMomentDegree );
-    //
-    // double x, y, z, w;
-    // Vector moment = testBase.ComputeSphericalBasis( 0, 1, 0 );
-    //// 9 basis moments if degree = 2
-    //// unsigned nTotalMoments = moment.size();
+    QMonteCarlo quad( 1000 );
+    int maxMomentDegree = 2;
+
+    // test quadrature;
+    std::cout << " Integration test:" << quad.Integrate( testFunc2 ) << "\n";
+
+    SphericalHarmonics testBase( maxMomentDegree );
+
+    double my, phi, w, resTest = 0.0;
+    Vector moment = testBase.ComputeSphericalBasis( 0, 1, 0 );
+    // 9 basis moments if degree = 2
+    // unsigned nTotalMoments = moment.size();
     //
-    // Vector results( moment.size(), 0.0 );
+    Vector results( moment.size(), 0.0 );
+
+    int idx_basis = 0, idx_basisM = 0;
+
+    for( unsigned idx_quad = 0; idx_quad < quad.GetNq(); idx_quad++ ) {
+        my     = quad.GetPointsSphere()[idx_quad][0];
+        phi    = quad.GetPointsSphere()[idx_quad][1];
+        w      = quad.GetWeights()[idx_quad];
+        moment = testBase.ComputeSphericalBasis( my, phi );
+
+        double funcEval = testFunc( my, phi );
+        //  std::cout << "f(my= " << my << ", phi= " << phi / M_PI << "*pi) =  " << funcEval << " w " << w << "\n";
+        resTest += w * funcEval;
+
+        std::cout << "my, phi " << my << " , " << phi << " ::\t mom5 " << moment[5] << " mom7 \t" << moment[7] << "\n";
+        // results[5] += w * moment[5] * moment[5];
+        // results[7] += w * moment[7] * moment[7];
+        //
+        for( int idx_l = 0; idx_l <= maxMomentDegree; idx_l++ ) {
+            idx_basis = testBase.GlobalIdxBasis( idx_l, 0 );
+
+            results[idx_basis] += w * moment[idx_basis] * moment[idx_basis];
+
+            for( int idx_k = 1; idx_k <= idx_l; idx_k++ ) {
+                idx_basis  = testBase.GlobalIdxBasis( idx_l, idx_k );
+                idx_basisM = testBase.GlobalIdxBasis( idx_l, -idx_k );
+
+                results[idx_basis] += w * moment[idx_basis] * moment[idx_basis];
+                results[idx_basisM] += w * moment[idx_basisM] * moment[idx_basisM];
+            }
+        }
+    }
+    std::cout << " testInt " << resTest << "\n";
+    std::cout << " Squared Integration:\n " << results << "\n";
     //
-    // int idx_basis = 0, idx_basisM = 0;
+    // Normalized integration
+    // Vector resultsNormal( moment.size(), 0.0 );
     //
     // for( unsigned idx_quad = 0; idx_quad < quad.GetNq(); idx_quad++ ) {
-    //    x      = quad.GetPoints()[idx_quad][0];
-    //    y      = quad.GetPoints()[idx_quad][1];
-    //    z      = quad.GetPoints()[idx_quad][2];
-    //    w      = quad.GetWeights()[idx_quad];
-    //    moment = testBase.ComputeSphericalBasis( x, y, z );
+    //     x      = quad.GetPoints()[idx_quad][0];
+    //     y      = quad.GetPoints()[idx_quad][1];
+    //     z      = quad.GetPoints()[idx_quad][2];
+    //     w      = quad.GetWeights()[idx_quad];
+    //     moment = testBase.ComputeSphericalBasis( x, y, z );
     //
-    //    for( int idx_l = 0; idx_l <= 2; idx_l++ ) {
-    //        idx_basis = testBase.GlobalIdxBasis( idx_l, 0 );
+    //     for( unsigned idx_sys = 0; idx_sys < 9; idx_sys++ ) {
     //
-    //        results[idx_basis] += w * moment[idx_basis] * moment[idx_basis];
-    //
-    //        for( int idx_k = 1; idx_k <= idx_l; idx_k++ ) {
-    //            idx_basis  = testBase.GlobalIdxBasis( idx_l, idx_k );
-    //            idx_basisM = testBase.GlobalIdxBasis( idx_l, -idx_k );
-    //            // std::cout << "(" << idx_basis << " | " << idx_basisM << ")\n";
-    //            results[idx_basis] += w * moment[idx_basis] * moment[idx_basis];
-    //            results[idx_basisM] += w * moment[idx_basisM] * moment[idx_basisM];
-    //        }
-    //    }
-    //}
-    // std::cout << "Normalization integration:\n " << results << "\n";
+    //         resultsNormal[idx_sys] += w * moment[idx_sys] * moment[idx_sys] / results[idx_sys];
+    //         // std::cout << idx_quad << ": " << results[0] << "\n";
+    //     }
+    // }
+    // std::cout << "moment integration:\n " << resultsNormal << "\n";
+
+    // Vector results2( moment.size(), 0.0 );
     //
-    //// Normalized integration
-    //// Vector resultsNormal( moment.size(), 0.0 );
-    ////
-    //// for( unsigned idx_quad = 0; idx_quad < quad.GetNq(); idx_quad++ ) {
-    ////     x      = quad.GetPoints()[idx_quad][0];
-    ////     y      = quad.GetPoints()[idx_quad][1];
-    ////     z      = quad.GetPoints()[idx_quad][2];
-    ////     w      = quad.GetWeights()[idx_quad];
-    ////     moment = testBase.ComputeSphericalBasis( x, y, z );
-    ////
-    ////     for( unsigned idx_sys = 0; idx_sys < 9; idx_sys++ ) {
-    ////
-    ////         resultsNormal[idx_sys] += w * moment[idx_sys] * moment[idx_sys] / results[idx_sys];
-    ////         // std::cout << idx_quad << ": " << results[0] << "\n";
-    ////     }
-    //// }
-    //// std::cout << "moment integration:\n " << resultsNormal << "\n";
+    // for( unsigned idx_quad = 0; idx_quad < quad.GetNq(); idx_quad++ ) {
+    //     x      = quad.GetPoints()[idx_quad][0];
+    //     y      = quad.GetPoints()[idx_quad][1];
+    //     z      = quad.GetPoints()[idx_quad][2];
+    //     w      = quad.GetWeights()[idx_quad];
+    //     moment = testBase.ComputeSphericalBasis( x, y, z );
     //
-    //// Vector results2( moment.size(), 0.0 );
-    ////
-    //// for( unsigned idx_quad = 0; idx_quad < quad.GetNq(); idx_quad++ ) {
-    ////     x      = quad.GetPoints()[idx_quad][0];
-    ////     y      = quad.GetPoints()[idx_quad][1];
-    ////     z      = quad.GetPoints()[idx_quad][2];
-    ////     w      = quad.GetWeights()[idx_quad];
-    ////     moment = testBase.ComputeSphericalBasis( x, y, z );
-    ////
-    ////     for( unsigned idx_sys = 1; idx_sys < 9; idx_sys++ ) {
-    ////         results2[idx_sys] += w * moment[idx_sys - 1] * moment[idx_sys];
-    ////         // std::cout << idx_quad << ": " << results[0] << "\n";
-    ////     }
-    //// }
-    //// std::cout << "moment integration:\n " << results2 << "\n";
+    //     for( unsigned idx_sys = 1; idx_sys < 9; idx_sys++ ) {
+    //         results2[idx_sys] += w * moment[idx_sys - 1] * moment[idx_sys];
+    //         // std::cout << idx_quad << ": " << results[0] << "\n";
+    //     }
+    // }
+    // std::cout << "moment integration:\n " << results2 << "\n";
     //
     // std::ofstream myfile;
     // myfile.open( "assLegendre.csv" );
@@ -106,20 +122,20 @@ int main( int argc, char** argv ) {
     // myfile.close();
     //
 
-    std::string filename = ParseArguments( argc, argv );
-
-    // CD  Load Settings from File
-    Config* config = new Config( filename );
-
-    // Print input file and run info to file
-    PrintLogHeader( filename );
-
-    // Build solver
-    Solver* solver = Solver::Create( config );
-
-    // Run solver and export
-    solver->Solve();
-    solver->Save();
+    //  std::string filename = ParseArguments( argc, argv );
+    //
+    //  // CD  Load Settings from File
+    //  Config* config = new Config( filename );
+    //
+    //  // Print input file and run info to file
+    //  PrintLogHeader( filename );
+    //
+    //  // Build solver
+    //  Solver* solver = Solver::Create( config );
+    //
+    //  // Run solver and export
+    //  solver->Solve();
+    //  solver->Save();
 
     MPI_Finalize();
     return EXIT_SUCCESS;

code/src/quadratures/qgausslegendretensorized.cpp

@@ -53,9 +53,14 @@ void QGaussLegendreTensorized::SetPointsAndWeights() {
 
     // resize points and weights
     _points.resize( _nq );
+    _pointsSphere.resize( _nq );
     for( auto& p : _points ) {
         p.resize( 3 );
     }
+    for( auto& p : _pointsSphere ) {
+        p.resize( 2 );
+    }
+
     _weights.resize( _nq );
 
     // transform tensorized (x,y,z)-grid to spherical grid points
@@ -64,7 +69,11 @@ void QGaussLegendreTensorized::SetPointsAndWeights() {
             _points[j * ( 2 * _order ) + i][0] = sqrt( 1 - nodes1D[j] * nodes1D[j] ) * std::cos( phi[i] );
             _points[j * ( 2 * _order ) + i][1] = sqrt( 1 - nodes1D[j] * nodes1D[j] ) * std::sin( phi[i] );
             _points[j * ( 2 * _order ) + i][2] = nodes1D[j];
-            _weights[j * ( 2 * _order ) + i]   = 2.0 * M_PI / _order * weights1D[j];
+
+            _pointsSphere[j * ( 2 * _order ) + i][0] = nodes1D[j];
+            _pointsSphere[j * ( 2 * _order ) + i][0] = phi[i];
+
+            _weights[j * ( 2 * _order ) + i] = 2.0 * M_PI / _order * weights1D[j];
         }
     }
 }

code/src/quadratures/qmontecarlo.cpp

@@ -10,24 +10,37 @@ QMonteCarlo::QMonteCarlo( unsigned order ) : QuadratureBase( order ) {
 void QMonteCarlo::SetPointsAndWeights() {
     // Nq random points on the sphere.
     _points.resize( GetNq() );
-
+    _pointsSphere.resize( GetNq() );
     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;
 
     for( unsigned i = 0; i < GetNq(); i++ ) {
         // https://mathworld.wolfram.com/SpherePointPicking.html describes how to generate random points on the sphere.
-        double x    = distribution( generator );
-        double y    = distribution( generator );
-        double z    = distribution( generator );
-        double norm = sqrt( x * x + y * y + z * z );
+        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;
+
+        // 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;
     }
 
     // Equal weights
     _weights = Vector( GetNq(), 4.0 * M_PI / GetNq() );
 }
 
+VectorVector QMonteCarlo::GetPointsSphere() const { return _pointsSphere; }
+
 void QMonteCarlo::SetConnectivity() {    // TODO
     // Not initialized for this quadrature.
     VectorVectorU connectivity;

code/src/quadratures/quadraturebase.cpp

@@ -31,6 +31,11 @@ double QuadratureBase::Integrate( double( f )( double x0, double x1, double x2 )
     return result;
 }
 
+VectorVector QuadratureBase::GetPointsSphere() const {
+    ErrorMessages::Error( "Quadrature points in spherical coordinates are not supported\nfor this quadrature. Exiting", CURRENT_FUNCTION );
+    return _pointsSphere;
+}
+
 std::vector<double> QuadratureBase::Integrate( std::vector<double>( f )( double x0, double x1, double x2 ), unsigned len ) {
     std::vector<double> result( len, 0.0 );
     std::vector<double> funcEval( len, 0.0 );

code/src/solvers/mnsolver.cpp

@@ -50,36 +50,33 @@ int MNSolver::GlobalIndex( int l, int k ) const {
 }
 
 void MNSolver::ComputeMoments() {
-    double x, y, z, w;
+    double my, phi, w;
 
     for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
-        x                  = _quadrature->GetPoints()[idx_quad][0];
-        y                  = _quadrature->GetPoints()[idx_quad][1];
-        z                  = _quadrature->GetPoints()[idx_quad][2];
-        _moments[idx_quad] = _basis.ComputeSphericalBasis( x, y, z );
+        phi = _quadrature->GetPointsSphere()[idx_quad][0];
+        my  = _quadrature->GetPointsSphere()[idx_quad][1];
+
+        _moments[idx_quad] = _basis.ComputeSphericalBasis( my, phi );
     }
 
     Vector normalizationFactor( _moments[0].size(), 0.0 );
 
     // Normalize Moments
 
-    for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
-        x = _quadrature->GetPoints()[idx_quad][0];
-        y = _quadrature->GetPoints()[idx_quad][1];
-        z = _quadrature->GetPoints()[idx_quad][2];
-        w = _quadrature->GetWeights()[idx_quad];
-
-        for( unsigned idx_sys = 0; idx_sys < _nTotalEntries; idx_sys++ ) {
-            normalizationFactor[idx_sys] += w * _moments[idx_quad][idx_sys] * _moments[idx_quad][idx_sys];
-        }
-    }
-    std::cout << "norm" << normalizationFactor << "\n";
-
-    for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
-        for( unsigned idx_sys = 0; idx_sys < _nTotalEntries; idx_sys++ ) {
-            _moments[idx_quad][idx_sys] /= sqrt( normalizationFactor[idx_sys] );
-        }
-    }
+    // for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
+    //    w = _quadrature->GetWeights()[idx_quad];
+    //
+    //    for( unsigned idx_sys = 0; idx_sys < _nTotalEntries; idx_sys++ ) {
+    //        normalizationFactor[idx_sys] += w * _moments[idx_quad][idx_sys] * _moments[idx_quad][idx_sys];
+    //    }
+    //}
+    // std::cout << "norm" << normalizationFactor << "\n";
+    //
+    // for( unsigned idx_quad = 0; idx_quad < _nq; idx_quad++ ) {
+    //    for( unsigned idx_sys = 0; idx_sys < _nTotalEntries; idx_sys++ ) {
+    //        _moments[idx_quad][idx_sys] /= sqrt( normalizationFactor[idx_sys] );
+    //    }
+    //}
 }
 
 Vector MNSolver::ConstructFlux( unsigned idx_cell ) {

code/tests/test_quadrature.cpp

@@ -92,6 +92,8 @@ double f( double x, double y, double z ) {
     return x * x + y * y + z * z;    // == 1
 }
 
+double g( double x, double y, double z ) { return x + y + z; }
+
 TEST_CASE( "Integrate a constant function.", "[quadrature]" ) {
     bool lowAccuracyTesting = false;
     for( auto quadraturename : quadraturenames ) {
@@ -111,6 +113,17 @@ TEST_CASE( "Integrate a constant function.", "[quadrature]" ) {
                 printf( "Computed result %.15f", Q->Integrate( f ) );
             }
             REQUIRE( approxequal( Q->Integrate( f ), 4.0 * M_PI, lowAccuracyTesting ) );
+
+            // non constant function
+            // if( !approxequal( Q->Integrate( g ), -0.0997858, lowAccuracyTesting ) ) {
+            //     printf( "Quadrature %d at order %d :  Error : %.15f (low accuracy testing was set to %d)\n",
+            //             quadraturename,
+            //             quadratureorder,
+            //             std::abs( Q->Integrate( g ) + 0.0997858 ),
+            //             lowAccuracyTesting );
+            //     printf( "Computed result %.15f", Q->Integrate( g ) );
+            // }
+            // REQUIRE( approxequal( Q->Integrate( g ), -0.0997858, lowAccuracyTesting ) );
         }
     }
 }