/*
 * vtkForCenos.cpp
 *
 *  Created on: Oct 11, 2017
 *      Author: vadims
 */

#include <boost/filesystem.hpp>
#include <chrono>
#include <complex>

#include <vtkVersion.h>
#include <vtkSmartPointer.h>
#include <vtkTetra.h>
#include <vtkCellArray.h>
#include <vtkXMLUnstructuredGridReader.h>
#include <vtkDataSetMapper.h>
#include <vtkActor.h>
#include <vtkRenderer.h>
#include <vtkRenderWindow.h>
#include <vtkRenderWindowInteractor.h>
#include <vtkPolyDataWriter.h>
#include <vtkXMLUnstructuredGridWriter.h>
#include <vtkUnstructuredGrid.h>
#include <vtkPointData.h>
#include <vtkVertexGlyphFilter.h>
#include <vtkHexahedron.h>
#include <vtkWedge.h>
#include <vtkPyramid.h>
#include <vtkQuad.h>
#include <vtkTriangle.h>
#include <vtkLine.h>
#include <vtkDoubleArray.h>
#include <vtkUnsignedIntArray.h>
#include <vtkCellData.h>
#include <vtkMultiBlockDataSet.h>
#include <vtkXMLMultiBlockDataWriter.h>
#include <vtkXMLCompositeDataWriter.h>
#include <vtkCollection.h>
#include <vtkDataSetCollection.h>
#include <vtkXMLDataSetWriter.h>
#include <vtkCompositeDataSet.h>
#include <vtkDataWriter.h>
#include <vtkEnSightWriter.h>
#include <vtkPointDataToCellData.h>
#include <vtkStringArray.h>
#include <vtkAppendFilter.h>
#include <vtkGeometryFilter.h>
#include <vtkPolyDataWriter.h>
#include <vtkInformation.h>
#include <vtkStreamingDemandDrivenPipeline.h>
#include <vtkXMLPolyDataWriter.h>
#include <vtkAppendPolyData.h>
#include <vtkGenericDataObjectReader.h>
#include <vtkDataArrayCollection.h>
#include <vtkDataArray.h>
#include <vtkCleanPolyData.h>
#include <vtkProbeFilter.h>
#include <vtkDelaunay3D.h>


#include <vtkUnstructuredGridAlgorithm.h>
#include <vtkExtractUnstructuredGrid.h>
#include <vtkBooleanOperationPolyDataFilter.h>



#include <limits>

//#include <vtkXMLDataObjectWriter.h>



#include "vtkForCenos.hpp"

void vtkForCenos::writeVTKlog(const std::string &text)
{
    
	char const* tmp = getenv( "TEMP" );
    std::string tempFolder( tmp );
	time_t t = time(0);
	struct tm * now = localtime( & t );
    std::ofstream log_file(
    	tempFolder + "\\cenos\\" + std::to_string(now->tm_year+1900) + "-" + std::to_string(now->tm_mon+1) +
		"-" + std::to_string(now->tm_mday) + ".cenos-vtk.log",
        std::ios_base::out | std::ios_base::app );
    std::string timeStr = std::to_string(now->tm_hour) + ":" + std::to_string(now->tm_min) +
			":" + std::to_string(now->tm_sec) + ". ";
    log_file << timeStr;
    log_file << text;
    log_file << std::endl;
    log_file.close();
    
}

vtkForCenos::vtkForCenos()
{
	this->multiBlockDataSet =
	     vtkSmartPointer<vtkMultiBlockDataSet>::New();
	this->multiPieceDataSet =
	     vtkSmartPointer<vtkMultiPieceDataSet>::New();
	this->_ens_writer =
	     vtkSmartPointer<cenosEnsightWriter>::New();
	this->_ens_reader =
	     vtkSmartPointer<vtkEnSightGoldBinaryReader >::New();
	this->unstructuredGrid =
	     vtkSmartPointer<vtkUnstructuredGrid>::New();
	nrOfSteps = 0;
}

int vtkForCenos::newTimeStep()
{
	for(auto &block: ugridVec)
	{
//		block->GetPointData()->Delete();
	}
	 return EXIT_SUCCESS;
}

int vtkForCenos::addBlockToGrid(std::vector<std::vector<float> > nodeArray,
		std::vector<std::vector<int> > elementArray, int entNr, float meshScale)
{
  vtkSmartPointer<vtkPoints> points =
    vtkSmartPointer<vtkPoints>::New();
  vtkSmartPointer<vtkUnstructuredGrid> ugrid =
	vtkSmartPointer<vtkUnstructuredGrid>::New();

  std::map<int, int> nodesPos;
  for (int j=0; j< nodeArray.size(); j++)
  {
	  nodesPos[nodeArray[j][0]] = j;
  }

   for (auto &p1: nodeArray)
  {
	  points->InsertNextPoint(p1[1]*meshScale, p1[2]*meshScale, p1[3]*meshScale);		//SetPoint will be faster!!!
  }

	ugrid->SetPoints(points);

  vtkSmartPointer<vtkIntArray> bID =
    vtkSmartPointer<vtkIntArray>::New();
  bID->SetNumberOfComponents(1);
  bID->SetName("BlockId");


  for (auto &e1: elementArray)
  {
	  if (e1[0] == 12)
	  {
		  //HEXA
		    vtkSmartPointer<vtkHexahedron> hexa =
		      vtkSmartPointer<vtkHexahedron>::New();
			for (int i=0; i<8; i++)
			{
				hexa->GetPointIds()->SetId(i, nodesPos[e1[i+1]]);
			}
			ugrid->InsertNextCell(hexa->GetCellType(),hexa->GetPointIds());
			bID->InsertNextValue(entNr);
	  }
	  if (e1[0] == 13)
	  {
		  //PRISM
		    vtkSmartPointer<vtkWedge> prism =
		      vtkSmartPointer<vtkWedge>::New();
			for (int i=0; i<6; i++)
			{
				prism->GetPointIds()->SetId(i, nodesPos[e1[i+1]]);
			}
			ugrid->InsertNextCell(prism->GetCellType(),prism->GetPointIds());
			bID->InsertNextValue(entNr);


	  }
	  if (e1[0] == 14)
	  {
		  //PYRAMID
		    vtkSmartPointer<vtkPyramid> pyra =
		      vtkSmartPointer<vtkPyramid>::New();
			for (int i=0; i<5; i++)
			{
				pyra->GetPointIds()->SetId(i,nodesPos[e1[i+1]]);
			}
			ugrid->InsertNextCell(pyra->GetCellType(),pyra->GetPointIds());
			bID->InsertNextValue(entNr);


	  }
	  if (e1[0] == 10)
	  {
		  //TETRA
		    vtkSmartPointer<vtkTetra> tetra =
		      vtkSmartPointer<vtkTetra>::New();
			for (int i=0; i<4; i++)
			{
				  tetra->GetPointIds()->SetId(i, nodesPos[e1[i+1]]);
			}
			ugrid->InsertNextCell(tetra->GetCellType(),tetra->GetPointIds());
			bID->InsertNextValue(entNr);

	  }
	  if (e1[0] == 9)
	  {
		  //QUAD
		    vtkSmartPointer<vtkQuad> quad =
		      vtkSmartPointer<vtkQuad>::New();
			for (int i=0; i<4; i++)
			{
				quad->GetPointIds()->SetId(i, nodesPos[e1[i+1]]);

			}
			ugrid->InsertNextCell(quad->GetCellType(),quad->GetPointIds());
			bID->InsertNextValue(entNr);

	  }
	  if (e1[0] == 5)
	  {
		  //TRIANGLE
		    vtkSmartPointer<vtkTriangle> triangle =
		      vtkSmartPointer<vtkTriangle>::New();
			for (int i=0; i<3; i++)
			{
				triangle->GetPointIds()->SetId(i, nodesPos[e1[i+1]]);

			}
			ugrid->InsertNextCell(triangle->GetCellType(),triangle->GetPointIds());
			bID->InsertNextValue(entNr);

	  }
	  if (e1[0] == 3)
	  {
		  //LINE
		    vtkSmartPointer<vtkLine> line =
		      vtkSmartPointer<vtkLine>::New();
			for (int i=0; i<2; i++)
			{
				line->GetPointIds()->SetId(i, nodesPos[e1[i+1]]);
			}
			ugrid->InsertNextCell(line->GetCellType(),line->GetPointIds());
			bID->InsertNextValue(entNr);
	  }
  }

  ugrid->GetCellData()->AddArray(bID);


  ugridVec.push_back(ugrid);

  indMap[entNr] = ugridVec.size()-1;
  int currentBlocks = multiBlockDataSet->GetNumberOfBlocks();
  multiBlockDataSet->SetNumberOfBlocks(currentBlocks+1);
//  multiBlockDataSet->SetBlock(currentBlocks, ugrid);

  return EXIT_SUCCESS;
}


int vtkForCenos::addUGrid(vtkSmartPointer<vtkUnstructuredGrid> ugrid, int entityId)
{
    ugridVec.push_back(ugrid);

    indMap[entityId] = ugridVec.size() - 1;
    int currentBlocks = multiBlockDataSet->GetNumberOfBlocks();
    multiBlockDataSet->SetNumberOfBlocks(currentBlocks + 1);
    return EXIT_SUCCESS;
}

vtkSmartPointer<vtkUnstructuredGrid> vtkForCenos::getGrid(int blockId)
{
    return ugridVec[indMap[blockId]];
}


int vtkForCenos::setBlockNames(std::vector<int> ids, std::vector<std::string> blockNames)
{
	_ens_writer->setNameList(ids, blockNames );
	ugridBlockNames = blockNames;
	  return EXIT_SUCCESS;
}

int vtkForCenos::writeTimeStep(double time)
{
	writeVTKlog("writing time step " + std::to_string(time));
    vtkSmartPointer<vtkAppendFilter> appender = vtkSmartPointer<vtkAppendFilter>::New();
    writeVTKlog("writing time step 1 " + std::to_string(ugridVec.size()));

    for (int i = 0; i < ugridVec.size(); ++i)
    {
        fillEmptySets(i);
    }

    for (int i = 0; i < ugridVec.size(); ++i)
    {
        multiBlockDataSet->SetBlock(i, ugridVec[i]);
    }

    for (int i = 0; i < multiBlockDataSet->GetNumberOfBlocks(); ++i)
    {
        appender->AddInputData(multiBlockDataSet->GetBlock(i));
    }

	appender->Update();

	_ens_writer->SetInputData( appender->GetOutput() );

	_ens_writer->SetTimeStep(nrOfSteps);

	tsteps.push_back(time);

	_ens_writer->Write();

	writeVTKlog("finished writing time step!");

	nrOfSteps++;

	return EXIT_SUCCESS;
}

void vtkForCenos::addTimeStepValue(double time)
{
	writeVTKlog("set time step " + std::to_string(time));
	_ens_writer->SetTimeStep(nrOfSteps);
	tsteps.push_back(time);
	nrOfSteps++;

}


int vtkForCenos::setOutput(std::string path, std::string basename)
{
	  boost::filesystem::path resdir(path);
		  if (!boost::filesystem::exists(resdir))
		  {
			  if (!boost::filesystem::create_directory(resdir))
			  {
				  std::cout << "Could not create directory " << path << "/results" << std::endl;
				  // ******************
				  return 0; 	// Needs proper return!!!!!
				  // ***************

			  }
		  }
	  _ens_writer->SetPath(path.c_str());
	  _ens_writer->SetBaseName(basename.c_str());
	  return EXIT_SUCCESS;
}

int vtkForCenos::finalize()
{
	writeVTKlog("finalizing case writing.");
	writeVTKlog("nr of timesteps. " + std::to_string(nrOfSteps));

	_ens_writer->WriteOCaseFile(nrOfSteps, tsteps);
	writeVTKlog("Finished case writing.");

	return EXIT_SUCCESS;
}

void vtkForCenos::enterDataSet(std::string dataset, int dim)
{
	bool found = false;
	for (auto &ds: dataSets)
	{
		if (ds.first == dataset)
		{
			found = true;
		}
	}
	if (!found)
	{
		dataSets.push_back(std::make_pair(dataset,dim));
	}
}

std::string vtkForCenos::mergeVTK(std::vector<std::string> fileList, std::vector<int> idList, std::string path, std::string domainName, std::vector<float>& boundVector)
{
	writeVTKlog("mergeVTK " );

	 vtkSmartPointer<vtkAppendPolyData> appendFilter =
	    vtkSmartPointer<vtkAppendPolyData>::New();
	 std::vector<vtkSmartPointer<vtkPolyData>> inputs;
	for (int fi=0; fi< fileList.size(); fi++)
	{
		std::string inpFileName = fileList[fi];
		std::replace(inpFileName.begin(), inpFileName.end(), '\\', '/');
		std::size_t pos = inpFileName.find(".vtk");
		inpFileName = inpFileName.substr(0,pos+4);

		 vtkSmartPointer<vtkPolyData> input =
		    vtkSmartPointer<vtkPolyData>::New();
	    vtkSmartPointer<vtkGenericDataObjectReader> reader =
	      vtkSmartPointer<vtkGenericDataObjectReader>::New();
	    reader->SetFileName(inpFileName.c_str());
	    reader->Update();

		 vtkSmartPointer<vtkCellArray> emptyCellArray =
		    vtkSmartPointer<vtkCellArray>::New();

	    vtkSmartPointer<vtkPolyData> polyData = reader->GetPolyDataOutput();
	    polyData->SetVerts(emptyCellArray);
	    polyData->SetLines(emptyCellArray);

		vtkSmartPointer<vtkUnsignedIntArray> dataN =
		    vtkSmartPointer<vtkUnsignedIntArray>::New();
		  dataN->SetNumberOfComponents(1);
		  dataN->SetName("boundaryId");

		 for(vtkIdType i = 0; i < polyData->GetNumberOfCells(); i++)
			{
				 dataN->InsertNextValue(idList[fi]);
			}

		 //polyData->GetCellData()->AddArray(dataN);
		 polyData->GetCellData()->SetScalars(dataN);

		input->ShallowCopy(polyData);

	    inputs.push_back(input);
		#if VTK_MAJOR_VERSION <= 5
		  appendFilter->AddInputConnection(inputs.back()->GetProducerPort());
		#else
		  appendFilter->AddInputData(inputs.back());
		#endif
	}

	  appendFilter->Update();

	  vtkSmartPointer<vtkCleanPolyData> cleanFilter =
	    vtkSmartPointer<vtkCleanPolyData>::New();
	  cleanFilter->SetInputConnection(appendFilter->GetOutputPort());
	  cleanFilter->Update();

		vtkSmartPointer<vtkPolyDataWriter> polywriter =
				vtkSmartPointer<vtkPolyDataWriter>::New();
		std::string filename = path + "\\vtk\\" + domainName +".vtk";
		writeVTKlog("writing VTK to file " + filename );

		polywriter->SetFileName(filename.c_str());
		#if VTK_MAJOR_VERSION <= 5
		  polywriter->SetInput(appendFilter->GetOutput());
		#else
		  polywriter->SetInputData(appendFilter->GetOutput());
		#endif

		for (int i=0; i<6; i++)
		{
			boundVector.push_back(polywriter->GetInput()->GetBounds()[i]);
		}

		polywriter->Write();

		return "vtk\\" + domainName +".vtk";
}


void vtkForCenos::fillEmptySets(int entNr)
{
	unsigned int numberOfArrays = ugridVec[entNr]->GetPointData()->GetNumberOfArrays();

	for (auto &ds: dataSets)
	{

		bool found = false;
		for(unsigned int i = 0; i < numberOfArrays; i++)
	    {
            auto t = ugridVec[entNr]->GetPointData()->GetArrayName(i);
			if (ugridVec[entNr]->GetPointData()->GetArrayName(i) == ds.first)
			{
				found = true;
			}
	    }
		if (!found)
		{
			if (ds.second==1)
				addZeroPointFieldToBlock_(ds.first, entNr);
			else if (ds.second==3)
				addZeroPointVectorFieldToBlock_(ds.first, entNr);
		}
	}
}

int vtkForCenos::addCurrentPointFieldToBlock(std::ifstream & resFile, std::string fieldName, int blockId)
{
	enterDataSet(fieldName, 1);
	vtkSmartPointer<vtkDoubleArray> dataN =
	    vtkSmartPointer<vtkDoubleArray>::New();
	  dataN->SetNumberOfComponents(1);
	  dataN->SetName(fieldName.c_str());

	std::string line;
	getline(resFile,line);
    int count = 0;
    std::istringstream cistr(line);
    cistr >> count;

	writeVTKlog("writing data " + fieldName + " in block " + std::to_string(blockId) + " in full");


	if (!resFile.is_open())
	{
		writeVTKlog("could not open file for read " + std::to_string(blockId) + ", assuming zeros for  " + fieldName );
        addZeroPointFieldToBlock(fieldName, blockId);
		return EXIT_SUCCESS;
	}
    std::vector<int> testvec;
    while (resFile.good())
    {
		getline(resFile,line);
		if (line.empty())
			break;
		std::istringstream lString(line);
		int nr;
		float v1;
		lString >> nr >> v1;
		dataN->InsertNextValue(v1);
        testvec.push_back(v1);
    }

    if (dataN->GetDataSize() != count)
    {
        std::cerr << "Field " << fieldName << " in entity nr. " << blockId << " does not contain all data points (" << dataN->GetDataSize() << "/" << count << ". This field will be ignored" << std::endl;
        return EXIT_SUCCESS;
    }

	ugridVec[indMap[blockId]]->GetPointData()->AddArray(dataN);
	writtenFields_Blocks.push_back(fieldName + "_" + std::to_string(blockId));
	return EXIT_SUCCESS;
}

int vtkForCenos::addZeroPointFieldToBlock(std::string fieldName, int blockId)
{
	vtkSmartPointer<vtkDoubleArray> dataN =
	    vtkSmartPointer<vtkDoubleArray>::New();
	  dataN->SetNumberOfComponents(1);
	  dataN->SetName(fieldName.c_str());

		writeVTKlog("writing data " + fieldName + " in block " + std::to_string(blockId) + " in 0");

	 for(vtkIdType i = 0; i < ugridVec[indMap[blockId]]->GetNumberOfPoints(); i++)
		{
		 	 dataN->InsertNextValue(0);
		}

	  ugridVec[indMap[blockId]]->GetPointData()->AddArray(dataN);
	  return EXIT_SUCCESS;
}

int vtkForCenos::addZeroPointFieldToBlock_(std::string fieldName, int blockId)
{
	vtkSmartPointer<vtkDoubleArray> dataN =
	    vtkSmartPointer<vtkDoubleArray>::New();
	  dataN->SetNumberOfComponents(1);
	  dataN->SetName(fieldName.c_str());

      writeVTKlog("writing data " + fieldName + " in block " + std::to_string(blockId) + " in 0, autofill");


	 for(vtkIdType i = 0; i < ugridVec[blockId]->GetNumberOfPoints(); i++)
		{
		 	 dataN->InsertNextValue(0);
		}

	  ugridVec[blockId]->GetPointData()->AddArray(dataN);
	  return EXIT_SUCCESS;
}

int vtkForCenos::addZeroPointVectorFieldToBlock(std::string fieldName, int blockId)
{
	vtkSmartPointer<vtkDoubleArray> dataN =
	    vtkSmartPointer<vtkDoubleArray>::New();
	  dataN->SetNumberOfComponents(3);
	  dataN->SetName(fieldName.c_str());

      writeVTKlog("writing data " + fieldName + " in block " + std::to_string(blockId) + " in 0, autofill");


	 for(vtkIdType i = 0; i < ugridVec[indMap[blockId]]->GetNumberOfPoints(); i++)
		{
		    float tuple[3];
		    tuple[0] = 0;
		    tuple[1] = 0;
		    tuple[2] = 0;
		    dataN->InsertNextTuple(tuple);
		}

	  ugridVec[indMap[blockId]]->GetPointData()->AddArray(dataN);
	  return EXIT_SUCCESS;
}

int vtkForCenos::addZeroPointVectorFieldToBlock_(std::string fieldName, int blockId)
{
	vtkSmartPointer<vtkDoubleArray> dataN =
	    vtkSmartPointer<vtkDoubleArray>::New();
	  dataN->SetNumberOfComponents(3);
	  dataN->SetName(fieldName.c_str());

      writeVTKlog("writing data " + fieldName + " in block " + std::to_string(blockId) + " in 0, autofill");

	 for(vtkIdType i = 0; i < ugridVec[blockId]->GetNumberOfPoints(); i++)
		{
		    float tuple[3];
		    tuple[0] = 0;
		    tuple[1] = 0;
		    tuple[2] = 0;
		    dataN->InsertNextTuple(tuple);
		}

	 ugridVec[blockId]->GetPointData()->AddArray(dataN);
	  return EXIT_SUCCESS;
}

int vtkForCenos::addCurrentCellFieldToBlock(std::ifstream & resFile, std::string fieldName, int blockId)
{
	enterDataSet(fieldName,1);
	vtkSmartPointer<vtkDoubleArray> dataN =
	    vtkSmartPointer<vtkDoubleArray>::New();
	  dataN->SetNumberOfComponents(1);
	  dataN->SetName(fieldName.c_str());

	std::string line;
	std::vector<int> elSize = getElementSizeArray(blockId);
	writeVTKlog( "entering reading cell field " + std::to_string(blockId) + " " + fieldName );

	int elNr=0;
	if (!resFile.is_open())
	{
		writeVTKlog( "could not open file for read " + std::to_string(blockId) + " " + fieldName );
		return 0;
	}

    while (resFile.good())
    {
		getline(resFile,line);
		if (line.empty())
			break;
		int ncount = elSize[elNr];
		std::size_t pos = line.find("   ");


		std::istringstream lString(line.substr(pos));

	//	writeVTKlog( "ncount " + std::to_string(ncount) );
/*
		for (int j = 0; j<ncount*3; j++)
		{
			lString	>> vl;
			writeVTKlog( "coord " + std::to_string(vl) );
		}
		*/

		float val1 = 0;
		if (ncount==4)
		{
			float v1, v2, v3, v4;
			lString	>> v1 >> v2 >> v3 >> v4;
			val1 = (v1 + v2 + v3 + v4)/4;
		}
		else
		{
			for (int j = 0; j<ncount; j++)
			{
				float v1;
				lString	>> v1;
				val1 = (val1*j + v1)/(j+1);
			}
		}

		dataN->InsertNextValue(val1);


		elNr++;
    }

	writeVTKlog( "done reading element field for block " + std::to_string(blockId) + " " + fieldName );

	  ugridVec[indMap[blockId]]->GetCellData()->AddArray(dataN);
	  writtenFields_Blocks.push_back(fieldName + "_" + std::to_string(blockId));
	  return EXIT_SUCCESS;
}


bool vtkForCenos::fieldIsAdded( std::string fieldName, int blockId)
{
	if (std::find(writtenFields_Blocks.begin(), writtenFields_Blocks.end(),fieldName + "_" + std::to_string(blockId))!=writtenFields_Blocks.end())
		return true;
	else
		return false;
}

int vtkForCenos::addZeroCellFieldToBlock( std::string fieldName, int blockId)
{
	vtkSmartPointer<vtkDoubleArray> dataN =
	    vtkSmartPointer<vtkDoubleArray>::New();
	  dataN->SetNumberOfComponents(1);
	  dataN->SetName(fieldName.c_str());

	std::vector<int> elSize = getElementSizeArray(blockId);


	 for(vtkIdType i = 0; i < ugridVec[indMap[blockId]]->GetNumberOfCells(); i++)
		{
		 	 dataN->InsertNextValue(0);
		}

	  ugridVec[indMap[blockId]]->GetCellData()->AddArray(dataN);
	  return EXIT_SUCCESS;
}

int vtkForCenos::addCurrentPointVectorFieldToBlock(std::ifstream & resFile, std::string fieldName, int blockId)
{
	enterDataSet(fieldName,3);
	vtkSmartPointer<vtkDoubleArray> dataN =
	    vtkSmartPointer<vtkDoubleArray>::New();
	  dataN->SetNumberOfComponents(3);
	  dataN->SetName(fieldName.c_str());

	std::string line;
	getline(resFile,line);
	writeVTKlog("writing data " + fieldName + " in block " + std::to_string(blockId) + " in full");


	if (!resFile.is_open())
	{
		std::cout << "could not open file for read " << blockId << " " << fieldName <<  std::endl;
		return 0;
	}
    while (resFile.good())
    {
		getline(resFile,line);
		if (line.empty())
			break;
		std::istringstream lString(line);
		int nr;
		float v1, v2, v3;
		lString >> nr >> v1 >> v2 >> v3;
		//while (lString	>> v1);
	    float tuple[3];
	    tuple[0] = v1;
	    tuple[1] = v2;
	    tuple[2] = v3;
	    dataN->InsertNextTuple(tuple);
    }

	  ugridVec[indMap[blockId]]->GetPointData()->AddArray(dataN);
	  writtenFields_Blocks.push_back(fieldName + "_" + std::to_string(blockId));
	  return EXIT_SUCCESS;
}


void vtkForCenos::writePolyDataToVTK(int entNr, std::string outFileName)
{

    vtkSmartPointer<vtkGeometryFilter> geometryFilter =
        vtkSmartPointer<vtkGeometryFilter>::New();

    geometryFilter->SetInputData(ugridVec[indMap[entNr]]);
    geometryFilter->Update();

    vtkSmartPointer<vtkPolyData> polydata = geometryFilter->GetOutput();

    vtkSmartPointer<vtkPolyDataWriter> writer = vtkSmartPointer<vtkPolyDataWriter>::New();
    writer->SetInputData(polydata);
    writer->SetFileName(outFileName.c_str());
    writer->Update();
    writer->Write();
}

std::vector<float> vtkForCenos::getPolyPoints(int entNr)
{

	  vtkSmartPointer<vtkGeometryFilter> geometryFilter =
	    vtkSmartPointer<vtkGeometryFilter>::New();

	  geometryFilter->SetInputData(ugridVec[indMap[entNr]]);

	  geometryFilter->Update();
      vtkSmartPointer<vtkPolyData> polydata = geometryFilter->GetOutput();

	  std::vector<float> outVec;
	  for(vtkIdType i = 0; i < polydata->GetNumberOfPoints(); i++)
	  {
	    double p[3];
	    polydata->GetPoint(i,p);

	    outVec.push_back(p[0]);

	    outVec.push_back(p[1]);

	    outVec.push_back(p[2]);
	  }

	  return outVec;
}

std::vector<int> vtkForCenos::getElementSizeArray(int blockId)
{
	  std::vector<int> outVec;
	  vtkSmartPointer<vtkIdList> cellPointIds =
	    vtkSmartPointer<vtkIdList>::New();
      vtkSmartPointer<vtkCellArray> cellArray = ugridVec[indMap[blockId]]->GetCells();

	  for(vtkIdType i = 0; i < cellArray->GetNumberOfCells(); i++)
	  {
		  cellArray->GetNextCell(cellPointIds);

		  outVec.push_back(cellPointIds->GetNumberOfIds());
	  }

	  return outVec;

}

std::vector<int> vtkForCenos::getPolyElements(int entNr)
{
	  vtkSmartPointer<vtkGeometryFilter> geometryFilter =
	    vtkSmartPointer<vtkGeometryFilter>::New();
	  geometryFilter->SetInputData(ugridVec[indMap[entNr]]);
	  geometryFilter->Update();
      vtkSmartPointer<vtkPolyData> polydata = geometryFilter->GetOutput();

	  std::vector<int> outVec;

	  vtkSmartPointer<vtkIdList> cellPointIds =
	    vtkSmartPointer<vtkIdList>::New();

	  for(vtkIdType i = 0; i < polydata->GetNumberOfCells(); i++)
	  {
		  polydata->GetCellPoints(i,cellPointIds);

		  outVec.push_back(cellPointIds->GetNumberOfIds());
		  for(vtkIdType j = 0; j < cellPointIds->GetNumberOfIds(); j++)
		  {
			  outVec.push_back(cellPointIds->GetId(j));
		  }
	  }

	  return outVec;

}

void vtkForCenos::writeLegacyVTK(float time, std::vector<std::string> postValuesForExport)
{
	std::string path = _ens_writer->GetPath();

	for (int i=0; i< ugridVec.size(); i++)
	{
		vtkSmartPointer<vtkGeometryFilter> geometryFilter =
				vtkSmartPointer<vtkGeometryFilter>::New();
		geometryFilter->SetInputData(ugridVec[i]);
		geometryFilter->Update();
        vtkSmartPointer<vtkPolyData> polydata = geometryFilter->GetOutput();

        vtkSmartPointer<vtkPointData> pointData = polydata->GetPointData();
		std::vector<std::string> pointDataNamesForRemove;
		for (int i = 0; i < pointData->GetNumberOfArrays(); i++)
		{
			std::string dataName = pointData->GetArray(i)->GetName();
		    if (std::find(postValuesForExport.begin(), postValuesForExport.end(), dataName) == postValuesForExport.end())
		    {
		    	pointDataNamesForRemove.push_back(dataName);
		    }
		}


		for (auto rmData: pointDataNamesForRemove)
		{
	    	polydata->GetPointData()->RemoveArray( rmData.c_str() );
		}


        if (pointData->GetNumberOfArrays() > 0)
        {
            vtkSmartPointer<vtkPolyDataWriter> polywriter =
                vtkSmartPointer<vtkPolyDataWriter>::New();
            std::string filename = path + "/" + ugridBlockNames[i] + "_" + std::to_string(time) + ".vtk";
            polywriter->SetFileName(filename.c_str());
            polywriter->SetInputData(polydata);
            polywriter->Write();
        }
	}
}

void vtkForCenos::setTransientGeometry(bool tg)
{
	_ens_writer->SetTransientGeometry(tg);
}

float vtkForCenos::getLastTimeInSet()
{
	writeVTKlog(" getLastTimeInSet");

	std::string path = _ens_writer->GetPath() ;
	path.append(std::string("/") + _ens_writer->GetBaseName() + std::string(".0.case") );
	_ens_reader->SetCaseFileName(path.c_str());
	_ens_reader->Update();

	float lastTimeInSet = 0;
	if (_ens_reader->GetTimeSets()->GetNumberOfItems() >0 )
	{
		vtkSmartPointer<vtkDataArray> timeSet = _ens_reader->GetTimeSets()->GetItem(0);

		for (int i=0; i< timeSet->GetNumberOfTuples(); i++)
		{
			lastTimeInSet = timeSet->GetTuple1(i);
		}
	}
	else
	{
		lastTimeInSet = 0;
	}

	return 	lastTimeInSet;
}

void vtkForCenos::setContinuation(double zerothTimeValue)
{
	writeVTKlog(" Setting continuation");
	_ens_writer->SetTransientGeometry(true);

	std::string path = _ens_writer->GetPath() ;
	path.append(std::string("/") + _ens_writer->GetBaseName() + std::string(".0.case") );
	_ens_reader->SetCaseFileName(path.c_str());
	_ens_reader->Update();


//	_ens_reader->GetOutput()->GetDataSet(0)->GetFieldData()->GetArray("BlockId");
    
	std::ostringstream stream;

	vtkSmartPointer<vtkAppendFilter> appender = vtkSmartPointer<vtkAppendFilter>::New();
	for(int i=0; i<_ens_reader->GetOutput()->GetNumberOfBlocks(); ++i)
	{
		appender->AddInputData(_ens_reader->GetOutput()->GetBlock(i));
	}

	appender->Update();

	_ens_writer->SetInputData( appender->GetOutput() );

	if (_ens_reader->GetTimeSets()->GetNumberOfItems() >0 )
	{
		vtkSmartPointer<vtkDataArray> timeSet = _ens_reader->GetTimeSets()->GetItem(0);
		for (int i=0; i< timeSet->GetNumberOfTuples(); i++)
		{
			double tuple;
			tuple = timeSet->GetTuple1(i);
			_ens_writer->SetTimeStep(i);
			tsteps.push_back(tuple);
			nrOfSteps++;
		}
	}
	else
	{
		_ens_writer->SetTimeStep(0);
		tsteps.push_back(zerothTimeValue);
		nrOfSteps++;
	}

	vtkDataArray *Blocks=appender->GetOutput()->GetCellData()->GetScalars("BlockId");
	if (Blocks==NULL) //|| strcmp(BlockData->GetName(),"BlockId"))
	  {
		writeVTKlog("Warning: No BlockIDs were found in source.");
	  }
	else
	{
		writeVTKlog("Blocks were found");
	}


//	_ens_writer->Write();

//	_ens_writer->WriteOCaseFile(nrOfSteps, tsteps);
}





void cenosEnsightWriter::WriteOCaseFile(int TotalTimeSteps, std::vector<double> tsteps)
{
  vtkUnstructuredGrid *input=this->GetInput();
  int i;

  this->ComputeNames();

  if (!this->BaseName)
  {
    vtkErrorMacro("A FileName or Path/BaseName must be specified.");
    return;
  }


  char charBuffer[512];
  snprintf(charBuffer,sizeof(charBuffer),"%s/%s.%d.case",this->Path,this->BaseName,this->ProcessNumber);

  //open the geometry file
  FILE *fd=nullptr;
  if (!(fd=OpenFile(charBuffer)))
  {
    return;
  }

  this->WriteTerminatedStringToFile("FORMAT\n",fd);
  this->WriteTerminatedStringToFile("type: ensight gold\n\n",fd);
  this->WriteTerminatedStringToFile("\nGEOMETRY\n",fd);

  //write the geometry file
  if (!this->TransientGeometry)
  {
    snprintf(charBuffer,sizeof(charBuffer),"model: %s.%d.00000.geo\n",this->BaseName,this->ProcessNumber);
    this->WriteTerminatedStringToFile(charBuffer,fd);
  }
  else
  {
    snprintf(charBuffer,sizeof(charBuffer),"model: 1 %s.%d.*****.geo\n",this->BaseName,this->ProcessNumber);
    this->WriteTerminatedStringToFile(charBuffer,fd);
  }

  this->WriteTerminatedStringToFile("\nVARIABLE\n",fd);

  char fileBuffer[512];

  //write the Node variable files
  for (i=0;i<input->GetPointData()->GetNumberOfArrays();i++)
  {

    strcpy(fileBuffer,input->GetPointData()->GetArray(i)->GetName());
    // skip arrays that were not written
    if (strcmp(fileBuffer,"GlobalElementId")==0)
    {
      continue;
    }
    if (strcmp(fileBuffer,"GlobalNodeId")==0)
    {
      continue;
    }
    if (strcmp(fileBuffer,"BlockId")==0)
    {
      continue;
    }
    // fileBuffer will be sanitized, and forward slashes removed.
    // forward slash is division sign in physics, thus cannot be removed!
    // fieldName will be used for writing variable names, and sanitized fileName for file names
    char fieldName[512];
    strncpy(fieldName, fileBuffer, 512);

    this->SanitizeFileName(fileBuffer);

    //figure out what kind of data it is
    char SmallBuffer[16];
    switch(input->GetPointData()->GetArray(i)->GetNumberOfComponents())
    {
    case(1):
      strcpy(SmallBuffer,"scalar");
      break;
    case(3):
      strcpy(SmallBuffer,"vector");
      break;
    case(6):
      strcpy(SmallBuffer,"tensor");
      break;
    case(9):
      strcpy(SmallBuffer,"tensor9");
      break;
    }
    if (TotalTimeSteps<=1)
    {
      snprintf(charBuffer,sizeof(charBuffer),"%s per node: %s %s.%d.00000_n.%s\n",
        SmallBuffer,
        fieldName, //fileBuffer,
        this->BaseName,
        this->ProcessNumber,
        fileBuffer);
    }
    else
    {
      snprintf(charBuffer,sizeof(charBuffer),"%s per node: 1 %s %s.%d.*****_n.%s\n",
        SmallBuffer,
        fieldName, //fileBuffer,
        this->BaseName,
        this->ProcessNumber,
        fileBuffer);


    }
    this->WriteTerminatedStringToFile(charBuffer,fd);
  }

  //write the cell variable files
  for (i=0;i<input->GetCellData()->GetNumberOfArrays();i++)
  {
    //figure out what kind of data it is
    char SmallBuffer[16];

    strcpy(fileBuffer,input->GetCellData()->GetArray(i)->GetName());
    // skip arrays that were not written
    if (strcmp(fileBuffer,"GlobalElementId")==0)
    {
      continue;
    }
    if (strcmp(fileBuffer,"GlobalNodeId")==0)
    {
      continue;
    }
    if (strcmp(fileBuffer,"BlockId")==0)
    {
      continue;
    }

// fileBuffer will be sanitized, and forward slashes removed.
// forward slash is division sign in physics, thus cannot be removed!
// fieldName will be used for writing variable names, and sanitized fileName for file names
    char fieldName[512];
    strncpy(fieldName, fileBuffer, 512);

    this->SanitizeFileName(fileBuffer);
    switch(input->GetCellData()->GetArray(i)->GetNumberOfComponents())
    {
    case(1):
      strcpy(SmallBuffer,"scalar");
      break;
    case(3):
      strcpy(SmallBuffer,"vector");
      break;
    case(6):
      strcpy(SmallBuffer,"tensor");
      break;
    case(9):
      strcpy(SmallBuffer,"tensor9");
      break;
    }
    if (TotalTimeSteps<=1)
    {
      snprintf(charBuffer,sizeof(charBuffer),"%s per element: %s_c %s.%d.00000_c.%s\n",
        SmallBuffer,
        fieldName, //fileBuffer,
        this->BaseName,
        this->ProcessNumber,
        fileBuffer);
    }
    else
    {
      snprintf(charBuffer,sizeof(charBuffer),"%s per element: 1 %s_c %s.%d.*****_c.%s\n",
        SmallBuffer,
        fieldName, //fileBuffer,
        this->BaseName,
        this->ProcessNumber,
        fileBuffer);
    }
    this->WriteTerminatedStringToFile(charBuffer,fd);
  }

  // Write time block always
  // Tested for result writing in AD.
  // TO DO test for IH
  //if (TotalTimeSteps>1)
  if (TotalTimeSteps>0)
  {
    this->WriteTerminatedStringToFile("\nTIME\n",fd);
    this->WriteTerminatedStringToFile("time set: 1\n",fd);
    snprintf(charBuffer,sizeof(charBuffer),"number of steps: %d\n",TotalTimeSteps);
    this->WriteTerminatedStringToFile(charBuffer,fd);
    this->WriteTerminatedStringToFile("filename start number: 00000\n",fd);
    this->WriteTerminatedStringToFile("filename increment: 00001\n",fd);
    this->WriteTerminatedStringToFile("time values: \n",fd);
    for (i=0;i<tsteps.size();i++)
    {
      double timestep=tsteps[i];

      snprintf(charBuffer,sizeof(charBuffer),"%f ",timestep);
      this->WriteTerminatedStringToFile(charBuffer,fd);
      if (i%6==0 && i>0)
      {
        this->WriteTerminatedStringToFile("\n",fd);
      }
    }
  }
  fclose(fd);
}

vtkStandardNewMacro(cenosEnsightWriter);

void cenosEnsightWriter::WriteData()
{
  int i;
  unsigned int ui;
  int blockCount=0;
  std::list<int>::iterator iter;

  if (this->TmpInput)
  {
    this->TmpInput->Delete();
    this->TmpInput = nullptr;
  }

  //figure out process ID

  this->ProcessNumber = 0;
  this->NumberOfProcesses = 1;

#ifdef VTK_USE_PARALLEL
  vtkMultiProcessController *c = vtkMultiProcessController::GetGlobalController();

  if (c != nullptr)
  {
    this->ProcessNumber=c->GetLocalProcessId();
    this->NumberOfProcesses = c->GetNumberOfProcesses();
  }
#endif

  vtkUnstructuredGrid *input=this->GetInput();
  vtkInformation* inInfo = this->GetInputInformation();

  if (this->GhostLevel >
      inInfo->Get(vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS()))
  {
    // re-execute pipeline if necessary to obtain ghost cells

    this->GetInputAlgorithm()->UpdateInformation();
    inInfo->Set(
      vtkStreamingDemandDrivenPipeline::UPDATE_NUMBER_OF_GHOST_LEVELS(),
      this->GhostLevel);
    this->GetInputAlgorithm()->Update();
  }

  //get the BlockID Cell Array
  vtkDataArray *BlockData=input->GetCellData()->GetScalars("BlockId");

  if (BlockData==nullptr || strcmp(BlockData->GetName(),"BlockId"))
  {
    BlockData=nullptr;
  }

  this->ComputeNames();

  if (!this->BaseName)
  {
    vtkErrorMacro("A FileName or Path/BaseName must be specified.");
    return;
  }

  this->SanitizeFileName(this->BaseName);

  char** blockNames=nullptr;
  int * elementIDs=nullptr;
  char charBuffer[1024];
  char fileBuffer[512];
  snprintf(charBuffer,sizeof(charBuffer),"%s/%s.%d.%05d.geo",
    this->Path,this->BaseName,this->ProcessNumber,
    this->TimeStep);

  //open the geometry file
  //only if timestep 0 and not transient geometry or transient geometry
  FILE *fd=nullptr;
  if (this->ShouldWriteGeometry())
  {
    if (!(fd=OpenFile(charBuffer)))
      return;
  }

  //Get the FILE's for Point Data Fields
  std::vector<FILE*> pointArrayFiles;
  int NumPointArrays=input->GetPointData()->GetNumberOfArrays();
  for (i=0;i<NumPointArrays;i++)
  {
    strcpy(fileBuffer,input->GetPointData()->GetArray(i)->GetName());
    this->SanitizeFileName(fileBuffer);
    snprintf(charBuffer,sizeof(charBuffer),"%s/%s.%d.%05d_n.%s",this->Path,this->BaseName,this->ProcessNumber,
      this->TimeStep,fileBuffer);
    FILE* ftemp=OpenFile(charBuffer);
    if (!ftemp)
    {
      fclose(fd);
      return;
    }
    pointArrayFiles.push_back(ftemp);

    //write the description line to the file
    WriteStringToFile(fileBuffer,ftemp);
  }

  //Get the FILE's for Cell Data Fields
  std::vector<FILE*> cellArrayFiles;
  int NumCellArrays=input->GetCellData()->GetNumberOfArrays();
  for (i=0;i<NumCellArrays;i++)
  {

    strcpy(fileBuffer,input->GetCellData()->GetArray(i)->GetName());
    this->SanitizeFileName(fileBuffer);
    snprintf(charBuffer,sizeof(charBuffer),"%s/%s.%d.%05d_c.%s",this->Path,this->BaseName,this->ProcessNumber,
      this->TimeStep,fileBuffer);
    FILE* ftemp=OpenFile(charBuffer);
    if (!ftemp)
    {
      fclose(fd);
      return;
    }
    cellArrayFiles.push_back(ftemp);

    //write the description line to the file
    WriteStringToFile(fileBuffer,ftemp);
  }

  //write the header information
  if (this->ShouldWriteGeometry())
  {
    this->WriteStringToFile("C Binary",fd);
    this->WriteStringToFile("Written by VTK EnSight Writer",fd);
    //if (this->Title)
    // this->WriteStringToFile(this->Title,fd);
    //else
    this->WriteStringToFile("No Title was Specified",fd);
    //we will specify node and element ID's
    this->WriteStringToFile("node id given\n",fd);
    this->WriteStringToFile("element id given\n",fd);
  }

  //get the Ghost Cell Array if it exists
  vtkDataArray *GhostData=input->GetCellData()->GetScalars(vtkDataSetAttributes::GhostArrayName());
  //if the strings are not the same then we did not get the ghostData array
  if (GhostData==nullptr || strcmp(GhostData->GetName(), vtkDataSetAttributes::GhostArrayName()))
  {
    GhostData=nullptr;
  }


  //data structure to get all the cells for a certain part
  //basically sort by part# and cell type
  std::map <int, std::vector <int> > CellsByPart;

  //just a list of part numbers
  std::list<int> partNumbers;

  //get all the part numbers in the unstructured grid and sort the cells
  //by part number
  for (i=0;i<input->GetNumberOfCells();i++)
  {
    int key=1;
    if (BlockData)
      key=(int)(BlockData->GetTuple(i)[0]);
    else
      cout << "No BlockID was found\n";
    if (CellsByPart.count(key)==0)
    {
      CellsByPart[key]=std::vector < int >() ;
    }
    CellsByPart[key].push_back(i);
    partNumbers.push_back(key);

  }

  //remove the duplicates from the partNumbers
  partNumbers.sort();
  partNumbers.unique();

  //write out each part
  for (iter=partNumbers.begin();iter!=partNumbers.end();++iter)
  {
    unsigned int j;
    std::list<int>::iterator iter2;
    int part=*iter;

    //write the part Header
    if (this->ShouldWriteGeometry())
    {
      blockCount+=1;
      this->WriteStringToFile("part",fd);
      this->WriteIntToFile(part,fd);
      //this->WriteStringToFile("VTK Part", fd); // this is original line
      this->WriteStringToFile(nameList[part].c_str(),fd);  // this is replaced line
      this->WriteStringToFile("coordinates",fd);
    }

    //write the part header for data files
    for (j=0;j<pointArrayFiles.size();j++)
    {
      this->WriteStringToFile("part",pointArrayFiles[j]);
      this->WriteIntToFile(part,pointArrayFiles[j]);
      this->WriteStringToFile("coordinates",pointArrayFiles[j]);
    }
    for (j=0;j<cellArrayFiles.size();j++)
    {
      this->WriteStringToFile("part",cellArrayFiles[j]);
      this->WriteIntToFile(part,cellArrayFiles[j]);
    }

    //list of VTK Node Indices per part
    std::list<int> NodesPerPart;

    //map that goes from NodeID to the order, used for element connectivity
    std::map<int, int> NodeIdToOrder;

    //get a list of all the nodes used for a particular part
    for (j=0;j<CellsByPart[part].size();j++)
    {
      vtkCell* cell=input->GetCell(CellsByPart[part][j]);
      vtkIdList* points=cell->GetPointIds();
      for (int k=0;k<points->GetNumberOfIds();k++)
      {
        NodesPerPart.push_back(points->GetId(k));
      }
    }

    //remove the duplicate Node ID's
    NodesPerPart.sort();
    NodesPerPart.unique();

    //write the number of nodes
    if (this->ShouldWriteGeometry())
    {
      this->WriteIntToFile(static_cast<int>(NodesPerPart.size()),fd);


      //write the Node ID's to the file
      //also set up the NodeID->order map
      int NodeCount=0;
      for (iter2=NodesPerPart.begin();iter2!=NodesPerPart.end();++iter2)
      {
        this->WriteIntToFile(*iter2,fd);
        NodeIdToOrder[*iter2]=NodeCount+1;
        NodeCount++;
      }


      //EnSight format requires all the X's, then all the Y's, then all the Z's
      //write the X Coordinates

      vtkPoints* inputPoints=input->GetPoints();
      for (iter2=NodesPerPart.begin();iter2!=NodesPerPart.end();++iter2)
      {
        this->WriteFloatToFile((float)(inputPoints->GetPoint(*iter2)[0]),fd);
      }
      for (iter2=NodesPerPart.begin();iter2!=NodesPerPart.end();++iter2)
      {
        this->WriteFloatToFile((float)(inputPoints->GetPoint(*iter2)[1]),fd);
      }
      for (iter2=NodesPerPart.begin();iter2!=NodesPerPart.end();++iter2)
      {
        this->WriteFloatToFile((float)(inputPoints->GetPoint(*iter2)[2]),fd);
      }

    }

    //write the Node Data for this part
    for (j=0;j<pointArrayFiles.size();j++)
    {
      vtkDataArray* DataArray=input->GetPointData()->GetArray(j);
      //figure out what type of data it is
      int DataSize=DataArray->GetNumberOfComponents();

      for (int CurrentDimension=0;
        CurrentDimension<DataSize;
        CurrentDimension++)
      {
        for (std::list<int>::iterator k=NodesPerPart.begin();
          k!=NodesPerPart.end();++k)
        {
          this->WriteFloatToFile((float)
            (DataArray->
             GetTuple(*k)[CurrentDimension]),
            pointArrayFiles[j]);
        }
      }
    }



    //now we need to sort the cell list by element type
    //map is indexed by cell type has a vector of cell ID's
    std::map<int, std::vector<int> > CellsByElement;
    for (j=0;j<CellsByPart[part].size();j++)
    {
      int CellType=input->GetCell(CellsByPart[part][j])->GetCellType();
      int ghostLevel=0;
      if (GhostData)
      {
        ghostLevel=(int)(GhostData->GetTuple(CellsByPart[part][j])[0]);
        if (ghostLevel & vtkDataSetAttributes::DUPLICATECELL)
        {
          ghostLevel=1;
        }
      }
      //we want to sort out the ghost cells from the normal cells
      //so the element type will be ghostMultiplier*ghostLevel+elementType
      CellType+=ghostLevel*GhostLevelMultiplier;
      if (CellsByElement.count(CellType)==0)
      {
        CellsByElement[CellType]=std::vector < int >() ;
      }
      CellsByElement[CellType].push_back(CellsByPart[part][j]);
    }

    //now we need to go through each element type that EnSight understands
    std::vector<int> elementTypes;

    //list the types that EnSight understands
    //the noticeable absences are the ones without a fixed number of Nodes
    //for the ghost cell types
    elementTypes.push_back(VTK_VERTEX);
    elementTypes.push_back(VTK_LINE);
    elementTypes.push_back(VTK_TRIANGLE);
    elementTypes.push_back(VTK_QUAD);
    elementTypes.push_back(VTK_POLYGON);
    elementTypes.push_back(VTK_TETRA);
    elementTypes.push_back(VTK_HEXAHEDRON);
    elementTypes.push_back(VTK_WEDGE);
    elementTypes.push_back(VTK_PYRAMID);
    elementTypes.push_back(VTK_CONVEX_POINT_SET);
    elementTypes.push_back(VTK_QUADRATIC_EDGE);
    elementTypes.push_back(VTK_QUADRATIC_TRIANGLE);
    elementTypes.push_back(VTK_QUADRATIC_QUAD);
    elementTypes.push_back(VTK_QUADRATIC_TETRA);
    elementTypes.push_back(VTK_QUADRATIC_HEXAHEDRON);
    elementTypes.push_back(VTK_QUADRATIC_WEDGE);
    elementTypes.push_back(VTK_QUADRATIC_PYRAMID);
    elementTypes.push_back(this->GhostLevelMultiplier+VTK_VERTEX);
    elementTypes.push_back(this->GhostLevelMultiplier+VTK_LINE);
    elementTypes.push_back(this->GhostLevelMultiplier+VTK_TRIANGLE);
    elementTypes.push_back(this->GhostLevelMultiplier+VTK_QUAD);
    elementTypes.push_back(this->GhostLevelMultiplier+VTK_POLYGON);
    elementTypes.push_back(this->GhostLevelMultiplier+VTK_TETRA);
    elementTypes.push_back(this->GhostLevelMultiplier+VTK_HEXAHEDRON);
    elementTypes.push_back(this->GhostLevelMultiplier+VTK_WEDGE);
    elementTypes.push_back(this->GhostLevelMultiplier+VTK_PYRAMID);
    elementTypes.push_back(this->GhostLevelMultiplier+VTK_CONVEX_POINT_SET);
    elementTypes.push_back(this->GhostLevelMultiplier+VTK_QUADRATIC_EDGE);
    elementTypes.push_back(this->GhostLevelMultiplier+VTK_QUADRATIC_TRIANGLE);
    elementTypes.push_back(this->GhostLevelMultiplier+VTK_QUADRATIC_QUAD);
    elementTypes.push_back(this->GhostLevelMultiplier+VTK_QUADRATIC_TETRA);
    elementTypes.push_back(this->GhostLevelMultiplier+VTK_QUADRATIC_HEXAHEDRON);
    elementTypes.push_back(this->GhostLevelMultiplier+VTK_QUADRATIC_WEDGE);
    elementTypes.push_back(this->GhostLevelMultiplier+VTK_QUADRATIC_PYRAMID);

    //write out each type of element
    if (this->ShouldWriteGeometry())
    {
      for (j=0;j<elementTypes.size();j++)
      {
        unsigned int k;
        int elementType=elementTypes[j];
        if (CellsByElement.count(elementType)>0)
        {
          //switch on element type to write correct type to file
          this->WriteElementTypeToFile(elementType,fd);

          //number of elements
          this->WriteIntToFile(
            static_cast<int>(CellsByElement[elementType].size()),fd);

          //element ID's
          for (k=0;k<CellsByElement[elementType].size();k++)
          {
            int CellId=CellsByElement[elementType][k];
            this->WriteIntToFile(CellId,fd);
          }

          //element conenctivity information
          for (k=0;k<CellsByElement[elementType].size();k++)
          {
            int CellId=CellsByElement[elementType][k];
            vtkIdList *PointIds=input->GetCell(CellId)->GetPointIds();
            for (int m=0;m<PointIds->GetNumberOfIds();m++)
            {
              int PointId=PointIds->GetId(m);
              this->WriteIntToFile(NodeIdToOrder[PointId],fd);
            }
          }
        }
      }
    }
    
    //write the Cell Data for this part
    for (j=0;j<cellArrayFiles.size();j++)
    {
      vtkDataArray* DataArray=input->GetCellData()->GetArray(j);
      //figure out what type of data it is
      int DataSize=DataArray->GetNumberOfComponents();

      for (unsigned int k=0;k<elementTypes.size();k++)
      {
        if (!CellsByElement[elementTypes[k]].empty())
        {
          this->WriteElementTypeToFile(elementTypes[k],
            cellArrayFiles[j]);
          for (unsigned int m=0;m<CellsByElement[elementTypes[k]].size();m++)
          {
            for ( int CurrentDimension = 0; CurrentDimension < DataSize;
                     CurrentDimension ++ )
            {
              this->WriteFloatToFile
                (  (float)
                   (   DataArray
                       ->GetTuple(  CellsByElement[ elementTypes[k] ][m]  )
                         [CurrentDimension]
                   ),
                   cellArrayFiles[j]
                );
            }
          }
        }
      }
    }
  }

  //now write the empty blocks
  //use the block list in the exodus model if it exists, otherwise
  //use the BlockID list if that exists.

  if (this->BlockIDs)
  {
    elementIDs=this->BlockIDs;
  }

  if (elementIDs)
  {
    for (i=0;i<this->NumberOfBlocks;i++)
    {
      unsigned int j;
      //figure out if the part was already written
      int part=elementIDs[i];
      if (   std::find(partNumbers.begin(), partNumbers.end(), part)
        == partNumbers.end() )
      {
        //no information about the part was written to the output files
        //so write some empty information
        if (this->ShouldWriteGeometry())
        {
          blockCount+=1;
          this->WriteStringToFile("part",fd);
          this->WriteIntToFile(part,fd);

          int exodusIndex =
            this->GetExodusModelIndex(elementIDs,this->NumberOfBlocks,part);

          if (exodusIndex!=-1 && blockNames)
          {
            snprintf(charBuffer,sizeof(charBuffer),"Exodus-%s-%d",blockNames[exodusIndex],part);
            this->WriteStringToFile(charBuffer,fd);
          }
          else
          {
            this->WriteStringToFile("VTK Part",fd);
          }
        }

        //write the part header for data files
        for (j=0;j<pointArrayFiles.size();j++)
        {
          this->WriteStringToFile("part",pointArrayFiles[j]);
          this->WriteIntToFile(part,pointArrayFiles[j]);
        }
        for (j=0;j<cellArrayFiles.size();j++)
        {
          this->WriteStringToFile("part",cellArrayFiles[j]);
          this->WriteIntToFile(part,cellArrayFiles[j]);
        }
      }
    }
  }
  
  if (this->TmpInput)
  {
    this->TmpInput->Delete();
    this->TmpInput = nullptr;
  }

  
  //close all the files
  if (fd)
  {
    fclose(fd);
  }

   for (ui = 0; ui < cellArrayFiles.size(); ui++)
   {
       fclose(cellArrayFiles[ui]);
   }

   for (ui = 0; ui < pointArrayFiles.size(); ui++)
   {
       fclose(pointArrayFiles[ui]);
   }
}

void cenosEnsightWriter::WriteStringToFile(const char* cstring, FILE* file)
{
    char cbuffer[81];

    unsigned long cstringLength = static_cast<unsigned long>(strlen(cstring));
    auto ts = vtkMath::Min(cstringLength, 80ul);
    memcpy(cbuffer, cstring, vtkMath::Min(cstringLength, 80ul));
    for (int i = cstringLength; i < 81; ++i)
    {
        cbuffer[i] = '\0';
    }

    // Write a constant 80 bytes to the file.
    fwrite(cbuffer, sizeof(char), 80, file);
}

void cenosEnsightWriter::WriteTerminatedStringToFile(const char* cstring, FILE* file)
{
    fwrite(cstring, sizeof(char), std::min(strlen(cstring), static_cast<size_t>(512)), file);
}


void cenosEnsightWriter::setNameList(std::vector<int> ids, std::vector<std::string> blockName)
{
	for (int i=0; i< ids.size(); i++)
	{
		nameList[ids[i]] = blockName[i];
	}
}

FILE* cenosEnsightWriter::OpenFile(char* name)
{
    FILE* fd = fopen(name, "wb");

    if (fd == nullptr)
    {
        vtkErrorMacro("Error opening " << name
            << ": " << strerror(errno));
        return nullptr;
    }
    return fd;
}

void cenosEnsightWriter::WriteIntToFile(const int i, FILE * file)
{
   fwrite(&i, sizeof(int), 1, file);
}

void cenosEnsightWriter::WriteFloatToFile(const float f, FILE * file)
{
  fwrite(&f, sizeof(float), 1, file);
}