#include "utilf.h"#include "markers.h"#include "extra1.h"#include "extra.h"#include "make_bc.h"#include "inout.h"Include dependency graph for extra1.c:
Go to the source code of this file.
Defines | |
| #define | NU 5 |
| #define | NV 5 |
| #define | NI 2 |
Functions | |
| int | aligned (real x0, real y0, real x1, real y1, real x2, real y2, real x3, real y3) |
| int | four_aligned (real *x, real *y, int n) |
| 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) |
| int | gauss3 (real A[12][12], real B[12], real C[12], real *det, int n) |
| 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 10 of file extra1.c. Referenced by extra_poly(). |
|
|
Definition at line 8 of file extra1.c. Referenced by extra_poly(). |
|
|
Definition at line 9 of file extra1.c. Referenced by extra_poly(). |
|
||||||||||||||||||||||||||||||||||||
|
Definition at line 13 of file extra1.c. References real. Referenced by four_aligned().
00015 {
00016 real dx, dy;
00017 dx = x1 - x0; dy = y1 - y0;
00018 if (dx*(y2 - y0) - dy*(x2 - x0) != 0.0 ||
00019 dx*(y3 - y0) - dy*(x3 - x0) != 0.0)
00020 return 0;
00021 return 1;
00022 }
|
|
||||||||||||||||||||||||||||||||||||||||
|
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 }
|
|
||||||||||||||||
|
Definition at line 25 of file extra1.c. References aligned(), and real. Referenced by extra_poly().
00026 {
00027 if (n < 4) return 0;
00028 if (n == 4) return aligned(x[0], y[0], x[1], y[1], x[2], y[2], x[3], y[3]);
00029 if (aligned(x[0], y[0], x[1], y[1], x[2], y[2], x[3], y[3]) ||
00030 aligned(x[1], y[1], x[2], y[2], x[3], y[3], x[4], y[4]) ||
00031 aligned(x[2], y[2], x[3], y[3], x[4], y[4], x[0], y[0]) ||
00032 aligned(x[3], y[3], x[4], y[4], x[0], y[0], x[1], y[1]) ||
00033 aligned(x[4], y[4], x[0], y[0], x[1], y[1], x[2], y[2]))
00034 return 1;
00035 return 0;
00036 }
|
|
||||||||||||||||||||||||
|
Definition at line 95 of file extra1.c. References real. Referenced by extra_poly().
00096 {
00097 int i, j, k, i1;
00098 real max, temp;
00099
00100 for (j = 0; j <= n - 2; j++) {
00101 i1 = j;
00102 max = fabs(A[j][j]);
00103 for (i = j + 1; i < n; i++) {
00104 if (fabs(A[i][j]) >= max) {
00105 i1 = i;
00106 max = fabs(A[i][j]);
00107 }
00108 }
00109 for (k = j; k < n; k++) {
00110 max = A[j][k]; A[j][k] = A[i1][k]; A[i1][k] = max;
00111 }
00112 max = B[j]; B[j] = B[i1]; B[i1] = max;
00113 for (i = j + 1; i < n; i++) {
00114 if (A[j][j] == 0.)
00115 return 1;
00116 else
00117 temp = A[i][j] / A[j][j];
00118 for (k = j; k < n; k++)
00119 A[i][k] -= temp * A[j][k];
00120 B[i] -= temp * B[j];
00121 }
00122 }
00123
00124 *det = 1.0;
00125 for (i = 0; i < n; i++)
00126 *det *= A[i][i];
00127 if (fabs(*det) < 1.e-10) return 1;
00128
00129 C[n - 1] = B[n - 1] / A[n - 1][n - 1];
00130 for (i = n - 2; i >= 0; i--) {
00131 C[i] = B[i];
00132 for (k = i + 1; k < n; k++)
00133 C[i] -= A[i][k] * C[k];
00134 C[i] /= A[i][i];
00135 }
00136 return 0;
00137 }
|
1.2.18