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