This graph shows which files directly or indirectly include this file:
Go to the source code of this file.
Functions | |
| void | interface_fluxes (interface in, real2D cc, real2D s1, real2D s2, real2D s3, real2D s4, int nx, int ny) |
| void | interface_hfrac (interface in, real2D fx, real2D c, real xo, real yo, int nx, int ny) |
| void | interface_vfrac (interface in, real2D fy, real2D c, real xo, real yo, int nx, int ny) |
|
||||||||||||||||||||||||||||||||||||
|
Definition at line 7 of file fluxes.c. References listinter::coord, fill0(), HOR_TYPE, INP, INTERPREC, MARK_TYPE, interface::n, polyroots(), real, real2D, splint, interface::spx, interface::spy, sq, listinter::t, interface::t, listinter::type, VER_TYPE, interface::x, and interface::y.
00010 {
00011 int i, j, k, l;
00012 int n, nin;
00013 real ti[14], t, t1, t2;
00014 real x, y, x1, y1, yo;
00015 polynom3 px, py;
00016 typedef struct {
00017 real t;
00018 int type, coord;
00019 } listinter;
00020 listinter listi[14], xl;
00021
00022 fill0(s1, nx, ny); fill0(s3, nx, ny);
00023 fill0(s2, nx, ny); fill0(s4, nx, ny);
00024 #if 0
00025 /* For each segment */
00026 for (l = 0; l < in.n - 1; l++) {
00027 i = (real)in.x[l];
00028 j = (real)in.y[l];
00029 t1 = in.t[l]; t2 = in.t[l+1];
00030 px = in.spx[l]; py = in.spy[l];
00031
00032 /* first point */
00033 listi[0].t = t1; listi[0].type = MARK_TYPE;
00034 nin = 1;
00035
00036 /* intersections with the horizontals */
00037 n = polyroots(py, t1, t2, ti, (real) j, INTERPREC);
00038 for (k = 0; k < n; k++) {
00039 listi[nin].t = ti[k]; listi[nin].type = HOR_TYPE;
00040 listi[nin++].coord = j;
00041 }
00042 n = polyroots(py, t1, t2, ti, (real) j + 1.0, INTERPREC);
00043 for (k = 0; k < n; k++) {
00044 listi[nin].t = ti[k]; listi[nin].type = HOR_TYPE;
00045 listi[nin++].coord = j + 1;
00046 }
00047 /* intersections with the verticals */
00048 n = polyroots(px, t1, t2, ti, (real) i, INTERPREC);
00049 for (k = 0; k < n; k++) {
00050 listi[nin].t = ti[k]; listi[nin].type = VER_TYPE;
00051 listi[nin++].coord = i;
00052 }
00053 n = polyroots(px, t1, t2, ti, (real) i + 1.0, INTERPREC);
00054 for (k = 0; k < n; k++) {
00055 listi[nin].t = ti[k]; listi[nin].type = VER_TYPE;
00056 listi[nin++].coord = i + 1;
00057 }
00058
00059 /* last point */
00060 listi[nin].t = t2; listi[nin++].type = MARK_TYPE;
00061
00062 /* sort the intersections in ascending order for t */
00063 for (k = 1; k < nin; k++) {
00064 xl = listi[k]; i = k - 1;
00065 while (i >= 0 && listi[i].t > xl.t) {
00066 listi[i + 1] = listi[i]; i--;
00067 }
00068 listi[i + 1] = xl;
00069 }
00070
00071 for (k = 0; k < nin - 1; k++) {
00072 t1 = listi[k].t; t2 = listi[k+1].t;
00073 t = 0.5*(t1 + t2);
00074 i = splint(px, t);
00075 j = splint(py, t);
00076
00077 x = splint(px, t1); x1 = splint(px, t2);
00078 y = splint(py, t1) - 2.; y1 = splint(py, t2) - 2.;
00079 yo = (real)j - 2.;
00080
00081 /* Intersection with the verticals */
00082 if (listi[k].type == VER_TYPE) {
00083 /* add the circulation of the vertical segment */
00084 if (listi[k].coord == i) {
00085 t = (1. + yo - y)*(1. + yo + y)/2.;
00086 s3[i][j] -= t;
00087 s1[i-1][j] += t;
00088 }
00089 else {
00090 t = (y - yo)*(y + yo)/2.;
00091 s1[i][j] += t;
00092 s3[i+1][j] -= t;
00093 }
00094 }
00095
00096 /* Intersection with the horizontals */
00097 if (listi[k].type == HOR_TYPE) {
00098 /* add the circulation of the horizontal segment */
00099 if (listi[k].coord == j) {
00100 t = (x - i)*yo;
00101 s4[i][j] -= t;
00102 s2[i][j-1] += t;
00103 }
00104 else {
00105 t = (1. + i - x)*(yo + 1.);
00106 s2[i][j] += t;
00107 s4[i][j+1] -= t;
00108 }
00109 }
00110
00111 /* add the circulation along the segment */
00112 x -= i; x1 -= i;
00113 t = (y1 - y)*((x1*y1 + x*y)/3. + (x1*y + x*y1)/6.);
00114 s1[i][j] += t;
00115 s3[i][j] -= t - (y1 - y)*(y1 + y)/2.;
00116 t = (sq(y1) + y1*y + sq(y))/3.;
00117 s2[i][j] += (x1 - x)*(yo*(y1 + y)/2. - t);
00118 s4[i][j] -= (x1 - x)*((yo + 1.)*(y1 + y)/2. - t);
00119 }
00120 }
00121 #endif
00122 for (i = 2; i < nx; i++)
00123 for (j = 2; j < ny; j++)
00124 if (INP(i,j)) {
00125 yo = (real)j - 2.;
00126 if (s1[i][j] <= 0.1*yo) s1[i][j] += yo + 0.5;
00127 else if (s1[i][j] > yo + 0.5) s1[i][j] -= yo + 0.5;
00128 if (s3[i][j] >= - 0.1*yo) s3[i][j] -= yo + 0.5;
00129 else if (s3[i][j] < - yo - 0.5) s3[i][j] += yo + 0.5;
00130 if (s2[i][j] <= 0.1*yo) s2[i][j] += yo + 1.;
00131 else if (s2[i][j] > yo + 1.) s2[i][j] -= yo + 1.;
00132 if (s4[i][j] >= - 0.1*yo) s4[i][j] -= yo;
00133 else if (s4[i][j] < -yo) s4[i][j] += yo;
00134
00135 if (s1[i][j] < 0.1*yo) {
00136 s1[i][j] += yo + 0.5;
00137 printf("s1[%d][%d]=%g <= 0.1\n", i, j, s1[i][j]);
00138 }
00139 else if (s1[i][j] > yo + 0.5) {
00140 s1[i][j] -= yo + 0.5;
00141 printf("s1[%d][%d]=%g >\n", i, j, s1[i][j]);
00142 }
00143 if (s3[i][j] > - 0.1*yo) {
00144 s3[i][j] -= yo + 0.5;
00145 printf("s3[%d][%d]=%g >= - 0.1\n", i, j, s3[i][j]);
00146 }
00147 else if (s3[i][j] < - yo - 0.5) {
00148 s3[i][j] += yo + 0.5;
00149 printf("s3[%d][%d]=%g <= \n", i, j, s3[i][j]);
00150 }
00151 if (s2[i][j] < 0.1*yo) {
00152 s2[i][j] += yo + 1.;
00153 printf("s2[%d][%d]=%g <= 0.1\n", i, j, s2[i][j]);
00154 }
00155 else if (s2[i][j] > yo + 1.) {
00156 s2[i][j] -= yo + 1.;
00157 printf("s2[%d][%d]=%g >\n", i, j, s2[i][j]);
00158 }
00159 if (s4[i][j] > - 0.1*yo) {
00160 s4[i][j] -= yo;
00161 printf("s4[%d][%d]=%g >= - 0.1\n", i, j, s4[i][j]);
00162 }
00163 else if (s4[i][j] < -yo) {
00164 s4[i][j] += yo;
00165 printf("s4[%d][%d]=%g <\n", i, j, s4[i][j]);
00166 }
00167 /* printf("%d,%d: s1: %g s2: %g s3: %g s4: %g\n", i, j,
00168 s1[i][j], s2[i][j], s3[i][j], s4[i][j]); */
00169 if (fabs(2.*(s1[i][j] + s3[i][j])/(s1[i][j] - s3[i][j])) > 1e-3) {
00170 fprintf(stderr, "interface_fluxes(): s1[%d][%d] = %g != - s3[%d][%d] = %g\n",
00171 i, j, s1[i][j], i, j, - s3[i][j]);
00172 exit(1);
00173 }
00174 }
00175 }
|
|
||||||||||||||||||||||||||||||||
|
Definition at line 178 of file fluxes.c. References polynom3::a, fill0(), INTERPREC, interface::n, polyroots(), real, real2D, splint, interface::spx, interface::spy, interface::t, interface::x, and interface::y.
00181 {
00182 int i, j, k, l;
00183 int n;
00184 real ti[14], t, t1, t2;
00185 real x;
00186 polynom3 px, py;
00187 FILE *fptr = fopen("hfrac", "wt");
00188
00189 fill0(fx, nx, ny);
00190
00191 /* For each segment */
00192 for (l = 0; l < in.n - 1; l++) {
00193 i = (real)in.x[l] - xo;
00194 j = (real)in.y[l] - yo;
00195 t1 = in.t[l]; t2 = in.t[l+1];
00196 px = in.spx[l]; py = in.spy[l];
00197 px.a -= xo; py.a -= yo;
00198
00199 /* intersections with the horizontals */
00200 n = polyroots(py, t1, t2, ti, (real) j, INTERPREC);
00201 for (k = 0; k < n; k++) {
00202 x = splint(px, ti[k]); i = x;
00203 fx[i][j] += (real)i + 1. - x;
00204 }
00205 n = polyroots(py, t1, t2, ti, (real) j + 1.0, INTERPREC);
00206 for (k = 0; k < n; k++) {
00207 x = splint(px, ti[k]); i = x;
00208 fx[i][j+1] += x - (real)i;
00209 }
00210 }
00211
00212 for (i = 2; i < nx; i++)
00213 for (j = 2; j < ny; j++) {
00214 if (fx[i][j] > 1.0) fx[i][j] -= 1.0;
00215 if (fx[i][j] == 0.0 && c[i][j] + c[i][j-1] >= 1.0) fx[i][j] = 1.0;
00216 if (fx[i][j] > 1.0 || fx[i][j] < 0.0) {
00217 fprintf(stderr, "interface_hfrac(): fx[%d][%d]=%g\n", i, j, fx[i][j]);
00218 exit(1);
00219 }
00220 fprintf(fptr, "%g %g\n%g %g\n\n",
00221 (real)i + xo + 0.5 + fx[i][j]/2., (real)j + yo,
00222 (real)i + yo + 0.5 - fx[i][j]/2., (real)j + yo);
00223 }
00224 fclose(fptr);
00225 }
|
|
||||||||||||||||||||||||||||||||
|
Definition at line 228 of file fluxes.c. References polynom3::a, fill0(), INTERPREC, interface::n, polyroots(), real, real2D, splint, interface::spx, interface::spy, sq, interface::t, interface::x, and interface::y.
00230 {
00231 int i, j, k, l;
00232 int n;
00233 real ti[14], t, t1, t2;
00234 real y;
00235 polynom3 px, py;
00236 FILE *fptr = fopen("vfrac", "wt");
00237
00238 fill0(fy, nx, ny);
00239
00240 /* For each segment */
00241 for (l = 0; l < in.n - 1; l++) {
00242 i = (real)in.x[l] - xo;
00243 j = (real)in.y[l] - yo;
00244 t1 = in.t[l]; t2 = in.t[l+1];
00245 px = in.spx[l]; py = in.spy[l];
00246 px.a -= xo; py.a -= yo;
00247
00248 /* intersections with the verticals */
00249 n = polyroots(px, t1, t2, ti, (real) i, INTERPREC);
00250 for (k = 0; k < n; k++) {
00251 y = splint(py, ti[k]); j = y;
00252 fy[i][j] += 0.5*(sq(y + yo - 2.) - sq((real)j + yo - 2.))/
00253 ((real)j + yo - 1.5);
00254 }
00255 n = polyroots(px, t1, t2, ti, (real) i + 1.0, INTERPREC);
00256 for (k = 0; k < n; k++) {
00257 y = splint(py, ti[k]); j = y;
00258 fy[i+1][j] += 0.5*(sq((real)j + yo - 1.) - sq(y + yo - 2.))/
00259 ((real)j + yo - 1.5);
00260 }
00261 }
00262
00263 for (i = 2; i < nx; i++)
00264 for (j = 2; j < ny; j++) {
00265 if (fy[i][j] > 1.0) fy[i][j] -= 1.0;
00266 if (fy[i][j] == 0.0 && c[i][j] + c[i-1][j] >= 1.0) fy[i][j] = 1.0;
00267 if (fy[i][j] > 1.0 || fy[i][j] < 0.0) {
00268 fprintf(stderr, "interface_vfrac(): fy[%d][%d]=%g\n", i, j, fy[i][j]);
00269 exit(1);
00270 }
00271 fprintf(fptr, "%g %g\n%g %g\n\n",
00272 (real)i + xo, (real)j + yo + 0.5 + fy[i][j]/2.,
00273 (real)i + yo, (real)j + yo + 0.5 - fy[i][j]/2.);
00274 }
00275 fclose(fptr);
00276 }
|
1.2.18