added monomial basis on unit sphere, added unit tests for this basis

code/include/common/config.h

@@ -83,6 +83,9 @@ class Config
     // Scattering Kernel
     KERNEL_NAME _kernelName; /*!< @brief Scattering Kernel Name*/
 
+    // Spherical Basis
+    SPHERICAL_BASIS_NAME _sphericalBasisName; /*!< @brief: Name of the basis on the unit sphere */
+
     // Optimizer
     OPTIMIZER_NAME _entropyOptimizerName; /*!< @brief Choice of optimizer */
     double _optimizerEpsilon;             /*!< @brief termination criterion epsilon for Newton Optmizer */
@@ -273,6 +276,8 @@ class Config
     // Scattering Kernel
     KERNEL_NAME inline GetKernelName() const { return _kernelName; }
 
+    // Basis name
+    SPHERICAL_BASIS_NAME inline GetSphericalBasisName() const { return _sphericalBasisName; }
     // Output Structure
     std::vector<VOLUME_OUTPUT> inline GetVolumeOutput() { return _volumeOutput; }
     unsigned short inline GetNVolumeOutput() { return _nVolumeOutput; }

code/include/common/globalconstants.h

@@ -141,4 +141,10 @@ enum SCALAR_OUTPUT { ITER, MASS, RMS_FLUX, VTK_OUTPUT, CSV_OUTPUT };
 inline std::map<std::string, SCALAR_OUTPUT> ScalarOutput_Map{
     { "ITER", ITER }, { "MASS", MASS }, { "RMS_FLUX", RMS_FLUX }, { "VTK_OUTPUT", VTK_OUTPUT }, { "CSV_OUTPUT", CSV_OUTPUT } };
 
+// Spherical Basis Name
+enum SPHERICAL_BASIS_NAME { SPHERICAL_HARMONICS, SPHERICAL_MONIMIALS };
+
+inline std::map<std::string, SPHERICAL_BASIS_NAME> SphericalBasis_Map{ { "SPHERICAL_HARMONICS", SPHERICAL_HARMONICS },
+                                                                       { "SPHERICAL_MONIMIALS", SPHERICAL_MONIMIALS } };
+
 #endif    // GLOBAL_CONSTANTS_H

code/include/toolboxes/sphericalharmonics.h

@@ -69,7 +69,7 @@ class SphericalHarmonics : public SphericalBasisBase
      */
     std::vector<double> _assLegendreP;
 
-    /*! @brief sperical harmonic basis functions of
+    /*! @brief: spherical harmonic basis functions of
      *         degree  0 <= l <= L and order -l <= k <= l
      *         length : N = L² +2L
      */

code/src/common/config.cpp

@@ -283,6 +283,10 @@ void Config::SetConfigOptions() {
     /*! @brief Scattering kernel \n DESCRIPTION: Describes used scattering kernel \n DEFAULT KERNEL_Isotropic \ingroup Config */
     AddEnumOption( "KERNEL", _kernelName, Kernel_Map, KERNEL_Isotropic );
 
+    /*! @brief Spherical Basis \n DESCRIPTION: Describes the chosen set of basis functions for on the unit sphere (e.g. for Moment methods) \n DEFAULT
+     * SPHERICAL_HARMONICS \ingroup Config */
+    AddEnumOption( "SPHERICAL_BASIS", _sphericalBasisName, SphericalBasis_Map, SPHERICAL_HARMONICS );
+
     // Output related options
     /*! @brief Volume output \n DESCRIPTION: Describes output groups to write to vtk \ingroup Config */
     AddEnumListOption( "VOLUME_OUTPUT", _nVolumeOutput, _volumeOutput, VolOutput_Map );

code/src/quadratures/quadraturebase.cpp

@@ -29,6 +29,7 @@ QuadratureBase* QuadratureBase::Create( Config* settings ) {
         case QUAD_Product: return new QProduct( settings );
         default: ErrorMessages::Error( "Creator for the chose quadrature does not yet exist. This is is the fault of the coder!", CURRENT_FUNCTION );
     }
+    return nullptr;
 }
 
 QuadratureBase* QuadratureBase::Create( QUAD_NAME name, unsigned quadOrder ) {
@@ -45,6 +46,7 @@ QuadratureBase* QuadratureBase::Create( QUAD_NAME name, unsigned quadOrder ) {
         case QUAD_Product: return new QProduct( quadOrder );
         default: ErrorMessages::Error( "Creator for the chose quadrature does not yet exist. This is is the fault of the coder!", CURRENT_FUNCTION );
     }
+    return nullptr;
 }
 
 double QuadratureBase::Integrate( double( f )( double x0, double x1, double x2 ) ) {

code/src/toolboxes/sphericalbasisbase.cpp

+/*!
+ * @file sphericalbasisbase.cpp
+ * @brief Base Class to handle basis classes on the unit sphere
+ * @author S. Schotthöfer
+ *
+ */
+
+#include "toolboxes/sphericalbasisbase.h"
+#include "common/config.h"
+#include "toolboxes/errormessages.h"
+#include "toolboxes/sphericalharmonics.h"
+#include "toolboxes/sphericalmonomials.h"
+
+SphericalBasisBase* SphericalBasisBase::Create( Config* settings ) {
+    SPHERICAL_BASIS_NAME name = settings->GetSphericalBasisName();
+    unsigned maxMomentDegree  = settings->GetMaxMomentDegree();
+
+    switch( name ) {
+        case SPHERICAL_HARMONICS: return new SphericalHarmonics( maxMomentDegree );
+        case SPHERICAL_MONIMIALS: return new SphericalMonomials( maxMomentDegree );
+
+        default: ErrorMessages::Error( "Creator for the chosen basis does not yet exist. This is is the fault of the coder!", CURRENT_FUNCTION );
+    }
+    return nullptr;
+}

code/src/toolboxes/sphericalmonomials.cpp

+/*!
+ * @file sphericalmonomials.cpp
+ * @brief Class for efficient computation of a spherical monomial basis
+ * @author S. Schotthöfer
+ */
+
+#include "toolboxes/sphericalmonomials.h"
+
+SphericalMonomials::SphericalMonomials( unsigned L_degree ) {
+    _LMaxDegree = L_degree;
+    _spatialDim = 3;    // Spatial dimension 2 is hardcoded
+
+    unsigned basisSize = ComputeBasisSize( L_degree );
+    _YBasis            = Vector( basisSize );
+}
+
+Vector SphericalMonomials::ComputeSphericalBasis( double my, double phi ) {
+    unsigned idx_vector = 0;
+    // go over all degrees of polynomials
+    for( unsigned idx_degree = 0; idx_degree <= _LMaxDegree; idx_degree++ ) {
+        // elem = Omega_x^a+ Omega_y^b +Omega_z^c : a+b+c = idx_degree
+        double omega_X = Omega_x( my, phi );
+        double omega_Y = Omega_y( my, phi );
+        double omega_Z = Omega_z( my );
+        double tmp     = .0;
+        for( unsigned a = 0; a <= idx_degree; a++ ) {
+            for( unsigned b = 0; b <= idx_degree - a; b++ ) {
+                unsigned c          = idx_degree - a - b;    // c uniquely defined
+                tmp                 = Power( omega_X, a ) * Power( omega_Y, b ) * Power( omega_Z, c );
+                _YBasis[idx_vector] = Power( omega_X, a ) * Power( omega_Y, b ) * Power( omega_Z, c );
+                idx_vector++;
+            }
+        }
+    }
+    return _YBasis;
+}
+
+Vector SphericalMonomials::ComputeSphericalBasis( double x, double y, double z ) {
+    // transform (x,y,z) into (my,phi)
+    double my  = z;
+    double phi = 0.0;
+
+    if( y >= 0 )
+        phi = acos( x );
+    else
+        phi = 2 * M_PI - acos( x );
+
+    return ComputeSphericalBasis( my, phi );
+}
+
+/*! @brief: Computes the length of the basis vector for a given max degree and
+ *          spatial dimension dim. len of a single oder: (degree + _spatialDim -1) over (degree)
+ *  @return: lenght of whole basis */
+unsigned SphericalMonomials::ComputeDimensionSize( unsigned degree ) {
+    return Factorial( degree + _spatialDim - 1 ) / ( Factorial( degree ) * Factorial( _spatialDim - 1 ) );
+}
+
+unsigned SphericalMonomials::ComputeBasisSize( unsigned degree ) {
+    unsigned basisLen = 0;
+    for( unsigned idx_degree = 0; idx_degree <= degree; idx_degree++ ) {
+        basisLen += ComputeDimensionSize( idx_degree );
+    }
+    return basisLen;
+}
+
+unsigned SphericalMonomials::Factorial( unsigned n ) { return ( n == 1 || n == 0 ) ? 1 : Factorial( n - 1 ) * n; }
+
+double SphericalMonomials::Omega_x( double my, double phi ) { return sqrt( 1 - my * my ) * sin( phi ); }
+double SphericalMonomials::Omega_y( double my, double phi ) { return sqrt( 1 - my * my ) * cos( phi ); }
+double SphericalMonomials::Omega_z( double my ) { return my; }
+
+double SphericalMonomials::Power( double basis, unsigned exponent ) {
+    if( exponent == 0 ) return 1.0;
+    double result = basis;
+    for( unsigned i = 1; i < exponent; i++ ) {
+        result = result * basis;
+    }
+    return result;
+}

code/tests/test_harmonics.cpp → code/tests/test_sphericalBasis.cpp

@@ -3,6 +3,7 @@
 #include "common/config.h"
 #include "quadratures/qgausslegendretensorized.h"
 #include "toolboxes/sphericalharmonics.h"
+#include "toolboxes/sphericalmonomials.h"
 
 #include <fstream>
 #include <iostream>
@@ -28,7 +29,7 @@ double P2_0( double my ) { return sqrt( 5 / ( 8 * M_PI ) ) * ( 3 * my * my - 1 )
 double P2_1( double my ) { return -1 * sqrt( 15 / ( 4 * M_PI ) ) * my * sqrt( 1 - my * my ); }
 double P2_2( double my ) { return sqrt( 15 / ( 16 * M_PI ) ) * ( 1 - my * my ); }
 
-TEST_CASE( "test the spherical harmonics basis computation", "[spherical_harmonics]" ) {
+TEST_CASE( "test  spherical harmonics basis ", "[spherical_harmonics]" ) {
 
     std::string filename = std::string( TESTS_PATH ) + "input/unit_tests/solvers/unit_harmonics.cfg";
 
@@ -221,3 +222,57 @@ TEST_CASE( "test the spherical harmonics basis computation", "[spherical_harmoni
         REQUIRE( errorWithinBounds );
     }
 }
+
+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 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
+
+TEST_CASE( "test spherical monomial basis", "[spherical_monomials]" ) {
+    unsigned maxMomentDegree = 2;    //==> 6+3+1 basis functions
+    SphericalMonomials testBase( maxMomentDegree );
+
+    SECTION( "Test against analytical solution" ) {
+        Vector moment;
+        std::vector<bool> validMoment( 10, true );
+        for( double my = -1.0; my < 1.0; my += 0.1 ) {
+
+            for( double phi = 0.0; phi < 2 * M_PI; phi += 0.1 ) {
+                moment = testBase.ComputeSphericalBasis( my, phi );
+
+                std::cout << "(algo|real): (" << moment[0] << " | " << SphericalMonomial_0( my, phi ) << ")" << std::endl;
+                std::cout << "(algo|real): (" << moment[1] << " | " << SphericalMonomial_1( my, phi ) << ")" << std::endl;
+                std::cout << "(algo|real): (" << moment[2] << " | " << SphericalMonomial_2( my, phi ) << ")" << std::endl;
+                std::cout << "(algo|real): (" << moment[3] << " | " << SphericalMonomial_3( my, phi ) << ")" << std::endl;
+                std::cout << "(algo|real): (" << moment[4] << " | " << SphericalMonomial_4( my, phi ) << ")" << std::endl;
+                std::cout << "(algo|real): (" << moment[5] << " | " << SphericalMonomial_5( my, phi ) << ")" << std::endl;
+                std::cout << "(algo|real): (" << moment[6] << " | " << SphericalMonomial_6( my, phi ) << ")" << std::endl;
+                std::cout << "(algo|real): (" << moment[7] << " | " << SphericalMonomial_7( my, phi ) << ")" << std::endl;
+                std::cout << "(algo|real): (" << moment[8] << " | " << SphericalMonomial_8( my, phi ) << ")" << std::endl;
+                std::cout << "(algo|real): (" << moment[9] << " | " << SphericalMonomial_9( my, phi ) << ")" << std::endl;
+
+                if( std::fabs( moment[0] - SphericalMonomial_0( my, phi ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validMoment[0] = false;
+                if( std::fabs( moment[1] - SphericalMonomial_1( my, phi ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validMoment[1] = false;
+                if( std::fabs( moment[2] - SphericalMonomial_2( my, phi ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validMoment[2] = false;
+                if( std::fabs( moment[3] - SphericalMonomial_3( my, phi ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validMoment[3] = false;
+                if( std::fabs( moment[4] - SphericalMonomial_4( my, phi ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validMoment[4] = false;
+                if( std::fabs( moment[5] - SphericalMonomial_5( my, phi ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validMoment[5] = false;
+                if( std::fabs( moment[6] - SphericalMonomial_6( my, phi ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validMoment[6] = false;
+                if( std::fabs( moment[7] - SphericalMonomial_7( my, phi ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validMoment[7] = false;
+                if( std::fabs( moment[8] - SphericalMonomial_8( my, phi ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validMoment[8] = false;
+                if( std::fabs( moment[9] - SphericalMonomial_9( my, phi ) ) > 1e2 * std::numeric_limits<double>::epsilon() ) validMoment[8] = false;
+            }
+        }
+        REQUIRE( std::all_of( validMoment.begin(), validMoment.end(), []( bool v ) { return v; } ) );
+    }
+}