/***************************************************************************************************
 THIS CODE ILLUSTRATES HOW TO CONVERT AN ANALYTICAL SOLUTION FROM POLAR COORDINATES TO
 CARTESIAN COORDINATES AND SUBSEQUENTLY BE USED TO INITIALIZE A FLUENT SIMULATION
 THIS ONLY WORKS FOR 2D SIMULATIONS.
 Coded by Tony Saad, © 2008
 downloaded from: http://pleasemakeanote.blogspot.com/
***************************************************************************************************/
#include "udf.h"
#define Pi 3.141592653

DEFINE_INIT(Initialize_Velocity_Field, theDomain)
{
 real xCentroid=0.0; /* abscissa of the cell centroid */
 real yCentroid = 0.0; /* ordinate of the cell centroid */
 real rCentroidMagnitude = 0.0; /* magnitude of the position vector */
 real theta = 0.0; /* angle made by the position vector and the x axis */
 real uR = 0.0; /* radial velocity (given by the analytical solution) */
 real uRx =0.0; /* x projection of the radial velocity */
 real uRy = 0.0; /* y projection of the radial velocity */
 real uTheta = 0.0; /* tangential velocity (given by the analytical solution) */
 real uThetax = 0.0; /* x projection of the tangential velocity */
 real uThetay = 0.0; /* y projection of the tangential velocity */
 /*declare thread variables*/
 cell_t theCell;
 real cellCentroidCoordinates[ND_ND];
 real origin[ND_ND];
 Thread* theThread;

 /*-------------------------------------------*/
 /* SET THE LOCATION OF THE ORIGIN */
 origin[0] = 0.0;
 origin[1] = 0.0;
 /*-------------------------------------------*/

 /* Loop over all cell threads in the domain*/
 thread_loop_c(theThread,theDomain) {
 /* loop over all cells in the domain thread */
 begin_c_loop_all (theCell,theThread) {

 /* get the centroid of the current cell */
 C_CENTROID(cellCentroidCoordinates,theCell,theThread);
 /* take into account the location of the origin */
 cellCentroidCoordinates[0] -= origin[0];
 cellCentroidCoordinates[1] -= origin[1];

 /* store the centroid location for ease of access */
 xCentroid = cellCentroidCoordinates[0];
 yCentroid = cellCentroidCoordinates[1];

 /* compute the radius of the current cell */
 rCentroidMagnitude = NV_MAG(cellCentroidCoordinates); 

 /*-------------------------------------------*/
 /* DEFINE THE RADIAL AND AXIAL VELOCITIES */
 /* ENTER YOUR ANALYTICAL SOLUTION HERE */
 uR = - rCentroidMagnitude;
 uTheta = 1/(rCentroidMagnitude + 1/rCentroidMagnitude);
 /*-------------------------------------------*/

 /* compute the angle made by the position vector and the x-axis */
 theta = atan(fabs(yCentroid/xCentroid));
 if (xCentroid<=0 && yCentroid>=0) {
 theta = Pi - theta;
 } else if (xCentroid<=0 && yCentroid<=0) {
 theta += Pi;
 } else if (xCentroid>0 && yCentroid<0) {
 theta =2*Pi - theta;
 }

 /* compute the x and y components of the radial velocity */
 uRx = uR*cos(theta);
 uRy = uR*sin(theta);

 /* compute the x and y components of the tangential velocity */
 uThetax = uTheta*cos(theta + Pi/2);
 uThetay = uTheta*sin(theta + Pi/2);

 /*set the x and y velocities in fluent*/
 C_U(theCell,theThread) = uRx + uThetax;
 C_V(theCell,theThread) = uRy + uThetay;
 }
 end_c_loop_all (theCell,theThread)
 }
}