New Angular Quadrature types added into OpenSn.

framework/math/quadratures/angular/angular_quadrature.h

@@ -11,27 +11,36 @@ namespace opensn
 {
 struct QuadraturePointPhiTheta;
 
+/// Angular quadrature type identifier.
 enum class AngularQuadratureType
 {
   ProductQuadrature = 1,
   SLDFEsq = 2,
   LebedevQuadrature = 3,
+  TriangularQuadrature = 4,
 };
 
+/// Quadrature point in spherical coordinates.
 struct QuadraturePointPhiTheta
 {
+  /// Azimuthal angle.
   double phi = 0.0;
+  /// Polar angle.
   double theta = 0.0;
 
   QuadraturePointPhiTheta(const double phi, const double theta) : phi(phi), theta(theta) {}
 };
 
+/// Base class for angular quadratures used in discrete ordinates transport calculations.
 class AngularQuadrature
 {
 public:
+  /// Spherical harmonic indices for moment-to-direction mappings.
   struct HarmonicIndices
   {
+    /// Degree of the spherical harmonic.
     unsigned int ell = 0;
+    /// Order of the spherical harmonic.
     int m = 0;
 
     HarmonicIndices() = default;
@@ -51,40 +60,46 @@ class AngularQuadrature
   {
   }
 
+  /// Discrete-to-moment operator matrix.
   std::vector<std::vector<double>> d2m_op_;
+  /// Moment-to-discrete operator matrix.
   std::vector<std::vector<double>> m2d_op_;
+  /// Mapping from moment index to spherical harmonic indices.
   std::vector<HarmonicIndices> m_to_ell_em_map_;
+  /// Quadrature type identifier.
   AngularQuadratureType type_;
+  /// Spatial dimension of the quadrature.
   unsigned int dimension_;
+  /// Maximum scattering order for moment calculations.
   unsigned int scattering_order_;
 
-  /// Populates a map of moment m to the Spherical Harmonic indices required.
+  /// Populate the map of moment index to spherical harmonic indices.
   void MakeHarmonicIndices();
 
 public:
   virtual ~AngularQuadrature() = default;
 
-  /// Computes the discrete to moment operator.
+  /// Compute the discrete-to-moment operator.
   void BuildDiscreteToMomentOperator();
 
-  /// Computes the moment to discrete operator.
+  /// Compute the moment-to-discrete operator.
   void BuildMomentToDiscreteOperator();
 
   /**
-   * Returns a reference to the precomputed d2m operator. The operator is accessed as [m][d], where
-   * m is the moment index and d is the direction index.
+   * Return a reference to the precomputed discrete-to-moment operator.
+   *
+   * The operator is accessed as [d][m], where m is the moment index and d is the direction index.
    */
   std::vector<std::vector<double>> const& GetDiscreteToMomentOperator() const;
 
   /**
-   * Returns a reference to the precomputed m2d operator. where m is the moment index and d is the
-   * direction index.
+   * Return a reference to the precomputed moment-to-discrete operator.
+   *
+   * The operator is accessed as [d][m], where m is the moment index and d is the direction index.
    */
   std::vector<std::vector<double>> const& GetMomentToDiscreteOperator() const;
 
-  /**
-   * Returns a reference to the precomputed harmonic index map.
-   */
+  /// Return a reference to the precomputed harmonic index map.
   const std::vector<HarmonicIndices>& GetMomentToHarmonicsIndexMap() const;
 
   unsigned int GetDimension() const { return dimension_; }
@@ -95,8 +110,11 @@ class AngularQuadrature
 
   AngularQuadratureType GetType() const { return type_; }
 
+  /// Quadrature point abscissae in spherical coordinates.
   std::vector<QuadraturePointPhiTheta> abscissae;
+  /// Quadrature weights.
   std::vector<double> weights;
+  /// Quadrature point direction vectors in Cartesian coordinates.
   std::vector<Vector3> omegas;
 };
 

framework/math/quadratures/angular/lebedev_quadrature.cc

@@ -101,4 +101,94 @@ LebedevQuadrature3DXYZ::LoadFromOrder(unsigned int quadrature_order, bool verbos
   }
 }
 
+LebedevQuadrature2DXY::LebedevQuadrature2DXY(unsigned int quadrature_order,
+                                             unsigned int scattering_order,
+                                             bool verbose)
+  : AngularQuadrature(AngularQuadratureType::LebedevQuadrature, 2, scattering_order)
+{
+  LoadFromOrder(quadrature_order, verbose);
+  MakeHarmonicIndices();
+  BuildDiscreteToMomentOperator();
+  BuildMomentToDiscreteOperator();
+}
+
+void
+LebedevQuadrature2DXY::LoadFromOrder(unsigned int quadrature_order, bool verbose)
+{
+  abscissae.clear();
+  weights.clear();
+  omegas.clear();
+
+  // Get points from LebedevOrders
+  const auto& points = LebedevOrders::GetOrderPoints(quadrature_order);
+
+  std::stringstream ostr;
+  double weight_sum = 0.0;
+
+  // Tolerance for determining if z is approximately zero (on the equator)
+  const double z_tolerance = 1.0e-12;
+
+  for (const auto& point : points)
+  {
+    const double x = point.x;
+    const double y = point.y;
+    const double z = point.z;
+    double w = point.weight;
+
+    // Skip points with z < 0 (lower hemisphere)
+    if (z < -z_tolerance)
+      continue;
+
+    // For points on the equator (z ≈ 0), halve the weight
+    if (std::fabs(z) <= z_tolerance)
+      w *= 0.5;
+
+    // Calculate phi and theta from x, y, z
+    const double r = std::sqrt(x * x + y * y + z * z);
+    const double theta = std::acos(z / r);
+    double phi = std::atan2(y, x);
+    if (phi < 0.0)
+      phi += 2.0 * M_PI;
+
+    // Create the point
+    QuadraturePointPhiTheta qpoint(phi, theta);
+    abscissae.push_back(qpoint);
+
+    // Create the direction vector
+    Vector3 omega{x / r, y / r, z / r};
+    omegas.push_back(omega);
+
+    // Store the weight (Doubled to renormalize from 0.5 to 1.0)
+    weights.push_back(2.0 * w);
+    weight_sum += 2.0 * w;
+
+    if (verbose)
+    {
+      ostr << "Varphi=" << std::fixed << std::setprecision(2) << qpoint.phi * 180.0 / M_PI
+           << " Theta=" << std::fixed << std::setprecision(2) << qpoint.theta * 180.0 / M_PI
+           << " Weight=" << std::scientific << std::setprecision(3) << w << '\n';
+    }
+  }
+
+  if (verbose)
+  {
+    log.Log() << "Loaded " << abscissae.size()
+              << " Lebedev 2D quadrature points (upper hemisphere) "
+              << "from quadrature order " << quadrature_order;
+    log.Log() << ostr.str() << "\n"
+              << "Weight sum=" << weight_sum;
+  }
+
+  // Check weight sum - should be 1.0 for upper hemisphere
+  const double expected_sum = 1.0;
+  if (std::fabs(weight_sum - expected_sum) > 1.0e-10)
+  {
+    if (verbose)
+    {
+      log.Log() << "Warning: Sum of weights differs from expected value 1.0.";
+      log.Log() << "Expected: " << expected_sum << ", Actual: " << weight_sum;
+    }
+  }
+}
+
 } // namespace opensn

framework/math/quadratures/angular/lebedev_quadrature.h

@@ -12,7 +12,7 @@ namespace opensn
 {
 
 /**
- * @brief Implementation of Lebedev quadrature for angular integration on the surface of a sphere.
+ * Lebedev quadrature for angular integration on the surface of a sphere.
  *
  * Lebedev quadrature provides a set of points on the surface of a sphere that allows for
  * symmetric and efficient angular integration. This implementation reads point data from
@@ -22,21 +22,55 @@ class LebedevQuadrature3DXYZ : public AngularQuadrature
 {
 public:
   /**
-   * @brief Constructor for Lebedev quadrature.
+   * Construct a Lebedev quadrature.
    *
-   * @param order The order of the Lebedev quadrature set to load
-   * @param verbose Flag to enable verbose output
+   * \param quadrature_order Order of the Lebedev quadrature set to load.
+   * \param scattering_order Scattering order for moment calculations.
+   * \param verbose Flag to enable verbose output.
    */
   LebedevQuadrature3DXYZ(unsigned int quadrature_order,
                          unsigned int scattering_order,
                          bool verbose = false);
 
 private:
   /**
-   * @brief Loads quadrature points for the specified order from predefined data.
+   * Load quadrature points for the specified order from predefined data.
    *
-   * @param order The order to load
-   * @param verbose Flag to enable verbose output
+   * \param quadrature_order Order to load.
+   * \param verbose Flag to enable verbose output.
+   */
+  void LoadFromOrder(unsigned int quadrature_order, bool verbose = false);
+};
+
+/**
+ * 2D Lebedev quadrature for angular integration on the upper hemisphere.
+ *
+ * This is a 2D version of the Lebedev quadrature that only includes points with z >= 0
+ * (i.e., polar angles theta <= pi/2). Points on the equator (z = 0) have their weights halved
+ * since they are shared between the upper and lower hemispheres.
+ */
+class LebedevQuadrature2DXY : public AngularQuadrature
+{
+public:
+  /**
+   * Construct a 2D Lebedev quadrature.
+   *
+   * \param quadrature_order Order of the Lebedev quadrature set to load.
+   * \param scattering_order Scattering order for moment calculations.
+   * \param verbose Flag to enable verbose output.
+   */
+  LebedevQuadrature2DXY(unsigned int quadrature_order,
+                        unsigned int scattering_order,
+                        bool verbose = false);
+
+private:
+  /**
+   * Load quadrature points for the specified order from predefined data.
+   *
+   * Keep only upper hemisphere points (z >= 0) and halve weights for z = 0.
+   *
+   * \param quadrature_order Order to load.
+   * \param verbose Flag to enable verbose output.
    */
   void LoadFromOrder(unsigned int quadrature_order, bool verbose = false);
 };