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