---

advect.h File Reference

This graph shows which files directly or indirectly include this file:

Go to the source code of this file.

Functions
void	advect (interface in, real2D u, real2D v, real2D ap, int nx, int ny, real *puaxis)
void	repulse (interface in, real cutoff, real coef, real lredis)

---

Function Documentation

void advect (  interface    in, real2D    u, real2D    v, real2D    ap, int    nx, int    ny, real *    puaxis )

Definition at line 81 of file advect.c. References bilinu(), bilinv(), extra_velocity_normals(), interface::n, real, real2D, splint1, interface::spx, interface::spy, sq, interface::t, UNDEFINED, interface::x, and interface::y. Referenced by timestep(). { int l; real x, y, dx, dy, n1, n2; real a11, a22, a12, a21; // FILE *fptr = fopen("markspeed", "at"); for (l = 0; l < in.n; l++) { x = in.x[l]; y = in.y[l]; if ((dx = bilinu(u, x, y)) == UNDEFINED || (dy = bilinv(v, x, y)) == UNDEFINED) { printf("advect(): extrapolating point %g,%g ... ", x, y); if (l == in.n - 1) { n1 = - splint1(in.spy[l-1], in.t[l]); n2 = splint1(in.spx[l-1], in.t[l]); } else { n1 = - splint1(in.spy[l], in.t[l]); n2 = splint1(in.spx[l], in.t[l]); } a11 = sqrt(sq(n1) + sq(n2)); n1 /= a11; n2 /= a11; printf("normal(%g,%g) ... ", n1, n2); extra_velocity_normals(x, y, n1, n2, x, y, u, v, ap, &dx, &dy, &a11, &a22, &a12, &a21, nx, ny); printf("Ok\n"); } #ifdef CORRECTOR x = in.x[l] + dx; y = in.y[l] + dy; if ((dx = bilinu(u, x, y)) == UNDEFINED || (dy = bilinv(v, x, y)) == UNDEFINED) { printf("advect(): extrapolating point %g,%g ... ", x, y); if (l == in.n - 1) { n1 = - splint1(in.spy[l-1], in.t[l]); n2 = splint1(in.spx[l-1], in.t[l]); } else { n1 = - splint1(in.spy[l], in.t[l]); n2 = splint1(in.spx[l], in.t[l]); } a11 = sqrt(sq(n1) + sq(n2)); n1 /= a11; n2 /= a11; printf("normal(%g,%g) ... ", n1, n2); extra_velocity_normals(x, y, n1, n2, x, y, u, v, ap, &dx, &dy, &a11, &a22, &a12, &a21, nx, ny); printf("Ok\n"); } #endif in.x[l] += dx; in.y[l] += dy; // printf("%d %g %g\n", l, dx, dy); /*cld*/ // if (l == 0) // *puaxis = dx; /*cld*/ // fprintf(fptr, "%g %g\n", in.t[l]/in.t[in.n-1], sqrt(sq(dx) + sq(dy))); // fprintf(fptr, "%g %g %g\n", in.t[l]/in.t[in.n-1], dx, dy); } // fprintf(fptr, "&\n"); // fclose(fptr); }

void repulse (  interface    in, real    cutoff, real    coef, real    lredis )

Definition at line 187 of file advect.c. References AXINATURAL, AXITRUE, interface::bc, interface::n, real, sq, interface::x, xi, interface::y, and yi. { int i, j, imin, imax, imin1, imax1, index; real xm, ym, xi, yi, t, r2; real *fx = (real *)calloc(in.n, sizeof(real)); real *fy = (real *)calloc(in.n, sizeof(real)); real *abx = (real *)malloc((in.n - 1)*sizeof(real)); real *aby = (real *)malloc((in.n - 1)*sizeof(real)); real *abn = (real *)malloc((in.n - 1)*sizeof(real)); #if 1 int *isx = (int *)malloc(in.n*sizeof(int)); int *isy = (int *)malloc(in.n*sizeof(int)); int *is; FILE *fptr = fopen("repulse", "wt"); cutoff *= cutoff; for (i = 0; i < in.n - 1; i++) { abx[i] = in.x[i+1] - in.x[i]; aby[i] = in.y[i+1] - in.y[i]; abn[i] = sqrt(sq(abx[i]) + sq(aby[i])); } /* sort the markers by ascending coordinates */ for (i = 0; i < in.n; i++) isx[i] = isy[i] = i; coord = in.x; qsort(isx, in.n, sizeof(int), cmp); coord = in.y; qsort(isy, in.n, sizeof(int), cmp); for (i = 0; i < in.n; i++) { xm = in.x[i]; ym = in.y[i]; /* find the points with an x-coordinate in [xm-lredis,xm+lredis] */ imax = find_smaller(in.x, isx, in.n, xm + lredis); imin = find_larger(in.x, isx, in.n, xm - lredis); /* find the points with an y-coordinate in [ym-lredis,ym+lredis] */ imax1 = find_smaller(in.y, isy, in.n, ym + lredis); imin1 = find_larger(in.y, isy, in.n, ym - lredis); if (imax - imin > imax1 - imin1) { imax = imax1; imin = imin1; is = isy; } else { is = isx; } /* for each point in the smallest interval */ for (j = imin; j <= imax; j++) { index = is[j]; if (index != i && index+1 != i) { if (index < in.n - 1) { t = ((xm - in.x[index])*abx[index] + (ym - in.y[index])*aby[index])/abn[index]; if (t <= 1.0) { if (t <= 0.0) { xi = in.x[index]; yi = in.y[index]; } else { xi = in.x[index] + t*abx[index]; yi = in.y[index] + t*aby[index]; } r2 = sq(xm - xi) + sq(ym - yi); if (r2 < cutoff) { r2 = coef*1./sq(r2); fx[i] -= r2*(xi - xm); fy[i] -= r2*(yi - ym); fprintf(fptr, "%g %g\n%g %g\n&\n\n", xi, yi, xm, ym); } } } } } } in.x[0] += fx[0]; in.x[in.n-1] += fx[in.n-1]; if (in.bc != AXITRUE && in.bc != AXINATURAL) { in.y[0] += fy[0]; in.y[in.n-1] += fy[in.n-1]; } for (i = 1; i < in.n - 1; i++) { in.x[i] += fx[i]; in.y[i] += fy[i]; } free(fx); free(fy); free(abx); free(aby); free(abn); free(isx); free(isy); fclose(fptr); #else FILE *fptr = fopen("repulse", "wt"); cutoff *= cutoff; for (i = 0; i < in.n - 1; i++) { abx[i] = in.x[i+1] - in.x[i]; aby[i] = in.y[i+1] - in.y[i]; abn[i] = sqrt(sq(abx[i]) + sq(aby[i])); } for (i = 0; i < in.n; i++) { xm = in.x[i]; ym = in.y[i]; for (j = 0; j < in.n - 1; j++) if (j != i && j+1 != i) { t = ((xm - in.x[j])*abx[j] + (ym - in.y[j])*aby[j])/abn[j]; if (t <= 1.0) { if (t < 0.0) { xi = in.x[j]; yi = in.y[j]; } else { xi = in.x[j] + t*abx[j]; yi = in.y[j] + t*aby[j]; } r2 = sq(xm - xi) + sq(ym - yi); if (r2 < cutoff) { r2 = coef*1./sq(r2); fx[i] -= r2*(xi - xm); fy[i] -= r2*(yi - ym); fprintf(fptr, "%g %g\n%g %g\n&\n\n", xi, yi, xm, ym); } } } } fclose(fptr); in.x[0] += fx[0]; in.x[in.n-1] += fx[in.n-1]; if (in.bc != AXITRUE && in.bc != AXINATURAL) { in.y[0] += fy[0]; in.y[in.n-1] += fy[in.n-1]; } for (i = 1; i < in.n - 1; i++) { in.x[i] += fx[i]; in.y[i] += fy[i]; } free(fx); free(fy); free(abx); free(aby); free(abn); #endif }

---