quadraturebase.cpp 4.16 KB
Newer Older
1
#include "quadratures/quadraturebase.h"
2
#include "common/config.h"
Jonas Kusch's avatar
Jonas Kusch committed
3
#include "quadratures/qgausslegendre1D.h"
4
5
6
#include "quadratures/qgausslegendretensorized.h"
#include "quadratures/qldfesa.h"
#include "quadratures/qlebedev.h"
7
8
#include "quadratures/qlevelsymmetric.h"
#include "quadratures/qmontecarlo.h"
9
#include "quadratures/qproduct.h"
steffen.schotthoefer's avatar
steffen.schotthoefer committed
10
#include "toolboxes/errormessages.h"
steffen.schotthoefer's avatar
steffen.schotthoefer committed
11

12
13
14
15
QuadratureBase::QuadratureBase( Config* settings ) {
    _settings = settings;
    _order    = settings->GetQuadOrder();
}
16

Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
17
18
QuadratureBase::QuadratureBase( unsigned order ) : _order( order ) { _settings = nullptr; }

19
20
QuadratureBase* QuadratureBase::CreateQuadrature( Config* settings ) {
    QUAD_NAME name = settings->GetQuadName();
21

22
    switch( name ) {
23
24
25
26
27
28
        case QUAD_MonteCarlo: return new QMonteCarlo( settings );
        case QUAD_GaussLegendreTensorized: return new QGaussLegendreTensorized( settings );
        case QUAD_GaussLegendre1D: return new QGaussLegendre1D( settings );
        case QUAD_LevelSymmetric: return new QLevelSymmetric( settings );
        case QUAD_LDFESA: return new QLDFESA( settings );
        case QUAD_Lebedev: return new QLebedev( settings );
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
29
30
31
32
        case QUAD_Product: return new ProductQuadrature( settings );
        default: ErrorMessages::Error( "Creator for the chose quadrature does not yet exist. This is is the fault of the coder!", CURRENT_FUNCTION );
    }
}
33

Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
34
QuadratureBase* QuadratureBase::CreateQuadrature( QUAD_NAME name, unsigned quadOrder ) {
35

36
    switch( name ) {
Steffen Schotthöfer's avatar
Steffen Schotthöfer committed
37
38
39
40
41
42
43
44
45
46
        case QUAD_MonteCarlo: return new QMonteCarlo( quadOrder );
        case QUAD_GaussLegendreTensorized:
            ErrorMessages::Error( "This quadrature must be initialized with a settings constructor!", CURRENT_FUNCTION );
            break;
        case QUAD_GaussLegendre1D: return new QGaussLegendre1D( quadOrder );
        case QUAD_LevelSymmetric: return new QLevelSymmetric( quadOrder );
        case QUAD_LDFESA: return new QLDFESA( quadOrder );
        case QUAD_Lebedev: return new QLebedev( quadOrder );
        case QUAD_Product: return new ProductQuadrature( quadOrder );
        default: ErrorMessages::Error( "Creator for the chose quadrature does not yet exist. This is is the fault of the coder!", CURRENT_FUNCTION );
47
48
    }
}
49

50
double QuadratureBase::Integrate( double( f )( double x0, double x1, double x2 ) ) {
51
    double result = 0.0;
52
53
54
55
56
57
58
59
60
    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];
        result += w * f( x, y, z );
    }
    return result;
}
61

62
63
64
65
66
VectorVector QuadratureBase::GetPointsSphere() const {
    ErrorMessages::Error( "Quadrature points in spherical coordinates are not supported\nfor this quadrature. Exiting", CURRENT_FUNCTION );
    return _pointsSphere;
}

67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
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 );

    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];
        funcEval = f( x, y, z );
        for( unsigned idx_len = 0; idx_len < len; idx_len++ ) {
            result[idx_len] += w * funcEval[idx_len];
        }
    }
    return result;
}

84
double QuadratureBase::SumUpWeights() { return sum( _weights ); }
85

86
void QuadratureBase::PrintWeights() {
87
    auto log = spdlog::get( "event" );
88
    for( unsigned i = 0; i < _nq; i++ ) {
89
        double w = _weights[i];
90
        log->info( w );
91
92
    }
}
93

94
void QuadratureBase::PrintPoints() {
95
    auto log = spdlog::get( "event" );
96
    for( unsigned i = 0; i < _nq; i++ ) {
97
98
99
        double x = _points[i][0];
        double y = _points[i][1];
        double z = _points[i][2];
100
        log->info( "{0}, {1}, {2}", x, y, z );
101
102
    }
}
103
void QuadratureBase::PrintPointsAndWeights() {
104
    auto log = spdlog::get( "event" );
105
    for( unsigned i = 0; i < _nq; i++ ) {
106
107
108
109
        double x = _points[i][0];
        double y = _points[i][1];
        double z = _points[i][2];
        double w = _weights[i];
110
        log->info( "{0}, {1}, {2}, {3}", x, y, z, w );
111
112
    }
}