field2gmsh_pos.cc
Go to the documentation of this file.
//
// Passage du format field au format bb de Bamg
//
// author: Pierre.Saramito@imag.fr
//
// date: 03/09/1997 ; 04/01/2018
//
/*Prog:field2bb
NAME: @code{field2gmsh_pos} - convert from .field file format to gmsh .pos one
@pindex field2gmsh_pos
@fiindex @file{.pos} gmsh field
@fiindex @file{.field} field
@toindex @code{gmsh}
SYNOPSIS:
@example
field2gmsh_pos < @var{input}.field > @var{output}.pos
@end example
DESCRIPTION:
Convert a @file{.field} file into a gmsh @file{.pos} one.
The output goes to standart output.
This command is useful for mesh adaptation with @code{gmsh}.
End:*/
#include "rheolef/compiler.h" // after index_set to avoid boost
#include "scatch.icc"
#include <iostream>
#include <iomanip>
#include <string>
#include <vector>
#include <array>
#include <limits>
#include <unistd.h>
using namespace std;
using namespace rheolef;
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 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'
};
struct geo_element: std::array<size_t,8> {
void setname (char name) { _variant = reference_element_variant(name); }
char name() const { return _reference_element_name [_variant]; }
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() : array<size_t,8>(), _variant(-1) {}
protected:
size_t _variant;
};
struct my_geo: vector<geo_element> {
// allocator:
my_geo();
// data:
size_t dimension;
size_t map_dimension;
size_t order;
string sys_coord;
array<size_t,reference_element__max_variant> size_by_variant;
array<size_t,3> size_by_dimension;
vector<point> node;
};
inline
my_geo::my_geo()
: dimension(0),
map_dimension(0),
order(1),
sys_coord("cartesian"),
size_by_variant(),
size_by_dimension(),
node()
{
}
} // namespace rheolef
void
get_geo (istream& in, my_geo& omega)
{
static const char* label_variant [] = {
"nodes", "edges", "triangles", "quadrangles", "tetrahedra", "prisms", "hexahedra" };
check_macro (scatch(in,"\nmesh",true), "input stream does not contains a geo.");
size_t version;
in >> version;
check_macro (version == 4, "mesh format version 4 expected, but format version " << version << " founded");
// get header
check_macro (scatch(in,"header",true), "geo file format version 4: \"header\" keyword not found");
string label;
omega.size_by_dimension.fill (0);
while (in.good()) {
size_t value;
in >> label;
if (label == "end") break;
if (label == "dimension") { in >> omega.dimension; }
else if (label == "order") { in >> omega.order; }
else if (label == "coordinate_system") { in >> omega.sys_coord; }
else {
size_t variant = 0;
for (; variant < reference_element__max_variant; variant++) {
if (label == label_variant[variant]) break;
}
check_macro (variant < reference_element__max_variant, "unexpected header member: \""<<label<<"\"");
in >> omega.size_by_variant [variant];
size_t dim = reference_element_dimension_by_variant(variant);
omega.size_by_dimension [dim] += omega.size_by_variant [variant];
}
}
check_macro (omega.order == 1, "unsupported geo order = "<<omega.order<<" > 1");
for (omega.map_dimension = 3; omega.map_dimension + 1 >= 1; omega.map_dimension--) {
if (omega.size_by_dimension [omega.map_dimension] != 0) break;
}
in >> label;
check_macro (label == "header", "geo file format version 4: \"end header\" keyword not found");
// get nodes
omega.node.resize (omega.size_by_variant[reference_element__p]);
for (size_t inod = 0; inod < omega.node.size(); ++inod) {
for (size_t i = 0; i < omega.dimension; ++i) {
in >> omega.node[inod][i];
}
}
// get elements
omega.resize (omega.size_by_dimension[omega.map_dimension]);
for (size_t ie = 0; ie < omega.size(); ++ie) {
char name;
in >> name;
omega[ie].setname(name);
for (size_t j = 0; j < omega[ie].size(); ++j) {
in >> omega[ie][j];
}
}
}
int main() {
//-------------------------
// input field
//-------------------------
scatch (cin, "field", true);
size_t version;
cin >> version;
check_macro (version == 1, "version " << version << " field format not yet supported");
size_t nbpts;
cin >> nbpts;
string geo_name;
string approx;
cin >> geo_name
>> approx;
check_macro (approx == "P1", "approx " << approx << " field not yet supported");
vector<double> uh (nbpts);
for (size_t i=0; i<nbpts; i++) {
cin >> uh[i];
}
//----------------------------
// input geo
//----------------------------
my_geo omega;
if (file_exists (geo_name + ".geo")) {
ifstream in ((geo_name + ".geo").c_str());
get_geo (in, omega);
} else if (file_exists (geo_name + ".geo.gz")) {
// use fifo:
// https://stackoverflow.com/questions/40740914/using-istream-to-read-from-named-pipe
bool do_verbose = true;
string fifo_name = tmpnam(0);
#ifdef TO_CLEAN
string command = "rm -f " + fifo_name;
if (do_verbose) cerr << "! " << command << endl;
check_macro (system (command.c_str()) == 0, "adapt: unix command failed");
#endif // TO_CLEAN
unlink (fifo_name.c_str()); // remove
mkfifo (fifo_name.c_str(), S_IWUSR | S_IRUSR); // | S_IRGRP | S_IROTH);
string command = "zcat " + geo_name + ".geo.gz > " + fifo_name + " &";
if (do_verbose) cerr << "! " << command << endl;
check_macro (system (command.c_str()) == 0, "adapt: unix command failed");
ifstream in (fifo_name.c_str());
get_geo (in, omega);
unlink (fifo_name.c_str()); // remove
#ifdef TO_CLEAN
command = "rm -f " + fifo_name;
if (do_verbose) cerr << "! " << command << endl;
check_macro (system (command.c_str()) == 0, "adapt: unix command failed");
#endif // TO_CLEAN
} else {
error_macro ("file \""<<geo_name<< ".geo[.gz]\" not found");
}
//----------------------------
// output pos
//----------------------------
/*
type #list-of-coords #list-of-values
--------------------------------------------------------------------
Scalar point SP 3 1 * nb-time-steps
Vector point VP 3 3 * nb-time-steps
Tensor point TP 3 9 * nb-time-steps
Scalar line SL 6 2 * nb-time-steps
Vector line VL 6 6 * nb-time-steps
Tensor line TL 6 18 * nb-time-steps
Scalar triangle ST 9 3 * nb-time-steps
Vector triangle VT 9 9 * nb-time-steps
Tensor triangle TT 9 27 * nb-time-steps
Scalar quadrangle SQ 12 4 * nb-time-steps
Vector quadrangle VQ 12 12 * nb-time-steps
Tensor quadrangle TQ 12 36 * nb-time-steps
Scalar tetrahedron SS 12 4 * nb-time-steps
Vector tetrahedron VS 12 12 * nb-time-steps
Tensor tetrahedron TS 12 36 * nb-time-steps
Scalar hexahedron SH 24 8 * nb-time-steps
Vector hexahedron VH 24 24 * nb-time-steps
Tensor hexahedron TH 24 72 * nb-time-steps
Scalar prism SI 18 6 * nb-time-steps
Vector prism VI 18 18 * nb-time-steps
Tensor prism TI 18 54 * nb-time-steps
Scalar pyramid SY 15 5 * nb-time-steps
Vector pyramid VY 15 15 * nb-time-steps
Tensor pyramid TY 15 45 * nb-time-steps
*/
cout << setprecision(numeric_limits<double>::digits10)
<< "View \"scalar\" {" << endl;
static char gmsh_pos_type [reference_element__max_variant];
gmsh_pos_type [reference_element__p] = 'P';
gmsh_pos_type [reference_element__e] = 'L';
gmsh_pos_type [reference_element__t] = 'T';
gmsh_pos_type [reference_element__q] = 'Q';
gmsh_pos_type [reference_element__T] = 'S';
gmsh_pos_type [reference_element__H] = 'H';
gmsh_pos_type [reference_element__P] = 'I';
#ifdef TODDO
switch (uh.get_space().valued_tag()) {
// ---------------------------
case space_constant::scalar: {
// ---------------------------
#endif // TODO
for (size_t ie = 0, ne = omega.size(); ie < ne; ie++) {
const geo_element& K = omega[ie];
cout << "S" << gmsh_pos_type[K.variant()] << "(";
for (size_t iloc = 0, nloc = K.size(); iloc < nloc; iloc++) {
const point& x = omega.node[K[iloc]];
for (size_t i_comp = 0; i_comp < 3; i_comp++) {
cout << x[i_comp];
if (i_comp != 2) cout << ",";
}
if (iloc != nloc-1) cout << ",";
}
cout << "){";
for (size_t iloc = 0, nloc = K.size(); iloc < nloc; iloc++) {
cout << uh[K[iloc]];
if (iloc != nloc-1) cout << ",";
}
cout << "};" << endl;
}
#ifdef TODO
}
// ---------------------------
case space_constant::vector: {
// ---------------------------
break;
}
// ---------------------------
case space_constant::tensor: {
// ---------------------------
field_component_const<T,sequential> uh_comp [3][3];
size_t d = uh.get_geo().dimension();
for (size_t i_comp = 0; i_comp < d; i_comp++) {
for (size_t j_comp = 0; j_comp < d; j_comp++) {
uh_comp[i_comp][j_comp].proxy_assign(uh(i_comp,j_comp));
}}
tensor_basic<T> value;
#define HAVE_ID_DEFAULT_TENSOR
#ifdef HAVE_ID_DEFAULT_TENSOR
// default (3,3) is 1 instead of 0
value (0,0) = value (1,1) = value (2,2) = 1;
#endif
for (size_t ie = 0, ne = omega.size(); ie < ne; ie++) {
const geo_element& K = omega[ie];
gmsh << "T" << gmsh_pos_type[K.variant()] << "(";
for (size_t iloc = 0, nloc = K.size(); iloc < nloc; iloc++) {
const point_basic<T>& x = omega.node(K[iloc]);
for (size_t i_comp = 0; i_comp < 3; i_comp++) {
gmsh << x[i_comp];
if (i_comp != 2) gmsh << ",";
}
if (iloc != nloc-1) gmsh << ",";
}
gmsh << "){";
for (size_t iloc = 0, nloc = K.size(); iloc < nloc; iloc++) {
size_t inod = K[iloc];
const point_basic<T>& x = omega.node(inod);
for (size_t i_comp = 0; i_comp < d; i_comp++) {
for (size_t j_comp = 0; j_comp < d; j_comp++) {
value(i_comp,j_comp) = uh_comp [i_comp][j_comp].dof (inod);
}}
for (size_t i_comp = 0; i_comp < 3; i_comp++) {
for (size_t j_comp = 0; j_comp < 3; j_comp++) {
gmsh << value(i_comp,j_comp);
if (i_comp != 2 || j_comp != 2 || iloc != nloc-1) gmsh << ",";
}}
}
gmsh << "};" << endl;
}
break;
}
default: error_macro ("put_gmsh: do not known how to print " << uh.valued() << "-valued field");
}
#endif // TODO
cout << "};" << endl;
}
point
see the point page for the full documentation
rheolef::field_component_const
Definition: field_component.h:143
rheolef::field_component_const::dof
const T & dof(size_type idof) const
Definition: field_component.h:169
rheolef::field_component_const::proxy_assign
field_component_const< T, M > & proxy_assign(const field_component_const< T, M > &uh_comp)
Definition: field_component.h:280
rheolef::geo_element
see the geo_element page for the full documentation
Definition: geo_element.h:102
rheolef::geo_element::variant
size_t variant() const
Definition: field2gmsh_pos.cc:91
rheolef::geo_element::size
size_t size() const
Definition: field2gmsh_pos.cc:94
rheolef::geo_element::size
size_type size() const
Definition: geo_element.h:168
rheolef::geo_element::_variant
size_t _variant
Definition: field2gmsh_pos.cc:97
rheolef::geo_element::dimension
size_t dimension() const
Definition: field2gmsh_pos.cc:92
rheolef::geo_element::setname
void setname(char name)
Definition: field2gmsh_pos.cc:89
rheolef::geo_element::variant
variant_type variant() const
Definition: geo_element.h:161
rheolef::geo_element::geo_element
geo_element()
Definition: field2gmsh_pos.cc:95
rheolef::geo_element::name
char name() const
Definition: field2gmsh_pos.cc:90
rheolef::geo_element::n_vertex
size_t n_vertex() const
Definition: field2gmsh_pos.cc:93
rheolef::point_basic
Definition: point.h:87
rheolef::point_basic::point_basic
point_basic(T x0=T(), T x1=T(), T x2=T())
Definition: field2gmsh_pos.cc:61
rheolef::point_basic::base
std::array< T, 3 > base
Definition: field2gmsh_pos.cc:60
rheolef::std
Definition: vec_expr_v2.h:402
rheolef::tensor_basic
Definition: tensor.h:90
error_macro
#define error_macro(message)
Definition: dis_macros.h:49
dimension
const size_t dimension
Definition: edge.icc:64
n_vertex
const size_t n_vertex
Definition: edge.icc:66
main
int main()
Definition: field2gmsh_pos.cc:180
get_geo
void get_geo(istream &in, my_geo &omega)
Definition: field2gmsh_pos.cc:126
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)")
gmsh
verbose clean transpose logscale grid shrink ball stereo iso volume skipvtk deformation fastfieldload lattice reader_on_stdin color format format format format format format gmsh
Definition: iorheo-members.h:83
mkgeo_ball.variant
int variant
Definition: mkgeo_ball.sh:149
mkgeo_ball.command
string command
Definition: mkgeo_ball.sh:136
mkgeo_ball.d
int d
Definition: mkgeo_ball.sh:154
mkgeo_ball.order
order
Definition: mkgeo_ball.sh:343
mkgeo_ball.dim
int dim
Definition: mkgeo_ball.sh:307
mkgeo_contraction.name
string name
Definition: mkgeo_contraction.sh:133
mkgeo_grid.sys_coord
string sys_coord
Definition: mkgeo_grid.sh:171
rheolef-config.version
version
Definition: rheolef-config.in:126
rheolef::space_constant::tensor
@ tensor
Definition: space_constant.h:138
rheolef::space_constant::scalar
@ scalar
Definition: space_constant.h:136
rheolef::space_constant::vector
@ vector
Definition: space_constant.h:137
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::file_exists
bool file_exists(const std::string &filename)
file_exists: see the rheostream page for the full documentation
Definition: scatch.icc:42
scatch.icc