vtk.cpp

Go to the documentation of this file.

//********************************************************************************************************
// code:                          ###  ###  #####   ####  ##  ##
//                       ##  ##   ########  ##     ##  ## ##  ##
//                       ##  ##   ## ## ##  ###    ##     ######
//                       ##  ##   ##    ##  ##     ##  ## ##  ##
//                       #####    ##    ##  #####   ####  ##  ##
//                       ##
//
// name:           vtk.cpp
// author(s):      Ladislav Svoboda
// language:       C, C++
// license:        This program is free software; you can redistribute it and/or modify
//                 it under the terms of the GNU Lesser General Public License as published by
//                 the Free Software Foundation; either version 2 of the License, or
//                 (at your option) any later version.
//
//                 This program is distributed in the hope that it will be useful,
//                 but WITHOUT ANY WARRANTY; without even the implied warranty of
//                 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
//                 GNU Lesser General Public License for more details.
//
//                 You should have received a copy of the GNU Lesser General Public License
//                 along with this program; if not, write to the Free Software
//                 Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
//********************************************************************************************************

#include "vtk.h"

#include "macros.h"
#include "librw.h"
#include "mathlib.h"

#include <time.h>
#include <math.h>


namespace mumech {


//*  ***   START   ***
void print_VTK_START (Stream *stream, const char *filename)
{
  if (stream->isFile()) {
      stream->open(STRM_file, "w", filename);
  }
  else {
    // xml document
    XMLDocument *doc = new XMLDocument(true, COLLAPSE_WHITESPACE);
    
    doc->InsertEndChild( doc->NewDeclaration() );
    
    XMLElement *VTKFile = doc->NewElement("VTKFile");   doc->InsertEndChild(VTKFile);
    VTKFile->SetAttribute("type", "UnstructuredGrid");
    VTKFile->SetAttribute("version", "0.1");
    VTKFile->SetAttribute("byte_order", "LittleEndian");
    
    // initiate stream
    stream->open(STRM_tixel, "w", filename, VTKFile -> InsertEndChild(doc->NewElement("UnstructuredGrid"))->ToElement()->InsertEndChild(doc->NewElement("Piece")) );
  }
}

//*  ***   HEAD   ***
void print_VTK_head (FILE *stream, const char *type)
{
  struct tm tmv;
  time_t today;
  
  //time and dates
  time( &today );
  tmv = * localtime( &today );
  
  // printing of the header
  fprintf( stream, "# vtk DataFile Version 3.0\n");
  fprintf( stream, "%s ", type);
  fprintf( stream, "VTK file created by muMECH on %d/%d/%d at %d:%d:%d\n",
       tmv.tm_year + 1900, tmv.tm_mon +1, tmv.tm_mday, tmv.tm_hour, tmv.tm_min, tmv.tm_sec );
  fprintf( stream, "ASCII\n" );
  fprintf( stream, "DATASET UNSTRUCTURED_GRID\n" );
}

//*  ***   NODES   ***
XMLElement* print_VTK_nodes_head (Stream *stream, long n)
{
  XMLElement *da = NULL;
  
  if (stream->isFile())
    fprintf (stream->file(), "POINTS %ld float\n", n);
  else {
    stream->tixel()->ToElement()->SetAttribute("NumberOfPoints", (int)n);
    //
    da = stream->tix_doc()->NewElement("DataArray");
    da->SetAttribute("format", "ascii");
    da->SetAttribute("type", "Float32");
    da->SetAttribute("NumberOfComponents", "3");
    
    stream->tixel()->InsertEndChild(stream->tix_doc()->NewElement("Points") );
    stream->tixel()->LastChild()->InsertEndChild(da);
  }
  
  return da;
}


//*  ***   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   ***
XMLElement* print_VTK_data_head (Stream *stream, const char *name, char type, char format, int slen)
{
  XMLElement *da = NULL;
  
  char flint[8];
  if      (format == 'i') { if (stream->isFile()) sprintf(flint,"int");    else sprintf(flint,"Int32");   }
  else if (format == 'f') { if (stream->isFile()) sprintf(flint,"float");  else sprintf(flint,"Float32"); }
  else errol;
  
  if (stream->isFile()) {
    if      (type == 's') { fprintf (stream->file(), "SCALARS %s %s %d\n", name, flint, slen);  fprintf (stream->file(), "LOOKUP_TABLE default\n");  }
    else if (type == 'v')   fprintf (stream->file(), "VECTORS %s %s\n", name, flint);
    else if (type == 't')   fprintf (stream->file(), "TENSORS %s %s\n", name, flint);
    else    _errorr ("error");
  }
  else {
    da = stream->tix_doc()->NewElement("DataArray");  stream->tixel()->InsertEndChild(da);
    da->SetAttribute("format", "ascii");
    da->SetAttribute("type", flint);
    da->SetAttribute("Name", name);
    if      (type == 's')  da->SetAttribute("NumberOfComponents", "1");
    else if (type == 'v')  da->SetAttribute("NumberOfComponents", "3");
    else if (type == 't')  da->SetAttribute("NumberOfComponents", "9");
    else    _errorr ("error");
  }
  
  return da;
}

void print_values (char *auxs, int n, const double *values, int precision)
{
  if (n==0) _errorr("");
  
  int n1 = 5 + precision + NUM_DIGITS_IN_PRINTED_EXPONENT;
  int pos = 0;
  
  sprintf (auxs+pos, "%*.*e", n1, precision, values[0]);  pos += n1;
  
  for (int i=1; i<n; i++) {
    sprintf (auxs+pos, " %*.*e", n1, precision, values[i]);
    pos += n1+1;
  }
  
  //sprintf (auxs+pos, "\n");
  if (auxs[pos] != '\0')
    errol;
}

void print_values_Scalar (char *auxs, const double *values, long len, int precision)
{
  int n1 = 5 + precision + NUM_DIGITS_IN_PRINTED_EXPONENT;
  
  for (long i=0; i<len; i++)
    sprintf (auxs+i*(n1+1), "%*.*e ", n1, precision, values[i]);

  auxs[len*(n1+1) - 1] = '\0';
}
void print_values_Vector (char *auxs, const double *values, bool twodim, int precision)
{
  int n1 = 5 + precision + NUM_DIGITS_IN_PRINTED_EXPONENT;
  
  if (twodim) {
    sprintf (auxs+0*(n1+1), "%*.*e ", n1, precision, values[0]);  sprintf (auxs+1*(n1+1), "%*.*e ", n1, precision, values[1]);  sprintf (auxs+2*(n1+1), "%*.*e", n1, precision, 0.0      );
  }
  else {
    sprintf (auxs+0*(n1+1), "%*.*e ", n1, precision, values[0]);  sprintf (auxs+1*(n1+1), "%*.*e ", n1, precision, values[1]);  sprintf (auxs+2*(n1+1), "%*.*e", n1, precision, values[2]);
  }
}
void print_values_Tensor (char *auxs, const double *values, bool twodim, int precision)
{
  int n1 = 5 + precision + NUM_DIGITS_IN_PRINTED_EXPONENT;
  
  if (twodim) {
    sprintf (auxs+0*(n1+1), "%*.*e ", n1, precision, values[0]);  sprintf (auxs+1*(n1+1), "%*.*e ", n1, precision, values[2]);  sprintf (auxs+2*(n1+1), "%*.*e ", n1, precision, 0.0      );
    sprintf (auxs+3*(n1+1), "%*.*e ", n1, precision, values[2]);  sprintf (auxs+4*(n1+1), "%*.*e ", n1, precision, values[1]);  sprintf (auxs+5*(n1+1), "%*.*e ", n1, precision, 0.0      );
    sprintf (auxs+6*(n1+1), "%*.*e ", n1, precision, 0.0      );  sprintf (auxs+7*(n1+1), "%*.*e ", n1, precision, 0.0      );  sprintf (auxs+8*(n1+1), "%*.*e",  n1, precision, 0.0      );
  }
  else {
    sprintf (auxs+0*(n1+1), "%*.*e ", n1, precision, values[0]);  sprintf (auxs+1*(n1+1), "%*.*e ", n1, precision, values[3]);  sprintf (auxs+2*(n1+1), "%*.*e ", n1, precision, values[5]);
    sprintf (auxs+3*(n1+1), "%*.*e ", n1, precision, values[3]);  sprintf (auxs+4*(n1+1), "%*.*e ", n1, precision, values[1]);  sprintf (auxs+5*(n1+1), "%*.*e ", n1, precision, values[4]);
    sprintf (auxs+6*(n1+1), "%*.*e ", n1, precision, values[5]);  sprintf (auxs+7*(n1+1), "%*.*e ", n1, precision, values[4]);  sprintf (auxs+8*(n1+1), "%*.*e",  n1, precision, values[2]);
  }
}

void print_auxs (bool lgc, Stream *stream, XMLElement *da, const char *auxs)
{
  if (auxs[20000] != '\0')
    _errorr("line length exceeds 20000, something is wrong");

  if (lgc)  fprintf (stream->file(), "%s\n", auxs);
  else      da->InsertEndChild( stream->tix_doc()->NewText(auxs) );
}

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


void scan_DATA_field_head (Stream *stream, Stream *strm, const char *str, char type, int nexc, char format, const char *name)
{
  char flint[8];
  if      (format == 'i') { if (stream->isFile()) sprintf(flint,"int");    else sprintf(flint,"Int32");   }
  else if (format == 'f') { if (stream->isFile()) sprintf(flint,"float");  else sprintf(flint,"Float32"); }
  else errol;
  
  char auxs[255];
  sprintf (auxs, "Error in section %s", name);
  
  if (stream->isFile()) {
    SP_scan_expected_word_exit   (str, flint, auxs, CASE);
    if (type == 'S') {
      SP_scan_expected_number_exit (str, nexc, auxs);
      FP_scan_expected_word_exit   (stream->file(), "LOOKUP_TABLE default", auxs, CASE);
      FP_skip_line (stream->file());
    }
    strm->redefine( stream->file() );
  }
  else {
    strm->redefine( XP_giveDAtext  (stream->tixel(), 0, false, "ascii", flint, name, &nexc));
  }
}



bool cmp_vtk_file (const char *file1, const char *file2, int verbose)
{
  double ABSZERO = 1.e-15;
  double RELZERO = 2.e-2;
  
  char *filename;
  
  // open stream with automatic detection VTK legaci/XML file format
  filename = new char[255];
  strcpy(filename, file1);
  Stream *stream1 = new Stream();
  stream1->open(STRM_void, "r", (const char*&)filename); // do not delete filename, it is deleted inside of the function
  filename = NULL;
  
  // open stream with automatic detection VTK legaci/XML file format
  filename = new char[255];
  strcpy(filename, file2);
  Stream *stream2 = new Stream();
  stream2->open(STRM_void, "r", (const char*&)filename); // do not delete filename, it is deleted inside of the function
  filename = NULL;
  
  
  bool status = false; //prnt = false;
  
  //
  char *WORD1 = NULL;
  char *WORD2 = NULL;
  //const XMLElement *xnel = NULL;
  //
  bool lgc = stream1->isFile();
  if (lgc != stream2->isFile())
    status++;
  
  //bool amd = false;
  //
  if (lgc) WORD1 = new char[255];
  if (lgc) WORD2 = new char[255];
  //
  
  // HEAD OD FILE
  if (lgc) {
    // 1. line
    FP_skip_line (stream1->file());
    FP_skip_line (stream2->file());
    
    // 2. line
    FP_scan_word (stream1->file(), WORD1);   FP_skip_line (stream1->file());
    FP_scan_word (stream2->file(), WORD2);   FP_skip_line (stream2->file());
    status +=  (strcmp(WORD1, WORD2) == 0) ? 0 : 1;
    
    // 3. line
    FP_scan_expected_word_exit (stream1->file(), "ASCII", "Invalid structure of the given VTK file", true);   FP_skip_line (stream1->file());
    FP_scan_expected_word_exit (stream2->file(), "ASCII", "Invalid structure of the given VTK file", true);   FP_skip_line (stream2->file());
    
    // 4. line
    FP_scan_expected_word_exit (stream1->file(), "DATASET UNSTRUCTURED_GRID", "Invalid structure of the given VTK file", true);   FP_skip_line (stream1->file());
    FP_scan_expected_word_exit (stream2->file(), "DATASET UNSTRUCTURED_GRID", "Invalid structure of the given VTK file", true);   FP_skip_line (stream2->file());
  }
  else {
    _errorr("not implemented");
  //  // check head of xml file and read possible ctrl section
  //  // return XMLElement "Piece"
  //  const XMLNode *parent = stream->tixnod();
  //  
  //  //* declaration - delete if present
  //  if (parent->FirstChild()->ToDeclaration() != NULL)  stream->tixnod()->DeleteChild( stream->tixnod()->FirstChild() );
  //  
  //  //* comment - delete if present
  //  while (true)
  //    if (parent->FirstChild()->ToComment() != NULL)  stream->tixnod()->DeleteChild( stream->tixnod()->FirstChild() );
  //    else                                            break;
  //  
  //  //
  //  parent = XP_give_unique_expected_elem (parent, "VTKFile");
  //  stream->relink_downF();
  //  
  //  ///
  //  XP_check_expected_attribute (parent, "type", "UnstructuredGrid");
  //  XP_check_expected_attribute (parent, "version", "0.1");
  //  XP_check_expected_attribute (parent, "byte_order", "LittleEndian");
  //  
  //  //
  //  parent = XP_give_unique_expected_elem (parent, "UnstructuredGrid", false);
  //  stream->relink_downF();
  //  
  //  // read AppendedData section
  //  parent = parent->NextSibling();
  //  if (! parent)
  //    _errorr("No AppendedData section");
  //  
  //  if (_STRCMP(parent->Value(), "AppendedData") != 0)     _errorr ("name of xelement is not \"AppendedData\"");
  //  
  //  if (! (parent = parent->FirstChild ()))  _errorr("No child in AppendedData");
  //  if (   parent->Value()[0] != '_'      )  _errorr("There is no character \'_\' at the start of section AppendedData");
  //  parent = parent->NextSibling();
  //  
  //  if (_STRCMP(parent->Value(), "Input_Type") != 0)  _errorr("name of xelement is not Input_Type");
  //  amd = this->file_type_s2e (parent->ToElement()->GetText());
  //  
  //  while (parent->NextSibling()) {
  //    parent = parent->NextSibling();
  //    const char *text = parent->Value();
  //    long auxl;
  //    
  //    if (_STRCMP(text, "SBA_mode") == 0) {
  //    sscanf (parent->ToElement()->GetText(), "%ld", &auxl);
  //    SBA_optimized = (bool)(auxl-1);
  //    }
  //    else if (_STRCMP(text, "Matrix_record") == 0) {
  //      double e, nu;
  //      sscanf (parent->ToElement()->GetText(), "%lf %lf", &e, &nu);
  //      this->set_matrixRecord (e, nu);
  //    }
  //    else if (_STRCMP(text, "Min_axes_diff") == 0) { double m;  sscanf (parent->ToElement()->GetText(), "%lf", &m);  this->set_semiaxes_min_difference(m); }
  //    else if (_STRCMP(text, "Max_axes_diff") == 0) { double m;  sscanf (parent->ToElement()->GetText(), "%lf", &m);  this->set_semiaxes_max_difference(m); }
  //  }
  //  
  //  // go to Piece
  //  XP_give_unique_expected_elem (stream->tixnod(), "Piece");
  //  stream->relink_downF();
  //  xnel = stream->tixel();
  }
  
  
  //Stream auxstrm;
  //int auxi, lc, les, ges;
  long nn, NN, NE, auxl1, auxl2, noel1, noel2; // , offset;
  double auxd1, auxd2;
  char ed='-';
  //const char *typs, *data=NULL, *offs=NULL;
  const char  *str1, *str2;
  char LINE1[1023], LINE2[1023];
  int lLINE=1023;
  VTKrange range = VTKR_void;
  //
  
  //cno = 0;
  //prnt = true;
  while (true) {
    //* key WORD
    if (lgc) {
      if (fgets (LINE1, lLINE, stream1->file()) == NULL)  { delete [] WORD1;  WORD1 = NULL; }
      if (fgets (LINE2, lLINE, stream2->file()) == NULL)  { delete [] WORD2;  WORD2 = NULL; }
      if (WORD1 == NULL  ||  WORD2 == NULL) {
    if (WORD1 != WORD2)
      status++;
    
    break;
      }
      
      str1 = LINE1;  SP_scan_word (str1, WORD1);
      str2 = LINE2;  SP_scan_word (str2, WORD2);
      status +=  (strcmp(WORD1, WORD2) == 0) ? 0 : 1;
    }
    else {
      // if (prnt) { if (! stream->relink_downF())  errol; else prnt = false; }
      // else      { if (! stream->relink_next())   break; }
      // WORD = (char *)stream->tixel()->Value();
    }
    
    //*  ******************
    //*  ***   POINTS   ***
    //*  ******************
    if      (_STRCMP(WORD1, "Points") == 0) {
      if (range == VTKR_void)  range = VTKR_points;  else errol;
      
      // number of nodes
      if (lgc) {
    SP_scan_number (str1, noel1);
    SP_scan_expected_word2_exit  (str1, "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));
      }
      
      // number of nodes
      if (lgc) {
    SP_scan_number (str2, noel2);
    SP_scan_expected_word2_exit  (str2, "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));
      }
      status +=  (noel1 == noel2) ? 0 : 1;
      NN = noel1;
      
      // coords
      noel1 *= 3;
      while (noel1) { noel1--;
    if (lgc) { status += ( fscanf (stream1->file(), "%lf", &auxd1) != 1 ); }
    // else     { status += ( sscanf (data,           "%lf %lf %lf", inclRec[noIncl]->origin, inclRec[noIncl]->origin+1, inclRec[noIncl]->origin+2) == 3 );  if (inclRec[noIncl]->origin[2] != 0.0)  _errorr("3. coordinates is not 0.0");  status += (SP_skip_word(data) && SP_skip_word(data) && SP_skip_word(data)); }
    if (lgc) { status += ( fscanf (stream2->file(), "%lf", &auxd2) != 1 ); }
    // else     { status += ( sscanf (data,           "%lf %lf %lf", inclRec[noIncl]->origin, inclRec[noIncl]->origin+1, inclRec[noIncl]->origin+2) == 3 );  if (inclRec[noIncl]->origin[2] != 0.0)  _errorr("3. coordinates is not 0.0");  status += (SP_skip_word(data) && SP_skip_word(data) && SP_skip_word(data)); }
        status += areSame (ABSZERO, RELZERO, auxd1, auxd2) ? 0 : 1;
    
    
      }
      
      if (lgc)  FP_skip_line (stream1->file());
      if (lgc)  FP_skip_line (stream2->file());
    }
    //*  *************************************
    //*  ***   UNSTRUCTURED GRID - CELLS   ***
    //*  *************************************
    else if (_STRCMP(WORD1, "CELLS") == 0) {
      if (range == VTKR_points)  range = VTKR_cells;  else errol;
      
      // number of cells
      if (lgc) {
        SP_scan_number (str1, noel1);
        SP_scan_number (str1, auxl1);
      }
      //    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);
      //    }
      // number of cells
      if (lgc) {
        SP_scan_number (str2, noel2);
        SP_scan_number (str2, auxl2);
      }
      status +=  (noel1 == noel2) ? 0 : 1;
      status +=  (auxl1 == auxl2) ? 0 : 1;
      NE = noel1;
      
      //if (!lgc) offset = 0;
      // elem nodes
      while (noel1) { noel1--;
        // nn
        if (lgc) { fscanf (stream1->file(),"%ld",&auxl1); }
        //else     { nn = offset;  status += SP_scan_number (offs, offset);  nn = offset - nn; }
        //
        if (lgc) { fscanf (stream2->file(),"%ld",&auxl2); }
    status +=  (auxl1 == auxl2) ? 0 : 1;
    
    nn = auxl1;
    while (nn) { nn--;
      if (lgc) fscanf (stream1->file(),"%ld", &auxl1);
      //else     status += SP_scan_number (data, auxl);
      if (lgc) fscanf (stream2->file(),"%ld", &auxl2);
      
      status +=  (auxl1 == auxl2) ? 0 : 1;
    }
      }
      
      if (lgc) FP_skip_line (stream1->file(), 1);
      if (lgc) FP_skip_line (stream2->file(), 1);
      
      // types of elems
      if (lgc) {
    fgets (LINE1, lLINE, stream1->file());  str1=LINE1;
    SP_scan_expected_word_exit   (str1, "CELL_TYPES", "Invalid structure of the given VTK file", CASE);
    SP_scan_expected_number_exit (str1,        noel2, "Invalid structure of the given VTK file, the number following CELL_TYPES");
      }
      // types of elems
      if (lgc) {
    fgets (LINE2, lLINE, stream2->file());  str2=LINE2;
    SP_scan_expected_word_exit   (str2, "CELL_TYPES", "Invalid structure of the given VTK file", CASE);
    SP_scan_expected_number_exit (str2,        noel2, "Invalid structure of the given VTK file, the number following CELL_TYPES");
      }
      
      while (noel2) { noel2--;
    if (lgc)  fscanf (stream1->file(),"%ld",&auxl1);
    //else      status += SP_scan_number (typs, auxl);
    if (lgc)  fscanf (stream2->file(),"%ld",&auxl2);
    //else      status += SP_scan_number (typs, auxl);
    status +=  (auxl1 == auxl2) ? 0 : 1;
      }
      
      if (lgc) FP_skip_line (stream1->file(), 1);
      if (lgc) FP_skip_line (stream2->file(), 1);
    }
    
    //*  *******************************
    //*  ***   POLYDATA - ELEMENTS   ***
    //*  *******************************
    else if (_STRCMP(WORD1, "LINES") == 0  ||  _STRCMP(WORD1, "POLYGONS" ) == 0  ||  _STRCMP(WORD1, "POLYS") == 0  ||  _STRCMP(WORD1, "Vertices") == 0  ||  _STRCMP(WORD1, "Verts") == 0  ||  _STRCMP(WORD1, "Strips") == 0) {
      _errorr2("Invalid structure of the given VTK file, the key word %s is not supported", WORD1);
    }
    //*  **********************
    //*  ***   POINT_DATA   ***
    //*  **********************
    else if (_STRCMP(WORD1, "POINT_DATA") == 0  ||  _STRCMP(WORD1, "PointData") == 0) {
      if (range == VTKR_cells  ||  range == VTKR_fields)  range = VTKR_pd;  else errol;
      ed = 'p';
      // head
      if (lgc) { SP_scan_expected_number_exit (str1, NN, "The number following POINT_DATA is not equal to number of inclusions"); }
      //else     { prnt = true; }
      if (lgc) { SP_scan_expected_number_exit (str2, NN, "The number following POINT_DATA is not equal to number of inclusions"); }
      //else     { prnt = true; }
    }
    //*  *********************
    //*  ***   CELL_DATA   ***
    //*  *********************
    else if (_STRCMP(WORD1, "CELL_DATA") == 0  ||  _STRCMP(WORD1, "CellData") == 0) {
      if (range == VTKR_cells || range == VTKR_data)  range = VTKR_cd;  else errol;
      ed = 'c';
      // head
      if (lgc) { SP_scan_expected_number_exit (str1, NE, "The number following CELL_DATA is not equal to number of inclusions"); }
      //else     { prnt = true; }
      if (lgc) { SP_scan_expected_number_exit (str2, NE, "The number following CELL_DATA is not equal to number of inclusions"); }
      //else     { prnt = true; }
    }
    //*  ******************
    //*  ***   FIELDS   ***
    //*  ******************
    else if (_STRCMP(WORD1, "FIELD") == 0) {
      if (range == VTKR_cells || range == VTKR_data)  range = VTKR_fields;  else errol;
      
      if (!lgc) _errorr("FIELD in nonlegacy");
      
      // head
      if (lgc)   SP_scan_word (str1, WORD1);
      if (lgc)   SP_scan_word (str2, WORD2);
      status +=  (strcmp(WORD1, WORD2) == 0) ? 0 : 1;
      
      // number of fields
      SP_scan_number (str1, noel1);
      SP_scan_number (str2, noel2);
      status +=  (noel1 == noel2) ? 0 : 1;
      
      for (int i=0; i<noel1; i++) {
    FP_scan_word (stream1->file(), WORD1);
    FP_scan_word (stream2->file(), WORD2);
    status +=  (strcmp(WORD1, WORD2) == 0) ? 0 : 1;
    
    fscanf (stream1->file(), "%ld", &auxl1);
    fscanf (stream2->file(), "%ld", &auxl2);
    status +=  (auxl1 == auxl2) ? 0 : 1;
    
    nn = auxl1;
    
    fscanf (stream1->file(), "%ld", &auxl1);
    fscanf (stream2->file(), "%ld", &auxl2);
    status +=  (auxl1 == auxl2) ? 0 : 1;
    
    nn *= auxl1;
    
    FP_scan_word (stream1->file(), WORD1);
    FP_scan_word (stream2->file(), WORD2);
    status +=  (strcmp(WORD1, WORD2) == 0) ? 0 : 1;
    
    if (_STRCMP(WORD1, "int") == 0  ||  _STRCMP(WORD1, "Int32") == 0)
      while (nn) { nn--;
        status += ( fscanf (stream1->file(), "%ld", &auxl1) != 1 );
        status += ( fscanf (stream2->file(), "%ld", &auxl2) != 1 );
        status += (auxl1 == auxl2) ? 0 : 1;
      }
    else if (_STRCMP(WORD1, "float") == 0  ||  _STRCMP(WORD1, "Float32") == 0)
      while (nn) { nn--;
        status += ( fscanf (stream1->file(), "%lf", &auxd1) != 1 );
        status += ( fscanf (stream2->file(), "%lf", &auxd2) != 1 );
            status += areSame (ABSZERO, RELZERO, auxd1, auxd2) ? 0 : 1;
      }
      }
    }
    //*  **************************
    //*  ***   DATA - SCALARS   ***
    //*  **************************
    else if (_STRCMP(WORD1, "SCALARS") == 0  ||  _STRCMP(WORD1, "VECTORS") == 0  ||  _STRCMP(WORD1, "TENSORS") == 0  ||  _STRCMP(WORD1, "DataArray") == 0) {
      if (range == VTKR_pd || range == VTKR_cd || range == VTKR_data)  range = VTKR_data;  else errol;
      
      char type = WORD1[0]; // save type of data - Scalar, Vector, Tensor, DataArray
      // head
      if (lgc)   SP_scan_word (str1, WORD1);
      //else     { WORD = (char *)stream->tixel()->ToElement()->Attribute("Name");  if (!WORD)  _errorr ("there is no attribute \"Name\""); }
      if (lgc)   SP_scan_word (str2, WORD2);
      //else     { WORD = (char *)stream->tixel()->ToElement()->Attribute("Name");  if (!WORD)  _errorr ("there is no attribute \"Name\""); }
      status +=  (strcmp(WORD1, WORD2) == 0) ? 0 : 1;
      
      // type
      if (lgc)   SP_scan_word (str1, WORD1);
      if (lgc)   SP_scan_word (str2, WORD2);
      status +=  (strcmp(WORD1, WORD2) == 0) ? 0 : 1;
      
      if (lgc) {
    if (type == 'S') {
      SP_scan_number (str1, auxl1);
      SP_scan_number (str2, auxl2);
      status +=  (auxl1 == auxl2) ? 0 : 1;
      nn = auxl1;
          status += FP_scan_expected_word (stream1->file(), "LOOKUP_TABLE default", CASE) ? 0 : 1;
          status += FP_scan_expected_word (stream2->file(), "LOOKUP_TABLE default", CASE) ? 0 : 1;
    }
    else if (type == 'V')  nn = 3;
    else if (type == 'T')  nn = 9;
    else if (type == 'F')  nn = 0;
    else errol;
      }
      //else     status += SP_scan_number (data, auxl);
      // strm->redefine( stream->file() );
      // else {
      // strm->redefine( XP_giveDAtext  (stream->tixel(), 0, false, "ascii", flint, name, &nexc));
      // }
      //FP_skip_line (stream->file());
      
      if      (ed == 'p')  nn = nn* NN;
      else if (ed == 'c')  nn = nn* NE;
      else errol;
      
      if (_STRCMP(WORD1, "int") == 0  ||  _STRCMP(WORD1, "Int32") == 0)
    while (nn) { nn--;
      if (lgc) { status += ( fscanf (stream1->file(), "%ld", &auxl1) != 1 ); }
      // else     { status += ( sscanf (data,           "%lf %lf %lf", inclRec[noIncl]->origin, inclRec[noIncl]->origin+1, inclRec[noIncl]->origin+2) == 3 );  if (inclRec[noIncl]->origin[2] != 0.0)  _errorr("3. coordinates is not 0.0");  status += (SP_skip_word(data) && SP_skip_word(data) && SP_skip_word(data)); }
      if (lgc) { status += ( fscanf (stream2->file(), "%ld", &auxl2) != 1 ); }
      // else     { status += ( sscanf (data,           "%lf %lf %lf", inclRec[noIncl]->origin, inclRec[noIncl]->origin+1, inclRec[noIncl]->origin+2) == 3 );  if (inclRec[noIncl]->origin[2] != 0.0)  _errorr("3. coordinates is not 0.0");  status += (SP_skip_word(data) && SP_skip_word(data) && SP_skip_word(data)); }
      status += (auxl1 == auxl2) ? 0 : 1;
    }
      else if (_STRCMP(WORD1, "float") == 0  ||  _STRCMP(WORD1, "Float32") == 0)
    while (nn) { nn--;
      if (lgc) { status += ( fscanf (stream1->file(), "%lf", &auxd1) != 1 ); }
      // else     { status += ( sscanf (data,           "%lf %lf %lf", inclRec[noIncl]->origin, inclRec[noIncl]->origin+1, inclRec[noIncl]->origin+2) == 3 );  if (inclRec[noIncl]->origin[2] != 0.0)  _errorr("3. coordinates is not 0.0");  status += (SP_skip_word(data) && SP_skip_word(data) && SP_skip_word(data)); }
      if (lgc) { status += ( fscanf (stream2->file(), "%lf", &auxd2) != 1 ); }
      // else     { status += ( sscanf (data,           "%lf %lf %lf", inclRec[noIncl]->origin, inclRec[noIncl]->origin+1, inclRec[noIncl]->origin+2) == 3 );  if (inclRec[noIncl]->origin[2] != 0.0)  _errorr("3. coordinates is not 0.0");  status += (SP_skip_word(data) && SP_skip_word(data) && SP_skip_word(data)); }
          status += areSame (ABSZERO, RELZERO, auxd1, auxd2) ? 0 : 1;
    }
      
      //*  *** DATA for CELLS ***
      if (lgc) FP_skip_line (stream1->file(), 1);
      //else { if (stream->tixel()->NextSibling() == NULL)  stream->relink_up(); }
      if (lgc) FP_skip_line (stream2->file(), 1);
    }
    else if (_STRCMP(WORD1, "") == 0) {;}
    else
      _errorr2 ("Invalid structure of the given VTK file, unexpected key word \"%s\"", WORD1);
    
    // 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 (lgc) delete [] WORD;
  
  stream1->close();
  stream2->close();
  delete stream1;
  delete stream2;
  //
  
  if (verbose > 0  &&  status >  0)  printf( " ERROR diffless %s %s \n", file1, file2);
  if (verbose > 1  &&  status == 0)  printf( " OK    diffless %s %s \n", file1, file2);
  
  return (status > 0);
} // end of function Problem :: read_vtk_file (const char *vtkTopologyFile)



} // end of namespace mumech

/*end of file*/