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