Logo Search packages:      
Sourcecode: ecs version File versions  Download package

ecs_vec_def_rc.c

/*============================================================================
 *  Définitions des fonctions
 *   associées aux structures `ecs_vec_int_t' et `ecs_vec_real_t' décrivant
 *   les vecteurs indexés entier et réel
 *   et propres aux vecteurs indexés
 *      liés aux champs principaux de type "definition"
 *  Ces fonctions participent à la fonctionnalité de "recollement conforme"
 *============================================================================*/

/*
  This file is part of the Code_Saturne Preprocessor, element of the
  Code_Saturne CFD tool.

  Copyright (C) 1999-2007 EDF S.A., France

  contact: saturne-support@edf.fr

  The Code_Saturne Preprocessor 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.

  The Code_Saturne Preprocessor 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 the Code_Saturne Preprocessor; if not, write to the
  Free Software Foundation, Inc.,
  51 Franklin St, Fifth Floor,
  Boston, MA  02110-1301  USA
*/


/*----------------------------------------------------------------------------
 *  Fichiers `include' librairie standard C
 *----------------------------------------------------------------------------*/

#include <assert.h>
#include <math.h>   /* sqrt() */
#include <stdlib.h>


/*----------------------------------------------------------------------------
 *  Fichiers `include' système ou BFT
 *----------------------------------------------------------------------------*/

#include <bft_error.h>
#include <bft_mem.h>
#include <bft_printf.h>


/*----------------------------------------------------------------------------
 *  Fichiers `include' visibles du  paquetage global "Utilitaire"
 *----------------------------------------------------------------------------*/

#include "ecs_chaine_glob.h"
#include "ecs_chrono.h"
#include "ecs_def.h"
#include "ecs_param_rc_glob.h"
#include "ecs_param_perio_glob.h"
#include "ecs_tab.h"


/*----------------------------------------------------------------------------
 *  Fichiers `include' visibles des paquetages visibles
 *----------------------------------------------------------------------------*/


/*----------------------------------------------------------------------------
 *  Fichiers `include' visibles du  paquetage courant
 *----------------------------------------------------------------------------*/

#include "ecs_vec_int.h"
#include "ecs_vec_real.h"
#include "ecs_vec_real_tri.h"
#include "ecs_vec_def.h"
#include "ecs_vec_def_perio.h"


/*----------------------------------------------------------------------------
 *  Fichier  `include' du  paquetage courant associe au fichier courant
 *----------------------------------------------------------------------------*/

#include "ecs_vec_def_rc.h"


/*----------------------------------------------------------------------------
 *  Fichiers `include' prives   du  paquetage courant
 *----------------------------------------------------------------------------*/

#include "ecs_vec_int_priv.h"
#include "ecs_vec_real_priv.h"


/*============================================================================
 *                       Macros globales au fichier
 *============================================================================*/

enum {
  X,
  Y,
  Z
} ;


#define ECS_LOC_RC_NBR_MOY_INTER_ARE  1
#define ECS_LOC_RC_NBR_INTER_ARE_INI  1024

#define ECS_LOC_DISTANCE(som1, som2)                                   \
  (sqrt(((*(som1    ) - *(som2    )) * (*(som1    ) - *(som2    ))) +  \
        ((*(som1 + 1) - *(som2 + 1)) * (*(som1 + 1) - *(som2 + 1))) +  \
        ((*(som1 + 2) - *(som2 + 2)) * (*(som1 + 2) - *(som2 + 2)))  )   )

#define ECS_LOC_MODULE_CARRE(vect) \
(vect[X] * vect[X] + vect[Y] * vect[Y] + vect[Z] * vect[Z])

#define ECS_LOC_PRODUIT_SCALAIRE(vect1, vect2)                        ( \
  vect1[X] * vect2[X] + vect1[Y] * vect2[Y] + vect1[Z] * vect2[Z] )

#define ECS_LOC_PRODUIT_VECTORIEL(prod_vect, vect1, vect2) ( \
prod_vect[X] = vect1[Y] * vect2[Z] - vect2[Y] * vect1[Z],    \
prod_vect[Y] = vect2[X] * vect1[Z] - vect1[X] * vect2[Z],    \
prod_vect[Z] = vect1[X] * vect2[Y] - vect2[X] * vect1[Y]   )


/*============================================================================
 *  Structures locales
 *============================================================================*/

/* Intersection arête-arête */

/*
  Afin de pouvoir parcourir les intersections par arête selon la coordonnée
  linéique croissante, on stocke les intersections calculées dans deux
  tableaux ; l'un contient le numéro d'arête et à la coordonnée linéique
  correspondante (on a donc deux entrées par intersection), et l'autre
  contient l'indice de chacune des deux entrées "arêtes" dans le premier.
  Cette association de structures permet le tri du tableau des arêtes par
  numéro d'arête et coordonnée linéique d'intersection, nécessaire à
  l'opérateur de découpage des arêtes selon leurs intersections et fusion des
  sommets confondus.
*/

typedef struct _ecs_rc_inter_are_are_t {
  ecs_size_t                 are ;            /* Arête                        */
  ecs_size_t                 i_inter ;        /* Indice de l'intersection     */
  ecs_real_t                 s_inter ;        /* Coordonnée s                 */
} ecs_rc_inter_are_are_t ;

typedef struct _ecs_rc_inter_are_inter_t {
  ecs_int_t                  ind_0 ;          /* Intersection 0               */
  ecs_int_t                  ind_1 ;          /* Intersection 1               */
} ecs_rc_inter_are_inter_t ;


typedef struct _ecs_rc_lis_inter_are_t {
  size_t                     nbr ;            /* Nombre d'éléments            */
  size_t                     nbr_max ;        /* Nombre d'éléments maximal    */
  ecs_rc_inter_are_are_t    *are ;            /* Coordonnées locales          */
  ecs_rc_inter_are_inter_t  *inter ;          /* Liste des intersections      */
} ecs_rc_lis_inter_are_t ;


/*============================================================================
 *                       Prototypes de fonctions privées
 *============================================================================*/

/*----------------------------------------------------------------------------
 *  Fonction qui renvoie l'ordonnancement des éléments
 *   en fonction de leur définition (sans modifier cette définition)
 *----------------------------------------------------------------------------*/

static ecs_tab_int_t ecs_loc_vec_def_rc__trie_elt
(
 const ecs_vec_real_t  *const this_vec_def_real,
 const ecs_real_t             tolerance
) ;


/*----------------------------------------------------------------------------
 * Fonction qui calcule les normales unitaires aux faces
 *----------------------------------------------------------------------------*/

void ecs_loc_vec_def_rc__normales
(
 const ecs_tab_int_t  *const tab_fac_select,
 const ecs_vec_int_t  *const vec_fac_def,
 const ecs_vec_int_t  *const vec_are_def,
 const ecs_vec_real_t *const vec_som_def,
       float          *const normale_fac
);


/*----------------------------------------------------------------------------
 *  Fonction qui découpe les arêtes
 *----------------------------------------------------------------------------*/

static ecs_vec_int_t * ecs_loc_vec_def_rc__decoup_are
(
 ecs_vec_int_t           *const vec_def_are,      /* <-> Déf. arêtes          */
 ecs_vec_real_t          *const vec_def_som,      /* <-> Déf. sommets         */
 ecs_rc_lis_inter_are_t         liste_inter_are,  /* <-> Déf. intersections   */
 ecs_tab_int_t           *const tab_equiv_som,    /* <-  Fusions de sommets   */
 ecs_tab_real_t          *const dist_max_som      /* <-> Tolérances           */
) ;


/*----------------------------------------------------------------------------
 *  Fonction qui détermine les sommets à fusionner pour le découpage
 * des arêtes ; On renvoie 0 si aucune transitivité de fusion n'a été
 * détectée, et 1 si une ou plusieurs transitivités ont du être éliminées
 * par réduction de tolérance ;
 *----------------------------------------------------------------------------*/

static ecs_int_t ecs_loc_vec_def_rc__fusion_inter
(
 ecs_vec_int_t           *const vec_def_are,             /*  -> Déf. arêtes  */
 ecs_vec_real_t          *const vec_def_som,             /*  -> Déf. sommets */
 ecs_rc_lis_inter_are_t         liste_inter_are,         /*  -> Intersects.  */
 ecs_int_t               *const som_inter,               /*  -> Sommets/int. */
 ecs_tab_real_t          *const dist_max_som,            /*  -> Tolérances   */
 ecs_int_t                      ind_are_inter_deb,       /*  -> Début arête  */
 size_t                         nbr_inter_are,           /*  -> Nb. inter    */
 ecs_int_t               *const equiv_som_are            /* <-> Fusions      */
);


/*----------------------------------------------------------------------------
 *  Fonction qui supprime les éventuelles arêtes dégénérées
 *----------------------------------------------------------------------------*/

static void  ecs_loc_vec_def_rc__nettoie_are
(
 ecs_vec_int_t           *const vec_def_are,      /* <-> Déf. arêtes          */
 ecs_vec_int_t           *const vec_are_old_new   /* <-> Renumérotation       */
) ;


/*----------------------------------------------------------------------------
 *  Fonction qui enlève les éventuels coins dégénérés des faces
 *----------------------------------------------------------------------------*/

static void  ecs_loc_vec_def_rc__nettoie_fac
(
 ecs_vec_int_t           *const vec_def_fac       /* <-> Déf. faces           */
) ;


/*----------------------------------------------------------------------------
 *  Construction des nouvelles arêtes decoupées par les intersections
 *
 *  La liste des sommets finaux issus d'une intersection ou d'une fusion
 *  est construite (et allouée) ici ; la structure liste_som_new
 *  correspondante est donc vide en entrée ;
 *----------------------------------------------------------------------------*/

static void ecs_loc_vec_def_rc__def_are_new
(
       ecs_vec_int_t         *const vec_def_fac,
       ecs_vec_int_t       * *const vec_def_are,
       ecs_vec_real_t        *const vec_def_som,
       ecs_vec_int_t       * *const vec_are_old_new,
       ecs_vec_int_t       * *const vec_som_old_new,
       ecs_vec_int_t         *const vec_fac_vis,
       ecs_vec_int_t       * *const vec_fac_perio,
       ecs_vec_int_t       * *const vec_are_perio,
       ecs_vec_int_t       * *const vec_som_perio,
 const ecs_tab_int_t         *const tab_fac_select,
       ecs_tab_int_t         *const liste_som_new,
 const ecs_param_rc_t               param_rc
) ;


/*----------------------------------------------------------------------------
 *  Fonction qui construit les nouvelles faces
 *----------------------------------------------------------------------------*/

static ecs_vec_int_t * ecs_loc_vec_def_rc__def_fac_new
(
 const ecs_vec_int_t    *const vec_def_fac,
 const ecs_vec_int_t    *const vec_def_are,
 const ecs_vec_real_t   *const vec_def_som,
 const ecs_vec_int_t    *const vec_fac_vis,
 const float            *const normale_fac,
       ecs_vec_int_t  * *const vec_fac_old_new,
 const ecs_tab_int_t    *const tab_fac_select,
       ecs_tab_int_t    *const liste_fac_err,
       size_t                  nbr_fac_ini,
 const ecs_real_t              epsilon_plan
) ;


/*----------------------------------------------------------------------------
 *  Fonction qui construit des nouvelles faces issues d'une face donnée
 *
 *  Renvoie :
 *    0 en cas de succès ;
 *    1 en cas de "cycle ouvert" (pas d'arête suivante acceptable) ;
 *    2 en cas d'arête parcourue deux fois pour une même face reconstruite ;
 *    3 en cas de boucle (nombre de sous-faces trop important et irréaliste).
 *----------------------------------------------------------------------------*/

static ecs_int_t ecs_loc_vec_def_rc__decoup_fac
(
 const ecs_vec_int_t    *const vec_def_fac,
 const ecs_vec_int_t    *const vec_def_are,
 const ecs_vec_real_t   *const vec_def_som,
 const float            *const normale_fac,
 const ecs_int_t        *const masque_fac_select,
 const ecs_vec_int_t    *const vec_fac_vis,
 const ecs_vec_int_t    *const vec_def_som_are,
 const ecs_vec_int_t    *const vec_def_are_fac,
       ecs_vec_int_t    *const vec_def_fac_new,
       ecs_vec_int_t    *const vec_fac_old_new,
 const ecs_int_t               ifac ,
 const ecs_real_t              epsilon_plan,
       size_t                  nbr_max_fac_dec,
       ecs_tab_int_t    *const liste_are_tete,
       ecs_tab_int_t    *const liste_are_ext,
       ecs_tab_int_t    *const liste_are_int,
       ecs_tab_int_t    *const liste_are_fac,
       size_t           *const max_pos_fac_new,
       size_t           *const max_val_fac_new
) ;


/*----------------------------------------------------------------------------
 *  Fonction qui détermine la liste des faces participant au recollement
 *  lorsqu'une définition des faces visibles par d'autres est déjà connue.
 *  Il n'est pas nécessaire que deux faces en vis-à-vis partiel se
 *  connaissent ; il suffit que l'une d'elles connaîsse l'autre (cas d'une
 *  filiation par exemple), même si une connaissance plus complète n'est
 *  pas gênante.
 *----------------------------------------------------------------------------*/

static ecs_tab_int_t  ecs_loc_vec_def_rc__lis_sel_vis
(
 const ecs_vec_int_t   *const vec_fac_vis
) ;


/*----------------------------------------------------------------------------
 *  Fonction qui construit les listes des faces nouvelles et modifiées
 *----------------------------------------------------------------------------*/

static void ecs_loc_vec_def_rc__lis_fac_new
(
 const ecs_vec_int_t    *const vec_def_fac,      /*  -> Déf. faces            */
 const ecs_vec_int_t    *const vec_def_are,      /*  -> Déf. arêtes           */
 const ecs_vec_real_t   *const vec_def_som,      /*  -> Déf. sommets          */
 const ecs_vec_int_t    *const vec_fac_old_new,  /*  -> Renumérotation faces  */
 const ecs_tab_int_t           liste_som_new,    /*  -> Soms. nouv. ou modif. */
       ecs_tab_int_t    *const liste_fac_new,    /* <-  Nouvelles faces       */
       ecs_tab_int_t    *const liste_fac_mod     /* <-  Faces modifiées       */
) ;


/*----------------------------------------------------------------------------
 *  Fonction qui determine les distances maximales associées aux sommets
 *----------------------------------------------------------------------------*/

static ecs_tab_real_t ecs_loc_vec_def_rc__dist_max
(
 const ecs_vec_int_t  *const vec_def_are,
       ecs_vec_real_t *const vec_def_som,
 const ecs_real_t            fraction_dist
) ;


/*----------------------------------------------------------------------------
 *  Fonction qui détermine les points d'intersection lorsqu'on ne connaît
 *  que le tableau de sélection des faces à recoller
 *----------------------------------------------------------------------------*/

static ecs_rc_lis_inter_are_t ecs_loc_vec_def_rc__int_are_sel
(
 const ecs_vec_int_t   *const vec_def_fac,
 const ecs_vec_int_t   *const vec_def_are,
       ecs_vec_real_t  *const vec_def_som,
 const ecs_tab_int_t   *const tab_fac_select,
 const ecs_tab_real_t         dist_max_som,
       size_t                 nbr_are_ini,
 const ecs_real_t             fraction_dist
) ;


/*----------------------------------------------------------------------------
 *  Fonction qui détermine les points d'intersection lorsque la liste des
 *  des faces visibles par d'autres est déjà connue. Il n'est pas
 *  nécessaire que deux faces en vis-à-vis partiel se connaissent ; il
 *  suffit que l'une d'elles connaîsse l'autre (cas d'une filiation par
 *  exemple), même si une connaîssance plus complète n'est pas gênante.
 *----------------------------------------------------------------------------*/

static ecs_rc_lis_inter_are_t ecs_loc_vec_def_rc__int_are_vis
(
 const ecs_vec_int_t   *const vec_def_fac,
 const ecs_vec_int_t   *const vec_def_are,
       ecs_vec_real_t  *const vec_def_som,
 const ecs_vec_int_t   *const vec_fac_vis,
 const ecs_tab_real_t         dist_max_som,
 const ecs_real_t             fraction_dist
) ;


/*----------------------------------------------------------------------------
 *  Fonction qui détermine les points d'intersection lorsque les
 *  arêtes à intersecter partagent au départ déjà un sommet.
 *----------------------------------------------------------------------------*/

#if defined(__HP_cc)   /* On limite le niveau d'optimisation max. */
#pragma OPT_LEVEL 2
#endif

static ecs_rc_lis_inter_are_t ecs_loc_vec_def_rc__int_are_som
(
 const ecs_vec_int_t   *const vec_def_fac,
 const ecs_vec_int_t   *const vec_def_are,
       ecs_vec_real_t  *const vec_def_som,
 const ecs_tab_int_t   *const tab_fac_select,
 const ecs_tab_real_t         dist_max_som,
 const ecs_real_t             fraction_dist
) ;


/*----------------------------------------------------------------------------
 *  Fonction qui ajoute une intersection arête-arête à la liste
 *  d'intersections arêtes-arêtes ;
 *----------------------------------------------------------------------------*/

static void ecs_loc_vec_def_rc__aj_inter
(
 ecs_int_t                 are_0,
 ecs_int_t                 are_1,
 ecs_real_t                s_inter_0,
 ecs_real_t                t_inter_1,
 ecs_rc_lis_inter_are_t  * liste_inter_are
) ;


/*----------------------------------------------------------------------------
 *  Fonction qui détermine le ou les points d'intersection de 2 arêtes
 * dans l'espace 3D. On cherche en premier lieu à établir une correspondance
 * entre sommets, puis une correspondance entre sommets, et enfin une
 * ``vraie intersection,'' le tout sous une tolérance donnée.
 * On renvoie les coordonnées locales s et t d'une intersection. Si s ou t
 * = 0 ou 1, on met sa valeur à -0.01 ou 1.01 respectivement, afin de
 * pouvoir remplacer un test d'égalité à 0 ou 1 par un test d'inégalité,
 * sachant que 0 <= s, t <= 1 pour toute intersection.
 *----------------------------------------------------------------------------*/

static void ecs_loc_vec_def_rc__int_are_3d
(
 const ecs_point_t   som_0_are_0,        /* --> Coordonnées des sommets       */
 const ecs_point_t   som_1_are_0,
 const ecs_point_t   som_0_are_1,
 const ecs_point_t   som_1_are_1,
 const ecs_real_t    d_max_s0_a0,        /* --> Tolérances sur les distances  */
 const ecs_real_t    d_max_s1_a0,
 const ecs_real_t    d_max_s0_a1,
 const ecs_real_t    d_max_s1_a1,
 const ecs_real_t    fraction_dist,
       ecs_real_t    s_inter_0[2],       /* <-- Coord. s intersections seg. 0 */
       ecs_real_t    t_inter_1[2],       /* <-- Coord. t intersections seg. 1 */
       size_t      * nbr_inter           /* <-- Nombre d'intersections        */
) ;


/*----------------------------------------------------------------------------
 *  Fonction qui renvoie le sommet associé à une intersection
 *  (en le créant si nécessaire)
 *----------------------------------------------------------------------------*/


static ecs_int_t  ecs_loc_vec_def_rc__som_inter
(
       ecs_vec_int_t           *const vec_def_are,
       ecs_vec_real_t          *const vec_def_som,
       size_t                  *const nbr_som_new,
       ecs_tab_real_t          *const dist_max_som,
 const ecs_rc_inter_are_are_t  *const inter_are
) ;


/*----------------------------------------------------------------------------
 *    Fonction de descente d'un arbre binaire pour le tri lexicographique
 *  d'un tableau d'intersections d'arêtes.
 *----------------------------------------------------------------------------*/

static void ecs_loc_rc__desc_arbre_inter
(
       ecs_int_t                        ltree,       /* --> Niveau arbre      */
 const ecs_int_t                        ntree,       /* --> Taille arbre      */
       ecs_rc_inter_are_are_t    *const table_are,   /* <-> Tableau à trier   */
       ecs_rc_inter_are_inter_t  *const table_inter  /* <-> Tableau associé   */
) ;


/*----------------------------------------------------------------------------
 *    Fonction de tri lexicographique d'une liste d'intersections d'arêtes
 *----------------------------------------------------------------------------*/

static void ecs_loc_vec_def_rc__trie_inter
(
 ecs_rc_lis_inter_are_t *const liste_inter
) ;


#if 0 && defined(DEBUG) && !defined(NDEBUG)

/*----------------------------------------------------------------------------
 *  Fonction qui affiche une structure ecs_vec_int_t
 *----------------------------------------------------------------------------*/

static void  ecs_dbg_vec_def_rc__imp_vec_int
(
 const char           *const titre,
 const char           *const nom,
 const ecs_vec_int_t  *const vec_int_dbg
) ;

#endif


#if 0 && defined(DEBUG) && !defined(NDEBUG)

/*----------------------------------------------------------------------------
 *  Fonction qui imprime une structure ecs_vec_real_t
 *----------------------------------------------------------------------------*/

static void  ecs_dbg_vec_def_rc__imp_vec_real
(
 const char            *const titre,
 const char            *const nom,
 const ecs_vec_real_t  *const vec_real_dbg
 ) ;

#endif


/*============================================================================
 *                             Fonctions publiques
 *============================================================================*/


/*----------------------------------------------------------------------------
 *  Fonction qui recolle les faces non conformes
 *
 *  Les listes des faces nouvelles ou modifiées sont construites (et allouées)
 *  ici ; les structures liste_fac_new et liste_fac_mod correspondantes sont
 *  donc vides en entrée ; idem pour liste_fac_err qui indiquera les indices
 *  des faces pour lesquelles le découpage en sous-faces a échoué
 *
 *  On prend en entrée soit une définition de type "visibilité" entre faces
 *  à recoller (par exemple une filiation), ou alors une simple liste de
 *  faces sélectionnées. L'un de vec_fac_vis et tab_fac_select doit donc
 *  être à NULL, et l'autre non.
 *----------------------------------------------------------------------------*/

void ecs_vec_def__recolle
(
       ecs_vec_int_t   * *const vec_def_fac,      /* <-> Déf. des faces       */
       ecs_vec_int_t   * *const vec_def_are,      /* <-> Déf. des arêtes      */
       ecs_vec_real_t  * *const vec_def_som,      /* <-> Déf. des sommets     */
       ecs_vec_int_t   * *const vec_fac_old_new,  /* <-  Renum. faces         */
       ecs_vec_int_t     *const vec_fac_vis,      /*  -> Visibilité faces     */
       ecs_vec_int_t   * *const vec_fac_perio,    /*  -> Correspondance perio */
       ecs_vec_int_t   * *const vec_are_perio,    /*  -> Correspondance perio */
       ecs_vec_int_t   * *const vec_som_perio,    /*  -> Correspondance perio */
 const ecs_tab_int_t     *const tab_fac_select,   /*  -> Faces sélectionnées  */
       ecs_tab_int_t     *const liste_fac_new,    /* <-  Nouvelles faces      */
       ecs_tab_int_t     *const liste_fac_mod,    /* <-  Faces modifiées      */
       ecs_tab_int_t     *const liste_fac_err,    /* <-  Faces non découpées  */
 const ecs_param_rc_t           param_rc          /*  -> Param. tolérance     */
)
{

  ecs_vec_int_t  * vec_def_fac_new ;

  ecs_vec_int_t  * vec_som_old_new ;
  ecs_vec_int_t  * vec_are_old_new ;

  ecs_tab_int_t    liste_som_new ;

  ecs_tab_int_t    tab_fac_select_tmp ;

  float        * normale_fac = NULL ;

  const size_t   nbr_fac_ini = (*vec_def_fac)->pos_nbr - 1 ;


  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  ECS_CHRONO("Joining: start                      ") ;


  assert((*vec_def_som)->pos_pas == ECS_DIM_3) ;

  assert (   (vec_fac_vis != NULL && tab_fac_select == NULL)
          || (vec_fac_vis == NULL && tab_fac_select != NULL)) ;


  bft_printf
    ("\n\n  %s\n", _("Joining parameter values:")) ;

  bft_printf
    (_("  - Fraction of the length of the shortest edge\n"
       "    connected to a vertex under which 2 edges\n"
       "    intersect and 2 vertices are merged:                %.10E\n"),
     param_rc.fraction_dist) ;

  bft_printf
    (_("  - Value of the dot product of the unit normals of\n"
       "    2 faces above which these faces are considered\n"
       "    to be coplanar:                                     %.10E\n"),
     param_rc.epsilon_plan) ;


  if (param_rc.fraction_dist >= 0.5) {
    ecs_warn() ;
    bft_printf(_("It is highly recommended that the joining parameter\n"
                 "representing the fraction of the minimum length of "
                 "the edges\n"
                 "incident to a vertex be in the interval ]0, 0.5[\n"
                 "--> Only use the current value if you have a good "
                 "understanding\n"
                 "    of the joining algorithm and you are sure that the gap\n"
                 "    between the faces to join is everywhere positive."));
  }

  if (param_rc.semi_conforme == ECS_TRUE)
    bft_printf(_("  - \"semi-conforming\" search algorithm for joining\n"
                 "    of faces already sharing some vertices.\n")) ;


  /*------------------------------------------*/
  /* Calcul des normales associées aux faces  */
  /*------------------------------------------*/

  if (tab_fac_select != NULL) {
    tab_fac_select_tmp.nbr = tab_fac_select->nbr;
    tab_fac_select_tmp.val = tab_fac_select->val;
  }
  else
    tab_fac_select_tmp = ecs_loc_vec_def_rc__lis_sel_vis(vec_fac_vis) ;

  BFT_MALLOC(normale_fac, tab_fac_select_tmp.nbr*3, float);

  ecs_loc_vec_def_rc__normales(&tab_fac_select_tmp,
                               *vec_def_fac,
                               *vec_def_are,
                               *vec_def_som,
                               normale_fac) ;

  tab_fac_select_tmp.nbr = 0;
  if (tab_fac_select != NULL)
    tab_fac_select_tmp.val = NULL;
  else
    BFT_FREE(tab_fac_select_tmp.val) ;


  /*---------------------------------------------------------*/
  /* Calcul éventuel des matrices de rotation en périodicité */
  /*---------------------------------------------------------*/

  if (param_rc.param_perio != NULL)
    ecs_vec_def_perio__mat_rot_3d((*(param_rc.param_perio)).angle,
                                  (*(param_rc.param_perio)).direction,
                                  (*(param_rc.param_perio)).matrice) ;


  /*==========================================================================*/
  /* 1ERE PHASE                                                               */
  /*==========================================================================*/
  /*                                                                          */
  /* DETERMINATION DES NOUVEAUX SOMMETS D'INTERSECTION DE FACES               */
  /*                                                                          */
  /*==========================================================================*/


  ECS_CHRONO("Joining: initialization done        ") ;


  ecs_loc_vec_def_rc__def_are_new(*vec_def_fac,
                                  vec_def_are,
                                  *vec_def_som,
                                  &vec_are_old_new,
                                  &vec_som_old_new,
                                  vec_fac_vis,
                                  vec_fac_perio,
                                  vec_are_perio,
                                  vec_som_perio,
                                  tab_fac_select,
                                  &liste_som_new,
                                  param_rc) ;


  ECS_CHRONO("Joining: new edges done             ") ;



  /*==========================================================================*/
  /* 2NDE PHASE                                                               */
  /*==========================================================================*/
  /*                                                                          */
  /* CONSTRUCTION DES NOUVELLES FACES                                         */
  /*                                                                          */
  /*==========================================================================*/


  /* Duplication des normales aux faces sélectionnées dans le cas périodique */

  if (param_rc.param_perio != NULL) {

    assert(tab_fac_select->nbr == ecs_vec_int__ret_val_nbr(*vec_fac_perio)) ;

    ecs_vec_def_perio__duplique_norm(*vec_fac_perio,
                                     *(param_rc.param_perio),
                                     &normale_fac) ;

  }


  /* Construction des nouvelles faces */

  if (tab_fac_select != NULL)

    vec_def_fac_new = ecs_loc_vec_def_rc__def_fac_new(*vec_def_fac,
                                                      *vec_def_are,
                                                      *vec_def_som,
                                                      NULL,
                                                      normale_fac,
                                                      vec_fac_old_new,
                                                      tab_fac_select,
                                                      liste_fac_err,
                                                      nbr_fac_ini,
                                                      param_rc.epsilon_plan) ;


  else {

    tab_fac_select_tmp = ecs_loc_vec_def_rc__lis_sel_vis(vec_fac_vis) ;

    vec_def_fac_new = ecs_loc_vec_def_rc__def_fac_new(*vec_def_fac,
                                                      *vec_def_are,
                                                      *vec_def_som,
                                                      vec_fac_vis,
                                                      normale_fac,
                                                      vec_fac_old_new,
                                                      &tab_fac_select_tmp,
                                                      liste_fac_err,
                                                      nbr_fac_ini,
                                                      param_rc.epsilon_plan) ;

    tab_fac_select_tmp.nbr = 0 ;
    BFT_FREE(tab_fac_select_tmp.val) ;

  }

  ecs_vec_int__detruit(*vec_def_fac) ;

  BFT_FREE(normale_fac) ;


  *vec_def_fac = vec_def_fac_new ;


  ECS_CHRONO("Joining: new faces done             ") ;



  if (param_rc.param_perio != NULL) {

    /* On nettoie le maillage des éléments périodiques n'étant pas intervenu
       au cours du recollement */

    ecs_vec_def_perio__nettoie(vec_def_fac,
                               vec_def_are,
                               vec_def_som,
                               *vec_fac_old_new,
                               vec_are_old_new,
                               *vec_fac_perio,
                               *vec_are_perio,
                               *vec_som_perio,
                               liste_fac_err,
                               &liste_som_new,
                               *param_rc.param_perio) ;

    ECS_CHRONO("Periodicity: elements pre-cleaning done  ") ;


    /* Traitement spécifique pour une périodicité de rotation de 180 degrés.
       (On teste la nature de la périodicité dans la routine) */

    ecs_vec_def_perio__trait_spec_sym(vec_def_som,
                                      vec_fac_perio,
                                      vec_are_perio,
                                      vec_som_perio,
                                      liste_som_new,
                                      *param_rc.param_perio) ;


    /* Héritage du recollement dans le cas d'une périodicité */

    ecs_vec_def_perio__herite(vec_def_fac,
                              vec_def_are,
                              vec_def_som,
                              *vec_fac_old_new,
                              vec_are_old_new,
                              vec_som_old_new,
                              vec_fac_perio,
                              vec_are_perio,
                              vec_som_perio,
                              liste_fac_err,
                              &liste_som_new,
                              *param_rc.param_perio) ;

    ECS_CHRONO("Periodicity: elements inheritance done  ") ;


    /* Nettoyage de certains éléments ayant servi à l'héritage */

    ecs_vec_def_perio__nettoie(vec_def_fac,
                               vec_def_are,
                               vec_def_som,
                               *vec_fac_old_new,
                               vec_are_old_new,
                               *vec_fac_perio,
                               *vec_are_perio,
                               *vec_som_perio,
                               liste_fac_err,
                               &liste_som_new,
                               *param_rc.param_perio) ;

    ECS_CHRONO("Periodicity: elements cleaning done      ") ;


    /* Marquage d'éventuelles faces périodiques "se voyant" elles-mêmes */

    ecs_vec_def_perio__lis_fac_err(*vec_fac_perio,
                                   liste_fac_err) ;


  }


  /* Suppression de tableaux qui ne sont plus utiles */

  vec_som_old_new = ecs_vec_int__detruit(vec_som_old_new) ;
  vec_are_old_new = ecs_vec_int__detruit(vec_are_old_new) ;


  /* Liste des faces dont au moins un sommet est issu d'une intersection
     d'arêtes ou d'une fusion de sommets (i.e. liste des faces nouvelles
     ou modifiées par le recollement)                                     */

  ecs_loc_vec_def_rc__lis_fac_new(*vec_def_fac,
                                  *vec_def_are,
                                  *vec_def_som,
                                  *vec_fac_old_new,
                                  liste_som_new,
                                  liste_fac_new,
                                  liste_fac_mod) ;


  /* Libération d'un tableau alloué dans ecs_loc_vec_def_rc__def_are_new() */

  liste_som_new.nbr = 0 ;
  BFT_FREE(liste_som_new.val) ;


  if (param_rc.param_perio != NULL) {

    /* On redimensionne finalement vec_fac_old_new */

    ecs_vec_int__redimensionne(*vec_fac_old_new,
                               nbr_fac_ini + 1,
                               (*vec_fac_old_new)->pos_tab[nbr_fac_ini] - 1) ;

    /* A ce stade de développement, on n'a plus besoin de la correspondance
       entre les sommets périodiques, ainsi que celle des arêtes périodiques.
       Il suffit de conserver vec_fac_perio.                                 */

    ecs_vec_int__detruit(*vec_som_perio) ;
    ecs_vec_int__detruit(*vec_are_perio) ;

    /* Vérification de l'existence faces périodiques (on peut avoir eu
       des intersections arête-arête sans obtenir de faces périodiques,
       on ne s'est alors pas arrêté au premier test). */

    if (ecs_vec_def__compte_elt_perio(*vec_fac_perio) == 0)

      bft_error(__FILE__, __LINE__, 0,
                _("Problem for periodicity %1d.\nNo correspondance "
                  "was found between periodic faces."),
                param_rc.param_perio->num_perio) ;


  }


#if 0 && defined(DEBUG) && !defined(NDEBUG)

  {

    ecs_int_t     iare ;
    ecs_int_t     ifac ;
    size_t        nbr_are_fac ;
    size_t        pos_are_fac ;

    printf("\n\nrubrique 16 ecs_vec_def__recollement\n") ;
    printf("definition finale des fac new\n") ;

    printf("\n(*vec_def_fac) :\n") ;
    for (ifac = 0 ; ifac < (*vec_def_fac)->pos_nbr - 1 ; ifac++) {
      pos_are_fac = (*vec_def_fac)->pos_tab[ifac]     - 1 ;
      nbr_are_fac = (*vec_def_fac)->pos_tab[ifac + 1] - 1 - pos_are_fac ;
      printf("  fac new %d :", ifac + 1) ;
      for (iare = 0 ; iare < nbr_are_fac ; iare++) {
        printf(" %d", (*vec_def_fac)->val_tab[pos_are_fac + iare]) ;
      }
      printf("\n") ;
    }
    fflush(stdout);

  }

#endif


#if 0 && defined(DEBUG) && !defined(NDEBUG)

  {

    ecs_int_t     iare ;
    ecs_int_t     ifac ;
    size_t        nbr_are_fac ;
    size_t        pos_are_fac ;

    printf("\n\nrubrique 17 ecs_vec_def__recollement\n") ;
    printf("definition finale des fac old en fonction des fac new\n") ;

    printf("\n(*vec_fac_old_new) :\n") ;
    for (ifac = 0 ; ifac < (*vec_fac_old_new)->pos_nbr - 1 ; ifac++) {
      pos_are_fac = (*vec_fac_old_new)->pos_tab[ifac]     - 1 ;
      nbr_are_fac = (*vec_fac_old_new)->pos_tab[ifac + 1] - 1 - pos_are_fac ;
      printf("  fac old %d :", ifac + 1) ;
      for (iare = 0 ; iare < nbr_are_fac ; iare++) {
        printf(" %d", (*vec_fac_old_new)->val_tab[pos_are_fac + iare]) ;
      }
      printf("\n") ;
    }
    fflush(stdout);

  }

#endif

}


/*============================================================================
 *                              Fonctions privées
 *============================================================================*/

/*----------------------------------------------------------------------------
 *  Fonction qui renvoie l'ordonnancement des éléments
 *   en fonction de leur définition (sans modifier cette définition)
 *----------------------------------------------------------------------------*/

static ecs_tab_int_t ecs_loc_vec_def_rc__trie_elt
(
 const ecs_vec_real_t  *const this_vec_def_real,
 const ecs_real_t             tolerance
)
{

  ecs_tab_int_t    vect_transf ;  /* Vecteur de transformation */

  ecs_int_t      nbr_elt ;        /* Nombre d'éléments */
  ecs_int_t      ielt ;           /* Indice de boucle sur les éléments */


  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  /*--------------------------------------------------------------------------*/
  /* Tri des éléments en fonction de leurs définitions                        */
  /*--------------------------------------------------------------------------*/

  nbr_elt = ecs_vec_real__ret_pos_nbr(this_vec_def_real) - 1 ;

  vect_transf.nbr = nbr_elt ;
  BFT_MALLOC(vect_transf.val, nbr_elt, ecs_int_t) ;

  for (ielt = 0 ; ielt < nbr_elt ; ielt++)
    vect_transf.val[ielt] = ielt ;

  ecs_vec_real__trie(this_vec_def_real,
                     &vect_transf,
                     tolerance) ;

  /*--------------------------------------------------------------------------*/
  /* Renvoi du vecteur de transformation                                      */
  /*--------------------------------------------------------------------------*/

  return vect_transf ;

}


/*----------------------------------------------------------------------------
 * Fonction qui calcule les normales unitaires aux faces
 *----------------------------------------------------------------------------*/

void ecs_loc_vec_def_rc__normales
(
 const ecs_tab_int_t  *const tab_fac_select,
 const ecs_vec_int_t  *const vec_fac_def,
 const ecs_vec_int_t  *const vec_are_def,
 const ecs_vec_real_t *const vec_som_def,
       float          *const normale_fac
)
{

  size_t     ifac_select ;

  ecs_int_t  icoo ;
  ecs_int_t  isom ;
  ecs_int_t  iare ;
  ecs_int_t  ifac ;
  ecs_int_t  itri ;

  ecs_int_t  icoo_inf ;
  ecs_int_t  iare_inf ;
  ecs_int_t  iare_sup ;

  ecs_int_t  nbr_max_som_fac ;
  ecs_int_t  nbr_som_fac ;
  ecs_int_t  ind_som ;
  ecs_int_t  ind_are ;

  ecs_point_t  vect1 ;
  ecs_point_t  vect2 ;

  ecs_point_t *coo_som_fac ;
  ecs_point_t  barycentre_fac ;
  ecs_point_t  normale_loc ;
  ecs_point_t  normale_tri ;

  double       inv_surface ;


  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  assert(vec_som_def->pos_pas == ECS_DIM_3) ;


  nbr_max_som_fac = 0 ;

  for (ifac_select = 0 ; ifac_select < tab_fac_select->nbr ; ifac_select++) {

    ifac = tab_fac_select->val[ifac_select] ;

    nbr_som_fac = vec_fac_def->pos_tab[ifac + 1] - vec_fac_def->pos_tab[ifac] ;

    if (nbr_max_som_fac <= nbr_som_fac) nbr_max_som_fac = nbr_som_fac ;

  }

  BFT_MALLOC(coo_som_fac, nbr_max_som_fac + 1, ecs_point_t) ;


  /*---------------------------------------------------------------------------
   *
   *                          Pi+1
   *              *---------*                   B  : barycentre du polygone P
   *             / .       .                                            \
   *            /   .     .   \                 Pi : sommet     du polygone P
   *           /     .   .                                      \
   *          /       . .  Ti   \               Ti : triangle
   *         *.........B.........* Pi
   *     Pn-1 \       . .       /
   *           \     .   .     /
   *            \   .     .   /
   *             \ .   T0  . /
   *              *---------*
   *            P0
   */


  /*=========================================================================*/
  /* Boucle sur les faces sélectionnées                                      */
  /*=========================================================================*/


  for (ifac_select = 0 ; ifac_select < tab_fac_select->nbr ; ifac_select++) {

    ifac = tab_fac_select->val[ifac_select] ;


    /*------------------------------------------------------------------------*/
    /*                                                                        */
    /* Definition du polygone P en fonction des sommets Pi de la face         */
    /*                                                                        */
    /*------------------------------------------------------------------------*/


    iare_inf = vec_fac_def->pos_tab[ifac]     - 1 ;
    iare_sup = vec_fac_def->pos_tab[ifac + 1] - 1 ;

    nbr_som_fac = 0 ;


    for (iare = iare_inf ; iare < iare_sup ; iare++) {

      ind_are  = vec_fac_def->val_tab[iare] ;

      if (ind_are > 0)
        isom = vec_are_def->pos_tab[ind_are          - 1]  - 1 ;
      else
        isom = vec_are_def->pos_tab[ECS_ABS(ind_are) - 1] ;

      ind_som  = vec_are_def->val_tab[isom] - 1 ;

      icoo_inf = vec_som_def->pos_pas * ind_som ;


      for (icoo = X ; icoo < ECS_DIM_3 ; icoo++) {

        coo_som_fac[nbr_som_fac][icoo] = vec_som_def->val_tab[icoo_inf + icoo] ;


      } /* Fin : boucle sur les coordonnees du sommet */


      nbr_som_fac++ ;



    } /* Fin : boucle sur les aretes de la face */


    for (icoo = X ; icoo < ECS_DIM_3 ; icoo++)
      coo_som_fac[nbr_som_fac][icoo] = coo_som_fac[0][icoo] ;



    /*------------------------------------------------------------------------*/
    /*                                                                        */
    /* Calcul des coordonnees du barycentre B du polygone P                   */
    /*                                                                        */
    /*                                                                        */
    /*  -->    1   n-1  -->                                                   */
    /*  OB  =  -  Somme OPi                                                   */
    /*         n   i=0                                                        */
    /*                                                                        */
    /*                                                                        */
    /*------------------------------------------------------------------------*/


    for (icoo = X ; icoo < ECS_DIM_3 ; icoo++) {

      barycentre_fac[icoo] = 0. ;

      for (isom = 0 ; isom < nbr_som_fac ; isom++)
        barycentre_fac[icoo] += coo_som_fac[isom][icoo] ;

      barycentre_fac[icoo] /= nbr_som_fac ;

    }


    for (icoo = X ; icoo < ECS_DIM_3 ; icoo++) {
      normale_loc[icoo] = 0. ;
    }


    /* Premiere boucle sur les triangles de la face (calcul des normales)     */
    /*========================================================================*/


    for (itri = 0 ; itri < nbr_som_fac ; itri++) {


      /*----------------------------------------------------------------------*/
      /*                                                                      */
      /* Calcul de la normale au triangle Ti :                                */
      /*                                                                      */
      /*  ->            -->   -->                                             */
      /*  N(Ti) = 1/2 ( BPi X BPi+1 )                                         */
      /*                                                                      */
      /*----------------------------------------------------------------------*/


      for (icoo = X ; icoo < ECS_DIM_3 ; icoo++) {

        vect1[icoo] = coo_som_fac[itri    ][icoo] - barycentre_fac[icoo] ;
        vect2[icoo] = coo_som_fac[itri + 1][icoo] - barycentre_fac[icoo] ;

      }

      ECS_LOC_PRODUIT_VECTORIEL(normale_tri, vect1, vect2) ;

      for (icoo = X ; icoo < ECS_DIM_3 ; icoo++)
        normale_tri[icoo] /= 2. ;


      /*----------------------------------------------------------------------*/
      /*                                                                      */
      /* Calcul de la normale au polygone                                     */
      /*  comme somme vectorielle des normales aux triangles :                */
      /*                                                                      */
      /*  ->      n-1   ->                                                    */
      /*  N(P) = Somme( N(Ti) )                                               */
      /*          i=0                                                         */
      /*                                                                      */
      /*----------------------------------------------------------------------*/


      for (icoo = X ; icoo < ECS_DIM_3 ; icoo++)
        normale_loc[icoo] += normale_tri[icoo] ;

    } /* Fin de la premiere boucle sur les triangles de la face */


    inv_surface = 1.0 / (sqrt(  normale_loc[0] * normale_loc[0]
                              + normale_loc[1] * normale_loc[1]
                              + normale_loc[2] * normale_loc[2])) ;


    /* Stockage du resultat dans la structure appropriee */
    /*---------------------------------------------------*/

    for (icoo = X ; icoo < ECS_DIM_3 ; icoo++) {

      normale_fac[ifac_select * ECS_DIM_3 + icoo]
        = normale_loc[icoo] * inv_surface ;

    }

  } /* Fin : boucle sur les faces */


  BFT_FREE(coo_som_fac) ;


}


/*----------------------------------------------------------------------------
 *  Construction des nouvelles arêtes decoupées par les intersections
 *
 *  La liste des sommets finaux issus d'une intersection ou d'une fusion
 *  est construite (et allouée) ici ; la structure liste_som_new
 *  correspondante est donc vide en entrée ;
 *----------------------------------------------------------------------------*/

static void ecs_loc_vec_def_rc__def_are_new
(
       ecs_vec_int_t         *const vec_def_fac,
       ecs_vec_int_t       * *const vec_def_are,
       ecs_vec_real_t        *const vec_def_som,
       ecs_vec_int_t       * *const vec_are_old_new,
       ecs_vec_int_t       * *const vec_som_old_new,
       ecs_vec_int_t         *const vec_fac_vis,
       ecs_vec_int_t       * *const vec_fac_perio,
       ecs_vec_int_t       * *const vec_are_perio,
       ecs_vec_int_t       * *const vec_som_perio,
 const ecs_tab_int_t         *const tab_fac_select,
       ecs_tab_int_t         *const liste_som_new,
 const ecs_param_rc_t               param_rc
)
{

  ecs_tab_int_t        tab_renum_are ;
  ecs_tab_int_t      * tab_signe_are ;

  ecs_rc_lis_inter_are_t   liste_inter_are ;

  ecs_tab_int_t        tab_equiv_som ;
  ecs_tab_real_t       dist_max_som ;

  const size_t         nbr_som_ini   = vec_def_som->pos_nbr - 1 ;
  const size_t         nbr_are_ini   = (*vec_def_are)->pos_nbr - 1 ;
  const ecs_real_t     fraction_dist = param_rc.fraction_dist ;


  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  dist_max_som = ecs_loc_vec_def_rc__dist_max(*vec_def_are,
                                              vec_def_som,
                                              fraction_dist) ;


  /* Duplication des éléments du maillage sélectionnés dans le cas périodique */

  if (param_rc.param_perio != NULL)

    ecs_vec_def_perio__duplique_def(vec_def_fac,
                                    *vec_def_are,
                                    vec_def_som,
                                    vec_fac_perio,
                                    vec_are_perio,
                                    vec_som_perio,
                                    tab_fac_select,
                                    &dist_max_som,
                                    *(param_rc.param_perio)) ;


  /* Recherche des intersections entre arêtes */

  if (vec_fac_vis == NULL && tab_fac_select != NULL) {

    if (param_rc.semi_conforme == ECS_FALSE)

      liste_inter_are = ecs_loc_vec_def_rc__int_are_sel(vec_def_fac,
                                                        *vec_def_are,
                                                        vec_def_som,
                                                        tab_fac_select,
                                                        dist_max_som,
                                                        nbr_are_ini,
                                                        fraction_dist) ;

    else

      liste_inter_are = ecs_loc_vec_def_rc__int_are_som(vec_def_fac,
                                                        *vec_def_are,
                                                        vec_def_som,
                                                        tab_fac_select,
                                                        dist_max_som,
                                                        fraction_dist) ;
  }
  else if (vec_fac_vis != NULL && tab_fac_select == NULL)

    liste_inter_are = ecs_loc_vec_def_rc__int_are_vis(vec_def_fac,
                                                      *vec_def_are,
                                                      vec_def_som,
                                                      vec_fac_vis,
                                                      dist_max_som,
                                                      fraction_dist) ;

  else
    assert(   (vec_fac_vis != NULL && tab_fac_select == NULL)
           || (vec_fac_vis == NULL && tab_fac_select != NULL)) ;


  /* Vérification de l'existence possible de recollement périodique */

  if (param_rc.param_perio != NULL && liste_inter_are.nbr == 0)

    bft_error(__FILE__, __LINE__, 0,
              _("Problem for periodicity %1d.\nNo correspondance "
                "was found between periodic faces."),
              param_rc.param_perio->num_perio) ;


  /* Tri du tableau des intersections */

  ecs_loc_vec_def_rc__trie_inter(&liste_inter_are) ;


  /* Découpages des arêtes */

  tab_equiv_som.nbr = 0 ;
  tab_equiv_som.val = NULL ;

  *vec_are_old_new = ecs_loc_vec_def_rc__decoup_are(*vec_def_are,
                                                    vec_def_som,
                                                    liste_inter_are,
                                                    &tab_equiv_som,
                                                    &dist_max_som) ;


  /* On en a terminé avec l'utilisation de la liste des intersections */

  BFT_FREE(liste_inter_are.are) ;
  BFT_FREE(liste_inter_are.inter) ;

  liste_inter_are.nbr     = 0 ;
  liste_inter_are.nbr_max = 0 ;
  liste_inter_are.are     = NULL ;
  liste_inter_are.inter   = NULL ;


  /* Fusion des sommets */

  *vec_som_old_new = ecs_vec_def__fusion_som(*vec_def_are,
                                             vec_def_som,
                                             &tab_equiv_som,
                                             &dist_max_som,
                                             liste_som_new) ;


#if 0 && defined(DEBUG) && !defined(NDEBUG)

  {

    printf("\n\ecs_loc_vec_def_rc__def_are_new\n") ;
    printf("nouvelles définitions des sommets\n") ;

    for (ind_som = 0 ;
         ind_som < vec_def_som->pos_nbr - 1 ;
         ind_som++)

      printf("%10d : %12.5e %12.5e %12.5e\n", ind_som + 1,
             vec_def_som->val_tab[ind_som*3],
             vec_def_som->val_tab[ind_som*3 + 1],
             vec_def_som->val_tab[ind_som*3 + 2]);

  }

#endif

#if 0 && defined(DEBUG) && !defined(NDEBUG)

  ecs_dbg_vec_def_rc__imp_vec_int("ecs_loc_vec_def_rc__def_are_new",
                                  "vec_som_old_new",
                                  vec_som_old_new) ;

#endif

#if 0 && defined(DEBUG) && !defined(NDEBUG)

  {
    ecs_int_t ind_are ;

    printf("\n\ecs_loc_vec_def_rc__def_are_new\n") ;
    printf("nouvelles définitions des arêtes en fonction des sommets\n") ;

    for (ind_are = 0 ;
         ind_are < vec_def_are->pos_nbr - 1 ;
         ind_are++)

      printf("%10d : %10d %10d\n", ind_are + 1,
             vec_def_are->val_tab[ind_are*2],
             vec_def_are->val_tab[ind_are*2 + 1]);

  }

#endif


  /* Libération de tableaux auxiliaires */

  tab_equiv_som.nbr = 0 ;
  BFT_FREE(tab_equiv_som.val) ;

  dist_max_som.nbr = 0 ;
  BFT_FREE(dist_max_som.val) ;


  /* Réactualisation du vecteur de correspondance des sommets
     dans le cas périodique */

  if (param_rc.param_perio != NULL)

    ecs_vec_def__maj_elt_perio(*vec_som_perio,
                               *vec_som_old_new,
                               vec_def_som->pos_nbr - 1) ;


  /*
    Pour la suite, on n'aura besoin dans vec_som_old_new que des relations
    entre les sommets initiaux (avant intersections) et les sommets
    fusionnés. Les sommets issus d'intersection ayant été rajoutés après
    la position "nbr_som_ini", il suffit de redimensionner à la baisse
    vec_som_old_new pour supprimer ces références.
  */

  ecs_vec_int__redimensionne(*vec_som_old_new,
                             nbr_som_ini + 1,
                             nbr_som_ini) ;


  /*--------------------------------------------------------*/
  /* Reactualisation des listes de sommets, arêtes et faces */
  /*--------------------------------------------------------*/


  bft_printf(_("\n  Merging of topologically identical edges:\n")) ;

  bft_printf(_("    Initial number of edges           : %10d\n"),
             (*vec_def_are)->pos_nbr - 1) ;

  /* Compactage des arêtes topologiquement identiques */
  /*--------------------------------------------------*/

  BFT_MALLOC(tab_signe_are, 1, ecs_tab_int_t) ;

  tab_renum_are = ecs_vec_def__fusionne(vec_def_are,
                                        tab_signe_are) ;


#if 0 && defined(DEBUG) && !defined(NDEBUG)

  {

    ecs_int_t          iare ;


    printf("\n\necs_loc_vec_def_rc__def_are_new\n") ;
    printf("définition des arêtes apres fusion\n") ;

    printf("\nvec_def_are :\n") ;
    for (iare = 0 ; iare < (*vec_def_are)->pos_nbr - 1 ; iare++)
      printf("  are %d : %d %d  (signe : %2d)\n", iare + 1,
             (*vec_def_are)->val_tab[2 * iare],
             (*vec_def_are)->val_tab[2 * iare + 1],
             tab_signe_are->val[iare]) ;
  }

#endif

  bft_printf(_("    Number of edges after merging     : %10d\n"),
             (*vec_def_are)->pos_nbr - 1) ;


  /* Réactualisation de la définition des anciennes arêtes */
  /*  en fonction des nouvelles                            */
  /*-------------------------------------------------------*/

  ecs_vec_int__renumerote_et_sgn(*vec_are_old_new,
                                 tab_renum_are,
                                 *tab_signe_are) ;

#if 0 && defined(DEBUG) && !defined(NDEBUG)
  {
    ecs_dbg_vec_def_rc__imp_vec_int
      ("Réactualisation de la définition des anciennes arêtes",
       "vec_are_old_new",
       *vec_are_old_new) ;
  }
#endif

  BFT_FREE(tab_renum_are.val) ;

  BFT_FREE(tab_signe_are->val) ;
  BFT_FREE(tab_signe_are) ;


  /* Suppression des arêtes dégénérées éventuelles */
  /*-----------------------------------------------*/

  ecs_loc_vec_def_rc__nettoie_are(*vec_def_are,
                                  *vec_are_old_new) ;


#if 0 && defined(DEBUG) && !defined(NDEBUG)

  ecs_dbg_vec_def_rc__imp_vec_int("ecs_loc_vec_def_rc__def_are_new",
                                  "vec_def_are (définition finale)",
                                  *vec_def_are) ;

#endif


  /* Remplacement dans la definition des faces,                           */
  /*  des anciennes arêtes par les (sous-)arêtes issues de leur decoupage */
  /*----------------------------------------------------------------------*/

  ecs_vec_int__remplace_ref(vec_def_fac,
                            *vec_are_old_new) ;


#if 0 && defined(DEBUG) && !defined(NDEBUG)

  ecs_dbg_vec_def_rc__imp_vec_int("ecs_loc_vec_def_rc__def_are_new",
                                  "vec_def_fac (avec nouvelles arêtes)",
                                  vec_def_fac) ;

#endif

  /* Suppression des coins dégénérés éventuels */
  /*-------------------------------------------*/

  ecs_loc_vec_def_rc__nettoie_fac(vec_def_fac);

}


/*----------------------------------------------------------------------------
 *  Fonction qui construit les nouvelles faces
 *----------------------------------------------------------------------------*/

static ecs_vec_int_t * ecs_loc_vec_def_rc__def_fac_new
(
 const ecs_vec_int_t    *const vec_def_fac,
 const ecs_vec_int_t    *const vec_def_are,
 const ecs_vec_real_t   *const vec_def_som,
 const ecs_vec_int_t    *const vec_fac_vis,
 const float            *const normale_fac,
       ecs_vec_int_t  * *const vec_fac_old_new,
 const ecs_tab_int_t    *const tab_fac_select,
       ecs_tab_int_t    *const liste_fac_err,
       size_t                  nbr_fac_ini,
 const ecs_real_t              epsilon_plan
)
{

  size_t         cpt_fac_select ;
  size_t         cpt_fac_new ;
  ecs_int_t      cpt_val_fac_new ;

  size_t         iare ;
  size_t         ifac ;
  size_t         ifac_select ;
  size_t         ipos ;
  size_t         ipos_deb ;
  size_t         isom ;
  size_t         ival ;

  size_t         max_pos_fac_new ;
  size_t         max_val_fac_new ;

  size_t         nbr_are ;
  size_t         nbr_are_fac ;
  size_t         nbr_fac ;
  size_t         nbr_fac_pb ;
  size_t         nbr_fac_pb_are_aucune ;
  size_t         nbr_fac_pb_are_double ;
  size_t         nbr_fac_pb_are_boucle ;
  size_t         nbr_max_fac_dec ;
  size_t         nbr_retour_pb_fac ;
  size_t         nbr_som ;
  size_t         nbr_val_fac ;

  ecs_int_t      num_are ;

  size_t         pos_are_fac ;

  ecs_int_t      retour ;

  ecs_tab_int_t  liste_are_tete ;
  ecs_tab_int_t  liste_are_ext ;
  ecs_tab_int_t  liste_are_fac ;
  ecs_tab_int_t  liste_are_int ;

  ecs_tab_int_t    masque_fac_select ;
  ecs_tab_int_t    tab_renum_fac ;

  ecs_tab_int_t  * tab_signe_fac ;

  ecs_vec_int_t  * vec_def_fac_new ;
  ecs_vec_int_t  * vec_som_are ;
  ecs_vec_int_t  * vec_are_fac ;

#if 0 && defined(DEBUG) && !defined(NDEBUG)
  ecs_int_t cpt_fac_dec = 0 ;
#endif


  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  /* Critère d'arrêt */

  nbr_max_fac_dec = 100 ;
  if (getenv("ECS_RC_MAX_FAC_DEC") != NULL) {
    nbr_max_fac_dec = atoi(getenv("ECS_RC_MAX_FAC_DEC")) ;
    nbr_max_fac_dec = ECS_MAX(1, nbr_max_fac_dec) ;
  }


  /* Initialisations */

  nbr_fac = vec_def_fac->pos_nbr - 1 ;
  nbr_are = vec_def_are->pos_nbr - 1 ;
  nbr_som = vec_def_som->pos_nbr - 1 ;
  nbr_val_fac = ecs_vec_int__ret_val_nbr(vec_def_fac) ;

  liste_are_tete.nbr = 4 ;
  liste_are_ext.nbr  = 4 ;
  liste_are_int.nbr  = 4 ;
  liste_are_fac.nbr  = 4 ;

  BFT_MALLOC(liste_are_tete.val, liste_are_tete.nbr, ecs_int_t) ;
  BFT_MALLOC(liste_are_ext.val,  liste_are_ext.nbr,  ecs_int_t) ;
  BFT_MALLOC(liste_are_int.val,  liste_are_int.nbr,  ecs_int_t) ;
  BFT_MALLOC(liste_are_fac.val,  liste_are_fac.nbr,  ecs_int_t) ;

  nbr_fac_pb = 0 ;
  nbr_fac_pb_are_aucune = 0 ;
  nbr_fac_pb_are_double = 0 ;
  nbr_fac_pb_are_boucle = 0 ;


  /* Construction de la nouvelle connectivité inverse "arête -> face" */

  vec_are_fac = ecs_vec_int__inverse(vec_def_fac,
                                     vec_def_are->pos_nbr - 1) ;


  vec_def_fac_new = ecs_vec_int__alloue(nbr_fac + 1,
                                        nbr_val_fac) ;

  vec_def_fac_new->pos_tab[0] = 1 ;

  max_pos_fac_new = nbr_fac + 1 ;
  max_val_fac_new = nbr_val_fac ;
  cpt_fac_select  = 0 ;
  cpt_fac_new     = 0 ;
  cpt_val_fac_new = 0 ;


  *vec_fac_old_new = ecs_vec_int__alloue(nbr_fac + 1,
                                         nbr_val_fac) ;

  (*vec_fac_old_new)->pos_tab[0] = 1 ;

  masque_fac_select.nbr = nbr_fac ;
  BFT_MALLOC(masque_fac_select.val, masque_fac_select.nbr, ecs_int_t) ;


  for (ifac = 0 ; ifac < nbr_fac ; ifac++)
    masque_fac_select.val[ifac] = -1 ;

  for (ifac_select = 0 ; ifac_select < tab_fac_select->nbr ; ifac_select++)
    masque_fac_select.val[tab_fac_select->val[ifac_select]] = cpt_fac_select++ ;


  /* On rajoute les faces dupliquées dans le cas périodique */

  for (ifac = nbr_fac_ini ; ifac < nbr_fac ; ifac++)
    masque_fac_select.val[ifac] = cpt_fac_select++ ;


  /*
    On compacte la connectivité arête -> faces de manière à ce qu'elle
    ne contienne plus que les références aux faces sélectionnées
    (-> gain CPU et mémoire)
  */

  for (ival = 0 ; ival < vec_are_fac->pos_tab[nbr_are] - 1 ; ival++)
    if (masque_fac_select.val[ECS_ABS(vec_are_fac->val_tab[ival]) - 1] < 0)
      vec_are_fac->val_tab[ival] = 0 ;

  ival = 0 ;
  for (iare = 0 ; iare < nbr_are ; iare++) {
    ipos_deb = ival ;
    for (ipos = vec_are_fac->pos_tab[iare    ] - 1 ;
         ipos < vec_are_fac->pos_tab[iare + 1] - 1 ;
         ipos++) {
      if (vec_are_fac->val_tab[ipos] != 0)
        vec_are_fac->val_tab[ival++] = vec_are_fac->val_tab[ipos] ;
    }
    vec_are_fac->pos_tab[iare] = ipos_deb + 1 ;
  }
  vec_are_fac->pos_tab[nbr_are] = ival + 1 ;

  ecs_vec_int__redimensionne(vec_are_fac, nbr_are + 1, ival) ;


  /* Construction de la nouvelle connectivité inverse "sommet -> arêtes" */

  vec_som_are = ecs_vec_int__inverse(vec_def_are,
                                     nbr_som) ;


  /*
    On compacte la connectivité sommet -> arêtes de manière à ce qu'elle
    ne contienne plus que les références aux arêtes des faces sélectionnées
    (-> gain CPU et mémoire)
  */

  for (ival = 0 ; ival < vec_som_are->pos_tab[nbr_som] - 1 ; ival++) {
    iare = vec_som_are->val_tab[ival] - 1 ;
    if (vec_are_fac->pos_tab[iare + 1] == vec_are_fac->pos_tab[iare])
      vec_som_are->val_tab[ival] = 0 ;
  }

  ival = 0 ;
  for (isom = 0 ; isom < nbr_som ; isom++) {
    ipos_deb = ival ;
    for (ipos = vec_som_are->pos_tab[isom    ] - 1 ;
         ipos < vec_som_are->pos_tab[isom + 1] - 1 ;
         ipos++) {
      if (vec_som_are->val_tab[ipos] != 0)
        vec_som_are->val_tab[ival++] = vec_som_are->val_tab[ipos] ;
    }
    vec_som_are->pos_tab[isom] = ipos_deb + 1 ;
  }
  vec_som_are->pos_tab[nbr_som] = ival + 1 ;

  ecs_vec_int__redimensionne(vec_som_are, nbr_som + 1, ival) ;


#if 0 && defined(DEBUG) && !defined(NDEBUG)

  {
    size_t   pos_som ;
    size_t   pos_som_sup ;

    printf("\n\nrubrique 14 ecs_loc_vec_def_rc__def_fac_new\n") ;
    printf("definition des som en fonction des are\n") ;


    printf("\nvec_som_are :\n") ;
    for (isom = 0 ; isom < vec_som_are->pos_nbr - 1 ; isom++) {
      pos_som = vec_som_are->pos_tab[isom]     - 1 ;
      pos_som_sup = vec_som_are->pos_tab[isom + 1] - 1 ;
      printf("  som %d :", isom + 1) ;
      for (iare = 0 ; iare < pos_som_sup - pos_som ; iare++) {
        num_are = vec_som_are->val_tab[pos_som + iare] ;
        printf(" %d", num_are) ;
      }
      printf("\n") ;
    }
    fflush(stdout);

  }

#endif


  /*======================*/
  /* Boucle sur les faces */
  /*======================*/


#if 0 && defined(DEBUG) && !defined(NDEBUG)
  printf("\n\nrubrique 15 ecs_loc_vec_def_rc__def_fac_new\n") ;
  printf("pour chaque fac old : definition des fac new\n") ;
#endif

  ifac = 0 ;

  while (ifac < nbr_fac) {

    pos_are_fac = vec_def_fac->pos_tab[ifac] - 1 ;

    nbr_are_fac
      = vec_def_fac->pos_tab[ifac + 1] - 1
      - pos_are_fac ;

    /*
      Si la face n'est pas sélectionnée ou que l'on n'est pas parvenu à la
      reconstruire, on la recopie simplement
    */

    if (masque_fac_select.val[ifac] < 0) {


      cpt_fac_new     = (*vec_fac_old_new)->pos_tab[ifac] - 1 ;
      cpt_val_fac_new = vec_def_fac_new->pos_tab[cpt_fac_new] - 1 ;


      if (max_val_fac_new <= cpt_val_fac_new + nbr_are_fac - 1) {
        max_val_fac_new = ECS_MAX(max_val_fac_new * 2,
                                  cpt_val_fac_new + nbr_are_fac) ;
        BFT_REALLOC(vec_def_fac_new->val_tab, max_val_fac_new, ecs_int_t) ;
      }

      for (iare = 0 ; iare < nbr_are_fac ; iare++) {
        num_are = vec_def_fac->val_tab[pos_are_fac + iare] ;
        vec_def_fac_new->val_tab[cpt_val_fac_new++] = num_are ;
      }

      if (max_pos_fac_new <= cpt_fac_new + 1) {
        max_pos_fac_new = ECS_MAX(max_pos_fac_new * 2, cpt_fac_new + 1) ;
        BFT_REALLOC(vec_def_fac_new->pos_tab, max_pos_fac_new, ecs_size_t) ;
        BFT_REALLOC((*vec_fac_old_new)->val_tab, max_pos_fac_new - 1,
                    ecs_int_t) ;
      }

      vec_def_fac_new->pos_tab[cpt_fac_new + 1]
        = vec_def_fac_new->pos_tab[cpt_fac_new] + nbr_are_fac ;

      (*vec_fac_old_new)->val_tab[cpt_fac_new] = cpt_fac_new + 1 ;

      (*vec_fac_old_new)->pos_tab[ifac + 1]
        = (*vec_fac_old_new)->pos_tab[ifac] + 1 ;

      cpt_fac_new += 1 ;


      ifac++ ;

      nbr_retour_pb_fac = 0 ;

      continue ;

    }

    /* Si la face est sélectionnée, on la traite */


    retour =  ecs_loc_vec_def_rc__decoup_fac(vec_def_fac,
                                             vec_def_are,
                                             vec_def_som,
                                             normale_fac,
                                             masque_fac_select.val,
                                             vec_fac_vis,
                                             vec_som_are,
                                             vec_are_fac,
                                             vec_def_fac_new,
                                             *vec_fac_old_new,
                                             ifac,
                                             epsilon_plan,
                                             nbr_max_fac_dec,
                                             &liste_are_tete,
                                             &liste_are_ext,
                                             &liste_are_int,
                                             &liste_are_fac,
                                             &max_pos_fac_new,
                                             &max_val_fac_new) ;

    if (retour != 0) {

      nbr_retour_pb_fac += 1 ;

      /*
        On rééssaie en permutant l'ordre des arêtes dans la définition
        de la face de manière à démarrer d'une arête différente
        (ce qui peut suffire dans quelques cas de faces déformées)
        on utilise liste_are_ext comme stockage temporaire, car il a
        été bien dimensionné lors des appels précédents
      */

      for (iare = 0 ; iare < nbr_are_fac ; iare++)
        liste_are_ext.val[iare] = vec_def_fac->val_tab[pos_are_fac + iare] ;

      vec_def_fac->val_tab[pos_are_fac] = liste_are_ext.val[nbr_are_fac - 1] ;
      for (iare = 1 ; iare < nbr_are_fac ; iare++)
        vec_def_fac->val_tab[pos_are_fac + iare]
          = liste_are_ext.val[iare - 1] ;

      /* Si on a épuisé les essais, on abandonne pour cette face */

      if (nbr_retour_pb_fac >= nbr_are_fac) {

        masque_fac_select.val[ifac] = -retour -1 ;
        nbr_fac_pb += 1 ;

        switch (retour) {
        case 1:
          nbr_fac_pb_are_aucune += 1 ;
          break ;
        case 2:
          nbr_fac_pb_are_double += 1 ;
          break ;
        case 3:
          nbr_fac_pb_are_boucle += 1 ;
          break ;
        }

      }

      /*
        On n'incrémente pas ifac ici afin de repasser dans le cas
        "masque_fac_select.val[ifac] < 0" de la boucle sur les faces
        ou dans le découpage de la face, selon le nombre de tentatives
        échouées
      */

    }
    else {

      ifac++ ;

      nbr_retour_pb_fac = 0 ;

    }


  }  /* Fin : boucle sur les faces */


  ecs_vec_int__detruit(vec_are_fac) ;


  /* Information utilisateur sur les faces qui n'ont pas pu être recollées */

  if (nbr_fac_pb > 0) {

    ecs_warn() ;
    bft_printf(_("Reconstruction problem for %lu face(s) during joining\n"
                 "(%lu open cycles, %lu edges traversed twice, %lu faces split\n"
                 "into more than ECS_RC_MAX_FAC_DEC = %u).\n\n"
                 "-> Eventually modify joining parameters"),
               (unsigned long)nbr_fac_pb, (unsigned long)nbr_fac_pb_are_aucune,
               (unsigned long)nbr_fac_pb_are_double,
               (unsigned long)nbr_fac_pb_are_boucle, (unsigned)nbr_max_fac_dec) ;

    if (liste_fac_err != NULL) {

      liste_fac_err->nbr = nbr_fac_pb ;
      BFT_MALLOC(liste_fac_err->val, liste_fac_err->nbr, ecs_int_t) ;

      nbr_fac_pb = 0 ;

      for (ifac = 0 ; ifac < nbr_fac ; ifac++) {
        if (masque_fac_select.val[ifac] < -1)
          liste_fac_err->val[nbr_fac_pb++] = ifac ;
      }
      assert(nbr_fac_pb == liste_fac_err->nbr) ;

    }

  }
  else if (liste_fac_err != NULL) {
    liste_fac_err->nbr = 0 ;
    liste_fac_err->val = NULL ;
  }

  masque_fac_select.nbr = 0 ;
  BFT_FREE(masque_fac_select.val) ;

  liste_are_fac.nbr = 0 ;
  BFT_FREE(liste_are_fac.val) ;
  liste_are_ext.nbr = 0 ;
  BFT_FREE(liste_are_ext.val) ;
  liste_are_int.nbr = 0 ;
  BFT_FREE(liste_are_int.val) ;
  liste_are_tete.nbr = 0 ;
  BFT_FREE(liste_are_tete.val) ;


  ecs_vec_int__detruit(vec_som_are) ;


  cpt_fac_new     = (*vec_fac_old_new)->pos_tab[ifac] - 1 ;
  cpt_val_fac_new = vec_def_fac_new->pos_tab[cpt_fac_new] - 1 ;


  BFT_REALLOC(vec_def_fac_new->pos_tab, cpt_fac_new + 1, ecs_size_t) ;

  BFT_REALLOC(vec_def_fac_new->val_tab, cpt_val_fac_new, ecs_int_t) ;

  vec_def_fac_new->pos_nbr = cpt_fac_new + 1 ;


  BFT_REALLOC((*vec_fac_old_new)->val_tab, cpt_val_fac_new, ecs_int_t) ;


  /* Certaines nouvelles faces intérieures à 2 faces initiales             */
  /*  sont définies 2 fois : il faut compacter la nouvelle liste des faces */
  /*                                                                       */
  /*                  .-----------------------.                            */
  /*          1ere    |                       |                            */
  /*          face    |                       |                            */
  /*          traitée |     .-----------------+----------.                 */
  /*                  |     |                 |          |                 */
  /*                  |     | Face intérieure |          |                 */
  /*                  |     |  définie 2 fois |          |                 */
  /*                  |     |                 |          |                 */
  /*                  `-----+-----------------'          | 2ème            */
  /*                        |                            | face            */
  /*                        |                            | traitée         */
  /*                        `----------------------------'                 */
  /*                                                                       */
  /*                                                                       */

  BFT_MALLOC(tab_signe_fac, 1, ecs_tab_int_t) ;

#if 0 && defined(DEBUG) && !defined(NDEBUG)

  {

    size_t     nbr_are_fac ;
    size_t     pos_are_fac ;


    printf("\n\nrubrique 15-1 ecs_vec_def__recollement\n") ;
    printf("definition des fac en fonction des are new avant fusion\n") ;


    printf("\nvec_def_fac :\n") ;
    for (ifac = 0 ; ifac < vec_def_fac_new->pos_nbr - 1 ; ifac++) {
      pos_are_fac = vec_def_fac_new->pos_tab[ifac]     - 1 ;
      nbr_are_fac
        = vec_def_fac_new->pos_tab[ifac + 1] - 1
        - pos_are_fac ;
      printf("  fac new %d :", ifac + 1) ;
      for (iare = 0 ; iare < nbr_are_fac ; iare++) {
        printf(" %d", vec_def_fac_new->val_tab[pos_are_fac + iare]) ;
      }
      printf("\n") ;
    }
    fflush(stdout);

  }

#endif


  tab_renum_fac = ecs_vec_def__fusionne(&vec_def_fac_new,
                                        tab_signe_fac) ;


#if 0 && defined(DEBUG) && !defined(NDEBUG)

  {

    size_t        ifac ;
    size_t        nbr_are_fac ;
    size_t        pos_are_fac ;

    printf("\n\nrubrique 15-2 ecs_vec_def__recollement\n") ;
    printf("definition des fac en fonction des are new apres fusion\n") ;


    printf("\ntab_renum_fac :\n") ;
    for (ifac = 0 ; ifac < tab_renum_fac.nbr ; ifac++)
      printf("  tab_renum_fac.val[%d] = %d\n", ifac, tab_renum_fac.val[ifac]) ;

    printf("\ntab_signe_fac :\n") ;
    for (ifac = 0 ; ifac < tab_signe_fac->nbr ; ifac++)
      printf("  tab_signe_fac->val[%d] = %d\n", ifac, tab_signe_fac->val[ifac]);


    printf("\nvec_def_fac :\n") ;
    for (ifac = 0 ; ifac < vec_def_fac_new->pos_nbr - 1 ; ifac++) {
      pos_are_fac = vec_def_fac_new->pos_tab[ifac] - 1 ;
      nbr_are_fac
        = vec_def_fac_new->pos_tab[ifac + 1] - 1
        - pos_are_fac ;
      printf("  fac new %d :", ifac + 1) ;
      for (iare = 0 ; iare < nbr_are_fac ; iare++) {
        printf(" %d", vec_def_fac_new->val_tab[pos_are_fac + iare]) ;
      }
      printf("\n") ;
    }
    fflush(stdout);

  }

#endif


  ecs_vec_int__renumerote_et_sgn(*vec_fac_old_new,
                                 tab_renum_fac,
                                 *tab_signe_fac) ;


  BFT_FREE(tab_signe_fac->val) ;
  BFT_FREE(tab_signe_fac) ;

  BFT_FREE(tab_renum_fac.val) ;


  /*
    Renumérotation de la liste des faces posant problème (non découpées,
    et ayant donc une entrée et une seule dans vec_fac_old_new).
   */

  if (liste_fac_err != NULL) {

    for (ifac = 0 ; ifac < liste_fac_err->nbr ; ifac++) {
      ipos = (*vec_fac_old_new)->pos_tab[liste_fac_err->val[ifac]] - 1 ;
      liste_fac_err->val[ifac] = ECS_ABS((*vec_fac_old_new)->val_tab[ipos]) - 1 ;
    }

  }


  return vec_def_fac_new ;


}


/*----------------------------------------------------------------------------
 *  Fonction qui construit des nouvelles faces issues d'une face donnée
 *
 *  Renvoie :
 *    0 en cas de succès ;
 *    1 en cas de "cycle ouvert" (pas d'arête suivante acceptable) ;
 *    2 en cas d'arête parcourue deux fois pour une même face reconstruite ;
 *    3 en cas de boucle (nombre de sous-faces trop important et irréaliste).
 *----------------------------------------------------------------------------*/

static ecs_int_t ecs_loc_vec_def_rc__decoup_fac
(
 const ecs_vec_int_t    *const vec_def_fac,
 const ecs_vec_int_t    *const vec_def_are,
 const ecs_vec_real_t   *const vec_def_som,
 const float            *const normale_fac,
 const ecs_int_t        *const masque_fac_select,
 const ecs_vec_int_t    *const vec_fac_vis,
 const ecs_vec_int_t    *const vec_def_som_are,
 const ecs_vec_int_t    *const vec_def_are_fac,
       ecs_vec_int_t    *const vec_def_fac_new,
       ecs_vec_int_t    *const vec_fac_old_new,
 const ecs_int_t               ifac,
 const ecs_real_t              epsilon_plan,
       size_t                  nbr_max_fac_dec,
       ecs_tab_int_t    *const liste_are_tete,
       ecs_tab_int_t    *const liste_are_ext,
       ecs_tab_int_t    *const liste_are_int,
       ecs_tab_int_t    *const liste_are_fac,
       size_t           *const max_pos_fac_new,
       size_t           *const max_val_fac_new
)
{

#define ECS_FCT_PRODUIT_VECTORIEL(prod_vect, vect1, vect2)   ( \
  prod_vect[0] = vect1[1] * vect2[2] - vect2[1] * vect1[2],    \
  prod_vect[1] = vect2[0] * vect1[2] - vect1[0] * vect2[2],    \
  prod_vect[2] = vect1[0] * vect2[1] - vect2[0] * vect1[1]   )

#define ECS_FCT_COSINUS(vect1, vect2)                         ( \
   ECS_LOC_PRODUIT_SCALAIRE(vect1, vect2)                       \
 / (sqrt(ECS_LOC_MODULE_CARRE(vect1)) * sqrt(ECS_LOC_MODULE_CARRE(vect2))) )


  ecs_bool_t     bool_fin ;
  ecs_bool_t     bool_test_are ;
  ecs_bool_t     bool_trouve ;

  size_t         cpt_are_fac_new ;
  size_t         cpt_are_int ;
  size_t         cpt_are_tete ;
  size_t         cpt_fac_new ;
  size_t         cpt_fac_new_old ;
  size_t         cpt_val_fac_new ;

  size_t         iare ;
  size_t         iare_ext ;
  size_t         iare_fac ;
  size_t         iare_int ;
  size_t         iare_int_bis ;
  size_t         iare_som ;
  size_t         iare_sui ;
  size_t         iare_sup_som ;
  size_t         iare_tete ;
  size_t         icoo ;
  size_t         ifac_are_connect ;
  size_t         ifac_sel ;
  size_t         ipos_vis ;

  size_t         nbr_are_ext ;
  size_t         nbr_are_fac ;
  size_t         nbr_som_are ;

  ecs_int_t      num_are ;
  ecs_int_t      num_are_connect ;
  ecs_int_t      num_are_suite ;
  ecs_int_t      num_are_suite_p_g ;
  ecs_int_t      num_are_suite_n_d ;
  ecs_int_t      num_are_tete ;
  ecs_int_t      num_fac ;
  ecs_int_t      num_som ;
  ecs_int_t      num_som_are ;
  ecs_int_t      num_som_are_prec ;
  ecs_int_t      num_som_extr ;
  ecs_int_t      num_som_extr_p_g ;
  ecs_int_t      num_som_extr_n_d ;
  ecs_int_t      num_som_init ;
  ecs_int_t      num_som_prec ;

  size_t         pos_are_fac ;
  size_t         pos_coo_som ;
  size_t         pos_coo_som_prec ;
  size_t         pos_fac_are_connect_inf ;
  size_t         pos_fac_are_connect_sup ;
  size_t         pos_are_som ;
  size_t         pos_som ;
  size_t         pos_som_are ;
  size_t         pos_som_sup ;

  ecs_int_t      val_som_are ;

  ecs_real_t     are_cos ;
  ecs_real_t     d_fac ;
  ecs_real_t     min_cos_p ;
  ecs_real_t     max_cos_n ;
  ecs_real_t     prod_scal ;
  ecs_real_t     prod_scal2 ;
  ecs_real_t     prod_norme ;

  ecs_point_t    vect_are_connect ;
  ecs_point_t    vect_are_trait ;
  ecs_point_t    vect_normale ;
  ecs_point_t    vect_normale_fac_are_connect ;
  ecs_point_t    prod_vect ;
  ecs_point_t    fac_coord_min ;
  ecs_point_t    fac_coord_max ;

  const ecs_real_t   borne_min_cos_c = -1.1 ;
  const ecs_real_t   borne_max_cos_c =  1.1 ;

  const ecs_real_t   eps_rc_produit_scalaire_min_c = 1.E-7 ;


  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  /* vec_fac_old_new->val_tab devrait contenir la suite 1, 2, 3, ... */

  cpt_fac_new         = vec_fac_old_new->pos_tab[ifac] - 1 ;
  cpt_val_fac_new     = vec_def_fac_new->pos_tab[cpt_fac_new] - 1 ;


  /* Traitement pour une face */
  /* ------------------------ */

  pos_are_fac = vec_def_fac->pos_tab[ifac] - 1 ;

  nbr_are_fac
    = vec_def_fac->pos_tab[ifac + 1] - 1
    - pos_are_fac ;


#if 0 && defined(DEBUG) && !defined(NDEBUG)
  printf("\nFace old %d :\n", ifac + 1) ;
#endif


  ifac_sel = masque_fac_select[ifac] ;

  for (icoo = 0 ; icoo < ECS_DIM_3 ; icoo++)
    vect_normale[icoo] = normale_fac[ifac_sel * ECS_DIM_3 + icoo] ;


  cpt_fac_new_old = 0 ;


  cpt_are_tete = 0 ;

  if (nbr_are_fac >= liste_are_tete->nbr) {
    liste_are_tete->nbr = 2 * (nbr_are_fac + 1) ;
    BFT_REALLOC(liste_are_tete->val, liste_are_tete->nbr, ecs_int_t) ;
  }
  if (nbr_are_fac >= liste_are_ext->nbr) {
    liste_are_ext->nbr = 2 * (nbr_are_fac + 1) ;
    BFT_REALLOC(liste_are_ext->val, liste_are_ext->nbr, ecs_int_t) ;
  }

  for (iare = 0 ; iare < nbr_are_fac ; iare++) {
    num_are = vec_def_fac->val_tab[pos_are_fac + iare] ;
    liste_are_tete->val[cpt_are_tete  ] = num_are ;
    liste_are_ext->val [cpt_are_tete++] = num_are ;
  }


  /* Coordonnées min et max de la face */

  num_are = vec_def_fac->val_tab[pos_are_fac] ;
  if (num_are > 0)
    pos_coo_som
      = (vec_def_are->val_tab[( num_are - 1) * 2    ] - 1) * ECS_DIM_3 ;
  else
    pos_coo_som
      = (vec_def_are->val_tab[(-num_are - 1) * 2 + 1] - 1) * ECS_DIM_3 ;
  for (icoo = 0 ; icoo < ECS_DIM_3 ; icoo++) {
    fac_coord_min[icoo] = vec_def_som->val_tab[pos_coo_som + icoo] ;
    fac_coord_max[icoo] = vec_def_som->val_tab[pos_coo_som + icoo] ;
  }

  for (iare = 1 ; iare < nbr_are_fac ; iare++) {

    num_are = vec_def_fac->val_tab[pos_are_fac + iare] ;
    if (num_are > 0)
      pos_coo_som
        = (vec_def_are->val_tab[( num_are - 1) * 2    ] - 1) * ECS_DIM_3 ;
    else
      pos_coo_som
        = (vec_def_are->val_tab[(-num_are - 1) * 2 + 1] - 1) * ECS_DIM_3 ;
    for (icoo = 0 ; icoo < ECS_DIM_3 ; icoo++) {
      fac_coord_min[icoo] = ECS_MIN(fac_coord_min[icoo],
                                    vec_def_som->val_tab[pos_coo_som + icoo]) ;
      fac_coord_max[icoo] = ECS_MAX(fac_coord_max[icoo],
                                    vec_def_som->val_tab[pos_coo_som + icoo]) ;
    }

  }

  d_fac = fac_coord_max[0] - fac_coord_min[0] ;
  for (icoo = 1 ; icoo < ECS_DIM_3 ; icoo++) {
    d_fac = ECS_MAX(d_fac, fac_coord_max[icoo] - fac_coord_min[icoo]) ;
  }

  for (icoo = 0 ; icoo < ECS_DIM_3 ; icoo++) {
    fac_coord_min[icoo] -= 0.5 * d_fac ;
    fac_coord_max[icoo] += 0.5 * d_fac ;
  }


  /* Tant qu'il existe une liste d'arêtes tête */

  while (cpt_are_tete != 0) {


    nbr_are_ext = cpt_are_tete ;
    cpt_are_int = 0 ;
    iare_tete = 0 ;


    /* Tant que la liste des arêtes tête n'est pas vide */

    while (iare_tete < cpt_are_tete) {

#if 0 && defined(DEBUG) && !defined(NDEBUG)
      printf("  face new %3d :", cpt_fac_new + 1) ;
#endif


      /* Initialisation de l'arête à traiter : 1ere arête tête */
      /*-------------------------------------------------------*/

      num_are_tete = liste_are_tete->val[iare_tete] ;

      num_are = num_are_tete ;


      /* Initialisation du sommet initial : 1er noeud de la 1ere arête tête */
      /*--------------------------------------------------------------------*/

      if (num_are > 0) {
        pos_som = vec_def_are->pos_tab[num_are          - 1]  - 1 ;
        pos_som_sup = pos_som + 1 ;
      }
      else if (num_are < 0) {
        pos_som = vec_def_are->pos_tab[ECS_ABS(num_are) - 1] ;
        pos_som_sup = pos_som - 1 ;
      }
      else {
        assert(num_are != 0) ;
      }

      num_som_init = vec_def_are->val_tab[pos_som] ;


      num_som_prec = num_som_init ;


      /* Initialisation du sommet à traiter : 2nd noeud de l'arête traitée */
      /*-------------------------------------------------------------------*/

      num_som = vec_def_are->val_tab[pos_som_sup] ;


      /* Initialisation de la définition de la nouvelle face */
      /*  avec la 1ère arête                                 */
      /*-----------------------------------------------------*/

      if (*max_val_fac_new <= cpt_val_fac_new) {
        *max_val_fac_new = ECS_MAX(*max_val_fac_new * 2, cpt_val_fac_new + 1) ;
        BFT_REALLOC(vec_def_fac_new->val_tab, *max_val_fac_new, ecs_int_t) ;
      }

      vec_def_fac_new->val_tab[cpt_val_fac_new++] = num_are ;

      cpt_are_fac_new = 0 ;
      liste_are_fac->val[cpt_are_fac_new++] = num_are ;

      if (cpt_are_fac_new >= liste_are_fac->nbr) {
        liste_are_fac->nbr = 2* (cpt_are_fac_new + 1) ;
        BFT_REALLOC(liste_are_fac->val, liste_are_fac->nbr, ecs_int_t) ;
      }


#if 0 && defined(DEBUG) && !defined(NDEBUG)
      printf("num are test %d\n", num_are) ;
#endif


      /* Tant que le sommet à traiter n'est pas le sommet initial */

      while (num_som != num_som_init) {


        /* Recherche des arêtes connectées au sommet à traiter, */
        /*  sauf l'arête en cours                               */
        /*------------------------------------------------------*/

        pos_are_som  = vec_def_som_are->pos_tab[num_som - 1] - 1 ;
        iare_sup_som = vec_def_som_are->pos_tab[num_som    ] - 1 ;

        nbr_som_are = iare_sup_som - pos_are_som ;


        if (nbr_som_are != 2) {


          /* Recherche parmi toutes les arêtes connectées au sommet traité, */
          /*  de l'arête suivante à prendre en compte                       */
          /*  pour la définition de la face                                 */
          /*  c'est l'arête la plus à gauche de l'arête traitée             */

          for (icoo = 0 ; icoo < ECS_DIM_3 ; icoo++)
            vect_are_trait[icoo]
              =   vec_def_som->val_tab[(num_som      - 1) * ECS_DIM_3 + icoo]
                - vec_def_som->val_tab[(num_som_prec - 1) * ECS_DIM_3 + icoo] ;

          min_cos_p = borne_max_cos_c ;
          max_cos_n = borne_min_cos_c ;

          bool_fin = ECS_FALSE ;


          for (iare_som = 0 ;
               iare_som < nbr_som_are && bool_fin == ECS_FALSE ; iare_som++) {

            num_are_connect
              = vec_def_som_are->val_tab[pos_are_som + iare_som];

            if (num_are_connect != ECS_ABS(num_are)) {

              pos_som_are = vec_def_are->pos_tab[num_are_connect - 1]  - 1 ;
              val_som_are = vec_def_are->val_tab[pos_som_are] ;

              if (val_som_are == num_som) {
                num_som_are      = vec_def_are->val_tab[pos_som_are + 1] ;
                num_som_are_prec = val_som_are ;
              }
              else {
                num_som_are      = val_som_are ;
                num_som_are_prec = vec_def_are->val_tab[pos_som_are + 1] ;
              }

              bool_test_are    = ECS_TRUE ;
              pos_coo_som      = (num_som_are - 1) * ECS_DIM_3 ;

              for (icoo = 0 ; icoo < ECS_DIM_3 ; icoo++) {
                if (   (vec_def_som->val_tab[pos_coo_som + icoo]
                        < fac_coord_min[icoo])
                    || (vec_def_som->val_tab[pos_coo_som + icoo]
                        > fac_coord_max[icoo]))
                  bool_test_are = ECS_FALSE ;
              }

              /* bool_test_are = ECS_TRUE ; */
              if (bool_test_are == ECS_TRUE) {

                pos_coo_som_prec = (num_som_are_prec-1)*ECS_DIM_3 ;

                for (icoo = 0 ; icoo < ECS_DIM_3 ; icoo++)
                  vect_are_connect[icoo]
                    =   vec_def_som->val_tab[pos_coo_som      + icoo]
                      - vec_def_som->val_tab[pos_coo_som_prec + icoo] ;


                /* L'arête connectée doit être dans le plan */
                /*   de la face à reconstruire :            */

                pos_fac_are_connect_inf
                  = vec_def_are_fac->pos_tab[num_are_connect - 1] - 1 ;

                pos_fac_are_connect_sup
                  = vec_def_are_fac->pos_tab[num_are_connect    ] - 1 ;

                bool_trouve = ECS_FALSE ;

                for (ifac_are_connect = pos_fac_are_connect_inf ;
                     ifac_are_connect < pos_fac_are_connect_sup
                       && bool_trouve == ECS_FALSE ;
                     ifac_are_connect++) {

                  num_fac
                    = ECS_ABS(vec_def_are_fac->val_tab[ifac_are_connect])
                    - 1 ;

                  if (num_fac != ifac) {

                    if (masque_fac_select[num_fac] > -1) {

                      ifac_sel = masque_fac_select[num_fac] ;

                      for (icoo = 0 ; icoo < ECS_DIM_3 ; icoo++)
                        vect_normale_fac_are_connect[icoo]
                          = normale_fac[ifac_sel * ECS_DIM_3 + icoo] ;

                      prod_scal = ECS_LOC_PRODUIT_SCALAIRE
                                    (vect_normale_fac_are_connect,
                                     vect_normale) ;

                      prod_scal2 = prod_scal * prod_scal ;

                      if (prod_scal2 > epsilon_plan) {

                        if (vec_fac_vis == NULL) {

                          bool_trouve = ECS_TRUE ;

                        }

                        /* Si l'on dispose d'informations sur les faces en
                           vis-à-vis, on s'assure aussi que l'arête appartient
                           bien à une de ces faces */

                        else {

                          for (ipos_vis = vec_fac_vis->pos_tab[ifac] - 1 ;
                               ipos_vis < vec_fac_vis->pos_tab[ifac+1] - 1 ;
                               ipos_vis++) {
                            if (vec_fac_vis->val_tab[ipos_vis] - 1 == num_fac) {
                              bool_trouve = ECS_TRUE;
                              break;
                            }
                          }

                          for (ipos_vis = vec_fac_vis->pos_tab[num_fac] - 1 ;
                               ipos_vis < vec_fac_vis->pos_tab[num_fac+1] - 1 ;
                               ipos_vis++) {
                            if (bool_trouve == ECS_TRUE)
                              break;
                            if (vec_fac_vis->val_tab[ipos_vis] - 1 == ifac)
                              bool_trouve = ECS_TRUE;
                          }

                        }

                      } /* if (prod_scal_2 > epsilon_plan) */

                    } /* if (masque_fac_select[num_fac] > -1) */

                  } /* if (num_fac != ifac) */

                  else {

                    bool_trouve = ECS_TRUE ;

                  }

                }


                if (bool_trouve == ECS_TRUE) {

                  /* L'arête connectée est bien dans le plan */
                  /*  de la face à reconstruire              */

                  /* Recherche parmi toutes les arêtes connectées           */
                  /*  au sommet traite,                                     */
                  /*  de l'arête suivante à prendre en compte               */
                  /*  pour la definition de la face :                       */
                  /*                                                        */
                  /*  (les arêtes sont supposées etre dirigées              */
                  /*   du sommet traité vers l'autre extrémité)             */
                  /*                                                        */
                  /*                            |                           */
                  /*                        \   |   /                       */
                  /*  arête suivante =====>  \  |  /                        */
                  /*                          \ | /      Partie gauche      */
                  /*                           \|/                          */
                  /*  arête traitée ==>  -->----o------  ............       */
                  /*                           /|\                          */
                  /*                          / | \      Partie droite      */
                  /*                         /  |  \                        */
                  /*                        /   |   \                       */
                  /*                                                        */
                  /*                                                        */
                  /* Une arête est dans la partie gauche de l'arête traitée */
                  /*  si elle fait un angle avec l'arête traitée            */
                  /*  dans l'intervalle [0, Pi[                             */
                  /*  c'est-a-dire si le sinus de l'angle est >= 0          */
                  /*  c'est-a-dire si le produit_mixte >= 0 :               */
                  /*                                                        */
                  /*                    ---->   ---->    ->                 */
                  /*                   (are_1 X are_2) . n   >= 0           */
                  /*     ->                                                 */
                  /*  ou n  est le vecteur normal à la face                 */
                  /*                                                        */
                  /* Parmi toutes les arêtes dans la partie gauche,         */
                  /*  l'arête la plus à gauche est celle dont le            */
                  /*  cosinus de l'angle est le plus petit.                 */
                  /*                                                        */
                  /* S'il n'y a aucune arête dans la partie gauche,         */
                  /*  l'arête la plus à gauche est celle dont le            */
                  /*  cosinus de l'angle est le plus grand.                 */
                  /*                                                        */

                  /* produit_mixte                                          */
                  /*  = PRODUIT_MIXTE(vect_are_trait, vect_are_connect);    */

                  ECS_FCT_PRODUIT_VECTORIEL(prod_vect,
                                            vect_are_trait,
                                            vect_are_connect);

                  are_cos = ECS_FCT_COSINUS(vect_are_trait, vect_are_connect);

                  prod_scal
                    = ECS_LOC_PRODUIT_SCALAIRE(prod_vect, vect_normale) ;

                  prod_norme
                    =   sqrt(ECS_LOC_MODULE_CARRE(prod_vect))
                      * sqrt(ECS_LOC_MODULE_CARRE(vect_normale)) ;


                  if (prod_scal
                      >= - eps_rc_produit_scalaire_min_c * prod_norme) {

                    /* L'arête suivante est l'arête la plus à gauche */

                    /* L'arête est dans la partie gauche                  */
                    /*  de l'arête traitée                                */

                    /* Parmi les arêtes à gauche de l'arête traitée       */
                    /*  on doit prendre celle ayant le plus petit cosinus */

                    if (are_cos < min_cos_p) {

                      min_cos_p = are_cos ;

                      if (val_som_are == num_som)
                        num_are_suite_p_g =   num_are_connect ;
                      else
                        num_are_suite_p_g = - num_are_connect ;

                      num_som_extr_p_g = num_som_are ;

                    }
                    /* else : rien à faire */


                  }
                  else {


                    /* L'arête est dans la partie droite                  */
                    /*  de l'arête traitée */

                    /* Parmi les arêtes à droite de l'arête traitée       */
                    /*  on doit prendre celle ayant le plus grand cosinus */

                    if (are_cos > max_cos_n) {

                      max_cos_n = are_cos ;

                      if (val_som_are == num_som)
                        num_are_suite_n_d =   num_are_connect ;
                      else
                        num_are_suite_n_d = - num_are_connect ;

                      num_som_extr_n_d = num_som_are ;

                    }
                    /* else : rien à faire */

                  }

                } /* Fin : si l'arête connectée est dans le meme plan */

              } /* Fin : si bool_test_are (arête ne sort pas de la boîte) */

            } /* Fin : si l'arête n'est pas l'arête traitée */
              /* else : rien à faire */

          } /* Fin : boucle sur les arêtes connectées au sommet traité */


          if (bool_fin == ECS_FALSE) {

            if (min_cos_p < borne_max_cos_c) {

              num_are_suite = num_are_suite_p_g ;
              num_som_extr  = num_som_extr_p_g ;

            }
            else if (max_cos_n > borne_min_cos_c) {

              num_are_suite = num_are_suite_n_d ;
              num_som_extr  = num_som_extr_n_d ;

            }
            else if (   min_cos_p >= borne_max_cos_c
                     && max_cos_n <= borne_min_cos_c) {

              /* Si l'on arrive ici, c'est que l'on a trouvé aucune
                 arête convenable */

              return 1 ;

            }


          } /* Fin : si on n'a pas encore atteint le sommet initial */


        } /* Fin : s'il y plus de 2 arêtes connectées au sommet */
        else {

          for (iare_som = 0 ; iare_som < nbr_som_are ; iare_som++) {

            num_are_connect
              = vec_def_som_are->val_tab[pos_are_som + iare_som];

            if (num_are_connect != ECS_ABS(num_are)) {

              pos_som_are = vec_def_are->pos_tab[num_are_connect - 1]  - 1 ;
              val_som_are = vec_def_are->val_tab[pos_som_are] ;

              if (val_som_are == num_som) {
                num_som_extr = vec_def_are->val_tab[pos_som_are + 1] ;
              }
              else {
                num_som_extr = val_som_are ;
              }

              if (val_som_are == num_som)
                num_are_suite =   num_are_connect ;
              else
                num_are_suite = - num_are_connect ;

            } /* Fin : si l'arête n'est pas l'arête traitée */

          } /* Fin : boucle sur les arêtes connectées au sommet traité */

        } /* Fin : s'il n'y a que deux arêtes connectées au sommet traité */


        /* Ajout à la définition de la face de l'arête la plus à gauche */
        /*--------------------------------------------------------------*/

        if (*max_val_fac_new <= cpt_val_fac_new) {
          *max_val_fac_new = ECS_MAX(*max_val_fac_new * 2,
                                     cpt_val_fac_new + 1) ;
          BFT_REALLOC(vec_def_fac_new->val_tab, *max_val_fac_new, ecs_int_t) ;
        }

        vec_def_fac_new->val_tab[cpt_val_fac_new++] = num_are_suite ;


        /* Afin de ne pas boucler indéfiniment,                        */
        /*  on vérifie que l'arête ne fait pas déja partie de la liste */
        /*  des arêtes définissant la face                             */

        iare_fac = 0 ;
        while (   iare_fac < cpt_are_fac_new
               && ECS_ABS(liste_are_fac->val[iare_fac])
                  != ECS_ABS(num_are_suite))
          iare_fac++ ;

        if (iare_fac != cpt_are_fac_new) {

          /* Problème de reconstruction */

          return 2 ;

        }


        liste_are_fac->val[cpt_are_fac_new++] = num_are_suite ;

        if (cpt_are_fac_new >= liste_are_fac->nbr) {
          liste_are_fac->nbr = 2 * (cpt_are_fac_new + 1) ;
          BFT_REALLOC(liste_are_fac->val, liste_are_fac->nbr, ecs_int_t) ;
        }


        /* Si l'arête suivante fait partie des arêtes têtes restantes */
        /*  on la supprime de la liste des arêtes tête                */
        /*------------------------------------------------------------*/

        iare = iare_tete + 1 ;
        while(iare < cpt_are_tete &&
              ECS_ABS(liste_are_tete->val[iare]) != ECS_ABS(num_are_suite))
          iare++ ;

        for (iare_sui = iare ;
             iare_sui < cpt_are_tete - 1 ;
             iare_sui++) {

          liste_are_tete->val[iare_sui] = liste_are_tete->val[iare_sui + 1] ;

          if (iare_sui >= liste_are_tete->nbr) {
            liste_are_tete->nbr = 2 * (iare_sui + 1) ;
            BFT_REALLOC(liste_are_tete->val, liste_are_tete->nbr, ecs_int_t) ;
          }

        }

        if (iare < cpt_are_tete)
          cpt_are_tete-- ;


        /* Le nouveau sommet à traiter devient le noeud de l'autre */
        /*  extrémite de l'arête la plus à gauche                  */
        /*---------------------------------------------------------*/

        num_som_prec = num_som ;
        num_som      = num_som_extr ;


        /* La nouvelle arête à traiter est l'arête la plus à gauche */
        /*----------------------------------------------------------*/

        assert(num_are_suite != num_are) ;

        num_are = num_are_suite ;


        /* L'arête est-elle une arête interne à la face initiale et          */
        /* participant seulement à la définition d'une seule nouvelle face ? */
        /*-------------------------------------------------------------------*/

        if (nbr_are_ext >= liste_are_ext->nbr) {
          liste_are_ext->nbr = 2 * (nbr_are_ext + 1) ;
          BFT_REALLOC(liste_are_ext->val, liste_are_ext->nbr, ecs_int_t) ;
        }

        iare_ext = 0 ;
        while (iare_ext < nbr_are_ext &&
               ECS_ABS(liste_are_ext->val[iare_ext]) != ECS_ABS(num_are))
          iare_ext++ ;

        if (   iare_ext != 0
            && iare_ext == nbr_are_ext) {

          /* L'arête est une arête interne à la face initiale :      */
          /* - si elle n'est pas deja stockée                        */
          /*   dans la liste des arêtes internes à la face initiale, */
          /*    on la stocke                                         */
          /* - si elle y est deja stockée,                           */
          /*    on la supprime de la liste                           */

          iare_int = 0 ;
          while (iare_int < cpt_are_int &&
                 ECS_ABS(liste_are_int->val[iare_int]) != ECS_ABS(num_are))
            iare_int++ ;

          if (iare_int == cpt_are_int) {

            /* L'arête n'est pas deja stockée                        */
            /*  dans la liste des arêtes internes à la face initiale */

            liste_are_int->val[cpt_are_int++] = num_are ;

            if (cpt_are_int >= liste_are_int->nbr) {
              liste_are_int->nbr = 2* (cpt_are_int + 1) ;
              BFT_REALLOC(liste_are_int->val, liste_are_int->nbr, ecs_int_t) ;
            }

          }
          else {

            /* L'arête est deja stockée :                                */
            /*  elle appartient donc à 2 nouvelles faces                 */
            /* On la supprime de la liste des arêtes internes solitaires */

            cpt_are_int-- ;

            for (iare_int_bis = iare_int ;
                 iare_int_bis < cpt_are_int ;
                 iare_int_bis++)
              liste_are_int->val[iare_int_bis]
                = liste_are_int->val[iare_int_bis + 1] ;

          }

        }
        /* else : rien à faire (l'arête est une arête externe) */


#if 0 && defined(DEBUG) && !defined(NDEBUG)
        printf(" %d", num_are) ;
#endif


      } /* Fin : tant que le sommet à traiter n'est pas le sommet initial */


      /* La construction de la nouvelle face est terminée */


#if 0 && defined(DEBUG) && !defined(NDEBUG)
      printf("\n") ;
#endif


      if (*max_pos_fac_new <= cpt_fac_new + 1) {
        *max_pos_fac_new = ECS_MAX(*max_pos_fac_new * 2, cpt_fac_new + 2) ;
        BFT_REALLOC(vec_def_fac_new->pos_tab, *max_pos_fac_new, ecs_size_t) ;
        BFT_REALLOC(vec_fac_old_new->val_tab, *max_pos_fac_new - 1, ecs_int_t) ;
      }

      vec_def_fac_new->pos_tab[cpt_fac_new + 1]
        = vec_def_fac_new->pos_tab[cpt_fac_new] + cpt_are_fac_new ;

      vec_fac_old_new->val_tab[cpt_fac_new] = cpt_fac_new + 1 ;

      cpt_fac_new_old++ ;

      if (cpt_fac_new_old > nbr_max_fac_dec) {

        /* Si une fac old se découpe en plus de `nbr_max_fac_dec'         */
        /*   fac new, on considère que le recollement de la face a échoué */

        return 3 ;

      }


      cpt_fac_new++ ;

      iare_tete++ ;


    } /* Fin : tant qu'il y a des arêtes tête */


    /* On forme une nouvelle liste d'arêtes tête                       */
    /*  à partir des arêtes internes à la face initiale                */
    /*  et ne participant qu'a la definition d'une seule nouvelle face */

    /* Ces arêtes doivent être prises dans le sens inverse pour */
    /*  la construction des nouvelles faces                     */
    /*  complètement intérieures à la face initiale             */

    cpt_are_tete = 0 ;

    if (cpt_are_int >= liste_are_ext->nbr) {
      liste_are_ext->nbr = 2 * (cpt_are_int + 1) ;
      BFT_REALLOC(liste_are_ext->val, liste_are_ext->nbr, ecs_int_t) ;
    }
    if (cpt_are_int >= liste_are_tete->nbr) {
      liste_are_tete->nbr = 2 * (cpt_are_int + 1) ;
      BFT_REALLOC(liste_are_tete->val, liste_are_tete->nbr, ecs_int_t) ;
    }

    for (iare_int = 0 ; iare_int < cpt_are_int ; iare_int++) {
      liste_are_tete->val[cpt_are_tete  ] = - liste_are_int->val[iare_int] ;
      liste_are_ext->val [cpt_are_tete++] = - liste_are_int->val[iare_int] ;
    }


  } /* Fin : tant qu'il y a une liste d'arêtes tête */


  vec_fac_old_new->pos_tab[ifac + 1]
    = vec_fac_old_new->pos_tab[ifac] + cpt_fac_new_old ;


  return 0 ;

}


/*----------------------------------------------------------------------------
 *  Fonction qui détermine la liste des faces participant au recollement
 *  lorsqu'une définition des faces visibles par d'autres est déjà connue.
 *  Il n'est pas nécessaire que deux faces en vis-à-vis partiel se
 *  connaissent ; il suffit que l'une d'elles connaîsse l'autre (cas d'une
 *  filiation par exemple), même si une connaîssance plus complète n'est
 *  pas gênante.
 *----------------------------------------------------------------------------*/

static ecs_tab_int_t  ecs_loc_vec_def_rc__lis_sel_vis
(
 const ecs_vec_int_t   *const vec_fac_vis
)
{

  ecs_int_t          cpt_fac_sel ;
  size_t             ifac ;

  size_t             nbr_fac ;
  size_t             ipos_vis ;
  ecs_int_t          ifac_vis ;

  ecs_int_t        * profil_fac ;

  ecs_tab_int_t      liste_fac_sel ;


  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  nbr_fac = vec_fac_vis->pos_nbr - 1 ;

  BFT_MALLOC(profil_fac, nbr_fac, ecs_int_t) ;

  for (ifac = 0 ; ifac < nbr_fac ; ifac++)
    profil_fac[ifac] = 0 ;


  /* Boucle principale sur les faces */

  for (ifac = 0 ; ifac < nbr_fac ; ifac++) {

    if (  vec_fac_vis->pos_tab[ifac + 1]
        > vec_fac_vis->pos_tab[ifac   ]) {

      profil_fac[ifac] = 1 ;

      /* Boucle sur les faces en vis-à-vis */

      for (ipos_vis = vec_fac_vis->pos_tab[ifac    ] - 1 ;
           ipos_vis < vec_fac_vis->pos_tab[ifac + 1] - 1 ;
           ipos_vis++                                     ) {

        ifac_vis = ECS_ABS(vec_fac_vis->val_tab[ipos_vis]) - 1 ;

        profil_fac[ifac_vis] = 1 ;

      }

    }

  }

  /* Le profil est maintenant crée */

  cpt_fac_sel = 0 ;

  for (ifac = 0 ; ifac < nbr_fac - 1 ; ifac++)
    if (profil_fac[ifac] != 0)
      cpt_fac_sel++ ;


  liste_fac_sel.nbr = cpt_fac_sel ;
  BFT_MALLOC(liste_fac_sel.val, liste_fac_sel.nbr, ecs_int_t) ;


  cpt_fac_sel = 0 ;

  for (ifac = 0 ; ifac < nbr_fac - 1 ; ifac++)
    if (profil_fac[ifac] != 0)
      liste_fac_sel.val[cpt_fac_sel++] = ifac ;


  /* Libération des tableaux de travail */

  BFT_FREE(profil_fac);


  return liste_fac_sel ;

}


/*----------------------------------------------------------------------------
 *  Fonction qui construit les listes des faces nouvelles et modifiées
 *----------------------------------------------------------------------------*/

static void ecs_loc_vec_def_rc__lis_fac_new
(
 const ecs_vec_int_t    *const vec_def_fac,      /*  -> Déf. faces            */
 const ecs_vec_int_t    *const vec_def_are,      /*  -> Déf. arêtes           */
 const ecs_vec_real_t   *const vec_def_som,      /*  -> Déf. sommets          */
 const ecs_vec_int_t    *const vec_fac_old_new,  /*  -> Renumérotation faces  */
 const ecs_tab_int_t           liste_som_new,    /*  -> Soms. nouv. ou modif. */
       ecs_tab_int_t    *const liste_fac_new,    /* <-  Nouvelles faces       */
       ecs_tab_int_t    *const liste_fac_mod     /* <-  Faces modifiées       */
)
{

  size_t          ifac ;
  size_t          isom ;
  size_t          ipos_mod ;
  size_t          ipos_new ;
  size_t          ival ;
  size_t          nbr_fac ;
  size_t          nbr_som ;
  size_t          pos_are_fac ;
  ecs_int_t       num_are ;
  size_t          pos_som ;

  ecs_int_t      *profil_fac ;
  ecs_bool_t     *profil_som ;


  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  liste_fac_new->nbr = 0 ;
  liste_fac_mod->nbr = 0 ;

  nbr_fac = vec_def_fac->pos_nbr - 1 ;
  nbr_som = vec_def_som->pos_nbr - 1 ;

  BFT_MALLOC(profil_som, nbr_som, ecs_bool_t) ;
  BFT_MALLOC(profil_fac, nbr_fac, ecs_int_t) ;


  for (isom = 0 ; isom < nbr_som ; isom++)
    profil_som[isom] = ECS_FALSE ;

  for (ipos_new = 0 ; ipos_new < liste_som_new.nbr ; ipos_new++)
    profil_som[liste_som_new.val[ipos_new]] = ECS_TRUE ;

  for (ifac = 0 ; ifac < nbr_fac ; ifac++)
    profil_fac[ifac] = 0 ;


  /* Détermination des faces issues de deux anciennes faces */

  for (ival = 0 ;
       ival < vec_fac_old_new->pos_tab[vec_fac_old_new->pos_nbr - 1] - 1 ;
       ival++)

    profil_fac[ECS_ABS(vec_fac_old_new->val_tab[ival]) - 1] += 1 ;


  for (ifac = 0 ; ifac < nbr_fac ; ifac++) {
    if (profil_fac[ifac] < 2)
      profil_fac[ifac] = 0 ;
    else
      liste_fac_new->nbr += 1 ;
  }


  /* Boucle sur les faces pour la détermination de faces modifiées */

  for (ifac = 0 ; ifac < nbr_fac ; ifac++) {


    for (pos_are_fac = vec_def_fac->pos_tab[ifac] - 1 ;
         pos_are_fac < vec_def_fac->pos_tab[ifac + 1] - 1 ;
         pos_are_fac++) {

      num_are = vec_def_fac->val_tab[pos_are_fac] ;


      if (num_are > 0)
        pos_som = vec_def_are->pos_tab[ num_are - 1] - 1 ;
      else if (num_are < 0)
        pos_som = vec_def_are->pos_tab[-num_are - 1] ;

      isom = vec_def_are->val_tab[pos_som] - 1 ;

      if (profil_som[isom] == ECS_TRUE &&
          profil_fac[ifac] == 0) {

        profil_fac[ifac] = 1 ;

        liste_fac_mod->nbr += 1 ;

      }

    }

  }


  BFT_FREE(profil_som) ;

  BFT_MALLOC(liste_fac_new->val, liste_fac_new->nbr, ecs_int_t) ;
  BFT_MALLOC(liste_fac_mod->val, liste_fac_mod->nbr, ecs_int_t) ;

  ipos_new = 0 ;
  ipos_mod = 0 ;

  for (ifac = 0 ; ifac < nbr_fac ; ifac++) {
    if (profil_fac[ifac] > 1)
      liste_fac_new->val[ipos_new++] = ifac ;
    else if (profil_fac[ifac] == 1)
      liste_fac_mod->val[ipos_mod++] = ifac ;
  }


  BFT_FREE(profil_fac) ;


}


/*----------------------------------------------------------------------------
 *  Fonction qui determine les distances maximales associées aux sommets
 *----------------------------------------------------------------------------*/

static ecs_tab_real_t ecs_loc_vec_def_rc__dist_max
(
 const ecs_vec_int_t  *const vec_def_are,
       ecs_vec_real_t *const vec_def_som,
 const ecs_real_t            fraction_dist
)
{

  size_t             iare ;
  size_t             isom ;
  size_t             ind_som_1 ;
  size_t             ind_som_2 ;
  ecs_real_t         lng_are ;

  ecs_tab_real_t     dist_max_som ;


  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  /* Construction des longueurs minimales associées à chaque sommet */


  /* À chaque sommet on associe la longueur de la plus petite   */
  /*  arête qui lui est connectée (arêtes intérieures comprises) */

  dist_max_som.nbr = vec_def_som->pos_nbr - 1 ;
  BFT_MALLOC(dist_max_som.val, dist_max_som.nbr, ecs_real_t) ;

  for (isom = 0 ; isom < dist_max_som.nbr ; isom++)
    dist_max_som.val[isom] = HUGE_VAL ;

  for (iare = 0 ; iare < vec_def_are->pos_nbr - 1 ; iare++) {

    ind_som_1 = vec_def_are->val_tab[iare * 2    ] - 1 ;
    ind_som_2 = vec_def_are->val_tab[iare * 2 + 1] - 1 ;

    lng_are
      = ECS_LOC_DISTANCE((vec_def_som->val_tab) + ind_som_1 * ECS_DIM_3,
                         (vec_def_som->val_tab) + ind_som_2 * ECS_DIM_3 ) ;

    dist_max_som.val[ind_som_1]
      = ECS_MIN(dist_max_som.val[ind_som_1], lng_are * fraction_dist);
    dist_max_som.val[ind_som_2]
      = ECS_MIN(dist_max_som.val[ind_som_2], lng_are * fraction_dist);

  }

  return dist_max_som ;

}


/*----------------------------------------------------------------------------
 *  Fonction qui détermine les points d'intersection lorsqu'on ne connaît
 *  que le tableau de sélection des faces à recoller
 *----------------------------------------------------------------------------*/

#if defined(__HP_cc)   /* On limite le niveau d'optimisation max. */
#pragma OPT_LEVEL 2
#endif

static ecs_rc_lis_inter_are_t ecs_loc_vec_def_rc__int_are_sel
(
 const ecs_vec_int_t   *const vec_def_fac,
 const ecs_vec_int_t   *const vec_def_are,
       ecs_vec_real_t  *const vec_def_som,
 const ecs_tab_int_t   *const tab_fac_select,
 const ecs_tab_real_t         dist_max_som,
       size_t                 nbr_are_ini,
 const ecs_real_t             fraction_dist
)
{

  size_t             iare ;
  size_t             iare_0 ;
  size_t             iare_1 ;
  size_t             iare_intrus ;
  size_t             iare_inf ;
  size_t             iare_mid ;
  size_t             iare_sup ;
  size_t             iare_min ;
  size_t             iare_max ;
  size_t             iare_select_0 ;
  size_t             iare_select_1 ;
  size_t             iare_select_min ;
  size_t             iare_select_max ;

  ecs_int_t          icoo ;
  ecs_int_t          icoo_perm ;
  ecs_int_t          iperm = 0 ;
  ecs_int_t          iloc;
  size_t             ifac_select ;
  size_t             i_inter ;
  size_t             isom ;

  size_t             isom_0[2] ;
  size_t             isom_1[2] ;

  ecs_int_t          ind_are ;

  size_t             nbr_are_fac ;
  size_t             nbr_inter;

  ecs_int_t          num_are ;
  ecs_int_t          num_fac ;

  size_t             pos_are_fac ;
  size_t             pos_som ;

  ecs_int_t        * liste_are_boite ;

  ecs_real_t         val_lim ;

  ecs_real_t         s_inter_0[2] ;
  ecs_real_t         t_inter_1[2] ;

  ecs_point_t      * boite_min = NULL ;
  ecs_point_t      * boite_max = NULL ;
  ecs_point_t        taille_boite_max = {0.0, 0.0, 0.0} ;

  size_t             cpt_are_boite ;
  size_t             nbr_max_are_boite = 10 ;

  ecs_point_t        coord_som_are_0[2] ;
  ecs_point_t        coord_som_are_1[2] ;

  ecs_bool_t         bool_aj_inter[2] ;

  ecs_tab_int_t      liste_are_select ;
  ecs_tab_int_t      tab_ord_are ;

  ecs_int_t          compte_test_inter ;
  ecs_int_t          compte_inter ;

  ecs_rc_lis_inter_are_t      liste_inter_are ;

  const size_t  nbr_are = vec_def_are->pos_nbr - 1 ;


  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/

  compte_test_inter    = 0 ;
  compte_inter         = 0 ;


  bft_printf(_("\n  Search for edge intersections:\n")) ;


  /* Initialisation de la liste des intersections */

  {
    liste_inter_are.nbr = 0 ;
    liste_inter_are.nbr_max = 0 ;
    liste_inter_are.are = NULL ;
    liste_inter_are.inter = NULL ;
  }

  /* Détermination des arêtes sur la zône sélectionnée (filtre) */

  {
    ecs_tab_int_t  profil_are_sel ;
    size_t         nbr_are_select  = 0 ;

    profil_are_sel.nbr = nbr_are ;

    BFT_MALLOC(profil_are_sel.val, profil_are_sel.nbr, ecs_int_t) ;

    for (iare = 0 ; iare < profil_are_sel.nbr ; iare++)
      profil_are_sel.val[iare] = 0 ;

    for (ifac_select = 0 ; ifac_select < tab_fac_select->nbr ; ifac_select++) {

      num_fac = tab_fac_select->val[ifac_select] ;

      pos_are_fac = vec_def_fac->pos_tab[num_fac]     - 1 ;
      nbr_are_fac
        = vec_def_fac->pos_tab[num_fac + 1] - 1
        - pos_are_fac ;

      for (iare = 0 ; iare < nbr_are_fac ; iare++) {

        num_are = vec_def_fac->val_tab[pos_are_fac + iare] ;

        ind_are = ECS_ABS(num_are) - 1 ;

        if (profil_are_sel.val[ind_are] == 0) {

          nbr_are_select += 1 ;
          profil_are_sel.val[ind_are] = 1 ;

        }

      }

    }


    /* dans le cas d'un recollement périodique, on sélectionne aussi
       les arêtes dupliquées se trouvant à la fin de vec_def_are     */

    for (iare = nbr_are_ini ; iare < nbr_are ; iare++) {

      nbr_are_select += 1 ;
      profil_are_sel.val[iare] = 1 ;

    }


    liste_are_select.nbr = nbr_are_select ;

    BFT_MALLOC(liste_are_select.val, liste_are_select.nbr, ecs_int_t) ;

    nbr_are_select = 0 ;

    for (iare = 0 ; iare < profil_are_sel.nbr ; iare++)

      if (profil_are_sel.val[iare] == 1)

        liste_are_select.val[nbr_are_select++] = iare ;

    BFT_FREE(profil_are_sel.val) ;
    profil_are_sel.nbr = 0 ;

  }

  /* Construction des boites englobantes pour les arêtes sélectionnées */

  {
    ecs_point_t taille_boite ;

    BFT_MALLOC(boite_min, liste_are_select.nbr, ecs_point_t) ;
    BFT_MALLOC(boite_max, liste_are_select.nbr, ecs_point_t) ;

    for (iare = 0 ; iare < liste_are_select.nbr ; iare++) {

      ind_are = liste_are_select.val[iare] ;

      isom = vec_def_are->val_tab[vec_def_are->pos_tab[ind_are] - 1] - 1 ;
      pos_som = vec_def_som->pos_pas * isom ;

      for (icoo = 0 ; icoo < ECS_DIM_3 ; icoo++) {
        icoo_perm = (icoo + iperm) % ((ecs_int_t) ECS_DIM_3) ;
        boite_min[iare][icoo] = vec_def_som->val_tab[pos_som + icoo_perm]
                                - (dist_max_som.val[isom] * 1.01) ;
        boite_max[iare][icoo] = vec_def_som->val_tab[pos_som + icoo_perm]
                                + (dist_max_som.val[isom] * 1.01) ;
      }

      isom = vec_def_are->val_tab[vec_def_are->pos_tab[ind_are]    ] - 1 ;
      pos_som = vec_def_som->pos_pas * isom ;

      for (icoo = 0 ; icoo < ECS_DIM_3 ; icoo++) {
        icoo_perm = (icoo + iperm) % ((ecs_int_t) ECS_DIM_3) ;
        boite_min[iare][icoo]
          = ECS_MIN(boite_min[iare][icoo],
                    vec_def_som->val_tab[pos_som + icoo_perm]
                    - (dist_max_som.val[isom] * 1.01)) ;
        boite_max[iare][icoo]
          = ECS_MAX(boite_max[iare][icoo],
                    vec_def_som->val_tab[pos_som + icoo_perm]
                    + (dist_max_som.val[isom] * 1.01)) ;
      }

      for (icoo = 0 ; icoo < ECS_DIM_3 ; icoo++) {
        taille_boite[icoo] = boite_max[iare][icoo] - boite_min[iare][icoo] ;
        if(taille_boite[icoo] > taille_boite_max[icoo])
          taille_boite_max[icoo] = taille_boite[icoo] ;
      }

    }

  }


  /* Ordonnancement des arêtes en fonction de leur boite englobante */

  {
    ecs_vec_real_t * vec_boite_min ;

    vec_boite_min = ecs_vec_real__alloue(liste_are_select.nbr + 1,
                                         liste_are_select.nbr * ECS_DIM_3) ;

    vec_boite_min->pos_pas = ECS_DIM_3 ;

    for (iare = 0 ; iare < liste_are_select.nbr ; iare++) {

      for (icoo = 0 ; icoo < ECS_DIM_3 ; icoo++)

        vec_boite_min->val_tab[iare * ECS_DIM_3 + icoo]
          = boite_min[iare][icoo] ;

    }

    tab_ord_are = ecs_loc_vec_def_rc__trie_elt(vec_boite_min,
                                               0.0) ;

    ecs_vec_real__detruit(vec_boite_min) ;

  }


  /*======================
   * Boucle sur les arêtes
   *======================*/

  BFT_MALLOC(liste_are_boite, nbr_max_are_boite, ecs_int_t) ;

  for (iare = 0 ; iare < liste_are_select.nbr - 1 ; iare++) {

    iare_select_0 = tab_ord_are.val[iare] ;
    iare_0        = liste_are_select.val[iare_select_0] ;

    iare_min = iare + 1 ;
    iare_max = liste_are_select.nbr - 1 ;

    cpt_are_boite = 0 ;

    for (iloc = 0 ; iloc < 2 ; iloc++) {

      isom_0[iloc]
        = vec_def_are->val_tab[vec_def_are->pos_tab[iare_0] + iloc - 1] - 1 ;
      pos_som = vec_def_som->pos_pas * isom_0[iloc] ;

      for (icoo = 0 ; icoo < ECS_DIM_3 ; icoo++)
        coord_som_are_0[iloc][icoo]
          = vec_def_som->val_tab[pos_som + icoo] ;

      for (icoo = 0 ; icoo < ECS_DIM_3 ; icoo++) {
        assert(coord_som_are_0[iloc][icoo] > boite_min[iare_select_0][icoo]);
        assert(coord_som_are_0[iloc][icoo] < boite_max[iare_select_0][icoo]);
      }

    }


    /*
      Recherche des intersections possibles entre arêtes ; l'intersection étant
      symétrique, il suffit de traiter les arêtes i et j pour j > i.
      On utilise un hybride dichotomie - parcours.
    */

    /* Dichotomie pour arête min. parcours */

    iare_inf = iare_min ;
    iare_sup = iare_max ;
    iare_mid = (iare_inf + iare_sup) / 2 ;

    val_lim  =   boite_min[iare_select_0][0] - taille_boite_max[0] ;

    while (iare_inf < iare_sup - 1) {
      if (val_lim < boite_min[tab_ord_are.val[iare_mid]][0])
        iare_sup = iare_mid ;
      else
        iare_inf = iare_mid ;
      iare_mid = (iare_inf + iare_sup) / 2 ;
    }

    iare_min = ECS_MIN(iare_inf, liste_are_select.nbr - 1) ;

    /* Dichotomie pour arête max. parcours */

    iare_inf = iare_min ;
    iare_sup = iare_max ;
    iare_mid = (iare_inf + iare_sup) / 2 ;

    val_lim  = boite_max[iare_select_0][0] ;

    while (iare_inf < iare_sup - 1) {
      if (val_lim < boite_min[tab_ord_are.val[iare_mid]][0])
        iare_sup = iare_mid ;
      else
        iare_inf = iare_mid ;
      iare_mid = (iare_inf + iare_sup) / 2 ;
    }

    iare_max = ECS_MIN(iare_sup, liste_are_select.nbr - 1) ;

    while (iare_min < iare_max) {

      iare_select_min = tab_ord_are.val[iare_min] ;
      iare_select_max = tab_ord_are.val[iare_max] ;

      if (   boite_min[iare_select_min][X] < boite_max[iare_select_0][X]
          && boite_min[iare_select_min][Y] < boite_max[iare_select_0][Y]
          && boite_min[iare_select_min][Z] < boite_max[iare_select_0][Z]
          && boite_max[iare_select_min][X] > boite_min[iare_select_0][X]
          && boite_max[iare_select_min][Y] > boite_min[iare_select_0][Y]
          && boite_max[iare_select_min][Z] > boite_min[iare_select_0][Z]) {

        if (nbr_max_are_boite <= cpt_are_boite) {
          nbr_max_are_boite *= 2 ;
          BFT_REALLOC(liste_are_boite, nbr_max_are_boite, ecs_int_t) ;
        }

        liste_are_boite[cpt_are_boite++] = iare_select_min ;

      }

      iare_min += 1 ;

      if (   boite_min[iare_select_max][X] < boite_max[iare_select_0][X]
          && boite_min[iare_select_max][Y] < boite_max[iare_select_0][Y]
          && boite_min[iare_select_max][Z] < boite_max[iare_select_0][Z]
          && boite_max[iare_select_max][X] > boite_min[iare_select_0][X]
          && boite_max[iare_select_max][Y] > boite_min[iare_select_0][Y]
          && boite_max[iare_select_max][Z] > boite_min[iare_select_0][Z]) {

        if (nbr_max_are_boite <= cpt_are_boite) {
          nbr_max_are_boite *= 2 ;
          BFT_REALLOC(liste_are_boite, nbr_max_are_boite, ecs_int_t) ;
        }

        liste_are_boite[cpt_are_boite++] = iare_select_max ;

      }

      iare_max -= 1 ;

    }

    /* Dernier élément possible (à ne pas oublier) */

    if (iare_min == iare_max) {

      iare_select_min = tab_ord_are.val[iare_min] ;

      if (   boite_min[iare_select_min][X] < boite_max[iare_select_0][X]
          && boite_min[iare_select_min][Y] < boite_max[iare_select_0][Y]
          && boite_min[iare_select_min][Z] < boite_max[iare_select_0][Z]
          && boite_max[iare_select_min][X] > boite_min[iare_select_0][X]
          && boite_max[iare_select_min][Y] > boite_min[iare_select_0][Y]
          && boite_max[iare_select_min][Z] > boite_min[iare_select_0][Z]) {

        if (nbr_max_are_boite <= cpt_are_boite) {
          nbr_max_are_boite *= 2 ;
          BFT_REALLOC(liste_are_boite, nbr_max_are_boite, ecs_int_t) ;
        }

        liste_are_boite[cpt_are_boite++] = iare_select_min ;

      }

      iare_min += 1 ;

    }

    /*
      Traitement des intersections possibles ;
    */

#if 0 && defined(DEBUG) && !defined(NDEBUG)
    printf("\n  arête (%d) de boucle principale : %d\n", iare, iare_0) ;
    printf("    boite_min = [%e ; %e ; %e]\n",
           boite_min[iare_select_0][X],
           boite_min[iare_select_0][Y],
           boite_min[iare_select_0][Z]) ;
    printf("    boite_max = [%e ; %e ; %e]\n",
           boite_max[iare_select_0][X],
           boite_max[iare_select_0][Y],
           boite_max[iare_select_0][Z]) ;
#endif

    for (iare_intrus = 0 ; iare_intrus < cpt_are_boite ; iare_intrus++) {

      iare_select_1 = liste_are_boite[iare_intrus] ;
      iare_1        = liste_are_select.val[iare_select_1] ;

      /*--------------------------------------------------*/
      /* Détermination de l'intersection entre les arêtes */
      /*--------------------------------------------------*/

      /* Test necéssaire en périodicité, les intersections ne sont
         possibles qu'entre une arête initiale et une arête dupliquée */

      if (   nbr_are_ini == nbr_are
          || (   nbr_are_ini < nbr_are
              && (   (iare_0 < nbr_are_ini && nbr_are_ini <= iare_1)
                  || (iare_1 < nbr_are_ini && nbr_are_ini <= iare_0)))) {


        for (iloc = 0 ; iloc < 2 ; iloc++)

          isom_1[iloc] =
            vec_def_are->val_tab[vec_def_are->pos_tab[iare_1] + iloc - 1] - 1 ;

#if 0 && defined(DEBUG) && !defined(NDEBUG)
        printf("\n    intersection possible avec arête (%d) : %d\n",
               iare_select_1, iare_1) ;
        printf("      boite_min = [%e ; %e ; %e]\n",
               boite_min[iare_select_1][X],
               boite_min[iare_select_1][Y],
               boite_min[iare_select_1][Z]) ;
        printf("      boite_max = [%e ; %e ; %e]\n",
               boite_max[iare_select_1][X],
               boite_max[iare_select_1][Y],
               boite_max[iare_select_1][Z]) ;
#endif

        for (iloc = 0 ; iloc < 2 ; iloc++) {

          pos_som = vec_def_som->pos_pas * isom_1[iloc] ;

          for (icoo = 0 ; icoo < ECS_DIM_3 ; icoo++)
            coord_som_are_1[iloc][icoo]
              = vec_def_som->val_tab[pos_som + icoo] ;
        }

        ecs_loc_vec_def_rc__int_are_3d(coord_som_are_0[0],
                                       coord_som_are_0[1],
                                       coord_som_are_1[0],
                                       coord_som_are_1[1],
                                       dist_max_som.val[isom_0[0]],
                                       dist_max_som.val[isom_0[1]],
                                       dist_max_som.val[isom_1[0]],
                                       dist_max_som.val[isom_1[1]],
                                       fraction_dist,
                                       s_inter_0,
                                       t_inter_1,
                                       &nbr_inter) ;
        compte_test_inter++ ;

        assert(nbr_inter <= 2) ;

#if 0 && defined(DEBUG) && !defined(NDEBUG)
        {
          for (i_inter = 0 ; i_inter < nbr_inter ; i_inter++) {
            if (   (0.0 < s_inter_0[i_inter] && s_inter_0[i_inter] < 1.0)
                || (0.0 < t_inter_1[i_inter] && t_inter_1[i_inter] < 1.0))
              printf("      intersection : s = %e ; t = %e\n",
                     s_inter_0[i_inter], t_inter_1[i_inter]) ;
            else
              printf("      recollement  : s = %d ; t = %d\n",
                     (int) s_inter_0[i_inter], (int) t_inter_1[i_inter]) ;
          }
        }
#endif

        /* Suppression intersections "triviales" éventuelles */

        for (i_inter = 0 ; i_inter < nbr_inter ; i_inter++) {
          bool_aj_inter[i_inter] = ECS_TRUE ;
          if (   (   s_inter_0[i_inter] < 0.0
                  && t_inter_1[i_inter] < 0.0
                  && isom_0[0] == isom_1[0])
              || (   s_inter_0[i_inter] < 0.0
                  && t_inter_1[i_inter] > 1.0
                  && isom_0[0] == isom_1[1])
              || (   s_inter_0[i_inter] > 1.0
                  && t_inter_1[i_inter] < 0.0
                  && isom_0[1] == isom_1[0])
              || (   s_inter_0[i_inter] > 1.0
                  && t_inter_1[i_inter] > 1.0
                  && isom_0[1] == isom_1[1]))
            bool_aj_inter[i_inter] = ECS_FALSE ;
        }

        if (nbr_inter == 2 && bool_aj_inter[1] == ECS_FALSE)
          nbr_inter = 1 ;
        else if (nbr_inter == 2 && bool_aj_inter[0] == ECS_FALSE) {
          nbr_inter = 1 ;
          s_inter_0[0] = s_inter_0[1] ;
          t_inter_1[0] = t_inter_1[1] ;
          bool_aj_inter[0] = bool_aj_inter[1] ;
        }
        if (nbr_inter == 1 && bool_aj_inter[0] == ECS_FALSE)
          nbr_inter = 0 ;

        compte_inter += nbr_inter ;

        {
          for (i_inter = 0 ; i_inter < nbr_inter ; i_inter++)

            /* Intersection arête-arête */

            ecs_loc_vec_def_rc__aj_inter(iare_0,
                                         iare_1,
                                         s_inter_0[i_inter],
                                         t_inter_1[i_inter],
                                         &liste_inter_are) ;

        }

      }

    } /* Fin test necéssaire en périodicité */

  }

  bft_printf(_("    Checked possibilities             : %10lu\n"),
             (unsigned long)compte_test_inter) ;
  bft_printf(_("    Intersections found               : %10lu\n"),
             (unsigned long)compte_inter) ;


  /* Libération des tableaux de travail */

  BFT_FREE(boite_min) ;
  BFT_FREE(boite_max) ;

  BFT_FREE(liste_are_boite) ;

  tab_ord_are.nbr = 0 ;
  BFT_FREE(tab_ord_are.val) ;

  liste_are_select.nbr = 0 ;
  BFT_FREE(liste_are_select.val) ;

  liste_inter_are.nbr_max = liste_inter_are.nbr ;
  BFT_REALLOC(liste_inter_are.are, liste_inter_are.nbr_max * 2,
              ecs_rc_inter_are_are_t);
  BFT_REALLOC(liste_inter_are.inter, liste_inter_are.nbr_max,
              ecs_rc_inter_are_inter_t);


  ECS_CHRONO("Joining: intersections done         ") ;

  return liste_inter_are ;

}


/*----------------------------------------------------------------------------
 *  Fonction qui détermine les points d'intersection lorsque la liste des
 *  des faces visibles par d'autres est déjà connue. Il n'est pas
 *  nécessaire que deux faces en vis-à-vis partiel se connaissent ; il
 *  suffit que l'une d'elles connaisse l'autre (cas d'une filiation par
 *  exemple), même si une connaissance plus complète n'est pas gênante.
 *----------------------------------------------------------------------------*/

#if defined(__HP_cc)   /* On limite le niveau d'optimisation max. */
#pragma OPT_LEVEL 2
#endif

static ecs_rc_lis_inter_are_t ecs_loc_vec_def_rc__int_are_vis
(
 const ecs_vec_int_t   *const vec_def_fac,
 const ecs_vec_int_t   *const vec_def_are,
       ecs_vec_real_t  *const vec_def_som,
 const ecs_vec_int_t   *const vec_fac_vis,
 const ecs_tab_real_t         dist_max_som,
 const ecs_real_t             fraction_dist
)
{

  size_t             ifac ;
  size_t             iare ;

  ecs_int_t          icoo ;
  size_t             i_inter ;

  size_t             taille_lst_are_tst ;
  ecs_int_t        * tab_val_are_tst ;
  ecs_int_t        * tab_ptr_are_tst ;
  size_t             pos_lst_are_tst ;
  size_t             ipos_vis ;
  size_t             ifac_vis ;
  size_t             pos_are_fac_vis ;
  size_t             nbr_are_fac_vis ;
  size_t             iare_vis ;
  ecs_int_t          num_are_vis ;
  ecs_int_t          ind_are_0 ;
  ecs_int_t          ind_are_1 ;
  ecs_int_t          ind_som_0_0 ;
  ecs_int_t          ind_som_1_0 ;
  ecs_int_t          ind_som_0_1 ;
  ecs_int_t          ind_som_1_1 ;
  ecs_int_t          ind_are_cmp ;

  size_t             pos_som_0_0 ;
  size_t             pos_som_0_1 ;
  size_t             pos_som_1_0 ;
  size_t             pos_som_1_1 ;

  ecs_point_t        coord_som_0_are_0 ;
  ecs_point_t        coord_som_1_are_0 ;
  ecs_point_t        coord_som_0_are_1 ;
  ecs_point_t        coord_som_1_are_1 ;
  ecs_point_t        boite_min_0 ;
  ecs_point_t        boite_max_0 ;
  ecs_point_t        boite_min_1 ;
  ecs_point_t        boite_max_1 ;

  size_t             nbr_are_fac ;
  size_t             nbr_inter ;

  ecs_int_t          num_are ;
  size_t             pos_are_fac ;

  ecs_real_t         s_inter_0[2] ;
  ecs_real_t         t_inter_1[2] ;

  ecs_bool_t         bool_aj_inter[2] ;

  ecs_int_t          compte_inter ;

  ecs_rc_lis_inter_are_t      liste_inter_are ;

  const size_t   nbr_are = vec_def_are->pos_nbr - 1 ;
  const size_t   nbr_fac = vec_def_fac->pos_nbr - 1 ;


  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/

  compte_inter         = 0 ;


  bft_printf(_("\n  Computation of edge intersections:\n")) ;


  /* Initialisation de la liste des intersections */

  {
    liste_inter_are.nbr = 0 ;
    liste_inter_are.nbr_max = 0 ;
    liste_inter_are.are = NULL ;
    liste_inter_are.inter = NULL ;

  }


  /*
    Initialisation de la liste chaînée d'intersections testées ;
    On surdimensionne, mais on n'initialise que la partie utile pour l'instant
 */

  taille_lst_are_tst = (nbr_are + 1) * 1.5 ;
  BFT_MALLOC(tab_val_are_tst, taille_lst_are_tst, ecs_int_t);
  BFT_MALLOC(tab_ptr_are_tst, taille_lst_are_tst, ecs_int_t) ;

  for (iare = 0 ; iare < nbr_are ; iare++) {
    tab_val_are_tst[iare] = -1 ;
    tab_ptr_are_tst[iare] = -1 ;
  }

  pos_lst_are_tst = nbr_are ;


  compte_inter = 0 ;


  /* Boucle principale sur les faces */
  /* ------------------------------- */

  for (ifac = 0 ; ifac < nbr_fac ; ifac++) {


    /* Boucle sur les faces en vis-à-vis */

    for (ipos_vis = vec_fac_vis->pos_tab[ifac    ] - 1 ;
         ipos_vis < vec_fac_vis->pos_tab[ifac + 1] - 1 ;
         ipos_vis++                                     ) {

      ifac_vis = ECS_ABS(vec_fac_vis->val_tab[ipos_vis]) - 1 ;


      /* On connaît les deux faces, on va boucler sur leurs arêtes */

      pos_are_fac = vec_def_fac->pos_tab[ifac] - 1 ;
      nbr_are_fac
        = vec_def_fac->pos_tab[ifac + 1] - 1
          - pos_are_fac ;

      pos_are_fac_vis = vec_def_fac->pos_tab[ifac_vis] - 1 ;
      nbr_are_fac_vis = vec_def_fac->pos_tab[ifac_vis + 1] - 1 - pos_are_fac_vis ;

      /* Boucle sur les arêtes de la face de boucle principale */

      for (iare = 0 ; iare < nbr_are_fac ; iare++) {

        num_are = ECS_ABS(vec_def_fac->val_tab[pos_are_fac + iare]) ;

        /* Boucle sur les arêtes de la face en vis-à-vis */

        for (iare_vis = 0 ; iare_vis < nbr_are_fac_vis ; iare_vis++) {

          num_are_vis
            = ECS_ABS(vec_def_fac->val_tab[pos_are_fac_vis + iare_vis]) ;

          ind_are_0 = ECS_MIN(num_are, num_are_vis) - 1 ;
          ind_are_1 = ECS_MAX(num_are, num_are_vis) - 1 ;

          ind_som_0_0   = vec_def_are->val_tab[ind_are_0 * 2    ] - 1 ;
          ind_som_1_0   = vec_def_are->val_tab[ind_are_0 * 2 + 1] - 1 ;
          ind_som_0_1   = vec_def_are->val_tab[ind_are_1 * 2    ] - 1 ;
          ind_som_1_1   = vec_def_are->val_tab[ind_are_1 * 2 + 1] - 1 ;

          /*
            On vérifie si une intersection entre ces deux arêtes
            n'a pas déjà été déterminée.
          */

          for (ind_are_cmp = ind_are_0 ;
               (   tab_val_are_tst[ind_are_cmp] != ind_are_1
                && tab_ptr_are_tst[ind_are_cmp] != -1) ;
               ind_are_cmp = tab_ptr_are_tst[ind_are_cmp]) ;

          if (ind_are_cmp != ind_are_1) {

            /* L'intersection entre ces arêtes n'est pas encore connue */

            pos_som_0_0   = ECS_DIM_3 * ind_som_0_0 ;
            pos_som_1_0   = ECS_DIM_3 * ind_som_1_0 ;
            pos_som_0_1   = ECS_DIM_3 * ind_som_0_1 ;
            pos_som_1_1   = ECS_DIM_3 * ind_som_1_1 ;

            for (icoo = 0 ; icoo < ECS_DIM_3 ; icoo++) {
              coord_som_0_are_0[icoo]
                = vec_def_som->val_tab[pos_som_0_0 + icoo] ;
              coord_som_1_are_0[icoo]
                = vec_def_som->val_tab[pos_som_1_0 + icoo] ;
              coord_som_0_are_1[icoo]
                = vec_def_som->val_tab[pos_som_0_1 + icoo] ;
              coord_som_1_are_1[icoo]
                = vec_def_som->val_tab[pos_som_1_1 + icoo] ;
            }

            /*
              On vérifie si les arêtes sont suffisamment voisines
              avant d'appeller le sous-programme de détection
              d'intersections plus coûteux (-> gain CPU estimé à
              environ 20% d'après quelques mesures).
            */

            for (icoo = 0 ; icoo < ECS_DIM_3 ; icoo++) {
              boite_min_0[icoo] = coord_som_0_are_0[icoo]
                                  - (dist_max_som.val[ind_som_0_0] * 1.01) ;
              boite_max_0[icoo] = coord_som_0_are_0[icoo]
                                  + (dist_max_som.val[ind_som_0_0] * 1.01) ;
              boite_min_1[icoo] = coord_som_0_are_1[icoo]
                                  - (dist_max_som.val[ind_som_0_1] * 1.01) ;
              boite_max_1[icoo] = coord_som_0_are_1[icoo]
                                  + (dist_max_som.val[ind_som_0_1] * 1.01) ;
            }

            for (icoo = 0 ; icoo < ECS_DIM_3 ; icoo++) {
              boite_min_0[icoo]
                = ECS_MIN(boite_min_0[icoo],
                          coord_som_1_are_0[icoo]
                          - (dist_max_som.val[ind_som_1_0] * 1.01)) ;
              boite_max_0[icoo]
                = ECS_MAX(boite_max_0[icoo],
                          coord_som_1_are_0[icoo]
                          + (dist_max_som.val[ind_som_1_0] * 1.01)) ;
              boite_min_1[icoo]
                = ECS_MIN(boite_min_1[icoo],
                          coord_som_1_are_1[icoo]
                          - (dist_max_som.val[ind_som_1_1] * 1.01)) ;
              boite_max_1[icoo]
                = ECS_MAX(boite_max_1[icoo],
                          coord_som_1_are_1[icoo]
                          + (dist_max_som.val[ind_som_1_1] * 1.01)) ;
            }

            /* Test si les boites englobantes s'intersectent */

            if (   boite_min_0[X] < boite_max_1[X]
                && boite_min_0[Y] < boite_max_1[Y]
                && boite_min_0[Z] < boite_max_1[Z]
                && boite_max_0[X] > boite_min_1[X]
                && boite_max_0[Y] > boite_min_1[Y]
                && boite_max_0[Z] > boite_min_1[Z]) {

              /*
                On ajoute l'arête 1 à la liste des arêtes testées
                en vis-à-vis de l'arête 1 et on recherche une vraie
                intersection.
              */

              if (tab_val_are_tst[ind_are_cmp] == -1)

                tab_val_are_tst[ind_are_cmp] = ind_are_1 ;

              else {

                if (pos_lst_are_tst == taille_lst_are_tst) {
                  taille_lst_are_tst *= 2 ;
                  BFT_REALLOC(tab_val_are_tst, taille_lst_are_tst,
                              ecs_int_t) ;
                  BFT_REALLOC(tab_ptr_are_tst, taille_lst_are_tst,
                              ecs_int_t) ;
                }

                tab_ptr_are_tst[ind_are_cmp] = pos_lst_are_tst ;
                tab_ptr_are_tst[pos_lst_are_tst] = -1 ;

                tab_val_are_tst[pos_lst_are_tst] = ind_are_1 ;

                pos_lst_are_tst += 1 ;

              }

              ecs_loc_vec_def_rc__int_are_3d(coord_som_0_are_0,
                                             coord_som_1_are_0,
                                             coord_som_0_are_1,
                                             coord_som_1_are_1,
                                             dist_max_som.val[ind_som_0_0],
                                             dist_max_som.val[ind_som_1_0],
                                             dist_max_som.val[ind_som_0_1],
                                             dist_max_som.val[ind_som_1_1],
                                             fraction_dist,
                                             s_inter_0,
                                             t_inter_1,
                                             &nbr_inter) ;

              assert(nbr_inter <= 2) ;

#if 0 && defined(DEBUG) && !defined(NDEBUG)
              {
                for (i_inter = 0 ; i_inter < nbr_inter ; i_inter++) {
                  if (   (   s_inter_0[i_inter] > 0.0
                          && s_inter_0[i_inter] < 1.0)
                      || (   t_inter_1[i_inter] > 0.0
                          && t_inter_1[i_inter] < 1.0))
                    printf("      intersection : s = %e ; t = %e\n",
                           s_inter_0[i_inter], t_inter_1[i_inter]) ;
                  else
                    printf("      recollement  : s = %d ; t = %d\n",
                           (int)s_inter_0[i_inter], (int)t_inter_1[i_inter]) ;
                }
              }
#endif


              /* Suppression intersections "triviales" éventuelles */

              for (i_inter = 0 ; i_inter < nbr_inter ; i_inter++) {

                bool_aj_inter[i_inter] = ECS_TRUE ;
                if (   (   s_inter_0[i_inter] < 0.0
                        && t_inter_1[i_inter] < 0.0
                        && ind_som_0_0 == ind_som_0_1)
                    || (   s_inter_0[i_inter] < 0.0
                        && t_inter_1[i_inter] > 1.0
                        && ind_som_0_0 == ind_som_1_1)
                    || (   s_inter_0[i_inter] > 1.0
                        && t_inter_1[i_inter] < 0.0
                        && ind_som_1_0 == ind_som_0_1)
                    || (   s_inter_0[i_inter] > 1.0
                        && t_inter_1[i_inter] > 1.0
                        && ind_som_1_0 == ind_som_1_1))
                  bool_aj_inter[i_inter] = ECS_FALSE ;

              }

              if (nbr_inter == 2 && bool_aj_inter[1] == ECS_FALSE)
                nbr_inter = 1 ;
              else if (nbr_inter == 2 && bool_aj_inter[0] == ECS_FALSE) {
                nbr_inter = 1 ;
                s_inter_0[0] = s_inter_0[1] ;
                t_inter_1[0] = t_inter_1[1] ;
                bool_aj_inter[0] = bool_aj_inter[1] ;
              }
              if (nbr_inter == 1 && bool_aj_inter[0] == ECS_FALSE)
                nbr_inter = 0 ;


              compte_inter += nbr_inter ;

              {
                for (i_inter = 0 ; i_inter < nbr_inter ; i_inter++)

                  /* Intersection arête-arête */

                  ecs_loc_vec_def_rc__aj_inter(ind_are_0,
                                               ind_are_1,
                                               s_inter_0[i_inter],
                                               t_inter_1[i_inter],
                                               &liste_inter_are) ;

              }

            } /* Fin test si les boites englobantes s'intersectent */

          } /*
              Fin du traitement du cas ou une intersection ind_are_0
              et ind_are_1 n'avait pas encore été testée.
            */

        } /* Fin boucle sur les arêtes de la face en vis-à-vis */

      } /* Fin boucle sur les arêtes de la face de boucle principale */

    } /* Fin boucle sur les faces en vis-à-vis */

  } /* Fin boucle principale sur les faces */

  bft_printf(_("    Intersections found               : %10lu\n"),
             (unsigned long)compte_inter) ;


  /* Libération des tableaux de travail */

  taille_lst_are_tst = 0 ;
  BFT_FREE(tab_val_are_tst);
  BFT_FREE(tab_ptr_are_tst) ;

  liste_inter_are.nbr_max = liste_inter_are.nbr ;
  BFT_REALLOC(liste_inter_are.are, liste_inter_are.nbr_max * 2,
              ecs_rc_inter_are_are_t);
  BFT_REALLOC(liste_inter_are.inter, liste_inter_are.nbr_max,
              ecs_rc_inter_are_inter_t);


  ECS_CHRONO("Joining: intersections done         ") ;

  return liste_inter_are ;

}


/*----------------------------------------------------------------------------
 *  Fonction qui détermine les points d'intersection lorsque les
 *  arêtes à intersecter partagent au départ déjà un sommet.
 *----------------------------------------------------------------------------*/

#if defined(__HP_cc)   /* On limite le niveau d'optimisation max. */
#pragma OPT_LEVEL 2
#endif

static ecs_rc_lis_inter_are_t ecs_loc_vec_def_rc__int_are_som
(
 const ecs_vec_int_t   *const vec_def_fac,
 const ecs_vec_int_t   *const vec_def_are,
       ecs_vec_real_t  *const vec_def_som,
 const ecs_tab_int_t   *const tab_fac_select,
 const ecs_tab_real_t         dist_max_som,
 const ecs_real_t             fraction_dist
)
{

  size_t             iare ;
  size_t             iare_0 ;
  size_t             iare_1 ;

  ecs_int_t          icoo ;
  ecs_int_t          iloc ;
  size_t             ifac_select ;
  size_t             i_inter ;
  size_t             isom ;
  size_t             isom_loc ;
  size_t             ival ;
  size_t             ipos ;
  size_t             ipos_deb ;

  ecs_int_t          isom_0[2] ;
  ecs_int_t          isom_1[2] ;

  ecs_int_t          ind_are ;

  size_t             nbr_are_fac ;
  size_t             nbr_inter;

  ecs_int_t          num_are ;
  ecs_int_t          num_fac ;

  size_t             pos_are_fac ;
  size_t             pos_som ;

  ecs_real_t         s_inter_0[2] ;
  ecs_real_t         t_inter_1[2] ;

  ecs_point_t        coord_som_are_0[2] ;
  ecs_point_t        coord_som_are_1[2] ;

  ecs_bool_t         bool_aj_inter[2] ;

  ecs_tab_int_t      liste_are_select ;

  size_t             compte_test_inter ;
  size_t             compte_inter ;

  ecs_vec_int_t    * vec_som_are ;

  ecs_rc_lis_inter_are_t      liste_inter_are ;

  const size_t    nbr_som = vec_def_som->pos_nbr - 1 ;
  const size_t    nbr_are = vec_def_are->pos_nbr - 1 ;


  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/

  compte_test_inter    = 0 ;
  compte_inter         = 0 ;


  bft_printf(_("\n  Search for edge intersections:\n")) ;


  /* Initialisation de la liste des intersections */

  {
    liste_inter_are.nbr = 0 ;
    liste_inter_are.nbr_max = 0 ;
    liste_inter_are.are = NULL ;
    liste_inter_are.inter = NULL ;
  }

  /* Détermination des arêtes sur la zône sélectionnée (filtre) */

  {
    ecs_tab_int_t  profil_are_sel ;
    size_t         nbr_are_select  = 0 ;

    profil_are_sel.nbr = nbr_are ;

    BFT_MALLOC(profil_are_sel.val, profil_are_sel.nbr, ecs_int_t) ;

    for (iare = 0 ; iare < profil_are_sel.nbr ; iare++)
      profil_are_sel.val[iare] = 0 ;

    for (ifac_select = 0 ; ifac_select < tab_fac_select->nbr ; ifac_select++) {

      num_fac = tab_fac_select->val[ifac_select] ;

      pos_are_fac = vec_def_fac->pos_tab[num_fac]     - 1 ;
      nbr_are_fac
        = vec_def_fac->pos_tab[num_fac + 1] - 1
        - pos_are_fac ;

      for (iare = 0 ; iare < nbr_are_fac ; iare++) {

        num_are = vec_def_fac->val_tab[pos_are_fac + iare] ;

        ind_are = ECS_ABS(num_are) - 1 ;

        if (profil_are_sel.val[ind_are] == 0) {

          nbr_are_select += 1 ;
          profil_are_sel.val[ind_are] = 1 ;

        }

      }

    }


    liste_are_select.nbr = nbr_are_select ;

    BFT_MALLOC(liste_are_select.val, liste_are_select.nbr, ecs_int_t) ;

    nbr_are_select = 0 ;

    for (iare = 0 ; iare < profil_are_sel.nbr ; iare++)

      if (profil_are_sel.val[iare] == 1)

        liste_are_select.val[nbr_are_select++] = iare ;


    /* Construction de la nouvelle connectivité inverse "sommet -> arêtes" */

    vec_som_are = ecs_vec_int__inverse(vec_def_are,
                                       nbr_som) ;


    /*
      On compacte la connectivité sommet -> arêtes de manière à ce qu'elle
      ne contienne plus que les références aux arêtes des faces sélectionnées
      (-> gain CPU et mémoire)
    */

    for (ival = 0 ; ival < vec_som_are->pos_tab[nbr_som] - 1 ; ival++) {
      iare = vec_som_are->val_tab[ival] - 1 ;
      if (profil_are_sel.val[iare] == 0)
        vec_som_are->val_tab[ival] = 0 ;
    }

    ival = 0 ;
    for (isom = 0 ; isom < nbr_som ; isom++) {
      ipos_deb = ival ;
      for (ipos = vec_som_are->pos_tab[isom    ] - 1 ;
           ipos < vec_som_are->pos_tab[isom + 1] - 1 ;
           ipos++) {
        if (vec_som_are->val_tab[ipos] != 0)
          vec_som_are->val_tab[ival++] = vec_som_are->val_tab[ipos] ;
      }
      vec_som_are->pos_tab[isom] = ipos_deb + 1 ;
    }
    vec_som_are->pos_tab[nbr_som] = ival + 1 ;

    ecs_vec_int__redimensionne(vec_som_are, nbr_som + 1, ival) ;


    /* On n'a plus besoin du profil */

    BFT_FREE(profil_are_sel.val) ;
    profil_are_sel.nbr = 0 ;

  }


  /*======================
   * Boucle sur les arêtes
   *======================*/

  for (iare = 0 ; iare < liste_are_select.nbr - 1 ; iare++) {

    iare_0        = liste_are_select.val[iare] ;

    for (iloc = 0 ; iloc < 2 ; iloc++) {

      isom_0[iloc]
        = vec_def_are->val_tab[vec_def_are->pos_tab[iare_0] + iloc - 1] - 1 ;
      pos_som = vec_def_som->pos_pas * isom_0[iloc] ;

      for (icoo = 0 ; icoo < ECS_DIM_3 ; icoo++)
        coord_som_are_0[iloc][icoo]
          = vec_def_som->val_tab[pos_som + icoo] ;

    }


    /*
      Recherche des intersections possibles entre arêtes ; l'intersection étant
      symétrique, il suffit de traiter les arêtes i et j pour j > i.
    */

    for (isom_loc = 0 ; isom_loc < 2 ; isom_loc++) {

      for (ipos = vec_som_are->pos_tab[isom_0[isom_loc]    ] - 1 ;
           ipos < vec_som_are->pos_tab[isom_0[isom_loc] + 1] - 1 ;
           ipos++) {

        iare_1 = vec_som_are->val_tab[ipos] - 1 ;

        if (iare_1 > iare_0) {

          /*--------------------------------------------------*/
          /* Détermination de l'intersection entre les arêtes */
          /*--------------------------------------------------*/


          for (iloc = 0 ; iloc < 2 ; iloc++) {

            isom_1[iloc] =
                vec_def_are->val_tab[vec_def_are->pos_tab[iare_1] + iloc - 1]
              - 1 ;

            pos_som = vec_def_som->pos_pas * isom_1[iloc] ;

            for (icoo = 0 ; icoo < ECS_DIM_3 ; icoo++)
              coord_som_are_1[iloc][icoo]
                = vec_def_som->val_tab[pos_som + icoo] ;

          }

          ecs_loc_vec_def_rc__int_are_3d(coord_som_are_0[0],
                                         coord_som_are_0[1],
                                         coord_som_are_1[0],
                                         coord_som_are_1[1],
                                         dist_max_som.val[isom_0[0]],
                                         dist_max_som.val[isom_0[1]],
                                         dist_max_som.val[isom_1[0]],
                                         dist_max_som.val[isom_1[1]],
                                         fraction_dist,
                                         s_inter_0,
                                         t_inter_1,
                                         &nbr_inter) ;
          compte_test_inter++ ;

          assert(nbr_inter <= 2) ;

#if 0 && defined(DEBUG) && !defined(NDEBUG)
          {
            for (i_inter = 0 ; i_inter < nbr_inter ; i_inter++) {
              if (   (0.0 < s_inter_0[i_inter] && s_inter_0[i_inter] < 1.0)
                  || (0.0 < t_inter_1[i_inter] && t_inter_1[i_inter] < 1.0))
                printf("      intersection : s = %e ; t = %e\n",
                       s_inter_0[i_inter], t_inter_1[i_inter]) ;
              else
                printf("      recollement  : s = %d ; t = %d\n",
                       (int) s_inter_0[i_inter], (int) t_inter_1[i_inter]) ;
            }
          }
#endif

        /* Suppression intersections "triviales" éventuelles */

          for (i_inter = 0 ; i_inter < nbr_inter ; i_inter++) {
            bool_aj_inter[i_inter] = ECS_TRUE ;
            if (   (   s_inter_0[i_inter] < 0.0
                    && t_inter_1[i_inter] < 0.0
                    && isom_0[0] == isom_1[0])
                || (   s_inter_0[i_inter] < 0.0
                    && t_inter_1[i_inter] > 1.0
                    && isom_0[0] == isom_1[1])
                || (   s_inter_0[i_inter] > 1.0
                    && t_inter_1[i_inter] < 0.0
                    && isom_0[1] == isom_1[0])
                || (   s_inter_0[i_inter] > 1.0
                    && t_inter_1[i_inter] > 1.0
                    && isom_0[1] == isom_1[1]))
              bool_aj_inter[i_inter] = ECS_FALSE ;
          }

          if (nbr_inter == 2 && bool_aj_inter[1] == ECS_FALSE)
            nbr_inter = 1 ;
          else if (nbr_inter == 2 && bool_aj_inter[0] == ECS_FALSE) {
            nbr_inter = 1 ;
            s_inter_0[0] = s_inter_0[1] ;
            t_inter_1[0] = t_inter_1[1] ;
            bool_aj_inter[0] = bool_aj_inter[1] ;
          }
          if (nbr_inter == 1 && bool_aj_inter[0] == ECS_FALSE)
            nbr_inter = 0 ;

          compte_inter += nbr_inter ;

          {
            for (i_inter = 0 ; i_inter < nbr_inter ; i_inter++)

              /* Intersection arête-arête */

              ecs_loc_vec_def_rc__aj_inter(iare_0,
                                           iare_1,
                                           s_inter_0[i_inter],
                                           t_inter_1[i_inter],
                                           &liste_inter_are) ;

          }

        }

      }

    } /* Fin boucle sur sommets arête */

  }  /* Fin boucle sur arêtes sélectionnées */

  bft_printf(_("    Checked possibilities             : %10lu\n"),
             (unsigned long)compte_test_inter) ;
  bft_printf(_("    Intersections found               : %10lu\n"),
             (unsigned long)compte_inter) ;


  /* Libération des tableaux de travail */

  liste_are_select.nbr = 0 ;
  BFT_FREE(liste_are_select.val) ;

  liste_inter_are.nbr_max = liste_inter_are.nbr ;
  BFT_REALLOC(liste_inter_are.are, liste_inter_are.nbr_max * 2,
              ecs_rc_inter_are_are_t);
  BFT_REALLOC(liste_inter_are.inter, liste_inter_are.nbr_max,
              ecs_rc_inter_are_inter_t);


  ECS_CHRONO("Joining: intersections doneions     ") ;

  return liste_inter_are ;

}


/*----------------------------------------------------------------------------
 *  Fonction qui ajoute une intersection arête-arête à la liste
 *  d'intersections arêtes-arêtes ;
 *----------------------------------------------------------------------------*/

static void ecs_loc_vec_def_rc__aj_inter
(
 ecs_int_t                 are_0,
 ecs_int_t                 are_1,
 ecs_real_t                s_inter_0,
 ecs_real_t                t_inter_1,
 ecs_rc_lis_inter_are_t  * liste_inter_are
)
{
  size_t  nbr = liste_inter_are->nbr ;
  size_t  ind_0 ;
  size_t  ind_1 ;


  /* Redimensionnement du tableau liste si nécessaire */

  if (liste_inter_are->nbr_max < liste_inter_are->nbr + 1) {

    if (liste_inter_are->nbr_max == 0)
      liste_inter_are->nbr_max += ECS_LOC_RC_NBR_INTER_ARE_INI ;

    else
      liste_inter_are->nbr_max *= 2 ;

    BFT_REALLOC(liste_inter_are->are,   liste_inter_are->nbr_max * 2,
                ecs_rc_inter_are_are_t);
    BFT_REALLOC(liste_inter_are->inter, liste_inter_are->nbr_max,
                ecs_rc_inter_are_inter_t);

  }

  /* Mise à jour de la structure */

  ind_0 = 2*nbr ;
  ind_1 = 2*nbr + 1 ;

  (liste_inter_are->are[ind_0]).are            = are_0 ;
  (liste_inter_are->are[ind_1]).are            = are_1 ;
  (liste_inter_are->are[ind_0]).i_inter        = nbr ;
  (liste_inter_are->are[ind_1]).i_inter        = nbr ;
  (liste_inter_are->are[ind_0]).s_inter        = s_inter_0 ;
  (liste_inter_are->are[ind_1]).s_inter        = t_inter_1 ;

  (liste_inter_are->inter[nbr]).ind_0          = ind_0 ;
  (liste_inter_are->inter[nbr]).ind_1          = ind_1 ;

  liste_inter_are->nbr = ++nbr ;

}


/*----------------------------------------------------------------------------
 *  Fonction qui découpe les arêtes
 *----------------------------------------------------------------------------*/

static ecs_vec_int_t * ecs_loc_vec_def_rc__decoup_are
(
 ecs_vec_int_t           *const vec_def_are,      /* <-> Déf. arêtes          */
 ecs_vec_real_t          *const vec_def_som,      /* <-> Déf. sommets         */
 ecs_rc_lis_inter_are_t         liste_inter_are,  /* <-> Déf. intersections   */
 ecs_tab_int_t           *const tab_equiv_som,    /* <-  Fusions de sommets   */
 ecs_tab_real_t          *const dist_max_som      /* <-> Tolérances           */
)
{

  size_t     nbr_are_old ;
  size_t     nbr_are_max ;
  size_t     nbr_are_new ;
  size_t     nbr_som_old ;
  size_t     nbr_som_max ;
  size_t     nbr_som_new ;

  size_t     ind_are ;
  size_t     ind_are_def ;
  size_t     ind_som_0 ;
  size_t     ind_som_1 ;
  size_t     ind_som_deb ;
  size_t     ind_som_fin ;
  size_t     nbr_inter_are ;
  size_t     ind_inter ;
  size_t     ind_are_inter ;
  size_t     ind_loc_inter ;
  size_t     ind_are_inter_deb ;
  size_t     ind_are_inter_fin ;

  size_t     cpt_ret_fusion ;

  ecs_int_t                 * som_inter ;
  ecs_rc_inter_are_are_t      inter_are_0 ;
  ecs_rc_inter_are_are_t      inter_are_1 ;
  ecs_rc_inter_are_inter_t    inter ;

  ecs_int_t                 * equiv_som_are ;
  size_t                      equiv_som_are_max ;

  ecs_vec_int_t             * vec_are_old_new ;


  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  bft_printf(_("\n  Subdivision of edges:\n")) ;

  /* Initialisations */

  cpt_ret_fusion = 0 ;

  nbr_are_old = vec_def_are->pos_nbr - 1 ;
  nbr_som_old = vec_def_som->pos_nbr - 1 ;

  nbr_are_max = nbr_are_old + (liste_inter_are.nbr * 2) ;
  nbr_som_max = nbr_som_old + (liste_inter_are.nbr * 2) ;

  nbr_are_new = nbr_are_old ;
  nbr_som_new = nbr_som_old ;

  ind_inter = 0 ;

  equiv_som_are_max = 0 ;
  equiv_som_are     = NULL ;

  bft_printf(_("    Initial number of vertices        : %10lu\n"),
             (unsigned long)nbr_som_old) ;
  bft_printf(_("    Initial number of edges           : %10lu\n"),
             (unsigned long)nbr_are_old) ;

  /* On redimensionne vec_def_are et vec_def_som (surdimensionnement) */

  ecs_vec_int__redimensionne(vec_def_are, nbr_are_max + 1, nbr_are_max * 2) ;

  ecs_vec_real__redimensionne(vec_def_som,
                              nbr_som_max + 1,
                              nbr_som_max * ECS_DIM_3) ;

  vec_are_old_new = ecs_vec_int__alloue(nbr_are_old + 1,
                                        nbr_are_old + liste_inter_are.nbr * 2) ;

  dist_max_som->nbr = nbr_som_old + (liste_inter_are.nbr * 2) ;
  BFT_REALLOC(dist_max_som->val, dist_max_som->nbr, ecs_real_t) ;


  /* Préparation du tableau d'équivalence des sommets */

  tab_equiv_som->nbr = nbr_som_max ;
  BFT_MALLOC(tab_equiv_som->val, tab_equiv_som->nbr, ecs_int_t) ;

  for (ind_som_0 = 0 ; ind_som_0 < nbr_som_max ; ind_som_0++)
    tab_equiv_som->val[ind_som_0] = -1 ;


  /* Tableau temporaire pour intersections : */

  BFT_MALLOC(som_inter, liste_inter_are.nbr * 2, ecs_int_t) ;


  /* Première boucle sur les intersections (-> sommets) */
  /*----------------------------------------------------*/

  for (ind_inter = 0 ; ind_inter < liste_inter_are.nbr ; ind_inter++) {

    inter = liste_inter_are.inter[ind_inter] ;

    inter_are_0 = liste_inter_are.are[inter.ind_0] ;

    inter_are_1 = liste_inter_are.are[inter.ind_1] ;


    /* intersection avec sommets existants ou non */

    ind_som_0 = ecs_loc_vec_def_rc__som_inter(vec_def_are,
                                              vec_def_som,
                                              &nbr_som_new,
                                              dist_max_som,
                                              &inter_are_0) ;

    ind_som_1 = ecs_loc_vec_def_rc__som_inter(vec_def_are,
                                              vec_def_som,
                                              &nbr_som_new,
                                              dist_max_som,
                                              &inter_are_1) ;

    som_inter[inter.ind_0] = ind_som_0 ;
    som_inter[inter.ind_1] = ind_som_1 ;

    ecs_vec_def__maj_equiv_som(ind_som_0, ind_som_1, tab_equiv_som) ;

  }

  /* On connaît maintenant le nombre de sommets après intersection */

  nbr_som_max = nbr_som_new ;

  ecs_vec_real__redimensionne(vec_def_som,
                              nbr_som_max + 1,
                              nbr_som_max * ECS_DIM_3) ;

  tab_equiv_som->nbr = nbr_som_max ;
  BFT_REALLOC(tab_equiv_som->val, tab_equiv_som->nbr, ecs_int_t) ;

  bft_printf(_("    Number of vertices after division : %10lu\n"),
             (unsigned long)nbr_som_new) ;


  /* Seconde boucle sur les intersections (-> arêtes) */
  /*--------------------------------------------------*/

  vec_are_old_new->pos_tab[0] = 1 ;

  ind_are = 0 ;
  ind_are_inter_deb = 0 ;

  while (ind_are_inter_deb < liste_inter_are.nbr * 2) {

    /* On passe à la prochaine arête découpée */

    while (ind_are < (liste_inter_are.are[ind_are_inter_deb]).are) {

      vec_are_old_new->pos_tab[ind_are + 1]
        = vec_are_old_new->pos_tab[ind_are] + 1 ;

      vec_are_old_new->val_tab[vec_are_old_new->pos_tab[ind_are] - 1]
        = ind_are + 1 ;

      ind_are++ ;

    }

    /* Initialisation du nombre de positions à 0 */

    vec_are_old_new->pos_tab[ind_are + 1]
      = vec_are_old_new->pos_tab[ind_are] ;

    /* On détermine la sous-liste correspondant à une même arête */

    for (ind_are_inter_fin = ind_are_inter_deb + 1 ;
             ind_are_inter_fin < liste_inter_are.nbr * 2
         && (   (liste_inter_are.are[ind_are_inter_deb]).are
             == (liste_inter_are.are[ind_are_inter_fin]).are) ;
         ind_are_inter_fin++) ;

    nbr_inter_are = ind_are_inter_fin - ind_are_inter_deb ;

    if (nbr_inter_are > equiv_som_are_max) {
      equiv_som_are_max = nbr_inter_are ;
      BFT_REALLOC(equiv_som_are, equiv_som_are_max, ecs_int_t) ;
    }

    cpt_ret_fusion += ecs_loc_vec_def_rc__fusion_inter(vec_def_are,
                                                       vec_def_som,
                                                       liste_inter_are,
                                                       som_inter,
                                                       dist_max_som,
                                                       ind_are_inter_deb,
                                                       nbr_inter_are,
                                                       equiv_som_are);


    /* Création effective des intersections et découpage de l'arête */

    ind_som_deb = vec_def_are->val_tab[ind_are * 2    ] - 1 ;
    ind_som_fin = vec_def_are->val_tab[ind_are * 2 + 1] - 1 ;

    ind_som_0 = ind_som_deb ;

    ind_are_def = ind_are ;

    for (ind_loc_inter = 0 ;
         ind_loc_inter < nbr_inter_are ;
         ind_loc_inter++) {

      ind_are_inter = ind_are_inter_deb + ind_loc_inter ;

      ind_inter = (liste_inter_are.are[ind_are_inter]).i_inter ;
      inter = liste_inter_are.inter[ind_inter] ;

      if ((liste_inter_are.are[inter.ind_0]).are == ind_are)
        ind_som_1 = som_inter[inter.ind_0] ;
      else
        ind_som_1 = som_inter[inter.ind_1] ;

      /* Traitement des équivalences. */

      if (equiv_som_are[ind_loc_inter] != -1) {

        if (ind_som_1 != (size_t)(equiv_som_are[ind_loc_inter]))

          ecs_vec_def__maj_equiv_som(ind_som_1,
                                     equiv_som_are[ind_loc_inter],
                                     tab_equiv_som) ;

        ind_som_1 = ECS_MIN(ind_som_1, (size_t)(equiv_som_are[ind_loc_inter])) ;

      }

      /* Construction d'une arête si le sommet courant n'a pas d'équivalent */

      else if (ind_som_0 != ind_som_1) {

        /* Renumérotation des arêtes */

        vec_are_old_new->val_tab[vec_are_old_new->pos_tab[ind_are + 1] - 1]
          = ind_are_def + 1 ;

        vec_are_old_new->pos_tab[ind_are + 1] += 1 ;

        /* position : (ind_are_def + 1)*2 + 1 = ind_are_def*2 + 3 */

        vec_def_are->pos_tab[ind_are_def   + 1] = ind_are_def*2 + 3 ;
        vec_def_are->val_tab[ind_are_def*2    ] = ind_som_0 + 1 ;
        vec_def_are->val_tab[ind_are_def*2 + 1] = ind_som_1 + 1 ;

        /*
          La première arête issue du découpage prend la place de l'arête
          d'origine (initialisation ind_are_def = ind_are ci-dessus) ;
          les autres sont ajoutées en fin de liste.
        */

        if (ind_are_def < nbr_are_old)
          ind_are_def = nbr_are_new ;
        else
          ind_are_def = ++nbr_are_new ;

        ind_som_0 = ind_som_1 ;

      }

    }

    /*
      On termine la définition si nécessaire (i.e. si pas d'intersection
      en s = 1)
    */

    if (ind_som_1 != ind_som_fin) {

      /* Renumérotation des arêtes */

      vec_are_old_new->val_tab[vec_are_old_new->pos_tab[ind_are + 1] - 1]
        = ind_are_def + 1 ;

      vec_are_old_new->pos_tab[ind_are + 1] += 1 ;

      /* position : (ind_are_def + 1)*2 + 1 = ind_are_def*2 + 3 */

      vec_def_are->pos_tab[ind_are_def   + 1] = ind_are_def*2 + 3 ;
      vec_def_are->val_tab[ind_are_def*2    ] = ind_som_1   + 1 ;
      vec_def_are->val_tab[ind_are_def*2 + 1] = ind_som_fin + 1 ;

      /*
        La première arête issue du découpage prend la place de l'arête
        d'origine (initialisation ind_are_def = ind_are ci-dessus) ;
        les autres sont ajoutées en fin de liste.
      */

      if (ind_are_def >= nbr_are_old)
        nbr_are_new++ ;

    }

    /* On passe au traitement de l'arête suivante */

    ind_are++ ;
    ind_are_inter_deb = ind_are_inter_fin ;

  }

  /* On passe à la dernière arête */

  while (ind_are < nbr_are_old) {

    vec_are_old_new->pos_tab[ind_are + 1]
      = vec_are_old_new->pos_tab[ind_are] + 1 ;

    vec_are_old_new->val_tab[vec_are_old_new->pos_tab[ind_are] - 1]
      = ind_are + 1 ;

    ind_are++ ;

  }

  assert (   vec_are_old_new->pos_tab[vec_are_old_new->pos_nbr - 1] - 1
          == nbr_are_new) ;

  bft_printf(_("    Number of edges after division    : %10lu\n"),
             (unsigned long)nbr_are_new) ;


  /* On redimensionne les tableaux (préalablement surdimensionnés) */

  ecs_vec_int__redimensionne(vec_def_are,
                             nbr_are_new + 1,
                             nbr_are_new * 2) ;

  ecs_vec_int__redimensionne(vec_are_old_new,
                             nbr_are_old + 1,
                             nbr_are_new) ;


#if 0 && defined(DEBUG) && !defined(NDEBUG)

  ecs_dbg_vec_def_rc__imp_vec_int("ecs_loc_vec_def_rc__decoup_are",
                                  "vec_are_old_new",
                                  vec_are_old_new) ;

#endif


  /* Libération mémoire et retour */

  BFT_FREE(som_inter) ;

  BFT_FREE(equiv_som_are) ;


  if (cpt_ret_fusion > 0) {
    ecs_warn() ;
    bft_printf
      (_("The intersection merging tolerance as been reduced for\n"
         "%ld edges so as to avoid crossing the local intersection\n"
         "tolerance through transitivity.\n"
         "-> Eventually reduce tolerance for this joining"),
       (long)cpt_ret_fusion) ;
  }


  return vec_are_old_new ;

}


/*----------------------------------------------------------------------------
 *  Fonction qui détermine les sommets à fusionner pour le découpage
 * des arêtes ; On renvoie 0 si aucune transitivité de fusion n'a été
 * détectée, et 1 si une ou plusieurs transitivités ont du être éliminées
 * par réduction de tolérance ;
 *----------------------------------------------------------------------------*/

static ecs_int_t ecs_loc_vec_def_rc__fusion_inter
(
 ecs_vec_int_t           *const vec_def_are,             /*  -> Déf. arêtes  */
 ecs_vec_real_t          *const vec_def_som,             /*  -> Déf. sommets */
 ecs_rc_lis_inter_are_t         liste_inter_are,         /*  -> Intersects.  */
 ecs_int_t               *const som_inter,               /*  -> Sommets/int. */
 ecs_tab_real_t          *const dist_max_som,            /*  -> Tolérances   */
 ecs_int_t                      ind_are_inter_deb,       /*  -> Début arête  */
 size_t                         nbr_inter_are,           /*  -> Nb. inter    */
 ecs_int_t               *const equiv_som_are            /* <-> Fusions      */
)
{

  size_t     ind_are ;
  size_t     ind_som_0 ;
  size_t     ind_som_1 ;
  size_t     ind_inter ;
  size_t     ind_are_inter ;
  size_t     ind_loc_inter ;

  ecs_real_t  d_min ;
  ecs_real_t  dist_min_inter ;
  ecs_real_t  dist_fusion_max ;
  ecs_real_t  dist_fusion_lim ;
  ecs_real_t  s_inter ;
  ecs_real_t  s_deb ;
  ecs_real_t  s_fin ;
  ecs_bool_t  reduire_dist_fusion ;
  ecs_real_t  lng_are ;

  ecs_rc_inter_are_inter_t    inter ;

  ecs_int_t  val_ret ;

  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  val_ret = 0 ;

#if 0 && defined(DEBUG) && !defined(NDEBUG)
  printf("arête %d ; intersections : %d\n",
         (liste_inter_are.are[ind_are_inter_deb]).are,
         nbr_inter_are) ;
#endif

  /* Initialisation pour la recherche d'intersections à fusionner */

  reduire_dist_fusion = ECS_FALSE ;

  ind_are = (liste_inter_are.are[ind_are_inter_deb]).are ;
  ind_som_0 = vec_def_are->val_tab[vec_def_are->pos_tab[ind_are] - 1] - 1 ;
  ind_som_1 = vec_def_are->val_tab[vec_def_are->pos_tab[ind_are]    ] - 1 ;
  lng_are
    = ECS_LOC_DISTANCE((vec_def_som->val_tab) + ind_som_0 * ECS_DIM_3,
                       (vec_def_som->val_tab) + ind_som_1 * ECS_DIM_3 ) ;

  equiv_som_are[0] = -1 ;

  dist_fusion_lim = lng_are ;

  /*
    On vérifie que la transitivité des recollements n'induise pas de
    recollement de sommets non voisins. Si c'est le cas, on baisse
    localement la tolérance.
  */

  do {

    reduire_dist_fusion = ECS_FALSE ;

    dist_fusion_max  = 0.0 ;

    /* Réinitialisation du tableau de fusion local de sommets */

    for (ind_loc_inter = 1 ;
         ind_loc_inter < nbr_inter_are ;
         ind_loc_inter++) {
      equiv_som_are[ind_loc_inter] = -1 ;
    }

    /* Autres initialisations */

    ind_inter = (liste_inter_are.are[ind_are_inter_deb]).i_inter ;
    inter = liste_inter_are.inter[ind_inter] ;

    ind_som_0 = som_inter[inter.ind_0] ;
    ind_som_1 = som_inter[inter.ind_1] ;

    dist_min_inter = ECS_MIN(dist_max_som->val[ind_som_0],
                             dist_max_som->val[ind_som_1]) ;
    dist_min_inter = ECS_MIN(dist_min_inter,
                             dist_fusion_lim) ;

    d_min = ECS_MIN(dist_min_inter, dist_fusion_lim) ;

    s_deb = liste_inter_are.are[ind_are_inter_deb].s_inter ;
    s_fin = s_deb ;

    /* Recherche des sommets à fusionner */

    for (ind_loc_inter = 1 ;
         ind_loc_inter < nbr_inter_are ;
         ind_loc_inter++) {

      ind_are_inter = ind_are_inter_deb + ind_loc_inter ;

      ind_inter = (liste_inter_are.are[ind_are_inter]).i_inter ;
      inter = liste_inter_are.inter[ind_inter] ;

      if ((liste_inter_are.are[inter.ind_0]).are == ind_are) {
        ind_som_0 = som_inter[inter.ind_0] ;
        ind_som_1 = som_inter[inter.ind_1] ;
      }
      else {
        ind_som_0 = som_inter[inter.ind_1] ;
        ind_som_1 = som_inter[inter.ind_0] ;
      }

      dist_min_inter = ECS_MIN(dist_max_som->val[ind_som_0],
                               dist_max_som->val[ind_som_1]) ;
      dist_min_inter = ECS_MIN(dist_min_inter,
                               dist_fusion_lim) ;

      d_min = ECS_MIN(d_min, dist_min_inter) ;

      s_inter = liste_inter_are.are[ind_are_inter].s_inter ;


      /* On peut fusionner l'intersection avec la précédente */

      if ((s_inter - s_deb) * lng_are <= d_min) {

        equiv_som_are[ind_loc_inter - 1] = ind_som_0 ;

        if ((s_inter - s_fin) * lng_are < dist_fusion_lim)
          dist_fusion_max = ECS_MAX(dist_fusion_max,
                                    (s_inter - s_fin) * lng_are) ;

#if 0 && defined(DEBUG) && !defined(NDEBUG)
        printf("inter %d : fusion avec précédente (s = %f, s = %f)\n",
               ind_loc_inter, s_fin, s_inter);
#endif

        s_fin = s_inter ;

      }

      /*
        On pourrait fusionner l'intersection avec la précédente,
        mais cela poserait un problème de fusions en chaîne
        (i.e. par transitivité) entre sommets trop distants
      */

      else if ((s_inter - s_fin) * lng_are <= dist_min_inter) {

        if ((s_inter - s_fin) * lng_are < dist_fusion_lim)
          dist_fusion_max = ECS_MAX(dist_fusion_max,
                                    (s_inter - s_fin) * lng_are) ;

        dist_fusion_lim = dist_fusion_max * 0.99 ;

#if 0 && defined(DEBUG) && !defined(NDEBUG)
        printf("inter %d : transitivité fusion (s = %f à s = %f)\n"
               " -> réduction tolérance à %f et reprise\n",
               ind_loc_inter, s_fin, s_inter, dist_fusion_lim);
#endif

        reduire_dist_fusion = ECS_TRUE ;

        val_ret = 1 ;

        break ;

      }

      /*
        On ne peut fusionner l'intersection avec la précédente ;
        on repart donc sur un nouveau groupe de sommets.
      */

      else {

        d_min = lng_are ;

        s_deb = s_inter ;

#if 0 && defined(DEBUG) && !defined(NDEBUG)
        printf("inter %d : pas de fusion avec précédente (s = %f, s = %f)\n",
               ind_loc_inter, s_fin, s_inter) ;
#endif

        s_fin = s_inter ;

      }

    }

  } while (reduire_dist_fusion == ECS_TRUE) ;


  /*
    Retour 0 si aucune transitivité détectée, et 1 si la tolérance à du
    être réduite localement pour éviter des transitivités.
  */

  return val_ret ;

}


/*----------------------------------------------------------------------------
 *  Fonction qui supprime les éventuelles arêtes dégénérées
 *----------------------------------------------------------------------------*/

static void  ecs_loc_vec_def_rc__nettoie_are
(
 ecs_vec_int_t           *const vec_def_are,      /* <-> Déf. arêtes          */
 ecs_vec_int_t           *const vec_are_old_new   /* <-> Renumérotation       */
)
{

  size_t     nbr_are_old ;
  size_t     nbr_are_ini ;
  size_t     nbr_are_new ;

  size_t     cpt_are ;
  size_t     ind_are ;
  size_t     ind_som_0 ;
  size_t     ind_som_1 ;

  size_t     cpt_pos_new ;
  size_t     ind_pos_deb ;
  size_t     ind_pos_fin ;
  size_t     ind_loc ;
  ecs_int_t  ind_are_new ;
  ecs_int_t  sgn ;

  ecs_int_t  * renum_are ;

  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  bft_printf(_("\n  Removal of degenerate edges if present      :\n")) ;

  /* Initialisations */

  nbr_are_ini = vec_def_are->pos_nbr - 1 ;
  nbr_are_old = vec_are_old_new->pos_nbr - 1 ;

  nbr_are_new = nbr_are_ini ;

  bft_printf(_("    Initial number of edges           : %10lu\n"),
             (unsigned long)nbr_are_ini) ;


  /* Boucle sur les nouvelles arêtes (renumérotation) */
  /* ------------------------------------------------ */

  BFT_MALLOC(renum_are, nbr_are_ini, ecs_int_t) ;

  cpt_are = 0 ;

  for (ind_are = 0 ; ind_are < nbr_are_ini ; ind_are++) {

    ind_som_0 = vec_def_are->val_tab[ind_are * 2    ] - 1 ;
    ind_som_1 = vec_def_are->val_tab[ind_are * 2 + 1] - 1 ;

    if (ind_som_0 != ind_som_1) {
      vec_def_are->val_tab[cpt_are * 2    ] = ind_som_0 + 1 ;
      vec_def_are->val_tab[cpt_are * 2 + 1] = ind_som_1 + 1 ;
      renum_are[ind_are] = cpt_are++ ;
    }
    else
      renum_are[ind_are] = -1 ;

  }

  nbr_are_new = cpt_are ;

  bft_printf(_("    Number of edges after processing  : %10lu\n"),
             (unsigned long)nbr_are_new) ;


  /* Si on n'a pas d'arêtes dégénérées, on peut sortir */

  if (nbr_are_new == nbr_are_ini) {

    BFT_FREE(renum_are) ;

    return ;

  }


  /* On redimensionne le tableau des arêtes */

  ecs_vec_int__redimensionne(vec_def_are,
                             nbr_are_new + 1,
                             nbr_are_new * 2) ;


  /* Boucle sur les anciennes arêtes */
  /* ------------------------------- */

  cpt_pos_new = 0 ;

  for (ind_are = 0 ; ind_are < nbr_are_old ; ind_are++) {

    ind_pos_deb = vec_are_old_new->pos_tab[ind_are    ] - 1 ;
    ind_pos_fin = vec_are_old_new->pos_tab[ind_are + 1] - 1 ;

    vec_are_old_new->pos_tab[ind_are] =  cpt_pos_new + 1 ;

    for (ind_loc = ind_pos_deb ; ind_loc < ind_pos_fin ; ind_loc++) {

      ind_are_new = renum_are[ECS_ABS(vec_are_old_new->val_tab[ind_loc]) - 1] ;

      if (ind_are_new != -1) {

        if (vec_are_old_new->val_tab[ind_loc] > 0)
          sgn =  1 ;
        else
          sgn = -1 ;
        vec_are_old_new->val_tab[cpt_pos_new++] =  sgn * (ind_are_new + 1) ;

      }

    }

    if (vec_are_old_new->pos_tab[ind_are] == cpt_pos_new + 1) {

      bft_error
        (__FILE__, __LINE__, 0,
         _("Error in conformal joining:\n"
           "all edges based on edge %ld are degenerate\n"
           "Eventually reduce the merging parameter (\"-fraction\")"),
         (long)(ind_are + 1)) ;

    }

  }

  vec_are_old_new->pos_tab[nbr_are_old] =  cpt_pos_new + 1 ;


  /* Libération et redimensionnement mémoire */

  BFT_FREE(renum_are) ;

  ecs_vec_int__redimensionne(vec_are_old_new,
                             nbr_are_old + 1,
                             cpt_pos_new) ;


#if 0 && defined(DEBUG) && !defined(NDEBUG)

  ecs_dbg_vec_def_rc__imp_vec_int("ecs_loc_vec_def_rc__nettoie_are",
                                  "vec_are_old_new",
                                  vec_are_old_new) ;

#endif

}


/*----------------------------------------------------------------------------
 *  Fonction qui enlève les éventuels coins dégénérés des faces
 *----------------------------------------------------------------------------*/

static void  ecs_loc_vec_def_rc__nettoie_fac
(
 ecs_vec_int_t           *const vec_def_fac       /* <-> Déf. faces           */
)
{


  size_t     nbr_are_cur ;
  size_t     nbr_are_ini ;
  size_t     nbr_are_fac_max ;

  size_t     ind_fac ;

  size_t     cpt_fac_val ;
  size_t     cpt_fac_smp ;
  size_t     cpt_are_cur ;
  size_t     cpt_are_del ;
  size_t     ind_pos_deb ;
  size_t     ind_pos_fin ;
  size_t     ind_loc ;

  ecs_int_t  * num_are_fac ;
  ecs_int_t  * del_are_fac ;

  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  bft_printf(_("\n  Removal of degenerate corners if present :\n")) ;

  /* Initialisations */

  nbr_are_fac_max = 0 ;
  num_are_fac     = NULL ;
  del_are_fac     = NULL ;

  cpt_fac_smp     = 0 ;
  cpt_fac_val     = 0 ;


  /*   Dans la définition des faces en fonction des nouvelles arêtes,  */
  /*  une même arête peut être référencée 2 fois en sens inverse,      */
  /*  ceci provenant de la fusion des sommets.                         */
  /*                                                                   */
  /*     x                                      x                      */
  /*     |\                                     |                      */
  /*     | \                                    |                      */
  /*   a2|  \a3                               A2|                      */
  /*     |   \                Fusion            |                      */
  /*     |    \      a4         des             |                      */
  /*  ---s1----s2------       sommets           x                      */
  /*     |      \             s1 et s2         / \                     */
  /*     |       \                            /   \                    */
  /*   a1|        \a4                      A1/     \A3                 */
  /*     |         \                        /       \                  */
  /*     |          \                      /         \                 */
  /*     x-----------x                    x-----------x                */
  /*          a5                                A4                     */
  /*                                                                   */
  /*                                                                   */
  /*  Face : a1 a2 a3 a4 a5           Face : A1 A2 -A2 A3 A4           */
  /*                                                                   */
  /*                                                                   */
  /* Attention : la configuration finale peut être                     */
  /*             A2 A1 A3 A4 -A2                                       */
  /*             ou les deux références d'arêtes à supprimer           */
  /*             sont en début et en fin de la définition de la face.  */
  /*                                                                   */
  /* Remarque :  plusieurs paires d'arêtes peuvent éventuellement      */
  /*             être référencées deux foix, sous la forme             */
  /*             ... A1 A2 -A2 -A1 ... (où la suppression de A2 fait   */
  /*             apparaître ... A1 -A1 ...) ; On effectue donc autant  */
  /*             de passages que nécessaire sur une face.              */
  /*                                                                   */


  /* Boucle sur les faces */
  /* --------------------- */

  for (ind_fac = 0 ; ind_fac < vec_def_fac->pos_nbr - 1 ; ind_fac++) {

    ind_pos_deb = vec_def_fac->pos_tab[ind_fac    ] - 1 ;
    ind_pos_fin = vec_def_fac->pos_tab[ind_fac + 1] - 1 ;

    nbr_are_ini = ind_pos_fin - ind_pos_deb ;

    if (nbr_are_fac_max < nbr_are_ini) {
      nbr_are_fac_max = nbr_are_ini ;
      BFT_REALLOC(num_are_fac, nbr_are_fac_max, ecs_int_t) ;
      BFT_REALLOC(del_are_fac, nbr_are_fac_max, ecs_int_t) ;
    }

    for (ind_loc = 0 ; ind_loc < nbr_are_ini ; ind_loc++) {
      num_are_fac[ind_loc] = vec_def_fac->val_tab[ind_pos_deb + ind_loc] ;
      del_are_fac[ind_loc] = 0 ;
    }

    /* Nettoyage de la définition locale */

    nbr_are_cur = nbr_are_ini ;

    do {

      /* détection dégénerescence */

      for (ind_loc = 0 ; ind_loc < nbr_are_cur ; ind_loc++)
        del_are_fac[ind_loc] = 0 ;

      for (ind_loc = 0 ; ind_loc < nbr_are_cur - 1 ; ind_loc++) {

        if (num_are_fac[ind_loc] == -num_are_fac[ind_loc + 1]) {
          del_are_fac[ind_loc    ] = 1 ;
          del_are_fac[ind_loc + 1] = 1 ;
        }

        if (num_are_fac[0] == -num_are_fac[nbr_are_cur - 1]) {
          del_are_fac[0              ] = 1 ;
          del_are_fac[nbr_are_cur - 1] = 1 ;
        }

      }

      /* On supprime les dégénérescences */

      cpt_are_del = 0 ;
      cpt_are_cur = 0 ;

      for (ind_loc = 0 ; ind_loc < nbr_are_cur ; ind_loc++) {

        if (del_are_fac[ind_loc] == 0)
          num_are_fac[cpt_are_cur++] = num_are_fac[ind_loc] ;
        else
          cpt_are_del++ ;

      }

      nbr_are_cur  = cpt_are_cur ;

    } while (cpt_are_del > 0) ;


    /* La définition locale est maintenant nettoyée */

    if (nbr_are_cur == 0)
      bft_error
        (__FILE__, __LINE__, 0,
         _("Error in conformal joining:\n"
           "face %lu is degenerate\n"
           "Eventually reduce the merging parameter (\"-fraction\")"),
         (unsigned long)(ind_fac + 1)) ;

    else if (nbr_are_cur != nbr_are_ini)
      cpt_fac_smp += 1 ;

    vec_def_fac->pos_tab[ind_fac] =  cpt_fac_val + 1 ;

    for (ind_loc = 0 ; ind_loc < nbr_are_cur ; ind_loc++)
      vec_def_fac->val_tab[cpt_fac_val++] = num_are_fac[ind_loc] ;

  }

  vec_def_fac->pos_tab[vec_def_fac->pos_nbr - 1] =  cpt_fac_val + 1 ;


  bft_printf(_("    Number of simplified faces        : %10lu\n"),
             (unsigned long)cpt_fac_smp ) ;


  /* Libération et redimensionnement mémoire */

  BFT_FREE(num_are_fac) ;
  BFT_FREE(del_are_fac) ;

  ecs_vec_int__redimensionne(vec_def_fac,
                             vec_def_fac->pos_nbr,
                             cpt_fac_val) ;


#if 0 && defined(DEBUG) && !defined(NDEBUG)

  ecs_dbg_vec_def_rc__imp_vec_int("ecs_loc_vec_def_rc__nettoie_fac",
                                  "vec_def_fac (after degen. corner removal)",
                                  vec_def_fac) ;

#endif


}


/*----------------------------------------------------------------------------
 *  Fonction qui détermine le ou les points d'intersection de 2 arêtes
 * dans l'espace 3D. On cherche en premier lieu à établir une correspondance
 * entre sommets, puis une correspondance entre sommets, et enfin une
 * ``vraie intersection,'' le tout sous une tolérance donnée.
 * On renvoie les coordonnées locales s et t d'une intersection. Si s ou t
 * = 0 ou 1, on met sa valeur à -0.01 ou 1.01 respectivement, afin de
 * pouvoir remplacer un test d'égalité à 0 ou 1 par un test d'inégalité,
 * sachant que 0 <= s, t <= 1 pour toute intersection.
 *----------------------------------------------------------------------------*/

static void ecs_loc_vec_def_rc__int_are_3d
(
 const ecs_point_t   som_0_are_0,        /* --> Coordonnées des sommets       */
 const ecs_point_t   som_1_are_0,
 const ecs_point_t   som_0_are_1,
 const ecs_point_t   som_1_are_1,
 const ecs_real_t    d_max_s0_a0,        /* --> Tolérances sur les distances  */
 const ecs_real_t    d_max_s1_a0,
 const ecs_real_t    d_max_s0_a1,
 const ecs_real_t    d_max_s1_a1,
 const ecs_real_t    fraction_dist,
       ecs_real_t    s_inter_0[2],       /* <-- Coord. s intersections seg. 0 */
       ecs_real_t    t_inter_1[2],       /* <-- Coord. t intersections seg. 1 */
       size_t      * nbr_inter           /* <-- Nombre d'intersections        */
)
{
  /*
   * Soient deux arêtes a0 et a1 d'extrémités (P00, P01) et (P10, P11).
   * Les coordonnées d'un point sur l'arête a0 sont de la forme :
   *   P00 + s.(P01 - P00), avec 0 <= s <= 1.
   * De même, les coordonées d'un point de a1 s'écrivent :
   *   P10 + t.(P11 - P10), avec 0 <= t <= 1.
   * La distance de ces deux points elevée au carré s'écrit donc :
   *   d(s, t) = || P10 - P00 + t.(P11 - P10) - s.(P01 - P00) ||^2
   * que l'on peut écrire sous la forme :
   *   d(s, t) = a.s^2 + 2b.s.t + c.t^2 + 2d.s + 2e.t + f
   * On cherche les valeurs de s et t (comprises entre 0 et 1) minimisant d.
   */

  ecs_int_t     cpt_inter_loc ;

  ecs_real_t    det ;
  ecs_real_t    inv_det ;
  ecs_real_t    dist2 ;
  ecs_real_t    a, b, c, d, e, f ;
  ecs_real_t    s, t ;

  ecs_real_t    d_max_s ;
  ecs_real_t    d_max_t ;

  ecs_point_t v0 ;
  ecs_point_t v1 ;
  ecs_point_t v2 ;

  ecs_bool_t    int_s0_a1 = ECS_FALSE ;
  ecs_bool_t    int_s1_a1 = ECS_FALSE ;


  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  /* Nombre d'intersections */

  *nbr_inter = 0 ;

  /* ->                      */
  /* v0 = vecteur (P00, P01) */

  v0[0] = som_1_are_0[0] - som_0_are_0[0];
  v0[1] = som_1_are_0[1] - som_0_are_0[1];
  v0[2] = som_1_are_0[2] - som_0_are_0[2];


  /* ->                      */
  /* v1 = vecteur (P10, P11) */

  v1[0] = som_1_are_1[0] - som_0_are_1[0];
  v1[1] = som_1_are_1[1] - som_0_are_1[1];
  v1[2] = som_1_are_1[2] - som_0_are_1[2];


  /* ->                      */
  /* v2 = vecteur (P00, P10) */

  v2[0] = som_0_are_1[0] - som_0_are_0[0];
  v2[1] = som_0_are_1[1] - som_0_are_0[1];
  v2[2] = som_0_are_1[2] - som_0_are_0[2];

  /*
   *       ->   ->
   * a =   v0 . v0                  ->
   *                                v0        P01
   *       ->   ->       P00 x-------------->
   * b = - v0 . v1            \
   *                           \               P11
   *       ->   ->              \             ^
   * c =   v1 . v1               \           /
   *                           -> \         /
   *       ->   ->             v2  \       / ->
   * d = - v0 . v2                  \     /  v1
   *                                 \   /
   *       ->   ->                    \ /
   * e =   v1 . v2                     x
   *                                  P10
   *
   *       ->   ->
   * f =   v2 . v2
   */

  a =   ECS_LOC_PRODUIT_SCALAIRE(v0, v0) ;
  b = - ECS_LOC_PRODUIT_SCALAIRE(v0, v1) ;
  c =   ECS_LOC_PRODUIT_SCALAIRE(v1, v1) ;
  d = - ECS_LOC_PRODUIT_SCALAIRE(v0, v2) ;
  e =   ECS_LOC_PRODUIT_SCALAIRE(v1, v2) ;
  f =   ECS_LOC_PRODUIT_SCALAIRE(v2, v2) ;

  assert(a > 0) ;
  assert(c > 0) ;

  assert(sqrt(a) * fraction_dist * 1.00001 >= d_max_s0_a0) ;
  assert(sqrt(a) * fraction_dist * 1.00001 >= d_max_s1_a0) ;
  assert(sqrt(c) * fraction_dist * 1.00001 >= d_max_s0_a1) ;
  assert(sqrt(c) * fraction_dist * 1.00001 >= d_max_s1_a1) ;


  /*-----------------------------------------------------------------
   * Première partie : une intersection de type "arêtes parallèles",
   * "extrémité-extrémité" ou "extrémité-intérieur" est-elle possible
   * sous la tolérance indiquée ?
   *-----------------------------------------------------------------*/


  /* Distance du sommet 0 de l'arête 0 à l'arête 1 */

  t = -e / c ;

  if (t > -fraction_dist && t < 1.0 + fraction_dist) { /* pré-filtrage */

    dist2 = t*(c*t + 2.0*e) + f ;

    assert(dist2 >= -1.e-8) ;


    if (dist2 < (d_max_s0_a0 * d_max_s0_a0) ) {

      /* On est sous la tolérance du sommet 0 */

      d_max_t = t*d_max_s0_a1 + (1.0-t)*d_max_s1_a1 ;

      if (dist2 < (d_max_t * d_max_t)) {

        /* On a intersection avec le sommet 0 ; on vérifie s'il s'agit
           d'une intersection de type "sommet-sommet" */

        cpt_inter_loc = 0 ;

        if (t <= fraction_dist + 0.01) {

          dist2 = f ;

          if (   dist2 <= (d_max_s0_a0 * d_max_s0_a0)
              && dist2 <= (d_max_s0_a1 * d_max_s0_a1)) {

            /* Intersection s = t = 0.0 sous la tolérance */

            int_s0_a1 = ECS_TRUE ;

            s_inter_0[*nbr_inter] = -0.01 ;
            t_inter_1[*nbr_inter] = -0.01 ;

            *nbr_inter += 1 ;

            assert(*nbr_inter <= 2) ;

            cpt_inter_loc = 1 ;

          }

        }
        else if (t >= 0.99 - fraction_dist) {

          dist2 = (c + 2.0*e) + f ;

          if (   dist2 <= (d_max_s0_a0 * d_max_s0_a0)
              && dist2 <= (d_max_s1_a1 * d_max_s1_a1)) {

            /* Intersection s = 0.0, t = 1.0 sous la tolérance */

            int_s1_a1 = ECS_TRUE ;

            s_inter_0[*nbr_inter] = -0.01 ;
            t_inter_1[*nbr_inter] =  1.01 ;

            *nbr_inter += 1 ;

            assert(*nbr_inter <= 2) ;

            cpt_inter_loc = 1 ;

          }

        }

        if (cpt_inter_loc == 0 && 0.0 < t && t < 1.0) {

          /* On a une vraie intersection en t */ ;

          s_inter_0[*nbr_inter] = -0.01 ;
          t_inter_1[*nbr_inter] = t ;

          *nbr_inter += 1 ;
          assert(*nbr_inter <= 2) ;

        }

      }

    }

  }

  /* Distance du sommet 1 de l'arête 0 à l'arête 1 */

  t = -(b + e) / c ; /* t = (v1.(v0-v2)) / c = (-b -e) /c */

  if (t > -fraction_dist && t < 1.0 + fraction_dist) { /* pré-filtrage */

    dist2 = ECS_ABS((a + 2.0*(b*t + d)) + t*(c*t + 2.0*e) + f) ;

    assert(dist2 >= -1.e-8) ;

    if (dist2 < (d_max_s1_a0 * d_max_s1_a0) ) {

      /* On est sous la tolérance du sommet 1 */

      d_max_t = t*d_max_s0_a1 + (1.0-t)*d_max_s1_a1 ;

      if (dist2 < (d_max_t * d_max_t)) {

        /* On a intersection avec le sommet 1 ; on vérifie s'il s'agit
           d'une intersection de type "sommet-sommet" */

        cpt_inter_loc = 0 ;

        if (t <= fraction_dist + 0.01) {

          dist2 = ECS_ABS(a + 2.0*d + f) ;

          if (   dist2 <= (d_max_s1_a0 * d_max_s1_a0)
              && dist2 <= (d_max_s0_a1 * d_max_s0_a1)) {

            /* Intersection s = 1.0, t = 0.0 sous la tolérance */

            int_s0_a1 = ECS_TRUE ;

            s_inter_0[*nbr_inter] =  1.01 ;
            t_inter_1[*nbr_inter] = -0.01 ;

            *nbr_inter += 1 ;

            cpt_inter_loc = 1 ;

          }

        }
        else if (t >= 0.99 - fraction_dist) {

          dist2 = ECS_ABS((a + 2.0*(b + d)) + (c + 2.0*e) + f) ;

          if (   dist2 <= (d_max_s1_a0 * d_max_s1_a0)
              && dist2 <= (d_max_s1_a1 * d_max_s1_a1)) {

            /* Intersection s = t = 1.0 sous la tolérance */

            int_s1_a1 = ECS_TRUE ;

            s_inter_0[*nbr_inter] =  1.01 ;
            t_inter_1[*nbr_inter] =  1.01 ;

            *nbr_inter += 1 ;

            cpt_inter_loc = 1 ;

          }

        }

        if (cpt_inter_loc == 0 && 0.0 < t && t < 1.0) {

          /* On a une vraie intersection en t */

          s_inter_0[*nbr_inter] = 1.01 ;
          t_inter_1[*nbr_inter] = t ;

          *nbr_inter += 1 ;
          assert(*nbr_inter <= 2) ;

        }

      }

    }

  }

  /* Si l'on a déjà deux intersections aux sommets, inutile de continuer */

  if (*nbr_inter == 2)

    return ;


  /*
    Si l'on n'a pas encore d'intersection, on teste pour les sommets
    de l'arête 1
  */

  if (int_s0_a1 == ECS_FALSE) {

    /* Distance du sommet 0 de l'arête 1 à l'arête 0 */

    s = -d / a ;

    if (s > 0.0 && s < 1.0) { /* pré-filtrage ; s = 0 et 1 déjà testés */

      dist2 = ECS_ABS(s*(a*s + 2.0*d) + f) ;

      if (dist2 < (d_max_s0_a1 * d_max_s0_a1) ) {

        /* On est sous la tolérance du sommet 0 */

        d_max_s = s*d_max_s0_a0 + (1.0-s)*d_max_s1_a0 ;

        if (dist2 < (d_max_s * d_max_s)) {

          /* On a intersection avec le sommet 0 */

          assert (0 < s && s < 1) ;

          s_inter_0[*nbr_inter] = s ;
          t_inter_1[*nbr_inter] = -0.01 ;

          *nbr_inter += 1 ;
          assert(*nbr_inter <= 2) ;

          if (*nbr_inter == 2)
            return ;

        }

      }

    }

  }

  if (int_s1_a1 == ECS_FALSE) {

    /* Distance du sommet 1 de l'arête 1 à l'arête 0 */

    s = -(b + d) / a ; /* s = (v0.(v1+v2)) / c = (-b -d) / a */

    if (s > 0.0 && s < 1.0) { /* pré-filtrage ; s = 0 et 1 déjà testés */

      dist2 = ECS_ABS(s*(a*s + 2.0*(b + d)) + (c + 2.0*e) + f) ;

      if (dist2 < (d_max_s1_a1 * d_max_s1_a1) ) {

        /* On est sous la tolérance du sommet 1 */

        d_max_s = s*d_max_s0_a0 + (1.0-s)*d_max_s1_a0 ;

        if (dist2 < (d_max_s * d_max_s)) {

          /* On a intersection avec le sommet 1 */

          assert (0 < s && s < 1) ;

          s_inter_0[*nbr_inter] = s ;
          t_inter_1[*nbr_inter] = 1.01 ;

          *nbr_inter += 1 ;
          assert(*nbr_inter <= 2) ;

          if (*nbr_inter == 2)
            return ;

        }

      }

    }

  }

  /* Si l'on au moins une intersection aux sommets, inutile de continuer */

  if (*nbr_inter > 0)

    return ;


  /*------------------------------------------------------------------
   * Seconde partie : aucune intersection de type "arêtes parallèles",
   * "extrémité-extrémité", ou "extrémité-intérieur" n'a été trouvée.
   * La seule possibilité restante est le cas où 0 < s < 1, 0 < t < 1.
   *------------------------------------------------------------------*/


  /*
   * Les arêtes sont parallèlles si a.c - b^2 = 0.
   * Cependant, ce terme est proportionnel à la quatrième puissance de la
   * longueur des arêtes. On préferera se fier à un angle, mais on saura
   * par la suite que si les arêtes ne sont pas parallèlles, alors
   * a.c - b^2 > 0.
   */

  /*
   * cos (theta) = v0.v1 / ||v0||.||v1||
   * cos^2 (theta) = (v0.v1)^2 / ((v0.v0).(v1.v1)) = (-b)^2 / (a.c)
   * On compare donc eps_paral^2 a (1 - b^2 / (a.c))
   */


  det = ECS_ABS(a * c - b * b);

  if (det <  0.001 * 0.001 * a * c)

    return ;


  /* Les arêtes ne sont pas parallèlles */


  s = b * e - c * d ;
  t = b * d - a * e ;

  if ( s >= 0. ) {

    if ( s <= det ) {

      if ( t >= 0. ) {

        if ( t <= det )  { /* région centrale (intérieure) */

          /* La distance minimale relie 2 points de la région centrale */
          inv_det = 1.0 / det ;
          s *= inv_det ;
          t *= inv_det ;
          dist2 = ECS_ABS(s*(a*s + 2.0*(b*t + d)) + t*(c*t + 2.0*e) + f) ;

          d_max_s = s*d_max_s0_a0 + (1.0-s)*d_max_s1_a0 ;
          d_max_t = t*d_max_s0_a1 + (1.0-t)*d_max_s1_a1 ;

          if (   dist2 <= (d_max_s * d_max_s)
              && dist2 <= (d_max_t * d_max_t)) {

            assert (0 <= s && s <= 1) ;
            assert (0 <= t && t <= 1) ;

            s_inter_0[*nbr_inter] = s ;
            t_inter_1[*nbr_inter] = t ;

            *nbr_inter += 1 ;
            assert(*nbr_inter <= 2) ;

          }

        }

      }

    }

  }

  return ;

}


/*----------------------------------------------------------------------------
 *  Fonction qui renvoie le sommet associé à une intersection
 *  (en le créant si nécessaire)
 *----------------------------------------------------------------------------*/


static ecs_int_t  ecs_loc_vec_def_rc__som_inter
(
       ecs_vec_int_t           *const vec_def_are,
       ecs_vec_real_t          *const vec_def_som,
       size_t                  *const nbr_som_new,
       ecs_tab_real_t          *const dist_max_som,
 const ecs_rc_inter_are_are_t  *const inter_are
)
{

  ecs_int_t  ind_som ;
  ecs_int_t  isom_0 ;
  ecs_int_t  isom_1 ;
  ecs_int_t  iloc ;

  ecs_real_t s = inter_are->s_inter ;


  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  /* Indices des sommets aux extrémités de l'arête */

  isom_0 = vec_def_are->val_tab[vec_def_are->pos_tab[inter_are->are] - 1] - 1 ;
  isom_1 = vec_def_are->val_tab[vec_def_are->pos_tab[inter_are->are]    ] - 1 ;


  /* intersection avec sommet existant ou non */

  if (inter_are->s_inter <= 0.0)

    ind_som = isom_0 ;

  else if (inter_are->s_inter >= 1.0)

    ind_som = isom_1 ;

  else {

    ind_som = *nbr_som_new ;

    assert (*nbr_som_new < vec_def_som->pos_nbr - 1) ;


    for (iloc = 0 ; iloc < ECS_DIM_3 ; iloc++) {

      vec_def_som->val_tab[(ind_som * ECS_DIM_3) + iloc]
        =   (1.0 - s) * vec_def_som->val_tab[(isom_0 * ECS_DIM_3) + iloc]
          +        s  * vec_def_som->val_tab[(isom_1 * ECS_DIM_3) + iloc] ;

    }

    dist_max_som->val[ind_som]
      = (1.0 - s) * dist_max_som->val[isom_0] + s * dist_max_som->val[isom_0] ;

    *nbr_som_new += 1 ;

  }

  return ind_som ;

}


/*----------------------------------------------------------------------------
 *    Fonction de descente d'un arbre binaire pour le tri lexicographique
 *  d'un tableau d'intersections d'arêtes.
 *----------------------------------------------------------------------------*/

static void ecs_loc_rc__desc_arbre_inter
(
       ecs_int_t                        ltree,       /* --> Niveau arbre      */
 const ecs_int_t                        ntree,       /* --> Taille arbre      */
       ecs_rc_inter_are_are_t    *const table_are,   /* <-> Tableau à trier   */
       ecs_rc_inter_are_inter_t  *const table_inter  /* <-> Tableau associé   */
)
{
  ecs_int_t               ktree, l_save;
  ecs_rc_inter_are_are_t  i_save, p1, p2;
  int                     isup;


  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  i_save = table_are[ltree];
  l_save = ltree;

  while (ltree <= (ntree/2)) {

    ktree = (2*ltree)+1;

    if (ktree < ntree - 1) {

      /* isup vrai si table_are[ktree+1] > table_are[ktree] */
      p1 = table_are[ktree + 1];
      p2 = table_are[ktree];
      isup = (p1.are > p2.are) ? 1 : 0;
      if (p1.are == p2.are)
        isup = (p1.s_inter > p2.s_inter) ? 1 : 0;

      /* si isup vrai, on incremente ktree */
      if (isup) ktree++;
    }

    if (ktree >= ntree) break;

    /* isup faux si table_are[ltree] (initial) < table_are[ktree] */
    p1 = i_save;
    p2 = table_are[ktree];
    isup = (p1.are < p2.are) ? 0 : 1;
    if (p1.are == p2.are)
      isup = (p1.s_inter < p2.s_inter) ? 0 : 1;

    /* si isup vrai */
    if (isup) break;

    table_are[ltree] = table_are[ktree];
    if ((table_inter[(table_are[ltree]).i_inter]).ind_0 == ktree)
      (table_inter[(table_are[ltree]).i_inter]).ind_0 = ltree ;
    else
      /* if ((table_inter[(table_are[ltree]).i_inter]).ind_1 == ktree) */
      (table_inter[(table_are[ltree]).i_inter]).ind_1 = ltree ;

    ltree = ktree;
  }

  table_are[ltree] = i_save;
  if ((table_inter[(table_are[ltree]).i_inter]).ind_0 == l_save)
    (table_inter[(table_are[ltree]).i_inter]).ind_0 = ltree ;
  else
    /* if ((table_inter[(table_are[ltree]).i_inter]).ind_1 == l_save) */
    (table_inter[(table_are[ltree]).i_inter]).ind_1 = ltree ;

}


/*----------------------------------------------------------------------------
 *    Fonction de tri lexicographique d'une liste d'intersections d'arêtes
 *----------------------------------------------------------------------------*/

static void ecs_loc_vec_def_rc__trie_inter
(
 ecs_rc_lis_inter_are_t *const liste_inter
)
{

  ecs_rc_inter_are_are_t  i_save;
  ecs_int_t               i;
  ecs_int_t               ind_0;
  ecs_int_t               ind_1;


  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/

  ECS_CHRONO("Joining: intersections sort start   ") ;


  /* creation de l'arbre binaire vec_renum->val_tab[liste_inter->nbr] */

  for (i = (liste_inter->nbr) - 1 ; i >= 0 ; i--)

    ecs_loc_rc__desc_arbre_inter (i,
                                  liste_inter->nbr*2,
                                  liste_inter->are,
                                  liste_inter->inter) ;

  /* tri de l'arbre binaire */

  for (i = (liste_inter->nbr*2) - 1 ; i > 0 ; i--) {

    i_save              = liste_inter->are[0] ;
    liste_inter->are[0] = liste_inter->are[i] ;
    liste_inter->are[i] = i_save ;

    if ((liste_inter->inter[(liste_inter->are[0]).i_inter]).ind_0 == i)
      (liste_inter->inter[(liste_inter->are[0]).i_inter]).ind_0 = 0 ;
    else
      /* if ((liste_inter->inter[(liste_inter->are[0]).i_inter]).ind_1 == i) */
      (liste_inter->inter[(liste_inter->are[0]).i_inter]).ind_1 = 0 ;

    if ((liste_inter->inter[i_save.i_inter]).ind_0 == 0)
      (liste_inter->inter[i_save.i_inter]).ind_0 = i ;
    else
      /* if ((liste_inter->inter[i_save.i_inter]).ind_1 == 0) */
      (liste_inter->inter[i_save.i_inter]).ind_1 = i ;

    ecs_loc_rc__desc_arbre_inter (0,
                                  i,
                                  liste_inter->are,
                                  liste_inter->inter) ;
  }


  for (i = 0 ; i < (ecs_int_t)(liste_inter->nbr) ; i++) {

    ind_0 = (liste_inter->inter[i]).ind_0 ;
    ind_1 = (liste_inter->inter[i]).ind_1 ;

    if (ind_0 > ind_1) {

      (liste_inter->inter[i]).ind_0 = ind_1 ;
      (liste_inter->inter[i]).ind_1 = ind_0 ;

    }

  }


#if 0 && defined(DEBUG) && !defined(NDEBUG)
  {
    printf("  intersection triée |   arête   |   inter   |    s_inter\n") ;
    for (i = 0 ; i < liste_inter->nbr*2 ; i++)
      printf("           %9d | %9d | %9d | %12.5e\n",
             i,
             (liste_inter->are[i]).are,
             (liste_inter->are[i]).i_inter,
             (liste_inter->are[i]).s_inter) ;
    for (i = 0 ; i < liste_inter->nbr-1 ; i++) {
      assert (   (liste_inter->are[i]).are
              <= (liste_inter->are[i+1]).are) ;
      assert (    (liste_inter->are[i]).are
               != (liste_inter->are[i+1]).are ||
              (liste_inter->are[i]).s_inter
              <= (liste_inter->are[i+1]).s_inter) ;
    }
    printf("\n  liens intersection |  indice_0 |  indice_1\n") ;
    for (i = 0 ; i < liste_inter->nbr ; i++)
      printf("           %9d | %9d | %9d\n",
             i,
             (liste_inter->inter[i]).ind_0,
             (liste_inter->inter[i]).ind_1) ;
  }
#endif

  ECS_CHRONO("Joining: intersections sort done    ") ;

}


#if 0 && defined(DEBUG) && !defined(NDEBUG)

/*----------------------------------------------------------------------------
 *  Fonction qui imprime une structure ecs_vec_int_t
 *----------------------------------------------------------------------------*/

static void  ecs_dbg_vec_def_rc__imp_vec_int
(
 const char           *const titre,
 const char           *const nom,
 const ecs_vec_int_t  *const vec_int_dbg
)
{

  ecs_int_t  ient ;
  ecs_int_t  ient_new ;
  size_t     nbr_old_new ;
  ecs_int_t  num_ent ;
  size_t     pos_ent_old ;



  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  {

    printf("\n\n%s\n", titre) ;

    printf("\n%s :\n", nom) ;
    for (ient = 0 ; ient < vec_int_dbg->pos_nbr - 1 ; ient++) {
      pos_ent_old = vec_int_dbg->pos_tab[ient] - 1 ;
      nbr_old_new
        = vec_int_dbg->pos_tab[ient + 1] - 1
        - pos_ent_old ;
      printf("  %d -> ", ient + 1) ;
      for (ient_new = 0 ; ient_new < nbr_old_new ; ient_new++) {
        num_ent = vec_int_dbg->val_tab[pos_ent_old + ient_new] ;
        printf(" %d", num_ent) ;
      }
      printf("\n") ;
    }
    fflush(stdout);

  }

}

#endif


#if 0 && defined(DEBUG) && !defined(NDEBUG)

/*----------------------------------------------------------------------------
 *  Fonction qui imprime une structure ecs_vec_real_t
 *----------------------------------------------------------------------------*/

static void  ecs_dbg_vec_def_rc__imp_vec_real
(
 const char            *const titre,
 const char            *const nom,
 const ecs_vec_real_t  *const vec_real_dbg
)
{

  ecs_int_t  ient ;
  size_t     ipos ;
  ecs_int_t  icoo ;

  ecs_real_t coo ;



  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  {

    printf("\n\n%s\n", titre) ;

    printf("\n%s :\n", nom) ;
    for (ient = 0 ; ient < vec_real_dbg->pos_nbr - 1 ; ient++) {
      ipos = vec_real_dbg->pos_pas * ient ;
      printf("  %d -> ", ient + 1) ;
      for (icoo = 0 ; icoo < vec_real_dbg->pos_pas ; icoo++) {
        coo = vec_real_dbg->val_tab[ipos + icoo] ;
        printf(" %f", coo) ;
      }
      printf("\n") ;
    }
    fflush(stdout);

  }

}

#endif

Generated by  Doxygen 1.6.0   Back to index