//  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: WriterMedit.cpp,v 1.4 2007/05/20 23:15:36 delpinux Exp $

#include <WriterMedit.hpp>
#include <Mesh.hpp>

#include <Structured3DMesh.hpp>
#include <SpectralMesh.hpp>
#include <MeshOfHexahedra.hpp>
#include <MeshOfTetrahedra.hpp>

#include <SurfaceMeshOfTriangles.hpp>
#include <SurfaceMeshOfQuadrangles.hpp>
#include <MeshOfTriangles.hpp>

#include <fstream>

#include <ScalarFunctionBase.hpp>
#include <FieldOfScalarFunction.hpp>

#include <FEMFunction.hpp>

void WriterMedit::
__fillCrossedComponent(const FieldOfScalarFunction& field,
		       Vector<real_t>& values) const
{
  const size_t numberOfComponents = field.numberOfComponents();
  ASSERT(values.size() == __mesh->numberOfVertices()*numberOfComponents);

  for (size_t i = 0; i<numberOfComponents; ++i) {
    const ScalarFunctionBase& function = *field.function(i);

    switch (function.type()) {
    case ScalarFunctionBase::femfunction: {
      const FEMFunctionBase& fem
	= static_cast<const FEMFunctionBase&>(function);
      if ((fem.baseMesh() == __mesh) and
	  ((fem.discretizationType() == DiscretizationType::lagrangianFEM1))) {
	for (size_t j=0; j<fem.values().size(); ++j) {
	  values[numberOfComponents*j+i] = fem[i];
	}
	break;
	} // if not continues the standard method
      }
    default: {
      for (size_t j=0; j<__mesh->numberOfVertices(); ++j) {
	const TinyVector<3,real_t>& X = __mesh->vertex(j);
	  values[numberOfComponents*j+i] = function(X);
      }
    }
    }
  }
}


void WriterMedit::
__proceedData() const
{
  std::string fileName = __fileName;
  fileName += ".bb";
  std::ofstream file(fileName.c_str());
  if (not(file)) {
    throw ErrorHandler(__FILE__,__LINE__,
		       "cannot open file '"+stringify(fileName)+"'",
		       ErrorHandler::normal);
  }

  if (__fieldList.size() + __scalarFunctionList.size() > 1) {
    throw ErrorHandler(__FILE__,__LINE__,
		       "writing file '"+stringify(__fileName.c_str())+
		       "': cannot save more than one field or function in Medit format!",
		       ErrorHandler::normal);
  }

  if (__scalarFunctionList.size()>0) {
    file << "3 1 " << __mesh->numberOfVertices() << " 2" << __CR;

    const ScalarFunctionBase& f = *__scalarFunctionList[0];

    switch (f.type()) {
    case ScalarFunctionBase::femfunction: {
      const FEMFunctionBase& fem = static_cast<const FEMFunctionBase&>(f);
      if ((fem.baseMesh() == __mesh) and
	  (fem.discretizationType() == DiscretizationType::lagrangianFEM1)) {
	for (size_t i=0; i<fem.values().size(); ++i) {
	  file << fem[i] << __CR;
	}
	break;
      } // if not continues the standard method
    }
    default: {
      for (size_t i=0; i<__mesh->numberOfVertices(); ++i) {
	const TinyVector<3,real_t>& X = __mesh->vertex(i);
	file << f(X) << __CR;
      }
    }
    }
  }

  if (__fieldList.size() > 0) {
    const FieldOfScalarFunction& field = *__fieldList[0];
    file << "3 " << field.numberOfComponents() << " " << __mesh->numberOfVertices() << " 2" << __CR;

    Vector<real_t> values(3*__mesh->numberOfVertices());

    this->__fillCrossedComponent(field, values);

    for (size_t i=0; i<values.size(); ++i) {
      file << values[i] << __CR;
    }
  }
}

template <typename MeshType>
void WriterMedit::
__saveElements(std::ostream& os,
	       const MeshType& mesh) const
{
  typedef typename MeshType::CellType CellType;

  os << mesh.numberOfCells() << __CR;

  for (size_t i = 0; i<mesh.numberOfCells(); ++i) {
    const CellType& c = mesh.cell(i);

    for (size_t j=0; j<CellType::NumberOfVertices; ++j) {
      os << mesh.vertexNumber(c(j))+1 << ' ';
    }

    os << c.reference() << __CR;
  }
}

void WriterMedit::
__proceedMesh() const
{
  std::string fileName = __fileName;
  fileName += ".mesh";
  std::ofstream file(fileName.c_str());
  if (not(file)) {
    throw ErrorHandler(__FILE__,__LINE__,
		       "cannot open file '"+stringify(fileName)+"'",
		       ErrorHandler::normal);
  }

  file << "MeshVersionFormatted 1" << __CR;
  file << "Dimension" << __CR;
  file << 3 << __CR;

  file << "Vertices " << __mesh->numberOfVertices() << __CR;

  for (size_t i=0; i<__mesh->numberOfVertices(); ++i) {
    const TinyVector<3>& X = __mesh->vertex(i);
    file << X[0] << ' ' <<  X[1] << ' ' << X[2]
	 << ' ' << __mesh->vertex(i).reference() << __CR;
  }

  switch (__mesh->type()) {
  case Mesh::cartesianHexahedraMesh: {
    const Structured3DMesh& mesh = dynamic_cast<const Structured3DMesh&>(*__mesh);
    file << "Hexahedra ";
    this->__saveElements(file, mesh);
    if (mesh.hasSurfaceMesh()) {
      file << "Quadrilaterals ";
      this->__saveElements(file, *mesh.surfaceMesh());
    }
    break;
  }
  case Mesh::spectralMesh: {
    const SpectralMesh& mesh = dynamic_cast<const SpectralMesh&>(*__mesh);
    file << "Hexahedra ";
    this->__saveElements(file, mesh);
    if (mesh.hasSurfaceMesh()) {
      file << "Quadrilaterals ";
      this->__saveElements(file, *mesh.surfaceMesh());
    }
    break;
  }
  case Mesh::hexahedraMesh: {
    const MeshOfHexahedra& mesh = dynamic_cast<const MeshOfHexahedra&>(*__mesh);
    file << "Hexahedra ";
    this->__saveElements(file, mesh);
    if (mesh.hasSurfaceMesh()) {
      file << "Quadrilaterals ";
      this->__saveElements(file, *mesh.surfaceMesh());
    }
    break;
  }
  case Mesh::tetrahedraMesh: {
    const MeshOfTetrahedra& mesh = dynamic_cast<const MeshOfTetrahedra&>(*__mesh);
    file << "Tetrahedra ";
    this->__saveElements(file, mesh);
    if (mesh.hasSurfaceMesh()) {
      file << "Triangles ";
      this->__saveElements<SurfaceMeshOfTriangles>(file, *mesh.surfaceMesh());
    }
    break;
  }
  case Mesh::trianglesMesh: {
    const MeshOfTriangles& mesh = dynamic_cast<const MeshOfTriangles&>(*__mesh);
    file << "Triangles ";
    this->__saveElements(file, mesh);
    break;
  }
  case Mesh::surfaceMeshTriangles: {
    const SurfaceMeshOfTriangles& mesh = dynamic_cast<const SurfaceMeshOfTriangles&>(*__mesh);
    file << "Triangles ";
    this->__saveElements(file, mesh);
    break;
  }
  case Mesh::surfaceMeshQuadrangles: {
    const SurfaceMeshOfQuadrangles& mesh = dynamic_cast<const SurfaceMeshOfQuadrangles&>(*__mesh);
    file << "Quadrilaterals ";
    this->__saveElements(file, mesh);
    break;
  }
  default: {
    throw ErrorHandler(__FILE__,__LINE__,
		       "unexpected mesh type",
		       ErrorHandler::unexpected);
  }
  }

  file << "End" << __CR;
}

void WriterMedit::
proceed() const
{
  if (__fieldList.size() + __scalarFunctionList.size() == 0) {
    this->__proceedMesh();
  } else {
    this->__proceedData();
  }
}


WriterMedit::
WriterMedit(ConstReferenceCounting<Mesh> mesh,
	    const std::string& fileName,
	    const std::string& CR)
  : WriterBase(mesh,fileName,CR)
{
  ;
}

WriterMedit::
~WriterMedit()
{
  ;
}