geo_seq_get_vtk.cc
Go to the documentation of this file.
//
// geo: vtk input
//
// author: Pierre.Saramito@imag.fr
//
// 5 march 2012
//
#include "rheolef/geo.h"
#include "rheolef/rheostream.h"
#include "rheolef/iorheo.h"
#include "vtk_cell_type.h"
namespace rheolef {
using namespace std;
template <class T>
idiststream&
geo_get_vtk (idiststream& ips, geo_basic<T,sequential>& omega)
{
using namespace std;
typedef typename geo_basic<T,sequential>::size_type size_type;
typedef typename geo_basic<T,sequential>::node_type node_type;
check_macro (ips.good(), "bad input stream for vtk");
istream& is = ips.is();
geo_header hdr;
hdr.order = 1;
// ------------------------------------------------------------------------
// 1) load the coordinates
// POINTS <np> float
// {xi yi zi} i=0..np-1
// ------------------------------------------------------------------------
check_macro (scatch(is,"POINTS"), "unexpected vtk input file");
size_type nnod;
is >> nnod;
if (nnod == 0) {
warning_macro("empty vtk mesh file");
return ips;
}
disarray<node_type, sequential> node (nnod);
scatch (is,"float");
node.get_values (ips);
check_macro (ips.good(), "bad input stream for vtk");
// ------------------------------------------------------------------------
// 1.2) compute dimension
// ------------------------------------------------------------------------
point xmin, xmax;
for (size_type k = 0; k < 3; ++k) {
xmin[k] = std::numeric_limits<T>::max();
xmax[k] = std::numeric_limits<T>::min();
}
for (size_t i = 0, n = node.size(); i < n; ++i) {
for (size_type k = 0; k < 3; ++k) {
xmin[k] = std::min (node[i][k], xmin[k]);
xmax[k] = std::max (node[i][k], xmax[k]);
}
}
hdr.dimension = 3;
for (size_type k = 0; k < 3; ++k) {
if (fabs(xmax[k] - xmin[k]) < std::numeric_limits<T>::epsilon())
hdr.dimension--;
}
// ------------------------------------------------------------------------
// 2) load the domain connectivity
// LINES <ne> <ncs>
// 2 {s1 s2} j=0..ne-1
// or
// POLYGONS <ne> <ncs>
// p {s1 .. sp} j=0..ne-1
// ------------------------------------------------------------------------
// 2.1) get connectivity header
string mark;
is >> ws >> mark;
if (mark == "VERTICES") { // skip unused paragraph in polydata v3.0
size_type n, n2, d, idx;
is >> n >> n2;
for (size_type i = 0; i < n; i++) {
is >> d >> idx;
}
}
do {
if (mark == "POLYGONS" || mark == "LINES" || mark == "CELLS") break;
} while (is >> ws >> mark);
bool have_vtk_cell_type = false;
string mesh_fmt = mark;
if (mesh_fmt == "LINES") hdr.map_dimension = 1;
else if (mesh_fmt == "POLYGONS") hdr.map_dimension = 2;
else if (mesh_fmt == "CELLS") { hdr.map_dimension = 0; have_vtk_cell_type = true; } // yet unknown
else error_macro ("geo: unexpected `" << mesh_fmt << "' in vtk input."
<< " Expect LINES, POLYGONS or CELLS.");
// note: when mesh_fmt == "CELLS", we will also get CELL_TYPES and deduce the map_dimension
size_type ne, ncs;
is >> ne >> ncs;
//
// 2.2) get all connectivity in a flat trash array
std::vector<size_type> cell_trash (ncs);
for (size_type ics = 0; ics < ncs; ics++) {
is >> cell_trash [ics];
}
std::vector<size_type> vtk_cell_type (ne);
if (have_vtk_cell_type) {
is >> ws >> mark;
check_macro (mark == "CELL_TYPES", "geo: unexpected `" << mark << "' in vtk input."
<< " Expect CELL_TYPES.");
size_type ne2;
is >> ne2;
check_macro (ne == ne2, "geo: unexpected CELL_TYPES size=`" << ne2 << "' , expect size="<<ne);
for (size_type ie = 0; ie < ne; ie++) {
is >> vtk_cell_type [ie];
}
}
typedef geo_element_auto<> element_t;
typedef disarray<element_t,sequential> disarray_t;
disarray_t elt (ne);
for (size_type ie = 0, ics = 0; ie < ne; ie++) {
size_type nloc = cell_trash [ics];
ics++;
size_type variant = (!have_vtk_cell_type) ?
reference_element::variant (nloc, hdr.map_dimension) :
vtk_cell_type2variant (vtk_cell_type [ie]);
elt[ie].reset (variant, hdr.order);
for (size_type iloc = 0 ; iloc < nloc; iloc++, ics++) {
elt[ie][iloc] = cell_trash [ics];
}
hdr.dis_size_by_dimension [elt[ie].dimension()]++;
hdr.dis_size_by_variant [elt[ie].variant()]++;
hdr.map_dimension = std::max (hdr.map_dimension, elt[ie].dimension());
}
hdr.dis_size_by_dimension [0] = nnod;
hdr.dis_size_by_variant [0] = nnod;
// ------------------------------------------------------------------------
// 2) split elements by variants
// ------------------------------------------------------------------------
std::array<disarray_t, reference_element::max_variant> tmp_geo_element;
if (hdr.map_dimension > 0) {
for (size_type variant = reference_element::first_variant_by_dimension(hdr.map_dimension);
variant < reference_element:: last_variant_by_dimension(hdr.map_dimension); variant++) {
geo_element_auto<> init_val (variant, hdr.order);
tmp_geo_element [variant] = disarray_t (hdr.dis_size_by_variant [variant], init_val);
}
std::array<size_type, 4> counter;
counter.fill(0);
for (size_type ie = 0; ie < ne; ie++) {
size_type variant = elt[ie].variant();
size_type ige = counter[variant];
tmp_geo_element [variant] [ige] = elt[ie];
counter[variant]++;
}
}
// ------------------------------------------------------------------------
// 3) build & upgrade
// ------------------------------------------------------------------------
omega.build_from_data (hdr, node, tmp_geo_element, true);
return ips;
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
template idiststream& geo_get_vtk<Float> (idiststream&, geo_basic<Float,sequential>&);
}// namespace rheolef
rheolef::reference_element::last_variant_by_dimension
static variant_type last_variant_by_dimension(size_type dim)
Definition: reference_element.h:150
vtk_cell_type.h
rheolef::geo_basic
generic mesh with rerefence counting
Definition: domain_indirect.h:64
warning_macro
#define warning_macro(message)
Definition: dis_macros.h:53
check_macro
check_macro(expr1.have_homogeneous_space(Xh1), "dual(expr1,expr2); expr1 should have homogeneous space. HINT: use dual(interpolate(Xh, expr1),expr2)")
rheolef::geo_header::dimension
size_type dimension
Definition: geo_header.h:39
rheolef::geo_get_vtk
idiststream & geo_get_vtk(idiststream &ips, geo_basic< T, sequential > &omega)
Definition: geo_seq_get_vtk.cc:39
rheolef::geo_element_auto
Definition: geo_element.h:99
rheolef::size_type
size_t size_type
Definition: basis_get.cc:76
rheolef::vtk_cell_type2variant
size_t vtk_cell_type2variant(size_t vtk_cell_type)
Definition: vtk_cell_type.cc:34
rheolef::reference_element::first_variant_by_dimension
static variant_type first_variant_by_dimension(size_type dim)
Definition: reference_element.h:148
mkgeo_ball.variant
variant
Definition: mkgeo_ball.sh:149
rheolef::geo_basic< T, sequential >
Definition: geo.h:1108
rheolef::reference_element::variant
variant_type variant() const
Definition: reference_element.h:99
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::geo_header::map_dimension
size_type map_dimension
Definition: geo_header.h:40
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37
error_macro
#define error_macro(message)
Definition: dis_macros.h:49
point
see the point page for the full documentation
rheolef::disarray< node_type, sequential >
mkgeo_ball.d
d
Definition: mkgeo_ball.sh:154
rheolef::space_numbering::nnod
size_type nnod(const basis_basic< T > &b, const geo_size &gs, size_type map_dim)
Definition: space_numbering.cc:54
rheolef::geo_header::dis_size_by_dimension
size_type dis_size_by_dimension[4]
Definition: geo_header.h:44
mkgeo_ball.n
n
Definition: mkgeo_ball.sh:150
rheolef::geo_header::order
size_type order
Definition: geo_header.h:42
size_type
field::size_type size_type
Definition: branch.cc:425
epsilon
Float epsilon
Definition: transmission_error.cc:25
rheolef::geo_get_vtk< Float >
template idiststream & geo_get_vtk< Float >(idiststream &, geo_basic< Float, sequential > &)
rheolef::std
Definition: vec_expr_v2.h:391
rheolef::geo_header
Definition: geo_header.h:32
rheolef::geo_header::dis_size_by_variant
size_type dis_size_by_variant[reference_element::max_variant]
Definition: geo_header.h:43
rheolef::geo_basic< T, sequential >::build_from_data
void build_from_data(const geo_header &hdr, const disarray< node_type, sequential > &node, std::array< disarray< geo_element_auto<>, sequential >, reference_element::max_variant > &tmp_geo_element, bool do_upgrade)