msh2geo.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 msh2geo unix command
// author: Pierre.Saramito@imag.fr
namespace rheolef {
} // namespace rheolef
#include "rheolef/compiler.h" // after index_set to avoid boost
#include "scatch.icc"
#include "index_set_header.icc"
#include "index_set_body.icc"
#include <array>
#include <cstring>
namespace rheolef {
template<class T>
struct point_basic: std::array<T,3> {
typedef std::array<T,3> base;
point_basic(T x0=T(), T x1=T(), T x2=T())
: base() {
base::operator[](0) = x0;
base::operator[](1) = x1;
base::operator[](2) = x2;
}
};
using point = point_basic<double>;
namespace edge {
#include "edge.icc"
} // namespace edge
namespace triangle {
#include "triangle.icc"
} // namespace triangle
namespace quadrangle {
#include "quadrangle.icc"
} // namespace quadrangle
namespace tetrahedron {
#include "tetrahedron.icc"
} // namespace tetrahedron
namespace prism {
#include "prism.icc"
} // namespace prism
namespace hexahedron {
#include "hexahedron.icc"
} // namespace hexahedron
} // namespace rheolef
#include "reference_element_aux.icc"
namespace rheolef {
// TODO: move somewhere in reference_element_xxx.cc
const char
_reference_element_name [reference_element__max_variant] = {
'p',
'e',
't',
'q',
'T',
'P',
'H'
};
} // namespace rheolef
namespace rheolef {
typedef int orientation_type;
typedef int shift_type;
// side (edge & face): element with reduced node list when order > 1, for sides (edges or faces)
struct geo_element_side: std::array<size_t,4> {
void setname (char name) { _variant = reference_element_variant(name); }
void setindex (size_t index) { _index = index; }
char name() const { return _reference_element_name [_variant]; }
size_t index() const { return _index; }
size_t variant() const { return _variant; }
size_t dimension()const { return reference_element_dimension_by_variant(variant()); }
size_t n_vertex() const { return reference_element_n_vertex (variant()); }
size_t size() const { return n_vertex(); }
geo_element_side()
: array<size_t,4>(), _variant(-1), _index(-1)
{}
protected:
size_t _variant;
size_t _index;
};
// side index with orient and shift for geo_element
struct geo_element_indirect {
typedef int orientation_type;
typedef int shift_type;
void setindex (size_t index) { _index = index; }
void setorientation (size_t orient) { _orient = orient; }
void setshift (size_t shift) { _shift = shift; }
size_t index() const { return _index; }
int orientation() const { return _orient; }
size_t shift() const { return _shift; }
geo_element_indirect()
: _index(-1), _orient(1), _shift(0) {}
protected:
size_t _index;
int _orient;
int _shift;
};
// element with full node list when order > 1
struct geo_element: vector<size_t> {
void setname (char name) { _name = name; }
void setorder (size_t order) { _order = order; }
void setindex (size_t index) { _index = index; }
void setgmshtype (size_t gmshtype) { _gmshtype = gmshtype; }
geo_element_indirect& edge_indirect (size_t i) { return _edge[i]; }
geo_element_indirect& face_indirect (size_t i) { return _face[i]; }
char name() const { return _name; }
size_t order() const { return _order; }
size_t index() const { return _index; }
size_t gmshtype() const { return _gmshtype; }
size_t variant() const { return reference_element_variant (name()); }
size_t dimension()const { return reference_element_dimension_by_name(name()); }
size_t n_vertex() const { return reference_element_n_vertex (variant()); }
const geo_element_indirect& edge_indirect (size_t i) const { return _edge[i]; }
const geo_element_indirect& face_indirect (size_t i) const { return _face[i]; }
size_t n_edge() const { return _reference_element_n_edge_by_variant [variant()]; }
size_t edge (size_t i) const { return _edge[i].index(); }
size_t edge_local_vertex(size_t iedg, size_t edg_iv) const;
size_t face (size_t i) const { return _face[i].index(); }
size_t n_face() const { return _reference_element_n_face_by_variant [variant()]; }
size_t face_n_vertex(size_t loc_ifac) const;
size_t face_local_vertex(size_t iedg, size_t edg_iv) const;
geo_element()
: vector<size_t>(), _name('z'), _order(-1), _index(-1), _gmshtype(-1), _edge(), _face()
{}
protected:
char _name; // TODO: store variant instead of name
size_t _order;
size_t _index;
size_t _gmshtype;
array<geo_element_indirect,12> _edge;
array<geo_element_indirect,6> _face;
};
size_t
geo_element::edge_local_vertex (size_t iedg, size_t edg_iv) const
{
switch (variant()) {
case reference_element__e: return vector<size_t>::operator[] (edg_iv);
case reference_element__t: return triangle::edge [iedg][edg_iv%2];
case reference_element__q: return quadrangle::edge [iedg][edg_iv%2];
case reference_element__T: return tetrahedron::edge [iedg][edg_iv%2];
case reference_element__P: return prism::edge [iedg][edg_iv%2];
case reference_element__H: return hexahedron::edge [iedg][edg_iv%2];
default: error_macro ("invalid variant"); return 0;
}
}
size_t
geo_element::face_local_vertex (size_t loc_ifac, size_t fac_iv) const
{
switch (variant()) {
case reference_element__t:
case reference_element__q: return vector<size_t>::operator[] (fac_iv);
case reference_element__T: return tetrahedron::face [loc_ifac][fac_iv%3];
case reference_element__P: return prism::face [loc_ifac][fac_iv%4];
case reference_element__H: return hexahedron::face [loc_ifac][fac_iv%4];
default: error_macro ("invalid variant"); return 0;
}
}
size_t
geo_element::face_n_vertex(size_t loc_ifac) const
{
switch (variant()) {
case reference_element__t: return 3;
case reference_element__q: return 4;
case reference_element__T: return 3;
case reference_element__H: return 4;
case reference_element__P: return (loc_ifac < 2) ? 3 : 4;
default: error_macro ("invalid variant"); return 0;
}
}
} // namespace rheolef
#include "msh2geo_defs.icc"
#include "msh2geo_node_renum.icc"
#include "geo_element_aux.icc"
using namespace std;
using namespace rheolef;
// necessite de numeroter les aretes
// puis de le stoker dans ts les triangles
// le num arete et son orient
// => doit etre appeler apres la num des aretes
void
numbering_Pk_dis_idof (
size_t order,
const array<size_t,reference_element__max_variant>& size_by_variant,
const array<size_t,reference_element__max_variant+1>& first_by_variant,
const array<size_t,reference_element__max_variant>& loc_ndof_by_variant,
const geo_element& K,
vector<size_t>& dis_idof1)
{
dis_idof1.resize (reference_element_n_node (K.variant(), order));
std::fill (dis_idof1.begin(), dis_idof1.end(), std::numeric_limits<size_t>::max());
for (size_t subgeo_variant = 0, variant = K.variant(); subgeo_variant <= variant; subgeo_variant++) {
size_t subgeo_dim = reference_element_dimension_by_variant (subgeo_variant);
size_t loc_sub_ndof = loc_ndof_by_variant [subgeo_variant];
for (size_t first_loc_idof = reference_element_first_inod_by_variant (variant, order, subgeo_variant),
last_loc_idof = reference_element_last_inod_by_variant (variant, order, subgeo_variant),
loc_idof = first_loc_idof;
loc_idof < last_loc_idof; loc_idof++) {
// 1) local loc_igev on subgeo
size_t loc_igev = (loc_idof - first_loc_idof) / loc_sub_ndof;
size_t loc_igev_j = (loc_idof - first_loc_idof) % loc_sub_ndof;
// 2) then compute ige; computation depends upon the subgeo dimension:
size_t ige;
switch (subgeo_dim) {
case 0: {
// convert node numbering to vertex numbering, for geo order > 1
size_t loc_inod = loc_idof; // only one dof per vertex
ige = K [loc_inod];
break;
}
case 1: {
loc_igev_j = geo_element_fix_edge_indirect (K, loc_igev, loc_igev_j, order);
size_t loc_ige = loc_igev;
ige = K.edge_indirect(loc_ige).index();
break;
}
case 2: {
size_t loc_ige = loc_igev;
if (subgeo_variant == reference_element__t) {
loc_igev_j = geo_element_fix_triangle_indirect (K, loc_igev, loc_igev_j, order);
} else {
size_t loc_ntri = (K.variant() == reference_element__P) ? 2 : 0;
loc_ige += loc_ntri;
loc_igev_j = geo_element_fix_quadrangle_indirect (K, loc_ige, loc_igev_j, order);
}
ige = K.face(loc_ige);
break;
}
case 3: {
ige = K.index();
break;
}
default: {
error_macro ("unexpected subgeo_dim="<<subgeo_dim);
}
}
// 3) then compute igev, by variant
size_t igev = ige;
for (size_t prev_subgeo_variant = reference_element_first_variant_by_dimension(subgeo_dim);
prev_subgeo_variant < subgeo_variant;
prev_subgeo_variant++) {
size_t nprev = size_by_variant [prev_subgeo_variant];
assert_macro (ige >= nprev, "invalid index");
igev -= nprev;
}
size_t dis_igev = igev;
// 4) finally compute dis_idof
size_t dis_idof = 0;
for (size_t prev_subgeo_variant = 0;
prev_subgeo_variant < subgeo_variant;
prev_subgeo_variant++) {
dis_idof += size_by_variant [prev_subgeo_variant]
*loc_ndof_by_variant [prev_subgeo_variant];
}
dis_idof += dis_igev
*loc_ndof_by_variant [subgeo_variant]
+ loc_igev_j;
assert_macro (loc_idof < dis_idof1.size(), "msh2geo: invalid index");
dis_idof1 [loc_idof] = dis_idof;
}
}
}
void
put_domain_noupgrade (
ostream& out,
size_t dim,
size_t dom_dim,
const map<size_t,list<size_t> >& domain_map,
map<size_t, string>& phys,
const vector<geo_element>& element)
{
if (dim == dom_dim && domain_map.size() == 1) return;
for (map<size_t,list<size_t> >::const_iterator
first = domain_map.begin(),
last = domain_map.end(); first != last; first++) {
size_t dom_idx = (*first).first;
const list<size_t>& dom = (*first).second;
string dom_name = phys[dom_idx];
if (dom_name == "") dom_name = "unnamed" + itos(dom_idx);
out << "domain" << endl
<< dom_name << endl
<< "1 " << dom_dim << " " << dom.size() << endl;
for (size_t variant = reference_element_first_variant_by_dimension(dom_dim);
variant < reference_element_last_variant_by_dimension(dom_dim); variant++) {
for (list<size_t>::const_iterator fe = dom.begin(), le = dom.end(); fe != le; fe++) {
size_t ie = *fe;
if (element[ie].variant() != variant) continue;
out << element[ie].name() << "\t";
if (element[ie].order() > 1) {
out << "p" << element[ie].order() << " ";
}
for (size_t j = 0; j < element[ie].size(); j++) {
out << element[ie][j] << " ";
}
out << endl;
}
}
out << endl;
}
}
void
put_domain_upgrade (
ostream& out,
size_t dim,
size_t dom_dim,
const map<size_t,list<size_t> >& domain_map,
map<size_t, string>& phys,
const vector<geo_element>& element,
const vector<geo_element_side>& edge,
const vector<index_set>& edge_ball,
const vector<geo_element_side>& face,
const vector<index_set>& face_ball)
{
if (dim == dom_dim && domain_map.size() == 1) return;
for (map<size_t,list<size_t> >::const_iterator
first = domain_map.begin(),
last = domain_map.end(); first != last; first++) {
size_t dom_idx = (*first).first;
const list<size_t>& dom = (*first).second;
string dom_name = phys[dom_idx];
if (dom_name == "") dom_name = "unnamed" + itos(dom_idx);
out << "domain" << endl
<< dom_name << endl
<< "2 " << dom_dim << " " << dom.size() << endl;
for (size_t variant = reference_element_first_variant_by_dimension(dom_dim);
variant < reference_element_last_variant_by_dimension(dom_dim); variant++) {
for (list<size_t>::const_iterator fe = dom.begin(), le = dom.end(); fe != le; fe++) {
size_t ie = *fe;
if (element[ie].variant() != variant) continue;
size_t isid = 0;
orientation_type orient = 1;
if (element[ie].name() == 'p') {
isid = element[ie][0];
orient = 1;
} else if (element[ie].name() == 'e') {
size_t inod0 = element[ie][0];
size_t inod1 = element[ie][1];
index_set iedge_set = edge_ball[inod0];
iedge_set.inplace_intersection (edge_ball[inod1]);
check_macro (iedge_set.size() == 1,
"msh2geo: error: connectivity problem (iedge.size="<<iedge_set.size()<<")");
isid = *(iedge_set.begin());
orient = (inod0 == edge[isid][0]) ? 1 : -1;
} else if (element[ie].dimension() == 2) {
size_t inod0 = element[ie][0];
index_set iface_set = face_ball[inod0];
for (size_t j = 1; j < element[ie].n_vertex(); ++j) {
size_t inodj = element[ie][j];
iface_set.inplace_intersection (face_ball[inodj]);
}
check_macro (iface_set.size() == 1,
"msh2geo: error: connectivity problem (iface.size="<<iface_set.size()<<")");
isid = *(iface_set.begin());
shift_type shift;
if (element[ie].name() == 't') {
geo_element_get_orientation_and_shift (
face[isid],
element[ie][0],
element[ie][1],
element[ie][2],
orient, shift);
} else {
geo_element_get_orientation_and_shift (
face[isid],
element[ie][0],
element[ie][1],
element[ie][2],
element[ie][3],
orient, shift);
}
} else { // 3d domain
isid = element[ie].index();
orient = 1;
}
out << orient*int(isid) << endl;
}
}
out << endl;
}
}
void msh2geo (istream& in, ostream& out, string sys_coord_name, bool do_upgrade)
{
// ----------------------------------
// 1. input gmsh
// ----------------------------------
// 1.0. preambule
check_macro (scatch(in,"$MeshFormat",true),
"msh2geo: input stream does not contains a gmsh mesh file ($MeshFormat not found).");
double gmsh_fmt_version;
size_t file_type, float_data_size;
in >> gmsh_fmt_version >> file_type >> float_data_size;
check_macro (gmsh_fmt_version <= 2.2,
"gmsh format version " << gmsh_fmt_version << " founded ; expect version 2.2 (HINT: use gmsh -format msh2)");
check_macro (file_type == 0, "msh2geo: unsupported gmsh non-ascii format");
check_macro (scatch(in,"$EndMeshFormat",true), "msh2geo: gmsh input error: $EndMeshFormat not found.");
//
// 1.1 optional domain names
//
bool full_match = true;
bool partial_match = !full_match;
check_macro (scatch(in,"$",partial_match), "msh2geo: gmsh input error: no more label found.");
size_t nphys;
map<size_t, string> phys;
string label;
in >> label;
size_t n_names = 0;
if (label == "PhysicalNames") {
in >> nphys;
for (size_t i = 0; i < nphys; i++) {
string name;
size_t dom_dim, dom_idx;
in >> dom_dim >> dom_idx;
// get name:
char c;
in >> std::ws >> c;
if (c != '"') name.push_back(c);
do {
in.get(c);
name.push_back(c);
} while (c != '"' && c != '\n');
// strip '"' in name
size_t start = 0, end = name.length();
if (name[start] == '"') start++;
if (name[end-1] == '"') end--;
name = name.substr(start,end-start);
// rename spaces and tabs
for (size_t i = 0; i < name.size(); i++) {
if (name[i] == ' ' || name[i] == '\t') name[i] = '_';
}
phys[dom_idx] = name.substr(start,end-start);
}
check_macro (scatch(in,"$EndPhysicalNames",true), "msh2geo: gmsh input error ($EndPhysicalNames not found).");
check_macro (scatch(in,"$",partial_match), "msh2geo: gmsh input error: no more label found.");
in >> label;
}
//
// 1.2. nodes
//
if (label != "Nodes") {
check_macro (scatch(in,"$Nodes",true), "msh2geo: $Nodes not found in gmsh file");
}
size_t nnode;
in >> nnode;
vector<point> node (nnode);
double infty = numeric_limits<double>::max();
point xmin ( infty, infty, infty);
point xmax (-infty, -infty, -infty);
for (size_t k = 0; k < node.size(); k++) {
size_t dummy;
in >> dummy;
for (size_t i = 0; i < 3; ++i) {
in >> node[k][i];
}
for (size_t j = 0 ; j < 3; j++) {
xmin[j] = min(node[k][j], xmin[j]);
xmax[j] = max(node[k][j], xmax[j]);
}
}
//
// dimension is deduced from bounding box
//
size_t dim = 3;
if (xmax[2] == xmin[2]) {
dim = 2;
if (xmax[1] == xmin[1]) {
dim = 1;
if (xmax[0] == xmin[0]) dim = 0;
}
}
//
// 1.3. elements
//
check_macro (scatch(in,"$Elements",true), "msh2geo: $Elements not found in gmsh file");
size_t n_element;
in >> n_element;
vector<geo_element> element (n_element);
vector<size_t> node_subgeo_variant (node.size(), std::numeric_limits<size_t>::max());
array<map<size_t, list<size_t> >,4> domain_map; // domain_map[dim][idom][ie]
size_t map_dim = 0;
size_t order = 0;
for (size_t i = 0; i < element.size(); i++) {
size_t id, dummy, gmshtype;
in >> id >> gmshtype;
size_t n_tag_gmsh;
in >> n_tag_gmsh;
size_t domain_idx = 0;
for (size_t j = 0 ; j < n_tag_gmsh; j++) {
// the first tag is the physical domain index
// the second tag is the object index, defined for all elements
// the third is zero (in all examples)
size_t tag_dummy;
in >> tag_dummy;
if (j == 0) {
domain_idx = tag_dummy;
}
}
check_macro (gmshtype < gmshtype_max,
"msh2geo: element #" << id << ": unexpected gmsh type '" << gmshtype << "'");
check_macro (gmsh_table[gmshtype].supported,
"msh2geo: element #" << id << ": unsupported gmsh type '" << gmshtype << "'");
element[i].setname (gmsh_table[gmshtype].name);
element[i].setorder(gmsh_table[gmshtype].order);
element[i].setgmshtype(gmshtype);
size_t nv = gmsh_table[gmshtype].nv;
size_t nn = gmsh_table[gmshtype].nn_tot;
size_t dim_i = element[i].dimension();
size_t variant = element[i].variant();
map_dim = max(map_dim,dim_i);
if (order == 0) {
order = element[i].order();
} else {
check_macro (order == element[i].order(), "msh2geo: unexpected order="<<element[i].order());
}
element[i].resize(nn);
for (size_t j = 0; j < nn; j++) {
in >> element[i][j];
element[i][j]--;
}
for (size_t subgeo_variant = 0; subgeo_variant <= variant; subgeo_variant++) { // set node dimension
for (size_t loc_inod = reference_element_first_inod_by_variant(variant,element[i].order(),subgeo_variant),
loc_nnod = reference_element_last_inod_by_variant(variant,element[i].order(),subgeo_variant);
loc_inod < loc_nnod; loc_inod++) {
node_subgeo_variant [element[i][loc_inod]] = subgeo_variant;
}
}
// from gmsh to rheolef/geo local node renumbering: element[i][j] are modified
msh2geo_node_renum (element[i], element[i].name(), element[i].order());
if (domain_idx != 0) {
domain_map[dim_i][domain_idx].push_back(i);
}
}
array<size_t,reference_element__max_variant> loc_ndof_by_variant;
reference_element_init_local_nnode_by_variant (order, loc_ndof_by_variant);
// set dis_ie: by increasing variant order, for map_dim element only
size_t last_index = 0;
for (size_t variant = reference_element_first_variant_by_dimension(map_dim);
variant < reference_element_last_variant_by_dimension(map_dim); variant++) {
for (size_t ie = 0; ie < element.size(); ++ie) {
if (element[ie].variant() != variant) continue;
element[ie].setindex (last_index);
switch (element[ie].dimension()) {
case 0: element[ie].resize(1); element[ie][0] = last_index; break;
case 1: element[ie].edge_indirect(0).setindex(last_index); break;
case 2: element[ie].face_indirect(0).setindex(last_index); break;
default: break;
}
last_index++;
}
}
// -------------------------------------------------
// 2. node reordering, first pass
// -------------------------------------------------
// permut: first vertex, then edge, faces and inernal nodes, by increasing subgeo_dim
vector<size_t> old2new_inode (node.size(), std::numeric_limits<size_t>::max());
if (true || !do_upgrade) {
// 2.1. when no upgrade
size_t new_inode = 0;
for (size_t subgeo_variant = 0; subgeo_variant < reference_element_last_variant_by_dimension(map_dim); subgeo_variant++) {
for (size_t old_inode = 0; old_inode < node.size(); old_inode++) {
if (node_subgeo_variant [old_inode] != subgeo_variant) continue;
old2new_inode[old_inode] = new_inode++;
}
}
} else {
// 2.2. when no upgrade: will be renumbered after side creation
for (size_t inode = 0; inode < node.size(); inode++) {
old2new_inode[inode] = inode;
}
}
for (size_t inode = 0; inode < node.size(); inode++) {
check_macro (old2new_inode[inode] != std::numeric_limits<size_t>::max(), "msh2geo: invalid permutation");
}
// inv permut
vector<size_t> new2old_inode (node.size(), std::numeric_limits<size_t>::max());
for (size_t inode = 0; inode < node.size(); inode++) {
new2old_inode[old2new_inode[inode]] = inode;
}
for (size_t inode = 0; inode < node.size(); inode++) {
check_macro (new2old_inode[inode] != std::numeric_limits<size_t>::max(), "msh2geo: invalid inverse permutation for inode="<<inode);
}
// for convenience, apply the new numbering to nodes and elements
{
vector<point> new_node (node.size());
for (size_t old_inode = 0; old_inode < node.size(); old_inode++) {
new_node [old2new_inode[old_inode]] = node [old_inode];
}
for (size_t inode = 0; inode < node.size(); inode++) {
node [inode] = new_node [inode];
}
}
for (size_t ie = 0; ie < element.size(); ++ie) {
for (size_t i = 0; i < element[ie].size(); ++i) {
element[ie][i] = old2new_inode[element[ie][i]];
}
}
// ------------------------------------------------------------------------
// 3.1) compute all edges
// ------------------------------------------------------------------------
vector<index_set> edge_ball (node.size());
vector<geo_element_side> edge;
if (do_upgrade && map_dim >= 2) {
list <geo_element_side> edge_list;
size_t nedg = 0;
for (size_t ie = 0; ie < element.size(); ++ie) {
if (element[ie].dimension() != map_dim) continue;
for (size_t iedg = 0; iedg < element[ie].n_edge(); ++iedg) {
size_t iv0 = element[ie].edge_local_vertex(iedg, 0);
size_t iv1 = element[ie].edge_local_vertex(iedg, 1);
size_t inod0 = element[ie][iv0];
size_t inod1 = element[ie][iv1];
index_set iedge_set = edge_ball[inod0];
iedge_set.inplace_intersection (edge_ball[inod1]);
check_macro (iedge_set.size() <= 1,
"msh2geo: error: connectivity problem (iedge.size="<<iedge_set.size()<<")");
if (iedge_set.size() == 1) continue; // edge already exists
edge_ball[inod0].insert (nedg);
edge_ball[inod1].insert (nedg);
geo_element_side new_edge;
new_edge.setname('e');
new_edge[0] = inod0;
new_edge[1] = inod1;
new_edge.setindex(nedg);
edge_list.push_back (new_edge);
nedg++;
}
}
edge.resize (edge_list.size());
std::copy (edge_list.begin(), edge_list.end(), edge.begin());
}
// propagate edge numbering inside elements
if (do_upgrade && map_dim >= 2) {
for (size_t ie = 0; ie < element.size(); ++ie) {
if (element[ie].dimension() != map_dim) continue;
for (size_t iedg = 0; iedg < element[ie].n_edge(); ++iedg) {
size_t iv0 = element[ie].edge_local_vertex(iedg, 0);
size_t iv1 = element[ie].edge_local_vertex(iedg, 1);
size_t inod0 = element[ie][iv0];
size_t inod1 = element[ie][iv1];
index_set iedge_set = edge_ball[inod0];
iedge_set.inplace_intersection (edge_ball[inod1]);
check_macro (iedge_set.size() == 1,
"msh2geo: error: connectivity problem (iedge.size="<<iedge_set.size()<<")");
size_t iedge = *(iedge_set.begin());
element[ie].edge_indirect(iedg).setindex (iedge);
element[ie].edge_indirect(iedg).setorientation((edge[iedge][0] == inod0) ? 1 : -1);
}
}
}
// ------------------------------------------------------------------------
// 3.2) compute all faces
// ------------------------------------------------------------------------
vector<index_set> face_ball (node.size());
vector<geo_element_side> face;
if (do_upgrade && map_dim >= 3) {
list <geo_element_side> face_list;
size_t nfac = 0;
for (size_t ie = 0; ie < element.size(); ++ie) {
if (element[ie].dimension() != map_dim) continue;
for (size_t loc_ifac = 0; loc_ifac < element[ie].n_face(); ++loc_ifac) {
size_t iv0 = element[ie].face_local_vertex(loc_ifac, 0);
size_t inod0 = element[ie][iv0];
index_set iface_set = face_ball[inod0];
for (size_t j = 1; j < element[ie].face_n_vertex(loc_ifac); ++j) {
size_t ivj = element[ie].face_local_vertex(loc_ifac, j);
size_t inodj = element[ie][ivj];
iface_set.inplace_intersection (face_ball[inodj]);
}
check_macro (iface_set.size() <= 1,
"msh2geo: error: connectivity problem (iface.size="<<iface_set.size()<<")");
if (iface_set.size() == 1) continue; // face already exists
geo_element_side new_face;
char face_name = (element[ie].face_n_vertex(loc_ifac) == 3) ? 't' : 'q';
new_face.setname(face_name);
new_face.setindex(nfac);
for (size_t j = 0; j < element[ie].face_n_vertex(loc_ifac); ++j) {
size_t ivj = element[ie].face_local_vertex(loc_ifac, j);
size_t inodj = element[ie][ivj];
face_ball[inodj].insert (nfac);
new_face[j] = inodj;
}
face_list.push_back (new_face);
nfac++;
}
}
face.resize (face_list.size());
std::copy (face_list.begin(), face_list.end(), face.begin());
}
// propagate face numbering inside elements
if (do_upgrade && map_dim >= 3) {
for (size_t ie = 0; ie < element.size(); ++ie) {
if (element[ie].dimension() != map_dim) continue;
for (size_t loc_ifac = 0; loc_ifac < element[ie].n_face(); ++loc_ifac) {
size_t iv0 = element[ie].face_local_vertex(loc_ifac, 0);
size_t inod0 = element[ie][iv0];
index_set iface_set = face_ball[inod0];
for (size_t j = 1; j < element[ie].face_n_vertex(loc_ifac); ++j) {
size_t ivj = element[ie].face_local_vertex(loc_ifac, j);
size_t inodj = element[ie][ivj];
iface_set.inplace_intersection (face_ball[inodj]);
}
check_macro (iface_set.size() == 1,
"msh2geo: error: connectivity problem (iface.size="<<iface_set.size()<<")");
size_t iface = *(iface_set.begin());
orientation_type orient;
shift_type shift;
if (element[ie].face_n_vertex(loc_ifac) == 3) {
geo_element_get_orientation_and_shift (
face[iface],
element[ie][element[ie].face_local_vertex(loc_ifac, 0)],
element[ie][element[ie].face_local_vertex(loc_ifac, 1)],
element[ie][element[ie].face_local_vertex(loc_ifac, 2)],
orient, shift);
} else {
geo_element_get_orientation_and_shift (
face[iface],
element[ie][element[ie].face_local_vertex(loc_ifac, 0)],
element[ie][element[ie].face_local_vertex(loc_ifac, 1)],
element[ie][element[ie].face_local_vertex(loc_ifac, 2)],
element[ie][element[ie].face_local_vertex(loc_ifac, 3)],
orient, shift);
}
element[ie].face_indirect(loc_ifac).setindex (iface);
element[ie].face_indirect(loc_ifac).setorientation (orient);
element[ie].face_indirect(loc_ifac).setshift (shift);
}
}
}
// ------------------------------------------------------------------------
// 3.4) count all elements
// ------------------------------------------------------------------------
array<size_t,reference_element__max_variant> size_by_variant;
array<size_t,reference_element__max_variant+1> first_by_variant;
size_by_variant.fill (0);
first_by_variant.fill (0);
if (true) {
size_t n_vertex = 0;
for (size_t inode = 0; inode < node.size(); inode++) {
if (node_subgeo_variant [inode] == reference_element__p) n_vertex++;
}
size_by_variant [reference_element__p] = n_vertex;
if (map_dim >= 2 && edge.size() != 0) {
size_by_variant [reference_element__e] = edge.size();
}
if (map_dim >= 3 && face.size() != 0) {
for (size_t loc_ifac = 0; loc_ifac < face.size(); ++loc_ifac) {
size_by_variant [face[loc_ifac].variant()]++;
}
}
for (size_t ie = 0; ie < element.size(); ++ie) {
if (element[ie].dimension() != map_dim) continue;
size_t variant = element[ie].variant();
size_by_variant [variant]++;
}
for (size_t dim = 0; dim <= 3; ++dim) {
for (size_t variant = reference_element_first_variant_by_dimension(dim);
variant < reference_element_last_variant_by_dimension(dim); variant++) {
first_by_variant [variant+1] = size_by_variant [variant];
}
}
}
// -------------------------------------------------
// 4. node reordering, second pass
// -------------------------------------------------
// need edges & faces with index & orient
if (do_upgrade) {
std::fill (old2new_inode.begin(), old2new_inode.end(), std::numeric_limits<size_t>::max());
std::fill (new2old_inode.begin(), new2old_inode.end(), std::numeric_limits<size_t>::max());
for (size_t variant = reference_element_first_variant_by_dimension(map_dim);
variant < reference_element_last_variant_by_dimension(map_dim); variant++) {
for (size_t ie = 0; ie < element.size(); ++ie) {
if (element[ie].variant() != variant) continue;
const geo_element& old_K = element[ie];
std::vector<size_t> new_K;
numbering_Pk_dis_idof (
order, size_by_variant, first_by_variant, loc_ndof_by_variant, old_K, new_K);
for (size_t loc_inod = 0; loc_inod < old_K.size(); loc_inod++) {
size_t old_inod = old_K [loc_inod];
size_t new_inod = new_K [loc_inod];
if (old2new_inode [old_inod] != std::numeric_limits<size_t>::max()) continue; // already done
old2new_inode [old_inod] = new_inod;
}
}
}
} else {
// when no upgrade:
for (size_t inode = 0; inode < node.size(); inode++) {
old2new_inode[inode] = inode;
}
}
for (size_t inode = 0; inode < node.size(); inode++) {
check_macro (old2new_inode[inode] != std::numeric_limits<size_t>::max(), "msh2geo: invalid permutation");
}
// inv permut
for (size_t inode = 0; inode < node.size(); inode++) {
new2old_inode[old2new_inode[inode]] = inode;
}
for (size_t inode = 0; inode < node.size(); inode++) {
check_macro (new2old_inode[inode] != std::numeric_limits<size_t>::max(), "msh2geo: invalid inverse permutation for inode="<<inode);
}
#ifdef TODO
#endif // TODO
// ----------------------------------
// 5. output geo
// ----------------------------------
// 5.1. output the mesh
size_t version = 4;
static const char* geo_variant_name [reference_element__max_variant] = {
"points",
"edges",
"triangles",
"quadrangles",
"tetrahedra",
"prisms",
"hexahedra"
};
out << setprecision(numeric_limits<double>::digits10)
<< "#!geo" << endl
<< endl
<< "mesh" << endl
<< version << endl
<< "header" << endl
<< " dimension\t" << dim << endl;
if (sys_coord_name != "cartesian") {
out << " coordinate_system " << sys_coord_name << endl;
}
if (order != 1) {
out << " order\t\t" << order << endl;
}
out << " nodes\t\t" << node.size() << endl;
if (map_dim > 0) {
for (size_t variant = reference_element_first_variant_by_dimension(map_dim);
variant < reference_element_last_variant_by_dimension(map_dim); variant++) {
if (size_by_variant[variant] > 0) {
out << " " << geo_variant_name [variant] << "\t" << size_by_variant[variant] << endl;
}
}
}
if (map_dim >= 3 && size_by_variant[reference_element__t] != 0) {
out << " triangles\t" << size_by_variant[reference_element__t] << endl;
}
if (map_dim >= 3 && size_by_variant[reference_element__q] != 0) {
out << " quadrangles\t" << size_by_variant[reference_element__q] << endl;
}
if (map_dim >= 2 && size_by_variant[reference_element__e] != 0) {
out << " edges\t" << size_by_variant[reference_element__e] << endl;
}
out << "end header" << endl
<< endl;
// nodes:
for (size_t inode = 0; inode < node.size(); inode++) {
for (size_t i = 0; i < dim; ++i) {
out << node[new2old_inode[inode]][i];
if (i+1 != dim) out << " ";
}
out << endl;
}
out << endl;
// elements:
for (size_t variant = reference_element_first_variant_by_dimension(map_dim);
variant < reference_element_last_variant_by_dimension(map_dim); variant++) {
for (size_t ie = 0; ie < element.size(); ie++) {
if (element[ie].variant() != variant) continue;
if (!do_upgrade) {
if (element[ie].name() != 'e' || element[ie].order() > 1) {
out << element[ie].name() << "\t";
}
if (element[ie].order() > 1) {
out << "p" << element[ie].order() << " ";
}
for (size_t iloc = 0, nloc = element[ie].size(); iloc < nloc; iloc++) {
out << element[ie][iloc];
if (iloc+1 != nloc) out << " ";
}
} else {
// upgrade: truncate inodes from pk to p1
out << element[ie].name() << "\t";
for (size_t iloc = 0, nloc = element[ie].n_vertex(); iloc < nloc; iloc++) {
out << element[ie][iloc];
if (iloc+1 != nloc) out << " ";
}
}
out << endl;
}
}
out << endl;
// faces: no-empty when upgrade & map_dim >= 3
for (size_t variant = reference_element_first_variant_by_dimension(2);
variant < reference_element_last_variant_by_dimension(2); variant++) {
for (size_t ie = 0; ie < face.size(); ie++) {
if (face[ie].variant() != variant) continue;
out << face[ie].name() << "\t";
for (size_t iloc = 0, nloc = face[ie].n_vertex(); iloc < nloc; iloc++) {
out << face[ie][iloc];
if (iloc+1 != nloc) out << " ";
}
out << endl;
}
}
// edges: no-empty when upgrade & map_dim >= 2
for (size_t ie = 0, ne = edge.size(); ie < ne; ++ie) {
out << "e\t" << old2new_inode[edge[ie][0]]
<< " " << old2new_inode[edge[ie][1]] << endl;
}
out << endl;
//
// 5.2. output domains
//
for (size_t d = 0; d <= 3; ++d) {
if (!do_upgrade) {
put_domain_noupgrade (out, map_dim, d, domain_map[d], phys, element);
} else {
put_domain_upgrade (out, map_dim, d, domain_map[d], phys, element,
edge, edge_ball, face, face_ball);
}
}
}
void usage()
{
cerr << "msh2geo: usage:" << endl
<< "msh2geo [-[no]upgrade][-rz|-zr] in.msh > out.geo" << endl
<< "msh2geo [-[no]upgrade][-rz|-zr] < in.msh > out.geo" << endl;
exit (1);
}
int main (int argc, char**argv) {
// ----------------------------
// scan the command line
// ----------------------------
string input_filename = "";
std::string sys_coord_name = "cartesian";
bool do_upgrade = true; // upgrade: put sides & truncate pk inodes to p1
for (int i = 1; i < argc; i++) {
if (strcmp (argv[i], "-rz") == 0) sys_coord_name = "rz";
else if (strcmp (argv[i], "-zr") == 0) sys_coord_name = "zr";
else if (strcmp (argv[i], "-cartesian") == 0) sys_coord_name = "cartesian";
else if (strcmp (argv[i], "-upgrade") == 0) do_upgrade = true;
else if (strcmp (argv[i], "-noupgrade") == 0) do_upgrade = false;
else if (argv[i][0] == '-') {
cerr << "field: unknown option `" << argv[i] << "'" << endl;
usage();
} else {
input_filename = argv[i];
}
}
if (input_filename == "") {
msh2geo (cin, cout, sys_coord_name, do_upgrade);
} else {
ifstream fin (input_filename.c_str());
if (!fin.good()) {
cerr << "msh2geo: unable to read file \""<<input_filename<<"\"" << endl; exit (1);
}
msh2geo (fin, cout, sys_coord_name, do_upgrade);
}
}
edge
see the edge page for the full documentation
hexahedron
see the hexahedron page for the full documentation
point
see the point page for the full documentation
prism
see the prism page for the full documentation
quadrangle
see the quadrangle page for the full documentation
rheolef::geo_element_indirect
Definition: geo_element_indirect.h:32
rheolef::geo_element_indirect::shift
shift_type shift() const
Definition: geo_element_indirect.h:53
rheolef::geo_element_indirect::setorientation
void setorientation(size_t orient)
Definition: msh2geo.cc:215
rheolef::geo_element_indirect::setindex
void setindex(size_t index)
Definition: msh2geo.cc:214
rheolef::geo_element_indirect::shift
size_t shift() const
Definition: msh2geo.cc:219
rheolef::geo_element_indirect::_index
size_t _index
Definition: msh2geo.cc:223
rheolef::geo_element_indirect::_shift
int _shift
Definition: msh2geo.cc:225
rheolef::geo_element_indirect::index
size_t index() const
Definition: msh2geo.cc:217
rheolef::geo_element_indirect::setshift
void setshift(size_t shift)
Definition: msh2geo.cc:216
rheolef::geo_element_indirect::geo_element_indirect
geo_element_indirect()
Definition: msh2geo.cc:220
rheolef::geo_element_indirect::shift_type
int shift_type
Definition: msh2geo.cc:213
rheolef::geo_element_indirect::orientation
int orientation() const
Definition: msh2geo.cc:218
rheolef::geo_element_indirect::_orient
int _orient
Definition: msh2geo.cc:224
rheolef::geo_element_indirect::orientation_type
int orientation_type
Definition: msh2geo.cc:212
rheolef::geo_element_indirect::index
size_type index() const
Definition: geo_element_indirect.h:51
rheolef::geo_element
see the geo_element page for the full documentation
Definition: geo_element.h:102
rheolef::geo_element::face
size_t face(size_t i) const
Definition: msh2geo.cc:247
rheolef::geo_element::_edge
array< geo_element_indirect, 12 > _edge
Definition: msh2geo.cc:259
rheolef::geo_element::edge_indirect
const geo_element_indirect & edge_indirect(size_t i) const
Definition: msh2geo.cc:242
rheolef::geo_element::variant
size_t variant() const
Definition: msh2geo.cc:239
rheolef::geo_element::edge_indirect
geo_element_indirect & edge_indirect(size_t i)
Definition: msh2geo.cc:233
rheolef::geo_element::setindex
void setindex(size_t index)
Definition: msh2geo.cc:231
rheolef::geo_element::face_indirect
geo_element_indirect & face_indirect(size_t i)
Definition: msh2geo.cc:234
rheolef::geo_element::edge
size_t edge(size_t i) const
Definition: msh2geo.cc:245
rheolef::geo_element::face
size_type face(size_type i) const
Definition: geo_element.h:210
rheolef::geo_element::_index
size_t _index
Definition: msh2geo.cc:257
rheolef::geo_element::size
size_type size() const
Definition: geo_element.h:168
rheolef::geo_element::index
size_t index() const
Definition: msh2geo.cc:237
rheolef::geo_element::setorder
void setorder(size_t order)
Definition: msh2geo.cc:230
rheolef::geo_element::dimension
size_t dimension() const
Definition: msh2geo.cc:240
rheolef::geo_element::setname
void setname(char name)
Definition: msh2geo.cc:229
rheolef::geo_element::_gmshtype
size_t _gmshtype
Definition: msh2geo.cc:258
rheolef::geo_element::edge_local_vertex
size_t edge_local_vertex(size_t iedg, size_t edg_iv) const
Definition: msh2geo.cc:263
rheolef::geo_element::order
size_t order() const
Definition: msh2geo.cc:236
rheolef::geo_element::n_face
size_t n_face() const
Definition: msh2geo.cc:248
rheolef::geo_element::setgmshtype
void setgmshtype(size_t gmshtype)
Definition: msh2geo.cc:232
rheolef::geo_element::_face
array< geo_element_indirect, 6 > _face
Definition: msh2geo.cc:260
rheolef::geo_element::face_local_vertex
size_t face_local_vertex(size_t iedg, size_t edg_iv) const
Definition: msh2geo.cc:276
rheolef::geo_element::_name
char _name
Definition: msh2geo.cc:255
rheolef::geo_element::face_indirect
const geo_element_indirect & face_indirect(size_t i) const
Definition: msh2geo.cc:243
rheolef::geo_element::face_n_vertex
size_t face_n_vertex(size_t loc_ifac) const
Definition: msh2geo.cc:288
rheolef::geo_element::variant
variant_type variant() const
Definition: geo_element.h:161
rheolef::geo_element::_order
size_t _order
Definition: msh2geo.cc:256
rheolef::geo_element::geo_element
geo_element()
Definition: msh2geo.cc:251
rheolef::geo_element::n_edge
size_t n_edge() const
Definition: msh2geo.cc:244
rheolef::geo_element::edge_indirect
const geo_element_indirect & edge_indirect(size_type i) const
Definition: geo_element.h:193
rheolef::geo_element::gmshtype
size_t gmshtype() const
Definition: msh2geo.cc:238
rheolef::geo_element::name
char name() const
Definition: geo_element.h:169
rheolef::geo_element::n_vertex
size_t n_vertex() const
Definition: msh2geo.cc:241
rheolef::geo_element::order
size_type order() const
Definition: geo_element.h:162
rheolef::index_set
Definition: index_set_header.icc:58
rheolef::index_set::inplace_intersection
void inplace_intersection(const index_set &b)
Definition: index_set_body.icc:129
rheolef::point_basic
Definition: point.h:87
rheolef::point_basic::point_basic
point_basic(T x0=T(), T x1=T(), T x2=T())
Definition: msh2geo.cc:136
rheolef::point_basic::base
std::array< T, 3 > base
Definition: msh2geo.cc:135
rheolef::std
Definition: vec_expr_v2.h:402
rheolef::reference_element__H
constexpr size_t reference_element__H
Definition: reference_element_aux.icc:39
rheolef::reference_element__e
constexpr size_t reference_element__e
Definition: reference_element_aux.icc:34
rheolef::reference_element__P
constexpr size_t reference_element__P
Definition: reference_element_aux.icc:38
rheolef::_reference_element_n_face_by_variant
static size_t _reference_element_n_face_by_variant[reference_element__max_variant]
Definition: reference_element_aux.icc:54
rheolef::reference_element__T
constexpr size_t reference_element__T
Definition: reference_element_aux.icc:37
rheolef::shift_type
int shift_type
Definition: msh2geo.cc:191
rheolef::orientation_type
int orientation_type
Definition: msh2geo.cc:190
rheolef::reference_element__t
constexpr size_t reference_element__t
Definition: reference_element_aux.icc:35
rheolef::reference_element__q
constexpr size_t reference_element__q
Definition: reference_element_aux.icc:36
rheolef::reference_element__max_variant
constexpr size_t reference_element__max_variant
Definition: reference_element_aux.icc:40
rheolef::_reference_element_n_edge_by_variant
static size_t _reference_element_n_edge_by_variant[reference_element__max_variant]
Definition: reference_element_aux.icc:53
rheolef::_reference_element_name
const char _reference_element_name[reference_element__max_variant]
Definition: field2gmsh_pos.cc:78
tetrahedron
see the tetrahedron page for the full documentation
triangle
see the triangle page for the full documentation
assert_macro
#define assert_macro(ok_condition, message)
Definition: dis_macros.h:113
error_macro
#define error_macro(message)
Definition: dis_macros.h:49
edge.icc
edge - reference element
dimension
const size_t dimension
Definition: edge.icc:64
n_vertex
const size_t n_vertex
Definition: edge.icc:66
T
Expr1::float_type T
Definition: field_expr.h:261
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
geo_element_aux.icc
hexahedron.icc
hexahedron - reference element
face
const size_t face[n_face][4]
Definition: hexahedron.icc:110
edge
const size_t edge[n_edge][2]
Definition: hexahedron.icc:118
index_set_body.icc
index_set_header.icc
numbering_Pk_dis_idof
void numbering_Pk_dis_idof(size_t order, const array< size_t, reference_element__max_variant > &size_by_variant, const array< size_t, reference_element__max_variant+1 > &first_by_variant, const array< size_t, reference_element__max_variant > &loc_ndof_by_variant, const geo_element &K, vector< size_t > &dis_idof1)
Definition: msh2geo.cc:313
usage
void usage()
Definition: msh2geo.cc:1051
main
int main(int argc, char **argv)
Definition: msh2geo.cc:1058
put_domain_noupgrade
void put_domain_noupgrade(ostream &out, size_t dim, size_t dom_dim, const map< size_t, list< size_t > > &domain_map, map< size_t, string > &phys, const vector< geo_element > &element)
Definition: msh2geo.cc:396
msh2geo
void msh2geo(istream &in, ostream &out, string sys_coord_name, bool do_upgrade)
Definition: msh2geo.cc:513
put_domain_upgrade
void put_domain_upgrade(ostream &out, size_t dim, size_t dom_dim, const map< size_t, list< size_t > > &domain_map, map< size_t, string > &phys, const vector< geo_element > &element, const vector< geo_element_side > &edge, const vector< index_set > &edge_ball, const vector< geo_element_side > &face, const vector< index_set > &face_ball)
Definition: msh2geo.cc:434
msh2geo_defs.icc
msh2geo_node_renum.icc
mkgeo_ball.variant
int variant
Definition: mkgeo_ball.sh:149
mkgeo_ball.c
int c
Definition: mkgeo_ball.sh:153
mkgeo_ball.d
int d
Definition: mkgeo_ball.sh:154
mkgeo_ball.order
order
Definition: mkgeo_ball.sh:343
mkgeo_ball.map_dim
map_dim
Definition: mkgeo_ball.sh:337
mkgeo_ball.dim
int dim
Definition: mkgeo_ball.sh:307
mkgeo_contraction.name
string name
Definition: mkgeo_contraction.sh:133
mkgeo_couette.nn
int nn
Definition: mkgeo_couette.sh:71
rheolef-config.version
version
Definition: rheolef-config.in:126
rheolef::io::out
@ out
Definition: rheostream.h:167
rheolef::space_constant::vector
@ vector
Definition: space_constant.h:137
rheolef::space_numbering::dis_idof
void dis_idof(const basis_basic< T > &b, const geo_size &gs, const geo_element &K, typename std::vector< size_type >::iterator dis_idof_tab)
Definition: space_numbering.cc:187
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
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
rheolef::itos
std::string itos(std::string::size_type i)
itos: see the rheostream page for the full documentation
rheolef::msh2geo_node_renum
void msh2geo_node_renum(vector< size_t > &element, size_t variant, size_t order)
Definition: msh2geo_node_renum.icc:447
prism.icc
prism - reference element
quadrangle.icc
quadrangle - reference element
scatch.icc
rheolef::const_iterator
Definition: field_expr_recursive.h:552
tetrahedron.icc
tetrahedron - reference element
triangle.icc
triangle - reference element