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