/**
#Surface tension effects using the compressible solver

*/

#include "tension.h"

event end_timestep (i++)
{
  face vector pf[];

  scalar skappa[];
  foreach ()
    skappa[] = 0.;
  if (f.sigma > 0.)
    curvature (f, skappa, f.sigma, add = true);
  boundary({skappa});

  foreach_face () {
    pf.x[] = 0.;
    if ((f[1] + f[])/2. > 0. && (f[1] + f[])/2. < 1. ) {
      if (f[1] > f[])
        pf.x[] = (1. - f[1])*skappa[1];
      else
        pf.x[] = (1. - f[])*skappa[];
    }
  }
  boundary((scalar *){pf});

  foreach () {
    double div1 = 0.;
    foreach_dimension () 
      div1 += (pf.x[1]*uf.x[1] - pf.x[]*uf.x[]);

    fE1[] -= f[]*div1/Delta*dt/cm[];
  }

  foreach_face () {
    pf.x[] = 0.;
    if ((f[1] + f[])/2. > 0. && (f[1] + f[])/2. < 1. ) {
      if (f[1] > f[])
        pf.x[] = - f[]*skappa[];
      else
        pf.x[] = - f[1] *skappa[1];
    }
  }
  boundary((scalar *){pf});

  foreach () {
    double div1 = 0.;
    foreach_dimension () 
      div1 += (pf.x[1]*uf.x[1] - pf.x[]*uf.x[]);

    fE2[] -= (1. - f[])*div1/Delta*dt/cm[];
  }

  boundary({fE1,fE2});
}