rheolef  6.5
geo_seq_get_bamg.cc
Go to the documentation of this file.
1 #include "rheolef/geo.h"
2 #include "rheolef/rheostream.h"
3 #include "rheolef/iorheo.h"
4 using namespace std;
5 namespace rheolef {
6 
7 template<class A>
8 static
9 void
10 bamg_load_element (std::istream& is, size_t variant, geo_element_auto<A>& K)
11 {
12  K.reset (variant, 1);
13  for (size_t iloc = 0, nloc = K.size(); iloc < nloc; iloc++) {
14  is >> K[iloc];
15  K[iloc]--;
16  }
17 }
18 template <class T, class A>
19 static
20 void
22  size_t dom_dim,
23  const array<geo_element_auto<A>, sequential, A>& elt,
24  const array<size_t, sequential, A>& elt_bamg_dom_id,
25  const std::vector<std::string>& dom_name,
26  const geo_basic<T,sequential>& omega,
27  vector<index_set>* ball)
28 {
29  using namespace std;
31  typedef pair<size_type,size_type> pair_t;
32  typedef geo_element_auto<A> element_t;
33  typedef array<element_t, sequential, A> e_array_t;
34  typedef array<size_type, sequential, A> i_array_t;
35  typedef map<size_type, pair_t, less<size_type>, A> map_t;
36 
37  A alloc = elt.get_allocator();
38  map_t bamg_id2idom (less<size_type>(), alloc);
39  size_type idom = 0;
40  for (size_type ie = 0, ne = elt_bamg_dom_id.size(); ie < ne; ie++) {
41  size_type bamg_id = elt_bamg_dom_id [ie];
42  typename map_t::iterator iter = bamg_id2idom.find (bamg_id);
43  if (iter != bamg_id2idom.end()) {
44  ((*iter).second.second)++;
45  continue;
46  }
47  bamg_id2idom.insert (pair<size_type,pair_t>(bamg_id, pair_t(idom,1)));
48  idom++;
49  }
50  size_type ndom = bamg_id2idom.size();
51  if (ndom != dom_name.size()) {
52  warning_macro ("geo::get_bamg: "
53  << dom_name.size() << " domain name(s) founded while "
54  << ndom << " bamg "<<dom_dim<<"d domain(s) are defined");
55  cerr << endl;
56  error_macro("HINT: check domain name file (.dmn)");
57  }
58  element_t dummy_elt (reference_element::p, 0, alloc);
59  vector<e_array_t, A> domain (ndom, e_array_t(0, dummy_elt, alloc), alloc);
60  vector<size_type, A> counter (ndom, 0, alloc);
61  for (size_type ie = 0, ne = elt_bamg_dom_id.size(); ie < ne; ie++) {
62  size_type bamg_dom_id = elt_bamg_dom_id [ie];
63  size_type idom = bamg_id2idom [bamg_dom_id].first;
64  if (domain[idom].size() == 0) {
65  size_type size = bamg_id2idom [bamg_dom_id].second;
66  domain[idom] = e_array_t (size, dummy_elt, alloc);
67  }
68  size_type dom_ie = counter[idom];
69  domain[idom][dom_ie] = elt[ie];
70  counter[idom]++;
71  }
72  for (typename map_t::const_iterator
73  iter = bamg_id2idom.begin(),
74  last = bamg_id2idom.end(); iter != last; ++iter) {
75  size_type bamg_id = (*iter).first;
76  size_type idom = (*iter).second.first;
77  size_type size = (*iter).second.second;
78  check_macro (idom < dom_name.size(), "invalid idom="<<idom<<" for domain name");
79  string name = dom_name [idom];
80 
81  domain_indirect_basic<sequential> dom (domain[idom], omega, ball);
82  dom.set_name (name);
83  dom.set_map_dimension (dom_dim);
84  omega.insert_domain_indirect (dom);
85  }
86  size_type ndom2 = omega.n_domain_indirect();
87 }
88 template <class T, class A>
89 static
90 void
92  const array<size_t, sequential, A>& edg_bdr_bamg_dom_id,
93  const array<size_t, sequential, A>& vert_bamg_dom_id,
94  const std::vector<std::string>& dom_name,
95  const geo_basic<T,sequential>& omega,
96  vector<index_set>* ball)
97 {
98  if (dom_name.size() == 0) return;
99  typedef geo_element_auto<A> element_t;
100  typedef array<element_t, sequential, A> v_array_t;
101  A alloc = edg_bdr_bamg_dom_id.get_allocator();
102  element_t dummy_elt (reference_element::p, 0, alloc);
103  std::set<size_t> vert_id;
104  for (size_t iv = 0, nv = vert_bamg_dom_id.size(); iv < nv; ++iv) {
105  size_t dom_id = vert_bamg_dom_id [iv];
106  if (dom_id == 0) continue;
107  vert_id.insert (dom_id);
108  }
109  for (size_t iedg_bdr = 0, nedg_bdr = edg_bdr_bamg_dom_id.size(); iedg_bdr < nedg_bdr; ++iedg_bdr) {
110  size_t dom_id = edg_bdr_bamg_dom_id [iedg_bdr];
111  vert_id.erase (dom_id);
112  }
113  check_macro (vert_id.size() == dom_name.size(),
114  "unexpected VerticeDomainNames with "<<dom_name.size()
115  <<" domain names while the mesh provides " << vert_id.size()
116  <<" vertex labels");
117  size_t idom = 0;
118  for (std::set<size_t>::const_iterator iter = vert_id.begin(); iter != vert_id.end(); ++iter, ++idom) {
119  size_t id = *iter;
120  string name = dom_name[idom];
121  size_t dom_size = 0;
122  for (size_t iv = 0, nv = vert_bamg_dom_id.size(); iv < nv; ++iv) {
123  if (vert_bamg_dom_id [iv] == id) dom_size++;
124  }
125  v_array_t vert_list (dom_size, dummy_elt, alloc);
126  for (size_t iv = 0, iv_dom = 0, nv = vert_bamg_dom_id.size(); iv < nv; ++iv) {
127  if (vert_bamg_dom_id [iv] != id) continue;
128  element_t& K = vert_list [iv_dom];
129  K.reset (reference_element::p, 1);
130  K[0] = iv;
131  iv_dom++;
132  }
133  domain_indirect_basic<sequential> dom (vert_list, omega, ball);
134  dom.set_name (name);
135  dom.set_map_dimension (0);
136  omega.insert_domain_indirect (dom);
137  }
138 }
139 template <class T>
142 {
143  using namespace std;
145  typedef typename geo_basic<T,sequential>::node_type node_type;
146 
147  check_macro (ips.good(), "bad input stream for bamg");
148  istream& is = ips.is();
149  geo_header hdr;
150  hdr.dimension = 2;
151  hdr.map_dimension = 2;
152  hdr.order = 1;
153  hdr.dis_size_by_dimension [2] = 0;
154 
155  typedef heap_allocator<size_type> alloc_t;
156  typedef geo_element_auto<alloc_t> element_t;
157  typedef array<node_type, sequential> n_array_t;
158  typedef array<element_t, sequential, alloc_t> e_array_t;
159  typedef array<size_type, sequential, alloc_t> i_array_t;
160 
161  alloc_t alloc;
162  n_array_t vertex;
163  i_array_t vert_bamg_dom_id;
164  i_array_t edg_bdr_bamg_dom_id;
165  i_array_t tri_bamg_dom_id;
166  i_array_t qua_bamg_dom_id;
167  e_array_t edg_bdr;
168  boost::array<e_array_t, reference_element::max_variant> tmp_geo_element;
169  string label;
170  while (is.good() && label != "End") {
171  is >> label;
172  if (label == "Vertices") {
173  // Vertices <np>
174  size_type nnod;
175  is >> nnod;
176  if (nnod == 0) {
177  warning_macro("empty bamg mesh file");
178  return ips;
179  }
180  hdr.dis_size_by_dimension [0] = nnod;
181  hdr.dis_size_by_variant [0] = nnod;
182  vertex.resize (nnod);
183  vert_bamg_dom_id = i_array_t (nnod, 0, alloc);
184  for (size_type inod = 0; inod < nnod; inod++) {
185  is >> vertex[inod][0] >> vertex[inod][1] >> vert_bamg_dom_id[inod];
186  }
187  check_macro (ips.good(), "bad input stream for bamg");
188  } else if (label == "Triangles") {
189  // Triangle <nt>
190  size_type nt;
191  is >> nt;
192  hdr.dis_size_by_dimension [2] += nt;
193  hdr.dis_size_by_variant [reference_element::t] = nt;
194  element_t init_tri (reference_element::t, hdr.order, alloc);
195  tmp_geo_element [reference_element::t] = e_array_t (nt, init_tri, alloc);
196  tri_bamg_dom_id = i_array_t (nt, 0, alloc);
197  for (size_type it = 0; it < nt; it++) {
198  bamg_load_element (is, reference_element::t,
199  tmp_geo_element[reference_element::t][it]);
200  is >> tri_bamg_dom_id[it];
201  }
202  } else if (label == "Quadrilaterals") {
203  // Quadrilaterals <nq>
204  size_type nq;
205  is >> nq;
206  hdr.dis_size_by_dimension [2] += nq;
207  hdr.dis_size_by_variant [reference_element::q] = nq;
208  element_t init_qua (reference_element::q, hdr.order, alloc);
209  tmp_geo_element [reference_element::q] = e_array_t (nq, init_qua, alloc);
210  qua_bamg_dom_id = i_array_t (nq, 0, alloc);
211  for (size_type iq = 0; iq < nq; iq++) {
212  bamg_load_element (is, reference_element::q,
213  tmp_geo_element[reference_element::q][iq]);
214  is >> qua_bamg_dom_id[iq];
215  }
216  } else if (label == "Edges") {
217  // Edges <nedg>
218  size_type nedg;
219  is >> nedg;
220  element_t init_edg (reference_element::e, hdr.order, alloc);
221  edg_bdr = e_array_t (nedg, init_edg, alloc);
222  edg_bdr_bamg_dom_id = i_array_t (nedg, 0, alloc);
223  for (size_type iedg = 0; iedg < nedg; iedg++) {
224  element_t& K = edg_bdr[iedg];
225  K.reset (reference_element::e, hdr.order);
226  for (size_type iloc = 0, nloc = 2; iloc < nloc; iloc++) {
227  is >> K[iloc];
228  K[iloc]--;
229  }
230  is >> edg_bdr_bamg_dom_id[iedg];
231  }
232  }
233  }
234  vector<string> vertice_domain_name;
235  vector<string> edg_dom_name;
236  vector<string> region_domain_name;
237  char c;
238  is >> ws >> c; // skip white and grab a char
239  while (c == 'E' || c == 'V' || c == 'R') {
240  is.unget(); // put char back
241  if (c == 'R') {
242  if (!scatch(is,"RegionDomainNames")) break;
243  size_type n_dom_region;
244  is >> n_dom_region;
245  region_domain_name.resize (n_dom_region);
246  for (size_type k = 0; k < n_dom_region; k++) {
247  is >> region_domain_name[k];
248  }
249  } else if (c == 'E') {
250  if (!scatch(is,"EdgeDomainNames")) break;
251  size_type n_dom_edge;
252  is >> n_dom_edge;
253  edg_dom_name.resize (n_dom_edge);
254  for (size_type k = 0; k < n_dom_edge; k++) {
255  is >> edg_dom_name[k];
256  }
257  } else if (c == 'V') {
258  if (!scatch(is,"VerticeDomainNames")) break;
259  size_type n_dom_vertice;
260  is >> n_dom_vertice;
261  vertice_domain_name.resize (n_dom_vertice);
262  for (size_type k = 0; k < n_dom_vertice; k++) {
263  is >> vertice_domain_name[k];
264  }
265  }
266  is >> ws >> c; // skip white and grab a char
267  }
268  bool do_upgrade = true;
269  omega.build_from_data (hdr, vertex, tmp_geo_element, do_upgrade);
270 
271  vector<index_set> ball[4];
272  build_domains (1, edg_bdr, edg_bdr_bamg_dom_id, edg_dom_name, omega, ball);
273  build_vertex_domains (edg_bdr_bamg_dom_id, vert_bamg_dom_id, vertice_domain_name, omega, ball);
274  return ips;
275 }
276 template idiststream& geo_get_bamg<Float> (idiststream&, geo_basic<Float,sequential>&);
277 
278 }// namespace rheolef
279