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 "PorousFlowCubic.h"
      11             : #include "libmesh/utility.h"
      12             : 
      13             : namespace PorousFlowCubic
      14             : {
      15             : Real
      16        6860 : cubic(Real x, Real x0, Real y0, Real y0p, Real x1, Real y1, Real y1p)
      17             : {
      18             :   mooseAssert(x0 != x1, "PorousFlowCubic: x0 cannot equal x1");
      19        6860 :   const Real d = x1 - x0;
      20             :   const Real d2 = Utility::pow<2>(d);
      21        6860 :   const Real mean = 0.5 * (x1 + x0);
      22        6860 :   const Real sq3 = 0.5 * std::sqrt(3.0) * d;
      23        6860 :   const Real term1 = y0p * (x - x0) * Utility::pow<2>(x - x1) /
      24        6860 :                      d2; // term1(x0) = term1(x1) = term1'(x1) = 0, term1'(x0) = y0p
      25        6860 :   const Real term2 = y1p * (x - x1) * Utility::pow<2>(x - x0) /
      26        6860 :                      d2; // term2(x0) = term2(x1) = term2'(x0) = 0, term2'(x1) = y1p
      27        6860 :   const Real term3 = (x - mean - sq3) * (x - mean) * (x - mean + sq3);
      28             :   // term3' = (x - mean) * (x - mean + sq3) + (x - mean - sq3) * (x - mean + sq3) + (x - mean - sq3)
      29             :   // * (x - mean)
      30             :   //        = 3 (x - mean)^2 - sq3^2
      31             :   // note term3' = 0 when x = mean +/- sq3/sqrt(3) = 0.5 * (x1 + x0) +/- 0.5 * (x1 - x0) = {x1, x0}
      32        6860 :   const Real term3_x0 = (x0 - mean - sq3) * (x0 - mean) * (x0 - mean + sq3);
      33        6860 :   const Real term3_x1 = (x1 - mean - sq3) * (x1 - mean) * (x1 - mean + sq3);
      34        6860 :   return (y0 * (term3 - term3_x1) + y1 * (term3_x0 - term3)) / (term3_x0 - term3_x1) + term1 +
      35        6860 :          term2;
      36             : }
      37             : 
      38             : Real
      39        6840 : dcubic(Real x, Real x0, Real y0, Real y0p, Real x1, Real y1, Real y1p)
      40             : {
      41             :   mooseAssert(x0 != x1, "PorousFlowCubic: x0 cannot equal x1");
      42        6840 :   const Real d = x1 - x0;
      43             :   const Real d2 = Utility::pow<2>(d);
      44        6840 :   const Real mean = 0.5 * (x1 + x0);
      45        6840 :   const Real sq3 = 0.5 * std::sqrt(3.0) * d;
      46        6840 :   const Real term1_prime = y0p * (Utility::pow<2>(x - x1) + (x - x0) * 2 * (x - x1)) / d2;
      47        6840 :   const Real term2_prime = y1p * (Utility::pow<2>(x - x0) + (x - x1) * 2 * (x - x0)) / d2;
      48             :   const Real term3_prime =
      49        6840 :       3.0 * Utility::pow<2>(mean) - 6 * mean * x - 0.75 * d2 + 3.0 * Utility::pow<2>(x);
      50        6840 :   const Real term3_x0 = (x0 - mean - sq3) * (x0 - mean) * (x0 - mean + sq3);
      51        6840 :   const Real term3_x1 = (x1 - mean - sq3) * (x1 - mean) * (x1 - mean + sq3);
      52        6840 :   return (y0 - y1) * term3_prime / (term3_x0 - term3_x1) + term1_prime + term2_prime;
      53             : }
      54             : }