Welcome to the NetCologne GmbH open source mirroring service!

This machine mirrors various open-source projects. 20 Gbit/s uplink.

If there are any issues or you want another project mirrored, please contact mirror-service -=AT=- netcologne DOT de !

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