1 /** @file gsFunctionExpr.hpp
2
3 @brief Provides implementation of FunctionExpr class.
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): A. Mantzaflaris
12 */
13
14 #pragma once
15
16 #include <gsCore/gsLinearAlgebra.h>
17
18 /* ExprTk options */
19
20 //This define will enable printing of debug information to stdout during
21 //the compilation process.
22 //#define exprtk_enable_debugging
23
24 // This define will disable the ability for expressions to have comments.
25 // Expressions that have comments when parsed with a build that has this
26 // option, will result in a compilation failure.
27 #define exprtk_disable_comments
28
29 // This define will disable the loop-wise 'break' and 'continue'
30 // capabilities. Any expression that contains those keywords will result
31 // in a compilation failure.
32 #define exprtk_disable_break_continue
33
34 // This define will disable the short-circuit '&' (and) and '|' (or)
35 // operators
36 #define exprtk_disable_sc_andor
37
38 // This define will disable all enhanced features such as strength
39 // reduction and special function optimisations and expression specific
40 // type instantiations. This feature will reduce compilation times and
41 // binary sizes but will also result in massive performance degradation
42 // of expression evaluations.
43 #if !defined(NDEBUG) || defined(__MINGW32__)
44 #define exprtk_disable_enhanced_features
45 #endif
46
47 // This define will disable all string processing capabilities. Any
48 // expression that contains a string or string related syntax will result
49 // in a compilation failure.
50 #define exprtk_disable_string_capabilities
51
52 #define exprtk_disable_rtl_io_file
53 #define exprtk_disable_rtl_vecops
54
55 // The order in which header files are included is essential.
56 //
57 // It is important that all forward declaration file
58 // "exprtk_X_forward.hpp"" are included BEFORE the header file
59 // "exprtk.hpp" so that specializations for is_true(), is_false(),
60 // etcetera are not known for ALL types that should be supported. All
61 // adaptor files "exprtk_X_adaptor.hpp" have to be included AFTER the
62 // file "exprtk.hpp".
63
64 #if defined(GISMO_WITH_ADIFF)
65 #define DScalar gismo::ad::DScalar2<real_t,-1>
66 #include <exprtk_ad_forward.hpp>
67 #endif
68
69 #if defined(GISMO_WITH_MPFR)
70 #include <exprtk_mpfr_forward.hpp>
71 #endif
72
73 #if defined(GISMO_WITH_GMP)
74 #include <exprtk_gmp_forward.hpp>
75 #endif
76
77 #if defined(GISMO_WITH_CODIPACK)
78 #include <gsCoDiPack/exprtk_codi_rf_forward.hpp>
79 #include <gsCoDiPack/exprtk_codi_rr_forward.hpp>
80 #endif
81
82 #if defined(GISMO_WITH_UNUM)
83 #include <gsUnum/exprtk_unum_posit_forward.hpp>
84 #endif
85
86 #include <exprtk.hpp>
87
88 #if defined(GISMO_WITH_ADIFF)
89 #include <exprtk_ad_adaptor.hpp>
90 #endif
91
92 #if defined(GISMO_WITH_MPFR)
93 #include <exprtk_mpfr_adaptor.hpp>
94 #endif
95
96 #if defined(GISMO_WITH_GMP)
97 #include <exprtk_gmp_adaptor.hpp>
98 #endif
99
100 #if defined(GISMO_WITH_CODIPACK)
101 #include <gsCoDiPack/exprtk_codi_rf_adaptor.hpp>
102 #include <gsCoDiPack/exprtk_codi_rr_adaptor.hpp>
103 #endif
104
105 #if defined(GISMO_WITH_UNUM)
106 #include <gsUnum/exprtk_unum_posit_adaptor.hpp>
107 #endif
108
109
110 #include <gsIO/gsXml.h>
111
112 namespace
113 {
114
115 // addition of mixed derivative for expressions
116 // see https://en.wikipedia.org/wiki/Finite_difference_coefficient
117 template <typename T>
mixed_derivative(const exprtk::expression<T> & e,T & x,T & y,const double & h=0.00001)118 T mixed_derivative(const exprtk::expression<T>& e,
119 T& x, T& y,
120 const double& h = 0.00001)
121 {
122 T num = T(0.0), tmp;
123 T x_init = x;
124 T y_init = y;
125
126 x = x_init + T(2.0) * h;
127 y = y_init + T(2.0) * h;
128 num += e.value();
129 y = y_init - T(2.0) * h;
130 num -= e.value();
131 x = x_init - T(2.0) * h;
132 num += e.value();
133 y = y_init + T(2.0) * h;
134 num -= e.value();
135
136 x = x_init + h;
137 y = y_init + h;
138 tmp = e.value();
139 y = y_init - h;
140 tmp -= e.value();
141 x = x_init - h;
142 tmp += e.value();
143 y = y_init + h;
144 tmp -= e.value();
145 num += 64* tmp;
146
147 x = x_init + T(2.0) * h;
148 y = y_init - h;
149 tmp = e.value();
150 y = y_init + h;
151 tmp -= e.value();
152 x = x_init - T(2.0) * h;
153 tmp += e.value();
154 y = y_init - h;
155 tmp -= e.value();
156
157 y = y_init + T(2.0) * h;
158 x = x_init - h;
159 tmp += e.value();
160 x = x_init + h;
161 tmp -= e.value();
162 y = y_init - T(2.0) * h;
163 tmp += e.value();
164 x = x_init - h;
165 tmp -= e.value();
166 num += 8* tmp;
167
168 x = x_init;
169 y = y_init;
170 return num / ( T(144.0)*h*h );
171 }
172
173 } //namespace
174
175 #define N_VARS 7
176
177 namespace gismo
178 {
179
180 template<typename T> class gsFunctionExpr<T>::gsFunctionExprPrivate
181 {
182 public:
183
184 #ifdef GISMO_WITH_ADIFF
185 typedef DScalar Numeric_t;
186 #else
187 typedef T Numeric_t;
188 #endif
189
190 typedef exprtk::symbol_table<Numeric_t> SymbolTable_t;
191 typedef exprtk::expression<Numeric_t> Expression_t;
192 typedef exprtk::parser<Numeric_t> Parser_t;
193
194 public:
195
gsFunctionExprPrivate(const short_t _dim)196 gsFunctionExprPrivate(const short_t _dim)
197 : vars(), dim(_dim)
198 {
199 GISMO_ENSURE( dim <= N_VARS, "The number of variables can be at most 7 (x,y,z,w,u,v,t)." );
200 init();
201 }
202
gsFunctionExprPrivate(const gsFunctionExprPrivate & other)203 gsFunctionExprPrivate(const gsFunctionExprPrivate & other)
204 : vars(), dim(other.dim)
205 {
206 GISMO_ASSERT ( string.size() == expression.size(), "Corrupted FunctionExpr");
207 init();
208 //copy_n(other.vars, N_VARS+1, vars);
209 string .reserve(string.size());
210 expression.reserve(string.size());
211 for (size_t i = 0; i!= other.string.size(); ++i)
212 addComponent(other.string[i]);
213 }
214
addComponent(const std::string & strExpression)215 void addComponent(const std::string & strExpression)
216 {
217 string.push_back( strExpression );// Keep string data
218 std::string & str = string.back();
219 str.erase(std::remove(str.begin(), str.end(),' '), str.end() );
220 gismo::util::string_replace(str, "**", "^");
221
222 // String expression
223 expression.push_back(Expression_t());
224 Expression_t & expr = expression.back();
225 //expr.release();
226 expr.register_symbol_table(symbol_table);
227
228 // Parser
229 Parser_t parser;
230 //Collect variable symbols
231 //parser.dec().collect_variables() = true;
232 bool success = parser.compile(str, expr);
233 if ( ! success )
234 gsWarn<<"gsFunctionExpr error: " <<parser.error() <<" while parsing "<<str<<"\n";
235 /*
236 typedef typename exprtk::parser_t::
237 dependent_entity_collector::symbol_t symbol_t;
238
239 std::deque<symbol_t> symbol_list;
240 parser.dec().symbols(symbol_list);
241 for (size_t i = 0; i != symbol_list.size(); ++i)
242 {
243 symbol_t& symbol = symbol_list[i];
244 // do something
245 }
246 //*/
247 }
248
init()249 void init()
250 {
251 //symbol_table.clear();
252 // Identify symbol table
253 symbol_table.add_variable("x",vars[0]);
254 symbol_table.add_variable("y",vars[1]);
255 symbol_table.add_variable("z",vars[2]);
256 symbol_table.add_variable("w",vars[3]);
257 symbol_table.add_variable("u",vars[4]);
258 symbol_table.add_variable("v",vars[5]);
259 symbol_table.add_variable("t",vars[6]);
260 //symbol_table.remove_variable("w",vars[3]);
261 symbol_table.add_pi();
262 //symbol_table.add_constant("C", 1);
263 }
264
265 public:
266 mutable Numeric_t vars[N_VARS];
267 SymbolTable_t symbol_table;
268 std::vector<Expression_t> expression;
269 std::vector<std::string> string;
270 short_t dim;
271
272 private:
273 gsFunctionExprPrivate();
274 gsFunctionExprPrivate operator= (const gsFunctionExprPrivate & other);
275 };
276
277 /* / /AM: under construction
278 template<typename T> class gsSymbolListPrivate
279 {
280 public:
281 exprtk::symbol_table<T> symbol_table;
282 std::vector<T> vars ;
283 std::vector<T> params;
284 };
285
286 template<typename T>
287 gsSymbolList<T>::gsSymbolList() : my(new gsSymbolListPrivate<T>) {}
288
289 template<typename T>
290 gsSymbolList<T>::~gsSymbolList()
291 {
292 delete my;
293 }
294
295 template<typename T>
296 void gsSymbolList<T>::setDefault()
297 {
298 my->symbol_table.clear();
299
300 // Identify symbol table
301 my->symbol_table.add_variable("x",my->vars[0]);
302 my->symbol_table.add_variable("y",my->vars[1]);
303 my->symbol_table.add_variable("z",my->vars[2]);
304 my->symbol_table.add_variable("w",my->vars[3]);
305 my->symbol_table.add_variable("u",my->vars[4]);
306 my->symbol_table.add_variable("v",my->vars[5]);
307 //my->symbol_table.remove_variable("w",my->vars[3]);
308 my->symbol_table.add_pi();
309 //my->symbol_table.add_constant("C", 1);
310 }
311
312 template<typename T>
313 bool gsSymbolList<T>::addConstant(const std::string & constant_name, const T& value)
314 {
315 return my->symbol_table.add_constant(constant_name, value);
316 }
317
318 template<typename T>
319 bool gsSymbolList<T>::addVariable(const std::string & variable_name)
320 {
321 my->vars.push_back(0.0);
322 return my->symbol_table.add_variable(variable_name, my->vars.back() );
323 }
324
325 template<typename T>
326 bool gsSymbolList<T>::addParameter(const std::string & variable_name)
327 {
328 my->params.push_back(0.0);
329 return my->symbol_table.add_variable(variable_name, my->params.back() );
330 }
331
332 template<typename T>
333 bool gsSymbolList<T>::hasSymbol(const std::string& symbol_name)
334 {
335 return true;
336 }
337 //*/
338
339 template<typename T>
gsFunctionExpr()340 gsFunctionExpr<T>::gsFunctionExpr() : my(new PrivateData_t(0))
341 { }
342
343 template<typename T>
gsFunctionExpr(const std::string & expression_string,short_t ddim)344 gsFunctionExpr<T>::gsFunctionExpr(const std::string & expression_string, short_t ddim)
345 : my(new PrivateData_t(ddim))
346 {
347 my->addComponent(expression_string);
348 }
349
350 template<typename T>
gsFunctionExpr(const std::string & expression_string1,const std::string & expression_string2,short_t ddim)351 gsFunctionExpr<T>::gsFunctionExpr(const std::string & expression_string1,
352 const std::string & expression_string2,
353 short_t ddim)
354 : my(new PrivateData_t(ddim))
355 {
356 my->addComponent(expression_string1);
357 my->addComponent(expression_string2);
358 }
359
360 template<typename T>
gsFunctionExpr(const std::string & expression_string1,const std::string & expression_string2,const std::string & expression_string3,short_t ddim)361 gsFunctionExpr<T>::gsFunctionExpr(const std::string & expression_string1,
362 const std::string & expression_string2,
363 const std::string & expression_string3,
364 short_t ddim)
365 : my(new PrivateData_t(ddim))
366 {
367 my->addComponent(expression_string1);
368 my->addComponent(expression_string2);
369 my->addComponent(expression_string3);
370 }
371
372 template<typename T>
gsFunctionExpr(const std::string & expression_string1,const std::string & expression_string2,const std::string & expression_string3,const std::string & expression_string4,short_t ddim)373 gsFunctionExpr<T>::gsFunctionExpr(const std::string & expression_string1,
374 const std::string & expression_string2,
375 const std::string & expression_string3,
376 const std::string & expression_string4,
377 short_t ddim)
378 : my(new PrivateData_t(ddim))
379 {
380 my->addComponent(expression_string1);
381 my->addComponent(expression_string2);
382 my->addComponent(expression_string3);
383 my->addComponent(expression_string4);
384 }
385
386 template<typename T>
gsFunctionExpr(const std::string & expression_string1,const std::string & expression_string2,const std::string & expression_string3,const std::string & expression_string4,const std::string & expression_string5,const std::string & expression_string6,const std::string & expression_string7,const std::string & expression_string8,const std::string & expression_string9,short_t ddim)387 gsFunctionExpr<T>::gsFunctionExpr(const std::string & expression_string1,
388 const std::string & expression_string2,
389 const std::string & expression_string3,
390 const std::string & expression_string4,
391 const std::string & expression_string5,
392 const std::string & expression_string6,
393 const std::string & expression_string7,
394 const std::string & expression_string8,
395 const std::string & expression_string9,
396 short_t ddim)
397 : my(new PrivateData_t(ddim))
398 {
399 my->addComponent(expression_string1);
400 my->addComponent(expression_string2);
401 my->addComponent(expression_string3);
402 my->addComponent(expression_string4);
403 my->addComponent(expression_string5);
404 my->addComponent(expression_string6);
405 my->addComponent(expression_string7);
406 my->addComponent(expression_string8);
407 my->addComponent(expression_string9);
408 }
409
410 template<typename T>
gsFunctionExpr(const std::vector<std::string> & expression_string,short_t ddim)411 gsFunctionExpr<T>::gsFunctionExpr(const std::vector<std::string> & expression_string,
412 short_t ddim)
413 : my(new PrivateData_t(ddim))
414 {
415 for (size_t i = 0; i!= expression_string.size(); ++i)
416 my->addComponent(expression_string[i]);
417 }
418
419 template<typename T>
gsFunctionExpr(const gsFunctionExpr & other)420 gsFunctionExpr<T>::gsFunctionExpr(const gsFunctionExpr & other)
421 {
422 my = new PrivateData_t(*other.my);
423 }
424 #if EIGEN_HAS_RVALUE_REFERENCES
425
426 template<typename T>
gsFunctionExpr(gsFunctionExpr && other)427 gsFunctionExpr<T>::gsFunctionExpr(gsFunctionExpr && other)
428 {
429 my = other.my; other.my = NULL;
430 }
431
432 template<typename T>
operator =(const gsFunctionExpr & other)433 gsFunctionExpr<T> & gsFunctionExpr<T>::operator=(const gsFunctionExpr& other)
434 {
435 if (this != &other)
436 {
437 delete my;
438 my = new PrivateData_t(*other.my);
439 }
440 return *this;
441 }
442
443 template<typename T>
operator =(gsFunctionExpr && other)444 gsFunctionExpr<T> & gsFunctionExpr<T>::operator=(gsFunctionExpr&& other)
445 {
446 if (this != &other)
447 {
448 delete my;
449 my = other.my; other.my = NULL;
450 }
451 return *this;
452 }
453
454 #else
455 template<typename T>
operator =(gsFunctionExpr other)456 gsFunctionExpr<T> & gsFunctionExpr<T>::operator=(gsFunctionExpr other)
457 {
458 std::swap(my,other.my);
459 return *this;
460 }
461 #endif
462
463 template<typename T>
~gsFunctionExpr()464 gsFunctionExpr<T>::~gsFunctionExpr()
465 {
466 delete my;
467 }
468
469 template<typename T>
domainDim() const470 short_t gsFunctionExpr<T>::domainDim() const
471 {
472 return my->dim;
473 }
474
475 template<typename T>
targetDim() const476 short_t gsFunctionExpr<T>::targetDim() const
477 {
478 return static_cast<short_t>(my->string.size());
479 }
480
481 template<typename T>
addComponent(const std::string & strExpression)482 void gsFunctionExpr<T>::addComponent(const std::string & strExpression)
483 {
484 my->addComponent(strExpression);
485 }
486
487 template<typename T>
expression(int i) const488 const std::string & gsFunctionExpr<T>::expression(int i) const
489 {
490 return my->string[i];
491 }
492
493 template<typename T>
set_x(T const & v) const494 void gsFunctionExpr<T>::set_x (T const & v) const { my->vars[0]= v; }
495
496 template<typename T>
set_y(T const & v) const497 void gsFunctionExpr<T>::set_y (T const & v) const { my->vars[1]= v; }
498
499 template<typename T>
set_z(T const & v) const500 void gsFunctionExpr<T>::set_z (T const & v) const { my->vars[2]= v; }
501
502 template<typename T>
set_w(T const & v) const503 void gsFunctionExpr<T>::set_w (T const & v) const { my->vars[3]= v; }
504
505 template<typename T>
set_u(T const & v) const506 void gsFunctionExpr<T>::set_u (T const & v) const { my->vars[4]= v; }
507
508 template<typename T>
set_v(T const & v) const509 void gsFunctionExpr<T>::set_v (T const & v) const { my->vars[5]= v; }
510
511 template<typename T>
set_t(T const & t) const512 void gsFunctionExpr<T>::set_t (T const & t) const { my->vars[6]= t; }
513
514 template<typename T>
eval_into(const gsMatrix<T> & u,gsMatrix<T> & result) const515 void gsFunctionExpr<T>::eval_into(const gsMatrix<T>& u, gsMatrix<T>& result) const
516 {
517 GISMO_ASSERT ( u.rows() == my->dim, "Inconsistent point dimension (expected: "
518 << my->dim <<", got "<< u.rows() <<")\n"<< *this);
519
520 const short_t n = targetDim();
521 result.resize(n, u.cols());
522
523 const PrivateData_t & expr =
524 # ifdef _OPENMP
525 omp_in_parallel() ? PrivateData_t(*my) :
526 # endif
527 *my;
528
529 for ( index_t p = 0; p!=u.cols(); p++ ) // for all evaluation points
530 {
531 copy_n(u.col(p).data(), expr.dim, expr.vars);
532
533 for (short_t c = 0; c!= n; ++c) // for all components
534 # ifdef GISMO_WITH_ADIFF
535 result(c,p) = expr.expression[c].value().getValue();
536 # else
537 result(c,p) = expr.expression[c].value();
538 # endif
539 }
540 }
541
542 template<typename T>
eval_component_into(const gsMatrix<T> & u,const index_t comp,gsMatrix<T> & result) const543 void gsFunctionExpr<T>::eval_component_into(const gsMatrix<T>& u, const index_t comp, gsMatrix<T>& result) const
544 {
545 GISMO_ASSERT ( u.rows() == my->dim, "Inconsistent point dimension (expected: "
546 << my->dim <<", got "<< u.rows() <<")");
547
548 GISMO_ASSERT (comp < targetDim(),
549 "Given component number is higher then number of components");
550
551 result.resize(1, u.cols());
552 for ( index_t p = 0; p!=u.cols(); ++p )
553 {
554 copy_n(u.col(p).data(), my->dim, my->vars);
555
556 # ifdef GISMO_WITH_ADIFF
557 result(0,p) = my->expression[comp].value().getValue();
558 # else
559 result(0,p) = my->expression[comp].value();
560 # endif
561 }
562 }
563
564 template<typename T>
deriv_into(const gsMatrix<T> & u,gsMatrix<T> & result) const565 void gsFunctionExpr<T>::deriv_into(const gsMatrix<T>& u, gsMatrix<T>& result) const
566 {
567 //gsDebug<< "Using finite differences (gsFunctionExpr::deriv_into) for derivatives.\n";
568 const short_t d = domainDim();
569 GISMO_ASSERT ( u.rows() == my->dim, "Inconsistent point dimension (expected: "
570 << my->dim <<", got "<< u.rows() <<")");
571
572 const short_t n = targetDim();
573 result.resize(d*n, u.cols());
574
575 const PrivateData_t & expr =
576 # ifdef _OPENMP
577 omp_in_parallel() ? PrivateData_t(*my) :
578 # endif
579 *my;
580
581 for ( index_t p = 0; p!=u.cols(); p++ ) // for all evaluation points
582 {
583 # ifdef GISMO_WITH_ADIFF
584 for (short_t k = 0; k!=d; ++k)
585 expr.vars[k].setVariable(k,d,u(k,p));
586 for (short_t c = 0; c!= n; ++c) // for all components
587 expr.expression[c].value().gradient_into(result.block(c*d,p,d,1));
588 //result.block(c*d,p,d,1) = expr.expression[c].value().getGradient(); //fails on constants
589 # else
590 copy_n(u.col(p).data(), expr.dim, expr.vars);
591 for (short_t c = 0; c!= n; ++c) // for all components
592 for ( short_t j = 0; j!=d; j++ ) // for all variables
593 result(c*d + j, p) =
594 exprtk::derivative<T>(expr.expression[c], expr.vars[j], 0.00001 ) ;
595 # endif
596 }
597 }
598
599 template<typename T>
deriv2_into(const gsMatrix<T> & u,gsMatrix<T> & result) const600 void gsFunctionExpr<T>::deriv2_into(const gsMatrix<T>& u, gsMatrix<T>& result) const
601 {
602 const short_t d = domainDim();
603 GISMO_ASSERT ( u.rows() == my->dim, "Inconsistent point dimension (expected: "
604 << my->dim <<", got "<< u.rows() <<")");
605
606 const short_t n = targetDim();
607 const index_t stride = d + d*(d-1)/2;
608 result.resize(stride*n, u.cols() );
609
610 const PrivateData_t & expr =
611 # ifdef _OPENMP
612 omp_in_parallel() ? PrivateData_t(*my) :
613 # endif
614 *my;
615
616 for ( index_t p = 0; p!=u.cols(); p++ ) // for all evaluation points
617 {
618 # ifndef GISMO_WITH_ADIFF
619 copy_n(u.col(p).data(), expr.dim, expr.vars);
620 # endif
621
622 for (short_t c = 0; c!= n; ++c) // for all components
623 {
624 # ifdef GISMO_WITH_ADIFF
625 for (index_t v = 0; v!=d; ++v)
626 expr.vars[v].setVariable(v,d,u(v,p));
627 const DScalar & ads = expr.expression[c].value();
628 const DScalar::Hessian_t & Hmat = ads.getHessian(); // note: can fail
629
630 for ( index_t k=0; k!=d; ++k)
631 {
632 result(k,p) = Hmat(k,k);
633 index_t m = d;
634 for ( index_t l=k+1; l<d; ++l)
635 result(m++,p) = Hmat(k,l);
636 }
637 # else
638 for (short_t k = 0; k!=d; ++k)
639 {
640 // H_{k,k}
641 result(k,p) = exprtk::
642 second_derivative<T>(expr.expression[c], expr.vars[k], 0.00001);
643
644 short_t m = d;
645 for (short_t l=k+1; l<d; ++l)
646 {
647 // H_{k,l}
648 result(m++,p) =
649 mixed_derivative<T>( expr.expression[c], expr.vars[k],
650 expr.vars[l], 0.00001 );
651 }
652 }
653 # endif
654 }
655 }
656 }
657
658 template<typename T>
659 gsMatrix<T>
hess(const gsMatrix<T> & u,unsigned coord) const660 gsFunctionExpr<T>::hess(const gsMatrix<T>& u, unsigned coord) const
661 {
662 //gsDebug<< "Using finite differences (gsFunctionExpr::hess) for Hessian.\n";
663 GISMO_ENSURE(coord == 0, "Error, function is real");
664 GISMO_ASSERT ( u.cols() == 1, "Need a single evaluation point." );
665 const index_t d = u.rows();
666 GISMO_ASSERT ( u.rows() == my->dim, "Inconsistent point dimension (expected: "
667 << my->dim <<", got "<< u.rows() <<")");
668
669 gsMatrix<T> res(d, d);
670
671 const PrivateData_t & expr =
672 # ifdef _OPENMP
673 omp_in_parallel() ? PrivateData_t(*my) :
674 # endif
675 *my;
676
677 # ifdef GISMO_WITH_ADIFF
678 for (index_t v = 0; v!=d; ++v)
679 expr.vars[v].setVariable(v, d, u(v,0) );
680 expr.expression[coord].value().hessian_into(res);
681 # else
682 copy_n(u.data(), expr.dim, expr.vars);
683 for( index_t j=0; j!=d; ++j )
684 {
685 res(j,j) = exprtk::
686 second_derivative<T>( expr.expression[coord], expr.vars[j], 0.00001);
687
688 for( index_t k = 0; k!=j; ++k )
689 res(k,j) = res(j,k) =
690 mixed_derivative<T>( expr.expression[coord], expr.vars[k],
691 expr.vars[j], 0.00001 );
692 }
693 # endif
694
695 return res;
696 }
697
698 template<typename T>
mderiv(const gsMatrix<T> & u,const index_t k,const index_t j) const699 gsMatrix<T> * gsFunctionExpr<T>::mderiv(const gsMatrix<T> & u,
700 const index_t k,
701 const index_t j) const
702 {
703 GISMO_ASSERT ( u.rows() == my->dim, "Inconsistent point size.");
704 const short_t n = targetDim();
705 gsMatrix<T> * res= new gsMatrix<T>(n,u.cols()) ;
706
707 const PrivateData_t & expr =
708 # ifdef _OPENMP
709 omp_in_parallel() ? PrivateData_t(*my) :
710 # endif
711 *my;
712
713 for( index_t p=0; p!=res->cols(); ++p )
714 {
715 # ifndef GISMO_WITH_ADIFF
716 copy_n(u.col(p).data(), expr.dim, expr.vars);
717 # endif
718
719 for (short_t c = 0; c!= n; ++c) // for all components
720 {
721 # ifdef GISMO_WITH_ADIFF
722 for (index_t v = 0; v!=expr.dim; ++v)
723 expr.vars[v].setVariable(v, expr.dim, u(v,p) );
724 (*res)(c,p) = expr.expression[c].value().getHessian()(k,j); //note: can fail
725 # else
726 (*res)(c,p) =
727 mixed_derivative<T>( expr.expression[c], expr.vars[k], expr.vars[j], 0.00001 ) ;
728 # endif
729 }
730 }
731 return res;
732 }
733
734 template<typename T>
laplacian(const gsMatrix<T> & u) const735 gsMatrix<T> gsFunctionExpr<T>::laplacian(const gsMatrix<T>& u) const
736 {
737 //gsDebug<< "Using finite differences (gsFunction::laplacian) for Laplacian.\n";
738 GISMO_ASSERT ( u.rows() == my->dim, "Inconsistent point size.");
739 const short_t n = targetDim();
740 gsMatrix<T> res(n,u.cols());
741
742 const PrivateData_t & expr =
743 # ifdef _OPENMP
744 omp_in_parallel() ? PrivateData_t(*my) :
745 # endif
746 *my;
747
748 for( index_t p = 0; p != res.cols(); ++p )
749 {
750 # ifndef GISMO_WITH_ADIFF
751 copy_n(u.col(p).data(), expr.dim, expr.vars);
752 # endif
753
754 for (short_t c = 0; c!= n; ++c) // for all components
755 {
756 # ifdef GISMO_WITH_ADIFF
757 for (index_t v = 0; v!=expr.dim; ++v)
758 expr.vars[v].setVariable(v, expr.dim, u(v,p) );
759 res(c,p) = expr.expression[c].value().getHessian().trace();
760 # else
761 T & val = res(c,p);
762 for ( index_t j = 0; j!=expr.dim; ++j )
763 val += exprtk::
764 second_derivative<T>( expr.expression[c], expr.vars[j], 0.00001 );
765 # endif
766 }
767 }
768 return res;
769 }
770
771 template<typename T>
print(std::ostream & os) const772 std::ostream & gsFunctionExpr<T>::print(std::ostream &os) const
773 {
774 os <<"[ ";
775 if( my->string.empty() )
776 os << "empty";
777 else
778 {
779 os << my->string[0];
780
781 for (short_t k = 1; k<targetDim(); ++k)
782 os <<", " << my->string[k];
783 }
784 os <<" ]";
785 return os;
786 }
787
788 namespace internal
789 {
790
791 /// @brief Get a FunctionsExpr from XML data
792 template<class T>
793 class gsXml< gsFunctionExpr<T> >
794 {
795 private:
gsXml()796 gsXml() { }
797 typedef gsFunctionExpr<T> Object;
798 public:
799 GSXML_COMMON_FUNCTIONS(Object);
tag()800 static std::string tag () { return "Function"; }
type()801 static std::string type () { return "FunctionExpr"; }
802
803 GSXML_GET_POINTER(Object);
804
get_into(gsXmlNode * node,Object & obj)805 static void get_into (gsXmlNode * node, Object & obj)
806 {
807 getFunctionFromXml(node, obj);
808 }
809
put(const Object & obj,gsXmlTree & data)810 static gsXmlNode * put (const Object & obj,
811 gsXmlTree & data )
812 {
813 gsXmlNode * func = makeNode("FunctionExpr", data);
814 func->append_attribute(makeAttribute("dim", obj.domainDim(), data));
815
816 const short_t tdim = obj.targetDim();
817
818 if ( tdim == 1)
819 {
820 func->value( makeValue(obj.expression(), data) );
821 }
822 else
823 {
824 gsXmlNode * cnode;
825 for (short_t c = 0; c!=tdim; ++c)
826 {
827 cnode = makeNode("c", obj.expression(c), data);
828 func->append_node(cnode);
829 }
830 }
831
832 return func;
833 }
834 };
835
836 } // internal
837
838 }; // namespace gismo
839