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

ecs_vec_def.c

/*============================================================================
 *  Définitions des fonctions
 *   associées aux structures `ecs_vec_int_t' et `ecs_vec_real_t' décrivant
 *   les vecteurs indexes entier et réel
 *   et propres aux vecteurs indexés
 *      liés aux champs principaux de type "définition"
 *============================================================================*/

/*
  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
*/


/*============================================================================
 *                                 Visibilité
 *============================================================================*/

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

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

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


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

#include "ecs_chrono.h"

#include "ecs_def.h"
#include "ecs_tab.h"

#include "ecs_elt_typ_liste.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_int_tri.h"
#include "ecs_vec_real_tri.h"


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

#include "ecs_vec_def.h"


/*----------------------------------------------------------------------------
 *  Fichiers `include' privés   du  paquetage courant
 *----------------------------------------------------------------------------*/

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


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

#if !defined(FLT_MAX)
#define FLT_MAX HUGE_VAL
#endif

enum {
  X,
  Y,
  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]   )


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


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


#define ECS_LOC_DETERMINANT(vect1, vect2, vect3) ( \
   ((vect1[Y] * vect2[Z] - vect2[Y] * vect1[Z]) * vect3[X]) \
 + ((vect2[X] * vect1[Z] - vect1[X] * vect2[Z]) * vect3[Y]) \
 + ((vect1[X] * vect2[Y] - vect2[X] * vect1[Y]) * vect3[Z]) )


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

/*----------------------------------------------------------------------------
 *  Fonction qui réalise l'ordonnancement des éléments
 *   en fonction de leur définition
 *----------------------------------------------------------------------------*/

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


/*----------------------------------------------------------------------------
 *  Fonction réalisant la liste compactée des éléments
 *         à partir de la liste ordonnée  des éléments
 *
 *  Le compactage consiste à supprimer les éléments géométriquement identiques
 *----------------------------------------------------------------------------*/

static ecs_vec_real_t * ecs_loc_vec_real_def__compacte
(                                                /* <-- Structure compactée   */
 const ecs_vec_real_t *const this_vec_real_ord,  /* --> Structure à compacter */
       ecs_tab_int_t  *const vect_transf         /* <-> Vecteur de transf.    */
) ;


/*----------------------------------------------------------------------------
 *  Fonction qui projette une face sur un plan parallèle à la face.
 *  Ce plan est ensuite assimilé au plan (Oxy).
 *----------------------------------------------------------------------------*/

static void ecs_loc_vec_def__plan_fac
(
 const ecs_int_t           nbr_som_fac,
       ecs_point_t  *const coo_som_fac
) ;


/*----------------------------------------------------------------------------
 * Fonction qui vérifie si un sommet d'un polygone (projeté en 2D) est
 *  convexe.
 *----------------------------------------------------------------------------*/

static ecs_bool_t ecs_loc_vec_def__test_convexe
(
 const ecs_int_t           som_prec,
 const ecs_int_t           som_cour,
 const ecs_int_t           som_suiv,
       ecs_point_t  *const coo_som_fac
) ;


/*----------------------------------------------------------------------------
 * Fonction qui vérifie si un sommet d'un polygone (projeté en 2D) est
 *  une "oreille".
 *----------------------------------------------------------------------------*/

static ecs_bool_t ecs_loc_vec_def__test_oreille
(
 const ecs_int_t           isom,
 const ecs_int_t           nbr_som_fac,
 const ecs_int_t    *const liste_prec,
 const ecs_int_t    *const liste_suiv,
 const ecs_bool_t   *const concave,
       ecs_point_t  *const coo_som_fac,
 const ecs_real_t          epsilon
) ;


/*----------------------------------------------------------------------------
 *  Fonction transformant une triangulation quelconque en triangulation
 *  de Delaunay grâce à un algorithme d'échange d'arêtes (méthode flip).
 *----------------------------------------------------------------------------*/

static void ecs_loc_vec_def__delaunay_flip
(
 const ecs_int_t           nbr_som_fac,
       ecs_int_t    *const liste_som_tria,
       ecs_int_t    *const liste_are_def,
       ecs_int_t    *const liste_voisinage_are,
       ecs_bool_t   *const liste_are_loc_delaunay,
       ecs_point_t  *const coo_som_fac
) ;


/*----------------------------------------------------------------------------
 *  Fonction renvoyant un booléen selon si l'arête est localement de
 *  Delaunay ou non
 *----------------------------------------------------------------------------*/

static ecs_bool_t ecs_loc_vec_def__are_loc_delaunay
(
 const ecs_int_t          isom_are_0,
 const ecs_int_t          isom_are_1,
 const ecs_int_t          isom_flip_0,
 const ecs_int_t          isom_flip_1,
       ecs_point_t *const coo_som_fac
) ;


/*----------------------------------------------------------------------------
 *  Fonction triant les indices des sommets définissant un triangle
 *----------------------------------------------------------------------------*/

static void ecs_loc_vec_def__trie_delaunay
(
       ecs_int_t   *const liste_som_tria
) ;


/*----------------------------------------------------------------------------
 *  Fonction réalisant le découpage des faces polygonales en triangles
 *----------------------------------------------------------------------------*/

static ecs_int_t ecs_loc_vec_def__dec_fac_tria
(
 const ecs_int_t            nbr_som_fac,
       ecs_int_t     *const liste_som_tria,
       ecs_int_t     *const liste_prec,
       ecs_int_t     *const liste_suiv,
       ecs_int_t     *const liste_are_def,
       ecs_int_t     *const liste_voisinage_are,
       ecs_bool_t    *const liste_are_loc_delaunay,
       ecs_bool_t    *const concave,
       ecs_point_t   *const coord_som
) ;


/*----------------------------------------------------------------------------
 *  Fonction qui transforme un tableau d'équivalence en liste chaînée simple
 *----------------------------------------------------------------------------*/

static void ecs_loc_vec_def__transf_equiv
(
 const ecs_vec_real_t    *const vec_def_som,    /* <-> Déf. sommets           */
       ecs_tab_int_t     *const tab_equiv_som   /* <-> Fusions de sommets     */
) ;


/*----------------------------------------------------------------------------
 *  Correction si nécessaire de l'orientation d'un quadrangle en connectivité
 *   nodale ; renvoie -1 si l'on ne parvient pas à corriger l'orientation
 *   ou que l'on demande une simple vérification (i.e. correc = ECS_FALSE),
 *   0 si l'orientation initiale est bonne, et 1 en cas de correction de
 *   l'orientation par une permutation de la connectivité locale.
 *----------------------------------------------------------------------------*/

static ecs_int_t ecs_loc_vec_def__orient_quad
(
 const ecs_real_t  coord[],
       ecs_int_t   connect[],
 const ecs_bool_t  correc
) ;


/*----------------------------------------------------------------------------
 *  Correction si nécessaire de l'orientation d'un tétraèdre en connectivité
 *   nodale ; renvoie -1 si l'on ne parvient pas à corriger l'orientation
 *   ou que l'on demande une simple vérification (i.e. correc = ECS_FALSE),
 *   0 si l'orientation initiale est bonne, et 1 en cas de correction de
 *   l'orientation par une permutation de la connectivité locale.
 *----------------------------------------------------------------------------*/

static ecs_int_t ecs_loc_vec_def__orient_tetra
(
 const ecs_real_t  coord[],
       ecs_int_t   connect[],
 const ecs_bool_t  correc
) ;


/*----------------------------------------------------------------------------
 *  Correction si nécessaire de l'orientation d'une pyramide en connectivité
 *   nodale ; renvoie -1 si l'on ne parvient pas à corriger l'orientation
 *   ou que l'on demande une simple vérification (i.e. correc = ECS_FALSE),
 *   0 si l'orientation initiale est bonne, et 1 en cas de correction de
 *   l'orientation par une permutation de la connectivité locale.
 *
 *  Tous les cas ne sont pas détectés ou traités : on suppose que le
 *   sommet est toujours en position 5, mais que la base peut être
 *   parcourue dans l'ordre 1 4 3 2, 1 2 4 3, ou 1 3 4 2 au lieu de 1 2 3 4,
 *   dans quel cas on la réordonne.
 *----------------------------------------------------------------------------*/

static ecs_int_t ecs_loc_vec_def__orient_pyram
(
 const ecs_real_t  coord[],
       ecs_int_t   connect[],
 const ecs_bool_t  correc
) ;


/*----------------------------------------------------------------------------
 *  Correction si nécessaire de l'orientation d'un prisme en connectivité
 *   nodale ; renvoie -1 si l'on ne parvient pas à corriger l'orientation
 *   ou que l'on demande une simple vérification (i.e. correc = ECS_FALSE),
 *   0 si l'orientation initiale est bonne, et 1 en cas de correction de
 *   l'orientation par une permutation de la connectivité locale.
 *
 *  Tous les cas ne sont pas détectés ou traités : on suppose que les
 *   bases peuvent être parcourues dans l'ordre 1 3 2 et 4 6 5 au lieu
 *   de 1 2 3 et 4 5 6, dans quel cas on les réordonne.
 *----------------------------------------------------------------------------*/

static ecs_int_t ecs_loc_vec_def__orient_prism
(
 const ecs_real_t  coord[],
       ecs_int_t   connect[],
 const ecs_bool_t  correc
) ;


/*----------------------------------------------------------------------------
 *  Correction si nécessaire de l'orientation d'un hexaèdre en connectivité
 *   nodale ; renvoie -1 si l'on ne parvient pas à corriger l'orientation
 *   ou que l'on demande une simple vérification (i.e. correc = ECS_FALSE),
 *   0 si l'orientation initiale est bonne, et 1 en cas de correction de
 *   l'orientation par une permutation de la connectivité locale.
 *
 *  Tous les cas ne sont pas détectés ou traités : on suppose que les
 *   bases peuvent être parcourues dans l'ordre soit 1 4 3 2 et 5 8 7 6,
 *   soit 1 2 4 3 et 5 6 8 7, soit 1 3 4 2 et 5 7 8 6 au lieu
 *   de 1 2 3 4 et 5 6 7 8, dans quel cas on les réordonne.
 *----------------------------------------------------------------------------*/

static ecs_int_t ecs_loc_vec_def__orient_hexa
(
 const ecs_real_t  coord[],
       ecs_int_t   connect[],
 const ecs_bool_t  correc
) ;


/*----------------------------------------------------------------------------
 *  Correction si nécessaire de l'orientation d'un polyedre en connectivité
 *   nodale ; renvoie -1 si l'on ne parvient pas à corriger l'orientation
 *   ou que l'on demande une simple vérification (i.e. correc = ECS_FALSE),
 *   0 si l'orientation initiale est bonne, et 1 en cas de correction de
 *   l'orientation par une permutation de la connectivité locale.
 *
 *  On suppose que toutes les faces polygonales sont bien définies, mais
 *   que certaines faces peuvent être définies avec une normale intérieure
 *   à la cellule et non pas extérieure. On suppose que la non-convexité
 *   du polyèdre est limitée.
 *----------------------------------------------------------------------------*/

static ecs_int_t ecs_loc_vec_def__orient_polyedre
(
 const ecs_real_t      coord[],
       ecs_int_t       connect[],
       ecs_int_t       size,
       ecs_bool_t      correc,
       ecs_tab_int_t  *face_index,
       ecs_tab_int_t  *face_marker,
       ecs_tab_int_t  *edges
);


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

/*----------------------------------------------------------------------------
 *  Fonction réalisant la décomposition des cellules en faces
 *   sur les vecteurs `ecs_vec_int_t' associes aux champs principaux
 *----------------------------------------------------------------------------*/

void ecs_vec_def__decompose_cel
(
 ecs_vec_int_t **vec_def_fac,        /* <-- Définitions des faces             */
 ecs_vec_int_t **vec_cel_def_fac,    /* <-- Définitions des cellules par face */
 ecs_vec_int_t  *vec_def_cel         /* --> Définitions des cellules          */
)
{

  size_t      nbr_cel ;
  size_t      nbr_def ;
  size_t      nbr_fac ;
  size_t      nbr_val_cel ;
  size_t      nbr_val_fac ;
  size_t      ind_pos_cel ;
  size_t      ind_pos_loc ;
  size_t      ind_pos_sui ;
  size_t      nbr_pos_loc ;
  ecs_int_t   num_def ;
  ecs_int_t   marqueur_fin ;

  size_t      isom ;
  size_t      icel ;
  int         idef ;
  size_t      cpt_fac ;
  int         ifac ;

  ecs_int_t   typ_geo_cel ;       /* Type géométrique élément en cours */

  ecs_int_t   typ_geo_base[9] = {ECS_ELT_TYP_NUL,
                                 ECS_ELT_TYP_NUL,
                                 ECS_ELT_TYP_NUL,
                                 ECS_ELT_TYP_NUL,
                                 ECS_ELT_TYP_CEL_TETRA,
                                 ECS_ELT_TYP_CEL_PYRAM,
                                 ECS_ELT_TYP_CEL_PRISM,
                                 ECS_ELT_TYP_NUL,
                                 ECS_ELT_TYP_CEL_HEXA};

  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  /*=================*/
  /* Initialisations */
  /*=================*/


  nbr_cel   = ecs_vec_int__ret_pos_nbr(vec_def_cel) - 1 ;


  /* Boucle de comptage pour l'allocation des sous-éléments */
  /*--------------------------------------------------------*/

  nbr_fac     = 0 ;
  nbr_val_fac = 0 ;
  nbr_def       = 0 ;

  for (icel = 0 ; icel < nbr_cel ; icel++ ) {

    ind_pos_cel = vec_def_cel->pos_tab[icel] - 1 ;

    nbr_val_cel = vec_def_cel->pos_tab[icel + 1] - 1 - ind_pos_cel ;

    /* Traitement des cellules "classiques" */
    /*--------------------------------------*/

    if (nbr_val_cel < 9) {

      typ_geo_cel = typ_geo_base[nbr_val_cel] ;

      /* Boucle sur les sous-éléments définissant la cellulle */

      for (ifac = 0 ;
           ifac < ecs_fic_elt_typ_liste_c[typ_geo_cel].nbr_sous_elt;
           ifac++ ) {

        const ecs_sous_elt_t  * sous_elt
          = &(ecs_fic_elt_typ_liste_c[typ_geo_cel].sous_elt[ifac]) ;

        nbr_fac++ ;

        for (idef = 0 ;
             idef < (ecs_fic_elt_typ_liste_c[sous_elt->elt_typ]).nbr_som ;
             idef++) ;

        nbr_val_fac += idef ;

      }

    }

    /* Traitement des éléments de type "polyèdre" */
    /*--------------------------------------------*/

    else {

      ind_pos_sui = vec_def_cel->pos_tab[icel + 1] - 1 ;
      nbr_pos_loc = ind_pos_sui - ind_pos_cel ;

      /* Convention : définition nodale cellule->sommets avec numéros de
         premiers sommets répétés en fin de liste pour marquer la fin
         de chaque face */

      marqueur_fin = -1 ;

      for (isom = ind_pos_cel ; isom < ind_pos_sui ; isom++) {

        if (vec_def_cel->val_tab[isom] != marqueur_fin) {
          nbr_val_fac += 1 ;
          if (marqueur_fin == -1)
            marqueur_fin = vec_def_cel->val_tab[isom] ;
        }
        else {
          marqueur_fin = -1 ;
          nbr_fac += 1 ;
        }

      }

    }

  } /* Fin de la boucle de comptage sur les cellules */


  /* Allocation et initialisation pour les faces            */
  /*  des vecteurs `ecs_vec_int_t' associés aux définitions */
  /*--------------------------------------------------------*/

  *vec_def_fac     = ecs_vec_int__alloue(nbr_fac + 1, nbr_val_fac) ;
  *vec_cel_def_fac = ecs_vec_int__alloue(nbr_cel + 1, nbr_fac) ;

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


  /*=======================================*/
  /* Boucle sur les cellules à transformer */
  /*=======================================*/


  cpt_fac  = 0 ;
  nbr_def    = 0 ;

  for (icel = 0 ; icel < nbr_cel ; icel++ ) {

    ind_pos_cel = vec_def_cel->pos_tab[icel] - 1 ;

    nbr_val_cel = vec_def_cel->pos_tab[icel + 1] - 1 - ind_pos_cel ;

    /*--------------------------------------*/
    /* Traitement des éléments "classiques" */
    /*--------------------------------------*/

    if (nbr_val_cel < 9) {

      typ_geo_cel = typ_geo_base[nbr_val_cel] ;

      /* Boucle sur les faces définissant la cellulle */
      /*==============================================*/

      for (ifac = 0 ;
           ifac < ecs_fic_elt_typ_liste_c[typ_geo_cel].nbr_sous_elt;
           ifac++ ) {

        const ecs_sous_elt_t  * sous_elt
          = &(ecs_fic_elt_typ_liste_c[typ_geo_cel].sous_elt[ifac]) ;


        /* Boucle sur les sommets définissant la face */
        /*--------------------------------------------*/

        for (idef = 0 ;
             idef < (ecs_fic_elt_typ_liste_c[sous_elt->elt_typ]).nbr_som ;
             idef++) {

          /* Définition de la face en fonction des sommets */

          num_def = sous_elt->som[idef] ;

          (*vec_def_fac)->val_tab[nbr_def++]
            = vec_def_cel->val_tab[ind_pos_cel + num_def - 1] ;

        }


        /* Position de la face dans sa définition en fonction des sommets */
        /*----------------------------------------------------------------*/

        (*vec_def_fac)->pos_tab[cpt_fac + 1]
          = (*vec_def_fac)->pos_tab[cpt_fac] + idef ;


        /* Détermination de la cellule en fonction des faces */
        /*---------------------------------------------------*/

        (*vec_cel_def_fac)->val_tab[cpt_fac] = cpt_fac + 1 ;

        cpt_fac++ ;

      } /* Fin de la boucle sur les faces d'une cellule */

    }


    /*--------------------------------------------*/
    /* Traitement des éléments de type "polyèdre" */
    /*--------------------------------------------*/

    else {

      /* Boucle sur les faces définissant le polyèdre */
      /*==============================================*/

      ifac = 0 ;
      marqueur_fin = -1 ;

      nbr_pos_loc = vec_def_cel->pos_tab[icel + 1] - 1 - ind_pos_cel ;

      for (ind_pos_loc = 0 ; ind_pos_loc < nbr_pos_loc ; ind_pos_loc++) {

        /* Définition de la face en fonction des sommets */

        if (vec_def_cel->val_tab[ind_pos_cel + ind_pos_loc] != marqueur_fin) {

          (*vec_def_fac)->val_tab[nbr_def++]
            = vec_def_cel->val_tab[ind_pos_cel + ind_pos_loc] ;

          if (marqueur_fin == -1)
            marqueur_fin = vec_def_cel->val_tab[ind_pos_cel + ind_pos_loc] ;

        }

        /* Position de la face dans sa définition en fonction des sommets */

        else {

          (*vec_def_fac)->pos_tab[cpt_fac + 1] = nbr_def + 1 ;

          marqueur_fin = -1 ;

          (*vec_cel_def_fac)->val_tab[cpt_fac] = cpt_fac + 1 ;

          /* Incrémentation du nombre de faces */

          ifac++ ;
          cpt_fac++ ;

        }

      } /* Fin de la boucle sur les faces du polyèdre */

    }


    /* A ce point, on a : ifac
       = ecs_fic_elt_typ_liste_c[vec_typ_geo_cel->val_tab[icel]].nbr_sous_elt
       pour des éléments classiques, et ifac est égal au nombre de faces pour
       des éléments polyédriques. */

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


  } /* Fin de la boucle sur les éléments */


  assert(ecs_vec_int__ret_val_nbr(*vec_def_fac) == nbr_val_fac) ;

}


/*----------------------------------------------------------------------------
 *  Fonction réalisant la décomposition des faces en arêtes
 *   sur les vecteurs `ecs_vec_int_t' associes aux champs principaux
 *----------------------------------------------------------------------------*/

void ecs_vec_def__decompose_fac
(
 ecs_vec_int_t **vec_def_are,        /* <-- Définitions des arêtes            */
 ecs_vec_int_t **vec_fac_def_are,    /* <-- Définitions des faces par arêtes  */
 ecs_vec_int_t  *vec_def_fac         /* --> Définitions des faces             */
)
{

  size_t      nbr_fac ;
  size_t      nbr_def ;
  size_t      nbr_are ;
  size_t      nbr_val_are ;
  size_t      ind_pos_fac ;
  size_t      ind_pos_loc ;
  size_t      nbr_pos_loc ;

  size_t      ifac ;
  size_t      cpt_are ;
  int         iare ;

  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  /*=================*/
  /* Initialisations */
  /*=================*/


  nbr_fac   = ecs_vec_int__ret_pos_nbr(vec_def_fac) - 1 ;

  /* Boucle de comptage pour l'allocation des sous-éléments */
  /*--------------------------------------------------------*/

  nbr_are     = 0 ;
  nbr_val_are = 0 ;
  nbr_def       = 0 ;

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

    nbr_pos_loc = vec_def_fac->pos_tab[ifac + 1] - vec_def_fac->pos_tab[ifac] ;

    nbr_are     += nbr_pos_loc ;
    nbr_val_are += nbr_pos_loc * 2 ;

  }


  /* Allocation et initialisation pour les arêtes,          */
  /*  des vecteurs `ecs_vec_int_t' associes aux définitions */
  /*--------------------------------------------------------*/

  *vec_def_are     = ecs_vec_int__alloue(nbr_are + 1, nbr_val_are) ;
  *vec_fac_def_are = ecs_vec_int__alloue(nbr_fac + 1, nbr_are) ;

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


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


  cpt_are  = 0 ;
  nbr_def    = 0 ;

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

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


    /* Boucle sur les arêtes définissant la face */
    /*===========================================*/

    nbr_pos_loc = vec_def_fac->pos_tab[ifac + 1] - 1 - ind_pos_fac ;

    iare = nbr_pos_loc ;

    for (ind_pos_loc = 0 ; ind_pos_loc < nbr_pos_loc ; ind_pos_loc++) {

      /* Définition de l'arête en fonction des sommets */

      (*vec_def_are)->val_tab[nbr_def++]
        = vec_def_fac->val_tab[ind_pos_fac + ((ind_pos_loc  )%nbr_pos_loc)] ;
      (*vec_def_are)->val_tab[nbr_def++]
        = vec_def_fac->val_tab[ind_pos_fac + ((ind_pos_loc+1)%nbr_pos_loc)] ;


      /* Position de l'arête dans sa définition en fonction des sommets */

      (*vec_def_are)->pos_tab[cpt_are + 1]
        = (*vec_def_are)->pos_tab[cpt_are] + 2 ;


      /* Détermination de la face en fonction des arêtes */

      (*vec_fac_def_are)->val_tab[cpt_are] = cpt_are + 1 ;

      cpt_are++ ;

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


    /* A ce point iare est égal au nombre d'arêtes de la face */

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


  } /* Fin de la boucle sur les faces */


  assert(ecs_vec_int__ret_val_nbr(*vec_def_are) == nbr_val_are) ;

}


/*----------------------------------------------------------------------------
 *  Fonction qui réalise la fusion des définitions des éléments
 *
 *  Les valeurs absolues sont triées si `bool_abs' est à `ECS_TRUE'
 *----------------------------------------------------------------------------*/

ecs_tab_int_t ecs_vec_def__fusionne
(
 ecs_vec_int_t * *const this_vec_def,
 ecs_tab_int_t   *      signe_elt
)
{

  size_t           cpt_sup_fin ;
  size_t           ind_cmp ;
  size_t           ind_inf ;
  size_t           ind_loc_sup ;
  size_t           ind_loc_cmp ;
  size_t           ind_pos ;
  size_t           ind_pos_loc ;
  size_t           ind_sup ;
  size_t           nbr_sup_ini ;
  size_t           nbr_inf ;
  ecs_int_t        num_inf ;
  ecs_int_t        num_inf_loc ;
  ecs_int_t        num_inf_min ;
  ecs_int_t        num_inf_min_cmp ;
  size_t           pos_cmp ;
  size_t           pos_cpt ;
  size_t           pos_sup ;
  int              sgn ;

  size_t           ind_pos_sup[3] ;
  size_t           ind_pos_cmp[3] ;

  ecs_tab_int_t    cpt_ref_inf ;

  ecs_vec_int_t  * vec_recherche ;

  ecs_vec_int_t  * vec_def_sup ;  /* Déf. entité supérieure (à compacter) */

  ecs_tab_int_t    tab_transf ;   /* Tableau de transformation */


  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  /* Initialisations */
  /*-----------------*/

  /* Initialisations avec précautions pour cas vide */

  tab_transf.nbr = 0 ;
  tab_transf.val = NULL ;

  if (this_vec_def == NULL)
    return tab_transf ;

  vec_def_sup = *this_vec_def ;

  if (vec_def_sup == NULL)
    return tab_transf ;

  nbr_sup_ini = vec_def_sup->pos_nbr - 1 ;
  cpt_sup_fin = 0;

  if (nbr_sup_ini < 1)
    return tab_transf ;

  if (signe_elt != NULL) {
    signe_elt->nbr = nbr_sup_ini ;
    BFT_MALLOC(signe_elt->val, nbr_sup_ini, ecs_int_t) ;
  }


  /* Comptage du nombre de sous-entités */

  nbr_inf = 0 ;

  for (ind_sup = 0 ;
       ind_sup < vec_def_sup->pos_tab[nbr_sup_ini] - 1 ;
       ind_sup++) {
    num_inf = vec_def_sup->val_tab[ind_sup];
    if ((size_t)(ECS_ABS(num_inf)) > nbr_inf)
      nbr_inf = ECS_ABS(num_inf);
  }

  /* Construction du tableau de recherche des sous-entités */
  /*-------------------------------------------------------*/

  vec_recherche = ecs_vec_int__alloue(nbr_inf + 1,
                                      vec_def_sup->pos_nbr - 1) ;

  /*
    Comptage pour chaque élément de l'entité inférieure du nombre de fois où il
    apparaît comme élement d'indice minimal d'un élément de l'entité supérieure
  */

  cpt_ref_inf = ecs_tab_int__cree_init(nbr_inf, 0) ;


  for (ind_sup = 0 ; ind_sup < vec_def_sup->pos_nbr - 1 ; ind_sup++) {

    ind_pos_loc = vec_def_sup->pos_tab[ind_sup] - 1 ;

    num_inf_min = ECS_ABS(vec_def_sup->val_tab[ind_pos_loc]) ;
    while (++ind_pos_loc < vec_def_sup->pos_tab[ind_sup + 1] -1) {
      num_inf_loc = ECS_ABS(vec_def_sup->val_tab[ind_pos_loc]) ;
      if (num_inf_loc < num_inf_min)
        num_inf_min = num_inf_loc ;
    }

    assert(num_inf_min > 0 && (size_t)num_inf_min <= nbr_inf) ;

    cpt_ref_inf.val[num_inf_min - 1] += 1 ;

  }


  /* Construction du vecteur */

  vec_recherche->pos_tab[0] = 1 ;

  for (ind_inf = 0 ; ind_inf < nbr_inf ; ind_inf++) {

    vec_recherche->pos_tab[ind_inf + 1]
      = vec_recherche->pos_tab[ind_inf] + cpt_ref_inf.val[ind_inf] ;

    cpt_ref_inf.val[ind_inf] = 0;

  }


  for (ind_sup = 0 ; ind_sup < vec_def_sup->pos_nbr - 1 ; ind_sup++) {

    ind_pos_loc = vec_def_sup->pos_tab[ind_sup] - 1 ;

    num_inf_min = ECS_ABS(vec_def_sup->val_tab[ind_pos_loc]) ;
    while (ind_pos_loc < vec_def_sup->pos_tab[ind_sup + 1] -1) {
      num_inf_loc = ECS_ABS(vec_def_sup->val_tab[ind_pos_loc]) ;
      ind_pos_loc++ ;
      if (num_inf_loc < num_inf_min)
        num_inf_min = num_inf_loc ;
    }

    ind_inf = num_inf_min - 1 ;

    ind_pos =   vec_recherche->pos_tab[ind_inf] - 1
              + cpt_ref_inf.val[ind_inf] ;

    cpt_ref_inf.val[ind_inf] += 1 ;

    vec_recherche->val_tab[ind_pos] = ind_sup + 1 ;

  }


  /* Libération du tableau auxiliaire */

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

  /* Allocation du tableau de transformation */

  tab_transf = ecs_tab_int__cree(nbr_sup_ini) ;

  for (ind_sup = 0 ; ind_sup < nbr_sup_ini ; ind_sup++)
    tab_transf.val[ind_sup] = ind_sup ;


  /* Boucle principale de recherche sur les éléments supérieurs */
  /*------------------------------------------------------------*/

  /*
    Le premier élément ne peut pas être fusionné avec un élément précédent,
    la boucle commence donc au deuxième
  */

  tab_transf.val[0] = 0 ;

  if (signe_elt != NULL)
    signe_elt->val[0] = 1 ;

  cpt_sup_fin       = 1 ;

  for (ind_sup = 1 ; ind_sup < vec_def_sup->pos_nbr - 1 ; ind_sup++) {

    /* Recherche élement entité inférieure de plus petit numéro référencé */

    ind_pos_sup[0] = vec_def_sup->pos_tab[ind_sup    ] - 1 ; /* début */
    ind_pos_sup[1] = vec_def_sup->pos_tab[ind_sup + 1] - 1 ; /* fin */
    ind_pos_sup[2] = vec_def_sup->pos_tab[ind_sup    ] - 1 ; /* plus petit */

    ind_pos_loc = ind_pos_sup[0] ;

    num_inf_min = ECS_ABS(vec_def_sup->val_tab[ind_pos_loc]) ;
    while (++ind_pos_loc < ind_pos_sup[1]) {
      num_inf_loc = ECS_ABS(vec_def_sup->val_tab[ind_pos_loc]) ;
      if (num_inf_loc < num_inf_min) {
        num_inf_min    = num_inf_loc ;
        ind_pos_sup[2] = ind_pos_loc ;
      }
    }

    /*
      On cherche des éléments de l'entité courante de plus petit numéro que
      l'entité courante ayant même plus petit élément de l'entité inférieure
      (recherche de candidats pour la fusion)
    */

    ind_inf = num_inf_min - 1 ;
    sgn     = 0 ;

    for (pos_cmp = vec_recherche->pos_tab[ind_inf]     - 1 ;
         pos_cmp < vec_recherche->pos_tab[ind_inf + 1] - 1 ;
         pos_cmp++) {

      ind_cmp = vec_recherche->val_tab[pos_cmp] - 1 ;

      /* Repérage point de départ pour comparaison */

      if (ind_cmp < ind_sup) {

        ind_pos_cmp[0] = vec_def_sup->pos_tab[ind_cmp    ] - 1 ; /* début */
        ind_pos_cmp[1] = vec_def_sup->pos_tab[ind_cmp + 1] - 1 ; /* fin */
        ind_pos_cmp[2] = vec_def_sup->pos_tab[ind_cmp    ] - 1 ; /* plus petit */

        assert(ind_pos_cmp[1] > ind_pos_cmp[0]);

        ind_pos_cmp[1] = vec_def_sup->pos_tab[ind_cmp + 1] - 1 ; /* fin */

        ind_pos_loc = ind_pos_cmp[0] ;

        num_inf_min_cmp = ECS_ABS(vec_def_sup->val_tab[ind_pos_loc]) ;
        while ((++ind_pos_loc) < ind_pos_cmp[1]) {
          num_inf_loc = ECS_ABS(vec_def_sup->val_tab[ind_pos_loc]) ;
          if (num_inf_loc < num_inf_min_cmp) {
            num_inf_min_cmp = num_inf_loc ;
            ind_pos_cmp[2]  = ind_pos_loc ;
          }
        }

        /* Comparaison des définitions */

        for (sgn = 1 ; sgn > -2 ; sgn -= 2) {

          ind_loc_sup = ind_pos_sup[2] ;
          ind_loc_cmp = ind_pos_cmp[2] ;

          do {

            ind_loc_sup++ ;
            if (ind_loc_sup == ind_pos_sup[1])
              ind_loc_sup = ind_pos_sup[0] ;

            ind_loc_cmp += sgn ;
            if (ind_loc_cmp == ind_pos_cmp[1])
              ind_loc_cmp = ind_pos_cmp[0] ;
            else if (ind_loc_cmp < ind_pos_cmp[0] || ind_loc_cmp > ind_pos_cmp[1])
              ind_loc_cmp = ind_pos_cmp[1] - 1 ;

          } while (   (   ECS_ABS(vec_def_sup->val_tab[ind_loc_sup])
                       == ECS_ABS(vec_def_sup->val_tab[ind_loc_cmp]))
                   && ind_loc_sup != ind_pos_sup[2]
                   && ind_loc_cmp != ind_pos_cmp[2]) ;

          if (   ind_loc_sup == ind_pos_sup[2]
              && ind_loc_cmp == ind_pos_cmp[2])
            break ; /* Sortie boucle sur signe parcours (1, -1, -3) */

        }

        /*
          Si sgn =  1, les entités sont confondues, de même sens ;
          Si sgn = -1, elles sont confondues, de sens inverse ;
          Sinon, elles ne sont pas confondues
        */

        if (sgn == 1 || sgn == -1)
          break ; /* Sortie boucle sur pos_cmp pour recherche de candidats */

      } /* Fin de la comparaison référence/candidat (if ind_cmp < ind_sup) */

    } /* Fin boucle sur pos_cmp recherche de candidats */


    /* Si on a trouvé une entité à fusionner */

    if (sgn == 1 || sgn == -1) {

      assert(pos_cmp < vec_recherche->pos_tab[ind_inf + 1] - 1) ;

      if (sgn == 1 && (ind_pos_cmp[1] - ind_pos_cmp[0] == 2)) {

        /*
          Si l'on a deux entités inférieures par entité supérieure,
          la permutation cyclique peut trouver une égalité quel
          que soit le signe ; on le corrige si nécessaire
        */

        if (   ECS_ABS(vec_def_sup->val_tab[ind_pos_sup[0]])
            != ECS_ABS(vec_def_sup->val_tab[ind_pos_cmp[0]]))
          sgn = -1 ;

      }

      tab_transf.val[ind_sup] = tab_transf.val[ind_cmp] ;

      if (signe_elt != NULL)
        signe_elt->val[ind_sup] = sgn ;

    }
    else {

      tab_transf.val[ind_sup] = cpt_sup_fin++ ;

      if (signe_elt != NULL)
        signe_elt->val[ind_sup] = 1 ;

    }


  } /* Fin boucle sur les éléments de l'entité supérieure (que l'on fusionne) */

  /* Libération du tableau de recherche */

  ecs_vec_int__detruit(vec_recherche) ;


  /* Compactage de la définition */
  /*-----------------------------*/

  /*
    Le premier élément ne peut pas être fusionné avec un élément précédent,
    la boucle commence donc au deuxième
  */

  cpt_sup_fin       = 1 ;

  for (ind_sup = 1 ; ind_sup < vec_def_sup->pos_nbr - 1 ; ind_sup++) {

    if (tab_transf.val[ind_sup] > (ecs_int_t)(cpt_sup_fin - 1)) {

      /* Recopie définition sur partie compactée du tableau */

      for (pos_cpt = vec_def_sup->pos_tab[cpt_sup_fin],
             pos_sup = vec_def_sup->pos_tab[ind_sup] ;
           pos_sup < vec_def_sup->pos_tab[ind_sup + 1] ;
           pos_cpt++, pos_sup++)
        vec_def_sup->val_tab[pos_cpt - 1] = vec_def_sup->val_tab[pos_sup - 1] ;

      cpt_sup_fin += 1 ;

      vec_def_sup->pos_tab[cpt_sup_fin] = pos_cpt ;

    }

  }

  /* Redimensionnement de l'entité compactée */

  ecs_vec_int__redimensionne(vec_def_sup,
                             cpt_sup_fin + 1,
                             vec_def_sup->pos_tab[cpt_sup_fin] - 1) ;

  /* Retour du tableau de transformation */

  return tab_transf ;

}


/*----------------------------------------------------------------------------
 *  Fonction qui fusionne des sommets équivalents
 *  Le pointeur sur un tableau de tolérances peut être à NULL, dans quel cas
 *  on prendra le même poids pour tous les sommets.
 *  Le pointeur sur le tableau liste_som_new peut aussi être à NULL, si l'on
 *  n'est pas intéressé par la liste des sommets fusionnés.
 *
 *  Remarque : le tableau d'équivalence (fusion) des sommets est construit de
 *             manière à ce qu'à un sommet ne comportant pas d'équivalent
 *             où de plus petit indice parmi ses équivalents corresponde la
 *             valeur -1, alors qu'un un sommet possédant des équivalents de
 *             plus petit indice corresponde le plus grand indice parmi ces
 *             équivalents (ce qui constitue une sorte de liste chaînée).
 *----------------------------------------------------------------------------*/

ecs_vec_int_t * ecs_vec_def__fusion_som
(
       ecs_vec_int_t     *const vec_def_are,    /* <-> Déf. arêtes            */
       ecs_vec_real_t    *const vec_def_som,    /* <-> Déf. sommets           */
       ecs_tab_int_t     *const tab_equiv_som,  /* <-> Fusions de sommets     */
 const ecs_tab_real_t    *const dist_max_som,   /*  -> Tolérances             */
       ecs_tab_int_t     *const liste_som_new   /* <-  Liste sommets modifiés */
)
{

  size_t     nbr_som_old ;
  size_t     nbr_som_new ;
  size_t     nbr_som_fus ;

  size_t     ind_som ;
  ecs_int_t  ind_som_loc ;
  ecs_int_t  ind_som_tmp ;

  ecs_int_t  icoo ;

  ecs_real_t  som_poids ;
  ecs_real_t  tmp_poids ;

  ecs_point_t  som_coord ;

  ecs_vec_int_t  * vec_som_old_new ;


  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  bft_printf(_("\n  Merging vertices:\n")) ;


  /* Étape préliminaire : transformation du tableau d'équivalence des */
  /*                      sommets en liste chaînée simple.            */

  ecs_loc_vec_def__transf_equiv(vec_def_som,
                                tab_equiv_som) ;

  /* Initialisations */
  /* --------------- */

  nbr_som_old = vec_def_som->pos_nbr - 1 ;

  bft_printf(_("    Initial number of vertices        : %10d\n"), nbr_som_old) ;

  vec_som_old_new = ecs_vec_int__alloue(nbr_som_old + 1,
                                        nbr_som_old) ;


  /* Initialisation du tableau de renumérotation des sommets */

  for (ind_som = 0 ; ind_som < nbr_som_old ; ind_som++) {

    vec_som_old_new->pos_tab[ind_som] = ind_som + 1 ;

    vec_som_old_new->val_tab[ind_som] = 0 ;

  }

  vec_som_old_new->pos_tab[nbr_som_old] = nbr_som_old + 1 ;


  /* Boucle principale sur les sommets */
  /*-----------------------------------*/

  nbr_som_new = 0 ;

  for (ind_som = 0 ; ind_som < nbr_som_old ; ind_som++) {

    /* Si le sommet n'a pas déjà été traité */

    if (vec_som_old_new->val_tab[ind_som] == 0) {

      /* Initialisation du sommet */

      som_poids   = 0.0 ;

      for (icoo = 0 ; icoo < ECS_DIM_3 ; icoo++)
        som_coord[icoo] = 0.0 ;

      ind_som_loc = ind_som ;

      /* Contribution (et marquage) des sommets équivalents */

      do {

        vec_som_old_new->val_tab[ind_som_loc] = nbr_som_new + 1 ;

        if (dist_max_som != NULL)
          tmp_poids = 1.0 / dist_max_som->val[ind_som_loc] ;

        else
          tmp_poids = 1.0 ;

        som_poids += tmp_poids ;

        for (icoo = 0 ; icoo < ECS_DIM_3 ; icoo++)
          som_coord[icoo]
            += (tmp_poids
                * vec_def_som->val_tab[ECS_DIM_3*ind_som_loc + icoo]) ;

        ind_som_loc = tab_equiv_som->val[ind_som_loc] ;

      } while (ind_som_loc != -1) ;

      /* Coordonnées finales */

      for (icoo = 0 ; icoo < ECS_DIM_3 ; icoo++)
        vec_def_som->val_tab[ECS_DIM_3*nbr_som_new + icoo]
          = som_coord[icoo] / som_poids ;

      /* On n'oublie pas d'incrémenter le nombre de sommets après fusion */

      nbr_som_new += 1 ;

    }

  }

  /* Fin de la boucle principale */
  /* --------------------------- */

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

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

  bft_printf(_("    Number of vertices after merging  : %10d\n"), nbr_som_new) ;


  /* Marquage des sommets issus d'une fusion ou intersection */
  /*---------------------------------------------------------*/

  /*
    On remarquera qu'une intersection entraîne toujours une fusion, les
    sommets étant crées sur chaque arête et marqués comme équivalents.
  */

  nbr_som_fus = 0 ;

  /*
    On n'aura plus besoin par la suite du tableau d'équivalence des sommets ;
    On le modifie donc de manière à ce que seul le premier sommet d'une liste
    chaînée de sommets équivalents (le sommet "récepteur") pointe sur son
    premier équivalent, pour s'en servir comme marqueur.
  */

  for (ind_som = 0 ; ind_som < tab_equiv_som->nbr ; ind_som++) {

    if (tab_equiv_som->val[ind_som] != -1) {

      nbr_som_fus += 1 ;

      ind_som_loc = tab_equiv_som->val[ind_som] ;

      while (ind_som_loc != -1) {

        ind_som_tmp = ind_som_loc ;
        ind_som_loc = tab_equiv_som->val[ind_som_loc] ;

        tab_equiv_som->val[ind_som_tmp] = -1 ;

      }

    }

  }

  /* Allocation de la liste des sommets issus d'une fusion */

  if (liste_som_new != NULL) {

    liste_som_new->nbr = nbr_som_fus ;

    BFT_MALLOC(liste_som_new->val, liste_som_new->nbr, ecs_int_t) ;


    /* Remplissage de la liste */

    nbr_som_fus = 0 ;

    for (ind_som = 0 ; ind_som < tab_equiv_som->nbr ; ind_som++) {

      if (tab_equiv_som->val[ind_som] != -1) {

        liste_som_new->val[nbr_som_fus]
          = vec_som_old_new->val_tab[ind_som] - 1 ;

        nbr_som_fus += 1 ;

      }

    }

    assert (nbr_som_fus == liste_som_new->nbr) ;

  }

  bft_printf(_("    Number of modified vertices       : %10d\n"), nbr_som_fus) ;


  /* Mise à jour de la définition des arêtes */
  /*-----------------------------------------*/

  for (ind_som = 0 ;
       ind_som < vec_def_are->pos_tab[vec_def_are->pos_nbr - 1] - 1 ;
       ind_som++)

    vec_def_are->val_tab[ind_som]
      = vec_som_old_new->val_tab[vec_def_are->val_tab[ind_som] - 1] ;


  /* Retour */

  return vec_som_old_new ;

}


/*----------------------------------------------------------------------------
 *  Suppression des sommets ne participant pas à la connectivité
 *   (renvoie le vecteur de renumérotation des sommets, ou NULL si aucune
 *   renumérotation n'a été nécessaire).
 *----------------------------------------------------------------------------*/

ecs_vec_int_t * ecs_vec_def__nettoie_nodal
(
       ecs_vec_real_t  *const vec_def_som,
 const ecs_vec_int_t   *const vec_def_are,
 const ecs_vec_int_t   *const vec_def_fac,
 const ecs_vec_int_t   *const vec_def_cel
)
{

  size_t          cpt_som ;
  size_t          iloc ;
  size_t          ipos_new ;
  size_t          ipos_old ;
  size_t          isom ;
  size_t          ival ;
  size_t          nbr_som ;
  size_t          nbr_val ;
  size_t          pas_dim ;

  ecs_tab_int_t   tab_ref_som ;

  ecs_vec_int_t  * vec_som_old_new ;


  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  assert(vec_def_som != NULL) ;


  /* marquage des sommets utilisés */

  nbr_som = vec_def_som->pos_nbr - 1 ;

  tab_ref_som = ecs_tab_int__cree_init(nbr_som, 0) ;


  /* arêtes */

  if (vec_def_are != NULL) {

    nbr_val = ecs_vec_int__ret_val_nbr(vec_def_are) ;

    for (ival = 0 ; ival < nbr_val ; ival++)
      tab_ref_som.val[vec_def_are->val_tab[ival] - 1] = 1 ;

  }


  /* faces */

  if (vec_def_fac != NULL) {

    nbr_val = ecs_vec_int__ret_val_nbr(vec_def_fac) ;

    for (ival = 0 ; ival < nbr_val ; ival++)
      tab_ref_som.val[vec_def_fac->val_tab[ival] - 1] = 1 ;

  }


  /* cellules */

  if (vec_def_cel != NULL) {

    nbr_val = ecs_vec_int__ret_val_nbr(vec_def_cel) ;

    for (ival = 0 ; ival < nbr_val ; ival++)
      tab_ref_som.val[vec_def_cel->val_tab[ival] - 1] = 1 ;

  }


  /* Préparation de la renumérotation */

  cpt_som = 0 ;

  for (isom = 0 ; isom < nbr_som ; isom++) {

    if (tab_ref_som.val[isom] != 0)
      cpt_som++ ;

  }


  /* Si tous les sommets sont référencés, rien à faire */

  if (cpt_som == nbr_som) {

    BFT_FREE(tab_ref_som.val) ;

    return NULL ;

  }


  /* Initialisation du vecteur de renumérotation */

  vec_som_old_new = ecs_vec_int__alloue(nbr_som + 1,
                                        cpt_som) ;

  vec_som_old_new->pos_tab[0] = 1 ;


  /* Compactage du tableau des sommets */

  cpt_som = 0 ;

  pas_dim = vec_def_som->pos_pas ;


  for (isom = 0 ; isom < nbr_som ; isom++) {

    if (tab_ref_som.val[isom] != 0) {

      ipos_new = pas_dim * (cpt_som) ;
      ipos_old = pas_dim * isom ;

      for (iloc = 0 ; iloc < pas_dim ; iloc++)
        vec_def_som->val_tab[ipos_new++] = vec_def_som->val_tab[ipos_old++] ;

      vec_som_old_new->val_tab[cpt_som] = cpt_som + 1 ;

      cpt_som++ ;

    }

    vec_som_old_new->pos_tab[isom + 1] = cpt_som + 1 ;

  }


  BFT_FREE(tab_ref_som.val) ;

  ecs_vec_real__redimensionne(vec_def_som,
                              cpt_som + 1,
                              cpt_som * pas_dim) ;


  return vec_som_old_new ;

}


/*----------------------------------------------------------------------------
 *  Correction si nécessaire de l'orientation des éléments en
 *   connectivité nodale.
 *----------------------------------------------------------------------------*/

void ecs_vec_def__orient_nodal
(
       ecs_vec_real_t  *const vec_def_som,    /*  -> Déf. sommets             */
 const ecs_vec_int_t   *const vec_def_fac,    /*  -> Déf. faces               */
 const ecs_vec_int_t   *const vec_def_cel,    /*  -> Déf. cellules            */
       ecs_tab_int_t   *const liste_cel_err,  /* <-  Cellules avec erreur     */
       ecs_tab_int_t   *const liste_cel_cor,  /* <-  Cellules corrigées       */
 const ecs_bool_t             correc_orient   /*  -> Correction ou non        */
)
{

  size_t      ipos_cel ;
  size_t      ipos_fac ;
  size_t      icel ;
  size_t      ifac ;
  size_t      nbr_som_loc ;
  size_t      nbr_fac ;
  size_t      nbr_cel ;
  ecs_int_t   typ_elt ;
  ecs_int_t   ret_orient ;

  ecs_int_t   cpt_orient_correc[ECS_ELT_TYP_FIN] ;
  ecs_int_t   cpt_orient_erreur[ECS_ELT_TYP_FIN] ;

  ecs_int_t   typ_cel_base[9] = {ECS_ELT_TYP_NUL,
                                 ECS_ELT_TYP_NUL,
                                 ECS_ELT_TYP_NUL,
                                 ECS_ELT_TYP_NUL,
                                 ECS_ELT_TYP_CEL_TETRA,
                                 ECS_ELT_TYP_CEL_PYRAM,
                                 ECS_ELT_TYP_CEL_PRISM,
                                 ECS_ELT_TYP_NUL,
                                 ECS_ELT_TYP_CEL_HEXA};

  ecs_int_t   cpt_cel_erreur = 0 ;
  ecs_int_t   cpt_cel_correc = 0 ;

  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  assert(vec_def_som != NULL) ;
  assert(vec_def_som->pos_pas == 3) ;

  for (typ_elt = 0 ; typ_elt < ECS_ELT_TYP_FIN ; typ_elt++) {
    cpt_orient_correc[typ_elt] = 0 ;
    cpt_orient_erreur[typ_elt] = 0 ;
  }


  bft_printf(_("\n  Element orientation check.\n\n")) ;

  /* faces */

  if (vec_def_fac != NULL) {

    nbr_fac = ecs_vec_int__ret_pos_nbr(vec_def_fac) - 1 ;

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

      ipos_fac    = vec_def_fac->pos_tab[ifac    ] - 1;
      nbr_som_loc = vec_def_fac->pos_tab[ifac + 1] - 1 - ipos_fac;

      if (nbr_som_loc == 4) {

        ret_orient
          = ecs_loc_vec_def__orient_quad(vec_def_som->val_tab,
                                         &(vec_def_fac->val_tab[ipos_fac]),
                                         correc_orient) ;

        if (ret_orient < 0)
          cpt_orient_erreur[ECS_ELT_TYP_FAC_QUAD] += 1 ;
        else if (ret_orient > 0)
          cpt_orient_correc[ECS_ELT_TYP_FAC_QUAD] += 1 ;

      }

    }

  }

  /* cellules */

  if (vec_def_cel != NULL) {

    ecs_tab_int_t face_index, face_marker, edges ;

    face_index.nbr = 0;
    face_index.val = NULL;

    face_marker.nbr = 0;
    face_marker.val = NULL;

    edges.nbr = 0;
    edges.val = NULL;

    nbr_cel = ecs_vec_int__ret_pos_nbr(vec_def_cel) - 1 ;

    for (icel = 0 ; icel < nbr_cel ; icel++) {

      ipos_cel = vec_def_cel->pos_tab[icel] - 1;

      nbr_som_loc = vec_def_cel->pos_tab[icel + 1] - 1 - ipos_cel ;

      if (nbr_som_loc < 9)
        typ_elt = typ_cel_base[nbr_som_loc] ;
      else
        typ_elt = ECS_ELT_TYP_CEL_POLY ;

      switch (typ_elt) {

      case ECS_ELT_TYP_CEL_TETRA:

        ret_orient
          = ecs_loc_vec_def__orient_tetra(vec_def_som->val_tab,
                                          &(vec_def_cel->val_tab[ipos_cel]),
                                          correc_orient) ;
        break ;

      case ECS_ELT_TYP_CEL_PYRAM:

        ret_orient
          = ecs_loc_vec_def__orient_pyram(vec_def_som->val_tab,
                                          &(vec_def_cel->val_tab[ipos_cel]),
                                          correc_orient) ;
        break ;

      case ECS_ELT_TYP_CEL_PRISM:

        ret_orient
          = ecs_loc_vec_def__orient_prism(vec_def_som->val_tab,
                                          &(vec_def_cel->val_tab[ipos_cel]),
                                          correc_orient) ;

        break ;

      case ECS_ELT_TYP_CEL_HEXA:

        ret_orient
          = ecs_loc_vec_def__orient_hexa(vec_def_som->val_tab,
                                         &(vec_def_cel->val_tab[ipos_cel]),
                                         correc_orient) ;

        break ;

      default: /* ECS_ELT_TYP_CEL_POLY */

        ret_orient
          = ecs_loc_vec_def__orient_polyedre(vec_def_som->val_tab,
                                             &(vec_def_cel->val_tab[ipos_cel]),
                                             (  vec_def_cel->pos_tab[icel + 1]
                                              - vec_def_cel->pos_tab[icel]),
                                             correc_orient,
                                             &face_index,
                                             &face_marker,
                                             &edges) ;

      } ;

      if (ret_orient < 0) {
        cpt_orient_erreur[typ_elt] += 1 ;
        if (liste_cel_err != NULL) {
          if (liste_cel_err->nbr == 0) {
            liste_cel_err->nbr = nbr_cel ;
            BFT_MALLOC(liste_cel_err->val, liste_cel_err->nbr, ecs_int_t) ;
          }
          liste_cel_err->val[cpt_cel_erreur] = icel ;
        }
        cpt_cel_erreur += 1 ;
      }

      else if (ret_orient > 0) {
        cpt_orient_correc[typ_elt] += 1 ;
        if (liste_cel_cor != NULL) {
          if (liste_cel_cor->nbr == 0) {
            liste_cel_cor->nbr = nbr_cel ;
            BFT_MALLOC(liste_cel_cor->val, liste_cel_cor->nbr, ecs_int_t) ;
          }
          liste_cel_cor->val[cpt_cel_correc] = icel ;
        }
        cpt_cel_correc += 1 ;
      }

    }

    if (face_index.nbr > 0)
      BFT_FREE(face_index.val);

    if (face_marker.nbr > 0)
      BFT_FREE(face_marker.val);

    if (edges.nbr > 0)
      BFT_FREE(edges.val);
  }

  /* Impression de messages d'avertissement */

  for (typ_elt = 0 ; typ_elt < ECS_ELT_TYP_FIN ; typ_elt++) {
    if (cpt_orient_correc[typ_elt] > 0) {
      ecs_warn() ;
      bft_printf(_("%d elements of type %s had to be re-oriented\n"),
                 (int)(cpt_orient_correc[typ_elt]),
                 _(ecs_fic_elt_typ_liste_c[typ_elt].nom)) ;
    }
    if (cpt_orient_erreur[typ_elt] > 0) {
      if (correc_orient == ECS_TRUE) {
        ecs_warn() ;
        bft_printf(_("%d elements of type %s were impossible to re-orient\n"),
                   (int)(cpt_orient_erreur[typ_elt]),
                   _(ecs_fic_elt_typ_liste_c[typ_elt].nom)) ;
      }
      else {
        ecs_warn() ;
        bft_printf(_("%d elements of type %s are mis-oriented\n"),
                   (int)(cpt_orient_erreur[typ_elt]),
                   _(ecs_fic_elt_typ_liste_c[typ_elt].nom)) ;
      }
    }
  }

  /* Redimensionnement de la liste des cellules toujours mal orientées */

  if (liste_cel_err != NULL) {
    if (liste_cel_err->nbr > 0) {
      liste_cel_err->nbr = cpt_cel_erreur ;
      BFT_REALLOC(liste_cel_err->val, liste_cel_err->nbr, ecs_int_t) ;
    }
  }
  if (liste_cel_cor != NULL) {
    if (liste_cel_cor->nbr > 0) {
      liste_cel_cor->nbr = cpt_cel_correc ;
      BFT_REALLOC(liste_cel_cor->val, liste_cel_cor->nbr, ecs_int_t) ;
    }
  }

}


/*----------------------------------------------------------------------------
 *  Fonction qui calcule les coordonnées minimales et maximales
 *----------------------------------------------------------------------------*/

void ecs_vec_def__calc_coo_ext
(
 ecs_vec_real_t * vec_som
)
{

  size_t      icoo ;
  size_t      ipos ;
  size_t      isom ;

  ecs_real_t  coo_min[3] ;
  ecs_real_t  coo_max[3] ;

  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  if (vec_som->pos_nbr < 2)
    return ;

  ipos = 0 ;

  for (icoo = 0 ; icoo < ECS_DIM_3 ; icoo++) {
    coo_min[icoo] = vec_som->val_tab[ipos + icoo] ;
    coo_max[icoo] = vec_som->val_tab[ipos + icoo] ;
  }

  for (isom = 1 ; isom < vec_som->pos_nbr - 1 ; isom++) {

    ipos = vec_som->pos_pas * isom ;

    for (icoo = 0 ; icoo < ECS_DIM_3 ; icoo++) {
      if (vec_som->val_tab[ipos + icoo] < coo_min[icoo])
        coo_min[icoo] = vec_som->val_tab[ipos + icoo] ;
      else if (vec_som->val_tab[ipos + icoo] > coo_max[icoo])
        coo_max[icoo] = vec_som->val_tab[ipos + icoo] ;
    }

  }

  bft_printf(_("\n  Domain coordinate extents:\n\n")) ;

  bft_printf("  [% 10.5e, % 10.5e, % 10.5e]\n",
             coo_min[0], coo_min[1], coo_min[2]) ;
  bft_printf("  [% 10.5e, % 10.5e, % 10.5e]\n",
             coo_max[0], coo_max[1], coo_max[2]) ;

}


/*----------------------------------------------------------------------------
 *  Fonction qui modifie les coordonnées du maillage
 *----------------------------------------------------------------------------*/

void ecs_vec_def__transf_coo
(
 ecs_vec_real_t  * vec_som,
 const double      matrice[3][4]
)
{

  size_t      icoo ;
  size_t      ipos ;
  size_t      isom ;

  ecs_real_t  coo_tmp[3] ;

  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  if (vec_som->pos_nbr < 2)
    return ;

  for (isom = 0; isom < vec_som->pos_nbr - 1 ; isom++) {

    ipos = vec_som->pos_pas * isom ;

    for (icoo = 0 ; icoo < ECS_DIM_3 ; icoo++)
      coo_tmp[icoo] = vec_som->val_tab[ipos + icoo] ;

    for (icoo = 0 ; icoo < ECS_DIM_3 ; icoo++)
      vec_som->val_tab[ipos + icoo]
        = (  matrice[icoo][0] * coo_tmp[0]
           + matrice[icoo][1] * coo_tmp[1]
           + matrice[icoo][2] * coo_tmp[2]
           + matrice[icoo][3] * 1.0) ;

  }

}


/*----------------------------------------------------------------------------
 *  Fonction qui réalise la fusion des définitions des éléments
 *----------------------------------------------------------------------------*/

ecs_tab_int_t ecs_vec_real_def__fusionne
(
 ecs_vec_real_t * *const this_vec_def_real
)
{

  ecs_vec_real_t * vec_elt_cpct ;  /* Vecteur compacte des éléments */

  ecs_tab_int_t    vect_transf ;   /* Vecteur de transformation */


  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  /* Ordonnancement des éléments en fonction de leur définition */
  /*------------------------------------------------------------*/


  vect_transf = ecs_loc_vec_real_def__trie_elt(this_vec_def_real,
                                               ECS_REAL_PRECISION) ;


  /* Compactage de la liste des éléments                      */
  /*  par suppression des éléments géométriquement identiques */
  /*----------------------------------------------------------*/


  vec_elt_cpct = ecs_loc_vec_real_def__compacte(*this_vec_def_real,
                                                &vect_transf ) ;

  ecs_vec_real__detruit(*this_vec_def_real) ;


  /* Affectation de la nouvelle liste compactée */
  /*--------------------------------------------*/


  *this_vec_def_real = vec_elt_cpct ;


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


  return vect_transf ;

}


/*----------------------------------------------------------------------------
 *  Fonction qui renvoie un tableau associant un type à chaque face, sous
 * forme de masque : 0 pour face isolée, 1 ou 2 pour face de bord (1 si
 * cellule avec cette face normale sortante, 2 si cellule avec cette face
 * normale entrante), 1+2 = 3 pour face interne, et 4 ou plus pour tous
 * les autres cas, correspondant à une erreur de connectivité (+4 pour faces
 * voyant au moins deux cellules avec face normale sortante, +8 pour faces
 * voyant au moins deux cellules avec face normale entrante).
 *
 *  Le type de chaque face pourra être modifié ultérieurement en fonction
 * des informations de périodicité.
 *----------------------------------------------------------------------------*/

ecs_tab_int_t ecs_vec_def__typ_fac_cel
(
 const ecs_vec_int_t  *const vec_def_cel,
       size_t                nbr_fac
)
{

  size_t      nbr_cel ;
  ecs_int_t   num_fac ;
  size_t      icel ;
  size_t      ifac ;
  size_t      ipos ;

  ecs_tab_int_t   typ_fac ;


  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  assert(vec_def_cel != NULL) ;


  /* Initialisations */

  nbr_cel   = ecs_vec_int__ret_pos_nbr(vec_def_cel) - 1 ;

  /* Type de face selon le nombre de cellules voisines : 0 pour face isolée,
     1 ou 2 pour face de bord, 3 pour faces internes, et 4 pour tous les
     autres cas (faces voyant au moins deux cellules sur un même côté) */

  typ_fac.nbr = nbr_fac;

  BFT_MALLOC(typ_fac.val, nbr_fac, ecs_int_t) ;

  for (ifac = 0 ; ifac < nbr_fac ; ifac++)
    typ_fac.val[ifac] = 0 ;


  /* Boucle sur les cellules : marquage */
  /*------------------------------------*/

  for (icel = 0 ; icel < nbr_cel ; icel++ ) {

    for (ipos = vec_def_cel->pos_tab[icel] - 1 ;
         ipos < vec_def_cel->pos_tab[icel+1] - 1 ;
         ipos++) {

      num_fac = vec_def_cel->val_tab[ipos] ;

      ifac = ECS_ABS(num_fac) - 1;

      if (num_fac > 0 && (typ_fac.val[ifac] & 1) == 0)
        typ_fac.val[ifac] += 1 ;

      else if (num_fac < 0 && (typ_fac.val[ifac] & 2) == 0)
        typ_fac.val[ifac] += 2 ;

      else {
        if (num_fac > 0 && (typ_fac.val[ifac] & 1) == 1)
          typ_fac.val[ifac] = typ_fac.val[ifac] | 4 ;
        else if (num_fac < 0 && (typ_fac.val[ifac] & 2) == 1)
          typ_fac.val[ifac] = typ_fac.val[ifac] | 8 ;
      }

    }

  }


  return typ_fac ;

}


/*----------------------------------------------------------------------------
 *  Fonction qui renvoie un tableau associant à chaque face les
 * numéros des cellules définies par cette face (normale sortante,
 * puis normale entrante). On affecte une valeur 0 lorsqu'il n'y a pas de
 * cellule correspondante directe (la périodicité n'est donc pas prise en
 * compte à ce niveau).
 *
 * On suppose que la cohérence du maillage a déjà été vérifiée et
 * qu'aucune face n'appartient à plus d'une cellule par côté.
 *----------------------------------------------------------------------------*/

ecs_tab_int_t ecs_vec_def__fac_cel
(
 const ecs_vec_int_t  *const vec_def_cel,
       size_t                nbr_fac
)
{

  size_t      nbr_cel ;
  ecs_int_t   num_fac ;
  size_t      icel ;
  size_t      ifac ;
  size_t      ipos ;

  ecs_tab_int_t   fac_cel ;


  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  assert(vec_def_cel != NULL) ;


  /* Initialisations */

  nbr_cel = ecs_vec_int__ret_pos_nbr(vec_def_cel) - 1 ;

  /* Allocation et mise à zéro des connectivités */

  fac_cel.nbr = nbr_fac * 2;

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

  for (ipos = 0 ; ipos < fac_cel.nbr ; ipos++)
    fac_cel.val[ipos] = 0 ;


  /* Boucle sur les cellules : marquage */
  /*------------------------------------*/

  for (icel = 0 ; icel < nbr_cel ; icel++ ) {

    for (ipos = vec_def_cel->pos_tab[icel] - 1 ;
         ipos < vec_def_cel->pos_tab[icel+1] - 1 ;
         ipos++) {

      num_fac = vec_def_cel->val_tab[ipos] ;

      ifac = ECS_ABS(num_fac) - 1;

      if (num_fac > 0) {
        assert(fac_cel.val[ifac*2] == 0) ;
        fac_cel.val[ifac*2] = icel + 1 ;
      }
      else {
        assert(fac_cel.val[ifac*2 + 1] == 0) ;
        fac_cel.val[ifac*2 + 1] = icel + 1 ;
      }

    }

  }


  return fac_cel ;

}


/*----------------------------------------------------------------------------
 *  Fonction qui renvoie un tableau associant à chaque cellule un code
 * en fonction des erreurs de connectivité éventuelles associées à cette
 * cellule (0 si pas d'erreur, 1 si une des faces définissant cette cellule
 * s'appuie sur plusieurs cellules du même côté).
 *----------------------------------------------------------------------------*/

ecs_tab_int_t ecs_vec_def__err_cel_connect
(
 const ecs_vec_int_t  *const vec_def_cel,
 const ecs_tab_int_t  *const typ_fac_cel
)
{

  size_t      nbr_cel ;
  ecs_int_t   num_fac ;
  size_t      icel ;
  size_t      ifac ;
  size_t      ipos ;

  ecs_tab_int_t   typ_cell_connect ;


  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  assert(vec_def_cel != NULL) ;


  /* Initialisations */

  nbr_cel = ecs_vec_int__ret_pos_nbr(vec_def_cel) - 1 ;

  /* Allocation et mise à zéro du tableau */

  typ_cell_connect.nbr = nbr_cel;

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

  for (icel = 0 ; icel < typ_cell_connect.nbr ; icel++)
    typ_cell_connect.val[icel] = 0 ;


  /* Boucle sur les cellules */
  /*-------------------------*/

  for (icel = 0 ; icel < nbr_cel ; icel++ ) {

    for (ipos = vec_def_cel->pos_tab[icel] - 1 ;
         ipos < vec_def_cel->pos_tab[icel+1] - 1 ;
         ipos++) {

      num_fac = vec_def_cel->val_tab[ipos] ;

      ifac = ECS_ABS(num_fac) - 1;

      if (num_fac > 0) {
        if (typ_fac_cel->val[ifac] & 4)
          typ_cell_connect.val[icel] = 1 ;
      }
      else {
        if (typ_fac_cel->val[ifac] & 8)
          typ_cell_connect.val[icel] = 1 ;
      }

    }

  }


  return typ_cell_connect ;

}


/*----------------------------------------------------------------------------
 *  Fonction qui construit la liste des cellules attachées à une liste
 *  de faces fournie en argument.
 *----------------------------------------------------------------------------*/

ecs_tab_int_t ecs_vec_def__liste_cel_fac
(
       size_t                nbr_fac,
 const ecs_vec_int_t  *const vec_def_cel,
 const ecs_tab_int_t         liste_fac
)
{

  size_t      cpt_cel ;
  size_t      icel ;
  size_t      ifac ;
  size_t      iloc ;
  size_t      ipos ;
  size_t      nbr_cel ;
  size_t      nbr_fac_cel ;
  ecs_int_t  *indic_cel ;
  ecs_int_t  *indic_fac ;

  ecs_tab_int_t  liste_cel ;

  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  assert(vec_def_cel != NULL) ;


  /* Initialisations */

  cpt_cel = 0 ;

  nbr_cel = vec_def_cel->pos_nbr - 1 ;

  /* Indicateurs */

  BFT_MALLOC(indic_cel, nbr_cel, ecs_int_t) ;
  BFT_MALLOC(indic_fac, nbr_fac, ecs_int_t) ;

  for (icel = 0 ; icel < nbr_cel ; icel++)
    indic_cel[icel] = 0 ;

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

  for (ifac = 0 ; ifac < liste_fac.nbr ; ifac++)
    indic_fac[liste_fac.val[ifac]] = 1 ;


  /* Première boucle sur les cellules : marquage */

  for (icel = 0 ; icel < nbr_cel ; icel++) {

    nbr_fac_cel
      = vec_def_cel->pos_tab[icel + 1]
      - vec_def_cel->pos_tab[icel] ;

    ipos = vec_def_cel->pos_tab[icel] - 1 ;

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

      ifac = ECS_ABS(vec_def_cel->val_tab[ipos]) - 1 ;

      ipos++ ;

      if (indic_fac[ifac] == 1 && indic_cel[icel] == 0) {
        indic_cel[icel] = 1 ;
        cpt_cel++ ;
      }

    }

  }

  BFT_FREE(indic_fac) ;

  /* Seconde boucle sur les cellules : remplissage */

  liste_cel.nbr = cpt_cel ;
  BFT_MALLOC(liste_cel.val, liste_cel.nbr, ecs_int_t) ;

  cpt_cel = 0 ;

  for (icel = 0 ; icel < nbr_cel ; icel++) {

    if (indic_cel[icel] == 1)
      liste_cel.val[cpt_cel++] = icel ;

  }

  BFT_FREE(indic_cel) ;

  return liste_cel ;

}


/*----------------------------------------------------------------------------
 *  Fonction qui construit la table de connectivité "faces -> sommets"
 *----------------------------------------------------------------------------*/

ecs_vec_int_t * ecs_vec_def__cree_fac_som
(
 const ecs_vec_int_t *const vec_def_fac,
 const ecs_vec_int_t *const vec_def_are
)
{

  ecs_vec_int_t * vec_connect_fac_som ;

  size_t     nbr_fac ;
  size_t     nbr_val ;

  size_t     isom ;
  size_t     iare ;
  size_t     ifac ;

  size_t     iare_inf ;
  size_t     iare_sup ;

  ecs_int_t  num_som ;
  ecs_int_t  num_are ;

  ecs_int_t  cpt_val ;


  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  nbr_fac = vec_def_fac->pos_nbr - 1 ;
  nbr_val = ecs_vec_int__ret_val_nbr(vec_def_fac) ;

  vec_connect_fac_som = ecs_vec_int__alloue(vec_def_fac->pos_nbr,
                                            nbr_val) ;

  vec_connect_fac_som->pos_tab[0] = 1 ;


  cpt_val = 0 ;


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

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


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


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

      num_are  = vec_def_fac->val_tab[iare] ;

      if (num_are > 0)
        isom = vec_def_are->pos_tab[num_are          - 1]  - 1 ;
      else
        isom = vec_def_are->pos_tab[ECS_ABS(num_are) - 1] ;

      num_som  = vec_def_are->val_tab[isom] ;


      vec_connect_fac_som->val_tab[cpt_val++] = num_som ;


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


  } /* Fin : boucle sur les faces */


  return vec_connect_fac_som ;


}


/*----------------------------------------------------------------------------
 *  Fonction préparant le découpage des faces  polygonales en triangles
 *
 *  Les faces doivent être définies en connectivité nodale.
 *
 *  Si un polygone à n sommets est decoupé en triangles,
 *  on obtient n-2 triangles.
 *----------------------------------------------------------------------------*/

void ecs_vec_def__cpt_poly_tria
(
 const ecs_vec_int_t  *const vec_def_fac,      /*  -> définition des faces    */
       size_t         *const fac_nbr_som_max,  /* <-  nbr max sommets/face    */
       size_t         *const nbr_fac_new,      /* <-  nbr faces decoupées     */
       size_t         *const nbr_fac_val_new   /* <-  dim val faces decoupées */
)
{

  size_t        nbr_fac ;
  size_t        ifac ;
  size_t        fac_nbr_som ;


  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  nbr_fac = vec_def_fac->pos_nbr - 1 ;


  *fac_nbr_som_max = 0 ;

  *nbr_fac_new     = 0 ;
  *nbr_fac_val_new = 0 ;


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


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


    fac_nbr_som =   vec_def_fac->pos_tab[ifac + 1]
                  - vec_def_fac->pos_tab[ifac] ;


    if (fac_nbr_som > *fac_nbr_som_max)

      *fac_nbr_som_max = fac_nbr_som ;


    if (fac_nbr_som <= 4) {

      /* La face ne doit pas être découpée */

      *nbr_fac_new     += 1 ;
      *nbr_fac_val_new += fac_nbr_som ;

    }
    else {

      *nbr_fac_new     += (fac_nbr_som - 2) ;
      *nbr_fac_val_new += (fac_nbr_som - 2) * 3 ;

    }


  } /* Fin : boucle sur les faces */


}


/*----------------------------------------------------------------------------
 *  Fonction réalisant le découpage des faces  polygonales en triangles
 *
 *  Les faces doivent être définies en connectivité nodale.
 *----------------------------------------------------------------------------*/

ecs_vec_int_t * ecs_vec_def__dec_poly_tria
(
       ecs_vec_int_t      *const vec_def_fac_new,
 const ecs_vec_int_t      *const vec_def_fac,
 const ecs_vec_real_t     *const vec_def_som,
       size_t                    fac_nbr_som_max,
       size_t                    nbr_fac_new
)
{

  size_t     cpt_fac_new ;
  size_t     cpt_val_new ;
  size_t     nbr_som ;
  size_t     nbr_fac ;

  ecs_int_t  icoo ;
  ecs_int_t  isom ;
  size_t     ifac ;
  ecs_int_t  itri ;

  ecs_int_t  icoo_inf ;
  ecs_int_t  isom_inf ;
  ecs_int_t  isom_sup ;

  ecs_int_t  num_som ;

  ecs_int_t  pos_fac_new ;
  ecs_int_t  pos_fac_dec ;

  ecs_int_t  fac_nbr_som ;
  ecs_int_t  fac_nbr_tria ;

  ecs_int_t     * som_fac ;

  ecs_int_t     * liste_som_tria ;
  ecs_int_t     * liste_prec ;
  ecs_int_t     * liste_suiv ;
  ecs_int_t     * liste_are_def ;
  ecs_int_t     * liste_voisinage_are ;
  ecs_bool_t    * liste_are_loc_delaunay ;
  ecs_bool_t    * concave ;

  ecs_point_t   * coord_som ;

  ecs_vec_int_t * vec_def_fac_dec ;

  const ecs_int_t  nbr_are_tot_max = fac_nbr_som_max*(fac_nbr_som_max - 1)/2 ;


#define ECS_FCT_CALCULE_DISTANCE(vect1, vect2)           \
  sqrt((vect2[X] - vect1[X]) * (vect2[X] - vect1[X]) +   \
       (vect2[Y] - vect1[Y]) * (vect2[Y] - vect1[Y]) +   \
       (vect2[Z] - vect1[Z]) * (vect2[Z] - vect1[Z])   )


  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  assert(vec_def_som->pos_pas == ECS_DIM_3) ;

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


  /* Définition des anciennes faces en fonction des nouvelles */

  vec_def_fac_dec  = ecs_vec_int__alloue(nbr_fac + 1,
                                         nbr_fac_new) ;

  vec_def_fac_dec->pos_tab[0] = 1 ;


  /* Définition des faces dans la nouvelle liste */

  cpt_fac_new = 0 ;
  cpt_val_new = 0 ;

  vec_def_fac_new->pos_tab[0] = 1 ;


  /* Définition des nouvelles faces issues du découpage */

  BFT_MALLOC(som_fac, fac_nbr_som_max, ecs_int_t) ;

  BFT_MALLOC(coord_som, fac_nbr_som_max, ecs_point_t) ;
  BFT_MALLOC(liste_prec, fac_nbr_som_max, ecs_int_t) ;
  BFT_MALLOC(liste_suiv, fac_nbr_som_max, ecs_int_t) ;
  BFT_MALLOC(concave, fac_nbr_som_max, ecs_bool_t) ;

  /* Le nombre de triangles après découpage doit être ou égal a
     (nombre de sommets - 2) */

  BFT_MALLOC(liste_som_tria, (fac_nbr_som_max-2)*3, ecs_int_t) ;

  /* On alloue aussi des tableaux de travail (de taille le nombre d'arêtes
     maximal, soit C(n,2)=n(n-1)/2 ) pour la triangulation de Delaunay */

  BFT_MALLOC(liste_are_def, nbr_are_tot_max*2, ecs_int_t) ;
  BFT_MALLOC(liste_voisinage_are, nbr_are_tot_max*2, ecs_int_t) ;
  BFT_MALLOC(liste_are_loc_delaunay, nbr_are_tot_max, ecs_bool_t) ;


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


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


    isom_inf = vec_def_fac->pos_tab[ifac]     - 1 ;
    isom_sup = vec_def_fac->pos_tab[ifac + 1] - 1 ;

    fac_nbr_som = isom_sup - isom_inf ;


    /* Préparation : boucle sur les sommets de la face */

    if (fac_nbr_som > 4) {

      for (isom = 0 ; isom < fac_nbr_som ; isom++) {

        num_som  = vec_def_fac->val_tab[isom_inf + isom] ;

        som_fac[isom] = num_som ;

        icoo_inf = vec_def_som->pos_pas * (num_som - 1) ;

        for (icoo = 0 ; icoo < ECS_DIM_3 ; icoo++)
          coord_som[isom][icoo] = vec_def_som->val_tab[icoo_inf + icoo] ;

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

    }


    if (fac_nbr_som <= 4) {

      /* La face est un triangle ou un quadrangle */
      /* =========================================*/

      pos_fac_dec = vec_def_fac_dec->pos_tab[ifac] ;

      vec_def_fac_dec->pos_tab[ifac + 1] = pos_fac_dec + 1 ;

      vec_def_fac_dec->val_tab[pos_fac_dec - 1] = cpt_fac_new + 1 ;


      /* Recopie de la définition de l'ancienne face */
      /*  dans la nouvelle liste des faces           */
      /*---------------------------------------------*/

      pos_fac_new = vec_def_fac_new->pos_tab[cpt_fac_new] ;

      vec_def_fac_new->pos_tab[cpt_fac_new + 1] = pos_fac_new + fac_nbr_som ;

      for (isom = 0 ; isom < fac_nbr_som ; isom++)
        vec_def_fac_new->val_tab[pos_fac_new - 1 + isom]
          = vec_def_fac->val_tab[isom_inf + isom] ;

      cpt_val_new += fac_nbr_som ;

      cpt_fac_new++ ;


    }
    else { /* if (fac_nbr_are > 4) */

      /* La face est un polygone et doit être découpée */
      /*===============================================*/


      /* Découpage effectif */
      /*--------------------*/

      fac_nbr_tria = ecs_loc_vec_def__dec_fac_tria(fac_nbr_som,
                                                   liste_som_tria,
                                                   liste_prec,
                                                   liste_suiv,
                                                   liste_are_def,
                                                   liste_voisinage_are,
                                                   liste_are_loc_delaunay,
                                                   concave,
                                                   coord_som) ;

      if (fac_nbr_tria != (fac_nbr_som - 2)) {

        bft_printf(_("  Vertices of face %d:\n"),
                   ifac + 1) ;

        for (isom = 0 ; isom < fac_nbr_som ; isom++) {

          num_som  = vec_def_fac->val_tab[isom_inf + isom] ;

          icoo_inf = vec_def_som->pos_pas * (num_som - 1) ;

          bft_printf("  %10d  [% 10.5e, % 10.5e, % 10.5e]\n",
                     num_som, vec_def_som->val_tab[icoo_inf],
                     vec_def_som->val_tab[icoo_inf + 1],
                     vec_def_som->val_tab[icoo_inf + 2]) ;

        }

        bft_printf(_("\n  Projected vertex coordinates of face %d:\n"),
                   ifac + 1) ;

        for (isom = 0 ; isom < fac_nbr_som ; isom++)
          bft_printf("  %10d  [% 10.5e, % 10.5e]\n",
                     isom + 1, coord_som[isom][0], coord_som[isom][1]) ;

        ecs_warn() ;
        bft_printf(_("Problem in triangulation of face <%ld> ;\n"
                     "we have %d vertices, so we should have %d triangles,\n"
                     "but be obtain %d.\n"),
                   (long)(ifac + 1), fac_nbr_som, fac_nbr_som-2, fac_nbr_tria) ;

      }


      /* Définition de l'ancienne face en fonction des triangles */
      /*---------------------------------------------------------*/


      pos_fac_dec = vec_def_fac_dec->pos_tab[ifac] ;

      vec_def_fac_dec->pos_tab[ifac + 1] = pos_fac_dec + fac_nbr_tria ;


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


        vec_def_fac_dec->val_tab[pos_fac_dec - 1 + itri] = cpt_fac_new + 1 ;


        /* Définition du triangle en fonction des sommets */
        /*------------------------------------------------*/

        pos_fac_new = vec_def_fac_new->pos_tab[cpt_fac_new] ;

        vec_def_fac_new->pos_tab[cpt_fac_new + 1] = pos_fac_new + 3 ;

        cpt_val_new += 3 ;

        for (isom = 0 ; isom < 3 ; isom++)

          vec_def_fac_new->val_tab[pos_fac_new + isom - 1]
            = vec_def_fac->val_tab[isom_inf + liste_som_tria[(itri*3) + isom]] ;

        cpt_fac_new++ ;

      }

    }

  } /* Fin : boucle sur les faces */


  /* Libération des tableaux de travail */

  BFT_FREE(som_fac) ;

  BFT_FREE(coord_som) ;
  BFT_FREE(liste_prec) ;
  BFT_FREE(liste_suiv) ;
  BFT_FREE(concave) ;
  BFT_FREE(liste_som_tria) ;

  BFT_FREE(liste_are_def) ;
  BFT_FREE(liste_voisinage_are) ;
  BFT_FREE(liste_are_loc_delaunay) ;


  /* Finalisation et vérification des dimensions */

  assert (nbr_fac_new >= cpt_fac_new) ;

  if (cpt_fac_new < nbr_fac_new) {

    ecs_warn() ;

    bft_printf(_("Triangulation problem for some faces;\n"
                 "probably due to excessive warping of original faces.\n")) ;

    bft_printf(_("Non critical for post-processing.\n")) ;

    ecs_vec_int__redimensionne(vec_def_fac_dec, nbr_fac + 1, cpt_fac_new) ;

    pos_fac_new = vec_def_fac_new->pos_tab[cpt_fac_new] ;

    ecs_vec_int__redimensionne(vec_def_fac_new,
                               cpt_fac_new + 1,
                               pos_fac_new - 1) ;

  }

  vec_def_fac_new->pos_nbr = cpt_fac_new + 1 ;


  return vec_def_fac_dec ;

}


/*----------------------------------------------------------------------------
 *  Fonction qui renvoie le nombre de faces de type "polygone"
 *  (basée sur la définition de faces par plus de 4 arêtes)
 *----------------------------------------------------------------------------*/

size_t ecs_vec_def__nbr_fac_poly
(
 const ecs_vec_int_t  *const vec_def_fac
)
{

  ecs_int_t  nbr_fac ;

  ecs_int_t  ifac ;

  size_t     nbr_fac_poly ;


  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  nbr_fac_poly = 0 ;

  nbr_fac = vec_def_fac->pos_nbr - 1 ;

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

    if (vec_def_fac->pos_tab[ifac + 1] - vec_def_fac->pos_tab[ifac] > 4)

      nbr_fac_poly++ ;

  }

  return nbr_fac_poly ;

}


/*----------------------------------------------------------------------------
 *  Fonction qui construit un tableau de booleens conforme à une liste
 *   de sous-éléments
 *  Un sous-élément est à `ECS_TRUE'
 *   s'il intervient dans la définition des éléments
 *----------------------------------------------------------------------------*/

void ecs_vec_def__cree_masque
(
 const ecs_tab_bool_t        bool_sselt_select,
 const ecs_vec_int_t  *const vec_def_elt
)
{

  size_t        nbr_elt ;
  size_t        num_sselt ;
  size_t        pos_inf ;
  size_t        pos_sup ;

  size_t        ielt ;
  size_t        ipos ;


  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  nbr_elt   = ecs_vec_int__ret_pos_nbr(vec_def_elt) - 1 ;


  for (ielt = 0 ; ielt < nbr_elt ; ielt++) {

    pos_inf = vec_def_elt->pos_tab[ielt    ] - 1 ;
    pos_sup = vec_def_elt->pos_tab[ielt + 1] - 1 ;


    for (ipos = pos_inf ; ipos < pos_sup ; ipos++) {

      num_sselt = ECS_ABS(vec_def_elt->val_tab[ipos]) - 1 ;

      if (bool_sselt_select.val[num_sselt] == ECS_FALSE) {

        /* Le sous-élément n'a pas encore été pris en compte */

        bool_sselt_select.val[num_sselt] = ECS_TRUE ;

      }
      /* else : rien à faire (le sous-élément a déjà été pris en compte) */

    } /* Fin boucle sur les positions de l'élément */

  } /* Fin boucle sur les éléments */


}


/*----------------------------------------------------------------------------
 *  Fusion des sommets confondus appartenant à des mêmes arêtes.
 *  La connectivité des arêtes est mise à jour.
 *----------------------------------------------------------------------------*/

ecs_vec_int_t * ecs_vec_def__nettoie_som_are
(
 ecs_vec_real_t  *const vec_def_som,
 ecs_vec_int_t   *const vec_def_are
)
{

  size_t     cpt_fusion ;

  size_t     icoo ;
  size_t     idim ;
  size_t     iare ;

  size_t     ipos_0 ;
  size_t     ipos_1 ;
  size_t     isom_0 ;
  size_t     isom_1 ;

  size_t     nbr_are ;
  size_t     nbr_som ;

  double      delta_coo ;
  double      lng_are_min ;
  double      lng_are_2 ;
  double      lng_min_2 ;

  ecs_tab_int_t  tab_equiv_som ;

  ecs_vec_int_t  *vec_som_old_new ;


  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  idim = vec_def_som->pos_pas ;

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


  lng_are_min = 1.0e-15 ;
  if (getenv("ECS_DEF_LNG_ARE_MIN") != NULL) {
    lng_are_min = atof(getenv("ECS_DEF_LNG_ARE_MIN")) ;
  }
  lng_min_2 = lng_are_min * ECS_ABS(lng_are_min) ;


  /* Marquage des sommets */
  /*----------------------*/

  tab_equiv_som = ecs_tab_int__cree_init(nbr_som, -1) ;

  cpt_fusion = 0 ;


  /* Boucle sur les arêtes pour détecter les sommets confondus */
  /*-----------------------------------------------------------*/

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

    ipos_0 = (vec_def_are->val_tab[iare*2    ] - 1) * idim ;
    ipos_1 = (vec_def_are->val_tab[iare*2 + 1] - 1) * idim ;

    lng_are_2 = 0.0 ;

    for (icoo = 0 ; icoo < idim ; icoo++) {

      delta_coo =   vec_def_som->val_tab[ipos_1 + icoo]
                  - vec_def_som->val_tab[ipos_0 + icoo] ;

      lng_are_2 += delta_coo * delta_coo ;

    }

    /* Sommets confondus détectés */

    if (lng_are_2 < lng_min_2) {

      isom_0 = vec_def_are->val_tab[iare*2    ] - 1 ;
      isom_1 = vec_def_are->val_tab[iare*2 + 1] - 1 ;

      ecs_vec_def__maj_equiv_som(isom_0,
                                 isom_1,
                                 &tab_equiv_som) ;

      cpt_fusion += 1 ;

    }

  }


  /* Fusion si nécessaire */

  if (cpt_fusion > 0) {

    bft_printf(_("\nMesh verification:\n\n"
                 "  %d vertices belong to edges of length less than %g\n"
                 "  and will be merged; (this tolerance may be modified\n"
                 "  using the ECS_DEF_LNG_ARE_MIN environment variable).\n"),
               (int)cpt_fusion, lng_are_min) ;

    vec_som_old_new = ecs_vec_def__fusion_som(vec_def_are,
                                              vec_def_som,
                                              &tab_equiv_som,
                                              NULL,
                                              NULL) ;

    bft_printf("\n") ;

  }
  else
    vec_som_old_new = NULL ;


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

  BFT_FREE(tab_equiv_som.val) ;

  return vec_som_old_new ;

}


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

ecs_vec_int_t  * ecs_vec_def__nettoie_are
(
 ecs_vec_int_t   *const vec_def_are
)
{

  ecs_int_t  nbr_are_old ;
  ecs_int_t  nbr_are_new ;

  ecs_int_t  cpt_are ;
  ecs_int_t  ind_are ;
  ecs_int_t  ind_som_0 ;
  ecs_int_t  ind_som_1 ;

  ecs_int_t  cpt_pos_new ;
  ecs_int_t  ind_are_new ;

  ecs_int_t  * renum_are ;

  ecs_vec_int_t  *vec_are_old_new ;


  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  /* Initialisations */

  nbr_are_old = vec_def_are->pos_nbr - 1 ;


  /* Boucle sur les arêtes (renumérotation) */
  /* -------------------------------------- */

  BFT_MALLOC(renum_are, nbr_are_old, ecs_int_t) ;

  cpt_are = 0 ;

  for (ind_are = 0 ; ind_are < nbr_are_old ; 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 ;



  /* Si on n'a pas d'arêtes dégénérées, on peut sortir */

  if (nbr_are_new == nbr_are_old) {

    BFT_FREE(renum_are) ;

    return NULL ;

  }

  bft_printf(_("\nMesh verification:\n\n"
               "  Removal of %d degenerate edges:\n"
               "    Initial number of edges           : %10d\n"
               "    Number of edges after processing  : %10d\n"),
             (nbr_are_old - nbr_are_new), nbr_are_old, nbr_are_new) ;


  /* 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 */
  /* ------------------------------- */

  vec_are_old_new = ecs_vec_int__alloue(nbr_are_old + 1,
                                        nbr_are_old) ;

  cpt_pos_new = 0 ;

  for (ind_are = 0 ; ind_are < nbr_are_old ; ind_are++) {

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

    ind_are_new = renum_are[ind_are] ;

    if (ind_are_new != -1)
        vec_are_old_new->val_tab[cpt_pos_new++] = ind_are_new + 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) ;

  return vec_are_old_new ;


}


/*----------------------------------------------------------------------------
 *  Fonction qui supprime les éventuelles faces dégénérées
 *----------------------------------------------------------------------------*/

ecs_vec_int_t  * ecs_vec_def__nettoie_fac
(
 ecs_vec_int_t   *const vec_def_fac
)
{

  ecs_int_t  nbr_fac_old ;
  ecs_int_t  nbr_fac_new ;

  ecs_int_t  ind_fac ;

  ecs_int_t  cpt_fac ;
  ecs_int_t  cpt_pos_new ;
  ecs_int_t  cpt_val ;
  ecs_int_t  ind_fac_new ;
  ecs_int_t  ind_val ;
  ecs_int_t  ind_val_deb ;
  ecs_int_t  ind_val_fin ;

  ecs_int_t  * renum_fac ;

  ecs_vec_int_t  *vec_fac_old_new ;


  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  /* Initialisations */

  nbr_fac_old = vec_def_fac->pos_nbr - 1 ;


  /* Boucle sur les faces (renumérotation) */
  /* ------------------------------------- */

  BFT_MALLOC(renum_fac, nbr_fac_old, ecs_int_t) ;

  cpt_fac = 0 ;
  cpt_val = 0 ;

  for (ind_fac = 0 ; ind_fac < nbr_fac_old ; ind_fac++) {

    ind_val_deb = vec_def_fac->pos_tab[ind_fac    ] - 1 ;
    ind_val_fin = vec_def_fac->pos_tab[ind_fac + 1] - 1 ;

    if (ind_val_fin - ind_val_deb > 2) {

      for (ind_val = ind_val_deb ; ind_val < ind_val_fin ; ind_val++)
        vec_def_fac->val_tab[cpt_val++] = vec_def_fac->val_tab[ind_val] ;

      renum_fac[ind_fac] = cpt_fac++;

      vec_def_fac->pos_tab[cpt_fac] = cpt_val + 1 ;

    }
    else

      renum_fac[ind_fac] = -1 ;

  }

  nbr_fac_new = cpt_fac ;



  /* Si on n'a pas de faces dégénérées, on peut sortir */

  if (nbr_fac_new == nbr_fac_old) {

    BFT_FREE(renum_fac) ;

    return NULL ;

  }

  bft_printf(_("\nMesh verification:\n\n"
               "  Removal of %d degenerate faces:\n"
               "    Initial number of faces           : %10d\n"
               "    Number of faces after processing  : %10d\n"),
             (nbr_fac_old - nbr_fac_new), nbr_fac_old, nbr_fac_new) ;


  /* On redimensionne le tableau des faces */

  ecs_vec_int__redimensionne(vec_def_fac,
                             nbr_fac_new + 1,
                             cpt_val) ;


  /* Boucle sur les anciennes faces */
  /* ------------------------------- */

  vec_fac_old_new = ecs_vec_int__alloue(nbr_fac_old + 1,
                                        nbr_fac_old) ;

  cpt_pos_new = 0 ;

  for (ind_fac = 0 ; ind_fac < nbr_fac_old ; ind_fac++) {

    vec_fac_old_new->pos_tab[ind_fac] =  cpt_pos_new + 1 ;

    ind_fac_new = renum_fac[ind_fac] ;

    if (ind_fac_new != -1)
        vec_fac_old_new->val_tab[cpt_pos_new++] = ind_fac_new + 1 ;

  }

  vec_fac_old_new->pos_tab[nbr_fac_old] =  cpt_pos_new + 1 ;


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

  BFT_FREE(renum_fac) ;

  ecs_vec_int__redimensionne(vec_fac_old_new,
                             nbr_fac_old + 1,
                             cpt_pos_new) ;

  return vec_fac_old_new ;


}


/*----------------------------------------------------------------------------
 *  Fonction qui met à jour le tableau des fusions de sommets
 *
 *  Remarque : le tableau d'équivalence (fusion) des sommets est construit de
 *             manière à ce qu'à un sommet ne comportant pas d'équivalent
 *             où de plus grand indice parmi ses équivalents corresponde la
 *             valeur -1, alors qu'un un sommet possédant des équivalents de
 *             plus grand indice corresponde le plus petit indice parmi ces
 *             équivalents (ce qui constitue une sorte de liste chaînée).
 *----------------------------------------------------------------------------*/

void ecs_vec_def__maj_equiv_som
(
 size_t                ind_som_0,
 size_t                ind_som_1,
 ecs_tab_int_t  *const tab_equiv_som
)
{

  ecs_int_t  ind_som_max ;
  ecs_int_t  ind_som_min ;
  ecs_int_t  ind_som_tmp ;

  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/

  ind_som_min = ECS_MIN(ind_som_0, ind_som_1) ;
  ind_som_max = ECS_MAX(ind_som_0, ind_som_1) ;

  /* On fusionne deux listes chaînées */
  /*----------------------------------*/

  while (   ind_som_max != ind_som_min
         && tab_equiv_som->val[ind_som_min] != ind_som_max) {

    /*
      On parcourt la liste chaînée correspondant à ind_som_min jusqu'au
      point d'insertion (si pas de point suivant dans la chaine,
      tab_equiv_som->val[ind_som_min] = -1).
    */

    while (   tab_equiv_som->val[ind_som_min] != -1
           && tab_equiv_som->val[ind_som_min] < ind_som_max)
      ind_som_min = tab_equiv_som->val[ind_som_min] ;

    ind_som_tmp = tab_equiv_som->val[ind_som_min] ;

    /*
      Si l'on est en fin de chaîne, on branche la liste chaînée correspondant
      à ind_som_max à la suite de celle correspondant à ind_som_min.
      Sinon, on doit reboucler
    */

    tab_equiv_som->val[ind_som_min] = ind_som_max ;

    if (ind_som_tmp != -1) {

      ind_som_min = ind_som_max ;
      ind_som_max = ind_som_tmp ;

    }

  }

}


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

/*----------------------------------------------------------------------------
 *  Fonction qui réalise l'ordonnancement des éléments
 *   en fonction de leur définition
 *----------------------------------------------------------------------------*/

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

  ecs_vec_real_t * vec_elt_ord ;  /* Vecteur ordonné des éléments */
  ecs_tab_int_t    vect_transf ;  /* Vecteur de transformation */

  size_t         nbr_elt ;        /* Nombre d'éléments */
  size_t         nbr_val ;
  size_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 ;
  nbr_val = ecs_vec_real__ret_val_nbr(*this_vec_def_real) ;

  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 ;

  vec_elt_ord = ecs_vec_real__alloue(nbr_elt + 1,
                                     nbr_val) ;

  ecs_vec_real__trie_et_renvoie(*this_vec_def_real,
                                vec_elt_ord,
                                &vect_transf,
                                tolerance) ;


  ecs_vec_real__detruit(*this_vec_def_real) ;


  /* Affectation de la nouvelle liste ordonnée */
  /*-------------------------------------------*/

  *this_vec_def_real = vec_elt_ord ;


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

  return vect_transf ;

}


/*----------------------------------------------------------------------------
 *  Fonction réalisant la liste compactée des éléments
 *         à partir de la liste ordonnée  des éléments
 *
 *  Le compactage consiste à supprimer les éléments géométriquement identiques
 *----------------------------------------------------------------------------*/

static ecs_vec_real_t * ecs_loc_vec_real_def__compacte
(
 const ecs_vec_real_t *const this_vec_real_ord,  /* --> Structure à compacter */
       ecs_tab_int_t  *const vect_transf         /* <-> Vecteur de transf.    */
)
{

  ecs_vec_real_t *vec_real_cpct ; /* Structure compactée des éléments */

  ecs_int_t    *vect_transf_cpct_val ;
  ecs_int_t    *vect_transf_val ;

  size_t        nbr_elt ;         /* Nombre d'éléments du vecteur à compacter */
  size_t        nbr_def ;         /* Nombre de définitions de l'ensemble des
                                     éléments du vecteur à compacter */

  ecs_real_t   *def_elt_ref ;

  ecs_int_t     nbr_elt_cpct ; /* Nombre d'éléments du vecteur compacte */
  ecs_int_t     cpt_val_cpct ; /* Compteur sur les valeurs des éléments
                                  du vecteur compacté */

  size_t        ielt ;         /* Indice sur les éléments à compacter */
  size_t        ipos ;         /* Indice sur les positions des définitions */
  size_t        ipos_ref ;

  ecs_bool_t      bool_diff ;


  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  /*=================*/
  /* Initialisations */
  /*=================*/


  nbr_elt = this_vec_real_ord->pos_nbr - 1 ;
  nbr_def = ecs_vec_real__ret_val_nbr(this_vec_real_ord) ;

  vec_real_cpct = ecs_vec_real__alloue(nbr_elt + 1,
                                       nbr_def) ;

  cpt_val_cpct = 0 ;


  BFT_MALLOC(vect_transf_cpct_val, nbr_elt, ecs_int_t) ;

  BFT_MALLOC(def_elt_ref, this_vec_real_ord->pos_pas, ecs_real_t) ;


  /*==============================================*/
  /* Repérage des éléments identiques             */
  /*  par comparaison de deux éléments successifs */
  /*==============================================*/


  /*---------------------------------------------*/
  /* Initialisation pour le tout premier élément */
  /*---------------------------------------------*/

  /* élément servant de référence au 1er tour */

  ipos_ref  = 0 ;

  for (ipos = 0 ; ipos < this_vec_real_ord->pos_pas ; ipos++) {
    def_elt_ref[ipos_ref++] = this_vec_real_ord->val_tab[ipos] ;
    vec_real_cpct->val_tab[cpt_val_cpct++] = this_vec_real_ord->val_tab[ipos];

  }

  vect_transf_cpct_val[0] = 0 ;

  nbr_elt_cpct = 1 ;


  /*--------------------------------*/
  /* Boucle sur les autres éléments */
  /*--------------------------------*/

  for (ielt = 1 ; ielt < nbr_elt ; ielt++) {

    ipos = this_vec_real_ord->pos_pas * ielt ;
    ipos_ref = 0 ;
    while (   ipos < this_vec_real_ord->pos_pas * (ielt + 1)
           && ipos_ref < this_vec_real_ord->pos_pas
           && ECS_ABS(  this_vec_real_ord->val_tab[ipos]
                      - def_elt_ref[ipos_ref]) < 2 * ECS_REAL_PRECISION) {
      ipos++ ;
      ipos_ref++ ;
    }

    if (ipos_ref == this_vec_real_ord->pos_pas) {

      /* Les deux éléments consécutifs dans la liste sont identiques */
      /*-------------------------------------------------------------*/

      vect_transf_cpct_val[ielt] = nbr_elt_cpct - 1 ;

      bool_diff = ECS_FALSE ;

    }
    else /* les deux éléments sont différents */

      bool_diff = ECS_TRUE ;


    if (bool_diff == ECS_TRUE) { /* Si les deux éléments sont différents */


      /* L'élément fait partie de la liste compactée des éléments */
      /*----------------------------------------------------------*/

      /* L'élément qui vient d'être compare       */
      /*  devient l'élément de référence          */
      /*  pour la comparaison du prochain élément */

      ipos_ref  = 0 ;

      for (ipos = this_vec_real_ord->pos_pas * ielt ;
           ipos < this_vec_real_ord->pos_pas * (ielt + 1) ;
           ipos++) {
        def_elt_ref[ipos_ref++] = this_vec_real_ord->val_tab[ipos] ;
        vec_real_cpct->val_tab[cpt_val_cpct++]
          = this_vec_real_ord->val_tab[ipos];
      }

      vect_transf_cpct_val[ielt] = nbr_elt_cpct ;

      nbr_elt_cpct++ ;


    } /* Fin : si les deux éléments sont différents */


  } /* Fin : boucle sur les autres éléments */


  BFT_FREE(def_elt_ref) ;


  vec_real_cpct->pos_nbr = nbr_elt_cpct + 1 ;
  vec_real_cpct->pos_pas = this_vec_real_ord->pos_pas ;


  BFT_REALLOC(vec_real_cpct->val_tab, cpt_val_cpct, ecs_real_t) ;


  BFT_MALLOC(vect_transf_val, nbr_elt, ecs_int_t) ;

  for (ielt = 0 ; ielt < nbr_elt ; ielt++)
    vect_transf_val[ielt] = vect_transf->val[ielt] ;


  for (ielt = 0 ; ielt < nbr_elt ; ielt++)
    vect_transf->val[ielt] = vect_transf_cpct_val[vect_transf_val[ielt]] ;


  BFT_FREE(vect_transf_val) ;
  BFT_FREE(vect_transf_cpct_val) ;


  /* Renvoi de la structure compactée des éléments */
  /*-----------------------------------------------*/

  return vec_real_cpct ;

}


/*----------------------------------------------------------------------------
 *  Fonction qui projette une face sur un plan parallèle à la face.
 *  Ce plan est ensuite assimilé au plan (Oxy).
 *----------------------------------------------------------------------------*/

static void ecs_loc_vec_def__plan_fac
(
 const ecs_int_t           nbr_som_fac,
       ecs_point_t  *const coo_som_fac
)
{

  ecs_int_t    icoo ;
  ecs_int_t    isom ;
  ecs_int_t    itri ;

  ecs_real_t   cost ;
  ecs_real_t   sint ;

  ecs_real_t   coo_tmp ;

  ecs_real_t   vect1[3] ;
  ecs_real_t   vect2[3] ;

  ecs_real_t   prod_vect[3] ;

  ecs_point_t  barycentre_fac ;
  ecs_point_t  normale_fac ;

  ecs_point_t *coo_som_fac_tmp ;


  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  /* Calcul des coordonnées du barycentre B du polygone P */
  /*======================================================*/

  for (icoo = 0 ; 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 = 0 ; icoo < ECS_DIM_3 ; icoo++)
    normale_fac[icoo] = 0. ;


  /* Calcul de la normale */
  /*======================*/

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

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

      vect1[icoo] = coo_som_fac[itri    ][icoo] - barycentre_fac[icoo] ;

      if (itri < nbr_som_fac - 1)
        vect2[icoo] = coo_som_fac[itri + 1][icoo] - barycentre_fac[icoo] ;
      else
        vect2[icoo] = coo_som_fac[0       ][icoo] - barycentre_fac[icoo] ;

    }

    normale_fac[0] += vect1[1] * vect2[2] - vect2[1] * vect1[2] ;
    normale_fac[1] += vect2[0] * vect1[2] - vect1[0] * vect2[2] ;
    normale_fac[2] += vect1[0] * vect2[1] - vect2[0] * vect1[1] ;

  }


  /* Projection dans un plan parallèle à la face */
  /*=============================================*/

  /* On ramène l'origine au centre de gravité de la fac */

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


  if (ECS_ABS(normale_fac[0]) > 1.e-12 || ECS_ABS(normale_fac[1]) > 1.e-12) {

    /* Première rotation d'axe (Oz) et d'angle (Ox, proj normale sur Oxy) */

    BFT_MALLOC(coo_som_fac_tmp, nbr_som_fac, ecs_point_t) ;

    vect1[0] = 1. ;
    vect1[1] = 0. ;
    vect1[2] = 0. ;

    vect2[0] = normale_fac[0] ;
    vect2[1] = normale_fac[1] ;
    vect2[2] = 0. ;

    ECS_LOC_PRODUIT_VECTORIEL(prod_vect, vect1, vect2) ;

    cost = ECS_LOC_PRODUIT_SCALAIRE(vect1, vect2) / ECS_LOC_MODULE(vect2) ;

    if (prod_vect[2] > 0.)
      sint =  ECS_LOC_MODULE(prod_vect) / ECS_LOC_MODULE(vect2) ;
    else
      sint = -ECS_LOC_MODULE(prod_vect) / ECS_LOC_MODULE(vect2) ;


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

      coo_som_fac_tmp[isom][0] =
         cost*coo_som_fac[isom][0] + sint*coo_som_fac[isom][1] ;
      coo_som_fac_tmp[isom][1] =
        -sint*coo_som_fac[isom][0] + cost*coo_som_fac[isom][1] ;
      coo_som_fac_tmp[isom][2] = coo_som_fac[isom][2] ;

    }


    /* Deuxième rotation d'axe (Oy) et d'angle (Oz', proj normale sur Ox'z) */

    vect1[0] =  0. ;
    vect1[1] =  0. ;
    vect1[2] =  1. ;

    vect2[0] =
      sqrt(normale_fac[0]*normale_fac[0] + normale_fac[1]*normale_fac[1]) ;
    vect2[1] = 0. ;
    vect2[2] = normale_fac[2] ;

    ECS_LOC_PRODUIT_VECTORIEL(prod_vect, vect1, vect2) ;

    cost = ECS_LOC_PRODUIT_SCALAIRE(vect1, vect2) / ECS_LOC_MODULE(vect2) ;

    if (prod_vect[2] > 0.)
      sint =  ECS_LOC_MODULE(prod_vect) / ECS_LOC_MODULE(vect2) ;
    else
      sint = -ECS_LOC_MODULE(prod_vect) / ECS_LOC_MODULE(vect2) ;


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

      coo_som_fac[isom][0] =
         cost*coo_som_fac_tmp[isom][0] + sint*coo_som_fac_tmp[isom][2] ;
      coo_som_fac[isom][1] = coo_som_fac_tmp[isom][1] ;
      coo_som_fac[isom][2] = 0. ;

    }


    BFT_FREE(coo_som_fac_tmp) ;


  }
  else {

    /* On écrase seulement la coordonnée z du sommet, en intervertissant
       éventuellement les coordonnées dans le plan de projection (Oxy).  */

    if (normale_fac[2] > 0.)
      for (isom = 0 ; isom < nbr_som_fac ; isom++)
        coo_som_fac[isom][2] = 0. ;

    else
      for (isom = 0 ; isom < nbr_som_fac ; isom++) {
        coo_tmp = coo_som_fac[isom][0] ;
        coo_som_fac[isom][0] = coo_som_fac[isom][1] ;
        coo_som_fac[isom][1] = coo_tmp ;
        coo_som_fac[isom][2] = 0. ;
      }

  }


}


/*----------------------------------------------------------------------------
 * Fonction qui vérifie si un sommet d'un polygone (projeté en 2D) est
 *  convexe.
 *----------------------------------------------------------------------------*/

static ecs_bool_t ecs_loc_vec_def__test_convexe
(
 const ecs_int_t           som_prec,
 const ecs_int_t           som_cour,
 const ecs_int_t           som_suiv,
       ecs_point_t  *const coo_som_fac
)
{

  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  /* sin(theta) = (v1 x v2) / (|v1| x |v2|), d'où le signe du sinus est celui
     de la composante z du produit vectoriel (v1 et v2 étant dans le plan) */

  if (((coo_som_fac[som_cour][0] - coo_som_fac[som_prec][0]) *
       (coo_som_fac[som_suiv][1] - coo_som_fac[som_cour][1])  ) -
      ((coo_som_fac[som_suiv][0] - coo_som_fac[som_cour][0]) *
       (coo_som_fac[som_cour][1] - coo_som_fac[som_prec][1])  )   > 0.0)

    return ECS_TRUE ;

  else

    return ECS_FALSE ;

}


/*----------------------------------------------------------------------------
 * Fonction qui vérifie si un sommet d'un polygone (projeté en 2D) est
 *  une "oreille".
 *----------------------------------------------------------------------------*/

static ecs_bool_t ecs_loc_vec_def__test_oreille
(
 const ecs_int_t           isom,
 const ecs_int_t           nbr_som_fac,
 const ecs_int_t    *const liste_prec,
 const ecs_int_t    *const liste_suiv,
 const ecs_bool_t   *const concave,
       ecs_point_t  *const coo_som_fac,
 const ecs_real_t          epsilon
)
{

  ecs_int_t iloc ;
  ecs_int_t som_prec ;
  ecs_int_t som_suiv ;

  ecs_real_t surf_2 ;
  ecs_real_t xisopa ;
  ecs_real_t yisopa ;

  ecs_real_t vect1[2] ;
  ecs_real_t vect2[2] ;
  ecs_real_t vect3[2] ;


  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  /* Si aucun sommet n'est concave, on a une "oreille" */

  for (iloc = 0 ; iloc < nbr_som_fac && concave[iloc] == ECS_FALSE ; iloc++) ;

  if (iloc == nbr_som_fac)

    return ECS_TRUE ;


  /* Si isom est convexe */

  else {

    if (concave[isom] == ECS_FALSE) {

      /* On vérifie si le triangle forme par liste_prec[isom], isom, et
         liste_suiv[isom] contient un sommet concave */

      som_prec = liste_prec[isom] ;
      som_suiv = liste_suiv[isom] ;

      vect2[0] = coo_som_fac[isom]    [0] - coo_som_fac[som_prec][0] ;
      vect2[1] = coo_som_fac[isom]    [1] - coo_som_fac[som_prec][1] ;
      vect3[0] = coo_som_fac[som_suiv][0] - coo_som_fac[som_prec][0] ;
      vect3[1] = coo_som_fac[som_suiv][1] - coo_som_fac[som_prec][1] ;

      surf_2 = vect2[0] * vect3[1] - vect3[0] * vect2[1] ;

      for (iloc  = liste_suiv[som_suiv] ;
           iloc != som_prec ;
           iloc  = liste_suiv[iloc]) {

        if (concave[iloc] == ECS_TRUE) {

          vect1[0] = coo_som_fac[iloc][0] - coo_som_fac[som_prec][0] ;
          vect1[1] = coo_som_fac[iloc][1] - coo_som_fac[som_prec][1] ;

          xisopa = (vect1[0]*vect3[1] - vect1[1]*vect3[0]) / surf_2 ;
          yisopa = (vect2[0]*vect1[1] - vect2[1]*vect1[0]) / surf_2 ;

          if ((1.0 - xisopa - yisopa > - epsilon) &&
              (      xisopa          > - epsilon) &&
              (               yisopa > - epsilon)   )

            return ECS_FALSE ;

        }

      }

      return ECS_TRUE ;

    }
    else

      return ECS_FALSE ;

  }

}


/*----------------------------------------------------------------------------
 *  Fonction transformant une triangulation quelconque en triangulation
 *  de Delaunay grâce à un algorithme d'échange d'arêtes (méthode flip).
 *----------------------------------------------------------------------------*/

static void ecs_loc_vec_def__delaunay_flip
(
 const ecs_int_t           nbr_som_fac,
       ecs_int_t    *const liste_som_tria,
       ecs_int_t    *const liste_are_def,
       ecs_int_t    *const liste_voisinage_are,
       ecs_bool_t   *const liste_are_loc_delaunay,
       ecs_point_t  *const coo_som_fac
)
{

  ecs_int_t    itri ;
  ecs_int_t    iare ;
  ecs_int_t    isom ;

  ecs_int_t    ind_tri ;
  ecs_int_t    ind_are ;
  ecs_int_t    ind_som ;

  ecs_int_t    ind_tri_0 ;
  ecs_int_t    ind_tri_1 ;

  ecs_int_t    ind_are_loc ;
  ecs_int_t    ind_are_flip ;

  ecs_int_t    isom_0 ;
  ecs_int_t    isom_1 ;

  ecs_int_t    isom_loc[2] ;
  ecs_int_t    isom_flip[2] ;
  ecs_int_t    isom_prec[2] ;
  ecs_int_t    isom_suiv[2] ;

  ecs_int_t    isom_inf ;
  ecs_int_t    isom_sup ;

  ecs_bool_t   bool_retour_debut ;
  ecs_bool_t   bool_fac_delaunay ;
  ecs_bool_t   bool_are_loc_delaunay ;
  ecs_bool_t   bool_quadrilatere_convexe ;

  const ecs_int_t    nbr_are_fac = nbr_som_fac*(nbr_som_fac - 1)/2 ;
  const ecs_int_t    nbr_tri_fac = nbr_som_fac - 2 ;


  /*
    Il y a (nbr_som_fac*(nbr_som_fac - 1)/2) combinaisons d'arête possibles.
    Le numéro d une arête peut être donné par :
    -> SOMME(k=0,k<isom_inf) { nbr_som_fac - k } + (isom_sup - isom_inf)
    i.e. nbr_som_fac*isom_inf - isom_inf*(isom_inf+1)/2 + (isom_sup - isom_inf)
    L'indice sera le numéro de l'arête moins un.

    Attention à l'ordre des arguments ! arg[0]=isom_inf, arg[1]=isom_sup
  */

#define ECS_LOC_IND_ARE(isom_inf, isom_sup) \
( nbr_som_fac*isom_inf - isom_inf*(isom_inf+1)/2 + isom_sup-isom_inf - 1 )


  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  /*-----------------*/
  /* Initialisations */
  /*-----------------*/

  for (isom_0 = 0 ; isom_0 < nbr_som_fac ; isom_0++) {
    for (isom_1 = isom_0 + 1 ; isom_1 < nbr_som_fac ; isom_1++) {

      ind_are = ECS_LOC_IND_ARE(isom_0, isom_1) ;

      /* Définition des arêtes */

      liste_are_def[2*ind_are    ] = isom_0 ;
      liste_are_def[2*ind_are + 1] = isom_1 ;

      /*
        Liste des triangles voisins partageant une arête :
        - (-1, -1) : l'arête n'existe pas
        - ( x, -1) : arête de bord ayant pour voisin le triangle x
        - ( x,  y) : arête interne ayant pour voisins les triangles x et y
      */

      liste_voisinage_are[2*ind_are    ] = -1 ;
      liste_voisinage_are[2*ind_are + 1] = -1 ;


      /* On initialise le tableau de booléens indiquant si l'arête est
         localement de Delaunay */

      liste_are_loc_delaunay[ind_are] = ECS_TRUE ;


    }

  }


  /* Premier parcours des triangles pour connaître le voisinage initial,
     ainsi que la liste des arêtes qui ne sont pas localement Delaunay */

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

    for (isom = 0 ; isom < 3 ; isom++) {

      isom_0 = liste_som_tria[(itri*3) +   isom      ] ;
      isom_1 = liste_som_tria[(itri*3) + ((isom+1)%3)] ;

      isom_inf = ECS_MIN(isom_0, isom_1) ;
      isom_sup = ECS_MAX(isom_0, isom_1) ;

      ind_are = ECS_LOC_IND_ARE(isom_inf, isom_sup) ;

      /* On met à jour le voisinage des arêtes */

      if (liste_voisinage_are[2*ind_are] == -1)
        liste_voisinage_are[2*ind_are] = itri ;
      else
        liste_voisinage_are[2*ind_are + 1] = itri ;

      /* Si ce n'est pas une arête de bord : */

      if (   !(isom_sup == isom_inf + 1)
          && !(isom_inf == 0 && isom_sup == nbr_som_fac - 1))
        liste_are_loc_delaunay[ind_are] = ECS_FALSE ;

    }

  }


  /*-----------------------------------*/
  /* Algorithme de flip proprement dit */
  /*-----------------------------------*/

  iare = 0 ;

  bool_retour_debut = ECS_FALSE ;
  bool_fac_delaunay = ECS_FALSE ;


  while(bool_fac_delaunay == ECS_FALSE) {

    if (liste_are_loc_delaunay[iare] == ECS_FALSE) {


      liste_are_loc_delaunay[iare] = ECS_TRUE ;


      isom_0 = liste_are_def[2*iare] ;
      isom_1 = liste_are_def[2*iare + 1] ;


      /* Premier travail sur les sommets des triangles */
      /*-----------------------------------------------*/

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

        ind_tri = liste_voisinage_are[2*iare + itri] ;


        for (isom = 0 ; isom < 3 ; isom++) {

          ind_som = liste_som_tria[3*ind_tri + isom] ;

          /* On cherche les sommets opposés */

          if (ind_som != isom_0 && ind_som !=isom_1)
            isom_flip[itri] = ind_som ;


          /* On cherche les sommets précedent et suivant */

          if (   ind_som == isom_0
                && liste_som_tria[3*ind_tri + ((isom + 1)%3)] == isom_1) {

            isom_prec[0] = liste_som_tria[3*ind_tri + ((isom + 2)%3)] ;
            isom_suiv[1] = isom_prec[0] ;

          }
          else if (   ind_som == isom_1
                   && liste_som_tria[3*ind_tri + ((isom + 1)%3)] == isom_0) {

            isom_suiv[0] = liste_som_tria[3*ind_tri + ((isom + 2)%3)] ;
            isom_prec[1] = isom_suiv[0] ;

          }

        } /* Fin parcours des sommets */

      } /* Fin parcours des triangles */


      /* Test sur la convexité du quadrilatère considéré */
      /*-------------------------------------------------*/

      if (   ecs_loc_vec_def__test_convexe(isom_prec[0],
                                           isom_0,
                                           isom_suiv[0],
                                           coo_som_fac) == ECS_TRUE

          && ecs_loc_vec_def__test_convexe(isom_prec[1],
                                           isom_1,
                                           isom_suiv[1],
                                           coo_som_fac) == ECS_TRUE) {

        bool_quadrilatere_convexe = ECS_TRUE ;

      }
      else {

        bool_quadrilatere_convexe = ECS_FALSE ;

      }


      /* On vérifie si l'arête est localement de Delaunay */
      /*--------------------------------------------------*/

      bool_are_loc_delaunay =
        ecs_loc_vec_def__are_loc_delaunay(isom_0,
                                          isom_1,
                                          isom_flip[0],
                                          isom_flip[1],
                                          coo_som_fac) ;


      /* Si l'arête n'est pas localement de Delaunay ... */
      /*-------------------------------------------------*/

      if (   bool_are_loc_delaunay == ECS_FALSE
          && bool_quadrilatere_convexe == ECS_TRUE) {


        isom_inf = ECS_MIN(isom_flip[0], isom_flip[1]) ;
        isom_sup = ECS_MAX(isom_flip[0], isom_flip[1]) ;

        ind_are_flip = ECS_LOC_IND_ARE(isom_inf, isom_sup) ;


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

          ind_tri_0 = liste_voisinage_are[2*iare + itri] ;
          ind_tri_1 = liste_voisinage_are[2*iare + ((itri + 1)%2)] ;


          /* On va changer la définition des triangles adjacents */
          /*-----------------------------------------------------*/

          for (isom = 0 ; isom < 3 ; isom++) {

            ind_som = liste_som_tria[3*ind_tri_0 + isom] ;

            if (ind_som == liste_are_def[2*iare + ((itri + 1)%2)])
              liste_som_tria[3*ind_tri_0 + isom] = isom_flip[(itri + 1)%2] ;

          }


          /* On renumérote le triangle dans l'ordre croissant de ses sommets */
          /*-----------------------------------------------------------------*/

          ecs_loc_vec_def__trie_delaunay(liste_som_tria + 3*ind_tri_0) ;


          /* On met à jour le voisinage et les arêtes non localament Delaunay */
          /*------------------------------------------------------------------*/

          for (isom = 0 ; isom < 3 ; isom++) {

            isom_loc[0] = liste_som_tria[3*ind_tri_0 + isom] ;
            isom_loc[1] = liste_som_tria[3*ind_tri_0 + ((isom + 1)%3)] ;

            isom_inf = ECS_MIN(isom_loc[0], isom_loc[1]) ;
            isom_sup = ECS_MAX(isom_loc[0], isom_loc[1]) ;

            ind_are_loc = ECS_LOC_IND_ARE(isom_inf, isom_sup) ;

            if (ind_are_loc < iare)
              bool_retour_debut = ECS_TRUE ;


            /* Si ce n'est pas la nouvelle arête ... */

            if (ind_are_loc != ind_are_flip) {

              if (liste_voisinage_are[2*ind_are_loc] == ind_tri_1)
                liste_voisinage_are[2*ind_are_loc] = ind_tri_0 ;
              else if (liste_voisinage_are[2*ind_are_loc + 1] == ind_tri_1)
                liste_voisinage_are[2*ind_are_loc + 1] = ind_tri_0 ;

              /* ... et que ce n'est pas une arête de bord */

              if (liste_voisinage_are[2*ind_are_loc + 1] != -1)
                liste_are_loc_delaunay[ind_are_loc] = ECS_FALSE ;

            }


          }


        } /* Fin parcours sur les triangles */


        ind_tri_0 = liste_voisinage_are[2*iare] ;
        ind_tri_1 = liste_voisinage_are[2*iare + 1] ;

        liste_voisinage_are[2*ind_are_flip] = ind_tri_0 ;
        liste_voisinage_are[2*ind_are_flip + 1] = ind_tri_1 ;

        liste_voisinage_are[2*iare] = -1 ;
        liste_voisinage_are[2*iare + 1] = -1 ;


      } /* Fin si l'arête n'était effectivement pas localement de Delaunay */


    } /* Fin si l'arête était initialement supposée non localement de Delaunay */


    if (iare == nbr_are_fac - 1 && bool_retour_debut == ECS_TRUE) {
      bool_retour_debut = ECS_FALSE ;
      iare = 0 ;
    }
    else if (iare == nbr_are_fac - 1 && bool_retour_debut == ECS_FALSE)
      bool_fac_delaunay = ECS_TRUE ;
    else
      iare++ ;


  } /* Fin de l'algorithme de flip */


}


/*----------------------------------------------------------------------------
 *  Fonction renvoyant un booléen selon si l'arête est localement de
 *  Delaunay ou non. On calcule la puissance d'un point par rapport au
 *  cercle circonscrit du triangle formé par les trois autres (dont deux
 *  définissant la diagonale considérée).
 *----------------------------------------------------------------------------*/

static ecs_bool_t ecs_loc_vec_def__are_loc_delaunay
(
 const ecs_int_t          isom_are_0,
 const ecs_int_t          isom_are_1,
 const ecs_int_t          isom_flip_0,
 const ecs_int_t          isom_flip_1,
       ecs_point_t *const coo_som_fac
)
{

  ecs_real_t   A ;
  ecs_real_t   B ;

  ecs_real_t   lambda[4] ;

  ecs_real_t   delta ;

  ecs_real_t   centre_x ;
  ecs_real_t   centre_y ;

  ecs_real_t   rayon ;

  ecs_real_t   puiss_point ;


  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  lambda[0] = 2*(coo_som_fac[isom_are_1][0] - coo_som_fac[isom_are_0][0]) ;
  lambda[1] = 2*(coo_som_fac[isom_are_1][1] - coo_som_fac[isom_are_0][1]) ;

  lambda[2] = 2*(coo_som_fac[isom_flip_0][0] - coo_som_fac[isom_are_0][0]) ;
  lambda[3] = 2*(coo_som_fac[isom_flip_0][1] - coo_som_fac[isom_are_0][1]) ;

  delta = lambda[1]*lambda[2] - lambda[0]*lambda[3] ;


  /* Si le triangle considéré est plat, on change de diagonale de manière
     automatique pour éviter une division par zéro. */

  if (ECS_ABS(delta) < 1.e-12)
    return ECS_FALSE ;


  A =   coo_som_fac[isom_are_1][0]*coo_som_fac[isom_are_1][0]
      - coo_som_fac[isom_are_0][0]*coo_som_fac[isom_are_0][0]
    +   coo_som_fac[isom_are_1][1]*coo_som_fac[isom_are_1][1]
      - coo_som_fac[isom_are_0][1]*coo_som_fac[isom_are_0][1] ;

  B =   coo_som_fac[isom_flip_0][0]*coo_som_fac[isom_flip_0][0]
      - coo_som_fac[isom_are_0][0]*coo_som_fac[isom_are_0][0]
    +   coo_som_fac[isom_flip_0][1]*coo_som_fac[isom_flip_0][1]
      - coo_som_fac[isom_are_0][1]*coo_som_fac[isom_are_0][1] ;


  /* Coordonnées du centre et rayon du cercle passant par les sommets
     de la diagonale et un troisième sommet. */

  centre_x = (lambda[1]*B - lambda[3]*A)/delta ;
  centre_y = (lambda[2]*A - lambda[0]*B)/delta ;

  rayon = sqrt(  (centre_x - coo_som_fac[isom_are_0][0])*
                 (centre_x - coo_som_fac[isom_are_0][0])
               + (centre_y - coo_som_fac[isom_are_0][1])*
                 (centre_y - coo_som_fac[isom_are_0][1])) ;


  /* On calcule la puissance du quatrième sommet du quadrilatère par rapport
     au cercle */

  puiss_point =   (coo_som_fac[isom_flip_1][0] - centre_x)*
                  (coo_som_fac[isom_flip_1][0] - centre_x)
                + (coo_som_fac[isom_flip_1][1] - centre_y)*
                  (coo_som_fac[isom_flip_1][1] - centre_y)
                - rayon*rayon ;


  /* On garde un peu de marge pour les cas où il n'y a pas de gain visible
     si l'on faisait un changement de diagonale */

  if (puiss_point > -1.e-12)
    return ECS_TRUE ;
  else
    return ECS_FALSE ;


}


/*----------------------------------------------------------------------------
 *  Fonction triant les indices des sommets définissant un triangle
 *----------------------------------------------------------------------------*/

static void ecs_loc_vec_def__trie_delaunay
(
 ecs_int_t   *const liste_som_tria
)
{

  ecs_int_t   imin ;
  ecs_int_t   imoy ;
  ecs_int_t   imax ;


  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  /* 1ière étape */

  if (*liste_som_tria < *(liste_som_tria + 1)) {

    imin = *liste_som_tria ;
    imoy = *(liste_som_tria + 1) ;
    imax = imoy ;

  }
  else {

    imin = *(liste_som_tria + 1) ;
    imoy = *liste_som_tria ;
    imax = imoy ;

  }


  /* 2ième étape */

  if (*(liste_som_tria + 2) < imin) {

    imoy = imin ;
    imin = *(liste_som_tria + 2) ;

  }
  else if (*(liste_som_tria + 2) > imax) {

    imax = *(liste_som_tria + 2) ;

  }
  else {

    imoy = *(liste_som_tria + 2) ;

  }


  /* Réaffectation */

  *liste_som_tria       = imin ;
  *(liste_som_tria + 1) = imoy ;
  *(liste_som_tria + 2) = imax ;


}


/*----------------------------------------------------------------------------
 *  Fonction réalisant le découpage des faces polygonales en triangles
 *----------------------------------------------------------------------------*/

static ecs_int_t ecs_loc_vec_def__dec_fac_tria
(
 const ecs_int_t            nbr_som_fac,
       ecs_int_t     *const liste_som_tria,
       ecs_int_t     *const liste_prec,
       ecs_int_t     *const liste_suiv,
       ecs_int_t     *const liste_are_def,
       ecs_int_t     *const liste_voisinage_are,
       ecs_bool_t    *const liste_are_loc_delaunay,
       ecs_bool_t    *const concave,
       ecs_point_t   *const coord_som
)
{

  ecs_int_t isom ;

  ecs_int_t  nbr_tria_fac = 0 ;
  ecs_int_t  num_essai    = 0 ;

  ecs_real_t epsilon[] = {1.0e-1, 1.0e-2, 0.0, -1.0e-2, -1.0e-1} ;


  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/



  /* On détermine dans quel plan on travaille (coord_som est écrasé) */

  ecs_loc_vec_def__plan_fac(nbr_som_fac, coord_som) ;


  /* Initialisations */

  while ((nbr_tria_fac != (nbr_som_fac - 2)) &&
         num_essai < 5) {

    nbr_tria_fac = 0 ; /* remise à zéro si essai précédent infructueux */

    for (isom = 0 ; isom < nbr_som_fac ; isom++) {
      liste_prec[isom] = isom - 1;
      liste_suiv[isom] = isom + 1;
    }
    liste_prec[0] = nbr_som_fac - 1 ;
    liste_suiv[nbr_som_fac - 1] = 0 ;

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

      if (ecs_loc_vec_def__test_convexe(liste_prec[isom],
                                        isom,
                                        liste_suiv[isom],
                                        coord_som) == ECS_TRUE)
        concave[isom] = ECS_FALSE ;
      else
        concave[isom] = ECS_TRUE ;

    }

    isom = 2 ;

    while (isom != 0 && isom != nbr_som_fac) {

      if (ecs_loc_vec_def__test_oreille(liste_prec[isom],
                                        nbr_som_fac,
                                        liste_prec,
                                        liste_suiv,
                                        concave,
                                        coord_som,
                                        epsilon[num_essai]) == ECS_TRUE) {

        /* On ajoute un triangle de sommets liste_prec[liste_prec[isom]],
           liste_prec[isom], isom */

        liste_som_tria[nbr_tria_fac*3    ] = liste_prec[liste_prec[isom]] ;
        liste_som_tria[nbr_tria_fac*3 + 1] = liste_prec[isom] ;
        liste_som_tria[nbr_tria_fac*3 + 2] = isom ;

        nbr_tria_fac += 1 ;


        /* On coupe l'oreille correspondant à liste_prec[isom] */

        liste_prec[isom] = liste_prec[liste_prec[isom]] ;
        liste_suiv[liste_prec[isom]] = isom ;

        if ((concave[isom] == ECS_TRUE) &&
            (ecs_loc_vec_def__test_convexe(liste_prec[isom],
                                           isom,
                                           liste_suiv[isom],
                                           coord_som) == ECS_TRUE))
          concave[isom] = ECS_FALSE ;

        if ((concave[liste_prec[isom]] == ECS_TRUE) &&
            (ecs_loc_vec_def__test_convexe(liste_prec[liste_prec[isom]],
                                           liste_prec[isom],
                                           isom,
                                           coord_som) == ECS_TRUE))
          concave[liste_prec[isom]] = ECS_FALSE ;

        if (liste_prec[isom] == 0)
          isom = liste_suiv[isom] ;

      }
      else
        isom = liste_suiv[isom] ;

    }

    num_essai++ ;

  }


  /* Algorithme de flip permettant d'avoir une triangulation de Delaunay
     à partir d'une triangulation quelconque */

  if (nbr_tria_fac == nbr_som_fac - 2)

    ecs_loc_vec_def__delaunay_flip(nbr_som_fac,
                                   liste_som_tria,
                                   liste_are_def,
                                   liste_voisinage_are,
                                   liste_are_loc_delaunay,
                                   coord_som) ;

  return nbr_tria_fac ;


}


/*----------------------------------------------------------------------------
 *  Fonction qui transforme un tableau d'équivalence en liste chaînée simple
 *----------------------------------------------------------------------------*/

static void ecs_loc_vec_def__transf_equiv
(
 const ecs_vec_real_t    *const vec_def_som,    /* <-> Déf. sommets           */
       ecs_tab_int_t     *const tab_equiv_som   /* <-> Fusions de sommets     */
)
{

  size_t     ind_som ;
  ecs_int_t  ind_som_min ;
  ecs_int_t  ind_som_max ;
  ecs_int_t  ind_som_tmp ;

  ecs_tab_int_t    tab_equiv_prec ;


  /*xxxxxxxxxxxxxxxxxxxxxxxxxxx Instructions xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx*/


  tab_equiv_prec.nbr = vec_def_som->pos_nbr - 1 ;
  BFT_MALLOC(tab_equiv_prec.val, tab_equiv_prec.nbr, ecs_int_t) ;

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

    tab_equiv_prec.val[ind_som] = -1 ;


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

    if (tab_equiv_som->val[ind_som] != -1) {

      ind_som_min = ind_som ;
      ind_som_max = tab_equiv_som->val[ind_som] ;

      assert (ind_som_min < ind_som_max) ;

      while (   ind_som_min != ind_som_max
             && tab_equiv_prec.val[ind_som_max] != ind_som_min) {

        /*
          On parcourt la liste inverse correspondant à ind_som_max jusqu'au
          point d'insertion (si pas de point precedent dans la chaine,
          tab_equiv_prec.val[ind_som_max] = -1).
        */

        while (tab_equiv_prec.val[ind_som_max] > ind_som_min)
          ind_som_max = tab_equiv_prec.val[ind_som_max] ;

        ind_som_tmp = tab_equiv_prec.val[ind_som_max] ;

        /*
          Si l'on est en début de chaîne, on branche la liste inverse
          correspondant à ind_som_min au début de celle correspondant
          à ind_som_max. Sinon, on doit reboucler.
        */

        tab_equiv_prec.val[ind_som_max] = ind_som_min ;

        if (ind_som_tmp != -1) {

          ind_som_max = ind_som_min ;
          ind_som_min = ind_som_tmp ;

        }

      }

    }

  }

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

    if (tab_equiv_prec.val[ind_som] != -1)

      tab_equiv_som->val[tab_equiv_prec.val[ind_som]] = ind_som ;

  }

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

}


/*----------------------------------------------------------------------------
 *  Correction si nécessaire de l'orientation d'un quadrangle en connectivité
 *   nodale ; renvoie -1 si l'on ne parvient pas à corriger l'orientation
 *   ou que l'on demande une simple vérification (i.e. correc = ECS_FALSE),
 *   0 si l'orientation initiale est bonne, et 1 en cas de correction de
 *   l'orientation par une permutation de la connectivité locale.
 *----------------------------------------------------------------------------*/

static ecs_int_t ecs_loc_vec_def__orient_quad
(
 const ecs_real_t  coord[],
       ecs_int_t   connect[],
 const ecs_bool_t  correc
)
{
  ecs_int_t   isom_tmp ;
  ecs_int_t   itri ;
  ecs_int_t   nb_angle_obtus ;
  ecs_real_t  sgn ;
  ecs_real_t  normale[3] ;
  ecs_real_t  vect1[3] ;
  ecs_real_t  vect2[3] ;
  ecs_real_t  prod_vect[3] ;

  ecs_int_t   passage = -1 ;


#define ECS_LOC_INIT_VECT(vect, i, j) ( \
  vect[0] =   coord[((connect[j-1] - 1) * 3)    ]  \
            - coord[((connect[i-1] - 1) * 3)    ], \
  vect[1] =   coord[((connect[j-1] - 1) * 3) + 1]  \
            - coord[((connect[i-1] - 1) * 3) + 1], \
  vect[2] =   coord[((connect[j-1] - 1) * 3) + 2]  \
            - coord[((connect[i-1] - 1) * 3) + 2] )


  /* Calcul d'une direction normale approchée de la face
     (peut être entrante ou sortante si la face est "croisée") */

  ECS_LOC_INIT_VECT(vect1, 1, 2) ;
  ECS_LOC_INIT_VECT(vect2, 1, 4) ;

  ECS_LOC_PRODUIT_VECTORIEL(normale, vect1, vect2) ;


  /* Boucle sur les renumérotations possibles */

  do {

    passage += 1 ;

    nb_angle_obtus = 0 ;

    /* Initialisation */

    switch(passage) {

    case 0:
      break ;

    case 1:
      if (correc == ECS_FALSE)
        return -1 ;
      isom_tmp = connect[2] ;
      connect[2] = connect[3] ;
      connect[3] = isom_tmp ;
      break ;

    default:
      return -1 ;

    }


    /* Boucle sur les coins */

    /* On compte les angles obtus, qui devraient être au nombre de 2 sur
       une face "croisée" (et 0 sur une face bien définie et convexe
       1 sur une face bien définie non convexe, soit 3 en apparaence
       sur une face bien définie convexe si la non-convexité se trouve
       au niveau du sommet 1 et que l'on a donc calculé une normale
       "à l'envers"). */

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

      ECS_LOC_INIT_VECT(vect1, ((itri+2) % 4) + 1, ((itri+1) % 4) + 1) ;
      ECS_LOC_INIT_VECT(vect2, ( itri    % 4) + 1, ((itri+1) % 4) + 1) ;

      ECS_LOC_PRODUIT_VECTORIEL(prod_vect, vect1, vect2) ;

      /* Angle obtus si produit mixte < 0, aigu sinon. */

      sgn = ECS_LOC_PRODUIT_SCALAIRE(prod_vect, normale) ;

      if (sgn < 0.0)
        nb_angle_obtus += 1;

    }

  } while (nb_angle_obtus == 2) ;

  return (ECS_MIN(passage, 1)) ;

#undef ECS_LOC_INIT_VECT

}


/*----------------------------------------------------------------------------
 *  Correction si nécessaire de l'orientation d'un tétraèdre en connectivité
 *   nodale ; renvoie -1 si l'on ne parvient pas à corriger l'orientation
 *   ou que l'on demande une simple vérification (i.e. correc = ECS_FALSE),
 *   0 si l'orientation initiale est bonne, et 1 en cas de correction de
 *   l'orientation par une permutation de la connectivité locale.
 *----------------------------------------------------------------------------*/

static ecs_int_t ecs_loc_vec_def__orient_tetra
(
 const ecs_real_t  coord[],
       ecs_int_t   connect[],
 const ecs_bool_t  correc
)
{
  ecs_int_t   isom_tmp ;
  ecs_real_t  det ;
  ecs_real_t  vect12[3] ;
  ecs_real_t  vect13[3] ;
  ecs_real_t  vect14[3] ;

  ecs_int_t   passage = -1 ;


#define ECS_LOC_INIT_VECT(vect, i, j) ( \
  vect[0] =   coord[((connect[j-1] - 1) * 3)    ]  \
            - coord[((connect[i-1] - 1) * 3)    ], \
  vect[1] =   coord[((connect[j-1] - 1) * 3) + 1]  \
            - coord[((connect[i-1] - 1) * 3) + 1], \
  vect[2] =   coord[((connect[j-1] - 1) * 3) + 2]  \
            - coord[((connect[i-1] - 1) * 3) + 2] )

#define ECS_LOC_PERMUTE(i, j) ( \
  isom_tmp = connect[i-1],     \
  connect[i-1] = connect[j-1], \
  connect[j-1] = isom_tmp      )


  /* Boucle sur les renumérotations possibles */

  do {

    passage += 1 ;

    /* Initialisation */

    switch(passage) {

    case 0:
      break ;

    case 1:
      if (correc == ECS_FALSE)
        return -1 ;
      ECS_LOC_PERMUTE(2, 3) ;
      break ;

    default: /* Retour connectivité d'origine et sortie */
      ECS_LOC_PERMUTE(2, 3) ;
      return -1 ;

    }

    ECS_LOC_INIT_VECT(vect12, 1, 2) ;
    ECS_LOC_INIT_VECT(vect13, 1, 3) ;
    ECS_LOC_INIT_VECT(vect14, 1, 4) ;

    det = ECS_LOC_DETERMINANT(vect12, vect13, vect14) ;

  } while (det < 0.0) ;

  return (ECS_MIN(passage, 1)) ;

#undef ECS_LOC_INIT_VECT
#undef ECS_LOC_PERMUTE

}


/*----------------------------------------------------------------------------
 *  Correction si nécessaire de l'orientation d'une pyramide en connectivité
 *   nodale ; renvoie -1 si l'on ne parvient pas à corriger l'orientation
 *   ou que l'on demande une simple vérification (i.e. correc = ECS_FALSE),
 *   0 si l'orientation initiale est bonne, et 1 en cas de correction de
 *   l'orientation par une permutation de la connectivité locale.
 *
 *  Tous les cas ne sont pas détectés ou traités : on suppose que le
 *   sommet est toujours en position 5, mais que la base peut être
 *   parcourue dans l'ordre 1 4 3 2, 1 2 4 3, ou 1 3 4 2 au lieu de 1 2 3 4,
 *   dans quel cas on la réordonne.
 *----------------------------------------------------------------------------*/

static ecs_int_t ecs_loc_vec_def__orient_pyram
(
 const ecs_real_t  coord[],
       ecs_int_t   connect[],
 const ecs_bool_t  correc
)
{
  ecs_int_t   isom ;
  ecs_int_t   isom_tmp ;
  ecs_int_t   connect_tmp[8] ;
  ecs_int_t   ret_base ;
  ecs_real_t  det ;
  ecs_real_t  vect1[3] ;
  ecs_real_t  vect2[3] ;
  ecs_real_t  vect3[3] ;

  ecs_int_t   retval = 0 ;
  ecs_int_t   passage = -1 ;


#define ECS_LOC_INIT_VECT(vect, i, j) ( \
  vect[0] =   coord[((connect_tmp[j-1] - 1) * 3)    ]  \
            - coord[((connect_tmp[i-1] - 1) * 3)    ], \
  vect[1] =   coord[((connect_tmp[j-1] - 1) * 3) + 1]  \
            - coord[((connect_tmp[i-1] - 1) * 3) + 1], \
  vect[2] =   coord[((connect_tmp[j-1] - 1) * 3) + 2]  \
            - coord[((connect_tmp[i-1] - 1) * 3) + 2] )

#define ECS_LOC_PERMUTE(i, j) ( \
  isom_tmp = connect[i-1],             \
  connect_tmp[i-1] = connect_tmp[j-1], \
  connect_tmp[j-1] = isom_tmp         )


  for (isom = 0 ; isom < 5 ; isom++)
    connect_tmp[isom] = connect[isom] ;

  /* Vérification et correction éventuelle de l'orientation de la base */

  ret_base = ecs_loc_vec_def__orient_quad(coord,
                                          connect_tmp,
                                          correc) ;

  retval = ECS_MAX(ret_base, retval) ;

  if ((correc == ECS_FALSE && ret_base != 0) || ret_base < 0)
    return - 1 ;


  /* Boucle sur les renumérotations possibles */

  do {

    passage += 1 ;

    /* Initialisation */

    switch(passage) {

    case 0:
      break ;

    case 1:
      if (correc == ECS_FALSE)
        return -1 ;
      else
        retval = 1 ;
      ECS_LOC_PERMUTE(2, 4) ;
      break ;

    default: /* Retour connectivité d'origine et sortie */
      ECS_LOC_PERMUTE(2, 4) ;
      return -1 ;

    }

    ECS_LOC_INIT_VECT(vect1, 1, 2) ;
    ECS_LOC_INIT_VECT(vect2, 1, 4) ;
    ECS_LOC_INIT_VECT(vect3, 1, 5) ;

    det = ECS_LOC_DETERMINANT(vect1, vect2, vect3) ;

  } while (det < 0.0) ;

  if (retval > 0) {
    for (isom = 0 ; isom < 5 ; isom++)
      connect[isom] = connect_tmp[isom] ;
  }

  return retval ;

#undef ECS_LOC_INIT_VECT
#undef ECS_LOC_PERMUTE

}


/*----------------------------------------------------------------------------
 *  Correction si nécessaire de l'orientation d'un prisme en connectivité
 *   nodale ; renvoie -1 si l'on ne parvient pas à corriger l'orientation
 *   ou que l'on demande une simple vérification (i.e. correc = ECS_FALSE),
 *   0 si l'orientation initiale est bonne, et 1 en cas de correction de
 *   l'orientation par une permutation de la connectivité locale.
 *
 *  Tous les cas ne sont pas détectés ou traités : on suppose que les
 *   bases peuvent être parcourues dans l'ordre 1 3 2 et 4 6 5 au lieu
 *   de 1 2 3 et 4 5 6, dans quel cas on les réordonne.
 *----------------------------------------------------------------------------*/

static ecs_int_t ecs_loc_vec_def__orient_prism
(
 const ecs_real_t  coord[],
       ecs_int_t   connect[],
 const ecs_bool_t  correc
)
{
  ecs_int_t   idim ;
  ecs_int_t   isom_tmp ;
  ecs_real_t  pscal ;
  ecs_real_t  vect1[3] ;
  ecs_real_t  vect2[3] ;
  ecs_real_t  vect3[3] ;

  ecs_int_t   passage    = -1 ;


#define ECS_LOC_INIT_VECT(vect, i, j) ( \
  vect[0] =   coord[((connect[j-1] - 1) * 3)    ]  \
            - coord[((connect[i-1] - 1) * 3)    ], \
  vect[1] =   coord[((connect[j-1] - 1) * 3) + 1]  \
            - coord[((connect[i-1] - 1) * 3) + 1], \
  vect[2] =   coord[((connect[j-1] - 1) * 3) + 2]  \
            - coord[((connect[i-1] - 1) * 3) + 2] )

#define ECS_LOC_PERMUTE(i, j) ( \
  isom_tmp = connect[i-1],     \
  connect[i-1] = connect[j-1], \
  connect[j-1] = isom_tmp      )


  /* Boucle sur les renumérotations possibles */

  do {

    passage += 1 ;

    /* Initialisation */

    switch(passage) {

    case 0:
      break ;

    case 1:
      if (correc == ECS_FALSE)
        return -1 ;
      ECS_LOC_PERMUTE(2, 3) ;
      ECS_LOC_PERMUTE(5, 6) ;
      break ;

    default: /* Retour connectivité d'origine et sortie */
      ECS_LOC_PERMUTE(2, 3) ;
      ECS_LOC_PERMUTE(5, 6) ;
      return -1 ;

    }

    ECS_LOC_INIT_VECT(vect2, 1, 2) ;
    ECS_LOC_INIT_VECT(vect3, 1, 3) ;

    ECS_LOC_PRODUIT_VECTORIEL(vect1, vect2, vect3) ;

    for (idim = 0 ; idim < 3 ; idim++)
      vect2[idim] = (  coord[((connect[4-1] - 1) * 3) + idim]
                     + coord[((connect[5-1] - 1) * 3) + idim]
                     + coord[((connect[6-1] - 1) * 3) + idim]
                     - coord[((connect[1-1] - 1) * 3) + idim]
                     - coord[((connect[2-1] - 1) * 3) + idim]
                     - coord[((connect[3-1] - 1) * 3) + idim]) ;

    pscal = ECS_LOC_PRODUIT_SCALAIRE(vect1, vect2) ;

  } while (pscal < 0.0) ;

  return (ECS_MIN(passage, 1)) ;

#undef ECS_LOC_INIT_VECT
#undef ECS_LOC_PERMUTE

}


/*----------------------------------------------------------------------------
 *  Correction si nécessaire de l'orientation d'un hexaèdre en connectivité
 *   nodale ; renvoie -1 si l'on ne parvient pas à corriger l'orientation
 *   ou que l'on demande une simple vérification (i.e. correc = ECS_FALSE),
 *   0 si l'orientation initiale est bonne, et 1 en cas de correction de
 *   l'orientation par une permutation de la connectivité locale.
 *
 *  Tous les cas ne sont pas détectés ou traités : on suppose que les
 *   bases peuvent être parcourues dans l'ordre soit 1 4 3 2 et 5 8 7 6,
 *   soit 1 2 4 3 et 5 6 8 7, soit 1 3 4 2 et 5 7 8 6 au lieu
 *   de 1 2 3 4 et 5 6 7 8, dans quel cas on les réordonne.
 *----------------------------------------------------------------------------*/

static ecs_int_t ecs_loc_vec_def__orient_hexa
(
 const ecs_real_t  coord[],
       ecs_int_t   connect[],
 const ecs_bool_t  correc
)
{
  ecs_int_t   idim ;
  ecs_int_t   isom ;
  ecs_int_t   isom_tmp ;
  ecs_int_t   connect_tmp[8] ;
  ecs_int_t   ret_base_1 ;
  ecs_int_t   ret_base_2 ;
  ecs_real_t  pscal ;
  ecs_real_t  vect1[3] ;
  ecs_real_t  vect2[3] ;
  ecs_real_t  vect3[3] ;

  ecs_int_t   retval = 0 ;
  ecs_int_t   passage = -1 ;


#define ECS_LOC_INIT_VECT(vect, i, j) ( \
  vect[0] =   coord[((connect_tmp[j-1] - 1) * 3)    ]  \
            - coord[((connect_tmp[i-1] - 1) * 3)    ], \
  vect[1] =   coord[((connect_tmp[j-1] - 1) * 3) + 1]  \
            - coord[((connect_tmp[i-1] - 1) * 3) + 1], \
  vect[2] =   coord[((connect_tmp[j-1] - 1) * 3) + 2]  \
            - coord[((connect_tmp[i-1] - 1) * 3) + 2] )

#define ECS_LOC_PERMUTE(i, j) ( \
  isom_tmp = connect[i-1],             \
  connect_tmp[i-1] = connect_tmp[j-1], \
  connect_tmp[j-1] = isom_tmp         )


  for (isom = 0 ; isom < 8 ; isom++)
    connect_tmp[isom] = connect[isom] ;

  /* Vérification et correction éventuelle de l'orientation des bases */

  ret_base_1 = ecs_loc_vec_def__orient_quad(coord,
                                            connect_tmp,
                                            correc) ;

  if ((correc == ECS_FALSE && ret_base_1 != 0) || ret_base_1 < 0)
    return - 1 ;

  else if (ret_base_1 > 0) {
    ret_base_2 = ecs_loc_vec_def__orient_quad(coord,
                                              connect_tmp + 4,
                                              correc) ;
    if (ret_base_2 != ret_base_1)
      return - 1 ;
    else
      retval = 1 ;
  }


  /* Boucle sur les renumérotations possibles */

  do {

    passage += 1 ;

    /* Initialisation */

    switch(passage) {

    case 0:
      break ;

    case 1:
      if (correc == ECS_FALSE)
        return -1 ;
      else
        retval = 1 ;
      ECS_LOC_PERMUTE(2, 4) ;
      ECS_LOC_PERMUTE(6, 8) ;
      break ;

    default:
      return -1 ;
      break ;

    }

    ECS_LOC_INIT_VECT(vect2, 1, 2) ;
    ECS_LOC_INIT_VECT(vect3, 1, 4) ;

    ECS_LOC_PRODUIT_VECTORIEL(vect1, vect2, vect3) ;

    for (idim = 0 ; idim < 3 ; idim++)
      vect2[idim] = (  coord[((connect_tmp[5-1] - 1) * 3) + idim]
                     + coord[((connect_tmp[6-1] - 1) * 3) + idim]
                     + coord[((connect_tmp[7-1] - 1) * 3) + idim]
                     + coord[((connect_tmp[8-1] - 1) * 3) + idim]
                     - coord[((connect_tmp[1-1] - 1) * 3) + idim]
                     - coord[((connect_tmp[2-1] - 1) * 3) + idim]
                     - coord[((connect_tmp[3-1] - 1) * 3) + idim]
                     - coord[((connect_tmp[4-1] - 1) * 3) + idim]) * 0.25 ;

    pscal = ECS_LOC_PRODUIT_SCALAIRE(vect1, vect2) ;

  } while (pscal < 0.0) ;

  if (retval > 0) {
    for (isom = 0 ; isom < 8 ; isom++)
      connect[isom] = connect_tmp[isom] ;
  }

  /* Vérification et correction éventuelle de l'orientation des cotés */

  connect_tmp[1-1] = connect[2-1] ;
  connect_tmp[2-1] = connect[3-1] ;
  connect_tmp[3-1] = connect[7-1] ;
  connect_tmp[4-1] = connect[6-1] ;

  ret_base_1 = ecs_loc_vec_def__orient_quad(coord,
                                            connect_tmp,
                                            correc) ;

  if (ret_base_1 != 0)
    return - 1 ;

  connect_tmp[1-1] = connect[1-1] ;
  connect_tmp[2-1] = connect[2-1] ;
  connect_tmp[3-1] = connect[6-1] ;
  connect_tmp[4-1] = connect[5-1] ;

  ret_base_1 = ecs_loc_vec_def__orient_quad(coord,
                                            connect_tmp,
                                            correc) ;

  if (ret_base_1 != 0)
    return - 1 ;


  return retval ;

#undef ECS_LOC_INIT_VECT
#undef ECS_LOC_PERMUTE

}


/*----------------------------------------------------------------------------
 * Compute contribution of a polygonal face to a cell's volume
 *
 * parameters:
 *   n_face_vertices -->  number of vertices
 *   vtx_coord       -->  vertex coordinates
 *   face_vtx_lst    -->  "face -> vertices" connectivity list
 *
 * returns:
 *   contribution of face to cell volume using Stoke's theorem
 *
 *                          Pi+1
 *              *---------*                   B  : barycentre of the polygon
 *             / .       . \
 *            /   .     .   \                 Pi : vertices of the polygon
 *           /     .   .     \
 *          /       . .  Ti   \               Ti : triangle
 *         *.........B.........* Pi
 *     Pn-1 \       . .       /
 *           \     .   .     /
 *            \   .     .   /
 *             \ .   T0  . /
 *              *---------*
 *            P0
 *----------------------------------------------------------------------------*/

static double
_compute_face_vol_contrib(ecs_int_t         n_face_vertices,
                          const ecs_real_t  vtx_coord[],
                          const ecs_int_t   face_vtx_lst[])
{
  ecs_int_t   i, tri_id, vtx_id;
  ecs_real_t  tri_center[3], tri_norm[3];
  ecs_real_t  vect1[3], vect2[3];

  ecs_real_t  face_barycentre[3] = {0., 0., 0.};

  double inv3 = 1./3.;

  double vol_contrib = 0.;

  /* Compute barycentre (B) coordinates for the polygon (P) */

  for (vtx_id = 0; vtx_id < n_face_vertices; vtx_id++) {
    size_t coord_id = face_vtx_lst[vtx_id] - 1;
    for (i = 0; i < 3; i++)
      face_barycentre[i] += vtx_coord[coord_id*3 + i];
  }

  for (i = 0; i < 3; i++)
    face_barycentre[i] /= n_face_vertices;

  /* Loop on triangles of the face */

  for (tri_id = 0 ; tri_id < n_face_vertices ; tri_id++) {

    size_t coord_id_1 = face_vtx_lst[tri_id % n_face_vertices] - 1;
    size_t coord_id_2 = face_vtx_lst[(tri_id + 1)% n_face_vertices] - 1;

    for (i = 0; i < 3; i++) {
      vect1[i] = vtx_coord[coord_id_1*3 + i] - face_barycentre[i];
      vect2[i] = vtx_coord[coord_id_2*3 + i] - face_barycentre[i];
    }

    ECS_LOC_PRODUIT_VECTORIEL(tri_norm, vect1, vect2);

    for (i = 0; i < 3; i++)
      tri_norm[i] *= 0.5;

    /* Computation of the center of the triangle Ti */

    for (i = 0; i < 3; i++)
      tri_center[i] = (  face_barycentre[i]
                       + vtx_coord[coord_id_1*3 + i]
                       + vtx_coord[coord_id_2*3 + i]) * inv3;

    /* Contribution to cell volume using Stoke's formula */

    for (i = 0; i < 3; i++)
      vol_contrib += (tri_norm[i] * tri_center[i]);

  } /* End of loop on triangles of the face */

  return vol_contrib;
}


/*----------------------------------------------------------------------------
 *  Correction si nécessaire de l'orientation d'un polyedre en connectivité
 *   nodale ; renvoie -1 si l'on ne parvient pas à corriger l'orientation
 *   ou que l'on demande une simple vérification (i.e. correc = ECS_FALSE),
 *   0 si l'orientation initiale est bonne, et 1 en cas de correction de
 *   l'orientation par une permutation de la connectivité locale.
 *
 *  On suppose que toutes les faces polygonales sont bien définies, mais
 *   que certaines faces peuvent être définies avec une normale intérieure
 *   à la cellule et non pas extérieure. On suppose que la non-convexité
 *   du polyèdre est limitée.
 *----------------------------------------------------------------------------*/

static ecs_int_t ecs_loc_vec_def__orient_polyedre
(
 const ecs_real_t      coord[],
       ecs_int_t       connect[],
       ecs_int_t       size,
       ecs_bool_t      correc,
       ecs_tab_int_t  *face_index,
       ecs_tab_int_t  *face_marker,
       ecs_tab_int_t  *edges
)
{
  size_t      i, j, face_id, edge_id ;
  ecs_int_t   ind_premier, ind_dernier ;

  double      cell_vol = 0.;
  size_t      n_unmarked_faces = 0 ;
  size_t      n_faces = 0 ;
  size_t      n_edges = 0 ;
  ecs_int_t   retval = 0 ;

  /* Ensure working arrays are large enough */

  /* Mark faces */

  {
    ecs_int_t   premier_som  = -1 ;

    for (i = 0 ; i < (size_t)size ; i++) {

      if (face_index->nbr < n_faces + 1) {
        face_index->nbr = ECS_MAX(n_faces + 1, face_index->nbr*2);
        BFT_REALLOC(face_index->val, face_index->nbr, ecs_int_t);
      }

      if (premier_som == -1) {
        face_index->val[n_faces] = i;
        premier_som = connect[i] ;
        ind_premier = i ;
      }
      else if (connect[i] == premier_som) {
        n_faces += 1;
        premier_som = -1 ;
        ind_dernier = i ;
      }
    }

    face_index->val[n_faces] = size;

    /* face_marker: 0 initially, 1 if oriented, -1 if inverted */

    if (face_marker->nbr < n_faces) {
      face_marker->nbr = ECS_MAX(n_faces, face_marker->nbr*2);
      BFT_REALLOC(face_marker->val, face_marker->nbr, ecs_int_t);
    }

    for (face_id = 0; face_id < n_faces; face_id++)
      face_marker->val[face_id] = 0;
  }

  if (n_faces == 0)
    return 0;

  /* Extract edges and build edge-face connectivity */
  /*------------------------------------------------*/

  for (face_id = 0; face_id < n_faces; face_id++) {

    size_t start_id = face_index->val[face_id];
    size_t end_id = face_index->val[face_id + 1] - 1;

    /* Loop on unassociated edges */

    for (i = start_id; i < end_id; i++) {

      ecs_int_t v_id_1 = connect[i], v_id_2 = connect[i + 1];
      ecs_int_t is_match = 0;

      /* Compare with other edges */

      for (edge_id = 0; edge_id < n_edges; edge_id++) {

        ecs_int_t v_id_1_cmp = edges->val[edge_id*4];
        ecs_int_t v_id_2_cmp = edges->val[edge_id*4 + 1];

        if ((v_id_1 == v_id_2_cmp) && (v_id_2 == v_id_1_cmp))
          is_match = 1;

        else if ((v_id_1 == v_id_1_cmp) && (v_id_2 == v_id_2_cmp))
          is_match = -1;

        /* Each edge should be shared by exactly 2 faces */

        if (is_match != 0) {
          if (edges->val[edge_id*4 + 3] == 0) {
            edges->val[edge_id*4 + 3] = is_match * (face_id + 1);
            break;
          }
          else
            return -1;
        }
      }

      /* If edge is unmatched, add it */

      if (is_match == 0) {

        if (edges->nbr < (n_edges + 1)*4) {
          edges->nbr = ECS_MAX(edges->nbr*2, (n_edges + 1)*4);
          BFT_REALLOC(edges->val, edges->nbr, ecs_int_t);
        }

        edges->val[n_edges*4] = v_id_1;
        edges->val[n_edges*4 + 1] = v_id_2;
        edges->val[n_edges*4 + 2] = face_id + 1;
        edges->val[n_edges*4 + 3] = 0;

        n_edges += 1;
      }

    }

  }

  /* Check if each edge is associated with 2 faces
     (i.e., that cell is closed) */

  for (edge_id = 0; edge_id < n_edges; edge_id++) {
    if (edges->val[edge_id*4 + 3] == 0)
      return -1;
  }

  /* Now check if face orientations are compatible */
  /*-----------------------------------------------*/

  face_marker->val[0] = 1; /* First face defines reference */
  n_unmarked_faces = n_faces - 1;

  while (n_unmarked_faces > 0) {

    for (edge_id = 0; edge_id < n_edges; edge_id++) {

      ecs_int_t face_num_2 = edges->val[edge_id*4 + 3];
      ecs_int_t face_id_1 = edges->val[edge_id*4 + 2] - 1;
      ecs_int_t face_id_2 = ECS_ABS(face_num_2) - 1;

      if (   face_marker->val[face_id_1] == 0
          && face_marker->val[face_id_2] != 0) {
        if (face_num_2 > 0)
          face_marker->val[face_id_1] = face_marker->val[face_id_2];
        else
          face_marker->val[face_id_1] = - face_marker->val[face_id_2];
        n_unmarked_faces -= 1;
      }

      else if (   face_marker->val[face_id_1] != 0
               && face_marker->val[face_id_2] == 0) {
        if (face_num_2 > 0)
          face_marker->val[face_id_2] = face_marker->val[face_id_1];
        else
          face_marker->val[face_id_2] = - face_marker->val[face_id_1];
        n_unmarked_faces -= 1;
      }
    }

  }

  for (edge_id = 0; edge_id < n_edges; edge_id++) {

    ecs_int_t face_num_2 = edges->val[edge_id*4 + 3];
    ecs_int_t face_id_1 = edges->val[edge_id*4 + 2] - 1;
    ecs_int_t face_id_2 = ECS_ABS(face_num_2) - 1;

    ecs_int_t marker_product
      = face_marker->val[face_id_1]*face_marker->val[face_id_2];

    if (   (face_num_2 < 0 && marker_product > 0)
        || (face_num_2 > 0 && marker_product < 0))
      return -1;

  }

  /* At this stage, topology is correct; check for inside-out cell */
  /*---------------------------------------------------------------*/

  for (face_id = 0; face_id < n_faces; face_id++) {

    size_t start_id = face_index->val[face_id];
    size_t end_id = face_index->val[face_id + 1] - 1;

    double vol_contrib
      = _compute_face_vol_contrib(end_id - start_id,
                                  coord,
                                  connect + start_id);
    cell_vol += vol_contrib * (double)(face_marker->val[face_id]);
  }

  /* Invert orientation if cell is inside_out */

  if (cell_vol < 0.) {
    for (face_id = 0; face_id < n_faces; face_id++)
      face_marker->val[face_id] *= -1;
  }

  /* Now correct connectivity if required */
  /*--------------------------------------*/

  if (correc == ECS_TRUE) {

    for (face_id = 0; face_id < n_faces; face_id++) {

      if (face_marker->val[face_id] == -1) {

        for (i = face_index->val[face_id] + 1,
               j = face_index->val[face_id + 1] - 2;
             i < j;
             i++, j--) {
          ecs_int_t connect_tmp = connect[i];
          connect[i] = connect[j];
          connect[j] = connect_tmp;
        }

        retval = 1;

      }
    }
  }

  else {

    for (face_id = 0; face_id < n_faces; face_id++) {
      if (face_marker->val[face_id] == -1)
        retval = -1;
    }

  }

  return retval ;
}


Generated by  Doxygen 1.6.0   Back to index