Rheolef  7.2
an efficient C++ finite element environment
space_constitution_assembly.cc
Go to the documentation of this file.
1 //
22 // helper for dof numbering with space product and component access
23 // continuation of space_constitution.cc
24 //
25 // contents:
26 // 8. dos_idof, for assembly
27 //
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"
34 #include <sstream>
35 
36 namespace rheolef {
37 
38 // ============================================================================
39 // 8. dis_idof, for assembly
40 // ============================================================================
41 template <class T, class M>
44 {
45  // lazy evaluated on the fly:
46  if ( _loc_ndof [hat_K.variant()] != std::numeric_limits<size_type>::max()) {
47  return _loc_ndof [hat_K.variant()];
48  }
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()];
54  }
55  // mixed multi-component:
56  size_type loc_ndof = 0;
57  for (const_iterator iter = _hier_constit.begin(), last = _hier_constit.end(); iter != last; ++iter) {
58  const space_constitution<T,M>& constit = *iter;
59  loc_ndof += constit.loc_ndof(hat_K); // recursive call
60  }
61  _loc_ndof [hat_K.variant()] = loc_ndof;
62  return loc_ndof;
63 }
64 // DG:
65 // integrate ("internal_sides", jump(u)*average(v));
66 // => in that case, have double dofs on internal sides
67 template <class T, class M>
68 static
69 bool is_dg_not_broken (
70  const space_constitution_rep<T,M>& constit,
71  const geo_basic<T,M>& omega_K,
72  const geo_element& K)
73 {
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());
77  // K is a side on omega_K and geo_space is_broken: two domains have the "is_broken" flag: "sides" and "internal_sides"
78  bool is_broken = constit.get_geo().is_broken(); // e.g. HDG lagrange multipliers
79  bool res = (is_dg && !is_broken);
80  return res;
81 }
82 // HDG:
83 // space Mh(omega["sides"],"P1d"); trial lambda(Mh);
84 // integrate (omega, on_local_sides(lambda));
85 // => in that case, have to loop on sides
86 template <class T, class M>
87 static
88 bool is_broken (
89  const space_constitution_rep<T,M>& constit,
90  const geo_basic<T,M>& omega_K)
91 {
92  // only two domains have the "is_broken" flag: "sides" and "internal_sides"
93  bool res = constit.get_geo().is_broken() && (omega_K.map_dimension() == constit.get_geo().map_dimension() + 1);
94  return res;
95 }
96 // HDG-POST:
97 // space Mhs(omega,"P1d[sides]"); trial lambda_s(Mhs);
98 // integrate (omega, on_local_sides(lambda_s));
99 // => in that case, have to loop on sides with bi-valued lambda_s
100 template <class T, class M>
101 static
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)
106 {
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;
110 }
111 // assembly loop on geo_elements K on a domain omega_K
112 // give its corresponding element K1 in space.geo
113 template <class T, class M>
114 const geo_element&
116  const geo_basic<T,M>& space_geo,
117  const geo_basic<T,M>& omega_K,
118  const geo_element& K_in)
119 {
120  bool is_broken = space_geo.is_broken() && (omega_K.map_dimension() == space_geo.map_dimension() + 1);
121  bool use_dom2bgd = (omega_K.variant() == geo_abstract_base_rep<T>::geo_domain &&
122  space_geo.variant() != geo_abstract_base_rep<T>::geo_domain);
123  bool use_bgd2dom = (omega_K.variant() != geo_abstract_base_rep<T>::geo_domain &&
124  space_geo.variant() == geo_abstract_base_rep<T>::geo_domain);
125  if (!is_broken && use_bgd2dom) { // get from background mesh:
126  return space_geo.bgd2dom_geo_element (K_in);
127  } else if (!is_broken && use_dom2bgd) { // get on background mesh:
128  return omega_K.dom2bgd_geo_element (K_in);
129  } else { // usual case: use K directly
130  return K_in;
131  }
132 }
133 // assembly_loc_ndof
134 // assembly_dis_idof : take into account DG on internal sides with doubled loc_ndof
135 template <class T, class M>
138  const geo_basic<T,M>& omega_K,
139  const geo_element& K_in) const
140 {
141  if (! is_hierarchical()) {
142 trace_macro("loc_ndof: omega_K="<<omega_K.name()<<", K_in="<<K_in.name()<<K_in.dis_ie()<<", space="<<name()<<"...");
143  const geo_element& K = assembly2space_geo_element (get_geo(), omega_K, K_in);
144  size_type loc_ndof = 0;
145  if (is_dg_not_broken (*this, omega_K, K_in)) {
146 trace_macro("loc_ndof: is_dg_not_broken...");
147  // dg: check for internal sides with doubled dofs
148  size_type map_d = K.dimension();
149  size_type dis_ie0 = K.master(0),
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()) {
155  // boundary side K
156  const geo_element& L = get_geo().dis_get_geo_element (map_d+1, dis_ie0);
157  loc_ndof = this->loc_ndof(L);
158  } else {
159  // internal side K
160  const geo_element& L0 = get_geo().dis_get_geo_element (map_d+1, dis_ie0);
161  const geo_element& L1 = get_geo().dis_get_geo_element (map_d+1, dis_ie1);
162  loc_ndof = this->loc_ndof(L0) + this->loc_ndof(L1);
163  }
164  } else if (is_dg_restricted_to_interface (*this, omega_K, K_in)) {
165 trace_macro("loc_ndof: is_dg_restricted_to_interface...");
166  // HDG-POST multiplier: check for boundary sides dofs
167  // example: K=e, b="P1d[sides]", omega_K="square" & get_geo="square"
168  size_type map_d = K.dimension();
169  size_type dis_ie0 = K.master(0),
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()) {
174  // boundary side K
175  const geo_element& L0 = get_geo().dis_get_geo_element (map_d+1, dis_ie0);
177  L0.get_side_informations (K, sid0);
178  loc_ndof = get_basis().local_ndof_on_side (L0, sid0);
179  } else {
180  // internal side K
181  const geo_element& L0 = get_geo().dis_get_geo_element (map_d+1, dis_ie0);
182  const geo_element& L1 = get_geo().dis_get_geo_element (map_d+1, dis_ie1);
183  side_information_type sid0, sid1;
184  L0.get_side_informations (K, sid0);
185  L1.get_side_informations (K, sid1);
186  loc_ndof = get_basis().local_ndof_on_side (L0, sid0)
187  + get_basis().local_ndof_on_side (L1, sid1);
188  // TODO: how to manage e.g. lambda_hs["interface"] = 1; it involves bi-valued side-based space (Pkd(S))^2
189  fatal_macro ("HDG-POST multipliplier: cannot be extracted on an internal interface (bi-valued)");
190  }
191  } else if (is_broken (*this, omega_K)) {
192 trace_macro("loc_ndof: is_broken...");
193  // HDG multiplier: check for boundary sides dofs
194  // example: omega_K="square" & get_geo="square[sides]"
195  size_type space_map_d = get_geo().map_dimension();
196  check_macro (space_map_d+1 == K.dimension(), "internal error: unexpected element dimension");
197  reference_element hat_K = K;
198  for (size_type isid = 0, nsid = hat_K.n_side(); isid < nsid; ++isid) {
199  reference_element hat_S = hat_K.side(isid);
200  loc_ndof += get_basis().ndof(hat_S);
201  }
202  } else {
203 trace_macro("loc_ndof: usual...");
204  // usual case:
205  loc_ndof = this->loc_ndof (K);
206  }
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");
208  return loc_ndof;
209  }
210  // hybrid multi-component:
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");
213  size_type loc_ndof = 0;
214  for (const_iterator iter = _hier_constit.begin(), last = _hier_constit.end(); iter != last; ++iter) {
215  const space_constitution<T,M>& constit = *iter;
216  loc_ndof += constit.assembly_loc_ndof(omega_K, K_in); // recursive call
217  }
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");
219  return loc_ndof;
220 }
221 template <class T, class M>
222 void
224  const geo_basic<T,M>& omega_K,
225  const geo_element& K_in,
226  typename std::vector<geo_element::size_type>::iterator& dis_idof_t,
227  const distributor& hier_ownership,
228  const std::vector<distributor>& start_by_flattened_component,
229  size_type& i_flat_comp) const
230 {
231  size_type space_map_d = get_geo().map_dimension();
232  size_type map_d = K_in.dimension();
233  if (! is_hierarchical()) {
234  const geo_element& K = assembly2space_geo_element (get_geo(), omega_K, K_in);
235  size_type loc_flat_comp_ndof = 0;
236  if (is_dg_not_broken (*this, omega_K, K_in)) {
237  // DG case:
238  size_type dis_ie0 = K.master(0),
239  dis_ie1 = K.master(1);
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()) {
243  // boundary side K with one adjacent element
244  const geo_element& L = get_geo().dis_get_geo_element (map_d+1, dis_ie0);
245  space_numbering::dis_idof (get_basis(), get_geo().sizes(), L, dis_idof_t);
246  loc_flat_comp_ndof = get_basis().ndof (L);
247  } else {
248  // internal side K with two adjacent elements
249  const geo_element& L0 = get_geo().dis_get_geo_element (map_d+1, dis_ie0);
250  const geo_element& L1 = get_geo().dis_get_geo_element (map_d+1, dis_ie1);
251  space_numbering::dis_idof (get_basis(), get_geo().sizes(), L0, dis_idof_t);
252  size_type loc_flat_comp_ndof0 = get_basis().ndof (L0);
253  space_numbering::dis_idof (get_basis(), get_geo().sizes(), L1, dis_idof_t + loc_flat_comp_ndof0);
254  loc_flat_comp_ndof = loc_flat_comp_ndof0 + get_basis().ndof (L1);
255  }
256  } else if (is_dg_restricted_to_interface (*this, omega_K, K_in)) {
257  // HDG-POST multiplier: check for boundary sides dofs
258  // example: K=e, b="P1d[sides]", omega_K="square" & get_geo="square"
259  size_type dis_ie0 = K.master(0),
260  dis_ie1 = K.master(1);
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()) {
264  // boundary side K with one adjacent element
265  const geo_element& L0 = get_geo().dis_get_geo_element (map_d+1, dis_ie0);
267  L0.get_side_informations (K, sid0);
268  std::vector<size_type> dis_idof_L0;
269  space_numbering::dis_idof (get_basis(), get_geo().sizes(), L0, 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];
275  }
276  loc_flat_comp_ndof = loc_idof_L0.size();
277  } else {
278  // internal side K with two adjacent elements
279  const geo_element& L0 = get_geo().dis_get_geo_element (map_d+1, dis_ie0);
280  const geo_element& L1 = get_geo().dis_get_geo_element (map_d+1, dis_ie1);
281  side_information_type sid0, sid1;
282  L0.get_side_informations (K, sid0);
283  L1.get_side_informations (K, sid1);
284  std::vector<size_type> dis_idof_L0, dis_idof_L1;
285  space_numbering::dis_idof (get_basis(), get_geo().sizes(), L0, dis_idof_L0);
286  space_numbering::dis_idof (get_basis(), get_geo().sizes(), L1, 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];
295  }
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)");
298  }
299  } else if (is_broken (*this, omega_K)) {
300  // HDG multiplier: check for boundary sides dofs
301  // example: omega_K="square" & get_geo="square[sides]"
302  check_macro (space_map_d+1 == K.dimension(), "internal error: unexpected element dimension");
303  for (size_type isid = 0, nsid = K_in.n_subgeo(space_map_d); isid < nsid; ++isid) {
304  size_type dis_isid_in = (K_in.dimension() == 1) ? K_in[isid] : (K_in.dimension() == 2) ? K_in.edge(isid) : K_in.face(isid);
305  const geo_element& S_in = omega_K.dis_get_geo_element (space_map_d, dis_isid_in);
306  bool have_bgd_dom = (get_geo().variant() == geo_abstract_base_rep<T>::geo_domain && get_geo().name() != omega_K.name());
307  const geo_element* S_ptr = (!have_bgd_dom) ? &S_in : &(get_geo().bgd2dom_geo_element(S_in));
308  const geo_element& S = *S_ptr;
309  space_numbering::dis_idof (get_basis(), get_geo().sizes(), S, dis_idof_t + loc_flat_comp_ndof);
310  loc_flat_comp_ndof += get_basis().ndof (S);
311  }
312  } else {
313  // non-broken cases: dg or usual
314  // usual: non-DG case:
315  space_numbering::dis_idof (get_basis(), get_geo().sizes(), K, dis_idof_t);
316  loc_flat_comp_ndof = get_basis().ndof (K);
317  }
318  // shift dis_idof depending on i_flat_comp
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];
321  distributor comp_ownership = ownership();
322  size_type iproc = comp_ownership.find_owner (dis_flat_comp_idof);
323  size_type first_dis_flat_comp_idof = comp_ownership.first_index(iproc);
324  size_type first_dis_idof = hier_ownership.first_index(iproc);
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;
329  size_type dis_idof = first_dis_idof + idof;
330  dis_idof_t [loc_flat_comp_idof] = dis_idof;
331  }
332  dis_idof_t += loc_flat_comp_ndof;
333  i_flat_comp++;
334  return;
335  }
336  // hierarchical case
337  for (size_type i_comp = 0, n_comp = size(); i_comp < n_comp; i_comp++) {
338  const space_constitution<T,M>& comp_constit = operator[] (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); // recursive call
340  }
341 }
342 template <class T, class M>
343 void
345  const geo_basic<T,M>& omega_K,
346  const geo_element& K_in,
347  std::vector<geo_element::size_type>& dis_idof_t) const
348 {
349  size_type i_flat_comp = 0;
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);
354 }
355 // ----------------------------------------------------------------------------
356 // instanciation in library
357 // ----------------------------------------------------------------------------
359 
360 #ifdef _RHEOLEF_HAVE_MPI
362 #endif // _RHEOLEF_HAVE_MPI
363 
364 } // namespace rheolef
see the distributor page for the full documentation
Definition: distributor.h:69
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 first_index(size_type iproc) const
global index range and local size owned by ip-th process
Definition: distributor.h:158
abstract base interface class
Definition: geo.h:248
see the geo_element page for the full documentation
Definition: geo_element.h:102
size_type edge(size_type i) const
Definition: geo_element.h:209
size_type face(size_type i) const
Definition: geo_element.h:210
size_type dimension() const
Definition: geo_element.h:167
size_type master(bool i) const
Definition: geo_element.h:165
orientation_type get_side_informations(const geo_element &S, size_type &loc_isid, size_type &shift) const
Definition: geo_element.cc:185
size_type n_subgeo(size_type subgeo_dim) const
Definition: geo_element.h:212
size_type dis_ie() const
Definition: geo_element.h:163
char name() const
Definition: geo_element.h:169
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)
Definition: dis_macros.h:111
#define fatal_macro(message)
Definition: dis_macros.h:33
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)