#include "utilf.h"
#include "markers.h"
#include "extra.h"
#include "make_bc.h"
#include "inout.h"
Include dependency graph for extra.c:
Go to the source code of this file.
Defines | |
#define | DIST 2.5 |
#define | SQR2 1.4142136 |
#define | SYMMETRY(val) (val < 2.0 ? 4. - val : val) |
#define | BC_YU(j, ny) (j <= 1 ? 3 - (j) : j >= ny ? 2*(ny) - (j) - 1 : j) |
#define | BC_YV(j, ny) (j < 2 ? 4 - (j) : j > ny ? 2*(ny) - (j) : j) |
#define | THETA1 0.463647609001 |
#define | THETA2 0.785398163398 |
#define | THETA3 1.10714871779 |
#define | THETA4 1.5707963268 |
#define | PREC 1e-3 |
#define | NCLOSEST 3 |
#define | SQCLOSEST sq(2*NCLOSEST + 1) |
#define | CLOSENESS(xu, yu, x, y, n1, n2) |
#define | GAUSSTOL 1.e-10 |
Functions | |
int | closest_normal (real x, real y, interface in, real *xmin, real *ymin, real *pmin, real *tmin, real *xn, real *yn) |
int | closest_u (real x, real y, real n1, real n2, real2D u, real2D ap, real *xu, real *yu, real *uu, real *xmax, real *xmin, real *ymax, real *ymin, int nx, int ny) |
int | closest_v (real x, real y, real n1, real n2, real2D v, real2D ap, real *xv, real *yv, real *vv, real *xmax, real *xmin, real *ymax, real *ymin, int nx, int ny) |
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 *dudx, real *dvdy, real *dudy, real *dvdx, 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 | gauss13 (real A[13][13], real B[13], real C[13], real *d, int n) |
int | gauss2 (real **A, real *B, real *C, int n) |
|
Definition at line 10 of file extra.c. Referenced by closest_u(). |
|
Definition at line 11 of file extra.c. Referenced by closest_v(). |
|
Value: (1.1*(sq(((xu) - (x))*(n1)) + \ sq(((yu) - (y))*(n2))) + \ 2.1*sq((n2)*((xu) - (x)) - n1*((yu) - (y)))) Definition at line 123 of file extra.c. Referenced by closest_u(), and closest_v(). |
|
|
|
|
|
Definition at line 121 of file extra.c. Referenced by closest_u(), and closest_v(). |
|
Definition at line 18 of file extra.c. Referenced by closest_normal(). |
|
Definition at line 122 of file extra.c. Referenced by extra_velocity_normals(), and extra_velocity_normals2(). |
|
|
|
|
|
|
|
|
|
|
|
|
|
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 130 of file extra.c. References BC_YU, CLOSENESS, INU, NCLOSEST, real, and real2D.
00134 { 00135 int i, j, k, i1, j1, j2, n = 0; 00136 real xl, yl, ul, dl; 00137 00138 *xmax = 0.0; *xmin = 2.*nx; *ymax = 0.0; *ymin = 2.*ny; 00139 00140 i = x; j = (int) ((real)y - 0.5); 00141 for (i1 = i - NCLOSEST; i1 <= i + NCLOSEST; i1++) 00142 for (j1 = j - NCLOSEST; j1 <= j + NCLOSEST; j1++) { 00143 j2 = BC_YU(j1,ny); 00144 if (INU(i1,j2)) { 00145 xu[n] = i1; yu[n] = j1 + 0.5; 00146 *xmax = xu[n] > *xmax ? xu[n] : *xmax; 00147 *ymax = yu[n] > *ymax ? yu[n] : *ymax; 00148 *xmin = xu[n] < *xmin ? xu[n] : *xmin; 00149 *ymin = yu[n] < *ymin ? yu[n] : *ymin; 00150 uu[n++] = u[i1][j2]; 00151 } 00152 } 00153 00154 /* sort the points */ 00155 for (k = 1; k < n; k++) { 00156 xl = xu[k]; yl = yu[k]; ul = uu[k]; 00157 i = k - 1; 00158 dl = CLOSENESS(xl, yl, x, y, n1, n2); 00159 while (i >= 0 && CLOSENESS(xu[i], yu[i], x, y, n1, n2) > dl) { 00160 xu[i+1] = xu[i]; yu[i+1] = yu[i]; uu[i+1] = uu[i]; 00161 i--; 00162 } 00163 xu[i+1] = xl; yu[i+1] = yl; uu[i+1] = ul; 00164 } 00165 return n; 00166 } |
|
Definition at line 169 of file extra.c. References BC_YV, CLOSENESS, INV, NCLOSEST, real, and real2D.
00173 { 00174 int i, j, k, i1, j1, j2, n = 0; 00175 real xl, yl, vl, dl; 00176 00177 *xmax = 0.0; *xmin = 2.*nx; *ymax = 0.0; *ymin = 2.*ny; 00178 00179 i = (int)((real)x - 0.5); j = y; 00180 for (i1 = i - NCLOSEST; i1 <= i + NCLOSEST; i1++) 00181 for (j1 = j - NCLOSEST; j1 <= j + NCLOSEST; j1++) { 00182 j2 = BC_YV(j1,ny); 00183 if (INV(i1,j2)) { 00184 xv[n] = i1 + 0.5; yv[n] = j1; 00185 *xmax = xv[n] > *xmax ? xv[n] : *xmax; 00186 *ymax = yv[n] > *ymax ? yv[n] : *ymax; 00187 *xmin = xv[n] < *xmin ? xv[n] : *xmin; 00188 *ymin = yv[n] < *ymin ? yv[n] : *ymin; 00189 if (j1 == j2) 00190 vv[n++] = v[i1][j1]; 00191 else 00192 vv[n++] = - v[i1][j2]; 00193 } 00194 } 00195 00196 /* sort the points */ 00197 for (k = 1; k < n; k++) { 00198 xl = xv[k]; yl = yv[k]; vl = vv[k]; 00199 i = k - 1; 00200 dl = CLOSENESS(xl, yl, x, y, n1, n2); 00201 while (i >= 0 && CLOSENESS(xv[i], yv[i], x, y, n1, n2) > dl) { 00202 xv[i+1] = xv[i]; yv[i+1] = yv[i]; vv[i+1] = vv[i]; 00203 i--; 00204 } 00205 xv[i+1] = xl; yv[i+1] = yl; vv[i+1] = vl; 00206 } 00207 return n; 00208 } |
|
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 802 of file extra.c. References real. Referenced by extra_velocity_normals2().
00803 { 00804 int i, j, k, i1; 00805 real max, temp, det = 1.0; 00806 00807 for (j = 0; j <= n - 2; j++) { 00808 i1 = j; 00809 max = fabs(A[j][j]); 00810 for (i = j + 1; i < n; i++) { 00811 if (fabs(A[i][j]) >= max) { 00812 i1 = i; 00813 max = fabs(A[i][j]); 00814 } 00815 } 00816 for (k = j; k < n; k++) { 00817 max = A[j][k]; A[j][k] = A[i1][k]; A[i1][k] = max; 00818 } 00819 max = B[j]; B[j] = B[i1]; B[i1] = max; 00820 for (i = j + 1; i < n; i++) { 00821 if (fabs(A[j][j]) < 1e-10) 00822 return 1; 00823 temp = A[i][j] / A[j][j]; 00824 for (k = j; k < n; k++) 00825 A[i][k] -= temp * A[j][k]; 00826 B[i] -= temp * B[j]; 00827 } 00828 } 00829 00830 for (i = 0; i < n; i++) 00831 det *= A[i][i]; 00832 *d = det; 00833 if (fabs(det) < 1.e-10) 00834 return 1; 00835 00836 if (fabs(A[n-1][n-1]) < 1.e-10) 00837 return 1; 00838 C[n - 1] = B[n - 1] / A[n - 1][n - 1]; 00839 for (i = n - 2; i >= 0; i--) { 00840 C[i] = B[i]; 00841 for (k = i + 1; k < n; k++) 00842 C[i] -= A[i][k] * C[k]; 00843 if (fabs(A[i][i]) < 1.e-10) 00844 return 1; 00845 C[i] /= A[i][i]; 00846 } 00847 return 0; 00848 } |
|
Definition at line 851 of file extra.c. References real.
00852 { 00853 int i, j, k, i1; 00854 real max, temp; 00855 00856 for (j = 0; j <= n - 2; j++) { 00857 i1 = j; 00858 max = fabs(A[j][j]); 00859 for (i = j + 1; i < n; i++) { 00860 if (fabs(A[i][j]) >= max) { 00861 i1 = i; 00862 max = fabs(A[i][j]); 00863 } 00864 } 00865 for (k = j; k < n; k++) { 00866 max = A[j][k]; A[j][k] = A[i1][k]; A[i1][k] = max; 00867 } 00868 max = B[j]; B[j] = B[i1]; B[i1] = max; 00869 for (i = j + 1; i < n; i++) { 00870 if (A[j][j] == 0.) 00871 return 1; 00872 else 00873 temp = A[i][j] / A[j][j]; 00874 for (k = j; k < n; k++) 00875 A[i][k] -= temp * A[j][k]; 00876 B[i] -= temp * B[j]; 00877 } 00878 } 00879 if (A[n - 1][n - 1] == 0.) 00880 return 1; 00881 else 00882 C[n - 1] = B[n - 1] / A[n - 1][n - 1]; 00883 for (i = n - 2; i >= 0; i--) { 00884 C[i] = B[i]; 00885 for (k = i + 1; k < n; k++) 00886 C[i] -= A[i][k] * C[k]; 00887 C[i] /= A[i][i]; 00888 } 00889 return 0; 00890 } |
|
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 } |