Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
space_constitution.cc
Go to the documentation of this file.
1
21//
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
43namespace rheolef {
44
45// ======================================================================
46// 1. allocators
47// ======================================================================
48template <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}
73template <class T, class M>
74void
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}
83template <class T, class M>
84void
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// -----------------------------------------------------------------------
102template <class T, class M>
103void
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}
175template <class T, class M>
176void
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}
194template <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// ======================================================================
207template <class T, class M>
208const 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}
232template <class T, class M>
233const 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}
245template <class T, class M>
246const 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// ======================================================================
255template <class T, class M>
256bool
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}
269template <class T, class M>
270bool
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
284template <class T, class M>
285inline
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}
306template <class T, class M>
307bool
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// =======================================================================
324template <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}
336template <class T, class M>
337void
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}
360template <class T, class M>
361void
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}
381template <class T, class M>
382void
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// ---------------------------------------------------------------,--------
400static
401void
402local_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}
429template <class T, class M>
430void
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();
450 variant < reference_element:: last_variant_by_dimension (map_d); ++variant) {
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}
459template <class T, class M>
460void
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}
478template <class T, class M>
479void
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// ======================================================================
490template<class T, class M>
491void
492space_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}
509template<class T, class M>
510std::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// ======================================================================
520template <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}
534template <class T, class M>
535void
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// ======================================================================
553template <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}
571template <class T, class M>
572void
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}
579template <class T, class M>
580void
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}
616template <class T, class M>
617void
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// ======================================================================
650template <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)
void set_dimension(size_type d)
size_type dimension() const
valued_type valued_tag() const
static std::string standard_naming(std::string family_name, size_t degree, const basis_option &sopt)
Definition basis_rep.cc:44
see the disarray page for the full documentation
Definition disarray.h:497
see the distributor page for the full documentation
Definition distributor.h:69
size_type last_index(size_type iproc) const
size_type find_owner(size_type dis_i) const
find iproc associated to a global index dis_i: CPU=log(nproc)
size_type dis_size() const
global and local sizes
size_type size(size_type iproc) const
size_type first_index(size_type iproc) const
global index range and local size owned by ip-th process
bool is_owned(size_type dis_i, size_type iproc) const
true when dis_i in [first_index(iproc):last_index(iproc)[
static const size_type decide
Definition distributor.h:83
std::allocator< int >::size_type size_type
Definition distributor.h:74
the finite element boundary domain
generic mesh with rerefence counting
Definition geo.h:1089
see the geo_element page for the full documentation
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
the finite element space
Definition space.h:382
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 basis_basic< T > & get_basis() const
const geo_basic< T, M > & get_geo() const
bool have_compact_support_inside_element() 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 geo_basic< T, M > & get_geo() const
const distributor & ownership() const
void do_act(const space_act &act)
rep::const_iterator const_iterator
Definition space_mult.h:61
#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.
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
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
Expr1::memory_type M