GetFEM  5.4.3
getfem_global_function.cc
1 /*===========================================================================
2 
3  Copyright (C) 2004-2022 Yves Renard
4  Copyright (C) 2016-2022 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 level_set &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 (const level_set &ls_ : lsets)
810  this->add_dependency(ls_);
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  this->add_dependency(ls);
819  }
820 
821  };
822 
823  pglobal_function
824  global_function_on_level_sets(const std::vector<level_set> &lsets,
825  const pxy_function &fn) {
826  return std::make_shared<global_function_on_levelsets_2D_>(lsets, fn);
827  }
828 
829  pglobal_function
830  global_function_on_level_set(const level_set &ls,
831  const pxy_function &fn) {
832  return std::make_shared<global_function_on_levelsets_2D_>(ls, fn);
833  }
834 
835 
836 
837 
838  // Global function for bspline basis
839  const scalar_type eps(1e-12);
840 
841  // order k = 3
842  scalar_type bsp3_01(scalar_type t) {
843  return (t >= -eps && t < 1) ? pow(1.-t,2)
844  : 0;
845  }
846  scalar_type bsp3_01_der(scalar_type t) {
847  return (t >= -eps && t < 1) ? 2.*t-2.
848  : 0;
849  }
850  scalar_type bsp3_01_der2(scalar_type t) {
851  return (t >= eps && t <= 1.-eps) ? 2.
852  : 0;
853  }
854  scalar_type bsp3_01_der2_with_hint(scalar_type t, scalar_type t_hint) {
855  return (t > -eps && t < 1.+eps && t_hint > 0 && t_hint < 1) ? 2.
856  : 0;
857  }
858  scalar_type bsp3_02(scalar_type t) {
859  if (t >= -eps) {
860  if (t < 1)
861  return (-1.5*t+2.)*t;
862  else if (t < 2)
863  return 0.5*pow(2.-t,2);
864  }
865  return 0;
866  }
867  scalar_type bsp3_02_der(scalar_type t) {
868  if (t >= -eps) {
869  if (t < 1)
870  return -3.*t+2.;
871  else if (t < 2)
872  return t-2.;
873  }
874  return 0;
875  }
876  scalar_type bsp3_02_der2(scalar_type t) {
877  if (t >= eps) {
878  if (t < 1.-eps)
879  return -3.;
880  else if (t < 1.+eps)
881  return 0;
882  else if (t <= 2.-eps)
883  return 1.;
884  }
885  return 0;
886  }
887  scalar_type bsp3_03(scalar_type t) {
888  if (t >= -eps) {
889  if (t < 1) {
890  return 0.5*t*t;
891  } else if (t < 2) {
892  t -= 1.;
893  return t*(1-t)+0.5;
894  } else if (t < 3) {
895  t = 3.-t;
896  return 0.5*t*t;
897  }
898  }
899  return 0;
900  }
901  scalar_type bsp3_03_der(scalar_type t) {
902  if (t >= -eps) {
903  if (t < 1) {
904  return t;
905  } else if (t < 2) {
906  t -= 1.;
907  return 1.-2.*t;
908  } else if (t < 3) {
909  return t-3.;
910  }
911  }
912  return 0;
913  }
914  scalar_type bsp3_03_der2(scalar_type t) {
915  if (t >= -eps) {
916  if (t < eps)
917  return 0;
918  else if (t <= 1.-eps)
919  return 1.;
920  else if (t < 1.+eps)
921  return 0;
922  else if (t <= 2.-eps)
923  return -2.;
924  else if (t < 2.+eps)
925  return 0;
926  else if (t <= 3.-eps)
927  return 1.;
928  else if (t < 3.+eps)
929  return 0;
930  }
931  return 0;
932  }
933 
934  // order k = 4
935  scalar_type bsp4_01(scalar_type t) {
936  return (t >= -eps && t < 1) ? pow(1.-t,3)
937  : 0;
938  }
939  scalar_type bsp4_01_der(scalar_type t) {
940  return (t >= -eps && t < 1) ? -3*pow(1.-t,2)
941  : 0;
942  }
943  scalar_type bsp4_01_der2(scalar_type t) {
944  return (t >= -eps && t < 1) ? 6*(1.-t)
945  : 0;
946  }
947  scalar_type bsp4_02(scalar_type t) {
948  if (t >= -eps) {
949  if (t < 1) {
950  return ((7./4.*t-9./2.)*t+3.)*t;
951  } else if (t < 2) {
952  return 1./4.*pow(2.-t,3);
953  }
954  }
955  return 0;
956  }
957  scalar_type bsp4_02_der(scalar_type t) {
958  if (t >= -eps) {
959  if (t < 1) {
960  return (21./4.*t-9.)*t+3.;
961  } else if (t < 2) {
962  return -3./4.*pow(2.-t,2);
963  }
964  }
965  return 0;
966  }
967  scalar_type bsp4_02_der2(scalar_type t) {
968  if (t >= -eps) {
969  if (t < 1) {
970  return 21./2.*t-9.;
971  } else if (t < 2) {
972  return 3./2.*(2.-t);
973  }
974  }
975  return 0;
976  }
977  scalar_type bsp4_03(scalar_type t) {
978  if (t >= -eps) {
979  if (t < 1) {
980  return (-11./12.*t+3./2.)*t*t;
981  } else if (t < 2) {
982  t -= 1;
983  return ((7./12.*t-5./4.)*t+1./4.)*t+7./12.;
984  } else if (t < 3) {
985  return 1./6.*pow(3.-t,3);
986  }
987  }
988  return 0;
989  }
990  scalar_type bsp4_03_der(scalar_type t) {
991  if (t >= -eps) {
992  if (t < 1) {
993  return (-11./4.*t+3.)*t;
994  } else if (t < 2) {
995  t -= 1;
996  return (7./4.*t-5./2.)*t+1./4.;
997  } else if (t < 3) {
998  return -1./2.*pow(3.-t,2);
999  }
1000  }
1001  return 0;
1002  }
1003  scalar_type bsp4_03_der2(scalar_type t) {
1004  if (t >= -eps) {
1005  if (t < 1) {
1006  return -11./2.*t+3.;
1007  } else if (t < 2) {
1008  t -= 1;
1009  return 7./2.*t-5./2.;
1010  } else if (t < 3) {
1011  return 3.-t;
1012  }
1013  }
1014  return 0;
1015  }
1016  scalar_type bsp4_04(scalar_type t) {
1017  if (t > 2) {
1018  t = 4.-t;
1019  }
1020  if (t >= -eps) {
1021  if (t < 1) {
1022  return 1./6.*pow(t,3);
1023  } else if (t <= 2) {
1024  t = t-1.;
1025  return ((-1./2.*t+1./2.)*t+1./2.)*t+1./6.;
1026  }
1027  }
1028  return 0;
1029  }
1030  scalar_type bsp4_04_der(scalar_type t) {
1031  scalar_type sgn(1);
1032  if (t > 2) {
1033  t = 4.-t;
1034  sgn = -1.;
1035  }
1036  if (t >= -eps) {
1037  if (t < 1) {
1038  return 1./2.*pow(t,2)*sgn;
1039  } else if (t < 2) {
1040  t = t-1.;
1041  return ((-3./2.*t+1.)*t+1./2.)*sgn;
1042  }
1043  }
1044  return 0;
1045  }
1046  scalar_type bsp4_04_der2(scalar_type t) {
1047  if (t > 2) {
1048  t = 4.-t;
1049  }
1050  if (t >= -eps) {
1051  if (t < 1) {
1052  return t;
1053  } else if (t < 2) {
1054  t = t-1.;
1055  return -3.*t+1.;
1056  }
1057  }
1058  return 0;
1059  }
1060 
1061 
1062  // order k = 5
1063  scalar_type bsp5_01(scalar_type t) {
1064  return (t >= -eps && t < 1) ? pow(1.-t,4)
1065  : 0;
1066  }
1067  scalar_type bsp5_01_der(scalar_type t) {
1068  return (t >= -eps && t < 1) ? -4.*pow(1.-t,3)
1069  : 0;
1070  }
1071  scalar_type bsp5_01_der2(scalar_type t) {
1072  return (t >= -eps && t < 1) ? 12.*pow(1.-t,2)
1073  : 0;
1074  }
1075  scalar_type bsp5_02(scalar_type t) {
1076  if (t >= -eps) {
1077  if (t < 1) {
1078  return (((-15./8.*t+7.)*t-9.)*t+4.)*t;
1079  } else if (t < 2) {
1080  return 1./8.*pow(2.-t,4);
1081  }
1082  }
1083  return 0;
1084  }
1085  scalar_type bsp5_02_der(scalar_type t) {
1086  if (t >= -eps) {
1087  if (t < 1) {
1088  return ((-15./2.*t+21.)*t-18.)*t+4.;
1089  } else if (t < 2) {
1090  return -1./2.*pow(2.-t,3);
1091  }
1092  }
1093  return 0;
1094  }
1095  scalar_type bsp5_02_der2(scalar_type t) {
1096  if (t >= -eps) {
1097  if (t < 1) {
1098  return (-45./2.*t+42.)*t-18.;
1099  } else if (t < 2) {
1100  return 3./2.*pow(2.-t,2);
1101  }
1102  }
1103  return 0;
1104  }
1105  scalar_type bsp5_03(scalar_type t) {
1106  if (t >= -eps) {
1107  if (t < 1) {
1108  return ((85./72.*t-11./3.)*t+3.)*t*t;
1109  } else if (t < 2) {
1110  t = 2.-t;
1111  return (((-23./72*t+2./9.)*t+1./3.)*t+2./9.)*t+1./18.;
1112  } else if (t < 3) {
1113  return 1./18.*pow(3.-t,4);
1114  }
1115  }
1116  return 0;
1117  }
1118  scalar_type bsp5_03_der(scalar_type t) {
1119  if (t >= -eps) {
1120  if (t < 1) {
1121  return ((85./18.*t-11.)*t+6.)*t;
1122  } else if (t < 2) {
1123  t = 2.-t;
1124  return -(((-23./18*t+2./3.)*t+2./3.)*t+2./9.);
1125  } else if (t < 3) {
1126  return -2./9.*pow(3.-t,3);
1127  }
1128  }
1129  return 0;
1130  }
1131  scalar_type bsp5_03_der2(scalar_type t) {
1132  if (t >= -eps) {
1133  if (t < 1) {
1134  return (85./6.*t-22.)*t+6.;
1135  } else if (t < 2) {
1136  t = 2.-t;
1137  return (-23./6*t+4./3.)*t+2./3.;
1138  } else if (t < 3) {
1139  return 2./3.*pow(3.-t,2);
1140  }
1141  }
1142  return 0;
1143  }
1144  scalar_type bsp5_04(scalar_type t) {
1145  if (t >= -eps) {
1146  if (t < 1) {
1147  return (-25./72.*t+2./3.)*t*t*t;
1148  } else if (t < 2) {
1149  t = 2.-t;
1150  return (((23./72.*t-5./9.)*t-1./3.)*t+4./9.)*t+4./9.;
1151  } else if (t < 3) {
1152  t = t-2.;
1153  return (((-13./72.*t+5./9.)*t-1./3.)*t-4./9.)*t+4./9.;
1154  } else if (t < 4) {
1155  return 1./24.*pow(4.-t,4);
1156  }
1157  }
1158  return 0;
1159  }
1160  scalar_type bsp5_04_der(scalar_type t) {
1161  if (t >= -eps) {
1162  if (t < 1) {
1163  return (-25./18.*t+2.)*t*t;
1164  } else if (t < 2) {
1165  t = 2.-t;
1166  return -(((23./18.*t-5./3.)*t-2./3.)*t+4./9.);
1167  } else if (t < 3) {
1168  t = t-2.;
1169  return ((-13./18.*t+5./3.)*t-2./3.)*t-4./9.;
1170  } else if (t < 4) {
1171  return -1./6.*pow(4.-t,3);
1172  }
1173  }
1174  return 0;
1175  }
1176  scalar_type bsp5_04_der2(scalar_type t) {
1177  if (t >= -eps) {
1178  if (t < 1) {
1179  return (-25./6.*t+4.)*t;
1180  } else if (t < 2) {
1181  t = 2.-t;
1182  return (23./6.*t-10./3.)*t-2./3.;
1183  } else if (t < 3) {
1184  t = t-2.;
1185  return -((-13./6.*t+10./3.)*t-2./3.);
1186  } else if (t < 4) {
1187  return 1./2.*pow(4.-t,2);
1188  }
1189  }
1190  return 0;
1191  }
1192  scalar_type bsp5_05(scalar_type t) {
1193  if (t > scalar_type(2.5)) {
1194  t = 5.-t;
1195  }
1196  if (t >= -eps) {
1197  if (t < 1) {
1198  return 1./24.*pow(t,4);
1199  } else if (t < 2) {
1200  t = t-1.;
1201  return (((-1./6.*t+1./6.)*t+1./4.)*t+1./6.)*t+1./24.;
1202  } else if (t < 3) {
1203  t = pow(t-2.5,2);
1204  return (1./4.*t-5./8.)*t+115./192.;
1205  }
1206  }
1207  return 0;
1208  }
1209  scalar_type bsp5_05_der(scalar_type t) {
1210  scalar_type sgn(1);
1211  if (t > scalar_type(2.5)) {
1212  t = 5.-t;
1213  sgn = -1.;
1214  }
1215  if (t >= -eps) {
1216  if (t < 1) {
1217  return 1./6.*pow(t,3)*sgn;
1218  } else if (t < 2) {
1219  t = t-1.;
1220  return (((-2./3.*t+1./2.)*t+1./2)*t+1./6.)*sgn;
1221  } else if (t < 3) {
1222  t = t-2.5;
1223  return (t*t-5./4.)*t*sgn;
1224  }
1225  }
1226  return 0;
1227  }
1228  scalar_type bsp5_05_der2(scalar_type t) {
1229  if (t > scalar_type(2.5)) {
1230  t = 5.-t;
1231  }
1232  if (t >= -eps) {
1233  if (t < 1) {
1234  return 1./2.*pow(t,2);
1235  } else if (t < 2) {
1236  t = t-1.;
1237  return (-2.*t+1.)*t+1./2;
1238  } else if (t < 3) {
1239  t = t-2.5;
1240  return 3*t*t-5./4.;
1241  }
1242  }
1243  return 0;
1244  }
1245 
1246 
1247 
1248  class global_function_x_bspline_
1249  : public global_function_simple, public context_dependencies {
1250  scalar_type xmin, xmax, xscale;
1251  std::function<scalar_type(scalar_type)> fx, fx_der, fx_der2;
1252  public:
1253  void update_from_context() const {}
1254 
1255  virtual scalar_type val(const base_node &pt) const
1256  {
1257  return fx(xscale*(pt[0]-xmin));
1258  }
1259  virtual void grad(const base_node &pt, base_small_vector &g) const
1260  {
1261  scalar_type dx = xscale*(pt[0]-xmin);
1262  g.resize(dim_);
1263  g[0] = xscale * fx_der(dx);
1264  }
1265  virtual void hess(const base_node &pt, base_matrix &h) const
1266  {
1267  scalar_type dx = xscale*(pt[0]-xmin);
1268  h.resize(dim_, dim_);
1269  gmm::clear(h);
1270  h(0,0) = xscale*xscale * fx_der2(dx);
1271  }
1272 
1273  virtual bool is_in_support(const base_node &pt) const {
1274  scalar_type dx = pt[0]-(xmin+xmax)/2;
1275  return (gmm::abs(dx)+1e-9 < gmm::abs(xmax-xmin)/2);
1276  }
1277 
1278  virtual void bounding_box(base_node &bmin, base_node &bmax) const {
1279  GMM_ASSERT1(bmin.size() == dim_ && bmax.size() == dim_,
1280  "Wrong dimensions");
1281  bmin[0] = std::min(xmin,xmax);
1282  bmax[0] = std::max(xmin,xmax);
1283  }
1284 
1285  global_function_x_bspline_(const scalar_type &xmin_, const scalar_type &xmax_,
1286  const size_type &order, const size_type &xtype)
1287  : global_function_simple(1), xmin(xmin_), xmax(xmax_),
1288  xscale(scalar_type(xtype)/(xmax-xmin))
1289  {
1290  GMM_ASSERT1(order >= 3 && order <= 5, "Only orders 3 to 5 are supported");
1291  GMM_ASSERT1(xtype >= 1 && xtype <= order, "Wrong input");
1292  if (order == 3) {
1293  if (xtype == 1) {
1294  fx = bsp3_01; fx_der = bsp3_01_der; fx_der2 = bsp3_01_der2;
1295  } else if (xtype == 2) {
1296  fx = bsp3_02; fx_der = bsp3_02_der; fx_der2 = bsp3_02_der2;
1297  } else if (xtype == 3) {
1298  fx = bsp3_03; fx_der = bsp3_03_der; fx_der2 = bsp3_03_der2;
1299  }
1300  } else if (order == 4) {
1301  if (xtype == 1) {
1302  fx = bsp4_01; fx_der = bsp4_01_der; fx_der2 = bsp4_01_der2;
1303  } else if (xtype == 2) {
1304  fx = bsp4_02; fx_der = bsp4_02_der; fx_der2 = bsp4_02_der2;
1305  } else if (xtype == 3) {
1306  fx = bsp4_03; fx_der = bsp4_03_der; fx_der2 = bsp4_03_der2;
1307  } else if (xtype == 4) {
1308  fx = bsp4_04; fx_der = bsp4_04_der; fx_der2 = bsp4_04_der2;
1309  }
1310  } else if (order == 5) {
1311  if (xtype == 1) {
1312  fx = bsp5_01; fx_der = bsp5_01_der; fx_der2 = bsp5_01_der2;
1313  } else if (xtype == 2) {
1314  fx = bsp5_02; fx_der = bsp5_02_der; fx_der2 = bsp5_02_der2;
1315  } else if (xtype == 3) {
1316  fx = bsp5_03; fx_der = bsp5_03_der; fx_der2 = bsp5_03_der2;
1317  } else if (xtype == 4) {
1318  fx = bsp5_04; fx_der = bsp5_04_der; fx_der2 = bsp5_04_der2;
1319  } else if (xtype == 5) {
1320  fx = bsp5_05; fx_der = bsp5_05_der; fx_der2 = bsp5_05_der2;
1321  }
1322  }
1323  }
1324 
1325  virtual ~global_function_x_bspline_()
1326  { DAL_STORED_OBJECT_DEBUG_DESTROYED(this, "Global function x bspline"); }
1327  };
1328 
1329 
1330 
1331  class global_function_xy_bspline_
1332  : public global_function_simple, public context_dependencies {
1333  scalar_type xmin, ymin, xmax, ymax, xscale, yscale;
1334  std::function<scalar_type(scalar_type)>
1335  fx, fy, fx_der, fy_der, fx_der2, fy_der2;
1336  public:
1337  void update_from_context() const {}
1338 
1339  virtual scalar_type val(const base_node &pt) const
1340  {
1341  return fx(xscale*(pt[0]-xmin))
1342  * fy(yscale*(pt[1]-ymin));
1343  }
1344  virtual void grad(const base_node &pt, base_small_vector &g) const
1345  {
1346  scalar_type dx = xscale*(pt[0]-xmin),
1347  dy = yscale*(pt[1]-ymin);
1348  g.resize(dim_);
1349  g[0] = xscale * fx_der(dx) * fy(dy);
1350  g[1] = fx(dx) * yscale * fy_der(dy);
1351  }
1352  virtual void hess(const base_node &pt, base_matrix &h) const
1353  {
1354  scalar_type dx = xscale*(pt[0]-xmin),
1355  dy = yscale*(pt[1]-ymin);
1356  h.resize(dim_, dim_);
1357  gmm::clear(h);
1358  h(0,0) = xscale*xscale * fx_der2(dx) * fy(dy);
1359  h(0,1) = h(1,0) = xscale * fx_der(dx) * yscale * fy_der(dy);
1360  h(1,1) = fx(dx) * yscale*yscale * fy_der2(dy);
1361  }
1362 
1363  virtual bool is_in_support(const base_node &pt) const {
1364  scalar_type dx = pt[0]-(xmin+xmax)/2,
1365  dy = pt[1]-(ymin+ymax)/2;
1366  return (gmm::abs(dx)+1e-9 < gmm::abs(xmax-xmin)/2) &&
1367  (gmm::abs(dy)+1e-9 < gmm::abs(ymax-ymin)/2);
1368  }
1369 
1370  virtual void bounding_box(base_node &bmin, base_node &bmax) const {
1371  GMM_ASSERT1(bmin.size() == dim_ && bmax.size() == dim_,
1372  "Wrong dimensions");
1373  bmin[0] = std::min(xmin,xmax);
1374  bmin[1] = std::min(ymin,ymax);
1375  bmax[0] = std::max(xmin,xmax);
1376  bmax[1] = std::max(ymin,ymax);
1377  }
1378 
1379  global_function_xy_bspline_(const scalar_type &xmin_, const scalar_type &xmax_,
1380  const scalar_type &ymin_, const scalar_type &ymax_,
1381  const size_type &order,
1382  const size_type &xtype, const size_type &ytype)
1383  : global_function_simple(2), xmin(xmin_), ymin(ymin_), xmax(xmax_), ymax(ymax_),
1384  xscale(scalar_type(xtype)/(xmax-xmin)), yscale(scalar_type(ytype)/(ymax-ymin))
1385  {
1386  GMM_ASSERT1(order >= 3 && order <= 5, "Wrong input");
1387  GMM_ASSERT1(xtype >= 1 && xtype <= order &&
1388  ytype >= 1 && ytype <= order, "Wrong input");
1389  if (order == 3) {
1390  if (xtype == 1) {
1391  fx = bsp3_01; fx_der = bsp3_01_der; fx_der2 = bsp3_01_der2;
1392  } else if (xtype == 2) {
1393  fx = bsp3_02; fx_der = bsp3_02_der; fx_der2 = bsp3_02_der2;
1394  } else if (xtype == 3) {
1395  fx = bsp3_03; fx_der = bsp3_03_der; fx_der2 = bsp3_03_der2;
1396  }
1397 
1398  if (ytype == 1) {
1399  fy = bsp3_01; fy_der = bsp3_01_der; fy_der2 = bsp3_01_der2;
1400  } else if (ytype == 2) {
1401  fy = bsp3_02; fy_der = bsp3_02_der; fy_der2 = bsp3_02_der2;
1402  } else if (ytype == 3) {
1403  fy = bsp3_03; fy_der = bsp3_03_der; fy_der2 = bsp3_03_der2;
1404  }
1405  } else if (order == 4) {
1406  if (xtype == 1) {
1407  fx = bsp4_01; fx_der = bsp4_01_der; fx_der2 = bsp4_01_der2;
1408  } else if (xtype == 2) {
1409  fx = bsp4_02; fx_der = bsp4_02_der; fx_der2 = bsp4_02_der2;
1410  } else if (xtype == 3) {
1411  fx = bsp4_03; fx_der = bsp4_03_der; fx_der2 = bsp4_03_der2;
1412  } else if (xtype == 4) {
1413  fx = bsp4_04; fx_der = bsp4_04_der; fx_der2 = bsp4_04_der2;
1414  }
1415 
1416  if (ytype == 1) {
1417  fy = bsp4_01; fy_der = bsp4_01_der; fy_der2 = bsp4_01_der2;
1418  } else if (ytype == 2) {
1419  fy = bsp4_02; fy_der = bsp4_02_der; fy_der2 = bsp4_02_der2;
1420  } else if (ytype == 3) {
1421  fy = bsp4_03; fy_der = bsp4_03_der; fy_der2 = bsp4_03_der2;
1422  } else if (ytype == 4) {
1423  fy = bsp4_04; fy_der = bsp4_04_der; fy_der2 = bsp4_04_der2;
1424  }
1425 
1426  } else if (order == 5) {
1427  if (xtype == 1) {
1428  fx = bsp5_01; fx_der = bsp5_01_der; fx_der2 = bsp5_01_der2;
1429  } else if (xtype == 2) {
1430  fx = bsp5_02; fx_der = bsp5_02_der; fx_der2 = bsp5_02_der2;
1431  } else if (xtype == 3) {
1432  fx = bsp5_03; fx_der = bsp5_03_der; fx_der2 = bsp5_03_der2;
1433  } else if (xtype == 4) {
1434  fx = bsp5_04; fx_der = bsp5_04_der; fx_der2 = bsp5_04_der2;
1435  } else if (xtype == 5) {
1436  fx = bsp5_05; fx_der = bsp5_05_der; fx_der2 = bsp5_05_der2;
1437  }
1438 
1439  if (ytype == 1) {
1440  fy = bsp5_01; fy_der = bsp5_01_der; fy_der2 = bsp5_01_der2;
1441  } else if (ytype == 2) {
1442  fy = bsp5_02; fy_der = bsp5_02_der; fy_der2 = bsp5_02_der2;
1443  } else if (ytype == 3) {
1444  fy = bsp5_03; fy_der = bsp5_03_der; fy_der2 = bsp5_03_der2;
1445  } else if (ytype == 4) {
1446  fy = bsp5_04; fy_der = bsp5_04_der; fy_der2 = bsp5_04_der2;
1447  } else if (ytype == 5) {
1448  fy = bsp5_05; fy_der = bsp5_05_der; fy_der2 = bsp5_05_der2;
1449  }
1450  }
1451  }
1452 
1453  virtual ~global_function_xy_bspline_()
1454  { DAL_STORED_OBJECT_DEBUG_DESTROYED(this, "Global function xy bspline"); }
1455  };
1456 
1457 
1458  class global_function_xyz_bspline_
1459  : public global_function_simple, public context_dependencies {
1460  scalar_type xmin, ymin, zmin, xmax, ymax, zmax, xscale, yscale, zscale;
1461  std::function<scalar_type(scalar_type)>
1462  fx, fy, fz, fx_der, fy_der, fz_der, fx_der2, fy_der2, fz_der2;
1463  public:
1464  void update_from_context() const {}
1465 
1466  virtual scalar_type val(const base_node &pt) const
1467  {
1468  return fx(xscale*(pt[0]-xmin))
1469  * fy(yscale*(pt[1]-ymin))
1470  * fz(zscale*(pt[2]-zmin));
1471  }
1472  virtual void grad(const base_node &pt, base_small_vector &g) const
1473  {
1474  scalar_type dx = xscale*(pt[0]-xmin),
1475  dy = yscale*(pt[1]-ymin),
1476  dz = zscale*(pt[2]-zmin);
1477  g.resize(dim_);
1478  g[0] = xscale * fx_der(dx) * fy(dy) * fz(dz);
1479  g[1] = fx(dx) * yscale * fy_der(dy) * fz(dz);
1480  g[2] = fx(dx) * fy(dy) * zscale * fz_der(dz);
1481  }
1482  virtual void hess(const base_node &pt, base_matrix &h) const
1483  {
1484  scalar_type dx = xscale*(pt[0]-xmin),
1485  dy = yscale*(pt[1]-ymin),
1486  dz = zscale*(pt[2]-zmin);
1487  h.resize(dim_, dim_);
1488  gmm::clear(h);
1489  h(0,0) = xscale*xscale * fx_der2(dx) * fy(dy) * fz(dz);
1490  h(1,1) = fx(dx) * yscale*yscale * fy_der2(dy) * fz(dz);
1491  h(2,2) = fx(dx) * fy(dy) * zscale*zscale * fz_der2(dz);
1492  h(0,1) = h(1,0) = xscale * fx_der(dx) * yscale * fy_der(dy) * fz(dz);
1493  h(0,2) = h(2,0) = xscale * fx_der(dx) * fy(dy) * zscale * fz_der(dz);
1494  h(1,2) = h(2,1) = fx(dx) * yscale * fy_der(dy) * zscale * fz_der(dz);
1495  }
1496 
1497  virtual bool is_in_support(const base_node &pt) const {
1498  scalar_type dx = pt[0]-(xmin+xmax)/2,
1499  dy = pt[1]-(ymin+ymax)/2,
1500  dz = pt[2]-(zmin+zmax)/2;
1501  return (gmm::abs(dx)+1e-9 < gmm::abs(xmax-xmin)/2) &&
1502  (gmm::abs(dy)+1e-9 < gmm::abs(ymax-ymin)/2) &&
1503  (gmm::abs(dz)+1e-9 < gmm::abs(zmax-zmin)/2);
1504  }
1505 
1506  virtual void bounding_box(base_node &bmin, base_node &bmax) const {
1507  GMM_ASSERT1(bmin.size() == dim_ && bmax.size() == dim_,
1508  "Wrong dimensions");
1509  bmin[0] = std::min(xmin,xmax);
1510  bmin[1] = std::min(ymin,ymax);
1511  bmin[2] = std::min(zmin,zmax);
1512  bmax[0] = std::max(xmin,xmax);
1513  bmax[1] = std::max(ymin,ymax);
1514  bmax[2] = std::max(zmin,zmax);
1515  }
1516 
1517  global_function_xyz_bspline_(const scalar_type &xmin_, const scalar_type &xmax_,
1518  const scalar_type &ymin_, const scalar_type &ymax_,
1519  const scalar_type &zmin_, const scalar_type &zmax_,
1520  const size_type &order,
1521  const size_type &xtype,
1522  const size_type &ytype,
1523  const size_type &ztype)
1524  : global_function_simple(3), xmin(xmin_), ymin(ymin_), zmin(zmin_),
1525  xmax(xmax_), ymax(ymax_), zmax(zmax_),
1526  xscale(scalar_type(xtype)/(xmax-xmin)),
1527  yscale(scalar_type(ytype)/(ymax-ymin)),
1528  zscale(scalar_type(ztype)/(zmax-zmin))
1529  {
1530  GMM_ASSERT1(order >= 3 && order <= 5, "Wrong input");
1531  GMM_ASSERT1(xtype >= 1 && xtype <= order &&
1532  ytype >= 1 && ytype <= order &&
1533  ztype >= 1 && ztype <= order, "Wrong input");
1534  if (order == 3) {
1535 
1536  if (xtype == 1) {
1537  fx = bsp3_01; fx_der = bsp3_01_der; fx_der2 = bsp3_01_der2;
1538  } else if (xtype == 2) {
1539  fx = bsp3_02; fx_der = bsp3_02_der; fx_der2 = bsp3_02_der2;
1540  } else if (xtype == 3) {
1541  fx = bsp3_03; fx_der = bsp3_03_der; fx_der2 = bsp3_03_der2;
1542  }
1543 
1544  if (ytype == 1) {
1545  fy = bsp3_01; fy_der = bsp3_01_der; fy_der2 = bsp3_01_der2;
1546  } else if (ytype == 2) {
1547  fy = bsp3_02; fy_der = bsp3_02_der; fy_der2 = bsp3_02_der2;
1548  } else if (ytype == 3) {
1549  fy = bsp3_03; fy_der = bsp3_03_der; fy_der2 = bsp3_03_der2;
1550  }
1551 
1552  if (ztype == 1) {
1553  fz = bsp3_01; fz_der = bsp3_01_der; fz_der2 = bsp3_01_der2;
1554  } else if (ztype == 2) {
1555  fz = bsp3_02; fz_der = bsp3_02_der; fz_der2 = bsp3_02_der2;
1556  } else if (ztype == 3) {
1557  fz = bsp3_03; fz_der = bsp3_03_der; fz_der2 = bsp3_03_der2;
1558  }
1559 
1560  } else if (order == 4) {
1561 
1562  if (xtype == 1) {
1563  fx = bsp4_01; fx_der = bsp4_01_der; fx_der2 = bsp4_01_der2;
1564  } else if (xtype == 2) {
1565  fx = bsp4_02; fx_der = bsp4_02_der; fx_der2 = bsp4_02_der2;
1566  } else if (xtype == 3) {
1567  fx = bsp4_03; fx_der = bsp4_03_der; fx_der2 = bsp4_03_der2;
1568  } else if (xtype == 4) {
1569  fx = bsp4_04; fx_der = bsp4_04_der; fx_der2 = bsp4_04_der2;
1570  }
1571 
1572  if (ytype == 1) {
1573  fy = bsp4_01; fy_der = bsp4_01_der; fy_der2 = bsp4_01_der2;
1574  } else if (ytype == 2) {
1575  fy = bsp4_02; fy_der = bsp4_02_der; fy_der2 = bsp4_02_der2;
1576  } else if (ytype == 3) {
1577  fy = bsp4_03; fy_der = bsp4_03_der; fy_der2 = bsp4_03_der2;
1578  } else if (ytype == 4) {
1579  fy = bsp4_04; fy_der = bsp4_04_der; fy_der2 = bsp4_04_der2;
1580  }
1581 
1582  if (ztype == 1) {
1583  fz = bsp4_01; fz_der = bsp4_01_der; fz_der2 = bsp4_01_der2;
1584  } else if (ztype == 2) {
1585  fz = bsp4_02; fz_der = bsp4_02_der; fz_der2 = bsp4_02_der2;
1586  } else if (ztype == 3) {
1587  fz = bsp4_03; fz_der = bsp4_03_der; fz_der2 = bsp4_03_der2;
1588  } else if (ztype == 4) {
1589  fz = bsp4_04; fz_der = bsp4_04_der; fz_der2 = bsp4_04_der2;
1590  }
1591 
1592  } else if (order == 5) {
1593 
1594  if (xtype == 1) {
1595  fx = bsp5_01; fx_der = bsp5_01_der; fx_der2 = bsp5_01_der2;
1596  } else if (xtype == 2) {
1597  fx = bsp5_02; fx_der = bsp5_02_der; fx_der2 = bsp5_02_der2;
1598  } else if (xtype == 3) {
1599  fx = bsp5_03; fx_der = bsp5_03_der; fx_der2 = bsp5_03_der2;
1600  } else if (xtype == 4) {
1601  fx = bsp5_04; fx_der = bsp5_04_der; fx_der2 = bsp5_04_der2;
1602  } else if (xtype == 5) {
1603  fx = bsp5_05; fx_der = bsp5_05_der; fx_der2 = bsp5_05_der2;
1604  }
1605 
1606  if (ytype == 1) {
1607  fy = bsp5_01; fy_der = bsp5_01_der; fy_der2 = bsp5_01_der2;
1608  } else if (ytype == 2) {
1609  fy = bsp5_02; fy_der = bsp5_02_der; fy_der2 = bsp5_02_der2;
1610  } else if (ytype == 3) {
1611  fy = bsp5_03; fy_der = bsp5_03_der; fy_der2 = bsp5_03_der2;
1612  } else if (ytype == 4) {
1613  fy = bsp5_04; fy_der = bsp5_04_der; fy_der2 = bsp5_04_der2;
1614  } else if (ytype == 5) {
1615  fy = bsp5_05; fy_der = bsp5_05_der; fy_der2 = bsp5_05_der2;
1616  }
1617 
1618  if (ztype == 1) {
1619  fz = bsp5_01; fz_der = bsp5_01_der; fz_der2 = bsp5_01_der2;
1620  } else if (ztype == 2) {
1621  fz = bsp5_02; fz_der = bsp5_02_der; fz_der2 = bsp5_02_der2;
1622  } else if (ztype == 3) {
1623  fz = bsp5_03; fz_der = bsp5_03_der; fz_der2 = bsp5_03_der2;
1624  } else if (ztype == 4) {
1625  fz = bsp5_04; fz_der = bsp5_04_der; fz_der2 = bsp5_04_der2;
1626  } else if (ztype == 5) {
1627  fz = bsp5_05; fz_der = bsp5_05_der; fz_der2 = bsp5_05_der2;
1628  }
1629 
1630  }
1631  }
1632 
1633  virtual ~global_function_xyz_bspline_()
1634  { DAL_STORED_OBJECT_DEBUG_DESTROYED(this, "Global function xyz bspline"); }
1635  };
1636 
1637 
1638  pglobal_function
1639  global_function_bspline(const scalar_type xmin, const scalar_type xmax,
1640  const size_type order, const size_type xtype) {
1641  return std::make_shared<global_function_x_bspline_>
1642  (xmin, xmax, order, xtype);
1643  }
1644 
1645  pglobal_function
1646  global_function_bspline(const scalar_type xmin, const scalar_type xmax,
1647  const scalar_type ymin, const scalar_type ymax,
1648  const size_type order,
1649  const size_type xtype, const size_type ytype) {
1650  return std::make_shared<global_function_xy_bspline_>
1651  (xmin, xmax, ymin, ymax, order, xtype, ytype);
1652  }
1653 
1654  pglobal_function
1655  global_function_bspline(const scalar_type xmin, const scalar_type xmax,
1656  const scalar_type ymin, const scalar_type ymax,
1657  const scalar_type zmin, const scalar_type zmax,
1658  const size_type order,
1659  const size_type xtype,
1660  const size_type ytype,
1661  const size_type ztype) {
1662  return std::make_shared<global_function_xyz_bspline_>
1663  (xmin, xmax, ymin, ymax, zmin, zmax, order,
1664  xtype, ytype, ztype);
1665  }
1666 
1667 
1668 
1669  // interpolator on mesh_fem
1670 
1671  interpolator_on_mesh_fem::interpolator_on_mesh_fem
1672  (const mesh_fem &mf_, const std::vector<scalar_type> &U_)
1673  : mf(mf_), U(U_) {
1674 
1675  if (mf.is_reduced()) {
1676  gmm::resize(U, mf.nb_basic_dof());
1677  gmm::mult(mf.extension_matrix(), U_, U);
1678  }
1679  base_node bmin, bmax;
1680  scalar_type EPS=1E-13;
1681  cv_stored = size_type(-1);
1682  boxtree.clear();
1683  for (dal::bv_visitor cv(mf.convex_index()); !cv.finished(); ++cv) {
1684  bounding_box(bmin, bmax, mf.linked_mesh().points_of_convex(cv),
1685  mf.linked_mesh().trans_of_convex(cv));
1686  for (auto&& val : bmin) val -= EPS;
1687  for (auto&& val : bmax) val += EPS;
1688  boxtree.add_box(bmin, bmax, cv);
1689  }
1690  boxtree.build_tree();
1691  }
1692 
1693  bool interpolator_on_mesh_fem::find_a_point(const base_node &pt,
1694  base_node &ptr,
1695  size_type &cv) const {
1696  bool gt_invertible;
1697  if (cv_stored != size_type(-1) && gic.invert(pt, ptr, gt_invertible)) {
1698  cv = cv_stored;
1699  if (gt_invertible)
1700  return true;
1701  }
1702 
1703  boxtree.find_boxes_at_point(pt, boxlst);
1704  for (const auto &box : boxlst) {
1705  gic = bgeot::geotrans_inv_convex
1706  (mf.linked_mesh().convex(box->id),
1707  mf.linked_mesh().trans_of_convex(box->id));
1708  cv_stored = box->id;
1709  if (gic.invert(pt, ptr, gt_invertible)) {
1710  cv = box->id;
1711  return true;
1712  }
1713  }
1714  return false;
1715  }
1716 
1717  bool interpolator_on_mesh_fem::eval(const base_node &pt, base_vector &val,
1718  base_matrix &grad) const {
1719  base_node ptref;
1720  size_type cv;
1721  base_vector coeff;
1722  dim_type q = mf.get_qdim(), N = mf.linked_mesh().dim();
1723  if (find_a_point(pt, ptref, cv)) {
1724  pfem pf = mf.fem_of_element(cv);
1725  bgeot::pgeometric_trans pgt =
1726  mf.linked_mesh().trans_of_convex(cv);
1727  base_matrix G;
1728  vectors_to_base_matrix(G, mf.linked_mesh().points_of_convex(cv));
1729  fem_interpolation_context fic(pgt, pf, ptref, G, cv, short_type(-1));
1730  slice_vector_on_basic_dof_of_element(mf, U, cv, coeff);
1731  // coeff.resize(mf.nb_basic_dof_of_element(cv));
1732  // gmm::copy(gmm::sub_vector
1733  // (U,gmm::sub_index(mf.ind_basic_dof_of_element(cv))), coeff);
1734  val.resize(q);
1735  pf->interpolation(fic, coeff, val, q);
1736  grad.resize(q, N);
1737  pf->interpolation_grad(fic, coeff, grad, q);
1738  return true;
1739  } else
1740  return false;
1741  }
1742 }
1743 
1744 /* 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
getfem_global_function.h
Definition of global functions to be used as base or enrichment functions in fem.
getfem_level_set.h
Define level-sets.
gmm::reshape
void reshape(M &v, size_type m, size_type n)
*‍/
Definition: gmm_blas.h:251
gmm::copy
void copy(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:978
gmm::fill
void fill(L &l, typename gmm::linalg_traits< L >::value_type x)
*‍/
Definition: gmm_blas.h:104
gmm::clear
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:59
gmm::resize
void resize(V &v, size_type n)
*‍/
Definition: gmm_blas.h:210
gmm::mult
void mult(const L1 &l1, const L2 &l2, L3 &l3)
*‍/
Definition: gmm_blas.h:1664
gmm::add
void add(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:1277
getfem::pfem
std::shared_ptr< const getfem::virtual_fem > pfem
type of pointer on a fem description
Definition: getfem_fem.h:244
bgeot::short_type
gmm::uint16_type short_type
used as the common short type integer in the library
Definition: bgeot_config.h:73
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
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
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:47
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:663