#include "udf.h"

int ncount;
int npost;

/*You need to modify these three lines*/
int POST_BOUNDARY_ZONE_ID[50]; 
int POST_CELL_ZONE_ID_1[50];
int POST_CELL_ZONE_ID_2[50]; 

float L_XX[50]; 
float L_YX[50]; 
float L_ZX[50]; 

float L_XY[50]; 
float L_YY[50]; 
float L_ZY[50]; 

float L_XZ[50]; 
float L_YZ[50]; 
float L_ZZ[50]; 


static void read_input()
{
	int i=0;
	FILE *fpin;
	
	if((fpin=fopen("input.txt","r"))==NULL)
	{
		printf("Input file does not exist!");
	}
	fscanf(fpin, "%d",&npost);
	for(i=0;i<npost;i++)
	{
		fscanf(fpin,"%d %d %d %f %f %f %f %f %f", &POST_BOUNDARY_ZONE_ID[i], &POST_CELL_ZONE_ID_1[i],
				&POST_CELL_ZONE_ID_2[i], &L_XX[i], &L_YX[i], &L_ZX[i], &L_XY[i], &L_YY[i], &L_ZY[i]);

		L_XZ[i]=L_YX[i]*L_ZY[i]-L_ZX[i]*L_YY[i];
		L_YZ[i]=L_ZX[i]*L_XY[i]-L_XX[i]*L_ZY[i];
		L_ZZ[i]=L_XX[i]*L_YY[i]-L_XY[i]*L_YX[i];
	}
	fclose(fpin);
}

static void init_udm()
{

  	Thread *tc;
  	cell_t c;

  	Domain *domain;
	domain = Get_Domain(1);  	

	thread_loop_c(tc,domain)
    	{
      		begin_c_loop(c,tc)
		{
	  		C_UDMI(c,tc,0)=0;
	  		C_UDMI(c,tc,1)=0;
	  		C_UDMI(c,tc,2)=0;
	  		C_UDMI(c,tc,3)=0;
		}
      		end_c_loop(c,tc)
	}
  	ncount=0;
}



static void cat_post(int flag)
{

  	float xf[ND_ND], x_vmax_global[ND_ND], x_vmax, y_vmax;
  	float x_face_local, y_face_local;
  	float xn[ND_ND],x_node_local,y_node_local;
  	float x_min, x_max, x_mid, y_min, y_max, y_mid;
  	float L_major, L_minor;
  	float eccent_x, eccent_y, eccent;
  	float vmag, vmax, vavg, vratio, vel_space; 
  	float A_face_vec[ND_ND], A_face_mag, A_tot;
  	float a_35, v_35, unif_35, a_65, v_65, unif_65, gamma;
  	int k;

  	face_t f, f_vmax;
  	cell_t c0, c1;
  	Thread *tf, *tc0, *tc1;
	
	FILE *fp;

  	Domain *domain;
	domain = Get_Domain(1);  	


	#define BIGNUM 1e20
	
	if((fp=fopen("output.txt", "r"))==NULL)
	{
		fp=fopen("output.txt","w+");
		fprintf(fp, "post-face    eccent     gamma     unif_35     unif_65     vratio     vmax \n");
		fclose(fp);
	} 

	fp=fopen("output.txt", "a");

	/*  find local y,z min/max extents based on node values */

  	x_min=BIGNUM;
  	y_min=BIGNUM;
  	x_max=-BIGNUM;
  	y_max=-BIGNUM;
  	tf=Lookup_Thread(domain, POST_BOUNDARY_ZONE_ID[flag]);
  	begin_f_loop(f,tf)
  	{
    		f_node_loop(f,tf,k)
    		{
      		xn[0]=NODE_X(F_NODE(f,tf,k));
      		xn[1]=NODE_Y(F_NODE(f,tf,k));
		      xn[2]=NODE_Z(F_NODE(f,tf,k));
		      x_node_local = NVD_DOT(xn,L_XX[flag],L_YX[flag],L_ZX[flag]);
		      y_node_local = NVD_DOT(xn,L_XY[flag],L_YY[flag],L_ZY[flag]);
		      x_min=MIN(x_node_local,x_min);
		      x_max=MAX(x_node_local,x_max);
		      y_min=MIN(y_node_local,y_min);
		      y_max=MAX(y_node_local,y_max);
		}
  	}
  	end_f_loop(f,tf)
  	x_mid=0.5*(x_min + x_max); 
  	y_mid=0.5*(y_min + y_max); 
  	L_major = x_max - x_min;
  	L_minor = y_max - y_min;


	Message("================================================================\n");
	Message("post-%d\n",flag);
	Message("================================================================\n");
  	Message("Geometry information:\n");
  	Message("(x_min, x_mid, x_max) = (%f, %f, %f) [m]\n", x_min, x_mid, x_max);
  	Message("(y_min, y_mid, y_max) = (%f, %f, %f) [m]\n", y_min, y_mid, y_max);
  	Message("(L_major, L_minor) = (%f, %f) [m]\n\n", L_major, L_minor);

	/* velocity index (eccentricity) computation */
  	vmax=0.0;
  	vavg=0.0;
  	A_tot=0.0;
  	tf=Lookup_Thread(domain, POST_BOUNDARY_ZONE_ID[flag]);
	tc0=Lookup_Thread(domain, POST_CELL_ZONE_ID_1[flag]);
  	tc1=Lookup_Thread(domain, POST_CELL_ZONE_ID_2[flag]);
  	begin_f_loop(f,tf)
  	{
    		F_AREA(A_face_vec,f,tf);
    		A_face_mag=NV_MAG(A_face_vec);
    		A_tot+=A_face_mag;
    		c0=F_C0(f,tf);
    		c1=F_C1(f,tf);
		vmag=0.5*(C_UDMI(c0,tc0,3)+C_UDMI(c1,tc1,3));
    		vavg+=vmag*A_face_mag;
    		if(vmag > vmax)
    		{
      		vmax=vmag;
      		f_vmax=f;
    		}
  	}
  	end_f_loop(f,tf)
  	vavg/=A_tot;
  	vratio=vmax/vavg;
  	F_CENTROID(x_vmax_global,f_vmax,tf);
  	x_vmax= NVD_DOT(x_vmax_global,L_XX[flag],L_YX[flag],L_ZX[flag]);
  	y_vmax= NVD_DOT(x_vmax_global,L_XY[flag],L_YY[flag],L_ZY[flag]);

  	eccent_x=2.0*(x_vmax-x_mid)/L_major;
  	eccent_y=2.0*(y_vmax-y_mid)/L_minor;
	eccent=sqrt(eccent_x*eccent_x+eccent_y*eccent_y);

  	Message("Eccentricity information:\n");
  	Message("(x_vmax, y_vmax) = (%f, %f) [m]\n", x_vmax, y_vmax);
  	Message("eccentricity = (%f, %f) ->  %f \n\n",eccent_x, eccent_y, eccent); 
  	Message("vmax = %f m/s occurs at face %d\n", vmax, f_vmax);

  	Message("vratio = %f\n", vratio);

	/* uniformity, gamma calculations */
  	v_35=0.35*vmax;
  	v_65=0.65*vmax;
  	a_35=0.0;
  	a_65=0.0;
  	gamma=0.0;
  	begin_f_loop(f,tf)
  	{
    		F_AREA(A_face_vec,f,tf);
    		A_face_mag=NV_MAG(A_face_vec);
    		c0=F_C0(f,tf);
    		c1=F_C1(f,tf);
     		vmag=0.5*(C_UDMI(c0,tc0,3)+C_UDMI(c1,tc1,3));
    		if(vmag >= v_35) a_35+=A_face_mag;
    		if(vmag >= v_65) a_65+=A_face_mag;
    		gamma+=ABS(vmag-vavg)*A_face_mag;
  	}
  	end_f_loop(f,tf)
  	gamma=1.0-gamma/(2.*vavg*A_tot);

  	unif_35=a_35/A_tot*100;
  	unif_65=a_65/A_tot*100;
  
  	Message("uniformity index = %f/%f\n",unif_65,unif_35);
  	Message("gamma = %f\n",gamma);

	
	fprintf(fp, "post%d    %9.3f  %9.3f  %9.3f   %9.3f  %9.3f  %9.3f \n", flag, eccent, gamma, unif_35, unif_65, vratio, vmax);
	fclose(fp);
	
}

static void set_udm(int flag)
{
	/* NOTE: THIS ROUTINE REQUIRES ALLOCATION OF */
	/* 4 USER DEFINED MEMORY LOCATIONS IN FLUENT */

	#define X_PRIME 0
	#define Y_PRIME 1
	#define Z_PRIME 2
	#define V_MAG 3

  	float xc[ND_ND];
  	Thread *tc1,*tc2;
  	cell_t c;

  	Domain *domain;
	domain = Get_Domain(1);  	


  	tc1=Lookup_Thread(domain, POST_CELL_ZONE_ID_1[flag]);
  	tc2=Lookup_Thread(domain, POST_CELL_ZONE_ID_2[flag]);

  	begin_c_loop(c,tc1)
  	{
    		C_CENTROID(xc,c,tc1);
    		C_UDMI(c,tc1,X_PRIME)=NVD_DOT(xc,L_XX[flag],L_YX[flag],L_ZX[flag]);
    		C_UDMI(c,tc1,Y_PRIME)=NVD_DOT(xc,L_XY[flag],L_YY[flag],L_ZY[flag]);
    		C_UDMI(c,tc1,Z_PRIME)=NVD_DOT(xc,L_XZ[flag],L_YZ[flag],L_ZZ[flag]);
    		C_UDMI(c,tc1,V_MAG)=sqrt(C_VMAG2(c,tc1));
  	}
  	end_c_loop(c,tc1)

  	if(POST_CELL_ZONE_ID_1[flag] != POST_CELL_ZONE_ID_2[flag])
  	{
    		begin_c_loop(c,tc2)
    		{
    			C_CENTROID(xc,c,tc2);
    			C_UDMI(c,tc2,X_PRIME)=NVD_DOT(xc,L_XX[flag],L_YX[flag],L_ZX[flag]);
    			C_UDMI(c,tc2,Y_PRIME)=NVD_DOT(xc,L_XY[flag],L_YY[flag],L_ZY[flag]);
    			C_UDMI(c,tc2,Z_PRIME)=NVD_DOT(xc,L_XZ[flag],L_YZ[flag],L_ZZ[flag]);
    			C_UDMI(c,tc2,V_MAG)=sqrt(C_VMAG2(c,tc2));
    		}
  		end_c_loop(c,tc2)
  	}
}



DEFINE_ON_DEMAND(Uniform_steady)
{
	int i=0;
	
	read_input();

	init_udm();
	for(i=0; i<npost;i++)
	{
		set_udm(i);
		cat_post(i);
	}
}
