---

surfaces.c

Go to the documentation of this file.

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

#define MAXROUND 0.00001


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;
    }
    
    /* 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];
    }
  }
  /* last point */
  xl.t = in.t[in.n-1]; xl.type = MARK_TYPE; xl.ip = in.n - 1;
  cut->list = (listinter *)realloc(cut->list,
                                   (cut->ncut + 1)*sizeof(listinter));
  cut->list[cut->ncut++] = xl;
}


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 x1, x2, y1, y2;

  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);

    /* add the circulations along the interface segment */
    px.a -= (real)i; py.a += cut.yo - 2.;
    x1 = splint(px, t1); x2 = splint(px, t2);
    y1 = splint(py, t1); y2 = splint(py, t2);

    wa[i][j] += (y2 - y1)*(x1 + x2)/2.;
    wc[i][j] += (y2 - y1)*((x2 - x1)*(y2 - y1)/3. + 
                           ((x2 - x1)*y1 + (y2 - y1)*x1)/2. + x1*y1);
    
    /* 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.);
    }
  }


#if 0
  /* 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) wa[i][j] -= 1.0;
      if (wa[i][j] <= a[i][j] - 0.5 && wa[i][j]<=0.0) wa[i][j] += 1.0;
      if (wa[i][j] <= a[i][j] - 0.5 && wa[i][j]<=0.0) 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) wc[i][j] -= cmax;
      if (wc[i][j] <= (a[i][j] - 0.5)*cmax && wc[i][j]<=0.0) wc[i][j] += cmax;
      if (wc[i][j] <= (a[i][j] - 0.5)*cmax && wc[i][j]<=0.0) wc[i][j] += cmax;

      if ((cut.yo != -0.5 || (j != 2 && j != ny)) && (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 && j != ny)) && (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] > cmax ? cmax : wc[i][j] < 0.0 ? 0.0 : wc[i][j];
      a[i][j] = wa[i][j] > 1.0 ? 1.0 : wa[i][j] < 0.0 ? 0.0 : wa[i][j];
    }

#else

  /* 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;

/*      wa[i][j] = wa[i][j] > 1.0 ? 1.0 : wa[i][j] < 0.0 ? 0.0 : wa[i][j]; */

      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;

/*      wc[i][j] = wc[i][j] > cmax ? cmax : wc[i][j] < 0.0 ? 0.0 : wc[i][j]; */

      if ((cut.yo != -0.5 || (j != 2 && j != ny)) && (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 && j != ny)) && (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] > cmax ? cmax : wc[i][j] < 0.0 ? 0.0 : wc[i][j];
      a[i][j] = wa[i][j] > 1.0 ? 1.0 : wa[i][j] < 0.0 ? 0.0 : wa[i][j];
    }
#endif
}


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 x1, y1, x2, y2;

  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);

    /* add the circulation along the interface segment */
    py.a += cut.yo - 2.; 
    px.a -= (real)i + 0.5 + za;
    x1 = splint(px, t1); x2 = splint(px, t2);
    y1 = splint(py, t1); y2 = splint(py, t2);
    rzdr[i][j] += (y2 - y1)*((x2 - x1)*(y2 - y1)/3. + 
                             ((x2 - x1)*y1 + (y2 - y1)*x1)/2. + x1*y1);
  }
}


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 x1, y1, x2, y2;

  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);
    
    /* add the circulation along the polynomial segment */
    py.a += cut.yo - 2.;
    x1 = splint(px, t1); x2 = splint(px, t2);
    y1 = splint(py, t1); y2 = splint(py, t2);
    t = cut.yo + (real)j - 1.5 + ra;
    rrdz[i][j] +=
      (x2 - x1)*(sq(y2 - y1)/3. + (2.*y1 - t)*(y2 - y1)/2. + y1*(y1 - t));
  }
}

---