#include "Mesh.h"
#include "Vector.h"
#include <iostream>
#include <map>
#include "ccmioutility.h"

using namespace std;

const int kMaxErrors = 10;

bool Mesh::FindFace(int id, FaceArray::iterator& fit)
{
    for (fit = this->internalFaces.begin();  fit != this->internalFaces.end();
									  ++fit)
	if (fit->id == id)
	    return(true);
    for (fit = this->boundaryFaces.begin();  fit != this->boundaryFaces.end();
									  ++fit)
	if (fit->id == id)
	    return(true);
    return(false);
}

bool Mesh::CompareInternalFace( const Mesh& mesh1, const MeshFace& face1,
				const Mesh& mesh2, const MeshFace& face2 )
{
    if (!mesh1.noCheck && !mesh2.noCheck &&
	(face1.cells[0] != face2.cells[0] ||
	 face1.cells[1] != face2.cells[1]))
	return(false);
    return(CompareFaceVertices(mesh1, face1, mesh2, face2));
}

bool Mesh::CompareBoundaryFace( const Mesh& mesh1, const MeshFace& face1,
				const Mesh& mesh2, const MeshFace& face2 )
{
    if (!mesh1.noCheck && !mesh2.noCheck &&
	(face1.cells[0] != face2.cells[0] /*||
         face1.boundary != face2.boundary*/))
	return(false);
    return(CompareFaceVertices(mesh1, face1, mesh2, face2));
}

bool Mesh::CompareFaceVertices( const Mesh& mesh1, const MeshFace& face1,
				const Mesh& mesh2, const MeshFace& face2 )
{
    int err = 0, index1, index2;
    const VertexArray& v1 = mesh1.vertices;
    const VertexArray& v2 = mesh2.vertices;

    if (face1.nVerts == face2.nVerts)
    {
	// Not only do the coordinates of the face matters, but the
	// order matters, too, since the direction of the winding should point
	// a certain way.
	for (int i = 0;  i < face1.nVerts;  ++i)
	{
	    index1 = face1.verts[i];
	    index2 = face2.verts[i];
	    if (v1[index1].coord[0] != v2[index2].coord[0] ||
		v1[index1].coord[1] != v2[index2].coord[1]||
		v1[index1].coord[2] != v2[index2].coord[2])
		err++;
	}
    }
    else
	err++;
    return(!err);
}

void Mesh::PrintFaceError(
    			   const Mesh& mesh1, const MeshFace& face1,
			   const Mesh& mesh2, const MeshFace& face2,
			   bool isInternal )
{
    if (isInternal)
	cout << "Internal face " << face1.id << "/" << face2.id << " differs:"
	     << endl;
    else
	cout << "Boundary face differs:" << endl;;
    cout << "Face1:  ";
    if (isInternal)
	cout << "cells { " << face1.cells[0] << ", "
	     << face1.cells[1] << " }     "
	     << face1.nVerts << " vertices";
    else
	cout << "cell " << face1.cells[0] << "   boundary "
	     << face1.boundary << "     "
	     << face1.nVerts << " vertices";
    cout << endl;
    for (int i = 0;  i < face1.nVerts;  ++i)
	cout << "\t" << face1.verts[i] << "\t("
	     << mesh1.vertices[face1.verts[i]].coord[0] << ", "
	     << mesh1.vertices[face1.verts[i]].coord[1] << ", "
	     << mesh1.vertices[face1.verts[i]].coord[2] << ") " << endl;

    cout << "Face2:  ";
    if (isInternal)
	cout << "cells { " << face2.cells[0] << ", "
	     << face2.cells[1] << " }     "
	     << face2.nVerts << " vertices";
    else
	cout << "cell " << face2.cells[0] << "   boundary "
	     << face2.boundary << "     "
	     << face2.nVerts << " vertices";
    cout << endl;
    for (int i = 0;  i < face2.nVerts;  ++i)
	cout << "\t" << face2.verts[i] << "\t("
	     << mesh2.vertices[face2.verts[i]].coord[0] << ", "
	     << mesh2.vertices[face2.verts[i]].coord[1] << ", "
	     << mesh2.vertices[face2.verts[i]].coord[2] << ") " << endl;
}

void Mesh::PrintInfo()
{
    cout << "\t" << this->nVertices << " vertices, " << this->nCells
	 << " cells, " << this->internalFaces.size() << " internal faces, "
	 << endl << "\t" << this->nBoundaries << " boundaries, "
	 << this->boundaryFaces.size() << " boundary faces." << endl;
}

bool Mesh::operator== ( const Mesh& that )
{
    if (this->err != kCCMIONoErr)
    {
	cout << "First file has error " << this->err << endl;
	return(false);
    }
    if (that.err != kCCMIONoErr)
    {
	cout << "Second file has error " << that.err << endl;
	return(false);
    }

    // Can only compare the number of cells because CP matches cause
    // additional faces and vertices to be created.  Also, internal faces
    // may become boundary faces if meshes are split into domains/subdomains.
    if (::abs(this->nCells - that.nCells) == 1)
	cout << "Warning:  cell count differs by 1 (" << this->nCells
	     << " to " << that.nCells
	     << "), possibly due to polyfile API bug" << endl;
    else if (this->nCells != that.nCells)
    {
	cout << "Warning:  Differing number of cells (" << this->nCells
	     << " to " << that.nCells << ")" << endl;
    }

    // The face and vertex IDs will probably be different between meshes,
    // but let's assume that the cell ID is the same.  Compare the meshes
    // by determining the cell's vertices.
    CellVertsArray cellVerts1, cellVerts2;
    CellVertsArray::iterator cell1, cell2;
    map< int, int >::iterator vi1, vi2;
    int nErrors = 0, cellNum = 0;

    CreateCellVertsArray(*this, cellVerts1);
    CreateCellVertsArray(that, cellVerts2);

    cell1 = cellVerts1.begin();
    cell2 = cellVerts2.begin();
    for(; cell1 != cellVerts1.end();  ++cell1, ++cell2, ++cellNum)
    {
	bool err = false;

	if (cell1->size() != cell2->size())
	    { err = true; goto error; }
	vi1 = cell1->begin();
	vi2 = cell2->begin();
	for (;  vi1 != cell1->end();  ++vi1, ++vi2)
	{
	    int v1 = vi1->second, v2 = vi2->second;
	    if (this->vertices[v1].coord[0] != that.vertices[v2].coord[0] ||
		this->vertices[v1].coord[1] != that.vertices[v2].coord[1] ||
		this->vertices[v1].coord[2] != that.vertices[v2].coord[2])
		{ err = true; goto error; }
	}
error:
	if (err && nErrors <= kMaxErrors)
	{
	    nErrors++;
	    if (nErrors == kMaxErrors)
	    {
		cout << "\t... more errors ..." << endl;
		continue;
	    }
	    cout << "Cell " << cellNum << " differs:" << endl;
	    vi1 = cell1->begin();
	    vi2 = cell2->begin();
	    for (;  vi1 != cell1->end() || vi2 != cell2->end(); )
	    {
		int v1 = vi1->second, v2 = vi2->second;
		if (vi1 != cell1->end())
		    cout << "\t" << v1 << "\t(" << this->vertices[v1].coord[0]
			 << ", " << this->vertices[v1].coord[1] << ", "
			 << this->vertices[v1].coord[2] << ")\t";
		else
		    cout << "\t\t\t\t\t";
		if (vi2 != cell2->end())
		    cout << "\t" << v1 << "\t(" << this->vertices[v2].coord[0]
			 << ", " << this->vertices[v2].coord[1] << ", "
			 << this->vertices[v2].coord[2] << ")";
		cout << endl;

		if (vi1 != cell1->end())
		    ++vi1;
		if (vi2 != cell2->end())
		    ++vi2;
	    }
	}
    }

    FaceArray::const_iterator fit1, fit2;
    
    if (this->boundaryFaces.size() != that.boundaryFaces.size())
    {
	cout << "Mesh " << this->filename << " has "
	     << this->boundaryFaces.size() << " boundary faces, but mesh "
	     << that.filename << " has " << that.boundaryFaces.size()
	     << "." << endl;
	return(false);
    }

    for (fit1 = this->boundaryFaces.begin(), fit2 = that.boundaryFaces.begin();
	 fit1 != this->boundaryFaces.end();  ++fit1, ++fit2)
    {
	if (!CompareBoundaryFace(*this, *fit1, that, *fit2))
	{
	    nErrors++;
	    if (nErrors == kMaxErrors)
	    {
		cout << "\t... more errors ..." << endl;
		continue;
	    }
	    else if (nErrors < kMaxErrors)
		PrintFaceError(*this, *fit1, that, *fit2, false);
	}
    }
    
    if (this->internalFaces.size() != that.internalFaces.size())
    {
	cout << "Mesh " << this->filename << " has "
	     << this->internalFaces.size() <<" internal faces, but mesh "
	     << that.filename << " has " << that.internalFaces.size()
	     << "." << endl;
	return(false);
    }

    for (fit1 = this->internalFaces.begin(), fit2 = that.internalFaces.begin();
	 fit1 != this->internalFaces.end();  ++fit1, ++fit2)
    {
	if (!CompareInternalFace(*this, *fit1, that, *fit2))
	{
	    nErrors++;
	    if (nErrors == kMaxErrors)
	    {
		cout << "\t... more errors ..." << endl;
		continue;
	    }
	    else if (nErrors < kMaxErrors)
		PrintFaceError(*this, *fit1, that, *fit2, true);
	}
    }
    
    return(nErrors == 0);
}

void Mesh::CheckFaceOrientation(CellVertsArray& cells)
{
    CellVertsArray::iterator cit;
    map<int, int>::iterator vit;
    vector<Vector> centroids;
    Vector v;

    centroids.reserve(cells.size());

    for (cit = cells.begin();  cit != cells.end();  ++cit)
    {
	float size;
	
	size = (float)cit->size();
	for (vit = cit->begin();  vit != cit->end();  ++vit)
	{
	    v[0] += this->vertices[vit->first].coord[0];
	    v[1] += this->vertices[vit->first].coord[1];
	    v[2] += this->vertices[vit->first].coord[2];
	}
	v[0] /= size;
	v[1] /= size;
	v[2] /= size;
	centroids.push_back(v);
    }

    FaceArray::iterator f;
    for (f = this->boundaryFaces.begin(); f != this->boundaryFaces.end(); ++f)
    {
	Vector v0(this->vertices[f->verts[0]].coord);
	Vector v1(this->vertices[f->verts[1]].coord);
	Vector v2(this->vertices[f->verts[f->nVerts - 1]].coord);
	Vector v, w, axis, center;
	
	v = v1 - v0;
	w = v2 - v0;
	axis = Vector::cross(v, w);
	center = centroids[f->cells[0]] - v0;
	if (Vector::dot(axis, center) > 0)	// axis should point outside
	{
	    cout << "Boundary " << f->boundary << ", face " << f->id
		 << " (cell " << f->cells[0]
		 << "):  face vertices are reversed." << endl;
	    cout << "\tv1: " << v1 << endl;
	    cout << "\tv2: " << v2 << endl;
	    cout << "\tcentroid: " << centroids[f->cells[0]] << endl;
	    cout << "\tv: " << v << endl;
	    cout << "\tw: " << w << endl;
	    cout << "\taxis: " << axis << endl;;
	    cout << "\tcenter: " << center << endl;
	    cout << "dot: " << Vector::dot(axis, center) << endl;
	}
    }

}

bool Mesh::CreateCellVertsArray( const Mesh& mesh, CellVertsArray& cells )
{
    int nCells = mesh.cells.size();
    FaceArray::const_iterator f;

    cells.clear();
    cells.resize(nCells + 1);
    for (f = mesh.internalFaces.begin();  f != mesh.internalFaces.end();  ++f)
    {
	if (f->cells[0] > nCells)
	    cout << "Internal error (file1):  not enough cells allocated!"
		 << endl;
	if (f->cells[1] > nCells)
	    cout << "Internal error (file2):  not enough cells allocated!"
		 << endl;
	for (int i = 0;  i < f->nVerts;  ++i)
	{
	    cells[f->cells[0]][f->verts[i]] = f->verts[i];
	    cells[f->cells[1]][f->verts[i]] = f->verts[i];
	}
    }

    for (f = mesh.boundaryFaces.begin();  f != mesh.boundaryFaces.end();  ++f)
    {
	for (int i = 0;  i < f->nVerts;  ++i)
	    cells[f->cells[0]][f->verts[i]] = f->verts[i];
    }

    return(true);
}

void Mesh::OffsetIDs( int offset )
{
    VertexArray::iterator v;
    for (v = this->vertices.begin();  v != this->vertices.end();  ++v)
	v->id += offset;
    for (v = this->packedVerts.begin();  v != this->packedVerts.end();  ++v)
	v->id += offset;

    CellArray::iterator c;
    for (c = this->cells.begin();  c != this->cells.end();  ++c)
	c->id += offset;

    FaceArray::iterator f;
    for (f = this->internalFaces.begin();  f != this->internalFaces.end();  ++f)
    {
	f->id += offset;
	for (int i = 0;  i < f->nVerts;  ++i)
	    f->verts[i] += offset;
	f->cells[0] += offset;
	f->cells[1] += offset;
    }

    for (f = this->boundaryFaces.begin();  f != this->boundaryFaces.end();  ++f)
    {
	f->id += offset;
	for (int i = 0;  i < f->nVerts;  ++i)
	    f->verts[i] += offset;
	f->cells[0] += offset;
	f->boundary += offset;
    }
}

template <class T>
CCMIOID Mesh::CreateIndexMap( vector<T>& items, CCMIOID& entity, CCMIOEntity type,
			    const char *description, bool isCell /*= false*/ )
{
    int max = -1;
    CCMIOID mapID;
    vector<int> ids;

    ids.reserve(items.size());
    for (unsigned int i = 0;  i < items.size();  ++i)
    {
	if (!isCell || (isCell && this->cellUsed[i]))
	{
	    ids.push_back(items[i].id);
	    if (max < items[i].id)
		max = items[i].id;
	}
    }

    CCMIONewEntity(&this->err, entity, type, description, &mapID);
    CCMIOWriteMap(&this->err, mapID, ids.size(), max, &ids[0],
		  kCCMIOStart, kCCMIOEnd);
    return(mapID);
}

void Mesh::WriteCCMIO( CCMIOID& state, CCMIOID& processor )
{
    CCMIOID verticesID, topologyID, mapID, id;
    CCMIOError *err = &this->err;
    map<int, bool> cellTypesUsed;

    CCMIOClearProcessor(err, state, processor, TRUE, TRUE, TRUE, TRUE, TRUE);

    // Write the vertices
    cout << "Writing " << this->packedVerts.size() << " vertices." << endl;
    mapID = CreateIndexMap<MeshVertex>(this->packedVerts, processor, kCCMIOMap,
				       "Vertex map");
    CCMIONewEntity(err, processor, kCCMIOVertices, "Vertices", &verticesID);
    {
	vector<float> v(this->packedVerts.size() * 3, 0);
	for (unsigned int i = 0;  i < this->packedVerts.size();  ++i)
	{
	    v[3 * i    ] = this->packedVerts[i].coord[0];
	    v[3 * i + 1] = this->packedVerts[i].coord[1];
	    v[3 * i + 2] = this->packedVerts[i].coord[2];
	}
	CCMIOWriteVerticesf(err, verticesID, 3, 1.0, mapID, &v[0],
			    kCCMIOStart, kCCMIOEnd);
    }

    // Write the cells
    cout << "Writing " << this->nCells << " cells." << endl;
    mapID = CreateIndexMap<MeshCell>(this->cells, processor, kCCMIOMap,
				     "Cell map", true);
    CCMIONewEntity(err, processor, kCCMIOTopology, "Topology", &topologyID);
    CCMIONewEntity(err, topologyID, kCCMIOCells, "Cells", &id);
    {
	unsigned int size = this->cells.size();
	vector<int> types;
	types.reserve(this->cells.size());
	for (unsigned int i = 0;  i < size;  ++i)
	{
	    if (this->cellUsed[i])
		types.push_back(this->cells[i].type);
	    cellTypesUsed[this->cells[i].type] = true;
	}
	CCMIOWriteCells(err, id, mapID, &types[0], kCCMIOStart, kCCMIOEnd);
    }

    // Write the faces
    FaceArray::iterator f;
    cout << "Writing " << this->internalFaces.size() << " internal faces." << endl;
    if (this->internalFaces.size() != 0)
    {
	vector<int> v, c;

	mapID = CreateIndexMap(this->internalFaces, processor, kCCMIOMap,
			       "Internal face map");
	CCMIONewEntity(err, topologyID, kCCMIOInternalFaces, "Internal faces", &id);
	for (f=this->internalFaces.begin(); f != this->internalFaces.end(); ++f)
	{
	    v.push_back(f->nVerts);
	    for (int i = 0;  i < f->nVerts;  ++i)
		v.push_back(f->verts[i]);
	    c.push_back(f->cells[0]);
	    c.push_back(f->cells[1]);
	}
	CCMIOWriteFaces(err, id, kCCMIOInternalFaces, mapID,
			v.size(), &v[0], kCCMIOStart, kCCMIOEnd);
	CCMIOWriteFaceCells(err, id, kCCMIOInternalFaces, mapID,
			    &c[0], kCCMIOStart, kCCMIOEnd);
    }
	
    // The boundary faces really need to be sorted into separate lists based
    // on boundary ID, but this is just a little test program, so we'll do
    // it the quick and inefficient way.
    cout << "Writing " << this->boundaryFaces.size() << " boundary faces." << endl;
    {
	int curID;
	vector<int> v, c;
	map<int, int> boundaryList;
	for (f = this->boundaryFaces.begin();f != this->boundaryFaces.end();++f)
	    boundaryList[f->boundary] = f->boundary;
	
	for (map<int, int>::iterator iter = boundaryList.begin();
	     iter != boundaryList.end();  ++iter)
	{
	    int pos = 0, max;
	    vector<int> ids;

	    v.resize(0);
	    c.resize(0);
	    curID = iter->first;
	    for (f=this->boundaryFaces.begin(); f != this->boundaryFaces.end();
									    ++f)
	    {
		if (f->boundary == curID)
		{
		    ids.push_back(f->id);
		    max = pos;
		    v.push_back(f->nVerts);
		    for (int i = 0;  i < f->nVerts;  ++i)
			v.push_back(f->verts[i]);
		    c.push_back(f->cells[0]);
		}
		pos++;
	    }
	    CCMIONewEntity(err, processor, kCCMIOMap, "Boundary faces", &mapID);
	    CCMIOWriteMap(err, mapID, ids.size(), max, &ids[0],
			  kCCMIOStart, kCCMIOEnd);
	    CCMIONewIndexedEntity(err, topologyID, kCCMIOBoundaryFaces,
				iter->first, "Boundary faces", &id);
	    CCMIOWriteFaces(err, id, kCCMIOBoundaryFaces, mapID, v.size(),
			    &v[0], kCCMIOStart, kCCMIOEnd);
	    CCMIOWriteFaceCells(err, id, kCCMIOBoundaryFaces, mapID,
				&c[0], kCCMIOStart, kCCMIOEnd);
	}
    }

    // Write out a dummy problem description
    CCMIOID problem;
    CCMIONode node;
    map<int, bool>::iterator iter;

    CCMIONewEntity(err, state, kCCMIOProblemDescription, "Dummy description",
		 &problem);
    for (iter = cellTypesUsed.begin();  iter != cellTypesUsed.end();  ++iter)
    {
	CCMIONewIndexedEntity(err, problem, kCCMIOCellType, iter->first,
			    "Dummy celltypes", &id);
	CCMIOGetEntityNode(err, id, &node);
	CCMIOWriteNodestr(err, node, "MaterialType", "plasma");
    }
    CCMIOWriteState(err, state, problem, NULL);

    // Write out the processor
    CCMIOWriteProcessor(err, processor, NULL, &verticesID, NULL, &topologyID,
		      NULL, NULL, NULL, NULL);
}
