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"
48template <
class T,
class M>
56 if (approx ==
"")
return;
57 const size_type unset = std::numeric_limits<basis_option::size_type>::max();
73template <
class T,
class M>
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);
83template <
class T,
class M>
88 if (! is_hierarchical()) {
90 terminal_constit.
do_act (act);
94 for (
iterator iter = _hier_constit.begin(), last = _hier_constit.end(); iter != last; ++iter) {
102template <
class T,
class M>
107 const distributor& start_by_flattened_component)
const
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();
116 for (
const_iterator iter = begin(), last = end(); iter != last; iter++) {
118 const std::string& dom_name = (*iter).get_domain_name();
119 size_type i_comp = (*iter).get_component_index();
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() <<
"\"");
125 if (prev_act != act) {
126 blocked_flag.dis_entry_assembly();
131 distributor ige_ownership = _omega.geo_element_ownership (dom_dim);
135 for (
size_type ioige = 0, noige = dom.size(); ioige < noige; ioige++) {
136 size_type ige = dom.oige (ioige).index();
138 const geo_element& S = _omega.get_geo_element (dom_dim, ige);
145 loc_idof_start = dim-1;
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];
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;
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;
164 basis scalar_basis = ...;
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];
169 normal.dis_entry (dis_scalar_idof) = copy_of_a_previous_computation;
175template <
class T,
class M>
179 const std::vector<distributor>& start_by_flattened_component,
182 if (! is_hierarchical()) {
184 terminal_constit.
data().build_blocked_flag (blocked_flag, ownership(), start_by_flattened_component [i_flat_comp]);
189 for (
const_iterator iter = _hier_constit.begin(), last = _hier_constit.end(); iter != last; ++iter) {
191 constit.
data().build_blocked_flag_recursive (blocked_flag, start_by_flattened_component, i_flat_comp);
194template <
class T,
class M>
200 build_blocked_flag_recursive (blocked_flag, _start_by_flattened_component, i_flat_comp);
201 blocked_flag.dis_entry_assembly();
207template <
class T,
class M>
211 if (! is_hierarchical()) {
212 return get_terminal().get_geo();
217 check_macro (_hier_constit.size() > 0,
"get_geo: empty space product");
219 for (
const_iterator iter = _hier_constit.begin(), last = _hier_constit.end(); iter != last; ++iter) {
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;
227 error_macro (
"get_geo: incompatible domains: \""<<ptr_dom->name()<<
"\" and \""<<ptr_comp_dom->name() <<
"\"");
232template <
class T,
class M>
236 if (! is_hierarchical()) {
237 return get_terminal().get_background_geo();
240 check_macro (_hier_constit.size() > 0,
"get_background_geo: empty space product");
245template <
class T,
class M>
249 check_macro(! is_hierarchical(),
"get_basis: undefined for heterogeneous space products");
250 return get_terminal().get_basis();
255template <
class T,
class M>
259 if (! is_hierarchical()) {
260 return get_terminal().get_basis().have_compact_support_inside_element();
263 for (
const_iterator iter = _hier_constit.begin(), last = _hier_constit.end(); iter != last; ++iter) {
269template <
class T,
class M>
273 if (! is_hierarchical()) {
274 return get_terminal().get_basis().is_discontinuous();
277 for (
const_iterator iter = _hier_constit.begin(), last = _hier_constit.end(); iter != last; ++iter) {
284template <
class T,
class M>
287 : _is_initialized(false),
289 _start_by_flattened_component(),
290 _start_by_component(),
292 _valued_tag(space_constant::mixed),
295 _hier_constit(expr.size()),
298 _loc_ndof.fill (std::numeric_limits<size_type>::max());
302 *iter_constit = Xi.get_constitution();
306template <
class T,
class M>
310 if (_is_hier != V2.
_is_hier) {
return false; }
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())) {
324template <
class T,
class M>
328 if (! is_hierarchical())
return 1;
330 for (
size_type i_comp = 0, n_comp = size(); i_comp < n_comp; i_comp++) {
332 flattened_size += comp_constit.
data()._init_flattened_size();
334 return flattened_size;
336template <
class T,
class M>
342 std::vector<distributor>& start_by_flattened_component)
const
344 if (! is_hierarchical()) {
345 if (! get_basis().is_initialized())
return;
346 start_by_flattened_component [i_flat_comp] =
distributor (dis_start_flat_comp_idof, comm(), start_flat_comp_idof);
350 start_flat_comp_idof += ndof;
351 dis_start_flat_comp_idof += dis_ndof;
355 for (
size_type i_comp = 0, n_comp = size(); i_comp < n_comp; 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);
360template <
class T,
class M>
364 if (! is_hierarchical()) {
365 if (! get_basis().is_initialized())
return;
366 _start_by_component.resize (1);
367 _start_by_component [0] =
distributor (0, comm(), 0);
371 _start_by_component.resize (size());
374 for (
size_type i_comp = 0, n_comp = size(); i_comp < n_comp; 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();
381template <
class T,
class M>
385 if (_is_initialized)
return;
386 _is_initialized =
true;
388 _flattened_size = _init_flattened_size();
389 _start_by_flattened_component.resize (_flattened_size);
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();
402local_append_external_dis_idof (
403 const std::vector<size_t>& comp_dis_idof_tab,
407 std::set<size_t>& ext_dof_set)
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)) {
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;
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);
429template <
class T,
class M>
433 std::set<size_type>& ext_dof_set,
435 const distributor& start_by_flattened_component)
const
440 std::vector<size_type> comp_dis_idof_tab;
442 for (
size_type ie = 0, ne = omega.size(); ie < ne; 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);
452 for (
auto iter : omega.get_external_geo_element_map(variant)) {
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);
459template <
class T,
class M>
462 std::set<size_type>& ext_dof_set,
464 const std::vector<distributor>& start_by_flattened_component,
467 if (! is_hierarchical()) {
469 append_external_dof (terminal_constit.
get_geo(), ext_dof_set, dof_ownership, start_by_flattened_component [i_flat_comp]);
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);
478template <
class T,
class M>
481 std::set<size_type>& ext_dof_set)
const
485 compute_external_dofs (ext_dof_set, ownership(), _start_by_flattened_component, i_flat_comp);
490template<
class T,
class M>
494 if (! is_hierarchical()) {
496 if (!terminal_constit.
get_basis().is_initialized())
return;
497 out << terminal_constit.
get_basis().name() <<
"{" << terminal_constit.
get_geo().name() <<
"}";
499 if (level > 0) out <<
"(";
501 for (
size_type i = 0, n = get_hierarchy().size(); i < n; i++) {
503 xi.
data().put (out, level+1);
504 if (i+1 < n) { out <<
"*"; }
506 if (level > 0) out <<
")";
509template<
class T,
class M>
513 std::ostringstream ostrstr;
515 return ostrstr.str();
520template <
class T,
class M>
524 if (! is_hierarchical()) {
525 return get_terminal().get_basis().degree();
528 for (
const_iterator iter = _hier_constit.begin(), last = _hier_constit.end(); iter != last; ++iter) {
530 degree_max = std::max(degree_max, constit.
degree_max());
534template <
class T,
class M>
538 if (! is_hierarchical()) {
541 get_terminal().neighbour_guard();
545 for (
const_iterator iter = _hier_constit.begin(), last = _hier_constit.end(); iter != last; ++iter) {
553template <
class T,
class M>
558 communicator comm =_terminal_constit.get_geo().comm();
562 return ios_ownership.
size();
565 for (
const_iterator iter = _hier_constit.begin(), last = _hier_constit.end(); iter != last; ++iter) {
571template <
class T,
class M>
579template <
class T,
class M>
592 _terminal_constit.set_ios_permutations (comp_idof2ios_dis_idof, comp_ios_idof2dis_idof);
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;
606 comp_start_idof += comp_ndof;
607 comp_start_dis_idof += comp_dis_ndof;
611 for (
const_iterator iter = _hier_constit.begin(), last = _hier_constit.end(); iter != last; ++iter) {
616template <
class T,
class M>
624 _terminal_constit.set_ios_permutations (idof2ios_dis_idof, ios_idof2dis_idof);
631 communicator comm1 = comm();
635 idof2ios_dis_idof.resize (dof_ownership, std::numeric_limits<size_type>::max());
638 set_ios_permutation_recursion (idof2ios_dis_idof, comp_start_idof, comp_start_dis_idof);
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);
650template <
class T,
class M>
656 if (! is_hierarchical()) {
658 size_type n_comp = get_terminal().get_basis().size();
659 size_type dis_idof = comp_dis_idof*n_comp + i_comp;
663 const distributor& comp_ownership = _hier_constit [i_comp].ownership();
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;
681#ifdef _RHEOLEF_HAVE_MPI
field::size_type size_type
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)
see the disarray page for the full documentation
see the distributor page for the full documentation
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
std::allocator< int >::size_type size_type
the finite element boundary domain
generic mesh with rerefence counting
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
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
hierarchy_type _hier_constit
bool have_compact_support_inside_element() 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
size_type ios_ndof() const
void _init_start_by_component() 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
size_type _init_flattened_size() const
bool is_discontinuous() 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
scalar_type _terminal_constit
size_type degree_max() 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 neighbour_guard() const
void do_act(const space_act &act)
size_type comp_dis_idof2dis_idof(size_type i_comp, size_type comp_dis_idof) const
space_constitution_terminal_rep()
container_type::size_type size_type
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
basis_basic< T > _fem_basis
void do_act(const space_act &act)
const basis_basic< T > & get_basis() const
const geo_basic< T, M > & get_geo() const
void do_act(const space_act &act)
bool have_compact_support_inside_element() const
size_type degree_max() const
const geo_basic< T, M > & get_background_geo() const
bool is_discontinuous() 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
size_type ios_ndof() const
const geo_basic< T, M > & get_geo() const
const distributor & ownership() const
void neighbour_guard() const
void do_act(const space_act &act)
rep::const_iterator const_iterator
#define assert_macro(ok_condition, message)
#define error_macro(message)
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)
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)
void basis_parse_from_string(const std::string &str, family_index_option_type &fio)