Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
field_lazy_terminal.h
Go to the documentation of this file.
1# ifndef _RHEOLEF_FIELD_LAZY_TERMINAL_H
2# define _RHEOLEF_FIELD_LAZY_TERMINAL_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// field_lazy = un-assembled field
24// as returned by lazy_integrate, lazy_interpolate or encapsulated field_basic
25// AUTHOR: Pierre.Saramito@imag.fr
26// DATE: 6 april 1920
27
28// SUMMARY:
29// 1. concept & base class : see field_lazy.h
30// 2. terminal
31// 2.1. field
32// 2.2. integrate
33// 2.3. integrate on a band
34// 2.4. interpolate
35// see also "field_lazy_node.h"
36//
37#include "rheolef/field_lazy.h"
38
39namespace rheolef {
40
41#ifdef TO_CLEAN
42// -------------------------------------------------------------------
43// 1. concept & base class
44// -------------------------------------------------------------------
45namespace details {
46// Define a trait type for detecting field expression valid arguments
47// template<class T> struct is_field_lazy: std::false_type {};
48// => should be defined in field.h for compilation reason,
49// otherwise Sfinae is always failed in field.h
50
51// Define a base class
52template<class Derived>
53class field_lazy_base {
54public:
55
56// no common methods yet
57
58protected:
59 Derived& derived() { return *static_cast< Derived*>(this); }
60 const Derived& derived() const { return *static_cast<const Derived*>(this); }
61};
62
63} // namespace details
64#endif // TO_CLEAN
65// -------------------------------------------------------------------
66// 2. terminal
67// -------------------------------------------------------------------
68// 2.1. field
69// -------------------------------------------------------------------
70// field_lazy_terminal_field encapsulates the field_basic class
71// It then appears as an unassembled field, for combining with
72// others unassembled fields in expressions
73//
74namespace details {
75
76template<class T, class M>
77class field_lazy_terminal_field: public field_lazy_base<field_lazy_terminal_field<T,M>> {
78public :
79// definitions:
80
83 using memory_type = M;
84 using scalar_type = T;
90 using vector_element_type = Eigen::Matrix<scalar_type,Eigen::Dynamic,1>;
91
92// allocator:
93
94 explicit field_lazy_terminal_field (const field_type& uh) : base(), _uh(uh) {}
95
96// accessors:
97
98 const geo_type& get_geo() const { return _uh.get_geo(); }
99 const space_type& get_space() const { return _uh.get_space(); }
100 bool is_on_band() const { return false; }
101 band_type get_band() const { return band_type(); }
102
103 void initialize (const geo_type& omega_K);
104
105 void evaluate (
106 const geo_type& omega_K,
107 const geo_element& K,
108 vector_element_type& uk) const;
109// data:
110protected:
112};
113// concept;
114template<class T, class M>
115struct is_field_lazy <field_lazy_terminal_field<T,M> > : std::true_type {};
116
117// inlined;
118template<class T, class M>
119void
121{
122 _uh.dis_dof_update();
123}
124template<class T, class M>
125void
127 const geo_type& omega_K,
128 const geo_element& K,
129 vector_element_type& uk) const
130{
131 std::vector<size_type> dis_idof;
132 _uh.get_space().get_constitution().assembly_dis_idof (omega_K, K, dis_idof);
133 size_type loc_ndof = dis_idof.size();
134 uk.resize (loc_ndof);
135 for (size_type loc_jdof = 0; loc_jdof < loc_ndof; ++loc_jdof) {
136 size_type dis_jdof = dis_idof[loc_jdof];
137 uk [loc_jdof] = _uh.dis_dof (dis_jdof);
138 }
139}
140
141}// namespace details
142// -------------------------------------------------------------------
143// 2.2. integrate
144// -------------------------------------------------------------------
145// field_lazy_terminal_integrate is returned by the integrate(domain) function
146// It contains the domain of integration and integration options (quadrature)
147// The sumitted elements K during evaluation belongs to a superset omega_K of "domain"
148// => the boolean is_on_domain[K_ie] behaves as an indicator function for "domain"
149//
150namespace details {
151
152template<class Expr>
154public :
155// definitions:
156
158 using memory_type = typename Expr::memory_type;
159 using scalar_type = typename Expr::scalar_type;
160 using space_type = typename Expr::space_type;
161 using geo_type = typename Expr::geo_type;
164 using vector_element_type = Eigen::Matrix<scalar_type,Eigen::Dynamic,1>;
165
166// allocators:
167
169 const Expr& expr,
170 const integrate_option& iopt);
171
173 const geo_type& domain,
174 const Expr& expr,
175 const integrate_option& iopt);
176
178 std::string domname,
179 const Expr& expr,
180 const integrate_option& iopt);
181
182// accessors:
183
184 const geo_type& get_geo() const { return _domain; }
185 const space_type& get_space() const { return _expr.get_vf_space(); }
186 bool is_on_band() const { return false; }
187 band_type get_band() const { return band_type(); }
188
189 void initialize (const geo_type& omega_K) const;
190
191 void evaluate (
192 const geo_type& omega_K,
193 const geo_element& K,
194 vector_element_type& uk) const;
195// data:
196protected:
198 mutable Expr _expr;
204// internal:
205 const geo_type& get_default_geo () const;
206};
207// inlined;
208template<class Expr>
210 const geo_type& domain,
211 const Expr& expr,
212 const integrate_option& iopt)
213: _domain(domain),
214 _expr(expr),
215 _iopt(iopt),
216 _is_on_domain(),
217 _prev_omega_K(),
218 _prev_K_dis_ie(std::numeric_limits<size_type>::max()),
219 _prev_uk()
220{
221}
222// missing domain:
223template<class Expr>
225 const Expr& expr,
226 const integrate_option& iopt)
227: _domain(),
228 _expr(expr),
229 _iopt(iopt),
230 _is_on_domain(),
231 _prev_omega_K(),
232 _prev_K_dis_ie(std::numeric_limits<size_type>::max()),
233 _prev_uk()
234{
236}
237// missing domain:
238template<class Expr>
240 std::string domname,
241 const Expr& expr,
242 const integrate_option& iopt)
243: _domain(),
244 _expr(expr),
245 _iopt(iopt),
246 _is_on_domain(),
247 _prev_omega_K(),
248 _prev_K_dis_ie(std::numeric_limits<size_type>::max()),
249 _prev_uk()
250{
251 _domain = get_default_geo() [domname];
252}
253template<class Expr>
256{
257 // TODO: improve the automatic determination of the domain ?
258 return get_space().get_constitution().get_geo();
259}
260template<class Expr>
261void
263{
265 _iopt._is_on_interface = false;
266 _iopt._is_inside_on_local_sides = false;
267 _expr.initialize (_domain, _iopt);
268 _prev_omega_K = omega_K;
269 _prev_K_dis_ie = std::numeric_limits<size_type>::max();
270 // --------------------------------
271 // initialize the domain indicator:
272 // --------------------------------
273 // TODO: sub-domain
274 // _domain could be a subdomain of omega_K:
275 // in that case evaluate(omega_K,K) should return a zero uk vector
276 // when K do not belongs to _domain
277 size_type K_map_d = omega_K.dimension();
278 _is_on_domain.resize (omega_K.geo_element_ownership(K_map_d));
279 std::fill (_is_on_domain.begin(), _is_on_domain.end(), false);
280 if (_domain.map_dimension() == omega_K.map_dimension()) {
281 if (_domain == omega_K || _domain.variant() != geo_abstract_base_rep<float_type>::geo_domain) {
282#ifdef TODO
283 // => _domain has the same element numbering as omega_K
284 // _domain.variant() is a geo or a geo_domain_indirect
285 trace_macro ("lazy_int(init): domain="<<_domain.name()<<": same numbering as "<<omega_K.name()<<"...");
286 size_type first_dis_ie = omega_K.geo_element_ownership (omega_K.map_dimension()).first_index();
287 trace_macro ("lazy_int(init): first_dis_ie="<<first_dis_ie);
288 trace_macro ("lazy_int(init): _is_on_domain.size="<<_is_on_domain.size());
289 for (size_type dom_ie = 0, dom_ne = _domain.size(); dom_ie < dom_ne; ++dom_ie) {
290 trace_macro ("lazy_int(init): dom_ie="<<dom_ie);
291 const geo_element& K = _domain [dom_ie];
292 size_type dis_ie = K.dis_ie();
293 trace_macro ("lazy_int(init): K.dis_ie="<<dis_ie);
294 check_macro (first_dis_ie <= dis_ie, "invalid index"); // PARANO
295 size_type ie = dis_ie - first_dis_ie;
296 trace_macro ("lazy_int(init): ie="<<ie);
297 check_macro (ie <= _is_on_domain.size(), "invalid index"); // PARANO
298 _is_on_domain [ie] = true;
299 }
300#endif // TODO
301 trace_macro ("lazy_int(init): domain="<<_domain.name()<<": same numbering as "<<omega_K.name()<<" done");
302 } else {
303 error_macro ("geo_domain="<<_domain.name()<<" subset of "<<omega_K.name()<<": not yet");
304 }
305 } else if (_domain.map_dimension()+1 == omega_K.map_dimension()) {
306 check_macro (_domain.get_background_geo() == omega_K, "unexpected domain \""
307 <<_domain.name()<<"\" as subset of \""<<omega_K.name()<<"\"");
308 if (_domain.variant() != geo_abstract_base_rep<float_type>::geo_domain) {
309 trace_macro ("lazy_int(init): domain(bdry or sides)="<<_domain.name()<<" subset of "<<omega_K.name()<<"...");
310 for (size_type dom_ie = 0, dom_ne = _domain.size(); dom_ie < dom_ne; ++dom_ie) {
311 const geo_element& S = _domain [dom_ie];
312 size_type dis_is = S.dis_ie();
313 // _domain="sides" or "boundary" or other d-1 sides domain
314 // for a side S, will access to its neighbours K0 & K1
315 size_type K_dis_ie = S.master(0); // TODO: domain = "interface" : orient could be -1 and then choose S.master(1) !!
316 check_macro (K_dis_ie != std::numeric_limits<size_type>::max(),
317 "unexpected isolated side S="<<S.name()<<S.dis_ie());
318 const geo_element& K = omega_K.dis_get_geo_element (K_map_d, K_dis_ie);
319 size_type dis_ie = K.dis_ie(); // not necessarily on the same proc as S
320 _is_on_domain.dis_entry (dis_ie) = true;
321 }
322 _is_on_domain.dis_entry_assembly();
323 trace_macro ("lazy_int(init): domain(bdry or sides)="<<_domain.name()<<" subset of "<<omega_K.name()<<" done");
324 } else {
325 error_macro ("geo_domain(bdry or sides)="<<_domain.name()<<" subset of "<<omega_K.name()<<": not yet");
326 }
327 } else {
328 error_macro ("unsupported domain \"" << _domain.name()
329 << "\" with map_dimension=" << _domain.map_dimension()
330 << " when integrated in geometry \"" << omega_K.name()
331 << "\" with map_dimension=" << omega_K.map_dimension());
332 }
333}
334/*
335 compute local idofs on S from the basis in K ?
336 all dofs in K are numberd from 0 to ndof(K)-1
337 for each S :
338 is_in_S[0:ndof(K)-1] = false
339 for each subgeo_dim=0 to S.dim :
340 for each subgeo T of S with T.dim=subgeo_dim
341 get T.first_idof and T.last_idof from basis infos
342 for idof=T.first_idof to T.last_idof-1
343 is_in_S[idof] = true
344 ndof_S = 0
345 for idof = 0 to ndof(K)-1
346 if is_in_S[idof] then ndof_S++
347 idof_S2idof.resize (ndof_S)
348 idof_S = 0
349 for idof = 0 to ndof(K)-1
350 if is_in_S[idof] then
351 idof_S2idof[idof_S++] = idof
352
353 finally:
354 for idof_S = 0 to ndof_S-1
355 uk[idof_S2idof[idof_S] = us[idof_S]
356 The computation of idof_S2idof could be done on the basis(hat_K)
357 one time for all for all sides S at initialization or on the fly
358 at the first request
359 - first step : we do it here, for all K
360 - second step : we move it on the basis class in an internal mutable data
361 For a space product eg Yh=Xh*Mh we have to concactenate the idof_S2idof arrays
362 => will be a member function of space class ?
363*/
364template <class T>
365void
367 const basis_basic<T>& b,
368 const reference_element& hat_K,
369 std::array<std::vector<size_t>,
371 std::array<size_t,
373 size_t& ndof_K)
374{
375 trace_macro ("compute_idof_S2idof_K: hat_K="<<hat_K.name()<<"...");
377 check_macro (hat_K.dimension() > 0, "idof_S2idof_K: invalid K.dimension=0");
378 size_type sid_dim = hat_K.dimension()-1;
379 trace_macro ("compute_idof_S2idof_K: sid_dim="<<sid_dim);
380 ndof_K = b.ndof (hat_K);
381 // ------------------------------
382 // 1) is_in_S [loc_isid] [idof_K]
383 // ------------------------------
384 std::array<std::vector<bool>, reference_element::max_side_by_variant> is_in_S;
385 std::array<std::map<size_type,size_type>, reference_element::max_side_by_variant> idof_S2idof_K_map;
386 for (size_type loc_isid = 0, loc_nsid = hat_K.n_subgeo(sid_dim); loc_isid < loc_nsid; ++loc_isid) {
387 is_in_S[loc_isid].resize (ndof_K);
388 std::fill (is_in_S[loc_isid].begin(), is_in_S[loc_isid].end(), false);
389 reference_element hat_S = hat_K.subgeo (sid_dim, loc_isid);
390 trace_macro ("compute_idof_S2idof_K: loc_isid="<<loc_isid<<", hat_S="<<hat_S.name());
391 size_type idof_S = 0;
392 for (size_type sub_sid_dim = 0; sub_sid_dim <= sid_dim; ++sub_sid_dim) {
393 size_type first_idof_K = b.first_idof_by_dimension (hat_K, sub_sid_dim);
394 trace_macro ("compute_idof_S2idof_K: first_idof_K(hat_K="<<hat_K.name()<<",sub_sid_dim="<<sub_sid_dim<<") = " <<first_idof_K);
395 for (size_type loc_isub_sid = 0; loc_isub_sid < hat_S.n_subgeo(sub_sid_dim); ++loc_isub_sid) {
396 reference_element hat_T = hat_S.subgeo (sub_sid_dim, loc_isub_sid);
397 size_type ndof_K_on_T = b.ndof_on_subgeo (hat_K.dimension(), hat_T.variant());
398 size_type loc_T_idx = 0;
399 switch (sub_sid_dim) {
400 case 0: {
401 loc_T_idx = hat_K.subgeo_local_vertex (sid_dim, loc_isid, loc_isub_sid);
402 break;
403 }
404 case 1: {
405 if (sid_dim == 1) {
406 loc_T_idx = loc_isid; // the side=edge in 2D itself
407 } else {
408 error_macro("compute_idof_S2idof_K: dof on edge in 3D: not yet");
409 }
410 break;
411 }
412 case 2:
413 default:{
414 if (sid_dim == 2) {
415 loc_T_idx = loc_isid; // the side=face in 3D itself
416 } else {
417 error_macro("compute_idof_S2idof_K: dof on face in 3D: not yet");
418 }
419 break;
420 }
421 }
422 trace_macro ("compute_idof_S2idof_K: ndof_K_on_T(hat_K.dim="<<hat_K.dimension()<<",hat_T="<<hat_T.name()<<") = " <<ndof_K_on_T);
423 trace_macro ("compute_idof_S2idof_K: T_idx="<<loc_T_idx);
424 size_type start_idof_K = first_idof_K + ndof_K_on_T*loc_T_idx;
425 size_type last_idof_K = first_idof_K + ndof_K_on_T*(loc_T_idx+1);
426 for (size_type idof_K = start_idof_K; idof_K < last_idof_K; ++idof_K) {
427 if (is_in_S [loc_isid] [idof_K]) continue;
428 is_in_S [loc_isid] [idof_K] = true;
429 idof_S2idof_K_map [loc_isid] [idof_S] = idof_K;
430 trace_macro ("idof_S2idof_K_map(isid="<<loc_isid<<",idof_S="<<idof_S<<") = "<<idof_K);
431 ++idof_S;
432 }
433 }
434 }
435 }
436#ifdef TO_CLEAN
437 // ------------------------------
438 // 2) count ndof on each S
439 // ------------------------------
440 for (size_type loc_isid = 0, loc_nsid = hat_K.n_subgeo(sid_dim); loc_isid < loc_nsid; ++loc_isid) {
441 ndof_S [loc_isid] = 0;
442 for (size_type idof_K = 0; idof_K < ndof_K; ++idof_K) {
443 if (is_in_S [loc_isid] [idof_K]) {
444 ndof_S [loc_isid]++;
445 }
446 }
447 trace_macro ("compute_idof_S2idof_K: ndof_S[isid="<<loc_isid<<"]="<<ndof_S[loc_isid]);
448 }
449#endif // TO_CLEAN
450 // ------------------------------
451 // 3) set idof_S2idof_K
452 // ------------------------------
453 for (size_type loc_isid = 0, loc_nsid = hat_K.n_subgeo(sid_dim); loc_isid < loc_nsid; ++loc_isid) {
454 ndof_S [loc_isid] = idof_S2idof_K_map [loc_isid].size();
455 idof_S2idof_K[loc_isid].resize (ndof_S [loc_isid]);
456 for (auto p : idof_S2idof_K_map [loc_isid]) {
457 size_type idof_S = p.first;
458 size_type idof_K = p.second;
459 if (is_in_S [loc_isid] [idof_K]) {
460 idof_S2idof_K [loc_isid] [idof_S] = idof_K;
461 trace_macro ("idof_K(isid="<<loc_isid<<",idof_S="<<idof_S<<") = "<< idof_K);
462 idof_S++;
463 }
464 }
465 }
466 trace_macro ("compute_idof_S2idof_K done");
467}
468template<class Expr>
469void
471 const geo_type& omega_K,
472 const geo_element& K,
473 vector_element_type& uk) const
474{
475 if (_prev_omega_K == omega_K && _prev_K_dis_ie == K.dis_ie()) {
476 uk = _prev_uk;
477 return;
478 }
479 if (_domain.map_dimension() == omega_K.map_dimension()) {
480 trace_macro ("lazy_int(eval): domain="<<_domain.name()<<": same numbering as "<<omega_K.name()<<"...");
481 _expr.evaluate (omega_K, K, uk);
482 trace_macro ("lazy_int(eval): domain="<<_domain.name()<<": same numbering as "<<omega_K.name()<<" done");
483 } else { // _domain.map_dimension()+1 == omega_K.map_dimension()
484 trace_macro ("lazy_int(eval): domain(bdry or sides)="<<_domain.name()<<" subset of "<<omega_K.name()<<"...");
485 std::array<std::vector<size_type>, reference_element::max_side_by_variant> idof_S2idof_K;
486 std::array<size_type, reference_element::max_side_by_variant> ndof_S;
487 size_type ndof_K = 0;
488 check_macro (get_space().size() == 0, "lazy_int(eval): "<<get_space().size()<<"-component space: not yet supported");
489 compute_idof_S2idof_K (get_space().get_basis(), K, idof_S2idof_K, ndof_S, ndof_K); // TODO: do it in space:: with concat
490 uk = vector_element_type::Zero (ndof_K, 1);
491 trace_macro ("lazy_int(eval): uk.size="<<uk.size());
492 check_macro (K.dimension() > 0, "lazy_int(eval): invalid K.dimension=0");
493 size_type sid_dim = K.dimension()-1;
495 for (size_type loc_isid = 0, nsid = K.n_subgeo(sid_dim); loc_isid < nsid; ++loc_isid) {
496 size_type dis_isid = K.subgeo_dis_index(sid_dim,loc_isid);
497 const geo_element& S = omega_K.dis_get_geo_element (sid_dim, dis_isid);
499 K.get_side_informations (S, sid);
500 _expr.evaluate (_domain, S, us);
501 trace_macro ("lazy_int(eval): us(loc_isid="<<loc_isid<<").size="<<us.size());
502 check_macro (size_type(us.size()) == ndof_S [loc_isid], "invalid us size");
503 for (size_type idof_S = 0; idof_S < ndof_S [loc_isid]; ++idof_S) {
504 size_type idof_K = idof_S2idof_K[loc_isid][idof_S];
505 trace_macro ("lazy_int(eval): uk[idof_K="<<idof_K<<"] += us[idof_S="<<idof_S<<"] = " << us[idof_S]);
506 uk[idof_K] += us[idof_S];
507 }
508 }
509 trace_macro ("lazy_int(eval): domain(bdry or sides)="<<_domain.name()<<" subset of "<<omega_K.name()<<" done");
510 }
511 _prev_uk = uk; // expensive to compute, so memorize it for common subexpressions
512 _prev_omega_K = omega_K;
513 _prev_K_dis_ie = K.dis_ie();
514 trace_macro("lazy_int(K="<<K.name()<<K.dis_ie()<<",prev="<<_prev_K_dis_ie<<"): compute");
515}
516
517template<class Expr>
518class field_lazy_terminal_integrate: public smart_pointer_nocopy<field_lazy_terminal_integrate_rep<Expr>>,
519 public field_lazy_base <field_lazy_terminal_integrate<Expr>> {
520public :
521// definitions:
522
526 using size_type = typename rep::size_type;
527 using memory_type = typename rep::memory_type;
528 using scalar_type = typename rep::scalar_type;
529 using space_type = typename rep::space_type;
530 using geo_type = typename rep::geo_type;
531 using band_type = typename rep::band_type;
532 using vector_element_type = typename rep::vector_element_type;
533
534// allocator:
535
537 const geo_type& domain,
538 const Expr& expr,
539 const integrate_option& iopt)
540 : base1(new_macro(rep(domain,expr,iopt))),
541 base2()
542 {}
543
545 const Expr& expr,
546 const integrate_option& iopt)
547 : base1(new_macro(rep(expr,iopt))),
548 base2()
549 {}
550
552 std::string domname,
553 const Expr& expr,
554 const integrate_option& iopt)
555 : base1(new_macro(rep(domname,expr,iopt))),
556 base2()
557 {}
558
559// accessors:
560
561 const geo_type& get_geo() const { return base1::data().get_geo(); }
562 const space_type& get_space() const { return base1::data().get_space(); }
563 bool is_on_band() const { return base1::data().is_on_band(); }
564 band_type get_band() const { return base1::data().get_band(); }
565
566 void initialize (const geo_type& omega_K) const { return base1::data().initialize (omega_K); }
567
568 void evaluate (
569 const geo_type& omega_K,
570 const geo_element& K,
571 vector_element_type& uk) const
572 { base1::data().evaluate (omega_K, K, uk); }
573};
574// concept;
575template<class Expr> struct is_field_lazy <field_lazy_terminal_integrate<Expr> > : std::true_type {};
576
577}// namespace details
578
579// 2.2.a) general call
580template <class Expr>
581inline
582typename
583std::enable_if<
585 ,details::field_lazy_terminal_integrate <Expr>
586>::type
589 const typename Expr::geo_type& domain,
590 const Expr& expr,
591 const integrate_option& iopt = integrate_option())
592{
593 return details::field_lazy_terminal_integrate<Expr> (domain, expr, iopt);
594}
595template <class Expr>
596inline
597typename
598std::enable_if<
599 details::is_field_expr_v2_variational_arg<Expr>::value
600 ,details::field_lazy_terminal_integrate <details::field_expr_quadrature_on_element<Expr>>
601>::type
605 const Expr& expr,
606 const integrate_option& iopt = integrate_option())
607{
609 return lazy_integrate (domain, expr_quad, iopt);
610}
611// 2.2.b) missing domain
612template <class Expr>
613inline
614typename
615std::enable_if<
616 details::is_field_expr_quadrature_arg<Expr>::value
617 ,details::field_lazy_terminal_integrate <Expr>
618>::type
621 const Expr& expr,
622 const integrate_option& iopt = integrate_option())
623{
625}
626template <class Expr>
627inline
628typename
629std::enable_if<
630 details::is_field_expr_v2_variational_arg<Expr>::value
631 ,details::field_lazy_terminal_integrate <details::field_expr_quadrature_on_element<Expr>>
632>::type
635 const Expr& expr,
636 const integrate_option& iopt = integrate_option())
637{
639 return lazy_integrate (expr_quad, iopt);
640}
641// 2.2.c) subdomain by its name
642template <class Expr>
643inline
644typename
645std::enable_if<
646 details::is_field_expr_quadrature_arg<Expr>::value
647 ,details::field_lazy_terminal_integrate <Expr>
648>::type
651 const std::string& domname,
652 const Expr& expr,
653 const integrate_option& iopt = integrate_option())
654{
655 return details::field_lazy_terminal_integrate<Expr> (domname, expr, iopt);
656}
657template <class Expr>
658inline
659typename
660std::enable_if<
661 details::is_field_expr_v2_variational_arg<Expr>::value
662 ,details::field_lazy_terminal_integrate <details::field_expr_quadrature_on_element<Expr>>
663>::type
666 const std::string& domname,
667 const Expr& expr,
668 const integrate_option& iopt = integrate_option())
669{
671 return lazy_integrate (domname, expr_quad, iopt);
672}
673// ----------------------------------------------
674// 2.3. integrate on a band
675// ----------------------------------------------
676namespace details {
677
678template<class Expr>
680public :
681// definitions:
682
684 using memory_type = typename Expr::memory_type;
685 using scalar_type = typename Expr::scalar_type;
686 using space_type = typename Expr::space_type;
687 using geo_type = typename Expr::geo_type;
690 using vector_element_type = Eigen::Matrix<scalar_type,Eigen::Dynamic,1>;
691
692// allocators:
693
695 const band_type& gh,
696 const Expr& expr,
697 const integrate_option& iopt);
698
699// accessors:
700
701 const geo_type& get_geo() const { return _gh.level_set(); }
702 const space_type& get_space() const { return _expr.get_vf_space(); }
703 bool is_on_band() const { return true; }
704 band_type get_band() const { return _gh; }
705
706 void initialize (const geo_type& omega_K) const;
707
708 void evaluate (
709 const geo_type& omega_K,
710 const geo_element& K,
711 vector_element_type& uk) const;
712// data:
713protected:
715 mutable Expr _expr;
720};
721// inlined;
722template<class Expr>
724 const band_type& gh,
725 const Expr& expr,
726 const integrate_option& iopt)
727: _gh(gh),
728 _expr(expr),
729 _iopt(iopt),
730 _prev_omega_K(),
731 _prev_K_dis_ie(std::numeric_limits<size_type>::max()),
732 _prev_uk()
733{
734}
735template<class Expr>
736void
738{
739 const geo_type& band = _gh.band(); // the 3D intersected tetras
740 const geo_type& bgd_omega = get_space().get_constitution().get_background_geo();
741 check_macro (omega_K == get_geo(),
742 "integrate on band: unexpected integration domain");
743 check_macro (band.get_background_geo() == bgd_omega,
744 "do_integrate: incompatible integration domain "<<_gh.level_set().name() << " and test function based domain "
745 << bgd_omega.name());
746 _expr.initialize (_gh, _iopt);
747 _prev_omega_K = omega_K;
748 _prev_K_dis_ie = std::numeric_limits<size_type>::max();
749}
750template<class Expr>
751void
753 const geo_type& omega_K,
754 const geo_element& K,
755 vector_element_type& uk) const
756{
757 if (_prev_omega_K == omega_K && _prev_K_dis_ie == K.dis_ie()) {
758 uk = _prev_uk;
759 return;
760 }
761 _expr.evaluate (omega_K, K, uk);
762 _prev_uk = uk; // expensive to compute, so memorize it for common subexpressions
763 _prev_omega_K = omega_K;
764 _prev_K_dis_ie = K.dis_ie();
765}
766template<class Expr>
767class field_lazy_terminal_integrate_band: public smart_pointer_nocopy<field_lazy_terminal_integrate_band_rep<Expr>>,
768 public field_lazy_base <field_lazy_terminal_integrate_band<Expr>> {
769public :
770// definitions:
771
775 using size_type = typename rep::size_type;
776 using memory_type = typename rep::memory_type;
777 using scalar_type = typename rep::scalar_type;
778 using space_type = typename rep::space_type;
779 using geo_type = typename rep::geo_type;
780 using band_type = typename rep::band_type;
781 using vector_element_type = typename rep::vector_element_type;
782
783// allocator:
784
786 const band_type& gh,
787 const Expr& expr,
788 const integrate_option& iopt)
789 : base1(new_macro(rep(gh,expr,iopt))),
790 base2()
791 {}
792
793// accessors:
794
795 const geo_type& get_geo() const { return base1::data().get_geo(); }
796 const space_type& get_space() const { return base1::data().get_space(); }
797 bool is_on_band() const { return base1::data().is_on_band(); }
798 band_type get_band() const { return base1::data().get_band(); }
799
800 void initialize (const geo_type& omega_K) const { return base1::data().initialize (omega_K); }
801
802 void evaluate (
803 const geo_type& omega_K,
804 const geo_element& K,
805 vector_element_type& uk) const
806 { base1::data().evaluate (omega_K, K, uk); }
807};
808// concept;
809template<class Expr> struct is_field_lazy <field_lazy_terminal_integrate_band<Expr> > : std::true_type {};
810
811}// namespace details
812
813template <class Expr>
814typename
815std::enable_if<
817 ,details::field_lazy_terminal_integrate_band <Expr>
818>::type
821 const band_basic<typename float_traits<typename Expr::scalar_type>::type,
822 typename Expr::memory_type>& gh,
823 const Expr& expr,
824 const integrate_option& iopt = integrate_option())
825{
827}
828template <class Expr>
829typename
830std::enable_if<
831 details::is_field_expr_v2_variational_arg<Expr>::value
832 ,details::field_lazy_terminal_integrate_band <details::field_expr_quadrature_on_element<Expr>>
833>::type
836 const band_basic<typename float_traits<typename Expr::scalar_type>::type,
837 typename Expr::memory_type>& gh,
838 const Expr& expr,
839 const integrate_option& iopt = integrate_option())
840{
841
843 return lazy_integrate (gh, expr_quad, iopt);
844}
845// ----------------------------------------------
846// 2.4. interpolate
847// ----------------------------------------------
848namespace details {
849
850template<class Expr>
852public :
853// definitions:
854
856 using memory_type = typename Expr::memory_type;
857 using value_type = typename Expr::value_type; // TODO: undeterminated_type<T>
858 using scalar_type = typename Expr::scalar_type;
860 using geo_type = geo_basic <float_type,memory_type>;
862 using band_type = band_basic <float_type,memory_type>;
863 using vector_element_type = Eigen::Matrix<scalar_type,Eigen::Dynamic,1>;
864
865// allocators:
866
868 const space_type& Xh,
869 const Expr& expr);
870
871// accessors:
872
873 const geo_type& get_geo() const { return _Xh.get_geo(); }
874 const space_type& get_space() const { return _Xh; }
875 bool is_on_band() const { return false; }
876 band_type get_band() const { return band_type(); }
877
878 void initialize (const geo_type& omega_K) const;
879
880 void evaluate (
881 const geo_type& omega_K,
882 const geo_element& K,
883 vector_element_type& uk) const;
884protected:
885// internals:
886
887 template <class Value>
888 typename std::enable_if<
890 ,void
891 >::type
892 evaluate_internal (
893 const geo_type& omega_K,
894 const geo_element& K,
895 vector_element_type& uk) const;
896
897 template <class Value>
898 typename std::enable_if<
900 ,void
901 >::type
903 const geo_type& omega_K,
904 const geo_element& K,
905 vector_element_type& uk) const;
906
907// data:
909 mutable Expr _expr;
915};
916// inlined;
917template<class Expr>
919 const space_type& Xh,
920 const Expr& expr)
921: _Xh(Xh),
922 _expr(expr),
923 _pops(),
924 _is_marked(),
925 _prev_omega_K(),
926 _prev_K_dis_ie(std::numeric_limits<size_type>::max()),
927 _prev_uk()
928{
929}
930template<class Expr>
931void
933{
934 check_macro (omega_K == _Xh.get_geo(),
935 "interpolate: incompatible space \""<<_Xh.name()<<"\" with geometry \""
936 << omega_K.name());
937 const basis_basic<float_type>& b = _Xh.get_basis();
938 const piola_fem<float_type>& pf = b.get_piola_fem();
939 integrate_option iopt;
940 _pops.initialize (omega_K.get_piola_basis(), b, iopt);
941 _expr.initialize (_Xh, _pops, iopt);
942 _expr.template valued_check<value_type>();
943 _prev_omega_K = omega_K;
944 _prev_K_dis_ie = std::numeric_limits<size_type>::max();
945 if (_Xh.get_basis().is_continuous()) {
946 std::set<size_type> ext_dis_idof;
947 _Xh.get_parent_subgeo_owner_dis_indexes (ext_dis_idof);
948 _is_marked.resize (_Xh.ownership());
949 std::fill (_is_marked.begin(), _is_marked.end(), false);
950 _is_marked.set_dis_indexes (ext_dis_idof);
951 } // if continuous
952}
953template<class Expr>
954template <class Value>
955typename std::enable_if<
957 ,void
958>::type
960 const geo_type& omega_K,
961 const geo_element& K,
962 vector_element_type& uk) const
963{
964 // 1) evaluate as a Value
965 Eigen::Matrix<Value,Eigen::Dynamic,1> value; // TODO: allocate it once on the mesh
966 _expr.evaluate (omega_K, K, value);
967 // 2) convert field "Values" at nodes of K to scalar "dofs"
968 const piola_fem<float_type>& pf = _Xh.get_basis().get_piola_fem();
969 if (pf.transform_need_piola()) {
970 const Eigen::Matrix<piola<float_type>,Eigen::Dynamic,1>& piola = _pops.get_piola (omega_K, K);
971 for (size_type loc_inod = 0, loc_nnod = value.size(); loc_inod < loc_nnod; ++loc_inod) {
972 // be carefull: piola_fem::inv_transform should support inplace call in the "value" arg
973 pf.inv_transform (piola[loc_inod], value[loc_inod], value[loc_inod]);
974 }
975 }
976 _Xh.get_basis().compute_dofs (K, value, uk);
977}
978template<class Expr>
979template <class Value>
980typename std::enable_if<
982 ,void
983>::type
985 const geo_type& omega_K,
986 const geo_element& K,
987 vector_element_type& uk) const
988{
989 using T = scalar_type;
990 // when valued_type is undeterminated at compile time e.g.
991 // field uh = lazy_interpolate (Xh, tau_h*vh);
992 // => undeterminated : could be scalar, vector or tensor valued
993 // so, ask to the Xh space at run-time:
994 switch (_Xh.valued_tag()) {
995#define _RHEOLEF_switch(VALUED,VALUE) \
996 case space_constant::VALUED: \
997 (*this).template evaluate_internal<VALUE> (omega_K, K, uk); break;
998_RHEOLEF_switch(scalar,T)
1001_RHEOLEF_switch(unsymmetric_tensor,tensor_basic<T>)
1002#ifdef TODO
1005#endif // TODO
1006#undef _RHEOLEF_switch
1007 default: error_macro ("unexpected argument valued tag="<<_Xh.valued_tag());
1008 }
1009}
1010template<class Expr>
1011void
1013 const geo_type& omega_K,
1014 const geo_element& K,
1015 vector_element_type& uk) const
1016{
1017 if (_prev_omega_K == omega_K && _prev_K_dis_ie == K.dis_ie()) {
1018 trace_macro("interpolate(K="<<K.name()<<K.dis_ie()<<",prev="<<_prev_K_dis_ie<<"): re-use");
1019 uk = _prev_uk;
1020 return;
1021 }
1022 trace_macro("interpolate(K="<<K.name()<<K.dis_ie()<<",prev="<<_prev_K_dis_ie<<"): compute");
1023 (*this).template evaluate_internal<value_type> (omega_K, K, uk);
1024
1025 // filter for continuous element: each dof value is transmitted only once
1026 // so the result could be summed as any others field_lazy<Expr>
1027 if (_Xh.get_basis().is_continuous()) {
1028 std::vector<size_type> dis_idx;
1029 _Xh.dis_idof (K, dis_idx);
1030 check_macro (dis_idx.size() == size_type(uk.size()),
1031 "invalid sizes: dis_idx.size="<<dis_idx.size()<<", uk.size="<<uk.size());
1032 size_type my_proc = _Xh.comm().rank();
1033 for (size_type loc_idof = 0, loc_ndof = dis_idx.size(); loc_idof < loc_ndof; ++loc_idof) {
1034 size_type dis_idof = dis_idx [loc_idof];
1035 size_type iproc = _Xh.get_parent_subgeo_owner (dis_idof);
1036 if (iproc != my_proc || _is_marked.dis_at (dis_idof)) {
1037 uk [loc_idof] = 0;
1038 continue;
1039 }
1040 // conserve the uk[loc_idof] value and remember it:
1041 _is_marked.dis_entry (dis_idof) = true;
1042 }
1043 }
1044 _prev_uk = uk; // expensive to compute, so memorize it for common subexpressions
1045 _prev_omega_K = omega_K;
1046 _prev_K_dis_ie = K.dis_ie();
1047}
1048template<class Expr>
1049class field_lazy_terminal_interpolate: public smart_pointer_nocopy<field_lazy_terminal_interpolate_rep<Expr>>,
1050 public field_lazy_base <field_lazy_terminal_interpolate<Expr>> {
1051public :
1052// definitions:
1053
1057 using size_type = typename rep::size_type;
1058 using memory_type = typename rep::memory_type;
1059 using scalar_type = typename rep::scalar_type;
1060 using space_type = typename rep::space_type;
1061 using geo_type = typename rep::geo_type;
1062 using band_type = typename rep::band_type;
1063 using vector_element_type = typename rep::vector_element_type;
1064
1065// allocator:
1066
1068 const space_type& Xh,
1069 const Expr& expr)
1070 : base1(new_macro(rep(Xh,expr))),
1071 base2()
1072 {}
1073
1074// accessors:
1075
1076 const geo_type& get_geo() const { return base1::data().get_geo(); }
1077 const space_type& get_space() const { return base1::data().get_space(); }
1078 bool is_on_band() const { return base1::data().is_on_band(); }
1079 band_type get_band() const { return base1::data().get_band(); }
1080
1081 void initialize (const geo_type& omega_K) const { return base1::data().initialize (omega_K); }
1082
1084 const geo_type& omega_K,
1085 const geo_element& K,
1086 vector_element_type& uk) const
1087 { base1::data().evaluate (omega_K, K, uk); }
1088};
1089// concept;
1090template<class Expr> struct is_field_lazy <field_lazy_terminal_interpolate<Expr> > : std::true_type {};
1091
1092}// namespace details
1093
1094// 2.4.a. general nonlinear expression
1096template<class Expr>
1097typename std::enable_if<
1099 && ! details::has_field_rdof_interface<Expr>::value // TODO: interpolate also rdof with an extended "evaluate" interface
1103 >
1104>::type
1106 const space_basic<
1108 ,typename Expr::memory_type>& Xh,
1109 const Expr& expr)
1110{
1112 return details::field_lazy_terminal_interpolate<wrap_t> (Xh, wrap_t(expr));
1113}
1114// 2.4.b. function & functor
1116template <class Expr>
1117inline
1118typename std::enable_if<
1119 details::is_field_function<Expr>::value
1120 ,details::field_lazy_terminal_interpolate<
1121 details::field_expr_v2_nonlinear_terminal_function<Expr>
1122 >
1123>::type
1133#ifdef TO_CLEAN
1134// 2.4.c. re-interpolation of fields and linear field expressions
1135// for change of mesh, of approx, ect
1136// not truly lazy, but here for the completeness of the interpolate interface
1138template<class T, class M>
1139inline
1140field_basic<T,M>
1141lazy_interpolate (const space_basic<T,M>& X2h, const field_basic<T,M>& u1h)
1142{
1143 return interpolate (X2h, u1h);
1144}
1145
1147template <class T, class M, class Expr>
1148inline
1149typename std::enable_if<
1150 details::is_field_wdof<Expr>::value
1151 && ! details::is_field<Expr>::value
1152 ,field_basic<T,M>
1153>::type
1154lazy_interpolate (const space_basic<T,M>& Xh, const Expr& expr)
1155{
1156 return interpolate (Xh, field_basic<T,M>(expr));
1157}
1158#endif // TO_CLEAN
1159
1160}// namespace rheolef
1161# endif /* _RHEOLEF_FIELD_LAZY_TERMINAL_H */
field::size_type size_type
Definition branch.cc:430
field gh(Float epsilon, Float t, const field &uh, const test &v)
see the band page for the full documentation
const geo_basic< T, M > & level_set() const
Definition band.h:109
band_basic< float_type, memory_type > band_type
void evaluate(const geo_type &omega_K, const geo_element &K, vector_element_type &uk) const
typename field_type::space_type space_type
Eigen::Matrix< scalar_type, Eigen::Dynamic, 1 > vector_element_type
typename float_traits< scalar_type >::type float_type
void evaluate(const geo_type &omega_K, const geo_element &K, vector_element_type &uk) const
Eigen::Matrix< scalar_type, Eigen::Dynamic, 1 > vector_element_type
typename float_traits< scalar_type >::type float_type
field_lazy_terminal_integrate_band_rep(const band_type &gh, const Expr &expr, const integrate_option &iopt)
field_lazy_terminal_integrate_band(const band_type &gh, const Expr &expr, const integrate_option &iopt)
void evaluate(const geo_type &omega_K, const geo_element &K, vector_element_type &uk) const
field_lazy_terminal_integrate_band_rep< Expr > rep
band_basic< float_type, memory_type > band_type
field_lazy_terminal_integrate_rep(const Expr &expr, const integrate_option &iopt)
void evaluate(const geo_type &omega_K, const geo_element &K, vector_element_type &uk) const
Eigen::Matrix< scalar_type, Eigen::Dynamic, 1 > vector_element_type
typename float_traits< scalar_type >::type float_type
void initialize(const geo_type &omega_K) const
void evaluate(const geo_type &omega_K, const geo_element &K, vector_element_type &uk) const
field_lazy_terminal_integrate_rep< Expr > rep
field_lazy_terminal_integrate(const Expr &expr, const integrate_option &iopt)
typename rep::vector_element_type vector_element_type
field_lazy_terminal_integrate(std::string domname, const Expr &expr, const integrate_option &iopt)
field_lazy_terminal_integrate(const geo_type &domain, const Expr &expr, const integrate_option &iopt)
void evaluate(const geo_type &omega_K, const geo_element &K, vector_element_type &uk) const
std::enable_if< is_undeterminated< Value >::value, void >::type evaluate_internal(const geo_type &omega_K, const geo_element &K, vector_element_type &uk) const
std::enable_if<!is_undeterminated< Value >::value, void >::type evaluate_internal(const geo_type &omega_K, const geo_element &K, vector_element_type &uk) const
Eigen::Matrix< scalar_type, Eigen::Dynamic, 1 > vector_element_type
typename float_traits< scalar_type >::type float_type
field_lazy_terminal_interpolate_rep(const space_type &Xh, const Expr &expr)
void evaluate(const geo_type &omega_K, const geo_element &K, vector_element_type &uk) const
field_lazy_terminal_interpolate_rep< Expr > rep
field_lazy_terminal_interpolate(const space_type &Xh, const Expr &expr)
typename rep::vector_element_type vector_element_type
see the disarray page for the full documentation
Definition disarray.h:497
geo_basic< float_type, memory_type > geo_type
Definition field.h:229
const geo_type & get_geo() const
Definition field.h:271
const space_type & get_space() const
Definition field.h:270
space_basic< float_type, memory_type > space_type
Definition field.h:230
abstract base interface class
Definition geo.h:248
generic mesh with rerefence counting
Definition geo.h:1089
see the geo_element page for the full documentation
reference_element::size_type size_type
size_type dimension() const
size_type master(bool i) const
orientation_type get_side_informations(const geo_element &S, size_type &loc_isid, size_type &shift) const
size_type n_subgeo(size_type subgeo_dim) const
size_type dis_ie() const
size_type subgeo_dis_index(size_type subgeo_dim, size_type i) const
see the integrate_option page for the full documentation
see the reference_element page for the full documentation
static const size_type max_side_by_variant
reference_element subgeo(size_type subgeo_dim, size_type loc_isid) const
variant_type variant() const
size_type subgeo_local_vertex(size_type subgeo_dim, size_type loc_isid, size_type loc_jsidvert) const
size_type n_subgeo(size_type subgeo_dim) const
see the tensor3 page for the full documentation
see the tensor4 page for the full documentation
see the tensor page for the full documentation
#define trace_macro(message)
Definition dis_macros.h:111
#define error_macro(message)
Definition dis_macros.h:49
void get_geo(istream &in, my_geo &omega)
Expr1::float_type T
Definition field_expr.h:230
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
#define _RHEOLEF_switch(VALUED, VALUE)
void compute_idof_S2idof_K(const basis_basic< T > &b, const reference_element &hat_K, std::array< std::vector< size_t >, reference_element::max_side_by_variant > &idof_S2idof_K, std::array< size_t, reference_element::max_side_by_variant > &ndof_S, size_t &ndof_K)
This file is part of Rheolef.
field_basic< T, M > lazy_interpolate(const space_basic< T, M > &X2h, const field_basic< T, M > &u1h)
see the interpolate page for the full documentation
Definition field.h:871
std::enable_if< details::is_field_expr_quadrature_arg< Expr >::value, details::field_lazy_terminal_integrate< Expr > >::type lazy_integrate(const typename Expr::geo_type &domain, const Expr &expr, const integrate_option &iopt=integrate_option())
see the integrate page for the full documentation
field_basic< T, M > interpolate(const space_basic< T, M > &V2h, const field_basic< T, M > &u1h)
see the interpolate page for the full documentation
STL namespace.
Definition sphere.icc:25
helper for std::complex<T>: get basic T type
Definition Float.h:93
Expr1::memory_type M