Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
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
37namespace rheolef {
38
39// -------------------------------------------------------------------
40// 1. concept & base class
41// -------------------------------------------------------------------
42namespace 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
49template<class Derived>
51public:
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>
71
72protected:
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)
86namespace details {
87
88template<class Expr>
90public :
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),
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:
132protected:
134 mutable Expr _expr;
139};
140// inlined;
141template<class Expr>
142void
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}
155template<class Expr>
156void
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
174template<class Expr>
175class form_lazy_terminal_integrate: public smart_pointer_nocopy<form_lazy_terminal_integrate_rep<Expr>>,
176 public form_lazy_base <form_lazy_terminal_integrate<Expr>> {
177public :
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;
189 using matrix_element_type = typename rep::matrix_element_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;
219template<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// ----------------------------------------------
225namespace details {
226
227template<class Expr>
229public :
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:
263protected:
265 mutable Expr _expr;
270};
271// inlined;
272template<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}
285template<class Expr>
286void
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}
300template<class Expr>
301void
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}
318template<class Expr>
319class 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>> {
321public :
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;
333 using matrix_element_type = typename rep::matrix_element_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;
362template<class Expr> struct is_form_lazy <form_lazy_terminal_integrate_band<Expr> > : std::true_type {};
363
364}// namespace details
365
366template <class Expr>
367typename
368std::enable_if<
370 ,details::form_lazy_terminal_integrate_band <Expr>
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}
381template <class Expr>
382typename
383std::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
std::enable_if< details::is_field_lazy< FieldExpr >::value, field_lazy_trans_mult_form< Derived, FieldExpr > >::type trans_mult(const FieldExpr &u_expr) const
field_lazy_trans_mult_form< Derived, field_lazy_terminal_field< T, M > > trans_mult(const field_basic< T, M > &uh) 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
form_lazy_terminal_integrate_band(const band_type &gh, const Expr &expr, const integrate_option &iopt)
form_lazy_terminal_integrate_band_rep< Expr > rep
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
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)
form_lazy_terminal_integrate_rep< Expr > rep
typename rep::matrix_element_type matrix_element_type
see the geo_element page for the full documentation
reference_element::size_type size_type
size_type dis_ie() const
see the integrate_option page for the full documentation
#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
STL namespace.