/**
# 3D meniscus

An initial perturbed free-surface is contained within a cube. The contact
angles on the faces of the cube are set to 45 degrees. A constant
surface tension is imposed and the viscous fluid relaxes toward its
equilibrium position: a 3D meniscus i.e. the intersection of a sphere
and a cube with the given volume and contact angles.

![Animation of the free-surface relaxation](meniscus3D/movie.mp4)

Time-implicit timestepping is used and greatly decreases the number of
timesteps required to reach convergence.

Both the maximum velocity and standard deviation on the curvature
decrease exponentially.

~~~gnuplot Convergence of the maximum velocity and standard deviation on curvature
set xlabel 'Time'
set logscale y
set xrange [0:5]
plot 'log' u 1:3 w l t 'Maximum velocity', \
     'log' u 1:5 w l t 'Standard deviation of curvature'
~~~
*/

#include "grid/multigrid.h"
#include "crobert/2_Implicit/hydro-tension.h"
#include "layered/implicit.h"
#include "layered/remap.h"
#include "view.h"

int main()
{
  L0 = 2.;
  origin (-L0/2., -L0/2.);

  G = 0.;
  nu = 1.;

  /**
  We use a large timestep and decrease the tolerance to get clear
  convergence of the maximum velocity and standard deviation of
  curvature. */
  
  TOLERANCE = 1e-8;
  CFL_H = 40;
  
  N = 64;
  nl = 1;

  run();
}

/**
## Contact angle boundary conditions 

This is necessary to balance the surface tension acceleration on the
side walls. This needs to match the contact angle boundary conditions
below. */

event half_advection (i++)
{
  foreach_dimension() {
    ha.n[left] = 0.;
    ha.n[right] = 0.;
    hu.n[left] = 0.;
    hu.n[right] = 0.;
    hf.n[left] = 0.;
    hf.n[right] = 0.;
  }
  boundary ((scalar *){ha, hu, hf});  
}

event init (i = 0)
{
  foreach_dimension() {
    eta[left]    = neumann (1./tan(45.*pi/180.));
    eta[right]   = neumann (1./tan(45.*pi/180.));
  }

  /**
  ## Initial conditions

  The initial free-surface shape. */

  foreach() {
    double H = 0.5 + 0.1*cos(2.*pi*x)*cos(2.*pi*y);
    zb[] = - 0.5;
    foreach_layer()
      h[] = H/nl;
  }
}

/**
## Outputs

We compute the interface curvature (times $\sigma$, which is
unity). This needs to be done in the `pressure()` event since the
$\sigma$ fields used by the `sigma_kappa()` macro are temporary. */

scalar kappa[];

event pressure (i++)
{
  foreach()
    kappa[] = sigma_kappa (eta, 0);
}
  
event logfile (i++)
{
  stats s = statsf (kappa);
  fprintf (stderr, "%g %g %g %g %g %d\n", t, dt, normf(u.x).max,
	   s.sum/s.volume, s.stddev, mgH.i);
}

event movie (i++; t <= 5)
{
  view (fov = 22, quat = {0.475152,0.161235,0.235565,0.832313},
	tx = 0.02, ty = 0., width = 800, height = 600);
  char s[80];
  sprintf (s, "t = %.2f", t);
  draw_string (s, size = 80);
  squares ("eta", linear = true, z = "eta", min = -0.3, max = 0.3);
  save ("movie.mp4");
}