//
//  define_init.c
//  initial_field_udf
//
//  This code reads in initial field information from the data file IC1.dat which contains information of the coordinates of the grids as well as the velocity of the flow in the corresponding grid. Then the code initializes the velocity by sweeping through the cells and defining the velocity at each cell according to the velocity values defined by the data file.
//
//  Created by Chiyu Max Jiang on 8/17/14.
//  Copyright (c) 2014 Chiyu Max Jiang. All rights reserved.
//

#include "udf.h"

DEFINE_INIT(define_init_velocity,d)
{
    int i, j, k;
	double a[36864][6]; //construct a matrix to store the data from the file
    FILE *fp = fopen("IC.dat", "r");
    
    cell_t c;
    Thread *t;
    real xc[ND_ND];
    real n1, n2, N1;
    int N2;

    /*scan the file and read the data into the matrix*/
    for (i = 0; i < 36864; i++) {
        for (j = 0; j < 6; j++) {
            fscanf(fp,"%lf",&a[i][j]);
        }
    }
	fclose(fp);
    
    /* loop over all cell threads in the domain */
    thread_loop_c(t,d)
    {
        n1 = 0;
        /* loop over all cells */
        begin_c_loop_all(c,t)
        {
			n1++;
            C_CENTROID(xc,c,t);
            if (n1-384*(ceil(n1/384)-1)<=192) {
            N1 = 192*(ceil(n1/384)-1)+n1-384*(ceil(n1/384)-1);
            }
            else {
            N1 = 192*(ceil(n1/384)-1)+n1-384*(ceil(n1/384)-1)-192;
            N1 = 192*(ceil(N1/192)-1)+193-(N1-192*(ceil(N1/192)-1));
            }
            n2 = (N1-192*(ceil(N1/192)-1)-1)*192+193-ceil(N1/192);
            N2 = n2;
            /*  define u v and w velocity components */
            C_U(c,t) = a[N2-1][3];
            C_V(c,t) = a[N2-1][4];
        }
        end_c_loop_all(c,t);
    }
    
}