GetFEM++  5.3
getfem_global_function.h
Go to the documentation of this file.
1 /* -*- c++ -*- (enables emacs c++ mode) */
2 /*===========================================================================
3 
4  Copyright (C) 2004-2017 Yves Renard
5  Copyright (C) 2016 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 <Yves.Renard@insa-lyon.fr>, 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"
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 
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 
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  */
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 
324 } /* end of namespace getfem. */
325 
326 #endif
base class for static stored objects
inherit from this class to define new global functions.
a general structure for interpolation of a function defined by a mesh_fem and a vector U at any point...
Define the getfem::mesh_fem class.
does the inversion of the geometric transformation for a given convex
Define a level-set.
size_t size_type
used as the common size type in the library
Definition: bgeot_poly.h:49
structure passed as the argument of fem interpolation functions.
Definition: getfem_fem.h:741
below a list of simple functions of (x,y) used for building the crack singular functions ...
A langage for generic assembly of pde boundary value problems.
GEneric Tool for Finite Element Methods.
Describe a finite element method linked to a mesh.
Balanced tree of n-dimensional rectangles.
Definition: bgeot_rtree.h:66
region-tree for window/point search on a set of rectangles.