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 "BSpline.h"
      11             : #include "libMeshReducedNamespace.h"
      12             : #include "MooseError.h"
      13             : #include "SplineUtils.h"
      14             : #include "Conversion.h"
      15             : 
      16             : using namespace libMesh;
      17             : 
      18             : namespace Moose
      19             : {
      20         272 : BSpline::BSpline(const unsigned int degree,
      21             :                  const libMesh::Point & start_point,
      22             :                  const libMesh::Point & end_point,
      23             :                  const libMesh::RealVectorValue & start_direction,
      24             :                  const libMesh::RealVectorValue & end_direction,
      25             :                  const unsigned int cps_per_half,
      26         272 :                  const libMesh::Real sharpness)
      27         272 :   : _degree(degree),
      28         272 :     _start_point(start_point),
      29         272 :     _end_point(end_point),
      30         272 :     _start_dir(start_direction),
      31         272 :     _end_dir(end_direction),
      32         272 :     _cps_per_half(cps_per_half),
      33         272 :     _sharpness(sharpness),
      34         272 :     _control_points(createControlPoints()),
      35         272 :     _knot_vector(buildKnotVector())
      36             : {
      37         272 : }
      38             : 
      39             : libMesh::Point
      40        5996 : BSpline::getPoint(const libMesh::Real t) const
      41             : {
      42             :   mooseAssert((t >= 0) && (t <= 1), "t should be in [0, 1]. t=" + std::to_string(t));
      43        5996 :   libMesh::Point returnPoint(0, 0, 0);
      44        5996 :   if (t == 1.0)
      45         272 :     return _control_points.back();
      46       49068 :   for (const auto i : index_range(_control_points))
      47       43344 :     returnPoint += BSpline::CdBBasis(t, i, _degree) * _control_points[i];
      48             : 
      49        5724 :   return returnPoint;
      50             : }
      51             : 
      52             : std::vector<libMesh::Point>
      53         272 : BSpline::createControlPoints() const
      54             : {
      55         272 :   std::vector<libMesh::Point> cps;
      56             :   /// call SplineUtils function
      57         272 :   cps = SplineUtils::bSplineControlPoints(
      58         272 :       _start_point, _end_point, _start_dir, _end_dir, _cps_per_half, _sharpness);
      59         272 :   return cps;
      60           0 : }
      61             : 
      62             : std::vector<libMesh::Real>
      63         272 : BSpline::buildKnotVector() const
      64             : {
      65         272 :   std::vector<libMesh::Real> knot_vector(_degree, 0); // initialize bottom to zeros
      66             : 
      67         272 :   const unsigned int num_points_left = _control_points.size() + 1 - _degree;
      68         272 :   if (_control_points.size() < _degree + 1)
      69           0 :     mooseError("Number of control points must be greater than or equal to degree + 1!");
      70             : 
      71        1730 :   for (const auto i : make_range(num_points_left))
      72             :   {
      73        1458 :     knot_vector.push_back(i); // count up
      74             :     mooseAssert(i < _control_points.size(),
      75             :                 "i must be less than the number of control points because of 0 indexing.");
      76             :   }
      77             : 
      78             :   // Filling with int(degree + 1) 0s at the beginning and 1s at the end makes the B-spline go
      79             :   // through the start and end point Note: the "+1" comes from the loop above which starts at 0 and
      80             :   // ends at num_points_left - 1
      81        1078 :   for (const auto i : make_range(_degree))
      82             :   {
      83         806 :     libmesh_ignore(i);
      84         806 :     knot_vector.push_back(num_points_left - 1);
      85             :   }
      86             : 
      87             :   // Normalize values s.t. the the max value is 1.
      88             :   // This must be done to have a constant domain of t to be between 0 and 1 (inclusive).
      89             :   // NOTE: the spacing is uniform, making this a uniform BSpline
      90        3342 :   for (auto & value : knot_vector)
      91             :   {
      92        3070 :     value /= (num_points_left - 1);
      93             :     mooseAssert(value <= 1.0 && value >= 0.0,
      94             :                 "knot_vector must be normalized to 1 with minimum value 0.");
      95             :   }
      96             : 
      97         272 :   return knot_vector;
      98           0 : }
      99             : 
     100             : libMesh::Real
     101      638448 : BSpline::CdBBasis(const libMesh::Real t, const unsigned int i, const unsigned int j) const
     102             : {
     103             :   mooseAssert(i + 1 < _knot_vector.size(), "Out of bounds access in knot vector");
     104             :   mooseAssert(_knot_vector[i] <= _knot_vector[i + 1],
     105             :               "Knot vector is ill-formatted. Ensure values are increasing." +
     106             :                   Moose::stringify(_knot_vector));
     107      638448 :   if (j == 0)
     108      340896 :     return ((t < _knot_vector[i + 1]) && (_knot_vector[i] <= t)) ? 1 : 0;
     109             :   else
     110      297552 :     return BSpline::firstCoeff(t, i, j) * BSpline::CdBBasis(t, i, j - 1) +
     111      297552 :            BSpline::secondCoeff(t, i, j) * BSpline::CdBBasis(t, i + 1, j - 1);
     112             : }
     113             : 
     114             : libMesh::Real
     115      297552 : BSpline::firstCoeff(const libMesh::Real t, const unsigned int i, const unsigned int j) const
     116             : {
     117             :   mooseAssert(i + j < _knot_vector.size(), "Out of bounds access in knot vector");
     118      297552 :   const auto ti = _knot_vector[i];
     119      297552 :   const auto tij = _knot_vector[i + j];
     120      297552 :   if (ti != tij)
     121      203124 :     return (t - ti) / (tij - ti);
     122             :   else
     123       94428 :     return 0;
     124             : }
     125             : 
     126             : libMesh::Real
     127      297552 : BSpline::secondCoeff(const libMesh::Real t, const unsigned int i, const unsigned int j) const
     128             : {
     129             :   mooseAssert(i + j + 1 < _knot_vector.size(), "Out of bounds access in knot vector");
     130      297552 :   libMesh::Real ti1 = _knot_vector[i + 1];
     131      297552 :   libMesh::Real tij1 = _knot_vector[i + j + 1];
     132      297552 :   if (ti1 != tij1)
     133      203124 :     return (tij1 - t) / (tij1 - ti1);
     134             :   else
     135       94428 :     return 0;
     136             : }
     137             : }