Rheolef  7.2
an efficient C++ finite element environment
field_expr_variational.h
Go to the documentation of this file.
1#ifndef _RHEOLEF_FIELD_EXPR_VARIATIONAL_H
2#define _RHEOLEF_FIELD_EXPR_VARIATIONAL_H
23//
24// 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. unary function
35// 1.1. unary node
36// 1.2. unary calls
37// 2. binary operators +- between two variational exprs
38// 2.1. binary node
39// 2.2. binary calls
40// 3. binary operators */ between a variational and a nonlinear expr
41// 3.1. helper
42// 3.2. binary node
43// 3.3. binary calls
44//
45#include "rheolef/field_expr_variational_terminal.h"
46
47namespace rheolef {
48
49// ---------------------------------------------------------------------------
50// 1. unary function
51// ---------------------------------------------------------------------------
52// 1.1. unary node
53// ---------------------------------------------------------------------------
54// ex: -v, 2*v, v/3
55namespace details {
56
57template <class UnaryFunction, class Expr>
59public:
60// typedefs:
61
63 typedef typename Expr::memory_type memory_type;
64 typedef typename details::generic_unary_traits<UnaryFunction>::template result_hint<typename Expr::value_type>::type
69 typedef typename Expr::vf_tag_type vf_tag_type;
75
76// alocators:
77
78 field_expr_v2_variational_unary (const UnaryFunction& f, const Expr& expr)
79 : _f(f), _expr(expr) {}
80
82 : _f(x._f), _expr(x._expr) {}
83
84// accessors:
85
86 static bool have_test_space() { return true; } // check !
87 const space_type& get_vf_space() const { return _expr.get_vf_space(); }
91 }
92 size_type n_derivative() const { return _expr.n_derivative(); }
93
94// mutable modifiers:
95
97 _expr.initialize (pops, iopt);
98 }
100 _expr.initialize (gh, pops, iopt);
101 }
102 // -------------
103 // evaluate
104 // -------------
105 // evaluate when all arg types are determinated
106 template<class Result, class Arg, class Status>
108 template<class M>
110 const self_type& obj,
111 const geo_basic<float_type,M>& omega_K,
112 const geo_element& K, Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value) const
113 {
114 fatal_macro ("invalid type resolution: Result="<<typename_macro(Result)
115 << ", Arg="<<typename_macro(Arg)
116 << ", UnaryFunction="<<typename_macro(UnaryFunction)
117 );
118 }
119 };
120 template<class Result, class Arg>
121 struct evaluate_call_check<Result,Arg,std::true_type> {
122 template<class M>
124 const self_type& obj,
125 const geo_basic<float_type,M>& omega_K,
126 const geo_element& K,
127 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value) const
128 {
129 Eigen::Matrix<Arg,Eigen::Dynamic,Eigen::Dynamic> arg_value;
130 obj._expr.evaluate (omega_K, K, arg_value);
131 value = arg_value.unaryExpr (obj._f);
132 }
133 template<class M>
135 const self_type& obj,
136 const geo_basic<float_type,M>& omega_K,
137 const geo_element& K,
138 const side_information_type& sid,
139 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value,
140 bool do_local_component_assembly) const
141 {
142 Eigen::Matrix<Arg,Eigen::Dynamic,Eigen::Dynamic> arg_value;
143 obj._expr.evaluate_on_side (omega_K, K, sid, arg_value, do_local_component_assembly);
144 value = arg_value.unaryExpr (obj._f);
145 }
146 };
147 template<class Result, class Arg, class M>
149 const geo_basic<float_type,M>& omega_K,
150 const geo_element& K,
151 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value) const
152 {
153 // check if Result is a valid return_type for this function:
154 typedef typename scalar_traits<Result>::type S;
155 typedef undeterminated_basic<S> undet;
156 typedef typename details::generic_unary_traits<UnaryFunction>::template hint<Arg,undet>::result_type result_type;
157 // TODO: instead of is_equal, could have compatible scalars T1,T2 ?
158 typedef typename details::is_equal<Result,result_type>::type status_t;
160 eval (*this, omega_K, K, value);
161 }
162 template<class Result, class Arg, class M>
164 const geo_basic<float_type,M>& omega_K,
165 const geo_element& K,
166 const side_information_type& sid,
167 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value,
168 bool do_local_component_assembly) const
169 {
170 // check if Result is a valid return_type for this function:
171 typedef typename scalar_traits<Result>::type S;
172 typedef undeterminated_basic<S> undet;
173 typedef typename details::generic_unary_traits<UnaryFunction>::template hint<Arg,undet>::result_type result_type;
174 // TODO: instead of is_equal, could have compatible scalars T1,T2 ?
175 typedef typename details::is_equal<Result,result_type>::type status_t;
177 eval (*this, omega_K, K, sid, value, do_local_component_assembly);
178 }
179 // when Arg is defined at compile time:
180 template<class This, class Result, class Arg, class Status>
182 template<class M>
184 const This& obj,
185 const geo_basic<float_type,M>& omega_K,
186 const geo_element& K,
187 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value) const
188 {
189 obj.template evaluate_internal<Result, Arg> (omega_K, K, value);
190 }
191 template<class M>
193 const This& obj,
194 const geo_basic<float_type,M>& omega_K,
195 const geo_element& K,
196 const side_information_type& sid,
197 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value,
198 bool do_local_component_assembly) const
199 {
200 obj.template evaluate_internal<Result, Arg> (omega_K, K, sid, value, do_local_component_assembly);
201 }
202 };
203 // when Arg is undeterminated at compile time
204 template<class This, class Result, class Arg>
205 struct evaluate_switch<This, Result, Arg, std::true_type> {
206 template<class M>
208 const This& obj,
209 const geo_basic<float_type,M>& omega_K,
210 const geo_element& K,
211 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value) const
212 {
213 typedef typename scalar_traits<Arg>::type T;
214 space_constant::valued_type arg_valued_tag = obj._expr.valued_tag();
215 switch (arg_valued_tag) {
216#define _RHEOLEF_switch(VALUED,VALUE) \
217 case space_constant::VALUED: \
218 obj.template evaluate_internal<Result, VALUE>(omega_K, K, value); break;
225#undef _RHEOLEF_switch
226 default: error_macro ("unexpected argument valued tag="<<arg_valued_tag);
227 }
228 }
229 template<class M>
231 const This& obj,
232 const geo_basic<float_type,M>& omega_K,
233 const geo_element& K,
234 const side_information_type& sid,
235 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value,
236 bool do_local_component_assembly) const
237 {
238 typedef typename scalar_traits<Arg>::type T;
239 space_constant::valued_type arg_valued_tag = obj._expr.valued_tag();
240 switch (arg_valued_tag) {
241#define _RHEOLEF_switch(VALUED,VALUE) \
242 case space_constant::VALUED: \
243 obj.template evaluate_internal<Result, VALUE>(omega_K, K, sid,value, do_local_component_assembly); break;
250#undef _RHEOLEF_switch
251 default: error_macro ("unexpected argument valued tag="<<arg_valued_tag);
252 }
253 }
254 };
255 template<class Result, class M>
256 void evaluate (
257 const geo_basic<float_type,M>& omega_K,
258 const geo_element& K,
259 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value) const
260 {
261 typedef typename details::generic_unary_traits<UnaryFunction>::template hint<typename Expr::value_type,Result>::argument_type
262 A1;
263 typedef typename is_undeterminated<A1>::type status_t;
265 eval (*this, omega_K, K, value);
266 }
267 template<class Result, class M>
269 const geo_basic<float_type,M>& omega_K,
270 const geo_element& K,
271 const side_information_type& sid,
272 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value,
273 bool do_local_component_assembly) const
274 {
275 typedef typename details::generic_unary_traits<UnaryFunction>::template hint<typename Expr::value_type,Result>::argument_type
276 A1;
277 typedef typename is_undeterminated<A1>::type status_t;
279 eval (*this, omega_K, K, sid, value, do_local_component_assembly);
280 }
281 template<class Result>
282 void valued_check() const {
283 typedef typename details::generic_unary_traits<UnaryFunction>::template hint<typename Expr::value_type,Result>::argument_type
284 A1;
285 if (! is_undeterminated<A1>::value) _expr.template valued_check<A1>();
286 }
287protected:
288// data:
289 UnaryFunction _f;
290 Expr _expr;
291};
292template<class F, class Expr> struct is_field_expr_v2_variational_arg <field_expr_v2_variational_unary<F,Expr> > : std::true_type {};
293
294} // namespace details
295
296// ---------------------------------------------------------------------------
297// 1.2. unary calls
298// ---------------------------------------------------------------------------
299
300#define _RHEOLEF_make_field_expr_v2_variational_unary_operator(FUNCTION,FUNCTOR) \
301template<class Expr> \
302inline \
303typename \
304std::enable_if< \
305 details::is_field_expr_v2_variational_arg<Expr>::value \
306 ,details::field_expr_v2_variational_unary< \
307 FUNCTOR \
308 ,Expr \
309 > \
310>::type \
311FUNCTION (const Expr& expr) \
312{ \
313 return details::field_expr_v2_variational_unary <FUNCTOR,Expr> (FUNCTOR(), expr); \
314}
315
320#undef _RHEOLEF_make_field_expr_v2_variational_unary_operator
321
322// ---------------------------------------------------------------------------
323// 2. binary operators +- between variationals
324// ---------------------------------------------------------------------------
325// 2.1. binary node
326// ---------------------------------------------------------------------------
327// ex: v+v, v-v
328
329namespace details {
330
331template<class BinaryFunction, class Expr1, class Expr2>
333public:
334// typedefs:
335
340 typename Expr1::value_type
341 ,typename Expr2::value_type>::type result_hint;
343 typename Expr1::value_type
344 ,typename Expr2::value_type
348 typedef space_basic<scalar_type,memory_type> space_type; // TODO: deduce from Exprs
349 typedef typename details::bf_vf_tag<BinaryFunction,
350 typename Expr1::vf_tag_type,
351 typename Expr2::vf_tag_type>::type vf_tag_type;
357
358// alocators:
359
361 const Expr1& expr1,
362 const Expr2& expr2)
363 : _f(f), _expr1(expr1), _expr2(expr2) {}
364
365// accessors:
366
367 static bool have_test_space() { return true; }
368 const space_type& get_vf_space() const { return _expr1.get_vf_space(); }
372 }
373 size_type n_derivative() const { return std::min(_expr1.n_derivative(), _expr2.n_derivative()); }
374
375// mutable modifiers:
376
377 // TODO: check that expr1 & expr2 have the same get_vf_space()
379 _expr1.initialize (pops, iopt);
380 _expr2.initialize (pops, iopt);
381 }
383 _expr1.initialize (gh, pops, iopt);
384 _expr2.initialize (gh, pops, iopt);
385 }
386 // evaluate when all arg types are determinated
387 template<class Result, class Arg1, class Arg2, class Status>
390 const self_type& obj,
392 const geo_element& K,
393 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value) const
394 {
395 fatal_macro ("invalid type resolution");
396 }
398 const self_type& obj,
400 const geo_element& K,
401 const side_information_type& sid,
402 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value) const
403 {
404 fatal_macro ("invalid type resolution");
405 }
406 };
407 template<class Result, class Arg1, class Arg2>
408 struct evaluate_call_check<Result,Arg1,Arg2,std::true_type> {
410 const self_type& obj,
412 const geo_element& K,
413 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value) const
414 {
415 Eigen::Matrix<Arg1,Eigen::Dynamic,Eigen::Dynamic> value1; obj._expr1.evaluate (omega_K, K, value1);
416 Eigen::Matrix<Arg2,Eigen::Dynamic,Eigen::Dynamic> value2; obj._expr2.evaluate (omega_K, K, value2);
417 // TODO: DVT_EIGEN_BLAS2
418 value.resize (value1.rows(), value1.cols());
419 for (size_type loc_inod = 0, loc_nnod = value1.rows(); loc_inod < loc_nnod; ++loc_inod) {
420 for (size_type loc_jdof = 0, loc_ndof = value1.cols(); loc_jdof < loc_ndof; ++loc_jdof) {
421 value(loc_inod,loc_jdof) = obj._f (value1(loc_inod,loc_jdof), value2(loc_inod,loc_jdof));
422 }}
423 }
425 const self_type& obj,
427 const geo_element& K,
428 const side_information_type& sid,
429 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value,
430 bool do_local_component_assembly) const
431 {
432 Eigen::Matrix<Arg1,Eigen::Dynamic,Eigen::Dynamic> value1; obj._expr1.evaluate_on_side (omega_K, K, sid, value1, do_local_component_assembly);
433 Eigen::Matrix<Arg2,Eigen::Dynamic,Eigen::Dynamic> value2; obj._expr2.evaluate_on_side (omega_K, K, sid, value2, do_local_component_assembly);
434 // TODO: DVT_EIGEN_BLAS2
435 value.resize (value1.rows(), value1.cols());
436 for (size_type loc_inod = 0, loc_nnod = value1.rows(); loc_inod < loc_nnod; ++loc_inod) {
437 for (size_type loc_jdof = 0, loc_ndof = value1.cols(); loc_jdof < loc_ndof; ++loc_jdof) {
438 value(loc_inod,loc_jdof) = obj._f (value1(loc_inod,loc_jdof), value2(loc_inod,loc_jdof));
439 }}
440 }
441 };
442 template<class Result, class Arg1, class Arg2>
445 const geo_element& K,
446 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value) const
447 {
449 // TODO: instead of is_equal, could have compatible scalars T1,T2 ?
450 typedef typename details::is_equal<Result,result_type>::type status_t;
452 eval (*this, omega_K, K, value);
453 }
454 template<class Result, class Arg1, class Arg2>
457 const geo_element& K,
458 const side_information_type& sid,
459 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value,
460 bool do_local_component_assembly) const
461 {
463 // TODO: instead of is_equal, could have compatible scalars T1,T2 ?
464 typedef typename details::is_equal<Result,result_type>::type status_t;
466 eval (*this, omega_K, K, sid, value, do_local_component_assembly);
467 }
468 template<class This, class Result, class Arg1, class Arg2, class Undet1, class Undet2>
470 // --------------------------------------------------------------------------
471 // when both args are defined at compile time:
472 // --------------------------------------------------------------------------
473 template<class This, class Result, class Arg1, class Arg2>
474 struct evaluate_switch<This, Result, Arg1, Arg2, std::false_type, std::false_type> {
476 const This& obj,
478 const geo_element& K,
479 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value) const
480 {
481 obj.template evaluate_internal<Result, Arg1, Arg2> (omega_K, K, value);
482 }
484 const This& obj,
486 const geo_element& K,
487 const side_information_type& sid,
488 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value,
489 bool do_local_component_assembly) const
490 {
491 obj.template evaluate_on_side_internal<Result, Arg1, Arg2> (omega_K, K, sid, value, do_local_component_assembly);
492 }
493 };
494 // --------------------------------------------------------------------------
495 // when first arg is undeterminated
496 // --------------------------------------------------------------------------
497 template<class This, class Result, class Arg1, class Arg2>
498 struct evaluate_switch<This, Result, Arg1, Arg2, std::true_type, std::false_type> {
500 const This& obj,
502 const geo_element& K,
503 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value) const
504 {
505 typedef typename scalar_traits<Arg1>::type T1;
506 space_constant::valued_type arg1_valued_tag = obj._expr1.valued_tag();
507 switch (arg1_valued_tag) {
508#define _RHEOLEF_switch1(VALUED1,VALUE1) \
509 case space_constant::VALUED1: \
510 obj.template evaluate_internal<Result, VALUE1, Arg2>(omega_K, K, value); break;
517#undef _RHEOLEF_switch1
518 default: error_macro ("unexpected first argument valued tag="<<arg1_valued_tag);
519 }
520 }
522 const This& obj,
524 const geo_element& K,
525 const side_information_type& sid,
526 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value,
527 bool do_local_component_assembly) const
528 {
529 typedef typename scalar_traits<Arg1>::type T1;
530 space_constant::valued_type arg1_valued_tag = obj._expr1.valued_tag();
531 switch (arg1_valued_tag) {
532#define _RHEOLEF_switch1(VALUED1,VALUE1) \
533 case space_constant::VALUED1: \
534 obj.template evaluate_on_side_internal<Result, VALUE1, Arg2>(omega_K, K, sid, value, do_local_component_assembly); break;
541#undef _RHEOLEF_switch1
542 default: error_macro ("unexpected first argument valued tag="<<arg1_valued_tag);
543 }
544 }
545 };
546 // --------------------------------------------------------------------------
547 // when second arg is undeterminated
548 // --------------------------------------------------------------------------
549 template<class This, class Result, class Arg1, class Arg2>
550 struct evaluate_switch<This, Result, Arg1, Arg2, std::false_type, std::true_type> {
552 const This& obj,
554 const geo_element& K,
555 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value) const
556 {
557 typedef typename scalar_traits<Arg2>::type T2;
558 space_constant::valued_type arg2_valued_tag = obj._expr2.valued_tag();
559 switch (arg2_valued_tag) {
560#define _RHEOLEF_switch2(VALUED2,VALUE2) \
561 case space_constant::VALUED2: \
562 obj.template evaluate_internal<Result, Arg1, VALUE2>(omega_K, K, value); break;
569#undef _RHEOLEF_switch2
570 default: error_macro ("unexpected second argument valued tag="<<arg2_valued_tag);
571 }
572 }
574 const This& obj,
576 const geo_element& K,
577 const side_information_type& sid,
578 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value,
579 bool do_local_component_assembly) const
580 {
581 typedef typename scalar_traits<Arg2>::type T2;
582 space_constant::valued_type arg2_valued_tag = obj._expr2.valued_tag();
583 switch (arg2_valued_tag) {
584#define _RHEOLEF_switch2(VALUED2,VALUE2) \
585 case space_constant::VALUED2: \
586 obj.template evaluate_on_side_internal<Result, Arg1, VALUE2>(omega_K, K, sid, value, do_local_component_assembly); break;
593#undef _RHEOLEF_switch2
594 default: error_macro ("unexpected second argument valued tag="<<arg2_valued_tag);
595 }
596 }
597 };
598 // --------------------------------------------------------------------------
599 // when both are undefined at compile time:
600 // --------------------------------------------------------------------------
601 template<class This, class Result, class Arg1, class Arg2>
602 struct evaluate_switch<This, Result, Arg1, Arg2, std::true_type, std::true_type> {
604 const This& obj,
606 const geo_element& K,
607 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value) const
608 {
609 typedef typename scalar_traits<Arg1>::type T1;
610 typedef typename scalar_traits<Arg2>::type T2;
611 space_constant::valued_type arg1_valued_tag = obj._expr1.valued_tag();
612 space_constant::valued_type arg2_valued_tag = obj._expr2.valued_tag();
613 switch (arg1_valued_tag) {
614#define _RHEOLEF_switch2(VALUE1,VALUED2,VALUE2) \
615 case space_constant::VALUED2: \
616 obj.template evaluate_internal<Result, VALUE1, VALUE2>(omega_K, K, value); break;
617#define _RHEOLEF_switch1(VALUED1,VALUE1) \
618 case space_constant::VALUED1: { \
619 switch (arg2_valued_tag) { \
620_RHEOLEF_switch2(VALUE1,scalar,T2) \
621_RHEOLEF_switch2(VALUE1,vector,point_basic<T2>) \
622_RHEOLEF_switch2(VALUE1,tensor,tensor_basic<T2>) \
623_RHEOLEF_switch2(VALUE1,unsymmetric_tensor,tensor_basic<T2>) \
624_RHEOLEF_switch2(VALUE1,tensor3,tensor3_basic<T2>) \
625_RHEOLEF_switch2(VALUE1,tensor4,tensor4_basic<T2>) \
626 default: error_macro ("unexpected second argument valued tag="<<arg2_valued_tag); \
627 } \
628 break; \
629 }
636#undef _RHEOLEF_switch2
637#undef _RHEOLEF_switch1
638 default: error_macro ("unexpected first argument valued tag="<<arg1_valued_tag);
639 }
640 }
642 const This& obj,
644 const geo_element& K,
645 const side_information_type& sid,
646 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value,
647 bool do_local_component_assembly) const
648 {
649 typedef typename scalar_traits<Arg1>::type T1;
650 typedef typename scalar_traits<Arg2>::type T2;
651 space_constant::valued_type arg1_valued_tag = obj._expr1.valued_tag();
652 space_constant::valued_type arg2_valued_tag = obj._expr2.valued_tag();
653 switch (arg1_valued_tag) {
654#define _RHEOLEF_switch2(VALUE1,VALUED2,VALUE2) \
655 case space_constant::VALUED2: \
656 obj.template evaluate_on_side_internal<Result, VALUE1, VALUE2>(omega_K, K, sid, value, do_local_component_assembly); break;
657#define _RHEOLEF_switch1(VALUED1,VALUE1) \
658 case space_constant::VALUED1: { \
659 switch (arg2_valued_tag) { \
660_RHEOLEF_switch2(VALUE1,scalar,T2) \
661_RHEOLEF_switch2(VALUE1,vector,point_basic<T2>) \
662_RHEOLEF_switch2(VALUE1,tensor,tensor_basic<T2>) \
663_RHEOLEF_switch2(VALUE1,unsymmetric_tensor,tensor_basic<T2>) \
664_RHEOLEF_switch2(VALUE1,tensor3,tensor3_basic<T2>) \
665_RHEOLEF_switch2(VALUE1,tensor4,tensor4_basic<T2>) \
666 default: error_macro ("unexpected second argument valued tag="<<arg2_valued_tag); \
667 } \
668 break; \
669 }
676#undef _RHEOLEF_switch2
677#undef _RHEOLEF_switch1
678 default: error_macro ("unexpected first argument valued tag="<<arg1_valued_tag);
679 }
680 }
681 };
682 // --------------------------------------------------------------------------
683 // main eval call:
684 // --------------------------------------------------------------------------
685 template<class Result>
686 struct hint {
688 typename Expr1::value_type
689 ,typename Expr2::value_type
690 ,Result>::first_argument_type A1;
692 typename Expr1::value_type
693 ,typename Expr2::value_type
694 ,Result>::second_argument_type A2;
695 };
696 template<class Result>
697 void evaluate (
699 const geo_element& K,
700 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value) const
701 {
702 typedef typename hint<Result>::A1 A1;
703 typedef typename hint<Result>::A2 A2;
704 typedef typename is_undeterminated<A1>::type undet_1;
705 typedef typename is_undeterminated<A2>::type undet_2;
707 eval (*this, omega_K, K, value);
708 }
709 template<class Result>
712 const geo_element& K,
713 const side_information_type& sid,
714 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value,
715 bool do_local_component_assembly) const
716 {
717 typedef typename hint<Result>::A1 A1;
718 typedef typename hint<Result>::A2 A2;
719 typedef typename is_undeterminated<A1>::type undet_1;
720 typedef typename is_undeterminated<A2>::type undet_2;
722 eval (*this, omega_K, K, sid, value, do_local_component_assembly);
723 }
724 template<class Result>
725 void valued_check() const {
726 typedef typename hint<Result>::A1 A1;
727 typedef typename hint<Result>::A2 A2;
728 if (! is_undeterminated<A1>::value) _expr1.template valued_check<A1>();
729 if (! is_undeterminated<A2>::value) _expr2.template valued_check<A2>();
730 }
731protected:
732// data:
735 Expr2 _expr2;
736};
737template <class F, class Expr1, class Expr2> struct is_field_expr_v2_variational_arg <field_expr_v2_variational_binary<F,Expr1,Expr2> > : std::true_type {};
738
739} // namespace details
740// ---------------------------------------------------------------------------
741// 2.2. binary call
742// ---------------------------------------------------------------------------
743namespace details {
744
745template<class Expr1, class Expr2, class Sfinae = void>
747
748template<class Expr1, class Expr2>
750 Expr1
751 ,Expr2
752 ,typename
753 std::enable_if<
754 is_field_expr_v2_variational_arg<Expr1>::value
755 && is_field_expr_v2_variational_arg<Expr2>::value
756 >::type
757>
758: and_type<
759 is_field_expr_v2_variational_arg<Expr1>
760 ,is_field_expr_v2_variational_arg<Expr2>
761 ,std::is_same <
762 typename Expr1::vf_tag_type
763 ,typename Expr2::vf_tag_type
764 >
765 >
766{};
767
768} // namespace details
769
770#define _RHEOLEF_make_field_expr_v2_variational_binary_operator_plus_minus(FUNCTION,FUNCTOR) \
771template<class Expr1, class Expr2> \
772inline \
773typename \
774std::enable_if< \
775 details::is_field_expr_v2_variational_binary_plus_minus <Expr1,Expr2>::value \
776 ,details::field_expr_v2_variational_binary< \
777 FUNCTOR \
778 ,Expr1 \
779 ,Expr2 \
780 > \
781>::type \
782FUNCTION (const Expr1& expr1, const Expr2& expr2) \
783{ \
784 return details::field_expr_v2_variational_binary <FUNCTOR, Expr1, Expr2> (FUNCTOR(), expr1, expr2); \
785}
786
789#undef _RHEOLEF_make_field_expr_v2_variational_binary_operator_plus_minus
790
791// ---------------------------------------------------------------------------
792// 3. binary operators */ between a variational and a nonlinear expr
793// ---------------------------------------------------------------------------
794// 3.1. helper
795// ---------------------------------------------------------------------------
796namespace details {
797
798 template<class This, class Arg1>
799 struct nl_switch {
800 typedef typename This::size_type size_type;
801 void element_initialize (const This& obj, const geo_element& K) const {
802 space_constant::valued_type nl_arg_valued_tag = obj._nl_expr.valued_tag();
803 switch (nl_arg_valued_tag) {
805 obj._nl_expr.evaluate (K, obj._scalar_nl_value_quad); break;
807 obj._nl_expr.evaluate (K, obj._vector_nl_value_quad); break;
810 obj._nl_expr.evaluate (K, obj._tensor_nl_value_quad); break;
812 obj._nl_expr.evaluate (K, obj._tensor3_nl_value_quad); break;
814 obj._nl_expr.evaluate (K, obj._tensor4_nl_value_quad); break;
815 default: error_macro ("unexpected first argument valued tag="<<nl_arg_valued_tag); // ICI
816 }
817 }
818 void element_initialize_on_side (const This& obj, const geo_element& K, const side_information_type& sid) const {
819 space_constant::valued_type nl_arg_valued_tag = obj._nl_expr.valued_tag();
820 switch (nl_arg_valued_tag) {
822 obj._nl_expr.evaluate_on_side (K, sid, obj._scalar_nl_value_quad); break;
824 obj._nl_expr.evaluate_on_side (K, sid, obj._vector_nl_value_quad); break;
827 obj._nl_expr.evaluate_on_side (K, sid, obj._tensor_nl_value_quad); break;
829 obj._nl_expr.evaluate_on_side (K, sid, obj._tensor3_nl_value_quad); break;
831 obj._nl_expr.evaluate_on_side (K, sid, obj._tensor4_nl_value_quad); break;
832 default: error_macro ("unexpected first argument valued tag="<<nl_arg_valued_tag);
833 }
834 }
835 Arg1 get_nl_value (const This& obj, size_type q) const {
836 // Arg1 may be solved at compile time for real evaluation
837 fatal_macro ("unexpected argument type="<<typename_macro(Arg1));
838 return Arg1();
839 }
840 };
841 template<class This> struct nl_switch<This,typename This::scalar_type> {
842 typedef typename This::size_type size_type;
843 typedef typename This::scalar_type scalar_type;
844 void element_initialize (const This& obj, const geo_element& K) const {
845 obj._nl_expr.evaluate (K, obj._scalar_nl_value_quad); }
846 void element_initialize_on_side (const This& obj, const geo_element& K, const side_information_type& sid) const {
847 obj._nl_expr.evaluate_on_side (K, sid, obj._scalar_nl_value_quad); }
848 const scalar_type& get_nl_value (const This& obj, size_type q) const {
849 return obj._scalar_nl_value_quad[q]; }
850 };
851 template<class This> struct nl_switch<This,point_basic<typename This::scalar_type> > {
852 typedef typename This::size_type size_type;
853 typedef typename This::scalar_type scalar_type;
854 void element_initialize (const This& obj, const geo_element& K) const {
855 obj._nl_expr.evaluate (K, obj._vector_nl_value_quad); }
856 void element_initialize_on_side (const This& obj, const geo_element& K, const side_information_type& sid) const {
857 obj._nl_expr.evaluate_on_side (K, sid, obj._vector_nl_value_quad); }
858 const point_basic<scalar_type>& get_nl_value (const This& obj, size_type q) const {
859 return obj._vector_nl_value_quad[q]; }
860 };
861 template<class This> struct nl_switch<This,tensor_basic<typename This::scalar_type> > {
862 typedef typename This::size_type size_type;
863 typedef typename This::scalar_type scalar_type;
864 void element_initialize (const This& obj, const geo_element& K) const {
865 obj._nl_expr.evaluate (K, obj._tensor_nl_value_quad); }
866 void element_initialize_on_side (const This& obj, const geo_element& K, const side_information_type& sid) const {
867 obj._nl_expr.evaluate_on_side (K, sid, obj._tensor_nl_value_quad); }
868 const tensor_basic<scalar_type>& get_nl_value (const This& obj, size_type q) const {
869 return obj._tensor_nl_value_quad[q]; }
870 };
871 template<class This> struct nl_switch<This,tensor3_basic<typename This::scalar_type> > {
872 typedef typename This::size_type size_type;
873 typedef typename This::scalar_type scalar_type;
874 void element_initialize (const This& obj, const geo_element& K) const {
875 obj._nl_expr.evaluate (K, obj._tensor3_nl_value_quad); }
876 void element_initialize_on_side (const This& obj, const geo_element& K, const side_information_type& sid) const {
877 obj._nl_expr.evaluate_on_side (K, sid, obj._tensor3_nl_value_quad); }
878 const tensor3_basic<scalar_type>& get_nl_value (const This& obj, size_type q) const {
879 return obj._tensor3_nl_value_quad[q]; }
880 };
881 template<class This> struct nl_switch<This,tensor4_basic<typename This::scalar_type> > {
882 typedef typename This::size_type size_type;
883 typedef typename This::scalar_type scalar_type;
884 void element_initialize (const This& obj, const geo_element& K) const {
885 obj._nl_expr.evaluate (K, obj._tensor4_nl_value_quad); }
886 void element_initialize_on_side (const This& obj, const geo_element& K, const side_information_type& sid) const {
887 obj._nl_expr.evaluate_on_side (K, sid, obj._tensor4_nl_value_quad); }
888 const tensor4_basic<scalar_type>& get_nl_value (const This& obj, size_type q) const {
889 return obj._tensor4_nl_value_quad[q]; }
890 };
891
892} // namespace details
893
894// ---------------------------------------------------------------------------
895// 3.2. binary node
896// ---------------------------------------------------------------------------
897//
898// function call: (f nl_expr vf_expr)
899// examples: f = operator*, operator/
900// eta_h*v
901// v/eta_h
902// dot(v,normal())
903// at any quadrature node xq, the compuation eta_h(xq) is performed
904// and then we loop on the basis functions for v :
905// eta_q = eta_h(xq);
906// for i=0..nk-1
907// value[i] = f (eta_q, v(xq)[i]);
908// since we can swap the two args (see the details::swap_fun<f> class),
909// we assume that the first argument is a field or a general field_nl_expr
910// and that the second argument is a test of a general field_vf_expr
911//
912// Implementation note: this operation do not reduces to field_expr_v2_variational_unary
913// with a class-function that contains eta_h since :
914// - the value of eta_h may be refreshed at each xq
915// (this could be achieved by replacing std::binder1st with an adequate extension)
916// - the valued category of eta_h is not always known at compile-time.
917// It is known in dot(eta_h,v) but not with eta_h*v
918// and the class-functions for field_expr_v2_variational_unary may have Arg1 and Result determined.
919// So we switch to a specific field_expr_v2_variational_binary_binded that is abble to solve the
920// valued type at run time. When it is possible, it is determined at compile-time.
921//
922namespace details {
923
924template<class BinaryFunction, class NLExpr, class VFExpr>
926public:
927// typedefs:
928
933 typename NLExpr::value_type
934 ,typename VFExpr::value_type>::type result_hint;
936 typename NLExpr::value_type
937 ,typename VFExpr::value_type
941 typedef space_basic<scalar_type,memory_type> space_type; // TODO: deduce from Exprs
942 typedef typename VFExpr::vf_tag_type vf_tag_type;
948
949// alocators:
950
952 const NLExpr& nl_expr,
953 const VFExpr& vf_expr)
954 : _f(f),
955 _nl_expr(nl_expr),
956 _vf_expr(vf_expr)
957 {}
958
960 : _f(x._f),
963 {}
964
965// accessors:
966
967 static bool have_test_space() { return true; } // deduce & check !
968 const space_type& get_vf_space() const { return _vf_expr.get_vf_space(); }
972 }
973 size_type n_derivative() const { return _vf_expr.n_derivative(); }
974
975// mutable modifiers:
976
978 _nl_expr.initialize (pops, iopt);
979 _vf_expr.initialize (pops, iopt);
980 }
982 _nl_expr.initialize ( pops, iopt);
983 _vf_expr.initialize (gh, pops, iopt);
984 }
985 // ---------------------------------------------
986 // evaluate
987 // ---------------------------------------------
988 // evaluate when all arg types are determinated
989 template<class Result, class Arg1, class Arg2, class Status>
992 const self_type& obj,
993 const geo_basic<float_type,memory_type>& omega_K,
994 const geo_element& K,
995 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value) const
996 {
997 fatal_macro ("invalid type resolution");
998 }
1000 const self_type& obj,
1001 const geo_basic<float_type,memory_type>& omega_K,
1002 const geo_element& K,
1003 const side_information_type& sid,
1004 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value,
1005 bool do_local_component_assembly) const
1006 {
1007 fatal_macro ("invalid type resolution");
1008 }
1009 };
1010 template<class Result, class Arg1, class Arg2>
1011 struct evaluate_call_check<Result,Arg1,Arg2,std::true_type> {
1013 const self_type& obj,
1014 const geo_basic<float_type,memory_type>& omega_K,
1015 const geo_element& K,
1016 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value) const
1017 {
1018 Eigen::Matrix<Arg1,Eigen::Dynamic,1> value1; obj._nl_expr.evaluate (omega_K, K, value1);
1019 Eigen::Matrix<Arg2,Eigen::Dynamic,Eigen::Dynamic> value2; obj._vf_expr.evaluate (omega_K, K, value2);
1020 check_macro (value1.size() == value2.rows(), "invalid sizes value1(nnod="<<value1.size()
1021 <<") and value2(nnod="<<value2.rows()<<",ndof="<<value2.cols()<<")");
1022 value.resize (value2.rows(), value2.cols());
1023 for (size_type loc_inod = 0, loc_nnod = value2.rows(); loc_inod < loc_nnod; ++loc_inod) {
1024 for (size_type loc_jdof = 0, loc_ndof = value2.cols(); loc_jdof < loc_ndof; ++loc_jdof) {
1025 value(loc_inod,loc_jdof) = obj._f (value1[loc_inod], value2(loc_inod,loc_jdof)); // TODO: DVT_EIGEN_BLAS2
1026 }}
1027 }
1029 const self_type& obj,
1030 const geo_basic<float_type,memory_type>& omega_K,
1031 const geo_element& K,
1032 const side_information_type& sid,
1033 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value,
1034 bool do_local_component_assembly) const
1035 {
1036 Eigen::Matrix<Arg1,Eigen::Dynamic,1> value1; obj._nl_expr.evaluate_on_side (omega_K, K, sid, value1);
1037 Eigen::Matrix<Arg2,Eigen::Dynamic,Eigen::Dynamic> value2; obj._vf_expr.evaluate_on_side (omega_K, K, sid, value2, do_local_component_assembly);
1038 check_macro (value1.size() == value2.rows(), "invalid sizes value1(nnod="<<value1.size()
1039 <<") and value2(nnod="<<value2.rows()<<",ndof="<<value2.cols()<<")");
1040 value.resize (value2.rows(), value2.cols());
1041 for (size_type loc_inod = 0, loc_nnod = value2.rows(); loc_inod < loc_nnod; ++loc_inod) {
1042 for (size_type loc_jdof = 0, loc_ndof = value2.cols(); loc_jdof < loc_ndof; ++loc_jdof) {
1043 value(loc_inod,loc_jdof) = obj._f (value1[loc_inod], value2(loc_inod,loc_jdof)); // TODO: DVT_EIGEN_BLAS2
1044 }}
1045 }
1046 };
1047 template<class Result, class Arg1, class Arg2>
1049 const geo_basic<float_type,memory_type>& omega_K,
1050 const geo_element& K,
1051 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value) const
1052 {
1054 // TODO: instead of is_equal, could have compatible scalars T1,T2 ?
1055 typedef typename details::is_equal<Result,result_type>::type status_t;
1057 eval (*this, omega_K, K, value);
1058 }
1059 template<class Result, class Arg1, class Arg2>
1061 const geo_basic<float_type,memory_type>& omega_K,
1062 const geo_element& K,
1063 const side_information_type& sid,
1064 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value,
1065 bool do_local_component_assembly) const
1066 {
1068 // TODO: instead of is_equal, could have compatible scalars T1,T2 ?
1069 typedef typename details::is_equal<Result,result_type>::type status_t;
1071 eval (*this, omega_K, K, sid, value, do_local_component_assembly);
1072 }
1073 template<class This, class Result, class Arg1, class Arg2, class Undet1, class Undet2>
1075 // --------------------------------------------------------------------------
1076 // when both args are defined at compile time:
1077 // --------------------------------------------------------------------------
1078 template<class This, class Result, class Arg1, class Arg2>
1079 struct evaluate_switch<This, Result, Arg1, Arg2, std::false_type, std::false_type> {
1081 const This& obj,
1082 const geo_basic<float_type,memory_type>& omega_K,
1083 const geo_element& K,
1084 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value) const
1085 {
1086 obj.template evaluate_internal<Result, Arg1, Arg2> (omega_K, K, value);
1087 }
1089 const This& obj,
1090 const geo_basic<float_type,memory_type>& omega_K,
1091 const geo_element& K,
1092 const side_information_type& sid,
1093 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value,
1094 bool do_local_component_assembly) const
1095 {
1096 obj.template evaluate_on_side_internal<Result, Arg1, Arg2> (omega_K, K, sid, value, do_local_component_assembly);
1097 }
1098 };
1099 // --------------------------------------------------------------------------
1100 // when first arg is undeterminated
1101 // --------------------------------------------------------------------------
1102 template<class This, class Result, class Arg1, class Arg2>
1103 struct evaluate_switch<This, Result, Arg1, Arg2, std::true_type, std::false_type> {
1105 const This& obj,
1106 const geo_basic<float_type,memory_type>& omega_K,
1107 const geo_element& K,
1108 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value) const
1109 {
1110 typedef typename scalar_traits<Arg1>::type T1;
1111 space_constant::valued_type arg1_valued_tag = obj._nl_expr.valued_tag();
1112 switch (arg1_valued_tag) {
1113#define _RHEOLEF_switch(A1_VALUED,A1_VALUE) \
1114 case space_constant::A1_VALUED: \
1115 obj.template evaluate_internal<Result, A1_VALUE, Arg2> (omega_K, K, value); break;
1122#undef _RHEOLEF_switch
1123 default: error_macro ("unexpected first argument valued tag="<<arg1_valued_tag);
1124 }
1125 }
1127 const This& obj,
1128 const geo_basic<float_type,memory_type>& omega_K,
1129 const geo_element& K,
1130 const side_information_type& sid,
1131 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value,
1132 bool do_local_component_assembly) const
1133 {
1134 typedef typename scalar_traits<Arg1>::type T1;
1135 space_constant::valued_type arg1_valued_tag = obj._nl_expr.valued_tag();
1136 switch (arg1_valued_tag) {
1137#define _RHEOLEF_switch(A1_VALUED,A1_VALUE) \
1138 case space_constant::A1_VALUED: \
1139 obj.template evaluate_on_side_internal<Result, A1_VALUE, Arg2> (omega_K, K, sid, value, do_local_component_assembly); break;
1146#undef _RHEOLEF_switch
1147 default: error_macro ("unexpected first argument valued tag="<<arg1_valued_tag);
1148 }
1149 }
1150 };
1151 // --------------------------------------------------------------------------
1152 // when second arg is undeterminated
1153 // --------------------------------------------------------------------------
1154 template<class This, class Result, class Arg1, class Arg2>
1155 struct evaluate_switch<This, Result, Arg1, Arg2, std::false_type, std::true_type> {
1157 const This& obj,
1158 const geo_basic<float_type,memory_type>& omega_K,
1159 const geo_element& K,
1160 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value) const
1161 {
1162 typedef typename scalar_traits<Arg2>::type T2;
1163 space_constant::valued_type arg2_valued_tag = obj._vf_expr.valued_tag();
1164 switch (arg2_valued_tag) {
1165#define _RHEOLEF_switch(A2_VALUED,A2_VALUE) \
1166 case space_constant::A2_VALUED: \
1167 obj.template evaluate_internal<Result, Arg1, A2_VALUE> (omega_K, K, value); break;
1174#undef _RHEOLEF_switch
1175 default: error_macro ("unexpected second argument valued tag="<<arg2_valued_tag);
1176 }
1177 }
1179 const This& obj,
1180 const geo_basic<float_type,memory_type>& omega_K,
1181 const geo_element& K,
1182 const side_information_type& sid,
1183 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value,
1184 bool do_local_component_assembly) const
1185 {
1186 typedef typename scalar_traits<Arg2>::type T2;
1187 space_constant::valued_type arg2_valued_tag = obj._vf_expr.valued_tag();
1188 switch (arg2_valued_tag) {
1189#define _RHEOLEF_switch(A2_VALUED,A2_VALUE) \
1190 case space_constant::A2_VALUED: \
1191 obj.template evaluate_on_side_internal<Result, Arg1, A2_VALUE> (omega_K, K, sid, value, do_local_component_assembly); break;
1198#undef _RHEOLEF_switch
1199 default: error_macro ("unexpected second argument valued tag="<<arg2_valued_tag);
1200 }
1201 }
1202 };
1203 // --------------------------------------------------------------------------
1204 // when both are undefined at compile time:
1205 // --------------------------------------------------------------------------
1206 template<class This, class Result, class Arg1, class Arg2>
1207 struct evaluate_switch<This, Result, Arg1, Arg2, std::true_type, std::true_type> {
1209 const This& obj,
1210 const geo_basic<float_type,memory_type>& omega_K,
1211 const geo_element& K,
1212 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value) const
1213 {
1214 typedef typename scalar_traits<Arg1>::type T1;
1215 typedef typename scalar_traits<Arg2>::type T2;
1216 space_constant::valued_type arg1_valued_tag = obj._nl_expr.valued_tag();
1217 space_constant::valued_type arg2_valued_tag = obj._vf_expr.valued_tag();
1218 switch (arg1_valued_tag) {
1219#define _RHEOLEF_switch_A2(A1_VALUE,A2_VALUED,A2_VALUE) \
1220 case space_constant::A2_VALUED: \
1221 obj.template evaluate_internal<Result, A1_VALUE, A2_VALUE> (omega_K, K, value); break;
1222
1223#define _RHEOLEF_switch(A1_VALUED,A1_VALUE) \
1224 case space_constant::A1_VALUED: { \
1225 switch (arg2_valued_tag) { \
1226_RHEOLEF_switch_A2(A1_VALUE,scalar,T2) \
1227_RHEOLEF_switch_A2(A1_VALUE,vector,point_basic<T2>) \
1228_RHEOLEF_switch_A2(A1_VALUE,tensor,tensor_basic<T2>) \
1229_RHEOLEF_switch_A2(A1_VALUE,unsymmetric_tensor,tensor_basic<T2>) \
1230_RHEOLEF_switch_A2(A1_VALUE,tensor3,tensor3_basic<T2>) \
1231_RHEOLEF_switch_A2(A1_VALUE,tensor4,tensor4_basic<T2>) \
1232 default: error_macro ("unexpected second argument valued tag="<<arg2_valued_tag); \
1233 } \
1234 break; \
1235 }
1242#undef _RHEOLEF_switch
1243#undef _RHEOLEF_switch_A2
1244 default: error_macro ("unexpected first argument valued tag="<<arg1_valued_tag);
1245 }
1246 }
1248 const This& obj,
1249 const geo_basic<float_type,memory_type>& omega_K,
1250 const geo_element& K,
1251 const side_information_type& sid,
1252 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value,
1253 bool do_local_component_assembly) const
1254 {
1255 typedef typename scalar_traits<Arg1>::type T1;
1256 typedef typename scalar_traits<Arg2>::type T2;
1257 space_constant::valued_type arg1_valued_tag = obj._nl_expr.valued_tag();
1258 space_constant::valued_type arg2_valued_tag = obj._vf_expr.valued_tag();
1259 switch (arg1_valued_tag) {
1260#define _RHEOLEF_switch_A2(A1_VALUE,A2_VALUED,A2_VALUE) \
1261 case space_constant::A2_VALUED: \
1262 obj.template evaluate_on_side_internal<Result, A1_VALUE, A2_VALUE> (omega_K, K, sid, value, do_local_component_assembly); break;
1263
1264#define _RHEOLEF_switch(A1_VALUED,A1_VALUE) \
1265 case space_constant::A1_VALUED: { \
1266 switch (arg2_valued_tag) { \
1267_RHEOLEF_switch_A2(A1_VALUE,scalar,T2) \
1268_RHEOLEF_switch_A2(A1_VALUE,vector,point_basic<T2>) \
1269_RHEOLEF_switch_A2(A1_VALUE,tensor,tensor_basic<T2>) \
1270_RHEOLEF_switch_A2(A1_VALUE,unsymmetric_tensor,tensor_basic<T2>) \
1271_RHEOLEF_switch_A2(A1_VALUE,tensor3,tensor3_basic<T2>) \
1272_RHEOLEF_switch_A2(A1_VALUE,tensor4,tensor4_basic<T2>) \
1273 default: error_macro ("unexpected second argument valued tag="<<arg2_valued_tag); \
1274 } \
1275 break; \
1276 }
1283#undef _RHEOLEF_switch
1284#undef _RHEOLEF_switch_A2
1285 default: error_macro ("unexpected first argument valued tag="<<arg1_valued_tag);
1286 }
1287 }
1288 };
1289 // --------------------------------------------------------------------------
1290 // hint fort argument's types when Result is known
1291 // --------------------------------------------------------------------------
1292 template<class Result>
1293 struct hint {
1294 typedef typename promote<
1295 typename NLExpr::value_type
1297 typename NLExpr::value_type
1298 ,typename VFExpr::value_type
1299 ,Result>::first_argument_type
1301 typedef typename promote<
1302 typename VFExpr::value_type
1304 typename NLExpr::value_type
1305 ,typename VFExpr::value_type
1306 ,Result>::second_argument_type
1308 };
1309 // --------------------------------------------------------------------------
1310 // evaluate : main call with possibly runtime type resolutiion
1311 // --------------------------------------------------------------------------
1312 template<class Result>
1314 const geo_basic<float_type,memory_type>& omega_K,
1315 const geo_element& K,
1316 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value) const
1317 {
1318 typedef typename hint<Result>::A1 A1;
1319 typedef typename hint<Result>::A2 A2;
1322 typedef typename is_undeterminated<A1>::type undet_1;
1323 typedef typename is_undeterminated<A2>::type undet_2;
1325 eval (*this, omega_K, K, value);
1326 }
1327 template<class Result>
1329 const geo_basic<float_type,memory_type>& omega_K,
1330 const geo_element& K,
1331 const side_information_type& sid,
1332 Eigen::Matrix<Result,Eigen::Dynamic,Eigen::Dynamic>& value,
1333 bool do_local_component_assembly) const
1334 {
1335 typedef typename hint<Result>::A1 A1;
1336 typedef typename hint<Result>::A2 A2;
1339 typedef typename is_undeterminated<A1>::type undet_1;
1340 typedef typename is_undeterminated<A2>::type undet_2;
1342 eval (*this, omega_K, K, sid, value, do_local_component_assembly);
1343 }
1344 template<class Value>
1346 const geo_basic<float_type,memory_type>& omega_K,
1347 const geo_element& S,
1348 const geo_element& K0,
1349 const geo_element& K1,
1350 const Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value0,
1351 const Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value1,
1352 Eigen::Matrix<Value,Eigen::Dynamic,Eigen::Dynamic>& value) const
1353 {
1354 _vf_expr.local_dg_merge_on_side (omega_K, S, K0, K1, value0, value1, value);
1355 }
1356 template<class Result>
1357 void valued_check() const {
1358 typedef typename hint<Result>::A1 A1;
1359 typedef typename hint<Result>::A2 A2;
1360 if (! is_undeterminated<A1>::value) _nl_expr.template valued_check<A1>();
1361 if (! is_undeterminated<A2>::value) _vf_expr.template valued_check<A2>();
1362 }
1363//protected:
1364// data:
1366 NLExpr _nl_expr;
1367 VFExpr _vf_expr;
1368};
1369template<class F, class Expr1, class Expr2> struct is_field_expr_v2_variational_arg <field_expr_v2_variational_binary_binded<F,Expr1,Expr2> > : std::true_type {};
1370
1371} // namespace details
1372// ---------------------------------------------------------------------------
1373// 3.3. binary calls
1374// ---------------------------------------------------------------------------
1375namespace details {
1376
1377template<class Expr1, class Expr2, class Sfinae = void>
1379
1380template<class Expr1, class Expr2>
1382 Expr1
1383 ,Expr2
1384 ,typename
1385 std::enable_if<
1386 is_field_expr_v2_nonlinear_arg <Expr1>::value
1387 && ! is_field_expr_v2_constant <Expr1>::value
1388 && is_field_expr_v2_variational_arg<Expr2>::value
1389 >::type
1390>
1391: std::true_type
1392{};
1393
1394template<class Expr1, class Expr2>
1397
1398} // namespace details
1399
1400#define _RHEOLEF_make_field_expr_v2_variational_binary_operator_multiplies_divides_left(FUNCTION,FUNCTOR) \
1401template<class Expr1, class Expr2> \
1402inline \
1403typename \
1404std::enable_if< \
1405 details::is_field_expr_v2_variational_binary_multiplies_divides_left <Expr1,Expr2>::value \
1406 ,details::field_expr_v2_variational_binary_binded< \
1407 FUNCTOR \
1408 ,typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr1>::type \
1409 ,Expr2 /* vf */ \
1410 > \
1411>::type \
1412FUNCTION (const Expr1& expr1, const Expr2& expr2) \
1413{ \
1414 typedef typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr1>::type wrap1_t; \
1415 return details::field_expr_v2_variational_binary_binded \
1416 <FUNCTOR, wrap1_t, Expr2> \
1417 (FUNCTOR(), wrap1_t(expr1), expr2); \
1418}
1419
1420#define _RHEOLEF_make_field_expr_v2_variational_binary_operator_multiplies_divides_right(FUNCTION,FUNCTOR) \
1421template<class Expr1, class Expr2> \
1422inline \
1423typename \
1424std::enable_if< \
1425 details::is_field_expr_v2_variational_binary_multiplies_divides_right <Expr1,Expr2>::value \
1426 ,details::field_expr_v2_variational_binary_binded< \
1427 details::swapper<FUNCTOR> \
1428 , typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr2>::type \
1429 ,Expr1 /* vf */ \
1430 > \
1431>::type \
1432FUNCTION (const Expr1& expr1, const Expr2& expr2) \
1433{ \
1434 typedef typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr2>::type wrap2_t; \
1435 return details::field_expr_v2_variational_binary_binded \
1436 <details::swapper<FUNCTOR>, wrap2_t, Expr1> \
1437 (details::swapper<FUNCTOR>(FUNCTOR()), wrap2_t(expr2), expr1); \
1438}
1439#define _RHEOLEF_make_field_expr_v2_variational_binary_operator_multiplies_divides(FUNCTION,FUNCTOR) \
1440 _RHEOLEF_make_field_expr_v2_variational_binary_operator_multiplies_divides_left (FUNCTION,FUNCTOR) \
1441 _RHEOLEF_make_field_expr_v2_variational_binary_operator_multiplies_divides_right (FUNCTION,FUNCTOR)
1442
1448#undef _RHEOLEF_make_field_expr_v2_variational_binary_operator_multiplies_divides_left
1449#undef _RHEOLEF_make_field_expr_v2_variational_binary_operator_multiplies_divides_right
1450#undef _RHEOLEF_make_field_expr_v2_variational_binary_operator_multiplies_divides
1451
1452// ---------------------------------------------------------------------------
1453// 5. binary operators */ between one variational and a constant
1454// ---------------------------------------------------------------------------
1455namespace details {
1456
1457template<class Expr1, class Expr2, class Sfinae = void>
1458struct is_field_expr_v2_variational_binary_multiplies_divides_constant_left : std::false_type {};
1459
1460template<class Expr1, class Expr2>
1461struct is_field_expr_v2_variational_binary_multiplies_divides_constant_left <
1462 Expr1
1463 ,Expr2
1464 ,typename
1465 std::enable_if<
1466 is_field_expr_v2_constant <Expr1>::value
1467 && is_field_expr_v2_variational_arg<Expr2>::value
1468 >::type
1469>
1470: std::true_type
1471{};
1472
1473template<class Expr1, class Expr2>
1474struct is_field_expr_v2_variational_binary_multiplies_divides_constant_right
1475: is_field_expr_v2_variational_binary_multiplies_divides_constant_left <Expr2,Expr1> {};
1476
1477} // namespace details
1478
1479#define _RHEOLEF_make_field_expr_v2_variational_binary_operator_multiplies_divides_constant_left(FUNCTION,FUNCTOR) \
1480template<class Expr1, class Expr2> \
1481inline \
1482typename \
1483std::enable_if< \
1484 details::is_field_expr_v2_variational_binary_multiplies_divides_constant_left <Expr1,Expr2>::value \
1485 ,details::field_expr_v2_variational_unary< \
1486 details::binder_first <FUNCTOR, Expr1> \
1487 ,Expr2 /* vf */ \
1488 > \
1489>::type \
1490FUNCTION (const Expr1& expr1, const Expr2& expr2) \
1491{ \
1492 return details::field_expr_v2_variational_unary \
1493 <details::binder_first <FUNCTOR,Expr1>, Expr2> \
1494 (details::binder_first <FUNCTOR,Expr1> (FUNCTOR(), expr1), expr2); \
1495}
1496
1497#define _RHEOLEF_make_field_expr_v2_variational_binary_operator_multiplies_divides_constant_right(FUNCTION,FUNCTOR) \
1498template<class Expr1, class Expr2> \
1499inline \
1500typename \
1501std::enable_if< \
1502 details::is_field_expr_v2_variational_binary_multiplies_divides_constant_right <Expr1,Expr2>::value \
1503 ,details::field_expr_v2_variational_unary< \
1504 details::binder_second <FUNCTOR, Expr2> \
1505 ,Expr1 /* vf */ \
1506 > \
1507>::type \
1508FUNCTION (const Expr1& expr1, const Expr2& expr2) \
1509{ \
1510 return details::field_expr_v2_variational_unary \
1511 <details::binder_second <FUNCTOR,Expr2>, Expr1> \
1512 (details::binder_second <FUNCTOR,Expr2> (FUNCTOR(), expr2), expr1); \
1513}
1514
1515#define _RHEOLEF_make_field_expr_v2_variational_binary_operator_multiplies_divides_constant(FUNCTION,FUNCTOR) \
1516 _RHEOLEF_make_field_expr_v2_variational_binary_operator_multiplies_divides_constant_left (FUNCTION,FUNCTOR) \
1517 _RHEOLEF_make_field_expr_v2_variational_binary_operator_multiplies_divides_constant_right (FUNCTION,FUNCTOR)
1518
1519
1525
1526#undef _RHEOLEF_make_field_expr_v2_variational_binary_operator_multiplies_divides_constant_right
1527#undef _RHEOLEF_make_field_expr_v2_variational_binary_operator_multiplies_divides_constant_left
1528#undef _RHEOLEF_make_field_expr_v2_variational_binary_operator_multiplies_divides_constant
1529
1530} // namespace rheolef
1531#endif // _RHEOLEF_FIELD_EXPR_VARIATIONAL_H
field::size_type size_type
Definition: branch.cc:430
field gh(Float epsilon, Float t, const field &uh, const test &v)
void evaluate_on_side_internal(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
promote_memory< typenameNLExpr::memory_type, typenameVFExpr::memory_type >::type memory_type
field_expr_v2_variational_binary_binded(const field_expr_v2_variational_binary_binded< BinaryFunction, NLExpr, VFExpr > &x)
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)
field_expr_v2_variational_binary_binded< BinaryFunction, NLExpr, VFExpr > self_type
field_expr_v2_variational_binary_binded(const BinaryFunction &f, const NLExpr &nl_expr, const VFExpr &vf_expr)
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
details::generic_binary_traits< BinaryFunction >::template hint< typenameNLExpr::value_type, typenameVFExpr::value_type, result_hint >::result_type value_type
details::generic_binary_traits< BinaryFunction >::template result_hint< typenameNLExpr::value_type, typenameVFExpr::value_type >::type result_hint
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_binary_binded< BinaryFunction, NLExpr, typename VFExpr::dual_self_type > dual_self_type
void evaluate_internal(const geo_basic< float_type, memory_type > &omega_K, const geo_element &K, Eigen::Matrix< Result, Eigen::Dynamic, Eigen::Dynamic > &value) const
void evaluate_on_side_internal(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
space_basic< scalar_type, memory_type > space_type
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)
details::generic_binary_traits< BinaryFunction >::template result_hint< typenameExpr1::value_type, typenameExpr2::value_type >::type result_hint
field_expr_v2_variational_binary< BinaryFunction, typename Expr1::dual_self_type, typename Expr2::dual_self_type > dual_self_type
field_expr_v2_variational_binary< BinaryFunction, Expr1, Expr2 > self_type
field_expr_v2_variational_binary(const BinaryFunction &f, const Expr1 &expr1, const Expr2 &expr2)
details::bf_vf_tag< BinaryFunction, typenameExpr1::vf_tag_type, typenameExpr2::vf_tag_type >::type vf_tag_type
static const space_constant::valued_type valued_hint
details::dual_vf_tag< vf_tag_type >::type vf_dual_tag_type
details::generic_binary_traits< BinaryFunction >::template hint< typenameExpr1::value_type, typenameExpr2::value_type, result_hint >::result_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
promote_memory< typenameExpr1::memory_type, typenameExpr2::memory_type >::type memory_type
void evaluate_internal(const geo_basic< float_type, memory_type > &omega_K, const geo_element &K, Eigen::Matrix< Result, Eigen::Dynamic, Eigen::Dynamic > &value) const
space_basic< scalar_type, memory_type > space_type
field_expr_v2_variational_unary< UnaryFunction, typename Expr::dual_self_type > dual_self_type
field_expr_v2_variational_unary< UnaryFunction, Expr > self_type
details::generic_unary_traits< UnaryFunction >::template result_hint< typenameExpr::value_type >::type value_type
void initialize(const piola_on_pointset< float_type > &pops, const integrate_option &iopt)
field_expr_v2_variational_unary(const field_expr_v2_variational_unary< UnaryFunction, Expr > &x)
void initialize(const band_basic< float_type, memory_type > &gh, const piola_on_pointset< float_type > &pops, const integrate_option &iopt)
void evaluate_internal(const geo_basic< float_type, M > &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_internal(const geo_basic< float_type, M > &omega_K, const geo_element &K, Eigen::Matrix< Result, Eigen::Dynamic, Eigen::Dynamic > &value) const
static const space_constant::valued_type valued_hint
details::dual_vf_tag< vf_tag_type >::type vf_dual_tag_type
field_expr_v2_variational_unary(const UnaryFunction &f, const Expr &expr)
void evaluate(const geo_basic< float_type, M > &omega_K, const geo_element &K, Eigen::Matrix< Result, Eigen::Dynamic, Eigen::Dynamic > &value) const
void evaluate_on_side(const geo_basic< float_type, M > &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
see the geo_element page for the full documentation
Definition: geo_element.h:102
reference_element::size_type size_type
Definition: geo_element.h:125
see the integrate_option page for the full documentation
the finite element space
Definition: space.h:382
rheolef::std type
typename details::generic_binary_traits< BinaryFunction >::template result_hint< typename Expr1::result_type, typename Expr2::result_type >::type result_type
rheolef::std BinaryFunction
rheolef::std value
rheolef::std Expr1
dot(x,y): see the expression page for the full documentation
see the tensor3 page for the full documentation
see the tensor4 page for the full documentation
see the tensor page for the full documentation
#define error_macro(message)
Definition: dis_macros.h:49
#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_binary_operator_multiplies_divides_right(FUNCTION, FUNCTOR)
#define _RHEOLEF_switch(VALUED, VALUE)
#define _RHEOLEF_make_field_expr_v2_variational_binary_operator_multiplies_divides_constant_right(FUNCTION, FUNCTOR)
#define _RHEOLEF_switch2(VALUED2, VALUE2)
#define _RHEOLEF_switch1(VALUED1, VALUE1)
This file is part of Rheolef.
T ddot(const tensor_basic< T > &a, const tensor_basic< T > &b)
ddot(x,y): see the expression page for the full documentation
Definition: tensor.cc:278
U tr(const tensor_basic< U > &a, size_t d=3)
T dddot(const tensor3_basic< T > &a, const tensor3_basic< T > &b)
Definition: tensor3.cc:94
_RHEOLEF_make_field_expr_v2_variational_binary_operator_multiplies_divides(operator*, details::multiplies) _RHEOLEF_make_field_expr_v2_variational_binary_operator_multiplies_divides_right(operator/
rheolef::std enable_if ::type dot const Expr1 expr1, const Expr2 expr2 dot(const Expr1 &expr1, const Expr2 &expr2)
Definition: vec_expr_v2.h:415
_RHEOLEF_make_field_expr_v2_variational_binary_operator_plus_minus(operator+, details::plus) _RHEOLEF_make_field_expr_v2_variational_binary_operator_plus_minus(operator-
csr< T, sequential > trans(const csr< T, sequential > &a)
trans(a): see the form page for the full documentation
Definition: csr.h:455
_RHEOLEF_make_field_expr_v2_variational_binary_operator_multiplies_divides_constant(operator*, details::multiplies) _RHEOLEF_make_field_expr_v2_variational_binary_operator_multiplies_divides_constant_right(operator/
_RHEOLEF_make_field_expr_v2_variational_unary_operator(operator+, details::unary_plus) _RHEOLEF_make_field_expr_v2_variational_unary_operator(operator-
t operator()(const t &a, const t &b)
Definition: space.cc:386
STL namespace.
Definition: cavity_dg.h:29
void operator()(const self_type &obj, const geo_basic< float_type, memory_type > &omega_K, const geo_element &K, Eigen::Matrix< Result, Eigen::Dynamic, Eigen::Dynamic > &value) const
details::generic_binary_traits< BinaryFunction >::template hint< typenameExpr1::value_type, typenameExpr2::value_type, Result >::second_argument_type A2
details::generic_binary_traits< BinaryFunction >::template hint< typenameExpr1::value_type, typenameExpr2::value_type, Result >::first_argument_type A1
void operator()(const self_type &obj, const geo_basic< float_type, memory_type > &omega_K, const geo_element &K, Eigen::Matrix< Result, Eigen::Dynamic, Eigen::Dynamic > &value) const
promote< typenameNLExpr::value_type, typenamedetails::generic_binary_traits< BinaryFunction >::templatehint< typenameNLExpr::value_type, typenameVFExpr::value_type, Result >::first_argument_type >::type A1
promote< typenameVFExpr::value_type, typenamedetails::generic_binary_traits< BinaryFunction >::templatehint< typenameNLExpr::value_type, typenameVFExpr::value_type, Result >::second_argument_type >::type A2
void operator()(const self_type &obj, const geo_basic< float_type, M > &omega_K, const geo_element &K, Eigen::Matrix< Result, Eigen::Dynamic, Eigen::Dynamic > &value) const
void operator()(const This &obj, const geo_basic< float_type, M > &omega_K, const geo_element &K, Eigen::Matrix< Result, Eigen::Dynamic, Eigen::Dynamic > &value) const
static space_constant::valued_type valued_tag(space_constant::valued_type, space_constant::valued_type)
Definition: expression.h:422
static space_constant::valued_type valued_tag(space_constant::valued_type)
Definition: expression.h:200
void element_initialize_on_side(const This &obj, const geo_element &K, const side_information_type &sid) const
const point_basic< scalar_type > & get_nl_value(const This &obj, size_type q) const
void element_initialize_on_side(const This &obj, const geo_element &K, const side_information_type &sid) const
const tensor3_basic< scalar_type > & get_nl_value(const This &obj, size_type q) const
void element_initialize_on_side(const This &obj, const geo_element &K, const side_information_type &sid) const
const tensor4_basic< scalar_type > & get_nl_value(const This &obj, size_type q) const
const tensor_basic< scalar_type > & get_nl_value(const This &obj, size_type q) const
void element_initialize_on_side(const This &obj, const geo_element &K, const side_information_type &sid) const
void element_initialize_on_side(const This &obj, const geo_element &K, const side_information_type &sid) const
const scalar_type & get_nl_value(const This &obj, size_type q) const
void element_initialize(const This &obj, const geo_element &K) const
void element_initialize_on_side(const This &obj, const geo_element &K, const side_information_type &sid) const
Arg1 get_nl_value(const This &obj, size_type q) const
void element_initialize(const This &obj, const geo_element &K) const
helper for generic field value_type: T, point_basic<T> or tensor_basic<T>