bamg2geo.cc
Go to the documentation of this file.
//
// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
//
// Rheolef is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// Rheolef is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with Rheolef; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
//
// =========================================================================
// the bamg2geo unix command
// author: Pierre.Saramito@imag.fr
// date: 4 january 2011
namespace rheolef {
} // namespace rheolef
#include "scatch.icc"
#include "rheolef/index_set_header.icc"
#include "rheolef/index_set_body.icc"
#include <iostream>
#include <fstream>
#include <iomanip>
#include <limits>
#include <string>
#include <cstring>
#include <vector>
#include <list>
#include <map>
#include <array>
using namespace std;
using namespace rheolef;
void usage() {
cerr << "bamg2geo: usage:" << endl
<< "bamg2geo "
<< "[-cartesian|-rz|-zr] "
<< "[input.bamg [input.dmn]]"
<< endl;
exit (1);
}
void bamg2geo(istream& bamg, istream& dmn, string syscoord) {
struct node_type: std::array<double,2> {
};
struct element_type: std::array<double,4> {
size_t nv;
size_t size() const { return nv; }
};
// ------------------------------------------------------------------------
// 1) load data
// ------------------------------------------------------------------------
typedef vector<node_type> n_array_t;
typedef vector<element_type> e_array_t;
typedef vector<size_t> i_array_t;
n_array_t node;
i_array_t node_bamg_dom_id;
e_array_t element;
i_array_t element_bamg_dom_id;
e_array_t edge_boundary;
i_array_t edge_boundary_bamg_dom_id;
list<element_type> edge_list;
e_array_t edge;
size_t ntri = 0, nqua = 0;
string label;
while (bamg.good() && label != "End") {
bamg >> label;
if (label == "Vertices") {
// ------------------------------------------------------------------------
// 1.1) load the coordinates
// Vertices <np>
// {xi yi dom_idx} i=0..np-1
// ------------------------------------------------------------------------
size_t nnod;
bamg >> nnod;
if (nnod == 0) {
cerr << "bamg2geo: error: empty bamg mesh file" << endl;
exit(1);
}
node.resize (nnod);
node_bamg_dom_id.resize (nnod);
for (size_t inod = 0; inod < nnod; inod++) {
bamg >> node[inod][0]
>> node[inod][1]
>> node_bamg_dom_id[inod];
}
} else if (label == "Triangles") {
// ------------------------------------------------------------------------
// 2.1) load the connectivity
// Triangle <nt>
// {s0 s1 s2 dom_idx} j=0..nt-1
// ------------------------------------------------------------------------
if (ntri != 0) {
cerr << "bamg2geo: error: triangles should precede quadrangles" << endl;
exit(1);
}
bamg >> ntri;
element.resize (ntri);
element_bamg_dom_id.resize (ntri);
for (size_t it = 0; it < ntri; it++) {
bamg >> element[it][0]
>> element[it][1]
>> element[it][2]
>> element_bamg_dom_id[it];
element[it][0]--;
element[it][1]--;
element[it][2]--;
element[it].nv = 3;
}
} else if (label == "Quadrilaterals") {
// ------------------------------------------------------------------------
// Quadrilaterals <nq>
// {s0 s1 s2 s3 dom_idx} j=0..nq-1
// ------------------------------------------------------------------------
bamg >> nqua;
cerr << "nqua = " << nqua << endl;
element.resize (ntri+nqua);
element_bamg_dom_id.resize (ntri+nqua);
for (size_t iq = 0; iq < nqua; iq++) {
bamg >> element[ntri+iq][0]
>> element[ntri+iq][1]
>> element[ntri+iq][2]
>> element[ntri+iq][3]
>> element_bamg_dom_id[ntri+iq];
element[ntri+iq][0]--;
element[ntri+iq][1]--;
element[ntri+iq][2]--;
element[ntri+iq][3]--;
element[ntri+iq].nv = 4;
}
} else if (label == "Edges") {
// ------------------------------------------------------------------------
// 2.3) load the boundary domains
// Edges <nedg>
// {s0 s1 dom_idx} j=0..nedg-1
// ------------------------------------------------------------------------
size_t nedg;
bamg >> nedg;
edge_boundary.resize (nedg);
edge_boundary_bamg_dom_id.resize (nedg);
for (size_t iedg = 0; iedg < nedg; iedg++) {
bamg >> edge_boundary[iedg][0]
>> edge_boundary[iedg][1]
>> edge_boundary_bamg_dom_id[iedg];
edge_boundary[iedg][0]--;
edge_boundary[iedg][1]--;
edge_boundary[iedg].nv = 2;
}
}
}
// ---------------------------------------------------------------
// 2) check rheolef extension to optional domain names
// ---------------------------------------------------------------
vector<string> node_domain_name;
vector<string> edge_domain_name;
vector<string> region_domain_name;
char c;
dmn >> ws >> c; // skip white and grab a char
// have "EdgeDomainNames" or "VertexDomainNames" ?
// bamg mesh may be followed by field data and such, so be carrefull...
while (c == 'E' || c == 'V' || c == 'R') {
dmn.unget(); // put char back
if (c == 'R') {
if (!scatch(dmn,"RegionDomainNames",true)) break;
size_t n_dom_region;
dmn >> n_dom_region;
region_domain_name.resize (n_dom_region);
for (size_t k = 0; k < n_dom_region; k++) {
dmn >> region_domain_name[k];
}
} else if (c == 'E') {
if (!scatch(dmn,"EdgeDomainNames",true)) break;
// ---------------------------------------------------------------
// get edge domain names
// ---------------------------------------------------------------
size_t n_dom_edge;
dmn >> n_dom_edge;
edge_domain_name.resize (n_dom_edge);
for (size_t k = 0; k < n_dom_edge; k++) {
dmn >> edge_domain_name[k];
}
} else if (c == 'V') {
if (!scatch(dmn,"VertexDomainNames",true)) break;
// ---------------------------------------------------------------
// get vertice domain names
// ---------------------------------------------------------------
size_t n_dom_vertice;
dmn >> n_dom_vertice;
node_domain_name.resize (n_dom_vertice);
for (size_t k = 0; k < n_dom_vertice; k++) {
dmn >> node_domain_name[k];
}
}
dmn >> ws >> c; // skip white and grab a char
}
// ------------------------------------------------------------------------
// 3) compute all edges
// ------------------------------------------------------------------------
vector<index_set> ball_edge (node.size());
size_t nedg = 0;
for (size_t ie = 0; ie < element.size(); ++ie) {
for (size_t iv0 = 0; iv0 < element[ie].nv; ++iv0) {
size_t iv1 = (iv0+1) % element[ie].nv;
size_t inod0 = element[ie][iv0];
size_t inod1 = element[ie][iv1];
index_set iedge_set = ball_edge[inod0];
iedge_set.inplace_intersection (ball_edge[inod1]);
if (iedge_set.size() > 1) {
cerr << "bamg2geo: error: connectivity problem (iedge.size="<<iedge_set.size()<<")" << endl;
exit(1);
}
if (iedge_set.size() == 1) continue; // edge already exists
ball_edge[inod0].insert (nedg);
ball_edge[inod1].insert (nedg);
element_type new_edge;
new_edge[0] = inod0;
new_edge[1] = inod1;
new_edge.nv = 2;
edge_list.push_back (new_edge);
nedg++;
}
}
edge.resize (nedg);
std::copy (edge_list.begin(), edge_list.end(), edge.begin());
// ------------------------------------------------------------------------
// 5) build 1d domains
// ------------------------------------------------------------------------
// 5.1) reduce the edge bamg domain id to [0:edge_ndom[ by counting domain id
typedef pair<size_t,size_t> pair_t;
typedef map<size_t, pair_t> map_t;
map_t edge_bamg_id2idom; // map[bamg_id] = {idom,size}
// map[bamg_id] will gives idom and the number of its elements
size_t edge_idom = 0;
for (size_t ieb = 0, neb = edge_boundary_bamg_dom_id.size(); ieb < neb; ieb++) {
size_t bamg_id = edge_boundary_bamg_dom_id [ieb];
typename map_t::iterator iter = edge_bamg_id2idom.find (bamg_id);
if (iter != edge_bamg_id2idom.end()) {
// here is a previous bamg_dom_id: increment #element counter
((*iter).second.second)++;
continue;
}
// here is a new bamg_dom_id: associated to idom and elt counter=1
edge_bamg_id2idom.insert (pair<size_t,pair_t>(bamg_id, pair_t(edge_idom,1)));
edge_idom++;
}
size_t edge_ndom = edge_bamg_id2idom.size();
// 5.2) check that edge_ndom matches the domain name disarray size
if (edge_ndom != edge_domain_name.size()) {
cerr << "bamg2geo: error: "
<< edge_domain_name.size() << " domain name(s) founded while "
<< edge_ndom << " bamg 1d domain(s) are defined" << endl
<< "HINT: check domain name file (.dmn)" << endl;
exit (1);
}
// 5.3) convert edge_boundary into an index in the edge[] table with orient
struct edge_orient_type { size_t index; int orient; };
vector<list<edge_orient_type> > edge_domain (edge_ndom);
for (size_t ieb = 0, neb = edge_boundary_bamg_dom_id.size(); ieb < neb; ieb++) {
size_t bamg_dom_id = edge_boundary_bamg_dom_id [ieb];
size_t edge_idom = edge_bamg_id2idom [bamg_dom_id].first;
size_t inod0 = edge_boundary[ieb][0];
size_t inod1 = edge_boundary[ieb][1];
index_set iedge_set = ball_edge[inod0];
iedge_set.inplace_intersection (ball_edge[inod1]);
if (iedge_set.size() != 1) {
cerr << "bamg2geo: error: connectivity problem (iedge.size="<<iedge_set.size()<<")" << endl;
exit(1);
}
size_t ie = *(iedge_set.begin());
edge_orient_type eo;
eo.index = ie;
eo.orient = (inod0 == edge[ie][0]) ? 1 : -1;
edge_domain[edge_idom].push_back (eo);
}
// ------------------------------------------------------------------------
// 6) build 0d domains
// ------------------------------------------------------------------------
vector<list<size_t> > node_domain;
if (node_domain_name.size() != 0) {
// 6.1) list all node domain id
std::set<size_t> node_id;
for (size_t iv = 0, nv = node_bamg_dom_id.size(); iv < nv; ++iv) {
size_t dom_id = node_bamg_dom_id [iv];
if (dom_id == 0) continue;
node_id.insert (dom_id);
}
cerr << "node_id1.size = " << node_id.size() << endl;
// 6.2) omit nodes that are marked with an edge id
for (size_t iedg_bdr = 0, nedg_bdr = edge_boundary_bamg_dom_id.size(); iedg_bdr < nedg_bdr; ++iedg_bdr) {
size_t dom_id = edge_boundary_bamg_dom_id [iedg_bdr];
node_id.erase (dom_id);
}
cerr << "node_id2.size = " << node_id.size() << endl;
if (node_id.size() != node_domain_name.size()) {
cerr << "bamg2geo: error: unexpected VertexDomainNames with "<<node_domain_name.size()
<<" domain names while the mesh provides " << node_id.size()
<<" node labels" << endl;
exit(1);
}
// 6.3) build 0d domain table
node_domain.resize (node_id.size());
size_t node_idom = 0;
for (set<size_t>::const_iterator iter = node_id.begin(); iter != node_id.end(); ++iter, ++node_idom) {
size_t bamg_id = *iter;
string name = node_domain_name[node_idom];
for (size_t i = 0, n = node_bamg_dom_id.size(); i < n; ++i) {
if (node_bamg_dom_id [i] != bamg_id) continue;
node_domain[node_idom].push_back (i);
}
}
}
// ------------------------------------------------------------------------
// 7) write .geo
// ------------------------------------------------------------------------
cout << "#!geo" << endl
<< endl
<< "mesh" << endl
<< "4" << endl
<< "header" << endl
<< " dimension 2" << endl;
if (syscoord != "cartesian") {
cout << " coordinate_system " << syscoord << endl;
}
cout << " nodes " << node.size() << endl;
if (ntri != 0) {
cout << " triangles " << ntri << endl;
}
if (nqua != 0) {
cout << " quadrangles " << nqua << endl;
}
if (edge.size() != 0) {
cout << " edges " << edge.size() << endl;
}
cout << "end header" << endl;
cout << setprecision(numeric_limits<double>::digits10);
for (size_t i = 0; i < node.size(); ++i) {
cout << node[i][0] << " "
<< node[i][1] << endl;
}
cout << endl;
for (size_t i = 0; i < element.size(); ++i) {
cout << ((element[i].nv == 3) ? 't' :'q')
<< "\t";
for (size_t iv = 0; iv < element[i].nv; ++iv) {
cout << element[i][iv];
if (iv+1 < element[i].nv) { cout << " "; }
}
cout << endl;
}
cout << endl;
for (auto e: edge) {
cout << "e\t"
<< e[0] << " "
<< e[1] << endl;
}
cout << endl;
for (size_t idom = 0; idom < edge_domain.size(); ++idom) {
cout << "domain" << endl
<< edge_domain_name[idom] << endl
<< "2 1 " << edge_domain[idom].size() << endl;
for (auto e: edge_domain[idom]) {
cout << e.index*e.orient << endl;
}
cout << endl;
}
for (size_t idom = 0; idom < node_domain.size(); ++idom) {
cout << "domain" << endl
<< node_domain_name[idom] << endl
<< "2 0 " << node_domain[idom].size() << endl;
for (auto i: node_domain[idom]) {
cout << i << endl;
}
cout << endl;
}
// TODO: 2d region domains
}
int main(int argc, char **argv) {
string syscoord = "cartesian";
string bamg_name, dmn_name;
for (int i = 1; i < argc; i++) {
if (strcmp (argv[i], "-cartesian") == 0) syscoord = "cartesian";
else if (strcmp (argv[i], "-rz") == 0) syscoord = "rz";
else if (strcmp (argv[i], "-zr") == 0) syscoord = "zr";
else if (argv[i][0] != '-') {
// input field on file
if (bamg_name == "") bamg_name = argv[i];
else if (dmn_name == "") dmn_name = argv[i];
else {
cerr << "bamg2geo: too many file names `" << argv[i] << "'" << endl;
usage();
}
} else {
cerr << "bamg2geo: unknown option `" << argv[i] << "'" << endl;
usage();
}
}
if (bamg_name == "") {
bamg2geo (cin,cin,syscoord);
return 0;
}
ifstream bamg (bamg_name.c_str());
if (!bamg.good()) {
cerr << "bamg2geo: unable to read file \""<<bamg_name<<"\"" << endl; exit (1);
}
if (dmn_name == "") {
bamg2geo (bamg,bamg,syscoord);
return 0;
}
ifstream dmn (dmn_name.c_str());
if (!dmn.good()) {
cerr << "bamg2geo: unable to read file \""<<dmn_name<<"\"" << endl; exit (1);
}
bamg2geo (bamg,dmn,syscoord);
}
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
mkgeo_ball.c
c
Definition: mkgeo_ball.sh:153
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
mkgeo_couette.dmn
dmn
Definition: mkgeo_couette.sh:122
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_ball.n
n
Definition: mkgeo_ball.sh:150
bamg2geo
void bamg2geo(istream &bamg, istream &dmn, string syscoord)
Definition: bamg2geo.cc:160
rheolef::std
Definition: vec_expr_v2.h:391
mkgeo_contraction.name
name
Definition: mkgeo_contraction.sh:133
edge
see the edge page for the full documentation