rheolef  6.5
geo_seq_get.cc
Go to the documentation of this file.
1 #include "rheolef/geo.h"
2 #include "rheolef/geo_domain.h"
3 #include "rheolef/rheostream.h"
4 #include "rheolef/iorheo.h"
5 #include "rheolef/rheostream.h"
6 
7 namespace rheolef {
8 
9 template <class T> idiststream& geo_get_vtk (idiststream& ips, geo_basic<T,sequential>& omega);
10 template <class T> idiststream& geo_get_bamg (idiststream& ips, geo_basic<T,sequential>& omega);
11 
12 template <class T>
15 {
17  if (format [iorheo::vtk]) { return geo_get_vtk (ips,*this); }
18  if (format [iorheo::bamg]) { return geo_get_bamg (ips,*this); }
19 
21  ptr->get (ips);
22  base::operator= (ptr);
23  return ips;
24 }
29 template <class T>
30 void
32 {
33  if (map_dimension() <= side_dim) return;
34  index_set empty_set; // TODO: add a global allocator to empty_set
35  array<index_set,sequential> ball (geo_element_ownership(0), empty_set);
36  {
37  size_type isid = 0;
38  for (const_iterator iter = begin(side_dim), last = end(side_dim); iter != last; iter++, isid++) {
39  const geo_element& S = *iter;
40  for (size_type iloc = 0, nloc = S.size(); iloc < nloc; iloc++) {
41  size_type iv = S[iloc];
42  ball [iv] += isid;
43  }
44  }
45  }
46  // => on numerote dis_iedg cette arete dans le geo_element K
47  for (size_type dim = side_dim+1; dim <= base::_map_dimension; dim++) {
48  for (iterator iter = begin(dim), last = end(dim); iter != last; iter++) {
49  geo_element& K = *iter;
50  for (size_type loc_isid = 0, loc_nsid = K.n_subgeo(side_dim); loc_isid < loc_nsid; loc_isid++) {
51  size_type loc_jv0 = K.subgeo_local_vertex (side_dim, loc_isid, 0);
52  size_type jv0 = K[loc_jv0];
53  index_set isid_set = ball [jv0]; // copy: will be intersected
54  for (size_type sid_jloc = 1, sid_nloc = K.subgeo_size (side_dim, loc_isid); sid_jloc < sid_nloc; sid_jloc++) {
55  size_type loc_jv = K.subgeo_local_vertex (side_dim, loc_isid, sid_jloc);
56  size_type jv = K[loc_jv];
57  const index_set& ball_jv = ball [jv];
58  isid_set.inplace_intersection (ball_jv);
59  }
60  check_macro (isid_set.size() == 1, "connectivity: side not found in the side set");
61  size_type isid = *(isid_set.begin());
62  const geo_element& S = get_geo_element(side_dim,isid);
63  if (side_dim == 1) {
64  size_type jv1 = K [K.subgeo_local_vertex (side_dim, loc_isid, 1)];
66  K.edge_indirect (loc_isid).set (orient, isid);
67  } else { // side_dim == 2
70  if (K.subgeo_size (side_dim, loc_isid) == 3) {
71  size_type jv1 = K [K.subgeo_local_vertex (side_dim, loc_isid, 1)];
72  size_type jv2 = K [K.subgeo_local_vertex (side_dim, loc_isid, 2)];
73  S.get_orientation_and_shift (jv0, jv1, jv2, orient, shift);
74  } else {
75  size_type jv1 = K [K.subgeo_local_vertex (side_dim, loc_isid, 1)];
76  size_type jv2 = K [K.subgeo_local_vertex (side_dim, loc_isid, 2)];
77  size_type jv3 = K [K.subgeo_local_vertex (side_dim, loc_isid, 3)];
78  S.get_orientation_and_shift (jv0, jv1, jv2, jv3, orient, shift);
79  }
80  K.face_indirect (loc_isid).set (orient, isid, shift);
81  }
82  }
83  }
84  }
85 }
86 template <class T>
89 {
90  using namespace std;
91  check_macro (ips.good(), "bad input stream for geo.");
92  istream& is = ips.is();
93  check_macro (dis_scatch(ips,ips.comm(),"\nmesh"), "input stream does not contains a geo.");
94  ips >> base::_version;
95  check_macro (base::_version == 4, "mesh format version 4 expected, but format version " << base::_version << " founded");
96  geo_header hdr;
97  ips >> hdr;
98  bool do_upgrade = iorheo::getupgrade(is);
99  if (do_upgrade || hdr.need_upgrade()) {
100  return get_upgrade (ips, hdr);
101  } else {
102  return get_standard (ips, hdr);
103  }
104 }
105 template <class T>
108 {
109  using namespace std;
110  check_macro (ips.good(), "bad input stream for geo.");
111  istream& is = ips.is();
112  base::_have_connectivity = true;
113  base::_name = "unnamed";
114  base::_dimension = hdr.dimension;
115  base::_map_dimension = hdr.map_dimension;
116  base::_sys_coord = hdr.sys_coord;
117  base::_numbering.set_degree (hdr.order);
118 
119  size_type nnod = hdr.dis_size_by_dimension [0];
120  size_type n_edg = hdr.dis_size_by_dimension [1];
121  size_type n_fac = hdr.dis_size_by_dimension [2];
122  size_type n_elt = hdr.dis_size_by_dimension [base::_map_dimension];
123  base::_node.resize (nnod);
124  if (base::_dimension > 0) {
125  base::_node.get_values (ips, _point_get<T>(geo_base_rep<T,sequential>::_dimension));
126  check_macro (ips.good(), "bad input stream for geo.");
127  }
128  base::_gs.node_ownership = base::_node.ownership();
129  if (base::_map_dimension > 0) {
130  for (size_type variant = reference_element::first_variant_by_dimension(base::_map_dimension);
131  variant < reference_element:: last_variant_by_dimension(base::_map_dimension); variant++) {
132  geo_element::parameter_type param (variant, 1);
133  base::_geo_element [variant].resize (hdr.dis_size_by_variant [variant], param);
134  base::_geo_element [variant].get_values (ips);
135  base::_gs.ownership_by_variant [variant] = base::_geo_element [variant].ownership();
136  }
137  base::_gs.ownership_by_dimension [base::_map_dimension] = distributor (n_elt, base::comm(), n_elt);
138  }
139  {
140  std::vector<size_type> node_subgeo_dim (nnod, size_type(-1));
141  size_type prev_variant = 0;
142  for (iterator iter = begin(base::_map_dimension), last = end(base::_map_dimension); iter != last; iter++) {
143  geo_element& K = *iter;
144  check_macro (prev_variant <= K.variant(), "elements should be numbered by increasing variants (petqTPH)");
145  prev_variant = K.variant();
146  for (size_type d = 0; d <= base::_map_dimension; d++) {
147  for (size_type loc_inod = K.first_inod(d), loc_nnod = K.last_inod(d); loc_inod < loc_nnod; loc_inod++) {
148  node_subgeo_dim [K[loc_inod]] = d;
149  }
150  }
151  }
152  size_type prev_node_dim = 0;
153  for (typename std::vector<size_type>::const_iterator iter = node_subgeo_dim.begin(), last = node_subgeo_dim.end();
154  iter != last; iter++) {
155  check_macro (prev_node_dim <= *iter, "nodes should be numbered by increasing subgeo dimension");
156  prev_node_dim = *iter;
157  }
158  }
159  // 5) compute n_vert (n_vert < nnod when order > 1) and set element indexes (K.dis_ie & K.ios_dis_ie)
160  size_type n_vert = 0;
161  if (base::_map_dimension == 0) {
162  n_vert = nnod;
163  } else {
164  std::vector<size_t> is_vertex (nnod, 0);
165  size_type ie = 0;
166  for (iterator iter = begin(base::_map_dimension), last = end(base::_map_dimension); iter != last; iter++, ie++) {
167  geo_element& K = *iter;
168  K.set_ios_dis_ie (ie);
169  K.set_dis_ie (ie);
170  if (base::order() > 1) {
171  for (size_type iloc = 0, nloc = K.size(); iloc < nloc; iloc++) {
172  is_vertex [K[iloc]] = 1;
173  }
174  }
175  }
176  if (base::order() == 1) {
177  n_vert = nnod;
178  } else {
179  n_vert = accumulate (is_vertex.begin(), is_vertex.end(), 0);
180  }
181  }
183  base::_geo_element [reference_element::p].resize (n_vert, param);
184  size_type first_iv = base::_node.ownership().first_index();
185  {
186  size_type iv = 0;
187  for (iterator iter = begin(0), last = end(0); iter != last; iter++, iv++) {
188  geo_element& P = *iter;
189  P[0] = first_iv + iv;
190  P.set_dis_ie (first_iv + iv); // TODO: P[0] & dis_ie redundant for `p'
191  P.set_ios_dis_ie (first_iv + iv);
192  }
193  }
194  base::_gs.ownership_by_variant [reference_element::p] = base::_geo_element [reference_element::p].ownership();
195  base::_gs.ownership_by_dimension [0] = base::_geo_element [reference_element::p].ownership();
196  if (base::_map_dimension > 0) {
197 
198  for (size_type side_dim = base::_map_dimension - 1; side_dim >= 1; side_dim--) {
200  variant < reference_element:: last_variant_by_dimension(side_dim); variant++) {
201  geo_element::parameter_type param (variant, 1);
202  base::_geo_element [variant].resize (hdr.dis_size_by_variant [variant], param);
203  base::_geo_element [variant].get_values (ips);
204  base::_gs.ownership_by_variant [variant] = base::_geo_element [variant].ownership();
205  }
206  size_type nsid = hdr.dis_size_by_dimension [side_dim];
207  base::_gs.ownership_by_dimension [side_dim] = distributor (nsid, base::comm(), nsid);
208  size_type isid = 0;
209  for (iterator iter = begin(side_dim), last = end(side_dim); iter != last; iter++, isid++) {
210  geo_element& S = *iter;
211  S.set_ios_dis_ie (isid);
212  S.set_dis_ie (isid);
213  }
214  }
215  }
216  vector<index_set> ball [4];
218  while (dom.get (ips, *this, ball)) {
219  base::_domains.push_back (dom);
220  }
221  set_element_side_index (1);
222  set_element_side_index (2);
224  return ips;
225 }
226 template <class T>
227 void
228 geo_rep<T,sequential>::dump (std::string name) const {
229  std::ofstream os ((name + "-dump.geo").c_str());
230  odiststream ods (os, base::_node.ownership().comm());
231  put_geo (ods);
232 }
233 template <class T>
234 void
235 geo_rep<T,sequential>::load (std::string filename, const communicator&)
236 {
237  idiststream ips;
238  ips.open (filename, "geo");
239  check_macro(ips.good(), "\"" << filename << "[.geo[.gz]]\" not found.");
240  get (ips);
241  std::string root_name = delete_suffix (delete_suffix(filename, "gz"), "geo");
242  std::string name = get_basename (root_name);
243  base::_name = name;
244 }
245 template class geo_rep<Float,sequential>;
247 
248 } // namespace rheolef
249