#include <RayTracingAngularQuadrature.h>

Public Member Functions
	RayTracingAngularQuadrature (const unsigned int dim, const unsigned int polar_order, const unsigned int azimuthal_order, const Real mu_min, const Real mu_max)
	Constructor.

void	rotate (const libMesh::Point &rotation_direction)
	Rotates the quadrature to a given direction.

const libMesh::Point &	getDirection (const unsigned int l) const
	Get the direction associated with direction `l`.

const std::vector< Real > &	getWeights (const unsigned int l) const
	Get the weights associated with the direction `l`.

Real	getTotalWeight (const unsigned int l) const
	Gets the total of the weights associated with the direction `l`.

const std::vector< Real > &	getPolarSins (const unsigned int l) const
	Gets the polar sins for the direction `l`.

std::size_t	numDirections () const
	Get the number of directions in the rotated and projected quadrature.

bool	hasDirection (const unsigned int l) const
	Whether or not the angular quadrature has direction `l`.

void	checkDirection (const unsigned int l) const
	Throws a MooseError if the angular quadrature does not have direction `l`.

std::size_t	numPolar (const unsigned int l) const
	The number of polar directions associated with the given direction.

unsigned int	polarOrder () const
	Get the polar order.

unsigned int	azimuthalOrder () const
	Get the azimuthal order.

Real	muMin () const
	Get the minimum mu.

Real	muMax () const
	Get the maximum mu.

const libMesh::Point &	currentRotationDirection () const
	Get the current rotation direction.

unsigned int	dim () const
	Get the quadrature dimension.

Static Public Member Functions
static void	gaussLegendre (const unsigned int order, std::vector< Real > &x, std::vector< Real > &w)
	Builds Gauss-Legendre quadrature on [0, 1] (symmetric about 0.5), with weights that sum to 1.

static void	chebyshev (const unsigned int order, std::vector< Real > &x, std::vector< Real > &w)
	Builds Chebyshev quadrature on [0, 2\pi] with weights that sum to 2\pi.

static void	rotationMatrix (const libMesh::Point &direction, libMesh::DenseMatrix< Real > &matrix)
	Builds the rotation matrix for direction `direction` into `matrix`.

static libMesh::Point	orthonormalVector (const libMesh::Point &v)
	Gets the vector that is orthonormal to `v`.

Private Member Functions
void	build ()
	Build the quadrature.

Private Attributes
const unsigned int	_dim
	The dimension.

const unsigned int	_polar_order
	The polar order.

const unsigned int	_azimuthal_order
	The azimuthal order.

const Real	_mu_min
	The minimum mu.

const Real	_mu_max
	The maximum mu.

std::vector< Real >	_phi
	Quadrature phi.

std::vector< Real >	_mu
	Quadrature mu.

std::vector< Real >	_polar_sin
	Quadrature polar sin.

std::vector< Real >	_w
	Quadrature combined weights.

libMesh::Point	_current_rotation_direction
	The current rotation direction.


std::vector< libMesh::Point >	_current_directions
	The current quadrature information.

std::vector< std::vector< Real > >	_current_weights

std::vector< std::vector< Real > >	_current_polar_sins

Detailed Description

Definition at line 19 of file RayTracingAngularQuadrature.h.

Constructor & Destructor Documentation

◆ RayTracingAngularQuadrature()

RayTracingAngularQuadrature::RayTracingAngularQuadrature	(	const unsigned int	dim,
		const unsigned int	polar_order,
		const unsigned int	azimuthal_order,
		const Real	mu_min,
		const Real	mu_max
	)

Constructor.

Parameters

dim	The mesh dimension
polar_order	Order for the polar quadrature
azimuthal_order	Azimuthal order for the quadrature
mu_min	Minimum mu for the quadrature
mu_max	Maximum mu for the quadrature

Definition at line 18 of file RayTracingAngularQuadrature.C.

  : _dim(dim),
    _polar_order(polar_order),
    _azimuthal_order(azimuthal_order),
    _mu_min(mu_min),
    _mu_max(mu_max)
{
  if (_polar_order == 0)
    mooseError("polar_order must be positive in RayTracingAngularQuadrature");
  if (_azimuthal_order == 0)
    mooseError("azimuthal_order must be positive in RayTracingAngularQuadrature");
  if (_mu_min >= _mu_max)
    mooseError("mu_min must be < mu_max in RayTracingAngularQuadrature");
  if (_mu_min < -1)
    mooseError("mu_min must be >= -1 in RayTracingAngularQuadrature");
  if (_mu_max > 1)
    mooseError("mu_max must be <= 1 in RayTracingAngularQuadrature");
  if (dim != 2 && dim != 3)
    mooseError("RayTracingAngularQuadrature only supports dimensions 2 and 3");
 
  // Build the quadrature
  build();
 
  // Default rotation is up
  rotate(libMesh::Point(0, 0, 1));
}

Member Function Documentation

◆ azimuthalOrder()

unsigned int RayTracingAngularQuadrature::azimuthalOrder ( ) const

inline

Get the azimuthal order.

Definition at line 146 of file RayTracingAngularQuadrature.h.

146{ return _azimuthal_order; }

Referenced by TEST().

◆ build()

void RayTracingAngularQuadrature::build ( )

private

Build the quadrature.

Definition at line 50 of file RayTracingAngularQuadrature.C.

{
  // Chebyshev quadrature on [0, 2\pi]
  std::vector<Real> chebyshev_x(_azimuthal_order);
  std::vector<Real> chebyshev_w(_azimuthal_order);
  chebyshev(_azimuthal_order, chebyshev_x, chebyshev_w);
 
  // Gauss-Legendre quadrature on [0, 1]
  std::vector<Real> gauss_legendre_x;
  std::vector<Real> gauss_legendre_w;
  gaussLegendre(_polar_order, gauss_legendre_x, gauss_legendre_w);
 
  _phi.resize(0);
  _mu.resize(0);
  _w.resize(0);
 
  // Build the product quadrature
  for (std::size_t i = 0; i < chebyshev_x.size(); ++i)
    for (std::size_t j = 0; j < gauss_legendre_x.size(); ++j)
    {
      // If 2D, throw away the downward mu
      if (_dim == 2 && gauss_legendre_x[j] < 0.5 - TOLERANCE * TOLERANCE)
        continue;
 
      _phi.push_back(chebyshev_x[i]);
      _mu.push_back(gauss_legendre_x[j] * (_mu_max - _mu_min) + _mu_min);
      _polar_sin.push_back(std::sin(std::acos(_mu.back())));
 
      const Real weight_factor =
          (_dim == 2 && !MooseUtils::absoluteFuzzyEqual(gauss_legendre_x[j], 0.5)) ? 2.0 : 1.0;
      _w.push_back(weight_factor * chebyshev_w[i] * gauss_legendre_w[j]);
    }
}

Referenced by RayTracingAngularQuadrature().

◆ chebyshev()

void RayTracingAngularQuadrature::chebyshev	(	const unsigned int	order,
		std::vector< Real > &	x,
		std::vector< Real > &	w
	)

static

Builds Chebyshev quadrature on [0, 2\pi] with weights that sum to 2\pi.

Parameters

order	The quadrature order
x	Vector to fill the points into
w	Vector to fill the weights into

Definition at line 247 of file RayTracingAngularQuadrature.C.

{
  x.resize(order);
  w.resize(order);
 
  for (std::size_t i = 0; i < order; ++i)
  {
    x[i] = 2 * (Real)i * M_PI / (Real)order;
    w[i] = 2 * M_PI / (Real)order;
  }
}

Referenced by build(), and TEST().

◆ checkDirection()

void RayTracingAngularQuadrature::checkDirection ( const unsigned int l ) const

Throws a MooseError if the angular quadrature does not have direction l.

Definition at line 240 of file RayTracingAngularQuadrature.C.

{
  if (!hasDirection(l))
    mooseError("RayTracingAngularQuadrature does not have direction ", l);
}

Referenced by getDirection(), getPolarSins(), getTotalWeight(), getWeights(), numPolar(), and RayTracingAngularQuadratureErrorTest::RayTracingAngularQuadratureErrorTest().

◆ currentRotationDirection()

const libMesh::Point & RayTracingAngularQuadrature::currentRotationDirection ( ) const

inline

Get the current rotation direction.

Definition at line 159 of file RayTracingAngularQuadrature.h.

159{ return _current_rotation_direction; }

RayTracingAngularQuadrature::_current_rotation_direction

libMesh::Point _current_rotation_direction

The current rotation direction.

Definition RayTracingAngularQuadrature.h:208

Referenced by TEST().

◆ dim()

unsigned int RayTracingAngularQuadrature::dim ( ) const

inline

Get the quadrature dimension.

Definition at line 164 of file RayTracingAngularQuadrature.h.

164{ return _dim; }

Referenced by RayTracingAngularQuadrature(), and testPolarSins2D().

◆ gaussLegendre()

void RayTracingAngularQuadrature::gaussLegendre	(	const unsigned int	order,
		std::vector< Real > &	x,
		std::vector< Real > &	w
	)

static

Builds Gauss-Legendre quadrature on [0, 1] (symmetric about 0.5), with weights that sum to 1.

Parameters

order	The quadrature order
x	Vector to fill the points into
w	Vector to fill the weights into

Definition at line 85 of file RayTracingAngularQuadrature.C.

{
  if (order == 0)
    mooseError("Order must be positive in gaussLegendre()");
 
  x.resize(order);
  w.resize(order);
  libMesh::DenseMatrix<Real> mat(order, order);
  libMesh::DenseVector<Real> lambda(order);
  libMesh::DenseVector<Real> lambdai(order);
  libMesh::DenseMatrix<Real> vec(order, order);
 
  for (unsigned int i = 1; i < order; ++i)
  {
    Real ri = i;
    mat(i, i - 1) = ri / std::sqrt(((2. * ri - 1.) * (2. * ri + 1.)));
    mat(i - 1, i) = mat(i, i - 1);
  }
  mat.evd_right(lambda, lambdai, vec);
 
  for (unsigned int i = 0; i < order; ++i)
  {
    x[i] = 0.5 * (lambda(i) + 1.0);
    w[i] = vec(0, i) * vec(0, i);
  }
 
  // Sort based on the points
  std::vector<std::size_t> sorted_indices(x.size());
  std::iota(sorted_indices.begin(), sorted_indices.end(), 0);
  std::stable_sort(sorted_indices.begin(),
                   sorted_indices.end(),
                   [&x](size_t i1, size_t i2) { return x[i1] < x[i2]; });
  const auto x_copy = x;
  const auto w_copy = w;
  for (std::size_t i = 0; i < x.size(); ++i)
  {
    x[i] = x_copy[sorted_indices[i]];
    w[i] = w_copy[sorted_indices[i]];
  }
}

Referenced by build(), RayTracingAngularQuadratureErrorTest::RayTracingAngularQuadratureErrorTest(), TEST(), and ViewFactorRayStudy::ViewFactorRayStudy().

◆ getDirection()

const libMesh::Point & RayTracingAngularQuadrature::getDirection ( const unsigned int l ) const

Get the direction associated with direction l.

This direction will be rotated per currentRotationDirection(), set by setRotation().

Definition at line 203 of file RayTracingAngularQuadrature.C.

{
  checkDirection(l);
  return _current_directions[l];
}

Referenced by ConeRayStudy::defineRays(), LotsOfRaysRayStudy::defineRays(), and testDirections().

◆ getPolarSins()

const std::vector< Real > & RayTracingAngularQuadrature::getPolarSins ( const unsigned int l ) const

Gets the polar sins for the direction l.

In the case of 2D, the 3D directions projected into the 2D plane may overlap. In this case, we combine the directions into a single direction with multiple polar sins. This is the reason for the vector return.

Definition at line 224 of file RayTracingAngularQuadrature.C.

{
  checkDirection(l);
  return _current_polar_sins[l];
}

Referenced by testPolarSins2D().

◆ getTotalWeight()

Real RayTracingAngularQuadrature::getTotalWeight ( const unsigned int l ) const

Gets the total of the weights associated with the direction l.

The weights across all directions sum to 2\pi.

In the case of 2D, the 3D directions projected into the 2D plane may overlap. In this case, we combine the directions into a single direction with multiple weights. This is the reason why a getter for the "total" weight for a single direction exists.

Definition at line 217 of file RayTracingAngularQuadrature.C.

{
  checkDirection(l);
  return std::accumulate(_current_weights[l].begin(), _current_weights[l].end(), (Real)0);
}

Referenced by ConeRayStudy::defineRays().

◆ getWeights()

const std::vector< Real > & RayTracingAngularQuadrature::getWeights ( const unsigned int l ) const

Get the weights associated with the direction l.

The weights across all directions sum to 2\pi.

In the case of 2D, the 3D directions projected into the 2D plane may overlap. In this case, we combine the directions into a single direction with multiple weights. This is the reason for the vector return.

Definition at line 210 of file RayTracingAngularQuadrature.C.

{
  checkDirection(l);
  return _current_weights[l];
}

Referenced by testWeights().

◆ hasDirection()

bool RayTracingAngularQuadrature::hasDirection ( const unsigned int l ) const

inline

Whether or not the angular quadrature has direction l.

Definition at line 121 of file RayTracingAngularQuadrature.h.

121{ return l < _current_directions.size(); }

Referenced by checkDirection(), and TEST().

◆ muMax()

Real RayTracingAngularQuadrature::muMax ( ) const

inline

Get the maximum mu.

Definition at line 154 of file RayTracingAngularQuadrature.h.

154{ return _mu_max; }

Referenced by TEST().

◆ muMin()

Real RayTracingAngularQuadrature::muMin ( ) const

inline

Get the minimum mu.

Definition at line 150 of file RayTracingAngularQuadrature.h.

150{ return _mu_min; }

Referenced by TEST().

◆ numDirections()

std::size_t RayTracingAngularQuadrature::numDirections ( ) const

inline

Get the number of directions in the rotated and projected quadrature.

Note that in the case of 2D, we overlap 3D directions projected into the 2D plane. Therefore, in 2D, this value may be less than the number of directions in the 3D quadrature.

Definition at line 117 of file RayTracingAngularQuadrature.h.

117{ return _current_directions.size(); }

Referenced by ConeRayStudy::defineRays(), LotsOfRaysRayStudy::defineRays(), TEST(), TEST(), testDirections(), testPolarSins2D(), and testWeights().

◆ numPolar()

std::size_t RayTracingAngularQuadrature::numPolar ( const unsigned int l ) const

The number of polar directions associated with the given direction.

In 2D, projecting the 3D directions into the 2D plane may result in directions that overlap, in which case we end up with one direction and multiple polars for said direction.

In 3D, this will always be 1.

Definition at line 231 of file RayTracingAngularQuadrature.C.

{
  checkDirection(l);
  if (_dim == 3)
    mooseAssert(_current_weights[l].size() == 1, "Should be 1 polar in 3D");
  return _current_polar_sins[l].size();
}

Referenced by testPolarSins2D().

◆ orthonormalVector()

libMesh::Point RayTracingAngularQuadrature::orthonormalVector ( const libMesh::Point & v )

static

Gets the vector that is orthonormal to v.

Definition at line 282 of file RayTracingAngularQuadrature.C.

{
  if (MooseUtils::absoluteFuzzyLessEqual(v.norm(), 0))
    ::mooseError("Vector v has norm close to 0 in orthonormalVector()");
 
  if (MooseUtils::absoluteFuzzyEqual(v(0), 0))
    return libMesh::Point(1, 0, 0);
  if (MooseUtils::absoluteFuzzyEqual(v(1), 0))
    return libMesh::Point(0, 1, 0);
  if (MooseUtils::absoluteFuzzyEqual(v(2), 0))
    return libMesh::Point(0, 0, 1);
 
  return libMesh::Point(-v(1), v(0), 0).unit();
}

Referenced by RayTracingAngularQuadratureErrorTest::RayTracingAngularQuadratureErrorTest(), and rotationMatrix().

◆ polarOrder()

unsigned int RayTracingAngularQuadrature::polarOrder ( ) const

inline

Get the polar order.

Definition at line 142 of file RayTracingAngularQuadrature.h.

142{ return _polar_order; }

Referenced by TEST().

◆ rotate()

void RayTracingAngularQuadrature::rotate ( const libMesh::Point & rotation_direction )

Rotates the quadrature to a given direction.

Definition at line 129 of file RayTracingAngularQuadrature.C.

{
  _current_rotation_direction = rotation_direction;
 
  libMesh::DenseMatrix<Real> rotation_matrix;
  rotationMatrix(rotation_direction.unit(), rotation_matrix);
 
  _current_directions.clear();
  _current_directions.reserve(_phi.size());
 
  _current_weights.clear();
  _current_weights.reserve(_w.size());
 
  _current_polar_sins.clear();
  _current_polar_sins.reserve(_polar_sin.size());
 
  libMesh::DenseVector<Real> omega(3);
  libMesh::DenseVector<Real> omega_p(3);
  Point direction;
 
  for (std::size_t q = 0; q < _phi.size(); ++q)
  {
    // Get direction as a DenseVector for rotation
    omega(0) = sqrt(1 - _mu[q] * _mu[q]) * cos(_phi[q]);
    omega(1) = sqrt(1 - _mu[q] * _mu[q]) * sin(_phi[q]);
    omega(2) = _mu[q];
 
    // Rotate and create the direction
    rotation_matrix.vector_mult(omega_p, omega);
    direction(0) = omega_p(0);
    direction(1) = omega_p(1);
    direction(2) = omega_p(2);
 
    // The index for the new direction in the "current" data
    // structures. By default - we add a new one. However,
    // in 2D we may have a duplicate projected direction
    // so we need to keep track of this in the case
    // that we are adding information to an already
    // existant direction
    std::size_t l = _current_directions.size();
 
    // If 2D, project to the xy plane and see if any other
    // projected directions are a duplicate
    if (_dim == 2)
    {
      direction(2) = 0;
      direction /= direction.norm();
 
      // If we loop through all current directions and find
      // no duplicates, this will set l back to _current_directions.size()
      for (l = 0; l < _current_directions.size(); ++l)
        if (_current_directions[l].absolute_fuzzy_equals(direction))
          break;
    }
 
    // If l is the current size of _current_directions, it means
    // that no duplicates were found and that we are inserting
    // a new direction
    if (l == _current_directions.size())
    {
      _current_directions.push_back(direction);
      _current_weights.resize(l + 1);
      _current_polar_sins.resize(l + 1);
    }
    else
      mooseAssert(_dim == 2, "Should only have duplicates in 2D");
 
    // Add to weights and sins for this direction
    _current_weights[l].push_back(_w[q]);
    _current_polar_sins[l].push_back(_polar_sin[q]);
  }
}

Referenced by ConeRayStudy::defineRays(), RayTracingAngularQuadrature(), TEST(), TEST(), TEST(), TEST(), and TEST().

◆ rotationMatrix()

void RayTracingAngularQuadrature::rotationMatrix	(	const libMesh::Point &	direction,
		libMesh::DenseMatrix< Real > &	matrix
	)

static

Builds the rotation matrix for direction direction into matrix.

Definition at line 262 of file RayTracingAngularQuadrature.C.

{
  matrix.resize(3, 3);
  matrix.zero();
 
  // Create a local coordinate system around direction
  const libMesh::Point tx = orthonormalVector(direction);
  const libMesh::Point ty = direction.cross(tx);
 
  // Create the rotation matrix
  for (unsigned int j = 0; j < 3; ++j)
  {
    matrix(j, 0) = tx(j);
    matrix(j, 1) = ty(j);
    matrix(j, 2) = direction(j);
  }
}

Referenced by rotate().

Member Data Documentation

◆ _azimuthal_order

const unsigned int RayTracingAngularQuadrature::_azimuthal_order

private

The azimuthal order.

Definition at line 177 of file RayTracingAngularQuadrature.h.

Referenced by azimuthalOrder(), build(), and RayTracingAngularQuadrature().

◆ _current_directions

std::vector<libMesh::Point> RayTracingAngularQuadrature::_current_directions

private

The current quadrature information.

We need "current" information because the quadrature can be rotated and because in 2D, there are cases where projected directions that are the same in the 2D plane are combined together with multiple weights and polar sins.

Definition at line 202 of file RayTracingAngularQuadrature.h.

Referenced by getDirection(), hasDirection(), numDirections(), and rotate().