GetFEM++  5.3
getfem_generic_assembly_functions_and_operators.cc
1 /*===========================================================================
2 
3  Copyright (C) 2013-2018 Yves Renard
4 
5  This file is a part of GetFEM++
6 
7  GetFEM++ is free software; you can redistribute it and/or modify it
8  under the terms of the GNU Lesser General Public License as published
9  by the Free Software Foundation; either version 3 of the License, or
10  (at your option) any later version along with the GCC Runtime Library
11  Exception either version 3.1 or (at your option) any later version.
12  This program is distributed in the hope that it will be useful, but
13  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
14  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
15  License and GCC Runtime Library Exception for more details.
16  You should have received a copy of the GNU Lesser General Public License
17  along with this program; if not, write to the Free Software Foundation,
18  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
19 
20 ===========================================================================*/
21 
22 #include "getfem/getfem_generic_assembly_semantic.h"
23 #include "getfem/getfem_generic_assembly_functions_and_operators.h"
24 #include "getfem/getfem_generic_assembly_compile_and_exec.h"
25 
26 /**
27  Providing for special Math functions unavailable on Intel or MSVS C++
28  compilers
29 */
30 
31 #if defined(_MSC_VER) && _MSC_VER < 1800
32 #include <boost/math/special_functions/acosh.hpp>
33 #include <boost/math/special_functions/asinh.hpp>
34 #include <boost/math/special_functions/atanh.hpp>
35 #include <boost/math/special_functions/erf.hpp>
36 typedef double (*BoostMathFunction)(double);
37 BoostMathFunction const acosh = boost::math::acosh<double>;
38 BoostMathFunction const asinh = boost::math::asinh<double>;
39 BoostMathFunction const atanh = boost::math::atanh<double>;
40 BoostMathFunction const erf = boost::math::erf<double>;
41 BoostMathFunction const erfc = boost::math::erfc<double>;
42 #endif
43 
44 namespace getfem {
45 
46  base_matrix& __mat_aux1()
47  {
48  DEFINE_STATIC_THREAD_LOCAL(base_matrix, m);
49  return m;
50  }
51 
52 
53 
54  //=========================================================================
55  // Structure dealing with predefined scalar functions.
56  //=========================================================================
57 
58  scalar_type ga_predef_function::operator()(scalar_type t_,
59  scalar_type u_) const {
60  switch(ftype_) {
61  case 0:
62  if (nbargs_ == 2)
63  return (*f2_)(t_, u_);
64  else
65  return (*f1_)(t_);
66  break;
67  case 1:
68  t.thrd_cast()[0] = t_; u.thrd_cast()[0] = u_;
69  workspace.thrd_cast().assembled_potential() = scalar_type(0);
70  ga_function_exec(*gis);
71  return workspace.thrd_cast().assembled_potential();
72  break;
73  }
74  return 0.;
75  }
76 
77  bool ga_predef_function::is_affine(const std::string &varname) const {
78  if (ftype_ == 1) {
79  for (size_type i = 0; i < workspace.thrd_cast().nb_trees(); ++i) {
80  const ga_workspace::tree_description &
81  td = workspace.thrd_cast().tree_info(i);
82  if (!(ga_is_affine(*(td.ptree), workspace, varname, "")))
83  return false;
84  }
85  return true;
86  }
87  return false;
88  }
89 
90 
91  static scalar_type ga_Heaviside(scalar_type t) { return (t >= 0.) ? 1.: 0.; }
92  static scalar_type ga_pos_part(scalar_type t) { return (t >= 0.) ? t : 0.; }
93  static scalar_type ga_reg_pos_part(scalar_type t, scalar_type eps)
94  { return (t >= eps) ? t-eps/2. : ((t <= 0) ? 0. : t*t/(2.*eps)); }
95  static scalar_type ga_der_reg_pos_part(scalar_type t, scalar_type eps)
96  { return (t >= eps) ? 1. : ((t <= 0) ? 0. : t/eps); }
97  static scalar_type ga_der2_reg_pos_part(scalar_type t, scalar_type eps)
98  { return (t >= eps) ? 0. : ((t <= 0) ? 0. : 1./eps); }
99 
100 
101  static scalar_type ga_half_sqr_pos_part(scalar_type t)
102  { return (t >= 0.) ? 0.5*t*t : 0.; }
103  static scalar_type ga_neg_part(scalar_type t) { return (t >= 0.) ? 0. : -t; }
104  static scalar_type ga_half_sqr_neg_part(scalar_type t)
105  { return (t >= 0.) ? 0. : 0.5*t*t; }
106  static scalar_type ga_sinc(scalar_type t) {// cardinal sine function sin(t)/t
107  if (gmm::abs(t) < 1E-4) {
108  scalar_type t2 = t*t;
109  return 1-t2/6.+ t2*t2/120.;
110  } else {
111  return sin(t)/t;
112  }
113  }
114  static scalar_type ga_sqr(scalar_type t) { return t*t; }
115  static scalar_type ga_max(scalar_type t, scalar_type u)
116  { return std::max(t,u); }
117  static scalar_type ga_min(scalar_type t, scalar_type u)
118  { return std::min(t,u); }
119  static scalar_type ga_abs(scalar_type t) { return gmm::abs(t); }
120  static scalar_type ga_sign(scalar_type t) { return (t >= 0.) ? 1.: -1.; }
121 
122  // Derivatives of predefined functions
123  static scalar_type ga_der_sinc(scalar_type t) {
124  if (gmm::abs(t) < 1E-4) {
125  scalar_type t2 = t*t;
126  return -t/3. + t*t2/30. -t*t2*t2/840.;
127  } else {
128  return (t*cos(t) - sin(t))/(t*t);
129  }
130  }
131  static scalar_type ga_der2_sinc(scalar_type t) {
132  if (gmm::abs(t) < 1E-4) {
133  scalar_type t2 = t*t;
134  return -1./3. + t2/10. -t2*t2/168.;
135  } else {
136  return ((2. - t*t)*sin(t) - 2.*t*cos(t))/(t*t*t);
137  }
138  }
139  static scalar_type ga_der_sqrt(scalar_type t) { return 0.5/sqrt(t); }
140  // static scalar_type ga_der_sqr(scalar_type t) { return 2*t; }
141  static scalar_type ga_der_t_pow(scalar_type t, scalar_type u)
142  { return u*pow(t,u-1.); }
143  static scalar_type ga_der_u_pow(scalar_type t, scalar_type u)
144  { return pow(t,u)*log(gmm::abs(t)); }
145  static scalar_type ga_der_log(scalar_type t) { return 1./t; }
146  static scalar_type ga_der_log10(scalar_type t) { return 1./(t*log(10.)); }
147  static scalar_type ga_der_tanh(scalar_type t)
148  { return 1.-gmm::sqr(tanh(t)); }
149  static scalar_type ga_der_asinh(scalar_type t)
150  { return 1./(sqrt(t*t+1.)); }
151  static scalar_type ga_der_acosh(scalar_type t)
152  { return 1./(sqrt(t*t-1.)); }
153  static scalar_type ga_der_atanh(scalar_type t)
154  { return 1./(1.-t*t); }
155  static scalar_type ga_der_cos(scalar_type t)
156  { return -sin(t); }
157  static scalar_type ga_der_tan(scalar_type t)
158  { return 1.+gmm::sqr(tan(t)); }
159  static scalar_type ga_der_asin(scalar_type t)
160  { return 1./(sqrt(1.-t*t)); }
161  static scalar_type ga_der_acos(scalar_type t)
162  { return -1./(sqrt(1.-t*t)); }
163  static scalar_type ga_der_atan(scalar_type t)
164  { return 1./(1.+t*t); }
165  static scalar_type ga_der_t_atan2(scalar_type t, scalar_type u)
166  { return u/(t*t+u*u); }
167  static scalar_type ga_der_u_atan2(scalar_type t, scalar_type u)
168  { return -t/(t*t+u*u); }
169  static scalar_type ga_der_erf(scalar_type t)
170  { return exp(-t*t)*2./sqrt(M_PI); }
171  static scalar_type ga_der_erfc(scalar_type t)
172  { return -exp(-t*t)*2./sqrt(M_PI); }
173  static scalar_type ga_der_neg_part(scalar_type t)
174  { return (t >= 0) ? 0. : -1.; }
175  static scalar_type ga_der_t_max(scalar_type t, scalar_type u)
176  { return (t-u >= 0) ? 1. : 0.; }
177  static scalar_type ga_der_u_max(scalar_type t, scalar_type u)
178  { return (u-t >= 0) ? 1. : 0.; }
179 
180 
181 
182  //=========================================================================
183  // Structure dealing with predefined operators.
184  //=========================================================================
185 
186  static void ga_init_scalar(bgeot::multi_index &mi) { mi.resize(0); }
187  static void ga_init_square_matrix(bgeot::multi_index &mi, size_type N)
188  { mi.resize(2); mi[0] = mi[1] = N; }
189 
190  // Norm Operator
191  struct norm_operator : public ga_nonlinear_operator {
192  bool result_size(const arg_list &args, bgeot::multi_index &sizes) const {
193  if (args.size() != 1 || args[0]->sizes().size() > 2) return false;
194  ga_init_scalar(sizes);
195  return true;
196  }
197 
198  void value(const arg_list &args, base_tensor &result) const
199  { result[0] = gmm::vect_norm2(args[0]->as_vector()); }
200 
201  // Derivative : u/|u|
202  void derivative(const arg_list &args, size_type,
203  base_tensor &result) const {
204  scalar_type no = gmm::vect_norm2(args[0]->as_vector());
205  if (no == scalar_type(0))
206  gmm::clear(result.as_vector());
207  else
208  gmm::copy(gmm::scaled(args[0]->as_vector(), scalar_type(1)/no),
209  result.as_vector());
210  }
211 
212  // Second derivative : (|u|^2 Id - u x u)/|u|^3
213  void second_derivative(const arg_list &args, size_type, size_type,
214  base_tensor &result) const {
215  const base_tensor &t = *args[0];
216  size_type N = t.size();
217  scalar_type no = gmm::vect_norm2(t.as_vector());
218  scalar_type no3 = no*no*no;
219 
220  if (no < 1E-25) no = 1E-25; // In order to avoid infinite values
221 
222  for (size_type i = 0; i < N; ++i)
223  for (size_type j = 0; j < N; ++j) {
224  result[j*N+i] = - t[i]*t[j] / no3;
225  if (i == j) result[j*N+i] += scalar_type(1)/no;
226  }
227  }
228  };
229 
230  // Norm_sqr Operator
231  struct norm_sqr_operator : public ga_nonlinear_operator {
232  bool result_size(const arg_list &args, bgeot::multi_index &sizes) const {
233  if (args.size() != 1 || args[0]->sizes().size() > 2) return false;
234  ga_init_scalar(sizes);
235  return true;
236  }
237 
238  void value(const arg_list &args, base_tensor &result) const
239  { result[0] = gmm::vect_norm2_sqr(args[0]->as_vector()); }
240 
241  // Derivative : 2*u
242  void derivative(const arg_list &args, size_type,
243  base_tensor &result) const {
244  gmm::copy(gmm::scaled(args[0]->as_vector(), scalar_type(2)),
245  result.as_vector());
246  }
247 
248  // Second derivative : Id
249  void second_derivative(const arg_list &args, size_type, size_type,
250  base_tensor &result) const {
251  const base_tensor &t = *args[0];
252  size_type N = t.size();
253  gmm::clear(result.as_vector());
254  for (size_type i = 0; i < N; ++i)
255  result[i*N+i] = scalar_type(2);
256  }
257  };
258 
259  // Det Operator
260  struct det_operator : public ga_nonlinear_operator {
261  bool result_size(const arg_list &args, bgeot::multi_index &sizes) const {
262  if (args.size() != 1 || args[0]->sizes().size() != 2
263  || args[0]->sizes()[0] != args[0]->sizes()[1]) return false;
264  ga_init_scalar(sizes);
265  return true;
266  }
267 
268  void value(const arg_list &args, base_tensor &result) const {
269  size_type N = args[0]->sizes()[0];
270  // base_matrix M(N, N);
271  // gmm::copy(args[0]->as_vector(), M.as_vector());
272  // result[0] = gmm::lu_det(M);
273  result[0] = bgeot::lu_det(&(*(args[0]->begin())), N);
274  }
275 
276  // Derivative : det(M)M^{-T}
277  void derivative(const arg_list &args, size_type,
278  base_tensor &result) const {
279  size_type N = args[0]->sizes()[0];
280  if (N) {
281  __mat_aux1().base_resize(N, N);
282  gmm::copy(args[0]->as_vector(), __mat_aux1().as_vector());
283  scalar_type det = bgeot::lu_inverse(__mat_aux1());
284  if (det == scalar_type(0))
285  gmm::clear(result.as_vector());
286  else {
287  auto it = result.begin();
288  auto ita = __mat_aux1().begin();
289  for (size_type j = 0; j < N; ++j, ++ita) {
290  auto itaa = ita;
291  *it = (*itaa) * det; ++it;
292  for (size_type i = 1; i < N; ++i, ++it)
293  { itaa += N; *it = (*itaa) * det; }
294  }
295  GA_DEBUG_ASSERT(it == result.end(), "Internal error");
296  }
297  }
298  }
299 
300  // Second derivative : det(M)(M^{-T}@M^{-T} - M^{-T}_{jk}M^{-T}_{li})
301  void second_derivative(const arg_list &args, size_type, size_type,
302  base_tensor &result) const {
303  size_type N = args[0]->sizes()[0];
304  __mat_aux1().base_resize(N, N);
305  gmm::copy(args[0]->as_vector(), __mat_aux1().as_vector());
306  scalar_type det = bgeot::lu_inverse(__mat_aux1());
307  if (det == scalar_type(0))
308  gmm::clear(result.as_vector());
309  else {
310  auto it = result.begin();
311  auto ita = __mat_aux1().begin(), ita_l = ita;
312  for (size_type l = 0; l < N; ++l, ++ita_l) {
313  auto ita_lk = ita_l, ita_jk = ita;
314  for (size_type k = 0; k < N; ++k, ita_lk += N, ita_jk += N) {
315  auto ita_j = ita;
316  for (size_type j = 0; j < N; ++j, ++ita_j, ++ita_jk) {
317  auto ita_ji = ita_j, ita_li = ita_l;
318  for (size_type i = 0; i < N; ++i, ++it, ita_ji += N, ita_li += N)
319  *it = ((*ita_ji) * (*ita_lk) - (*ita_jk) * (*ita_li)) * det;
320  }
321  }
322  }
323  GA_DEBUG_ASSERT(it == result.end(), "Internal error");
324  }
325  }
326  };
327 
328  // Inverse Operator (for square matrices)
329  struct inverse_operator : public ga_nonlinear_operator {
330  bool result_size(const arg_list &args, bgeot::multi_index &sizes) const {
331  if (args.size() != 1 || args[0]->sizes().size() != 2
332  || args[0]->sizes()[0] != args[0]->sizes()[1]) return false;
333  ga_init_square_matrix(sizes, args[0]->sizes()[0]);
334  return true;
335  }
336 
337  // Value : M^{-1}
338  void value(const arg_list &args, base_tensor &result) const {
339  size_type N = args[0]->sizes()[0];
340  __mat_aux1().base_resize(N, N);
341  gmm::copy(args[0]->as_vector(), __mat_aux1().as_vector());
342  bgeot::lu_inverse(__mat_aux1());
343  gmm::copy(__mat_aux1().as_vector(), result.as_vector());
344  }
345 
346  // Derivative : -M^{-1}{ik}M^{-1}{lj} (comes from H -> M^{-1}HM^{-1})
347  void derivative(const arg_list &args, size_type,
348  base_tensor &result) const { // to be verified
349  size_type N = args[0]->sizes()[0];
350  __mat_aux1().base_resize(N, N);
351  gmm::copy(args[0]->as_vector(), __mat_aux1().as_vector());
352  bgeot::lu_inverse(__mat_aux1());
353  auto it = result.begin();
354  auto ita = __mat_aux1().begin(), ita_l = ita;
355  for (size_type l = 0; l < N; ++l, ++ita_l) {
356  auto ita_k = ita;
357  for (size_type k = 0; k < N; ++k, ita_k += N) {
358  auto ita_lj = ita_l;
359  for (size_type j = 0; j < N; ++j, ++ita_lj) {
360  auto ita_ik = ita_k;
361  for (size_type i = 0; i < N; ++i, ++it, ++ita_ik)
362  *it = -(*ita_ik) * (*ita_lj);
363  }
364  }
365  }
366  GA_DEBUG_ASSERT(it == result.end(), "Internal error");
367  }
368 
369  // Second derivative :
370  // M^{-1}{ik}M^{-1}{lm}M^{-1}{nj} + M^{-1}{im}M^{-1}{mk}M^{-1}{lj}
371  // comes from (H,K) -> M^{-1}HM^{-1}KM^{-1} + M^{-1}KM^{-1}HM^{-1}
372  void second_derivative(const arg_list &args, size_type, size_type,
373  base_tensor &result) const { // to be verified
374  size_type N = args[0]->sizes()[0];
375  __mat_aux1().base_resize(N, N);
376  gmm::copy(args[0]->as_vector(), __mat_aux1().as_vector());
377  bgeot::lu_inverse(__mat_aux1());
378  base_tensor::iterator it = result.begin();
379  for (size_type n = 0; n < N; ++n)
380  for (size_type m = 0; m < N; ++m)
381  for (size_type l = 0; l < N; ++l)
382  for (size_type k = 0; k < N; ++k)
383  for (size_type j = 0; j < N; ++j)
384  for (size_type i = 0; i < N; ++i, ++it)
385  *it = __mat_aux1()(i,k)*__mat_aux1()(l,m)*__mat_aux1()(n,j)
386  + __mat_aux1()(i,m)*__mat_aux1()(m,k)*__mat_aux1()(l,j);
387  GA_DEBUG_ASSERT(it == result.end(), "Internal error");
388  }
389  };
390 
391  //=========================================================================
392  // Initialization of predefined functions and operators.
393  //=========================================================================
394 
395  ga_predef_function::ga_predef_function()
396  : expr_(""), derivative1_(""), derivative2_(""), gis(nullptr) {}
397 
398  ga_predef_function::ga_predef_function(pscalar_func_onearg f,
399  size_type dtype__,
400  const std::string &der)
401  : ftype_(0), dtype_(dtype__), nbargs_(1), f1_(f), expr_(""),
402  derivative1_(der), derivative2_("") {}
403  ga_predef_function::ga_predef_function(pscalar_func_twoargs f,
404  size_type dtype__,
405  const std::string &der1,
406  const std::string &der2)
407  : ftype_(0), dtype_(dtype__), nbargs_(2), f2_(f),
408  expr_(""), derivative1_(der1), derivative2_(der2), gis(nullptr) {}
409  ga_predef_function::ga_predef_function(const std::string &expr__)
410  : ftype_(1), dtype_(3), nbargs_(1), expr_(expr__),
411  derivative1_(""), derivative2_(""), t(1, 0.), u(1, 0.), gis(nullptr) {}
412 
413 
414  ga_predef_function_tab::ga_predef_function_tab() {
415 
416  ga_predef_function_tab &PREDEF_FUNCTIONS = *this;
417 
418  // Power functions and their derivatives
419  PREDEF_FUNCTIONS["sqrt"] = ga_predef_function(sqrt, 1, "DER_PDFUNC_SQRT");
420  PREDEF_FUNCTIONS["sqr"] = ga_predef_function(ga_sqr, 2, "2*t");
421  PREDEF_FUNCTIONS["pow"] = ga_predef_function(pow, 1, "DER_PDFUNC1_POW",
422  "DER_PDFUNC2_POW");
423 
424  PREDEF_FUNCTIONS["DER_PDFUNC_SQRT"] =
425  ga_predef_function(ga_der_sqrt, 2, "-0.25/(t*sqrt(t))");
426  PREDEF_FUNCTIONS["DER_PDFUNC1_POW"] =
427  ga_predef_function(ga_der_t_pow, 2, "u*(u-1)*pow(t,u-2)",
428  "pow(t,u-1)*(u*log(t)+1)");
429  PREDEF_FUNCTIONS["DER_PDFUNC2_POW"] =
430  ga_predef_function(ga_der_u_pow, 2, "pow(t,u-1)*(u*log(t)+1)",
431  "pow(t,u)*sqr(log(t))");
432 
433  // Hyperbolic functions
434  PREDEF_FUNCTIONS["exp"] = ga_predef_function(exp, 1, "exp");
435  PREDEF_FUNCTIONS["log"] = ga_predef_function(log, 1, "DER_PDFUNC_LOG");
436  PREDEF_FUNCTIONS["log10"] =
437  ga_predef_function(log10, 1, "DER_PDFUNC_LOG10");
438  PREDEF_FUNCTIONS["sinh"] = ga_predef_function(sinh, 1, "cosh");
439  PREDEF_FUNCTIONS["cosh"] = ga_predef_function(cosh, 1, "sinh");
440  PREDEF_FUNCTIONS["tanh"] = ga_predef_function(tanh, 1, "DER_PDFUNC_TANH");
441  PREDEF_FUNCTIONS["asinh"] =
442  ga_predef_function(asinh, 1, "DER_PDFUNC_ASINH");
443  PREDEF_FUNCTIONS["acosh"] =
444  ga_predef_function(acosh, 1, "DER_PDFUNC_ACOSH");
445  PREDEF_FUNCTIONS["atanh"] =
446  ga_predef_function(atanh, 1, "DER_PDFUNC_ATANH");
447 
448 
449  PREDEF_FUNCTIONS["DER_PDFUNC_LOG"] =
450  ga_predef_function(ga_der_log, 2, "-1/sqr(t)");
451  PREDEF_FUNCTIONS["DER_PDFUNC_LOG10"] =
452  ga_predef_function(ga_der_log10, 2, "-1/(sqr(t)*log(10))");
453  PREDEF_FUNCTIONS["DER_PDFUNC_TANH"] =
454  ga_predef_function(ga_der_tanh, 2, "2*tanh(t)*(sqr(tanh(t))-1)");
455  PREDEF_FUNCTIONS["DER_PDFUNC_ASINH"] =
456  ga_predef_function(ga_der_asinh, 2, "-t/(pow(t*t+1,1.5))");
457  PREDEF_FUNCTIONS["DER_PDFUNC_ACOSH"] =
458  ga_predef_function(ga_der_acosh, 2, "-t/(pow(t*t-1,1.5))");
459  PREDEF_FUNCTIONS["DER_PDFUNC_ATANH"] =
460  ga_predef_function(ga_der_atanh, 2, "2*t/sqr(1-t*t)");
461 
462 
463  // Trigonometric functions
464  PREDEF_FUNCTIONS["sin"] = ga_predef_function(sin, 1, "cos");
465  PREDEF_FUNCTIONS["cos"] = ga_predef_function(cos, 1, "DER_PDFUNC_COS");
466  PREDEF_FUNCTIONS["tan"] = ga_predef_function(tan, 1, "DER_PDFUNC_TAN");
467  PREDEF_FUNCTIONS["asin"] = ga_predef_function(asin, 1, "DER_PDFUNC_ASIN");
468  PREDEF_FUNCTIONS["acos"] = ga_predef_function(acos, 1, "DER_PDFUNC_ACOS");
469  PREDEF_FUNCTIONS["atan"] = ga_predef_function(atan, 1, "DER_PDFUNC_ATAN");
470  PREDEF_FUNCTIONS["atan2"] = ga_predef_function(atan2,1,"DER_PDFUNC1_ATAN2",
471  "DER_PDFUNC2_ATAN2");
472  PREDEF_FUNCTIONS["sinc"] = ga_predef_function(ga_sinc, 1,
473  "DER_PDFUNC_SINC");
474  PREDEF_FUNCTIONS["DER_PDFUNC_SINC"] =ga_predef_function(ga_der_sinc, 1,
475  "DER2_PDFUNC_SINC");
476  PREDEF_FUNCTIONS["DER2_PDFUNC_SINC"] = ga_predef_function(ga_der2_sinc);
477 
478 
479  PREDEF_FUNCTIONS["DER_PDFUNC_COS"] =
480  ga_predef_function(ga_der_cos, 2, "-cos(t)");
481  PREDEF_FUNCTIONS["DER_PDFUNC_TAN"] =
482  ga_predef_function(ga_der_tan, 2, "2*tan(t)/sqr(cos(t))");
483  // PREDEF_FUNCTIONS["DER_PDFUNC_TAN"] =
484  // ga_predef_function(ga_der_tan, 2, "2*tan(t)*(1+sqr(tan(t)))");
485  PREDEF_FUNCTIONS["DER_PDFUNC_ASIN"] =
486  ga_predef_function(ga_der_asin, 2, "t/(pow(1-t*t,1.5))");
487  PREDEF_FUNCTIONS["DER_PDFUNC_ACOS"] =
488  ga_predef_function(ga_der_acos, 2, "-t/(pow(1-t*t,1.5))");
489  PREDEF_FUNCTIONS["DER_PDFUNC_ATAN"] =
490  ga_predef_function(ga_der_atan, 2, "-2*t/sqr(1+t*t)");
491  PREDEF_FUNCTIONS["DER_PDFUNC1_ATAN2"] =
492  ga_predef_function(ga_der_t_atan2, 2, "-2*t*u/sqr(sqr(u)+sqr(t))",
493  "(sqrt(t)-sqr(u))/sqr(sqr(u)+sqr(t))");
494  PREDEF_FUNCTIONS["DER_PDFUNC2_ATAN2"] =
495  ga_predef_function(ga_der_u_atan2, 2,
496  "(sqrt(t)-sqr(u))/sqr(sqr(u)+sqr(t))",
497  "2*t*u/sqr(sqr(u)+sqr(t))");
498 
499 
500  // Error functions
501  PREDEF_FUNCTIONS["erf"]
502  = ga_predef_function(erf, 1, "DER_PDFUNC_ERF");
503  PREDEF_FUNCTIONS["erfc"]
504  = ga_predef_function(erfc, 1, "DER_PDFUNC_ERFC");
505 
506  PREDEF_FUNCTIONS["DER_PDFUNC_ERF"] =
507  ga_predef_function(ga_der_erf, 2, "exp(-t*t)*2/sqrt(pi)");
508  PREDEF_FUNCTIONS["DER_PDFUNC_ERFC"] =
509  ga_predef_function(ga_der_erfc, 2, "-exp(-t*t)*2/sqrt(pi)");
510 
511 
512 
513  // Miscellaneous functions
514  PREDEF_FUNCTIONS["Heaviside"] = ga_predef_function(ga_Heaviside); // ga_predef_function(ga_Heaviside, 2, "(0)");
515  PREDEF_FUNCTIONS["sign"] = ga_predef_function(ga_sign);
516  PREDEF_FUNCTIONS["abs"] = ga_predef_function(ga_abs, 1, "sign");
517  PREDEF_FUNCTIONS["pos_part"]
518  = ga_predef_function(ga_pos_part, 1, "Heaviside");
519  PREDEF_FUNCTIONS["half_sqr_pos_part"]
520  = ga_predef_function(ga_half_sqr_pos_part, 1, "pos_part");
521  PREDEF_FUNCTIONS["neg_part"]
522  = ga_predef_function(ga_neg_part, 1, "DER_PDFUNC_NEG_PART");
523  PREDEF_FUNCTIONS["half_sqr_neg_part"]
524  = ga_predef_function(ga_half_sqr_neg_part, 2, "-neg_part(t)");
525  PREDEF_FUNCTIONS["reg_pos_part"]
526  = ga_predef_function(ga_reg_pos_part, 1, "DER_REG_POS_PART", "");
527  PREDEF_FUNCTIONS["DER_REG_POS_PART"]
528  = ga_predef_function(ga_der_reg_pos_part, 1, "DER2_REG_POS_PART", "");
529  PREDEF_FUNCTIONS["DER_REG_POS_PART"]
530  = ga_predef_function(ga_der2_reg_pos_part);
531 
532  PREDEF_FUNCTIONS["max"]
533  = ga_predef_function(ga_max, 1, "DER_PDFUNC1_MAX", "DER_PDFUNC2_MAX");
534  PREDEF_FUNCTIONS["min"]
535  = ga_predef_function(ga_min, 1, "DER_PDFUNC2_MAX", "DER_PDFUNC1_MAX");
536 
537  PREDEF_FUNCTIONS["DER_PDFUNC_NEG_PART"]
538  = ga_predef_function(ga_der_neg_part, 2, "-Heaviside(-t)");
539  PREDEF_FUNCTIONS["DER_PDFUNC1_MAX"] = ga_predef_function(ga_der_t_max);
540  PREDEF_FUNCTIONS["DER_PDFUNC2_MAX"] = ga_predef_function(ga_der_u_max);
541 
542  }
543 
544  ga_spec_function_tab::ga_spec_function_tab() {
545  // Predefined special functions
546  ga_spec_function_tab &SPEC_FUNCTIONS = *this;
547 
548  SPEC_FUNCTIONS.insert("pi");
549  SPEC_FUNCTIONS.insert("meshdim");
550  SPEC_FUNCTIONS.insert("timestep");
551  SPEC_FUNCTIONS.insert("qdim");
552  SPEC_FUNCTIONS.insert("qdims");
553  SPEC_FUNCTIONS.insert("Id");
554  }
555 
556  ga_spec_op_tab::ga_spec_op_tab() {
557  // Predefined special operators
558  ga_spec_op_tab &SPEC_OP = *this;
559  SPEC_OP.insert("X");
560  SPEC_OP.insert("element_size");
561  SPEC_OP.insert("element_K");
562  SPEC_OP.insert("element_B");
563  SPEC_OP.insert("Normal");
564  SPEC_OP.insert("Sym");
565  SPEC_OP.insert("Skew");
566  SPEC_OP.insert("Def");
567  SPEC_OP.insert("Trace");
568  SPEC_OP.insert("Deviator");
569  SPEC_OP.insert("Interpolate");
570  SPEC_OP.insert("Interpolate_filter");
571  SPEC_OP.insert("Elementary_transformation");
572  SPEC_OP.insert("Xfem_plus");
573  SPEC_OP.insert("Xfem_minus");
574  SPEC_OP.insert("Print");
575  SPEC_OP.insert("Reshape");
576  SPEC_OP.insert("Swap_indices");
577  SPEC_OP.insert("Index_move_last");
578  SPEC_OP.insert("Contract");
579  SPEC_OP.insert("Diff");
580  SPEC_OP.insert("Grad");
581  }
582 
583 
584  ga_predef_operator_tab::ga_predef_operator_tab(void) {
585  // Predefined operators
586  ga_predef_operator_tab &PREDEF_OPERATORS = *this;
587 
588  PREDEF_OPERATORS.add_method("Norm", std::make_shared<norm_operator>());
589  PREDEF_OPERATORS.add_method("Norm_sqr",
590  std::make_shared<norm_sqr_operator>());
591  PREDEF_OPERATORS.add_method("Det", std::make_shared<det_operator>());
592  PREDEF_OPERATORS.add_method("Inv", std::make_shared<inverse_operator>());
593  }
594 
595 
596 
597  bool ga_function_exists(const std::string &name) {
598  const ga_predef_function_tab &PREDEF_FUNCTIONS
599  = dal::singleton<ga_predef_function_tab>::instance(0);
600  return PREDEF_FUNCTIONS.find(name) != PREDEF_FUNCTIONS.end();
601  }
602 
603 
604  void ga_define_function(const std::string &name, size_type nbargs,
605  const std::string &expr, const std::string &der1,
606  const std::string &der2) {
607  auto guard = omp_guard{};
608 
609  auto &PREDEF_FUNCTIONS = dal::singleton<ga_predef_function_tab>::instance(0);
610  if(PREDEF_FUNCTIONS.find(name) != PREDEF_FUNCTIONS.end()) return;
611  GMM_ASSERT1(nbargs >= 1 && nbargs <= 2, "Generic assembly only allows "
612  "the definition of scalar function with one or two arguments");
613  { // Only for syntax analysis
614  base_vector t(1);
615  ga_workspace workspace;
616  workspace.add_fixed_size_variable("t", gmm::sub_interval(0,1), t);
617  if (nbargs == 2)
618  workspace.add_fixed_size_variable("u", gmm::sub_interval(0,1), t);
619  workspace.add_function_expression(expr);
620  }
621 
622  PREDEF_FUNCTIONS[name] = ga_predef_function(expr);
623  ga_predef_function &F = PREDEF_FUNCTIONS[name];
624  F.gis = std::make_unique<instruction_set>();
625  for (size_type thread = 0; thread < num_threads(); ++thread)
626  {
627  F.workspace(thread).add_fixed_size_variable("t", gmm::sub_interval(0,1),
628  F.t(thread));
629  if (nbargs == 2)
630  F.workspace(thread).add_fixed_size_variable("u", gmm::sub_interval(0,1),
631  F.u(thread));
632  F.workspace(thread).add_function_expression(expr);
633  ga_compile_function(F.workspace(thread), (*F.gis)(thread), true);
634  }
635  F.nbargs_ = nbargs;
636  if (nbargs == 1) {
637  if (der1.size()) { F.derivative1_ = der1; F.dtype_ = 2; }
638  } else {
639  if (der1.size() && der2.size()) {
640  F.derivative1_ = der1; F.derivative2_ = der2; F.dtype_ = 2;
641  }
642  }
643  }
644 
645  void ga_define_function(const std::string &name, pscalar_func_onearg f,
646  const std::string &der) {
647  ga_predef_function_tab &PREDEF_FUNCTIONS
648  = dal::singleton<ga_predef_function_tab>::instance(0);
649  PREDEF_FUNCTIONS[name] = ga_predef_function(f, 1, der);
650  ga_predef_function &F = PREDEF_FUNCTIONS[name];
651  if (der.size() == 0) F.dtype_ = 0;
652  else if (!(ga_function_exists(der))) F.dtype_ = 2;
653  }
654 
655  void ga_define_function(const std::string &name, pscalar_func_twoargs f,
656  const std::string &der1, const std::string &der2) {
657  ga_predef_function_tab &PREDEF_FUNCTIONS
658  = dal::singleton<ga_predef_function_tab>::instance(0);
659  PREDEF_FUNCTIONS[name] = ga_predef_function(f, 1, der1, der2);
660  ga_predef_function &F = PREDEF_FUNCTIONS[name];
661  if (der1.size() == 0 || der2.size() == 0)
662  F.dtype_ = 0;
663  else if (!(ga_function_exists(der1)) || !(ga_function_exists(der2)))
664  F.dtype_ = 2;
665  }
666 
667  void ga_undefine_function(const std::string &name) {
668  ga_predef_function_tab &PREDEF_FUNCTIONS
669  = dal::singleton<ga_predef_function_tab>::instance(0);
670  ga_predef_function_tab::iterator it = PREDEF_FUNCTIONS.find(name);
671  if (it != PREDEF_FUNCTIONS.end()) {
672  PREDEF_FUNCTIONS.erase(name);
673  std::string name0 = "DER_PDFUNC_" + name;
674  ga_undefine_function(name0);
675  std::string name1 = "DER_PDFUNC1_" + name;
676  ga_undefine_function(name1);
677  std::string name2 = "DER_PDFUNC2_" + name;
678  ga_undefine_function(name2);
679  }
680  }
681 
682  //=========================================================================
683  // User defined functions
684  //=========================================================================
685 
686  ga_function::ga_function(const ga_workspace &workspace_,
687  const std::string &e)
688  : local_workspace(true, workspace_), expr(e), gis(0) {}
689 
690  ga_function::ga_function(const model &md, const std::string &e)
691  : local_workspace(md), expr(e), gis(0) {}
692 
693  ga_function::ga_function(const std::string &e)
694  : local_workspace(), expr(e), gis(0) {}
695 
696  ga_function::ga_function(const ga_function &gaf)
697  : local_workspace(gaf.local_workspace), expr(gaf.expr), gis(0)
698  { if (gaf.gis) compile(); }
699 
700  void ga_function::compile() const {
701  if (gis) delete gis;
702  gis = new ga_instruction_set;
703  local_workspace.clear_expressions();
704  local_workspace.add_function_expression(expr);
705  ga_compile_function(local_workspace, *gis, false);
706  }
707 
708  ga_function &ga_function::operator =(const ga_function &gaf) {
709  if (gis) delete gis;
710  gis = 0;
711  local_workspace = gaf.local_workspace;
712  expr = gaf.expr;
713  if (gaf.gis) compile();
714  return *this;
715  }
716 
717  ga_function::~ga_function() { if (gis) delete gis; gis = 0; }
718 
719  const base_tensor &ga_function::eval() const {
720  GMM_ASSERT1(gis, "Uncompiled function");
721  gmm::clear(local_workspace.assembled_tensor().as_vector());
722  ga_function_exec(*gis);
723  return local_workspace.assembled_tensor();
724  }
725 
726  void ga_function::derivative(const std::string &var) {
727  GMM_ASSERT1(gis, "Uncompiled function");
728  if (local_workspace.nb_trees()) {
729  ga_tree tree = *(local_workspace.tree_info(0).ptree);
730  ga_derivative(tree, local_workspace, *((const mesh *)(0)), var, "", 1);
731  if (tree.root) {
732  ga_semantic_analysis(tree, local_workspace, *((const mesh *)(0)),
733  1, false, true);
734  // To be improved to suppress test functions in the expression ...
735  // ga_replace_test_by_cte do not work in all operations like
736  // vector components x(1)
737  // ga_replace_test_by_cte(tree.root, false);
738  // ga_semantic_analysis(tree, local_workspace, *((const mesh *)(0)), 1,
739  // false, true);
740  }
741  expr = ga_tree_to_string(tree);
742  }
743  if (gis) delete gis;
744  gis = 0;
745  compile();
746  }
747 
748 } /* end of namespace */
gmm::vect_norm2_sqr
number_traits< typename linalg_traits< V >::value_type >::magnitude_type vect_norm2_sqr(const V &v)
squared Euclidean norm of a vector.
Definition: gmm_blas.h:544
getfem_generic_assembly_semantic.h
Semantic analysis of assembly trees and semantic manipulations.
dal::singleton::instance
static T & instance()
Instance from the current thread.
Definition: dal_singleton.h:146
gmm::vect_norm2
number_traits< typename linalg_traits< V >::value_type >::magnitude_type vect_norm2(const V &v)
Euclidean norm of a vector.
Definition: gmm_blas.h:557
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_assembling.h:48
gmm::clear
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:59