// 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/. */

#ifndef __vtk_writer_spamesh_h
#define __vtk_writer_spamesh_h

#include "DofMap.hpp"
#include "FELib.hpp"
#include "BasisFunction.hpp"
#include "CollectiveBasisFunction.hpp"
#include "CollectiveMappingFunction.hpp"
#include "VTKHelperSpaMesh.hpp"
#include "ReferenceSpaceFiniteElement.hpp"
#include "MappingBasisFunction.hpp"

template<unsigned int Dim, typename FEManagerT>
class VTKWriterSpaMesh 
{
public:
  typedef FEManagerT FEManagerType;
  typedef typename FEManagerType::FE FiniteElementSpaceType;
  typedef typename FEManagerType::MappingFunctionType MappingFunctionType;
  typedef typename FEManagerType::BasisFunctionType BasisFunctionType;
  typedef SpaMesh::Connectivity::TopoPointVector TopoPointVector;
  
  VTKWriterSpaMesh() : 
    _totNumberOfNodesToPlot(0),
    _totNumberOfElementsToPlot(0),
    _connectivitySize(0),
    _initialized(false)
  {
    ;
  }

  ~VTKWriterSpaMesh() {}

  template<typename TopologyType, unsigned int dim> void init(const SpaMesh::Mesh<dim>& mesh, TopologyType& feDataAccessor, const DofMap& dofMap)
  {
    this->initSizes(mesh,feDataAccessor,dofMap);
    this->buildShapeMatrix(mesh,feDataAccessor,dofMap);
    this->buildConnectivity(mesh,feDataAccessor,dofMap);
  }

  template<typename TopologyType, typename Derived, unsigned int dim> 
    void writeOnce(const SpaMesh::Mesh<dim>& mesh, TopologyType& feDataAccessor, const DofMap& dofMap, const Eigen::MatrixBase<Derived>& local_solution)
  {
    this->initSizes(mesh,feDataAccessor,dofMap);
    this->build_discontinuous_solution_vector(mesh,feDataAccessor, dofMap, local_solution);
    this->buildConnectivity(mesh,feDataAccessor,dofMap);
    this->write(mesh,dofMap);
  }
  
  template<typename TopologyType, typename Derived, unsigned int dim> 
    void writeStored(const SpaMesh::Mesh<dim>& mesh, TopologyType& feDataAccessor, const DofMap& dofMap, const Eigen::MatrixBase<Derived>& local_solution)
  {
     this->update_solution_vector(mesh, feDataAccessor, dofMap, local_solution);
     this->write(mesh, dofMap);
  }
  
  template<typename TopologyType, typename Derived, unsigned int dim> void writeOnTheFly(const SpaMesh::Mesh<dim>& mesh,
                                                                                         TopologyType& feDataAccessor, 
                                                                                         const DofMap& dofMap, 
				 				                         const Eigen::MatrixBase<Derived>& local_solution,
								                         const unsigned int step = 0,
								                         const double time = 0)
  {
    const unsigned int n_vars = dofMap.getNVars();
    const unsigned int p_order = dofMap.getPOrder();
    const typename DofMap::VarNamesType var_names = dofMap.getVarNames();
    const typename DofMap::VarScalingsType var_scalings = dofMap.getVarScalings();
    Quadrature::PointVector referenceSpaceNodes;
    ElemType elemTypeOld = NODEELEM;
    BasisFunctionType basisFunction;
    MappingFunctionType mappingFunction;
    typename BasisFunctionType::ShapeMatrix& shapeMatrix = basisFunction.getShapeFunction();
    basisFunction.doComputeShapeDeriv(false);
    Quadrature::PointVector& nodes = mappingFunction.getQpointX();
    this->initSizes(mesh,feDataAccessor,dofMap);

    //write pvtu header
    if (FELib::processor_id() == 0)
    {
      char filename[1024];
      sprintf(filename,"%s%d%s","output_",step,".pvtu");
      FILE* ppFile=fopen(filename,"w");
      fprintf (ppFile,"%s","<VTKFile type=\"PUnstructuredGrid\"");
      fprintf (ppFile,"%s"," version=\"0.1\"");
      fprintf (ppFile,"%s\n"," byte_order=\"LittleEndian\">");
      fprintf (ppFile,"%s\n","   <PUnstructuredGrid GhostLevel=\"0\">");
      fprintf (ppFile,"%s\n","         <FieldData>");
      fprintf (ppFile,"%s%s%s\n","            <DataArray type=\"Float32\" Name=\"","time","\" NumberOfComponents=\"1\" NumberOfTuples=\"1\" format=\"ascii\">");
      fprintf (ppFile,"%f\n",              time);
      fprintf (ppFile,"%s\n","            </DataArray>");
      fprintf (ppFile,"%s\n","         </FieldData>");
      fprintf (ppFile,"%s\n","         <PPoints>");
      fprintf (ppFile,"%s%s%s\n","            <PDataArray type=\"Float32\" NumberOfComponents=\"","3","\">");
      fprintf (ppFile,"%s\n","            </PDataArray>");
      fprintf (ppFile,"%s\n","         </PPoints>");
      fprintf (ppFile,"%s\n","         <PPointData>");
      for(unsigned int v = 0; v < n_vars; v++)
      {
        fprintf (ppFile,"%s%s%s\n","            <PDataArray type=\"Float32\" Name=\"",var_names(v).c_str(),"\" NumberOfComponents=\"1\">");
        fprintf (ppFile,"%s\n","            </PDataArray>");
      }
      fprintf (ppFile,"%s\n","         </PPointData>");
      for (unsigned int i = 0; i < FELib::n_processors(); i++)
      {
        fprintf (ppFile,"%s%d%s%02d%s\n","      <Piece  Source=\"output_",step,"_",i,".vtu\">");
        fprintf (ppFile,"%s\n","      </Piece>  ");        
      }
      fprintf (ppFile,"%s\n","   </PUnstructuredGrid>");
      fprintf (ppFile,"%s\n","</VTKFile>");
      fclose (ppFile);       
    }

    // write header
    char filename[2048];
    sprintf(filename,"output_%d_%02d.vtu",step,FELib::processor_id());
    FILE* pFile=fopen(filename,"w");
    fprintf (pFile,"%s","<VTKFile type=\"UnstructuredGrid\"");
    fprintf (pFile,"%s"," version=\"0.1\"");
    fprintf (pFile,"%s\n"," byte_order=\"LittleEndian\">");
    fprintf (pFile,"%s\n","   <UnstructuredGrid>");
    fprintf (pFile,"%s","      <Piece  ");
    fprintf (pFile,"%s%d%s","NumberOfPoints=\"",_totNumberOfNodesToPlot,"\"  ");
    fprintf (pFile,"%s%d%s\n","NumberOfCells=\"",_totNumberOfElementsToPlot,"\">");
    
    // write mesh nodes  
    fprintf (pFile,"%s\n","         <Points>");
    fprintf (pFile,"%s%s%s\n","            <DataArray type=\"Float32\" NumberOfComponents=\"","3","\" format=\"ascii\">");
    elemTypeOld = NODEELEM;
    for (int elem = 0; elem < mesh.elemCard(); elem++)
    {
      feDataAccessor.initData(mesh,elem);
      const ElemType elemType = feDataAccessor.getElemType();
      const TopoPointVector& elem_coord = feDataAccessor.getElemCoord();
      if (elemType != elemTypeOld)
      {
        VTKHelperSpaMesh::getSubElementsNodesE(elemType, p_order, referenceSpaceNodes);
        elemTypeOld = elemType;
      }
      FiniteElementSpaceType::FiniteElementMappingType::computePhysicalSpacePoints(mappingFunction,referenceSpaceNodes,elemType,elem_coord,ReferenceFiniteElement::getGeoOrder(elemType));
      for (unsigned int n = 0; n < nodes.rows(); n++)
        for(unsigned int i = 0; i < 3; i++)
          fprintf (pFile,"%f ",nodes(n,i));
    }
    fprintf (pFile,"%s\n","");
    fprintf (pFile,"%s\n","            </DataArray>");
    fprintf (pFile,"%s\n","         </Points>");

    //  write solutions 
    elemTypeOld = NODEELEM;
    fprintf (pFile,"%s\n","         <PointData>");
    for(unsigned int j=0;j<n_vars;++j)
    {
      fprintf (pFile,"%s%s%s\n","            <DataArray type=\"Float32\" Name=\"",var_names(j).c_str(),"\" NumberOfComponents=\"1\" format=\"ascii\">");
      for (int elem = 0; elem < mesh.elemCard(); elem++)
      {
        feDataAccessor.initData(mesh,elem);
        const ElemType elemType = feDataAccessor.getElemType();
        if (elemType != elemTypeOld)
        {
          VTKHelperSpaMesh::getSubElementsNodesE(elemType, p_order, referenceSpaceNodes);
          elemTypeOld = elemType;
        }
        FiniteElementSpaceType::initFENodalPoints(basisFunction,p_order,mappingFunction,ReferenceFiniteElement::getGeoOrder(elemType),feDataAccessor,referenceSpaceNodes);
        for (unsigned int i = 0; i < referenceSpaceNodes.rows(); i++)
	{
          //double sol = (local_solution.row(elem_local_id)).segment(j*shapeMatrix.rows(),shapeMatrix.rows()) * shapeMatrix.col(i);
	  //the new sol computation above is cleaner but does not work for p0
          double sol = local_solution.block(elem,j*shapeMatrix.rows(),1,shapeMatrix.rows()).row(0) * (shapeMatrix.col(i));
	  sol *= var_scalings[j];
          fprintf (pFile,"%f ", sol);
	}
      }
      fprintf (pFile,"%s\n","");
      fprintf (pFile,"%s\n","            </DataArray>");
    }
    fprintf (pFile,"%s\n","         </PointData>");

    // write elements connectivity
    unsigned int numberOfSubNodes;
    unsigned int numberOfSubElementNodes;
    Eigen::Matrix<unsigned int, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> numbering;
    fprintf (pFile,"%s\n","         <Cells>");
    fprintf (pFile,"%s\n","            <DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">");
    ElemTypesMesh itemTypeOld = SpaMesh::node;
    unsigned int nodesCounter = 0;
    for (int elem = 0; elem < mesh.elemCard(); elem++)
    {
      ElemTypesMesh itemType = static_cast<ElemTypesMesh>(mesh.elemMesh.elemTypeFromDimType(mesh.elemMesh.typeFromLID(elem)));
      if (itemType != itemTypeOld)
      {
        VTKHelperSpaMesh::getSubElementsNodesNumbering(itemType,p_order,numbering);
        numberOfSubNodes = VTKHelperSpaMesh::getSubElementsTotalNodesNumber(itemType,p_order,numberOfSubElementNodes);
        itemTypeOld = itemType;
      }
      for (unsigned int i = 0; i < numbering.rows(); i++)
        for (unsigned int j = 0; j < numbering.cols(); j++)
          fprintf (pFile, "%d ", numbering(i,j) + nodesCounter);
      nodesCounter += numberOfSubNodes;
    }
    fprintf (pFile,"%s\n","");
    fprintf (pFile,"%s\n","            </DataArray>");
    
    // write elements offset 
    itemTypeOld = SpaMesh::node;
    unsigned int offset = 0;
    unsigned int numberOfSubElements;
    fprintf (pFile,"%s\n","            <DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">");
    for (int elem = 0; elem < mesh.elemCard(); elem++)
    {
      ElemTypesMesh itemType = static_cast<ElemTypesMesh>(mesh.elemMesh.elemTypeFromDimType(mesh.elemMesh.typeFromLID(elem)));
      if (itemType != itemTypeOld)
      {
        numberOfSubNodes = VTKHelperSpaMesh::getSubElementsTotalNodesNumber(itemType,p_order,numberOfSubElementNodes);
        numberOfSubElements = VTKHelperSpaMesh::getNumberOfSubElements(itemType,p_order);
        itemTypeOld = itemType;
      }
      for (unsigned int i = 0; i < numberOfSubElements; i++)	  
      {
        offset += numberOfSubElementNodes;
        fprintf (pFile, "%d ", offset);
      }
    }
//    for(unsigned int i = 0; i < _totNumberOfElementsToPlot; i++)
    fprintf (pFile,"%s\n","");
    fprintf (pFile,"%s\n","            </DataArray>");

    // write elements types 
    unsigned int vtkCellType;
    itemTypeOld = SpaMesh::node;
    fprintf (pFile,"%s\n","            <DataArray type=\"Int32\" Name=\"types\" format=\"ascii\">");
    for (int elem = 0; elem < mesh.elemCard(); elem++)
    {
      ElemTypesMesh itemType = static_cast<ElemTypesMesh>(mesh.elemMesh.elemTypeFromDimType(mesh.elemMesh.typeFromLID(elem)));
      if (itemType != itemTypeOld)
      {
        numberOfSubElements = VTKHelperSpaMesh::getNumberOfSubElements(itemType,p_order);
        vtkCellType = VTKHelperSpaMesh::vtkItemType(itemType);
        itemTypeOld = itemType;
      }
      for (unsigned int i = 0; i < numberOfSubElements; i++)	  
      {
        fprintf (pFile, "%d ", vtkCellType);
      }
    }
    fprintf (pFile,"%s\n","");
    fprintf (pFile,"%s\n","            </DataArray>");
    fprintf (pFile,"%s\n","         </Cells>");

    fprintf (pFile,"%s\n","      </Piece>  ");        
    fprintf (pFile,"%s\n","   </UnstructuredGrid>");
    fprintf (pFile,"%s\n","</VTKFile>");
    fclose (pFile);       
  }

private:  
  template<typename TopologyType, unsigned int dim> void initSizes(const SpaMesh::Mesh<dim>& mesh, TopologyType& feDataAccessor, const DofMap& dofMap)
  {
    const unsigned int p_order = dofMap.getPOrder();
    _totNumberOfElementsToPlot = 0;
    _totNumberOfNodesToPlot = 0;
    _connectivitySize = 0;
  
    ElemTypesMesh itemTypeOld = SpaMesh::node;
    unsigned int numberOfSubElements;
    unsigned int numberOfSubNodes;
    unsigned int numberOfSubElementNodes;
    for (int elem = 0; elem < mesh.elemCard(); elem++)
    {
      ElemTypesMesh itemType = static_cast<ElemTypesMesh>(mesh.elemMesh.elemTypeFromDimType(mesh.elemMesh.typeFromLID(elem)));
      if (itemType != itemTypeOld)
      {
        numberOfSubElements = VTKHelperSpaMesh::getNumberOfSubElements(itemType,p_order);
        numberOfSubNodes = VTKHelperSpaMesh::getSubElementsTotalNodesNumber(itemType,p_order,numberOfSubElementNodes);
        itemTypeOld = itemType;
      }
      _totNumberOfElementsToPlot += numberOfSubElements; 
      _totNumberOfNodesToPlot += numberOfSubNodes; 
      _connectivitySize += numberOfSubElementNodes*numberOfSubElements;
    }
    _initialized = true;
  }
  
  template<typename TopologyType, unsigned int dim> void buildConnectivity(const SpaMesh::Mesh<dim>& mesh, TopologyType& feDataAccessor, const DofMap& dofMap)
  {
    assert(_initialized == true);
    _elementConnOffset.resize(_totNumberOfElementsToPlot);
    _elementTypes.resize(_totNumberOfElementsToPlot);
    _elementConn.resize(_connectivitySize);
    const unsigned int p_order = dofMap.getPOrder();
  
    Eigen::Matrix<unsigned int, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> numbering;
    unsigned int subElementsCounter = 0;
    unsigned int nodesCounter = 0;
    unsigned int connCounter = 0;
    unsigned int offset = 0;
    unsigned int vtkCellType = 0;
    unsigned int numberOfSubNodes;
    unsigned int numberOfSubElementNodes;
    ElemTypesMesh itemTypeOld = SpaMesh::node;
    for (int elem = 0; elem < mesh.elemCard(); elem++)
    {
      ElemTypesMesh itemType = mesh.elemMesh.elemTypeFromDimType(mesh.elemMesh.typeFromLID(elem));
      if (itemType != itemTypeOld)
      {
        VTKHelperSpaMesh::getSubElementsNodesNumbering(itemType,p_order,numbering);
        numberOfSubNodes = VTKHelperSpaMesh::getSubElementsTotalNodesNumber(itemType,p_order,numberOfSubElementNodes);
        vtkCellType = VTKHelperSpaMesh::vtkItemType(itemType);
        itemTypeOld = itemType;
      }
      for (unsigned int i = 0; i < numbering.rows(); i++)
        for (unsigned int j = 0; j < numbering.cols(); j++)
          _elementConn(connCounter++) = numbering(i,j) + nodesCounter; 
      nodesCounter += numberOfSubNodes;
      for (unsigned int i = 0; i < numbering.rows(); i++)	  
      {
        offset += numberOfSubElementNodes;
        _elementConnOffset[subElementsCounter] = offset;
        _elementTypes[subElementsCounter] = vtkCellType;
        subElementsCounter++;
      }
    }
  }
  
  template<typename TopologyType, unsigned int dim> void buildShapeMatrix(const SpaMesh::Mesh<dim>& mesh, TopologyType& feDataAccessor, const DofMap& dofMap)
  {
    assert(_initialized == true);
    const unsigned int p_order = dofMap.getPOrder();
    Quadrature::PointVector referenceSpaceNodes;
    ElemType elemTypeOld = NODEELEM;
    
    const unsigned int dofsNumber = FiniteElementSpaceType::FiniteElementBasisType::dofsNumber(Dim, p_order);
    _collBasisFunction.resizeAllCollective(dofsNumber, _totNumberOfNodesToPlot);
    _mappingFunction.resizeQpointsXCollective(_totNumberOfNodesToPlot);
    BasisFunctionType auxBasisFunction;
    MappingFunctionType auxMappingFunction;
    unsigned int qpIndex = 0;
    _collQPIndices.clear();
    _collQPIndices.reserve(mesh.elemCard()+1);
    _collQPIndices.push_back(qpIndex);
    for (int elem = 0; elem < mesh.elemCard(); elem++)
    {
      feDataAccessor.initData(mesh,elem);
      const ElemType elemType = feDataAccessor.getElemType();
      if (elemType != elemTypeOld)
      {
        VTKHelperSpaMesh::getSubElementsNodesE(elemType, p_order, referenceSpaceNodes);
        elemTypeOld = elemType;
      }
      FiniteElementSpaceType::initFENodalPoints(auxBasisFunction,p_order,auxMappingFunction,ReferenceFiniteElement::getGeoOrder(elemType),feDataAccessor,referenceSpaceNodes);
     
      unsigned int nodesSize = referenceSpaceNodes.rows();
      assert(nodesSize==auxBasisFunction.getShapeFunction().cols());
      assert(nodesSize==auxMappingFunction.getQpointX().rows());

      this->_collBasisFunction.getShapeFunction().block(0,qpIndex,dofsNumber,nodesSize) = auxBasisFunction.getShapeFunction();
      this->_mappingFunction.getQpointX().block(qpIndex,0,nodesSize,3) = auxMappingFunction.getQpointX();
      qpIndex += nodesSize;
      _collQPIndices.push_back(qpIndex);
    }
  }
 
  template<unsigned int dim>
  void write(const SpaMesh::Mesh<dim>& mesh, const DofMap& dofMap)
  {
    const unsigned int n_vars = dofMap.getNVars();
    const typename DofMap::VarNamesType var_names = dofMap.getVarNames();
    //write pvtu header
    if (FELib::processor_id() == 0)
    {
      char filename[1024];
      sprintf(filename,"output.pvtu");
      FILE* pFile=fopen(filename,"w");
      fprintf (pFile,"%s","<VTKFile type=\"PUnstructuredGrid\"");
      fprintf (pFile,"%s"," version=\"0.1\"");
      fprintf (pFile,"%s\n"," byte_order=\"LittleEndian\">");
      fprintf (pFile,"%s\n","   <PUnstructuredGrid GhostLevel=\"0\">");
      fprintf (pFile,"%s\n","         <PPoints>");
      fprintf (pFile,"%s%s%s\n","            <PDataArray type=\"Float32\" NumberOfComponents=\"","3","\">");
      fprintf (pFile,"%s\n","            </PDataArray>");
      fprintf (pFile,"%s\n","         </PPoints>");
      fprintf (pFile,"%s\n","         <PPointData>");
      for(unsigned int v = 0; v < n_vars; v++)
      {
        fprintf (pFile,"%s%s%s\n","            <PDataArray type=\"Float32\" Name=\"",var_names(v).c_str(),"\" NumberOfComponents=\"1\">");
        fprintf (pFile,"%s\n","            </PDataArray>");
      }
      fprintf (pFile,"%s\n","         </PPointData>");
      for (unsigned int i = 0; i < FELib::n_processors(); i++)
      {
        fprintf (pFile,"%s%02d%s\n","      <Piece  Source=\"output_",i,".vtu\">");
        fprintf (pFile,"%s\n","      </Piece>  ");        
      }
      fprintf (pFile,"%s\n","   </PUnstructuredGrid>");
      fprintf (pFile,"%s\n","</VTKFile>");
    }
  
    // write header
    char filename[1024];
    sprintf(filename,"output_%02d.vtu",FELib::processor_id());
    FILE* pFile=fopen(filename,"w");
    fprintf (pFile,"%s","<VTKFile type=\"UnstructuredGrid\"");
    fprintf (pFile,"%s"," version=\"0.1\"");
    fprintf (pFile,"%s\n"," byte_order=\"LittleEndian\">");
    fprintf (pFile,"%s\n","   <UnstructuredGrid>");
    fprintf (pFile,"%s","      <Piece  ");
    fprintf (pFile,"%s%d%s","NumberOfPoints=\"",_totNumberOfNodesToPlot,"\"  ");
    fprintf (pFile,"%s%d%s\n","NumberOfCells=\"",_totNumberOfElementsToPlot,"\">");
    
    // write mesh nodes  
    fprintf (pFile,"%s\n","         <Points>");
    fprintf (pFile,"%s%s%s\n","            <DataArray type=\"Float32\" NumberOfComponents=\"","3","\" format=\"ascii\">");
    Quadrature::PointVector& nodes = _mappingFunction.getQpointX();
    assert(nodes.rows() == _totNumberOfNodesToPlot);
    for (unsigned int n = 0; n < _totNumberOfNodesToPlot; n++)
      for(unsigned int i = 0; i < 3; i++)
         fprintf (pFile,"%f ",nodes(n,i));
    fprintf (pFile,"%s\n","");
    fprintf (pFile,"%s\n","            </DataArray>");
    fprintf (pFile,"%s\n","         </Points>");
  
    //  write solutions 
    fprintf (pFile,"%s\n","         <PointData>");
    for(unsigned int j=0;j<n_vars;++j)
    {
      fprintf (pFile,"%s%s%s\n","            <DataArray type=\"Float32\" Name=\"",var_names(j).c_str(),"\" NumberOfComponents=\"1\" format=\"ascii\">");
      assert(_nodal_solution.rows() == _totNumberOfNodesToPlot && _nodal_solution.cols() == n_vars);
      for(unsigned int i = 0; i < _totNumberOfNodesToPlot; i++)
        fprintf (pFile,"%f ", _nodal_solution(i,j));
      fprintf (pFile,"%s\n","");
      fprintf (pFile,"%s\n","            </DataArray>");
    }
    fprintf (pFile,"%s\n","         </PointData>");
  
    // write elements connectivity
    fprintf (pFile,"%s\n","         <Cells>");
    fprintf (pFile,"%s\n","            <DataArray type=\"Int32\" Name=\"connectivity\" format=\"ascii\">");
    for(unsigned int i = 0; i < _connectivitySize; i++)
      fprintf (pFile, "%d ", _elementConn(i));
    fprintf (pFile,"%s\n","");
    fprintf (pFile,"%s\n","            </DataArray>");
    
    // write elements offset 
    fprintf (pFile,"%s\n","            <DataArray type=\"Int32\" Name=\"offsets\" format=\"ascii\">");
    assert(_elementConnOffset.size() == _totNumberOfElementsToPlot);
    for(unsigned int i = 0; i < _totNumberOfElementsToPlot; i++)
      fprintf (pFile, "%d ", _elementConnOffset(i));
    fprintf (pFile,"%s\n","");
    fprintf (pFile,"%s\n","            </DataArray>");
  
    // write elements types 
    fprintf (pFile,"%s\n","            <DataArray type=\"Int32\" Name=\"types\" format=\"ascii\">");
    for(unsigned int i = 0; i < _totNumberOfElementsToPlot; i++)
      fprintf (pFile, "%d ", _elementTypes(i));
    fprintf (pFile,"%s\n","");
    fprintf (pFile,"%s\n","            </DataArray>");
    fprintf (pFile,"%s\n","         </Cells>");
  
    fprintf (pFile,"%s\n","      </Piece>  ");        
    fprintf (pFile,"%s\n","   </UnstructuredGrid>");
    fprintf (pFile,"%s\n","</VTKFile>");
    fclose (pFile);       
  }

  template<typename TopologyType, typename Derived, unsigned int dim> 
       void build_discontinuous_solution_vector(const SpaMesh::Mesh<dim>& mesh, TopologyType& feDataAccessor, const DofMap& dofMap, const Eigen::MatrixBase<Derived>& local_solution)
  {
    assert(_initialized == true);
    const unsigned int n_vars = dofMap.getNVars();
    const unsigned int p_order = dofMap.getPOrder();
    const typename DofMap::VarScalingsType var_scalings = dofMap.getVarScalings();
    _nodal_solution.resize(_totNumberOfNodesToPlot,n_vars);
    _mappingFunction.resizeQpointsXCollective(_totNumberOfNodesToPlot);
  
    Quadrature::PointVector referenceSpaceNodes;
    ElemType elemTypeOld = NODEELEM;
    BasisFunctionType basisFunction;
    MappingFunctionType auxMappingFunction;
    typename BasisFunctionType::ShapeMatrix& shapeMatrix = basisFunction.getShapeFunction();
    basisFunction.doComputeShapeDeriv(false);
 
    unsigned int nodesCounter = 0;
    for (int elem = 0; elem < mesh.elemCard(); elem++)
    {
      feDataAccessor.initData(mesh,elem);
      const ElemType elemType = feDataAccessor.getElemType();
      if (elemType != elemTypeOld)
      {
        VTKHelperSpaMesh::getSubElementsNodesE(elemType, p_order, referenceSpaceNodes);
        elemTypeOld = elemType;
      }

      FiniteElementSpaceType::initFENodalPoints(basisFunction,p_order,auxMappingFunction,ReferenceFiniteElement::getGeoOrder(elemType),feDataAccessor,referenceSpaceNodes);
      _mappingFunction.getQpointX().block(nodesCounter,0,referenceSpaceNodes.rows(),3) = auxMappingFunction.getQpointX();
      for (unsigned int var = 0; var < n_vars; var++)
        for (unsigned int j = 0; j < referenceSpaceNodes.rows(); j++)
	{
          _nodal_solution(j+nodesCounter,var) = (local_solution.row(elem)).segment(var*shapeMatrix.rows(),shapeMatrix.rows()) * shapeMatrix.col(j);
          _nodal_solution(j+nodesCounter,var) *= var_scalings[var];
	}
      nodesCounter+=referenceSpaceNodes.rows();
    }
    assert(nodesCounter == _totNumberOfNodesToPlot);
  }
  
  template<typename TopologyType, typename Derived, unsigned int dim> 
         void update_solution_vector(const SpaMesh::Mesh<dim>& mesh, TopologyType& feDataAccessor, const DofMap& dofMap, const Eigen::MatrixBase<Derived>& local_solution)
  {
    assert(_initialized == true);
    const unsigned int n_vars = dofMap.getNVars();
    const unsigned int p_order = dofMap.getPOrder();
    const typename DofMap::VarScalingsType var_scalings = dofMap.getVarScalings();
    _nodal_solution.resize(_totNumberOfNodesToPlot,n_vars);
  
    typename BasisFunctionType::ShapeMatrix& shapeMatrix = _collBasisFunction.getShapeFunction();
    
    unsigned int firstIndex, lastIndex;
    for (int elem = 0; elem < mesh.elemCard(); elem++)
    {
      firstIndex = _collQPIndices[elem];
      lastIndex = _collQPIndices[elem+1];
      for (unsigned int var = 0; var < n_vars; var++)
        for (unsigned int j = firstIndex; j < lastIndex; j++)
	{
          _nodal_solution(j,var) = (local_solution.row(elem)).segment(var*shapeMatrix.rows(),shapeMatrix.rows()) * shapeMatrix.col(j);
          _nodal_solution(j,var) *= var_scalings[var];
	}
    }
  }
  
  template<typename TopologyType, typename Derived, unsigned int dim> 
    void build_solution_vector(const SpaMesh::Mesh<dim>& mesh, TopologyType& feDataAccessor, const DofMap& dofMap, const Eigen::MatrixBase<Derived>& local_solution)
  {
  }
  
private:   
  CollectiveBasisFunction<Dim> _collBasisFunction;

  std::vector<unsigned int> _collQPIndices;
    
  CollectiveMappingFunction<MappingBasisFunction<Dim,Dim> > _mappingFunction;

  unsigned int _totNumberOfNodesToPlot;
  
  unsigned int _totNumberOfElementsToPlot;

  unsigned int _connectivitySize;

  Eigen::Matrix<double,Eigen::Dynamic,Eigen::Dynamic> _nodal_solution; 

  Eigen::Matrix<unsigned int,Eigen::Dynamic,1> _elementConnOffset; 
  
  Eigen::Matrix<unsigned int,Eigen::Dynamic,1> _elementConn; 
  
  Eigen::Matrix<unsigned int,Eigen::Dynamic,1> _elementTypes; 

  bool _initialized;
};



#endif
