/**********************************************************************
   UDF to calculate temperature field function and store in           
   user-defined memory. Also print min, max, avg temperatures.             
***********************************************************************/
#include "udf.h"

DEFINE_ON_DEMAND(on_demand2_calc)
{
   Domain *d;
   /* declare domain pointer since it is not passed as an  
                 //argument to the DEFINE macro  */
   real velocity;
   Thread *t;
   cell_t c;
   Node *nod;
   int n;
   double pos_x;
   double pos_y;
   double xvel;
   double yvel;
   
   double xprev=0;
   double yprev=0;
   
	
   double G_theta = 0;
   double G_theta2 = 0;
   double viz;	//water = 1.00481e-6;
   double viz2 = 1.00481e-6;
   double Cu = 1.9; //ke = 0.09
   double G = 0;
   double theta = 0;
   double cellNo = 0;
   double maxViz = 0;

   d = Get_Domain(1);     /* Get the domain using Fluent utility */

   /* Loop over all cell threads in the domain */
   thread_loop_c(t,d)
     {

     /* Compute max, min, volume-averaged temperature */

     /* Loop over all cells  */
     begin_c_loop(c,t)
       {
		   double delta_x = 0;
			double delta_y = 0;
			int run_number = 0;
		   for (n=0; n < cell_type_nnodes[(int)C_TYPE(c,t)]; n++) 	 
		   
			   //c_node_loop(c,t,n)
				{
					nod = C_NODE(c, t, n); 
					pos_x=NODE_X(nod);
					pos_y=NODE_Y(nod);
					
				if(run_number>0)
				{
					if(delta_x<fabs(pos_x-xprev))
						delta_x= fabs(pos_x-xprev);
					
					if(delta_y < fabs(pos_y-yprev))
						delta_y = fabs(pos_y-yprev);
					
				}
				run_number = run_number +1;
				xprev = pos_x;
				yprev = pos_y;

				}
		// kinematic viscosity, v  = Cu*k^2/e
		   viz = Cu*C_K(c,t)*C_K(c,t)/C_D(c,t) + viz2;
			
			theta = theta + delta_x*delta_y/(fabs(C_U(c,t)*delta_y) + fabs(C_V(c,t)*delta_x));
		   G = G + sqrt(C_D(c,t)/viz);

		   //Gtheta2 is sum of all cells and take average
		   G_theta2 = G_theta2 + sqrt(C_D(c,t)/viz)*delta_x*delta_y/(fabs(C_U(c,t)*delta_y) + fabs(C_V(c,t)*delta_x));
			cellNo = cellNo + 1;
			   
       }
     end_c_loop(c,t)
	G = G/cellNo;
	 //sum of all theta
	 printf("\n theta = %g\n", theta);
	 //sum of all cell Gtheta
	 printf("\n Gtheta2 = %g\n", G_theta2);
	 theta = theta/cellNo;
	 G_theta2 = G_theta2/cellNo;
	 G_theta = G*theta; //average G times average time

     

     }
	
   printf("\n theta = %g\n", theta);
   printf("\n G = %g\n", G);
   printf("\n Gtheta = %g\n", G_theta);
	printf("\n Gtheta2 = %g\n", G_theta2);
   
	// viscosity (ut + u)
	printf("\n Viz = %g\n", viz);
   printf("\n cellNo = %g\n", cellNo);

}