#ifndef POST_CPP
#define POST_CPP

#include "Post.h"
#include <math.h>
#include <iostream>

using namespace std;

const int kMaxError = 10;
const char kFieldSetName[] = "Default field set";

class Int
{
 public:
    int val;

    Int() : val(0) {};
    inline bool operator==(int that) const { return(this->val == that); };
    inline bool operator!=(int that) const { return(!(*this == that)); };
    inline void operator++(int) { this->val++; };
};

std::ostream& operator<<(std::ostream& out, BoundaryData& data)
{
    cout << "{ " << data.id << ", " << data.priority << ", (";
    for (int j = 0;  j < 9;  ++j)
	cout << data.value[j] << ", ";
    cout << ")}";
    return(out);
};

bool Post::operator== (const Post& that) const
{
    int nErrors = 0;
    FieldArray::const_iterator f1, f2;

    if (this->fields.size() != that.fields.size())
    {
	cout << "Files do not have the same number of data fields ("
	     << this->fields.size() << " and " << that.fields.size() << ")"
	     << endl;
	nErrors++;
    }

    for (f1 = this->fields.begin(); f1 != this->fields.end();  ++f1)
    {
	f2 = that.fields.find(String(f1->second->name));
	if (f2 == that.fields.end())
	{
	    cout << "Both files do not have field '" << f1->second->name
		 << "'" << endl;
	    nErrors++;
	    continue;
	}

	// Make sure that we are comparing scalar data to scalar data, etc.
	if(f1->second->GetDims() != f2->second->GetDims())
	{
	    cout << "Field '" << f1->second->name
		 << "' has differing dimensionality ("
		 << f1->second->GetDims() << " and " << f2->second->GetDims()
		 << ")" << endl;
	    nErrors++;
	    continue;
	}

	cout << "Comparing '" << f1->second->name << "'" << endl;

	// Check that the same cells are marked for both fields
	vector<bool>::iterator cm1, cm2;
	int count = 0;	// Start at zero because the array starts at zero
			// (although element zero is unused)

	if (f1->second->cellmap.size() != f2->second->cellmap.size())
	{
	    cout << "Warning:  Field '" << f1->second->name
		 << "' has differing sizes ("
		 << f1->second->cellmap.size() << " and "
		 << f2->second->cellmap.size() << ")" << endl;
	    nErrors++;
	}
	cm1 = f1->second->cellmap.begin();
	cm2 = f2->second->cellmap.begin();
	for(; cm1 != f1->second->cellmap.end() &&
		cm2 != f2->second->cellmap.end(); ++cm1,++cm2)
	{
	    if (*cm1 != *cm2)
	    {
		if (nErrors < kMaxError)
		    cout << "Cell " << count
			 << " is marked as used by only "
			 << ((*cm1) ? "file1" : "file2") << "." << endl;
		nErrors++;
	    }
	    count++;

	}

	// Check that the data is the same
	vector<float>::iterator d1, d2;
	int inc = 1;
	if (f1->second->cellmap.size() != f2->second->cellmap.size())
	{
	    if (nErrors < kMaxError)
		cout << "Field '" << f1->second->name
		     << "' has differing amounts of data ("
		     << f1->second->cellmap.size() << " and "
		     << f2->second->cellmap.size() << ")" << endl;
	    nErrors++;
	    continue;
	}
	d1 = f1->second->data.begin();
	d2 = f2->second->data.begin();
	count = 0;
	if (f1->second->GetDims() == 2)
	    inc = 3;
	for(;  d1 != f1->second->data.end() && d2 != f2->second->data.end();
	    						       d1+=inc, d2+=inc)
	{
	    if (count++ == 0)	// Nothing store at zero;  fortran indexing
		continue;
	    if (f1->second->GetDims() == 1)
	    {
		float debug1 = *d1, debug2 = *d2;
//		if (*d1 != *d2)
		if (debug1 != debug2)
		{
		    if (nErrors < kMaxError)
			cout << "Field '" << f1->second->name << "' cell "
			     << count << " has differing values (" << *d1
			     << " and " << *d2 << ")" << endl;
		    nErrors++;
		}
	    }
	    else if (f1->second->GetDims() == 2)
	    {
		if (*d1 != *d2 || *(d1+1) != *(d2+1) || *(d1+2) != *(d2+2))
		{
		    if (nErrors < kMaxError)
			cout << "Field '" << f1->second->name << "' cell "
			     << count << " has differing values: "
			     << "(" << *d1 << "," << *(d1+1) << "," << *(d1+2)
			     << ")" << " and "
			     << "(" << *d2 << "," << *(d2+1) << "," << *(d2+2)
			     << ")" << endl;
		    nErrors++;
		}
	    }
	}

	// Compare boundary data
	PostField::BoundaryArray::iterator bd;
	if (f1->second->boundaryData.size() != f2->second->boundaryData.size())
	    cout << "Warning:  Field '" << f1->second->name
		 << "' has boundary data of "
		 << "differing sizes (" << f1->second->boundaryData.size()
		 << " and " << f2->second->boundaryData.size() << ")" << endl;
	map<int, Int> tmp;
	map<int, Int>::iterator tmpIter;
	for (bd = f1->second->boundaryData.begin();
	     			     bd != f1->second->boundaryData.end(); ++bd)
	    tmp[bd->second.id]++;
	for (bd = f2->second->boundaryData.begin();
	     			     bd != f2->second->boundaryData.end(); ++bd)
	    tmp[bd->second.id]++;

	for (tmpIter = tmp.begin();  tmpIter != tmp.end();  ++tmpIter)
	{
	    if (tmpIter->second == 2 &&
		(f1->second->boundaryData[tmpIter->first]!=
				     f2->second->boundaryData[tmpIter->first]))
	    {
		if (nErrors < kMaxError)
		    cout << "Field '" << f1->second->name << "' boundary "
			 << tmpIter->first << " have different values:" << endl
			 << "file1: "<< f1->second->boundaryData[tmpIter->first]
			 << "file1: "<< f1->second->boundaryData[tmpIter->first]
			 << endl;
		nErrors++;
	    }
	    else if (tmpIter->second == 1)
		{
		    BoundaryData b;
		    b.id = tmpIter->first;

		    if ((f1->second->boundaryData.find(tmpIter->first) !=
			 f1->second->boundaryData.end() &&
			 f1->second->boundaryData[tmpIter->first] != b) ||
			(f2->second->boundaryData.find(tmpIter->first) !=
			 f2->second->boundaryData.end() &&
			 f2->second->boundaryData[tmpIter->first] != b))
		    {
			if (nErrors < kMaxError)
			    cout << "Field '" << f1->second->name
				 << "' boundary "
				 << tmpIter->first << " have different values."
				 << endl;
			nErrors++;
		    }
		}
	}
    }

    if (nErrors > kMaxError)
	cout << "\t... " << nErrors - kMaxError << " more error"
	     << ((nErrors - kMaxError > 1) ? "s" : "") << " ..." << endl;
    return(nErrors == 0);
}

bool Post::WriteCCMIO( CCMIOID& processor )
{
    int i, dimSizes[] = { 0, 1, 3, 9 };
    unsigned int n;
    CCMIOID topology, cell, cellMap, fieldSet, field, data, bMap, bound, bfMap;
    CCMIOError error = kCCMIONoErr, *err = &error;

    // The map for the cell data is going to be exactly the same as for the
    // cells in the mesh, so find that map and reuse it.  This only works
    // because I happen to know that I'm writing the mesh and solution to
    // the same file.
    CCMIOReadProcessor(err, processor, NULL, &topology, NULL, NULL);
    CCMIOGetEntity(err, topology, kCCMIOCells, 0, &cell);
    CCMIOReadCells(err, cell, &cellMap, NULL, kCCMIOStart, kCCMIOEnd);

    // Get/make the first field set
    if (CCMIONextEntity(err, processor, kCCMIOFieldSet, &i, &fieldSet) != kCCMIONoErr)
    {
	*err = kCCMIONoErr;
	CCMIONewEntity(err, processor, kCCMIOFieldSet, kFieldSetName, &fieldSet);
    }

    FieldArray::iterator iter;
    for (iter = this->fields.begin();  iter != this->fields.end();  ++iter)
    {
	int nDims = iter->second->GetDims();

	// Write the cell data
	cout << "Writing post data field '" << iter->first.c_str() << "'"
	     << endl;
	CCMIONewField(err, fieldSet, iter->first.c_str(), iter->first.c_str(),
		      (CCMIODimensionality)nDims, &field);
	CCMIONewEntity(err, field, kCCMIOFieldData, NULL, &data);
	if (nDims == 1)
	    CCMIOWriteFieldDataf(err, data, cellMap, kCCMIOCell,
				 &iter->second->data[dimSizes[nDims]],
				 kCCMIOStart, kCCMIOEnd);
	else if (nDims == 2)
	{
	    int size = iter->second->data.size() / 3 - 1;
	    CCMIOID vectorField[3], fd;
	    char suffix[][5] = { " (U)", " (V)", " (W)" };

	    for (int i = 0;  i < 3;  ++i)
	    {
		vector<float> data(size);
		for (int j = 0;  j < size;  ++j)
		    data[j] = iter->second->data[3 + 3 * j + i];
		CCMIONewField(err, fieldSet, (iter->first + suffix[i]).c_str(),
			      (iter->first + suffix[i]).c_str(), kCCMIOScalar,
			      &vectorField[i]);
		CCMIONewEntity(err, vectorField[i], kCCMIOFieldData, NULL, &fd);
		CCMIOWriteFieldDataf(err, fd, cellMap, kCCMIOCell, &data[0],
				     kCCMIOStart, kCCMIOEnd);
		CCMIOWriteMultiDimensionalFieldData(err, field,
						    (CCMIOComponent)i,
						    vectorField[i]);
	    }
	}

	if (iter->second->boundaryData.empty())
	    continue;

	// Since the NPOST format doesn't save the boundary face data
	// individually, we need to get the list of faces for this piece of
	// data.
	int max = -1, nPts = (int)pow(3.0, static_cast < int > (nDims - 1));
	vector<int> boundaryIDs;
	vector<float> boundaryVals;
	PostField::BoundaryArray::iterator bIter;
	for (bIter = iter->second->boundaryData.begin();
	     bIter != iter->second->boundaryData.end();  ++bIter)
	{
	    CCMIOGetEntity(err, topology, kCCMIOBoundaryFaces, bIter->first,
			 &bound);
	    CCMIOReadFaces(err, bound, kCCMIOBoundaryFaces, &bfMap, NULL,
			   NULL, kCCMIOStart, kCCMIOEnd);
	    CCMIOEntitySize(err, bound, &n, NULL);
	    vector<int> ids(n);
	    CCMIOReadMap(err, bfMap, &ids[0], kCCMIOStart, kCCMIOEnd);
	    for (unsigned int j = 0;  j < n;  ++j)
	    {
		boundaryIDs.push_back(ids[j]);
		if (max < ids[j])
		    max = ids[j];
		for (int k = 0;  k < nPts;  ++k)
		    boundaryVals.push_back(bIter->second.value[k]);
	    }
	}

	// Write the boundary data
	CCMIONewEntity(err, field, kCCMIOFieldData, NULL, &data);
	CCMIONewEntity(err, data, kCCMIOMap, "Boundary data map", &bMap);
	CCMIOWriteMap(err, bMap, boundaryIDs.size(), max, &boundaryIDs[0],
		      kCCMIOStart, kCCMIOEnd);
	CCMIOWriteFieldDataf(err, data, bMap, kCCMIOFace,
			     &boundaryVals[0], kCCMIOStart, kCCMIOEnd);
    }
    CCMIOWriteProcessor(err, processor, NULL, NULL, NULL, NULL, NULL, NULL, 
		      NULL, &fieldSet);
    CCMIOCloseFile(err, topology);
    if (error != kCCMIONoErr)
	throw(IOException("Error writing solution data"));
    return(true);
}

#endif // POST_CPP
