Rheolef  7.2
an efficient C++ finite element environment
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 
39 namespace rheolef {
40 
41 #ifdef TO_CLEAN
42 // -------------------------------------------------------------------
43 // 1. concept & base class
44 // -------------------------------------------------------------------
45 namespace 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
52 template<class Derived>
53 class field_lazy_base {
54 public:
55 
56 // no common methods yet
57 
58 protected:
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 //
74 namespace details {
75 
76 template<class T, class M>
77 class field_lazy_terminal_field: public field_lazy_base<field_lazy_terminal_field<T,M>> {
78 public :
79 // definitions:
80 
83  using memory_type = M;
84  using scalar_type = T;
87  using geo_type = typename field_type::geo_type;
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:
110 protected:
112 };
113 // concept;
114 template<class T, class M>
115 struct is_field_lazy <field_lazy_terminal_field<T,M> > : std::true_type {};
116 
117 // inlined;
118 template<class T, class M>
119 void
121 {
122  _uh.dis_dof_update();
123 }
124 template<class T, class M>
125 void
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 //
150 namespace details {
151 
152 template<class Expr>
154 public :
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:
196 protected:
198  mutable Expr _expr;
204 // internal:
205  const geo_type& get_default_geo () const;
206 };
207 // inlined;
208 template<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:
223 template<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:
238 template<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 }
253 template<class Expr>
256 {
257  // TODO: improve the automatic determination of the domain ?
258  return get_space().get_constitution().get_geo();
259 }
260 template<class Expr>
261 void
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 */
364 template <class T>
365 void
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 }
468 template<class Expr>
469 void
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 
517 template<class Expr>
518 class field_lazy_terminal_integrate: public smart_pointer_nocopy<field_lazy_terminal_integrate_rep<Expr>>,
519  public field_lazy_base <field_lazy_terminal_integrate<Expr>> {
520 public :
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;
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;
575 template<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
580 template <class Expr>
581 inline
582 typename
583 std::enable_if<
586 >::type
589  const typename Expr::geo_type& domain,
590  const Expr& expr,
591  const integrate_option& iopt = integrate_option())
592 {
594 }
595 template <class Expr>
596 inline
597 typename
598 std::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
612 template <class Expr>
613 inline
614 typename
615 std::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 }
626 template <class Expr>
627 inline
628 typename
629 std::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
642 template <class Expr>
643 inline
644 typename
645 std::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 }
657 template <class Expr>
658 inline
659 typename
660 std::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 // ----------------------------------------------
676 namespace details {
677 
678 template<class Expr>
680 public :
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:
713 protected:
715  mutable Expr _expr;
720 };
721 // inlined;
722 template<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 }
735 template<class Expr>
736 void
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 }
750 template<class Expr>
751 void
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 }
766 template<class Expr>
767 class 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>> {
769 public :
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;
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;
809 template<class Expr> struct is_field_lazy <field_lazy_terminal_integrate_band<Expr> > : std::true_type {};
810 
811 }// namespace details
812 
813 template <class Expr>
814 typename
815 std::enable_if<
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 }
828 template <class Expr>
829 typename
830 std::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 // ----------------------------------------------
848 namespace details {
849 
850 template<class Expr>
852 public :
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;
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;
884 protected:
885 // internals:
886 
887  template <class Value>
888  typename std::enable_if<
890  ,void
891  >::type
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;
917 template<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 }
930 template<class Expr>
931 void
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 }
953 template<class Expr>
954 template <class Value>
955 typename 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 }
978 template<class Expr>
979 template <class Value>
980 typename 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;
1002 #ifdef TODO
1005 #endif // TODO
1006 #undef _RHEOLEF_switch
1007  default: error_macro ("unexpected argument valued tag="<<_Xh.valued_tag());
1008  }
1009 }
1010 template<class Expr>
1011 void
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 }
1048 template<class Expr>
1049 class field_lazy_terminal_interpolate: public smart_pointer_nocopy<field_lazy_terminal_interpolate_rep<Expr>>,
1050  public field_lazy_base <field_lazy_terminal_interpolate<Expr>> {
1051 public :
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;
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 
1083  void evaluate (
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;
1090 template<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
1096 template<class Expr>
1097 typename 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 {
1113 }
1114 // 2.4.b. function & functor
1116 template <class Expr>
1117 inline
1118 typename 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
1125  const space_basic<
1128  const Expr& expr)
1129 {
1132 }
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
1138 template<class T, class M>
1139 inline
1140 field_basic<T,M>
1141 lazy_interpolate (const space_basic<T,M>& X2h, const field_basic<T,M>& u1h)
1142 {
1143  return interpolate (X2h, u1h);
1144 }
1145 
1147 template <class T, class M, class Expr>
1148 inline
1149 typename std::enable_if<
1150  details::is_field_wdof<Expr>::value
1151  && ! details::is_field<Expr>::value
1152  ,field_basic<T,M>
1153 >::type
1154 lazy_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
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(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)
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
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
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
typename float_traits< scalar_type >::type float_type
field_lazy_terminal_interpolate_rep(const space_type &Xh, const Expr &expr)
space_basic< float_type, memory_type > space_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_interpolate(const space_type &Xh, const Expr &expr)
typename rep::vector_element_type vector_element_type
const space_type & get_space() const
Definition: field.h:270
const geo_type & get_geo() const
Definition: field.h:271
geo_basic< float_type, memory_type > geo_type
Definition: field.h:229
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
Definition: geo_element.h:102
reference_element::size_type size_type
Definition: geo_element.h:125
size_type dimension() const
Definition: geo_element.h:167
size_type master(bool i) const
Definition: geo_element.h:165
orientation_type get_side_informations(const geo_element &S, size_type &loc_isid, size_type &shift) const
Definition: geo_element.cc:185
size_type n_subgeo(size_type subgeo_dim) const
Definition: geo_element.h:212
size_type dis_ie() const
Definition: geo_element.h:163
char name() const
Definition: geo_element.h:169
size_type subgeo_dis_index(size_type subgeo_dim, size_type i) const
Definition: geo_element.h:214
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
typename scalar_traits< value_type >::type scalar_type
rheolef::std type
size_t size_type
Definition: basis_get.cc:76
rheolef::std value
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)
void dis_idof(const basis_basic< T > &b, const geo_size &gs, const geo_element &K, typename std::vector< size_type >::iterator dis_idof_tab)
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
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
Definition: interpolate.cc:233
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
Definition: sphere.icc:25
helper for std::complex<T>: get basic T type
Definition: Float.h:93
Expr1::memory_type M
Definition: vec_expr_v2.h:416