Rheolef  7.1
an efficient C++ finite element environment
integrate.h
Go to the documentation of this file.
1 #ifndef _RHEO_INTEGRATE_H
2 #define _RHEO_INTEGRATE_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 
24 namespace rheolef {
155 } // namespace rheolef
156 
157 // Implementation note
158 // -------------------
159 // SUMMARY:
160 // 1. numeric integration
161 // 1.1. general integration of a nonlinear expression
162 // 1.2. measure of the domain
163 // 1.3. when the valued result type is undetermined
164 // 2. field-result integration of a variational expression
165 // 2.1. general call
166 // 2.2. missing domain
167 // 2.3. subdomain by its name
168 // 3. form-result integration of a variational expression
169 // 3.1. general call
170 // 3.2. missing domain
171 // 3.3. subdomain by its name
172 // 4. variational integration: on a band
173 //
174 #include "rheolef/field_expr.h"
175 #include "rheolef/field_expr_variational.h"
176 #include "rheolef/form_expr_variational.h"
177 
178 #include "rheolef/field_expr_value_assembly.h"
179 #include "rheolef/field_vf_assembly.h"
180 #include "rheolef/form_vf_assembly.h"
181 #include "rheolef/form_expr_quadrature.h"
182 #include "rheolef/field_expr_quadrature.h"
183 
184 #include "rheolef/functor.h" // used to convert functions to functors
185 
186 namespace rheolef {
187 
188 // ---------------------------------------------------
189 // 1. numeric integration
190 // ---------------------------------------------------
191 // 1.1. general integration of a nonlinear expression
192 // ---------------------------------------------------
193 template <class T, class M, class Expr,
194  class Result = typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr>::type::value_type>
195 inline
196 typename std::enable_if<
197  details::is_field_expr_v2_nonlinear_arg<Expr>::value
198  && ! is_undeterminated<Result>::value,
199  Result
200 >::type
202 integrate (const geo_basic<T,M>& omega, const Expr& expr, const integrate_option& iopt,
203  Result dummy = Result())
204 {
206  if (omega.map_dimension() < omega.get_background_geo().map_dimension()) {
207  omega.get_background_geo().neighbour_guard();
208  }
209  Result result(0);
210  field_expr_v2_value_assembly (omega, wrap_t(expr), iopt, result);
211  return result;
212 }
213 // ---------------------------------------------------
214 // 1.2. measure of the domain
215 // ---------------------------------------------------
216 template <class T, class M>
217 T
220 {
221  if (iopt.get_order() == std::numeric_limits<integrate_option::size_type>::max()) {
222  iopt.set_order(0);
223  }
225  return integrate (omega, one, iopt);
226 }
227 // ---------------------------------------------------
228 // 1.3. when the valued result type is undetermined
229 // ---------------------------------------------------
230 // TODO: return a overdetermined<T> value that is an union of all possibilities with a valued_tag
231 template<class T, class M, class Expr>
232 inline
233 typename std::enable_if<
234  details::is_field_expr_v2_nonlinear_arg<Expr>::value
235  && is_undeterminated<typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr>::type::value_type>::value,
236  typename scalar_traits<typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr>::type::value_type>::type
237 >::type
239 integrate (const geo_basic<T,M>& omega, const Expr& expr, const integrate_option& iopt)
240 {
241  typedef typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr>::type::value_type undef_t;
242  typedef typename scalar_traits<undef_t>::type scalar_type;
243  switch (expr.valued_tag()) {
244  case space_constant::scalar: {
245  return integrate (omega, expr, iopt, scalar_type());
246  }
247  // others type: problem on how to return a run-type type ?
248  // TODO: return an overdetermined union type that convert to one of scalar, point, tensor, etc ?
249  default:
250  warning_macro ("Expr="<<pretty_typename_macro(Expr));
251  error_macro ("integrate: not yet for `"
252  << space_constant::valued_name (expr.valued_tag())
253  << "' valued expression");
254  return 0;
255  }
256 }
257 // -------------------------------------------------------
258 // 2. field-result integration of a variational expression
259 // -------------------------------------------------------
260 // 2.1. general call
261 // -------------------------------------------------------
262 template <class T, class M, class Expr>
263 inline
264 typename
265 std::enable_if<
266  details::is_field_expr_quadrature_arg<Expr>::value
267  ,field_basic<T,M>
268 >::type
271  const geo_basic<T,M>& domain,
272  const Expr& expr,
273  const integrate_option& iopt = integrate_option())
274 {
276  lh.assembly (domain, expr, iopt);
277  return lh;
278 }
279 template <class T, class M, class Expr>
280 inline
281 typename
282 std::enable_if<
283  details::is_field_expr_v2_variational_arg<Expr>::value
284  ,field_basic<T,M>
285 >::type
288  const geo_basic<T,M>& domain,
289  const Expr& expr,
290  const integrate_option& fopt = integrate_option())
291 {
293  return integrate (domain, expr_quad, fopt);
294 }
295 // ----------------------------------------------
296 // 2.2. missing domain
297 // ----------------------------------------------
298 template <class Expr>
299 inline
300 typename
301 std::enable_if<
302  details::is_field_expr_quadrature_arg<Expr>::value
303  ,field_basic <typename Expr::scalar_type, typename Expr::memory_type>
304 >::type
307  const Expr& expr,
308  const integrate_option& iopt = integrate_option())
309 {
312  dom = expr.get_vf_space().get_constitution().get_geo();
313  lh.assembly (dom, expr, iopt);
314  return lh;
315 }
316 template <class Expr>
317 inline
318 typename
319 std::enable_if<
320  details::is_field_expr_v2_variational_arg<Expr>::value
321  ,field_basic <typename Expr::scalar_type, typename Expr::memory_type>
322 >::type
325  const Expr& expr,
326  const integrate_option& fopt = integrate_option())
327 {
329  return integrate (expr_quad, fopt);
330 }
331 // ----------------------------------------------
332 // 2.3. subdomain by its name
333 // ----------------------------------------------
334 template <class Expr>
335 inline
336 typename
337 std::enable_if<
338  details::is_field_expr_quadrature_arg<Expr>::value
339  ,field_basic <typename Expr::scalar_type, typename Expr::memory_type>
340 >::type
343  const std::string& domname,
344  const Expr& expr,
345  const integrate_option& iopt = integrate_option())
346 {
349  dom = expr.get_vf_space().get_constitution().get_geo() [domname];
350  lh.assembly (dom, expr, iopt);
351  return lh;
352 }
353 template <class Expr>
354 inline
355 typename
356 std::enable_if<
357  details::is_field_expr_v2_variational_arg<Expr>::value
358  ,field_basic <typename Expr::scalar_type, typename Expr::memory_type>
359 >::type
362  const std::string& domname,
363  const Expr& expr,
364  const integrate_option& fopt = integrate_option())
365 {
367  return integrate (domname, expr_quad, fopt);
368 }
369 // -------------------------------------------------------
370 // 3. form-result integration of a variational expression
371 // -------------------------------------------------------
372 // 3.1. general call
373 // -------------------------------------------------------
374 template <class T, class M, class Expr>
375 inline
376 typename
377 std::enable_if<
378  details::is_form_expr_quadrature_arg<Expr>::value
379  ,form_basic <typename Expr::scalar_type, typename Expr::memory_type>
380 >::type
383  const geo_basic<T,M>& domain,
384  const Expr& expr,
385  const integrate_option& fopt = integrate_option())
386 {
388  a.assembly (domain, expr, fopt);
389  return a;
390 }
391 template <class T, class M, class Expr>
392 inline
393 typename
394 std::enable_if<
395  details::is_form_expr_v2_variational_arg<Expr>::value
396  ,form_basic <typename Expr::scalar_type, typename Expr::memory_type>
397 >::type
400  const geo_basic<T,M>& domain,
401  const Expr& expr,
402  const integrate_option& fopt = integrate_option())
403 {
405  return integrate (domain, expr_quad, fopt);
406 }
407 // ----------------------------------------------
408 // 3.2. missing domain
409 // ----------------------------------------------
410 template <class Expr>
411 inline
412 typename
413 std::enable_if<
414  details::is_form_expr_quadrature_arg<Expr>::value
415  ,form_basic <typename Expr::scalar_type, typename Expr::memory_type>
416 >::type
419  const Expr& expr,
420  const integrate_option& fopt = integrate_option())
421 {
424  dom_trial = expr.get_trial_space().get_constitution().get_geo(),
425  dom_test = expr.get_test_space().get_constitution().get_geo(),
426  dom;
427  // dom = intersection of trial & test domain definition
428  if (dom_trial.is_broken() && dom_test.is_broken() &&
429  expr.get_trial_space().get_constitution().is_hierarchical() &&
430  expr.get_test_space().get_constitution().is_hierarchical() ) {
431  dom = dom_test.get_background_geo();
432  } else if (dom_trial.name() == dom_test.name() ||
433  dom_trial.name() == dom_test.get_background_geo().name() ||
434  dom_trial.is_broken()) {
435  dom = dom_test;
436  } else if (dom_test.name() == dom_trial.get_background_geo().name() ||
437  dom_test.is_broken()) {
438  dom = dom_trial;
439  } else {
440  error_macro("integrate: incompatible domains: trial \""<<dom_trial.name()
441  << "\" and \"" << dom_test.name() << "\"");
442  }
443  trace_macro ("dom_trial="<<dom_trial.name()<<" dom_test="<<dom_test.name()<<" -> dom="<<dom.name());
444  a.assembly (dom, expr, fopt);
445  return a;
446 }
447 template <class Expr>
448 inline
449 typename
450 std::enable_if<
451  details::is_form_expr_v2_variational_arg<Expr>::value
452  ,form_basic <typename Expr::scalar_type, typename Expr::memory_type>
453 >::type
456  const Expr& expr,
457  const integrate_option& fopt = integrate_option())
458 {
460  return integrate (expr_quad, fopt);
461 }
462 // ----------------------------------------------
463 // 3.3. subdomain by its name
464 // ----------------------------------------------
465 template <class Expr>
466 inline
467 typename
468 std::enable_if<
469  details::is_form_expr_quadrature_arg<Expr>::value
470  ,form_basic <typename Expr::scalar_type, typename Expr::memory_type>
471 >::type
474  const std::string& domname,
475  const Expr& expr,
476  const integrate_option& fopt = integrate_option())
477 {
480  dom = expr.get_trial_space().get_constitution().get_background_geo()[domname];
481  a.assembly (dom, expr, fopt);
482  return a;
483 }
484 template <class Expr>
485 inline
486 typename
487 std::enable_if<
488  details::is_form_expr_v2_variational_arg<Expr>::value
489  ,form_basic <typename Expr::scalar_type, typename Expr::memory_type>
490 >::type
493  const std::string& domname,
494  const Expr& expr,
495  const integrate_option& fopt = integrate_option())
496 {
498  return integrate (domname, expr_quad, fopt);
499 }
500 // ----------------------------------------------
501 // 4. variational integration: on a band
502 // ----------------------------------------------
503 template <class T, class M, class Expr>
504 inline
505 typename
506 std::enable_if<
507  details::is_field_expr_quadrature_arg<Expr>::value
508  ,field_basic <typename Expr::scalar_type, typename Expr::memory_type>
509 >::type
512  const band_basic<T,M>& gh,
513  const Expr& expr,
514  const integrate_option& iopt = integrate_option())
515 {
517  lh.assembly (gh, expr, iopt);
518  return lh;
519 }
520 template <class T, class M, class Expr>
521 inline
522 typename
523 std::enable_if<
524  details::is_field_expr_v2_variational_arg<Expr>::value
525  ,field_basic <typename Expr::scalar_type, typename Expr::memory_type>
526 >::type
529  const band_basic<T,M>& gh,
530  const Expr& expr,
531  const integrate_option& iopt = integrate_option())
532 {
533 
535  return integrate (gh, expr_quad, iopt);
536 }
537 template <class T, class M, class Expr>
538 inline
539 typename
540 std::enable_if<
541  details::is_form_expr_quadrature_arg<Expr>::value
542  ,form_basic <typename Expr::scalar_type, typename Expr::memory_type>
543 >::type
546  const band_basic<T,M>& gh,
547  const Expr& expr,
548  const integrate_option& fopt = integrate_option())
549 {
551  a.assembly (gh, expr, fopt);
552  return a;
553 }
554 template <class T, class M, class Expr>
555 inline
556 typename
557 std::enable_if<
558  details::is_form_expr_v2_variational_arg<Expr>::value
559  ,form_basic <typename Expr::scalar_type, typename Expr::memory_type>
560 >::type
563  const band_basic<T,M>& gh,
564  const Expr& expr,
565  const integrate_option& fopt = integrate_option())
566 {
568  return integrate (gh, expr_quad, fopt);
569 }
570 
571 }// namespace rheolef
572 #endif // _RHEO_INTEGRATE_H
rheolef::geo_basic
generic mesh with rerefence counting
Definition: geo.h:1089
rheolef::domain
geo domain
Definition: geo_domain_indirect.h:58
rheolef::integrate_option::get_order
size_t get_order() const
Definition: integrate_option.h:242
rheolef::details::f_constant
Definition: space_constant.h:297
warning_macro
#define warning_macro(message)
Definition: dis_macros.h:53
gh
field gh(Float epsilon, Float t, const field &uh, const test &v)
Definition: burgers_diffusion_operators.icc:37
mkgeo_ball.expr
expr
Definition: mkgeo_ball.sh:361
rheolef::integrate_option::set_order
void set_order(size_t r)
Definition: integrate_option.h:254
rheolef::dummy
static iorheo::force_initialization dummy
Definition: iorheo.cc:147
rheolef::integrate
std::enable_if< details::is_field_expr_v2_nonlinear_arg< Expr >::value &&! is_undeterminated< Result >::value, Result >::type integrate(const geo_basic< T, M > &omega, const Expr &expr, const integrate_option &iopt, Result dummy=Result())
see the integrate page for the full documentation
Definition: integrate.h:202
rheolef::details::field_expr_v2_nonlinear_terminal_wrapper_traits::type
Expr type
Definition: field_expr_terminal.h:98
rheolef::value
rheolef::std value
rheolef::integrate_option
see the integrate_option page for the full documentation
Definition: integrate_option.h:125
rheolef::space_constant::scalar
@ scalar
Definition: space_constant.h:136
a
Definition: diffusion_isotropic.h:25
rheolef::type
rheolef::std type
rheolef::field_basic
Definition: field.h:235
lh
field lh(Float epsilon, Float t, const test &v)
Definition: burgers_diffusion_operators.icc:25
rheolef::space_constant::valued_name
const std::string & valued_name(valued_type valued_tag)
Definition: space_constant.cc:43
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
error_macro
#define error_macro(message)
Definition: dis_macros.h:49
rheolef::scalar_type
typename scalar_traits< value_type >::type scalar_type
Definition: field_expr_recursive.h:500
rheolef::details::field_expr_v2_value_assembly
void field_expr_v2_value_assembly(const geo_basic< T, M > &omega, const Expr &expr, const integrate_option &iopt, Result &result)
Definition: field_expr_value_assembly.h:55
rheolef::form_basic
Definition: form.h:144
rheolef::scalar_traits::type
T type
Definition: point.h:324
rheolef::details::field_expr_quadrature_on_element
Definition: field_expr_quadrature.h:67
mkgeo_ball.a
int a
Definition: mkgeo_ball.sh:151
rheolef::details::form_expr_quadrature_on_element
Definition: form_expr_quadrature.h:65
M
Expr1::memory_type M
Definition: vec_expr_v2.h:416
T
Expr1::float_type T
Definition: field_expr.h:261
trace_macro
#define trace_macro(message)
Definition: dis_macros.h:111