G+Smo  24.08.0
Geometry + Simulation Modules
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
gsBasis.hpp
Go to the documentation of this file.
1 
14 #pragma once
15 
16 #include <gsCore/gsBasisFun.h>
17 #include <gsCore/gsDomainIterator.h>
18 #include <gsCore/gsBoundary.h>
19 #include <gsCore/gsGeometry.h>
20 
21 namespace gismo
22 {
23 
24 template<class T>
25 gsBasis<T>::~gsBasis()
26 { }
27 
28 template<class T>
29 gsBasisFun<T> gsBasis<T>::function(index_t i) const
30 {
31  return gsBasisFun<T>(*this,i);
32 }
33 
34 
35 // Evaluates a linear combination of basis functions (default implementation)
36 template<class T>
37 void gsBasis<T>::evalFunc_into(const gsMatrix<T> &u,
38  const gsMatrix<T> & coefs,
39  gsMatrix<T>& result) const
40 {
41  gsMatrix<T> B ;
42  gsMatrix<index_t> actives;
43 
44  // compute function values
45  this->eval_into(u,B);
46  // compute active functions
47  this->active_into(u,actives);
48 
49  // compute result as linear combination of
50  // "coefs(actives)" and B
51  linearCombination_into( coefs, actives, B, result );
52 }
53 
54 
55 // Evaluates the Jacobian of the function given by coefs (default implementation)
56 // For each point, result contains a geomDim x parDim matrix block containing the Jacobian matrix
57 template<class T>
58 void gsBasis<T>::jacobianFunc_into(const gsMatrix<T> &u, const gsMatrix<T> & coefs, gsMatrix<T>& result) const
59 {
60  const index_t n = coefs.cols();
61  const index_t numPts = u.cols(); // at how many points to evaluate the gradients
62  const index_t pardim = this->dim();
63 
64  result.setZero( n, pardim * numPts );
65  gsMatrix<T> B;
66  gsMatrix<index_t> ind;
67 
68  this->deriv_into(u,B);
69  this->active_into(u,ind);
70  const index_t numAct=ind.rows();
71 
72  for (index_t p = 0; p < numPts; ++p) // p = point
73  for (index_t c=0; c<n; ++c ) // c = component
74  for ( index_t a=0; a< numAct ; ++a ) // a = active function
75  {
76  result.block(c,p*pardim, 1, pardim).noalias() +=
77  coefs(ind(a,p), c) * B.block(a*pardim, p, pardim, 1).transpose();
78  }
79 }
80 
81 // Evaluates the partial derivatives of the function given by coefs (default implementation)
82 // For each point, result contains one column with stacked gradients of the components
83 template<class T>
84 void gsBasis<T>::derivFunc_into(const gsMatrix<T> &u, const gsMatrix<T> & coefs, gsMatrix<T>& result) const
85 {
86  gsMatrix<T> B;
87  gsMatrix<index_t> actives;
88 
89  // compute first derivatives
90  this->deriv_into(u,B);
91  // compute active functions
92  this->active_into(u,actives);
93 
94  // compute result as linear combination of
95  // "coefs(actives)" and B
96  linearCombination_into( coefs, actives, B, result );
97 }
98 
99 // Evaluates the second derivatives of the function given by coefs (default implementation)
100 template<class T>
101 void gsBasis<T>::deriv2Func_into(const gsMatrix<T> & u,
102  const gsMatrix<T> & coefs,
103  gsMatrix<T>& result) const
104 {
105  gsMatrix<T> B;
106  gsMatrix<index_t> actives;
107 
108  // compute second derivatives
109  this->deriv2_into(u,B);
110  // compute active functions
111  this->active_into(u,actives);
112 
113  // compute result as linear combination of
114  // "coefs(actives)" and B
115  linearCombination_into( coefs, actives, B, result );
116 }
117 
118 template<class T>
119 void gsBasis<T>::evalAllDersFunc_into(const gsMatrix<T> &u,
120  const gsMatrix<T> & coefs,
121  const unsigned n,
122  std::vector< gsMatrix<T> >& result,
123  bool sameElement) const
124 {
125  // resize result so that it will hold
126  // function values and up to the n-th derivatives
127  result.resize(n+1);
128  if ( 0 == u.cols() ) return;
129 
130  // B will contain the derivatives up to order n
131  std::vector< gsMatrix<T> >B;
132  // actives will contain the indices of the basis functions
133  // which are active at the evaluation points
134  gsMatrix<index_t> actives;
135 
136  this->evalAllDers_into(u,n,B,sameElement);
137  if (sameElement)
138  this->active_into(u.col(0), actives);
139  else
140  this->active_into(u, actives);
141 
142  // for derivatives 0 to n, evaluate the function by linear combination
143  // of coefficients with the respective function values/derivatives
144  for( unsigned i = 0; i <= n; i++)
145  linearCombination_into( coefs, actives, B[i], result[i], sameElement);
146 }
147 
148 
149 template<class T>
150 void gsBasis<T>::linearCombination_into(const gsMatrix<T> & coefs,
151  const gsMatrix<index_t> & actives,
152  const gsMatrix<T> & values,
153  gsMatrix<T> & result, bool sameElement)
154 {
155  const index_t numPts = values.cols() ;
156  const index_t tarDim = coefs.cols() ;
157  const index_t stride = values.rows() / actives.rows();
158 
159  GISMO_ASSERT( actives.rows() * stride == values.rows(),
160  "Number of values "<<values.rows()<< " and actives "<<actives.rows()<<" does not fit together");
161 
162  result.resize( tarDim * stride, numPts );
163  result.setZero();
164 
165  if (sameElement)
166  {
167  for ( index_t i = 0; i < actives.rows(); ++i ) // for all nonzero basis functions
168  for ( index_t c = 0; c < tarDim; ++c ) // for all components of the geometry
169  result.middleRows( stride * c, stride).noalias() +=
170  coefs( actives.at(i), c) * values.middleRows(stride * i, stride);
171  }
172  else
173  {
174  for ( index_t pt = 0; pt < numPts; ++pt ) // For pt, i.e., for every column of u
175  for ( index_t i = 0; i < actives.rows(); ++i ) // for all nonzero basis functions
176  for ( index_t c = 0; c < tarDim; ++c ) // for all components of the geometry
177  result.block( stride * c, pt, stride, 1).noalias() +=
178  coefs( actives(i,pt), c) * values.block( stride * i, pt, stride, 1);
179  }
180 }
181 
182 
183 template<class T>
184 inline gsMatrix<T> gsBasis<T>::laplacian(const gsMatrix<T> & u ) const
185 {
186  gsMatrix<T> tmp;
187  this->deriv2_into(u,tmp);
188  return tmp.colwise().sum();
189 }
190 
191 template<class T> inline
192 gsSparseMatrix<T> gsBasis<T>::collocationMatrix(const gsMatrix<T> & u) const
193 {
194  gsSparseMatrix<T> result( u.cols(), this->size() );
195  gsMatrix<T> ev;
196  gsMatrix<index_t> act;
197 
198  eval_into (u.col(0), ev);
199  active_into(u.col(0), act);
200  result.reservePerColumn( act.rows() );
201  for (index_t i=0; i!=act.rows(); ++i)
202  result.insert(0, act.at(i) ) = ev.at(i);
203 
204  for (index_t k=1; k!=u.cols(); ++k)
205  {
206  eval_into (u.col(k), ev );
207  active_into(u.col(k), act);
208  for (index_t i=0; i!=act.rows(); ++i)
209  result.insert(k, act.at(i) ) = ev.at(i);
210  }
211 
212  result.makeCompressed();
213  return result;
214 }
215 
216 template<class T> inline
217 memory::unique_ptr<gsGeometry<T> > gsBasis<T>::interpolateData( gsMatrix<T> const& vals,
218  gsMatrix<T> const& pts) const
219 {
220  GISMO_ASSERT (dim() == pts.rows() , "Wrong dimension of the points("<<
221  pts.rows()<<", expected "<<dim() <<").");
222  GISMO_ASSERT (this->size() == pts.cols() , "Expecting as many points as the basis functions." );
223  GISMO_ASSERT (this->size() == vals.cols(), "Expecting as many values as the number of points." );
224 
225  gsSparseMatrix<T> Cmat = collocationMatrix(pts);
226  gsMatrix<T> x ( this->size(), vals.rows());
227 
228  // typename gsSparseSolver<T>::BiCGSTABIdentity solver( Cmat );
229  // typename gsSparseSolver<T>::BiCGSTABDiagonal solver( Cmat );
230  // typename gsSparseSolver<T>::QR solver( Cmat );
231  typename gsSparseSolver<T>::BiCGSTABILUT solver( Cmat );
232 
233  // Solves for many right hand side columns
234  x = solver.solve( vals.transpose() );
235 
236  // gsDebug <<"gs Interpolate error : " << solver.error() << std::"\n";
237  // gsDebug <<"gs Interpolate iters : " << solver.iterations() << std::"\n";
238  // gsDebug <<"intpl sol : " << x.transpose() << std::"\n";
239 
240  return makeGeometry( give(x) );
241 }
242 
243 template<class T> inline
244 memory::unique_ptr<gsGeometry<T> > gsBasis<T>::interpolateAtAnchors(gsMatrix<T> const & vals) const
245 {
246  GISMO_ASSERT (this->size() == vals.cols(),
247  "Expecting as many values as the number of basis functions." );
248  gsMatrix<T> pts;
249  anchors_into(pts);
250  return interpolateData(vals, pts);
251 }
252 
253 /*
254 template<class T> inline
255 gsGeometry<T> * gsBasis<T>::interpolateAtAnchors(gsFunction<T> const & func) const
256 {
257  gsMatrix<T> pts, vals;
258  anchors_into(pts);
259  func.eval_into(pts, vals);
260  return interpolateData(vals, pts);
261 }
262 
263 template<class T>
264 gsGeometry<T> * gsBasis<T>::projectL2(gsFunction<T> const & func) const
265 {
266 
267  return NULL;
268 }
269 */
270 
271 template<class T> inline
272 void gsBasis<T>::anchors_into(gsMatrix<T>&) const
273 { GISMO_NO_IMPLEMENTATION }
274 
275 template<class T> inline
276 void gsBasis<T>::anchor_into(index_t, gsMatrix<T>&) const
277 { GISMO_NO_IMPLEMENTATION }
278 
279 
280 template<class T> inline
281 void gsBasis<T>::connectivityAtAnchors(gsMesh<T> & mesh) const
282 {
283  gsMatrix<T> nodes = anchors();
284  nodes.transposeInPlace();// coefficient vectors have ctrl points at rows
285  connectivity(nodes, mesh);
286 }
287 
288 template<class T>
289 void gsBasis<T>::connectivity(const gsMatrix<T> &, gsMesh<T> &) const
290 { GISMO_NO_IMPLEMENTATION }
291 
292 template<class T>
293 void gsBasis<T>::active_into(const gsMatrix<T> &, gsMatrix<index_t>&) const
294 { GISMO_NO_IMPLEMENTATION }
295 
296 template <class T>
297 bool gsBasis<T>::isActive(const index_t, const gsVector<T> &) const
298 { GISMO_NO_IMPLEMENTATION }
299 
300 template<class T>
301 void gsBasis<T>::numActive_into(const gsMatrix<T> &, gsVector<index_t>&) const
302 { GISMO_NO_IMPLEMENTATION }
303 
304 template<class T>
305 void gsBasis<T>::activeCoefs_into(const gsVector<T> &, const gsMatrix<T> &,
306  gsMatrix<T>&) const
307 { GISMO_NO_IMPLEMENTATION }
308 
309 template<class T>
310 gsMatrix<index_t>
311 gsBasis<T>::allBoundary() const
312 { GISMO_NO_IMPLEMENTATION }
313 
314 template<class T>
315 gsMatrix<index_t>
316 gsBasis<T>::boundaryOffset(boxSide const &,index_t) const
317 { GISMO_NO_IMPLEMENTATION }
318 
319 template<class T>
320 index_t
321 gsBasis<T>::functionAtCorner(boxCorner const &) const
322 { GISMO_NO_IMPLEMENTATION }
323 
325 template<class T>
326 gsBasis<T>* gsBasis<T>::boundaryBasis_impl(boxSide const &) const
327 { GISMO_NO_IMPLEMENTATION }
329 
330 template<class T>
331 typename gsBasis<T>::uPtr gsBasis<T>::componentBasis(boxComponent b) const
332 {
333  GISMO_ASSERT( b.totalDim() == this->dim(), "The dimensions do not agree." );
334 
335  const short_t dim = this->dim();
336 
337  uPtr result;
338  short_t d=0;
339  for (short_t i=0; i<dim; ++i)
340  {
341  boxComponent::location loc = b.locationForDirection(i);
342  if (loc)
343  {
344  if (result)
345  result = result->boundaryBasis( boxSide(loc+2*d) );
346  else
347  result = this->boundaryBasis( boxSide(loc+2*d) );
348  }
349  else
350  ++d;
351  }
352 
353  if (!result)
354  result = clone();
355 
356  return result;
357 }
358 
359 template<class T>
360 typename gsBasis<T>::uPtr gsBasis<T>::componentBasis_withIndices(boxComponent b, gsMatrix<index_t>& indices, bool noBoundary) const
361 {
362  GISMO_ASSERT( b.totalDim() == this->dim(), "The dimensions do not agree." );
363  const short_t dim = this->dim();
364 
365  uPtr result;
366  short_t d=0;
367  for (short_t i=0; i<dim; ++i)
368  {
369  boxComponent::location loc = b.locationForDirection(i);
370  if (loc)
371  {
372  if (result)
373  {
374  gsMatrix<index_t> tmp = result->boundary( boxSide(loc+2*d) );
375  for (index_t j=0; j<tmp.size(); ++j)
376  tmp(j,0) = indices(tmp(j,0),0);
377  tmp.swap(indices);
378  result = result->boundaryBasis( boxSide(loc+2*d) );
379  }
380  else
381  {
382  indices = this->boundary( boxSide(loc+2*d) );
383  result = this->boundaryBasis( boxSide(loc+2*d) );
384  }
385  }
386  else
387  ++d;
388  }
389 
390  if (!result)
391  {
392  result = clone();
393  const index_t sz = this->size();
394  indices.resize(sz,1);
395  for (index_t i=0;i<sz;++i)
396  indices(i,0) = i;
397  }
398 
399  if (noBoundary && d > 0)
400  {
401 
402  gsMatrix<index_t> bdy_indices = result->allBoundary();
403 
404  const index_t indices_sz = indices.rows();
405  const index_t bdy_indices_sz = bdy_indices.rows();
406 
407  // Copy all entries from indices to indices_cleaned except
408  // those with indices in bdy_indices
409 
410  gsMatrix<index_t> indices_cleaned(indices_sz - bdy_indices_sz, 1);
411  index_t j = 0, t = 0;
412  for (index_t i = 0; i < indices_sz; ++i)
413  {
414  if (util::greater(i, bdy_indices(j, 0)) && j < bdy_indices_sz)
415  ++j;
416  if (util::less(i, bdy_indices(j, 0)) || j == bdy_indices_sz)
417  {
418  indices_cleaned(t, 0) = indices(i, 0);
419  ++t;
420  }
421  }
422  GISMO_ASSERT(t == indices_cleaned.rows(), "Internal error.");
423  indices.swap(indices_cleaned);
424  }
425 
426  return result;
427 
428 }
429 
430 template<class T>
431 gsMatrix<T> gsBasis<T>::support() const
432 { GISMO_NO_IMPLEMENTATION }
433 
434 template<class T>
435 gsMatrix<T> gsBasis<T>::support(const index_t &) const
436 { GISMO_NO_IMPLEMENTATION }
437 
438 template<class T>
439 gsMatrix<T> gsBasis<T>::supportInterval(index_t dir) const
440 { return support().row(dir); }
441 
442 template<class T>
443 void gsBasis<T>::eval_into(const gsMatrix<T> &, gsMatrix<T>&) const
444 { GISMO_NO_IMPLEMENTATION }
445 
446 template<class T>
447 void gsBasis<T>::evalSingle_into(index_t, const gsMatrix<T> &, gsMatrix<T>&) const
448 { GISMO_NO_IMPLEMENTATION }
449 
450 template<class T>
451 void gsBasis<T>::deriv_into(const gsMatrix<T> &, gsMatrix<T>&) const
452 { GISMO_NO_IMPLEMENTATION }
453 
454 template<class T>
455 void gsBasis<T>::derivSingle_into(index_t,
456  const gsMatrix<T> &,
457  gsMatrix<T>&) const
458 { GISMO_NO_IMPLEMENTATION }
459 
460 template<class T>
461 void gsBasis<T>::deriv2_into(const gsMatrix<T> &, gsMatrix<T>&) const
462 { GISMO_NO_IMPLEMENTATION }
463 
464 template<class T>
465 void gsBasis<T>::deriv2Single_into(index_t,
466  const gsMatrix<T> &,
467  gsMatrix<T>&) const
468 { GISMO_NO_IMPLEMENTATION }
469 
470 template<class T>
471 void gsBasis<T>::evalAllDersSingle_into(index_t, const gsMatrix<T> &,
472  int, gsMatrix<T>&) const
473 { GISMO_NO_IMPLEMENTATION }
474 
475 template<class T>
476 void gsBasis<T>::evalDerSingle_into(index_t, const
477  gsMatrix<T> &, int,
478  gsMatrix<T>&) const
479 { GISMO_NO_IMPLEMENTATION }
480 
481 
482 template<class T>
483 typename gsBasis<T>::uPtr gsBasis<T>::create() const
484 { GISMO_NO_IMPLEMENTATION }
485 
486 
487 template<class T>
488 typename gsBasis<T>::uPtr gsBasis<T>::tensorize(const gsBasis &) const
489 { GISMO_NO_IMPLEMENTATION }
490 
491 template<class T>
492 typename gsBasis<T>::domainIter
493 gsBasis<T>::makeDomainIterator() const
494 { GISMO_NO_IMPLEMENTATION }
495 
496 template<class T>
497 typename gsBasis<T>::domainIter
498 gsBasis<T>::makeDomainIterator(const boxSide &) const
499 { GISMO_NO_IMPLEMENTATION }
500 
501 template<class T>
502 size_t gsBasis<T>::numElements(boxSide const &) const
503 { GISMO_NO_IMPLEMENTATION }
504 
505 template<class T>
506 size_t gsBasis<T>::elementIndex(const gsVector<T> &) const
507 { GISMO_NO_IMPLEMENTATION }
508 
509 template<class T>
510 gsMatrix<T> gsBasis<T>::elementInSupportOf(index_t) const
511 { GISMO_NO_IMPLEMENTATION }
512 
513 template<class T>
514 const gsBasis<T>& gsBasis<T>::component(short_t) const
515 { GISMO_NO_IMPLEMENTATION }
516 
517 template<class T>
518 gsBasis<T>& gsBasis<T>::component(short_t i)
519 { return const_cast<gsBasis<T>&>(const_cast<const gsBasis<T>*>(this)->component(i));}
520 
521 template<class T>
522 std::vector<index_t> gsBasis<T>::asElements(gsMatrix<T> const &, int) const
523 { GISMO_NO_IMPLEMENTATION }
524 
525 template<class T>
526 std::vector<index_t> gsBasis<T>::asElementsUnrefine(gsMatrix<T> const &, int) const
527 { GISMO_NO_IMPLEMENTATION }
528 
529 template<class T>
530 void gsBasis<T>::refine(gsMatrix<T> const &, int)
531 { GISMO_NO_IMPLEMENTATION }
532 
533 template<class T>
534 void gsBasis<T>::unrefine(gsMatrix<T> const &, int)
535 { GISMO_NO_IMPLEMENTATION }
536 
537 template<class T>
538 void gsBasis<T>::refineElements(std::vector<index_t> const &)
539 { GISMO_NO_IMPLEMENTATION }
540 
541 template<class T>
542 void gsBasis<T>::unrefineElements(std::vector<index_t> const &)
543 { GISMO_NO_IMPLEMENTATION }
544 
545 template<class T>
546 void gsBasis<T>::refineElements_withCoefs(gsMatrix<T> &,std::vector<index_t> const &)
547 { GISMO_NO_IMPLEMENTATION }
548 
549 template<class T>
550 void gsBasis<T>::unrefineElements_withCoefs(gsMatrix<T> &,std::vector<index_t> const &)
551 { GISMO_NO_IMPLEMENTATION }
552 
553 template<class T>
554 void gsBasis<T>::uniformRefine(int, int, int)
555 { GISMO_NO_IMPLEMENTATION }
556 
557 template<class T>
558 void gsBasis<T>::uniformRefine_withCoefs(gsMatrix<T>& , int , int , int )
559 { GISMO_NO_IMPLEMENTATION }
560 
561 template<class T>
562 void gsBasis<T>::uniformRefine_withTransfer(gsSparseMatrix<T,RowMajor> &,
563  int, int)
564 { GISMO_NO_IMPLEMENTATION }
565 
566 template<class T>
567 void gsBasis<T>::uniformCoarsen(int)
568 { GISMO_NO_IMPLEMENTATION }
569 
570 template<class T>
571 void gsBasis<T>::uniformCoarsen_withCoefs(gsMatrix<T>& , int )
572 { GISMO_NO_IMPLEMENTATION }
573 
574 template<class T>
575 void gsBasis<T>::uniformCoarsen_withTransfer(gsSparseMatrix<T,RowMajor> &,
576  int)
577 { GISMO_NO_IMPLEMENTATION }
578 
579 template<class T>
580 void gsBasis<T>::degreeElevate(short_t const &, short_t const)
581 { GISMO_NO_IMPLEMENTATION }
582 
583 template<class T>
584 void gsBasis<T>::degreeReduce(short_t const &, short_t const)
585 { GISMO_NO_IMPLEMENTATION }
586 
587 template<class T>
588 void gsBasis<T>::degreeIncrease(short_t const &, short_t const)
589 { GISMO_NO_IMPLEMENTATION }
590 
591 template<class T>
592 void gsBasis<T>::degreeDecrease(short_t const &, short_t const)
593 { GISMO_NO_IMPLEMENTATION }
594 
595 template<class T>
596 void gsBasis<T>::setDegree(short_t const& i)
597 {
598  const short_t dm = this->dim();
599  for (short_t k = 0; k!=dm; ++k)
600  {
601  const short_t p = this->degree(k);
602  if ( i > p )
603  {
604  this->degreeElevate(i-p, k);
605  }
606  else if ( i < p )
607  {
608  this->degreeReduce(p-i, k);
609  }
610  }
611 }
612 
613 template<class T>
614 void gsBasis<T>::setDegreePreservingMultiplicity(short_t const& i)
615 {
616  for ( short_t d = 0; d < dim(); ++ d )
617  {
618  if ( i > degree(d) )
619  degreeIncrease(i-degree(d),d);
620  else if ( i < degree(d) )
621  degreeDecrease(-i+degree(d),d);
622  }
623 }
624 
625 template<class T>
626 void gsBasis<T>::elevateContinuity(int const &)
627 { GISMO_NO_IMPLEMENTATION }
628 
629 template<class T>
630 void gsBasis<T>::reduceContinuity(int const &)
631 { GISMO_NO_IMPLEMENTATION }
632 
633 template<class T>
634 gsDomain<T> * gsBasis<T>::domain() const
635 { GISMO_NO_IMPLEMENTATION }
636 
637 template<class T>
638 short_t gsBasis<T>::maxDegree() const
639 { GISMO_NO_IMPLEMENTATION }
640 
641 template<class T>
642 short_t gsBasis<T>::minDegree() const
643 { GISMO_NO_IMPLEMENTATION }
644 
645 template<class T>
646 short_t gsBasis<T>::totalDegree() const
647 { GISMO_NO_IMPLEMENTATION }
648 
649 template<class T>
650 short_t gsBasis<T>::degree(short_t) const
651 { GISMO_NO_IMPLEMENTATION }
652 
653 template<class T>
654 void gsBasis<T>::reverse()
655 { GISMO_NO_IMPLEMENTATION }
656 
657 template<class T>
658 void gsBasis<T>::matchWith(const boundaryInterface &, const gsBasis<T> &,
659  gsMatrix<index_t> &, gsMatrix<index_t> &, index_t) const
660 { GISMO_NO_IMPLEMENTATION }
661 
662 template<class T>
663 T gsBasis<T>::getMinCellLength() const
664 {
665  const domainIter it = this->makeDomainIterator();
666  T h = 0;
667  for (; it->good(); it->next() )
668  {
669  const T sz = it->getMinCellLength();
670  if ( sz < h || h == 0 ) h = sz;
671  }
672  return h;
673 }
674 
675 template<class T>
676 T gsBasis<T>::getMaxCellLength() const
677 {
678  const domainIter it = this->makeDomainIterator();
679  T h = 0;
680  for (; it->good(); it->next() )
681  {
682  const T sz = it->getMaxCellLength();
683  if ( sz > h ) h = sz;
684  }
685  return h;
686 }
687 
688 
689 template<class T> inline
690 std::vector<gsSparseMatrix<T> >
691 gsBasis<T>::collocationMatrixWithDeriv(const gsBasis<T> & b, const gsMatrix<T> & u)
692 {
693  int dim = b.domainDim();
694  std::vector<gsSparseMatrix<T>> result(dim+1, gsSparseMatrix<T>( u.cols(), b.size() ));
695  std::vector<gsMatrix<T>> ev;
696  gsMatrix<index_t> act;
697 
698  b.evalAllDers_into (u.col(0), 1, ev);
699  b.active_into(u.col(0), act);
700  result[0].reservePerColumn( act.rows() );
701  result[1].reservePerColumn( act.rows() );
702  if (dim==2)
703  result[2].reservePerColumn( act.rows() );
704  for (index_t i=0; i!=act.rows(); ++i)
705  {
706  result[0].insert(0, act.at(i) ) = ev[0].at(i);
707  result[1].insert(0, act.at(i) ) = ev[1].at(dim*i);
708  if (dim == 2)
709  result[2].insert(0, act.at(i) ) = ev[1].at(dim*i+1);
710  }
711  for (index_t k=1; k!=u.cols(); ++k)
712  {
713  b.evalAllDers_into (u.col(k), 1, ev );
714  b.active_into(u.col(k), act);
715  for (index_t i=0; i!=act.rows(); ++i)
716  {
717  result[0].insert(k, act.at(i) ) = ev[0].at(i);
718  result[1].insert(k, act.at(i) ) = ev[1].at(dim*i);
719  if (dim == 2)
720  result[2].insert(k, act.at(i) ) = ev[1].at(dim*i +1);
721  }
722  }
723 
724  result[0].makeCompressed();
725  result[1].makeCompressed();
726  if (dim == 2)
727  result[2].makeCompressed();
728  return result;
729 }
730 
731 
732 
733 
734 // gsBasis<T>::linearComb(active, evals, m_tmpCoefs, result);
735 // gsBasis<T>::jacobianFromGradients(active, grads, m_tmpCoefs, result);
736 
737 /*
738 template<class T>
739 void gsBasis<T>::linearComb(const gsMatrix<index_t> & actives,
740  const gsMatrix<T> & basisVals,
741  const gsMatrix<T> & coefs,
742  gsMatrix<T>& result )
743 {
744  // basisVals.rows()==1 (or else basisVals.rows() == coefs.cols() and .cwiseProd)
745  result.resize(coefs.cols(), basisVals.cols()) ;
746 
747  for ( index_t j=0; j!=basisVals.cols() ; j++ ) // for all basis function values
748  {
749  //todo grab result.col(j)
750  result.col(j) = basisVals(0,j) * coefs.row( actives(0,j) ) ;//transpose ?
751  for ( index_t i=1; i< actives.rows() ; i++ )
752  result.col(j) += basisVals(i,j) * coefs.row( actives(i,j) ) ;
753  }
754 }
755 */
756 
757 }; // namespace gismo
gismo::gsDomain
Class representing a domain. i.e. a collection of elements (triangles, rectangles, cubes, simplices.
Definition: gsDomain.h:31
gismo::boxComponent::totalDim
short_t totalDim() const
Dimension of the computational domain (the box itself)
Definition: gsBoundary.h:506
gismo::gsBasis::evalAllDersFunc_into
virtual void evalAllDersFunc_into(const gsMatrix< T > &u, const gsMatrix< T > &coefs, const unsigned n, std::vector< gsMatrix< T > > &result, bool sameElement=false) const
Evaluates all derivatives up to order n of the function described by coefs at points u...
Definition: gsBasis.hpp:119
gismo::gsBasis::evalAllDersSingle_into
virtual void evalAllDersSingle_into(index_t i, const gsMatrix< T > &u, int n, gsMatrix< T > &result) const
Evaluate the basis function i and its derivatives up to order n at points u into result.
Definition: gsBasis.hpp:471
gismo::gsBasis::getMaxCellLength
virtual T getMaxCellLength() const
Get the maximum mesh size, as expected for approximation error estimates.
Definition: gsBasis.hpp:676
GISMO_NO_IMPLEMENTATION
#define GISMO_NO_IMPLEMENTATION
Definition: gsDebug.h:129
gismo::gsBasis::create
virtual gsBasis::uPtr create() const
Create an empty basis of the derived type and return a pointer to it.
Definition: gsBasis.hpp:483
gismo::util::less
bool less(T1 t1, T2 t2)
Definition: gsTemplateTools.h:247
gismo::gsBasis::derivSingle_into
virtual void derivSingle_into(index_t i, const gsMatrix< T > &u, gsMatrix< T > &result) const
Evaluates the (partial) derivatives of the i-th basis function at points u into result.
Definition: gsBasis.hpp:455
gismo::gsBasis::degreeIncrease
virtual void degreeIncrease(short_t const &i=1, short_t const dir=-1)
Elevate the degree of the basis by the given amount, preserve knots multiplicity. ...
Definition: gsBasis.hpp:588
gismo::gsBasis::supportInterval
gsMatrix< T > supportInterval(index_t dir) const
Returns an interval that contains the parameter values in direction dir.
Definition: gsBasis.hpp:439
short_t
#define short_t
Definition: gsConfig.h:35
gismo::gsBasis::refineElements
virtual void refineElements(std::vector< index_t > const &boxes)
Refinement function, with different sytax for different basis.
Definition: gsBasis.hpp:538
gsBoundary.h
Provides structs and classes related to interfaces and boundaries.
gismo::gsFunctionSet::evalAllDers_into
virtual void evalAllDers_into(const gsMatrix< T > &u, int n, std::vector< gsMatrix< T > > &result, bool sameElement=false) const
Evaluate the nonzero functions and their derivatives up to order n at points u into result...
Definition: gsFunctionSet.hpp:81
gismo::gsBasis::connectivity
virtual void connectivity(const gsMatrix< T > &nodes, gsMesh< T > &mesh) const
Definition: gsBasis.hpp:289
gismo::gsBasis::degree
virtual short_t degree(short_t i) const
Degree with respect to the i-th variable. If the basis is a tensor product of (piecewise) polynomial ...
Definition: gsBasis.hpp:650
gismo::gsBasis::minDegree
virtual short_t minDegree() const
If the basis is of polynomial or piecewise polynomial type, then this function returns the minimum po...
Definition: gsBasis.hpp:642
gismo::boundary
Struct that defines the boundary sides and corners and types of a geometric object.
Definition: gsBoundary.h:55
gismo::gsBasis::elementInSupportOf
virtual gsMatrix< T > elementInSupportOf(index_t j) const
Returns (the coordinates of) an element in the support of basis function j.
Definition: gsBasis.hpp:510
gismo::gsBasis::evalDerSingle_into
virtual void evalDerSingle_into(index_t i, const gsMatrix< T > &u, int n, gsMatrix< T > &result) const
Evaluate the (partial) derivative(s) of order n the i-th basis function at points u into result...
Definition: gsBasis.hpp:476
gsBasisFun.h
Provides definition of the BasisFun class.
gismo::util::greater
bool greater(T1 t1, T2 t2)
Definition: gsTemplateTools.h:313
gismo::gsBasis::derivFunc_into
virtual void derivFunc_into(const gsMatrix< T > &u, const gsMatrix< T > &coefs, gsMatrix< T > &result) const
Evaluate the derivatives of the function described by coefs at points u.
Definition: gsBasis.hpp:84
gismo::gsBasis::jacobianFunc_into
virtual void jacobianFunc_into(const gsMatrix< T > &u, const gsMatrix< T > &coefs, gsMatrix< T > &result) const
Evaluate the Jacobian of the function described by coefs at points u. Jacobian matrices are stacked i...
Definition: gsBasis.hpp:58
gismo::gsFunctionSet::size
virtual index_t size() const
size
Definition: gsFunctionSet.h:578
gismo::give
S give(S &x)
Definition: gsMemory.h:266
gsGeometry.h
Provides declaration of Geometry abstract interface.
gismo::gsBasis::elementIndex
virtual size_t elementIndex(const gsVector< T > &u) const
Returns an index for the element which contains point u.
Definition: gsBasis.hpp:506
index_t
#define index_t
Definition: gsConfig.h:32
gismo::gsBasis::numElements
virtual size_t numElements(boxSide const &s=0) const
The number of elements on side s.
Definition: gsBasis.hpp:502
gismo::gsBasis::componentBasis_withIndices
virtual uPtr componentBasis_withIndices(boxComponent b, gsMatrix< index_t > &indices, bool noBoundary=true) const
Returns the basis that corresponds to the component.
Definition: gsBasis.hpp:360
gismo::gsBasis::component
virtual const gsBasis< T > & component(short_t i) const
For a tensor product basis, return the (const) 1-d basis for the i-th parameter component.
Definition: gsBasis.hpp:514
gismo::gsBasis::domain
virtual gsDomain< T > * domain() const
Definition: gsBasis.hpp:634
GISMO_ASSERT
#define GISMO_ASSERT(cond, message)
Definition: gsDebug.h:89
gismo::gsBasis::deriv2_into
virtual void deriv2_into(const gsMatrix< T > &u, gsMatrix< T > &result) const
Evaluate the second derivatives of all active basis function at points u.
Definition: gsBasis.hpp:461
gismo::gsBasis::uniformRefine
virtual void uniformRefine(int numKnots=1, int mul=1, int dir=-1)
Refine the basis uniformly by inserting numKnots new knots with multiplicity mul on each knot span...
Definition: gsBasis.hpp:554
gismo::gsBasis::getMinCellLength
virtual T getMinCellLength() const
Get the minimum mesh size, as expected for inverse inequalities.
Definition: gsBasis.hpp:663
gismo::boxComponent::locationForDirection
location locationForDirection(index_t direction) const
Definition: gsBoundary.cpp:154
gismo::gsMesh
Class Representing a triangle mesh with 3D vertices.
Definition: gsMesh.h:31
gismo::gsBasis::uniformCoarsen_withCoefs
virtual void uniformCoarsen_withCoefs(gsMatrix< T > &coefs, int numKnots=1)
Coarsen the basis uniformly.
Definition: gsBasis.hpp:571
gismo::gsBasis::interpolateData
memory::unique_ptr< gsGeometry< T > > interpolateData(gsMatrix< T > const &vals, gsMatrix< T > const &pts) const
Applies interpolation given the parameter values pts and values vals.
Definition: gsBasis.hpp:217
gismo::gsMatrix::at
T at(index_t i) const
Returns the i-th element of the vectorization of the matrix.
Definition: gsMatrix.h:211
gismo::gsBasis::degreeDecrease
virtual void degreeDecrease(short_t const &i=1, short_t const dir=-1)
Lower the degree of the basis by the given amount, preserving knots multiplicity. ...
Definition: gsBasis.hpp:592
gismo::gsBasis::uniformCoarsen
virtual void uniformCoarsen(int numKnots=1)
Coarsen the basis uniformly by removing groups of numKnots consecutive knots, each knot removed mul t...
Definition: gsBasis.hpp:567
gismo::gsBasis::componentBasis
virtual uPtr componentBasis(boxComponent b) const
Returns the basis that corresponds to the component.
Definition: gsBasis.hpp:331
gismo::gsBasis::anchors_into
virtual void anchors_into(gsMatrix< T > &result) const
Returns the anchor points that represent the members of the basis in result. There is exactly one anc...
Definition: gsBasis.hpp:272
gismo::gsBasis::activeCoefs_into
virtual void activeCoefs_into(const gsVector< T > &u, const gsMatrix< T > &coefs, gsMatrix< T > &result) const
Returns the matrix result of active coefficients at points u, each row being one coefficient. The order of the rows is the same as active_into and eval_into functions.
Definition: gsBasis.hpp:305
gismo::gsBasis::totalDegree
virtual short_t totalDegree() const
If the basis is of polynomial or piecewise polynomial type, then this function returns the total poly...
Definition: gsBasis.hpp:646
gismo::gsBasis::uniformCoarsen_withTransfer
virtual void uniformCoarsen_withTransfer(gsSparseMatrix< T, RowMajor > &transfer, int numKnots=1)
Coarsen the basis uniformly and produce a sparse matrix which maps coarse coefficient vectors to refi...
Definition: gsBasis.hpp:575
gismo::gsBasis::uniformRefine_withCoefs
virtual void uniformRefine_withCoefs(gsMatrix< T > &coefs, int numKnots=1, int mul=1, int dir=-1)
Refine the basis uniformly.
Definition: gsBasis.hpp:558
gismo::gsBasis::refineElements_withCoefs
virtual void refineElements_withCoefs(gsMatrix< T > &coefs, std::vector< index_t > const &boxes)
Refine basis and geometry coefficients to levels.
Definition: gsBasis.hpp:546
gismo::gsBasis::setDegreePreservingMultiplicity
void setDegreePreservingMultiplicity(short_t const &i)
Set the degree of the basis (either increase or decrecee) in order to have degree equal to i...
Definition: gsBasis.hpp:614
gsDomainIterator.h
Provides declaration of DomainIterator abstract interface.
gismo::gsBasis::connectivityAtAnchors
virtual void connectivityAtAnchors(gsMesh< T > &mesh) const
Definition: gsBasis.hpp:281
gismo::gsBasis::refine
virtual void refine(gsMatrix< T > const &boxes, int refExt=0)
Refine the basis on the area defined by the matrix boxes.
Definition: gsBasis.hpp:530
gismo::gsBasis::function
gsBasisFun< T > function(index_t i) const
Returns the i-th basis function as a gsFunction.
Definition: gsBasis.hpp:29
gismo::gsFunctionSet::domainDim
virtual short_t domainDim() const =0
Dimension of the (source) domain.
gismo::gsBasis::linearCombination_into
static void linearCombination_into(const gsMatrix< T > &coefs, const gsMatrix< index_t > &actives, const gsMatrix< T > &values, gsMatrix< T > &result, bool sameElement=false)
Computes the linear combination coefs * values( actives )
Definition: gsBasis.hpp:150
gismo::boxSide
Struct which represents a certain side of a box.
Definition: gsBoundary.h:84
gismo::gsBasis::evalSingle_into
virtual void evalSingle_into(index_t i, const gsMatrix< T > &u, gsMatrix< T > &result) const
Evaluate the i-th basis function at points u into result.
Definition: gsBasis.hpp:447
gismo::gsBasisFun
Represents an individual function in a function set, or a certain component of a vector-valued functi...
Definition: gsBasisFun.h:36
gismo::gsBasis::numActive_into
virtual void numActive_into(const gsMatrix< T > &u, gsVector< index_t > &result) const
Returns the number of active (nonzero) basis functions at points u in result.
Definition: gsBasis.hpp:301
gismo::gsBasis::eval_into
virtual void eval_into(const gsMatrix< T > &u, gsMatrix< T > &result) const
Evaluates nonzero basis functions at point u into result.
Definition: gsBasis.hpp:443
gismo::gsBasis::interpolateAtAnchors
virtual memory::unique_ptr< gsGeometry< T > > interpolateAtAnchors(gsMatrix< T > const &vals) const
Applies interpolation of values pts using the anchors as parameter points. May be reimplemented in de...
Definition: gsBasis.hpp:244
gismo::gsBasis::uPtr
memory::unique_ptr< gsBasis > uPtr
Unique pointer for gsBasis.
Definition: gsBasis.h:89
gismo::gsBasis::tensorize
virtual gsBasis::uPtr tensorize(const gsBasis &other) const
Return a tensor basis of this and other.
Definition: gsBasis.hpp:488
gismo::boxComponent::location
location
Represents a location.
Definition: gsBoundary.h:530
gismo::gsBasis::boundaryOffset
virtual gsMatrix< index_t > boundaryOffset(boxSide const &s, index_t offset) const
Definition: gsBasis.hpp:316
gismo::gsBasis::reduceContinuity
virtual void reduceContinuity(int const &i=1)
Reduces the continuity of the basis along element boundaries.
Definition: gsBasis.hpp:630
gismo::gsBasis::allBoundary
virtual gsMatrix< index_t > allBoundary() const
Returns the indices of the basis functions that are nonzero at the domain boundary.
Definition: gsBasis.hpp:311
gismo::gsBasis::matchWith
virtual void matchWith(const boundaryInterface &bi, const gsBasis< T > &other, gsMatrix< index_t > &bndThis, gsMatrix< index_t > &bndOther, index_t offset=0) const
Computes the indices of DoFs that match on the interface bi. The interface is assumed to be a common ...
Definition: gsBasis.hpp:658
gismo::gsBasis::reverse
virtual void reverse()
Reverse the basis.
Definition: gsBasis.hpp:654
gismo::gsBasis::degreeElevate
virtual void degreeElevate(short_t const &i=1, short_t const dir=-1)
Elevate the degree of the basis by the given amount, preserve smoothness.
Definition: gsBasis.hpp:580
gismo::gsBasis::isActive
virtual bool isActive(const index_t i, const gsVector< T > &u) const
Returns true if there the point u with non-zero value or derivatives when evaluated at the basis func...
Definition: gsBasis.hpp:297
gismo::gsBasis::maxDegree
virtual short_t maxDegree() const
If the basis is of polynomial or piecewise polynomial type, then this function returns the maximum po...
Definition: gsBasis.hpp:638
gismo::boundaryInterface
Struct which represents an interface between two patches.
Definition: gsBoundary.h:649
gismo::gsBasis::evalFunc_into
virtual void evalFunc_into(const gsMatrix< T > &u, const gsMatrix< T > &coefs, gsMatrix< T > &result) const
Evaluate the function described by coefs at points u, i.e., evaluates a linear combination of coefs x...
Definition: gsBasis.hpp:37
gismo::gsBasis::makeDomainIterator
virtual domainIter makeDomainIterator() const
Create a domain iterator for the computational mesh of this basis, that points to the first element o...
Definition: gsBasis.hpp:493
gismo::gsBasis::uniformRefine_withTransfer
virtual void uniformRefine_withTransfer(gsSparseMatrix< T, RowMajor > &transfer, int numKnots=1, int mul=1)
Refine the basis uniformly.
Definition: gsBasis.hpp:562
gismo::gsBasis::deriv_into
virtual void deriv_into(const gsMatrix< T > &u, gsMatrix< T > &result) const
Evaluates the first partial derivatives of the nonzero basis function.
Definition: gsBasis.hpp:451
gismo::boxCorner
Struct which represents a certain corner of a hyper-cube.
Definition: gsBoundary.h:291
gismo::gsBasis::collocationMatrix
gsSparseMatrix< T > collocationMatrix(gsMatrix< T > const &u) const
Computes the collocation matrix w.r.t. points u.
Definition: gsBasis.hpp:192
gismo::gsBasis::anchor_into
virtual void anchor_into(index_t i, gsMatrix< T > &result) const
Returns the anchor point for member i of the basis.
Definition: gsBasis.hpp:276
gismo::gsBasis::support
virtual gsMatrix< T > support() const
Returns (a bounding box for) the domain of the whole basis.
Definition: gsBasis.hpp:431
gismo::gsBasis::laplacian
virtual gsMatrix< T > laplacian(const gsMatrix< T > &u) const
Compute the Laplacian of all nonzero basis functions at points u.
Definition: gsBasis.hpp:184
gismo::gsBasis
A basis represents a family of scalar basis functions defined over a common parameter domain...
Definition: gsBasis.h:78
gismo::gsBasis::setDegree
void setDegree(short_t const &i)
Set the degree of the basis (either elevate or reduce) in order to have degree equal to i wrt to each...
Definition: gsBasis.hpp:596
gismo::gsBasis::deriv2Single_into
virtual void deriv2Single_into(index_t i, const gsMatrix< T > &u, gsMatrix< T > &result) const
Evaluate the (partial) derivatives of the i-th basis function at points u into result.
Definition: gsBasis.hpp:465
gismo::gsBasis::deriv2Func_into
virtual void deriv2Func_into(const gsMatrix< T > &u, const gsMatrix< T > &coefs, gsMatrix< T > &result) const
Evaluates the second derivatives of the function described by coefs at points u.
Definition: gsBasis.hpp:101
gismo::gsBasis::elevateContinuity
virtual void elevateContinuity(int const &i=1)
Elevates the continuity of the basis along element boundaries.
Definition: gsBasis.hpp:626
gismo::gsBasis::active_into
virtual void active_into(const gsMatrix< T > &u, gsMatrix< index_t > &result) const
Returns the indices of active basis functions at points u, as a list of indices, in result...
Definition: gsBasis.hpp:293
gismo::boxComponent
Struct which represents a certain component (interior, face, egde, corner).
Definition: gsBoundary.h:445
gismo::gsBasis::degreeReduce
virtual void degreeReduce(short_t const &i=1, short_t const dir=-1)
Reduce the degree of the basis by the given amount, preserve smoothness.
Definition: gsBasis.hpp:584