Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
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
27namespace 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
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
204namespace rheolef {
205
206// -------------------------------------------------------------------
207// 3. unary expressions
208// -------------------------------------------------------------------
209// 3.1. unary_minus & unary_plus
210// -------------------------------------------------------------------
211namespace details {
212
213template<class Unop, class Expr>
214class form_lazy_unop: public form_lazy_base<form_lazy_unop<Unop,Expr>> {
215public :
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:
252protected:
253 Unop _unop;
254 Expr _expr;
255};
256// concept;
257template<class Unop, class Expr>
258struct is_form_lazy <form_lazy_unop<Unop,Expr> > : std::true_type {};
259
260// inlined:
261template<class Unop, class Expr>
262void
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) \
277template<class Expr, class Sfinae = typename std::enable_if<details::is_form_lazy<Expr>::value, Expr>::type> \
278details::form_lazy_unop<NAME,Expr> \
279operator 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// -------------------------------------------------------------------
289namespace details {
290
291template<class Expr>
293public :
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),
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:
330protected:
331 mutable Expr _expr;
335};
336// inlined:
337template<class Expr>
338void
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}
353template<class Expr>
354void
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
376template<class Expr>
377class form_lazy_invert: public smart_pointer_nocopy<form_lazy_invert_rep<Expr>>,
378 public form_lazy_base <form_lazy_invert<Expr>> {
379public :
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;
391 using matrix_element_type = typename rep::matrix_element_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;
418template<class Expr> struct is_form_lazy <form_lazy_invert<Expr> > : std::true_type {};
419
420}// namespace details
421
423template<class Expr, class Sfinae = typename std::enable_if<details::is_form_lazy<Expr>::value, Expr>::type>
425inv (const Expr& a)
426{
428}
429// -------------------------------------------------------------------
430// 3.3. transpose
431// -------------------------------------------------------------------
432namespace details {
433
434template<class Expr>
435class form_lazy_transpose: public form_lazy_base<form_lazy_transpose<Expr>> {
436public :
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:
472protected:
473 Expr _expr;
474};
475// concept;
476template<class Expr> struct is_form_lazy <form_lazy_transpose<Expr> > : std::true_type {};
477
478// inlined:
479template<class Expr>
480void
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
493template<class Expr, class Sfinae = typename std::enable_if<details::is_form_lazy<Expr>::value, Expr>::type>
495trans (const Expr& a)
496{
498}
499// -------------------------------------------------------------------
500// 4. binary expressions
501// -------------------------------------------------------------------
502// 4.1. add & minus
503// -------------------------------------------------------------------
504namespace details {
505
506template<class Binop, class Expr1, class Expr2>
507class form_lazy_add: public form_lazy_base<form_lazy_add<Binop,Expr1,Expr2>> {
508public :
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:
546protected:
547 Binop _binop;
548 Expr1 _expr1;
549 Expr2 _expr2;
550};
551// concept;
552template<class Binop, class Expr1, class Expr2>
553struct is_form_lazy <form_lazy_add<Binop,Expr1,Expr2> > : std::true_type {};
554
555// inlined:
556template<class Binop, class Expr1, class Expr2>
557void
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}
573template<class Binop, class Expr1, class Expr2>
576{
577 // TODO: subdomains e.g. robin
578 return _expr1.get_geo();
579}
580template<class Binop, class Expr1, class Expr2>
581void
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) \
600template<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> \
603details::form_lazy_add<details::NAME,Expr1,Expr2> \
604operator OP (const Expr1& a, const Expr2& b) \
605{ \
606 return details::form_lazy_add<details::NAME,Expr1,Expr2> (a,b); \
607}
610#undef _RHEOLEF_form_lazy_add
611// -------------------------------------------------------------------
612// 4.2. multiply
613// -------------------------------------------------------------------
614// TODO: shared on common subexpressions
615namespace details {
616
617template<class Expr1, class Expr2>
619public :
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),
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:
658protected:
659 mutable Expr1 _expr1;
660 mutable Expr2 _expr2;
664};
665// inlined:
666template<class Expr1, class Expr2>
667void
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}
694template<class Expr1, class Expr2>
697{
698 // TODO: subdomains e.g. robin
699 return _expr1.get_geo();
700}
701template<class Expr1, class Expr2>
702void
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}
726template<class Expr1, class Expr2>
727class form_lazy_multiply: public smart_pointer_nocopy<form_lazy_multiply_rep<Expr1,Expr2>>,
728 public form_lazy_base <form_lazy_multiply<Expr1,Expr2>> {
729public :
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;
741 using matrix_element_type = typename rep::matrix_element_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;
769template<class Expr1, class Expr2> struct is_form_lazy <form_lazy_multiply<Expr1,Expr2> > : std::true_type {};
770
771}// namespace details
772
774template<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>
778operator* (const Expr1& a, const Expr2& b)
779{
781}
782// -------------------------------------------------------
783// 5. form allocator
784// -------------------------------------------------------
785template<class T, class M>
786template<class Expr, class Sfinae>
787inline
788form_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}
796template<class T, class M>
797template<class Expr, class Sfinae>
798inline
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
typename Expr1::memory_type memory_type
const space_type & get_test_space() const
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
typename float_traits< scalar_type >::type float_type
typename Expr1::space_type space_type
const space_type & get_trial_space() const
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_test_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
typename float_traits< scalar_type >::type float_type
const space_type & get_trial_space() const
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_test_space() const
typename rep::memory_type memory_type
void initialize(const geo_type &omega_K) const
typename rep::size_type size_type
typename rep::band_type band_type
typename rep::space_type space_type
const space_type & get_trial_space() const
const geo_type & get_geo() const
smart_pointer_nocopy< rep > base1
typename rep::scalar_type scalar_type
form_lazy_invert_rep< Expr > rep
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
typename Expr1::memory_type memory_type
const space_type & get_test_space() const
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
typename float_traits< scalar_type >::type float_type
typename Expr1::space_type space_type
const space_type & get_trial_space() const
void evaluate(const geo_type &omega_K, const geo_element &K, matrix_element_type &ak) const
const space_type & get_test_space() const
typename rep::memory_type memory_type
void initialize(const geo_type &omega_K) const
typename rep::space_type space_type
form_lazy_multiply_rep< Expr1, Expr2 > rep
const space_type & get_trial_space() const
smart_pointer_nocopy< rep > base1
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_test_space() 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
typename float_traits< scalar_type >::type float_type
const space_type & get_trial_space() const
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_test_space() 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)
typename float_traits< scalar_type >::type float_type
const space_type & get_trial_space() const
const geo_type & get_geo() const
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
reference_element::size_type size_type
size_type dis_ie() const
#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)")
#define _RHEOLEF_form_lazy_unop(OP, NAME)
-a, +a: see the form page for the full documentation
#define _RHEOLEF_form_lazy_add(OP, NAME)
a+b, a-b: see the form page for the full documentation
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
STL namespace.