Rheolef  7.2
an efficient C++ finite element environment
form_expr_variational.h
Go to the documentation of this file.
1 #ifndef _RHEOLEF_FORM_EXPR_VARIATIONAL_H
2 #define _RHEOLEF_FORM_EXPR_VARIATIONAL_H
23 //
24 // variational expressions are used for form assembly
25 //
26 // author: Pierre.Saramito@imag.fr
27 //
28 // date: 23 september 2015
29 //
30 // Notes: use template expressions and SFINAE techniques
31 //
32 // SUMMARY:
33 // 1. concept
34 // 2. unary function
35 // 2.1. unary node
36 // 2.2. unary calls
37 // 3. binary operators +- between two variational forms
38 // 3.1. binary node
39 // 3.2. binary calls
40 // 4. binary operators * between two variational fields that returns a form
41 // 4.1. unary node
42 // 4.2. unary calls
43 // 5. binary operators */ between a variational form and a nonlinear expr
44 // 5.1. binary node
45 // 5.2. binary calls
46 // 6. binary operators */ between a variational form and a constant
47 //
48 #include "rheolef/field_expr_variational.h"
49 
50 namespace rheolef {
51 
52 // -------------------------------------------------------------------
53 // 1. concept
54 // -------------------------------------------------------------------
55 namespace details {
56 
57 // Define a trait type for detecting form expression valid arguments
58 template<class T> struct is_form_expr_v2_variational_arg : std::false_type {};
59 
60 } // namespace details
61 // ---------------------------------------------------------------------------
62 // 2. unary function
63 // example: -(u*v), 2*(u*v), (u*v)/2
64 // ---------------------------------------------------------------------------
65 // 2.1. unary node
66 // ---------------------------------------------------------------------------
67 namespace details {
68 
69 template<class UnaryFunction, class Expr>
71 public:
72 // typedefs:
73 
75  typedef typename Expr::memory_type memory_type;
77  typename Expr::value_type>::type result_hint;
79  typename Expr::value_type
84  typedef typename Expr::vf_tag_type vf_tag_type;
90  typedef typename Expr::maybe_symmetric::type maybe_symmetric;
91 
93 
94 // alocators:
95 
96  form_expr_v2_variational_unary (const UnaryFunction& f, const Expr& expr)
97  : _f(f), _expr(expr) {}
98 
100  : _f(x._f), _expr(x._expr) {}
101 
104  _f = x._f; _expr = x._expr; return *this; }
105 
106 // accessors:
107 
108  const space_type& get_trial_space() const { return _expr.get_trial_space(); }
109  const space_type& get_test_space() const { return _expr.get_test_space(); }
110  size_type n_derivative() const { return _expr.n_derivative(); }
111 
112 // mutable modifiers:
113 
114  void initialize (const piola_on_pointset<float_type>& pops, const integrate_option& iopt) {
115  _expr.initialize (pops, iopt);
116  }
118  _expr.initialize (gh, pops, iopt);
119  }
120  template<class Result>
121  void evaluate (
122  const geo_basic<float_type,memory_type>& omega_K,
123  const geo_element& K,
124  Eigen::Tensor<Result,3>& value) const
125  {
126  // Arg=Result=float_type for f in general: elementary matrix,
127  // so use Result-valued "value" for arg of f:
128  _expr.evaluate (omega_K, K, value);
129  for (size_type i = 0, ni = value.dimension(0); i < ni; ++i) {
130  for (size_type j = 0, nj = value.dimension(1); j < nj; ++j) {
131  for (size_type k = 0, nk = value.dimension(2); k < nk; ++k) {
132  value(i,j,k) = _f (value(i,j,k));
133  }}}
134  }
135  template<class Result>
136  bool valued_check() const {
137  typedef Result A1;
138  if (! is_undeterminated<A1>::value) return _expr.template valued_check<A1>();
139  return true;
140  }
141 protected:
142 // data:
143  UnaryFunction _f;
144  Expr _expr;
145 };
146 template<class F, class Expr> struct is_form_expr_v2_variational_arg <form_expr_v2_variational_unary<F,Expr> > : std::true_type {};
147 
148 } // namespace details
149 // ---------------------------------------------------------------------------
150 // 2.2. unary calls
151 // ---------------------------------------------------------------------------
152 
153 #define _RHEOLEF_make_form_expr_v2_variational_unary(FUNCTION,FUNCTOR) \
154 template<class Expr> \
155 inline \
156 typename \
157 std::enable_if< \
158  details::is_form_expr_v2_variational_arg<Expr>::value \
159  ,details::form_expr_v2_variational_unary< \
160  FUNCTOR \
161  ,Expr \
162  > \
163 >::type \
164 FUNCTION (const Expr& expr) \
165 { \
166  return details::form_expr_v2_variational_unary <FUNCTOR,Expr> (FUNCTOR(), expr); \
167 }
168 
171 #undef _RHEOLEF_make_form_expr_v2_variational_unary
172 
173 // ---------------------------------------------------------------------------
174 // 3. binary operators +- between two variational forms
175 // ---------------------------------------------------------------------------
176 // 3.1. binary node
177 // ---------------------------------------------------------------------------
178 // example: operator+ between two forms as in
179 // (u*v) + dot(grad(u),grad(v))
180 
181 namespace details {
182 
183 template<class BinaryFunction, class Expr1, class Expr2>
185 public:
186 // typedefs:
187 
192  typename Expr1::value_type
193  ,typename Expr2::value_type>::type result_hint;
195  typename Expr1::value_type
196  ,typename Expr2::value_type
200  typedef space_basic<scalar_type,memory_type> space_type; // TODO: deduce from Exprs
201  typedef typename details::bf_vf_tag<BinaryFunction,
202  typename Expr1::vf_tag_type,
203  typename Expr2::vf_tag_type>::type vf_tag_type;
207  typedef form_expr_v2_variational_binary<BinaryFunction,typename Expr1::dual_self_type,
208  typename Expr2::dual_self_type>
210  typedef typename and_type<typename Expr1::maybe_symmetric::type,
211  typename Expr2::maybe_symmetric::type>::type
213 
215 
216 // alocators:
217 
219  const Expr1& expr1,
220  const Expr2& expr2)
221  : _f(f), _expr1(expr1), _expr2(expr2) {}
222 
224  : _f(x._f), _expr1(x._expr1), _expr2(x._expr2) {}
225 
228  { _f = x._f; _expr1 = x._expr1; _expr2 = x._expr2; return *this; }
229 
230 // accessors:
231 
232  const space_type& get_trial_space() const { return _expr1.get_trial_space(); }
233  const space_type& get_test_space() const { return _expr1.get_test_space(); }
234  size_type n_derivative() const { return std::min(_expr1.n_derivative(), _expr2.n_derivative()); }
235 
236 // mutable modifiers:
237 
238  // TODO: at init, check that exp1 & expr2 has the same test & trial spaces
239  void initialize (const piola_on_pointset<float_type>& pops, const integrate_option& iopt) {
240  _expr1.initialize (pops, iopt);
241  _expr2.initialize (pops, iopt);
242  }
244  _expr1.initialize (gh, pops, iopt);
245  _expr2.initialize (gh, pops, iopt);
246  }
247  // Result is float_type in general: elementary matrix
248  template<class Result>
249  void evaluate (
250  const geo_basic<float_type,memory_type>& omega_K,
251  const geo_element& K,
252  Eigen::Tensor<Result,3>& value) const
253  {
254  // for f=operator+ => sum of two elementary matrix of the same type
255  // TODO: otherwise Result and 2 could be obtained from the hint<> helper
256  typedef Result A1;
257  typedef Result A2;
258  Eigen::Tensor<float_type,3> value1; _expr1.evaluate (omega_K, K, value1);
259  Eigen::Tensor<float_type,3> value2; _expr2.evaluate (omega_K, K, value2);
260  value.resize (value1.dimension(0), value1.dimension(1), value1.dimension(2));
261  check_macro (value1.dimension(0) == value2.dimension(0) &&
262  value1.dimension(1) == value2.dimension(1) &&
263  value1.dimension(2) == value2.dimension(2),
264  "invalid sizes value1("
265  << value1.dimension(0) <<","<< value1.dimension(1) <<","<< value1.dimension(2) << ") and value2("
266  << value2.dimension(0) <<","<< value2.dimension(1) <<","<< value2.dimension(2) <<")");
267  for (size_type i = 0, ni = value.dimension(0); i < ni; ++i) {
268  for (size_type j = 0, nj = value.dimension(1); j < nj; ++j) {
269  for (size_type k = 0, nk = value.dimension(2); k < nk; ++k) {
270  value(i,j,k) = _f (value1(i,j,k), value2(i,j,k));
271  }}}
272  }
273  template<class Result>
275  const geo_basic<float_type,memory_type>& omega_K,
276  const geo_element& K,
277  const side_information_type& sid,
278  Eigen::Tensor<Result,3>& value) const
279  {
280  // for f=operator+ => sum of two elementary matrix of the same type
281  // TODO: otherwise Result and 2 could be obtained from the hint<> helper
282  typedef Result A1;
283  typedef Result A2;
284  Eigen::Tensor<float_type,3> value1; _expr1.evaluate_on_side (omega_K, K, sid, value1);
285  Eigen::Tensor<float_type,3> value2; _expr2.evaluate_on_side (omega_K, K, sid, value2);
286  value.resize (value1.dimension(0), value1.dimension(1), value1.dimension(2));
287  check_macro (value1.dimension(0) == value2.dimension(0) &&
288  value1.dimension(1) == value2.dimension(1) &&
289  value1.dimension(2) == value2.dimension(2),
290  "invalid sizes value1("
291  << value1.dimension(0) <<","<< value1.dimension(1) <<","<< value1.dimension(2) << ") and value2("
292  << value2.dimension(0) <<","<< value2.dimension(1) <<","<< value2.dimension(2) <<")");
293  for (size_type i = 0, ni = value.dimension(0); i < ni; ++i) {
294  for (size_type j = 0, nj = value.dimension(1); j < nj; ++j) {
295  for (size_type k = 0, nk = value.dimension(2); k < nk; ++k) {
296  value(i,j,k) = _f (value1(i,j,k), value2(i,j,k));
297  }}}
298  }
299  template<class Result>
300  bool valued_check() const {
301  typedef Result A1;
302  typedef Result A2;
303  bool status = true;
304  if (! is_undeterminated<A1>::value) status &= _expr1.template valued_check<A1>();
305  if (! is_undeterminated<A2>::value) status &= _expr2.template valued_check<A2>();
306  return status;
307  }
308 protected:
309 // data:
312  Expr2 _expr2;
313 };
314 template<class F, class Expr1, class Expr2> struct is_form_expr_v2_variational_arg <form_expr_v2_variational_binary<F,Expr1,Expr2> > : std::true_type {};
315 
316 } // namespace details
317 
318 // ---------------------------------------------------------------------------
319 // 3.2. binary calls
320 // ---------------------------------------------------------------------------
321 
322 #define _RHEOLEF_form_expr_v2_variational_binary(FUNCTION,FUNCTOR) \
323 template <class Expr1, class Expr2> \
324 inline \
325 typename \
326 std::enable_if< \
327  details::is_form_expr_v2_variational_arg <Expr1>::value \
328  && details::is_form_expr_v2_variational_arg <Expr2>::value \
329  ,details::form_expr_v2_variational_binary< \
330  FUNCTOR \
331  ,Expr1 \
332  ,Expr2 \
333  > \
334 >::type \
335 FUNCTION (const Expr1& expr1, const Expr2& expr2) \
336 { \
337  return details::form_expr_v2_variational_binary \
338  <FUNCTOR, Expr1, Expr2> \
339  (FUNCTOR(), expr1, expr2); \
340 }
341 
344 
345 #undef _RHEOLEF_form_expr_v2_variational_binary
346 
347 // ---------------------------------------------------------------------------
348 // 4. binary operators * between two variational fields that returns a form
349 // example: integrate(u*v)
350 // ---------------------------------------------------------------------------
351 // 4.1. binary node
352 // ---------------------------------------------------------------------------
353 namespace details {
354 
355 template<class BinaryFunction, class Expr1, class Expr2>
357 public:
358 // typedefs:
359 
364  typename Expr1::value_type
365  ,typename Expr2::value_type>::type result_hint;
367  typename Expr1::value_type
368  ,typename Expr2::value_type
372  typedef space_basic<scalar_type,memory_type> space_type; // TODO: deduce from Exprs
373  typedef typename details::bf_vf_tag<BinaryFunction,
374  typename Expr1::vf_tag_type,
375  typename Expr2::vf_tag_type>::type vf_tag_type;
382  typename details::is_equal<
383  Expr1,
384  typename Expr2::dual_self_type>::type>::type
386 
388 
389 // alocators:
390 
392  const Expr1& expr1,
393  const Expr2& expr2)
394  : _f(f), _expr1(expr1), _expr2(expr2) {}
397  _f = x._f; _expr1 = x._expr1; _expr2 = x._expr2; return *this; }
398 
399 // accessors:
400 
401  const space_type& get_test_space() const {
403  _expr1.get_vf_space() : _expr2.get_vf_space();
404  }
405  const space_type& get_trial_space() const {
407  _expr1.get_vf_space() : _expr2.get_vf_space();
408  }
409 
410  size_type n_derivative() const { return _expr1.n_derivative() + _expr2.n_derivative(); }
411 
412 // mutable modifiers:
413 
414  void initialize (
415  const piola_on_pointset<float_type>& pops,
416  const integrate_option& iopt)
417  {
418  _expr1.initialize (pops, iopt);
419  _expr2.initialize (pops, iopt);
420  }
421  void initialize (
423  const piola_on_pointset<float_type>& pops,
424  const integrate_option& iopt)
425  {
426  _expr1.initialize (gh, pops, iopt);
427  _expr2.initialize (gh, pops, iopt);
428  }
429  // TODO: DVT_EXPR_SPARSE_MATRIX return array[q] of SparseMatrix(i,j)
430  template<class ValueType, class Arg1, class Arg2>
432  const geo_basic<float_type,memory_type>& omega_K,
433  const geo_element& K,
434  Eigen::Tensor<ValueType,3>& value) const
435  {
436  typedef long int eig_idx_t;
438  // expr1 is trial and expr2 is test
439  Eigen::Matrix<Arg2,Eigen::Dynamic,Eigen::Dynamic> v_test; _expr2.evaluate (omega_K, K, v_test);
440  Eigen::Matrix<Arg1,Eigen::Dynamic,Eigen::Dynamic> u_trial; _expr1.evaluate (omega_K, K, u_trial);
441  check_macro(u_trial.rows() == v_test.rows(), "mult: invalid sizes u_trial("
442  <<u_trial.rows()<<","<<u_trial.cols() <<") and v_test("
443  <<v_test.rows() <<","<<v_test.cols()<<")");
444  eig_idx_t nq = u_trial.rows();
445  eig_idx_t ni = v_test.cols();
446  eig_idx_t nj = u_trial.cols();
447  value.resize(nq,ni,nj);
448  for (eig_idx_t q = 0; q < nq; ++q) {
449  for (eig_idx_t i = 0; i < ni; ++i) {
450  for (eig_idx_t j = 0; j < nj; ++j) {
451  value(q,i,j) = _f (u_trial(q,j), v_test(q,i));
452  }}}
453  } else {
454  // expr2 is trial and expr1 is test
455  Eigen::Matrix<Arg2,Eigen::Dynamic,Eigen::Dynamic> v_test; _expr1.evaluate (omega_K, K, v_test);
456  Eigen::Matrix<Arg1,Eigen::Dynamic,Eigen::Dynamic> u_trial; _expr2.evaluate (omega_K, K, u_trial);
457  check_macro(u_trial.rows() == v_test.rows(), "binary: invalid sizes");
458  eig_idx_t nq = u_trial.rows();
459  eig_idx_t ni = v_test.cols();
460  eig_idx_t nj = u_trial.cols();
461  value.resize(nq,ni,nj);
462  for (eig_idx_t q = 0; q < nq; ++q) {
463  for (eig_idx_t i = 0; i < ni; ++i) {
464  for (eig_idx_t j = 0; j < nj; ++j) {
465  value(q,i,j) = _f (u_trial(q,j), v_test(q,i));
466  }}}
467  }
468  }
469  template<class ValueType, class Arg1, class Arg2>
471  const geo_basic<float_type,memory_type>& omega_K,
472  const geo_element& K,
473  const side_information_type& sid,
474  Eigen::Tensor<ValueType,3>& value) const
475  {
476  typedef long int eig_idx_t;
477  bool do_local_component_assembly = true;
479  // expr1 is trial and expr2 is test
480  Eigen::Matrix<Arg2,Eigen::Dynamic,Eigen::Dynamic> v_test; _expr2.evaluate_on_side (omega_K, K, sid, v_test, do_local_component_assembly);
481  Eigen::Matrix<Arg1,Eigen::Dynamic,Eigen::Dynamic> u_trial; _expr1.evaluate_on_side (omega_K, K, sid, u_trial, do_local_component_assembly);
482  check_macro(u_trial.rows() == v_test.rows(), "invalid sizes u_trial("
483  <<u_trial.rows()<<","<<u_trial.cols() <<") and v_test("
484  <<v_test.rows() <<","<<v_test.cols()<<")");
485  eig_idx_t nq = u_trial.rows();
486  eig_idx_t ni = v_test.cols();
487  eig_idx_t nj = u_trial.cols();
488  value.resize(nq,ni,nj);
489  for (eig_idx_t q = 0; q < nq; ++q) {
490  for (eig_idx_t i = 0; i < ni; ++i) {
491  for (eig_idx_t j = 0; j < nj; ++j) {
492  value(q,i,j) = _f (u_trial(q,j), v_test(q,i));
493  }}}
494  } else {
495  // expr2 is trial and expr1 is test
496  Eigen::Matrix<Arg2,Eigen::Dynamic,Eigen::Dynamic> v_test; _expr1.evaluate_on_side (omega_K, K, sid, v_test, do_local_component_assembly);
497  Eigen::Matrix<Arg1,Eigen::Dynamic,Eigen::Dynamic> u_trial; _expr2.evaluate_on_side (omega_K, K, sid, u_trial, do_local_component_assembly);
498  check_macro(u_trial.rows() == v_test.rows(), "binary: invalid sizes");
499  eig_idx_t nq = u_trial.rows();
500  eig_idx_t ni = v_test.cols();
501  eig_idx_t nj = u_trial.cols();
502  value.resize(nq,ni,nj);
503  for (eig_idx_t q = 0; q < nq; ++q) {
504  for (eig_idx_t i = 0; i < ni; ++i) {
505  for (eig_idx_t j = 0; j < nj; ++j) {
506  value(q,i,j) = _f (u_trial(q,j), v_test(q,i));
507  }}}
508  }
509  }
510  // when both args are defined at compile time:
511  template<class This, class ValueType,
512  class Arg1, space_constant::valued_type Arg1Tag,
513  class Arg2, space_constant::valued_type Arg2Tag>
515  typedef typename This::float_type float_type;
516  typedef typename This::memory_type memory_type;
518  const This& obj,
519  const geo_basic<float_type,memory_type>& omega_K,
520  const geo_element& K,
521  Eigen::Tensor<ValueType,3>& value) const
522  {
523  obj.template evaluate_internal<ValueType, Arg1, Arg2> (omega_K, K, value);
524  }
526  const This& obj,
527  const geo_basic<float_type,memory_type>& omega_K,
528  const geo_element& K,
529  const side_information_type& sid,
530  Eigen::Tensor<ValueType,3>& value) const
531  {
532  obj.template evaluate_on_side_internal<ValueType, Arg1, Arg2> (omega_K, K, sid, value);
533  }
534  };
535  template<class ValueType>
536  void evaluate (
537  const geo_basic<float_type,memory_type>& omega_K,
538  const geo_element& K,
539  Eigen::Tensor<ValueType,3>& value) const
540  {
542  typename Expr1::value_type
543  ,typename Expr2::value_type
544  ,ValueType>::first_argument_type first_argument_type;
546  typename Expr1::value_type
547  ,typename Expr2::value_type
548  ,ValueType>::second_argument_type second_argument_type;
551  evaluate_switch <self_type, ValueType,
552  first_argument_type, first_argument_tag,
553  second_argument_type, second_argument_tag> eval;
554  eval (*this, omega_K, K, value);
555  }
556  template<class ValueType>
558  const geo_basic<float_type,memory_type>& omega_K,
559  const geo_element& K,
560  const side_information_type& sid,
561  Eigen::Tensor<ValueType,3>& value) const
562  {
564  typename Expr1::value_type
565  ,typename Expr2::value_type
566  ,ValueType>::first_argument_type first_argument_type;
568  typename Expr1::value_type
569  ,typename Expr2::value_type
570  ,ValueType>::second_argument_type second_argument_type;
573  evaluate_switch <self_type, ValueType,
574  first_argument_type, first_argument_tag,
575  second_argument_type, second_argument_tag> eval;
576  eval (*this, omega_K, K, sid, value);
577  }
578  template<class ValueType>
579  bool valued_check() const {
581  typename Expr1::value_type
582  ,typename Expr2::value_type
583  ,ValueType>::first_argument_type A1;
585  typename Expr1::value_type
586  ,typename Expr2::value_type
587  ,ValueType>::second_argument_type A2;
588  if (! is_undeterminated<A1>::value) _expr1.template valued_check<A1>();
589  if (! is_undeterminated<A2>::value) _expr2.template valued_check<A2>();
590  return true;
591  }
592 protected:
593 // data:
596  Expr2 _expr2;
597 };
598 template<class F, class Expr1, class Expr2> struct is_form_expr_v2_variational_arg <form_expr_v2_variational_binary_field<F,Expr1,Expr2> > : std::true_type {};
599 
600 } // namespace details
601 // ---------------------------------------------------------------------------
602 // 4.2. binary calls
603 // ---------------------------------------------------------------------------
604 namespace details {
605 
606 template<class Expr1, class Expr2, class Sfinae = void>
607 struct is_form_expr_v2_variational_binary_field : std::false_type {};
608 
609 template <class Expr1, class Expr2>
611  Expr1
612  ,Expr2
613  ,typename
614  std::enable_if<
615  is_field_expr_v2_variational_arg<Expr1>::value
616  && is_field_expr_v2_variational_arg<Expr2>::value
617  >::type
618 >
619 : and_type<
620  is_field_expr_v2_variational_arg<Expr1>
621  ,is_field_expr_v2_variational_arg<Expr2>
622  ,std::is_same <
623  typename Expr1::vf_tag_type
624  ,typename dual_vf_tag<typename Expr2::vf_tag_type>::type
625  >
626  >
627 {};
628 
629 } // namespace details
630 
631 #define _RHEOLEF_form_expr_v2_variational_binary_field(FUNCTION,FUNCTOR) \
632 template <class Expr1, class Expr2> \
633 inline \
634 typename \
635 std::enable_if< \
636  details::is_form_expr_v2_variational_binary_field <Expr1,Expr2>::value \
637  ,details::form_expr_v2_variational_binary_field< \
638  FUNCTOR \
639  ,Expr1 \
640  ,Expr2 \
641  > \
642 >::type \
643 FUNCTION (const Expr1& expr1, const Expr2& expr2) \
644 { \
645  return details::form_expr_v2_variational_binary_field \
646  <FUNCTOR, Expr1, Expr2> \
647  (FUNCTOR(), expr1, expr2); \
648 }
649 
654 
655 #undef _RHEOLEF_form_expr_v2_variational_binary_field
656 
657 // ---------------------------------------------------------------------------
658 // 5. binary operators */ between a variational form and a nonlinear expr
659 // ---------------------------------------------------------------------------
660 // 5.1. binary node
661 // ---------------------------------------------------------------------------
662 // example: integrate(eta_h*(u*v))
663 
664 namespace details {
665 
666 template<class BinaryFunction, class NLExpr, class VFExpr>
668 public:
669 // typedefs:
670 
675  typename NLExpr::value_type
676  ,typename VFExpr::value_type>::type result_hint;
678  typename NLExpr::value_type
679  ,typename VFExpr::value_type
683  typedef space_basic<scalar_type,memory_type> space_type; // TODO: deduce from Exprs
684  typedef typename details::bf_vf_tag<BinaryFunction,
686  typename VFExpr::vf_tag_type>::type vf_tag_type;
692  typedef typename VFExpr::maybe_symmetric::type maybe_symmetric;
693  // TODO: symmetry: works only when eta_h is scalar
694  // TODO: problem when ddot(eta_h,otimes(u,v)) when eta_h is unsymmetric tensor
695  // and "unsymmetric tensor" is not known at compile time
696 
698 
699 // alocators:
700 
702  const NLExpr& nl_expr,
703  const VFExpr& vf_expr)
704  : _f(f),
705  _nl_expr(nl_expr),
706  _vf_expr(vf_expr)
707  {}
708 
709 // accessors:
710 
711  const space_type& get_trial_space() const { return _vf_expr.get_trial_space(); }
712  const space_type& get_test_space() const { return _vf_expr.get_test_space(); }
713  size_type n_derivative() const { return _vf_expr.n_derivative(); }
714 
715 // mutable modifiers:
716 
717  void initialize (const piola_on_pointset<float_type>& pops, const integrate_option& iopt) {
718  _nl_expr.initialize (pops, iopt);
719  _vf_expr.initialize (pops, iopt);
720  }
722  _nl_expr.initialize ( pops, iopt);
723  _vf_expr.initialize (gh, pops, iopt);
724  }
725  // ---------------------------------------------
726  // element initialize: evaluate nl_expr
727  // ---------------------------------------------
728  template<class Result>
729  void evaluate (
730  const geo_basic<float_type,memory_type>& omega_K,
731  const geo_element& K,
732  Eigen::Tensor<Result,3>& value) const
733  {
734  typedef Result Arg1; // TODO: switch
735  typedef Result Arg2; // TODO: switch ; is float_type in general, as elementary matrix
736  Eigen::Matrix<Arg1,Eigen::Dynamic,1> value1; _nl_expr.evaluate (omega_K, K, value1);
737  Eigen::Tensor<float_type,3> value2; _vf_expr.evaluate (omega_K, K, value2);
738  check_macro (value1.size() == value2.dimension(0), "invalid sizes");
739  size_type loc_nnod = value2.dimension(0);
740  size_type loc_ndof1 = value2.dimension(1);
741  size_type loc_ndof2 = value2.dimension(2);
742  value.resize(loc_nnod,loc_ndof1,loc_ndof2);
743  for (size_type loc_inod = 0; loc_inod < loc_nnod; ++loc_inod) {
744  for (size_type loc_jdof = 0; loc_jdof < loc_ndof1; ++loc_jdof) {
745  for (size_type loc_kdof = 0; loc_kdof < loc_ndof2; ++loc_kdof) {
746  value(loc_inod,loc_jdof,loc_kdof)
747  = _f (value1(loc_inod,loc_jdof), value2(loc_inod,loc_jdof,loc_kdof));
748  }}}
749  }
750  template<class Result>
751  bool valued_check() const {
752  typedef Result A1;
753  typedef Result A2;
754  bool status = true;
755  if (! is_undeterminated<A1>::value) status &= _nl_expr.template valued_check<A1>();
756  if (! is_undeterminated<A2>::value) status &= _vf_expr.template valued_check<A2>();
757  return status;
758  }
759 //protected:
760 // data:
762  NLExpr _nl_expr;
763  VFExpr _vf_expr;
764 };
765 template<class F, class Expr1, class Expr2> struct is_form_expr_v2_variational_arg <form_expr_v2_variational_binary_binded<F,Expr1,Expr2> > : std::true_type {};
766 
767 } // namespace details
768 
769 // ---------------------------------------------------------------------------
770 // 5.2. binary calls
771 // ---------------------------------------------------------------------------
772 namespace details {
773 
774 template<class Expr1, class Expr2, class Sfinae = void>
776 
777 template<class Expr1, class Expr2>
779  Expr1
780  ,Expr2
781  ,typename
782  std::enable_if<
783  is_field_expr_v2_nonlinear_arg <Expr1>::value
784  && ! is_rheolef_arithmetic <Expr1>::value
785  && is_form_expr_v2_variational_arg<Expr2>::value
786  >::type
787 >
788 : std::true_type
789 {};
790 
791 template<class Expr1, class Expr2>
794 
795 } // namespace details
796 
797 #define _RHEOLEF_make_form_expr_v2_variational_binary_operator_multiplies_divides_left(FUNCTION,FUNCTOR) \
798 template<class Expr1, class Expr2> \
799 inline \
800 typename \
801 std::enable_if< \
802  details::is_form_expr_v2_variational_binary_multiplies_divides_left <Expr1,Expr2>::value \
803  ,details::form_expr_v2_variational_binary_binded< \
804  FUNCTOR \
805  ,typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr1>::type \
806  ,Expr2 /* vf */ \
807  > \
808 >::type \
809 FUNCTION (const Expr1& expr1, const Expr2& expr2) \
810 { \
811  typedef typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr1>::type wrap1_t; \
812  return details::form_expr_v2_variational_binary_binded \
813  <FUNCTOR, wrap1_t, Expr2> \
814  (FUNCTOR(), wrap1_t(expr1), expr2); \
815 }
816 
817 #define _RHEOLEF_make_form_expr_v2_variational_binary_operator_multiplies_divides_right(FUNCTION,FUNCTOR) \
818 template<class Expr1, class Expr2> \
819 inline \
820 typename \
821 std::enable_if< \
822  details::is_form_expr_v2_variational_binary_multiplies_divides_right <Expr1,Expr2>::value \
823  ,details::form_expr_v2_variational_binary_binded< \
824  details::swapper<FUNCTOR> \
825  , typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr2>::type \
826  ,Expr1 /* vf */ \
827  > \
828 >::type \
829 FUNCTION (const Expr1& expr1, const Expr2& expr2) \
830 { \
831  typedef typename details::field_expr_v2_nonlinear_terminal_wrapper_traits<Expr2>::type wrap2_t; \
832  return details::form_expr_v2_variational_binary_binded \
833  <details::swapper<FUNCTOR>, wrap2_t, Expr1> \
834  (details::swapper<FUNCTOR>(FUNCTOR()), wrap2_t(expr2), expr1); \
835 }
836 #define _RHEOLEF_make_form_expr_v2_variational_binary_operator_multiplies_divides(FUNCTION,FUNCTOR) \
837  _RHEOLEF_make_form_expr_v2_variational_binary_operator_multiplies_divides_left (FUNCTION,FUNCTOR) \
838  _RHEOLEF_make_form_expr_v2_variational_binary_operator_multiplies_divides_right (FUNCTION,FUNCTOR)
839 
845 #undef _RHEOLEF_make_form_expr_v2_variational_binary_operator_multiplies_divides_left
846 #undef _RHEOLEF_make_form_expr_v2_variational_binary_operator_multiplies_divides_right
847 #undef _RHEOLEF_make_form_expr_v2_variational_binary_operator_multiplies_divides
848 
849 // ---------------------------------------------------------------------------
850 // 6. binary operators */ between a variational form and a constant
851 // ---------------------------------------------------------------------------
852 namespace details {
853 
854 template<class Expr1, class Expr2, class Sfinae = void>
855 struct is_form_expr_v2_variational_binary_multiplies_divides_constant_left : std::false_type {};
856 
857 template<class Expr1, class Expr2>
858 struct is_form_expr_v2_variational_binary_multiplies_divides_constant_left <
859  Expr1
860  ,Expr2
861  ,typename
862  std::enable_if<
863  is_rheolef_arithmetic <Expr1>::value
864  && is_form_expr_v2_variational_arg<Expr2>::value
865  >::type
866 >
867 : std::true_type
868 {};
869 
870 template<class Expr1, class Expr2>
871 struct is_form_expr_v2_variational_binary_multiplies_divides_constant_right
872 : is_form_expr_v2_variational_binary_multiplies_divides_constant_left <Expr2,Expr1> {};
873 
874 } // namespace details
875 
876 #define _RHEOLEF_make_form_expr_v2_variational_binary_operator_multiplies_divides_constant_left(FUNCTION,FUNCTOR) \
877 template<class Expr1, class Expr2> \
878 inline \
879 typename \
880 std::enable_if< \
881  details::is_form_expr_v2_variational_binary_multiplies_divides_constant_left <Expr1,Expr2>::value \
882  ,details::form_expr_v2_variational_unary< \
883  details::binder_first <FUNCTOR, Expr1> \
884  ,Expr2 /* vf */ \
885  > \
886 >::type \
887 FUNCTION (const Expr1& expr1, const Expr2& expr2) \
888 { \
889  return details::form_expr_v2_variational_unary \
890  <details::binder_first <FUNCTOR,Expr1>, Expr2> \
891  (details::binder_first <FUNCTOR,Expr1> (FUNCTOR(), expr1), expr2); \
892 }
893 
894 #define _RHEOLEF_make_form_expr_v2_variational_binary_operator_multiplies_divides_constant_right(FUNCTION,FUNCTOR) \
895 template<class Expr1, class Expr2> \
896 inline \
897 typename \
898 std::enable_if< \
899  details::is_form_expr_v2_variational_binary_multiplies_divides_constant_right <Expr1,Expr2>::value \
900  ,details::form_expr_v2_variational_unary< \
901  details::binder_second <FUNCTOR, Expr2> \
902  ,Expr1 /* vf */ \
903  > \
904 >::type \
905 FUNCTION (const Expr1& expr1, const Expr2& expr2) \
906 { \
907  return details::form_expr_v2_variational_unary \
908  <details::binder_second <FUNCTOR,Expr2>, Expr1> \
909  (details::binder_second <FUNCTOR,Expr2> (FUNCTOR(), expr2), expr1); \
910 }
911 
912 #define _RHEOLEF_make_form_expr_v2_variational_binary_operator_multiplies_divides_constant(FUNCTION,FUNCTOR) \
913  _RHEOLEF_make_form_expr_v2_variational_binary_operator_multiplies_divides_constant_left (FUNCTION,FUNCTOR) \
914  _RHEOLEF_make_form_expr_v2_variational_binary_operator_multiplies_divides_constant_right (FUNCTION,FUNCTOR)
915 
916 
922 
923 #undef _RHEOLEF_make_form_expr_v2_variational_binary_operator_multiplies_divides_constant_right
924 #undef _RHEOLEF_make_form_expr_v2_variational_binary_operator_multiplies_divides_constant_left
925 #undef _RHEOLEF_make_form_expr_v2_variational_binary_operator_multiplies_divides_constant
926 
927 } // namespace rheolef
928 #endif // _RHEOLEF_FORM_EXPR_VARIATIONAL_H
field gh(Float epsilon, Float t, const field &uh, const test &v)
details::generic_binary_traits< BinaryFunction >::template result_hint< typename NLExpr::value_type,typename VFExpr::value_type >::type result_hint
form_expr_v2_variational_binary_binded< BinaryFunction, NLExpr, typename VFExpr::dual_self_type > dual_self_type
details::bf_vf_tag< BinaryFunction, details::vf_tag_00, typename VFExpr::vf_tag_type >::type vf_tag_type
void initialize(const piola_on_pointset< float_type > &pops, const integrate_option &iopt)
form_expr_v2_variational_binary_binded< BinaryFunction, NLExpr, VFExpr > self_type
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 hint< typename NLExpr::value_type,typename VFExpr::value_type,result_hint >::result_type value_type
form_expr_v2_variational_binary_binded(const BinaryFunction &f, const NLExpr &nl_expr, const VFExpr &vf_expr)
static const space_constant::valued_type valued_hint
details::dual_vf_tag< vf_tag_type >::type vf_dual_tag_type
void evaluate(const geo_basic< float_type, memory_type > &omega_K, const geo_element &K, Eigen::Tensor< Result, 3 > &value) const
promote_memory< typename NLExpr::memory_type, typename VFExpr::memory_type >::type memory_type
details::generic_binary_traits< BinaryFunction >::template result_hint< typename Expr1::value_type,typename Expr2::value_type >::type result_hint
form_expr_v2_variational_binary_field(const BinaryFunction &f, const Expr1 &expr1, const Expr2 &expr2)
promote_memory< typename Expr1::memory_type, typename Expr2::memory_type >::type memory_type
void initialize(const piola_on_pointset< float_type > &pops, const integrate_option &iopt)
void evaluate(const geo_basic< float_type, memory_type > &omega_K, const geo_element &K, Eigen::Tensor< ValueType, 3 > &value) const
void initialize(const band_basic< float_type, memory_type > &gh, const piola_on_pointset< float_type > &pops, const integrate_option &iopt)
void evaluate_on_side_internal(const geo_basic< float_type, memory_type > &omega_K, const geo_element &K, const side_information_type &sid, Eigen::Tensor< ValueType, 3 > &value) const
details::generic_binary_traits< BinaryFunction >::template hint< typename Expr1::value_type,typename Expr2::value_type,result_hint >::result_type value_type
void evaluate_internal(const geo_basic< float_type, memory_type > &omega_K, const geo_element &K, Eigen::Tensor< ValueType, 3 > &value) const
details::bf_vf_tag< BinaryFunction, typename Expr1::vf_tag_type, typename Expr2::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
form_expr_v2_variational_binary_field< BinaryFunction, Expr1, Expr2 > self_type
form_expr_v2_variational_binary_field< BinaryFunction, Expr1, Expr2 > & operator=(const form_expr_v2_variational_binary_field< BinaryFunction, Expr1, Expr2 > &x)
form_expr_v2_variational_binary_field< BinaryFunction, typename Expr1::dual_self_type, typename Expr2::dual_self_type > dual_self_type
void evaluate_on_side(const geo_basic< float_type, memory_type > &omega_K, const geo_element &K, const side_information_type &sid, Eigen::Tensor< ValueType, 3 > &value) const
and_type< typename details::generic_binary_traits< BinaryFunction >::is_symmetric::type, typename details::is_equal< Expr1, typename Expr2::dual_self_type >::type >::type maybe_symmetric
details::generic_binary_traits< BinaryFunction >::template result_hint< typename Expr1::value_type,typename Expr2::value_type >::type result_hint
space_basic< scalar_type, memory_type > space_type
promote_memory< typename Expr1::memory_type, typename Expr2::memory_type >::type memory_type
void initialize(const piola_on_pointset< float_type > &pops, const integrate_option &iopt)
form_expr_v2_variational_binary< BinaryFunction, Expr1, Expr2 > self_type
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 hint< typename Expr1::value_type,typename Expr2::value_type,result_hint >::result_type value_type
form_expr_v2_variational_binary(const form_expr_v2_variational_binary< BinaryFunction, Expr1, Expr2 > &x)
form_expr_v2_variational_binary< BinaryFunction, Expr1, Expr2 > & operator=(const form_expr_v2_variational_binary< BinaryFunction, Expr1, Expr2 > &x)
details::bf_vf_tag< BinaryFunction, typename Expr1::vf_tag_type, typename Expr2::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
form_expr_v2_variational_binary(const BinaryFunction &f, const Expr1 &expr1, const Expr2 &expr2)
void evaluate_on_side(const geo_basic< float_type, memory_type > &omega_K, const geo_element &K, const side_information_type &sid, Eigen::Tensor< Result, 3 > &value) const
form_expr_v2_variational_binary< BinaryFunction, typename Expr1::dual_self_type, typename Expr2::dual_self_type > dual_self_type
and_type< typename Expr1::maybe_symmetric::type, typename Expr2::maybe_symmetric::type >::type maybe_symmetric
void evaluate(const geo_basic< float_type, memory_type > &omega_K, const geo_element &K, Eigen::Tensor< Result, 3 > &value) const
form_expr_v2_variational_unary(const UnaryFunction &f, const Expr &expr)
space_basic< scalar_type, memory_type > space_type
form_expr_v2_variational_unary(const form_expr_v2_variational_unary< UnaryFunction, Expr > &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)
form_expr_v2_variational_unary< UnaryFunction, typename Expr::dual_self_type > dual_self_type
details::generic_unary_traits< UnaryFunction >::template hint< typename Expr::value_type,result_hint >::result_type value_type
static const space_constant::valued_type valued_hint
details::dual_vf_tag< vf_tag_type >::type vf_dual_tag_type
details::generic_unary_traits< UnaryFunction >::template result_hint< typename Expr::value_type >::type result_hint
form_expr_v2_variational_unary< UnaryFunction, Expr > & operator=(const form_expr_v2_variational_unary< UnaryFunction, Expr > &x)
void evaluate(const geo_basic< float_type, memory_type > &omega_K, const geo_element &K, Eigen::Tensor< Result, 3 > &value) const
form_expr_v2_variational_unary< UnaryFunction, Expr > self_type
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
rheolef::std BinaryFunction
rheolef::std value
rheolef::std Expr1
dot(x,y): see the expression page for the full documentation
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_form_expr_v2_variational_binary_operator_multiplies_divides_right(FUNCTION, FUNCTOR)
#define _RHEOLEF_make_form_expr_v2_variational_binary_operator_multiplies_divides_constant_right(FUNCTION, FUNCTOR)
std::pair< std::false_type, std::false_type > vf_tag_00
This file is part of Rheolef.
T dddot(const tensor3_basic< T > &a, const tensor3_basic< T > &b)
Definition: tensor3.cc:94
rheolef::std enable_if ::type dot const Expr1 expr1, const Expr2 expr2 dot(const Expr1 &expr1, const Expr2 &expr2)
dot(x,y): see the expression page for the full documentation
Definition: vec_expr_v2.h:415
_RHEOLEF_form_expr_v2_variational_binary(operator+, details::plus) _RHEOLEF_form_expr_v2_variational_binary(operator-
_RHEOLEF_form_expr_v2_variational_binary_field(operator*, details::multiplies) _RHEOLEF_form_expr_v2_variational_binary_field(dot
_RHEOLEF_make_form_expr_v2_variational_binary_operator_multiplies_divides(operator*, details::multiplies) _RHEOLEF_make_form_expr_v2_variational_binary_operator_multiplies_divides_right(operator/
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
_RHEOLEF_make_form_expr_v2_variational_binary_operator_multiplies_divides_constant(operator*, details::multiplies) _RHEOLEF_make_form_expr_v2_variational_binary_operator_multiplies_divides_constant_right(operator/
_RHEOLEF_make_form_expr_v2_variational_unary(operator+, details::unary_plus) _RHEOLEF_make_form_expr_v2_variational_unary(operator-
Definition: cavity_dg.h:29
void operator()(const This &obj, const geo_basic< float_type, memory_type > &omega_K, const geo_element &K, Eigen::Tensor< ValueType, 3 > &value) const