Rheolef  7.2
an efficient C++ finite element environment
geo_seq_get.cc
Go to the documentation of this file.
1 #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&
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 }
51 template <class T>
52 void
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  disarray<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::_gs._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)];
99  K.edge_indirect (loc_isid).set (orient, isid);
100  } else { // side_dim == 2
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 // =========================================================================
125 template <class T>
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>
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::_gs._map_dimension = hdr.map_dimension;
161  base::_sys_coord = hdr.sys_coord;
162  base::_piola_basis.reset_family_index (hdr.order);
164  size_type n_edg = hdr.dis_size_by_dimension [1];
165  size_type n_fac = hdr.dis_size_by_dimension [2];
166  size_type n_elt = hdr.dis_size_by_dimension [base::_gs._map_dimension];
167  // ------------------------------------------------------------------------
168  // 2) get coordinates
169  // ------------------------------------------------------------------------
170  base::_node.resize (nnod);
171  if (base::_dimension > 0) {
172  base::_node.get_values (ips, _point_get<T>(geo_base_rep<T,sequential>::_dimension));
173  check_macro (ips.good(), "bad input stream for geo.");
174  }
175  base::_gs.node_ownership = base::_node.ownership();
176  // ------------------------------------------------------------------------
177  // 3) get elements
178  // ------------------------------------------------------------------------
179  if (base::_gs._map_dimension > 0) {
180  for (size_type variant = reference_element::first_variant_by_dimension(base::_gs._map_dimension);
181  variant < reference_element:: last_variant_by_dimension(base::_gs._map_dimension); variant++) {
183  base::_geo_element [variant].resize (hdr.dis_size_by_variant [variant], param);
184  base::_geo_element [variant].get_values (ips);
185  base::_gs.ownership_by_variant [variant] = base::_geo_element [variant].ownership();
186  }
187  base::_gs.ownership_by_dimension [base::_gs._map_dimension] = distributor (n_elt, base::comm(), n_elt);
188  }
189  // ------------------------------------------------------------------------
190  // 4) check that nodes are numbered by increasing node_subgeo_dim
191  // ------------------------------------------------------------------------
192  // ICI va devenir obsolete car les noeuds seront numerotes par _numbering=Pk_numbering
193  {
194  std::vector<size_type> node_subgeo_dim (nnod, size_type(-1));
195  size_type prev_variant = 0;
196  for (iterator iter = begin(base::_gs._map_dimension), last = end(base::_gs._map_dimension); iter != last; iter++) {
197  geo_element& K = *iter;
198  check_macro (prev_variant <= K.variant(), "elements should be numbered by increasing variants (petqTPH)");
199  prev_variant = K.variant();
200  for (size_type d = 0; d <= base::_gs._map_dimension; d++) {
201  for (size_type loc_inod = K.first_inod(d), loc_nnod = K.last_inod(d); loc_inod < loc_nnod; loc_inod++) {
202  node_subgeo_dim [K[loc_inod]] = d;
203  }
204  }
205  }
206  size_type prev_node_dim = 0;
207  for (typename std::vector<size_type>::const_iterator iter = node_subgeo_dim.begin(), last = node_subgeo_dim.end();
208  iter != last; iter++) {
209  check_macro (prev_node_dim <= *iter, "nodes should be numbered by increasing subgeo dimension");
210  prev_node_dim = *iter;
211  }
212  }
213  // ------------------------------------------------------------------------
214  // 5) compute n_vert (n_vert < nnod when order > 1) and set element indexes (K.dis_ie & K.ios_dis_ie)
215  // ------------------------------------------------------------------------
216  size_type n_vert = 0;
217  if (base::_gs._map_dimension == 0) {
218  n_vert = nnod;
219  } else {
220  std::vector<size_t> is_vertex (nnod, 0);
221  size_type ie = 0;
222  for (iterator iter = begin(base::_gs._map_dimension), last = end(base::_gs._map_dimension); iter != last; iter++, ie++) {
223  geo_element& K = *iter;
224  K.set_ios_dis_ie (ie);
225  K.set_dis_ie (ie);
226  if (base::order() > 1) {
227  for (size_type iloc = 0, nloc = K.size(); iloc < nloc; iloc++) {
228  is_vertex [K[iloc]] = 1;
229  }
230  }
231  }
232  if (base::order() == 1) {
233  n_vert = nnod;
234  } else {
235  n_vert = accumulate (is_vertex.begin(), is_vertex.end(), 0);
236  }
237  }
238  // ------------------------------------------------------------------------
239  // 6) create vertex-element (0d elements)
240  // ------------------------------------------------------------------------
242  base::_geo_element [reference_element::p].resize (n_vert, param);
243  size_type first_iv = base::_node.ownership().first_index();
244  {
245  size_type iv = 0;
246  for (iterator iter = begin(0), last = end(0); iter != last; iter++, iv++) {
247  geo_element& P = *iter;
248  P[0] = first_iv + iv;
249  P.set_dis_ie (first_iv + iv); // TODO: P[0] & dis_ie redundant for `p'
250  P.set_ios_dis_ie (first_iv + iv);
251  }
252  }
253  // ownership_by_dimension[0]: used by connectivity
254  base::_gs.ownership_by_variant [reference_element::p] = base::_geo_element [reference_element::p].ownership();
255  base::_gs.ownership_by_dimension [0] = base::_geo_element [reference_element::p].ownership();
256  // ------------------------------------------------------------------------
257  // 7) get faces & edge
258  // ------------------------------------------------------------------------
259  if (base::_gs._map_dimension > 0) {
260 
261  for (size_type side_dim = base::_gs._map_dimension - 1; side_dim >= 1; side_dim--) {
265  base::_geo_element [variant].resize (hdr.dis_size_by_variant [variant], param);
266  base::_geo_element [variant].get_values (ips);
267  base::_gs.ownership_by_variant [variant] = base::_geo_element [variant].ownership();
268  }
269  size_type nsid = hdr.dis_size_by_dimension [side_dim];
270  base::_gs.ownership_by_dimension [side_dim] = distributor (nsid, base::comm(), nsid);
271  size_type isid = 0;
272  for (iterator iter = begin(side_dim), last = end(side_dim); iter != last; iter++, isid++) {
273  geo_element& S = *iter;
274  S.set_ios_dis_ie (isid);
275  S.set_dis_ie (isid);
276  }
277  }
278  }
279  // ------------------------------------------------------------------------
280  // 8) get domain, until end-of-file
281  // ------------------------------------------------------------------------
282  vector<index_set> ball [4];
284  while (dom.get (ips, *this, ball)) {
285  base::_domains.push_back (dom);
286  }
287  // ------------------------------------------------------------------------
288  // 9) set indexes on faces and edges of elements, for P2 approx
289  // ------------------------------------------------------------------------
290  set_element_side_index (1);
291  set_element_side_index (2);
292  // ------------------------------------------------------------------------
293  // 10) bounding box: _xmin, _xmax
294  // ------------------------------------------------------------------------
296  return ips;
297 }
298 // ----------------------------------------------------------------------------
299 // dump
300 // ----------------------------------------------------------------------------
301 template <class T>
302 void
303 geo_rep<T,sequential>::dump (std::string name) const {
304  std::ofstream os ((name + "-dump.geo").c_str());
305  odiststream ods (os, base::_node.ownership().comm());
306  put_geo (ods);
307 }
308 // ----------------------------------------------------------------------------
309 // read from file
310 // ----------------------------------------------------------------------------
311 template <class T>
312 void
313 geo_rep<T,sequential>::load (std::string filename, const communicator&)
314 {
315  idiststream ips;
316  ips.open (filename, "geo");
317  check_macro(ips.good(), "\"" << filename << "[.geo[.gz]]\" not found.");
318  get (ips);
319  std::string root_name = delete_suffix (delete_suffix(filename, "gz"), "geo");
320  std::string name = get_basename (root_name);
321  base::_name = name;
322 }
323 // ----------------------------------------------------------------------------
324 // instanciation in library
325 // ----------------------------------------------------------------------------
326 template class geo_rep<Float,sequential>;
328 
329 } // namespace rheolef
see the disarray page for the full documentation
Definition: disarray.h:497
see the distributor page for the full documentation
Definition: distributor.h:69
idiststream & get(idiststream &ips, const geo_rep< T, sequential > &omega, std::vector< index_set > *ball)
geo_element_hack::size_type size_type
Definition: geo.h:260
base class for M=sequential or distributed meshes representations
Definition: geo.h:528
generic mesh with rerefence counting
Definition: geo.h:1089
void set(orientation_type orient, size_type ige, size_type shift=0)
see the geo_element page for the full documentation
Definition: geo_element.h:102
geo_element_indirect::orientation_type orientation_type
Definition: geo_element.h:134
bool get_orientation_and_shift(const geo_element &S, orientation_type &orient, shift_type &shift) const
return orientation and shift between *this element and S
Definition: geo_element.cc:114
size_type subgeo_local_vertex(size_type subgeo_dim, size_type i_subgeo, size_type i_subgeo_vertex) const
Definition: geo_element.h:223
size_type first_inod(size_type subgeo_dim) const
Definition: geo_element.h:225
size_type size() const
Definition: geo_element.h:168
void set_ios_dis_ie(size_type ios_dis_ie)
Definition: geo_element.h:173
size_type subgeo_size(size_type subgeo_dim, size_type loc_isid) const
Definition: geo_element.h:221
variant_type variant() const
Definition: geo_element.h:161
orientation_type get_edge_orientation(size_type dis_iv0, size_type dis_iv1) const
Definition: geo_element.cc:155
geo_element_indirect::shift_type shift_type
Definition: geo_element.h:135
size_type n_subgeo(size_type subgeo_dim) const
Definition: geo_element.h:212
const geo_element_indirect & face_indirect(size_type i) const
Definition: geo_element.h:197
const geo_element_indirect & edge_indirect(size_type i) const
Definition: geo_element.h:193
size_type last_inod(size_type subgeo_dim) const
Definition: geo_element.h:227
void set_dis_ie(size_type dis_ie)
Definition: geo_element.h:172
idiststream & get(idiststream &)
io for geo
Definition: geo_seq_get.cc:127
sequential mesh representation
Definition: geo.h:778
idiststream: see the diststream page for the full documentation
Definition: diststream.h:336
std::istream & is()
Definition: diststream.h:400
void open(std::string filename, std::string suffix="", const communicator &comm=communicator())
This routine opens a physical input file.
Definition: diststream.cc:85
bool good() const
Definition: diststream.cc:124
const communicator & comm() const
Definition: diststream.h:356
void inplace_intersection(const index_set &b)
odiststream: see the diststream page for the full documentation
Definition: diststream.h:137
static variant_type last_variant_by_dimension(size_type dim)
static variant_type first_variant_by_dimension(size_type dim)
static const variant_type p
static std::set< size_t > empty_set
Definition: space_seq.cc:59
size_t size_type
Definition: basis_get.cc:76
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
verbose clean transpose logscale grid shrink ball stereo iso volume skipvtk deformation fastfieldload lattice reader_on_stdin color format format bamg
verbose clean transpose logscale grid shrink ball stereo iso volume skipvtk deformation fastfieldload lattice reader_on_stdin color format format format format format format format format format format format format format format format format vtk
verbose clean transpose logscale grid shrink ball stereo iso volume skipvtk deformation fastfieldload lattice reader_on_stdin color format format format format format format format format format format format format format format format format format format dump
size_type nnod(const basis_basic< T > &b, const geo_size &gs, size_type map_dim)
This file is part of Rheolef.
void compute_bbox(const geo_base_rep< T, M > &omega, const geo_element &K, point_basic< T > &xmin, point_basic< T > &xmax)
Definition: geo_locate.cc:50
string delete_suffix(const string &name, const string &suffix)
delete_suffix: see the rheostream page for the full documentation
Definition: rheostream.cc:226
string get_basename(const string &name)
get_basename: see the rheostream page for the full documentation
Definition: rheostream.cc:258
bool dis_scatch(idiststream &ips, const communicator &comm, std::string ch)
distributed version of scatch(istream&,string)
Definition: diststream.cc:44
idiststream & geo_get_vtk(idiststream &ips, geo_basic< T, sequential > &omega)
idiststream & geo_get_bamg(idiststream &ips, geo_basic< T, sequential > &omega)
void load(idiststream &in, Float &p, field &uh)
point input helper
Definition: geo.h:155
bool need_upgrade() const
Definition: geo_header.cc:79
size_type map_dimension
Definition: geo_header.h:40
size_type dis_size_by_dimension[4]
Definition: geo_header.h:44
coordinate_type sys_coord
Definition: geo_header.h:41
size_type dimension
Definition: geo_header.h:39
size_type dis_size_by_variant[reference_element::max_variant]
Definition: geo_header.h:43
geo iterator
Definition: geo.h:193