Line data    Source code

       1             : /********************************************************************/
       2             : /*                  SOFTWARE COPYRIGHT NOTIFICATION                 */
       3             : /*                             Cardinal                             */
       4             : /*                                                                  */
       5             : /*                  (c) 2021 UChicago Argonne, LLC                  */
       6             : /*                        ALL RIGHTS RESERVED                       */
       7             : /*                                                                  */
       8             : /*                 Prepared by UChicago Argonne, LLC                */
       9             : /*               Under Contract No. DE-AC02-06CH11357               */
      10             : /*                With the U. S. Department of Energy               */
      11             : /*                                                                  */
      12             : /*             Prepared by Battelle Energy Alliance, LLC            */
      13             : /*               Under Contract No. DE-AC07-05ID14517               */
      14             : /*                With the U. S. Department of Energy               */
      15             : /*                                                                  */
      16             : /*                 See LICENSE for full restrictions                */
      17             : /********************************************************************/
      18             : 
      19             : #include "SymmetryPointGenerator.h"
      20             : #include "GeometryUtils.h"
      21             : #include "UserErrorChecking.h"
      22             : 
      23             : registerMooseObject("CardinalApp", SymmetryPointGenerator);
      24             : 
      25             : InputParameters
      26         160 : SymmetryPointGenerator::validParams()
      27             : {
      28         160 :   InputParameters params = GeneralUserObject::validParams();
      29         320 :   params.addRequiredParam<Point>("normal", "Normal of the symmetry plane");
      30         320 :   params.addParam<Point>("rotation_axis",
      31             :                          "If rotationally symmetric, the axis about which to rotate. "
      32             :                          "If not specified, then the geometry is mirror-symmetric.");
      33         320 :   params.addRangeCheckedParam<Real>(
      34             :       "rotation_angle",
      35             :       "rotation_angle > 0 & rotation_angle <= 180",
      36             :       "If rotationally symmetric, the angle (degrees) from the 'normal' plane "
      37             :       "through which to rotate to form the symmetric wedge. If not specified, then the"
      38             :       "geometry is mirror-symmetric.");
      39         160 :   params.addClassDescription("Maps from a point (x, y, z) to a new point that is either "
      40             :                              "mirror-symmetric or rotationally-symmetric from the point.");
      41         160 :   return params;
      42           0 : }
      43             : 
      44          83 : SymmetryPointGenerator::SymmetryPointGenerator(const InputParameters & params)
      45         166 :   : ThreadedGeneralUserObject(params), _rotational_symmetry(isParamValid("rotation_axis"))
      46             : {
      47         332 :   checkJointParams(params, {"rotation_axis", "rotation_angle"}, "specifying rotational symmetry");
      48             : 
      49         249 :   _normal = geom_utils::unitVector(getParam<Point>("normal"), "normal");
      50             : 
      51          83 :   if (_rotational_symmetry)
      52             :   {
      53          40 :     const auto & angle = getParam<Real>("rotation_angle");
      54             : 
      55         120 :     _rotational_axis = geom_utils::unitVector(getParam<Point>("rotation_axis"), "rotation_axis");
      56             : 
      57             :     // the symmetry axis needs to be perpendicular to the plane normal
      58          40 :     if (!MooseUtils::absoluteFuzzyEqual(_rotational_axis * _normal, 0.0))
      59           3 :       paramError("rotation_axis", "The 'rotation_axis' must be perpendicular to the 'normal'!");
      60             : 
      61             :     // unit circle must be divisible by angle
      62          37 :     if (!MooseUtils::absoluteFuzzyEqual(fmod(360.0, angle), 0))
      63           3 :       paramError("rotation_angle",
      64             :                  "The unit circle must be evenly divisible by the 'rotation_angle'!");
      65             : 
      66          34 :     _angle = angle * M_PI / 180.0;
      67             : 
      68          34 :     _zero_theta = _normal.cross(_rotational_axis);
      69          34 :     _zero_theta = _zero_theta / _zero_theta.norm();
      70             : 
      71          34 :     _reflection_normal = geom_utils::rotatePointAboutAxis(_normal, -_angle / 2.0, _rotational_axis);
      72          34 :     _reflection_normal = _reflection_normal / _reflection_normal.norm();
      73             :   }
      74         243 : }
      75             : 
      76             : bool
      77      164643 : SymmetryPointGenerator::onPositiveSideOfPlane(const Point & p, const Point & normal) const
      78             : {
      79      164643 :   return p * normal > 0;
      80             : }
      81             : 
      82             : Point
      83       84134 : SymmetryPointGenerator::reflectPointAcrossPlane(const Point & p, const Point & normal) const
      84             : {
      85             :   Real coeff = -normal * p;
      86       84134 :   return p + 2.0 * coeff * normal;
      87             : }
      88             : 
      89             : int
      90       62668 : SymmetryPointGenerator::sector(const Point & p) const
      91             : {
      92       62668 :   Real theta = acos(p * _zero_theta / p.norm());
      93       62668 :   if (onPositiveSideOfPlane(p, _normal))
      94       31666 :     theta = 2.0 * M_PI - theta;
      95             : 
      96       62668 :   return theta / _angle;
      97             : }
      98             : 
      99             : Point
     100      164631 : SymmetryPointGenerator::transformPoint(const Point & p) const
     101             : {
     102      164631 :   Point pt = p;
     103             : 
     104      164631 :   if (_rotational_symmetry)
     105             :   {
     106             :     // first, find the closest point on the plane
     107             :     Real coeff = -_rotational_axis * p;
     108             :     Point vec_to_pt = p + coeff * _rotational_axis;
     109             : 
     110             :     // get the sector - we only need to do a transformation if not in the first sector
     111       62656 :     int s = sector(vec_to_pt);
     112       62656 :     if (s != 0)
     113             :     {
     114       51790 :       pt = geom_utils::rotatePointAboutAxis(p, s * _angle, _rotational_axis);
     115             : 
     116             :       // if the sector was odd, we also need to reflect the point about an axis
     117             :       // halfway between the symmetry plane and the zero-theta line
     118       51790 :       if (s % 2 != 0)
     119       31613 :         pt = reflectPointAcrossPlane(pt, _reflection_normal);
     120             :     }
     121             :   }
     122             :   else
     123             :   {
     124      101975 :     if (onPositiveSideOfPlane(pt, _normal))
     125       52521 :       pt = reflectPointAcrossPlane(pt, _normal);
     126             :   }
     127             : 
     128      164631 :   return pt;
     129             : }