GetFEM  5.5
getfem_global_function.cc
1 /*===========================================================================
2 
3  Copyright (C) 2004-2026 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, see https://www.gnu.org/licenses/.
19 
20 ===========================================================================*/
21 
22 #include <getfem/getfem_global_function.h>
23 #include <getfem/getfem_level_set.h>
24 
25 namespace getfem {
26 
27 
28  // Partial implementation of abstract class global_function_simple
29 
30  scalar_type global_function_simple::val
31  (const fem_interpolation_context &c) const {
32  base_node pt = c.xreal();
33  GMM_ASSERT1(pt.size() == dim_, "Point of wrong size (" << pt.size() << ") "
34  << "passed to a global function of dim = "<< dim_ <<".");
35  return this->val(pt);
36  }
37 
38  void global_function_simple::grad
39  (const fem_interpolation_context &c, base_small_vector &g) const {
40  base_node pt = c.xreal();
41  GMM_ASSERT1(pt.size() == dim_, "Point of wrong size (" << pt.size() << ") "
42  << "passed to a global function of dim = "<< dim_ <<".");
43  this->grad(pt, g);
44  }
45 
46  void global_function_simple::hess
47  (const fem_interpolation_context &c, base_matrix &h) const {
48  base_node pt = c.xreal();
49  GMM_ASSERT1(pt.size() == dim_, "Point of wrong size (" << pt.size() << ") "
50  << "passed to a global function of dim = "<< dim_ <<".");
51  this->hess(pt, h);
52  }
53 
54  // Implementation of global_function_parser
55 
56  scalar_type global_function_parser::val(const base_node &pt) const {
57  const bgeot::base_tensor &t = tensor_val(pt);
58  GMM_ASSERT1(t.size() == 1, "Wrong size of expression result "
59  << f_val.expression());
60  return t[0];
61  }
62 
63  const base_tensor &global_function_parser::tensor_val(const base_node &pt) const {
64  gmm::copy(pt, pt_);
65  return f_val.eval();
66  }
67 
68  void global_function_parser::grad(const base_node &pt, base_small_vector &g) const {
69  g.resize(dim_);
70  gmm::copy(pt, pt_);
71  const bgeot::base_tensor &t = f_grad.eval();
72  GMM_ASSERT1(t.size() == dim_, "Wrong size of expression result "
73  << f_grad.expression());
74  gmm::copy(t.as_vector(), g);
75 
76  }
77 
78  void global_function_parser::hess(const base_node &pt, base_matrix &h) const {
79  h.resize(dim_, dim_);
80  gmm::copy(pt, pt_);
81  const bgeot::base_tensor &t = f_hess.eval();
82  GMM_ASSERT1(t.size() == size_type(dim_*dim_),
83  "Wrong size of expression result " << f_hess.expression());
84  gmm::copy(t.as_vector(), h.as_vector());
85  }
86 
87  global_function_parser::global_function_parser(dim_type dim__,
88  const std::string &sval,
89  const std::string &sgrad,
90  const std::string &shess)
91  : global_function_simple(dim__),
92  f_val(gw, sval), f_grad(gw, sgrad), f_hess(gw, shess) {
93 
94  size_type N(dim_);
95  pt_.resize(N);
96  gmm::fill(pt_, scalar_type(0));
97  gw.add_fixed_size_variable("X", gmm::sub_interval(0, N), pt_);
98  if (N >= 1) gw.add_macro("x", "X(1)");
99  if (N >= 2) gw.add_macro("y", "X(2)");
100  if (N >= 3) gw.add_macro("z", "X(3)");
101  if (N >= 4) gw.add_macro("w", "X(4)");
102 
103  f_val.compile();
104  f_grad.compile();
105  f_hess.compile();
106  }
107 
108  // Implementation of global_function_sum
109 
110  scalar_type global_function_sum::val
111  (const fem_interpolation_context &c) const {
112  scalar_type res(0);
113  for (const auto &f : functions)
114  res += f->val(c);
115  return res;
116  }
117 
118  void global_function_sum::grad
119  (const fem_interpolation_context &c, base_small_vector &g) const {
120  g.resize(dim_);
121  gmm::clear(g);
122  base_small_vector gg(dim_);
123  for (const auto &f : functions) {
124  f->grad(c, gg);
125  gmm::add(gg, g);
126  }
127  }
128 
129  void global_function_sum::hess
130  (const fem_interpolation_context &c, base_matrix &h) const {
131  h.resize(dim_, dim_);
132  gmm::clear(h);
133  base_matrix hh(dim_, dim_);
134  for (const auto &f : functions) {
135  f->hess(c, hh);
136  gmm::add(hh.as_vector(), h.as_vector());
137  }
138  }
139 
140  bool global_function_sum::is_in_support(const base_node &p) const {
141  for (const auto &f : functions)
142  if (f->is_in_support(p)) return true;
143  return false;
144  }
145 
146  void global_function_sum::bounding_box
147  (base_node &bmin_, base_node &bmax_) const {
148  if (functions.size() > 0)
149  functions[0]->bounding_box(bmin_, bmax_);
150  base_node bmin0(dim()), bmax0(dim());
151  for (const auto &f : functions) {
152  f->bounding_box(bmin0, bmax0);
153  for (size_type i=0; i < dim(); ++i) {
154  if (bmin0[i] < bmin_[i]) bmin_[i] = bmin0[i];
155  if (bmax0[i] > bmax_[i]) bmax_[i] = bmax0[i];
156  }
157  }
158  }
159 
160  global_function_sum::global_function_sum(const std::vector<pglobal_function> &funcs)
161  : global_function((funcs.size() > 0) ? funcs[0]->dim() : 0), functions(funcs) {
162  for (const auto &f : functions)
163  GMM_ASSERT1(f->dim() == dim(), "Incompatible dimensions among the provided"
164  " global functions");
165  }
166 
167  global_function_sum::global_function_sum(pglobal_function f1, pglobal_function f2)
168  : global_function(f1->dim()), functions(2) {
169  functions[0] = f1;
170  functions[1] = f2;
171  GMM_ASSERT1(f1->dim() == dim() && f2->dim() == dim(),
172  "Incompatible dimensions between the provided global functions");
173  }
174 
175  global_function_sum::global_function_sum(pglobal_function f1, pglobal_function f2,
176  pglobal_function f3)
177  : global_function(f1->dim()), functions(3) {
178  functions[0] = f1;
179  functions[1] = f2;
180  functions[2] = f3;
181  GMM_ASSERT1(f1->dim() == dim() && f2->dim() == dim() && f3->dim() == dim(),
182  "Incompatible dimensions between the provided global functions");
183  }
184 
185  global_function_sum::global_function_sum(pglobal_function f1, pglobal_function f2,
186  pglobal_function f3, pglobal_function f4)
187  : global_function(f1->dim()), functions(4) {
188  functions[0] = f1;
189  functions[1] = f2;
190  functions[2] = f3;
191  functions[3] = f4;
192  GMM_ASSERT1(f1->dim() == dim() && f2->dim() == dim() && f3->dim() == dim(),
193  "Incompatible dimensions between the provided global functions");
194  }
195 
196 
197  // Implementation of global_function_product
198 
199  scalar_type global_function_product::val
200  (const fem_interpolation_context &c) const {
201  return f1->val(c) * f2->val(c);
202  }
203 
204  void global_function_product::grad
205  (const fem_interpolation_context &c, base_small_vector &g) const {
206  g.resize(dim_);
207  base_small_vector gg(dim_);
208  f1->grad(c, gg);
209  gmm::copy(gmm::scaled(gg, f2->val(c)), g);
210  f2->grad(c, gg);
211  gmm::add(gmm::scaled(gg, f1->val(c)), g);
212  }
213 
214  void global_function_product::hess
215  (const fem_interpolation_context &c, base_matrix &h) const {
216  h.resize(dim_, dim_);
217  gmm::clear(h);
218  base_matrix hh(dim_, dim_);
219  f1->hess(c, hh);
220  gmm::copy(gmm::scaled(hh, f2->val(c)), h);
221  f2->hess(c, hh);
222  gmm::add(gmm::scaled(hh, f1->val(c)), h);
223  base_small_vector g1(dim_), g2(dim_);
224  f1->grad(c, g1);
225  f2->grad(c, g2);
226  gmm::rank_one_update(h, g1, g2);
227  gmm::rank_one_update(h, g2, g1);
228  }
229 
230  bool global_function_product::is_in_support(const base_node &p) const {
231  return f1->is_in_support(p) && f2->is_in_support(p);
232  }
233 
234  void global_function_product::bounding_box
235  (base_node &bmin_, base_node &bmax_) const {
236  base_node bmin0(dim()), bmax0(dim());
237  f1->bounding_box(bmin0, bmax0);
238  f2->bounding_box(bmin_, bmax_);
239  for (size_type i=0; i < dim(); ++i) {
240  if (bmin0[i] > bmin_[i]) bmin_[i] = bmin0[i];
241  if (bmax0[i] < bmax_[i]) bmax_[i] = bmax0[i];
242  if (bmin_[i] > bmax_[i])
243  GMM_WARNING1("Global function product with vanishing basis function");
244  }
245  }
246 
247  global_function_product::global_function_product(pglobal_function f1_, pglobal_function f2_)
248  : global_function(f1_->dim()), f1(f1_), f2(f2_) {
249  GMM_ASSERT1(f2->dim() == dim(), "Incompatible dimensions between the provided"
250  " global functions");
251  }
252 
253 
254  // Implementation of global_function_bounded
255 
256  bool global_function_bounded::is_in_support(const base_node &pt) const {
257  if (has_expr) {
258  gmm::copy(pt, pt_);
259  const bgeot::base_tensor &t = f_val.eval();
260  GMM_ASSERT1(t.size() == 1, "Wrong size of expression result "
261  << f_val.expression());
262  return (t[0] > scalar_type(0));
263  }
264  return true;
265  }
266 
267  global_function_bounded::global_function_bounded(pglobal_function f_,
268  const base_node &bmin_,
269  const base_node &bmax_,
270  const std::string &is_in_expr)
271  : global_function(f_->dim()), f(f_), bmin(bmin_), bmax(bmax_),
272  f_val(gw, is_in_expr) {
273 
274  has_expr = !is_in_expr.empty();
275  if (has_expr) {
276  size_type N(dim_);
277  pt_.resize(N);
278  gmm::fill(pt_, scalar_type(0));
279  gw.add_fixed_size_variable("X", gmm::sub_interval(0, N), pt_);
280  if (N >= 1) gw.add_macro("x", "X(1)");
281  if (N >= 2) gw.add_macro("y", "X(2)");
282  if (N >= 3) gw.add_macro("z", "X(3)");
283  if (N >= 4) gw.add_macro("w", "X(4)");
284  f_val.compile();
285  }
286  }
287 
288  // Implementation of some useful xy functions
289 
290  parser_xy_function::parser_xy_function(const std::string &sval,
291  const std::string &sgrad,
292  const std::string &shess)
293  : f_val(gw, sval), f_grad(gw, sgrad), f_hess(gw, shess),
294  ptx(1), pty(1), ptr(1), ptt(1) {
295 
296  gw.add_fixed_size_constant("x", ptx);
297  gw.add_fixed_size_constant("y", pty);
298  gw.add_fixed_size_constant("r", ptr);
299  gw.add_fixed_size_constant("theta", ptt);
300 
301  f_val.compile();
302  f_grad.compile();
303  f_hess.compile();
304  }
305 
306  scalar_type
307  parser_xy_function::val(scalar_type x, scalar_type y) const {
308  ptx[0] = double(x); // x
309  pty[0] = double(y); // y
310  ptr[0] = double(sqrt(fabs(x*x+y*y))); // r
311  ptt[0] = double(atan2(y,x)); // theta
312 
313  const bgeot::base_tensor &t = f_val.eval();
314  GMM_ASSERT1(t.size() == 1, "Wrong size of expression result "
315  << f_val.expression());
316  return t[0];
317  }
318 
319  base_small_vector
320  parser_xy_function::grad(scalar_type x, scalar_type y) const {
321  ptx[0] = double(x); // x
322  pty[0] = double(y); // y
323  ptr[0] = double(sqrt(fabs(x*x+y*y))); // r
324  ptt[0] = double(atan2(y,x)); // theta
325 
326  base_small_vector res(2);
327  const bgeot::base_tensor &t = f_grad.eval();
328  GMM_ASSERT1(t.size() == 2, "Wrong size of expression result "
329  << f_grad.expression());
330  gmm::copy(t.as_vector(), res);
331  return res;
332  }
333 
334  base_matrix
335  parser_xy_function::hess(scalar_type x, scalar_type y) const {
336  ptx[0] = double(x); // x
337  pty[0] = double(y); // y
338  ptr[0] = double(sqrt(fabs(x*x+y*y))); // r
339  ptt[0] = double(atan2(y,x)); // theta
340 
341  base_matrix res(2,2);
342  const bgeot::base_tensor &t = f_hess.eval();
343  GMM_ASSERT1(t.size() == 4, "Wrong size of expression result "
344  << f_hess.expression());
345  gmm::copy(t.as_vector(), res.as_vector());
346  return res;
347  }
348 
349  /* the basic singular functions for 2D cracks */
350  scalar_type
351  crack_singular_xy_function::val(scalar_type x, scalar_type y) const {
352  scalar_type sgny = (y < 0 ? -1.0 : 1.0);
353  scalar_type r = sqrt(x*x + y*y);
354  if (r < 1e-10) return 0;
355 
356  /* The absolute value is unfortunately necessary, otherwise, sqrt(-1e-16)
357  can be required ...
358  */
359  scalar_type sin2 = sqrt(gmm::abs(.5-x/(2*r))) * sgny;
360  scalar_type cos2 = sqrt(gmm::abs(.5+x/(2*r)));
361 
362  scalar_type res = 0;
363  switch(l){
364 
365  /* First order enrichement displacement field (linear elasticity) */
366 
367  case 0 : res = sqrt(r)*sin2; break;
368  case 1 : res = sqrt(r)*cos2; break;
369  case 2 : res = sin2*y/sqrt(r); break;
370  case 3 : res = cos2*y/sqrt(r); break;
371 
372  /* Second order enrichement of displacement field (linear elasticity) */
373 
374  case 4 : res = sqrt(r)*r*sin2; break;
375  case 5 : res = sqrt(r)*r*cos2; break;
376  case 6 : res = sin2*cos2*cos2*r*sqrt(r); break;
377  case 7 : res = cos2*cos2*cos2*r*sqrt(r); break;
378 
379  /* First order enrichement of pressure field (linear elasticity) mixed formulation */
380 
381  case 8 : res = -sin2/sqrt(r); break;
382  case 9 : res = cos2/sqrt(r); break;
383 
384  /* Second order enrichement of pressure field (linear elasticity) mixed formulation */
385 
386  case 10 : res = sin2*sqrt(r); break; // cos2*cos2
387  case 11 : res = cos2*sqrt(r); break;
388 
389  /* First order enrichement of displacement field (nonlinear elasticity)[Rodney Stephenson Journal of elasticity VOL.12 No. 1, January 1982] */
390 
391  case 12 : res = r*sin2*sin2; break;
392  case 13 : res = sqrt(r)*sin2; break;
393 
394 /* First order enrichement of pressure field (nonlinear elasticity) */
395 
396  case 14 : res = sin2/r; break;
397  case 15 : res = cos2/r; break;
398 
399  default: GMM_ASSERT2(false, "arg");
400  }
401  return res;
402  }
403 
404 
405  base_small_vector
406  crack_singular_xy_function::grad(scalar_type x, scalar_type y) const {
407  scalar_type sgny = (y < 0 ? -1.0 : 1.0);
408  scalar_type r = sqrt(x*x + y*y);
409 
410  if (r < 1e-10) {
411  GMM_WARNING0("Warning, point close to the singularity (r=" << r << ")");
412  }
413 
414  /* The absolute value is unfortunately necessary, otherwise, sqrt(-1e-16)
415  can be required ...
416  */
417  scalar_type sin2 = sqrt(gmm::abs(.5-x/(2*r))) * sgny;
418  scalar_type cos2 = sqrt(gmm::abs(.5+x/(2*r)));
419 
420  base_small_vector res(2);
421  switch(l){
422  /* First order enrichement displacement field (linear elasticity) */
423  case 0 :
424  res[0] = -sin2/(2*sqrt(r));
425  res[1] = cos2/(2*sqrt(r));
426  break;
427  case 1 :
428  res[0] = cos2/(2*sqrt(r));
429  res[1] = sin2/(2*sqrt(r));
430  break;
431  case 2 :
432  res[0] = cos2*((-5*cos2*cos2) + 1. + 4*(cos2*cos2*cos2*cos2))/sqrt(r);
433  res[1] = sin2*((-3*cos2*cos2) + 1. + 4*(cos2*cos2*cos2*cos2))/sqrt(r);
434  break;
435  case 3 :
436  res[0] = -cos2*cos2*sin2*((4*cos2*cos2) - 3.)/sqrt(r);
437  res[1] = cos2*((4*cos2*cos2*cos2*cos2) + 2. - (5*cos2*cos2))/sqrt(r);
438  break;
439  /* Second order enrichement of displacement field (linear elasticity) */
440 
441  case 4 :
442  res[0] = sin2 *((4*cos2*cos2)-3.)*sqrt(r)/2.;
443  res[1] = cos2*(5. - (4*cos2*cos2))*sqrt(r)/2. ;
444  break;
445  case 5 :
446  res[0] = cos2*((4*cos2*cos2)-1.)*sqrt(r)/2.;
447  res[1] = sin2*((4*cos2*cos2)+1.)*sqrt(r)/2. ;
448  break;
449 
450  case 6 :
451  res[0] = sin2*cos2*cos2*sqrt(r)/2.;
452  res[1] = cos2*(2. - (cos2*cos2))*sqrt(r)/2.;
453  break;
454  case 7 :
455  res[0] = 3*cos2*cos2*cos2*sqrt(r)/2.;
456  res[1] = 3*sin2*cos2*cos2*sqrt(r)/2.;
457  break;
458 
459  /* First order enrichement of pressure field (linear elasticity) mixed formulation */
460 
461  case 8 :
462  res[0] =sin2*((4*cos2*cos2)-1.)/(2*sqrt(r)*r);
463  res[1] =-cos2*((4*cos2*cos2)-3.)/(2*sqrt(r)*r);
464  break;
465  case 9 :
466  res[0] =-cos2*((2*cos2*cos2) - 3.)/(2*sqrt(r)*r);
467  res[1] =-sin2*((4*cos2*cos2)-1.)/(2*sqrt(r)*r);
468  break;
469 
470  /* Second order enrichement of pressure field (linear elasticity) mixed formulation */
471  case 10 :
472  res[0] = -sin2/(2*sqrt(r));
473  res[1] = cos2/(2*sqrt(r));
474  break;
475  case 11 :
476  res[0] = cos2/(2*sqrt(r));
477  res[1] = sin2/(2*sqrt(r));
478  break;
479 
480  /* First order enrichement of displacement field (nonlinear elasticity)[Rodney Stephenson Journal of elasticity VOL.12 No. 1, January 1982] */
481 
482  case 12 :
483  res[0] = sin2*sin2;
484  res[1] = 0.5*cos2*sin2;
485  break;
486  case 13 :
487  res[0] = -sin2/(2*sqrt(r));
488  res[1] = cos2/(2*sqrt(r));
489  break;
490 
491 /* First order enrichement of pressure field (****Nonlinear elasticity*****) */
492 
493 
494  case 14 :
495  res[0] = -sin2/r;
496  res[1] = cos2/(2*r);
497  break;
498  case 15 :
499  res[0] = -cos2/r;
500  res[1] = -sin2/(2*r);
501  break;
502 
503 
504  default: GMM_ASSERT2(false, "oups");
505  }
506  return res;
507  }
508 
509  base_matrix
510  crack_singular_xy_function::hess(scalar_type x, scalar_type y) const {
511  scalar_type sgny = (y < 0 ? -1.0 : 1.0);
512  scalar_type r = sqrt(x*x + y*y);
513 
514  if (r < 1e-10) {
515  GMM_WARNING0("Warning, point close to the singularity (r=" << r << ")");
516  }
517 
518  /* The absolute value is unfortunately necessary, otherwise, sqrt(-1e-16)
519  can be required ...
520  */
521  scalar_type sin2 = sqrt(gmm::abs(.5-x/(2*r))) * sgny;
522  scalar_type cos2 = sqrt(gmm::abs(.5+x/(2*r)));
523 
524  base_matrix res(2,2);
525  switch(l){
526  case 0 :
527  res(0,0) = (-sin2*x*x + 2.0*cos2*x*y + sin2*y*y) / (4*pow(r, 3.5));
528  res(0,1) = (-cos2*x*x - 2.0*sin2*x*y + cos2*y*y) / (4*pow(r, 3.5));
529  res(1,0) = res(0, 1);
530  res(1,1) = (sin2*x*x - 2.0*cos2*x*y - sin2*y*y) / (4*pow(r, 3.5));
531  break;
532  case 1 :
533  res(0,0) = (-cos2*x*x - 2.0*sin2*x*y + cos2*y*y) / (4*pow(r, 3.5));
534  res(0,1) = (sin2*x*x - 2.0*cos2*x*y - sin2*y*y) / (4*pow(r, 3.5));
535  res(1,0) = res(0, 1);
536  res(1,1) = (cos2*x*x + 2.0*sin2*x*y - cos2*y*y) / (4*pow(r, 3.5));
537  break;
538  case 2 :
539  res(0,0) = 3.0*y*(sin2*x*x + 2.0*cos2*x*y - sin2*y*y) / (4*pow(r, 4.5));
540  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)
541  / (4*pow(r, 4.5));
542  res(1,0) = res(0, 1);
543  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)
544  / (4*pow(r, 4.5));
545  break;
546  case 3 :
547  res(0,0) = 3.0*y*(cos2*x*x - 2.0*sin2*x*y - cos2*y*y) / (4*pow(r, 4.5));
548  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)
549  / (4*pow(r, 4.5));
550  res(1,0) = res(0, 1);
551  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)
552  / (4*pow(r, 4.5));
553  break;
554  default: GMM_ASSERT2(false, "oups");
555  }
556  return res;
557  }
558 
559  scalar_type
560  cutoff_xy_function::val(scalar_type x, scalar_type y) const {
561  scalar_type res = 1;
562  switch (fun) {
563  case EXPONENTIAL_CUTOFF: {
564  res = (a4>0) ? exp(-a4 * gmm::sqr(x*x+y*y)) : 1;
565  break;
566  }
567  case POLYNOMIAL_CUTOFF: {
568  assert(r0 > r1);
569  scalar_type r = gmm::sqrt(x*x+y*y);
570 
571  if (r <= r1) res = 1;
572  else if (r >= r0) res = 0;
573  else {
574  // scalar_type c = 6./(r0*r0*r0 - r1*r1*r1 + 3*r1*r0*(r1-r0));
575  // scalar_type k = -(c/6.)*(-pow(r0,3.) + 3*r1*pow(r0,2.));
576  // res = (c/3.)*pow(r,3.) - (c*(r0 + r1)/2.)*pow(r,2.) +
577  // c*r0*r1*r + k;
578  scalar_type c = 1./pow(r0-r1,3.0);
579  res = c*(r*(r*(2.0*r-3.0*(r0+r1))+6.0*r1*r0) + r0*r0*(r0-3.0*r1));
580  }
581  break;
582  }
583  case POLYNOMIAL2_CUTOFF: {
584  assert(r0 > r1);
585  scalar_type r = gmm::sqrt(x*x+y*y);
586  if (r <= r1) res = scalar_type(1);
587  else if (r >= r0) res = scalar_type(0);
588  else {
589 // scalar_type c = 1./((-1./30.)*(pow(r1,5) - pow(r0,5))
590 // + (1./6.)*(pow(r1,4)*r0 - r1*pow(r0,4))
591 // - (1./3.)*(pow(r1,3)*pow(r0,2) -
592 // pow(r1,2)*pow(r0,3)));
593 // scalar_type k = 1. - c*((-1./30.)*pow(r1,5) +
594 // (1./6.)*pow(r1,4)*r0 -
595 // (1./3.)*pow(r1,3)*pow(r0,2));
596 // res = c*( (-1./5.)*pow(r,5) + (1./2.)* (r1+r0)*pow(r,4) -
597 // (1./3.)*(pow(r1,2)+pow(r0,2) + 4.*r0*r1)*pow(r,3) +
598 // r0*r1*(r0+r1)*pow(r,2) - pow(r0,2)*pow(r1,2)*r) + k;
599  res = (r*(r*(r*(r*(-6.0*r + 15.0*(r0+r1))
600  - 10.0*(r0*r0 + 4.0*r1*r0 + r1*r1))
601  + 30.0 * r0*r1*(r0+r1)) - 30.0*r1*r1*r0*r0)
602  + r0*r0*r0*(r0*r0-5.0*r1*r0+10*r1*r1)) / pow(r0-r1, 5.0);
603  }
604  break;
605  }
606  default : res = 1;
607  }
608  return res;
609  }
610 
611  base_small_vector
612  cutoff_xy_function::grad(scalar_type x, scalar_type y) const {
613  base_small_vector res(2);
614  switch (fun) {
615  case EXPONENTIAL_CUTOFF: {
616  scalar_type r2 = x*x+y*y, ratio = -4.*exp(-a4*r2*r2)*a4*r2;
617  res[0] = ratio*x;
618  res[1] = ratio*y;
619  break;
620  }
621  case POLYNOMIAL_CUTOFF: {
622  scalar_type r = gmm::sqrt(x*x+y*y);
623  scalar_type ratio = 0;
624 
625  if ( r > r1 && r < r0 ) {
626  // scalar_type c = 6./(pow(r0,3.) - pow(r1,3.) + 3*r1*r0*(r1-r0));
627  // ratio = c*(r - r0)*(r - r1);
628  ratio = 6.*(r - r0)*(r - r1)/pow(r0-r1, 3.);
629  }
630  res[0] = ratio*x/r;
631  res[1] = ratio*y/r;
632  break;
633  }
634  case POLYNOMIAL2_CUTOFF: {
635  scalar_type r = gmm::sqrt(x*x+y*y);
636  scalar_type ratio = 0;
637  if (r > r1 && r < r0) {
638 // scalar_type c = 1./((-1./30.)*(pow(r1,5) - pow(r0,5))
639 // + (1./6.)*(pow(r1,4)*r0 - r1*pow(r0,4))
640 // - (1./3.)*(pow(r1,3)*pow(r0,2)
641 // - pow(r1,2)*pow(r0,3)));
642 // ratio = - c*gmm::sqr(r-r0)*gmm::sqr(r-r1);
643  ratio = -30.0*gmm::sqr(r-r0)*gmm::sqr(r-r1) / pow(r0-r1, 5.0);
644  }
645  res[0] = ratio*x/r;
646  res[1] = ratio*y/r;
647  break;
648  }
649  default :
650  res[0] = 0;
651  res[1] = 0;
652  }
653  return res;
654  }
655 
656  base_matrix
657  cutoff_xy_function::hess(scalar_type x, scalar_type y) const {
658  base_matrix res(2,2);
659  switch (fun) {
660  case EXPONENTIAL_CUTOFF: {
661  scalar_type r2 = x*x+y*y, r4 = r2*r2;
662  res(0,0) = 4.0*a4*(-3.0*x*x - y*y + 4.0*a4*x*x*r4)*exp(-a4*r4);
663  res(1,0) = 8.0*a4*x*y*(-1.0 + 2.0*a4*r4)*exp(-a4*r4);
664  res(0,1) = res(1,0);
665  res(1,1) = 4.0*a4*(-3.0*y*y - x*x + 4.0*a4*y*y*r4)*exp(-a4*r4);
666  break;
667  }
668  case POLYNOMIAL_CUTOFF: {
669  scalar_type r2 = x*x+y*y, r = gmm::sqrt(r2), c=6./(pow(r0-r1,3.)*r*r2);
670  if ( r > r1 && r < r0 ) {
671  res(0,0) = c*(x*x*r2 + r1*r0*y*y - r*r2*(r0+r1) + r2*r2);
672  res(1,0) = c*x*y*(r2 - r1*r0);
673  res(0,1) = res(1,0);
674  res(1,1) = c*(y*y*r2 + r1*r0*x*x - r*r2*(r0+r1) + r2*r2);
675  }
676  break;
677  }
678  case POLYNOMIAL2_CUTOFF: {
679  scalar_type r2 = x*x+y*y, r = gmm::sqrt(r2), r3 = r*r2;
680  if (r > r1 && r < r0) {
681  scalar_type dp = -30.0*(r1-r)*(r1-r)*(r0-r)*(r0-r) / pow(r0-r1, 5.0);
682  scalar_type ddp = 60.0*(r1-r)*(r0-r)*(r0+r1-2.0*r) / pow(r0-r1, 5.0);
683  scalar_type rx= x/r, ry= y/r, rxx= y*y/r3, rxy= -x*y/r3, ryy= x*x/r3;
684  res(0,0) = ddp*rx*rx + dp*rxx;
685  res(1,0) = ddp*rx*ry + dp*rxy;
686  res(0,1) = res(1,0);
687  res(1,1) = ddp*ry*ry + dp*ryy;
688  }
689  break;
690  }
691  }
692  return res;
693  }
694 
695  cutoff_xy_function::cutoff_xy_function(int fun_num, scalar_type r,
696  scalar_type r1_, scalar_type r0_)
697  {
698  fun = fun_num;
699  r1 = r1_; r0 = r0_;
700  a4 = 0;
701  if (r > 0) a4 = pow(2.7/r,4.);
702  }
703 
704 
705  struct global_function_on_levelsets_2D_ :
706  public global_function, public context_dependencies {
707  const std::vector<level_set> dummy_lsets;
708  const std::vector<level_set> &lsets;
709  const level_set &ls;
710  mutable pmesher_signed_distance mls_x, mls_y;
711  mutable size_type cv;
712 
713  pxy_function fn;
714 
715  void update_mls(size_type cv_, size_type n) const {
716  if (cv_ != cv) {
717  cv=cv_;
718  if (lsets.size() == 0) {
719  mls_x = ls.mls_of_convex(cv, 1);
720  mls_y = ls.mls_of_convex(cv, 0);
721  } else {
722  base_node pt(n);
723  scalar_type d = scalar_type(-2);
724  for (const level_set &ls_ : lsets) {
725  pmesher_signed_distance mls_xx, mls_yy;
726  mls_xx = ls_.mls_of_convex(cv, 1);
727  mls_yy = ls_.mls_of_convex(cv, 0);
728  scalar_type x = (*mls_xx)(pt), y = (*mls_yy)(pt);
729  scalar_type d2 = gmm::sqr(x) + gmm::sqr(y);
730  if (d < scalar_type(-1) || d2 < d) {
731  d = d2;
732  mls_x = mls_xx;
733  mls_y = mls_yy;
734  }
735  }
736  }
737  }
738  }
739 
740  virtual scalar_type val(const fem_interpolation_context& c) const {
741  update_mls(c.convex_num(), c.xref().size());
742  scalar_type x = (*mls_x)(c.xref());
743  scalar_type y = (*mls_y)(c.xref());
744  if (c.xfem_side() > 0 && y <= 1E-13) y = 1E-13;
745  if (c.xfem_side() < 0 && y >= -1E-13) y = -1E-13;
746  return fn->val(x,y);
747  }
748  virtual void grad(const fem_interpolation_context& c,
749  base_small_vector &g) const {
750  size_type P = c.xref().size();
751  base_small_vector dx(P), dy(P), dfr(2);
752 
753  update_mls(c.convex_num(), P);
754  scalar_type x = mls_x->grad(c.xref(), dx);
755  scalar_type y = mls_y->grad(c.xref(), dy);
756  if (c.xfem_side() > 0 && y <= 0) y = 1E-13;
757  if (c.xfem_side() < 0 && y >= 0) y = -1E-13;
758 
759  base_small_vector gfn = fn->grad(x,y);
760  gmm::mult(c.B(), gfn[0]*dx + gfn[1]*dy, g);
761  }
762  virtual void hess(const fem_interpolation_context&c,
763  base_matrix &h) const {
764  size_type P = c.xref().size(), N = c.N();
765  base_small_vector dx(P), dy(P), dfr(2), dx_real(N), dy_real(N);
766 
767  update_mls(c.convex_num(), P);
768  scalar_type x = mls_x->grad(c.xref(), dx);
769  scalar_type y = mls_y->grad(c.xref(), dy);
770  if (c.xfem_side() > 0 && y <= 0) y = 1E-13;
771  if (c.xfem_side() < 0 && y >= 0) y = -1E-13;
772 
773  base_small_vector gfn = fn->grad(x,y);
774  base_matrix hfn = fn->hess(x,y);
775 
776  base_matrix hx, hy, hx_real(N*N, 1), hy_real(N*N, 1);
777  mls_x->hess(c.xref(), hx);
778  mls_x->hess(c.xref(), hy);
779  gmm::reshape(hx, P*P, 1);
780  gmm::reshape(hy, P*P, 1);
781 
782  gmm::mult(c.B3(), hx, hx_real);
783  gmm::mult(c.B32(), gmm::scaled(dx, -1.0), gmm::mat_col(hx_real, 0));
784  gmm::mult(c.B3(), hy, hy_real);
785  gmm::mult(c.B32(), gmm::scaled(dy, -1.0), gmm::mat_col(hy_real, 0));
786  gmm::mult(c.B(), dx, dx_real);
787  gmm::mult(c.B(), dy, dy_real);
788 
789  for (size_type i = 0; i < N; ++i)
790  for (size_type j = 0; j < N; ++j) {
791  h(i, j) = hfn(0,0) * dx_real[i] * dx_real[j]
792  + hfn(0,1) * dx_real[i] * dy_real[j]
793  + hfn(1,0) * dy_real[i] * dx_real[j]
794  + hfn(1,1) * dy_real[i] * dy_real[j]
795  + gfn[0] * hx_real(i * N + j, 0) + gfn[1] * hy_real(i*N+j,0);
796  }
797  }
798 
799  void update_from_context() const { cv = size_type(-1); }
800 
801  global_function_on_levelsets_2D_(const std::vector<level_set> &lsets_,
802  const pxy_function &fn_)
803  : global_function(2), dummy_lsets(0, dummy_level_set()),
804  lsets(lsets_), ls(dummy_level_set()), fn(fn_) {
805  GMM_ASSERT1(lsets.size() > 0, "The list of level sets should"
806  " contain at least one level set.");
807  cv = size_type(-1);
808  for (const level_set &ls_ : lsets)
809  this->add_dependency(ls_);
810  }
811 
812  global_function_on_levelsets_2D_(const level_set &ls_,
813  const pxy_function &fn_)
814  : global_function(2), dummy_lsets(0, dummy_level_set()),
815  lsets(dummy_lsets), ls(ls_), fn(fn_) {
816  cv = size_type(-1);
817  this->add_dependency(ls);
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 
835 
836 
837  // Global function for bspline basis
838  const scalar_type eps(1e-12);
839 
840  // order k = 3
841  scalar_type bsp3_01(scalar_type t) {
842  return (t >= -eps && t < 1) ? pow(1.-t,2)
843  : 0;
844  }
845  scalar_type bsp3_01_der(scalar_type t) {
846  return (t >= -eps && t < 1) ? 2.*t-2.
847  : 0;
848  }
849  scalar_type bsp3_01_der2(scalar_type t) {
850  return (t >= eps && t <= 1.-eps) ? 2.
851  : 0;
852  }
853  scalar_type bsp3_01_der2_with_hint(scalar_type t, scalar_type t_hint) {
854  return (t > -eps && t < 1.+eps && t_hint > 0 && t_hint < 1) ? 2.
855  : 0;
856  }
857  scalar_type bsp3_02(scalar_type t) {
858  if (t >= -eps) {
859  if (t < 1)
860  return (-1.5*t+2.)*t;
861  else if (t < 2)
862  return 0.5*pow(2.-t,2);
863  }
864  return 0;
865  }
866  scalar_type bsp3_02_der(scalar_type t) {
867  if (t >= -eps) {
868  if (t < 1)
869  return -3.*t+2.;
870  else if (t < 2)
871  return t-2.;
872  }
873  return 0;
874  }
875  scalar_type bsp3_02_der2(scalar_type t) {
876  if (t >= eps) {
877  if (t < 1.-eps)
878  return -3.;
879  else if (t < 1.+eps)
880  return 0;
881  else if (t <= 2.-eps)
882  return 1.;
883  }
884  return 0;
885  }
886  scalar_type bsp3_03(scalar_type t) {
887  if (t >= -eps) {
888  if (t < 1) {
889  return 0.5*t*t;
890  } else if (t < 2) {
891  t -= 1.;
892  return t*(1-t)+0.5;
893  } else if (t < 3) {
894  t = 3.-t;
895  return 0.5*t*t;
896  }
897  }
898  return 0;
899  }
900  scalar_type bsp3_03_der(scalar_type t) {
901  if (t >= -eps) {
902  if (t < 1) {
903  return t;
904  } else if (t < 2) {
905  t -= 1.;
906  return 1.-2.*t;
907  } else if (t < 3) {
908  return t-3.;
909  }
910  }
911  return 0;
912  }
913  scalar_type bsp3_03_der2(scalar_type t) {
914  if (t >= -eps) {
915  if (t < eps)
916  return 0;
917  else if (t <= 1.-eps)
918  return 1.;
919  else if (t < 1.+eps)
920  return 0;
921  else if (t <= 2.-eps)
922  return -2.;
923  else if (t < 2.+eps)
924  return 0;
925  else if (t <= 3.-eps)
926  return 1.;
927  else if (t < 3.+eps)
928  return 0;
929  }
930  return 0;
931  }
932 
933  // order k = 4
934  scalar_type bsp4_01(scalar_type t) {
935  return (t >= -eps && t < 1) ? pow(1.-t,3)
936  : 0;
937  }
938  scalar_type bsp4_01_der(scalar_type t) {
939  return (t >= -eps && t < 1) ? -3*pow(1.-t,2)
940  : 0;
941  }
942  scalar_type bsp4_01_der2(scalar_type t) {
943  return (t >= -eps && t < 1) ? 6*(1.-t)
944  : 0;
945  }
946  scalar_type bsp4_02(scalar_type t) {
947  if (t >= -eps) {
948  if (t < 1) {
949  return ((7./4.*t-9./2.)*t+3.)*t;
950  } else if (t < 2) {
951  return 1./4.*pow(2.-t,3);
952  }
953  }
954  return 0;
955  }
956  scalar_type bsp4_02_der(scalar_type t) {
957  if (t >= -eps) {
958  if (t < 1) {
959  return (21./4.*t-9.)*t+3.;
960  } else if (t < 2) {
961  return -3./4.*pow(2.-t,2);
962  }
963  }
964  return 0;
965  }
966  scalar_type bsp4_02_der2(scalar_type t) {
967  if (t >= -eps) {
968  if (t < 1) {
969  return 21./2.*t-9.;
970  } else if (t < 2) {
971  return 3./2.*(2.-t);
972  }
973  }
974  return 0;
975  }
976  scalar_type bsp4_03(scalar_type t) {
977  if (t >= -eps) {
978  if (t < 1) {
979  return (-11./12.*t+3./2.)*t*t;
980  } else if (t < 2) {
981  t -= 1;
982  return ((7./12.*t-5./4.)*t+1./4.)*t+7./12.;
983  } else if (t < 3) {
984  return 1./6.*pow(3.-t,3);
985  }
986  }
987  return 0;
988  }
989  scalar_type bsp4_03_der(scalar_type t) {
990  if (t >= -eps) {
991  if (t < 1) {
992  return (-11./4.*t+3.)*t;
993  } else if (t < 2) {
994  t -= 1;
995  return (7./4.*t-5./2.)*t+1./4.;
996  } else if (t < 3) {
997  return -1./2.*pow(3.-t,2);
998  }
999  }
1000  return 0;
1001  }
1002  scalar_type bsp4_03_der2(scalar_type t) {
1003  if (t >= -eps) {
1004  if (t < 1) {
1005  return -11./2.*t+3.;
1006  } else if (t < 2) {
1007  t -= 1;
1008  return 7./2.*t-5./2.;
1009  } else if (t < 3) {
1010  return 3.-t;
1011  }
1012  }
1013  return 0;
1014  }
1015  scalar_type bsp4_04(scalar_type t) {
1016  if (t > 2) {
1017  t = 4.-t;
1018  }
1019  if (t >= -eps) {
1020  if (t < 1) {
1021  return 1./6.*pow(t,3);
1022  } else if (t <= 2) {
1023  t = t-1.;
1024  return ((-1./2.*t+1./2.)*t+1./2.)*t+1./6.;
1025  }
1026  }
1027  return 0;
1028  }
1029  scalar_type bsp4_04_der(scalar_type t) {
1030  scalar_type sgn(1);
1031  if (t > 2) {
1032  t = 4.-t;
1033  sgn = -1.;
1034  }
1035  if (t >= -eps) {
1036  if (t < 1) {
1037  return 1./2.*pow(t,2)*sgn;
1038  } else if (t < 2) {
1039  t = t-1.;
1040  return ((-3./2.*t+1.)*t+1./2.)*sgn;
1041  }
1042  }
1043  return 0;
1044  }
1045  scalar_type bsp4_04_der2(scalar_type t) {
1046  if (t > 2) {
1047  t = 4.-t;
1048  }
1049  if (t >= -eps) {
1050  if (t < 1) {
1051  return t;
1052  } else if (t < 2) {
1053  t = t-1.;
1054  return -3.*t+1.;
1055  }
1056  }
1057  return 0;
1058  }
1059 
1060 
1061  // order k = 5
1062  scalar_type bsp5_01(scalar_type t) {
1063  return (t >= -eps && t < 1) ? pow(1.-t,4)
1064  : 0;
1065  }
1066  scalar_type bsp5_01_der(scalar_type t) {
1067  return (t >= -eps && t < 1) ? -4.*pow(1.-t,3)
1068  : 0;
1069  }
1070  scalar_type bsp5_01_der2(scalar_type t) {
1071  return (t >= -eps && t < 1) ? 12.*pow(1.-t,2)
1072  : 0;
1073  }
1074  scalar_type bsp5_02(scalar_type t) {
1075  if (t >= -eps) {
1076  if (t < 1) {
1077  return (((-15./8.*t+7.)*t-9.)*t+4.)*t;
1078  } else if (t < 2) {
1079  return 1./8.*pow(2.-t,4);
1080  }
1081  }
1082  return 0;
1083  }
1084  scalar_type bsp5_02_der(scalar_type t) {
1085  if (t >= -eps) {
1086  if (t < 1) {
1087  return ((-15./2.*t+21.)*t-18.)*t+4.;
1088  } else if (t < 2) {
1089  return -1./2.*pow(2.-t,3);
1090  }
1091  }
1092  return 0;
1093  }
1094  scalar_type bsp5_02_der2(scalar_type t) {
1095  if (t >= -eps) {
1096  if (t < 1) {
1097  return (-45./2.*t+42.)*t-18.;
1098  } else if (t < 2) {
1099  return 3./2.*pow(2.-t,2);
1100  }
1101  }
1102  return 0;
1103  }
1104  scalar_type bsp5_03(scalar_type t) {
1105  if (t >= -eps) {
1106  if (t < 1) {
1107  return ((85./72.*t-11./3.)*t+3.)*t*t;
1108  } else if (t < 2) {
1109  t = 2.-t;
1110  return (((-23./72*t+2./9.)*t+1./3.)*t+2./9.)*t+1./18.;
1111  } else if (t < 3) {
1112  return 1./18.*pow(3.-t,4);
1113  }
1114  }
1115  return 0;
1116  }
1117  scalar_type bsp5_03_der(scalar_type t) {
1118  if (t >= -eps) {
1119  if (t < 1) {
1120  return ((85./18.*t-11.)*t+6.)*t;
1121  } else if (t < 2) {
1122  t = 2.-t;
1123  return -(((-23./18*t+2./3.)*t+2./3.)*t+2./9.);
1124  } else if (t < 3) {
1125  return -2./9.*pow(3.-t,3);
1126  }
1127  }
1128  return 0;
1129  }
1130  scalar_type bsp5_03_der2(scalar_type t) {
1131  if (t >= -eps) {
1132  if (t < 1) {
1133  return (85./6.*t-22.)*t+6.;
1134  } else if (t < 2) {
1135  t = 2.-t;
1136  return (-23./6*t+4./3.)*t+2./3.;
1137  } else if (t < 3) {
1138  return 2./3.*pow(3.-t,2);
1139  }
1140  }
1141  return 0;
1142  }
1143  scalar_type bsp5_04(scalar_type t) {
1144  if (t >= -eps) {
1145  if (t < 1) {
1146  return (-25./72.*t+2./3.)*t*t*t;
1147  } else if (t < 2) {
1148  t = 2.-t;
1149  return (((23./72.*t-5./9.)*t-1./3.)*t+4./9.)*t+4./9.;
1150  } else if (t < 3) {
1151  t = t-2.;
1152  return (((-13./72.*t+5./9.)*t-1./3.)*t-4./9.)*t+4./9.;
1153  } else if (t < 4) {
1154  return 1./24.*pow(4.-t,4);
1155  }
1156  }
1157  return 0;
1158  }
1159  scalar_type bsp5_04_der(scalar_type t) {
1160  if (t >= -eps) {
1161  if (t < 1) {
1162  return (-25./18.*t+2.)*t*t;
1163  } else if (t < 2) {
1164  t = 2.-t;
1165  return -(((23./18.*t-5./3.)*t-2./3.)*t+4./9.);
1166  } else if (t < 3) {
1167  t = t-2.;
1168  return ((-13./18.*t+5./3.)*t-2./3.)*t-4./9.;
1169  } else if (t < 4) {
1170  return -1./6.*pow(4.-t,3);
1171  }
1172  }
1173  return 0;
1174  }
1175  scalar_type bsp5_04_der2(scalar_type t) {
1176  if (t >= -eps) {
1177  if (t < 1) {
1178  return (-25./6.*t+4.)*t;
1179  } else if (t < 2) {
1180  t = 2.-t;
1181  return (23./6.*t-10./3.)*t-2./3.;
1182  } else if (t < 3) {
1183  t = t-2.;
1184  return -((-13./6.*t+10./3.)*t-2./3.);
1185  } else if (t < 4) {
1186  return 1./2.*pow(4.-t,2);
1187  }
1188  }
1189  return 0;
1190  }
1191  scalar_type bsp5_05(scalar_type t) {
1192  if (t > scalar_type(2.5)) {
1193  t = 5.-t;
1194  }
1195  if (t >= -eps) {
1196  if (t < 1) {
1197  return 1./24.*pow(t,4);
1198  } else if (t < 2) {
1199  t = t-1.;
1200  return (((-1./6.*t+1./6.)*t+1./4.)*t+1./6.)*t+1./24.;
1201  } else if (t < 3) {
1202  t = pow(t-2.5,2);
1203  return (1./4.*t-5./8.)*t+115./192.;
1204  }
1205  }
1206  return 0;
1207  }
1208  scalar_type bsp5_05_der(scalar_type t) {
1209  scalar_type sgn(1);
1210  if (t > scalar_type(2.5)) {
1211  t = 5.-t;
1212  sgn = -1.;
1213  }
1214  if (t >= -eps) {
1215  if (t < 1) {
1216  return 1./6.*pow(t,3)*sgn;
1217  } else if (t < 2) {
1218  t = t-1.;
1219  return (((-2./3.*t+1./2.)*t+1./2)*t+1./6.)*sgn;
1220  } else if (t < 3) {
1221  t = t-2.5;
1222  return (t*t-5./4.)*t*sgn;
1223  }
1224  }
1225  return 0;
1226  }
1227  scalar_type bsp5_05_der2(scalar_type t) {
1228  if (t > scalar_type(2.5)) {
1229  t = 5.-t;
1230  }
1231  if (t >= -eps) {
1232  if (t < 1) {
1233  return 1./2.*pow(t,2);
1234  } else if (t < 2) {
1235  t = t-1.;
1236  return (-2.*t+1.)*t+1./2;
1237  } else if (t < 3) {
1238  t = t-2.5;
1239  return 3*t*t-5./4.;
1240  }
1241  }
1242  return 0;
1243  }
1244 
1245 
1246 
1247  class global_function_x_bspline_
1248  : public global_function_simple, public context_dependencies {
1249  scalar_type xmin, xmax, xscale;
1250  std::function<scalar_type(scalar_type)> fx, fx_der, fx_der2;
1251  public:
1252  void update_from_context() const {}
1253 
1254  virtual scalar_type val(const base_node &pt) const
1255  {
1256  return fx(xscale*(pt[0]-xmin));
1257  }
1258  virtual void grad(const base_node &pt, base_small_vector &g) const
1259  {
1260  scalar_type dx = xscale*(pt[0]-xmin);
1261  g.resize(dim_);
1262  g[0] = xscale * fx_der(dx);
1263  }
1264  virtual void hess(const base_node &pt, base_matrix &h) const
1265  {
1266  scalar_type dx = xscale*(pt[0]-xmin);
1267  h.resize(dim_, dim_);
1268  gmm::clear(h);
1269  h(0,0) = xscale*xscale * fx_der2(dx);
1270  }
1271 
1272  virtual bool is_in_support(const base_node &pt) const {
1273  scalar_type dx = pt[0]-(xmin+xmax)/2;
1274  return (gmm::abs(dx)+1e-9 < gmm::abs(xmax-xmin)/2);
1275  }
1276 
1277  virtual void bounding_box(base_node &bmin, base_node &bmax) const {
1278  GMM_ASSERT1(bmin.size() == dim_ && bmax.size() == dim_,
1279  "Wrong dimensions");
1280  bmin[0] = std::min(xmin,xmax);
1281  bmax[0] = std::max(xmin,xmax);
1282  }
1283 
1284  global_function_x_bspline_(const scalar_type &xmin_, const scalar_type &xmax_,
1285  const size_type &order, const size_type &xtype)
1286  : global_function_simple(1), xmin(xmin_), xmax(xmax_),
1287  xscale(scalar_type(xtype)/(xmax-xmin))
1288  {
1289  GMM_ASSERT1(order >= 3 && order <= 5, "Only orders 3 to 5 are supported");
1290  GMM_ASSERT1(xtype >= 1 && xtype <= order, "Wrong input");
1291  if (order == 3) {
1292  if (xtype == 1) {
1293  fx = bsp3_01; fx_der = bsp3_01_der; fx_der2 = bsp3_01_der2;
1294  } else if (xtype == 2) {
1295  fx = bsp3_02; fx_der = bsp3_02_der; fx_der2 = bsp3_02_der2;
1296  } else if (xtype == 3) {
1297  fx = bsp3_03; fx_der = bsp3_03_der; fx_der2 = bsp3_03_der2;
1298  }
1299  } else if (order == 4) {
1300  if (xtype == 1) {
1301  fx = bsp4_01; fx_der = bsp4_01_der; fx_der2 = bsp4_01_der2;
1302  } else if (xtype == 2) {
1303  fx = bsp4_02; fx_der = bsp4_02_der; fx_der2 = bsp4_02_der2;
1304  } else if (xtype == 3) {
1305  fx = bsp4_03; fx_der = bsp4_03_der; fx_der2 = bsp4_03_der2;
1306  } else if (xtype == 4) {
1307  fx = bsp4_04; fx_der = bsp4_04_der; fx_der2 = bsp4_04_der2;
1308  }
1309  } else if (order == 5) {
1310  if (xtype == 1) {
1311  fx = bsp5_01; fx_der = bsp5_01_der; fx_der2 = bsp5_01_der2;
1312  } else if (xtype == 2) {
1313  fx = bsp5_02; fx_der = bsp5_02_der; fx_der2 = bsp5_02_der2;
1314  } else if (xtype == 3) {
1315  fx = bsp5_03; fx_der = bsp5_03_der; fx_der2 = bsp5_03_der2;
1316  } else if (xtype == 4) {
1317  fx = bsp5_04; fx_der = bsp5_04_der; fx_der2 = bsp5_04_der2;
1318  } else if (xtype == 5) {
1319  fx = bsp5_05; fx_der = bsp5_05_der; fx_der2 = bsp5_05_der2;
1320  }
1321  }
1322  }
1323 
1324  virtual ~global_function_x_bspline_()
1325  { DAL_STORED_OBJECT_DEBUG_DESTROYED(this, "Global function x bspline"); }
1326  };
1327 
1328 
1329 
1330  class global_function_xy_bspline_
1331  : public global_function_simple, public context_dependencies {
1332  scalar_type xmin, ymin, xmax, ymax, xscale, yscale;
1333  std::function<scalar_type(scalar_type)>
1334  fx, fy, fx_der, fy_der, fx_der2, fy_der2;
1335  public:
1336  void update_from_context() const {}
1337 
1338  virtual scalar_type val(const base_node &pt) const
1339  {
1340  return fx(xscale*(pt[0]-xmin))
1341  * fy(yscale*(pt[1]-ymin));
1342  }
1343  virtual void grad(const base_node &pt, base_small_vector &g) const
1344  {
1345  scalar_type dx = xscale*(pt[0]-xmin),
1346  dy = yscale*(pt[1]-ymin);
1347  g.resize(dim_);
1348  g[0] = xscale * fx_der(dx) * fy(dy);
1349  g[1] = fx(dx) * yscale * fy_der(dy);
1350  }
1351  virtual void hess(const base_node &pt, base_matrix &h) const
1352  {
1353  scalar_type dx = xscale*(pt[0]-xmin),
1354  dy = yscale*(pt[1]-ymin);
1355  h.resize(dim_, dim_);
1356  gmm::clear(h);
1357  h(0,0) = xscale*xscale * fx_der2(dx) * fy(dy);
1358  h(0,1) = h(1,0) = xscale * fx_der(dx) * yscale * fy_der(dy);
1359  h(1,1) = fx(dx) * yscale*yscale * fy_der2(dy);
1360  }
1361 
1362  virtual bool is_in_support(const base_node &pt) const {
1363  scalar_type dx = pt[0]-(xmin+xmax)/2,
1364  dy = pt[1]-(ymin+ymax)/2;
1365  return (gmm::abs(dx)+1e-9 < gmm::abs(xmax-xmin)/2) &&
1366  (gmm::abs(dy)+1e-9 < gmm::abs(ymax-ymin)/2);
1367  }
1368 
1369  virtual void bounding_box(base_node &bmin, base_node &bmax) const {
1370  GMM_ASSERT1(bmin.size() == dim_ && bmax.size() == dim_,
1371  "Wrong dimensions");
1372  bmin[0] = std::min(xmin,xmax);
1373  bmin[1] = std::min(ymin,ymax);
1374  bmax[0] = std::max(xmin,xmax);
1375  bmax[1] = std::max(ymin,ymax);
1376  }
1377 
1378  global_function_xy_bspline_(const scalar_type &xmin_, const scalar_type &xmax_,
1379  const scalar_type &ymin_, const scalar_type &ymax_,
1380  const size_type &order,
1381  const size_type &xtype, const size_type &ytype)
1382  : global_function_simple(2), xmin(xmin_), ymin(ymin_), xmax(xmax_), ymax(ymax_),
1383  xscale(scalar_type(xtype)/(xmax-xmin)), yscale(scalar_type(ytype)/(ymax-ymin))
1384  {
1385  GMM_ASSERT1(order >= 3 && order <= 5, "Wrong input");
1386  GMM_ASSERT1(xtype >= 1 && xtype <= order &&
1387  ytype >= 1 && ytype <= order, "Wrong input");
1388  if (order == 3) {
1389  if (xtype == 1) {
1390  fx = bsp3_01; fx_der = bsp3_01_der; fx_der2 = bsp3_01_der2;
1391  } else if (xtype == 2) {
1392  fx = bsp3_02; fx_der = bsp3_02_der; fx_der2 = bsp3_02_der2;
1393  } else if (xtype == 3) {
1394  fx = bsp3_03; fx_der = bsp3_03_der; fx_der2 = bsp3_03_der2;
1395  }
1396 
1397  if (ytype == 1) {
1398  fy = bsp3_01; fy_der = bsp3_01_der; fy_der2 = bsp3_01_der2;
1399  } else if (ytype == 2) {
1400  fy = bsp3_02; fy_der = bsp3_02_der; fy_der2 = bsp3_02_der2;
1401  } else if (ytype == 3) {
1402  fy = bsp3_03; fy_der = bsp3_03_der; fy_der2 = bsp3_03_der2;
1403  }
1404  } else if (order == 4) {
1405  if (xtype == 1) {
1406  fx = bsp4_01; fx_der = bsp4_01_der; fx_der2 = bsp4_01_der2;
1407  } else if (xtype == 2) {
1408  fx = bsp4_02; fx_der = bsp4_02_der; fx_der2 = bsp4_02_der2;
1409  } else if (xtype == 3) {
1410  fx = bsp4_03; fx_der = bsp4_03_der; fx_der2 = bsp4_03_der2;
1411  } else if (xtype == 4) {
1412  fx = bsp4_04; fx_der = bsp4_04_der; fx_der2 = bsp4_04_der2;
1413  }
1414 
1415  if (ytype == 1) {
1416  fy = bsp4_01; fy_der = bsp4_01_der; fy_der2 = bsp4_01_der2;
1417  } else if (ytype == 2) {
1418  fy = bsp4_02; fy_der = bsp4_02_der; fy_der2 = bsp4_02_der2;
1419  } else if (ytype == 3) {
1420  fy = bsp4_03; fy_der = bsp4_03_der; fy_der2 = bsp4_03_der2;
1421  } else if (ytype == 4) {
1422  fy = bsp4_04; fy_der = bsp4_04_der; fy_der2 = bsp4_04_der2;
1423  }
1424 
1425  } else if (order == 5) {
1426  if (xtype == 1) {
1427  fx = bsp5_01; fx_der = bsp5_01_der; fx_der2 = bsp5_01_der2;
1428  } else if (xtype == 2) {
1429  fx = bsp5_02; fx_der = bsp5_02_der; fx_der2 = bsp5_02_der2;
1430  } else if (xtype == 3) {
1431  fx = bsp5_03; fx_der = bsp5_03_der; fx_der2 = bsp5_03_der2;
1432  } else if (xtype == 4) {
1433  fx = bsp5_04; fx_der = bsp5_04_der; fx_der2 = bsp5_04_der2;
1434  } else if (xtype == 5) {
1435  fx = bsp5_05; fx_der = bsp5_05_der; fx_der2 = bsp5_05_der2;
1436  }
1437 
1438  if (ytype == 1) {
1439  fy = bsp5_01; fy_der = bsp5_01_der; fy_der2 = bsp5_01_der2;
1440  } else if (ytype == 2) {
1441  fy = bsp5_02; fy_der = bsp5_02_der; fy_der2 = bsp5_02_der2;
1442  } else if (ytype == 3) {
1443  fy = bsp5_03; fy_der = bsp5_03_der; fy_der2 = bsp5_03_der2;
1444  } else if (ytype == 4) {
1445  fy = bsp5_04; fy_der = bsp5_04_der; fy_der2 = bsp5_04_der2;
1446  } else if (ytype == 5) {
1447  fy = bsp5_05; fy_der = bsp5_05_der; fy_der2 = bsp5_05_der2;
1448  }
1449  }
1450  }
1451 
1452  virtual ~global_function_xy_bspline_()
1453  { DAL_STORED_OBJECT_DEBUG_DESTROYED(this, "Global function xy bspline"); }
1454  };
1455 
1456 
1457  class global_function_xyz_bspline_
1458  : public global_function_simple, public context_dependencies {
1459  scalar_type xmin, ymin, zmin, xmax, ymax, zmax, xscale, yscale, zscale;
1460  std::function<scalar_type(scalar_type)>
1461  fx, fy, fz, fx_der, fy_der, fz_der, fx_der2, fy_der2, fz_der2;
1462  public:
1463  void update_from_context() const {}
1464 
1465  virtual scalar_type val(const base_node &pt) const
1466  {
1467  return fx(xscale*(pt[0]-xmin))
1468  * fy(yscale*(pt[1]-ymin))
1469  * fz(zscale*(pt[2]-zmin));
1470  }
1471  virtual void grad(const base_node &pt, base_small_vector &g) const
1472  {
1473  scalar_type dx = xscale*(pt[0]-xmin),
1474  dy = yscale*(pt[1]-ymin),
1475  dz = zscale*(pt[2]-zmin);
1476  g.resize(dim_);
1477  g[0] = xscale * fx_der(dx) * fy(dy) * fz(dz);
1478  g[1] = fx(dx) * yscale * fy_der(dy) * fz(dz);
1479  g[2] = fx(dx) * fy(dy) * zscale * fz_der(dz);
1480  }
1481  virtual void hess(const base_node &pt, base_matrix &h) const
1482  {
1483  scalar_type dx = xscale*(pt[0]-xmin),
1484  dy = yscale*(pt[1]-ymin),
1485  dz = zscale*(pt[2]-zmin);
1486  h.resize(dim_, dim_);
1487  gmm::clear(h);
1488  h(0,0) = xscale*xscale * fx_der2(dx) * fy(dy) * fz(dz);
1489  h(1,1) = fx(dx) * yscale*yscale * fy_der2(dy) * fz(dz);
1490  h(2,2) = fx(dx) * fy(dy) * zscale*zscale * fz_der2(dz);
1491  h(0,1) = h(1,0) = xscale * fx_der(dx) * yscale * fy_der(dy) * fz(dz);
1492  h(0,2) = h(2,0) = xscale * fx_der(dx) * fy(dy) * zscale * fz_der(dz);
1493  h(1,2) = h(2,1) = fx(dx) * yscale * fy_der(dy) * zscale * fz_der(dz);
1494  }
1495 
1496  virtual bool is_in_support(const base_node &pt) const {
1497  scalar_type dx = pt[0]-(xmin+xmax)/2,
1498  dy = pt[1]-(ymin+ymax)/2,
1499  dz = pt[2]-(zmin+zmax)/2;
1500  return (gmm::abs(dx)+1e-9 < gmm::abs(xmax-xmin)/2) &&
1501  (gmm::abs(dy)+1e-9 < gmm::abs(ymax-ymin)/2) &&
1502  (gmm::abs(dz)+1e-9 < gmm::abs(zmax-zmin)/2);
1503  }
1504 
1505  virtual void bounding_box(base_node &bmin, base_node &bmax) const {
1506  GMM_ASSERT1(bmin.size() == dim_ && bmax.size() == dim_,
1507  "Wrong dimensions");
1508  bmin[0] = std::min(xmin,xmax);
1509  bmin[1] = std::min(ymin,ymax);
1510  bmin[2] = std::min(zmin,zmax);
1511  bmax[0] = std::max(xmin,xmax);
1512  bmax[1] = std::max(ymin,ymax);
1513  bmax[2] = std::max(zmin,zmax);
1514  }
1515 
1516  global_function_xyz_bspline_(const scalar_type &xmin_, const scalar_type &xmax_,
1517  const scalar_type &ymin_, const scalar_type &ymax_,
1518  const scalar_type &zmin_, const scalar_type &zmax_,
1519  const size_type &order,
1520  const size_type &xtype,
1521  const size_type &ytype,
1522  const size_type &ztype)
1523  : global_function_simple(3), xmin(xmin_), ymin(ymin_), zmin(zmin_),
1524  xmax(xmax_), ymax(ymax_), zmax(zmax_),
1525  xscale(scalar_type(xtype)/(xmax-xmin)),
1526  yscale(scalar_type(ytype)/(ymax-ymin)),
1527  zscale(scalar_type(ztype)/(zmax-zmin))
1528  {
1529  GMM_ASSERT1(order >= 3 && order <= 5, "Wrong input");
1530  GMM_ASSERT1(xtype >= 1 && xtype <= order &&
1531  ytype >= 1 && ytype <= order &&
1532  ztype >= 1 && ztype <= order, "Wrong input");
1533  if (order == 3) {
1534 
1535  if (xtype == 1) {
1536  fx = bsp3_01; fx_der = bsp3_01_der; fx_der2 = bsp3_01_der2;
1537  } else if (xtype == 2) {
1538  fx = bsp3_02; fx_der = bsp3_02_der; fx_der2 = bsp3_02_der2;
1539  } else if (xtype == 3) {
1540  fx = bsp3_03; fx_der = bsp3_03_der; fx_der2 = bsp3_03_der2;
1541  }
1542 
1543  if (ytype == 1) {
1544  fy = bsp3_01; fy_der = bsp3_01_der; fy_der2 = bsp3_01_der2;
1545  } else if (ytype == 2) {
1546  fy = bsp3_02; fy_der = bsp3_02_der; fy_der2 = bsp3_02_der2;
1547  } else if (ytype == 3) {
1548  fy = bsp3_03; fy_der = bsp3_03_der; fy_der2 = bsp3_03_der2;
1549  }
1550 
1551  if (ztype == 1) {
1552  fz = bsp3_01; fz_der = bsp3_01_der; fz_der2 = bsp3_01_der2;
1553  } else if (ztype == 2) {
1554  fz = bsp3_02; fz_der = bsp3_02_der; fz_der2 = bsp3_02_der2;
1555  } else if (ztype == 3) {
1556  fz = bsp3_03; fz_der = bsp3_03_der; fz_der2 = bsp3_03_der2;
1557  }
1558 
1559  } else if (order == 4) {
1560 
1561  if (xtype == 1) {
1562  fx = bsp4_01; fx_der = bsp4_01_der; fx_der2 = bsp4_01_der2;
1563  } else if (xtype == 2) {
1564  fx = bsp4_02; fx_der = bsp4_02_der; fx_der2 = bsp4_02_der2;
1565  } else if (xtype == 3) {
1566  fx = bsp4_03; fx_der = bsp4_03_der; fx_der2 = bsp4_03_der2;
1567  } else if (xtype == 4) {
1568  fx = bsp4_04; fx_der = bsp4_04_der; fx_der2 = bsp4_04_der2;
1569  }
1570 
1571  if (ytype == 1) {
1572  fy = bsp4_01; fy_der = bsp4_01_der; fy_der2 = bsp4_01_der2;
1573  } else if (ytype == 2) {
1574  fy = bsp4_02; fy_der = bsp4_02_der; fy_der2 = bsp4_02_der2;
1575  } else if (ytype == 3) {
1576  fy = bsp4_03; fy_der = bsp4_03_der; fy_der2 = bsp4_03_der2;
1577  } else if (ytype == 4) {
1578  fy = bsp4_04; fy_der = bsp4_04_der; fy_der2 = bsp4_04_der2;
1579  }
1580 
1581  if (ztype == 1) {
1582  fz = bsp4_01; fz_der = bsp4_01_der; fz_der2 = bsp4_01_der2;
1583  } else if (ztype == 2) {
1584  fz = bsp4_02; fz_der = bsp4_02_der; fz_der2 = bsp4_02_der2;
1585  } else if (ztype == 3) {
1586  fz = bsp4_03; fz_der = bsp4_03_der; fz_der2 = bsp4_03_der2;
1587  } else if (ztype == 4) {
1588  fz = bsp4_04; fz_der = bsp4_04_der; fz_der2 = bsp4_04_der2;
1589  }
1590 
1591  } else if (order == 5) {
1592 
1593  if (xtype == 1) {
1594  fx = bsp5_01; fx_der = bsp5_01_der; fx_der2 = bsp5_01_der2;
1595  } else if (xtype == 2) {
1596  fx = bsp5_02; fx_der = bsp5_02_der; fx_der2 = bsp5_02_der2;
1597  } else if (xtype == 3) {
1598  fx = bsp5_03; fx_der = bsp5_03_der; fx_der2 = bsp5_03_der2;
1599  } else if (xtype == 4) {
1600  fx = bsp5_04; fx_der = bsp5_04_der; fx_der2 = bsp5_04_der2;
1601  } else if (xtype == 5) {
1602  fx = bsp5_05; fx_der = bsp5_05_der; fx_der2 = bsp5_05_der2;
1603  }
1604 
1605  if (ytype == 1) {
1606  fy = bsp5_01; fy_der = bsp5_01_der; fy_der2 = bsp5_01_der2;
1607  } else if (ytype == 2) {
1608  fy = bsp5_02; fy_der = bsp5_02_der; fy_der2 = bsp5_02_der2;
1609  } else if (ytype == 3) {
1610  fy = bsp5_03; fy_der = bsp5_03_der; fy_der2 = bsp5_03_der2;
1611  } else if (ytype == 4) {
1612  fy = bsp5_04; fy_der = bsp5_04_der; fy_der2 = bsp5_04_der2;
1613  } else if (ytype == 5) {
1614  fy = bsp5_05; fy_der = bsp5_05_der; fy_der2 = bsp5_05_der2;
1615  }
1616 
1617  if (ztype == 1) {
1618  fz = bsp5_01; fz_der = bsp5_01_der; fz_der2 = bsp5_01_der2;
1619  } else if (ztype == 2) {
1620  fz = bsp5_02; fz_der = bsp5_02_der; fz_der2 = bsp5_02_der2;
1621  } else if (ztype == 3) {
1622  fz = bsp5_03; fz_der = bsp5_03_der; fz_der2 = bsp5_03_der2;
1623  } else if (ztype == 4) {
1624  fz = bsp5_04; fz_der = bsp5_04_der; fz_der2 = bsp5_04_der2;
1625  } else if (ztype == 5) {
1626  fz = bsp5_05; fz_der = bsp5_05_der; fz_der2 = bsp5_05_der2;
1627  }
1628 
1629  }
1630  }
1631 
1632  virtual ~global_function_xyz_bspline_()
1633  { DAL_STORED_OBJECT_DEBUG_DESTROYED(this, "Global function xyz bspline"); }
1634  };
1635 
1636 
1637  pglobal_function
1638  global_function_bspline(const scalar_type xmin, const scalar_type xmax,
1639  const size_type order, const size_type xtype) {
1640  return std::make_shared<global_function_x_bspline_>
1641  (xmin, xmax, order, xtype);
1642  }
1643 
1644  pglobal_function
1645  global_function_bspline(const scalar_type xmin, const scalar_type xmax,
1646  const scalar_type ymin, const scalar_type ymax,
1647  const size_type order,
1648  const size_type xtype, const size_type ytype) {
1649  return std::make_shared<global_function_xy_bspline_>
1650  (xmin, xmax, ymin, ymax, order, xtype, ytype);
1651  }
1652 
1653  pglobal_function
1654  global_function_bspline(const scalar_type xmin, const scalar_type xmax,
1655  const scalar_type ymin, const scalar_type ymax,
1656  const scalar_type zmin, const scalar_type zmax,
1657  const size_type order,
1658  const size_type xtype,
1659  const size_type ytype,
1660  const size_type ztype) {
1661  return std::make_shared<global_function_xyz_bspline_>
1662  (xmin, xmax, ymin, ymax, zmin, zmax, order,
1663  xtype, ytype, ztype);
1664  }
1665 
1666 
1667 
1668  // interpolator on mesh_fem
1669 
1670  interpolator_on_mesh_fem::interpolator_on_mesh_fem
1671  (const mesh_fem &mf_, const std::vector<scalar_type> &U_)
1672  : mf(mf_), U(U_) {
1673 
1674  if (mf.is_reduced()) {
1675  gmm::resize(U, mf.nb_basic_dof());
1676  gmm::mult(mf.extension_matrix(), U_, U);
1677  }
1678  base_node bmin, bmax;
1679  scalar_type EPS=1E-13;
1680  cv_stored = size_type(-1);
1681  boxtree.clear();
1682  for (dal::bv_visitor cv(mf.convex_index()); !cv.finished(); ++cv) {
1683  bounding_box(bmin, bmax, mf.linked_mesh().points_of_convex(cv),
1684  mf.linked_mesh().trans_of_convex(cv));
1685  for (auto&& val : bmin) val -= EPS;
1686  for (auto&& val : bmax) val += EPS;
1687  boxtree.add_box(bmin, bmax, cv);
1688  }
1689  boxtree.build_tree();
1690  }
1691 
1692  bool interpolator_on_mesh_fem::find_a_point(const base_node &pt,
1693  base_node &ptr,
1694  size_type &cv) const {
1695  bool gt_invertible;
1696  if (cv_stored != size_type(-1) && gic.invert(pt, ptr, gt_invertible)) {
1697  cv = cv_stored;
1698  if (gt_invertible)
1699  return true;
1700  }
1701 
1702  boxtree.find_boxes_at_point(pt, boxlst);
1703  for (const auto &box : boxlst) {
1704  gic = bgeot::geotrans_inv_convex
1705  (mf.linked_mesh().convex(box->id),
1706  mf.linked_mesh().trans_of_convex(box->id));
1707  cv_stored = box->id;
1708  if (gic.invert(pt, ptr, gt_invertible)) {
1709  cv = box->id;
1710  return true;
1711  }
1712  }
1713  return false;
1714  }
1715 
1716  bool interpolator_on_mesh_fem::eval(const base_node &pt, base_vector &val,
1717  base_matrix &grad) const {
1718  base_node ptref;
1719  size_type cv;
1720  base_vector coeff;
1721  dim_type q = mf.get_qdim(), N = mf.linked_mesh().dim();
1722  if (find_a_point(pt, ptref, cv)) {
1723  pfem pf = mf.fem_of_element(cv);
1724  bgeot::pgeometric_trans pgt =
1725  mf.linked_mesh().trans_of_convex(cv);
1726  base_matrix G;
1727  vectors_to_base_matrix(G, mf.linked_mesh().points_of_convex(cv));
1728  fem_interpolation_context fic(pgt, pf, ptref, G, cv, short_type(-1));
1729  slice_vector_on_basic_dof_of_element(mf, U, cv, coeff);
1730  // coeff.resize(mf.nb_basic_dof_of_element(cv));
1731  // gmm::copy(gmm::sub_vector
1732  // (U,gmm::sub_index(mf.ind_basic_dof_of_element(cv))), coeff);
1733  val.resize(q);
1734  pf->interpolation(fic, coeff, val, q);
1735  grad.resize(q, N);
1736  pf->interpolation_grad(fic, coeff, grad, q);
1737  return true;
1738  } else
1739  return false;
1740  }
1741 }
1742 
1743 /* end of namespace getfem */
bgeot::geotrans_inv_convex
does the inversion of the geometric transformation for a given convex
Definition: bgeot_geotrans_inv.h:63
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:250
gmm::copy
void copy(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:976
gmm::fill
void fill(L &l, typename gmm::linalg_traits< L >::value_type x)
*‍/
Definition: gmm_blas.h:103
gmm::clear
void clear(L &l)
clear (fill with zeros) a vector or matrix.
Definition: gmm_blas.h:58
gmm::resize
void resize(V &v, size_type n)
*‍/
Definition: gmm_blas.h:209
gmm::mult
void mult(const L1 &l1, const L2 &l2, L3 &l3)
*‍/
Definition: gmm_blas.h:1663
gmm::add
void add(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:1275
getfem::pfem
std::shared_ptr< const getfem::virtual_fem > pfem
type of pointer on a fem description
Definition: getfem_fem.h:243
bgeot::short_type
gmm::uint16_type short_type
used as the common short type integer in the library
Definition: bgeot_config.h:72
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:48
bgeot::pgeometric_trans
std::shared_ptr< const bgeot::geometric_trans > pgeometric_trans
pointer type for a geometric transformation
Definition: bgeot_geometric_trans.h:185
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:229
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:46
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:662