Rheolef  7.2
an efficient C++ finite element environment
space.cc
Go to the documentation of this file.
1 
22 #include "rheolef/space.h"
23 #include "rheolef/space_numbering.h"
24 #include "rheolef/space_mult.h"
25 #include "rheolef/piola_util.h"
26 #include "rheolef/piola_on_pointset.h"
27 
28 namespace rheolef {
29 
30 // ======================================================================
31 // space_base_rep
32 // ======================================================================
33 template <class T, class M>
35  const geo_basic<T,M>& omega_in,
36  std::string approx,
37  std::string valued)
38 : _constit(omega_in,
39  (approx == "" || valued == "scalar") ? approx : valued+"("+approx+")"),
40  _xdof(),
41  _have_freezed(false),
42  _idof2blk_dis_iub(),
43  _has_nt_basis(),
44  _normal(),
45  _iu_ownership(),
46  _ib_ownership(),
47  _parent_subgeo_owner_initialized(false),
48  _parent_subgeo_owner_reattribution_required(false),
49  _parent_subgeo_owner()
50 {
51  // omega_in is compressed by constitution_terminal() allocator
52  // when it is a domain: then it becomes a geo_domain, with compressed numbering
53  // so, omega_in can be different from omega:
54  if (approx == "") return; // empty element => default cstor
55  _idof2blk_dis_iub.resize (_constit.ownership());
56  _has_nt_basis.resize (_constit.ownership()); // TODO: _constit.scalar_ownership()
57  _normal.resize (_constit.ownership());
58  init_xdof();
59 }
60 template <class T, class M>
61 space_base_rep<T,M>::space_base_rep (
62  const geo_basic<T,M>& omega_in,
63  const basis_basic<T>& b)
64 : _constit(omega_in, b.name()),
65  _xdof(),
66  _have_freezed(false),
67  _idof2blk_dis_iub(),
68  _has_nt_basis(),
69  _normal(),
70  _iu_ownership(),
71  _ib_ownership(),
72  _parent_subgeo_owner_initialized(false),
73  _parent_subgeo_owner_reattribution_required(false),
74  _parent_subgeo_owner()
75 {
76  // omega_in is compressed by constitution_terminal() allocator
77  // when it is a domain: then it becomes a geo_domain, with compressed numbering
78  // so, omega_in can be different from omega:
79  if (! b.is_initialized()) return; // empty element => default cstor
80  _idof2blk_dis_iub.resize (_constit.ownership());
81  _has_nt_basis.resize (_constit.ownership());
82  _normal.resize (_constit.ownership());
83  init_xdof();
84 }
85 template <class T, class M>
86 space_base_rep<T,M>::space_base_rep (const space_constitution<T,M>& constit)
87 : _constit(constit),
88  _xdof(),
89  _have_freezed(false),
90  _idof2blk_dis_iub(),
91  _has_nt_basis(),
92  _normal(),
93  _iu_ownership(),
94  _ib_ownership(),
95  _parent_subgeo_owner_initialized(false),
96  _parent_subgeo_owner_reattribution_required(false),
97  _parent_subgeo_owner()
98 {
99  distributor idof_ownership = _constit.ownership();
100  _idof2blk_dis_iub.resize (idof_ownership);
101  _has_nt_basis.resize (idof_ownership);
102  _normal.resize (idof_ownership);
103  init_xdof();
104 }
105 template <class T, class M>
106 space_base_rep<T,M>::space_base_rep (const space_mult_list<T,M>& expr)
107 : _constit(expr),
108  _xdof(),
109  _have_freezed(false),
110  _idof2blk_dis_iub(),
111  _has_nt_basis(),
112  _normal(),
113  _iu_ownership(),
114  _ib_ownership(),
115  _parent_subgeo_owner_initialized(false),
116  _parent_subgeo_owner_reattribution_required(false),
117  _parent_subgeo_owner()
118 {
119  _idof2blk_dis_iub.resize (_constit.ownership());
120  _has_nt_basis.resize (_constit.ownership());
121  _normal.resize (_constit.ownership());
122  init_xdof();
123 }
124 template <class T, class M>
125 void
126 space_base_rep<T,M>::init_xdof()
127 {
128  if (_constit.is_hierarchical()) {
129  trace_macro ("init_xdof: hierarchical space: xdof undefined");
130  return;
131  }
132  const geo_basic<T,M>& omega = get_geo();
133  const basis_basic<T>& b = get_basis();
134  size_type dis_nnod = space_numbering::dis_nnod (b, omega.sizes(), omega.map_dimension());
135  size_type nnod = space_numbering:: nnod (b, omega.sizes(), omega.map_dimension());
136  communicator comm = ownership().comm();
137  distributor nod_ownership (dis_nnod, comm, nnod);
138  _xdof.resize (nod_ownership);
139  piola_on_pointset<T> pops;
140  pops.initialize (omega.get_piola_basis(), b, integrate_option());
141  std::vector<size_type> dis_inod_tab;
142  for (size_type ie = 0, ne = omega.size(); ie < ne; ie++) {
143  const geo_element& K = omega [ie];
144  const Eigen::Matrix<piola<T>,Eigen::Dynamic,1>& piola = pops.get_piola (omega, K);
145  space_numbering::dis_inod (b, omega.sizes(), K, dis_inod_tab);
146  for (size_type loc_inod = 0, loc_nnod = piola.size(); loc_inod < loc_nnod; ++loc_inod) {
147  _xdof.dis_entry (dis_inod_tab[loc_inod]) = piola [loc_inod].F;
148  }
149  }
150  _xdof.dis_entry_assembly();
151  // add external nodes on current element:
152  if (is_distributed<M>::value && comm.size() > 1) {
153  index_set dis_inod_set;
154  for (size_type ie = 0, ne = omega.size(); ie < ne; ++ie) {
155  const geo_element& K = omega [ie];
156  space_numbering::dis_inod (b, omega.sizes(), K, dis_inod_tab);
157  for (size_type loc_inod = 0, loc_nnod = dis_inod_tab.size(); loc_inod < loc_nnod; ++loc_inod) {
158  size_type dis_inod = dis_inod_tab[loc_inod];
159  if (nod_ownership.is_owned(dis_inod)) continue;
160  check_macro (dis_inod < dis_nnod, "invalid dis_inod="<<dis_inod<<" out of range [0:"<<dis_nnod<<"[");
161  dis_inod_set.insert (dis_inod);
162  }
163  }
164  _xdof.set_dis_indexes (dis_inod_set);
165  }
166 }
167 template <class T, class M>
168 void
169 space_base_rep<T,M>::base_freeze_body () const
170 {
171  _constit.neighbour_guard(); // geo could uses S.master(i), a neighbour connnectivity
172  // -----------------------------------------------------------------------
173  // 1) loop on domains: mark blocked dofs
174  // -----------------------------------------------------------------------
175  disarray<size_type,M> blocked_flag = _constit.build_blocked_flag();
176 
177  // copy the blocked_flag into the dis_iub disarray, as the "blocked" bit:
178  for (size_type idof = 0, ndof = blocked_flag.size(); idof < ndof; idof++) {
179  _idof2blk_dis_iub [idof].set_blocked (blocked_flag[idof]);
180  }
181  // -----------------------------------------------------------------------
182  // 2) init numbering
183  // -----------------------------------------------------------------------
184  size_type n_unknown = 0;
185  size_type n_blocked = 0;
186  for (size_type idof = 0, ndof = _idof2blk_dis_iub.size(); idof < ndof; idof++) {
187  bool blk = _idof2blk_dis_iub [idof].is_blocked();
188  if (! blk) {
189  _idof2blk_dis_iub[idof].set_dis_iub (n_unknown);
190  n_unknown++;
191  } else {
192  _idof2blk_dis_iub[idof].set_dis_iub (n_blocked);
193  n_blocked++;
194  }
195  }
196  size_type dis_n_unknown = n_unknown;
197  size_type dis_n_blocked = n_blocked;
198 #ifdef _RHEOLEF_HAVE_MPI
199  if (is_distributed<M>::value) {
200  dis_n_unknown = mpi::all_reduce (comm(), dis_n_unknown, std::plus<T>());
201  dis_n_blocked = mpi::all_reduce (comm(), dis_n_blocked, std::plus<T>());
202  }
203 #endif // // _RHEOLEF_HAVE_MPI
204  _iu_ownership = distributor (dis_n_unknown, comm(), n_unknown);
205  _ib_ownership = distributor (dis_n_blocked, comm(), n_blocked);
206 }
207 // ----------------------------------------------------------------------------
208 // still used by geo_subdivide.cc and internally in space.cc
209 // ----------------------------------------------------------------------------
210 template <class T, class M>
211 void
212 space_base_rep<T,M>::dis_idof (const geo_element& K, std::vector<size_type>& dis_idof) const
213 {
214  freeze_guard();
215  _constit.assembly_dis_idof (get_geo(), K, dis_idof);
216 }
217 // ----------------------------------------------------------------------------
218 // name : e.g. "P1(square)", for field_expr<Expr> checks
219 // ----------------------------------------------------------------------------
220 template <class T, class M>
221 std::string
223 {
224  return _constit.name();
225 }
226 // ----------------------------------------------------------------------------
227 // u["left"]
228 // => build here the requiresd temporary indirect disarray
229 // ----------------------------------------------------------------------------
239 template <class T, class M>
241 space_base_rep<T,M>::build_dom_dis_idof2bgd_dis_idof (
242  const space_base_rep<T,M>& Wh, // Wh = space on domain
243  const std::string& dom_name // redundant: contained in Wh.get_geo()
244  ) const
245 {
246  const geo_basic<T,M>& bgd_gamma = get_geo()[dom_name];
247  return build_dom_dis_idof2bgd_dis_idof (Wh, bgd_gamma);
248 }
249 template <class T, class M>
251 space_base_rep<T,M>::build_dom_dis_idof2bgd_dis_idof (
252  const space_base_rep<T,M>& Wh, // Wh = space on domain
253  const geo_basic<T,M>& bgd_gamma2
254  ) const
255 {
256  // TODO: move it to the call ?
257  const geo_basic<T,M>& bgd_gamma = bgd_gamma2.get_background_domain();
258 
259  // TODO: verifier que le domaine bgd_gamma est compatible:
260  // => il doit y avoir un meme maillage de base a bgd_gamma & Wh.get_geo
261  const space_base_rep<T,M>& Vh = *this; // Vh = space on background mesh
262  Vh.freeze_guard();
263  Wh.freeze_guard();
264  const geo_basic<T,M>& dom_gamma = Wh.get_geo();
265  size_type map_dim = dom_gamma.map_dimension();
266  check_macro (dom_gamma.size() == bgd_gamma.size(), "incompatible domains");
267  distributor dom_ownership = Wh.ownership();
268  distributor bgd_ownership = Vh.ownership();
269  size_type my_proc = dom_ownership.comm().rank();
270 trace_macro("Vh="<<Vh.name()<<", Wh="<<Wh.name()<<", dom_dis_idof2bgd_dis_idof.size=Wh.size="<<dom_ownership.size()<<"...");
271  size_type first_dom_dis_idof = dom_ownership.first_index();
272  size_type first_bgd_dis_idof = bgd_ownership.first_index();
273  std::vector<size_type> dom_dis_idofs, bgd_dis_idofs;
274  disarray<size_type, M> dom_dis_idof2bgd_dis_idof (dom_ownership, std::numeric_limits<size_type>::max());
275  std::set<size_type> ext_dom_dis_idof;
276  for (size_type ige = 0, nge = dom_gamma.size(); ige < nge; ige++) {
277  const geo_element& dom_S = dom_gamma[ige];
278  const geo_element& bgd_S = bgd_gamma[ige];
279  Wh.dis_idof (dom_S, dom_dis_idofs);
280  Vh.dis_idof (bgd_S, bgd_dis_idofs);
281 trace_macro("ige="<<ige<<": dom_S="<<dom_S.name()<<dom_S.dis_ie()<<", bgd_S="<<bgd_S.name()<<bgd_S.dis_ie()
282  << ": dom_dis_idofs.size="<<dom_dis_idofs.size()
283  << ", bgd_dis_idofs.size="<<bgd_dis_idofs.size());
284  for (size_type loc_idof = 0, loc_ndof = dom_dis_idofs.size(); loc_idof < loc_ndof; loc_idof++) {
285  size_type dom_dis_idof = dom_dis_idofs [loc_idof];
286  size_type bgd_dis_idof = bgd_dis_idofs [loc_idof];
287  dom_dis_idof2bgd_dis_idof.dis_entry (dom_dis_idof) = bgd_dis_idof;
288  trace_macro("ige="<<ige<<": dom_dis_idof2bgd_dis_idof["<<dom_dis_idof<<"] = "<<bgd_dis_idof);
289  if (! dom_ownership.is_owned (dom_dis_idof, my_proc)) {
290  ext_dom_dis_idof.insert (dom_dis_idof);
291  }
292  }
293  }
294  trace_macro ("ext_dom_dis_idof.size="<<ext_dom_dis_idof.size());
295  dom_dis_idof2bgd_dis_idof.dis_entry_assembly();
296  dom_dis_idof2bgd_dis_idof.set_dis_indexes (ext_dom_dis_idof); // doit etre en deuxieme...?
297 #ifdef TO_CLEAN
298  // move to local numbering:
299  for (size_type dom_idof = 0, dom_ndof = dom_dis_idof2bgd_dis_idof.size(); dom_idof < dom_ndof; dom_idof++) {
300  size_type bgd_dis_idof = dom_dis_idof2bgd_dis_idof [dom_idof];
301  size_type dom_dis_idof = dom_idof + first_dom_dis_idof;
302  check_macro (bgd_ownership.is_owned (bgd_dis_idof), "bgd_dis_idof="<<bgd_dis_idof<<" is out of range ["<<first_bgd_dis_idof<<":"<< bgd_ownership.last_index()<<"[");
303  size_type bgd_idof = bgd_dis_idof - first_bgd_dis_idof;
304  dom_dis_idof2bgd_dis_idof [dom_idof] = bgd_idof;
305  }
306 trace_macro("Vh="<<Vh.name()<<", Wh="<<Wh.name()<<", dom_dis_idof2bgd_dis_idof.size=Wh.size="<<dom_ownership.size()<<" done");
307 #endif // TO_CLEAN
308  return dom_dis_idof2bgd_dis_idof;
309 }
310 #define _RHEOLEF_space_real(M) \
311 template <class T> \
312 space_basic<T,M> \
313 space_basic<T,M>::real() \
314 { \
315  return space_basic<T,M> (geo_basic<T,M>(details::zero_dimension()), "P1"); \
316 }
318 #ifdef _RHEOLEF_HAVE_MPI
320 #endif // _RHEOLEF_HAVE_MPI
321 #undef _RHEOLEF_space_real
322 // ======================================================================
323 // parent_subgeo_owner
324 // ======================================================================
325 /* IMPLEMENTATION NOTES:
326 
327  Consider the following situation:
328 
329  P0 S0 top
330  +-----+------+--
331  |\ K0 |
332  S1 | \ |
333  | K1 \|
334  P1 +-----+
335  |
336  left |
337  +
338 
339  vertex P0 is shared by edges S0 and S1 in triangle K0
340  edge S0 belongs to domain "top"
341  edge S1 belongs to domain "left"
342  assume that the mesh partition attributes
343  K0 on partition iproc=0
344  K1 on partition iproc=1
345  then, sides and vertex are also attributed to a proc
346  via a recursive descent on subgeos:
347  iproc=0 K0, S0, P0
348  iproc=1 K1, S1, P1
349  when explorating S1, its vertex P0 is already attributed to iproc=0
350 
351  Then, consider:
352  space Wh (omega["left"], "P1");
353  The element edge S1 contains an **orphean** dof associated to vertex P0 :
354  it do not belong to S1 owner iproc=1 and is not referenced by any others sides
355  of "left" that belongs to its owner iproc=0
356  So, in order to interpolate in a lazy element-by-element loop, we have to
357  re-attribute this dof to S0 owner, iproc=1.
358 
359  notes:
360  - this idof iproc owner reattribution required only when
361  * nproc > 1
362  * continuous elements
363  * omega is a geo_domain, i.e. some parent elements are missing
364  otherwise could return usual dis_idof owner
365 
366  - all dofs are associated to a subgeo
367  - when continuous approx, this subgeo corresponds to a mesh subgeo
368  - when discontinuous, this subgeo de-multiplicated from mesh:
369  * verticies 1D, faces 3D or edges 2D : doubled when internals
370  * edges 3D and verticies 2D and 3D : de-multiplied
371  in all cases, dofs conceptually conserve the memory of the parent subgeo
372 */
373 template <class T, class M>
376 {
377  parent_subgeo_owner_check();
378  if (! _parent_subgeo_owner_reattribution_required) {
379  return ownership().find_owner (dis_idof);
380  }
381  return _parent_subgeo_owner.dis_at (dis_idof);
382 }
383 namespace details {
384 struct pair_dof_parent_owners_min_op {
385  using t = std::pair<size_t,size_t>; // pair: first=parent_owner & second=dof_owner
386  t operator() (const t& a, const t& b) {
387  // if "a" or "b" parent owner match the dof owner: select it
388  // otherwise, choose the parent owner arbitrarily
389  // avoid unset=max(size_t) by tacking the min
390  if (a.first == a.second) return a;
391  if (b.first == b.second) return b;
392  size_t parent_owner = std::min (a.first, b.first); // avoid unset=max
393  return std::make_pair (parent_owner, a.second);
394  }
395 };
396 } // namespace details
397 
398 template <class T, class M>
399 bool
400 space_base_rep<T,M>::parent_subgeo_owner_guard() const
401 {
402  if (_parent_subgeo_owner_initialized) {
403  return _parent_subgeo_owner_reattribution_required;
404  }
405  _parent_subgeo_owner_initialized = true;
406  freeze_guard();
407 
408  check_macro (! get_constitution().is_hierarchical(), "parent_subgeo_owner_reattribution: invalid product-space");
409  const space_constitution_terminal<T,M>& terminal_constit = get_constitution().get_terminal();
410  size_type nproc = ownership().comm().size();
411  size_type my_proc = ownership().comm().rank();
412  const geo_basic<T,M>& omega = terminal_constit.get_geo();
413  _parent_subgeo_owner_reattribution_required
414  = terminal_constit.get_basis().is_continuous()
415  && nproc > 1
416  && omega.variant() == geo_abstract_base_rep<T>::geo_domain;
417  if (! _parent_subgeo_owner_reattribution_required) {
418  return _parent_subgeo_owner_reattribution_required;
419  }
420  // try to reduce comms: when a parent owner is the dof owner, prefer it
421  // --> see pair_dof_parent_owners_min_op
422  std::vector<size_type> dis_idx;
423  size_type unset = std::numeric_limits<size_type>::max();
424  std::pair<size_type,size_type> unset_pair (unset, unset);
425  disarray<std::pair<size_type,size_type>> dof_parent_choice (ownership(), unset_pair);
426  size_type parent_proc = my_proc;
427  std::set<size_type> ext_dis_idof;
428  for (size_type ie = 0, ne = omega.size(); ie < ne; ++ie) {
429  const geo_element& K = omega [ie];
430  get_constitution().assembly_dis_idof (omega, K, dis_idx);
431  for (size_type loc_idof = 0, loc_ndof = dis_idx.size(); loc_idof < loc_ndof; loc_idof++) {
432  size_type dis_idof = dis_idx[loc_idof];
433  size_type dof_proc = ownership().find_owner (dis_idof);
434  dof_parent_choice.dis_entry (dis_idof) = std::make_pair (parent_proc, dof_proc);
435  if (dof_proc != my_proc) {
436  ext_dis_idof.insert (dis_idof);
437  }
438  }
439  }
440  dof_parent_choice.dis_entry_assembly (details::pair_dof_parent_owners_min_op());
441  dof_parent_choice.set_dis_indexes (ext_dis_idof);
442 
443  _parent_subgeo_owner.resize (ownership());
444  std::fill (_parent_subgeo_owner.begin(), _parent_subgeo_owner.end(), std::numeric_limits<size_type>::max());
445  for (size_type ie = 0, ne = omega.size(); ie < ne; ++ie) {
446  const geo_element& K = omega [ie];
447  get_constitution().assembly_dis_idof (omega, K, dis_idx);
448  for (size_type loc_idof = 0, loc_ndof = dis_idx.size(); loc_idof < loc_ndof; loc_idof++) {
449  size_type dis_idof = dis_idx[loc_idof];
450  size_type chosen_parent_proc = dof_parent_choice.dis_at(dis_idof).first;
451  // when chosen_parent_proc == unset : it means that this dof is orphean
452  // ie there are no element K in the geo_domain containing this dof and such that K belongs to the dof_owner
453  size_type reattributed_dof_proc = (chosen_parent_proc != unset) ? chosen_parent_proc : my_proc;
454  _parent_subgeo_owner.dis_entry (dis_idof) = reattributed_dof_proc;
455 #ifdef TO_CLEAN
456  size_type dof_proc = ownership().find_owner (dis_idof);
457  if (dof_proc != my_proc) {
458  warning_macro ("K="<<K.name()<<K.dis_ie()<<", dis_idof="<<dis_idof<<": dof_proc="<<dof_proc<<" != my_proc="<<my_proc<<": chosen_parent_proc="<<chosen_parent_proc<<", reattributed_dof_proc="<<reattributed_dof_proc);
459  }
460 #endif // TO_CLEAN
461  }
462  }
463  // choose for each dis_idof the proc owner that has the min index:
464  _parent_subgeo_owner.dis_entry_assembly (details::generic_set_min_op());
465  _parent_subgeo_owner.set_dis_indexes (ext_dis_idof);
466  return _parent_subgeo_owner_reattribution_required;
467 }
468 // ----------------------------------------------------------------------------
469 // instanciation in library
470 // ----------------------------------------------------------------------------
471 #define _RHEOLEF_instanciation(T,M) \
472 template class space_base_rep<T,M>; \
473 template space_basic<T,M> space_basic<T,M>::real();
474 
475 _RHEOLEF_instanciation(Float,sequential)
476 #ifdef _RHEOLEF_HAVE_MPI
477 _RHEOLEF_instanciation(Float,distributed)
478 #endif // _RHEOLEF_HAVE_MPI
479 #undef _RHEOLEF_instanciation
480 
481 } // namespace rheolef
field::size_type size_type
Definition: branch.cc:430
see the Float page for the full documentation
see the communicator page for the full documentation
space_pair_type::size_type size_type
Definition: space.h:164
#define trace_macro(message)
Definition: dis_macros.h:111
#define warning_macro(message)
Definition: dis_macros.h:53
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 dis_inod(const basis_basic< T > &b, const geo_size &gs, const geo_element &K, typename std::vector< size_type >::iterator dis_inod_tab)
size_type nnod(const basis_basic< T > &b, const geo_size &gs, size_type map_dim)
size_type dis_nnod(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.
t operator()(const t &a, const t &b)
Definition: space.cc:386
#define _RHEOLEF_instanciation(T, M)
Definition: space.cc:471
#define _RHEOLEF_space_real(M)
Definition: space.cc:310
Expr1::memory_type M
Definition: vec_expr_v2.h:416