GetFEM  5.4.2
getfem_global_function.cc
1 /*===========================================================================
2 
3  Copyright (C) 2004-2020 Yves Renard
4  Copyright (C) 2016-2020 Konstantinos Poulios
5 
6  This file is a part of GetFEM
7 
8  GetFEM is free software; you can redistribute it and/or modify it
9  under the terms of the GNU Lesser General Public License as published
10  by the Free Software Foundation; either version 3 of the License, or
11  (at your option) any later version along with the GCC Runtime Library
12  Exception either version 3.1 or (at your option) any later version.
13  This program is distributed in the hope that it will be useful, but
14  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16  License and GCC Runtime Library Exception for more details.
17  You should have received a copy of the GNU Lesser General Public License
18  along with this program; if not, write to the Free Software Foundation,
19  Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
20 
21 ===========================================================================*/
22 
23 #include <getfem/getfem_global_function.h>
24 #include <getfem/getfem_level_set.h>
25 
26 namespace getfem {
27 
28 
29  // Partial implementation of abstract class global_function_simple
30 
31  scalar_type global_function_simple::val
32  (const fem_interpolation_context &c) const {
33  base_node pt = c.xreal();
34  GMM_ASSERT1(pt.size() == dim_, "Point of wrong size (" << pt.size() << ") "
35  << "passed to a global function of dim = "<< dim_ <<".");
36  return this->val(pt);
37  }
38 
39  void global_function_simple::grad
40  (const fem_interpolation_context &c, base_small_vector &g) const {
41  base_node pt = c.xreal();
42  GMM_ASSERT1(pt.size() == dim_, "Point of wrong size (" << pt.size() << ") "
43  << "passed to a global function of dim = "<< dim_ <<".");
44  this->grad(pt, g);
45  }
46 
47  void global_function_simple::hess
48  (const fem_interpolation_context &c, base_matrix &h) const {
49  base_node pt = c.xreal();
50  GMM_ASSERT1(pt.size() == dim_, "Point of wrong size (" << pt.size() << ") "
51  << "passed to a global function of dim = "<< dim_ <<".");
52  this->hess(pt, h);
53  }
54 
55  // Implementation of global_function_parser
56 
57  scalar_type global_function_parser::val(const base_node &pt) const {
58  const bgeot::base_tensor &t = tensor_val(pt);
59  GMM_ASSERT1(t.size() == 1, "Wrong size of expression result "
60  << f_val.expression());
61  return t[0];
62  }
63 
64  const base_tensor &global_function_parser::tensor_val(const base_node &pt) const {
65  gmm::copy(pt, pt_);
66  return f_val.eval();
67  }
68 
69  void global_function_parser::grad(const base_node &pt, base_small_vector &g) const {
70  g.resize(dim_);
71  gmm::copy(pt, pt_);
72  const bgeot::base_tensor &t = f_grad.eval();
73  GMM_ASSERT1(t.size() == dim_, "Wrong size of expression result "
74  << f_grad.expression());
75  gmm::copy(t.as_vector(), g);
76 
77  }
78 
79  void global_function_parser::hess(const base_node &pt, base_matrix &h) const {
80  h.resize(dim_, dim_);
81  gmm::copy(pt, pt_);
82  const bgeot::base_tensor &t = f_hess.eval();
83  GMM_ASSERT1(t.size() == size_type(dim_*dim_),
84  "Wrong size of expression result " << f_hess.expression());
85  gmm::copy(t.as_vector(), h.as_vector());
86  }
87 
88  global_function_parser::global_function_parser(dim_type dim__,
89  const std::string &sval,
90  const std::string &sgrad,
91  const std::string &shess)
92  : global_function_simple(dim__),
93  f_val(gw, sval), f_grad(gw, sgrad), f_hess(gw, shess) {
94 
95  size_type N(dim_);
96  pt_.resize(N);
97  gmm::fill(pt_, scalar_type(0));
98  gw.add_fixed_size_variable("X", gmm::sub_interval(0, N), pt_);
99  if (N >= 1) gw.add_macro("x", "X(1)");
100  if (N >= 2) gw.add_macro("y", "X(2)");
101  if (N >= 3) gw.add_macro("z", "X(3)");
102  if (N >= 4) gw.add_macro("w", "X(4)");
103 
104  f_val.compile();
105  f_grad.compile();
106  f_hess.compile();
107  }
108 
109  // Implementation of global_function_sum
110 
111  scalar_type global_function_sum::val
112  (const fem_interpolation_context &c) const {
113  scalar_type res(0);
114  for (const auto &f : functions)
115  res += f->val(c);
116  return res;
117  }
118 
119  void global_function_sum::grad
120  (const fem_interpolation_context &c, base_small_vector &g) const {
121  g.resize(dim_);
122  gmm::clear(g);
123  base_small_vector gg(dim_);
124  for (const auto &f : functions) {
125  f->grad(c, gg);
126  gmm::add(gg, g);
127  }
128  }
129 
130  void global_function_sum::hess
131  (const fem_interpolation_context &c, base_matrix &h) const {
132  h.resize(dim_, dim_);
133  gmm::clear(h);
134  base_matrix hh(dim_, dim_);
135  for (const auto &f : functions) {
136  f->hess(c, hh);
137  gmm::add(hh.as_vector(), h.as_vector());
138  }
139  }
140 
141  bool global_function_sum::is_in_support(const base_node &p) const {
142  for (const auto &f : functions)
143  if (f->is_in_support(p)) return true;
144  return false;
145  }
146 
147  void global_function_sum::bounding_box
148  (base_node &bmin_, base_node &bmax_) const {
149  if (functions.size() > 0)
150  functions[0]->bounding_box(bmin_, bmax_);
151  base_node bmin0(dim()), bmax0(dim());
152  for (const auto &f : functions) {
153  f->bounding_box(bmin0, bmax0);
154  for (size_type i=0; i < dim(); ++i) {
155  if (bmin0[i] < bmin_[i]) bmin_[i] = bmin0[i];
156  if (bmax0[i] > bmax_[i]) bmax_[i] = bmax0[i];
157  }
158  }
159  }
160 
161  global_function_sum::global_function_sum(const std::vector<pglobal_function> &funcs)
162  : global_function((funcs.size() > 0) ? funcs[0]->dim() : 0), functions(funcs) {
163  for (const auto &f : functions)
164  GMM_ASSERT1(f->dim() == dim(), "Incompatible dimensions among the provided"
165  " global functions");
166  }
167 
168  global_function_sum::global_function_sum(pglobal_function f1, pglobal_function f2)
169  : global_function(f1->dim()), functions(2) {
170  functions[0] = f1;
171  functions[1] = f2;
172  GMM_ASSERT1(f1->dim() == dim() && f2->dim() == dim(),
173  "Incompatible dimensions between the provided global functions");
174  }
175 
176  global_function_sum::global_function_sum(pglobal_function f1, pglobal_function f2,
177  pglobal_function f3)
178  : global_function(f1->dim()), functions(3) {
179  functions[0] = f1;
180  functions[1] = f2;
181  functions[2] = f3;
182  GMM_ASSERT1(f1->dim() == dim() && f2->dim() == dim() && f3->dim() == dim(),
183  "Incompatible dimensions between the provided global functions");
184  }
185 
186  global_function_sum::global_function_sum(pglobal_function f1, pglobal_function f2,
187  pglobal_function f3, pglobal_function f4)
188  : global_function(f1->dim()), functions(4) {
189  functions[0] = f1;
190  functions[1] = f2;
191  functions[2] = f3;
192  functions[3] = f4;
193  GMM_ASSERT1(f1->dim() == dim() && f2->dim() == dim() && f3->dim() == dim(),
194  "Incompatible dimensions between the provided global functions");
195  }
196 
197 
198  // Implementation of global_function_product
199 
200  scalar_type global_function_product::val
201  (const fem_interpolation_context &c) const {
202  return f1->val(c) * f2->val(c);
203  }
204 
205  void global_function_product::grad
206  (const fem_interpolation_context &c, base_small_vector &g) const {
207  g.resize(dim_);
208  base_small_vector gg(dim_);
209  f1->grad(c, gg);
210  gmm::copy(gmm::scaled(gg, f2->val(c)), g);
211  f2->grad(c, gg);
212  gmm::add(gmm::scaled(gg, f1->val(c)), g);
213  }
214 
215  void global_function_product::hess
216  (const fem_interpolation_context &c, base_matrix &h) const {
217  h.resize(dim_, dim_);
218  gmm::clear(h);
219  base_matrix hh(dim_, dim_);
220  f1->hess(c, hh);
221  gmm::copy(gmm::scaled(hh, f2->val(c)), h);
222  f2->hess(c, hh);
223  gmm::add(gmm::scaled(hh, f1->val(c)), h);
224  base_small_vector g1(dim_), g2(dim_);
225  f1->grad(c, g1);
226  f2->grad(c, g2);
227  gmm::rank_one_update(h, g1, g2);
228  gmm::rank_one_update(h, g2, g1);
229  }
230 
231  bool global_function_product::is_in_support(const base_node &p) const {
232  return f1->is_in_support(p) && f2->is_in_support(p);
233  }
234 
235  void global_function_product::bounding_box
236  (base_node &bmin_, base_node &bmax_) const {
237  base_node bmin0(dim()), bmax0(dim());
238  f1->bounding_box(bmin0, bmax0);
239  f2->bounding_box(bmin_, bmax_);
240  for (size_type i=0; i < dim(); ++i) {
241  if (bmin0[i] > bmin_[i]) bmin_[i] = bmin0[i];
242  if (bmax0[i] < bmax_[i]) bmax_[i] = bmax0[i];
243  if (bmin_[i] > bmax_[i])
244  GMM_WARNING1("Global function product with vanishing basis function");
245  }
246  }
247 
248  global_function_product::global_function_product(pglobal_function f1_, pglobal_function f2_)
249  : global_function(f1_->dim()), f1(f1_), f2(f2_) {
250  GMM_ASSERT1(f2->dim() == dim(), "Incompatible dimensions between the provided"
251  " global functions");
252  }
253 
254 
255  // Implementation of global_function_bounded
256 
257  bool global_function_bounded::is_in_support(const base_node &pt) const {
258  if (has_expr) {
259  gmm::copy(pt, pt_);
260  const bgeot::base_tensor &t = f_val.eval();
261  GMM_ASSERT1(t.size() == 1, "Wrong size of expression result "
262  << f_val.expression());
263  return (t[0] > scalar_type(0));
264  }
265  return true;
266  }
267 
268  global_function_bounded::global_function_bounded(pglobal_function f_,
269  const base_node &bmin_,
270  const base_node &bmax_,
271  const std::string &is_in_expr)
272  : global_function(f_->dim()), f(f_), bmin(bmin_), bmax(bmax_),
273  f_val(gw, is_in_expr) {
274 
275  has_expr = !is_in_expr.empty();
276  if (has_expr) {
277  size_type N(dim_);
278  pt_.resize(N);
279  gmm::fill(pt_, scalar_type(0));
280  gw.add_fixed_size_variable("X", gmm::sub_interval(0, N), pt_);
281  if (N >= 1) gw.add_macro("x", "X(1)");
282  if (N >= 2) gw.add_macro("y", "X(2)");
283  if (N >= 3) gw.add_macro("z", "X(3)");
284  if (N >= 4) gw.add_macro("w", "X(4)");
285  f_val.compile();
286  }
287  }
288 
289  // Implementation of some useful xy functions
290 
291  parser_xy_function::parser_xy_function(const std::string &sval,
292  const std::string &sgrad,
293  const std::string &shess)
294  : f_val(gw, sval), f_grad(gw, sgrad), f_hess(gw, shess),
295  ptx(1), pty(1), ptr(1), ptt(1) {
296 
297  gw.add_fixed_size_constant("x", ptx);
298  gw.add_fixed_size_constant("y", pty);
299  gw.add_fixed_size_constant("r", ptr);
300  gw.add_fixed_size_constant("theta", ptt);
301 
302  f_val.compile();
303  f_grad.compile();
304  f_hess.compile();
305  }
306 
307  scalar_type
308  parser_xy_function::val(scalar_type x, scalar_type y) const {
309  ptx[0] = double(x); // x
310  pty[0] = double(y); // y
311  ptr[0] = double(sqrt(fabs(x*x+y*y))); // r
312  ptt[0] = double(atan2(y,x)); // theta
313 
314  const bgeot::base_tensor &t = f_val.eval();
315  GMM_ASSERT1(t.size() == 1, "Wrong size of expression result "
316  << f_val.expression());
317  return t[0];
318  }
319 
320  base_small_vector
321  parser_xy_function::grad(scalar_type x, scalar_type y) const {
322  ptx[0] = double(x); // x
323  pty[0] = double(y); // y
324  ptr[0] = double(sqrt(fabs(x*x+y*y))); // r
325  ptt[0] = double(atan2(y,x)); // theta
326 
327  base_small_vector res(2);
328  const bgeot::base_tensor &t = f_grad.eval();
329  GMM_ASSERT1(t.size() == 2, "Wrong size of expression result "
330  << f_grad.expression());
331  gmm::copy(t.as_vector(), res);
332  return res;
333  }
334 
335  base_matrix
336  parser_xy_function::hess(scalar_type x, scalar_type y) const {
337  ptx[0] = double(x); // x
338  pty[0] = double(y); // y
339  ptr[0] = double(sqrt(fabs(x*x+y*y))); // r
340  ptt[0] = double(atan2(y,x)); // theta
341 
342  base_matrix res(2,2);
343  const bgeot::base_tensor &t = f_hess.eval();
344  GMM_ASSERT1(t.size() == 4, "Wrong size of expression result "
345  << f_hess.expression());
346  gmm::copy(t.as_vector(), res.as_vector());
347  return res;
348  }
349 
350  /* the basic singular functions for 2D cracks */
351  scalar_type
352  crack_singular_xy_function::val(scalar_type x, scalar_type y) const {
353  scalar_type sgny = (y < 0 ? -1.0 : 1.0);
354  scalar_type r = sqrt(x*x + y*y);
355  if (r < 1e-10) return 0;
356 
357  /* The absolute value is unfortunately necessary, otherwise, sqrt(-1e-16)
358  can be required ...
359  */
360  scalar_type sin2 = sqrt(gmm::abs(.5-x/(2*r))) * sgny;
361  scalar_type cos2 = sqrt(gmm::abs(.5+x/(2*r)));
362 
363  scalar_type res = 0;
364  switch(l){
365 
366  /* First order enrichement displacement field (linear elasticity) */
367 
368  case 0 : res = sqrt(r)*sin2; break;
369  case 1 : res = sqrt(r)*cos2; break;
370  case 2 : res = sin2*y/sqrt(r); break;
371  case 3 : res = cos2*y/sqrt(r); break;
372 
373  /* Second order enrichement of displacement field (linear elasticity) */
374 
375  case 4 : res = sqrt(r)*r*sin2; break;
376  case 5 : res = sqrt(r)*r*cos2; break;
377  case 6 : res = sin2*cos2*cos2*r*sqrt(r); break;
378  case 7 : res = cos2*cos2*cos2*r*sqrt(r); break;
379 
380  /* First order enrichement of pressure field (linear elasticity) mixed formulation */
381 
382  case 8 : res = -sin2/sqrt(r); break;
383  case 9 : res = cos2/sqrt(r); break;
384 
385  /* Second order enrichement of pressure field (linear elasticity) mixed formulation */
386 
387  case 10 : res = sin2*sqrt(r); break; // cos2*cos2
388  case 11 : res = cos2*sqrt(r); break;
389 
390  /* First order enrichement of displacement field (nonlinear elasticity)[Rodney Stephenson Journal of elasticity VOL.12 No. 1, January 1982] */
391 
392  case 12 : res = r*sin2*sin2; break;
393  case 13 : res = sqrt(r)*sin2; break;
394 
395 /* First order enrichement of pressure field (nonlinear elasticity) */
396 
397  case 14 : res = sin2/r; break;
398  case 15 : res = cos2/r; break;
399 
400  default: GMM_ASSERT2(false, "arg");
401  }
402  return res;
403  }
404 
405 
406  base_small_vector
407  crack_singular_xy_function::grad(scalar_type x, scalar_type y) const {
408  scalar_type sgny = (y < 0 ? -1.0 : 1.0);
409  scalar_type r = sqrt(x*x + y*y);
410 
411  if (r < 1e-10) {
412  GMM_WARNING0("Warning, point close to the singularity (r=" << r << ")");
413  }
414 
415  /* The absolute value is unfortunately necessary, otherwise, sqrt(-1e-16)
416  can be required ...
417  */
418  scalar_type sin2 = sqrt(gmm::abs(.5-x/(2*r))) * sgny;
419  scalar_type cos2 = sqrt(gmm::abs(.5+x/(2*r)));
420 
421  base_small_vector res(2);
422  switch(l){
423  /* First order enrichement displacement field (linear elasticity) */
424  case 0 :
425  res[0] = -sin2/(2*sqrt(r));
426  res[1] = cos2/(2*sqrt(r));
427  break;
428  case 1 :
429  res[0] = cos2/(2*sqrt(r));
430  res[1] = sin2/(2*sqrt(r));
431  break;
432  case 2 :
433  res[0] = cos2*((-5*cos2*cos2) + 1. + 4*(cos2*cos2*cos2*cos2))/sqrt(r);
434  res[1] = sin2*((-3*cos2*cos2) + 1. + 4*(cos2*cos2*cos2*cos2))/sqrt(r);
435  break;
436  case 3 :
437  res[0] = -cos2*cos2*sin2*((4*cos2*cos2) - 3.)/sqrt(r);
438  res[1] = cos2*((4*cos2*cos2*cos2*cos2) + 2. - (5*cos2*cos2))/sqrt(r);
439  break;
440  /* Second order enrichement of displacement field (linear elasticity) */
441 
442  case 4 :
443  res[0] = sin2 *((4*cos2*cos2)-3.)*sqrt(r)/2.;
444  res[1] = cos2*(5. - (4*cos2*cos2))*sqrt(r)/2. ;
445  break;
446  case 5 :
447  res[0] = cos2*((4*cos2*cos2)-1.)*sqrt(r)/2.;
448  res[1] = sin2*((4*cos2*cos2)+1.)*sqrt(r)/2. ;
449  break;
450 
451  case 6 :
452  res[0] = sin2*cos2*cos2*sqrt(r)/2.;
453  res[1] = cos2*(2. - (cos2*cos2))*sqrt(r)/2.;
454  break;
455  case 7 :
456  res[0] = 3*cos2*cos2*cos2*sqrt(r)/2.;
457  res[1] = 3*sin2*cos2*cos2*sqrt(r)/2.;
458  break;
459 
460  /* First order enrichement of pressure field (linear elasticity) mixed formulation */
461 
462  case 8 :
463  res[0] =sin2*((4*cos2*cos2)-1.)/(2*sqrt(r)*r);
464  res[1] =-cos2*((4*cos2*cos2)-3.)/(2*sqrt(r)*r);
465  break;
466  case 9 :
467  res[0] =-cos2*((2*cos2*cos2) - 3.)/(2*sqrt(r)*r);
468  res[1] =-sin2*((4*cos2*cos2)-1.)/(2*sqrt(r)*r);
469  break;
470 
471  /* Second order enrichement of pressure field (linear elasticity) mixed formulation */
472  case 10 :
473  res[0] = -sin2/(2*sqrt(r));
474  res[1] = cos2/(2*sqrt(r));
475  break;
476  case 11 :
477  res[0] = cos2/(2*sqrt(r));
478  res[1] = sin2/(2*sqrt(r));
479  break;
480 
481  /* First order enrichement of displacement field (nonlinear elasticity)[Rodney Stephenson Journal of elasticity VOL.12 No. 1, January 1982] */
482 
483  case 12 :
484  res[0] = sin2*sin2;
485  res[1] = 0.5*cos2*sin2;
486  break;
487  case 13 :
488  res[0] = -sin2/(2*sqrt(r));
489  res[1] = cos2/(2*sqrt(r));
490  break;
491 
492 /* First order enrichement of pressure field (****Nonlinear elasticity*****) */
493 
494 
495  case 14 :
496  res[0] = -sin2/r;
497  res[1] = cos2/(2*r);
498  break;
499  case 15 :
500  res[0] = -cos2/r;
501  res[1] = -sin2/(2*r);
502  break;
503 
504 
505  default: GMM_ASSERT2(false, "oups");
506  }
507  return res;
508  }
509 
510  base_matrix
511  crack_singular_xy_function::hess(scalar_type x, scalar_type y) const {
512  scalar_type sgny = (y < 0 ? -1.0 : 1.0);
513  scalar_type r = sqrt(x*x + y*y);
514 
515  if (r < 1e-10) {
516  GMM_WARNING0("Warning, point close to the singularity (r=" << r << ")");
517  }
518 
519  /* The absolute value is unfortunately necessary, otherwise, sqrt(-1e-16)
520  can be required ...
521  */
522  scalar_type sin2 = sqrt(gmm::abs(.5-x/(2*r))) * sgny;
523  scalar_type cos2 = sqrt(gmm::abs(.5+x/(2*r)));
524 
525  base_matrix res(2,2);
526  switch(l){
527  case 0 :
528  res(0,0) = (-sin2*x*x + 2.0*cos2*x*y + sin2*y*y) / (4*pow(r, 3.5));
529  res(0,1) = (-cos2*x*x - 2.0*sin2*x*y + cos2*y*y) / (4*pow(r, 3.5));
530  res(1,0) = res(0, 1);
531  res(1,1) = (sin2*x*x - 2.0*cos2*x*y - sin2*y*y) / (4*pow(r, 3.5));
532  break;
533  case 1 :
534  res(0,0) = (-cos2*x*x - 2.0*sin2*x*y + cos2*y*y) / (4*pow(r, 3.5));
535  res(0,1) = (sin2*x*x - 2.0*cos2*x*y - sin2*y*y) / (4*pow(r, 3.5));
536  res(1,0) = res(0, 1);
537  res(1,1) = (cos2*x*x + 2.0*sin2*x*y - cos2*y*y) / (4*pow(r, 3.5));
538  break;
539  case 2 :
540  res(0,0) = 3.0*y*(sin2*x*x + 2.0*cos2*x*y - sin2*y*y) / (4*pow(r, 4.5));
541  res(0,1) = (-2.0*sin2*x*x*x - 5.0*cos2*y*x*x + 4.0*sin2*y*y*x + cos2*y*y*y)
542  / (4*pow(r, 4.5));
543  res(1,0) = res(0, 1);
544  res(1,1) = (4.0*cos2*x*x*x - 7.0*sin2*y*x*x - 2.0*cos2*y*y*x - sin2*y*y*y)
545  / (4*pow(r, 4.5));
546  break;
547  case 3 :
548  res(0,0) = 3.0*y*(cos2*x*x - 2.0*sin2*x*y - cos2*y*y) / (4*pow(r, 4.5));
549  res(0,1) = (-2.0*cos2*x*x*x + 5.0*sin2*y*x*x + 4.0*cos2*y*y*x - sin2*y*y*y)
550  / (4*pow(r, 4.5));
551  res(1,0) = res(0, 1);
552  res(1,1) = (-4.0*sin2*x*x*x - 7.0*cos2*y*x*x + 2.0*sin2*y*y*x - cos2*y*y*y)
553  / (4*pow(r, 4.5));
554  break;
555  default: GMM_ASSERT2(false, "oups");
556  }
557  return res;
558  }
559 
560  scalar_type
561  cutoff_xy_function::val(scalar_type x, scalar_type y) const {
562  scalar_type res = 1;
563  switch (fun) {
564  case EXPONENTIAL_CUTOFF: {
565  res = (a4>0) ? exp(-a4 * gmm::sqr(x*x+y*y)) : 1;
566  break;
567  }
568  case POLYNOMIAL_CUTOFF: {
569  assert(r0 > r1);
570  scalar_type r = gmm::sqrt(x*x+y*y);
571 
572  if (r <= r1) res = 1;
573  else if (r >= r0) res = 0;
574  else {
575  // scalar_type c = 6./(r0*r0*r0 - r1*r1*r1 + 3*r1*r0*(r1-r0));
576  // scalar_type k = -(c/6.)*(-pow(r0,3.) + 3*r1*pow(r0,2.));
577  // res = (c/3.)*pow(r,3.) - (c*(r0 + r1)/2.)*pow(r,2.) +
578  // c*r0*r1*r + k;
579  scalar_type c = 1./pow(r0-r1,3.0);
580  res = c*(r*(r*(2.0*r-3.0*(r0+r1))+6.0*r1*r0) + r0*r0*(r0-3.0*r1));
581  }
582  break;
583  }
584  case POLYNOMIAL2_CUTOFF: {
585  assert(r0 > r1);
586  scalar_type r = gmm::sqrt(x*x+y*y);
587  if (r <= r1) res = scalar_type(1);
588  else if (r >= r0) res = scalar_type(0);
589  else {
590 // scalar_type c = 1./((-1./30.)*(pow(r1,5) - pow(r0,5))
591 // + (1./6.)*(pow(r1,4)*r0 - r1*pow(r0,4))
592 // - (1./3.)*(pow(r1,3)*pow(r0,2) -
593 // pow(r1,2)*pow(r0,3)));
594 // scalar_type k = 1. - c*((-1./30.)*pow(r1,5) +
595 // (1./6.)*pow(r1,4)*r0 -
596 // (1./3.)*pow(r1,3)*pow(r0,2));
597 // res = c*( (-1./5.)*pow(r,5) + (1./2.)* (r1+r0)*pow(r,4) -
598 // (1./3.)*(pow(r1,2)+pow(r0,2) + 4.*r0*r1)*pow(r,3) +
599 // r0*r1*(r0+r1)*pow(r,2) - pow(r0,2)*pow(r1,2)*r) + k;
600  res = (r*(r*(r*(r*(-6.0*r + 15.0*(r0+r1))
601  - 10.0*(r0*r0 + 4.0*r1*r0 + r1*r1))
602  + 30.0 * r0*r1*(r0+r1)) - 30.0*r1*r1*r0*r0)
603  + r0*r0*r0*(r0*r0-5.0*r1*r0+10*r1*r1)) / pow(r0-r1, 5.0);
604  }
605  break;
606  }
607  default : res = 1;
608  }
609  return res;
610  }
611 
612  base_small_vector
613  cutoff_xy_function::grad(scalar_type x, scalar_type y) const {
614  base_small_vector res(2);
615  switch (fun) {
616  case EXPONENTIAL_CUTOFF: {
617  scalar_type r2 = x*x+y*y, ratio = -4.*exp(-a4*r2*r2)*a4*r2;
618  res[0] = ratio*x;
619  res[1] = ratio*y;
620  break;
621  }
622  case POLYNOMIAL_CUTOFF: {
623  scalar_type r = gmm::sqrt(x*x+y*y);
624  scalar_type ratio = 0;
625 
626  if ( r > r1 && r < r0 ) {
627  // scalar_type c = 6./(pow(r0,3.) - pow(r1,3.) + 3*r1*r0*(r1-r0));
628  // ratio = c*(r - r0)*(r - r1);
629  ratio = 6.*(r - r0)*(r - r1)/pow(r0-r1, 3.);
630  }
631  res[0] = ratio*x/r;
632  res[1] = ratio*y/r;
633  break;
634  }
635  case POLYNOMIAL2_CUTOFF: {
636  scalar_type r = gmm::sqrt(x*x+y*y);
637  scalar_type ratio = 0;
638  if (r > r1 && r < r0) {
639 // scalar_type c = 1./((-1./30.)*(pow(r1,5) - pow(r0,5))
640 // + (1./6.)*(pow(r1,4)*r0 - r1*pow(r0,4))
641 // - (1./3.)*(pow(r1,3)*pow(r0,2)
642 // - pow(r1,2)*pow(r0,3)));
643 // ratio = - c*gmm::sqr(r-r0)*gmm::sqr(r-r1);
644  ratio = -30.0*gmm::sqr(r-r0)*gmm::sqr(r-r1) / pow(r0-r1, 5.0);
645  }
646  res[0] = ratio*x/r;
647  res[1] = ratio*y/r;
648  break;
649  }
650  default :
651  res[0] = 0;
652  res[1] = 0;
653  }
654  return res;
655  }
656 
657  base_matrix
658  cutoff_xy_function::hess(scalar_type x, scalar_type y) const {
659  base_matrix res(2,2);
660  switch (fun) {
661  case EXPONENTIAL_CUTOFF: {
662  scalar_type r2 = x*x+y*y, r4 = r2*r2;
663  res(0,0) = 4.0*a4*(-3.0*x*x - y*y + 4.0*a4*x*x*r4)*exp(-a4*r4);
664  res(1,0) = 8.0*a4*x*y*(-1.0 + 2.0*a4*r4)*exp(-a4*r4);
665  res(0,1) = res(1,0);
666  res(1,1) = 4.0*a4*(-3.0*y*y - x*x + 4.0*a4*y*y*r4)*exp(-a4*r4);
667  break;
668  }
669  case POLYNOMIAL_CUTOFF: {
670  scalar_type r2 = x*x+y*y, r = gmm::sqrt(r2), c=6./(pow(r0-r1,3.)*r*r2);
671  if ( r > r1 && r < r0 ) {
672  res(0,0) = c*(x*x*r2 + r1*r0*y*y - r*r2*(r0+r1) + r2*r2);
673  res(1,0) = c*x*y*(r2 - r1*r0);
674  res(0,1) = res(1,0);
675  res(1,1) = c*(y*y*r2 + r1*r0*x*x - r*r2*(r0+r1) + r2*r2);
676  }
677  break;
678  }
679  case POLYNOMIAL2_CUTOFF: {
680  scalar_type r2 = x*x+y*y, r = gmm::sqrt(r2), r3 = r*r2;
681  if (r > r1 && r < r0) {
682  scalar_type dp = -30.0*(r1-r)*(r1-r)*(r0-r)*(r0-r) / pow(r0-r1, 5.0);
683  scalar_type ddp = 60.0*(r1-r)*(r0-r)*(r0+r1-2.0*r) / pow(r0-r1, 5.0);
684  scalar_type rx= x/r, ry= y/r, rxx= y*y/r3, rxy= -x*y/r3, ryy= x*x/r3;
685  res(0,0) = ddp*rx*rx + dp*rxx;
686  res(1,0) = ddp*rx*ry + dp*rxy;
687  res(0,1) = res(1,0);
688  res(1,1) = ddp*ry*ry + dp*ryy;
689  }
690  break;
691  }
692  }
693  return res;
694  }
695 
696  cutoff_xy_function::cutoff_xy_function(int fun_num, scalar_type r,
697  scalar_type r1_, scalar_type r0_)
698  {
699  fun = fun_num;
700  r1 = r1_; r0 = r0_;
701  a4 = 0;
702  if (r > 0) a4 = pow(2.7/r,4.);
703  }
704 
705 
706  struct global_function_on_levelsets_2D_ :
707  public global_function, public context_dependencies {
708  const std::vector<level_set> dummy_lsets;
709  const std::vector<level_set> &lsets;
710  const level_set &ls;
711  mutable pmesher_signed_distance mls_x, mls_y;
712  mutable size_type cv;
713 
714  pxy_function fn;
715 
716  void update_mls(size_type cv_, size_type n) const {
717  if (cv_ != cv) {
718  cv=cv_;
719  if (lsets.size() == 0) {
720  mls_x = ls.mls_of_convex(cv, 1);
721  mls_y = ls.mls_of_convex(cv, 0);
722  } else {
723  base_node pt(n);
724  scalar_type d = scalar_type(-2);
725  for (const auto &ls_ : lsets) {
726  pmesher_signed_distance mls_xx, mls_yy;
727  mls_xx = ls_.mls_of_convex(cv, 1);
728  mls_yy = ls_.mls_of_convex(cv, 0);
729  scalar_type x = (*mls_xx)(pt), y = (*mls_yy)(pt);
730  scalar_type d2 = gmm::sqr(x) + gmm::sqr(y);
731  if (d < scalar_type(-1) || d2 < d) {
732  d = d2;
733  mls_x = mls_xx;
734  mls_y = mls_yy;
735  }
736  }
737  }
738  }
739  }
740 
741  virtual scalar_type val(const fem_interpolation_context& c) const {
742  update_mls(c.convex_num(), c.xref().size());
743  scalar_type x = (*mls_x)(c.xref());
744  scalar_type y = (*mls_y)(c.xref());
745  if (c.xfem_side() > 0 && y <= 1E-13) y = 1E-13;
746  if (c.xfem_side() < 0 && y >= -1E-13) y = -1E-13;
747  return fn->val(x,y);
748  }
749  virtual void grad(const fem_interpolation_context& c,
750  base_small_vector &g) const {
751  size_type P = c.xref().size();
752  base_small_vector dx(P), dy(P), dfr(2);
753 
754  update_mls(c.convex_num(), P);
755  scalar_type x = mls_x->grad(c.xref(), dx);
756  scalar_type y = mls_y->grad(c.xref(), dy);
757  if (c.xfem_side() > 0 && y <= 0) y = 1E-13;
758  if (c.xfem_side() < 0 && y >= 0) y = -1E-13;
759 
760  base_small_vector gfn = fn->grad(x,y);
761  gmm::mult(c.B(), gfn[0]*dx + gfn[1]*dy, g);
762  }
763  virtual void hess(const fem_interpolation_context&c,
764  base_matrix &h) const {
765  size_type P = c.xref().size(), N = c.N();
766  base_small_vector dx(P), dy(P), dfr(2), dx_real(N), dy_real(N);
767 
768  update_mls(c.convex_num(), P);
769  scalar_type x = mls_x->grad(c.xref(), dx);
770  scalar_type y = mls_y->grad(c.xref(), dy);
771  if (c.xfem_side() > 0 && y <= 0) y = 1E-13;
772  if (c.xfem_side() < 0 && y >= 0) y = -1E-13;
773 
774  base_small_vector gfn = fn->grad(x,y);
775  base_matrix hfn = fn->hess(x,y);
776 
777  base_matrix hx, hy, hx_real(N*N, 1), hy_real(N*N, 1);
778  mls_x->hess(c.xref(), hx);
779  mls_x->hess(c.xref(), hy);
780  gmm::reshape(hx, P*P, 1);
781  gmm::reshape(hy, P*P, 1);
782 
783  gmm::mult(c.B3(), hx, hx_real);
784  gmm::mult(c.B32(), gmm::scaled(dx, -1.0), gmm::mat_col(hx_real, 0));
785  gmm::mult(c.B3(), hy, hy_real);
786  gmm::mult(c.B32(), gmm::scaled(dy, -1.0), gmm::mat_col(hy_real, 0));
787  gmm::mult(c.B(), dx, dx_real);
788  gmm::mult(c.B(), dy, dy_real);
789 
790  for (size_type i = 0; i < N; ++i)
791  for (size_type j = 0; j < N; ++j) {
792  h(i, j) = hfn(0,0) * dx_real[i] * dx_real[j]
793  + hfn(0,1) * dx_real[i] * dy_real[j]
794  + hfn(1,0) * dy_real[i] * dx_real[j]
795  + hfn(1,1) * dy_real[i] * dy_real[j]
796  + gfn[0] * hx_real(i * N + j, 0) + gfn[1] * hy_real(i*N+j,0);
797  }
798  }
799 
800  void update_from_context() const { cv = size_type(-1); }
801 
802  global_function_on_levelsets_2D_(const std::vector<level_set> &lsets_,
803  const pxy_function &fn_)
804  : global_function(2), dummy_lsets(0, dummy_level_set()),
805  lsets(lsets_), ls(dummy_level_set()), fn(fn_) {
806  GMM_ASSERT1(lsets.size() > 0, "The list of level sets should"
807  " contain at least one level set.");
808  cv = size_type(-1);
809  for (size_type i = 0; i < lsets.size(); ++i)
810  this->add_dependency(lsets[i]);
811  }
812 
813  global_function_on_levelsets_2D_(const level_set &ls_,
814  const pxy_function &fn_)
815  : global_function(2), dummy_lsets(0, dummy_level_set()),
816  lsets(dummy_lsets), ls(ls_), fn(fn_) {
817  cv = size_type(-1);
818  }
819 
820  };
821 
822  pglobal_function
823  global_function_on_level_sets(const std::vector<level_set> &lsets,
824  const pxy_function &fn) {
825  return std::make_shared<global_function_on_levelsets_2D_>(lsets, fn);
826  }
827 
828  pglobal_function
829  global_function_on_level_set(const level_set &ls,
830  const pxy_function &fn) {
831  return std::make_shared<global_function_on_levelsets_2D_>(ls, fn);
832  }
833 
834  // interpolator on mesh_fem
835 
836  interpolator_on_mesh_fem::interpolator_on_mesh_fem
837  (const mesh_fem &mf_, const std::vector<scalar_type> &U_)
838  : mf(mf_), U(U_) {
839 
840  if (mf.is_reduced()) {
841  gmm::resize(U, mf.nb_basic_dof());
842  gmm::mult(mf.extension_matrix(), U_, U);
843  }
844  base_node bmin, bmax;
845  scalar_type EPS=1E-13;
846  cv_stored = size_type(-1);
847  boxtree.clear();
848  for (dal::bv_visitor cv(mf.convex_index()); !cv.finished(); ++cv) {
849  bounding_box(bmin, bmax, mf.linked_mesh().points_of_convex(cv),
850  mf.linked_mesh().trans_of_convex(cv));
851  for (auto&& val : bmin) val -= EPS;
852  for (auto&& val : bmax) val += EPS;
853  boxtree.add_box(bmin, bmax, cv);
854  }
855  boxtree.build_tree();
856  }
857 
858  bool interpolator_on_mesh_fem::find_a_point(const base_node &pt,
859  base_node &ptr,
860  size_type &cv) const {
861  bool gt_invertible;
862  if (cv_stored != size_type(-1) && gic.invert(pt, ptr, gt_invertible)) {
863  cv = cv_stored;
864  if (gt_invertible)
865  return true;
866  }
867 
868  boxtree.find_boxes_at_point(pt, boxlst);
869  for (const auto &box : boxlst) {
870  gic = bgeot::geotrans_inv_convex
871  (mf.linked_mesh().convex(box->id),
872  mf.linked_mesh().trans_of_convex(box->id));
873  cv_stored = box->id;
874  if (gic.invert(pt, ptr, gt_invertible)) {
875  cv = box->id;
876  return true;
877  }
878  }
879  return false;
880  }
881 
882  bool interpolator_on_mesh_fem::eval(const base_node &pt, base_vector &val,
883  base_matrix &grad) const {
884  base_node ptref;
885  size_type cv;
886  base_vector coeff;
887  dim_type q = mf.get_qdim(), N = mf.linked_mesh().dim();
888  if (find_a_point(pt, ptref, cv)) {
889  pfem pf = mf.fem_of_element(cv);
890  bgeot::pgeometric_trans pgt =
891  mf.linked_mesh().trans_of_convex(cv);
892  base_matrix G;
893  vectors_to_base_matrix(G, mf.linked_mesh().points_of_convex(cv));
894  fem_interpolation_context fic(pgt, pf, ptref, G, cv, short_type(-1));
895  slice_vector_on_basic_dof_of_element(mf, U, cv, coeff);
896  // coeff.resize(mf.nb_basic_dof_of_element(cv));
897  // gmm::copy(gmm::sub_vector
898  // (U,gmm::sub_index(mf.ind_basic_dof_of_element(cv))), coeff);
899  val.resize(q);
900  pf->interpolation(fic, coeff, val, q);
901  grad.resize(q, N);
902  pf->interpolation_grad(fic, coeff, grad, q);
903  return true;
904  } else
905  return false;
906  }
907 }
908 
909 /* end of namespace getfem */
bgeot::geotrans_inv_convex
does the inversion of the geometric transformation for a given convex
Definition: bgeot_geotrans_inv.h:64
gmm::resize
void resize(M &v, size_type m, size_type n)
*‍/
Definition: gmm_blas.h:231
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_global_function.h
Definition of global functions to be used as base or enrichment functions in fem.
bgeot::short_type
gmm::uint16_type short_type
used as the common short type integer in the library
Definition: bgeot_config.h:72
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:46
gmm::reshape
void reshape(M &v, size_type m, size_type n)
*‍/
Definition: gmm_blas.h:250
getfem::pfem
std::shared_ptr< const getfem::virtual_fem > pfem
type of pointer on a fem description
Definition: getfem_fem.h:244
getfem_level_set.h
Define level-sets.
getfem::slice_vector_on_basic_dof_of_element
void slice_vector_on_basic_dof_of_element(const mesh_fem &mf, const VEC1 &vec, size_type cv, VEC2 &coeff, size_type qmult1=size_type(-1), size_type qmult2=size_type(-1))
Given a mesh_fem.
Definition: getfem_mesh_fem.h:659
bgeot::pgeometric_trans
std::shared_ptr< const bgeot::geometric_trans > pgeometric_trans
pointer type for a geometric transformation
Definition: bgeot_geometric_trans.h:186
dal::add_dependency
void add_dependency(pstatic_stored_object o1, pstatic_stored_object o2)
Add a dependency, object o1 will depend on object o2.
Definition: dal_static_stored_objects.cc:230