///
/// This file is part of Rheolef.
///
/// Copyright (C) 2000-2009 Pierre Saramito <Pierre.Saramito@imag.fr>
///
/// Rheolef is free software; you can redistribute it and/or modify
/// it under the terms of the GNU General Public License as published by
/// the Free Software Foundation; either version 2 of the License, or
/// (at your option) any later version.
///
/// Rheolef is distributed in the hope that it will be useful,
/// but WITHOUT ANY WARRANTY; without even the implied warranty of
/// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
/// GNU General Public License for more details.
///
/// You should have received a copy of the GNU General Public License
/// along with Rheolef; if not, write to the Free Software
/// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307  USA
/// 
/// =========================================================================
/// Banded level set routines 
///
/// Authors: Lara Aborm, Jocelyn Etienne, Pierre Saramito
///
#include "rheolef/level_set.h"
#include "rheolef/field.h"
namespace rheolef {

// --------------------------------------------------------------------------------
// gestion de la numerotation locale
// --------------------------------------------------------------------------------
// TODO: move in reference_element
static
size_t
edge_t_iloc (size_t l, size_t m)
{
  static const int edge_t_iloc_table [3][3] = { 
			    {-1, 0, 2},
			    { 0,-1, 1},
			    { 2, 1,-1}};
  return size_t(edge_t_iloc_table [l][m]);
}
static
size_t
edge_T_iloc (size_t l, size_t m)
{
  static const int edge_T_iloc_table [4][4] = {
                            {-1, 0, 2, 3},
                            { 0,-1, 1, 4},
                            { 2, 1,-1, 5},
                            { 3, 4, 5,-1}};
  return size_t(edge_T_iloc_table [l][m]);
}
static
size_t
face_T_iloc (size_t l, size_t m, size_t n)
{
  static size_t face_T_iloc_table [4][4][4];
  bool static initialized = false;
  if (!initialized) {
    for (size_t i = 0; i < 4; i++) 
      for (size_t j = 0; j < 4; j++) 
        for (size_t k = 0; k < 4; k++) 
          face_T_iloc_table [i][j][k] = size_t(-1);
    reference_element hat_K (reference_element::T);
    for (size_t i_face = 0; i_face < hat_K.n_face(); i_face++) {
      size_t p[3];
      for (size_t k = 0; k < 3; k++) p[k] = hat_K.subgeo_local_vertex(2,i_face,k);
      face_T_iloc_table [p[0]][p[1]][p[2]] = i_face;
      face_T_iloc_table [p[0]][p[2]][p[1]] = i_face;
      face_T_iloc_table [p[1]][p[0]][p[2]] = i_face;
      face_T_iloc_table [p[1]][p[2]][p[0]] = i_face;
      face_T_iloc_table [p[2]][p[0]][p[1]] = i_face;
      face_T_iloc_table [p[2]][p[1]][p[0]] = i_face;
    }
  }
  return face_T_iloc_table [l][m][n];
}
// --------------------------------------------------------------------------------
// gestion de la precision
// --------------------------------------------------------------------------------
static Float level_set_epsilon = 100*std::numeric_limits<Float>::epsilon();

template <class T>
static
bool
is_zero (const T& x) {
  // TODO: control by level_set_option ?
  return fabs(x) <= level_set_epsilon;
}
template <class T>
static
bool
have_same_sign (const T& x, const T& y) {
  using namespace details;
  return !is_zero(x) && !is_zero(y) && x*y > 0;
}
template <class T>
static
bool
have_opposite_sign (const T& x, const T& y) {
  using namespace details;
  return !is_zero(x) && !is_zero(y) && x*y < 0;
}
// --------------------------------------------------------------------------------
// 2D: fonctions locales : sur un seul element triangle
// --------------------------------------------------------------------------------
// appele lors du 1er passage qui liste les elements de la bande cas de dimension 2
template <class T>
static
bool
belongs_to_band_t (const std::vector<T>& f) {
  using namespace details;
  if (have_same_sign(f[0],f[1]) && have_same_sign(f[0],f[2])) return false;
  // on rejette le triangle dans tous les sommets de meme signe :
  if (is_zero(f[0]) && is_zero(f[1]) && is_zero(f[2])) return false;
  // on rejette les triangles dont un sommet :
  if (is_zero(f[0]) && have_same_sign(f[1],f[2])) return false;
  if (is_zero(f[1]) && have_same_sign(f[0],f[2])) return false;
  if (is_zero(f[2]) && have_same_sign(f[0],f[1])) return false;
  return true;
}
// apellee lors du calcul des matrices M_K et A_K pour K dans la bande cas de dimension 2
template <class T>
static
size_t
isolated_vertex_t (const std::vector<T>& f)
{
  using namespace details;
  /* on retourne le sommet isole a chaque fois */
  if (have_same_sign    (f[0],f[1]) && have_opposite_sign(f[0],f[2]))  return 2;
  if (have_opposite_sign(f[0],f[1]) && have_same_sign    (f[0],f[2]))  return 1;
  if (have_opposite_sign(f[0],f[1]) && have_opposite_sign(f[0],f[2]))  return 0;
  if (is_zero(f[0])                 && have_opposite_sign(f[1],f[2]))  return 1; /* on peut retourner 2 de meme*/

  if (is_zero(f[1]) && have_opposite_sign(f[0],f[2]))  return 0; /* on peut retourner 2 */
  if (is_zero(f[2]) && have_opposite_sign(f[0],f[1]))  return 0; /* on peut retourner 1 */
  if (is_zero(f[0]) && is_zero(f[1]) && ! is_zero(f[2])) return 2;
  if (is_zero(f[0]) && is_zero(f[2]) && ! is_zero(f[1])) return 1;
  return 0; /* f1 == 0 et f2 == 0 et f0 != 0 */
}
template <class T>
static
void
subcompute_matrix_t (
    const std::vector<point_basic<T> >& x,
    const std::vector<T>&               f,
    std::vector<size_t>&                j,
    point_basic<T>&                     a,
    point_basic<T>&                     b,
    T&                                  S)
{
  using namespace details;
  j.resize (3);
  j[0] = isolated_vertex_t (f);
  j[1] = (j[0]+1) % 3;
  j[2] = (j[1]+1) % 3;
  // edge {j1,j2} has normal oriented as grad(f), in f>0 direction:
  if (! is_zero(f[j[0]]) && f[j[0]] < 0) std::swap (j[1], j[2]);
  T theta_1= f[j[1]]/(f[j[1]]-f[j[0]]);
  T theta_2= f[j[2]]/(f[j[2]]-f[j[0]]);
  // calcul des coordonnes d'intersection
  a = theta_1*x[j[0]]+(1-theta_1)*x[j[1]];
  b = theta_2*x[j[0]]+(1-theta_2)*x[j[2]];
  S = sqrt(pow(a[0]-b[0],2)+pow(a[1]-b[1],2));
  if (is_zero(f[j[1]]) && is_zero(f[j[2]])) {
      S /= 2;
  }
}
// --------------------------------------------------------------------------------
// 3D: fonctions locales : sur un seul element tetraedre
// --------------------------------------------------------------------------------

class quadruplet {
public:
    quadruplet (size_t a=0, size_t b=0, size_t c=0, size_t d=0) {
        q[0]=a;
        q[1]=b;
        q[2]=c;
        q[3]=d;
    }
    size_t  operator[] (size_t i) const {
        return q[i%4];
    }
    size_t& operator[] (size_t i)       {
        return q[i%4];
    }
    friend std::ostream& operator<< (std::ostream& out, const quadruplet& q) {
        out << "((" << q[0] << "," << q[1] << "), (" << q[2] << "," << q[3] << "))";
        return out;
    }
protected:
    size_t q[4];
};
// appele lors du 1er passage qui liste les elements de la bande cas de dimension 3
template <class T>
static
bool
belongs_to_band_T (const std::vector<T>& f)
{
  using namespace details;
  if (have_same_sign(f[0],f[1]) && have_same_sign(f[0],f[2]) && have_same_sign(f[2],f[3])) return false;
  // cas ou 4 points s'annulent en dimension 3 est degenere
  if (is_zero(f[0]) && is_zero(f[1]) && is_zero(f[2]) && is_zero(f[3])) return false;

  if (is_zero(f[0]) && have_same_sign(f[1],f[2]) && have_same_sign(f[1],f[3])) return false;
  if (is_zero(f[1]) && have_same_sign(f[0],f[2]) && have_same_sign(f[0],f[3])) return false;
  if (is_zero(f[2]) && have_same_sign(f[0],f[1]) && have_same_sign(f[0],f[3])) return false;
  if (is_zero(f[3]) && have_same_sign(f[0],f[1]) && have_same_sign(f[0],f[2])) return false;
  // cas ou f s'annule sur 2 sommets et garde le meme signe sur les 2 autres sommets est exclu
  if (is_zero(f[0]) && is_zero(f[1]) && have_same_sign(f[2],f[3])) return false;
  if (is_zero(f[0]) && is_zero(f[2]) && have_same_sign(f[1],f[3])) return false;
  if (is_zero(f[0]) && is_zero(f[3]) && have_same_sign(f[1],f[2])) return false;
  if (is_zero(f[1]) && is_zero(f[2]) && have_same_sign(f[0],f[3])) return false;
  if (is_zero(f[1]) && is_zero(f[3]) && have_same_sign(f[0],f[2])) return false;
  if (is_zero(f[2]) && is_zero(f[3]) && have_same_sign(f[1],f[0])) return false;
  return true;
}
// apellee lors du calcul des matrices M_K et A_K pour T dans la bande cas de dimension 
template <class T>
bool
intersection_is_quadrilateral_T (const std::vector<T>& f, quadruplet& q)
{
  if (have_same_sign(f[0],f[1]) && have_opposite_sign(f[0],f[2]) && have_same_sign(f[2],f[3])) {
    if (f[0] > 0) q = quadruplet(0,1, 2,3);
    else          q = quadruplet(2,3, 0,1);
    return true;
  }
  if (have_opposite_sign(f[0],f[1]) && have_same_sign(f[0],f[2]) && have_opposite_sign(f[2],f[3])) {
    if (f[0] < 0) q = quadruplet(0,2, 1,3);
    else          q = quadruplet(1,3, 0,2);
    return true;
  }
  if (have_opposite_sign(f[0],f[1]) && have_opposite_sign(f[0],f[2]) && have_opposite_sign(f[2],f[3])) {
    if (f[0] > 0) q = quadruplet(0,3, 1,2);
    else          q = quadruplet(1,2, 0,3);
    return true;
  }
  return false;
}
// cas d'une intersection triangle:
template <class T>
static
size_t
isolated_vertex_T (const std::vector<T>& f)
{
  using namespace details;
  // cas ou l'intersection est un triangle
  if (have_opposite_sign(f[0],f[1]) && have_opposite_sign(f[0],f[2]) && have_same_sign    (f[2],f[3])) return 0;
  if (have_same_sign    (f[0],f[1]) && have_opposite_sign(f[0],f[2]) && have_opposite_sign(f[2],f[3])) return 2;
  if (have_same_sign    (f[0],f[1]) && have_same_sign    (f[0],f[2]) && have_opposite_sign(f[2],f[3])) return 3;
  // cas ou f s'annule sur un sommet et change de signe sur les 2 autres sommets
  if (have_opposite_sign(f[0],f[1]) && have_same_sign(f[0],f[2]) && have_same_sign(f[2],f[3])) return 1;
  if (is_zero(f[0]) && have_same_sign    (f[1],f[2]) && have_opposite_sign(f[1],f[3])) return 3;
  if (is_zero(f[0]) && have_opposite_sign(f[1],f[2]) && have_same_sign    (f[1],f[3])) return 2;
  if (is_zero(f[0]) && have_opposite_sign(f[1],f[2]) && have_opposite_sign(f[1],f[3])) return 1;
  if (is_zero(f[1]) && have_opposite_sign(f[0],f[2]) && have_same_sign    (f[0],f[3])) return 2;
  if (is_zero(f[1]) && have_same_sign    (f[0],f[2]) && have_opposite_sign(f[0],f[3])) return 3;
  if (is_zero(f[1]) && have_opposite_sign(f[0],f[2]) && have_opposite_sign(f[0],f[3])) return 0;
  if (is_zero(f[2]) && have_opposite_sign(f[0],f[1]) && have_same_sign    (f[0],f[3])) return 1;
  if (is_zero(f[2]) && have_same_sign    (f[0],f[1]) && have_opposite_sign(f[0],f[3])) return 3;
  if (is_zero(f[2]) && have_opposite_sign(f[0],f[1]) && have_opposite_sign(f[0],f[3])) return 0;
  if (is_zero(f[3]) && have_opposite_sign(f[0],f[1]) && have_same_sign    (f[0],f[2])) return 1;
  if (is_zero(f[3]) && have_same_sign    (f[0],f[1]) && have_opposite_sign(f[0],f[2])) return 2;
  if (is_zero(f[3]) && have_opposite_sign(f[0],f[1]) && have_opposite_sign(f[0],f[2])) return 0;
  // cas ou f s'annule en 2 sommets et change de signe sur les 2 autres sommets
  if (is_zero(f[0]) && is_zero(f[1]) && have_opposite_sign(f[2],f[3])) return 2; // ou 3
  if (is_zero(f[0]) && is_zero(f[2]) && have_opposite_sign(f[1],f[3])) return 1;
  if (is_zero(f[0]) && is_zero(f[3]) && have_opposite_sign(f[1],f[2])) return 1;
  if (is_zero(f[1]) && is_zero(f[2]) && have_opposite_sign(f[0],f[3])) return 0;
  if (is_zero(f[1]) && is_zero(f[3]) && have_opposite_sign(f[0],f[2])) return 0;
  if (is_zero(f[2]) && is_zero(f[3]) && have_opposite_sign(f[0],f[1])) return 0;

  // le triangle d'intersection est la face du tetradre ou f s'annule sur les 3 sommets
  if (is_zero(f[0]) && is_zero(f[1]) && is_zero(f[2]) && ! is_zero(f[3])) return 3;
  if (is_zero(f[0]) && is_zero(f[1]) && is_zero(f[3]) && ! is_zero(f[2])) return 2;
  if (is_zero(f[1]) && is_zero(f[2]) && is_zero(f[3]) && ! is_zero(f[0])) return 0;
  return 1;
}
template <class T>
static
void
subcompute_matrix_triangle_T (
  const std::vector<point_basic<T> >& x,
  const std::vector<T>&               f,
  std::vector<size_t>&                j,
  point_basic<T>&                     a,
  point_basic<T>&                     b,
  point_basic<T>&                     c,
  T&                                  aire)
{
  using namespace details;
  j.resize(4);
  j[0] = isolated_vertex_T (f);
  j[1] = (j[0]+1) % 4;
  j[2] = (j[1]+1) % 4;
  j[3] = (j[2]+1) % 4;
  // orient
  if (! is_zero(f[j[0]]) && ((f[j[0]] > 0 && j[0] % 2 == 0) || (f[j[0]] < 0 && j[0] % 2 == 1))) 
    std::swap (j[1], j[2]);
  T theta_1 = f[j[1]]/(f[j[1]]-f[j[0]]);
  T theta_2 = f[j[2]]/(f[j[2]]-f[j[0]]);
  T theta_3 = f[j[3]]/(f[j[3]]-f[j[0]]);

  /* calcul des coordonnees d'intersection */
  a = theta_1*x[j[0]]+(1-theta_1)*x[j[1]];
  b = theta_2*x[j[0]]+(1-theta_2)*x[j[2]];
  c = theta_3*x[j[0]]+(1-theta_3)*x[j[3]];
  aire = 0.5* norm (vect( b-a, c-a));
  if (is_zero(f[j[1]]) && is_zero(f[j[2]]) && is_zero(f[j[3]])) {
    aire /= 2;
  }
}
template <class T>
static
void
subcompute_matrix_quadrilateral_T (
    const std::vector<point_basic<T> >& x,
    const std::vector<T>&               f,
    const quadruplet&                   q,
    point_basic<T>&                     a,
    point_basic<T>&                     b,
    point_basic<T>&                     c,
    point_basic<T>&                     d,
    T&                                  aire_Q)
{
    // intersection:
    // a = segment {x(q0) x(q2)} inter {f=0}
    // b = segment {x(q1) x(q2)} inter {f=0}
    // d = segment {x(q1) x(q3)} inter {f=0}
    // c = segment {x(q0) x(q3)} inter {f=0}
    // quadrilatere abdc = triangle(abd) union triangle(adc)
    T theta_1 = f[q[2]]/(f[q[2]]-f[q[0]]);
    T theta_2 = f[q[2]]/(f[q[2]]-f[q[1]]);
    T theta_3 = f[q[3]]/(f[q[3]]-f[q[0]]);
    T theta_4 = f[q[3]]/(f[q[3]]-f[q[1]]);
    /* calcul des coordonnees d'intersection */
    a = theta_1*x[q[0]]+(1-theta_1)*x[q[2]];
    b = theta_2*x[q[1]]+(1-theta_2)*x[q[2]];
    c = theta_3*x[q[0]]+(1-theta_3)*x[q[3]];
    d = theta_4*x[q[1]]+(1-theta_4)*x[q[3]];
    aire_Q = 0.5*norm(vect(a-c,a-b)) + 0.5*norm(vect(d-c,d-b));
}
// --------------------------------------------------------------------------------
// level_set: compte the intersection mesh
// --------------------------------------------------------------------------------
typedef geo_element_auto<heap_allocator<geo_element::size_type> > element_type;

template <class T, class M>
void
gamma_list2disarray (
  const std::list<point_basic<T> >&             gamma_node_list,
  std::array<std::list<std::pair<element_type,size_t> >, 
               reference_element::max_variant>  gamma_side_list,
  const communicator&                           comm,
  size_t                                        d,
  disarray<point_basic<T>, M>&                  gamma_node,
  std::array<disarray<element_type,M>,
	       reference_element::max_variant>& gamma_side,
  disarray<size_t,M>&                           sid_ie2bnd_ie)
{
  typedef geo_element::size_type size_type;
  // 1) nodes:
  size_type nnod = gamma_node_list.size();
  distributor gamma_node_ownership (distributor::decide, comm, nnod);
  gamma_node.resize (gamma_node_ownership);
  typename disarray<point_basic<T>, M>::iterator node_iter = gamma_node.begin();
  for (typename std::list<point_basic<T> >::const_iterator
	iter = gamma_node_list.begin(),
	last = gamma_node_list.end();
	iter != last; iter++, node_iter++) {
    *node_iter = *iter;
  }
  // 2) sides:
  heap_allocator<size_type> alloc;
  size_type map_dim = d-1;
  size_type order = 1;
  size_type nsid = 0;
  for (size_type variant = reference_element::first_variant_by_dimension(map_dim);
                 variant < reference_element:: last_variant_by_dimension(map_dim); variant++) {
    size_type nsidv = gamma_side_list [variant].size();
    distributor gamma_sidv_ownership (distributor::decide, comm, nsidv);
    nsid += nsidv;
    element_type element_init (variant, order, alloc);
    gamma_side[variant].resize (gamma_sidv_ownership, element_init);
    typename disarray<element_type, M>::iterator side_iter = gamma_side[variant].begin();
    for (typename std::list<std::pair<element_type,size_type> >::const_iterator
	  iter = gamma_side_list[variant].begin(),
	  last = gamma_side_list[variant].end();
	  iter != last; iter++, side_iter++) {
      *side_iter = (*iter).first;
    }
  }
  // 3) side2band correspondance
  distributor gamma_sid_ownership (distributor::decide, comm, nsid);
  const size_type undef      = std::numeric_limits<size_t>::max();
  sid_ie2bnd_ie.resize (gamma_sid_ownership, undef);
  typename disarray<size_type, M>::iterator idx_iter = sid_ie2bnd_ie.begin();
  for (size_type variant = reference_element::first_variant_by_dimension(map_dim);
                 variant < reference_element:: last_variant_by_dimension(map_dim); variant++) {
    for (typename std::list<std::pair<element_type,size_type> >::const_iterator
	  iter = gamma_side_list[variant].begin(),
	  last = gamma_side_list[variant].end();
	  iter != last; iter++, idx_iter++) {
      *idx_iter = (*iter).second;
    }
  }
}
struct to_solve {
  size_t variant, S_ige, k, dis_i;
  to_solve (size_t variant1, size_t S_ige1=0, size_t k1=0, size_t dis_i1=0) 
   : variant(variant1), S_ige(S_ige1), k(k1), dis_i(dis_i1) {}
};
template <class T, class M>
geo_basic<T,M>
level_set_internal (
    const field_basic<T,M>&      fh,
    const level_set_option& opt,
    std::vector<size_t>&         bnd_dom_ie_list,
    disarray<size_t,M>&          sid_ie2bnd_ie)
{
  using namespace std;
  using namespace details;
  typedef geo_element::size_type size_type;
  level_set_epsilon = opt.epsilon; // set global variable
  fh.dis_dof_update();
  const geo_basic<T,M>&   lambda = fh.get_geo();
  const space_basic<T,M>& Xh     = fh.get_space();
  check_macro(lambda.order() == 1,        "Only order=1 level set mesh supported");
  check_macro(fh.get_approx() == "P1", "Only P1 level set function supported");
  size_type order = 1;
  std::vector<size_type> dis_idof;
  std::vector<T>     f;
  size_type d = lambda.dimension();
  size_type map_dim = d-1;
  heap_allocator<size_type> alloc;
  std::array<list<pair<element_type,size_type> >, 
               reference_element::max_variant> gamma_side_list;
  list<point_basic<T> >                        gamma_node_list;
  std::vector<size_type>       j(d+1);
  std::vector<point_basic<T> > x(d+1);
  const size_type not_marked = numeric_limits<size_t>::max();
  const size_type undef      = numeric_limits<size_t>::max();

  distributor                 node_ownership = lambda.sizes().node_ownership;
  size_type                   first_dis_inod = node_ownership.first_index();
  disarray<size_type>         marked_node (node_ownership, not_marked);
  disarray<size_type>         extern_node (node_ownership, not_marked);
  set<size_type>              ext_marked_node_idx;
  list<to_solve>              node_to_solve;

  distributor                 edge_ownership = lambda.sizes().ownership_by_dimension[1];
  size_type                   first_dis_iedg = edge_ownership.first_index();
  disarray<size_type>         marked_edge      (edge_ownership, not_marked);
  disarray<std::pair<size_type,point_basic<T> > > 
                              extern_edge      (edge_ownership, std::make_pair(not_marked,point_basic<T>()));
  set<size_type>              ext_marked_edge_idx;
  list<to_solve>              edge_to_solve;

  distributor                 face_ownership = lambda.sizes().ownership_by_dimension[2];
  size_type                   first_dis_ifac = face_ownership.first_index();
  disarray<size_type>         marked_face (face_ownership, not_marked);

  communicator comm = node_ownership.comm();
  size_type my_proc = comm.rank();

  bnd_dom_ie_list.resize(0);

  // ------------------------------------------------------------
  // 1) loop on lambda & build intersection sides
  // ------------------------------------------------------------
  for (size_type ie = 0, ne = lambda.size(), bnd_ie = 0; ie < ne; ie++) {
    // ---------------------------------------------------------
    // 1.1) fast check if there is an intersection:
    // ---------------------------------------------------------
    const geo_element& K = lambda [ie];
    Xh.dis_idof (K, dis_idof);
    f.resize (dis_idof.size());
    for (size_type loc_idof = 0, loc_ndof = dis_idof.size(); loc_idof < loc_ndof; loc_idof++) {
      f [loc_idof] = fh.dis_dof (dis_idof[loc_idof]);
    }
    bool do_intersect = false;
    switch (K.variant()) {
      case reference_element::t:
        do_intersect = belongs_to_band_t (f);
        break;
      case reference_element::T: {
        do_intersect = belongs_to_band_T (f);
        break;
      }
      default :
        error_macro("level set: element type `" << K.name() << "' not yet supported");
    }
    if (! do_intersect) continue;
    bnd_dom_ie_list.push_back (ie);
    // ---------------------------------------------------------
    // 1.2) compute the intersection
    // ---------------------------------------------------------
    x.resize   (dis_idof.size());
    for (size_type loc_idof = 0, loc_ndof = dis_idof.size(); loc_idof < loc_ndof; loc_idof++) {
      size_type loc_inod = loc_idof; // assume here that fh has an isoparametric approx
      size_type dis_inod = K [loc_inod];
      x [loc_idof] = lambda.dis_node(dis_inod);
    }
    element_type S (alloc);
    switch (K.variant()) {
      case reference_element::t: {
        // ---------------------------------------------------------
        // 1.2.a) triangle -> 1 edge
        // ---------------------------------------------------------
        point_basic<T> a, b;
        T length;
        subcompute_matrix_t (x, f, j, a, b, length);
        if (is_zero(f[j[1]]) && is_zero(f[j[2]])) {
          // ---------------------------------------------------------
          // 1.2.a.1) edge {j1,j2} is included in the bbox mesh
          // ---------------------------------------------------------
          for (size_type k = 0; k < 2; k++) {
            size_type dis_inod = K[j[k+1]];
	    if (node_ownership.is_owned (dis_inod)) {
              size_type inod = dis_inod - first_dis_inod;
	      if (marked_node [inod] == not_marked) {
	        marked_node [inod] = gamma_node_list.size();
                gamma_node_list.push_back (lambda.node(inod));
              }
            } else {
              // dis_inod is owned by another partition jproc:
              // if there is another neighbour element K' from the same partition jproc
              //     then jproc will insert dis_inod in gamma
	      //     so, there is nothing to do here
              // otherwise, dis_inod will be orphan in jproc and will be re-affected to my_proc
  	      extern_node.dis_entry (dis_inod) = my_proc;
            }
          }
	  size_type loc_iedg = edge_t_iloc (j[1], j[2]);
	  size_type dis_iedg = K.edge (loc_iedg);
	  if (edge_ownership.is_owned (dis_iedg)) {
            size_type iedg = dis_iedg - first_dis_iedg;
	    if (marked_edge [iedg] == not_marked) {
              S.reset(reference_element::e, order);
	      marked_edge [iedg] = gamma_side_list [S.variant()].size();
	      size_type S_ige = marked_edge [iedg];
              for (size_type k = 0; k < S.n_node(); k++) {
                size_type dis_inod = K[j[k+1]];
	        if (node_ownership.is_owned (dis_inod)) {
                  size_type inod = dis_inod - first_dis_inod;
                  S[k] = marked_node [inod];
                } else {
	          S[k] = undef;
		  node_to_solve.push_back (to_solve (S.variant(), S_ige, k, dis_inod));
                  ext_marked_node_idx.insert (dis_inod);
                }
              }
              gamma_side_list [S.variant()].push_back (make_pair(S,bnd_ie));
            }
          } else {
            // intersection is an edge of lambda that is owned by another partition jproc
            // the neighbour element K' across this e"dge belongs to the same partition jproc
            // and will insert dis_iedg in gamma
	    // so, there is nothing to do here
          }
        } else {
          // ---------------------------------------------------------
          // 1.2.a.2) edge {j1,j2} is interior to the triangle
          // ---------------------------------------------------------
          S.reset(reference_element::e, order);
	  size_type S_ige = gamma_side_list [S.variant()].size ();
	  point_basic<T> xx[2] = {a,b};
          for (size_type k = 0; k < 2; k++) {
            if (! is_zero(f[j[k+1]]) && ! is_zero(f[j[0]])) {
              // xk is inside edge {j0,j[k+1]} of triangle K:
              size_type loc_iedg = edge_t_iloc (j[0], j[k+1]);
              size_type dis_iedg = K.edge (loc_iedg);
	      if (edge_ownership.is_owned (dis_iedg)) {
                size_type iedg = dis_iedg - first_dis_iedg;
		if (marked_edge [iedg] == not_marked) {
	          marked_edge [iedg] = gamma_node_list.size();
                  gamma_node_list.push_back (xx[k]);
                }
	        S[k] = marked_edge [iedg];
              } else {
	        S[k] = undef;
		edge_to_solve.push_back (to_solve (S.variant(), S_ige, k, dis_iedg));
                ext_marked_edge_idx.insert (dis_iedg);
  	        extern_edge.dis_entry (dis_iedg) = std::make_pair (my_proc, xx[k]);
              }
            } else { // xk is at edge boundary: a node of the 2d mesh
              size_type dis_inod = (!is_zero(f[j[0]])) ? K[j[k+1]] : K[j[0]];
	      if (node_ownership.is_owned (dis_inod)) {
                size_type inod = dis_inod - first_dis_inod;
	        if (marked_node [inod] == not_marked) {
	          marked_node [inod] = gamma_node_list.size();
                  gamma_node_list.push_back (lambda.node(inod));
                }
	        S[k] = marked_node [inod];
              } else {
	        S[k] = undef;
		node_to_solve.push_back (to_solve (S.variant(), S_ige, k, dis_inod));
                ext_marked_node_idx.insert (dis_inod);
  	        extern_node.dis_entry (dis_inod) = my_proc;
              }
            }
          }
          // S[0] == S[1] when is_zero(f[j[0]]) but f[j[0]] != 0, i.e. precision pbs
          check_macro (S[0] != S[1] || S[0] == undef, "degenerate 2d intersection");
          gamma_side_list [S.variant()].push_back (make_pair(S,bnd_ie));
        }
        break;
      }
      case reference_element::T: {
        // ---------------------------------------------------------
        // 1.2.b) tetrahedron -> 1 triangle or 1 quadrangle
        // ---------------------------------------------------------
        quadruplet q;
        point_basic<T> a, b, c, d;
        T aire;
        if (!intersection_is_quadrilateral_T (f, q)) {
          subcompute_matrix_triangle_T (x, f, j, a, b, c, aire);
          if (is_zero(f[j[1]]) && is_zero(f[j[2]]) && is_zero(f[j[3]])) {
            // the full face {j1,j2,j3} is included in the surface mesh:
            for (size_type k = 0; k < 3; k++) {
              size_type dis_inod = K[j[k+1]];
	      if (node_ownership.is_owned (dis_inod)) {
                size_type inod = dis_inod - first_dis_inod;
	        if (marked_node [inod] == not_marked) {
	          marked_node [inod] = gamma_node_list.size();
                  gamma_node_list.push_back (lambda.node(inod));
                }
              } else {
                // dis_inod is owned by another partition jproc:
                // if there is another neighbour element K' from the same partition jproc
                //     then jproc will insert dis_inod in gamma
	        //     so, there is nothing to do here
                // otherwise, dis_inod will be orphan in jproc and will be re-affected to my_proc
  	        extern_node.dis_entry (dis_inod) = my_proc;
              }
            }
	    size_type loc_ifac = face_T_iloc (j[1], j[2], j[3]);
	    size_type dis_ifac = K.face (loc_ifac);
	    if (face_ownership.is_owned (dis_ifac)) {
              size_type ifac = dis_ifac - first_dis_ifac;
	      if (marked_face [ifac] == not_marked) {
                S.reset(reference_element::t, order);
	        marked_face [ifac] = gamma_side_list [S.variant()].size();
	        size_type S_ige = marked_face [ifac];
                for (size_type k = 0; k < S.n_node(); k++) {
                  size_type dis_inod = K[j[k+1]];
	          if (node_ownership.is_owned (dis_inod)) {
                    size_type inod = dis_inod - first_dis_inod;
                    S[k] = marked_node [inod];
                  } else {
	            S[k] = undef;
		    node_to_solve.push_back (to_solve (S.variant(), S_ige, k, dis_inod));
                    ext_marked_node_idx.insert (dis_inod);
                  }
                }
                gamma_side_list [S.variant()].push_back (make_pair(S,bnd_ie));
              } else {
	        // the side will be inserted by the neighbour of K in another partition
	        // so, nothing to do
              }
            } else {
              // intersection is a face of lambda that is owned by another partition jproc
              // the neighbour element K' across this face belongs to the same partition jproc
              // and will insert dis_ifac in gamma
	      // so, there is nothing to do here
            }
          } else {
            // create the new face {j1,j2,j3} by intersections:
            S.reset(reference_element::t, order);
	    size_type S_ige = gamma_side_list [S.variant()].size();
	    point_basic<T> xx[3] = {a,b,c};
            for (size_type k = 0; k < 3; k++) {
              if (! is_zero(f[j[k+1]]) && ! is_zero(f[j[0]])) {
                // xk is inside edge {j0,j[k+1]} of triangle K:
                size_type loc_iedg = edge_T_iloc (j[0], j[k+1]);
                size_type dis_iedg = K.edge (loc_iedg);
	        if (edge_ownership.is_owned (dis_iedg)) {
                  size_type iedg = dis_iedg - first_dis_iedg;
	          if (marked_edge [iedg] == not_marked) {
	            marked_edge [iedg] = gamma_node_list.size();
                    gamma_node_list.push_back (xx[k]);
                  }
	          S[k] = marked_edge [iedg];
                } else {
	          S[k] = undef;
		  edge_to_solve.push_back (to_solve (S.variant(), S_ige, k, dis_iedg));
                  ext_marked_edge_idx.insert (dis_iedg);
  	          extern_edge.dis_entry (dis_iedg) = std::make_pair (my_proc, xx[k]);
                }
              } else { // xk is at edge boundary: a node of the 2d mesh
                size_type dis_inod = (!is_zero(f[j[0]])) ? K[j[k+1]] : K[j[0]];
	        if (node_ownership.is_owned (dis_inod)) {
                  size_type inod = dis_inod - first_dis_inod;
	          if (marked_node [inod] == not_marked) {
	            marked_node [inod] = gamma_node_list.size();
                    gamma_node_list.push_back (lambda.node(inod));
                  }
	          S[k] = marked_node [inod];
                } else {
	          S[k] = undef;
		  node_to_solve.push_back (to_solve (S.variant(), S_ige, k, dis_inod));
                  ext_marked_node_idx.insert (dis_inod);
  	          extern_node.dis_entry (dis_inod) = my_proc;
                }
              }
            }
            // S[0] == S[j] when is_zero(f[j[0]]) but f[j[0]] != 0, i.e. precision pbs
            check_macro ((S[0] != S[1] || S[0] == undef) &&
                         (S[1] != S[2] || S[1] == undef) &&
                         (S[2] != S[0] || S[2] == undef), "degenerate 3d intersection");
            gamma_side_list [S.variant()].push_back (make_pair(S,bnd_ie));
          }
        } else {
          // create the new quadri face by intersections:
          subcompute_matrix_quadrilateral_T (x, f, q, a, b, c, d, aire);
          S.reset (reference_element::q, order);
	  size_type S_ige  = gamma_side_list[reference_element::q].size();
	  size_type S1_ige = gamma_side_list[reference_element::t].size(); // or S1+S2=two 't'
	  size_type S2_ige = S1_ige + 1;
          size_type iloc_S1[4] = {0, 1, 2, undef}; // the way to slip a q into 2 t
          size_type iloc_S2[4] = {0, undef, 1, 2};
	  point_basic<T>  xx[4] = {a,b,d,c};
	  size_type s[4] = {q[0],q[2],q[1],q[3]};
          for (size_type k = 0; k < 4; k++) {
            size_type k1 = (k+1) % 4;
            if (! is_zero(f[s[k]]) && ! is_zero(f[s[k1]])) {
              // xk is inside edge {j0,j[k+1]} of triangle K:
              size_type loc_iedg = edge_T_iloc (s[k], s[k1]);
              size_type dis_iedg = K.edge (loc_iedg);
	      if (edge_ownership.is_owned (dis_iedg)) {
                size_type iedg = dis_iedg - first_dis_iedg;
	        if (marked_edge [iedg] == not_marked) {
	          marked_edge [iedg] = gamma_node_list.size();
                  gamma_node_list.push_back (xx[k]);
                }
	        S[k] = marked_edge [iedg]; 
              } else {
	        S[k] = undef;
	        if (!opt.split_to_triangle) {
		  edge_to_solve.push_back (to_solve (reference_element::q, S_ige, k, dis_iedg));
		} else {
		  if (iloc_S1[k] != undef) edge_to_solve.push_back (to_solve (reference_element::t, S1_ige, iloc_S1[k], dis_iedg));
		  if (iloc_S2[k] != undef) edge_to_solve.push_back (to_solve (reference_element::t, S2_ige, iloc_S2[k], dis_iedg));
                }
                ext_marked_edge_idx.insert (dis_iedg);
  	        extern_edge.dis_entry (dis_iedg) = std::make_pair (my_proc, xx[k]);
              }
            } else { // xk is at edge boundary: a node of the 2d mesh
              size_type dis_inod = is_zero(f[s[k]]) ? K[s[k]] : K[s[k1]];
	      if (node_ownership.is_owned (dis_inod)) {
                size_type inod = dis_inod - first_dis_inod;
	        if (marked_node [inod] == not_marked) {
	          marked_node [inod] = gamma_node_list.size();
                  gamma_node_list.push_back (lambda.node(inod));
                }
	        S[k] = marked_node [inod];
              } else {
	        S[k] = undef;
	        if (!opt.split_to_triangle) {
		  node_to_solve.push_back (to_solve (reference_element::q, S_ige, k, dis_inod));
		} else {
		  if (iloc_S1[k] != undef) node_to_solve.push_back (to_solve (reference_element::t, S1_ige, iloc_S1[k], dis_inod));
		  if (iloc_S2[k] != undef) node_to_solve.push_back (to_solve (reference_element::t, S2_ige, iloc_S2[k], dis_inod));
                }
                ext_marked_node_idx.insert (dis_inod);
              }
            }
          }
	  if (!opt.split_to_triangle) {
            check_macro ((S[0] != S[1] || S[0] == undef) &&
                         (S[1] != S[2] || S[1] == undef) &&
                         (S[2] != S[3] || S[2] == undef) &&
                         (S[3] != S[0] || S[3] == undef), "degenerate 3d intersection");
            // S[0] == S[j] when is_zero(f[j[0]]) but f[j[0]] != 0, i.e. precision pbs
            gamma_side_list [S.variant()].push_back (make_pair(S,bnd_ie));
	  } else {
	    // split quadri into 2 triangles
	    // one K -> two (S1,S2) faces: table element2face may return a pair of size_t
	    // but S1-> and S2->K only is required during assembly
            element_type S1 (alloc);
            S1.reset (reference_element::t, order);
	    for (size_type k = 0; k < 4; k++) {
	      if (iloc_S1[k] != undef) S1 [iloc_S1[k]] = S[k];
            }
            check_macro ((S1[0] != S1[1] || S1[0] == undef) &&
                         (S1[1] != S1[2] || S1[1] == undef) &&
                         (S1[2] != S1[0] || S1[2] == undef), "degenerate 3d intersection");
            gamma_side_list [S1.variant()].push_back (make_pair(S1,bnd_ie));

            element_type S2 (alloc);
            S2.reset (reference_element::t, order);
	    for (size_type k = 0; k < 4; k++) {
	      if (iloc_S2[k] != undef) S2 [iloc_S2[k]] = S[k];
            }
            check_macro ((S2[0] != S2[1] || S2[0] == undef) &&
                         (S2[1] != S2[2] || S2[1] == undef) &&
                         (S2[2] != S2[0] || S2[2] == undef), "degenerate 3d intersection");
            gamma_side_list [S2.variant()].push_back (make_pair(S2,bnd_ie));
          }
        } // if-else
        break;
      }
      default : {
        error_macro("level set intersection: element not yet implemented: " << K.name());
      }
    }
    bnd_ie++;
  }
  extern_node.dis_entry_assembly();
  extern_edge.dis_entry_assembly();
  // ------------------------------------------------------------
  // 2) solve orphan nodes, if any
  // ------------------------------------------------------------
  // 2.1.a) re-affect orphan node to another process where gamma use it
  distributor comm_ownership (comm.size(), comm, 1);
  disarray<index_set,M> orphan_node (comm_ownership, index_set());
  for (size_type inod = 0, nnod = node_ownership.size(); inod < nnod; inod++) {
    if (!(extern_node [inod] != not_marked && marked_node [inod] == not_marked)) continue;
    // inod is orphan in this proc: not used for gamma (but used for lambda)
    // re-affect it to a process that use it for gamma
    size_type iproc = extern_node[inod];
    size_type dis_inod = first_dis_inod + inod;
    index_set dis_inod_set;
    dis_inod_set.insert (dis_inod);
    orphan_node.dis_entry (iproc) += dis_inod_set;
  }
  orphan_node.dis_entry_assembly();

  // 2.1.b) re-affect orphan edge to another process where gamma use it
  // there could be orphan edge in 2d if lambda is a part of a regular mesh:
  // => a boundary edge can belong to another proc. This is not yet handled
  disarray<index_set,M> orphan_edge (comm_ownership, index_set());
  for (size_type iedg = 0, nedg = edge_ownership.size(); iedg < nedg; iedg++) {
    if (!(extern_edge [iedg].first != not_marked && marked_edge [iedg] == not_marked)) continue;
    // iedg is orphan in this proc: not used for gamma (but used for lambda)
    // re-affect it to a process that use it for gamma
    size_type iproc = extern_edge[iedg].first;
    size_type dis_iedg = first_dis_iedg + iedg;
    index_set dis_iedg_set;
    dis_iedg_set.insert (dis_iedg);
    orphan_edge.dis_entry (iproc) += dis_iedg_set;
  }
  orphan_edge.dis_entry_assembly ();
  check_macro (orphan_edge[0].size() == 0, "unexpected orphan edges");

  // 2.2) count total nodes used by gamma
  size_type  orphan_gamma_nnod = orphan_node[0].size();
  size_type regular_gamma_nnod = gamma_node_list.size();
  size_type         gamma_nnod = regular_gamma_nnod + orphan_gamma_nnod;
  distributor gamma_node_ownership (distributor::decide, comm, gamma_nnod);
  // 2.3) shift marked_node & marked_edge from gamma_inod local count to gamma_dis_inod one
  size_type gamma_first_dis_inod = gamma_node_ownership.first_index();
  for (size_type inod = 0, nnod = node_ownership.size(); inod < nnod; inod++) {
    if (marked_node [inod] == not_marked) continue;
    marked_node [inod] += gamma_first_dis_inod;
  }
  for (size_type iedg = 0, nedg = edge_ownership.size(); iedg < nedg; iedg++) {
    if (marked_edge [iedg] == not_marked) continue;
    marked_edge [iedg] += gamma_first_dis_inod;
  }
  // 2.4) append orphan node to regular one in gamma_node_list and set marked_node
  for (index_set::const_iterator
        iter = orphan_node[0].begin(),
        last = orphan_node[0].end(); iter != last; ++iter) {
    size_type dis_inod = *iter;
    marked_node.dis_entry(dis_inod) = gamma_first_dis_inod + gamma_node_list.size();
    gamma_node_list.push_back (lambda.dis_node(dis_inod));
  }
  marked_node.dis_entry_assembly();
  marked_edge.dis_entry_assembly();
  // ------------------------------------------------------------
  // 3) convert lists to fixed size distributed arrays
  // ------------------------------------------------------------
  disarray<point_basic<T>, M>                                            gamma_node;
  std::array<disarray<element_type,M>, reference_element::max_variant> gamma_side;
  gamma_list2disarray (gamma_node_list, gamma_side_list, comm, d, gamma_node, gamma_side, sid_ie2bnd_ie);
  // ------------------------------------------------------------
  // 4) replace inod to dis_inod in element lists
  // ------------------------------------------------------------
  for (size_type variant = reference_element::first_variant_by_dimension(map_dim);
                 variant < reference_element:: last_variant_by_dimension(map_dim); variant++) {
    for (size_type ige = 0, nge = gamma_side[variant].size(); ige < nge; ige++) {
      element_type& S = gamma_side[variant][ige];
      for (size_type loc_inod = 0, loc_nnod = S.n_node(); loc_inod < loc_nnod; loc_inod++) {
        if (S[loc_inod] == undef) continue; // external node, will be solved later
        S[loc_inod] += gamma_first_dis_inod;
      }
    }
  }
  // ----------------------------------------------------------------
  // 5) solve intersection that are located on external edges & nodes
  // ----------------------------------------------------------------
  marked_node.set_dis_indexes (ext_marked_node_idx);
  for (list<to_solve>::const_iterator iter = node_to_solve.begin(),
                                      last = node_to_solve.end(); iter != last; iter++) {
    const to_solve& x = *iter;
    element_type& S = gamma_side[x.variant][x.S_ige];
    check_macro (S[x.k] == undef, "external index already solved");
    size_type dis_inod = x.dis_i;
    size_type iproc = node_ownership.find_owner(dis_inod);
    size_type gamma_dis_inod = marked_node.dis_at (dis_inod);
    S[x.k] = gamma_dis_inod;
  }

  marked_edge.set_dis_indexes (ext_marked_edge_idx);
  for (list<to_solve>::const_iterator iter = edge_to_solve.begin(),
                                      last = edge_to_solve.end(); iter != last; iter++) {
    const to_solve& x = *iter;
    element_type& S = gamma_side[x.variant][x.S_ige];
    check_macro (S[x.k] == undef, "external index already solved");
    size_type dis_iedg = x.dis_i;
    size_type iproc = edge_ownership.find_owner(dis_iedg);
    size_type gamma_dis_inod = marked_edge.dis_at (dis_iedg);
    S[x.k] = gamma_dis_inod;
  }
  // ------------------------------------------------------------
  // 6) convert intersection to geo
  // ------------------------------------------------------------
  geo_basic<T,M> gamma (lambda, gamma_node, gamma_side);
  return gamma;
}
template <class T, class M>
geo_basic<T,M>
level_set (
    const field_basic<T,M>&      fh,
    const level_set_option& opt)
{
  typedef geo_element::size_type size_type;
  std::vector<size_type> bnd_dom_ie_list;
  disarray<size_type,M>     sid_ie2bnd_ie;
  return level_set_internal (fh, opt, bnd_dom_ie_list, sid_ie2bnd_ie);
}
// ----------------------------------------------------------------------------
// instanciation in library
// ----------------------------------------------------------------------------
#define _RHEOLEF_instanciation(T,M)                                                     \
template geo_basic<T,M> level_set_internal (						\
    const field_basic<T,M>&,								\
    const level_set_option&,							\
    std::vector<size_t>&,								\
    disarray<size_t,M>&);								\
template geo_basic<T,M> level_set (							\
    const field_basic<T,M>&,								\
    const level_set_option&);

_RHEOLEF_instanciation(Float,sequential)
#ifdef _RHEOLEF_HAVE_MPI
_RHEOLEF_instanciation(Float,distributed)
#endif // _RHEOLEF_HAVE_MPI

} // namespace