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         164 : CubicTransition::CubicTransition(const Real & x_center, const Real & transition_width)
      17             :   : SmoothTransition(x_center, transition_width),
      18             : 
      19         164 :     _A(0.0),
      20         164 :     _B(0.0),
      21         164 :     _C(0.0),
      22         164 :     _D(0.0),
      23             : 
      24         164 :     _initialized(false)
      25             : {
      26         164 : }
      27             : 
      28             : void
      29         164 : 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         164 :   DenseMatrix<Real> mat(4, 4);
      37             : 
      38         164 :   mat(0, 0) = std::pow(_x1, 3);
      39         164 :   mat(0, 1) = std::pow(_x1, 2);
      40         164 :   mat(0, 2) = _x1;
      41         164 :   mat(0, 3) = 1.0;
      42             : 
      43         164 :   mat(1, 0) = std::pow(_x2, 3);
      44         164 :   mat(1, 1) = std::pow(_x2, 2);
      45         164 :   mat(1, 2) = _x2;
      46         164 :   mat(1, 3) = 1.0;
      47             : 
      48         164 :   mat(2, 0) = 3.0 * std::pow(_x1, 2);
      49         164 :   mat(2, 1) = 2.0 * _x1;
      50         164 :   mat(2, 2) = 1.0;
      51         164 :   mat(2, 3) = 0.0;
      52             : 
      53         164 :   mat(3, 0) = 3.0 * std::pow(_x2, 2);
      54         164 :   mat(3, 1) = 2.0 * _x2;
      55         164 :   mat(3, 2) = 1.0;
      56         164 :   mat(3, 3) = 0.0;
      57             : 
      58         164 :   DenseVector<Real> rhs(4);
      59         164 :   rhs(0) = f1_end_value;
      60         164 :   rhs(1) = f2_end_value;
      61         164 :   rhs(2) = df1dx_end_value;
      62         164 :   rhs(3) = df2dx_end_value;
      63             : 
      64         164 :   DenseVector<Real> coefs(4);
      65         164 :   mat.lu_solve(rhs, coefs);
      66             : 
      67         164 :   _A = coefs(0);
      68         164 :   _B = coefs(1);
      69         164 :   _C = coefs(2);
      70         164 :   _D = coefs(3);
      71             : 
      72         164 :   _initialized = true;
      73         164 : }
      74             : 
      75             : Real
      76        3242 : CubicTransition::value(const Real & x, const Real & f1, const Real & f2) const
      77             : {
      78             :   mooseAssert(_initialized, "initialize() must be called.");
      79             : 
      80        3242 :   if (x <= _x1)
      81         696 :     return f1;
      82        2546 :   else if (x >= _x2)
      83        1424 :     return f2;
      84             :   else
      85        1122 :     return _A * std::pow(x, 3) + _B * std::pow(x, 2) + _C * x + _D;
      86             : }
      87             : 
      88             : Real
      89        1270 : CubicTransition::derivative(const Real & x, const Real & df1dx, const Real & df2dx) const
      90             : {
      91             :   mooseAssert(_initialized, "initialize() must be called.");
      92             : 
      93        1270 :   if (x <= _x1)
      94         460 :     return df1dx;
      95         810 :   else if (x >= _x2)
      96         450 :     return df2dx;
      97             :   else
      98         360 :     return 3.0 * _A * std::pow(x, 2) + 2.0 * _B * x + _C;
      99             : }