qmontecarlo.cpp 1.67 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
#include "qmontecarlo.h"

QMonteCarlo::QMonteCarlo( int order ) : Quadrature( order ) {
    SetName( ComputeName() );
    SetNq( ComputeNq() );
    SetPoints( ComputePoints() );
    SetWeights( ComputeWeights() );
    SetConnectivity( ComputeConnectivity() );
};

std::string QMonteCarlo::ComputeName() {
    // Human readable description.
    return "Monte Carlo Quadrature";
};

int QMonteCarlo::ComputeNq() {
    // For Monte Carlo Quadrature, nq = order.
    return GetOrder();
};

blaze::DynamicVector<blaze::DynamicVector<double>> QMonteCarlo::ComputePoints() {
    // Nq random points on the sphere.
    int nq = GetNq();
    blaze::DynamicVector<blaze::DynamicVector<double>> points( nq );

    std::default_random_engine generator;
    std::normal_distribution<double> distribution( 0.0, 1.0 );

    for( int i = 0; i < nq; 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 );
        blaze::DynamicVector<double> point( {x / norm, y / norm, z / norm} );
        points[i] = point;
    }
    return points;
};

blaze::DynamicVector<double> QMonteCarlo::ComputeWeights() {
    // Equal weights
    int nq = GetNq();
    blaze::DynamicVector<double> weights( nq, 4.0 * M_PI / nq );
    return weights;
};

blaze::DynamicVector<blaze::DynamicVector<int>> QMonteCarlo::ComputeConnectivity() {
    // Not initialized for this quadrature.
    blaze::DynamicVector<blaze::DynamicVector<int>> connectivity;
    return connectivity;
};