Rheolef  7.1
an efficient C++ finite element environment
geo_seq_get_bamg.cc
Go to the documentation of this file.
1 //
22 // geo: bamg input
23 //
24 // author: Pierre.Saramito@imag.fr
25 //
26 // 5 march 2012
27 //
28 #include "rheolef/geo.h"
29 #include "rheolef/rheostream.h"
30 #include "rheolef/iorheo.h"
31 using namespace std;
32 namespace rheolef {
33 
34 static
35 void
36 bamg_load_element (std::istream& is, size_t variant, geo_element_auto<>& K)
37 {
38  K.reset (variant, 1);
39  for (size_t iloc = 0, nloc = K.size(); iloc < nloc; iloc++) {
40  is >> K[iloc];
41  K[iloc]--;
42  }
43 }
44 template <class T>
45 static
46 void
47 build_domains (
48  size_t dom_dim,
49  const disarray<geo_element_auto<>, sequential>& elt,
50  const disarray<size_t, sequential>& elt_bamg_dom_id,
51  const std::vector<std::string>& dom_name,
52  const geo_basic<T,sequential>& omega,
53  vector<index_set>* ball)
54 {
55  using namespace std;
57  typedef pair<size_type,size_type> pair_t;
58  typedef geo_element_auto<> element_t;
59  typedef disarray<element_t, sequential> e_array_t;
60  typedef disarray<size_type, sequential> i_array_t;
61  typedef map<size_type, pair_t, less<size_type> > map_t;
62 
63  // -------------------------------------------------------------------
64  // 1) reduce the bamg domain id to [0:ndom[ by counting domain id
65  // -------------------------------------------------------------------
66  // map[bamg_id] will gives idom and the number of its elements
67  map_t bamg_id2idom; // TODO (less<size_type>());
68  size_type idom = 0;
69  for (size_type ie = 0, ne = elt_bamg_dom_id.size(); ie < ne; ie++) {
70  size_type bamg_id = elt_bamg_dom_id [ie];
71  typename map_t::iterator iter = bamg_id2idom.find (bamg_id);
72  if (iter != bamg_id2idom.end()) {
73  // here is a previous bamg_dom_id: increment counter
74  ((*iter).second.second)++;
75  continue;
76  }
77  // here is a new bamg_dom_id: associated to idom and elt counter=1
78  bamg_id2idom.insert (pair<size_type,pair_t>(bamg_id, pair_t(idom,1)));
79  idom++;
80  }
81  size_type ndom = bamg_id2idom.size();
82  // -------------------------------------------------------------------
83  // 2) check that ndom matches the domain name disarray size
84  // -------------------------------------------------------------------
85  if (ndom != dom_name.size()) {
86  warning_macro ("geo::get_bamg: "
87  << dom_name.size() << " domain name(s) founded while "
88  << ndom << " bamg "<<dom_dim<<"d domain(s) are defined");
89  cerr << endl;
90  error_macro("HINT: check domain name file (.dmn)");
91  }
92  // -------------------------------------------------------------------
93  // 3) create domain disarray: loop on elements
94  // -------------------------------------------------------------------
95  element_t dummy_elt (reference_element::p, 0);
96  vector<e_array_t> domain (ndom, e_array_t(0, dummy_elt));
97  vector<size_type> counter (ndom, 0);
98  for (size_type ie = 0, ne = elt_bamg_dom_id.size(); ie < ne; ie++) {
99  size_type bamg_dom_id = elt_bamg_dom_id [ie];
100  size_type idom = bamg_id2idom [bamg_dom_id].first;
101  if (domain[idom].size() == 0) {
102  size_type size = bamg_id2idom [bamg_dom_id].second;
103  domain[idom] = e_array_t (size, dummy_elt);
104  }
105  size_type dom_ie = counter[idom];
106  domain[idom][dom_ie] = elt[ie];
107  counter[idom]++;
108  }
109  // -------------------------------------------------------------------
110  // 3) insert domains in the mesh: loop on domains
111  // -------------------------------------------------------------------
112  for (typename map_t::const_iterator
113  iter = bamg_id2idom.begin(),
114  last = bamg_id2idom.end(); iter != last; ++iter) {
115  size_type bamg_id = (*iter).first;
116  size_type idom = (*iter).second.first;
117  size_type size = (*iter).second.second;
118  check_macro (idom < dom_name.size(), "invalid idom="<<idom<<" for domain name");
119  string name = dom_name [idom];
120 
121  domain_indirect_basic<sequential> dom (domain[idom], omega, ball);
122  dom.set_name (name);
123  dom.set_map_dimension (dom_dim);
124  omega.insert_domain_indirect (dom);
125  }
126  size_type ndom2 = omega.n_domain_indirect();
127 }
128 template <class T>
129 static
130 void
131 build_vertex_domains (
132  const disarray<size_t, sequential>& edg_bdr_bamg_dom_id,
133  const disarray<size_t, sequential>& vert_bamg_dom_id,
134  const std::vector<std::string>& dom_name,
135  const geo_basic<T,sequential>& omega,
136  vector<index_set>* ball)
137 {
138  if (dom_name.size() == 0) return;
139  typedef geo_element_auto<> element_t;
140  typedef disarray<element_t, sequential> v_array_t;
141  element_t dummy_elt (reference_element::p, 0);
142  // -------------------------------------------------------------------
143  // 1) list all vertex domain id
144  // -------------------------------------------------------------------
145  std::set<size_t> vert_id;
146  for (size_t iv = 0, nv = vert_bamg_dom_id.size(); iv < nv; ++iv) {
147  size_t dom_id = vert_bamg_dom_id [iv];
148  if (dom_id == 0) continue;
149  vert_id.insert (dom_id);
150  }
151  // -------------------------------------------------------------------
152  // 2) omit vertex that are marked with an edge id
153  // -------------------------------------------------------------------
154  for (size_t iedg_bdr = 0, nedg_bdr = edg_bdr_bamg_dom_id.size(); iedg_bdr < nedg_bdr; ++iedg_bdr) {
155  size_t dom_id = edg_bdr_bamg_dom_id [iedg_bdr];
156  vert_id.erase (dom_id);
157  }
158  check_macro (vert_id.size() == dom_name.size(),
159  "unexpected VertexDomainNames with "<<dom_name.size()
160  <<" domain names while the mesh provides " << vert_id.size()
161  <<" vertex labels");
162  // -------------------------------------------------------------------
163  // 3) loop on vertex domain and insert it in mesh
164  // -------------------------------------------------------------------
165  size_t idom = 0;
166  for (std::set<size_t>::const_iterator iter = vert_id.begin(); iter != vert_id.end(); ++iter, ++idom) {
167  size_t id = *iter;
168  string name = dom_name[idom];
169  size_t dom_size = 0;
170  for (size_t iv = 0, nv = vert_bamg_dom_id.size(); iv < nv; ++iv) {
171  if (vert_bamg_dom_id [iv] == id) dom_size++;
172  }
173  v_array_t vert_list (dom_size, dummy_elt);
174  for (size_t iv = 0, iv_dom = 0, nv = vert_bamg_dom_id.size(); iv < nv; ++iv) {
175  if (vert_bamg_dom_id [iv] != id) continue;
176  element_t& K = vert_list [iv_dom];
177  K.reset (reference_element::p, 1);
178  K[0] = iv;
179  iv_dom++;
180  }
181  domain_indirect_basic<sequential> dom (vert_list, omega, ball);
182  dom.set_name (name);
183  dom.set_map_dimension (0);
184  omega.insert_domain_indirect (dom);
185  }
186 }
187 template <class T>
188 idiststream&
189 geo_get_bamg (idiststream& ips, geo_basic<T,sequential>& omega)
190 {
191  using namespace std;
193  typedef typename geo_basic<T,sequential>::node_type node_type;
194 
195  check_macro (ips.good(), "bad input stream for bamg");
196  istream& is = ips.is();
197  geo_header hdr;
198  hdr.dimension = 2;
199  hdr.map_dimension = 2;
200  hdr.order = 1;
201  hdr.dis_size_by_dimension [2] = 0;
202 
203  // ------------------------------------------------------------------------
204  // 1) load data
205  // ------------------------------------------------------------------------
206  typedef geo_element_auto<> element_t;
207  typedef disarray<node_type, sequential> n_array_t;
208  typedef disarray<element_t, sequential> e_array_t;
209  typedef disarray<size_type, sequential> i_array_t;
210 
211  n_array_t vertex;
212  i_array_t vert_bamg_dom_id;
213  i_array_t edg_bdr_bamg_dom_id;
214  i_array_t tri_bamg_dom_id;
215  i_array_t qua_bamg_dom_id;
216  e_array_t edg_bdr;
217  std::array<e_array_t, reference_element::max_variant> tmp_geo_element;
218  string label;
219  while (is.good() && label != "End") {
220  is >> label;
221  if (label == "Vertices") {
222  // ------------------------------------------------------------------------
223  // 1.1) load the coordinates
224  // Vertices <np>
225  // {xi yi dom_idx} i=0..np-1
226  // ------------------------------------------------------------------------
227  size_type nnod;
228  is >> nnod;
229  if (nnod == 0) {
230  warning_macro("empty bamg mesh file");
231  return ips;
232  }
233  hdr.dis_size_by_dimension [0] = nnod;
234  hdr.dis_size_by_variant [0] = nnod;
235  vertex.resize (nnod);
236  vert_bamg_dom_id = i_array_t (nnod, 0);
237  for (size_type inod = 0; inod < nnod; inod++) {
238  is >> vertex[inod][0] >> vertex[inod][1] >> vert_bamg_dom_id[inod];
239  }
240  check_macro (ips.good(), "bad input stream for bamg");
241  } else if (label == "Triangles") {
242  // ------------------------------------------------------------------------
243  // 2.1) load the connectivity
244  // Triangle <nt>
245  // {s0 s1 s2 dom_idx} j=0..nt-1
246  // ------------------------------------------------------------------------
247  size_type nt;
248  is >> nt;
249  hdr.dis_size_by_dimension [2] += nt;
250  hdr.dis_size_by_variant [reference_element::t] = nt;
251  element_t init_tri (reference_element::t, hdr.order);
252  tmp_geo_element [reference_element::t] = e_array_t (nt, init_tri);
253  tri_bamg_dom_id = i_array_t (nt, 0);
254  for (size_type it = 0; it < nt; it++) {
255  bamg_load_element (is, reference_element::t,
256  tmp_geo_element[reference_element::t][it]);
257  is >> tri_bamg_dom_id[it];
258  }
259  } else if (label == "Quadrilaterals") {
260  // ------------------------------------------------------------------------
261  // Quadrilaterals <nq>
262  // {s0 s1 s2 s3 dom_idx} j=0..nq-1
263  // ------------------------------------------------------------------------
264  size_type nq;
265  is >> nq;
266  hdr.dis_size_by_dimension [2] += nq;
267  hdr.dis_size_by_variant [reference_element::q] = nq;
268  element_t init_qua (reference_element::q, hdr.order);
269  tmp_geo_element [reference_element::q] = e_array_t (nq, init_qua);
270  qua_bamg_dom_id = i_array_t (nq, 0);
271  for (size_type iq = 0; iq < nq; iq++) {
272  bamg_load_element (is, reference_element::q,
273  tmp_geo_element[reference_element::q][iq]);
274  is >> qua_bamg_dom_id[iq];
275  }
276  } else if (label == "Edges") {
277  // ------------------------------------------------------------------------
278  // 2.3) load the boundary domains
279  // Edges <nedg>
280  // {s0 s1 dom_idx} j=0..nedg-1
281  // ------------------------------------------------------------------------
282  size_type nedg;
283  is >> nedg;
284  element_t init_edg (reference_element::e, hdr.order);
285  edg_bdr = e_array_t (nedg, init_edg);
286  edg_bdr_bamg_dom_id = i_array_t (nedg, 0);
287  for (size_type iedg = 0; iedg < nedg; iedg++) {
288  element_t& K = edg_bdr[iedg];
289  K.reset (reference_element::e, hdr.order);
290  for (size_type iloc = 0, nloc = 2; iloc < nloc; iloc++) {
291  is >> K[iloc];
292  K[iloc]--;
293  }
294  is >> edg_bdr_bamg_dom_id[iedg];
295  }
296  }
297  }
298  // ---------------------------------------------------------------
299  // 2) check rheolef extension to optional domain names
300  // ---------------------------------------------------------------
301  vector<string> vertice_domain_name;
302  vector<string> edg_dom_name;
303  vector<string> region_domain_name;
304  char c;
305  is >> ws >> c; // skip white and grab a char
306  // have "EdgeDomainNames" or "VertexDomainNames" ?
307  // bamg mesh may be followed by field data and such, so be carrefull...
308  while (c == 'E' || c == 'V' || c == 'R') {
309  is.unget(); // put char back
310  if (c == 'R') {
311  if (!scatch(is,"RegionDomainNames")) break;
312  size_type n_dom_region;
313  is >> n_dom_region;
314  region_domain_name.resize (n_dom_region);
315  for (size_type k = 0; k < n_dom_region; k++) {
316  is >> region_domain_name[k];
317  }
318  } else if (c == 'E') {
319  if (!scatch(is,"EdgeDomainNames")) break;
320  // ---------------------------------------------------------------
321  // get edge domain names
322  // ---------------------------------------------------------------
323  size_type n_dom_edge;
324  is >> n_dom_edge;
325  edg_dom_name.resize (n_dom_edge);
326  for (size_type k = 0; k < n_dom_edge; k++) {
327  is >> edg_dom_name[k];
328  }
329  } else if (c == 'V') {
330  if (!scatch(is,"VertexDomainNames")) break;
331  // ---------------------------------------------------------------
332  // get vertice domain names
333  // ---------------------------------------------------------------
334  size_type n_dom_vertice;
335  is >> n_dom_vertice;
336  vertice_domain_name.resize (n_dom_vertice);
337  for (size_type k = 0; k < n_dom_vertice; k++) {
338  is >> vertice_domain_name[k];
339  }
340  }
341  is >> ws >> c; // skip white and grab a char
342  }
343  // ------------------------------------------------------------------------
344  // 3) build & upgrade
345  // ------------------------------------------------------------------------
346  bool do_upgrade = true;
347  omega.build_from_data (hdr, vertex, tmp_geo_element, do_upgrade);
348 
349  // ------------------------------------------------------------------------
350  // 4) get domain, until end-of-file
351  // ------------------------------------------------------------------------
352  vector<index_set> ball[4];
353  build_domains (1, edg_bdr, edg_bdr_bamg_dom_id, edg_dom_name, omega, ball);
354  build_vertex_domains (edg_bdr_bamg_dom_id, vert_bamg_dom_id, vertice_domain_name, omega, ball);
355  //TODO: region(d=2) domains
356  return ips;
357 }
358 // ----------------------------------------------------------------------------
359 // instanciation in library
360 // ----------------------------------------------------------------------------
361 template idiststream& geo_get_bamg<Float> (idiststream&, geo_basic<Float,sequential>&);
362 
363 }// namespace rheolef
field::size_type size_type
Definition: branch.cc:425
see the disarray page for the full documentation
Definition: disarray.h:459
void build_from_data(const geo_header &hdr, const disarray< node_type, sequential > &node, std::array< disarray< geo_element_auto<>, sequential >, reference_element::max_variant > &tmp_geo_element, bool do_upgrade)
generic mesh with rerefence counting
Definition: geo.h:1089
#define error_macro(message)
Definition: dis_macros.h:49
#define warning_macro(message)
Definition: dis_macros.h:53
const point vertex[n_vertex]
Definition: edge.icc:67
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
size_type nnod(const basis_basic< T > &b, const geo_size &gs, size_type map_dim)
This file is part of Rheolef.
template idiststream & geo_get_bamg< Float >(idiststream &, geo_basic< Float, sequential > &)
bool scatch(std::istream &in, const std::string &ch, bool full_match=true)
scatch: see the rheostream page for the full documentation
Definition: scatch.icc:52
idiststream & geo_get_bamg(idiststream &ips, geo_basic< T, sequential > &omega)
size_type map_dimension
Definition: geo_header.h:40
size_type dis_size_by_dimension[4]
Definition: geo_header.h:44
size_type dimension
Definition: geo_header.h:39
size_type dis_size_by_variant[reference_element::max_variant]
Definition: geo_header.h:43