qmontecarlo.cpp 1.67 KB
 Thomas committed Mar 30, 2020 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> QMonteCarlo::ComputePoints() { // Nq random points on the sphere. int nq = GetNq(); blaze::DynamicVector> points( nq ); std::default_random_engine generator; std::normal_distribution 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 point( {x / norm, y / norm, z / norm} ); points[i] = point; } return points; }; blaze::DynamicVector QMonteCarlo::ComputeWeights() { // Equal weights int nq = GetNq(); blaze::DynamicVector weights( nq, 4.0 * M_PI / nq ); return weights; }; blaze::DynamicVector> QMonteCarlo::ComputeConnectivity() { // Not initialized for this quadrature. blaze::DynamicVector> connectivity; return connectivity; };``````