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        1258 : 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             :   // ensure that vec is normalised
      28        1258 :   vec /= std::sqrt(vec * vec);
      29             : 
      30             :   // construct v0 and v1 to be orthonormal to vec
      31             :   // and form a RH basis, that is, so v1 x vec = v0
      32             : 
      33             :   // Use Gram-Schmidt method to find v1.
      34        1258 :   GenericRealVectorValue<is_ad> v1;
      35             :   // Need a prototype for v1 first, and this is done by looking at the smallest component of vec
      36        1258 :   GenericRealVectorValue<is_ad> w(std::abs(vec(0)), std::abs(vec(1)), std::abs(vec(2)));
      37        1258 :   if ((w(2) >= w(1) && w(1) >= w(0)) || (w(1) >= w(2) && w(2) >= w(0)))
      38             :     // vec(0) is the smallest component
      39        1063 :     v1(0) = 1;
      40         195 :   else if ((w(2) >= w(0) && w(0) >= w(1)) || (w(0) >= w(2) && w(2) >= w(1)))
      41             :     // vec(1) is the smallest component
      42         185 :     v1(1) = 1;
      43             :   else
      44             :     // vec(2) is the smallest component
      45          10 :     v1(2) = 1;
      46             :   // now Gram-Schmidt
      47        1258 :   v1 -= (v1 * vec) * vec;
      48        1258 :   v1 /= std::sqrt(v1 * v1);
      49             : 
      50             :   // now use v0 = v1 x vec
      51        1258 :   GenericRealVectorValue<is_ad> v0;
      52        1258 :   v0(0) = v1(1) * vec(2) - v1(2) * vec(1);
      53        1258 :   v0(1) = v1(2) * vec(0) - v1(0) * vec(2);
      54        1258 :   v0(2) = v1(0) * vec(1) - v1(1) * vec(0);
      55             : 
      56             :   // the desired rotation matrix is just
      57        1258 :   GenericRealTensorValue<is_ad> rot(
      58        1258 :       v0(0), v0(1), v0(2), v1(0), v1(1), v1(2), vec(0), vec(1), vec(2));
      59        2516 :   return rot;
      60             : }
      61             : 
      62             : /// provides a rotation matrix that will rotate the vector vec1 to vec2
      63             : template <bool is_ad = false>
      64             : GenericRealTensorValue<is_ad>
      65         170 : rotVec1ToVec2(GenericRealVectorValue<is_ad> vec1, GenericRealVectorValue<is_ad> vec2)
      66             : // provides a rotation matrix that will rotate the vector vec1 to the vector vec2
      67             : {
      68         170 :   GenericRealTensorValue<is_ad> rot1_to_z = rotVecToZ<is_ad>(vec1);
      69         170 :   GenericRealTensorValue<is_ad> rot2_to_z = rotVecToZ<is_ad>(vec2);
      70         340 :   return rot2_to_z.transpose() * rot1_to_z;
      71         170 : }
      72             : 
      73             : /// provides a rotation matrix that will rotate the vector vec1 to the [1,0,0], assuming vec1[2]==0
      74             : template <bool is_ad = false>
      75             : GenericRealTensorValue<is_ad>
      76          12 : rotVec2DToX(const GenericRealVectorValue<is_ad> & vec)
      77             : // provides a rotation matrix that will rotate the vector `vec` to the [1,0,0], assuming vec[2]==0
      78             : {
      79          12 :   const GenericReal<is_ad> theta = std::atan2(vec(1), vec(0));
      80          12 :   const GenericReal<is_ad> st = std::sin(theta);
      81          12 :   const GenericReal<is_ad> ct = std::cos(theta);
      82          12 :   return GenericRealTensorValue<is_ad>(ct, st, 0., -st, ct, 0., 0., 0., 1.);
      83             : }
      84             : }