~ubuntu-branches/ubuntu/raring/rheolef/raring-proposed

« back to all changes in this revision

Viewing changes to skit/ptst2/mpi_boost_polymorphe_unsrt_tst.cc

  • Committer: Package Import Robot
  • Author(s): Pierre Saramito
  • Date: 2012-04-06 09:12:21 UTC
  • mfrom: (1.1.5)
  • Revision ID: package-import@ubuntu.com-20120406091221-m58me99p1nxqui49
Tags: 6.0-1
* New upstream release 6.0 (major changes):
  - massively distributed and parallel support
  - full FEM characteristic method (Lagrange-Gakerkin method) support
  - enhanced users documentation 
  - source code supports g++-4.7 (closes: #667356)
* debian/control: dependencies for MPI distributed solvers added
* debian/rules: build commands simplified
* debian/librheolef-dev.install: man1/* to man9/* added
* debian/changelog: package description rewritted (closes: #661689)

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 distributed 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
 
// Maquette d'une classe geo avec type elements varies (polymorphe)
23
 
//  2eme version
24
 
//
25
 
// MPI: on envoie les elements par groupes de type HPTqtep homogenes
26
 
// ce qui fait que ca reste performant du point de vue MPI
27
 
//
28
 
// limitation:
29
 
//  * Les elements ont supposes ici consecutifs.
30
 
//  * Pour s'etendre au cas des variantes curvilignes (Hc,...,qc,tc)
31
 
//    cela pose un probleme, car les elements touchant la frontiere 
32
 
//    ne sont pas consecutifs par numeros.
33
 
//    On pourrait renumeroter tous les elements apres l'identification
34
 
//    de la frontiere, mais on peut eviter cette renumerotation tardive.
35
 
//  * Pour faire sauter l'hypothese des consecutif,
36
 
//    on peut aussi envoyer une paire (index,element) puis ranger x(index) = element
37
 
//  => c'est l'etape suivante
38
 
//
39
 
// test: on construit sur le proc=0 et envoie au proc=1 qui depaquette
40
 
//
41
 
// Author: Pierre.Saramito@imag.fr
42
 
//
43
 
// Date: 14 dec 2010
44
 
//
45
 
#include "rheolef/config.h"
46
 
#ifndef _RHEOLEF_HAVE_MPI
47
 
int main() { return 0; }
48
 
#else // _RHEOLEF_HAVE_MPI
49
 
 
50
 
#include "rheolef/array.h"
51
 
#include <boost/array.hpp>
52
 
#include <boost/ptr_container/ptr_container.hpp>
53
 
namespace rheolef {
54
 
// =====================================================================
55
 
class geo_element : boost::noncopyable {
56
 
public:
57
 
    typedef std::vector<int>::size_type size_type;
58
 
    typedef enum {
59
 
        p = 0, 
60
 
        e = 1, 
61
 
        t = 2, 
62
 
        q = 3, 
63
 
        T = 4, 
64
 
        P = 5, 
65
 
        H = 6, 
66
 
        enum_type_max = 7 
67
 
    } enum_type; 
68
 
 
69
 
    static char tab_name [enum_type_max];
70
 
 
71
 
    virtual ~geo_element () {}
72
 
    virtual enum_type type() const = 0;
73
 
    virtual char      name() const = 0;
74
 
    virtual size_type size() const = 0;
75
 
    virtual size_type dimension() const = 0;
76
 
    virtual size_type  operator[] (size_type i) const = 0;
77
 
    virtual size_type& operator[] (size_type i)       = 0;
78
 
 
79
 
    virtual geo_element* do_clone() const = 0;
80
 
    geo_element* clone() const { return do_clone(); }
81
 
    geo_element* new_clone (const geo_element& K) { return K.clone(); }
82
 
};
83
 
char geo_element::tab_name [geo_element::enum_type_max] = {'p', 'e', 't', 'q', 'T', 'P', 'H'};
84
 
inline
85
 
std::ostream&
86
 
operator<< (std::ostream& os, const geo_element& K) {
87
 
        os << K.name() << "\t";
88
 
        for (geo_element::size_type iloc = 0; iloc < K.size(); iloc++) {
89
 
          os << K[iloc];
90
 
          if (iloc < K.size() - 1) os << " ";
91
 
        }
92
 
        return os;
93
 
}
94
 
template <class Array>
95
 
void
96
 
dump (const geo_element& K1, Array& K2) {
97
 
        check_macro (K1.size() == K2.size(), "invalid element sizes: K1.size="<<K1.size()<<" and K2.size="<<K2.size());
98
 
        for (geo_element::size_type iloc = 0; iloc < K1.size(); iloc++) {
99
 
            K2[iloc] = K1[iloc];
100
 
        }
101
 
102
 
class geo_element_t : public geo_element {
103
 
    static const enum_type _type =  t;
104
 
    static const char      _name = 't';
105
 
    static const size_type _size =  3;
106
 
    static const size_type _dim  =  2;
107
 
public:
108
 
    typedef boost::array<size_type,_size> array_type;
109
 
 
110
 
    geo_element_t () /* : geo_element() */ {
111
 
        for (size_type iloc = 0; iloc < size(); iloc++) {
112
 
            _vertex[iloc] = std::numeric_limits<size_type>::max();
113
 
        }
114
 
    }
115
 
    geo_element_t (size_type a, size_type b, size_type c) : geo_element() {
116
 
        _vertex[0] = a;
117
 
        _vertex[1] = b;
118
 
        _vertex[2] = c;
119
 
    } 
120
 
    geo_element_t (const array_type& K) {  _vertex = K; }
121
 
    geo_element* do_clone() const {
122
 
        return new geo_element_t (_vertex[0], _vertex[1], _vertex[2]);
123
 
    }
124
 
    enum_type type() const      { return _type; };
125
 
    char      name() const      { return _name; };
126
 
    size_type size() const      { return _size; }
127
 
    size_type dimension() const { return _dim; };
128
 
    size_type  operator[] (size_type i) const { return _vertex[i]; }
129
 
    size_type& operator[] (size_type i)       { return _vertex[i]; }
130
 
    template<class Archive>
131
 
    void serialize (Archive& ar, const unsigned int version) {
132
 
        for (size_type iloc = 0; iloc < size(); iloc++) {
133
 
          ar & _vertex[iloc];
134
 
        }
135
 
    }
136
 
// data:
137
 
protected:
138
 
    array_type _vertex;
139
 
};
140
 
// =====================================================================
141
 
class geo_element_q : public geo_element {
142
 
    static const enum_type _type =  q;
143
 
    static const char      _name = 'q';
144
 
    static const size_type _size =  4;
145
 
    static const size_type _dim  =  2;
146
 
public:
147
 
    typedef boost::array<size_type,_size> array_type;
148
 
    geo_element_q () : geo_element() { 
149
 
        for (size_type iloc = 0; iloc < size(); iloc++) {
150
 
            _vertex[iloc] = std::numeric_limits<size_type>::max();
151
 
        }
152
 
    }
153
 
    geo_element_q (size_type a, size_type b, size_type c, size_type d) : geo_element() {
154
 
        _vertex[0] = a;
155
 
        _vertex[1] = b;
156
 
        _vertex[2] = c;
157
 
        _vertex[3] = d;
158
 
    } 
159
 
    geo_element_q (const array_type& K) {  _vertex = K; }
160
 
    geo_element* do_clone() const {
161
 
        return new geo_element_q (_vertex[0], _vertex[1], _vertex[2], _vertex[3]);
162
 
    }
163
 
    enum_type type() const      { return _type; };
164
 
    char      name() const      { return _name; };
165
 
    size_type size() const      { return _size; }
166
 
    size_type dimension() const { return _dim; };
167
 
    size_type  operator[] (size_type i) const { return _vertex[i]; }
168
 
    size_type& operator[] (size_type i)       { return _vertex[i]; }
169
 
    template<class Archive> 
170
 
    void serialize (Archive& ar, const unsigned int version) {
171
 
        for (size_type iloc = 0; iloc < size(); iloc++) {
172
 
          ar & _vertex[iloc];
173
 
        }
174
 
    }
175
 
// data:
176
 
protected:
177
 
    array_type _vertex;
178
 
};
179
 
} // namespace rheolef
180
 
// =====================================================================
181
 
// Some serializable types, like geo_element, have a fixed amount of data stored at fixed field positions.
182
 
// When this is the case, boost::mpi can optimize their serialization and transmission to avoid extraneous 
183
 
// copy operations.
184
 
// To enable this optimization, we specialize the type trait is_mpi_datatype, e.g.:
185
 
namespace boost {
186
 
 namespace mpi {
187
 
  template <> struct is_mpi_datatype<rheolef::geo_element_t> : mpl::true_ { };
188
 
  template <> struct is_mpi_datatype<rheolef::geo_element_q> : mpl::true_ { };
189
 
 } // namespace mpi
190
 
} // namespace boost
191
 
 
192
 
// =====================================================================
193
 
// classe polymorphe & MPI
194
 
// =====================================================================
195
 
 
196
 
using namespace rheolef;
197
 
using namespace std;
198
 
 
199
 
struct my_polymorphic_geo {
200
 
  typedef boost::ptr_vector<geo_element>::size_type size_type;
201
 
  typedef boost::ptr_vector<geo_element>::iterator iterator;
202
 
  typedef boost::ptr_vector<geo_element>::const_iterator const_iterator;
203
 
  my_polymorphic_geo () : _x(), _comm() {
204
 
    fill (_count_by_type, _count_by_type+geo_element::enum_type_max, 0);
205
 
    if (_comm.rank() == 0) {
206
 
      _x.push_back (new geo_element_t (1,2,3));
207
 
      _x.push_back (new geo_element_q (11,12,13,14));
208
 
      _x.push_back (new geo_element_t (4,5,6));
209
 
      _x.push_back (new geo_element_q (15,16,17,18));
210
 
      _x.push_back (new geo_element_t (7,8,9));
211
 
      _count_by_type [geo_element::t] = 3;
212
 
      _count_by_type [geo_element::q] = 2;
213
 
    }
214
 
  }
215
 
  void put () const {
216
 
    cout << "geo " << _x.size() << endl;
217
 
    for (size_t i = 0; i < _x.size(); i++) {
218
 
      const geo_element& K = _x[i];
219
 
#ifndef TO_CLEAN
220
 
      const std::type_info& rtti = typeid(K);
221
 
      cerr << "type " << rtti.name() << ": " << K << endl;
222
 
#endif // TO_CLEAN
223
 
      cout << K << endl;
224
 
    }
225
 
  }
226
 
  void test_mpi () {
227
 
    pair<size_type,size_type> count_tag_pair [size_type(geo_element::enum_type_max)];
228
 
    int tag_count = 0;
229
 
    if (_comm.rank() == 0) {
230
 
      for (size_type i = 0; i < size_type(geo_element::enum_type_max); i++) {
231
 
        size_type tag = 100+i;
232
 
        count_tag_pair [i] = make_pair(_count_by_type[i], tag);
233
 
      }
234
 
      _comm.send (1, tag_count, count_tag_pair, size_type(geo_element::enum_type_max));
235
 
      // envoi
236
 
      vector<pair<size_type,geo_element_t::array_type> > buffer_t (_count_by_type[geo_element::t]);
237
 
      vector<pair<size_type,geo_element_q::array_type> > buffer_q (_count_by_type[geo_element::q]);
238
 
      vector<size_type> last_by_type (geo_element::enum_type_max, 0);
239
 
      for (size_type ie = 0; ie < _x.size(); ie++) {
240
 
        const geo_element& K = _x[ie];
241
 
        size_type ibuf = last_by_type[K.type()];
242
 
        last_by_type[K.type()]++;
243
 
        switch (K.type()) {
244
 
          // -----------------
245
 
          // dump triangles:
246
 
          // -----------------
247
 
          case geo_element::t: {
248
 
            buffer_t[ibuf].first = ie;
249
 
            dump (_x[ie], buffer_t[ibuf].second);
250
 
            break;
251
 
          }
252
 
          // -----------------
253
 
          // dump quadrangles:
254
 
          // -----------------
255
 
          case geo_element::q: {
256
 
            buffer_q[ibuf].first = ie;
257
 
            dump (_x[ie], buffer_q[ibuf].second);
258
 
            break;
259
 
          }
260
 
          default : {
261
 
            error_macro ("unexpected element " << geo_element::tab_name[K.type()]);
262
 
          }
263
 
        }
264
 
      }
265
 
      // -----------------
266
 
      // send triangles:
267
 
      // -----------------
268
 
      if (count_tag_pair[geo_element::t].first != 0) {
269
 
        _comm.send (1,
270
 
                count_tag_pair[geo_element::t].second,
271
 
                buffer_t.begin().operator->(),  
272
 
                count_tag_pair[geo_element::t].first);
273
 
      }
274
 
      // -----------------
275
 
      // send quadrangles:
276
 
      // -----------------
277
 
      if (count_tag_pair[geo_element::q].first != 0) {
278
 
        _comm.send (1,
279
 
                count_tag_pair[geo_element::q].second,
280
 
                buffer_q.begin().operator->(),  
281
 
                count_tag_pair[geo_element::q].first);
282
 
      }
283
 
    } else if (_comm.rank() == 1) {
284
 
      _comm.recv (0, tag_count, count_tag_pair, size_type(geo_element::enum_type_max));
285
 
      vector<pair<size_type,geo_element_t::array_type> > buffer_t (count_tag_pair[geo_element::t].first);
286
 
      vector<pair<size_type,geo_element_q::array_type> > buffer_q (count_tag_pair[geo_element::q].first);
287
 
      size_type total_size = 0;
288
 
      for (size_type i = 0; i < size_type(geo_element::enum_type_max); i++) {
289
 
        _count_by_type[i] = count_tag_pair [i].first;
290
 
        total_size += _count_by_type[i];
291
 
      }
292
 
      _x.reserve (total_size);
293
 
      // resize...
294
 
#ifdef TODO
295
 
      _x.resize (total_size);
296
 
#else // TODO
297
 
      for (size_type j = 0; j < total_size; j++) {
298
 
        _x.push_back (new geo_element_t);
299
 
      }
300
 
#endif // TODO
301
 
      // -----------------
302
 
      // recv triangles:
303
 
      // -----------------
304
 
      if (count_tag_pair[geo_element::t].first != 0) {
305
 
        _comm.recv (0,
306
 
                count_tag_pair[geo_element::t].second,
307
 
                buffer_t.begin().operator->(),  
308
 
                count_tag_pair[geo_element::t].first);
309
 
        for (size_type ibuf = 0; ibuf < buffer_t.size(); ibuf++) {
310
 
          size_type ie = buffer_t[ibuf].first;
311
 
          //_x.insert (_x.begin()+ie, new geo_element_t (buffer_t[ibuf].second));
312
 
          //warning_macro ("replace["<<ie<<"]="<<buffer_t[ibuf].second);
313
 
          check_macro (!_x.is_null(ie), "is_null "<<ie);
314
 
          _x.replace (_x.begin()+ie, new geo_element_t (buffer_t[ibuf].second));
315
 
        }
316
 
      }
317
 
      // -----------------
318
 
      // recv quadrangles:
319
 
      // -----------------
320
 
      if (count_tag_pair[geo_element::q].first != 0) {
321
 
        _comm.recv (0,
322
 
                count_tag_pair[geo_element::q].second,
323
 
                buffer_q.begin().operator->(),  
324
 
                count_tag_pair[geo_element::q].first);
325
 
        for (size_type ibuf = 0; ibuf < buffer_q.size(); ibuf++) {
326
 
          size_type ie = buffer_q[ibuf].first;
327
 
          //_x.insert (_x.begin()+ie, new geo_element_q (buffer_q[ibuf].second));
328
 
          //warning_macro ("replace["<<ie<<"]="<<buffer_q[ibuf].second);
329
 
          check_macro (!_x.is_null(ie), "is_null "<<ie);
330
 
          _x.replace (_x.begin()+ie, new geo_element_q(buffer_q[ibuf].second));
331
 
        }
332
 
      }
333
 
      put();
334
 
    }
335
 
  }
336
 
// data:
337
 
  boost::ptr_vector<geo_element> _x;
338
 
  mpi::communicator              _comm;
339
 
  size_type _count_by_type [geo_element::enum_type_max];
340
 
};
341
 
int main(int argc, char**argv) {
342
 
  environment distributed(argc, argv);
343
 
  mpi::communicator comm;
344
 
  if (comm.size() == 1) {
345
 
        cerr << "not sequential : use mpirun -np 2 !" << endl;
346
 
        return 0;
347
 
  }
348
 
  my_polymorphic_geo omega;
349
 
  omega.test_mpi();
350
 
}
351
 
#endif // _RHEOLEF_HAVE_MPI