Rheolef  7.2
an efficient C++ finite element environment
field_lazy_form_mult.h
Go to the documentation of this file.
1 # ifndef _RHEOLEF_FIELD_LAZY_FORM_MULT_H
2 # define _RHEOLEF_FIELD_LAZY_FORM_MULT_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*field_lazy and such
24 // AUTHOR: Pierre.Saramito@imag.fr
25 // DATE: 6 april 1920
26 //
27 // SUMMARY:
28 // 1. form_lazy*field_lazy -> field_lazy
29 // 2. form_lazy.trans_mult (field_lazy) -> field_lazy
30 //
31 #include "rheolef/field_lazy_node.h"
32 #include "rheolef/form.h"
33 
34 namespace rheolef {
35 
36 // -------------------------------------------------------------------
37 // 1. form_lazy*field_lazy -> field_lazy
38 // -------------------------------------------------------------------
39 namespace details {
40 
41 template<class FormExpr, class FieldExpr>
43 public :
44 // definitions:
45 
47  using memory_type = typename FieldExpr::memory_type;
48  using scalar_type = typename FieldExpr::scalar_type;
49  using space_type = typename FieldExpr::space_type;
50  using geo_type = typename FieldExpr::geo_type;
53  using vector_element_type = Eigen::Matrix<scalar_type,Eigen::Dynamic,1>;
54  using matrix_element_type = Eigen::Matrix<scalar_type,Eigen::Dynamic,Eigen::Dynamic>;
55 
56 // allocator:
57 
58  field_lazy_mult_form_rep (const FormExpr& a_expr, const FieldExpr& u_expr);
59 
60 // accessors:
61 
62  const geo_type& get_geo() const { return _a_expr.get_geo(); }
63  const space_type& get_space() const { return _a_expr.get_test_space(); }
64  bool is_on_band() const { return _u_expr.is_on_band(); }
65  band_type get_band() const { return _u_expr.get_band(); }
66 
67  void initialize (const geo_type& omega_K) const;
68 
69  void evaluate (
70  const geo_type& omega_K,
71  const geo_element& K,
72  vector_element_type& uk) const;
73 // data:
74 protected:
75  mutable FormExpr _a_expr;
76  mutable FieldExpr _u_expr;
80 };
81 // inlined;
82 template<class FormExpr, class FieldExpr>
84  const FormExpr& a_expr,
85  const FieldExpr& u_expr)
86  : _a_expr(a_expr),
87  _u_expr(u_expr),
88  _prev_omega_K(),
89  _prev_K_dis_ie(std::numeric_limits<size_type>::max()),
90  _prev_vk()
91 {
92 }
93 template<class FormExpr, class FieldExpr>
94 void
96 {
97  // TODO: subdomains e.g. robin
98  check_macro (_a_expr.get_trial_space().get_geo() == _u_expr.get_geo(),
99  "lazy_multiply: different domain not yet supported");
100 
101  check_macro (_a_expr.get_trial_space() == _u_expr.get_space(),
102  "lazy_multiply: incompatible spaces \""
103  <<_a_expr.get_trial_space().name()<<"\" and \""
104  <<_u_expr.get_space().name()<<"\"");
105 
106  if (! _a_expr. get_test_space().get_constitution().have_compact_support_inside_element() ||
107  ! _a_expr.get_trial_space().get_constitution().have_compact_support_inside_element()) {
108  warning_macro("lazy_multiply: requires compact support inside elements (e.g. discontinuous or bubble)");
109  warning_macro("lazy_multiply: form spaces was "
110  << "[\"" <<_a_expr.get_trial_space().name()<<"\", \"" <<_a_expr.get_test_space().name()<<"\"]");
111  fatal_macro("lazy_multiply: HINT: convert to \"form\" before to do the product");
112  }
113  _a_expr.initialize (omega_K);
114  _u_expr.initialize (omega_K);
115  _prev_omega_K = omega_K;
116  _prev_K_dis_ie = std::numeric_limits<size_type>::max();
117 }
118 template<class FormExpr, class FieldExpr>
119 void
121  const geo_type& omega_K,
122  const geo_element& K,
123  vector_element_type& vk) const
124 {
125  if (_prev_omega_K == omega_K && _prev_K_dis_ie == K.dis_ie()) {
126  trace_macro("amulx(K="<<K.name()<<K.dis_ie()<<",prev="<<_prev_K_dis_ie<<"): re-use");
127  vk = _prev_vk;
128  return;
129  }
130  trace_macro("amulx(K="<<K.name()<<K.dis_ie()<<",prev="<<_prev_K_dis_ie<<"): compute");
133  _a_expr.evaluate (omega_K, K, ak);
134  _u_expr.evaluate (omega_K, K, uk);
135 #ifdef _RHEOLEF_PARANO
136  check_macro (ak.cols() == uk.size(), "a*u: invalid sizes");
137 #endif // _RHEOLEF_PARANO
138  vk = ak*uk;
139  _prev_omega_K = omega_K;
140  _prev_vk = vk; // expensive to compute, so memorize it for common subexpressions
141  _prev_K_dis_ie = K.dis_ie();
142 }
143 template<class FormExpr, class FieldExpr>
144 class field_lazy_mult_form: public smart_pointer_nocopy<field_lazy_mult_form_rep<FormExpr,FieldExpr>>,
145  public field_lazy_base <field_lazy_mult_form <FormExpr,FieldExpr>> {
146 public :
147 // definitions:
148 
152  using size_type = typename rep::size_type;
153  using memory_type = typename rep::memory_type;
154  using scalar_type = typename rep::scalar_type;
155  using space_type = typename rep::space_type;
156  using geo_type = typename rep::geo_type;
157  using band_type = typename rep::band_type;
159 
160 // allocator:
161 
162  field_lazy_mult_form (const FormExpr& a_expr, const FieldExpr& u_expr)
163  : base1(new_macro(rep(a_expr,u_expr))),
164  base2()
165  {}
166 
167 // accessors:
168 
169  const geo_type& get_geo() const { return base1::data().get_geo(); }
170  const space_type& get_space() const { return base1::data().get_space(); }
171  bool is_on_band() const { return base1::data().is_on_band(); }
172  band_type get_band() const { return base1::data().get_band(); }
173 
174  void initialize (const geo_type& omega_K) const { base1::data().initialize (omega_K); }
175 
176  void evaluate (
177  const geo_type& omega_K,
178  const geo_element& K,
179  vector_element_type& uk) const
180  { base1::data().evaluate (omega_K, K, uk); }
181 };
182 // concept;
183 template<class FormExpr, class FieldExpr>
184 struct is_field_lazy <field_lazy_mult_form<FormExpr,FieldExpr> > : std::true_type {};
185 
186 }// namespace details
187 
189 template<class FormExpr, class FieldExpr,
190  class Sfinae1 = typename std::enable_if<details:: is_form_lazy<FormExpr> ::value, FormExpr>::type,
191  class Sfinae2 = typename std::enable_if<details::is_field_lazy<FieldExpr>::value, FieldExpr>::type>
193 operator* (const FormExpr& a, const FieldExpr& u)
194 {
196 }
198 template<
199  class FormExpr
200  ,class Sfinae = typename std::enable_if<details::is_form_lazy<FormExpr>::value, FormExpr>::type
201 >
202 details::field_lazy_mult_form<
203  FormExpr
204  ,details::field_lazy_terminal_field<
205  typename FormExpr::scalar_type
206  ,typename FormExpr::memory_type
207  >
208 >
210  const FormExpr& a,
212 {
213  using wrap_type
215  typename FormExpr::scalar_type
216  ,typename FormExpr::memory_type>;
218 }
219 // -------------------------------------------------------------------
220 // 2. form_lazy.trans_mult (field_lazy) -> field_lazy
221 // -------------------------------------------------------------------
222 namespace details {
223 
224 template<class FormExpr, class FieldExpr>
226 public :
227 // definitions:
228 
230  using memory_type = typename FieldExpr::memory_type;
231  using scalar_type = typename FieldExpr::scalar_type;
232  using space_type = typename FieldExpr::space_type;
233  using geo_type = typename FieldExpr::geo_type;
236  using vector_element_type = Eigen::Matrix<scalar_type,Eigen::Dynamic,1>;
237  using matrix_element_type = Eigen::Matrix<scalar_type,Eigen::Dynamic,Eigen::Dynamic>;
238 
239 // allocator:
240 
241  field_lazy_trans_mult_form_rep (const FormExpr& a_expr, const FieldExpr& u_expr);
242 
243 // accessors:
244 
245  const geo_type& get_geo() const { return _a_expr.get_geo(); }
246  const space_type& get_space() const { return _a_expr.get_trial_space(); }
247  bool is_on_band() const { return _u_expr.is_on_band(); }
248  band_type get_band() const { return _u_expr.get_band(); }
249 
250  void initialize (const geo_type& omega_K) const;
251 
252  void evaluate (
253  const geo_type& omega_K,
254  const geo_element& K,
255  vector_element_type& uk) const;
256 // data:
257 protected:
258  mutable FormExpr _a_expr;
259  mutable FieldExpr _u_expr;
263 };
264 // inlined;
265 template<class FormExpr, class FieldExpr>
267  const FormExpr& a_expr,
268  const FieldExpr& u_expr)
269  : _a_expr(a_expr),
270  _u_expr(u_expr),
271  _prev_omega_K(),
272  _prev_K_dis_ie(std::numeric_limits<size_type>::max()),
273  _prev_vk()
274 {
275 }
276 template<class FormExpr, class FieldExpr>
277 void
279 {
280  // TODO: subdomains e.g. robin
281  check_macro (_a_expr.get_test_space().get_geo() == _u_expr.get_geo(),
282  "lazy_trans_mult: different domain not yet supported");
283 
284  check_macro (_a_expr.get_test_space() == _u_expr.get_space(),
285  "lazy_trans_mult: incompatible spaces \""
286  <<_a_expr.get_test_space().name()<<"\" and \""
287  <<_u_expr.get_space().name()<<"\"");
288 
289  if (! _a_expr. get_test_space().get_constitution().have_compact_support_inside_element() ||
290  ! _a_expr.get_trial_space().get_constitution().have_compact_support_inside_element()) {
291  warning_macro("lazy_trans_mult: requires compact support inside elements (e.g. discontinuous or bubble)");
292  warning_macro("lazy_trans_mult: form spaces was "
293  << "[\"" <<_a_expr.get_trial_space().name()<<"\", \"" <<_a_expr.get_test_space().name()<<"\"]");
294  fatal_macro("lazy_trans_mult: HINT: convert to \"form\" before to do the product");
295  }
296  _a_expr.initialize (omega_K);
297  _u_expr.initialize (omega_K);
298  _prev_omega_K = omega_K;
299  _prev_K_dis_ie = std::numeric_limits<size_type>::max();
300 }
301 template<class FormExpr, class FieldExpr>
302 void
304  const geo_type& omega_K,
305  const geo_element& K,
306  vector_element_type& vk) const
307 {
308  if (_prev_K_dis_ie == K.dis_ie()) {
309  trace_macro("atansx(K="<<K.name()<<K.dis_ie()<<",prev="<<_prev_K_dis_ie<<"): re-use");
310  vk = _prev_vk;
311  return;
312  }
313  trace_macro("atansx(K="<<K.name()<<K.dis_ie()<<",prev="<<_prev_K_dis_ie<<"): compute");
316  _a_expr.evaluate (omega_K, K, ak);
317  _u_expr.evaluate (omega_K, K, uk);
318 #ifdef _RHEOLEF_PARANO
319  check_macro (ak.rows() == uk.size(), "a*u: invalid sizes");
320 #endif // _RHEOLEF_PARANO
321  vk = uk.transpose()*ak;
322  _prev_vk = vk; // expensive to compute, so memorize it for common subexpressions
323  _prev_omega_K = omega_K;
324  _prev_K_dis_ie = K.dis_ie();
325 }
326 template<class FormExpr, class FieldExpr>
327 class field_lazy_trans_mult_form: public smart_pointer_nocopy<field_lazy_trans_mult_form_rep<FormExpr,FieldExpr>>,
328  public field_lazy_base <field_lazy_trans_mult_form <FormExpr,FieldExpr>> {
329 public :
330 // definitions:
331 
335  using size_type = typename rep::size_type;
336  using memory_type = typename rep::memory_type;
337  using scalar_type = typename rep::scalar_type;
338  using space_type = typename rep::space_type;
339  using geo_type = typename rep::geo_type;
340  using band_type = typename rep::band_type;
342 
343 // allocator:
344 
345  field_lazy_trans_mult_form (const FormExpr& a_expr, const FieldExpr& u_expr)
346  : base1(new_macro(rep(a_expr,u_expr))),
347  base2()
348  {}
349 
350 // accessors:
351 
352  const geo_type& get_geo() const { return base1::data().get_geo(); }
353  const space_type& get_space() const { return base1::data().get_space(); }
354  bool is_on_band() const { return base1::data().is_on_band(); }
355  band_type get_band() const { return base1::data().get_band(); }
356 
357  void initialize (const geo_type& omega_K) const { base1::data().initialize (omega_K); }
358 
359  void evaluate (
360  const geo_type& omega_K,
361  const geo_element& K,
362  vector_element_type& uk) const
363  { base1::data().evaluate (omega_K, K, uk); }
364 };
365 // concept;
366 template<class FormExpr, class FieldExpr>
367 struct is_field_lazy <field_lazy_trans_mult_form<FormExpr,FieldExpr> > : std::true_type {};
368 
369 }// namespace details
370 }// namespace rheolef
371 # endif /* _RHEOLEF_FIELD_LAZY_FORM_MULT_H */
band_basic< float_type, memory_type > band_type
typename FieldExpr::space_type space_type
typename FieldExpr::scalar_type scalar_type
void initialize(const geo_type &omega_K) const
field_lazy_mult_form_rep(const FormExpr &a_expr, const FieldExpr &u_expr)
void evaluate(const geo_type &omega_K, const geo_element &K, vector_element_type &uk) const
typename FieldExpr::memory_type memory_type
Eigen::Matrix< scalar_type, Eigen::Dynamic, Eigen::Dynamic > matrix_element_type
Eigen::Matrix< scalar_type, Eigen::Dynamic, 1 > vector_element_type
typename float_traits< scalar_type >::type float_type
void initialize(const geo_type &omega_K) const
void evaluate(const geo_type &omega_K, const geo_element &K, vector_element_type &uk) const
typename rep::vector_element_type vector_element_type
field_lazy_mult_form(const FormExpr &a_expr, const FieldExpr &u_expr)
band_basic< float_type, memory_type > band_type
field_lazy_trans_mult_form_rep(const FormExpr &a_expr, const FieldExpr &u_expr)
void initialize(const geo_type &omega_K) const
void evaluate(const geo_type &omega_K, const geo_element &K, vector_element_type &uk) const
Eigen::Matrix< scalar_type, Eigen::Dynamic, Eigen::Dynamic > matrix_element_type
Eigen::Matrix< scalar_type, Eigen::Dynamic, 1 > vector_element_type
typename float_traits< scalar_type >::type float_type
void initialize(const geo_type &omega_K) const
void evaluate(const geo_type &omega_K, const geo_element &K, vector_element_type &uk) const
field_lazy_trans_mult_form(const FormExpr &a_expr, const FieldExpr &u_expr)
typename rep::vector_element_type vector_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
rheolef::std type
rheolef::std value
#define trace_macro(message)
Definition: dis_macros.h:111
#define fatal_macro(message)
Definition: dis_macros.h:33
#define warning_macro(message)
Definition: dis_macros.h:53
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.
csr< T, sequential > operator*(const T &lambda, const csr< T, sequential > &a)
Definition: csr.h:437
Definition: leveque.h:25