G+Smo  25.01.0
Geometry + Simulation Modules
 
Loading...
Searching...
No Matches
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
34namespace gismo{
35namespace expr{
36
40class unitVec_expr : public _expr<unitVec_expr >
41{
42public:
43 typedef real_t Scalar;
44private:
45 index_t _index;
46 index_t _dim;
47
48public:
49 unitVec_expr(const index_t index, const index_t dim) : _index(index), _dim(dim) { }
50
51public:
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
79template<class E>
80class var1_expr : public _expr<var1_expr<E> >
81{
82public:
83 typedef typename E::Scalar Scalar;
84
85private:
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;
92public:
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
125private:
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
247template<class E1, class E2, class E3>
248class var2dot_expr : public _expr<var2dot_expr<E1,E2,E3> >
249{
250public:
251 typedef typename E1::Scalar Scalar;
252
253private:
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
259public:
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
374template<class E1, class E2, class E3>
375class var2deriv2dot_expr : public _expr<var2deriv2dot_expr<E1,E2,E3> >
376{
377public:
378 typedef typename E1::Scalar Scalar;
379
380private:
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
386public:
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
498template<class E1, class E2>
499class var2_expr : public _expr<var2_expr<E1,E2> >
500{
501public:
502 typedef typename E1::Scalar Scalar;
503
504private:
505
506 typename E1::Nested_t _u;
507 typename E2::Nested_t _v;
508 typename gsGeometryMap<Scalar>::Nested_t _G;
509
510public:
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
569private:
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
658template<class E>
659class tvar1_expr : public _expr<tvar1_expr<E> >
660{
661public:
662 typedef typename E::Scalar Scalar;
663
664private:
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
672public:
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
704private:
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
852template<class E>
853class ovar1_expr : public _expr<ovar1_expr<E> >
854{
855public:
856 typedef typename E::Scalar Scalar;
857
858private:
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
866public:
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
899private:
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
1000template<class E1, class E2, class E3>
1001class ovar2dot_expr : public _expr<ovar2dot_expr<E1,E2,E3> >
1002 {
1003public:
1004 typedef typename E1::Scalar Scalar;
1005
1006private:
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;
1017public:
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
1198template<class E1, class E2>
1199class deriv2dot_expr : public _expr<deriv2dot_expr<E1, E2> >
1200{
1201 typename E1::Nested_t _u;
1202 typename E2::Nested_t _v;
1203
1204public:
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
1255private:
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
1440template<class E>
1441class deriv2_expr : public _expr<deriv2_expr<E> >
1442{
1443 typename E::Nested_t _u;
1444
1445public:
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 GISMO_UNUSED(k); //should this use k??
1504 /*
1505 Here, we compute the hessian of the geometry map.
1506 The hessian of the geometry map c has the form: hess(c)
1507 [d11 c1, d11 c2, d11 c3]
1508 [d22 c1, d22 c2, d22 c3]
1509 [d12 c1, d12 c2, d12 c3]
1510
1511 The geometry map has components c=[c1,c2,c3]
1512 */
1513 // evaluate the geometry map of U
1514 tmp = u.data().values[2];
1515 res.resize(rows(),cols());
1516 for (index_t comp = 0; comp != u.source().targetDim(); comp++)
1517 res.col(comp) = tmp.block(comp*rows(),0,rows(),1); //star,length
1518 return res;
1519 }
1520
1521 template<class U> inline
1522 typename util::enable_if< util::is_same<U,gsFeSolution<Scalar> >::value, const gsMatrix<Scalar> & >::type
1523 eval_impl(const U & u, const index_t k) const
1524 {
1525 /*
1526 Here, we compute the hessian of the geometry map.
1527 The hessian of the geometry map c has the form: hess(c)
1528 [d11 c1, d11 c2, d11 c3]
1529 [d22 c1, d22 c2, d22 c3]
1530 [d12 c1, d12 c2, d12 c3]
1531
1532 The geometry map has components c=[c1,c2,c3]
1533 */
1534 // evaluate the geometry map of U
1535 hess_expr<gsFeSolution<Scalar>> sHess = hess_expr<gsFeSolution<Scalar>>(_u);
1536 res = sHess.eval(k).transpose();
1537 return res;
1538 }
1539
1541 template<class U> inline
1542 typename util::enable_if<util::is_same<U,gsFeSpace<Scalar> >::value, const gsMatrix<Scalar> & >::type
1543 eval_impl(const U & u, const index_t k) const
1544 {
1545 /*
1546 Here, we compute the hessian of the basis with n actives.
1547 The hessian of the basis u has the form: hess(u)
1548 active 1 active 2 active n = cardinality
1549 [d11 u1, d11 u2, d11 u3] [d11 u1, d11 u2, d11 u3] ... [d11 u1, d11 u2, d11 u3]
1550 [d22 u1, d22 u2, d22 u3] [d22 u1, d22 u2, d22 u3] ... [d22 u1, d22 u2, d22 u3]
1551 [d12 u1, d12 u2, d12 u3] [d12 u1, d12 u2, d12 u3] ... [d12 u1, d12 u2, d12 u3]
1552
1553 Here, the basis function has components u = [u1,u2,u3]. Since they are evaluated for scalars
1554 we use blockDiag to make copies for all components ui
1555
1556 active 1 active 2 active k = cardinality/dim active 1 active 2k active 1 active 2k
1557 [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]
1558 [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]
1559 [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]
1560
1561 */
1562 const index_t numAct = u.data().values[0].rows(); // number of actives of a basis function
1563 const index_t cardinality = u.cardinality(); // total number of actives (=3*numAct)
1564
1565 res.resize(rows(), _u.dim() *_u.cardinality()); // (3 x 3*cardinality)
1566 res.setZero();
1567
1568 tmp = _u.data().values[2].reshapeCol(k, cols(), numAct );
1569 for (index_t d = 0; d != cols(); ++d)
1570 {
1571 const index_t s = d*(cardinality + 1);
1572 for (index_t i = 0; i != numAct; ++i)
1573 res.col(s+i*_u.cols()) = tmp.col(i);
1574 }
1575
1576 return res;
1577 }
1578
1579 template<class U> inline
1580 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
1581 rows_impl(const U & u) const
1582 {
1583 return u.source().domainDim() * ( u.source().domainDim() + 1 ) / 2;
1584 }
1585
1586 template<class U> inline
1587 typename util::enable_if< util::is_same<U,gsFeSolution<Scalar> >::value, index_t >::type
1588 rows_impl(const U & u) const
1589 {
1590 return u.parDim() * ( u.parDim() + 1 ) / 2;
1591 }
1592
1594 template<class U> inline
1595 typename util::enable_if< util::is_same<U,gsFeVariable<Scalar> >::value || util::is_same<U,gsGeometryMap<Scalar> >::value, index_t >::type
1596 cols_impl(const U & u) const
1597 {
1598 return u.source().targetDim();
1599 }
1600
1601 template<class U> inline
1602 typename util::enable_if< util::is_same<U,gsFeSolution<Scalar> >::value || util::is_same<U,gsFeSpace<Scalar> >::value, index_t >::type
1603 cols_impl(const U & u) const
1604 {
1605 return u.dim();
1606 }
1607
1608};
1609
1617template<class E1, class E2, class E3>
1618class flatdot_expr : public _expr<flatdot_expr<E1,E2,E3> >
1619{
1620public:
1621 typedef typename E1::Scalar Scalar;
1622
1623private:
1624 typename E1::Nested_t _A;
1625 typename E2::Nested_t _B;
1626 typename E3::Nested_t _C;
1627 mutable gsMatrix<Scalar> eA, eB, eC, tmp, res;
1628
1629public:
1630 enum {Space = 3, ScalarValued = 0, ColBlocks = 0};
1631
1632public:
1633
1634 flatdot_expr(_expr<E1> const& A, _expr<E2> const& B, _expr<E3> const& C) : _A(A),_B(B),_C(C)
1635 {
1636 }
1637
1638 const gsMatrix<Scalar> & eval(const index_t k) const
1639 {
1640 const index_t Ac = _A.cols();
1641 const index_t An = _A.cardinality();
1642 const index_t Bc = _B.cols();
1643 const index_t Bn = _B.cardinality();
1644
1645 eA = _A.eval(k);
1646 eB = _B.eval(k);
1647 eC = _C.eval(k);
1648
1649 GISMO_ASSERT(Bc==_A.rows(), "Dimensions: "<<Bc<<","<< _A.rows()<< "do not match");
1650 GISMO_STATIC_ASSERT(E1::Space==1, "First entry should be rowSpan" );
1651 GISMO_STATIC_ASSERT(E2::Space==2, "Second entry should be colSpan.");
1652
1653 res.resize(An, Bn);
1654
1655 // to do: evaluate the jacobians internally for the basis functions and loop over dimensions as well, since everything is same for all dimensions.
1656 for (index_t i = 0; i!=An; ++i) // for all actives u
1657 for (index_t j = 0; j!=Bn; ++j) // for all actives v
1658 {
1659 tmp.noalias() = eB.middleCols(i*Bc,Bc) * eA.middleCols(j*Ac,Ac);
1660
1661 tmp(0,0) *= eC.at(0);
1662 tmp(0,1) *= eC.at(2);
1663 tmp(1,0) *= eC.at(2);
1664 tmp(1,1) *= eC.at(1);
1665 res(i,j) = tmp.sum();
1666 }
1667 return res;
1668 }
1669
1670 index_t rows() const { return 1; }
1671 index_t cols() const { return 1; }
1672
1673 void parse(gsExprHelper<Scalar> & evList) const
1674 { _A.parse(evList);_B.parse(evList);_C.parse(evList); }
1675
1676 const gsFeSpace<Scalar> & rowVar() const { return _A.rowVar(); }
1677 const gsFeSpace<Scalar> & colVar() const { return _B.colVar(); }
1678 index_t cardinality_impl() const { return _A.cardinality_impl(); }
1679
1680 void print(std::ostream &os) const { os << "flatdot("; _A.print(os);_B.print(os);_C.print(os); os<<")"; }
1681};
1682
1692template<class E1, class E2, class E3>
1693class flatdot2_expr : public _expr<flatdot2_expr<E1,E2,E3> >
1694{
1695public:
1696 typedef typename E1::Scalar Scalar;
1697
1698private:
1699 typename E1::Nested_t _A;
1700 typename E2::Nested_t _B;
1701 typename E3::Nested_t _C;
1702 mutable gsMatrix<Scalar> eA, eB, eC, res, tmp;
1703
1704public:
1705 enum {Space = E1::Space, ScalarValued = 0, ColBlocks = 0};
1706
1707 flatdot2_expr(_expr<E1> const& A, _expr<E2> const& B, _expr<E3> const& C) : _A(A),_B(B),_C(C)
1708 {
1709 }
1710
1711 const gsMatrix<Scalar> & eval(const index_t k) const
1712 {
1713 const index_t Ac = _A.cols();
1714 const index_t An = _A.cardinality();
1715 const index_t Bn = _B.cardinality();
1716
1717 eA = _A.eval(k);
1718 eB = _B.eval(k);
1719 eC = _C.eval(k);
1720
1721 GISMO_ASSERT(_B.rows()==_A.cols(), "Dimensions: "<<_B.rows()<<","<< _A.cols()<< "do not match");
1722 GISMO_STATIC_ASSERT(E1::Space==1, "First entry should be rowSpan");
1723 GISMO_STATIC_ASSERT(E2::Space==2, "Second entry should be colSpan.");
1724 GISMO_ASSERT(_C.cols()==_B.rows(), "Dimensions: "<<_C.rows()<<","<< _B.rows()<< "do not match");
1725
1726 res.resize(An, Bn);
1727 for (index_t i = 0; i!=An; ++i)
1728 for (index_t j = 0; j!=Bn; ++j)
1729 {
1730 tmp = eA.middleCols(i*Ac,Ac) * eB.col(j); // E_f_der2
1731 tmp.row(2) *= 2.0; // multiply the third row of E_f_der2 by 2 for voight notation
1732 res(i,j) = (eC.row(0) * tmp.col(0)).value(); // E_f^T * mm * E_f_der2
1733 }
1734 return res;
1735 }
1736
1737 index_t rows() const { return 1; }
1738 index_t cols() const { return 1; }
1739
1740 void parse(gsExprHelper<Scalar> & evList) const
1741 { _A.parse(evList);_B.parse(evList);_C.parse(evList); }
1742
1743 const gsFeSpace<Scalar> & rowVar() const { return _A.rowVar(); }
1744 const gsFeSpace<Scalar> & colVar() const { return _B.colVar(); }
1745 index_t cardinality_impl() const { return _A.cardinality_impl(); }
1746
1747 void print(std::ostream &os) const { os << "flatdot2("; _A.print(os);_B.print(os);_C.print(os); os<<")"; }
1748};
1749
1754template<class T> class cartcovinv_expr ;
1755
1756template<class T>
1757class cartcov_expr : public _expr<cartcov_expr<T> >
1758{
1759 typename gsGeometryMap<T>::Nested_t _G;
1760
1761public:
1762 typedef T Scalar;
1763
1764 enum {Space = 0, ScalarValued = 0, ColBlocks = 0};
1765
1766 cartcov_expr(const gsGeometryMap<T> & G) : _G(G) { }
1767
1768 mutable gsMatrix<Scalar> covBasis;
1769 mutable gsMatrix<Scalar,3,3> result;
1770 mutable gsVector<Scalar> normal, tmp;
1771 mutable gsVector<Scalar> e1, e2, a1, a2;
1772
1773 gsMatrix<Scalar> eval(const index_t k) const
1774 {
1775 if (_G.targetDim()==3)
1776 {
1777 // Compute covariant bases in deformed and undeformed configuration
1778 normal = _G.data().normals.col(k);
1779 normal.normalize();
1780 covBasis.resize(3,3);
1781 covBasis.leftCols(2) = _G.data().jacobian(k);
1782 covBasis.col(2) = normal;
1783 }
1784 else if (_G.targetDim()==2)
1785 {
1786 // Compute covariant bases in deformed and undeformed configuration
1787 covBasis.resize(2,2);
1788 covBasis = _G.data().jacobian(k);
1789 }
1790 else
1791 GISMO_ERROR("Not implemented");
1792
1793 a1 = covBasis.col(0);
1794 a2 = covBasis.col(1);
1795
1796 e1 = a1.normalized();
1797 e2 = (a2-(a2.dot(e1)*e1)).normalized();
1798
1799 result(0,0) = (e1.dot(a1))*(a1.dot(e1)); // 1111
1800 result(0,1) = (e1.dot(a2))*(a2.dot(e1)); // 1122
1801 result(0,2) = 2*(e1.dot(a1))*(a2.dot(e1)); // 1112
1802 // Row 1
1803 result(1,0) = (e2.dot(a1))*(a1.dot(e2)); // 2211
1804 result(1,1) = (e2.dot(a2))*(a2.dot(e2)); // 2222
1805 result(1,2) = 2*(e2.dot(a1))*(a2.dot(e2)); // 2212
1806 // Row 2
1807 result(2,0) = (e1.dot(a1))*(a1.dot(e2)); // 1211
1808 result(2,1) = (e1.dot(a2))*(a2.dot(e2)); // 1222
1809 result(2,2) = (e1.dot(a1))*(a2.dot(e2)) + (e1.dot(a2))*(a1.dot(e2)); // 1212
1810 return result;
1811 }
1812
1813 cartcovinv_expr<T> inv() const
1814 {
1815 GISMO_ASSERT(rows() == cols(), "The Jacobian matrix is not square");
1816 return cartcovinv_expr<T>(_G);
1817 }
1818
1819 index_t rows() const { return 3; }
1820
1821 index_t cols() const { return 3; }
1822
1823 static const gsFeSpace<Scalar> & rowVar() { return gsNullExpr<Scalar>::get(); }
1824 static const gsFeSpace<Scalar> & colVar() { return gsNullExpr<Scalar>::get(); }
1825
1826 void parse(gsExprHelper<Scalar> & evList) const
1827 {
1828 //GISMO_ASSERT(NULL!=m_fd, "FeVariable: FuncData member not registered");
1829 evList.add(_G);
1830 _G.data().flags |= NEED_NORMAL|NEED_DERIV;
1831 }
1832
1833 void print(std::ostream &os) const { os << "cartcov("; _G.print(os); os <<")"; }
1834};
1835
1836template<class T>
1837class cartcovinv_expr : public _expr<cartcovinv_expr<T> >
1838{
1839 typename gsGeometryMap<T>::Nested_t _G;
1840
1841public:
1842 typedef T Scalar;
1843
1844 enum {Space = 0, ScalarValued = 0, ColBlocks = 0};
1845
1846 cartcovinv_expr(const gsGeometryMap<T> & G) : _G(G) { }
1847
1848 mutable gsMatrix<T> temp;
1849
1850 gsMatrix<T> eval(const index_t k) const
1851 {
1852 cartcov_expr<Scalar> cartcov = cartcov_expr<Scalar>(_G);
1853 temp = (cartcov.eval(k)).reshape(3,3).inverse();
1854 return temp;
1855 }
1856
1857 index_t rows() const { return 3; }
1858
1859 index_t cols() const { return 3; }
1860
1861 static const gsFeSpace<Scalar> & rowVar() { return gsNullExpr<Scalar>::get(); }
1862 static const gsFeSpace<Scalar> & colVar() { return gsNullExpr<Scalar>::get(); }
1863
1864 void parse(gsExprHelper<Scalar> & evList) const
1865 {
1866 //GISMO_ASSERT(NULL!=m_fd, "FeVariable: FuncData member not registered");
1867 cartcov(_G).parse(evList); //
1868
1869 evList.add(_G);
1870 _G.data().flags |= NEED_NORMAL|NEED_DERIV;
1871 }
1872
1873 void print(std::ostream &os) const { os << "cartcovinv("; _G.print(os); os <<")"; }
1874};
1875
1879template<class T> class cartconinv_expr ;
1880
1881template<class T>
1882class cartcon_expr : public _expr<cartcon_expr<T> >
1883{
1884 typename gsGeometryMap<T>::Nested_t _G;
1885
1886public:
1887 typedef T Scalar;
1888
1889 enum {Space = 0, ScalarValued = 0, ColBlocks = 0};
1890
1891 cartcon_expr(const gsGeometryMap<T> & G) : _G(G) { }
1892
1893 mutable gsMatrix<Scalar> covBasis, conBasis, covMetric, conMetric, cartBasis;
1894 mutable gsMatrix<Scalar,3,3> result;
1895 mutable gsVector<Scalar> normal, tmp;
1896 mutable gsVector<Scalar> e1, e2, ac1, ac2, a1, a2;
1897
1898 gsMatrix<Scalar> eval(const index_t k) const
1899 {
1900 if (_G.targetDim()==3)
1901 {
1902 // Compute covariant bases in deformed and undeformed configuration
1903 normal = _G.data().normals.col(k);
1904 normal.normalize();
1905
1906 covBasis.resize(3,3);
1907 conBasis.resize(3,3);
1908 covBasis.leftCols(2) = _G.data().jacobian(k);
1909 covBasis.col(2) = normal;
1910 covMetric = covBasis.transpose() * covBasis;
1911
1912 conMetric = covMetric.inverse();
1913
1914 conBasis.col(0) = conMetric(0,0)*covBasis.col(0)+conMetric(0,1)*covBasis.col(1)+conMetric(0,2)*covBasis.col(2);
1915 conBasis.col(1) = conMetric(1,0)*covBasis.col(0)+conMetric(1,1)*covBasis.col(1)+conMetric(1,2)*covBasis.col(2);
1916 }
1917 else if (_G.targetDim()==2)
1918 {
1919 // Compute covariant bases in deformed and undeformed configuration
1920 normal.resize(3);
1921 normal<<0,0,1;
1922 normal.normalize();
1923 covBasis.resize(2,2);
1924 conBasis.resize(2,2);
1925 // Compute covariant bases in deformed and undeformed configuration
1926 covBasis = _G.data().jacobian(k);
1927 covMetric = covBasis.transpose() * covBasis;
1928
1929 conMetric = covMetric.inverse();
1930
1931 conBasis.col(0) = conMetric(0,0)*covBasis.col(0)+conMetric(0,1)*covBasis.col(1);
1932 conBasis.col(1) = conMetric(1,0)*covBasis.col(0)+conMetric(1,1)*covBasis.col(1);
1933 }
1934 else
1935 GISMO_ERROR("Not implemented");
1936
1937 a1 = covBasis.col(0);
1938 a2 = covBasis.col(1);
1939
1940 e1 = a1.normalized();
1941 e2 = (a2-(a2.dot(e1)*e1)).normalized();
1942
1943 ac1 = conBasis.col(0);
1944 ac2 = conBasis.col(1);
1945
1946 result(0,0) = (e1.dot(ac1))*(ac1.dot(e1)); // 1111
1947 result(0,1) = (e1.dot(ac2))*(ac2.dot(e1)); // 1122
1948 result(0,2) = 2*(e1.dot(ac1))*(ac2.dot(e1)); // 1112
1949 // Row 1
1950 result(1,0) = (e2.dot(ac1))*(ac1.dot(e2)); // 2211
1951 result(1,1) = (e2.dot(ac2))*(ac2.dot(e2)); // 2222
1952 result(1,2) = 2*(e2.dot(ac1))*(ac2.dot(e2)); // 2212
1953 // Row 2
1954 result(2,0) = (e1.dot(ac1))*(ac1.dot(e2)); // 1211
1955 result(2,1) = (e1.dot(ac2))*(ac2.dot(e2)); // 1222
1956 result(2,2) = (e1.dot(ac1))*(ac2.dot(e2)) + (e1.dot(ac2))*(ac1.dot(e2)); // 1212
1957 return result;
1958
1959 }
1960
1961 cartconinv_expr<T> inv() const
1962 {
1963 GISMO_ASSERT(rows() == cols(), "The Jacobian matrix is not square");
1964 return cartconinv_expr<T>(_G);
1965 }
1966
1967 index_t rows() const { return 3; }
1968
1969 index_t cols() const { return 3; }
1970
1971 static const gsFeSpace<Scalar> & rowVar() { return gsNullExpr<Scalar>::get(); }
1972 static const gsFeSpace<Scalar> & colVar() { return gsNullExpr<Scalar>::get(); }
1973
1974 void parse(gsExprHelper<Scalar> & evList) const
1975 {
1976 //GISMO_ASSERT(NULL!=m_fd, "FeVariable: FuncData member not registered");
1977 evList.add(_G);
1978 _G.data().flags |= NEED_NORMAL|NEED_DERIV;
1979 }
1980
1981 void print(std::ostream &os) const { os << "cartcon("; _G.print(os); os <<")"; }
1982};
1983
1984template<class T>
1985class cartconinv_expr : public _expr<cartconinv_expr<T> >
1986{
1987private:
1988 typename gsGeometryMap<T>::Nested_t _G;
1989 mutable gsMatrix<T> temp;
1990
1991public:
1992 typedef T Scalar;
1993
1994 enum {Space = 0, ScalarValued = 0, ColBlocks = 0};
1995
1996 cartconinv_expr(const gsGeometryMap<T> & G) : _G(G) { }
1997
1998 gsMatrix<T> eval(const index_t k) const
1999 {
2000 cartcon_expr<Scalar> cartcon = cartcon_expr<Scalar>(_G);
2001 temp = (cartcon.eval(k)).reshape(3,3).inverse();
2002 return temp;
2003 }
2004
2005 index_t rows() const { return 3; }
2006
2007 index_t cols() const { return 3; }
2008
2009 static const gsFeSpace<Scalar> & rowVar() { return gsNullExpr<Scalar>::get(); }
2010 static const gsFeSpace<Scalar> & colVar() { return gsNullExpr<Scalar>::get(); }
2011
2012 void parse(gsExprHelper<Scalar> & evList) const
2013 {
2014 cartcon(_G).parse(evList); //
2015
2016 //GISMO_ASSERT(NULL!=m_fd, "FeVariable: FuncData member not registered");
2017 evList.add(_G);
2018 _G.data().flags |= NEED_NORMAL|NEED_DERIV;
2019 }
2020
2021 void print(std::ostream &os) const { os << "cartconinv("; _G.print(os); os <<")"; }
2022};
2023
2024EIGEN_STRONG_INLINE
2025unitVec_expr uv(const index_t index, const index_t dim) { return unitVec_expr(index,dim); }
2026
2027template<class E> EIGEN_STRONG_INLINE
2028var1_expr<E> var1(const E & u, const gsGeometryMap<typename E::Scalar> & G) { return var1_expr<E>(u, G); }
2029
2030template<class E1, class E2> EIGEN_STRONG_INLINE
2031var2_expr<E1,E2> var2(const E1 & u, const E2 & v, const gsGeometryMap<typename E1::Scalar> & G)
2032{ return var2_expr<E1,E2>(u,v, G); }
2033
2034template<class E1, class E2, class E3> EIGEN_STRONG_INLINE
2035var2dot_expr<E1,E2,E3> var2(const E1 & u, const E2 & v, const gsGeometryMap<typename E1::Scalar> & G, const E3 & Ef)
2036{ return var2dot_expr<E1,E2,E3>(u,v, G, Ef); }
2037
2038template<class E1, class E2, class E3> EIGEN_STRONG_INLINE
2039var2deriv2dot_expr<E1,E2,E3> var2deriv2(const E1 & u, const E2 & v, const gsGeometryMap<typename E1::Scalar> & G, const E3 & Ef)
2040{ return var2deriv2dot_expr<E1,E2,E3>(u,v, G, Ef); }
2041
2042template<class E> EIGEN_STRONG_INLINE
2043tvar1_expr<E> tvar1(const E & u, const gsGeometryMap<typename E::Scalar> & G) { return tvar1_expr<E>(u, G); }
2044
2045template<class E> EIGEN_STRONG_INLINE
2046ovar1_expr<E> ovar1(const E & u, const gsGeometryMap<typename E::Scalar> & G) { return ovar1_expr<E>(u, G); }
2047
2048template<class E1, class E2, class E3> EIGEN_STRONG_INLINE
2049ovar2dot_expr<E1,E2,E3> ovar2(const E1 & u, const E2 & v, const gsGeometryMap<typename E1::Scalar> & G, const E3 & C)
2050{ return ovar2dot_expr<E1,E2,E3>(u,v, G, C); }
2051
2052// template<class E1, class E2> EIGEN_STRONG_INLINE
2053// hessdot_expr<E1,E2> hessdot(const E1 & u, const E2 & v) { return hessdot_expr<E1,E2>(u, v); }
2054
2055template<class E> EIGEN_STRONG_INLINE
2056deriv2_expr<E> deriv2(const E & u) { return deriv2_expr<E>(u); }
2057
2058template<class E1, class E2> EIGEN_STRONG_INLINE
2059deriv2dot_expr<E1, E2> deriv2(const E1 & u, const E2 & v) { return deriv2dot_expr<E1, E2>(u,v); }
2060
2061template<class E1, class E2, class E3> EIGEN_STRONG_INLINE
2062flatdot_expr<E1,E2,E3> flatdot(const E1 & u, const E2 & v, const E3 & w)
2063{ return flatdot_expr<E1,E2,E3>(u, v, w); }
2064
2065template<class E1, class E2, class E3> EIGEN_STRONG_INLINE
2066flatdot2_expr<E1,E2,E3> flatdot2(const E1 & u, const E2 & v, const E3 & w)
2067{ return flatdot2_expr<E1,E2,E3>(u, v, w); }
2068
2069template<class E> EIGEN_STRONG_INLINE
2070cartcov_expr<E> cartcov(const gsGeometryMap<E> & G) { return cartcov_expr<E>(G); }
2071
2072template<class E> EIGEN_STRONG_INLINE
2073cartcon_expr<E> cartcon(const gsGeometryMap<E> & G) { return cartcon_expr<E>(G); }
2074
2075}
2076}
gismo::expr::cartconinv_expr
Definition gsThinShellUtils.h:1986
gismo::expr::cartcovinv_expr
Definition gsThinShellUtils.h:1838
gismo::expr::deriv2_expr
Computes the second derivative of an expression.
Definition gsThinShellUtils.h:1442
gismo::expr::deriv2_expr::eval_impl
util::enable_if< util::is_same< U, gsFeSpace< Scalar > >::value, constgsMatrix< Scalar > & >::type eval_impl(const U &u, const index_t k) const
Spexialization for a space.
Definition gsThinShellUtils.h:1543
gismo::expr::deriv2dot_expr
Expression that takes the second derivative of an expression and multiplies it with a row vector.
Definition gsThinShellUtils.h:1200
gismo::expr::flatdot2_expr
Computes the product of expressions E1 and E2 and multiplies with a vector E3 in voight notation.
Definition gsThinShellUtils.h:1694
gismo::expr::flatdot_expr
Computes the product of expressions E1 and E2 and multiplies with a vector E3 in voight notation.
Definition gsThinShellUtils.h:1619
gismo::expr::gsFeSpace
Definition gsExpressions.h:973
gismo::expr::ovar1_expr
Expression for the first variation of the outer normal.
Definition gsThinShellUtils.h:854
gismo::expr::ovar2dot_expr
Expression for the second variation of the outer normal times a vector.
Definition gsThinShellUtils.h:1002
gismo::expr::tangent_expr
Definition gsExpressions.h:3081
gismo::expr::tvar1_expr
Expression for the first variation of the outer tangent.
Definition gsThinShellUtils.h:660
gismo::expr::unitVec_expr
Simple expression for the unit vector of length dim and with value 1 on index.
Definition gsThinShellUtils.h:41
gismo::expr::var1_expr
Expression for the first variation of the surface normal.
Definition gsThinShellUtils.h:81
gismo::expr::var2_expr
Second variation of the normal.
Definition gsThinShellUtils.h:500
gismo::expr::var2deriv2dot_expr
Second variation of the surface normal times the second derivative of the geometry map times a vector...
Definition gsThinShellUtils.h:376
gismo::expr::var2dot_expr
Second variation of the surface normal times a vector.
Definition gsThinShellUtils.h:249
gismo::gsExprHelper
Definition gsExprHelper.h:27
gismo::gsMatrix
A matrix with arbitrary coefficient type and fixed or dynamic size.
Definition gsMatrix.h:41
gismo::gsMatrix::reshapeCol
gsAsMatrix< T, Dynamic, Dynamic > reshapeCol(index_t c, index_t n, index_t m)
Returns column c of the matrix resized to n x m matrix This function assumes that the matrix is size ...
Definition gsMatrix.h:231
gismo::gsMatrix::col3d
Col3DType col3d(index_t c)
Returns column c as a fixed-size 3D vector.
Definition gsMatrix.h:244
gismo::gsMatrix::at
T at(index_t i) const
Returns the i-th element of the vectorization of the matrix.
Definition gsMatrix.h:211
gismo::gsVector
A vector with arbitrary coefficient type and fixed or dynamic size.
Definition gsVector.h:37
gismo::gsVector::at
T at(index_t i) const
Returns the i-th element of the vector.
Definition gsVector.h:177
gsAssemblerOptions.h
Provides assembler and solver options.
gsBasis.h
Provides declaration of Basis abstract interface.
gsBoundaryConditions.h
Provides gsBoundaryConditions class.
short_t
#define short_t
Definition gsConfig.h:35
index_t
#define index_t
Definition gsConfig.h:32
GISMO_NO_IMPLEMENTATION
#define GISMO_NO_IMPLEMENTATION
Definition gsDebug.h:129
GISMO_ERROR
#define GISMO_ERROR(message)
Definition gsDebug.h:118
gsWarn
#define gsWarn
Definition gsDebug.h:50
GISMO_UNUSED
#define GISMO_UNUSED(x)
Definition gsDebug.h:112
gsInfo
#define gsInfo
Definition gsDebug.h:43
GISMO_ASSERT
#define GISMO_ASSERT(cond, message)
Definition gsDebug.h:89
gsDofMapper.h
Provides the gsDofMapper class for re-indexing DoFs.
gsExpressions.h
Defines different expressions.
gsFuncData.h
This object is a cache for computed values from an evaluator.
gsLinearAlgebra.h
This is the main header file that collects wrappers of Eigen for linear algebra.
gsPointGrid.h
Provides functions to generate structured point data.
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::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::grad
EIGEN_STRONG_INLINE grad_expr< E > grad(const E &u)
The gradient of a variable.
Definition gsExpressions.h:4490
gismo
The G+Smo namespace, containing all definitions for the library.
gismo::NEED_VALUE
@ NEED_VALUE
Value of the object.
Definition gsForwardDeclarations.h:72
gismo::NEED_DERIV2
@ NEED_DERIV2
Second derivatives.
Definition gsForwardDeclarations.h:80
gismo::NEED_DERIV
@ NEED_DERIV
Gradient of the object.
Definition gsForwardDeclarations.h:73
gismo::NEED_NORMAL
@ NEED_NORMAL
Normal vector of the object.
Definition gsForwardDeclarations.h:85
gismo::NEED_OUTER_NORMAL
@ NEED_OUTER_NORMAL
Outward normal on the boundary.
Definition gsForwardDeclarations.h:86
gismo::NEED_ACTIVE
@ NEED_ACTIVE
Active function ids.
Definition gsForwardDeclarations.h:84