www.mooseframework.org
TricrystalTripleJunctionIC.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 
9 #include "MooseRandom.h"
10 #include "MooseMesh.h"
11 #include "MathUtils.h"
12 
13 template <>
14 InputParameters
16 {
17  InputParameters params = validParams<InitialCondition>();
18  params.addClassDescription("Tricrystal with a triple junction");
19  params.addRequiredParam<unsigned int>("op_num", "Number of grain order parameters");
20  params.addRequiredParam<unsigned int>("op_index", "Index for the current grain order parameter");
21  params.addParam<Real>("theta1", 135.0, "Angle of first grain at triple junction in degrees");
22  params.addParam<Real>("theta2", 135.0, "Angle of second grain at triple junction in degrees");
23  params.addParam<Point>(
24  "junction",
25  "The point where the triple junction is located. Default is the center of the mesh");
26  return params;
27 }
28 
29 TricrystalTripleJunctionIC::TricrystalTripleJunctionIC(const InputParameters & parameters)
30  : InitialCondition(parameters),
31  _mesh(_fe_problem.mesh()),
32  _op_num(getParam<unsigned int>("op_num")),
33  _op_index(getParam<unsigned int>("op_index")),
34  _theta1(getParam<Real>("theta1")),
35  _theta2(getParam<Real>("theta2"))
36 {
37  if (_op_num != 3)
38  mooseError("Tricrystal ICs must have op_num = 3");
39 
40  if (_theta1 + _theta2 >= 360.0)
41  mooseError("Sum of the angles must total less than 360 degrees");
42 
43  // Make sure that _junction is in the domain
44  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
45  {
46  if ((_mesh.getMinInDimension(i) > _junction(i)) || (_mesh.getMaxInDimension(i) < _junction(i)))
47  mooseError("Triple junction out of bounds");
48  }
49 
50  // Default junction point is the center
51  if (!parameters.isParamValid("junction"))
52  {
53  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
54  _junction(i) = (_mesh.getMaxInDimension(i) - _mesh.getMinInDimension(i)) / 2.0;
55  }
56  else
57  _junction = getParam<Point>("junction");
58 
59  // Convert the angles to radians
60  _theta1 = _theta1 * libMesh::pi / 180.0;
61  _theta2 = _theta2 * libMesh::pi / 180.0;
62 
63  // Change the first angle to be measured from the +x axis
64  _theta1 = 3.0 * libMesh::pi / 2.0 - _theta1;
65 
66  // Change the third angle to be measured from the +x axis
67  _theta2 = _theta2 - libMesh::pi / 2.0;
68 
69  // Only compute the tangent once for computational efficiency
70  _tan_theta1 = std::tan(_theta1);
71  _tan_theta2 = std::tan(_theta2);
72 }
73 
74 Real
76 {
77  /*
78  * This does all the work to create a triple junction that looks like the letter Y
79  */
80  Real dist_left; // The distance from the point to the line specified by _theta1
81  Real dist_right; // The distance from the point to the line specified by _theta2
82 
83  // Check if the point is above or below the left-most line
84  // Function to use is y = _junction(1) + (x - _junction(0)) * std::tan(libMesh::pi/2.0 - _theta1)
85  if (_theta1 == 0) // Handle tan(0) case
86  dist_left = p(1) - _junction(1);
87  else
88  dist_left = p(1) - (_junction(1) + (p(0) - _junction(0)) * _tan_theta1);
89 
90  // Check if the point is above or below the right-most line
91  // Function to use is y = _junction(1) + (x - _junction(0))*std::tan(-(libMesh::pi/2.0 - _theta2))
92  if (_theta2 == 0) // Handle tan(0) classes
93  dist_right = p(1) - _junction(1);
94  else
95  dist_right = p(1) - (_junction(1) + (p(0) - _junction(0)) * _tan_theta2);
96 
97  // Check if the point is to the left or right of the middle line
98  Real dist_center = p(0) - _junction(0); // Negative value if the point is to the left
99 
100  if (_tan_theta1 > 0 && _theta1 <= libMesh::pi / 2.0) // Case for large left grain
101  {
102  /*
103  * There's a lot going on here. The first statement tells MOOSE to check and
104  * see if the current point is above the line defined by the first angle, only
105  * if it is past the center line. All other points to the left of the center
106  * line are going to be part of the 1st grain (_op_index == 1)
107  * The second statement defines the second grain by the line defined by _theta2
108  * and marks everything below that line, and to the right of the center line.
109  * The third statement takes care of everything in between.
110  */
111  if ((((dist_left >= 0 && dist_center >= 0) || (dist_center < 0)) && _op_index == 1) ||
112  (dist_right <= 0 && dist_center > 0 && _op_index == 2) ||
113  (dist_left < 0 && dist_center > 0 && dist_right > 0 && _op_index == 3))
114  return 1.0;
115  else
116  return 0.0;
117  }
118 
119  // This does a similar thing as the above case, but switches it for the right and left grains
120  else if (_tan_theta2 < 0 && _theta2 >= libMesh::pi / 2.0) // Case for large right grain
121  {
122  if ((dist_left <= 0 && dist_center <= 0 && _op_index == 1) ||
123  (((dist_right >= 0 && dist_center <= 0) || (dist_center > 0)) && _op_index == 2) ||
124  (dist_left > 0 && dist_right < 0 && dist_center < 0 && _op_index == 3))
125  return 1.0;
126  else
127  return 0.0;
128  }
129 
130  else // all other cases
131  {
132  if ((dist_left <= 0 && dist_center <= 0 && _op_index == 1) || // First grain
133  (dist_right <= 0 && dist_center > 0 && _op_index == 2) || // Second grain
134  (((dist_left > 0 && dist_center < 0) || (dist_right > 0 && dist_center >= 0)) &&
135  _op_index == 3)) // Third grain
136  return 1.0;
137  else
138  return 0.0;
139  }
140 }
const unsigned int _op_num
Number of order parameters.
virtual Real value(const Point &p)
Real _theta2
Angle of third grain at triple junction in radians.
Real _tan_theta1
tangent of the first angle after a shift of pi/2
Real _theta1
Angle of first grain at triple junction in radians.
InputParameters validParams< TricrystalTripleJunctionIC >()
Real _tan_theta2
tangent of the second angle after a shift of pi/2
Point _junction
Point where the triple junction occurs.
TricrystalTripleJunctionIC(const InputParameters &parameters)