2
/// This file is part of Rheolef.
4
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
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.
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.
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
20
/// =========================================================================
23
NAME: @code{msh2geo} -- msh2geo - convert gmsh mesh in geo format
27
@fiindex @file{.msh} gmsh mesh
28
@fiindex @file{.geo} mesh
32
msh2geo [-zr|-rz] < @var{input}[.msh] > @var{output}.geo
36
Convert a gmsh @file{.msh} into @file{.geo} one.
37
The output goes to standart output.
38
See the @code{gmsh} documentation for a detailed description
39
of the @file{.mshcad} input file for @code{gmsh}.
42
gmsh -2 toto.mshcad -o toto.msh
43
msh2geo < toto.msh > toto.geo
44
gmsh -2 -order 2 toto.mshcad -o toto2.msh
45
msh2geo < toto2.msh > toto2.geo
47
COORDINATE SYSTEM OPTION:
48
Most of rheolef codes are coordinate-system independant.
49
The coordinate system is specified in the geometry file @file{.geo}.
51
@cindex axisymmetric coordinate system
54
the 2d mesh is axisymmetric: @code{zr} (resp. @code{rz}) stands when the symmetry is related
55
to the first (resp. second) coordinate.
58
Pk triangle, when k>=5, may have internal nodes renumbered: from the
61
The nodes of a curved element are numbered in the following order:
63
the element principal vertices;
64
the internal nodes for each edge;
65
the internal nodes for each face;
66
the volume internal nodes.
68
The numbering for face and volume internal nodes is recursive,
69
i.e., the numbering follows that of the nodes of an embedded face/volume.
70
The higher order nodes are assumed to be equispaced on the element.
73
In rheolef, internal triangle nodes are numbered from left to right and then
74
from bottom to top. The numbering differ for triangle when k >= 5.
75
Thus, @code{msh2geo} fix the corresponding internal nodes numbering during the
78
Pk tetrahedrons and hexaedrons in gmsh and rheolef has not the same edge-node order
80
E.g. for tetrahedrons, edges 13 and 23 should be swaped
81
and reoriented as 32 and 31.
82
Thus, @code{msh2geo} fix the corresponding internal nodes numbering.
85
Fix for P3-tetra: swap edges orientations for 3,4,5
86
and swap faces 1 and 2. Check P4(T) for face orientation.
87
Perform face visualisation with gnuplot face fill.
89
See also hexa edges orient and faces numbers and orient.
91
Check that node are numbered by vertex-node, then edge-node, then face(tri,qua)-node and then volume(T,P,H)-node.
92
Otherwise, renumber all nodes.
94
Support for high order >= 6 element ? not documented in gmsh, but gmsh supports it at run
97
#include "rheolef/point.h"
98
#include "rheolef/index_set.h"
99
#include "rheolef/reference_element.h"
100
#include "rheolef/rheostream.h"
103
using namespace rheolef;
105
#include "msh2geo_defs.icc"
106
void msh2geo_node_renum (vector<size_t>& element, size_type variant, size_type order);
113
const map<size_t,list<size_t> >& domain_map,
114
map<size_t, string>& phys,
115
const vector<char>& element_name,
116
const vector<size_t>& element_order,
117
const vector<vector<size_t> >& element,
118
const vector<size_t>& old2new_inode)
120
if (dim == dom_dim && domain_map.size() == 1) return;
121
for (map<size_t,list<size_t> >::const_iterator
122
first = domain_map.begin(),
123
last = domain_map.end(); first != last; first++) {
124
size_type dom_idx = (*first).first;
125
const list<size_t>& dom = (*first).second;
126
string dom_name = phys[dom_idx];
127
if (dom_name == "") dom_name = "unamed" + itos(dom_idx);
128
out << "domain" << endl
130
<< "1 " << dom_dim << " " << dom.size() << endl;
131
for (size_type variant = reference_element::first_variant_by_dimension(dom_dim);
132
variant < reference_element:: last_variant_by_dimension(dom_dim); variant++) {
133
for (list<size_t>::const_iterator fe = dom.begin(), le = dom.end(); fe != le; fe++) {
135
if (reference_element::variant(element_name[i]) != variant) continue;
136
out << element_name[i] << "\t";
137
if (element_order[i] > 1) {
138
out << "p" << element_order[i] << " ";
140
for (size_type j = 0; j < element[i].size(); j++) {
141
out << old2new_inode[element[i][j]] << " ";
149
void msh2geo (istream& in, ostream& out, std::string sys_coord_name = "cartesian")
151
// ----------------------------------
153
// ----------------------------------
155
check_macro (scatch(in,"$MeshFormat"),
156
"input stream does not contains a gmsh mesh file ($MeshFormat not found).");
157
Float gmsh_fmt_version;
158
size_type file_type, float_data_size;
159
in >> gmsh_fmt_version >> file_type >> float_data_size;
160
if (gmsh_fmt_version > 2.2) {
161
warning_macro("gmsh format version " << gmsh_fmt_version << " founded ; expect version 2.2");
163
check_macro (file_type == 0, "unsupported gmsh non-ascii format");
164
check_macro (scatch(in,"$EndMeshFormat"), "gmsh input error: $EndMeshFormat not found.");
166
// 1.1 optional domain names
168
check_macro (scatch(in,"$"), "gmsh input error: no more label found.");
170
map<size_t, string> phys;
173
size_type n_names = 0;
174
if (label == "PhysicalNames") {
176
for (size_type i = 0; i < nphys; i++) {
178
size_type dom_dim, dom_idx;
179
in >> dom_dim >> dom_idx;
183
if (c != '"') name.push_back(c);
187
} while (c != '"' && c != '\n');
189
size_type start = 0, end = name.length();
190
if (name[start] == '"') start++;
191
if (name[end-1] == '"') end--;
192
name = name.substr(start,end-start);
193
// rename spaces and tabs
194
for (size_t i = 0; i < name.size(); i++) {
195
if (name[i] == ' ' || name[i] == '\t') name[i] = '_';
197
phys[dom_idx] = name.substr(start,end-start);
199
check_macro (scatch(in,"$EndPhysicalNames"), "gmsh input error ($EndPhysicalNames not found).");
200
check_macro (scatch(in,"$"), "gmsh input error: no more label found.");
206
check_macro (label == "Nodes", "$Nodes not found in gmsh file");
209
vector<point> node(nnode);
210
Float infty = numeric_limits<Float>::max();
211
point xmin ( infty, infty, infty);
212
point xmax (-infty, -infty, -infty);
213
for (size_type k = 0; k < nnode; k++) {
215
in >> dummy >> node[k];
216
for (size_type j = 0 ; j < 3; j++) {
217
xmin[j] = min(node[k][j], xmin[j]);
218
xmax[j] = max(node[k][j], xmax[j]);
222
// dimension is deduced from bounding box
225
if (xmax[2] == xmin[2]) {
227
if (xmax[1] == xmin[1]) {
229
if (xmax[0] == xmin[0]) dim = 0;
235
check_macro (scatch(in,"$Elements"), "$Elements not found in gmsh file");
238
vector<size_t> element_gmshtype(nelement);
239
vector<char> element_name (nelement);
240
vector<size_t> element_order (nelement);
241
vector<vector<size_t> > element (nelement);
242
vector<size_t> node_subgeo_variant (nnode, size_t(-1));
243
map<size_t, list<size_t> > point_domain_map, edge_domain_map, face_domain_map, volume_domain_map;
244
size_type map_dim = 0;
246
size_type size_by_dim [4] = {0,0,0,0};
247
size_type size_by_variant [reference_element::max_variant] = {0,0,0,0,0,0,0};
248
for (size_type i = 0; i < nelement; i++) {
249
size_type id, dummy, gmshtype;
250
in >> id >> gmshtype;
251
size_type n_tag_gmsh;
253
size_type domain_idx = 0;
254
for (size_type j = 0 ; j < n_tag_gmsh; j++) {
255
// the first tag is the physical domain index
256
// the second tag is the object index, defined for all elements
257
// the third is zero (in all examples)
261
domain_idx = tag_dummy;
264
check_macro (gmshtype < gmshtype_max,
265
"element #" << id << ": unexpected gmsh type '" << gmshtype << "'");
266
check_macro (gmsh_table[gmshtype].supported,
267
"element #" << id << ": unsupported gmsh type '" << gmshtype << "'");
269
element_name [i] = gmsh_table[gmshtype].variant;
270
element_order[i] = gmsh_table[gmshtype].order;
271
size_type nv = gmsh_table[gmshtype].nv;
272
size_type nn = gmsh_table[gmshtype].nn_tot;
273
size_type dim_i = element_dimension(element_name[i]);
274
size_by_dim[dim_i]++;
275
size_t variant = reference_element::variant(element_name[i]);
276
size_by_variant [variant]++;
277
map_dim = max(map_dim,dim_i);
278
order = max(order,element_order[i]);
279
element[i].resize(nn);
280
for (size_type j = 0; j < nn; j++) {
284
for (size_type subgeo_variant = 0; subgeo_variant <= variant; subgeo_variant++) { // set node dimension
285
for (size_type loc_inod = reference_element::first_inod_by_variant(variant,element_order[i],subgeo_variant),
286
loc_nnod = reference_element:: last_inod_by_variant(variant,element_order[i],subgeo_variant);
287
loc_inod < loc_nnod; loc_inod++) {
288
node_subgeo_variant [element[i][loc_inod]] = subgeo_variant;
291
msh2geo_node_renum (element[i], element_name[i], element_order[i]);
292
if (domain_idx != 0) {
294
case 0: point_domain_map[domain_idx].push_back(i); break;
295
case 1: edge_domain_map[domain_idx].push_back(i); break;
296
case 2: face_domain_map[domain_idx].push_back(i); break;
297
default: volume_domain_map[domain_idx].push_back(i); break;
301
// -------------------------------------------------
303
// -------------------------------------------------
304
// permut: first vertex, then edge, faces and inernal nodes, by increasing subgeo_dim
305
vector<size_t> old2new_inode (nnode, size_t(-1));
307
for (size_type subgeo_variant = 0; subgeo_variant < reference_element::last_variant_by_dimension(map_dim); subgeo_variant++) {
308
for (size_type k = 0; k < nnode; k++) {
309
if (node_subgeo_variant [k] != subgeo_variant) continue;
310
old2new_inode[k] = inode++;
314
vector<size_t> new2old_inode (nnode, size_t(-1));
315
for (size_type k = 0; k < nnode; k++) {
316
new2old_inode[old2new_inode[k]] = k;
318
// ----------------------------------
320
// ----------------------------------
321
size_type version = 4;
322
static const char* geo_variant_name [reference_element::max_variant] = {
331
out << setprecision(numeric_limits<Float>::digits10)
337
<< " dimension\t" << dim << endl;
338
if (sys_coord_name != "cartesian") {
339
out << " coordinate_system " << sys_coord_name << endl;
342
out << " order\t\t" << order << endl;
344
out << " nodes\t\t" << nnode << endl;
347
for (size_type variant = reference_element::first_variant_by_dimension(map_dim);
348
variant < reference_element:: last_variant_by_dimension(map_dim); variant++) {
349
if (size_by_variant[variant] > 0) {
350
out << " " << geo_variant_name [variant] << "\t" << size_by_variant[variant] << endl;
354
out << "end header" << endl
356
for (size_type k = 0; k < nnode; k++) {
357
node[new2old_inode[k]].put (out, dim);
361
for (size_type variant = reference_element::first_variant_by_dimension(map_dim);
362
variant < reference_element:: last_variant_by_dimension(map_dim); variant++) {
363
for (size_type i = 0; i < nelement; i++) {
364
if (reference_element::variant(element_name[i]) != variant) continue;
365
if (element_name[i] != 'e' || element_order[i] > 1) {
366
out << element_name[i] << "\t";
368
if (element_order[i] > 1) {
369
out << "p" << element_order[i] << " ";
371
for (size_type iloc = 0, nloc = element[i].size(); iloc < nloc; iloc++) {
372
out << old2new_inode[element[i][iloc]];
373
if (iloc != nloc - 1) out << " ";
379
// ----------------------------------
381
// ----------------------------------
382
put_domain (out, map_dim, 0, point_domain_map, phys, element_name, element_order, element, old2new_inode);
383
put_domain (out, map_dim, 1, edge_domain_map, phys, element_name, element_order, element, old2new_inode);
384
put_domain (out, map_dim, 2, face_domain_map, phys, element_name, element_order, element, old2new_inode);
385
put_domain (out, map_dim, 3, volume_domain_map, phys, element_name, element_order, element, old2new_inode);
389
cerr << "msh2geo: usage:" << endl
390
<< "msh2geo [-rz|-zr] < in.msh | geo -upgrade - > out.geo" << endl;
393
int main (int argc, char**argv) {
394
// ----------------------------
395
// scan the command line
396
// ----------------------------
397
std::string sys_coord_name = "cartesian";
398
for (int i = 1; i < argc; i++) {
399
if (strcmp (argv[i], "-rz") == 0) sys_coord_name = "rz";
400
else if (strcmp (argv[i], "-zr") == 0) sys_coord_name = "zr";
402
cerr << "field: unknown option `" << argv[i] << "'" << endl;
406
msh2geo (cin, cout, sys_coord_name);
added
removed
nfem/pbin/msh2geo.cc
