/*==========================================================================

  Program:   Reading and Visualizing FEM data with VTK
  Fraunhofer IWU: Sandra Scherer

===========================================================================*/

#include "vtkDoubleArray.h"
#include "vtkPointData.h"
#include "vtkPoints.h"
#include "vtkUnstructuredGrid.h"
#include "vtkPoints.h"
#include "vtkTetra.h"
#include "vtkHexahedron.h"
#include "vtkLookupTable.h"
#include "vtkScalarBarActor.h"
#include "vtkScalarsToColors.h"
#include "vtkCellArray.h"
#include "vtkSmartPointer.h"

#include "ansysFEMresults.h"
#include "stdio.h"
#include <iostream>
#include <fstream>
#include <string>
#include <vector>

using namespace std;

//NR ifstream Referenz, Zuweisung durch Initialisierungsliste
FEresult::FEresult(std::ifstream &file) : istr(file) {
	
	nodeCoordinates = vtkPoints::New();
	uGrid = vtkUnstructuredGrid::New();
	lut = vtkLookupTable::New();
	scalars = vtkDoubleArray::New();

	readInputFile();
	buildGeometry(elementType);	
	buildLookUpTable();	
}


void FEresult::readInputFile()
{
	istr>>d;
	numNodes=(long)d;
	setNumNodes(numNodes);	
	numbers.resize(numNodes);
	
	istr>>d;
	numCells=(long)d;
	setNumCells(numCells);
	cells = new int(numCells);
	hexVector.resize(numCells);	
	tetraVector.resize(numCells);
	uGrid->Allocate(numCells); //NR allocate verschoben, numCells im Konstruktor noch nicht bekannt, daher Programmabsturz

	istr>>d;
	elementType=(long)d;
	setElementType(elementType);

	cerr << "Reading ANSYS Input File..." << endl;
	cerr << "Number of Nodes:  " << FEresult::getNumNodes()<< endl;
	cerr << "Number of Cells:  " << FEresult::getNumCells() << endl;
	cerr << "ANSYS Element type: " << FEresult::getElementType() << endl;
	
	n=0;
	stop=0;	
	
	float* xCoordinates = new float[numNodes];
	float* yCoordinates = new float[numNodes];
    float* zCoordinates = new float[numNodes];
	
	double x, y, z;	
	int* numbers = new int[numNodes];

	cerr << "Reading node coordinates..." << endl;

	while ((!stop) && (n < numNodes))
	{
		istr>>nodeNumber;
		//NR STL map um Lücken zu schließen
		pointIndexLUT.insert(lmap::value_type(nodeNumber, n));
		
	// read the coordinates and build the next node
		istr >> x;		
		istr >> y;
		istr >> z;
	
		xCoordinates[n] = (float)x; yCoordinates[n] = (float)y; zCoordinates[n] = (float)z;

		numbers[n] = n; //NR Laufindex statt nodeNumber
		
		nodeCoordinates->InsertPoint(n, x, y, z);
		
        ++n;		
	}	

	std::cout << "number of nodes / vtkPoints size: " << numNodes << "/" << nodeCoordinates->GetNumberOfPoints() << std::endl;
	

}			

void FEresult::buildGeometry(int elementType)
{
	if (elementType == 185)
	{
		buildHexa();
	}
	else if (elementType ==187)
	{
		buildTetra();
	}
}

//Element type is 8-node hexaeder
void FEresult::buildHexa()
{

	cerr << "Build hexaeder mesh..." << endl;
	long index;
	for(long i = 0; i < numCells; i++)
	{
		vtkSmartPointer<vtkHexahedron> theHex = vtkSmartPointer<vtkHexahedron>::New();	
		hexa = new hexaCell();

		istr>>cellNumber;	

		//NR wie bei den Punkten, map zur Erhaltung der Reihenfolge
		hexaIndexLUT.insert(lmap::value_type(cellNumber,i));

		hexa->id = i; 

		theHex->GetPointIds()->SetNumberOfIds(8);  

		istr>>index;
		index = (pointIndexLUT.find(index))->second; //NR suche Knoten in map und gibt neuen Index zurück, analog für die anderen Knoten
		theHex->GetPointIds()->SetId(4,index); 
		hexa->node1 = index;

		istr>>index;
		index = (pointIndexLUT.find(index))->second;
		theHex->GetPointIds()->SetId(5,index);
		hexa->node2 = index;

		istr>>index;	
		index = (pointIndexLUT.find(index))->second;
		theHex->GetPointIds()->SetId(1,index); 
		hexa->node3 = index;

		istr>>index;	
		index = (pointIndexLUT.find(index))->second;
		theHex->GetPointIds()->SetId(0,index); 
		hexa->node4 = index;

		istr>>index;	
		index = (pointIndexLUT.find(index))->second;
		theHex->GetPointIds()->SetId(7,index); 
		hexa->node5 = index;

		istr>>index;	
		index = (pointIndexLUT.find(index))->second;
		theHex->GetPointIds()->SetId(6,index); 
		hexa->node6 = index;

		istr>>index;	
		index = (pointIndexLUT.find(index))->second;
		theHex->GetPointIds()->SetId(2,index);
		hexa->node7 = index;

		istr>>index;	
		index = (pointIndexLUT.find(index))->second;
		theHex->GetPointIds()->SetId(3,index); 
		hexa->node8 = index;

		hexVector.at(i) = hexa;		
	
		vtkSmartPointer<vtkCellArray> hexas = vtkSmartPointer<vtkCellArray>::New();
		hexas->InsertNextCell(theHex);

		uGrid->InsertNextCell(theHex->GetCellType(), theHex->GetPointIds());
	}
	std::cout << "hexa element map size: " << hexaIndexLUT.size() << std::endl;
}

//-------------------------------------------------------------------------------------------

//Element type is 4-node tetraeder
void FEresult::buildTetra()
{
			cerr << "Build tetraeder mesh..." << endl;
			int index;
			for(int i = 0; i < numCells; i++)
			{
				istr>>cellNumber;
				tetraIndexLUT.insert(lmap::value_type(cellNumber, i)); //NR das gleiche wie bei Hexa
		
				tetra = new tetraCell;
				tetra->id = i;	

				vtkSmartPointer<vtkTetra> theTetra = vtkSmartPointer<vtkTetra>::New();		
				theTetra->GetPointIds()->SetNumberOfIds(4);  

				istr>>index;
				index = (pointIndexLUT.find(index))->second;
				theTetra->GetPointIds()->SetId(0,index); 
				tetra->node1 = index;

				istr>>index;
				index = (pointIndexLUT.find(index))->second;
				theTetra->GetPointIds()->SetId(1,index);
				tetra->node2 = index;

				istr>>index;	
				index = (pointIndexLUT.find(index))->second;
				theTetra->GetPointIds()->SetId(2,index); 
				tetra->node3 = index;

				istr>>index;	
				index = (pointIndexLUT.find(index))->second;
				theTetra->GetPointIds()->SetId(3,index); 
				tetra->node4 = index;

				vtkSmartPointer<vtkCellArray> tetras = vtkSmartPointer<vtkCellArray>::New();
				tetras->InsertNextCell(theTetra);

				tetraVector.at(i) = tetra;
				 
				uGrid->InsertNextCell(theTetra->GetCellType(), theTetra->GetPointIds());				
			}
			std::cout << "tetra element map size: " << tetraIndexLUT.size() << std::endl;
}

void FEresult::buildLookUpTable()
{ 
	cerr << "Build LookUpTable... " << endl;
	lut->SetNumberOfTableValues(9);  

	// Ansys color table
	lut->SetTableValue(0,	0,	0,	1, 1);  
	lut->SetTableValue(1,	0,	0.7,1, 1); 
	lut->SetTableValue(2,	0,	1,	1, 1); 
	lut->SetTableValue(3,	0,	1,	0.7, 1); 
	lut->SetTableValue(4,	0,	1,	0, 1); 
	lut->SetTableValue(5,	0.7,1,	0, 1); 
	lut->SetTableValue(6,	1,	1,	0, 1); 
	lut->SetTableValue(7,	1,	0.7,0, 1); 
	lut->SetTableValue(8,	1,	0,	0, 1); 

	lut->SetRange(0,0.1);	
	lut->Build();
}

void FEresult::buildScalars()
{	
	vtkIdType numberOfPoints;	
    numberOfPoints = uGrid->GetNumberOfPoints();
    cerr << "number of grid points:" << numberOfPoints << endl;

	istr>>numStressValues;
		
	scalars->SetNumberOfValues(numberOfPoints);
	for (vtkIdType i = 0; i < scalars->GetNumberOfTuples(); ++i) {
		scalars->SetValue(i, 0.0);
	}


	//NR ein if unterbricht den Programmablauf, daher ifs in Schleifen minimieren und nach außen ziehen
	if (elementType == 187) {
		for (n = 0; n < numStressValues; ++n) {
			istr>>cellNr;

			//NR Suche Zellennummer in map, nutzt 0 wenn Zellennummer nicht vorhanden z.B. Scalarwert
			lmap::iterator it = tetraIndexLUT.find(cellNr);
			if (it != tetraIndexLUT.end())
				cellNr = it->second;
			else {
				std::cerr << "cell number not found("<<cellNr<<"), fall back to zero" << std::endl;
				cellNr = 0;
			}
			setCellNumber(cellNr);
			istr>>value;
			buildTetraScalars();
		}
	}
	if (elementType == 185) {
		for (n = 0; n < numStressValues; ++n) {
			istr>>cellNr;
			lmap::iterator it = hexaIndexLUT.find(cellNr);
			if (it != hexaIndexLUT.end())
				cellNr = it->second;
			else {
				//std::cerr << "cell number not found("<<cellNr<<"), fall back to zero" << std::endl;
				cellNr = 0;
			}
			setCellNumber(cellNr);
			istr>>value;
			buildHexaScalars();
		}
	}	
}

void FEresult::buildHexaScalars()
{				
					long k = cellNr;
					scalars->SetValue(hexVector.at(k)->node1,value);	
					scalars->SetValue(hexVector.at(k)->node2,value);
					scalars->SetValue(hexVector.at(k)->node3,value);
					scalars->SetValue(hexVector.at(k)->node4,value);	
					scalars->SetValue(hexVector.at(k)->node5,value);
					scalars->SetValue(hexVector.at(k)->node6,value);
					scalars->SetValue(hexVector.at(k)->node7,value);
					scalars->SetValue(hexVector.at(k)->node8,value);	
}

void FEresult::buildTetraScalars()
{
	//NR durch die Schließung der Lücken kann der Index direkt benutzt anstatt erst gesucht werden
	//NR netter Nebeneffekt: nun ist das Programm ca. 20 mal schneller
			long k = cellNr;
			scalars->SetValue(tetraVector.at(k)->node1,value);	
			scalars->SetValue(tetraVector.at(k)->node2,value);
			scalars->SetValue(tetraVector.at(k)->node3,value);
			scalars->SetValue(tetraVector.at(k)->node4,value);		
}

void FEresult::setNumNodes(int numberOfNodes)
{
	numNodes = numberOfNodes;
}

 void FEresult::setNumCells(int numberOfCells)
{
	numCells = numberOfCells;
}

 void FEresult::setElementType(int element){

	elementType = element;
}

 void FEresult::setCellNumber(int num)
{
	cellNr = num;
}

 int FEresult::getCellNumber()
{
	return cellNr;
}
long FEresult::getNumNodes()
{
	return numNodes;
}

long FEresult::getNumCells()
{
	return numCells;
}

long FEresult::getElementType()
{
	 return elementType;
}