G+Smo  24.08.0
Geometry + Simulation Modules
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
gsThinShellUtils.h
Go to the documentation of this file.
1 
16 #pragma once
17 
19 #include <gsCore/gsLinearAlgebra.h>
20 #include <gsCore/gsBasis.h>
21 #include <gsCore/gsFuncData.h>
22 #include <gsCore/gsDofMapper.h>
23 
24 #include <gsPde/gsBoundaryConditions.h>
25 
26 #include <gsUtils/gsPointGrid.h>
27 
28 #include <gsAssembler/gsAssemblerOptions.h>
29 #include <gsAssembler/gsExpressions.h>
30 
31 
32 # define MatExprType auto
33 
34 namespace gismo{
35 namespace expr{
36 
40 class unitVec_expr : public _expr<unitVec_expr >
41 {
42 public:
43  typedef real_t Scalar;
44 private:
45  index_t _index;
46  index_t _dim;
47 
48 public:
49  unitVec_expr(const index_t index, const index_t dim) : _index(index), _dim(dim) { }
50 
51 public:
52  enum{ Space = 0, ScalarValued= 0, ColBlocks= 0};
53 
54  gsMatrix<Scalar> eval(const index_t) const
55  {
56  gsMatrix<Scalar> vec = gsMatrix<Scalar>::Zero(_dim,1);
57  // vec.setZero();
58  vec(_index,0) = 1;
59  return vec;
60  }
61 
62  index_t rows() const { return _dim; }
63  index_t cols() const { return 1; }
64  void parse(gsExprHelper<Scalar> & evList) const
65  { }
66 
67  const gsFeSpace<Scalar> & rowVar() const {return gsNullExpr<Scalar>::get();}
68  const gsFeSpace<Scalar> & colVar() const {return gsNullExpr<Scalar>::get();}
69 
70  void print(std::ostream &os) const { os << "uv("<<_dim <<")";}
71 };
72 
73 
79 template<class E>
80 class var1_expr : public _expr<var1_expr<E> >
81 {
82 public:
83  typedef typename E::Scalar Scalar;
84 
85 private:
86  typename E::Nested_t _u;
87  typename gsGeometryMap<Scalar>::Nested_t _G;
88 
89  mutable gsMatrix<Scalar> res;
90  mutable gsMatrix<Scalar> bGrads, cJac;
91  mutable gsVector<Scalar,3> m_v, normal;
92 public:
93  enum{ Space = E::Space, ScalarValued= 0, ColBlocks= 0};
94 
95  var1_expr(const E & u, const gsGeometryMap<Scalar> & G) : _u(u), _G(G) { }
96 
97 # define Eigen gsEigen
98  EIGEN_MAKE_ALIGNED_OPERATOR_NEW
99 # undef Eigen
100 
101  // helper function
102  static inline gsVector<Scalar,3> vecFun(index_t pos, Scalar val)
103  {
104  gsVector<Scalar,3> result = gsVector<Scalar,3>::Zero();
105  result[pos] = val;
106  return result;
107  }
108 
109  const gsMatrix<Scalar> & eval(const index_t k) const {return eval_impl(_u,k); }
110 
111  index_t rows() const { return 1; }
112  index_t cols() const { return cols_impl(_u); }
113 
114  void parse(gsExprHelper<Scalar> & evList) const
115  {
116  parse_impl<E>(evList);
117  }
118 
119  const gsFeSpace<Scalar> & rowVar() const { return _u.rowVar(); }
120  const gsFeSpace<Scalar> & colVar() const {return gsNullExpr<Scalar>::get();}
121  index_t cardinality_impl() const { return _u.cardinality_impl(); }
122 
123  void print(std::ostream &os) const { os << "var1("; _u.print(os); os <<")"; }
124 
125 private:
126  template<class U> inline
127  typename util::enable_if< util::is_same<U,gsFeSpace<Scalar> >::value,void>::type
128  parse_impl(gsExprHelper<Scalar> & evList) const
129  {
130  evList.add(_u);
131  _u.data().flags |= NEED_ACTIVE | NEED_GRAD; // need actives for cardinality
132  evList.add(_G);
133  _G.data().flags |= NEED_NORMAL | NEED_DERIV;
134  }
135 
136  template<class U> inline
137  typename util::enable_if< util::is_same<U,gsFeSolution<Scalar>>::value || util::is_same<U,gsFeVariable<Scalar>>::value,void>::type
138  parse_impl(gsExprHelper<Scalar> & evList) const
139  {
140  evList.add(_G);
141  _G.data().flags |= NEED_NORMAL | NEED_DERIV;
142 
143  grad(_u).parse(evList); //
144 
145  _u.parse(evList);
146  }
147 
148  template<class U> inline
149  typename util::enable_if< util::is_same<U,gsFeSpace<Scalar> >::value, const gsMatrix<Scalar> & >::type
150  eval_impl(const U & u, const index_t k) const
151  {
152  const index_t A = _u.cardinality()/_u.dim(); // _u.data().actives.rows()
153  res.resize(_u.cardinality(), cols()); // rows()*
154 
155  normal = _G.data().normal(k);// not normalized to unit length
156  normal.normalize();
157 
158  bGrads = _u.data().values[1].col(k);
159  cJac = _G.data().values[1].reshapeCol(k, _G.data().dim.first, _G.data().dim.second).transpose();
160  const Scalar measure = (cJac.col3d(0).cross( cJac.col3d(1) )).norm();
161 
162  for (index_t d = 0; d!= cols(); ++d) // for all basis function components
163  {
164  const short_t s = d*A;
165  for (index_t j = 0; j!= A; ++j) // for all actives
166  {
167  // Jac(u) ~ Jac(G) with alternating signs ?..
168  m_v.noalias() = (vecFun(d, bGrads.at(2*j ) ).cross( cJac.col3d(1) )
169  - vecFun(d, bGrads.at(2*j+1) ).cross( cJac.col3d(0) )) / measure;
170 
171  // --------------- First variation of the normal
172  res.row(s+j).noalias() = (m_v - ( normal*m_v.transpose() ) * normal).transpose(); // outer-product version
173  }
174  }
175  return res;
176  }
177 
178  template<class U> inline
179  typename util::enable_if< util::is_same<U,gsFeVariable<Scalar> >::value, const gsMatrix<Scalar> & >::type
180  eval_impl(const U & u, const index_t k) const
181  {
182  grad_expr<U> vGrad = grad_expr<U>(u);
183  bGrads = vGrad.eval(k);
184 
185  return make_vector(k);
186  }
187 
188  // Specialisation for a solution
189  template<class U> inline
190  typename util::enable_if< util::is_same<U,gsFeSolution<Scalar> >::value, const gsMatrix<Scalar> & >::type
191  eval_impl(const U & u, const index_t k) const
192  {
193  GISMO_ASSERT(1==_u.data().actives.cols(), "Single actives expected");
194  grad_expr<gsFeSolution<Scalar>> sGrad = grad_expr<gsFeSolution<Scalar>>(_u);
195  bGrads = sGrad.eval(k);
196 
197  return make_vector(k);
198  }
199 
200  const gsMatrix<Scalar> & make_vector(const index_t k) const
201  {
202  res.resize(rows(), cols()); // rows()*
203  normal = _G.data().normal(k);// not normalized to unit length
204  normal.normalize();
205 
206  cJac = _G.data().values[1].reshapeCol(k, _G.data().dim.first, _G.data().dim.second).transpose();
207  const Scalar measure = (cJac.col3d(0).cross( cJac.col3d(1) )).norm();
208 
209  m_v.noalias() = ( ( bGrads.col3d(0) ).cross( cJac.col3d(1) )
210  - ( bGrads.col3d(1) ).cross( cJac.col3d(0) ) ) / measure;
211 
212  // --------------- First variation of the normal
213  res = (m_v - ( normal*m_v.transpose() ) * normal).transpose(); // outer-product version
214  return res;
215  }
216 
217  template<class U> inline
218  typename util::enable_if< util::is_same<U,gsGeometryMap<Scalar> >::value, index_t >::type
219  cols_impl(const U & u) const
220  {
221  return _u.data().dim.second;
222  }
223 
224  template<class U> inline
225  typename util::enable_if<util::is_same<U,gsFeSpace<Scalar> >::value || util::is_same<U,gsFeSolution<Scalar> >::value, index_t >::type
226  cols_impl(const U & u) const
227  {
228  return _u.dim();
229  }
230 
231  template<class U> inline
232  typename util::enable_if<util::is_same<U,gsFeVariable<Scalar> >::value, index_t >::type
233  cols_impl(const U & u) const
234  {
235  return _u.source().targetDim();
236  }
237 
238 };
239 
247 template<class E1, class E2, class E3>
248 class var2dot_expr : public _expr<var2dot_expr<E1,E2,E3> >
249 {
250 public:
251  typedef typename E1::Scalar Scalar;
252 
253 private:
254  typename E1::Nested_t _u;
255  typename E2::Nested_t _v;
256  typename gsGeometryMap<Scalar>::Nested_t _G;
257  typename E3::Nested_t _Ef;
258 
259 public:
260  enum{ Space = 3, ScalarValued= 0, ColBlocks= 0 };
261 
262  var2dot_expr( const E1 & u, const E2 & v, const gsGeometryMap<Scalar> & G, _expr<E3> const& Ef) : _u(u),_v(v), _G(G), _Ef(Ef) { }
263 
264  mutable gsMatrix<Scalar> res;
265 
266  mutable gsMatrix<Scalar> uGrads, vGrads, cJac, cDer2, evEf, result;
267  mutable gsVector<Scalar> m_u, m_v, normal, m_uv, m_u_der, n_der, n_der2, tmp; // memomry leaks when gsVector<T,3>, i.e. fixed dimension
268 # define Eigen gsEigen
269  EIGEN_MAKE_ALIGNED_OPERATOR_NEW
270 # undef Eigen
271 
272  // helper function
273  static inline gsVector<Scalar,3> vecFun(index_t pos, Scalar val)
274  {
275  gsVector<Scalar,3> result = gsVector<Scalar,3>::Zero();
276  result[pos] = val;
277  return result;
278  }
279 
280  const gsMatrix<Scalar> & eval(const index_t k) const
281  {
282  res.resize(_u.cardinality(), _u.cardinality());
283 
284  normal = _G.data().normal(k);
285  normal.normalize();
286  uGrads = _u.data().values[1].col(k);
287  vGrads = _v.data().values[1].col(k);
288  cJac = _G.data().values[1].reshapeCol(k, _G.data().dim.first, _G.data().dim.second).transpose();
289 
290  const index_t cardU = _u.data().values[0].rows(); // number of actives per component of u
291  const index_t cardV = _v.data().values[0].rows(); // number of actives per component of v
292  const Scalar measure = (cJac.col3d(0).cross( cJac.col3d(1) )).norm();
293 
294  evEf = _Ef.eval(k);
295 
296  for (index_t j = 0; j!= cardU; ++j) // for all basis functions u (1)
297  {
298  for (index_t i = 0; i!= cardV; ++i) // for all basis functions v (1)
299  {
300  for (index_t d = 0; d!= _u.dim(); ++d) // for all basis functions u (2)
301  {
302  m_u.noalias() = ( vecFun(d, uGrads.at(2*j ) ).cross( cJac.col3d(1) )
303  -vecFun(d, uGrads.at(2*j+1) ).cross( cJac.col3d(0) ))
304  / measure;
305 
306  const short_t s = d*cardU;
307 
308  for (index_t c = 0; c!= _v.dim(); ++c) // for all basis functions v (2)
309  {
310  const short_t r = c*cardV;
311  m_v.noalias() = ( vecFun(c, vGrads.at(2*i ) ).cross( cJac.col3d(1) )
312  -vecFun(c, vGrads.at(2*i+1) ).cross( cJac.col3d(0) ))
313  / measure;
314 
315  // n_der.noalias() = (m_v - ( normal.dot(m_v) ) * normal);
316  n_der.noalias() = (m_v - ( normal*m_v.transpose() ) * normal); // outer-product version
317 
318  m_uv.noalias() = ( vecFun(d, uGrads.at(2*j ) ).cross( vecFun(c, vGrads.at(2*i+1) ) )
319  +vecFun(c, vGrads.at(2*i ) ).cross( vecFun(d, uGrads.at(2*j+1) ) ))
320  / measure; //check
321 
322  m_u_der.noalias() = (m_uv - ( normal.dot(m_v) ) * m_u);
323  // m_u_der.noalias() = (m_uv - ( normal*m_v.transpose() ) * m_u); // outer-product version TODO
324 
325  // --------------- Second variation of the normal
326  tmp = m_u_der - (m_u.dot(n_der) + normal.dot(m_u_der) ) * normal - (normal.dot(m_u) ) * n_der;
327  // tmp = m_u_der - (m_u.dot(n_der) + normal.dot(m_u_der) ) * normal - (normal.dot(m_u) ) * n_der;
328 
329  // Evaluate the product
330  result = evEf * tmp;
331 
332  res(s + j, r + i ) = result(0,0);
333  }
334  }
335  }
336  }
337  return res;
338  }
339 
340  index_t rows() const
341  {
342  return 1; // because the resulting matrix has scalar entries for every combination of active basis functions
343  }
344 
345  index_t cols() const
346  {
347  return 1; // because the resulting matrix has scalar entries for every combination of active basis functions
348  }
349 
350  void parse(gsExprHelper<Scalar> & evList) const
351  {
352  evList.add(_u);
353  _u.data().flags |= NEED_ACTIVE | NEED_VALUE | NEED_GRAD;
354  evList.add(_v);
355  _v.data().flags |= NEED_ACTIVE | NEED_VALUE | NEED_GRAD;
356  evList.add(_G);
357  _G.data().flags |= NEED_NORMAL | NEED_DERIV;
358  _Ef.parse(evList);
359  }
360 
361  const gsFeSpace<Scalar> & rowVar() const { return _u.rowVar(); }
362  const gsFeSpace<Scalar> & colVar() const { return _v.rowVar(); }
363 
364  void print(std::ostream &os) const { os << "var2("; _u.print(os); os <<")"; }
365 };
366 
374 template<class E1, class E2, class E3>
375 class var2deriv2dot_expr : public _expr<var2deriv2dot_expr<E1,E2,E3> >
376 {
377 public:
378  typedef typename E1::Scalar Scalar;
379 
380 private:
381  typename E1::Nested_t _u;
382  typename E2::Nested_t _v;
383  typename gsGeometryMap<Scalar>::Nested_t _G;
384  typename E3::Nested_t _Ef;
385 
386 public:
387  enum{ Space = 3, ScalarValued= 0, ColBlocks= 0 };
388 
389  var2deriv2dot_expr( const E1 & u, const E2 & v, const gsGeometryMap<Scalar> & G, _expr<E3> const& Ef) : _u(u),_v(v), _G(G), _Ef(Ef) { }
390 
391  mutable gsMatrix<Scalar> res;
392 
393  mutable gsMatrix<Scalar> uGrads, vGrads, cJac, cDer2, evEf, result;
394  mutable gsVector<Scalar> m_u, m_v, normal, m_uv, m_u_der, n_der, n_der2, tmp; // memomry leaks when gsVector<T,3>, i.e. fixed dimension
395 # define Eigen gsEigen
396  EIGEN_MAKE_ALIGNED_OPERATOR_NEW
397 # undef Eigen
398 
399  // helper function
400  static inline gsVector<Scalar,3> vecFun(index_t pos, Scalar val)
401  {
402  gsVector<Scalar,3> result = gsVector<Scalar,3>::Zero();
403  result[pos] = val;
404  return result;
405  }
406 
407  const gsMatrix<Scalar> & eval(const index_t k) const
408  {
409  res.resize(_u.cardinality(), _u.cardinality());
410 
411  normal = _G.data().normal(k);
412  normal.normalize();
413  uGrads = _u.data().values[1].col(k);
414  vGrads = _v.data().values[1].col(k);
415  cJac = _G.data().values[1].reshapeCol(k, _G.data().dim.first, _G.data().dim.second).transpose();
416  cDer2 = _G.data().values[2].reshapeCol(k, _G.data().dim.second, _G.data().dim.second);
417 
418  const index_t cardU = _u.data().values[0].rows(); // number of actives per component of u
419  const index_t cardV = _v.data().values[0].rows(); // number of actives per component of v
420  const Scalar measure = (cJac.col3d(0).cross( cJac.col3d(1) )).norm();
421 
422  evEf = _Ef.eval(k);
423 
424  for (index_t j = 0; j!= cardU; ++j) // for all basis functions u (1)
425  {
426  for (index_t i = 0; i!= cardV; ++i) // for all basis functions v (1)
427  {
428  for (index_t d = 0; d!= _u.dim(); ++d) // for all basis functions u (2)
429  {
430  m_u.noalias() = ( vecFun(d, uGrads.at(2*j ) ).cross( cJac.col3d(1) )
431  -vecFun(d, uGrads.at(2*j+1) ).cross( cJac.col3d(0) ))
432  / measure;
433 
434  const short_t s = d*cardU;
435 
436  for (index_t c = 0; c!= _v.dim(); ++c) // for all basis functions v (2)
437  {
438  const short_t r = c*cardV;
439  m_v.noalias() = ( vecFun(c, vGrads.at(2*i ) ).cross( cJac.col3d(1) )
440  -vecFun(c, vGrads.at(2*i+1) ).cross( cJac.col3d(0) ))
441  / measure;
442 
443  // n_der.noalias() = (m_v - ( normal.dot(m_v) ) * normal);
444  n_der.noalias() = (m_v - ( normal*m_v.transpose() ) * normal); // outer-product version
445 
446  m_uv.noalias() = ( vecFun(d, uGrads.at(2*j ) ).cross( vecFun(c, vGrads.at(2*i+1) ) )
447  +vecFun(c, vGrads.at(2*i ) ).cross( vecFun(d, uGrads.at(2*j+1) ) ))
448  / measure; //check
449 
450  m_u_der.noalias() = (m_uv - ( normal.dot(m_v) ) * m_u);
451  // m_u_der.noalias() = (m_uv - ( normal*m_v.transpose() ) * m_u); // outer-product version TODO
452 
453  // --------------- Second variation of the normal
454  tmp = m_u_der - (m_u.dot(n_der) + normal.dot(m_u_der) ) * normal - (normal.dot(m_u) ) * n_der;
455  // tmp = m_u_der - (m_u.dot(n_der) + normal.dot(m_u_der) ) * normal - (normal.dot(m_u) ) * n_der;
456 
457  // Evaluate the product
458  tmp = cDer2 * tmp; // E_f_der2, last component
459  tmp.row(2) *= 2.0;
460  result = evEf * tmp;
461 
462  res(s + j, r + i ) = result(0,0);
463  }
464  }
465  }
466  }
467  return res;
468  }
469 
470  index_t rows() const
471  {
472  return 1; // because the resulting matrix has scalar entries for every combination of active basis functions
473  }
474 
475  index_t cols() const
476  {
477  return 1; // because the resulting matrix has scalar entries for every combination of active basis functions
478  }
479 
480  void parse(gsExprHelper<Scalar> & evList) const
481  {
482  evList.add(_u);
483  _u.data().flags |= NEED_ACTIVE | NEED_VALUE | NEED_GRAD;
484  evList.add(_v);
485  _v.data().flags |= NEED_ACTIVE | NEED_VALUE | NEED_GRAD;
486  evList.add(_G);
487  _G.data().flags |= NEED_NORMAL | NEED_DERIV | NEED_2ND_DER;
488  _Ef.parse(evList);
489  }
490 
491  const gsFeSpace<Scalar> & rowVar() const { return _u.rowVar(); }
492  const gsFeSpace<Scalar> & colVar() const { return _v.rowVar(); }
493 
494  void print(std::ostream &os) const { os << "var2dot("; _u.print(os), _v.print(os), _G.print(os), _Ef.print(os); os <<")"; }
495 };
496 
498 template<class E1, class E2>
499 class var2_expr : public _expr<var2_expr<E1,E2> >
500 {
501 public:
502  typedef typename E1::Scalar Scalar;
503 
504 private:
505 
506  typename E1::Nested_t _u;
507  typename E2::Nested_t _v;
508  typename gsGeometryMap<Scalar>::Nested_t _G;
509 
510 public:
511  enum{ Space = E1::Space, ScalarValued= 0, ColBlocks= 0 };
512 
513  var2_expr( const E1 & u, const E2 & v, const gsGeometryMap<Scalar> & G) : _u(u),_v(v), _G(G) { }
514 
515  mutable gsMatrix<Scalar> res;
516 
517  mutable gsMatrix<Scalar> uGrads, vGrads, cJac, cDer2, result;
518  mutable gsVector<Scalar> m_u, m_v, normal, m_uv, m_u_der, n_der, n_der2, tmp; // memomry leaks when gsVector<T,3>, i.e. fixed dimension
519 
520 # define Eigen gsEigen
521  EIGEN_MAKE_ALIGNED_OPERATOR_NEW
522 # undef Eigen
523 
524  // helper function
525  static inline gsVector<Scalar,3> vecFun(index_t pos, Scalar val)
526  {
527  gsVector<Scalar,3> result = gsVector<Scalar,3>::Zero();
528  result[pos] = val;
529  return result;
530  }
531 
532  const gsMatrix<Scalar> & eval(const index_t k) const {return eval_impl(_u,_v,k); }
533 
534  index_t rows() const
535  {
536  return 3;
537  }
538 
539  index_t cols() const
540  {
541  return 1;
542  }
543 
544  void parse(gsExprHelper<Scalar> & evList) const
545  {
546  grad(_u).parse(evList); //
547  grad(_v).parse(evList); //
548 
549  _u.parse(evList);
550  evList.add(_u);
551  _u.data().flags |= NEED_GRAD | NEED_ACTIVE;
552 
553  _v.parse(evList);
554  evList.add(_v);
555  _v.data().flags |= NEED_GRAD | NEED_ACTIVE;
556 
557  _G.parse(evList);
558  evList.add(_G);
559  _G.data().flags |= NEED_NORMAL | NEED_DERIV | NEED_DERIV2 ;
560  }
561 
562  const gsFeSpace<Scalar> & rowVar() const { return _u.rowVar(); }
563  const gsFeSpace<Scalar> & colVar() const { return _v.rowVar(); }
564 
565  enum{rowSpan = 0, colSpan = 0};
566 
567  void print(std::ostream &os) const { os << "var2("; _u.print(os), _v.print(os), _G.print(os); os <<")"; }
568 
569 private:
570 // Specialisation for a space
571  template<class U, class V> inline
572  typename util::enable_if< !(util::is_same<U,gsFeVariable<Scalar> >::value && util::is_same<V,gsFeVariable<Scalar> >::value), const gsMatrix<Scalar> & >::type
573  eval_impl(const U & u, const V & v, const index_t k) const
574  {
575  GISMO_NO_IMPLEMENTATION;
576  }
577 
578  template<class U, class V> inline
579  typename util::enable_if< util::is_same<U,gsFeVariable<Scalar> >::value && util::is_same<V,gsFeVariable<Scalar> >::value, const gsMatrix<Scalar> & >::type
580  eval_impl(const U & u, const V & v, const index_t k) const
581  {
582  grad_expr<U> uGrad = grad_expr<U>(u);
583  uGrads = uGrad.eval(k);
584  grad_expr<V> vGrad = grad_expr<V>(v);
585  vGrads = vGrad.eval(k);
586 
587  res.resize(rows(), cols()); // rows()*
588 
589  normal = _G.data().normal(k);
590  normal.normalize();
591 
592  cJac = _G.data().values[1].reshapeCol(k, _G.data().dim.first, _G.data().dim.second).transpose();
593  cDer2 = _G.data().values[2].reshapeCol(k, _G.data().dim.second, _G.data().dim.second);
594 
595  const Scalar measure = (cJac.col3d(0).cross( cJac.col3d(1) )).norm();
596 
597  m_u.noalias() = ( ( uGrads.col(0).template head<3>() ).cross( cJac.col(1).template head<3>() )
598  - ( uGrads.col(1).template head<3>() ).cross( cJac.col(0).template head<3>() ) )
599  / measure;
600 
601  m_v.noalias() = ( ( vGrads.col(0).template head<3>() ).cross( cJac.col(1).template head<3>() )
602  - ( vGrads.col(1).template head<3>() ).cross( cJac.col(0).template head<3>() ) )
603  / measure;
604 
605  // n_der.noalias() = (m_v - ( normal.dot(m_v) ) * normal);
606  n_der.noalias() = (m_v - ( normal*m_v.transpose() ) * normal); // outer-product version
607 
608  m_uv.noalias() = ( ( uGrads.col(0).template head<3>() ).cross( vGrads.col(1).template head<3>() )
609  - ( vGrads.col(1).template head<3>() ).cross( uGrads.col(0).template head<3>() ) )
610  / measure;
611 
612  m_u_der.noalias() = (m_uv - ( normal.dot(m_v) ) * m_u);
613  // m_u_der.noalias() = (m_uv - ( normal*m_v.transpose() ) * m_u); // outer-product version TODO
614 
615  // --------------- Second variation of the normal
616  tmp = m_u_der - (m_u.dot(n_der) + normal.dot(m_u_der) ) * normal - (normal.dot(m_u) ) * n_der;
617  // tmp = m_u_der - (m_u.dot(n_der) + normal.dot(m_u_der) ) * normal - (normal.dot(m_u) ) * n_der;
618 
619  // Evaluate the product
620  tmp = cDer2 * tmp; // E_f_der2, last component
621  tmp.row(2) *= 2.0;
622  res = tmp;
623  return res;
624  }
625 
626  // // Specialisation for a solution
627  // template<class U> inline
628  // typename util::enable_if< util::is_same<U,gsFeSolution<Scalar> >::value, const gsMatrix<Scalar> & >::type
629  // eval_impl(const U & u, const index_t k) const
630  // {
631  // GISMO_ASSERT(1==_u.data().actives.cols(), "Single actives expected");
632  // solGrad_expr<Scalar> sGrad = solGrad_expr<Scalar>(_u);
633  // res.resize(rows(), cols()); // rows()*
634 
635  // normal = _G.data().normal(k);// not normalized to unit length
636  // normal.normalize();
637  // bGrads = sGrad.eval(k);
638  // cJac = _G.data().values[1].reshapeCol(k, _G.data().dim.first, _G.data().dim.second).transpose();
639  // const Scalar measure = _G.data().measures.at(k);
640 
641  // m_v.noalias() = ( ( bGrads.col(0).template head<3>() ).cross( cJac.col(1).template head<3>() )
642  // - ( bGrads.col(1).template head<3>() ).cross( cJac.col(0).template head<3>() ) ) / measure;
643 
644  // // --------------- First variation of the normal
645  // // res.row(s+j).noalias() = (m_v - ( normal.dot(m_v) ) * normal).transpose();
646  // res = (m_v - ( normal*m_v.transpose() ) * normal).transpose(); // outer-product version
647  // return res;
648  // }
649 
650 };
651 
652 
658 template<class E>
659 class tvar1_expr : public _expr<tvar1_expr<E> >
660 {
661 public:
662  typedef typename E::Scalar Scalar;
663 
664 private:
665  typename E::Nested_t _u;
666  typename gsGeometryMap<Scalar>::Nested_t _G;
667 
668  mutable gsVector<Scalar,3> onormal, tangent, utangent, dtan;
669  mutable gsVector<Scalar> tmp;
670  mutable gsMatrix<Scalar> bGrads, cJac, res;
671 
672 public:
673  enum{ Space = E::Space, ScalarValued= 0, ColBlocks= 0};
674 
675  tvar1_expr(const E & u, const gsGeometryMap<Scalar> & G) : _u(u), _G(G) { }
676 
677 # define Eigen gsEigen
678  EIGEN_MAKE_ALIGNED_OPERATOR_NEW
679 # undef Eigen
680 
681  // helper function
682  static inline gsVector<Scalar,3> vecFun(index_t pos, Scalar val)
683  {
684  gsVector<Scalar,3> result = gsVector<Scalar,3>::Zero();
685  result[pos] = val;
686  return result;
687  }
688  const gsMatrix<Scalar> & eval(const index_t k) const {return eval_impl(_u,k); }
689 
690  index_t rows() const { return 1; }
691  index_t cols() const { return _u.dim(); }
692 
693  void parse(gsExprHelper<Scalar> & evList) const
694  {
695  parse_impl<E>(evList);
696  }
697 
698  const gsFeSpace<Scalar> & rowVar() const { return _u.rowVar(); }
699  const gsFeSpace<Scalar> & colVar() const {return gsNullExpr<Scalar>::get();}
700  index_t cardinality_impl() const { return _u.cardinality_impl(); }
701 
702  void print(std::ostream &os) const { os << "tvar("; _G.print(os); os <<")"; }
703 
704 private:
705  template<class U> inline
706  typename util::enable_if< !util::is_same<U,gsFeSolution<Scalar> >::value,void>::type
707  parse_impl(gsExprHelper<Scalar> & evList) const
708  {
709  evList.add(_u);
710  _u.data().flags |= NEED_GRAD | NEED_ACTIVE;
711  evList.add(_G);
712  _G.data().flags |= NEED_NORMAL | NEED_OUTER_NORMAL | NEED_DERIV;
713  }
714 
715  template<class U> inline
716  typename util::enable_if< util::is_same<U,gsFeSolution<Scalar> >::value,void>::type
717  parse_impl(gsExprHelper<Scalar> & evList) const
718  {
719  evList.add(_G);
720  _G.data().flags |= NEED_NORMAL | NEED_OUTER_NORMAL | NEED_DERIV;
721 
722  grad(_u).parse(evList); //
723 
724  _u.parse(evList);
725  }
726 
727  template<class U> inline
728  typename util::enable_if< !util::is_same<U,gsFeSolution<Scalar> >::value, const gsMatrix<Scalar> & >::type
729  eval_impl(const U & u, const index_t k) const
730  {
731  GISMO_ASSERT(_G.data().dim.second==3,"Domain dimension should be 3, is "<<_G.data().dim.second);
732 
733  const index_t A = _u.cardinality()/_u.dim();
734  res.resize(_u.cardinality(), cols()); // rows()*
735  res.setZero();
736  cJac = _G.data().jacobian(k);
737  cJac.colwise().normalize();
738 
739  onormal = _G.data().outNormal(k);
740  onormal.normalize();
741  tmp = cJac.transpose() * onormal;
742  tmp.colwise().normalize(); //normalize the inner product for fair comparison
743 
744  Scalar tol = 1e-8;
745  /*
746  We can check which column of the Jacobian corresponds to the outer normal vector or to the tangent.
747  The tangent is a covariant vector and hence the column of the Jacobian should be equal to the tangent.
748  The normal is a contravariant vector and hence the corresponding column of the Jacobian times the outward normal should give 1. We use this property.
749  */
750  index_t colIndex = -1;
751  if ( (math::abs(tmp.at(0)) < tol) && (math::abs(tmp.at(1)) > 1-tol ) ) // then the normal is vector 2 and the tangent vector 1
752  colIndex = 0;
753  else if ( (math::abs(tmp.at(1)) < tol) && (math::abs(tmp.at(0)) > 1-tol ) ) // then the normal is vector 1 and the tangent vector 2
754  colIndex = 1;
755  else // then the normal is unknown??
756  gsInfo<<"warning: choice unknown\n";
757 
758  // tangent = cJac.col(colIndex);
759  tangent_expr<Scalar> tan_expr = tangent_expr<Scalar>(_G);
760  tangent = tan_expr.eval(k);
761  utangent = tangent.normalized();
762 
763  index_t sign = 0;
764  if (colIndex!=-1)
765  {
766  Scalar dot = tangent.dot(cJac.col(colIndex));
767  sign = (Scalar(0) < dot) - (dot < Scalar(0));
768  }
769  else
770  {
771  gsWarn<<"No suitable tangent and outer normal vector found for point "<<_G.data().values[0].transpose()<<"\n";
772  return res;
773  }
774 
775  // Now we will compute the derivatives of the basis functions
776  bGrads = _u.data().values[1].col(k);
777  for (index_t d = 0; d!= cols(); ++d) // for all basis function components
778  {
779  const short_t s = d*A;
780  for (index_t j = 0; j!= A; ++j) // for all actives
781  {
782  // The tangent vector is in column colIndex in cJac and thus in 2*j+colIndex in bGrads.
783  // Furthermore, as basis function for dimension d, it has a nonzero in entry d, and zeros elsewhere
784  dtan = sign*vecFun(d, bGrads.at(2*j+colIndex));
785  res.row(s+j).noalias() = (1 / tangent.norm() * ( dtan - ( utangent.transpose() * dtan ) * utangent ) ).transpose();
786  }
787  }
788  return res;
789  }
790 
791  template<class U> inline
792  typename util::enable_if< util::is_same<U,gsFeSolution<Scalar> >::value, const gsMatrix<Scalar> & >::type
793  eval_impl(const U & u, const index_t k) const
794  {
795  GISMO_ASSERT(1==_u.data().actives.cols(), "Single actives expected");
796  res.resize(rows(), cols());
797  res.setZero();
798 
799  cJac = _G.data().jacobian(k);
800  cJac.colwise().normalize();
801 
802  onormal = _G.data().outNormal(k);
803  onormal.normalize();
804  tmp = cJac.transpose() * onormal;
805  tmp.colwise().normalize(); //normalize the inner product for fair comparison
806 
807  Scalar tol = 1e-8;
808 
809  /*
810  We can check which column of the Jacobian corresponds to the outer normal vector or to the tangent.
811  The tangent is a covariant vector and hence the column of the Jacobian should be equal to the tangent.
812  The normal is a contravariant vector and hence the corresponding column of the Jacobian times the outward normal should give 1. We use this property.
813  */
814  index_t colIndex = -1;
815  if ( (math::abs(tmp.at(0)) < tol) && (math::abs(tmp.at(1)) > 1-tol ) ) // then the normal is vector 2 and the tangent vector 1
816  colIndex = 0;
817  else if ( (math::abs(tmp.at(1)) < tol) && (math::abs(tmp.at(0)) > 1-tol ) ) // then the normal is vector 1 and the tangent vector 2
818  colIndex = 1;
819  else // then the normal is unknown??
820  gsInfo<<"warning: choice unknown\n";
821 
822  // tangent = cJac.col(colIndex);
823  tangent_expr<Scalar> tan_expr = tangent_expr<Scalar>(_G);
824  tangent = tan_expr.eval(k);
825  // utangent = tangent.normalized();
826 
827  index_t sign = 0;
828  if (colIndex!=-1)
829  {
830  Scalar dot = tangent.dot(cJac.col(colIndex));
831  sign = (Scalar(0) < dot) - (dot < Scalar(0));
832  }
833  else
834  {
835  gsWarn<<"No suitable tangent and outer normal vector found for point "<<_G.data().values[0].transpose()<<"\n";
836  return res;
837  }
838 
839  bGrads = _u.data().values[1].col(k);
840  dtan = sign*bGrads.col(colIndex);
841  res.noalias() = (1 / tangent.norm() * ( dtan - ( tangent * dtan ) * tangent / (tangent.norm() * tangent.norm()) )).transpose();
842  return res;
843  }
844 
845 };
846 
852 template<class E>
853 class ovar1_expr : public _expr<ovar1_expr<E> >
854 {
855 public:
856  typedef typename E::Scalar Scalar;
857 
858 private:
859  typename E::Nested_t _u;
860  typename gsGeometryMap<Scalar>::Nested_t _G;
861 
862  mutable gsVector<Scalar,3> tangent, normal, tvar, snvar;
863 
864  mutable gsMatrix<Scalar> tvarMat, snvarMat, res;
865 
866 public:
867  enum{ Space = E::Space, ScalarValued= 0, ColBlocks= 0};
868 
869  ovar1_expr(const E & u, const gsGeometryMap<Scalar> & G) : _u(u), _G(G) { }
870 
871 # define Eigen gsEigen
872  EIGEN_MAKE_ALIGNED_OPERATOR_NEW
873 # undef Eigen
874 
875  // helper function
876  static inline gsVector<Scalar,3> vecFun(index_t pos, Scalar val)
877  {
878  gsVector<Scalar,3> result = gsVector<Scalar,3>::Zero();
879  result[pos] = val;
880  return result;
881  }
882 
883  const gsMatrix<Scalar> & eval(const index_t k) const { return eval_impl(_u,k); }
884 
885  index_t rows() const { return 1; }
886  index_t cols() const { return _u.dim(); }
887 
888  void parse(gsExprHelper<Scalar> & evList) const
889  {
890  parse_impl<E>(evList);
891  }
892 
893  const gsFeSpace<Scalar> & rowVar() const { return _u.rowVar(); }
894  const gsFeSpace<Scalar> & colVar() const {return gsNullExpr<Scalar>::get();}
895  index_t cardinality_impl() const { return _u.cardinality_impl(); }
896 
897  void print(std::ostream &os) const { os << "ovar("; _G.print(os); os <<")"; }
898 
899 private:
900  template<class U> inline
901  typename util::enable_if< !util::is_same<U,gsFeSolution<Scalar> >::value,void>::type
902  parse_impl(gsExprHelper<Scalar> & evList) const
903  {
904  evList.add(_u);
905  _u.data().flags |= NEED_GRAD | NEED_ACTIVE;
906  evList.add(_G);
907  _G.data().flags |= NEED_NORMAL | NEED_OUTER_NORMAL | NEED_DERIV; // all needed?
908 
909  tv(_G).parse(evList);
910  tvar1(_u,_G).parse(evList);
911  var1(_u,_G).parse(evList);
912  }
913 
914  template<class U> inline
915  typename util::enable_if< util::is_same<U,gsFeSolution<Scalar> >::value,void>::type
916  parse_impl(gsExprHelper<Scalar> & evList) const
917  {
918  evList.add(_G);
919  _G.data().flags |= NEED_NORMAL | NEED_OUTER_NORMAL | NEED_DERIV; // all needed?
920 
921  grad(_u).parse(evList); //
922  tv(_G).parse(evList);
923 
924  _u.parse(evList);
925  }
926 
927  template<class U> inline
928  typename util::enable_if< !util::is_same<U,gsFeSolution<Scalar> >::value, const gsMatrix<Scalar> & >::type
929  eval_impl(const U & u, const index_t k) const
930  {
931  GISMO_ASSERT(_G.data().dim.second==3,"Domain dimension should be 3, is "<<_G.data().dim.second);
932  const index_t A = _u.cardinality()/_u.dim(); // _u.data().actives.rows()
933  res.resize(_u.cardinality(), cols()); // rows()*
934 
935  tangent_expr<Scalar> tan_expr = tangent_expr<Scalar>(_G);
936  tangent = tan_expr.eval(k);
937  tangent.normalize();
938 
939  // For the normal vector variation
940  normal = _G.data().normal(k);
941  normal.normalize();
942 
943  tvar1_expr<E> tvar_expr = tvar1_expr<E>(_u,_G);
944  tvarMat = tvar_expr.eval(k);
945 
946  var1_expr<E> snvar_expr = var1_expr<E>(_u,_G);
947  snvarMat = snvar_expr.eval(k);
948 
949  for (index_t d = 0; d!= cols(); ++d) // for all basis function components
950  {
951  const short_t s = d*A;
952  for (index_t j = 0; j!= A; ++j) // for all actives
953  {
954  tvar = tvarMat.row(s+j);
955  snvar= snvarMat.row(s+j);
956 
957  // VARIATION OF THE OUTER NORMAL
958  res.row(s+j).noalias() = tvar.cross(normal) + tangent.cross(snvar);
959  }
960  }
961  return res;
962  }
963 
964  template<class U> inline
965  typename util::enable_if< util::is_same<U,gsFeSolution<Scalar> >::value, const gsMatrix<Scalar> & >::type
966  eval_impl(const U & u, const index_t k) const
967  {
968  GISMO_ASSERT(_G.data().dim.second==3,"Domain dimension should be 3, is "<<_G.data().dim.second);
969  res.resize(rows(), cols());
970 
971  tangent_expr<Scalar> tan_expr = tangent_expr<Scalar>(_G);
972  tangent = tan_expr.eval(k);
973  tangent.normalize();
974 
975  // For the normal vector variation
976  normal = _G.data().normal(k);
977  normal.normalize();
978 
979  tvar1_expr<E> tvar_expr = tvar1_expr<E>(_u,_G);
980  tvar = tvar_expr.eval(k);
981 
982  var1_expr<E> snvar_expr = var1_expr<E>(_u,_G);
983  snvar = snvar_expr.eval(k);
984 
985 
986  // VARIATION OF THE OUTER NORMAL
987  res.noalias() = tvar.cross(normal) + tangent.cross(snvar);
988 
989  return res;
990  }
991 };
992 
1000 template<class E1, class E2, class E3>
1001 class ovar2dot_expr : public _expr<ovar2dot_expr<E1,E2,E3> >
1002  {
1003 public:
1004  typedef typename E1::Scalar Scalar;
1005 
1006 private:
1007  typename E1::Nested_t _u;
1008  typename E2::Nested_t _v;
1009  typename gsGeometryMap<Scalar>::Nested_t _G;
1010  typename E3::Nested_t _C;
1011 
1012  mutable gsVector<Scalar,3> normal, onormal, tangent, utangent, dtanu, dtanv, tvaru, tvarv, tvar2,
1013  mu, mv, muv, mu_der, snvaru, snvarv, snvar2, nvar2;
1014  mutable gsMatrix<Scalar> uGrads, vGrads, cJac, res, eC;
1015 
1016  mutable gsMatrix<Scalar> tvaruMat, tvarvMat, tmp;
1017 public:
1018  enum{ Space = 3, ScalarValued= 0, ColBlocks= 0 };
1019 
1020  ovar2dot_expr(const E1 & u, const E2 & v, const gsGeometryMap<Scalar> & G, _expr<E3> const& C) : _u(u), _v(v), _G(G), _C(C) { }
1021 
1022 # define Eigen gsEigen
1023  EIGEN_MAKE_ALIGNED_OPERATOR_NEW
1024 # undef Eigen
1025 
1026  // helper function
1027  static inline gsVector<Scalar,3> vecFun(index_t pos, Scalar val)
1028  {
1029  gsVector<Scalar,3> result = gsVector<Scalar,3>::Zero();
1030  result[pos] = val;
1031  return result;
1032  }
1033 
1034  const gsMatrix<Scalar> & eval(const index_t k) const
1035  {
1036  GISMO_ASSERT(_G.data().dim.second==3,"Domain dimension should be 3, is "<<_G.data().dim.second);
1037  GISMO_ASSERT(_C.cols()*_C.rows()==3, "Size of vector is incorrect");
1038  res.resize(_u.cardinality(), _u.cardinality());
1039  res.setZero();
1040 
1041  eC = _C.eval(k);
1042  eC.resize(1,3);
1043 
1044  cJac = _G.data().values[1].reshapeCol(k, _G.data().dim.first, _G.data().dim.second).transpose();
1045  cJac.colwise().normalize();
1046 
1047  onormal = _G.data().outNormal(k);
1048  onormal.normalize();
1049  tmp = cJac.transpose() * onormal;
1050  tmp.colwise().normalize(); //normalize the inner product for fair comparison
1051 
1052  Scalar tol = 1e-8;
1053 
1054  /*
1055  We can check which column of the Jacobian corresponds to the outer normal vector or to the tangent.
1056  The tangent is a covariant vector and hence the column of the Jacobian should be equal to the tangent.
1057  The normal is a contravariant vector and hence the corresponding column of the Jacobian times the outward normal should give 1. We use this property.
1058  */
1059  index_t colIndex = -1;
1060  if ( (math::abs(tmp.at(0)) < tol) && (math::abs(tmp.at(1)) > 1-tol ) ) // then the normal is vector 2 and the tangent vector 1
1061  colIndex = 0;
1062  else if ( (math::abs(tmp.at(1)) < tol) && (math::abs(tmp.at(0)) > 1-tol ) ) // then the normal is vector 1 and the tangent vector 2
1063  colIndex = 1;
1064  else // then the normal is unknown??
1065  gsInfo<<"warning: choice unknown\n";
1066 
1067  // tangent = cJac.col(colIndex);
1068  // utangent = tangent / tangent.norm();
1069 
1070  // Required for the normal vector variation
1071  normal = _G.data().normal(k);
1072  normal.normalize();
1073  uGrads = _u.data().values[1].col(k);
1074  vGrads = _v.data().values[1].col(k);
1075 
1076  const index_t cardU = _u.data().values[0].rows(); // number of actives per component of u
1077  const index_t cardV = _v.data().values[0].rows(); // number of actives per component of v
1078  const Scalar measure = (cJac.col3d(0).cross( cJac.col3d(1) )).norm();
1079 
1080  tangent_expr<Scalar> tan_expr = tangent_expr<Scalar>(_G);
1081  tangent = tan_expr.eval(k);
1082  utangent = tangent.normalized();
1083 
1084  index_t sign = 0;
1085  if (colIndex!=-1)
1086  {
1087  Scalar dot = tangent.dot(cJac.col(colIndex));
1088  sign = (Scalar(0) < dot) - (dot < Scalar(0));
1089  }
1090  else
1091  {
1092  gsWarn<<"No suitable tangent and outer normal vector found for point "<<_G.data().values[0].transpose()<<"\n";
1093  return res;
1094  }
1095 
1096  // // For the normal vector variation
1097  // normal = _G.data().normal(k);
1098  // normal.normalize();
1099 
1100  // tvar1_expr<E1> tvaru_expr = tvar1_expr<E1>(_u,_G);
1101  // tvaruMat = tvaru_expr.eval(k);
1102  // tvar1_expr<E2> tvarv_expr = tvar1_expr<E2>(_v,_G);
1103  // tvarvMat = tvarv_expr.eval(k);
1104 
1105  for (index_t j = 0; j!= cardU; ++j) // for all basis functions u (1)
1106  {
1107  for (index_t i = 0; i!= cardV; ++i) // for all basis functions v (1)
1108  {
1109  for (index_t d = 0; d!= _u.dim(); ++d) // for all basis functions u (2)
1110  {
1111  const short_t s = d*cardU;
1112 
1113  // // first variation of the tangent
1114  // tvaru = tvaruMat.row(s+j);
1115  // first variation of the tangent (colvector)
1116  dtanu = sign*vecFun(d, uGrads.at(2*j+colIndex));
1117  tvaru = 1 / tangent.norm() * ( dtanu - ( utangent.dot(dtanu) ) * utangent );
1118 
1119  // first variation of the surface normal (colvector)
1120  mu.noalias() = ( vecFun(d, uGrads.at(2*j ) ).cross( cJac.col3d(1) )
1121  -vecFun(d, uGrads.at(2*j+1) ).cross( cJac.col3d(0) ))
1122  / measure;
1123  snvaru.noalias() = (mu - ( normal.dot(mu) ) * normal);
1124 
1125  for (index_t c = 0; c!= _v.dim(); ++c) // for all basis functions v (2)
1126  {
1127  const short_t r = c*cardV;
1128 
1129  // // first variation of the tangent
1130  // tvarv = tvarvMat.row(r+i);
1131  // first variation of the tangent (colvector)
1132  dtanv = sign*vecFun(c, vGrads.at(2*i+colIndex));
1133  tvarv = 1 / tangent.norm() * ( dtanv - ( utangent.dot(dtanv) ) * utangent );
1134 
1135  // first variation of the surface normal (colvector)
1136  mv.noalias() = ( vecFun(c, vGrads.at(2*i ) ).cross( cJac.col3d(1) )
1137  -vecFun(c, vGrads.at(2*i+1) ).cross( cJac.col3d(0) ))
1138  / measure;
1139  snvarv.noalias() = (mv - ( normal.dot(mv) ) * normal);
1140 
1141 
1142  // See point 2 of
1143  // Herrema, Austin J. et al. 2021. “Corrigendum to ‘Penalty Coupling of Non-Matching Isogeometric Kirchhoff–Love Shell Patches with Application to Composite Wind Turbine Blades’ [Comput. Methods Appl. Mech. Engrg. 346 (2019) 810–840].” Computer Methods in Applied Mechanics and Engineering 373: 113488.
1144  // Second variation of the tangent (colvector)
1145  tvar2 = -1 / tangent.norm() * (
1146  ( tvarv.dot(dtanu) * utangent )
1147  + ( utangent.dot(dtanv) * tvaru )
1148  + ( utangent.dot(dtanu) * tvarv )
1149  );
1150 
1151  // Second variation of the surface normal (colvector)
1152  muv.noalias() = ( vecFun(d, uGrads.at(2*j ) ).cross( vecFun(c, vGrads.at(2*i+1) ) )
1153  +vecFun(c, vGrads.at(2*i ) ).cross( vecFun(d, uGrads.at(2*j+1) ) ))
1154  / measure;
1155 
1156  mu_der.noalias() = (muv - ( normal.dot(mv) ) * mu);
1157 
1158  snvar2 = mu_der - (mu.dot(snvarv) + normal.dot(mu_der) ) * normal - (normal.dot(mu) ) * snvarv;
1159 
1160  // Second variation of the outer normal (colvector)
1161  nvar2 = tvar2.cross(normal) + tvaru.cross(snvarv) + tvarv.cross(snvaru) + utangent.cross(snvar2);
1162 
1163  res(s + j, r + i ) = (eC * nvar2).value();
1164  }
1165  }
1166  }
1167  }
1168  return res;
1169  }
1170 
1171  index_t rows() const { return 1; }
1172 
1173  index_t cols() const { return 1; }
1174 
1175  void parse(gsExprHelper<Scalar> & evList) const
1176  {
1177  evList.add(_u);
1178  _u.data().flags |= NEED_ACTIVE | NEED_VALUE | NEED_GRAD;
1179  evList.add(_v);
1180  _v.data().flags |= NEED_ACTIVE | NEED_VALUE | NEED_GRAD;
1181  evList.add(_G);
1182  _G.data().flags |= NEED_NORMAL | NEED_OUTER_NORMAL | NEED_DERIV;
1183  _C.parse(evList);
1184  }
1185 
1186  const gsFeSpace<Scalar> & rowVar() const { return _u.rowVar(); }
1187  const gsFeSpace<Scalar> & colVar() const { return _v.rowVar(); }
1188 
1189  void print(std::ostream &os) const { os << "nvar2("; _u.print(os); os <<")"; }
1190 };
1191 
1198 template<class E1, class E2>
1199 class deriv2dot_expr : public _expr<deriv2dot_expr<E1, E2> >
1200 {
1201  typename E1::Nested_t _u;
1202  typename E2::Nested_t _v;
1203 
1204 public:
1205  enum{ Space = (E1::Space == 1 || E2::Space == 1) ? 1 : 0,
1206  ScalarValued= 0,
1207  ColBlocks= 0
1208  };
1209 
1210  typedef typename E1::Scalar Scalar;
1211 
1212  deriv2dot_expr(const E1 & u, const E2 & v) : _u(u), _v(v) { }
1213 
1214  mutable gsMatrix<Scalar> res,tmp, vEv;
1215 
1216  const gsMatrix<Scalar> & eval(const index_t k) const {return eval_impl(_u,k); }
1217 
1218  index_t rows() const
1219  {
1220  return 1; //since the product with another vector is computed
1221  }
1222 
1223  index_t cols() const
1224  {
1225  return cols_impl(_u);
1226  }
1227 
1228  void parse(gsExprHelper<Scalar> & evList) const
1229  {
1230  parse_impl<E1>(evList);
1231  }
1232 
1233  const gsFeSpace<Scalar> & rowVar() const
1234  {
1235  if (E1::Space == 1 && E2::Space == 0)
1236  return _u.rowVar();
1237  else if (E1::Space == 0 && E2::Space == 1)
1238  return _v.rowVar();
1239  else
1240  return gsNullExpr<Scalar>::get();
1241  }
1242 
1243  const gsFeSpace<Scalar> & colVar() const
1244  {
1245  if (E1::Space == 1 && E2::Space == 0)
1246  return _v.colVar();
1247  else if (E1::Space == 0 && E2::Space == 1)
1248  return _u.colVar();
1249  else
1250  return gsNullExpr<Scalar>::get();
1251  }
1252 
1253  void print(std::ostream &os) const { os << "deriv2("; _u.print(os); _v.print(os); os <<")"; }
1254 
1255 private:
1256  template<class U> inline
1257  typename util::enable_if< util::is_same<U,gsGeometryMap<Scalar> >::value,void>::type
1258  parse_impl(gsExprHelper<Scalar> & evList) const
1259  {
1260  _u.parse(evList);
1261  evList.add(_u); // We manage the flags of _u "manually" here (sets data)
1262  _u.data().flags |= NEED_DERIV2; // define flags
1263 
1264  _v.parse(evList); // We need to evaluate _v (_v.eval(.) is called)
1265 
1266  // Note: evList.parse(.) is called only in exprAssembler for the global expression
1267  }
1268 
1269  template<class U> inline
1270  typename util::enable_if< util::is_same<U,gsFeSpace<Scalar> >::value,void>::type
1271  parse_impl(gsExprHelper<Scalar> & evList) const
1272  {
1273  _u.parse(evList);
1274  evList.add(_u); // We manage the flags of _u "manually" here (sets data)
1275  _u.data().flags |= NEED_DERIV2; // define flags
1276 
1277  _v.parse(evList); // We need to evaluate _v (_v.eval(.) is called)
1278 
1279  // Note: evList.parse(.) is called only in exprAssembler for the global expression
1280  }
1281 
1282  template<class U> inline
1283  typename util::enable_if< util::is_same<U,gsFeVariable<Scalar> >::value,void>::type
1284  parse_impl(gsExprHelper<Scalar> & evList) const
1285  {
1286  _u.parse(evList);
1287  evList.add(_u); // We manage the flags of _u "manually" here (sets data)
1288  _u.data().flags |= NEED_DERIV2; // define flags
1289 
1290  _v.parse(evList); // We need to evaluate _v (_v.eval(.) is called)
1291 
1292  // Note: evList.parse(.) is called only in exprAssembler for the global expression
1293  }
1294 
1295  template<class U> inline
1296  typename util::enable_if< util::is_same<U,gsFeSolution<Scalar> >::value,void>::type
1297  parse_impl(gsExprHelper<Scalar> & evList) const
1298  {
1299  _u.parse(evList); //
1300  hess(_u).parse(evList); //
1301 
1302  // evList.add(_u); // We manage the flags of _u "manually" here (sets data)
1303  _v.parse(evList); // We need to evaluate _v (_v.eval(.) is called)
1304  }
1305 
1306  template<class U> inline
1307  typename util::enable_if< util::is_same<U,gsGeometryMap<Scalar> >::value, const gsMatrix<Scalar> & >::type
1308  eval_impl(const U & u, const index_t k) const
1309  {
1310  /*
1311  Here, we multiply the hessian of the geometry map by a vector, which possibly has multiple actives.
1312  The hessian of the geometry map c has the form: hess(c)
1313  [d11 c1, d11 c2, d11 c3]
1314  [d22 c1, d22 c2, d22 c3]
1315  [d12 c1, d12 c2, d12 c3]
1316  And we want to compute [d11 c .v; d22 c .v; d12 c .v] ( . denotes a dot product and c and v are both vectors)
1317  So we simply evaluate for every active basis function v_k the product hess(c).v_k
1318  */
1319 
1320  // evaluate the geometry map of U
1321  tmp =_u.data().values[2].reshapeCol(k, cols(), _u.data().dim.second );
1322  vEv = _v.eval(k);
1323  res = vEv * tmp.transpose();
1324  return res;
1325  }
1326 
1327  template<class U> inline
1328  typename util::enable_if<util::is_same<U,gsFeSpace<Scalar> >::value, const gsMatrix<Scalar> & >::type
1329  eval_impl(const U & u, const index_t k) const
1330  {
1331  /*
1332  We assume that the basis has the form v*e_i where e_i is the unit vector with 1 on index i and 0 elsewhere
1333  This implies that hess(v) = [hess(v_1), hess(v_2), hess(v_3)] only has nonzero entries in column i. Hence,
1334  hess(v) . normal = hess(v_i) * n_i (vector-scalar multiplication. The result is then of the form
1335  [hess(v_1)*n_1 .., hess(v_2)*n_2 .., hess(v_3)*n_3 ..]. Here, the dots .. represent the active basis functions.
1336  */
1337  const index_t numAct = u.data().values[0].rows(); // number of actives of a basis function
1338  const index_t cardinality = u.cardinality(); // total number of actives (=3*numAct)
1339  res.resize(rows()*cardinality, cols() );
1340  tmp.transpose() =_u.data().values[2].reshapeCol(k, cols(), numAct );
1341  vEv = _v.eval(k);
1342 
1343  for (index_t i = 0; i!=_u.dim(); i++)
1344  res.block(i*numAct, 0, numAct, cols() ) = tmp * vEv.at(i);
1345  return res;
1346  }
1347 
1348  template<class U> inline
1349  typename util::enable_if<util::is_same<U,gsFeVariable<Scalar> >::value, const gsMatrix<Scalar> & >::type
1350  eval_impl(const U & u, const index_t k) const
1351  {
1352  /*
1353  Here, we multiply the hessian of the geometry map by a vector, which possibly has multiple actives.
1354  The laplacian of the variable has the form: hess(v)
1355  [d11 c1, d22 c1, d12 c1]
1356  [d11 c2, d22 c2, d12 c2]
1357  [d11 c3, d22 c3, d12 c3]
1358  And we want to compute [d11 c .v; d22 c .v; d12 c .v] ( . denotes a dot product and c and v are both vectors)
1359  So we simply evaluate for every active basis function v_k the product hess(c).v_k
1360  */
1361 
1362  gsMatrix<> tmp2;
1363  tmp = u.data().values[2].col(k);
1364  index_t nDers = _u.source().domainDim() * (_u.source().domainDim() + 1) / 2;
1365  index_t dim = _u.source().targetDim();
1366  tmp2.resize(nDers, dim);
1367  for (index_t comp = 0; comp != u.source().targetDim(); comp++)
1368  tmp2.col(comp) = tmp.block(comp * nDers, 0, nDers, 1); //star,length
1369 
1370  vEv = _v.eval(k);
1371  res = vEv * tmp2.transpose();
1372  return res;
1373  }
1374 
1375  template<class U> inline
1376  typename util::enable_if<util::is_same<U,gsFeSolution<Scalar> >::value, const gsMatrix<Scalar> & >::type
1377  eval_impl(const U & u, const index_t k) const
1378  {
1379  /*
1380  Here, we multiply the hessian of the geometry map by a vector, which possibly has multiple actives.
1381  The hessian of the geometry map c has the form: hess(c)
1382  [d11 c1, d22 c1, d12 c1]
1383  [d11 c2, d22 c2, d12 c2]
1384  [d11 c3, d22 c3, d12 c3]
1385  And we want to compute [d11 c .v; d22 c .v; d12 c .v] ( . denotes a dot product and c and v are both vectors)
1386  So we simply evaluate for every active basis function v_k the product hess(c).v_k
1387  */
1388 
1389  hess_expr<gsFeSolution<Scalar>> sHess = hess_expr<gsFeSolution<Scalar>>(_u); // NOTE: This does not parse automatically!
1390  tmp = sHess.eval(k); //.transpose();
1391  vEv = _v.eval(k);
1392  res = vEv * tmp;
1393  return res;
1394  }
1395 
1396  template<class U> inline
1397  typename util::enable_if< util::is_same<U,gsGeometryMap<Scalar> >::value ||
1398  util::is_same<U,gsFeVariable<Scalar> >::value, index_t >::type
1399  cols_impl(const U & u) const
1400  {
1401  return _u.targetDim();
1402  }
1403 
1404  template<class U> inline
1405  typename util::enable_if< util::is_same<U,gsFeSpace<Scalar> >::value ||
1406  util::is_same<U,gsFeSolution<Scalar> >::value, index_t >::type
1407  cols_impl(const U & u) const
1408  {
1409  return _u.dim();
1410  }
1411 };
1412 
1413 
1414 
1440 template<class E>
1441 class deriv2_expr : public _expr<deriv2_expr<E> >
1442 {
1443  typename E::Nested_t _u;
1444 
1445 public:
1446  // enum {ColBlocks = E::rowSpan }; // ????
1447  enum{ Space = E::Space, ScalarValued= 0, ColBlocks = (1==E::Space?1:0) };
1448 
1449  typedef typename E::Scalar Scalar;
1450 
1451  deriv2_expr(const E & u) : _u(u) { }
1452 
1453  mutable gsMatrix<Scalar> res, tmp;
1454 
1455  const gsMatrix<Scalar> & eval(const index_t k) const {return eval_impl(_u,k); }
1456 
1457  index_t rows() const //(components)
1458  {
1459  return rows_impl(_u); // _u.dim() for space or targetDim() for geometry
1460  }
1461 
1462  index_t cols() const
1463  {
1464  return cols_impl(_u);
1465  }
1466 
1467  void parse(gsExprHelper<Scalar> & evList) const
1468  {
1469  _u.parse(evList);
1470  _u.data().flags |= NEED_DERIV2;
1471  }
1472 
1473  const gsFeSpace<Scalar> & rowVar() const { return _u.rowVar(); }
1474  const gsFeSpace<Scalar> & colVar() const { return _u.colVar(); }
1475 
1476  index_t cardinality_impl() const { return _u.cardinality_impl(); }
1477 
1478  void print(std::ostream &os) const { os << "deriv2("; _u.print(os); os <<")"; }
1479 
1480  private:
1481  template<class U> inline
1482  typename util::enable_if< util::is_same<U,gsGeometryMap<Scalar> >::value, const gsMatrix<Scalar> & >::type
1483  eval_impl(const U & u, const index_t k) const
1484  {
1485  /*
1486  Here, we compute the hessian of the geometry map.
1487  The hessian of the geometry map c has the form: hess(c)
1488  [d11 c1, d11 c2, d11 c3]
1489  [d22 c1, d22 c2, d22 c3]
1490  [d12 c1, d12 c2, d12 c3]
1491 
1492  The geometry map has components c=[c1,c2,c3]
1493  */
1494  // evaluate the geometry map of U
1495  res = _u.data().values[2].reshapeCol(k, cols(), _u.data().dim.second );
1496  return res;
1497  }
1498 
1499  template<class U> inline
1500  typename util::enable_if< util::is_same<U,gsFeVariable<Scalar> >::value, const gsMatrix<Scalar> & >::type
1501  eval_impl(const U & u, const index_t k) const
1502  {
1503  /*
1504  Here, we compute the hessian of the geometry map.
1505  The hessian of the geometry map c has the form: hess(c)
1506  [d11 c1, d11 c2, d11 c3]
1507  [d22 c1, d22 c2, d22 c3]
1508  [d12 c1, d12 c2, d12 c3]
1509 
1510  The geometry map has components c=[c1,c2,c3]
1511  */
1512  // evaluate the geometry map of U
1513  tmp = _u.data().values[2];
1514  res.resize(rows(),cols());
1515  for (index_t comp = 0; comp != _u.source().targetDim(); comp++)
1516  res.col(comp) = tmp.block(comp*rows(),0,rows(),1); //star,length
1517  return res;
1518  }
1519 
1520  template<class U> inline
1521  typename util::enable_if< util::is_same<U,gsFeSolution<Scalar> >::value, const gsMatrix<Scalar> & >::type
1522  eval_impl(const U & u, const index_t k) const
1523  {
1524  /*
1525  Here, we compute the hessian of the geometry map.
1526  The hessian of the geometry map c has the form: hess(c)
1527  [d11 c1, d11 c2, d11 c3]
1528  [d22 c1, d22 c2, d22 c3]
1529  [d12 c1, d12 c2, d12 c3]
1530 
1531  The geometry map has components c=[c1,c2,c3]
1532  */
1533  // evaluate the geometry map of U
1534  hess_expr<gsFeSolution<Scalar>> sHess = hess_expr<gsFeSolution<Scalar>>(_u);
1535  res = sHess.eval(k).transpose();
1536  return res;
1537  }
1538 
1540  template<class U> inline
1541  typename util::enable_if<util::is_same<U,gsFeSpace<Scalar> >::value, const gsMatrix<Scalar> & >::type
1542  eval_impl(const U & u, const index_t k) const
1543  {
1544  /*
1545  Here, we compute the hessian of the basis with n actives.
1546  The hessian of the basis u has the form: hess(u)
1547  active 1 active 2 active n = cardinality
1548  [d11 u1, d11 u2, d11 u3] [d11 u1, d11 u2, d11 u3] ... [d11 u1, d11 u2, d11 u3]
1549  [d22 u1, d22 u2, d22 u3] [d22 u1, d22 u2, d22 u3] ... [d22 u1, d22 u2, d22 u3]
1550  [d12 u1, d12 u2, d12 u3] [d12 u1, d12 u2, d12 u3] ... [d12 u1, d12 u2, d12 u3]
1551 
1552  Here, the basis function has components u = [u1,u2,u3]. Since they are evaluated for scalars
1553  we use blockDiag to make copies for all components ui
1554 
1555  active 1 active 2 active k = cardinality/dim active 1 active 2k active 1 active 2k
1556  [d11 u, 0, 0] [d11 u, 0, 0] ... [d11 u, 0, 0] [0, d11 u, 0] ... [0, d11 u, 0] [0, d11 u, 0] ... [0, d11 u, 0]
1557  [d22 u, 0, 0] [d22 u, 0, 0] ... [d22 u, 0, 0] [0, d22 u, 0] ... [0, d22 u, 0] [0, d22 u, 0] ... [0, d22 u, 0]
1558  [d12 u, 0, 0] [d12 u, 0, 0] ... [d12 u, 0, 0] [0, d12 u, 0] ... [0, d12 u, 0] [0, d12 u, 0] ... [0, d12 u, 0]
1559 
1560  */
1561  const index_t numAct = u.data().values[0].rows(); // number of actives of a basis function
1562  const index_t cardinality = u.cardinality(); // total number of actives (=3*numAct)
1563 
1564  res.resize(rows(), _u.dim() *_u.cardinality()); // (3 x 3*cardinality)
1565  res.setZero();
1566 
1567  tmp = _u.data().values[2].reshapeCol(k, cols(), numAct );
1568  for (index_t d = 0; d != cols(); ++d)
1569  {
1570  const index_t s = d*(cardinality + 1);
1571  for (index_t i = 0; i != numAct; ++i)
1572  res.col(s+i*_u.cols()) = tmp.col(i);
1573  }
1574 
1575  return res;
1576  }
1577 
1578  template<class U> inline
1579  typename util::enable_if< util::is_same<U,gsFeVariable<Scalar> >::value || util::is_same<U,gsGeometryMap<Scalar> >::value || util::is_same<U,gsFeSpace<Scalar> >::value, index_t >::type
1580  rows_impl(const U & u) const
1581  {
1582  return _u.source().domainDim() * ( _u.source().domainDim() + 1 ) / 2;
1583  }
1584 
1585  template<class U> inline
1586  typename util::enable_if< util::is_same<U,gsFeSolution<Scalar> >::value, index_t >::type
1587  rows_impl(const U & u) const
1588  {
1589  return _u.parDim() * ( _u.parDim() + 1 ) / 2;
1590  }
1591 
1593  template<class U> inline
1594  typename util::enable_if< util::is_same<U,gsFeVariable<Scalar> >::value || util::is_same<U,gsGeometryMap<Scalar> >::value, index_t >::type
1595  cols_impl(const U & u) const
1596  {
1597  return u.source().targetDim();
1598  }
1599 
1600  template<class U> inline
1601  typename util::enable_if< util::is_same<U,gsFeSolution<Scalar> >::value || util::is_same<U,gsFeSpace<Scalar> >::value, index_t >::type
1602  cols_impl(const U & u) const
1603  {
1604  return _u.dim();
1605  }
1606 
1607 };
1608 
1616 template<class E1, class E2, class E3>
1617 class flatdot_expr : public _expr<flatdot_expr<E1,E2,E3> >
1618 {
1619 public:
1620  typedef typename E1::Scalar Scalar;
1621 
1622 private:
1623  typename E1::Nested_t _A;
1624  typename E2::Nested_t _B;
1625  typename E3::Nested_t _C;
1626  mutable gsMatrix<Scalar> eA, eB, eC, tmp, res;
1627 
1628 public:
1629  enum {Space = 3, ScalarValued = 0, ColBlocks = 0};
1630 
1631 public:
1632 
1633  flatdot_expr(_expr<E1> const& A, _expr<E2> const& B, _expr<E3> const& C) : _A(A),_B(B),_C(C)
1634  {
1635  }
1636 
1637  const gsMatrix<Scalar> & eval(const index_t k) const
1638  {
1639  const index_t Ac = _A.cols();
1640  const index_t An = _A.cardinality();
1641  const index_t Bc = _B.cols();
1642  const index_t Bn = _B.cardinality();
1643 
1644  eA = _A.eval(k);
1645  eB = _B.eval(k);
1646  eC = _C.eval(k);
1647 
1648  GISMO_ASSERT(Bc==_A.rows(), "Dimensions: "<<Bc<<","<< _A.rows()<< "do not match");
1649  GISMO_STATIC_ASSERT(E1::Space==1, "First entry should be rowSpan" );
1650  GISMO_STATIC_ASSERT(E2::Space==2, "Second entry should be colSpan.");
1651 
1652  res.resize(An, Bn);
1653 
1654  // to do: evaluate the jacobians internally for the basis functions and loop over dimensions as well, since everything is same for all dimensions.
1655  for (index_t i = 0; i!=An; ++i) // for all actives u
1656  for (index_t j = 0; j!=Bn; ++j) // for all actives v
1657  {
1658  tmp.noalias() = eB.middleCols(i*Bc,Bc) * eA.middleCols(j*Ac,Ac);
1659 
1660  tmp(0,0) *= eC.at(0);
1661  tmp(0,1) *= eC.at(2);
1662  tmp(1,0) *= eC.at(2);
1663  tmp(1,1) *= eC.at(1);
1664  res(i,j) = tmp.sum();
1665  }
1666  return res;
1667  }
1668 
1669  index_t rows() const { return 1; }
1670  index_t cols() const { return 1; }
1671 
1672  void parse(gsExprHelper<Scalar> & evList) const
1673  { _A.parse(evList);_B.parse(evList);_C.parse(evList); }
1674 
1675  const gsFeSpace<Scalar> & rowVar() const { return _A.rowVar(); }
1676  const gsFeSpace<Scalar> & colVar() const { return _B.colVar(); }
1677  index_t cardinality_impl() const { return _A.cardinality_impl(); }
1678 
1679  void print(std::ostream &os) const { os << "flatdot("; _A.print(os);_B.print(os);_C.print(os); os<<")"; }
1680 };
1681 
1691 template<class E1, class E2, class E3>
1692 class flatdot2_expr : public _expr<flatdot2_expr<E1,E2,E3> >
1693 {
1694 public:
1695  typedef typename E1::Scalar Scalar;
1696 
1697 private:
1698  typename E1::Nested_t _A;
1699  typename E2::Nested_t _B;
1700  typename E3::Nested_t _C;
1701  mutable gsMatrix<Scalar> eA, eB, eC, res, tmp;
1702 
1703 public:
1704  enum {Space = E1::Space, ScalarValued = 0, ColBlocks = 0};
1705 
1706  flatdot2_expr(_expr<E1> const& A, _expr<E2> const& B, _expr<E3> const& C) : _A(A),_B(B),_C(C)
1707  {
1708  }
1709 
1710  const gsMatrix<Scalar> & eval(const index_t k) const
1711  {
1712  const index_t Ac = _A.cols();
1713  const index_t An = _A.cardinality();
1714  const index_t Bn = _B.cardinality();
1715 
1716  eA = _A.eval(k);
1717  eB = _B.eval(k);
1718  eC = _C.eval(k);
1719 
1720  GISMO_ASSERT(_B.rows()==_A.cols(), "Dimensions: "<<_B.rows()<<","<< _A.cols()<< "do not match");
1721  GISMO_STATIC_ASSERT(E1::Space==1, "First entry should be rowSpan");
1722  GISMO_STATIC_ASSERT(E2::Space==2, "Second entry should be colSpan.");
1723  GISMO_ASSERT(_C.cols()==_B.rows(), "Dimensions: "<<_C.rows()<<","<< _B.rows()<< "do not match");
1724 
1725  res.resize(An, Bn);
1726  for (index_t i = 0; i!=An; ++i)
1727  for (index_t j = 0; j!=Bn; ++j)
1728  {
1729  tmp = eA.middleCols(i*Ac,Ac) * eB.col(j); // E_f_der2
1730  tmp.row(2) *= 2.0; // multiply the third row of E_f_der2 by 2 for voight notation
1731  res(i,j) = (eC.row(0) * tmp.col(0)).value(); // E_f^T * mm * E_f_der2
1732  }
1733  return res;
1734  }
1735 
1736  index_t rows() const { return 1; }
1737  index_t cols() const { return 1; }
1738 
1739  void parse(gsExprHelper<Scalar> & evList) const
1740  { _A.parse(evList);_B.parse(evList);_C.parse(evList); }
1741 
1742  const gsFeSpace<Scalar> & rowVar() const { return _A.rowVar(); }
1743  const gsFeSpace<Scalar> & colVar() const { return _B.colVar(); }
1744  index_t cardinality_impl() const { return _A.cardinality_impl(); }
1745 
1746  void print(std::ostream &os) const { os << "flatdot2("; _A.print(os);_B.print(os);_C.print(os); os<<")"; }
1747 };
1748 
1753 template<class T> class cartcovinv_expr ;
1754 
1755 template<class T>
1756 class cartcov_expr : public _expr<cartcov_expr<T> >
1757 {
1758  typename gsGeometryMap<T>::Nested_t _G;
1759 
1760 public:
1761  typedef T Scalar;
1762 
1763  enum {Space = 0, ScalarValued = 0, ColBlocks = 0};
1764 
1765  cartcov_expr(const gsGeometryMap<T> & G) : _G(G) { }
1766 
1767  mutable gsMatrix<Scalar> covBasis;
1768  mutable gsMatrix<Scalar,3,3> result;
1769  mutable gsVector<Scalar> normal, tmp;
1770  mutable gsVector<Scalar> e1, e2, a1, a2;
1771 
1772  gsMatrix<Scalar> eval(const index_t k) const
1773  {
1774  if (_G.targetDim()==3)
1775  {
1776  // Compute covariant bases in deformed and undeformed configuration
1777  normal = _G.data().normals.col(k);
1778  normal.normalize();
1779  covBasis.resize(3,3);
1780  covBasis.leftCols(2) = _G.data().jacobian(k);
1781  covBasis.col(2) = normal;
1782  }
1783  else if (_G.targetDim()==2)
1784  {
1785  // Compute covariant bases in deformed and undeformed configuration
1786  covBasis.resize(2,2);
1787  covBasis = _G.data().jacobian(k);
1788  }
1789  else
1790  GISMO_ERROR("Not implemented");
1791 
1792  a1 = covBasis.col(0);
1793  a2 = covBasis.col(1);
1794 
1795  e1 = a1.normalized();
1796  e2 = (a2-(a2.dot(e1)*e1)).normalized();
1797 
1798  result(0,0) = (e1.dot(a1))*(a1.dot(e1)); // 1111
1799  result(0,1) = (e1.dot(a2))*(a2.dot(e1)); // 1122
1800  result(0,2) = 2*(e1.dot(a1))*(a2.dot(e1)); // 1112
1801  // Row 1
1802  result(1,0) = (e2.dot(a1))*(a1.dot(e2)); // 2211
1803  result(1,1) = (e2.dot(a2))*(a2.dot(e2)); // 2222
1804  result(1,2) = 2*(e2.dot(a1))*(a2.dot(e2)); // 2212
1805  // Row 2
1806  result(2,0) = (e1.dot(a1))*(a1.dot(e2)); // 1211
1807  result(2,1) = (e1.dot(a2))*(a2.dot(e2)); // 1222
1808  result(2,2) = (e1.dot(a1))*(a2.dot(e2)) + (e1.dot(a2))*(a1.dot(e2)); // 1212
1809  return result;
1810  }
1811 
1812  cartcovinv_expr<T> inv() const
1813  {
1814  GISMO_ASSERT(rows() == cols(), "The Jacobian matrix is not square");
1815  return cartcovinv_expr<T>(_G);
1816  }
1817 
1818  index_t rows() const { return 3; }
1819 
1820  index_t cols() const { return 3; }
1821 
1822  static const gsFeSpace<Scalar> & rowVar() { return gsNullExpr<Scalar>::get(); }
1823  static const gsFeSpace<Scalar> & colVar() { return gsNullExpr<Scalar>::get(); }
1824 
1825  void parse(gsExprHelper<Scalar> & evList) const
1826  {
1827  //GISMO_ASSERT(NULL!=m_fd, "FeVariable: FuncData member not registered");
1828  evList.add(_G);
1829  _G.data().flags |= NEED_NORMAL|NEED_DERIV;
1830  }
1831 
1832  void print(std::ostream &os) const { os << "cartcov("; _G.print(os); os <<")"; }
1833 };
1834 
1835 template<class T>
1836 class cartcovinv_expr : public _expr<cartcovinv_expr<T> >
1837 {
1838  typename gsGeometryMap<T>::Nested_t _G;
1839 
1840 public:
1841  typedef T Scalar;
1842 
1843  enum {Space = 0, ScalarValued = 0, ColBlocks = 0};
1844 
1845  cartcovinv_expr(const gsGeometryMap<T> & G) : _G(G) { }
1846 
1847  mutable gsMatrix<T> temp;
1848 
1849  gsMatrix<T> eval(const index_t k) const
1850  {
1851  cartcov_expr<Scalar> cartcov = cartcov_expr<Scalar>(_G);
1852  temp = (cartcov.eval(k)).reshape(3,3).inverse();
1853  return temp;
1854  }
1855 
1856  index_t rows() const { return 3; }
1857 
1858  index_t cols() const { return 3; }
1859 
1860  static const gsFeSpace<Scalar> & rowVar() { return gsNullExpr<Scalar>::get(); }
1861  static const gsFeSpace<Scalar> & colVar() { return gsNullExpr<Scalar>::get(); }
1862 
1863  void parse(gsExprHelper<Scalar> & evList) const
1864  {
1865  //GISMO_ASSERT(NULL!=m_fd, "FeVariable: FuncData member not registered");
1866  cartcov(_G).parse(evList); //
1867 
1868  evList.add(_G);
1869  _G.data().flags |= NEED_NORMAL|NEED_DERIV;
1870  }
1871 
1872  void print(std::ostream &os) const { os << "cartcovinv("; _G.print(os); os <<")"; }
1873 };
1874 
1878 template<class T> class cartconinv_expr ;
1879 
1880 template<class T>
1881 class cartcon_expr : public _expr<cartcon_expr<T> >
1882 {
1883  typename gsGeometryMap<T>::Nested_t _G;
1884 
1885 public:
1886  typedef T Scalar;
1887 
1888  enum {Space = 0, ScalarValued = 0, ColBlocks = 0};
1889 
1890  cartcon_expr(const gsGeometryMap<T> & G) : _G(G) { }
1891 
1892  mutable gsMatrix<Scalar> covBasis, conBasis, covMetric, conMetric, cartBasis;
1893  mutable gsMatrix<Scalar,3,3> result;
1894  mutable gsVector<Scalar> normal, tmp;
1895  mutable gsVector<Scalar> e1, e2, ac1, ac2, a1, a2;
1896 
1897  gsMatrix<Scalar> eval(const index_t k) const
1898  {
1899  if (_G.targetDim()==3)
1900  {
1901  // Compute covariant bases in deformed and undeformed configuration
1902  normal = _G.data().normals.col(k);
1903  normal.normalize();
1904 
1905  covBasis.resize(3,3);
1906  conBasis.resize(3,3);
1907  covBasis.leftCols(2) = _G.data().jacobian(k);
1908  covBasis.col(2) = normal;
1909  covMetric = covBasis.transpose() * covBasis;
1910 
1911  conMetric = covMetric.inverse();
1912 
1913  conBasis.col(0) = conMetric(0,0)*covBasis.col(0)+conMetric(0,1)*covBasis.col(1)+conMetric(0,2)*covBasis.col(2);
1914  conBasis.col(1) = conMetric(1,0)*covBasis.col(0)+conMetric(1,1)*covBasis.col(1)+conMetric(1,2)*covBasis.col(2);
1915  }
1916  else if (_G.targetDim()==2)
1917  {
1918  // Compute covariant bases in deformed and undeformed configuration
1919  normal.resize(3);
1920  normal<<0,0,1;
1921  normal.normalize();
1922  covBasis.resize(2,2);
1923  conBasis.resize(2,2);
1924  // Compute covariant bases in deformed and undeformed configuration
1925  covBasis = _G.data().jacobian(k);
1926  covMetric = covBasis.transpose() * covBasis;
1927 
1928  conMetric = covMetric.inverse();
1929 
1930  conBasis.col(0) = conMetric(0,0)*covBasis.col(0)+conMetric(0,1)*covBasis.col(1);
1931  conBasis.col(1) = conMetric(1,0)*covBasis.col(0)+conMetric(1,1)*covBasis.col(1);
1932  }
1933  else
1934  GISMO_ERROR("Not implemented");
1935 
1936  a1 = covBasis.col(0);
1937  a2 = covBasis.col(1);
1938 
1939  e1 = a1.normalized();
1940  e2 = (a2-(a2.dot(e1)*e1)).normalized();
1941 
1942  ac1 = conBasis.col(0);
1943  ac2 = conBasis.col(1);
1944 
1945  result(0,0) = (e1.dot(ac1))*(ac1.dot(e1)); // 1111
1946  result(0,1) = (e1.dot(ac2))*(ac2.dot(e1)); // 1122
1947  result(0,2) = 2*(e1.dot(ac1))*(ac2.dot(e1)); // 1112
1948  // Row 1
1949  result(1,0) = (e2.dot(ac1))*(ac1.dot(e2)); // 2211
1950  result(1,1) = (e2.dot(ac2))*(ac2.dot(e2)); // 2222
1951  result(1,2) = 2*(e2.dot(ac1))*(ac2.dot(e2)); // 2212
1952  // Row 2
1953  result(2,0) = (e1.dot(ac1))*(ac1.dot(e2)); // 1211
1954  result(2,1) = (e1.dot(ac2))*(ac2.dot(e2)); // 1222
1955  result(2,2) = (e1.dot(ac1))*(ac2.dot(e2)) + (e1.dot(ac2))*(ac1.dot(e2)); // 1212
1956  return result;
1957 
1958  }
1959 
1960  cartconinv_expr<T> inv() const
1961  {
1962  GISMO_ASSERT(rows() == cols(), "The Jacobian matrix is not square");
1963  return cartconinv_expr<T>(_G);
1964  }
1965 
1966  index_t rows() const { return 3; }
1967 
1968  index_t cols() const { return 3; }
1969 
1970  static const gsFeSpace<Scalar> & rowVar() { return gsNullExpr<Scalar>::get(); }
1971  static const gsFeSpace<Scalar> & colVar() { return gsNullExpr<Scalar>::get(); }
1972 
1973  void parse(gsExprHelper<Scalar> & evList) const
1974  {
1975  //GISMO_ASSERT(NULL!=m_fd, "FeVariable: FuncData member not registered");
1976  evList.add(_G);
1977  _G.data().flags |= NEED_NORMAL|NEED_DERIV;
1978  }
1979 
1980  void print(std::ostream &os) const { os << "cartcon("; _G.print(os); os <<")"; }
1981 };
1982 
1983 template<class T>
1984 class cartconinv_expr : public _expr<cartconinv_expr<T> >
1985 {
1986 private:
1987  typename gsGeometryMap<T>::Nested_t _G;
1988  mutable gsMatrix<T> temp;
1989 
1990 public:
1991  typedef T Scalar;
1992 
1993  enum {Space = 0, ScalarValued = 0, ColBlocks = 0};
1994 
1995  cartconinv_expr(const gsGeometryMap<T> & G) : _G(G) { }
1996 
1997  gsMatrix<T> eval(const index_t k) const
1998  {
1999  cartcon_expr<Scalar> cartcon = cartcon_expr<Scalar>(_G);
2000  temp = (cartcon.eval(k)).reshape(3,3).inverse();
2001  return temp;
2002  }
2003 
2004  index_t rows() const { return 3; }
2005 
2006  index_t cols() const { return 3; }
2007 
2008  static const gsFeSpace<Scalar> & rowVar() { return gsNullExpr<Scalar>::get(); }
2009  static const gsFeSpace<Scalar> & colVar() { return gsNullExpr<Scalar>::get(); }
2010 
2011  void parse(gsExprHelper<Scalar> & evList) const
2012  {
2013  cartcon(_G).parse(evList); //
2014 
2015  //GISMO_ASSERT(NULL!=m_fd, "FeVariable: FuncData member not registered");
2016  evList.add(_G);
2017  _G.data().flags |= NEED_NORMAL|NEED_DERIV;
2018  }
2019 
2020  void print(std::ostream &os) const { os << "cartconinv("; _G.print(os); os <<")"; }
2021 };
2022 
2023 EIGEN_STRONG_INLINE
2024 unitVec_expr uv(const index_t index, const index_t dim) { return unitVec_expr(index,dim); }
2025 
2026 template<class E> EIGEN_STRONG_INLINE
2027 var1_expr<E> var1(const E & u, const gsGeometryMap<typename E::Scalar> & G) { return var1_expr<E>(u, G); }
2028 
2029 template<class E1, class E2> EIGEN_STRONG_INLINE
2030 var2_expr<E1,E2> var2(const E1 & u, const E2 & v, const gsGeometryMap<typename E1::Scalar> & G)
2031 { return var2_expr<E1,E2>(u,v, G); }
2032 
2033 template<class E1, class E2, class E3> EIGEN_STRONG_INLINE
2034 var2dot_expr<E1,E2,E3> var2(const E1 & u, const E2 & v, const gsGeometryMap<typename E1::Scalar> & G, const E3 & Ef)
2035 { return var2dot_expr<E1,E2,E3>(u,v, G, Ef); }
2036 
2037 template<class E1, class E2, class E3> EIGEN_STRONG_INLINE
2038 var2deriv2dot_expr<E1,E2,E3> var2deriv2(const E1 & u, const E2 & v, const gsGeometryMap<typename E1::Scalar> & G, const E3 & Ef)
2039 { return var2deriv2dot_expr<E1,E2,E3>(u,v, G, Ef); }
2040 
2041 template<class E> EIGEN_STRONG_INLINE
2042 tvar1_expr<E> tvar1(const E & u, const gsGeometryMap<typename E::Scalar> & G) { return tvar1_expr<E>(u, G); }
2043 
2044 template<class E> EIGEN_STRONG_INLINE
2045 ovar1_expr<E> ovar1(const E & u, const gsGeometryMap<typename E::Scalar> & G) { return ovar1_expr<E>(u, G); }
2046 
2047 template<class E1, class E2, class E3> EIGEN_STRONG_INLINE
2048 ovar2dot_expr<E1,E2,E3> ovar2(const E1 & u, const E2 & v, const gsGeometryMap<typename E1::Scalar> & G, const E3 & C)
2049 { return ovar2dot_expr<E1,E2,E3>(u,v, G, C); }
2050 
2051 // template<class E1, class E2> EIGEN_STRONG_INLINE
2052 // hessdot_expr<E1,E2> hessdot(const E1 & u, const E2 & v) { return hessdot_expr<E1,E2>(u, v); }
2053 
2054 template<class E> EIGEN_STRONG_INLINE
2055 deriv2_expr<E> deriv2(const E & u) { return deriv2_expr<E>(u); }
2056 
2057 template<class E1, class E2> EIGEN_STRONG_INLINE
2058 deriv2dot_expr<E1, E2> deriv2(const E1 & u, const E2 & v) { return deriv2dot_expr<E1, E2>(u,v); }
2059 
2060 template<class E1, class E2, class E3> EIGEN_STRONG_INLINE
2061 flatdot_expr<E1,E2,E3> flatdot(const E1 & u, const E2 & v, const E3 & w)
2062 { return flatdot_expr<E1,E2,E3>(u, v, w); }
2063 
2064 template<class E1, class E2, class E3> EIGEN_STRONG_INLINE
2065 flatdot2_expr<E1,E2,E3> flatdot2(const E1 & u, const E2 & v, const E3 & w)
2066 { return flatdot2_expr<E1,E2,E3>(u, v, w); }
2067 
2068 template<class E> EIGEN_STRONG_INLINE
2069 cartcov_expr<E> cartcov(const gsGeometryMap<E> & G) { return cartcov_expr<E>(G); }
2070 
2071 template<class E> EIGEN_STRONG_INLINE
2072 cartcon_expr<E> cartcon(const gsGeometryMap<E> & G) { return cartcon_expr<E>(G); }
2073 
2074 }
2075 }
gismo::expr::deriv2_expr
Computes the second derivative of an expression.
Definition: gsThinShellUtils.h:1441
GISMO_NO_IMPLEMENTATION
#define GISMO_NO_IMPLEMENTATION
Definition: gsDebug.h:129
gismo::gsExprHelper
Definition: gsExpressions.h:96
gismo::NEED_DERIV
Gradient of the object.
Definition: gsForwardDeclarations.h:73
short_t
#define short_t
Definition: gsConfig.h:35
gismo::gsMatrix::col3d
Col3DType col3d(index_t c)
Returns column c as a fixed-size 3D vector.
Definition: gsMatrix.h:244
gsExpressions.h
Defines different expressions.
gismo::expr::unitVec_expr
Simple expression for the unit vector of length dim and with value 1 on index.
Definition: gsThinShellUtils.h:40
gismo::NEED_OUTER_NORMAL
Outward normal on the boundary.
Definition: gsForwardDeclarations.h:86
gismo::expr::var2deriv2dot_expr
Second variation of the surface normal times the second derivative of the geometry map times a vector...
Definition: gsThinShellUtils.h:375
gismo::expr::tv
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
gismo::expr::var1_expr
Expression for the first variation of the surface normal.
Definition: gsThinShellUtils.h:80
gsBasis.h
Provides declaration of Basis abstract interface.
gismo::NEED_DERIV2
Second derivatives.
Definition: gsForwardDeclarations.h:80
index_t
#define index_t
Definition: gsConfig.h:32
gsAssemblerOptions.h
Provides assembler and solver options.
gismo::expr::flatdot_expr
Computes the product of expressions E1 and E2 and multiplies with a vector E3 in voight notation...
Definition: gsThinShellUtils.h:1617
GISMO_ASSERT
#define GISMO_ASSERT(cond, message)
Definition: gsDebug.h:89
gsDofMapper.h
Provides the gsDofMapper class for re-indexing DoFs.
gismo::gsMatrix::at
T at(index_t i) const
Returns the i-th element of the vectorization of the matrix.
Definition: gsMatrix.h:211
gsWarn
#define gsWarn
Definition: gsDebug.h:50
gismo::NEED_NORMAL
Normal vector of the object.
Definition: gsForwardDeclarations.h:85
gismo::expr::var2_expr
Second variation of the normal.
Definition: gsThinShellUtils.h:499
gismo::expr::deriv2_expr::eval_impl
util::enable_if< util::is_same< U, gsFeSpace< Scalar > >::value, const gsMatrix< Scalar > & >::type eval_impl(const U &u, const index_t k) const
Spexialization for a space.
Definition: gsThinShellUtils.h:1542
gismo::expr::ovar1_expr
Expression for the first variation of the outer normal.
Definition: gsThinShellUtils.h:853
gismo::expr::deriv2dot_expr
Expression that takes the second derivative of an expression and multiplies it with a row vector...
Definition: gsThinShellUtils.h:1199
gsInfo
#define gsInfo
Definition: gsDebug.h:43
gismo::NEED_ACTIVE
Active function ids.
Definition: gsForwardDeclarations.h:84
gismo::expr::flatdot2_expr
Computes the product of expressions E1 and E2 and multiplies with a vector E3 in voight notation...
Definition: gsThinShellUtils.h:1692
gismo::gsVector::at
T at(index_t i) const
Returns the i-th element of the vector.
Definition: gsVector.h:177
gismo::expr::ovar2dot_expr
Expression for the second variation of the outer normal times a vector.
Definition: gsThinShellUtils.h:1001
gismo::expr::reshape
EIGEN_STRONG_INLINE reshape_expr< E > const reshape(E const &u, index_t n, index_t m)
Reshape an expression.
Definition: gsExpressions.h:1925
gismo::expr::cartconinv_expr
Definition: gsThinShellUtils.h:1878
gismo::NEED_VALUE
Value of the object.
Definition: gsForwardDeclarations.h:72
gismo::expr::grad
EIGEN_STRONG_INLINE grad_expr< E > grad(const E &u)
The gradient of a variable.
Definition: gsExpressions.h:4490
gsPointGrid.h
Provides functions to generate structured point data.
gsBoundaryConditions.h
Provides gsBoundaryConditions class.
gismo::expr::cartcovinv_expr
Definition: gsThinShellUtils.h:1753
GISMO_ERROR
#define GISMO_ERROR(message)
Definition: gsDebug.h:118
gsLinearAlgebra.h
This is the main header file that collects wrappers of Eigen for linear algebra.
gismo::expr::abs
EIGEN_STRONG_INLINE abs_expr< E > abs(const E &u)
Absolute value.
Definition: gsExpressions.h:4486
gismo::expr::var2dot_expr
Second variation of the surface normal times a vector.
Definition: gsThinShellUtils.h:248
gismo::expr::tvar1_expr
Expression for the first variation of the outer tangent.
Definition: gsThinShellUtils.h:659
gsFuncData.h
This object is a cache for computed values from an evaluator.
gismo::expr::tangent_expr
Definition: gsExpressions.h:3080