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             : #include "Moose.h"
      13             : #include "libmesh/vector_value.h"
      14             : #include "libmesh/tensor_value.h"
      15             : #include "MooseTypes.h"
      16             : /**
      17             :  * Utility functions to return rotations matrics
      18             :  */
      19             : namespace RotationMatrix
      20             : {
      21             : /// provides a rotation matrix that will rotate the vector vec to the z axis (the "2" direction)
      22             : template <bool is_ad = false>
      23             : GenericRealTensorValue<is_ad>
      24        1432 : rotVecToZ(GenericRealVectorValue<is_ad> vec)
      25             : // provides a rotation matrix that will rotate the vector vec to the z axis (the "2" direction)
      26             : {
      27             :   using std::sqrt, std::abs;
      28             :   // ensure that vec is normalised
      29        1432 :   vec /= sqrt(vec * vec);
      30             : 
      31             :   // construct v0 and v1 to be orthonormal to vec
      32             :   // and form a RH basis, that is, so v1 x vec = v0
      33             : 
      34             :   // Use Gram-Schmidt method to find v1.
      35        1432 :   GenericRealVectorValue<is_ad> v1;
      36             :   // Need a prototype for v1 first, and this is done by looking at the smallest component of vec
      37        1432 :   GenericRealVectorValue<is_ad> w(abs(vec(0)), abs(vec(1)), abs(vec(2)));
      38        1432 :   if ((w(2) >= w(1) && w(1) >= w(0)) || (w(1) >= w(2) && w(2) >= w(0)))
      39             :     // vec(0) is the smallest component
      40        1198 :     v1(0) = 1;
      41         234 :   else if ((w(2) >= w(0) && w(0) >= w(1)) || (w(0) >= w(2) && w(2) >= w(1)))
      42             :     // vec(1) is the smallest component
      43         214 :     v1(1) = 1;
      44             :   else
      45             :     // vec(2) is the smallest component
      46          20 :     v1(2) = 1;
      47             :   // now Gram-Schmidt
      48        1432 :   v1 -= (v1 * vec) * vec;
      49        1432 :   v1 /= sqrt(v1 * v1);
      50             : 
      51             :   // now use v0 = v1 x vec
      52        1432 :   GenericRealVectorValue<is_ad> v0;
      53        1432 :   v0(0) = v1(1) * vec(2) - v1(2) * vec(1);
      54        1432 :   v0(1) = v1(2) * vec(0) - v1(0) * vec(2);
      55        1432 :   v0(2) = v1(0) * vec(1) - v1(1) * vec(0);
      56             : 
      57             :   // the desired rotation matrix is just
      58        1432 :   GenericRealTensorValue<is_ad> rot(
      59        1432 :       v0(0), v0(1), v0(2), v1(0), v1(1), v1(2), vec(0), vec(1), vec(2));
      60        2864 :   return rot;
      61             : }
      62             : 
      63             : /// provides a rotation matrix that will rotate the vector vec1 to vec2
      64             : template <bool is_ad = false>
      65             : GenericRealTensorValue<is_ad>
      66         206 : rotVec1ToVec2(GenericRealVectorValue<is_ad> vec1, GenericRealVectorValue<is_ad> vec2)
      67             : // provides a rotation matrix that will rotate the vector vec1 to the vector vec2
      68             : {
      69         206 :   GenericRealTensorValue<is_ad> rot1_to_z = rotVecToZ<is_ad>(vec1);
      70         206 :   GenericRealTensorValue<is_ad> rot2_to_z = rotVecToZ<is_ad>(vec2);
      71         412 :   return rot2_to_z.transpose() * rot1_to_z;
      72         206 : }
      73             : 
      74             : /// provides a rotation matrix that will rotate the vector vec1 to the [1,0,0], assuming vec1[2]==0
      75             : template <bool is_ad = false>
      76             : GenericRealTensorValue<is_ad>
      77          24 : rotVec2DToX(const GenericRealVectorValue<is_ad> & vec)
      78             : // provides a rotation matrix that will rotate the vector `vec` to the [1,0,0], assuming vec[2]==0
      79             : {
      80             :   using std::atan2, std::sin, std::cos;
      81          24 :   const GenericReal<is_ad> theta = atan2(vec(1), vec(0));
      82          24 :   const GenericReal<is_ad> st = sin(theta);
      83          24 :   const GenericReal<is_ad> ct = cos(theta);
      84          24 :   return GenericRealTensorValue<is_ad>(ct, st, 0., -st, ct, 0., 0., 0., 1.);
      85             : }
      86             : }