806 lines
27 KiB
C
806 lines
27 KiB
C
/* routines dealing with hashgrid for lennardjones.c */
|
|
/* Note: a possible improvement would be to use pointers instead of integers */
|
|
/* when referencing other hashgrid cells or particles */
|
|
|
|
int bc_grouped(int bc)
|
|
/* regroup boundary conditions by type: rectangular = 0, periodic = 1, Klein = 2, Boy = 3, L shape = 4, ... */
|
|
{
|
|
switch (bc) {
|
|
case (BC_SCREEN): return(0);
|
|
case (BC_RECTANGLE): return(0);
|
|
case (BC_CIRCLE): return(0);
|
|
case (BC_PERIODIC): return(1);
|
|
case (BC_PERIODIC_CIRCLE): return(1);
|
|
case (BC_EHRENFEST): return(0);
|
|
case (BC_PERIODIC_FUNNEL): return(1);
|
|
case (BC_RECTANGLE_LID): return(0);
|
|
case (BC_RECTANGLE_WALL): return(0);
|
|
case (BC_PERIODIC_TRIANGLE): return(1);
|
|
case (BC_KLEIN): return(2);
|
|
case (BC_SCREEN_BINS): return(0);
|
|
case (BC_BOY): return(3);
|
|
case (BC_GENUS_TWO): return(4);
|
|
case (BC_ABSORBING): return(0);
|
|
case (BC_REFLECT_ABS): return(0);
|
|
case (BC_REFLECT_ABS_BOTTOM): return(0);
|
|
default:
|
|
{
|
|
printf("Warning: Hashgrid will not be properly initialised, update bc_grouped()\n\n");
|
|
return(-1);
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
int mhash(int i, int j)
|
|
{
|
|
return(i*HASHY + j);
|
|
}
|
|
|
|
|
|
int hash_cell(double x, double y)
|
|
/* compute hash grid position of particle at (x,y) */
|
|
/* returns number of hash cell */
|
|
{
|
|
static int first = 1;
|
|
static double lx, ly, padding;
|
|
int i, j;
|
|
|
|
if (first)
|
|
{
|
|
if (!NO_WRAP_BC) padding = 0.0;
|
|
else padding = HASHGRID_PADDING;
|
|
lx = BCXMAX - BCXMIN + 2.0*padding;
|
|
ly = BCYMAX - (BCYMIN) + 2.0*padding;
|
|
first = 0;
|
|
}
|
|
|
|
if (CENTER_VIEW_ON_OBSTACLE) x -= xshift;
|
|
|
|
i = (int)((double)HASHX*(x - BCXMIN + padding)/lx);
|
|
j = (int)((double)HASHY*(y - BCYMIN + padding)/ly);
|
|
|
|
if (i<0) i = 0;
|
|
else if (i>=HASHX) i = HASHX-1;
|
|
if (j<0) j = 0;
|
|
else if (j>=HASHY) j = HASHY-1;
|
|
|
|
return(mhash(i,j));
|
|
// printf("Mapped (%.3f,%.3f) to (%i, %i)\n", x, y, ij[0], ij[1]);
|
|
}
|
|
|
|
|
|
void init_hashgrid(t_hashgrid hashgrid[HASHX*HASHY])
|
|
/* initialise table of neighbouring cells for each hashgrid cell, depending on boundary condition */
|
|
{
|
|
int i, j, k, p, q, m, i1, j1;
|
|
|
|
printf("Initializing hash grid\n");
|
|
|
|
/* bulk of the table */
|
|
for (i=0; i<HASHX-1; i++)
|
|
for (j=0; j<HASHY-1; j++)
|
|
{
|
|
m = mhash(i, j);
|
|
hashgrid[m].nneighb = 9;
|
|
|
|
for (p=-1; p<2; p++)
|
|
for (q=-1; q<2; q++)
|
|
hashgrid[m].neighbour[3*p+q+4] = mhash(i+p, j+q);
|
|
}
|
|
|
|
/* different boundary conditions */
|
|
switch (bc_grouped(BOUNDARY_COND)) {
|
|
case (0): /* rectangular b.c. */
|
|
{
|
|
/* left/right boundaries */
|
|
for (j=1; j<HASHY-1; j++)
|
|
{
|
|
for (i=0; i<HASHX; i+=HASHX-1) hashgrid[mhash(i, j)].nneighb = 6;
|
|
|
|
for (q=-1; q<2; q++)
|
|
{
|
|
for (p=0; p<2; p++) hashgrid[mhash(0, j)].neighbour[2*(q+1)+p] = mhash(p, j+q);
|
|
for (p=-1; p<1; p++) hashgrid[mhash(HASHX-1, j)].neighbour[2*(q+1)+p+1] = mhash(HASHX-1+p, j+q);
|
|
}
|
|
}
|
|
|
|
/* top/down boundaries */
|
|
for (i=1; i<HASHX-1; i++)
|
|
{
|
|
for (j=0; j<HASHY; j+=HASHY-1) hashgrid[mhash(i, j)].nneighb = 6;
|
|
|
|
for (p=-1; p<2; p++)
|
|
{
|
|
for (q=0; q<2; q++) hashgrid[mhash(i, 0)].neighbour[2*(p+1)+q] = mhash(i+p, q);
|
|
for (q=-1; q<1; q++) hashgrid[mhash(i, HASHY-1)].neighbour[2*(p+1)+q+1] = mhash(i+p, HASHY-1+q);
|
|
}
|
|
}
|
|
|
|
/* corners */
|
|
for (i=0; i<HASHX; i+=HASHX-1) for (j=0; j<HASHY; j+=HASHY-1) hashgrid[mhash(i, j)].nneighb = 4;
|
|
for (p=0; p<2; p++) for (q=0; q<2; q++)
|
|
{
|
|
hashgrid[mhash(0,0)].neighbour[2*p+q] = mhash(p, q);
|
|
hashgrid[mhash(HASHX-1,0)].neighbour[2*p+q] = mhash(HASHX-1-p, q);
|
|
hashgrid[mhash(0,HASHY-1)].neighbour[2*p+q] = mhash(p, HASHY-1-q);
|
|
hashgrid[mhash(HASHX-1,HASHY-1)].neighbour[2*p+q] = mhash(HASHX-1-p, HASHY-1-q);
|
|
}
|
|
break;
|
|
}
|
|
|
|
case(1): /* periodic b.c. */
|
|
{
|
|
/* left/right boundaries */
|
|
for (j=0; j<HASHY; j++) for (i=0; i<HASHX; i+=HASHX-1)
|
|
{
|
|
hashgrid[mhash(i, j)].nneighb = 9;
|
|
|
|
for (p=-1; p<2; p++) for (q=-1; q<2; q++)
|
|
{
|
|
i1 = (i+p+HASHX)%HASHX;
|
|
j1 = (j+q+HASHY)%HASHY;
|
|
hashgrid[mhash(i,j)].neighbour[3*(p+1)+q+1] = mhash(i1, j1);
|
|
}
|
|
}
|
|
|
|
/* top/bottom boundaries */
|
|
for (i=1; i<HASHX-1; i++) for (j=0; j<HASHY; j+=HASHY-1)
|
|
{
|
|
hashgrid[mhash(i, j)].nneighb = 9;
|
|
|
|
for (p=-1; p<2; p++) for (q=-1; q<2; q++)
|
|
{
|
|
i1 = (i+p+HASHX)%HASHX;
|
|
j1 = (j+q+HASHY)%HASHY;
|
|
hashgrid[mhash(i,j)].neighbour[3*(p+1)+q+1] = mhash(i1, j1);
|
|
}
|
|
}
|
|
break;
|
|
}
|
|
|
|
case(2): /* Klein bottle b.c. */
|
|
{
|
|
/* left/right boundaries */
|
|
for (j=0; j<HASHY; j++) for (i=0; i<HASHX; i+=HASHX-1)
|
|
{
|
|
hashgrid[mhash(i, j)].nneighb = 9;
|
|
|
|
for (p=-1; p<2; p++) for (q=-1; q<2; q++)
|
|
{
|
|
i1 = (i+p+HASHX)%HASHX;
|
|
if (((i == 0)&&(p >= 0))||((i == HASHX-1)&&(p <= 0))) j1 = (j+q+HASHY)%HASHY;
|
|
else j1 = (2*HASHY-1-j-q)%HASHY;
|
|
hashgrid[mhash(i,j)].neighbour[3*(p+1)+q+1] = mhash(i1, j1);
|
|
}
|
|
}
|
|
|
|
/* top/bottom boundaries */
|
|
for (i=1; i<HASHX-1; i++) for (j=0; j<HASHY; j+=HASHY-1)
|
|
{
|
|
hashgrid[mhash(i, j)].nneighb = 9;
|
|
|
|
for (p=-1; p<2; p++) for (q=-1; q<2; q++)
|
|
{
|
|
i1 = (i+p+HASHX)%HASHX;
|
|
j1 = (j+q+HASHY)%HASHY;
|
|
hashgrid[mhash(i,j)].neighbour[3*(p+1)+q+1] = mhash(i1, j1);
|
|
}
|
|
}
|
|
break;
|
|
}
|
|
|
|
case(3): /* Boy surface/projective plane b.c. */
|
|
{
|
|
/* left/right boundaries */
|
|
for (j=1; j<HASHY-1; j++) for (i=0; i<HASHX; i+=HASHX-1)
|
|
{
|
|
hashgrid[mhash(i, j)].nneighb = 9;
|
|
|
|
for (p=-1; p<2; p++) for (q=-1; q<2; q++)
|
|
{
|
|
i1 = (i+p+HASHX)%HASHX;
|
|
if (((i == 0)&&(p >= 0))||((i == HASHX-1)&&(p <= 0))) j1 = (j+q+HASHY)%HASHY;
|
|
else j1 = (2*HASHY-1-j-q)%HASHY;
|
|
hashgrid[mhash(i,j)].neighbour[3*(p+1)+q+1] = mhash(i1, j1);
|
|
}
|
|
}
|
|
|
|
/* top/bottom boundaries */
|
|
for (i=1; i<HASHX-1; i++) for (j=0; j<HASHY; j+=HASHY-1)
|
|
{
|
|
hashgrid[mhash(i, j)].nneighb = 9;
|
|
|
|
for (p=-1; p<2; p++) for (q=-1; q<2; q++)
|
|
{
|
|
if (((j == 0)&&(q >= 0))||((j == HASHY-1)&&(q <= 0))) i1 = (i+p+HASHX)%HASHX;
|
|
else i1 = (2*HASHX-1-i-p)%HASHX;
|
|
j1 = (j+q+HASHY)%HASHY;
|
|
hashgrid[mhash(i,j)].neighbour[3*(p+1)+q+1] = mhash(i1, j1);
|
|
}
|
|
}
|
|
|
|
/* corners */
|
|
for (i=0; i<HASHX; i+=HASHX-1) for (j=0; j<HASHY; j+=HASHY-1) hashgrid[mhash(i, j)].nneighb = 7;
|
|
|
|
for (p=0; p<2; p++) for (q=0; q<2; q++) hashgrid[mhash(0,0)].neighbour[2*p+q] = mhash(p,q);
|
|
hashgrid[mhash(0,0)].neighbour[4] = mhash(HASHX-1,HASHY-1);
|
|
hashgrid[mhash(0,0)].neighbour[5] = mhash(HASHX-2,HASHY-1);
|
|
hashgrid[mhash(0,0)].neighbour[6] = mhash(HASHX-1,HASHY-2);
|
|
|
|
for (p=0; p<2; p++) for (q=0; q<2; q++) hashgrid[mhash(HASHX-1,0)].neighbour[2*p+q] = mhash(HASHX-1-p,q);
|
|
hashgrid[mhash(HASHX-1,0)].neighbour[4] = mhash(0,HASHY-1);
|
|
hashgrid[mhash(HASHX-1,0)].neighbour[5] = mhash(1,HASHY-1);
|
|
hashgrid[mhash(HASHX-1,0)].neighbour[6] = mhash(0,HASHY-2);
|
|
|
|
for (p=0; p<2; p++) for (q=0; q<2; q++) hashgrid[mhash(0,HASHY-1)].neighbour[2*p+q] = mhash(p,HASHY-1-q);
|
|
hashgrid[mhash(0,HASHY-1)].neighbour[4] = mhash(HASHX-1,0);
|
|
hashgrid[mhash(0,HASHY-1)].neighbour[5] = mhash(HASHX-2,1);
|
|
hashgrid[mhash(0,HASHY-1)].neighbour[6] = mhash(HASHX-1,2);
|
|
|
|
for (p=0; p<2; p++) for (q=0; q<2; q++) hashgrid[mhash(HASHX-1,HASHY-1)].neighbour[2*p+q] = mhash(HASHX-1-p,HASHY-1-q);
|
|
hashgrid[mhash(HASHX-1,HASHY-1)].neighbour[4] = mhash(0,1);
|
|
hashgrid[mhash(HASHX-1,HASHY-1)].neighbour[5] = mhash(1,1);
|
|
hashgrid[mhash(HASHX-1,HASHY-1)].neighbour[6] = mhash(0,2);
|
|
break;
|
|
}
|
|
|
|
case(4): /* genus 2 (L-shape) b.c. */
|
|
{
|
|
if ((HASHX%2 == 1)||(HASHY%2 == 1))
|
|
{
|
|
printf("Error: HASHX and HASHY must be even for this boundary condition\n");
|
|
exit(0);
|
|
}
|
|
|
|
/* dummy value for unused cells */
|
|
for (i=HASHX/2; i<HASHX; i++) for (j=HASHY/2; j<HASHY; j++) hashgrid[mhash(i, j)].nneighb = 0;
|
|
|
|
/* left boundary */
|
|
for (j=0; j<HASHY; j++)
|
|
{
|
|
hashgrid[mhash(0, j)].nneighb = 9;
|
|
|
|
for (p=-1; p<2; p++) for (q=-1; q<2; q++)
|
|
{
|
|
j1 = (j+q+HASHY)%HASHY;
|
|
if (j1 < HASHY/2) i1 = (p+HASHX)%HASHX;
|
|
else i1 = (p+HASHX/2)%(HASHX/2);
|
|
hashgrid[mhash(0, j)].neighbour[3*(p+1)+q+1] = mhash(i1, j1);
|
|
}
|
|
}
|
|
|
|
/* right boundary */
|
|
for (j=0; j<HASHY; j++)
|
|
{
|
|
if (j < HASHY/2)
|
|
{
|
|
hashgrid[mhash(HASHX-1, j)].nneighb = 9;
|
|
|
|
for (p=-1; p<2; p++) for (q=-1; q<2; q++)
|
|
{
|
|
j1 = (j+q+HASHY)%HASHY;
|
|
if (j1 < HASHY/2) i1 = (HASHX-1+p)%HASHX;
|
|
else i1 = (HASHX/2-1+p)%(HASHX/2);
|
|
hashgrid[mhash(HASHX-1,j)].neighbour[3*(p+1)+q+1] = mhash(i1, j1);
|
|
}
|
|
}
|
|
else
|
|
{
|
|
hashgrid[mhash(HASHX/2-1, j)].nneighb = 9;
|
|
// hashgrid[mhash(HASHX-1, j)].nneighb = 0; /* dummy value for unused boundary */
|
|
|
|
for (p=-1; p<2; p++) for (q=-1; q<2; q++)
|
|
{
|
|
j1 = (j+q+HASHY)%HASHY;
|
|
if (j1 < HASHY/2) i1 = (HASHX-1+p)%HASHX;
|
|
else i1 = (HASHX/2-1+p)%(HASHX/2);
|
|
hashgrid[mhash(HASHX/2-1, j)].neighbour[3*(p+1)+q+1] = mhash(i1, j1);
|
|
}
|
|
}
|
|
|
|
}
|
|
|
|
/* bottom boundary */
|
|
for (i=1; i<HASHX-1; i++)
|
|
{
|
|
hashgrid[mhash(i, 0)].nneighb = 9;
|
|
|
|
for (p=-1; p<2; p++) for (q=-1; q<2; q++)
|
|
{
|
|
i1 = (i+p+HASHX)%HASHX;
|
|
if (i1 < HASHX/2) j1 = (q+HASHY)%HASHY;
|
|
else j1 = (q+HASHY/2)%(HASHY/2);
|
|
hashgrid[mhash(i,0)].neighbour[3*(p+1)+q+1] = mhash(i1, j1);
|
|
}
|
|
}
|
|
|
|
/* top boundary */
|
|
for (i=1; i<HASHX-1; i++)
|
|
{
|
|
if (i < HASHX/2)
|
|
{
|
|
hashgrid[mhash(i, HASHY-1)].nneighb = 9;
|
|
|
|
for (p=-1; p<2; p++) for (q=-1; q<2; q++)
|
|
{
|
|
i1 = (i+p+HASHX/2)%(HASHX/2);
|
|
if (i1 < HASHX/2) j1 = (HASHY-1+q)%HASHY;
|
|
else j1 = (HASHY/2-1+q)%(HASHY/2);
|
|
hashgrid[mhash(i, HASHY-1)].neighbour[3*(p+1)+q+1] = mhash(i1, j1);
|
|
}
|
|
}
|
|
else
|
|
{
|
|
hashgrid[mhash(i, HASHY/2-1)].nneighb = 9;
|
|
// hashgrid[mhash(i, HASHY-1)].nneighb = 0; /* dummy value for unused boundary */
|
|
|
|
for (p=-1; p<2; p++) for (q=-1; q<2; q++)
|
|
{
|
|
i1 = (i+p+HASHX)%HASHX;
|
|
if (i1 < HASHX/2) j1 = (HASHY-1+q)%HASHY;
|
|
else j1 = (HASHY/2-1+q)%(HASHY/2);
|
|
hashgrid[mhash(i, HASHY/2-1)].neighbour[3*(p+1)+q+1] = mhash(i1, j1);
|
|
}
|
|
}
|
|
}
|
|
|
|
/* TO DO : add more cells for "corners" ? */
|
|
break;
|
|
}
|
|
|
|
default: /* do nothing */;
|
|
}
|
|
|
|
// for (i=0; i<HASHX; i++)
|
|
// {
|
|
// for (j=0; j<HASHY; j++)
|
|
// {
|
|
// for (k=0; k<hashgrid[mhash(i,j)].nneighb; k++)
|
|
// {
|
|
// m = hashgrid[mhash(i,j)].neighbour[k];
|
|
// p = m/HASHY;
|
|
// q = m%HASHY;
|
|
// if (vabs((double)(p-i)) + vabs((double)(q-j)) > 2.0)
|
|
// printf("Grid cell (%i, %i) - neighbour %i = %i = (%i, %i)\n", i, j, k, m, p, q);
|
|
// }
|
|
// // sleep(1);
|
|
// }
|
|
// // sleep(1);
|
|
// }
|
|
sleep(1);
|
|
}
|
|
|
|
void update_hashgrid(t_particle* particle, t_hashgrid* hashgrid, int verbose)
|
|
{
|
|
int i, j, k, n, m, max = 0, hashcell;
|
|
|
|
// printf("Updating hashgrid_number\n");
|
|
for (i=0; i<HASHX*HASHY; i++) hashgrid[i].number = 0;
|
|
// printf("Updated hashgrid_number\n");
|
|
|
|
/* place each particle in hash grid */
|
|
for (k=0; k<ncircles; k++)
|
|
if (particle[k].active)
|
|
{
|
|
// printf("placing circle %i\t", k);
|
|
hashcell = hash_cell(particle[k].xc, particle[k].yc);
|
|
n = hashgrid[hashcell].number;
|
|
if (n < HASHMAX) hashgrid[hashcell].particles[n] = k;
|
|
else printf("Too many particles in hash cell (%i, %i), try increasing HASHMAX\n", i, j);
|
|
hashgrid[hashcell].number++;
|
|
particle[k].hashcell = hashcell;
|
|
|
|
if (n > max) max = n;
|
|
}
|
|
|
|
if(verbose) printf("Maximal number of particles per hash cell: %i\n", max);
|
|
}
|
|
|
|
|
|
int wrap_particle(t_particle* particle, double *px, double *py)
|
|
/* relocate particles in case of periodic and similar boundary conditions */
|
|
{
|
|
double x, y, x1, y1, x2, y2;
|
|
int move = 0;
|
|
|
|
x = particle->xc;
|
|
y = particle->yc;
|
|
x1 = x;
|
|
y1 = y;
|
|
|
|
switch (bc_grouped(BOUNDARY_COND)) {
|
|
case (0): /* rectangular b.c. */
|
|
{
|
|
/* do nothing */
|
|
return(0);
|
|
break;
|
|
}
|
|
case (1): /* periodic b.c. */
|
|
{
|
|
if (x < BCXMIN)
|
|
{
|
|
x1 += BCXMAX - BCXMIN;
|
|
move++;
|
|
}
|
|
else if (x > BCXMAX)
|
|
{
|
|
x1 += BCXMIN - BCXMAX;
|
|
move++;
|
|
}
|
|
if (y > BCYMAX)
|
|
{
|
|
y1 += BCYMIN - BCYMAX;
|
|
move++;
|
|
}
|
|
else if (y < BCYMIN)
|
|
{
|
|
y1 += BCYMAX - BCYMIN;
|
|
move++;
|
|
}
|
|
particle->xc = x1;
|
|
particle->yc = y1;
|
|
return(move);
|
|
break;
|
|
}
|
|
case (2): /* Klein bottle b.c. */
|
|
{
|
|
if (y > BCYMAX)
|
|
{
|
|
y1 += BCYMIN - BCYMAX;
|
|
move++;
|
|
}
|
|
else if (y < BCYMIN)
|
|
{
|
|
y1 += BCYMAX - BCYMIN;
|
|
move++;
|
|
}
|
|
if (x < BCXMIN)
|
|
{
|
|
x1 += BCXMAX - BCXMIN;
|
|
y1 = -y1;
|
|
*py *= -1.0;
|
|
particle->angle = DPI - particle->angle;
|
|
move++;
|
|
}
|
|
else if (x > BCXMAX)
|
|
{
|
|
x1 += BCXMIN - BCXMAX;
|
|
y1 = -y1;
|
|
particle->angle = DPI - particle->angle;
|
|
*py *= -1.0;
|
|
move++;
|
|
}
|
|
particle->xc = x1;
|
|
particle->yc = y1;
|
|
return(move);
|
|
break;
|
|
}
|
|
case (3): /* Boy surface b.c. */
|
|
{
|
|
if ((y < BCYMAX)&&(y > BCYMIN))
|
|
{
|
|
if (x > BCXMAX)
|
|
{
|
|
x1 += BCXMIN - BCXMAX;
|
|
y1 = -y1;
|
|
particle->angle *= -1.0;
|
|
particle->vy *= -1.0;
|
|
*py *= -1.0;
|
|
move++;
|
|
}
|
|
else if (x < BCXMIN)
|
|
{
|
|
x1 += BCXMAX - BCXMIN;
|
|
y1 = -y1;
|
|
particle->angle *= -1.0;
|
|
particle->vy *= -1.0;
|
|
*py *= -1.0;
|
|
move++;
|
|
}
|
|
}
|
|
// x = x1;
|
|
// y = y1;
|
|
if ((x < BCXMAX)&&(x > BCXMIN))
|
|
{
|
|
if (y > BCYMAX)
|
|
{
|
|
y1 += BCYMIN - BCYMAX;
|
|
x1 = -x1;
|
|
particle->angle = PI - particle->angle;
|
|
particle->vx *= -1.0;
|
|
*px *= -1.0;
|
|
move++;
|
|
}
|
|
else if (y < BCYMIN)
|
|
{
|
|
y1 += BCYMAX - BCYMIN;
|
|
x1 = -x1;
|
|
particle->angle = PI - particle->angle;
|
|
particle->vx *= -1.0;
|
|
*px *= -1.0;
|
|
move++;
|
|
}
|
|
}
|
|
if (((x >= BCXMAX)||(x <= BCXMIN))&&((y >= BCYMAX)||(y <= BCYMIN)))
|
|
{
|
|
/* This case can lead to numerical instabilities, and can be avoided by putting obstacles in the corners */
|
|
printf("Double wrap!\n");
|
|
if (x >= BCXMAX) x1 = BCXMAX - BCXMIN - x1;
|
|
else x1 = -BCXMAX + BCXMIN - x1;
|
|
if (y >= BCXMAX) y1 = BCYMAX - BCYMIN - y1;
|
|
else y1 = -BCYMAX + BCYMIN - y1;
|
|
|
|
particle->vx *= -1.0;
|
|
*px *= -1.0;
|
|
particle->vy *= -1.0;
|
|
*py *= -1.0;
|
|
|
|
|
|
// if (x1 >= BCXMAX) x1 = BCXMAX - 1.0e-5;
|
|
// else if (x1 <= BCXMIN) x1 = BCXMIN + 1.0e-5;
|
|
// if (y1 >= BCXMAX) y1 = BCYMAX - 1.0e-5;
|
|
// else if (y1 <= BCYMIN) y1 = BCYMIN + 1.0e-5;
|
|
move++;
|
|
}
|
|
|
|
particle->xc = x1;
|
|
particle->yc = y1;
|
|
return(move);
|
|
break;
|
|
}
|
|
case (4): /* genus two (L-shaped domain) b.c. */
|
|
{
|
|
if ((x > 0.0)&&(y > 0.0))
|
|
{
|
|
if (x > y) y1 -= 0.5*(BCYMAX - BCYMIN);
|
|
else x1 -= 0.5*(BCXMAX - BCXMIN);
|
|
move++;
|
|
}
|
|
|
|
else
|
|
{
|
|
|
|
if (x < BCXMIN)
|
|
{
|
|
if (y < 0.0) x1 += BCXMAX - BCXMIN;
|
|
else x1 += 0.5*(BCXMAX - BCXMIN);
|
|
move++;
|
|
}
|
|
else
|
|
{
|
|
if ((y < 0.0)&&(x > BCXMAX))
|
|
{
|
|
x1 += BCXMIN - BCXMAX;
|
|
move++;
|
|
}
|
|
else if ((y >= 0.0)&&(x > 0.0)&&(x < OBSTACLE_RADIUS))
|
|
{
|
|
x1 += 0.5*(BCXMIN - BCXMAX);
|
|
move++;
|
|
}
|
|
}
|
|
|
|
if (y < BCYMIN)
|
|
{
|
|
if (x1 < 0.0) y1 += BCYMAX - BCYMIN;
|
|
else y1 += 0.5*(BCYMAX - BCYMIN);
|
|
move++;
|
|
}
|
|
else
|
|
{
|
|
if ((x1 < 0.0)&&(y > BCYMAX))
|
|
{
|
|
y1 += BCYMIN - BCYMAX;
|
|
move++;
|
|
}
|
|
else if ((x1 >= 0.0)&&(y > 0.0)&&(y < OBSTACLE_RADIUS))
|
|
{
|
|
y1 += 0.5*(BCYMIN - BCYMAX);
|
|
move++;
|
|
}
|
|
}
|
|
}
|
|
|
|
// if (move > 0) printf("Moved particle from (%.3lg, %.3lg) to (%.3lg, %.3lg)\n", x, y, x1, y1);
|
|
|
|
particle->xc = x1;
|
|
particle->yc = y1;
|
|
return(move);
|
|
break;
|
|
}
|
|
default:
|
|
{
|
|
/* do nothing */
|
|
break;
|
|
}
|
|
}
|
|
|
|
}
|
|
|
|
|
|
void wrap_relative_positions(double x1, double y1, double *x2, double *y2)
|
|
/* computes relative positions of particles, taking boundary conditions into account */
|
|
{
|
|
double dxhalf = 0.5*(BCXMAX - BCXMIN), dyhalf = 0.5*(BCYMAX - BCYMIN);
|
|
double dx, dy, x3, y3, x4, y4, dh, dv, wrap_x, wrap_y;
|
|
int verbose = 0;
|
|
|
|
dx = *x2 - x1;
|
|
dy = *y2 - y1;
|
|
|
|
switch (bc_grouped(BOUNDARY_COND)) {
|
|
case (0): /* rectangular b.c. */
|
|
{
|
|
/* do nothing */
|
|
break;
|
|
}
|
|
case (1): /* periodic b.c. */
|
|
{
|
|
if (dx > dxhalf) *x2 -= BCXMAX - BCXMIN;
|
|
else if (-dx > dxhalf) *x2 += BCXMAX - BCXMIN;
|
|
if (dy > dyhalf) *y2 -= BCYMAX - BCYMIN;
|
|
else if (-dy > dyhalf) *y2 += BCYMAX - BCYMIN;
|
|
// printf("(x2,y2) = (%.3lg, %.3lg)\n", *x2, *y2);
|
|
break;
|
|
}
|
|
case (2): /* Klein bottle b.c. */
|
|
{
|
|
if (dx > dxhalf)
|
|
{
|
|
*x2 -= BCXMAX - BCXMIN;
|
|
*y2 *= -1.0;
|
|
}
|
|
else if (-dx > dxhalf)
|
|
{
|
|
*x2 += BCXMAX - BCXMIN;
|
|
*y2 *= -1.0;
|
|
}
|
|
dy = *y2 - y1;
|
|
if (dy > dyhalf) *y2 -= BCYMAX - BCYMIN;
|
|
else if (-dy > dyhalf) *y2 += BCYMAX - BCYMIN;
|
|
break;
|
|
}
|
|
case (3): /* Boy surface - TO FIX */
|
|
{
|
|
// wrap_x = 2.0*(BCXMAX - BCXMIN)/(double)HASHX;
|
|
// wrap_y = 2.0*(BCYMAX - BCYMIN)/(double)HASHY;
|
|
wrap_x = dxhalf;
|
|
wrap_y = dyhalf;
|
|
|
|
/* find which wrapped point is closest to (x1, y1) */
|
|
dh = 100.0;
|
|
if (dx > wrap_x)
|
|
{
|
|
x3 = *x2 - BCXMAX + BCXMIN;
|
|
y3 = -*y2;
|
|
dh = (x3-x1)*(x3-x1) + (y3-y1)*(y3-y1);
|
|
}
|
|
else if (-dx > wrap_x)
|
|
{
|
|
x3 = *x2 + BCXMAX - BCXMIN;
|
|
y3 = -*y2;
|
|
dh = (x3-x1)*(x3-x1) + (y3-y1)*(y3-y1);
|
|
}
|
|
if ((verbose)&&(dh < 100.0)&&(dh > 0.0))
|
|
{
|
|
printf("Case 1: (x1,y1) = (%.3lg, %.3lg)\t", x1, y1);
|
|
printf("(x2,y2) = (%.3lg, %.3lg) -> (%.3lg, %.3lg)\n", *x2, *y2, x3, y3);
|
|
}
|
|
|
|
dv = 100.0;
|
|
if (dy > wrap_y)
|
|
{
|
|
x4 = -*x2;
|
|
y4 = *y2 - BCYMAX + BCYMIN;
|
|
dv = (x4-x1)*(x4-x1) + (y4-y1)*(y4-y1);
|
|
}
|
|
else if (-dy > wrap_y)
|
|
{
|
|
x4 = -*x2;
|
|
y4 = *y2 + BCYMAX - BCYMIN;
|
|
dv = (x4-x1)*(x4-x1) + (y4-y1)*(y4-y1);
|
|
}
|
|
if ((verbose)&&(dv < 100.0)&&(dv > 0.0))
|
|
{
|
|
printf("Case 2: (x1,y1) = (%.3lg, %.3lg)\t", x1, y1);
|
|
printf("(x2,y2) = (%.3lg, %.3lg) -> (%.3lg, %.3lg)\n", *x2, *y2, x4, y4);
|
|
}
|
|
|
|
if (dh < 100.0)
|
|
{
|
|
*x2 = x3;
|
|
*y2 = y3;
|
|
}
|
|
if (dv < dh)
|
|
{
|
|
*x2 = x4;
|
|
*y2 = y4;
|
|
}
|
|
break;
|
|
}
|
|
case (4): /* genus 2 (L-shaped) b.c. */
|
|
{
|
|
x3 = *x2;
|
|
y3 = *y2;
|
|
if ((x1 < 0.0)&&(y1 < 0.0))
|
|
{
|
|
if (dx > dxhalf) *x2 -= (BCXMAX - BCXMIN);
|
|
if (dy > dyhalf) *y2 -= (BCYMAX - BCYMIN);
|
|
}
|
|
else if ((x1 >= 0.0)&&(y1 < 0.0))
|
|
{
|
|
if (dx < -dxhalf) *x2 += (BCXMAX - BCXMIN);
|
|
if (dy > 0.5*dyhalf) *y2 -= 0.5*(BCYMAX - BCYMIN);
|
|
else if (dy < -0.5*dyhalf) *y2 += 0.5*(BCYMAX - BCYMIN);
|
|
}
|
|
else if ((x1 < 0.0)&&(y1 >= 0.0))
|
|
{
|
|
if (dy < -dyhalf) *y2 += BCYMAX - BCYMIN;
|
|
if (dx > 0.5*dxhalf) *x2 -= 0.5*(BCXMAX - BCXMIN);
|
|
else if (dx < -0.5*dxhalf) *x2 += 0.5*(BCXMAX - BCXMIN);
|
|
}
|
|
if ((verbose)&&(module2(*x2 - x1, *y2 - y1) > dyhalf))
|
|
{
|
|
printf("(x1,y1) = (%.3lg, %.3lg)\n", x1, y1);
|
|
printf("(x2,y2) = (%.3lg, %.3lg) ", x3, y3);
|
|
printf("-> (%.3lg, %.3lg)\n", *x2, *y2);
|
|
}
|
|
// printf("(x2,y2) = (%.3lg, %.3lg)\n", *x2, *y2);
|
|
break;
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
void compute_relative_positions(t_particle particle[NMAXCIRCLES], t_hashgrid hashgrid[HASHX*HASHY])
|
|
{
|
|
int j, i0, j0, m0, k, m, p, q, n = 0;
|
|
double x1, x2, y1, y2, xtemp, ytemp;
|
|
|
|
for (j=0; j < ncircles; j++) if (particle[j].active)
|
|
{
|
|
// i0 = particle[j].hashx;
|
|
// j0 = particle[j].hashy;
|
|
// m0 = mhash(i0, j0);
|
|
m0 = particle[j].hashcell;
|
|
x1 = particle[j].xc;
|
|
y1 = particle[j].yc;
|
|
n = 0;
|
|
|
|
for (q=0; q<hashgrid[m0].nneighb; q++)
|
|
{
|
|
m = hashgrid[m0].neighbour[q];
|
|
for (k=0; k<hashgrid[m].number; k++)
|
|
if ((hashgrid[m].particles[k]!=j)&&(particle[hashgrid[m].particles[k]].active))
|
|
{
|
|
if (n < 9*HASHMAX)
|
|
{
|
|
p = hashgrid[m].particles[k];
|
|
|
|
x2 = particle[p].xc;
|
|
y2 = particle[p].yc;
|
|
xtemp = x2;
|
|
ytemp = y2;
|
|
|
|
if (bc_grouped(BOUNDARY_COND) != 0) wrap_relative_positions(x1, y1, &x2, &y2);
|
|
|
|
particle[j].hashneighbour[n] = p;
|
|
particle[j].deltax[n] = x2 - x1;
|
|
particle[j].deltay[n] = y2 - y1;
|
|
// if ((j%50 == 0)&&((vabs(x2-x1)>1.0)||(vabs(y2-y1)>1.0)))
|
|
// if (((vabs(x2-x1)>0.7)||(vabs(y2-y1)>0.7)))
|
|
// {
|
|
// printf("(x1, y1) = (%.3lg, %.3lg), (x2, y2) = (%.3lg, %.3lg) -> (%.3lg, %.3lg)\n", x1, y1, xtemp, ytemp, x2, y2);
|
|
// printf("particle[%i].delta[%i] = (%.4lg, %.4lg)\n", j, n, particle[j].deltax[n], particle[j].deltay[n]);
|
|
// }
|
|
n++;
|
|
}
|
|
else printf("Not enough memory in particle.deltax, particle.deltay\n");
|
|
}
|
|
}
|
|
particle[j].hash_nneighb = n;
|
|
}
|
|
}
|
|
|