Rheolef  7.1
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&
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 }
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>
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::_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>;
327 template idiststream& geo_basic<Float,sequential>::get (idiststream&);
328 
329 } // namespace rheolef
rheolef::reference_element::last_variant_by_dimension
static variant_type last_variant_by_dimension(size_type dim)
Definition: reference_element.h:150
rheolef::compute_bbox
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
rheolef::geo_basic
generic mesh with rerefence counting
Definition: domain_indirect.h:64
rheolef::geo_element_indirect::set
void set(orientation_type orient, size_type ige, size_type shift=0)
Definition: geo_element_indirect.h:69
rheolef::geo_element::subgeo_size
size_type subgeo_size(size_type subgeo_dim, size_type loc_isid) const
Definition: geo_element.h:221
rheolef::geo_element::get_orientation_and_shift
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
rheolef::delete_suffix
string delete_suffix(const string &name, const string &suffix)
delete_suffix: see the rheostream page for the full documentation
Definition: rheostream.cc:222
rheolef::geo_element::last_inod
size_type last_inod(size_type subgeo_dim) const
Definition: geo_element.h:227
rheolef::domain_indirect_basic< sequential >::get
idiststream & get(idiststream &ips, const geo_rep< T, sequential > &omega, std::vector< index_set > *ball)
Definition: domain_indirect.h:535
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
rheolef::index_set::inplace_intersection
void inplace_intersection(const index_set &b)
Definition: index_set_body.icc:126
rheolef::geo_element::face_indirect
const geo_element_indirect & face_indirect(size_type i) const
Definition: geo_element.h:197
rheolef::geo_header::dimension
size_type dimension
Definition: geo_header.h:39
rheolef::geo_abstract_base_rep::size_type
geo_element_hack::size_type size_type
Definition: geo.h:260
rheolef::geo_get_vtk
idiststream & geo_get_vtk(idiststream &ips, geo_basic< T, sequential > &omega)
Definition: geo_seq_get_vtk.cc:39
rheolef::geo_rep< T, sequential >
Definition: geo.h:786
rheolef::geo_base_rep
base class for M=sequential or distributed meshes representations
Definition: geo.h:528
dump
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
Definition: iorheo-members.h:119
rheolef::get_basename
string get_basename(const string &name)
get_basename: see the rheostream page for the full documentation
Definition: rheostream.cc:254
rheolef::distributor
see the distributor page for the full documentation
Definition: distributor.h:62
mkgeo_ball.order
order
Definition: mkgeo_ball.sh:343
rheolef::geo_element::n_subgeo
size_type n_subgeo(size_type subgeo_dim) const
Definition: geo_element.h:212
rheolef::geo_element::size
size_type size() const
Definition: geo_element.h:168
rheolef::geo_element
see the geo_element page for the full documentation
Definition: geo_element.h:102
rheolef::geo_rep
sequential mesh representation
Definition: domain_indirect.h:63
rheolef::geo_element::shift_type
geo_element_indirect::shift_type shift_type
Definition: geo_element.h:135
rheolef::index_set
Definition: index_set_header.icc:58
rheolef::size_type
size_t size_type
Definition: basis_get.cc:76
rheolef::geo_element::variant
variant_type variant() const
Definition: geo_element.h:161
rheolef::reference_element::first_variant_by_dimension
static variant_type first_variant_by_dimension(size_type dim)
Definition: reference_element.h:148
rheolef::empty_set
static std::set< size_t > empty_set
Definition: space_seq.cc:59
bamg
verbose clean transpose logscale grid shrink ball stereo iso volume skipvtk deformation fastfieldload lattice reader_on_stdin color format format bamg
Definition: iorheo-members.h:71
rheolef::geo_element::set_dis_ie
void set_dis_ie(size_type dis_ie)
Definition: geo_element.h:172
mkgeo_ball.variant
variant
Definition: mkgeo_ball.sh:149
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 vtk
Definition: iorheo-members.h:113
rheolef::geo_element::subgeo_local_vertex
size_type subgeo_local_vertex(size_type subgeo_dim, size_type i_subgeo, size_type i_subgeo_vertex) const
Definition: geo_element.h:223
rheolef::geo_header::map_dimension
size_type map_dimension
Definition: geo_header.h:40
rheolef::geo_iterator
geo iterator
Definition: geo.h:193
rheolef::geo_element::edge_indirect
const geo_element_indirect & edge_indirect(size_type i) const
Definition: geo_element.h:193
rheolef::geo_element::orientation_type
geo_element_indirect::orientation_type orientation_type
Definition: geo_element.h:134
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
rheolef::geo_header::sys_coord
coordinate_type sys_coord
Definition: geo_header.h:41
rheolef::geo_element::parameter_type
Definition: geo_element.h:136
rheolef::odiststream
odiststream: see the diststream page for the full documentation
Definition: diststream.h:126
rheolef::disarray
see the disarray page for the full documentation
Definition: disarray.h:459
mkgeo_ball.d
d
Definition: mkgeo_ball.sh:154
rheolef::geo_element::first_inod
size_type first_inod(size_type subgeo_dim) const
Definition: geo_element.h:225
rheolef::domain_indirect_basic< sequential >
Definition: domain_indirect.h:342
rheolef::geo_element::get_edge_orientation
orientation_type get_edge_orientation(size_type dis_iv0, size_type dis_iv1) const
Definition: geo_element.cc:155
rheolef::reference_element::p
static const variant_type p
Definition: reference_element.h:75
rheolef::space_numbering::nnod
size_type nnod(const basis_basic< T > &b, const geo_size &gs, size_type map_dim)
Definition: space_numbering.cc:54
rheolef::geo_header::dis_size_by_dimension
size_type dis_size_by_dimension[4]
Definition: geo_header.h:44
rheolef::_point_get
point input helper
Definition: geo.h:155
load
void load(idiststream &in, Float &p, field &uh)
Definition: p_laplacian_post.cc:64
rheolef::geo_header::order
size_type order
Definition: geo_header.h:42
rheolef::geo_rep< T, sequential >::get
idiststream & get(idiststream &)
io for geo
Definition: geo_seq_get.cc:127
mkgeo_ball.dim
dim
Definition: mkgeo_ball.sh:307
rheolef::dis_scatch
bool dis_scatch(idiststream &ips, const communicator &comm, std::string ch)
distributed version of scatch(istream&,string)
Definition: diststream.cc:43
rheolef::std
Definition: vec_expr_v2.h:391
rheolef::geo_header
Definition: geo_header.h:32
mkgeo_contraction.name
name
Definition: mkgeo_contraction.sh:133
rheolef::geo_header::dis_size_by_variant
size_type dis_size_by_variant[reference_element::max_variant]
Definition: geo_header.h:43
rheolef::geo_get_bamg
idiststream & geo_get_bamg(idiststream &ips, geo_basic< T, sequential > &omega)
Definition: geo_seq_get_bamg.cc:189
rheolef::geo_element::set_ios_dis_ie
void set_ios_dis_ie(size_type ios_dis_ie)
Definition: geo_element.h:173