161void bamg2geo(istream& bamg, istream& dmn,
string syscoord) {
162 struct node_type: std::array<double,2> {
164 struct element_type: std::array<double,4> {
166 size_t size()
const {
return nv; }
171 typedef vector<node_type> n_array_t;
172 typedef vector<element_type> e_array_t;
173 typedef vector<size_t> i_array_t;
176 i_array_t node_bamg_dom_id;
178 i_array_t element_bamg_dom_id;
179 e_array_t edge_boundary;
180 i_array_t edge_boundary_bamg_dom_id;
181 list<element_type> edge_list;
183 size_t ntri = 0, nqua = 0;
185 while (
bamg.good() && label !=
"End") {
187 if (label ==
"Vertices") {
196 cerr <<
"bamg2geo: error: empty bamg mesh file" << endl;
200 node_bamg_dom_id.resize (nnod);
201 for (
size_t inod = 0; inod < nnod; inod++) {
202 bamg >> node[inod][0]
204 >> node_bamg_dom_id[inod];
206 }
else if (label ==
"Triangles") {
213 cerr <<
"bamg2geo: error: triangles should precede quadrangles" << endl;
217 element.resize (ntri);
218 element_bamg_dom_id.resize (ntri);
219 for (
size_t it = 0; it < ntri; it++) {
220 bamg >> element[it][0]
223 >> element_bamg_dom_id[it];
229 }
else if (label ==
"Quadrilaterals") {
235 cerr <<
"nqua = " << nqua << endl;
236 element.resize (ntri+nqua);
237 element_bamg_dom_id.resize (ntri+nqua);
238 for (
size_t iq = 0; iq < nqua; iq++) {
239 bamg >> element[ntri+iq][0]
240 >> element[ntri+iq][1]
241 >> element[ntri+iq][2]
242 >> element[ntri+iq][3]
243 >> element_bamg_dom_id[ntri+iq];
244 element[ntri+iq][0]--;
245 element[ntri+iq][1]--;
246 element[ntri+iq][2]--;
247 element[ntri+iq][3]--;
248 element[ntri+iq].nv = 4;
250 }
else if (label ==
"Edges") {
258 edge_boundary.resize (nedg);
259 edge_boundary_bamg_dom_id.resize (nedg);
260 for (
size_t iedg = 0; iedg < nedg; iedg++) {
261 bamg >> edge_boundary[iedg][0]
262 >> edge_boundary[iedg][1]
263 >> edge_boundary_bamg_dom_id[iedg];
264 edge_boundary[iedg][0]--;
265 edge_boundary[iedg][1]--;
266 edge_boundary[iedg].nv = 2;
273 vector<string> node_domain_name;
274 vector<string> edge_domain_name;
275 vector<string> region_domain_name;
280 while (c ==
'E' || c ==
'V' || c ==
'R') {
283 if (!
scatch(dmn,
"RegionDomainNames",
true))
break;
286 region_domain_name.resize (n_dom_region);
287 for (
size_t k = 0; k < n_dom_region; k++) {
288 dmn >> region_domain_name[k];
290 }
else if (c ==
'E') {
291 if (!
scatch(dmn,
"EdgeDomainNames",
true))
break;
297 edge_domain_name.resize (n_dom_edge);
298 for (
size_t k = 0; k < n_dom_edge; k++) {
299 dmn >> edge_domain_name[k];
301 }
else if (c ==
'V') {
302 if (!
scatch(dmn,
"VertexDomainNames",
true))
break;
306 size_t n_dom_vertice;
307 dmn >> n_dom_vertice;
308 node_domain_name.resize (n_dom_vertice);
309 for (
size_t k = 0; k < n_dom_vertice; k++) {
310 dmn >> node_domain_name[k];
318 vector<index_set> ball_edge (node.size());
320 for (
size_t ie = 0; ie < element.size(); ++ie) {
321 for (
size_t iv0 = 0; iv0 < element[ie].nv; ++iv0) {
322 size_t iv1 = (iv0+1) % element[ie].nv;
323 size_t inod0 = element[ie][iv0];
324 size_t inod1 = element[ie][iv1];
327 if (iedge_set.size() > 1) {
328 cerr <<
"bamg2geo: error: connectivity problem (iedge.size="<<iedge_set.size()<<
")" << endl;
331 if (iedge_set.size() == 1)
continue;
332 ball_edge[inod0].insert (nedg);
333 ball_edge[inod1].insert (nedg);
334 element_type new_edge;
338 edge_list.push_back (new_edge);
343 std::copy (edge_list.begin(), edge_list.end(),
edge.begin());
348 typedef pair<size_t,size_t> pair_t;
349 typedef map<size_t, pair_t> map_t;
350 map_t edge_bamg_id2idom;
352 size_t edge_idom = 0;
353 for (
size_t ieb = 0, neb = edge_boundary_bamg_dom_id.size(); ieb < neb; ieb++) {
354 size_t bamg_id = edge_boundary_bamg_dom_id [ieb];
355 typename map_t::iterator iter = edge_bamg_id2idom.find (bamg_id);
356 if (iter != edge_bamg_id2idom.end()) {
358 ((*iter).second.second)++;
362 edge_bamg_id2idom.insert (pair<size_t,pair_t>(bamg_id, pair_t(edge_idom,1)));
365 size_t edge_ndom = edge_bamg_id2idom.size();
367 if (edge_ndom != edge_domain_name.size()) {
368 cerr <<
"bamg2geo: error: "
369 << edge_domain_name.size() <<
" domain name(s) founded while "
370 << edge_ndom <<
" bamg 1d domain(s) are defined" << endl
371 <<
"HINT: check domain name file (.dmn)" << endl;
375 struct edge_orient_type {
long index;
int orient; };
376 vector<list<edge_orient_type> > edge_domain (edge_ndom);
377 for (
size_t ieb = 0, neb = edge_boundary_bamg_dom_id.size(); ieb < neb; ieb++) {
378 size_t bamg_dom_id = edge_boundary_bamg_dom_id [ieb];
379 size_t edge_idom = edge_bamg_id2idom [bamg_dom_id].first;
380 size_t inod0 = edge_boundary[ieb][0];
381 size_t inod1 = edge_boundary[ieb][1];
384 if (iedge_set.size() != 1) {
385 cerr <<
"bamg2geo: error: connectivity problem (iedge.size="<<iedge_set.size()<<
")" << endl;
388 size_t ie = *(iedge_set.begin());
391 eo.orient = (inod0 ==
edge[ie][0]) ? 1 : -1;
392 edge_domain[edge_idom].push_back (eo);
397 vector<list<size_t> > node_domain;
398 if (node_domain_name.size() != 0) {
400 std::set<size_t> node_id;
401 for (
size_t iv = 0, nv = node_bamg_dom_id.size(); iv < nv; ++iv) {
402 size_t dom_id = node_bamg_dom_id [iv];
403 if (dom_id == 0)
continue;
404 node_id.insert (dom_id);
406 cerr <<
"node_id1.size = " << node_id.size() << endl;
409 for (
size_t iedg_bdr = 0, nedg_bdr = edge_boundary_bamg_dom_id.size(); iedg_bdr < nedg_bdr; ++iedg_bdr) {
410 size_t dom_id = edge_boundary_bamg_dom_id [iedg_bdr];
411 node_id.erase (dom_id);
413 cerr <<
"node_id2.size = " << node_id.size() << endl;
414 if (node_id.size() != node_domain_name.size()) {
415 cerr <<
"bamg2geo: error: unexpected VertexDomainNames with "<<node_domain_name.size()
416 <<
" domain names while the mesh provides " << node_id.size()
417 <<
" node labels" << endl;
421 node_domain.resize (node_id.size());
422 size_t node_idom = 0;
423 for (set<size_t>::const_iterator iter = node_id.begin(); iter != node_id.end(); ++iter, ++node_idom) {
424 size_t bamg_id = *iter;
425 string name = node_domain_name[node_idom];
426 for (
size_t i = 0, n = node_bamg_dom_id.size(); i < n; ++i) {
427 if (node_bamg_dom_id [i] != bamg_id)
continue;
428 node_domain[node_idom].push_back (i);
435 cout <<
"#!geo" << endl
440 <<
" dimension 2" << endl;
441 if (syscoord !=
"cartesian") {
442 cout <<
" coordinate_system " << syscoord << endl;
444 cout <<
" nodes " << node.size() << endl;
446 cout <<
" triangles " << ntri << endl;
449 cout <<
" quadrangles " << nqua << endl;
451 if (
edge.size() != 0) {
452 cout <<
" edges " <<
edge.size() << endl;
454 cout <<
"end header" << endl;
456 cout << setprecision(numeric_limits<double>::digits10);
457 for (
size_t i = 0; i < node.size(); ++i) {
458 cout << node[i][0] <<
" "
459 << node[i][1] << endl;
462 for (
size_t i = 0; i < element.size(); ++i) {
463 cout << ((element[i].nv == 3) ?
't' :
'q')
465 for (
size_t iv = 0; iv < element[i].nv; ++iv) {
466 cout << element[i][iv];
467 if (iv+1 < element[i].nv) { cout <<
" "; }
478 for (
size_t idom = 0; idom < edge_domain.size(); ++idom) {
479 cout <<
"domain" << endl
480 << edge_domain_name[idom] << endl
481 <<
"2 1 " << edge_domain[idom].size() << endl;
482 for (
auto e: edge_domain[idom]) {
483 cout << e.index*e.orient << endl;
487 for (
size_t idom = 0; idom < node_domain.size(); ++idom) {
488 cout <<
"domain" << endl
489 << node_domain_name[idom] << endl
490 <<
"2 0 " << node_domain[idom].size() << endl;
491 for (
auto i: node_domain[idom]) {