Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "RayTracingAngularQuadrature.h"
11 :
12 : // MOOSE includes
13 : #include "MooseUtils.h"
14 :
15 : // libMesh includes
16 : #include "libmesh/dense_vector.h"
17 :
18 697 : RayTracingAngularQuadrature::RayTracingAngularQuadrature(const unsigned int dim,
19 : const unsigned int polar_order,
20 : const unsigned int azimuthal_order,
21 : const Real mu_min,
22 697 : const Real mu_max)
23 697 : : _dim(dim),
24 697 : _polar_order(polar_order),
25 697 : _azimuthal_order(azimuthal_order),
26 697 : _mu_min(mu_min),
27 697 : _mu_max(mu_max)
28 : {
29 697 : if (_polar_order == 0)
30 2 : mooseError("polar_order must be positive in RayTracingAngularQuadrature");
31 695 : if (_azimuthal_order == 0)
32 2 : mooseError("azimuthal_order must be positive in RayTracingAngularQuadrature");
33 693 : if (_mu_min >= _mu_max)
34 2 : mooseError("mu_min must be < mu_max in RayTracingAngularQuadrature");
35 691 : if (_mu_min < -1)
36 2 : mooseError("mu_min must be >= -1 in RayTracingAngularQuadrature");
37 689 : if (_mu_max > 1)
38 2 : mooseError("mu_max must be <= 1 in RayTracingAngularQuadrature");
39 687 : if (dim != 2 && dim != 3)
40 2 : mooseError("RayTracingAngularQuadrature only supports dimensions 2 and 3");
41 :
42 : // Build the quadrature
43 685 : build();
44 :
45 : // Default rotation is up
46 685 : rotate(libMesh::Point(0, 0, 1));
47 685 : }
48 :
49 : void
50 685 : RayTracingAngularQuadrature::build()
51 : {
52 : // Chebyshev quadrature on [0, 2\pi]
53 685 : std::vector<Real> chebyshev_x(_azimuthal_order);
54 685 : std::vector<Real> chebyshev_w(_azimuthal_order);
55 685 : chebyshev(_azimuthal_order, chebyshev_x, chebyshev_w);
56 :
57 : // Gauss-Legendre quadrature on [0, 1]
58 : std::vector<Real> gauss_legendre_x;
59 : std::vector<Real> gauss_legendre_w;
60 685 : gaussLegendre(_polar_order, gauss_legendre_x, gauss_legendre_w);
61 :
62 685 : _phi.resize(0);
63 685 : _mu.resize(0);
64 685 : _w.resize(0);
65 :
66 : // Build the product quadrature
67 2847 : for (std::size_t i = 0; i < chebyshev_x.size(); ++i)
68 10392 : for (std::size_t j = 0; j < gauss_legendre_x.size(); ++j)
69 : {
70 : // If 2D, throw away the downward mu
71 8230 : if (_dim == 2 && gauss_legendre_x[j] < 0.5 - TOLERANCE * TOLERANCE)
72 1449 : continue;
73 :
74 6781 : _phi.push_back(chebyshev_x[i]);
75 6781 : _mu.push_back(gauss_legendre_x[j] * (_mu_max - _mu_min) + _mu_min);
76 6781 : _polar_sin.push_back(std::sin(std::acos(_mu.back())));
77 :
78 : const Real weight_factor =
79 6781 : (_dim == 2 && !MooseUtils::absoluteFuzzyEqual(gauss_legendre_x[j], 0.5)) ? 2.0 : 1.0;
80 6781 : _w.push_back(weight_factor * chebyshev_w[i] * gauss_legendre_w[j]);
81 : }
82 685 : }
83 :
84 : void
85 688 : RayTracingAngularQuadrature::gaussLegendre(const unsigned int order,
86 : std::vector<Real> & x,
87 : std::vector<Real> & w)
88 : {
89 688 : if (order == 0)
90 2 : mooseError("Order must be positive in gaussLegendre()");
91 :
92 686 : x.resize(order);
93 686 : w.resize(order);
94 686 : libMesh::DenseMatrix<Real> mat(order, order);
95 686 : libMesh::DenseVector<Real> lambda(order);
96 686 : libMesh::DenseVector<Real> lambdai(order);
97 686 : libMesh::DenseMatrix<Real> vec(order, order);
98 :
99 3513 : for (unsigned int i = 1; i < order; ++i)
100 : {
101 2827 : Real ri = i;
102 2827 : mat(i, i - 1) = ri / std::sqrt(((2. * ri - 1.) * (2. * ri + 1.)));
103 2827 : mat(i - 1, i) = mat(i, i - 1);
104 : }
105 686 : mat.evd_right(lambda, lambdai, vec);
106 :
107 4199 : for (unsigned int i = 0; i < order; ++i)
108 : {
109 3513 : x[i] = 0.5 * (lambda(i) + 1.0);
110 3513 : w[i] = vec(0, i) * vec(0, i);
111 : }
112 :
113 : // Sort based on the points
114 686 : std::vector<std::size_t> sorted_indices(x.size());
115 : std::iota(sorted_indices.begin(), sorted_indices.end(), 0);
116 686 : std::stable_sort(sorted_indices.begin(),
117 : sorted_indices.end(),
118 5801 : [&x](size_t i1, size_t i2) { return x[i1] < x[i2]; });
119 686 : const auto x_copy = x;
120 686 : const auto w_copy = w;
121 4199 : for (std::size_t i = 0; i < x.size(); ++i)
122 : {
123 3513 : x[i] = x_copy[sorted_indices[i]];
124 3513 : w[i] = w_copy[sorted_indices[i]];
125 : }
126 1372 : }
127 :
128 : void
129 24414 : RayTracingAngularQuadrature::rotate(const libMesh::Point & rotation_direction)
130 : {
131 24414 : _current_rotation_direction = rotation_direction;
132 :
133 24414 : libMesh::DenseMatrix<Real> rotation_matrix;
134 24414 : rotationMatrix(rotation_direction.unit(), rotation_matrix);
135 :
136 24414 : _current_directions.clear();
137 24414 : _current_directions.reserve(_phi.size());
138 :
139 24414 : _current_weights.clear();
140 24414 : _current_weights.reserve(_w.size());
141 :
142 24414 : _current_polar_sins.clear();
143 24414 : _current_polar_sins.reserve(_polar_sin.size());
144 :
145 24414 : libMesh::DenseVector<Real> omega(3);
146 24414 : libMesh::DenseVector<Real> omega_p(3);
147 : Point direction;
148 :
149 262815 : for (std::size_t q = 0; q < _phi.size(); ++q)
150 : {
151 : // Get direction as a DenseVector for rotation
152 238401 : omega(0) = sqrt(1 - _mu[q] * _mu[q]) * cos(_phi[q]);
153 238401 : omega(1) = sqrt(1 - _mu[q] * _mu[q]) * sin(_phi[q]);
154 238401 : omega(2) = _mu[q];
155 :
156 : // Rotate and create the direction
157 238401 : rotation_matrix.vector_mult(omega_p, omega);
158 238401 : direction(0) = omega_p(0);
159 238401 : direction(1) = omega_p(1);
160 238401 : direction(2) = omega_p(2);
161 :
162 : // The index for the new direction in the "current" data
163 : // structures. By default - we add a new one. However,
164 : // in 2D we may have a duplicate projected direction
165 : // so we need to keep track of this in the case
166 : // that we are adding information to an already
167 : // existant direction
168 : std::size_t l = _current_directions.size();
169 :
170 : // If 2D, project to the xy plane and see if any other
171 : // projected directions are a duplicate
172 238401 : if (_dim == 2)
173 : {
174 12384 : direction(2) = 0;
175 12384 : direction /= direction.norm();
176 :
177 : // If we loop through all current directions and find
178 : // no duplicates, this will set l back to _current_directions.size()
179 47180 : for (l = 0; l < _current_directions.size(); ++l)
180 36034 : if (_current_directions[l].absolute_fuzzy_equals(direction))
181 : break;
182 : }
183 :
184 : // If l is the current size of _current_directions, it means
185 : // that no duplicates were found and that we are inserting
186 : // a new direction
187 238401 : if (l == _current_directions.size())
188 : {
189 237163 : _current_directions.push_back(direction);
190 237163 : _current_weights.resize(l + 1);
191 237163 : _current_polar_sins.resize(l + 1);
192 : }
193 : else
194 : mooseAssert(_dim == 2, "Should only have duplicates in 2D");
195 :
196 : // Add to weights and sins for this direction
197 238401 : _current_weights[l].push_back(_w[q]);
198 238401 : _current_polar_sins[l].push_back(_polar_sin[q]);
199 : }
200 24414 : }
201 :
202 : const libMesh::Point &
203 418611 : RayTracingAngularQuadrature::getDirection(const unsigned int l) const
204 : {
205 418611 : checkDirection(l);
206 418611 : return _current_directions[l];
207 : }
208 :
209 : const std::vector<Real> &
210 210 : RayTracingAngularQuadrature::getWeights(const unsigned int l) const
211 : {
212 210 : checkDirection(l);
213 210 : return _current_weights[l];
214 : }
215 :
216 : Real
217 1050 : RayTracingAngularQuadrature::getTotalWeight(const unsigned int l) const
218 : {
219 1050 : checkDirection(l);
220 1050 : return std::accumulate(_current_weights[l].begin(), _current_weights[l].end(), (Real)0);
221 : }
222 :
223 : const std::vector<Real> &
224 32 : RayTracingAngularQuadrature::getPolarSins(const unsigned int l) const
225 : {
226 32 : checkDirection(l);
227 32 : return _current_polar_sins[l];
228 : }
229 :
230 : std::size_t
231 20 : RayTracingAngularQuadrature::numPolar(const unsigned int l) const
232 : {
233 20 : checkDirection(l);
234 : if (_dim == 3)
235 : mooseAssert(_current_weights[l].size() == 1, "Should be 1 polar in 3D");
236 20 : return _current_polar_sins[l].size();
237 : }
238 :
239 : void
240 419925 : RayTracingAngularQuadrature::checkDirection(const unsigned int l) const
241 : {
242 419925 : if (!hasDirection(l))
243 2 : mooseError("RayTracingAngularQuadrature does not have direction ", l);
244 419923 : }
245 :
246 : void
247 686 : RayTracingAngularQuadrature::chebyshev(const unsigned int order,
248 : std::vector<Real> & x,
249 : std::vector<Real> & w)
250 : {
251 686 : x.resize(order);
252 686 : w.resize(order);
253 :
254 2948 : for (std::size_t i = 0; i < order; ++i)
255 : {
256 2262 : x[i] = 2 * (Real)i * M_PI / (Real)order;
257 2262 : w[i] = 2 * M_PI / (Real)order;
258 : }
259 686 : }
260 :
261 : void
262 24414 : RayTracingAngularQuadrature::rotationMatrix(const libMesh::Point & direction,
263 : libMesh::DenseMatrix<Real> & matrix)
264 : {
265 24414 : matrix.resize(3, 3);
266 : matrix.zero();
267 :
268 : // Create a local coordinate system around direction
269 24414 : const libMesh::Point tx = orthonormalVector(direction);
270 24414 : const libMesh::Point ty = direction.cross(tx);
271 :
272 : // Create the rotation matrix
273 97656 : for (unsigned int j = 0; j < 3; ++j)
274 : {
275 73242 : matrix(j, 0) = tx(j);
276 73242 : matrix(j, 1) = ty(j);
277 73242 : matrix(j, 2) = direction(j);
278 : }
279 24414 : }
280 :
281 : libMesh::Point
282 24416 : RayTracingAngularQuadrature::orthonormalVector(const libMesh::Point & v)
283 : {
284 24416 : if (MooseUtils::absoluteFuzzyLessEqual(v.norm(), 0))
285 2 : ::mooseError("Vector v has norm close to 0 in orthonormalVector()");
286 :
287 24414 : if (MooseUtils::absoluteFuzzyEqual(v(0), 0))
288 : return libMesh::Point(1, 0, 0);
289 7744 : if (MooseUtils::absoluteFuzzyEqual(v(1), 0))
290 : return libMesh::Point(0, 1, 0);
291 17 : if (MooseUtils::absoluteFuzzyEqual(v(2), 0))
292 : return libMesh::Point(0, 0, 1);
293 :
294 1 : return libMesh::Point(-v(1), v(0), 0).unit();
295 : }
|
