Rheolef  7.2
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. scalar-result 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 // 5. lazy-field-result integration of a variational expression
174 // 5.1. general call
175 // 5.2. missing domain
176 // 5.3. subdomain by its name
177 // 6. lazy-form-result integration of a variational expression
178 // 6.1. general call
179 // 6.2. missing domain
180 // 6.3. subdomain by its name
181 //
182 #include "rheolef/field_expr.h"
183 #include "rheolef/field_expr_variational.h"
184 #include "rheolef/form_expr_variational.h"
185 
186 #include "rheolef/field_expr_value_assembly.h"
187 #include "rheolef/field_vf_assembly.h"
188 #include "rheolef/form_vf_assembly.h"
189 #include "rheolef/form_expr_quadrature.h"
190 #include "rheolef/field_expr_quadrature.h"
191 #include "rheolef/form_lazy_expr.h"
192 
193 #include "rheolef/functor.h" // used to convert functions to functors
194 
195 namespace rheolef {
196 
197 // ---------------------------------------------------
198 // 1. scalar-result integration
199 // ---------------------------------------------------
200 // 1.1. general integration of a nonlinear expression
201 // ---------------------------------------------------
202 template <class T, class M, class Expr,
203  class Result = typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr>::type::value_type>
204 inline
205 typename std::enable_if<
206  details::is_field_expr_v2_nonlinear_arg<Expr>::value
207  && ! is_undeterminated<Result>::value,
208  Result
209 >::type
211 integrate (const geo_basic<T,M>& omega, const Expr& expr, const integrate_option& iopt,
212  Result dummy = Result())
213 {
215  if (omega.map_dimension() < omega.get_background_geo().map_dimension()) {
216  omega.get_background_geo().neighbour_guard();
217  }
218  Result result(0);
219  field_expr_v2_value_assembly (omega, wrap_t(expr), iopt, result);
220  return result;
221 }
222 // ---------------------------------------------------
223 // 1.2. measure of the domain
224 // ---------------------------------------------------
225 template <class T, class M>
226 T
228 integrate (const geo_basic<T,M>& omega, integrate_option&& iopt = integrate_option())
229 {
230  if (iopt.get_order() == std::numeric_limits<integrate_option::size_type>::max()) {
231  iopt.set_order(0);
232  }
234  return integrate (omega, one, iopt);
235 }
236 // ---------------------------------------------------
237 // 1.3. when the valued result type is undetermined
238 // ---------------------------------------------------
239 // TODO: return a overdetermined<T> value that is an union of all possibilities with a valued_tag
240 template<class T, class M, class Expr>
241 inline
242 typename std::enable_if<
243  details::is_field_expr_v2_nonlinear_arg<Expr>::value
244  && is_undeterminated<typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr>::type::value_type>::value,
245  typename scalar_traits<typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr>::type::value_type>::type
246 >::type
248 integrate (const geo_basic<T,M>& omega, const Expr& expr, const integrate_option& iopt)
249 {
250  typedef typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr>::type::value_type undef_t;
251  typedef typename scalar_traits<undef_t>::type scalar_type;
252  switch (expr.valued_tag()) {
253  case space_constant::scalar: {
254  return integrate (omega, expr, iopt, scalar_type());
255  }
256  // others type: problem on how to return a run-type type ?
257  // TODO: return an overdetermined union type that convert to one of scalar, point, tensor, etc ?
258  default:
259  warning_macro ("Expr="<<pretty_typename_macro(Expr));
260  error_macro ("integrate: not yet for `"
261  << space_constant::valued_name (expr.valued_tag())
262  << "' valued expression");
263  return 0;
264  }
265 }
266 // -------------------------------------------------------
267 // 2. field-result integration of a variational expression
268 // -------------------------------------------------------
269 // 2.1. general call
270 // -------------------------------------------------------
271 template <class T, class M, class Expr>
272 inline
273 typename
274 std::enable_if<
275  details::is_field_expr_quadrature_arg<Expr>::value
276  ,field_basic<T,M>
277 >::type
280  const geo_basic<T,M>& domain,
281  const Expr& expr,
282  const integrate_option& iopt = integrate_option())
283 {
285  lh.do_integrate (domain, expr, iopt);
286  return lh;
287 }
288 template <class T, class M, class Expr>
289 inline
290 typename
291 std::enable_if<
292  details::is_field_expr_v2_variational_arg<Expr>::value
293  ,field_basic<T,M>
294 >::type
297  const geo_basic<T,M>& domain,
298  const Expr& expr,
299  const integrate_option& fopt = integrate_option())
300 {
302  return integrate (domain, expr_quad, fopt);
303 }
304 // ----------------------------------------------
305 // 2.2. missing domain
306 // ----------------------------------------------
307 template <class Expr>
308 inline
309 typename
310 std::enable_if<
311  details::is_field_expr_quadrature_arg<Expr>::value
312  ,field_basic <typename Expr::scalar_type, typename Expr::memory_type>
313 >::type
316  const Expr& expr,
317  const integrate_option& iopt = integrate_option())
318 {
321  dom = expr.get_vf_space().get_constitution().get_geo();
322  lh.do_integrate (dom, expr, iopt);
323  return lh;
324 }
325 template <class Expr>
326 inline
327 typename
328 std::enable_if<
329  details::is_field_expr_v2_variational_arg<Expr>::value
330  ,field_basic <typename Expr::scalar_type, typename Expr::memory_type>
331 >::type
334  const Expr& expr,
335  const integrate_option& fopt = integrate_option())
336 {
338  return integrate (expr_quad, fopt);
339 }
340 // ----------------------------------------------
341 // 2.3. subdomain by its name
342 // ----------------------------------------------
343 template <class Expr>
344 inline
345 typename
346 std::enable_if<
347  details::is_field_expr_quadrature_arg<Expr>::value
348  ,field_basic <typename Expr::scalar_type, typename Expr::memory_type>
349 >::type
352  const std::string& domname,
353  const Expr& expr,
354  const integrate_option& iopt = integrate_option())
355 {
358  dom = expr.get_vf_space().get_constitution().get_geo() [domname];
359  lh.do_integrate (dom, expr, iopt);
360  return lh;
361 }
362 template <class Expr>
363 inline
364 typename
365 std::enable_if<
366  details::is_field_expr_v2_variational_arg<Expr>::value
367  ,field_basic <typename Expr::scalar_type, typename Expr::memory_type>
368 >::type
371  const std::string& domname,
372  const Expr& expr,
373  const integrate_option& fopt = integrate_option())
374 {
376  return integrate (domname, expr_quad, fopt);
377 }
378 // -------------------------------------------------------
379 // 3. form-result integration of a variational expression
380 // -------------------------------------------------------
381 // 3.1. general call
382 // -------------------------------------------------------
383 template <class T, class M, class Expr>
384 inline
385 typename
386 std::enable_if<
387  details::is_form_expr_quadrature_arg<Expr>::value
388  ,form_basic <typename Expr::scalar_type, typename Expr::memory_type>
389 >::type
392  const geo_basic<T,M>& domain,
393  const Expr& expr,
394  const integrate_option& fopt = integrate_option())
395 {
397  a.do_integrate (domain, expr, fopt);
398  return a;
399 }
400 template <class T, class M, class Expr>
401 inline
402 typename
403 std::enable_if<
404  details::is_form_expr_v2_variational_arg<Expr>::value
405  ,form_basic <typename Expr::scalar_type, typename Expr::memory_type>
406 >::type
409  const geo_basic<T,M>& domain,
410  const Expr& expr,
411  const integrate_option& fopt = integrate_option())
412 {
414  return integrate (domain, expr_quad, fopt);
415 }
416 // ----------------------------------------------
417 // 3.2. missing domain
418 // ----------------------------------------------
419 template <class Expr>
420 inline
421 typename
422 std::enable_if<
423  details::is_form_expr_quadrature_arg<Expr>::value
424  ,form_basic <typename Expr::scalar_type, typename Expr::memory_type>
425 >::type
428  const Expr& expr,
429  const integrate_option& fopt = integrate_option())
430 {
433  dom_trial = expr.get_trial_space().get_constitution().get_geo(),
434  dom_test = expr.get_test_space().get_constitution().get_geo(),
435  dom;
436  // dom = intersection of trial & test domain definition
437  if (dom_trial.is_broken() && dom_test.is_broken() &&
438  expr.get_trial_space().get_constitution().is_hierarchical() &&
439  expr.get_test_space().get_constitution().is_hierarchical() ) {
440  dom = dom_test.get_background_geo();
441  } else if (dom_trial.name() == dom_test.name() ||
442  dom_trial.name() == dom_test.get_background_geo().name() ||
443  dom_trial.is_broken()) {
444  dom = dom_test;
445  } else if (dom_test.name() == dom_trial.get_background_geo().name() ||
446  dom_test.is_broken()) {
447  dom = dom_trial;
448  } else {
449  error_macro("integrate: incompatible domains: trial \""<<dom_trial.name()
450  << "\" and \"" << dom_test.name() << "\"");
451  }
452  trace_macro ("dom_trial="<<dom_trial.name()<<" dom_test="<<dom_test.name()<<" -> dom="<<dom.name());
453  a.do_integrate (dom, expr, fopt);
454  return a;
455 }
456 template <class Expr>
457 inline
458 typename
459 std::enable_if<
460  details::is_form_expr_v2_variational_arg<Expr>::value
461  ,form_basic <typename Expr::scalar_type, typename Expr::memory_type>
462 >::type
465  const Expr& expr,
466  const integrate_option& fopt = integrate_option())
467 {
469  return integrate (expr_quad, fopt);
470 }
471 // ----------------------------------------------
472 // 3.3. subdomain by its name
473 // ----------------------------------------------
474 template <class Expr>
475 inline
476 typename
477 std::enable_if<
478  details::is_form_expr_quadrature_arg<Expr>::value
479  ,form_basic <typename Expr::scalar_type, typename Expr::memory_type>
480 >::type
483  const std::string& domname,
484  const Expr& expr,
485  const integrate_option& fopt = integrate_option())
486 {
489  dom = expr.get_trial_space().get_constitution().get_background_geo()[domname];
490  a.do_integrate (dom, expr, fopt);
491  return a;
492 }
493 template <class Expr>
494 inline
495 typename
496 std::enable_if<
497  details::is_form_expr_v2_variational_arg<Expr>::value
498  ,form_basic <typename Expr::scalar_type, typename Expr::memory_type>
499 >::type
502  const std::string& domname,
503  const Expr& expr,
504  const integrate_option& fopt = integrate_option())
505 {
507  return integrate (domname, expr_quad, fopt);
508 }
509 // ----------------------------------------------
510 // 4. variational integration: on a band
511 // ----------------------------------------------
512 template <class T, class M, class Expr>
513 inline
514 typename
515 std::enable_if<
516  details::is_field_expr_quadrature_arg<Expr>::value
517  ,field_basic <typename Expr::scalar_type, typename Expr::memory_type>
518 >::type
521  const band_basic<T,M>& gh,
522  const Expr& expr,
523  const integrate_option& iopt = integrate_option())
524 {
526  lh.do_integrate (gh, expr, iopt);
527  return lh;
528 }
529 template <class T, class M, class Expr>
530 inline
531 typename
532 std::enable_if<
533  details::is_field_expr_v2_variational_arg<Expr>::value
534  ,field_basic <typename Expr::scalar_type, typename Expr::memory_type>
535 >::type
538  const band_basic<T,M>& gh,
539  const Expr& expr,
540  const integrate_option& iopt = integrate_option())
541 {
542 
544  return integrate (gh, expr_quad, iopt);
545 }
546 template <class T, class M, class Expr>
547 inline
548 typename
549 std::enable_if<
550  details::is_form_expr_quadrature_arg<Expr>::value
551  ,form_basic <typename Expr::scalar_type, typename Expr::memory_type>
552 >::type
555  const band_basic<T,M>& gh,
556  const Expr& expr,
557  const integrate_option& fopt = integrate_option())
558 {
560  a.do_integrate (gh, expr, fopt);
561  return a;
562 }
563 template <class T, class M, class Expr>
564 inline
565 typename
566 std::enable_if<
567  details::is_form_expr_v2_variational_arg<Expr>::value
568  ,form_basic <typename Expr::scalar_type, typename Expr::memory_type>
569 >::type
572  const band_basic<T,M>& gh,
573  const Expr& expr,
574  const integrate_option& fopt = integrate_option())
575 {
577  return integrate (gh, expr_quad, fopt);
578 }
579 #ifdef TO_CLEAN
580 // -------------------------------------------------------
581 // 5. field-result integration of a variational expression
582 // -------------------------------------------------------
583 // 5.1. general call
584 // -------------------------------------------------------
585 template <class Expr>
586 inline
587 typename
588 std::enable_if<
589  details::is_field_expr_quadrature_arg<Expr>::value
590  ,details::field_lazy_terminal_integrate <Expr>
591 >::type
594  const typename Expr::geo_type& domain,
595  const Expr& expr,
596  const integrate_option& iopt = integrate_option())
597 {
598  return details::field_lazy_terminal_integrate<Expr> (domain, expr, iopt);
599 }
600 template <class Expr>
601 inline
602 typename
603 std::enable_if<
604  details::is_field_expr_v2_variational_arg<Expr>::value
605  ,details::field_lazy_terminal_integrate <details::field_expr_quadrature_on_element<Expr>>
606 >::type
609  const geo_basic<typename Expr::scalar_type, typename Expr::memory_type>& domain,
610  const Expr& expr,
611  const integrate_option& iopt = integrate_option())
612 {
613  details::field_expr_quadrature_on_element<Expr> expr_quad(expr);
614  return lazy_integrate (domain, expr_quad, iopt);
615 }
616 // ----------------------------------------------
617 // 5.2. missing domain
618 // ----------------------------------------------
619 template <class Expr>
620 inline
621 typename
622 std::enable_if<
623  details::is_field_expr_quadrature_arg<Expr>::value
624  ,details::field_lazy_terminal_integrate <Expr>
625 >::type
628  const Expr& expr,
629  const integrate_option& fopt = integrate_option())
630 {
631  return details::field_lazy_terminal_integrate<Expr> (expr, iopt);
632 }
633 template <class Expr>
634 inline
635 typename
636 std::enable_if<
637  details::is_field_expr_v2_variational_arg<Expr>::value
638  ,details::field_lazy_terminal_integrate <details::field_expr_quadrature_on_element<Expr>>
639 >::type
642  const Expr& expr,
643  const integrate_option& fopt = integrate_option())
644 {
645  details::field_expr_quadrature_on_element<Expr> expr_quad(expr);
646  return lazy_integrate (expr_quad, fopt);
647 }
648 // ----------------------------------------------
649 // 5.3. subdomain by its name
650 // ----------------------------------------------
651 template <class Expr>
652 inline
653 typename
654 std::enable_if<
655  details::is_field_expr_quadrature_arg<Expr>::value
656  ,details::field_lazy_terminal_integrate <Expr>
657 >::type
660  const std::string& domname,
661  const Expr& expr,
662  const integrate_option& fopt = integrate_option())
663 {
664  return details::field_lazy_terminal_integrate<Expr> (domname, expr, iopt);
665 }
666 template <class Expr>
667 inline
668 typename
669 std::enable_if<
670  details::is_field_expr_v2_variational_arg<Expr>::value
671  ,details::field_lazy_terminal_integrate <details::field_expr_quadrature_on_element<Expr>>
672 >::type
675  const std::string& domname,
676  const Expr& expr,
677  const integrate_option& fopt = integrate_option())
678 {
679  details::field_expr_quadrature_on_element<Expr> expr_quad(expr);
680  return lazy_integrate (domname, expr_quad, fopt);
681 }
682 #endif // TO_CLEAN
683 // -------------------------------------------------------
684 // 6. lazy_form-result integration of a variational expression
685 // -------------------------------------------------------
686 // 6.1. general call
687 // -------------------------------------------------------
688 template <class Expr>
689 inline
690 typename
691 std::enable_if<
692  details::is_form_expr_quadrature_arg<Expr>::value
693  ,details::form_lazy_terminal_integrate <Expr>
694 >::type
698  const Expr& expr,
699  const integrate_option& iopt = integrate_option())
700 {
702 }
703 template <class Expr>
704 inline
705 typename
706 std::enable_if<
707  details::is_form_expr_v2_variational_arg<Expr>::value
708  ,details::form_lazy_terminal_integrate <details::form_expr_quadrature_on_element<Expr>>
709 >::type
713  const Expr& expr,
714  const integrate_option& iopt = integrate_option())
715 {
717  return lazy_integrate (domain, expr_quad, iopt);
718 }
719 // ----------------------------------------------
720 // 6.2. missing domain
721 // ----------------------------------------------
722 template <class Expr>
723 inline
724 typename
725 std::enable_if<
726  details::is_form_expr_quadrature_arg<Expr>::value
727  ,details::form_lazy_terminal_integrate <Expr>
728 >::type
731  const Expr& expr,
732  const integrate_option& fopt = integrate_option())
733 {
734  // TODO: improve the automatic determination of the domain
735  // by an union that operates recrusively
736  constexpr bool has_on_local_sides = details::is_form_expr_quadrature_on_side_arg<Expr>::value;
738  dom_trial = expr.get_trial_space().get_constitution().get_geo(),
739  dom_test = expr. get_test_space().get_constitution().get_geo(),
740  dom;
741  // dom = intersection of trial & test domain definition
742  if (dom_trial.is_broken() && dom_test.is_broken() &&
743  expr.get_trial_space().get_constitution().is_hierarchical() &&
744  expr. get_test_space().get_constitution().is_hierarchical() ) {
745  dom = dom_test.get_background_geo();
746  } else if (dom_trial.name() == dom_test.name() ||
747  dom_trial.name() == dom_test.get_background_geo().name() ||
748  dom_trial.is_broken()) {
749  dom = has_on_local_sides ? dom_trial : dom_test;
750  } else if (dom_test.name() == dom_trial.get_background_geo().name() ||
751  dom_test.is_broken()) {
752  dom = has_on_local_sides ? dom_test : dom_trial;
753  } else {
754  error_macro("integrate: incompatible domains: trial \""<<dom_trial.name()
755  << "\" and \"" << dom_test.name() << "\"");
756  }
757  trace_macro("Expr="<<pretty_typename_macro(Expr));
758  trace_macro ("space_trial="<<expr.get_trial_space().name()<<" space_test="<<expr. get_test_space().name());
759  trace_macro ("dom_trial="<<dom_trial.name()<<" dom_test="<<dom_test.name()<<" -> dom="<<dom.name());
760  return lazy_integrate (dom, expr, fopt);
761 }
762 template <class Expr>
763 inline
764 typename
765 std::enable_if<
766  details::is_form_expr_v2_variational_arg<Expr>::value
767  ,details::form_lazy_terminal_integrate <details::form_expr_quadrature_on_element<Expr>>
768 >::type
771  const Expr& expr,
772  const integrate_option& fopt = integrate_option())
773 {
775  return lazy_integrate (expr_quad, fopt);
776 }
777 // ----------------------------------------------
778 // 6.3. subdomain by its name
779 // ----------------------------------------------
780 template <class Expr>
781 inline
782 typename
783 std::enable_if<
784  details::is_form_expr_quadrature_arg<Expr>::value
785  ,details::form_lazy_terminal_integrate <Expr>
786 >::type
789  const std::string& domname,
790  const Expr& expr,
791  const integrate_option& fopt = integrate_option())
792 {
795  dom = expr.get_trial_space().get_constitution().get_background_geo()[domname];
796  return lazy_integrate (dom, expr, fopt);
797 }
798 template <class Expr>
799 inline
800 typename
801 std::enable_if<
802  details::is_form_expr_v2_variational_arg<Expr>::value
803  ,details::form_lazy_terminal_integrate <details::form_expr_quadrature_on_element<Expr>>
804 >::type
807  const std::string& domname,
808  const Expr& expr,
809  const integrate_option& fopt = integrate_option())
810 {
812  return lazy_integrate (domname, expr_quad, fopt);
813 }
814 
815 }// namespace rheolef
816 #endif // _RHEO_INTEGRATE_H
field lh(Float epsilon, Float t, const test &v)
field gh(Float epsilon, Float t, const field &uh, const test &v)
typename scalar_traits< value_type >::type scalar_type
rheolef::std type
rheolef::std value
static iorheo::force_initialization dummy
Definition: iorheo.cc:147
#define trace_macro(message)
Definition: dis_macros.h:111
#define error_macro(message)
Definition: dis_macros.h:49
#define warning_macro(message)
Definition: dis_macros.h:53
Expr1::float_type T
Definition: field_expr.h:230
void field_expr_v2_value_assembly(const geo_basic< T, M > &omega, const Expr &expr0, const integrate_option &iopt, Result &result)
const std::string & valued_name(valued_type valued_tag)
This file is part of Rheolef.
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:211
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
Expr1::memory_type M
Definition: vec_expr_v2.h:416