---

extra.c

Go to the documentation of this file.

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

#define DIST 2.5
#define SQR2 1.4142136
#define SYMMETRY(val) (val < 2.0 ? 4. - val : val)
#define BC_YU(j,ny) (j <= 1 ? 3 - (j) : j >= ny ? 2*(ny) - (j) - 1 : j)
#define BC_YV(j,ny) (j < 2 ? 4 - (j) : j > ny ? 2*(ny) - (j) : j)

#define THETA1 0.463647609001   /* atan(0.5) */
#define THETA2 0.785398163398   /* PI/4 */
#define THETA3 1.10714871779    /* atan(2) */
#define THETA4 1.5707963268     /* PI/2 */

#define PREC 1e-3

static int gauss(real A[7][7], real B[7], real C[7], int n);


/* normal distance between a point and a segment. if the foot of the normal
   is not in the segment, the minimum distance between each of the end points
   is returned.
*/
static real normal_dist(real xa, real ya, real xb, real yb, real x, real y,
                        int *which_end)
{
  real ab, abac, ac, bc;
  ab = sqrt(sq(xb - xa) + sq(yb - ya));
  if (ab == 0.0) {
    *which_end = 1;
    return sqrt(sq(x - xa) + sq(y - ya));
  }
  abac = (xb - xa)*(x - xa) + (yb - ya)*(y - ya);
  if (abac < 0.0 || abac > sq(ab)) {
    ac = sqrt(sq(x - xa) + sq(y - ya));
    bc = sqrt(sq(x - xb) + sq(y - yb));
    if (ac < bc) { *which_end = 1; return ac; }
    *which_end = 2;
    return bc;
  }
  *which_end = 0;
  ac = sqrt(sq(x - xa) + sq(y - ya));
  bc = sqrt(sq(x - xb) + sq(y - yb));
  return fabs((xb - xa)*(y - ya) - (yb - ya)*(x - xa))/ab;
}


int closest_normal(real x, real y, interface in,
                   real *xmin, real *ymin, real *pmin, real *tmin,
                   real *xn, real *yn)
{
  int i, span;
  real d, dmin, dl;
  int imin, found = 0, which;

  /* find the closest segment */
#if 0
  dmin = normal_dist(in.x[0], in.y[0], in.x[1], in.y[1], x, y, &which);
  imin = 0;
  for (i = 1; i < in.n - 1; i++) {
    d = normal_dist(in.x[i], in.y[i], in.x[i+1], in.y[i+1], x, y, &which);
    if (d < dmin) { 
      dmin = d;
      if (which == 0 || which == 1) imin = i;
      else imin = i + 1;
    }
  }
#else
  /* find the closest marker */
  dmin = sq(x - in.x[0]) + sq(y - in.y[0]); imin = 0;
  for (i = 1; i < in.n; i++) {
    d = sq(x - in.x[i]) + sq(y - in.y[i]);
    if (d < dmin) { dmin = d; imin = i; }
  }  
#endif

  /* find the closest point on the interface */
  i = imin; span = 1;
  if (i < in.n - 1)
    found = distmin(in.spx[i], in.spy[i], x, y, tmin,
                    in.t[i], in.t[i+1], PREC);
  while (!found && (i + span < in.n - 1 || i - span >= 0)) {
    if (i + span < in.n - 1 &&
        (found = distmin(in.spx[i+span], in.spy[i+span], x, y, tmin,
                         in.t[i+span], in.t[i+span+1], PREC)))
      imin = i + span;
    else if (i - span >= 0 &&
             (found = distmin(in.spx[i-span], in.spy[i-span], x, y, tmin,
                              in.t[i-span], in.t[i-span+1], PREC)))
      imin = i - span;    
    span++;
  }
  if (!found) { 
    FILE *fptr = fopen("closest_splines", "wt");
    interface_write_splines(in, fptr, 1000);
    fclose(fptr);
    printf("&&&& %g %g %d %g\n", x, y, imin, dmin); 
      printf("%g %g\n%g %g\n%g %g\n", in.x[imin-1], in.y[imin-1],
             in.x[imin], in.y[imin],
             in.x[imin+1], in.y[imin+1]);
      return 0; 
  }

  *xn = - splint1(in.spy[imin], *tmin);
  *yn = splint1(in.spx[imin], *tmin);

  dl = sqrt(sq(*xn) + sq(*yn));
  *xn /= dl; *yn /= dl;
  *xmin = splint(in.spx[imin], *tmin);
  *ymin = splint(in.spy[imin], *tmin);
  *pmin = (*tmin - in.t[imin])/(in.t[imin+1] - in.t[imin])
    *(in.p[imin+1] - in.p[imin]) + in.p[imin];

  return 1;
}


#define NCLOSEST 3
#define SQCLOSEST sq(2*NCLOSEST + 1)
 #define CLOSENESS(xu, yu, x, y, n1, n2)  (1.1*(sq(((xu) - (x))*(n1)) + \
    sq(((yu) - (y))*(n2))) + \
    2.1*sq((n2)*((xu) - (x)) - n1*((yu) - (y))))
/* #define CLOSENESS(xu, yu, x, y, n1, n2)  MAX(sq(((xu) - (x))*(n1)) + \
   sq(((yu) - (y))*(n2)), \
   2.1*sq((n2)*((xu) - (x)) - n1*((yu) - (y)))) */

int closest_u(real x, real y, real n1, real n2, real2D u, real2D ap,
              real *xu, real *yu, real *uu, 
              real *xmax, real *xmin, real *ymax, real *ymin,
              int nx, int ny)
{
  int i, j, k, i1, j1, j2, n = 0;
  real xl, yl, ul, dl;

  *xmax = 0.0; *xmin = 2.*nx; *ymax = 0.0; *ymin = 2.*ny;

  i = x; j = (int) ((real)y - 0.5);
  for (i1 = i - NCLOSEST; i1 <= i + NCLOSEST; i1++)
    for (j1 = j - NCLOSEST; j1 <= j + NCLOSEST; j1++) {
      j2 = BC_YU(j1,ny);
      if (INU(i1,j2)) {
        xu[n] = i1; yu[n] = j1 + 0.5;
        *xmax = xu[n] > *xmax ? xu[n] : *xmax;
        *ymax = yu[n] > *ymax ? yu[n] : *ymax;
        *xmin = xu[n] < *xmin ? xu[n] : *xmin;
        *ymin = yu[n] < *ymin ? yu[n] : *ymin;
        uu[n++] = u[i1][j2];
      }
    }

  /* sort the points */
  for (k = 1; k < n; k++) {
    xl = xu[k]; yl = yu[k]; ul = uu[k];
    i = k - 1;
    dl = CLOSENESS(xl, yl, x, y, n1, n2);
    while (i >= 0 && CLOSENESS(xu[i], yu[i], x, y, n1, n2) > 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, real n1, real n2, real2D v, real2D ap,
              real *xv, real *yv, real *vv, 
              real *xmax, real *xmin, real *ymax, real *ymin,
              int nx, int ny)
{
  int i, j, k, i1, j1, j2, n = 0;
  real xl, yl, vl, dl;

  *xmax = 0.0; *xmin = 2.*nx; *ymax = 0.0; *ymin = 2.*ny;

  i = (int)((real)x - 0.5); j = y;
  for (i1 = i - NCLOSEST; i1 <= i + NCLOSEST; i1++)
    for (j1 = j - NCLOSEST; j1 <= j + NCLOSEST; j1++) {
      j2 = BC_YV(j1,ny);
      if (INV(i1,j2)) {
        xv[n] = i1 + 0.5; yv[n] = j1;
        *xmax = xv[n] > *xmax ? xv[n] : *xmax;
        *ymax = yv[n] > *ymax ? yv[n] : *ymax;
        *xmin = xv[n] < *xmin ? xv[n] : *xmin;
        *ymin = yv[n] < *ymin ? yv[n] : *ymin;
        if (j1 == j2)
          vv[n++] = v[i1][j1];
        else
          vv[n++] = - v[i1][j2];
      }
    }

  /* sort the points */
  for (k = 1; k < n; k++) {
    xl = xv[k]; yl = yv[k]; vl = vv[k]; 
    i = k - 1;
    dl = CLOSENESS(xl, yl, x, y, n1, n2);
    while (i >= 0 && CLOSENESS(xv[i], yv[i], x, y, n1, n2) > 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;
}

/* #define DEBUG_EXTRA */

void extra_velocity_normals(real x, real y,
                            real n1, real n2, real xn, real yn,
                            real2D u, real2D v, real2D ap,
                            real *eu, real *ev,
                            real *dudx, real *dvdy, real *dudy, real *dvdx,
                            int nx, int ny)
{
  int i, j, k, l, l1, kmin, kmax, lmin, lmax;
  int nsx = 0, nsy = 0, dx, dy;
  real xSx = 0., xSxx = 0., xSxy = 0., xSy = 0., xSyy = 0.;
  real ySx = 0., ySxx = 0., ySxy = 0., ySy = 0., ySyy = 0.;
  real xSu = 0., xSxu = 0., xSyu = 0.;
  real ySv = 0., ySxv = 0., ySyv = 0.;
  real sys[7][7] = {{0.,0.,0.,0.,0.,0.,0.},
                    {0.,0.,0.,0.,0.,0.,0.},
                    {0.,0.,0.,0.,0.,0.,0.},
                    {0.,0.,0.,0.,0.,0.,0.},
                    {0.,0.,0.,0.,0.,0.,0.},
                    {0.,0.,0.,0.,0.,0.,0.},
                    {0.,0.,0.,0.,0.,0.,0.}};
  real rhs[7] = {0.,0.,0.,0.,0.,0.,0.}, sol[7], dummy;
#ifdef DEBUG_EXTRA  
FILE *fptr1 = fopen("extra_upoints", "at");
FILE *fptr2 = fopen("extra_vpoints", "at");
FILE *fptr3 = fopen("extra_points", "at");
FILE *fptr4 = fopen("extra_normals", "at");
#endif

#if 0
  /* find the velocity points using a 5x5 or 4x4 stencil */
  i = x + 0.5; j = y;
  dx = dy = 1;

  do {
    nsx = 0;
    xSx = xSxx = xSxy = xSy = xSyy = 0.;
    xSu = xSxu = xSyu = 0.;
    
    if (modf(x + 0.5, &dummy) == 0.0) { kmin = i - dx; kmax = i + dx - 1; }
    else { kmin = i - dx; kmax = i + dx; }

    if (modf(y, &dummy) == 0.0) { lmin = j - dy; lmax = j + dy - 1; }
    else { lmin = j - dy; lmax = j + dy; }

    for (k = kmin; k <= kmax; k++)
      for (l = lmin; l <= lmax; l++) {
        l1 = l < 2 ? 3 - l : l;
        if (INU(k,l1)) {
          xSx += (real)k; xSxx += sq((real)k);
          xSxy += (0.5 + (real)l)*(real)k;
          xSy += 0.5 + (real)l; xSyy += sq(0.5 + (real)l);
          xSu += u[k][l1]; xSxu += u[k][l1]*(real)k;
          xSyu += u[k][l1]*(0.5 + (real)l);
#ifdef DEBUG_EXTRA
          fprintf(fptr1, "%g %g %g\n", (real)k, 0.5 + (real)l, u[k][l1]);
#endif
          nsx++;
        }
      }
    dx++; dy++;
  } while((nsx <=  kmax - kmin + 1 || nsx <= lmax - lmin + 1) && dx <= 3);
  
  i = x; j = y + 0.5;
  dx = dy = 1;
  do {
    nsy = 0;
    ySx = ySxx = ySxy = ySy = ySyy = 0.;
    ySv = ySxv = ySyv = 0.;

    if (modf(x, &dummy) == 0.0) { kmin = i - dx; kmax = i + dx - 1; }
    else { kmin = i - dx; kmax = i + dx; }

    if (modf(y + 0.5, &dummy) == 0.0) { lmin = j - dy; lmax = j + dy - 1; }
    else { lmin = j - dy; lmax = j + dy; }

    for (k = kmin; k <= kmax; k++)
      for (l = lmin; l <= lmax; l++) {
        l1 = l < 2 ? 4 - l : l;
        if (INV(k,l1)) {
          ySx += 0.5 + (real)k; ySxx += sq(0.5 + (real)k); 
          ySxy += (0.5 + (real)k)*(real)l;
          ySy += (real)l; ySyy += sq((real)l);
          dummy = l1 == l ? v[k][l1] : - v[k][l1];
          ySv += dummy; ySxv += dummy*(0.5 + (real)k);
          ySyv += dummy*(real)l;
#ifdef DEBUG_EXTRA
          fprintf(fptr2, "%g %g %g\n", 0.5+(real)k, (real)l, dummy);
#endif
          nsy++;
        }
      }
    dx++; dy++;
  } while((nsy <= kmax - kmin + 1 || nsy <= lmax - lmin + 1) && dx <= 3);
#else
  real xu[SQCLOSEST], yu[SQCLOSEST], uu[SQCLOSEST];
  real xv[SQCLOSEST], yv[SQCLOSEST], vv[SQCLOSEST];
  real xmaxu, xminu, ymaxu, yminu;
  real xmaxv, xminv, ymaxv, yminv;
  int n;

  nsx = nsy = 5;
  if ((n = closest_u(xn, yn, n1, n2, u, ap, xu, yu, uu, 
                     &xmaxu, &xminu, &ymaxu, &yminu, nx, ny)) < nsx) {
    printf("extra_velocity_normals(%g,%g): not enough u's\n", x, y);
    printf("xn: %g yn: %g n1: %g n2: %g\n", xn, yn, n1, n2);
    for (i = 0; i < n; i++)
      printf("  %g %g\n", xu[i], yu[i]);
    exit(1);
  }
  if ((n = closest_v(xn, yn, n1, n2, v, ap, xv, yv, vv, 
                     &xmaxv, &xminv, &ymaxv, &yminv, nx, ny)) < nsy) {
    printf("extra_velocity_normals(%g,%g): not enough v's\n", x, y);
    printf("xn: %g yn: %g n1: %g n2: %g\n", xn, yn, n1, n2);
    for (i = 0; i < n; i++)
      printf("  %g %g\n", xv[i], yv[i]);
    exit(1);
  }

  xSx = xSxx = xSxy = xSy = xSyy = 0.;
  xSu = xSxu = xSyu = 0.;
  for (i = 0; i < nsx; i++) {
    xSx += xu[i]; xSxx += sq(xu[i]);
    xSxy += xu[i]*yu[i];
    xSy += yu[i]; xSyy += sq(yu[i]);
    xSu += uu[i]; xSxu += uu[i]*xu[i];
    xSyu += uu[i]*yu[i];
#ifdef DEBUG_EXTRA
    fprintf(fptr1, "%g %g %g\n", xu[i], yu[i], uu[i]);
#endif
  }

  ySx = ySxx = ySxy = ySy = ySyy = 0.;
  ySv = ySxv = ySyv = 0.;
  for (i = 0; i < nsy; i++) {
    ySx += xv[i]; ySxx += sq(xv[i]); 
    ySxy += xv[i]*yv[i];
    ySy += yv[i]; ySyy += sq(yv[i]);
    ySv += vv[i]; ySxv += vv[i]*xv[i];
    ySyv += yv[i]*vv[i];
#ifdef DEBUG_EXTRA
    fprintf(fptr2, "%g %g %g\n", xv[i], yv[i], vv[i]);
#endif
  }
#endif

  /* build the system of equations
     sys[i][0]*u0+sys[i][1]*v0+sys[i][2]*A11+sys[i][3]*A22+
     sys[i][4]*A12+sys[i][5]*A21+sys[i][6]*lambda=rhs[i] */
  sys[0][2] = n1*n2; sys[0][3] = -n1*n2;
  sys[0][4] = (sq(n2) - sq(n1))/2.0; sys[0][5] = sys[0][4];
  sys[1][0] = nsx; sys[1][2] = xSx; sys[1][4] = xSy;
  rhs[1] = xSu;
  sys[2][1] = nsy; sys[2][3] = ySy; sys[2][5] = ySx;
  rhs[2] = ySv;
  sys[3][0] = xSx; sys[3][2] = xSxx; sys[3][4] = xSxy;
  sys[3][6] = n1*n2/2.;
  rhs[3] = xSxu;
  sys[4][1] = ySy; sys[4][3] = ySyy; sys[4][5] = ySxy;
  sys[4][6] = - n1*n2/2.;
  rhs[4] = ySyv;
  sys[5][0] = xSy; sys[5][2] = xSxy; sys[5][4] = xSyy;
  sys[5][6] = (sq(n2) - sq(n1))/4.0;
  rhs[5] = xSyu;
  sys[6][1] = ySx; sys[6][3] = ySxy; sys[6][5] = ySxx;
  sys[6][6] = (sq(n2) - sq(n1))/4.0;
  rhs[6] = ySxv;

#ifdef DEBUG_EXTRA
  fprintf(fptr1, "&\n\n");
  fprintf(fptr2, "&\n\n");
  fprintf(fptr3, "%g %g\n&\n\n", x, y);
  fprintf(fptr4, "%g %g\n%g %g\n&\n\n", xn, yn, xn + n1, yn + n2);
  fclose(fptr1); fclose(fptr2); fclose(fptr3); fclose(fptr4);
#endif

  if (gauss(sys, rhs, sol, 7)) {
    printf("Error: extra_velocity(%g, %g): no solution\n", x, y);
    printf("       nsx: %d nsy: %d\n", nsx, nsy);
    for (i = 0; i < nsx; i++)
      printf("%g %g %g\n", xu[i], yu[i], uu[i]);
    printf("\n");
    for (i = 0; i < nsy; i++)
      printf("%g %g %g\n", xv[i], yv[i], vv[i]);
    printf("\n");
    for (i = 0; i < 7; i++) {
      for (j = 0; j < 7; j++)
        printf("%10g ", sys[i][j]);
      printf("\n");
    }
    exit(1);
  }

  *dudx = sol[2]; *dvdy = sol[3]; *dudy = sol[4]; *dvdx = sol[5];
  *eu = sol[0] + sol[2]*x + sol[4]*y;
  *ev = sol[1] + sol[5]*x + sol[3]*y;

  if (*eu != *eu || *ev != *ev) {
    printf("Error: extra_velocity(%g, %g): NaN\n", x, y);
    printf("       eu: %g ev: %g\n", *eu, *ev);
    printf("       nsx: %d nsy: %d\n", nsx, nsy);
    for (i = 0; i < nsx; i++)
      printf("%g %g %g\n", xu[i], yu[i], uu[i]);
    printf("\n");
    for (i = 0; i < nsy; i++)
      printf("%g %g %g\n", xv[i], yv[i], vv[i]);
    printf("\n");
    for (i = 0; i < 7; i++) {
      for (j = 0; j < 7; j++)
        printf("%10g ", sys[i][j]);
      printf("\n");
    }
    exit(1);
  }
}


void extra_velocity_normals2(real x, real y,
                             real n1, real n2, real xn, real yn,
                             real2D u, real2D v, 
                             real2D ap,
                             real *eu, real *ev,
                             real *dudx, real *dudy, real *dvdx, real *dvdy,
                             int nx, int ny)
{
  int i, j, k, l, kmin, kmax, lmin, lmax;
  int nsx = 0, nsy = 0, dx, dy;
  real x1, y1;
  real xSx, xSxx, xSxy, xSy, xSyy, xSxxy, xSxyy, xSyyy, xSxxx,
    xSxxxy, xSxxyy, xSxyyy, xSyyyy, xSxxxx, xSu, xSxu, xSyu,
    xSxxu, xSyyu, xSxyu;
  real ySx, ySxx, ySxy, ySy, ySyy, ySxxy, ySxyy, ySyyy, ySxxx,
    ySxxxy, ySxxyy, ySxyyy, ySyyyy, ySxxxx, ySv, ySxv, ySyv,
    ySxxv, ySyyv, ySxyv;
  real sys[13][13] = {{0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},
                      {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},
                      {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},
                      {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},
                      {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},
                      {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},
                      {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},
                      {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},
                      {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},
                      {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},
                      {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},
                      {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.},
                      {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}};
  real rhs[13] = {0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.,0.}, sol[13], dummy;
  real n1n2 = n1*n2, sqn = (sq(n2) - sq(n1))/2.;
  real xu[SQCLOSEST], yu[SQCLOSEST], uu[SQCLOSEST];
  real xv[SQCLOSEST], yv[SQCLOSEST], vv[SQCLOSEST];
  real xmaxu, ymaxu, xminu, yminu, xmaxv, ymaxv, xminv, yminv, det;
  /* #define DEBUG_EXTRA */
#ifdef DEBUG_EXTRA  
  static int cas;
  FILE *fptr1, *fptr2, *fptr3, *fptr4;
  char s[256];

  sprintf(s, "extra_upoints%d", cas);
  fptr1 = fopen(s, "wt");
  sprintf(s, "extra_vpoints%d", cas);
  fptr2 = fopen(s, "wt");
  sprintf(s, "extra_u%d", cas);
  fptr3 = fopen(s, "wt");
  sprintf(s, "extra_v%d", cas++);
  fptr4 = fopen(s, "wt");
#endif

  nsx = nsy = 12;
  if (closest_u(xn, yn, n1, n2, u, ap, xu, yu, uu, 
                &xmaxu, &xminu, &ymaxu, &yminu, nx, ny) < nsx) {
    fprintf(stderr, "extra_velocity_normals2(%g,%g): not enough u's\n", x, y);
    exit(1);
  }
  if (closest_v(xn, yn, n1, n2, v, ap, xv, yv, vv, 
                &xmaxv, &xminv, &ymaxv, &yminv, nx, ny) < nsy) {
    fprintf(stderr, "extra_velocity_normals2(%g,%g): not enough v's\n", x, y);
    exit(1);
  }

  xSx = xSxx = xSxy = xSy = xSyy = xSxxy = xSxyy = xSyyy = xSxxx =
    xSxxxy = xSxxyy = xSxyyy = xSyyyy = xSxxxx = xSu = xSxu = xSyu =
    xSxxu = xSyyu = xSxyu = 0.0;
  for (i = 0; i < nsx; i++) {
    x1 = xu[i]; y1 = yu[i];
    
    xSx += x1; xSxx += sq(x1);
    xSxy += x1*y1;
    xSy += y1; xSyy += sq(y1);
    
    xSxxy += sq(x1)*y1;
    xSxyy += x1*sq(y1);
    xSyyy += sq(y1)*y1;
    xSxxx += sq(x1)*x1;
    
    xSxxxy += sq(x1)*x1*y1;
    xSxxyy += sq(x1)*sq(y1);
    xSxyyy += x1*sq(y1)*y1;
    
    xSyyyy += sq(y1)*sq(y1);
    
    xSxxxx += sq(x1)*sq(x1);
    
    xSu += uu[i]; 
    xSxu += x1*uu[i];
    xSyu += y1*uu[i];
    xSxxu += sq(x1)*uu[i];
    xSyyu += sq(y1)*uu[i];
    xSxyu += x1*y1*uu[i];
  }

  ySx = ySxx = ySxy = ySy = ySyy = ySxxy = ySxyy = ySyyy = ySxxx =
    ySxxxy = ySxxyy = ySxyyy = ySyyyy = ySxxxx = ySv = ySxv = ySyv =
    ySxxv = ySyyv = ySxyv = 0.0;
  for (i = 0; i < nsy; i++) {
    x1 = xv[i]; y1 = yv[i];
    
    ySx += x1; ySxx += sq(x1);
    ySxy += x1*y1;
    ySy += y1; ySyy += sq(y1);
    
    ySxxy += sq(x1)*y1;
    ySxyy += x1*sq(y1);
    ySyyy += sq(y1)*y1;
    ySxxx += sq(x1)*x1;
    
    ySxxxy += sq(x1)*x1*y1;
    ySxxyy += sq(x1)*sq(y1);
    ySxyyy += x1*sq(y1)*y1;
    
    ySyyyy += sq(y1)*sq(y1);
    
    ySxxxx += sq(x1)*sq(x1);
    
    ySv += vv[i]; 
    ySxv += x1*vv[i];
    ySyv += y1*vv[i];
    ySxxv += sq(x1)*vv[i];
    ySyyv += sq(y1)*vv[i];
    ySxyv += x1*y1*vv[i];
  }

  sys[0][0] = 2.*n1n2*xn; sys[0][1] = 2.*sqn*yn; 
  sys[0][2] = n1n2*yn + sqn*xn; sys[0][3] = n1n2; sys[0][4] = sqn;
  sys[0][6] = 2.*sqn*xn; sys[0][7] = - 2.*n1n2*yn;
  sys[0][8] = sqn*yn - n1n2*xn; sys[0][9] = sqn; sys[0][10] = - n1n2;  

  sys[1][0] = xSxx; sys[1][1] = xSyy; sys[1][2] = xSxy;
  sys[1][3] = xSx; sys[1][4] = xSy; sys[1][5] = nsx;  
  rhs[1] = xSu;
  
  sys[2][0] = xSxxy; sys[2][1] = xSyyy; sys[2][2] = xSxyy;
  sys[2][3] = xSxy; sys[2][4] = xSyy; sys[2][5] = xSy;
  sys[2][12] = sqn/2.;
  rhs[2] = xSyu;

  sys[3][0] = xSxxx; sys[3][1] = xSxyy; sys[3][2] = xSxxy;
  sys[3][3] = xSxx; sys[3][4] = xSxy; sys[3][5] = xSx;
  sys[3][12] = n1n2/2.;
  rhs[3] = xSxu;
  
  sys[4][0] = xSxxxy; sys[4][1] = xSxyyy; sys[4][2] = xSxxyy;
  sys[4][3] = xSxxy; sys[4][4] = xSxyy; sys[4][5] = xSxy;
  sys[4][12] = (n1n2*yn + sqn*xn)/2.;
  rhs[4] = xSxyu;

  sys[5][0] = xSxxyy; sys[5][1] = xSyyyy; sys[5][2] = xSxyyy;
  sys[5][3] = xSxyy; sys[5][4] = xSyyy; sys[5][5] = xSyy;
  sys[5][12] = sqn*yn;
  rhs[5] = xSyyu;

  sys[6][0] = xSxxxx; sys[6][1] = xSxxyy; sys[6][2] = xSxxxy;
  sys[6][3] = xSxxx; sys[6][4] = xSxxy; sys[6][5] = xSxx;
  sys[6][12] = n1n2*xn;
  rhs[6] = xSxxu;


  sys[7][6] = ySxx; sys[7][7] = ySyy; sys[7][8] = ySxy;
  sys[7][9] = ySx; sys[7][10] = ySy; sys[7][11] = nsy;  
  rhs[7] = ySv;
  
  sys[8][6] = ySxxy; sys[8][7] = ySyyy; sys[8][8] = ySxyy;
  sys[8][9] = ySxy; sys[8][10] = ySyy; sys[8][11] = ySy;
  sys[8][12] = - n1n2/2.;
  rhs[8] = ySyv;

  sys[9][6] = ySxxx; sys[9][7] = ySxyy; sys[9][8] = ySxxy;
  sys[9][9] = ySxx; sys[9][10] = ySxy; sys[9][11] = ySx;
  sys[9][12] = sqn/2.;
  rhs[9] = ySxv;
  
  sys[10][6] = ySxxxy; sys[10][7] = ySxyyy; sys[10][8] = ySxxyy;
  sys[10][9] = ySxxy; sys[10][10] = ySxyy; sys[10][11] = ySxy;
  sys[10][12] = (- n1n2*xn + sqn*yn)/2.;
  rhs[10] = ySxyv;

  sys[11][6] = ySxxyy; sys[11][7] = ySyyyy; sys[11][8] = ySxyyy;
  sys[11][9] = ySxyy; sys[11][10] = ySyyy; sys[11][11] = ySyy;
  sys[11][12] = - n1n2*yn;
  rhs[11] = ySyyv;

  sys[12][6] = ySxxxx; sys[12][7] = ySxxyy; sys[12][8] = ySxxxy;
  sys[12][9] = ySxxx; sys[12][10] = ySxxy; sys[12][11] = ySxx;
  sys[12][12] = sqn*xn;
  rhs[12] = ySxxv;

#ifdef DEBUG_EXTRA
  for (i = 0; i < nsx; i++)
    fprintf(fptr1, "%g %g %f\n", xu[i], yu[i], uu[i]);
  for (i = 0; i < nsy; i++)
    fprintf(fptr2, "%g %g %f\n", xv[i], yv[i], vv[i]);
  fclose(fptr1); fclose(fptr2);
#endif

  if (gauss13(sys, rhs, sol, &det, 13)) {
    printf("Error: extra_velocity2(%g, %g): no solution\n", x, y);
    printf("       nsx: %d nsy: %d\n", nsx, nsy);
    for (i = 0; i < 13; i++) {
      for (j = 0; j < 13; j++)
        printf("%5g ", sys[i][j]);
      printf("\n");
    }
    exit(1);
  }

#ifdef DEBUG_EXTRA
  for (x1 = xminu; x1 <= xmaxu; x1 += 1.0) {
    for (y1 = yminu; y1 <= ymaxu; y1 += 1.0)
      fprintf(fptr3, "%g %g %g\n", x1, y1, 
              sol[0]*sq(x1) + sol[1]*sq(y1) + sol[2]*x1*y1 + 
              sol[3]*x1 + sol[4]*y1 + sol[5]);
    fprintf(fptr3, "\n");
  }
  fclose(fptr3);
  for (x1 = xminv; x1 <= xmaxv; x1 += 1.0) {
    for (y1 = yminv; y1 <= ymaxv; y1 += 1.0)
      fprintf(fptr4, "%g %g %g\n", x1, y1, 
              sol[6]*sq(x1) + sol[7]*sq(y1) + sol[8]*x1*y1 + 
              sol[9]*x1 + sol[10]*y1 + sol[11]);
    fprintf(fptr4, "\n");
  }
  fclose(fptr4);
  x1 = 0.0;
  for (i = 0; i < nsx; i++) {
    y1 = sol[0]*sq(xu[i]) + sol[1]*sq(yu[i]) + sol[2]*xu[i]*yu[i] + 
      sol[3]*xu[i] + sol[4]*yu[i] + sol[5];
    x1 += sq(y1 - uu[i]);
  }
  for (i = 0; i < nsy; i++) {
    y1 = sol[6]*sq(xv[i]) + sol[7]*sq(yv[i]) + sol[8]*xv[i]*yv[i] + 
      sol[9]*xv[i] + sol[10]*yv[i] + sol[11];
    x1 += sq(y1 - vv[i]);
  }
  printf("%g %g %g\n", x, y, x1);
#endif

  *eu = sol[0]*sq(x) + sol[1]*sq(y) + sol[2]*x*y + 
    sol[3]*x + sol[4]*y + sol[5];
  *ev = sol[6]*sq(x) + sol[7]*sq(y) + sol[8]*x*y + 
    sol[9]*x + sol[10]*y + sol[11];
  *dudx = 2.*sol[0]*x + sol[2]*y + sol[3];
  *dudy = 2.*sol[1]*y + sol[2]*x + sol[4];
  *dvdx = 2.*sol[6]*x + sol[8]*y + sol[9];
  *dvdy = 2.*sol[7]*y + sol[8]*x + sol[10];
}


void extra_velocity(real x, real y,
                    real2D u, real2D v, real2D ap,
                    interface in,
                    real *eu, real *ev,
                    int nx, int ny)
{
  real n1, n2;
  real a11, a22, a12, a21;
  real xmin, ymin, pmin, tmin;
  if (!closest_normal(x, y, in, &xmin, &ymin, &pmin, &tmin, &n1, &n2)) {
    fprintf(stderr, "extra_velocity(%g, %g): no closest_normal\n", x, y);
    exit(1);
  }
#ifdef EXTRA_LINEAR
  extra_velocity_normals(x, y, n1, n2, xmin, ymin,
                         u, v, ap, 
                         eu, ev, &a11, &a22, &a12, &a21,
                         nx, ny);
#else
  extra_velocity_normals2(x, y, n1, n2, xmin, ymin,
                          u, v, ap, 
                          eu, ev, &a11, &a22, &a12, &a21,
                          nx, ny);
#endif
}


real kappa_normals(polynom3 px, polynom3 py, real t, real *nx, real *ny)
{
  real xp, yp, y, dl, kappa;

  xp = splint1(px, t);
  yp = splint1(py, t);
  dl = sqrt(sq(xp) + sq(yp));
  kappa = (xp*splint2(py, t) - yp*splint2(px, t))/(sq(dl)*dl);
  /* axisymmetric term */
  y = splint(py, t) - 2.0;
  if (fabs(y) < 1e-2) /* second order approximation near the axis */
    kappa -= splint2(px, t)/splint1(py, t)/dl;
  else
    kappa -= xp/y/dl;
  *nx = - yp/dl;
  *ny = xp/dl;
  return kappa;
}


real avg_kappa_normals(polynom3 px, polynom3 py, real t1, real t2,
                       real *nx, real *ny)
{
  static real coef[7] = {1., 3., 3., 2., 3., 3., 1.}; 
  int i;
  real xp, yp, y, h = (t2 - t1)/6., t, kappasum = 0.0, lsum = 0.0, dl2;

  for (i = 0; i < 7; i++) {
    t = t1 + h*i;
    xp = splint1(px, t);
    yp = splint1(py, t);
    dl2 = sq(xp) + sq(yp);
    kappasum += coef[i]*((xp*splint2(py, t) - yp*splint2(px, t))/dl2);
    /* axisymmetric term */
    y = splint(py, t) - 2.0;
    if (fabs(y) < 1e-2) /* second order approximation near the axis */
      kappasum -= coef[i]*splint2(px, t)/splint1(py, t);
    else
      kappasum -= coef[i]*xp/y;
    lsum += coef[i]*sqrt(dl2);
  }
  *nx = 8.*(splint(py, t2) - splint(py, t1))/(3.*h*lsum);
  *ny = - 8.*(splint(px, t2) - splint(px, t1))/(3.*h*lsum);
  return kappasum/lsum;
}


/* -------------------------------------------------------------------------
                Resoud un systeme d'equations lineaires
                du type A * C = B par la methode du pivot de Gauss
   ------------------------------------------------------------------------- */
#define GAUSSTOL 1.e-10

static int gauss(real A[7][7], real B[7], real C[7], 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 (fabs(A[j][j]) < GAUSSTOL)
        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];
    }
  }

  if (fabs(A[n-1][n-1]) < GAUSSTOL)
    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];
    if (fabs(A[i][i]) < GAUSSTOL)
      return 1;
    C[i] /= A[i][i];
  }
  return 0;
}


int gauss13(real A[13][13], real B[13], real C[13], real *d, int n)
{
  int i, j, k, i1;
  real max, temp, det = 1.0;
  
  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 (fabs(A[j][j]) < 1e-10)
        return 1;
      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];
    }
  }

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

  if (fabs(A[n-1][n-1]) < 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];
    if (fabs(A[i][i]) < 1.e-10)
      return 1;
    C[i] /= A[i][i];
  }
  return 0;
}


int gauss2(real **A, real *B, real *C, 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];
    }
  }
  if (A[n - 1][n - 1] == 0.)
    return 1;
  else
    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;
}

---