This graph shows which files directly or indirectly include this file:
Go to the source code of this file.
Functions | |
| int | closest_u (real x, real y, real2D u, real2D cu, real *xu, real *yu, real *uu, int nx, int ny) |
| int | closest_v (real x, real y, real2D v, real2D cv, real *xv, real *yv, real *vv, int nx, int ny) |
| void | extra_poly (real x, real y, real2D u, real2D cu, real2D v, real2D cv, interface in, real *eu, real *ev, int nx, int ny) |
|
||||||||||||||||||||||||||||||||||||||||
|
Definition at line 39 of file extra1.c. References INU, real, real2D, and sq. Referenced by extra_poly(), extra_velocity_normals(), and extra_velocity_normals2().
00041 {
00042 int i, j, k, i1, j1, n = 0;
00043 real xl, yl, ul, dl;
00044
00045 i = x; j = (int) ((real)y - 0.5);
00046 for (i1 = i - 2; i1 <= i + 2; i1++)
00047 for (j1 = j - 2; j1 <= j + 2; j1++)
00048 if (INU(i1,j1)) {
00049 xu[n] = i1; yu[n] = j1 + 0.5; uu[n++] = u[i1][j1];
00050 }
00051
00052 /* sort the points */
00053 for (k = 1; k < n; k++) {
00054 xl = xu[k]; yl = yu[k]; ul = uu[k];
00055 i = k - 1;
00056 dl = sq(x - xl) + sq(y - yl);
00057 while (i >= 0 && sq(xu[i] - x) + sq(yu[i] - y) > dl) {
00058 xu[i+1] = xu[i]; yu[i+1] = yu[i]; uu[i+1] = uu[i];
00059 i--;
00060 }
00061 xu[i+1] = xl; yu[i+1] = yl; uu[i+1] = ul;
00062 }
00063 return n;
00064 }
|
|
||||||||||||||||||||||||||||||||||||||||
|
Definition at line 67 of file extra1.c. References INV, real, real2D, and sq. Referenced by extra_poly(), extra_velocity_normals(), and extra_velocity_normals2().
00069 {
00070 int i, j, k, i1, j1, n = 0;
00071 real xl, yl, vl, dl;
00072
00073 i = (int)((real)x - 0.5); j = y;
00074 for (i1 = i - 2; i1 <= i + 2; i1++)
00075 for (j1 = j - 2; j1 <= j + 2; j1++)
00076 if (INV(i1,j1)) {
00077 xv[n] = i1 + 0.5; yv[n] = j1; vv[n++] = v[i1][j1];
00078 }
00079
00080 /* sort the points */
00081 for (k = 1; k < n; k++) {
00082 xl = xv[k]; yl = yv[k]; vl = vv[k];
00083 i = k - 1;
00084 dl = sq(x - xl) + sq(y - yl);
00085 while (i >= 0 && sq(xv[i] - x) + sq(yv[i] - y) > dl) {
00086 xv[i+1] = xv[i]; yv[i+1] = yv[i]; vv[i+1] = vv[i];
00087 i--;
00088 }
00089 xv[i+1] = xl; yv[i+1] = yl; vv[i+1] = vl;
00090 }
00091 return n;
00092 }
|
|
||||||||||||||||||||||||||||||||||||||||||||||||
|
Definition at line 140 of file extra1.c. References closest_normal(), closest_u(), closest_v(), four_aligned(), gauss3(), interface_arc(), interface::n, NI, NU, NV, real, real2D, splint, splint1, interface::spx, interface::spy, sq, interface::t, xi, and yi.
00144 {
00145 real xu[25], yu[25], uu[25];
00146 real xv[25], yv[25], vv[25];
00147 real xi[NI], yi[NI], xni[NI], yni[NI];
00148 real xmin, ymin, pmin, tmin, xn, yn, t, t1, dl;
00149 int i, j, n;
00150 real a[12][12], b[12], c[12], tt, tp, det;
00151 FILE *fptr;
00152 static int graph, graph1;
00153
00154 if (!closest_normal(x, y, in, &xmin, &ymin, &pmin, &tmin, &xn, &yn)) {
00155 fprintf(stderr, "extra_poly(%g,%g): no closest normal\n", x, y);
00156 exit(1);
00157 }
00158 for (n = 0, t = tmin - 0.5; n < NI; n++, t += 1.0)
00159 if (t < 0.0 || t > in.t[in.n-1]) {
00160 t1 = t < 0.0 ? -t : 2.*in.t[in.n-1] - t;
00161 i = interface_arc(in, t1);
00162 xi[n] = splint(in.spx[i], t1);
00163 yi[n] = 4. - splint(in.spy[i], t1);
00164 xni[n] = - splint1(in.spy[i], t1);
00165 yni[n] = - splint1(in.spx[i], t1);
00166 dl = sqrt(sq(xni[n]) + sq(yni[n]));
00167 xni[n] /= dl; yni[n] /= dl;
00168 }
00169 else {
00170 i = interface_arc(in, t);
00171 xi[n] = splint(in.spx[i], t);
00172 yi[n] = splint(in.spy[i], t);
00173 xni[n] = - splint1(in.spy[i], t);
00174 yni[n] = splint1(in.spx[i], t);
00175 dl = sqrt(sq(xni[n]) + sq(yni[n]));
00176 xni[n] /= dl; yni[n] /= dl;
00177 }
00178
00179 fptr = fopen("ipoints", "at");
00180 for (i = 0; i < NI; i++)
00181 fprintf(fptr, "%g %g\n%g %g\n%g %g\n\n", x, y, xi[i], yi[i],
00182 xi[i] + xni[i], yi[i] + yni[i]);
00183 fclose(fptr);
00184
00185 if (closest_u(xmin, ymin, u, cu, xu, yu, uu, nx, ny) < NU) {
00186 fprintf(stderr, "extra_poly(%g,%g): not enough u's\n", x, y);
00187 exit(1);
00188 }
00189 if (closest_v(xmin, ymin, v, cv, xv, yv, vv, nx, ny) < NV) {
00190 fprintf(stderr, "extra_poly(%g,%g): not enough v's\n", x, y);
00191 exit(1);
00192 }
00193
00194 fptr = fopen("upoints", "at");
00195 for (i = 0; i < NU; i++)
00196 fprintf(fptr, "%g %g\n%g %g\n\n", x, y, xu[i], yu[i]);
00197 fclose(fptr);
00198 fptr = fopen("vpoints", "at");
00199 for (i = 0; i < NV; i++)
00200 fprintf(fptr, "%g %g\n%g %g\n\n", x, y, xv[i], yv[i]);
00201 fclose(fptr);
00202
00203 if (four_aligned(xu, yu, 4)) {
00204 xu[3] = xu[4]; yu[3] = yu[4]; uu[3] = uu[4];
00205 xu[4] = xu[5]; yu[4] = yu[5]; uu[4] = uu[5];
00206 }
00207 if (four_aligned(xu, yu, 5)) {
00208 xu[4] = xu[5]; yu[4] = yu[5]; uu[4] = uu[5];
00209 }
00210 if (four_aligned(xv, yv, 4)) {
00211 xv[3] = xv[4]; yv[3] = yv[4]; vv[3] = vv[4];
00212 xv[4] = xv[5]; yv[4] = yv[5]; vv[4] = vv[5];
00213 }
00214 if (four_aligned(xv, yv, 5)) {
00215 xv[4] = xv[5]; yv[4] = yv[5]; vv[4] = vv[5];
00216 }
00217
00218 /* build the system of equations */
00219 for (i = 0; i < NU; i++) {
00220 a[i][0] = sq(xu[i]); a[i][1] = sq(yu[i]); a[i][2] = xu[i]*yu[i];
00221 a[i][3] = xu[i]; a[i][4] = yu[i]; a[i][5] = 1.0;
00222 a[i][6] = a[i][7] = a[i][8] = a[i][9] = a[i][10] = a[i][11] = 0.0;
00223 c[i] = uu[i];
00224 }
00225 for (i = 0; i < NV; i++) {
00226 a[i+NU][0] = a[i+NU][1] = a[i+NU][2] =
00227 a[i+NU][3] = a[i+NU][4] = a[i+NU][5] = 0.0;
00228 a[i+NU][6] = sq(xv[i]); a[i+NU][7] = sq(yv[i]); a[i+NU][8] = xv[i]*yv[i];
00229 a[i+NU][9] = xv[i]; a[i+NU][10] = yv[i]; a[i+NU][11] = 1.0;
00230 c[i+NU] = vv[i];
00231 }
00232 for (i = 0; i < NI; i++) {
00233 tt = -xni[i]*yni[i]; tp = 0.5*(sq(yni[i]) - sq(xni[i]));
00234 a[i+NU+NV][0] = -2.*tt*xi[i]; a[i+NU+NV][1] = 2.*tp*yi[i];
00235 a[i+NU+NV][2] = tp*xi[i] - tt*yi[i];
00236 a[i+NU+NV][3] = -tt; a[i+NU+NV][4] = tp; a[i+NU+NV][5] = 0.0;
00237 a[i+NU+NV][6] = 2.*tp*xi[i]; a[i+NU+NV][7] = 2.*tt*yi[i];
00238 a[i+NU+NV][8] = tt*xi[i] + tp*yi[i];
00239 a[i+NU+NV][9] = tp; a[i+NU+NV][10] = tt; a[i+NU+NV][11] = 0.0;
00240 c[i+NU+NV] = 0.0;
00241 }
00242
00243 /*
00244 for (i = 0; i < 12; i++) {
00245 for (j = 0; j < 12; j++)
00246 printf("%g\t", a[i][j]);
00247 printf("\n");
00248 }
00249 printf("\n");
00250 */
00251
00252 if (gauss3(a, c, b, &det, 12)) {
00253 FILE *fptr;
00254 char s[256];
00255 sprintf(s, "nosolution.%d.xmgr", graph++);
00256 fptr = fopen(s, "wt");
00257 fprintf(stderr, "extra_poly(%g,%g): no solution det = %g\n", x, y, det);
00258 fprintf(fptr, "%g %g\n&\n", x, y);
00259 for (i = 0; i < NU; i++)
00260 fprintf(fptr, "%g %g\n", xu[i], yu[i]);
00261 fprintf(fptr, "&\n");
00262 for (i = 0; i < NV; i++)
00263 fprintf(fptr, "%g %g\n", xv[i], yv[i]);
00264 fprintf(fptr, "&\n");
00265 for (i = 0; i < NI; i++)
00266 fprintf(fptr, "%g %g\n", xi[i], yi[i]);
00267 fprintf(fptr, "&\n");
00268 fclose(fptr);
00269 /* exit(1); */
00270 }
00271 /*
00272 for (i = 0; i < 12; i++)
00273 printf("%g\t", b[i]);
00274 printf("\n");
00275 */
00276 *eu = b[0]*sq(x) + b[1]*sq(y) + b[2]*x*y + b[3]*x + b[4]*y + b[5];
00277 *ev = b[6]*sq(x) + b[7]*sq(y) + b[8]*x*y + b[9]*x + b[10]*y + b[11];
00278
00279 {
00280 real maxu = 0, minu = nx, x1, y1, maxv = 0, minv = ny;
00281 char s[256];
00282 static int truc;
00283 sprintf(s, "toto.%d", truc++);
00284 fptr = fopen(s, "wt");
00285
00286 fprintf(fptr, "# %g %g\n# &\n", x, y);
00287 for (i = 0; i < NU; i++)
00288 fprintf(fptr, "# %g %g\n", xu[i], yu[i]);
00289 fprintf(fptr, "# &\n");
00290 for (i = 0; i < NV; i++)
00291 fprintf(fptr, "# %g %g\n", xv[i], yv[i]);
00292 fprintf(fptr, "# &\n");
00293 for (i = 0; i < NI; i++)
00294 fprintf(fptr, "# %g %g\n", xi[i], yi[i]);
00295 fprintf(fptr, "# &\n");
00296
00297 for (i = 0; i < NU; i++) {
00298 maxu = xu[i] > maxu ? xu[i] : maxu;
00299 minu = xu[i] < minu ? xu[i] : minu;
00300 maxv = yu[i] > maxv ? yu[i] : maxv;
00301 minv = yu[i] < minv ? yu[i] : minv;
00302 }
00303 for (x1 = minu; x1 <= maxu; x1 += (maxu - minu)/10.) {
00304 for (y1 = minv; y1 <= maxv; y1 += (maxv - minv)/10.)
00305 fprintf(fptr, "%g %g %g\n", x1, y1,
00306 b[0]*sq(x1) + b[1]*sq(y1) + b[2]*x1*y1 + b[3]*x1 + b[4]*y1 + b[5]);
00307 fprintf(fptr, "\n");
00308 }
00309 fclose(fptr);
00310 }
00311 }
|
1.2.18