This graph shows which files directly or indirectly include this file:
Go to the source code of this file.
Defines | |
| #define | INTERPLATE_FULL 1 |
| #define | INTERPLATE_INTER 2 |
Functions | |
| void | extra_velocity_normals (real x, real y, real n1, real n2, real xn, real yn, real2D u, real2D v, real2D ap, real *eu, real *ev, real *a11, real *a22, real *a12, real *a21, int nx, int ny) |
| void | extra_velocity_normals2 (real x, real y, real n1, real n2, real xn, real yn, real2D u, real2D v, real2D ap, real *eu, real *ev, real *dudx, real *dudy, real *dvdx, real *dvdy, int nx, int ny) |
| void | extra_velocity (real x, real y, real2D u, real2D v, real2D ap, interface in, real *eu, real *ev, int nx, int ny) |
| real | kappa_normals (polynom3 px, polynom3 py, real t, real *nx, real *ny) |
| real | avg_kappa_normals (polynom3 px, polynom3 py, real t1, real t2, real *nx, real *ny) |
| int | closest_normal (real x, real y, interface in, real *xmin, real *ymin, real *pmin, real *tmin, real *xn, real *yn) |
|
|
Definition at line 1 of file extra.h. Referenced by bc_vector(). |
|
|
|
|
||||||||||||||||||||||||||||
|
Definition at line 725 of file extra.c. References real, splint, splint1, splint2, and sq.
00727 {
00728 static real coef[7] = {1., 3., 3., 2., 3., 3., 1.};
00729 int i;
00730 real xp, yp, y, h = (t2 - t1)/6., t, kappasum = 0.0, lsum = 0.0, dl2;
00731
00732 for (i = 0; i < 7; i++) {
00733 t = t1 + h*i;
00734 xp = splint1(px, t);
00735 yp = splint1(py, t);
00736 dl2 = sq(xp) + sq(yp);
00737 kappasum += coef[i]*((xp*splint2(py, t) - yp*splint2(px, t))/dl2);
00738 /* axisymmetric term */
00739 y = splint(py, t) - 2.0;
00740 if (fabs(y) < 1e-2) /* second order approximation near the axis */
00741 kappasum -= coef[i]*splint2(px, t)/splint1(py, t);
00742 else
00743 kappasum -= coef[i]*xp/y;
00744 lsum += coef[i]*sqrt(dl2);
00745 }
00746 *nx = 8.*(splint(py, t2) - splint(py, t1))/(3.*h*lsum);
00747 *ny = - 8.*(splint(px, t2) - splint(px, t1))/(3.*h*lsum);
00748 return kappasum/lsum;
00749 }
|
|
||||||||||||||||||||||||||||||||||||||||
|
Definition at line 51 of file extra.c. References distmin(), interface_write_splines(), interface::n, PREC, real, splint, splint1, interface::spx, interface::spy, sq, interface::t, interface::x, and interface::y. Referenced by extra_poly(), and extra_velocity().
00054 {
00055 int i, span;
00056 real d, dmin, dl;
00057 int imin, found = 0, which;
00058
00059 /* find the closest segment */
00060 #if 0
00061 dmin = normal_dist(in.x[0], in.y[0], in.x[1], in.y[1], x, y, &which);
00062 imin = 0;
00063 for (i = 1; i < in.n - 1; i++) {
00064 d = normal_dist(in.x[i], in.y[i], in.x[i+1], in.y[i+1], x, y, &which);
00065 if (d < dmin) {
00066 dmin = d;
00067 if (which == 0 || which == 1) imin = i;
00068 else imin = i + 1;
00069 }
00070 }
00071 #else
00072 /* find the closest marker */
00073 dmin = sq(x - in.x[0]) + sq(y - in.y[0]); imin = 0;
00074 for (i = 1; i < in.n; i++) {
00075 d = sq(x - in.x[i]) + sq(y - in.y[i]);
00076 if (d < dmin) { dmin = d; imin = i; }
00077 }
00078 #endif
00079
00080 /* find the closest point on the interface */
00081 i = imin; span = 1;
00082 if (i < in.n - 1)
00083 found = distmin(in.spx[i], in.spy[i], x, y, tmin,
00084 in.t[i], in.t[i+1], PREC);
00085 while (!found && (i + span < in.n - 1 || i - span >= 0)) {
00086 if (i + span < in.n - 1 &&
00087 (found = distmin(in.spx[i+span], in.spy[i+span], x, y, tmin,
00088 in.t[i+span], in.t[i+span+1], PREC)))
00089 imin = i + span;
00090 else if (i - span >= 0 &&
00091 (found = distmin(in.spx[i-span], in.spy[i-span], x, y, tmin,
00092 in.t[i-span], in.t[i-span+1], PREC)))
00093 imin = i - span;
00094 span++;
00095 }
00096 if (!found) {
00097 FILE *fptr = fopen("closest_splines", "wt");
00098 interface_write_splines(in, fptr, 1000);
00099 fclose(fptr);
00100 printf("&&&& %g %g %d %g\n", x, y, imin, dmin);
00101 printf("%g %g\n%g %g\n%g %g\n", in.x[imin-1], in.y[imin-1],
00102 in.x[imin], in.y[imin],
00103 in.x[imin+1], in.y[imin+1]);
00104 return 0;
00105 }
00106
00107 *xn = - splint1(in.spy[imin], *tmin);
00108 *yn = splint1(in.spx[imin], *tmin);
00109
00110 dl = sqrt(sq(*xn) + sq(*yn));
00111 *xn /= dl; *yn /= dl;
00112 *xmin = splint(in.spx[imin], *tmin);
00113 *ymin = splint(in.spy[imin], *tmin);
00114 *pmin = (*tmin - in.t[imin])/(in.t[imin+1] - in.t[imin])
00115 *(in.p[imin+1] - in.p[imin]) + in.p[imin];
00116
00117 return 1;
00118 }
|
|
||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 678 of file extra.c. References closest_normal(), extra_velocity_normals(), extra_velocity_normals2(), real, and real2D. Referenced by bc_vector_div().
00683 {
00684 real n1, n2;
00685 real a11, a22, a12, a21;
00686 real xmin, ymin, pmin, tmin;
00687 if (!closest_normal(x, y, in, &xmin, &ymin, &pmin, &tmin, &n1, &n2)) {
00688 fprintf(stderr, "extra_velocity(%g, %g): no closest_normal\n", x, y);
00689 exit(1);
00690 }
00691 #ifdef EXTRA_LINEAR
00692 extra_velocity_normals(x, y, n1, n2, xmin, ymin,
00693 u, v, ap,
00694 eu, ev, &a11, &a22, &a12, &a21,
00695 nx, ny);
00696 #else
00697 extra_velocity_normals2(x, y, n1, n2, xmin, ymin,
00698 u, v, ap,
00699 eu, ev, &a11, &a22, &a12, &a21,
00700 nx, ny);
00701 #endif
00702 }
|
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 212 of file extra.c. References closest_u(), closest_v(), INU, INV, real, real2D, sq, and SQCLOSEST. Referenced by advect(), and extra_velocity().
00218 {
00219 int i, j, k, l, l1, kmin, kmax, lmin, lmax;
00220 int nsx = 0, nsy = 0, dx, dy;
00221 real xSx = 0., xSxx = 0., xSxy = 0., xSy = 0., xSyy = 0.;
00222 real ySx = 0., ySxx = 0., ySxy = 0., ySy = 0., ySyy = 0.;
00223 real xSu = 0., xSxu = 0., xSyu = 0.;
00224 real ySv = 0., ySxv = 0., ySyv = 0.;
00225 real sys[7][7] = {{0.,0.,0.,0.,0.,0.,0.},
00226 {0.,0.,0.,0.,0.,0.,0.},
00227 {0.,0.,0.,0.,0.,0.,0.},
00228 {0.,0.,0.,0.,0.,0.,0.},
00229 {0.,0.,0.,0.,0.,0.,0.},
00230 {0.,0.,0.,0.,0.,0.,0.},
00231 {0.,0.,0.,0.,0.,0.,0.}};
00232 real rhs[7] = {0.,0.,0.,0.,0.,0.,0.}, sol[7], dummy;
00233 #ifdef DEBUG_EXTRA
00234 FILE *fptr1 = fopen("extra_upoints", "at");
00235 FILE *fptr2 = fopen("extra_vpoints", "at");
00236 FILE *fptr3 = fopen("extra_points", "at");
00237 FILE *fptr4 = fopen("extra_normals", "at");
00238 #endif
00239
00240 #if 0
00241 /* find the velocity points using a 5x5 or 4x4 stencil */
00242 i = x + 0.5; j = y;
00243 dx = dy = 1;
00244
00245 do {
00246 nsx = 0;
00247 xSx = xSxx = xSxy = xSy = xSyy = 0.;
00248 xSu = xSxu = xSyu = 0.;
00249
00250 if (modf(x + 0.5, &dummy) == 0.0) { kmin = i - dx; kmax = i + dx - 1; }
00251 else { kmin = i - dx; kmax = i + dx; }
00252
00253 if (modf(y, &dummy) == 0.0) { lmin = j - dy; lmax = j + dy - 1; }
00254 else { lmin = j - dy; lmax = j + dy; }
00255
00256 for (k = kmin; k <= kmax; k++)
00257 for (l = lmin; l <= lmax; l++) {
00258 l1 = l < 2 ? 3 - l : l;
00259 if (INU(k,l1)) {
00260 xSx += (real)k; xSxx += sq((real)k);
00261 xSxy += (0.5 + (real)l)*(real)k;
00262 xSy += 0.5 + (real)l; xSyy += sq(0.5 + (real)l);
00263 xSu += u[k][l1]; xSxu += u[k][l1]*(real)k;
00264 xSyu += u[k][l1]*(0.5 + (real)l);
00265 #ifdef DEBUG_EXTRA
00266 fprintf(fptr1, "%g %g %g\n", (real)k, 0.5 + (real)l, u[k][l1]);
00267 #endif
00268 nsx++;
00269 }
00270 }
00271 dx++; dy++;
00272 } while((nsx <= kmax - kmin + 1 || nsx <= lmax - lmin + 1) && dx <= 3);
00273
00274 i = x; j = y + 0.5;
00275 dx = dy = 1;
00276 do {
00277 nsy = 0;
00278 ySx = ySxx = ySxy = ySy = ySyy = 0.;
00279 ySv = ySxv = ySyv = 0.;
00280
00281 if (modf(x, &dummy) == 0.0) { kmin = i - dx; kmax = i + dx - 1; }
00282 else { kmin = i - dx; kmax = i + dx; }
00283
00284 if (modf(y + 0.5, &dummy) == 0.0) { lmin = j - dy; lmax = j + dy - 1; }
00285 else { lmin = j - dy; lmax = j + dy; }
00286
00287 for (k = kmin; k <= kmax; k++)
00288 for (l = lmin; l <= lmax; l++) {
00289 l1 = l < 2 ? 4 - l : l;
00290 if (INV(k,l1)) {
00291 ySx += 0.5 + (real)k; ySxx += sq(0.5 + (real)k);
00292 ySxy += (0.5 + (real)k)*(real)l;
00293 ySy += (real)l; ySyy += sq((real)l);
00294 dummy = l1 == l ? v[k][l1] : - v[k][l1];
00295 ySv += dummy; ySxv += dummy*(0.5 + (real)k);
00296 ySyv += dummy*(real)l;
00297 #ifdef DEBUG_EXTRA
00298 fprintf(fptr2, "%g %g %g\n", 0.5+(real)k, (real)l, dummy);
00299 #endif
00300 nsy++;
00301 }
00302 }
00303 dx++; dy++;
00304 } while((nsy <= kmax - kmin + 1 || nsy <= lmax - lmin + 1) && dx <= 3);
00305 #else
00306 real xu[SQCLOSEST], yu[SQCLOSEST], uu[SQCLOSEST];
00307 real xv[SQCLOSEST], yv[SQCLOSEST], vv[SQCLOSEST];
00308 real xmaxu, xminu, ymaxu, yminu;
00309 real xmaxv, xminv, ymaxv, yminv;
00310 int n;
00311
00312 nsx = nsy = 5;
00313 if ((n = closest_u(xn, yn, n1, n2, u, ap, xu, yu, uu,
00314 &xmaxu, &xminu, &ymaxu, &yminu, nx, ny)) < nsx) {
00315 printf("extra_velocity_normals(%g,%g): not enough u's\n", x, y);
00316 printf("xn: %g yn: %g n1: %g n2: %g\n", xn, yn, n1, n2);
00317 for (i = 0; i < n; i++)
00318 printf(" %g %g\n", xu[i], yu[i]);
00319 exit(1);
00320 }
00321 if ((n = closest_v(xn, yn, n1, n2, v, ap, xv, yv, vv,
00322 &xmaxv, &xminv, &ymaxv, &yminv, nx, ny)) < nsy) {
00323 printf("extra_velocity_normals(%g,%g): not enough v's\n", x, y);
00324 printf("xn: %g yn: %g n1: %g n2: %g\n", xn, yn, n1, n2);
00325 for (i = 0; i < n; i++)
00326 printf(" %g %g\n", xv[i], yv[i]);
00327 exit(1);
00328 }
00329
00330 xSx = xSxx = xSxy = xSy = xSyy = 0.;
00331 xSu = xSxu = xSyu = 0.;
00332 for (i = 0; i < nsx; i++) {
00333 xSx += xu[i]; xSxx += sq(xu[i]);
00334 xSxy += xu[i]*yu[i];
00335 xSy += yu[i]; xSyy += sq(yu[i]);
00336 xSu += uu[i]; xSxu += uu[i]*xu[i];
00337 xSyu += uu[i]*yu[i];
00338 #ifdef DEBUG_EXTRA
00339 fprintf(fptr1, "%g %g %g\n", xu[i], yu[i], uu[i]);
00340 #endif
00341 }
00342
00343 ySx = ySxx = ySxy = ySy = ySyy = 0.;
00344 ySv = ySxv = ySyv = 0.;
00345 for (i = 0; i < nsy; i++) {
00346 ySx += xv[i]; ySxx += sq(xv[i]);
00347 ySxy += xv[i]*yv[i];
00348 ySy += yv[i]; ySyy += sq(yv[i]);
00349 ySv += vv[i]; ySxv += vv[i]*xv[i];
00350 ySyv += yv[i]*vv[i];
00351 #ifdef DEBUG_EXTRA
00352 fprintf(fptr2, "%g %g %g\n", xv[i], yv[i], vv[i]);
00353 #endif
00354 }
00355 #endif
00356
00357 /* build the system of equations
00358 sys[i][0]*u0+sys[i][1]*v0+sys[i][2]*A11+sys[i][3]*A22+
00359 sys[i][4]*A12+sys[i][5]*A21+sys[i][6]*lambda=rhs[i] */
00360 sys[0][2] = n1*n2; sys[0][3] = -n1*n2;
00361 sys[0][4] = (sq(n2) - sq(n1))/2.0; sys[0][5] = sys[0][4];
00362 sys[1][0] = nsx; sys[1][2] = xSx; sys[1][4] = xSy;
00363 rhs[1] = xSu;
00364 sys[2][1] = nsy; sys[2][3] = ySy; sys[2][5] = ySx;
00365 rhs[2] = ySv;
00366 sys[3][0] = xSx; sys[3][2] = xSxx; sys[3][4] = xSxy;
00367 sys[3][6] = n1*n2/2.;
00368 rhs[3] = xSxu;
00369 sys[4][1] = ySy; sys[4][3] = ySyy; sys[4][5] = ySxy;
00370 sys[4][6] = - n1*n2/2.;
00371 rhs[4] = ySyv;
00372 sys[5][0] = xSy; sys[5][2] = xSxy; sys[5][4] = xSyy;
00373 sys[5][6] = (sq(n2) - sq(n1))/4.0;
00374 rhs[5] = xSyu;
00375 sys[6][1] = ySx; sys[6][3] = ySxy; sys[6][5] = ySxx;
00376 sys[6][6] = (sq(n2) - sq(n1))/4.0;
00377 rhs[6] = ySxv;
00378
00379 #ifdef DEBUG_EXTRA
00380 fprintf(fptr1, "&\n\n");
00381 fprintf(fptr2, "&\n\n");
00382 fprintf(fptr3, "%g %g\n&\n\n", x, y);
00383 fprintf(fptr4, "%g %g\n%g %g\n&\n\n", xn, yn, xn + n1, yn + n2);
00384 fclose(fptr1); fclose(fptr2); fclose(fptr3); fclose(fptr4);
00385 #endif
00386
00387 if (gauss(sys, rhs, sol, 7)) {
00388 printf("Error: extra_velocity(%g, %g): no solution\n", x, y);
00389 printf(" nsx: %d nsy: %d\n", nsx, nsy);
00390 for (i = 0; i < nsx; i++)
00391 printf("%g %g %g\n", xu[i], yu[i], uu[i]);
00392 printf("\n");
00393 for (i = 0; i < nsy; i++)
00394 printf("%g %g %g\n", xv[i], yv[i], vv[i]);
00395 printf("\n");
00396 for (i = 0; i < 7; i++) {
00397 for (j = 0; j < 7; j++)
00398 printf("%10g ", sys[i][j]);
00399 printf("\n");
00400 }
00401 exit(1);
00402 }
00403
00404 *dudx = sol[2]; *dvdy = sol[3]; *dudy = sol[4]; *dvdx = sol[5];
00405 *eu = sol[0] + sol[2]*x + sol[4]*y;
00406 *ev = sol[1] + sol[5]*x + sol[3]*y;
00407
00408 if (*eu != *eu || *ev != *ev) {
00409 printf("Error: extra_velocity(%g, %g): NaN\n", x, y);
00410 printf(" eu: %g ev: %g\n", *eu, *ev);
00411 printf(" nsx: %d nsy: %d\n", nsx, nsy);
00412 for (i = 0; i < nsx; i++)
00413 printf("%g %g %g\n", xu[i], yu[i], uu[i]);
00414 printf("\n");
00415 for (i = 0; i < nsy; i++)
00416 printf("%g %g %g\n", xv[i], yv[i], vv[i]);
00417 printf("\n");
00418 for (i = 0; i < 7; i++) {
00419 for (j = 0; j < 7; j++)
00420 printf("%10g ", sys[i][j]);
00421 printf("\n");
00422 }
00423 exit(1);
00424 }
00425 }
|
|
||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 428 of file extra.c. References closest_u(), closest_v(), gauss13(), real, real2D, sq, and SQCLOSEST. Referenced by extra_velocity().
00435 {
00436 int i, j, k, l, kmin, kmax, lmin, lmax;
00437 int nsx = 0, nsy = 0, dx, dy;
00438 real x1, y1;
00439 real xSx, xSxx, xSxy, xSy, xSyy, xSxxy, xSxyy, xSyyy, xSxxx,
00440 xSxxxy, xSxxyy, xSxyyy, xSyyyy, xSxxxx, xSu, xSxu, xSyu,
00441 xSxxu, xSyyu, xSxyu;
00442 real ySx, ySxx, ySxy, ySy, ySyy, ySxxy, ySxyy, ySyyy, ySxxx,
00443 ySxxxy, ySxxyy, ySxyyy, ySyyyy, ySxxxx, ySv, ySxv, ySyv,
00444 ySxxv, ySyyv, ySxyv;
00445 real sys[13][13] = {{0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},
00446 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},
00447 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},
00448 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},
00449 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},
00450 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},
00451 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},
00452 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},
00453 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},
00454 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},
00455 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},
00456 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},
00457 {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}};
00458 real rhs[13] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}, sol[13], dummy;
00459 real n1n2 = n1*n2, sqn = (sq(n2) - sq(n1))/2.;
00460 real xu[SQCLOSEST], yu[SQCLOSEST], uu[SQCLOSEST];
00461 real xv[SQCLOSEST], yv[SQCLOSEST], vv[SQCLOSEST];
00462 real xmaxu, ymaxu, xminu, yminu, xmaxv, ymaxv, xminv, yminv, det;
00463 /* #define DEBUG_EXTRA */
00464 #ifdef DEBUG_EXTRA
00465 static int cas;
00466 FILE *fptr1, *fptr2, *fptr3, *fptr4;
00467 char s[256];
00468
00469 sprintf(s, "extra_upoints%d", cas);
00470 fptr1 = fopen(s, "wt");
00471 sprintf(s, "extra_vpoints%d", cas);
00472 fptr2 = fopen(s, "wt");
00473 sprintf(s, "extra_u%d", cas);
00474 fptr3 = fopen(s, "wt");
00475 sprintf(s, "extra_v%d", cas++);
00476 fptr4 = fopen(s, "wt");
00477 #endif
00478
00479 nsx = nsy = 12;
00480 if (closest_u(xn, yn, n1, n2, u, ap, xu, yu, uu,
00481 &xmaxu, &xminu, &ymaxu, &yminu, nx, ny) < nsx) {
00482 fprintf(stderr, "extra_velocity_normals2(%g,%g): not enough u's\n", x, y);
00483 exit(1);
00484 }
00485 if (closest_v(xn, yn, n1, n2, v, ap, xv, yv, vv,
00486 &xmaxv, &xminv, &ymaxv, &yminv, nx, ny) < nsy) {
00487 fprintf(stderr, "extra_velocity_normals2(%g,%g): not enough v's\n", x, y);
00488 exit(1);
00489 }
00490
00491 xSx = xSxx = xSxy = xSy = xSyy = xSxxy = xSxyy = xSyyy = xSxxx =
00492 xSxxxy = xSxxyy = xSxyyy = xSyyyy = xSxxxx = xSu = xSxu = xSyu =
00493 xSxxu = xSyyu = xSxyu = 0.0;
00494 for (i = 0; i < nsx; i++) {
00495 x1 = xu[i]; y1 = yu[i];
00496
00497 xSx += x1; xSxx += sq(x1);
00498 xSxy += x1*y1;
00499 xSy += y1; xSyy += sq(y1);
00500
00501 xSxxy += sq(x1)*y1;
00502 xSxyy += x1*sq(y1);
00503 xSyyy += sq(y1)*y1;
00504 xSxxx += sq(x1)*x1;
00505
00506 xSxxxy += sq(x1)*x1*y1;
00507 xSxxyy += sq(x1)*sq(y1);
00508 xSxyyy += x1*sq(y1)*y1;
00509
00510 xSyyyy += sq(y1)*sq(y1);
00511
00512 xSxxxx += sq(x1)*sq(x1);
00513
00514 xSu += uu[i];
00515 xSxu += x1*uu[i];
00516 xSyu += y1*uu[i];
00517 xSxxu += sq(x1)*uu[i];
00518 xSyyu += sq(y1)*uu[i];
00519 xSxyu += x1*y1*uu[i];
00520 }
00521
00522 ySx = ySxx = ySxy = ySy = ySyy = ySxxy = ySxyy = ySyyy = ySxxx =
00523 ySxxxy = ySxxyy = ySxyyy = ySyyyy = ySxxxx = ySv = ySxv = ySyv =
00524 ySxxv = ySyyv = ySxyv = 0.0;
00525 for (i = 0; i < nsy; i++) {
00526 x1 = xv[i]; y1 = yv[i];
00527
00528 ySx += x1; ySxx += sq(x1);
00529 ySxy += x1*y1;
00530 ySy += y1; ySyy += sq(y1);
00531
00532 ySxxy += sq(x1)*y1;
00533 ySxyy += x1*sq(y1);
00534 ySyyy += sq(y1)*y1;
00535 ySxxx += sq(x1)*x1;
00536
00537 ySxxxy += sq(x1)*x1*y1;
00538 ySxxyy += sq(x1)*sq(y1);
00539 ySxyyy += x1*sq(y1)*y1;
00540
00541 ySyyyy += sq(y1)*sq(y1);
00542
00543 ySxxxx += sq(x1)*sq(x1);
00544
00545 ySv += vv[i];
00546 ySxv += x1*vv[i];
00547 ySyv += y1*vv[i];
00548 ySxxv += sq(x1)*vv[i];
00549 ySyyv += sq(y1)*vv[i];
00550 ySxyv += x1*y1*vv[i];
00551 }
00552
00553 sys[0][0] = 2.*n1n2*xn; sys[0][1] = 2.*sqn*yn;
00554 sys[0][2] = n1n2*yn + sqn*xn; sys[0][3] = n1n2; sys[0][4] = sqn;
00555 sys[0][6] = 2.*sqn*xn; sys[0][7] = - 2.*n1n2*yn;
00556 sys[0][8] = sqn*yn - n1n2*xn; sys[0][9] = sqn; sys[0][10] = - n1n2;
00557
00558 sys[1][0] = xSxx; sys[1][1] = xSyy; sys[1][2] = xSxy;
00559 sys[1][3] = xSx; sys[1][4] = xSy; sys[1][5] = nsx;
00560 rhs[1] = xSu;
00561
00562 sys[2][0] = xSxxy; sys[2][1] = xSyyy; sys[2][2] = xSxyy;
00563 sys[2][3] = xSxy; sys[2][4] = xSyy; sys[2][5] = xSy;
00564 sys[2][12] = sqn/2.;
00565 rhs[2] = xSyu;
00566
00567 sys[3][0] = xSxxx; sys[3][1] = xSxyy; sys[3][2] = xSxxy;
00568 sys[3][3] = xSxx; sys[3][4] = xSxy; sys[3][5] = xSx;
00569 sys[3][12] = n1n2/2.;
00570 rhs[3] = xSxu;
00571
00572 sys[4][0] = xSxxxy; sys[4][1] = xSxyyy; sys[4][2] = xSxxyy;
00573 sys[4][3] = xSxxy; sys[4][4] = xSxyy; sys[4][5] = xSxy;
00574 sys[4][12] = (n1n2*yn + sqn*xn)/2.;
00575 rhs[4] = xSxyu;
00576
00577 sys[5][0] = xSxxyy; sys[5][1] = xSyyyy; sys[5][2] = xSxyyy;
00578 sys[5][3] = xSxyy; sys[5][4] = xSyyy; sys[5][5] = xSyy;
00579 sys[5][12] = sqn*yn;
00580 rhs[5] = xSyyu;
00581
00582 sys[6][0] = xSxxxx; sys[6][1] = xSxxyy; sys[6][2] = xSxxxy;
00583 sys[6][3] = xSxxx; sys[6][4] = xSxxy; sys[6][5] = xSxx;
00584 sys[6][12] = n1n2*xn;
00585 rhs[6] = xSxxu;
00586
00587
00588 sys[7][6] = ySxx; sys[7][7] = ySyy; sys[7][8] = ySxy;
00589 sys[7][9] = ySx; sys[7][10] = ySy; sys[7][11] = nsy;
00590 rhs[7] = ySv;
00591
00592 sys[8][6] = ySxxy; sys[8][7] = ySyyy; sys[8][8] = ySxyy;
00593 sys[8][9] = ySxy; sys[8][10] = ySyy; sys[8][11] = ySy;
00594 sys[8][12] = - n1n2/2.;
00595 rhs[8] = ySyv;
00596
00597 sys[9][6] = ySxxx; sys[9][7] = ySxyy; sys[9][8] = ySxxy;
00598 sys[9][9] = ySxx; sys[9][10] = ySxy; sys[9][11] = ySx;
00599 sys[9][12] = sqn/2.;
00600 rhs[9] = ySxv;
00601
00602 sys[10][6] = ySxxxy; sys[10][7] = ySxyyy; sys[10][8] = ySxxyy;
00603 sys[10][9] = ySxxy; sys[10][10] = ySxyy; sys[10][11] = ySxy;
00604 sys[10][12] = (- n1n2*xn + sqn*yn)/2.;
00605 rhs[10] = ySxyv;
00606
00607 sys[11][6] = ySxxyy; sys[11][7] = ySyyyy; sys[11][8] = ySxyyy;
00608 sys[11][9] = ySxyy; sys[11][10] = ySyyy; sys[11][11] = ySyy;
00609 sys[11][12] = - n1n2*yn;
00610 rhs[11] = ySyyv;
00611
00612 sys[12][6] = ySxxxx; sys[12][7] = ySxxyy; sys[12][8] = ySxxxy;
00613 sys[12][9] = ySxxx; sys[12][10] = ySxxy; sys[12][11] = ySxx;
00614 sys[12][12] = sqn*xn;
00615 rhs[12] = ySxxv;
00616
00617 #ifdef DEBUG_EXTRA
00618 for (i = 0; i < nsx; i++)
00619 fprintf(fptr1, "%g %g %f\n", xu[i], yu[i], uu[i]);
00620 for (i = 0; i < nsy; i++)
00621 fprintf(fptr2, "%g %g %f\n", xv[i], yv[i], vv[i]);
00622 fclose(fptr1); fclose(fptr2);
00623 #endif
00624
00625 if (gauss13(sys, rhs, sol, &det, 13)) {
00626 printf("Error: extra_velocity2(%g, %g): no solution\n", x, y);
00627 printf(" nsx: %d nsy: %d\n", nsx, nsy);
00628 for (i = 0; i < 13; i++) {
00629 for (j = 0; j < 13; j++)
00630 printf("%5g ", sys[i][j]);
00631 printf("\n");
00632 }
00633 exit(1);
00634 }
00635
00636 #ifdef DEBUG_EXTRA
00637 for (x1 = xminu; x1 <= xmaxu; x1 += 1.0) {
00638 for (y1 = yminu; y1 <= ymaxu; y1 += 1.0)
00639 fprintf(fptr3, "%g %g %g\n", x1, y1,
00640 sol[0]*sq(x1) + sol[1]*sq(y1) + sol[2]*x1*y1 +
00641 sol[3]*x1 + sol[4]*y1 + sol[5]);
00642 fprintf(fptr3, "\n");
00643 }
00644 fclose(fptr3);
00645 for (x1 = xminv; x1 <= xmaxv; x1 += 1.0) {
00646 for (y1 = yminv; y1 <= ymaxv; y1 += 1.0)
00647 fprintf(fptr4, "%g %g %g\n", x1, y1,
00648 sol[6]*sq(x1) + sol[7]*sq(y1) + sol[8]*x1*y1 +
00649 sol[9]*x1 + sol[10]*y1 + sol[11]);
00650 fprintf(fptr4, "\n");
00651 }
00652 fclose(fptr4);
00653 x1 = 0.0;
00654 for (i = 0; i < nsx; i++) {
00655 y1 = sol[0]*sq(xu[i]) + sol[1]*sq(yu[i]) + sol[2]*xu[i]*yu[i] +
00656 sol[3]*xu[i] + sol[4]*yu[i] + sol[5];
00657 x1 += sq(y1 - uu[i]);
00658 }
00659 for (i = 0; i < nsy; i++) {
00660 y1 = sol[6]*sq(xv[i]) + sol[7]*sq(yv[i]) + sol[8]*xv[i]*yv[i] +
00661 sol[9]*xv[i] + sol[10]*yv[i] + sol[11];
00662 x1 += sq(y1 - vv[i]);
00663 }
00664 printf("%g %g %g\n", x, y, x1);
00665 #endif
00666
00667 *eu = sol[0]*sq(x) + sol[1]*sq(y) + sol[2]*x*y +
00668 sol[3]*x + sol[4]*y + sol[5];
00669 *ev = sol[6]*sq(x) + sol[7]*sq(y) + sol[8]*x*y +
00670 sol[9]*x + sol[10]*y + sol[11];
00671 *dudx = 2.*sol[0]*x + sol[2]*y + sol[3];
00672 *dudy = 2.*sol[1]*y + sol[2]*x + sol[4];
00673 *dvdx = 2.*sol[6]*x + sol[8]*y + sol[9];
00674 *dvdy = 2.*sol[7]*y + sol[8]*x + sol[10];
00675 }
|
|
||||||||||||||||||||||||
|
Definition at line 705 of file extra.c. References real, splint, splint1, splint2, and sq.
00706 {
00707 real xp, yp, y, dl, kappa;
00708
00709 xp = splint1(px, t);
00710 yp = splint1(py, t);
00711 dl = sqrt(sq(xp) + sq(yp));
00712 kappa = (xp*splint2(py, t) - yp*splint2(px, t))/(sq(dl)*dl);
00713 /* axisymmetric term */
00714 y = splint(py, t) - 2.0;
00715 if (fabs(y) < 1e-2) /* second order approximation near the axis */
00716 kappa -= splint2(px, t)/splint1(py, t)/dl;
00717 else
00718 kappa -= xp/y/dl;
00719 *nx = - yp/dl;
00720 *ny = xp/dl;
00721 return kappa;
00722 }
|
1.2.18