#include "utilf.h"#include "markers.h"#include "advect.h"#include "extra.h"#include "inout.h"Include dependency graph for advect.c:
Go to the source code of this file.
Defines | |
| #define | OUTU(i, j) (u[i][j] == UNDEFINED) |
| #define | OUTV(i, j) (v[i][j] == UNDEFINED) |
Functions | |
| real | bilinu (real2D u, real x, real y) |
| real | bilinv (real2D v, real x, real y) |
| 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) |
|
|
Definition at line 7 of file advect.c. Referenced by bilinu(). |
|
|
Definition at line 8 of file advect.c. Referenced by bilinv(). |
|
||||||||||||||||||||||||||||||||
|
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().
00083 {
00084 int l;
00085 real x, y, dx, dy, n1, n2;
00086 real a11, a22, a12, a21;
00087 // FILE *fptr = fopen("markspeed", "at");
00088
00089 for (l = 0; l < in.n; l++) {
00090 x = in.x[l]; y = in.y[l];
00091 if ((dx = bilinu(u, x, y)) == UNDEFINED ||
00092 (dy = bilinv(v, x, y)) == UNDEFINED) {
00093 printf("advect(): extrapolating point %g,%g ... ", x, y);
00094 if (l == in.n - 1) {
00095 n1 = - splint1(in.spy[l-1], in.t[l]);
00096 n2 = splint1(in.spx[l-1], in.t[l]);
00097 }
00098 else {
00099 n1 = - splint1(in.spy[l], in.t[l]);
00100 n2 = splint1(in.spx[l], in.t[l]);
00101 }
00102 a11 = sqrt(sq(n1) + sq(n2));
00103 n1 /= a11; n2 /= a11;
00104 printf("normal(%g,%g) ... ", n1, n2);
00105 extra_velocity_normals(x, y, n1, n2, x, y,
00106 u, v, ap, &dx, &dy,
00107 &a11, &a22, &a12, &a21, nx, ny);
00108 printf("Ok\n");
00109 }
00110 #ifdef CORRECTOR
00111 x = in.x[l] + dx; y = in.y[l] + dy;
00112 if ((dx = bilinu(u, x, y)) == UNDEFINED ||
00113 (dy = bilinv(v, x, y)) == UNDEFINED) {
00114 printf("advect(): extrapolating point %g,%g ... ", x, y);
00115 if (l == in.n - 1) {
00116 n1 = - splint1(in.spy[l-1], in.t[l]);
00117 n2 = splint1(in.spx[l-1], in.t[l]);
00118 }
00119 else {
00120 n1 = - splint1(in.spy[l], in.t[l]);
00121 n2 = splint1(in.spx[l], in.t[l]);
00122 }
00123 a11 = sqrt(sq(n1) + sq(n2));
00124 n1 /= a11; n2 /= a11;
00125 printf("normal(%g,%g) ... ", n1, n2);
00126 extra_velocity_normals(x, y, n1, n2, x, y,
00127 u, v, ap, &dx, &dy,
00128 &a11, &a22, &a12, &a21, nx, ny);
00129 printf("Ok\n");
00130 }
00131 #endif
00132 in.x[l] += dx;
00133 in.y[l] += dy;
00134 // printf("%d %g %g\n", l, dx, dy);
00135
00136 /*cld*/
00137 // if (l == 0)
00138 // *puaxis = dx;
00139 /*cld*/
00140 // fprintf(fptr, "%g %g\n", in.t[l]/in.t[in.n-1], sqrt(sq(dx) + sq(dy)));
00141 // fprintf(fptr, "%g %g %g\n", in.t[l]/in.t[in.n-1], dx, dy);
00142 }
00143
00144 // fprintf(fptr, "&\n");
00145 // fclose(fptr);
00146
00147 }
|
|
||||||||||||||||
|
Definition at line 13 of file advect.c. References OUTU, real, real2D, and UNDEFINED. Referenced by advect().
00014 {
00015 int i, j;
00016 i = x; j = (int) ((real)y - 0.5);
00017 x -= (real)i;
00018 y -= 0.5 + (real)j;
00019 if (OUTU(i,j)) {
00020 x = 1. - x; y = 1. - y;
00021 if (x + y > 1.0 || OUTU(i,j+1) || OUTU(i+1,j) || OUTU(i+1,j+1))
00022 return UNDEFINED;
00023 return x*u[i][j+1] + y*u[i+1][j] + (1. - x - y)*u[i+1][j+1];
00024 }
00025 if (OUTU(i+1,j+1)) {
00026 if (x + y > 1.0 || OUTU(i+1,j) || OUTU(i,j+1))
00027 return UNDEFINED;
00028 return x*u[i+1][j] + y*u[i][j+1] + (1. - x - y)*u[i][j];
00029 }
00030 if (OUTU(i+1,j)) {
00031 y = 1. - y;
00032 if (x + y > 1.0 || OUTU(i,j+1))
00033 return UNDEFINED;
00034 return x*u[i+1][j+1] + y*u[i][j] + (1. - x - y)*u[i][j+1];
00035 }
00036 if (OUTU(i,j+1)) {
00037 x = 1. - x;
00038 if (x + y > 1.0)
00039 return UNDEFINED;
00040 return x*u[i][j] + y*u[i+1][j+1] + (1. - x - y)*u[i+1][j];
00041 }
00042 return u[i][j]*(1. - x - y + x*y) + u[i+1][j]*x*(1. - y) +
00043 u[i][j+1]*y*(1. - x) + u[i+1][j+1]*x*y;
00044 }
|
|
||||||||||||||||
|
Definition at line 47 of file advect.c. References OUTV, real, real2D, and UNDEFINED. Referenced by advect().
00048 {
00049 int i, j;
00050 i = (int)((real) x - 0.5); j = y;
00051 x -= 0.5 + (real)i;
00052 y -= (real) j;
00053 if (OUTV(i,j)) {
00054 x = 1. - x; y = 1. - y;
00055 if (x + y > 1.0 || OUTV(i,j+1) || OUTV(i+1,j) || OUTV(i+1,j+1))
00056 return UNDEFINED;
00057 return x*v[i][j+1] + y*v[i+1][j] + (1. - x - y)*v[i+1][j+1];
00058 }
00059 if (OUTV(i+1,j+1)) {
00060 if (x + y > 1.0 || OUTV(i,j+1) || OUTV(i+1,j))
00061 return UNDEFINED;
00062 return x*v[i+1][j] + y*v[i][j+1] + (1. - x - y)*v[i][j];
00063 }
00064 if (OUTV(i+1,j)) {
00065 y = 1. - y;
00066 if (x + y > 1.0 || OUTV(i,j+1))
00067 return UNDEFINED;
00068 return x*v[i+1][j+1] + y*v[i][j] + (1. - x - y)*v[i][j+1];
00069 }
00070 if (OUTV(i,j+1)) {
00071 x = 1. - x;
00072 if (x + y > 1.0)
00073 return UNDEFINED;
00074 return x*v[i][j] + y*v[i+1][j+1] + (1. - x - y)*v[i+1][j];
00075 }
00076 return v[i][j]*(1. - x - y + x*y) + v[i+1][j]*x*(1. - y) +
00077 v[i][j+1]*y*(1. - x) + v[i+1][j+1]*x*y;
00078 }
|
|
||||||||||||||||||||
|
Definition at line 187 of file advect.c. References AXINATURAL, AXITRUE, interface::bc, interface::n, real, sq, interface::x, xi, interface::y, and yi.
00188 {
00189 int i, j, imin, imax, imin1, imax1, index;
00190 real xm, ym, xi, yi, t, r2;
00191 real *fx = (real *)calloc(in.n, sizeof(real));
00192 real *fy = (real *)calloc(in.n, sizeof(real));
00193 real *abx = (real *)malloc((in.n - 1)*sizeof(real));
00194 real *aby = (real *)malloc((in.n - 1)*sizeof(real));
00195 real *abn = (real *)malloc((in.n - 1)*sizeof(real));
00196 #if 1
00197 int *isx = (int *)malloc(in.n*sizeof(int));
00198 int *isy = (int *)malloc(in.n*sizeof(int));
00199 int *is;
00200 FILE *fptr = fopen("repulse", "wt");
00201
00202 cutoff *= cutoff;
00203 for (i = 0; i < in.n - 1; i++) {
00204 abx[i] = in.x[i+1] - in.x[i];
00205 aby[i] = in.y[i+1] - in.y[i];
00206 abn[i] = sqrt(sq(abx[i]) + sq(aby[i]));
00207 }
00208
00209 /* sort the markers by ascending coordinates */
00210 for (i = 0; i < in.n; i++)
00211 isx[i] = isy[i] = i;
00212 coord = in.x;
00213 qsort(isx, in.n, sizeof(int), cmp);
00214 coord = in.y;
00215 qsort(isy, in.n, sizeof(int), cmp);
00216
00217 for (i = 0; i < in.n; i++) {
00218 xm = in.x[i]; ym = in.y[i];
00219 /* find the points with an x-coordinate in [xm-lredis,xm+lredis] */
00220 imax = find_smaller(in.x, isx, in.n, xm + lredis);
00221 imin = find_larger(in.x, isx, in.n, xm - lredis);
00222 /* find the points with an y-coordinate in [ym-lredis,ym+lredis] */
00223 imax1 = find_smaller(in.y, isy, in.n, ym + lredis);
00224 imin1 = find_larger(in.y, isy, in.n, ym - lredis);
00225 if (imax - imin > imax1 - imin1) { imax = imax1; imin = imin1; is = isy; }
00226 else { is = isx; }
00227 /* for each point in the smallest interval */
00228 for (j = imin; j <= imax; j++) {
00229 index = is[j];
00230 if (index != i && index+1 != i) {
00231 if (index < in.n - 1) {
00232 t = ((xm - in.x[index])*abx[index] +
00233 (ym - in.y[index])*aby[index])/abn[index];
00234 if (t <= 1.0) {
00235 if (t <= 0.0) { xi = in.x[index]; yi = in.y[index]; }
00236 else {
00237 xi = in.x[index] + t*abx[index];
00238 yi = in.y[index] + t*aby[index];
00239 }
00240 r2 = sq(xm - xi) + sq(ym - yi);
00241 if (r2 < cutoff) {
00242 r2 = coef*1./sq(r2);
00243 fx[i] -= r2*(xi - xm);
00244 fy[i] -= r2*(yi - ym);
00245 fprintf(fptr, "%g %g\n%g %g\n&\n\n", xi, yi, xm, ym);
00246 }
00247 }
00248 }
00249 }
00250 }
00251 }
00252 in.x[0] += fx[0]; in.x[in.n-1] += fx[in.n-1];
00253 if (in.bc != AXITRUE && in.bc != AXINATURAL) {
00254 in.y[0] += fy[0]; in.y[in.n-1] += fy[in.n-1];
00255 }
00256 for (i = 1; i < in.n - 1; i++) {
00257 in.x[i] += fx[i]; in.y[i] += fy[i];
00258 }
00259 free(fx); free(fy);
00260 free(abx); free(aby); free(abn);
00261 free(isx); free(isy);
00262 fclose(fptr);
00263 #else
00264 FILE *fptr = fopen("repulse", "wt");
00265
00266 cutoff *= cutoff;
00267 for (i = 0; i < in.n - 1; i++) {
00268 abx[i] = in.x[i+1] - in.x[i];
00269 aby[i] = in.y[i+1] - in.y[i];
00270 abn[i] = sqrt(sq(abx[i]) + sq(aby[i]));
00271 }
00272 for (i = 0; i < in.n; i++) {
00273 xm = in.x[i]; ym = in.y[i];
00274 for (j = 0; j < in.n - 1; j++)
00275 if (j != i && j+1 != i) {
00276 t = ((xm - in.x[j])*abx[j] + (ym - in.y[j])*aby[j])/abn[j];
00277 if (t <= 1.0) {
00278 if (t < 0.0) {
00279 xi = in.x[j]; yi = in.y[j];
00280 }
00281 else {
00282 xi = in.x[j] + t*abx[j];
00283 yi = in.y[j] + t*aby[j];
00284 }
00285 r2 = sq(xm - xi) + sq(ym - yi);
00286 if (r2 < cutoff) {
00287 r2 = coef*1./sq(r2);
00288 fx[i] -= r2*(xi - xm);
00289 fy[i] -= r2*(yi - ym);
00290 fprintf(fptr, "%g %g\n%g %g\n&\n\n", xi, yi, xm, ym);
00291 }
00292 }
00293 }
00294 }
00295 fclose(fptr);
00296
00297 in.x[0] += fx[0]; in.x[in.n-1] += fx[in.n-1];
00298 if (in.bc != AXITRUE && in.bc != AXINATURAL) {
00299 in.y[0] += fy[0]; in.y[in.n-1] += fy[in.n-1];
00300 }
00301 for (i = 1; i < in.n - 1; i++) {
00302 in.x[i] += fx[i]; in.y[i] += fy[i];
00303 }
00304 free(fx); free(fy);
00305 free(abx); free(aby); free(abn);
00306 #endif
00307 }
|
1.2.18