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 "BimodalInverseSuperellipsoidsIC.h"
      11             : 
      12             : // MOOSE includes
      13             : #include "MooseMesh.h"
      14             : #include "MooseVariable.h"
      15             : 
      16             : registerMooseObject("PhaseFieldApp", BimodalInverseSuperellipsoidsIC);
      17             : 
      18             : InputParameters
      19          11 : BimodalInverseSuperellipsoidsIC::validParams()
      20             : {
      21          11 :   InputParameters params = BimodalSuperellipsoidsIC::validParams();
      22          11 :   params.addClassDescription("Bimodal size distribution of large particles (specified in input "
      23             :                              "file, value invalue) and small particles (placed randomly inside the "
      24             :                              "larger particles, value outvalue)");
      25          11 :   return params;
      26           0 : }
      27             : 
      28           6 : BimodalInverseSuperellipsoidsIC::BimodalInverseSuperellipsoidsIC(const InputParameters & parameters)
      29           6 :   : BimodalSuperellipsoidsIC(parameters)
      30             : {
      31           6 : }
      32             : 
      33             : void
      34           6 : BimodalInverseSuperellipsoidsIC::initialSetup()
      35             : {
      36           6 :   if (_size_variation_type == 2 && _size_variation > 0.0)
      37           0 :     paramError("size_variation",
      38             :                "If size_variation > 0.0, you must pass in a size_variation_type in "
      39             :                "BimodalInverseSuperellipsoidsIC");
      40             : 
      41           6 :   BimodalSuperellipsoidsIC::initialSetup();
      42           6 : }
      43             : 
      44             : Real
      45       10000 : BimodalInverseSuperellipsoidsIC::value(const Point & p)
      46             : {
      47       10000 :   Real value = _outvalue;
      48             :   Real val2 = 0.0;
      49             : 
      50             :   // First loop over the specified superellipsoids
      51       20000 :   for (unsigned int ellip = 0; ellip < _x_positions.size() && value != _invalue; ++ellip)
      52             :   {
      53       10000 :     val2 = computeSuperellipsoidValue(
      54             :         p, _centers[ellip], _as[ellip], _bs[ellip], _cs[ellip], _ns[ellip]);
      55       10000 :     if ((val2 > value && _invalue > _outvalue) || (val2 < value && _outvalue > _invalue))
      56             :       value = val2;
      57             :   }
      58             : 
      59             :   // Then loop over the randomly positioned particles and set value inside them back to outvalue
      60       90000 :   for (unsigned int ellip = _x_positions.size(); ellip < _x_positions.size() + _npart; ++ellip)
      61             :   {
      62       80000 :     val2 = computeSuperellipsoidInverseValue(
      63             :         p, _centers[ellip], _as[ellip], _bs[ellip], _cs[ellip], _ns[ellip]);
      64       80000 :     if ((val2 < value && _invalue > _outvalue) || (val2 > value && _outvalue > _invalue))
      65             :       value = val2;
      66             :   }
      67             : 
      68       10000 :   return value;
      69             : }
      70             : 
      71             : void
      72           6 : BimodalInverseSuperellipsoidsIC::computeSuperellipsoidCenters()
      73             : {
      74           6 :   _centers.resize(_x_positions.size() + _npart);
      75             : 
      76             :   // First place the specified (large) particles from the input file
      77          12 :   for (unsigned int i = 0; i < _x_positions.size(); ++i)
      78             :   {
      79           6 :     _centers[i](0) = _x_positions[i];
      80           6 :     _centers[i](1) = _y_positions[i];
      81           6 :     _centers[i](2) = _z_positions[i];
      82             :   }
      83             : 
      84             :   // Next place the randomly positioned (small) particles
      85          54 :   for (unsigned int i = _x_positions.size(); i < _x_positions.size() + _npart; ++i)
      86             :   {
      87             :     unsigned int num_tries = 0;
      88             : 
      89         612 :     while (num_tries < _max_num_tries)
      90             :     {
      91         612 :       num_tries++;
      92             : 
      93             :       RealTensorValue ran;
      94         612 :       ran(0, 0) = _random.rand(_tid);
      95         612 :       ran(1, 1) = _random.rand(_tid);
      96         612 :       ran(2, 2) = _random.rand(_tid);
      97             : 
      98         612 :       _centers[i] = _bottom_left + ran * _range;
      99             : 
     100             :       // check for collisions with the specified (large) and randomly placed small particles
     101        1086 :       for (unsigned int j = 0; j < i; ++j)
     102             :       {
     103             :         // Compute the distance r1 from the center of each specified superellipsoid to its
     104             :         // outside edge along the vector between the specified superellipsoid and the current
     105             :         // randomly positioned one.
     106             :         // This uses the equation for a superellipse in polar coordinates and substitutes
     107             :         // distances for sin, cos functions.
     108        1038 :         Point dist_vec = _mesh.minPeriodicVector(_var.number(), _centers[i], _centers[j]);
     109        1038 :         const Real dist = dist_vec.norm();
     110             : 
     111             :         // First calculate rmn1 = r1^(-n), replacing sin, cos functions with distances
     112        1038 :         Real rmn1 = (std::pow(std::abs(dist_vec(0) / dist / _as[j]), _ns[j]) +
     113        1038 :                      std::pow(std::abs(dist_vec(1) / dist / _bs[j]), _ns[j]) +
     114        1038 :                      std::pow(std::abs(dist_vec(2) / dist / _cs[j]), _ns[j]));
     115             :         // Then calculate r1 from rmn1
     116        1038 :         const Real r1 = std::pow(rmn1, (-1.0 / _ns[j]));
     117             : 
     118             :         // Now calculate the distance r2 from the center of the randomly placed
     119             :         // superellipsoid to its outside edge in the same manner
     120        1038 :         Real rmn2 = (std::pow(std::abs(dist_vec(0) / dist / _as[i]), _ns[i]) +
     121        1038 :                      std::pow(std::abs(dist_vec(1) / dist / _bs[i]), _ns[i]) +
     122        1038 :                      std::pow(std::abs(dist_vec(2) / dist / _cs[i]), _ns[i]));
     123        1038 :         const Real r2 = std::pow(rmn2, (-1.0 / _ns[i]));
     124             : 
     125        1038 :         if (j < _x_positions.size())
     126             :         {
     127         612 :           if (r1 - dist - r2 < _large_spac)
     128         564 :             goto fail;
     129             :         }
     130             :         else
     131             :         {
     132         426 :           if (dist - r1 - r2 < _small_spac)
     133         114 :             goto fail;
     134             :         }
     135             :       }
     136             : 
     137             :       // accept the position of the new center
     138          48 :       goto accept;
     139             : 
     140             :     // retry a new position until tries are exhausted
     141             :     fail:
     142             :       continue;
     143             :     }
     144             : 
     145           0 :     if (num_tries == _max_num_tries)
     146           0 :       mooseError("Too many tries in MultiSmoothCircleIC");
     147             : 
     148           0 :   accept:
     149             :     continue;
     150          48 :   }
     151           6 : }