Rheolef  7.1
an efficient C++ finite element environment
space_constitution.cc
Go to the documentation of this file.
1 //
22 // helper for dof numbering with space product and component access
23 //
24 // contents:
25 // 1. allocators
26 // 2. accessors
27 // 3. inquire
28 // 4. initialize
29 // 5. output
30 // 6. do_act: block & unblock on a domain
31 // 7. ios accessors
32 // 8. dis_idof, for assembly: moved to space_constitution_assembly.cc
33 //
34 #include "rheolef/geo_domain.h"
35 #include "rheolef/space.h"
36 #include "rheolef/space_numbering.h"
37 #include "rheolef/space_mult.h"
38 #include "rheolef/piola_util.h"
39 #include "rheolef/basis_get.h"
40 #include <sstream>
41 
42 namespace rheolef {
43 
44 // ======================================================================
45 // 1. allocators
46 // ======================================================================
47 template <class T, class M>
49  const geo_basic<T,M>& omega,
50  std::string approx)
51  : _acts(),
52  _omega(compact(omega)),
53  _fem_basis()
54 {
55  if (approx == "") return;
56  const size_type unset = std::numeric_limits<basis_option::size_type>::max();
58  basis_parse_from_string (approx, fio);
59  if (fio.option.valued_tag() != space_constant::scalar && fio.option.dimension() == unset) {
60  // when vector/tensor on surface geometry: fix the number of components by forcing the physical dimension
61  // e.g. 3D-vector on a surface (map_dimension = 2)
62  fio.option.set_dimension (omega.dimension());
63  }
65  omega.coordinate_system() != space_constant::cartesian) {
66  // when tensor on an axisymmetric geometry: should add a theta-theta tensor component
67  fio.option.set_coordinate_system (omega.coordinate_system());
68  }
69  std::string approx_d = basis_rep<T>::standard_naming (fio.family, fio.index, fio.option);
70  _fem_basis = basis_basic<T>(approx_d);
71 }
72 template <class T, class M>
73 void
75 {
76  check_macro (get_geo().map_dimension() < get_geo().dimension() ||
77  !get_basis().have_compact_support_inside_element(),
78  "try to [un]block domain `" << act.get_domain_name() << "' on discontinuous or bubble space("
79  << get_basis().name()<<"(" <<_omega.name()<<"))");
80  _acts.push_back (act);
81 }
82 template <class T, class M>
83 void
85 {
87  if (! is_hierarchical()) {
88  space_scalar_constitution<T,M>& scalar_constit = get_scalar();
89  scalar_constit.do_act (act);
90  return;
91  }
92  // hierarchical case:
93  for (iterator iter = _hier_constit.begin(), last = _hier_constit.end(); iter != last; ++iter) {
94  space_constitution<T,M>& constit = *iter;
95  constit.do_act (act); // recursive call
96  }
97 }
98 // -----------------------------------------------------------------------
99 // loop on domains: mark blocked dofs (called by space::freeze)
100 // -----------------------------------------------------------------------
101 template <class T, class M>
102 void
104  disarray<size_type,M>& blocked_flag, // disarray<bool,M> not supported
105  const distributor& comp_ownership,
106  const distributor& start_by_flattened_component) const
107 {
108  check_macro (get_basis().is_initialized(), "space with undefined finite element method cannot be used");
109  std::vector<size_type> comp_dis_idof_t;
110  distributor dof_ownership = blocked_flag.ownership();
111  size_type dis_ndof = dof_ownership.dis_size();
112  bool prev_act = space_act::block;
113  size_type n_comp = get_basis().size();
114  size_type dim = get_geo().dimension();
115  for (const_iterator iter = begin(), last = end(); iter != last; iter++) {
116  typename space_act::act_type act = (*iter).get_act();
117  const std::string& dom_name = (*iter).get_domain_name();
118  size_type i_comp = (*iter).get_component_index();
119  check_macro (i_comp == space_act::unset_index || i_comp < n_comp,
120  "invalid blocked domain for the "<<i_comp << "-th component in a "
121  <<n_comp<<"-component non-hierarchical space with basis=\""
122  << get_basis().name() << "\"");
123  // sync all blocked_flags when there is an alternance of block and unblock (bug fixed)
124  if (prev_act != act) {
125  blocked_flag.dis_entry_assembly();
126  prev_act = act;
127  }
128  const domain_indirect_basic<M>& dom = _omega.get_domain_indirect(dom_name);
129  size_type dom_dim = dom.map_dimension();
130  distributor ige_ownership = _omega.geo_element_ownership (dom_dim);
131  size_type first_dis_ige = ige_ownership.first_index();
132  size_type last_dis_ige = ige_ownership.last_index();
133  bool blk = (act == space_act::block || act == space_act::block_n) ? true : false;
134  for (size_type ioige = 0, noige = dom.size(); ioige < noige; ioige++) {
135  size_type ige = dom.oige (ioige).index();
136  size_type dis_ige = ige + first_dis_ige;
137  const geo_element& S = _omega.get_geo_element (dom_dim, ige);
138  space_numbering::dis_idof (get_basis(), _omega.sizes(), S, comp_dis_idof_t);
139  size_type loc_idof_start, loc_idof_incr;
140  if (act == space_act::block || act == space_act::unblock) {
141  loc_idof_start = (i_comp == space_act::unset_index) ? 0 : i_comp;
142  loc_idof_incr = (i_comp == space_act::unset_index) ? 1 : n_comp;
143  } else { // block_n or unblock_n
144  loc_idof_start = dim-1;
145  loc_idof_incr = dim;
146  }
147  // vector-valued : do not work yet with RTk, but only with (Pk)^d
148  // note that RTk Hdiv-conform could block_n in theory
149  for (size_type loc_idof = loc_idof_start, loc_ndof = comp_dis_idof_t.size(); loc_idof < loc_ndof; loc_idof += loc_idof_incr) {
150  size_type comp_dis_idof = comp_dis_idof_t [loc_idof];
151  size_type iproc = comp_ownership.find_owner (comp_dis_idof);
152  size_type first_comp_dis_idof = comp_ownership.first_index (iproc);
153  assert_macro (comp_dis_idof >= first_comp_dis_idof, "unexpected comp_dis_idof");
154  size_type comp_idof = comp_dis_idof - first_comp_dis_idof;
155  size_type comp_start_idof = start_by_flattened_component.size(iproc);
156  size_type idof = comp_start_idof + comp_idof;
157  size_type first_dis_idof = dof_ownership.first_index(iproc);
158  size_type dis_idof = first_dis_idof + idof;
159  assert_macro (dis_idof < dis_ndof, "unexpected dis_idof");
160  blocked_flag.dis_entry (dis_idof) = blk;
161  }
162 #ifdef TODO
163  basis scalar_basis = ...; // scalar Pk
164  space_numbering::dis_idof (scalar_basis, _omega.sizes(), S, scalar_dis_idof_t);
165  for (size_type loc_scalar_idof = 0, loc_scalar_ndof = comp_dis_idof_t.size()/n_comp; loc_scalar_idof < loc_scalar_ndof; ++loc_scalar_idof) {
166  size_type dis_scalar_idof = scalar_dis_idof_t [loc_idof]; // as if we do a scalar Pk element
167  has_nt_basis.dis_entry (dis_scalar_idof) = (act == space_act::block_n || act == space_act::unblock_n);
168  normal.dis_entry (dis_scalar_idof) = copy_of_a_previous_computation;
169  }
170 #endif // TODO
171  }
172  }
173 }
174 template <class T, class M>
175 void
177  disarray<size_type,M>& blocked_flag, // disarray<bool,M> not supported
178  const std::vector<distributor>& start_by_flattened_component,
179  size_type& i_flat_comp) const
180 {
181  if (! is_hierarchical()) {
182  const space_scalar_constitution<T,M>& scalar_constit = get_scalar();
183  scalar_constit.data().build_blocked_flag (blocked_flag, ownership(), start_by_flattened_component [i_flat_comp]);
184  i_flat_comp++;
185  return;
186  }
187  // hierarchical case:
188  for (const_iterator iter = _hier_constit.begin(), last = _hier_constit.end(); iter != last; ++iter) {
189  const space_constitution<T,M>& constit = *iter;
190  constit.data().build_blocked_flag_recursive (blocked_flag, start_by_flattened_component, i_flat_comp); // recursive call
191  }
192 }
193 template <class T, class M>
196 {
197  disarray<size_type,M> blocked_flag (ownership(), 0); // disarray<bool,M> not supported
198  size_type i_flat_comp = 0;
199  build_blocked_flag_recursive (blocked_flag, _start_by_flattened_component, i_flat_comp);
200  blocked_flag.dis_entry_assembly();
201  return blocked_flag;
202 }
203 // ======================================================================
204 // 2. accessors
205 // ======================================================================
206 template <class T, class M>
207 const geo_basic<T,M>&
209 {
210  if (! is_hierarchical()) {
211  return get_scalar().get_geo();
212  }
213  // heterogeneous: get as geo the maximum domain of definition
214  // => it is the intersection of all domain definition
215  // implementation note: use "const geo*" ptr_dom for returning a "const geo&"
216  check_macro (_hier_constit.size() > 0, "get_geo: empty space product");
217  const geo_basic<T,M>* ptr_dom = &(_hier_constit[0].get_geo()); // recursive call
218  for (const_iterator iter = _hier_constit.begin(), last = _hier_constit.end(); iter != last; ++iter) {
219  const space_constitution<T,M>& constit = *iter;
220  const geo_basic<T,M>* ptr_comp_dom = &(constit.get_geo()); // recursive call
221  if (ptr_comp_dom->name() == ptr_dom->name()) continue;
222  if (ptr_dom->get_background_geo().name() == ptr_comp_dom->name()) continue;
223  if (ptr_comp_dom->get_background_geo().name() == ptr_dom->name()) {
224  ptr_dom = ptr_comp_dom;
225  } else {
226  error_macro ("get_geo: incompatible domains: \""<<ptr_dom->name()<<"\" and \""<<ptr_comp_dom->name() << "\"");
227  }
228  }
229  return *(ptr_dom);
230 }
231 template <class T, class M>
232 const geo_basic<T,M>&
234 {
235  if (! is_hierarchical()) {
236  return get_scalar().get_background_geo();
237  }
238  // component have the same background mesh: assume it without check (TODO: constit.check_have_same_bgd_geo())
239  check_macro (_hier_constit.size() > 0, "get_background_geo: empty space product");
240  const_iterator iter = _hier_constit.begin();
241  const space_constitution<T,M>& constit = (*iter);
242  return constit.get_background_geo(); // recursive call
243 }
244 template <class T, class M>
245 const basis_basic<T>&
247 {
248  check_macro(! is_hierarchical(), "get_basis: undefined for heterogeneous space products");
249  return get_scalar().get_basis();
250 }
251 // ======================================================================
252 // 3. inquire
253 // ======================================================================
254 template <class T, class M>
255 bool
257 {
258  if (! is_hierarchical()) {
259  return get_scalar().get_basis().have_compact_support_inside_element();
260  }
261  bool has = true;
262  for (const_iterator iter = _hier_constit.begin(), last = _hier_constit.end(); iter != last; ++iter) {
263  const space_constitution<T,M>& constit = *iter;
264  has = has && constit.have_compact_support_inside_element(); // recursive call
265  }
266  return has;
267 }
268 template <class T, class M>
269 bool
271 {
272  if (! is_hierarchical()) {
273  return get_scalar().get_basis().is_discontinuous();
274  }
275  bool is = true;
276  for (const_iterator iter = _hier_constit.begin(), last = _hier_constit.end(); iter != last; ++iter) {
277  const space_constitution<T,M>& constit = *iter;
278  is = is && constit.is_discontinuous(); // recursive call
279  }
280  return is;
281 }
282 template <class T, class M>
284  const geo_basic<T,M>& omega,
285  std::string approx)
286  : _is_initialized(false),
287  _flattened_size(0),
288  _start_by_flattened_component(),
289  _start_by_component(),
290  _ownership(),
291  _valued_tag(space_constant::scalar),
292  _is_hier(false),
293  _scalar_constit(),
294  _hier_constit(),
295  _loc_ndof()
296 {
297  _loc_ndof.fill (std::numeric_limits<size_type>::max());
298  _scalar_constit = scalar_type (omega, approx);
299  if (_scalar_constit.get_basis().is_initialized()) {
300  _valued_tag = _scalar_constit.get_basis().valued_tag();
301  }
302  initialize();
303 }
304 // build a space::constitution hierrachy from a list of spaces
305 template <class T, class M>
307  : _is_initialized(false),
308  _flattened_size(0),
309  _start_by_flattened_component(),
310  _start_by_component(),
311  _ownership(),
312  _valued_tag(space_constant::mixed),
313  _is_hier(true),
314  _scalar_constit(),
315  _hier_constit(expr.size()),
316  _loc_ndof()
317 {
318  _loc_ndof.fill (std::numeric_limits<size_type>::max());
320  for (typename space_mult_list<T,M>::const_iterator iter = expr.begin(), last = expr.end(); iter != last; ++iter, ++iter_constit) {
321  const space_basic<T,M>& Xi = *iter;
322  *iter_constit = Xi.get_constitution();
323  }
324  initialize();
325 }
326 template <class T, class M>
327 bool
329 {
330  if (_is_hier != V2._is_hier) { return false; }
331  if (!_is_hier) { return (_scalar_constit == V2._scalar_constit); }
332  // here, two hierarchies:
333  if (_hier_constit.size() != V2._hier_constit.size()) return false;
334  for (const_iterator iter1 = _hier_constit.begin(), iter2 = V2._hier_constit.begin(), last1 = _hier_constit.end(); iter1 != last1; ++iter1, ++iter2) {
335  if (! (*iter1).data().operator==((*iter2).data())) { // recursive call
336  return false;
337  }
338  }
339  return true;
340 }
341 // =======================================================================
342 // 4. initialize
343 // =======================================================================
344 template <class T, class M>
347 {
348  if (! is_hierarchical()) return 1;
349  size_type flattened_size = 0;
350  for (size_type i_comp = 0, n_comp = size(); i_comp < n_comp; i_comp++) {
351  const space_constitution<T,M>& comp_constit = operator[] (i_comp);
352  flattened_size += comp_constit.data()._init_flattened_size(); // recursive call
353  }
354  return flattened_size;
355 }
356 template <class T, class M>
357 void
359  size_type& i_flat_comp,
360  size_type& start_flat_comp_idof,
361  size_type& dis_start_flat_comp_idof,
362  std::vector<distributor>& start_by_flattened_component) const
363 {
364  if (! is_hierarchical()) {
365  if (! get_basis().is_initialized()) return; // empty numbering => empty space
366  start_by_flattened_component [i_flat_comp] = distributor (dis_start_flat_comp_idof, comm(), start_flat_comp_idof);
367  size_type ndof = space_numbering:: ndof (get_basis(), get_geo().sizes(), get_geo().map_dimension());
368  size_type dis_ndof = space_numbering::dis_ndof (get_basis(), get_geo().sizes(), get_geo().map_dimension());
369  i_flat_comp++;
370  start_flat_comp_idof += ndof;
371  dis_start_flat_comp_idof += dis_ndof;
372  return;
373  }
374  // hierarchical case:
375  for (size_type i_comp = 0, n_comp = size(); i_comp < n_comp; i_comp++) {
376  const space_constitution<T,M>& comp_constit = operator[] (i_comp);
377  comp_constit.data()._init_start_by_flattened_component (i_flat_comp, start_flat_comp_idof, dis_start_flat_comp_idof, start_by_flattened_component); // recursive call
378  }
379 }
380 template <class T, class M>
381 void
383 {
384  if (! is_hierarchical()) {
385  if (! get_basis().is_initialized()) return; // empty numbering => empty space
386  _start_by_component.resize (1);
387  _start_by_component [0] = distributor (0, comm(), 0);
388  return;
389  }
390  // hierarchical case:
391  _start_by_component.resize (size());
392  size_type comp_start_idof = 0;
393  size_type comp_start_dis_idof = 0;
394  for (size_type i_comp = 0, n_comp = size(); i_comp < n_comp; i_comp++) {
395  const space_constitution<T,M>& comp_constit = operator[] (i_comp);
396  _start_by_component [i_comp] = distributor (comp_start_dis_idof, comm(), comp_start_idof);
397  comp_start_idof += comp_constit.ownership(). size();
398  comp_start_dis_idof += comp_constit.ownership().dis_size();
399  }
400 }
401 template <class T, class M>
402 void
404 {
405  if (_is_initialized) return;
406  _is_initialized = true;
407 
408  _flattened_size = _init_flattened_size();
409  _start_by_flattened_component.resize (_flattened_size);
410  size_type i_flat_comp = 0;
411  size_type start_flat_comp_idof = 0;
412  size_type dis_start_flat_comp_idof = 0;
413  _init_start_by_flattened_component (i_flat_comp, start_flat_comp_idof, dis_start_flat_comp_idof, _start_by_flattened_component);
414  _ownership = distributor (dis_start_flat_comp_idof, comm(), start_flat_comp_idof);
415  _init_start_by_component();
416 }
417 // ----------------------------------------------------------------------
418 // get external idofs: used by space<distributed>::freeze
419 // ---------------------------------------------------------------,--------
420 static
421 void
422 local_append_external_dis_idof (
423  const std::vector<size_t>& comp_dis_idof_tab,
424  const distributor& comp_ownership,
425  const distributor& dof_ownership,
426  const distributor& start_by_flattened_component,
427  std::set<size_t>& ext_dof_set)
428 {
429  using size_type = size_t;
430  size_type comp_dis_ndof = comp_ownership.dis_size();
431  for (size_type loc_idof = 0, loc_ndof = comp_dis_idof_tab.size(); loc_idof < loc_ndof; loc_idof++) {
432  size_type comp_dis_idof = comp_dis_idof_tab [loc_idof];
433  assert_macro (comp_dis_idof < comp_dis_ndof, "idof " << comp_dis_idof_tab [loc_idof] << " out of range[0:"<< comp_dis_ndof << "[");
434  if (! comp_ownership.is_owned (comp_dis_idof)) {
435  size_type iproc = comp_ownership.find_owner (comp_dis_idof);
436  size_type comp_first_dis_idof = comp_ownership.first_index(iproc);
437  assert_macro (comp_dis_idof >= comp_first_dis_idof, "unexpected comp_dis_idof");
438  size_type comp_idof = comp_dis_idof - comp_first_dis_idof;
439  size_type comp_start_idof = start_by_flattened_component.size(iproc);
440  size_type idof = comp_start_idof + comp_idof;
441  size_type first_dis_idof = dof_ownership.first_index(iproc);
442  size_type dis_idof = first_dis_idof + idof;
443  assert_macro (! dof_ownership.is_owned (dis_idof), "unexpected dis_idof="<<dis_idof<<" in range ["
444  << first_dis_idof << ":" << dof_ownership.last_index() << "[");
445  ext_dof_set.insert (dis_idof);
446  }
447  }
448 }
449 template <class T, class M>
450 void
452  const geo_basic<T,M>& omega,
453  std::set<size_type>& ext_dof_set,
454  const distributor& dof_ownership,
455  const distributor& start_by_flattened_component) const
456 {
457  // assume here a scalar (non-multi valued) case:
458  distributor comp_ownership = ownership();
459  size_type comp_dis_ndof = comp_ownership.dis_size();
460  std::vector<size_type> comp_dis_idof_tab;
461  // loop on internal elements
462  for (size_type ie = 0, ne = omega.size(); ie < ne; ie++) {
463  const geo_element& K = omega [ie];
464  assembly_dis_idof (omega, K, comp_dis_idof_tab);
465  local_append_external_dis_idof (comp_dis_idof_tab, comp_ownership, dof_ownership, start_by_flattened_component, ext_dof_set);
466  }
467  // loop also on external elements
468  size_type map_d = omega.map_dimension();
471  // for (auto iter : omega._geo_element[variant].get_dis_map_entries()) {
472  for (auto iter : omega.get_external_geo_element_map(variant)) {
473  const geo_element& K = iter.second;
474  assembly_dis_idof (omega, K, comp_dis_idof_tab);
475  local_append_external_dis_idof (comp_dis_idof_tab, comp_ownership, dof_ownership, start_by_flattened_component, ext_dof_set);
476  }
477  }
478 }
479 template <class T, class M>
480 void
482  std::set<size_type>& ext_dof_set,
483  const distributor& dof_ownership,
484  const std::vector<distributor>& start_by_flattened_component,
485  size_type& i_flat_comp) const
486 {
487  if (! is_hierarchical()) {
488  const space_scalar_constitution<T,M>& scalar_constit = get_scalar();
489  append_external_dof (scalar_constit.get_geo(), ext_dof_set, dof_ownership, start_by_flattened_component [i_flat_comp]);
490  i_flat_comp++;
491  return;
492  }
493  // hierarchical case:
494  for (const_iterator iter = _hier_constit.begin(), last = _hier_constit.end(); iter != last; ++iter) {
495  (*iter).data().compute_external_dofs (ext_dof_set, dof_ownership, start_by_flattened_component, i_flat_comp);
496  }
497 }
498 template <class T, class M>
499 void
501  std::set<size_type>& ext_dof_set) const
502 {
503  ext_dof_set.clear();
504  size_type i_flat_comp = 0;
505  compute_external_dofs (ext_dof_set, ownership(), _start_by_flattened_component, i_flat_comp);
506 }
507 // ======================================================================
508 // 5. output
509 // ======================================================================
510 template<class T, class M>
511 void
512 space_constitution_rep<T,M>::put (std::ostream& out, size_type level) const
513 {
514  if (! is_hierarchical()) {
515  const space_scalar_constitution<T,M>& scalar_constit = get_scalar();
516  if (!scalar_constit.get_basis().is_initialized()) return;
517  out << scalar_constit.get_basis().name() << "{" << scalar_constit.get_geo().name() << "}";
518  } else {
519  if (level > 0) out << "(";
520  typename space_constitution<T,M>::hierarchy_type::const_iterator x = get_hierarchy().begin();
521  for (size_type i = 0, n = get_hierarchy().size(); i < n; i++) {
522  const space_constitution<T,M>& xi = x[i];
523  xi.data().put (out, level+1); // recursive call
524  if (i+1 < n) { out << "*"; }
525  }
526  if (level > 0) out << ")";
527  }
528 }
529 template<class T, class M>
530 std::string
532 {
533  std::ostringstream ostrstr;
534  put (ostrstr, 0);
535  return ostrstr.str();
536 }
537 // ======================================================================
538 // 6. do_act: block & unblock on a domain
539 // ======================================================================
540 template <class T, class M>
543 {
544  if (! is_hierarchical()) {
545  return get_scalar().get_basis().degree();
546  }
547  size_type degree_max = 0;
548  for (const_iterator iter = _hier_constit.begin(), last = _hier_constit.end(); iter != last; ++iter) {
549  const space_constitution<T,M>& constit = *iter;
550  degree_max = std::max(degree_max, constit.degree_max()); // recursive call
551  }
552  return degree_max;
553 }
554 template <class T, class M>
555 void
557 {
558  if (! is_hierarchical()) {
559  if (get_geo().dimension() == get_geo().map_dimension()) {
560  // 3d surface mesh or 2d lineique mesh are not yet treated for neighbours
561  get_scalar().neighbour_guard();
562  }
563  return;
564  }
565  for (const_iterator iter = _hier_constit.begin(), last = _hier_constit.end(); iter != last; ++iter) {
566  const space_constitution<T,M>& constit = *iter;
567  constit.neighbour_guard(); // recursive call
568  }
569 }
570 // ======================================================================
571 // 7. ios accessors
572 // ======================================================================
573 template <class T, class M>
576 {
577  if (!_is_hier) {
578  communicator comm =_scalar_constit.get_geo().comm();
579  size_type dis_ndof = space_numbering::dis_ndof (_scalar_constit.get_basis(), _scalar_constit.get_geo().sizes(), _scalar_constit.get_geo().map_dimension());
580  // TODO: memorize ios_ownership : avoid comms & risks of errors
581  distributor ios_ownership (dis_ndof, comm, distributor::decide);
582  return ios_ownership.size();
583  }
584  size_type ios_ndof = 0;
585  for (const_iterator iter = _hier_constit.begin(), last = _hier_constit.end(); iter != last; ++iter) {
586  const space_constitution<T,M>& constit = *iter;
587  ios_ndof += constit.ios_ndof(); // recursive call
588  }
589  return ios_ndof;
590 }
591 template <class T, class M>
592 void
594  disarray<size_type,M>& idof2ios_dis_idof,
595  disarray<size_type,M>& ios_idof2dis_idof) const
596 {
597  space_numbering::set_ios_permutations (get_basis(), get_geo(), idof2ios_dis_idof, ios_idof2dis_idof);
598 }
599 template <class T, class M>
600 void
602  disarray<size_type,M>& idof2ios_dis_idof,
603  size_type& comp_start_idof,
604  size_type& comp_start_dis_idof) const
605 {
606  if (!_is_hier) {
607  // non-hierarchical case:
608  disarray<size_type,M> comp_idof2ios_dis_idof;
609  disarray<size_type,M> comp_ios_idof2dis_idof;
610  // TODO: avoid allocating comp_idof2ios_dis_idof: numbering supports directly comp_start_{dis_}idof
611  // TODO: avoid invert permut here: numbering support no invert perm arg
612  _scalar_constit.set_ios_permutations (comp_idof2ios_dis_idof, comp_ios_idof2dis_idof);
613  // then shift
614  size_type comp_ndof = comp_idof2ios_dis_idof.size();
615  size_type comp_dis_ndof = comp_idof2ios_dis_idof.dis_size();
616  size_type ndof = idof2ios_dis_idof.size();
617  size_type dis_ndof = idof2ios_dis_idof.dis_size();
618  for (size_type comp_idof = 0; comp_idof < comp_ndof; comp_idof++) {
619  size_type comp_ios_dis_idof = comp_idof2ios_dis_idof [comp_idof];
620  size_type ios_dis_idof = comp_start_dis_idof + comp_ios_dis_idof;
621  size_type idof = comp_start_idof + comp_idof;
622  assert_macro (idof < ndof, "unexpected idof="<<idof<<" out of range [0:"<<ndof<<"[");
623  assert_macro (ios_dis_idof < dis_ndof, "unexpected ios_dis_idof="<<ios_dis_idof<<" out of range [0:"<<dis_ndof<<"[");
624  idof2ios_dis_idof [idof] = ios_dis_idof;
625  }
626  comp_start_idof += comp_ndof;
627  comp_start_dis_idof += comp_dis_ndof;
628  return;
629  }
630  // hierarchical case:
631  for (const_iterator iter = _hier_constit.begin(), last = _hier_constit.end(); iter != last; ++iter) {
632  const space_constitution<T,M>& constit = *iter;
633  constit.set_ios_permutation_recursion (idof2ios_dis_idof, comp_start_idof, comp_start_dis_idof); // recursive call
634  }
635 }
636 template <class T, class M>
637 void
639  disarray<size_type,M>& idof2ios_dis_idof,
640  disarray<size_type,M>& ios_idof2dis_idof) const
641 {
642  if (!_is_hier) {
643  // non-hierarchical case: direct
644  _scalar_constit.set_ios_permutations (idof2ios_dis_idof, ios_idof2dis_idof);
645  return;
646  }
647  // hierarchical case:
648  // ----------------------------------------------
649  // 1) compute permutation
650  // ----------------------------------------------
651  communicator comm1 = comm();
652  distributor dof_ownership = ownership();
653  size_type ndof1 = dof_ownership.size();
654  size_type dis_ndof1 = dof_ownership.dis_size();
655  idof2ios_dis_idof.resize (dof_ownership, std::numeric_limits<size_type>::max());
656  size_type comp_start_idof = 0;
657  size_type comp_start_dis_idof = 0;
658  set_ios_permutation_recursion (idof2ios_dis_idof, comp_start_idof, comp_start_dis_idof);
659  // ----------------------------------------------
660  // 2) invert permutation into ios_idof2dis_idof
661  // ----------------------------------------------
662  size_type ios_ndof1 = ios_ndof();
663  distributor ios_dof_ownership (dis_ndof1, idof2ios_dis_idof.comm(), ios_ndof1);
664  ios_idof2dis_idof.resize (ios_dof_ownership, std::numeric_limits<size_type>::max());
665  idof2ios_dis_idof.reverse_permutation (ios_idof2dis_idof);
666 }
667 // ----------------------------------------------------------------------------
668 // instanciation in library
669 // ----------------------------------------------------------------------------
673 
674 #ifdef _RHEOLEF_HAVE_MPI
678 #endif // _RHEOLEF_HAVE_MPI
679 
680 } // namespace rheolef
rheolef::space_constitution_rep::_init_flattened_size
size_type _init_flattened_size() const
Definition: space_constitution.cc:346
rheolef::space_numbering::ndof
size_type ndof(const basis_basic< T > &b, const geo_size &gs, size_type map_dim)
Definition: space_numbering.cc:28
rheolef::reference_element::last_variant_by_dimension
static variant_type last_variant_by_dimension(size_type dim)
Definition: reference_element.h:150
rheolef::compact
geo_basic< T, M > compact(const geo_basic< T, M > &gamma)
Definition: geo.cc:219
rheolef::space_scalar_constitution
Definition: space_constitution.h:148
rheolef::geo_basic
generic mesh with rerefence counting
Definition: domain_indirect.h:64
rheolef::family_index_option_type::family
std::string family
Definition: basis_get.h:34
rheolef::space_constitution_rep::build_blocked_flag_recursive
void build_blocked_flag_recursive(disarray< size_type, M > &blocked_flag, const std::vector< distributor > &start_by_component, size_type &i_comp) const
Definition: space_constitution.cc:176
rheolef::space_constitution_rep::set_ios_permutations
void set_ios_permutations(disarray< size_type, M > &idof2ios_dis_idof, disarray< size_type, M > &ios_idof2dis_idof) const
Definition: space_constitution.cc:638
rheolef::space_constitution_rep::is_discontinuous
bool is_discontinuous() const
Definition: space_constitution.cc:270
rheolef::space_constitution_rep::set_ios_permutation_recursion
void set_ios_permutation_recursion(disarray< size_type, M > &idof2ios_dis_idof, size_type &comp_start_idof, size_type &comp_start_dis_idof) const
Definition: space_constitution.cc:601
rheolef::io::out
Definition: rheostream.h:167
rheolef::distributor::last_index
size_type last_index(size_type iproc) const
Definition: distributor.h:157
rheolef::space_constitution_rep::have_compact_support_inside_element
bool have_compact_support_inside_element() const
Definition: space_constitution.cc:256
rheolef::space_numbering::dis_ndof
size_type dis_ndof(const basis_basic< T > &b, const geo_size &gs, size_type map_dim)
Definition: space_numbering.cc:41
rheolef::space_constitution_rep::neighbour_guard
void neighbour_guard() const
Definition: space_constitution.cc:556
mkgeo_ball.expr
expr
Definition: mkgeo_ball.sh:361
rheolef::space_numbering::set_ios_permutations
void set_ios_permutations(const basis_basic< T > &b, const geo_basic< T, distributed > &omega, disarray< size_type, distributed > &idof2ios_dis_idof, disarray< size_type, distributed > &ios_idof2dis_idof)
Definition: space_numbering.cc:275
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
rheolef::space_constitution_rep::degree_max
size_type degree_max() const
Definition: space_constitution.cc:542
rheolef::put
void put(std::ostream &out, std::string name, const tiny_matrix< T > &a)
Definition: tiny_lu.h:155
rheolef::space_constitution::is_discontinuous
bool is_discontinuous() const
Definition: space_constitution.h:437
rheolef::family_index_option_type::index
size_t index
Definition: basis_get.h:35
rheolef::space_scalar_constitution_rep::set_ios_permutations
void set_ios_permutations(disarray< size_type, M > &idof2ios_dis_idof, disarray< size_type, M > &ios_idof2dis_idof) const
Definition: space_constitution.cc:593
rheolef::normal
details::field_expr_v2_nonlinear_terminal_function< details::normal_pseudo_function< Float > > normal()
normal: see the expression page for the full documentation
Definition: field_expr_terminal.h:439
rheolef::space_constitution_rep::get_background_geo
const geo_basic< T, M > & get_background_geo() const
Definition: space_constitution.cc:233
rheolef::space_scalar_constitution_rep::const_iterator
container_type::const_iterator const_iterator
Definition: space_constitution.h:85
rheolef::space_scalar_constitution_rep
Definition: space_constitution.h:81
rheolef::distributor::find_owner
size_type find_owner(size_type dis_i) const
find iproc associated to a global index dis_i: CPU=log(nproc)
Definition: distributor.cc:106
rheolef::space_constitution_rep::ios_ndof
size_type ios_ndof() const
Definition: space_constitution.cc:575
rheolef::space_scalar_constitution::do_act
void do_act(const space_act &act)
Definition: space_constitution.h:197
rheolef::space_act::get_domain_name
const std::string & get_domain_name() const
Definition: space_constitution.h:66
rheolef::space_constant::cartesian
Definition: space_constant.h:122
rheolef::space_constitution_rep::_init_start_by_component
void _init_start_by_component() const
Definition: space_constitution.cc:382
rheolef::space_scalar_constitution_rep::space_scalar_constitution_rep
space_scalar_constitution_rep()
Definition: space_constitution.h:89
rheolef::space_constitution_rep::_init_start_by_flattened_component
void _init_start_by_flattened_component(size_type &i_flat_comp, size_type &start_flat_comp_idof, size_type &dis_start_flat_comp_idof, std::vector< distributor > &start_by_flattened_component) const
Definition: space_constitution.cc:358
rheolef::space_basic
the finite element space
Definition: space.h:352
rheolef::distributor
see the distributor page for the full documentation
Definition: distributor.h:62
rheolef::space_act::act_type
act_type
Definition: space_constitution.h:46
rheolef::geo_element
see the geo_element page for the full documentation
Definition: geo_element.h:102
rheolef::space_act
Definition: space_constitution.h:38
rheolef::space_constant::mixed
Definition: space_constant.h:142
rheolef::space_act::block
Definition: space_constitution.h:47
rheolef::distributor::decide
static const size_type decide
Definition: distributor.h:76
rheolef::space_constitution_rep::put
void put(std::ostream &out, size_type level=0) const
Definition: space_constitution.cc:512
rheolef::size_type
size_t size_type
Definition: basis_get.cc:76
rheolef::space_numbering::dis_idof
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)
Definition: space_numbering.cc:187
rheolef::space_constitution::do_act
void do_act(const space_act &act)
Definition: space_constitution.h:424
rheolef::reference_element::first_variant_by_dimension
static variant_type first_variant_by_dimension(size_type dim)
Definition: reference_element.h:148
rheolef::basis_basic
Definition: basis.h:206
rheolef::space_constant::tensor
Definition: space_constant.h:138
mkgeo_ball.variant
variant
Definition: mkgeo_ball.sh:149
rheolef::space_constitution_rep::const_iterator
hierarchy_type::const_iterator const_iterator
Definition: space_constitution.h:216
rheolef::space_constant::scalar
Definition: space_constant.h:136
rheolef::space_constitution_rep::iterator
hierarchy_type::iterator iterator
Definition: space_constitution.h:215
rheolef::space_constitution_rep::get_geo
const geo_basic< T, M > & get_geo() const
Definition: space_constitution.cc:208
rheolef::space_constitution_rep::_valued_tag
valued_type _valued_tag
Definition: space_constitution.h:334
rheolef::space_scalar_constitution::get_basis
const basis_basic< T > & get_basis() const
Definition: space_constitution.h:170
rheolef::space_constitution_rep::get_basis
const basis_basic< T > & get_basis() const
Definition: space_constitution.cc:246
basis
see the basis page for the full documentation
rheolef::space_constitution::get_background_geo
const geo_basic< T, M > & get_background_geo() const
Definition: space_constitution.h:417
rheolef::space_constitution::neighbour_guard
void neighbour_guard() const
Definition: space_constitution.h:439
rheolef::space_scalar_constitution_rep::build_blocked_flag
void build_blocked_flag(disarray< size_type, M > &blocked_flag, const distributor &comp_ownership, const distributor &start_by_component) const
Definition: space_constitution.cc:103
assert_macro
#define assert_macro(ok_condition, message)
Definition: dis_macros.h:113
rheolef::space_scalar_constitution::get_geo
const geo_basic< T, M > & get_geo() const
Definition: space_constitution.h:168
rheolef::space_constitution_rep::do_act
void do_act(const space_act &act)
Definition: space_constitution.cc:84
rheolef::space_constitution::have_compact_support_inside_element
bool have_compact_support_inside_element() const
Definition: space_constitution.h:436
rheolef::space_scalar_constitution_rep::size_type
container_type::size_type size_type
Definition: space_constitution.h:84
rheolef::space_constitution_rep::size_type
hierarchy_type::size_type size_type
Definition: space_constitution.h:214
rheolef::space_constitution::ownership
const distributor & ownership() const
Definition: space_constitution.h:410
rheolef::space_constitution::degree_max
size_type degree_max() const
Definition: space_constitution.h:438
rheolef::family_index_option_type
Definition: basis_get.h:31
dimension
const size_t dimension
Definition: edge.icc:64
rheolef::space_mult_list
Definition: space.h:111
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::distributor::dis_size
size_type dis_size() const
global and local sizes
Definition: distributor.h:207
rheolef::basis_option::set_dimension
void set_dimension(size_type d)
Definition: basis_option.h:150
error_macro
#define error_macro(message)
Definition: dis_macros.h:49
rheolef::space_scalar_constitution_rep::do_act
void do_act(const space_act &act)
Definition: space_constitution.cc:74
rheolef::basis_option::set_coordinate_system
void set_coordinate_system(coordinate_type sc)
Definition: basis_option.h:151
rheolef::smart_pointer_base< space_scalar_constitution_rep< T, M >, details::constructor_copy< space_scalar_constitution_rep< T, M > > >::data
const space_scalar_constitution_rep< T, M > & data() const
Definition: smart_pointer.h:266
rheolef::basis_option::dimension
size_type dimension() const
Definition: basis_option.h:138
rheolef::space_constitution::get_geo
const geo_basic< T, M > & get_geo() const
Definition: space_constitution.h:416
rheolef::space_constitution::set_ios_permutation_recursion
void set_ios_permutation_recursion(disarray< size_type, M > &idof2ios_dis_idof, size_type &comp_start_idof, size_type &comp_start_dis_idof) const
Definition: space_constitution.h:474
rheolef::space_act::unset_index
static const size_type unset_index
Definition: space_constitution.h:44
rheolef::disarray< size_type, M >
rheolef::space_constitution_rep::_is_hier
bool _is_hier
Definition: space_constitution.h:335
rheolef::space_constitution_rep::space_constitution_rep
space_constitution_rep()
Definition: space_constitution.h:346
rheolef::space_constitution_rep::operator==
bool operator==(const space_constitution_rep< T, M > &V2) const
Definition: space_constitution.cc:328
rheolef::family_index_option_type::option
basis_option option
Definition: basis_get.h:36
rheolef::space_act::block_n
Definition: space_constitution.h:49
rheolef::space_constitution_rep::scalar_type
space_scalar_constitution< T, M > scalar_type
Definition: space_constitution.h:213
rheolef::space_constitution_rep::_loc_ndof
std::array< size_type, reference_element::max_variant > _loc_ndof
Definition: space_constitution.h:338
rheolef::space_scalar_constitution_rep::_fem_basis
basis_basic< T > _fem_basis
Definition: space_constitution.h:145
rheolef::basis_rep::standard_naming
static std::string standard_naming(std::string family_name, size_t degree, const basis_option &sopt)
Definition: basis_rep.cc:44
rheolef::space_constitution_rep::build_blocked_flag
disarray< size_type, M > build_blocked_flag() const
Definition: space_constitution.cc:195
mkgeo_ball.n
n
Definition: mkgeo_ball.sh:150
rheolef::space_constitution_rep::compute_external_dofs
void compute_external_dofs(std::set< size_type > &ext_dof_set) const
Definition: space_constitution.cc:500
rheolef::basis_parse_from_string
void basis_parse_from_string(const std::string &str, family_index_option_type &fio)
Definition: basis_get.cc:142
size_type
field::size_type size_type
Definition: branch.cc:425
rheolef::space_constitution_rep::_hier_constit
hierarchy_type _hier_constit
Definition: space_constitution.h:337
rheolef::space_constitution_rep::_scalar_constit
scalar_type _scalar_constit
Definition: space_constitution.h:336
rheolef::distributor::is_owned
bool is_owned(size_type dis_i, size_type iproc) const
true when dis_i in [first_index(iproc):last_index(iproc)[
Definition: distributor.h:213
rheolef::distributor::first_index
size_type first_index(size_type iproc) const
global index range and local size owned by ip-th process
Definition: distributor.h:151
get_geo
void get_geo(istream &in, my_geo &omega)
Definition: field2gmsh_pos.cc:126
rheolef::space_mult_list::const_iterator
rep::const_iterator const_iterator
Definition: space_mult.h:61
rheolef::space_constitution_rep::append_external_dof
void append_external_dof(const geo_basic< T, M > &dom, std::set< size_type > &ext_dof_set, const distributor &dof_ownership, const distributor &start_by_component) const
Definition: space_constitution.cc:451
rheolef::space_constitution::ios_ndof
size_type ios_ndof() const
Definition: space_constitution.h:414
rheolef::space_constitution_rep::name
std::string name() const
Definition: space_constitution.cc:531
rheolef::space_act::unblock
Definition: space_constitution.h:48
mkgeo_ball.dim
dim
Definition: mkgeo_ball.sh:307
rheolef::space_constitution_rep::initialize
void initialize() const
Definition: space_constitution.cc:403
rheolef::domain_indirect_basic
the finite element boundary domain
Definition: domain_indirect.h:335
rheolef::distributor::size_type
std::allocator< int >::size_type size_type
Definition: distributor.h:67
mkgeo_contraction.name
name
Definition: mkgeo_contraction.sh:133
M
Expr1::memory_type M
Definition: vec_expr_v2.h:385
rheolef::distributor::size
size_type size(size_type iproc) const
Definition: distributor.h:163
rheolef::space_constitution
Definition: space_constitution.h:32
rheolef::space_act::unblock_n
Definition: space_constitution.h:50
rheolef::basis_option::valued_tag
valued_type valued_tag() const
Definition: basis_option.h:245
rheolef::space_constitution_rep
Definition: space_constitution.h:208