---

fluxes.c

Go to the documentation of this file.

#include "utilf.h"
#include "markers.h"
#include "fluxes.h"
#include "inout.h"


void interface_fluxes(interface in, real2D cc,
                      real2D s1, real2D s2, real2D s3, real2D s4,
                      int nx, int ny)
{
  int i, j, k, l;
  int n, nin;
  real ti[14], t, t1, t2;
  real x, y, x1, y1, yo;
  polynom3 px, py;
  typedef struct {
    real t;
    int type, coord;
  } listinter;
  listinter listi[14], xl;

  fill0(s1, nx, ny);  fill0(s3, nx, ny);  
  fill0(s2, nx, ny);  fill0(s4, nx, ny);
#if 0
  /* For each segment */  
  for (l = 0; l < in.n - 1; l++) {
    i = (real)in.x[l]; 
    j = (real)in.y[l];
    t1 = in.t[l]; t2 = in.t[l+1];
    px = in.spx[l]; py = in.spy[l];
    
    /* 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;
    }
    
    for (k = 0; k < nin - 1; k++)  {
      t1 = listi[k].t; t2 = listi[k+1].t;
      t = 0.5*(t1 + t2);
      i = splint(px, t);
      j = splint(py, t);
      
      x = splint(px, t1); x1 = splint(px, t2); 
      y = splint(py, t1) - 2.; y1 = splint(py, t2) - 2.;
      yo = (real)j - 2.;
      
      /* Intersection with the verticals */
      if (listi[k].type == VER_TYPE) {
        /* add the circulation of the vertical segment */
        if (listi[k].coord == i) {
          t = (1. + yo - y)*(1. + yo + y)/2.;
          s3[i][j] -= t;
          s1[i-1][j] += t;
        } 
        else {
          t = (y - yo)*(y + yo)/2.;
          s1[i][j] += t;
          s3[i+1][j] -= t;
        }
      }

      /* Intersection with the horizontals */
      if (listi[k].type == HOR_TYPE) {
        /* add the circulation of the horizontal segment */
        if (listi[k].coord == j) {
          t = (x - i)*yo;
          s4[i][j] -= t;
          s2[i][j-1] += t;
        } 
        else {
          t = (1. + i - x)*(yo + 1.);
          s2[i][j] += t;
          s4[i][j+1] -= t;
        }
      }

      /* add the circulation along the segment */
      x -= i; x1 -= i;
      t = (y1 - y)*((x1*y1 + x*y)/3. + (x1*y + x*y1)/6.);
      s1[i][j] += t;
      s3[i][j] -= t - (y1 - y)*(y1 + y)/2.;
      t = (sq(y1) + y1*y + sq(y))/3.;
      s2[i][j] += (x1 - x)*(yo*(y1 + y)/2. - t);
      s4[i][j] -= (x1 - x)*((yo + 1.)*(y1 + y)/2. - t);
    }
  }
#endif
  for (i = 2; i < nx; i++)
    for (j = 2; j < ny; j++)
      if (INP(i,j)) {
        yo = (real)j - 2.;
        if (s1[i][j] <= 0.1*yo) s1[i][j] += yo + 0.5;
        else if (s1[i][j] > yo + 0.5) s1[i][j] -= yo + 0.5;
        if (s3[i][j] >= - 0.1*yo) s3[i][j] -= yo + 0.5;
        else if (s3[i][j] < - yo - 0.5) s3[i][j] += yo + 0.5;
        if (s2[i][j] <= 0.1*yo) s2[i][j] += yo + 1.;
        else if (s2[i][j] > yo + 1.) s2[i][j] -= yo + 1.;
        if (s4[i][j] >= - 0.1*yo) s4[i][j] -= yo;
        else if (s4[i][j] < -yo) s4[i][j] += yo;

        if (s1[i][j] < 0.1*yo) { 
          s1[i][j] += yo + 0.5; 
          printf("s1[%d][%d]=%g <= 0.1\n", i, j, s1[i][j]);
        }
        else if (s1[i][j] > yo + 0.5) {
          s1[i][j] -= yo + 0.5;
          printf("s1[%d][%d]=%g >\n", i, j, s1[i][j]);
        }
        if (s3[i][j] > - 0.1*yo) { 
          s3[i][j] -= yo + 0.5;
          printf("s3[%d][%d]=%g >= - 0.1\n", i, j, s3[i][j]);
        }
        else if (s3[i][j] < - yo - 0.5) {
          s3[i][j] += yo + 0.5;
          printf("s3[%d][%d]=%g <= \n", i, j, s3[i][j]);
        }
        if (s2[i][j] < 0.1*yo) {
          s2[i][j] += yo + 1.;
          printf("s2[%d][%d]=%g <= 0.1\n", i, j, s2[i][j]);
        }
        else if (s2[i][j] > yo + 1.) {
          s2[i][j] -= yo + 1.;
          printf("s2[%d][%d]=%g >\n", i, j, s2[i][j]);
        }
        if (s4[i][j] > - 0.1*yo) {
          s4[i][j] -= yo;
          printf("s4[%d][%d]=%g >= - 0.1\n", i, j, s4[i][j]);
        }
        else if (s4[i][j] < -yo) {
          s4[i][j] += yo;
          printf("s4[%d][%d]=%g <\n", i, j, s4[i][j]);
        }
        /*      printf("%d,%d: s1: %g s2: %g s3: %g s4: %g\n", i, j, 
                s1[i][j], s2[i][j], s3[i][j], s4[i][j]); */
        if (fabs(2.*(s1[i][j] + s3[i][j])/(s1[i][j] - s3[i][j])) > 1e-3) {
          fprintf(stderr, "interface_fluxes(): s1[%d][%d] = %g != - s3[%d][%d] = %g\n",
                  i, j, s1[i][j], i, j, - s3[i][j]);
          exit(1);
        }
      }
}


void interface_hfrac(interface in, real2D fx, real2D c,
                     real xo, real yo, 
                     int nx, int ny)
{
  int i, j, k, l;
  int n;
  real ti[14], t, t1, t2;
  real x;
  polynom3 px, py;
FILE *fptr = fopen("hfrac", "wt");

 fill0(fx, nx, ny);

  /* For each segment */  
  for (l = 0; l < in.n - 1; l++) {
    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;
    
    /* intersections with the horizontals */
    n = polyroots(py, t1, t2, ti, (real) j, INTERPREC);
    for (k = 0; k < n; k++) {
      x = splint(px, ti[k]); i = x;
      fx[i][j] += (real)i + 1. - x;
    }
    n = polyroots(py, t1, t2, ti, (real) j + 1.0, INTERPREC);
    for (k = 0; k < n; k++) {
      x = splint(px, ti[k]); i = x;
      fx[i][j+1] += x - (real)i;
    }
  }

  for (i = 2; i < nx; i++)
    for (j = 2; j < ny; j++) {
      if (fx[i][j] > 1.0) fx[i][j] -= 1.0;
      if (fx[i][j] == 0.0 && c[i][j] + c[i][j-1] >= 1.0) fx[i][j] = 1.0;
      if (fx[i][j] > 1.0 || fx[i][j] < 0.0) {
        fprintf(stderr, "interface_hfrac(): fx[%d][%d]=%g\n", i, j, fx[i][j]);
        exit(1);
      }
      fprintf(fptr, "%g %g\n%g %g\n\n", 
              (real)i + xo + 0.5 + fx[i][j]/2., (real)j + yo,
              (real)i + yo + 0.5 - fx[i][j]/2., (real)j + yo);
    }
  fclose(fptr);
}


void interface_vfrac(interface in, real2D fy, real2D c,
                     real xo, real yo, int nx, int ny)
{
  int i, j, k, l;
  int n;
  real ti[14], t, t1, t2;
  real y;
  polynom3 px, py;
FILE *fptr = fopen("vfrac", "wt");

 fill0(fy, nx, ny);

  /* For each segment */  
  for (l = 0; l < in.n - 1; l++) {
    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;
    
    /* intersections with the verticals */
    n = polyroots(px, t1, t2, ti, (real) i, INTERPREC);
    for (k = 0; k < n; k++) {
      y = splint(py, ti[k]); j = y;
      fy[i][j] += 0.5*(sq(y + yo - 2.) - sq((real)j + yo - 2.))/
        ((real)j + yo - 1.5);
    }
    n = polyroots(px, t1, t2, ti, (real) i + 1.0, INTERPREC);
    for (k = 0; k < n; k++) {
      y = splint(py, ti[k]); j = y;
      fy[i+1][j] += 0.5*(sq((real)j + yo - 1.) - sq(y + yo - 2.))/
        ((real)j + yo - 1.5);
    }
  }

  for (i = 2; i < nx; i++)
    for (j = 2; j < ny; j++) {
      if (fy[i][j] > 1.0) fy[i][j] -= 1.0;
      if (fy[i][j] == 0.0 && c[i][j] + c[i-1][j] >= 1.0) fy[i][j] = 1.0;
      if (fy[i][j] > 1.0 || fy[i][j] < 0.0) {
        fprintf(stderr, "interface_vfrac(): fy[%d][%d]=%g\n", i, j, fy[i][j]);
        exit(1);
      }
      fprintf(fptr, "%g %g\n%g %g\n\n", 
              (real)i + xo, (real)j + yo + 0.5 + fy[i][j]/2.,
              (real)i + yo, (real)j + yo + 0.5 - fy[i][j]/2.);
    }
  fclose(fptr);
}

---