G+Smo  24.08.0
Geometry + Simulation Modules
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
gsExpressions.h
Go to the documentation of this file.
1 
14 #pragma once
15 
16 #include <gsCore/gsFuncData.h>
19 
20 
21 namespace gismo
22 {
23 
24 // Adaptor to compute Hessian
25 template <typename Derived>
26 void secDerToHessian(const gsEigen::DenseBase<Derived> & secDers,
27  const index_t dim,
28  gsMatrix<typename Derived::Scalar> & hessian)
29 {
30  const index_t sz = dim*(dim+1)/2;
31  const gsAsConstMatrix<typename Derived::Scalar>
32  ders(secDers.derived().data(), sz, secDers.size() / sz );
33  hessian.resize(dim*dim, ders.cols() );
34 
35  switch ( dim )
36  {
37  case 1:
38  hessian = secDers.transpose(); //==ders
39  break;
40  case 2:
41  hessian.row(0)=ders.row(0);//0,0
42  hessian.row(1)=//1,0
43  hessian.row(2)=ders.row(2);//0,1
44  hessian.row(3)=ders.row(1);//1,1
45  break;
46  case 3:
47  hessian.row(0)=ders.row(0);//0,0
48  hessian.row(3)=//0,1
49  hessian.row(1)=ders.row(3);//1,0
50  hessian.row(6)=//0,2
51  hessian.row(2)=ders.row(4);//2,0
52  hessian.row(4)=ders.row(1);//1,1
53  hessian.row(7)=//1,2
54  hessian.row(5)=ders.row(5);//2,1
55  hessian.row(8)=ders.row(2);//2,2
56  break;
57  default:
58  break;
59  }
60 }
61 
63 template<class T> struct gsFeSpaceData
64 {
65  gsFeSpaceData(const gsFunctionSet<T> & _fs, index_t _dim, index_t _id):
66  fs(&_fs), dim(give(_dim)), id(give(_id)) { }
67 
68  const gsFunctionSet<T> * fs;
69  index_t dim, id;
70  gsDofMapper mapper;
71  gsMatrix<T> fixedDofs;
72  index_t cont; //int. coupling
73 
74  bool valid() const
75  {
76  GISMO_ASSERT(nullptr!=fs, "Invalid pointer.");
77  return static_cast<size_t>(fs->size()*dim)==mapper.mapSize();
78  }
79 
80  void init()
81  {
82  GISMO_ASSERT(nullptr!=fs, "Invalid pointer.");
83  if (const gsMultiBasis<T> * mb =
84  dynamic_cast<const gsMultiBasis<T>*>(fs) )
85  mapper = gsDofMapper(*mb, dim );
86  else if (const gsBasis<T> * b =
87  dynamic_cast<const gsBasis<T>*>(fs) )
88  mapper = gsDofMapper(*b, dim );
89  mapper.finalize();
90  fixedDofs.clear();
91  cont = -1;
92  }
93 };
94 
95 // Forward declaration in gismo namespace
96 template<class T> class gsExprHelper;
97 
105 namespace expr
106 {
107 
108 template <class E> struct is_arithmetic{enum{value=0};};
109 template <> struct is_arithmetic<real_t>{enum{value=1};};
110 template <typename E, bool = is_arithmetic<E>::value >
111 class _expr {using E::GISMO_ERROR_expr;};
112 
113 template<class T> class gsFeSpace;
114 template<class T> class gsFeVariable;
115 template<class T> class gsFeSolution;
116 template<class E> class symm_expr;
117 template<class E> class symmetrize_expr;
118 template<class E> class normalized_expr;
119 template<class E> class trace_expr;
120 template<class E> class integral_expr;
121 template<class E> class adjugate_expr;
122 template<class E> class norm_expr;
123 template<class E> class sqNorm_expr;
124 template<class E> class det_expr;
125 template<class E> class value_expr;
126 template<class E> class asdiag_expr;
127 template<class E> class max_expr;
128 template<class E> class rowsum_expr;
129 template<class E> class colsum_expr;
130 template<class E> class col_expr;
131 template<class T> class meas_expr;
132 template<class E> class inv_expr;
133 template<class E, bool cw = false> class tr_expr;
134 template<class E> class cb_expr;
135 template<class E> class abs_expr;
136 template<class E> class pow_expr;
137 template<class E> class sign_expr;
138 template<class E> class ppart_expr;
139 template<class E> class exp_expr;
140 template<class E> class ppartval_expr;
141 template<class T> class cdiam_expr;
142 template<class E> class temp_expr;
143 template<class E1, class E2, bool = E1::ColBlocks && !E1::ScalarValued && !E2::ScalarValued> class mult_expr
144 {using E1::GISMO_ERROR_mult_expr_has_invalid_template_arguments;};
145 
146 // Call as pow(a,b)
147 template<class E> pow_expr<E>
148 pow(_expr<E> const& u, real_t q) { return pow_expr<E>(u,q); }
149 
150 /*
151  Traits class for expressions
152 */
153 template <typename E> struct expr_traits
154 {
155 public:
156 // typedef typename E::Scalar Scalar;
157  typedef real_t Scalar;//todo
158  typedef const E Nested_t;
159 };
160 
161 # define Temporary_t typename util::conditional<ScalarValued,Scalar, \
162  typename gsMatrix<Scalar>::Base >::type
163 #if __cplusplus >= 201402L || _MSVC_LANG >= 201402L // c++14
164 # define MatExprType auto
165 # define AutoReturn_t auto
166 //note: in c++11 auto-return requires -> decltype(.)
167 #else // 199711L, 201103L
168 # define MatExprType typename gsMatrix<real_t>::constRef
169 # define AutoReturn_t typename util::conditional<ScalarValued,real_t,MatExprType>::type
170 #endif
171 
175 template <typename E>
176 class _expr<E, false>
177 {
178 protected://private:
179  _expr(){}
180  _expr(const _expr&) { }
181 public:
182  // Defined in derived classes: enum { Space, ScalarValued, ColBlocks }
183  // - ScalarValued: 0 is a scalar (must have Space=0),1 one denotes gsMatrix
184  // - ColBlocks: the expression stacks matrices per basis function
185  // - Space: 0: not a trial nor a test object (eg. normal vector, force function)
186  // 1: a test object (essentially a right-hand side vector expression)
187  // 2: a trial object
188  // 3: a trial+trial object (essentially a matrix expression)
189 
190  typedef typename expr_traits<E>::Nested_t Nested_t;
191  typedef typename expr_traits<E>::Scalar Scalar;
192 
194  void print(std::ostream &os) const
195  {
196  //gsInfo<<"\n Space="<<E::Space<<", ScV="<<E::ScalarValued<<", ColBlocks="<<E::ColBlocks<<"\n";
197  static_cast<E const&>(*this).print(os);
198  os<<"\n";
199  /*
200  std::string tmp(__PRETTY_FUNCTION__);
201  tmp.erase(0,74);
202  tmp.erase(tmp.size()-42,42);
203  size_t pos = 0;
204  while((pos=tmp.find(", false",0))!=std::string::npos) tmp.erase(pos,7);
205  while((pos=tmp.find(", true",0))!=std::string::npos) tmp.erase(pos,6);
206  while((pos=tmp.find("gismo::expr::",0))!=std::string::npos) tmp.erase(pos,13);
207  while((pos=tmp.find("_expr",0))!=std::string::npos) tmp.erase(pos,5);
208  while((pos=tmp.find("<double>",0))!=std::string::npos) tmp.erase(pos,8);
209  // while((pos=tmp.find("<long double>",0))!=std::string::npos) tmp.erase(pos,13);
210  // while((pos=tmp.find("<float>",0))!=std::string::npos) tmp.erase(pos,7);
211  tmp.erase(std::remove_if(tmp.begin(),tmp.end(),::isspace),tmp.end());
212  os<<tmp<<"\n";
213  */
214  }
215 
216  std::ostream & printDetail(std::ostream &os) const
217  {
218  os << (isVectorTr() ? "VectorTr " :
219  (isVector() ? "Vector " :
220  (isMatrix() ? "Matrix " :
221  "Scalar ") ) )
222  <<"expression of size "<< rows() // bug: this might be invalid if unparsed
223  << " x "<<cols()<<"\n";
224  print(os);
225  return os;
226  }
227 
229  MatExprType eval(const index_t k) const
230  { return static_cast<E const&>(*this).eval(k); }
231 
233  tr_expr<E> tr() const
234  { return tr_expr<E,false>(static_cast<E const&>(*this)); }
235 
237  tr_expr<E,true> cwisetr() const
238  { return tr_expr<E,true>(static_cast<E const&>(*this)); }
239 
241  cb_expr<E> cb() const
242  { return cb_expr<E>(static_cast<E const&>(*this)); }
243 
245  sign_expr<E> sgn(Scalar tolerance=0) const
246  { return sign_expr<E>(static_cast<E const&>(*this), tolerance); }
247 
249  exp_expr<E> exp() const
250  { return exp_expr<E>(static_cast<E const&>(*this)); }
251 
254  { return ppart_expr<E>(static_cast<E const&>(*this)); }
255  ppartval_expr<E> ppartval() const
256  { return ppartval_expr<E>(static_cast<E const&>(*this)); }
257 
259  mult_expr<real_t, ppart_expr<mult_expr<double,E,false>> , false>
260  npart() const { return -1* ( -(*this) ).ppart() ; }
261 
263  temp_expr<E> temp() const
264  { return temp_expr<E>(static_cast<E const&>(*this)); }
265 
267  inv_expr<E> const inv() const
268  { return inv_expr<E>(static_cast<E const&>(*this)); }
269 
271  trace_expr<E> trace() const
272  { return trace_expr<E>(static_cast<E const&>(*this)); }
273 
275  adjugate_expr<E> adj() const
276  { return adjugate_expr<E>(static_cast<E const&>(*this)); }
277 
279  norm_expr<E> norm() const
280  { return norm_expr<E>(static_cast<E const&>(*this)); }
281 
283  normalized_expr<E> normalized() const
284  { return normalized_expr<E>(static_cast<E const&>(*this)); }
285 
287  det_expr<E> det() const
288  { return det_expr<E>(static_cast<E const&>(*this)); }
289 
291  sqNorm_expr<E> sqNorm() const
292  { return sqNorm_expr<E>(static_cast<E const&>(*this)); }
293 
295  mult_expr<E,E,0> (sqr)() const { return (*this)*(*this); }
296 
297  symm_expr<E> symm() const
298  { return symm_expr<E>(static_cast<E const&>(*this)); }
299 
300  symmetrize_expr<E> symmetrize() const
301  { return symmetrize_expr<E>(static_cast<E const&>(*this)); }
302 
305  value_expr<E> val() const
306  { return value_expr<E>(static_cast<E const&>(*this)); }
307 
310  { return asdiag_expr<E>(static_cast<E const&>(*this)); }
311 
313  max_expr<E> max() const
314  { return max_expr<E>(static_cast<E const&>(*this)); }
315 
317  rowsum_expr<E> rowSum() const
318  { return rowsum_expr<E>(static_cast<E const&>(*this)); }
319 
321  colsum_expr<E> colSum() const
322  { return colsum_expr<E>(static_cast<E const&>(*this)); }
323 
324  col_expr<E> operator[](const index_t i) const
325  { return col_expr<E>(static_cast<E const&>(*this),i); }
326 
328  index_t rows() const
329  { return static_cast<E const&>(*this).rows(); }
330 
332  index_t cols() const
333  { return static_cast<E const&>(*this).cols(); }
334 
335  index_t cardinality() const
336  { return static_cast<E const&>(*this).cardinality_impl(); }
337 
338  static index_t cardinality_impl() { return 1; }
339 
342  bool isScalar() const { return rows()*cols()<=1; }
343 
344  static bool isVector () { return 1==E::Space; }
345  static bool isVectorTr() { return 2==E::Space; }
346  static bool isMatrix () { return 3==E::Space; }
347 
350  void parse(gsExprHelper<Scalar> & evList) const
351  { static_cast<E const&>(*this).parse(evList); }
352 
353  template<class op> void apply(op & _op) const
354  { static_cast<E const&>(*this).apply(_op); }
355 
358  const gsFeSpace<Scalar> & rowVar() const
359  {
360  // assert ValueType!=0
361  return static_cast<E const&>(*this).rowVar();
362  }
363 
366  const gsFeSpace<Scalar> & colVar() const
367  {
368  // assert ValueType==2
369  return static_cast<E const&>(*this).colVar();
370  }
371 
372  // Overload conversions, eg. converts _expr<mult_expr> to
373  // mult_expr.
374  operator E&() { return static_cast< E&>(*this); }
375  operator E const&() const { return static_cast<const E&>(*this); }
376 
377  E const & derived() const { return static_cast<const E&>(*this); }
378 };
379 
381 template <typename E>
382 std::ostream &operator<<(std::ostream &os, const _expr<E> & b)
383 {b.print(os); return os; }
384 
385 }
386 }
387 
389 
390 namespace gismo
391 {
392 namespace expr
393 {
394 
395 /*
396  Null expression is a compatibility expression invalid at runtime
397 */
398 template<class T>
399 class gsNullExpr : public _expr<gsNullExpr<T> >
400 {
401 public:
402 
403  operator const gsFeSpace<T> & () const
404  {
405  static gsFeSpace<T> vv(-1);
406  return vv;
407  }
408 
409  typedef T Scalar;
410  gsMatrix<T> eval(const index_t) const { GISMO_ERROR("gsNullExpr"); }
411  inline index_t rows() const { GISMO_ERROR("gsNullExpr"); }
412  inline index_t cols() const { GISMO_ERROR("gsNullExpr"); }
413  void parse(gsExprHelper<T> &) const { }
414 
415  const gsFeSpace<T> & rowVar() const { GISMO_ERROR("gsNullExpr"); }
416  const gsFeSpace<T> & colVar() const { GISMO_ERROR("gsNullExpr"); }
417 
418  void print(std::ostream &os) const { os << "NullExpr"; }
419 
420  static const gsNullExpr & get()
421  {
422  static gsNullExpr o;
423  return o;
424  }
425 //private:
426  gsNullExpr() {}
427 };
428 
429 
430 template<class E>
431 class symbol_expr : public _expr<E>
432 {
433 public:
434  typedef typename expr_traits<E>::Scalar Scalar;
435 
436  friend class gismo::gsExprHelper<Scalar>;
437 protected:
438  const gsFunctionSet<Scalar> * m_fs;
439  const gsFuncData<Scalar> * m_fd;
440  index_t m_d;
441  bool m_isAcross;
442 
443 public:
444 
446  bool isAcross() const { return m_isAcross; }
447 
448  E right() const
449  {
450  E ac(this->derived());
451  ac.m_fs = m_fs;//needed?
452  ac.m_isAcross = true;
453  return ac;
454  }
455 
456  E left() const
457  {
458  E ac(this->derived());
459  ac.m_fs = m_fs;
460  ac.m_isAcross = false;
461  return ac;
462  }
463 
465  const gsFunctionSet<Scalar> & source() const {return *m_fs;}
466 
468  const gsFuncData<Scalar> & data() const
469  {
470  GISMO_ASSERT(NULL!=m_fd, "FuncData member not registered "<<this<<"/"<< m_fs);
471  return *m_fd;
472  }
473 
474  // used by FeSpace, FeVariable, ..
475  void parse(gsExprHelper<Scalar> & evList) const
476  {
477  evList.add(*this);
478  this->m_fd->flags |= NEED_VALUE | NEED_ACTIVE;
479  }
480 
481  index_t cardinality_impl() const
482  {
483  GISMO_ASSERT(this->data().actives.rows()!=0,"Cardinality depends on the NEED_ACTIVE flag");
484  return m_d * this->data().actives.rows();
485  }
486 
487  //public for now due to Bc
488  void setSource(const gsFunctionSet<Scalar> & fs) { m_fs = &fs;}
489 
490 private:
491  void setData(const gsFuncData<Scalar> & val) { m_fd = &val;}
492  void setDim(index_t _d) { m_d = _d; }
493  void clear() { m_fs = NULL; }
494 
495 protected:
496  explicit symbol_expr(index_t _d)
497  : m_fs(NULL), m_fd(NULL), m_d(_d), m_isAcross(false) { }
498 
499 public:
500  bool isValid() const { return NULL!=m_fd && NULL!=m_fs; }
501 
502  // component
503  // expr comp(const index_t i) const { return comp_expr<Scalar>(*this,i); }
504  // eval(k).col(i)
505 
506  // The evaluation return rows for (basis) functions and columns
507  // for (coordinate) components
508  MatExprType eval(const index_t k) const
509  { return m_fd->values[0].col(k).blockDiag(m_d); }
510  //{ return m_fd->values[0].col(k); }
511 
512  const gsFeSpace<Scalar> & rowVar() const {return gsNullExpr<Scalar>::get();}
513  const gsFeSpace<Scalar> & colVar() const {return gsNullExpr<Scalar>::get();}
514 
515  index_t rows() const
516  {
517  GISMO_ASSERT(NULL!=m_fs, "FeVariable: Function source not registered");
518  return m_fs->targetDim();
519  }
520 
521  index_t cols() const { return m_d; }
522 
523  void print(std::ostream &os) const { os << "u"; }
524 
525 public:
526 
528  index_t dim() const { return m_d;}
529 
532  index_t targetDim() const { return m_fs->targetDim(); }
533 
535  index_t parDim() const { return m_fs->domainDim(); }
536 
537  index_t cSize() const
538  {
539  GISMO_ASSERT(0!=m_fd->values[0].size(),"Probable error.");
540  return m_fd->values[0].rows();
541  } // coordinate size
542 };
543 
544 /*
545  Column expression
546 */
547 template<class E>
548 class col_expr : public _expr<col_expr<E> >
549 {
550  typename E::Nested_t _c;
551  const index_t _i;
552 public:
553  typedef typename E::Scalar Scalar;
554  typedef const col_expr<E> Nested_t;
555 
556  enum { Space = E::Space, ScalarValued = 0, ColBlocks = 0 };
557 
558  col_expr(const E & c, const index_t i) : _c(c), _i(i) { }
559 
560 public:
561 
562  //ConstColXpr
563  inline MatExprType eval(const index_t k) const { return _c.eval(k).col(_i); }
564 
565  index_t rows() const { return _c.rows(); }
566  index_t cols() const { return 1; }
567  void parse(gsExprHelper<Scalar> & evList) const { _c.parse(evList); }
568 
569  const gsFeSpace<Scalar> & rowVar() const { return _c.rowVar(); }
570  const gsFeSpace<Scalar> & colVar() const { return _c.colVar(); }
571 
572  void print(std::ostream &os) const { os<<_c<<"["<<_i<<"]"; }
573 };
574 
575 /*
576  Expression for a constant value
577 */
578 template<class T>
579 class _expr<T, true> : public _expr<_expr<T> >
580 {
581  const T _c;
582 public:
583  typedef T Scalar;
584  typedef const _expr<T> Nested_t;
585 
586  explicit _expr(Scalar c) : _c(give(c)) { }
587 
588 public:
589  enum {Space = 0, ScalarValued = 1, ColBlocks= 0};
590 
591  inline Scalar eval(const index_t ) const { return _c; }
592 
593  inline _expr val() const { return *this; }
594  index_t rows() const { return 0; }
595  index_t cols() const { return 0; }
596  void parse(gsExprHelper<Scalar> &) const { }
597 
598  const gsFeSpace<T> & rowVar() const { return gsNullExpr<T>::get(); }
599  const gsFeSpace<T> & colVar() const { return gsNullExpr<T>::get(); }
600 
601  void print(std::ostream &os) const { os<<_c; }
602 };
603 
604 /*
605  Geometry map expression
606 */
607 template<class T>
608 class gsGeometryMap : public _expr<gsGeometryMap<T> >
609 {
610  const gsFunctionSet<T> * m_fs;
611  const gsMapData<T> * m_fd;
612  //index_t d, n;
613 
614  bool m_isAcross;
615 
616 public:
617  enum {Space = 0, ScalarValued= 0, ColBlocks= 0};
618 
619  bool isAcross() const { return m_isAcross; }
620 
621  gsGeometryMap right() const
622  {
623  gsGeometryMap ac;
624  ac.m_fs = m_fs;
625  ac.m_isAcross = true;
626  return ac;
627  }
628 
629  gsGeometryMap left() const
630  {
631  gsGeometryMap ac;
632  ac.m_fs = m_fs;
633  ac.m_isAcross = false;
634  return ac;
635  }
636 
638  const gsFunctionSet<T> & source() const {return *m_fs;}
639 
641  const gsMapData<T> & data() const
642  {
643  GISMO_ASSERT(NULL!=m_fd, "gsGeometryMap: invalid data "<< m_fs <<","<<m_fd);
644  return *m_fd;
645  }
646 
647  index_t targetDim() const { return m_fs->targetDim();}
648  index_t domainDim() const { return m_fs->domainDim();}
649 
651  void copyCoefs( const gsGeometryMap<T> & other) const
652  {
653  GISMO_ASSERT( dynamic_cast<const gsMultiPatch<T>*>( this->m_fs ), "error");
654  const gsMultiPatch<T> & thisMP = static_cast<const gsMultiPatch<T>&>(*this->m_fs );
655  GISMO_ASSERT( dynamic_cast<const gsMultiPatch<T>*>( other.m_fs ), "error");
656  const gsMultiPatch<T> & otherMP = static_cast<const gsMultiPatch<T>&>(*other.m_fs );
657  GISMO_ASSERT( (thisMP.domainDim()==otherMP.domainDim())&&
658  (thisMP.geoDim()==otherMP.geoDim())&&
659  (thisMP.coefsSize() == otherMP.coefsSize())&&
660  (thisMP.nPatches()==otherMP.nPatches()),
661  "The geometryMaps are not compatible!");
662 
663  // For every patch of the MultiPatch
664  for ( index_t p=0; p < thisMP.nPatches(); p++ )
665  {
666  // Copy coeffs of the other MultiPatch
667  thisMP.patch(p).coefs() = otherMP.patch(p).coefs();
668  }
669 
670  } //end copyCoffs
671 
672  void deformBy( const gsFeSolution<T> & displacement) const
673  {
674  const index_t dim = m_fs->domainDim();
675 
676  const gsMultiBasis<T> & mb = static_cast<const gsMultiBasis<T>&>(displacement.space().source());
677  const gsMultiPatch<T> & mp = static_cast<const gsMultiPatch<T>&>(*this->m_fs );
678  GISMO_ASSERT( dynamic_cast<const gsMultiBasis<T>*>(&displacement.space().source()), "error");
679  GISMO_ASSERT( dynamic_cast<const gsMultiPatch<T>*>( this->m_fs), "error");
680 
681  // For every patch of the MultiPatch
682  for ( size_t p=0; p < mp.nPatches(); p++ )
683  {
684  // Get the patch's coefficients
685  gsMatrix<T> &result = mp.patch(p).coefs();
686 
687  // Number of basis functions of patch with index p
688  const index_t sz = mb[p].size();
689 
690  // For all components
691  for (index_t c = 0; c!=dim; c++)
692  {
693  // loop over all basis functions (even the eliminated ones)
694  for (index_t i = 0; i < sz; ++i)
695  {
696  const int ii = displacement.mapper().index(i, p, c);
697  if ( displacement.mapper().is_free_index(ii) ) // DoF value is in the defVector
698  {
699  result(i,c) += displacement.coefs().at(ii);
700  }
701  else
702  {
703  result(i,c) += displacement.fixedPart().at( displacement.mapper().global_to_bindex(ii));
704  }
705  }
706  }
707  }
708  } // end deformBy
709 
710 public:
711  typedef T Scalar;
712 
713  friend class gismo::gsExprHelper<Scalar>;
714 
715  void print(std::ostream &os) const { os << "G"; }
716 
717  auto eval(const index_t k) const -> decltype(m_fd->values[0].col(k))
718  { return m_fd->values[0].col(k); }
719 
720 protected:
721 
722  gsGeometryMap() : m_fs(NULL), m_fd(NULL), m_isAcross(false) { }
723 
724  void setSource(const gsFunctionSet<Scalar> & fs) { m_fs = &fs;}
725  void setData(const gsMapData<Scalar> & val) { m_fd = &val;}
726 
727 public:
728 
729  index_t rows() const { return m_fs->targetDim(); }
730  index_t cols() const { return 1; }
731 
732  const gsFeSpace<T> & rowVar() const { return gsNullExpr<T>::get(); }
733  const gsFeSpace<T> & colVar() const { return gsNullExpr<T>::get(); }
734 
735  void parse(gsExprHelper<Scalar> & evList) const
736  {
737  evList.add(*this);
738  m_fd->flags |= NEED_VALUE;
739  }
740 };
741 
742 // Traits for gsGeometryMap
743 template <typename T> struct expr_traits<gsGeometryMap<T> >
744 {
745  typedef T Scalar;
746  typedef const gsGeometryMap<T> Nested_t; // nesting without ref!
747 };
748 
749 /*
750  Expression for the measure of a geometry map
751 */
752 template<class T>
753 class meas_expr : public _expr<meas_expr<T> >
754 {
755  typename gsGeometryMap<T>::Nested_t _G;
756 
757 public:
758  enum {Space = 0, ScalarValued = 1, ColBlocks = 0};
759 
760  typedef T Scalar;
761 
762  meas_expr(const gsGeometryMap<T> & G) : _G(G) { }
763 
764  T eval(const index_t k) const
765  {
766  return _G.data().measures.at(k);
767  }
768 
769  index_t rows() const { return 0; }
770  index_t cols() const { return 0; }
771  void parse(gsExprHelper<Scalar> & evList) const
772  {
773  evList.add(_G);
774  _G.data().flags |= NEED_MEASURE;
775  }
776 
777  const gsFeSpace<T> & rowVar() const { return gsNullExpr<T>::get(); }
778  const gsFeSpace<T> & colVar() const { return gsNullExpr<T>::get(); }
779 
780  void print(std::ostream &os) const { os << "meas("; _G.print(os); os <<")"; }
781 };
782 
783 
784 /*
785  An element object collecting relevant expressions
786 */
787 template<class T>
788 class gsFeElement
789 {
790  friend class cdiam_expr<T>;
791 
792  const gsDomainIterator<T> * m_di;
793 
794  const gsVector<T> * m_weights;
795  //const gsMatrix<T> * m_points;
796 
797  gsFeElement(const gsFeElement &);
798 public:
799  typedef T Scalar;
800 
801  gsFeElement() : m_di(NULL), m_weights(nullptr) { }
802 
803  void set(const gsDomainIterator<T> & di, const gsVector<T> & weights)
804  { m_di = &di, m_weights = &weights; }
805 
806  bool isValid() const { return nullptr!=m_weights; }
807 
808  const gsVector<T> & weights() const {return *m_weights;}
809 
810  template<class E>
811  integral_expr<E> integral(const _expr<E>& ff) const
812  { return integral_expr<E>(*this,ff); }
813 
814  typedef integral_expr<T> AreaRetType;
815  AreaRetType area() const
816  { return integral(_expr<T,true>(1)); }
817 
818  typedef integral_expr<meas_expr<T> > PHAreaRetType;
820  PHAreaRetType area(const gsGeometryMap<Scalar> & _G) const
821  { return integral(meas_expr<T>(_G)); }
822 
823  typedef pow_expr<integral_expr<T> > DiamRetType;
825  DiamRetType diam() const //-> int(1)^(1/d)
826  { return pow(integral(_expr<T,true>(1)),(T)(1)/(T)(2)); }
827 
828  typedef pow_expr<integral_expr<meas_expr<T> > > PHDiamRetType;
830  PHDiamRetType diam(const gsGeometryMap<Scalar> & _G) const
831  { return pow(integral(meas_expr<T>(_G)),(T)(1)/(T)(2)); }
832 
833  //const gsMatrix<T> points() const {return pts;}
834 
835  //index_t dim() { return di->
836 
837  void print(std::ostream &os) const { os << "e"; }
838 
839  void parse(gsExprHelper<T> & evList) const
840  {
841  GISMO_ERROR("EL");
842  }
843 };
844 
848 template<class E>
849 class integral_expr : public _expr<integral_expr<E> >
850 {
851 public:
852  //typedef typename E::Scalar Scalar;
853  typedef real_t Scalar;
854  mutable Scalar m_val;
855 private:
856  const gsFeElement<Scalar> & _e;
857  typename _expr<E>::Nested_t _ff;
858 public:
859  enum {Space= 0, ScalarValued= 1, ColBlocks = 0};
860 
861  integral_expr(const gsFeElement<Scalar> & el, const _expr<E> & u)
862  : m_val(-1), _e(el), _ff(u) { }
863 
864  const Scalar & eval(const index_t k) const
865  {
866  GISMO_ENSURE(_e.isValid(), "Element is valid within integrals only.");
867  // if (0==k)
868  {
869  const Scalar * w = _e.weights().data();
870  m_val = (*w) * _ff.val().eval(0);
871  for (index_t j = 1; j != _e.weights().rows(); ++j)
872  m_val += (*(++w)) * _ff.val().eval(j);
873  }
874  return m_val;
875  }
876 
877  inline const integral_expr<E> & val() const { return *this; }
878  inline index_t rows() const { return 0; }
879  inline index_t cols() const { return 0; }
880  void parse(gsExprHelper<Scalar> & evList) const
881  {
882  _ff.parse(evList);
883  }
884 
885  const gsFeSpace<Scalar> & rowVar() const { return gsNullExpr<Scalar>::get(); }
886  const gsFeSpace<Scalar> & colVar() const { return gsNullExpr<Scalar>::get(); }
887 
888  void print(std::ostream &os) const
889  {
890  os << "integral(";
891  _ff.print(os);
892  os <<")";
893  }
894 };
895 
896 /*
897  template<class T>
898  class parNv_expr : public _expr<parNv_expr<T> >
899  {
900  const gsFeElement<T> & _e;
901  public:
902  typedef T Scalar;
903 
904  enum {ScalarValued = 1};
905 
906  explicit parNv_expr(const gsFeElement<T> & el) : _e(el) { }
907 
908  T eval(const index_t k) const { return _e.m_di->getCellSize(); }
909 
910  inline cdiam_expr<T> val() const { return *this; }
911  inline index_t rows() const { return 0; }
912  inline index_t cols() const { return 0; }
913  void parse(gsExprHelper<Scalar> &) const { }
914  const gsFeVariable<T> & rowVar() const { gsNullExpr<Scalar>::get(); }
915  const gsFeVariable<T> & colVar() const { gsNullExpr<Scalar>::get(); }
916 
917  void print(std::ostream &os) const
918  { os << "diam(e)"; }
919  };
920 */
921 
922 template <typename T>
923 struct expr_traits<gsFeVariable<T> >
924 {
925  typedef T Scalar;
926  typedef const gsFeVariable<T> Nested_t;
927 };
928 
933 template<class T>
934 class gsFeVariable : public symbol_expr< gsFeVariable<T> >
935 {
936  friend class gismo::gsExprHelper<T>;
937  typedef symbol_expr< gsFeVariable<T> > Base;
938 protected:
939  explicit gsFeVariable(index_t _d = 1) : Base(_d) { }
940 public:
941  enum {Space = 0, ScalarValued = 0, ColBlocks = 0};
942 };
943 
944 template<class T>
945 class gsComposition : public symbol_expr< gsComposition<T> >
946 { //comp(f,G)
947  friend class gismo::gsExprHelper<T>;
948  typedef symbol_expr< gsComposition<T> > Base;
949  typename gsGeometryMap<T>::Nested_t _G;
950 protected:
951  explicit gsComposition(const gsGeometryMap<T> & G, index_t _d = 1)
952  : Base(_d), _G(G) { }
953 public:
954  enum {Space = 0, ScalarValued= 0, ColBlocks= 0};
955 
956  typename gsMatrix<T>::constColumn
957  eval(const index_t k) const { return this->m_fd->values[0].col(k); }
958 
959  const gsGeometryMap<T> & inner() const { return _G;};
960 
961  void parse(gsExprHelper<T> & evList) const
962  {
963  //evList.add(_G); //done in gsExprHelper
964  evList.add(*this);
965  this->data().flags |= NEED_VALUE|NEED_ACTIVE;
966  //_G.data().flags |= NEED_VALUE; //done in gsExprHelper
967  }
968 };
969 
970 
978 template<class T>
979 class gsFeSpace :public symbol_expr< gsFeSpace<T> >
980 {
981  friend class gsNullExpr<T>;
982 
983 protected:
984  typedef symbol_expr< gsFeSpace<T> > Base;
985 
986  // contains id, mapper, fixedDofs, etc
987  gsFeSpaceData<T> * m_sd;
988 
989 public:
990  enum{Space = 1, ScalarValued=0, ColBlocks=0};// test space
991 
992  typedef const gsFeSpace Nested_t; //no ref
993 
994  typedef T Scalar;
995 
996  const gsFeSpace<T> & rowVar() const {return *this;}
997 
998  gsDofMapper & mapper()
999  {
1000  GISMO_ASSERT(NULL!=m_sd, "Space/mapper not properly initialized.");
1001  return m_sd->mapper;
1002  }
1003 
1004  const gsDofMapper & mapper() const
1005  {return const_cast<gsFeSpace*>(this)->mapper();}
1006 
1007  inline const gsMatrix<T> & fixedPart() const {return m_sd->fixedDofs;}
1008  gsMatrix<T> & fixedPart() {return m_sd->fixedDofs;}
1009 
1010  index_t id() const { return (m_sd ? m_sd->id : -101); }
1011  void setSpaceData(gsFeSpaceData<T>& sd) {m_sd = &sd;}
1012 
1013  index_t interfaceCont() const {return m_sd->cont;}
1014  index_t & setInterfaceCont(const index_t _r) const
1015  {
1016  GISMO_ASSERT(_r>-2 && _r<1, "Invalid or not implemented (r="<<_r<<").");
1017  return m_sd->cont = _r;
1018  }
1019 
1020  gsFeSolution<T> function(const gsMatrix<T>& solVector) const
1021  { return gsFeSolution<T>(*this); }
1022 
1023  void getCoeffs(const gsMatrix<T>& solVector, gsMatrix<T> & result,
1024  const index_t p = 0) const
1025  {
1026  const index_t dim = this->dim();
1027 
1028  const gsMultiBasis<T> & mb = static_cast<const gsMultiBasis<T>&>(this->source());
1029  GISMO_ASSERT( dynamic_cast<const gsMultiBasis<T>*>(&this->source()), "error");
1030 
1031  // Reconstruct solution coefficients on patch p
1032  const index_t sz = mb[p].size();
1033  result.resize(sz, dim!=1 ? dim : solVector.cols()); // (!)
1034 
1035  for (index_t c = 0; c!=dim; c++) // for all components
1036  {
1037  // loop over all basis functions (even the eliminated ones)
1038  for (index_t i = 0; i < sz; ++i)
1039  {
1040  const int ii = m_sd->mapper.index(i, p, c);
1041  if ( m_sd->mapper.is_free_index(ii) ) // DoF value is in the solVector
1042  for(index_t s = 0; s != solVector.cols(); ++s )
1043  result(i,c+s) = solVector(ii,s); //assume dim==1 xor solVector.cols()==1
1044  else // eliminated DoF: fill with Dirichlet data
1045  result(i,c) = m_sd->fixedDofs.at( m_sd->mapper.global_to_bindex(ii) );
1046  }
1047  }
1048  }
1049 
1050  // space restrictTo(boundaries);
1051  // space restrictTo(bcRefList domain);
1052 
1053  void setupMapper(gsDofMapper dofsMapper) const
1054  {
1055  GISMO_ASSERT( dofsMapper.isFinalized(), "The provided dof-mapper is not finalized.");
1056  GISMO_ASSERT( dofsMapper.mapSize()==static_cast<size_t>(this->source().size()*dofsMapper.numComponents()), "The dof-mapper is not consistent: mapSize()="<<dofsMapper.mapSize()<<"!="<<static_cast<size_t>(this->source().size())<<"=this->source().size()");
1057  m_sd->mapper = give(dofsMapper);
1058  }
1059 
1060  void setup(const index_t _icont = -1) const
1061  {
1062  this->setInterfaceCont(_icont);
1063  m_sd->mapper = gsDofMapper();
1064 
1065  if (const gsMultiBasis<T> * mb =
1066  dynamic_cast<const gsMultiBasis<T>*>(&this->source()) )
1067  {
1068  m_sd->mapper = gsDofMapper(*mb, this->dim() );
1069  //m_mapper.init(*mb, this->dim()); //bug
1070  if ( 0==this->interfaceCont() ) // Conforming boundaries ?
1071  {
1072  for ( gsBoxTopology::const_iiterator it = mb->topology().iBegin();
1073  it != mb->topology().iEnd(); ++it )
1074  {
1075  mb->matchInterface(*it, m_sd->mapper);
1076  }
1077  }
1078  }
1079 
1080  if (const gsMappedBasis<2,T> * mb =
1081  dynamic_cast<const gsMappedBasis<2,T>*>(&this->source()) )
1082  {
1083  m_sd->mapper.setIdentity(mb->nPatches(), mb->size() , this->dim());
1084  }
1085 
1086  m_sd->mapper.finalize();
1087  }
1088 
1089  void setup(const gsBoundaryConditions<T> & bc, const index_t dir_values,
1090  const index_t _icont = -1) const
1091  {
1092  this->setInterfaceCont(_icont);
1093  m_sd->mapper = gsDofMapper();
1094  const index_t dim = this->dim();
1095  const gsMultiBasis<T> *mb = dynamic_cast<const gsMultiBasis<T> *>(&this->source());
1096  if (mb != nullptr)
1097  {
1098  m_sd->mapper = gsDofMapper(*mb, this->dim());
1099  //m_mapper.init(*mb, this->dim()); //bug
1100  if (0 == this->interfaceCont()) // Conforming boundaries ?
1101  {
1102  for (gsBoxTopology::const_iiterator it = mb->topology().iBegin();
1103  it != mb->topology().iEnd(); ++it) {
1104  mb->matchInterface(*it, m_sd->mapper);
1105  }
1106  }
1107 
1108  // Strong Dirichlet conditions
1109  gsMatrix<index_t> bnd;
1110  for (typename gsBoundaryConditions<T>::const_iterator
1111  it = bc.begin("Dirichlet"); it != bc.end("Dirichlet"); ++it)
1112  {
1113  const index_t cc = it->unkComponent();
1114  GISMO_ASSERT(static_cast<size_t>(it->ps.patch) < this->mapper().numPatches(),
1115  "Problem: a boundary condition is set on a patch id which does not exist.");
1116 
1117  bnd = mb->basis(it->ps.patch).boundary(it->ps.side());
1118  m_sd->mapper.markBoundary(it->ps.patch, bnd, cc);
1119  }
1120  // Clamped boundary condition (per DoF)
1121  gsMatrix<index_t> bnd1;
1122  for (typename gsBoundaryConditions<T>::const_iterator
1123  it = bc.begin("Clamped"); it != bc.end("Clamped"); ++it)
1124  {
1125  const index_t cc = it->unkComponent();
1126 
1127  GISMO_ASSERT(static_cast<size_t>(it->ps.patch) < this->mapper().numPatches(),
1128  "Problem: a boundary condition is set on a patch id which does not exist.");
1129 
1130  bnd = mb->basis(it->ps.patch).boundaryOffset(it->ps.side(), 0);
1131  bnd1 = mb->basis(it->ps.patch).boundaryOffset(it->ps.side(), 1);
1132  // Cast to tensor b-spline basis
1133  if (!it->ps.parameter())
1134  bnd.swap(bnd1);
1135  for (index_t c = 0; c!=dim; c++) // for all components
1136  {
1137  if (c==cc || cc==-1 )
1138  for (index_t k = 0; k < bnd.size(); ++k)
1139  m_sd->mapper.matchDof(it->ps.patch, (bnd)(k, 0),
1140  it->ps.patch, (bnd1)(k, 0), c);
1141  }
1142 
1143  }
1144 
1145  // Collapsed
1146  for (typename gsBoundaryConditions<T>::const_iterator
1147  it = bc.begin("Collapsed"); it != bc.end("Collapsed"); ++it)
1148  {
1149  const index_t cc = it->unkComponent();
1150 
1151  GISMO_ASSERT(static_cast<size_t>(it->ps.patch) < this->mapper().numPatches(),
1152  "Problem: a boundary condition is set on a patch id which does not exist.");
1153 
1154  bnd = mb->basis(it->ps.patch).boundary(it->ps.side());
1155 
1156  // match all DoFs to the first one of the side
1157  for (index_t c = 0; c!=dim; c++) // for all components
1158  {
1159  if (c==cc || cc==-1)
1160  for (index_t k = 0; k < bnd.size() - 1; ++k)
1161  m_sd->mapper.matchDof(it->ps.patch, (bnd)(0, 0),
1162  it->ps.patch, (bnd)(k + 1, 0), c);
1163  }
1164  }
1165 
1166  // Coupled
1167  for (typename gsBoundaryConditions<T>::const_cpliterator
1168  it = bc.coupledBegin(); it != bc.coupledEnd(); ++it)
1169  {
1170  const index_t cc = it->component;
1171 
1172  GISMO_ASSERT(static_cast<size_t>(it->ifc.first().patch) < this->mapper().numPatches(),
1173  "Problem: a boundary condition is set on a patch id which does not exist.");
1174  GISMO_ASSERT(static_cast<size_t>(it->ifc.second().patch) < this->mapper().numPatches(),
1175  "Problem: a boundary condition is set on a patch id which does not exist.");
1176 
1177 
1178  bnd = mb->basis(it->ifc.first().patch).boundary(it->ifc.first().side());
1179  bnd1 = mb->basis(it->ifc.second().patch).boundary(it->ifc.second().side());
1180 
1181  // match all DoFs to the first one of the side
1182  for (index_t c = 0; c!=dim; c++) // for all components
1183  {
1184  if (c==cc || cc==-1)
1185  {
1186  for (index_t k = 0; k < bnd.size() - 1; ++k)
1187  m_sd->mapper.matchDof(it->ifc.first() .patch, (bnd)(0, 0),
1188  it->ifc.first() .patch, (bnd)(k + 1, 0), c);
1189  for (index_t k = 0; k < bnd1.size(); ++k)
1190  m_sd->mapper.matchDof(it->ifc.first() .patch, (bnd)(0, 0),
1191  it->ifc.second().patch, (bnd1)(k, 0), c);
1192  }
1193  }
1194  }
1195 
1196  // corners
1197  for (typename gsBoundaryConditions<T>::const_citerator
1198  it = bc.cornerBegin(); it != bc.cornerEnd(); ++it)
1199  {
1200  for (index_t r = 0; r!=this->dim(); ++r)
1201  {
1202  if (it->component!=-1 && r!=it->component) continue;
1203 
1204  //assumes (unk == -1 || it->unknown == unk)
1205  GISMO_ASSERT(static_cast<size_t>(it->patch) < mb->nBases(),
1206  "Problem: a corner boundary condition is set on a patch id which does not exist.");
1207  m_sd->mapper.eliminateDof(mb->basis(it->patch).functionAtCorner(it->corner),
1208  it->patch, it->component);
1209  }
1210  }
1211 
1212  } else if (const gsBasis<T> *b =
1213  dynamic_cast<const gsBasis<T> *>(&this->source()))
1214  {
1215  m_sd->mapper = gsDofMapper(*b, this->dim() );
1216  gsMatrix<index_t> bnd;
1217  for (typename gsBoundaryConditions<T>::const_iterator
1218  it = bc.begin("Dirichlet"); it != bc.end("Dirichlet"); ++it) {
1219  GISMO_ASSERT(it->ps.patch == 0,
1220  "Problem: a boundary condition is set on a patch id which does not exist.");
1221 
1222  bnd = b->boundary(it->ps.side());
1223  m_sd->mapper.markBoundary(0, bnd, it->unkComponent());
1224  }
1225 
1226  for (typename gsBoundaryConditions<T>::const_iterator
1227  it = bc.begin("Clamped"); it != bc.end("Clamped"); ++it) {
1228  GISMO_ASSERT(it->ps.patch == 0,
1229  "Problem: a boundary condition is set on a patch id which does not exist.");
1230 
1231  bnd = b->boundary(it->ps.side());
1232  //const index_t cc = it->unkComponent();
1233  // m_sd->mapper.markBoundary(0, bnd, 0);
1234  }
1235 
1236  m_sd->mapper = gsDofMapper(*b);
1237  for (typename gsBoundaryConditions<T>::const_iterator
1238  it = bc.begin("Collapsed"); it != bc.end("Collapsed"); ++it) {
1239  GISMO_ASSERT(it->ps.patch == 0,
1240  "Problem: a boundary condition is set on a patch id which does not exist.");
1241 
1242  bnd = b->boundary(it->ps.side());
1243  //const index_t cc = it->unkComponent();
1244  // m_sd->mapper.markBoundary(0, bnd, 0);
1245  }
1246  } else if (const gsMappedBasis<2, T> *mapb =
1247  dynamic_cast<const gsMappedBasis<2, T> *>(&this->source()))
1248  {
1249  m_sd->mapper.setIdentity(mapb->nPatches(), mapb->size(), this->dim());
1250 
1251  if (0 == this->interfaceCont()) // C^0 matching interface
1252  {
1253  gsMatrix<index_t> int1, int2;
1254  for (gsBoxTopology::const_iiterator it = mapb->getTopol().iBegin();
1255  it != mapb->getTopol().iEnd(); ++it) {
1256  int1 = mapb->basis(it->first().patch).boundaryOffset(it->first().side(), 0);
1257  int2 = mapb->basis(it->second().patch).boundaryOffset(it->second().side(), 0);
1258 
1259  m_sd->mapper.matchDofs(it->first().patch, int1, it->second().patch, int2);
1260  }
1261  }
1262  if (1 == this->interfaceCont()) // C^1 matching interface
1263  {
1264  GISMO_ERROR("Boundary offset function is not implemented for gsMappedBasis in general.");
1265  }
1266 
1267  gsMatrix<index_t> bnd;
1268  for (typename gsBoundaryConditions<T>::const_iterator
1269  it = bc.begin("Dirichlet"); it != bc.end("Dirichlet"); ++it) {
1270  const index_t cc = it->unkComponent();
1271  GISMO_ASSERT(static_cast<size_t>(it->ps.patch) < this->mapper().numPatches(),
1272  "Problem: a boundary condition is set on a patch id which does not exist.");
1273 
1274  bnd = mapb->basis(it->ps.patch).boundary(it->ps.side());
1275  m_sd->mapper.markBoundary(it->ps.patch, bnd, cc);
1276  }
1277 
1278  // Clamped boundary condition (per DoF)
1279  gsMatrix<index_t> bnd1;
1280  for (typename gsBoundaryConditions<T>::const_iterator
1281  it = bc.begin("Clamped"); it != bc.end("Clamped"); ++it) {
1282  const index_t cc = it->unkComponent();
1283 
1284  GISMO_ASSERT(static_cast<size_t>(it->ps.patch) < this->mapper().numPatches(),
1285  "Problem: a boundary condition is set on a patch id which does not exist.");
1286 
1287  bnd = mapb->basis(it->ps.patch).boundaryOffset(it->ps.side(), 0);
1288  bnd1 = mapb->basis(it->ps.patch).boundaryOffset(it->ps.side(), 1);
1289 
1290  // Cast to tensor b-spline basis
1291  if (mapb != NULL) // clamp adjacent dofs
1292  {
1293  if (!it->ps.parameter())
1294  bnd.swap(bnd1);
1295  for (index_t c = 0; c!=dim; c++) // for all components
1296  {
1297  if (c==cc || cc==-1 )
1298  for (index_t k = 0; k < bnd.size() - 1; ++k)
1299  m_sd->mapper.matchDof( it->ps.patch, (bnd)(k, 0),
1300  it->ps.patch, (bnd1)(k, 0), c);
1301  }
1302  } else
1303  gsWarn << "Unable to apply clamped condition.\n";
1304  }
1305 
1306  // COLLAPSED
1307  for (typename gsBoundaryConditions<T>::const_iterator
1308  it = bc.begin("Collapsed"); it != bc.end("Collapsed"); ++it) {
1309  const index_t cc = it->unkComponent();
1310 
1311  GISMO_ASSERT(static_cast<size_t>(it->ps.patch) < this->mapper().numPatches(),
1312  "Problem: a boundary condition is set on a patch id which does not exist.");
1313 
1314  bnd = mapb->basis(it->ps.patch).boundary(it->ps.side());
1315 
1316  // Cast to tensor b-spline basis
1317  if (mapb != NULL) // clamp adjacent dofs
1318  {
1319  // match all DoFs to the first one of the side
1320  for (index_t c = 0; c!=dim; c++) // for all components
1321  {
1322  if (c==cc || cc==-1)
1323  for (index_t k = 0; k < bnd.size() - 1; ++k)
1324  m_sd->mapper.matchDof(it->ps.patch, (bnd)(0, 0),
1325  it->ps.patch, (bnd)(k + 1, 0), c);
1326  }
1327  }
1328  }
1329 
1330  // corners
1331  for (typename gsBoundaryConditions<T>::const_citerator
1332  it = bc.cornerBegin(); it != bc.cornerEnd(); ++it)
1333  {
1334  //assumes (unk == -1 || it->unknown == unk)
1335  GISMO_ASSERT(it->patch < mapb->nPieces(),
1336  "Problem: a corner boundary condition is set on a patch id which does not exist.");
1337  m_sd->mapper.eliminateDof(mapb->basis(it->patch).functionAtCorner(it->corner), it->patch, it->component);
1338  }
1339  } else
1340  {
1341  GISMO_ASSERT(0 == bc.size(), "Problem: BCs are ignored.");
1342  m_sd->mapper.setIdentity(this->source().nPieces(), this->source().size());
1343  }
1344 
1345  m_sd->mapper.finalize();
1346 
1347  // Compute Dirichlet node values
1348  gsDirichletValues(bc, dir_values, *this);
1349  }
1350 
1351  void print(std::ostream &os) const { os << "u"; }
1352 
1353 protected:
1354  friend class gismo::gsExprHelper<Scalar>;
1355  friend class symbol_expr<gsFeSpace>;
1356  explicit gsFeSpace(index_t _d = 1) : Base(_d), m_sd(nullptr) { }
1357 };
1358 
1359 template<class T> inline bool
1360 operator== (const gsFeSpace<T> & a, const gsFeSpace<T> & b)
1361 { return a.id()== b.id() && a.isAcross()==b.isAcross(); }
1362 
1363 /*
1364  Expression representing a function given by a vector of
1365  coefficients in a gsFeSpace.
1366 
1367  Typically it used for accessing the solution of a boundary-value
1368  problem.
1369 */
1370 template<class T>
1371 class gsFeSolution : public _expr<gsFeSolution<T> >
1372 {
1373 protected:
1374  const gsFeSpace<T> _u;
1375  gsMatrix<T> * _Sv;
1376  bool m_isAcross;
1377 
1378 public:
1379  typedef T Scalar;
1380  enum {Space = 0, ScalarValued= 0, ColBlocks= 0};
1381 
1382  bool isAcross() const { return m_isAcross; }
1383 
1384  gsFeSolution right() const
1385  {
1386  gsFeSolution ac(*this);
1387  ac.m_isAcross = true;
1388  return ac;
1389  }
1390 
1391  gsFeSolution left() const { return gsFeSolution(*this); }
1392 
1393  explicit gsFeSolution(const gsFeSpace<T> & u) : _u(u), _Sv(NULL) { }
1394 
1395  gsFeSolution(const gsFeSpace<T> & u, gsMatrix<T> & Sv) : _u(u), _Sv(&Sv) { }
1396 
1397  const gsFeSpace<T> & space() const {return _u;};
1398 
1399  mutable gsMatrix<T> res;
1400  const gsMatrix<T> & eval(index_t k) const
1401  {
1402  bool singleActives = (1 == _u.data().actives.cols());
1403 
1404  res.setZero(_u.dim(), 1);
1405  const gsDofMapper & map = _u.mapper();
1406  GISMO_ASSERT(_Sv->size()==map.freeSize(), "The solution vector has wrong dimensions: "<<_Sv->size()<<" != "<<map.freeSize());
1407 
1408  for (index_t c = 0; c!=_u.dim(); c++) // for all components
1409  {
1410  for (index_t i = 0; i!=_u.data().actives.rows(); ++i)
1411  {
1412  const index_t ii = map.index(_u.data().actives(i, singleActives ? 0 : k), _u.data().patchId, c);
1413  if ( map.is_free_index(ii) ) // DoF value is in the solVector
1414  res.at(c) += _Sv->at(ii) * _u.data().values[0](i,k);
1415  else
1416  res.at(c) += _u.data().values[0](i,k) *
1417  _u.fixedPart().at( map.global_to_bindex(ii) );
1418  }
1419  }
1420  return res;
1421  }
1422 
1423  //template<class U>
1424  //linearComb(U & ie){ sum up ie[_u] times the _Sv }
1425  // ie.eval(k), _u.data().actives(), fixedPart() - see lapl_expr
1426 
1427  const gsFeSpace<Scalar> & rowVar() const {return gsNullExpr<Scalar>::get();}
1428  const gsFeSpace<Scalar> & colVar() const {return gsNullExpr<Scalar>::get();}
1429  index_t rows() const {return _u.dim(); }
1430 
1431  static index_t cols() {return 1; }
1432 
1433  void parse(gsExprHelper<Scalar> & evList) const
1434  {
1435  evList.add(_u);
1436  _u.data().flags |= NEED_VALUE | NEED_ACTIVE;
1437  }
1438 
1439  void print(std::ostream &os) const { os << "s"; }
1440 
1441 public:
1442  index_t dim() const { return _u.dim();}
1443 
1444  index_t parDim() const
1445  { return _u.source().domainDim(); }
1446 
1447  //gsDofMapper & mapper() {return _u.mapper();}
1448  const gsDofMapper & mapper() const {return _u.mapper();}
1449 
1450  inline const gsMatrix<T> & fixedPart() const {return _u.fixedPart();}
1451  gsMatrix<T> & fixedPart() {return _u.fixedPart();}
1452 
1453  //gsFuncData<T> & data() {return _u.data();}
1454  const gsFuncData<T> & data() const {return _u.data();}
1455 
1456  void setSolutionVector(gsMatrix<T>& solVector)
1457  { _Sv = & solVector; }
1458 
1464  void setComponent(index_t component, real_t value, index_t patch=-1)
1465  {
1466  gsMatrix<T> & solVector = *_Sv;
1467  const gsDofMapper & mapper = _u.mapper();
1468 
1469  index_t patchStart, patchEnd;
1470  if (patch==-1){
1471  patchStart = 0;
1472  patchEnd = _u.mapper().numPatches();
1473  }
1474  else{
1475  patchStart = patch;
1476  patchEnd = patch + 1;
1477  }
1478 
1479  for (index_t p=patchStart; p!=patchEnd; ++p)
1480  {
1481  for (size_t i = 0; i != mapper.patchSize(p, component); ++i)
1482  {
1483  const index_t ii = mapper.index(i, p, component);
1484  if ( mapper.is_free_index(ii) ) // DoF value is in the solVector
1485  solVector.at(ii) = value;
1486  }
1487  }
1488  }
1489 
1490  const gsMatrix<T> & coefs() const { return *_Sv; }
1491  //gsMatrix<T> & coefs() { return *_Sv; } // wd4702 ?
1492 
1493  //const gsMatrix<T> & coefs(component, patch) const { return *_Sv; }
1494 
1496  void perturbLocal(T val, index_t j, index_t p = 0)
1497  {
1498  // GISMO_ASSERT(1==_u.data().actives.cols(), "Single actives expected");
1499  //if (_u.mapper().is_free_index(j) )
1500  //{
1501  GISMO_ASSERT(j<_Sv->size(), "Solution vector is not initialized/allocated, sz="<<_Sv->size() );
1502  _Sv->at(j) += val;
1503  //}
1504  //else
1505  // _u.fixedPart().at( _u.mapper().global_to_bindex(j) ) += val;
1506  }
1507 
1509  void extract(gsMatrix<T> & result, const index_t p = 0) const
1510  { _u.getCoeffs(*_Sv, result, p); }
1511 
1514  void extractFull(gsMatrix<T> & result) const
1515  {
1516  index_t offset;
1517  const index_t dim = _u.dim();
1518  const size_t totalSz = _u.mapper().mapSize();
1519  result.resize(totalSz, 1);
1520  for (size_t p=0; p!=_u.mapper().numPatches(); ++p)
1521  {
1522  offset = _u.mapper().offset(p);
1523  // Reconstruct solution coefficients on patch p
1524 
1525  for (index_t c = 0; c!=dim; c++) // for all components
1526  {
1527  const index_t sz = _u.mapper().patchSize(p,c);
1528 
1529  // loop over all basis functions (even the eliminated ones)
1530  for (index_t i = 0; i < sz; ++i)
1531  {
1532  //gsDebugVar(i);
1533  const int ii = _u.mapper().index(i, p, c);
1534  //gsDebugVar(ii);
1535  if ( _u.mapper().is_free_index(ii) ) // DoF value is in the solVector
1536  {
1537  result(i+offset,0) = _Sv->at(ii);
1538  }
1539  else // eliminated DoF: fill with Dirichlet data
1540  result(i+offset,0) = _u.fixedPart().at( _u.mapper().global_to_bindex(ii) );
1541  }
1542  offset += sz;
1543  }
1544  }
1545  }
1546 
1548  void extract(gsMultiPatch<T> & result) const
1549  {
1550  result.clear();
1551 
1552  if( const gsMultiBasis<T>* basis = dynamic_cast<const gsMultiBasis<T>* >(&_u.source()) )
1553  for (size_t i = 0; i != basis->nBases(); ++i)
1554  {
1555  memory::unique_ptr<gsGeometry<T> > p(this->extractPiece(i));
1556  result.addPatch(*p);
1557  }
1558  }
1559 
1561  void extract(gsMappedSpline<2,T> & result) const
1562  {
1563  if( const gsMappedBasis<2,T>* basis = dynamic_cast<const gsMappedBasis<2,T>* >(&_u.source()) )
1564  {
1565  gsMatrix<T> coefs;
1566  this->extractFull(coefs);
1567  coefs.resize(coefs.rows()/_u.dim(),_u.dim());
1568  result.init(*basis,coefs);
1569  }
1570  }
1571 
1573  memory::unique_ptr<gsGeometry<T> > extractPiece(const index_t p) const
1574  {
1575  if ( const gsBasis<T> * b = dynamic_cast<const gsBasis<T>*>(&_u.source().piece(p)) )
1576  {
1577  gsMatrix<T> cf;
1578  extract(cf, p);
1579  return b->makeGeometry(give(cf));
1580  }
1581  GISMO_ERROR("gsFeSolution: Extraction error");
1582  }
1583 
1584  // insert g-coefficients to the solution vector
1585  void insert(const gsGeometry<T> & g, const index_t p = 0) const
1586  {
1587  const gsMatrix<T> & cf = g.coefs();
1588  gsMatrix<T> & sol = *_Sv;
1589  //gsMatrix<T> & fixedPart = _u.fixedPart();
1590  const gsDofMapper & mapper = _u.mapper();
1591  for (index_t c = 0; c!=_u.dim(); c++) // for all components
1592  {
1593  for (index_t i = 0; i != cf.rows(); ++i)
1594  {
1595  const index_t ii = mapper.index(i, p, c);
1596  if ( mapper.is_free_index(ii) ) // DoF value is in the solVector
1597  sol.at(ii) = cf(i, c);
1598  /*
1599  else
1600  {
1601  fixedPart.row(m_sd->mapper.global_to_bindex(ii)) = cf.row(i);
1602  }
1603  */
1604  }
1605  }
1606  }
1607 };
1608 
1609 /*
1610  Expression for the transpose of an expression
1611 */
1612 template<class E, bool cw>
1613 class tr_expr : public _expr<tr_expr<E,cw> >
1614 {
1615  typename E::Nested_t _u;
1616 
1617 public:
1618 
1619  typedef typename E::Scalar Scalar;
1620 
1621  tr_expr(_expr<E> const& u)
1622  : _u(u) { }
1623 
1624 public:
1625  enum {ColBlocks = E::ColBlocks, ScalarValued=E::ScalarValued};
1626  enum {Space = cw?E::Space:(E::Space==1?2:(E::Space==2?1:E::Space))};
1627 
1628  mutable Temporary_t res;
1629  const Temporary_t & eval(const index_t k) const
1630  {
1631  if (E::ColBlocks)
1632  res = _u.eval(k).blockTranspose( _u.cardinality() );
1633  else
1634  res = _u.eval(k).transpose();
1635  return res;
1636  }
1637 
1638  index_t rows() const { return _u.cols(); }
1639 
1640  index_t cols() const { return _u.rows(); }
1641 
1642  void parse(gsExprHelper<Scalar> & evList) const
1643  { _u.parse(evList); }
1644 
1645  const gsFeSpace<Scalar> & rowVar() const { return cw?_u.rowVar():_u.colVar(); }
1646  const gsFeSpace<Scalar> & colVar() const { return cw?_u.colVar():_u.rowVar(); }
1647 
1648  index_t cardinality_impl() const { return _u.cardinality_impl(); }
1649 
1650  void print(std::ostream &os) const { os<<"("; _u.print(os); os <<")\u1D40"; }
1651 private:
1652 /*
1653  template<class U> EIGEN_STRONG_INLINE MatExprType
1654  eval_impl(const U k, typename util::enable_if<1==ColBlocks,U>::type* = nullptr)
1655  { return _u.eval(k).blockTranspose(_u.cols()/_u.rows()); }
1656 
1657  template<class U> EIGEN_STRONG_INLINE MatExprType
1658  eval_impl(const U k, typename util::enable_if<0==ColBlocks,U>::type* = nullptr)
1659  { return _u.eval(k).transpose(); }
1660 */
1661 };
1662 
1663 /*
1664  Expression to make an expression colblocks
1665 */
1666 template<class E>
1667 class cb_expr : public _expr<cb_expr<E> >
1668 {
1669  typename E::Nested_t _u;
1670 
1671 public:
1672 
1673  typedef typename E::Scalar Scalar;
1674 
1675  cb_expr(_expr<E> const& u)
1676  : _u(u) { }
1677 
1678 public:
1679  enum {ColBlocks = 1, ScalarValued=E::ScalarValued};
1680  enum {Space = E::Space};
1681 
1682  mutable gsMatrix<Scalar> ev, res;
1683 
1684  const gsMatrix<Scalar> & eval(const index_t k) const
1685  {
1686  //return _u.eval(k).transpose();
1687  // /*
1688  ev = _u.eval(k);
1689  if (E::ColBlocks)
1690  {
1691  return ev;
1692  }
1693  else
1694  {
1695  res = _u.eval(k);
1696 
1697  GISMO_ASSERT(res.rows() % _u.rows() == 0 && res.cols() % _u.cols() == 0,"Result is not a multiple of the space dimensions?");
1698 
1699  index_t cardinality;
1700  if ( (cardinality = res.rows() / _u.rows()) >= 1 && res.cols() / _u.cols() == 1 ) // stored in rows
1701  {
1702  res.resize(_u.rows(), cardinality * _u.cols());
1703  for (index_t r = 0; r!=cardinality; r++)
1704  res.block(0 , r * _u.cols(), _u.rows(), _u.cols()) = ev.block( r * _u.rows(), 0, _u.rows(), _u.cols() );
1705  }
1706  else if ( (cardinality = res.rows() / _u.rows()) == 1 && res.cols() / _u.cols() >= 1 ) // stored in cols ----->>>> This is already colBlocks???
1707  {
1708  res.resize(_u.rows(), cardinality * _u.cols());
1709  for (index_t r = 0; r!=cardinality; r++)
1710  res.block(0 , r * _u.cols(), _u.rows(), _u.cols()) = ev.block( 0, r * _u.cols(), _u.rows(), _u.cols() );
1711  }
1712  }
1713  return res;
1714  }
1715 
1716  index_t rows() const { return _u.rows(); }
1717 
1718  index_t cols() const { return _u.cols(); }
1719 
1720  void parse(gsExprHelper<Scalar> & evList) const
1721  { _u.parse(evList); }
1722 
1723  const gsFeSpace<Scalar> & rowVar() const { return _u.rowVar(); }
1724  const gsFeSpace<Scalar> & colVar() const { return _u.colVar(); }
1725 
1726  index_t cardinality_impl() const
1727  {
1728  res = _u.eval(0);
1729  index_t cardinality;
1730  if ( res.rows() / _u.rows() >= 1 && res.cols() / _u.cols() == 1 ) // stored in rows
1731  cardinality = res.rows() / _u.rows();
1732  else if ( res.rows() / _u.rows() == 1 && res.cols() / _u.cols() >= 1 )
1733  cardinality = res.cols() / _u.cols();
1734  else
1735  GISMO_ERROR("Cardinality for cb_expr cannot be determined.");
1736 
1737  return cardinality;
1738  }
1739 
1740  void print(std::ostream &os) const { os<<"{"; _u.print(os); os <<"}"; }
1741 };
1742 
1743 
1744 /*
1745  Expression for an evaluation of the (sub-)expression in temporary memory
1746 */
1747 template<class E>
1748 class temp_expr : public _expr<temp_expr<E> >
1749 {
1750  typename E::Nested_t _u;
1751  typedef typename E::Scalar Scalar;
1752  mutable gsMatrix<Scalar> tmp;
1753 
1754 public:
1755  temp_expr(_expr<E> const& u)
1756  : _u(u) { }
1757 
1758 public:
1759  enum {Space = E::Space, ScalarValued = E::ScalarValued,
1760  ColBlocks = E::ColBlocks};
1761 
1762  // template<bool S = ColBlocks>
1763  // typename util::enable_if<S,MatExprType>::type
1764  const gsMatrix<Scalar> & eval(const index_t k) const
1765  {
1766  tmp = _u.eval(k);
1767  return tmp;
1768  }
1769 
1770  index_t rows() const { return _u.rows(); }
1771  index_t cols() const { return _u.cols(); }
1772  void parse(gsExprHelper<Scalar> & evList) const { _u.parse(evList); }
1773  const gsFeSpace<Scalar> & rowVar() const { return _u.rowVar(); }
1774  const gsFeSpace<Scalar> & colVar() const { return _u.colVar(); }
1775 
1776  void print(std::ostream &os) const { _u.print(os); }
1777 };
1778 
1779 /*
1780  Expression for the trace of a (matrix) expression
1781 */
1782 template<class E>
1783 class trace_expr : public _expr<trace_expr<E> >
1784 {
1785 public:
1786  typedef typename E::Scalar Scalar;
1787  enum {ScalarValued = 0, Space = E::Space, ColBlocks= 0};
1788 
1789 private:
1790  typename E::Nested_t _u;
1791  mutable gsMatrix<Scalar> res;
1792 
1793 public:
1794  trace_expr(_expr<E> const& u) : _u(u)
1795  {
1796  // gcc 4.8.4: invalid read due to _u.rows() using gsFuncData
1797  //GISMO_ASSERT(0== _u.cols()%_u.rows(), "Expecting square-block expression, got " << _u.rows() <<" x "<< _u.cols() );
1798  }
1799 
1800  // choose if ColBlocks
1801  const gsMatrix<Scalar> & eval(const index_t k) const
1802  {
1803  auto tmp = _u.eval(k);
1804  const index_t cb = _u.rows();
1805  const index_t r = _u.cardinality();
1806  if (Space==1)
1807  res.resize(r, 1);
1808  else
1809  res.resize(1, r);
1810 
1811  for (index_t i = 0; i!=r; ++i)
1812  res.at(i) = tmp.middleCols(i*cb,cb).trace();
1813  return res;
1814  }
1815 
1816  // choose if !ColBlocks
1817  //todo: Scalar eval(const index_t k) const
1818 
1819  index_t rows() const { return _u.cols() / _u.rows(); } //_u.cardinality()?
1820  index_t cols() const { return 1; }
1821 
1822  index_t cardinality_impl() const { return _u.cardinality(); }
1823 
1824  void parse(gsExprHelper<Scalar> & evList) const
1825  { _u.parse(evList); }
1826 
1827  const gsFeSpace<Scalar> & rowVar() const { return _u.rowVar(); }
1828  const gsFeSpace<Scalar> & colVar() const { return _u.colVar(); }
1829 
1830  void print(std::ostream &os) const { os << "trace("; _u.print(os); os<<")"; }
1831 };
1832 
1833 /*
1834  Expression for the adjugate of a (matrix) expression
1835 */
1836 template<class E>
1837 class adjugate_expr : public _expr<adjugate_expr<E> >
1838 {
1839 public:
1840  typedef typename E::Scalar Scalar;
1841  enum {ScalarValued = 0, ColBlocks = E::ColBlocks};
1842  enum {Space = E::Space};
1843 private:
1844  typename E::Nested_t _u;
1845  mutable gsMatrix<Scalar> res;
1846 
1847 public:
1848  adjugate_expr(_expr<E> const& u) : _u(u)
1849  {
1850  // gcc 4.8.4: invalid read due to _u.rows() using gsFuncData
1851  //GISMO_ASSERT(0== _u.cols()%_u.rows(), "Expecting square-block expression, got " << _u.rows() <<" x "<< _u.cols() );
1852  }
1853 
1854  // choose if ColBlocks
1855  const gsMatrix<Scalar> & eval(const index_t k) const
1856  {
1857  auto tmp = _u.eval(k);
1858  const index_t cb = _u.rows();
1859  const index_t r = _u.cols() / cb;
1860  res.resize(_u.rows(),_u.cols());
1861  for (index_t i = 0; i!=r; ++i){
1862  res.middleCols(i*cb,cb) = tmp.middleCols(i*cb,cb).adjugate();
1863  }
1864  return res;
1865  }
1866 
1867  // choose if !ColBlocks
1868  //todo: Scalar eval(const index_t k) const
1869 
1870  index_t rows() const { return _u.rows(); }
1871  index_t cols() const { return _u.cols(); }
1872 
1873  void parse(gsExprHelper<Scalar> & evList) const
1874  { _u.parse(evList); }
1875 
1876  const gsFeSpace<Scalar> & rowVar() const { return _u.rowVar(); }
1877  const gsFeSpace<Scalar> & colVar() const { return _u.colVar(); }
1878 
1879  void print(std::ostream &os) const { os << "adj("; _u.print(os); os<<")"; }
1880 };
1881 
1882 template<class E>
1883 class reshape_expr : public _expr<reshape_expr<E> >
1884 {
1885 public:
1886  typedef typename E::Scalar Scalar;
1887  enum {ScalarValued = 0, ColBlocks = E::ColBlocks};
1888  enum {Space = E::Space};
1889 private:
1890  typename E::Nested_t _u;
1891  index_t _n, _m;
1892  mutable gsMatrix<Scalar> tmp;
1893 
1894 public:
1895 
1896  //the reshaping is done column-wise
1897  reshape_expr(_expr<E> const& u, index_t n, index_t m) : _u(u), _n(n), _m(m)
1898  {
1899  //GISMO_ASSERT( _u.rows()*_u.cols() == _n*_m, "Wrong dimension"); //
1900  }
1901 
1902  const gsAsConstMatrix<Scalar> eval(const index_t k) const
1903  {
1904  // Note: this assertion would fail in the constructor!
1905  GISMO_ASSERT( _u.rows()*_u.cols() == _n*_m, "Wrong dimension "
1906  << _u.rows() << " x "<<_u.cols() << "!=" << _n << " * "<< _m );
1907  tmp = _u.eval(k);
1908  return gsAsConstMatrix<Scalar>(tmp.data(),_n,_m);
1909  }
1910 
1911  index_t rows() const { return _n; }
1912  index_t cols() const { return _m; }
1913 
1914  void parse(gsExprHelper<Scalar> & evList) const
1915  { _u.parse(evList); }
1916 
1917  const gsFeSpace<Scalar> & rowVar() const { return _u.rowVar(); }
1918  const gsFeSpace<Scalar> & colVar() const { return _u.colVar(); }
1919 
1920  void print(std::ostream &os) const { os << "reshape("; _u.print(os); os<<","<<_n<<","<<_m<<")"; }
1921 };
1922 
1924 template <typename E> EIGEN_STRONG_INLINE
1925 reshape_expr<E> const reshape(E const & u, index_t n, index_t m)
1926 { return reshape_expr<E>(u, n, m); }
1927 
1928 template<class E>
1929 class replicate_expr : public _expr<replicate_expr<E> >
1930 {
1931 public:
1932  typedef typename E::Scalar Scalar;
1933  enum {ScalarValued = 0, Space = E::Space, ColBlocks= E::ColBlocks};
1934 private:
1935  typename E::Nested_t _u;
1936  index_t _n, _m;
1937  mutable gsMatrix<Scalar> tmp;
1938 
1939 public:
1940 
1941  //the replicate is done nxm times
1942  replicate_expr(_expr<E> const& u, index_t n, index_t m) : _u(u), _n(n), _m(m)
1943  {
1944  }
1945 
1946  auto eval(const index_t k) const -> decltype(tmp.replicate(_n,_m))
1947  {
1948  tmp = _u.eval(k);
1949  return tmp.replicate(_n,_m);
1950  }
1951 
1952  index_t rows() const { return _n*_u.rows(); }
1953  index_t cols() const { return _m*_u.cols(); }
1954 
1955  void parse(gsExprHelper<Scalar> & evList) const
1956  { _u.parse(evList); }
1957 
1958  const gsFeSpace<Scalar> & rowVar() const { return _u.rowVar(); }
1959  const gsFeSpace<Scalar> & colVar() const { return _u.colVar(); }
1960  index_t cardinality_impl() const { return _u.cardinality_impl(); }
1961 
1962  void print(std::ostream &os) const { os << "replicate("; _u.print(os); os<<","<<_n<<","<<_m<<")"; }
1963 };
1964 
1966 template <typename E> EIGEN_STRONG_INLINE
1967 replicate_expr<E> const replicate(E const & u, index_t n, index_t m = 1)
1968 { return replicate_expr<E>(u, n, m); }
1969 
1977 template<class E>
1978 class flat_expr : public _expr<flat_expr<E> >
1979 {
1980 public:
1981  typedef typename E::Scalar Scalar;
1982  enum {ScalarValued = 0, Space = E::Space, ColBlocks= 0}; // to do: ColBlocks
1983 private:
1984  typename E::Nested_t _u;
1985  mutable gsMatrix<Scalar> tmp;
1986 
1987 public:
1988 
1989  flat_expr(_expr<E> const& u) : _u(u)
1990  {
1991  //GISMO_ASSERT( _u.rows()*_u.cols() == _n*_m, "Wrong dimension"); //
1992  }
1993 
1994  const gsMatrix<Scalar> & eval(const index_t k) const
1995  {
1996  tmp = _u.eval(k);
1997  const index_t numActives = _u.cardinality();
1998 
1999  for (index_t i = 0; i<numActives; ++i)
2000  {
2001  tmp(0,2*i+1) += tmp(1,2*i);
2002  std::swap(tmp(1,2*i), tmp(1,2*i+1));
2003  }
2004 
2005  tmp.resize(4,numActives);
2006  tmp.conservativeResize(3,numActives);
2007 
2008  if ( 1==Space )
2009  tmp.transposeInPlace();
2010  else if (2!=Space) // if not colSpan and not rowSpan
2011  tmp.transposeInPlace();
2012 
2013  return tmp;
2014  }
2015 
2016  index_t rows() const { return 1; }
2017  index_t cols() const { return 3; }
2018 
2019  void parse(gsExprHelper<Scalar> & evList) const
2020  { _u.parse(evList); }
2021 
2022  const gsFeSpace<Scalar> & rowVar() const { return _u.rowVar(); }
2023  const gsFeSpace<Scalar> & colVar() const { return _u.colVar(); }
2024  index_t cardinality_impl() const { return _u.cardinality_impl(); }
2025 
2026  void print(std::ostream &os) const { os << "flat("; _u.print(os); os<<")"; }
2027 };
2028 
2030 template <typename E> EIGEN_STRONG_INLINE
2031 flat_expr<E> const flat(E const & u)
2032 { return flat_expr<E>(u); }
2033 
2034 
2035 /*
2036  Expression for the diagonal(s) of a (matrix) expression
2037 */
2038  template<class E>
2039  class diag_expr : public _expr<diag_expr<E> >
2040  {
2041  public:
2042  typedef typename E::Scalar Scalar;
2043  enum {Space=0, ColBlocks=E::ColBlocks, ScalarValued = 0};
2044  private:
2045  typename E::Nested_t _u;
2046  mutable gsMatrix<Scalar> res;
2047 
2048  public:
2049  diag_expr(_expr<E> const& u) : _u(u)
2050  {
2051  GISMO_ASSERT(0== _u.cols()%_u.rows(), "Expecting square-block expression, got "
2052  << _u.rows() <<" x "<< _u.cols() );
2053  }
2054 
2055  const gsMatrix<Scalar> & eval(const index_t k) const
2056  {
2057  // Assume mat ??
2058  MatExprType tmp = _u.eval(k);
2059  const index_t cb = _u.rows();
2060  const index_t r = _u.cols() / cb;
2061  res.resize(r, cb);
2062  for (index_t i = 0; i!=r; ++i)
2063  res.row(i) = tmp.middleCols(i*cb,cb).diagonal();
2064  return res;
2065  }
2066 
2067  index_t rows() const { return _u.cols() / _u.rows(); }
2068  index_t cols() const { return _u.rows(); }
2069 
2070  void parse(gsExprHelper<Scalar> & evList) const
2071  { _u.parse(evList); }
2072 
2073  const gsFeSpace<Scalar> & rowVar() const { return _u.rowVar(); }
2074  const gsFeSpace<Scalar> & colVar() const { return _u.colVar(); }
2075 
2076 
2077  void print(std::ostream &os) const { os << "diag("; _u.print(os); os<<")"; }
2078  };
2079 
2081 template <typename E> EIGEN_STRONG_INLINE
2082 diag_expr<E> const diagonal(E const & u)
2083 { return diag_expr<E>(u); }
2084 
2085 #define GISMO_EXPR_VECTOR_EXPRESSION(name, mname, isSv) \
2086  template<class E> class name##_##expr : public _expr<name##_##expr<E> > { \
2087  typename E::Nested_t _u; \
2088  public: \
2089  typedef typename E::Scalar Scalar; \
2090  enum {Space= E::Space, ScalarValued= isSv, ColBlocks= E::ColBlocks}; \
2091  name##_##expr(_expr<E> const& u) : _u(u) { } \
2092  mutable Temporary_t tmp; \
2093  const Temporary_t & eval(const index_t k) const { \
2094  tmp = _u.eval(k).mname(); return tmp; } \
2095  index_t rows() const { return isSv ? 0 : _u.rows(); } \
2096  index_t cols() const { return isSv ? 0 : _u.cols(); } \
2097  void parse(gsExprHelper<Scalar> & evList) const { _u.parse(evList); } \
2098  const gsFeSpace<Scalar> & rowVar() const {return gsNullExpr<Scalar>::get();} \
2099  const gsFeSpace<Scalar> & colVar() const {return gsNullExpr<Scalar>::get();} \
2100  void print(std::ostream &os) const \
2101  { os << #name <<"("; _u.print(os); os <<")"; } \
2102  };
2103 
2105 GISMO_EXPR_VECTOR_EXPRESSION(norm,norm,1);
2107 GISMO_EXPR_VECTOR_EXPRESSION(sqNorm,squaredNorm,1);
2109 GISMO_EXPR_VECTOR_EXPRESSION(normalized,normalized,0);
2112 // GISMO_EXPR_VECTOR_EXPRESSION(cwSqr,array().square,0)
2113 // GISMO_EXPR_VECTOR_EXPRESSION(sum,array().sum,1)
2114 // GISMO_EXPR_VECTOR_EXPRESSION(sqrt,array().sqrt,0)
2115 //GISMO_EXPR_VECTOR_EXPRESSION(abs,array().abs,0)
2116 
2117 //Determinant
2118 GISMO_EXPR_VECTOR_EXPRESSION(det,determinant,1);
2119 
2120 //GISMO_EXPR_VECTOR_EXPRESSION(replicate,replicate,0);
2121 
2122 #undef GISMO_EXPR_VECTOR_EXPRESSION
2123 
2127 template<class E>
2128 class asdiag_expr : public _expr<asdiag_expr<E> >
2129 {
2130 public:
2131  typedef typename E::Scalar Scalar;
2132 private:
2133  typename E::Nested_t _u;
2134  mutable gsMatrix<Scalar> res;
2135 
2136 public:
2137  enum{Space = E::Space, ScalarValued= 0, ColBlocks= E::ColBlocks};
2138 
2139  asdiag_expr(_expr<E> const& u) : _u(u) { }
2140 
2141 public:
2142 
2143  const gsMatrix<Scalar> & eval(const index_t k) const
2144  {
2145  auto m = _u.eval(k);
2146  const index_t r = m.rows();
2147  const index_t c = m.cols();
2148  res.resize(r,r*c);
2149  for (index_t i = 0; i!=c; ++i)
2150  res.middleCols(i*r,r) = m.col(i).asDiagonal();
2151  return res;
2152  }
2153 
2154  const gsFeSpace<Scalar> & rowVar() const { return _u.rowVar(); }
2155  const gsFeSpace<Scalar> & colVar() const { return _u.colVar(); }
2156 
2157  index_t cardinality_impl() const { return _u.cardinality_impl(); }
2158 
2159  index_t rows() const { return _u.rows(); }
2160  index_t cols() const { return _u.rows() * _u.cols(); }
2161  void parse(gsExprHelper<Scalar> & evList) const
2162  { _u.parse(evList); }
2163 
2164  void print(std::ostream &os) const { os << "diag("; _u.print(os); os <<")";}
2165 };
2166 
2167 // Takes the max of a vector
2168 template<class E>
2169 class max_expr : public _expr<max_expr<E> >
2170 {
2171 public:
2172  typedef typename E::Scalar Scalar;
2173  enum {ScalarValued = 0, Space = E::Space, ColBlocks = 1};
2174 private:
2175  typename E::Nested_t _u;
2176  mutable gsMatrix<Scalar> tmp;
2177  mutable gsMatrix<Scalar> res;
2178 
2179 public:
2180 
2181  max_expr(_expr<E> const& u) : _u(u)
2182  {
2183  //GISMO_ASSERT( _u.rows()*_u.cols() == _n*_m, "Wrong dimension"); //
2184  }
2185 
2186  const gsMatrix<Scalar> & eval(const index_t k) const {return eval_impl(_u,k); }
2187 
2188  index_t rows() const { return 1; }
2189  index_t cols() const { return 1; }
2190  void setFlag() const { _u.setFlag(); }
2191 
2192  void parse(gsExprHelper<Scalar> & evList) const
2193  { _u.parse(evList); }
2194 
2195  const gsFeSpace<Scalar> & rowVar() const { return _u.rowVar(); }
2196  const gsFeSpace<Scalar> & colVar() const { return _u.colVar(); }
2197  index_t cardinality_impl() const { return _u.cardinality_impl(); }
2198 
2199  void print(std::ostream &os) const { os << "max("; _u.print(os); os<<")"; }
2200 private:
2201  template<class U> inline
2202  typename util::enable_if< util::is_same<U,gsFeSpace<Scalar> >::value, const gsMatrix<Scalar> & >::type
2203  eval_impl(const U & u, const index_t k) const
2204  {
2205  tmp = u.eval(k);
2206 
2207  res.resize(1,u.cardinality());
2208  if (E::ColBlocks)
2209  for (index_t c=0; c!=_u.cardinality(); c++)
2210  res(0,c) = tmp.block(0,c*u.cols(),u.rows(),u.cols()).maxCoeff();
2211  else
2212  for (index_t c=0; c!=_u.rows(); c++)
2213  res(0,c) = tmp.block(c*u.rows(),0,u.rows(),u.cols()).maxCoeff();
2214  return res;
2215  }
2216 
2217 
2218  template<class U> inline
2219  typename util::enable_if< !util::is_same<U,gsFeSpace<Scalar> >::value, const gsMatrix<Scalar> & >::type
2220  eval_impl(const U & u, const index_t k) const
2221  {
2222  res = u.eval(k).colwise().maxCoeff();
2223  return res;
2224  }
2225 };
2226 
2227 template<class E>
2228 class rowsum_expr : public _expr<rowsum_expr<E> >
2229 {
2230 public:
2231  typedef typename E::Scalar Scalar;
2232  enum {ScalarValued = 0, Space = E::Space, ColBlocks = E::ColBlocks};
2233 private:
2234  typename E::Nested_t _u;
2235  mutable gsMatrix<Scalar> tmp;
2236 
2237 public:
2238 
2239  rowsum_expr(_expr<E> const& u) : _u(u)
2240  {
2241  //GISMO_ASSERT( _u.rows()*_u.cols() == _n*_m, "Wrong dimension"); //
2242  }
2243 
2244  const gsMatrix<Scalar> & eval(const index_t k) const
2245  {
2246  tmp = _u.eval(k).rowwise().sum();
2247  return tmp;
2248  }
2249 
2250  index_t rows() const { return _u.rows(); }
2251  index_t cols() const { return 1; }
2252  void setFlag() const { _u.setFlag(); }
2253 
2254  void parse(gsExprHelper<Scalar> & evList) const
2255  { _u.parse(evList); }
2256 
2257  const gsFeSpace<Scalar> & rowVar() const { return _u.rowVar(); }
2258  const gsFeSpace<Scalar> & colVar() const { return _u.colVar(); }
2259  index_t cardinality_impl() const { return _u.cardinality_impl(); }
2260 
2261  void print(std::ostream &os) const { os << "rowsum("; _u.print(os); os<<")"; }
2262 };
2263 
2264 template<class E>
2265 class colsum_expr : public _expr<colsum_expr<E> >
2266 {
2267 public:
2268  typedef typename E::Scalar Scalar;
2269  enum {ScalarValued = 0, Space = E::Space, ColBlocks = E::ColBlocks};
2270 private:
2271  typename E::Nested_t _u;
2272  mutable gsMatrix<Scalar> tmp;
2273 
2274 public:
2275 
2276  colsum_expr(_expr<E> const& u) : _u(u)
2277  {
2278  //GISMO_ASSERT( _u.rows()*_u.cols() == _n*_m, "Wrong dimension"); //
2279  }
2280 
2281  const gsMatrix<Scalar> & eval(const index_t k) const
2282  {
2283  tmp = _u.eval(k).colwise().sum();
2284  return tmp;
2285  }
2286 
2287  index_t rows() const { return _u.rows(); }
2288  index_t cols() const { return 1; }
2289  void setFlag() const { _u.setFlag(); }
2290 
2291  void parse(gsExprHelper<Scalar> & evList) const
2292  { _u.parse(evList); }
2293 
2294  const gsFeSpace<Scalar> & rowVar() const { return _u.rowVar(); }
2295  const gsFeSpace<Scalar> & colVar() const { return _u.colVar(); }
2296  index_t cardinality_impl() const { return _u.cardinality_impl(); }
2297 
2298  void print(std::ostream &os) const { os << "colsum("; _u.print(os); os<<")"; }
2299 };
2300 
2304 class idMat_expr : public _expr<idMat_expr >
2305 {
2306 public:
2307  typedef real_t Scalar;
2308  enum {Space = 0, ScalarValued = 0, ColBlocks = 0};
2309 private:
2310  index_t _dim;
2311 
2312 public:
2313  idMat_expr(const index_t dim) : _dim(dim) { }
2314 
2315 public:
2316 
2318  {
2319  return gsMatrix<Scalar>::Identity(_dim,_dim);
2320  }
2321 
2322  index_t rows() const { return _dim; }
2323  index_t cols() const { return _dim; }
2324  void parse(gsExprHelper<Scalar> & ) const { }
2325 
2326  const gsFeSpace<Scalar> & rowVar() const {return gsNullExpr<Scalar>::get();}
2327  const gsFeSpace<Scalar> & colVar() const {return gsNullExpr<Scalar>::get();}
2328 
2329  void print(std::ostream &os) const { os << "id("<<_dim <<")";}
2330 };
2331 
2335 class constMat_expr : public _expr<constMat_expr >
2336 {
2337 public:
2338  typedef real_t Scalar;
2339  enum {Space = 0, ScalarValued = 0, ColBlocks = 0};
2340 private:
2341  gsMatrix<Scalar> _mat;
2342 
2343 public:
2344  constMat_expr(const gsMatrix<Scalar> mat) : _mat(mat) { }
2345 
2346 public:
2347 
2348  gsMatrix<Scalar> eval(const index_t) const
2349  {
2350  return _mat;
2351  }
2352 
2353  index_t rows() const { return _mat.rows(); }
2354  index_t cols() const { return _mat.cols(); }
2355  void parse(gsExprHelper<Scalar> & ) const { }
2356 
2357  const gsFeSpace<Scalar> & rowVar() const {return gsNullExpr<Scalar>::get();}
2358  const gsFeSpace<Scalar> & colVar() const {return gsNullExpr<Scalar>::get();}
2359 
2360  void print(std::ostream &os) const { os << "constMat";}
2361 };
2362 
2366 template<class E>
2367 class sign_expr : public _expr<sign_expr<E> >
2368 {
2369  typename E::Nested_t _u;
2370  typename E::Scalar _tol;
2371 public:
2372  typedef typename E::Scalar Scalar;
2373  enum {ScalarValued = 1, Space = E::Space, ColBlocks= 0};
2374 
2375  sign_expr(_expr<E> const& u, Scalar tolerance = 0.0) : _u(u),_tol(tolerance){
2376  GISMO_ASSERT( _tol >= 0, "Tolerance for sign_expr should be a positive number.");
2377  }
2378 
2379  Scalar eval(const index_t k) const
2380  {
2381  const Scalar v = _u.val().eval(k);
2382  return ( v>_tol ? 1 : ( v<-_tol ? -1 : 0 ) );
2383  }
2384 
2385  static index_t rows() { return 0; }
2386  static index_t cols() { return 0; }
2387 
2388  void parse(gsExprHelper<Scalar> & el) const
2389  { _u.parse(el); }
2390 
2391  static bool isScalar() { return true; }
2392 
2393  const gsFeSpace<Scalar> & rowVar() const {return gsNullExpr<Scalar>::get();}
2394  const gsFeSpace<Scalar> & colVar() const {return gsNullExpr<Scalar>::get();}
2395 
2396  void print(std::ostream &os) const { os<<"sgn("; _u.print(os); os <<")"; }
2397 };
2398 
2402 template<class E>
2403 class exp_expr : public _expr<exp_expr<E> >
2404 {
2405  typename E::Nested_t _u;
2406  public:
2407  typedef typename E::Scalar Scalar;
2408  enum {ScalarValued = 1, Space = E::Space, ColBlocks= 0};
2409 
2410  exp_expr(_expr<E> const& u) : _u(u) { }
2411 
2412  Scalar eval(const index_t k) const
2413  {
2414  const Scalar v = _u.val().eval(k);
2415  return math::exp(v);
2416  }
2417 
2418  static index_t rows() { return 0; }
2419  static index_t cols() { return 0; }
2420 
2421  void parse(gsExprHelper<Scalar> & el) const
2422  { _u.parse(el); }
2423 
2424  static bool isScalar() { return true; }
2425 
2426  const gsFeSpace<Scalar> & rowVar() const {return gsNullExpr<Scalar>::get();}
2427  const gsFeSpace<Scalar> & colVar() const {return gsNullExpr<Scalar>::get();}
2428 
2429  void print(std::ostream &os) const { os<<"exp("; _u.print(os); os <<")"; }
2430 };
2431 
2435 template<class E>
2436 class ppart_expr : public _expr<ppart_expr<E> >
2437 {
2438 public:
2439  typedef typename E::Scalar Scalar;
2440  enum {ScalarValued = E::ScalarValued, Space = E::Space, ColBlocks= E::ColBlocks};
2441 private:
2442  typename E::Nested_t _u;
2443  mutable gsMatrix<Scalar> res;
2444 public:
2445 
2446  ppart_expr(_expr<E> const& u) : _u(u) { }
2447 
2448  const gsMatrix<Scalar> & eval(index_t k) const
2449  {
2450  res = _u.eval(k).cwiseMax(0.0); // component-wise maximum with zero
2451  return res;
2452  }
2453 
2454 
2455  index_t rows() const { return _u.rows(); }
2456  index_t cols() const { return _u.cols(); }
2457 
2458  void parse(gsExprHelper<Scalar> & el) const
2459  { _u.parse(el); }
2460 
2461  const gsFeSpace<Scalar> & rowVar() const {return _u.rowVar();}
2462  const gsFeSpace<Scalar> & colVar() const {return _u.colVar();}
2463 
2464  void print(std::ostream &os) const { os<<"posPart("; _u.print(os); os <<")"; }
2465 };
2466 
2470 template<class E>
2471 class ppartval_expr : public _expr<ppartval_expr<E> >
2472 {
2473  typename E::Nested_t _u;
2474  public:
2475  typedef typename E::Scalar Scalar;
2476  enum {ScalarValued = 1, Space = 0, ColBlocks= 0};
2477  mutable Scalar res;
2478  public:
2479 
2480  ppartval_expr(_expr<E> const& u) : _u(u) { }
2481 
2482  Scalar & eval(index_t k) const
2483  {
2484  res = std::max(0.0,_u.eval(k));
2485  return res; // component-wise maximum with zero
2486  }
2487 
2488  index_t rows() const { return 0; }
2489  index_t cols() const { return 0; }
2490 
2491  void parse(gsExprHelper<Scalar> & evList) const
2492  { _u.parse(evList); }
2493 
2494  const gsFeSpace<Scalar> & rowVar() const {return gsNullExpr<Scalar>::get();}
2495  const gsFeSpace<Scalar> & colVar() const {return gsNullExpr<Scalar>::get();}
2496 
2497  void print(std::ostream &os) const { os<<"posPart("; _u.print(os); os <<")"; }
2498 };
2499 
2503 template<class E>
2504 class pow_expr : public _expr<pow_expr<E> >
2505 {
2506  typename E::Nested_t _u;
2507 
2508 public:
2509  typedef typename E::Scalar Scalar;
2510  enum {ScalarValued = 1, Space = E::Space, ColBlocks= E::ColBlocks};
2511 
2512  Scalar _q;// power
2513 
2514  pow_expr(_expr<E> const& u, Scalar q) : _u(u), _q(q) { }
2515 
2516  Scalar eval(const index_t k) const
2517  {
2518  const Scalar v = _u.val().eval(k);
2519  return math::pow(v,_q);
2520  }
2521 
2522  static index_t rows() { return 0; }
2523  static index_t cols() { return 0; }
2524 
2525  void parse(gsExprHelper<Scalar> & el) const
2526  { _u.parse(el); }
2527 
2528  static bool isScalar() { return true; }
2529 
2530  const gsFeSpace<Scalar> & rowVar() const {return gsNullExpr<Scalar>::get();}
2531  const gsFeSpace<Scalar> & colVar() const {return gsNullExpr<Scalar>::get();}
2532 
2533  void print(std::ostream &os) const { os<<"pow("; _u.print(os); os <<")"; }
2534 };
2535 
2541 template <typename E1, typename E2>
2542 class matrix_by_space_expr : public _expr<matrix_by_space_expr<E1,E2> >
2543 {
2544 public:
2545  typedef typename E1::Scalar Scalar;
2546  enum {ScalarValued = 0, ColBlocks = 1};
2547  enum {Space = E2::Space};
2548 
2549 private:
2550  typename E1::Nested_t _u;
2551  typename E2::Nested_t _v;
2552  mutable gsMatrix<Scalar> res;
2553 
2554 public:
2555  matrix_by_space_expr(E1 const& u, E2 const& v) : _u(u), _v(v) { }
2556 
2557 
2558  // choose if ColBlocks
2559  const gsMatrix<Scalar> & eval(const index_t k) const
2560  {
2561  const index_t r = _u.rows();
2562  const index_t N = _v.cols() / (r*r);
2563 
2564  const auto uEv = _u.eval(k);
2565  const auto vEv = _v.eval(k);
2566 
2567  res.resize(r, N*r*r);
2568  // gsDebugVar(res.cols());
2569  for (index_t s = 0; s!=r; ++s)
2570  for (index_t i = 0; i!=N; ++i)
2571  {
2572  res.middleCols((s*N + i)*r,r).noalias() =
2573  uEv.col(s) * vEv.middleCols((s*N + i)*r,r).row(s);
2574  //uEv*vEv.middleCols((s*N + i)*r,r);
2575  }
2576  //meaning: [Jg Jg Jg] * Jb ..
2577  return res;
2578  }
2579 
2580  index_t rows() const { return _u.cols(); }
2581  index_t cols() const { return _v.cols(); }
2582 
2583  void parse(gsExprHelper<Scalar> & evList) const
2584  { _u.parse(evList); _v.parse(evList); }
2585 
2586  const gsFeSpace<Scalar> & rowVar() const { return _v.rowVar(); }
2587  const gsFeSpace<Scalar> & colVar() const { return _v.colVar(); }
2588 
2589  void print(std::ostream &os) const { os << "matrix_by_space("; _u.print(os); os<<")"; }
2590 };
2591 
2597 template <typename E1, typename E2>
2598 class matrix_by_space_expr_tr : public _expr<matrix_by_space_expr_tr<E1,E2> >
2599 {
2600 public:
2601  typedef typename E1::Scalar Scalar;
2602  enum {ScalarValued = 0, ColBlocks = 1};
2603  enum {Space = E2::Space};
2604 
2605 private:
2606  typename E1::Nested_t _u;
2607  typename E2::Nested_t _v;
2608  mutable gsMatrix<Scalar> res;
2609 
2610 public:
2611  matrix_by_space_expr_tr(E1 const& u, E2 const& v) : _u(u), _v(v) { }
2612 
2613 
2614  // choose if ColBlocks
2615  const gsMatrix<Scalar> & eval(const index_t k) const
2616  {
2617  const index_t r = _u.rows();
2618  const index_t N = _v.cols() / (r*r);
2619 
2620  const auto uEv = _u.eval(k);
2621  const auto vEv = _v.eval(k);
2622 
2623  res.resize(r, N*r*r);
2624  for (index_t s = 0; s!=r; ++s)
2625  for (index_t i = 0; i!=N; ++i)
2626  {
2627  res.middleCols((s*N + i)*r,r).noalias() =
2628  uEv.transpose()*vEv.middleCols((s*N + i)*r,r).transpose();
2629  }
2630  //meaning: [Jg Jg Jg] * Jb ..
2631  return res;
2632  }
2633 
2634  index_t rows() const { return _u.cols(); }
2635  index_t cols() const { return _v.cols(); }
2636 
2637  void parse(gsExprHelper<Scalar> & evList) const
2638  { _u.parse(evList); _v.parse(evList); }
2639 
2640  const gsFeSpace<Scalar> & rowVar() const { return _v.rowVar(); }
2641  const gsFeSpace<Scalar> & colVar() const { return _v.colVar(); }
2642 
2643  void print(std::ostream &os) const { os << "matrix_by_space_tr("; _u.print(os); os<<")"; }
2644 };
2645 
2646 /*
2647  Adaptor for scalar-valued expression
2648 */
2649 template<class E>
2650 class value_expr : public _expr<value_expr<E> >
2651 {
2652  typename E::Nested_t _u;
2653 
2654 public:
2655  typedef typename E::Scalar Scalar;
2656  value_expr(_expr<E> const& u) : _u(u)
2657  {
2658  // rows/cols not known at construction time
2659  //GISMO_ASSERT(u.rows()*u.cols()<=1, "Expression\n"<<u<<"is not a scalar.");
2660  }
2661 
2662 public:
2663  enum {Space= 0, ScalarValued= 1, ColBlocks= 0};
2664 
2665  Scalar eval(const index_t k) const { return eval_impl(_u,k); }
2666 
2667  // enables expr.val().val()
2668  inline value_expr<E> val() const { return *this; }
2669  index_t rows() const { return 0; }
2670  index_t cols() const { return 0; }
2671  void parse(gsExprHelper<Scalar> & evList) const
2672  { _u.parse(evList); }
2673 
2674  static bool isScalar() { return true; }
2675 
2676  const gsFeSpace<Scalar> & rowVar() const {return gsNullExpr<Scalar>::get();}
2677  const gsFeSpace<Scalar> & colVar() const {return gsNullExpr<Scalar>::get();}
2678 
2679  void print(std::ostream &os) const { _u.print(os); }
2680 
2681  // Math functions. eg
2682  // sqrt_mexpr<T> sqrt() { return sqrt_mexpr<T>(*this); }
2683 private:
2684  template<class U> static inline
2685  typename util::enable_if<U::ScalarValued,Scalar>::type
2686  eval_impl(const U & u, const index_t k) { return u.eval(k); }
2687 
2688  template<class U> static inline
2689  typename util::enable_if<!U::ScalarValued,Scalar>::type
2690  eval_impl(const U & u, const index_t k) { return u.eval(k).value(); }
2691 };
2692 
2693 template<class E>
2694 class abs_expr : public _expr<abs_expr<E> >
2695 {
2696  typename E::Nested_t _u;
2697 
2698 public:
2699  typedef typename E::Scalar Scalar;
2700  explicit abs_expr(_expr<E> const& u) : _u(u) { }
2701 
2702 public:
2703  enum {Space= 0, ScalarValued= 1, ColBlocks= 0};
2704 
2705  Scalar eval(const index_t k) const { return abs_expr::eval_impl(_u,k); }
2706 
2707  index_t rows() const { return _u.rows(); }
2708  index_t cols() const { return _u.cols(); }
2709  void parse(gsExprHelper<Scalar> & evList) const
2710  { _u.parse(evList); }
2711 
2712  static bool isScalar() { return true; }
2713 
2714  const gsFeSpace<Scalar> & rowVar() const {return gsNullExpr<Scalar>::get();}
2715  const gsFeSpace<Scalar> & colVar() const {return gsNullExpr<Scalar>::get();}
2716 
2717  void print(std::ostream &os) const { _u.print(os); }
2718 
2719  // Math functions. eg
2720  // sqrt_mexpr<T> sqrt() { return sqrt_mexpr<T>(*this); }
2721 private:
2722  template<class U> static inline
2723  typename util::enable_if<U::ScalarValued,Scalar>::type
2724  eval_impl(const U & u, const index_t k) {return math::abs(u.eval(k)); }
2725  template<class U> static inline
2726  typename util::enable_if<!U::ScalarValued,gsMatrix<Scalar> >::type
2727  eval_impl(const U & u, const index_t k) { return u.eval(k).cwiseAbs(); }
2728 };
2729 
2730 /*
2731  Expression for the gradient of a finite element variable
2732 
2733  Transposed gradient vectors are returned as a matrix
2734 */
2735 template<class E>
2736 class grad_expr : public _expr<grad_expr<E> >
2737 {
2738  typename E::Nested_t _u;
2739 public:
2740  enum {Space = E::Space, ScalarValued= 0, ColBlocks= 0};
2741 
2742  typedef typename E::Scalar Scalar;
2743  mutable gsMatrix<Scalar> tmp;
2744 
2745  grad_expr(const E & u) : _u(u)
2746  { GISMO_ASSERT(1==u.dim(),"grad(.) requires 1D variable, use jac(.) instead.");}
2747 
2748  const gsMatrix<Scalar> & eval(const index_t k) const
2749  {
2750  // Assumes: derivatives are in _u.data().values[1]
2751  // gsExprHelper acounts for compositions/physical expressions
2752  // so that derivs are directly computed
2753  tmp = _u.data().values[1].reshapeCol(k, cols(), cardinality_impl()).transpose();
2754  return tmp;
2755  }
2756 
2757  index_t rows() const { return 1 /*==u.dim()*/; }
2758 
2759  index_t cols() const { return _u.source().domainDim(); }
2760 
2761  index_t cardinality_impl() const
2762  { return _u.data().values[1].rows() / cols(); }
2763 
2764  void parse(gsExprHelper<Scalar> & evList) const
2765  {
2766  evList.add(_u);
2767  _u.data().flags |= NEED_GRAD;
2768  }
2769 
2770  const gsFeSpace<Scalar> & rowVar() const { return _u.rowVar(); }
2771  const gsFeSpace<Scalar> & colVar() const
2772  {return gsNullExpr<Scalar>::get();}
2773 
2774  void print(std::ostream &os) const { os << "\u2207("; _u.print(os); os <<")"; }
2775 private:
2776 
2777  template<class U> static inline
2778  typename util::enable_if<util::is_same<U,gsComposition<Scalar> >::value,
2779  const gsMatrix<Scalar> &>::type
2780  eval_impl(const U & u, const index_t k)
2781  {
2782  return u.eval(k);
2783  }
2784 };
2785 
2786 /*
2787  \brief Expression for the gradient of a finite element variable
2788 
2789  Transposed gradient vectors are returned as a matrix.
2790  This specialization is for a gsFeSolution object
2791 */
2792 
2793 template<class T>
2794 class grad_expr<gsFeSolution<T> > : public _expr<grad_expr<gsFeSolution<T> > >
2795 {
2796 protected:
2797  const gsFeSolution<T> _u;
2798 
2799 public:
2800  typedef T Scalar;
2801  enum {Space = 0, ScalarValued= 0, ColBlocks= 0};
2802 
2803  explicit grad_expr(const gsFeSolution<T> & u) : _u(u) { }
2804 
2805  mutable gsMatrix<T> res;
2806  const gsMatrix<T> & eval(index_t k) const
2807  {
2808  GISMO_ASSERT(1==_u.data().actives.cols(), "Single actives expected");
2809 
2810  res.setZero(_u.dim(), _u.parDim());
2811  const gsDofMapper & map = _u.mapper();
2812  for (index_t c = 0; c!= _u.dim(); c++)
2813  {
2814  for (index_t i = 0; i!=_u.data().actives.size(); ++i)
2815  {
2816  const index_t ii = map.index(_u.data().actives.at(i), _u.data().patchId,c);
2817  if ( map.is_free_index(ii) ) // DoF value is in the solVector
2818  {
2819  res.row(c) += _u.coefs().at(ii) *
2820  _u.data().values[1].col(k).segment(i*_u.parDim(), _u.parDim()).transpose();
2821  }
2822  else
2823  {
2824  res.row(c) +=
2825  _u.fixedPart().at( map.global_to_bindex(ii) ) *
2826  _u.data().values[1].col(k).segment(i*_u.parDim(), _u.parDim()).transpose();
2827  }
2828  }
2829  }
2830  return res;
2831  }
2832 
2833  index_t rows() const {return _u.dim();}
2834  index_t cols() const {return _u.parDim(); }
2835 
2836  const gsFeSpace<Scalar> & rowVar() const
2837  {return gsNullExpr<Scalar>::get();}
2838  const gsFeSpace<Scalar> & colVar() const
2839  {return gsNullExpr<Scalar>::get();}
2840 
2841  void parse(gsExprHelper<Scalar> & evList) const
2842  {
2843  _u.parse(evList); // add symbol
2844  evList.add(_u.space());
2845  _u.data().flags |= NEED_GRAD|NEED_ACTIVE; // define flags
2846  }
2847 
2848  void print(std::ostream &os) const { os << "\u2207(s)"; }
2849 };
2850 
2851 /*
2852  Expression for the derivative of the jacobian of a spline geometry map,
2853  with respect to the coordinate c.
2854 
2855  It returns a matrix with the gradient of u in row d.
2856 */
2857 template<class E>
2858 class dJacdc_expr : public _expr<dJacdc_expr<E> >
2859 {
2860  typename E::Nested_t _u;
2861 public:
2862  enum{ Space = E::Space, ScalarValued = 0, ColBlocks = (1==E::Space?1:0)};
2863 
2864  typedef typename E::Scalar Scalar;
2865 
2866  mutable gsMatrix<Scalar> res;
2867  index_t _c;
2868 
2869  dJacdc_expr(const E & u, index_t c) : _u(u), _c(c)
2870  { GISMO_ASSERT(1==u.dim(),"grad(.) requires 1D variable, use jac(.) instead.");}
2871 
2872  const gsMatrix<Scalar> & eval(const index_t k) const
2873  {
2874  index_t dd = _u.source().domainDim();
2875  index_t n = _u.rows();
2876  res.setZero(dd, dd*n);
2877 
2878  gsMatrix<Scalar> grad = _u.data().values[1].reshapeCol(k, dd, n);
2879  for(index_t i = 0; i < n; i++){
2880  res.row(_c).segment(i*dd,dd) = grad.col(i);
2881  }
2882  return res;
2883  }
2884 
2885  index_t rows() const { return _u.source().domainDim(); }
2886 
2887  index_t cols() const { return _u.source().domainDim()*_u.rows(); }
2888 
2889  void parse(gsExprHelper<Scalar> & evList) const
2890  {
2891  evList.add(_u);
2892  _u.data().flags |= NEED_GRAD;
2893  }
2894 
2895  const gsFeSpace<Scalar> & rowVar() const { return _u.rowVar(); }
2896  const gsFeSpace<Scalar> & colVar() const
2897  {return gsNullExpr<Scalar>::get();}
2898 
2899  void print(std::ostream &os) const { os << "dJacdc("; _u.print(os); os <<")"; }
2900 };
2901 
2902 
2903 /*
2904  Expression for the nabla (\f$\nabla\f$) of a finite element variable,
2905 */
2906 template<class T>
2907 class nabla_expr : public _expr<nabla_expr<T> >
2908 {
2909  typename gsFeVariable<T>::Nested_t u;
2910 
2911 public:
2912  typedef T Scalar;
2913  enum{Space = 1};
2914 
2915  /* // todo
2916  nabla_expr(const gsGeometryMap<T> & G)
2917  : m_data(G.data()) { }
2918  */
2919 
2920  nabla_expr(const gsFeVariable<T> & _u) : u(_u)
2921  {
2922  //GISMO_ASSERT(u.parDim()==u.dim(),"nabla(.) requires tarDim==parDim:"
2923  // << u.parDim() <<"!="<< u.dim() <<"\n" );
2924  }
2925 
2926  mutable gsMatrix<Scalar> res;
2927 
2928  const gsMatrix<Scalar> & eval(const index_t k) const
2929  {
2930  const index_t d = u.cols();
2931  const index_t n = rows() / d;
2932  res.setZero(rows(), d);
2933 
2934  for (index_t i = 0; i!=d; ++i)
2935  res.col(i).segment(i*n,n) = u.data().values[1].reshapeCol(k, d, n).row(i);
2936  return res;
2937  }
2938 
2939  index_t rows() const { return u.rows(); }
2940  index_t cols() const { return u.cols(); }
2941 
2942  void parse(gsExprHelper<Scalar> & evList) const
2943  {
2944  evList.add(u);
2945  u.data().flags |= NEED_GRAD;
2946  }
2947 
2948  const gsFeSpace<T> & rowVar() const { return u.rowVar(); }
2949  const gsFeSpace<T> & colVar() const
2950  {return gsNullExpr<Scalar>::get();}
2951 
2952  void print(std::ostream &os) const { os << "nabla("; u.print(os); os <<")"; }
2953 };
2954 
2955 /*
2956  Expression for the nabla2 (\f$\nabla^2\f$ or Del2) of a finite element variable,
2957  see also https://en.wikipedia.org/wiki/Del
2958 
2959  Transposed pure second derivatives are returned as a matrix
2960 */
2961 template<class T>
2962 class nabla2_expr : public _expr<nabla2_expr<T> >
2963 {
2964  typename gsFeVariable<T>::Nested_t u;
2965 
2966 public:
2967  typedef T Scalar;
2968  enum{Space = 1};
2969 
2970  /* // todo
2971  nabla2_expr(const gsGeometryMap<T> & G)
2972  : m_data(G.data()) { }
2973  */
2974 
2975  nabla2_expr(const gsFeVariable<T> & _u)
2976  : u(_u)
2977  { }
2978 
2979  MatExprType eval(const index_t k) const
2980  {
2981  // numActive x parDim
2982  return u.data().values[2]
2983  .reShapeCol(k, u.data().values[2].rows()/u.cSize(), u.cSize() )
2984  .topRows(u.parDim()).transpose();
2985  }
2986 
2987  index_t rows() const { return u.rows(); }
2988  index_t cols() const { return u.parDim(); }
2989 
2990  void parse(gsExprHelper<Scalar> & evList) const
2991  {
2992  evList.add(u);
2993  u.data().flags |= NEED_DERIV2;
2994  }
2995 
2996  const gsFeSpace<T> & rowVar() const { return u.rowVar(); }
2997  const gsFeSpace<T> & colVar() const
2998  {return gsNullExpr<T>::get();}
2999 };
3000 
3002 template<class T>
3003 nabla2_expr<T> nabla2(const gsFeVariable<T> & u) { return nabla2_expr<T>(u); }
3004 // #define lapl(x) nabla2(x).sum() // assume tarDim==1
3005 
3010 template<class T>
3011 class onormal_expr : public _expr<onormal_expr<T> >
3012 {
3013  typename gsGeometryMap<T>::Nested_t _G;
3014 
3015 public:
3016  typedef T Scalar;
3017  enum {Space = 0, ScalarValued= 0, ColBlocks= 0};
3018 
3019  explicit onormal_expr(const gsGeometryMap<T> & G) : _G(G) { }
3020 
3021  auto eval(const index_t k) const -> decltype(_G.data().outNormals.col(k))
3022  { return _G.data().outNormals.col(k); }
3023 
3024  index_t rows() const { return _G.source().targetDim(); }
3025  index_t cols() const { return 1; }
3026 
3027  const gsFeSpace<T> & rowVar() const {return gsNullExpr<T>::get();}
3028  const gsFeSpace<T> & colVar() const {return gsNullExpr<T>::get();}
3029 
3030  void parse(gsExprHelper<Scalar> & evList) const
3031  {
3032  evList.add(_G);
3033  _G.data().flags |= NEED_OUTER_NORMAL;
3034  }
3035 
3036  void print(std::ostream &os) const { os << "nv("; _G.print(os); os <<")"; }
3037 };
3038 
3043 template<class T>
3044 class normal_expr : public _expr<normal_expr<T> >
3045 {
3046  typename gsGeometryMap<T>::Nested_t _G;
3047 
3048 public:
3049  typedef T Scalar;
3050  enum {Space = 0, ScalarValued= 0, ColBlocks= 0};
3051 
3052  normal_expr(const gsGeometryMap<T> & G) : _G(G)
3053  {
3054  GISMO_ENSURE( _G.source().domainDim()+1 == _G.source().targetDim(), "Surface normal requires codimension 1");
3055  }
3056 
3057  auto eval(const index_t k) const -> decltype(_G.data().normals.col(k))
3058  { return _G.data().normals.col(k); }
3059 
3060  index_t rows() const { return _G.source().targetDim(); }
3061  index_t cols() const { return 1; }
3062 
3063  const gsFeSpace<T> & rowVar() const {return gsNullExpr<T>::get();}
3064  const gsFeSpace<T> & colVar() const {return gsNullExpr<T>::get();}
3065 
3066  void parse(gsExprHelper<Scalar> & evList) const
3067  {
3068  evList.add(_G);
3069  _G.data().flags |= NEED_NORMAL;
3070  }
3071 
3072  void print(std::ostream &os) const { os << "sn("; _G.print(os); os <<")"; }
3073 };
3074 
3079 template<class T>
3080 class tangent_expr : public _expr<tangent_expr<T> >
3081 {
3082  typename gsGeometryMap<T>::Nested_t _G;
3083 
3084 public:
3085  typedef T Scalar;
3086  enum {Space = 0, ScalarValued= 0, ColBlocks= 0};
3087 
3088  tangent_expr(const gsGeometryMap<T> & G) : _G(G) { }
3089 
3090  mutable gsVector<Scalar> res;
3091  const gsVector<Scalar> & eval(const index_t k) const
3092  {
3093  if (_G.targetDim()==2)
3094  {
3095  res = _G.data().outNormals.col(k);//2x1
3096  std::swap( res(0,0), res(1,0) );
3097  res(0,0) *= -1;
3098  return res;
3099  }
3100  else if (_G.targetDim()==3)
3101  {
3102  res.resize(3);
3103  res.col3d(0) = _G.data().normals.col3d(k)
3104  .cross( _G.data().outNormals.col3d(k) );
3105  return res;
3106  }
3107  else
3108  GISMO_ERROR("Function not implemented for dimension"<<_G.targetDim());
3109 
3110  }
3111 
3112  index_t rows() const { return _G.source().targetDim(); }
3113  index_t cols() const { return 1; }
3114 
3115  static const gsFeSpace<Scalar> & rowVar() {return gsNullExpr<Scalar>::get();}
3116  static const gsFeSpace<Scalar> & colVar() {return gsNullExpr<Scalar>::get();}
3117 
3118  void parse(gsExprHelper<Scalar> & evList) const
3119  {
3120  evList.add(_G);
3121  _G.data().flags |= NEED_NORMAL;
3122  _G.data().flags |= NEED_OUTER_NORMAL;
3123  }
3124 
3125  void print(std::ostream &os) const { os << "tv("; _G.print(os); os <<")"; }
3126 };
3127 
3131 template<class E>
3132 class lapl_expr : public _expr<lapl_expr<E> >
3133 {
3134  typename E::Nested_t _u;
3135 
3136 public:
3137  typedef typename E::Scalar Scalar;
3138  enum {Space = E::Space, ScalarValued= 0, ColBlocks= 0};
3139 
3140  lapl_expr(const E & u) : _u(u) { }
3141 
3142  auto eval(const index_t k) const -> decltype(_u.data().laplacians.col(k))
3143  {
3144  // numActive x 1
3145  return _u.data().laplacians.col(k);
3146  //todo: replace by
3147  // NEED_DERIV2
3148  // ..nabla2.sum()
3149  }
3150 
3151  index_t rows() const { return _u.data().laplacians.rows(); }
3152  index_t cols() const { return 1; }
3153 
3154  index_t cardinality_impl() const { return _u.cardinality_impl(); }
3155 
3156  void parse(gsExprHelper<Scalar> & evList) const
3157  {
3158  evList.add(_u);
3159  _u.data().flags |= NEED_LAPLACIAN;
3160  }
3161 
3162  static const gsFeSpace<Scalar> & rowVar() {return E::rowVar();}
3163  static const gsFeSpace<Scalar> & colVar() {return gsNullExpr<Scalar>::get();}
3164 
3165  void print(std::ostream &os) const { os << "\u2206("; _u.print(os); os <<")"; } //or \u0394
3166 };
3167 
3168 /*
3169  Expression for the Laplacian of a finite element solution
3170 */
3171 template<class T>
3172 class lapl_expr<gsFeSolution<T> > : public _expr<lapl_expr<gsFeSolution<T> > >
3173 {
3174 protected:
3175  const gsFeSolution<T> _u;
3176 
3177 public:
3178  typedef T Scalar;
3179  enum {Space = 0, ScalarValued= 0, ColBlocks= 0};
3180 
3181  lapl_expr(const gsFeSolution<T> & u) : _u(u) { }
3182 
3183  mutable gsMatrix<T> res;
3184  const gsMatrix<T> & eval(const index_t k) const
3185  {
3186  GISMO_ASSERT(1==_u.data().actives.cols(), "Single actives expected");
3187 
3188  res.setZero(_u.dim(), 1); // scalar, but per component
3189  const gsDofMapper & map = _u.mapper();
3190 
3191  index_t numActs = _u.data().values[0].rows();
3192  index_t numDers = _u.parDim() * (_u.parDim() + 1) / 2;
3193  gsMatrix<T> deriv2;
3194 
3195  for (index_t c = 0; c!= _u.dim(); c++)
3196  for (index_t i = 0; i!=numActs; ++i)
3197  {
3198  const index_t ii = map.index(_u.data().actives.at(i), _u.data().patchId,c);
3199  deriv2 = _u.data().values[2].block(i*numDers,k,_u.parDim(),1); // this only takes d11, d22, d33 part. For all the derivatives [d11, d22, d33, d12, d13, d23]: col.block(i*numDers,k,numDers,1)
3200  if ( map.is_free_index(ii) ) // DoF value is in the solVector
3201  res.at(c) += _u.coefs().at(ii) * deriv2.sum();
3202  else
3203  res.at(c) +=_u.fixedPart().at( map.global_to_bindex(ii) ) * deriv2.sum();
3204  }
3205  return res;
3206  }
3207 
3208  index_t rows() const { return _u.dim(); }
3209  index_t cols() const { return 1; }
3210 
3211  void parse(gsExprHelper<Scalar> & evList) const
3212  {
3213  evList.add(_u.space());
3214  _u.data().flags |= NEED_ACTIVE | NEED_DERIV2;
3215  }
3216 
3217  const gsFeSpace<Scalar> & rowVar() const {return gsNullExpr<T>::get();}
3218  const gsFeSpace<Scalar> & colVar() const {return gsNullExpr<T>::get();}
3219 
3220  void print(std::ostream &os) const { os << "\u2206(s)"; }
3221 };
3222 
3223 /*
3224  Expression for the (precomputed) second fundamental form of a surface
3225 */
3226 template<class T>
3227 class fform2nd_expr : public _expr<fform2nd_expr<T> >
3228 {
3229  typename gsGeometryMap<T>::Nested_t _G;
3230 public:
3231  typedef T Scalar;
3232  enum {Space = 0, ScalarValued = 0, ColBlocks = 0};
3233 
3234  fform2nd_expr(const gsGeometryMap<T> & G) : _G(G) { }
3235 
3236  const gsAsConstMatrix<Scalar> eval(const index_t k) const
3237  {
3238  return gsAsConstMatrix<Scalar>(_G.data().fundForms.col(k).data(),rows(),cols());
3239  }
3240 
3241  index_t rows() const { return _G.source().domainDim() ; }
3242  index_t cols() const { return _G.source().domainDim() ; }
3243 
3244  void parse(gsExprHelper<Scalar> & evList) const
3245  {
3246  evList.add(_G);
3247  _G.data().flags |= NEED_2ND_FFORM;
3248  }
3249 
3250  const gsFeSpace<Scalar> & rowVar() const {return gsNullExpr<T>::get();}
3251  const gsFeSpace<Scalar> & colVar() const {return gsNullExpr<T>::get();}
3252 
3253  void print(std::ostream &os) const { os << "fform2nd("; _G.print(os); os <<")"; }
3254 };
3255 
3256 /*
3257  Expression for the (precomputed) inverse of the Jacobian matrix of
3258  a geometry map
3259 */
3260 template<class T>
3261 class jacInv_expr : public _expr<jacInv_expr<T> >
3262 {
3263  typename gsGeometryMap<T>::Nested_t _G;
3264 public:
3265  typedef T Scalar;
3266  enum {Space = 0, ScalarValued = 0, ColBlocks = 0};
3267 
3268  jacInv_expr(const gsGeometryMap<T> & G) : _G(G)
3269  {
3270  // Note: for non-square Jacobian matrices, generalized inverse, i.e.: (J^t J)^{-t} J^t
3271  //GISMO_ASSERT(rows() == cols(), "The Jacobian matrix is not square");
3272  }
3273 
3274  MatExprType eval(const index_t k) const { return _G.data().jacInvTr.reshapeCol(k,cols(),rows()).transpose(); }
3275 
3276  index_t rows() const { return _G.source().domainDim(); }
3277  index_t cols() const { return _G.source().targetDim(); }
3278 
3279  void parse(gsExprHelper<Scalar> & evList) const
3280  {
3281  evList.add(_G);
3282  _G.data().flags |= NEED_GRAD_TRANSFORM;
3283  }
3284 
3285  const gsFeSpace<Scalar> & rowVar() const {return gsNullExpr<T>::get();}
3286  const gsFeSpace<Scalar> & colVar() const {return gsNullExpr<T>::get();}
3287 
3288  // todo mat_expr ?
3289  // tr() const --> _G.data().fundForm(k)
3290 
3291  void print(std::ostream &os) const { os << "jacInv("; _G.print(os); os <<")"; }
3292 };
3293 
3294 /*
3295  Expression for the Jacobian matrix of a FE variable
3296 */
3297 template<class E>
3298 class jac_expr : public _expr<jac_expr<E> >
3299 {
3300  typename E::Nested_t _u;
3301 public:
3302  enum {ColBlocks = (1==E::Space?1:0) };
3303  enum {Space = E::Space, ScalarValued= 0 };
3304 
3305  typedef typename E::Scalar Scalar;
3306 
3307  mutable gsMatrix<Scalar> res;
3308 
3309  jac_expr(const E & _u_)
3310  : _u(_u_) { }
3311 
3312  MatExprType eval(const index_t k) const
3313  {
3314  if (0!=Space)
3315  {
3316  // Dim x (numActive*Dim)
3317  res = _u.data().values[1].col(k).transpose().blockDiag(_u.dim());
3318  }
3319  else
3320  {
3321  res = _u.data().values[1]
3322  .reshapeCol(k, _u.parDim(), _u.targetDim()).transpose()
3323  .blockDiag(_u.dim());
3324  }
3325  return res;
3326  }
3327 
3328  const gsFeSpace<Scalar> & rowVar() const { return _u.rowVar(); }
3329  //const gsFeSpace<Scalar> & rowVar() const { return rowVar_impl<E>(); }
3330  const gsFeSpace<Scalar> & colVar() const { return gsNullExpr<Scalar>::get(); }
3331 
3332  index_t rows() const { return rows_impl(_u); }
3333  index_t cols() const { return cols_impl(_u); }
3334 
3335  // index_t rows() const { return _u.dim(); }
3336  // index_t cols() const { return _u.source().domainDim(); }
3337 
3338  index_t cardinality_impl() const
3339  {
3340  return _u.dim() * _u.data().actives.rows();
3341  }
3342 
3343  void parse(gsExprHelper<Scalar> & evList) const
3344  {
3345  evList.add(_u);
3346  _u.data().flags |= NEED_DERIV|NEED_ACTIVE;
3347  //note: cardinality() depends on actives
3348  }
3349 
3350  void print(std::ostream &os) const { os << "\u2207("; _u.print(os);os <<")"; }
3351 
3352 private:
3353 
3354  // The jacobian is different for gsFeVariable, gsFeSolution and gsFeSpace
3355  // gsFeSolution: Does not work
3356  // gsFeVariable: dim()=1 and source().targetDim()=d
3357  // gsFeSpace: dim()=d and source().targetDim()=1
3358  template<class U> inline
3359  typename util::enable_if<!(util::is_same<U,gsFeSpace<Scalar> >::value), index_t >::type // What about solution??
3360  rows_impl(const U & u) const
3361  {
3362  return u.source().targetDim();
3363  }
3364 
3365  template<class U> inline
3366  typename util::enable_if< (util::is_same<U,gsFeSpace<Scalar> >::value), index_t >::type
3367  rows_impl(const U & u) const
3368  {
3369  return u.dim();
3370  }
3371 
3372  template<class U> inline
3373  typename util::enable_if<!(util::is_same<U,gsFeSpace<Scalar> >::value), index_t >::type
3374  cols_impl(const U & u) const
3375  {
3376  return u.source().domainDim();
3377  }
3378 
3379  template<class U> inline
3380  typename util::enable_if< (util::is_same<U,gsFeSpace<Scalar> >::value), index_t >::type
3381  cols_impl(const U & u) const
3382  {
3383  return u.source().domainDim();
3384  }
3385 
3386  // The jacobian is different for gsFeVariable, gsFeSolution and gsFeSpace
3387  // gsFeSolution: Does not work
3388  // gsFeVariable: rowVar = NULL
3389  // gsFeSpace: rowVar = u
3390  template<class U> inline
3391  typename util::enable_if<!(util::is_same<U,gsFeSpace<Scalar> >::value), const gsFeSpace<Scalar> & >::type
3392  rowVar_impl() const
3393  {
3394  return gsNullExpr<Scalar>::get();
3395  }
3396 
3397  template<class U> inline
3398  typename util::enable_if<(util::is_same<U,gsFeSpace<Scalar> >::value), const gsFeSpace<Scalar> & >::type
3399  rowVar_impl() const
3400  {
3401  return _u;
3402  }
3403 };
3404 
3405 /*
3406  Expression for the Jacobian matrix of a geometry map
3407 */
3408 template<class T>
3409 class jac_expr<gsGeometryMap<T> > : public _expr<jac_expr<gsGeometryMap<T> > >
3410 {
3411  typename gsGeometryMap<T>::Nested_t _G;
3412 
3413 public:
3414  typedef T Scalar;
3415  enum {Space = 0, ScalarValued= 0, ColBlocks= 0};
3416 
3417  jac_expr(const gsGeometryMap<T> & G) : _G(G) { }
3418  MatExprType eval(const index_t k) const
3419  {
3420  // TarDim x ParDim
3421  return _G.data().values[1]
3422  .reshapeCol(k, _G.data().dim.first, _G.data().dim.second).transpose();
3423  }
3424 
3425  index_t rows() const { return _G.source().targetDim(); }
3426 
3427  index_t cols() const { return _G.source().domainDim(); }
3428 
3429  static const gsFeSpace<Scalar> & rowVar() { return gsNullExpr<Scalar>::get(); }
3430  static const gsFeSpace<Scalar> & colVar() { return gsNullExpr<Scalar>::get(); }
3431 
3432  void parse(gsExprHelper<Scalar> & evList) const
3433  {
3434  evList.add(_G);
3435  _G.data().flags |= NEED_DERIV;
3436  }
3437 
3438  meas_expr<T> absDet() const
3439  {
3440  GISMO_ASSERT(rows() == cols(), "The Jacobian matrix is not square");
3441  return meas_expr<T>(_G);
3442  }
3443 
3444  jacInv_expr<T> inv() const
3445  {
3446  GISMO_ASSERT(rows() == cols(), "The Jacobian matrix is not square");
3447  return jacInv_expr<T>(_G);
3448  }
3449 
3451  jacInv_expr<T> ginv() const { return jacInv_expr<T>(_G); }
3452 
3453  void print(std::ostream &os) const { os << "\u2207("; _G.print(os); os <<")"; }
3454 };
3455 
3456 template<class E>
3457 class hess_expr : public _expr<hess_expr<E> >
3458 {
3459 public:
3460  typedef typename E::Scalar Scalar;
3461 private:
3462  typename E::Nested_t _u;
3463  mutable gsMatrix<Scalar> res;
3464 public:
3465  enum {ScalarValued = 0, ColBlocks = (1==E::Space?1:0) };
3466  enum {Space = E::Space };
3467 
3468 public:
3469  hess_expr(const E & u) : _u(u)
3470  {
3471  //gsInfo << "\n-expression is space ? "<<E::Space <<"\n"; _u.print(gsInfo);
3472  //GISMO_ASSERT(1==_u.dim(),"hess(.) requires 1D variable");
3473  }
3474 
3475  const gsMatrix<Scalar> & eval(const index_t k) const
3476  {
3477  const gsFuncData<Scalar> & dd = _u.data();
3478  const index_t sz = cardinality_impl();
3479  res.resize(dd.dim.first, sz*dd.dim.first);
3480  secDerToHessian(dd.values[2].col(k), dd.dim.first, res);
3481  res.resize(dd.dim.first, res.cols()*dd.dim.first);
3482  // Note: auto returns by value here,
3483  // in C++11 we may add in -> decltype(res) &
3484  return res;
3485  }
3486 
3487  index_t rows() const { return _u.data().dim.first; }
3488  index_t cols() const
3489  { return rows();
3490  //return 2*_u.data().values[2].rows() / (1+_u.data().dim.first);
3491  }
3492 
3493  index_t cardinality_impl() const
3494  {
3495  return 2*_u.data().values[2].rows()/
3496  (_u.data().dim.first*(1+_u.data().dim.first));
3497  //gsDebugVar(_u.data().values.front().rows());//empty!
3498  }
3499 
3500  void parse(gsExprHelper<Scalar> & evList) const
3501  {
3502  evList.add(_u);
3503  _u.data().flags |= NEED_2ND_DER;
3504  }
3505 
3506  const gsFeSpace<Scalar> & rowVar() const { return _u.rowVar(); }
3507  const gsFeSpace<Scalar> & colVar() const { return gsNullExpr<Scalar>::get(); }
3508 
3509  void print(std::ostream &os) const
3510  // { os << "hess("; _u.print(os);os <<")"; }
3511  { os << "\u210D(U)"; }
3512 };
3513 
3514 template<class T>
3515 class hess_expr<gsFeSolution<T> > : public _expr<hess_expr<gsFeSolution<T> > >
3516 {
3517 protected:
3518  const gsFeSolution<T> _u;
3519 
3520 public:
3521  typedef T Scalar;
3522  enum{Space = 0, ScalarValued = 0, ColBlocks = 0 };
3523 
3524  hess_expr(const gsFeSolution<T> & u) : _u(u) { }
3525 
3526  mutable gsMatrix<T> res;
3527  const gsMatrix<T> & eval(const index_t k) const
3528  {
3529  GISMO_ASSERT(1==_u.data().actives.cols(), "Single actives expected. Actives: \n"<<_u.data().actives);
3530 
3531  const gsDofMapper & map = _u.mapper();
3532 
3533  const index_t numActs = _u.data().values[0].rows();
3534  const index_t pdim = _u.parDim();
3535  index_t numDers = pdim*(pdim+1)/2;
3536  gsMatrix<T> deriv2;
3537 
3538  // In the scalar case, the hessian is returned as a pdim x pdim matrix
3539  if (1==_u.dim())
3540  {
3541  res.setZero(numDers,1);
3542  for (index_t i = 0; i!=numActs; ++i)
3543  {
3544  const index_t ii = map.index(_u.data().actives.at(i), _u.data().patchId,0);
3545  deriv2 = _u.data().values[2].block(i*numDers,k,numDers,1);
3546  if ( map.is_free_index(ii) ) // DoF value is in the solVector
3547  res += _u.coefs().at(ii) * deriv2;
3548  else
3549  res +=_u.fixedPart().at( map.global_to_bindex(ii) ) * deriv2;
3550  }
3551  secDerToHessian(res, pdim, deriv2);
3552  res.swap(deriv2);
3553  res.resize(pdim,pdim);
3554  }
3555  // In the vector case, the hessian is returned as a matrix where each row corresponds to the component of the solution and contains the derivatives in the columns
3556  else
3557  {
3558  res.setZero(rows(), numDers);
3559  for (index_t c = 0; c != _u.dim(); c++)
3560  for (index_t i = 0; i != numActs; ++i) {
3561  const index_t ii = map.index(_u.data().actives.at(i), _u.data().patchId, c);
3562  deriv2 = _u.space().data().values[2].block(i * numDers, k, numDers,
3563  1).transpose(); // start row, start col, rows, cols
3564  if (map.is_free_index(ii)) // DoF value is in the solVector
3565  res.row(c) += _u.coefs().at(ii) * deriv2;
3566  else
3567  res.row(c) += _u.fixedPart().at(map.global_to_bindex(ii)) * deriv2;
3568  }
3569  }
3570  return res;
3571  }
3572 
3573  index_t rows() const
3574  {
3575  if (1==_u.dim())
3576  return _u.parDim();
3577  else
3578  return _u.dim(); // number of components
3579  }
3580  index_t cols() const
3581  {
3582  if (1==_u.dim())
3583  return _u.parDim();
3584  // second derivatives in the columns; i.e. [d11, d22, d33, d12, d13, d23]
3585  else
3586  return _u.parDim() * (_u.parDim() + 1) / 2;
3587  }
3588 
3589  const gsFeSpace<Scalar> & rowVar() const { return gsNullExpr<Scalar>::get(); }
3590  const gsFeSpace<Scalar> & colVar() const { return gsNullExpr<Scalar>::get(); }
3591 
3592  void parse(gsExprHelper<Scalar> & evList) const
3593  {
3594  _u.parse(evList); // add symbol
3595  evList.add(_u.space());
3596  _u.data().flags |= NEED_ACTIVE | NEED_VALUE | NEED_DERIV2;
3597  }
3598 
3599  void print(std::ostream &os) const { os << "\u210D(s)"; }
3600 };
3601 
3602 
3603 /*
3604  Expression for the partial derivative (matrices) of the Jacobian
3605  matrix of the geometry map
3606 */
3607 template<class T>
3608 class dJacG_expr : public _expr<dJacG_expr<T> >
3609 {
3610  typename gsGeometryMap<T>::Nested_t _G;
3611 
3612  mutable gsMatrix<T> res;
3613 public:
3614  typedef T Scalar;
3615 
3616  dJacG_expr(const gsGeometryMap<T> & G) : _G(G) { }
3617 
3618  MatExprType eval(const index_t k) const
3619  {
3620  const index_t sz = _G.data().values[0].rows();
3621  const index_t s = _G.data().derivSize(); //dim.first*(_G.data().dim.first+1)/2;
3622  (void)s;
3623  res.resize(_G.data().dim.second, sz*_G.data().dim.first);
3624  res.setOnes();//_G.data().values[2].segment(i*k,k); // todo
3625  return res;
3626  }
3627 
3628  index_t rows() const { return _G.source().targetDim(); }
3629  index_t cols() const { return _G.source().domainDim(); }
3630 
3631  void parse(gsExprHelper<Scalar> & evList) const
3632  {
3633  evList.add(_G);
3634  _G.data().flags |= NEED_2ND_DER;
3635  }
3636 };
3637 
3638 
3639 /*
3640  Expression for the curl
3641 */
3642 template<class T>
3643 class curl_expr : public _expr<curl_expr<T> >
3644 {
3645 public:
3646  typedef T Scalar;
3647 private:
3648  typename gsFeVariable<T>::Nested_t _u;
3649  mutable gsMatrix<Scalar> res;
3650 public:
3651  enum{ Space = 1, ScalarValued= 0, ColBlocks= 0};
3652 
3653  curl_expr(const gsFeVariable<T> & u) : _u(u)
3654  { GISMO_ASSERT(3==u.dim(),"curl(.) requires 3D variable."); }
3655 
3656  const gsMatrix<T> & eval(const index_t k) const
3657  {
3658  res.setZero( rows(), _u.dim());
3659  const index_t na = _u.data().values[0].rows();
3660  gsAsConstMatrix<T, Dynamic, Dynamic> pd =
3661  _u.data().values[1].reshapeCol(k, cols(), na);
3662 
3663  res.col(0).segment(na ,na) = -pd.row(2);
3664  res.col(0).segment(2*na,na) = pd.row(1);
3665  res.col(1).segment(0 ,na) = pd.row(2);
3666  res.col(1).segment(2*na,na) = -pd.row(0);
3667  res.col(2).segment(0 ,na) = -pd.row(1);
3668  res.col(2).segment(na ,na) = pd.row(0);
3669  return res;
3670  }
3671 
3672  index_t rows() const { return _u.dim() * _u.data().values[0].rows(); }
3673  index_t cols() const { return _u.data().dim.first; }
3674 
3675  void parse(gsExprHelper<Scalar> & evList) const
3676  {
3677  evList.add(_u);
3678  _u.data().flags |= NEED_GRAD;
3679  }
3680 
3681  const gsFeSpace<T> & rowVar() const { return _u.rowVar(); }
3682  const gsFeSpace<T> & colVar() const {return gsNullExpr<T>::get();}
3683 
3684  void print(std::ostream &os) const { os << "curl("; _u.print(os); os <<")"; }
3685 };
3686 
3687 /*
3688  Expression for multiplication operation (first version)
3689 
3690  First argument E1 has ColBlocks = false
3691 
3692  Partial specialization for (right) blockwise multiplication
3693 
3694  B * [A1 A2 A3] = [B*A1 B*A2 B*A3]
3695 
3696 */
3697 template <typename E1, typename E2>
3698 class mult_expr<E1,E2,false> : public _expr<mult_expr<E1, E2, false> >
3699 {
3700  typename E1::Nested_t _u;
3701  typename E2::Nested_t _v;
3702 
3703 public:
3704  enum {ScalarValued = E1::ScalarValued && E2::ScalarValued,
3705  ColBlocks = E2::ColBlocks};
3706  enum {Space = (int)E1::Space + (int)E2::Space };
3707 
3708  typedef typename E1::Scalar Scalar;
3709 
3710  mult_expr(_expr<E1> const& u,
3711  _expr<E2> const& v)
3712  : _u(u), _v(v) { }
3713 
3714  mutable Temporary_t tmp;
3715  const Temporary_t & eval(const index_t k) const
3716  {
3717  GISMO_ASSERT(0==_u.cols()*_v.rows() || _u.cols() == _v.rows(),
3718  "Wrong dimensions "<<_u.cols()<<"!="<<_v.rows()<<" in * operation:\n"
3719  << _u <<" times \n" << _v );
3720 
3721  // Note: a * b * c --> (a*b).eval()*c
3722  tmp = _u.eval(k) * _v.eval(k);
3723  return tmp; // assumes result is not scalarvalued
3724  }
3725 
3726  index_t rows() const { return E1::ScalarValued ? _v.rows() : _u.rows(); }
3727  index_t cols() const { return E2::ScalarValued ? _u.cols() : _v.cols(); }
3728  void parse(gsExprHelper<Scalar> & evList) const
3729  { _u.parse(evList); _v.parse(evList); }
3730 
3731 
3732  index_t cardinality_impl() const
3733  { return 0==E1::Space ? _v.cardinality(): _u.cardinality(); }
3734 
3735  const gsFeSpace<Scalar> & rowVar() const
3736  { return 0==E1::Space ? _v.rowVar() : _u.rowVar(); }
3737  const gsFeSpace<Scalar> & colVar() const
3738  { return 0==E2::Space ? _u.colVar() : _v.colVar(); }
3739 
3740  void print(std::ostream &os) const { _u.print(os); os<<"*"; _v.print(os); }
3741 };
3742 
3743 /*
3744  Expression for multiplication operation (second version)
3745 
3746  First argument E1 has ColBlocks = true
3747 
3748  Partial specialization for (right) blockwise multiplication
3749  [A1 A2 A3] * B = [A1*B A2*B A3*B]
3750 
3751  as well as
3752 
3753  both are ColBlocks: [A1 A2 A3] * [B1 B2 B3] = [A1*B1 A2*B2 A3*B3]
3754  [A2*B1 .. ]
3755  [ ]
3756 */
3757 template <typename E1, typename E2>
3758 class mult_expr<E1, E2, true> : public _expr<mult_expr<E1, E2, true> >
3759 {
3760 public:
3761  typedef typename E2::Scalar Scalar;
3762 private:
3763  typename E1::Nested_t _u;
3764  typename E2::Nested_t _v;
3765 
3766  mutable gsMatrix<Scalar> res;
3767 public:
3768  enum {ScalarValued = 0, ColBlocks = E1::ColBlocks}; //(!)
3769  enum {Space = (int)E1::Space + (int)E2::Space };
3770 
3771  mult_expr(_expr<E1> const& u,
3772  _expr<E2> const& v)
3773  : _u(u), _v(v)
3774  {
3775 
3776  }
3777 
3778  const gsMatrix<Scalar> & eval(const index_t k) const
3779  {
3780  const index_t uc = _u.cols();
3781  const index_t ur = _u.rows();
3782  const index_t nb = _u.cardinality();
3783  const auto tmpA = _u.eval(k);
3784  const auto tmpB = _v.eval(k);
3785 
3786  const index_t vc = _v.cols();
3787 
3788  // either _v.cardinality()==1 or _v.cardinality()==_u.cardinality()
3789  if ( 1 == _v.cardinality() ) //second is not ColBlocks
3790  {
3791  res.resize(ur, vc*nb);
3792  GISMO_ASSERT(tmpA.cols()==uc*nb, "Dimension error.. "<< tmpA.cols()<<"!="<<uc*nb );
3793  GISMO_ASSERT(1==_v.cardinality(), "Dimension error");
3794  //gsInfo<<"cols = "<<res.cols()<<"; rows = "<<res.rows()<<"\n";
3795  for (index_t i = 0; i!=nb; ++i)
3796  res.middleCols(i*vc,vc).noalias()
3797  = tmpA.middleCols(i*uc,uc) * tmpB;
3798  }
3799  // both are ColBlocks: [A1 A2 A3] * [B1 B2 B3] = [A1*B1 A2*B2 A3*B3]
3800  // [A2*B1 .. ]
3801  // [ ]
3802  else
3803  {
3804  const index_t nbv = _v.cardinality();
3805  res.resize(ur*nb, vc*nbv);
3806  for (index_t i = 0; i!=nb; ++i)
3807  for (index_t j = 0; j!=nbv; ++j)
3808  {
3809  res.block(i*ur,j*vc,ur,vc).noalias() =
3810  tmpA.middleCols(i*uc,uc) * tmpB.middleCols(j*vc,vc);
3811  // res.middleCols(i*vc,vc).noalias()
3812  // = tmpA.middleCols(i*uc,uc) * tmpB.middleCols(i*vc,vc);
3813  }
3814  }
3815  return res;
3816  }
3817 
3818  index_t rows() const {
3819  return _u.rows();
3820  }
3821  index_t cols() const {
3822  return _v.cols();
3823  }
3824 
3825  void parse(gsExprHelper<Scalar> & evList) const
3826  { _u.parse(evList); _v.parse(evList); }
3827 
3828 
3829  index_t cardinality_impl() const { return _u.cardinality(); }
3830 
3831  const gsFeSpace<Scalar> & rowVar() const { return _u.rowVar(); }
3832  const gsFeSpace<Scalar> & colVar() const
3833  {
3834  if ( 1 == _v.cardinality() )
3835  return _u.colVar();
3836  else
3837  return _v.colVar();
3838  }
3839 
3840  void print(std::ostream &os) const
3841  { os << "("; _u.print(os);os <<"*";_v.print(os);os << ")"; }
3842 };
3843 
3844 /*
3845  Expression for multiplication operation (third version)
3846 
3847  Scalar multiplication
3848 */
3849 template <typename E2>
3850 class mult_expr<typename E2::Scalar, E2, false>
3851  : public _expr<mult_expr<typename E2::Scalar, E2, false> >
3852 // template <typename E> class scmult_expr : public _expr<scmult_expr<E> >
3853 {
3854 public:
3855  typedef typename E2::Scalar Scalar;
3856 private:
3857  Scalar const _c;
3858  typename E2::Nested_t _v;
3859 
3860  //mult_expr(const mult_expr&);
3861 public:
3862  enum {ScalarValued = E2::ScalarValued, ColBlocks = E2::ColBlocks};
3863  enum {Space = E2::Space};
3864 
3865  mult_expr(Scalar const & c, _expr<E2> const& v)
3866  : _c(c), _v(v) { }
3867 
3868  EIGEN_STRONG_INLINE AutoReturn_t eval(const index_t k) const
3869  {
3870  return ( _c * _v.eval(k) );
3871 
3872  }
3873 
3874  index_t rows() const { return _v.rows(); }
3875  index_t cols() const { return _v.cols(); }
3876 
3877  void parse(gsExprHelper<Scalar> & evList) const
3878  { _v.parse(evList); }
3879 
3880  index_t cardinality_impl() const
3881  { return _v.cardinality(); }
3882 
3883  const gsFeSpace<Scalar> & rowVar() const { return _v.rowVar(); }
3884  const gsFeSpace<Scalar> & colVar() const { return _v.colVar(); }
3885 
3886  void print(std::ostream &os) const { os << _c <<"*";_v.print(os); }
3887 };
3888 
3889 
3890 template <typename E1, typename E2>
3891 class collapse_expr : public _expr<collapse_expr<E1, E2> >
3892 {
3893  typename E1::Nested_t _u;
3894  typename E2::Nested_t _v;
3895 
3896 public:
3897  enum {ScalarValued = 0, ColBlocks = 0};
3898  enum { Space = (int)E1::Space + (int)E2::Space };
3899 
3900  typedef typename E1::Scalar Scalar;
3901 
3902  mutable gsMatrix<Scalar> res;
3903 
3904  collapse_expr(_expr<E1> const& u,
3905  _expr<E2> const& v)
3906  : _u(u), _v(v) { }
3907 
3908  //EIGEN_STRONG_INLINE MatExprType
3909  const gsMatrix<Scalar> &
3910  eval(const index_t k) const
3911  {
3912  const index_t nb = rows();
3913  const auto tmpA = _u.eval(k);
3914  const auto tmpB = _v.eval(k);
3915 
3916  if (E1::ColBlocks)
3917  {
3918  const index_t ur = _v.rows();
3919  res.resize(nb, ur);
3920  for (index_t i = 0; i!=nb; ++i)
3921  {
3922  res.row(i).transpose().noalias() = tmpA.middleCols(i*ur,ur) * tmpB;
3923  }
3924  }
3925  else if (E2::ColBlocks)
3926  {
3927  const index_t ur = _u.cols();
3928  res.resize(nb, ur);
3929  for (index_t i = 0; i!=nb; ++i)
3930  {
3931  res.row(i).noalias() = tmpA * tmpB.middleCols(i*ur,ur);
3932  }
3933  }
3934 
3935  return res;
3936  }
3937 
3938  index_t rows() const { return E1::ColBlocks ? _u.cols() / _v.rows() : _v.cols() / _u.cols() ; }
3939  index_t cols() const { return E1::ColBlocks ? _v.rows() : _u.cols(); }
3940 
3941  void parse(gsExprHelper<Scalar> & evList) const
3942  { _u.parse(evList); _v.parse(evList); }
3943 
3944  const gsFeSpace<Scalar> & rowVar() const
3945  { return E1::ColBlocks ? _u.rowVar() : _v.rowVar(); }
3946  const gsFeSpace<Scalar> & colVar() const
3947  {
3948  GISMO_ERROR("none");
3949  }
3950 
3951  void print(std::ostream &os) const { _u.print(os); os<<"~"; _v.print(os); }
3952 };
3953 
3954 // Multi-matrix collapsed by a vector
3955 template <typename E1, typename E2> //EIGEN_STRONG_INLINE
3956 //collapse_expr<E1,E2> const operator&(<E1> const& u, _expr<E2> const& v)
3957 collapse_expr<E1,E2> collapse( _expr<E1> const& u, _expr<E2> const& v)
3958 { return collapse_expr<E1, E2>(u, v); }
3959 
3960 
3961 /*
3962  Expression for the Frobenius matrix (or double dot) product (first
3963  version) Also block-wise
3964 
3965  [A1 A2 A3] . [B1 B2 B3]
3966  =
3967  [ A1.B1 A1.B2 A1.B3 ]
3968  [ A2.B1 A2.B2 A2.B3 ]
3969  [ A3.B1 A3.B2 A3.B3 ]
3970 */
3971 template <typename E1, typename E2, bool = E2::ColBlocks>
3972 class frprod_expr : public _expr<frprod_expr<E1, E2> >
3973 {
3974 public:
3975  typedef typename E1::Scalar Scalar;
3976  enum {ScalarValued = 0, ColBlocks=E2::ColBlocks};
3977  enum { Space = (int)E1::Space + (int)E2::Space };
3978  // E1 E2 this (16 cases..)
3979  // 0 0 0
3980  // 1 1
3981  // 2 2
3982  // 3 3
3983 
3984 private:
3985  typename E1::Nested_t _u;
3986  typename E2::Nested_t _v;
3987 
3988  mutable gsMatrix<Scalar> res;
3989 
3990 public:
3991 
3992  frprod_expr(_expr<E1> const& u, _expr<E2> const& v)
3993  : _u(u), _v(v)
3994  {
3995  //todo: add check() functions, which will evaluate expressions on an empty matrix (no points) to setup initial dimensions ???
3996  //GISMO_ASSERT(_u.rows() == _v.rows(),
3997  // "Wrong dimensions "<<_u.rows()<<"!="<<_v.rows()<<" in % operation");
3998  //GISMO_ASSERT(_u.cols() == _v.cols(),
3999  // "Wrong dimensions "<<_u.cols()<<"!="<<_v.cols()<<" in % operation");
4000  }
4001 
4002  const gsMatrix<Scalar> & eval(const index_t k) const //todo: specialize for nb==1
4003  {
4004  // assert _u.size()==_v.size()
4005  const index_t rb = _u.rows();
4006  const index_t nb = _u.cardinality();
4007  auto A = _u.eval(k);
4008  auto B = _v.eval(k);
4009  res.resize(nb, nb);
4010  for (index_t i = 0; i!=nb; ++i) // all with all
4011  for (index_t j = 0; j!=nb; ++j)
4012  res(i,j) =
4013  (A.middleCols(i*rb,rb).array() * B.middleCols(j*rb,rb).array()).sum();
4014  return res;
4015  }
4016 
4017  index_t rows() const { return _u.cols() / _u.rows(); }
4018  index_t cols() const { return _u.cols() / _u.rows(); }
4019 
4020  void parse(gsExprHelper<Scalar> & evList) const
4021  { _u.parse(evList); _v.parse(evList); }
4022 
4023  const gsFeSpace<Scalar> & rowVar() const { return _u.rowVar(); }
4024  const gsFeSpace<Scalar> & colVar() const { return _v.colVar(); }
4025 
4026  void print(std::ostream &os) const
4027  { os << "("; _u.print(os); os<<" % "; _v.print(os); os<<")";}
4028 };
4029 
4030 /*
4031 
4032  Expression for the Frobenius matrix (or double dot) product (second
4033  version), When left hand only side is block-wise
4034 
4035  [A1 A2 A3] : B = [A1:B A2:B A3:B]
4036 */
4037 template <typename E1, typename E2>
4038 class frprod_expr<E1,E2,false> : public _expr<frprod_expr<E1, E2,false> >
4039 {
4040 public:
4041  typedef typename E1::Scalar Scalar;
4042  enum {ScalarValued = 0, Space = E1::Space, ColBlocks= E1::ColBlocks};
4043 
4044 private:
4045  typename E1::Nested_t _u;
4046  typename E2::Nested_t _v;
4047 
4048  mutable gsMatrix<Scalar> res;
4049 
4050 public:
4051 
4052  frprod_expr(_expr<E1> const& u, _expr<E2> const& v)
4053  : _u(u), _v(v)
4054  {
4055  // gsInfo << "expression is space ? "<<E1::Space <<"\n"; _u.print(gsInfo);
4056  // GISMO_ASSERT(_u.rows() == _v.rows(),
4057  // "Wrong dimensions "<<_u.rows()<<"!="<<_v.rows()<<" in % operation");
4058  // GISMO_ASSERT(_u.cols() == _v.cols(),
4059  // "Wrong dimensions "<<_u.cols()<<"!="<<_v.cols()<<" in % operation");
4060  }
4061 
4062  const gsMatrix<Scalar> & eval(const index_t k) const //todo: specialize for nb==1
4063  {
4064  // assert _u.size()==_v.size()
4065  auto A = _u.eval(k);
4066  auto B = _v.eval(k);
4067  const index_t rb = A.rows(); //==cb
4068  const index_t nb = _u.cardinality();
4069  res.resize(nb, 1);
4070  for (index_t i = 0; i!=nb; ++i) // all with all
4071  res(i,0) =
4072  (A.middleCols(i*rb,rb).array() * B.array()).sum();
4073  return res;
4074  }
4075 
4076  index_t rows() const { return _u.cols() / _u.rows(); }
4077  index_t cols() const { return 1; }
4078 
4079  void parse(gsExprHelper<Scalar> & evList) const
4080  { _u.parse(evList); _v.parse(evList); }
4081 
4082 
4083  const gsFeSpace<Scalar> & rowVar() const { return _u.rowVar(); }
4084  const gsFeSpace<Scalar> & colVar() const { return _v.rowVar(); }
4085 
4086  void print(std::ostream &os) const
4087  { os << "("; _u.print(os); os<<" % "; _v.print(os); os<<")";}
4088 };
4089 
4090 /*
4091  Expression for scalar division operation (first version)
4092 */
4093 template <typename E1, typename E2>
4094 class divide_expr : public _expr<divide_expr<E1,E2> >
4095 {
4096  typename E1::Nested_t _u;
4097  typename E2::Nested_t _v;
4098 
4099 public:
4100  typedef typename E1::Scalar Scalar;
4101 
4102 public:
4103  enum {ScalarValued = E1::ScalarValued, ColBlocks= E2::ColBlocks};
4104  enum {Space = E1::Space}; // The denominator E2 has to be scalar.
4105 
4106  divide_expr(_expr<E1> const& u, _expr<E2> const& v)
4107  : _u(u), _v(v)
4108  {
4109  GISMO_STATIC_ASSERT(E2::ScalarValued, "The denominator needs to be scalar valued.");
4110  }
4111 
4112  AutoReturn_t eval(const index_t k) const
4113  { return ( _u.eval(k) / _v.eval(k) ); }
4114 
4115  index_t rows() const { return _u.rows(); }
4116  index_t cols() const { return _u.cols(); }
4117 
4118  void parse(gsExprHelper<Scalar> & evList) const
4119  { _u.parse(evList); _v.parse(evList); }
4120 
4121 
4122  const gsFeSpace<Scalar> & rowVar() const { return _u.rowVar(); }
4123  const gsFeSpace<Scalar> & colVar() const { return _u.colVar(); }
4124 
4125  void print(std::ostream &os) const
4126  { os << "("; _u.print(os);os <<" / ";_v.print(os);os << ")"; }
4127 };
4128 
4129 /*
4130  Division specialization (second version) for constant value denominator
4131 */
4132 template <typename E1>
4133 class divide_expr<E1,typename E1::Scalar>
4134  : public _expr<divide_expr<E1,typename E1::Scalar> >
4135 {
4136 public:
4137  typedef typename E1::Scalar Scalar;
4138 
4139 private:
4140  typename E1::Nested_t _u;
4141  Scalar const _c;
4142 
4143 public:
4144  enum {Space= E1::Space, ScalarValued = E1::ScalarValued, ColBlocks= E1::ColBlocks};
4145 
4146  divide_expr(_expr<E1> const& u, Scalar const c)
4147  : _u(u), _c(c) { }
4148 
4149  AutoReturn_t eval(const index_t k) const
4150  { return ( _u.eval(k) / _c ); }
4151 
4152  index_t rows() const { return _u.rows(); }
4153  index_t cols() const { return _u.cols(); }
4154 
4155  void parse(gsExprHelper<Scalar> & evList) const
4156  { _u.parse(evList); }
4157 
4158 
4159  const gsFeSpace<Scalar> & rowVar() const { return _u.rowVar(); }
4160  const gsFeSpace<Scalar> & colVar() const { return _u.colVar(); }
4161 
4162  void print(std::ostream &os) const
4163  { os << "("; _u.print(os);os <<"/"<< _c << ")"; }
4164 };
4165 
4166 /*
4167  Division specialization (third version) for constant value
4168  numerator
4169 */
4170 template <typename E2>
4171 class divide_expr<typename E2::Scalar,E2>
4172  : public _expr<divide_expr<typename E2::Scalar,E2> >
4173 {
4174 public:
4175  typedef typename E2::Scalar Scalar;
4176 
4177 private:
4178  Scalar const _c;
4179  typename E2::Nested_t _u;
4180 public:
4181  enum {Space= 0, ScalarValued = 1, ColBlocks= 0};
4182 
4183  divide_expr(Scalar const c, _expr<E2> const& u)
4184  : _c(c), _u(u)
4185  { GISMO_STATIC_ASSERT(E2::ScalarValued, "The denominator needs to be scalar valued."); }
4186 
4187  Scalar eval(const index_t k) const
4188  { return ( _c / _u.val().eval(k) ); }
4189 
4190  index_t rows() const { return 0; }
4191  index_t cols() const { return 0; }
4192 
4193  void parse(gsExprHelper<Scalar> & evList) const
4194  { _u.parse(evList); }
4195 
4196 
4197  const gsFeSpace<Scalar> & rowVar() const { return _u.rowVar(); }
4198  const gsFeSpace<Scalar> & colVar() const { return _u.colVar(); }
4199 
4200  void print(std::ostream &os) const
4201  { os << "("<< _c <<"/";_u.print(os);os << ")";}
4202 };
4203 
4204 /*
4205  Expression for addition operation
4206 */
4207 template <typename E1, typename E2>
4208 class add_expr : public _expr<add_expr<E1, E2> >
4209 {
4210  typename E1::Nested_t _u;
4211  typename E2::Nested_t _v;
4212 
4213 public:
4214  enum {ScalarValued = E1::ScalarValued && E2::ScalarValued,
4215  ColBlocks = E1::ColBlocks || E2::ColBlocks };
4216  enum {Space = E1::Space}; // == E2::Space
4217 
4218  typedef typename E1::Scalar Scalar;
4219 
4220  add_expr(_expr<E1> const& u, _expr<E2> const& v)
4221  : _u(u), _v(v)
4222  {
4223  GISMO_ENSURE((int)E1::Space == (int)E2::Space &&
4224  _u.rowVar()==_v.rowVar() && _u.colVar()==_v.colVar(),
4225  "Error: adding apples and oranges (use comma instead),"
4226  " namely:\n" << _u <<"\n"<<_v<<
4227  " \nvars:\n" << _u.rowVar().id()<<"!="<<_v.rowVar().id() <<", "<< _u.colVar().id()<<"!="<<_v.colVar().id()<<
4228  " \nspaces:\n" << (int)E1::Space<< "!="<< (int)E2::Space
4229  );
4230  }
4231 
4232  mutable Temporary_t res;
4233  const Temporary_t & eval(const index_t k) const
4234  {
4235  GISMO_ASSERT(_u.rows() == _v.rows(),
4236  "Wrong dimensions "<<_u.rows()<<"!="<<_v.rows()<<" in + operation:\n"
4237  << _u <<" plus \n" << _v );
4238  GISMO_ASSERT(_u.cols() == _v.cols(),
4239  "Wrong dimensions "<<_u.cols()<<"!="<<_v.cols()<<" in + operation:\n"
4240  << _u <<" plus \n" << _v );
4241  res = _u.eval(k) + _v.eval(k);
4242  return res;
4243  }
4244 
4245  index_t rows() const { return _u.rows(); }
4246  index_t cols() const { return _u.cols(); }
4247 
4248  void parse(gsExprHelper<Scalar> & evList) const
4249  { _u.parse(evList); _v.parse(evList); }
4250 
4251 
4252  index_t cardinality_impl() const { return _u.cardinality_impl(); }
4253 
4254  const gsFeSpace<Scalar> & rowVar() const { return _u.rowVar(); }
4255  const gsFeSpace<Scalar> & colVar() const { return _v.colVar(); }
4256 
4257  void print(std::ostream &os) const
4258  { os << "("; _u.print(os);os <<" + ";_v.print(os);os << ")"; }
4259 };
4260 
4261 /*
4262  lincom_expr (lc) ?
4263  Expression for (square) matrix summation operation
4264 
4265  M [r x r*k] is a list of matrices
4266  Summation is done over k,
4267  [M1 M2 .. Mk]
4268 
4269  u [s x k] is a list of vectors
4270 
4271  Computed quantity is of size [r x r*s] and contains
4272  [ ... sum_k(Mk * u(s,k) ) ... ]_s
4273 */
4274 template <typename E1, typename E2>
4275 class summ_expr : public _expr<summ_expr<E1,E2> >
4276 {
4277 public:
4278  typedef typename E1::Scalar Scalar;
4279 
4280  enum {Space = E1::Space, ScalarValued= 0, ColBlocks= E2::ColBlocks};
4281 
4282  summ_expr(E1 const& u, E2 const& M) : _u(u), _M(M) { }
4283 
4284  const gsMatrix<Scalar> & eval(const index_t k) const
4285  {
4286  auto sl = _u.eval(k);
4287  const index_t sr = sl.rows();
4288  auto ml = _M.eval(k);
4289  const index_t mr = ml.rows();
4290  const index_t mb = _M.cardinality();
4291 
4292  GISMO_ASSERT(_M.cols()==_M.rows(),"Matrix must be square: "<< _M.rows()<<" x "<< _M.cols() << " expr: "<< _M );
4293  GISMO_ASSERT(mb==_u.cols(),"cardinality must match vector, but card(M)="<<_M.cardinality()<<" and cols(u)="<<_u.cols());
4294 
4295  res.setZero(mr, sr * mr);
4296  for (index_t i = 0; i!=sr; ++i)
4297  for (index_t t = 0; t!=mb; ++t) // lc
4298  res.middleCols(i*mr,mr) += sl(i,t) * ml.middleCols(t*mr,mr);
4299  return res;
4300  }
4301 
4302  index_t rows() const { return _M.rows(); }
4303  index_t cols() const { return _M.rows(); } //_u.rows()
4304 
4305  void parse(gsExprHelper<Scalar> & evList) const
4306  { _u.parse(evList); _M.parse(evList); }
4307 
4308  const gsFeSpace<Scalar> & rowVar() const { return _u.rowVar(); }
4309  const gsFeSpace<Scalar> & colVar() const { return gsNullExpr<Scalar>::get(); }
4310 
4311  index_t cardinality_impl() const
4312  { GISMO_ERROR("Something went terribly wrong"); }
4313 
4314  void print(std::ostream &os) const
4315  { os << "sum("; _M.print(os); os<<","; _u.print(os); os<<")"; }
4316 
4317 private:
4318  typename E1::Nested_t _u;
4319  typename E2::Nested_t _M;
4320 
4321  mutable gsMatrix<Scalar> res;
4322 };
4323 
4324 
4325 /*
4326  Expression for subtraction operation
4327 */
4328 template <typename E1, typename E2>
4329 class sub_expr : public _expr<sub_expr<E1, E2> >
4330 {
4331  typename E1::Nested_t _u;
4332  typename E2::Nested_t _v;
4333 
4334 public:
4335  enum {ScalarValued = E1::ScalarValued && E2::ScalarValued,
4336  ColBlocks = E1::ColBlocks || E2::ColBlocks };
4337  enum {Space = E1::Space}; // == E2::Space
4338 
4339  typedef typename E1::Scalar Scalar;
4340 
4341  sub_expr(_expr<E1> const& u, _expr<E2> const& v)
4342  : _u(u), _v(v)
4343  {
4344  GISMO_ENSURE((int)E1::Space == (int)E2::Space &&
4345  _u.rowVar()==_v.rowVar() && _u.colVar()==_v.colVar(),
4346  "Error: substracting apples from oranges (use comma instead),"
4347  " namely:\n" << _u <<"\n"<<_v);
4348  }
4349 
4350  mutable Temporary_t res;
4351  const Temporary_t & eval(const index_t k) const
4352  {
4353  // GISMO_ASSERT(_u.rowVar().id()==_v.rowVar().id() && _u.rowVar().isAcross()==_v.rowVar().isAcross(),
4354  // "The row spaces are not split compatibly.");
4355  // GISMO_ASSERT(_u.colVar().id()==_v.colVar().id() && _u.colVar().isAcross()==_v.colVar().isAcross(),
4356  // "The col spaces are not split compatibly.");
4357  GISMO_ASSERT(_u.rows() == _v.rows(),
4358  "Wrong dimensions "<<_u.rows()<<"!="<<_v.rows()<<" in - operation:\n" << _u <<" minus \n" << _v );
4359  GISMO_ASSERT(_u.cols() == _v.cols(),
4360  "Wrong dimensions "<<_u.cols()<<"!="<<_v.cols()<<" in - operation:\n" << _u <<" minus \n" << _v );
4361  GISMO_ASSERT(_u.cardinality() == _u.cardinality(),
4362  "Cardinality "<< _u.cardinality()<<" != "<< _v.cardinality());
4363  //return (_u.eval(k) - _v.eval(k) ).eval();
4364  //return (_u.eval(k) - _v.eval(k) ); // any temporary matrices eval(.) will leak mem.
4365  res = _u.eval(k) - _v.eval(k);
4366  return res;
4367  }
4368 
4369  index_t rows() const { return _u.rows(); }
4370  index_t cols() const { return _u.cols(); }
4371 
4372  void parse(gsExprHelper<Scalar> & evList) const
4373  { _u.parse(evList); _v.parse(evList); }
4374 
4375  const gsFeSpace<Scalar> & rowVar() const { return _u.rowVar(); }
4376  const gsFeSpace<Scalar> & colVar() const { return _v.colVar(); }
4377 
4378  index_t cardinality_impl() const
4379  {
4380  GISMO_ASSERT(_u.cardinality() == _u.cardinality(),
4381  "Cardinality "<< _u.cardinality()<<" != "<< _v.cardinality());
4382  return _u.cardinality();
4383  }
4384 
4385  void print(std::ostream &os) const
4386  { os << "("; _u.print(os); os<<" - ";_v.print(os); os << ")";}
4387 };
4388 
4389 /*
4390  Expression for symmetrization operation
4391 */
4392 template <typename E>
4393 class symm_expr : public _expr<symm_expr<E> >
4394 {
4395  typename E::Nested_t _u;
4396 
4397  mutable gsMatrix<typename E::Scalar> tmp;
4398 public:
4399  typedef typename E::Scalar Scalar;
4400 
4401  enum { Space = (0==E::Space ? 0 : E::Space), ScalarValued= E::ScalarValued, ColBlocks= E::ColBlocks };
4402 
4403  symm_expr(_expr<E> const& u)
4404  : _u(u) { }
4405 
4406  MatExprType eval(const index_t k) const
4407  {
4408  //const MatExprType tmp = _u.eval(k);
4409  tmp = _u.eval(k);
4410  // todo: avoid temporary or not ?
4411  return tmp * tmp.transpose();
4412  }
4413 
4414  index_t rows() const { return _u.rows(); }
4415  index_t cols() const { return _u.rows(); }
4416 
4417  void parse(gsExprHelper<Scalar> & evList) const
4418  { _u.parse(evList); }
4419 
4420  const gsFeSpace<Scalar> & rowVar() const { return _u.rowVar(); }
4421  const gsFeSpace<Scalar> & colVar() const { return _u.rowVar(); }
4422 
4423  void print(std::ostream &os) const { os << "symm("; _u.print(os); os <<")"; }
4424 };
4425 
4426 template <typename E>
4427 class symmetrize_expr : public _expr<symmetrize_expr<E> >
4428 {
4429  typename E::Nested_t _u;
4430 
4431  mutable gsMatrix<typename E::Scalar> tmp;
4432 public:
4433  enum { Space = (0==E::Space ? 0 : 3), ScalarValued=E::ScalarValued, ColBlocks= E::ColBlocks };
4434  typedef typename E::Scalar Scalar;
4435 
4436  symmetrize_expr(_expr<E> const& u)
4437  : _u(u) { }
4438 
4439  MatExprType eval(const index_t k) const
4440  {
4441  //const MatExprType tmp = _u.eval(k);
4442  tmp = _u.eval(k);
4443  // todo: avoid temporary or not ?
4444  return tmp + tmp.transpose();
4445  }
4446 
4447  index_t rows() const { return _u.rows(); }
4448  index_t cols() const { return _u.rows(); }
4449 
4450  void parse(gsExprHelper<Scalar> & evList) const
4451  { _u.parse(evList); }
4452 
4453  const gsFeSpace<Scalar> & rowVar() const { return _u.rowVar(); }
4454  const gsFeSpace<Scalar> & colVar() const { return _u.rowVar(); }
4455  index_t cardinality_impl() const { return _u.cardinality_impl(); }
4456 
4457  void print(std::ostream &os) const { os << "symmetrize("; _u.print(os); os <<")"; }
4458 };
4459 
4460 /* Symmetrization operation
4461  template <typename E> symm_expr<E> const
4462  symm(_expr<E> const& u) { return symm_expr<E>(u);}
4463 */
4464 
4465 #undef MatExprType
4466 #undef AutoReturn_t
4467 //----------------------------------------------------------------------------------
4468 
4470 EIGEN_STRONG_INLINE idMat_expr id(const index_t dim) { return idMat_expr(dim); }
4471 
4472 EIGEN_STRONG_INLINE constMat_expr ones(const index_t dim) {
4473  gsMatrix<real_t> ones(dim, dim);
4474  ones.fill(1);
4475  return constMat_expr(ones);
4476  }
4477 
4478 EIGEN_STRONG_INLINE constMat_expr mat(const gsMatrix<real_t> mat) { return constMat_expr(mat); }
4479 
4480 // Returns the unit as an expression
4481 //EIGEN_STRONG_INLINE _expr<real_t> one() { return _expr<real_t,true>(1); }
4482 
4483 
4485 template<class E> EIGEN_STRONG_INLINE
4486 abs_expr<E> abs(const E & u) { return abs_expr<E>(u); }
4487 
4489 template<class E> EIGEN_STRONG_INLINE
4490 grad_expr<E> grad(const E & u) { return grad_expr<E>(u); }
4491 
4493 template<class E> EIGEN_STRONG_INLINE
4494 dJacdc_expr<E> dJacdc(const E & u, index_t c) { return dJacdc_expr<E>(u,c); }
4495 
4497 template<class T> EIGEN_STRONG_INLINE
4498 curl_expr<T> curl(const gsFeVariable<T> & u) { return curl_expr<T>(u); }
4499 
4501 template<class T> EIGEN_STRONG_INLINE
4502 nabla_expr<T> nabla(const gsFeVariable<T> & u) { return nabla_expr<T>(u); }
4503 
4505 template<class T> EIGEN_STRONG_INLINE
4506 onormal_expr<T> nv(const gsGeometryMap<T> & u) { return onormal_expr<T>(u); }
4507 
4508 template<class T> EIGEN_STRONG_INLINE
4509 normal_expr<T> sn(const gsGeometryMap<T> & u) { return normal_expr<T>(u); }
4510 
4512 template<class T> EIGEN_STRONG_INLINE
4513 tangent_expr<T> tv(const gsGeometryMap<T> & u) { return tangent_expr<T>(u); }
4514 
4515 template<class E> EIGEN_STRONG_INLINE
4516 lapl_expr<E> lapl(const symbol_expr<E> & u) { return lapl_expr<E>(u); }
4517 
4518 template<class T> EIGEN_STRONG_INLINE
4519 lapl_expr<gsFeSolution<T> > lapl(const gsFeSolution<T> & u)
4520 { return lapl_expr<gsFeSolution<T> >(u); }
4521 
4523 template<class T> EIGEN_STRONG_INLINE fform2nd_expr<T> fform2nd(const gsGeometryMap<T> & G)
4524 { return fform2nd_expr<T>(G); }
4525 
4527 template<class E> EIGEN_STRONG_INLINE
4528 jac_expr<E> jac(const symbol_expr<E> & u) { return jac_expr<E>(u); }
4529 
4531 template<class T> EIGEN_STRONG_INLINE
4532 jac_expr<gsGeometryMap<T> > jac(const gsGeometryMap<T> & G) {return jac_expr<gsGeometryMap<T> >(G);}
4533 
4535 template<class T> EIGEN_STRONG_INLINE
4536 grad_expr<gsFeSolution<T> > jac(const gsFeSolution<T> & s) {return grad_expr<gsFeSolution<T> >(s);}
4537 
4538 template<class E> EIGEN_STRONG_INLINE
4539 hess_expr<E> hess(const symbol_expr<E> & u) { return hess_expr<E>(u); }
4540 
4542 template<class T> EIGEN_STRONG_INLINE
4543 hess_expr<gsGeometryMap<T> > hess(const gsGeometryMap<T> & u) { return hess_expr<gsGeometryMap<T> >(u); }
4544 
4546 template<class T> EIGEN_STRONG_INLINE
4547 hess_expr<gsFeSolution<T> > hess(const gsFeSolution<T> & u) { return hess_expr<gsFeSolution<T> >(u); }
4548 
4550 template<class T> EIGEN_STRONG_INLINE
4551 dJacG_expr<T> dJac(const gsGeometryMap<T> & G) { return dJacG_expr<T>(G); }
4552 
4554 template<class T> EIGEN_STRONG_INLINE
4555 meas_expr<T> meas(const gsGeometryMap<T> & G) { return meas_expr<T>(G); }
4556 
4558 template <typename E1, typename E2> EIGEN_STRONG_INLINE
4559 mult_expr<E1,E2> const operator*(_expr<E1> const& u, _expr<E2> const& v)
4560 { return mult_expr<E1, E2>(u, v); }
4561 
4562 template <typename E2> EIGEN_STRONG_INLINE
4563 mult_expr<typename E2::Scalar,E2,false> const
4564 operator*(typename E2::Scalar const& u, _expr<E2> const& v)
4565 { return mult_expr<typename E2::Scalar, E2, false>(u, v); }
4566 
4567 template <typename E1> EIGEN_STRONG_INLINE
4568 mult_expr<typename E1::Scalar,E1,false> const
4569 operator*(_expr<E1> const& v, typename E1::Scalar const& u)
4570 { return mult_expr<typename E1::Scalar,E1, false>(u, v); }
4571 
4572 template <typename E1> EIGEN_STRONG_INLINE
4573 mult_expr<typename E1::Scalar,E1,false> const
4574 operator-(_expr<E1> const& u)
4575 { return mult_expr<typename E1::Scalar,E1, false>(-1, u); }
4576 
4577 template <typename E> mult_expr<constMat_expr, E> const
4578 operator*( gsMatrix<typename E::Scalar> const& u, _expr<E> const& v)
4579 { return mult_expr<constMat_expr, E>(mat(u), v); }
4580 
4581 template <typename E> mult_expr<E, constMat_expr> const
4582 operator*(_expr<E> const& u, gsMatrix<typename E::Scalar> const& v)
4583 { return mult_expr<E, constMat_expr>(u, mat(v) ); }
4584 
4586 template <typename E1, typename E2> EIGEN_STRONG_INLINE
4587 frprod_expr<E1,E2> const operator%(_expr<E1> const& u, _expr<E2> const& v)
4588 { return frprod_expr<E1, E2>(u, v); }
4589 
4591 template <typename E1, typename E2> EIGEN_STRONG_INLINE
4592 divide_expr<E1,E2> const operator/(_expr<E1> const& u, _expr<E2> const& v)
4593 { return divide_expr<E1,E2>(u, v); }
4594 
4595 template <typename E> EIGEN_STRONG_INLINE
4596 divide_expr<E,typename E::Scalar> const
4597 operator/(_expr<E> const& u, const typename E::Scalar v)
4598 { return divide_expr<E,typename E::Scalar>(u, v); }
4599 
4600 template <typename E> EIGEN_STRONG_INLINE
4601 divide_expr<typename E::Scalar,E> const
4602 operator/(const typename E::Scalar u, _expr<E> const& v)
4603 { return divide_expr<typename E::Scalar,E>(u, v); }
4604 
4606 template <typename E1, typename E2> EIGEN_STRONG_INLINE
4607 add_expr<E1,E2> const operator+(_expr<E1> const& u, _expr<E2> const& v)
4608 { return add_expr<E1, E2>(u, v); }
4609 
4611 template <typename E> EIGEN_STRONG_INLINE
4612 add_expr< E, _expr<typename E::Scalar, true> >
4613 operator+(_expr<E> const& u, const typename E::Scalar v)
4614 { return add_expr<E,_expr<typename E::Scalar>>(u, _expr<typename E::Scalar,true>(v)); }
4615 
4617 template <typename E> EIGEN_STRONG_INLINE
4618 add_expr< E, _expr<typename E::Scalar, true> >
4619 operator+(const typename E::Scalar v, _expr<E> const& u)
4620 { return add_expr<E,_expr<typename E::Scalar>>(u, _expr<typename E::Scalar,true>(v)); }
4621 
4623 template <typename E1, typename E2> EIGEN_STRONG_INLINE
4624 summ_expr<E1,E2> const summ(E1 const & u, E2 const& M)
4625 { return summ_expr<E1,E2>(u, M); }
4626 
4629 template <typename E1, typename E2> EIGEN_STRONG_INLINE
4630 matrix_by_space_expr<E1,E2> const matrix_by_space(E1 const & u, E2 const& v)
4631 { return matrix_by_space_expr<E1,E2>(u, v); }
4632 
4635 template <typename E1, typename E2> EIGEN_STRONG_INLINE
4636 matrix_by_space_expr_tr<E1,E2> const matrix_by_space_tr(E1 const & u, E2 const& v)
4637 { return matrix_by_space_expr_tr<E1,E2>(u, v); }
4638 
4640 template <typename E1, typename E2> EIGEN_STRONG_INLINE
4641 sub_expr<E1,E2> const operator-(_expr<E1> const& u, _expr<E2> const& v)
4642 { return sub_expr<E1, E2>(u, v); }
4643 
4644 template <typename E2> EIGEN_STRONG_INLINE
4645 sub_expr<_expr<typename E2::Scalar>,E2> const
4646 operator-(typename E2::Scalar const& s, _expr<E2> const& v)
4647 {
4648  // assert E2::ScalarValued
4649  return sub_expr<_expr<typename E2::Scalar>, E2>(_expr<typename E2::Scalar>(s), v);
4650 }
4651 
4652 
4653 // Shortcuts for common quantities, for instance function
4654 // transformations by the geometry map \a G
4655 #define GISMO_SHORTCUT_VAR_EXPRESSION(name,impl) template<class E> EIGEN_STRONG_INLINE \
4656  auto name(const E & u) -> decltype(impl) { return impl; }
4657 #define GISMO_SHORTCUT_MAP_EXPRESSION(name,impl) template<class T> EIGEN_STRONG_INLINE \
4658  auto name(const gsGeometryMap<T> & G) -> decltype(impl) { return impl; }
4659 #define GISMO_SHORTCUT_PHY_EXPRESSION(name,impl) template<class E> EIGEN_STRONG_INLINE \
4660  auto name(const E & u, const gsGeometryMap<typename E::Scalar> & G) -> decltype(impl) { return impl; }
4661 
4662 // Divergence
4663 GISMO_SHORTCUT_VAR_EXPRESSION( div, jac(u).trace() )
4664 GISMO_SHORTCUT_PHY_EXPRESSION( idiv, ijac(u,G).trace() )
4665 
4666 // The unit (normalized) boundary (outer pointing) normal
4667 GISMO_SHORTCUT_MAP_EXPRESSION(unv, nv(G).normalized() )
4668 // The unit (normalized) boundary (surface) normal
4669 GISMO_SHORTCUT_MAP_EXPRESSION(usn, sn(G).normalized() )
4670 
4671 GISMO_SHORTCUT_PHY_EXPRESSION(igrad, grad(u)*jac(G).ginv() ) // transpose() problem ??
4672 GISMO_SHORTCUT_VAR_EXPRESSION(igrad, grad(u) ) // u is presumed to be defined over G
4673 
4674 GISMO_SHORTCUT_PHY_EXPRESSION( ijac, jac(u) * jac(G).ginv())
4675 
4676 // note and todo: does this work for non-scalar solutions?
4677 GISMO_SHORTCUT_PHY_EXPRESSION(ihess,
4678  jac(G).ginv().tr()*( hess(u) - summ(igrad(u,G),hess(G)) ) * jac(G).ginv() )
4679 GISMO_SHORTCUT_VAR_EXPRESSION(ihess, hess(u) )
4680 
4681 GISMO_SHORTCUT_PHY_EXPRESSION(ilapl, ihess(u,G).trace() )
4682 GISMO_SHORTCUT_VAR_EXPRESSION(ilapl, hess(u).trace() )
4683 
4684 GISMO_SHORTCUT_VAR_EXPRESSION(fform, jac(u).tr()*jac(u) )
4685 GISMO_SHORTCUT_VAR_EXPRESSION(shapeop, fform(u).inv() * fform2nd(u) )
4686 
4687 #undef GISMO_SHORTCUT_PHY_EXPRESSION
4688 #undef GISMO_SHORTCUT_VAR_EXPRESSION
4689 #undef GISMO_SHORTCUT_MAP_EXPRESSION
4690 
4691 } // namespace expr
4692 
4693 } //namespace gismo
EIGEN_STRONG_INLINE matrix_by_space_expr< E1, E2 > const matrix_by_space(E1 const &u, E2 const &v)
Definition: gsExpressions.h:4630
asdiag_expr< E > asDiag() const
Returns a diagonal matrix expression of the vector expression.
Definition: gsExpressions.h:309
#define MatExprType
[Include namespace]
Definition: gsThinShellUtils.h:32
mult_expr< real_t, ppart_expr< mult_expr< double, E, false > >, false > npart() const
Returns the expression&#39;s negative part.
Definition: gsExpressions.h:260
Definition: gsExpressions.h:120
Definition: gsExpressions.h:96
Definition: gsExpressions.h:137
size_t mapSize() const
Returns the total number of patch-local degrees of freedom that are being mapped. ...
Definition: gsDofMapper.h:473
Gradient of the object.
Definition: gsForwardDeclarations.h:73
EIGEN_STRONG_INLINE diag_expr< E > const diagonal(E const &u)
Get diagonal elements of matrix as a vector.
Definition: gsExpressions.h:2082
Matrix< Scalar, Dynamic, Dynamic > cramerInverse() const
Inversion for small matrices using Cramer&#39;s Rule.
Definition: gsMatrixAddons.h:55
Definition: gsExpressions.h:3044
normalized_expr< E > normalized() const
Returns the vector normalized to unit length.
Definition: gsExpressions.h:283
The functions compute Dirichlet degrees of freedom using various methods.
EIGEN_STRONG_INLINE matrix_by_space_expr_tr< E1, E2 > const matrix_by_space_tr(E1 const &u, E2 const &v)
Definition: gsExpressions.h:4636
sqNorm_expr< E > sqNorm() const
Returns the squared Euclidean norm of the expression.
Definition: gsExpressions.h:291
index_t rows() const
Returns the row-size of the expression.
Definition: gsExpressions.h:328
Definition: gsExpressions.h:2304
ppart_expr< E > ppart() const
Returns the expression&#39;s positive part.
Definition: gsExpressions.h:253
const gsFeSpace< Scalar > & rowVar() const
Definition: gsExpressions.h:358
Outward normal on the boundary.
Definition: gsForwardDeclarations.h:86
EIGEN_STRONG_INLINE tangent_expr< T > tv(const gsGeometryMap< T > &u)
The tangent boundary vector of a geometry map in 2D.
Definition: gsExpressions.h:4513
norm_expr< E > norm() const
Returns the Euclidean norm of the expression.
Definition: gsExpressions.h:279
adjugate_expr< E > adj() const
Returns the adjugate of the expression (for matrix-valued expressions)
Definition: gsExpressions.h:275
Definition: gsExpressions.h:3011
S give(S &x)
Definition: gsMemory.h:266
EIGEN_STRONG_INLINE dJacG_expr< T > dJac(const gsGeometryMap< T > &G)
The partial derivatives of the Jacobian matrix of a geometry map.
Definition: gsExpressions.h:4551
Second derivatives.
Definition: gsForwardDeclarations.h:80
#define index_t
Definition: gsConfig.h:32
#define GISMO_ENSURE(cond, message)
Definition: gsDebug.h:102
const gsFeElement< Scalar > & _e
Reference to the element.
Definition: gsExpressions.h:856
EIGEN_STRONG_INLINE fform2nd_expr< T > fform2nd(const gsGeometryMap< T > &G)
The second fundamental form of G.
Definition: gsExpressions.h:4523
EIGEN_STRONG_INLINE frprod_expr< E1, E2 > const operator%(_expr< E1 > const &u, _expr< E2 > const &v)
Frobenious product (also known as double dot product) operator for expressions.
Definition: gsExpressions.h:4587
det_expr< E > det() const
Returns the determinant of the expression.
Definition: gsExpressions.h:287
#define GISMO_ASSERT(cond, message)
Definition: gsDebug.h:89
EIGEN_STRONG_INLINE summ_expr< E1, E2 > const summ(E1 const &u, E2 const &M)
Matrix-summation operator for expressions.
Definition: gsExpressions.h:4624
Maintains a mapping from patch-local dofs to global dof indices and allows the elimination of individ...
Definition: gsDofMapper.h:68
Definition: gsExpressions.h:3132
Struct containing information for matrix assembly.
Definition: gsExpressions.h:63
colsum_expr< E > colSum() const
Returns the colSum of a matrix.
Definition: gsExpressions.h:321
Laplacian.
Definition: gsForwardDeclarations.h:83
const gsFeSpace< Scalar > & colVar() const
Definition: gsExpressions.h:366
void finalize()
Must be called after all boundaries and interfaces have been marked to set up the dof numbering...
Definition: gsDofMapper.cpp:240
value_expr< E > val() const
Definition: gsExpressions.h:305
temp_expr< E > temp() const
Returns an evaluation of the (sub-)expression in temporary memory.
Definition: gsExpressions.h:263
exp_expr< E > exp() const
Returns exp(expression)
Definition: gsExpressions.h:249
#define gsWarn
Definition: gsDebug.h:50
Definition: gsExpressions.h:2335
Holds a set of patch-wise bases and their topology information.
Definition: gsMultiBasis.h:36
Normal vector of the object.
Definition: gsForwardDeclarations.h:85
inv_expr< E > const inv() const
Returns the inverse of the expression (for matrix-valued expressions)
Definition: gsExpressions.h:267
Provides declaration of Basis abstract interface.
EIGEN_STRONG_INLINE mult_expr< E1, E2 > const operator*(_expr< E1 > const &u, _expr< E2 > const &v)
Multiplication operator for expressions.
Definition: gsExpressions.h:4559
EIGEN_STRONG_INLINE onormal_expr< T > nv(const gsGeometryMap< T > &u)
The (outer pointing) boundary normal of a geometry map.
Definition: gsExpressions.h:4506
Definition: gsExpressions.h:139
Definition: gsExpressions.h:140
cb_expr< E > cb() const
Returns the puts the expression to colBlocks.
Definition: gsExpressions.h:241
tr_expr< E > tr() const
Returns the transpose of the expression.
Definition: gsExpressions.h:233
void print(std::ostream &os) const
Prints the expression as a string to os.
Definition: gsExpressions.h:194
MatExprType eval(const index_t k) const
Evaluates the expression at evaluation point indexed by k.
Definition: gsExpressions.h:229
Active function ids.
Definition: gsForwardDeclarations.h:84
Interface for the set of functions defined on a domain (the total number of functions in the set equa...
Definition: gsFuncData.h:23
EIGEN_STRONG_INLINE nabla_expr< T > nabla(const gsFeVariable< T > &u)
The nabla ( ) of a finite element variable.
Definition: gsExpressions.h:4502
Definition: gsExpressions.h:138
Second fundamental form.
Definition: gsForwardDeclarations.h:87
Definition: gsExpressions.h:126
The density of the measure pull back.
Definition: gsForwardDeclarations.h:76
trace_expr< E > trace() const
Returns the trace of the expression (for matrix-valued expressions)
Definition: gsExpressions.h:271
tr_expr< E, true > cwisetr() const
Returns the coordinate-wise transpose of the expression.
Definition: gsExpressions.h:237
EIGEN_STRONG_INLINE curl_expr< T > curl(const gsFeVariable< T > &u)
The curl of a finite element variable.
Definition: gsExpressions.h:4498
Definition: gsExpressions.h:2598
size_t numPatches() const
Returns the number of patches present underneath the mapper.
Definition: gsDofMapper.h:469
EIGEN_STRONG_INLINE reshape_expr< E > const reshape(E const &u, index_t n, index_t m)
Reshape an expression.
Definition: gsExpressions.h:1925
EIGEN_STRONG_INLINE dJacdc_expr< E > dJacdc(const E &u, index_t c)
The derivative of the jacobian of a geometry map with respect to a coordinate.
Definition: gsExpressions.h:4494
nabla2_expr< T > nabla2(const gsFeVariable< T > &u)
The nabla2 ( ) of a finite element variable.
Definition: gsExpressions.h:3003
Definition: gsExpressions.h:136
EIGEN_STRONG_INLINE flat_expr< E > const flat(E const &u)
Make a matrix 2x2 expression &quot;flat&quot;.
Definition: gsExpressions.h:2031
Value of the object.
Definition: gsForwardDeclarations.h:72
EIGEN_STRONG_INLINE grad_expr< E > grad(const E &u)
The gradient of a variable.
Definition: gsExpressions.h:4490
index_t cols() const
Returns the column-size of the expression.
Definition: gsExpressions.h:332
EIGEN_STRONG_INLINE idMat_expr id(const index_t dim)
The identity matrix of dimension dim.
Definition: gsExpressions.h:4470
#define GISMO_ERROR(message)
Definition: gsDebug.h:118
EIGEN_STRONG_INLINE meas_expr< T > meas(const gsGeometryMap< T > &G)
The measure of a geometry map.
Definition: gsExpressions.h:4555
EIGEN_STRONG_INLINE abs_expr< E > abs(const E &u)
Absolute value.
Definition: gsExpressions.h:4486
EIGEN_STRONG_INLINE divide_expr< E1, E2 > const operator/(_expr< E1 > const &u, _expr< E2 > const &v)
Scalar division operator for expressions.
Definition: gsExpressions.h:4592
Gradient transformation matrix.
Definition: gsForwardDeclarations.h:77
static bool isVector()
rowSpan &amp;&amp; !colSpan
Definition: gsExpressions.h:344
max_expr< E > max() const
Returns the rowSum of a matrix.
Definition: gsExpressions.h:313
GISMO_EXPR_VECTOR_EXPRESSION(norm, norm, 1)
Eucledian Norm.
void parse(gsExprHelper< Scalar > &evList) const
Parse the expression and discover the list of evaluation sources, also sets the required evaluation f...
Definition: gsExpressions.h:350
This object is a cache for computed values from an evaluator.
EIGEN_STRONG_INLINE add_expr< E1, E2 > const operator+(_expr< E1 > const &u, _expr< E2 > const &v)
Addition operator for expressions.
Definition: gsExpressions.h:4607
rowsum_expr< E > rowSum() const
Returns the rowSum of a matrix.
Definition: gsExpressions.h:317
Definition: gsExpressions.h:1978
EIGEN_STRONG_INLINE replicate_expr< E > const replicate(E const &u, index_t n, index_t m=1)
Replicate an expression.
Definition: gsExpressions.h:1967
Definition: gsExpressions.h:114
Defines expression precomputed_expr.
A basis represents a family of scalar basis functions defined over a common parameter domain...
Definition: gsBasis.h:78
sign_expr< E > sgn(Scalar tolerance=0) const
Returns the sign of the expression.
Definition: gsExpressions.h:245
Definition: gsExpressions.h:3080
bool isScalar() const
Returns true iff the expression is scalar-valued.
Definition: gsExpressions.h:342
Definition: gsExpressions.h:2542
EIGEN_STRONG_INLINE jac_expr< E > jac(const symbol_expr< E > &u)
The Jacobian matrix of a FE variable.
Definition: gsExpressions.h:4528