field_seq_put_gmsh_pos.cc

Go to the documentation of this file.

//
// gmsh pos output, for mesh adaptation
//
// author: Pierre.Saramito@imag.fr
//
// date: 12 may 1997  update: 23 oct 2011
//
#include "rheolef/field.h"
#include "rheolef/piola_util.h"
#include "rheolef/rheostream.h"
#include "rheolef/iorheo.h"
#include "rheolef/field_evaluate.h"
#include "rheolef/space_component.h"
#include "rheolef/field_expr.h"
 
namespace rheolef { 
using namespace std;
 
template <class T>
odiststream&
field_put_gmsh_pos (odiststream& ods, const field_basic<T,sequential>& uh, std::string name)
{
  typedef typename field_basic<T,sequential>::size_type size_type;
  if (name == "") { name = uh.get_space().valued(); }
  ostream& gmsh = ods.os();
  gmsh << setprecision(numeric_limits<T>::digits10);
  check_macro (uh.get_space().get_basis().name() == "P1",
    "gmsh: unsupported " << uh.get_space().get_basis().name() << " approximation");
  check_macro (uh.get_geo().order() == 1,
    "gmsh: unsupported geo order > 1");
  /*
          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
  */
  gmsh << "View \"" << name << "\" {" << 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';
  const geo_basic<T,sequential>& omega = uh.get_geo();
  switch (uh.get_space().valued_tag()) {
    // ---------------------------
    case space_constant::scalar: {
    // ---------------------------
      for (size_type ie = 0, ne = omega.size(); ie < ne; ie++) {
        const geo_element& K = omega[ie];
        gmsh << "S" << gmsh_pos_type[K.variant()] << "(";
        for (size_type iloc = 0, nloc = K.size(); iloc < nloc; iloc++) {
          const point_basic<T>& x = omega.node(K[iloc]);
          for (size_type 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_type iloc = 0, nloc = K.size(); iloc < nloc; iloc++) {
          gmsh << uh.dof(K[iloc]);
          if (iloc != nloc-1) gmsh << ",";
        }
        gmsh << "};" << endl;
      }
      break;
    }
#ifdef TODO
    // ---------------------------
    case space_constant::vector: {
    // ---------------------------
      break;
    }
#endif // TODO
    // ---------------------------
    case space_constant::tensor: {
    // ---------------------------
      size_type d = uh.get_geo().dimension();
#ifdef TO_CLEAN
      details::field_rdof_sliced_const<field_basic<T,sequential>> uh_comp [3][3];
      for (size_type i_comp = 0; i_comp < d; i_comp++) {
      for (size_type j_comp = 0; j_comp < d; j_comp++) {
        uh_comp[i_comp][j_comp].proxy_assign(uh(i_comp,j_comp));
      }}
#endif // TO_CLEAN
      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_type ie = 0, ne = omega.size(); ie < ne; ie++) {
        const geo_element& K = omega[ie];
        gmsh << "T" << gmsh_pos_type[K.variant()] << "(";
        for (size_type iloc = 0, nloc = K.size(); iloc < nloc; iloc++) {
          const point_basic<T>& x = omega.node(K[iloc]);
          for (size_type 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_type iloc = 0, nloc = K.size(); iloc < nloc; iloc++) {
          size_type inod = K[iloc];
          const point_basic<T>& x = omega.node(inod);
          for (size_type i_comp = 0; i_comp < d; i_comp++) {
          for (size_type j_comp = 0; j_comp < d; j_comp++) {
            value(i_comp,j_comp) = uh (i_comp,j_comp).dof (inod);
          }}
          for (size_type i_comp = 0; i_comp < 3; i_comp++) {
          for (size_type 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");
  }
  gmsh << "};" << endl;
  return ods;
}
template <class T>
odiststream&
field_put_gmsh_pos (odiststream& ods, const field_basic<T,sequential>& uh)
{
  return field_put_gmsh_pos (ods, uh, "");
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
template odiststream& field_put_gmsh_pos<Float>  (odiststream&, const field_basic<Float,sequential>&, std::string);
template odiststream& field_put_gmsh_pos<Float>  (odiststream&, const field_basic<Float,sequential>&);
 
}// namespace rheolef