/*
*				example3.c:
*
*		User initialization example for Front Package:
*
*	Copyright 1999 by The University at Stony Brook, All rights reserved.
*
*	This is example of three circles all moving the a normal velocity.
*	Bifurcation occurs when they meet each other. FronTier solves
*	the bifurcation automatically.
*
*	This example shows how to convert the interface into arrays
*	of floats. Functions introduced in this example include:
*/

#include <FronTier.h>

	/*  Function Declarations */

static float multi_circle_func(POINTER,float*);
static void map_output_interface(Front*);

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

typedef struct {
        int num_cir;
        float **cen;
        float *rad;
} TMC_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;
	TMC_PARAMS mc_params;

	/* Initialize basic computational data */

	f_basic.dim = 3;	
	f_basic.L[0] = 0.0;	f_basic.L[1] = 0.0;	f_basic.L[2] = 0.0;
	f_basic.U[0] = 1.0;	f_basic.U[1] = 1.0;	f_basic.U[2] = 1.0;
	f_basic.gmax[0] = 40;	f_basic.gmax[1] = 40;	f_basic.gmax[2] = 40;
	f_basic.boundary[0][0] = f_basic.boundary[0][1] = DIRICHLET_BOUNDARY;
	f_basic.boundary[1][0] = f_basic.boundary[1][1] = DIRICHLET_BOUNDARY;
	f_basic.boundary[2][0] = f_basic.boundary[2][1] = DIRICHLET_BOUNDARY;
	f_basic.size_of_intfc_state = 0;
	

	FrontStartUp(&front,&f_basic);

	/* Start testing */

	front.step = 0;


	/* Initialize interface through level function */

	mc_params.num_cir = 1;
        vector(&mc_params.rad,mc_params.num_cir,FLOAT);
        matrix(&mc_params.cen,mc_params.num_cir,3,FLOAT);
	mc_params.cen[0][0] = 0.4;
	mc_params.cen[0][1] = 0.4;
	mc_params.cen[0][2] = 0.4;
	mc_params.rad[0] = 0.3;

	level_func_pack.neg_component = 1;
	level_func_pack.pos_component = 2;
	level_func_pack.func_params = (POINTER)&mc_params;
	level_func_pack.func = multi_circle_func;
	
	FrontInitIntfc(&front,&level_func_pack);

	map_output_interface(&front);

	clean_up(0);
}

static float multi_circle_func(
        POINTER func_params,
        float *coords)
{
        TMC_PARAMS *mc_params = (TMC_PARAMS*)func_params;
        float **cen = mc_params->cen;
        float *rad = mc_params->rad;
        float num_cir = mc_params->num_cir;
        float dist,dmin;
        int   i,imin;
                                                                                
        dmin = HUGE;
        for (i = 0; i < num_cir; ++i)
        {
            dist = sqrt(sqr(coords[0] - cen[i][0]) +
                        sqr(coords[1] - cen[i][1]) +
			sqr(coords[2] - cen[i][2]));
            if (dist - rad[i] < dmin)
            {
                dmin = dist - rad[i];
                imin = i;
            }
        }
        dist = sqrt(sqr(coords[0] - cen[imin][0]) +
                    sqr(coords[1] - cen[imin][1]) +
		    sqr(coords[2] - cen[imin][2])) - rad[imin];
        return dist;
}       /* end ellipse_func */

static void SaveAsTecplot(
	char *filename, 
	int nP, 
	float *coords, 
	int nTri, 
	int *vIndices)
{
	FILE *hfile = fopen(filename, "w");
	fprintf(hfile, "TITLE = %s \n", filename);
	fprintf(hfile, "VARIABLES = \"X\" \"Y\" \"Z\" \n");
	fprintf(hfile, "ZONE N=%d, E=%d, F=FEPOINT, ET=TETRAHEDRON \n", 
					nP, nTri);
	for(int i=0; i<nP; i++)
	    fprintf(hfile,"%f %f %f \n",coords[i*3+0],coords[i*3+1],
	    				coords[i*3+2]);
	for(int i=0; i<nTri; i++)
	    fprintf(hfile, "%d %d %d %d", vIndices[i*3+0], vIndices[i*3+1],
					vIndices[i*3+2], vIndices[i*3+0]);
	fclose(hfile);
		
}

static void map_output_interface(
	Front *front)
{
	INTERFACE *intfc = front->interf;
        char filename[100];
	int i,j,step;
	step = front->step;
			
	int num_points = NumOfIntfcPoints(intfc);
	int num_tris = NumOfIntfcTris(intfc);
	int dim = 3;
	float *coords = (float*)malloc(num_points*dim*sizeof(float));
	int *vertex_indices = (int*)malloc(num_tris*3*sizeof(int));
	TRI **tris = (TRI**)malloc(num_tris*sizeof(TRI*));
	char *out_name = "output";
	
	// the interface
	printf("Number of interface points = %d\n",num_points);
	printf("Number of interface triangles = %d\n\n",num_tris);

	ArrayOfIntfcTris_FT(intfc,tris);	
	printf("Print first 10 triangles:\n\n");
	for (i = 0; i < 10; i++)
	    print_tri(tris[i],intfc);
	printf("\n\n");

	sprintf(filename,"%s-%d.intfc.plt",out_name,step);
	ArrayOfIntfcTris(intfc,coords,vertex_indices);	
	SaveAsTecplot(filename,num_points,coords,num_tris,vertex_indices);
	
	int num_surfaces = NumOfSurfaces(intfc);
	printf("Number of surfaces = %d\n",num_surfaces);
	SURFACE **surfaces;
	surfaces = (SURFACE**)malloc(num_surfaces*sizeof(SURFACE*));
	
	ArrayOfSurfaces(intfc, surfaces);
	// print the tris of each surface
	for(i = 0; i < num_surfaces; i++)
	{
	    sprintf(filename,"%s-%d.surface-%d.plt",out_name,step,i);
	    num_points = NumOfSurfPoints(surfaces[i]);
	    num_tris   = NumOfSurfTris(surfaces[i]);
	    printf("Number of points on surface %d = %d\n",
	    			i+1,num_points);
	    printf("Number of tris on surface %d = %d\n",
	    			i+1,num_tris);

	    ArrayOfSurfTris_FT(surfaces[i],tris);	
	    printf("Print first 10 triangles:\n\n");
	    for (j = 0; j < 10; j++)
	        print_tri(tris[i],intfc);
	    printf("\n\n");

	    ArrayOfSurfTris(surfaces[i],coords,vertex_indices);
	    SaveAsTecplot(filename,num_points,coords,num_tris,vertex_indices);
	}
}	/* end map_output */