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             : #pragma once
      11             : 
      12             : // MOOSE includes
      13             : #include "MooseTypes.h"
      14             : 
      15             : // libMesh includes
      16             : #include "libmesh/dense_matrix.h"
      17             : #include "libmesh/point.h"
      18             : 
      19             : class RayTracingAngularQuadrature
      20             : {
      21             : public:
      22             :   /**
      23             :    * Constructor.
      24             :    *
      25             :    * @param dim The mesh dimension
      26             :    * @param polar_order Order for the polar quadrature
      27             :    * @param azimuthal_order Azimuthal order for the quadrature
      28             :    * @param mu_min Minimum mu for the quadrature
      29             :    * @param mu_max Maximum mu for the quadrature
      30             :    */
      31             :   RayTracingAngularQuadrature(const unsigned int dim,
      32             :                               const unsigned int polar_order,
      33             :                               const unsigned int azimuthal_order,
      34             :                               const Real mu_min,
      35             :                               const Real mu_max);
      36             : 
      37             :   /**
      38             :    * Builds Gauss-Legendre quadrature on [0, 1] (symmetric about 0.5), with
      39             :    * weights that sum to 1.
      40             :    *
      41             :    * @param order The quadrature order
      42             :    * @param x Vector to fill the points into
      43             :    * @param w Vector to fill the weights into
      44             :    */
      45             :   static void gaussLegendre(const unsigned int order, std::vector<Real> & x, std::vector<Real> & w);
      46             : 
      47             :   /**
      48             :    * Builds Chebyshev quadrature on [0, 2\pi] with weights that sum to 2\pi.
      49             :    *
      50             :    * @param order The quadrature order
      51             :    * @param x Vector to fill the points into
      52             :    * @param w Vector to fill the weights into
      53             :    */
      54             :   static void chebyshev(const unsigned int order, std::vector<Real> & x, std::vector<Real> & w);
      55             : 
      56             :   /**
      57             :    * Builds the rotation matrix for direction \p direction into \p matrix.
      58             :    */
      59             :   static void rotationMatrix(const libMesh::Point & direction, libMesh::DenseMatrix<Real> & matrix);
      60             :   /**
      61             :    * Gets the vector that is orthonormal to \p v.
      62             :    */
      63             :   static libMesh::Point orthonormalVector(const libMesh::Point & v);
      64             : 
      65             :   /**
      66             :    * Rotates the quadrature to a given direction.
      67             :    */
      68             :   void rotate(const libMesh::Point & rotation_direction);
      69             : 
      70             :   /**
      71             :    * Get the direction associated with direction \p l.
      72             :    *
      73             :    * This direction will be rotated per currentRotationDirection(),
      74             :    * set by setRotation().
      75             :    */
      76             :   const libMesh::Point & getDirection(const unsigned int l) const;
      77             :   /**
      78             :    * Get the weights associated with the direction \p l.
      79             :    *
      80             :    * The weights across all directions sum to 2\pi.
      81             :    *
      82             :    * In the case of 2D, the 3D directions projected into the
      83             :    * 2D plane may overlap. In this case, we combine the directions
      84             :    * into a single direction with multiple weights. This is the
      85             :    * reason for the vector return.
      86             :    */
      87             :   const std::vector<Real> & getWeights(const unsigned int l) const;
      88             :   /**
      89             :    * Gets the total of the weights associated with the direction \p l.
      90             :    *
      91             :    * The weights across all directions sum to 2\pi.
      92             :    *
      93             :    * In the case of 2D, the 3D directions projected into the
      94             :    * 2D plane may overlap. In this case, we combine the directions
      95             :    * into a single direction with multiple weights. This is the
      96             :    * reason why a getter for the "total" weight for a single
      97             :    * direction exists.
      98             :    */
      99             :   Real getTotalWeight(const unsigned int l) const;
     100             :   /**
     101             :    * Gets the polar sins for the direction \p l.
     102             :    *
     103             :    * In the case of 2D, the 3D directions projected into the
     104             :    * 2D plane may overlap. In this case, we combine the directions
     105             :    * into a single direction with multiple polar sins. This is the
     106             :    * reason for the vector return.
     107             :    */
     108             :   const std::vector<Real> & getPolarSins(const unsigned int l) const;
     109             : 
     110             :   /**
     111             :    * Get the number of directions in the rotated and projected quadrature.
     112             :    *
     113             :    * Note that in the case of 2D, we overlap 3D directions projected
     114             :    * into the 2D plane. Therefore, in 2D, this value may be less
     115             :    * than the number of directions in the 3D quadrature.
     116             :    */
     117             :   std::size_t numDirections() const { return _current_directions.size(); }
     118             :   /**
     119             :    * Whether or not the angular quadrature has direction \p l.
     120             :    */
     121      419925 :   bool hasDirection(const unsigned int l) const { return l < _current_directions.size(); }
     122             : 
     123             :   /**
     124             :    * Throws a MooseError if the angular quadrature does not have direction \p l.
     125             :    */
     126             :   void checkDirection(const unsigned int l) const;
     127             : 
     128             :   /**
     129             :    * The number of polar directions associated with the given direction.
     130             :    *
     131             :    * In 2D, projecting the 3D directions into the 2D plane may result
     132             :    * in directions that overlap, in which case we end up with one
     133             :    * direction and multiple polars for said direction.
     134             :    *
     135             :    * In 3D, this will always be 1.
     136             :    */
     137             :   std::size_t numPolar(const unsigned int l) const;
     138             : 
     139             :   /**
     140             :    * Get the polar order
     141             :    */
     142           1 :   unsigned int polarOrder() const { return _polar_order; }
     143             :   /**
     144             :    * Get the azimuthal order
     145             :    */
     146           1 :   unsigned int azimuthalOrder() const { return _azimuthal_order; }
     147             :   /**
     148             :    * Get the minimum mu
     149             :    */
     150           1 :   Real muMin() const { return _mu_min; }
     151             :   /**
     152             :    * Get the maximum mu
     153             :    */
     154           1 :   Real muMax() const { return _mu_max; }
     155             : 
     156             :   /**
     157             :    * Get the current rotation direction
     158             :    */
     159             :   const libMesh::Point & currentRotationDirection() const { return _current_rotation_direction; }
     160             : 
     161             :   /**
     162             :    * Get the quadrature dimension
     163             :    */
     164           5 :   unsigned int dim() const { return _dim; }
     165             : 
     166             : private:
     167             :   /**
     168             :    * Build the quadrature
     169             :    */
     170             :   void build();
     171             : 
     172             :   /// The dimension
     173             :   const unsigned int _dim;
     174             :   /// The polar order
     175             :   const unsigned int _polar_order;
     176             :   /// The azimuthal order
     177             :   const unsigned int _azimuthal_order;
     178             :   /// The minimum mu
     179             :   const Real _mu_min;
     180             :   /// The maximum mu
     181             :   const Real _mu_max;
     182             : 
     183             :   /// Quadrature phi
     184             :   std::vector<Real> _phi;
     185             :   /// Quadrature mu
     186             :   std::vector<Real> _mu;
     187             :   /// Quadrature polar sin
     188             :   std::vector<Real> _polar_sin;
     189             :   /// Quadrature combined weights
     190             :   std::vector<Real> _w;
     191             : 
     192             :   /**
     193             :    * The current quadrature information.
     194             :    *
     195             :    * We need "current" information because the quadrature
     196             :    * can be rotated and because in 2D, there are cases where
     197             :    * projected directions that are the same in the 2D plane
     198             :    * are combined together with multiple weights and polar
     199             :    * sins.
     200             :    */
     201             :   ///@{
     202             :   std::vector<libMesh::Point> _current_directions;
     203             :   std::vector<std::vector<Real>> _current_weights;
     204             :   std::vector<std::vector<Real>> _current_polar_sins;
     205             :   ///@}
     206             : 
     207             :   /// The current rotation direction
     208             :   libMesh::Point _current_rotation_direction;
     209             : };