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 "CubicTransition.h"
      11             : #include "MooseError.h"
      12             : 
      13             : #include "libmesh/dense_matrix.h"
      14             : #include "libmesh/dense_vector.h"
      15             : 
      16          98 : CubicTransition::CubicTransition(const Real & x_center, const Real & transition_width)
      17             :   : SmoothTransition(x_center, transition_width),
      18             : 
      19          98 :     _A(0.0),
      20          98 :     _B(0.0),
      21          98 :     _C(0.0),
      22          98 :     _D(0.0),
      23             : 
      24          98 :     _initialized(false)
      25             : {
      26          98 : }
      27             : 
      28             : void
      29          98 : CubicTransition::initialize(const Real & f1_end_value,
      30             :                             const Real & f2_end_value,
      31             :                             const Real & df1dx_end_value,
      32             :                             const Real & df2dx_end_value)
      33             : {
      34             :   // compute cubic polynomial coefficients
      35             : 
      36          98 :   DenseMatrix<Real> mat(4, 4);
      37             : 
      38          98 :   mat(0, 0) = std::pow(_x1, 3);
      39          98 :   mat(0, 1) = std::pow(_x1, 2);
      40          98 :   mat(0, 2) = _x1;
      41          98 :   mat(0, 3) = 1.0;
      42             : 
      43          98 :   mat(1, 0) = std::pow(_x2, 3);
      44          98 :   mat(1, 1) = std::pow(_x2, 2);
      45          98 :   mat(1, 2) = _x2;
      46          98 :   mat(1, 3) = 1.0;
      47             : 
      48          98 :   mat(2, 0) = 3.0 * std::pow(_x1, 2);
      49          98 :   mat(2, 1) = 2.0 * _x1;
      50          98 :   mat(2, 2) = 1.0;
      51          98 :   mat(2, 3) = 0.0;
      52             : 
      53          98 :   mat(3, 0) = 3.0 * std::pow(_x2, 2);
      54          98 :   mat(3, 1) = 2.0 * _x2;
      55          98 :   mat(3, 2) = 1.0;
      56          98 :   mat(3, 3) = 0.0;
      57             : 
      58          98 :   DenseVector<Real> rhs(4);
      59          98 :   rhs(0) = f1_end_value;
      60          98 :   rhs(1) = f2_end_value;
      61          98 :   rhs(2) = df1dx_end_value;
      62          98 :   rhs(3) = df2dx_end_value;
      63             : 
      64          98 :   DenseVector<Real> coefs(4);
      65          98 :   mat.lu_solve(rhs, coefs);
      66             : 
      67          98 :   _A = coefs(0);
      68          98 :   _B = coefs(1);
      69          98 :   _C = coefs(2);
      70          98 :   _D = coefs(3);
      71             : 
      72          98 :   _initialized = true;
      73          98 : }
      74             : 
      75             : Real
      76        2342 : CubicTransition::value(const Real & x, const Real & f1, const Real & f2) const
      77             : {
      78             :   mooseAssert(_initialized, "initialize() must be called.");
      79             : 
      80        2342 :   if (x <= _x1)
      81         533 :     return f1;
      82        1809 :   else if (x >= _x2)
      83        1019 :     return f2;
      84             :   else
      85         790 :     return _A * std::pow(x, 3) + _B * std::pow(x, 2) + _C * x + _D;
      86             : }
      87             : 
      88             : Real
      89        1042 : CubicTransition::derivative(const Real & x, const Real & df1dx, const Real & df2dx) const
      90             : {
      91             :   mooseAssert(_initialized, "initialize() must be called.");
      92             : 
      93        1042 :   if (x <= _x1)
      94         381 :     return df1dx;
      95         661 :   else if (x >= _x2)
      96         375 :     return df2dx;
      97             :   else
      98         286 :     return 3.0 * _A * std::pow(x, 2) + 2.0 * _B * x + _C;
      99             : }