geo_seq_put_gmsh.cc
Go to the documentation of this file.
//
// gmsh output
//
// author: Pierre.Saramito@imag.fr
//
// date: 26 mars 2012
//
#include "rheolef/geo.h"
#include "rheolef/space_numbering.h"
#include "rheolef/piola_util.h"
#include "rheolef/rheostream.h"
#include "rheolef/iorheo.h"
using namespace std;
namespace rheolef {
// ----------------------------------------------------------------------------
// one element puts
// ----------------------------------------------------------------------------
template <class T>
static
void
put_edge (ostream& gmsh, const geo_element& K, const basis_basic<T>& my_numb, const geo_basic<T,sequential>& omega)
{
typedef typename geo_basic<T,sequential>::size_type size_type;
static size_type order2gmsh_type [11] = {0, 1, 8, 26, 27, 28, 62, 63, 64, 65, 66 };
size_type my_order = my_numb.degree();
check_macro (my_order <= 10, "gmsh output: element 'e' order > 10 not yet supported");
std::vector<size_type> inod;
space_numbering::dis_idof (my_numb, omega.sizes(), K, inod);
gmsh << K.dis_ie()+1 << " " << order2gmsh_type [my_order] << " 2 99 2"; // TODO: domains
for (size_type iloc = 0, nloc = inod.size(); iloc < nloc; iloc++) {
gmsh << " " << inod[iloc]+1;
}
gmsh << endl;
}
template <class T>
static
void
put_triangle (ostream& gmsh, const geo_element& K, const basis_basic<T>& my_numb, const geo_basic<T,sequential>& omega)
{
typedef typename geo_basic<T,sequential>::size_type size_type;
static size_type order2gmsh_type [11] = {0, 2, 9, 21, 23, 25, 42, 43, 44, 45, 46};
size_type my_order = my_numb.degree();
// TODO: permutations of internal nodes for order >= 4
check_macro (my_order <= 3, "gmsh output: element 't' order > 10 not yet supported");
std::vector<size_type> inod;
space_numbering::dis_idof (my_numb, omega.sizes(), K, inod);
gmsh << K.dis_ie()+1 << " " << order2gmsh_type [my_order] << " 2 99 2"; // TODO: domains
for (size_type iloc = 0, nloc = inod.size(); iloc < nloc; iloc++) {
gmsh << " " << inod[iloc]+1;
}
gmsh << endl;
}
#ifdef TODO
template <class T>
static
void
put_quadrangle (ostream& gmsh, const geo_element& K, const basis_basic<T>& my_numb, const geo_basic<T,sequential>& omega)
{
typedef typename geo_basic<T,sequential>::size_type size_type;
typedef point_basic<size_type> ilat;
std::vector<size_type> inod;
space_numbering::dis_idof (my_numb, omega.sizes(), K, inod);
size_type my_order = my_numb.degree();
for (size_type i = 0; i < my_order; i++) {
for (size_type j = 0; j < my_order; j++) {
size_type loc_inod00 = reference_element_q::ilat2loc_inod (my_order, ilat(i, j));
size_type loc_inod10 = reference_element_q::ilat2loc_inod (my_order, ilat(i+1, j));
size_type loc_inod11 = reference_element_q::ilat2loc_inod (my_order, ilat(i+1, j+1));
size_type loc_inod01 = reference_element_q::ilat2loc_inod (my_order, ilat(i, j+1));
gmsh << "4\t" << inod[loc_inod00] << " "
<< inod[loc_inod10] << " "
<< inod[loc_inod11] << " "
<< inod[loc_inod01] << endl;
}
}
}
template <class T>
static
void
put_tetrahedron (ostream& gmsh, const geo_element& K, const basis_basic<T>& my_numb, const geo_basic<T,sequential>& omega)
{
typedef typename geo_basic<T,sequential>::size_type size_type;
typedef point_basic<size_type> ilat;
std::vector<size_type> inod;
space_numbering::dis_idof (my_numb, omega.sizes(), K, inod);
size_type my_order = my_numb.degree();
for (size_type i = 0; i < my_order; i++) {
for (size_type j = 0; i+j < my_order; j++) {
for (size_type k = 0; i+j+k < my_order; k++) {
size_type loc_inod000 = reference_element_T::ilat2loc_inod (my_order, ilat(i, j, k));
size_type loc_inod100 = reference_element_T::ilat2loc_inod (my_order, ilat(i+1, j, k));
size_type loc_inod010 = reference_element_T::ilat2loc_inod (my_order, ilat(i, j+1, k));
size_type loc_inod001 = reference_element_T::ilat2loc_inod (my_order, ilat(i, j, k+1));
gmsh << "4\t" << inod[loc_inod000] << " "
<< inod[loc_inod100] << " "
<< inod[loc_inod010] << " "
<< inod[loc_inod001] << endl;
if (i+j+k+2 > my_order) continue;
// complete the ijk-th cube: 4 more tetras
size_type loc_inod110 = reference_element_T::ilat2loc_inod (my_order, ilat(i+1, j+1, k));
size_type loc_inod101 = reference_element_T::ilat2loc_inod (my_order, ilat(i+1, j, k+1));
size_type loc_inod011 = reference_element_T::ilat2loc_inod (my_order, ilat(i, j+1, k+1));
gmsh << "4\t" << inod[loc_inod100] << " " // face in x0 & x2 direction
<< inod[loc_inod101] << " "
<< inod[loc_inod010] << " "
<< inod[loc_inod001] << endl
<< "4\t" << inod[loc_inod010] << " " // face in x1 & x2 direction
<< inod[loc_inod011] << " "
<< inod[loc_inod001] << " "
<< inod[loc_inod101] << endl
<< "4\t" << inod[loc_inod100] << " "
<< inod[loc_inod101] << " "
<< inod[loc_inod110] << " "
<< inod[loc_inod010] << endl
<< "4\t" << inod[loc_inod010] << " "
<< inod[loc_inod110] << " "
<< inod[loc_inod011] << " "
<< inod[loc_inod101] << endl;
// the last 6th sub-tetra that fully fills the ijk-th cube
if (i+j+k+3 > my_order) continue;
size_type loc_inod111 = reference_element_T::ilat2loc_inod (my_order, ilat(i+1, j+1, k+1));
gmsh << "4\t" << inod[loc_inod111] << " " // face in x0 & x2 direction
<< inod[loc_inod101] << " "
<< inod[loc_inod011] << " "
<< inod[loc_inod110] << endl;
}
}
}
}
static
void
raw_put_prism (ostream& gmsh,
size_t i000, size_t i100, size_t i010,
size_t i001, size_t i101, size_t i011)
{
// gmsh prism has swaped x & y axis order: 00z 10z 01z replaced by 00z 01z 10z
gmsh << "6\t" << i000 << " "
<< i010 << " "
<< i100 << " "
<< i001 << " "
<< i011 << " "
<< i101 << endl;
}
template <class T>
static
void
put_prism (ostream& gmsh, const geo_element& K, const basis_basic<T>& my_numb, const geo_basic<T,sequential>& omega, const disarray<point_basic<Float>,sequential>& my_node)
{
typedef typename geo_basic<T,sequential>::size_type size_type;
typedef point_basic<size_type> ilat;
std::vector<size_type> inod;
space_numbering::dis_idof (my_numb, omega.sizes(), K, inod);
size_type my_order = my_numb.degree();
for (size_type k = 0; k < my_order; k++) {
for (size_type j = 0; j < my_order; j++) {
for (size_type i = 0; i+j < my_order; i++) {
size_type loc_inod000 = reference_element_P::ilat2loc_inod (my_order, ilat(i, j, k));
size_type loc_inod100 = reference_element_P::ilat2loc_inod (my_order, ilat(i+1, j, k));
size_type loc_inod010 = reference_element_P::ilat2loc_inod (my_order, ilat(i, j+1, k));
size_type loc_inod001 = reference_element_P::ilat2loc_inod (my_order, ilat(i, j, k+1));
size_type loc_inod101 = reference_element_P::ilat2loc_inod (my_order, ilat(i+1, j, k+1));
size_type loc_inod011 = reference_element_P::ilat2loc_inod (my_order, ilat(i, j+1, k+1));
raw_put_prism (gmsh,
inod[loc_inod000],
inod[loc_inod100],
inod[loc_inod010],
inod[loc_inod001],
inod[loc_inod101],
inod[loc_inod011]);
if (i+j+1 >= my_order) continue;
size_type loc_inod110 = reference_element_P::ilat2loc_inod (my_order, ilat(i+1, j+1, k));
size_type loc_inod111 = reference_element_P::ilat2loc_inod (my_order, ilat(i+1, j+1, k+1));
raw_put_prism (gmsh,
inod[loc_inod100],
inod[loc_inod110],
inod[loc_inod010],
inod[loc_inod101],
inod[loc_inod111],
inod[loc_inod011]);
}
}
}
}
template <class T>
static
void
put_hexahedron (ostream& gmsh, const geo_element& K, const basis_basic<T>& my_numb, const geo_basic<T,sequential>& omega)
{
typedef typename geo_basic<T,sequential>::size_type size_type;
typedef point_basic<size_type> ilat;
std::vector<size_type> inod;
space_numbering::dis_idof (my_numb, omega.sizes(), K, inod);
size_type my_order = my_numb.degree();
for (size_type i = 0; i < my_order; i++) {
for (size_type j = 0; j < my_order; j++) {
for (size_type k = 0; k < my_order; k++) {
size_type loc_inod000 = reference_element_H::ilat2loc_inod (my_order, ilat(i, j, k));
size_type loc_inod100 = reference_element_H::ilat2loc_inod (my_order, ilat(i+1, j, k));
size_type loc_inod110 = reference_element_H::ilat2loc_inod (my_order, ilat(i+1, j+1, k));
size_type loc_inod010 = reference_element_H::ilat2loc_inod (my_order, ilat(i, j+1, k));
size_type loc_inod001 = reference_element_H::ilat2loc_inod (my_order, ilat(i, j, k+1));
size_type loc_inod101 = reference_element_H::ilat2loc_inod (my_order, ilat(i+1, j, k+1));
size_type loc_inod011 = reference_element_H::ilat2loc_inod (my_order, ilat(i, j+1, k+1));
size_type loc_inod111 = reference_element_H::ilat2loc_inod (my_order, ilat(i+1, j+1, k+1));
gmsh << "8\t" << inod[loc_inod000] << " "
<< inod[loc_inod100] << " "
<< inod[loc_inod110] << " "
<< inod[loc_inod010] << " "
<< inod[loc_inod001] << " "
<< inod[loc_inod101] << " "
<< inod[loc_inod111] << " "
<< inod[loc_inod011] << endl;
}
}
}
}
#endif // TODO
template <class T>
static
void
put (ostream& gmsh, const geo_element& K, const basis_basic<T>& my_numb, const geo_basic<T,sequential>& omega, const disarray<point_basic<Float>,sequential>& my_node)
{
switch (K.variant()) {
#ifdef TODO
case reference_element::p: gmsh << "1\t" << K[0] << endl; break;
#endif // TODO
case reference_element::e: put_edge (gmsh, K, my_numb, omega); break;
case reference_element::t: put_triangle (gmsh, K, my_numb, omega); break;
#ifdef TODO
case reference_element::q: put_quadrangle (gmsh, K, my_numb, omega); break;
case reference_element::T: put_tetrahedron (gmsh, K, my_numb, omega); break;
case reference_element::P: put_prism (gmsh, K, my_numb, omega, my_node); break;
case reference_element::H: put_hexahedron (gmsh, K, my_numb, omega); break;
#endif // TODO
default: error_macro ("unsupported element variant `" << K.name() <<"'");
}
}
// ----------------------------------------------------------------------------
// geo puts
// ----------------------------------------------------------------------------
template <class T>
odiststream&
geo_put_gmsh (odiststream& ops, const geo_basic<T,sequential>& omega, const basis_basic<T>& my_numb, const disarray<point_basic<T>,sequential>& my_node)
{
//
// 0) pre-requises
//
typedef typename geo_basic<T,sequential>::size_type size_type;
size_type my_order = my_numb.degree();
ostream& gmsh = ops.os();
//
// 1) put header
//
gmsh << setprecision (std::numeric_limits<T>::digits10)
<< "$MeshFormat" << endl
<< "2.2 0 8" << endl
<< "$EndMeshFormat" << endl;
// TODO: add domains: scan by domain and add for earch element to a domain list
//
// 2) put nodes
//
gmsh << "$Nodes" << endl
<< my_node.size() << endl;
for (size_type inod = 0, nnod = my_node.size(); inod < nnod; inod++) {
gmsh << inod+1 << " " << my_node[inod] << endl;
}
gmsh << "$EndNodes" << endl;
//
// 3) put elements
//
size_type map_dim = omega.map_dimension();
gmsh << "$Elements" << endl
<< omega.size() << endl;
for (size_type ie = 0, ne = omega.size(); ie < ne; ie++) {
const geo_element& K = omega.get_geo_element (map_dim, ie);
put (gmsh, K, my_numb, omega, my_node);
}
gmsh << "$EndElements" << endl;
return ops;
}
template <class T>
odiststream&
geo_put_gmsh (odiststream& ops, const geo_basic<T,sequential>& omega)
{
basis_basic<T> my_numb ("P" + itos(omega.order()));
return geo_put_gmsh (ops, omega, my_numb, omega.get_nodes());
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
template odiststream& geo_put_gmsh<Float> (odiststream&, const geo_basic<Float,sequential>&, const basis_basic<Float>&, const disarray<point_basic<Float>,sequential>&);
template odiststream& geo_put_gmsh<Float> (odiststream&, const geo_basic<Float,sequential>&);
}// namespace rheolef
size_type
field::size_type size_type
Definition: branch.cc:425
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.order
order
Definition: mkgeo_ball.sh:343
rheolef
This file is part of Rheolef.
Definition: compiler_eigen.h:37