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 "CrossIC.h"
      11             : 
      12             : // Portions of this code Copyright 2007-2009 Roy Stogner
      13             : //
      14             : // Permission is hereby granted, free of charge, to any person obtaining
      15             : // a copy of this software and associated documentation files (the
      16             : //"Software"), to deal in the Software without restriction, including
      17             : // without limitation the rights to use, copy, modify, merge, publish,
      18             : // distribute, sublicense, and/or sell copies of the Software, and to
      19             : // permit persons to whom the Software is furnished to do so, subject to
      20             : // the following conditions:
      21             : //
      22             : // THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
      23             : // EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
      24             : // MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
      25             : // IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
      26             : // CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
      27             : // TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
      28             : // SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
      29             : 
      30             : registerMooseObject("PhaseFieldApp", CrossIC);
      31             : 
      32             : InputParameters
      33         349 : CrossIC::validParams()
      34             : {
      35         349 :   InputParameters params = C1ICBase::validParams();
      36         349 :   params.addClassDescription("Cross-shaped initial condition");
      37         698 :   params.addParam<Real>("x1", 0.0, "The x coordinate of the lower left-hand corner of the box");
      38         698 :   params.addParam<Real>("y1", 0.0, "The y coordinate of the lower left-hand corner of the box");
      39         698 :   params.addParam<Real>("x2", 1.0, "The x coordinate of the upper right-hand corner of the box");
      40         698 :   params.addParam<Real>("y2", 1.0, "The y coordinate of the upper right-hand corner of the box");
      41         349 :   return params;
      42           0 : }
      43             : 
      44         187 : CrossIC::CrossIC(const InputParameters & parameters)
      45             :   : C1ICBase(parameters),
      46         187 :     _x1(parameters.get<Real>("x1")),
      47         187 :     _y1(parameters.get<Real>("y1")),
      48         187 :     _x2(parameters.get<Real>("x2")),
      49         374 :     _y2(parameters.get<Real>("y2"))
      50             : {
      51         187 : }
      52             : 
      53             : Real
      54      813717 : CrossIC::value(const Point & p)
      55             : {
      56      813717 :   const Real x = (p(0) - _x1) / (_x2 - _x1);
      57      813717 :   const Real y = (p(1) - _y1) / (_y2 - _y1);
      58             : 
      59      813717 :   const Real cmax = _average + _amplitude;
      60      813717 :   const Real cmin = _average - _amplitude;
      61             : 
      62      813717 :   if (x > (0.5 + _length / 2.0 + _buffer + _interface))
      63             :   {
      64             :     return cmin;
      65             :   }
      66      735617 :   else if (x > 0.5 + _length / 2.0)
      67             :   {
      68       45822 :     if (y > 0.5 + _width / 2.0 + _buffer + _interface)
      69             :     {
      70             :       return cmin;
      71             :     }
      72       27936 :     else if (y > 0.5 + _width / 2.0)
      73             :     {
      74        2616 :       Real xd = x - .5 - _length / 2.0;
      75        2616 :       Real yd = y - .5 - _width / 2.0;
      76        2616 :       Real r = std::sqrt(xd * xd + yd * yd);
      77        2616 :       return interfaceValue(r);
      78             :     }
      79       25320 :     else if (y > 0.5 - _width / 2.0)
      80             :     {
      81        4818 :       return interfaceValue(x - .5 - _length / 2.0);
      82             :     }
      83       20502 :     else if (y > 0.5 - _width / 2. - _buffer - _interface)
      84             :     {
      85        2476 :       Real xd = x - .5 - _length / 2.0;
      86        2476 :       Real yd = y - .5 + _width / 2.0;
      87        2476 :       Real r = std::sqrt(xd * xd + yd * yd);
      88        2476 :       return interfaceValue(r);
      89             :     }
      90             :     else
      91             :     {
      92             :       return cmin;
      93             :     }
      94             :   }
      95      689795 :   else if (x > 0.5 + _width / 2. + 2. * _buffer + _interface)
      96             :   {
      97      136755 :     if (y > 0.5 + _width / 2 + _buffer + _interface)
      98             :       return cmin;
      99       83492 :     else if (y > 0.5 + _width / 2.0)
     100        7964 :       return interfaceValue(y - .5 - _width / 2.0);
     101       75528 :     else if (y > 0.5 - _width / 2.0)
     102             :       return cmax;
     103       61227 :     else if (y > 0.5 - _width / 2. - _buffer - _interface)
     104        7374 :       return interfaceValue(.5 - _width / 2. - y);
     105             :     else
     106             :       return cmin;
     107             :   }
     108      553040 :   else if (x > 0.5 + _width / 2.0)
     109             :   {
     110       67418 :     if (y > 0.5 + _length / 2. + _buffer + _interface)
     111             :     {
     112             :       return cmin;
     113             :     }
     114       60136 :     else if (y > 0.5 + _length / 2.0)
     115             :     {
     116        4244 :       Real xd = x - (.5 + _width / 2.0);
     117        4244 :       Real yd = y - (.5 + _length / 2.0);
     118        4244 :       Real r = std::sqrt(xd * xd + yd * yd);
     119        4244 :       return interfaceValue(r);
     120             :     }
     121       55892 :     else if (y > 0.5 + _width / 2. + 2. * _buffer + _interface)
     122             :     {
     123       12442 :       return interfaceValue(x - .5 - _width / 2.0);
     124             :     }
     125       43450 :     else if (y > 0.5 + _width / 2.0)
     126             :     {
     127        6440 :       Real xd = x - (.5 + _width / 2. + 2. * _buffer + _interface);
     128        6440 :       Real yd = y - (.5 + _width / 2. + 2. * _buffer + _interface);
     129        6440 :       Real r = std::sqrt(xd * xd + yd * yd);
     130        6440 :       return interfaceValue(2. * _buffer + _interface - r);
     131             :     }
     132       37010 :     else if (y > 0.5 - _width / 2.0)
     133             :     {
     134             :       return cmax;
     135             :     }
     136       30408 :     else if (y > 0.5 - _width / 2. - 2. * _buffer - _interface)
     137             :     {
     138        5672 :       Real xd = x - (.5 + _width / 2. + 2. * _buffer + _interface);
     139        5672 :       Real yd = y - (.5 - _width / 2. - 2. * _buffer - _interface);
     140        5672 :       Real r = std::sqrt(xd * xd + yd * yd);
     141        5672 :       return interfaceValue(2. * _buffer + _interface - r);
     142             :     }
     143       24736 :     else if (y > 0.5 - _length / 2.0)
     144             :     {
     145       12950 :       return interfaceValue(x - .5 - _width / 2.0);
     146             :     }
     147       11786 :     else if (y > 0.5 - _length / 2. - _buffer - _interface)
     148             :     {
     149        4504 :       Real xd = x - (.5 + _width / 2.0);
     150        4504 :       Real yd = y - (.5 - _length / 2.0);
     151        4504 :       Real r = std::sqrt(xd * xd + yd * yd);
     152        4504 :       return interfaceValue(r);
     153             :     }
     154             :     else
     155             :     {
     156             :       return cmin;
     157             :     }
     158             :   }
     159      485622 :   else if (x > 0.5 - _width / 2.0)
     160             :   {
     161       82197 :     if (y > 0.5 + _length / 2 + _buffer + _interface)
     162             :       return cmin;
     163       75185 :     else if (y > 0.5 + _length / 2.0)
     164        4818 :       return interfaceValue(y - .5 - _length / 2.0);
     165       70367 :     else if (y > 0.5 - _length / 2.0)
     166             :       return cmax;
     167       14093 :     else if (y > 0.5 - _length / 2. - _buffer - _interface)
     168        6244 :       return interfaceValue(.5 - _length / 2. - y);
     169             :     else
     170             :       return cmin;
     171             :   }
     172      403425 :   else if (x > 0.5 - _width / 2. - 2. * _buffer - _interface)
     173             :   {
     174       79038 :     if (y > 0.5 + _length / 2. + _buffer + _interface)
     175             :     {
     176             :       return cmin;
     177             :     }
     178       72556 :     else if (y > 0.5 + _length / 2.0)
     179             :     {
     180        3924 :       Real xd = x - (.5 - _width / 2.0);
     181        3924 :       Real yd = y - (.5 + _length / 2.0);
     182        3924 :       Real r = std::sqrt(xd * xd + yd * yd);
     183        3924 :       return interfaceValue(r);
     184             :     }
     185       68632 :     else if (y > 0.5 + _width / 2. + 2. * _buffer + _interface)
     186             :     {
     187       11738 :       return interfaceValue(.5 - _width / 2. - x);
     188             :     }
     189       56894 :     else if (y > 0.5 + _width / 2.0)
     190             :     {
     191        5864 :       Real xd = x - (.5 - _width / 2. - 2. * _buffer - _interface);
     192        5864 :       Real yd = y - (.5 + _width / 2. + 2. * _buffer + _interface);
     193        5864 :       Real r = std::sqrt(xd * xd + yd * yd);
     194        5864 :       return interfaceValue(2. * _buffer + _interface - r);
     195             :     }
     196       51030 :     else if (y > 0.5 - _width / 2.0)
     197             :     {
     198             :       return cmax;
     199             :     }
     200       43074 :     else if (y > 0.5 - _width / 2. - 2. * _buffer - _interface)
     201             :     {
     202        8636 :       Real xd = x - (.5 - _width / 2. - 2. * _buffer - _interface);
     203        8636 :       Real yd = y - (.5 - _width / 2. - 2. * _buffer - _interface);
     204        8636 :       Real r = std::sqrt(xd * xd + yd * yd);
     205        8636 :       return interfaceValue(2. * _buffer + _interface - r);
     206             :     }
     207       34438 :     else if (y > 0.5 - _length / 2.0)
     208             :     {
     209       19866 :       return interfaceValue(.5 - _width / 2. - x);
     210             :     }
     211       14572 :     else if (y > 0.5 - _length / 2. - _buffer - _interface)
     212             :     {
     213        6416 :       Real xd = x - (.5 - _width / 2.0);
     214        6416 :       Real yd = y - (.5 - _length / 2.0);
     215        6416 :       Real r = std::sqrt(xd * xd + yd * yd);
     216        6416 :       return interfaceValue(r);
     217             :     }
     218             :     else
     219             :     {
     220             :       return cmin;
     221             :     }
     222             :   }
     223      324387 :   else if (x > 0.5 - _length / 2.0)
     224             :   {
     225      177685 :     if (y > 0.5 + _width / 2 + _buffer + _interface)
     226             :       return cmin;
     227      123612 :     else if (y > 0.5 + _width / 2.0)
     228        8084 :       return interfaceValue(y - .5 - _width / 2.0);
     229      115528 :     else if (y > 0.5 - _width / 2.0)
     230             :       return cmax;
     231       97311 :     else if (y > 0.5 - _width / 2. - _buffer - _interface)
     232       12762 :       return interfaceValue(.5 - _width / 2. - y);
     233             :     else
     234             :       return cmin;
     235             :   }
     236      146702 :   else if (x > (.5 - _length / 2. - _buffer - _interface))
     237             :   {
     238       60232 :     if (y > 0.5 + _width / 2. + _buffer + _interface)
     239             :     {
     240             :       return cmin;
     241             :     }
     242       41076 :     else if (y > 0.5 + _width / 2.0)
     243             :     {
     244        2816 :       Real xd = x - (.5 - _length / 2.0);
     245        2816 :       Real yd = y - .5 - _width / 2.0;
     246        2816 :       Real r = std::sqrt(xd * xd + yd * yd);
     247        2816 :       return interfaceValue(r);
     248             :     }
     249       38260 :     else if (y > 0.5 - _width / 2.0)
     250             :     {
     251        6244 :       return interfaceValue(.5 - _length / 2. - x);
     252             :     }
     253       32016 :     else if (y > 0.5 - _width / 2. - _buffer - _interface)
     254             :     {
     255        4104 :       Real xd = x - (.5 - _length / 2.0);
     256        4104 :       Real yd = y - .5 + _width / 2.0;
     257        4104 :       Real r = std::sqrt(xd * xd + yd * yd);
     258        4104 :       return interfaceValue(r);
     259             :     }
     260             :     else
     261             :       return cmin;
     262             :   }
     263             :   else
     264             :     return cmin;
     265             : }
     266             : 
     267             : RealGradient
     268      100152 : CrossIC::gradient(const Point & p)
     269             : {
     270      100152 :   Point pxminus = p, pxplus = p, pyminus = p, pyplus = p;
     271             : 
     272      100152 :   pxminus(0) -= TOLERANCE;
     273      100152 :   pyminus(1) -= TOLERANCE;
     274      100152 :   pxplus(0) += TOLERANCE;
     275      100152 :   pyplus(1) += TOLERANCE;
     276             : 
     277      100152 :   Number uxminus = value(pxminus), uxplus = value(pxplus), uyminus = value(pyminus),
     278      100152 :          uyplus = value(pyplus);
     279             : 
     280      100152 :   return Gradient((uxplus - uxminus) / 2.0 / TOLERANCE, (uyplus - uyminus) / 2.0 / TOLERANCE);
     281             : }