// SpaFEDTe, a Template based C++ library for creating 
// Discontinuous Finite Element Spaces,
// Copyright (C) 2012 Lorenzo Alessio Botti

/* This library 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 3.0 of the License, or (at your option) any later version. */

/* This software 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 software; if not, a copy of the full */
/* GNU Lesser General Public License can be found at */
/* http://www.gnu.org/licenses/. */

// This implementation is mostly borrowed from MOAB readers

#ifndef __read_mesh_vtu__h
#define __read_mesh_vtu__h

#include "VtuMeshUtil.hpp"
#include "rapidxml.hpp"
#include "rapidxml_utils.hpp"
#include <string>
#include <map>
#include "Mesh.hpp"

using namespace SpaMesh;
using namespace Connectivity;

class ReadMeshVtu 
{
   
public:

  ReadMeshVtu();

  ~ReadMeshVtu();

  template<unsigned int dim>
  void load_file( const char *filename, Mesh<dim>& mesh);
  
  template<unsigned int dim>
  void load_mesh( const char *filename, Mesh<dim>& mesh);

private:
  template<unsigned int dim>
  unsigned long computeNodesSize(const Mesh<dim>& mesh)
  {
    return mesh.nodesCoords.cols();
  }

  template<unsigned int dim>
  unsigned long computeCellsSize(const Mesh<dim>& mesh)
  {
    unsigned long size = 0;
    for (unsigned int i = 0; i < mesh.elemMesh.elemNodesConnectivity.size(); i++)
      size += mesh.elemMesh.elemNodesConnectivity[i].cols();
    for (unsigned int i = 0; i < mesh.faceMesh.elemNodesConnectivity.size(); i++)
      size += mesh.faceMesh.elemNodesConnectivity[i].cols();
    return size;
  }
 
  template<unsigned int dim>
  void vtu_read_unstructured_piece (const rapidxml::xml_node<> *& piece, 
			                 Mesh<dim>& mesh, std::vector<ElemTypesMesh>& elemTypes);

  template<unsigned int dim>
  void read_vertices( const char *& pointCoords,
                           const long numberOfPoints,
                 	   const long numberOfComponents,
			   Mesh<dim>& mesh);
 
  template<unsigned int dim>
  void read_elements( std::map<std::string, const char *>& cellData,
                           const long numberOfCells,
			   Mesh<dim>& mesh, std::vector<ElemTypesMesh>& elemTypes);

  template<unsigned int dim>
  void read_element_faces_data(const char * dataArrayString, 
                                    long numberOfComponents,
        		            const std::vector<ElemTypesMesh>& elemTypes,
                                    ElemMesh<dim>& mesh);
                          
  template<unsigned int dim>
  void read_face_neighbors_data(const char * dataArrayString, 
        		             const std::vector<ElemTypesMesh>& elemTypes,
                                     FaceMesh<dim>& mesh);
                          
private:
  std::string mPartitionTagName;
};

template<unsigned int dim>
void ReadMeshVtu::load_mesh(const char *filename, Mesh<dim>& mesh)
{
  std::string fileNameString = std::string(filename);
  int dotindex = fileNameString.find_last_of("."); 
  std::string extension = fileNameString.substr(dotindex+1); 

  if (strcmp(extension.c_str(),"vtu") == 0)
    this->load_file(filename, mesh);
  else if (strcmp(extension.c_str(),"pvtu") == 0)
  {
    std::cout<<"Vtu reader: Wrong extensions "<<extension<<std::endl;
  }
}

template<unsigned int dim>
void ReadMeshVtu::load_file(const char *filename, Mesh<dim>& mesh)
{
  rapidxml::file<> fileToRead(filename);
  rapidxml::xml_document<> doc;    
  doc.parse<rapidxml::parse_trim_whitespace|rapidxml::parse_validate_closing_tags>(fileToRead.data());    

  const rapidxml::xml_node<>* vtkFile = doc.first_node();
  if (!vtkFile)
    std::cerr<<"VtuReader error: no xml nodes in file "<<filename;
  // Check header
  if (strcmp(vtkFile->name(),"VTKFile") != 0)
    std::cerr<<"VtuReader error in file "<<filename
             <<": not reading a VTKFile. First node is "<<vtkFile->name();
  if (vtkFile->first_attribute("type"))
  {
    if (strcmp(vtkFile->first_attribute("type")->value(),"UnstructuredGrid") != 0)
      std::cerr<<"VtuReader error in file "<<filename
               <<": not reading an unstructured grid. File type is "
               <<vtkFile->first_attribute("type")->value();
  }
  else
    std::cerr<<"VtuReader error in file "<<filename
             <<": missing type of VTKFile, continue reading anyway";

  if (vtkFile->first_attribute("version"))
  {
    if (strcmp(vtkFile->first_attribute("version")->value(),"0.1") != 0 &&
        strcmp(vtkFile->first_attribute("version")->value(),"1.0") != 0)
      std::cerr<<"VtuReader error in file "<<filename
               <<": version mismatch. Version is "
               <<vtkFile->first_attribute("version")->value();
  }
  else
    std::cerr<<"VtuReader error in file "<<filename
             <<": missing version of VTKFile, continue reading anyway";
  
 
  // Read a piece of unstructured grid 
  if (!vtkFile->first_node("UnstructuredGrid"))
    std::cerr<<"VtuReader error: no UnstructuredGrid in VTKfile "<<filename;
  const rapidxml::xml_node<>* piece = vtkFile->first_node("UnstructuredGrid")->first_node("Piece");
  if (!piece)
    std::cerr<<"VtuReader error: no UnstructuredGrid/Piece in VTKfile "<<filename;
  std::vector<ElemTypesMesh> elemTypes(0);
  this->vtu_read_unstructured_piece (piece, mesh, elemTypes);
    
  long elem_count = this->computeCellsSize(mesh);
  long vertex_count = this->computeNodesSize(mesh);
  
  // Read CellData
  const rapidxml::xml_node<>* cellData = piece->first_node("CellData");
  if (cellData)
    for(rapidxml::xml_node<> *dataArray=cellData->first_node("DataArray"); dataArray; dataArray=dataArray->next_sibling("DataArray"))
    {
      long nComponents = 0;
      if (dataArray->first_attribute("NumberOfComponents"))
        nComponents = atol(dataArray->first_attribute("NumberOfComponents")->value());
      else
        nComponents = 1;
      if (!dataArray->first_attribute("Name"))
      {
        std::cerr<<"VtuReader error: missing name of CellData DataArray, aborting.";
	std::abort();
      }
      if (!dataArray->first_attribute("type"))
      {
        std::cerr<<"VtuReader error: missing type of CellData DataArray "
                 <<dataArray->first_attribute("Name")->value()<<", aborting.";
	std::abort();
      }
      if (strcmp(dataArray->value(), "\0") == 0)
      {
        std::cerr<<"VtuReader error: missing data of CellData DataArray "
                 <<dataArray->first_attribute("Name")->value()<<", aborting.";
	std::abort();
      }
      if (strcmp(dataArray->first_attribute("Name")->value(), "ELEMENT_FACES_LID") == 0)
      {
        if (strcmp(dataArray->first_attribute("type")->value(), "Int32") == 0 ||
            strcmp(dataArray->first_attribute("type")->value(), "Int64") == 0)
	{
          this->read_element_faces_data(dataArray->value(), nComponents, elemTypes, mesh.elemMesh);
	}
        else
	{
          std::cerr<<"VtuReader error: wrong type for element faces local ids, aborting.";
   	  std::abort();
        }
      }
      if (strcmp(dataArray->first_attribute("Name")->value(), "ELEMENT_NEIGHBORS_BY_FACE_LID") == 0)
      {
	if (nComponents != 2)
	{
          std::cerr<<"VtuReader error: wrong number of components for face neighbors local ids, aborting.";
	  std::abort();
	}
        if (strcmp(dataArray->first_attribute("type")->value(), "Int32") == 0 ||
            strcmp(dataArray->first_attribute("type")->value(), "Int64") == 0)
	{
          this->read_face_neighbors_data(dataArray->value(), elemTypes, mesh.faceMesh);
	}
        else
	{
          std::cerr<<"VtuReader error: wrong type for face neighbors local ids, aborting.";
	  std::abort();
	}
      }
    }
//// Read PointData
//  const rapidxml::xml_node<>* pointData = piece->first_node("PointData");
//  if (pointData)
//  {
//    std::vector<Range> vertex_list(1);
//    vertex_list[0] = vertices;
//    for(rapidxml::xml_node<> *dataArray=pointData->first_node("DataArray"); dataArray; dataArray=dataArray->next_sibling("DataArray"))
//    {
//      long nComponents = 0;
//      if (dataArray->first_attribute("NumberOfComponents"))
//        nComponents = atol(dataArray->first_attribute("NumberOfComponents")->value());
//      else
//        nComponents = 1;
//      if (!dataArray->first_attribute("Name"))
//        std::cerr<<"VtuReader error: missing name of CellData DataArray, aborting.");
//      if (!dataArray->first_attribute("type"))
//        std::cerr<<"VtuReader error: missing type of CellData DataArray "
//                 <<dataArray->first_attribute("Name")->value()<<", aborting.");
//      if (strcmp(dataArray->value(), "\0") == 0)
//        std::cerr<<"VtuReader error: missing data of CellData DataArray "
//                 <<dataArray->first_attribute("Name")->value()<<", aborting.");
//      this->read_tag_data(dataArray->value(), nComponents, vertex_count, vertex_list, 
//                          dataArray->first_attribute("Name")->value(), dataArray->first_attribute("type")->value());
//    }
//  }
}

template<unsigned int dim>
void ReadMeshVtu::vtu_read_unstructured_piece (const rapidxml::xml_node<> *& piece, 
                                                    Mesh<dim>& mesh,
						    std::vector<ElemTypesMesh>& elemTypes)
{
  if (!piece->first_attribute("NumberOfPoints")) 
  {
    std::cerr<<"VtuReader error: missing number of points, aborting.";
    std::abort();
  }
  long numberOfPoints = atol(piece->first_attribute("NumberOfPoints")->value());

  const rapidxml::xml_node<> * pointCoords = piece->first_node("Points")->first_node("DataArray");
  if (!pointCoords)
  {
    std::cerr<<"VtuReader error: missing point coordinates node, aborting.";
    std::abort();
  }
  if (!pointCoords->first_attribute("NumberOfComponents"))
  {
    std::cerr<<"VtuReader error: missing number of Points coordinates, aborting.";
    std::abort();
  }
  long numberOfComponents = atol(pointCoords->first_attribute("NumberOfComponents")->value());
  if (strcmp(pointCoords->first_attribute("NumberOfComponents")->value(),"3") != 0)
    std::cerr<<"VtuReader error: reading points with "
             <<pointCoords->first_attribute("NumberOfComponents")->value()
             <<" coordinates. 3 coordinates required!";
  if (!pointCoords->first_attribute("type"))
  {
    std::cerr<<"VtuReader error: missing coordinates data type, aborting.";
    std::abort();
  }
  const char * dataType = pointCoords->first_attribute("type")->value();
  if (!VtuMeshUtil::check_vtu_data_type(dataType,"Float")) 
    std::cerr<<"VtuReader error: wrong coordinates data type "<<dataType;
  if (!pointCoords->first_attribute("format"))
  {
    std::cerr<<"VtuReader error: missing coordinates data type, aborting.";
    std::abort();
  }
  if (strcmp(pointCoords->first_attribute("format")->value(),"ascii") != 0)
  {
    std::cerr<<"VtuReader error: data should be in ascii format. Format is "
             <<pointCoords->first_attribute("format")->value();
    std::abort();
  }
    
  // Create vertices and read coordinates
  const char * coordsString = pointCoords->value();
  if (strcmp(coordsString, "\0") == 0)
  {
    std::cerr<<"VtuReader error: missing value of point coordinates array, aborting.";
    std::abort();
  }
  read_vertices( coordsString, numberOfPoints, numberOfComponents, mesh);
  
  if (!piece->first_attribute("NumberOfCells"))
  {
    std::cerr<<"VtuReader error: missing number of cells, aborting.";
    std::abort();
  }
  // Read element connectivity for all elements in this piece
  long numberOfCells = atol(piece->first_attribute("NumberOfCells")->value());
  rapidxml::xml_node<> * cells = piece->first_node("Cells");
  if (!cells)
  {
    std::cerr<<"VtuReader error: missing cells node, aborting.";
    std::abort();
  }
  std::map<std::string, const char *> cellData;
  std::map<std::string, const char *> cellDataTypes;
  for(rapidxml::xml_node<> *dataArray=cells->first_node("DataArray"); dataArray; dataArray=dataArray->next_sibling("DataArray"))
  {
    if (!dataArray->first_attribute("Name"))
    {
      std::cerr<<"VtuReader error: missing name of cells DataArray, aborting.";
      std::abort();
    }
    if (strcmp(dataArray->value(), "\0") == 0)
    {
      std::cerr<<"VtuReader error: missing data of cells DataArray "
	       <<dataArray->first_attribute("Name")->value()<<", aborting.";
      std::abort();
    }
    cellData[std::string(dataArray->first_attribute("Name")->value())] 
                                    = dataArray->value();
    if (!dataArray->first_attribute("type"))
    {
      std::cerr<<"VtuReader error: missing type of cells DataArray "
	       <<dataArray->first_attribute("Name")->value()<<", aborting.";
      std::abort();
    }
    cellDataTypes[std::string(dataArray->first_attribute("Name")->value())] 
                                    = dataArray->first_attribute("type")->value();
    if (!dataArray->first_attribute("format"))
    {
      std::cerr<<"VtuReader error: missing format of cells DataArray "
	       <<dataArray->first_attribute("Name")->value()<<", aborting.";
      std::abort();
    }
    if (strcmp(dataArray->first_attribute("format")->value(),"ascii") != 0)
    {
      std::cerr<<"VtuReader error: Cell data should be in ascii format. Format is "
               <<dataArray->first_attribute("format")->value();
      std::abort();
    }
  }
  if (cellData.find("connectivity") == cellData.end() || 
      cellData.find("offsets") == cellData.end() || 
      cellData.find("types") == cellData.end())
  {
    std::cerr<<"VtuReader error: Cell data missing.";
    std::abort();
  }
  if (!VtuMeshUtil::check_vtu_data_type(cellDataTypes["connectivity"],"Int") || 
      !VtuMeshUtil::check_vtu_data_type(cellDataTypes["offsets"],"Int") || 
      !VtuMeshUtil::check_vtu_data_type(cellDataTypes["types"],"UInt"))
  {
    std::cerr<<"VtuReader error: wrong data type in Cells section";
    std::abort();
  }
  read_elements( cellData, numberOfCells, mesh, elemTypes);
}

template <unsigned int dim>
void ReadMeshVtu::read_elements( std::map<std::string, const char *>& cellData,
                                      const long numberOfCells,
				      Mesh<dim>& mesh, std::vector<ElemTypesMesh>& elemTypes)
{
  const char * offsets = cellData["offsets"];
  const char * connectivity = cellData["connectivity"];
  const char * types = cellData["types"];
  if (!VtuMeshUtil::checkStringValuesSize(offsets,numberOfCells))
  {
    std::cerr<<"VtuReader error: cell offsets array is too short";
    std::abort();
  }
  if (!VtuMeshUtil::checkStringValuesSize(types,numberOfCells))
  {
    std::cerr<<"VtuReader error: cell types array is too short";
    std::abort();
  }
  if (!VtuMeshUtil::checkConnectivityStringSize(connectivity,offsets))
  {
    std::cerr<<"VtuReader error: connectivity array is too short";
    std::abort();
  }
  long connectivityStringOffset = 0;
  long offsetStringOffset = 0;
  long typeStringOffset = 0;
  char * endOffset;
  char * endConnectivity;
  char * endType;
  // Create elements in blocks of the same type
  // Preserve the order in which the elements 
  // were read for later reading attribute data.
  long offset_bkp = 0;
  long totElems = 0;
  while (*(types+typeStringOffset) != '\0')
  {
    long vertSum = 0;
    long blockElems = 0;
    unsigned vtk_type_bkp = strtoul (types+typeStringOffset, NULL, 10);
    unsigned vtk_type = vtk_type_bkp;
    long offset = strtol (offsets + offsetStringOffset, NULL, 10);
    long nVerts = offset-offset_bkp;
    for(;;)
    {
      vtk_type = strtoul (types + typeStringOffset, &endType, 10);
      if (vtk_type == vtk_type_bkp && *(types+typeStringOffset) != '\0') 
      {
        offset = strtol (offsets + offsetStringOffset, &endOffset, 10);
        typeStringOffset = endType - types;
        offsetStringOffset = endOffset - offsets;
        assert (offsetStringOffset <= strlen(offsets));
        assert (typeStringOffset <= strlen(types));
	assert (nVerts == (offset-offset_bkp)/(blockElems+1));
	vertSum += nVerts;
        blockElems++;
        totElems++;
        assert (totElems <= numberOfCells);
      }
      else 
	break;
    }
    long blockConnSize = offset-offset_bkp;
    offset_bkp = offset;
    
    const ElemTypesMesh type = VtuMeshUtil::vtuElemTypes[vtk_type_bkp].spa_type;
    std::size_t firstColIndex;
    mesh.resizeElemNodesConnConservative(type, blockElems, firstColIndex);
    ElemNodesConnMatrixDynamic& connMatrix = mesh.getElemNodesConnForType(type);
    long blockConn = 0;
    for (long e = 0; e < blockElems; e++)
    {
      elemTypes.push_back(type);
      for (long i = 0; i < nVerts; i++)
      {
        connMatrix(i,e+firstColIndex) =  strtoul (connectivity+connectivityStringOffset, &endConnectivity, 10);
        connectivityStringOffset = endConnectivity - connectivity;
        assert (connectivityStringOffset<=strlen(connectivity));
        blockConn++;
        assert(blockConn <= blockConnSize);
      }
      const unsigned* order = VtuMeshUtil::vtuElemTypes[vtk_type_bkp].node_order;
      if ( order )
      {
	Connectivity::TopoIndexVector tmp_conn_list = connMatrix.col(e+firstColIndex);
        for (int j = 0; j < nVerts; ++j)
          connMatrix(order[j],e+firstColIndex) = tmp_conn_list[j]; 
      }
    }
    assert(blockConn == blockConnSize);
  }
  assert (totElems == numberOfCells);
  assert (*(connectivity+connectivityStringOffset) == '\0');
  assert (*(offsets+offsetStringOffset) == '\0');
  assert (*(types+typeStringOffset) == '\0');
}

template<unsigned int dim>
void ReadMeshVtu::read_vertices( const char *& coordsString,
                                  long numberOfPoints,
				  long numberOfComponents,
				  Mesh<dim>& mesh)
{
  long totalCoords = numberOfPoints * numberOfComponents;
  assert(numberOfComponents == 3);

  if (!VtuMeshUtil::checkStringValuesSize(coordsString,numberOfPoints*numberOfComponents))
  {
    std::cerr<<"VtuReader error: points coordinates array is too short";
    std::abort();
  }

  mesh.nodesCoords.resize(Eigen::NoChange,numberOfPoints);

  char * endToken;
  long offset = 0;
  for (long vtx = 0; vtx < numberOfPoints; ++vtx)
    for (unsigned int d = 0; d < 3; d++)
    {
      mesh.nodesCoords(d,vtx) = strtod (coordsString+offset, &endToken);
      offset = endToken - coordsString;
      assert (offset>0);
      assert (offset<=strlen(coordsString));
    }
  
}

template<unsigned int dim>
void ReadMeshVtu::read_element_faces_data(const char * dataArrayString, 
                                               long numberOfComponents,
			                       const std::vector<ElemTypesMesh>& elemTypes,
  		                               ElemMesh<dim>& mesh)
{
  long numberOfValues = elemTypes.size();
  if (!VtuMeshUtil::checkStringValuesSize(dataArrayString,numberOfValues*numberOfComponents))
  {
    std::cerr<<"VtuReader error: data array is too short"<<std::endl;
    std::abort();
  }
  size_t size = 0;

  std::cout<<"Reading element faces"<<std::endl;
  char * endToken;
  long offset = 0;
  ElemTypesMesh elemTypeMesh_bkp = elemTypes[0];
  long blockElems = 0;
  for (std::vector<ElemTypesMesh>::const_iterator iter = elemTypes.begin(); ; iter++)
  {
    if (iter != elemTypes.end() && *iter == elemTypeMesh_bkp)
      blockElems++;
    else
    {
      unsigned int elemDim = Connectivity::dimensionFromType(elemTypeMesh_bkp);
      if (elemDim != dim)
      {
        for (long i = 0; i < numberOfComponents*blockElems; i++)
        {
          int elementFace = strtol (dataArrayString+offset, &endToken, 10);
          offset = endToken - dataArrayString;
          assert (offset<=strlen(dataArrayString));
	  if (elementFace != -1)
            std::cout<<"Error, VtuReader found a face for an element of dimension "<<elemDim<<std::endl;
        }
      }
      else 
      {
        typename MeshTypes<dim>::Types elemMeshType = static_cast<typename MeshTypes<dim>::Types>(mesh.dimMeshElemType(elemTypeMesh_bkp));
        std::size_t firstColIndex;
        mesh.resizeElemSidesConnConservative(elemMeshType, blockElems, firstColIndex);
	int numberOfCompTrue = Connectivity::elemTypesSideNumberAll[elemTypeMesh_bkp];
        ElemSidesConnMatrixDynamic& connMatrix = mesh.getElemSidesConnForType(elemMeshType);
        for (long e = 0; e < blockElems; e++)
	{
          for (long i = 0; i < numberOfCompTrue; i++)
	  {
            connMatrix(i,e+firstColIndex) = strtol (dataArrayString+offset, &endToken, 10);
            offset = endToken - dataArrayString;
            assert (offset<=strlen(dataArrayString));
	  }
	  for (long i = numberOfCompTrue; i < numberOfComponents; i++)
	  {
            int elementFace = strtol (dataArrayString+offset, &endToken, 10);
            offset = endToken - dataArrayString;
            assert (offset<=strlen(dataArrayString));
	    if (elementFace != -1)
	    {
              std::cout<<"Error, VtuReader found "<<i+1 <<" faces for an element of type "<<elemTypeMesh_bkp<<std::endl;
	      std::cout<<"Element number is "<<e<<" among "<<blockElems<<" elements of type "<<elemTypeMesh_bkp<<std::endl;
	    }
	  }
	}
      }
      if (iter == elemTypes.end())
	break;
      elemTypeMesh_bkp = *iter;
      blockElems = 1;
    }
  }
  assert (*(dataArrayString+offset) == '\0');
}

template<unsigned int dim>
void ReadMeshVtu::read_face_neighbors_data(const char * dataArrayString, 
			                        const std::vector<ElemTypesMesh>& elemTypes,
  		                                FaceMesh<dim>& mesh)
{
  long numberOfValues = elemTypes.size();
  long numberOfComponents = 2;
  if (!VtuMeshUtil::checkStringValuesSize(dataArrayString,numberOfValues*2))
    std::cerr<<"VtuReader error: data array is too short"<<std::endl;
  size_t size = 0;

  std::cout<<"Reading face neighbors"<<std::endl;
  char * endToken;
  long offset = 0;
  ElemTypesMesh elemTypeMesh_bkp = elemTypes[0];
  long blockElems = 0;
  for (std::vector<ElemTypesMesh>::const_iterator iter = elemTypes.begin(); ; iter++)
  {
    if (iter != elemTypes.end() && *iter == elemTypeMesh_bkp)
      blockElems++;
    else
    {
      unsigned int elemDim = Connectivity::dimensionFromType(elemTypeMesh_bkp);
      if (elemDim != dim)
      {
        for (long i = 0; i < numberOfComponents*blockElems; i++)
        {
          int elementFace = strtol (dataArrayString+offset, &endToken, 10);
          offset = endToken - dataArrayString;
          assert (offset<=strlen(dataArrayString));
	  if (elementFace != -1)
            std::cout<<"Error, VtuReader found a neighbor for an element of dimension "<<elemDim<<std::endl;
        }
      }
      else 
      {
        typename MeshTypes<dim>::Types elemMeshType = static_cast<typename MeshTypes<dim>::Types>(mesh.dimMeshElemType(elemTypeMesh_bkp));
        std::size_t firstColIndex;
        mesh.resizeSideElemsConnConservative(elemMeshType, blockElems, firstColIndex);
        SideElemsConnMatrix& connMatrix = mesh.getSideElemsConnForType(elemMeshType);
        for (long e = 0; e < blockElems; e++)
	{
          for (long i = 0; i < 2; i++)
	  {
            connMatrix(i,e+firstColIndex) = strtol (dataArrayString+offset, &endToken, 10);
            offset = endToken - dataArrayString;
            assert (offset<=strlen(dataArrayString));
	  }
	}
      }
      if (iter == elemTypes.end())
	break;
      elemTypeMesh_bkp = *iter;
      blockElems = 1;
    }
  }
  assert (*(dataArrayString+offset) == '\0');
}
#endif




