From 92549a619d41a910d1c6b020725b1e1610aba9bd Mon Sep 17 00:00:00 2001 From: Arthur Date: Mon, 23 Dec 2024 13:10:11 +0100 Subject: [PATCH] Template creation 2 store GaussLegendre weights --- src/demeter/quadrature/: | 73 +++++++++++++++++++++ src/demeter/quadrature/:! | 20 ++++++ src/demeter/quadrature/gauss_legendre.cpp | 30 ++++----- src/demeter/quadrature/gauss_legendre.hpp | 22 ++++--- src/demeter/quadrature/w | 77 +++++++++++++++++++++++ 5 files changed, 197 insertions(+), 25 deletions(-) create mode 100644 src/demeter/quadrature/: create mode 100644 src/demeter/quadrature/:! create mode 100644 src/demeter/quadrature/w diff --git a/src/demeter/quadrature/: b/src/demeter/quadrature/: new file mode 100644 index 0000000..94c3f9e --- /dev/null +++ b/src/demeter/quadrature/: @@ -0,0 +1,73 @@ +#include "gauss_legendre.hpp" + +#include +#include +#include + +// To do +// 1 - Definir une structure qui permettent d'aller chercher +// les bons tableaux selon le nombre d'angles polaires +// 2 - Simplifier la structure car plusieurs fois les memes +// valeurs a un signe pres +// 3 - Automatiser le calcul des racines et des poids pour +// aller jusqua 48 +// Tous les poids sont a diviser par deux +// teta est dans [0,pi] pas besoin des angles teta +// dans [pi,2pi] on fait teta_k+pi avec +// w_{k+pi} = w_k + +namespace Demeter { + + +// N = 2 +template <> +const auto Eigen::ArrayXd::GaussLegendre<2>::test_ << 2.000000; +// nteta = 1 + //ArrayXd costeta(1); + //ArrayXd weights(1); + +// costeta << 0.00000; +// weights << 2.00000; + +//// nteta = 2 +// ArrayXd costeta(2); +// ArrayXd weights(2); +// +// costeta << -0.57735, 0.57735; +// weights << 1.00000, 1.00000; +// +//// nteta = 4 +// ArrayXd costeta(3); +// ArrayXd weights(3); +// +// costeta << -0.339981, 0.339981, +// -0.861136, 0.861136; +// weights << 0.652145, 0.652145, +// 0.347855, 0.347855; + + +// nteta = 6 +// nteta = 8 +// nteta = 10 +// nteta = 12 +// nteta = 14 +// nteta = 16 +// nteta = 18 +// nteta = 20 +// nteta = 22 +// nteta = 24 +// nteta = 26 +// nteta = 28 +// nteta = 30 +// nteta = 32 +// nteta = 34 +// nteta = 36 +// nteta = 38 +// nteta = 40 +// nteta = 42 +// nteta = 44 +// nteta = 46 +// nteta = 48 + +} // namespace Demeter + diff --git a/src/demeter/quadrature/:! b/src/demeter/quadrature/:! new file mode 100644 index 0000000..6a09ca6 --- /dev/null +++ b/src/demeter/quadrature/:! @@ -0,0 +1,20 @@ +#pragma once + + +#include +#include + +#include +#include + +namespace Demeter { + +template +class GaussLegendre { + public: + std::span wgt() const { return {weights_.begin(), weights_.end()}; } + + private: + static const std::array weights_; +}; +} // namespace Demeter diff --git a/src/demeter/quadrature/gauss_legendre.cpp b/src/demeter/quadrature/gauss_legendre.cpp index 4eeaff6..eda4954 100644 --- a/src/demeter/quadrature/gauss_legendre.cpp +++ b/src/demeter/quadrature/gauss_legendre.cpp @@ -5,12 +5,8 @@ #include // To do -// 1 - Definir une structure qui permettent d'aller chercher -// les bons tableaux selon le nombre d'angles polaires -// 2 - Simplifier la structure car plusieurs fois les memes -// valeurs a un signe pres -// 3 - Automatiser le calcul des racines et des poids pour -// aller jusqua 48 +// 1 - Automatiser le calcul des racines et des poids pour +// aller jusqua 48 (python) // Tous les poids sont a diviser par deux // teta est dans [0,pi] pas besoin des angles teta // dans [pi,2pi] on fait teta_k+pi avec @@ -18,15 +14,20 @@ namespace Demeter { -void GaussLegendre::compute_angles() { - // to do -} +// nteta = 2 +template <> +Eigen::ArrayXd GaussLegendre<2>::weights_(1); +template <> +Eigen::ArrayXd GaussLegendre<2>::costeta_(1); +template <> +Eigen::ArrayXd GaussLegendre<2>::teta_(1); +template <> +void GaussLegendre<2>::initialization() { + weights_ << 2.00000; + costeta_ << 0.00000; + teta_ << 0.00000; +} -void GaussLegendre::compute_weights() { - // to do -} - - // nteta = 1 //ArrayXd costeta(1); //ArrayXd weights(1); @@ -75,4 +76,3 @@ void GaussLegendre::compute_weights() { // nteta = 48 } // namespace Demeter - diff --git a/src/demeter/quadrature/gauss_legendre.hpp b/src/demeter/quadrature/gauss_legendre.hpp index ec1d267..5d56ee9 100644 --- a/src/demeter/quadrature/gauss_legendre.hpp +++ b/src/demeter/quadrature/gauss_legendre.hpp @@ -1,24 +1,26 @@ #pragma once - -#include #include namespace Demeter { +template class GaussLegendre { public: + // Define array initializer + static void initialization(); - GaussLegendre(const int& nphi, const int& nteta) - : nphi_(nphi), nteta_(nteta) { - } + // Define getter methods + const Eigen::ArrayXd& getweights() const { return weights_; } + const Eigen::ArrayXd& getcosteta() const { return costeta_; } + const Eigen::ArrayXd& getteta() const { return teta_; } - void compute_angles(); - void compute_weights(); private: - size_t nphi_; // nber of azimutal angles - size_t nteta_; // nber of polar angles - + static Eigen::ArrayXd weights_; + static Eigen::ArrayXd costeta_; + static Eigen::ArrayXd teta_; }; + + } // namespace Demeter diff --git a/src/demeter/quadrature/w b/src/demeter/quadrature/w new file mode 100644 index 0000000..997ee6a --- /dev/null +++ b/src/demeter/quadrature/w @@ -0,0 +1,77 @@ +#include "gauss_legendre.hpp" + +#include +#include +#include + +// To do +// 1 - Definir une structure qui permettent d'aller chercher +// les bons tableaux selon le nombre d'angles polaires +// 2 - Simplifier la structure car plusieurs fois les memes +// valeurs a un signe pres +// 3 - Automatiser le calcul des racines et des poids pour +// aller jusqua 48 +// Tous les poids sont a diviser par deux +// teta est dans [0,pi] pas besoin des angles teta +// dans [pi,2pi] on fait teta_k+pi avec +// w_{k+pi} = w_k + +namespace Demeter { + + +// N = 2 +template <> + +const auto GaussLegendre<2>::weights_ = flux_ = ArrayXd::Ones(2); +const std::array GaussLegendre<2>::wgt_ = { + 1.00000000000000000000000000000000000e+00}; + +// nteta = 1 + //ArrayXd costeta(1); + //ArrayXd weights(1); + +// costeta << 0.00000; +// weights << 2.00000; + +//// nteta = 2 +// ArrayXd costeta(2); +// ArrayXd weights(2); +// +// costeta << -0.57735, 0.57735; +// weights << 1.00000, 1.00000; +// +//// nteta = 4 +// ArrayXd costeta(3); +// ArrayXd weights(3); +// +// costeta << -0.339981, 0.339981, +// -0.861136, 0.861136; +// weights << 0.652145, 0.652145, +// 0.347855, 0.347855; + + +// nteta = 6 +// nteta = 8 +// nteta = 10 +// nteta = 12 +// nteta = 14 +// nteta = 16 +// nteta = 18 +// nteta = 20 +// nteta = 22 +// nteta = 24 +// nteta = 26 +// nteta = 28 +// nteta = 30 +// nteta = 32 +// nteta = 34 +// nteta = 36 +// nteta = 38 +// nteta = 40 +// nteta = 42 +// nteta = 44 +// nteta = 46 +// nteta = 48 + +} // namespace Demeter +