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