Float u_init (const point& x) {
    if (x[0] < 0.6) return 1;
    if (x[0] < 0.8) return 5*(0.8-x[0]);
    return 0;
}
point f (const Float& u) { return point(sqr(u)/2); }
Float g = 1;
Float phi (const point& nu, Float a, Float b) {
  if ((nu[0] >= 0 && a <= b) || (nu[0] <= 0 && a >= b)) 
    return nu[0]*min(sqr(a),sqr(b))/2;
  else
    return nu[0]*max(sqr(a),sqr(b))/2;
}
