Rheolef  7.2
an efficient C++ finite element environment
 
Loading...
Searching...
No Matches
bamg2geo.cc
Go to the documentation of this file.
1
3//
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
25namespace 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
151using namespace std;
152using namespace rheolef;
153void usage() {
154 cerr << "bamg2geo: usage:" << endl
155 << "bamg2geo "
156 << "[-cartesian|-rz|-zr] "
157 << "[input.bamg [input.dmn]]"
158 << endl;
159 exit (1);
160}
161void 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}
498int 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
see the edge page for the full documentation
void inplace_intersection(const index_set &b)
int main()
Definition field2bb.cc:58
verbose clean transpose logscale grid shrink ball stereo iso volume skipvtk deformation fastfieldload lattice reader_on_stdin color format format bamg
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
STL namespace.