www.mooseframework.org
EulerAngleUpdater.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
7 #include "EulerAngleUpdater.h"
10 #include "RotationTensor.h"
11 
12 template <>
13 InputParameters
15 {
16  InputParameters params = validParams<EulerAngleProvider>();
17  params.addClassDescription(
18  "Provide updated euler angles after rigid body rotation of the grains.");
19  params.addRequiredParam<UserObjectName>("grain_tracker_object",
20  "The FeatureFloodCount UserObject to get values from.");
21  params.addRequiredParam<UserObjectName>("euler_angle_provider",
22  "Name of Euler angle provider user object");
23  params.addRequiredParam<UserObjectName>("grain_torques_object",
24  "Name of Euler angle provider user object");
25  params.addRequiredParam<VectorPostprocessorName>("grain_volumes",
26  "The feature volume VectorPostprocessorValue.");
27  params.addParam<Real>("rotation_constant", 1.0, "Constant value characterizing grain rotation");
28  return params;
29 }
30 
31 EulerAngleUpdater::EulerAngleUpdater(const InputParameters & params)
32  : EulerAngleProvider(params),
33  _grain_tracker(getUserObject<GrainTrackerInterface>("grain_tracker_object")),
34  _euler(getUserObject<EulerAngleProvider>("euler_angle_provider")),
35  _grain_torque(getUserObject<GrainForceAndTorqueInterface>("grain_torques_object")),
36  _grain_volumes(getVectorPostprocessorValue("grain_volumes", "feature_volumes")),
37  _mr(getParam<Real>("rotation_constant")),
38  _first_time(true)
39 {
40 }
41 
42 void
44 {
45  const auto grain_num = _grain_tracker.getTotalFeatureCount();
46 
47  if (_first_time)
48  {
49  _angles.resize(grain_num);
50  _angles_old.resize(grain_num);
51  for (unsigned int i = 0; i < grain_num; ++i)
52  _angles[i] = _euler.getEulerAngles(i); // Read initial euler angles
53  }
54 
55  unsigned int angle_size = _angles.size();
56  for (unsigned int i = angle_size; i < grain_num; ++i) // if new grains are created
57  _angles.push_back(_euler.getEulerAngles(i)); // Assign initial euler angles
58 
59  for (unsigned int i = 0; i < grain_num; ++i)
60  {
61  if (!_first_time && !_fe_problem.converged())
62  _angles[i] = _angles_old[i];
63 
64  RealGradient torque = _grain_torque.getTorqueValues()[i];
65 
66  if (i <= angle_size) // if new grains are created
67  _angles_old[i] = _angles[i];
68  else
69  _angles_old.push_back(_angles[i]);
70 
71  RotationTensor R0(_angles_old[i]); // RotationTensor as per old euler angles
72  RealVectorValue torque_rotated =
73  R0 * torque; // Applied torque is rotated to allign with old grain axes
74  RealVectorValue omega =
75  _mr / _grain_volumes[i] * torque_rotated; // Angular velocity as per old grain axes
85  RealVectorValue angle_change;
86  angle_change(0) = omega(2) * _dt;
87  angle_change(1) =
88  (omega(0) * std::cos(angle_change(0)) + omega(1) * std::sin(angle_change(0))) * _dt;
89  angle_change(2) = (omega(0) * std::sin(angle_change(0)) * std::sin(angle_change(1)) -
90  omega(1) * std::cos(angle_change(0)) * std::sin(angle_change(1)) +
91  omega(2) * std::cos(angle_change(1))) *
92  _dt;
93  angle_change *= (180.0 / libMesh::pi);
94 
95  RotationTensor R1(angle_change); // Rotation matrix due to torque
102  RealTensorValue R = R1 * R0;
103 
104  if (R(2, 2) != 1.0 && R(2, 2) != -1.0) // checks if cos(Phi) = 1 or -1
105  {
106  _angles[i].phi1 = std::atan2(R(2, 0), -R(2, 1)) * (180.0 / libMesh::pi);
107  _angles[i].Phi = std::acos(R(2, 2)) * (180.0 / libMesh::pi);
108  _angles[i].phi2 = std::atan2(R(0, 2), R(1, 2)) * (180.0 / libMesh::pi);
109  }
110  else if (R(2, 2) == 1.0) // special case for Phi = 0.0
111  {
112  if (R0(2, 2) == 1.0)
113  // when Phi_old = 0.0; all the rotations are about z axis and angles accumulates after each
114  // rotation
115  _angles[i].phi1 = _angles_old[i].phi1 + _angles_old[i].phi2 + angle_change(0);
116  else
117  _angles[i].phi1 = angle_change(0);
118  // Comply with bunge euler angle definitions, 0.0 <= phi1 <= 360.0
119  if (std::abs(_angles[i].phi1) > 360.0)
120  {
121  int laps = _angles[i].phi1 / 360.0;
122  _angles[i].phi1 -= laps * 360.0;
123  }
124  _angles[i].Phi = 0.0;
125  _angles[i].phi2 = -_angles[i].phi1 + std::atan2(R(0, 1), R(1, 1)) * (180.0 / libMesh::pi);
126  }
127  else
128  {
129  if (R0(2, 2) == 1.0)
130  _angles[i].phi1 = _angles_old[i].phi1 + _angles_old[i].phi2 + angle_change(0);
131  else
132  _angles[i].phi1 = angle_change(0);
133  // Comply with bunge euler angle definitions, 0.0 <= phi1 <= 360.0
134  if (std::abs(_angles[i].phi1) > 360.0)
135  {
136  int laps = _angles[i].phi1 / 360.0;
137  _angles[i].phi1 -= laps * 360.0;
138  }
139  _angles[i].Phi = 180.0;
140  _angles[i].phi2 = _angles[i].phi1 - std::atan2(-R(0, 1), -R(1, 1)) * (180.0 / libMesh::pi);
141  }
142 
143  // Following checks and updates are done only to comply with bunge euler angle definitions, 0.0
144  // <= phi1/phi2 <= 360.0
145  if (_angles[i].phi1 < 0.0)
146  _angles[i].phi1 += 360.0;
147  if (_angles[i].phi2 < 0.0)
148  _angles[i].phi2 += 360.0;
149  if (_angles[i].Phi < 0.0)
150  mooseError("Euler angle out of range.");
151  }
152 
153  _first_time = false;
154 }
155 
156 unsigned int
158 {
159  return _angles.size();
160 }
161 
162 const EulerAngles &
164 {
165  mooseAssert(i < getGrainNum(), "Requesting Euler angles for an invalid grain id");
166  return _angles[i];
167 }
168 
169 const EulerAngles &
171 {
172  mooseAssert(i < getGrainNum(), "Requesting Euler angles for an invalid grain id");
173  return _angles_old[i];
174 }
virtual const EulerAngles & getEulerAnglesOld(unsigned int) const
const GrainTrackerInterface & _grain_tracker
const GrainForceAndTorqueInterface & _grain_torque
This class defines the interface for the GrainTracking objects.
This class provides interface for extracting the forces and torques computed in other UserObjects...
virtual std::size_t getTotalFeatureCount() const =0
Returns a number large enough to contain the largest ID for all grains in use.
InputParameters validParams< EulerAngleUpdater >()
const VectorPostprocessorValue & _grain_volumes
virtual unsigned int getGrainNum() const override
InputParameters validParams< EulerAngleProvider >()
const EulerAngleProvider & _euler
This is a RealTensor version of a rotation matrix It is instantiated with the Euler angles...
EulerAngleUpdater(const InputParameters &parameters)
virtual void initialize() override
std::vector< EulerAngles > _angles_old
Euler angle triplet.
Definition: EulerAngles.h:20
virtual const std::vector< RealGradient > & getTorqueValues() const =0
Abstract base class for user objects that implement the Euler Angle provider interface.
std::vector< EulerAngles > _angles
virtual const EulerAngles & getEulerAngles(unsigned int) const override
virtual const EulerAngles & getEulerAngles(unsigned int) const =0