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