---

advect.c

Go to the documentation of this file.

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

#define OUTU(i,j) (u[i][j] == UNDEFINED)
#define OUTV(i,j) (v[i][j] == UNDEFINED)

/*cld 
static real bilinu(real2D u, real x, real y)
cld*/
real bilinu(real2D u, real x, real y)
{
  int i, j;
  i = x; j = (int) ((real)y - 0.5);
  x -= (real)i;
  y -= 0.5 + (real)j;
  if (OUTU(i,j)) {
    x = 1. - x; y = 1. - y;
    if (x + y > 1.0 || OUTU(i,j+1) || OUTU(i+1,j) || OUTU(i+1,j+1))
      return UNDEFINED;
    return x*u[i][j+1] + y*u[i+1][j] + (1. - x - y)*u[i+1][j+1];
  }
  if (OUTU(i+1,j+1)) {
    if (x + y > 1.0 || OUTU(i+1,j) || OUTU(i,j+1))
      return UNDEFINED;
    return x*u[i+1][j] + y*u[i][j+1] + (1. - x - y)*u[i][j];
  }
  if (OUTU(i+1,j)) {
    y = 1. - y;
    if (x + y > 1.0 || OUTU(i,j+1))
      return UNDEFINED;
    return x*u[i+1][j+1] + y*u[i][j] + (1. - x - y)*u[i][j+1];
  }
  if (OUTU(i,j+1)) {
    x = 1. - x;
    if (x + y > 1.0)
      return UNDEFINED;
    return x*u[i][j] + y*u[i+1][j+1] + (1. - x - y)*u[i+1][j];
  }
  return u[i][j]*(1. - x - y + x*y) + u[i+1][j]*x*(1. - y) +
    u[i][j+1]*y*(1. - x) + u[i+1][j+1]*x*y;
}


real bilinv(real2D v, real x, real y)
{
  int i, j;
  i = (int)((real) x - 0.5); j = y;
  x -= 0.5 + (real)i;
  y -= (real) j;  
  if (OUTV(i,j)) {
    x = 1. - x; y = 1. - y;
    if (x + y > 1.0 || OUTV(i,j+1) || OUTV(i+1,j) || OUTV(i+1,j+1))
      return UNDEFINED;
    return x*v[i][j+1] + y*v[i+1][j] + (1. - x - y)*v[i+1][j+1];
  }
  if (OUTV(i+1,j+1)) {
    if (x + y > 1.0 || OUTV(i,j+1) || OUTV(i+1,j))
      return UNDEFINED;      
    return x*v[i+1][j] + y*v[i][j+1] + (1. - x - y)*v[i][j];
  }
  if (OUTV(i+1,j)) {
    y = 1. - y;
    if (x + y > 1.0 || OUTV(i,j+1))
      return UNDEFINED;
    return x*v[i+1][j+1] + y*v[i][j] + (1. - x - y)*v[i][j+1];
  }
  if (OUTV(i,j+1)) {
    x = 1. - x;
    if (x + y > 1.0)
      return UNDEFINED;
    return x*v[i][j] + y*v[i+1][j+1] + (1. - x - y)*v[i+1][j];
  }
  return v[i][j]*(1. - x - y + x*y) + v[i+1][j]*x*(1. - y) +
    v[i][j+1]*y*(1. - x) + v[i+1][j+1]*x*y;
}


void advect(interface in, real2D u, real2D v, real2D ap, int nx, int ny, 
            real *puaxis)
{
  int l;
  real x, y, dx, dy, n1, n2;
  real a11, a22, a12, a21;
//  FILE *fptr = fopen("markspeed", "at"); 

  for (l = 0; l < in.n; l++) {
    x = in.x[l]; y = in.y[l];
    if ((dx = bilinu(u, x, y)) == UNDEFINED ||
        (dy = bilinv(v, x, y)) == UNDEFINED) {
      printf("advect(): extrapolating point %g,%g ... ", x, y);
      if (l == in.n - 1) {
        n1 = - splint1(in.spy[l-1], in.t[l]);
        n2 = splint1(in.spx[l-1], in.t[l]);
      } 
      else {
        n1 = - splint1(in.spy[l], in.t[l]);
        n2 = splint1(in.spx[l], in.t[l]);
      }
      a11 = sqrt(sq(n1) + sq(n2));
      n1 /= a11; n2 /= a11;
      printf("normal(%g,%g) ... ", n1, n2);
      extra_velocity_normals(x, y, n1, n2, x, y,
                             u, v, ap, &dx, &dy,
                             &a11, &a22, &a12, &a21, nx, ny);
      printf("Ok\n");
    }
#ifdef CORRECTOR
    x = in.x[l] + dx; y = in.y[l] + dy;
    if ((dx = bilinu(u, x, y)) == UNDEFINED ||
        (dy = bilinv(v, x, y)) == UNDEFINED) {
      printf("advect(): extrapolating point %g,%g ... ", x, y);
      if (l == in.n - 1) {
        n1 = - splint1(in.spy[l-1], in.t[l]);
        n2 = splint1(in.spx[l-1], in.t[l]);
      } 
      else {
        n1 = - splint1(in.spy[l], in.t[l]);
        n2 = splint1(in.spx[l], in.t[l]);
      }
      a11 = sqrt(sq(n1) + sq(n2));
      n1 /= a11; n2 /= a11;
      printf("normal(%g,%g) ... ", n1, n2);
      extra_velocity_normals(x, y, n1, n2, x, y,
                             u, v, ap, &dx, &dy,
                             &a11, &a22, &a12, &a21, nx, ny);
      printf("Ok\n");
    }
#endif
    in.x[l] += dx;
    in.y[l] += dy;
//    printf("%d %g %g\n", l, dx, dy);

/*cld*/
//    if (l == 0)
//          *puaxis = dx;
/*cld*/
//        fprintf(fptr, "%g %g\n", in.t[l]/in.t[in.n-1], sqrt(sq(dx) + sq(dy))); 
//        fprintf(fptr, "%g %g %g\n", in.t[l]/in.t[in.n-1], dx, dy); 
  }  
  
//  fprintf(fptr, "&\n");
//  fclose(fptr);
  
}


static real *coord;
static int cmp(const void *i1, const void *i2)
{
  if (coord[*((int *)i1)] > coord[*((int *)i2)]) return 1;
  if (coord[*((int *)i1)] < coord[*((int *)i2)]) return -1;
  return 0;
}


static int find_smaller(real *coord, int *index, int n, real c) {
  int i, imin = 0, imax = n - 1;
  if (c < coord[index[0]]) return -1;
  if (c > coord[index[n-1]]) return n-1;
  i = (imin + imax)/2;
  while (coord[index[i]] > c || coord[index[i+1]] < c) {
    if (coord[index[i]] > c) imax = i;
    if (coord[index[i]] < c) imin = i;
    i = (imin + imax)/2;
  }
  return i;
}


static int find_larger(real *coord, int *index, int n, real c) {
  int i, imin = 0, imax = n - 1;
  if (c > coord[index[n-1]]) return -1;
  if (c < coord[index[0]]) return 0;
  i = (imin + imax)/2;
  while (coord[index[i]] < c || coord[index[i-1]] > c) {
    if (coord[index[i]] > c) imax = i;
    if (coord[index[i]] < c) imin = i;
    i = (imin + imax)/2;
  }
  return i;
}


void repulse(interface in, real cutoff, real coef, real lredis)
{
  int i, j, imin, imax, imin1, imax1, index;
  real xm, ym, xi, yi, t, r2;
  real *fx = (real *)calloc(in.n, sizeof(real));
  real *fy = (real *)calloc(in.n, sizeof(real));
  real *abx = (real *)malloc((in.n - 1)*sizeof(real));
  real *aby = (real *)malloc((in.n - 1)*sizeof(real));
  real *abn = (real *)malloc((in.n - 1)*sizeof(real));
#if 1
  int *isx = (int *)malloc(in.n*sizeof(int));
  int *isy = (int *)malloc(in.n*sizeof(int));
  int *is;
FILE *fptr = fopen("repulse", "wt");

  cutoff *= cutoff;
  for (i = 0; i < in.n - 1; i++) {
    abx[i] = in.x[i+1] - in.x[i]; 
    aby[i] = in.y[i+1] - in.y[i];
    abn[i] = sqrt(sq(abx[i]) + sq(aby[i]));
  }

  /* sort the markers by ascending coordinates */
  for (i = 0; i < in.n; i++)
    isx[i] = isy[i] = i;
  coord = in.x;
  qsort(isx, in.n, sizeof(int), cmp);
  coord = in.y;
  qsort(isy, in.n, sizeof(int), cmp);

  for (i = 0; i < in.n; i++) {
    xm = in.x[i]; ym = in.y[i];
    /* find the points with an x-coordinate in [xm-lredis,xm+lredis] */
    imax = find_smaller(in.x, isx, in.n, xm + lredis);
    imin = find_larger(in.x, isx, in.n, xm - lredis);
    /* find the points with an y-coordinate in [ym-lredis,ym+lredis] */
    imax1 = find_smaller(in.y, isy, in.n, ym + lredis);
    imin1 = find_larger(in.y, isy, in.n, ym - lredis);
    if (imax - imin > imax1 - imin1) { imax = imax1; imin = imin1; is = isy; }
    else { is = isx; }
    /* for each point in the smallest interval */
    for (j = imin; j <= imax; j++) {
      index = is[j];
      if (index != i && index+1 != i) {
        if (index < in.n - 1) {
          t = ((xm - in.x[index])*abx[index] +
               (ym - in.y[index])*aby[index])/abn[index];
          if (t <= 1.0) {
            if (t <= 0.0) { xi = in.x[index]; yi = in.y[index]; }
            else {
              xi = in.x[index] + t*abx[index];
              yi = in.y[index] + t*aby[index];
            }
            r2 = sq(xm - xi) + sq(ym - yi);
            if (r2 < cutoff) {
              r2 = coef*1./sq(r2);
              fx[i] -= r2*(xi - xm);
              fy[i] -= r2*(yi - ym);
              fprintf(fptr, "%g %g\n%g %g\n&\n\n", xi, yi, xm, ym);
            }
          }
        }
      }
    }
  }
  in.x[0] += fx[0]; in.x[in.n-1] += fx[in.n-1];
  if (in.bc != AXITRUE && in.bc != AXINATURAL) {
    in.y[0] += fy[0]; in.y[in.n-1] += fy[in.n-1];
  }
  for (i = 1; i < in.n - 1; i++) {
    in.x[i] += fx[i]; in.y[i] += fy[i];
  }
  free(fx); free(fy); 
  free(abx); free(aby); free(abn);
  free(isx); free(isy);
  fclose(fptr);
#else
FILE *fptr = fopen("repulse", "wt");

  cutoff *= cutoff;
  for (i = 0; i < in.n - 1; i++) {
    abx[i] = in.x[i+1] - in.x[i]; 
    aby[i] = in.y[i+1] - in.y[i];
    abn[i] = sqrt(sq(abx[i]) + sq(aby[i]));
  }
  for (i = 0; i < in.n; i++) {
    xm = in.x[i]; ym = in.y[i];
    for (j = 0; j < in.n - 1; j++) 
      if (j != i && j+1 != i) {
        t = ((xm - in.x[j])*abx[j] + (ym - in.y[j])*aby[j])/abn[j];
        if (t <= 1.0) {
          if (t < 0.0) {
            xi = in.x[j]; yi = in.y[j];
          }
          else {
            xi = in.x[j] + t*abx[j];
            yi = in.y[j] + t*aby[j];
          }
          r2 = sq(xm - xi) + sq(ym - yi);
          if (r2 < cutoff) {
            r2 = coef*1./sq(r2);
            fx[i] -= r2*(xi - xm);
            fy[i] -= r2*(yi - ym);
fprintf(fptr, "%g %g\n%g %g\n&\n\n", xi, yi, xm, ym);
          }
        }
      }
  }
  fclose(fptr);

  in.x[0] += fx[0]; in.x[in.n-1] += fx[in.n-1];
  if (in.bc != AXITRUE && in.bc != AXINATURAL) {
    in.y[0] += fy[0]; in.y[in.n-1] += fy[in.n-1];
  }
  for (i = 1; i < in.n - 1; i++) {
    in.x[i] += fx[i]; in.y[i] += fy[i];
  }
  free(fx); free(fy); 
  free(abx); free(aby); free(abn);
#endif
}

---