Rheolef  7.2
an efficient C++ finite element environment
bamg2geo.cc
Go to the documentation of this file.
1 //
4 // Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
5 //
6 // Rheolef is free software; you can redistribute it and/or modify
7 // it under the terms of the GNU General Public License as published by
8 // the Free Software Foundation; either version 2 of the License, or
9 // (at your option) any later version.
10 //
11 // Rheolef is distributed in the hope that it will be useful,
12 // but WITHOUT ANY WARRANTY; without even the implied warranty of
13 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 // GNU General Public License for more details.
15 //
16 // You should have received a copy of the GNU General Public License
17 // along with Rheolef; if not, write to the Free Software
18 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
19 //
20 // =========================================================================
21 // the bamg2geo unix command
22 // author: Pierre.Saramito@imag.fr
23 // date: 4 january 2011
24 
25 namespace rheolef {
134 } // namespace rheolef
135 
136 #include "scatch.icc"
137 #include "rheolef/compiler.h"
138 #include "rheolef/index_set_header.icc"
139 #include "rheolef/index_set_body.icc"
140 #include <iostream>
141 #include <fstream>
142 #include <iomanip>
143 #include <limits>
144 #include <string>
145 #include <cstring>
146 #include <vector>
147 #include <list>
148 #include <map>
149 #include <array>
150 
151 using namespace std;
152 using namespace rheolef;
153 void usage() {
154  cerr << "bamg2geo: usage:" << endl
155  << "bamg2geo "
156  << "[-cartesian|-rz|-zr] "
157  << "[input.bamg [input.dmn]]"
158  << endl;
159  exit (1);
160 }
161 void bamg2geo(istream& bamg, istream& dmn, string syscoord) {
162  struct node_type: std::array<double,2> {
163  };
164  struct element_type: std::array<double,4> {
165  size_t nv;
166  size_t size() const { return nv; }
167  };
168  // ------------------------------------------------------------------------
169  // 1) load data
170  // ------------------------------------------------------------------------
171  typedef vector<node_type> n_array_t;
172  typedef vector<element_type> e_array_t;
173  typedef vector<size_t> i_array_t;
174 
175  n_array_t node;
176  i_array_t node_bamg_dom_id;
177  e_array_t element;
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;
182  e_array_t edge;
183  size_t ntri = 0, nqua = 0;
184  string label;
185  while (bamg.good() && label != "End") {
186  bamg >> label;
187  if (label == "Vertices") {
188  // ------------------------------------------------------------------------
189  // 1.1) load the coordinates
190  // Vertices <np>
191  // {xi yi dom_idx} i=0..np-1
192  // ------------------------------------------------------------------------
193  size_t nnod;
194  bamg >> nnod;
195  if (nnod == 0) {
196  cerr << "bamg2geo: error: empty bamg mesh file" << endl;
197  exit(1);
198  }
199  node.resize (nnod);
200  node_bamg_dom_id.resize (nnod);
201  for (size_t inod = 0; inod < nnod; inod++) {
202  bamg >> node[inod][0]
203  >> node[inod][1]
204  >> node_bamg_dom_id[inod];
205  }
206  } else if (label == "Triangles") {
207  // ------------------------------------------------------------------------
208  // 2.1) load the connectivity
209  // Triangle <nt>
210  // {s0 s1 s2 dom_idx} j=0..nt-1
211  // ------------------------------------------------------------------------
212  if (ntri != 0) {
213  cerr << "bamg2geo: error: triangles should precede quadrangles" << endl;
214  exit(1);
215  }
216  bamg >> ntri;
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]
221  >> element[it][1]
222  >> element[it][2]
223  >> element_bamg_dom_id[it];
224  element[it][0]--;
225  element[it][1]--;
226  element[it][2]--;
227  element[it].nv = 3;
228  }
229  } else if (label == "Quadrilaterals") {
230  // ------------------------------------------------------------------------
231  // Quadrilaterals <nq>
232  // {s0 s1 s2 s3 dom_idx} j=0..nq-1
233  // ------------------------------------------------------------------------
234  bamg >> nqua;
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;
249  }
250  } else if (label == "Edges") {
251  // ------------------------------------------------------------------------
252  // 2.3) load the boundary domains
253  // Edges <nedg>
254  // {s0 s1 dom_idx} j=0..nedg-1
255  // ------------------------------------------------------------------------
256  size_t nedg;
257  bamg >> nedg;
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;
267  }
268  }
269  }
270  // ---------------------------------------------------------------
271  // 2) check rheolef extension to optional domain names
272  // ---------------------------------------------------------------
273  vector<string> node_domain_name;
274  vector<string> edge_domain_name;
275  vector<string> region_domain_name;
276  char c;
277  dmn >> ws >> c; // skip white and grab a char
278  // have "EdgeDomainNames" or "VertexDomainNames" ?
279  // bamg mesh may be followed by field data and such, so be carrefull...
280  while (c == 'E' || c == 'V' || c == 'R') {
281  dmn.unget(); // put char back
282  if (c == 'R') {
283  if (!scatch(dmn,"RegionDomainNames",true)) break;
284  size_t n_dom_region;
285  dmn >> n_dom_region;
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];
289  }
290  } else if (c == 'E') {
291  if (!scatch(dmn,"EdgeDomainNames",true)) break;
292  // ---------------------------------------------------------------
293  // get edge domain names
294  // ---------------------------------------------------------------
295  size_t n_dom_edge;
296  dmn >> n_dom_edge;
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];
300  }
301  } else if (c == 'V') {
302  if (!scatch(dmn,"VertexDomainNames",true)) break;
303  // ---------------------------------------------------------------
304  // get vertice domain names
305  // ---------------------------------------------------------------
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];
311  }
312  }
313  dmn >> ws >> c; // skip white and grab a char
314  }
315  // ------------------------------------------------------------------------
316  // 3) compute all edges
317  // ------------------------------------------------------------------------
318  vector<index_set> ball_edge (node.size());
319  size_t nedg = 0;
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];
325  index_set iedge_set = ball_edge[inod0];
326  iedge_set.inplace_intersection (ball_edge[inod1]);
327  if (iedge_set.size() > 1) {
328  cerr << "bamg2geo: error: connectivity problem (iedge.size="<<iedge_set.size()<<")" << endl;
329  exit(1);
330  }
331  if (iedge_set.size() == 1) continue; // edge already exists
332  ball_edge[inod0].insert (nedg);
333  ball_edge[inod1].insert (nedg);
334  element_type new_edge;
335  new_edge[0] = inod0;
336  new_edge[1] = inod1;
337  new_edge.nv = 2;
338  edge_list.push_back (new_edge);
339  nedg++;
340  }
341  }
342  edge.resize (nedg);
343  std::copy (edge_list.begin(), edge_list.end(), edge.begin());
344  // ------------------------------------------------------------------------
345  // 5) build 1d domains
346  // ------------------------------------------------------------------------
347  // 5.1) reduce the edge bamg domain id to [0:edge_ndom[ by counting domain id
348  typedef pair<size_t,size_t> pair_t;
349  typedef map<size_t, pair_t> map_t;
350  map_t edge_bamg_id2idom; // map[bamg_id] = {idom,size}
351  // map[bamg_id] will gives idom and the number of its elements
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()) {
357  // here is a previous bamg_dom_id: increment #element counter
358  ((*iter).second.second)++;
359  continue;
360  }
361  // here is a new bamg_dom_id: associated to idom and elt counter=1
362  edge_bamg_id2idom.insert (pair<size_t,pair_t>(bamg_id, pair_t(edge_idom,1)));
363  edge_idom++;
364  }
365  size_t edge_ndom = edge_bamg_id2idom.size();
366  // 5.2) check that edge_ndom matches the domain name disarray 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;
372  exit (1);
373  }
374  // 5.3) convert edge_boundary into an index in the edge[] table with orient
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];
382  index_set iedge_set = ball_edge[inod0];
383  iedge_set.inplace_intersection (ball_edge[inod1]);
384  if (iedge_set.size() != 1) {
385  cerr << "bamg2geo: error: connectivity problem (iedge.size="<<iedge_set.size()<<")" << endl;
386  exit(1);
387  }
388  size_t ie = *(iedge_set.begin());
389  edge_orient_type eo;
390  eo.index = ie;
391  eo.orient = (inod0 == edge[ie][0]) ? 1 : -1;
392  edge_domain[edge_idom].push_back (eo);
393  }
394  // ------------------------------------------------------------------------
395  // 6) build 0d domains
396  // ------------------------------------------------------------------------
397  vector<list<size_t> > node_domain;
398  if (node_domain_name.size() != 0) {
399  // 6.1) list all node domain id
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);
405  }
406  cerr << "node_id1.size = " << node_id.size() << endl;
407 
408  // 6.2) omit nodes that are marked with an edge id
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);
412  }
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;
418  exit(1);
419  }
420  // 6.3) build 0d domain table
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);
429  }
430  }
431  }
432  // ------------------------------------------------------------------------
433  // 7) write .geo
434  // ------------------------------------------------------------------------
435  cout << "#!geo" << endl
436  << endl
437  << "mesh" << endl
438  << "4" << endl
439  << "header" << endl
440  << " dimension 2" << endl;
441  if (syscoord != "cartesian") {
442  cout << " coordinate_system " << syscoord << endl;
443  }
444  cout << " nodes " << node.size() << endl;
445  if (ntri != 0) {
446  cout << " triangles " << ntri << endl;
447  }
448  if (nqua != 0) {
449  cout << " quadrangles " << nqua << endl;
450  }
451  if (edge.size() != 0) {
452  cout << " edges " << edge.size() << endl;
453  }
454  cout << "end header" << endl;
455 
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;
460  }
461  cout << endl;
462  for (size_t i = 0; i < element.size(); ++i) {
463  cout << ((element[i].nv == 3) ? 't' :'q')
464  << "\t";
465  for (size_t iv = 0; iv < element[i].nv; ++iv) {
466  cout << element[i][iv];
467  if (iv+1 < element[i].nv) { cout << " "; }
468  }
469  cout << endl;
470  }
471  cout << endl;
472  for (auto e: edge) {
473  cout << "e\t"
474  << e[0] << " "
475  << e[1] << endl;
476  }
477  cout << endl;
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;
484  }
485  cout << endl;
486  }
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]) {
492  cout << i << endl;
493  }
494  cout << endl;
495  }
496  // TODO: 2d region domains
497 }
498 int main(int argc, char **argv) {
499  string syscoord = "cartesian";
500  string bamg_name, dmn_name;
501  for (int i = 1; i < argc; i++) {
502  if (strcmp (argv[i], "-cartesian") == 0) syscoord = "cartesian";
503  else if (strcmp (argv[i], "-rz") == 0) syscoord = "rz";
504  else if (strcmp (argv[i], "-zr") == 0) syscoord = "zr";
505  else if (argv[i][0] != '-') {
506  // input field on file
507  if (bamg_name == "") bamg_name = argv[i];
508  else if (dmn_name == "") dmn_name = argv[i];
509  else {
510  cerr << "bamg2geo: too many file names `" << argv[i] << "'" << endl;
511  usage();
512  }
513  } else {
514  cerr << "bamg2geo: unknown option `" << argv[i] << "'" << endl;
515  usage();
516  }
517  }
518  if (bamg_name == "") {
519  bamg2geo (cin,cin,syscoord);
520  return 0;
521  }
522  ifstream bamg (bamg_name.c_str());
523  if (!bamg.good()) {
524  cerr << "bamg2geo: unable to read file \""<<bamg_name<<"\"" << endl; exit (1);
525  }
526  if (dmn_name == "") {
527  bamg2geo (bamg,bamg,syscoord);
528  return 0;
529  }
530  ifstream dmn (dmn_name.c_str());
531  if (!dmn.good()) {
532  cerr << "bamg2geo: unable to read file \""<<dmn_name<<"\"" << endl; exit (1);
533  }
534  bamg2geo (bamg,dmn,syscoord);
535 }
void bamg2geo(istream &bamg, istream &dmn, string syscoord)
Definition: bamg2geo.cc:161
void usage()
Definition: bamg2geo.cc:153
int main(int argc, char **argv)
Definition: bamg2geo.cc:498
see the edge page for the full documentation
void inplace_intersection(const index_set &b)
const size_t edge[n_edge][2]
Definition: hexahedron.icc:118
verbose clean transpose logscale grid shrink ball stereo iso volume skipvtk deformation fastfieldload lattice reader_on_stdin color format format bamg
size_type nnod(const basis_basic< T > &b, const geo_size &gs, size_type map_dim)
This file is part of Rheolef.
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:44