//  This file is part of ff3d - http://www.freefem.org/ff3d
//  Copyright (C) 2001, 2002, 2003 Stphane Del Pino

//  This program is free software; you can redistribute it and/or modify
//  it under the terms of the GNU General Public License as published by
//  the Free Software Foundation; either version 2, 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 General Public License for more details.

//  You should have received a copy of the GNU General Public License
//  along with this program; if not, write to the Free Software Foundation,
//  Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.  

//  $Id: GmshFormatReader.cpp,v 1.6 2004/12/31 16:38:58 delpinux Exp $

#include <GmshFormatReader.hpp>

#include <SurfaceMeshOfTriangles.hpp>
#include <MeshOfTetrahedra.hpp>

void
GmshFormatReader::__readVertices()
{
  const int numberOfVerices = this->__getInteger();
  ffout(3) << "- Number Of Vertices: " << numberOfVerices << '\n';
  if (numberOfVerices<0) {
    throw ErrorHandler(__FILE__,__LINE__,
		       "reading file '"+this->__fileName
		       +"': number of vertices is negative",
		       ErrorHandler::normal);
  }

  __verticesNumbers.resize(numberOfVerices);
  __vertices = new VerticesSet(numberOfVerices);
  VerticesSet& vertices = *__vertices;

  for (int i=0; i<numberOfVerices; ++i) {
    __verticesNumbers[i] = this->__getInteger();
    const real_t x = this->__getReal();
    const real_t y = this->__getReal();
    const real_t z = this->__getReal();
    vertices[i] = TinyVector<3, real_t>(x,y,z);
  }
}

void
GmshFormatReader::__readElements()
{
  const int numberOfElements = this->__getInteger();
  ffout(3) << "- Number Of Elements: " << numberOfElements << '\n';
  if (numberOfElements<0) {
    throw ErrorHandler(__FILE__,__LINE__,
		       "reading file '"+__fileName
		       +"': number of elements is negative",
		       ErrorHandler::normal);
  }

  __elementType.resize(numberOfElements);
  __references.resize(numberOfElements);
  __elementVertices.resize(numberOfElements);

  for (int i=0; i<numberOfElements; ++i) {
    this->__getInteger(); // drops element number
    const int elementType = this->__getInteger()-1;

    if ((elementType < 0) or (elementType > 14)) {
      throw ErrorHandler(__FILE__,__LINE__,
			 "reading file '"+__fileName
			 +"': unknown element type '"+stringify(elementType)+"'",
			 ErrorHandler::normal);
    }
    __elementType[i] = elementType;
    __references[i] = this->__getInteger(); // physical reference
    this->__getInteger(); // drops entity reference

    const int numberOfVertices = this->__getInteger();
    if (numberOfVertices != __numberOfPrimitiveNodes[elementType]) {
      std::stringstream errorMsg;
      errorMsg << "reading file '" <<__fileName
	       << "':number of vertices (" << numberOfVertices
	       << ") does not match expected ("
	       << __numberOfPrimitiveNodes[elementType] << ")" << std::ends;

      throw ErrorHandler(__FILE__,__LINE__,
			 errorMsg.str(),
			 ErrorHandler::normal);
    }
    __elementVertices[i].resize(numberOfVertices);
    for (int j=0; j<numberOfVertices; ++j) {
      __elementVertices[i][j] = this->__getInteger();
    }
  }  
}

void
GmshFormatReader::__proceedData()
{
  TinyVector<15,int> elementNumber = 0;
  std::vector<std::set<size_t> > elementReferences(15);

  for (size_t i=0; i<__elementType.size(); ++i) {
    const int elementType = __elementType[i];
    elementNumber[elementType]++;
    elementReferences[elementType].insert(__references[i]);
  }

  for (size_t i=0; i<elementNumber.size(); ++i) {
    if (elementNumber[i]>0) {
      ffout(3) << "  - Number of "
	       << __primitivesNames[i]
	       << ": " << elementNumber[i];
      if (not(__supportedPrimitives[i])) {
	ffout(3) << " [IGNORED]";
      } else {
	ffout(3) << " references: ";
	for(std::set<size_t>::const_iterator iref
	      = elementReferences[i].begin();
	    iref != elementReferences[i].end(); ++iref) {
	  if (iref != elementReferences[i].begin()) ffout(3) << ',';
	  ffout(3) << *iref;
	}

	switch (i) {
	  // Supported entities
	case 1: { // triangles
	  __triangles = new Vector<Triangle>(elementNumber[i]);
	  break;
	}
	case 3: { // tetrahedra
	  __tetrahedra = new Vector<Tetrahedron>(elementNumber[i]);
	  break;
	}
	  // Unsupported entities
	case  0: // edge
	case  2: // quadrangle
	case  4: // hexaredron
	case  5: // prism
	case  6: // pyramid
	case  7: // second order edge
	case  8: // second order triangle
	case  9: // second order quadrangle
	case 10: // second order tetrahedron
	case 11: // second order hexahedron
	case 12: // second order prism
	case 13: // second order pyramid
	case 14:{// point
	  break; 
	}
	default: {
	  throw ErrorHandler(__FILE__,__LINE__,
			     "reading file '"+__fileName
			     +"': ff3d cannot read this kind of element",
			     ErrorHandler::normal);
	}
	}
      }
      ffout(3) << '\n';
    }
  }

  ffout(3) << "- Building correspondance table\n";
  int minNumber = std::numeric_limits<int>::max();
  int maxNumber = std::numeric_limits<int>::min();
  for (size_t i=0; i<__verticesNumbers.size(); ++i) {
    const int& vertexNumber = __verticesNumbers[i];
    minNumber = std::min(minNumber,vertexNumber);
    maxNumber = std::max(maxNumber,vertexNumber);
  }

  if (minNumber<0) {
    throw ErrorHandler(__FILE__,__LINE__,
		       "reading file '"+__fileName
		       +"': ff3d does not implement negative vertices number",
		       ErrorHandler::normal);
  }

 // A value of -1 means that the vertex is unknown
  __verticesCorrepondance.resize(maxNumber+1,-1);

  for (size_t i=0; i<__verticesNumbers.size(); ++i) {
    __verticesCorrepondance[__verticesNumbers[i]] = i;
  }

  // reset element number to count them while filling structures
  elementNumber = 0;  

  VerticesSet& vertices = *__vertices;

  // Element structures filling
  for (size_t i=0; i<__elementType.size(); ++i) {
    std::vector<int>& elementVertices = __elementVertices[i];
    switch (__elementType[i]) {
      // Supported entities
    case 1: { // triangles
      int& triangleNumber = elementNumber[1];

      const size_t a = __verticesCorrepondance[elementVertices[0]];
      const size_t b = __verticesCorrepondance[elementVertices[1]];
      const size_t c = __verticesCorrepondance[elementVertices[2]];
      (*__triangles)[triangleNumber]
	= Triangle(vertices[a], vertices[b], vertices[c],
		   __references[i]);
      triangleNumber++; // one more triangle
      break;
    }
    case 3: { // tetrahedra
      int& tetrahedronNumber = elementNumber[3];

      const int a = __verticesCorrepondance[elementVertices[0]];
      const int b = __verticesCorrepondance[elementVertices[1]];
      const int c = __verticesCorrepondance[elementVertices[2]];
      const int d = __verticesCorrepondance[elementVertices[3]];
      if ((a<0) or (b<0) or (c<0) or (d<0)) {
	throw ErrorHandler(__FILE__,__LINE__,
			   "reading file '"+__fileName
			   +"': error reading element "+stringify(i)
			   +" [bad vertices definition]",
			   ErrorHandler::normal);
      }
      (*__tetrahedra)[tetrahedronNumber]
	= Tetrahedron(vertices[a], vertices[b], vertices[c], vertices[d],
		      __references[i]);
      tetrahedronNumber++; // one more triangle
      break;
    }
      // Unsupported entities
    case  0: // edge
    case  2: // quadrangle
    case  4: // hexaredron
    case  5: // prism
    case  6: // pyramid
    case  7: // second order edge
    case  8: // second order triangle
    case  9: // second order quadrangle
    case 10: // second order tetrahedron
    case 11: // second order hexahedron
    case 12: // second order prism
    case 13: // second order pyramid
    case 14:{// point
      break; 
    }
    default: {
	  throw ErrorHandler(__FILE__,__LINE__,
			     "reading file '"+__fileName
			     +"': ff3d cannot read this kind of element",
			     ErrorHandler::normal);
    }
    }
  }
}

GmshFormatReader::Keyword
GmshFormatReader::__nextKeyword()
{
  GmshFormatReader::Keyword kw;

  char pKeyword[256];
  int retval = fscanf(__ifh,"%255s",pKeyword);
  const std::string aKeyword = pKeyword;

  if (retval < 0)  {
    kw.second = EndOfFile;
    return kw;
  } else if (retval == 0) {
    kw.second = Unknown;
    return kw;
  }

  KeywordList::iterator i = __keywordList.find(aKeyword.c_str());

  if (i != __keywordList.end()) {
    kw.first = (*i).first;
    kw.second = (*i).second;
    return kw;
  }

  throw ErrorHandler(__FILE__,__LINE__,
		     "reading file '"+__fileName
		     +"': unknown keyword "+aKeyword,
		     ErrorHandler::normal);

  kw.first = aKeyword;
  kw.second = Unknown;
  return kw;
}

GmshFormatReader::GmshFormatReader(const std::string & s)
  : MeshReader(s)
{
  __keywordList["$NOD"]    = NOD;
  __keywordList["$ENDNOD"] = ENDNOD;
  __keywordList["$ELM"]    = ELM;
  __keywordList["$ENDELM"] = ENDELM;

  __numberOfPrimitiveNodes.resize(16);
  __numberOfPrimitiveNodes[ 0] =  2; // edge
  __numberOfPrimitiveNodes[ 1] =  3; // triangle
  __numberOfPrimitiveNodes[ 2] =  4; // quadrangle
  __numberOfPrimitiveNodes[ 3] =  4; // Tetrahedron
  __numberOfPrimitiveNodes[ 4] =  8; // Hexaredron
  __numberOfPrimitiveNodes[ 5] =  6; // Prism
  __numberOfPrimitiveNodes[ 6] =  5; // Pyramid
  __numberOfPrimitiveNodes[ 7] =  3; // second order edge
  __numberOfPrimitiveNodes[ 8] =  6; // second order triangle
  __numberOfPrimitiveNodes[ 9] =  9; // second order quadrangle
  __numberOfPrimitiveNodes[10] = 10; // second order tetrahedron
  __numberOfPrimitiveNodes[11] = 27; // second order hexahedron
  __numberOfPrimitiveNodes[12] = 18; // second order prism
  __numberOfPrimitiveNodes[13] = 14; // second order pyramid
  __numberOfPrimitiveNodes[14] =  1; // point                    

  __primitivesNames[0]     = "edges";
  __supportedPrimitives[0] = false;
  __primitivesNames[1]     = "triangles";
  __supportedPrimitives[1] = true;
  __primitivesNames[2]     = "quadrangles";
  __supportedPrimitives[2] = true;
  __primitivesNames[3]     = "tetrahedra";
  __supportedPrimitives[3] = true;
  __primitivesNames[4]     = "hexahedra";
  __supportedPrimitives[4] = true;
  __primitivesNames[5]     = "prisms";
  __supportedPrimitives[5] = false;
  __primitivesNames[6]     = "pyramids";
  __supportedPrimitives[6] = false;
  __primitivesNames[7]     = "second order edges";
  __supportedPrimitives[7] = false;
  __primitivesNames[8]     = "second order triangles";
  __supportedPrimitives[8] = false;
  __primitivesNames[9]     = "second order quadrangles";
  __supportedPrimitives[9] = false;
  __primitivesNames[10]    = "second order tetrahedra";
  __supportedPrimitives[10]= false;
  __primitivesNames[11]    = "second order hexahedra";
  __supportedPrimitives[11]= false;
  __primitivesNames[12]    = "second order prisms";
  __supportedPrimitives[12]= false;
  __primitivesNames[13]    = "second order pyramids";
  __supportedPrimitives[13]= false;
  __primitivesNames[14]    = "point";
  __supportedPrimitives[14]= false;

  ffout(3) << "\n";
  ffout(3) << "Reading file " << s << '\n';

  // Getting vertices list
  GmshFormatReader::Keyword kw = this->__nextKeyword();
  if (kw.second != NOD) {
    throw ErrorHandler(__FILE__,__LINE__,
		       "reading file '"+__fileName
		       +"': expecting $NOD, '"+kw.first+"' was found",
		       ErrorHandler::normal);
  }
  
  this->__readVertices();

  kw = this->__nextKeyword();
  if (kw.second != ENDNOD) {
    throw ErrorHandler(__FILE__,__LINE__,
		       "reading file '"+__fileName
		       +"': expecting $ENDNOD, '"+kw.first+"' was found",
		       ErrorHandler::normal);
  }

  // Getting elements list
  kw = this->__nextKeyword();
  if (kw.second != ELM) {
    throw ErrorHandler(__FILE__,__LINE__,
		       "reading file '"+__fileName
		       +"': expecting $ELM, '"+kw.first+"' was found",
		       ErrorHandler::normal);
  }

  this->__readElements();

  kw = this->__nextKeyword();
  if (kw.second != ENDELM) {
    throw ErrorHandler(__FILE__,__LINE__,
		       "reading file '"+__fileName
		       +"': expecting $ENDELM, '"+kw.first+"' was found",
		       ErrorHandler::normal);
  }

  this->__proceedData();

  this->__createMesh();
}