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