1 2 /** @file gsBSplineSolver.h 3 4 @brief Provides classes and functions to solve equations involving B-splines. 5 6 This file is part of the G+Smo library. 7 8 This Source Code Form is subject to the terms of the Mozilla Public 9 License, v. 2.0. If a copy of the MPL was not distributed with this 10 file, You can obtain one at http://mozilla.org/MPL/2.0/. 11 12 Author(s): A. Mantzaflaris, A. Bressan 13 */ 14 15 16 17 #pragma once 18 19 #include <gsCore/gsLinearAlgebra.h> 20 #include <gsCore/gsForwardDeclarations.h> 21 22 namespace gismo { 23 24 25 enum RootType 26 { 27 odd, 28 even, 29 odd_interval, 30 even_interval, 31 invalid 32 }; 33 34 35 36 37 template <typename T=real_t> 38 struct Root 39 { 40 // requires cleanup 41 RootType type; 42 T parameter; 43 T begParameter; 44 T endParameter; 45 gsVector<T> point; 46 gsVector<T> begPoint; 47 gsVector<T> endPoint; 48 RootRoot49 Root(bool is_odd, T par, const gsVector<T> &mpoint) 50 : type(is_odd?odd:even), begParameter(0), endParameter(0) 51 { 52 parameter=par; 53 point=mpoint; 54 } RootRoot55 Root(bool is_odd, T parB, T parE, const gsVector<T> &pointB, const gsVector<T> &pointE) 56 : type(is_odd?odd_interval:even_interval), parameter(0) 57 { 58 begParameter=parB; 59 endParameter=parE; 60 begPoint=pointB; 61 endPoint=pointE; 62 } RootRoot63 Root() : type(invalid){} 64 }; 65 66 67 /** 68 \brief find intersections of a BSpline curve with an hyperplane 69 70 This function tries to be robust and to report correctly intersections 71 of the curve with the given hyperplane of thickness 2*tolerance. 72 Intersections are roots of a B-Spline curve and they can be of four types: 73 74 -odd points, i.e. the position with respect to the hyperplane changes 75 at the intersection 76 77 -even points, i.e. the position with respect to the hyperplane is the same 78 before and after the intersection 79 80 -odd intervals, i.e. the curve stays for a full parametric interval in the 81 hyperplane and then exit on the other side of it 82 83 -even intervals, i.e. the curve stays for a full parametric interval in the 84 hyperplane and then exit on the side it comes from 85 86 The intersections are reported as encountered while following the curve 87 in the direction given by increasing parameter. 88 89 This function assumes an open knot vector: the first and last control points 90 are in the curve. 91 If the start or the end of the the curve is on the hyperplane it is reported 92 as an odd intersection. This means that this special case must be handled 93 outside of this function when used to determine if a point is inside a 2D area 94 bounded by a closed curve. 95 96 \ingroup Nurbs 97 **/ 98 template <typename T> 99 unsigned findHyperPlaneIntersections ( 100 const gsBSpline<T> &curve, 101 const gsVector<T> &normal, 102 T reference, 103 T tolerance, 104 std::vector<Root<T> > &roots 105 ); 106 107 108 109 /** 110 A univariate root solver for B-spline curves. 111 */ 112 template<class T> 113 class gsBSplineSolver 114 { 115 typedef gsKnotVector<T> knotVectorType; 116 117 public: 118 /// Default empty constructor gsBSplineSolver()119 gsBSplineSolver() : m_n(1),m_k(1), m_d(1), eps(1e-7) { } 120 121 /// Destructor ~gsBSplineSolver()122 ~gsBSplineSolver() { } 123 124 public: 125 /// Return true if the first root exist, the value of the root is in this->value() 126 bool firstRoot(gsBSpline<T> const & bsp, int const & coord = 0, 127 T const & tr = 0, T const & tol = 1e-7, unsigned const & N=100) 128 { 129 initSolver(bsp,coord,tr,tol,N); 130 return nextRoot(); 131 } 132 133 /// Next root (requires that first root has been called before) 134 bool nextRoot (); 135 136 /// The value of the current root value()137 inline T value () { return x;} 138 139 // Return a vector with all the roots 140 void allRoots (gsBSpline<T> const & bsp, std::vector<T> & result, 141 int const & coord = 0, 142 T const & tr = 0, T const & tol = 1e-7, unsigned const & N=100) 143 { 144 result.clear(); 145 initSolver(bsp,coord,tr,tol,N); 146 147 while ( nextRoot() ) 148 { 149 result.push_back( value() ) ; 150 } 151 //gsLog<< "gsBSplineSolver: found "<< result.size() <<" roots.\n "; 152 } 153 154 private: 155 /// Initialize the solver with B-spline data initSolver(gsBSpline<T> const & bsp,int const & coord,T const & tr,T const & tol,unsigned const & N)156 void initSolver(gsBSpline<T> const & bsp , int const & coord, 157 T const & tr, T const & tol, unsigned const &N) 158 { 159 m_n = bsp.coefsSize(); 160 m_d = bsp.degree(); 161 m_k = 1; 162 eps = tol; 163 x = 0.0; 164 maxn = N + m_n; 165 166 const knotVectorType & kv = bsp.knots(); 167 168 m_t.resize(m_n+N+m_d+1); 169 m_c.resize(m_n+N); 170 171 for ( unsigned i = 0; i < m_n; ++i ) 172 { 173 m_c[i]= bsp.coef(i, coord) - tr; 174 m_t[i]= kv[i]; 175 } 176 177 for ( unsigned i = m_n; i < m_n + m_d +1 ; ++i ) 178 m_t[i]= kv[i]; 179 } 180 181 /// insert knot x in interval mu by Boehms algorithm 182 /// Note: t,c must have size at least n+1, n+d+2 respectively 183 /// require that x>=t[mu] 184 int insertKnot (int mu) ; 185 186 // Data members 187 private: 188 189 // m_t: knot vector, m_c: coefficients 190 std::vector<T> m_c, m_t; 191 //gsVector<T> m_c, m_t; 192 193 // m_n: coefficient size, maxn: maximum coeff. size // m_k span of root 194 unsigned m_n, maxn, m_k ; 195 196 int m_d; // m_d: degree 197 198 T eps; // tolerance 199 200 T x; // root value 201 202 }; // class gsClass 203 204 } //namespace gismo 205 206 #ifndef GISMO_BUILD_LIB 207 #include GISMO_HPP_HEADER(gsBSplineSolver.hpp) 208 #endif 209 210