/******************************************************************************\
*
*  MODULE:      PatchedMesh.cp
*
*  PURPOSE:     To provide a Demo ImportPIE template. 
*
*  FUNCTIONS:  Read a set of rectangular patches form a contour layer and create a quardilateral mesh
*
*  COMMENTS:   
*
\******************************************************************************/

#include <math.h>
#include <strstream.h>
#include <stdlib.h>
#include <string.h>
#include <stdio.h>

#include "ImportPIE.h"
#include "ANECB.h"

#pragma lib_export on
extern "C" void GetANEFunctions(long* numNames, struct ANEPIEDesc*** descriptors);
#pragma lib_export off

#define FALSE 0
#define TRUE 1

/*-------------------------------------------------------------------------
// sReserve
//--------------------------------------------------------------------------
// this function is used to handle dynamic array
//--------------------------------------------------------------------------*/
 
static void sReserve(char** buf, long* bufLen, long minLength)
{
	long newLen;

	if (*bufLen > minLength)
		return;

	newLen = *bufLen*1.5;
	if (newLen <= minLength)
		newLen = minLength + 32;

	if (*buf)
		*buf = (char*)realloc(*buf, newLen);
	else
		*buf = (char*)malloc(newLen);
	*bufLen = newLen;

} // sReserve

static char* sOutputBuffer;
static long sOutputBufferMaxSize;
static long sOutputBufferLen;

static void sInitOutputBuffer()
{
	sOutputBuffer = 0;
	sOutputBufferMaxSize = 0;
	sOutputBufferLen = 0;

	sReserve(&sOutputBuffer, &sOutputBufferMaxSize, 100);
	sOutputBufferLen = 0;
	sOutputBuffer[0] = '\0'; 
}

static void sWriteString(char* str)
{
	sReserve(&sOutputBuffer, &sOutputBufferMaxSize, sOutputBufferLen+strlen(str)+1);
	strcat(sOutputBuffer+sOutputBufferLen, str);
	sOutputBufferLen += strlen(str);
}

static void sWriteInt(long number)
{
	static char numberBuf[100];
	sprintf(numberBuf, "%ld", number);
	sWriteString(numberBuf);
}

static void sWriteDouble(double number)
{
	static char numberBuf[100];
	sprintf(numberBuf, "%lg", number);
	sWriteString(numberBuf);
}

static char* sGetOutputBuffer()
{
	return sOutputBuffer;
}

static void sTerminateOutputBuffer()
{
	free(sOutputBuffer);
}

//==============================================================================
// ImportPIE definitions
//==============================================================================

//
// some struct declerations that we will use

struct SPoint							// 2D point struct
	{
		double x;
		double y;

	}; // SPoint

struct SPatch						// a quadrilateral patch struct
	{
		SPoint point[4];			// 4 points of the patch
		long nodeNumber[4];		// 4 node numbers of the 4 points
		long neighbour[4];			// patch number of neighbour pathces
		int backindex[4];			// number of neighbour pathces
		double density[2];			// x, y densitys
		int densWasCorrected[2];		// true after correcting the density

	}; // SPatch

typedef SPatch* SPatchPtr;

struct SElement					// a quadrilateral element
	{
		long nodeNumber[4];		// 4 node numbers of the 4 points

	}; // SElement


	//
	// definition for patches array

	SPatchPtr*	patches = 0;
	long patchesArraySize = 0;
	long numPatches = 0;

	SPoint* nodes = 0;
	long nodesArraySize = 0;
	long currentNodeNumber = 0;


//------------------------------------------------------------------------------
// MIN/MAX:
//------------------------------------------------------------------------------
// Utility macro to get the min/max value of 2 numbers
//------------------------------------------------------------------------------

#define MIN(a,b) ((a)<(b)?(a):(b))
#define MAX(a,b) ((a)>(b)?(a):(b))


//------------------------------------------------------------------------------
// PointsEqual:
//------------------------------------------------------------------------------
// Utility function to compare two points for equality.
// Usually the two points are not perfectly equal up to the last significant digit,
// so compare up to 6 digits
//------------------------------------------------------------------------------

int PointsEqual(SPoint* pt1, SPoint* pt2)
{
	return fabs( (pt1->x - pt2->x) / pt1->x ) < 1E-6 &&  fabs( (pt1->y - pt2->y) / pt1->y ) < 1E-6;

} // PointsEqual


//------------------------------------------------------------------------------
// GetDensityIndex:
//------------------------------------------------------------------------------
// figure out if the 'neighbour' is on the x direction (return 0)
// or on the y direction (return 1)
//------------------------------------------------------------------------------

int GetDensityIndex(SPatch* patch, int neighbour)
{
	double maxdy = 0;
	long index;

	for (int i=0; i<4; i++)
	{
		double y1 = patch->point[i].y;
		double y2 = patch->point[(i+1)%4].y;

		if (y1 > y2)																		// check that this edge go down
		{
			if ( y1-y2 > maxdy)														// look for the leftmost edge
			{																					// this edge is in the y direction!
				index = i;																	// save it's index
				maxdy = y1-y2;
			}
		}
	} // for i

	return neighbour%2==index%2?1:0;										// if the nieghbour modulo 2 is the same as the found
																								// edge index modulo 2, it's the y direction (return 1)

} // GetDensityIndex


//------------------------------------------------------------------------------
// GetPatchNeighbourIndex:
//------------------------------------------------------------------------------
// Utility function to find the index of 'neighbour' patch within
// the neighbours of 'patch'
//------------------------------------------------------------------------------

int GetPatchNeighbourIndex(SPatch* patch, SPatch* neighbour)
{
	for (int i=0; i<4; i++)
	{
		if (patches[patch->neighbour[i]]==neighbour)
			return i;
	}

	return -1;

} // GetPatchNeighbourIndex


//------------------------------------------------------------------------------
// GetMaxPatchDensity:
//------------------------------------------------------------------------------
//
//------------------------------------------------------------------------------

long GetMaxPatchDensity(long patchIndex, long neighbourIndex)
{
	if (patchIndex==-1)
		return 0;

	SPatch* patch = patches[patchIndex];

	int newIndex = (neighbourIndex+2)%4;
	long d = GetMaxPatchDensity(patch->neighbour[newIndex], patch->backindex[newIndex]);
	int densIndex = GetDensityIndex(patch, neighbourIndex);
	
	return MAX(d, patch->density[densIndex]);

} // GetMaxPatchDensity


//------------------------------------------------------------------------------
// SetPatchDensity:
//------------------------------------------------------------------------------
//
//------------------------------------------------------------------------------

void SetPatchDensity(long patchIndex, long neighbourIndex, double dens)
{
	if (patchIndex==-1)
		return;

	SPatch* patch = patches[patchIndex];

	int newIndex = (neighbourIndex+2)%4;
	SetPatchDensity(patch->neighbour[newIndex], patch->backindex[newIndex], dens);

	int densIndex = GetDensityIndex(patch, neighbourIndex);
	
	//
	// should check for plonters here!@#$%:
	// if (patch->densWasCorrected[densIndex]) { do something, quick! }

	if (patch->densWasCorrected[densIndex])
	{
		patch->densWasCorrected[densIndex] = FALSE;
	}

	patch->density[densIndex] = dens;
	patch->densWasCorrected[densIndex] = TRUE;

} // SetPatchDensity


//------------------------------------------------------------------------------
// SearchNode:
//------------------------------------------------------------------------------
//
//------------------------------------------------------------------------------

long SearchNode(SPoint pt, double epsilon)
{
	for (long i=0; i<currentNodeNumber; i++)
	{
		if ( fabs(pt.x-nodes[i].x)<=epsilon && fabs(pt.y-nodes[i].y)<=epsilon )
			return i;
	}

	return -1;

} // SearchNode


//------------------------------------------------------------------------------
// GPatchedMeshImportProc:
//------------------------------------------------------------------------------
// Import and mesh:
// 	Import a set of (hopefuly) rectangular patches from the source information layer
//	and build a quad-mesh
//------------------------------------------------------------------------------

static void
GPatchedMeshImportProc(ANE_PTR aneHandle, const ANE_STR fileName, ANE_PTR layerHandle)
{
	//
	// prepare array of patches

	patches = 0;
	patchesArraySize = 0;
	numPatches = 0;
	sReserve((char**)&patches, &patchesArraySize, sizeof(SPatch*)*10);			// reserve room for at least 10 patches

	char* ibuff;
	char* obuff;

	//
	// get a string from the source layer describing the pathces.
	// the string is in the information layer export form:
	//
	// <num-of-points> <parameter1> <parameter2>
	// <x1> <y1> 
	// <x2> <y2> 
	// <x3> <y3> 
	// <x4> <y4> 
	// <x5> <y5> 
	//
	// Notes:
	//	<num-of-points> should be 5
	//	<parameter1> will be used as the density in the x direction
	//	<parameter2> will be used as the density in the y direction
	//	<x5> and <y5> should be the same as <x1> and <y1>!

	ANE_ExportTextFromOtherLayer(aneHandle, layerHandle, &ibuff);

	//
	// prepare to read the patches

	istrstream importText(ibuff);
	long numpoints;
	double x, y;
	int problemEncountered = FALSE;
	SPatch* patch;

	//
	// loop and read all the patches

	while(TRUE)
	{
		//
		// skip comments and empty lines

		do
		{
			importText.ignore(256, '\n');
		} while (importText.peek() == '#' || importText.peek() == '\n');

		if (!importText || importText.eof())
			break;

		//
		// read header line

		importText >> numpoints;					// should be 5 points
		if (numpoints != 5 || importText.eof())
		{
			problemEncountered = TRUE;
			break;
		}

		patch = new SPatch;
		importText >> patch->density[0];		if (importText.eof()) break;
		importText >> patch->density[1];

		patch->densWasCorrected[0] = FALSE;
		patch->densWasCorrected[1] = FALSE;

		//
		// skip comment and empty lines

		do
		{
			importText.ignore(256, '\n');
		} while (importText.peek() == '#' || importText.peek() == '\n');

		if (!importText || importText.eof())
			break;

		//
		// read x and y values of the 5 points

		importText >> patch->point[0].x;
		importText >> patch->point[0].y;
		patch->nodeNumber[0] = -1;				// reset node number
		patch->neighbour[0] = -1;					// reset neighbour number

		importText >> patch->point[1].x;
		importText >> patch->point[1].y;
		patch->nodeNumber[1] = -1;				// reset node number
		patch->neighbour[1] = -1;					// reset neighbour number

		importText >> patch->point[2].x;
		importText >> patch->point[2].y;
		patch->nodeNumber[2] = -1;				// reset node number
		patch->neighbour[2] = -1;					// reset neighbour number

		importText >> patch->point[3].x;
		importText >> patch->point[3].y;
		patch->nodeNumber[3] = -1;				// reset node number
		patch->neighbour[3] = -1;					// reset neighbour number

		importText >> x;								// should be equal to patch->point[0].x
		importText >> y;								// should be equal to patch->point[0].y

		if (!importText)								// check for a problem
			break;

		//
		// add the new patch to the list

		sReserve((char**)&patches, &patchesArraySize, sizeof(SPatch*)*(numPatches+1));		// reserve room for at least 'numPatches' patches
		patches[numPatches++] = patch;																						// insert patch into array
	};

//???	ANE_MemDelete(aneHandle, ibuff);							// clear the block of memory returned from 'ANE_ExportTextFromOtherLayer()'

	//
	// prepare an array of node positions
	// max possible nodes are 4 times the patches count

	nodes = 0;
	nodesArraySize = 0;
	currentNodeNumber = 0;
	sReserve((char**)&nodes, &nodesArraySize, sizeof(SPoint)*(numPatches*4+1));										// reserve room for nodes

	//
	// prepare an array of elements

	SElement* elements = 0;
	long numElements = 0;
	long elemsArraySize = 0;
	sReserve((char**)&elements, &elemsArraySize, sizeof(SElement)*(numPatches*4+1));				// reserve room for elements

	//
	// loop for all pathces and calc the node numbers:

	long i, j, k, l;

	for (i=0; i<numPatches; i++)																									// for all patches
	{
		patch = patches[i];																											// save a pointer to the current patch
		
		//
		// loop for all nodes and assign them the proper node-numbers

		for (j=0; j<4; j++)																											// loop for patch's 9 nodes
		{
			//
			// search all known nodes, is it equal to our vertex?

			int nodeFound = FALSE;
			for (k=0; k<currentNodeNumber && !nodeFound; k++)													// check all known nodes
				if (PointsEqual(&nodes[k], &patch->point[j]))															// test equality
					nodeFound = TRUE;																								// a node was found
			
			if (nodeFound)
			{
				// a node was found
				// take it's number
				
				patch->nodeNumber[j] = k-1;																					// take the found node number
			}
			else
			{
				// a node was NOT found
				// create a new node and
				// take a new number
				
				patch->nodeNumber[j] = currentNodeNumber;															// take the found node number
				nodes[currentNodeNumber++] = patch->point[j];														// save the node's position
			}

		} // for j

	} // for i

	//
	// figure out patch neighbours

	for (i=0; i<numPatches; i++)																									// for all patches
	{
		patch = patches[i];																											// save a pointer to the current patch
		
		//
		// loop for all nodes and assign them the proper node-numbers

		for (j=0; j<4; j++)																											// loop for patch's 4 nodes
			if (patch->neighbour[j]==-1)
			{
				long n1 = patch->nodeNumber[j];
				long n2 = patch->nodeNumber[(j+1)%4];

				// search neighbour patch
				
				int found = FALSE;
				for (k=0; k<numPatches && !found; k++)
				{
					SPatch* p =  patches[k];
					
					for (l=0; l<4 && !found; l++)
						if (p->nodeNumber[l]==n2 && p->nodeNumber[(l+1)%4]==n1 )
						{
							found = TRUE;
							patch->neighbour[j]=k;
							patch->backindex[j]=l;

							p->neighbour[l]=i;
							p->backindex[l]=j;
						} // for l, if neighbour found

				} // for k

			} // for j

	} // for i

	//
	// prepare meshing densities

	for (i=0; i<numPatches; i++)																									// for all patches
	{
		long dens1, dens2, dens, densIndex;

		patch = patches[i];																											// save a pointer to the current patch

		densIndex = GetDensityIndex(patch, 0);

		if (!patch->densWasCorrected[densIndex])
		{
			dens1 = GetMaxPatchDensity(patch->neighbour[0], patch->backindex[0] );
			dens2 = GetMaxPatchDensity(patch->neighbour[2], patch->backindex[2] );

			dens  = MAX(patch->density[densIndex], MAX(dens1, dens2));

			patch->density[densIndex] = dens;
			patch->densWasCorrected[densIndex] = TRUE;

			SetPatchDensity(patch->neighbour[0], patch->backindex[0], dens);
			SetPatchDensity(patch->neighbour[2], patch->backindex[2], dens);
		
		}

		densIndex = (densIndex+1)%2;

		if (!patch->densWasCorrected[densIndex])
		{
			dens1 = GetMaxPatchDensity(patch->neighbour[1], patch->backindex[1] );
			dens2 = GetMaxPatchDensity(patch->neighbour[3], patch->backindex[3] );

			dens  = MAX(patch->density[densIndex], MAX(dens1, dens2));

			patch->density[densIndex] = dens;
			patch->densWasCorrected[densIndex] = TRUE;

			SetPatchDensity(patch->neighbour[1], patch->backindex[1], dens);
			SetPatchDensity(patch->neighbour[3], patch->backindex[3], dens);
		}

	} // for i


	//**************
	// do the meshing here
	//**************


	for (i=0; i<numPatches; i++)																									// for all patches
	{
		patch = patches[i];																											// save a pointer to the current patch

		double maxx = MAX(patch->point[0].x, MAX(patch->point[1].x, MAX(patch->point[2].x, patch->point[3].x)));
		double minx = MIN(patch->point[0].x, MIN(patch->point[1].x, MIN(patch->point[2].x, patch->point[3].x)));
		double maxy = MAX(patch->point[0].y, MAX(patch->point[1].y, MAX(patch->point[2].y, patch->point[3].y)));
		double miny = MIN(patch->point[0].y, MIN(patch->point[1].y, MIN(patch->point[2].y, patch->point[3].y)));
		double epsilon = 1e-6*((maxx-minx)+(maxy-miny));

		long dens1, dens2, densIndex;
		densIndex = GetDensityIndex(patch, 0);
		dens1 = patch->density[densIndex];
		dens2 = patch->density[(densIndex+1)%2];

		long numNodes = (dens1+2)*(dens2+2);
		long* nodeNumber = 0;
		long nodeNumberArraySize = 0;
		sReserve((char**)&nodeNumber, &nodeNumberArraySize, sizeof(long)*(numNodes+1));		// reserve room for node numbers

		for (j=0; j<numNodes; j++)
			nodeNumber[j] = -1;

		//
		// save patch's 4 nodes

		nodeNumber[0]										= patch->nodeNumber[0];
		nodeNumber[dens1+1]							= patch->nodeNumber[1];
		nodeNumber[(dens1+2)*(dens2+2)-1]	= patch->nodeNumber[2];
		nodeNumber[(dens1+2)*(dens2+1)]		= patch->nodeNumber[3];

		for (j=0; j<dens2+2; j++)
			for (k=0; k<dens1+2; k++)
			{
				long n = j*(dens1+2)+k;
				if (nodeNumber[n] == -1)
				{
					sReserve((char**)&nodes, &nodesArraySize, sizeof(SPoint)*(currentNodeNumber+2));										// reserve room for nodes

					SPoint p1, p2, pt;
					double ratio;
					
					// pt1 is between node[0] and node[3]
					
					ratio = (double)(j)/(dens2+1);
					p1.x = patch->point[0].x*(1-ratio) + patch->point[3].x*ratio;
					p1.y = patch->point[0].y*(1-ratio) + patch->point[3].y*ratio;
	
					// pt2 is between node[1] and node[2]
					
					p2.x = patch->point[1].x*(1-ratio) + patch->point[2].x*ratio;
					p2.y = patch->point[1].y*(1-ratio) + patch->point[2].y*ratio;
					
					// our point should be between p1 and p2

					ratio = (double)(k)/(dens1+1);
					pt.x = p1.x*(1-ratio) + p2.x*ratio;
					pt.y = p1.y*(1-ratio) + p2.y*ratio;
					
					long foundNode = SearchNode(pt, epsilon);
					if (foundNode == -1)
					{
						nodeNumber[n] = currentNodeNumber;
						nodes[currentNodeNumber++] = pt;
					}
					else
					{
						nodeNumber[n] = foundNode;
					}
				}
			}

		//
		// reserve space for new elements

		long newElements = (dens1+1) * (dens2+1);
		sReserve((char**)&elements, &elemsArraySize, sizeof(SElement)*(numElements+newElements+1));		// reserve room for elements

		//
		// fill new elements with the proper node numbers
		
		for (j=0; j<dens2+1; j++)
			for (k=0; k<dens1+1; k++)
			{
				SElement element;

				element.nodeNumber[0] = nodeNumber[j*(dens1+2)+k];
				element.nodeNumber[1] = nodeNumber[j*(dens1+2)+k+1];
				element.nodeNumber[2] = nodeNumber[(j+1)*(dens1+2)+k+1];
				element.nodeNumber[3] = nodeNumber[(j+1)*(dens1+2)+k];
	
				elements[numElements++] = element;																											// insert patch into array
			}

		delete nodeNumber;								// clear the block of memory used for the node-numbers

	} // for i

	delete patches;										// clear the block of memory used for the patches

	//
	// prepare output string
	// the string should be in the mesh layer's import form:
	//
	// <num-of-elements> <num-of-nodes> <num-of-element-parameters> <num-of-node-parameters>
	// N <node-number1> <x1> <y1> 
	// N <node-number2> <x2> <y2> 
	// .
	// .
	// E <elem-number1> <node1> <node2> <node3> <node4>
	// E <elem-number2> <node1> <node2> <node3> <node4>
	// .
	// .
	//
	// Notes:
	//	<num-of-element-parameters> and <num-of-node-parameters>
	//	 are both 0 in our example

	sInitOutputBuffer();
//	ostrstream exportText;
	
	// export header line

	sWriteInt(numElements);
	sWriteString(" ");
	sWriteInt(currentNodeNumber);
	sWriteString(" 0 0\n");
//	exportText << numElements << " " << currentNodeNumber << " 0 0\n";

	// export node lines

	for (i=0; i<currentNodeNumber; i++)
	{
		sWriteString("N ");
		sWriteInt(i+1);
		sWriteString(" ");
		sWriteDouble(nodes[i].x);
		sWriteString(" ");
		sWriteDouble(nodes[i].y);
		sWriteString("\n");
		// exportText << "N " << i+1 << " " << nodes[i].x << " " << nodes[i].y << "\n";
	}

	// export element lines

	for (i=0; i<numElements; i++)
	{
		sWriteString("E ");
		sWriteInt(i+1);
		sWriteString(" ");
		sWriteInt(elements[i].nodeNumber[0]+1);
		sWriteString(" ");
		sWriteInt(elements[i].nodeNumber[1]+1);
		sWriteString(" ");
		sWriteInt(elements[i].nodeNumber[2]+1);
		sWriteString(" ");
		sWriteInt(elements[i].nodeNumber[3]+1);
		sWriteString("\n");
//		exportText << "E " << i+1 << " " << elements[i].nodeNumber[0]+1
//												 << " " << elements[i].nodeNumber[1]+1
//												 << " " << elements[i].nodeNumber[2]+1
//												 << " " << elements[i].nodeNumber[3]+1 << "\n";
	}

	obuff = sGetOutputBuffer();

	ANE_ImportTextToLayer(aneHandle, obuff);

	sTerminateOutputBuffer();
	//
	// cleanup
	
	delete nodes;										// clear the block of memory used for the nodes
	delete elements;								// clear the block of memory used for the elements

} // GPatchedMeshImportProc



//=================================================================
//	PIE handling functions and declerations
//=================================================================

static struct ANEPIEDesc* gPIEDescroptors[1];									// list of PIE descriptors for all parts

void
GetANEFunctions(long* numNames, struct ANEPIEDesc*** descriptors)
{ 
	*numNames = 0;																			// initialize the number of PIEs

	//
	// prepare Import descriptor for "GPatchedMeshImportProc" ImportPIE

	ANEPIEDesc*		meshDesc				= new ANEPIEDesc;				// create PIE descriptor
	ImportPIEDesc*	meshImportDesc	= new ImportPIEDesc;			// create ImportPIE additional information record

	//
	// fill in the PIE descriptor

	meshDesc->name			= "MeshPatches";										// PIE name
	meshDesc->type			= kImportPIE;											// type of PIE: ImportPIE
	meshDesc->descriptor	= meshImportDesc;									// address of additional information record

	//
	// fill in the ImportPIE additional information record

	meshImportDesc->version					= IMPORT_PIE_VERSION;		// descriptor version, use constant here...
	meshImportDesc->name						= "Mesh Patches...";				// title of menu
	meshImportDesc->importFlags			= kImportFromLayer;			// import from another layer
	meshImportDesc->toLayerTypes		= kPIEQuadMeshLayer;			// import only TO quad-mesh layers
	meshImportDesc->fromLayerTypes	= (EPIELayerType)(kPIEInformationLayer |
																kPIEDomainLayer);				// import FROM information (contour) and domain layers
	meshImportDesc->doImportProc		= GPatchedMeshImportProc;	// the address of the import function

	gPIEDescroptors[(*numNames)++] = meshDesc;							// put descriptor into array

	*descriptors = gPIEDescroptors;

} // GetANEFunctions