rheolef  6.5
geo_mpi_get.cc
Go to the documentation of this file.
1 #include "rheolef/config.h"
2 #ifdef _RHEOLEF_HAVE_MPI
3 #include "rheolef/geo.h"
4 #include "rheolef/geo_domain.h"
5 #include "rheolef/dis_macros.h"
6 #include "rheolef/rheostream.h"
7 #include "rheolef/iorheo.h"
8 #include "rheolef/index_set.h"
9 
10 namespace rheolef {
11 
12 extern
13 array<size_t>
15  hack_array<geo_element_hack> ios_geo_element [reference_element::max_variant],
16  const distributor& ownership_by_dimension,
17  size_t map_dim,
18  size_t dis_nv);
19 
20 #ifdef TO_CLEAN
21 std::streambuf* strm_buffer = 0;
22 std::ofstream file;
23 void start_redirect (std::ostream& strm, std::string filename) {
24  file.open (filename.c_str());
25  strm_buffer = strm.rdbuf(); // save output buffer of the stream
26  strm.rdbuf (file.rdbuf()); // redirect ouput into the file
27 }
28 void end_redirect (std::ostream& strm, std::string="") {
29  strm.rdbuf (strm_buffer);
30 }
31 #endif // TO_CLEAN
32 
33 static
34 distributor
37  size_t side_dim)
38 {
39  using namespace std;
43  if (first_variant >= last_variant) return distributor();
44  const communicator& comm = ios_geo_element [first_variant].ownership().comm();
45 
46  size_type nge = 0;
47  size_type dis_nge = 0;
49  variant < reference_element:: last_variant_by_dimension(side_dim); variant++) {
50  nge += ios_geo_element [variant].size();
51  dis_nge += ios_geo_element [variant].dis_size();
52  }
53  return distributor (dis_nge, comm, nge);
54 }
55 static
56 void
59  size_t side_dim,
60  distributor& ios_ownership,
61  distributor ios_ownership_by_variant [reference_element::max_variant])
62 {
63  using namespace std;
67  if (first_variant >= last_variant) return;
68  const communicator& comm = ios_geo_element [first_variant].ownership().comm();
69 
70  size_type dis_nge = 0;
72  variant < reference_element:: last_variant_by_dimension(side_dim); variant++) {
73  dis_nge += ios_geo_element [variant].dis_size();
74  }
75  ios_ownership = distributor (dis_nge, comm, distributor::decide);
76 
77  size_type first_dis_v = 0;
79  variant < reference_element:: last_variant_by_dimension(side_dim); variant++) {
80 
81  size_type dis_ngev = ios_geo_element [variant].dis_size();
82  size_type last_dis_v = first_dis_v + dis_ngev;
83  size_type first_ios_dis_nge = ios_ownership.first_index();
84  size_type last_ios_dis_nge = ios_ownership. last_index();
85  size_type first_ios_dis_ngev = min (max (first_ios_dis_nge, first_dis_v), last_ios_dis_nge);
86  size_type last_ios_dis_ngev = max (min ( last_ios_dis_nge, last_dis_v), first_ios_dis_nge);
87  size_type ios_ngev = last_ios_dis_ngev - first_ios_dis_ngev;
88  ios_ownership_by_variant [variant] = distributor (dis_ngev, comm, ios_ngev);
89  first_dis_v = last_dis_v;
90  }
91 }
92 // => a side S belongs to the same proc as at least a super-side K that contains it
93 void
96  const geo_size& ios_gs, // TODO: not used...
97  size_t S_dim,
98  boost::array<std::vector<size_t>, 4>& massive_partition_by_dimension,
99  boost::array<array<size_t>, reference_element::max_variant>& partition_by_variant)
100 {
101  typedef size_t size_type;
102  distributor ios_vertex_ownership = ios_geo_element [reference_element::p].ownership();
103  size_type dis_nv = ios_vertex_ownership.dis_size();
104  size_type first_ios_dis_iv = ios_vertex_ownership.first_index();
105  index_set empty_set; // TODO: add a global allocator to empty_set
106  array<index_set,distributed> K_ball (ios_vertex_ownership, empty_set);
107  size_type K_dim = S_dim + 1;
108  size_type K_ios_ige = 0;
109  distributor true_K_ios_ge_ownership = build_true_ios_ge_ownership_by_dimension (ios_geo_element, K_dim);
110  size_type first_K_ios_dis_ige = true_K_ios_ge_ownership.first_index();
111  for (size_type K_variant = reference_element::first_variant_by_dimension(K_dim);
112  K_variant < reference_element:: last_variant_by_dimension(K_dim); K_variant++) {
113  distributor K_ios_gev_ownership = ios_geo_element [K_variant].ownership();
114  for (size_type K_ios_igev = 0, K_ios_ngev = K_ios_gev_ownership.size(); K_ios_igev < K_ios_ngev; K_ios_igev++, K_ios_ige++) {
115  const geo_element& K = ios_geo_element [K_variant] [K_ios_igev];
116  size_type K_ios_dis_ige = first_K_ios_dis_ige + K_ios_ige;
117  for (size_type iloc = 0, nloc = K.size(); iloc < nloc; iloc++) {
118  size_type ios_dis_iv = K[iloc];
119  assert_macro (ios_dis_iv < dis_nv, "K={"<<K<<"}: "<<iloc<<"-vertex index " << ios_dis_iv << " out of range [0:"<<dis_nv<<"[");
120  if (ios_vertex_ownership.is_owned(ios_dis_iv)) {
121  size_type ios_iv = ios_dis_iv - first_ios_dis_iv;
122  K_ball[ios_iv] += K_ios_dis_ige;
123  } else {
124  index_set K_ios_dis_ige_set;
125  K_ios_dis_ige_set += K_ios_dis_ige;
126  K_ball.dis_entry (ios_dis_iv) += K_ios_dis_ige_set; // not so efficient: union with {dis_iegv}
127  }
128  }
129  }
130  }
131  K_ball.dis_entry_assembly();
132  index_set ext_ios_dis_iv;
133  for (size_type S_variant = reference_element::first_variant_by_dimension(S_dim);
134  S_variant < reference_element:: last_variant_by_dimension(S_dim); S_variant++) {
135  distributor ios_gev_ownership = ios_geo_element [S_variant].ownership();
136  for (size_type ios_igev = 0, ios_ngev = ios_gev_ownership.size(); ios_igev < ios_ngev; ios_igev++) {
137  geo_element& S = ios_geo_element [S_variant] [ios_igev];
138  for (size_type iloc = 0, nloc = S.size(); iloc < nloc; iloc++) {
139  size_type ios_dis_iv = S[iloc];
140  if (! ios_vertex_ownership.is_owned (ios_dis_iv)) {
141  ext_ios_dis_iv += ios_dis_iv;
142  }
143  }
144  }
145  }
146  K_ball.set_dis_indexes (ext_ios_dis_iv);
147  distributor true_S_ios_ge_ownership = build_true_ios_ge_ownership_by_dimension (ios_geo_element, S_dim);
148  size_type S_dis_nge = true_S_ios_ge_ownership.dis_size();
149  size_type first_S_ios_dis_ige = true_S_ios_ge_ownership.first_index();
150  std::vector<size_type> tmp_S_massive_partition (S_dis_nge, 0); // TODO: massive memory area
151  size_type S_ios_ige = 0;
152  for (size_type S_variant = reference_element::first_variant_by_dimension(S_dim);
153  S_variant < reference_element:: last_variant_by_dimension(S_dim); S_variant++) {
154  distributor ios_gev_ownership = ios_geo_element [S_variant].ownership();
155  size_type first_ios_dis_igev = ios_gev_ownership.first_index();
156  for (size_type ios_igev = 0, ios_ngev = ios_gev_ownership.size(); ios_igev < ios_ngev; ios_igev++, S_ios_ige++) {
157  geo_element& S = ios_geo_element [S_variant] [ios_igev];
158  index_set K_ios_ige_set = K_ball.dis_at (S[0]);
159  for (size_type iloc = 0, nloc = S.size(); iloc < nloc; iloc++) {
160  K_ios_ige_set.inplace_intersection (K_ball.dis_at(S[iloc]));
161  }
162  check_macro (K_ios_ige_set.size() > 0, "connectivity: S={"<<S<<"} not found in the side set");
163  size_type S_owner = 0;
164  for (index_set::const_iterator iter = K_ios_ige_set.begin(), last = K_ios_ige_set.end(); iter != last; iter++) {
165  size_type K_ios_ige = *iter;
166  S_owner = std::max(S_owner, massive_partition_by_dimension[K_dim][K_ios_ige]);
167  }
168  size_type S_ios_dis_ige = first_S_ios_dis_ige + S_ios_ige;
169  tmp_S_massive_partition [S_ios_dis_ige] = S_owner;
170  }
171  }
172  massive_partition_by_dimension[S_dim].resize (tmp_S_massive_partition.size(), 0);
173  communicator comm = ios_vertex_ownership.comm();
174  mpi::all_reduce (
175  comm,
176  tmp_S_massive_partition.begin().operator->(),
177  tmp_S_massive_partition.size(),
178  massive_partition_by_dimension[S_dim].begin().operator->(),
179  mpi::maximum<size_type>());
180  size_type S_ios_dis_ige = first_S_ios_dis_ige;
181  for (size_type S_variant = reference_element::first_variant_by_dimension(S_dim);
182  S_variant < reference_element:: last_variant_by_dimension(S_dim); S_variant++) {
183  partition_by_variant [S_variant].resize (ios_geo_element [S_variant].ownership());
184  for (size_type S_ios_igev = 0, S_ios_negv = partition_by_variant [S_variant].size();
185  S_ios_igev < S_ios_negv; S_ios_igev++, S_ios_dis_ige++) {
186  partition_by_variant [S_variant][S_ios_igev] = massive_partition_by_dimension[S_dim] [S_ios_dis_ige];
187  }
188  }
189 }
190 void
193  const geo_size& ios_gs,
194  size_t dis_nv,
195  size_t side_dim,
196  hack_array<geo_element_hack> geo_element [reference_element::max_variant],
197  geo_size& gs,
198  array<size_t> igev2ios_dis_igev [reference_element::max_variant],
199  array<size_t> ios_igev2dis_igev [reference_element::max_variant],
200  array<size_t> ios_ige2dis_ige [4],
201  boost::array<array<size_t>, reference_element::max_variant>& partition_by_variant)
202 {
204  typedef hack_array<geo_element_hack>::iterator iterator_by_variant;
205  typedef hack_array<geo_element_hack>::const_iterator const_iterator_by_variant;
206 
207  communicator comm = ios_geo_element[0].ownership().comm();
208  size_type nge = 0;
209  for (size_type variant = reference_element::first_variant_by_dimension(side_dim);
210  variant < reference_element:: last_variant_by_dimension(side_dim); variant++) {
211 
212  ios_geo_element [variant].repartition (
213  partition_by_variant [variant],
214  geo_element [variant],
215  igev2ios_dis_igev [variant],
216  ios_igev2dis_igev [variant]);
217 
218  gs.ownership_by_variant [variant] = geo_element [variant].ownership();
219  nge += geo_element [variant].size();
220  }
221  gs.ownership_by_dimension [side_dim] = distributor (distributor::decide, comm, nge);
222  {
223  size_type first_v = 0;
224  for (size_type variant = reference_element::first_variant_by_dimension(side_dim);
225  variant < reference_element:: last_variant_by_dimension(side_dim); variant++) {
226 
227  gs.first_by_variant[variant] = distributor (distributor::decide, comm, first_v);
228  first_v += geo_element [variant].size();
229  }
230  }
231  ios_ige2dis_ige [side_dim].resize (ios_gs.ownership_by_dimension [side_dim]); // OK
232  {
233  size_type first_dis_ige = gs.ownership_by_dimension [side_dim].first_index();
234  size_type ige = 0;
235  for (size_type variant = reference_element::first_variant_by_dimension(side_dim);
236  variant < reference_element:: last_variant_by_dimension(side_dim); variant++) {
237  size_type first_v = gs.first_by_variant [variant].size();
238  for (size_type igev = 0, ngev = geo_element [variant].size(); igev < ngev; igev++, ige++) {
239  ::rheolef::geo_element& S = geo_element [variant] [igev]; // GNU C++ 4.5 BUG : add ::rheolef::
240  size_type dis_ige = first_dis_ige + ige;
241  S.set_dis_ie (dis_ige);
242  size_type ios_dis_ige = S.ios_dis_ie();
243  ios_ige2dis_ige [side_dim].dis_entry (ios_dis_ige) = dis_ige;
244  }
245  }
246  ios_ige2dis_ige [side_dim].dis_entry_assembly();
247  }
248 }
249 void
251  const std::vector<geo_element::size_type>& new_global_node_num,
252  size_t dis_nnod,
254 {
255  using namespace std;
257  size_type first_dis_igev = gev.ownership().first_index();
258  for (size_type igev = 0, ngev = gev.size(); igev < ngev; igev++) {
259  geo_element& S = gev [igev];
260  for (size_type iloc = 0, nloc = S.n_node(); iloc < nloc; iloc++) {
261  assert_macro (S[iloc] < dis_nnod, "node index "<<S[iloc]<<" out of range [0:"<<dis_nnod<<"[");
262  assert_macro (new_global_node_num[S[iloc]] < dis_nnod, "new node index "<<new_global_node_num[S[iloc]] <<" out of range [0:"<<dis_nnod<<"[");
263  S[iloc] = new_global_node_num [S[iloc]];
264  }
265  }
266 }
274 template <class T>
275 void
277 {
278  distributor vertex_ownership = base::_geo_element[reference_element::p].ownership();
279  size_type first_dis_iv = vertex_ownership.first_index();
280  size_type last_dis_iv = vertex_ownership.last_index();
281  size_type dis_nv = vertex_ownership.dis_size();
282 
283  distributor node_ownership = base::_node.ownership();
284  size_type first_dis_inod = node_ownership.first_index();
285  size_type last_dis_inod = node_ownership.last_index();
286  size_type dis_nnod = node_ownership.dis_size();
287 
288  std::set<size_type> ext_vertex_set;
289  std::set<size_type> ext_node_set;
290  std::vector<size_type> dis_inod1;
291  for (size_type dim = 1; dim <= base::_map_dimension; dim++) {
292  for (size_type ige = 0, nge = base::_gs.ownership_by_dimension[dim].size(); ige < nge; ige++) {
293  const geo_element& K = get_geo_element(dim,ige);
294  base::_numbering.dis_idof (base::_gs, K, dis_inod1);
295  for (size_type loc_inod = 0, loc_nnod = dis_inod1.size(); loc_inod < loc_nnod; loc_inod++) {
296  size_type dis_inod = dis_inod1 [loc_inod];
297  assert_macro (dis_inod < dis_nnod, "node index "<< dis_inod <<" out of range [0:"<<dis_nnod<<"[");
298  if (node_ownership.is_owned (dis_inod)) continue;
299  ext_node_set.insert (dis_inod);
300  if (loc_inod >= K.size()) continue;
301  size_type dis_iv = base::dis_inod2dis_iv (dis_inod);
302  check_macro (dis_iv < dis_nv, "vertex index "<< dis_iv <<" out of range [0:"<<dis_nv<<"[");
303  check_macro (!vertex_ownership.is_owned (dis_iv), "strange bug: not owned (nod) but is_owned(ver)");
304  ext_vertex_set.insert (dis_iv);
305  }
306  }
307  }
308  base::_node.set_dis_indexes (ext_node_set);
309  base::_geo_element [reference_element::p].set_dis_indexes (ext_vertex_set);
310 }
315 template <class T>
316 void
318 {
319  distributor vertex_ownership = base::_gs.ownership_by_dimension [0];
320  distributor node_ownership = base::_node.ownership();
321  size_type nv = vertex_ownership.size();
322  size_type first_dis_iv = vertex_ownership.first_index();
323  size_type last_dis_iv = vertex_ownership. last_index();
324  size_type dis_nnod = node_ownership.dis_size(); // for bound checks
325  index_set ext_iv_set; // size=O((N/nproc)^((d-1)/d))
326  for (size_type dim = side_dim; dim <= base::_map_dimension; dim++) {
329  for (size_type igev = 0, ngev = base::_geo_element[var].size(); igev < ngev; igev++) {
330  const geo_element& K = base::_geo_element [var] [igev];
331  for (size_type iloc = 0, nloc = K.size(); iloc < nloc; iloc++) {
332  size_type dis_inod = K[iloc];
333  if (! node_ownership.is_owned(dis_inod)) { // vertex ownership follows node ownership
334  size_type dis_iv = base::dis_inod2dis_iv (dis_inod);
335  ext_iv_set.insert (dis_iv);
336  }
337  }
338  }
339  }
340  }
341  index_set empty_set; // TODO: add a global allocator to empty_set
342  array<index_set,distributed> ball [reference_element::max_variant];
344  variant < reference_element:: last_variant_by_dimension(side_dim); variant++) {
345  ball [variant].resize (vertex_ownership, empty_set);
346  distributor gev_ownership = base::_gs.ownership_by_variant [variant];
347  size_type first_dis_igev = gev_ownership.first_index();
348  for (size_type igev = 0, ngev = base::_geo_element [variant].size(); igev < ngev; igev++) {
349  const geo_element& S = base::_geo_element [variant] [igev];
350  size_type dis_igev = first_dis_igev + igev;
351  for (size_type iloc = 0, nloc = S.size(); iloc < nloc; iloc++) {
352  size_type dis_inod = S[iloc];
353  assert_macro (dis_inod < dis_nnod, "S={"<<S<<"}: "<<iloc<<"-node index " << dis_inod << " out of range [0:"<<dis_nnod<<"[");
354  size_type dis_iv = base::dis_inod2dis_iv (dis_inod);
355  if (vertex_ownership.is_owned(dis_iv)) {
356  size_type iv = dis_iv - first_dis_iv;
357  ball [variant][iv] += dis_igev;
358  } else {
359  index_set dis_igev_set;
360  dis_igev_set += dis_igev;
361  ball [variant].dis_entry (dis_iv) += dis_igev_set; // not so efficient: union with {dis_iegv}
362  }
363  }
364  }
365  ball [variant].dis_entry_assembly();
366  // -> ceux issus des triangles, tetra, ect, et externes a la partition ?
367  ball [variant].set_dis_indexes (ext_iv_set);
368  }
370  variant < reference_element:: last_variant_by_dimension(side_dim); variant++) {
371 
372  std::set<size_type> ext_igev_set; // size=O((N/nproc)^((d-1)/d))
373  distributor gev_ownership = base::_gs.ownership_by_variant [variant];
374  size_type dis_ngev = gev_ownership.dis_size();
375  for (size_type iv = 0, nv = vertex_ownership.size(); iv < nv; iv++) {
376 
377  const index_set& ball_iv = ball [variant] [iv];
378  for (index_set::const_iterator iter = ball_iv.begin(), last = ball_iv.end(); iter != last; iter++) {
379  size_type dis_igev = *iter;
380  assert_macro (dis_igev < dis_ngev, "index "<<dis_igev<<" of element variant="<<variant<<" is out of range [0:"<<dis_ngev<<"[");
381  if (! gev_ownership.is_owned (dis_igev)) {
382  ext_igev_set.insert (dis_igev);
383  }
384  }
385  }
386  typedef array<index_set>::scatter_map_type map_type;
387  const map_type& ball_dis = ball [variant].get_dis_map_entries();
388  for (typename map_type::const_iterator iter_b = ball_dis.begin(), last_b = ball_dis.end(); iter_b != last_b; iter_b++) {
389  size_type dis_iv = (*iter_b).first;
390  const index_set& ball_dis_iv = (*iter_b).second;
391  for (index_set::const_iterator iter = ball_dis_iv.begin(), last = ball_dis_iv.end(); iter != last; iter++) {
392  size_type dis_igev = *iter;
393  assert_macro (dis_igev < dis_ngev, "index "<<dis_igev<<" of element variant="<<variant<<" is out of range [0:"<<dis_ngev<<"[");
394  if (! gev_ownership.is_owned (dis_igev)) {
395  ext_igev_set.insert (dis_igev);
396  }
397  }
398  }
399  base::_geo_element[variant].set_dis_indexes (ext_igev_set);
400  }
401  // => on numerote dis_iedg cette arete dans le geo_element K
402  for (size_type dim = side_dim+1; dim <= base::_map_dimension; dim++) {
405  for (size_type igev = 0, ngev = base::_geo_element [var].size(); igev < ngev; igev++) {
406  geo_element& K = base::_geo_element [var] [igev];
407  size_type S_iv [4];
408  size_type S_inod[4];
409  for (size_type loc_isid = 0, loc_nsid = K.n_subgeo(side_dim); loc_isid < loc_nsid; loc_isid++) {
410  size_type S_size = K.subgeo_size (side_dim, loc_isid);
411  size_type S_variant = reference_element::variant (S_size, side_dim);
412  size_type loc_jv0 = K.subgeo_local_vertex (side_dim, loc_isid, 0);
413  S_inod [0] = K[loc_jv0];
414  S_iv [0] = base::dis_inod2dis_iv (S_inod[0]);
415  const array<index_set,distributed>& ball_S_variant = ball [S_variant];
416  index_set igev_set = ball_S_variant.dis_at (S_iv[0]);
417  for (size_type sid_jloc = 1, sid_nloc = S_size; sid_jloc < sid_nloc; sid_jloc++) {
418  size_type loc_jv = K.subgeo_local_vertex (side_dim, loc_isid, sid_jloc);
419  S_inod [sid_jloc] = K[loc_jv];
420  S_iv [sid_jloc] = base::dis_inod2dis_iv (S_inod[sid_jloc]);
421  const index_set& ball_jv = ball_S_variant.dis_at (S_iv[sid_jloc]);
422  igev_set.inplace_intersection (ball_jv);
423  }
424  check_macro (igev_set.size() > 0, "connectivity: "<<S_size<<"-side ("
425  <<S_iv[0]<<","<<S_iv[1]<<","<<S_iv[2]<<","<<S_iv[3]<<") not found in the side set");
426  check_macro (igev_set.size() == 1, "connectivity: the same side is multiply represented");
427  size_type dis_igev = *(igev_set.begin());
428  const geo_element& S = base::_geo_element[S_variant].dis_at(dis_igev);
429  size_type dis_isid = S.dis_ie();
430  if (side_dim == 1) {
431  geo_element::orientation_type orient = S.get_edge_orientation (S_inod[0], S_inod[1]);
432  K.edge_indirect (loc_isid).set (orient, dis_isid);
433  } else { // side_dim == 2
436  if (K.subgeo_size (side_dim, loc_isid) == 3) {
437  S.get_orientation_and_shift (S_inod[0], S_inod[1], S_inod[2], orient, shift);
438  } else {
439  S.get_orientation_and_shift (S_inod[0], S_inod[1], S_inod[2], S_inod[3], orient, shift);
440  }
441  K.face_indirect (loc_isid).set (orient, dis_isid, shift);
442  }
443  }
444  }
445  }
446  }
447 }
448 template <class T>
449 void
451  boost::array<size_type,reference_element::max_variant>& local_ndof_by_variant,
452  array<size_type,distributed>& idof2ios_dis_idof) const
453 {
454  size_type ndof = 0;
455  for (size_type variant = 0;
456  variant < reference_element::last_variant_by_dimension(base::map_dimension());
457  variant++) {
458  ndof += base::_gs.ownership_by_variant [variant].size() * local_ndof_by_variant [variant];
459  }
460  distributor dof_ownership (distributor::decide, base::comm(), ndof);
461  idof2ios_dis_idof.resize (dof_ownership);
462  size_type first_dis_idof= dof_ownership.first_index();
463  size_type idof = 0;
464  size_type first_ios_dis_v = 0;
465  for (size_type variant = 0;
466  variant < reference_element::last_variant_by_dimension(base::map_dimension());
467  variant++) {
468  size_type loc_ndof = local_ndof_by_variant [variant];
469  for (size_type igev = 0, ngev = base::_gs.ownership_by_variant [variant].size(); igev < ngev; igev++) {
470  size_type ios_dis_igev = _igev2ios_dis_igev [variant] [igev];
471  for (size_type loc_idof = 0; loc_idof < loc_ndof; loc_idof++, idof++) {
472  size_type ios_dis_idof = first_ios_dis_v + ios_dis_igev*loc_ndof + loc_idof;
473  idof2ios_dis_idof [idof] = ios_dis_idof;
474  }
475  }
476  first_ios_dis_v += base::_gs.ownership_by_variant [variant].dis_size() * loc_ndof;
477  }
478 }
479 template <class T>
480 void
482 {
483  boost::array<size_type,reference_element::max_variant> loc_ndof_by_variant;
484  reference_element::init_local_nnode_by_variant (base::order(), loc_ndof_by_variant);
485  set_ios_permutation (loc_ndof_by_variant, _inod2ios_dis_inod);
486 }
487 template <class T>
490 {
491  using namespace std;
492  check_macro (ips.good(), "bad input stream for geo.");
493  communicator comm = base::_geo_element[reference_element::p].ownership().comm();
494 
496  size_type my_proc = comm.rank();
497  check_macro (dis_scatch(ips, ips.comm(), "\nmesh"), "input stream does not contains a geo.");
498  ips >> base::_version;
499  check_macro (base::_version == 4, "geo version < 4 not supported (HINT: see geo -upgrade)");
500 
501  geo_header hdr;
502  ips >> hdr;
503  check_macro (! hdr.need_upgrade(),
504  "unsupported geo without connectivity in the distributed version: HINT: use geo_upgrade");
505  base::_have_connectivity = true;
506  base::_dimension = hdr.dimension;
507  base::_map_dimension = hdr.map_dimension;
508  base::_sys_coord = hdr.sys_coord;
509  base::_name = "unnamed";
510  base::_numbering.set_degree (hdr.order);
511 
512  size_type dis_nnod = hdr.dis_size_by_dimension [0];
513  size_type dis_nedg = hdr.dis_size_by_dimension [1];
514  size_type dis_nfac = hdr.dis_size_by_dimension [2];
515  size_type dis_ne = hdr.dis_size_by_dimension [base::_map_dimension];
516  size_type ios_size_by_variant [reference_element::max_variant];
517  std::fill (ios_size_by_variant, ios_size_by_variant+reference_element::max_variant, 0);
518  array<node_type> ios_node (dis_nnod);
519  ios_node.get_values (ips, _point_get<T>(base::_dimension));
520  check_macro (ips.good(), "bad input stream for geo.");
522  size_type ios_nv = 0;
523  if (base::_map_dimension == 0) {
524  ios_nv = ios_node.size();
525  } else {
526  // set elt ios index and compute map_dim & ios_nv (ios_nv < ios_nnod when order > 1)
527  size_type ios_ne = 0;
528  size_type dis_ne = 0;
529  for (size_type variant = reference_element::first_variant_by_dimension(base::_map_dimension);
530  variant < reference_element:: last_variant_by_dimension(base::_map_dimension); variant++) {
531  geo_element::parameter_type param (variant, 1);
532  ios_geo_element [variant].resize (hdr.dis_size_by_variant [variant], param);
533  ios_geo_element [variant].get_values (ips);
534  _ios_gs.first_by_variant [variant] = distributor (distributor::decide, base::comm(), ios_ne);
535  ios_ne += ios_geo_element [variant].size();
536  dis_ne += ios_geo_element [variant].dis_size();
537  }
539  ios_geo_element,
540  base::_map_dimension,
541  _ios_gs.ownership_by_dimension [base::_map_dimension],
542  _ios_gs.ownership_by_variant);
543 
544  vector<size_type> local_is_vertex (dis_nnod, 0);
545  for (size_type variant = reference_element::first_variant_by_dimension(base::_map_dimension);
546  variant < reference_element:: last_variant_by_dimension(base::_map_dimension); variant++) {
547  size_type first_ios_dis_v = _ios_gs.first_by_variant [variant].dis_size();
548  size_type first_ios_dis_igev = ios_geo_element [variant].ownership().first_index();
549  size_type ios_igev = 0;
550  for (iterator_by_variant iter = ios_geo_element [variant].begin(), last = ios_geo_element [variant].end();
551  iter != last; iter++, ios_igev++) {
552  geo_element& K = *iter;
553  size_type ios_dis_ie = first_ios_dis_v + first_ios_dis_igev + ios_igev;
554  K.set_ios_dis_ie (ios_dis_ie); // OK
555  ios_size_by_variant [K.variant()]++;
556  if (base::order() > 1) {
557  for (size_type iloc = 0, nloc = K.size(); iloc < nloc; iloc++) {
558  local_is_vertex [K[iloc]] = 1;
559  }
560  }
561  }
562  }
563  if (base::order() == 1) {
564  ios_nv = ios_node.size();
565  } else {
566  vector<size_type> global_is_vertex (dis_nnod, 0);
567  mpi::all_reduce (
568  comm,
569  local_is_vertex.begin().operator->(),
570  local_is_vertex.size(),
571  global_is_vertex.begin().operator->(),
572  mpi::maximum<size_type>());
573  ios_nv = accumulate (global_is_vertex.begin() + ios_node.ownership().first_index(),
574  global_is_vertex.begin() + ios_node.ownership().last_index(), 0);
575  }
576  }
577  distributor ios_vertex_ownership;
578  if (base::order() == 1) {
579  ios_vertex_ownership = ios_node.ownership();
580  } else {
581  ios_vertex_ownership = distributor (distributor::decide, base::comm(), ios_nv); // not sure...
582  }
583  size_type dis_nv = ios_vertex_ownership.dis_size();
584  {
586  ios_geo_element [reference_element::p].resize (ios_vertex_ownership, param);
587  _ios_gs.ownership_by_variant [reference_element::p] = ios_vertex_ownership;
588  _ios_gs.ownership_by_dimension [0] = ios_vertex_ownership;
589  size_type first_ios_dis_iv = ios_vertex_ownership.first_index();
590  for (size_type ios_iv = 0, ios_nv = ios_vertex_ownership.size(); ios_iv < ios_nv; ios_iv++) {
591  geo_element& P = ios_geo_element [reference_element::p] [ios_iv];
592  size_type ios_dis_iv = first_ios_dis_iv + ios_iv;
593  P [0] = ios_dis_iv;
594  P.set_ios_dis_ie (ios_dis_iv);
595  ios_size_by_variant [P.variant()]++;
596  }
597  }
598  if (base::_map_dimension > 0) {
599  for (size_type side_dim = base::_map_dimension - 1; side_dim >= 1; side_dim--) {
600  size_type ios_ngev = 0;
601  size_type dis_ngev = 0;
603  variant < reference_element:: last_variant_by_dimension(side_dim); variant++) {
604  geo_element::parameter_type param (variant, 1);
605  ios_geo_element [variant].resize (hdr.dis_size_by_variant [variant], param);
606  ios_geo_element [variant].get_values (ips);
607  _ios_gs.first_by_variant [variant] = distributor (distributor::decide, base::comm(), ios_ngev);
608  ios_ngev += ios_geo_element [variant].size();
609  dis_ngev += ios_geo_element [variant].dis_size();
610  }
612  ios_geo_element,
613  side_dim,
614  _ios_gs.ownership_by_dimension [side_dim],
615  _ios_gs.ownership_by_variant);
616 
618  variant < reference_element:: last_variant_by_dimension(side_dim); variant++) {
619  size_type ios_igev = 0;
620  size_type first_dis_v = _ios_gs.first_by_variant [variant].dis_size();
621  size_type first_ios_dis_igev = ios_geo_element [variant].ownership().first_index();
622  for (iterator_by_variant iter = ios_geo_element [variant].begin(), last = ios_geo_element [variant].end();
623  iter != last; iter++, ios_igev++) {
624  geo_element& S = *iter;
625  size_type ios_dis_igev = first_dis_v + first_ios_dis_igev + ios_igev;
626  S.set_ios_dis_ie (ios_dis_igev); // OK
627  ios_size_by_variant [S.variant()]++;
628  }
629  }
630  }
631  }
632  distributor true_ios_ownership_by_dimension [4];
633  true_ios_ownership_by_dimension [base::_map_dimension]
634  = build_true_ios_ge_ownership_by_dimension (ios_geo_element, base::_map_dimension);
635 
636  array<size_type> partition = geo_mpi_partition (
637  ios_geo_element,
638  true_ios_ownership_by_dimension [base::_map_dimension],
639  base::_map_dimension,
640  dis_nv);
641 
642  boost::array<array<size_t>, reference_element::max_variant> partition_by_variant;
643  {
644  typename array<size_type>::const_iterator iter = partition.begin();
645  for (size_type variant = reference_element::first_variant_by_dimension(base::_map_dimension);
646  variant < reference_element:: last_variant_by_dimension(base::_map_dimension); variant++) {
647  partition_by_variant [variant].resize (ios_geo_element [variant].ownership());
648  for (typename array<size_type>::iterator iter_by_var = partition_by_variant [variant].begin(),
649  last_by_var = partition_by_variant [variant].end();
650  iter_by_var != last_by_var; iter_by_var++, iter++) {
651  *iter_by_var = *iter;
652  }
653  }
654  }
656  ios_geo_element,
657  _ios_gs,
658  dis_nv,
659  base::_map_dimension,
660  base::_geo_element,
661  base::_gs,
662  _igev2ios_dis_igev,
663  _ios_igev2dis_igev,
664  _ios_ige2dis_ige,
665  partition_by_variant);
666 
667 
668  boost::array<std::vector<size_t>, 4> massive_partition_by_dimension;
669  {
670  std::vector<size_t> tmp_massive_partition (dis_ne, 0);
671  size_type true_first_dis_ige = true_ios_ownership_by_dimension [base::_map_dimension].first_index();
672  for (size_type ios_ie = 0, ios_ne = partition.size(); ios_ie < ios_ne; ios_ie++) {
673  size_type ios_dis_ie = true_first_dis_ige + ios_ie;
674  tmp_massive_partition [ios_dis_ie] = partition [ios_ie];
675  }
676  massive_partition_by_dimension [base::_map_dimension].resize (dis_ne);
677  mpi::all_reduce (
678  comm,
679  tmp_massive_partition.begin().operator->(),
680  tmp_massive_partition.size(),
681  massive_partition_by_dimension[base::_map_dimension].begin().operator->(),
682  mpi::maximum<size_type>());
683  }
684  for (size_type supergeo_dim = base::_map_dimension; supergeo_dim > 0; supergeo_dim--) {
686  ios_geo_element,
687  _ios_gs,
688  supergeo_dim-1,
689  massive_partition_by_dimension,
690  partition_by_variant);
691  }
692  partition_by_variant [reference_element::p].resize (ios_vertex_ownership);
693  size_type first_ios_dis_iv = ios_vertex_ownership.first_index();
694  for (size_type ios_iv = 0, ios_nv = ios_vertex_ownership.size(); ios_iv < ios_nv; ios_iv++) {
695  size_type ios_dis_iv = first_ios_dis_iv + ios_iv;
696  partition_by_variant [reference_element::p] [ios_iv] = massive_partition_by_dimension[0] [ios_dis_iv];
697  }
698  array<size_type> iv2ios_dis_iv;
699  ios_geo_element [reference_element::p].repartition (
700  partition_by_variant [reference_element::p],
701  base::_geo_element [reference_element::p],
702  iv2ios_dis_iv,
703  _ios_ige2dis_ige[0]);
704 
705  distributor vertex_ownership = base::_geo_element[reference_element::p].ownership();
706  base::_gs.ownership_by_dimension [0] = vertex_ownership;
707  base::_gs.ownership_by_variant [reference_element::p] = vertex_ownership;
708  base::_gs.first_by_variant [reference_element::p] = distributor (0, base::comm(), 0);
709  size_type first_dis_iv = vertex_ownership.first_index();
710  size_type last_dis_iv = vertex_ownership.last_index();
711  _igev2ios_dis_igev [reference_element::p].resize (vertex_ownership);
712  for (size_type iv = 0, nv = base::_geo_element[reference_element::p].size(); iv < nv; iv++) {
713  geo_element& P = base::_geo_element[reference_element::p] [iv];
714  size_type dis_iv = first_dis_iv + iv;
715  P[0] = dis_iv;
716  P.set_dis_ie (dis_iv);
717  _igev2ios_dis_igev [reference_element::p] [iv] = P.ios_dis_ie();
718  }
719  if (base::_map_dimension > 0) {
720  for (size_type side_dim = base::_map_dimension-1; side_dim > 0; side_dim--) {
721 
723  ios_geo_element,
724  _ios_gs,
725  dis_nv,
726  side_dim,
727  base::_geo_element,
728  base::_gs,
729  _igev2ios_dis_igev,
730  _ios_igev2dis_igev,
731  _ios_ige2dis_ige,
732  partition_by_variant);
733  }
734  }
735  do {
737  bool status = dom.get (ips, *this);
738  if (!status) break;
739  base::_domains.push_back (dom);
740  } while (true);
741  node_renumbering (ios_node.ownership());
742  distributor node_ownership = _inod2ios_dis_inod.ownership();
743  base::_gs.node_ownership = node_ownership;
744  _ios_inod2dis_inod.resize (ios_node.ownership(), std::numeric_limits<size_type>::max());
745  _inod2ios_dis_inod.reverse_permutation (_ios_inod2dis_inod);
746 
747  base::_node.resize (node_ownership);
748  ios_node.permutation_apply (_ios_inod2dis_inod, base::_node);
749  vector<size_type> new_local_node_num (dis_nnod, 0);
750  size_type first_ios_inod = _ios_inod2dis_inod.ownership().first_index();
751  size_type last_ios_inod = _ios_inod2dis_inod.ownership(). last_index();
752  for (size_type dis_ios_inod = first_ios_inod; dis_ios_inod < last_ios_inod; dis_ios_inod++) {
753  size_type ios_inod = dis_ios_inod - first_ios_inod;
754  new_local_node_num [dis_ios_inod] = _ios_inod2dis_inod[ios_inod];
755  }
756  vector<size_type> new_global_node_num (dis_nnod, 0);
757  mpi::all_reduce (
758  comm,
759  new_local_node_num.begin().operator->(),
760  new_local_node_num.size(),
761  new_global_node_num.begin().operator->(),
762  mpi::maximum<size_type>());
763  for (size_type dim = base::_map_dimension; dim > 0; dim--) {
765  variant < reference_element:: last_variant_by_dimension(dim); variant++) {
767  new_global_node_num,
768  dis_nnod,
769  base::_geo_element [variant]);
770  }
771  }
772  for (size_type dim = base::_map_dimension - 1; base::_map_dimension > 0 && dim > 0; dim--) {
773  set_element_side_index (dim);
774  }
775  build_external_entities ();
777  return ips;
778 }
779 template <class T>
780 void
782  std::string filename,
783  const communicator& comm)
784 {
785  idiststream ips;
786  ips.open (filename, "geo", comm);
787  check_macro(ips.good(), "\"" << filename << "[.geo[.gz]]\" not found.");
788  get (ips);
789  std::string root_name = delete_suffix (delete_suffix(filename, "gz"), "geo");
790  std::string name = get_basename (root_name);
791  base::_name = name;
792 }
793 template class geo_rep<Float,distributed>;
794 
795 } // namespace rheolef
796 #endif // _RHEOLEF_HAVE_MPI
797