Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
field_expr_variational_terminal.h
Go to the documentation of this file.
1#ifndef _RHEOLEF_FIELD_EXPR_VARIATIONAL_TERMINAL_H
2#define _RHEOLEF_FIELD_EXPR_VARIATIONAL_TERMINAL_H
23//
24// terminals variational expressions are used for field assembly
25// e.g. for right-hand-side of linear systems
26//
27// author: Pierre.Saramito@imag.fr
28//
29// date: 21 september 2015
30//
31// Notes: use template expressions and SFINAE techniques
32//
33// SUMMARY:
34// 1. concept
35// 2. grad, grad_s, D, etc
36// 3. div, div_s
37// 4. curl
38// 5. jump, average, inner, outer
39//
40#include "rheolef/field_expr.h"
41#include "rheolef/test.h"
42#include "rheolef/test_component.h"
43
44namespace rheolef {
45
46// -------------------------------------------------------------------
47// 1. concept
48// -------------------------------------------------------------------
49namespace details {
50
51// Define a trait type for detecting field expression valid arguments
52template<class T> struct is_field_expr_v2_variational_arg : std::false_type {};
53template<class T, class M, class VfTag> struct is_field_expr_v2_variational_arg<test_basic<T,M,VfTag> > : std::true_type {};
54template<class T, class M, class VfTag> struct is_field_expr_v2_variational_arg<test_component<T,M,VfTag> > : std::true_type {};
55
56} // namespace details
57// ---------------------------------------------------------------------------
58// 2. grad, grad_s, D, etc
59// ---------------------------------------------------------------------------
60namespace details {
61
62template<class Expr>
64public:
65// typedefs:
66
68 typedef typename Expr::memory_type memory_type;
75 typedef typename Expr::vf_tag_type vf_tag_type;
81
82// alocators:
83
85 : _expr(expr),
86 _gopt(gopt)
87 {
88 check_macro (gopt.broken || get_vf_space().get_basis().is_continuous(),
89 "grad(.): unexpected " << get_vf_space().get_basis().name()
90 << " discontinuous approximation (HINT: consider grad_h(.))");
91 }
96
97// accessors:
98
99 const space_type& get_vf_space() const { return _expr.get_vf_space(); }
102 space_constant::valued_type v = _expr.valued_tag();
103 switch (v) {
107 default:
108 fatal_macro ("unexpected " << space_constant::valued_name(v) << "-valued argument for grad() operator");
110 }
111 }
112 size_type n_derivative() const { return _expr.n_derivative() + 1; }
113
114// mutable modifiers:
115
118 const integrate_option& iopt)
119 {
120 _expr.initialize (pops, iopt);
121 }
125 const integrate_option& iopt)
126 {
127 _expr.initialize (gh, pops, iopt);
128 }
129 template<class Result>
130 void
133 const geo_element& K,
134 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value) const
135 {
136 _expr.template evaluate<Result,details::differentiate_option::gradient> (omega_K, K, _gopt, value);
137 }
138 template<class Result>
139 void
142 const geo_element& K,
143 const side_information_type& sid,
144 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value,
145 bool do_local_component_assembly) const
146 {
147 _expr.template evaluate_on_side<Result,details::differentiate_option::gradient> (omega_K, K, sid, _gopt, value, do_local_component_assembly);
148 }
149 template<class Value>
151 const geo_basic<float_type,memory_type>& omega_K,
152 const geo_element& S,
153 const geo_element& K0,
154 const geo_element& K1,
155 const Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value0,
156 const Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value1,
157 Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value) const
158 {
159 _expr.local_dg_merge_on_side (omega_K, S, K0, K1, value0, value1, value);
160 }
161 template<class Result>
162 void valued_check() const {
163 _expr.template grad_valued_check<Result>();
164 }
165protected:
166// data:
167 Expr _expr;
169};
170template<class Expr> struct is_field_expr_v2_variational_arg <field_expr_v2_variational_grad<Expr> > : std::true_type {};
171
172} // namespace details
173
174// grad(v)
175template<class Expr>
176inline
177typename
178std::enable_if<
181>::type
182grad (const Expr& expr)
183{
184 return details::field_expr_v2_variational_grad <Expr> (expr);
185}
186// grad_s(v)
187template<class Expr>
188inline
189typename
190std::enable_if<
191 details::is_field_expr_v2_variational_arg<Expr>::value
192 ,details::field_expr_v2_variational_grad<Expr>
193>::type
194grad_s (const Expr& expr)
195{
196 details::differentiate_option opt;
197 opt.surfacic = true;
198 return details::field_expr_v2_variational_grad <Expr> (expr, opt);
199}
200// grad_h(v)
201template<class Expr>
202inline
203typename
204std::enable_if<
205 details::is_field_expr_v2_variational_arg<Expr>::value
206 ,details::field_expr_v2_variational_grad<Expr>
207>::type
208grad_h (const Expr& expr)
209{
210 details::differentiate_option opt;
211 opt.broken = true;
212 return details::field_expr_v2_variational_grad <Expr> (expr, opt);
213}
214// D(v)
215template<class Expr>
216inline
217typename
218std::enable_if<
219 details::is_field_expr_v2_variational_arg<Expr>::value
220 ,details::field_expr_v2_variational_grad<Expr>
221>::type
222D (const Expr& expr)
223{
224 details::differentiate_option opt;
225 opt.symmetrized = true;
226 return details::field_expr_v2_variational_grad <Expr> (expr, opt);
227}
228// Ds(v)
229template<class Expr>
230inline
231typename
232std::enable_if<
233 details::is_field_expr_v2_variational_arg<Expr>::value
234 ,details::field_expr_v2_variational_grad<Expr>
235>::type
236Ds (const Expr& expr)
237{
238 details::differentiate_option opt;
239 opt.symmetrized = true;
240 opt.surfacic = true;
241 return details::field_expr_v2_variational_grad <Expr> (expr, opt);
242}
243// Dh(v)
244template<class Expr>
245inline
246typename
247std::enable_if<
248 details::is_field_expr_v2_variational_arg<Expr>::value
249 ,details::field_expr_v2_variational_grad<Expr>
250>::type
251Dh (const Expr& expr)
252{
253 details::differentiate_option opt;
254 opt.symmetrized = true;
255 opt.broken = true;
256 return details::field_expr_v2_variational_grad <Expr> (expr, opt);
257}
258// ---------------------------------------------------------------------------
259// 3. div, div_s
260// ---------------------------------------------------------------------------
261namespace details {
262
263template<class Expr>
265public:
266// typedefs:
267
269 typedef typename Expr::memory_type memory_type;
276 typedef typename Expr::vf_tag_type vf_tag_type;
282
283// alocators:
284
286 : _expr(expr),
287 _gopt(gopt)
288 {
289 check_macro (gopt.broken || get_vf_space().get_basis().is_continuous(),
290 "div(.): unexpected " << get_vf_space().get_basis().name()
291 << " discontinuous approximation (HINT: consider div_h(.))");
292 }
293
294// accessors:
295
296 const space_type& get_vf_space() const { return _expr.get_vf_space(); }
299 space_constant::valued_type v = _expr.valued_tag();
300 switch (v) {
304 default:
305 fatal_macro ("unexpected " << space_constant::valued_name(v) << "-valued argument for div() operator");
307 }
308 }
309 size_type n_derivative() const { return _expr.n_derivative() + 1; }
310
311// mutable modifiers:
312
314 _expr.initialize (pops, iopt);
315 }
317 _expr.initialize (gh, pops, iopt);
318 }
319 template<class Result>
320 void evaluate (
322 const geo_element& K,
323 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value) const
324 {
325 _expr.template evaluate<Result,details::differentiate_option::divergence> (omega_K, K, _gopt, value);
326 }
327 template<class Result>
330 const geo_element& K,
331 const side_information_type& sid,
332 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value,
333 bool do_local_component_assembly) const
334 {
335 _expr.template evaluate_on_side<Result,details::differentiate_option::divergence> (omega_K, K, sid, _gopt, value, do_local_component_assembly);
336 }
337 template<class Result>
338 void valued_check() const {
339 _expr.template div_valued_check<Result>();
340 }
341 template<class Value>
343 const geo_basic<float_type,memory_type>& omega_K,
344 const geo_element& S,
345 const geo_element& K0,
346 const geo_element& K1,
347 const Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value0,
348 const Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value1,
349 Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value) const
350 {
351 _expr.local_dg_merge_on_side (omega_K, S, K0, K1, value0, value1, value);
352 }
353protected:
354// data:
355 Expr _expr;
357};
358template<class Expr> struct is_field_expr_v2_variational_arg <field_expr_v2_variational_div<Expr> > : std::true_type {};
359
360} // namespace details
361
362// div(v)
363template<class Expr>
364inline
365typename
366std::enable_if<
369>::type
370div (const Expr& expr)
371{
372 return details::field_expr_v2_variational_div <Expr> (expr);
373}
374// div_s(v)
375template<class Expr>
376inline
377typename
378std::enable_if<
379 details::is_field_expr_v2_variational_arg<Expr>::value
380 ,details::field_expr_v2_variational_div<Expr>
381>::type
382div_s (const Expr& expr)
383{
384 details::differentiate_option opt;
385 opt.surfacic = true;
386 return details::field_expr_v2_variational_div <Expr> (expr, opt);
387}
388// div_h(v)
389template<class Expr>
390inline
391typename
392std::enable_if<
393 details::is_field_expr_v2_variational_arg<Expr>::value
394 ,details::field_expr_v2_variational_div<Expr>
395>::type
396div_h (const Expr& expr)
397{
398 details::differentiate_option opt;
399 opt.broken = true;
400 return details::field_expr_v2_variational_div <Expr> (expr, opt);
401}
402// ---------------------------------------------------------------------------
403// 4. curl
404// ---------------------------------------------------------------------------
405namespace details {
406
407template<class Expr>
409public:
410// typedefs:
411
413 typedef typename Expr::memory_type memory_type;
416 // value_type = vctor when d=2 and Expr is scalar or when d=3
417 // = scalar when d=2 and Expr is vector
418 // thus is undeterminated at compile-time
422 typedef typename Expr::vf_tag_type vf_tag_type;
428
429// alocators:
430
432 : _expr(expr),
433 _gopt(gopt)
434 {
435 check_macro (gopt.broken || get_vf_space().get_basis().is_continuous(),
436 "curl(.): unexpected " << get_vf_space().get_basis().name()
437 << " discontinuous approximation (HINT: consider curl_h(.))");
438 }
439
440// accessors:
441
442 const space_type& get_vf_space() const { return _expr.get_vf_space(); }
445 space_constant::valued_type arg_v = _expr.valued_tag();
446 switch (arg_v) {
449 size_type d = get_vf_space().get_geo().dimension();
451 }
452 default:
453 fatal_macro ("unexpected " << space_constant::valued_name(arg_v) << "-valued argument for curl() operator");
455 }
456 }
457 size_type n_derivative() const { return _expr.n_derivative() + 1; }
458
459// mutable modifiers:
460
462 _expr.initialize (pops, iopt);
463 }
465 _expr.initialize (gh, pops, iopt);
466 }
467 template<class Result>
468 void evaluate (
470 const geo_element& K,
471 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value) const
472 {
473 _expr.template evaluate<Result,details::differentiate_option::curl> (omega_K, K, _gopt, value);
474 }
475 template<class Value>
477 const geo_basic<float_type,memory_type>& omega_K,
478 const geo_element& S,
479 const geo_element& K0,
480 const geo_element& K1,
481 const Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value0,
482 const Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value1,
483 Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value) const
484 {
485 _expr.local_dg_merge_on_side (omega_K, S, K0, K1, value0, value1, value);
486 }
488 _expr.template valued_check<point_basic<scalar_type> >();
489 size_type d = get_vf_space().get_geo().dimension();
490 check_macro (d==2, "unexpected "<<d<<"D physical dimension for the scalar-valued curl() operator");
491 }
493 size_type d = get_vf_space().get_geo().dimension();
494 check_macro (d==2 || d==3, "unexpected "<<d<<"D physical dimension for the vector-valued curl() operator");
495 if (d == 2) {
496 _expr.template valued_check<scalar_type>();
497 } else {
498 _expr.template valued_check<point_basic<scalar_type> >();
499 }
500 }
502 fatal_macro ("unexpected tensor-valued result for the curl() operator");
503 }
504 template<class Result>
505 void valued_check() const {
506 _valued_check_internal(Result());
507 }
508protected:
509// data:
510 Expr _expr;
512};
513template<class Expr> struct is_field_expr_v2_variational_arg <field_expr_v2_variational_curl<Expr> > : std::true_type {};
514
515} // namespace details
516
517// curl(v)
518template<class Expr>
519inline
520typename
521std::enable_if<
524>::type
525curl (const Expr& expr)
526{
527 return details::field_expr_v2_variational_curl <Expr> (expr);
528}
529// bcurl(v) = Batchelor curl
530template<class Expr>
531inline
532typename
533std::enable_if<
534 details::is_field_expr_v2_variational_arg<Expr>::value
535 ,details::field_expr_v2_variational_curl<Expr>
536>::type
537bcurl (const Expr& expr)
538{
539 details::differentiate_option opt;
540 opt.batchelor_curl = true;
541 return details::field_expr_v2_variational_curl <Expr> (expr, opt);
542}
543// ---------------------------------------------------------------------------
544// 5. jump, average, inner, outer
545// ---------------------------------------------------------------------------
546// discontinuous Galerkin operators
547// in expressions templates for variationnal formulations
548//
549/*
550
551SPECIFICATION: a first example
552
553 Let v be a function defined over Omega
554 and discontinuous across internal sides (e.g. Pkd).
555 Let f be a function defined on Oemga.
556
557 We want to assembly
558 l(v) = int_{internal sides} f(x) [v](x) ds
559 where [v] is the jump of v across internal sides.
560
561 => l(v) = sum_{K is internal sides} int_K f(x) [v](x) ds
562
563 Let K be an internal side of the mesh of Omega.
564
565 int_K f(x) [v](x) ds
566 = int_{hat_K} f(F(hat_x))
567 [v](F(hat_x))
568 det(DF(hat_x)) d hat_s
569
570 where F is the piola transformation from the reference element hat_K to K:
571 F : hat_K ---> K
572 hat_x |--> x = F(hat_x)
573
574 The fonction v is not defined on a basis over internal sides K but over
575 elements L of the mesh of Omega.
576 Let L0 and L1 the two elements such that K is the common side of L0 and L1
577 and K is oriented from L0 to L1:
578 [v] = v0 - v1 on K, where v0=v/L0 and v1=v/L1.
579
580 Let G0 the piola transformation from the reference element tilde_L to L0:
581 G0 : tilde_L ---> L0
582 tilde_x |--> x = G0(tilde_x)
583
584 Conversely, let G1 the piola transformation from the reference element tilde_L to L1.
585
586 int_K f(x) [v](x) ds
587 = int_{hat_K} f(F(hat_x))
588 (v0-v1)(F(hat_x))
589 det(DF(hat_x)) d hat_s
590
591 The the basis fonction v0 and v1 are defined by using tilde_v, on the reference element tilde_L:
592 v0(x) = tilde_v (G0^{-1}(x))
593 v1(x) = tilde_v (G1^{-1}(x))
594 and with x=F(hat_x):
595 v0(F(hat_x)) = tilde_v (G0^{-1}(F(hat_x)))
596 v1(F(hat_x)) = tilde_v (G1^{-1}(F(hat_x)))
597
598 Thus:
599 int_K f(x) [v](x) ds
600 = int_{hat_K} f(F(hat_x))
601 ( tilde_v (G0^{-1}(F(hat_x)))
602 - tilde_v (G1^{-1}(F(hat_x))) )
603 det(DF(hat_x)) ds
604
605 Observe that H0=G0^{-1}oF is linear:
606 H0 : hat_K ---> tilde0_K subset tilde_L
607 hat_x ---> tilde0_x = H0(hat_x)
608 Conversely:
609 H1 : hat_K ---> tilde1_K subset tilde_L
610 hat_x ---> tilde1_x = H1(hat_x)
611 Thus, K linearly transforms by H0 into a side tilde0_K of the reference element tilde_L
612 and, by H1, into another side tilde1_K of tilde_L.
613
614 int_K f(x) [v](x) ds
615 = int_{hat_K} f(F(hat_x))
616 ( tilde_v (H0(hat_x))
617 - tilde_v (H1(hat_x)) )
618 det(DF(hat_x)) ds
619
620 Let (hat_xq, hat_wq)_{q=0...} a quadrature formulae over hat_K.
621 The integral becomes:
622
623 int_K f(x) [v](x) ds
624 = sum_q f(F(hat_xq))
625 ( tilde_v (H0(hat_xq))
626 - tilde_v (H1(hat_xq)) )
627 det(DF(hat_xq)) hat_wq
628
629 Then, the basis functions tilde_v can be computed one time for all
630 over all the sides tilde(i)_K of the reference element tilde_L, i=0..nsides(tilde_L)
631 at the quadratures nodes tilde(i)_xq = Hi(hat_xq):
632 tilde_v (Hi(hat_xq)), i=0..nsides(tilde_L), q=0..nq(hat_K)
633
634SPECIFICATION: a second example
635
636 We want to assembly
637 l(v) = int_{internal sides} f(x) [grad(v).n](x) ds
638 where [grad(v).n] is the jump of the normal derivative of v across internal sides.
639
640 Let K be an internal side of the mesh of Omega.
641
642 int_K f(x) [grad(v).n](x) ds
643 = int_{hat_K} f(F(hat_x))
644 [grad(v).n](F(hat_x))
645 det(DF(hat_x)) d hat_s
646
647 = int_{hat_K} f(F(hat_x))
648 (grad(v0).n)(F(hat_x))
649 det(DF(hat_x)) d hat_s
650 - int_{hat_K} f(F(hat_x))
651 (grad(v1).n)(F(hat_x))
652 det(DF(hat_x)) d hat_s
653
654 where v0=v/L0 and v1=v/L1 and Li are the two elements containing the side K.
655 Let us fix one of the Li and omits the i subscript.
656 The computation reduces to evaluate:
657
658 int_K f(x) grad(v).n(x) ds
659 = sum_q f(F(hat_xq))
660 (grad(v).n)(F(hat_xq))
661 det(DF(hat_xq)) hat_wq
662
663 From the gradient transformation:
664
665 grad(v)(F(hat_xq)) = DG^{-T}(H(hat_xq)) * tilde_grad(tilde_v)(H(hat_xq))
666
667 where H = G^{-1}oF is linear from hat_K to tilde_K subset tilde_L.
668
669 int_K f(x) grad(v).n(x) ds
670 = sum_q f(F(hat_xq))
671 DG^{-T}(H(hat_xq))*tilde_grad(tilde_v)(H(hat_xq))
672 .n(F(hat_xq))
673 det(DF(hat_xq)) hat_wq
674
675 We can evaluate one time for all the gradients of basis functions tilde_v
676 on the quadrature nodes of each sides tilde_K of tilde_L :
677 tilde_grad(tilde_v)(H(hat_xq))
678 The piola basis functions and their derivatives are also evaluated one time for all on these nodes :
679 DG^{-T}(H(hat_xq))
680
681 The normal vector
682 n(xq), xq=F(hat_xq), q=...
683 should be evaluated on K, not on L that has no normal vector.
684
685IMPLEMENTATION: bassis evaluation => test.cc
686
687 The basis_on_pointset class extends to the case of an integration over a side of
688
689 test_rep<T,M>::initialize (const geo_basic<float_type,M>& dom, const quadrature<T>& quad, const integrate_option& iopt) const {
690 _basis_on_quad.set (quad, get_vf_space().get_basis());
691 _piola_on_quad.set (quad, get_vf_space().get_geo().get_piola_basis());
692 => inchange'
693 }
694 test_rep<T,M>::element_initialize (const geo_element& L, size_type loc_isid=-1) const {
695 if (loc_isid != -1) {
696 basis_on_quad.restrict_on_side (tilde_L, loc_isid);
697 piola_on_quad.restrict_on_side (tilde_L, loc_isid);
698 }
699 }
700 test_rep<T,M>::basis_evaluate (...) {
701 // Then, a subsequent call to
702 basis_on_quad.evaluate (tilde_L, q);
703 // will restrict to the loc_isid part.
704 }
705
706IMPLEMENTATION: normal vector => field_vf_expr.h & field_nl_expr_terminal.h
707 on propage des vf_expr aux nl_expr le fait qu'on travaille sur une face :
708 class nl_helper {
709 void element_initialize (const This& obj, const geo_element& L, size_type loc_isid=-1) const {
710 obj._nl_expr.evaluate (L, isid, obj._vector_nl_value_quad);
711 }
712 };
713 pour la classe normal :
714 field_expr_terminal_normal::evaluate (L, loc_isid, value) {
715 if (loc_isid != -1) K=side(L,loc_isid); else K=L;
716 puis inchange.
717 }
718 pour la classe terminal_field: si on evalue un field uh qui est discontinu :
719 on sait sur quelle face il se restreint :
720 field_expr_terminal_field::evaluate (L, loc_isid, value) {
721 if (loc_isid != -1) {
722 _basis_on_quad.restrict_on_side (tilde_L, loc_isid);
723 }
724 for (q..) {
725 general_field_evaluate (_uh, _basis_on_quad, tilde_L, _dis_idof, q, value[q]);
726 }
727 }
728IMPLEMENTATION: bassis evaluation => basis_on_pointset.cc
729 c'est la que se fait le coeur du travail :
730 basis_on_pointset::restrict_on_side (tilde_L, loc_isid)
731 => initialise
732
733 a l'initialisation, on evalue une fois pour tte
734 sur toutes les faces en transformant la quadrature via
735 tilde(i)_xq = Hi(hat_xq)
736 tilde(i)_wq = ci*hat_wq
737 avec
738 ci = meas(tilde(i)_K)/meas(hat_K)
739
740 puis :
741 basis_on_pointset::evaluate (tilde_L, q)
742 on se baladera dans la tranche [loc_isid*nq, (loc_isid+1)*nq[
743 du coup, on positionne un pointeur de debut q_start = loc_isid*nq
744 et une taille q_size = nq
745 si les faces sont differentes (tri,qua) dans un prisme, il faudra
746 un tableau de pointeurs pour gerer cela :
747 q_start [loc_nsid+1]
748 q_size [loc_isid] = q_start[loc_isid+1] - q_start[loc_isid]
749
750 basis_on_pointset::begin() { return _val[_curr_K_variant][_curr_q].begin() + q_start[_curr_K_variant][loc_isid]; }
751 basis_on_pointset::begin() { return _val[_curr_K_variant][_curr_q].begin() + q_start[_curr_K_variant][loc_isid+1]; }
752
753 et le tour est joue' !
754
755PLAN DE DEVELOPPEMENT:
756 1) DG transport
757 basis_on_pointset.cc
758 test.cc
759 essais :
760 lh = integrate(jump(v)*f);
761 convect_dg2.cc
762 2) DG diffusion : avec normale et gradient
763 field_vf_expr.h
764 class nl_helper
765 field_nl_expr_terminal.h
766 field_expr_terminal_normal::evaluate (L, loc_isid, value)
767 field_expr_terminal_field ::evaluate (L, loc_isid, value)
768
769*/
770
771namespace details {
772// ---------------------------------------------------------------------------
773// 5.1 class dg
774// ---------------------------------------------------------------------------
775template<class Expr>
777public:
778// typedefs:
779
780 using size_type = typename Expr::size_type;
781 using memory_type = typename Expr::memory_type;
782 using value_type = typename Expr::value_type;
783 using scalar_type = typename Expr::scalar_type;
784 using float_type = typename Expr::float_type;
785 using space_type = typename Expr::space_type;
786 using vf_tag_type = typename Expr::vf_tag_type;
790 using is_elementwise = std::false_type; // assembles two adjacents elements
791
792// alocators:
793
794 field_expr_v2_variational_dg (const Expr& expr, const float_type& c0, const float_type& c1);
795
796// accessors:
797
798 const space_type& get_vf_space() const { return _expr0.get_vf_space(); }
799 static const space_constant::valued_type valued_hint = Expr::valued_hint;
800 space_constant::valued_type valued_tag() const { return _expr0.valued_tag(); }
801 size_type n_derivative() const { return _expr0.n_derivative(); }
802
804 _expr0.initialize (pops, iopt);
805 _expr1.initialize (pops, iopt);
806 _bgd_omega = get_vf_space().get_constitution().get_background_geo();
807 }
809 _expr0.initialize (gh, pops, iopt);
810 _expr1.initialize (gh, pops, iopt);
811 fatal_macro ("unsupported discontinuous Galerkin on a band"); // how to define background mesh _bgd_omega ?
812 }
813 template<class Result>
814 void evaluate (
816 const geo_element& K,
817 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value) const;
818
819 template<class Result>
820 void valued_check() const {
821 check_macro (get_vf_space().get_constitution().is_discontinuous(),
822 "unexpected continuous test-function in space " << get_vf_space().name()
823 << " for jump or average operator (HINT: omit jump or average)");
824 _expr0.template valued_check<Result>();
825 }
826protected:
827// data:
828 mutable Expr _expr0;
829 mutable Expr _expr1;
835};
836// concept:
837template<class Expr> struct is_field_expr_v2_variational_arg <field_expr_v2_variational_dg<Expr> > : std::true_type {};
838
839// allocator:
840template<class Expr>
841inline
843 const Expr& expr,
844 const float_type& c0,
845 const float_type& c1)
846: _expr0(expr),
847 _expr1(expr),
848 _c0(c0),
849 _c1(c1),
850 _tilde0_L0(),
851 _tilde1_L1(),
852 _bgd_omega()
853{
854}
855// ---------------------------------------------------------------------------
856// 5.2. evaluate
857// ---------------------------------------------------------------------------
858template<class Expr>
859template<class Result>
860void
863 const geo_element& K,
864 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value) const
865{
866trace_macro ("evaluate_dg: "<<K.name()<<K.dis_ie()<<"...")
867 check_macro (_bgd_omega == omega_K.get_background_geo().get_background_geo(),
868 "discontinuous Galerkin: incompatible integration domain "<<omega_K.name() << " and test function in "
869 << get_vf_space().name());
870 size_type L_map_d = K.dimension() + 1;
871 size_type L_dis_ie0 = K.master(0);
872 size_type L_dis_ie1 = K.master(1);
873 check_macro (L_dis_ie0 != std::numeric_limits<size_type>::max(),
874 "unexpected isolated mesh side "<<K.name()<<K.dis_ie()<<" dim="<<K.dimension());
875 const geo_element& L0 = _bgd_omega.dis_get_geo_element (L_map_d, L_dis_ie0);
877 L0.get_side_informations (K, sid0);
878 if (L_dis_ie1 == std::numeric_limits<size_type>::max()) {
879 // K is a boundary side
880 bool do_local_component_assembly = true;
881 // average(v)=jump(v)=inner(v)=outer(v)=v on the boundary
882 _expr0.evaluate_on_side (omega_K, L0, sid0, value, do_local_component_assembly);
883 } else {
884 // K is an internal side
885 const geo_element& L1 = _bgd_omega.dis_get_geo_element (L_map_d, L_dis_ie1);
887 L1.get_side_informations (K, sid1);
888 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic> value0, value1;
889 bool do_local_component_assembly = false;
890 _expr0.evaluate_on_side (omega_K, L0, sid0, value0, do_local_component_assembly);
891 _expr1.evaluate_on_side (omega_K, L1, sid1, value1, do_local_component_assembly);
892 for (size_type loc_inod = 0, loc_nnod = value0.rows(); loc_inod < loc_nnod; ++loc_inod) {
893 for (size_type loc_idof = 0, loc_ndof = value0.cols(); loc_idof < loc_ndof; ++loc_idof) {
894 value0 (loc_inod,loc_idof) = _c0*value0 (loc_inod,loc_idof); // TODO: DVT_EIGEN_BLAS2
895 }}
896 for (size_type loc_inod = 0, loc_nnod = value1.rows(); loc_inod < loc_nnod; ++loc_inod) {
897 for (size_type loc_idof = 0, loc_ndof = value1.cols(); loc_idof < loc_ndof; ++loc_idof) {
898 value1 (loc_inod,loc_idof) = _c1*value1 (loc_inod,loc_idof); // TODO: DVT_EIGEN_BLAS2
899 }}
900 // merge loc_dofs at the component level
901 // assume that space is not hierarchical, otherwise see dg_merge
902 _expr0.local_dg_merge_on_side (omega_K, K, L0, L1, value0, value1, value);
903 }
904trace_macro ("evaluate_dg: "<<K.name()<<K.dis_ie()<<" done")
905}
906
907} // namespace details
908// ---------------------------------------------------------------------------
909// 5.4. user-level interface
910// ---------------------------------------------------------------------------
911#define _RHEOLEF_make_field_expr_v2_variational_dg(op,c0,c1) \
912template<class Expr> \
913inline \
914typename \
915std::enable_if< \
916 details::is_field_expr_v2_variational_arg<Expr>::value \
917 ,details::field_expr_v2_variational_dg<Expr> \
918>::type \
919op (const Expr& expr) \
920{ \
921 return details::field_expr_v2_variational_dg <Expr> (expr, c0, c1); \
922}
923
928#undef _RHEOLEF_make_field_expr_v2_variational_dg
929
930} // namespace rheolef
931#endif // _RHEOLEF_FIELD_EXPR_VARIATIONAL_TERMINAL_H
field gh(Float epsilon, Float t, const field &uh, const test &v)
field_expr_v2_variational_curl(const Expr &expr, const differentiate_option &gopt=differentiate_option())
void _valued_check_internal(const tensor_basic< scalar_type > &) const
void initialize(const piola_on_pointset< float_type > &pops, const integrate_option &iopt)
void initialize(const band_basic< float_type, memory_type > &gh, const piola_on_pointset< float_type > &pops, const integrate_option &iopt)
scalar_traits< typenameExpr::value_type >::type scalar_type
field_expr_v2_variational_curl< typename Expr::dual_self_type > dual_self_type
void _valued_check_internal(const point_basic< scalar_type > &) const
details::dual_vf_tag< vf_tag_type >::type vf_dual_tag_type
void local_dg_merge_on_side(const geo_basic< float_type, memory_type > &omega_K, const geo_element &S, const geo_element &K0, const geo_element &K1, const Eigen::Matrix< Value, Eigen::Dynamic, Eigen::Dynamic > &value0, const Eigen::Matrix< Value, Eigen::Dynamic, Eigen::Dynamic > &value1, Eigen::Matrix< Value, Eigen::Dynamic, Eigen::Dynamic > &value) const
void evaluate(const geo_basic< float_type, memory_type > &omega_K, const geo_element &K, Eigen::Matrix< Result, Eigen::Dynamic, Eigen::Dynamic > &value) const
void initialize(const piola_on_pointset< float_type > &pops, const integrate_option &iopt)
typename details::dual_vf_tag< vf_tag_type >::type vf_dual_tag_type
field_expr_v2_variational_dg(const Expr &expr, const float_type &c0, const float_type &c1)
void initialize(const band_basic< float_type, memory_type > &gh, const piola_on_pointset< float_type > &pops, const integrate_option &iopt)
field_expr_v2_variational_dg< typename Expr::dual_self_type > dual_self_type
void evaluate(const geo_basic< float_type, memory_type > &omega_K, const geo_element &K, Eigen::Matrix< Result, Eigen::Dynamic, Eigen::Dynamic > &value) const
void initialize(const piola_on_pointset< float_type > &pops, const integrate_option &iopt)
space_constant::rank_down< typenameExpr::value_type >::type value_type
void initialize(const band_basic< float_type, memory_type > &gh, const piola_on_pointset< float_type > &pops, const integrate_option &iopt)
scalar_traits< typenameExpr::value_type >::type scalar_type
field_expr_v2_variational_div(const Expr &expr, const differentiate_option &gopt=differentiate_option())
details::dual_vf_tag< vf_tag_type >::type vf_dual_tag_type
void local_dg_merge_on_side(const geo_basic< float_type, memory_type > &omega_K, const geo_element &S, const geo_element &K0, const geo_element &K1, const Eigen::Matrix< Value, Eigen::Dynamic, Eigen::Dynamic > &value0, const Eigen::Matrix< Value, Eigen::Dynamic, Eigen::Dynamic > &value1, Eigen::Matrix< Value, Eigen::Dynamic, Eigen::Dynamic > &value) const
void evaluate_on_side(const geo_basic< float_type, memory_type > &omega_K, const geo_element &K, const side_information_type &sid, Eigen::Matrix< Result, Eigen::Dynamic, Eigen::Dynamic > &value, bool do_local_component_assembly) const
void evaluate(const geo_basic< float_type, memory_type > &omega_K, const geo_element &K, Eigen::Matrix< Result, Eigen::Dynamic, Eigen::Dynamic > &value) const
field_expr_v2_variational_div< typename Expr::dual_self_type > dual_self_type
field_expr_v2_variational_grad(const field_expr_v2_variational_grad< Expr > &x)
void initialize(const piola_on_pointset< float_type > &pops, const integrate_option &iopt)
field_expr_v2_variational_grad(const Expr &expr, const differentiate_option &gopt=differentiate_option())
field_expr_v2_variational_grad< typename Expr::dual_self_type > dual_self_type
void initialize(const band_basic< float_type, memory_type > &gh, const piola_on_pointset< float_type > &pops, const integrate_option &iopt)
scalar_traits< typenameExpr::value_type >::type scalar_type
details::dual_vf_tag< vf_tag_type >::type vf_dual_tag_type
void local_dg_merge_on_side(const geo_basic< float_type, memory_type > &omega_K, const geo_element &S, const geo_element &K0, const geo_element &K1, const Eigen::Matrix< Value, Eigen::Dynamic, Eigen::Dynamic > &value0, const Eigen::Matrix< Value, Eigen::Dynamic, Eigen::Dynamic > &value1, Eigen::Matrix< Value, Eigen::Dynamic, Eigen::Dynamic > &value) const
space_constant::rank_up< typenameExpr::value_type >::type value_type
void evaluate_on_side(const geo_basic< float_type, memory_type > &omega_K, const geo_element &K, const side_information_type &sid, Eigen::Matrix< Result, Eigen::Dynamic, Eigen::Dynamic > &value, bool do_local_component_assembly) const
void evaluate(const geo_basic< float_type, memory_type > &omega_K, const geo_element &K, Eigen::Matrix< Result, Eigen::Dynamic, Eigen::Dynamic > &value) const
generic mesh with rerefence counting
Definition geo.h:1089
see the geo_element page for the full documentation
reference_element::size_type size_type
size_type dimension() const
size_type master(bool i) const
orientation_type get_side_informations(const geo_element &S, size_type &loc_isid, size_type &shift) const
size_type dis_ie() const
see the integrate_option page for the full documentation
see the reference_element page for the full documentation
the finite element space
Definition space.h:382
#define trace_macro(message)
Definition dis_macros.h:111
#define fatal_macro(message)
Definition dis_macros.h:33
Expr1::float_type T
Definition field_expr.h:230
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
#define _RHEOLEF_make_field_expr_v2_variational_dg(op, c0, c1)
const std::string & valued_name(valued_type valued_tag)
This file is part of Rheolef.
std::enable_if< details::has_field_rdof_interface< Expr >::value, details::field_expr_v2_nonlinear_terminal_field< typenameExpr::scalar_type, typenameExpr::memory_type, details::differentiate_option::gradient > >::type grad_s(const Expr &expr)
grad_s(uh): see the expression page for the full documentation
std::enable_if< details::has_field_rdof_interface< Expr >::value, details::field_expr_v2_nonlinear_terminal_field< typenameExpr::scalar_type, typenameExpr::memory_type, details::differentiate_option::gradient > >::type grad(const Expr &expr)
grad(uh): see the expression page for the full documentation
std::enable_if< details::has_field_rdof_interface< Expr >::value, details::field_expr_v2_nonlinear_terminal_field< typenameExpr::scalar_type, typenameExpr::memory_type, details::differentiate_option::gradient > >::type D(const Expr &expr)
D(uh): see the expression page for the full documentation.
std::enable_if< details::has_field_rdof_interface< Expr >::value, details::field_expr_v2_nonlinear_terminal_field< typenameExpr::scalar_type, typenameExpr::memory_type, details::differentiate_option::divergence > >::type div(const Expr &expr)
div(uh): see the expression page for the full documentation
std::enable_if< details::has_field_rdof_interface< Expr >::value, details::field_expr_v2_nonlinear_terminal_field< typenameExpr::scalar_type, typenameExpr::memory_type, details::differentiate_option::divergence > >::type div_h(const Expr &expr)
div_h(uh): see the expression page for the full documentation
std::enable_if< details::is_field_expr_v2_variational_arg< Expr >::value, details::field_expr_v2_variational_curl< Expr > >::type bcurl(const Expr &expr)
std::enable_if< details::has_field_rdof_interface< Expr >::value, details::field_expr_v2_nonlinear_terminal_field< typenameExpr::scalar_type, typenameExpr::memory_type, details::differentiate_option::curl > >::type curl(const Expr &expr)
curl(uh): see the expression page for the full documentation
std::enable_if< details::has_field_rdof_interface< Expr >::value, details::field_expr_v2_nonlinear_terminal_field< typenameExpr::scalar_type, typenameExpr::memory_type, details::differentiate_option::gradient > >::type grad_h(const Expr &expr)
grad_h(uh): see the expression page for the full documentation
std::enable_if< details::has_field_rdof_interface< Expr >::value, details::field_expr_v2_nonlinear_terminal_field< typenameExpr::scalar_type, typenameExpr::memory_type, details::differentiate_option::gradient > >::type Dh(const Expr &expr)
Dh(uh): see the expression page for the full documentation.
std::enable_if< details::has_field_rdof_interface< Expr >::value, details::field_expr_v2_nonlinear_terminal_field< typenameExpr::scalar_type, typenameExpr::memory_type, details::differentiate_option::divergence > >::type div_s(const Expr &expr)
div_s(uh): see the expression page for the full documentation
std::enable_if< details::has_field_rdof_interface< Expr >::value, details::field_expr_v2_nonlinear_terminal_field< typenameExpr::scalar_type, typenameExpr::memory_type, details::differentiate_option::gradient > >::type Ds(const Expr &expr)
Ds(uh): see the expression page for the full documentation.
undeterminated_basic< typename scalar_traits< T >::type > type
point_basic< typename scalar_traits< T >::type > type
helper for generic field value_type: T, point_basic<T> or tensor_basic<T>
Expr1::memory_type M