Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
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
24namespace 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
195namespace rheolef {
196
197// ---------------------------------------------------
198// 1. scalar-result integration
199// ---------------------------------------------------
200// 1.1. general integration of a nonlinear expression
201// ---------------------------------------------------
202template <class T, class M, class Expr,
203 class Result = typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr>::type::value_type>
204inline
205typename std::enable_if<
206 details::is_field_expr_v2_nonlinear_arg<Expr>::value
207 && ! is_undeterminated<Result>::value,
208 Result
209>::type
211integrate (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// ---------------------------------------------------
225template <class T, class M>
226T
229{
230 if (iopt.get_order() == std::numeric_limits<integrate_option::size_type>::max()) {
231 iopt.set_order(0);
232 }
233 details::f_constant <point_basic<T>,T> one(1);
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
240template<class T, class M, class Expr>
241inline
242typename 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
248integrate (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;
252 switch (expr.valued_tag()) {
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// -------------------------------------------------------
271template <class T, class M, class Expr>
272inline
273typename
274std::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}
288template <class T, class M, class Expr>
289inline
290typename
291std::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// ----------------------------------------------
307template <class Expr>
308inline
309typename
310std::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{
319 field_basic <typename Expr::scalar_type, typename Expr::memory_type> lh;
320 const geo_basic <typename Expr::scalar_type, typename Expr::memory_type>&
321 dom = expr.get_vf_space().get_constitution().get_geo();
322 lh.do_integrate (dom, expr, iopt);
323 return lh;
324}
325template <class Expr>
326inline
327typename
328std::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// ----------------------------------------------
343template <class Expr>
344inline
345typename
346std::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{
356 field_basic <typename Expr::scalar_type, typename Expr::memory_type> lh;
357 const geo_basic <typename Expr::scalar_type, typename Expr::memory_type>&
358 dom = expr.get_vf_space().get_constitution().get_geo() [domname];
359 lh.do_integrate (dom, expr, iopt);
360 return lh;
361}
362template <class Expr>
363inline
364typename
365std::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// -------------------------------------------------------
383template <class T, class M, class Expr>
384inline
385typename
386std::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}
400template <class T, class M, class Expr>
401inline
402typename
403std::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// ----------------------------------------------
419template <class Expr>
420inline
421typename
422std::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{
431 form_basic <typename Expr::scalar_type, typename Expr::memory_type> a;
432 geo_basic <typename Expr::scalar_type, typename Expr::memory_type>
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}
456template <class Expr>
457inline
458typename
459std::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// ----------------------------------------------
474template <class Expr>
475inline
476typename
477std::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{
487 form_basic <typename Expr::scalar_type, typename Expr::memory_type> a;
488 const geo_basic <typename Expr::scalar_type, typename Expr::memory_type>&
489 dom = expr.get_trial_space().get_constitution().get_background_geo()[domname];
490 a.do_integrate (dom, expr, fopt);
491 return a;
492}
493template <class Expr>
494inline
495typename
496std::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// ----------------------------------------------
512template <class T, class M, class Expr>
513inline
514typename
515std::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{
525 field_basic <typename Expr::scalar_type, typename Expr::memory_type> lh;
526 lh.do_integrate (gh, expr, iopt);
527 return lh;
528}
529template <class T, class M, class Expr>
530inline
531typename
532std::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}
546template <class T, class M, class Expr>
547inline
548typename
549std::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{
559 form_basic <typename Expr::scalar_type, typename Expr::memory_type> a;
560 a.do_integrate (gh, expr, fopt);
561 return a;
562}
563template <class T, class M, class Expr>
564inline
565typename
566std::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// -------------------------------------------------------
585template <class Expr>
586inline
587typename
588std::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}
600template <class Expr>
601inline
602typename
603std::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// ----------------------------------------------
619template <class Expr>
620inline
621typename
622std::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}
633template <class Expr>
634inline
635typename
636std::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// ----------------------------------------------
651template <class Expr>
652inline
653typename
654std::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}
666template <class Expr>
667inline
668typename
669std::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// -------------------------------------------------------
688template <class Expr>
689inline
690typename
691std::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{
701 return details::form_lazy_terminal_integrate<Expr> (domain, expr, iopt);
702}
703template <class Expr>
704inline
705typename
706std::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// ----------------------------------------------
722template <class Expr>
723inline
724typename
725std::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;
737 geo_basic <typename Expr::scalar_type, typename Expr::memory_type>
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}
762template <class Expr>
763inline
764typename
765std::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// ----------------------------------------------
780template <class Expr>
781inline
782typename
783std::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{
793 form_basic <typename Expr::scalar_type, typename Expr::memory_type> a;
794 const geo_basic <typename Expr::scalar_type, typename Expr::memory_type>&
795 dom = expr.get_trial_space().get_constitution().get_background_geo()[domname];
796 return lazy_integrate (dom, expr, fopt);
797}
798template <class Expr>
799inline
800typename
801std::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)
generic mesh with rerefence counting
Definition geo.h:1089
see the integrate_option page for the full documentation
typename scalar_traits< value_type >::type scalar_type
#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
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