~ubuntu-branches/ubuntu/saucy/rheolef/saucy

« back to all changes in this revision

Viewing changes to nfem/plib/domain_indirect_mpi.cc

  • Committer: Bazaar Package Importer
  • Author(s): Pierre Saramito
  • Date: 2011-03-23 11:14:26 UTC
  • mfrom: (1.1.4 upstream)
  • Revision ID: james.westby@ubuntu.com-20110323111426-cjvhey7lxt6077ty
Tags: 5.93-1
* New upstream release (minor changes):
  - some extra warning message deleted in heap_allocator
  - graphic output with mayavi2 fixed
  - add doc refman in .info and .pdf format

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
///
 
2
/// This file is part of Rheolef.
 
3
///
 
4
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
 
5
///
 
6
/// Rheolef is free software; you can redistribute it and/or modify
 
7
/// it under the terms of the GNU General Public License as published by
 
8
/// the Free Software Foundation; either version 2 of the License, or
 
9
/// (at your option) any later version.
 
10
///
 
11
/// Rheolef is sequential in the hope that it will be useful,
 
12
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
 
13
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 
14
/// GNU General Public License for more details.
 
15
///
 
16
/// You should have received a copy of the GNU General Public License
 
17
/// along with Rheolef; if not, write to the Free Software
 
18
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
 
19
/// 
 
20
/// =========================================================================
 
21
 
 
22
#include "rheolef/config.h"
 
23
 
 
24
#ifdef _RHEOLEF_HAVE_MPI
 
25
#include "rheolef/domain_indirect.h"
 
26
#include "rheolef/geo.h"
 
27
#include "rheolef/geo_connectivity.h"
 
28
 
 
29
namespace rheolef {
 
30
 
 
31
// ----------------------------------------------------------------------------
 
32
/// @brief domain_indirect_rep<mpi> i/o in format version 2
 
33
// ----------------------------------------------------------------------------
 
34
/**
 
35
 Implementation notes:
 
36
 
 
37
 The domain is implemented as an array of elements, e.g. edges or faces,
 
38
 but it could also be vertices or volume elements, for volumic domains.
 
39
 Says a domain of "geo_element" for simplicity. 
 
40
 The orientation is important for boundaries or interfaces
 
41
 and the direction of the normal has a sense in domains.
 
42
 It is described as an array of a pair: an orientation
 
43
 and an index of geo_element, denoted in brief as "oige".
 
44
 The index part of the pair, in distributed environment, is denoted
 
45
 as "dis_ige".
 
46
 An index that navigates in this array of pair is denoted
 
47
 as ioige, a short name for an index of (a pair of) orientation and
 
48
 index of a geo_element.
 
49
 
 
50
 1) Initial numbering ("ini")
 
51
 
 
52
 The array of oige pairs that constitutes the domain contains "dis_noige" entries:
 
53
 this array is denoted as "ini_oige".
 
54
 Its distributor is denoted as ini_ioige_ownership,
 
55
 and constructed simply from dis_noige and the communicator obtained from
 
56
 the geometry omega.
 
57
        ini_dis_ige = ini_oige[ini_ioige].index
 
58
 
 
59
 2) First renumbering ("ios")
 
60
 
 
61
 Indexes of geo_element (the dis_ige part of the pair oige), that are read from the
 
62
 input file, may be changed to a new numbering:
 
63
 the geo_elements has been renumbered after the partition of the geometry.
 
64
 The new geo_element numbering is also used for computations, see e.g. the "space" class.
 
65
 The initial numbering is refered in the geometry omega as,
 
66
 e.g. ios_dis_iedg or ios_dis_ifac, says here ios_dis_ige.
 
67
 The new side numbering is obtained by:
 
68
        dis_ige = omega.ios_ige2dis_ige (ios_ige)
 
69
 There is a first drawback in distributed environment:
 
70
 this call is valid when ios_ige is local to the current processor,
 
71
 while only the global ios_dis_ige index is valid. The ios_ige are
 
72
 distributed according to ios_ige_ownership (e.g. ios_edge_ownership
 
73
 or ios_face_ownership), as specified in the geometry omega.
 
74
 Thus, the array of osig pairs may be fisrt re-distributed from
 
75
 ini_ioige_ownership to another distributor, denoted as ios_ioige_ownership,
 
76
 so that the index part of the pair matches the geo_element ios_ige_ownership
 
77
 distribution. After redistribution of the array of pairs, these pairs
 
78
 have been reordered, i.e. renumbered. The permutation and inverse
 
79
 permutation arrays are denoted by:
 
80
        ini_ioige2ios_dis_ioige,        according to ini_ioige_ownership
 
81
        ios_ioige2ini_dis_ioige,        according to ios_ioige_ownership
 
82
 Also, the initial "ini_oige" array elements are reordered and re-distributed
 
83
 in a new array denoted as "ios_oige". The index of geo_element part of an
 
84
 element writes:
 
85
        ios_dis_ige = ios_oige [ios_ioige].index
 
86
 It can be converted into the
 
87
 local geo_element numbering as "ios_ige" and is then suitable for the obtention
 
88
 of the new geo_element numbering used for computations
 
89
        dis_ige = omega.ios_ige2dis_ige (ios_ige)
 
90
 
 
91
 3) Second renumbering (no prefix)
 
92
 
 
93
 The dis_ige index obtained by this way refers to a geo_element that is not owned in
 
94
 general by the current processor.
 
95
 So, the array of pairs may be a second time redistributed
 
96
 accordingly to the renumbered geo_element, i.e. the "ige_ownership" distributor
 
97
 (e.g. edge_ownership or face_ownership), as specified by the geometry omega.
 
98
 
 
99
 An array of pairs, containing the orientation and the new geo_element index "dis_ige"
 
100
 is created and denoted as "tmp_oige" :
 
101
        tmp_oige[ios_ioige] = pair(orient, dis_ige)
 
102
 Notice that this array is based on the ios_ige_ownership distributor.
 
103
 A new distribution follows the distribution of geo_elements and is denoted as "oige":
 
104
        oige[ioige] = pair(orient,dis_ige)
 
105
 such that the locally owned part of oige contains indexes of geo_elements
 
106
 that are also locally owned by the geometry.
 
107
 The permutation and inverse permutation arrays are denoted by:
 
108
        ios_ioige2dis_ioige,    according to ios_ioige_ownership
 
109
        ioige2ios_dis_ioige,    according to     ioige_ownership
 
110
 
 
111
 The "oige" array is stored directly in the class domain: domain is 
 
112
 implemented as a derived class of the array class which represents "oige".
 
113
 Thus "oige" is unamed in the domain class.
 
114
 
 
115
 4) Back to "ini" renumbering
 
116
 
 
117
 In order to write domains into files with the initial ordering and indexes
 
118
 of geo_element, we also build the direct permutation and inverse permutation arrays:
 
119
 
 
120
 for (ios_ioige=0..) {
 
121
    ini_dis_ioige = ios_ioige2ini_dis_ioige [ios_ioige]
 
122
    dis_ioige     = ios_ioige2dis_ioige     [ios_ioige]
 
123
    ini_ioige2dis_ioige.dis_entry (ini_dis_ioige) = dis_ioige
 
124
    ioige2ini_dis_ioige.dis_entry (dis_ioige      = ini_dis_ioige
 
125
 }
 
126
 ioige2ini_dis_ioige.dis_entry_assembly()
 
127
 ini_ioige2dis_ioige.dis_entry_assembly()
 
128
 
 
129
 5) Output domains
 
130
 
 
131
 for (ioige=0...)
 
132
    ini_dis_ioige = ioige2ini_dis_ioige [ioige]
 
133
    orient = oige [ioige].orient
 
134
    ige    = oige [ioige].index
 
135
    ios_dis_ige = omega.ige2ios_dis_ige (ige)
 
136
    ini_oige.dis_entry (ini_dis_ioige) = pair(orient, ios_dis_ige)
 
137
 }
 
138
 ini_oige.dis_entry_assembly
 
139
 ini_oige.put
 
140
 
 
141
*/
 
142
template<class U>
 
143
idiststream&
 
144
domain_indirect_rep<distributed>::get (
 
145
    idiststream&                   ips,
 
146
    const geo_rep<U,distributed>& omega)
 
147
{
 
148
  const size_type unset = std::numeric_limits<size_type>::max();
 
149
  // ---------------------
 
150
  // 1) get initial array
 
151
  // ---------------------
 
152
  // 1.1) get header
 
153
  communicator comm = omega.comm();
 
154
  if ( ! dis_scatch (ips, comm, "\ndomain")) {
 
155
    ips.is().setstate (std::ios::badbit);
 
156
    return ips;
 
157
  }
 
158
  size_type version, dis_noige;
 
159
  ips >> base::_name
 
160
      >> version >> base::_map_dim >> dis_noige;
 
161
  check_macro (version == 2, "unsupported version="<<version<<" domain format");
 
162
  std::fill (_ini_size_by_variant, _ini_size_by_variant+reference_element::max_size, 0);
 
163
 
 
164
  // 1.2) get data
 
165
  distributor ini_ioige_ownership (dis_noige, comm);
 
166
  array<domain_pair_type,distributed> ini_oige (ini_ioige_ownership);
 
167
  ini_oige.get_values (ips);
 
168
  // ---------------------------
 
169
  // 2) first renumbering (ios)
 
170
  // ---------------------------
 
171
  // 2.1) compute ios_owner for each oriented-side
 
172
  distributor ios_ige_ownership = omega.geo_element_ios_ownership (base::_map_dim);
 
173
  array<size_type>  ios_owner (ini_ioige_ownership, unset);
 
174
  for (size_type ini_ioige = 0, ini_noige = ini_oige.size();
 
175
        ini_ioige < ini_noige; ini_ioige++) {
 
176
    size_type ios_ige = ini_oige [ini_ioige].index();
 
177
    ios_owner [ini_ioige] = ios_ige_ownership.find_owner (ios_ige);
 
178
  }
 
179
  // 2.2) ios repartition
 
180
  array<domain_pair_type>  ios_oige;
 
181
  array<size_type>  ini_ioige2ios_dis_ioige;
 
182
  array<size_type>  ios_ioige2ini_dis_ioige;
 
183
  ini_oige.repartition (
 
184
        ios_owner,
 
185
        ios_oige,
 
186
        ios_ioige2ini_dis_ioige,
 
187
        ini_ioige2ios_dis_ioige);
 
188
 
 
189
  // ---------------------
 
190
  // 3) first renumbering
 
191
  // ---------------------
 
192
  // 3.1) geo_element renumbering
 
193
  distributor ios_ioige_ownership = ios_oige.ownership();
 
194
  array<domain_pair_type>  tmp_oige (ios_ioige_ownership);
 
195
  for (size_type ios_ioige = 0, ios_noige = ios_ioige_ownership.size();
 
196
        ios_ioige < ios_noige; ios_ioige++) {
 
197
    orientation_type orient      = ios_oige [ios_ioige].orientation();
 
198
    size_type        ios_dis_ige = ios_oige [ios_ioige].index();
 
199
    size_type        ios_ige     = ios_dis_ige - ios_ige_ownership.first_index();
 
200
    size_type        dis_ige     = omega.ios_ige2dis_ige (base::_map_dim, ios_ige);
 
201
    tmp_oige [ios_ioige].set (orient, dis_ige);
 
202
  }
 
203
  // 3.2) compute ownership for each oriented-side
 
204
  distributor ige_ownership = omega.geo_element_ownership (base::_map_dim);
 
205
  array<size_type>  partition (ios_ioige_ownership, unset);
 
206
  for (size_type ios_ioige = 0, ios_noige = ios_ioige_ownership.size();
 
207
        ios_ioige < ios_noige; ios_ioige++) {
 
208
    size_type ige = tmp_oige [ios_ioige].index();
 
209
    partition [ios_ioige] = ige_ownership.find_owner (ige);
 
210
  }
 
211
  // 3.3) repartition
 
212
  array<size_type>  ios_ioige2dis_ioige;
 
213
  array<size_type>  ioige2ios_dis_ioige;
 
214
  tmp_oige.repartition (
 
215
        partition,
 
216
        *this,
 
217
        ioige2ios_dis_ioige,
 
218
        ios_ioige2dis_ioige);
 
219
  
 
220
  // 3.4) shift from "dis_ioige" to a "ioige" numbering local to each process:
 
221
  distributor ioige_ownership = array<domain_pair_type>::ownership();
 
222
  for (size_type ioige = 0, noige = ioige_ownership.size(); ioige < noige; ioige++) {
 
223
    size_type        dis_ige = operator[] (ioige).index();
 
224
    size_type        ige = dis_ige - ige_ownership.first_index();
 
225
    operator[] (ioige).set_index (ige);
 
226
  }
 
227
  // ----------------------------------------------
 
228
  // 4) Back to "ini" renumbering: set permutations
 
229
  // ----------------------------------------------
 
230
  _ini_ioige2dis_ioige.resize (ini_ioige_ownership, unset);
 
231
  _ioige2ini_dis_ioige.resize (    ioige_ownership, unset);
 
232
  for (size_type ios_ioige = 0, ios_noige = ios_ioige_ownership.size();
 
233
        ios_ioige < ios_noige; ios_ioige++) {
 
234
    size_type ini_dis_ioige = ios_ioige2ini_dis_ioige [ios_ioige];
 
235
    size_type dis_ioige     = ios_ioige2dis_ioige     [ios_ioige];
 
236
    _ini_ioige2dis_ioige.dis_entry (ini_dis_ioige) = dis_ioige;
 
237
    _ioige2ini_dis_ioige.dis_entry (dis_ioige)     = ini_dis_ioige;
 
238
 }
 
239
 _ioige2ini_dis_ioige.dis_entry_assembly();
 
240
 _ini_ioige2dis_ioige.dis_entry_assembly();
 
241
 
 
242
  return ips;
 
243
}
 
244
template<class U>
 
245
odiststream&
 
246
domain_indirect_rep<distributed>::put (
 
247
    odiststream&                   ops,
 
248
    const geo_rep<U,distributed>& omega) const
 
249
{
 
250
  using namespace std;
 
251
  ops << "domain" << endl
 
252
      << base::_name << endl
 
253
      << "2 " << base::_map_dim << " " << dis_size() << endl;
 
254
  distributor ioige_ownership = array<domain_pair_type>::ownership();
 
255
  distributor ini_ioige_ownership = _ini_ioige2dis_ioige.ownership();
 
256
  distributor ige_ownership = omega.geo_element_ownership (base::_map_dim);
 
257
  array<domain_pair_type> ini_oige (ini_ioige_ownership);
 
258
  for (size_type ioige = 0, noige = ioige_ownership.size(); ioige < noige; ioige++) {
 
259
    size_type ini_dis_ioige = _ioige2ini_dis_ioige [ioige];
 
260
    orientation_type orient = operator[] (ioige).orientation();
 
261
    size_type        ige    = operator[] (ioige).index();
 
262
    size_type   ios_dis_ige = omega.ige2ios_dis_ige (base::_map_dim, ige);
 
263
    ini_oige.dis_entry (ini_dis_ioige) = domain_pair_type (orient, ios_dis_ige);
 
264
  }
 
265
  ini_oige.dis_entry_assembly();
 
266
  ini_oige.put_values (ops);
 
267
  return ops;
 
268
}
 
269
// ----------------------------------------------------------------------------
 
270
// instanciation in library
 
271
// ----------------------------------------------------------------------------
 
272
template
 
273
idiststream&
 
274
domain_indirect_rep<distributed>::get (
 
275
    idiststream&                       ips,
 
276
    const geo_rep<Float,distributed>& omega);
 
277
 
 
278
template
 
279
odiststream&
 
280
domain_indirect_rep<distributed>::put (
 
281
    odiststream&                       ops,
 
282
    const geo_rep<Float,distributed>& omega) const;
 
283
 
 
284
} // namespace rheolef
 
285
#endif // _RHEOLEF_HAVE_MPI