Rheolef  7.2
an efficient C++ finite element environment
form_lazy_terminal.h
Go to the documentation of this file.
1 # ifndef _RHEOLEF_FORM_LAZY_TERMINAL_H
2 # define _RHEOLEF_FORM_LAZY_TERMINAL_H
3 //
4 // This file is part of Rheolef.
5 //
6 // Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
7 //
8 // Rheolef is free software; you can redistribute it and/or modify
9 // it under the terms of the GNU General Public License as published by
10 // the Free Software Foundation; either version 2 of the License, or
11 // (at your option) any later version.
12 //
13 // Rheolef is distributed in the hope that it will be useful,
14 // but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 // GNU General Public License for more details.
17 //
18 // You should have received a copy of the GNU General Public License
19 // along with Rheolef; if not, write to the Free Software
20 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
21 //
22 // =========================================================================
23 // form_lazy_expr = result of a lazy integrate() that are combined in exprs
24 // AUTHOR: Pierre.Saramito@imag.fr
25 // DATE: 30 march 1920
26 //
27 // SUMMARY:
28 // 1. concept & base class
29 // 2. terminal
30 // 2.1. integrate
31 // 2.2. integrate on a band
32 // see also "form_lazy_expr.h"
33 
34 #include "rheolef/form.h"
35 #include "rheolef/field_lazy_form_mult.h"
36 
37 namespace rheolef {
38 
39 // -------------------------------------------------------------------
40 // 1. concept & base class
41 // -------------------------------------------------------------------
42 namespace details {
43 // Define a trait type for detecting form expression valid arguments
44 // template<class T> struct is_form_lazy: std::false_type {};
45 // => should be defined in form.h for compilation reason,
46 // otherwise Sfinae is always failed in form.h
47 
48 // Define a base class
49 template<class Derived>
51 public:
52 // a.trans_mult(u)
53 
54  template<class FieldExpr>
55  typename std::enable_if<
58  >::type
59  trans_mult (const FieldExpr& u_expr) const {
61  }
62  // TODO: how to obtain T and M from Derived ?
63  // typename Derived::scalar_type => compilation failed
64  // see eg EigenBase in eigen library
65  template<class T, class M>
67  trans_mult (const field_basic<T,M>& uh) const {
70  }
71 
72 protected:
73  Derived& derived() { return *static_cast< Derived*>(this); }
74  const Derived& derived() const { return *static_cast<const Derived*>(this); }
75 };
76 
77 } // namespace details
78 // -------------------------------------------------------------------
79 // 2. terminal
80 // -------------------------------------------------------------------
81 // 2.1. integrate
82 // -------------------------------------------------------------------
83 // form_lazy_terminal_integrate is returned by the integrate() function
84 // It contains the domain of integration
85 // and integration options (quadrature)
86 namespace details {
87 
88 template<class Expr>
90 public :
91 // definitions:
92 
94  using memory_type = typename Expr::memory_type;
95  using scalar_type = typename Expr::scalar_type;
96  using space_type = typename Expr::space_type;
97  using geo_type = typename Expr::geo_type;
100  using matrix_element_type = Eigen::Matrix<scalar_type,Eigen::Dynamic,Eigen::Dynamic>;
101 
102 // allocator:
103 
105  const geo_type& domain,
106  const Expr& expr,
107  const integrate_option& iopt)
108  : _domain(domain),
109  _expr(expr),
110  _iopt(iopt),
111  _prev_omega_K(),
112  _prev_K_dis_ie(std::numeric_limits<size_type>::max()),
113  _prev_ak()
114  {}
115 
116 // accessors:
117 
118  const geo_type& get_geo() const { return _domain; }
119  const space_type& get_trial_space() const { return _expr.get_trial_space(); }
120  const space_type& get_test_space() const { return _expr.get_test_space(); }
121  bool is_on_band() const { return false; }
122  band_type get_band() const { return band_type(); }
123 
124  void initialize (const geo_type& omega_K) const;
125  bool is_diagonal() const { return false; }
126 
127  void evaluate (
128  const geo_type& omega_K,
129  const geo_element& K,
130  matrix_element_type& ak) const;
131 // data:
132 protected:
134  mutable Expr _expr;
139 };
140 // inlined;
141 template<class Expr>
142 void
144 {
145  // TODO: sub-domain
146  // _domain could be a subdomain of omega_K:
147  // in that case evaluate(omega_K,K) should return a zero ak matrix
148  // when K do not belongs to _domain
149  _iopt._is_on_interface = false;
150  _iopt._is_inside_on_local_sides = false;
151  _expr.initialize (_domain, _iopt);
152  _prev_omega_K = omega_K;
153  _prev_K_dis_ie = std::numeric_limits<size_type>::max();
154 }
155 template<class Expr>
156 void
158  const geo_type& omega_K,
159  const geo_element& K,
160  matrix_element_type& ak) const
161 {
162  if (_prev_omega_K == omega_K && _prev_K_dis_ie == K.dis_ie()) {
163  trace_macro("interpolate(K="<<K.name()<<K.dis_ie()<<",prev="<<_prev_K_dis_ie<<"): re-use");
164  ak = _prev_ak;
165  return;
166  }
167  _expr.evaluate (omega_K, K, ak);
168  _prev_ak = ak; // expensive to compute, so memorize it for common subexpressions
169  _prev_omega_K = omega_K;
170  _prev_K_dis_ie = K.dis_ie();
171  trace_macro("interpolate(K="<<K.name()<<K.dis_ie()<<",prev="<<_prev_K_dis_ie<<"): compute");
172 }
173 
174 template<class Expr>
175 class form_lazy_terminal_integrate: public smart_pointer_nocopy<form_lazy_terminal_integrate_rep<Expr>>,
176  public form_lazy_base <form_lazy_terminal_integrate<Expr>> {
177 public :
178 // definitions:
179 
183  using size_type = typename rep::size_type;
184  using memory_type = typename rep::memory_type;
185  using scalar_type = typename rep::scalar_type;
186  using space_type = typename rep::space_type;
187  using geo_type = typename rep::geo_type;
188  using band_type = typename rep::band_type;
190 
191 // allocator:
192 
194  const geo_type& domain,
195  const Expr& expr,
196  const integrate_option& iopt)
197  : base1(new_macro(rep(domain,expr,iopt))),
198  base2()
199  {}
200 
201 // accessors:
202 
203  const geo_type& get_geo() const { return base1::data().get_geo(); }
204  const space_type& get_trial_space() const { return base1::data().get_trial_space(); }
205  const space_type& get_test_space() const { return base1::data().get_test_space(); }
206  bool is_on_band() const { return base1::data().is_on_band(); }
207  band_type get_band() const { return base1::data().get_band(); }
208 
209  void initialize (const geo_type& omega_K) const { return base1::data().initialize (omega_K); }
210  bool is_diagonal() const { return base1::data().is_diagonal(); }
211 
212  void evaluate (
213  const geo_type& omega_K,
214  const geo_element& K,
215  matrix_element_type& ak) const
216  { base1::data().evaluate (omega_K, K, ak); }
217 };
218 // concept;
219 template<class Expr> struct is_form_lazy <form_lazy_terminal_integrate<Expr> > : std::true_type {};
220 
221 }// namespace details
222 // ----------------------------------------------
223 // 2.2. integrate on a band
224 // ----------------------------------------------
225 namespace details {
226 
227 template<class Expr>
229 public :
230 // definitions:
231 
233  using memory_type = typename Expr::memory_type;
234  using scalar_type = typename Expr::scalar_type;
235  using space_type = typename Expr::space_type;
236  using geo_type = typename Expr::geo_type;
239  using matrix_element_type = Eigen::Matrix<scalar_type,Eigen::Dynamic,Eigen::Dynamic>;
240 
241 // allocators:
242 
244  const band_type& gh,
245  const Expr& expr,
246  const integrate_option& iopt);
247 
248 // accessors:
249 
250  const geo_type& get_geo() const { return _gh.level_set(); }
251  const space_type& get_trial_space() const { return _expr.get_trial_space(); }
252  const space_type& get_test_space() const { return _expr.get_test_space(); }
253  bool is_on_band() const { return true; }
254  band_type get_band() const { return _gh; }
255 
256  void initialize (const geo_type& omega_K) const;
257 
258  void evaluate (
259  const geo_type& omega_K,
260  const geo_element& K,
261  matrix_element_type& ak) const;
262 // data:
263 protected:
265  mutable Expr _expr;
270 };
271 // inlined;
272 template<class Expr>
274  const band_type& gh,
275  const Expr& expr,
276  const integrate_option& iopt)
277 : _gh(gh),
278  _expr(expr),
279  _iopt(iopt),
280  _prev_omega_K(),
281  _prev_K_dis_ie(std::numeric_limits<size_type>::max()),
282  _prev_ak()
283 {
284 }
285 template<class Expr>
286 void
288 {
289  const geo_type& band = _gh.band(); // the 3D intersected tetras
290  const geo_type& bgd_omega = get_test_space().get_constitution().get_background_geo();
291  check_macro (omega_K == get_geo(),
292  "integrate on band: unexpected integration domain");
293  check_macro (band.get_background_geo() == bgd_omega,
294  "do_integrate: incompatible integration domain "<<_gh.level_set().name() << " and test function based domain "
295  << bgd_omega.name());
296  _expr.initialize (_gh, _iopt);
297  _prev_omega_K = omega_K;
298  _prev_K_dis_ie = std::numeric_limits<size_type>::max();
299 }
300 template<class Expr>
301 void
303  const geo_type& omega_K,
304  const geo_element& K,
305  matrix_element_type& ak) const
306 {
307  if (_prev_omega_K == omega_K && _prev_K_dis_ie == K.dis_ie()) {
308  trace_macro("interpolate(K="<<K.name()<<K.dis_ie()<<",prev="<<_prev_K_dis_ie<<"): re-use");
309  ak = _prev_ak;
310  return;
311  }
312  // TODO: memorize _prev
313  _expr.evaluate (omega_K, K, ak);
314  _prev_ak = ak; // expensive to compute, so memorize it for common subexpressions
315  _prev_omega_K = omega_K;
316  _prev_K_dis_ie = K.dis_ie();
317 }
318 template<class Expr>
319 class form_lazy_terminal_integrate_band: public smart_pointer_nocopy<form_lazy_terminal_integrate_band_rep<Expr>>,
320  public form_lazy_base <form_lazy_terminal_integrate_band<Expr>> {
321 public :
322 // definitions:
323 
327  using size_type = typename rep::size_type;
328  using memory_type = typename rep::memory_type;
329  using scalar_type = typename rep::scalar_type;
330  using space_type = typename rep::space_type;
331  using geo_type = typename rep::geo_type;
332  using band_type = typename rep::band_type;
334 
335 // allocator:
336 
338  const band_type& gh,
339  const Expr& expr,
340  const integrate_option& iopt)
341  : base1(new_macro(rep(gh,expr,iopt))),
342  base2()
343  {}
344 
345 // accessors:
346 
347  const geo_type& get_geo() const { return base1::data().get_geo(); }
348  const space_type& get_trial_space() const { return base1::data().get_trial_space(); }
349  const space_type& get_test_space() const { return base1::data().get_test_space(); }
350  bool is_on_band() const { return base1::data().is_on_band(); }
351  band_type get_band() const { return base1::data().get_band(); }
352 
353  void initialize (const geo_type& omega_K) const { return base1::data().initialize (omega_K); }
354 
355  void evaluate (
356  const geo_type& omega_K,
357  const geo_element& K,
358  matrix_element_type& ak) const
359  { base1::data().evaluate (omega_K, K, ak); }
360 };
361 // concept;
362 template<class Expr> struct is_form_lazy <form_lazy_terminal_integrate_band<Expr> > : std::true_type {};
363 
364 }// namespace details
365 
366 template <class Expr>
367 typename
368 std::enable_if<
371 >::type
374  const band_basic<typename float_traits<typename Expr::scalar_type>::type,
375  typename Expr::memory_type>& gh,
376  const Expr& expr,
377  const integrate_option& iopt = integrate_option())
378 {
380 }
381 template <class Expr>
382 typename
383 std::enable_if<
384  details::is_form_expr_v2_variational_arg<Expr>::value
385  ,details::form_lazy_terminal_integrate_band <details::form_expr_quadrature_on_element<Expr>>
386 >::type
389  const band_basic<typename float_traits<typename Expr::scalar_type>::type,
390  typename Expr::memory_type>& gh,
391  const Expr& expr,
392  const integrate_option& iopt = integrate_option())
393 {
394 
396  return lazy_integrate (gh, expr_quad, iopt);
397 }
398 
399 }// namespace rheolef
400 # endif /* _RHEOLEF_FORM_LAZY_TERMINAL_H */
field gh(Float epsilon, Float t, const field &uh, const test &v)
see the band page for the full documentation
const geo_basic< T, M > & level_set() const
Definition: band.h:109
field_lazy_trans_mult_form< Derived, field_lazy_terminal_field< T, M > > trans_mult(const field_basic< T, M > &uh) const
const Derived & derived() const
std::enable_if< details::is_field_lazy< FieldExpr >::value,field_lazy_trans_mult_form< Derived, FieldExpr > >::type trans_mult(const FieldExpr &u_expr) const
void evaluate(const geo_type &omega_K, const geo_element &K, matrix_element_type &ak) const
Eigen::Matrix< scalar_type, Eigen::Dynamic, Eigen::Dynamic > matrix_element_type
form_lazy_terminal_integrate_band_rep(const band_type &gh, const Expr &expr, const integrate_option &iopt)
typename float_traits< scalar_type >::type float_type
void evaluate(const geo_type &omega_K, const geo_element &K, matrix_element_type &ak) const
void initialize(const geo_type &omega_K) const
form_lazy_terminal_integrate_band(const band_type &gh, const Expr &expr, const integrate_option &iopt)
typename rep::matrix_element_type matrix_element_type
band_basic< float_type, memory_type > band_type
void evaluate(const geo_type &omega_K, const geo_element &K, matrix_element_type &ak) const
void initialize(const geo_type &omega_K) const
form_lazy_terminal_integrate_rep(const geo_type &domain, const Expr &expr, const integrate_option &iopt)
Eigen::Matrix< scalar_type, Eigen::Dynamic, Eigen::Dynamic > matrix_element_type
typename float_traits< scalar_type >::type float_type
void evaluate(const geo_type &omega_K, const geo_element &K, matrix_element_type &ak) const
void initialize(const geo_type &omega_K) const
form_lazy_terminal_integrate(const geo_type &domain, const Expr &expr, const integrate_option &iopt)
typename rep::matrix_element_type matrix_element_type
see the geo_element page for the full documentation
Definition: geo_element.h:102
reference_element::size_type size_type
Definition: geo_element.h:125
size_type dis_ie() const
Definition: geo_element.h:163
char name() const
Definition: geo_element.h:169
see the integrate_option page for the full documentation
rheolef::std type
#define trace_macro(message)
Definition: dis_macros.h:111
void get_geo(istream &in, my_geo &omega)
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
This file is part of Rheolef.
std::enable_if< details::is_field_expr_quadrature_arg< Expr >::value,details::field_lazy_terminal_integrate< Expr >>::type lazy_integrate(const typename Expr::geo_type &domain, const Expr &expr, const integrate_option &iopt=integrate_option())
see the integrate page for the full documentation