// Wsimplex.cpp: implementation of the Wsimplex class.
//
//////////////////////////////////////////////////////////////////////

#include "dxstdafx.h"
#include "Wsimplex.h"
#include "Vector.hpp"
#include "Wbucket.h"
#include "WsubMesh.h"

//////////////////////////////////////////////////////////////////////
// ruction/Destruction
//////////////////////////////////////////////////////////////////////

Wsimplex::Wsimplex(int dim)
{
	deleted = 0;
	circumcentre = new Vector(0,0,0);
	dimension = dim;

	vertex = new Vector*[dimension+1];
	neighbour = new Wsimplex*[dimension+1];
}

Wsimplex::Wsimplex(Vector *q, Vector *w, int order)
{	
	dimension = 1;
	deleted = 0;
	vertex = new Vector*[dimension+1];
	neighbour = new Wsimplex*[dimension+1];

	circumcentre = new Vector(0,0,0);

	neighbour[0] = NULL;
	neighbour[1] = NULL;
	inFrontOf = NULL;
	behind = NULL;

	if(!order)
	{
	vertex[0] = q;
	vertex[1] = w;
	}
	else
	{
		if(q<w)
		{
			vertex[0]=q;
			vertex[1]=w;
		}
		else
		{
			vertex[0] = w;
			vertex[1] = q;
		}

	}
//WARNING correct ordering by pointer value assumed if(order == 0)
	updateCentre();
}


Wsimplex::Wsimplex(Vector *q, Vector *w, Vector *e, int order)
{	
	dimension = 2;
	deleted = 0;
	vertex = new Vector*[dimension+1];
	neighbour = new Wsimplex*[dimension+1];

	circumcentre = new Vector(0,0,0);
	neighbour[0] = NULL;
	neighbour[1] = NULL;
	neighbour[2] = NULL;
	inFrontOf = NULL;
	behind = NULL;

	if(!order)
	{
	vertex[0] = q;
	vertex[1] = w;
	vertex[2] = e;
	}
	else
	{
		if(q<w)
		{
			if(q<e)
			{
				vertex[0]=q;
				if(w<e)
				{
					vertex[1]=w;
					vertex[2]=e;
				}
				else
				{
					vertex[1]=e;
					vertex[2]=w;
				}
			}
			else
			{
				vertex[0] = e;
				vertex[1] = q;
				vertex[2] = w;
			}
		}
		else
		{
			if(q<e)
			{
				vertex[0] = w;
				vertex[1] = q;
				vertex[2] = e;
			}
			else
			{
				vertex[2] = q;
				if(w<e)
				{
					vertex[0]=w;
					vertex[1]=e;
				}
				else
				{
					vertex[0]=e;
					vertex[1]=w;
				}
			}
		}
	}

//WARNING correct ordering by pointer value assumed if(order == 0)
	updateCentre();
}

Wsimplex::Wsimplex(Vector *sharp, Wsimplex *flat)
{
	dimension = flat->getDimension()+1;
	if(dimension == 2)
		printf("brek \n");

	vertex = new Vector*[dimension+1];
	neighbour = new Wsimplex*[dimension+1];
	
	deleted = 0;
	circumcentre = new Vector(0,0,0);
	
	for(int shg = 0;shg < dimension+1;shg++)
		neighbour[shg] = NULL;
	inFrontOf = NULL;
	behind = NULL;

 	Vector** tri = flat->getVertices();

	int sindex;
	for(sindex = 0;sindex<dimension;sindex++)
		if(sharp < tri[sindex])	break;;

	vertex[sindex] = sharp;
	neighbour[sindex] = flat->getBehind();

	switch(dimension)
	{
	case 3: vertex[(sindex>2)?2:3] = tri[2];
			neighbour[(sindex>2)?2:3] = flat->getNeighboursFront(2);
	case 2: vertex[(sindex>1)?1:2] = tri[1];
			neighbour[(sindex>1)?1:2] = flat->getNeighboursFront(1);
	case 1: vertex[sindex?0:1] = tri[0];
			neighbour[sindex?0:1] = flat->getNeighboursFront(0);
	}
	
//tell the neighbours to link here
	for(int i=0;i<dimension+1;i++)
		if(neighbour[i])
			neighbour[i]->helloNeighbour(this);

 	updateCentre();
}


Wsimplex::~Wsimplex()
{
	if(!backlist.empty())
	{
		::std::vector<Wbucket*>::const_iterator a = backlist.begin();
		while(a != backlist.end())
		{
			(*a)->eliminate(this);
			a++;
		}
	}
//WARNING cc is possibly used in refine
	delete circumcentre;
}


Vector* Wsimplex::getCircumcentre() 
{
	return circumcentre;
}

double Wsimplex::getRadius() 
{
	return radius;
}

int Wsimplex::inCircle(Vector *guess) 
{
	Vector dist = *circumcentre - *guess;
	return sqrad - 0.0000001> dist.norm2();
}

void Wsimplex::linkback(Wbucket *container)
{
	backlist.push_back(container);
}


void Wsimplex::updateCentre()
{
	if(dimension == 3)
	{
	double bax = vertex[1]->x - vertex[0]->x;
	double cax = vertex[2]->x - vertex[0]->x;
	double dax = vertex[3]->x - vertex[0]->x;
	double bay = vertex[1]->y - vertex[0]->y;
	double cay = vertex[2]->y - vertex[0]->y;
	double day = vertex[3]->y - vertex[0]->y;
	double baz = vertex[1]->z - vertex[0]->z;
	double caz = vertex[2]->z - vertex[0]->z;
	double daz = vertex[3]->z - vertex[0]->z;

	double bxcx = bay * caz - baz * cay;
	double bxcy = baz * cax - bax * caz;
	double bxcz = bax * cay - bay * cax;

	double dxbx = day * baz - daz * bay;
	double dxby = daz * bax - dax * baz;
	double dxbz = dax * bay - day * bax;

	double cxdx = cay * daz - caz * day;
	double cxdy = caz * dax - cax * daz;
	double cxdz = cax * day - cay * dax;

	double da = dax*dax + day*day + daz*daz;
	double ca = cax*cax + cay*cay + caz*caz;
	double ba = bax*bax + bay*bay + baz*baz;

	double nx = da * bxcx + ca * dxbx + ba * cxdx;
	double ny = da * bxcy + ca * dxby + ba * cxdy;
	double nz = da * bxcz + ca * dxbz + ba * cxdz;

	double det = 2 * ( bax * cay * daz + bay * caz * dax + baz * cax * day
		    - bax * caz * day - bay * cax * daz - baz * cay * dax);

	sqrad = (nx*nx + ny*ny + nz*nz) / (det * det);
	radius = sqrt(sqrad);
	circumcentre->x = vertex[0]->x + nx / det;
	circumcentre->y = vertex[0]->y + ny / det;
	circumcentre->z = vertex[0]->z + nz / det;
	
	}
	else
	if(dimension == 2)
	{
/*
	double els[9];
	els[0] = vertex[1]->x - vertex[0]->x;
	els[3] = vertex[1]->y - vertex[0]->y;
	els[6] = vertex[1]->z - vertex[0]->z;
	els[1] = vertex[2]->x - vertex[0]->x;
	els[4] = vertex[2]->y - vertex[0]->y;
	els[7] = vertex[2]->z - vertex[0]->z;
	els[2] = els[3]*els[7] - els[6]*els[4];
	els[5] = els[6]*els[1] - els[0]*els[7];
	els[8] = els[0]*els[4] - els[3]*els[1];
	Matrix M(3,3,els);

	double vavg[3];
	vavg[0] = (	els[0]*(vertex[1]->x + vertex[0]->x) +
				els[3]*(vertex[1]->y + vertex[0]->y) +
				els[6]*(vertex[1]->z + vertex[0]->z)
				)/2;
	vavg[1] = (	els[1]*(vertex[2]->x + vertex[0]->x) +
				els[4]*(vertex[2]->y + vertex[0]->y) +
				els[7]*(vertex[2]->z + vertex[0]->z)
				)/2;
	vavg[2] =	els[2]*vertex[0]->x +
				els[5]*vertex[0]->y +
				els[8]*vertex[0]->z;
	Vector V(3,vavg);

	M.invert();
	Vector ResX = M*V;
	
	circumcentre->x = ResX(1);
	circumcentre->y = ResX(2);
	circumcentre->z = ResX(3);
	double r1 = circumcentre->x - vertex[0]->x;
	double r2 = circumcentre->y - vertex[0]->y;
	double r3 = circumcentre->z - vertex[0]->z;
	sqrad = r1*r1 + r2*r2 + r3*r3;
	radius = sqrt(sqrad);
*/	
	}
	//WARN 'dim=1' case missing
}


WsubMesh* Wsimplex::cavity(Vector *proxy)
{
	if(deleted)
		return (WsubMesh*)0x00000001;
	if(inCircle(proxy))
	{
		deleted = 1;

		//WARN out of chain here, or undertaker later

		//into deleted chain
		WsubMesh* collector = new WsubMesh();
		
		collector->setInsider(this);

		for(int i=0;i<dimension+1;i++)
		{
			WsubMesh* a = NULL;
			if(neighbour[i])
			{
				a = neighbour[i]->cavity(proxy);
			}
//else it is outer space, notice a == NULL then
			if(a)
			{
//WARNING n[i] returned, therefore all his neighbours are dead,
// he will never be called again, but cannot be deleted yet!!!
// 0x1 return indicates a tetra earlier in call chain
				if(a != (WsubMesh*)0x00000001)
				{
//					delete neighbour[i];
					collector->merge(a);
				}
			}
			else
			{
				Wsimplex* n;
				if(dimension == 3)
					n = new Wsimplex(
					vertex[(i==0)?1:0],
					vertex[(i<2)?2:1],
					vertex[(i<3)?3:2]);
				else
					n= new Wsimplex(
					vertex[(i==0)?1:0],
					vertex[(i==2)?1:2]);
				n->behind = neighbour[i];
				collector->add(n);
			}

		}

		return collector;
		
	}
	else
		return NULL;
}

void Wsimplex::helloNeighbour(Wsimplex* theguynext)
{
	Vector** ott = theguynext->getVertices();
	if(!(vertex[0] == ott[0] || vertex[0] == ott[1]))
	{
		neighbour[0] = theguynext;
		return;
	}
	
	if(!(vertex[dimension] == ott[dimension] || vertex[dimension] == ott[dimension-1]))
	{
		neighbour[dimension] = theguynext;
		return;
	}

	for(int t=1;t<dimension;t++)
		if(!(vertex[t] == ott[t] || vertex[t] == ott[t+1] || vertex[t] == ott[t-1]))
		{
			neighbour[t] = theguynext;
			return;
		}
	printf("ERROR: helloNeighbour found no matching nodes!\n");
}

int Wsimplex::getDimension()
{
	return dimension;
}

Vector** Wsimplex::getVertices()
{
	return vertex;
}

Wsimplex** Wsimplex::getNeighbours()
{
	return neighbour;
}


Wsimplex::Wsimplex(Vector *q, Vector *w, Vector *e, Vector *r)
{

//WARNING position is not checked for qualit
//quality NOT SET!!
//this constructor was marked obsolete
//but it is used for constructing from file

	dimension = 3;
	deleted = 0;
	for(int i=0;i<4;i++)
		neighbour[i] = NULL;
	vertex[0] = q;
	vertex[1] = w;
	vertex[2] = e;
	vertex[3] = r;
	circumcentre = new Vector(0,0,0);
	updateCentre(); 
}


void Wsimplex::print()
{
//	cout << "Wsimplex:\n";
//	for(int i=0;i<=dimension;i++)
//		vertex[i]->print();
}

int Wsimplex::adjacentToNothing()
{
	int c=0;
	for(int y=0;y<=dimension;y++)
		if(!neighbour[y])
			c++;
	return c;
}

int Wsimplex::verifyDelaunay(Vector* proxy)
{
	if(proxy)
	{
		if(inSmallSphere(proxy))
			return 0;
		else
			return 1;
	}
	else
	return 1;

}

int Wsimplex::inSmallSphere(Vector *guess)
{
	Vector xx = (*circumcentre) - (*guess);

	return sqrad-0.1 > xx.norm2();
}

void Wsimplex::cutlinkback(Wbucket *container)
{
	if(!backlist.empty())
	{
		::std::vector<Wbucket*>::iterator a = backlist.begin();
		while(a != backlist.end())
		{
			if(*a == container)
			{
				backlist.erase(a);
				break;
			}
			else
				a++;
		}
	}
}

Vector* Wsimplex::getInteriorPoint()
{
	double ix = 0, iy = 0, iz = 0;
	for(int c=0;c<=dimension;c++)
	{
		ix += vertex[c]->x;
		iy += vertex[c]->y;
		iz += vertex[c]->z;
	}
	return new Vector(ix/(double)(dimension+1),iy/(double)(dimension+1),iz/(double)(dimension+1));
}

int Wsimplex::nearEdge()
{
//WARNING no id in Vector, should identify boundary nodes using geometry
	return 0;

}

Wsimplex* Wsimplex::getNeighboursFront(int which)
{
	if(neighbour[which])
		return neighbour[which]->getFront();
	else
		return 0x0;
}

Wsimplex* Wsimplex::getFront()
{
	return inFrontOf;
}

Wsimplex* Wsimplex::getBehind()
{
	return behind;
}

int Wsimplex::isComplete()
{
	Wsimplex** a = neighbour;
	int com = 0;
	for(int c=0; c<dimension; c++)
		com = com && *(a++);
	return com;
}

int Wsimplex::join(Wsimplex *other)
{
	Vector** ott = other->getVertices();

	if(dimension == 1)
	{
		if(vertex[0] == ott[0])
		{
			neighbour[1] = other;
			other->neighbour[1] = this;
		}
		else
			if(vertex[0] == ott[1])
			{
				neighbour[1] = other;
				other->neighbour[0] = this;
			}
		else
			if(vertex[1] == ott[0])
			{
				neighbour[0] = other;
				other->neighbour[1] = this;
			}
		else
			if(vertex[1] == ott[1])
			{
				neighbour[0] = other;
				other->neighbour[0] = this;
			}
		return ((neighbour[0] && neighbour[1])?1:0) |
			((other->neighbour[0] && other->neighbour[1])?2:0);
	}

	if(vertex[0] == ott[0])
	{
		if(vertex[1] == ott[1])
		{
			neighbour[2] = other;
			other->neighbour[2] = this;
		}
		else
			if(vertex[1] == ott[2])
			{
				neighbour[2] = other;
				other->neighbour[1] = this;
			}
			else
				if(vertex[2] == ott[1])
				{
					neighbour[1] = other;
					other->neighbour[2] = this;
				}
				else
					if(vertex[2] == ott[2])
					{
						neighbour[1] = other;
						other->neighbour[1] = this;
					}
	}
	else
		if(vertex[0] == ott[1])
		{
			if(vertex[1] == ott[2])
			{
				neighbour[2] = other;
				other->neighbour[0] = this;
			}
			else
				if(vertex[2] == ott[2])
				{
					neighbour[1] = other;
					other->neighbour[0] = this;
				}
		}
		else
			if(vertex[1] == ott[0])
			{
				if(vertex[2] == ott[1])
				{
				neighbour[0] = other;
				other->neighbour[2] = this;
				}
				else
					if(vertex[2] == ott[2])
					{
					neighbour[0] = other;
					other->neighbour[1] = this;
					}
			}
			else
				if((vertex[1] == ott[1]) && (vertex[2] == ott[2]))
					{
					neighbour[0] = other;
					other->neighbour[0] = this;
					}
				else
					return -1; //join impossible;
//check if complete, return proper value
	return ((neighbour[0] && neighbour[1] && neighbour[2])?1:0) |
			((other->neighbour[0] && other->neighbour[1] && other->neighbour[2])?2:0);
}