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

« back to all changes in this revision

Viewing changes to nfem/plib/geo_seq.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 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/geo_connectivity.h"
24
 
#include "rheolef/rheostream.h"
25
 
 
26
 
namespace rheolef {
27
 
 
28
 
// =========================================================================
29
 
// geo_base_rep members
30
 
// =========================================================================
31
 
template <class T, class M>
32
 
geo_base_rep<T,M>::~geo_base_rep()
33
 
{
34
 
  warning_macro ("dstor geo_base_rep: name="<<name()<<", dim="<<dimension()<<", map_dim="<<map_dimension()<<" size="<<size());
35
 
}
36
 
template <class T, class M>
37
 
geo_base_rep<T,M>::geo_base_rep (const geo_base_rep<T,M>& omega2)
38
 
{
39
 
  error_macro ("physical copy of geo_base_rep: dim="<<omega2.dimension()<<", dis_size="<<omega2.dis_size());
40
 
}
41
 
template <class T, class M>
42
 
const domain_indirect_basic<M>&
43
 
geo_base_rep<T,M>::get_domain_indirect (const std::string& name) const
44
 
{
45
 
        for (typename std::vector<domain_indirect_basic<M> >::const_iterator
46
 
                        iter = _domains.begin(),
47
 
                        last = _domains.end(); iter != last; iter++) {
48
 
                const domain_indirect_basic<M>& dom = *iter;
49
 
                if (name == dom.name()) return dom;
50
 
        }
51
 
        error_macro ("undefined domain \""<<name<<"\" in mesh \"" << _name << "\"");
52
 
        return *(_domains.begin()); // not reached
53
 
}
54
 
// =========================================================================
55
 
// geo_rep<seq> members
56
 
// =========================================================================
57
 
/// @brief io for geo
58
 
template <class T>
59
 
        idiststream&
60
 
geo_rep<T,sequential>::get (idiststream& ips)
61
 
{
62
 
warning_macro ("get...");
63
 
  using namespace std;
64
 
  istream& is = ips.is();
65
 
warning_macro ("get [0]");
66
 
  check_macro (is.good(), "bad input stream for geo.");
67
 
 
68
 
warning_macro ("get [0.b]");
69
 
  if (!scatch(is,"\nmesh"))
70
 
        error_macro("input stream does not contains a geo.");
71
 
  //
72
 
  // 1) get header
73
 
  //
74
 
warning_macro ("get [1]");
75
 
  size_type n_vert;
76
 
  size_type n_elt;
77
 
  size_type n_edg = 0;
78
 
  size_type n_fac = 0;
79
 
 
80
 
  base::_name = "unnamed";
81
 
 
82
 
  is >> base::_version
83
 
        >> base::_dimension
84
 
        >> n_vert
85
 
        >> n_elt;
86
 
  if (base::_version < 2) {
87
 
        warning_macro ("mesh version < 2 is obsolete");
88
 
  } else {
89
 
        if (base::_dimension >= 3) {
90
 
                is >> n_fac;
91
 
        }
92
 
        if (base::_dimension >= 2) {
93
 
                is >> n_edg;
94
 
        }
95
 
  }
96
 
  //
97
 
  // 2) get coordinates
98
 
  //
99
 
warning_macro ("get [2]");
100
 
  base::_vertex.resize (n_vert);
101
 
  base::_vertex.get_values (ips, _point_get<T>(geo_base_rep<T,sequential>::_dimension));
102
 
  base::_geo_element[0].resize (base::_vertex.ownership());
103
 
  size_type first_iv = base::_vertex.ownership().first_index();
104
 
  for (size_type iv = 0; iv < n_vert; iv++) {
105
 
        base::_geo_element[0][iv] = geo_element_p(first_iv + iv);
106
 
        geo_element& P = base::_geo_element[0][iv];
107
 
        P.set_ios_dis_ie (first_iv + iv);
108
 
        P.set_dis_ie     (first_iv + iv);
109
 
  }
110
 
  check_macro (is.good(), "bad input stream for geo.");
111
 
  //
112
 
  // 3) get elements
113
 
  //
114
 
warning_macro ("get [3]");
115
 
  polymorphic_array<geo_element,sequential> elt (n_elt);
116
 
  elt.get_values (ips);
117
 
  base::_map_dimension = 0;
118
 
  for (size_type ie = 0; ie < n_elt; ie++) {
119
 
        geo_element& K = elt[ie];
120
 
        base::_map_dimension = std::max (K.dimension(), base::_map_dimension);
121
 
        K.set_ios_dis_ie (ie);
122
 
        K.set_dis_ie (ie);
123
 
  }
124
 
  if (base::_map_dimension == base::_dimension) {
125
 
        base::_geo_element[base::_map_dimension] = elt;
126
 
  }
127
 
  //
128
 
  // 4) get faces & edges
129
 
  //
130
 
warning_macro ("get [4]");
131
 
  if   (base::_version  >= 2 && base::_dimension >= 3) {
132
 
        base::_geo_element[2].resize (n_fac);
133
 
        base::_geo_element[2].get_values (ips);
134
 
        for (size_type ifac = 0; ifac < n_fac; ifac++) {
135
 
                geo_element& F = base::_geo_element[2] [ifac];
136
 
                F.set_ios_dis_ie (ifac);
137
 
                F.set_dis_ie     (ifac);
138
 
        }
139
 
  }
140
 
  if   (base::_version  >= 2 && base::_dimension >= 2) {
141
 
        base::_geo_element[1].resize (n_edg);
142
 
        base::_geo_element[1].get_values (ips);
143
 
        for (size_type iedg = 0; iedg < n_edg; iedg++) {
144
 
                geo_element& E = base::_geo_element[1] [iedg];
145
 
                E.set_ios_dis_ie (iedg);
146
 
                E.set_dis_ie     (iedg);
147
 
        }
148
 
  }
149
 
  //
150
 
  // 5) get domain, until end-of-file
151
 
  //
152
 
warning_macro ("get [5]");
153
 
  vector<set<size_type> > ball [4];
154
 
  domain_indirect_basic<sequential> dom;
155
 
  while (dom.get (ips, *this, ball)) {
156
 
warning_macro ("get [5] domain...");
157
 
     base::_domains.push_back (dom);
158
 
warning_macro ("get [5] domain done: " << dom.name());
159
 
  }
160
 
  //
161
 
  // 6) size by
162
 
  //
163
 
warning_macro ("get [6]");
164
 
  base::reset_size_by();
165
 
warning_macro ("get done");
166
 
  return ips;
167
 
}
168
 
template <class T, class M>
169
 
void
170
 
geo_base_rep<T,M>::reset_size_by()
171
 
{
172
 
  std::fill (_size_by_dimension, _size_by_dimension+4, 0);
173
 
  std::fill (_size_by_variant,   _size_by_variant+reference_element::max_size, 0);
174
 
  for (size_type dim = 0; dim <= _map_dimension; dim++) {
175
 
    _size_by_dimension [dim] = _geo_element [dim].ownership().size();
176
 
    if (dim <= 1) {
177
 
      _size_by_variant [dim] = _size_by_dimension [dim];
178
 
      continue;
179
 
    }
180
 
    const polymorphic_array<geo_element,M>& ge = _geo_element[dim]; // otherwise bug with physical copy
181
 
    for (typename polymorphic_array<geo_element,M>::const_iterator iter = ge.begin(), last = ge.end();
182
 
                iter != last; iter++) {
183
 
        const geo_element& K = (*iter);
184
 
        _size_by_variant [K.variant()]++;
185
 
    }
186
 
  }
187
 
}
188
 
template <class T>
189
 
odiststream&
190
 
geo_rep<T,sequential>::put (odiststream& ops)  const {
191
 
  using namespace std;
192
 
  size_type n_vert = base::_vertex.dis_size ();
193
 
  size_type n_elt  = base::_geo_element[base::_map_dimension].dis_size ();
194
 
  //
195
 
  // put header
196
 
  //
197
 
  ops << "#!geo" << endl
198
 
      << endl
199
 
      << "mesh" << endl
200
 
      << geo_base_rep<T,sequential>::_version << " " 
201
 
      << geo_base_rep<T,sequential>::_dimension << " " 
202
 
      << n_vert << " "
203
 
      << n_elt;
204
 
  if   (base::_version >= 2) {
205
 
    if (base::_dimension >= 3) {
206
 
      ops << " " << base::_geo_element[2].dis_size();
207
 
    }
208
 
    if (base::_dimension >= 2) {
209
 
      ops << " " << base::_geo_element[1].dis_size();
210
 
    }
211
 
  }
212
 
  ops << endl << endl;
213
 
  //
214
 
  // put vertices & elements
215
 
  //
216
 
  geo_base_rep<T,sequential>::_vertex.put_values (ops, _point_put<T>(geo_base_rep<T,sequential>::_dimension));
217
 
  ops << endl;
218
 
  base::_geo_element[base::_map_dimension].put_values (ops);
219
 
  ops << endl;
220
 
  //
221
 
  // put faces & edges
222
 
  //
223
 
  if   (base::_version >= 2 && base::_dimension >= 3) {
224
 
        base::_geo_element[2].put_values (ops); 
225
 
  }
226
 
  if   (base::_version >= 2 && base::_dimension >= 2) {
227
 
        base::_geo_element[1].put_values (ops); 
228
 
  }
229
 
  //
230
 
  // put domains
231
 
  //
232
 
  for (typename std::vector<domain_indirect_basic<sequential> >::const_iterator
233
 
        iter = base::_domains.begin(), last = base::_domains.end();
234
 
        iter != last; ++iter) {
235
 
    ops << endl;
236
 
    (*iter).put (ops);
237
 
  }
238
 
  return ops;
239
 
}
240
 
template <class T>
241
 
void
242
 
geo_rep<T,sequential>::dump (std::string name)  const {
243
 
  geo_base_rep<T,sequential>::_vertex.dump        (name + "-vert");
244
 
  base::_geo_element[base::_map_dimension].dump(name + "-elem");
245
 
}
246
 
// ----------------------------------------------------------------------------
247
 
// read from file
248
 
// ----------------------------------------------------------------------------
249
 
template <class T>
250
 
void
251
 
geo_rep<T,sequential>::load (std::string filename) 
252
 
{
253
 
  idiststream ips;
254
 
  ips.open (filename, "geo");
255
 
  check_macro(ips.good(), "\"" << filename << "[.geo[.gz]]\" not found.");
256
 
  get (ips);
257
 
  std::string root_name = delete_suffix (delete_suffix(filename, "gz"), "geo");
258
 
  std::string name = get_basename (root_name);
259
 
  base::_name = name;
260
 
}
261
 
// ----------------------------------------------------------------------------
262
 
// instanciation in library
263
 
// ----------------------------------------------------------------------------
264
 
template class geo_base_rep<Float,sequential>;
265
 
#ifdef _RHEOLEF_HAVE_MPI
266
 
template class geo_base_rep<Float,distributed>;
267
 
#endif // _RHEOLEF_HAVE_MPI
268
 
 
269
 
template class geo_rep<Float,sequential>;
270
 
 
271
 
template geo_basic<Float,sequential> compact (const geo_basic<Float,sequential>&);
272
 
 
273
 
} // namespace rheolef