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"
41 template <
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;
67 template <
class T,
class M>
69 bool 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);
86 template <
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);
100 template <
class T,
class M>
102 bool 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_restricted_to_sides()
108 && constit.get_geo().get_background_geo().map_dimension() == omega_K.map_dimension()
109 && omega_K.map_dimension() == K.dimension() + 1;
113 template <
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);
135 template <
class T,
class M>
141 if (! is_hierarchical()) {
145 if (is_dg_not_broken (*
this, omega_K, K_in)) {
150 dis_ie1 = K.master(1);
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()
153 <<
" associated to K="<<K.name()<<K.dis_ie() <<
" in space_geo=" <<
get_geo().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)) {
165 trace_macro(
"loc_ndof: is_dg_restricted_to_interface...");
170 dis_ie1 = K.master(1);
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)) {
196 check_macro (space_map_d+1 == K.dimension(),
"internal error: unexpected element dimension");
198 for (
size_type isid = 0, nsid = hat_K.n_side(); isid < nsid; ++isid) {
200 loc_ndof += get_basis().ndof(hat_S);
205 loc_ndof = this->loc_ndof (K);
207 trace_macro(
"loc_ndof: omega_K="<<omega_K.name()<<
", K_in="<<K_in.
name()<<K_in.
dis_ie()<<
", space="<<
name()<<
": result="<<loc_ndof<<
" done");
211 trace_macro(
"loc_ndof_hier: omega_K="<<omega_K.name()<<
", K_in="<<K_in.
name()<<K_in.
dis_ie()<<
", space="<<
name()<<
"...");
212 check_macro (_is_hier,
"unexpected hierarchical space");
214 for (
const_iterator iter = _hier_constit.begin(), last = _hier_constit.end(); iter != last; ++iter) {
218 trace_macro(
"loc_ndof_hier: omega_K="<<omega_K.name()<<
", K_in="<<K_in.
name()<<K_in.
dis_ie()<<
", space="<<
name()<<
": result="<<loc_ndof<<
" done");
221 template <
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);
342 template <
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
362 #endif // _RHEOLEF_HAVE_MPI