//	-*- Mode: C++ -*-
//	navi.h	 Copyright (C) 1990-1997,1999 Yukinori NAGATANI
//	
//	The Program for culculation
//	Navier-Stokes equation
//	of fuild motion
//
//	Class Module
//	
//	Aug,24,1990	Navi2
//	Mar, 3,1995	Navi on X Windows
//	Nov,10,1997	C++ version
//	Nov,10,1999	C++ version (ppm)

#include <stdio.h>
#include <iostream.h>
#include <fstream.h>
#include <stdlib.h>	// atoi atof
#include <math.h>

typedef	double	Real;
typedef	char	*String;

#ifndef Bool
typedef	int	Bool;
const Bool	True = 1,	False = 0;
#endif

String	szMessage =
	"navi:Illegal option error!\n"
	"Available options : \n"
	"    [-xsize ##] [-ysize ##] [-load  FILE_NAME] [-save FILE_NAME]\n";

void Error (String szMessage)
{
	cerr << szMessage << endl;
	exit (0);
}

#define	TIME_RECORD
#ifdef TIME_RECORD
#include <time.h>	// time related functions
// Get Formated Time String
String getTimeString (void)
{
	static char	szTimeString [256];
	time_t		t;
	
	time (&t);
	strftime (szTimeString, 256,
		  "%Y/%m/%d %a %H:%M:%S",
		  localtime (&t));
	
	return	szTimeString;
}
#endif

class Point {
public:
	enum BoundaryType {
		FREE		= 0,
		IN		= 1,
		OUT		= 2,
		FREE_SURFACE	= 3,
		SOLID		= 4,
		UP_SURFACE	= 5,
		DOWN_SURFACE	= 6,
		RIGHT_SURFACE	= 7,
		LEFT_SURFACE	= 8
	};
	
	Real		z;	// rotation
	Real		f;	// stream function
	BoundaryType	b;	// boundary condition
	
	Point (void)
	{
		z = 0.0;
		f = 0.0;
		b = FREE;
	}
	friend ostream & operator << (ostream& s, const Point &p)
	{
		s << p.z << '	'
		  << p.f << '	'
		  << p.b << endl;
		return	s;
	}
	friend istream & operator >> (istream& s, Point &p)
	{
		s >> p.z
		  >> p.f
		  >> (int &) p.b;
		return	s;
	}
};

class Field {
public:
	int	xMax,	yMax;	// Matrix Size
	Point	*field;
	
	Field (void)
	{
		xMax = yMax = 0;
		field = NULL;
	}
	inline Point & value (int x, int y) const
	{
		// Caution !! NO DOMAIN-CHECKING of (x,y) for speed-up.
		return	field [x + y*xMax];
	}
	Field & setSize (int xMax, int yMax)
	{
		if (this->xMax == xMax &&
		    this->yMax == yMax &&
		    field != NULL)
			return	*this;
		
		delete field;
		this->xMax = xMax;
		this->yMax = yMax;
		field = new Point [yMax * xMax];

		return	*this;
	}
	friend ostream & operator << (ostream& s, const Field &f)
	{
		s << ";; Field Data" << endl
		  << "size	" << f.xMax << "	" << f.yMax << endl
		  << "# main data" << endl;
		
		for (int y = 0 ; y < f.yMax ; y++) {
			for (int x = 0 ; x < f.xMax ; x++) {
				s << f.value (x,y);
			}
		}
		return	s;
	}
	friend istream & operator >> (istream& s, Field &f)
	{
		char	szBuffer[256];
		for (;;) {
			if (s.getline (szBuffer, 256) == NULL) {
				// s.close ();
				return s;	// End of File
			}
			if (szBuffer[0] == ';')	// Comment
				continue;
			if (szBuffer[0] == 's')
				sscanf (szBuffer, "%*s %d %d",
					&(f.xMax), &(f.yMax));
			if (szBuffer[0] == '#')	// data start
				break;
		}
		f.setSize (f.xMax, f.yMax);
		for (int y = 0 ; y < f.yMax ; y++) {
			for (int x = 0 ; x < f.xMax ; x++) {
				s >> f.value (x,y);
			}
		}
		return	s;
	}
	Field & operator = (const Field & f)	// substitution
	{
		if (this == &f)
			return	*this;
		if (xMax != f.xMax || yMax != f.yMax)
			Error ("Size Unmach in Operator Field & = (const Field &)");
	
		for (int y = 0 ; y < yMax ; y++) {
			for (int x = 0 ; x < xMax ; x++) {
				value (x,y) = f.value (x,y);
			}
		}
		return	*this;
	}
	void swap (Field & f)
	{
		if (xMax != f.xMax || yMax != f.yMax)
			Error ("Size Unmach in swap");
		
		Point	* swapPoint = field;
		field	= f.field;
		f.field	= swapPoint;
	}
};

class NavierStokes {
public:
	Field		fldCurrent,	fldNext;
	
	int		xMax,	yMax;	// Matrix Size
	
	Real		dx,	dy;	// X_axis and Y_axis differential
	Real		dt;		// Time differential
	long int	nTime;		// Time
	
	Real		vx;		// Velocity ( X_axis )
	Real		nu;		// stickiness coefficient
	
	String		szLoadFileName;	// Load File Name
	ifstream	fsInput;
	Bool		bfsInputOpened;
	
	String		szSaveFileName;	// Periodic Save File Name
	int		nPeriod;	// Periode of data file save
	
	NavierStokes (void)
	{
		xMax = 64 * 2;
		yMax = 64;
		
		dx = dy = 1.0 / yMax;	// Height of Wind-Tunnel is 1.0
		vx = 1.0;		// Flow Velocity
		nu = 0.001;		// stickiness coefficient
		// L := Object Size in Wind-Tunnel
		// Reynolds number := vx * L / nu = L * 1000
		
		dt = 0.0005;		// Time Step
		nTime = 0;
		
		szLoadFileName	= "navi.dat";
		bfsInputOpened	= False;
		
		szSaveFileName	= "navi.dat";
		nPeriod = 200;
	}
	setParams (int nArgc, String *pszArgv)
	{
		pszArgv++;
		for (; nArgc > 1  && (*pszArgv)[0] == '-' ; nArgc--, pszArgv++) {
			if (!strcmp (*pszArgv, "-xsize")) {
				nArgc--;
				pszArgv++;
				if (nArgc <= 0)
					Error (szMessage);
				xMax = atoi (*pszArgv);
				continue;
			}
			if (!strcmp (*pszArgv, "-ysize")) {
				nArgc--;
				pszArgv++;
				if (nArgc <= 0)
					Error (szMessage);
				yMax = atoi (*pszArgv);
				continue;
			}
			if (!strcmp (*pszArgv, "-load")) {
				nArgc--;
				pszArgv++;
				if (nArgc <= 0)
					Error (szMessage);
				szLoadFileName = *pszArgv;
				continue;
			}
			if (!strcmp (*pszArgv, "-save")) {
				nArgc--;
				pszArgv++;
				if (nArgc <= 0)
					Error (szMessage);
				szSaveFileName = *pszArgv;
				continue;
			}
			Error (szMessage);
		}
		
		fldCurrent.setSize (xMax, yMax);
		fldNext   .setSize (xMax, yMax);
		
		setInitialize (fldCurrent);
		setInitialize (fldNext);
		
	}
	Bool SaveToFile (const String szSaveFileName = NULL) const
	{
		String	szFileName = 
			(szSaveFileName != NULL ?
			 szSaveFileName : this->szSaveFileName);
		
		// Add data to Text file
		ofstream	fsOut (szFileName, ios::app);
		if (fsOut == NULL) {
			cerr << "navi:Cannot Open File <" << szFileName 
			     <<  "> for Save." << endl;
			return False;
		}
		fsOut << ";; NAVI_DATA" << endl;
#ifdef TIME_RECORD
		fsOut << ";; Date " << getTimeString () << endl;
#endif
		fsOut <<  "time	" << nTime << endl
		      <<  "dt	" << dt	   << endl
		      << fldCurrent
		      << endl;
		fsOut.close ();
		return	True;
	}
	Bool ReadFileOneUnit (void)
	{
		if (!bfsInputOpened) {
			// Open Text file
			fsInput.open (szLoadFileName, ios::in);
			if (fsInput == NULL) {
				cerr << "navi:Cannot Open File <"
				     << szLoadFileName
				     << "> for Load." << endl;
				return False;
			}
			bfsInputOpened = True;
		}
		
		char		szBuffer[256];
		for (;;) {
			if (fsInput.getline (szBuffer, 256) == NULL) {
				// End of File
				fsInput.close ();
				bfsInputOpened = False;
				return False;
			}
			if (szBuffer[0] == ';')	// Comment
				continue;
			if (szBuffer[0] == 't')
				sscanf (szBuffer, "%*s %ld", &nTime);
			if (szBuffer[0] == 'd') {
				sscanf (szBuffer, "%*s %lf", &dt);
				break;
			}
		}
		fsInput >> fldCurrent;
		xMax = fldCurrent.xMax;
		yMax = fldCurrent.yMax;
		return True;
	}
	int nReadFileLastUnit ()
	{
		int	nUnit;
		for (nUnit = 0 ; ReadFileOneUnit () ; nUnit++);
		return	nUnit;
	}
	
	void setInitialize (Field & field) const
	{
		int	xMax = field.xMax,
			yMax = field.yMax;
		for (int y = 0 ; y < yMax ; y++){
		    for (int x = 0 ; x < xMax ; x++) {
			
			Point	& p = field.value (x,y);
			    
			p.b  = Point::FREE;
			p.f  = y*vx*dy;
			p.z  = 0;
			
			int	xB1 = (int)(0.225*xMax), yB1 = (int)(0.45*yMax),
				xB2 = (int)(0.275*xMax), yB2 = (int)(0.55*yMax);
			Real	SurfaceFlow = (0.5*(yMax-1))*vx*dy;
			
			// Flow Input & Output
			if (x == 0) {
				p.b = Point::IN;
				p.f = y*vx*dy;
				p.z = 0;
			}
			if (x == xMax-1) {
				p.b = Point::OUT;
				p.f = y*vx*dy;
				p.z = 0;
			}
			// Flow Guid Surface
			if (y == 0 || y == yMax-1) {
				p.b = Point::FREE_SURFACE;
			}
			// Solid
			if (xB1 <  x && x <  xB2 &&
			    yB1 <  y && y <  yB2 ){
				p.b = Point::SOLID;
				p.f = SurfaceFlow;
				p.z = 0;
			}
			// Solid Surface
			if (xB1 <  x  && x <  xB2 &&
			    yB1 == y              ) {
				p.b = Point::DOWN_SURFACE;
				p.f = SurfaceFlow;
			}
			if (xB1 <  x && x <  xB2 &&
			                y == yB2 ) {
				p.b = Point::UP_SURFACE;
				p.f = SurfaceFlow;
			}
			if (xB1 == x &&
			    yB1 <  y && y < yB2) {
				p.b = Point::LEFT_SURFACE;
				p.f = SurfaceFlow;
			}
			if (            x == xB2 &&
			    yB1 <  y && y <  yB2 ) {
				p.b = Point::RIGHT_SURFACE;
				p.f = SurfaceFlow;
			}
		    }
		}
	
	}
	// Get Rotation from Stream Function on Boundaries
	void boundary (Field & field) const
	{
		Real	dxdx = dx*dx,
			dydy = dy*dy;
		
		for (int y = 0 ; y < field.yMax ; y++) {
		    for (int x = 0 ; x < field.xMax ; x++) {
			
			Point	& p = field.value (x,y);
			
			switch (p.b) {
			case Point::DOWN_SURFACE:
				p.z = 2 * ( p.f - field.value (x  ,y-1).f ) / dydy;
				break;
			case Point::UP_SURFACE:
				p.z = 2 * ( p.f - field.value (x  ,y+1).f ) / dydy;
				break;
			case Point::LEFT_SURFACE:	
				p.z = 2 * ( p.f - field.value (x-1,y  ).f ) / dxdx;
				break;
			case Point::RIGHT_SURFACE:
				p.z = 2 * ( p.f - field.value (x+1,y  ).f ) / dxdx;
				break;
			}
		    }
		}
	}
	
	// solve  Poisson equation of stream function
	void SubPoison (Field & fIn, Field & fOut) const
	{
		Real	dxdx = dx*dx,
			dydy = dy*dy,
			dxdxdydy = dxdx*dydy,
			dr2 = 1/(2*(dxdx+dydy));
		
		boundary (fIn);
		for (int y = 0 ; y < fOut.yMax ; y++) {
		    for (int x = 0 ; x < fOut.xMax ; x++) {
			if (fIn.value(x,y).b == Point::FREE) {
				
				fOut.value(x,y).f = dr2 * (
					  dxdx * (   fIn .value (x  ,y+1).f
						   + fOut.value (x  ,y-1).f )
					+ dydy * (   fOut.value (x-1,y  ).f
						   + fIn .value (x+1,y  ).f )
					+ dxdxdydy * fIn .value (x  ,y  ).z  );
			}
		    }
		}
	}
	void poison (Field & fCurr, Field & fWork) const
	{
		fWork = fCurr;	// copy
		
		for (int k = 0 ; k < 20 ; k++) {
			SubPoison (fCurr, fWork);
			SubPoison (fWork, fCurr);
		}
	}
	void CalculationStepSub (Field & fIn, Field & fOut) const
	{
		Real	dtp4dxdy  = dt / (4*dx*dy),
			dtnupdxdx = dt * nu / (dx*dx),
			dtnupdydy = dt * nu / (dy*dy);
		
		// rotation boundary condition
		boundary (fIn);
		
		// solve rotation transforn equation
		for (int y = 0 ; y < fOut.yMax ; y++) {
		    for (int x = 0 ; x < fOut.xMax ; x++) {
			if (fIn.value (x,y).b == Point::FREE) {
			
				fOut.value (x,y).z
				    = fIn.value (x,y).z 
				    -dtp4dxdy*(( fIn.value (x  ,y+1).f
						-fIn.value (x  ,y-1).f )
					      *( fIn.value (x+1,y  ).z
						-fIn.value (x-1,  y).z )
					     - ( fIn.value (x+1,y  ).f
						-fIn.value (x-1,  y).f )
					      *( fIn.value (x  ,y+1).z
						-fIn.value (x  ,y-1).z ) )
				    +dtnupdxdx * (    fIn.value (x-1,y  ).z
						   -2*fIn.value (x  ,y  ).z
						   +  fIn.value (x+1,y  ).z )
				    +dtnupdydy * (    fIn.value (x  ,y-1).z
						   -2*fIn.value (x,y    ).z
						   +  fIn.value (x  ,y+1).z );
			}
		    }
		}
		// solve Poisson equation of stream function
		// both sorce and resul are in fOut,
		// fIn is only work area.
		poison (fOut, fIn);
	}
	void CalculationStep (void)
	{
		CalculationStepSub (fldCurrent, fldNext);
		fldCurrent .swap (fldNext);
		nTime ++;
	}
	
	Bool bJustOnPeriod (void) const
	{
		return	(nTime % (long int) nPeriod == 0);
	}
	Bool SaveStep (void)
	{
		Bool	bSaved = False;
		if (bJustOnPeriod()) {
			SaveToFile ();
			bSaved = True;
		}
		
		// CalculationStep ();
		
		return	bSaved;
	}
};