Main Page   Data Structures   File List   Data Fields   Globals  

surfaces.c File Reference

#include "utilf.h"
#include "markers.h"
#include "surfaces.h"

Include dependency graph for surfaces.c:

Include dependency graph

Go to the source code of this file.

Defines

#define MAXROUND   0.00001

Functions

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)

Define Documentation

#define MAXROUND   0.00001
 

Definition at line 5 of file surfaces.c.

Referenced by interface_surfaces().

Function Documentation

void cut_interface (  interface    in,
interface_cut *    cut,
real    xo,
real    yo
) 
 

Definition at line 8 of file surfaces.c.

References polynom3::a, listinter::coord, HOR_TYPE, INTERPREC, listinter::ip, interface_cut::list, MARK_TYPE, interface::n, interface_cut::ncut, polyroots(), real, interface::spx, interface::spy, listinter::t, interface::t, listinter::type, VER_TYPE, interface::x, interface_cut::xo, interface::y, and interface_cut::yo. 

00009 {
00010   int i, j, k, l;
00011   int n, nin;
00012   real ti[14], t1, t2;
00013   polynom3 px, py;
00014   listinter listi[14], xl;
00015 
00016   if (cut->list == NULL) cut->list = (listinter *)malloc(sizeof(listinter));
00017   cut->ncut = 0;
00018   cut->xo = xo; cut->yo = yo;
00019 
00020   /* For each segment */
00021   for (l = 0; l < in.n - 1; l++) {
00022     /* the positions are shifted (-xo, -yo) */
00023     i = (real)(in.x[l] - xo); 
00024     j = (real)(in.y[l] - yo);
00025     t1 = in.t[l]; t2 = in.t[l+1];
00026     px = in.spx[l]; py = in.spy[l];
00027     px.a -= xo; py.a -= yo;
00028     
00029     /* first point */
00030     listi[0].t = t1; listi[0].type = MARK_TYPE;
00031     nin = 1;
00032     
00033     /* intersections with the horizontals */
00034     n = polyroots(py, t1, t2, ti, (real) j, INTERPREC);
00035     for (k = 0; k < n; k++) {
00036       listi[nin].t = ti[k]; listi[nin].type = HOR_TYPE; 
00037       listi[nin++].coord = j;
00038     }
00039     n = polyroots(py, t1, t2, ti, (real) j + 1.0, INTERPREC);
00040     for (k = 0; k < n; k++) {
00041       listi[nin].t = ti[k]; listi[nin].type = HOR_TYPE; 
00042       listi[nin++].coord = j + 1;
00043     }
00044     /* intersections with the verticals */
00045     n = polyroots(px, t1, t2, ti, (real) i, INTERPREC);
00046     for (k = 0; k < n; k++) {
00047       listi[nin].t = ti[k]; listi[nin].type = VER_TYPE; 
00048       listi[nin++].coord = i;
00049     }
00050     n = polyroots(px, t1, t2, ti, (real) i + 1.0, INTERPREC);
00051     for (k = 0; k < n; k++) {
00052       listi[nin].t = ti[k]; listi[nin].type = VER_TYPE; 
00053       listi[nin++].coord = i + 1;
00054     }
00055     
00056     /* sort the intersections in ascending order for t */
00057     for (k = 1; k < nin; k++) {
00058       xl = listi[k]; i = k - 1;
00059       while (i >= 0 && listi[i].t > xl.t) {
00060         listi[i + 1] = listi[i]; i--;
00061       }
00062       listi[i + 1] = xl;
00063     }
00064 
00065     /* add the points to the list */
00066     cut->list = (listinter *)realloc(cut->list,
00067                                      (cut->ncut + nin)*sizeof(listinter));
00068     for (k = 0; k < nin; k++) {
00069       listi[k].ip = l;
00070       cut->list[cut->ncut++] = listi[k];
00071     }
00072   }
00073   /* last point */
00074   xl.t = in.t[in.n-1]; xl.type = MARK_TYPE; xl.ip = in.n - 1;
00075   cut->list = (listinter *)realloc(cut->list,
00076                                    (cut->ncut + 1)*sizeof(listinter));
00077   cut->list[cut->ncut++] = xl;
00078 }

void interface_cut_write (  interface    in,
interface_cut    cut,
char *    file
) 
 

Definition at line 81 of file surfaces.c.

References HOR_TYPE, listinter::ip, interface_cut::list, MARK_TYPE, interface_cut::ncut, real, splint, interface::spx, interface::spy, listinter::t, listinter::type, and VER_TYPE. 

00082 {
00083   int i, ip;
00084   polynom3 px, py;
00085   real t;
00086   FILE *fptr = fopen(file, "wt");
00087 
00088   for (i = 0; i < cut.ncut; i++)
00089     if (cut.list[i].type == MARK_TYPE) {
00090       ip = cut.list[i].ip;
00091       px = in.spx[ip]; py = in.spy[ip];
00092       t = cut.list[i].t;
00093       fprintf(fptr, "%g %g\n", splint(px, t), splint(py, t));
00094     }
00095   fprintf(fptr, "&\n");
00096   for (i = 0; i < cut.ncut; i++)
00097     if (cut.list[i].type == HOR_TYPE) {
00098       ip = cut.list[i].ip;
00099       px = in.spx[ip]; py = in.spy[ip];
00100       t = cut.list[i].t;
00101       fprintf(fptr, "%g %g\n", splint(px, t), splint(py, t));
00102     }
00103   fprintf(fptr, "&\n");
00104   for (i = 0; i < cut.ncut; i++)
00105     if (cut.list[i].type == VER_TYPE) {
00106       ip = cut.list[i].ip;
00107       px = in.spx[ip]; py = in.spy[ip];
00108       t = cut.list[i].t;
00109       fprintf(fptr, "%g %g\n", splint(px, t), splint(py, t));
00110     }
00111   fclose(fptr);
00112 }

void interface_rrdz (  interface    in,
interface_cut    cut,
real    ra,
real2D    rrdz,
int    nx,
int    ny
) 
 

Definition at line 349 of file surfaces.c.

References polynom3::a, fill0(), listinter::ip, interface_cut::list, interface_cut::ncut, real, real2D, splint, interface::spx, interface::spy, sq, listinter::t, interface_cut::xo, and interface_cut::yo. 

00352 {
00353   int i, j, l;
00354   real t, t1, t2;
00355   polynom3 px, py;
00356   real x1, y1, x2, y2;
00357 
00358   fill0(rrdz, nx, ny);
00359 
00360   /* For each intersection */
00361   for (l = 0; l < cut.ncut - 1; l++) {
00362     px = in.spx[cut.list[l].ip];
00363     py = in.spy[cut.list[l].ip];
00364     px.a -= cut.xo;
00365     py.a -= cut.yo;
00366 
00367     t1 = cut.list[l].t; t2 = cut.list[l+1].t;
00368     t = 0.5*(t1 + t2);
00369     i = splint(px, t);
00370     j = splint(py, t);
00371     
00372     /* add the circulation along the polynomial segment */
00373     py.a += cut.yo - 2.;
00374     x1 = splint(px, t1); x2 = splint(px, t2);
00375     y1 = splint(py, t1); y2 = splint(py, t2);
00376     t = cut.yo + (real)j - 1.5 + ra;
00377     rrdz[i][j] +=
00378       (x2 - x1)*(sq(y2 - y1)/3. + (2.*y1 - t)*(y2 - y1)/2. + y1*(y1 - t));
00379   }
00380 }

void interface_rzdr (  interface    in,
interface_cut    cut,
real    za,
real2D    rzdr,
int    nx,
int    ny
) 
 

Definition at line 316 of file surfaces.c.

References polynom3::a, fill0(), listinter::ip, interface_cut::list, interface_cut::ncut, real, real2D, splint, interface::spx, interface::spy, listinter::t, interface_cut::xo, and interface_cut::yo. 

00319 {
00320   int i, j, l;
00321   real t, t1, t2;
00322   polynom3 px, py;
00323   real x1, y1, x2, y2;
00324 
00325   fill0(rzdr, nx, ny);
00326   /* For each intersection */
00327   for (l = 0; l < cut.ncut - 1; l++) {
00328     px = in.spx[cut.list[l].ip];
00329     py = in.spy[cut.list[l].ip];
00330     px.a -= cut.xo;
00331     py.a -= cut.yo;
00332 
00333     t1 = cut.list[l].t; t2 = cut.list[l+1].t;
00334     t = 0.5*(t1 + t2);
00335     i = splint(px, t);
00336     j = splint(py, t);
00337 
00338     /* add the circulation along the interface segment */
00339     py.a += cut.yo - 2.; 
00340     px.a -= (real)i + 0.5 + za;
00341     x1 = splint(px, t1); x2 = splint(px, t2);
00342     y1 = splint(py, t1); y2 = splint(py, t2);
00343     rzdr[i][j] += (y2 - y1)*((x2 - x1)*(y2 - y1)/3. + 
00344                              ((x2 - x1)*y1 + (y2 - y1)*x1)/2. + x1*y1);
00345   }
00346 }

void interface_surfaces (  interface    in,
interface_cut    cut,
real2D    sx,
real2D    sy,
real2D    a,
real2D    c,
real2D    wa,
real2D    wc,
int    nx,
int    ny
) 
 

Definition at line 115 of file surfaces.c.

References polynom3::a, listinter::coord, fill0(), HOR_TYPE, listinter::ip, interface_cut::list, MAXROUND, interface_cut::ncut, real, real2D, splint, interface::spx, interface::spy, sq, listinter::t, listinter::type, VER_TYPE, interface_cut::xo, and interface_cut::yo. 

00119 {
00120   int i, j, l;
00121   real t, t1, t2;
00122   real x, y, cmax;
00123   polynom3 px, py;
00124   real x1, x2, y1, y2;
00125 
00126   fill0(wa, nx, ny); fill0(wc, nx, ny);
00127   fill0(sx, nx, ny); fill0(sy, nx, ny);
00128 
00129   /* For each intersection */
00130   for (l = 0; l < cut.ncut - 1; l++) {
00131     px = in.spx[cut.list[l].ip];
00132     py = in.spy[cut.list[l].ip];
00133     px.a -= cut.xo;
00134     py.a -= cut.yo;
00135 
00136     t1 = cut.list[l].t; t2 = cut.list[l+1].t;
00137     t = 0.5*(t1 + t2);
00138     i = splint(px, t);
00139     j = splint(py, t);
00140 
00141     x = splint(px, t1);
00142     y = splint(py, t1);
00143 
00144     /* add the circulations along the interface segment */
00145     px.a -= (real)i; py.a += cut.yo - 2.;
00146     x1 = splint(px, t1); x2 = splint(px, t2);
00147     y1 = splint(py, t1); y2 = splint(py, t2);
00148 
00149     wa[i][j] += (y2 - y1)*(x1 + x2)/2.;
00150     wc[i][j] += (y2 - y1)*((x2 - x1)*(y2 - y1)/3. + 
00151                            ((x2 - x1)*y1 + (y2 - y1)*x1)/2. + x1*y1);
00152     
00153     /* Intersection with the verticals */
00154     if (cut.list[l].type == VER_TYPE) {
00155       /* deal with roundoff errors */
00156       if (y > (real)j + 1.0) y = (real)j + 1.0; 
00157       else if (y < (real)j) y = (real)j;
00158       /* add the circulation of the vertical segment */
00159       if (cut.list[l].coord == i) {
00160         wa[i-1][j] += (real)j + 1.0 - y;
00161         t = 0.5*(sq((real)j + cut.yo - 1.) - sq(y + cut.yo - 2.));
00162         wc[i-1][j] += t;
00163         sy[i][j] += t;
00164       } 
00165       else {
00166         wa[i][j] += y - (real)j;
00167         t = 0.5*(sq(y + cut.yo - 2.) - sq((real)j + cut.yo - 2.));
00168         wc[i][j] += t;
00169         sy[i+1][j] += t;
00170       }
00171     }
00172     /* Intersection with the horizontals */
00173     else if (cut.list[l].type == HOR_TYPE) {
00174       /* deal with roundoff errors */
00175       if (x > (real)i + 1.0) x = (real)i + 1.0; 
00176       else if (x < (real)i) x = (real)i;
00177       /* add the circulation of the horizontal segment */
00178       if (cut.list[l].coord == j)
00179         sx[i][j] += (x - (real)i)*((real)j + cut.yo - 2.);
00180       else
00181         sx[i][j+1] += ((real)i + 1. - x)*((real)j + cut.yo - 1.);
00182     }
00183   }
00184 
00185 
00186 #if 0
00187   /* merdic formula */
00188   for (i = 1; i <= nx; i++)
00189     for (j = 2; j <= ny; j++) {
00190       t = wa[i][j] = modf(wa[i][j], &y);
00191       if (wa[i][j] >= a[i][j] + 0.5 && wa[i][j]>=1.0) wa[i][j] -= 1.0;
00192       if (wa[i][j] <= a[i][j] - 0.5 && wa[i][j]<=0.0) wa[i][j] += 1.0;
00193       if (wa[i][j] <= a[i][j] - 0.5 && wa[i][j]<=0.0) wa[i][j] += 1.0;
00194 
00195       cmax = (real)j + cut.yo - 1.5;
00196       wc[i][j] = cmax == 0.0 ? 0.0 : cmax*modf(wc[i][j]/cmax, &y);
00197       if (wc[i][j] >= (a[i][j] + 0.5)*cmax && wc[i][j]>=cmax) wc[i][j] -= cmax;
00198       if (wc[i][j] <= (a[i][j] - 0.5)*cmax && wc[i][j]<=0.0) wc[i][j] += cmax;
00199       if (wc[i][j] <= (a[i][j] - 0.5)*cmax && wc[i][j]<=0.0) wc[i][j] += cmax;
00200 
00201       if ((cut.yo != -0.5 || (j != 2 && j != ny)) && (wa[i][j] < - MAXROUND
00202                                          || wa[i][j] > 1.0 + MAXROUND))
00203         printf("WARNING wa[%d][%d]: %g %g %g xo: %g yo: %g\n", 
00204                i, j, wa[i][j], a[i][j], t, cut.xo, cut.yo);
00205       if ((cut.yo != -0.5 || (j != 2 && j != ny)) && (wc[i][j] < 
00206                                                       - MAXROUND*cmax 
00207                                          || wc[i][j] > cmax*(1. + MAXROUND)))
00208         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);
00209       c[i][j] = wc[i][j] > cmax ? cmax : wc[i][j] < 0.0 ? 0.0 : wc[i][j];
00210       a[i][j] = wa[i][j] > 1.0 ? 1.0 : wa[i][j] < 0.0 ? 0.0 : wa[i][j];
00211     }
00212 
00213 #else
00214 
00215   /* merdic formula */
00216   for (i = 1; i <= nx; i++)
00217     for (j = 2; j <= ny; j++) {
00218       t = wa[i][j] = modf(wa[i][j], &y);
00219       if (wa[i][j] >= a[i][j] + 0.5) wa[i][j] -= 1.0;
00220       if (wa[i][j] <= a[i][j] - 0.5) wa[i][j] += 1.0;
00221       if (wa[i][j] <= a[i][j] - 0.5) wa[i][j] += 1.0;
00222 
00223 /*      wa[i][j] = wa[i][j] > 1.0 ? 1.0 : wa[i][j] < 0.0 ? 0.0 : wa[i][j]; */
00224 
00225       cmax = (real)j + cut.yo - 1.5;
00226       wc[i][j] = cmax == 0.0 ? 0.0 : cmax*modf(wc[i][j]/cmax, &y);
00227       if (wc[i][j] >= (a[i][j] + 0.5)*cmax) wc[i][j] -= cmax;
00228       if (wc[i][j] <= (a[i][j] - 0.5)*cmax) wc[i][j] += cmax;
00229       if (wc[i][j] <= (a[i][j] - 0.5)*cmax) wc[i][j] += cmax;
00230 
00231 /*      wc[i][j] = wc[i][j] > cmax ? cmax : wc[i][j] < 0.0 ? 0.0 : wc[i][j]; */
00232 
00233       if ((cut.yo != -0.5 || (j != 2 && j != ny)) && (wa[i][j] < - MAXROUND
00234                                          || wa[i][j] > 1.0 + MAXROUND))
00235         printf("WARNING wa[%d][%d]: %g %g %g xo: %g yo: %g\n",
00236                i, j, wa[i][j], a[i][j], t, cut.xo, cut.yo);
00237       if ((cut.yo != -0.5 || (j != 2 && j != ny)) && (wc[i][j] <
00238                                                       - MAXROUND*cmax
00239                                          || wc[i][j] > cmax*(1. + MAXROUND)))
00240         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);
00241       c[i][j] = wc[i][j] > cmax ? cmax : wc[i][j] < 0.0 ? 0.0 : wc[i][j];
00242       a[i][j] = wa[i][j] > 1.0 ? 1.0 : wa[i][j] < 0.0 ? 0.0 : wa[i][j];
00243     }
00244 #endif
00245 }

void sx_fill (  real2D    sx,
real2D    a,
real    xo,
real    yo,
interface    in,
int    nx,
int    ny
) 
 

Definition at line 248 of file surfaces.c.

References interface_inside(), real, and real2D. 

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

void sx_write (  real2D    sx,
real    xo,
real    yo,
int    nx,
int    ny,
FILE *    fptr
) 
 

Definition at line 290 of file surfaces.c.

References real, and real2D. 

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

void sy_fill (  real2D    sy,
real2D    a,
real    xo,
real    yo,
interface    in,
int    nx,
int    ny
) 
 

Definition at line 269 of file surfaces.c.

References interface_inside(), real, and real2D. 

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

void sy_write (  real2D    sy,
real    xo,
real    yo,
int    nx,
int    ny,
FILE *    fptr
) 
 

Definition at line 303 of file surfaces.c.

References real, and real2D. 

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

Generated on Wed Feb 19 22:28:40 2003 for Markers by doxygen1.2.18