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

90.91

/framework/math/quadratures/angular/triangular_quadrature.cc

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

#include "framework/math/quadratures/angular/triangular_quadrature.h"
#include "framework/math/quadratures/gausslegendre_quadrature.h"
#include "framework/logging/log.h"
#include "framework/runtime.h"
#include <cmath>

namespace opensn
{

void
GLCTriangularQuadrature3DXYZ::AssembleTriangularCosines(
  const std::vector<std::vector<double>>& azimuthal_per_polar,
  const std::vector<double>& polar,
  const std::vector<std::vector<double>>& wts_per_polar,
  bool verbose)
{
  const size_t Np = polar.size();

  polar_ang_ = polar;
  azimuthal_per_polar_ = azimuthal_per_polar;

  if (verbose)
  {
    log.Log() << "Polar angles:";
    for (const auto& ang : polar_ang_)
      log.Log() << ang;
  }

  // Create angle pairs with varying azimuthal angles per polar level
  map_directions_.clear();
  for (unsigned int j = 0; j < Np; ++j)
    map_directions_.emplace(j, std::vector<unsigned int>());

  abscissae_.clear();
  weights_.clear();
  weight_sum_ = 0.0;

  unsigned int direction_index = 0;
  for (unsigned int j = 0; j < Np; ++j)
  {
    const size_t Na_j = azimuthal_per_polar[j].size();

    if (verbose)
    {
      log.Log() << "Polar level " << j << " (theta=" << polar[j] * 180.0 / M_PI << " deg): " << Na_j
                << " azimuthal angles";
    }

    for (unsigned int i = 0; i < Na_j; ++i)
    {
      map_directions_[j].emplace_back(direction_index);

      const auto abscissa = QuadraturePointPhiTheta(azimuthal_per_polar[j][i], polar[j]);
      abscissae_.emplace_back(abscissa);

      const double weight = wts_per_polar[j][i];
      weights_.emplace_back(weight);
      weight_sum_ += weight;

      ++direction_index;
    }
  }

  // Create omega list
  omegas_.clear();
  for (const auto& qpoint : abscissae_)
  {
    Vector3 new_omega;
    new_omega.x = std::sin(qpoint.theta) * std::cos(qpoint.phi);
    new_omega.y = std::sin(qpoint.theta) * std::sin(qpoint.phi);
    new_omega.z = std::cos(qpoint.theta);

    omegas_.emplace_back(new_omega);

    if (verbose)
      log.Log() << "Quadrature angle=" << new_omega.PrintStr();
  }

  // Normalize weights to 1.0
  for (auto& w : weights_)
    w /= weight_sum_;

  // Compute and store sum of weights after normalization
  weight_sum_ = 0.0;
  for (auto& w : weights_)
    weight_sum_ += w;
}

GLCTriangularQuadrature3DXYZ::GLCTriangularQuadrature3DXYZ(unsigned int Npolar,
                                                           unsigned int scattering_order,
                                                           bool verbose,
                                                           OperatorConstructionMethod method)
  : TriangularQuadrature(3, scattering_order, method)
{
  if (Npolar % 2 != 0)
    throw std::invalid_argument("GLCTriangularQuadrature3DXYZ: Npolar must be even.");

  // Compute maximum number of azimuthal angles at the equator
  // For the triangular reduction pattern with 4 fewer azimuthal angles per level from equator,
  // this ensures exactly 4 azimuthal angles at the poles (1 per octant)
  const unsigned int Nazimuthal = 2 * Npolar;

  n_polar_ = Npolar;
  n_azimuthal_ = Nazimuthal;

  GaussLegendreQuadrature gl_polar(Npolar);

  // Create polar angles
  std::vector<double> polar;
  polar.reserve(Npolar);
  for (unsigned int j = 0; j < Npolar; ++j)
    polar.emplace_back(M_PI - std::acos(gl_polar.qpoints[j][0]));

  // Calculate number of azimuthal angles per polar level
  // At equator: Nazimuthal angles
  // Each level away from equator: 4 less azimuthal angles (1 less per octant)
  // This maintains octant symmetry in the triangular pattern
  std::vector<unsigned int> num_azimuthal_per_level(Npolar);
  for (unsigned int j = 0; j < Npolar; ++j)
  {
    // Distance from equator in terms of polar levels
    // For even Npolar, the equatorial levels are at Npolar/2 - 1 and Npolar/2
    // Levels from equator: how many steps from the nearest equatorial level
    unsigned int levels_from_equator = 0;
    if (j < Npolar / 2)
    {
      // Upper hemisphere: distance from equatorial level (Npolar/2 - 1)
      levels_from_equator = (Npolar / 2 - 1) - j;
    }
    else
    {
      // Lower hemisphere: distance from equatorial level (Npolar/2)
      levels_from_equator = j - (Npolar / 2);
    }
    // Reduce by 4 azimuthal angles per level (1 per octant) to maintain symmetry
    num_azimuthal_per_level[j] = Nazimuthal - 4 * levels_from_equator;
  }

  // Create azimuthal angles and weights per polar level
  std::vector<std::vector<double>> azimuthal_per_polar(Npolar);
  std::vector<std::vector<double>> wts_per_polar(Npolar);

  for (unsigned int j = 0; j < Npolar; ++j)
  {
    const unsigned int Na_j = num_azimuthal_per_level[j];

    // Create azimuthal angles for this polar level (Chebyshev distribution)
    for (unsigned int i = 0; i < Na_j; ++i)
      azimuthal_per_polar[j].emplace_back(M_PI * (2 * (i + 1) - 1) / Na_j);

    // Create weights: Gauss-Legendre weight (polar) * Gauss-Chebyshev weight (azimuthal)
    // Gauss-Chebyshev weight for integration over [0, 2π) is 2π/N
    const double gl_weight = gl_polar.weights[j];
    const double gc_weight = 2.0 * M_PI / Na_j;
    for (unsigned int i = 0; i < Na_j; ++i)
      wts_per_polar[j].emplace_back(gl_weight * gc_weight);
  }

  // Initialize
  AssembleTriangularCosines(azimuthal_per_polar, polar, wts_per_polar, verbose);
  MakeHarmonicIndices();
  BuildDiscreteToMomentOperator();
  BuildMomentToDiscreteOperator();
}

void
GLCTriangularQuadrature2DXY::AssembleTriangularCosines(
  const std::vector<std::vector<double>>& azimuthal_per_polar,
  const std::vector<double>& polar,
  const std::vector<std::vector<double>>& wts_per_polar,
  bool verbose)
{
  const size_t Np = polar.size();

  polar_ang_ = polar;
  azimuthal_per_polar_ = azimuthal_per_polar;

  if (verbose)
  {
    log.Log() << "Polar angles:";
    for (const auto& ang : polar_ang_)
      log.Log() << ang;
  }

  // Create angle pairs with varying azimuthal angles per polar level
  map_directions_.clear();
  for (unsigned int j = 0; j < Np; ++j)
    map_directions_.emplace(j, std::vector<unsigned int>());

  abscissae_.clear();
  weights_.clear();
  weight_sum_ = 0.0;

  unsigned int direction_index = 0;
  for (unsigned int j = 0; j < Np; ++j)
  {
    const size_t Na_j = azimuthal_per_polar[j].size();

    if (verbose)
    {
      log.Log() << "Polar level " << j << " (theta=" << polar[j] * 180.0 / M_PI << " deg): " << Na_j
                << " azimuthal angles";
    }

    for (unsigned int i = 0; i < Na_j; ++i)
    {
      map_directions_[j].emplace_back(direction_index);

      const auto abscissa = QuadraturePointPhiTheta(azimuthal_per_polar[j][i], polar[j]);
      abscissae_.emplace_back(abscissa);

      const double weight = wts_per_polar[j][i];
      weights_.emplace_back(weight);
      weight_sum_ += weight;

      ++direction_index;
    }
  }

  // Create omega list
  omegas_.clear();
  for (const auto& qpoint : abscissae_)
  {
    Vector3 new_omega;
    new_omega.x = std::sin(qpoint.theta) * std::cos(qpoint.phi);
    new_omega.y = std::sin(qpoint.theta) * std::sin(qpoint.phi);
    new_omega.z = std::cos(qpoint.theta);

    omegas_.emplace_back(new_omega);

    if (verbose)
      log.Log() << "Quadrature angle=" << new_omega.PrintStr();
  }

  // Normalize weights to 1.0
  for (auto& w : weights_)
    w /= weight_sum_;

  // Compute and store sum of weights after normalization
  weight_sum_ = 0.0;
  for (auto& w : weights_)
    weight_sum_ += w;
}

GLCTriangularQuadrature2DXY::GLCTriangularQuadrature2DXY(unsigned int Npolar,
                                                         unsigned int scattering_order,
                                                         bool verbose,
                                                         OperatorConstructionMethod method)
  : TriangularQuadrature(2, scattering_order, method)
{
  if (Npolar % 2 != 0)
    throw std::invalid_argument("GLCTriangularQuadrature2DXY: Npolar must be even.");

  // For 2D, we only use the upper hemisphere (half of the polar levels)
  const unsigned int half = Npolar / 2;

  // Compute maximum number of azimuthal angles at the equator
  // For the triangular reduction pattern with 4 fewer azimuthal angles per level from equator,
  // this ensures exactly 4 azimuthal angles at the pole (1 per quadrant)
  const unsigned int Nazimuthal = 2 * Npolar;

  n_polar_ = Npolar;
  n_azimuthal_ = Nazimuthal;

  GaussLegendreQuadrature gl_polar(Npolar);

  // Create polar angles (only upper hemisphere - first half of GL points)
  std::vector<double> polar;
  polar.reserve(half);
  for (unsigned int j = 0; j < half; ++j)
    polar.emplace_back(M_PI - acos(gl_polar.qpoints[j][0]));

  // Calculate number of azimuthal angles per polar level
  // At equator (j = half - 1): Nazimuthal angles
  // Each level toward the pole: 4 less azimuthal angles (1 less per quadrant)
  std::vector<unsigned int> num_azimuthal_per_level(half);
  for (unsigned int j = 0; j < half; ++j)
  {
    // Distance from equator: equator is at j = half - 1
    const unsigned int levels_from_equator = (half - 1) - j;
    // Reduce by 4 azimuthal angles per level (1 per quadrant) to maintain symmetry
    num_azimuthal_per_level[j] = Nazimuthal - 4 * levels_from_equator;
  }

  // Create azimuthal angles and weights per polar level
  std::vector<std::vector<double>> azimuthal_per_polar(half);
  std::vector<std::vector<double>> wts_per_polar(half);

  for (unsigned int j = 0; j < half; ++j)
  {
    const unsigned int Na_j = num_azimuthal_per_level[j];

    // Create azimuthal angles for this polar level (Chebyshev distribution)
    for (unsigned int i = 0; i < Na_j; ++i)
      azimuthal_per_polar[j].emplace_back(M_PI * (2 * (i + 1) - 1) / Na_j);

    // Create weights: Gauss-Legendre weight (polar) * Gauss-Chebyshev weight (azimuthal)
    // Gauss-Chebyshev weight for integration over [0, 2π) is 2π/N
    const double gl_weight = gl_polar.weights[j];
    const double gc_weight = 2.0 * M_PI / Na_j;
    for (unsigned int i = 0; i < Na_j; ++i)
      wts_per_polar[j].emplace_back(gl_weight * gc_weight);
  }

  // Initialize
  AssembleTriangularCosines(azimuthal_per_polar, polar, wts_per_polar, verbose);
  MakeHarmonicIndices();
  BuildDiscreteToMomentOperator();
  BuildMomentToDiscreteOperator();
}

} // namespace opensn

1	// SPDX-FileCopyrightText: 2026 The OpenSn Authors <https://open-sn.github.io/opensn/>
2	// SPDX-License-Identifier: MIT
3
4	#include "framework/math/quadratures/angular/triangular_quadrature.h"
5	#include "framework/math/quadratures/gausslegendre_quadrature.h"
6	#include "framework/logging/log.h"
7	#include "framework/runtime.h"
8	#include <cmath>
9
10	namespace opensn
11	{
12
13	void
14	GLCTriangularQuadrature3DXYZ::AssembleTriangularCosines(	1✔
15	const std::vector<std::vector<double>>& azimuthal_per_polar,
16	const std::vector<double>& polar,
17	const std::vector<std::vector<double>>& wts_per_polar,
18	bool verbose)
19	{
20	const size_t Np = polar.size();	1✔
21
22	polar_ang_ = polar;	1✔
23	azimuthal_per_polar_ = azimuthal_per_polar;	1✔
24
25	if (verbose)	1✔
26	{
27	log.Log() << "Polar angles:";	×
NEW 28	for (const auto& ang : polar_ang_)	×
29	log.Log() << ang;	×
30	}
31
32	// Create angle pairs with varying azimuthal angles per polar level
33	map_directions_.clear();	1✔
34	for (unsigned int j = 0; j < Np; ++j)	5✔
35	map_directions_.emplace(j, std::vector<unsigned int>());	8✔
36
37	abscissae_.clear();	1✔
38	weights_.clear();	1✔
39	weight_sum_ = 0.0;	1✔
40
41	unsigned int direction_index = 0;	1✔
42	for (unsigned int j = 0; j < Np; ++j)	5✔
43	{
44	const size_t Na_j = azimuthal_per_polar[j].size();	4✔
45
46	if (verbose)	4✔
47	{
48	log.Log() << "Polar level " << j << " (theta=" << polar[j] * 180.0 / M_PI << " deg): " << Na_j	×
49	<< " azimuthal angles";	×
50	}
51
52	for (unsigned int i = 0; i < Na_j; ++i)	28✔
53	{
54	map_directions_[j].emplace_back(direction_index);	24✔
55
56	const auto abscissa = QuadraturePointPhiTheta(azimuthal_per_polar[j][i], polar[j]);	24✔
57	abscissae_.emplace_back(abscissa);	24✔
58
59	const double weight = wts_per_polar[j][i];	24✔
60	weights_.emplace_back(weight);	24✔
61	weight_sum_ += weight;	24✔
62
63	++direction_index;	24✔
64	}
65	}
66
67	// Create omega list
68	omegas_.clear();	1✔
69	for (const auto& qpoint : abscissae_)	25✔
70	{
71	Vector3 new_omega;	24✔
72	new_omega.x = std::sin(qpoint.theta) * std::cos(qpoint.phi);	24✔
73	new_omega.y = std::sin(qpoint.theta) * std::sin(qpoint.phi);	24✔
74	new_omega.z = std::cos(qpoint.theta);	24✔
75
76	omegas_.emplace_back(new_omega);	24✔
77
78	if (verbose)	24✔
79	log.Log() << "Quadrature angle=" << new_omega.PrintStr();	×
80	}
81
82	// Normalize weights to 1.0
83	for (auto& w : weights_)	25✔
84	w /= weight_sum_;	24✔
85
86	// Compute and store sum of weights after normalization
87	weight_sum_ = 0.0;	1✔
88	for (auto& w : weights_)	25✔
89	weight_sum_ += w;	24✔
90	}	1✔
91
92	GLCTriangularQuadrature3DXYZ::GLCTriangularQuadrature3DXYZ(unsigned int Npolar,	1✔
93	unsigned int scattering_order,
94	bool verbose,
95	OperatorConstructionMethod method)	1✔
96	: TriangularQuadrature(3, scattering_order, method)	1✔
97	{
98	if (Npolar % 2 != 0)	1✔
99	throw std::invalid_argument("GLCTriangularQuadrature3DXYZ: Npolar must be even.");	×
100
101	// Compute maximum number of azimuthal angles at the equator
102	// For the triangular reduction pattern with 4 fewer azimuthal angles per level from equator,
103	// this ensures exactly 4 azimuthal angles at the poles (1 per octant)
104	const unsigned int Nazimuthal = 2 * Npolar;	1✔
105
106	n_polar_ = Npolar;	1✔
107	n_azimuthal_ = Nazimuthal;	1✔
108
109	GaussLegendreQuadrature gl_polar(Npolar);	1✔
110
111	// Create polar angles
112	std::vector<double> polar;	1✔
113	polar.reserve(Npolar);	1✔
114	for (unsigned int j = 0; j < Npolar; ++j)	5✔
115	polar.emplace_back(M_PI - std::acos(gl_polar.qpoints[j][0]));	4✔
116
117	// Calculate number of azimuthal angles per polar level
118	// At equator: Nazimuthal angles
119	// Each level away from equator: 4 less azimuthal angles (1 less per octant)
120	// This maintains octant symmetry in the triangular pattern
121	std::vector<unsigned int> num_azimuthal_per_level(Npolar);	1✔
122	for (unsigned int j = 0; j < Npolar; ++j)	5✔
123	{
124	// Distance from equator in terms of polar levels
125	// For even Npolar, the equatorial levels are at Npolar/2 - 1 and Npolar/2
126	// Levels from equator: how many steps from the nearest equatorial level
127	unsigned int levels_from_equator = 0;	4✔
128	if (j < Npolar / 2)	4✔
129	{
130	// Upper hemisphere: distance from equatorial level (Npolar/2 - 1)
131	levels_from_equator = (Npolar / 2 - 1) - j;	2✔
132	}
133	else
134	{
135	// Lower hemisphere: distance from equatorial level (Npolar/2)
136	levels_from_equator = j - (Npolar / 2);	2✔
137	}
138	// Reduce by 4 azimuthal angles per level (1 per octant) to maintain symmetry
139	num_azimuthal_per_level[j] = Nazimuthal - 4 * levels_from_equator;	4✔
140	}
141
142	// Create azimuthal angles and weights per polar level
143	std::vector<std::vector<double>> azimuthal_per_polar(Npolar);	1✔
144	std::vector<std::vector<double>> wts_per_polar(Npolar);	1✔
145
146	for (unsigned int j = 0; j < Npolar; ++j)	5✔
147	{
148	const unsigned int Na_j = num_azimuthal_per_level[j];	4✔
149
150	// Create azimuthal angles for this polar level (Chebyshev distribution)
151	for (unsigned int i = 0; i < Na_j; ++i)	28✔
152	azimuthal_per_polar[j].emplace_back(M_PI * (2 * (i + 1) - 1) / Na_j);	24✔
153
154	// Create weights: Gauss-Legendre weight (polar) * Gauss-Chebyshev weight (azimuthal)
155	// Gauss-Chebyshev weight for integration over [0, 2π) is 2π/N
156	const double gl_weight = gl_polar.weights[j];	4✔
157	const double gc_weight = 2.0 * M_PI / Na_j;	4✔
158	for (unsigned int i = 0; i < Na_j; ++i)	28✔
159	wts_per_polar[j].emplace_back(gl_weight * gc_weight);	24✔
160	}
161
162	// Initialize
163	AssembleTriangularCosines(azimuthal_per_polar, polar, wts_per_polar, verbose);	1✔
164	MakeHarmonicIndices();	1✔
165	BuildDiscreteToMomentOperator();	1✔
166	BuildMomentToDiscreteOperator();	1✔
167	}	2✔
168
169	void
170	GLCTriangularQuadrature2DXY::AssembleTriangularCosines(	2✔
171	const std::vector<std::vector<double>>& azimuthal_per_polar,
172	const std::vector<double>& polar,
173	const std::vector<std::vector<double>>& wts_per_polar,
174	bool verbose)
175	{
176	const size_t Np = polar.size();	2✔
177
178	polar_ang_ = polar;	2✔
179	azimuthal_per_polar_ = azimuthal_per_polar;	2✔
180
181	if (verbose)	2✔
182	{
183	log.Log() << "Polar angles:";	×
NEW 184	for (const auto& ang : polar_ang_)	×
UNCOV 185	log.Log() << ang;	×
186	}
187
188	// Create angle pairs with varying azimuthal angles per polar level
189	map_directions_.clear();	2✔
190	for (unsigned int j = 0; j < Np; ++j)	6✔
191	map_directions_.emplace(j, std::vector<unsigned int>());	8✔
192
193	abscissae_.clear();	2✔
194	weights_.clear();	2✔
195	weight_sum_ = 0.0;	2✔
196
197	unsigned int direction_index = 0;	2✔
198	for (unsigned int j = 0; j < Np; ++j)	6✔
199	{
200	const size_t Na_j = azimuthal_per_polar[j].size();	4✔
201
202	if (verbose)	4✔
203	{
204	log.Log() << "Polar level " << j << " (theta=" << polar[j] * 180.0 / M_PI << " deg): " << Na_j	×
UNCOV 205	<< " azimuthal angles";	×
206	}
207
208	for (unsigned int i = 0; i < Na_j; ++i)	28✔
209	{
210	map_directions_[j].emplace_back(direction_index);	24✔
211
212	const auto abscissa = QuadraturePointPhiTheta(azimuthal_per_polar[j][i], polar[j]);	24✔
213	abscissae_.emplace_back(abscissa);	24✔
214
215	const double weight = wts_per_polar[j][i];	24✔
216	weights_.emplace_back(weight);	24✔
217	weight_sum_ += weight;	24✔
218
219	++direction_index;	24✔
220	}
221	}
222
223	// Create omega list
224	omegas_.clear();	2✔
225	for (const auto& qpoint : abscissae_)	26✔
226	{
227	Vector3 new_omega;	24✔
228	new_omega.x = std::sin(qpoint.theta) * std::cos(qpoint.phi);	24✔
229	new_omega.y = std::sin(qpoint.theta) * std::sin(qpoint.phi);	24✔
230	new_omega.z = std::cos(qpoint.theta);	24✔
231
232	omegas_.emplace_back(new_omega);	24✔
233
234	if (verbose)	24✔
UNCOV 235	log.Log() << "Quadrature angle=" << new_omega.PrintStr();	×
236	}
237
238	// Normalize weights to 1.0
239	for (auto& w : weights_)	26✔
240	w /= weight_sum_;	24✔
241
242	// Compute and store sum of weights after normalization
243	weight_sum_ = 0.0;	2✔
244	for (auto& w : weights_)	26✔
245	weight_sum_ += w;	24✔
246	}	2✔
247
248	GLCTriangularQuadrature2DXY::GLCTriangularQuadrature2DXY(unsigned int Npolar,	2✔
249	unsigned int scattering_order,
250	bool verbose,
251	OperatorConstructionMethod method)	2✔
252	: TriangularQuadrature(2, scattering_order, method)	2✔
253	{
254	if (Npolar % 2 != 0)	2✔
UNCOV 255	throw std::invalid_argument("GLCTriangularQuadrature2DXY: Npolar must be even.");	×
256
257	// For 2D, we only use the upper hemisphere (half of the polar levels)
258	const unsigned int half = Npolar / 2;	2✔
259
260	// Compute maximum number of azimuthal angles at the equator
261	// For the triangular reduction pattern with 4 fewer azimuthal angles per level from equator,
262	// this ensures exactly 4 azimuthal angles at the pole (1 per quadrant)
263	const unsigned int Nazimuthal = 2 * Npolar;	2✔
264
265	n_polar_ = Npolar;	2✔
266	n_azimuthal_ = Nazimuthal;	2✔
267
268	GaussLegendreQuadrature gl_polar(Npolar);	2✔
269
270	// Create polar angles (only upper hemisphere - first half of GL points)
271	std::vector<double> polar;	2✔
272	polar.reserve(half);	2✔
273	for (unsigned int j = 0; j < half; ++j)	6✔
274	polar.emplace_back(M_PI - acos(gl_polar.qpoints[j][0]));	4✔
275
276	// Calculate number of azimuthal angles per polar level
277	// At equator (j = half - 1): Nazimuthal angles
278	// Each level toward the pole: 4 less azimuthal angles (1 less per quadrant)
279	std::vector<unsigned int> num_azimuthal_per_level(half);	2✔
280	for (unsigned int j = 0; j < half; ++j)	6✔
281	{
282	// Distance from equator: equator is at j = half - 1
283	const unsigned int levels_from_equator = (half - 1) - j;	4✔
284	// Reduce by 4 azimuthal angles per level (1 per quadrant) to maintain symmetry
285	num_azimuthal_per_level[j] = Nazimuthal - 4 * levels_from_equator;	4✔
286	}
287
288	// Create azimuthal angles and weights per polar level
289	std::vector<std::vector<double>> azimuthal_per_polar(half);	2✔
290	std::vector<std::vector<double>> wts_per_polar(half);	2✔
291
292	for (unsigned int j = 0; j < half; ++j)	6✔
293	{
294	const unsigned int Na_j = num_azimuthal_per_level[j];	4✔
295
296	// Create azimuthal angles for this polar level (Chebyshev distribution)
297	for (unsigned int i = 0; i < Na_j; ++i)	28✔
298	azimuthal_per_polar[j].emplace_back(M_PI * (2 * (i + 1) - 1) / Na_j);	24✔
299
300	// Create weights: Gauss-Legendre weight (polar) * Gauss-Chebyshev weight (azimuthal)
301	// Gauss-Chebyshev weight for integration over [0, 2π) is 2π/N
302	const double gl_weight = gl_polar.weights[j];	4✔
303	const double gc_weight = 2.0 * M_PI / Na_j;	4✔
304	for (unsigned int i = 0; i < Na_j; ++i)	28✔
305	wts_per_polar[j].emplace_back(gl_weight * gc_weight);	24✔
306	}
307
308	// Initialize
309	AssembleTriangularCosines(azimuthal_per_polar, polar, wts_per_polar, verbose);	2✔
310	MakeHarmonicIndices();	2✔
311	BuildDiscreteToMomentOperator();	2✔
312	BuildMomentToDiscreteOperator();	2✔
313	}	4✔
314
315	} // 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