/* 
A "chunk" represents a portion of the FEM mesh, managing and
displaying it.

Orion Sky Lawlor, olawlor@acm.org, 12/1/2001
*/
#include <string.h>
#include <stdlib.h>
#include "main.h"

//List of chunks (a very low-level API)
chunk *chunks::chunkList;
int chunks::n; //Number of chunks
void chunks::init(int nChunks)
{
	n=nChunks;
	chunkList=new chunk[n];
	for (int i=0;i<n;i++)
		chunkList[i].setNumber(i);
}

//------------ chunk ----------------
chunk::chunk(void) {
	u=NULL;
	norm=vtkPolyDataNormals::New();
	map=vtkPolyDataMapper::New();
	map->SetInput(norm->GetOutput());
}

chunk::~chunk() {
	delete u;
	norm->Delete();
	map->Delete();
}

void chunk::setUpdate(NetFEM_update *uNew)
{
	delete u;
	u=uNew;
}

void chunk::addGeometry(outputFilter *o)
{
	vtkPolyData *poly=o->addGeometry(*this);
	//norm->SetInput(poly);
	map->SetInput(poly);
	poly->Delete();
}
void chunk::updateRange(outputFilter *o)
{
	o->updateRange(*this).apply(map);
}

void floatRange::apply(vtkPolyDataMapper *map) const
{
	map->SetScalarRange((float *)r);
}

//-------------------- Displacement Filter: ---------------
dispFilter *dispFilter::cur;
dispFilter::~dispFilter() {}

void nodeDispFilter::setUpdate(const NetFEM_update *u)
{
	const NetFEM_doubleField &f=u->getNodes().getField(fieldNo);
	stride=f.getStride();
	data=f.getData(0);
}

//------------ geometry filters ---------
geomFilter *geomFilter::cur;

class pointSource {
protected:
	const double *coor;
	dispFilter *df;
	double scale;
public:
	pointSource(const double *coor_,dispFilter *df_) 
	  :coor(coor_), df(df_) {scale=df->getScale();}
	virtual ~pointSource() {}
	virtual double *get(int pointNo,double *sto) =0;
};
class pointSource2D : public pointSource {
public:
	pointSource2D(const double *coor_,dispFilter *df_)
	   :pointSource(coor_,df_) {}
	virtual double *get(int pointNo,double *sto) {
		const double *disp=df->getDisp(pointNo);
		sto[0]=disp[0]*scale+coor[2*pointNo+0];
		sto[1]=disp[1]*scale+coor[2*pointNo+1];
		sto[2]=0.0;
		return sto;
	}
};
class pointSource3D : public pointSource {
public:
	pointSource3D(const double *coor_,dispFilter *df_)
	   :pointSource(coor_,df_) {}
	virtual double *get(int pointNo,double *sto) {
		const double *disp=df->getDisp(pointNo);
		sto[0]=disp[0]*scale+coor[3*pointNo+0];
		sto[1]=disp[1]*scale+coor[3*pointNo+1];
		sto[2]=disp[2]*scale+coor[3*pointNo+2];
		return sto;
	}
};
pointSource *buildPointSource(int nDim,const double *coor,dispFilter *df)
{
	if (nDim==2) return new pointSource2D(coor,df);
	else if (nDim==3) return new pointSource3D(coor,df);
	else bad("Only 2 and 3 dimentional computations supported!\n");
	return NULL;
}

vtkPolyData *nodeGeomFilter::addGeometry(chunk &c,int &nCopies)
{
	const NetFEM_update *u=c.getUpdate();
	const NetFEM_nodes &n=u->getNodes();
	int nNodes=n.getItems();
	
	vtkPolyData *poly=vtkPolyData::New();
	nCopies=1;
	
	//Build the polygon points
	vtkPoints *pts=vtkPoints::New();
	pts->Allocate(nNodes);
	dispFilter *df=dispFilter::get(); df->setUpdate(u);
	pointSource *p=buildPointSource(u->getDim(),n.getCoord().getData(0),df);
	double sto[3];
	for (int i=0;i<nNodes;i++)
		pts->InsertPoint(i,p->get(i,sto));
	delete p;
	poly->SetPoints(pts);
	pts->Delete();
	
	addConnectivity(poly,nNodes,c);
	
	return poly;
}

void cloudGeomFilter::addConnectivity(vtkPolyData *poly,int nNodes,chunk &c)
{
	//Build the node vertices
	vtkCellArray *faces=vtkCellArray::New();
	faces->Allocate(faces->EstimateSize(nNodes,1));
	for (int i=0;i<nNodes;i++)
		faces->InsertNextCell(1,&i);
	poly->SetVerts(faces);
	faces->Delete();
}

void triGeomFilter::addConnectivity(vtkPolyData *poly,int nNodes,chunk &c)
{
	//Build the polygon faces
	vtkCellArray *faces=vtkCellArray::New();
	const NetFEM_elems &e=c.getUpdate()->getElem(elemType);
	int nEl=e.getItems();
	int nPer=e.getNodesPer();
	const int *conn=e.getNodes(0);
//	if (nPer!=3) bad("Can only plot triangular elements!\n");
	faces->Allocate(faces->EstimateSize(nEl,nPer));
	for (int i=0;i<nEl;i++)
		faces->InsertNextCell(nPer,(int *)&conn[nPer*i]);
	poly->SetPolys(faces);
	faces->Delete();
}

void outlineGeomFilter::addConnectivity(vtkPolyData *poly,int nNodes,chunk &c)
{
	//Build the polygon faces
	vtkCellArray *faces=vtkCellArray::New();
	const NetFEM_elems &e=c.getUpdate()->getElem(elemType);
	int nEl=e.getItems();
	int nPer=e.getNodesPer();
	const int *conn=e.getNodes(0);
	faces->Allocate(faces->EstimateSize(nEl*nPer,2));
	for (int i=0;i<nEl;i++) {
		const int *iconn=&conn[nPer*i];
		int oconn[2];
		for (int j=0;j<nPer;j++) {
			oconn[0]=iconn[j];
			oconn[1]=iconn[(j+1)%nPer];
			faces->InsertNextCell(2,oconn);
		}
	}
	poly->SetLines(faces);
	faces->Delete();
}


vtkPolyData *elemGeomFilter::addGeometry(chunk &c,int &nCopies)
{
	const NetFEM_update *u=c.getUpdate();
	const NetFEM_nodes &n=u->getNodes();
	const NetFEM_elems &e=u->getElem(elemType);
	int nElem=e.getItems();
	int nPer=e.getNodesPer();
	
	vtkPolyData *poly=vtkPolyData::New();	
//Silly VTK doesn't support polygon-valued functions, so we simulate them
// by building a separate copy of each element's surrounding points,
// and then setting all their values equal.
	nCopies=nPer;
	
	//Build the polygon points and connectivity simultaniously
	vtkPoints *pts=vtkPoints::New();
	pts->Allocate(nElem*nPer);
	vtkCellArray *faces=vtkCellArray::New();
	faces->Allocate(faces->EstimateSize(nElem,nPer));
	dispFilter *df=dispFilter::get(); df->setUpdate(u);
	pointSource *p=buildPointSource(u->getDim(),n.getCoord().getData(0),df);
	double sto[3];
	for (int ne=0;ne<nElem;ne++) {
		const int *iconn=e.getNodes(ne);
		int oconn[20]; //Assumes nPer <= 20
		for (int i=0;i<nPer;i++)
			oconn[i]=pts->InsertNextPoint(p->get(iconn[i],sto));
		faces->InsertNextCell(nPer,oconn);
	}
	delete p;
	poly->SetPoints(pts); pts->Delete();
	poly->SetPolys(faces); faces->Delete();
	
	return poly;
}


//----------- field locators ----------
fieldLoc::~fieldLoc() {}

const NetFEM_doubleField &nodeFieldLoc::getField(const NetFEM_update *u) const
{
	return u->getNodes().getField(fieldNo);
}

const NetFEM_doubleField &elemFieldLoc::getField(const NetFEM_update *u) const
{
	return u->getElem(elemType).getField(fieldNo);
}

//----------- field filters ------------
fieldFilter::~fieldFilter() {}

void scalarFieldFilter::setScalars(vtkScalars *vColor,int n,int copies,chunk &c)
{
	int stride=loc->getStride(c.getUpdate());
	const double *srcVal=loc->getData(c.getUpdate());
	
	for (int i=0;i<n;i++) 
	for (int j=0;j<copies;j++)
	{
		vColor->InsertNextScalar(srcVal[stride*i]);
	}
}
void magFieldFilter::setScalars(vtkScalars *vColor,int n,int copies,chunk &c)
{
	int stride=loc->getStride(c.getUpdate());
	const double *srcVal=loc->getData(c.getUpdate());
	
	if (stride==2) //2D vector magnitude
	  for (int i=0;i<n;i++) 
	  for (int j=0;j<copies;j++)
	  {
		const double *v=&srcVal[stride*i];
		vColor->InsertNextScalar((float)sqrt(v[0]*v[0]+v[1]*v[1]));
	  }
	else
	if (stride==3) //3D vector magnitude
	  for (int i=0;i<n;i++) 
	  for (int j=0;j<copies;j++)
	  {
		const double *v=&srcVal[stride*i];
		vColor->InsertNextScalar((float)sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2]));
	  }
}

//Colors chunks by their chunk number
void patchworkFieldFilter::setScalars(vtkScalars *vColor,int n,int copies,chunk &c)
{
	float color=1+c.getNumber();
	for (int i=0;i<n;i++) 
	for (int j=0;j<copies;j++)
	{
		vColor->InsertNextScalar(color);
	}
}


//----------- output filters -----------
// Abstract base class:
float floatRange::BIG=1.0e30;
outputFilter *outputFilter::cur=NULL;

outputFilter::~outputFilter() {}

void outputFilter::beginRendering(void)
{
	r.reset();
}
const floatRange &outputFilter::updateRange(chunk &c)
{
	return r;
}

vtkPolyData *nodeOutputFilter::addGeometry(chunk &c)
{
//Build geometry from update:
	int nCopies=1;
	vtkPolyData *poly=geomFilter::get()->addGeometry(c,nCopies);
	
	//Build the polygon scalars
	vtkScalars *skl=vtkScalars::New();
	int nNodes=c.getUpdate()->getNodes().getItems();
	skl->Allocate(nNodes);
	f->setScalars(skl,nNodes,nCopies,c);
	poly->GetPointData()->SetScalars(skl);
	float rng[2];
	skl->GetRange(rng);
	addToRange(rng);	
	skl->Delete();
	
	return poly;
}

vtkPolyData *elemOutputFilter::addGeometry(chunk &c)
{
//Build geometry from update:
	int nCopies=1;
	vtkPolyData *poly=geomFilter::get()->addGeometry(c,nCopies);
	
	//Build the polygon scalars
	vtkScalars *skl=vtkScalars::New();
	int n=c.getUpdate()->getElem(elemType).getItems();
	skl->Allocate(n);
	f->setScalars(skl,n,nCopies,c);
	poly->GetPointData()->SetScalars(skl);
	float rng[2];
	skl->GetRange(rng);
	addToRange(rng);	
	skl->Delete();
	
	return poly;
}


