Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
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
34namespace rheolef {
35
36// -------------------------------------------------------------------
37// 1. form_lazy*field_lazy -> field_lazy
38// -------------------------------------------------------------------
39namespace details {
40
41template<class FormExpr, class FieldExpr>
43public :
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:
74protected:
75 mutable FormExpr _a_expr;
76 mutable FieldExpr _u_expr;
80};
81// inlined;
82template<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}
93template<class FormExpr, class FieldExpr>
94void
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}
118template<class FormExpr, class FieldExpr>
119void
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}
143template<class FormExpr, class FieldExpr>
144class 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>> {
146public :
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;
158 using vector_element_type = typename rep::vector_element_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;
183template<class FormExpr, class FieldExpr>
184struct is_field_lazy <field_lazy_mult_form<FormExpr,FieldExpr> > : std::true_type {};
185
186}// namespace details
187
189template<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>
193operator* (const FormExpr& a, const FieldExpr& u)
194{
196}
198template<
199 class FormExpr
200 ,class Sfinae = typename std::enable_if<details::is_form_lazy<FormExpr>::value, FormExpr>::type
201>
202details::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// -------------------------------------------------------------------
222namespace details {
223
224template<class FormExpr, class FieldExpr>
226public :
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:
257protected:
258 mutable FormExpr _a_expr;
259 mutable FieldExpr _u_expr;
263};
264// inlined;
265template<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}
276template<class FormExpr, class FieldExpr>
277void
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}
301template<class FormExpr, class FieldExpr>
302void
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}
326template<class FormExpr, class FieldExpr>
327class 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>> {
329public :
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;
341 using vector_element_type = typename rep::vector_element_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;
366template<class FormExpr, class FieldExpr>
367struct 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 */
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
field_lazy_mult_form_rep< FormExpr, FieldExpr > rep
typename rep::vector_element_type vector_element_type
field_lazy_mult_form(const FormExpr &a_expr, const FieldExpr &u_expr)
field_lazy_trans_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
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
field_lazy_trans_mult_form_rep< FormExpr, FieldExpr > rep
see the geo_element page for the full documentation
reference_element::size_type size_type
size_type dis_ie() const
#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
STL namespace.
Definition leveque.h:25