/************************************************************************************
FronTier is a set of libraries that implements differnt types of Front Traking algorithms.
Front Tracking is a numerical method for the solution of partial differential equations 
whose solutions have discontinuities.  


Copyright (C) 1999 by The University at Stony Brook. 
 

This library is free software; you can redistribute it and/or
modify it under the terms of the GNU Lesser General Public
License as published by the Free Software Foundation; either
version 2.1 of the License, or (at your option) any later version.

This library is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
Lesser General Public License for more details.

You should have received a copy of the GNU Lesser General Public
License along with this library; if not, write to the Free Software
Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA

******************************************************************************/


/*
*		ex_mvortex.c
*
*	I made this example by example00.c
*		
*
*
*
*/

#include <FronTier.h>

	/*  Function Declarations */
int     init_io( int,char**);
static void test_propagate(Front*);
static float level_circle_func(POINTER,float*);
static int test_vortex_vel(POINTER,Front*,POINT*,HYPER_SURF_ELEMENT*,
	                       HYPER_SURF*,float*);

char in_name[100],sc_name[100],out_name[100];

/********************************************************************
 *	Level function parameters for the initial interface 	    *
 ********************************************************************/

typedef struct {
	        /* equation for line is x^2/a^2 + y^2/b^2 = 1 */
        float x0;
        float y0;        
        float R; 
} CIRCLE_PARAMS;


/********************************************************************
 *	Velocity function parameters for the front	 	    *
 ********************************************************************/
/* typedef struct {
	float i1,i2;
        float cen1[2],cen2[2];
} DOUBLE_VORTEX_PARAMS;
*/

/*
typedef struct {
	int dim;
	float coeff;
	float time;
	float cos_time;
	char type[10];
	float cen[MAXD];
	float rad;
} VORTEX_PARAMS;
*/

int main(int argc, char **argv)
{
	static Front front;
	static RECT_GRID comp_grid;
	F_BASIC_DATA f_basic;
	LEVEL_FUNC_PACK level_func_pack;
	VELO_FUNC_PACK velo_func_pack;
	CIRCLE_PARAMS circle_params;	/* level function parameters */
	VORTEX_PARAMS vortex_params; /* velocity function parameters */
	Locstate  sl;

	init_io(argc,argv);

	/* Initialize basic computational data */

	f_basic.dim = 2;	
	f_basic.L[0] = 0.0;	f_basic.L[1] = 0.0;
	f_basic.U[0] = 1.0;	f_basic.U[1] = 1.0;
	f_basic.gmax[0] = 100;	f_basic.gmax[1] = 100;		//myex grid size
	f_basic.boundary[0][0] = f_basic.boundary[0][1] = PERIODIC_BOUNDARY;
	f_basic.boundary[1][0] = f_basic.boundary[1][1] = PERIODIC_BOUNDARY;
	f_basic.size_of_intfc_state = 0;

	FrontStartUp(&front,&f_basic);

	/* Start testing */

	front.step = 0;


	/* Initialize interface through level function */

	//circle_params.x0 = 0.5;
	//circle_params.y0 = 0.25;
	//circle_params.R = 0.15;

	level_func_pack.neg_component = 1;
	level_func_pack.pos_component = 2;
	//level_func_pack.func_params = (POINTER)&circle_params;
	//level_func_pack.func = level_circle_func;

	level_func_pack.func_params = NULL;
	level_func_pack.func = NULL;

	level_func_pack.num_points = 250;		//myex num points
	level_func_pack.is_closed_curve = YES;

	matrix(&level_func_pack.point_array,level_func_pack.num_points,2,sizeof(float));
	int i;
	for (i = 0; i < level_func_pack.num_points; ++i)
	{
	    float phi = i*2.0*PI/(float)level_func_pack.num_points;
	    level_func_pack.point_array[i][0] = 0.5 + 0.15*cos(phi);
	    level_func_pack.point_array[i][1] = 0.25 + 0.15*sin(phi);
	}

	FrontInitIntfc(&front,&level_func_pack);

	/* Initialize velocity field function */

	vortex_params.dim=2;
	vortex_params.type[0]='M';
	vortex_params.coeff=1.0;
	vortex_params.time=0;
	vortex_params.cos_time=1.5;
	vortex_params.cen[0]=0.5;
	vortex_params.cen[1]=0.25;
	vortex_params.rad=0.15;

	velo_func_pack.func_params = (POINTER)&vortex_params;
	velo_func_pack.func = test_vortex_vel;

	FrontInitVelo(&front,&velo_func_pack);

	/* Propagate the front */

	test_propagate(&front);

	clean_up(0);
	return 0;
}


/*ARGSUSED*/
	int	init_io(
	int		argc,
	char		**argv)
{
	argc--;
	argv++;
	strcpy(out_name,"intfc");
	while (argc >= 1)
	{
	    if (argv[0][0] != '-')
	    {
		printf("Usage: example -o output\n");
		exit(1);
	    }
	    switch(argv[0][1]) {
	    case 'o':
		zero_scalar(out_name,100);
		strcpy(out_name,argv[1]);
		freopen(out_name,"w",stdout);
		argc -= 2;
		argv += 2;
		break;
	    }
	}
}

#define is_output_time(front,i,interval)                        \
        ((front)->time+front->dt >= (i)*(interval)) 

static  void test_propagate(
        Front *front)
{
        int ip,im,status,count;
        Front *newfront;
        float max_time,dt,dt_frac,CFL;
        float print_time_interval,movie_frame_interval;
        int max_step;
        bool is_print_time, is_movie_time, time_limit_reached;
        char s[10];
        float fcrds[MAXD];
        int  dim = front->rect_grid->dim;

	max_time = 1.5;
	max_step = 400;
	print_time_interval = 1.0;
	movie_frame_interval = 0.02;

        redistribute(front,YES,NO);

        printf("dim = %d\n", dim);
        front->time = 0.0;

        print_front_output(front,out_name);
        show_front_output(front,out_name);

        front->dt = 0.0;
        status = FrontAdvance(front->dt,&dt_frac,front,&newfront,
                                (POINTER)NULL);
        ip = im = 1;
        CFL = Time_step_factor(front);
	printf("CFL = %f\n",CFL);

	Frequency_of_redistribution(front,GENERAL_WAVE)=2; // myex: Frequency of redistribution

	printf("Frequency_of_redistribution(front,GENERAL_WAVE) = %d\n",
		Frequency_of_redistribution(front,GENERAL_WAVE));
        time_limit_reached = NO;

        for (;;)
        {
            front->dt = CFL*(*front->max_front_time_step)(front,fcrds);

            is_print_time = is_movie_time = NO;
            if (is_output_time(front,im,movie_frame_interval))
            {
                front->dt = im*movie_frame_interval - front->time;
                is_movie_time = YES;
            }
            if (is_output_time(front,ip,print_time_interval))
            {
                front->dt = ip*print_time_interval - front->time;
                is_print_time = YES;
            }
            if (front->time+front->dt > max_time)
            {
                    front->dt = max_time - front->time;
                    time_limit_reached = YES;
            }
            status = FrontAdvance(front->dt,&dt_frac,front,&newfront,
                                (POINTER)NULL);
            count = 0;
            while (status == MODIFY_TIME_STEP)
            {
                    front->dt *= CFL;
                    status = FrontAdvance(front->dt,&dt_frac,front,&newfront,
                                (POINTER)NULL);
                    count++;
                    if (count > 10) clean_up(ERROR);
            }
            assign_interface_and_free_front(front,newfront);

            ++front->step;
            front->time += front->dt;

            printf("\ntime = %f   step = %5d   dt = %f\n",
                        front->time,front->step,front->dt);
            if (is_print_time ||
                time_limit_reached || front->step >= max_step)
            {
                    print_front_output(front,out_name);
                    ++ip;
            }
            if (is_movie_time ||
                time_limit_reached || front->step >= max_step)
            {
                    show_front_output(front,out_name);
                    ++im;
            }

            fflush(stdout);
            if (time_limit_reached || front->step >= max_step)
                    break;
        }
        (void) delete_interface(front->interf);
}       /* end test_propagate */

/********************************************************************
 *	Sample (circle) level function for the initial interface    *
 ********************************************************************/

static float level_circle_func(
        POINTER func_params,
        float *coords)
{
	CIRCLE_PARAMS *circle_params = (CIRCLE_PARAMS*)func_params;
	float x0,y0,R,dist;

	x0 = circle_params->x0;
	y0 = circle_params->y0;
	R  = circle_params->R;

	dist = sqrt(sqr(coords[0] - x0) + sqr(coords[1] - y0)) - R;
	return dist;
}	/* end level_circle_func */


/********************************************************************
 *	Sample (circle) velocity function for the front    *
 ********************************************************************/

/*
static int test_vortex_vel(
	POINTER params,
	Front *front,
	POINT *p,
	HYPER_SURF_ELEMENT *hse,
	HYPER_SURF *hs,
	float *vel)
{
	DOUBLE_VORTEX_PARAMS *dv_params = (DOUBLE_VORTEX_PARAMS*)params;
	float *coords = Coords(p);
	float d1,d2;
	float s1,s2;
	float *cen1 = dv_params->cen1;
	float *cen2 = dv_params->cen2;
	float dx1,dy1;
	float dx2,dy2;

	dx1 = coords[0] - cen1[0]; 
	dy1 = coords[1] - cen1[1];
	dx2 = coords[0] - cen2[0]; 
	dy2 = coords[1] - cen2[1];

	d1 = sqrt(sqr(dx1) + sqr(dy1));
	d2 = sqrt(sqr(dx2) + sqr(dy2));

	s1 = dv_params->i1/2.0/PI/d1;
	s2 = dv_params->i2/2.0/PI/d2;

	vel[0] =  s1*dy1/d1 + s2*dy2/d2;
	vel[1] = -s1*dx1/d1 - s2*dx2/d2;
}
*/

static int test_vortex_vel(
        POINTER params,
        Front *front,
        POINT *p,
        HYPER_SURF_ELEMENT *hse,
        HYPER_SURF *hs,
        float *vel)
{
        VORTEX_PARAMS *vortex_params;
        int i, dim;
        float coeff,coeff2,xtemp,ytemp;
        char *type;
        float *coords = Coords(p);
        vortex_params = (VORTEX_PARAMS*)params;
	float x,y,z;

	dim = vortex_params->dim;
        type = vortex_params->type;

        coeff2 = 1.0;
        if(vortex_params->time > 0 && front->time >= vortex_params->time)
            coeff2 = -1.0;
        if(vortex_params->cos_time > 0 )
            coeff2 = cos((front->time*PI)/vortex_params->cos_time);
	x = coords[0];
	y = coords[1];
	if (dim == 3)
	    z = coords[2];

	if (dim == 2)
	{
	    switch(type[0])
	    {
	        case 'm':
		case 'M':
		    //four vortex work
		    xtemp = 4*PI*(x+0.5);
		    ytemp = 4*PI*(y+0.5);
		    vel[0] = coeff2*sin(xtemp)*sin(ytemp); 
		    vel[1] = coeff2*cos(xtemp)*cos(ytemp); 
		    //end four vortex work
		    break;
		case 's':
		case 'S':
			//2D single vortex motion
		    vel[0] = -coeff2*sin(PI*x)*sin(PI*x)*sin(2*PI*y);
		    vel[1] = coeff2*sin(2*PI*x)*sin(PI*y)*sin(PI*y);
		    break;
		default:
		    screen("Undefined vortex type!\n");
		    clean_up(ERROR);
	    }
	}
	else if (dim == 3)
	{
	    switch(type[0])
	    {
	        case 'm':
	        case 'M':
		    //double vortex work
	            vel[0] = coeff2*2*sin(PI*x)*sin(PI*x)*sin(2*PI*y)*sin(2*PI*z);
	            vel[1] = -coeff2*sin(2*PI*x)*sin(PI*y)*sin(PI*y)*sin(2*PI*z);
	            vel[2] = -coeff2*sin(2*PI*x)*sin(2*PI*y)*sin(PI*z)*sin(PI*z);
		    break;
	        case 's':
	        case 'S':
			//shearing flow motion
		    vel[0] = coeff2*sin(PI*x)*sin(PI*x)*sin(2*PI*y);
		    vel[1] = -coeff2*sin(2*PI*x)*sin(PI*y)*sin(PI*y);
		    vel[2] = coeff2*sqr(1-2*sqrt((x-0.5)*(x-0.5)+(y-0.5)*(y-0.5)));
		    /*
			 * 3D single vortex motion
		    vel[0] = -coeff2*sin(PI*x)*sin(PI*x)*sin(2*PI*y);
		    vel[1] = coeff2*sin(2*PI*x)*sin(PI*y)*sin(PI*y);
		    vel[2] = 0;*/
		    break;
	        default:
	            screen("Undefined time dependency type!\n");
	            clean_up(ERROR);
	    }
	}    
}       /* end vortex_vel */