#ifndef _POLY_WRITER_H
#define _POLY_WRITER_H
#include "head.h"

class poly_writer
{
public:
	poly_writer();
	~poly_writer();
	int readdat(string);
	int poly_plane(string);
	int poly_topo(string);
	int run(string,string,int,double,double*);
	
private:
	int M,N;
	double bottom_depth;
	double xmin,xmax,dx,ymin,ymax,dy,zmin,zmax;
	double *topo;
};

poly_writer::poly_writer()
{
	topo = NULL;
}

poly_writer::~poly_writer()
{
	if(topo!=NULL) delete[] topo;
}

int poly_writer::run(string inname,string outname,int bottom_type,double depth,double* Range)
{
	bottom_depth = depth;
	xmin = *Range; dx = *(Range+1); xmax = *(Range+2);
	ymin = *(Range+3); dy = *(Range+4); ymax = *(Range+5);
	M = (xmax-xmin)/dx+1;
	N = (ymax-ymin)/dy+1;
	xmax = xmin+(M-1)*dx;
	ymax = ymin+(N-1)*dy;
	if (readdat(inname)) return 1;
	if (bottom_type==1)
	{
		if (poly_topo(outname)) return 1;
	}
	else if (bottom_type==2)
	{
		if (poly_plane(outname)) return 1;
	}
	return 0;
}

int poly_writer::readdat(string filename)
{
	const char* openname = filename.c_str();
	ifstream infile(openname);
	if (!infile)
	{
		cout<<filename<<" not found, program stopped..."<<endl;
		return 1;
	}
	else
	{
		double temp;
		topo = new double [M*N];
		zmin = 1e20; zmax = -1e20;
		for (int i = 0; i < M*N; i++)
		{
			infile>>temp>>temp>>topo[i];
			if (zmin>topo[i]) zmin = topo[i];
			if (zmax<topo[i]) zmax = topo[i];
		}
		infile.close();
		return 0;
	}
}

int poly_writer::poly_plane(string filename)
{
	const char* savename = filename.c_str();
	ofstream outfile(savename);
	if (!outfile)
	{
		cout<<"can not create "<<filename<<endl;
		return 1;
	}
	else
	{
		outfile<<"# Part 1 - node list"<<endl;
		outfile<<"# node count, 3 dim, no attribute, no boundary marker"<<endl;
		outfile<<M*N+4<<" "<<3<<" "<<0<<" "<<0<<endl;
		outfile<<"# Node index, node coordinates"<<endl;
		for(int i=0;i<M*N;i++)
		{
			outfile<<i+1<<" "<<xmin+dx*(i%M)<<" "<<ymin+dy*(i/M)<<" "<<topo[i]<<endl;
		}
		outfile<<M*N+1<<" "<<xmin<<" "<<ymin<<" "<<zmin-bottom_depth<<endl;
		outfile<<M*N+2<<" "<<xmin<<" "<<ymax<<" "<<zmin-bottom_depth<<endl;
		outfile<<M*N+3<<" "<<xmax<<" "<<ymin<<" "<<zmin-bottom_depth<<endl;
		outfile<<M*N+4<<" "<<xmax<<" "<<ymax<<" "<<zmin-bottom_depth<<endl<<endl;

		outfile<<"# Part 2 - facet list"<<endl;
	    outfile<<"# facet count, no boundary marker"<<endl;
		outfile<<(M-1)*(N-1)*2+5<<" "<<0<<endl;
		outfile<<"# 1 polygon, no hole, no boundary marker"<<endl;
		for(int i=1;i<=N-1;i++)
		{
			for(int j=1;j<=M-1;j++)
			{
				outfile<<1<<endl<<3<<" "<<(i-1)*M+j<<" "<<(i-1)*M+j+1<<" "<<i*M+j+1<<endl;
			}
		}
		for(int i=N;i>=2;i--)
		{
			for(int j=M;j>=2;j--)
			{
				outfile<<1<<endl<<3<<" "<<(i-1)*M+j<<" "<<(i-1)*M+j-1<<" "<<(i-2)*M+j-1<<endl;
			}
		}

		outfile<<1<<endl<<4<<" "<<M*N+1<<" "<<M*N+2<<" "<<M*N+4<<" "<<M*N+3<<endl;
		
		outfile<<1<<endl<<N+2<<" ";
		for(int i=1;i<=N;i++)
		{
			outfile<<(i-1)*M+1<<" ";
		}
		outfile<<M*N+2<<" "<<M*N+1<<endl;

		outfile<<1<<endl<<N+2<<" ";
		for(int i=1;i<=N;i++)
		{
			outfile<<i*M<<" ";
		}
		outfile<<M*N+4<<" "<<M*N+3<<endl;

		outfile<<1<<endl<<M+2<<" ";
		for(int i=1;i<=M;i++)
		{
			outfile<<i<<" ";
		}
		outfile<<M*N+3<<" "<<M*N+1<<endl;

		outfile<<1<<endl<<M+2<<" ";
		for(int i=1;i<=M;i++)
		{
			outfile<<(N-1)*M+i<<" ";
		}
		outfile<<M*N+4<<" "<<M*N+2<<endl<<endl;

		outfile<<"# Part 3 - hole list"<<endl;
    	outfile<<"0 # no hole"<<endl;
    	outfile<<"# Part 4 - region list"<<endl;
    	outfile<<"0 # no region"<<endl;
    	outfile<<"#generated by xyz2poly "<<filename;
    	outfile.close();
    	return 0;
	}
}

int poly_writer::poly_topo(string filename)
{
	const char* savename = filename.c_str();
	ofstream outfile(savename);
	if (!outfile)
	{
		cout<<"can not create "<<filename<<endl;
		return 1;
	}
	else
	{
		outfile<<"# Part 1 - node list"<<endl;
		outfile<<"# node count, 3 dim, no attribute, no boundary marker"<<endl;
		outfile<<2*M*N<<" "<<3<<" "<<0<<" "<<0<<endl;
		outfile<<"# Node index, node coordinates"<<endl;
		for(int i=0;i<M*N;i++)
		{
			outfile<<i+1<<" "<<xmin+dx*(i%M)<<" "<<ymin+dy*(i/M)<<" "<<topo[i]<<endl;
		}
		for(int i=0;i<M*N;i++)
		{
			outfile<<M*N+i+1<<" "<<xmin+dx*(i%M)<<" "<<ymin+dy*(i/M)<<" "<<topo[i]-bottom_depth<<endl;
		}

		outfile<<endl<<"# Part 2 - facet list"<<endl;
	    outfile<<"# facet count, no boundary marker"<<endl;
		outfile<<(M-1)*(N-1)*4+4<<" "<<0<<endl;
		outfile<<"# 1 polygon, no hole, no boundary marker"<<endl;
		for(int i=1;i<=N-1;i++)
		{
			for(int j=1;j<=M-1;j++)
			{
				outfile<<1<<endl<<3<<" "<<(i-1)*M+j<<" "<<(i-1)*M+j+1<<" "<<i*M+j+1<<endl;
			}
		}
		for(int i=N;i>=2;i--)
		{
			for(int j=M;j>=2;j--)
			{
				outfile<<1<<endl<<3<<" "<<(i-1)*M+j<<" "<<(i-1)*M+j-1<<" "<<(i-2)*M+j-1<<endl;
			}
		}

		for(int i=1;i<=N-1;i++)
		{
			for(int j=1;j<=M-1;j++)
			{
				outfile<<1<<endl<<3<<" "<<M*N+(i-1)*M+j<<" "<<M*N+(i-1)*M+j+1<<" "<<M*N+i*M+j+1<<endl;
			}
		}
		for(int i=N;i>=2;i--)
		{
			for(int j=M;j>=2;j--)
			{
				outfile<<1<<endl<<3<<" "<<M*N+(i-1)*M+j<<" "<<M*N+(i-1)*M+j-1<<" "<<M*N+(i-2)*M+j-1<<endl;
			}
		}
		
		outfile<<1<<endl<<N*2<<" ";
		for(int i=1;i<=N;i++)
		{
			outfile<<(i-1)*M+1<<" ";
		}
		for(int i=N;i>=1;i--)
		{
			outfile<<M*N+(i-1)*M+1<<" ";
		}
		outfile<<endl;

		outfile<<1<<endl<<N*2<<" ";
		for(int i=1;i<=N;i++)
		{
			outfile<<i*M<<" ";
		}
		for(int i=N;i>=1;i--)
		{
			outfile<<M*N+i*M<<" ";
		}
		outfile<<endl;

		outfile<<1<<endl<<M*2<<" ";
		for(int i=1;i<=M;i++)
		{
			outfile<<i<<" ";
		}
		for(int i=M;i>=1;i--)
		{
			outfile<<M*N+i<<" ";
		}
		outfile<<endl;

		outfile<<1<<endl<<M*2<<" ";
		for(int i=1;i<=M;i++)
		{
			outfile<<(N-1)*M+i<<" ";
		}
		for(int i=M;i>=1;i--)
		{
			outfile<<M*N+(N-1)*M+i<<" ";
		}
		outfile<<endl<<endl;

		outfile<<"# Part 3 - hole list"<<endl;
    	outfile<<"0 # no hole"<<endl;
    	outfile<<"# Part 4 - region list"<<endl;
    	outfile<<"0 # no region"<<endl;
    	outfile<<"#generated by xyz2poly "<<filename;
    	outfile.close();
    	return 0;
	}
}

#endif