surfaces.h File Reference

00048 {
00049   int i, j, k, l;
00050   int n, nin;
00051   real ti[14], t1, t2;
00052   polynom3 px, py;
00053   listinter listi[14], xl;
00054 
00055   if (cut->list == NULL) cut->list = (listinter *)malloc(sizeof(listinter));
00056   cut->ncut = 0;
00057   cut->xo = xo; cut->yo = yo;
00058 
00059   /* For each segment */
00060   for (l = 0; l < in.n - 1; l++) {
00061     /* the positions are shifted (-xo, -yo) */
00062     i = (real)(in.x[l] - xo); 
00063     j = (real)(in.y[l] - yo);
00064     t1 = in.t[l]; t2 = in.t[l+1];
00065     px = in.spx[l]; py = in.spy[l];
00066     px.a -= xo; py.a -= yo;
00067     
00068     /* first point */
00069     listi[0].t = t1; listi[0].type = MARK_TYPE;
00070     nin = 1;
00071     
00072     /* intersections with the horizontals */
00073     n = polyroots(py, t1, t2, ti, (real) j, INTERPREC);
00074     for (k = 0; k < n; k++) {
00075       listi[nin].t = ti[k]; listi[nin].type = HOR_TYPE; 
00076       listi[nin++].coord = j;
00077     }
00078     n = polyroots(py, t1, t2, ti, (real) j + 1.0, INTERPREC);
00079     for (k = 0; k < n; k++) {
00080       listi[nin].t = ti[k]; listi[nin].type = HOR_TYPE; 
00081       listi[nin++].coord = j + 1;
00082     }
00083     /* intersections with the verticals */
00084     n = polyroots(px, t1, t2, ti, (real) i, INTERPREC);
00085     for (k = 0; k < n; k++) {
00086       listi[nin].t = ti[k]; listi[nin].type = VER_TYPE; 
00087       listi[nin++].coord = i;
00088     }
00089     n = polyroots(px, t1, t2, ti, (real) i + 1.0, INTERPREC);
00090     for (k = 0; k < n; k++) {
00091       listi[nin].t = ti[k]; listi[nin].type = VER_TYPE; 
00092       listi[nin++].coord = i + 1;
00093     }
00094     
00095     /* last point */
00096     listi[nin].t = t2; listi[nin++].type = MARK_TYPE;
00097     
00098     /* sort the intersections in ascending order for t */
00099     for (k = 1; k < nin; k++) {
00100       xl = listi[k]; i = k - 1;
00101       while (i >= 0 && listi[i].t > xl.t) {
00102         listi[i + 1] = listi[i]; i--;
00103       }
00104       listi[i + 1] = xl;
00105     }
00106 
00107     /* add the points to the list */
00108     cut->list = (listinter *)realloc(cut->list,
00109                                      (cut->ncut + nin)*sizeof(listinter));
00110     for (k = 0; k < nin; k++) {
00111       listi[k].ip = l;
00112       cut->list[cut->ncut++] = listi[k];
00113     }
00114   }
00115 }

00119 {
00120   int i, ip;
00121   polynom3 px, py;
00122   real t;
00123   FILE *fptr = fopen(file, "wt");
00124 
00125   for (i = 0; i < cut.ncut; i++)
00126     if (cut.list[i].type == MARK_TYPE) {
00127       ip = cut.list[i].ip;
00128       px = in.spx[ip]; py = in.spy[ip];
00129       t = cut.list[i].t;
00130       fprintf(fptr, "%g %g\n", splint(px, t), splint(py, t));
00131     }
00132   fprintf(fptr, "&\n");
00133   for (i = 0; i < cut.ncut; i++)
00134     if (cut.list[i].type == HOR_TYPE) {
00135       ip = cut.list[i].ip;
00136       px = in.spx[ip]; py = in.spy[ip];
00137       t = cut.list[i].t;
00138       fprintf(fptr, "%g %g\n", splint(px, t), splint(py, t));
00139     }
00140   fprintf(fptr, "&\n");
00141   for (i = 0; i < cut.ncut; i++)
00142     if (cut.list[i].type == VER_TYPE) {
00143       ip = cut.list[i].ip;
00144       px = in.spx[ip]; py = in.spy[ip];
00145       t = cut.list[i].t;
00146       fprintf(fptr, "%g %g\n", splint(px, t), splint(py, t));
00147     }
00148   fclose(fptr);
00149 }

00010 {
00011   int i;
00012   real a, d, dmin = 1e30;
00013   for (i = 0; i < in.n - 1; i++) {
00014     d = sqrt(sq(x - 0.5*(in.x[i] + in.x[i+1])) + 
00015              sq(y - 0.5*(in.y[i] + in.y[i+1])));
00016     a = (in.x[i+1] - in.x[i])*(y - in.y[i]) - 
00017       (x - in.x[i])*(in.y[i+1] - in.y[i]);
00018     if (d < fabs(dmin)) dmin = a >= 0.0 ? d : -d;
00019   }
00020   if (dmin >= 0.0) return dmin;
00021   return 0.0;
00022 }

00360 {
00361   int i, j, l;
00362   real t, t1, t2;
00363   polynom3 px, py;
00364   real dz[3], rdz[6], rrdzp[9], irrdz[10];
00365 
00366   fill0(rrdz, nx, ny);
00367 
00368   /* For each intersection */
00369   for (l = 0; l < cut.ncut - 1; l++) {
00370     px = in.spx[cut.list[l].ip];
00371     py = in.spy[cut.list[l].ip];
00372     px.a -= cut.xo;
00373     py.a -= cut.yo;
00374 
00375     t1 = cut.list[l].t; t2 = cut.list[l+1].t;
00376     t = 0.5*(t1 + t2);
00377     i = splint(px, t);
00378     j = splint(py, t);
00379     
00380     py.a += cut.yo - 2.; px.a += cut.xo;
00381     polyderiv((real *)&px, 4, dz);
00382     polymul((real *)&py, 4, dz, 3, rdz);
00383     py.a -= cut.yo + (real)j - 1.5 + ra;
00384     polymul((real *)&py, 4, rdz, 6, rrdzp);
00385     polyint(rrdzp, 9, irrdz);
00386     
00387     /* add the circulation along the polynomial segment */
00388     rrdz[i][j] += polyeval(irrdz, 10, t2) - polyeval(irrdz, 10, t1);
00389   }
00390 }

00324 {
00325   int i, j, l;
00326   real t, t1, t2;
00327   polynom3 px, py;
00328   real dr[3], rdr[6], rzdrp[9], irzdr[10];
00329 
00330   fill0(rzdr, nx, ny);
00331 
00332   /* For each intersection */
00333   for (l = 0; l < cut.ncut - 1; l++) {
00334     px = in.spx[cut.list[l].ip];
00335     py = in.spy[cut.list[l].ip];
00336     px.a -= cut.xo;
00337     py.a -= cut.yo;
00338 
00339     t1 = cut.list[l].t; t2 = cut.list[l+1].t;
00340     t = 0.5*(t1 + t2);
00341     i = splint(px, t);
00342     j = splint(py, t);
00343 
00344     py.a += cut.yo - 2.; 
00345     polyderiv((real *)&py, 4, dr);
00346     polymul((real *)&py, 4, dr, 3, rdr);
00347     px.a -= (real)i + 0.5 + za;
00348     polymul((real *)&px, 4, rdr, 6, rzdrp);
00349     polyint(rzdrp, 9, irzdr);
00350     
00351     /* add the circulation along the polynomial segment */
00352     rzdr[i][j] += polyeval(irzdr, 10, t2) - polyeval(irzdr, 10, t1);
00353   }
00354 }

00156 {
00157   int i, j, l;
00158   real t, t1, t2;
00159   real x, y, cmax;
00160   polynom3 px, py;
00161   real dy[3], xdy[6], ixdy[7], xydy[9], ixydy[10];
00162 
00163   fill0(wa, nx, ny); fill0(wc, nx, ny);
00164   fill0(sx, nx, ny); fill0(sy, nx, ny);
00165 
00166   /* For each intersection */
00167   for (l = 0; l < cut.ncut - 1; l++) {
00168     px = in.spx[cut.list[l].ip];
00169     py = in.spy[cut.list[l].ip];
00170     px.a -= cut.xo;
00171     py.a -= cut.yo;
00172 
00173     t1 = cut.list[l].t; t2 = cut.list[l+1].t;
00174     t = 0.5*(t1 + t2);
00175     i = splint(px, t);
00176     j = splint(py, t);
00177     
00178     x = splint(px, t1);
00179     y = splint(py, t1);
00180 
00181     px.a -= (real)i;
00182     polyderiv((real *)&py, 4, dy);
00183     polymul((real *)&px, 4, dy, 3, xdy);
00184     polyint(xdy, 6, ixdy);
00185 
00186     py.a += cut.yo - 2.;
00187     polymul((real *)&py, 4, xdy, 6, xydy);
00188     polyint(xydy, 9, ixydy);
00189     
00190     /* Intersection with the verticals */
00191     if (cut.list[l].type == VER_TYPE) {
00192       /* deal with roundoff errors */
00193       if (y > (real)j + 1.0) y = (real)j + 1.0; 
00194       else if (y < (real)j) y = (real)j;
00195       /* add the circulation of the vertical segment */
00196       if (cut.list[l].coord == i) {
00197         wa[i-1][j] += (real)j + 1.0 - y;
00198         t = 0.5*(sq((real)j + cut.yo - 1.) - sq(y + cut.yo - 2.));
00199         wc[i-1][j] += t;
00200         sy[i][j] += t;
00201       } 
00202       else {
00203         wa[i][j] += y - (real)j;
00204         t = 0.5*(sq(y + cut.yo - 2.) - sq((real)j + cut.yo - 2.));
00205         wc[i][j] += t;
00206         sy[i+1][j] += t;
00207       }
00208     }
00209     /* Intersection with the horizontals */
00210     else if (cut.list[l].type == HOR_TYPE) {
00211       /* deal with roundoff errors */
00212       if (x > (real)i + 1.0) x = (real)i + 1.0; 
00213       else if (x < (real)i) x = (real)i;
00214       /* add the circulation of the horizontal segment */
00215       if (cut.list[l].coord == j)
00216         sx[i][j] += (x - (real)i)*((real)j + cut.yo - 2.);
00217       else
00218         sx[i][j+1] += ((real)i + 1. - x)*((real)j + cut.yo - 1.);
00219     }
00220     
00221     /* add the circulation along the polynomial segment */
00222     wa[i][j] += polyeval(ixdy, 7, t2) - polyeval(ixdy, 7, t1);
00223     wc[i][j] += polyeval(ixydy, 10, t2) - polyeval(ixydy, 10, t1);
00224   }
00225 
00226   /* merdic formula */
00227   for (i = 1; i <= nx; i++)
00228     for (j = 2; j <= ny; j++) {
00229       t = wa[i][j] = modf(wa[i][j], &y);
00230       if (wa[i][j] >= a[i][j] + 0.5) wa[i][j] -= 1.0;
00231       if (wa[i][j] <= a[i][j] - 0.5) wa[i][j] += 1.0;
00232       if (wa[i][j] <= a[i][j] - 0.5) wa[i][j] += 1.0;
00233 
00234       cmax = (real)j + cut.yo - 1.5;
00235       wc[i][j] = cmax == 0.0 ? 0.0 : cmax*modf(wc[i][j]/cmax, &y);
00236       if (wc[i][j] >= (a[i][j] + 0.5)*cmax) wc[i][j] -= cmax;
00237       if (wc[i][j] <= (a[i][j] - 0.5)*cmax) wc[i][j] += cmax;
00238       if (wc[i][j] <= (a[i][j] - 0.5)*cmax) wc[i][j] += cmax;
00239 
00240       if ((cut.yo != -0.5 || j != 2) && (wa[i][j] < - MAXROUND
00241                                          || wa[i][j] > 1.0 + MAXROUND))
00242         printf("WARNING wa[%d][%d]: %g %g %g xo: %g yo: %g\n", 
00243                i, j, wa[i][j], a[i][j], t, cut.xo, cut.yo);
00244       if ((cut.yo != -0.5 || j != 2) && (wc[i][j] < - MAXROUND*cmax 
00245                                          || wc[i][j] > cmax*(1. + MAXROUND)))
00246         printf("WARNING wc[%d][%d]: %g c[%d][%d]:%g cmax: %g xo,yo: %g,%g\n", i, j, wc[i][j], i, j, c[i][j], cmax, cut.xo, cut.yo);
00247       c[i][j] = wc[i][j];
00248       a[i][j] = wa[i][j];
00249     }
00250 }

00255 {
00256   int i, j;
00257   real sxmax;
00258 
00259   for (i = 2; i <= nx; i++)
00260     for (j = 3; j <= ny; j++) {
00261       sxmax = (real)j + yo - 2.;
00262       if (sx[i][j] > sxmax) sx[i][j] -= sxmax;
00263       if (sx[i][j] == 0.0 && (a[i][j] > 0.0 || a[i][j-1] > 0.0))
00264         if (a[i][j] == 1.0 || a[i][j-1] == 1.0 || 
00265             interface_inside(in, (real)i + xo + 0.5, (real)j + yo))
00266           sx[i][j] = sxmax;
00267       if (sx[i][j] < 0.0 || sx[i][j] > sxmax)
00268         printf("WARNING sx[%d][%d]:%g sxmax: %g xo,yo: %g,%g\n", 
00269                i, j, sx[i][j], sxmax, xo, yo);
00270     }
00271 }

00296 {
00297   int i, j;
00298   for (i = 1; i <= nx; i++)
00299     for (j = 3; j <= ny; j++)
00300       fprintf(fptr, "%g %g\n%g %g\n\n", 
00301               (real)i + 0.5 + xo - 0.5*sx[i][j]/((real)j + yo - 1.5), 
00302               (real)j + yo,
00303               (real)i + 0.5 + xo + 0.5*sx[i][j]/((real)j + yo - 1.5), 
00304               (real)j + yo);
00305 }

00276 {
00277   int i, j;
00278   real symax;
00279 
00280   for (i = 2; i <= nx; i++)
00281     for (j = 2; j <= ny; j++) {
00282       symax = (real)j + yo - 1.5;
00283       if (sy[i][j] > symax) sy[i][j] -= symax;
00284       if (sy[i][j] == 0.0 && (a[i][j] > 0.0 || a[i-1][j] > 0.0))
00285         if (a[i][j] == 1.0 || a[i-1][j] == 1.0 || 
00286             interface_inside(in, (real)i + xo, (real)j + yo + 0.5))
00287           sy[i][j] = symax;
00288       if (symax > 0.0 && (sy[i][j] < 0.0 || sy[i][j] > symax))
00289         printf("WARNING sy[%d][%d]:%g symax: %g xo,yo: %g,%g\n", 
00290                i, j, sy[i][j], symax, xo, yo);
00291     }
00292 }

00309 {
00310   int i, j;
00311   for (i = 1; i <= nx; i++)
00312     for (j = 2; j <= ny; j++)
00313       fprintf(fptr, "%g %g\n%g %g\n\n", 
00314               (real)i + xo, 
00315               (real)j + 0.5 + yo - 0.5*sy[i][j]/((real)j + yo - 1.5), 
00316               (real)i + xo,
00317               (real)j + 0.5 + yo + 0.5*sy[i][j]/((real)j + yo - 1.5));
00318 }


Data Structures
struct	interface_cut
struct	listinter
Functions
real	interface_inside (interface in, real x, real y)
void	cut_interface (interface in, interface_cut *cut, real xo, real yo)
void	interface_cut_write (interface in, interface_cut cut, char *file)
void	interface_surfaces (interface in, interface_cut cut, real2D sx, real2D sy, real2D a, real2D c, real2D wa, real2D wc, int nx, int ny)
void	sx_fill (real2D sx, real2D a, real xo, real yo, interface in, int nx, int ny)
void	sy_fill (real2D sy, real2D a, real xo, real yo, interface in, int nx, int ny)
void	sx_write (real2D sx, real xo, real yo, int nx, int ny, FILE *fptr)
void	sy_write (real2D sy, real xo, real yo, int nx, int ny, FILE *fptr)
void	interface_rzdr (interface in, interface_cut cut, real za, real2D rzdr, int nx, int ny)
void	interface_rrdz (interface in, interface_cut cut, real ra, real2D rrdz, int nx, int ny)

surfaces.h File Reference

Data Structures

Functions

Function Documentation