← Back to branch summary

~ubuntu-branches/ubuntu/trusty/rheolef/trusty

« back to all changes in this revision

Viewing changes to nfem/plib/geo_seq_get.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
http://bugs.debian.org/667356
http://bugs.debian.org/661689
* 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:
nfem/plib/geo_seq_get.cc
 
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
#include "rheolef/geo.h"
 
22
#include "rheolef/geo_domain.h"
 
23
#include "rheolef/rheostream.h"
 
24
#include "rheolef/iorheo.h"
 
25
#include "rheolef/rheostream.h"
 
26
 
 
27
namespace rheolef {
 
28
 
 
29
template <class T> idiststream& geo_get_vtk  (idiststream& ips, geo_basic<T,sequential>& omega);
 
30
template <class T> idiststream& geo_get_bamg (idiststream& ips, geo_basic<T,sequential>& omega);
 
31
 
 
32
template <class T>
 
33
idiststream&
 
34
geo_basic<T,sequential>::get (idiststream& ips)
 
35
{
 
36
  iorheo::flag_type format = iorheo::flags(ips.is()) & iorheo::format_field;
 
37
  if (format [iorheo::vtk])     { return geo_get_vtk     (ips,*this); }
 
38
  if (format [iorheo::bamg])    { return geo_get_bamg    (ips,*this); }
 
39
 
 
40
  // else: standard .geo format:
 
41
  // allocate a new geo_rep object (TODO: do a dynamic_cast ?)
 
42
  geo_rep<T,sequential>* ptr = new_macro((geo_rep<T,sequential>));
 
43
  ptr->get (ips);
 
44
  base::operator= (ptr);
 
45
  return ips;
 
46
}
 
47
/** ------------------------------------------------------------------------
 
48
 * on any 3d geo_element K, set K.dis_iface(iloc) number
 
49
 * ------------------------------------------------------------------------
 
50
 */
 
51
template <class T>
 
52
void
 
53
geo_rep<T,sequential>::set_element_side_index (size_type side_dim)
 
54
{
 
55
  if (map_dimension() <= side_dim) return;
 
56
  // ------------------------------------------------------------------------
 
57
  // 1) ball(X) := { E; X is a vertex of E }
 
58
  // ------------------------------------------------------------------------
 
59
  index_set empty_set; // TODO: add a global allocator to empty_set
 
60
  array<index_set,sequential> ball (geo_element_ownership(0), empty_set);
 
61
  {
 
62
    size_type isid = 0;
 
63
    for (const_iterator iter = begin(side_dim), last = end(side_dim); iter != last; iter++, isid++) {
 
64
      const geo_element& S = *iter;
 
65
      for (size_type iloc = 0, nloc = S.size(); iloc < nloc; iloc++) {
 
66
        size_type iv = S[iloc];
 
67
        ball [iv] += isid;
 
68
      }
 
69
    }
 
70
  }
 
71
  // ------------------------------------------------------------------------
 
72
  // 2) pour K dans partition(iproc)
 
73
  //      pour (dis_A,dis_B) arete de K
 
74
  //        set = dis_ball(dis_A) inter dis_ball(dis_B) = {dis_iedg}
 
75
  //        E = dis_edges(dis_iedg)
 
76
  //        => on numerote dis_iedg cette arete dans le geo_element K
 
77
  //        et on indique son orient en comparant a E, arete qui definit l'orient
 
78
  // ------------------------------------------------------------------------
 
79
  for (size_type dim = side_dim+1; dim <= base::_map_dimension; dim++) {
 
80
    for (iterator iter = begin(dim), last = end(dim); iter != last; iter++) {
 
81
      geo_element& K = *iter;
 
82
      for (size_type loc_isid = 0, loc_nsid = K.n_subgeo(side_dim); loc_isid < loc_nsid; loc_isid++) {
 
83
        size_type loc_jv0 = K.subgeo_local_vertex (side_dim, loc_isid, 0);
 
84
        size_type     jv0 = K[loc_jv0];
 
85
        index_set isid_set = ball [jv0]; // copy: will be intersected
 
86
        for (size_type sid_jloc = 1, sid_nloc = K.subgeo_size (side_dim, loc_isid); sid_jloc < sid_nloc; sid_jloc++) { 
 
87
          size_type loc_jv = K.subgeo_local_vertex (side_dim, loc_isid, sid_jloc);
 
88
          size_type     jv = K[loc_jv];
 
89
          const index_set& ball_jv = ball [jv];
 
90
          isid_set.inplace_intersection (ball_jv);
 
91
        }
 
92
        check_macro (isid_set.size() == 1, "connectivity: side not found in the side set");
 
93
        size_type isid = *(isid_set.begin());
 
94
        const geo_element& S = get_geo_element(side_dim,isid);
 
95
        if (side_dim == 1) {
 
96
          // side: edge
 
97
          size_type jv1 = K [K.subgeo_local_vertex (side_dim, loc_isid, 1)];
 
98
          geo_element::orientation_type orient = S.get_edge_orientation (jv0, jv1);
 
99
          K.edge_indirect (loc_isid).set (orient, isid);
 
100
        } else { // side_dim == 2 
 
101
          geo_element::orientation_type orient;
 
102
          geo_element::shift_type       shift;
 
103
          if (K.subgeo_size (side_dim, loc_isid) == 3) {
 
104
            // side: triangle
 
105
            size_type jv1 = K [K.subgeo_local_vertex (side_dim, loc_isid, 1)];
 
106
            size_type jv2 = K [K.subgeo_local_vertex (side_dim, loc_isid, 2)];
 
107
            S.get_orientation_and_shift (jv0, jv1, jv2, orient, shift);
 
108
          } else {
 
109
            // side: quadrangle
 
110
            size_type jv1 = K [K.subgeo_local_vertex (side_dim, loc_isid, 1)];
 
111
            size_type jv2 = K [K.subgeo_local_vertex (side_dim, loc_isid, 2)];
 
112
            size_type jv3 = K [K.subgeo_local_vertex (side_dim, loc_isid, 3)];
 
113
            S.get_orientation_and_shift (jv0, jv1, jv2, jv3, orient, shift);
 
114
          }
 
115
          K.face_indirect (loc_isid).set (orient, isid, shift);
 
116
        }
 
117
      }
 
118
    }
 
119
  }
 
120
}
 
121
// =========================================================================
 
122
// get
 
123
// =========================================================================
 
124
/// @brief io for geo
 
125
template <class T>
 
126
idiststream&
 
127
geo_rep<T,sequential>::get (idiststream& ips)
 
128
{
 
129
  using namespace std;
 
130
  check_macro (ips.good(), "bad input stream for geo.");
 
131
  istream& is = ips.is();
 
132
  // ------------------------------------------------------------------------
 
133
  // 0) get header
 
134
  // ------------------------------------------------------------------------
 
135
  check_macro (dis_scatch(ips,ips.comm(),"\nmesh"), "input stream does not contains a geo.");
 
136
  ips >> base::_version;
 
137
  check_macro (base::_version == 4, "mesh format version 4 expected, but format version " << base::_version << " founded");
 
138
  geo_header hdr;
 
139
  ips >> hdr;
 
140
  bool do_upgrade = iorheo::getupgrade(is);
 
141
  if (do_upgrade || hdr.need_upgrade()) {
 
142
    return get_upgrade  (ips, hdr);
 
143
  } else {
 
144
    return get_standard (ips, hdr);
 
145
  }
 
146
}
 
147
template <class T>
 
148
idiststream&
 
149
geo_rep<T,sequential>::get_standard (idiststream& ips, const geo_header& hdr)
 
150
{
 
151
  using namespace std;
 
152
  check_macro (ips.good(), "bad input stream for geo.");
 
153
  istream& is = ips.is();
 
154
  // ------------------------------------------------------------------------
 
155
  // 1) store header infos in geo
 
156
  // ------------------------------------------------------------------------
 
157
  base::_have_connectivity = true;
 
158
  base::_name          = "unnamed";
 
159
  base::_dimension     = hdr.dimension;
 
160
  base::_map_dimension = hdr.map_dimension;
 
161
  base::_sys_coord     = hdr.sys_coord;
 
162
  base::_numbering.set_degree (hdr.order);
 
163
 
 
164
  // compute map_dimension = max dimension with non-empty geo_element set
 
165
  size_type nnod  = hdr.dis_size_by_dimension [0];
 
166
  size_type n_edg = hdr.dis_size_by_dimension [1];
 
167
  size_type n_fac = hdr.dis_size_by_dimension [2];
 
168
  size_type n_elt = hdr.dis_size_by_dimension [base::_map_dimension];
 
169
  // ------------------------------------------------------------------------
 
170
  // 2) get coordinates
 
171
  // ------------------------------------------------------------------------
 
172
  base::_node.resize (nnod);
 
173
  if (base::_dimension > 0) {
 
174
    base::_node.get_values (ips, _point_get<T>(geo_base_rep<T,sequential>::_dimension));
 
175
    check_macro (ips.good(), "bad input stream for geo.");
 
176
  }
 
177
  base::_gs.node_ownership = base::_node.ownership();
 
178
  // ------------------------------------------------------------------------
 
179
  // 3) get elements
 
180
  // ------------------------------------------------------------------------
 
181
  if (base::_map_dimension > 0) {
 
182
    for (size_type variant = reference_element::first_variant_by_dimension(base::_map_dimension);
 
183
                   variant < reference_element:: last_variant_by_dimension(base::_map_dimension); variant++) {
 
184
      geo_element::parameter_type param (variant, 1);
 
185
      base::_geo_element [variant].resize (hdr.dis_size_by_variant [variant], param);
 
186
      base::_geo_element [variant].get_values (ips);
 
187
      base::_gs.ownership_by_variant [variant] = base::_geo_element [variant].ownership();
 
188
    }
 
189
    base::_gs.ownership_by_dimension [base::_map_dimension] = distributor (n_elt, base::comm(), n_elt);
 
190
  }
 
191
  // ------------------------------------------------------------------------
 
192
  // 4) check that nodes are numbered by increasing node_subgeo_dim
 
193
  // ------------------------------------------------------------------------
 
194
  // ICI va devenir obsolete car les noeuds seront numerotes par _numbering=Pk_numbering
 
195
  {
 
196
    std::vector<size_type> node_subgeo_dim (nnod, size_type(-1));
 
197
    size_type prev_variant = 0;
 
198
    for (iterator iter = begin(base::_map_dimension), last = end(base::_map_dimension); iter != last; iter++) {
 
199
      geo_element& K = *iter;
 
200
      check_macro (prev_variant <= K.variant(), "elements should be numbered by increasing variants (petqTPH)");
 
201
      prev_variant = K.variant();
 
202
      for (size_type d = 0; d <= base::_map_dimension; d++) { 
 
203
        for (size_type loc_inod = K.first_inod(d), loc_nnod = K.last_inod(d); loc_inod < loc_nnod; loc_inod++) {
 
204
          node_subgeo_dim [K[loc_inod]] = d;
 
205
        }
 
206
      }
 
207
    }
 
208
    size_type prev_node_dim = 0;
 
209
    for (typename std::vector<size_type>::const_iterator iter = node_subgeo_dim.begin(), last = node_subgeo_dim.end();
 
210
                iter != last; iter++) {
 
211
      check_macro (prev_node_dim <= *iter, "nodes should be numbered by increasing subgeo dimension");
 
212
      prev_node_dim = *iter;
 
213
    }
 
214
  }
 
215
  // ------------------------------------------------------------------------
 
216
  // 5) compute n_vert (n_vert < nnod when order > 1) and set element indexes (K.dis_ie & K.ios_dis_ie)
 
217
  // ------------------------------------------------------------------------
 
218
  size_type n_vert = 0;
 
219
  if (base::_map_dimension == 0) {
 
220
    n_vert = nnod;
 
221
  } else {
 
222
    std::vector<size_t> is_vertex (nnod, 0);
 
223
    size_type ie = 0;
 
224
    for (iterator iter = begin(base::_map_dimension), last = end(base::_map_dimension); iter != last; iter++, ie++) {
 
225
        geo_element& K = *iter;
 
226
        K.set_ios_dis_ie (ie);
 
227
        K.set_dis_ie (ie);
 
228
        if (base::order() > 1) {
 
229
          for (size_type iloc = 0, nloc = K.size(); iloc < nloc; iloc++) {
 
230
            is_vertex [K[iloc]] = 1;
 
231
          }
 
232
        }
 
233
    }
 
234
    if (base::order() == 1) {
 
235
      n_vert = nnod;
 
236
    } else {
 
237
      n_vert = accumulate (is_vertex.begin(), is_vertex.end(), 0);
 
238
    }
 
239
  }
 
240
  // ------------------------------------------------------------------------
 
241
  // 6) create vertex-element (0d elements)
 
242
  // ------------------------------------------------------------------------
 
243
  geo_element::parameter_type param (reference_element::p, 1);
 
244
  base::_geo_element [reference_element::p].resize (n_vert, param);
 
245
  size_type first_iv = base::_node.ownership().first_index();
 
246
  {
 
247
    size_type iv = 0;
 
248
    for (iterator iter = begin(0), last = end(0); iter != last; iter++, iv++) {
 
249
      geo_element& P = *iter;
 
250
      P[0]            = first_iv + iv;
 
251
      P.set_dis_ie     (first_iv + iv); // TODO: P[0] & dis_ie redundant for `p'
 
252
      P.set_ios_dis_ie (first_iv + iv);
 
253
    }
 
254
  }
 
255
  // ownership_by_dimension[0]: used by connectivity
 
256
  base::_gs.ownership_by_variant   [reference_element::p] = base::_geo_element [reference_element::p].ownership();
 
257
  base::_gs.ownership_by_dimension [0] = base::_geo_element [reference_element::p].ownership();
 
258
  // ------------------------------------------------------------------------
 
259
  // 7) get faces & edge
 
260
  // ------------------------------------------------------------------------
 
261
  if (base::_map_dimension > 0) {
 
262
    
 
263
    for (size_type side_dim = base::_map_dimension - 1; side_dim >= 1; side_dim--) {
 
264
      for (size_type variant = reference_element::first_variant_by_dimension(side_dim);
 
265
                     variant < reference_element:: last_variant_by_dimension(side_dim); variant++) {
 
266
        geo_element::parameter_type param (variant, 1);
 
267
        base::_geo_element [variant].resize (hdr.dis_size_by_variant [variant], param);
 
268
        base::_geo_element [variant].get_values (ips);
 
269
        base::_gs.ownership_by_variant [variant] = base::_geo_element [variant].ownership();
 
270
      }
 
271
      size_type nsid = hdr.dis_size_by_dimension [side_dim];
 
272
      base::_gs.ownership_by_dimension [side_dim] = distributor (nsid, base::comm(), nsid);
 
273
      size_type isid = 0;
 
274
      for (iterator iter = begin(side_dim), last = end(side_dim); iter != last; iter++, isid++) {
 
275
        geo_element& S = *iter;
 
276
        S.set_ios_dis_ie (isid);
 
277
        S.set_dis_ie     (isid);
 
278
      }
 
279
    }
 
280
  }
 
281
  // ------------------------------------------------------------------------
 
282
  // 8) get domain, until end-of-file
 
283
  // ------------------------------------------------------------------------
 
284
  vector<index_set> ball [4];
 
285
  domain_indirect_basic<sequential> dom;
 
286
  while (dom.get (ips, *this, ball)) {
 
287
     base::_domains.push_back (dom);
 
288
  }
 
289
  // ------------------------------------------------------------------------
 
290
  // 9) set indexes on faces and edges of elements, for P2 approx
 
291
  // ------------------------------------------------------------------------
 
292
  set_element_side_index (1);
 
293
  set_element_side_index (2);
 
294
  // ------------------------------------------------------------------------
 
295
  // 10) bounding box: _xmin, _xmax
 
296
  // ------------------------------------------------------------------------
 
297
  base::compute_bbox();
 
298
  return ips;
 
299
}
 
300
// ----------------------------------------------------------------------------
 
301
// dump
 
302
// ----------------------------------------------------------------------------
 
303
template <class T>
 
304
void
 
305
geo_rep<T,sequential>::dump (std::string name)  const {
 
306
  std::ofstream os ((name + "-dump.geo").c_str());
 
307
  odiststream ods (os, base::_node.ownership().comm());
 
308
  put_geo (ods);
 
309
}
 
310
// ----------------------------------------------------------------------------
 
311
// read from file
 
312
// ----------------------------------------------------------------------------
 
313
template <class T>
 
314
void
 
315
geo_rep<T,sequential>::load (std::string filename) 
 
316
{
 
317
  idiststream ips;
 
318
  ips.open (filename, "geo");
 
319
  check_macro(ips.good(), "\"" << filename << "[.geo[.gz]]\" not found.");
 
320
  get (ips);
 
321
  std::string root_name = delete_suffix (delete_suffix(filename, "gz"), "geo");
 
322
  std::string name = get_basename (root_name);
 
323
  base::_name = name;
 
324
}
 
325
// ----------------------------------------------------------------------------
 
326
// instanciation in library
 
327
// ----------------------------------------------------------------------------
 
328
template class geo_rep<Float,sequential>;
 
329
template idiststream& geo_basic<Float,sequential>::get (idiststream&);
 
330
 
 
331
} // namespace rheolef