Rheolef  7.1
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/index_set_header.icc"
138 #include "rheolef/index_set_body.icc"
139 #include <iostream>
140 #include <fstream>
141 #include <iomanip>
142 #include <limits>
143 #include <string>
144 #include <cstring>
145 #include <vector>
146 #include <list>
147 #include <map>
148 #include <array>
149 
150 using namespace std;
151 using namespace rheolef;
152 void usage() {
153  cerr << "bamg2geo: usage:" << endl
154  << "bamg2geo "
155  << "[-cartesian|-rz|-zr] "
156  << "[input.bamg [input.dmn]]"
157  << endl;
158  exit (1);
159 }
160 void bamg2geo(istream& bamg, istream& dmn, string syscoord) {
161  struct node_type: std::array<double,2> {
162  };
163  struct element_type: std::array<double,4> {
164  size_t nv;
165  size_t size() const { return nv; }
166  };
167  // ------------------------------------------------------------------------
168  // 1) load data
169  // ------------------------------------------------------------------------
170  typedef vector<node_type> n_array_t;
171  typedef vector<element_type> e_array_t;
172  typedef vector<size_t> i_array_t;
173 
174  n_array_t node;
175  i_array_t node_bamg_dom_id;
176  e_array_t element;
177  i_array_t element_bamg_dom_id;
178  e_array_t edge_boundary;
179  i_array_t edge_boundary_bamg_dom_id;
180  list<element_type> edge_list;
181  e_array_t edge;
182  size_t ntri = 0, nqua = 0;
183  string label;
184  while (bamg.good() && label != "End") {
185  bamg >> label;
186  if (label == "Vertices") {
187  // ------------------------------------------------------------------------
188  // 1.1) load the coordinates
189  // Vertices <np>
190  // {xi yi dom_idx} i=0..np-1
191  // ------------------------------------------------------------------------
192  size_t nnod;
193  bamg >> nnod;
194  if (nnod == 0) {
195  cerr << "bamg2geo: error: empty bamg mesh file" << endl;
196  exit(1);
197  }
198  node.resize (nnod);
199  node_bamg_dom_id.resize (nnod);
200  for (size_t inod = 0; inod < nnod; inod++) {
201  bamg >> node[inod][0]
202  >> node[inod][1]
203  >> node_bamg_dom_id[inod];
204  }
205  } else if (label == "Triangles") {
206  // ------------------------------------------------------------------------
207  // 2.1) load the connectivity
208  // Triangle <nt>
209  // {s0 s1 s2 dom_idx} j=0..nt-1
210  // ------------------------------------------------------------------------
211  if (ntri != 0) {
212  cerr << "bamg2geo: error: triangles should precede quadrangles" << endl;
213  exit(1);
214  }
215  bamg >> ntri;
216  element.resize (ntri);
217  element_bamg_dom_id.resize (ntri);
218  for (size_t it = 0; it < ntri; it++) {
219  bamg >> element[it][0]
220  >> element[it][1]
221  >> element[it][2]
222  >> element_bamg_dom_id[it];
223  element[it][0]--;
224  element[it][1]--;
225  element[it][2]--;
226  element[it].nv = 3;
227  }
228  } else if (label == "Quadrilaterals") {
229  // ------------------------------------------------------------------------
230  // Quadrilaterals <nq>
231  // {s0 s1 s2 s3 dom_idx} j=0..nq-1
232  // ------------------------------------------------------------------------
233  bamg >> nqua;
234  cerr << "nqua = " << nqua << endl;
235  element.resize (ntri+nqua);
236  element_bamg_dom_id.resize (ntri+nqua);
237  for (size_t iq = 0; iq < nqua; iq++) {
238  bamg >> element[ntri+iq][0]
239  >> element[ntri+iq][1]
240  >> element[ntri+iq][2]
241  >> element[ntri+iq][3]
242  >> element_bamg_dom_id[ntri+iq];
243  element[ntri+iq][0]--;
244  element[ntri+iq][1]--;
245  element[ntri+iq][2]--;
246  element[ntri+iq][3]--;
247  element[ntri+iq].nv = 4;
248  }
249  } else if (label == "Edges") {
250  // ------------------------------------------------------------------------
251  // 2.3) load the boundary domains
252  // Edges <nedg>
253  // {s0 s1 dom_idx} j=0..nedg-1
254  // ------------------------------------------------------------------------
255  size_t nedg;
256  bamg >> nedg;
257  edge_boundary.resize (nedg);
258  edge_boundary_bamg_dom_id.resize (nedg);
259  for (size_t iedg = 0; iedg < nedg; iedg++) {
260  bamg >> edge_boundary[iedg][0]
261  >> edge_boundary[iedg][1]
262  >> edge_boundary_bamg_dom_id[iedg];
263  edge_boundary[iedg][0]--;
264  edge_boundary[iedg][1]--;
265  edge_boundary[iedg].nv = 2;
266  }
267  }
268  }
269  // ---------------------------------------------------------------
270  // 2) check rheolef extension to optional domain names
271  // ---------------------------------------------------------------
272  vector<string> node_domain_name;
273  vector<string> edge_domain_name;
274  vector<string> region_domain_name;
275  char c;
276  dmn >> ws >> c; // skip white and grab a char
277  // have "EdgeDomainNames" or "VertexDomainNames" ?
278  // bamg mesh may be followed by field data and such, so be carrefull...
279  while (c == 'E' || c == 'V' || c == 'R') {
280  dmn.unget(); // put char back
281  if (c == 'R') {
282  if (!scatch(dmn,"RegionDomainNames",true)) break;
283  size_t n_dom_region;
284  dmn >> n_dom_region;
285  region_domain_name.resize (n_dom_region);
286  for (size_t k = 0; k < n_dom_region; k++) {
287  dmn >> region_domain_name[k];
288  }
289  } else if (c == 'E') {
290  if (!scatch(dmn,"EdgeDomainNames",true)) break;
291  // ---------------------------------------------------------------
292  // get edge domain names
293  // ---------------------------------------------------------------
294  size_t n_dom_edge;
295  dmn >> n_dom_edge;
296  edge_domain_name.resize (n_dom_edge);
297  for (size_t k = 0; k < n_dom_edge; k++) {
298  dmn >> edge_domain_name[k];
299  }
300  } else if (c == 'V') {
301  if (!scatch(dmn,"VertexDomainNames",true)) break;
302  // ---------------------------------------------------------------
303  // get vertice domain names
304  // ---------------------------------------------------------------
305  size_t n_dom_vertice;
306  dmn >> n_dom_vertice;
307  node_domain_name.resize (n_dom_vertice);
308  for (size_t k = 0; k < n_dom_vertice; k++) {
309  dmn >> node_domain_name[k];
310  }
311  }
312  dmn >> ws >> c; // skip white and grab a char
313  }
314  // ------------------------------------------------------------------------
315  // 3) compute all edges
316  // ------------------------------------------------------------------------
317  vector<index_set> ball_edge (node.size());
318  size_t nedg = 0;
319  for (size_t ie = 0; ie < element.size(); ++ie) {
320  for (size_t iv0 = 0; iv0 < element[ie].nv; ++iv0) {
321  size_t iv1 = (iv0+1) % element[ie].nv;
322  size_t inod0 = element[ie][iv0];
323  size_t inod1 = element[ie][iv1];
324  index_set iedge_set = ball_edge[inod0];
325  iedge_set.inplace_intersection (ball_edge[inod1]);
326  if (iedge_set.size() > 1) {
327  cerr << "bamg2geo: error: connectivity problem (iedge.size="<<iedge_set.size()<<")" << endl;
328  exit(1);
329  }
330  if (iedge_set.size() == 1) continue; // edge already exists
331  ball_edge[inod0].insert (nedg);
332  ball_edge[inod1].insert (nedg);
333  element_type new_edge;
334  new_edge[0] = inod0;
335  new_edge[1] = inod1;
336  new_edge.nv = 2;
337  edge_list.push_back (new_edge);
338  nedg++;
339  }
340  }
341  edge.resize (nedg);
342  std::copy (edge_list.begin(), edge_list.end(), edge.begin());
343  // ------------------------------------------------------------------------
344  // 5) build 1d domains
345  // ------------------------------------------------------------------------
346  // 5.1) reduce the edge bamg domain id to [0:edge_ndom[ by counting domain id
347  typedef pair<size_t,size_t> pair_t;
348  typedef map<size_t, pair_t> map_t;
349  map_t edge_bamg_id2idom; // map[bamg_id] = {idom,size}
350  // map[bamg_id] will gives idom and the number of its elements
351  size_t edge_idom = 0;
352  for (size_t ieb = 0, neb = edge_boundary_bamg_dom_id.size(); ieb < neb; ieb++) {
353  size_t bamg_id = edge_boundary_bamg_dom_id [ieb];
354  typename map_t::iterator iter = edge_bamg_id2idom.find (bamg_id);
355  if (iter != edge_bamg_id2idom.end()) {
356  // here is a previous bamg_dom_id: increment #element counter
357  ((*iter).second.second)++;
358  continue;
359  }
360  // here is a new bamg_dom_id: associated to idom and elt counter=1
361  edge_bamg_id2idom.insert (pair<size_t,pair_t>(bamg_id, pair_t(edge_idom,1)));
362  edge_idom++;
363  }
364  size_t edge_ndom = edge_bamg_id2idom.size();
365  // 5.2) check that edge_ndom matches the domain name disarray size
366  if (edge_ndom != edge_domain_name.size()) {
367  cerr << "bamg2geo: error: "
368  << edge_domain_name.size() << " domain name(s) founded while "
369  << edge_ndom << " bamg 1d domain(s) are defined" << endl
370  << "HINT: check domain name file (.dmn)" << endl;
371  exit (1);
372  }
373  // 5.3) convert edge_boundary into an index in the edge[] table with orient
374  struct edge_orient_type { size_t index; int orient; };
375  vector<list<edge_orient_type> > edge_domain (edge_ndom);
376  for (size_t ieb = 0, neb = edge_boundary_bamg_dom_id.size(); ieb < neb; ieb++) {
377  size_t bamg_dom_id = edge_boundary_bamg_dom_id [ieb];
378  size_t edge_idom = edge_bamg_id2idom [bamg_dom_id].first;
379  size_t inod0 = edge_boundary[ieb][0];
380  size_t inod1 = edge_boundary[ieb][1];
381  index_set iedge_set = ball_edge[inod0];
382  iedge_set.inplace_intersection (ball_edge[inod1]);
383  if (iedge_set.size() != 1) {
384  cerr << "bamg2geo: error: connectivity problem (iedge.size="<<iedge_set.size()<<")" << endl;
385  exit(1);
386  }
387  size_t ie = *(iedge_set.begin());
388  edge_orient_type eo;
389  eo.index = ie;
390  eo.orient = (inod0 == edge[ie][0]) ? 1 : -1;
391  edge_domain[edge_idom].push_back (eo);
392  }
393  // ------------------------------------------------------------------------
394  // 6) build 0d domains
395  // ------------------------------------------------------------------------
396  vector<list<size_t> > node_domain;
397  if (node_domain_name.size() != 0) {
398  // 6.1) list all node domain id
399  std::set<size_t> node_id;
400  for (size_t iv = 0, nv = node_bamg_dom_id.size(); iv < nv; ++iv) {
401  size_t dom_id = node_bamg_dom_id [iv];
402  if (dom_id == 0) continue;
403  node_id.insert (dom_id);
404  }
405  cerr << "node_id1.size = " << node_id.size() << endl;
406 
407  // 6.2) omit nodes that are marked with an edge id
408  for (size_t iedg_bdr = 0, nedg_bdr = edge_boundary_bamg_dom_id.size(); iedg_bdr < nedg_bdr; ++iedg_bdr) {
409  size_t dom_id = edge_boundary_bamg_dom_id [iedg_bdr];
410  node_id.erase (dom_id);
411  }
412  cerr << "node_id2.size = " << node_id.size() << endl;
413  if (node_id.size() != node_domain_name.size()) {
414  cerr << "bamg2geo: error: unexpected VertexDomainNames with "<<node_domain_name.size()
415  <<" domain names while the mesh provides " << node_id.size()
416  <<" node labels" << endl;
417  exit(1);
418  }
419  // 6.3) build 0d domain table
420  node_domain.resize (node_id.size());
421  size_t node_idom = 0;
422  for (set<size_t>::const_iterator iter = node_id.begin(); iter != node_id.end(); ++iter, ++node_idom) {
423  size_t bamg_id = *iter;
424  string name = node_domain_name[node_idom];
425  for (size_t i = 0, n = node_bamg_dom_id.size(); i < n; ++i) {
426  if (node_bamg_dom_id [i] != bamg_id) continue;
427  node_domain[node_idom].push_back (i);
428  }
429  }
430  }
431  // ------------------------------------------------------------------------
432  // 7) write .geo
433  // ------------------------------------------------------------------------
434  cout << "#!geo" << endl
435  << endl
436  << "mesh" << endl
437  << "4" << endl
438  << "header" << endl
439  << " dimension 2" << endl;
440  if (syscoord != "cartesian") {
441  cout << " coordinate_system " << syscoord << endl;
442  }
443  cout << " nodes " << node.size() << endl;
444  if (ntri != 0) {
445  cout << " triangles " << ntri << endl;
446  }
447  if (nqua != 0) {
448  cout << " quadrangles " << nqua << endl;
449  }
450  if (edge.size() != 0) {
451  cout << " edges " << edge.size() << endl;
452  }
453  cout << "end header" << endl;
454 
455  cout << setprecision(numeric_limits<double>::digits10);
456  for (size_t i = 0; i < node.size(); ++i) {
457  cout << node[i][0] << " "
458  << node[i][1] << endl;
459  }
460  cout << endl;
461  for (size_t i = 0; i < element.size(); ++i) {
462  cout << ((element[i].nv == 3) ? 't' :'q')
463  << "\t";
464  for (size_t iv = 0; iv < element[i].nv; ++iv) {
465  cout << element[i][iv];
466  if (iv+1 < element[i].nv) { cout << " "; }
467  }
468  cout << endl;
469  }
470  cout << endl;
471  for (auto e: edge) {
472  cout << "e\t"
473  << e[0] << " "
474  << e[1] << endl;
475  }
476  cout << endl;
477  for (size_t idom = 0; idom < edge_domain.size(); ++idom) {
478  cout << "domain" << endl
479  << edge_domain_name[idom] << endl
480  << "2 1 " << edge_domain[idom].size() << endl;
481  for (auto e: edge_domain[idom]) {
482  cout << e.index*e.orient << endl;
483  }
484  cout << endl;
485  }
486  for (size_t idom = 0; idom < node_domain.size(); ++idom) {
487  cout << "domain" << endl
488  << node_domain_name[idom] << endl
489  << "2 0 " << node_domain[idom].size() << endl;
490  for (auto i: node_domain[idom]) {
491  cout << i << endl;
492  }
493  cout << endl;
494  }
495  // TODO: 2d region domains
496 }
497 int main(int argc, char **argv) {
498  string syscoord = "cartesian";
499  string bamg_name, dmn_name;
500  for (int i = 1; i < argc; i++) {
501  if (strcmp (argv[i], "-cartesian") == 0) syscoord = "cartesian";
502  else if (strcmp (argv[i], "-rz") == 0) syscoord = "rz";
503  else if (strcmp (argv[i], "-zr") == 0) syscoord = "zr";
504  else if (argv[i][0] != '-') {
505  // input field on file
506  if (bamg_name == "") bamg_name = argv[i];
507  else if (dmn_name == "") dmn_name = argv[i];
508  else {
509  cerr << "bamg2geo: too many file names `" << argv[i] << "'" << endl;
510  usage();
511  }
512  } else {
513  cerr << "bamg2geo: unknown option `" << argv[i] << "'" << endl;
514  usage();
515  }
516  }
517  if (bamg_name == "") {
518  bamg2geo (cin,cin,syscoord);
519  return 0;
520  }
521  ifstream bamg (bamg_name.c_str());
522  if (!bamg.good()) {
523  cerr << "bamg2geo: unable to read file \""<<bamg_name<<"\"" << endl; exit (1);
524  }
525  if (dmn_name == "") {
526  bamg2geo (bamg,bamg,syscoord);
527  return 0;
528  }
529  ifstream dmn (dmn_name.c_str());
530  if (!dmn.good()) {
531  cerr << "bamg2geo: unable to read file \""<<dmn_name<<"\"" << endl; exit (1);
532  }
533  bamg2geo (bamg,dmn,syscoord);
534 }
mkgeo_ball.n
int n
Definition: mkgeo_ball.sh:150
edge
const size_t edge[n_edge][2]
Definition: hexahedron.icc:118
rheolef::index_set::inplace_intersection
void inplace_intersection(const index_set &b)
Definition: index_set_body.icc:126
usage
void usage()
Definition: bamg2geo.cc:152
rheolef::index_set
Definition: index_set_header.icc:58
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::scatch
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
scatch.icc
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
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
main
int main(int argc, char **argv)
Definition: bamg2geo.cc:497
mkgeo_couette.dmn
string dmn
Definition: mkgeo_couette.sh:122
bamg2geo
void bamg2geo(istream &bamg, istream &dmn, string syscoord)
Definition: bamg2geo.cc:160
rheolef::std
Definition: vec_expr_v2.h:402
mkgeo_contraction.name
string name
Definition: mkgeo_contraction.sh:133
mkgeo_ball.c
int c
Definition: mkgeo_ball.sh:153
edge
see the edge page for the full documentation