sandbox/popinet/meanflow.c
Evolution by mean curvature flow
Solve \displaystyle \mathbf{n}\cdot\partial_\tau\mathbf{X} = -\kappa See e.g Deckelnick et al, 2005.
set size ratio -1
plot 'out' w l t ''
#include "grid/multigrid.h"
#include "advection.h"
#include "vof.h"
#include "curvature.h"
scalar f[];
scalar * tracers = NULL, * interfaces = {f};
int main()
{
L0 = 1 [0];
origin (-0.5, -0.5);
N = 64;
// this is the diffusive timestep limit (with diffusion coeff unity)
DT = sq(L0/N)/2.;
run();
}
double circle (double x, double y)
{
double theta = atan2(y,x);
double r = sqrt(sq(x) + sq(y));
return r - 0.4*(1. + 0.2*cos(8.*theta));
}
event init (i = 0)
{
fraction (f, circle(x, y));
}
coord normal (Point point, scalar c)
{
coord n = mycs (point, c);
double nn = 0.;
foreach_dimension()
nn += sq(n.x);
nn = sqrt(nn);
foreach_dimension()
n.x /= nn;
return n;
}
scalar kappa[];
event velocity (i++) {
curvature (f, kappa);
foreach_face() {
coord n;
foreach_dimension()
n.x = 0.;
double k = 0.;
if (f[] > 0. && f[] < 1.) {
n = normal (point, f);
k = kappa[];
if (f[-1] > 0. && f[-1] < 1.) {
coord np = normal (neighborp(-1), f);
foreach_dimension()
n.x = (n.x + np.x)/2.;
k = (k + kappa[-1])/2.;
}
}
else if (f[-1] > 0. && f[-1] < 1.) {
k = kappa[-1];
n = normal (neighborp(-1), f);
}
u.x[] = - (k /*+ 2.5*/)*n.x;
}
#if 0
FILE * fp = fopen("velo", "w");
foreach()
fprintf (fp, "%g %g %g %g\n", x, y,
(uf.x[] + uf.x[1])/2.,
(uf.y[] + uf.y[0,1])/2.);
fclose (fp);
#endif
}
event logfile (t += 0.001; t <= sq(0.1)) {
// stats s = statsf (kappa);
// fprintf (stderr, "# %f %d %.12f %g %g\n", t, i, s.sum, s.min, s.max);
output_facets (f, stdout);
}
#if 0
event gfsview (i++) {
static FILE * fp = popen ("gfsview2D -s mean.gfv", "w");
output_gfs (fp, t = t);
}
#endif