2
/// This file is part of Rheolef.
4
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
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.
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.
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
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"
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);
34
geo_basic<T,sequential>::get (idiststream& ips)
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); }
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>));
44
base::operator= (ptr);
47
/** ------------------------------------------------------------------------
48
* on any 3d geo_element K, set K.dis_iface(iloc) number
49
* ------------------------------------------------------------------------
53
geo_rep<T,sequential>::set_element_side_index (size_type side_dim)
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);
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];
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);
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);
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) {
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);
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);
115
K.face_indirect (loc_isid).set (orient, isid, shift);
121
// =========================================================================
123
// =========================================================================
124
/// @brief io for geo
127
geo_rep<T,sequential>::get (idiststream& ips)
130
check_macro (ips.good(), "bad input stream for geo.");
131
istream& is = ips.is();
132
// ------------------------------------------------------------------------
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");
140
bool do_upgrade = iorheo::getupgrade(is);
141
if (do_upgrade || hdr.need_upgrade()) {
142
return get_upgrade (ips, hdr);
144
return get_standard (ips, hdr);
149
geo_rep<T,sequential>::get_standard (idiststream& ips, const geo_header& hdr)
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);
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.");
177
base::_gs.node_ownership = base::_node.ownership();
178
// ------------------------------------------------------------------------
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();
189
base::_gs.ownership_by_dimension [base::_map_dimension] = distributor (n_elt, base::comm(), n_elt);
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
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;
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;
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) {
222
std::vector<size_t> is_vertex (nnod, 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);
228
if (base::order() > 1) {
229
for (size_type iloc = 0, nloc = K.size(); iloc < nloc; iloc++) {
230
is_vertex [K[iloc]] = 1;
234
if (base::order() == 1) {
237
n_vert = accumulate (is_vertex.begin(), is_vertex.end(), 0);
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();
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);
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) {
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();
271
size_type nsid = hdr.dis_size_by_dimension [side_dim];
272
base::_gs.ownership_by_dimension [side_dim] = distributor (nsid, base::comm(), nsid);
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);
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);
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();
300
// ----------------------------------------------------------------------------
302
// ----------------------------------------------------------------------------
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());
310
// ----------------------------------------------------------------------------
312
// ----------------------------------------------------------------------------
315
geo_rep<T,sequential>::load (std::string filename)
318
ips.open (filename, "geo");
319
check_macro(ips.good(), "\"" << filename << "[.geo[.gz]]\" not found.");
321
std::string root_name = delete_suffix (delete_suffix(filename, "gz"), "geo");
322
std::string name = get_basename (root_name);
325
// ----------------------------------------------------------------------------
326
// instanciation in library
327
// ----------------------------------------------------------------------------
328
template class geo_rep<Float,sequential>;
329
template idiststream& geo_basic<Float,sequential>::get (idiststream&);
331
} // namespace rheolef