/*============================================================================
*
*                    Code_Saturne version 1.3
*                    ------------------------
*
*
*     This file is part of the Code_Saturne Kernel, element of the
*     Code_Saturne CFD tool.
*
*     Copyright (C) 1998-2007 EDF S.A., France
*
*     contact: saturne-support@edf.fr
*
*     The Code_Saturne Kernel 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 Kernel 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 Kernel; if not, write to the
*     Free Software Foundation, Inc.,
*     51 Franklin St, Fifth Floor,
*     Boston, MA  02110-1301  USA
*
*============================================================================*/

/*============================================================================
 * Grandeurs associées à un maillage
 *============================================================================*/

/* includes système */

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <assert.h>


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

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


/*----------------------------------------------------------------------------
 *  Fichiers `include' locaux
 *----------------------------------------------------------------------------*/

#include "cs_base.h"
#include "cs_maillage.h"
#include "cs_maillage_connect.h"
#include "cs_parallel.h"
#include "cs_perio.h"


/*----------------------------------------------------------------------------
 *  Fichiers `include' associés au fichier courant
 *----------------------------------------------------------------------------*/

#include "cs_maillage_grd.h"


#ifdef __cplusplus
extern "C" {
#endif /* __cplusplus */


/*=============================================================================
 * Définitions de macros locales
 *============================================================================*/

  enum {X, Y, Z} ;

#define _CS_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 _CS_PRODUIT_SCALAIRE(vect1, vect2) \
  (vect1[X] * vect2[X] + vect1[Y] * vect2[Y] + vect1[Z] * vect2[Z])

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


/*============================================================================
 *  Variables globales statiques
 *============================================================================*/

/* Pointeur sur les grandeurs du maillage principal */
cs_maillage_grd_t  *cs_glob_maillage_grd = NULL;

/* Choix de l'algorithme de calcul des centres de gravité des cellules */
static int cs_glob_maillage_grd_cdg_cel = 0;


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

/*----------------------------------------------------------------------------
 * Calcul des grandeurs associées à des faces (intérieures ou de bord)
 *
 *                          Pi+1
 *              *---------*                   B  : barycentre du polygone P
 *             / .       . \
 *            /   .     .   \                 Pi : sommet     du polygone P
 *           /     .   .     \
 *          /       . .  Ti   \               Ti : triangle
 *         *.........B.........* Pi
 *     Pn-1 \       . .       /
 *           \     .   .     /
 *            \   .     .   /
 *             \ .   T0  . /
 *              *---------*
 *            P0
 *----------------------------------------------------------------------------*/

static void _cs_maillage_grd_calc_faces
(
 const cs_int_t            dim,             /* --> Dimension de l'espace      */
 const cs_int_t            nbr_fac,         /* --> Nombre de faces            */
 const cs_real_t           coo_som[],       /* --> Coordonnées des sommets    */
 const cs_int_t            pos_fac_som[],   /* --> Positions de départ dans
                                                   connect. faces -> sommets  */
 const cs_int_t            val_fac_som[],   /* --> Connect. faces -> sommets  */
       cs_real_t           cdg_fac[],       /* <-- CDG faces                  */
       cs_real_t           surf_fac[]       /* <-- Normales faces             */
);


/*----------------------------------------------------------------------------
 * Calcul des centres des cellules C                                          *
 *                                                                            *
 * à partir de leurs n sommets S(i) où i=0,n-1                                *
 *                                                                            *
 *                                                                            *
 * -->      1    n-1   -->                                                    *
 * OB(C) = ---  Somme  OSi                                                    *
 *          n    i=0                                                          *
 *----------------------------------------------------------------------------*/

static void _cs_maillage_grd_cdg_cel_som
(
 const cs_maillage_t  *const maillage,      /* --> Maillage associé           */
       cs_real_t             cen_cel[]      /* <-- Barycentres des cellules   */
);


/*----------------------------------------------------------------------------
 * Calcul des centres G des cellules C                                        *
 *                                                                            *
 * à partir de leurs n faces F(i) où i=0,n-1                                  *
 *                                                                            *
 *                                                                            *
 *           n-1                                                              *
 *          Somme  Surf(Fi) G(Fi)                                             *
 *           i=0                                                              *
 *  G(C) = -----------------------                                            *
 *           n-1                                                              *
 *          Somme  Surf(Fi)                                                   *
 *           i=0                                                              *
 *----------------------------------------------------------------------------*/

static void _cs_maillage_grd_cdg_cel_fac
(
 const cs_maillage_t  *const maillage,   /* --> Maillage associé              */
 const cs_real_t             surf_fac[], /* --> Surfaces des faces internes   */
 const cs_real_t             cdg_fac[],  /* --> Centres de gravité des
                                                faces internes                */
 const cs_real_t             surf_fbr[], /* --> Surfaces des faces de bord    */
 const cs_real_t             cdg_fbr[],  /* --> Centres de gravité des
                                                faces de bord                 */
       cs_real_t             cen_cel[]   /* <-- Centres de gravité des
                                                cellules                      */
);


/*----------------------------------------------------------------------------
 * Calcul du volume des cellules C à partir de leurs n faces F(i)             *
 * et de leur centre de gravité G(Fi) où i=0,n-1                              *
 *                                                                            *
 *         1    n-1                                                           *
 *  G(C) = - .  Somme  Surf(Fi) G(Fi)                                         *
 *         3    i = 0                                                         *
 *----------------------------------------------------------------------------*/

static void _cs_maillage_grd_vol_cel
(
 const cs_maillage_t  *const maillage,   /* --> Maillage associé              */
 const cs_real_t             surf_fac[], /* --> Surfaces des faces internes   */
 const cs_real_t             cdg_fac[],  /* --> Centres de gravité des
                                                faces internes                */
 const cs_real_t             surf_fbr[], /* --> Surfaces des faces de bord    */
 const cs_real_t             cdg_fbr[],  /* --> Centres de gravité des
                                                faces de bord                 */
       cs_real_t             vol_cel[]   /* <-- Volume des cellules           */
);


/*============================================================================
 *  Fonctions publiques pour API Fortran
 *============================================================================*/

/*----------------------------------------------------------------------------
 * Interrogation ou modification du choix de l'algorithme de calcul
 * des centres de gravité des cellules ;
 * Au retour, IOPT contient la valeur (1 à 2) correspondant à l'algorithme
 * choisi
 *
 * Interface Fortran :
 *
 * SUBROUTINE ALGCEN (IOPT)
 * *****************
 *
 * INTEGER          IOPT        : <-> : Choix de l'algorithme
 *                                      < 0 : interrogation
 *                                        0 : calcul basé sur les centres et
 *                                            surfaces des faces (par défaut)
 *                                        1 : calcul basé sur les sommets
 *----------------------------------------------------------------------------*/

void CS_PROCF (algcen, ALGCEN)
(
 cs_int_t  *const iopt
)
{
  int iopt_ret = cs_maillage_grd_choix_cdg_cel((int)(*iopt));

  *iopt = iopt_ret;
}


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

/*----------------------------------------------------------------------------
 * Interrogation ou modification du choix de l'algorithme de calcul
 * des centres de gravité des cellules ;
 *  < 0 : interrogation
 *    0 : calcul basé sur les centres et surfaces des faces (choix par défaut)
 *    1 : calcul basé sur les sommets
 * Renvoie la valeur (1 à 2) correspondant à l'algorithme choisi
 *----------------------------------------------------------------------------*/

int cs_maillage_grd_choix_cdg_cel
(
 const int choix_algo
)
{
  if (choix_algo > 1)
    bft_error(__FILE__, __LINE__,0,
              _("L'indicateur de choix d'algorithme de calcul du centre "
                "de gravité\n"
                "des cellules peut prendre les valeurs suivantes :\n"
                "  0 : calcul basé sur les centres et surfaces des faces\n"
                "  1 : calcul basé sur les sommets\n"
                "et non %d."), cs_glob_maillage_grd_cdg_cel);

  else if (choix_algo >= 0)
    cs_glob_maillage_grd_cdg_cel = choix_algo;

  return cs_glob_maillage_grd_cdg_cel;
}

/*----------------------------------------------------------------------------
 * Création d'une structure grandeurs maillage
 *----------------------------------------------------------------------------*/

cs_maillage_grd_t  * cs_maillage_grd_cree
(
 void
)
{
  cs_maillage_grd_t  * maillage_grd;

  BFT_MALLOC (maillage_grd, 1, cs_maillage_grd_t);

  maillage_grd->cen_cel = NULL;
  maillage_grd->vol_cel = NULL;
  maillage_grd->surf_fac = NULL;
  maillage_grd->surf_fbr = NULL;
  maillage_grd->cdg_fac = NULL;
  maillage_grd->cdg_fbr = NULL;

  return (maillage_grd);
}


/*----------------------------------------------------------------------------
 * Destruction d'une structure grandeurs maillage
 *----------------------------------------------------------------------------*/

cs_maillage_grd_t  * cs_maillage_grd_detruit
(
 cs_maillage_grd_t  * maillage_grd
)
{

  BFT_FREE(maillage_grd->cen_cel);
  BFT_FREE(maillage_grd->vol_cel);
  BFT_FREE(maillage_grd->surf_fac);
  BFT_FREE(maillage_grd->surf_fbr);
  BFT_FREE(maillage_grd->cdg_fac);
  BFT_FREE(maillage_grd->cdg_fbr);

  BFT_FREE(maillage_grd);

  return (maillage_grd);
}

/*----------------------------------------------------------------------------
 * Calcul des grandeurs associées à un maillage
 *----------------------------------------------------------------------------*/

void cs_maillage_grd_calc
(
 const cs_maillage_t      *const maillage,      /* -->  maillage utilisé      */
       cs_maillage_grd_t  *const maillage_grd   /* <-> grandeurs associées    */
)
{
  cs_int_t  dim = maillage->dim;
  cs_int_t  nbr_fac = maillage->nbr_fac;
  cs_int_t  nbr_fbr = maillage->nbr_fbr;
  cs_int_t  nbr_cel = maillage->nbr_cel;
  cs_int_t  nbr_cel_etendu = maillage->nbr_cel_etendu;


  /* Allocation des structures s'il ne s'agit pas d'une mise à jour */

  if (maillage_grd->surf_fac == NULL)
    BFT_MALLOC(maillage_grd->surf_fac, nbr_fac*dim, cs_real_t);

  if (maillage_grd->cdg_fac == NULL)
    BFT_MALLOC(maillage_grd->cdg_fac, nbr_fac*dim, cs_real_t);

  if (maillage_grd->surf_fbr == NULL)
    BFT_MALLOC(maillage_grd->surf_fbr, nbr_fbr*dim, cs_real_t);

  if (maillage_grd->cdg_fbr == NULL)
    BFT_MALLOC(maillage_grd->cdg_fbr, nbr_fbr*dim, cs_real_t);

  if (maillage_grd->cen_cel == NULL)
    BFT_MALLOC(maillage_grd->cen_cel, nbr_cel_etendu*dim, cs_real_t);

  if (maillage_grd->vol_cel == NULL)
    BFT_MALLOC(maillage_grd->vol_cel, nbr_cel_etendu, cs_real_t);


  /* Calcul des centres de gravité et des surfaces des faces internes */

  _cs_maillage_grd_calc_faces(dim,
                              nbr_fac,
                              maillage->coo_som,
                              maillage->pos_fac_som,
                              maillage->val_fac_som,
                              maillage_grd->cdg_fac,
                              maillage_grd->surf_fac);


  /* Calcul des centres de gravité et des surfaces des faces de bord */

  _cs_maillage_grd_calc_faces(dim,
                              nbr_fbr,
                              maillage->coo_som,
                              maillage->pos_fbr_som,
                              maillage->val_fbr_som,
                              maillage_grd->cdg_fbr,
                              maillage_grd->surf_fbr);


  /* Calcul des centres des cellules basé sur le barycentre des faces
     ou sur les sommets */

  switch (cs_glob_maillage_grd_cdg_cel) {

  case 0:
    _cs_maillage_grd_cdg_cel_fac(maillage,
                                 maillage_grd->surf_fac,
                                 maillage_grd->cdg_fac,
                                 maillage_grd->surf_fbr,
                                 maillage_grd->cdg_fbr,
                                 maillage_grd->cen_cel);
    break;

  case 1:
    _cs_maillage_grd_cdg_cel_som(maillage,
                                 maillage_grd->cen_cel);
    break;

  default:
    assert(0);

  }

  /* Calcul du volume des cellules */

  _cs_maillage_grd_vol_cel(maillage,
                           maillage_grd->surf_fac,
                           maillage_grd->cdg_fac,
                           maillage_grd->surf_fbr,
                           maillage_grd->cdg_fbr,
                           maillage_grd->vol_cel);

  if (cs_glob_base_nbr > 1) {

    cs_int_t  idim, icel;
    cs_real_t  *coord_tmp = NULL;

    /* Synchronisation des coordonnées des centres des cellules
       entre processeurs dans les zones halo */

    BFT_MALLOC(coord_tmp, nbr_cel_etendu, cs_real_t);

    for (idim = 0 ; idim < dim ; idim++) {
      for (icel = 0 ; icel < nbr_cel ; icel++)
        coord_tmp[icel] = maillage_grd->cen_cel[icel*dim + idim];
      for (icel = nbr_cel ; icel < nbr_cel_etendu ; icel++)
        coord_tmp[icel] = 0.;
      CS_PROCF (parcom, PARCOM) (coord_tmp);
      for (icel = nbr_cel ; icel < nbr_cel_etendu ; icel++)
        maillage_grd->cen_cel[icel*dim + idim] = coord_tmp[icel];
    }

    BFT_FREE(coord_tmp);

    CS_PROCF (pargve, PARGVE) (&dim, maillage_grd->cen_cel);

    /* Synchronisation du volume des cellules entre processeurs
       dans les zones halo */

    for (icel = nbr_cel ; icel < nbr_cel_etendu ; icel++)
      maillage_grd->vol_cel[icel] = 0.;

    CS_PROCF (parcom, PARCOM) (maillage_grd->vol_cel);

  }

  /* En cas de halo monoprocesseur (périodicité), mise à zéro du halo */

  else if (nbr_cel_etendu > nbr_cel) {

    cs_int_t  idim, icel;

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

      /* Traitement du volume des cellules */

      maillage_grd->vol_cel[icel] = 0.;

      /* Traitement des centres cellule */

      for (idim = 0 ; idim < dim ; idim++)
        maillage_grd->cen_cel[icel*dim + idim] = 0.;
    }

  }

  /* En cas de périodicité, mise à jour du halo */

  if (maillage->nbr_per > 0)
    cs_perio_sync_geo() ;

}


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

/*----------------------------------------------------------------------------
 * Calcul des grandeurs associées à des faces (intérieures ou de bord)
 *
 *                          Pi+1
 *              *---------*                   B  : barycentre du polygone P
 *             / .       . \
 *            /   .     .   \                 Pi : sommet     du polygone P
 *           /     .   .     \
 *          /       . .  Ti   \               Ti : triangle
 *         *.........B.........* Pi
 *     Pn-1 \       . .       /
 *           \     .   .     /
 *            \   .     .   /
 *             \ .   T0  . /
 *              *---------*
 *            P0
 *----------------------------------------------------------------------------*/

static void _cs_maillage_grd_calc_faces
(
 const cs_int_t            dim,             /* --> Dimension de l'espace      */
 const cs_int_t            nbr_fac,         /* --> Nombre de faces            */
 const cs_real_t           coo_som[],       /* --> Coordonnées des sommets    */
 const cs_int_t            pos_fac_som[],   /* --> Positions de départ dans
                                                   connect. faces -> sommets  */
 const cs_int_t            val_fac_som[],   /* --> Connect. faces -> sommets  */
       cs_real_t           cdg_fac[],       /* <-- CDG faces                  */
       cs_real_t           surf_fac[]       /* <-- Normales faces             */
)
{
  cs_int_t     ifac;
  cs_int_t     nbr_som_fac, nbr_max_som_fac;
  cs_int_t     icoo, icoo_inf;
  cs_int_t     isom, isom_inf, isom_sup;
  cs_int_t     itri;

  cs_point_t  *coo_som_fac;
  cs_point_t   barycentre_fac;
  cs_point_t   normale_fac;
  cs_point_t   centre_fac;
  cs_point_t   centre_tri;
  cs_point_t  *normale_tri;
  cs_point_t   vect1, vect2;

  cs_real_t   surface_fac;
  cs_real_t   surface_tri;
  cs_real_t   contrib_vol_fac;
  cs_real_t   contrib_vol_tri;
  cs_real_t   correc_cdg;


  /* Retour si l'on n'a pas les informations suffisantes
     (cas SolCom hors radiatif ou Enveloppe pré 1.2.d sans option "-n") */

  if (pos_fac_som == NULL && val_fac_som == NULL)
    return;


  /* Vérifications */
  /*---------------*/

  if (dim != 3)
    bft_error(__FILE__, __LINE__,0,
               _("Le calcul des grandeurs géométriques associées aux faces\n"
                 "ne sont implantés qu'en 3D."));

  assert(cdg_fac != NULL || nbr_fac == 0);
  assert(surf_fac != NULL || nbr_fac == 0);


  /* Comptages et allocations */
  /*--------------------------*/

  nbr_max_som_fac = 0;

  for (ifac = 0 ; ifac < nbr_fac ; ifac++) {
    nbr_som_fac = pos_fac_som[ifac + 1] - pos_fac_som[ifac];
    if (nbr_max_som_fac <= nbr_som_fac) nbr_max_som_fac = nbr_som_fac;
  }

  BFT_MALLOC(coo_som_fac, nbr_max_som_fac + 1, cs_point_t);
  BFT_MALLOC(normale_tri, nbr_max_som_fac    , cs_point_t);


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

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

    contrib_vol_tri = 0. ;
    surface_fac = 0.0 ;

    /*------------------------------------------------------------------------*/
    /* Définition du polygone P en fonction des sommets Pi de la face         */
    /*------------------------------------------------------------------------*/

    isom_inf = pos_fac_som[ifac]     - 1;
    isom_sup = pos_fac_som[ifac + 1] - 1;

    nbr_som_fac = 0;

    for (isom = isom_inf ; isom < isom_sup ; isom++) {

      icoo_inf = 3 * (val_fac_som[isom] - 1);

      for (icoo = X ; icoo < 3 ; icoo++)
        coo_som_fac[nbr_som_fac][icoo] = coo_som[icoo_inf + icoo];

      nbr_som_fac++;

    }

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

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

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

      barycentre_fac[icoo] = 0.0;

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

      barycentre_fac[icoo] /= nbr_som_fac ;

    }

    for (icoo = X ; icoo < 3 ; icoo++) {
      normale_fac[icoo] = 0.0;
      centre_fac [icoo] = 0.0;
    }


    /* Première boucle sur les triangles de la face (calcul des normales)     */
    /*========================================================================*/

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

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

      for (icoo = X ; icoo < 3 ; icoo++) {
        vect1[icoo] = coo_som_fac[itri    ][icoo] - barycentre_fac[icoo];
        vect2[icoo] = coo_som_fac[itri + 1][icoo] - barycentre_fac[icoo];
      }

      _CS_PRODUIT_VECTORIEL(normale_tri[itri], vect1, vect2);

      for (icoo = X ; icoo < 3 ; icoo++)
        normale_tri[itri][icoo] *= 0.5;

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

      for (icoo = X ; icoo < 3 ; icoo++)
        normale_fac[icoo] += normale_tri[itri][icoo];

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


    /* Seconde boucle sur les triangles de la face (pour le barycentre) */
    /*==================================================================*/

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

      /*----------------------------------------------------------------------*/
      /* Calcul du centre de gravite G(Ti) du triangle Ti :                   */
      /*                                                                      */
      /*  -->            -->  -->   -->                                       */
      /*  OG(Ti) = 1/3 ( OB + OPi + OPi+1 )                                   */
      /*                                                                      */
      /* Et de la contribution au volume du triangle Ti :                     */
      /*                                                                      */
      /*  -->    ->                                                           */
      /*  OG(Ti).N(Ti)                                                        */
      /*----------------------------------------------------------------------*/

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

        centre_tri[icoo] =   barycentre_fac[icoo]
                           + coo_som_fac[itri    ][icoo]
                           + coo_som_fac[itri + 1][icoo];

        centre_tri[icoo] /= 3.0;

        contrib_vol_tri += (centre_tri[icoo] * normale_tri[itri][icoo]);

      }

      /*----------------------------------------------------------------------*/
      /* Calcul de la surface du triangle Ti comme module de la normale       */
      /*                                                                      */
      /*               ->                                                     */
      /*  Surf(Ti) = | N(Ti) |                                                */
      /*----------------------------------------------------------------------*/

      surface_tri = _CS_MODULE(normale_tri[itri]);

      if (_CS_PRODUIT_SCALAIRE(normale_tri[itri], normale_fac) < 0.0)
        surface_tri *= -1.0;

      surface_fac += surface_tri;

      /*----------------------------------------------------------------------*/
      /*   n-1                                                                */
      /*  Somme  Surf(Ti) G(Ti)                                               */
      /*   i=0                                                                */
      /*----------------------------------------------------------------------*/

      for (icoo = X ; icoo < 3 ; icoo++)
        centre_fac[icoo] += surface_tri * centre_tri[icoo];

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


    /*------------------------------------------------------------------------*/
    /* Calcul du centre de gravite G(P) du polygone P :                       */
    /*                                                                        */
    /*           n-1                                                          */
    /*          Somme  Surf(Ti) G(Ti)                                         */
    /*           i=0                                                          */
    /*  G(P) = -----------------------                                        */
    /*           n-1                                                          */
    /*          Somme  Surf(Ti)                                               */
    /*           i=0                                                          */
    /*                                                                        */
    /* Calcul de la contribution au volume du polygone (avant correction) :   */
    /*                                                                        */
    /*  -->    ->                                                             */
    /*  OG(P).N(P)                                                            */
    /*------------------------------------------------------------------------*/

    contrib_vol_fac = 0.0;

    for (icoo = X ; icoo < 3 ; icoo++) {
      centre_fac[icoo] = centre_fac[icoo] / surface_fac;
      contrib_vol_fac += (centre_fac[icoo] * normale_fac[icoo]);
    }

    correc_cdg = (contrib_vol_tri - contrib_vol_fac)
                 / (surface_fac * surface_fac);

    for (icoo = X ; icoo < 3 ; icoo++)
      centre_fac[icoo] += correc_cdg * normale_fac[icoo];

    /* Stockage du résultat dans les structures appropriées */
    /*------------------------------------------------------*/

    for (icoo = X ; icoo < 3 ; icoo++) {
      cdg_fac [ifac*3 + icoo] = centre_fac[icoo];
      surf_fac[ifac*3 + icoo] = normale_fac[icoo];
    }

  } /* Fin : boucle sur les faces */

  BFT_FREE(normale_tri);
  BFT_FREE(coo_som_fac);

}


/*----------------------------------------------------------------------------
 * Calcul des centres G des cellules C                                        *
 *                                                                            *
 * à partir de leurs n faces F(i) où i=0,n-1                                  *
 *                                                                            *
 *                                                                            *
 *           n-1                                                              *
 *          Somme  Surf(Fi) G(Fi)                                             *
 *           i=0                                                              *
 *  G(C) = -----------------------                                            *
 *           n-1                                                              *
 *          Somme  Surf(Fi)                                                   *
 *           i=0                                                              *
 *----------------------------------------------------------------------------*/

static void _cs_maillage_grd_cdg_cel_fac
(
 const cs_maillage_t  *const maillage,   /* --> Maillage associé              */
 const cs_real_t             surf_fac[], /* --> Surfaces des faces internes   */
 const cs_real_t             cdg_fac[],  /* --> Centres de gravité des
                                                faces internes                */
 const cs_real_t             surf_fbr[], /* --> Surfaces des faces de bord    */
 const cs_real_t             cdg_fbr[],  /* --> Centres de gravité des
                                                faces de bord                 */
       cs_real_t             cen_cel[]   /* <-- Centres de gravité des
                                                cellules                      */
)
{
  cs_int_t   ifac, icoo, ival;
  cs_int_t   icel, icel1, icel2;

  cs_real_t  *surf_tot;
  cs_real_t   surf;
  cs_real_t   vect_loc_surf[3];


  /* Connectivité du maillage */

  const  cs_int_t   dim = maillage->dim;
  const  cs_int_t   nbr_fac = maillage->nbr_fac;
  const  cs_int_t   nbr_fbr = maillage->nbr_fbr;
  const  cs_int_t   nbr_cel = maillage->nbr_cel;
  const  cs_int_t   nbr_cel_etendu = maillage->nbr_cel_etendu;
  const  cs_int_t  *fac_cel = maillage->fac_cel;
  const  cs_int_t  *fbr_cel = maillage->fbr_cel;


  /* Retour si l'on n'a pas les informations suffisantes
     (cas SolCom hors radiatif ou Enveloppe pré 1.2.d sans option "-n") */

  if (maillage->val_fac_som == NULL && maillage->val_fbr_som == NULL)
    return;


  /* Vérifications */
  /* ------------- */

  if (dim != 3)
    bft_error(__FILE__, __LINE__,0,
               _("Le calcul des centres des cellules\n"
                 "ne sont implantés qu'en 3D."));

  assert(cen_cel != NULL);


  /* Initialisation */
  /* -------------- */

  BFT_MALLOC(surf_tot, nbr_cel_etendu, cs_real_t);

  for (ival = 0 ; ival < nbr_cel_etendu ; ival++) {
    surf_tot[ival] = 0.;

    for (icoo = 0; icoo < dim ; icoo++) {
      cen_cel[dim*ival + icoo] = 0.;
    }

  }


  /* ----------------------------- */
  /* Boucle sur les faces internes */
  /* ----------------------------- */

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

    /* ----------------------------------------------------------- */
    /* Pour chacune des deux cellules voisines à la face interne,  */
    /* on met à jour le numérateur de cen_cel et surf_tot          */
    /* ----------------------------------------------------------- */

    icel1 = fac_cel[2*ifac    ] - 1;
    icel2 = fac_cel[2*ifac + 1] - 1;

    /* Calcul de la surface de la face */

    for (icoo = 0 ; icoo < dim ; icoo++)
      vect_loc_surf[icoo] = surf_fac[dim*ifac + icoo];

    surf = _CS_MODULE(vect_loc_surf);

    surf_tot[icel1] += surf;
    surf_tot[icel2] += surf;

    /* Calcul du numérateur */

    for (icoo = 0 ; icoo < dim ; icoo++) {
      cen_cel[dim*icel1 + icoo] += cdg_fac[dim*ifac + icoo]*surf;
      cen_cel[dim*icel2 + icoo] += cdg_fac[dim*ifac + icoo]*surf;
    }

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


  /* ----------------------------- */
  /* Boucle sur les faces de bord  */
  /* ----------------------------- */

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

    /* ----------------------------------------------------------- */
    /* Pour chaque cellule voisine de la face bord,                */
    /* on met à jour le numérateur de cen_cel et surf_tot          */
    /* ----------------------------------------------------------- */

    icel1 = fbr_cel[ifac] - 1;

    /* Calcul de la surface de la face */

    for (icoo = 0 ; icoo < dim ; icoo++)
      vect_loc_surf[icoo] = surf_fbr[dim*ifac + icoo];

    surf = _CS_MODULE(vect_loc_surf);
    surf_tot[icel1] += surf;

    /* Calcul du numérateur */

    for (icoo = 0 ; icoo < dim ; icoo++)
      cen_cel[dim*icel1 + icoo] += cdg_fbr[dim*ifac + icoo]*surf;

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


  /* ------------------------------------------------------------------*/
  /* Boucle sur les cellules pour finaliser la détermination du centre */
  /* ------------------------------------------------------------------*/

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

    for (icoo = 0 ; icoo < dim ; icoo++ )
      cen_cel[icel*dim + icoo] /= surf_tot[icel];

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


  /* Libération de la mémoire des pointeurs */
  BFT_FREE(surf_tot);

}


/*----------------------------------------------------------------------------
 * Calcul des centres des cellules C                                          *
 *                                                                            *
 * à partir de leurs n sommets S(i) où i=0,n-1                                *
 *                                                                            *
 *                                                                            *
 * -->      1    n-1   -->                                                    *
 * OB(C) = ---  Somme  OSi                                                    *
 *          n    i=0                                                          *
 *----------------------------------------------------------------------------*/

static void _cs_maillage_grd_cdg_cel_som
(
 const cs_maillage_t  *const maillage,      /* --> Maillage associé           */
       cs_real_t             cen_cel[]      /* <-- Barycentres des cellules   */
)
{
  cs_int_t     icel, ifac, isom, icoo;
  cs_int_t     ipos_fac, ipos_som;
  cs_int_t     num_fac;

  cs_int_t     cpt_som_cel;

  cs_int_t    *ind_cel_som;

  cs_int_t    *pos_fac_som_loc;
  cs_int_t    *val_fac_som_loc;

  cs_int_t    *pos_cel_fac = NULL;
  cs_int_t    *val_cel_fac = NULL;


  /* Retour si l'on n'a pas les informations suffisantes
     (cas SolCom hors radiatif ou Enveloppe pré 1.2.d sans option "-n") */

  if (maillage->val_fac_som == NULL && maillage->val_fbr_som == NULL)
    return;


  /* Vérifications */
  /* ------------- */

  if (maillage->dim != 3)
    bft_error(__FILE__, __LINE__,0,
               _("Le calcul des centres des cellules\n"
                 "ne sont implantés qu'en 3D."));

  assert(cen_cel != NULL);


  /* Allocations et Initialisation */
  /* ----------------------------- */

  BFT_MALLOC(ind_cel_som, maillage->nbr_som, cs_int_t);

  for (isom = 0 ; isom < maillage->nbr_som ; isom++)
    ind_cel_som[isom] = -1;

  /* Extraction de la connectivité "cellules -> faces" */

  cs_maillage_ret_cel_fac(maillage,
                          maillage->nbr_cel,
                          NULL,
                          &pos_cel_fac,
                          &val_cel_fac);


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

  for (icel = 0 ; icel < maillage->nbr_cel ; icel++) {

    cpt_som_cel = 0;

    for (icoo = 0 ; icoo < 3 ; icoo++)
      cen_cel[icel*3 + icoo] = 0.;

    /* Boucle sur les faces de la cellule */

    for (ipos_fac = pos_cel_fac[icel] ;
         ipos_fac < pos_cel_fac[icel + 1] ;
         ipos_fac++) {

      num_fac = val_cel_fac[ipos_fac - 1];

      /* Face interne ou face de bord */

      if (num_fac > maillage->nbr_fbr) {
        ifac = num_fac - maillage->nbr_fbr - 1;
        pos_fac_som_loc = maillage->pos_fac_som;
        val_fac_som_loc = maillage->val_fac_som;
      }
      else {
        ifac = num_fac - 1;
        pos_fac_som_loc = maillage->pos_fbr_som;
        val_fac_som_loc = maillage->val_fbr_som;
      }

      /* Boucle sur les sommets de la face */

      for (ipos_som = pos_fac_som_loc[ifac] ;
           ipos_som < pos_fac_som_loc[ifac + 1] ;
           ipos_som++) {

        isom = val_fac_som_loc[ipos_som] - 1;

        if (ind_cel_som[isom] < icel) {
          for (icoo = 0 ; icoo < 3 ; icoo++)
            cen_cel[icel*3 + icoo] += maillage->coo_som[isom*3 + icoo];
          cpt_som_cel += 1;
          ind_cel_som[isom] = icel;
        }

      }

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

    for (icoo = 0 ; icoo < 3 ; icoo++)
      cen_cel[icel*3 + icoo] /= (double)cpt_som_cel;

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


  /* Libération de la mémoire allouée */

  BFT_FREE(ind_cel_som);

  BFT_FREE(pos_cel_fac);
  BFT_FREE(val_cel_fac);

}


/*----------------------------------------------------------------------------
 * Calcul du volume des cellules C à partir de leurs n faces F(i)             *
 * et de leur centre de gravité G(Fi) où i=0,n-1                              *
 *                                                                            *
 *         1    n-1                                                           *
 *  G(C) = - .  Somme  Surf(Fi) G(Fi)                                         *
 *         3    i = 0                                                         *
 *----------------------------------------------------------------------------*/

static void _cs_maillage_grd_vol_cel
(
 const cs_maillage_t  *const maillage,   /* --> Maillage associé              */
 const cs_real_t             surf_fac[], /* --> Surfaces des faces internes   */
 const cs_real_t             cdg_fac[],  /* --> Centres de gravité des
                                                faces internes                */
 const cs_real_t             surf_fbr[], /* --> Surfaces des faces de bord    */
 const cs_real_t             cdg_fbr[],  /* --> Centres de gravité des
                                                faces de bord                 */
       cs_real_t             vol_cel[]   /* <-- Volume des cellules           */
)
{
  cs_int_t  i_elt1, i_elt2;
  cs_int_t  i_fac, i_fbr, i_cel;
  cs_int_t  i_dim;

  cs_real_t  flux = 0;

  const cs_real_t un_tiers = 1.0/3.0;
  const cs_int_t dim = maillage->dim;

  /* Initialisation */
  /*----------------*/

  for (i_cel = 0; i_cel < maillage->nbr_cel_etendu; i_cel++)
    vol_cel[i_cel] = 0;

  /* Boucle sur les faces internes */
  /* ----------------------------- */

  for (i_fac = 0; i_fac < maillage->nbr_fac; i_fac++) {

    i_elt1 = maillage->fac_cel[2*i_fac] - 1;
    i_elt2 = maillage->fac_cel[2*i_fac + 1] - 1;

    flux = 0;
    for (i_dim = 0; i_dim < dim; i_dim++)
      flux += surf_fac[dim*i_fac + i_dim] * cdg_fac[dim*i_fac + i_dim];

    vol_cel[i_elt1] += flux;
    vol_cel[i_elt2] -= flux;

  }

  /* Boucle sur les faces de bord */
  /* ---------------------------- */

  for (i_fbr = 0; i_fbr < maillage->nbr_fbr; i_fbr++) {

    i_elt1 = maillage->fbr_cel[i_fbr] - 1;

    flux = 0;
    for (i_dim = 0; i_dim < dim; i_dim++)
      flux += surf_fbr[dim*i_fbr + i_dim] * cdg_fbr[dim*i_fbr + i_dim];

    vol_cel[i_elt1] += flux;

  }

  /* Calcul du volume */

  for (i_cel = 0; i_cel < maillage->nbr_cel; i_cel++)
    vol_cel[i_cel] *= un_tiers;

}


/*----------------------------------------------------------------------------*/

#if 0 /* Morceau de test bonne orientation faces */

  cs_int_t   idim, ifac, icel;
  cs_real_t  cdgfac[3];
  cs_real_t  cdgcel[3];
  cs_real_t  normale[3];
  cs_real_t  pscal;

  for (ifac = 0 ; ifac < maillage->nbr_fbr ; ifac++) {

    icel = maillage->fbr_cel[ifac] - 1;
    pscal = 0;

    for (idim = 0 ; idim < 3 ; idim++) {
      cdgcel[idim]  = cs_glob_maillage_grd->cen_cel[icel*3 + idim];
      cdgfac[idim]  = cs_glob_maillage_grd->cdg_fbr[ifac*3 + idim];
      normale[idim] = cs_glob_maillage_grd->surf_fbr[ifac*3 + idim];
      pscal += normale[idim] * (cdgfac[idim] - cdgcel[idim]);
    }

    if (pscal < 0.0)
      printf("num_fac_brd = %d, num_cel = %d, pscal = %f\n",
             ifac + 1, icel + 1, pscal);
  }

#endif


/*----------------------------------------------------------------------------*/

/* Suppression des macros locales */

#undef _CS_PRODUIT_VECTORIEL
#undef _CS_PRODUIT_SCALAIRE
#undef _CS_MODULE


#ifdef __cplusplus
}
#endif /* __cplusplus */
