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 "libmesh/ignore_warnings.h"
      11             : #include "EulerAngles.h"
      12             : #include "MooseRandom.h"
      13             : 
      14         498 : EulerAngles::EulerAngles()
      15             : {
      16         498 :   phi1 = 0.0;
      17         498 :   Phi = 0.0;
      18         498 :   phi2 = 0.0;
      19         498 : }
      20             : 
      21      169016 : EulerAngles::EulerAngles(const Eigen::Quaternion<Real> & q)
      22             : {
      23      169016 :   phi1 = std::atan2((q.x() * q.z() + q.w() * q.y()), -(-q.w() * q.x() + q.y() * q.z())) *
      24             :          (180.0 / libMesh::pi);
      25             : 
      26             :   // avoid NAN value from sqrt(val) while computing Phi
      27      169016 :   auto val = 1.0 - std::pow(q.w() * q.w() - q.x() * q.x() - q.y() * q.y() + q.z() * q.z(), 2.0);
      28      169016 :   if (val < 0.0 && !MooseUtils::absoluteFuzzyEqual(val, 0.0))
      29           0 :     mooseError("Euler angle conversion is not successful due to invalid quaternion value.");
      30             : 
      31      169016 :   Phi = std::atan2(std::sqrt(std::abs(val)),
      32      169016 :                    q.w() * q.w() - q.x() * q.x() - q.y() * q.y() + q.z() * q.z()) *
      33             :         (180.0 / libMesh::pi);
      34      169016 :   phi2 = std::atan2((q.x() * q.z() - q.w() * q.y()), (q.w() * q.x() + q.y() * q.z())) *
      35             :          (180.0 / libMesh::pi);
      36             : 
      37             :   // Following checks and updates are done only to comply with bunge euler angle definitions, 0.0
      38             :   // <= phi1/phi2 <= 360.0
      39      169016 :   if (phi1 < 0.0)
      40       37410 :     phi1 += 360.0;
      41      169016 :   if (phi2 < 0.0)
      42           0 :     phi2 += 360.0;
      43      169016 :   if (Phi < 0.0)
      44           0 :     mooseError("Euler angle out of range.");
      45      169016 : }
      46             : 
      47         377 : EulerAngles::EulerAngles(const Real & v0, const Real & v1, const Real & v2)
      48             : {
      49         377 :   phi1 = v0;
      50         377 :   Phi = v1;
      51         377 :   phi2 = v2;
      52         377 : }
      53             : 
      54             : Eigen::Quaternion<Real>
      55           0 : EulerAngles::toQuaternion()
      56             : {
      57             :   Eigen::Quaternion<Real> q;
      58             : 
      59             :   Real cPhi1PlusPhi2, cphi, cPhi1MinusPhi2;
      60             :   Real sPhi1PlusPhi2, sphi, sPhi1MinusPhi2;
      61             : 
      62             :   /**
      63             :    * "NASA Mission Planning and Analysis Division.
      64             :    * "Euler Angles, Quaternions, and Transformation Matrices". NASA.
      65             :    */
      66             : 
      67           0 :   cPhi1PlusPhi2 = std::cos((phi1 * libMesh::pi / 180.0 + phi2 * libMesh::pi / 180.0) / 2.0);
      68           0 :   cphi = std::cos(Phi * libMesh::pi / 360.0);
      69           0 :   cPhi1MinusPhi2 = std::cos((phi1 * libMesh::pi / 180.0 - phi2 * libMesh::pi / 180.0) / 2.0);
      70             : 
      71           0 :   sPhi1PlusPhi2 = std::sin((phi1 * libMesh::pi / 180.0 + phi2 * libMesh::pi / 180.0) / 2.0);
      72           0 :   sphi = std::sin(Phi * libMesh::pi / 360.0);
      73           0 :   sPhi1MinusPhi2 = std::sin((phi1 * libMesh::pi / 180.0 - phi2 * libMesh::pi / 180.0) / 2.0);
      74             : 
      75           0 :   q.w() = cphi * cPhi1PlusPhi2;
      76           0 :   q.x() = sphi * cPhi1MinusPhi2;
      77           0 :   q.y() = sphi * sPhi1MinusPhi2;
      78           0 :   q.z() = cphi * sPhi1PlusPhi2;
      79             : 
      80           0 :   return q;
      81             : }
      82             : 
      83             : void
      84           0 : EulerAngles::random()
      85             : {
      86           0 :   phi1 = MooseRandom::rand() * 360.0;
      87           0 :   Phi = std::acos(1.0 - 2.0 * MooseRandom::rand()) / libMesh::pi * 180.0;
      88           0 :   phi2 = MooseRandom::rand() * 360;
      89           0 : }
      90             : 
      91             : void
      92           0 : EulerAngles::random(MooseRandom & random)
      93             : {
      94           0 :   phi1 = random.rand(0) * 360.0;
      95           0 :   Phi = std::acos(1.0 - 2.0 * random.rand(0)) / libMesh::pi * 180.0;
      96           0 :   phi2 = random.rand(0) * 360;
      97           0 : }
      98             : #include "libmesh/restore_warnings.h"