RayTracingAngularQuadrature Class Reference

#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. More...
void	rotate (const libMesh::Point &rotation_direction)
	Rotates the quadrature to a given direction. More...
const libMesh::Point &	getDirection (const unsigned int l) const
	Get the direction associated with direction l. More...
const std::vector< Real > &	getWeights (const unsigned int l) const
	Get the weights associated with the direction l. More...
Real	getTotalWeight (const unsigned int l) const
	Gets the total of the weights associated with the direction l. More...
const std::vector< Real > &	getPolarSins (const unsigned int l) const
	Gets the polar sins for the direction l. More...
std::size_t	numDirections () const
	Get the number of directions in the rotated and projected quadrature. More...
bool	hasDirection (const unsigned int l) const
	Whether or not the angular quadrature has direction l. More...
void	checkDirection (const unsigned int l) const
	Throws a MooseError if the angular quadrature does not have direction l. More...
std::size_t	numPolar (const unsigned int l) const
	The number of polar directions associated with the given direction. More...
unsigned int	polarOrder () const
	Get the polar order. More...
unsigned int	azimuthalOrder () const
	Get the azimuthal order. More...
Real	muMin () const
	Get the minimum mu. More...
Real	muMax () const
	Get the maximum mu. More...
const libMesh::Point &	currentRotationDirection () const
	Get the current rotation direction. More...
unsigned int	dim () const
	Get the quadrature dimension. More...

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, with weights that sum to 1. More...
static void	chebyshev (const unsigned int order, std::vector< Real > &x, std::vector< Real > &w)
	Builds Chebyshev quadrature on [0, 2] with weights that sum to 2. More...
static void	rotationMatrix (const libMesh::Point &direction, libMesh::DenseMatrix< Real > &matrix)
	Builds the rotation matrix for direction direction into matrix. More...
static libMesh::Point	orthonormalVector (const libMesh::Point &v)
	Gets the vector that is orthonormal to v. More...

Private Member Functions
void	build ()
	Build the quadrature. More...

Private Attributes
const unsigned int	_dim
	The dimension. More...
const unsigned int	_polar_order
	The polar order. More...
const unsigned int	_azimuthal_order
	The azimuthal order. More...
const Real	_mu_min
	The minimum mu. More...
const Real	_mu_max
	The maximum mu. More...
std::vector< Real >	_phi
	Quadrature phi. More...
std::vector< Real >	_mu
	Quadrature mu. More...
std::vector< Real >	_polar_sin
	Quadrature polar sin. More...
std::vector< Real >	_w
	Quadrature combined weights. More...
libMesh::Point	_current_rotation_direction
	The current rotation direction. More...
std::vector< libMesh::Point >	_current_directions
	The current quadrature information. More...
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.

Referenced by TEST().

{ return _azimuthal_order; }

build()

void RayTracingAngularQuadrature::build ( )

private

Build the quadrature.

Definition at line 50 of file RayTracingAngularQuadrature.C.

Referenced by RayTracingAngularQuadrature().

{
  // 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]);
    }
}

chebyshev()

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

static

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

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.

Referenced by build(), and TEST().

{
  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;
  }
}

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.

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

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

currentRotationDirection()

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

inline

Get the current rotation direction.

Definition at line 159 of file RayTracingAngularQuadrature.h.

Referenced by TEST().

{ return _current_rotation_direction; }

dim()

unsigned int RayTracingAngularQuadrature::dim ( ) const

inline

Get the quadrature dimension.

Definition at line 164 of file RayTracingAngularQuadrature.h.

Referenced by RayTracingAngularQuadrature(), and testPolarSins2D().

{ return _dim; }

gaussLegendre()

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

static

Builds Gauss-Legendre quadrature on 0, 1, 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.

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

{
  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]];
  }
}

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.

Referenced by testDirections().

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

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.

Referenced by testPolarSins2D().

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

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.

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);
}

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.

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.

Referenced by testWeights().

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

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.

Referenced by checkDirection(), and TEST().

{ return l < _current_directions.size(); }

muMax()

Real RayTracingAngularQuadrature::muMax ( ) const

inline

Get the maximum mu.

Definition at line 154 of file RayTracingAngularQuadrature.h.

Referenced by TEST().

{ return _mu_max; }

muMin()

Real RayTracingAngularQuadrature::muMin ( ) const

inline

Get the minimum mu.

Definition at line 150 of file RayTracingAngularQuadrature.h.

Referenced by TEST().

{ return _mu_min; }

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.

Referenced by TEST(), testDirections(), testPolarSins2D(), and testWeights().

{ return _current_directions.size(); }

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.

Referenced by testPolarSins2D().

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

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.

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

{
  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();
}

polarOrder()

unsigned int RayTracingAngularQuadrature::polarOrder ( ) const

inline

Get the polar order.

Definition at line 142 of file RayTracingAngularQuadrature.h.

Referenced by TEST().

{ return _polar_order; }

rotate()

void RayTracingAngularQuadrature::rotate

(

const libMesh::Point &

rotation_direction

)

Rotates the quadrature to a given direction.

Definition at line 129 of file RayTracingAngularQuadrature.C.

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

{
  _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]);
  }
}

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.

Referenced by rotate().

{
  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);
  }
}

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().

_current_polar_sins

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

private

Definition at line 204 of file RayTracingAngularQuadrature.h.

Referenced by getPolarSins(), numPolar(), and rotate().

_current_rotation_direction

libMesh::Point RayTracingAngularQuadrature::_current_rotation_direction

private

The current rotation direction.

Definition at line 208 of file RayTracingAngularQuadrature.h.

Referenced by currentRotationDirection(), and rotate().

_current_weights

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

private

Definition at line 203 of file RayTracingAngularQuadrature.h.

Referenced by getTotalWeight(), getWeights(), numPolar(), and rotate().

_dim

const unsigned int RayTracingAngularQuadrature::_dim

private

The dimension.

Definition at line 173 of file RayTracingAngularQuadrature.h.

Referenced by build(), dim(), numPolar(), and rotate().

_mu

std::vector<Real> RayTracingAngularQuadrature::_mu

private

Quadrature mu.

Definition at line 186 of file RayTracingAngularQuadrature.h.

Referenced by build(), and rotate().

_mu_max

const Real RayTracingAngularQuadrature::_mu_max

private

The maximum mu.

Definition at line 181 of file RayTracingAngularQuadrature.h.

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

_mu_min

const Real RayTracingAngularQuadrature::_mu_min

private

The minimum mu.

Definition at line 179 of file RayTracingAngularQuadrature.h.

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

_phi

std::vector<Real> RayTracingAngularQuadrature::_phi

private

Quadrature phi.

Definition at line 184 of file RayTracingAngularQuadrature.h.

Referenced by build(), and rotate().

_polar_order

const unsigned int RayTracingAngularQuadrature::_polar_order

private

The polar order.

Definition at line 175 of file RayTracingAngularQuadrature.h.

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

_polar_sin

std::vector<Real> RayTracingAngularQuadrature::_polar_sin

private

Quadrature polar sin.

Definition at line 188 of file RayTracingAngularQuadrature.h.

Referenced by build(), and rotate().

_w

std::vector<Real> RayTracingAngularQuadrature::_w

private

Quadrature combined weights.

Definition at line 190 of file RayTracingAngularQuadrature.h.

Referenced by build(), and rotate().

---

The documentation for this class was generated from the following files:

• ray_tracing/include/utils/RayTracingAngularQuadrature.h
• ray_tracing/src/utils/RayTracingAngularQuadrature.C