25654184358

Committed 09 May 2026 03:58AM UTC coverage: 75.687%. Remained the same

Build # 25654184358

Build Type

push

github

Committed by

web-flow

Commit Message

Merge pull request #1042 from wdhawkins/ang_quadrature_refactor

Refactoring angular quadrature classes.

Coverage Stats

230 of 278 new or added lines in 27 files covered. (82.73%)

181 existing lines in 12 files now uncovered.

22258 of 29408 relevant lines covered (75.69%)

64739266.89 hits per line

Source File
Press 'n' to go to next uncovered line, 'b' for previous

75.49

/framework/math/quadratures/angular/lebedev_quadrature.cc

// SPDX-FileCopyrightText: 2025 The OpenSn Authors <https://open-sn.github.io/opensn/>
// SPDX-License-Identifier: MIT

#include "framework/math/quadratures/angular/lebedev_quadrature.h"
#include "framework/math/quadratures/angular/lebedev_orders.h"
#include "framework/logging/log.h"
#include "framework/runtime.h"
#include <fstream>
#include <sstream>
#include <iostream>
#include <iomanip>
#include <cmath>
#include <stdexcept>

namespace opensn
{

LebedevQuadrature3DXYZ::LebedevQuadrature3DXYZ(unsigned int quadrature_order,
                                               unsigned int scattering_order,
                                               bool verbose,
                                               OperatorConstructionMethod method)
  : AngularQuadrature(AngularQuadratureType::LEBEDEV_QUADRATURE, 3, scattering_order, method),
    quadrature_order_(quadrature_order)
{
  LoadFromOrder(quadrature_order, verbose);
  MakeHarmonicIndices();
  BuildDiscreteToMomentOperator();
  BuildMomentToDiscreteOperator();
}

void
LebedevQuadrature3DXYZ::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;

  for (const auto& point : points)
  {
    const double x = point.x;
    const double y = point.y;
    const double z = point.z;
    const double w = point.weight;

    // 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
    weights_.push_back(w);
    weight_sum += 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 " << points.size() << " Lebedev quadrature points from quadrature order "
              << quadrature_order;
    log.Log() << ostr.str() << "\n"
              << "Weight sum=" << weight_sum;
  }

  // Check weight sum
  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.";
      log.Log() << "Expected: " << expected_sum << ", Actual: " << weight_sum;
    }

    // Normalize weights
    const double scale_factor = expected_sum / weight_sum;
    for (auto& w : weights_)
      w *= scale_factor;

    if (verbose)
      log.Log() << "Weights have been normalized to sum to 1.";
  }
}

LebedevQuadrature2DXY::LebedevQuadrature2DXY(unsigned int quadrature_order,
                                             unsigned int scattering_order,
                                             bool verbose,
                                             OperatorConstructionMethod method)
  : AngularQuadrature(AngularQuadratureType::LEBEDEV_QUADRATURE, 2, scattering_order, method),
    quadrature_order_(quadrature_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 " << GetNumAngles() << " 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

1	// SPDX-FileCopyrightText: 2025 The OpenSn Authors <https://open-sn.github.io/opensn/>
2	// SPDX-License-Identifier: MIT
3
4	#include "framework/math/quadratures/angular/lebedev_quadrature.h"
5	#include "framework/math/quadratures/angular/lebedev_orders.h"
6	#include "framework/logging/log.h"
7	#include "framework/runtime.h"
8	#include <fstream>
9	#include <sstream>
10	#include <iostream>
11	#include <iomanip>
12	#include <cmath>
13	#include <stdexcept>
14
15	namespace opensn
16	{
17
18	LebedevQuadrature3DXYZ::LebedevQuadrature3DXYZ(unsigned int quadrature_order,	6✔
19	unsigned int scattering_order,
20	bool verbose,
21	OperatorConstructionMethod method)	6✔
22	: AngularQuadrature(AngularQuadratureType::LEBEDEV_QUADRATURE, 3, scattering_order, method),
23	quadrature_order_(quadrature_order)	6✔
24	{
25	LoadFromOrder(quadrature_order, verbose);	6✔
26	MakeHarmonicIndices();	6✔
27	BuildDiscreteToMomentOperator();	6✔
28	BuildMomentToDiscreteOperator();	6✔
29	}	6✔
30
31	void
32	LebedevQuadrature3DXYZ::LoadFromOrder(unsigned int quadrature_order, bool verbose)	6✔
33	{
34	abscissae_.clear();	6✔
35	weights_.clear();	6✔
36	omegas_.clear();	6✔
37
38	// Get points from LebedevOrders
39	const auto& points = LebedevOrders::GetOrderPoints(quadrature_order);	6✔
40
41	std::stringstream ostr;	6✔
42	double weight_sum = 0.0;	6✔
43
44	for (const auto& point : points)	378✔
45	{
46	const double x = point.x;	372✔
47	const double y = point.y;	372✔
48	const double z = point.z;	372✔
49	const double w = point.weight;	372✔
50
51	// Calculate phi and theta from x, y, z
52	const double r = std::sqrt(x * x + y * y + z * z);	372✔
53	const double theta = std::acos(z / r);	372✔
54	double phi = std::atan2(y, x);	372✔
55	if (phi < 0.0)	372✔
56	phi += 2.0 * M_PI;	162✔
57
58	// Create the point
59	QuadraturePointPhiTheta qpoint(phi, theta);	372✔
60	abscissae_.push_back(qpoint);	372✔
61
62	// Create the direction vector
63	Vector3 omega{x / r, y / r, z / r};	372✔
64	omegas_.push_back(omega);	372✔
65
66	// Store the weight
67	weights_.push_back(w);	372✔
68	weight_sum += w;	372✔
69
70	if (verbose)	372✔
71	{
UNCOV 72	ostr << "Varphi=" << std::fixed << std::setprecision(2) << qpoint.phi * 180.0 / M_PI	×
73	<< " Theta=" << std::fixed << std::setprecision(2) << qpoint.theta * 180.0 / M_PI	×
74	<< " Weight=" << std::scientific << std::setprecision(3) << w << '\n';	×
75	}
76	}
77
78	if (verbose)	6✔
79	{
UNCOV 80	log.Log() << "Loaded " << points.size() << " Lebedev quadrature points from quadrature order "	×
81	<< quadrature_order;	×
82	log.Log() << ostr.str() << "\n"	×
83	<< "Weight sum=" << weight_sum;	×
84	}
85
86	// Check weight sum
87	const double expected_sum = 1.0;	6✔
88	if (std::fabs(weight_sum - expected_sum) > 1.0e-10)	6✔
89	{
UNCOV 90	if (verbose)	×
91	{
UNCOV 92	log.Log() << "Warning: Sum of weights differs from expected value 1.";	×
93	log.Log() << "Expected: " << expected_sum << ", Actual: " << weight_sum;	×
94	}
95
96	// Normalize weights
UNCOV 97	const double scale_factor = expected_sum / weight_sum;	×
NEW 98	for (auto& w : weights_)	×
99	w *= scale_factor;	×
100
UNCOV 101	if (verbose)	×
102	log.Log() << "Weights have been normalized to sum to 1.";	×
103	}
104	}	6✔
105
106	LebedevQuadrature2DXY::LebedevQuadrature2DXY(unsigned int quadrature_order,	1✔
107	unsigned int scattering_order,
108	bool verbose,
109	OperatorConstructionMethod method)	1✔
110	: AngularQuadrature(AngularQuadratureType::LEBEDEV_QUADRATURE, 2, scattering_order, method),
111	quadrature_order_(quadrature_order)	1✔
112	{
113	LoadFromOrder(quadrature_order, verbose);	1✔
114	MakeHarmonicIndices();	1✔
115	BuildDiscreteToMomentOperator();	1✔
116	BuildMomentToDiscreteOperator();	1✔
117	}	1✔
118
119	void
120	LebedevQuadrature2DXY::LoadFromOrder(unsigned int quadrature_order, bool verbose)	1✔
121	{
122	abscissae_.clear();	1✔
123	weights_.clear();	1✔
124	omegas_.clear();	1✔
125
126	// Get points from LebedevOrders
127	const auto& points = LebedevOrders::GetOrderPoints(quadrature_order);	1✔
128
129	std::stringstream ostr;	1✔
130	double weight_sum = 0.0;	1✔
131
132	// Tolerance for determining if z is approximately zero (on the equator)
133	const double z_tolerance = 1.0e-12;	1✔
134
135	for (const auto& point : points)	7✔
136	{
137	const double x = point.x;	6✔
138	const double y = point.y;	6✔
139	const double z = point.z;	6✔
140	double w = point.weight;	6✔
141
142	// Skip points with z < 0 (lower hemisphere)
143	if (z < -z_tolerance)	6✔
144	continue;	1✔
145
146	// For points on the equator (z ≈ 0), halve the weight
147	if (std::fabs(z) <= z_tolerance)	5✔
148	w *= 0.5;	4✔
149
150	// Calculate phi and theta from x, y, z
151	const double r = std::sqrt(x * x + y * y + z * z);	5✔
152	const double theta = std::acos(z / r);	5✔
153	double phi = std::atan2(y, x);	5✔
154	if (phi < 0.0)	5✔
155	phi += 2.0 * M_PI;	1✔
156
157	// Create the point
158	QuadraturePointPhiTheta qpoint(phi, theta);	5✔
159	abscissae_.push_back(qpoint);	5✔
160
161	// Create the direction vector
162	Vector3 omega{x / r, y / r, z / r};	5✔
163	omegas_.push_back(omega);	5✔
164
165	// Store the weight (Doubled to renormalize from 0.5 to 1.0)
166	weights_.push_back(2.0 * w);	5✔
167	weight_sum += 2.0 * w;	5✔
168
169	if (verbose)	5✔
170	{
UNCOV 171	ostr << "Varphi=" << std::fixed << std::setprecision(2) << qpoint.phi * 180.0 / M_PI	×
UNCOV 172	<< " Theta=" << std::fixed << std::setprecision(2) << qpoint.theta * 180.0 / M_PI	×
173	<< " Weight=" << std::scientific << std::setprecision(3) << w << '\n';	×
174	}
175	}
176
177	if (verbose)	1✔
178	{
NEW 179	log.Log() << "Loaded " << GetNumAngles() << " Lebedev 2D quadrature points (upper hemisphere) "	×
UNCOV 180	<< "from quadrature order " << quadrature_order;	×
181	log.Log() << ostr.str() << "\n"	×
182	<< "Weight sum=" << weight_sum;	×
183	}
184
185	// Check weight sum - should be 1.0 for upper hemisphere
186	const double expected_sum = 1.0;	1✔
187	if (std::fabs(weight_sum - expected_sum) > 1.0e-10)	1✔
188	{
UNCOV 189	if (verbose)	×
190	{
UNCOV 191	log.Log() << "Warning: Sum of weights differs from expected value 1.0.";	×
192	log.Log() << "Expected: " << expected_sum << ", Actual: " << weight_sum;	×
193	}
194	}
195	}	1✔
196
197	} // namespace opensn

Open-Sn / opensn / 25654184358

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous