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 "EulerAngleUpdater.h"
      11             : #include "GrainTrackerInterface.h"
      12             : #include "GrainForceAndTorqueInterface.h"
      13             : #include "RotationTensor.h"
      14             : #include "SystemBase.h"
      15             : 
      16             : registerMooseObject("PhaseFieldApp", EulerAngleUpdater);
      17             : 
      18             : InputParameters
      19          28 : EulerAngleUpdater::validParams()
      20             : {
      21          28 :   InputParameters params = EulerAngleProvider::validParams();
      22          28 :   params.addClassDescription(
      23             :       "Provide updated euler angles after rigid body rotation of the grains.");
      24          56 :   params.addRequiredParam<UserObjectName>("grain_tracker_object",
      25             :                                           "The FeatureFloodCount UserObject to get values from.");
      26          56 :   params.addRequiredParam<UserObjectName>("euler_angle_provider",
      27             :                                           "Name of Euler angle provider user object");
      28          56 :   params.addRequiredParam<UserObjectName>("grain_torques_object",
      29             :                                           "Name of Euler angle provider user object");
      30          56 :   params.addRequiredParam<VectorPostprocessorName>("grain_volumes",
      31             :                                                    "The feature volume VectorPostprocessorValue.");
      32          56 :   params.addParam<Real>("rotation_constant", 1.0, "Constant value characterizing grain rotation");
      33          28 :   return params;
      34           0 : }
      35             : 
      36          14 : EulerAngleUpdater::EulerAngleUpdater(const InputParameters & params)
      37             :   : EulerAngleProvider(params),
      38          14 :     _grain_tracker(getUserObject<GrainTrackerInterface>("grain_tracker_object")),
      39          14 :     _euler(getUserObject<EulerAngleProvider>("euler_angle_provider")),
      40          14 :     _grain_torque(getUserObject<GrainForceAndTorqueInterface>("grain_torques_object")),
      41          14 :     _grain_volumes(getVectorPostprocessorValue("grain_volumes", "feature_volumes")),
      42          28 :     _mr(getParam<Real>("rotation_constant")),
      43          28 :     _first_time(declareRestartableData<bool>("first_time_euler_update", true)),
      44          14 :     _first_time_recovered(_app.isRecovering()),
      45          28 :     _t_step_old(declareRestartableData<int>("euler_update_tstep_old", -1)),
      46          42 :     _angles_old(declareRestartableData<std::vector<EulerAngles>>("euler_angles_old"))
      47             : {
      48          14 : }
      49             : 
      50             : void
      51          44 : EulerAngleUpdater::initialize()
      52             : {
      53          44 :   const auto grain_num = _grain_tracker.getTotalFeatureCount();
      54             : 
      55          44 :   if (_first_time)
      56             :   {
      57          12 :     _angles.resize(grain_num);
      58          12 :     _angles_old.resize(grain_num);
      59          24 :     for (unsigned int i = 0; i < grain_num; ++i)
      60          12 :       _angles[i] = _euler.getEulerAngles(i); // Read initial euler angles
      61             :   }
      62             : 
      63          44 :   unsigned int angle_size = _angles.size();
      64          44 :   for (unsigned int i = angle_size; i < grain_num; ++i) // if new grains are created
      65           0 :     _angles.push_back(_euler.getEulerAngles(i));        // Assign initial euler angles
      66             : 
      67          88 :   for (unsigned int i = 0; i < grain_num; ++i)
      68             :   {
      69             :     // We dont want to use old angles on recovery
      70          44 :     if (!_first_time && !_first_time_recovered && !_fe_problem.converged(_sys.number()))
      71           0 :       _angles[i] = _angles_old[i];
      72          44 :     _first_time_recovered = false;
      73             : 
      74          44 :     RealGradient torque = _grain_torque.getTorqueValues()[i];
      75             : 
      76          44 :     if (i <= angle_size && _t_step > _t_step_old) // if new grains are created
      77          42 :       _angles_old[i] = _angles[i];
      78           2 :     else if (i > angle_size)
      79           0 :       _angles_old.push_back(_angles[i]);
      80             : 
      81          44 :     RotationTensor R0(_angles_old[i]); // RotationTensor as per old euler angles
      82             :     RealVectorValue torque_rotated =
      83          44 :         R0 * torque; // Applied torque is rotated to allign with old grain axes
      84             :     RealVectorValue omega =
      85          44 :         _mr / _grain_volumes[i] * torque_rotated; // Angular velocity as per old grain axes
      86             :     /**
      87             :      * Change in euler angles are obtained from the torque & angular velocities about the material
      88             :      * axes.
      89             :      * Change in phi1, Phi and phi2 are caused by rotation about z axis, x' axis & z'' axis,
      90             :      * respectively.
      91             :      * Components of the angular velocities across z, x' and z'' axes are obtained from the torque
      92             :      * values.
      93             :      * This yields change in euler angles due to grain rotation.
      94             :      */
      95             :     RealVectorValue angle_change;
      96          44 :     angle_change(0) = omega(2) * _dt;
      97          44 :     angle_change(1) =
      98          44 :         (omega(0) * std::cos(angle_change(0)) + omega(1) * std::sin(angle_change(0))) * _dt;
      99          44 :     angle_change(2) = (omega(0) * std::sin(angle_change(0)) * std::sin(angle_change(1)) -
     100          44 :                        omega(1) * std::cos(angle_change(0)) * std::sin(angle_change(1)) +
     101          44 :                        omega(2) * std::cos(angle_change(1))) *
     102          44 :                       _dt;
     103             :     angle_change *= (180.0 / libMesh::pi);
     104             : 
     105          44 :     RotationTensor R1(angle_change); // Rotation matrix due to torque
     106             :     /**
     107             :      * Final RotationMatrix = RotationMatrix due to applied torque X old RotationMatrix
     108             :      * Updated Euler angles are obtained by back-tracking the angles from the rotation matrix
     109             :      * For details about the componenets of the rotation matrix please refer to RotationTensor.C
     110             :      * Phi = acos(R33); phi1 = atan2(R31,-R32); phi2 = atan2(R13,R23) for phi != 0.0 por 180.0
     111             :      */
     112          44 :     RealTensorValue R = R1 * R0;
     113             : 
     114          44 :     if (R(2, 2) != 1.0 && R(2, 2) != -1.0) // checks if cos(Phi) = 1 or -1
     115             :     {
     116          44 :       _angles[i].phi1 = std::atan2(R(2, 0), -R(2, 1)) * (180.0 / libMesh::pi);
     117          44 :       _angles[i].Phi = std::acos(R(2, 2)) * (180.0 / libMesh::pi);
     118          44 :       _angles[i].phi2 = std::atan2(R(0, 2), R(1, 2)) * (180.0 / libMesh::pi);
     119             :     }
     120           0 :     else if (R(2, 2) == 1.0) // special case for Phi = 0.0
     121             :     {
     122           0 :       if (R0(2, 2) == 1.0)
     123             :         // when Phi_old = 0.0; all the rotations are about z axis and angles accumulates after each
     124             :         // rotation
     125           0 :         _angles[i].phi1 = _angles_old[i].phi1 + _angles_old[i].phi2 + angle_change(0);
     126             :       else
     127           0 :         _angles[i].phi1 = angle_change(0);
     128             :       // Comply with bunge euler angle definitions, 0.0 <= phi1 <= 360.0
     129           0 :       if (std::abs(_angles[i].phi1) > 360.0)
     130             :       {
     131           0 :         int laps = _angles[i].phi1 / 360.0;
     132           0 :         _angles[i].phi1 -= laps * 360.0;
     133             :       }
     134           0 :       _angles[i].Phi = 0.0;
     135           0 :       _angles[i].phi2 = -_angles[i].phi1 + std::atan2(R(0, 1), R(1, 1)) * (180.0 / libMesh::pi);
     136             :     }
     137             :     else
     138             :     {
     139           0 :       if (R0(2, 2) == 1.0)
     140           0 :         _angles[i].phi1 = _angles_old[i].phi1 + _angles_old[i].phi2 + angle_change(0);
     141             :       else
     142           0 :         _angles[i].phi1 = angle_change(0);
     143             :       // Comply with bunge euler angle definitions, 0.0 <= phi1 <= 360.0
     144           0 :       if (std::abs(_angles[i].phi1) > 360.0)
     145             :       {
     146           0 :         int laps = _angles[i].phi1 / 360.0;
     147           0 :         _angles[i].phi1 -= laps * 360.0;
     148             :       }
     149           0 :       _angles[i].Phi = 180.0;
     150           0 :       _angles[i].phi2 = _angles[i].phi1 - std::atan2(-R(0, 1), -R(1, 1)) * (180.0 / libMesh::pi);
     151             :     }
     152             : 
     153             :     // Following checks and updates are done only to comply with bunge euler angle definitions, 0.0
     154             :     // <= phi1/phi2 <= 360.0
     155          44 :     if (_angles[i].phi1 < 0.0)
     156          31 :       _angles[i].phi1 += 360.0;
     157          44 :     if (_angles[i].phi2 < 0.0)
     158          44 :       _angles[i].phi2 += 360.0;
     159          44 :     if (_angles[i].Phi < 0.0)
     160           0 :       mooseError("Euler angle out of range.");
     161             :   }
     162             : 
     163          44 :   _first_time = false;
     164          44 :   _t_step_old = _t_step;
     165          44 : }
     166             : 
     167             : const EulerAngles &
     168          13 : EulerAngleUpdater::getEulerAnglesOld(unsigned int i) const
     169             : {
     170             :   mooseAssert(i < getGrainNum(), "Requesting Euler angles for an invalid grain id");
     171          13 :   return _angles_old[i];
     172             : }