geometry.cpp

Go to the documentation of this file.

#include "geometry.h"
#include "point.h"
#include "cell.h"
#include "module_t3d.h"
#include "module_oofem.h"
#include "dupl.h"
#ifdef __MEER_MODULE
#include "module_meer.h"
#endif

#include "mathlib.h"
#include "gelib.h"
#include "librw.h"

#include <time.h>


namespace midaspace {

//*  ********************************************************************************************
//*  *** *** *** ***                      CLASS MESHGEOMETRY                      *** *** *** ***
//*  ********************************************************************************************
Geometry :: Geometry (long gid, const Problem *pd, MMprocessing mmp) : ProblemSubject (gid, pd)
{
  //*  ***   GENERAL INFO   ***
  gp = mmp;
  //* STATE
  nodeprop = elemprop = false;
  
  //*  ***   GEOMETRY DATA   ***
  origPjntsNum = origElemsNum = 0;
  Pjnts.initialize(0,10);  //* !!!!!!! toto zrusit
  Edges.initialize(0,10);
  Faces.initialize(0,10);
  Elems.initialize(0,10);
  
  //*  ***   GEOMETRY STATS   ***
  connectivity_assembled = false;
  VTKCTcl = VTKCTct = VTKCTch = -1;
  aver_elem_circ = -1.0;
  anyBeamElem = -1;
  boubox_assembled = false;
}
Geometry :: ~Geometry ()
{
  Pjnts.delete_objects();
  Edges.delete_objects();
  Faces.delete_objects();
  Elems.delete_objects();
}

void Geometry :: initialize (void)
{
  origPjntsNum = Pjnts();
  origElemsNum = Elems();
  
  long i;
  
  //* set property 1 for all elements
  if (Pd->give_PDBO(PDBO_allElemsProp1))
    for (i=0; i<Elems(); i++)  Elems[i]->set_mprop(1);
  
  //*  *** initialization of points ***
  for (i=0; i<Pjnts(); i++)  Pjnts[i]->initialize();
  //*  *** initialization of elements ***
  for (i=0; i<Elems(); i++)  Elems[i]->initialize();
}
void Geometry :: finitialize (void)
{
  long i;
  
  //*  *** initialization of points ***
  for (i=0; i<Pjnts(); i++)  Pjnts[i]->finitialize();
  //*  *** initialization of elements ***
  for (i=0; i<Elems(); i++)  Elems[i]->finitialize();
}
void Geometry :: checkConsistency (void) const
{
  long i;
  
  // *** Nodes/Elements ***
  for (i=0; i<Pjnts(); i++) Pjnts[i]->checkConsistency();
  for (i=0; i<Elems(); i++) Elems[i]->checkConsistency();
  
  // *** Node/Edge/Face/Elem[id]->gid = id ***
  for (i=0; i<Pjnts(); i++)  if ( !Pjnts[i]  &&  Pjnts[i]->give_ID() != i )  _errorr3 ("Node %ld is on position %ld, it is forbidden", Pjnts[i]->give_ID(), i);
  for (i=0; i<Edges(); i++)  if ( !Edges[i]  ||  Edges[i]->give_ID() != i )  _errorr3 ("Edge %ld is on position %ld, it is forbidden", Edges[i]->give_ID(), i);
  for (i=0; i<Faces(); i++)  if ( !Faces[i]  &&  Faces[i]->give_ID() != i )  _errorr3 ("Face %ld is on position %ld, it is forbidden", Faces[i]->give_ID(), i);
  for (i=0; i<Elems(); i++)  if ( !Elems[i]  &&  Elems[i]->give_ID() != i )  _errorr3 ("Elem %ld is on position %ld, it is forbidden", Elems[i]->give_ID(), i);
  
# ifdef DEBUG
  long j;
  // *** check right node numbers of elements ***
  for (i=0; i<Elems(); i++)
    for (j=0; j<Elems[i]->give_nno(); j++)
      if (Elems[i]->give_point(j) == NULL)
    _errorr3 ( "The %d-th node on element %d is not set", j, i);
# endif
}



//*  *************************
//*  ***   GEOMETRY DATA   ***
//*  *************************
// ###
// ###   manipulation with arrays of components (Nodes, Vertices, Elements....
// ###

// vem kod z  void LocalCoordSystem :: initialize (Geometry *geom)
const Point* Geometry :: give_Pjnt_with_prop (long prp, bool unique) const
{
  Point *p=NULL;
  for (long i=0; i<Pjnts(); i++)
    if (Pjnts[i]->has_mproperty (prp)) {
      if (unique && p)  _errorr2("node with property %d is not unique", prp);
      p = Pjnts[i];
      if (!unique)  break;
    }
  return p;
}



//*  \todo             dodelat
//*  !!!!!!!!!!!!!!!!  tady musi byt zohledneno ze se menei STATS
//*  !!!!!!!!!!!!!!!!  kdyz pridam element tak musim resetnout statistiky a opravit connectivity


template <class T>
void TF_GPA_delete_damned (T &list)
{
  for (long i=0; i<list(); i++)
    if (list[i]->give_delete_flag())
      list.delete_object(i);
  
  TF_GPA_shake_down_reidoid(list);
}

template <class T>
void TF_GPA_add_another (GPA<T> &list, T *comp, const Geometry *mg)           { comp->reset_Geom(mg);  comp->set_ID(list());  list.add(comp); }

template <class T>
void TF_GPA_wdg_another (GPA<T> &list, T *comp, const Geometry *mg, long nid) { comp->set_Geom(mg);  comp->set_ID(nid);     list.wedge(nid, comp); TF_GPA_reidoid (list); }

void Geometry :: add_another_Pjnt(Point   *comp)           { TF_GPA_add_another (Pjnts, comp, this);       this->geom_stats_assembling(comp); }
void Geometry :: add_another_Edge(Edge    *comp)           { TF_GPA_add_another (Edges, comp, this);       this->geom_stats_assembling(comp); }
void Geometry :: add_another_Face(Face    *comp)           { TF_GPA_add_another (Faces, comp, this);       this->geom_stats_assembling(comp); }
void Geometry :: add_another_Elem(Element *comp)           { TF_GPA_add_another (Elems, comp, this);       this->geom_stats_assembling(comp); }

void Geometry :: wdg_another_Elem(Element *comp, long nid) { TF_GPA_wdg_another (Elems, comp, this, nid);  this->geom_stats_assembling(comp); }

void Geometry :: replace_Pjnt_by (Point *Pold, Point *Pnew)
{
  Pold->make_invisible(Pnew, false);
  Pjnts.assign(Pold->give_ID(), Pnew);
  Pnew->set_ID(Pold->give_ID());
  delete Pold;
}


//
//  NEni dodelana, aktivuj damn_orphaned
//
void Geometry :: kill_damned_elems (void)
{
  long i;
  
  for (i=0; i<Elems(); i++)
    if ( Elems[i]->give_delete_flag() )
      Elems[i]->connectivity_removing ();
  
  //for (i=0; i<Nn ; i++)  Nodes[i]->damn_if_orphan();
  //for (i=0; i<Ned; i++)  Edges[i]->damn_if_orphan();
  //for (i=0; i<Nfa; i++)  Faces[i]->damn_if_orphan();
  //
  // tady killnout osirele nodes a edges a elems...
  //
  
  TF_GPA_delete_damned (Elems);
}

void Geometry :: letoff_aa (void)
{
  Pjnts.resize(0);
  Edges.resize(0);
  Faces.resize(0);
  Elems.resize(0);
}


//*  **************************
//*  ***   GEOMETRY STATS   ***
//*  **************************
void Geometry :: connectivity_allocation (void)
{
  if (2*Pjnts() > Edges.give_asize())  Edges.reallocup(2*Pjnts());
  if (2*Pjnts() > Faces.give_asize())  Faces.reallocup(2*Pjnts());
}

void Geometry :: connectivity_assembling (void)
{
  if (connectivity_assembled) return;
  connectivity_assembled = true;
  
  this->connectivity_allocation();
  
  for (long i=0; i<Elems(); i++)
    Elems[i]->connectivity_assembling(false);
}

void Geometry :: VTKCT_assembling (long indx)
{
  long i;
  
  if (indx > -1) {
    i = indx;
    indx++;
  }
  else {
    VTKCTcl = VTKCTct = VTKCTch = 0;
    i=0;
    indx = Elems();
  }
  
  for (; i<indx; i++)
    switch (CellGeometry_e2i_VTK(Elems[i]->give_cellGeom())) {
    case 3:   VTKCTcl++;     break;
    case 5:   VTKCTct++;     break;
    case 9:   VTKCTct++;     break;
    case 12:  VTKCTch++;     break;
    case 0:   _errorr2 ("give_VTKcelltype is not implemented on element %ld", i);
    default:  _errorr ("unsuported VTK cell type");
    }
}

void Geometry :: give_vtk_polydata_counts (long &cl, long &cp) const
{
  if (! this->VTKCT_is_assembled())  ((Mesh *)this)->VTKCT_assembling ();
  
  cl = VTKCTcl;
  cp = VTKCTct;
}

bool Geometry :: only_vtk_polydata (void) const
{
  if (! this->VTKCT_is_assembled())  ((Mesh *)this)->VTKCT_assembling ();
  
  if (VTKCTch)  return false;
  ;             return true;
}

bool Geometry :: only_1d_elements (void) const
{
  if (! this->VTKCT_is_assembled())  ((Mesh *)this)->VTKCT_assembling ();
  
  if (VTKCTcl > 0  &&  VTKCTct == 0  &&  VTKCTch == 0)  return true;
  ;                                                     return false;
}

void Geometry :: aver_elem_circ_assembling (long indx)
{
  if (indx > -1) {
    aver_elem_circ = ( aver_elem_circ * (Elems() - 1) + Elems[indx]->give_circum() ) / Elems();
  }
  else {
    aver_elem_circ = 0.0;
    for (long i=0; i<Elems(); i++)
      aver_elem_circ += Elems[i]->give_circum();
    
    aver_elem_circ /= Elems();
  }
}
void Geometry :: anyBeamElem_assembling (long indx)
{
  if (indx > -1) {
    if (Elems[indx]->give_classid() == classBeam)  anyBeamElem++;
  }
  else {
    anyBeamElem = 0;
    for (long i=0; i<Elems(); i++)
      if (Elems[i]->give_classid() == classBeam)  anyBeamElem++;
  }
}
double Geometry :: give_aver_elem_circ (void) const { if (aver_elem_circ < 0.0)  ((Mesh *)this)->aver_elem_circ_assembling();  return aver_elem_circ; }
long   Geometry :: give_anyBeamElem    (void) const { if (anyBeamElem    < 0  )  ((Mesh *)this)->anyBeamElem_assembling();     return anyBeamElem;    }

void Geometry :: boubox_assembling (long indx)
{
  const PoinT* coo;
  
  if (indx > -1) {
    coo = Pjnts[indx]->give_coords();
    
    for (long j=0; j<3; j++)
      if      (boubox_min[j] > (*coo)[j])  boubox_min[j] = (*coo)[j];
      else if (boubox_max[j] < (*coo)[j])  boubox_max[j] = (*coo)[j];
  }
  else {
    long i, j;
    const PoinT* coo;
    // set boundary box
    coo = Pjnts[0]->give_coords();
    for (i=0; i<3; i++)
      boubox_min[i] = boubox_max[i] = (*coo)[i];
    
    for (i=1; i<Pjnts(); i++) {
      coo = Pjnts[i]->give_coords();
      for (j=0; j<3; j++)
    if      (boubox_min[j] > (*coo)[j])  boubox_min[j] = (*coo)[j];
    else if (boubox_max[j] < (*coo)[j])  boubox_max[j] = (*coo)[j];
    }
    
    boubox_assembled = true;
  }
  
  boubox_diff.copy(&boubox_max);
  boubox_diff.sub (&boubox_min);
}
const PoinT* Geometry :: give_boubox_diff (void) const
{
  if (boubox_assembled == false)
    ((Mesh *)this)->boubox_assembling();
  
  return &boubox_diff;
}

void Geometry :: geom_stats_assembling (const GeometryComponent *comp)
{
  if (dynamic_cast<const Point *>(comp)) {
    if (this->boubox_assembled)  this->boubox_assembling(comp->give_ID());
  }
  else if (dynamic_cast<const Cell *> (comp)) {
    if (this->connectivity_assembled)  ((Cell*)comp)->connectivity_assembling(true);
    
    if (dynamic_cast<const Element *> (comp)) {
      if (this->VTKCT_is_assembled())  this->VTKCT_assembling          (comp->give_ID());
      if (this->aver_elem_circ > 0.0)  this->aver_elem_circ_assembling (comp->give_ID());
      if (this->anyBeamElem    >= 0 )  this->anyBeamElem_assembling    (comp->give_ID());
    }
  }
  else errol;
}


VTKdataset Geometry :: giveVTKDS (void)
{
  for (long i=0; i<Elems(); i++)
    if (Elems[i]->give_cellGeom() == CG_Hexahedron)  return VTKDS_unstruct;
  
  return VTKDS_polydata;
}

VTKPDtopology Geometry :: give_VTKPDtopology_of_elem (long i) const
{
  return Elems[i]->give_VTKPDtopology();
}

void Geometry :: print_info (void) const
{
  //* model geom info
  this->print_geom_info ();
  
  //* mesh quality
  const Mesh *m = dynamic_cast<const Mesh *>(this);
  if (m)
    m->mesh_quality();
  
}

void Geometry :: print_geom_info (void) const
{
  Pd->sodriver()->print_info_message (SODE_start_blue_arrow, "%s%ld] - geometry stats:", (this->isModel() ? "Model[" : "Mesh["), this->give_ID());
  
  //virtual CellGeometry give_cellGeom (void) const { return CG_Line; }
  long i;
  long cBeam,cTriangle,cQuads,cPolygon,cTetra,cBrick;
  ;    cBeam=cTriangle=cQuads=cPolygon=cTetra=cBrick = 0;
  long Ne = Elems();
  
  for (i=0; i<Ne; i++)
    switch (Elems[i]->give_cellGeom()) {
    case CG_Line:        cBeam++;     break;
    case CG_Triangle:    cTriangle++; break;
    case CG_Quadrangle:  cQuads++;    break;
    case CG_Polygon:     cPolygon++;  break;
    case CG_Tetrahedron: cTetra++;    break;
    case CG_Hexahedron:  cBrick++;    break;
    default: _errorr2 ("Unknown element geom type (element number %d)", i+1);
    }
  
  ;              Pd->sodriver()->print_info_message (SODE_row,"\n     Number of nodes:       %6ld",                 Pjnts());
  if (cBeam)     Pd->sodriver()->print_info_message (SODE_row,  "     Number of beams:       %6ld (%2ld%% of %ld)", cBeam,     (100*cBeam    /Ne), Ne);
  if (cTriangle) Pd->sodriver()->print_info_message (SODE_row,  "     Number of triangles:   %6ld (%2ld%% of %ld)", cTriangle, (100*cTriangle/Ne), Ne);
  if (cQuads)    Pd->sodriver()->print_info_message (SODE_row,  "     Number of quadrangles: %6ld (%2ld%% of %ld)", cQuads,    (100*cQuads   /Ne), Ne);
  if (cPolygon)  Pd->sodriver()->print_info_message (SODE_row,  "     Number of polygons:    %6ld (%2ld%% of %ld)", cPolygon,  (100*cPolygon /Ne), Ne);
  if (cTetra)    Pd->sodriver()->print_info_message (SODE_row,  "     Number of tetras:      %6ld (%2ld%% of %ld)", cTetra,    (100*cTetra   /Ne), Ne);
  if (cBrick)    Pd->sodriver()->print_info_message (SODE_row,  "     Number of bricks:      %6ld (%2ld%% of %ld)", cBrick,    (100*cBrick   /Ne), Ne);
  
  //* volumes
  int dim;
  double lav, cs, LAV[3], VOL[2];
  cs = LAV[0] = LAV[1] = LAV[2] =  VOL[0] = VOL[1] = 0.0;
  //
  bool novol = false;
  for (i=0; i<Ne; i++) {
    dim = Elems[i]->give_dimension();
    lav = Elems[i]->give_lav();
    
    if (Pd->give_fulldata()) {
      //* toto muzu predelat na volani fce element->give_volume()
      if      (dim == 1 &&  (cs = Elems[i]->give_elemAttribs()->give_cs()->give_area())      <= 0.0)  novol = true;
      else if (dim == 2 &&  (cs = Elems[i]->give_elemAttribs()->give_cs()->give_thickness()) <= 0.0)  novol = true;
    }
    
    if      (dim == 1) { LAV[0] += lav;  VOL[0] += cs*lav; }
    else if (dim == 2) { LAV[1] += lav;  VOL[1] += cs*lav; }
    else                 LAV[2] += lav;
  }
  if (novol)
    VOL[0] = VOL[1] = 0.0;
  
  if (LAV[0] > 0.0)  fprintf (stdout, "     Suma of lengths/volumes:    %6g / %6g\n", LAV[0], VOL[0]);
  if (LAV[1] > 0.0)  fprintf (stdout, "     Suma of areas  /volumes:    %6g / %6g\n", LAV[1], VOL[1]);
  if (LAV[2] > 0.0)  fprintf (stdout, "     Suma of volumes/volumes:    %6g / %6g\n", LAV[2], LAV[2]);
}


//*  ************************************
//*  ***   READING OF GEOMETRY MESH   ***
//*  ************************************
FILE* scan_DATA_head (FILE *stream, const char *str, const char *expnumtype, int *noc, const char *descript)
{
  int n;
  char auxs[255], LINE[1023];
  sprintf (auxs, "Error in section %s", descript);
  
  SP_scan_expected_word_exit   (str, expnumtype, auxs, CASE);
  SP_scan_number               (str, n);
  if (*noc && *noc != n)  _errorr3 ("in SCALAR, number of components is not \"%d\" but %d", *noc, n);
  else *noc = n;
  fgets (LINE, 1023, stream);  str=LINE;
  SP_scan_expected_word_exit   (str, "LOOKUP_TABLE", auxs, CASE);
  
  return stream;
}

void scan_DATA_field (Stream *stream, const char *str, int nexc, long m, Lmtrx *dataL, Dmtrx *dataD, Lvctr *sprs, const char *name)
{
  char flint[8];
  if      ( dataL && !dataD) { if (stream->isFile()) sprintf(flint,"int");    else sprintf(flint,"Int32");   }
  else if (!dataL &&  dataD) { if (stream->isFile()) sprintf(flint,"float");  else sprintf(flint,"Float32"); }
  else errol;
  
  Stream strm;
  int n=0;
  if (stream->isFile()) strm.redefine( scan_DATA_head (stream->file(),  str,               flint, &n, name));
  else                  strm.redefine( XP_giveDAtext  (stream->tixel(), 0, false, "ascii", flint, name, &n));
  
  if (nexc && nexc != n)  _errorr4 ("in DataArray %s, number of components is not \"%d\" but %d", name, nexc, n);
  
  if (sprs)  sprs->resize_ignore_vals(m);
  long i,j;
  if (dataL) { dataL->resize_ignore_vals(m, n);  for (i=0; i<m; i++) { if (sprs)  ST_scan_number (&strm, (*sprs)[i]);  for (j=0; j<n; j++) if (ST_scan_number (&strm, (*dataL)(i,j)) != true)  _errorr3 ("There is few \"%s\" numbers in section %s", flint, name); } }
  else       { dataD->resize_ignore_vals(m, n);  for (i=0; i<m; i++) { if (sprs)  ST_scan_number (&strm, (*sprs)[i]);  for (j=0; j<n; j++) if (ST_scan_number (&strm, (*dataD)(i,j)) != true)  _errorr3 ("There is few \"%s\" numbers in section %s", flint, name); } }
  
  if (stream->isFile())  FP_skip_line (stream->file());
}


Point* Model :: allocate_point (long gid, const PoinT *coo, char attribs_alloc)
{
  return (Point*)(new Vertex (this, gid, coo, attribs_alloc));
}

Point* Mesh  :: allocate_point (long gid, const PoinT *coo, char attribs_alloc)
{
  return (Point*)(new Node   (this, gid, coo, attribs_alloc) );
}


Element* allocate_element (bool mdl, CellGeometry eg, long gid, long oid, Geometry *geom, int np=0, bool attribs_alloc=true, long dom=0, long lid=-1)
{
  if (lid == -1) lid = gid;
  
  if (mdl) {
    switch (eg) {
    case CG_Line:        if (np>2)  return new PolyLine   (gid, oid, geom, np);
                         else       return new Line       (gid, oid, geom);  break;
    case CG_Polygon:                return new PolygonMdl (gid, oid, geom);  break;
    default: errol;
    }
  }
  else {
    switch (eg) {
    case CG_Line:        if (np>2)  errol;
                         else       return new Beam       (gid, oid, geom, attribs_alloc, dom, lid);  break;
    case CG_Triangle:               return new Triangle   (gid, oid, geom, attribs_alloc, dom, lid);  break;
    case CG_Quadrangle:             return new Quadrangle (gid, oid, geom, attribs_alloc, dom, lid);  break;
    case CG_Hexahedron:             return new Brick      (gid, oid, geom, attribs_alloc, dom, lid);  break;
    default: errol;
    }
  }
  
  return NULL;
}

void Geometry :: read_addata_VTK (const char *filename, bool sparse)
{
  char *file = new char[255];
  sprintf (file, "%s%s", Pd->give_IN_Path(), filename);
  
  Stream *stream = new Stream();
  stream->open(STRM_void, "r", (const char*&)file);
  
  if (stream->isFile())  FP_skip_line (stream->file());
  
  //* reading
  this->read_VTK ('-', stream, true, sparse);
  
  //* close
  stream->close();
  delete stream;
}

enum VTKrange { VTKR_void, VTKR_points, VTKR_cells, VTKR_lines, VTKR_polygs, VTKR_pred, VTKR_pd, VTKR_cd, VTKR_scals };
void Geometry :: read_VTK (char DATASET, Stream *stream, bool addata, bool sparseflag)
{
  bool lgc = stream->isFile();
  bool mdl = this->isModel();
  
  // variable declaration
  bool status = true, parent;
  int auxi;
  long i, j, id, nn, auxl, noel, cno, offset, nSprsItms=0;
  PoinT coord;
  char ed='-';
  const char *data=NULL, *offs=NULL, *typs, *str;
  char LINE[1023], *WORD=NULL;
  int lLINE=1023;
  Element *elem = NULL;
  VTKrange range;
  const XMLElement *xnel = NULL;
  Stream strm;
  Lmtrx dataL( 1 , ( Pjnts() > Elems() ? Pjnts() : Elems() ) );
  Dmtrx dataD( 1 , ( Pjnts() > Elems() ? Pjnts() : Elems() ) );
  Lvctr *sparse = NULL;
  if (sparseflag) sparse = new Lvctr(Pjnts());
  
  //
  cno = 0;
  
  if (lgc) WORD = new char[255];
  else     xnel = stream->tixel();
  
  if (sparse && !addata)  errol;
  if (addata)  range = VTKR_pred;
  else         range = VTKR_void;
  parent = true;
  while (true) {
    //* key WORD
    if (lgc) {
      if (fgets (LINE, lLINE, stream->file()) == NULL)  break;
      str = LINE;
      SP_scan_word (str, WORD);
    }
    else {
      if (parent) { if (! stream->relink_downF())  errol; else parent = false; }
      else        { if (! stream->relink_next())   break; }
      WORD = (char *)stream->tixel()->Value();
    }
    
    //*  ******************
    //*  ***   POINTS   ***
    //*  ******************
    if      (strcasecmp(WORD, "Points") == 0) {
      if (range == VTKR_void)  range = VTKR_points;  else errol;
      
      // number of nodes
      if (lgc) {
    SP_scan_number (str, noel);
    SP_scan_expected_word2_exit  (str, "double", "float", "Invalid structure of the given VTK file, unsupported datatype of POINTS", CASE);
      }
      else {
    xnel->ToElement()->QueryLongAttribute("NumberOfPoints", &noel);
    data = XP_giveDAtext (stream->tixel(), 1, true, "ascii", "Float32", NULL, &(auxi=3));
      }
      
      Pjnts.reallocup(noel);
      
      // coords
      while (noel) { noel--;
    if (lgc) status += coord.scan_xyz(stream->file());
    else     status += coord.scan_xyz(data);
    Pjnts.add( allocate_point (Pjnts(), &coord) );
      }
      if (lgc)  FP_skip_line (stream->file());
    }
    else if (strcmp(WORD, "Vertices") == 0  ||  strcmp(WORD, "Verts") == 0) {
      if (!lgc)  xnel->ToElement()->QueryLongAttribute("NumberOfVerts" , &noel);
      else       noel = 1;
      if (!noel)  continue;  else  _errorr("error v");
    }
    //*  *************************************
    //*  ***   UNSTRUCTURED GRID - CELLS   ***
    //*  *************************************
    else if (strcasecmp(WORD, "CELLS") == 0) {
      if (DATASET != 'U') errol;
      if (range == VTKR_points)  range = VTKR_cells;  else errol;
      
      // number of cells
      if (lgc) {
    SP_scan_number (str, noel);
    SP_scan_number (str, cno);
      }
      else {
    xnel->ToElement()->QueryLongAttribute("NumberOfCells" , &noel);
    data = XP_giveDAtext (stream->tixel(), 1, false, "ascii", "Int32", "connectivity", NULL);
    offs = XP_giveDAtext (stream->tixel(), 2, false, "ascii", "Int32", "offsets", NULL);
    typs = XP_giveDAtext (stream->tixel(), 3, true,  "ascii", "UInt8", "types", NULL);
      }
      
      Elems.reallocup(noel);
      
      // types of elems
      if (lgc) {
    FP_skip_line (stream->file(), noel);
    fgets (LINE, lLINE, stream->file());  str=LINE;
    SP_scan_expected_word_exit   (str, "CELL_TYPES", "Invalid structure of the given VTK file", CASE);
    SP_scan_expected_number_exit (str,         noel, "Invalid structure of the given VTK file, the number following CELL_TYPES");
      }
      
      while (noel) { noel--;
    if (lgc)  fscanf (stream->file(),"%ld",&auxl);
    else      status += SP_scan_number (typs, auxl);
    switch (auxl) {
    case  3:      Elems.add( allocate_element(mdl, CG_Line,       Elems(), Elems(), this) );  break;
    case  5:      Elems.add( allocate_element(mdl, CG_Triangle,   Elems(), Elems(), this) );  break;
    case  7:      Elems.add( allocate_element(mdl, CG_Polygon,    Elems(), Elems(), this) );  break;
    case  9:      Elems.add( allocate_element(mdl, CG_Quadrangle, Elems(), Elems(), this) );  break;
    case 12:      Elems.add( allocate_element(mdl, CG_Hexahedron, Elems(), Elems(), this) );  break;
    default: _errorr3 ("Unknown %ld-th CELL_TYPE: %ld", Elems()+1, auxl);
    }
      }
      if (lgc) {
    rewind (stream->file());
    FP_skip_line (stream->file(), 5+Pjnts()+1);
      }
      else offset = 0;
      // elem nodes
      for (i=0; i<Elems(); i++) {
    // nn
    if (lgc) { fscanf (stream->file(),"%ld",&nn); cno = cno - nn - 1; }
    else     { nn = offset;  status += SP_scan_number (offs, offset);  nn = offset - nn; }
    //
    if (Elems[i]->give_classid() == classPolygonMdl)  ((PolygonMdl*)Elems[i])->allocNnNe(nn);
    if (nn != Elems[i]->give_nno())
      _errorr4 ("The number of element(%ld) nodes %ld is not equal to %ld, wrong array offsets", i, auxl, nn);
    // loop
    for (j=0; j<nn; j++) {
      if (lgc) fscanf (stream->file(),"%ld", &auxl);
      else     status += SP_scan_number (data, auxl);
      Elems[i]->set_point(j, auxl);
    }
      }
      if (lgc) FP_skip_line (stream->file(), 1+1+Elems());
    }
    //*  *******************************
    //*  ***   POLYDATA - ELEMENTS   ***
    //*  *******************************
    else if (strcasecmp(WORD, "LINES") == 0  ||  strcasecmp(WORD, "POLYGONS" ) == 0  ||  strcasecmp(WORD, "POLYS") == 0) {
      if (DATASET != 'P') errol;
      if (WORD[0] == 'L')  { if (range == VTKR_points)                         range = VTKR_lines;  else errol; }
      if (WORD[0] == 'P')  { if (range == VTKR_points || range == VTKR_lines)  range = VTKR_polygs; else errol; }
      
      // number of elements
      if (lgc) {
    SP_scan_number (str, noel);
    SP_scan_number (str, cno);
      }
      else {
    if (WORD[0] == 'L')  xnel->ToElement()->QueryLongAttribute("NumberOfLines" , &noel);
    if (WORD[0] == 'P')  xnel->ToElement()->QueryLongAttribute("NumberOfPolys" , &noel);
    if (!noel)  continue;
    data = XP_giveDAtext (stream->tixel(), 1, false, "ascii", "Int32", "connectivity", NULL);
    offs = XP_giveDAtext (stream->tixel(), 2, true, "ascii", "Int32", "offsets", NULL);
    offset = 0;
      }
      
      Elems.reallocup(Elems() + noel);
      
      while (noel) { (noel)--;
    // count of nodes
    if (lgc)   fscanf (stream->file(),"%ld", &nn);
    else     { nn = offset;  status += SP_scan_number (offs, offset);  nn = offset - nn; }
    if (WORD[0] == 'L') {
      if      (nn == 2)  elem = allocate_element(mdl, CG_Line,       Elems(), Elems(), this);
      else if (nn >= 3)  elem = allocate_element(mdl, CG_Line,       Elems(), Elems(), this, nn);
      else    _errorr2 ("Unsupported number of nodes on element \"%ld\"", Elems()+1);
    }
    if (WORD[0] == 'P') {
      if      (nn == 3)  elem = allocate_element(mdl, CG_Triangle,   Elems(), Elems(), this);
      else if (nn == 4)  elem = allocate_element(mdl, CG_Quadrangle, Elems(), Elems(), this);
      else    _errorr2 ("Unsupported number of nodes on element \"%ld\"", Elems()+1);
    }
    // nodes
    for (j=0; j<nn; j++) { 
      if (lgc) fscanf (stream->file(),"%ld", &auxl);
      else     status += SP_scan_number (data, auxl);
      elem->set_point(j, auxl);
    }
    
    Elems.add(elem);
    
    if (lgc)  cno = cno - nn - 1;
      }
      if (lgc)  FP_skip_line (stream->file());
    }
    else if (strcmp(WORD, "Strips") == 0) {
      if (!lgc)  xnel->ToElement()->QueryLongAttribute("NumberOfStrips", &noel);
      else       noel = 1;
      if (!noel) continue;  else  _errorr("error s");
    }
    //*  **********************
    //*  ***   POINT_DATA   ***
    //*  **********************
    else if (strcmp(WORD, "POINT_DATA") == 0  ||  strcmp(WORD, "PointData") == 0) {
      if (range == VTKR_cells || range == VTKR_lines || range == VTKR_polygs || range == VTKR_pred)  range = VTKR_pd;  else errol;
      ed = 'p';
      // head
      if (lgc) { if (sparse) SP_scan_number               (str, nSprsItms);
                 else        SP_scan_expected_number_exit (str, Pjnts(), "The number following POINT_DATA is not equal to Nn"); }
      else     { if (sparse) stream->tixel()->QueryLongAttribute("NumberOfPoints", &nSprsItms);
             parent = true; }
      if (!sparse) nSprsItms = Pjnts();
    }
    //*  *********************
    //*  ***   CELL_DATA   ***
    //*  *********************
    else if (strcmp(WORD, "CELL_DATA") == 0  ||  strcmp(WORD, "CellData") == 0) {
      if (range == VTKR_cells || range == VTKR_lines || range == VTKR_polygs || range == VTKR_pred || range == VTKR_scals)  range = VTKR_cd;  else errol;
      ed = 'c';
      // head
      if (lgc) { if (sparse) SP_scan_number               (str, nSprsItms);
             else        SP_scan_expected_number_exit (str, Elems(), "The number following  CELL_DATA is not equal to Ne"); }
      else     { if (sparse) stream->tixel()->QueryLongAttribute("NumberOfCells", &nSprsItms);
             parent = true; }
      if (!sparse) nSprsItms = Elems();
    }
    //*  **************************
    //*  ***   DATA - SCALARS   ***
    //*  **************************
    else if (strcmp(WORD, "SCALARS") == 0  ||  strcmp(WORD, "DataArray") == 0) {
      if (range == VTKR_pd || range == VTKR_cd || range == VTKR_scals)  range = VTKR_scals;  else errol;
      
      // head
      if (lgc)   SP_scan_word (str, WORD);
      else     { WORD = (char *)stream->tixel()->ToElement()->Attribute("Name");  if (!WORD)  _errorr ("there is no attribute \"Name\""); }
      
      //*  *** DATA for POINTS ***
      if (ed == 'p') {
    if      (strcmp(WORD, "IDS_Prescribed_Values"       ) == 0) { scan_DATA_field (stream, str, 0, nSprsItms, &dataL,   NULL, sparse, WORD);  for (i=0; i<nSprsItms; i++) { id = (sparse ? (*sparse)[i] : i);                   ((Node*)Pjnts[id])->give_pointAttribs()->set_nDOFs_BC(dataL.give_ccols(), dataL.give_ptr2val(i,0)); } }
    else if (strcmp(WORD, "ID_SET_IDS_Prescribed_Values") == 0) { scan_DATA_field (stream, str, 1, nSprsItms, &dataL,   NULL, sparse, WORD);  for (i=0; i<nSprsItms; i++) { id = (sparse ? (*sparse)[i] : i);  if (dataL(i,0))  ((Node*)Pjnts[id])->give_pointAttribs()->set_dofbc_with_ID(                  dataL(i,0)-1 );             } }
    else if (strcmp(WORD, "ID_Boundary_Condition"       ) == 0) { scan_DATA_field (stream, str, 1, nSprsItms, &dataL,   NULL, sparse, WORD);  for (i=0; i<nSprsItms; i++) { id = (sparse ? (*sparse)[i] : i);  if (dataL(i,0))  ((Node*)Pjnts[id])->give_pointAttribs()->add_load    (                  dataL(i,0)-1) ;             } }
    else if (strcmp(WORD, "Full_Hinge"                  ) == 0) { scan_DATA_field (stream, str, 1, nSprsItms, &dataL,   NULL, sparse, WORD);  for (i=0; i<nSprsItms; i++) { id = (sparse ? (*sparse)[i] : i);                   ((Node*)Pjnts[id])->give_pointAttribs()->set_hinge      (               dataL(i,0)  ) ;             } }
    else if (strcmp(WORD, "ada_refined_derivatives"     ) == 0) { scan_DATA_field (stream, str, 0, nSprsItms,   NULL, &dataD, sparse, WORD);  for (i=0; i<nSprsItms; i++) { id = (sparse ? (*sparse)[i] : i);                   ((Node*)Pjnts[id])->set_resultN               (dataD.give_ccols(), dataD.give_ptr2val(i,0), 0, RTN_ada_refder); }      ((Mesh*)this)->set_RTN(RTN_ada_refder, true);    }
    else if (strcmp(WORD, "Property"                    ) == 0) { scan_DATA_field (stream, str, 1, nSprsItms, &dataL,   NULL, sparse, WORD);  for (i=0; i<nSprsItms; i++) { id = (sparse ? (*sparse)[i] : i);  if (dataL(i,0))  ((Node*)Pjnts[id])->set_mprop                 (                  dataL(i,0)  ) ;             } }
    else  _errorr2 ("Invalid name of SCALARS/DataArray for POINT_DATA/PointData \"%s\"", WORD);
      }
      
      //*  *** DATA for CELLS ***
      if (ed == 'c') {
    if      (strcmp(WORD, "Virtual"                     ) == 0) { if (!mdl) errol;  scan_DATA_field (stream, str, 1, nSprsItms, &dataL,   NULL, sparse, WORD);  for (i=0; i<nSprsItms; i++) { id = (sparse ? (*sparse)[i] : i);                        ((GelemAttribs*)Elems[id]->give_elemAttribs())->set_virtual     (                      dataL             (i,0)  ) ; } }
    else if (strcmp(WORD, "Property"                    ) == 0) { elemprop = true;  scan_DATA_field (stream, str, 1, nSprsItms, &dataL,   NULL, sparse, WORD);  for (i=0; i<nSprsItms; i++) { id = (sparse ? (*sparse)[i] : i);  if (dataL(i,0))                       Elems[id]->set_mprop                      (                      dataL             (i,0)  ) ; } }
    else if (strcmp(WORD, "IDS_Prescribed_Values"       ) == 0) { if (!mdl) errol;  scan_DATA_field (stream, str, 0, nSprsItms, &dataL,   NULL, sparse, WORD);  for (i=0; i<nSprsItms; i++) { id = (sparse ? (*sparse)[i] : i);                        ((GelemAttribs*)Elems[id]->give_elemAttribs())->set_nDOFs_BC    (dataL.give_ccols(),   dataL.give_ptr2val(i,0)  ) ; } }
    else if (strcmp(WORD, "ID_SET_IDS_Prescribed_Values") == 0) { if (!mdl) errol;  scan_DATA_field (stream, str, 1, nSprsItms, &dataL,   NULL, sparse, WORD);  for (i=0; i<nSprsItms; i++) { id = (sparse ? (*sparse)[i] : i);  if (dataL(i,0))       ((GelemAttribs*)Elems[id]->give_elemAttribs())->set_dofbc_with_ID(                     dataL             (i,0)-1) ; } }
    else if (strcmp(WORD, "ID_Cross-Section"            ) == 0) {                   scan_DATA_field (stream, str, 1, nSprsItms, &dataL,   NULL, sparse, WORD);  for (i=0; i<nSprsItms; i++) { id = (sparse ? (*sparse)[i] : i);  if (dataL(i,0))                       Elems[id]->give_elemAttribs() ->set_cs          (                      dataL             (i,0)-1) ; } }
    else if (strcmp(WORD, "ID_Material"                 ) == 0) {                   scan_DATA_field (stream, str, 1, nSprsItms, &dataL,   NULL, sparse, WORD);  for (i=0; i<nSprsItms; i++) { id = (sparse ? (*sparse)[i] : i);  if (dataL(i,0))                       Elems[id]->give_elemAttribs() ->set_mat         (                      dataL             (i,0)-1) ; } }
    else if (strncmp(WORD, "ID_Boundary_Condition", 21  ) == 0) {                   scan_DATA_field (stream, str, 1, nSprsItms, &dataL,   NULL, sparse, WORD);  for (i=0; i<nSprsItms; i++) { id = (sparse ? (*sparse)[i] : i);  if (dataL(i,0))                       Elems[id]                     ->set_load        (                      dataL             (i,0)-1) ; } }
    else if (strcmp(WORD, "symStiffMtrxOfMat"           ) == 0) {                   scan_DATA_field (stream, str, 0, nSprsItms,   NULL, &dataD, sparse, WORD);  for (i=0; i<nSprsItms; i++) { id = (sparse ? (*sparse)[i] : i);                        ((GelemAttribs*)Elems[id]->give_elemAttribs())->set_MatStiffMtrx(dataD.give_ccols(),   dataD.give_ptr2val(i,0)  ) ; } }
    else if (strcmp(WORD, "LCS_xy_vector"               ) == 0) {                   scan_DATA_field (stream, str, 3, nSprsItms,   NULL, &dataD, sparse, WORD);  for (i=0; i<nSprsItms; i++) { id = (sparse ? (*sparse)[i] : i);                                        Elems[id]->give_elemAttribs() ->set_lcs         (LCSd_xy, LCSt_Vector, dataD.give_ptr2val(i,0)  ) ; } }
    else if (strcmp(WORD, "LCS_xz_vector"               ) == 0) {                   scan_DATA_field (stream, str, 3, nSprsItms,   NULL, &dataD, sparse, WORD);  for (i=0; i<nSprsItms; i++) { id = (sparse ? (*sparse)[i] : i);                                        Elems[id]->give_elemAttribs() ->set_lcs         (LCSd_xz, LCSt_Vector, dataD.give_ptr2val(i,0)  ) ; } }
    else if (strcmp(WORD, "MeshGen_elemSize"            ) == 0) { if (!mdl) errol;  scan_DATA_field (stream, str, 1, nSprsItms,   NULL, &dataD, sparse, WORD);  for (i=0; i<nSprsItms; i++) { id = (sparse ? (*sparse)[i] : i);  if (dataD(i,0) > 0.0) ((Gelement*)    Elems[id]                    )->set_elemSize    (                      dataD             (i,0)  ) ; } }
    else if (strcmp(WORD, "MeshGen_elemCount"           ) == 0) { if (!mdl) errol;  scan_DATA_field (stream, str, 1, nSprsItms, &dataL,   NULL, sparse, WORD);  for (i=0; i<nSprsItms; i++) { id = (sparse ? (*sparse)[i] : i);  if (dataL(i,0) > 0)   ((Gelement*)    Elems[id]                    )->set_elemCount   (                      dataL             (i,0)  ) ; } }
    else if (strcmp(WORD, "ada_error_pct_l"             ) == 0) {                   scan_DATA_field (stream, str, 0, nSprsItms,   NULL, &dataD, sparse, WORD);  for (i=0; i<nSprsItms; i++) { id = (sparse ? (*sparse)[i] : i);                        ((FElement*)    Elems[id]                    )->set_result      (                      dataD             (i,0), 0, RTE_ada_error_pct_l); } ((Mesh*)this)->set_RTE(RTE_ada_error_pct_l, true);  }
    else if (strcmp(WORD, "ada_refsizel_relative"       ) == 0) {                   scan_DATA_field (stream, str, 0, nSprsItms,   NULL, &dataD, sparse, WORD);  for (i=0; i<nSprsItms; i++) { id = (sparse ? (*sparse)[i] : i);                        ((FElement*)    Elems[id]                    )->set_result      (                      dataD             (i,0), 0, RTE_ada_h_new); } ((Mesh*)this)->set_RTE(RTE_ada_h_new, true);  }
    else _errorr2 ("Invalid name of SCALARS/DataArray in CELL_DATA/CellData of the given VTK XML file, \"%s\"", WORD);
      }
      if (!lgc)  if (stream->tixel()->NextSibling() == NULL)  stream->relink_up();
    }
    else if (strcmp(WORD, "") == 0) {;}
    else _errorr2 ("Invalid structure of the given VTK file, unexpected key word \"%s\"", WORD);
    
    if (lgc) { if (cno) _errorr2 ("cno(%ld) is not zero", cno); }
    else {
      if (data && data[0] != '\0')  _errorr2("CellData of name \"%s\" has too many numbers", WORD);
      if (offs && offs[0] != '\0')  _errorr("error");
    }
  }
  
  if (Pjnts() == 0)  _errorr ("Invalid data in the given VTK file, number of points is zero");
  if (Elems() == 0)  _errorr ("Invalid data in the given VTK file, number of cells is zero");
  
  if (lgc) delete [] WORD;
  if (sparseflag) delete sparse;
  if (status == false)  _errorr("error");
}

void Geometry :: join_geom (Geometry *geom)
{
  Pjnts.reallocplus(Pjnts() + geom->give_Npts());
  Edges.reallocplus(Edges() + geom->give_Neds());
  Faces.reallocplus(Faces() + geom->give_Nfas());
  Elems.reallocplus(Elems() + geom->give_Nels());
  
  long i;
  for (i=0; i<geom->give_Npts(); i++)  this->add_another_Pjnt(const_cast<Point*  >(geom->give_Pjnt(i)));
  for (i=0; i<geom->give_Neds(); i++)  this->add_another_Edge(const_cast<Edge*   >(geom->give_Edge(i)));
  for (i=0; i<geom->give_Nfas(); i++)  this->add_another_Face(const_cast<Face*   >(geom->give_Face(i)));
  for (i=0; i<geom->give_Nels(); i++)  this->add_another_Elem(const_cast<Element*>(geom->give_Elem(i)));
}


bool Geometry :: give_3Delement_with_point_inside (const PoinT *point, const GeometryComponent *&comp, PoinT *nc) const
{
  if ( ! connectivity_is_assembled() )  _errorr("connectivity is not initiated");
  
  long i, j, mind, numsuperelem;
  double dist, distance, *distances;
  const Point *node = NULL;
  const Element *elem;
  
  // finds the most closes node to the point
  distance = 1.0e100;
  for (i=0; i<Pjnts(); i++) {
    dist = point->dist2_to(Pjnts[i]->give_coords());
    if (dist < distance && Pjnts[i]->is_on_3d_element()) {
      distance = dist;
      node = Pjnts[i];
    }
  }
  if (!node) _errorr0;
  //
  numsuperelem = node->give_numsuperelem();
  
  distances = new double [numsuperelem];
  
  // computing of suma of distances startpoint-node[i] for every element superior to closednode
  for (i=0; i<numsuperelem; i++) {
    distances[i] = 0.0;
    elem = node->give_superelem(i);
    for (j=0; j<elem->give_nno(); j++)
      distances[i] += point->dist2_to (elem->give_point(j)->give_coords());
  }
  
  for (i=0; i<numsuperelem; i++) {
    // finding the least distance
    mind = 0;
    for (j=1; j<numsuperelem; j++)
      if (distances[mind] > distances[j])  mind = j;
    
    distances[mind] = 1.0e100;
    elem = node->give_superelem(mind);
    
    if (elem->give_dimension() != 3)  continue;
    
    if ( ((FElement*)elem)->is_point_on (point, comp, nc) ) {
      delete [] distances;
      return true;
    }
  }
  
  delete [] distances;
  return false;
}




//*  ********************************************************************************************
//*  *** *** *** ***                         CLASS MODEL                          *** *** *** ***
//*  ********************************************************************************************
Model :: Model (long gid, const Problem *p, MMprocessing mmp) : Geometry (gid, p, mmp)
{

}
Model :: ~Model (void)
{

}


//*  **********************************************
//*  ***   READING OF REINFORCEMENT POLYLINES   ***
//*  **********************************************
void Model :: read_model_polylines (const char *path, const char *filename)
{
  FILE *in = _openFilePN("r", "PolyLine", path, filename);
  
  elemprop = true;
  
  long i, j, id, nv, ne, prop;
  PoinT coords;
  
  fscanf (in, "%ld", &ne);
  Elems.reallocup(ne);
  
  PolyLine *pl;
  
  for (i=0; i<ne; i++) {
    fscanf (in, "%ld", &id);  id--; if (i != id)  _errorr("mismatch in polyline file");
    fscanf (in, "%ld", &nv);
    fscanf (in, "%ld", &prop);
    
    pl = new PolyLine(id, id, this, nv);
    Pjnts.reallocup(Pjnts()+nv);
    pl->set_mprop(prop);
    
    for (j=0; j<nv; j++) {
      coords.scan_xyz(in);
      Pjnts.add( new Vertex (this, Pjnts(), &coords, 'n') );
      Pjnts.last()->set_mprop(Pjnts());
      pl->set_point(j, Pjnts.last());
    }
    
    Elems.add(pl);
  }
  
  fclose (in);
}



//*  **********************************************
//*  ***   READING OF REINFORCEMENT POLYLINES   ***
//*  **********************************************
struct Predp {
  char name[15];
  int id, nl, nk;
  GPA<PoinT> points;

  double length (void) { return points.last()->x; }
  
  int is_my_ch(char *chn) {
    return (strncmp(name, chn, strlen(name)) == 0 ? strlen(name) : 0);
  }
};
void Model :: read_model_MELNIK (void)
{
  FILE *in;
  PolyLine *pl;
  char filename[255] = "";
  char name[15], c;
  long auxl, id, i, nno;
  double auxd;
  //long i, j, id, nv, ne, prop;
  PoinT coords, *auxP;
  GPA<Predp> prdi;
  Predp *prd;
  
  int quart;
  int right;
  bool propsp;
  
  switch (Pd->give_melnik()) {
  case 10:  quart = 4;  propsp = false;  break;  // celek
  case 20:  quart = 2;  propsp = false;  break;  // polovina
  case 21:  quart = 2;  propsp = true;   break;  // polovina + specialni property
  case 40:  quart = 1;  propsp = false;  break;  // ctvrtina
  default: _errorr("\"-P_melnik x\" missing");
  }
  
  elemprop = true;
  
  for (int f=0; f<7; f++) {
    if (f == 0) { if (1             ) sprintf (filename, "this.model.1_dolni_krajni_L_"); else continue; }
    if (f == 1) { if (1 && quart > 1) sprintf (filename, "this.model.1_dolni_krajni_P_"); else continue; }
    if (f == 2) { if (1             ) sprintf (filename, "this.model.2_dolni_stredni_");  else continue; }
    if (f == 3) { if (1             ) sprintf (filename, "this.model.3_horni_konzol_L_"); else continue; }
    if (f == 4) { if (1 && quart > 1) sprintf (filename, "this.model.3_horni_konzol_P_"); else continue; }
    if (f == 5) { if (1             ) sprintf (filename, "this.model.4_horni_zmonol_L_"); else continue; }
    if (f == 6) { if (1 && quart > 1) sprintf (filename, "this.model.4_horni_zmonol_P_"); else continue; }
    
    //* predpeti
    in = _openFilePNS("r", "Melnik file", "", filename, "P.txt");
    
    while (true) {
      prd = new Predp;
      
      if (fscanf (in, "%d", &prd->id) == EOF)  break;
      FP_scan_word (in, prd->name);
      fscanf (in, "%d", &prd->nl);
      fscanf (in, "%d", &prd->nk);
      
      while (true) {
    fscanf (in, "%lf", &auxd);
    
    auxP = new PoinT;
    auxP->x = auxd;
      
    prd->points.add(auxP);
    
    FP_skip_space (in);
    if ((c = fgetc(in)) == '\n')  break;
    else ungetc(c,in);
      }
      
      for (i=0; i<prd->points()-1; i++) {
    prd->points[i]->scan_y(in);
    prd->points[i]->y *= -1.0;
      }
      FP_skip_space (in);
      if (fgetc(in) != '\n')  errol;
      
      prdi.add(prd);
    }
    
    fclose (in);
    
    
    //* kabely
    for (int kz=0; kz<2; kz++) {
      if (kz==0)                         in = _openFilePNS("r", "Melnik file", "", filename, "K.txt");
      else if (quart < 4) continue; else in = _openFilePNS("r", "Melnik file", "", filename, "Z.txt");
      
      pl = new PolyLine(Elems(), Elems(), this, 0);
      id = -1;
      while (true) {
    c = fscanf (in, "%ld", &auxl);  auxl--;
    
    if (id == -1) { if (auxl != 0)  errol; }
    else if (auxl == 0 || c == EOF) {
      auxl = 0;
      for (i=0; i<prdi(); i++)
        if (prdi[i]->is_my_ch(name) > auxl) {
          prd = prdi[i];
          auxl = prdi[i]->is_my_ch(name);
        }
      
      if (!propsp) pl->set_mprop(prd->nl);                  // pro testovani
      else         pl->set_mprop(prd->nl*1000+prd->id*10);  // aktualni vystup
      //pl->set_prop_node_inher(true);
      prd->nk--;
      
      //* kontrola
      if (0) {
        nno = pl->give_nno();
        if (nno != prd->points()-1 && !(f==2 && quart==1) )
          _errorr6("f=%ld kz=%d  %s %ld != %ld", f, kz, name, nno, prd->points()-1);
        auxd = (pl->give_lav() - prd->length());
        if ( (auxd > 0.1  ||  auxd < -0.1) && !(f==2 && quart==1) )
          _warningg5("f=%ld  %s length %lf != %lf", f, name, pl->give_lav(), prd->length());
        //if ( (pl->give_lav()/1000.0 - prd->length()) > 0.1 )  errol;
        if (f==1 || f==4 || f==6) {
          if (quart==1)  errol;
          for (i=0; i<nno; i++) {
        if ( fabs(pl->give_point(i)->give_coord(0) - prd->points[nno-1-i][0].x) > 0.001)
          errol;
        if ( fabs(pl->give_point(i)->give_coord(2) - prd->points[nno-1-i][0].y) > 0.001)
          _warningg5("f=%ld  %s z: %lf != %lf", f, name, pl->give_point(i)->give_coord(2), prd->points[nno-1-i][0].y);
          }
        }
        else {
          for (i=0; i<nno; i++) {
        if ( fabs(pl->give_point(i)->give_coord(0) - prd->points[i][0].x) > 0.001 && !(f==2 && quart==1))
          errol;
        if ( fabs(pl->give_point(i)->give_coord(2) - prd->points[i][0].y) > 0.001)
          _warningg5("f=%ld  %s z: %lf != %lf", f, name, pl->give_point(i)->give_coord(2), prd->points[i][0].y);
          }
        }
      }
      
      Elems.add(pl);
      
      if (c == EOF) { pl = NULL; break; }
      
      pl = new PolyLine(Elems(), Elems(), this, 0);
      id = -1;
    }
    id ++;
    if (id == 0)  right = 0;
    
    coords.scan_xyz(in);
    coords.z *= -1.0;
    coords.dvd(1000.0);
    if (quart==1 && coords.x > 146.0) {
      coords.x = 146.0;
      right++;
    }
    
    // DOCASNE prekryvani vyztuze L a P pro provnani a zjisteni symetricnosti
    //if (coords.x > 146.0) {
    //  coords.x = -1 * ( coords.x - 292.0);
    //}
    
    if (right < 2) {
      Pjnts.add( new Vertex (this, Pjnts(), &coords, 'n'));
      
      pl->add_point(Pjnts.last());
    }
    
    if (id)  FP_scan_expected_word (in, name);
    else     FP_scan_word          (in, name);
    
    FP_skip_space (in);
    if (fgetc(in) != '\n')  errol;
    
    //if (right == 1)  right = 2;
      }
      
      fclose (in);
    }
    
    if (quart==4)
      for (i=0; i<prdi(); i++)
    if (prdi[i]->nk != 0)
      _warningg4("f=%ld predp skup %ld, %d != 0", f, i, prdi[i]->nk);
    
    prdi.delete_objects();
    prdi.resize(0);
  }
}

//void Model :: findRANofSkewLines (void)
//{
//  double p1[3], p2[3];
//  
//  int a = give_edge_points_commperpend_skewlines ( Elems[0]->give_point(0)->give_coords(),
//                         Elems[0]->give_point(1)->give_coords(),
//                         Elems[1]->give_point(0)->give_coords(),
//                         Elems[1]->give_point(1)->give_coords(),
//                         p1, p2);
//  
//  printf ("p1 %lf %lf %lf\n", p1[0], p1[1], p1[2]);
//  printf ("p2 %lf %lf %lf\n", p2[0], p2[1], p2[2]);
//  printf ("a = %d\n", a);
//}

void Model :: transform_to_mesh (void)
{
  MMprocessing gp = this->give_gp();
  
  switch (gp) {
  case MMP_primary:
  case MMP_secondary: {
    //* new mesh
    Mesh *mesh = new Mesh(0, Pd, 1, gp);
    mesh->set_Master_model(this);
    
    //* transformation
    if (Pd->give_IN_meshGen_elemSize() == 0.0)
      mesh->read_mesh_equal_to_model();
    else {
      this->generate_mesh_primary();
      mesh->read_mesh_T3d(NULL, NULL);
    }
    
    if (gp == MMP_primary && Pd->give_cMesh())  ((Problem*)Pd)->wedge_mesh(0, mesh);
    else                                        ((Problem*)Pd)->add_mesh(mesh);
    
    mesh->initialize();
    
    break;
  }
  case MMP_RFbyHN:
    this->generate_mesh_RFbyHN();
    break;
  default: errol;
  }
}


void Model :: generate_mesh_RFbyHN (void)
{
  GPA<const Node> hnodes;        // (hanging)nodes on polyline
  const PolyLine *pol;
  PolyLine *newpol;
  Mesh* pm = Pd->give_Mesh(0);   // parent mesh geometry
  
  pm->connectivity_assembling ();
 
  //* loop over elements = RF bars
  for (long i=0; i<Elems(); i++) {
    //for (long i=37; i<38; i++) {
    pol = dynamic_cast<const PolyLine *>(Elems[i]);
    if (!pol)  errol;
    
    newpol = pol->generate_mesh_RFbyHN (pm, hnodes);
    
    if (newpol) {
      if (i+1 == Elems())  this->add_another_Elem(newpol);
      else                 this->wdg_another_Elem(newpol, i+1);
      
      hnodes.resize(1);
    }
    else
      hnodes.resize(0);
  }
}


void Model :: generate_mesh_primary (void)
{
  //* info print
  double sec = clock();
  Pd->sodriver()->print_info_message (SODE_start_blue_arrow, "Model[%ld] - mesh generation by T3d", this->give_ID());
  
  this->connectivity_assembling ();
  
  switch (Pd->give_P_mesher()) {
  case MGT_T3d: {
    femFileFormat fff = FFF_T3d;
    
    //* T3D INPUT FILE
    FILE *out = _openFilePNS ("w", "T3d input", Pd->give_OUT_moFILE()->give_path(), Pd->give_OUT_moFILE()->give_base(), ".T3d.in");
    
    long i;
    for (i=0; i<Pjnts(); i++)  Pjnts[i]->print_row(out, fff);  fprintf (out,"\n");
    for (i=0; i<Edges(); i++)  Edges[i]->print_row(out, fff);  fprintf (out,"\n");
    for (i=0; i<Faces(); i++)  Faces[i]->print_row(out, fff);  fprintf (out,"\n");
    for (i=0; i<Elems(); i++)  Elems[i]->print_row(out, fff);
    
    fclose (out);
    
    
    const char *basename = Pd->give_OUT_moFILE()->give_name();
    char cmd[255] = "";
    
    //if (IN_resultfile != NULL)  errol;
    //IN_resultfile = new FiLe(FFF_OOFEM, this->OUT_moFILE->give_path(), this->OUT_moFILE->give_base(), ".oofem.out");
    
    //sprintf (cmd, "-i %s.T3d.in -o %s.T3d.out -d %.2lf -p 8 -g %s.T3d.native.log",   basename, basename, Pd->give_IN_elemsize(), basename);
    sprintf (cmd, "-i %s.T3d.in -o %s.T3d.out -d %.2lf -p 8",   basename, basename, Pd->give_IN_meshGen_elemSize());
    //sprintf (cmd, "-i %s.T3d.in -o %s.T3d.out -p 8",   basename, basename);
    
    //sprintf (cmd, "-i %s.T3d.in -o %s.T3d.out -d %.2lf -p 264", basename, basename, Pd->give_IN_elemsize());
    
    if (Pd->give_IN_file_bgm())
      sprintf (cmd, "%s -m %s",  cmd, Pd->give_IN_file_bgm()->give_name());
    
    if (Pd->give_P_mesherbinary())  mesh_generate_T3d (Pd, Pd->give_P_mesherbinary()->give_name_or_null(), cmd);
    else                            mesh_generate_T3d (Pd, NULL                                          , cmd);
    
    //fprintf (stdout, "%s\n", cmd);
    break;
  }
  default: errol;
  }
  
  Pd->sodriver()->print_info_time_green_ok (SODE_end_green_ok, (clock() - sec) / (double)CLOCKS_PER_SEC);
}


//*  ********************************************************************************************
//*  *** *** *** ***                         CLASS FEMESH                         *** *** *** ***
//*  ********************************************************************************************
Mesh :: Mesh (long gid, const Problem *p, long numdom, MMprocessing mmp) : Geometry (gid, p, mmp)
{
  Master = NULL;
  
  //*  *** SET sequent or parallel version ***
  NumDomains = numdom;
  if      (NumDomains==1) Parallel = false;
  else if (NumDomains>1)  Parallel = true;
  else    _errorr ("Wrong number of domains");
  
  if (Parallel) {
    NumNodes = new long[NumDomains];
    NumElems = new long[NumDomains];
  }
  else { NumNodes = NumElems = NULL; }
  
  cSubDomains = 0;
  //* RESULTS
  rslts_solver = SOL_Void;
  rslts_nsteps = 0;
  rslts_step = NULL;
  rslts_loadlevel= NULL;
  memset(RTN, 0, cRTN*sizeof(bool));
  memset(RTE, 0, cRTE*sizeof(bool));
  critical_coeff = 0.0;
  //* Adaptivity
  nRegions = 0;
}
Mesh :: ~Mesh (void)
{
  if (Parallel) {
    deallocateCheck (NumNodes);
    deallocateCheck (NumElems);
  }
  
  delete rslts_step;
  delete rslts_loadlevel;
  
  Lata.delete_objects();
  Data.delete_objects();
}

void Mesh :: initialize (void)
{
  Geometry :: initialize ();
  
  //* divide quadrangles into two triangles
  if (Pd->give_PDBO(PDBO_divideQuads)) {
    Triangle *tr1, *tr2;
    Quadrangle *quad;
    long ne = Elems();
    for (long i=0; i<ne; i++) {
      quad = dynamic_cast<Quadrangle *> (Elems[i]);
      if (quad == NULL)  continue;
      
      bool fds = quad->is_first_diag_short();
      tr1 = new Triangle(quad, true,  fds);
      tr2 = new Triangle(quad, false, fds);
      
      quad->connectivity_assembling();
      tr1 ->connectivity_assembling();
      tr2 ->connectivity_assembling();
      
      ((Face*)tr1->give_face(0))->give_facedgeAttribs()->set_new (quad->give_face(0)->give_facedgeAttribs());
      ((Face*)tr2->give_face(0))->give_facedgeAttribs()->set_new (quad->give_face(0)->give_facedgeAttribs());
      
      quad->connectivity_removing();
      delete quad;
      Elems.assign(i, tr1);
      this->add_another_Element (tr2, 0);
    }
    if (Elems() - ne) {
      fprintf (stdout,"\n    Divided quadrangles   %ld\n", Elems() - ne);
      // toto smaz: print_info_message ('S',ETC_DEFAULT,"");
    }
  }
}

void Mesh :: checkConsistency (void) const
{
  Geometry :: checkConsistency();
}

void Mesh :: add_another_Element (FElement *elem, long dom)
{
  elem->set_lid(Parallel ? NumElems[dom]++ : Elems());
  
  this->add_another_Elem (elem);
}

void Mesh :: mesh_quality (void) const
{
  ((Mesh*)this)->connectivity_assembling ();
  
  if ( ! Pd->give_PDBO (PDBO_OUT_MeshQual) )  return;
  
  Pd->sodriver()->print_info_message (SODE_start_blue_arrow, "Compute element quality");
  
  long i;
  double q;
  long Ne = Elems();
  Lvctr Q(7), Qi(Ne);
  
  Q.zero();
  
  for (i=0; i<Ne; i++) {
    q = Elems[i]->give_quality();
    
    if      (q < 0.0 )  { Q[0]++; Qi[i] = 7; }
    else if (q < 0.01)  { Q[1]++; Qi[i] = 6; }
    else if (q < 0.2 )  { Q[2]++; Qi[i] = 5; }
    else if (q < 0.4 )  { Q[3]++; Qi[i] = 4; }
    else if (q < 0.6 )  { Q[4]++; Qi[i] = 3; }
    else if (q < 0.8 )  { Q[5]++; Qi[i] = 2; }
    else if (q < 1.0 )  { Q[6]++; Qi[i] = 1; }
    else _errorr2 ("Quality %e > 1.0", q);
  }
  
  // check
  long sum = 0;
  for (i=0; i<7; i++)  sum += Q[i];
  if (sum != Ne) _errorr ("sum != Ne");
  // check
  
  if (true) {
    fprintf (stdout,"\n    Element quality, %ld%% of %ld\n", (100*(Ne-Q[0]))/Ne, Ne);
    
    if (Q[0]<Ne) {
      fprintf (stdout,"    %5ld (%3ld%%) [1] fine           (0.8 - 1.0)\n", Q[6], (100*Q[6])/Ne);
      fprintf (stdout,"    %5ld (%3ld%%) [2] good           (0.6 - 0.8)\n", Q[5], (100*Q[5])/Ne);
      fprintf (stdout,"    %5ld (%3ld%%) [3] low            (0.4 - 0.6)\n", Q[4], (100*Q[4])/Ne);
      fprintf (stdout,"    %5ld (%3ld%%) [4] wrong          (0.2 - 0.4)\n", Q[3], (100*Q[3])/Ne);
      fprintf (stdout,"    %5ld (%3ld%%) [5] ugly           (0.01- 0.2)\n", Q[2], (100*Q[2])/Ne);
      fprintf (stdout,"    %5ld (%3ld%%) [6] non-acceptable (0.0 -0.01)\n", Q[1], (100*Q[1])/Ne);
      fprintf (stdout,"    %5ld (%3ld%%) [7] --             (    < 0.0)\n", Q[0], (100*Q[0])/Ne);
    }
    // toto smaz: print_info_message ('S',ETC_DEFAULT,"");
  }
  
  if      (Q[0] == Ne) Pd->sodriver()->print_info_message        (SODE_end,              "--");
  else if (Q[1] != 0)  Pd->sodriver()->print_info_message_colour (SODE_end, ETC_D_RED,   "NON-ACCEPTABLY LOW QUALITY");
  else                 Pd->sodriver()->print_info_message_colour (SODE_end, ETC_D_GREEN, "ok");
  
  //this->print_VTK_CD_scal (true, ".mesh_quality", "mesh_quality", (const Xvctr**)&Qi);
}


void Mesh :: RIGIDmatToRAN (void)
{
  this->connectivity_assembling ();
  
  Pd->sodriver()->print_info_message (SODE_start_blue_arrow, "%s%ld] - RIGID elems to RigidArm switch", (this->isModel() ? "Model[" : "Mesh["), this->give_ID());
  
  long i, nd=0;
  
  for (i=0; i<Elems(); i++) {
    if (Elems[i]->give_elemAttribs()->give_mat() == NULL  ||  Elems[i]->give_elemAttribs()->give_mat()->give_type() != MAT_RIGID)  continue;
    
    //* check
    if (Elems[i]->give_classid() != classBeam)  _errorr("RigidBody material is not at beam element");
    if (Elems[i]->give_nno()     != 2)          _errorr("RigidBody material is not at beam element with 2 nodes");
    
    if ( this->FElems(i)->give_mdl_masterel() && 
     this->FElems(i)->give_mdl_masterel()->give_lav() != this->Elems[i]->give_lav() )
      _errorr("Element with rigid body mat. is not same length as master");
    
    const Node *n1 = dynamic_cast<const Node *> (Elems[i]->give_point(0));
    const Node *n2 = dynamic_cast<const Node *> (Elems[i]->give_point(1));
    
    if (n1->give_classid() != classNode)  _errorr("First node at element with RigidBody material is not Node");
    if (n2->give_classid() != classNode)  _errorr("First node at element with RigidBody material is not Node");
    
    if (n1->give_pointAttribs()->give_nDOFs() != n2->give_pointAttribs()->give_nDOFs())  _errorr("Not same Number of DOFS at nodes of first and second node at element with RigidBody material");
    
    //* delete element
    Elems[i]->connectivity_removing();
    Elems.delete_object(i);   nd++;
    
    //* make rigid arm node
    RigidArmNode *ran = new RigidArmNode (n1);
    ((RANAttribs*)ran->give_pointAttribs())->complete_setup_yourself(n2);
    
    //* smazat n1 a nahradit ran
    this->replace_Pjnt_by ((Node*)n1, ran); //
  }
  
  // shake
  TF_GPA_shake_down_reidoid(Elems);
  
  // rm RIGID mat
  ((Problem*)Pd)->rm_MAT_RIGID();
  
  Pd->sodriver()->print_info_message (SODE_end, "%d", nd);
}


//*  ************************************
//*  ***   READING OF GEOMETRY MESH   ***
//*  ************************************
void Mesh :: read_mesh_ANSYS (char *path, const char *filename)
{
  // open file
  FILE *in = _openFilePN("r", "Mesh", path, filename);
  
  // variable declaration
  char c;
  long i,id,nl,domid;
  PoinT coord;
  FElement *elem;
  //
  domid=0;
  
  
  // GET NUMBER OF LINES
  nl = FP_number_of_lines (in, true);
  
  Pjnts.reallocup(nl);
  Elems.reallocup(nl);
  
  
  //*  ***  HEAD  ***
  FP_skip_line (in,4);
  
  // main loop
  while ((c = getc(in)) != EOF) {
    if (c=='\n') continue;
    
    switch (c){
    case 'K':{ // node
      //while ( getc(in) != ',') ;
      ;                      if ( getc(in) != ',')  _errorr3 ("There is no separation char ',' nn=%ld, nel=%ld", Pjnts(), Elems());
      fscanf (in,"%ld",&id); if ( getc(in) != ',')  _errorr3 ("There is no separation char ',' nn=%ld, nel=%ld", Pjnts(), Elems());
      coord.scan_x(in);      if ( getc(in) != ',')  _errorr3 ("There is no separation char ',' nn=%ld, nel=%ld", Pjnts(), Elems());
      coord.scan_y(in);      if ( getc(in) != ',')  _errorr3 ("There is no separation char ',' nn=%ld, nel=%ld", Pjnts(), Elems());
      coord.scan_z(in);
      
      // CTENI COORDS HODIT DO POINT.cpp read input
      Pjnts.add( allocate_point (id-1, &coord) );
      break;
    }
    case 'A':{ // element - triangle
      elem = new Triangle (Elems(),Elems(), this, true, domid, Elems());
      
      for (i=0; i<elem->give_nno(); i++) {
    if ( getc(in) != ',')  _errorr3 ("There is no separation char ',' nn=%ld, nel=%ld", Pjnts(), Elems());
    fscanf (in,"%ld", &id);
    elem->set_node(i, id-1);
      }
      Elems.add(elem);
      break;
    }
    default:
      FP_skip_line(in);
    }
  }
  
  fclose (in);
}

//*  ************************************
//*  ***   READING OF GEOMETRY MESH   ***
//*  ************************************
void Mesh :: read_mesh_JKTK (char *path, const char *filename)
{
  // open file
  FILE *in = _openFilePN("r", "Mesh", path, filename);
  
  // variable declaration
  long i, auxl, n;
  const char *str;
  char LINE[1023];
  femFileFormat fff = FFF_JKTK;
  
  //* nodes
  fscanf (in, "%ld", &n);
  Pjnts.reallocup(n);
  for (i=0; i<n; i++) {
    FP_scan_line_skip_emptyORcommented (in, LINE);  str=LINE;
    
    SP_scan_expected_number_exit (str, i+1, "Invalid structure of file with JKTK format: wrong node id");
    Pjnts.add((Point*)new Node (this, i, NULL));
    Nodes(i)->read_input (str, fff);
  }
  
  //* connectivity will be assembled
  this->connectivity_allocation();
  
  //* elements
  fscanf (in, "%ld", &n);
  Elems.reallocup(n);
  for (i=0; i<n; i++) {
    FP_scan_line_skip_emptyORcommented (in, LINE);  str=LINE;
    
    SP_scan_expected_number_exit (str, i+1, "Invalid structure of file with JKTK format: wrong node id");
    SP_scan_number (str, auxl);
    Elems.add( allocate_element(false, CellGeometry_i2e_JKTK(auxl), Elems(), Elems(), this) );
    
    // if quadratic element -> edit attributes  !!!!!
    
    FElems(i)->read_input (str, fff);
  }
  
  if (EOF != fscanf (in, "%ld", &auxl))  errol;
  
  fclose (in);
}


//*  ************************************
//*  ***   READING OF GEOMETRY MESH   ***
//*  ************************************
void Mesh :: read_mesh_T3d (char *path, const char *filename)
{
  FILE *in;
  long i, id, lid, deg, auxl, domid, trash, el[7], cND, outspec, nno;
  bool auxb;
  PoinT coord;
  int t3dout;
  CellGeometry cg;
  
  const char *str;
  char LINE[1023];
  
  
  //* connectivity will be assembled
  this->connectivity_allocation();
  
  for (domid=0; domid<NumDomains; domid++) {
    if (filename) {
      char domsuff[8] = "";
      if (Parallel)  sprintf (domsuff,".%ld",domid);
      in = _openFilePNS ("r", "Mesh", path, filename, domsuff);
    }
    else {
      if (Parallel)  errol;
      in = _openFilePNS ("r", "Mesh", NULL, Pd->give_OUT_moFILE()->give_name(), ".T3d.out");
    }
    
    //*  ***  HEAD  ***
    fscanf (in, "%d",          &t3dout);                if (t3dout!=3 && t3dout!=7)     _errorr("error in t3d mesh file");
    fscanf (in, "%ld %ld %ld", &deg, &auxl, &outspec);  if (outspec!=8 && outspec!=264) _errorr("error in t3d mesh file");
    if (Parallel) {  fscanf (in,"%ld",&auxl);       if (auxl!=domid+1)              _errorr("error in t3d mesh file");  }
    
    if (t3dout!=7 && Parallel) errol;
    
    // Nn
    fscanf (in, "%ld", &cND);
    if (domid==0) Pjnts.resizeup(cND);
    else          if (Pjnts() != cND)  _errorr ("error in t3d mesh file");
    
    // Ne
    for (i=0; i<7; i++) {
      if (i<t3dout) fscanf (in, "%ld", el+i);
      else el[i] = 0;
      if (i) el[i] += el[i-1];
    }
    if (domid==0) Elems.resizeup(el[6]);
    else          if (Elems() != el[6])  _errorr ("error in t3d mesh file");
    
    if (Parallel) {
      fscanf (in,"%ld", NumNodes+domid);  cND = NumNodes[domid];
      NumElems[domid] = 0;
      for (i=0;i<7;i++) { fscanf (in,"%ld",el+i);  NumElems[domid] += el[i]; if (i) el[i] += el[i-1]; }
    }
    
    
    //*  ***  NODES  ***
    for (i=0; i<cND; i++) {
      if (!Parallel) { fscanf (in, "%ld",               &id);          lid = id; }
      else           { fscanf (in, "%ld %ld %ld", &lid, &id, &trash);            }
      
      if (id<0) { auxb = true;  id = -id-1;}
      else      { auxb = false; id =  id-1;}
      
      FP_scan_line_skip_emptyORcommented (in, LINE);  str=LINE;
      
      if (Pjnts[id]==NULL) { Pjnts.assign(id, allocate_point (id, NULL, '-'));  Nodes(id)->read_input (str, FFF_T3d);  if (Parallel) Nodes(id)->initialize_domli (auxb,domid,lid-1); }
      else                   Nodes(id)->add_domain           (id,                                      str, FFF_T3d,                                              auxb,domid,lid-1);
      
      //* RAN a HN se naalokuji podle genatu, viz Point :: attributes_allocation
      //* pak ale musim z Node udelat RAN nebo HN
      //* !! v t3d se attributes stejne nenaplnujou, tak bych je nemusel vubec allokovat, pak bych si usetril pretahavani na
      //* zatim to nebude implementovat, az budu do vtk zadelave HN tak tam se bude cist property az v PointData
      //* tak bych to mohl nejak sloucit
      if (Pjnts[id]->give_pointAttribs()->give_classid() == classHNAttribs) {
    PointAttribs* na = Nodes(id)->release_attributes();
    Node *newnode = new HangingNode(Nodes(id), false);
    newnode->assign_attributes(na);
    delete Pjnts[id];
    Pjnts.assign(id, newnode);
      }
      if (Pjnts[id]->give_pointAttribs()->give_classid() == classRANAttribs) {
    PointAttribs* na = Nodes(id)->release_attributes();
    Node *newnode = new RigidArmNode(Nodes(id), false);
    newnode->assign_attributes(na);
    delete Pjnts[id];
    Pjnts.assign(id, newnode);
      }
    }
    
    
    //*  ***  ELEMENTS  ***
    if (deg != 1)  errol; // jinak bude jiny pocet uzlu
    //
    for (i=0; i<el[6]; i++) {
      if      (i < el[0]) {                  nno = 2;  cg = CG_Line;     }
      else if (i < el[1]) {                  nno = 3;  cg = CG_Triangle; }
      else if (i < el[2]) { if (t3dout==3) { nno = 4;  cg = CG_Tetrahedron; } else { nno = 4;  cg = CG_Quadrangle; } }
      else if (i < el[3]) {                                                          nno = 4;  cg = CG_Tetrahedron;  }
      else if (i < el[4]) {                                                          nno = 6;  cg = CG_Wedge;        }
      else if (i < el[5]) {                                                          nno = 5;  cg = CG_Pyramid;      }
      else if (i < el[6]) {                                                          nno = 8;  cg = CG_Hexahedron;   }
      else errol;
      
      FP_scan_line_skip_emptyORcommented (in, LINE);  str=LINE;
      
      SP_scan_number (str, id); id--;
      
      Elems.assign(id, allocate_element(false, cg, id, id, this, nno, false, domid, i) );
      
      FElems(id)->read_input (str, FFF_T3d);
    }
    
    if (EOF != fscanf (in, "%ld", &auxl))  errol;
    
    fclose (in);
  }// end of loop over domains
}


//*  ************************************
//*  ***   READING OF GEOMETRY MESH   ***
//*  ************************************
void Mesh :: read_mesh_equal_to_model (void)
{
  long i, j;
  const Vertex  *masterex;
  const Element *masterel;
  FElement *elem = NULL;
  
  
  //* nodes
  Pjnts.resizeup(Master->give_Npts());
  
  for (i=0; i<Pjnts(); i++) {
    masterex = dynamic_cast<const Vertex *>(Master->give_Pjnt(i));   if (!masterex)  errol;
    
    Pjnts.assign(i, new Node(masterex, 'n'));
    Pjnts[i]->reset_Geom(this);
    Pjnts[i]->reset_property(0, i+1);
    Nodes(i)->set_master_component (i+1, Master, 1);
  }
  
  
  //* elements
  Elems.reallocup(Master->give_Nels());
  
  for (i=0; i<Master->give_Nels(); i++) {
    masterel = Master->give_Elem(i);
    const GelemAttribs *a = dynamic_cast<const GelemAttribs*>(masterel->give_elemAttribs());  if (!a) errol;
    
    if (masterel->give_classid() == classPolyLine  ||  masterel->give_classid() == classLine) {
      if (masterel->give_nno()-2)
    Elems.reallocplus(masterel->give_nno()-2);
      
      for (j=0; j<masterel->give_nno()-1; j++) {
    elem = new Beam (Elems(), Elems(), this, false, 0, Elems());
    ((FElement*)elem)->attributes_allocation(a);
    
    // nodes
    if (elem->give_nno() != 2)  errol;
    elem->set_node(0, masterel->give_point(j  )->give_ID());
    elem->set_node(1, masterel->give_point(j+1)->give_ID());
    
    elem->set_model_prop(i+1, Master, false);
    Elems.add(elem);
      }
    }
    else  errol;
  }
  
  elemprop = Master->give_elemprop();
  nodeprop = Master->give_nodeprop();
}


//*  ************************************
//*  ***   READING OF GEOMETRY MESH   ***
//*  ************************************
void Mesh :: read_mesh_OOFEM (FILE *stream, long nn, long ne)
{
  // variable declaration
  long i;
  const char *str;
  char LINE[1023];
  femFileFormat fff = FFF_OOFEM;
  char WORD[255];
  
  
  //* NODES
  //
  Pjnts.reallocup(nn);
  for (i=0; i<nn; i++) {
    FP_scan_line_skip_emptyORcommented (stream, LINE);  str=LINE;
    SP_scan_word (str, WORD);
    SP_scan_expected_number_exit (str, i+1, "Invalid structure of OOFEM input file");
    
    if      (STRCMP(WORD, "Node"        ) == 0)  Pjnts.add((Point*)new Node         (this, i, NULL));
    else if (STRCMP(WORD, "RigidArmNode") == 0)  Pjnts.add((Point*)new RigidArmNode (this, i, NULL));
    else if (STRCMP(WORD, "HangingNode" ) == 0)  Pjnts.add((Point*)new HangingNode  (this, i, NULL));
    else _errorr2("Unknown OOFEM node type %s", WORD);
    
    Nodes(i)->read_input (str, fff);
  }
  
  //* ELEMENTS
  FiniteElementTypeSet fets;
  //
  Elems.reallocup(ne);
  for (i=0; i<ne; i++) {
    FP_scan_line_skip_emptyORcommented (stream, LINE);  str=LINE;
    SP_scan_word (str, WORD);
    SP_scan_expected_number_exit (str, i+1, "Invalid structure of OOFEM input file");
    //
    
    FETSet_si2set (WORD, 0, SOL_OOFEM, Pd->give_analgroup(), &fets);
    
    Elems.add( allocate_element(false, fets.cg, Elems(), Elems(), this) );
    
    FElems(i)->give_elemAttribs()->set_FETS(&fets);
    FElems(i)->read_input (str, fff);
  }
}

//*  ************************************
//*  ***   READING OF GEOMETRY MESH   ***
//*  ************************************
void Mesh :: read_mesh_SIFEL (FILE *stream)
{
  // variable declaration
  long i, nn, ne;
  const char *str;
  char LINE[1023];
  femFileFormat fff = FFF_SIFEL;
  int auxi;
  
  
  //* NODES
  FP_skip_comment(stream);  fscanf (stream, "%ld", &nn);
  Pjnts.reallocup(nn);
  for (i=0; i<nn; i++) {
    FP_scan_line_skip_emptyORcommented (stream, LINE);  str=LINE;
    SP_skip_word (str, 4); SP_scan_number (str, auxi);  str=LINE;
    SP_scan_expected_number_exit (str, i+1, "Invalid structure of SIFEL input file");
    
    if      (auxi > 0)  Pjnts.add((Point*)new Node         (this, i, NULL));
    else if (auxi < 0)  Pjnts.add((Point*)new HangingNode  (this, i, NULL));
    else _errorr("Unknown SIFEL node");
    
    Nodes(i)->read_input (str, fff);
  }
  
  //* CONSTRAINED NODES
  long id;
  FP_skip_comment(stream);  fscanf (stream, "%ld", &nn);
  for (i=0; i<nn; i++) {
    fscanf (stream, "%ld", &id); id--;
    FP_scan_line (stream, LINE);
    Nodes(id)->give_pointAttribs()->set_dofbc(LINE, fff);
  }
  
  //* ELEMENTS
  FiniteElementTypeSet fets;
  FP_skip_comment(stream);  fscanf (stream, "%ld", &ne);
  Elems.reallocup(ne);
  for (i=0; i<ne; i++) {
    FP_scan_line_skip_emptyORcommented (stream, LINE);  str=LINE;
    //
    SP_scan_expected_number_exit (str, i+1, "Invalid structure of OOFEM input file");
    SP_scan_number               (str, auxi);
    
    FETSet_si2set (NULL, auxi, SOL_SIFEL, Pd->give_analgroup(), &fets);
    
    Elems.add( allocate_element(false, fets.cg, Elems(), Elems(), this) );
    
    FElems(i)->give_elemAttribs()->set_FETS(&fets);
    FElems(i)->read_input (str, fff);
  }
}


//*  ************************************
//*  ***   READING OF GEOMETRY MESH   ***
//*  ************************************
void scan_expect_num (FILE *stream, long expctd)
{
  long aux;
  fscanf (stream, "%ld", &aux);
  if ( aux != expctd)
    _errorr3 ("Invalid structure of the given UNV file: actual != expected -- %ld != %ld", aux, expctd);
}

void Mesh :: read_mesh_UNV (char *path, const char *filename)
{
  // open file
  FILE *in = _openFilePN("r", "Mesh", path, filename);
  
  // variable declaration
  long auxl, nn, te, id, j, domid;
  PoinT coord;
  FElement *elem = NULL;
  //
  domid=0;
  
  
  // GET NUMBER OF LINES
  auxl = FP_number_of_lines (in, true);
  
  Pjnts.reallocup(auxl);
  Elems.reallocup(auxl);
  
  //* nodes
  scan_expect_num (in,   -1);
  scan_expect_num (in, 2411);
  
  id=0;
  while (true) {
    fscanf (in, "%ld", &auxl);
    if (auxl == -1) break;
    if (auxl != id+1)  _errorr ("Invalid structure of the given UNV file");
    
    scan_expect_num (in,  1);
    scan_expect_num (in,  1);
    scan_expect_num (in, 11);
    
    coord.scan_xyz (in);
    Pjnts.add( allocate_point (id, &coord) );
    
    id++;
  }
  
  //* elements
  scan_expect_num (in,   -1);
  scan_expect_num (in, 2412);
  
  id=0;
  while (true) {
    fscanf (in, "%ld", &auxl);
    if (auxl == -1) break;
    //if (auxl != id-1)  _errorr ("Invalid structure of the given UNV file");
    
    fscanf (in,"%ld %ld %ld %ld %ld", &te, &auxl, &auxl, &auxl, &nn);
    switch (te) {
    case 21:  elem = new Beam     (id, id, this, true, domid, id);  fscanf (in,"%ld %ld %ld", &auxl, &auxl, &auxl);  break;
    case 91:  elem = new Triangle (id, id, this, true, domid, id);  break;
    default: _errorr3 ("Unknown %ld-th \"fe descriptor id\": %ld", id+1, te);
    }
    
    // nodes
    if (elem->give_nno() != nn)  errol;
    for (j=0; j<nn; j++) { fscanf (in, "%ld", &auxl);  elem->set_node(j, auxl-1); }
    
    Elems.add(elem);
    
    id++;
  }
  
  if (EOF != fscanf (in, "%ld", &auxl))  errol;
  
  fclose (in);
}



//*  *****************************************************************************************
//*  ***                                READING OF RESULTS                                 ***
//*  *****************************************************************************************
void Mesh :: read_structural_analysis_output (void)
{
  Pd->sodriver()->print_info_message (SODE_start_blue_arrow, "Reading of analysis output file");
  
  const FiLe *file;
  
  //* MAIN OUTPUT
  file = Pd->give_IN_file_results();
  
  if (file == NULL) errol;
  
  switch (file->give_ff()) {
    case FFF_OOFEM:  this->read_output_OOFEM ();  break;
    case FFF_SIFEL:  this->read_output_SIFEL ();  break;
    default:  _errorr2 ("Unsuppoted type of output file from structural analysis package, type %ld", file->give_ff());
  }
  
  //* ADDITIONAL OUTPUT
  file = Pd->give_IN_file_results_addataVTK();
  
  if (file)
    this->read_addata_VTK (file->give_name(), false);
  
  
  Pd->sodriver()->print_info_message_core (SODE_green_ok);
}


long give_nstepsOut (FILE *stream, femFileFormat ff, PAType pt, const Problem *Pd)
{
  char c;
  long answer = 0;
  long pos = 100;
  
  if (ff == FFF_OOFEM) {
    while (answer == 0) {
      if (Pd->solver_file_ptr)  return 1;
      
      //if (Pd->solver_disk_write) {
      //    if (fseek(stream,                       -pos, SEEK_END) != 0)  _errorr ("There is no complete data of one time step in OOFEM output file");
      //}
      //else                        {
      //    return 1;
      //    // fseek nejak nefunguje nebo co
      //    if (fseek(stream, pos, SEEK_CUR) != 0)
      //      _errorr ("There is no complete data of one time step in OOFEM output file");
      //}
      //
      
      if (fseek(stream,                       -pos, SEEK_END) != 0)  _errorr ("There is no complete data of one time step in OOFEM output file");
      
      FP_skip_line(stream);
      
      while (true) {
    c = fgetc(stream);
    if ( c == 'U' &&
         FP_scan_expected_word (stream, "ser") &&
         FP_scan_expected_word (stream, "time") &&
         FP_scan_expected_word (stream, "consumed") &&
         FP_scan_expected_word (stream, "by") &&
         FP_scan_expected_word (stream, "solution") &&
         FP_scan_expected_word (stream, "step"))
      fscanf (stream, "%ld", &answer);
    
    if (c == '\n')  continue;
    if (! FP_skip_line(stream))  break;
      }
      
      pos = 1.9 * pos;
    }
  }
  else if (ff == FFF_SIFEL) {
    switch (pt) {
    case PAT_static_linear:  answer = 1;  break;
    case PAT_stationary:     answer = 1;  break;
    case PAT_nonstationary:  
      while (answer == 0) {
    if (fseek(stream, -pos, SEEK_END) != 0)  _errorr ("There is no complete data of one time step in SIFEL output file");
    
    while (true) {
      // go to new line
      FP_skip_line(stream);
      // go to the line starting with '*'
      if (! FP_skip_to_line_starting_with(stream, '*'))  break;
      
      FP_skip_line(stream);
      if ( FP_scan_expected_word (stream, "Step") &&
           FP_scan_expected_word (stream, "number=") )
        fscanf (stream, "%ld", &answer);
    }
    
    pos = 1.9 * pos;
      }
      answer++;
      break;
    default: _errorr2("Unknown problem type %d", (int)pt);
    }
  }
  return answer;
}

void Mesh :: read_output_OOFEM (void)
{
  if (rslts_solver != SOL_Void)  errol;
  rslts_solver = SOL_OOFEM;
  
  double auxd;
  PAType pt = Pd->give_analloctype();
  
  // stability in oofem runs only for 1d elements
  if (pt == PAT_stability_linear  &&  !Pd->give_primary_geometry()->only_1d_elements())
    pt = PAT_static_linear;
  
  // open file
  const char* name = Pd->give_IN_file_results()->give_name_or_null();
  if (name == NULL)  name = Pd->give_output_file();
  FILE *in;
  //if (Pd->solver_file_ptr)  in = fmemopen(Pd->solver_file_ptr, Pd->solver_file_size, "r");
  //else
  if (name != NULL)  in = _openFilePN  ("r", "Mesh", Pd->give_OUT_Path(), name);
  else               in = _openFilePNS ("r", "Mesh", Pd->give_OUT_Path(), Pd->give_OUT_moFILE()->give_name(), ".oofem.out");
  
  if (Pd->give_nsteps() == 1)  rslts_nsteps = 1;
  else {
    rslts_nsteps = give_nstepsOut (in, FFF_OOFEM, Pd->give_analtype(), Pd);
    rewind(in);
  }
  
  rslts_step = new Lvctr(rslts_nsteps);
  if (pt == PAT_static_nonlinear) rslts_loadlevel = new Dvctr(rslts_nsteps);
  
  
  //* START READING
  FP_skip_line (in, 5);
  FP_scan_expected_word_exit(in, "#    #  #    #  #       #       #    #   OOFEM ver.", "Invalid structure of OOFEM output file", false);
  fscanf (in, "%lf", &auxd);
  ((Problem*)Pd)->set_OOFEM_ver((int)(auxd*10));
  
  FP_skip_line (in, 8);
  
  long i, s;
  for (s=0; s<rslts_nsteps; s++) {
    if (pt == PAT_static_nonlinear) {
      if (! FP_skip_to_line_starting_with (in, "Output for time",        true))  _errorr("Invalid structure of OOFEM output file");
      fscanf (in, "%lf", &auxd);
      if (! FP_skip_expected_string       (in, ", solution step number", true))  _errorr("Invalid structure of OOFEM output file");
      fscanf (in, "%ld", &((*rslts_step)[s]));
      FP_skip_line (in, 1);
      if (! FP_skip_expected_string       (in, "Reached load level :",   true))  _errorr("Invalid structure of OOFEM output file");
      fscanf (in, "%lf", &((*rslts_loadlevel)[s]));
    }
    
    
    // ***   NODES   ***
    if (pt == PAT_stability_linear) { if (! FP_skip_behind_line_starting_with (in, "Linear static:", true))      _errorr("Invalid structure of OOFEM output file"); }
    else                           { if (! FP_skip_behind_line_starting_with (in, "DofManager output:", true))  _errorr("Invalid structure of OOFEM output file"); }
    FP_skip_line (in, 1);
    RTN[RTN_displacement] = true;
    
    for (i=0; i<Pjnts(); i++)  Nodes(i)->read_output_OOFEM (in, s, RTN_displacement);
    
    
    // ***   ELEMENTS   ***
    if (pt != PAT_stability_linear) {
      if (! FP_skip_behind_line_starting_with (in, "Element output:", true))  _errorr("Invalid structure of OOFEM output file");
      FP_skip_line (in, 1);
    }
    RTE[RTE_global_strain] = true;
    RTE[RTE_global_stress] = true;
    
    for (i=0; i<Elems(); i++)  FElems(i)->read_output_OOFEM (in, s);
    
    
    // ***   REACTIONS   ***
    if (! FP_skip_behind_line_starting_with (in, "\tR E A C T I O N S  O U T P U T:", true))  _errorr("Invalid structure of OOFEM output file");
    FP_skip_line (in, 3);
    RTN[RTN_reaction] = true;
    
    for (i=0; i<Pjnts(); i++)  Nodes(i)->read_output_OOFEM (in, s, RTN_reaction);
  }
  
  if (pt == PAT_stability_linear) {
    if (! FP_skip_behind_line_starting_with (in, "Eigen Values are:", true)) _errorr("Invalid structure of OOFEM output file");
    FP_skip_line (in);
    fscanf (in, "%lf", &critical_coeff);
  }
  
  //
  fclose (in);
}

void Mesh :: read_output_SIFEL (void)
{
  if (rslts_solver != SOL_Void)  errol;
  rslts_solver = SOL_SIFEL;
  
  //double auxd;
  //* temporary, see read_output_OOFEM
  bool nonlin = false;
  if (Pd->give_analtype() == PAT_nonstationary)  nonlin = true;
  
  // open file
  const char* name = Pd->give_IN_file_results()->give_name_or_null();
  if (name == NULL)  name = Pd->give_output_file();
  if (name == NULL)  errol;
  FILE *in = _openFilePN("r", "Mesh", Pd->give_OUT_Path(), name);
  
  // nsteps
  rslts_nsteps = give_nstepsOut (in, FFF_SIFEL, Pd->give_analtype(), Pd);
  rslts_step   = new Lvctr(rslts_nsteps);
  if (nonlin) rslts_loadlevel = new Dvctr(rslts_nsteps);
  rewind(in);
  
  // READ
  bool skip = false;
  long i, s, auxl;
  char LINE[1023];
  if (nonlin)  s = -1; // nonlin
  else         s = 0;  // lin
  while (true) {
    if (! FP_skip_to_line_starting_with (in, "**", true))  break;
    
    if (fgetc(in) == '*') {
      FP_skip_line (in);
      if (! FP_scan_expected_word (in, "Step")   )  _errorr("Invalid structure of SIFEL output file");
      if (! FP_scan_expected_word (in, "number="))  _errorr("Invalid structure of SIFEL output file");
      fscanf (in, "%ld", &auxl);
      skip = (s == auxl) ? true : false;
      s = auxl;
      FP_skip_line (in, 2);
      continue;
    }
    if (skip) continue;
    if (s==-1)  _errorr("Invalid structure of SIFEL output file");
    
    FP_scan_line (in, LINE);
    if (LINE[strlen(LINE)-1] != ':')  _errorr("Invalid structure of SIFEL output file");
    
    if      (strcmp(LINE, "Nodal displacements:") == 0) { RTN[RTN_displacement]   = true;  FP_skip_line (in, 1);  for (i=0; i<Pjnts(); i++)   Nodes(i)->read_output_SIFEL (in, s, RTN_displacement);   }
    else if (strcmp(LINE, "Reactions:"          ) == 0) { RTN[RTN_reaction]       = true;  FP_skip_line (in, 1);  for (i=0; i<Pjnts(); i++)   Nodes(i)->read_output_SIFEL (in, s, RTN_reaction);       }
    else if (strcmp(LINE, "Element strains:"    ) == 0) { RTE[RTE_global_strain]  = true;  FP_skip_line (in, 1);  for (i=0; i<Elems(); i++)  FElems(i)->read_output_SIFEL (in, s, RTE_global_strain);         }
    else if (strcmp(LINE, "Element stresses:"   ) == 0) { RTE[RTE_global_stress]  = true;  FP_skip_line (in, 1);  for (i=0; i<Elems(); i++)  FElems(i)->read_output_SIFEL (in, s, RTE_global_stress);         }
    else if (strcmp(LINE, "Nodal unknowns :"    ) == 0) { RTN[RTN_transpunknowns] = true;  FP_skip_line (in, 1);  for (i=0; i<Pjnts(); i++)   Nodes(i)->read_output_SIFEL (in, s, RTN_transpunknowns); }
    else if (strcmp(LINE, "Element gradients :" ) == 0) { RTE[RTE_gradient]       = true;  FP_skip_line (in, 1);  for (i=0; i<Elems(); i++)  FElems(i)->read_output_SIFEL (in, s, RTE_gradient);       }
    else if (strcmp(LINE, "Element fluxes :"    ) == 0) { RTE[RTE_fluxes]         = true;  FP_skip_line (in, 1);  for (i=0; i<Elems(); i++)  FElems(i)->read_output_SIFEL (in, s, RTE_fluxes);         }
    else    _errorr2("Unsuported section in SIFEL output file, \"%s\"", LINE);
  }
  
  //
  fclose (in);
}

// //*  ************************************
// //*  ***   READING OF GEOMETRY MESH   ***
// //*  ************************************
// /**
//    cte VLNA soubory, zatim omezene
// */
// void Mesh :: read_mesh_VLNA (char *path, char *filename)
// {
//   // open file
//   FILE *in = _openFilePN("r", "Mesh", path, filename);
//
//   // variable declaration
//   long i,domid,auxl,nb;
//   double coord[3];
//   //
//   domid=0;
//
//
//   //*  ***  POINTS  ***
//   // head
//   fscanf (in,"%ld",&Nn);   FP_skip_line (in);
//   Pjnts = new Node* [Nn];
//   // coords
//   for (i=0; i<Nn; i++) {
//     FP_scan_array (in, coord, 3);
//     Pjnts[i] = allocate_Point (i, 0, coord);
//   }
//
//   //*  ***  ELEMENTS  ***
//   fscanf (in,"%ld",&Ne);
//   Elems = new FElement* [Ne];
//   for (i=0; i<Ne; i++) {
//     Elems[i] = new Beam (i,this,domid,i);
//     // nodes
//     fscanf (in,"%ld", &auxl);  Elems[i]->set_node(0, auxl);
//     fscanf (in,"%ld", &auxl);  Elems[i]->set_node(1, auxl);
//   }
//
//   //*  ***  BOUNDARY  ***
//   int ndofs = 6;
//   int *bc0 = new int[ndofs];  fill_all_by (bc0, (long)ndofs, 0);
//   int *bc1 = new int[ndofs];  fill_all_by (bc1, (long)ndofs, 1);
//
//   //
//   fscanf (in,"%ld",&nb);
//   //
//   for (i=0; i<nb; i++) {
//     fscanf (in,"%ld",&auxl);
//     Pjnts[auxl]->initialize_BC(ndofs, bc1);
//   }
//
//   for (i=0; i<Nn; i++)
//     if (NULL == Pjnts[i]->give_BC() )
//       Pjnts[i]->initialize_BC(ndofs, bc0);
//
//
//   //*  ***  SNIH - zatizeni  ***
//   fscanf (in,"%ld",&nb);
//   //
//   for (i=0; i<nb; i++) {
//     fscanf (in,"%ld",&auxl);
//     Pjnts[auxl]->set_prop(1);
//   }
//
//
//   //
//   fclose (in);
// }


void Mesh :: find_slaves (void)
{
  long i, j;
  const GeometryComponent *comp;
  
  PoinT nc;
  
  for (i=0; i<Elems(); i++) {
    if (! Elems[i]->give_elemAttribs()->give_HNmstr()) continue;
    
    for (j=0; j<Pjnts(); j++) {
      if (! Elems[i]->is_point_on (Pjnts[j]->give_coords(), comp, &nc) )  continue;
      
      switch (comp->give_classid()) {
      case classNode:  errol;  break;
      case classEdge:
    //this->node_to_hn (j, comp, nc);
    errol;
    
    break;
      default: errol;
      }
      
    }
  }
}


template <class T>
void TF_GPA_cleanup_duplicity (const StdoutDriver *sodriver, T &list, const char *name, char t, bool shake)
{
  sodriver->print_info_message (SODE_start_blue_arrow, "Cleaning of %s duplicities", name);
  
  long i, nd=0;
  for (i=list()-1; i>=0; i--)
    if ( list[i] && list[i]->invisible_duplicated(t) ) {
      list.delete_object(i);
      nd++;
    }
  
  // shake
  if (shake)  TF_GPA_shake_down_reidoid(list);
  
  sodriver->print_info_message (SODE_end, "%d", nd);
}

void Mesh :: cleanup_duplicities (void)
{
  this->connectivity_assembling ();
  
  if (Pd->give_PDBO (PDBO_MultipNode)) {
    this->check_duplicity_nodes();
    TF_GPA_cleanup_duplicity (Pd->sodriver(), Pjnts, "nodes",                 'a', true );
    
    ((Problem*)Pd)->set_PDBO (PDBO_MultipNode, false);
  }
  
  if (Pd->give_PDBO (PDBO_MultipElem)) {
    //* do not change the order of following 
    TF_GPA_cleanup_duplicity (Pd->sodriver(), Elems, "elements - collapsed",  'c', false);
    TF_GPA_cleanup_duplicity (Pd->sodriver(), Elems, "elements - duplicated", 'r', true );
    TF_GPA_cleanup_duplicity (Pd->sodriver(), Faces, "faces",                 'a', true );
    TF_GPA_cleanup_duplicity (Pd->sodriver(), Edges, "edges",                 'a', true );
    
    ((Problem*)Pd)->set_PDBO (PDBO_MultipElem, false);
  }
}


void Mesh :: check_duplicity_nodes (void)
{
  Pd->sodriver()->print_info_message (SODE_start_blue_arrow, "Checking of node duplicities");
  
  //
  long i,j;
  double circum;
  
  long Nn = Pjnts();
  long Ne = Elems();
  
  // compute average circum
  circum = 0.0;
  for (i=0; i<Ne; i++)
    circum += Elems[i]->give_circum();
  circum = circum/Ne;
  
  // ALLOC DP
  DuplicatePoints *dp = new DuplicatePoints(Nn, circum, ZERO_NC);
  
  // set up coords in dp
  for (i=0; i<Nn; i++)
    dp->set_point (i, Pjnts[i]->give_coords() );
  
  // allocate duplicity arrays
  long nd;
  long *Cduplicity, **duplicity;
  
  Cduplicity = new long[Nn/2];
  duplicity = new long*[Nn/2];  for (i=0; i<Nn/2; i++) duplicity[i] = NULL;
  
  dp->assign_cellpoints();
  dp->find_duplicitys(nd, Cduplicity, duplicity);
  
  // save duplicities on nodes
  for (i=0; i<nd; i++) {
    // check the nodes in duplicity squad are not hanging/rigidarm type
    for (j=0; j<Cduplicity[i]; j++)
      if (Pjnts[duplicity[i][j]]->give_classid() != classNode)
    _errorr ("The duplicated node %ld is not pure node, is hanging node or what");
    
    // save of particular nodes
    for (j=1; j<Cduplicity[i]; j++)
      Nodes(duplicity[i][j]) -> setup_duplicity_master( Nodes(duplicity[i][0]) );
    
  }
  
  // DELETE DP
  delete dp;
  deallocateCheck (Cduplicity);
  deallocateCheck ( duplicity, Nn/2, false);
  
  Pd->sodriver()->print_info_message (SODE_end, "found nnd = %d", nd);
}


void Mesh :: find_subdomains (void)
{
  this->connectivity_assembling ();
  
  Pd->sodriver()->print_info_message (SODE_start_blue_arrow, "Subdomain finding");
  
  long Nn = Pjnts();
  long i, nn, level;
  long *cnisd = new long[Nn];
  
  cSubDomains = 0;
  cnisd[cSubDomains] = 0;
  
  level = 0;
  while (true) {
    nn = 0;
    for (i=0; i<Nn; i++) {
      if (level  &&  Pjnts[i]->give_subdom() != -2)  continue;
      
      Nodes(i)->find_parent_subdom (cSubDomains, &nn, &level);
      
      if (nn) {
    cnisd[cSubDomains] += nn;
    nn = 0;
    if (level == 0)
      cnisd[++cSubDomains] = 0;
      }
    }
    //printf ("level %ld;", level);
    if (level) {
      if (--level)  level = 1;
      else          cnisd[++cSubDomains] = 0;
    }

    for (i=0; i<cSubDomains; i++)
      nn += cnisd[i];

    //printf ("level %ld; nn = %ld(%ld); cSubDomains = %ld \n", level, nn, Nn, cSubDomains);

    if (nn == Nn) break;
  }

  cNodesInSD = new long[cSubDomains];
  for (i=0; i<cSubDomains; i++)
    cNodesInSD[i] = cnisd[i];
  
  delete [] cnisd;
  
  Pd->sodriver()->print_info_message (SODE_end, "%ld", cSubDomains);
}
void Mesh :: delete_subdomains_except (long leave)
{
  long i,j;

  for (i=0; i<Pjnts(); i++)
    if ( Pjnts[i]->give_subdom() != leave ) {
      for (j=0; j<Nodes(i)->give_numsuperedge(); j++)  ((Edge*   )Nodes(i)->give_superedge(j))->set_delete_flag(true);
      for (j=0; j<Nodes(i)->give_numsuperface(); j++)  ((Face*   )Nodes(i)->give_superface(j))->set_delete_flag(true);
      for (j=0; j<Nodes(i)->give_numsuperelem(); j++)  ((FElement*)Nodes(i)->give_superelem(j))->set_delete_flag(true);
      Nodes(i)->set_delete_flag(true);
    }
  
  TF_GPA_delete_damned (Pjnts);
  TF_GPA_delete_damned (Edges);
  TF_GPA_delete_damned (Faces);
  TF_GPA_delete_damned (Elems);
}


//void Mesh :: delete_cmfrs (void)
//{
//  Pd->sodriver()->print_info_message (SODE_start_blue_arrow, "Camfours finding");
//
//  this->init_connectivity ();
//
//  long i;
//  long nCmfrs;
//
//  nCmfrs = 0;
//  for (i=0; i<Elems(); i++)
//    // zjisti zda aspon jeden konec prutu se nedotyka pouze jednoho dalsiho
//    if ( Elems[i]->is_cmfr() ) {
//      nCmfrs++;
//    }
//
//  print_info_message ('e',ETC_DEFAULT,"%ld form %ld", nCmfrs, Elems());
//}


//void Mesh :: vlna (void)
//{
//  //
//  if (! Pd->give_PDBO(PDBO_Vlna)) return;
//  //
//  printf ("VLNA\n");
//
//  //
//  long i;
//  long Nn = Pjnts();
//  Lvctr *nodevars = new Lvctr(Nn);
//
//  // *** property == snih ***
//  if (0) {
//    for (i=0; i<Nn; i++)
//      nodevars[0](i) = Pjnts[i]->give_property();
//    //
//    this->print_VTK_PD_scal (false, ".snih", "snih", (const Xvctr**)&nodevars);
//    exit (0);
//  }
//
//  delete nodevars;
//
//  ///*  *** nalezeni subdomen a tisk trojbarevne ***
//  //if (0) {
//  //long *nsups;
//  //  this->find_subdomains();
//  //  // nalezeni podepreni
//  //  allocate    (nsups, cSubDomains);
//  //  fill_all_by (nsups, cSubDomains, (long)0);
//  //  for (i=0; i<Nn; i++)
//  //    if ( Pjnts[i]->give_BC()[0] )
//  //  nsups[Pjnts[i]->give_subdom()] ++;
//  //
//  //  long nonsupSD = 0;
//  //  for (i=0; i<cSubDomains; i++)
//  //    if (nsups[i] == 0)
//  //  nonsupSD++;
//  //
//  //
//  //  // info tisk
//  //  printf ("Pocet               segmentu: %ld\n", cSubDomains);
//  //  printf ("Pocet nepodeprenych segmentu: %ld\n", nonsupSD);
//  //
//  //  //* aux print
//  //  FILE *out = fopen ("aux.txt","w");   openFileTest (out,"txt file");
//  //  for (i=0; i<cSubDomains; i++)  fprintf (out,"%ld\n", cNodesInSD[i]);
//  //  fclose (out);
//  //
//  //  // tisk trojbarevne
//  //  for (i=0; i<Nn; i++)
//  //    if ( Pjnts[i]->give_subdom() == 0)     nodevars[i] = 1;   // klenbova subdomena
//  //    else {
//  //  if (nsups[Pjnts[i]->give_subdom()])  nodevars[i] = 2;   // dalsi SD - podeprene
//  //  else                                 nodevars[i] = 3;   // dalsi SD - nepodeprene
//  //    }
//  //
//  //  //
//  //  this->print_VTK (".trojbarevnej",1, &nodevars);
//  //  exit (0);
//  //}
//
//
//
//
//  // *** nalezeni camfouru ***
//  if (0) {
//    while (true) {
//      this->delete_cmfrs();
//      this->kill_damned_elems();
//
//      this->find_subdomains();
//      if (cSubDomains == 1) break;
//
//      // najde nejvetsi subdomenu
//      long ind, max = 0;
//      for (i=0; i<cSubDomains; i++) if (cNodesInSD[i] > max) { max = cNodesInSD[i]; ind = i; }
//
//      // smaze nejvetsi SD
//      this->delete_subdomains_except (ind);
//      // vynuluje SD
//      for (i=0; i<Nn; i++)
//  Pjnts[i]->set_subdom(-1);
//      cSubDomains = 0;
//      delete [] cNodesInSD;
//    }
//    printf ("\n Nn = %ld\n",Nn);
//    this->print_VTK (false, ".in");
//  }
//
//
//  // *** nalezeni camfouru ***
//  if (Pd->give_PDBO(PDBO_uncfrd)) {
//    this->delete_cmfrs();
//    this->kill_damned_elems();
//
//    // tisk
//    Lvctr *elemvars = new Lvctr(Elems());
//    for (i=0; i<Elems(); i++)
//      if ( Elems[i]->give_delete_flag() )   elemvars[0](i) = 1;
//      else                                  elemvars[0](i) = 0;
//
//    this->print_VTK_CD_scal (false, ".camfour", "camfours", (const Xvctr**)&elemvars);
//    delete elemvars;
//    //exit (0);
//  }
//
//  // *** nalezeni subdomen a smazani mimoklenbovych ***
//  if (1) {
//    this->find_subdomains();
//    
//    // najde nejvetsi subdomenu
//    long ind=0, max = 0;
//    for (i=0; i<cSubDomains; i++)
//      if (cNodesInSD[i] > max) {
//  max = cNodesInSD[i];
//  ind = i;
//      }
//    
//    // smaze nejvetsi SD
//    this->delete_subdomains_except (ind);
//    this->print_VTK (false, ".in");
//  }
//  
//  //
//  printf ("KONEC\n");
//  //exit (0);
//}


void Mesh :: sort_polydata (Xvctr *data) const
{
  long li, pi;
  const_cast<Mesh*>(this)->give_vtk_polydata_counts (pi, li);
  li = 0;
  
  Xvctr *aux = NULL;
  if (dynamic_cast<Lvctr *>(data))  aux = new Lvctr((Lvctr *)data);
  if (dynamic_cast<Dvctr *>(data))  aux = new Dvctr((Dvctr *)data);
  
  for (long i=0; i<Elems(); i++)
    switch (this->give_VTKPDtopology_of_elem(i)) {
    case VTKPD_LINE:     data->cpat(li++, aux, i);  break;
    case VTKPD_POLYGON:  data->cpat(pi++, aux, i);  break;
    default:  _errorr ("unsuported VTK cell type");
    }
  
  delete aux;
}

//*  ***
//*  ***   BC assign   ***
//*  ***
void Geometry :: assign_NodeLoad_by_point_ID    (const BoundaryCond* bc, long id)
{
  Pjnts[id]->give_attributes()->add_load(bc);
}
void Geometry :: assign_NodeLoad_by_point_prop  (const BoundaryCond* bc, long prop)
{
  if (prop < 1)  errol;
  
  for (long i=0; i<Pjnts(); i++)
    if (Pjnts[i]->give_mpropertyORzero() == prop)
      Pjnts[i]->give_attributes()->add_load(bc);
}
void Geometry :: assign_NodeLoad_by_point_coords(const BoundaryCond* bc, const PoinT *coords)
{
  // 0.1 je tam proto aby norm byl mensi nez vzdalenost pro duplicities
  double norm = 0.1 * this->give_aver_elem_circ ();
  
  for (long i=0; i<Pjnts(); i++)
    if (Pjnts[i]->is_identical_to (norm, coords))
      Pjnts[i]->give_attributes()->add_load(bc);
}

void Geometry :: assign_BodyLoad_by_elem_ID     (const BoundaryCond* bc, long id)
{
  Elems[id]->give_attributes()->add_load(bc);
}
void Geometry :: assign_BodyLoad_by_elem_prop   (const BoundaryCond* bc, long prop)
{
  for (long i=0; i<Elems(); i++)
    if (Elems[i]->give_mpropertyORzero() == prop)
      Elems[i]->give_attributes()->add_load(bc);
}

void Geometry :: assign_EdgeLoad_by_elem_ID     (const BoundaryCond* bc, long id, int indx)
{
  if (indx == -1) { if (Elems[id]->give_ned() == 1)  indx = 0;  else _errorr("ConstantEdgeLoad at element with more than 1 edge"); }
  ((Edge*)Elems[id]->give_edge(indx))->give_attributes()->add_load(bc);
}
void Geometry :: assign_EdgeLoad_by_elem_prop   (const BoundaryCond* bc, long prop)
{
  for (long i=0; i<Elems(); i++)
    if (Elems[i]->give_mpropertyORzero() == prop) {
      if (Elems[i]->give_ned() != 1)  _errorr("ConstantEdgeLoad at element with more than 1 edge");
      ((Edge*)Elems[i]->give_edge(0))->give_attributes()->add_load(bc);
    }
}

void Geometry :: assign_FaceLoad_by_elem_ID     (const BoundaryCond* bc, long id, int indx)
{
  if (indx == -1) { if (Elems[id]->give_nfa() == 1)  indx = 0;  else _errorr("ConstantFaceLoad at element with more than 1 face"); }
  ((Face*)Elems[id]->give_face(indx))->give_attributes()->add_load(bc);
}
void Geometry :: assign_FaceLoad_by_elem_prop   (const BoundaryCond* bc, long prop)
{
  for (long i=0; i<Elems(); i++)
    if (Elems[i]->give_mpropertyORzero() == prop) {
      if (Elems[i]->give_nfa() != 1)  _errorr("ConstantFaceLoad at element with more than 1 face");
      ((Face*)Elems[i]->give_face(0))->give_attributes()->add_load(bc);
    }
}
void Geometry :: assign_FaceLoad_by_face_prop   (const BoundaryCond* bc, long prop)
{
  for (long i=0; i<Faces(); i++)
    if (Faces[i]->give_mpropertyORzero() == prop)
      Faces[i]->give_attributes()->add_load(bc);
}


DOFsPerNode Geometry :: give_DOFsPerNode (void)
{
  DOFsPerNode answer, aux;
  answer = Elems[0]->give_elemAttribs()->give_dpn();
  
  for (long i=1; i<Elems() ; i++) {
    aux = Elems[i]->give_elemAttribs()->give_dpn();
    
    if (aux != answer)
      _errorr ("mismatch in DOFsPerNode");
  }
  
  return answer;
}

//void Mesh :: shake_nodes (void)
//{
//  Node* node;
//
//  node = Pjnts[0];
//
//  for (long i=0; i<Nn-1; i++)
//    Pjnts[i] = Pjnts[i+1];
//
//  Pjnts[Nn-1] = node;
//}
//
//void Mesh :: shake_elements (void)
//{
//  FElement* elem;
//
//  elem = Elems[0];
//
//  for (long i=0; i<Ne-1; i++)
//    Elems[i] = Elems[i+1];
//
//  Elems[Ne-1] = elem;
//}





//void Mesh :: print_geom_info (void)
//  long i;
//  long nebe,ntrs,nshl,nhmt,nebr,nrbr,nqbr;
//  ;    nebe=ntrs=nshl=nhmt=nebr=nrbr=nqbr=0;
//
//  for (i=0; i<Ne; i++)
//    switch (Elems[i]->te){
//    case LiBeam:                     //  1
//    case Beam:      nebe++; break;   //  2
//    case Truss:     ntrs++; break;   //  3
//    case Shell:     nshl++; break;   // 31
//    case Quad1hmt:  nhmt++; break;   // 35
//    case  Brick:    nebr++; break;   // 53
//    case RBrick:    nrbr++; break;   // 54
//    case QBrick:    nqbr++; break;   // 55
//    default: _errorr2 ("Unknown element type (element number %d) in function output_print", i+1);
//    }
//
//  for (i=0;i<Np;i++){
//    ntrs += PolyLines[i]->give_nhn() - 1;
//  }
//
//  ;          fprintf (stdout,"\n Number of elems:    %6ld\n",Ne);
//  if (nebe)  fprintf (stdout,  " Number of beams:    %6ld\n",nebe);
//  if (ntrs)  fprintf (stdout,  " Number of trusses:  %6ld\n",ntrs);
//  if (nshl)  fprintf (stdout,  " Number of Shells:   %6ld\n",nshl);
//  if (nhmt)  fprintf (stdout,  " Number of Quad1hmt: %6ld\n",nhmt);
//  if (nebr)  fprintf (stdout,  " Number of bricks:   %6ld\n",nebr);
//  if (nrbr)  fprintf (stdout,  " Number of Rbricks:  %6ld\n",nrbr);
//  if (nqbr)  fprintf (stdout,  " Number of Qbricks:  %6ld\n",nqbr);
//}


//*  ********************************************************************************************
//*  *** *** *** ***                       OUTPUT FUNCTIONS                       *** *** *** ***
//*  ********************************************************************************************

//*  *************************************************
//*  *** ***          BLOCKS PRINTING          *** ***
//*  *************************************************

void Mesh :: print_block_nodes (FILE *stream, femFileFormat fff, int did)
{
  for (long i=0; i<Pjnts(); i++)
    if (Nodes(i)->give_lid_id(did) != -1)
      Nodes(i)->print_row (stream, fff, did);
  
}
void Mesh :: print_block_supported_dofs (FILE *stream, femFileFormat fff, int did)
{
  if (fff == FFF_SIFEL) {
    long n = 0;
    for (long i=0; i<Pjnts(); i++)
      if (Nodes(i)->give_lid_id(did) != -1)
    if (Nodes(i)->give_pointAttribs()->is_supported())
      n++;
    
    fprintf (stream,"\n%ld\n", n);
  }
  
  if (fff == FFF_SIFEL  ||  fff == FFF_ANSYS) {
    for (long i=0; i<Pjnts(); i++)
      if (Nodes(i)->give_lid_id(did) != -1)
    Nodes(i)->give_pointAttribs()->print_supported_BC_to_line(stream, fff);
  }
}
void Mesh :: print_block_elems (FILE *stream, femFileFormat fff, int did)
{
  for (long i=0; i<Elems(); i++)
    if (FElems(i)->give_domain() == did)
      FElems(i)->print_row (stream, fff);
  
}
void Mesh :: print_block_nodal_load (FILE *stream, femFileFormat fff, int did)
{
  if (fff == FFF_SIFEL) {
    long n = 0;
    for (long i=0; i<Pjnts(); i++)
      if (Nodes(i)->give_pointAttribs()->is_loaded())
    n++;
    
    fprintf (stream, "%ld\n", n);
  }
  
  if (fff == FFF_SIFEL  ||  fff == FFF_ANSYS) {
    for (long i=0; i<Pjnts(); i++)
      Nodes(i)->give_pointAttribs()->print_bc(stream, fff);
  }
}
void Mesh :: print_block_gravity_load_ANSYS (FILE *stream) const
{
  const GPA<BoundaryCond>* bcs = Pd->give_BCs();
  
  for (long i=0; i<(*bcs)(); i++) {
    switch ((*bcs)[i]->give_loctype()) {
    case BC_PV:  break;
    case BC_NL:  break;
    case BC_DW:  
      //if ((*bcs)[i]->give_ncomp() != 3)  errol;
      
      //* all elements have to have dead weight load
      long j, k;
      for (j=0; j<Elems(); j++) {
    const GPA<const BoundaryCond>* bbc = FElems(j)->give_elemAttribs()->give_loads();
    for (k=0; k<(*bbc)(); k++)
      if ((*bbc)[k] == (*bcs)[i])
        break;
    
    if (k == (*bbc)()) _errorr2("element %ld has no dead weight load assign", FElems(j)->give_ID()+1);
      }
      
      fprintf (stream, "ACEL");
      for (j=0; j<(*bcs)[i]->give_ncomp(); j++)
    fprintf (stream, ", %.6g", 9.81 * (*bcs)[i]->give_component(j));
      
      fprintf (stream, "\n\n");
      break;
    default: _errorr2("Unknown boundary condition type %d", (*bcs)[i]->give_loctype());
    }
  }
}



template <class T>
void TF_GPA_print_list (FILE *stream, const char *name, const char *name2, T &list, long id)
{
  fprintf (stream,"\n  %s %5ld  %s ", name, id+1, name2);
  for (long i=0; i<list(); i++)
    fprintf (stream,"  %ld", list[i]->give_ID()+1);
}

void Mesh :: print_control (FILE *conout) const
{
  long i;
  
  fprintf (conout,"\n\n\n\n");
  fprintf (conout,"  N    N   OOOO   DDDDD   EEEEE   SSSS  \n");
  fprintf (conout,"  NN   N  O    O  D    D  E      S    S \n");
  fprintf (conout,"  N N  N  O    O  D    D  EEE     SS    \n");
  fprintf (conout,"  N  N N  O    O  D    D  EEE       SS  \n");
  fprintf (conout,"  N   NN  O    O  D    D  E      S    S \n");
  fprintf (conout,"  N    N   OOOO   DDDDD   EEEEE   SSSS  \n");
  
  fprintf (conout,"\n\n EDGES SUPERIOR TO NODE \n "   );  for (i=0;i<Pjnts();i++)  TF_GPA_print_list (conout, "node", "superior edges ",    *give_Node(i)->give_superedges(), i);
  fprintf (conout,"\n\n ELEMENTS SUPERIOR TO NODE \n ");  for (i=0;i<Pjnts();i++)  TF_GPA_print_list (conout, "node", "superior elements ", *give_Node(i)->give_superelems(), i);
  
  
  fprintf (conout,"\n\n\n");
  fprintf (conout,"  EEEEE  DDDDD    GGGG   EEEEE   SSSS  \n");
  fprintf (conout,"  E      D    D  G    G  E      S    S \n");
  fprintf (conout,"  EEE    D    D  G       EEE     SS    \n");
  fprintf (conout,"  EEE    D    D  G  GGG  EEE       SS  \n");
  fprintf (conout,"  E      D    D  G    G  E      S    S \n");
  fprintf (conout,"  EEEEE  DDDDD    GGGG   EEEEE   SSSS  \n");
  
  fprintf (conout,"\n\n NODES ON EDGE \n "            );  for (i=0;i<Edges();i++)  TF_GPA_print_list (conout, "edge", "nodes",              *Edges[i]->give_points(),   i);
  fprintf (conout,"\n\n FACES SUPERIOR TO EDGE \n "   );  for (i=0;i<Edges();i++)  TF_GPA_print_list (conout, "edge", "superior faces ",    *Edges[i]->give_superfaces(), i);
  fprintf (conout,"\n\n ELEMENTS SUPERIOR TO EDGE \n ");  for (i=0;i<Edges();i++)  TF_GPA_print_list (conout, "edge", "superior elements ", *Edges[i]->give_superelems(), i);
  
  
  fprintf (conout,"\n\n\n");
  fprintf (conout,"  FFFFF   AAAA    CCCC   EEEEE   SSSS  \n");
  fprintf (conout,"  F      A    A  C    C  E      S    S \n");
  fprintf (conout,"  FFF    A    A  C       EEE     SS    \n");
  fprintf (conout,"  FFF    AAAAAA  C       EEE       SS  \n");
  fprintf (conout,"  F      A    A  C    C  E      S    S \n");
  fprintf (conout,"  F      A    A   CCCC   EEEEE   SSSS  \n");
  
  fprintf (conout,"\n\n NODES ON FACE \n "            );  for (i=0;i<Faces();i++)  TF_GPA_print_list (conout, "face", "nodes",             *Faces[i]->give_points(),     i);
  fprintf (conout,"\n\n EDGES ON FACE \n "            );  for (i=0;i<Faces();i++)  TF_GPA_print_list (conout, "face", "edges",             *Faces[i]->give_edges(),     i);
  fprintf (conout,"\n\n ELEMENTS SUPERIOR TO FACE \n ");  for (i=0;i<Faces();i++)  TF_GPA_print_list (conout, "face", "superior elements ", *Faces[i]->give_superelems(), i);
  
  
  fprintf (conout,"\n\n\n");
  fprintf (conout,"  EEEEE  L      EEEEE  M    M  EEEEE  N    N  TTTTTT   SSSS  \n");
  fprintf (conout,"  E      L      E      MM  MM  E      NN   N    TT    S    S \n");
  fprintf (conout,"  EEE    L      EEE    M MM M  EEE    N N  N    TT     SS    \n");
  fprintf (conout,"  EEE    L      EEE    M    M  EEE    N  N N    TT       SS  \n");
  fprintf (conout,"  E      L      E      M    M  E      N   NN    TT    S    S \n");
  fprintf (conout,"  EEEEE  LLLLL  EEEEE  M    M  EEEEE  N    N    TT     SSSS  \n");
  
  fprintf (conout,"\n\n NODES ON ELEMENT \n "); for (i=0;i<Elems();i++) TF_GPA_print_list (conout, "element", "nodes", *Elems[i]->give_points(), i);
  fprintf (conout,"\n\n EDGES ON ELEMENT \n "); for (i=0;i<Elems();i++) TF_GPA_print_list (conout, "element", "edges", *Elems[i]->give_edges(), i);
  fprintf (conout,"\n\n FACES ON ELEMENT \n "); for (i=0;i<Elems();i++) TF_GPA_print_list (conout, "element", "faces", *Elems[i]->give_faces(), i);
}



//*  ************************************************
//*  *** ***        PRINTING VTK FILE         *** ***
//*  ************************************************

void Mesh :: resize_Lata (long newsize)
{
  if (Lata.give_size() >= newsize)  return;
  
  long oldsize = Lata.give_size();
  Lata.resize(newsize);
  
  for (; oldsize<newsize; oldsize++)
    Lata.assign(oldsize, new Lvctr);
  
}

void Mesh :: resize_Lata_members (long newsize)
{
  if (Lata[0]->give_size() != newsize) {
    long i = Lata.give_size();
    while (i--)
      Lata[i]->resize_ignore_vals(newsize);
  }
}

void Mesh :: resize_Data (long newsize)
{
  if (Data.give_size() >= newsize)  return;
  
  long oldsize = Data.give_size();
  Data.resize(newsize);
  
  for (; oldsize<newsize; oldsize++)
    Data.assign(oldsize, new Dvctr);
  
}

void Mesh :: resize_Data_members (long newsize)
{
  if (Data[0]->give_size() != newsize) {
    long i = Data.give_size();
    while (i--)
      Data[i]->resize_ignore_vals(newsize);
  }
}



//*  ***   START   ***
void print_VTK_START (Stream *stream, const Mesh *mg, const char *name, VTKdataset ds)
{
  char *file = new char[1023];
  
  if (stream->isFile()) {
    if (mg->give_Pd()->give_OUT_Path()) sprintf (file, "%s%s%s%s", mg->give_Pd()->give_OUT_Path(), mg->give_Pd()->give_OUT_moFILE()->give_name(), name, ".vtk");
    else                                sprintf (file,   "%s%s%s",                                 mg->give_Pd()->give_OUT_moFILE()->give_name(), name, ".vtk");
    stream->open(STRM_file, "w", (const char *&)file);
  }
  else {
    // name
    if (mg->give_Pd()->give_OUT_Path())  sprintf (file, "%s%s%s", mg->give_Pd()->give_OUT_Path(), mg->give_Pd()->give_OUT_moFILE()->give_name(), name);
    else                                 sprintf (file,   "%s%s",                                 mg->give_Pd()->give_OUT_moFILE()->give_name(), name);
    
    if      (ds == VTKDS_polydata)  strcat(file, ".vtp");
    else if (ds == VTKDS_unstruct)  strcat(file, ".vtu");
    else _errorr ("unsuported DataSet");
    
    // xml document
    XMLDocument *doc = new XMLDocument(true, COLLAPSE_WHITESPACE);
    
    char auxs[20];
    if (ds == VTKDS_polydata)  sprintf(auxs, "%s", "PolyData");
    else                       sprintf(auxs, "%s", "UnstructuredGrid");
    
    doc->InsertEndChild( doc->NewDeclaration() );
    
    XMLElement *VTKFile = doc->NewElement("VTKFile");   doc->InsertEndChild(VTKFile);
    VTKFile->SetAttribute("type", auxs);
    VTKFile->SetAttribute("version", "0.1");
    VTKFile->SetAttribute("byte_order", "LittleEndian");
    
    // initiate stream
    stream->open(STRM_tixel, "w", (const char *&)file, VTKFile -> InsertEndChild(doc->NewElement(auxs))->ToElement()->InsertEndChild(doc->NewElement("Piece")) );
  }
  
}

//*  ***    END    ***
void print_VTK_FINISH (Stream *stream)
{
  stream->close();
}

//*  ***   HEAD   ***
void print_VTK_head (FILE *stream, VTKdataset ds)
{
  fprintf (stream,"# vtk DataFile Version 3.0\n");
  fprintf (stream,"VTK file written by MIDAS\n");
  fprintf (stream,"ASCII\n");
  if      (ds == VTKDS_polydata)  fprintf (stream,"DATASET POLYDATA\n");
  else if (ds == VTKDS_unstruct)  fprintf (stream,"DATASET UNSTRUCTURED_GRID\n");
  else _errorr ("unsuported dataset");
}


//*  ***   NODES   ***
void Mesh :: print_VTK_nodes (Stream *stream) const
{
  XMLElement *da = NULL;
  
  if (stream->isFile())
    fprintf (stream->file(), "POINTS %ld float\n", Pjnts());
  else {
    stream->tixel()->ToElement()->SetAttribute("NumberOfPoints", (int)Pjnts());
    //
    da = stream->tix_doc()->NewElement("DataArray");
    da->SetAttribute("type", "Float32");
    da->SetAttribute("NumberOfComponents", "3");
    da->SetAttribute("format", "ascii");
    
    stream->tixel()->InsertEndChild(stream->tix_doc()->NewElement("Points") );
    stream->tixel()->LastChild()->InsertEndChild(da);
  }
  
  char auxs[255];
  for (long i=0; i<Pjnts(); i++) {
    Pjnts[i]->print_row_VTX (auxs);
    
    if (stream->isFile())  fprintf (stream->file(), "%s\n",auxs);
    else                   da->InsertEndChild( stream->tix_doc()->NewText(auxs) );
  }
}
//*  ***   ELEMENTS   ***
void Mesh :: print_VTK_elems (Stream *stream, VTKdataset ds) const
{
  if (ds == VTKDS_polydata) {
    long cl, cp;
    FElems(0)->Msh()->give_vtk_polydata_counts(cl, cp);
    
    if (stream->isFile()) {
      if (cl) {
    fprintf (stream->file(), "LINES %ld %ld\n", cl, cl*3);
    for (long i=0; i<Elems(); i++)
      if (VTKPD_LINE == Elems[i]->give_VTKPDtopology())
        Elems[i]->print_row_VTK (stream->file());
      }
      
      if (cp) {
    long nnumbs = 0;
    for (long i=0; i<Elems(); i++)
      if (VTKPD_POLYGON == Elems[i]->give_VTKPDtopology())
        nnumbs += Elems[i]->give_nno();
    fprintf (stream->file(), "POLYGONS %ld %ld\n", cp, cp + nnumbs);
    for (long i=0; i<Elems(); i++)
      if (VTKPD_POLYGON == Elems[i]->give_VTKPDtopology())
        Elems[i]->print_row_VTK (stream->file());
      }
    }
    else {
      stream->tixel()->SetAttribute("NumberOfVerts",  0 );
      stream->tixel()->SetAttribute("NumberOfLines",  (int)cl);
      stream->tixel()->SetAttribute("NumberOfStrips", 0 );
      stream->tixel()->SetAttribute("NumberOfPolys",  (int)cp);
      
      XMLElement *lines = stream->tix_doc()->NewElement("Lines");   stream->tixel()->InsertEndChild(lines);
      XMLElement *polys = stream->tix_doc()->NewElement("Polys");   stream->tixel()->InsertEndChild(polys);

      XMLElement *dalc = stream->tix_doc()->NewElement("DataArray");  lines->InsertEndChild(dalc);
      dalc->SetAttribute("format", "ascii");
      dalc->SetAttribute("type", "Int32");
      dalc->SetAttribute("Name", "connectivity");

      XMLElement *dalo = (XMLElement*)dalc->ShallowClone(NULL);   lines->InsertEndChild(dalo);
      dalo->SetAttribute("Name", "offsets");

      XMLElement *dapc = (XMLElement*)dalc->ShallowClone(NULL);   polys->InsertEndChild(dapc);
      XMLElement *dapo = (XMLElement*)dalo->ShallowClone(NULL);   polys->InsertEndChild(dapo);
      
      char auxs[255];
      long offl=0, offp=0;
      for (long i=0; i<Elems(); i++) {
    Elems[i]->print_row_VTX (auxs);
    
    switch (Elems[i]->give_VTKPDtopology()) {
    case VTKPD_LINE:    dalc->InsertEndChild(stream->tix_doc()->NewText(auxs));  sprintf (auxs," %ld", offl += Elems[i]->give_nno());  dalo->InsertEndChild(stream->tix_doc()->NewText(auxs)); break;
    case VTKPD_POLYGON: dapc->InsertEndChild(stream->tix_doc()->NewText(auxs));  sprintf (auxs," %ld", offp += Elems[i]->give_nno());  dapo->InsertEndChild(stream->tix_doc()->NewText(auxs)); break;
    default:
      _errorr ("unsuported VTKPDtopology");
    }
      }
    }
  }
  else if (ds == VTKDS_unstruct) {
    if (stream->isFile()) {
      _errorr ("The output to legacy VTK for UNSTRUCTURED dataset is not implemented, use VTU please");
    }
    else {
      stream->tixel()->SetAttribute("NumberOfCells",  (int)Elems());
      
      XMLElement *cells = stream->tix_doc()->NewElement("Cells");   stream->tixel()->InsertEndChild(cells);
      
      XMLElement *dac = stream->tix_doc()->NewElement("DataArray");  cells->InsertEndChild(dac);
      dac->SetAttribute("format", "ascii");
      dac->SetAttribute("type", "Int32");
      dac->SetAttribute("Name", "connectivity");
      
      XMLElement *dao = (XMLElement*)dac->ShallowClone(NULL);   cells->InsertEndChild(dao);
      dao->SetAttribute("Name", "offsets");
      
      XMLElement *dat = (XMLElement*)dac->ShallowClone(NULL);   cells->InsertEndChild(dat);
      dat->SetAttribute("type", "UInt8");
      dat->SetAttribute("Name", "types");
      
      char auxs[255];
      long off=0;
      for (long i=0; i<Elems(); i++) {
    Elems[i]->print_row_VTX (auxs);                                         dac->InsertEndChild(stream->tix_doc()->NewText(auxs));
    sprintf (auxs," %ld", off += Elems[i]->give_nno());                     dao->InsertEndChild(stream->tix_doc()->NewText(auxs));
    sprintf (auxs," %d", CellGeometry_e2i_VTK(Elems[i]->give_cellGeom()));  dat->InsertEndChild(stream->tix_doc()->NewText(auxs));
      }
    }
  }
  else
    _errorr ("unsuported dataset");
}


//*  ***   HEAD + NODES + ELEMENTS   ***
void Mesh :: print_VTK_body (Stream *stream, VTKdataset ds) const
{
  if (stream->isFile())
    print_VTK_head (stream->file(), ds);
  
  print_VTK_nodes (stream);
  print_VTK_elems (stream, ds);
}


//*  ***
//*  ***   DATA - beru long a double a tisknu to jako int a float   ***
//*  ***

//*  ***   HEAD   ***
void print_VTK_init_point_data (Stream *stream, long n) {
  if (stream->isFile())
    fprintf (stream->file(), "POINT_DATA %ld\n", n);
  else {
    stream->tixel()->InsertEndChild( stream->tix_doc()->NewElement("PointData") );
    stream->relink_downL();
  }
}
void print_VTK_init_cell_data (Stream *stream, long n) {
  if (stream->isFile())
    fprintf (stream->file(), "CELL_DATA %ld\n", n);
  else {
    stream->tixel()->InsertEndChild( stream->tix_doc()->NewElement("CellData") );
    stream->relink_downL();
  }
}

//  ***   SCALAR / VECTOR / TENSOR  body   ***
void print_VTK_data_body (Stream *stream, long n, const char *name, const GPA<Array1d> *data, char type, bool r2z=false)
{
  XMLElement *da = NULL;
  
  // check integer x double
  bool I = false;
  if      ((*data)[0]->give_arrayTypedef() == ATlong  ) I = true;
  else if ((*data)[0]->give_arrayTypedef() == ATdouble) I = false;
  else _errorr ("error");
  
  //
  char tn[8];
  if (I) { if (stream->isFile()) sprintf(tn,"int");   else sprintf(tn,"Int32");   }
  else   { if (stream->isFile()) sprintf(tn,"float"); else sprintf(tn,"Float32"); }
  
  long len = (*data)[0]->give_size();
  
  if (stream->isFile()) {
    if      (type == 's') { fprintf (stream->file(), "SCALARS %s %s %ld\n", name, tn, len);  fprintf (stream->file(), "LOOKUP_TABLE default\n");  }
    else if (type == 'v')   fprintf (stream->file(), "VECTORS %s %s\n",    name, tn);
    else if (type == 't')   fprintf (stream->file(), "TENSORS %s %s\n",    name, tn);
    else    _errorr ("error");
  }
  else {
    da = stream->tix_doc()->NewElement("DataArray");  stream->tixel()->InsertEndChild(da);
    da->SetAttribute("format", "ascii");
    da->SetAttribute("type", tn);
    da->SetAttribute("Name", name);
    if (type == 's')  da->SetAttribute("NumberOfComponents", (int)len);
    if (type == 'v')  da->SetAttribute("NumberOfComponents", "3");
    if (type == 't')  da->SetAttribute("NumberOfComponents", "9");
  }
  
  int prec = 6;
  double zero, relzero;
  if (r2z) { zero = 1.e-14;  relzero = 1.e-5; }
  else     { zero = 0.0;     relzero = 1.e-5; }
  
  char *auxs;
  if      (type == 's')  auxs = new char[         (*data)[0] ->length_printed       (prec)];
  else if (type == 'v')  auxs = new char[((Xvctr*)(*data)[0])->length_printed_vector(prec)];
  else                   auxs = new char[((Xvctr*)(*data)[0])->length_printed_tensor(prec)];
  
  for (long i=0; i<n; i++) {
    if      (type == 's')           (*data)[i] ->print           (auxs, prec, zero);
    else if (type == 'v')  ((Xvctr*)(*data)[i])->print_vector    (auxs, prec, zero, true);
    else                   ((Xvctr*)(*data)[i])->print_symtensor (auxs, prec, ((Xvctr*)(*data)[i])->give_zero(zero, relzero), true);
    
    if (stream->isFile())  fprintf (stream->file(), "%s\n",auxs);
    else                   da->InsertEndChild( stream->tix_doc()->NewText(auxs) );
  }
  
  delete [] auxs;
}

//*  ***
//*  ***   FULL VTK FILE   ***
//*  ***
//* VTK file with one point data scalar
//  !!!!!! musim to predelat, aby se tam posilalo pole vektoru [Nn,1] a ne [1,Nn] !!!!!!!!!!
//void Mesh :: print_VTK_PD_scal (bool xml, const char *suff, const char *name, const Xvctr **data)
//{
//  Stream *stream = NULL;
//  VTKdataset ds;
//  
//  if (this->only_vtk_polydata())  ds = VTKDS_polydata;
//  else                            ds = VTKDS_unstruct;
//  
//  if (!xml)  stream = new Stream(STRM_file);
//  else       stream = new Stream(STRM_tixel);
//  
//  print_VTK_START            (stream, this     , suff, ds);
//  print_VTK_body             (stream,                  ds);
//  print_VTK_init_point_data  (stream, this->give_Nn()    );
//  print_VTK_data_body        (stream, this->give_Nn(), name, data, 's');
//  print_VTK_FINISH           (stream);
//  delete stream;
//}

//void Mesh :: print_VTK_CD_scal (bool xml, const char *suff, const char *name, const Xvctr **data) const
//{
//  Stream *stream = NULL;
//  VTKdataset ds;
//  
//  if (this->only_vtk_polydata())  ds = VTKDS_polydata;
//  else                            ds = VTKDS_unstruct;
//  
//  if (!xml)  stream = new Stream(STRM_file);
//  else       stream = new Stream(STRM_tixel);
//  
//  print_VTK_START            (stream, this     , suff, ds);
//  print_VTK_body             (stream                 , ds);
//  print_VTK_init_cell_data   (stream, this->give_Nels()    );
//  
//  // pro polydata: je to treba sortnout aby byly lines a polys za sebou
//  // !!! sort meni/realokuje data ale nevraci spravnej pointer kterej se pak spatne deletuje !!!
//  //
//  if (ds == VTKDS_polydata)  this->sort_polydata ((Xvctr*)data[0]);
//  //if (!xml && ds == VTKDS_polydata)  _warningg("tady je zakomentovanej sort");
//  
//  print_VTK_data_body        (stream, this->give_Nels(), name, data, 's');
//  print_VTK_FINISH           (stream);
//  delete stream;
//}

void Mesh :: print_VTK (femFileFormat format, const char *suff)
{
  Stream *stream = NULL;
  VTKdataset ds = VTKDS_Void;
  
  switch (format) {
  case FFF_VTK:  stream = new Stream(STRM_file);   ds =                            this->giveVTKDS();  break;
  case FFF_VTP:  stream = new Stream(STRM_tixel);  ds = VTKDS_polydata;  if (ds != this->giveVTKDS())  _errorr ("The output to xml VTK for POLYDATA dataset is not supported for this mesh (probably because of brick elements), use VTU please");  break;
  case FFF_VTU:  stream = new Stream(STRM_tixel);  ds = VTKDS_unstruct;                                break;
  default: _errorr("File format for output of results is not specified ");
  }
  
  print_VTK_START            (stream, this, suff, ds);
  print_VTK_body             (stream            , ds);
  
  //* PROPERTY
  if (Pd->give_PDBO(PDBO_OUT_print_property)) {
    this->resize_Lata (Elems());
    this->resize_Lata_members(1);
    
    print_VTK_init_cell_data (stream, Elems());
    
    for (long i=0; i<Elems(); i++)  (*Lata[i])[0] = FElems(i)->give_mpropertyORzero();
    
    print_VTK_data_body (stream, Elems(), "Property", (GPA<Array1d>*)&Lata, 's', true);
  }
  //* PROPERTY
  if (Pd->give_PDBO(PDBO_OUT_print_parent_property)) {
    this->resize_Lata (Elems());
    this->resize_Lata_members(1);
    
    print_VTK_init_cell_data (stream, Elems());
    
    for (long i=0; i<Elems(); i++)  (*Lata[i])[0] = FElems(i)->give_parent_prop();
    
    print_VTK_data_body (stream, Elems(), "Property", (GPA<Array1d>*)&Lata, 's', true);
  }
  
  print_VTK_FINISH           (stream);
  delete stream;
}


void Mesh :: print_JKTK (femFileFormat format, const char *suff)
{
  Pd->sodriver()->print_info_message (SODE_start_blue_arrow, "JKTK file print");
  
  FILE *out = _openFilePNS ("w", "JKTK mesh", this->give_Pd()->give_OUT_Path(), this->give_Pd()->give_OUT_moFILE()->give_name(), ".jktk.top");
  
  
  //*  ***  NODES  ***
  fprintf (out, "%ld\n", Pjnts());
  this->print_block_nodes(out, FFF_JKTK, 0);
  
  //*  ***  ELEMENTS  ***
  fprintf (out, "\n%ld\n", Elems());
  this->print_block_elems(out, FFF_JKTK, 0);
  
  
  fclose (out);
  
  //* file with hanging nodes
  long i, n;
  for (n=i=0; i<Pjnts(); i++)
    if (Pjnts[i]->give_classid() == classHangingNode)
      n++;
  
  if (n) {
    FILE *out = _openFilePNS ("w", "JKTK mesh", this->give_Pd()->give_OUT_Path(), this->give_Pd()->give_OUT_moFILE()->give_name(), ".jktk.hn");
    
    fprintf (out, "%ld\n", n);
    
    for (i=0; i<Pjnts(); i++)
      if (Pjnts[i]->give_classid() == classHangingNode) {
    fprintf (out, "%ld", Pjnts[i]->give_ID()+1);
    Pjnts[i]->give_pointAttribs()->print_row(out, FFF_JKTK, 0);
    fprintf (out, "\n");
      }
    
    fclose (out);
  }
  
  Pd->sodriver()->print_info_message_core (SODE_green_ok);
}


void print_dat_file (const char *file, long n, const Dvctr *data1, const Dvctr *data2, bool zeroline)
{
  if (data1 == NULL  &&  data2 == NULL)  errol;
  
  FILE *in = _openFilePN("w", "Dat file", NULL, file);
  
  //Xvctr *aux = NULL;
  //if (dynamic_cast<Lvctr *>(data))  aux = new Lvctr((Lvctr *)data);
  //if (dynamic_cast<Dvctr *>(data))  aux = new Dvctr((Dvctr *)data);
  
  if (data1 && data1->give_size() < n)  errol;
  if (data2 && data2->give_size() < n)  errol;
  
  if (zeroline)  fprintf (in, "%9.6e %9.6e\n", 0.0, 0.0);
  
  long i;
  if      (data1 == NULL)  for (i=0; i<n; i++)  fprintf (in, "%9.6e %9.6e\n", (double)i+1, (*data2)[i]);
  else if (data2 == NULL)  for (i=0; i<n; i++)  fprintf (in, "%9.6e %9.6e\n", (*data1)[i], (double)i+1);
  else                     for (i=0; i<n; i++)  fprintf (in, "%9.6e %9.6e\n", (*data1)[i], (*data2)[i]);
  
  fclose (in);
}

//*  ************************************************
//*  *** ***             RESULTS              *** ***
//*  ************************************************
void Mesh :: print_results (void)
{
  Pd->sodriver()->print_info_message (SODE_start_blue_arrow, "Results printing");
  
  //* check out the geometry is really preserved
  if (Pd->give_PDBO(PDBO_preserveGeom)) {
    if (origPjntsNum != Pjnts())  errol;
    if (origElemsNum != Elems())  errol;
    for (long i=0; i<Pjnts(); i++)
      if (Pjnts[i]->give_ID() != Pjnts[i]->give_origID()) errol;
  }
  
  //* output format
  bool       xml = true;
  VTKdataset ds  = VTKDS_Void;
  
  switch (Pd->give_OUT_moFILE()->give_ff()) {
  case FFF_VTK: xml = false; ds =                            this->giveVTKDS();  break;
  case FFF_VTP: xml = true;  ds = VTKDS_polydata;  if (ds != this->giveVTKDS())  _errorr ("The output to xml VTK for POLYDATA dataset is not supported for this mesh (probably because of brick elements), use VTU please");  break;
  case FFF_VTU: xml = true;  ds = VTKDS_unstruct;                                break;
  default: _errorr("File format for output of results is not specified ");
  }
  
  Dvctr steps(rslts_nsteps);
  long numsteps = 0;
  //int globnDOFs = Pd->give_global_nDOFs();
  
  //* variables
  long i, step;
  Stream *stream;
  char suff[12];
  SStype SST;
  
  long Nn = Pjnts();
  long Ne = Elems();
  long NnNe = (Nn>Ne ? Nn : Ne);
  long dstep = Pd->give_OUT_printStep();
  
  //* ALLOCATION
  this->resize_Data (NnNe);
  this->resize_Lata (NnNe);
  pData.free_ptrs();  pData.resize(NnNe);
  
  //* type of stress/strain results at elements
  //* docasne vyhodit node priznaky varianty !!!!
  //* S - sigma, T - tau
  //*                       <--        OOFEM        -->         Pozn.
  //* SST_truss         ... Loc   N  0  0  0 0 0   [N??]
  //* SST_beam
  //* SST_3d
  //* SST_3dshell
  //* SST_3d_nodes
  //* SST_3dshell_nodes
  //* SST_2dstrain      ... Glob  Sx Sy 0  0 0 Txy [Pa]     OOFEM   dava stress_z, ale ja ho nezobrazuju
  //* SST_2dstress      ... Glob  Sx Sy 0  0 0 Txy [Pa]     OOFEM nedava strain_z
  //* SST_transp2d      ...
  //* 
  int sst[cSSTypes];   memset(sst, 0, cSSTypes*sizeof(int));
  for (i=0; i<Ne; i++)
    sst[FElems(i)->give_elemAttribs()->give_sst()]++;
  
  //* plasticity
  bool plasticityMaterial = false;
  for (i=0; i<Ne; i++)
    if (FElems(i)->give_elemAttribs()->give_mat() && FElems(i)->give_elemAttribs()->give_mat()->isOOFEMplast()) {
      plasticityMaterial = true;
      break;
    }
  
  Dvctr *yielding=NULL;
  if (plasticityMaterial)  yielding = new Dvctr(rslts_nsteps);
  
  //* LOOP OVER STEPS
  //* na elementech ss prumeruju na jednu hodnotu
  for (step=0; step<rslts_nsteps; step++) {
    // skip steps
    if (step!=1000 && step!=rslts_nsteps-1) {
      // skip non yielding steps
      //if (plasticityMaterial) {
      //    int sum = 0;
      //    for (i=0; i<Ne; i++)  sum += ( FElems(i)->give_results(step, RTE_plaststrain) != NULL ) ? 1 : 0;
      //    if (sum == 0 || sum == Ne)   continue;
      //}
      // skip steps
      if (step%dstep!=0)  continue;
    }
    
    steps[numsteps] = step+1;
    
    if (xml)  stream = new Stream(STRM_tixel);
    else      stream = new Stream(STRM_file);
    
    if (rslts_nsteps == 1)   sprintf (suff, ".rslts");
    else                     sprintf (suff, ".rslts.%03ld", step+1);
    
    //*  ***   BODY   ***
    print_VTK_START (stream, this, suff, ds);
    print_VTK_body  (stream            , ds);
    
    
    //*  **********************
    //*  ***   POINT DATA   ***
    //*  **********************
    print_VTK_init_point_data (stream, Nn);
    
    
    //* DISPLACEMENT
    if (RTN[RTN_displacement]) {   for (i=0; i<Nn; i++)  Nodes(i)->setup_full_alloc_DOFvals_at(Data[i], RTN_displacement, step);        print_VTK_data_body (stream, Nn, "uknw_displacement", (GPA<Array1d>*)&Data, 'v');
                                   for (i=0; i<Nn; i++)  Data[i]->shift(3);                                                             print_VTK_data_body (stream, Nn, "uknw_rotation",     (GPA<Array1d>*)&Data, 'v');
                   for (i=0; i<Nn; i++)  Data[i]->deshift();                                                                                                                                               }
    
    //* REACTIONS
    if (RTN[RTN_reaction]) {       for (i=0; i<Nn; i++)  Nodes(i)->setup_full_alloc_DOFvals_at(Data[i], RTN_reaction, step);            print_VTK_data_body (stream, Nn, "reactions_forces",  (GPA<Array1d>*)&Data,  'v');
                                   for (i=0; i<Nn; i++)  Data[i]->shift(3);                                                             print_VTK_data_body (stream, Nn, "reactions_moments", (GPA<Array1d>*)&Data,  'v');
                   for (i=0; i<Nn; i++)  Data[i]->deshift();                                                                                                                                               }
    
    //* TRANSPORT UNKNOWNS
    if (RTN[RTN_transpunknowns]) { for (i=0; i<Nn; i++)  pData.assign(i, Nodes(i)->give_resultsN_dv(step, RTN_transpunknowns));  print_VTK_data_body (stream, Nn, "uknw_transp_media", (GPA<Array1d>*)&pData, 's'); }
    
    this->resize_Data_members(1);
    
    //* ADAPTIVITY, PERCENTAGE ERROR AT ELEMENT
    if (RTN[RTN_ada_h_act]) {
      for (i=0; i<Nn; i++)
    (*Data[i])[0] = Nodes(i)->give_resultsN_ds(step, RTN_ada_h_act)[0]();
      print_VTK_data_body (stream, Nn, "ada_h_act",         (GPA<Array1d>*)&Data,  's');
    }
    if (RTN[RTN_ada_h_new]) {      for (i=0; i<Nn; i++)  (*Data[i])[0] = Nodes(i)->give_resultsN_ds(step, RTN_ada_h_new)[0]();           print_VTK_data_body (stream, Nn, "ada_h_new",         (GPA<Array1d>*)&Data,  's'); }
    
    //* ADAPTIVITY, REFINED DERIVATIVES
    //if (RTN[RTN_ada_refder]) {
    //  for (i=0; i<Nn; i++)  datap[i] = (Dvctr*)Nodes(i)->give_resultsN_d(0   , RTN_ada_refder);      print_VTK_data_body (stream, Nn, "othr_ada_refder", (GPA<Array1d>*)&pData, 's');
    //}
    if (0) {
      char label[24];
      
      // normala pro vsechny regiony
      if (0) {
     for (long regid=0; regid<this->nRegions; regid++) {
       sprintf (label, "ada_normal_reg%02ld", regid);  for (i=0; i<Nn; i++)  pData.assign(i, Nodes(i)->give_resultsE_dv(regid, step, RTE_ada_normal));          print_VTK_data_body (stream, Nn, label,        (GPA<Array1d>*)&pData, 'v');
     }
      }
      
      // ada otocena rovina
      if (0) {
    this->resize_Data_members(3);
    
    // pouze jeden region, ten posledni
    long regid = this->nRegions-1;
    sprintf (label, "ada_osax_reg%02ld", regid);  for (i=0; i<Nn; i++)  Data[i]->beCopyOf( Nodes(i)->give_resultsE_dm(regid, step, RTE_ada_normal)->give_ptr2val(0,0));  print_VTK_data_body (stream, Nn, label,  (GPA<Array1d>*)&Data, 'v');
    sprintf (label, "ada_osay_reg%02ld", regid);  for (i=0; i<Nn; i++)  Data[i]->beCopyOf( Nodes(i)->give_resultsE_dm(regid, step, RTE_ada_normal)->give_ptr2val(1,0));  print_VTK_data_body (stream, Nn, label,  (GPA<Array1d>*)&Data, 'v');
    sprintf (label, "ada_osaz_reg%02ld", regid);  for (i=0; i<Nn; i++)  Data[i]->beCopyOf( Nodes(i)->give_resultsE_dm(regid, step, RTE_ada_normal)->give_ptr2val(2,0));  print_VTK_data_body (stream, Nn, label,  (GPA<Array1d>*)&Data, 'v');
      }
    }
    if (Pd->give_PDBO(PDBO_P_refder)) {
      long regid;
      const Dvctr* auxD;
      char label[24];
      
      this->resize_Data_members(3);
      
      for (regid=0; regid<this->nRegions; regid++) {
    
    for (i=0; i<Nn; i++) {
      auxD = Nodes(i)->give_resultsE_dv(regid, step, RTE_global_stress);
      Data[i][0][0] = Data[i][0][1] = 0.0;
      if (auxD == NULL) Data[i][0][2] = 0.0;
      else              Data[i][0][2] = auxD[0][6]; // sigma x
    }
    
    sprintf (label, "othr_ada_refder_%02ld", regid);
    print_VTK_data_body (stream, Nn, label, (GPA<Array1d>*)&Data, 'v');
      }
    }
    
    
    //* VARIABLES - D x R  strain x stress
    //if (RTE[RTE_strain]  &&  sst[SST = SST_3d]) { for (i=0; i<Nn; i++)  Nodes(i)->give_ssstate (data[i], SST, RVTstrn_displ, 'g', step);  print_VTK_data_body (stream, Nn, "Dstrain", (const Array1d**)data, 't', true); }
    //if (RTE[RTE_stress]  &&  sst[SST = SST_3d]) { for (i=0; i<Nn; i++)  Nodes(i)->give_ssstate (data[i], SST, RVTstrs_displ, 'g', step);  print_VTK_data_body (stream, Nn, "Dstress", (const Array1d**)data, 't', true); }
    
    //if (RTE[RTE_strain]  &&  sst[SST = SST_3dshell]) { for (i=0; i<Nn; i++)  Nodes(i)->give_ssstate (data[i], SST, RVTstrn_displ, 'g', step);  print_VTK_data_body (stream, Nn, "Dstrain", (const Array1d**)data, 't', true);
    //                                                       for (i=0; i<Nn; i++)  Nodes(i)->give_ssstate (data[i], SST, RVTstrn_rotat, 'g', step);  print_VTK_data_body (stream, Nn, "Rstrain", (const Array1d**)data, 't', true); }
    //if (RTE[RTE_stress]  &&  sst[SST = SST_3dshell]) { for (i=0; i<Nn; i++)  Nodes(i)->give_ssstate (data[i], SST, RVTstrs_displ, 'g', step);  print_VTK_data_body (stream, Nn, "Dstress", (const Array1d**)data, 't', true);
    //                                                       for (i=0; i<Nn; i++)  Nodes(i)->give_ssstate (data[i], SST, RVTstrs_rotat, 'g', step);  print_VTK_data_body (stream, Nn, "Rstress", (const Array1d**)data, 't', true); }
    
    //if (sst[SST_beam])                            {  for (i=0; i<Nn; i++)  Nodes(i)->give_ssstate (data[i], RVT_endDispl, step); print_VTK_data_body (stream, Nn, "endDispl", (const Array1d**)data, 't', true);  }
    //if (sst[SST_beam])                            {  for (i=0; i<Nn; i++)  Nodes(i)->give_ssstate (data[i], RVT_endForce, step); print_VTK_data_body (stream, Nn, "endForce", (const Array1d**)data, 't', true);  }
    
    //* 2dstress - 3 nenulove slozky ve vektoru v globalnim SS, predpokladam rovinu xy
    //* do uzlu budu hazet podle nejblizsiho ip, ale zatim to implementovat nebudu
    
    if (xml)  stream->relink_up();
    
    
    //*  *********************
    //*  ***   CELL DATA   ***
    //*  *********************
    //* head
    print_VTK_init_cell_data (stream, Ne);
    
    //*
    //* ***   NON-RESULS VALUES   ***
    //*
    
    //* MODELA PARENT ELEMENT ID and PROPERTY
    this->resize_Lata_members(1);
    if (Pd->give_PDBO(PDBO_OUT_print_model_parent)) {        for (i=0; i<Ne; i++)  (*Lata[i])[0] = FElems(i)->give_parent_id();                                 print_VTK_data_body (stream, Ne, "ID_model_parent",      (GPA<Array1d>*)&Lata, 's', true); }
    if (Pd->give_PDBO(PDBO_OUT_print_property)) {            for (i=0; i<Ne; i++)  (*Lata[i])[0] = FElems(i)->give_parent_prop();                               print_VTK_data_body (stream, Ne, "Property",             (GPA<Array1d>*)&Lata, 's', true); }
    
    
    //*
    //* ***   STRAIN & STRESS   ***
    //*
    
    //* 1d ELEMENTS
    
    //* TRUSS ELEMENTS - strain x force
    if (RTE[RTE_global_strain] && sst[SST = SST_truss]) {    for (i=0; i<Ne; i++)  FElems(i)->give_ssstate (Data[i], SST, RVTstrn, 'g', step);                  print_VTK_data_body (stream, Ne, "strain_truss",         (GPA<Array1d>*)&Data, 's', true); }
    if (RTE[RTE_global_stress] && sst[SST = SST_truss]) {    for (i=0; i<Ne; i++)  FElems(i)->give_ssstate (Data[i], SST, RVTstrs, 'g', step);                  print_VTK_data_body (stream, Ne, "stress_truss",         (GPA<Array1d>*)&Data, 's', true); }
    
    //* BEAM ELEMENTS
    if (RTE[RTE_global_strain] && sst[SST = SST_beam]) {     for (i=0; i<Ne; i++)  FElems(i)->give_ssstate (Data[i], SST, RVTstrn, 'g', step);                  print_VTK_data_body (stream, Ne, "strain_beam3D_NVM",    (GPA<Array1d>*)&Data, 'v', true);
                                                             for (i=0; i<Ne; i++)  Data[i]->shift(3);                                                           print_VTK_data_body (stream, Ne, "strain_beam3D_NVyVz",  (GPA<Array1d>*)&Data, 'v', true);
                                 for (i=0; i<Ne; i++)  Data[i]->shift(3);                                                           print_VTK_data_body (stream, Ne, "strain_beam3D_MxMyMz", (GPA<Array1d>*)&Data, 'v', true);
                                 for (i=0; i<Ne; i++)  Data[i]->shift(-6);                                                                                                                                                     }
    if (RTE[RTE_global_stress] && sst[SST = SST_beam]) {     for (i=0; i<Ne; i++)  FElems(i)->give_ssstate (Data[i], SST, RVTstrs, 'g', step);                  print_VTK_data_body (stream, Ne, "stress_beam3D_NVM",    (GPA<Array1d>*)&Data, 'v', true);
                                                         for (i=0; i<Ne; i++)  Data[i]->shift(3);                                                           print_VTK_data_body (stream, Ne, "stress_beam3D_NVyVz",  (GPA<Array1d>*)&Data, 'v', true);
                                 for (i=0; i<Ne; i++)  Data[i]->shift(3);                                                           print_VTK_data_body (stream, Ne, "stress_beam3D_MxMyMz", (GPA<Array1d>*)&Data, 'v', true);
                                 for (i=0; i<Ne; i++)  Data[i]->shift(-6);                                                                                                                                                     }
    
    //* 2d ELEMENTS
    
    //* planestrain - 3 nenulove slozky ve vektoru v globalnim SS, predpokladam rovinu xy
    if (RTE[RTE_global_strain] && sst[SST = SST_plstrain]) { for (i=0; i<Ne; i++)  FElems(i)->give_ssstate (Data[i], SST, RVTstrn_displ, 'g', step);            print_VTK_data_body (stream, Ne, "strain_planestrain",   (GPA<Array1d>*)&Data, 'v', true); }
    if (RTE[RTE_global_stress] && sst[SST = SST_plstrain]) { for (i=0; i<Ne; i++)  FElems(i)->give_ssstate (Data[i], SST, RVTstrs_displ, 'g', step);            print_VTK_data_body (stream, Ne, "stress_planestrain",   (GPA<Array1d>*)&Data, 'v', true); }
    
    //* planestress - 3 nenulove slozky ve vektoru v globalnim SS, predpokladam rovinu xy
    if (RTE[RTE_global_strain] && sst[SST = SST_plstress]) { for (i=0; i<Ne; i++)  FElems(i)->give_ssstate (Data[i], SST, RVTstrn_displ, 'g', step);            print_VTK_data_body (stream, Ne, "strain_planestress",   (GPA<Array1d>*)&Data, 'v', true); }
    if (RTE[RTE_global_stress] && sst[SST = SST_plstress]) { for (i=0; i<Ne; i++)  FElems(i)->give_ssstate (Data[i], SST, RVTstrs_displ, 'g', step);            print_VTK_data_body (stream, Ne, "stress_planestress",   (GPA<Array1d>*)&Data, 'v', true); }
    
    //* shell stra
    if (RTE[RTE_global_strain] && sst[SST = SST_3dshell])  { for (i=0; i<Ne; i++)  FElems(i)->give_ssstate (Data[i], SST, RVTstrn_displ, 'g', step);            print_VTK_data_body (stream, Ne, "strainD_shell_glob",   (GPA<Array1d>*)&Data, 't', true);
                                                             for (i=0; i<Ne; i++)  FElems(i)->give_ssstate (Data[i], SST, RVTstrn_rotat, 'g', step);            print_VTK_data_body (stream, Ne, "strainR_shell_glob",   (GPA<Array1d>*)&Data, 't', true); }
    if (RTE[RTE_global_stress] && sst[SST = SST_3dshell])  { for (i=0; i<Ne; i++)  FElems(i)->give_ssstate (Data[i], SST, RVTstrs_displ, 'g', step);            print_VTK_data_body (stream, Ne, "stressD_shell_glob",   (GPA<Array1d>*)&Data, 't', true);
                                                             for (i=0; i<Ne; i++)  FElems(i)->give_ssstate (Data[i], SST, RVTstrs_rotat, 'g', step);            print_VTK_data_body (stream, Ne, "stressR_shell_glob",   (GPA<Array1d>*)&Data, 't', true); }
    
    //* 3d ELEMENTS
    
    //*
    if (RTE[RTE_global_strain] && sst[SST = SST_3d]) {       for (i=0; i<Ne; i++)  FElems(i)->give_ssstate (Data[i], SST, RVTstrn_displ, 'g', step);            print_VTK_data_body (stream, Ne, "strain_3D",            (GPA<Array1d>*)&Data, 't', true); }
    if (RTE[RTE_global_stress] && sst[SST = SST_3d]) {       for (i=0; i<Ne; i++)  FElems(i)->give_ssstate (Data[i], SST, RVTstrs_displ, 'g', step);            print_VTK_data_body (stream, Ne, "stress_3D",            (GPA<Array1d>*)&Data, 't', true); }
     
    //* SECONDARY values
    
    //* plasticityMaterial
    if (plasticityMaterial) {                                for (i=0; i<Ne; i++) (*Lata[i])[0] = ( FElems(i)->give_results(step, RTE_plaststrain) != NULL );   print_VTK_data_body (stream, Ne, "yielding",             (GPA<Array1d>*)&Lata, 's');
      
      for (i=0; i<Ne; i++)  if ( (*Lata[i])[0] )  break;
      if (i==Ne)  (*yielding)[numsteps] = 0;
      else        (*yielding)[numsteps] = 1;
    }
    
    //* CROSS SECTION usage/leverage/utilization/level of usage
    if (Pd->give_PDBO(PDBO_OUT_print_CSusage)  &&  this->give_rslts_nsteps()) {      this->resize_Data_members(1);
                                                             for (i=0; i<Ne; i++)  (*Data[i])[0] = FElems(i)->give_CSusage_elast();                             print_VTK_data_body (stream, Ne, "CSusage_elast",        (GPA<Array1d>*)&Data, 's', true);
                                     for (i=0; i<Ne; i++)  (*Data[i])[0] = FElems(i)->give_CSusage_elast_rel();                         print_VTK_data_body (stream, Ne, "CSusage_elast_rel",    (GPA<Array1d>*)&Data, 's', true);
                                     for (i=0; i<Ne; i++)  (*Lata[i])[0] = FElems(i)->give_CSusage_elast_bool();                        print_VTK_data_body (stream, Ne, "CSusage_elast_bool",   (GPA<Array1d>*)&Lata, 's');       }
    
    //* ADAPTIVITY, PERCENTAGE ERROR AT ELEMENT
    this->resize_Data_members(1);
    this->resize_Lata_members(1);
    //
    if (RTE[RTE_ada_error_pct_l]) {                          for (i=0; i<Ne; i++)  (*Data[i])[0] = FElems(i)->give_results_ds(step, RTE_ada_error_pct_l)[0]();  print_VTK_data_body (stream, Ne, "ada_error_pct_l",      (GPA<Array1d>*)&Data, 's');       }
    if (RTE[RTE_ada_error_pct_g]) {                          for (i=0; i<Ne; i++)  (*Data[i])[0] = FElems(i)->give_results_ds(step, RTE_ada_error_pct_g)[0]();  print_VTK_data_body (stream, Ne, "ada_error_pct_g",      (GPA<Array1d>*)&Data, 's');       }
    if (RTE[RTE_ada_h_act]) {                                for (i=0; i<Ne; i++)  (*Data[i])[0] = FElems(i)->give_results_ds(step, RTE_ada_h_act)[0]();        print_VTK_data_body (stream, Ne, "ada_h_act",            (GPA<Array1d>*)&Data, 's');       }
    if (RTE[RTE_ada_h_new]) {                                for (i=0; i<Ne; i++)  (*Data[i])[0] = FElems(i)->give_results_ds(step, RTE_ada_h_new)[0]();        print_VTK_data_body (stream, Ne, "ada_h_new",            (GPA<Array1d>*)&Data, 's');       }
    //if (RTE[RTE_ada_eia])         {                          for (i=0; i<Ne; i++)  (*Data[i])[0] = FElems(i)->give_results_ds(step, RTE_ada_eia)[0]();          print_VTK_data_body (stream, Ne, "ada_error_aver",       (GPA<Array1d>*)&Data, 's');       }
    //if (RTE[RTE_ada_normal]) {                             for (i=0; i<Ne; i++)   pData.assign(i,FElems(i)->give_results_dv(step, RTE_ada_normal));           print_VTK_data_body (stream, Ne, "ada_normal",           (GPA<Array1d>*)&pData, 'v');      }
    if (this->nRegions) {                                    for (i=0; i<Ne; i++)  (*Lata[i])[0] = FElems(i)->give_regid() + 1;                                 print_VTK_data_body (stream, Ne, "ada_regions",          (GPA<Array1d>*)&Lata, 's');       }
    
    
    if (xml)  stream->relink_up(); // end of CELL DATA
    if (xml)  stream->relink_up(); // end of Piece
    if (xml)  stream->relink_up(); // end of UnstructuredGrid
    
    
    //*  *************************
    //*  ***   APPENDED DATA   ***
    //*  *************************
    //* 
    if (xml) {
      char str[1023];
      XMLElement *txl;
      
      txl = stream->tix_doc()->NewElement("AppendedData");    stream->tixel()->InsertEndChild(txl);
      txl->InsertEndChild( stream->tix_doc()->NewText("_") );
      stream->relink_downL();
      
      if (1)                                              { txl = stream->tix_doc()->NewElement("Volume");                         stream->tixel()->InsertEndChild(txl);  sprintf (str, " %lf ", this->volume_core());             txl->InsertEndChild( stream->tix_doc()->NewText(str) ); }
      if (Pd->give_analloctype() == PAT_stability_linear) { txl = stream->tix_doc()->NewElement("Stability_critical_coefficient"); stream->tixel()->InsertEndChild(txl);  sprintf (str, " %lf ", this->critical_coeff);            txl->InsertEndChild( stream->tix_doc()->NewText(str) ); }
      if (Pd->give_adpa())                                { txl = stream->tix_doc()->NewElement("Error");                          stream->tixel()->InsertEndChild(txl);  sprintf (str, " %.2lf ", Pd->give_adpa()->give_error()); txl->InsertEndChild( stream->tix_doc()->NewText(str) ); }
    }
    
    //*  ***   END   ***
    print_VTK_FINISH (stream);
    delete stream;
    numsteps++;
  }
  
  if (rslts_nsteps > 1) {
    if (plasticityMaterial) {
      print_dat_file ("this.dat.yielding.dat", numsteps, &steps, yielding, true);
      delete yielding;
    }
    
    //const char* anal = Pd->give_OOanalysis();
    //if (anal && SP_scan_expected_word (anal, "NonLinearStatic", CASE)) {
    long idnode = 1;
    int index = 0;
    const Node* node = this->give_Node(idnode);
    
    Dvctr displ(rslts_nsteps);
    
    node->give_displcmnt_in_all_steps (&displ, index);
    
    //if (! rslts_loadlevel)errol;
    
    print_dat_file ("this.dat", rslts_nsteps, &displ, rslts_loadlevel, true);
  }
  
  
  Pd->sodriver()->print_info_message_core (SODE_green_ok);
}

//void Mesh :: print_cabletensileonly (bool mft)
//{
//  if (this->give_anyBeamElem() == 0)  return;
//  
//  Lvctr *data = new Lvctr(Elems());
//  long i, step;
//  char str[16];
//  const Beam *beam;
//  
//  for (step=0; step<nstepsOut; step++) {
//    data->zero();
//    for (i=0; i<Elems(); i++) {
//      beam = dynamic_cast<const Beam *>(Elems[i]);
//      if (beam  &&  0.0 > beam->give_normal_force(step))  data[0](i) = 1;
//    }
//    
//    if (nstepsOut == 1)   sprintf (str, ".cbltens");
//    else                  sprintf (str, ".cbltens.%03ld", step);
//    
//    this->print_VTK_CD_scal (mft, str, "CableTensileStress", (const Array1d**)&data);
//  }
//  
//  delete data;
//}


void Mesh :: print_characteristics_to_VTK (bool mft)
{
  long i;
  Stream *stream = new Stream(STRM_tixel);
  VTKdataset ds = this->giveVTKDS();
  
  long Nn = Pjnts();
  long Ne = Elems();
  long NnNe = (Nn>Ne ? Nn : Ne);
  
  //*  ALLOCATION  *//
  this->resize_Lata (NnNe);
  
  //*  ************************
  //*  ***   START + BODY   ***
  //*  ************************
  print_VTK_START (stream, this, ".ctrl", ds);
  print_VTK_body  (stream      , ds);
  
  
  //*  **********************
  //*  ***   POINT DATA   ***
  //*  **********************
  print_VTK_init_point_data (stream, Nn);
  
  //* support
  for (i=0; i<Nn; i++)  Nodes(i)->give_pointAttribs()->assemble_full_constrained (Lata[i]);      print_VTK_data_body (stream, Nn, "supports_forces",  (GPA<Array1d>*)&Lata, 'v');
  for (i=0; i<Nn; i++)  Lata[i]->shift(3);                                                       print_VTK_data_body (stream, Nn, "supports_moments", (GPA<Array1d>*)&Lata, 'v');
  for (i=0; i<Nn; i++)  Lata[i]->deshift();
  
  this->resize_Lata_members(1);
  for (i=0; i<Nn; i++)  Lata[i][0][0] = Nodes(i)->give_pointAttribs()->is_supported();           print_VTK_data_body (stream, Nn, "supports_bool",     (GPA<Array1d>*)&Lata, 's');
  
  //* end point data
  stream->relink_up();
  
  
  //* *********************
  //* ***   CELL DATA   ***
  //* *********************
  print_VTK_init_cell_data (stream, Ne);
  
  this->resize_Lata_members(1);
  
  //* material
  for (i=0; i<Ne; i++)  Lata[i][0][0] = Elems[i]->give_elemAttribs()->give_mat_ID();       print_VTK_data_body (stream, Ne, "materials",      (GPA<Array1d>*)&Lata, 's');
  for (i=0; i<Ne; i++)  Lata[i][0][0] = Elems[i]->give_elemAttribs()->give_cs ()->give_ID();    print_VTK_data_body (stream, Ne, "cross-sections", (GPA<Array1d>*)&Lata, 's');
  
  //*  *********************
  //*  ***      END      ***
  //*  *********************
  print_VTK_FINISH (stream);
  delete stream;
}


//*  *****************************************************************************************
//*  ***                                    MADAPTIVITY                                    ***
//*  *****************************************************************************************
void Mesh :: adaptivity (void)
{
#ifdef __MEER_MODULE
  
  Pd->sodriver()->print_info_message (SODE_start_blue_arrow, "Adaptivity - relative stress error estimate");
  
  // allocate
  MEER_interface *ada = new MEER_interface;
  
  // run meer
  ada->local_initialize(this);
  ada->initialize();
  double err = ada->solve(Pd->give_adpa()->give_required_error());
  
  Pd->give_adpa()->set_error(err);
  
  // results at elem
  const double* errori_l = ada->give_error_at_elements_local();
  const double* errori_g = ada->give_error_at_elements_global();
  const double* h_old_elem = ada->give_h_old_at_elems();
  const double* h_new_elem = ada->give_h_new_at_elems();
  Dvctr*** rnvals = ada->give_refined_node_values_2();
  
  // docas
  const double* eia = ada->give_abs_error_at_elements();
  
  long i, j, istep = 0;
  long nn = Pjnts();
  long ne = Elems();
  
  
  //* save results to element
  this->set_RTE(RTE_ada_error_pct_l, true);
  this->set_RTE(RTE_ada_error_pct_g, true);
  this->set_RTE(RTE_ada_h_act,       true);
  this->set_RTE(RTE_ada_h_new,       true);
  this->set_RTE(RTE_ada_normal,      true);
  
  // docas
  this->set_RTE(RTE_ada_eia,      true);
  
  for (i=0; i<ne; i++) {
    // save results
    FElems(i)->set_result(errori_l[i],                 istep, RTE_ada_error_pct_l);
    FElems(i)->set_result(errori_g[i],                 istep, RTE_ada_error_pct_g);
    FElems(i)->set_result(h_old_elem[i],               istep, RTE_ada_h_act);
    FElems(i)->set_result(h_new_elem[i],               istep, RTE_ada_h_new);
    FElems(i)->set_result(ada->give_normal_at_elem(i), istep, RTE_ada_normal);
    
    // docas
    FElems(i)->set_result(eia[i],                      istep, RTE_ada_eia);
  }
  
  //* save results to node
  Dvctr h_old_node(nn);
  Dvctr h_new_node(nn);
  double sizei, size, hold, hnew;
  
  // zprumeruju h_* do uzlu
  for (i=0; i<nn; i++) {
    const GPA<const Element>* supelems = Nodes(i)->give_superelems();
    
    if (supelems[0]()) {
      //* vazeny prumer
      size = hold = hnew = 0;
      for (j=0; j<supelems[0](); j++) {
    sizei = supelems[0][j]->give_lav();
    
    size += sizei;
    hold += sizei * h_old_elem[supelems[0][j]->give_ID()];
    hnew += sizei * h_new_elem[supelems[0][j]->give_ID()];
      }
      
      h_old_node[i] = hold / size;
      h_new_node[i] = hnew / size;
      
      // //* nejmensi
      // for (j=0; j<supelems[0](); j++) {
      //   hnew = h_new_elem[supelems[0][j]->give_ID()];
      //   if (h_new_node[i] > hnew)
      //    h_new_node[i] = hnew;
      // }
    }
    else {
      h_old_node[i] = h_new_node[i] = 0.0;
    }
  }
  
  // tisk ada_refder, normal, local coord system
  if (0) {
    // save results to node
    this->set_RTN(RTN_ada_refder, true);
    
    //const VectoR *v;
    //VectoR v0;  v0.x = v0.y = v0.z = 0.0;
    const Dmtrx *m;
    
    for (i=0; i<this->nRegions; i++)
      for (j=0; j<nn; j++) {
    Nodes(j)->add_resultE (rnvals[i][j], i, istep, RTE_global_stress);
    rnvals[i][j] = NULL;
    
    //v = ada->give_normal_at_node_and_region(j, i);
    //if (v == NULL)  v = &v0;
    //Nodes(j)->set_resultE (v, i, istep, RTE_ada_normal);
    m = ada->give_rotmat_at_node(j);
    Nodes(j)->set_resultE (m, i, istep, RTE_ada_normal);
      }
  }
  
  // save results
  this->set_RTN(RTN_ada_h_act, true);
  this->set_RTN(RTN_ada_h_new, true);
  
  for (i=0; i<nn; i++) {
    Nodes(i)->set_resultN(h_old_node[i], istep, RTN_ada_h_act);
    Nodes(i)->set_resultN(h_new_node[i], istep, RTN_ada_h_new);
  }
  
  
  //* print BackGroundMesh
  GPA<const FElement> Beams(0,ne), Triangles(0,ne), Quads(0,ne), Tetras(0,ne), Bricks(0,ne);
  
  for (i=0; i<ne; i++)
    switch (Elems[i]->give_cellGeom()) {
    case CG_Line:        Beams.add(FElems(i));      break;
    case CG_Triangle:    Triangles.add(FElems(i));  break;
    case CG_Quadrangle:  Quads.add(FElems(i));      break;
    case CG_Tetrahedron: Tetras.add(FElems(i));     break;
    case CG_Hexahedron:  Bricks.add(FElems(i));     break;
    default: _errorr2 ("Unknown element geom type (element number %d)", i+1);
    }
  
  FILE *bgm = _openFilePNS ("w", "T3d input", Pd->give_OUT_moFILE()->give_path(), Pd->give_OUT_moFILE()->give_base(), ".T3d.bgm");
  
  fprintf(bgm, "7 1\n");
  fprintf(bgm, "%ld %ld %ld %ld %ld %d %d %ld\n", Pjnts(), Beams(), Triangles(), Quads(), Tetras(), 0, 0, Bricks()); // 0 = pyramids, wedges
  
  // nodes
  for (i=0; i<nn; i++)
    fprintf(bgm, "%ld %e %e %e  %e\n", i+1, Pjnts[i]->give_coord(0), Pjnts[i]->give_coord(1), Pjnts[i]->give_coord(2), h_new_node[i]);
  
  //
  long ii = 1;
  
  if (Triangles()) {
    for (i=0; i<Triangles(); i++) {
      fprintf(bgm, "%ld", i+1);
      for (j=0; j<3; j++)  fprintf(bgm, " %ld", Triangles[i]->give_node(j)->give_ID()+1);
      fprintf(bgm, "\n");
    }
  }
  
  if (Quads()) {
    for (i=0; i<Quads(); i++) {
      fprintf(bgm, "%ld", i+1);
      for (j=0; j<4; j++)  fprintf(bgm, " %ld", Quads[i]->give_node(j)->give_ID()+1);
      fprintf(bgm, "\n");
    }
  }
  
  if (Beams())  errol;
  if (Tetras()) errol;
  if (Bricks()) errol;
  
  fclose (bgm);
  
  
  delete ada;
  
  //Pd->sodriver()->print_info_message_core (SODE_green_ok);
  Pd->sodriver()->print_info_message (SODE_end, "%.2lf", err);
  
#else
  _errorr("MEER is not linked.");
#endif
}

#ifdef __MEER_MODULE
void Mesh :: regions_init (void)
{
  //* smyce pres vsechny elementy - zatim jsou vsechny v regionu 0
  this->nRegions = Pd->give_cCS();
  for (long i=0; i<Elems(); i++)
    FElems(i)->set_regid (FElems(i)->give_elemAttribs()->give_cs()->give_ID());
  
}
#endif

double Mesh :: volume_core (void) const
{
  //
  if (! Pd->give_fulldata()) errol;
  
  long i;
  double cs, vol = 0.0;
  for (i=0; i<Elems(); i++) {
    switch (Elems[i]->give_dimension()) {
    case 1:  if ((cs = Elems[i]->give_elemAttribs()->give_cs()->give_area())      <= 0.0) errol;        break;
    case 2:  cs = Elems[i]->give_elemAttribs()->give_cs()->give_thickness();  if (cs < 0.0)  cs = 0.0;  break;
    case 3:  cs = 1.0;                                                                                  break;
    default: errol;
    }
    
    vol += cs * Elems[i]->give_lav();
  }
  
  return vol;
}

} // namespace midaspace