G+Smo  24.08.0
Geometry + Simulation Modules
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
gsBSplineAlgorithms.h
Go to the documentation of this file.
1 
14 #pragma once
15 
16 
17 namespace gismo
18 {
19 
27 namespace bspline
28 {
33  template <class T, typename KnotIterator, typename Derived>
34  void evalBasis( T u,
35  KnotIterator knot,
36  int deg,
37  gsEigen::MatrixBase<Derived> const & result )
38  {
39  STACK_ARRAY(T, left, deg + 1);
40  STACK_ARRAY(T, right, deg + 1);
41  evalBasis( u, knot, deg, left, right, result);
42  }
43 
44  template <class T, typename KnotIterator, typename Derived>
45  void evalBasis( T u,
46  KnotIterator knot,
47  int deg,
48  T left[], T right[],
49  gsEigen::MatrixBase<Derived> const & result )
50  {
51  gsEigen::MatrixBase<Derived>& res = const_cast<gsEigen::MatrixBase<Derived>&>(result);
52 
53  res(0,0)= (T)(1.0); // 0-th degree function value
54 
55  for(int j=1; j<= deg; j++) // For all degrees ( ndu column)
56  {
57  left[j] = u - *(knot+1-j);
58  right[j] = *(knot+j) - u;
59  T saved = (T)(0) ;
60  for(int r=0; r<j ; r++) // For all (except the last) basis functions of degree j ( ndu row)
61  {
62  const T temp = res(r,0) / ( right[r+1] + left[j-r] );
63  res(r,0) = saved + right[r+1] * temp ;// r-th function value of degree j
64  saved = left[j-r] * temp ;
65  }
66  res(j,0) = saved;
67  }
68  }
69 
71  template <class T, typename KnotIterator, typename Derived>
72  void evalDeg1Basis( const T & u,
73  const KnotIterator knot,
74  gsEigen::MatrixBase<Derived> const & result )
75  {
76  result(0,0) = u - (*knot+1);
77  result(1,0) = (*knot) - u ;
78  const T dn = (*knot) - (*knot+1) ;
79  if( dn != 0.0 )
80  result.array /= dn;
81  }
82 
84  template <class T, typename KnotIterator, typename Derived>
85  void evalDeg2Basis( const T & u,
86  const KnotIterator knot,
87  gsEigen::MatrixBase<Derived> const & result )
88  {
89 
90  }
91 
93  template <class T, typename KnotIterator, typename MatrixType>
94  void evalDeg3Basis( const T & u,
95  const KnotIterator knot,
96  MatrixType & result )
97  {
98 
99  }
100 
101 
106  template <class T, typename KnotIterator>
107  void deBoorTriangle( T u,
108  KnotIterator knot,
109  int deg,
110  T N[] )
111  {
112  const int table_size = deg+1;
113  T saved;
114 
115  for( int j = 0; j <= deg; j++ ) // Initialize zeroth-degree functions
116  if( u >= *(knot+j) && u < *(knot+j+1) )
117  N[ j*table_size ] = 1;
118  else
119  N[ j*table_size ] = 0;
120  for( int k = 1; k <= deg; k++ ) // Compute full triangular table
121  {
122  if( N[(k-1)] == 0 )
123  saved = 0;
124  else
125  saved = ( (u-*(knot) )*N[ k-1 ])/( *(knot+k)- *(knot));
126  for( int j = 0; j < deg-k+1; j++)
127  {
128  if( N[ (j+1)*table_size + (k-1) ] == 0 )
129  {
130  N[ j*table_size + k ] = saved;
131  saved = 0;
132  }
133  else
134  {
135  const T Uleft = *(knot+j+1);
136  const T Uright = *(knot+j+k+1);
137  const T temp = N[ (j+1)*table_size + (k-1) ]/(Uright-Uleft);
138  N[ j*table_size + k ] = saved + (Uright-u)*temp;
139  saved = (u-Uleft)*temp;
140  }
141  }
142  }
143  }
144 
145 
151  template <class T, typename KnotIterator, typename MatrixType>
152  void evalGeo( const T & u,
153  const KnotIterator & knot,
154  int deg,
155  MatrixType coefs, // makes a local copy of coefs
156  MatrixType & result )
157  {
158  result.resize(deg+1, 1 ) ;
159  T tmp;
160 
161  for ( int r=0; r!= deg; r++ ) // triangle step
162  //for ( int r=s; r< deg; r++ ) // TO DO: improve using multiplicity s
163  for ( int i=deg; i>r; i-- ) // recursive computation
164  {
165  tmp= ( u - *(knot+i) ) / ( *(knot+i+deg-r) - *(knot+i) );
166  coefs.row(i) = (T(1)-tmp)*coefs.row(i-1) + tmp* coefs.row(i) ;
167  }
168  result.noalias() = coefs.row(deg);
169  }
170 
175  template <class T, typename KnotIterator, typename MatrixType>
176  void evalBasis2( const T & u,
177  const KnotIterator & knot,
178  int deg,
179  MatrixType & result )
180  {
181  MatrixType pseudo_coefs;
182  pseudo_coefs.setIdentity(deg+1,deg+1);
183  evalGeoAlg(u, knot, deg, pseudo_coefs, result);
184  }
185 
186 
189  template <class T, typename KnotIterator, typename MatrixType>
191  {
192 
193 
194  }
195 
198  template <class T, typename KnotIterator, typename MatrixType>
199  void derivBasis( )
200  { }
201 
204  template <class T, typename KnotIterator, typename MatrixType>
206  { }
207 
210  template <class T, typename KnotIterator, typename MatrixType>
211  void allDersBasis( )
212  { }
213 
216  template <class T, typename KnotIterator, typename MatrixType>
217  void allDersSingle( )
218  { }
219 
220 
221 
223 template<class Basis_t>
224 void degreeElevateBSpline(Basis_t &basis,
226  short_t m)
227 {
228  typedef typename Basis_t::Scalar_t T;
229 
230  GISMO_ASSERT(m >= 0 && m<512, "Can only elevate degree by a positive (not enormous) amount.");
231  GISMO_ASSERT(basis.size() == coefs.rows(), "Invalid coefficients");
232 
233  if (m==0) return;
234 
235  const short_t p = basis.degree();
236  const index_t ncoefs = coefs.rows();
237  const index_t n = coefs.cols();
238  const gsKnotVector<T> & knots = basis.knots();
239 
240  // compute original derivative coefficients P (recurrence)
241  // original derivative coefficients
242 // # if defined(__GNUC__)
243 // STACK_ARRAY(gsMatrix<T>,P,p+1);
244 // # else // Note: No non-POD VLAs in clang/MSVC compiler
245  std::vector<gsMatrix<T> > P(p+1);
246 // # endif
247  for(short_t i=0;i<p+1;i++)
248  P[i].setZero(ncoefs - i, n);
249 
250  // insert first row
251  P[0].swap(coefs);
252  // fill table of derivative coefficients
253  for(short_t j=1; j<=p;j++)
254  for(index_t i=0; i<ncoefs-j; i++)
255  {
256  if(knots[i+p+1]>knots[i+j])
257  P[j].row(i).noalias() =
258  ( P[j-1].row(i+1) - P[j-1].row(i) ) / ( knots[i+p+1] - knots[i+j] );
259  }
260 
261  // loop over unique knot values (exept last one)
262  const typename gsKnotVector<T>::multContainer &
263  mult = knots.multiplicities(); // vector of mulitplicities
264 
265  // degree elevate basis
266  basis.degreeElevate(m);
267  const index_t ncoefs_new = basis.size();
268  const short_t p_new = basis.degree();
269 
270  // new (elevated) derivative coefficients
271 // # if defined(__GNUC__)
272 // STACK_ARRAY(gsMatrix<T>,Q,p_new+1);
273 // # else // Note: No non-POD VLAs in clang/MSVC compiler
274  std::vector<gsMatrix<T> > Q(p_new+1);
275 // # endif
276  for(short_t i=0; i<p_new+1; i++)
277  Q[i].setZero(ncoefs_new - i, n);
278 
279  // loop over knot intervals (with positive measure):
280  // prescibe known coefficients (see Huang Theorem 2) and fill
281  // corresponding triangular table using backwards recurrence
282 
283  // precompute factors
284  gsVector<T> factor = gsVector<T>::Ones(p+1);
285  for(short_t j=0; j<=p; j++)
286  for(short_t l=1; l<=j; l++)
287  factor[j] *= static_cast<T>(p+1-l) / static_cast<T>(p_new+1-l);
288 
289  //set known coefficients
290  // for(int j=0; j<=p; ++j)
291  // Q[j].row(0)=factor[j]*P[j].row(0);
292 
293  int betak = 0; // sum of interior mulitplicities
294  for(size_t k=0; k<mult.size()-1; k++)
295  {
296  //set known coefficients
297  for(short_t j=p+1-mult[k]; j<=p; ++j)
298  Q[j].row(betak+k*m) = factor[j] * P[j].row(betak);
299 
300  for(short_t j=p_new-1; j>=0;j--)
301  {
302  // fill triangular table
303  for(short_t i=1; i<=p_new-j; i++)
304  {
305  const int ik= i+betak+k*m; // update index i for the considered knot value
306  if(knots[ik+p_new]>knots[ik+j])
307  Q[j].row(ik).noalias() =
308  Q[j].row(ik-1) + Q[j+1].row(ik-1) * (knots[ik-1+p_new+1]-knots[ik-1+j+1]);
309 
310  }
311  }
312 
313  betak += mult[k+1];
314  }
315 
316  // Insert new coefficients into coefs
317  coefs.swap( Q[0] );
318 }
319 
320 
321 } // namespace bspline
322 
323 }// namespace gismo
void evalBasis2(const T &u, const KnotIterator &knot, int deg, MatrixType &result)
Definition: gsBSplineAlgorithms.h:176
void allDersSingle()
Definition: gsBSplineAlgorithms.h:217
multContainer multiplicities() const
Returns vector of multiplicities of the knots.
Definition: gsKnotVector.hpp:714
#define short_t
Definition: gsConfig.h:35
void evalDeg3Basis(const T &u, const KnotIterator knot, MatrixType &result)
Evaluation for degree 3 B-spline basis.
Definition: gsBSplineAlgorithms.h:94
void evalDeg2Basis(const T &u, const KnotIterator knot, gsEigen::MatrixBase< Derived > const &result)
Evaluation for degree 2 B-spline basis.
Definition: gsBSplineAlgorithms.h:85
size_t size() const
Number of knots (including repetitions).
Definition: gsKnotVector.h:242
#define index_t
Definition: gsConfig.h:32
A matrix with arbitrary coefficient type and fixed or dynamic size.
Definition: gsMatrix.h:38
#define GISMO_ASSERT(cond, message)
Definition: gsDebug.h:89
void degreeElevateBSpline(Basis_t &basis, gsMatrix< typename Basis_t::Scalar_t > &coefs, short_t m)
Increase the degree of a 1D B-spline from degree p to degree p + m.
Definition: gsBSplineAlgorithms.h:224
void allDersBasis()
Definition: gsBSplineAlgorithms.h:211
void evalGeo(const T &u, const KnotIterator &knot, int deg, MatrixType coefs, MatrixType &result)
Definition: gsBSplineAlgorithms.h:152
void evalBasis(T u, KnotIterator knot, int deg, gsEigen::MatrixBase< Derived > const &result)
Definition: gsBSplineAlgorithms.h:34
void deBoorTriangle(T u, KnotIterator knot, int deg, T N[])
Definition: gsBSplineAlgorithms.h:107
void evalBasisSingle()
Definition: gsBSplineAlgorithms.h:190
Class for representing a knot vector.
Definition: gsKnotVector.h:79
void derivBasis()
Definition: gsBSplineAlgorithms.h:199
void derivBasisSingle()
Definition: gsBSplineAlgorithms.h:205
void evalDeg1Basis(const T &u, const KnotIterator knot, gsEigen::MatrixBase< Derived > const &result)
Evaluation for degree 1 B-spline basis.
Definition: gsBSplineAlgorithms.h:72