1 /** @file gsBSplineBasis.h
2 
3     @brief Provides declaration of BSplineBasis class
4 
5     This file is part of the G+Smo library.
6 
7     This Source Code Form is subject to the terms of the Mozilla Public
8     License, v. 2.0. If a copy of the MPL was not distributed with this
9     file, You can obtain one at http://mozilla.org/MPL/2.0/.
10 
11     Author(s): A. Mantzaflaris, D. Mokris
12 */
13 
14 
15 #pragma once
16 
17 #include <gsCore/gsForwardDeclarations.h>
18 #include <gsCore/gsConstantBasis.h>
19 
20 #include <gsTensor/gsTensorBasis.h>
21 #include <gsTensor/gsTensorDomainIterator.h>
22 #include <gsTensor/gsTensorDomainBoundaryIterator.h>
23 
24 #include <gsNurbs/gsKnotVector.h>
25 
26 namespace gismo
27 {
28 
29 /// @brief Traits for BSplineBasis in more dimensions
30 template<short_t d, class T>
31 struct gsBSplineTraits
32 {
33     typedef gsKnotVector<T> KnotVectorType;
34 
35     typedef gsTensorBSplineBasis<d,T> Basis;
36     typedef gsTensorNurbsBasis<d,T>   RatBasis;
37     typedef gsTensorBSpline<d,T>      Geometry;
38     typedef gsTensorNurbs<d,T>        RatGeometry;
39 };
40 template<class T>
41 struct gsBSplineTraits<1,T>
42 {
43     typedef gsKnotVector<T> KnotVectorType;
44     typedef gsBSplineBasis<T>         Basis;
45     typedef gsNurbsBasis<T>           RatBasis;
46     typedef gsBSpline<T>              Geometry;
47     typedef gsNurbs<T>                RatGeometry;
48 };
49 template<class T>
50 struct gsBSplineTraits<0,T>
51 {
52     typedef gsKnotVector<T> KnotVectorType;
53     typedef gsConstantBasis<T>                       Basis;
54     typedef gsConstantBasis<T>                       RatBasis;
55     typedef gsConstantFunction<T>                    Geometry;
56     typedef gsConstantFunction<T>                    RatGeometry;
57 };
58 
59 /** \brief
60     A univariate B-spline basis.
61 
62     \tparam T coefficient type
63     \tparam KnotVectorType the type of knot vector to use
64 
65     \ingroup basis
66     \ingroup Nurbs
67 */
68 template<class T>
69 class gsTensorBSplineBasis<1,T> : public gsTensorBasis<1,T>
70 {
71 public:
72     typedef gsKnotVector<T> KnotVectorType;
73 
74     typedef gsTensorBasis<1,T> Base;
75 
76     typedef gsBSplineBasis<T> Self_t;
77 
78     typedef gsTensorBSplineBasis TensorSelf_t;
79     typedef gsTensorBSplineBasis<1,T> TensorSelf_tt;
80 
81     /// @brief Coefficient type
82     typedef T Scalar_t;
83 
84     /// @brief Associated geometry type
85     typedef typename gsBSplineTraits<1,T>::Geometry GeometryType;
86 
87     /// @brief Associated Boundary basis type
88     typedef typename gsBSplineTraits<0,T>::Basis BoundaryBasisType;
89 
90     /// @brief Dimension of the parameter domain
91     static const short_t Dim = 1;
92 
93     /// @brief Smart pointer for gsTensorBSplineBasis
94     typedef memory::shared_ptr< Self_t > Ptr;
95 
96     /// @brief Smart pointer for gsTensorBSplineBasis
97     typedef memory::unique_ptr< Self_t > uPtr;
98 
99 public:
100 
101     // Look at gsBasis class for a description
102     // Note: Specializing pointer type at return
103     //GISMO_UPTR_FUNCTION_PURE(TensorSelf_tt, clone)
104     private: virtual gsTensorBSplineBasis * clone_impl() const = 0;
105     public: uPtr clone() const { return uPtr(dynamic_cast<Self_t*>(clone_impl())); }
106 
107     // gsTensorBSplineBasis( const Base & o)
108     // {
109     //     const gsTensorBSplineBasis * a;
110     //     if ( ( a = dynamic_cast<const gsTensorBSplineBasis *>( &o )) )
111     //     {
112     //         m_p        = a->degree() ;
113     //         m_knots    = KnotVectorType( a->knots() );
114     //         m_periodic = a->m_periodic;
115     //     }
116     //     else
117     //         GISMO_ERROR("Cannot convert "<<o<<" to gsTensorBSplineBasis\n");
118     // }
119 
120     static Self_t * New(std::vector<gsBasis<T>*> & bb );
121 
122     static Self_t * New(std::vector<Self_t*> & bb )
123     {
124         return new Self_t(*bb.front());
125     }
126 
127     static uPtr make(std::vector<gsBasis<T>*> & bb )
128     {
129         return uPtr(New(bb));
130     }
131 
132     static uPtr make(std::vector<Self_t*> & bb )
133     {
134         return uPtr(new Self_t(*bb.front()));
135     }
136 
137     static uPtr make( const KnotVectorType & KV )
138     {
139         return uPtr(new Self_t(KV));
140     }
141 
142     // Note: these casts can be dangerous
143     // operator Self_t &() { return dynamic_cast<Self_t&>(*this);}
144     // operator const Self_t &() const { return dynamic_cast<const Self_t&>(*this);}
145 
146 public:
147 
148     void swap(gsTensorBSplineBasis& other)
149     {
150         std::swap(m_p, other.m_p);
151         std::swap(m_periodic, other.m_periodic);
152         m_knots.swap(other.m_knots);
153     }
154 
155 /* Virtual member functions required by the base class */
156 
157     // Look at gsBasis class for a description
158     short_t domainDim() const { return Dim; }
159 
160     // Unhide/forward gsTensorBasis<1,T>::size(k), since the overload
161     // with size() automatically hides it in this class
162     using Base::size;
163 
164     // Look at gsBasis class for a description
165     index_t size() const { return m_knots.size() - m_p - 1 - m_periodic; }
166 
167     // Look at gsBasis class for a description
168     size_t numElements() const { return m_knots.numElements(); }
169     using Base::numElements; //unhide
170 
171     // Look at gsBasis class for a description
172     size_t elementIndex(const gsVector<T> & u ) const;
173 
174     // Same as gsBasis::elementIndex but argument is a value instead of a vector
175     size_t elementIndex(T u ) const;
176 
177     // Look at gsBasis class for a description
178     gsMatrix<T> elementInSupportOf(index_t j) const;
179 
180     /// @brief Returns span (element) indices of the beginning and end
181     /// of the support of the i-th basis function.
182     void elementSupport_into(const index_t i, gsMatrix<index_t,1,2>& result) const
183     {
184         gsMatrix<index_t> tmp_vec;
185         m_knots.supportIndex_into(i, tmp_vec);
186         result = tmp_vec.cwiseMax(0).cwiseMin(m_knots.numElements());
187     }
188 
189     template <int _Rows>
190     gsMatrix<T> elementDom(const gsMatrix<index_t,_Rows,2> & box) const
191     {
192         gsMatrix<T> rvo(1,2);
193         rvo.at(0) = m_knots.uValue(box.at(0));
194         rvo.at(1) = m_knots.uValue(box.at(1));
195         return rvo;
196     }
197 
198     // Look at gsBasis class for a description
199     const TensorSelf_t & component(short_t i) const = 0;
200 
201     // Look at gsBasis class for a description
202     TensorSelf_t & component(short_t i) = 0;
203 
204     /// @brief Returns the anchors (greville points) of the basis
205     void anchors_into(gsMatrix<T> & result) const
206     {
207         m_knots.greville_into(result);
208     }
209 
210     /// @brief Returns the anchors (greville points) of the basis
211     void anchor_into(index_t i, gsMatrix<T> & result) const
212     {
213         result.resize(1,1);
214         result(0,0) = m_knots.greville(i);
215     }
216 
217     // Look at gsBasis class for a description
218     void connectivity(const gsMatrix<T> & nodes,
219                       gsMesh<T> & mesh) const;
220 
221     // Look at gsBasis class for a description
222     void active_into(const gsMatrix<T> & u, gsMatrix<index_t>& result) const;
223 
224     // Look at gsBasis class for a description
225     bool isActive(const index_t i, const gsVector<T> & u) const;
226 
227     // Look at gsBasis class for a description
228     gsMatrix<index_t> allBoundary( ) const ;
229 
230     // Look at gsBasis class for a description
231     gsMatrix<index_t> boundaryOffset(boxSide const & s,index_t offset) const;
232 
233 #ifdef __DOXYGEN__
234     /// @brief Gives back the boundary basis at boxSide s
235     typename BoundaryBasisType::uPtr boundaryBasis(boxSide const & s);
236 #endif
237     GISMO_UPTR_FUNCTION_DEC(gsConstantBasis<T>, boundaryBasis, boxSide const &)
238 
239     // Look at gsBasis class for a description
240     gsMatrix<T> support() const ;
241 
242     // Look at gsBasis class for a description
243     gsMatrix<T> support(const index_t & i ) const ;
244 
245     /// @brief Only meaningfull for periodic basis: For basis members that have
246     /// a twin, this function returns the other twin index, otherwise it
247     /// returns the same index as the argument
248     index_t twin(index_t i) const ;
249 
250     // Look at gsBasis class for a description
251     virtual void eval_into(const gsMatrix<T> & u, gsMatrix<T>& result) const;
252 
253     // Look at gsBasis class for a description
254     virtual void evalSingle_into(index_t i, const gsMatrix<T> & u, gsMatrix<T>& result) const;
255 
256     // Look at gsBasis class for a description
257     virtual void evalFunc_into(const gsMatrix<T> & u, const gsMatrix<T> & coefs, gsMatrix<T>& result) const;
258 
259     // Look at gsBasis class for a description
260     void deriv_into(const gsMatrix<T> & u, gsMatrix<T>& result ) const ;
261 
262     // Look at gsBasis class for a description
263     void derivSingle_into(index_t i, const gsMatrix<T> & u, gsMatrix<T>& result ) const ;
264 
265     // Look at gsBasis class for a description
266     void deriv_into(const gsMatrix<T> & u, const gsMatrix<T> & coefs, gsMatrix<T>& result ) const ;
267 
268     // Look at gsBasis class for a description
269     void deriv2_into(const gsMatrix<T> & u, gsMatrix<T>& result ) const ;
270 
271     // Look at gsBasis class for a description
272     void deriv2Single_into(index_t i, const gsMatrix<T> & u, gsMatrix<T>& result ) const ;
273 
274     // Look at gsBasis class for a description
275     void deriv2_into(const gsMatrix<T> & u, const gsMatrix<T> & coefs, gsMatrix<T>& result ) const ;
276 
277     // Look at gsBasis class for a description
278     gsMatrix<T> laplacian(const gsMatrix<T> & u ) const ;
279 
280     // Look at gsBasis class for a description
281     typename gsBasis<T>::uPtr tensorize(const gsBasis<T> & other) const;
282 
283     /// @brief Check the BSplineBasis for consistency
284     bool check() const
285     {
286         if ( m_periodic > 0 )
287         {
288             // Periodicity check wrt knot values
289             return (
290                 m_knots.degree()      == m_p &&
291                 (int)m_knots.size()   >  2*m_p+1
292                 );
293         }
294         else
295         {
296             return (
297                 m_knots.degree()      == m_p &&
298                 (int)m_knots.size()   >  2*m_p+1
299                 );
300         }
301     }
302 
303     /// @brief Prints the object as a string.
304     std::ostream &print(std::ostream &os) const;
305 
306     /// @brief Return a string with detailed information on the basis.
307     std::string detail() const;
308 
309     // Look at gsBasis class for a description
310     virtual void evalDerSingle_into(index_t i, const gsMatrix<T> & u,
311                                     int n, gsMatrix<T>& result) const;
312 
313     // Look at gsBasis class for a description
314     virtual void evalAllDers_into(const gsMatrix<T> & u, int n,
315                                   std::vector<gsMatrix<T> >& result) const;
316 
317     // Look at gsBasis class for a description
318     virtual void evalAllDersSingle_into(index_t i, const gsMatrix<T> & u,
319                                         int n, gsMatrix<T>& result) const;
320 
321     // Look at gsBasis class for a description
322     short_t degree(short_t i) const
323     {
324         GISMO_UNUSED(i);
325         GISMO_ASSERT(i==0,"Asked for degree(i) in 1D basis.");
326         return m_p;
327     }
328 
329     short_t degree() const {return m_p;}
330 
331     // Look at gsBasis class for a description
332     short_t maxDegree()   const { return m_p; }
333 
334     // Look at gsBasis class for a description
335     short_t minDegree()   const { return m_p; }
336 
337     // Look at gsBasis class for a description
338     short_t totalDegree() const { return m_p; }
339 
340     /// @brief Returns the order of the B-spline  basis
341     inline unsigned order() const { return m_p+1; }
342 
343     /// @brief True iff the point \a pp is in the domain of the basis
344     inline bool inDomain(T const & pp) const
345     { return ( (pp >= *(m_knots.begin()+m_p)) &&  (pp <= *(m_knots.end()-m_p-1) ) ); }
346 
347     /// @brief Returns the starting value of the domain of the basis
348     T domainStart() const { return *(m_knots.begin()+m_p); }
349 
350     /// @brief Returns the ending value of the domain of the basis
351     T domainEnd() const { return *(m_knots.end()-m_p-1); }
352 
353     /// @brief Returns length of the ``active" part of the knot vector.
354     T _activeLength() const { return domainEnd() - domainStart(); }
355 
356     /// @brief Returns the index of the first active (ie. non-zero) basis function at point u
357     /// Takes into account non-clamped knots.
358     inline index_t firstActive(T u) const {
359         return ( inDomain(u) ? (m_knots.iFind(u)-m_knots.begin()) - m_p : 0 );
360     }
361 
362     // Number of active functions at any point of the domain
363     inline index_t numActive() const { return m_p + 1; }
364 
365     // Look at gsBasis class for a description
366     gsDomain<T> * domain() const { return const_cast<KnotVectorType *>(&m_knots); }
367 
368     /// @brief Returns the knot vector of the basis
369     const KnotVectorType & knots (int i  = 0) const
370     {
371         GISMO_ENSURE(i==0, "Invalid knots requested");
372         return m_knots;
373     }
374     KnotVectorType & knots (int i  = 0)
375     {
376         GISMO_ENSURE(i==0, "Invalid knots requested");
377         return m_knots;
378     }
379 
380     T knot(index_t const & i) const { return m_knots[i];}
381 
382     /// @brief Inserts the knot \em knot in the underlying knot vector.
383     void insertKnot(T knot, int mult=1)
384     { m_knots.insert( knot, mult); }
385 
386     // compatibility with tensor-bsplines
387     void insertKnots(const std::vector< std::vector<T> >& refineKnots)
388     {
389         GISMO_ASSERT( refineKnots.size() == 1, "refineKnots vector has wrong size" );
390         this->knots().insert(refineKnots.front());
391     }
392 
393     // Look at gsBasis class for a description
394     void refineElements(std::vector<index_t> const & elements)
395     { m_knots.refineSpans(elements); }
396 
397     // Look at gsBasis class for a description
398     void uniformRefine(int numKnots = 1, int mul=1)
399     { m_knots.uniformRefine(numKnots,mul); }
400 
401     // Look at gsBasis class for a description
402     void uniformRefine_withCoefs(gsMatrix<T>& coefs, int numKnots = 1, int mul=1);
403 
404     // Look at gsBasis class for a description
405     void uniformRefine_withTransfer(gsSparseMatrix<T,RowMajor> & transfer, int numKnots = 1, int mul=1);
406 
407     // Look at gsBasis class for a description
408     void uniformCoarsen(int numKnots = 1)
409     { m_knots.coarsen(numKnots); }
410 
411     // Look at gsBasis class for a description
412     void uniformCoarsen_withTransfer(gsSparseMatrix<T,RowMajor> & transfer, int numKnots = 1);
413 
414     /// @brief Refine the basis by inserting the given knots and perform knot
415     /// refinement for the given coefficient matrix.
416     void refine_withCoefs(gsMatrix<T>& coefs, const std::vector<T>& knots);
417 
418     // compatibility with tensor-bsplines
419     void refine_withCoefs(gsMatrix<T> & coefs,
420                           const std::vector< std::vector<T> >& refineKnots)
421     {
422         refine_withCoefs(coefs,refineKnots.front() );
423     }
424 
425     /// @brief Refine the basis by inserting the given knots and produce a sparse matrix which maps coarse coefficient vectors to refined ones.
426     void refine_withTransfer(gsSparseMatrix<T,RowMajor> & transfer, const std::vector<T>& knots);
427 
428     void refine_withTransfer(gsSparseMatrix<T,RowMajor> & transfer,
429                              const std::vector<std::vector<T> >& knots)
430     {
431         GISMO_ASSERT(knots.size()==1, "the knots you want to insert do not have the right dimension");
432         refine_withTransfer(transfer,knots.front());
433         //GISMO_NO_IMPLEMENTATION
434     }
435 
436     /// @brief Increases the degree without adjusting the smoothness at inner
437     /// knots, except from the knot values in \a knots (constrained
438     /// knots of initial geometry)
439     ///
440     /// This type of refinement is known as k-refinement.  Note that
441     /// this type of refinement is ment to be performed after one (or
442     /// more) h-refinement steps the parent mesh (\a other),
443     /// otherwise this function is equivalent to p-refinement of \a other.
444     ///
445     /// @param other parent/reference mesh determining the
446     /// smoothness at the inner knots.
447     /// @param i number of k-refinement steps to perform
448     ///
449     /// @remarks Not tested yet!
450     void refine_k(const TensorSelf_t & other, int const & i = 1)
451     {
452         GISMO_ASSERT( m_p >= other.m_p, "Degree of other knot-vector should be lower.");
453         //for (typename std::vector<T>::iterator it =
454         //         m_knots.begin(); it != m_knots.end(); ++it)
455         //    GISMO_ASSERT( has(*it), "Knot "<< *it<<" is not in the knot vector.");
456 
457         // grab unique knots
458         const std::vector<T> knots = other.m_knots.unique();
459 
460         // Increase the degree without adjusting any knot
461         m_p += i;
462         m_knots.set_degree(m_p);
463         // Adjust (reduce) smoothness to satisfy initial constraint knots
464         m_knots.insert(knots,i);
465     }
466 
467     /// @brief p-refinement (essentially degree elevation)
468     void refine_p(short_t const & i = 1)
469     { degreeElevate(i); }
470 
471     /// @brief Uniform h-refinement (placing \a i new knots inside each knot-span
472     void refine_h(short_t const & i = 1)
473     { uniformRefine(i); }
474 
475     /// @brief Elevate the degree of the basis and preserve the smoothness
476     void degreeElevate(short_t const & i = 1, short_t const dir = -1)
477     {
478         GISMO_UNUSED(dir);
479         GISMO_ASSERT( dir == -1 || dir == 0, "Invalid direction");
480         m_p+=i; m_knots.degreeElevate(i);
481     }
482 
483     // Look at gsBasis for documentation
484     void degreeReduce (short_t const & i = 1, short_t const dir = -1)
485     {
486         GISMO_UNUSED(dir);
487         GISMO_ASSERT( dir == -1 || dir == 0, "Invalid direction");
488         GISMO_ASSERT( i<=m_p, "Cannot reduce degree to negative");
489         m_p-=i; m_knots.degreeReduce(i);
490         //m_periodic =
491     }
492 
493     // Look at gsBasis for documentation
494     void degreeIncrease(short_t const & i = 1, short_t const dir = -1)
495     {
496         GISMO_UNUSED(dir);
497         GISMO_ASSERT( dir == -1 || dir == 0, "Invalid direction");
498         m_p+=i; m_knots.degreeIncrease(i);
499     }
500 
501     // Look at gsBasis for documentation
502     void degreeDecrease(short_t const & i = 1, short_t const dir = -1)
503     {
504         GISMO_UNUSED(dir);
505         GISMO_ASSERT( dir == -1 || dir == 0, "Invalid direction");
506         m_p-=i; m_knots.degreeDecrease(i);
507     }
508 
509     /// @brief Elevate spline continuity at interior knots by \a i
510     void elevateContinuity(int const & i = 1)
511     {
512         GISMO_ASSERT( i>=0 && ( m_knots.size()>static_cast<size_t>(2*(m_p+1)) || i<=m_p ),
513                       "Cannot achieve continuity less than C^{-1} at interior knots.");
514         m_knots.reduceMultiplicity(i);
515     }
516 
517     /// @brief Reduces spline continuity at interior knots by \a i
518     void reduceContinuity(int const & i = 1)
519     {
520         GISMO_ASSERT( i>=0 && ( m_knots.size()>static_cast<size_t>(2*(m_p+1)) || i<=m_p ),
521                       "Cannot achieve continuity less than C^{-1} at interior knots.");
522         // TODO check: max interior mult + i <= m_p+1
523         m_knots.increaseMultiplicity(i);
524     }
525 
526     /// @brief Tells, whether the basis is periodic.
527     bool isPeriodic() const
528     {
529         return (m_periodic > 0);
530     }
531 
532     // Look at gsBasis class for a description
533     index_t functionAtCorner(boxCorner const & c) const;
534 
535     /// @brief Returns number of functions crossing the boundary of the knot vector.
536     int numCrossingFunctions () const
537     {
538         return m_periodic;
539     }
540 
541     /// @brief Checks, if both endknots have multiplicity m_p + 1.
542     bool isClamped() const
543     {
544         if( m_knots[0] != m_knots[m_p])
545             return false;
546 
547         else if( m_knots[m_knots.size() - m_p -1] != m_knots[m_knots.size()-1])
548             return false;
549 
550         else
551             return true;
552     }
553 
554     /// @brief If flag is true, tries to convert the basis to periodic
555     /// (succeeds only if the knot vector is suitable).
556     void setPeriodic(bool flag = true)
557     {
558         if ( flag )
559             _convertToPeriodic();
560         else
561             m_periodic = 0;
562     }
563 
564     // Compatible with tensor B-spline basis
565     void setPeriodic(int dir)
566     {
567         GISMO_UNUSED(dir);
568         GISMO_ASSERT(dir==0, "Invalid direction");
569             _convertToPeriodic();
570     }
571 
572     /// @brief Returns the multiplicity of the first ``significant" knot
573     /// (i.e., the m_p+1st). If it is different from the multiplicity
574     /// of the corresponding knot at the end, returns zero.
575     int borderKnotMult() const;
576 
577     typename gsBasis<T>::domainIter makeDomainIterator() const
578     {
579         return typename gsBasis<T>::domainIter(new gsTensorDomainIterator<T,1>(*this));
580     }
581 
582     typename gsBasis<T>::domainIter makeDomainIterator(const boxSide & s) const
583     {
584         return ( s == boundary::none ?
585                  typename gsBasis<T>::domainIter(new gsTensorDomainIterator<T,1>(*this)) :
586                  typename gsBasis<T>::domainIter(new gsTensorDomainBoundaryIterator<T,1>(*this, s))
587                 );
588     }
589 
590     /// @brief Moves the knot vectors to enforce periodicity.
591     void enforceOuterKnotsPeriodic();
592 
593     // Look at gsBasis class for a description
594     void reverse() { m_knots.reverse(); }
595 
596     void matchWith(const boundaryInterface & bi,
597                    const gsBasis<T> & other,
598                    gsMatrix<index_t> & bndThis,
599                    gsMatrix<index_t> & bndOther) const;
600 
601 protected:
602 
603     /// @brief Tries to convert the basis into periodic
604     void _convertToPeriodic();
605 
606     /// @brief Adjusts endknots so that the knot vector can be made periodic.
607     void _stretchEndKnots();
608 
609 public:
610 
611     /// @brief Helper function for evaluation with periodic basis.
612     ///
613     /// @param coefs coefficients (control points, one per row) before
614     /// converting the basis into periodic.
615     ///
616     /// @return copy of coefs with the first m_periodic rows copied to
617     /// the last m_periodic rows.
618     gsMatrix<T> perCoefs( const gsMatrix<T>& coefs ) const
619     {
620         gsMatrix<T> per_coefs = coefs;
621         per_coefs.bottomRows( m_periodic ) = coefs.topRows( m_periodic );
622         return per_coefs;
623     }
624 
625     gsMatrix<T> perCoefs(const gsMatrix<T>& coefs, int dir) const
626     {
627         GISMO_UNUSED(dir);
628         GISMO_ASSERT(dir==0, "Error");
629         return perCoefs(coefs);
630     }
631 
632     /// @brief Helper function for transforming periodic coefficients
633     /// to full coefficients
634     void expandCoefs(gsMatrix<T> & coefs) const
635     {
636         const index_t sz = coefs.rows();
637         coefs.conservativeResize(sz+m_periodic, Eigen::NoChange);
638         coefs.bottomRows( m_periodic ) = coefs.topRows( m_periodic );
639     }
640 
641     /// @brief Helper function for transforming full coefficients to
642     /// periodic coefficients
643     void trimCoefs(gsMatrix<T> & coefs) const
644     {
645         const index_t sz = coefs.rows();
646         coefs.conservativeResize(sz-m_periodic, Eigen::NoChange);
647     }
648 
649     /// @brief Returns the size of the basis ignoring the bureaucratic way of
650     /// turning the basis into periodic.
651     int trueSize() const
652     { return this->size() + m_periodic; }
653 
654 // Data members
655 protected:
656 
657     /// @brief Degree
658     short_t m_p;
659 
660     /// @brief Knot vector
661     KnotVectorType m_knots;
662 
663     /// @brief Denotes whether the basis is periodic, ( 0 -- non-periodic, >0 -- number of ``crossing" functions)
664     int m_periodic;
665 
666     /*/// @brief Multiplicity of the p+1st knot from the beginning and from the end.
667       int m_bordKnotMulti;*/
668 
669 }; // class gsTensorBSplineBasis<1>
670 
671 
672 //Using C++11 alias:
673 // template<class T, class KnotVectorType>
674 // using gsBSplineBasis = gsTensorBSplineBasis<1,T>
675 
676 /** \brief
677     A univariate B-spline basis.
678 
679     \tparam T coefficient type
680     \tparam KnotVectorType the type of knot vector to use
681 
682     \ingroup basis
683     \ingroup Nurbs
684 */
685 template<class T>
686 class gsBSplineBasis : public gsTensorBSplineBasis<1,T>
687 {
688 public:
689     typedef memory::shared_ptr< gsBSplineBasis > Ptr;
690     typedef memory::unique_ptr< gsBSplineBasis > uPtr;
691 
692     typedef gsKnotVector<T> KnotVectorType;
693     typedef gsTensorBSplineBasis<1,T> Base;
694     typedef gsBSplineBasis<T> Self_t;
695 
696     /// @brief Associated geometry type
697     typedef typename gsBSplineTraits<1,T>::Geometry GeometryType;
698 
699     /// @brief Associated Boundary basis type
700     typedef typename gsBSplineTraits<0,T>::Basis BoundaryBasisType;
701 
702 public:
703 
704     // Default empty constructor
705     explicit gsBSplineBasis(const bool periodic = false )
706     {
707         m_p = 0;
708         m_knots.initClamped(0);
709         m_periodic = 0;
710 
711         if( periodic )
712             this->_convertToPeriodic();
713 
714         if( ! this->check() )
715             gsWarn << "Warning: Inconsistent "<< *this<< "\n";
716     }
717 
718     /// @brief Construct BSpline basis of a knot vector
719     explicit gsBSplineBasis(KnotVectorType KV, const bool periodic = false)
720     {
721         m_p        = KV.degree();
722         m_knots.swap(KV);
723         m_periodic = 0;
724 
725         if( periodic )
726             this->_convertToPeriodic();
727 
728         if( ! this->check() )
729             gsWarn << "Warning: Insconsistent "<< *this<< "\n";
730     }
731 
732     /// @brief Compatibility constructor with input an std::vector containing
733     /// a single knotvector
734     explicit gsBSplineBasis(std::vector<KnotVectorType> KV)
735     {
736         GISMO_ASSERT(1 == KV.size(), "Expecting a single knotvector." );
737 
738         m_p        = KV.front().degree();
739         m_knots.swap(KV.front());
740         m_periodic = 0;
741     }
742 
743     /// @brief Construct a BSpline basis
744     /// @param u0 starting parameter
745     /// @param u1 end parameter parameter
746     /// @param interior number of interior knots
747     /// @param degree degree of the spline space
748     /// @param mult_interior multiplicity at the interior knots
749     /// @param periodic specifies if basis is periodic or not
750     gsBSplineBasis(const T u0, const T u1, const unsigned interior,
751                    const int degree, const unsigned mult_interior=1,
752                    const bool periodic = false )
753     {
754         m_p = degree;
755         m_knots.initUniform(u0, u1, interior, m_p+1, mult_interior, m_p);
756         m_periodic = 0;
757 
758         if( periodic )
759             this->_convertToPeriodic();
760 
761         if( ! this->check() )
762             gsWarn << "Warning: Insconsistent "<< *this<< "\n";
763     }
764 
765 /*
766     /// @brief Copy Constructor
767     gsBSplineBasis( const gsBSplineBasis & o)
768     : Base(o)
769     { }
770 */
771 
772     // Look at gsBasis class for a description
773     GISMO_CLONE_FUNCTION(gsBSplineBasis)
774 
775     // Look at gsBasis class for a description
776     Self_t & component(short_t i);
777 
778     // Look at gsBasis class for a description
779     const Self_t & component(short_t i) const;
780 
781     memory::unique_ptr<gsGeometry<T> > makeGeometry( gsMatrix<T> coefs ) const;
782 
783 private:
784 
785     using Base::m_p;
786     using Base::m_knots;
787     using Base::m_periodic;
788 };
789 
790 
791 } // namespace gismo
792 
793 
794 // *****************************************************************
795 #ifndef GISMO_BUILD_LIB
796 #include GISMO_HPP_HEADER(gsBSplineBasis.hpp)
797 #else
798 #ifdef gsBSplineBasis_EXPORT
799 #include GISMO_HPP_HEADER(gsBSplineBasis.hpp) // needed on thebrain?
800 #undef  EXTERN_CLASS_TEMPLATE
801 #define EXTERN_CLASS_TEMPLATE CLASS_TEMPLATE_INST
802 #endif
803 namespace gismo
804 {
805 template<class T> const short_t gsTensorBSplineBasis<1,T>::Dim; //-O3 (SLE11) fix
806 EXTERN_CLASS_TEMPLATE gsTensorBSplineBasis<1,real_t>;
807 EXTERN_CLASS_TEMPLATE gsBSplineBasis<real_t>;
808 }
809 #endif
810 // *****************************************************************
811