/************************************************************************************
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
******************************************************************************/
/*
* test_map_3d.c:
*
* User initialization example for Front Package:
*
* Copyright 1999 by The University at Stony Brook, All rights reserved.
*
* This example shows how to use an array of points to make a curve.
*
*/
#include
/* Function Declarations */
static int init_io( int, char**);
static void test_propagate(Front*);
static int test_vortex_vel(POINTER,Front*,POINT*,HYPER_SURF_ELEMENT*,
HYPER_SURF*,float*);
char in_name[100],sc_name[100],out_name[100];
/********************************************************************
* Velocity function parameters for the front *
********************************************************************/
typedef struct {
float i1,i2;
float cen1[2],cen2[2];
} DOUBLE_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;
DOUBLE_VORTEX_PARAMS dv_params; /* velocity function parameters */
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; //ex15: grid size
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.size_of_intfc_state = 0;
FrontStartUp(&front,&f_basic);
/* Start testing */
front.step = 0;
/* Initialize interface through level function */
level_func_pack.neg_component = 1;
level_func_pack.pos_component = 2;
level_func_pack.func_params = NULL;
level_func_pack.func = NULL;
level_func_pack.num_points = 250; //ex15: num points
level_func_pack.is_closed_curve = YES;
// ---original--- matrix(&level_func_pack.point_array,500,2,sizeof(float));
matrix(&level_func_pack.point_array,level_func_pack.num_points,2,sizeof(float));
int i;
// ---original--- for (i = 0; i < 500; ++i)
for (i = 0; i < level_func_pack.num_points; ++i)
{
// --- original --- float phi = i*2.0*PI/500.0;
float phi = i*2.0*PI/(float)level_func_pack.num_points;
level_func_pack.point_array[i][0] = 0.5 + 0.3*cos(phi);
level_func_pack.point_array[i][1] = 0.5 + 0.3*sin(phi);
}
FrontInitIntfc(&front,&level_func_pack);
/* Initialize velocity field function */
dv_params.cen1[0] = 0.25;
dv_params.cen1[1] = 0.50;
dv_params.cen2[0] = 0.75;
dv_params.cen2[1] = 0.50;
dv_params.i1 = -0.5;
dv_params.i2 = 0.5;
velo_func_pack.func_params = (POINTER)&dv_params;
velo_func_pack.func = test_vortex_vel;
FrontInitVelo(&front,&velo_func_pack);
/* Propagate the front */
test_propagate(&front);
clean_up(0);
clean_up(0);
}
static 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 = 2.0;
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; // ex15: 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) 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;
} /* end test_vortex_vel */