28#include "rheolef/geo_domain.h"
29#include "rheolef/space.h"
30#include "rheolef/space_numbering.h"
31#include "rheolef/space_mult.h"
32#include "rheolef/piola_util.h"
33#include "rheolef/basis_get.h"
41template <
class T,
class M>
46 if ( _loc_ndof [hat_K.
variant()] != std::numeric_limits<size_type>::max()) {
47 return _loc_ndof [hat_K.
variant()];
49 if (! is_hierarchical()) {
50 size_type loc_comp_ndof = get_basis().ndof (hat_K);
51 size_type n_comp = (!_is_hier ? 1 : size());
52 _loc_ndof [hat_K.
variant()] = n_comp*loc_comp_ndof;
53 return _loc_ndof [hat_K.
variant()];
57 for (
const_iterator iter = _hier_constit.begin(), last = _hier_constit.end(); iter != last; ++iter) {
61 _loc_ndof [hat_K.
variant()] = loc_ndof;
67template <
class T,
class M>
69bool is_dg_not_broken (
74 bool is_dg = (!constit.
get_basis().is_continuous()
75 && omega_K.map_dimension() + 1 == constit.
get_geo().map_dimension()
76 && omega_K.map_dimension() == K.
dimension());
78 bool is_broken = constit.
get_geo().is_broken();
79 bool res = (is_dg && !is_broken);
86template <
class T,
class M>
89 const space_constitution_rep<T,M>& constit,
90 const geo_basic<T,M>& omega_K)
93 bool res = constit.get_geo().is_broken() && (omega_K.map_dimension() == constit.get_geo().map_dimension() + 1);
100template <
class T,
class M>
102bool is_dg_restricted_to_interface (
103 const space_constitution_rep<T,M>& constit,
104 const geo_basic<T,M>& omega_K,
105 const geo_element& K)
107 return constit.get_basis().option().is_trace_n()
108 && constit.get_geo().get_background_geo().map_dimension() == omega_K.map_dimension()
109 && omega_K.map_dimension() == K.dimension() + 1;
113template <
class T,
class M>
120 bool is_broken = space_geo.is_broken() && (omega_K.map_dimension() == space_geo.map_dimension() + 1);
125 if (!is_broken && use_bgd2dom) {
126 return space_geo.bgd2dom_geo_element (K_in);
127 }
else if (!is_broken && use_dom2bgd) {
128 return omega_K.dom2bgd_geo_element (K_in);
135template <
class T,
class M>
141 if (! is_hierarchical()) {
145 if (is_dg_not_broken (*
this, omega_K, K_in)) {
151 check_macro (dis_ie0 != std::numeric_limits<size_type>::max(),
152 "unexpected isolated side K_in="<<K_in.
name()<<K_in.
dis_ie() <<
" in omega_K=" << omega_K.name()
154 if (dis_ie1 == std::numeric_limits<size_type>::max()) {
157 loc_ndof = this->loc_ndof(
L);
162 loc_ndof = this->loc_ndof(L0) + this->loc_ndof(L1);
164 }
else if (is_dg_restricted_to_interface (*
this, omega_K, K_in)) {
165trace_macro(
"loc_ndof: is_dg_restricted_to_interface...");
171 check_macro (dis_ie0 != std::numeric_limits<size_type>::max(),
172 "unexpected isolated mesh side K.dis_ie="<<K.
dis_ie());
173 if (dis_ie1 == std::numeric_limits<size_type>::max()) {
178 loc_ndof = get_basis().local_ndof_on_side (L0, sid0);
186 loc_ndof = get_basis().local_ndof_on_side (L0, sid0)
187 + get_basis().local_ndof_on_side (L1, sid1);
189 fatal_macro (
"HDG-POST multipliplier: cannot be extracted on an internal interface (bi-valued)");
191 }
else if (is_broken (*
this, omega_K)) {
200 loc_ndof += get_basis().ndof(hat_S);
205 loc_ndof = this->loc_ndof (K);
207trace_macro(
"loc_ndof: omega_K="<<omega_K.name()<<
", K_in="<<K_in.
name()<<K_in.
dis_ie()<<
", space="<<
name()<<
": result="<<loc_ndof<<
" done");
212 check_macro (_is_hier,
"unexpected hierarchical space");
214 for (
const_iterator iter = _hier_constit.begin(), last = _hier_constit.end(); iter != last; ++iter) {
218trace_macro(
"loc_ndof_hier: omega_K="<<omega_K.name()<<
", K_in="<<K_in.
name()<<K_in.
dis_ie()<<
", space="<<
name()<<
": result="<<loc_ndof<<
" done");
221template <
class T,
class M>
226 typename std::vector<geo_element::size_type>::iterator& dis_idof_t,
228 const std::vector<distributor>& start_by_flattened_component,
233 if (! is_hierarchical()) {
236 if (is_dg_not_broken (*
this, omega_K, K_in)) {
240 check_macro (dis_ie0 != std::numeric_limits<size_type>::max(),
241 "unexpected isolated mesh side K.dis_ie="<<K.
dis_ie());
242 if (dis_ie1 == std::numeric_limits<size_type>::max()) {
246 loc_flat_comp_ndof = get_basis().ndof (
L);
252 size_type loc_flat_comp_ndof0 = get_basis().ndof (L0);
254 loc_flat_comp_ndof = loc_flat_comp_ndof0 + get_basis().ndof (L1);
256 }
else if (is_dg_restricted_to_interface (*
this, omega_K, K_in)) {
261 check_macro (dis_ie0 != std::numeric_limits<size_type>::max(),
262 "unexpected isolated mesh side K.dis_ie="<<K.
dis_ie());
263 if (dis_ie1 == std::numeric_limits<size_type>::max()) {
268 std::vector<size_type> dis_idof_L0;
270 Eigen::Matrix<size_type,Eigen::Dynamic,1> loc_idof_L0;
271 get_basis().local_idof_on_side (L0, sid0, loc_idof_L0);
272 for (
size_type loc_sid_idof = 0, loc_sid_ndof = loc_idof_L0.size(); loc_sid_idof < loc_sid_ndof; ++loc_sid_idof) {
273 size_type loc_idof = loc_idof_L0 [loc_sid_idof];
274 dis_idof_t [loc_sid_idof] = dis_idof_L0 [loc_idof];
276 loc_flat_comp_ndof = loc_idof_L0.size();
284 std::vector<size_type> dis_idof_L0, dis_idof_L1;
287 Eigen::Matrix<size_type,Eigen::Dynamic,1> loc_idof_L0, loc_idof_L1;
288 get_basis().local_idof_on_side (L0, sid0, loc_idof_L0);
289 get_basis().local_idof_on_side (L1, sid1, loc_idof_L1);
290 check_macro (loc_idof_L0.size() == loc_idof_L1.size(),
"unexpected bi-valued side: unmatch sizes");
291 for (
size_type loc_sid_idof = 0, loc_sid_ndof = loc_idof_L0.size(); loc_sid_idof < loc_sid_ndof; ++loc_sid_idof) {
292 size_type loc_idof = loc_idof_L0 [loc_sid_idof];
293 dis_idof_t [loc_sid_idof] = dis_idof_L0 [loc_idof];
294 dis_idof_t [loc_sid_ndof+loc_sid_idof] = dis_idof_L1 [loc_idof];
296 loc_flat_comp_ndof = loc_idof_L0.size() + loc_idof_L1.size();
297 fatal_macro (
"HDG-POST multipliplier: cannot be extracted on an internal interface (bi-valued)");
299 }
else if (is_broken (*
this, omega_K)) {
303 for (
size_type isid = 0, nsid = K_in.
n_subgeo(space_map_d); isid < nsid; ++isid) {
305 const geo_element& S_in = omega_K.dis_get_geo_element (space_map_d, dis_isid_in);
307 const geo_element* S_ptr = (!have_bgd_dom) ? &S_in : &(
get_geo().bgd2dom_geo_element(S_in));
310 loc_flat_comp_ndof += get_basis().ndof (S);
316 loc_flat_comp_ndof = get_basis().ndof (K);
319 for (
size_type loc_flat_comp_idof = 0; loc_flat_comp_idof < loc_flat_comp_ndof; ++loc_flat_comp_idof) {
320 size_type dis_flat_comp_idof = dis_idof_t [loc_flat_comp_idof];
325 check_macro (dis_flat_comp_idof >= first_dis_flat_comp_idof,
"unexpected dis_comp_idof");
326 size_type flat_comp_idof = dis_flat_comp_idof - first_dis_flat_comp_idof;
327 size_type start_flat_comp_idof = start_by_flattened_component [i_flat_comp].size(iproc);
328 size_type idof = start_flat_comp_idof + flat_comp_idof;
330 dis_idof_t [loc_flat_comp_idof] =
dis_idof;
332 dis_idof_t += loc_flat_comp_ndof;
337 for (
size_type i_comp = 0, n_comp = size(); i_comp < n_comp; i_comp++) {
339 comp_constit.
data()._assembly_dis_idof_recursive (omega_K, K_in, dis_idof_t, hier_ownership, start_by_flattened_component, i_flat_comp);
342template <
class T,
class M>
347 std::vector<geo_element::size_type>& dis_idof_t)
const
350 size_type loc_ndof = assembly_loc_ndof (omega_K, K_in);
351 dis_idof_t.resize (loc_ndof);
352 typename std::vector<geo_element::size_type>::iterator first_dis_idof_t = dis_idof_t.begin();
353 _assembly_dis_idof_recursive (omega_K, K_in, first_dis_idof_t, ownership(), _start_by_flattened_component, i_flat_comp);
360#ifdef _RHEOLEF_HAVE_MPI
see the distributor page for the full documentation
size_type find_owner(size_type dis_i) const
find iproc associated to a global index dis_i: CPU=log(nproc)
size_type first_index(size_type iproc) const
global index range and local size owned by ip-th process
abstract base interface class
see the geo_element page for the full documentation
size_type edge(size_type i) const
size_type face(size_type i) const
size_type dimension() const
size_type master(bool i) const
orientation_type get_side_informations(const geo_element &S, size_type &loc_isid, size_type &shift) const
size_type n_subgeo(size_type subgeo_dim) const
see the reference_element page for the full documentation
reference_element side(size_type loc_isid) const
variant_type variant() const
const basis_basic< T > & get_basis() const
void assembly_dis_idof(const geo_basic< T, M > &dom, const geo_element &bgd_K, std::vector< geo_element::size_type > &dis_idof) const
void _assembly_dis_idof_recursive(const geo_basic< T, M > &dom, const geo_element &bgd_K, typename std::vector< geo_element::size_type >::iterator &dis_idof_t, const distributor &hier_ownership, const std::vector< distributor > &start_by_flattened_component, size_type &i_flat_comp) const
size_type loc_ndof(const reference_element &hat_K) const
hierarchy_type::const_iterator const_iterator
hierarchy_type::size_type size_type
const geo_basic< T, M > & get_geo() const
size_type assembly_loc_ndof(const geo_basic< T, M > &dom, const geo_element &bgd_K) const
size_type loc_ndof(const reference_element &hat_K) const
size_type assembly_loc_ndof(const geo_basic< T, M > &dom, const geo_element &bgd_K) const
#define trace_macro(message)
#define fatal_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)
This file is part of Rheolef.
const geo_element & assembly2space_geo_element(const geo_basic< T, M > &space_geo, const geo_basic< T, M > &omega_K, const geo_element &K_in)