---

extra1.c

Go to the documentation of this file.

#include "utilf.h"
#include "markers.h"
#include "extra1.h"
#include "extra.h"
#include "make_bc.h"
#include "inout.h"

#define NU 5
#define NV 5
#define NI 2


int aligned(real x0, real y0, real x1, real y1, 
            real x2, real y2, real x3, real y3)
{
  real dx, dy;
  dx = x1 - x0; dy = y1 - y0;
  if (dx*(y2 - y0) - dy*(x2 - x0) != 0.0 ||
      dx*(y3 - y0) - dy*(x3 - x0) != 0.0)
    return 0;
  return 1;
}


int four_aligned(real *x, real *y, int n)
{
  if (n < 4) return 0;
  if (n == 4) return aligned(x[0], y[0], x[1], y[1], x[2], y[2], x[3], y[3]);
  if (aligned(x[0], y[0], x[1], y[1], x[2], y[2], x[3], y[3]) ||
      aligned(x[1], y[1], x[2], y[2], x[3], y[3], x[4], y[4]) ||
      aligned(x[2], y[2], x[3], y[3], x[4], y[4], x[0], y[0]) ||
      aligned(x[3], y[3], x[4], y[4], x[0], y[0], x[1], y[1]) ||
      aligned(x[4], y[4], x[0], y[0], x[1], y[1], x[2], y[2]))
    return 1;
  return 0;
}


int closest_u(real x, real y, real2D u, real2D cu,
              real *xu, real *yu, real *uu, int nx, int ny)
{
  int i, j, k, i1, j1, n = 0;
  real xl, yl, ul, dl;

  i = x; j = (int) ((real)y - 0.5);
  for (i1 = i - 2; i1 <= i + 2; i1++)
    for (j1 = j - 2; j1 <= j + 2; j1++)
      if (INU(i1,j1)) {
        xu[n] = i1; yu[n] = j1 + 0.5; uu[n++] = u[i1][j1];
      }

  /* sort the points */
  for (k = 1; k < n; k++) {
    xl = xu[k]; yl = yu[k]; ul = uu[k];
    i = k - 1;
    dl = sq(x - xl) + sq(y - yl);
    while (i >= 0 && sq(xu[i] - x) + sq(yu[i] - y) > dl) {
      xu[i+1] = xu[i]; yu[i+1] = yu[i]; uu[i+1] = uu[i]; 
      i--;
    }
    xu[i+1] = xl; yu[i+1] = yl; uu[i+1] = ul;
  }
  return n;
}


int closest_v(real x, real y, real2D v, real2D cv,
              real *xv, real *yv, real *vv, int nx, int ny)
{
  int i, j, k, i1, j1, n = 0;
  real xl, yl, vl, dl;

  i = (int)((real)x - 0.5); j = y;
  for (i1 = i - 2; i1 <= i + 2; i1++)
    for (j1 = j - 2; j1 <= j + 2; j1++)
      if (INV(i1,j1)) {
        xv[n] = i1 + 0.5; yv[n] = j1; vv[n++] = v[i1][j1];
      }

  /* sort the points */
  for (k = 1; k < n; k++) {
    xl = xv[k]; yl = yv[k]; vl = vv[k]; 
    i = k - 1;
    dl = sq(x - xl) + sq(y - yl);
    while (i >= 0 && sq(xv[i] - x) + sq(yv[i] - y) > dl) {
      xv[i+1] = xv[i]; yv[i+1] = yv[i]; vv[i+1] = vv[i]; 
      i--;
    }
    xv[i+1] = xl; yv[i+1] = yl; vv[i+1] = vl;
  }
  return n;
}


int gauss3(real A[12][12], real B[12], real C[12], real *det, int n)
{
  int i, j, k, i1;
  real max, temp;
  
  for (j = 0; j <= n - 2; j++) {
    i1 = j;
    max = fabs(A[j][j]);
    for (i = j + 1; i < n; i++) {
      if (fabs(A[i][j]) >= max) {
        i1 = i;
        max = fabs(A[i][j]);
      }
    }
    for (k = j; k < n; k++) {
      max = A[j][k]; A[j][k] = A[i1][k]; A[i1][k] = max;
    }
    max = B[j]; B[j] = B[i1]; B[i1] = max;
    for (i = j + 1; i < n; i++) {
      if (A[j][j] == 0.)
        return 1;
      else
        temp = A[i][j] / A[j][j];
      for (k = j; k < n; k++)
        A[i][k] -= temp * A[j][k];
      B[i] -= temp * B[j];
    }
  }

  *det = 1.0;
  for (i = 0; i < n; i++)
    *det *= A[i][i];
  if (fabs(*det) < 1.e-10) return 1;

  C[n - 1] = B[n - 1] / A[n - 1][n - 1];
  for (i = n - 2; i >= 0; i--) {
    C[i] = B[i];
    for (k = i + 1; k < n; k++)
      C[i] -= A[i][k] * C[k];
    C[i] /= A[i][i];
  }
  return 0;
}


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)
{
  real xu[25], yu[25], uu[25];
  real xv[25], yv[25], vv[25];
  real xi[NI], yi[NI], xni[NI], yni[NI];
  real xmin, ymin, pmin, tmin, xn, yn, t, t1, dl;
  int i, j, n;
  real a[12][12], b[12], c[12], tt, tp, det;
FILE *fptr;
static int graph, graph1;

  if (!closest_normal(x, y, in, &xmin, &ymin, &pmin, &tmin, &xn, &yn)) {
    fprintf(stderr, "extra_poly(%g,%g): no closest normal\n", x, y);
    exit(1);
  }
  for (n = 0, t = tmin - 0.5; n < NI; n++, t += 1.0)
    if (t < 0.0 || t > in.t[in.n-1]) {
      t1 = t < 0.0 ? -t : 2.*in.t[in.n-1] - t;
      i = interface_arc(in, t1);
      xi[n] = splint(in.spx[i], t1);
      yi[n] = 4. - splint(in.spy[i], t1);
      xni[n] = - splint1(in.spy[i], t1);
      yni[n] = - splint1(in.spx[i], t1);
      dl = sqrt(sq(xni[n]) + sq(yni[n]));
      xni[n] /= dl; yni[n] /= dl;
    }
    else {
      i = interface_arc(in, t);
      xi[n] = splint(in.spx[i], t);
      yi[n] = splint(in.spy[i], t);
      xni[n] = - splint1(in.spy[i], t);
      yni[n] = splint1(in.spx[i], t);
      dl = sqrt(sq(xni[n]) + sq(yni[n]));
      xni[n] /= dl; yni[n] /= dl;
    }

  fptr = fopen("ipoints", "at");
  for (i = 0; i < NI; i++)
    fprintf(fptr, "%g %g\n%g %g\n%g %g\n\n", x, y, xi[i], yi[i],
            xi[i] + xni[i], yi[i] + yni[i]);
  fclose(fptr);
 
  if (closest_u(xmin, ymin, u, cu, xu, yu, uu, nx, ny) < NU) {
    fprintf(stderr, "extra_poly(%g,%g): not enough u's\n", x, y);
    exit(1);
  }
  if (closest_v(xmin, ymin, v, cv, xv, yv, vv, nx, ny) < NV) {
    fprintf(stderr, "extra_poly(%g,%g): not enough v's\n", x, y);
    exit(1);
  }
  
  fptr = fopen("upoints", "at");
  for (i = 0; i < NU; i++)
    fprintf(fptr, "%g %g\n%g %g\n\n", x, y, xu[i], yu[i]);
  fclose(fptr);
  fptr = fopen("vpoints", "at");
  for (i = 0; i < NV; i++)
    fprintf(fptr, "%g %g\n%g %g\n\n", x, y, xv[i], yv[i]);
  fclose(fptr);

  if (four_aligned(xu, yu, 4)) {
    xu[3] = xu[4]; yu[3] = yu[4]; uu[3] = uu[4];
    xu[4] = xu[5]; yu[4] = yu[5]; uu[4] = uu[5];
  }
  if (four_aligned(xu, yu, 5)) {
    xu[4] = xu[5]; yu[4] = yu[5]; uu[4] = uu[5];
  }
  if (four_aligned(xv, yv, 4)) {
    xv[3] = xv[4]; yv[3] = yv[4]; vv[3] = vv[4];
    xv[4] = xv[5]; yv[4] = yv[5]; vv[4] = vv[5];
  }
  if (four_aligned(xv, yv, 5)) {
    xv[4] = xv[5]; yv[4] = yv[5]; vv[4] = vv[5];
  }
  
  /* build the system of equations */
  for (i = 0; i < NU; i++) {
    a[i][0] = sq(xu[i]); a[i][1] = sq(yu[i]); a[i][2] = xu[i]*yu[i];
    a[i][3] = xu[i]; a[i][4] = yu[i]; a[i][5] = 1.0;
    a[i][6] = a[i][7] = a[i][8] = a[i][9] = a[i][10] = a[i][11] = 0.0;
    c[i] = uu[i];
  }
  for (i = 0; i < NV; i++) {
    a[i+NU][0] = a[i+NU][1] = a[i+NU][2] = 
      a[i+NU][3] = a[i+NU][4] = a[i+NU][5] = 0.0;
    a[i+NU][6] = sq(xv[i]); a[i+NU][7] = sq(yv[i]); a[i+NU][8] = xv[i]*yv[i];
    a[i+NU][9] = xv[i]; a[i+NU][10] = yv[i]; a[i+NU][11] = 1.0;
    c[i+NU] = vv[i];
  }
  for (i = 0; i < NI; i++) {
    tt = -xni[i]*yni[i]; tp = 0.5*(sq(yni[i]) - sq(xni[i]));
    a[i+NU+NV][0] = -2.*tt*xi[i]; a[i+NU+NV][1] = 2.*tp*yi[i];
    a[i+NU+NV][2] = tp*xi[i] - tt*yi[i];
    a[i+NU+NV][3] = -tt; a[i+NU+NV][4] = tp; a[i+NU+NV][5] = 0.0;
    a[i+NU+NV][6] = 2.*tp*xi[i]; a[i+NU+NV][7] = 2.*tt*yi[i];
    a[i+NU+NV][8] = tt*xi[i] + tp*yi[i];
    a[i+NU+NV][9] = tp; a[i+NU+NV][10] = tt; a[i+NU+NV][11] = 0.0;
    c[i+NU+NV] = 0.0;
  }

  /*
    for (i = 0; i < 12; i++) {
    for (j = 0; j < 12; j++)
    printf("%g\t", a[i][j]);
    printf("\n");
    }
    printf("\n");
  */

  if (gauss3(a, c, b, &det, 12)) {
    FILE *fptr;
    char s[256];
    sprintf(s, "nosolution.%d.xmgr", graph++);
    fptr = fopen(s, "wt");
    fprintf(stderr, "extra_poly(%g,%g): no solution det = %g\n", x, y, det);
    fprintf(fptr, "%g %g\n&\n", x, y);
    for (i = 0; i < NU; i++)
      fprintf(fptr, "%g %g\n", xu[i], yu[i]);
    fprintf(fptr, "&\n");
    for (i = 0; i < NV; i++)
      fprintf(fptr, "%g %g\n", xv[i], yv[i]);
    fprintf(fptr, "&\n");
    for (i = 0; i < NI; i++)
      fprintf(fptr, "%g %g\n", xi[i], yi[i]);
    fprintf(fptr, "&\n");
    fclose(fptr);
    /*    exit(1); */
  }
  /*
  for (i = 0; i < 12; i++)
    printf("%g\t", b[i]);
  printf("\n");
  */
  *eu = b[0]*sq(x) + b[1]*sq(y) + b[2]*x*y + b[3]*x + b[4]*y + b[5];
  *ev = b[6]*sq(x) + b[7]*sq(y) + b[8]*x*y + b[9]*x + b[10]*y + b[11];

  {
    real maxu = 0, minu = nx, x1, y1, maxv = 0, minv = ny;
    char s[256];
    static int truc;
    sprintf(s, "toto.%d", truc++);
    fptr = fopen(s, "wt");

    fprintf(fptr, "# %g %g\n# &\n", x, y);
    for (i = 0; i < NU; i++)
      fprintf(fptr, "# %g %g\n", xu[i], yu[i]);
    fprintf(fptr, "# &\n");
    for (i = 0; i < NV; i++)
      fprintf(fptr, "# %g %g\n", xv[i], yv[i]);
    fprintf(fptr, "# &\n");
    for (i = 0; i < NI; i++)
      fprintf(fptr, "# %g %g\n", xi[i], yi[i]);
    fprintf(fptr, "# &\n");

    for (i = 0; i < NU; i++) {
      maxu = xu[i] > maxu ? xu[i] : maxu;
      minu = xu[i] < minu ? xu[i] : minu;
      maxv = yu[i] > maxv ? yu[i] : maxv;
      minv = yu[i] < minv ? yu[i] : minv;
    }
    for (x1 = minu; x1 <= maxu; x1 += (maxu - minu)/10.) {
      for (y1 = minv; y1 <= maxv; y1 += (maxv - minv)/10.)
        fprintf(fptr, "%g %g %g\n", x1, y1,
                b[0]*sq(x1) + b[1]*sq(y1) + b[2]*x1*y1 + b[3]*x1 + b[4]*y1 + b[5]);
      fprintf(fptr, "\n");
    }
    fclose(fptr);
  }
}

---