Commit 8f055f77 authored by Thomas's avatar Thomas
Browse files

Initial dummy quadrature implementation.

parent 048e360f
#ifndef QUADRATURE_H
#define QUADRATURE_H
#include <blaze/Math.h>
#include <math.h>
#include <random>
#include <string>
class Quadrature
{
public:
Quadrature( int order ) {
SetOrder( order );
SetName( ComputeName() );
SetNq( ComputeNq() );
SetPoints( ComputePoints() );
SetWeights( ComputeWeights() );
SetConnectivity( ComputeConnectivity() );
};
virtual ~Quadrature() = 0;
// Virtual methods that every quadrature has to implement
virtual std::string ComputeName();
virtual int ComputeNq();
virtual blaze::DynamicVector<blaze::DynamicVector<double>> ComputePoints();
virtual blaze::DynamicVector<double> ComputeWeights();
virtual blaze::DynamicVector<blaze::DynamicVector<int>> ComputeConnectivity();
// Setter
void SetName( std::string name ) { _name = name; };
void SetOrder( int order ) { _order = order; };
void SetNq( int nq ) { _nq = nq; };
void SetPoints( blaze::DynamicVector<blaze::DynamicVector<double>> points ) { _points = points; };
void SetWeights( blaze::DynamicVector<double> weights ) { _weights = weights; };
void SetConnectivity( blaze::DynamicVector<blaze::DynamicVector<int>> connectivity ) { _connectivity = connectivity; };
// Getter
std::string GetName() { return _name; };
int GetOrder() { return _order; };
int GetNq() { return _nq; };
blaze::DynamicVector<blaze::DynamicVector<double>> GetPoints() { return _points; };
blaze::DynamicVector<double> GetWeights() { return _weights; };
blaze::DynamicVector<blaze::DynamicVector<int>> GetConnectivity() { return _connectivity; };
protected:
std::string _name;
int _order;
int _nq;
blaze::DynamicVector<blaze::DynamicVector<double>> _points;
blaze::DynamicVector<double> _weights;
blaze::DynamicVector<blaze::DynamicVector<int>> _connectivity;
};
#endif // QUADRATURE_H
#include "quadrature.h"
#include <iostream>
int main( int argc, char** argv ) {
std::cout << "Hello world!" << std::endl;
Quadrature* Q = GetQuadrature( "montecarlo", 10 );
Q->PrintWeights();
return EXIT_SUCCESS;
}
#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;
};
#ifndef QMONTECARLO_H
#define QMONTECARLO_H
#include "quadrature.h"
class QMonteCarlo:public Quadrature
{
public:
QMonteCarlo( int order );
virtual ~QMonteCarlo() {}
std::string ComputeName();
int ComputeNq();
blaze::DynamicVector<blaze::DynamicVector<double>> ComputePoints();
blaze::DynamicVector<double> ComputeWeights();
blaze::DynamicVector<blaze::DynamicVector<int>> ComputeConnectivity();
};
#endif // QMONTECARLO_H
#include "quadrature.h"
#include "qmontecarlo.h"
Quadrature::Quadrature( int order ) : _order( order ) {}
Quadrature* GetQuadrature( std::string name, int order ) {
if( name == "montecarlo" ) {
return new QMonteCarlo( order );
}
return new QMonteCarlo( order );
}
#ifndef QUADRATURE_H
#define QUADRATURE_H
#include <blaze/Math.h>
#include <iostream>
#include <string>
#include <vector>
using namespace std;
class Quadrature
{
public:
Quadrature( int order );
virtual ~Quadrature(){};
virtual std::string ComputeName() = 0;
virtual int ComputeNq() = 0;
virtual blaze::DynamicVector<blaze::DynamicVector<double>> ComputePoints() = 0;
virtual blaze::DynamicVector<double> ComputeWeights() = 0;
virtual blaze::DynamicVector<blaze::DynamicVector<int>> ComputeConnectivity() = 0;
void PrintWeights() {
for( int i = 0; i < _nq; i++ ) {
std::cout << _weights[i] << std::endl;
}
}
// Setter
void SetName( std::string name ) { _name = name; };
void SetOrder( int order ) { _order = order; };
void SetNq( int nq ) { _nq = nq; };
void SetPoints( blaze::DynamicVector<blaze::DynamicVector<double>> points ) { _points = points; };
void SetWeights( blaze::DynamicVector<double> weights ) { _weights = weights; };
void SetConnectivity( blaze::DynamicVector<blaze::DynamicVector<int>> connectivity ) { _connectivity = connectivity; };
// Getter
std::string GetName() { return _name; };
int GetOrder() { return _order; };
int GetNq() { return _nq; };
blaze::DynamicVector<blaze::DynamicVector<double>> GetPoints() { return _points; };
blaze::DynamicVector<double> GetWeights() { return _weights; };
blaze::DynamicVector<blaze::DynamicVector<int>> GetConnectivity() { return _connectivity; };
protected:
std::string _name;
int _order;
int _nq;
blaze::DynamicVector<blaze::DynamicVector<double>> _points;
blaze::DynamicVector<double> _weights;
blaze::DynamicVector<blaze::DynamicVector<int>> _connectivity;
};
// Quadrature Hub
Quadrature* GetQuadrature( std::string name, int order );
#endif // QUADRATURE_H
#include "catch.hpp"
#include "factorial_example.h"
TEST_CASE( "Factorials are computed", "[factorial]" ) {
REQUIRE( Factorial( 0 ) == 1 );
REQUIRE( Factorial( 1 ) == 1 );
REQUIRE( Factorial( 2 ) == 2 );
REQUIRE( Factorial( 3 ) == 6 );
REQUIRE( Factorial( 10 ) == 3628800 );
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment