Improving and testing Legendre quadrature

kennethassogba · Dec 26, 2024 · 86af7c4 · 86af7c4
1 parent 93577b9
commit 86af7c4
Show file tree

Hide file tree

Showing 4 changed files with 78 additions and 31 deletions.
diff --git a/examples/main.cpp b/examples/main.cpp
@@ -5,13 +5,29 @@
 
 #include "demeter/model.hpp"
 #include "demeter/solve.hpp"
+#include "demeter/quadrature.hpp"
 
 int main() {
   using namespace Demeter;
   using namespace Eigen;
 
   CheckOpenMP();
 
+  // Define variables
+  constexpr size_t nteta   = 48;
+  constexpr size_t ngroups = 2;
+  constexpr size_t naniso  = 0;
+  constexpr size_t nmoms   = (naniso+1)*(naniso+1);
+
+  // Test sum of gaussian weights
+  GaussLegendre<nteta> gauss;
+  gauss.getweights();
+
+  ArrayXd weights(nteta/2);
+  weights = gauss.getweights();
+  std::cout << weights;
+  std::cout << "\n";
+
   // Test linear solvers
   Eigen::Matrix2f A, b;
   A << 2, -1, -1, 3; 
@@ -24,11 +40,7 @@ int main() {
   std::cout << '\n';
 
 
-  // Define variables
-  size_t ngroups = 2;
-  size_t naniso  = 0;
-  size_t nmoms   = (naniso+1)*(naniso+1);
-
+
   // Define cross sections
   ArrayXd  sigma_t(ngroups);
   ArrayXXd sigma_s(ngroups,ngroups);

diff --git a/src/demeter/quadrature.hpp b/src/demeter/quadrature.hpp
@@ -1,3 +1,3 @@
 #pragma once
 
-#include "solve/gauss_legendre.hpp"
+#include "quadrature/gauss_legendre.hpp"
diff --git a/src/demeter/quadrature/gauss_legendre.cpp b/src/demeter/quadrature/gauss_legendre.cpp
@@ -7,6 +7,7 @@
 namespace Demeter {
 
 // nteta = 2
+//template class GaussLegendre<2>;
 template <>
 Eigen::ArrayXd GaussLegendre<2>::weights_(1);
 template <>
@@ -17,9 +18,10 @@ template <>
 void GaussLegendre<2>::initialization() {
     weights_ << 1.0000000;
     costeta_ << 0.57735027;
-    teta_ << 0.00000;
+    teta_ = costeta_.acos() * 180.0 / M_PI;
 }
 // nteta = 4
+template class GaussLegendre<4>;
 template <>
 Eigen::ArrayXd GaussLegendre<4>::weights_(2);
 template <>
@@ -30,9 +32,10 @@ template <>
 void GaussLegendre<4>::initialization() {
     weights_ << 0.65214515, 0.34785485;
     costeta_ << 0.33998104, 0.86113631;
-    teta_ << 0.00000;
+    teta_ = costeta_.acos() * 180.0 / M_PI;
 }
 // nteta = 6
+template class GaussLegendre<6>;
 template <>
 Eigen::ArrayXd GaussLegendre<6>::weights_(3);
 template <>
@@ -43,9 +46,10 @@ template <>
 void GaussLegendre<6>::initialization() {
     weights_ << 0.46791393, 0.36076157, 0.17132449;
     costeta_ << 0.23861919, 0.66120939, 0.93246951;
-    teta_ << 0.00000;
+    teta_ << 0.23861919, 0.66120939, 0.93246951;
 }
 // nteta = 8
+template class GaussLegendre<8>;
 template <>
 Eigen::ArrayXd GaussLegendre<8>::weights_(4);
 template <>
@@ -56,9 +60,10 @@ template <>
 void GaussLegendre<8>::initialization() {
     weights_ << 0.36268378, 0.31370665, 0.22238103, 0.10122854;
     costeta_ << 0.18343464, 0.52553241, 0.79666648, 0.96028986;
-    teta_ << 0.00000;
+    teta_ = costeta_.acos() * 180.0 / M_PI;
 }
 // nteta = 10
+template class GaussLegendre<10>;
 template <>
 Eigen::ArrayXd GaussLegendre<10>::weights_(5);
 template <>
@@ -69,9 +74,10 @@ template <>
 void GaussLegendre<10>::initialization() {
     weights_ << 0.29552422, 0.26926672, 0.21908636, 0.14945135, 0.06667134;
     costeta_ << 0.14887434, 0.43339539, 0.67940957, 0.86506337, 0.97390653;
-    teta_ << 0.00000;
+    teta_ = costeta_.acos() * 180.0 / M_PI;
 }
 // nteta = 12
+template class GaussLegendre<12>;
 template <>
 Eigen::ArrayXd GaussLegendre<12>::weights_(6);
 template <>
@@ -84,9 +90,10 @@ void GaussLegendre<12>::initialization() {
                 0.04717534;
     costeta_ << 0.12523341, 0.3678315,  0.58731795, 0.76990267, 0.90411726,
                 0.98156063;
-    teta_ << 0.00000;
+    teta_ = costeta_.acos() * 180.0 / M_PI;
 }
 // nteta = 14
+template class GaussLegendre<14>;
 template <>
 Eigen::ArrayXd GaussLegendre<14>::weights_(7);
 template <>
@@ -99,9 +106,10 @@ void GaussLegendre<14>::initialization() {
                 0.08015809, 0.03511946;
     costeta_ << 0.10805495, 0.31911237, 0.51524864, 0.6872929, 0.82720132, 
                 0.92843488, 0.98628381;
-    teta_ << 0.00000;
+    teta_ = costeta_.acos() * 180.0 / M_PI;
 }
 // nteta = 16
+template class GaussLegendre<16>;
 template <>
 Eigen::ArrayXd GaussLegendre<16>::weights_(8);
 template <>
@@ -114,9 +122,10 @@ void GaussLegendre<16>::initialization() {
                 0.09515851, 0.06225352, 0.02715246;
     costeta_ << 0.09501251, 0.28160355, 0.45801678, 0.61787624, 0.75540441, 
                 0.8656312, 0.94457502, 0.98940093;
-    teta_ << 0.00000;
+    teta_ = costeta_.acos() * 180.0 / M_PI;
 }
 // nteta = 18
+template class GaussLegendre<18>;
 template <>
 Eigen::ArrayXd GaussLegendre<18>::weights_(9);
 template <>
@@ -128,9 +137,10 @@ void GaussLegendre<18>::initialization() {
     weights_ << 0.16914238, 0.16427648, 0.15468468, 0.14064291, 0.12255521, 
                 0.10094204, 0.07642573, 0.04971455, 0.02161601;
     costeta_ << 0.00000;
-    teta_ << 0.00000;
+    teta_ = costeta_.acos() * 180.0 / M_PI;
 }
 // nteta = 20
+template class GaussLegendre<20>;
 template <>
 Eigen::ArrayXd GaussLegendre<20>::weights_(10);
 template <>
@@ -143,9 +153,10 @@ void GaussLegendre<20>::initialization() {
                 0.10193012, 0.08327674, 0.06267205, 0.04060143, 0.01761401;
     costeta_ << 0.07652652, 0.22778585, 0.37370609, 0.510867,   0.63605368,
                 0.74633191, 0.83911697, 0.91223443, 0.96397193, 0.9931286;
-    teta_ << 0.00000;
+    teta_ = costeta_.acos() * 180.0 / M_PI;
 }
 // nteta = 22
+template class GaussLegendre<22>;
 template <>
 Eigen::ArrayXd GaussLegendre<22>::weights_(11);
 template <>
@@ -160,9 +171,10 @@ void GaussLegendre<22>::initialization() {
     costeta_ << 0.06973927, 0.20786043, 0.34193582, 0.46935584, 0.5876404,  
                 0.69448726, 0.78781681, 0.86581258, 0.92695677, 0.9700605,  
                 0.99429459;
-    teta_ << 0.00000;
+    teta_ = costeta_.acos() * 180.0 / M_PI;
 }
 // nteta = 24
+template class GaussLegendre<24>;
 template <>
 Eigen::ArrayXd GaussLegendre<24>::weights_(12);
 template <>
@@ -177,9 +189,10 @@ void GaussLegendre<24>::initialization() {
     costeta_ << 0.06405689, 0.19111887, 0.31504268, 0.43379351, 0.54542147, 
                 0.64809365, 0.74012419, 0.82000199, 0.88641553, 0.93827455, 
                 0.97472856, 0.99518722;
-    teta_ << 0.00000;
+    teta_ = costeta_.acos() * 180.0 / M_PI;
 }
 // nteta = 26
+template class GaussLegendre<26>;
 template <>
 Eigen::ArrayXd GaussLegendre<26>::weights_(13);
 template <>
@@ -194,9 +207,10 @@ void GaussLegendre<26>::initialization() {
     costeta_ << 0.05923009, 0.17685882, 0.29200484, 0.40305176, 0.50844071,
                 0.60669229, 0.69642726, 0.77638595, 0.84544594, 0.90263786,
                 0.94715907, 0.97838545, 0.9958857;
-    teta_ << 0.00000;
+    teta_ = costeta_.acos() * 180.0 / M_PI;
 }
 // nteta = 28
+template class GaussLegendre<28>;
 template <>
 Eigen::ArrayXd GaussLegendre<28>::weights_(14);
 template <>
@@ -211,9 +225,10 @@ void GaussLegendre<28>::initialization() {
     costeta_ << 0.05507929, 0.16456928, 0.27206163, 0.37625152, 0.47587422, 
                 0.56972047, 0.65665109, 0.73561088, 0.80564137, 0.86589252, 
                 0.91563303, 0.95425928, 0.98130317, 0.9964425;
-    teta_ << 0.00000;
+    teta_ = costeta_.acos() * 180.0 / M_PI;
 }
 // nteta = 30
+template class GaussLegendre<30>;
 template <>
 Eigen::ArrayXd GaussLegendre<30>::weights_(15);
 template <>
@@ -228,9 +243,10 @@ void GaussLegendre<30>::initialization() {
     costeta_ << 0.05147184, 0.15386991, 0.25463693, 0.35270473, 0.44703377, 
                 0.53662415, 0.62052618, 0.69785049, 0.76777743, 0.82956576, 
                 0.88256054, 0.92620005, 0.96002186, 0.98366812, 0.99689348;
-    teta_ << 0.00000;
+    teta_ = costeta_.acos() * 180.0 / M_PI;
 }
 // nteta = 32
+template class GaussLegendre<32>;
 template <>
 Eigen::ArrayXd GaussLegendre<32>::weights_(16);
 template <>
@@ -247,9 +263,10 @@ void GaussLegendre<32>::initialization() {
                 0.50689991, 0.58771576, 0.66304427, 0.73218212, 0.7944838 , 
                 0.84936761, 0.89632116, 0.93490608, 0.96476226, 0.98561151, 
                 0.99726386;
-    teta_ << 0.00000;
+    teta_ = costeta_.acos() * 180.0 / M_PI;
 }
 // nteta = 34
+template class GaussLegendre<34>;
 template <>
 Eigen::ArrayXd GaussLegendre<34>::weights_(17);
 template <>
@@ -266,9 +283,10 @@ void GaussLegendre<34>::initialization() {
                 0.48010655, 0.5578755 , 0.63102173, 0.69893911, 0.76106488,
                 0.81688423, 0.86593464, 0.90780968, 0.9421624 , 0.96870826,
                 0.98722782, 0.99757175;
-    teta_ << 0.00000;
+    teta_ = costeta_.acos() * 180.0 / M_PI;
 }
 // nteta = 36
+template class GaussLegendre<36>;
 template <>
 Eigen::ArrayXd GaussLegendre<36>::weights_(18);
 template <>
@@ -285,9 +303,10 @@ void GaussLegendre<36>::initialization() {
                 0.4558639 , 0.53068029, 0.60156766, 0.66800124, 0.72948917, 
                 0.78557623, 0.83584717, 0.8799298 , 0.91749777, 0.94827298, 
                 0.97202769, 0.98858648, 0.99783046;
-    teta_ << 0.00000;
+    teta_ = costeta_.acos() * 180.0 / M_PI;
 }
 // nteta = 38
+template class GaussLegendre<38>;
 template <>
 Eigen::ArrayXd GaussLegendre<38>::weights_(19);
 template <>
@@ -304,9 +323,10 @@ void GaussLegendre<38>::initialization() {
                 0.43384717, 0.50583472, 0.57445602, 0.63925442, 0.69979868, 
                 0.7556859 , 0.80654417, 0.85203502, 0.89185574, 0.92574133, 
                 0.95346633, 0.97484633, 0.98973945, 0.99804993;
-    teta_ << 0.00000;
+    teta_ = costeta_.acos() * 180.0 / M_PI;
 }
 // nteta = 40
+template class GaussLegendre<40>;
 template <>
 Eigen::ArrayXd GaussLegendre<40>::weights_(20);
 template <>
@@ -323,9 +343,10 @@ void GaussLegendre<40>::initialization() {
                 0.4137792 , 0.4830758 , 0.54946713, 0.61255389, 0.67195668, 
                 0.72731826, 0.77830565, 0.82461223, 0.8659595 , 0.90209881, 
                 0.93281281, 0.95791682, 0.97725995, 0.99072624, 0.99823771;
-    teta_ << 0.00000;
+    teta_ = costeta_.acos() * 180.0 / M_PI;
 }
 // nteta = 42
+template class GaussLegendre<42>;
 template <>
 Eigen::ArrayXd GaussLegendre<42>::weights_(21);
 template <>
@@ -344,9 +365,10 @@ void GaussLegendre<42>::initialization() {
                 0.70049459, 0.75127994, 0.79796205, 0.84028598, 0.87802057, 
                 0.91095972, 0.93892356, 0.96175937, 0.97934251, 0.99157729, 
                 0.99839962;
-    teta_ << 0.00000;
+    teta_ = costeta_.acos() * 180.0 / M_PI;
 }
 // nteta = 44
+template class GaussLegendre<44>;
 template <>
 Eigen::ArrayXd GaussLegendre<44>::weights_(22);
 template <>
@@ -365,9 +387,10 @@ void GaussLegendre<44>::initialization() {
                 0.67518607, 0.72553105, 0.77226148, 0.81514454, 0.8539666 , 
                 0.88853424, 0.91867526, 0.94423951, 0.96509965, 0.98115183, 
                 0.99231639, 0.9985402;
-    teta_ << 0.00000;
+    teta_ = costeta_.acos() * 180.0 / M_PI;
 }
 // nteta = 46
+template class GaussLegendre<46>;
 template <>
 Eigen::ArrayXd GaussLegendre<46>::weights_(23);
 template <>
@@ -386,9 +409,10 @@ void GaussLegendre<46>::initialization() {
                 0.65133485, 0.70106951, 0.74760536, 0.79073006, 0.83024684, 
                 0.86597539, 0.89775271, 0.9254338 , 0.94889236, 0.96802139, 
                 0.98273367, 0.99296235, 0.99866304;
-    teta_ << 0.00000;
+    teta_ = costeta_.acos() * 180.0 / M_PI;
 }
 // nteta = 48
+template class GaussLegendre<48>;
 template <>
 Eigen::ArrayXd GaussLegendre<48>::weights_(24);
 template <>
@@ -407,6 +431,6 @@ void GaussLegendre<48>::initialization() {
                 0.6288674 , 0.67787238, 0.72403413, 0.76715903, 0.8070662 , 
                 0.84358826, 0.87657202, 0.90587914, 0.93138669, 0.9529877 , 
                 0.97059159, 0.98412458, 0.99353017, 0.99877101;
-    teta_ << 0.00000;
+    teta_ = costeta_.acos() * 180.0 / M_PI;
 }
 } // namespace Demeter
diff --git a/src/demeter/quadrature/gauss_legendre.hpp b/src/demeter/quadrature/gauss_legendre.hpp
@@ -2,6 +2,13 @@
 
 #include <Eigen/Dense>
 
+// to do
+
+// fonction qui permet assurer 
+// somme des poids vaut 1
+// permet assurer angles dans 0 pi
+
+
 namespace Demeter {
 
  /**
@@ -19,10 +26,14 @@ template <std::size_t nteta>
 class GaussLegendre {
  public:
 
+  GaussLegendre() {
+      initialization();
+  }
+
   /**
   * @brief Allocate the arrays to store the data.
   */
-  static void initialization();
+  void initialization();
 
   // Define getter methods
   const Eigen::ArrayXd& getweights() const { return weights_; }