Loading...
Searching...
No Matches
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" + std::to_string(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" + std::to_string(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);
}
}
point
see the point 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::shift_type
short int shift_type
Definition geo_element_indirect.h:39
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::orientation
int orientation() const
Definition msh2geo.cc:218
rheolef::geo_element_indirect::orientation_type
short int orientation_type
Definition geo_element_indirect.h:38
rheolef::geo_element_indirect::_orient
int _orient
Definition msh2geo.cc:224
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::face_indirect
const geo_element_indirect & face_indirect(size_t i) const
Definition msh2geo.cc:243
rheolef::geo_element::variant
size_t variant() const
Definition msh2geo.cc:239
rheolef::geo_element::setindex
void setindex(size_t index)
Definition msh2geo.cc:231
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::edge_indirect
geo_element_indirect & edge_indirect(size_t i)
Definition msh2geo.cc:233
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::edge_indirect
const geo_element_indirect & edge_indirect(size_t i) const
Definition msh2geo.cc:242
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::edge_indirect
const geo_element_indirect & edge_indirect(size_type i) const
Definition geo_element.h:193
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_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::face_indirect
geo_element_indirect & face_indirect(size_t i)
Definition msh2geo.cc:234
rheolef::geo_element::n_edge
size_t n_edge() const
Definition msh2geo.cc:244
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:126
rheolef::point_basic
Definition point.h:86
rheolef::point_basic::point_basic
point_basic(T x0=T(), T x1=T(), T x2=T())
Definition msh2geo.cc:136
rheolef::point_basic< FT >::base
std::array< FT, 3 > base
Definition field2gmsh_pos.cc:60
rheolef::point
point_basic< Float > point
Definition point.h:163
rheolef::shift_type
int shift_type
Definition msh2geo.cc:191
rheolef::orientation_type
int orientation_type
Definition msh2geo.cc:190
rheolef::reference_element__max_variant
constexpr size_t reference_element__max_variant
Definition reference_element_aux.icc:40
rheolef::_reference_element_name
const char _reference_element_name[reference_element__max_variant]
Definition field2gmsh_pos.cc:78
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
main
int main()
Definition field2bb.cc:58
T
Expr1::float_type T
Definition field_expr.h:230
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
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
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
rheolef::edge
Definition reference_element.cc:29
rheolef::hexahedron
Definition reference_element.cc:49
rheolef::hexahedron::face
const size_t face[n_face][4]
Definition reference_element.cc:111
rheolef::hexahedron::edge
const size_t edge[n_edge][2]
Definition reference_element.cc:119
rheolef::prism
Definition reference_element.cc:45
rheolef::prism::face
const size_t face[n_face][4]
Definition reference_element.cc:87
rheolef::prism::edge
const size_t edge[n_edge][2]
Definition reference_element.cc:94
rheolef::quadrangle
Definition reference_element.cc:37
rheolef::quadrangle::edge
const size_t edge[n_edge][2]
Definition reference_element.cc:84
rheolef::tetrahedron
Definition reference_element.cc:41
rheolef::tetrahedron::face
const size_t face[n_face][3]
Definition reference_element.cc:117
rheolef::tetrahedron::edge
const size_t edge[n_edge][2]
Definition reference_element.cc:123
rheolef::triangle
Definition reference_element.cc:33
rheolef::triangle::edge
const size_t edge[n_edge][2]
Definition reference_element.cc:87
rheolef
This file is part of Rheolef.
Definition compiler_eigen.h:39
rheolef::msh2geo_node_renum
void msh2geo_node_renum(vector< size_t > &element, size_t variant, size_t order)
Definition msh2geo_node_renum.icc:447
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:44
std
STL namespace.
d
size_t d
Definition neumann-laplace-lambda.cc:28
prism.icc
prism - reference element
quadrangle.icc
quadrangle - reference element
scatch.icc
rheolef::const_iterator
Definition field_expr_recursive.h:593
tetrahedron.icc
tetrahedron - reference element
triangle.icc
triangle - reference element