1 /** @file gsThinShellUtils.h
2
3 @brief Utilities for gsThinShellAssembler. Mainly expressions
4
5 This file is part of the G+Smo library.
6
7 This Source Code Form is subject to the terms of the Mozilla Public
8 License, v. 2.0. If a copy of the MPL was not distributed with this
9 file, You can obtain one at http://mozilla.org/MPL/2.0/.
10
11 Author(s):
12 H.M. Verhelst (2019-..., TU Delft)
13 A. Mantzaflaris (2019-..., Inria)
14 */
15
16 //! [Include namespace]
17 #include <gsCore/gsLinearAlgebra.h>
18 #include <gsCore/gsBasis.h>
19 #include <gsCore/gsFuncData.h>
20 #include <gsCore/gsDofMapper.h>
21
22 #include <gsPde/gsBoundaryConditions.h>
23
24 #include <gsUtils/gsPointGrid.h>
25
26 #include <gsAssembler/gsAssemblerOptions.h>
27 #include <gsAssembler/gsExpressions.h>
28
29
30 # define MatExprType auto
31
32 namespace gismo{
33 namespace expr{
34
35 /**
36 * @brief Simple expression for the unit vector of length \a dim and with value 1 on \a index
37 */
38 class unitVec_expr : public _expr<unitVec_expr >
39 {
40 public:
41 typedef real_t Scalar;
42 private:
43 index_t _index;
44 index_t _dim;
45
46 public:
unitVec_expr(const index_t index,const index_t dim)47 unitVec_expr(const index_t index, const index_t dim) : _index(index), _dim(dim) { }
48
49 public:
50 enum{ Space = 0, ScalarValued= 0, ColBlocks= 0};
51
eval(const index_t)52 gsMatrix<Scalar> eval(const index_t) const
53 {
54 gsMatrix<Scalar> vec = gsMatrix<Scalar>::Zero(_dim,1);
55 // vec.setZero();
56 vec(_index,0) = 1;
57 return vec;
58 }
59
rows()60 index_t rows() const { return _dim; }
cols()61 index_t cols() const { return 1; }
parse(gsSortedVector<const gsFunctionSet<Scalar> * > &)62 void parse(gsSortedVector<const gsFunctionSet<Scalar>*> & ) const { }
63
rowVar()64 const gsFeSpace<Scalar> & rowVar() const {return gsNullExpr<Scalar>::get();}
colVar()65 const gsFeSpace<Scalar> & colVar() const {return gsNullExpr<Scalar>::get();}
66
print(std::ostream & os)67 void print(std::ostream &os) const { os << "uv("<<_dim <<")";}
68 };
69
70
71 /**
72 * @brief Expression for the first variation of the surface normal
73 *
74 * @tparam E Object type
75 */
76 template<class E>
77 class var1_expr : public _expr<var1_expr<E> >
78 {
79 public:
80 typedef typename E::Scalar Scalar;
81
82 private:
83 typename E::Nested_t _u;
84 typename gsGeometryMap<Scalar>::Nested_t _G;
85
86 mutable gsMatrix<Scalar> res;
87 mutable gsMatrix<Scalar> bGrads, cJac;
88 mutable gsVector<Scalar,3> m_v, normal;
89 public:
90 enum{ Space = E::Space, ScalarValued= 0, ColBlocks= 0};
91
var1_expr(const E & u,const gsGeometryMap<Scalar> & G)92 var1_expr(const E & u, const gsGeometryMap<Scalar> & G) : _u(u), _G(G) { }
93
94 EIGEN_MAKE_ALIGNED_OPERATOR_NEW
95
96 // helper function
vecFun(index_t pos,Scalar val)97 static inline gsVector<Scalar,3> vecFun(index_t pos, Scalar val)
98 {
99 gsVector<Scalar,3> result = gsVector<Scalar,3>::Zero();
100 result[pos] = val;
101 return result;
102 }
103
eval(const index_t k)104 const gsMatrix<Scalar> & eval(const index_t k) const {return eval_impl(_u,k); }
105
rows()106 index_t rows() const { return 1; }
cols()107 index_t cols() const { return _u.dim(); }
108
parse(gsExprHelper<Scalar> & evList)109 void parse(gsExprHelper<Scalar> & evList) const
110 {
111 parse_impl<E>(evList);
112 }
113
rowVar()114 const gsFeSpace<Scalar> & rowVar() const { return _u.rowVar(); }
colVar()115 const gsFeSpace<Scalar> & colVar() const {return gsNullExpr<Scalar>::get();}
cardinality_impl()116 index_t cardinality_impl() const { return _u.cardinality_impl(); }
117
print(std::ostream & os)118 void print(std::ostream &os) const { os << "var1("; _u.print(os); os <<")"; }
119
120 private:
121 template<class U> inline
122 typename util::enable_if< !util::is_same<U,gsFeSolution<Scalar> >::value,void>::type
parse_impl(gsExprHelper<Scalar> & evList)123 parse_impl(gsExprHelper<Scalar> & evList) const
124 {
125 evList.add(_u);
126 _u.data().flags |= NEED_GRAD;
127 evList.add(_G);
128 _G.data().flags |= NEED_NORMAL | NEED_DERIV | NEED_MEASURE;
129 }
130
131 template<class U> inline
132 typename util::enable_if< util::is_same<U,gsFeSolution<Scalar> >::value,void>::type
parse_impl(gsExprHelper<Scalar> & evList)133 parse_impl(gsExprHelper<Scalar> & evList) const
134 {
135 evList.add(_G);
136 _G.data().flags |= NEED_NORMAL | NEED_DERIV | NEED_MEASURE;
137
138 grad(_u).parse(evList); //
139
140 _u.parse(evList);
141 }
142
143 template<class U> inline
144 typename util::enable_if< util::is_same<U,gsFeSpace<Scalar> >::value, const gsMatrix<Scalar> & >::type
eval_impl(const U & u,const index_t k)145 eval_impl(const U & u, const index_t k) const
146 {
147 const index_t A = _u.cardinality()/_u.dim(); // _u.data().actives.rows()
148 res.resize(_u.cardinality(), cols()); // rows()*
149
150 normal = _G.data().normal(k);// not normalized to unit length
151 normal.normalize();
152 bGrads = _u.data().values[1].col(k);
153 cJac = _G.data().values[1].reshapeCol(k, _G.data().dim.first, _G.data().dim.second).transpose();
154 const Scalar measure = _G.data().measures.at(k);
155
156 for (index_t d = 0; d!= cols(); ++d) // for all basis function components
157 {
158 const short_t s = d*A;
159 for (index_t j = 0; j!= A; ++j) // for all actives
160 {
161 // Jac(u) ~ Jac(G) with alternating signs ?..
162 m_v.noalias() = (vecFun(d, bGrads.at(2*j ) ).cross( cJac.col(1).template head<3>() )
163 - vecFun(d, bGrads.at(2*j+1) ).cross( cJac.col(0).template head<3>() )) / measure;
164
165 // --------------- First variation of the normal
166 // res.row(s+j).noalias() = (m_v - ( normal.dot(m_v) ) * normal).transpose();
167 res.row(s+j).noalias() = (m_v - ( normal*m_v.transpose() ) * normal).transpose(); // outer-product version
168 }
169 }
170 return res;
171 }
172
173 template<class U> inline
174 typename util::enable_if< util::is_same<U,gsFeSolution<Scalar> >::value, const gsMatrix<Scalar> & >::type
eval_impl(const U & u,const index_t k)175 eval_impl(const U & u, const index_t k) const
176 {
177 GISMO_ASSERT(1==_u.data().actives.cols(), "Single actives expected");
178 grad_expr<gsFeSolution<Scalar>> sGrad = grad_expr<gsFeSolution<Scalar>>(_u);
179 res.resize(rows(), cols()); // rows()*
180
181 normal = _G.data().normal(k);// not normalized to unit length
182 normal.normalize();
183 bGrads = sGrad.eval(k);
184 cJac = _G.data().values[1].reshapeCol(k, _G.data().dim.first, _G.data().dim.second).transpose();
185 const Scalar measure = _G.data().measures.at(k);
186
187 m_v.noalias() = ( ( bGrads.col(0).template head<3>() ).cross( cJac.col(1).template head<3>() )
188 - ( bGrads.col(1).template head<3>() ).cross( cJac.col(0).template head<3>() ) ) / measure;
189
190 // --------------- First variation of the normal
191 // res.row(s+j).noalias() = (m_v - ( normal.dot(m_v) ) * normal).transpose();
192 res = (m_v - ( normal*m_v.transpose() ) * normal).transpose(); // outer-product version
193 return res;
194 }
195 };
196
197
198 /**
199 * @brief Second variation of the surface normal times a vector.
200 *
201 * @tparam E1 Type of u
202 * @tparam E2 Type of v
203 * @tparam E3 Type of the vector
204 */
205 template<class E1, class E2, class E3>
206 class var2_expr : public _expr<var2_expr<E1,E2,E3> >
207 {
208 public:
209 typedef typename E1::Scalar Scalar;
210
211 private:
212 typename E1::Nested_t _u;
213 typename E2::Nested_t _v;
214 typename gsGeometryMap<Scalar>::Nested_t _G;
215 typename E3::Nested_t _Ef;
216
217 public:
218 enum{ Space = 3, ScalarValued= 0, ColBlocks= 0 };
219
var2_expr(const E1 & u,const E2 & v,const gsGeometryMap<Scalar> & G,_expr<E3> const & Ef)220 var2_expr( const E1 & u, const E2 & v, const gsGeometryMap<Scalar> & G, _expr<E3> const& Ef) : _u(u),_v(v), _G(G), _Ef(Ef) { }
221
222 mutable gsMatrix<Scalar> res;
223
224 mutable gsMatrix<Scalar> uGrads, vGrads, cJac, cDer2, evEf, result;
225 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
226 EIGEN_MAKE_ALIGNED_OPERATOR_NEW
227
228 // helper function
vecFun(index_t pos,Scalar val)229 static inline gsVector<Scalar,3> vecFun(index_t pos, Scalar val)
230 {
231 gsVector<Scalar,3> result = gsVector<Scalar,3>::Zero();
232 result[pos] = val;
233 return result;
234 }
235
eval(const index_t k)236 const gsMatrix<Scalar> & eval(const index_t k) const
237 {
238 res.resize(_u.cardinality(), _u.cardinality());
239
240 normal = _G.data().normal(k);
241 normal.normalize();
242 uGrads = _u.data().values[1].col(k);
243 vGrads = _v.data().values[1].col(k);
244 cJac = _G.data().values[1].reshapeCol(k, _G.data().dim.first, _G.data().dim.second).transpose();
245 cDer2 = _G.data().values[2].reshapeCol(k, _G.data().dim.second, _G.data().dim.second);
246
247 const index_t cardU = _u.data().values[0].rows(); // number of actives per component of u
248 const index_t cardV = _v.data().values[0].rows(); // number of actives per component of v
249 const Scalar measure = _G.data().measures.at(k);
250
251 evEf = _Ef.eval(k);
252
253 for (index_t j = 0; j!= cardU; ++j) // for all basis functions u (1)
254 {
255 for (index_t i = 0; i!= cardV; ++i) // for all basis functions v (1)
256 {
257 for (index_t d = 0; d!= _u.dim(); ++d) // for all basis functions u (2)
258 {
259 m_u.noalias() = ( vecFun(d, uGrads.at(2*j ) ).cross( cJac.col(1).template head<3>() )
260 -vecFun(d, uGrads.at(2*j+1) ).cross( cJac.col(0).template head<3>() ))
261 / measure;
262
263 const short_t s = d*cardU;
264
265 for (index_t c = 0; c!= _v.dim(); ++c) // for all basis functions v (2)
266 {
267 const short_t r = c*cardV;
268 m_v.noalias() = ( vecFun(c, vGrads.at(2*i ) ).cross( cJac.col(1).template head<3>() )
269 -vecFun(c, vGrads.at(2*i+1) ).cross( cJac.col(0).template head<3>() ))
270 / measure;
271
272 // n_der.noalias() = (m_v - ( normal.dot(m_v) ) * normal);
273 n_der.noalias() = (m_v - ( normal*m_v.transpose() ) * normal); // outer-product version
274
275 m_uv.noalias() = ( vecFun(d, uGrads.at(2*j ) ).cross( vecFun(c, vGrads.at(2*i+1) ) )
276 +vecFun(c, vGrads.at(2*i ) ).cross( vecFun(d, uGrads.at(2*j+1) ) ))
277 / measure; //check
278
279 m_u_der.noalias() = (m_uv - ( normal.dot(m_v) ) * m_u);
280 // m_u_der.noalias() = (m_uv - ( normal*m_v.transpose() ) * m_u); // outer-product version TODO
281
282 // --------------- Second variation of the normal
283 tmp = m_u_der - (m_u.dot(n_der) + normal.dot(m_u_der) ) * normal - (normal.dot(m_u) ) * n_der;
284 // tmp = m_u_der - (m_u.dot(n_der) + normal.dot(m_u_der) ) * normal - (normal.dot(m_u) ) * n_der;
285
286 // Evaluate the product
287 tmp = cDer2 * tmp; // E_f_der2, last component
288 tmp.row(2) *= 2.0;
289 result = evEf * tmp;
290
291 res(s + j, r + i ) = result(0,0);
292 }
293 }
294 }
295 }
296 return res;
297 }
298
rows()299 index_t rows() const
300 {
301 return 1; // because the resulting matrix has scalar entries for every combination of active basis functions
302 }
303
cols()304 index_t cols() const
305 {
306 return 1; // because the resulting matrix has scalar entries for every combination of active basis functions
307 }
308
parse(gsExprHelper<Scalar> & evList)309 void parse(gsExprHelper<Scalar> & evList) const
310 {
311 evList.add(_u);
312 _u.data().flags |= NEED_ACTIVE | NEED_VALUE | NEED_GRAD;
313 evList.add(_v);
314 _v.data().flags |= NEED_ACTIVE | NEED_VALUE | NEED_GRAD;
315 evList.add(_G);
316 _G.data().flags |= NEED_NORMAL | NEED_DERIV | NEED_2ND_DER | NEED_MEASURE;
317 _Ef.parse(evList);
318 }
319
rowVar()320 const gsFeSpace<Scalar> & rowVar() const { return _u.rowVar(); }
colVar()321 const gsFeSpace<Scalar> & colVar() const { return _v.rowVar(); }
322
print(std::ostream & os)323 void print(std::ostream &os) const { os << "var2("; _u.print(os); os <<")"; }
324 };
325
326 /**
327 * @brief Expression for the first variation of the outer tangent
328 *
329 * @tparam E Object type
330 */
331 template<class E>
332 class tvar1_expr : public _expr<tvar1_expr<E> >
333 {
334 public:
335 typedef typename E::Scalar Scalar;
336
337 private:
338 typename E::Nested_t _u;
339 typename gsGeometryMap<Scalar>::Nested_t _G;
340
341 mutable gsVector<Scalar,3> onormal, tangent, dtan;
342 mutable gsVector<Scalar> tmp;
343 mutable gsMatrix<Scalar> bGrads, cJac, res;
344
345 public:
346 enum{ Space = E::Space, ScalarValued= 0, ColBlocks= 0};
347
tvar1_expr(const E & u,const gsGeometryMap<Scalar> & G)348 tvar1_expr(const E & u, const gsGeometryMap<Scalar> & G) : _u(u), _G(G) { }
349
350 EIGEN_MAKE_ALIGNED_OPERATOR_NEW
351
352 // helper function
vecFun(index_t pos,Scalar val)353 static inline gsVector<Scalar,3> vecFun(index_t pos, Scalar val)
354 {
355 gsVector<Scalar,3> result = gsVector<Scalar,3>::Zero();
356 result[pos] = val;
357 return result;
358 }
359
eval(const index_t k)360 const gsMatrix<Scalar> & eval(const index_t k) const {return eval_impl(_u,k); }
361
rows()362 index_t rows() const { return 1; }
cols()363 index_t cols() const { return _u.dim(); }
364
parse(gsExprHelper<Scalar> & evList)365 void parse(gsExprHelper<Scalar> & evList) const
366 {
367 parse_impl<E>(evList);
368 }
369
rowVar()370 const gsFeSpace<Scalar> & rowVar() const { return _u.rowVar(); }
colVar()371 const gsFeSpace<Scalar> & colVar() const {return gsNullExpr<Scalar>::get();}
cardinality_impl()372 index_t cardinality_impl() const { return _u.cardinality_impl(); }
373
print(std::ostream & os)374 void print(std::ostream &os) const { os << "tvar("; _G.print(os); os <<")"; }
375
376 private:
377 template<class U> inline
378 typename util::enable_if< !util::is_same<U,gsFeSolution<Scalar> >::value,void>::type
parse_impl(gsExprHelper<Scalar> & evList)379 parse_impl(gsExprHelper<Scalar> & evList) const
380 {
381 evList.add(_u);
382 _u.data().flags |= NEED_GRAD | NEED_ACTIVE;
383 evList.add(_G);
384 _G.data().flags |= NEED_NORMAL | NEED_OUTER_NORMAL | NEED_DERIV | NEED_MEASURE;
385 }
386
387 template<class U> inline
388 typename util::enable_if< util::is_same<U,gsFeSolution<Scalar> >::value,void>::type
parse_impl(gsExprHelper<Scalar> & evList)389 parse_impl(gsExprHelper<Scalar> & evList) const
390 {
391 evList.add(_G);
392 _G.data().flags |= NEED_NORMAL | NEED_OUTER_NORMAL | NEED_DERIV | NEED_MEASURE;
393
394 grad(_u).parse(evList); //
395
396 _u.parse(evList);
397 }
398
399 template<class U> inline
400 typename util::enable_if< !util::is_same<U,gsFeSolution<Scalar> >::value, const gsMatrix<Scalar> & >::type
eval_impl(const U & u,const index_t k)401 eval_impl(const U & u, const index_t k) const
402 {
403 GISMO_ASSERT(_G.data().dim.second==3,"Domain dimension should be 3, is "<<_G.data().dim.second);
404
405 const index_t A = _u.cardinality()/_u.dim();
406 res.resize(_u.cardinality(), cols()); // rows()*
407 cJac = _G.data().jacobian(k);
408
409 onormal = _G.data().outNormal(k);
410 tmp = cJac.transpose() * onormal;
411 Scalar tol = 1e-8;
412
413 /*
414 We can check which column of the Jacobian corresponds to the outer normal vector or to the tangent.
415 The tangent is a covariant vector and hence the column of the Jacobian should be equal to the tangent.
416 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.
417 */
418 index_t colIndex = -1;
419 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
420 colIndex = 0;
421 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
422 colIndex = 1;
423 else // then the normal is unknown??
424 gsInfo<<"warning: choice unknown\n";
425
426 tangent = cJac.col(colIndex);
427
428 // Now we will compute the derivatives of the basis functions
429 bGrads = _u.data().values[1].col(k);
430 for (index_t d = 0; d!= cols(); ++d) // for all basis function components
431 {
432 const short_t s = d*A;
433 for (index_t j = 0; j!= A; ++j) // for all actives
434 {
435 // The tangent vector is in column colIndex in cJac and thus in 2*j+colIndex in bGrads.
436 // Furthermore, as basis function for dimension d, it has a nonzero in entry d, and zeros elsewhere
437 dtan = vecFun(d, bGrads.at(2*j+colIndex));
438 res.row(s+j).noalias() = (1 / tangent.norm() * ( dtan - ( tangent.transpose() * dtan ) * tangent / (tangent.norm() * tangent.norm()) )).transpose();
439 }
440 }
441 return res;
442 }
443
444 template<class U> inline
445 typename util::enable_if< util::is_same<U,gsFeSolution<Scalar> >::value, const gsMatrix<Scalar> & >::type
eval_impl(const U & u,const index_t k)446 eval_impl(const U & u, const index_t k) const
447 {
448 GISMO_ASSERT(1==_u.data().actives.cols(), "Single actives expected");
449 res.resize(rows(), cols());
450
451 cJac = _G.data().jacobian(k);
452 onormal = _G.data().outNormal(k);
453 tmp = cJac.transpose() * onormal;
454 Scalar tol = 1e-8;
455
456 /*
457 We can check which column of the Jacobian corresponds to the outer normal vector or to the tangent.
458 The tangent is a covariant vector and hence the column of the Jacobian should be equal to the tangent.
459 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.
460 */
461 index_t colIndex;
462 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
463 colIndex = 0;
464 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
465 colIndex = 1;
466 else // then the normal is unknown??
467 gsInfo<<"warning: choice unknown\n";
468
469 tangent = cJac.col(colIndex);
470 bGrads = _u.data().values[1].col(k);
471 dtan = bGrads.col(colIndex);
472 res.noalias() = (1 / tangent.norm() * ( dtan - ( tangent * dtan ) * tangent / (tangent.norm() * tangent.norm()) )).transpose();
473 return res;
474 }
475
476 };
477
478 /**
479 * @brief Expression for the first variation of the outer normal
480 *
481 * @tparam E Object type
482 */
483 template<class E>
484 class ovar1_expr : public _expr<ovar1_expr<E> >
485 {
486 public:
487 typedef typename E::Scalar Scalar;
488
489 private:
490 typename E::Nested_t _u;
491 typename gsGeometryMap<Scalar>::Nested_t _G;
492
493 mutable gsVector<Scalar,3> onormal, tangent, utangent, normal, dtan, tvar, snvar, m_v;
494 mutable gsVector<Scalar> tmp;
495 mutable gsMatrix<Scalar> bGrads, cJac, res;
496
497 public:
498 enum{ Space = E::Space, ScalarValued= 0, ColBlocks= 0};
499
ovar1_expr(const E & u,const gsGeometryMap<Scalar> & G)500 ovar1_expr(const E & u, const gsGeometryMap<Scalar> & G) : _u(u), _G(G) { }
501
502 EIGEN_MAKE_ALIGNED_OPERATOR_NEW
503
504 // helper function
vecFun(index_t pos,Scalar val)505 static inline gsVector<Scalar,3> vecFun(index_t pos, Scalar val)
506 {
507 gsVector<Scalar,3> result = gsVector<Scalar,3>::Zero();
508 result[pos] = val;
509 return result;
510 }
511
eval(const index_t k)512 const gsMatrix<Scalar> & eval(const index_t k) const { return eval_impl(_u,k); }
513
rows()514 index_t rows() const { return 1; }
cols()515 index_t cols() const { return _u.dim(); }
516
parse(gsExprHelper<Scalar> & evList)517 void parse(gsExprHelper<Scalar> & evList) const
518 {
519 parse_impl<E>(evList);
520 }
521
rowVar()522 const gsFeSpace<Scalar> & rowVar() const { return _u.rowVar(); }
colVar()523 const gsFeSpace<Scalar> & colVar() const {return gsNullExpr<Scalar>::get();}
cardinality_impl()524 index_t cardinality_impl() const { return _u.cardinality_impl(); }
525
print(std::ostream & os)526 void print(std::ostream &os) const { os << "ovar("; _G.print(os); os <<")"; }
527
528 private:
529 template<class U> inline
530 typename util::enable_if< !util::is_same<U,gsFeSolution<Scalar> >::value,void>::type
parse_impl(gsExprHelper<Scalar> & evList)531 parse_impl(gsExprHelper<Scalar> & evList) const
532 {
533 evList.add(_u);
534 _u.data().flags |= NEED_GRAD | NEED_ACTIVE;
535 evList.add(_G);
536 _G.data().flags |= NEED_NORMAL | NEED_OUTER_NORMAL | NEED_DERIV | NEED_MEASURE;
537 }
538
539 template<class U> inline
540 typename util::enable_if< util::is_same<U,gsFeSolution<Scalar> >::value,void>::type
parse_impl(gsExprHelper<Scalar> & evList)541 parse_impl(gsExprHelper<Scalar> & evList) const
542 {
543 evList.add(_G);
544 _G.data().flags |= NEED_NORMAL | NEED_OUTER_NORMAL | NEED_DERIV | NEED_MEASURE;
545
546 grad(_u).parse(evList); //
547
548 _u.parse(evList);
549 }
550
551 template<class U> inline
552 typename util::enable_if< !util::is_same<U,gsFeSolution<Scalar> >::value, const gsMatrix<Scalar> & >::type
eval_impl(const U & u,const index_t k)553 eval_impl(const U & u, const index_t k) const
554 {
555 GISMO_ASSERT(_G.data().dim.second==3,"Domain dimension should be 3, is "<<_G.data().dim.second);
556 const index_t A = _u.cardinality()/_u.dim(); // _u.data().actives.rows()
557 res.resize(_u.cardinality(), cols()); // rows()*
558
559 onormal = _G.data().outNormal(k);
560 cJac = _G.data().values[1].reshapeCol(k, _G.data().dim.first, _G.data().dim.second).transpose();
561 tmp = cJac.transpose() * onormal;
562 Scalar tol = 1e-8;
563
564 /*
565 We can check which column of the Jacobian corresponds to the outer normal vector or to the tangent.
566 The tangent is a covariant vector and hence the column of the Jacobian should be equal to the tangent.
567 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.
568 */
569 index_t colIndex = -1;
570 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
571 colIndex = 0;
572 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
573 colIndex = 1;
574 else // then the normal is unknown??
575 gsInfo<<"warning: choice unknown\n";
576
577 tangent = cJac.col(colIndex);
578 utangent = tangent / tangent.norm();
579
580 // For the normal vector variation
581 normal = _G.data().normal(k);
582 normal.normalize();
583 bGrads = _u.data().values[1].col(k);
584 const Scalar measure = _G.data().measures.at(k);
585
586 for (index_t d = 0; d!= cols(); ++d) // for all basis function components
587 {
588 const short_t s = d*A;
589 for (index_t j = 0; j!= A; ++j) // for all actives
590 {
591 // VARIATION OF THE TANGENT
592 // The tangent vector is in column colIndex in cJac and thus in 2*j+colIndex in bGrads.
593 // Furthermore, as basis function for dimension d, it has a nonzero in entry d, and zeros elsewhere
594 dtan = vecFun(d, bGrads.at(2*j+colIndex));
595 tvar.noalias() = 1 / tangent.norm() * ( dtan - ( utangent.dot(dtan) ) * utangent);
596
597 // VARIATION OF THE NORMAL
598 // Jac(u) ~ Jac(G) with alternating signs ?..
599 m_v.noalias() = (vecFun(d, bGrads.at(2*j ) ).cross( cJac.col(1).template head<3>() )
600 - vecFun(d, bGrads.at(2*j+1) ).cross( cJac.col(0).template head<3>() )) / measure;
601
602 // --------------- First variation of the normal
603 snvar.noalias() = m_v - ( normal*m_v.transpose() ) * normal;
604
605 // VARIATION OF THE OUTER NORMAL
606 res.row(s+j).noalias() = tvar.cross(normal) + utangent.cross(snvar);
607 }
608 }
609 return res;
610 }
611
612 template<class U> inline
613 typename util::enable_if< util::is_same<U,gsFeSolution<Scalar> >::value, const gsMatrix<Scalar> & >::type
eval_impl(const U & u,const index_t k)614 eval_impl(const U & u, const index_t k) const
615 {
616 GISMO_ASSERT(_G.data().dim.second==3,"Domain dimension should be 3, is "<<_G.data().dim.second);
617 res.resize(rows(), cols());
618
619 cJac = _G.data().jacobian(k);
620 onormal = _G.data().outNormal(k);
621 normal = _G.data().normal(k);
622 tmp = cJac.transpose() * onormal;
623 Scalar tol = 1e-8;
624
625 /*
626 We can check which column of the Jacobian corresponds to the outer normal vector or to the tangent.
627 The tangent is a covariant vector and hence the column of the Jacobian should be equal to the tangent.
628 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.
629 */
630 index_t colIndex;
631 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
632 colIndex = 0;
633 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
634 colIndex = 1;
635 else // then the normal is unknown??
636 gsInfo<<"warning: choice unknown\n";
637
638 tangent = cJac.col(colIndex);
639 utangent = tangent / tangent.norm();
640
641 // Now we will compute the derivatives of the basis functions
642 bGrads = _u.data().values[1].col(k);
643 const Scalar measure = _G.data().measures.at(k);
644
645 // VARIATION OF THE TANGENT
646 // The tangent vector is in column colIndex in cJac and colIndex in bGrads.
647 dtan = bGrads.col(colIndex);
648 tvar.noalias() = 1 / tangent.norm() * ( dtan - ( utangent.dot(dtan) ) * utangent );
649
650 // VARIATION OF THE NORMAL
651 m_v.noalias() = ( ( bGrads.col(0).template head<3>() ).cross( cJac.col(1).template head<3>() )
652 - ( bGrads.col(1).template head<3>() ).cross( cJac.col(0).template head<3>() ) ) / measure;
653
654 // --------------- First variation of the normal
655 snvar.noalias() = m_v - ( normal.dot(m_v) ) * normal;
656
657 // VARIATION OF THE OUTER NORMAL
658 res.noalias() = tvar.cross(normal) + utangent.cross(snvar);
659
660 return res;
661 }
662 };
663
664 /**
665 * @brief Expression for the second variation of the outer normal times a vector
666 *
667 * @tparam E1 Type of u
668 * @tparam E2 Type of v
669 * @tparam E3 Type of the vector
670 */
671 template<class E1, class E2, class E3>
672 class ovar2_expr : public _expr<ovar2_expr<E1,E2,E3> >
673 {
674 public:
675 typedef typename E1::Scalar Scalar;
676
677 private:
678 typename E1::Nested_t _u;
679 typename E2::Nested_t _v;
680 typename gsGeometryMap<Scalar>::Nested_t _G;
681 typename E3::Nested_t _C;
682
683 mutable gsVector<Scalar,3> onormal, normal, tangent, utangent, dtanu, dtanv, tvaru, tvarv, tvar2,
684 mu, mv, muv, mu_der, snvaru, snvarv, snvar2, nvar2;
685 mutable gsVector<Scalar> tmp;
686 mutable gsMatrix<Scalar> uGrads, vGrads, cJac, res, eC;
687
688 public:
689 enum{ Space = 3, ScalarValued= 0, ColBlocks= 0 };
690
ovar2_expr(const E1 & u,const E2 & v,const gsGeometryMap<Scalar> & G,_expr<E3> const & C)691 ovar2_expr(const E1 & u, const E2 & v, const gsGeometryMap<Scalar> & G, _expr<E3> const& C) : _u(u), _v(v), _G(G), _C(C) { }
692
693 EIGEN_MAKE_ALIGNED_OPERATOR_NEW
694
695 // helper function
vecFun(index_t pos,Scalar val)696 static inline gsVector<Scalar,3> vecFun(index_t pos, Scalar val)
697 {
698 gsVector<Scalar,3> result = gsVector<Scalar,3>::Zero();
699 result[pos] = val;
700 return result;
701 }
702
eval(const index_t k)703 const gsMatrix<Scalar> & eval(const index_t k) const
704 {
705 GISMO_ASSERT(_G.data().dim.second==3,"Domain dimension should be 3, is "<<_G.data().dim.second);
706 GISMO_ASSERT(_C.cols()*_C.rows()==3, "Size of vector is incorrect");
707 res.resize(_u.cardinality(), _u.cardinality());
708
709 eC = _C.eval(k);
710 eC.resize(1,3);
711
712 cJac = _G.data().values[1].reshapeCol(k, _G.data().dim.first, _G.data().dim.second).transpose();
713
714 onormal = _G.data().outNormal(k);
715 tmp = cJac.transpose() * onormal;
716 Scalar tol = 1e-8;
717
718 /*
719 We can check which column of the Jacobian corresponds to the outer normal vector or to the tangent.
720 The tangent is a covariant vector and hence the column of the Jacobian should be equal to the tangent.
721 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.
722 */
723 index_t colIndex = -1;
724 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
725 colIndex = 0;
726 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
727 colIndex = 1;
728 else // then the normal is unknown??
729 gsInfo<<"warning: choice unknown\n";
730
731 tangent = cJac.col(colIndex);
732 utangent = tangent / tangent.norm();
733
734 // Required for the normal vector variation
735 normal = _G.data().normal(k);
736 normal.normalize();
737 uGrads = _u.data().values[1].col(k);
738 vGrads = _v.data().values[1].col(k);
739
740 const index_t cardU = _u.data().values[0].rows(); // number of actives per component of u
741 const index_t cardV = _v.data().values[0].rows(); // number of actives per component of v
742 const Scalar measure = _G.data().measures.at(k);
743
744 for (index_t j = 0; j!= cardU; ++j) // for all basis functions u (1)
745 {
746 for (index_t i = 0; i!= cardV; ++i) // for all basis functions v (1)
747 {
748 for (index_t d = 0; d!= _u.dim(); ++d) // for all basis functions u (2)
749 {
750 const short_t s = d*cardU;
751
752 // first variation of the tangent (colvector)
753 dtanu = vecFun(d, uGrads.at(2*j+colIndex));
754 tvaru = 1 / tangent.norm() * ( dtanu - ( utangent.dot(dtanu) ) * utangent );
755
756 // first variation of the surface normal (colvector)
757 mu.noalias() = ( vecFun(d, uGrads.at(2*j ) ).cross( cJac.col(1).template head<3>() )
758 -vecFun(d, uGrads.at(2*j+1) ).cross( cJac.col(0).template head<3>() ))
759 / measure;
760 snvaru.noalias() = (mu - ( normal.dot(mu) ) * normal);
761
762 for (index_t c = 0; c!= _v.dim(); ++c) // for all basis functions v (2)
763 {
764 const short_t r = c*cardV;
765
766 // first variation of the tangent (colvector)
767 dtanv = vecFun(c, vGrads.at(2*i+colIndex));
768 tvarv = 1 / tangent.norm() * ( dtanv - ( utangent.dot(dtanv) ) * utangent );
769
770 // first variation of the surface normal (colvector)
771 mv.noalias() = ( vecFun(c, vGrads.at(2*i ) ).cross( cJac.col(1).template head<3>() )
772 -vecFun(c, vGrads.at(2*i+1) ).cross( cJac.col(0).template head<3>() ))
773 / measure;
774 snvarv.noalias() = (mv - ( normal.dot(mv) ) * normal);
775
776
777 // Second variation of the tangent (colvector)
778 // tvar2 = 1 / tangent.norm() * ( tvarv.dot(dtanu) * tangent )
779 // + 1 / (tangent.norm()*tangent.norm())
780 // * ( 2*( (tangent.dot(dtanu))*(tangent.dot(dtanv))*tangent )
781 // - ( tangent.dot(dtanu)*dtanv ) - ( tangent.dot(dtanv)*dtanu) );
782 tvar2 = -1 / tangent.norm() * (
783 ( tvarv.dot(dtanu) * utangent )
784 + ( utangent.dot(dtanv) * tvaru )
785 + ( utangent.dot(dtanu) * tvarv )
786 );
787
788 // Second variation of the surface normal (colvector)
789 muv.noalias() = ( vecFun(d, uGrads.at(2*j ) ).cross( vecFun(c, vGrads.at(2*i+1) ) )
790 +vecFun(c, vGrads.at(2*i ) ).cross( vecFun(d, uGrads.at(2*j+1) ) ))
791 / measure;
792
793 mu_der.noalias() = (muv - ( normal.dot(mv) ) * mu);
794
795 snvar2 = mu_der - (mu.dot(snvarv) + normal.dot(mu_der) ) * normal - (normal.dot(mu) ) * snvarv;
796
797 // Second variation of the outer normal (colvector)
798 nvar2 = tvar2.cross(normal) + tvaru.cross(snvarv) + tvarv.cross(snvaru) + utangent.cross(snvar2);
799
800 res(s + j, r + i ) = (eC * nvar2)(0,0);
801 }
802 }
803 }
804 }
805 return res;
806 }
807
rows()808 index_t rows() const { return 1; }
809
cols()810 index_t cols() const { return 1; }
811
parse(gsExprHelper<Scalar> & evList)812 void parse(gsExprHelper<Scalar> & evList) const
813 {
814 evList.add(_u);
815 _u.data().flags |= NEED_ACTIVE | NEED_VALUE | NEED_GRAD;
816 evList.add(_v);
817 _v.data().flags |= NEED_ACTIVE | NEED_VALUE | NEED_GRAD;
818 evList.add(_G);
819 _G.data().flags |= NEED_NORMAL | NEED_OUTER_NORMAL | NEED_DERIV | NEED_MEASURE;
820 _C.parse(evList);
821 }
822
rowVar()823 const gsFeSpace<Scalar> & rowVar() const { return _u.rowVar(); }
colVar()824 const gsFeSpace<Scalar> & colVar() const { return _v.rowVar(); }
825
print(std::ostream & os)826 void print(std::ostream &os) const { os << "nvar2("; _u.print(os); os <<")"; }
827 };
828
829 /**
830 * @brief Expression that takes the second derivative of an expression and multiplies it with a row vector
831 *
832 * @tparam E1 Expression
833 * @tparam E2 Row vector
834 */
835 template<class E1, class E2>
836 class deriv2dot_expr : public _expr<deriv2dot_expr<E1, E2> >
837 {
838 typename E1::Nested_t _u;
839 typename E2::Nested_t _v;
840
841 public:
842 enum{ Space = E1::Space, ScalarValued= 0, ColBlocks= 0 };
843
844 typedef typename E1::Scalar Scalar;
845
deriv2dot_expr(const E1 & u,const E2 & v)846 deriv2dot_expr(const E1 & u, const E2 & v) : _u(u), _v(v) { }
847
848 mutable gsMatrix<Scalar> res,tmp, vEv;
849
eval(const index_t k)850 const gsMatrix<Scalar> & eval(const index_t k) const {return eval_impl(_u,k); }
851
rows()852 index_t rows() const
853 {
854 return 1; //since the product with another vector is computed
855 }
856
cols()857 index_t cols() const
858 {
859 return cols_impl(_u);
860 }
861
parse(gsExprHelper<Scalar> & evList)862 void parse(gsExprHelper<Scalar> & evList) const
863 {
864 parse_impl<E1>(evList);
865 }
866
rowVar()867 const gsFeSpace<Scalar> & rowVar() const
868 {
869 // Note: what happens if E2 is a space? The following can fix it:
870 // if (E1::Space == 1 && E2::Space == 0)
871 // return _u.rowVar();
872 // else if (E1::Space == 0 && E2::Space == 1)
873 // return _v.rowVar();
874 // else
875
876 return _u.rowVar();
877 }
878
colVar()879 const gsFeSpace<Scalar> & colVar() const
880 {
881 // Note: what happens if E2 is a space? The following can fix it:
882 // if (E1::Space == 1 && E2::Space == 0)
883 // return _v.rowVar();
884 // else if (E1::Space == 0 && E2::Space == 1)
885 // return _u.rowVar();
886 // else
887
888 return _v.rowVar();
889 }
890
print(std::ostream & os)891 void print(std::ostream &os) const { os << "deriv2("; _u.print(os); _v.print(os); os <<")"; }
892
893 private:
894 template<class U> inline
895 typename util::enable_if< !util::is_same<U,gsFeSolution<Scalar> >::value,void>::type
parse_impl(gsExprHelper<Scalar> & evList)896 parse_impl(gsExprHelper<Scalar> & evList) const
897 {
898 evList.add(_u); // We manage the flags of _u "manually" here (sets data)
899 _u.data().flags |= NEED_DERIV2; // define flags
900
901 _v.parse(evList); // We need to evaluate _v (_v.eval(.) is called)
902
903 // Note: evList.parse(.) is called only in exprAssembler for the global expression
904 }
905
906 template<class U> inline
907 typename util::enable_if< util::is_same<U,gsFeSolution<Scalar> >::value,void>::type
parse_impl(gsExprHelper<Scalar> & evList)908 parse_impl(gsExprHelper<Scalar> & evList) const
909 {
910 _u.parse(evList); //
911 hess(_u).parse(evList); //
912
913 // evList.add(_u); // We manage the flags of _u "manually" here (sets data)
914 _v.parse(evList); // We need to evaluate _v (_v.eval(.) is called)
915 }
916
917 template<class U> inline
918 typename util::enable_if< util::is_same<U,gsGeometryMap<Scalar> >::value, const gsMatrix<Scalar> & >::type
eval_impl(const U & u,const index_t k)919 eval_impl(const U & u, const index_t k) const
920 {
921 /*
922 Here, we multiply the hessian of the geometry map by a vector, which possibly has multiple actives.
923 The hessian of the geometry map c has the form: hess(c)
924 [d11 c1, d11 c2, d11 c3]
925 [d22 c1, d22 c2, d22 c3]
926 [d12 c1, d12 c2, d12 c3]
927 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)
928 So we simply evaluate for every active basis function v_k the product hess(c).v_k
929 */
930
931
932 // evaluate the geometry map of U
933 tmp =_u.data().values[2].reshapeCol(k, cols(), _u.data().dim.second );
934 vEv = _v.eval(k);
935 res = vEv * tmp.transpose();
936 return res;
937 }
938
939 template<class U> inline
940 typename util::enable_if<util::is_same<U,gsFeSpace<Scalar> >::value, const gsMatrix<Scalar> & >::type
eval_impl(const U & u,const index_t k)941 eval_impl(const U & u, const index_t k) const
942 {
943 /*
944 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
945 This implies that hess(v) = [hess(v_1), hess(v_2), hess(v_3)] only has nonzero entries in column i. Hence,
946 hess(v) . normal = hess(v_i) * n_i (vector-scalar multiplication. The result is then of the form
947 [hess(v_1)*n_1 .., hess(v_2)*n_2 .., hess(v_3)*n_3 ..]. Here, the dots .. represent the active basis functions.
948 */
949 const index_t numAct = u.data().values[0].rows(); // number of actives of a basis function
950 const index_t cardinality = u.cardinality(); // total number of actives (=3*numAct)
951 res.resize(rows()*cardinality, cols() );
952 tmp.transpose() =_u.data().values[2].reshapeCol(k, cols(), numAct );
953 vEv = _v.eval(k);
954
955 for (index_t i = 0; i!=_u.dim(); i++)
956 res.block(i*numAct, 0, numAct, cols() ) = tmp * vEv.at(i);
957
958 return res;
959 }
960
961 template<class U> inline
962 typename util::enable_if<util::is_same<U,gsFeSolution<Scalar> >::value, const gsMatrix<Scalar> & >::type
eval_impl(const U & u,const index_t k)963 eval_impl(const U & u, const index_t k) const
964 {
965 /*
966 Here, we multiply the hessian of the geometry map by a vector, which possibly has multiple actives.
967 The hessian of the geometry map c has the form: hess(c)
968 [d11 c1, d22 c1, d12 c1]
969 [d11 c2, d22 c2, d12 c2]
970 [d11 c3, d22 c3, d12 c3]
971 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)
972 So we simply evaluate for every active basis function v_k the product hess(c).v_k
973 */
974
975 hess_expr<gsFeSolution<Scalar>> sHess = hess_expr<gsFeSolution<Scalar>>(_u); // NOTE: This does not parse automatically!
976 tmp = sHess.eval(k);
977 vEv = _v.eval(k);
978 res = vEv * tmp;
979 return res;
980 }
981
982 template<class U> inline
983 typename util::enable_if< util::is_same<U,gsGeometryMap<Scalar> >::value, index_t >::type
cols_impl(const U & u)984 cols_impl(const U & u) const
985 {
986 return _u.data().dim.second;
987 }
988
989 template<class U> inline
990 typename util::enable_if< !util::is_same<U,gsGeometryMap<Scalar> >::value, index_t >::type
cols_impl(const U & u)991 cols_impl(const U & u) const
992 {
993 return _u.dim();
994 }
995
996 // template<class U> inline
997 // typename util::enable_if<util::is_same<U,gsFeSolution<Scalar> >::value, index_t >::type
998 // cols_impl(const U & u) const
999 // {
1000 // return _u.dim();
1001 // }
1002
1003 };
1004
1005
1006
1007 /**
1008 * @brief Computes the second derivative of an expression
1009 *
1010 * It assumes that the vector of basis functions is of the form v = u*e_i where u
1011 * is the scalar basis function u: [0,1]^3 -> R^1 and e_i is the unit vector with a 1 on index i and a 0 elsewhere.
1012 * Let us define the following blocks
1013 * hess1(u) = hess2(u) = hess3(u) =
1014 * [d11 u , 0 , 0 ] | [0 , d11 u , 0 ] | [0 , 0 , d11 u ]
1015 * [d22 u , 0 , 0 ] | [0 , d22 u , 0 ] | [0 , 0 , d22 u ]
1016 * [d12 u , 0 , 0 ] | [0 , d12 u , 0 ] | [0 , 0 , d12 u ]
1017 *
1018 * Then the deriv2(u) is defined as follows (for k number of actives)
1019 * [hess1(u)_1]
1020 * ...
1021 * [hess1(u)_k]
1022 * [hess2(u)_1]
1023 * ...
1024 * [hess2(u)_k]
1025 * [hess3(u)_1]
1026 * ...
1027 * [hess3(u)_k]
1028 *
1029 *
1030 * @tparam E Expression type
1031 */
1032 template<class E>
1033 class deriv2_expr : public _expr<deriv2_expr<E> >
1034 {
1035 typename E::Nested_t _u;
1036
1037 public:
1038 // enum {ColBlocks = E::rowSpan }; // ????
1039 enum{ Space = E::Space, ScalarValued= 0, ColBlocks = (1==E::Space?1:0) };
1040
1041 typedef typename E::Scalar Scalar;
1042
deriv2_expr(const E & u)1043 deriv2_expr(const E & u) : _u(u) { }
1044
1045 mutable gsMatrix<Scalar> res, tmp;
1046
eval(const index_t k)1047 const gsMatrix<Scalar> & eval(const index_t k) const {return eval_impl(_u,k); }
1048
rows()1049 index_t rows() const //(components)
1050 {
1051 return 3; // _u.dim() for space or targetDim() for geometry
1052 }
1053
cols()1054 index_t cols() const
1055 {
1056 return _u.source().domainDim() * ( _u.source().domainDim() + 1 ) / 2;
1057 }
1058
parse(gsExprHelper<Scalar> & evList)1059 void parse(gsExprHelper<Scalar> & evList) const
1060 {
1061 _u.parse(evList);
1062 _u.data().flags |= NEED_DERIV2;
1063 }
1064
rowVar()1065 const gsFeSpace<Scalar> & rowVar() const { return _u.rowVar(); }
colVar()1066 const gsFeSpace<Scalar> & colVar() const { return _u.colVar(); }
1067
cardinality_impl()1068 index_t cardinality_impl() const { return _u.cardinality_impl(); }
1069
print(std::ostream & os)1070 void print(std::ostream &os) const { os << "deriv2("; _u.print(os); os <<")"; }
1071
1072 private:
1073 template<class U> inline
1074 typename util::enable_if< util::is_same<U,gsGeometryMap<Scalar> >::value, const gsMatrix<Scalar> & >::type
eval_impl(const U & u,const index_t k)1075 eval_impl(const U & u, const index_t k) const
1076 {
1077 /*
1078 Here, we compute the hessian of the geometry map.
1079 The hessian of the geometry map c has the form: hess(c)
1080 [d11 c1, d11 c2, d11 c3]
1081 [d22 c1, d22 c2, d22 c3]
1082 [d12 c1, d12 c2, d12 c3]
1083
1084 The geometry map has components c=[c1,c2,c3]
1085 */
1086 // evaluate the geometry map of U
1087 res = _u.data().values[2].reshapeCol(k, cols(), _u.data().dim.second );
1088 return res;
1089 }
1090
1091 template<class U> inline
1092 typename util::enable_if< util::is_same<U,gsFeVariable<Scalar> >::value, const gsMatrix<Scalar> & >::type
eval_impl(const U & u,const index_t k)1093 eval_impl(const U & u, const index_t k) const
1094 {
1095 /*
1096 Here, we compute the hessian of the geometry map.
1097 The hessian of the geometry map c has the form: hess(c)
1098 [d11 c1, d11 c2, d11 c3]
1099 [d22 c1, d22 c2, d22 c3]
1100 [d12 c1, d12 c2, d12 c3]
1101
1102 The geometry map has components c=[c1,c2,c3]
1103 */
1104 // evaluate the geometry map of U
1105 tmp = _u.data().values[2];
1106 res.resize(rows(),cols());
1107 for (index_t comp = 0; comp != _u.source().targetDim(); comp++)
1108 res.col(comp) = tmp.block(comp*rows(),0,rows(),1); //star,length
1109 return res;
1110 }
1111
1112 template<class U> inline
1113 typename util::enable_if< util::is_same<U,gsFeSolution<Scalar> >::value, const gsMatrix<Scalar> & >::type
eval_impl(const U & u,const index_t k)1114 eval_impl(const U & u, const index_t k) const
1115 {
1116 /*
1117 Here, we compute the hessian of the geometry map.
1118 The hessian of the geometry map c has the form: hess(c)
1119 [d11 c1, d11 c2, d11 c3]
1120 [d22 c1, d22 c2, d22 c3]
1121 [d12 c1, d12 c2, d12 c3]
1122
1123 The geometry map has components c=[c1,c2,c3]
1124 */
1125 // evaluate the geometry map of U
1126 hess_expr<gsFeSolution<Scalar>> sHess = hess_expr<gsFeSolution<Scalar>>(_u);
1127 res = sHess.eval(k).transpose();
1128 return res;
1129 }
1130
1131 /// Spexialization for a space
1132 template<class U> inline
1133 typename util::enable_if<util::is_same<U,gsFeSpace<Scalar> >::value, const gsMatrix<Scalar> & >::type
eval_impl(const U & u,const index_t k)1134 eval_impl(const U & u, const index_t k) const
1135 {
1136 /*
1137 Here, we compute the hessian of the basis with n actives.
1138 The hessian of the basis u has the form: hess(u)
1139 active 1 active 2 active n = cardinality
1140 [d11 u1, d11 u2, d11 u3] [d11 u1, d11 u2, d11 u3] ... [d11 u1, d11 u2, d11 u3]
1141 [d22 u1, d22 u2, d22 u3] [d22 u1, d22 u2, d22 u3] ... [d22 u1, d22 u2, d22 u3]
1142 [d12 u1, d12 u2, d12 u3] [d12 u1, d12 u2, d12 u3] ... [d12 u1, d12 u2, d12 u3]
1143
1144 Here, the basis function has components u = [u1,u2,u3]. Since they are evaluated for scalars
1145 we use blockDiag to make copies for all components ui
1146
1147 active 1 active 2 active k = cardinality/dim active 1 active 2k active 1 active 2k
1148 [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]
1149 [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]
1150 [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]
1151
1152 */
1153 const index_t numAct = u.data().values[0].rows(); // number of actives of a basis function
1154 const index_t cardinality = u.cardinality(); // total number of actives (=3*numAct)
1155
1156 res.resize(rows(), _u.dim() *_u.cardinality()); // (3 x 3*cardinality)
1157 res.setZero();
1158
1159 tmp = _u.data().values[2].reshapeCol(k, cols(), numAct );
1160 for (index_t d = 0; d != cols(); ++d)
1161 {
1162 const index_t s = d*(cardinality + 1);
1163 for (index_t i = 0; i != numAct; ++i)
1164 res.col(s+i*_u.cols()) = tmp.col(i);
1165 }
1166
1167 return res;
1168 }
1169
1170 };
1171
1172 /**
1173 * @brief Computes the product of expressions \a E1 and \a E2 and multiplies with a vector \a E3 in voight notation
1174 *
1175 * @tparam E1 Type of expression 1
1176 * @tparam E2 Type of expression 2
1177 * @tparam E3 Type of the vector expression
1178 */
1179 template<class E1, class E2, class E3>
1180 class flatdot_expr : public _expr<flatdot_expr<E1,E2,E3> >
1181 {
1182 public:
1183 typedef typename E1::Scalar Scalar;
1184
1185 private:
1186 typename E1::Nested_t _A;
1187 typename E2::Nested_t _B;
1188 typename E3::Nested_t _C;
1189 mutable gsMatrix<Scalar> eA, eB, eC, tmp, res;
1190
1191 public:
1192 enum {Space = E1::Space, ScalarValued = 0, ColBlocks = 0};
1193
1194 public:
1195
flatdot_expr(_expr<E1> const & A,_expr<E2> const & B,_expr<E3> const & C)1196 flatdot_expr(_expr<E1> const& A, _expr<E2> const& B, _expr<E3> const& C) : _A(A),_B(B),_C(C)
1197 {
1198 }
1199
eval(const index_t k)1200 const gsMatrix<Scalar> & eval(const index_t k) const
1201 {
1202 const index_t Ac = _A.cols();
1203 const index_t An = _A.cardinality();
1204 const index_t Bc = _B.cols();
1205 const index_t Bn = _B.cardinality();
1206
1207 eA = _A.eval(k);
1208 eB = _B.eval(k);
1209 eC = _C.eval(k);
1210
1211 GISMO_ASSERT(Bc==_A.rows(), "Dimensions: "<<Bc<<","<< _A.rows()<< "do not match");
1212 GISMO_STATIC_ASSERT(E1::Space==1, "First entry should be rowSpan" );
1213 GISMO_STATIC_ASSERT(E2::Space==2, "Second entry should be colSpan.");
1214
1215 res.resize(An, Bn);
1216
1217 // to do: evaluate the jacobians internally for the basis functions and loop over dimensions as well, since everything is same for all dimensions.
1218 for (index_t i = 0; i!=An; ++i) // for all actives u
1219 for (index_t j = 0; j!=Bn; ++j) // for all actives v
1220 {
1221 tmp.noalias() = eB.middleCols(i*Bc,Bc) * eA.middleCols(j*Ac,Ac);
1222
1223 tmp(0,0) *= eC.at(0);
1224 tmp(0,1) *= eC.at(2);
1225 tmp(1,0) *= eC.at(2);
1226 tmp(1,1) *= eC.at(1);
1227 res(i,j) = tmp.sum();
1228 }
1229 return res;
1230 }
1231
rows()1232 index_t rows() const { return 1; }
cols()1233 index_t cols() const { return 1; }
1234
parse(gsExprHelper<Scalar> & evList)1235 void parse(gsExprHelper<Scalar> & evList) const
1236 { _A.parse(evList);_B.parse(evList);_C.parse(evList); }
1237
rowVar()1238 const gsFeSpace<Scalar> & rowVar() const { return _A.rowVar(); }
colVar()1239 const gsFeSpace<Scalar> & colVar() const { return _B.colVar(); }
cardinality_impl()1240 index_t cardinality_impl() const { return _A.cardinality_impl(); }
1241
print(std::ostream & os)1242 void print(std::ostream &os) const { os << "flatdot("; _A.print(os);_B.print(os);_C.print(os); os<<")"; }
1243 };
1244
1245 /**
1246 * @brief Computes the product of expressions \a E1 and \a E2 and multiplies with a vector \a E3 in voight notation
1247 *
1248 * NOTE: SPECIALIZED FOR THE SECOND VARIATION OF THE CURVATURE!
1249 *
1250 * @tparam E1 Type of expression 1
1251 * @tparam E2 Type of expression 2
1252 * @tparam E3 Type of the vector expression
1253 */
1254 template<class E1, class E2, class E3>
1255 class flatdot2_expr : public _expr<flatdot2_expr<E1,E2,E3> >
1256 {
1257 public:
1258 typedef typename E1::Scalar Scalar;
1259
1260 private:
1261 typename E1::Nested_t _A;
1262 typename E2::Nested_t _B;
1263 typename E3::Nested_t _C;
1264 mutable gsMatrix<Scalar> eA, eB, eC, res, tmp;
1265
1266 public:
1267 enum {Space = E1::Space, ScalarValued = 0, ColBlocks = 0};
1268
flatdot2_expr(_expr<E1> const & A,_expr<E2> const & B,_expr<E3> const & C)1269 flatdot2_expr(_expr<E1> const& A, _expr<E2> const& B, _expr<E3> const& C) : _A(A),_B(B),_C(C)
1270 {
1271 }
1272
eval(const index_t k)1273 const gsMatrix<Scalar> & eval(const index_t k) const
1274 {
1275 const index_t Ac = _A.cols();
1276 const index_t An = _A.cardinality();
1277 const index_t Bn = _B.cardinality();
1278
1279 eA = _A.eval(k);
1280 eB = _B.eval(k);
1281 eC = _C.eval(k);
1282
1283 GISMO_ASSERT(_B.rows()==_A.cols(), "Dimensions: "<<_B.rows()<<","<< _A.cols()<< "do not match");
1284 GISMO_STATIC_ASSERT(E1::Space==1, "First entry should be rowSpan");
1285 GISMO_STATIC_ASSERT(E2::Space==2, "Second entry should be colSpan.");
1286 GISMO_ASSERT(_C.cols()==_B.rows(), "Dimensions: "<<_C.rows()<<","<< _B.rows()<< "do not match");
1287
1288 res.resize(An, Bn);
1289 for (index_t i = 0; i!=An; ++i)
1290 for (index_t j = 0; j!=Bn; ++j)
1291 {
1292 tmp = eA.middleCols(i*Ac,Ac) * eB.col(j); // E_f_der2
1293 tmp.row(2) *= 2.0; // multiply the third row of E_f_der2 by 2 for voight notation
1294 res(i,j) = (eC.row(0) * tmp.col(0)).value(); // E_f^T * mm * E_f_der2
1295 }
1296 return res;
1297 }
1298
rows()1299 index_t rows() const { return 1; }
cols()1300 index_t cols() const { return 1; }
1301
parse(gsExprHelper<Scalar> & evList)1302 void parse(gsExprHelper<Scalar> & evList) const
1303 { _A.parse(evList);_B.parse(evList);_C.parse(evList); }
1304
rowVar()1305 const gsFeSpace<Scalar> & rowVar() const { return _A.rowVar(); }
colVar()1306 const gsFeSpace<Scalar> & colVar() const { return _B.colVar(); }
cardinality_impl()1307 index_t cardinality_impl() const { return _A.cardinality_impl(); }
1308
print(std::ostream & os)1309 void print(std::ostream &os) const { os << "flatdot2("; _A.print(os);_B.print(os);_C.print(os); os<<")"; }
1310 };
1311
1312 /// Expression for the transformation matrix FROM local covariant TO local cartesian bases, based on a geometry map
1313 template<class T> class cartcovinv_expr ;
1314
1315 template<class T>
1316 class cartcov_expr : public _expr<cartcov_expr<T> >
1317 {
1318 typename gsGeometryMap<T>::Nested_t _G;
1319
1320 public:
1321 typedef T Scalar;
1322
1323 enum {Space = 0, ScalarValued = 0, ColBlocks = 0};
1324
cartcov_expr(const gsGeometryMap<T> & G)1325 cartcov_expr(const gsGeometryMap<T> & G) : _G(G) { }
1326
1327 mutable gsMatrix<Scalar> covBasis, conBasis, covMetric, conMetric, cartBasis;
1328 mutable gsMatrix<Scalar,3,3> result;
1329 mutable gsVector<Scalar> normal, tmp;
1330 mutable gsVector<Scalar> e1, e2, a1, a2;
1331
eval(const index_t k)1332 gsMatrix<Scalar> eval(const index_t k) const
1333 {
1334 if (_G.targetDim()==3)
1335 {
1336 // Compute covariant bases in deformed and undeformed configuration
1337 normal = _G.data().normals.col(k);
1338 normal.normalize();
1339 covBasis.resize(3,3);
1340 conBasis.resize(3,3);
1341 covBasis.leftCols(2) = _G.data().jacobian(k);
1342 covBasis.col(2) = normal;
1343 covMetric = covBasis.transpose() * covBasis;
1344
1345 conMetric = covMetric.inverse();
1346
1347 conBasis.col(1) = conMetric(1,0)*covBasis.col(0)+conMetric(1,1)*covBasis.col(1)+conMetric(1,2)*covBasis.col(2);
1348
1349 e1 = covBasis.col(0); e1.normalize();
1350 e2 = conBasis.col(1); e2.normalize();
1351
1352 a1 = covBasis.col(0);
1353 a2 = covBasis.col(1);
1354
1355 result(0,0) = (e1.dot(a1))*(a1.dot(e1));
1356 result(0,1) = (e1.dot(a2))*(a2.dot(e2));
1357 result(0,2) = 2*(e1.dot(a1))*(a2.dot(e1));
1358 // Row 1
1359 result(1,0) = (e2.dot(a1))*(a1.dot(e2));
1360 result(1,1) = (e2.dot(a2))*(a2.dot(e2));
1361 result(1,2) = 2*(e2.dot(a1))*(a2.dot(e2));
1362 // Row 2
1363 result(2,0) = (e1.dot(a1))*(a1.dot(e2));
1364 result(2,1) = (e1.dot(a2))*(a2.dot(e2));
1365 result(2,2) = (e1.dot(a1))*(a2.dot(e2)) + (e1.dot(a2))*(a1.dot(e2));
1366
1367 // return result.inverse(); // !!!!
1368 return result;
1369 }
1370 else if (_G.targetDim()==2)
1371 {
1372 // Compute covariant bases in deformed and undeformed configuration
1373 covBasis.resize(2,2);
1374 conBasis.resize(2,2);
1375 covBasis = _G.data().jacobian(k);
1376 covMetric = covBasis.transpose() * covBasis;
1377 conMetric = covMetric.inverse();
1378 conBasis.col(1) = conMetric(1,0)*covBasis.col(0)+conMetric(1,1)*covBasis.col(1);
1379
1380 e1 = covBasis.col(0); e1.normalize();
1381 e2 = conBasis.col(1); e2.normalize();
1382 // e3 = normal;
1383
1384 a1 = covBasis.col(0);
1385 a2 = covBasis.col(1);
1386
1387 result(0,0) = (e1.dot(a1))*(a1.dot(e1));
1388 result(0,1) = (e1.dot(a2))*(a2.dot(e2));
1389 result(0,2) = 2*(e1.dot(a1))*(a2.dot(e1));
1390 // Row 1
1391 result(1,0) = (e2.dot(a1))*(a1.dot(e2));
1392 result(1,1) = (e2.dot(a2))*(a2.dot(e2));
1393 result(1,2) = 2*(e2.dot(a1))*(a2.dot(e2));
1394 // Row 2
1395 result(2,0) = (e1.dot(a1))*(a1.dot(e2));
1396 result(2,1) = (e1.dot(a2))*(a2.dot(e2));
1397 result(2,2) = (e1.dot(a1))*(a2.dot(e2)) + (e1.dot(a2))*(a1.dot(e2));
1398
1399 return result;
1400 }
1401 else
1402 GISMO_ERROR("Not implemented");
1403 }
1404
inv()1405 cartcovinv_expr<T> inv() const
1406 {
1407 GISMO_ASSERT(rows() == cols(), "The Jacobian matrix is not square");
1408 return cartcovinv_expr<T>(_G);
1409 }
1410
rows()1411 index_t rows() const { return 3; }
1412
cols()1413 index_t cols() const { return 3; }
1414
rowVar()1415 static const gsFeSpace<Scalar> & rowVar() { return gsNullExpr<Scalar>::get(); }
colVar()1416 static const gsFeSpace<Scalar> & colVar() { return gsNullExpr<Scalar>::get(); }
1417
parse(gsExprHelper<Scalar> & evList)1418 void parse(gsExprHelper<Scalar> & evList) const
1419 {
1420 //GISMO_ASSERT(NULL!=m_fd, "FeVariable: FuncData member not registered");
1421 evList.add(_G);
1422 _G.data().flags |= NEED_NORMAL|NEED_DERIV;
1423 }
1424
print(std::ostream & os)1425 void print(std::ostream &os) const { os << "cartcov("; _G.print(os); os <<")"; }
1426 };
1427
1428 template<class T>
1429 class cartcovinv_expr : public _expr<cartcovinv_expr<T> >
1430 {
1431 typename gsGeometryMap<T>::Nested_t _G;
1432
1433 public:
1434 typedef T Scalar;
1435
1436 enum {Space = 0, ScalarValued = 0, ColBlocks = 0};
1437
cartcovinv_expr(const gsGeometryMap<T> & G)1438 cartcovinv_expr(const gsGeometryMap<T> & G) : _G(G) { }
1439
1440 mutable gsMatrix<T> temp;
1441
eval(const index_t k)1442 gsMatrix<T> eval(const index_t k) const
1443 {
1444 temp = (cartcov_expr<gsGeometryMap<T> >(_G).eval(k)).reshape(3,3).inverse();
1445 return temp;
1446 }
1447
rows()1448 index_t rows() const { return 3; }
1449
cols()1450 index_t cols() const { return 3; }
1451
rowVar()1452 static const gsFeSpace<Scalar> & rowVar() { return gsNullExpr<Scalar>::get(); }
colVar()1453 static const gsFeSpace<Scalar> & colVar() { return gsNullExpr<Scalar>::get(); }
1454
parse(gsExprHelper<Scalar> & evList)1455 void parse(gsExprHelper<Scalar> & evList) const
1456 {
1457 //GISMO_ASSERT(NULL!=m_fd, "FeVariable: FuncData member not registered");
1458 evList.add(_G);
1459 _G.data().flags |= NEED_NORMAL|NEED_DERIV;
1460 }
1461
print(std::ostream & os)1462 void print(std::ostream &os) const { os << "cartcovinv("; _G.print(os); os <<")"; }
1463 };
1464
1465
1466 /// Expression for the transformation matrix FROM local contravariant TO local cartesian bases, based on a geometry map
1467 template<class T> class cartconinv_expr ;
1468
1469 template<class T>
1470 class cartcon_expr : public _expr<cartcon_expr<T> >
1471 {
1472 typename gsGeometryMap<T>::Nested_t _G;
1473
1474 public:
1475 typedef T Scalar;
1476
1477 enum {Space = 0, ScalarValued = 0, ColBlocks = 0};
1478
cartcon_expr(const gsGeometryMap<T> & G)1479 cartcon_expr(const gsGeometryMap<T> & G) : _G(G) { }
1480
1481 mutable gsMatrix<Scalar> covBasis, conBasis, covMetric, conMetric, cartBasis;
1482 mutable gsMatrix<Scalar,3,3> result;
1483 mutable gsVector<Scalar> normal, tmp;
1484 mutable gsVector<Scalar> e1, e2, ac1, ac2;
1485
eval(const index_t k)1486 gsMatrix<Scalar> eval(const index_t k) const
1487 {
1488 if (_G.targetDim()==3)
1489 {
1490 // Compute covariant bases in deformed and undeformed configuration
1491 normal = _G.data().normals.col(k);
1492 normal.normalize();
1493 covBasis.resize(3,3);
1494 conBasis.resize(3,3);
1495 // Compute covariant bases in deformed and undeformed configuration
1496 normal = _G.data().normals.col(k);
1497 normal.normalize();
1498 covBasis.leftCols(2) = _G.data().jacobian(k);
1499 covBasis.col(2) = normal;
1500 covMetric = covBasis.transpose() * covBasis;
1501
1502 conMetric = covMetric.inverse();
1503
1504 conBasis.col(0) = conMetric(0,0)*covBasis.col(0)+conMetric(0,1)*covBasis.col(1)+conMetric(0,2)*covBasis.col(2);
1505 conBasis.col(1) = conMetric(1,0)*covBasis.col(0)+conMetric(1,1)*covBasis.col(1)+conMetric(1,2)*covBasis.col(2);
1506
1507 e1 = covBasis.col(0); e1.normalize();
1508 e2 = conBasis.col(1); e2.normalize();
1509 // e3 = normal;
1510
1511 ac1 = conBasis.col(0);
1512 ac2 = conBasis.col(1);
1513
1514 result(0,0) = (e1.dot(ac1))*(ac1.dot(e1));
1515 result(0,1) = (e1.dot(ac2))*(ac2.dot(e2));
1516 result(0,2) = 2*(e1.dot(ac1))*(ac2.dot(e1));
1517 // Row 1
1518 result(1,0) = (e2.dot(ac1))*(ac1.dot(e2));
1519 result(1,1) = (e2.dot(ac2))*(ac2.dot(e2));
1520 result(1,2) = 2*(e2.dot(ac1))*(ac2.dot(e2));
1521 // Row 2
1522 result(2,0) = (e1.dot(ac1))*(ac1.dot(e2));
1523 result(2,1) = (e1.dot(ac2))*(ac2.dot(e2));
1524 result(2,2) = (e1.dot(ac1))*(ac2.dot(e2)) + (e1.dot(ac2))*(ac1.dot(e2));
1525
1526 return result;
1527 }
1528 else if (_G.targetDim()==2)
1529 {
1530 // Compute covariant bases in deformed and undeformed configuration
1531 normal = _G.data().normals.col(k);
1532 normal.normalize();
1533 covBasis.resize(2,2);
1534 conBasis.resize(2,2);
1535 // Compute covariant bases in deformed and undeformed configuration
1536 covBasis = _G.data().jacobian(k);
1537 covMetric = covBasis.transpose() * covBasis;
1538
1539 conMetric = covMetric.inverse();
1540
1541 conBasis.col(0) = conMetric(0,0)*covBasis.col(0)+conMetric(0,1)*covBasis.col(1);
1542 conBasis.col(1) = conMetric(1,0)*covBasis.col(0)+conMetric(1,1)*covBasis.col(1);
1543
1544 e1 = covBasis.col(0); e1.normalize();
1545 e2 = conBasis.col(1); e2.normalize();
1546
1547 ac1 = conBasis.col(0);
1548 ac2 = conBasis.col(1);
1549
1550 result(0,0) = (e1.dot(ac1))*(ac1.dot(e1));
1551 result(0,1) = (e1.dot(ac2))*(ac2.dot(e2));
1552 result(0,2) = 2*(e1.dot(ac1))*(ac2.dot(e1));
1553 // Row 1
1554 result(1,0) = (e2.dot(ac1))*(ac1.dot(e2));
1555 result(1,1) = (e2.dot(ac2))*(ac2.dot(e2));
1556 result(1,2) = 2*(e2.dot(ac1))*(ac2.dot(e2));
1557 // Row 2
1558 result(2,0) = (e1.dot(ac1))*(ac1.dot(e2));
1559 result(2,1) = (e1.dot(ac2))*(ac2.dot(e2));
1560 result(2,2) = (e1.dot(ac1))*(ac2.dot(e2)) + (e1.dot(ac2))*(ac1.dot(e2));
1561
1562 return result;
1563 }
1564 else
1565 GISMO_ERROR("Not implemented");
1566
1567 }
1568
inv()1569 cartconinv_expr<T> inv() const
1570 {
1571 GISMO_ASSERT(rows() == cols(), "The Jacobian matrix is not square");
1572 return cartconinv_expr<T>(_G);
1573 }
1574
rows()1575 index_t rows() const { return 3; }
1576
cols()1577 index_t cols() const { return 3; }
1578
rowVar()1579 static const gsFeSpace<Scalar> & rowVar() { return gsNullExpr<Scalar>::get(); }
colVar()1580 static const gsFeSpace<Scalar> & colVar() { return gsNullExpr<Scalar>::get(); }
1581
parse(gsExprHelper<Scalar> & evList)1582 void parse(gsExprHelper<Scalar> & evList) const
1583 {
1584 //GISMO_ASSERT(NULL!=m_fd, "FeVariable: FuncData member not registered");
1585 evList.add(_G);
1586 _G.data().flags |= NEED_NORMAL|NEED_DERIV;
1587 }
1588
print(std::ostream & os)1589 void print(std::ostream &os) const { os << "cartcon("; _G.print(os); os <<")"; }
1590 };
1591
1592 template<class T>
1593 class cartconinv_expr : public _expr<cartconinv_expr<T> >
1594 {
1595 private:
1596 typename gsGeometryMap<T>::Nested_t _G;
1597 mutable gsMatrix<T> temp;
1598
1599 public:
1600 typedef T Scalar;
1601
cartconinv_expr(const gsGeometryMap<T> & G)1602 cartconinv_expr(const gsGeometryMap<T> & G) : _G(G) { }
1603
1604
eval(const index_t k)1605 gsMatrix<T> eval(const index_t k) const
1606 {
1607 temp = (cartcon_expr<gsGeometryMap<T> >(_G).eval(k)).reshape(3,3).inverse();
1608 return temp;
1609 }
1610
rows()1611 index_t rows() const { return 3; }
1612
cols()1613 index_t cols() const { return 3; }
1614
rowVar()1615 static const gsFeSpace<Scalar> & rowVar() { return gsNullExpr<Scalar>::get(); }
colVar()1616 static const gsFeSpace<Scalar> & colVar() { return gsNullExpr<Scalar>::get(); }
1617
parse(gsExprHelper<Scalar> & evList)1618 void parse(gsExprHelper<Scalar> & evList) const
1619 {
1620 //GISMO_ASSERT(NULL!=m_fd, "FeVariable: FuncData member not registered");
1621 evList.add(_G);
1622 _G.data().flags |= NEED_NORMAL|NEED_DERIV;
1623 }
1624
print(std::ostream & os)1625 void print(std::ostream &os) const { os << "cartconinv("; _G.print(os); os <<")"; }
1626 };
1627
1628 EIGEN_STRONG_INLINE
uv(const index_t index,const index_t dim)1629 unitVec_expr uv(const index_t index, const index_t dim) { return unitVec_expr(index,dim); }
1630
1631 template<class E> EIGEN_STRONG_INLINE
var1(const E & u,const gsGeometryMap<typename E::Scalar> & G)1632 var1_expr<E> var1(const E & u, const gsGeometryMap<typename E::Scalar> & G) { return var1_expr<E>(u, G); }
1633
1634 template<class E1, class E2, class E3> EIGEN_STRONG_INLINE
var2(const E1 & u,const E2 & v,const gsGeometryMap<typename E1::Scalar> & G,const E3 & Ef)1635 var2_expr<E1,E2,E3> var2(const E1 & u, const E2 & v, const gsGeometryMap<typename E1::Scalar> & G, const E3 & Ef)
1636 { return var2_expr<E1,E2,E3>(u,v, G, Ef); }
1637
1638 template<class E> EIGEN_STRONG_INLINE
tvar1(const E & u,const gsGeometryMap<typename E::Scalar> & G)1639 tvar1_expr<E> tvar1(const E & u, const gsGeometryMap<typename E::Scalar> & G) { return tvar1_expr<E>(u, G); }
1640
1641 template<class E> EIGEN_STRONG_INLINE
ovar1(const E & u,const gsGeometryMap<typename E::Scalar> & G)1642 ovar1_expr<E> ovar1(const E & u, const gsGeometryMap<typename E::Scalar> & G) { return ovar1_expr<E>(u, G); }
1643
1644 template<class E1, class E2, class E3> EIGEN_STRONG_INLINE
ovar2(const E1 & u,const E2 & v,const gsGeometryMap<typename E1::Scalar> & G,const E3 & C)1645 ovar2_expr<E1,E2,E3> ovar2(const E1 & u, const E2 & v, const gsGeometryMap<typename E1::Scalar> & G, const E3 & C)
1646 { return ovar2_expr<E1,E2,E3>(u,v, G, C); }
1647
1648 // template<class E1, class E2> EIGEN_STRONG_INLINE
1649 // hessdot_expr<E1,E2> hessdot(const E1 & u, const E2 & v) { return hessdot_expr<E1,E2>(u, v); }
1650
1651 template<class E> EIGEN_STRONG_INLINE
deriv2(const E & u)1652 deriv2_expr<E> deriv2(const E & u) { return deriv2_expr<E>(u); }
1653
1654 template<class E1, class E2> EIGEN_STRONG_INLINE
deriv2(const E1 & u,const E2 & v)1655 deriv2dot_expr<E1, E2> deriv2(const E1 & u, const E2 & v) { return deriv2dot_expr<E1, E2>(u,v); }
1656
1657 template<class E1, class E2, class E3> EIGEN_STRONG_INLINE
flatdot(const E1 & u,const E2 & v,const E3 & w)1658 flatdot_expr<E1,E2,E3> flatdot(const E1 & u, const E2 & v, const E3 & w)
1659 { return flatdot_expr<E1,E2,E3>(u, v, w); }
1660
1661 template<class E1, class E2, class E3> EIGEN_STRONG_INLINE
flatdot2(const E1 & u,const E2 & v,const E3 & w)1662 flatdot2_expr<E1,E2,E3> flatdot2(const E1 & u, const E2 & v, const E3 & w)
1663 { return flatdot2_expr<E1,E2,E3>(u, v, w); }
1664
1665 template<class E> EIGEN_STRONG_INLINE
cartcov(const gsGeometryMap<E> & G)1666 cartcov_expr<E> cartcov(const gsGeometryMap<E> & G) { return cartcov_expr<E>(G); }
1667
1668 template<class E> EIGEN_STRONG_INLINE
cartcon(const gsGeometryMap<E> & G)1669 cartcon_expr<E> cartcon(const gsGeometryMap<E> & G) { return cartcon_expr<E>(G); }
1670
1671 }
1672 }
1673