GetFEM  5.5
getfem_global_function.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2004-2026 Yves Renard
5  Copyright (C) 2016-2020 Konstantinos Poulios
6 
7  This file is a part of GetFEM
8 
9  GetFEM is free software; you can redistribute it and/or modify it
10  under the terms of the GNU Lesser General Public License as published
11  by the Free Software Foundation; either version 3 of the License, or
12  (at your option) any later version along with the GCC Runtime Library
13  Exception either version 3.1 or (at your option) any later version.
14  This program is distributed in the hope that it will be useful, but
15  WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
16  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
17  License and GCC Runtime Library Exception for more details.
18  You should have received a copy of the GNU Lesser General Public License
19  along with this program. If not, see https://www.gnu.org/licenses/.
20 
21  As a special exception, you may use this file as it is a part of a free
22  software library without restriction. Specifically, if other files
23  instantiate templates or use macros or inline functions from this file,
24  or you compile this file and link it with other files to produce an
25  executable, this file does not by itself cause the resulting executable
26  to be covered by the GNU Lesser General Public License. This exception
27  does not however invalidate any other reasons why the executable file
28  might be covered by the GNU Lesser General Public License.
29 
30 ===========================================================================*/
31 
32 /**@file getfem_global_function.h
33  @author Yves Renard <Yves.Renard@insa-lyon.fr>, J. Pommier
34  @date March, 2005.
35  @brief Definition of global functions to be used as base or enrichment
36  functions in fem.
37 */
38 #ifndef GETFEM_GLOBAL_FUNCTION_H__
39 #define GETFEM_GLOBAL_FUNCTION_H__
40 
41 #include "bgeot_rtree.h"
42 #include "getfem_mesh_fem.h"
43 #include "getfem_generic_assembly.h"
44 
45 namespace getfem {
46 
47  /// inherit from this class to define new global functions.
48  class global_function : virtual public dal::static_stored_object {
49  protected:
50  const dim_type dim_;
51  public:
52  dim_type dim() const { return dim_; }
53 
54  virtual scalar_type val(const fem_interpolation_context&) const
55  { GMM_ASSERT1(false, "this global_function has no value"); }
56  virtual void grad(const fem_interpolation_context&, base_small_vector&) const
57  { GMM_ASSERT1(false, "this global_function has no gradient"); }
58  virtual void hess(const fem_interpolation_context&, base_matrix&) const
59  { GMM_ASSERT1(false, "this global_function has no hessian"); }
60 
61  virtual bool is_in_support(const base_node & /* pt */ ) const { return true; }
62  virtual void bounding_box(base_node &bmin, base_node &bmax) const {
63  GMM_ASSERT1(bmin.size() == dim_ && bmax.size() == dim_,
64  "Wrong dimensions");
65  for (auto&& xx : bmin) xx = -1e+25;
66  for (auto&& xx : bmax) xx = 1e+25;
67  }
68 
69  global_function(dim_type dim__) : dim_(dim__)
70  { DAL_STORED_OBJECT_DEBUG_CREATED(this, "Global function");}
71  virtual ~global_function()
72  { DAL_STORED_OBJECT_DEBUG_DESTROYED(this, "Global function"); }
73  };
74 
75  typedef std::shared_ptr<const global_function> pglobal_function;
76 
77 
78  class global_function_simple : public global_function {
79  public:
80  // These virtual functions can not be further overriden in derived classes
81  virtual scalar_type val(const fem_interpolation_context&) const final;
82  virtual void grad(const fem_interpolation_context&, base_small_vector&) const final;
83  virtual void hess(const fem_interpolation_context&, base_matrix&) const final;
84  // You have to override these instead
85  virtual scalar_type val(const base_node &pt) const = 0;
86  virtual void grad(const base_node &pt, base_small_vector&) const = 0;
87  virtual void hess(const base_node &pt, base_matrix&) const = 0;
88 
89  global_function_simple(dim_type dim__) : global_function(dim__)
90  { DAL_STORED_OBJECT_DEBUG_CREATED(this, "Global function simple"); }
91  virtual ~global_function_simple()
92  { DAL_STORED_OBJECT_DEBUG_DESTROYED(this, "Global function simple"); }
93  };
94 
95  class global_function_parser : public global_function_simple {
96  ga_workspace gw;
97  ga_function f_val, f_grad, f_hess;
98  mutable model_real_plain_vector pt_;
99  public:
100  virtual scalar_type val(const base_node &pt) const;
101  virtual const base_tensor &tensor_val(const base_node &pt) const;
102  virtual void grad(const base_node &pt, base_small_vector &g) const;
103  virtual void hess(const base_node &pt, base_matrix &h) const;
104 
105  global_function_parser(dim_type dim_,
106  const std::string &sval,
107  const std::string &sgrad="",
108  const std::string &shess="");
109  virtual ~global_function_parser()
110  { DAL_STORED_OBJECT_DEBUG_DESTROYED(this, "Global function parser"); }
111  };
112 
113 
114  class global_function_sum : public global_function {
115  std::vector<pglobal_function> functions;
116  public:
117  virtual scalar_type val(const fem_interpolation_context&) const;
118  virtual void grad(const fem_interpolation_context&, base_small_vector&) const;
119  virtual void hess(const fem_interpolation_context&, base_matrix&) const;
120  virtual bool is_in_support(const base_node &p) const;
121  virtual void bounding_box(base_node &bmin_, base_node &bmax_) const;
122  global_function_sum(const std::vector<pglobal_function> &funcs);
123  global_function_sum(pglobal_function f1, pglobal_function f2);
124  global_function_sum(pglobal_function f1, pglobal_function f2,
125  pglobal_function f3);
126  global_function_sum(pglobal_function f1, pglobal_function f2,
127  pglobal_function f3, pglobal_function f4);
128  virtual ~global_function_sum()
129  { DAL_STORED_OBJECT_DEBUG_DESTROYED(this, "Global function sum"); }
130  };
131 
132  class global_function_product : public global_function {
133  pglobal_function f1, f2;
134  public:
135  virtual scalar_type val(const fem_interpolation_context&) const;
136  virtual void grad(const fem_interpolation_context&, base_small_vector&) const;
137  virtual void hess(const fem_interpolation_context&, base_matrix&) const;
138  virtual bool is_in_support(const base_node &p) const;
139  virtual void bounding_box(base_node &bmin_, base_node &bmax_) const;
140  global_function_product(pglobal_function f1_, pglobal_function f2_);
141  virtual ~global_function_product()
142  { DAL_STORED_OBJECT_DEBUG_DESTROYED(this, "Global function product"); }
143  };
144 
145  class global_function_bounded : public global_function {
146  pglobal_function f;
147  const base_node bmin, bmax;
148  bool has_expr;
149  ga_workspace gw;
150  ga_function f_val;
151  mutable model_real_plain_vector pt_;
152  public:
153  virtual scalar_type val(const fem_interpolation_context &c) const
154  { return f->val(c); }
155  virtual void grad(const fem_interpolation_context &c, base_small_vector &g) const
156  { f->grad(c, g); }
157  virtual void hess(const fem_interpolation_context &c, base_matrix &h) const
158  { f->hess(c, h); }
159 
160  virtual bool is_in_support(const base_node &) const;
161  virtual void bounding_box(base_node &bmin_, base_node &bmax_) const {
162  bmin_ = bmin;
163  bmax_ = bmax;
164  }
165  // is_in_expr should evaluate to a positive result inside the support of the
166  // function and negative elsewhere
167  global_function_bounded(pglobal_function f_,
168  const base_node &bmin_, const base_node &bmax_,
169  const std::string &is_in_expr="");
170  virtual ~global_function_bounded()
171  { DAL_STORED_OBJECT_DEBUG_DESTROYED(this, "Global function bounded"); }
172  };
173 
174 
175  /** a general structure for interpolation of a function defined
176  by a mesh_fem and a vector U at any point
177  (interpolation of value and gradient).
178  */
179 
180  struct interpolator_on_mesh_fem {
181  const mesh_fem &mf;
182  std::vector<scalar_type> U;
183 
184  mutable bgeot::rtree boxtree;
185  mutable size_type cv_stored;
186  mutable bgeot::rtree::pbox_set boxlst;
187  mutable bgeot::geotrans_inv_convex gic;
188 
189 
190  interpolator_on_mesh_fem(const mesh_fem &mf_,
191  const std::vector<scalar_type> &U_);
192  bool find_a_point(const base_node &pt, base_node &ptr,
193  size_type &cv) const;
194  bool eval(const base_node &pt, base_vector &val, base_matrix &grad) const;
195  };
196 
197  typedef std::shared_ptr<const interpolator_on_mesh_fem>
198  pinterpolator_on_mesh_fem;
199 
200 
201  /** below a list of simple functions of (x,y)
202  used for building the crack singular functions
203  */
204  class abstract_xy_function : virtual public dal::static_stored_object {
205  public:
206  virtual scalar_type val(scalar_type x, scalar_type y) const = 0;
207  virtual base_small_vector grad(scalar_type x, scalar_type y) const = 0;
208  virtual base_matrix hess(scalar_type x, scalar_type y) const = 0;
209  virtual ~abstract_xy_function() {}
210  };
211 
212  typedef std::shared_ptr<const abstract_xy_function> pxy_function;
213 
214  struct parser_xy_function : public abstract_xy_function {
215  ga_workspace gw;
216  ga_function f_val, f_grad, f_hess;
217  mutable model_real_plain_vector ptx, pty, ptr, ptt;
218 
219  virtual scalar_type val(scalar_type x, scalar_type y) const;
220  virtual base_small_vector grad(scalar_type x, scalar_type y) const;
221  virtual base_matrix hess(scalar_type x, scalar_type y) const;
222 
223  parser_xy_function(const std::string &sval,
224  const std::string &sgrad="0;0;",
225  const std::string &shess="0;0;0;0;");
226  virtual ~parser_xy_function() {}
227  };
228 
229  struct crack_singular_xy_function : public abstract_xy_function {
230  unsigned l; /* 0 <= l <= 6 */
231  virtual scalar_type val(scalar_type x, scalar_type y) const;
232  virtual base_small_vector grad(scalar_type x, scalar_type y) const;
233  virtual base_matrix hess(scalar_type x, scalar_type y) const;
234  crack_singular_xy_function(unsigned l_) : l(l_) {}
235  virtual ~crack_singular_xy_function() {}
236  };
237 
238  struct cutoff_xy_function : public abstract_xy_function {
239  enum { NOCUTOFF = -1,
240  EXPONENTIAL_CUTOFF = 0,
241  POLYNOMIAL_CUTOFF = 1,
242  POLYNOMIAL2_CUTOFF= 2 };
243  int fun;
244  scalar_type a4, r1, r0;
245  virtual scalar_type val(scalar_type x, scalar_type y) const;
246  virtual base_small_vector grad(scalar_type x, scalar_type y) const;
247  virtual base_matrix hess(scalar_type x, scalar_type y) const;
248  cutoff_xy_function(int fun_num, scalar_type r,
249  scalar_type r1, scalar_type r0);
250  virtual ~cutoff_xy_function() {}
251  };
252 
253  struct interpolated_xy_function : public abstract_xy_function {
254  pinterpolator_on_mesh_fem itp;
255  size_type component;
256  virtual scalar_type val(scalar_type x, scalar_type y) const {
257  base_vector v; base_matrix g;
258  itp->eval(base_node(x,y), v, g);
259  return v[component];
260  }
261  virtual base_small_vector grad(scalar_type x, scalar_type y) const {
262  base_vector v; base_matrix g;
263  itp->eval(base_node(x,y), v, g);
264  return base_small_vector(g(component,0), g(component,1));
265  }
266  virtual base_matrix hess(scalar_type, scalar_type) const
267  { GMM_ASSERT1(false, "Sorry, to be done ..."); }
268  interpolated_xy_function(const pinterpolator_on_mesh_fem &itp_,
269  size_type c) :
270  itp(itp_), component(c) {}
271  virtual ~interpolated_xy_function() {}
272  };
273 
274  struct product_of_xy_functions : public abstract_xy_function {
275  pxy_function fn1, fn2;
276  scalar_type val(scalar_type x, scalar_type y) const {
277  return fn1->val(x,y) * fn2->val(x,y);
278  }
279  base_small_vector grad(scalar_type x, scalar_type y) const {
280  return fn1->grad(x,y)*fn2->val(x,y) + fn1->val(x,y)*fn2->grad(x,y);
281  }
282  virtual base_matrix hess(scalar_type x, scalar_type y) const {
283  base_matrix h = fn1->hess(x,y);
284  gmm::scale(h, fn2->val(x,y));
285  gmm::add(gmm::scaled(fn2->hess(x,y), fn1->val(x,y)), h);
286  gmm::rank_two_update(h, fn1->grad(x,y), fn2->grad(x,y));
287  return h;
288  }
289  product_of_xy_functions(pxy_function &fn1_, pxy_function &fn2_)
290  : fn1(fn1_), fn2(fn2_) {}
291  virtual ~product_of_xy_functions() {}
292  };
293 
294  struct add_of_xy_functions : public abstract_xy_function {
295  pxy_function fn1, fn2;
296  scalar_type val(scalar_type x, scalar_type y) const {
297  return fn1->val(x,y) + fn2->val(x,y);
298  }
299  base_small_vector grad(scalar_type x, scalar_type y) const {
300  return fn1->grad(x,y) + fn2->grad(x,y);
301  }
302  virtual base_matrix hess(scalar_type x, scalar_type y) const {
303  base_matrix h = fn1->hess(x,y);
304  gmm::add(fn2->hess(x,y), h);
305  return h;
306  }
307  add_of_xy_functions(const pxy_function &fn1_, const pxy_function &fn2_)
308  : fn1(fn1_), fn2(fn2_) {}
309  virtual ~add_of_xy_functions() {}
310  };
311 
312 
313  /** some useful global functions */
314  class level_set;
315 
316  pglobal_function
317  global_function_on_level_set(const level_set &ls, const pxy_function &fn);
318 
319  pglobal_function
320  global_function_on_level_sets(const std::vector<level_set> &lsets,
321  const pxy_function &fn);
322 
323  pglobal_function
324  global_function_bspline(const scalar_type xmin, const scalar_type xmax,
325  const size_type order, const size_type xtype);
326 
327  pglobal_function
328  global_function_bspline(const scalar_type xmin, const scalar_type xmax,
329  const scalar_type ymin, const scalar_type ymax,
330  const size_type order,
331  const size_type xtype, const size_type ytype);
332 
333  pglobal_function
334  global_function_bspline(const scalar_type xmin, const scalar_type xmax,
335  const scalar_type ymin, const scalar_type ymax,
336  const scalar_type zmin, const scalar_type zmax,
337  const size_type order,
338  const size_type xtype, const size_type ytype,
339  const size_type ztype);
340 } /* end of namespace getfem. */
341 
342 #endif
bgeot_rtree.h
region-tree for window/point search on a set of rectangles.
bgeot::geotrans_inv_convex
does the inversion of the geometric transformation for a given convex
Definition: bgeot_geotrans_inv.h:63
bgeot::rtree
Balanced tree of n-dimensional rectangles.
Definition: bgeot_rtree.h:97
dal::static_stored_object
base class for static stored objects
Definition: dal_static_stored_objects.h:206
getfem::abstract_xy_function
below a list of simple functions of (x,y) used for building the crack singular functions
Definition: getfem_global_function.h:204
getfem::fem_interpolation_context
structure passed as the argument of fem interpolation functions.
Definition: getfem_fem.h:749
getfem::global_function
inherit from this class to define new global functions.
Definition: getfem_global_function.h:48
getfem::mesh_fem
Describe a finite element method linked to a mesh.
Definition: getfem_mesh_fem.h:151
getfem_generic_assembly.h
A language for generic assembly of pde boundary value problems.
getfem_mesh_fem.h
Define the getfem::mesh_fem class.
gmm::add
void add(const L1 &l1, L2 &l2)
*‍/
Definition: gmm_blas.h:1275
bgeot::size_type
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:48
getfem
GEneric Tool for Finite Element Methods.
Definition: getfem_accumulated_distro.h:46
getfem::interpolator_on_mesh_fem
a general structure for interpolation of a function defined by a mesh_fem and a vector U at any point...
Definition: getfem_global_function.h:180