Rheolef  7.2
an efficient C++ finite element environment
form_lazy_expr.h
Go to the documentation of this file.
1 # ifndef _RHEOLEF_FORM_LAZY_EXPR_H
2 # define _RHEOLEF_FORM_LAZY_EXPR_H
3 //
4 // This file is part of Rheolef.
5 //
6 // Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
7 //
8 // Rheolef is free software; you can redistribute it and/or modify
9 // it under the terms of the GNU General Public License as published by
10 // the Free Software Foundation; either version 2 of the License, or
11 // (at your option) any later version.
12 //
13 // Rheolef is distributed in the hope that it will be useful,
14 // but WITHOUT ANY WARRANTY; without even the implied warranty of
15 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
16 // GNU General Public License for more details.
17 //
18 // You should have received a copy of the GNU General Public License
19 // along with Rheolef; if not, write to the Free Software
20 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
21 //
22 // =========================================================================
23 // form_lazy_expr = result of a lazy integrate() that are combined in exprs
24 // AUTHOR: Pierre.Saramito@imag.fr
25 // DATE: 30 march 1920
26 
27 namespace rheolef {
116 } // namespace rheolef
117 
118 //
119 // SUMMARY: see also "form_lazy_terminal.h"
120 // 3. unary expressions
121 // 3.1. unary_minus & unary_plus
122 // 3.2. invert
123 // 3.3. transpose
124 // 4. binary expressions
125 // 4.1. plus & minus
126 // 4.2. multiply
127 // 5. form allocator
128 //
129 
130 /*
131  Question:
132  What is the difference between
133  field_lazy_xxx and field_expr_xxx ?
134  form_lazy_xxx and form_expr_xxx ?
135  Response:
136  "lazy" means un-assembled matrix and vectors, here fields and forms.
137  The idea is to compute on the fly fields and forms during an iteration
138  on the element of the mesh.
139 
140  exemple 1:
141  wh = uh + vh;
142  1) When both uh and vh are assembled (ie field)
143  the loop is unrolled at the "dof" level:
144  wh(i)= uh(i) + vh(i)
145  => field_expr unroll at the "dof" level
146  2) Here, either uh or vh are un-assembled (ie field_lazy)
147  the loop is unrolled at the element level:
148  wh/K = uh/K + vh/K
149  => field_lazy unrolls at the element level
150 
151  exemple 2:
152  wh = a*uh + b*vh;
153  1) When both uh and vh are assembled (ie field)
154  the loop is unrolled at the "dof" level:
155  wh(i)= sum_j a_ij*uh(j) + sum_k b_ik*vh(k)
156  => field_expr unroll at the dof level
157  2) Here, either uh or vh are un-assembled (ie field_lazy)
158  the loop is unrolled at the element level:
159  wh/K = a/K*uh/K + b/K*vh/K
160  assuming that we use discontiuous elements (eg "Pkd")
161  otherwise a and/or b are not diagonal versus K
162  => field_lazy unrolls at the element level
163 
164  exemple 3:
165  auto a = ...
166  auto b = ...
167  auto c = ...
168  form S = inv(a + b*inv(c)*trans(b));
169  The two imbricated inversions are possible here
170  based on an unassembled form_lazy a, b, c.
171  Otherwise, with a, b, c as forms, this is not possible.
172  Such an idom is typical in the HHO method, see dirichlet_hho.cc
173 
174  IMPLEMENTATION NOTES:
175  * Expressions operates at the elementy level:
176  element matrix are: inverted, transposed, mult, etc.
177  Effective computations are delayed until the full
178  expr is converted to form_basic<T,M>
179  * A base class defines the common methods for all form_lazy_xxx
180  it uses the CRTP C++ idiom:
181  https://en.wikipedia.org/wiki/Curiously_recurring_template_pattern
182  * Some computations are expensive :
183  - for field: O(n^p) with p > 1 where n is the typical element-vector size
184  - for forms: O(n^p) with p > 2 where n is the typical element-matrix size
185  Are converned:
186  - field: a*u, u=integrate(...), u=interpolate(...)
187  - forms: inv(a), a*b, a=integrate(...)
188  For these expensive computations, it is possible to share the result
189  of computations by speciying where is the common sub-expresion:
190  example:
191  auto a1 = inv(m)*b - b*inv(m); // inv(m) is computed 2 times
192  auto inv_m = inv(m); // the common sub-expr
193  auto a2 = inv_m*b - b*inv_m; // inv(m) computed once on the fly
194  requirements:
195  - The initialize() method should be "const" for the smart_pointer
196  to avoid generating a true copy when calling initialize()
197  - This also avoid the true-copy semantic of the imbedded eigen::matrix
198  that stores the shared result
199 
200 
201 */
202 #include "rheolef/form_lazy_terminal.h"
203 
204 namespace rheolef {
205 
206 // -------------------------------------------------------------------
207 // 3. unary expressions
208 // -------------------------------------------------------------------
209 // 3.1. unary_minus & unary_plus
210 // -------------------------------------------------------------------
211 namespace details {
212 
213 template<class Unop, class Expr>
214 class form_lazy_unop: public form_lazy_base<form_lazy_unop<Unop,Expr>> {
215 public :
216 // definitions:
217 
220  using memory_type = typename Expr::memory_type;
221  using scalar_type = typename Expr::scalar_type;
222  using space_type = typename Expr::space_type;
223  using geo_type = typename Expr::geo_type;
226  using matrix_element_type = Eigen::Matrix<scalar_type,Eigen::Dynamic,Eigen::Dynamic>;
227 
228 // allocator:
229 
230  form_lazy_unop (const Unop& unop, const Expr& expr)
231  : base(),
232  _unop(unop),
233  _expr(expr)
234  {}
235 
236 // accessors:
237 
238  const geo_type& get_geo() const { return _expr.get_geo(); }
239  const space_type& get_trial_space() const { return _expr.get_trial_space(); }
240  const space_type& get_test_space() const { return _expr.get_test_space(); }
241  bool is_on_band() const { return false; }
242  band_type get_band() const { return band_type(); }
243 
244  void initialize (const geo_type& omega_K) { _expr.initialize (omega_K); }
245  bool is_diagonal() const { return _expr.is_diagonal(); }
246 
247  void evaluate (
248  const geo_type& omega_K,
249  const geo_element& K,
250  matrix_element_type& ak) const;
251 // data:
252 protected:
253  Unop _unop;
254  Expr _expr;
255 };
256 // concept;
257 template<class Unop, class Expr>
258 struct is_form_lazy <form_lazy_unop<Unop,Expr> > : std::true_type {};
259 
260 // inlined:
261 template<class Unop, class Expr>
262 void
264  const geo_type& omega_K,
265  const geo_element& K,
266  matrix_element_type& bk) const
267 {
269  _expr.evaluate (omega_K, K, ak);
270  bk = ak.unaryExpr(_unop);
271 }
272 
273 }// namespace details
274 
276 #define _RHEOLEF_form_lazy_unop(OP,NAME) \
277 template<class Expr, class Sfinae = typename std::enable_if<details::is_form_lazy<Expr>::value, Expr>::type> \
278 details::form_lazy_unop<NAME,Expr> \
279 operator OP (const Expr& a) \
280 { \
281  return details::form_lazy_unop<NAME,Expr> (NAME(),a); \
282 }
285 #undef _RHEOLEF_form_lazy_unop
286 // -------------------------------------------------------------------
287 // 3.2. invert
288 // -------------------------------------------------------------------
289 namespace details {
290 
291 template<class Expr>
293 public :
294 // definitions:
295 
297  using memory_type = typename Expr::memory_type;
298  using scalar_type = typename Expr::scalar_type;
299  using space_type = typename Expr::space_type;
300  using geo_type = typename Expr::geo_type;
303  using matrix_element_type = Eigen::Matrix<scalar_type,Eigen::Dynamic,Eigen::Dynamic>;
304 
305 // allocator:
306 
307  form_lazy_invert_rep (const Expr& expr)
308  : _expr(expr),
309  _prev_omega_K(),
310  _prev_K_dis_ie(std::numeric_limits<size_type>::max()),
311  _prev_ak()
312  {}
313 
314 // accessors:
315 
316  const geo_type& get_geo() const { return _expr.get_geo(); }
317  const space_type& get_trial_space() const { return _expr.get_trial_space(); }
318  const space_type& get_test_space() const { return _expr.get_test_space(); }
319  bool is_on_band() const { return false; }
320  band_type get_band() const { return band_type(); }
321 
322  void initialize (const geo_type& omega_K) const;
323  bool is_diagonal() const { return _expr.is_diagonal(); }
324 
325  void evaluate (
326  const geo_type& omega_K,
327  const geo_element& K,
328  matrix_element_type& ak) const;
329 // data:
330 protected:
331  mutable Expr _expr;
335 };
336 // inlined:
337 template<class Expr>
338 void
340 {
341  check_macro (get_trial_space() == get_test_space(),
342  "inv(form): spaces "
343  << "\"" <<get_trial_space().name()<<"\" and \"" <<get_test_space().name()<<"\""
344  << " should be equal");
345 
346  check_macro (get_trial_space().get_constitution().have_compact_support_inside_element() &&
347  get_test_space().get_constitution().have_compact_support_inside_element(),
348  "inv(form): requires compact support inside elements (e.g. discontinuous or bubble)");
349  _expr.initialize (omega_K);
350  _prev_omega_K = omega_K;
351  _prev_K_dis_ie = std::numeric_limits<size_type>::max();
352 }
353 template<class Expr>
354 void
356  const geo_type& omega_K,
357  const geo_element& K,
358  matrix_element_type& ak) const
359 {
360  if (_prev_omega_K == omega_K && _prev_K_dis_ie == K.dis_ie()) {
361  trace_macro("inv(K="<<K.name()<<K.dis_ie()<<",prev="<<_prev_K_dis_ie<<"): re-use");
362  ak = _prev_ak;
363  return;
364  }
365  trace_macro("inv(K="<<K.name()<<K.dis_ie()<<",prev="<<_prev_K_dis_ie<<"): compute");
366  _expr.evaluate (omega_K, K, ak);
367 #ifdef _RHEOLEF_PARANO
368  check_macro (ak.rows() == ak.cols(), "inv: matrix should be square");
369 #endif // _RHEOLEF_PARANO
370  details::local_invert (ak, _expr.is_diagonal());
371  _prev_ak = ak; // expensive to compute, so memorize it for common subexpressions
372  _prev_omega_K = omega_K;
373  _prev_K_dis_ie = K.dis_ie();
374 }
375 
376 template<class Expr>
377 class form_lazy_invert: public smart_pointer_nocopy<form_lazy_invert_rep<Expr>>,
378  public form_lazy_base <form_lazy_invert<Expr>> {
379 public :
380 // definitions:
381 
385  using size_type = typename rep::size_type;
386  using memory_type = typename rep::memory_type;
387  using scalar_type = typename rep::scalar_type;
388  using space_type = typename rep::space_type;
389  using geo_type = typename rep::geo_type;
390  using band_type = typename rep::band_type;
392 
393 // allocator:
394 
395  form_lazy_invert (const Expr& expr)
396  : base1(new_macro(rep(expr))),
397  base2()
398  {}
399 
400 // accessors:
401 
402  const geo_type& get_geo() const { return base1::data().get_geo(); }
403  const space_type& get_trial_space() const { return base1::data().get_trial_space(); }
404  const space_type& get_test_space() const { return base1::data().get_test_space(); }
405  bool is_on_band() const { return base1::data().is_on_band(); }
406  band_type get_band() const { return base1::data().get_band(); }
407 
408  void initialize (const geo_type& omega_K) const { return base1::data().initialize (omega_K); }
409  bool is_diagonal() const { return base1::data().is_diagonal(); }
410 
411  void evaluate (
412  const geo_type& omega_K,
413  const geo_element& K,
414  matrix_element_type& ak) const
415  { return base1::data().evaluate (omega_K, K, ak); }
416 };
417 // concept;
418 template<class Expr> struct is_form_lazy <form_lazy_invert<Expr> > : std::true_type {};
419 
420 }// namespace details
421 
423 template<class Expr, class Sfinae = typename std::enable_if<details::is_form_lazy<Expr>::value, Expr>::type>
425 inv (const Expr& a)
426 {
428 }
429 // -------------------------------------------------------------------
430 // 3.3. transpose
431 // -------------------------------------------------------------------
432 namespace details {
433 
434 template<class Expr>
435 class form_lazy_transpose: public form_lazy_base<form_lazy_transpose<Expr>> {
436 public :
437 // definitions:
438 
441  using memory_type = typename Expr::memory_type;
442  using scalar_type = typename Expr::scalar_type;
443  using space_type = typename Expr::space_type;
444  using geo_type = typename Expr::geo_type;
447  using matrix_element_type = Eigen::Matrix<scalar_type,Eigen::Dynamic,Eigen::Dynamic>;
448 
449 // allocator:
450 
451  form_lazy_transpose (const Expr& expr)
452  : base(),
453  _expr(expr)
454  {}
455 
456 // accessors:
457 
458  const geo_type& get_geo() const { return _expr.get_geo(); }
459  const space_type& get_trial_space() const { return _expr.get_test_space(); } // swapped!
460  const space_type& get_test_space() const { return _expr.get_trial_space(); } // swapped!
461  bool is_on_band() const { return false; }
462  band_type get_band() const { return band_type(); }
463 
464  void initialize (const geo_type& omega_K) { _expr.initialize (omega_K); }
465  bool is_diagonal() const { return _expr.is_diagonal(); }
466 
467  void evaluate (
468  const geo_type& omega_K,
469  const geo_element& K,
470  matrix_element_type& ak) const;
471 // data:
472 protected:
473  Expr _expr;
474 };
475 // concept;
476 template<class Expr> struct is_form_lazy <form_lazy_transpose<Expr> > : std::true_type {};
477 
478 // inlined:
479 template<class Expr>
480 void
482  const geo_type& omega_K,
483  const geo_element& K,
484  matrix_element_type& ak) const
485 {
486  _expr.evaluate (omega_K, K, ak);
487  ak.transposeInPlace();
488 }
489 
490 }// namespace details
491 
493 template<class Expr, class Sfinae = typename std::enable_if<details::is_form_lazy<Expr>::value, Expr>::type>
495 trans (const Expr& a)
496 {
498 }
499 // -------------------------------------------------------------------
500 // 4. binary expressions
501 // -------------------------------------------------------------------
502 // 4.1. add & minus
503 // -------------------------------------------------------------------
504 namespace details {
505 
506 template<class Binop, class Expr1, class Expr2>
507 class form_lazy_add: public form_lazy_base<form_lazy_add<Binop,Expr1,Expr2>> {
508 public :
509 // definitions:
510 
513  using memory_type = typename Expr1::memory_type;
514  using scalar_type = typename Expr1::scalar_type;
515  using space_type = typename Expr1::space_type;
516  using geo_type = typename Expr1::geo_type;
519  using matrix_element_type = Eigen::Matrix<scalar_type,Eigen::Dynamic,Eigen::Dynamic>;
520 
521 // allocator:
522 
523  form_lazy_add (const Expr1& expr1, const Expr2& expr2)
524  : base(),
525  _binop(),
526  _expr1(expr1),
527  _expr2(expr2) {}
528 
529 // accessors:
530 
531  const geo_type& get_geo() const;
532  const space_type& get_trial_space() const { return _expr1.get_trial_space(); }
533  const space_type& get_test_space() const { return _expr1.get_test_space(); }
534  bool is_on_band() const { return false; }
535  band_type get_band() const { return band_type(); }
536 
537  void initialize (const geo_type& omega_K);
538  bool is_diagonal() const { return _expr1.is_diagonal() &&
539  _expr2.is_diagonal(); }
540 
541  void evaluate (
542  const geo_type& omega_K,
543  const geo_element& K,
544  matrix_element_type& ak) const;
545 // data:
546 protected:
547  Binop _binop;
549  Expr2 _expr2;
550 };
551 // concept;
552 template<class Binop, class Expr1, class Expr2>
553 struct is_form_lazy <form_lazy_add<Binop,Expr1,Expr2> > : std::true_type {};
554 
555 // inlined:
556 template<class Binop, class Expr1, class Expr2>
557 void
559 {
560  // TODO: subdomains e.g. robin
561  check_macro (_expr1.get_geo() == _expr2. get_geo(),
562  "lazy_add: different domain not yet supported");
563 
564  check_macro (_expr1.get_trial_space() == _expr2.get_trial_space() &&
565  _expr1. get_test_space() == _expr2. get_test_space(),
566  "lazy_add: incompatible spaces "
567  << "[\"" <<_expr1.get_trial_space().name()<<"\", \"" <<_expr1.get_test_space().name()<<"\"] and "
568  << "[\"" <<_expr2.get_trial_space().name()<<"\", \"" <<_expr2.get_test_space().name()<<"\"]");
569 
570  _expr1.initialize (omega_K);
571  _expr2.initialize (omega_K);
572 }
573 template<class Binop, class Expr1, class Expr2>
576 {
577  // TODO: subdomains e.g. robin
578  return _expr1.get_geo();
579 }
580 template<class Binop, class Expr1, class Expr2>
581 void
583  const geo_type& omega_K,
584  const geo_element& K,
585  matrix_element_type& ck) const
586 {
587  matrix_element_type ak, bk;
588  _expr1.evaluate (omega_K, K, ak);
589  _expr2.evaluate (omega_K, K, bk);
590 #ifdef _RHEOLEF_PARANO
591  check_macro (ak.rows() == bk.rows() && ak.cols() == bk.cols(), "a+b: invalid sizes");
592 #endif // _RHEOLEF_PARANO
593  ck = ak.binaryExpr (bk, _binop);
594 }
595 
596 }// namespace details
597 
599 #define _RHEOLEF_form_lazy_add(OP,NAME) \
600 template<class Expr1, class Expr2, \
601  class Sfinae1 = typename std::enable_if<details::is_form_lazy<Expr1>::value, Expr1>::type, \
602  class Sfinae2 = typename std::enable_if<details::is_form_lazy<Expr2>::value, Expr2>::type> \
603 details::form_lazy_add<details::NAME,Expr1,Expr2> \
604 operator OP (const Expr1& a, const Expr2& b) \
605 { \
606  return details::form_lazy_add<details::NAME,Expr1,Expr2> (a,b); \
607 }
609 _RHEOLEF_form_lazy_add(-,minus)
610 #undef _RHEOLEF_form_lazy_add
611 // -------------------------------------------------------------------
612 // 4.2. multiply
613 // -------------------------------------------------------------------
614 // TODO: shared on common subexpressions
615 namespace details {
616 
617 template<class Expr1, class Expr2>
619 public :
620 // definitions:
621 
623  using memory_type = typename Expr1::memory_type;
624  using scalar_type = typename Expr1::scalar_type;
625  using space_type = typename Expr1::space_type;
626  using geo_type = typename Expr1::geo_type;
629  using matrix_element_type = Eigen::Matrix<scalar_type,Eigen::Dynamic,Eigen::Dynamic>;
630 
631 // allocator:
632 
633  form_lazy_multiply_rep (const Expr1& expr1, const Expr2& expr2)
634  : _expr1(expr1),
635  _expr2(expr2),
636  _prev_omega_K(),
637  _prev_K_dis_ie(std::numeric_limits<size_type>::max()),
638  _prev_ck()
639  {}
640 
641 // accessors:
642 
643  const geo_type& get_geo() const;
644  const space_type& get_trial_space() const { return _expr2.get_trial_space(); }
645  const space_type& get_test_space() const { return _expr1.get_test_space(); }
646  bool is_on_band() const { return false; }
647  band_type get_band() const { return band_type(); }
648 
649  void initialize (const geo_type& omega_K) const;
650  bool is_diagonal() const { return _expr1.is_diagonal() &&
651  _expr2.is_diagonal(); }
652 
653  void evaluate (
654  const geo_type& omega_K,
655  const geo_element& K,
656  matrix_element_type& ak) const;
657 // data:
658 protected:
659  mutable Expr1 _expr1;
660  mutable Expr2 _expr2;
664 };
665 // inlined:
666 template<class Expr1, class Expr2>
667 void
669 {
670  // TODO: subdomains e.g. robin
671  check_macro (_expr1.get_geo() == _expr2. get_geo(),
672  "lazy_multiply: different domain not yet supported");
673 
674  check_macro (_expr1.get_trial_space() == _expr2. get_test_space(),
675  "lazy_multiply: incompatible spaces \""
676  <<_expr1.get_trial_space().name()<<"\" and \""
677  <<_expr2. get_test_space().name()<<"\"");
678 
679  if (! _expr1. get_test_space().get_constitution().have_compact_support_inside_element() ||
680  ! _expr1.get_trial_space().get_constitution().have_compact_support_inside_element() ||
681  ! _expr2.get_trial_space().get_constitution().have_compact_support_inside_element() ) {
682  warning_macro("lazy_multiply: requires compact support inside elements (e.g. discontinuous or bubble)");
683  warning_macro("lazy_multiply: space was "
684  << "[\"" <<_expr1.get_trial_space().name()<<"\", \"" <<_expr1.get_test_space().name()<<"\"] and "
685  << "[\"" <<_expr2.get_trial_space().name()<<"\", \"" <<_expr2.get_test_space().name()<<"\"]");
686  fatal_macro("lazy_multiply: HINT: convert to \"form\" before to do the product");
687  }
688  _expr1.initialize (omega_K);
689  _expr2.initialize (omega_K);
690  _prev_omega_K = omega_K;
691  _prev_K_dis_ie = std::numeric_limits<size_type>::max();
692  trace_macro("mult(omega_K="<<omega_K.name()<<"): init prev:="<<_prev_K_dis_ie<<" for this="<<(this));
693 }
694 template<class Expr1, class Expr2>
697 {
698  // TODO: subdomains e.g. robin
699  return _expr1.get_geo();
700 }
701 template<class Expr1, class Expr2>
702 void
704  const geo_type& omega_K,
705  const geo_element& K,
706  matrix_element_type& ck) const
707 {
708  if (_prev_omega_K == omega_K && _prev_K_dis_ie == K.dis_ie()) {
709  trace_macro("mult(K="<<K.name()<<K.dis_ie()<<",prev="<<_prev_K_dis_ie<<"): re-use this="<<(this));
710  ck = _prev_ck;
711  return;
712  }
713  trace_macro("mult(K="<<K.name()<<K.dis_ie()<<",prev="<<_prev_K_dis_ie<<"): compute this="<<(this));
714  matrix_element_type ak, bk;
715  _expr1.evaluate (omega_K, K, ak);
716  _expr2.evaluate (omega_K, K, bk);
717 #ifdef _RHEOLEF_PARANO
718  check_macro (ak.cols() == bk.rows(), "a*b: invalid sizes");
719 #endif // _RHEOLEF_PARANO
720  ck = ak*bk;
721  _prev_ck = ck; // expensive to compute, so memorize it for common subexpressions
722  _prev_omega_K = omega_K;
723  _prev_K_dis_ie = K.dis_ie();
724  trace_macro("mult(K="<<K.name()<<K.dis_ie()<<"): prev:="<<_prev_K_dis_ie<<" for this="<<(this));
725 }
726 template<class Expr1, class Expr2>
727 class form_lazy_multiply: public smart_pointer_nocopy<form_lazy_multiply_rep<Expr1,Expr2>>,
728  public form_lazy_base <form_lazy_multiply<Expr1,Expr2>> {
729 public :
730 // definitions:
731 
735  using size_type = typename rep::size_type;
736  using memory_type = typename rep::memory_type;
737  using scalar_type = typename rep::scalar_type;
738  using space_type = typename rep::space_type;
739  using geo_type = typename rep::geo_type;
740  using band_type = typename rep::band_type;
742 
743 // allocator:
744 
745  form_lazy_multiply (const Expr1& expr1, const Expr2& expr2)
746  : base1(new_macro(rep(expr1,expr2))),
747  base2()
748  {}
749 
750 // accessors:
751 
752  const geo_type& get_geo() const { return base1::data().get_geo(); }
753  const space_type& get_trial_space() const { return base1::data().get_trial_space(); }
754  const space_type& get_test_space() const { return base1::data().get_test_space(); }
755  bool is_on_band() const { return base1::data().is_on_band(); }
756  band_type get_band() const { return base1::data().get_band(); }
757 
758  void initialize (const geo_type& omega_K) const { base1::data().initialize (omega_K); }
759  bool is_diagonal() const { return base1::data().is_diagonal(); }
760 
761  void evaluate (
762  const geo_type& omega_K,
763  const geo_element& K,
764  matrix_element_type& ak) const
765  { base1::data().evaluate (omega_K, K, ak); }
766 
767 };
768 // concept;
769 template<class Expr1, class Expr2> struct is_form_lazy <form_lazy_multiply<Expr1,Expr2> > : std::true_type {};
770 
771 }// namespace details
772 
774 template<class Expr1, class Expr2,
775  class Sfinae1 = typename std::enable_if<details::is_form_lazy<Expr1>::value, Expr1>::type,
776  class Sfinae2 = typename std::enable_if<details::is_form_lazy<Expr2>::value, Expr2>::type>
778 operator* (const Expr1& a, const Expr2& b)
779 {
781 }
782 // -------------------------------------------------------
783 // 5. form allocator
784 // -------------------------------------------------------
785 template<class T, class M>
786 template<class Expr, class Sfinae>
787 inline
788 form_basic<T,M>&
790 {
791  // here is the main call to the effective computation
792  // of all sparse matrix asssembly from element matrix:
793  convert_from_form_lazy (a);
794  return *this;
795 }
796 template<class T, class M>
797 template<class Expr, class Sfinae>
798 inline
800 : _X(), _Y(), _uu(), _ub(), _bu(), _bb()
801 {
802  operator= (a);
803 }
804 
805 }// namespace rheolef
806 # endif /* _RHEOLEF_FORM_LAZY_EXPR_H */
const geo_type & get_geo() const
band_basic< float_type, memory_type > band_type
geo_element::size_type size_type
typename Expr1::scalar_type scalar_type
void evaluate(const geo_type &omega_K, const geo_element &K, matrix_element_type &ak) const
const space_type & get_trial_space() const
typename Expr1::memory_type memory_type
form_lazy_add(const Expr1 &expr1, const Expr2 &expr2)
void initialize(const geo_type &omega_K)
Eigen::Matrix< scalar_type, Eigen::Dynamic, Eigen::Dynamic > matrix_element_type
const space_type & get_test_space() const
typename float_traits< scalar_type >::type float_type
typename Expr1::space_type space_type
typename Expr1::geo_type geo_type
band_basic< float_type, memory_type > band_type
void evaluate(const geo_type &omega_K, const geo_element &K, matrix_element_type &ak) const
const space_type & get_trial_space() const
void initialize(const geo_type &omega_K) const
typename Expr::scalar_type scalar_type
typename Expr::memory_type memory_type
Eigen::Matrix< scalar_type, Eigen::Dynamic, Eigen::Dynamic > matrix_element_type
const space_type & get_test_space() const
typename float_traits< scalar_type >::type float_type
typename Expr::space_type space_type
void evaluate(const geo_type &omega_K, const geo_element &K, matrix_element_type &ak) const
const space_type & get_trial_space() const
typename rep::memory_type memory_type
void initialize(const geo_type &omega_K) const
const geo_type & get_geo() const
typename rep::size_type size_type
typename rep::band_type band_type
typename rep::space_type space_type
typename rep::geo_type geo_type
const space_type & get_test_space() const
typename rep::scalar_type scalar_type
typename rep::matrix_element_type matrix_element_type
band_basic< float_type, memory_type > band_type
typename Expr1::scalar_type scalar_type
void evaluate(const geo_type &omega_K, const geo_element &K, matrix_element_type &ak) const
const space_type & get_trial_space() const
typename Expr1::memory_type memory_type
void initialize(const geo_type &omega_K) const
form_lazy_multiply_rep(const Expr1 &expr1, const Expr2 &expr2)
Eigen::Matrix< scalar_type, Eigen::Dynamic, Eigen::Dynamic > matrix_element_type
const space_type & get_test_space() const
typename float_traits< scalar_type >::type float_type
typename Expr1::space_type space_type
void evaluate(const geo_type &omega_K, const geo_element &K, matrix_element_type &ak) const
const space_type & get_trial_space() const
typename rep::memory_type memory_type
void initialize(const geo_type &omega_K) const
const geo_type & get_geo() const
typename rep::size_type size_type
typename rep::band_type band_type
typename rep::space_type space_type
const space_type & get_test_space() const
typename rep::scalar_type scalar_type
form_lazy_multiply(const Expr1 &expr1, const Expr2 &expr2)
typename rep::matrix_element_type matrix_element_type
band_basic< float_type, memory_type > band_type
void evaluate(const geo_type &omega_K, const geo_element &K, matrix_element_type &ak) const
const space_type & get_trial_space() const
const geo_type & get_geo() const
void initialize(const geo_type &omega_K)
typename Expr::scalar_type scalar_type
typename Expr::memory_type memory_type
Eigen::Matrix< scalar_type, Eigen::Dynamic, Eigen::Dynamic > matrix_element_type
const space_type & get_test_space() const
typename float_traits< scalar_type >::type float_type
typename Expr::space_type space_type
band_basic< float_type, memory_type > band_type
geo_element::size_type size_type
void evaluate(const geo_type &omega_K, const geo_element &K, matrix_element_type &ak) const
const space_type & get_trial_space() const
const geo_type & get_geo() const
void initialize(const geo_type &omega_K)
typename Expr::geo_type geo_type
typename Expr::scalar_type scalar_type
typename Expr::memory_type memory_type
Eigen::Matrix< scalar_type, Eigen::Dynamic, Eigen::Dynamic > matrix_element_type
form_lazy_unop(const Unop &unop, const Expr &expr)
const space_type & get_test_space() const
typename float_traits< scalar_type >::type float_type
typename Expr::space_type space_type
form_basic< T, M > & operator=(const form_basic< T, M > &)
Definition: form.h:329
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
size_type dis_ie() const
Definition: geo_element.h:163
char name() const
Definition: geo_element.h:169
rheolef::std type
rheolef::std Expr1
dot(x,y): see the expression page for the full documentation
#define trace_macro(message)
Definition: dis_macros.h:111
#define fatal_macro(message)
Definition: dis_macros.h:33
#define warning_macro(message)
Definition: dis_macros.h:53
void get_geo(istream &in, my_geo &omega)
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
void local_invert(Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &m, bool is_diag)
This file is part of Rheolef.
tensor_basic< T > inv(const tensor_basic< T > &a, size_t d)
Definition: tensor.cc:219
csr< T, sequential > trans(const csr< T, sequential > &a)
trans(a): see the form page for the full documentation
Definition: csr.h:455
csr< T, sequential > operator*(const T &lambda, const csr< T, sequential > &a)
Definition: csr.h:437
_RHEOLEF_form_lazy_unop(+, details::unary_plus) _RHEOLEF_form_lazy_unop(-
_RHEOLEF_form_lazy_add(+, plus) _RHEOLEF_form_lazy_add(-