---

surfaces.poly.c

Go to the documentation of this file.

#include "utilf.h"
#include "markers.h"
#include "surfaces.h"
#include "polynomial.h"

#define MAXROUND 0.00001


real interface_inside(interface in, real x, real y)
{
  int i;
  real a, d, dmin = 1e30;
  for (i = 0; i < in.n - 1; i++) {
    d = sqrt(sq(x - 0.5*(in.x[i] + in.x[i+1])) + 
             sq(y - 0.5*(in.y[i] + in.y[i+1])));
    a = (in.x[i+1] - in.x[i])*(y - in.y[i]) - 
      (x - in.x[i])*(in.y[i+1] - in.y[i]);
    if (d < fabs(dmin)) dmin = a >= 0.0 ? d : -d;
  }
  if (dmin >= 0.0) return dmin;
  return 0.0;
}


void interface_fill(interface in, real2D ap, int nx, int ny)
{
  int i, j;

  fill0(ap, nx, ny);
  for (i = 2; i <= nx; i++)
    for (j = 2; j <= ny; j++)
      if (interface_inside(in, (real)i, (real)j)) {
        ap[i][j] += 0.25; ap[i-1][j] += 0.25;
        ap[i][j-1] += 0.25; ap[i-1][j-1] += 0.25;
      }
  for (j = 2; j <= ny; j++)
    if (interface_inside(in, 1.0, (real)j)) {
      ap[1][j] += 0.25; ap[1][j-1] += 0.25;
    }
  for (i = 2; i <= nx; i++)
    if (interface_inside(in, (real)i, 1.0)) {
      ap[i][1] += 0.25; ap[i-1][1] += 0.25;
    }
}


void cut_interface(interface in, interface_cut *cut, real xo, real yo)
{
  int i, j, k, l;
  int n, nin;
  real ti[14], t1, t2;
  polynom3 px, py;
  listinter listi[14], xl;

  if (cut->list == NULL) cut->list = (listinter *)malloc(sizeof(listinter));
  cut->ncut = 0;
  cut->xo = xo; cut->yo = yo;

  /* For each segment */
  for (l = 0; l < in.n - 1; l++) {
    /* the positions are shifted (-xo, -yo) */
    i = (real)(in.x[l] - xo); 
    j = (real)(in.y[l] - yo);
    t1 = in.t[l]; t2 = in.t[l+1];
    px = in.spx[l]; py = in.spy[l];
    px.a -= xo; py.a -= yo;
    
    /* first point */
    listi[0].t = t1; listi[0].type = MARK_TYPE;
    nin = 1;
    
    /* intersections with the horizontals */
    n = polyroots(py, t1, t2, ti, (real) j, INTERPREC);
    for (k = 0; k < n; k++) {
      listi[nin].t = ti[k]; listi[nin].type = HOR_TYPE; 
      listi[nin++].coord = j;
    }
    n = polyroots(py, t1, t2, ti, (real) j + 1.0, INTERPREC);
    for (k = 0; k < n; k++) {
      listi[nin].t = ti[k]; listi[nin].type = HOR_TYPE; 
      listi[nin++].coord = j + 1;
    }
    /* intersections with the verticals */
    n = polyroots(px, t1, t2, ti, (real) i, INTERPREC);
    for (k = 0; k < n; k++) {
      listi[nin].t = ti[k]; listi[nin].type = VER_TYPE; 
      listi[nin++].coord = i;
    }
    n = polyroots(px, t1, t2, ti, (real) i + 1.0, INTERPREC);
    for (k = 0; k < n; k++) {
      listi[nin].t = ti[k]; listi[nin].type = VER_TYPE; 
      listi[nin++].coord = i + 1;
    }
    
    /* last point */
    listi[nin].t = t2; listi[nin++].type = MARK_TYPE;
    
    /* sort the intersections in ascending order for t */
    for (k = 1; k < nin; k++) {
      xl = listi[k]; i = k - 1;
      while (i >= 0 && listi[i].t > xl.t) {
        listi[i + 1] = listi[i]; i--;
      }
      listi[i + 1] = xl;
    }

    /* add the points to the list */
    cut->list = (listinter *)realloc(cut->list,
                                     (cut->ncut + nin)*sizeof(listinter));
    for (k = 0; k < nin; k++) {
      listi[k].ip = l;
      cut->list[cut->ncut++] = listi[k];
    }
  }
}


void interface_cut_write(interface in, interface_cut cut, char *file)
{
  int i, ip;
  polynom3 px, py;
  real t;
  FILE *fptr = fopen(file, "wt");

  for (i = 0; i < cut.ncut; i++)
    if (cut.list[i].type == MARK_TYPE) {
      ip = cut.list[i].ip;
      px = in.spx[ip]; py = in.spy[ip];
      t = cut.list[i].t;
      fprintf(fptr, "%g %g\n", splint(px, t), splint(py, t));
    }
  fprintf(fptr, "&\n");
  for (i = 0; i < cut.ncut; i++)
    if (cut.list[i].type == HOR_TYPE) {
      ip = cut.list[i].ip;
      px = in.spx[ip]; py = in.spy[ip];
      t = cut.list[i].t;
      fprintf(fptr, "%g %g\n", splint(px, t), splint(py, t));
    }
  fprintf(fptr, "&\n");
  for (i = 0; i < cut.ncut; i++)
    if (cut.list[i].type == VER_TYPE) {
      ip = cut.list[i].ip;
      px = in.spx[ip]; py = in.spy[ip];
      t = cut.list[i].t;
      fprintf(fptr, "%g %g\n", splint(px, t), splint(py, t));
    }
  fclose(fptr);
}


void interface_surfaces(interface in, interface_cut cut,
                        real2D sx, real2D sy, real2D a, real2D c,
                        real2D wa, real2D wc,
                        int nx, int ny)
{
  int i, j, l;
  real t, t1, t2;
  real x, y, cmax;
  polynom3 px, py;
  real dy[3], xdy[6], ixdy[7], xydy[9], ixydy[10];

  fill0(wa, nx, ny); fill0(wc, nx, ny);
  fill0(sx, nx, ny); fill0(sy, nx, ny);

  /* For each intersection */
  for (l = 0; l < cut.ncut - 1; l++) {
    px = in.spx[cut.list[l].ip];
    py = in.spy[cut.list[l].ip];
    px.a -= cut.xo;
    py.a -= cut.yo;

    t1 = cut.list[l].t; t2 = cut.list[l+1].t;
    t = 0.5*(t1 + t2);
    i = splint(px, t);
    j = splint(py, t);
    
    x = splint(px, t1);
    y = splint(py, t1);

    px.a -= (real)i;
    polyderiv((real *)&py, 4, dy);
    polymul((real *)&px, 4, dy, 3, xdy);
    polyint(xdy, 6, ixdy);

    py.a += cut.yo - 2.;
    polymul((real *)&py, 4, xdy, 6, xydy);
    polyint(xydy, 9, ixydy);
    
    /* Intersection with the verticals */
    if (cut.list[l].type == VER_TYPE) {
      /* deal with roundoff errors */
      if (y > (real)j + 1.0) y = (real)j + 1.0; 
      else if (y < (real)j) y = (real)j;
      /* add the circulation of the vertical segment */
      if (cut.list[l].coord == i) {
        wa[i-1][j] += (real)j + 1.0 - y;
        t = 0.5*(sq((real)j + cut.yo - 1.) - sq(y + cut.yo - 2.));
        wc[i-1][j] += t;
        sy[i][j] += t;
      } 
      else {
        wa[i][j] += y - (real)j;
        t = 0.5*(sq(y + cut.yo - 2.) - sq((real)j + cut.yo - 2.));
        wc[i][j] += t;
        sy[i+1][j] += t;
      }
    }
    /* Intersection with the horizontals */
    else if (cut.list[l].type == HOR_TYPE) {
      /* deal with roundoff errors */
      if (x > (real)i + 1.0) x = (real)i + 1.0; 
      else if (x < (real)i) x = (real)i;
      /* add the circulation of the horizontal segment */
      if (cut.list[l].coord == j)
        sx[i][j] += (x - (real)i)*((real)j + cut.yo - 2.);
      else
        sx[i][j+1] += ((real)i + 1. - x)*((real)j + cut.yo - 1.);
    }
    
    /* add the circulation along the polynomial segment */
    wa[i][j] += polyeval(ixdy, 7, t2) - polyeval(ixdy, 7, t1);
    wc[i][j] += polyeval(ixydy, 10, t2) - polyeval(ixydy, 10, t1);
  }

  /* merdic formula */
  for (i = 1; i <= nx; i++)
    for (j = 2; j <= ny; j++) {
      t = wa[i][j] = modf(wa[i][j], &y);
      if (wa[i][j] >= a[i][j] + 0.5) wa[i][j] -= 1.0;
      if (wa[i][j] <= a[i][j] - 0.5) wa[i][j] += 1.0;
      if (wa[i][j] <= a[i][j] - 0.5) wa[i][j] += 1.0;

      cmax = (real)j + cut.yo - 1.5;
      wc[i][j] = cmax == 0.0 ? 0.0 : cmax*modf(wc[i][j]/cmax, &y);
      if (wc[i][j] >= (a[i][j] + 0.5)*cmax) wc[i][j] -= cmax;
      if (wc[i][j] <= (a[i][j] - 0.5)*cmax) wc[i][j] += cmax;
      if (wc[i][j] <= (a[i][j] - 0.5)*cmax) wc[i][j] += cmax;

      if ((cut.yo != -0.5 || j != 2) && (wa[i][j] < - MAXROUND
                                         || wa[i][j] > 1.0 + MAXROUND))
        printf("WARNING wa[%d][%d]: %g %g %g xo: %g yo: %g\n", 
               i, j, wa[i][j], a[i][j], t, cut.xo, cut.yo);
      if ((cut.yo != -0.5 || j != 2) && (wc[i][j] < - MAXROUND*cmax 
                                         || wc[i][j] > cmax*(1. + MAXROUND)))
        printf("WARNING wc[%d][%d]: %g c[%d][%d]:%g cmax: %g xo,yo: %g,%g\n", i, j, wc[i][j], i, j, c[i][j], cmax, cut.xo, cut.yo);
      c[i][j] = wc[i][j];
      a[i][j] = wa[i][j];
    }
}


void sx_fill(real2D sx, real2D a, real xo, real yo, 
             interface in, int nx, int ny)
{
  int i, j;
  real sxmax;

  for (i = 2; i <= nx; i++)
    for (j = 3; j <= ny; j++) {
      sxmax = (real)j + yo - 2.;
      if (sx[i][j] > sxmax) sx[i][j] -= sxmax;
      if (sx[i][j] == 0.0 && (a[i][j] > 0.0 || a[i][j-1] > 0.0))
        if (a[i][j] == 1.0 || a[i][j-1] == 1.0 || 
            interface_inside(in, (real)i + xo + 0.5, (real)j + yo))
          sx[i][j] = sxmax;
      if (sx[i][j] < 0.0 || sx[i][j] > sxmax)
        printf("WARNING sx[%d][%d]:%g sxmax: %g xo,yo: %g,%g\n", 
               i, j, sx[i][j], sxmax, xo, yo);
    }
}


void sy_fill(real2D sy, real2D a, real xo, real yo, 
             interface in, int nx, int ny)
{
  int i, j;
  real symax;

  for (i = 2; i <= nx; i++)
    for (j = 2; j <= ny; j++) {
      symax = (real)j + yo - 1.5;
      if (sy[i][j] > symax) sy[i][j] -= symax;
      if (sy[i][j] == 0.0 && (a[i][j] > 0.0 || a[i-1][j] > 0.0))
        if (a[i][j] == 1.0 || a[i-1][j] == 1.0 || 
            interface_inside(in, (real)i + xo, (real)j + yo + 0.5))
          sy[i][j] = symax;
      if (symax > 0.0 && (sy[i][j] < 0.0 || sy[i][j] > symax))
        printf("WARNING sy[%d][%d]:%g symax: %g xo,yo: %g,%g\n", 
               i, j, sy[i][j], symax, xo, yo);
    }
}


void sx_write(real2D sx, real xo, real yo, int nx, int ny, FILE *fptr)
{
  int i, j;
  for (i = 1; i <= nx; i++)
    for (j = 3; j <= ny; j++)
      fprintf(fptr, "%g %g\n%g %g\n\n", 
              (real)i + 0.5 + xo - 0.5*sx[i][j]/((real)j + yo - 1.5), 
              (real)j + yo,
              (real)i + 0.5 + xo + 0.5*sx[i][j]/((real)j + yo - 1.5), 
              (real)j + yo);
}


void sy_write(real2D sy, real xo, real yo, int nx, int ny, FILE *fptr)
{
  int i, j;
  for (i = 1; i <= nx; i++)
    for (j = 2; j <= ny; j++)
      fprintf(fptr, "%g %g\n%g %g\n\n", 
              (real)i + xo, 
              (real)j + 0.5 + yo - 0.5*sy[i][j]/((real)j + yo - 1.5), 
              (real)i + xo,
              (real)j + 0.5 + yo + 0.5*sy[i][j]/((real)j + yo - 1.5));
}


void interface_rzdr(interface in, interface_cut cut, real za,
                    real2D rzdr,
                    int nx, int ny)
{
  int i, j, l;
  real t, t1, t2;
  polynom3 px, py;
  real dr[3], rdr[6], rzdrp[9], irzdr[10];

  fill0(rzdr, nx, ny);

  /* For each intersection */
  for (l = 0; l < cut.ncut - 1; l++) {
    px = in.spx[cut.list[l].ip];
    py = in.spy[cut.list[l].ip];
    px.a -= cut.xo;
    py.a -= cut.yo;

    t1 = cut.list[l].t; t2 = cut.list[l+1].t;
    t = 0.5*(t1 + t2);
    i = splint(px, t);
    j = splint(py, t);

    py.a += cut.yo - 2.; 
    polyderiv((real *)&py, 4, dr);
    polymul((real *)&py, 4, dr, 3, rdr);
    px.a -= (real)i + 0.5 + za;
    polymul((real *)&px, 4, rdr, 6, rzdrp);
    polyint(rzdrp, 9, irzdr);
    
    /* add the circulation along the polynomial segment */
    rzdr[i][j] += polyeval(irzdr, 10, t2) - polyeval(irzdr, 10, t1);
  }
}


void interface_rrdz(interface in, interface_cut cut, real ra,
                    real2D rrdz,
                    int nx, int ny)
{
  int i, j, l;
  real t, t1, t2;
  polynom3 px, py;
  real dz[3], rdz[6], rrdzp[9], irrdz[10];

  fill0(rrdz, nx, ny);

  /* For each intersection */
  for (l = 0; l < cut.ncut - 1; l++) {
    px = in.spx[cut.list[l].ip];
    py = in.spy[cut.list[l].ip];
    px.a -= cut.xo;
    py.a -= cut.yo;

    t1 = cut.list[l].t; t2 = cut.list[l+1].t;
    t = 0.5*(t1 + t2);
    i = splint(px, t);
    j = splint(py, t);
    
    py.a += cut.yo - 2.; px.a += cut.xo;
    polyderiv((real *)&px, 4, dz);
    polymul((real *)&py, 4, dz, 3, rdz);
    py.a -= cut.yo + (real)j - 1.5 + ra;
    polymul((real *)&py, 4, rdz, 6, rrdzp);
    polyint(rrdzp, 9, irrdz);
    
    /* add the circulation along the polynomial segment */
    rrdz[i][j] += polyeval(irrdz, 10, t2) - polyeval(irrdz, 10, t1);
  }
}

---