Main Page   Data Structures   File List   Data Fields   Globals  

markers.c File Reference

#include <malloc.h>
#include "utilf.h"
#include "markers.h"

Include dependency graph for markers.c:

Include dependency graph

Go to the source code of this file.

Defines

#define PERIODIC(i, nx, ip)   {ip=(i-2)%(nx-2)+2; if(ip<2) ip+=nx-2;}
#define PERIODIC_DL(dl, nx)

Functions

int interface_arc (interface in, real t)
void interface_write_markers (interface in, FILE *fptr)
void interface_write_markers_speed (real time, interface in, real2D u, real2D v, real tau, real du, int nx, int ny, FILE *fptr)
void interface_write_splines (interface in, FILE *fptr, int nb)
void curvature_write (interface in, FILE *fptr)
real interfaces_area (interface *in, int n)
real interfaces_axivolume (interface *in, int n)
real interfaces_xmoment (interface *in, int n)
real interfaces_ymoment (interface *in, int n)
real interfaces_x_amp (interface *in, int n)
real interfaces_y_amp (interface *in, int n)
void interface_splines (interface in, int nx, int ny)
void interface_redis (interface *in, real space)
void interface_redis_scale (interface *in, interface src, real space, real scale)
void interface_filter (interface *in, int nx, int ny)

Define Documentation

#define PERIODIC (  i,
nx,
ip    )     {ip=(i-2)%(nx-2)+2; if(ip<2) ip+=nx-2;}
 

Definition at line 6 of file markers.c.

#define PERIODIC_DL (  dl,
nx    ) 
 

Value:

{if(fabs(dl)>0.5*(nx-2))\
                              dl=dl>0.0?dl-nx+2.:dl+nx-2.;}

Definition at line 7 of file markers.c.

Referenced by interface_splines().

Function Documentation

void curvature_write (  interface    in,
FILE *    fptr
) 
 

Definition at line 68 of file markers.c.

References interface::n, real, splint1, splint2, interface::spx, interface::spy, sq, and interface::t.

Referenced by initialize(), and timestep(). 

00069 {
00070   int i;
00071   real x1, y1, dl, kappa;
00072 
00073   fprintf(fptr, "\n");
00074   for (i = 0; i < in.n - 1; i++) {
00075     x1 = splint1(in.spx[i], in.t[i]);
00076     y1 = splint1(in.spy[i], in.t[i]);    
00077     dl = sqrt(sq(x1) + sq(y1));
00078     kappa = (x1*splint2(in.spy[i], in.t[i]) - y1*splint2(in.spx[i], in.t[i]))
00079             / (dl*dl*dl);
00080     fprintf(fptr, "%g %g\n", in.t[i]/in.t[in.n-1], kappa);
00081   }
00082 }

int interface_arc (  interface    in,
real    t
) 
 

Definition at line 10 of file markers.c.

References interface::n, real, and interface::t.

Referenced by extra_poly(). 

00011 {
00012   int imin = 0, imax = in.n - 1, i;
00013   
00014   if (t < 0.0 || t > in.t[in.n - 1]) return -1;
00015   
00016   while(imax - imin > 1) {
00017     i = (imin + imax)/2;
00018     if (t == in.t[i]) return i;
00019     if (t > in.t[i]) imin = i;
00020     else imax = i;
00021   }
00022   return imin;
00023 }

void interface_filter (  interface *    in,
int    nx,
int    ny
) 
 

Definition at line 412 of file markers.c.

References interface_splines(), interface::n, real, splint, interface::spx, interface::spy, interface::t, interface::x, and interface::y. 

00413 {
00414   int i;
00415   real t;
00416 
00417   for (i = 0; i < in->n - 1; i++) {
00418     t = (in->t[i] + in->t[i+1])/2.;
00419     in->x[i] = splint(in->spx[i], t);
00420     in->y[i] = splint(in->spy[i], t);
00421   }
00422   in->x[in->n-1] = in->x[0];
00423   in->y[in->n-1] = in->y[0];
00424   
00425   interface_splines(*in, nx, ny);
00426 }

void interface_redis (  interface *    in,
real    space
) 
 

Definition at line 327 of file markers.c.

References interface::n, real, splint, interface::spx, interface::spy, interface::t, interface::x, and interface::y.

Referenced by initialize(), and timestep(). 

00328 {
00329   int i, j, nb, k = 0;
00330   real l, dl, xe, ye;
00331 
00332   nb = (int)(in->t[in->n-1] / space) + 2;
00333   dl = in->t[in->n-1] / (real)(nb - 1);
00334   if (dl >= 1.0)
00335     printf("WARNING interface_redis(): dl=%g\n", dl);
00336 
00337   xe = in->x[in->n-1];  /* saving the coordinates of the last point */
00338   ye = in->y[in->n-1];  
00339 
00340   in->x = (real *)realloc(in->x, nb * sizeof(real));
00341   in->y = (real *)realloc(in->y, nb * sizeof(real));
00342 
00343   l = dl;               /* Starting at i=1 the first point is unchanged */
00344   for (i = 1; i < nb - 1; i++, l += dl) {
00345     /* find the next interval [j;j+1] such as t[j] <= l <= t[j+1] */
00346     while(in->t[k+1] < l)
00347       k++;
00348 
00349 #ifdef DEBUG_REDIS
00350     if (k >= in->n-1) {
00351       printf("interface_redis: Error i=%d/%d k too large: %d/%d\n",
00352               i, nb - 1, k, in->n - 1);
00353       exit(1);
00354     }
00355     if (in->t[k] > l || in->t[k+1] < l) {
00356       printf("interface_redis: Error i=%d/%d k=%d, %g <= %g <= %g\n",
00357               i, nb - 1, k, in->t[k], l, in->t[k+1]);
00358       exit(1);
00359     }
00360 #endif    
00361     
00362     in->x[i] = splint(in->spx[k], l);
00363     in->y[i] = splint(in->spy[k], l);
00364   }
00365 
00366   in->x[nb - 1] = xe;  /* The last point is unchanged */
00367   in->y[nb - 1] = ye;  /* Copying the saved coordinates */
00368   in->t = (real *)realloc(in->t, nb * sizeof(real));
00369   in->spx = (polynom3 *)realloc(in->spx, (nb - 1) * sizeof(polynom3));
00370   in->spy = (polynom3 *)realloc(in->spy, (nb - 1) * sizeof(polynom3));
00371 #ifdef FREE
00372   in->p = (real *)realloc(in->p, nb * sizeof(real));
00373 #endif
00374   in->n = nb;
00375 }

void interface_redis_scale (  interface *    in,
interface    src,
real    space,
real    scale
) 
 

Definition at line 378 of file markers.c.

References interface::n, real, splint, interface::spx, interface::spy, interface::t, interface::x, and interface::y.

Referenced by mgsolve(). 

00380 {
00381   int i, j, nb, k = 0;
00382   real l, dl, xe, ye;
00383 
00384   nb = (int)(src.t[src.n-1] / space) + 2;
00385   dl = src.t[src.n-1] / (real)(nb - 1);
00386 
00387   in->x = (real *)realloc(in->x, nb * sizeof(real));
00388   in->y = (real *)realloc(in->y, nb * sizeof(real));
00389   in->t = (real *)realloc(in->t, nb * sizeof(real));
00390   in->spx = (polynom3 *)realloc(in->spx, (nb - 1) * sizeof(polynom3));
00391   in->spy = (polynom3 *)realloc(in->spy, (nb - 1) * sizeof(polynom3));
00392 #ifdef FREE
00393   in->p = (real *)realloc(in->p, nb * sizeof(real));
00394 #endif
00395   in->n = nb;
00396 
00397   l = dl;               /* Starting at i=1 the first point is unchanged */
00398   for (i = 1; i < nb - 1; i++, l += dl) {
00399     /* find the next interval [j;j+1] such as t[j] <= l <= t[j+1] */
00400     while(src.t[k+1] < l)
00401       k++;
00402     in->x[i] = scale*(splint(src.spx[k], l) - 2.) + 2.;
00403     in->y[i] = scale*(splint(src.spy[k], l) - 2.) + 2.;
00404   }
00405   in->x[0] = scale*(src.x[0] - 2.) + 2.;
00406   in->y[0] = scale*(src.y[0] - 2.) + 2.;
00407   in->x[nb - 1] = scale*(src.x[src.n - 1] - 2.) + 2.;
00408   in->y[nb - 1] = scale*(src.y[src.n - 1] - 2.) + 2.;
00409 }

void interface_splines (  interface    in,
int    nx,
int    ny
) 
 

Definition at line 246 of file markers.c.

References APPROX_PERIODIC, AXIBOTH, AXINATURAL, AXITRUE, interface::bc, INTERPREC, interface::n, PERIODIC_DL, Pspline(), real, SP_DERIVATIVE, SP_NATURAL, spline(), interface::spx, interface::spy, sq, interface::t, TRUE_PERIODIC, interface::x, and interface::y.

Referenced by initialize(), interface_filter(), mgsolve(), and timestep(). 

00247 {
00248   real *d2y, l = 0.0, d0, d1, det, dy0, dy1, yp;
00249   real frac;
00250   int i;
00251 
00252   d2y = (real *)malloc(in.n * sizeof(real));
00253 
00254   in.t[0] = 0.0;
00255   for (i = 1; i < in.n; i++) {
00256     /* Move slightly the points when they are near the grid lines */
00257     frac = fabs(modf(in.x[i], &d0));
00258     if (frac < INTERPREC || 1. - frac < INTERPREC || 
00259         fabs(0.5 - frac) < INTERPREC) {
00260       printf("WARNING: interfaces_splines(): moving x[%d]: %.10f to %.10f\n",
00261              i, in.x[i], in.x[i] + 10.*INTERPREC);
00262       in.x[i] += 10.*INTERPREC;
00263     }
00264     frac = fabs(modf(in.y[i], &d0));
00265     if (frac < INTERPREC || 1. - frac < INTERPREC || 
00266         fabs(0.5 - frac) < INTERPREC) {
00267       printf("WARNING: interfaces_splines(): moving y[%d]: %.10f to %.10f\n", 
00268              i, in.y[i], in.y[i] + 10.*INTERPREC);
00269       in.y[i] += 10.*INTERPREC;
00270     }
00271 
00272     d0 = sqrt(sq(in.x[i] - in.x[i-1]) + sq(in.y[i] - in.y[i-1]));
00273 /*
00274     if (d0 >= 1.0)
00275       printf("WARNING: interfaces_splines(): d0=%f >= 1.0\n", d0);
00276 */
00277     l += d0;
00278     in.t[i] = l;
00279   }
00280 
00281   switch(in.bc) {
00282   case TRUE_PERIODIC:
00283     Pspline(in.t, in.x, d2y, in.t[in.n - 1], in.n - 1, in.spx);
00284     Pspline(in.t, in.y, d2y, in.t[in.n - 1], in.n - 1, in.spy);
00285     break;
00286   case APPROX_PERIODIC:
00287     /* We need an estimate of the derivative for the periodic ends 
00288        We suppose a 2nd order interpolation for the points 1, 0, n-2
00289        will do */
00290     d0 = in.t[in.n-1] - in.t[in.n-2];  
00291     d1 = d0 + in.t[1] - in.t[0];
00292     det = d1*d0*(d0 - d1);
00293     
00294     dy0 = in.x[0] - in.x[in.n-2];
00295     dy1 = in.x[1] - in.x[in.n-2];
00296     PERIODIC_DL(dy0, nx);
00297     PERIODIC_DL(dy1, nx);
00298     yp = (2.*(dy0*d1 - dy1*d0)*d0 + dy1*d0*d0 - dy0*d1*d1) / det;
00299     spline(in.t, in.x, in.n, yp, yp, d2y, in.spx, SP_DERIVATIVE);
00300     
00301     dy0 = in.y[0] - in.y[in.n-2];
00302     dy1 = in.y[1] - in.y[in.n-2];
00303     yp = (2.*(dy0*d1 - dy1*d0)*d0 + dy1*d0*d0 - dy0*d1*d1) / det;
00304     spline(in.t, in.y, in.n, yp, yp, d2y, in.spy, SP_DERIVATIVE);
00305     break;
00306   case AXITRUE:
00307     spline(in.t, in.x, in.n, 0.0, 0.0, d2y, in.spx, SP_DERIVATIVE);
00308     spline(in.t, in.y, in.n, 1.0, -1.0, d2y, in.spy, SP_DERIVATIVE);
00309     break;
00310   case AXINATURAL:
00311     spline(in.t, in.x, in.n, 0.0, 0.0, d2y, in.spx, SP_NATURAL);
00312     spline(in.t, in.y, in.n, 0.0, 0.0, d2y, in.spy, SP_NATURAL);
00313     break;
00314   case AXIBOTH:
00315     spline(in.t, in.x, in.n, 0.0, 0.0, d2y, in.spx, SP_DERIVATIVE);
00316     spline(in.t, in.y, in.n, 1.0, 1.0, d2y, in.spy, SP_DERIVATIVE);
00317     break;
00318   default:
00319     printf("interface_splines(): bad value for the bc parameter: bc = %d\n", in.bc);
00320     exit(1);
00321   }
00322     
00323   free(d2y);
00324 }

void interface_write_markers (  interface    in,
FILE *    fptr
) 
 

Definition at line 26 of file markers.c.

References interface::n, interface::x, and interface::y.

Referenced by initialize(), pressure(), and timestep(). 

00027 {
00028   int i;
00029 
00030   fprintf(fptr, "\n");
00031   for (i = 0; i <= in.n-1; i++)
00032     fprintf(fptr, "%g %g\n", in.x[i], in.y[i]); 
00033   fprintf(fptr, "&\n"); 
00034 }

void interface_write_markers_speed (  real    time,
interface    in,
real2D    u,
real2D    v,
real    tau,
real    du,
int    nx,
int    ny,
FILE *    fptr
) 
 

Definition at line 36 of file markers.c.

References interface::n, real, real2D, interface::x, and interface::y. 

00039 {
00040 /* Physical velocity*/
00041   int i;
00042   real h=1./(nx-2);
00043 
00044   fprintf(fptr, "\n");
00045   for (i = 0; i <= in.n-1; i++)
00046     fprintf(fptr, "%g %g %g %g %g\n", time, in.x[i], in.y[i], 
00047                     bilinu(u, in.x[i], in.y[i])*du*h/tau, 
00048                     bilinu(u, in.x[i], in.y[i])*du*h/tau); 
00049   fprintf(fptr, "&\n"); 
00050 }

void interface_write_splines (  interface    in,
FILE *    fptr,
int    nb
) 
 

Definition at line 52 of file markers.c.

References interface::n, real, splint_find(), interface::spx, interface::spy, and interface::t.

Referenced by closest_normal(), initialize(), and timestep(). 

00053 {
00054   int i;
00055   real t;
00056 
00057   fprintf(fptr, "\n");
00058   for (i = 0; i < nb; i++) {
00059     t = (real) i * in.t[in.n-1] / (real) (nb - 1);
00060     fprintf(fptr, "%g %g\n", 
00061             splint_find(in.t, in.spx, in.n, t),
00062             splint_find(in.t, in.spy, in.n, t));
00063   }
00064   fprintf(fptr, "&\n");
00065 }

real interfaces_area (  interface *    in,
int    n
) 
 

Definition at line 85 of file markers.c.

References polynom3::a, polynom3::b, polynom3::c, polynom3::d, interface::n, real, interface::spx, interface::spy, and interface::t. 

00086 {
00087   int i, j;
00088   real ar = 0.0, a1, a2, a3, a4, a5, a6, t2, t1;
00089   interface in1;
00090 
00091   for (j = 0; j < n; j++) {
00092     in1 = in[j];
00093     for (i = 0; i < in1.n - 1; i++) 
00094       {
00095         a1 = in1.spx[i].a*in1.spy[i].b; 
00096         a2 = 0.5*(in1.spx[i].b*in1.spy[i].b + 2.*in1.spx[i].a*in1.spy[i].c);
00097         a3 = (in1.spx[i].c*in1.spy[i].b + 2.*in1.spx[i].b*in1.spy[i].c 
00098               + 3.*in1.spx[i].a*in1.spy[i].d)/3.;
00099         a4 = 0.25*(in1.spx[i].d*in1.spy[i].b + 2.*in1.spx[i].c*in1.spy[i].c 
00100                    + 3.*in1.spx[i].b*in1.spy[i].d);
00101         a5 = (3.*in1.spx[i].c*in1.spy[i].d + 2.*in1.spx[i].d*in1.spy[i].c)/5.;
00102         a6 = 0.5*in1.spx[i].d*in1.spy[i].d;
00103         t1 = in1.t[i]; t2 = in1.t[i+1];
00104         ar += t2*(a1 + t2*(a2 + t2*(a3 + t2*(a4 + t2*(a5 + t2*a6))))) -
00105           t1*(a1 + t1*(a2 + t1*(a3 + t1*(a4 + t1*(a5 + t1*a6)))));
00106       }
00107   }
00108   return ar;
00109 }

real interfaces_axivolume (  interface *    in,
int    n
) 
 

Definition at line 112 of file markers.c.

References polynom3::a, polynom3::b, polynom3::c, polynom3::d, interface::n, real, interface::spx, interface::spy, interface::t, interface::x, and interface::y.

Referenced by timestep(). 

00113 {
00114   int i, j;
00115   real ar = 0.0, a1, a2, a3, a4, a5, a6, a7, a8, a9, t2, t1;
00116   real ax, bx, cx, dx, ay, by, cy, dy;
00117   real x1, y1, x2, y2;
00118   interface in1;
00119 
00120   for (j = 0; j < n; j++) {
00121     in1 = in[j];
00122     for (i = 0; i < in1.n - 1; i++) {
00123 #ifdef POLYVOLUME
00124       /* integration of xy dy along the polynomial segment */
00125       ax = in1.spx[i].a; bx = in1.spx[i].b;
00126       cx = in1.spx[i].c; dx = in1.spx[i].d;
00127       ay = in1.spy[i].a - 2.0; by = in1.spy[i].b;
00128       cy = in1.spy[i].c; dy = in1.spy[i].d;
00129       
00130       a1 = ax*ay*by;
00131       a2 = ax*ay*cy+(ax*by+bx*ay)*by/2.;
00132       a3 = ax*ay*dy+2.0/3.0*(ax*by+bx*ay)*cy+(ax*cy+bx*by+cx*ay)*by/3.;
00133       a4 = 3.0/4.0*(ax*by+bx*ay)*dy+(ax*cy+bx*by+cx*ay)*cy/2.
00134         +(ax*dy+bx*cy+cx*by+dx*ay)*by/4.;
00135       a5 = 3.0/5.0*(ax*cy+bx*by+cx*ay)*dy
00136         +2.0/5.0*(ax*dy+bx*cy+cx*by+dx*ay)*cy+(bx*dy+cx*cy+dx*by)*by/5.;
00137       a6 = (ax*dy+bx*cy+cx*by+dx*ay)*dy/2+(bx*dy+cx*cy+dx*by)*cy/3
00138         +(cx*dy+dx*cy)*by/6.;
00139       a7 = 3.0/7.0*(bx*dy+cx*cy+dx*by)*dy+2.0/7.0*(cx*dy+dx*cy)*cy
00140         +dx*dy*by/7.;
00141       a8 = 3.0/8.0*(cx*dy+dx*cy)*dy+dx*dy*cy/4.;
00142       a9 = dx*dy*dy/3.;
00143       
00144       t1 = in1.t[i]; t2 = in1.t[i+1];
00145       ar += t2*(a1 + t2*(a2 + t2*(a3 + t2*(a4 + t2*(a5 + t2*(a6 + t2*(a7 + t2*(a8 + t2*a9)))))))) -
00146         t1*(a1 + t1*(a2 + t1*(a3 + t1*(a4 + t1*(a5 + t1*(a6 + t1*(a7 + t1*(a8 + t1*a9))))))));
00147 #else
00148       /* integration of xy dy along the segment */
00149       x1 = in1.x[i]; y1 = in1.y[i] - 2.;
00150       x2 = in1.x[i+1]; y2 = in1.y[i+1] - 2.;
00151       ar += ((x2 - x1)*(y2 - y1)/3. + 
00152              (y1*(x2 - x1) + x1*(y2 - y1))/2. + x1*y1)*(y2 - y1);
00153 #endif
00154     }
00155   }
00156   return 2.*PI*ar;
00157 }

real interfaces_x_amp (  interface *    in,
int    n
) 
 

Definition at line 218 of file markers.c.

References interface::n, real, and interface::x. 

00219 {
00220   int i, j;
00221   real max = -1e6, min = 1e6;
00222 
00223   for (j = 0; j < n; j++)
00224     for (i = 0; i < in[j].n; i++) {
00225       max = in[j].x[i] > max ? in[j].x[i] : max;
00226       min = in[j].x[i] < min ? in[j].x[i] : min;
00227     }
00228   return fabs(max - min);
00229 }

real interfaces_xmoment (  interface *    in,
int    n
) 
 

Definition at line 160 of file markers.c.

References interface::n, real, splint, splint1, interface::spx, interface::spy, and interface::t. 

00161 {
00162   /* compute Sum{x^2.dx.dy} */
00163   /* coefficients of Newton-Cotes formula */
00164   static real coef[7] = {1., 3., 3., 2., 3., 3., 1.}; 
00165   int i, j, k;
00166   real x, yp, h, t1, t2, t, sum, total = 0.0;
00167   interface in1;
00168   polynom3 px, py;
00169 
00170   for (k = 0; k < n; k++) {
00171     in1 = in[k];
00172     for (j = 0; j < in1.n - 1; j++) {
00173       t1 = in1.t[j]; t2 = in1.t[j+1];
00174       px = in1.spx[j]; py = in1.spy[j];
00175       sum = 0.0; h = (t2 - t1)/6.;
00176       for (i = 0; i < 7; i++) {
00177         t = t1 + h*i;
00178         x = splint(px, t) - 2.;
00179         yp = splint1(py, t);
00180         sum += coef[i]*x*x*x*yp;
00181       }
00182       total += 3.*h*sum/8.;
00183     }
00184   }
00185   return total/3.;
00186 }

real interfaces_y_amp (  interface *    in,
int    n
) 
 

Definition at line 232 of file markers.c.

References interface::n, real, and interface::y. 

00233 {
00234   int i, j;
00235   real max = -1e6, min = 1e6;
00236 
00237   for (j = 0; j < n; j++)
00238     for (i = 0; i < in[j].n; i++) {
00239       max = in[j].y[i] > max ? in[j].y[i] : max;
00240       min = in[j].y[i] < min ? in[j].y[i] : min;
00241     }
00242   return fabs(max - min);
00243 }

real interfaces_ymoment (  interface *    in,
int    n
) 
 

Definition at line 189 of file markers.c.

References interface::n, real, splint, splint1, interface::spx, interface::spy, and interface::t. 

00190 {
00191   /* compute Sum{y^2.dx.dy} */
00192   /* coefficients of Newton-Cotes formula */
00193   static real coef[7] = {1., 3., 3., 2., 3., 3., 1.}; 
00194   int i, j, k;
00195   real y, xp, h, t1, t2, t, sum, total = 0.0;
00196   interface in1;
00197   polynom3 px, py;
00198 
00199   for (k = 0; k < n; k++) {
00200     in1 = in[k];
00201     for (j = 0; j < in1.n - 1; j++) {
00202       t1 = in1.t[j]; t2 = in1.t[j+1];
00203       px = in1.spx[j]; py = in1.spy[j];
00204       sum = 0.0; h = (t2 - t1)/6.;
00205       for (i = 0; i < 7; i++) {
00206         t = t1 + h*i;
00207         y = splint(py, t) - 2.;
00208         xp = splint1(px, t);
00209         sum += coef[i]*y*y*y*xp;
00210       }
00211       total += 3.*h*sum/8.;
00212     }
00213   }
00214   return total/3.;
00215 }

Generated on Wed Feb 19 22:27:49 2003 for Markers by doxygen1.2.18