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

/*============================================================================
 * Définitions, variables globales, et fonctions de base
 *============================================================================*/


/* includes système */

#include <assert.h>
#include <errno.h>
#include <stdio.h>
#include <stdarg.h>
#include <string.h>
#include <math.h>

/* Includes librairie BFT et FVM */

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

/* Includes librairie */

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

#include "cs_voiset.h"

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

/*============================================================================
 * Variables statiques locales
 *============================================================================*/

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

/*============================================================================
 * Prototypes de fonctions FORTRAN appelées
 *============================================================================*/

/*----------------------------------------------------------------------------
 *  Calcul des gradients par moindres carrés (standard ou avec voisinage étendu)
 *----------------------------------------------------------------------------*/

void CS_PROCF (gradmc, GRADMC)
(
 const cs_int_t   *const idbia0,      /* --> numéro 1ère case libre dans IA   */
 const cs_int_t   *const idbra0,      /* --> numéro 1ère case libre dans RA   */
 const cs_int_t   *const ncelet,      /* --> nombre de cellules étendu        */
 const cs_int_t   *const ncel,        /* --> nombre de cellules               */
 const cs_int_t   *const nfac,        /* --> nombre de faces internes         */
 const cs_int_t   *const nfabor,      /* --> nombre de faces de bord          */
 const cs_int_t   *const ncelbr,      /* --> nombre de cellules de bord       */
 const cs_int_t   *const nituse,      /* --> longueur du tableau ituser[]     */
 const cs_int_t   *const nrtuse,      /* --> longueur du tableau rtuser[]     */
 const cs_int_t   *const inc,         /* --> 0 ou 1 : incrément ou non        */
 const cs_int_t   *const iccocg,      /* --> 1 ou 0 : recalcul de COCG ou non */
 const cs_int_t   *const nswrgp,      /* --> >1 : avec reconstruction         */
 const cs_int_t   *const idimte,      /* --> 0, 1, 2 : variable
                                             scalaire, vecteur, tenseur       */
 const cs_int_t   *const itenso,      /* --> pour la périodicité de rotation  */
 const cs_int_t   *const iphydp,      /* --> prise en compte de P hydrostatiqu*/
 const cs_int_t   *const imrgra,      /* --> mode de calcul du gradient       */
 const cs_int_t   *const iwarnp,      /* --> niveau d'impression              */
 const cs_int_t   *const nfecra,      /* --> unité sortie standard            */
 const cs_real_t  *const epsrgp,      /* --> précision pour calcul du gradient
                                             itératif                         */
 const cs_real_t  *const extrap,      /* --> extrapolation des grad. au bord  */
 const cs_int_t          ifacel[],    /* --> liste des faces internes         */
 const cs_int_t          ifabor[],    /* --> liste des faces de bord          */
 const cs_int_t          icelbr[],    /* --> numéro des cellules ayant au
                                             moins une face de bord           */
 const cs_int_t          ipcvse[],    /* --> pour chaque cellule, position
                                             dans ielvse de la première cellule
                                             du voisinage étendu              */
 const cs_int_t          ielvse[],    /* --> numéro des cellules du voisinage
                                             étendu                           */
       cs_int_t          ituser[],    /* <-> tab. complémentaire utilisateur  */
       cs_int_t          ia[],        /* <-> macro-tableau entier             */
 const cs_real_t         volume[],    /* --> volumes des cellules             */
 const cs_real_t         surfac[],    /* --> surfaces des faces internes      */
 const cs_real_t         surfbo[],    /* --> surfaces des faces de bord       */
 const cs_real_t         surfbn[],    /* --> norme de surfbo                  */
 const cs_real_t         pond[],      /* --> pondération géométrique faces int*/
 const cs_real_t         dist[],      /* --> distance entre I' et J' faces int*/
 const cs_real_t         distbr[],    /* --> distance du cdgfbo à I' fac. brd */
 const cs_real_t         dijpf[],     /* --> vecteur I'J' aux faces internes  */
 const cs_real_t         diipb[],     /* --> vecteur II' aux faces de bord    */
 const cs_real_t         fextx[],     /* --> composantes de la force          */
 const cs_real_t         fexty[],     /*     extérieure générant la pression  */
 const cs_real_t         fextz[],     /*     hydrostatique                    */
 const cs_real_t         xyzcen[],    /* --> c.d.g. des cellules              */
 const cs_real_t         cdgfac[],    /* --> c.d.g. des faces internes        */
 const cs_real_t         cdgfbo[],    /* --> c.d.g. des faces de bord         */
 const cs_real_t         coefap[],    /* --> cond. limites aux faces de bord  */
 const cs_real_t         coefbp[],    /* --> cond. limites aux faces de bord  */
 const cs_real_t         pvar[],      /* --> variable dont on calcule le grad */
       cs_real_t         cocgb[],     /* <-> contribution à COCG des faces
                                             internes des cellules de bord    */
       cs_real_t         cocg[],      /* <->  contribution à COCG des faces
                                             internes des cellules de bord    */
       cs_real_t         rtuser[],    /* <-> tab. complémentaire utilisateur  */
       cs_real_t         dpdx[],      /* <-- composante x du gradient         */
       cs_real_t         dpdy[],      /* <-- composante y du gradient         */
       cs_real_t         dpdz[],      /* <-- composante z du gradient         */
       cs_real_t         bx[],        /*  -  tab. de travail local            */
       cs_real_t         by[],        /*  -  tab. de travail local            */
       cs_real_t         bz[],        /*  -  tab. de travail local            */
       cs_real_t         ra[]         /* <-> macro-tableau réel               */
 );


/*----------------------------------------------------------------------------
 *  Limitation des gradients (standard ou avec voisinage étendu)
 *----------------------------------------------------------------------------*/

void CS_PROCF (limgrd, LIMGRD)
(
 const cs_int_t   *const ndim,        /* --> dimension de l'espace            */
 const cs_int_t   *const ncelet,      /* --> nombre de cellules étendu        */
 const cs_int_t   *const ncel,        /* --> nombre de cellules               */
 const cs_int_t   *const nfac,        /* --> nombre de faces internes         */
 const cs_int_t   *const imrgra,      /* --> mode de calcul du gradient       */
 const cs_int_t   *const imligp,      /* --> mode de limitation du gradient   */
 const cs_int_t   *const idimte,      /* --> 0, 1, 2 : variable
                                             scalaire, vecteur, tenseur       */
 const cs_int_t   *const itenso,      /* --> pour la périodicité de rotation  */
 const cs_int_t   *const iwarnp,      /* --> niveau d'impression              */
 const cs_int_t   *const nfecra,      /* --> unité sortie standard            */
 const cs_real_t  *const climgp,      /* --> coefficient de limitation du grad*/
 const cs_int_t          ifacel[],    /* --> liste des faces internes         */
 const cs_int_t          ipcvse[],    /* --> pour chaque cellule, position
                                             dans ielvse de la première cellule
                                             du voisinage étendu              */
 const cs_int_t          ielvse[],    /* --> numéro des cellules du voisinage
                                             étendu                           */
 const cs_real_t         xyzcen[],    /* --> c.d.g. des cellules              */
 const cs_real_t         pvar[],      /* --> variable dont on calcule le grad */
       cs_real_t         dpdx[],      /* <-> composante x du gradient         */
       cs_real_t         dpdy[],      /* <-> composante y du gradient         */
       cs_real_t         dpdz[],      /* <-> composante z du gradient         */
       cs_real_t         faclim[],    /*  -  tab. de travail local            */
       cs_real_t         denom[],     /*  -  tab. de travail local            */
       cs_real_t         denum[]      /*  -  tab. de travail local            */
);

void CS_PROCF (filtre, FILTRE)
(
 const cs_int_t   *const ncelet,      /* --> nombre de cellules étendu        */
 const cs_int_t   *const ncel,        /* --> nombre de cellules               */
 const cs_int_t   *const nfac,        /* --> nombre de faces internes         */
 const cs_int_t   *const nfabor,      /* --> nombre de faces de bord          */
 const cs_int_t          ifacel[],    /* --> liste des faces internes         */
 const cs_int_t          ifabor[],    /* --> liste des faces de bord          */
 const cs_int_t          ipcvse[],    /* --> pour chaque cellule, position
                                             dans ielvse de la première cellule
                                             du voisinage étendu              */
 const cs_int_t          ielvse[],    /* --> numéro des cellules du voisinage
                                             étendu                           */
 const cs_real_t         volume[],    /* --> volumes des cellules             */
 const cs_real_t         pvar[],      /* --> variable a filtrer               */
       cs_real_t         varfil[],    /* --> resultat de la variable filtree  */
       cs_real_t         w1[],        /*  -  tab. de travail local            */
       cs_real_t         w2[]         /*  -  tab. de travail local            */
);


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

/*----------------------------------------------------------------------------
 * Definition du voisinage étendu pour le calcul des gradients par moindres
 *   carrés : au premier passage, on choisit les cellules et on modifie la
 *   structure maillage qui stocke le voisinage étendu.
 *            ensuite, on se contente de retourner les pointeurs
 *
 * Attention :  le voisinage étendu permet de toucher toutes les cellules
 *              qui partagent un noeud avec une cellule donnée, mais il ne
 *              comporte pas les cellules visibles d'une cellule au travers
 *              des faces (elles sont accessibles par IFACEL)
 *
 * Interface Fortran :
 *
 * SUBROUTINE DEFVSE
 * *****************
 *    & ( NDIM   , NCEL   , NFAC   , NNOD   ,
 *    &   IMRGRA , ANOMAX ,
 *    &   IFACEL , IPNFAC , NODFAC , XYZCEN , SURFAC )
 *
 *----------------------------------------------------------------------------*/

void CS_PROCF (defvse, DEFVSE)
(
 const cs_int_t   *const ndim,    /* --> Dimension = 3                        */
 const cs_int_t   *const ncel,    /* --> Nombre de cellules                   */
 const cs_int_t   *const nfac,    /* --> Nombre de faces internes             */
 const cs_int_t   *const nnod,    /* --> Nombre de sommets                    */
 const cs_int_t   *const imrgra,  /* --> Mode de reconstruction               */
 const cs_real_t  *const anomax,  /* --> Angle de non orthogonalité en radian
                                         au-delà duquel on retient dans le
                                         voisinage étendu de chaque cellule
                                         voisine les cellules dont un sommet
                                         est sur la face                      */
 const cs_int_t   *const ifacel,  /* --> Connectivité face interne -> cellules*/
 const cs_int_t   *const ipnfac,  /* --> Position du premier noeud d'une face
                                         dans nodfac                          */
 const cs_int_t   *const nodfac,  /* --> Numéro des noeuds des faces internes */
 const cs_real_t  *const xyzcen,  /* --> Coord. du centre des cellules        */
 const cs_real_t  *const surfac   /* --> Vecteur surface des faces internes   */
)
{

  cs_maillage_t         *maillage = cs_glob_maillage;

  cs_int_t               ind_som;
  cs_int_t               ind_fac;
  cs_int_t               ind_cel;
  cs_int_t               ind_cel_i;
  cs_int_t               ind_cel_j;
  cs_int_t               ind_bcl_cel;
  cs_int_t               ind_bcl_cel_vse;
  cs_int_t               ind_bcl_fac;
  cs_int_t               num_cel;
  cs_int_t               num_cel_vse;
  cs_int_t               nbr_val_connect_som_cel;
  cs_int_t               nbr_val_connect_som_cel_estime;
  cs_int_t               nbr_cel_elimine;
  cs_int_t               ind_val_last;
  cs_int_t               ind_pos_last;
  cs_int_t               icoo;

  cs_int_t              *ielvse;
  cs_int_t              *ipcvse;
  cs_int_t              *pos_connect_som_fac;
  cs_int_t              *val_connect_som_fac;
  cs_int_t              *pos_connect_som_cel;
  cs_int_t              *val_connect_som_cel;

  cs_real_t              vij[CS_DIM_3];
  cs_real_t              vnf[CS_DIM_3];
  cs_real_t              normeij;
  cs_real_t              normenf;
  cs_real_t              cos_ij_nf;
  cs_real_t              cos_ij_nf_min;

  cs_bool_t              bool_cel_vu;

  const cs_real_t             *xyzcen_avec_voiset;

  static  cs_int_t cs_loc_premier_appel    = 0;

  enum {X, Y, Z};

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





  /* Au premier passage, on sélectionne les cellules */

  if (cs_loc_premier_appel == 0) {

    /* Pour la suite, on note que l'on est passé une fois */

    cs_loc_premier_appel = 1;


    /* Si on n'a pas de voisinage étendu, on s'arrête */
    ipcvse = maillage->pos_connect_voiset_cel_cel_dom;
    ielvse = maillage->val_connect_voiset_cel_cel_dom;

    if (   ipcvse == NULL || ielvse == NULL
        || maillage->ind_type_voiset == CS_MAILLAGE_TYPE_VOISET_SANS) {
      bft_error (__FILE__, __LINE__, 0,
                 _("Lorsque la méthode de calcul des gradients par moindres\n"
                   "carrés sur support étendu est activée (IMRGRA = <%d>), \n"
                   "il est indispensable de générer un voisinage étendu    \n"
                   "en utilisant l'Enveloppe en pré-traitement avec        \n"
                   "l'option  -ve  (COMMANDE_VOISET dans le lanceur)."),
                 (int)(*imrgra));
    }


    /* Selon imrgra, on sélectionne le voisinage de diverses façons */

    /*
         Par défaut : on conserve tout
         =============================
    */

    /*
         imrgra = 3 : sélection sur critère de non orthogonalité
         =======================================================
                    pour chaque face interne, au-dessus d'un angle de non
                    orthogonalité donné, on retient dans le voisinage étendu
                    des deux cellules voisines de la face toutes les cellules
                    qui s'appuient sur un des sommets de la face
    */

    if ((*imrgra) == 3) {

      /*
          On construit d'abord la connectivité sommet -> cellules
          ------------------------------------------------------
          On a déjà face -> sommets, il faut l'inverser
          on utilisera ensuite face -> cellules
      */



      /*
          On construit la connectivité sommet -> faces
          en inversant face -> sommets
      */

      /* Allocation de pos_connect_som_fac */
      BFT_MALLOC(pos_connect_som_fac, (*nnod)+1, cs_int_t);
      /* On compte les faces de chaque sommet, on utilise pos_connect_som_fac*/
      for (ind_som = 0 ; ind_som < (*nnod)+1 ; ind_som++)
        pos_connect_som_fac[ind_som] = 0;
      for (ind_fac = 0 ; ind_fac < (*nfac) ; ind_fac++)
        for (ind_bcl_fac = ipnfac[ind_fac  ] - 1 ;
             ind_bcl_fac < ipnfac[ind_fac+1] - 1 ; ind_bcl_fac++) {
          ind_som = nodfac[ind_bcl_fac] - 1;
          pos_connect_som_fac[ind_som+1]++;
        }
      /* Construction de pos_connect_som_fac */
      pos_connect_som_fac[0] = 1;
      for (ind_som = 0 ; ind_som < (*nnod) ; ind_som++)
        pos_connect_som_fac[ind_som+1] += pos_connect_som_fac[ind_som];
      /* Allocations de val_connect_som_fac */
      BFT_MALLOC(val_connect_som_fac, pos_connect_som_fac[(*nnod)]-1, cs_int_t);
      /* Construction de val_connect_som_fac : pour cela, on utilise
         pos_connect_som_fac pour repérer la case à remplir dans
         val_connect_som_fac : on l'incrémente donc à chaque nouvelle valeur
         ajoutée. A la fin, on a donc décalé les valeurs de pos_connect_som_fac
         d'une case et il faut les remettre en place. */
      for (ind_fac = 0 ; ind_fac < (*nfac) ; ind_fac++)
        for (ind_bcl_fac = ipnfac[ind_fac  ] - 1 ;
             ind_bcl_fac < ipnfac[ind_fac+1] - 1 ; ind_bcl_fac++) {
          ind_som = nodfac[ind_bcl_fac] - 1;
          val_connect_som_fac[pos_connect_som_fac[ind_som]-1] = ind_fac + 1;
          pos_connect_som_fac[ind_som]++;
        }
      for (ind_som = (*nnod) ; ind_som > 0 ; ind_som--)
        pos_connect_som_fac[ind_som] = pos_connect_som_fac[ind_som-1];
      pos_connect_som_fac[0] = 1;

#if 0
      for (ind_som = 0 ; ind_som < (*nnod) ; ind_som++) {
        printf(" sommet %d :: ", ind_som+1) ;
        for (ind_bcl_cel = pos_connect_som_fac[ind_som  ] - 1 ;
             ind_bcl_cel < pos_connect_som_fac[ind_som+1] - 1 ;
             ind_bcl_cel++)
          printf(" %d ", val_connect_som_fac[ind_bcl_cel]) ;
        printf("\n") ;
      }
#endif

      /*
          On construit la connectivité sommet -> cellules
          avec sommet -> faces et face -> cellules
      */

      /* Allocations */
      nbr_val_connect_som_cel_estime = 3*(*nnod);
      BFT_MALLOC(val_connect_som_cel, nbr_val_connect_som_cel_estime, cs_int_t);
      BFT_MALLOC(pos_connect_som_cel, (*nnod)+1, cs_int_t);
      /* Initialisations */
      nbr_val_connect_som_cel = 0;
      pos_connect_som_cel[0] = 1;
      /* Pour tous les sommets */
      for (ind_som = 0 ; ind_som < (*nnod) ; ind_som++) {
        /* On balaye les faces */
        for (ind_bcl_fac = pos_connect_som_fac[ind_som  ] - 1 ;
             ind_bcl_fac < pos_connect_som_fac[ind_som+1] - 1 ; ind_bcl_fac++) {
          ind_fac = val_connect_som_fac[ind_bcl_fac] - 1;
          /* Et on note les cellules (il faut absolument
             qu'elles ne soient pas en plusieurs exemplaires */
          /* Voisin 1 */
          num_cel = ifacel[2*ind_fac  ];
          /* Test si deja vu */
          bool_cel_vu = CS_FALSE;
          ind_bcl_cel = pos_connect_som_cel[ind_som] - 1;
          while (   (bool_cel_vu == CS_FALSE)
                 && (ind_bcl_cel < nbr_val_connect_som_cel)) {
            if (num_cel == val_connect_som_cel[ind_bcl_cel])
              bool_cel_vu = CS_TRUE;
            ind_bcl_cel++;
          }
          /* Si pas encore vu, on le note */
          if (bool_cel_vu == CS_FALSE) {
            if (nbr_val_connect_som_cel + 1 > nbr_val_connect_som_cel_estime) {
              nbr_val_connect_som_cel_estime *= 2;
              BFT_REALLOC(val_connect_som_cel,
                          nbr_val_connect_som_cel_estime, cs_int_t);
            }
            val_connect_som_cel[nbr_val_connect_som_cel] = num_cel;
            nbr_val_connect_som_cel++;
          }
          /* Voisin 2 */
          num_cel = ifacel[2*ind_fac+1];
          /* Test si deja vu */
          bool_cel_vu = CS_FALSE;
          ind_bcl_cel = pos_connect_som_cel[ind_som] - 1;
          while (   (bool_cel_vu == CS_FALSE)
                 && (ind_bcl_cel < nbr_val_connect_som_cel)) {
            if (num_cel == val_connect_som_cel[ind_bcl_cel])
              bool_cel_vu = CS_TRUE;
            ind_bcl_cel++;
          }
          /* Si pas encore vu, on le note */
          if (bool_cel_vu == CS_FALSE) {
            if (nbr_val_connect_som_cel + 1 > nbr_val_connect_som_cel_estime) {
              nbr_val_connect_som_cel_estime *=2;
              BFT_REALLOC(val_connect_som_cel,
                          nbr_val_connect_som_cel_estime, cs_int_t);
            }
            val_connect_som_cel[nbr_val_connect_som_cel] = num_cel;
            nbr_val_connect_som_cel++;
          }
        }
        pos_connect_som_cel[ind_som+1] = nbr_val_connect_som_cel + 1;
      }
      /* Réallocation */
      BFT_REALLOC(val_connect_som_cel, nbr_val_connect_som_cel, cs_int_t);

      /* Libération de sommet -> faces */
      BFT_FREE(val_connect_som_fac);
      BFT_FREE(pos_connect_som_fac);

#if 0
      for (ind_som = 0 ; ind_som < (*nnod) ; ind_som++) {
        printf(" sommet %d :: ", ind_som+1) ;
        for (ind_bcl_cel = pos_connect_som_cel[ind_som  ] - 1 ;
             ind_bcl_cel < pos_connect_som_cel[ind_som+1] - 1 ;
             ind_bcl_cel++)
          printf(" %d ", val_connect_som_cel[ind_bcl_cel]) ;
        printf("\n") ;
      }
#endif

      /*
          On sélectionne enfin les cellules à conserver dans le voisinage
          ---------------------------------------------------------------
          On les marque d'abord, on élimine celles qui sont inutiles ensuite
      */

      /*
          Marquage des cellules à éliminer : en négatif
      */

      assert((*ndim) == CS_DIM_3);

      /* Coordonnées du centre */
      if (   maillage->ind_type_voiset == CS_MAILLAGE_TYPE_VOISET_SEPARE
          && maillage->nbr_cel_fac_par_voiset > 0) {
        xyzcen_avec_voiset = maillage->coord_cel_avec_voiset;
      }
      else
        xyzcen_avec_voiset = xyzcen;


      /* Sélection par les faces internes */
      for (ind_fac = 0 ; ind_fac < (*nfac) ; ind_fac++) {

      /* On calcule le cosinus de l'angle de non orthogonalité
         des faces internes (angle entre la normale et la droite
         passant par le centre des cellules voisines */

        /* Vecteur IJ */
        ind_cel_i = ifacel[2*ind_fac  ] - 1;
        ind_cel_j = ifacel[2*ind_fac+1] - 1;
        icoo = 0;
        vij[icoo]
          = xyzcen_avec_voiset[CS_DIM_3*ind_cel_j+icoo]
          - xyzcen_avec_voiset[CS_DIM_3*ind_cel_i+icoo];
        icoo = 1;
        vij[icoo]
          = xyzcen_avec_voiset[CS_DIM_3*ind_cel_j+icoo]
          - xyzcen_avec_voiset[CS_DIM_3*ind_cel_i+icoo];
        icoo = 2;
        vij[icoo]
          = xyzcen_avec_voiset[CS_DIM_3*ind_cel_j+icoo]
          - xyzcen_avec_voiset[CS_DIM_3*ind_cel_i+icoo];
        normeij = CS_LOC_MODULE(vij);
        assert (normeij > 0.);
        /* Vecteur normal à la face */
        icoo = 0;
        vnf[icoo] = surfac[CS_DIM_3*ind_fac+icoo];
        icoo = 1;
        vnf[icoo] = surfac[CS_DIM_3*ind_fac+icoo];
        icoo = 2;
        vnf[icoo] = surfac[CS_DIM_3*ind_fac+icoo];
        normenf = CS_LOC_MODULE(vnf);
        assert (normenf > 0.);
        /* Produit scalaire vij.vnf */
        cos_ij_nf = (vij[0]*vnf[0]+vij[1]*vnf[1]+vij[2]*vnf[2])
          / (normeij*normenf);

        /* Comparaison à une limite prédéfinie,
           En dessous de cette limite, c'est non orthogonal et
           on garde la cellule dans le voisinage étendu des deux voisines
           de la face (on la marque en négatif, puis on changera tous les
           signes, pour que les cellules à éliminer soient en négatif */

        cos_ij_nf_min = cos((*anomax));

        /* Test sur l'angle de non orthogonalité */
        if (cos_ij_nf < cos_ij_nf_min) {
          /* Pour chaque cellule voisine de la face :
             intersection du groupe des cellules du voisinage étendu
             que l'on n'a pas encore décidé de garder avec
             le groupe de cellules possédant un sommet sur la face */
          /* Cellule voisine 1 (uniquement si elle est réelle) */
          if (ind_cel_i < (*ncel)) {
            ind_cel = ind_cel_i;
            /* Voisinage étendu non encore retenu */
            for (ind_bcl_cel_vse = ipcvse[ind_cel  ] - 1 ;
                 ind_bcl_cel_vse < ipcvse[ind_cel+1] - 1 ;
                 ind_bcl_cel_vse++) {
              num_cel_vse = ielvse[ind_bcl_cel_vse];
              if (num_cel_vse > 0)
                /* Cellules partageant un sommet avec la face */
                for (ind_bcl_fac = ipnfac[ind_fac  ] - 1 ;
                     ind_bcl_fac < ipnfac[ind_fac+1] - 1 ; ind_bcl_fac++) {
                  ind_som = nodfac[ind_bcl_fac] - 1;
                  for (ind_bcl_cel = pos_connect_som_cel[ind_som  ] - 1 ;
                       ind_bcl_cel < pos_connect_som_cel[ind_som+1] - 1 ;
                       ind_bcl_cel++) {
                    num_cel = val_connect_som_cel[ind_bcl_cel];
                    /* Comparaison et sélection */
                    if (num_cel == num_cel_vse && ielvse[ind_bcl_cel_vse] > 0)
                      ielvse[ind_bcl_cel_vse] = - ielvse[ind_bcl_cel_vse];
                  }
                }
            }
          }
          /* Cellule voisine 2 (uniquement si elle est réelle) */
          if (ind_cel_j < (*ncel)) {
            ind_cel = ind_cel_j;
            /* Voisinage étendu non encore retenu */
            for (ind_bcl_cel_vse = ipcvse[ind_cel  ] - 1 ;
                 ind_bcl_cel_vse < ipcvse[ind_cel+1] - 1 ;
                 ind_bcl_cel_vse++) {
              num_cel_vse = ielvse[ind_bcl_cel_vse];
              if (num_cel_vse > 0)
                /* Cellules partageant un sommet avec la face */
                for (ind_bcl_fac = ipnfac[ind_fac  ] - 1 ;
                     ind_bcl_fac < ipnfac[ind_fac+1] - 1 ; ind_bcl_fac++) {
                  ind_som = nodfac[ind_bcl_fac] - 1;
                  for (ind_bcl_cel = pos_connect_som_cel[ind_som  ] - 1 ;
                       ind_bcl_cel < pos_connect_som_cel[ind_som+1] - 1 ;
                       ind_bcl_cel++) {
                    num_cel = val_connect_som_cel[ind_bcl_cel];
                    /* Comparaison et sélection */
                  if (num_cel == num_cel_vse && ielvse[ind_bcl_cel_vse] > 0)
                      ielvse[ind_bcl_cel_vse] = - ielvse[ind_bcl_cel_vse] ;
                  }
                }
            }
          }
        }
      }

      /* Libération de sommet -> cellules */
      BFT_FREE(val_connect_som_cel);
      BFT_FREE(pos_connect_som_cel);

      /* On change tous les signes de ielvse, pour que les cellules à éliminer
         soient en négatif */
      for (ind_bcl_cel = 0 ; ind_bcl_cel < ipcvse[(*ncel)]-1 ; ind_bcl_cel++)
        ielvse[ind_bcl_cel] = - ielvse[ind_bcl_cel];

      /* On élimine les cellules marquées négativement */
      nbr_cel_elimine = 0;
      ind_val_last = -1;
      ind_pos_last =  0;
      for (ind_cel = 0 ; ind_cel < (*ncel) ; ind_cel++) {
        for (ind_bcl_cel = ind_pos_last ;
             ind_bcl_cel < ipcvse[ind_cel+1] - 1 ; ind_bcl_cel++) {
          if (ielvse[ind_bcl_cel] > 0) {
            ind_val_last++;
            ielvse[ind_val_last] = ielvse[ind_bcl_cel];
          }
          else {
            nbr_cel_elimine++;
          }
        }
        ind_pos_last = ipcvse[ind_cel+1] - 1;
        ipcvse[ind_cel+1] = ipcvse[ind_cel+1] - nbr_cel_elimine;
      }

      maillage->nbr_val_connect_voiset_cel_cel_dom
        = maillage->nbr_val_connect_voiset_cel_cel_dom - nbr_cel_elimine;
      /* Réallocation avec la dimension exacte */
      BFT_REALLOC(maillage->val_connect_voiset_cel_cel_dom,
                  ipcvse[(*ncel)]-1, cs_int_t);

      bft_printf(_("\n %d cellules conservées dans le voisinage étendu\n"),
                 maillage->nbr_val_connect_voiset_cel_cel_dom);

#if 0
      for (ind_cel = 0 ; ind_cel < (*ncel) ; ind_cel++) {
        printf(" cellule %d :: ", ind_cel+1);
        for (ind_bcl_cel = ipcvse[ind_cel  ] - 1 ;
             ind_bcl_cel < ipcvse[ind_cel+1] - 1 ;
             ind_bcl_cel++)
          printf(" %d ", ielvse[ind_bcl_cel]);
        printf("\n");
      }
#endif

    }

  }

}


/*----------------------------------------------------------------------------
 * Appel de la routine FORTRAN GRADMC pour le calcul des gradients par moindres
 *   carrés. On ajoute aux arguments le voisinage étendu de la structure
 *   maillage.
 *
 * Interface Fortran :
 *
 * SUBROUTINE CGRDMC
 * *****************
 *
 *    & ( IDBIA0 , IDBRA0 ,
 *    &   NCELET , NCEL   , NFAC   , NFABOR , NCELBR , NITUSE , NRTUSE ,
 *    &   INC    , ICCOCG , NSWRGP , IDIMTE , ITENSO , IPHYDP , IMRGRA ,
 *    &   IWARNP , NFECRA , EPSRGP , EXTRAP ,
 *    &   IFACEL , IFABOR , ICELBR , ITUSER , IA     ,
 *    &   VOLUME , SURFAC , SURFBO , SURFBN , POND   ,
 *    &   DIST   , DISTBR , DIJPF  , DIIPB  ,
 *    &   FEXTX  , FEXTY  , FEXTZ  ,
 *    &   XYZCEN , CDGFAC , CDGFBO , COEFAP , COEFBP , PVAR   ,
 *    &   COCGB  , COCG   , RTUSER ,
 *    &   DPDX   , DPDY   , DPDZ   ,
 *    &   BX     , BY     , BZ     , RA     )
 *
 *----------------------------------------------------------------------------*/

void CS_PROCF (cgrdmc, CGRDMC)
(
 const cs_int_t   *const idbia0,      /* --> numéro 1ère case libre dans IA   */
 const cs_int_t   *const idbra0,      /* --> numéro 1ère case libre dans RA   */
 const cs_int_t   *const ncelet,      /* --> nombre de cellules étendu        */
 const cs_int_t   *const ncel,        /* --> nombre de cellules               */
 const cs_int_t   *const nfac,        /* --> nombre de faces internes         */
 const cs_int_t   *const nfabor,      /* --> nombre de faces de bord          */
 const cs_int_t   *const ncelbr,      /* --> nombre de cellules de bord       */
 const cs_int_t   *const nituse,      /* --> longueur du tableau ituser[]     */
 const cs_int_t   *const nrtuse,      /* --> longueur du tableau rtuser[]     */
 const cs_int_t   *const inc,         /* --> 0 ou 1 : incrément ou non        */
 const cs_int_t   *const iccocg,      /* --> 1 ou 0 : recalcul de COCG ou non */
 const cs_int_t   *const nswrgp,      /* --> >1 : avec reconstruction         */
 const cs_int_t   *const idimte,      /* --> 0, 1, 2 : variable
                                             scalaire, vecteur, tenseur       */
 const cs_int_t   *const itenso,      /* --> pour la périodicité de rotation  */
 const cs_int_t   *const iphydp,      /* --> prise en compte de P hydrostatiqu*/
 const cs_int_t   *const imrgra,      /* --> mode de calcul du gradient       */
 const cs_int_t   *const iwarnp,      /* --> niveau d'impression              */
 const cs_int_t   *const nfecra,      /* --> unité sortie standard            */
 const cs_real_t  *const epsrgp,      /* --> précision pour calcul du gradient
                                             itératif                         */
 const cs_real_t  *const extrap,      /* --> extrapolation des grad. au bord  */
 const cs_int_t          ifacel[],    /* --> liste des faces internes         */
 const cs_int_t          ifabor[],    /* --> liste des faces de bord          */
 const cs_int_t          icelbr[],    /* --> numéro des cellules ayant au
                                             moins une face de bord           */
       cs_int_t          ituser[],    /* <-> tab. complémentaire utilisateur  */
       cs_int_t          ia[],        /* <-> macro-tableau entier             */
 const cs_real_t         volume[],    /* --> volumes des cellules             */
 const cs_real_t         surfac[],    /* --> surfaces des faces internes      */
 const cs_real_t         surfbo[],    /* --> surfaces des faces de bord       */
 const cs_real_t         surfbn[],    /* --> norme de surfbo                  */
 const cs_real_t         pond[],      /* --> pondération géométrique faces int*/
 const cs_real_t         dist[],      /* --> distance entre I' et J' faces int*/
 const cs_real_t         distbr[],    /* --> distance du cdgfbo à I' fac. brd */
 const cs_real_t         dijpf[],     /* --> vecteur I'J' aux faces internes  */
 const cs_real_t         diipb[],     /* --> vecteur II' aux faces de bord    */
 const cs_real_t         fextx[],     /* --> composantes de la force          */
 const cs_real_t         fexty[],     /*     extérieure générant la pression  */
 const cs_real_t         fextz[],     /*     hydrostatique                    */
 const cs_real_t         xyzcen[],    /* --> c.d.g. des cellules              */
 const cs_real_t         cdgfac[],    /* --> c.d.g. des faces internes        */
 const cs_real_t         cdgfbo[],    /* --> c.d.g. des faces de bord         */
 const cs_real_t         coefap[],    /* --> cond. limites aux faces de bord  */
 const cs_real_t         coefbp[],    /* --> cond. limites aux faces de bord  */
 const cs_real_t         pvar[],      /* --> variable dont on calcule le grad */
       cs_real_t         cocgb[],     /* <-> contribution à COCG des faces
                                             internes des cellules de bord    */
       cs_real_t         cocg[],      /* <->  contribution à COCG des faces
                                             internes des cellules de bord    */
       cs_real_t         rtuser[],    /* <-> tab. complémentaire utilisateur  */
       cs_real_t         dpdx[],      /* <-- composante x du gradient         */
       cs_real_t         dpdy[],      /* <-- composante y du gradient         */
       cs_real_t         dpdz[],      /* <-- composante z du gradient         */
       cs_real_t         bx[],        /*  -  tab. de travail local            */
       cs_real_t         by[],        /*  -  tab. de travail local            */
       cs_real_t         bz[],        /*  -  tab. de travail local            */
       cs_real_t         ra[]         /* <-> macro-tableau réel               */
)
{
#if defined(_CS_HAVE_MPI)
        cs_int_t         iel;
        cs_int_t         nceltot;
#endif /*_CS_HAVE_MPI */
        cs_int_t        *ipcvse;
        cs_int_t        *ielvse;
  const cs_real_t       *xyzcen_avec_voiset;
  const cs_real_t       *pvar_avec_voiset;
  const cs_real_t       *fextx_avec_voiset;
  const cs_real_t       *fexty_avec_voiset;
  const cs_real_t       *fextz_avec_voiset;
#if defined(_CS_HAVE_MPI)
        cs_real_t       *pvar_voiset;
        cs_real_t       *fextx_voiset;
        cs_real_t       *fexty_voiset;
        cs_real_t       *fextz_voiset;
        cs_real_t       *fextx_avec_voiset_tmp;
        cs_real_t       *fexty_avec_voiset_tmp;
        cs_real_t       *fextz_avec_voiset_tmp;
#endif /*_CS_HAVE_MPI */

  cs_maillage_t         *maillage = cs_glob_maillage;
#if defined(_CS_HAVE_MPI)
  cs_maillage_tmp_t     *maillage_tmp = cs_glob_maillage_tmp;
#endif /*_CS_HAVE_MPI */


  /*
     Remarque sur les tests faisant intervenir _CS_HAVE_MPI

     Dans cs_maillage.h, la structure cs_maillage_tmp_t n'est definie
       definie que si on _CS_HAVE_MPI (si le parallélisme est possible)
     Ici, on utilise cs_maillage_tmp_t s'il y a un voisinage separe. C'est
       correct car il ne peut y avoir un voisinage separé que en parallele
       etant donne le choix fait dans l'enveloppe (un voisinage separe
       en non parallele n'apporterait rien car il ne permettrait pas de
       reduire le volume des echanges).
     On ne risque donc pas de faire appel à cs_maillage_tmp_t si _CS_HAVE_MPI
       n'est pas defini, mais il faut cependant proteger cs_maillage_tmp_t
       par des tests  sur _CS_HAVE_MPI pour permettre la compilation
       lorsque _CS_HAVE_MPI n'est pas defini
     Par coherence, on fait de meme avec les variables pvar_voiset
        fext?_voiset et fext?_avec_voiset_tmp, nceltot, iel
  */


  /* La complexite du traitement dans cette fonction est liée à la prise en
     compte d'un voisinage étendu séparé du halo classique */


  /* On va transmettre au gradient par moindre carrés une variable
     pvar_avec_voiset qui sera :
       maillage_tmp->var_cel_avec_voiset lorsque le voisinage étendu existe
                                         et qu'il est traité séparément
       pvar                              sinon (cas classiques)

     La variable  maillage_tmp->var_cel_avec_voiset devra contenir
       pvar dans les ncelet premieres cases et les valeurs actualisees
       du voisinage étendu séparé dans le reste, noté pvar_voiset

     Il faut de même transmettre la variable allongée correcte pour les
       centres de gravité des cellules et,
       lorsque la pression hydrostatique est prise en compte, transmettre
       également des variables allongées pour la fext.

     L'utilisation des pointeurs fext?_avec_voiset_tmp au lieu de
       fext?_avec_voiset directement permet d'éviter un warning et de conserver
       le caractere const des fext? en argument.
  */



  /*
     Echanges pour mise à jour du voisinage étendu séparé
     ----------------------------------------------------
  */


#if defined(_CS_HAVE_MPI)
  if (maillage->ind_type_voiset == CS_MAILLAGE_TYPE_VOISET_SEPARE) {

    /* Identification des pointeurs sur le voisinage étendu séparé */

    if (maillage->nbr_cel_fac_par_voiset > 0) {
      pvar_voiset = maillage_tmp->var_cel_avec_voiset + (*ncelet);
    }
    else {
      pvar_voiset = NULL;
    }


    if (*iphydp != 0) {

      if (maillage->nbr_cel_fac_par_voiset > 0) {

        nceltot = (*ncelet)+maillage->nbr_cel_fac_par_voiset;
        BFT_MALLOC(fextx_avec_voiset_tmp,nceltot,cs_real_t);
        BFT_MALLOC(fexty_avec_voiset_tmp,nceltot,cs_real_t);
        BFT_MALLOC(fextz_avec_voiset_tmp,nceltot,cs_real_t);

        fextx_voiset = fextx_avec_voiset_tmp + (*ncelet);
        fexty_voiset = fexty_avec_voiset_tmp + (*ncelet);
        fextz_voiset = fextz_avec_voiset_tmp + (*ncelet);
      }
      else {

        fextx_avec_voiset_tmp = NULL;
        fexty_avec_voiset_tmp = NULL;
        fextz_avec_voiset_tmp = NULL;

        fextx_voiset = NULL;
        fexty_voiset = NULL;
        fextz_voiset = NULL;
      }
    }


    /* Echanges parallèles et périodiques
       Attention, un processeur peut n'avoir rien à recevoir
       (maillage->nbr_cel_fac_par_voiset = 0), mais avoir cependant
       des éléments à envoyer : on doit donc entrer dans parcve
       pour tous les processeurs lorsque l'on est en parallèle   */

    if (maillage->nbr_dom > 1)
      cs_parallel_parcve(pvar, pvar_voiset);

    if (maillage->nbr_per > 0)
      cs_perio_percve   (pvar, pvar_voiset);


    if (*iphydp != 0) {

      if (maillage->nbr_dom > 1) {
        cs_parallel_parcve(fextx, fextx_voiset);
        cs_parallel_parcve(fexty, fexty_voiset);
        cs_parallel_parcve(fextz, fextz_voiset);
      }

      if (maillage->nbr_per > 0) {
        cs_perio_percve   (fextx, fextx_voiset);
        cs_perio_percve   (fexty, fexty_voiset);
        cs_perio_percve   (fextz, fextz_voiset);
      }

    }


  }
#endif /*_CS_HAVE_MPI */


  /*
     Appel du gradient par moindres carres
     -------------------------------------
  */


  /* Identification du pointeur sur la variable allongée
     (de longueur ncelet + voisinage étendu séparé éventuel)
     et copie des valeurs 1-ncelet si nécessaire */

#if defined(_CS_HAVE_MPI)
  if (   maillage->ind_type_voiset == CS_MAILLAGE_TYPE_VOISET_SEPARE
      && maillage->nbr_cel_fac_par_voiset > 0) {

    for (iel = 0 ; iel < (*ncelet) ; iel++)
      maillage_tmp->var_cel_avec_voiset[iel] = pvar[iel];
    pvar_avec_voiset = maillage_tmp->var_cel_avec_voiset;


    if (*iphydp != 0) {

      for (iel = 0 ; iel < (*ncelet) ; iel++)
        fextx_avec_voiset_tmp[iel] = fextx[iel];
      fextx_avec_voiset = fextx_avec_voiset_tmp;
      for (iel = 0 ; iel < (*ncelet) ; iel++)
        fexty_avec_voiset_tmp[iel] = fexty[iel];
      fexty_avec_voiset = fexty_avec_voiset_tmp;
      for (iel = 0 ; iel < (*ncelet) ; iel++)
        fextz_avec_voiset_tmp[iel] = fextz[iel];
      fextz_avec_voiset = fextz_avec_voiset_tmp;

    }
    else {

      fextx_avec_voiset = fextx;
      fexty_avec_voiset = fexty;
      fextz_avec_voiset = fextz;

    }
  }

  else {
#endif /*_CS_HAVE_MPI */

    pvar_avec_voiset = pvar;


    fextx_avec_voiset = fextx;
    fexty_avec_voiset = fexty;
    fextz_avec_voiset = fextz;

#if defined(_CS_HAVE_MPI)
  }
#endif /*_CS_HAVE_MPI */

  /* Identification du pointeur sur les coordonnées allongées
     (de longueur ncelet + voisinage étendu séparé éventuel)*/

  if (   maillage->ind_type_voiset == CS_MAILLAGE_TYPE_VOISET_SEPARE
      && maillage->nbr_cel_fac_par_voiset > 0) {
    xyzcen_avec_voiset = maillage->coord_cel_avec_voiset;
  }
  else
    xyzcen_avec_voiset = xyzcen;


  /*  Pointeurs sur le voisinage étendu */

  ipcvse = maillage->pos_connect_voiset_cel_cel_dom;
  ielvse = maillage->val_connect_voiset_cel_cel_dom;

  /*  Appel du calcul de gradient */

  CS_PROCF(gradmc, GRADMC)
    (idbia0, idbra0,
     ncelet, ncel  , nfac  , nfabor, ncelbr, nituse, nrtuse,
     inc   , iccocg, nswrgp, idimte, itenso, iphydp, imrgra,
     iwarnp, nfecra, epsrgp, extrap,
     ifacel, ifabor, icelbr, ipcvse, ielvse, ituser, ia    ,
     volume, surfac, surfbo, surfbn, pond  ,
     dist  , distbr, dijpf , diipb ,
     fextx_avec_voiset, fexty_avec_voiset, fextz_avec_voiset,
     xyzcen_avec_voiset,
     cdgfac, cdgfbo, coefap, coefbp, pvar_avec_voiset,
     cocgb , cocg  , rtuser,
     dpdx  , dpdy  , dpdz  ,
     bx    , by    , bz    , ra   );


  /* Libération si on a réservé quelque chose */

#if defined(_CS_HAVE_MPI)
  if (   (maillage->ind_type_voiset == CS_MAILLAGE_TYPE_VOISET_SEPARE)
      && (*iphydp != 0)
      && (maillage->nbr_cel_fac_par_voiset > 0)) {
    BFT_FREE(fextx_avec_voiset_tmp);
    BFT_FREE(fexty_avec_voiset_tmp);
    BFT_FREE(fextz_avec_voiset_tmp);
  }
#endif /*_CS_HAVE_MPI */
}


/*----------------------------------------------------------------------------
 * Appel de la routine FORTRAN LIMGRD pour la limitation éventuelle des
 *   gradients. On ajoute aux arguments le voisinage étendu de la structure
 *   maillage.
 *
 * Interface Fortran :
 *
 * SUBROUTINE CLMGRD
 * *****************
 *
 *    & ( NDIM   , NCELET , NCEL   , NFAC   ,
 *    &   IMRGRA , IMLIGP , IDIMTE , ITENSO , IWARNP , NFECRA , CLIMGP ,
 *    &   IFACEL ,
 *    &   XYZCEN , PVAR   ,
 *    &   DPDX   , DPDY   , DPDZ   ,
 *    &   FACLIM , DENOM  , DENUM  )
 *
 *----------------------------------------------------------------------------*/

void CS_PROCF (clmgrd, CLMGRD)
(
 const cs_int_t   *const ndim,        /* --> dimension de l'espace            */
 const cs_int_t   *const ncelet,      /* --> nombre de cellules étendu        */
 const cs_int_t   *const ncel,        /* --> nombre de cellules               */
 const cs_int_t   *const nfac,        /* --> nombre de faces internes         */
 const cs_int_t   *const imrgra,      /* --> mode de calcul du gradient       */
 const cs_int_t   *const imligp,      /* --> mode de limitation du gradient   */
 const cs_int_t   *const idimte,      /* --> 0, 1, 2 : variable
                                             scalaire, vecteur, tenseur       */
 const cs_int_t   *const itenso,      /* --> pour la périodicité de rotation  */
 const cs_int_t   *const iwarnp,      /* --> niveau d'impression              */
 const cs_int_t   *const nfecra,      /* --> unité sortie standard            */
 const cs_real_t  *const climgp,      /* --> coefficient de limitation du grad*/
 const cs_int_t          ifacel[],    /* --> liste des faces internes         */
 const cs_real_t         xyzcen[],    /* --> c.d.g. des cellules              */
 const cs_real_t         pvar[],      /* --> variable dont on calcule le grad */
       cs_real_t         dpdx[],      /* <-> composante x du gradient         */
       cs_real_t         dpdy[],      /* <-> composante y du gradient         */
       cs_real_t         dpdz[],      /* <-> composante z du gradient         */
       cs_real_t         faclim[],    /*  -  tab. de travail local            */
       cs_real_t         denom[],     /*  -  tab. de travail local            */
       cs_real_t         denum[]      /*  -  tab. de travail local            */
)
{
#if defined(_CS_HAVE_MPI)
        cs_int_t         iel;
        cs_int_t         nceltot;
#endif /*_CS_HAVE_MPI */
        cs_int_t        *ipcvse;
        cs_int_t        *ielvse;
  const cs_real_t       *xyzcen_avec_voiset;
  const cs_real_t       *pvar_avec_voiset;
        cs_real_t       *dpdx_avec_voiset;
        cs_real_t       *dpdy_avec_voiset;
        cs_real_t       *dpdz_avec_voiset;
        cs_real_t       *denum_avec_voiset;
        cs_real_t       *denom_avec_voiset;
#if defined(_CS_HAVE_MPI)
        cs_real_t       *pvar_voiset;
        cs_real_t       *dpdx_voiset;
        cs_real_t       *dpdy_voiset;
        cs_real_t       *dpdz_voiset;
        cs_real_t       *dpdx_avec_voiset_tmp;
        cs_real_t       *dpdy_avec_voiset_tmp;
        cs_real_t       *dpdz_avec_voiset_tmp;
        cs_real_t       *denum_avec_voiset_tmp;
        cs_real_t       *denom_avec_voiset_tmp;
#endif /*_CS_HAVE_MPI */

  cs_maillage_t         *maillage = cs_glob_maillage;
#if defined(_CS_HAVE_MPI)
  cs_maillage_tmp_t     *maillage_tmp = cs_glob_maillage_tmp;
#endif /*_CS_HAVE_MPI */


  /*
     Remarque sur les tests faisant intervenir _CS_HAVE_MPI

     Dans cs_maillage.h, la structure cs_maillage_tmp_t n'est definie
       definie que si on _CS_HAVE_MPI (si le parallélisme est possible)
     Ici, on utilise cs_maillage_tmp_t s'il y a un voisinage separe. C'est
       correct car il ne peut y avoir un voisinage separé que en parallele
       etant donne le choix fait dans l'enveloppe (un voisinage separe
       en non parallele n'apporterait rien car il ne permettrait pas de
       reduire le volume des echanges).
     On ne risque donc pas de faire appel à cs_maillage_tmp_t si _CS_HAVE_MPI
       n'est pas defini, mais il faut cependant proteger cs_maillage_tmp_t
       par des tests  sur _CS_HAVE_MPI pour permettre la compilation
       lorsque _CS_HAVE_MPI n'est pas defini
     Par coherence, on fait de meme avec les variables pvar_voiset
        fext?_voiset et fext?_avec_voiset_tmp, nceltot, iel
  */


  /* La complexite du traitement dans cette fonction est liée à la prise en
     compte d'un voisinage étendu séparé du halo classique */


  /* On va transmettre au gradient par moindre carrés une variable
     pvar_avec_voiset qui sera :
       maillage_tmp->var_cel_avec_voiset lorsque le voisinage étendu existe
                                         et qu'il est traité séparément
       pvar                              sinon (cas classiques)

     La variable  maillage_tmp->var_cel_avec_voiset devra contenir
       pvar dans les ncelet premieres cases et les valeurs actualisees
       du voisinage étendu séparé dans le reste, noté pvar_voiset

     Il faut de même transmettre la variable allongée correcte pour les
       centres de gravité des cellules et,
       pour le mode de limitation imligp = 1, transmettre également
       une variable allongée pour le gradient et les tableaux de travail
       denom et denum.
  */



  /*
     Echanges pour mise à jour du voisinage étendu séparé
     ----------------------------------------------------
  */


#if defined(_CS_HAVE_MPI)
  if (maillage->ind_type_voiset == CS_MAILLAGE_TYPE_VOISET_SEPARE) {

    /* Identification des pointeurs sur le voisinage étendu séparé */

    if (maillage->nbr_cel_fac_par_voiset > 0) {
      pvar_voiset =   maillage_tmp->var_cel_avec_voiset
	            + maillage->nbr_cel_etendu;
    }
    else {
      pvar_voiset = NULL;
    }


    /* Pour le gradient et deux tableaux de travail denom et denum.
       L'utilisation d'un pointeur _tmp pour le gradient (au lieu d'un
       pointeur dpdx_avec_voiset directement) provient du fait qu'on a recopié
       la fonction CGRDMC. Ici ce n'était pas indispensable, dpdx en entrée
       n'étant pas const.
       Pour denom et denum, pas nécessaire d'initialiser les tableaux
       on n'a donc pas besoin de denom_voiset et denum_voiset */

    if (*imligp == 1) {

      if (maillage->nbr_cel_fac_par_voiset > 0) {

        nceltot = maillage->nbr_cel_etendu + maillage->nbr_cel_fac_par_voiset;

        BFT_MALLOC(dpdx_avec_voiset_tmp, nceltot, cs_real_t);
        BFT_MALLOC(dpdy_avec_voiset_tmp, nceltot, cs_real_t);
        BFT_MALLOC(dpdz_avec_voiset_tmp, nceltot, cs_real_t);

        BFT_MALLOC(denum_avec_voiset_tmp, nceltot, cs_real_t);
        BFT_MALLOC(denom_avec_voiset_tmp, nceltot, cs_real_t);

        dpdx_voiset = dpdx_avec_voiset_tmp + maillage->nbr_cel_etendu;
        dpdy_voiset = dpdy_avec_voiset_tmp + maillage->nbr_cel_etendu;
        dpdz_voiset = dpdz_avec_voiset_tmp + maillage->nbr_cel_etendu;

      }
      else {

        dpdx_avec_voiset_tmp = NULL;
        dpdy_avec_voiset_tmp = NULL;
        dpdz_avec_voiset_tmp = NULL;

        denum_avec_voiset_tmp = NULL;
        denom_avec_voiset_tmp = NULL;

        dpdx_voiset = NULL;
        dpdy_voiset = NULL;
        dpdz_voiset = NULL;

      }
    }


    /* Echanges parallèles et périodiques
       Attention, un processeur peut n'avoir rien à recevoir
       (maillage->nbr_cel_fac_par_voiset = 0), mais avoir cependant
       des éléments à envoyer : on doit donc entrer dans parcve
       pour tous les processeurs lorsque l'on est en parallèle   */

    if (maillage->nbr_dom > 1)
      cs_parallel_parcve(pvar, pvar_voiset);

    if (maillage->nbr_per > 0)
      cs_perio_percve(pvar, pvar_voiset);

    /* Echange pour les gradients
       (Pour les tableaux de travail, échange inutile) */
    if (*imligp == 1) {

      if (maillage->nbr_dom > 1) {
        cs_parallel_parcve(dpdx, dpdx_voiset);
        cs_parallel_parcve(dpdy, dpdy_voiset);
        cs_parallel_parcve(dpdz, dpdz_voiset);
      }

      if (maillage->nbr_per > 0) {
        cs_perio_percve(dpdx, dpdx_voiset);
        cs_perio_percve(dpdy, dpdy_voiset);
        cs_perio_percve(dpdz, dpdz_voiset);
      }

    }


  }
#endif /*_CS_HAVE_MPI */

  /*
     Appel de la limitation du gradient
     ----------------------------------
  */


  /* Identification du pointeur sur la variable allongée
     (de longueur ncelet + voisinage étendu séparé éventuel)
     et copie des valeurs 1-ncelet si nécessaire */

#if defined(_CS_HAVE_MPI)
  if (   maillage->ind_type_voiset == CS_MAILLAGE_TYPE_VOISET_SEPARE
      && maillage->nbr_cel_fac_par_voiset > 0) {

    for (iel = 0 ; iel < (*ncelet) ; iel++)
       maillage_tmp->var_cel_avec_voiset[iel] = pvar[iel];
    pvar_avec_voiset = maillage_tmp->var_cel_avec_voiset;


    if (*imligp == 1) {

      for (iel = 0 ; iel < (*ncelet) ; iel++)
        dpdx_avec_voiset_tmp[iel] = dpdx[iel];
      dpdx_avec_voiset = dpdx_avec_voiset_tmp;
      for (iel = 0 ; iel < (*ncelet) ; iel++)
        dpdy_avec_voiset_tmp[iel] = dpdy[iel];
      dpdy_avec_voiset = dpdy_avec_voiset_tmp;
      for (iel = 0 ; iel < (*ncelet) ; iel++)
        dpdz_avec_voiset_tmp[iel] = dpdz[iel];
      dpdz_avec_voiset = dpdz_avec_voiset_tmp;

      denum_avec_voiset = denum_avec_voiset_tmp;
      denom_avec_voiset = denom_avec_voiset_tmp;

    }
    else {

      dpdx_avec_voiset = dpdx;
      dpdy_avec_voiset = dpdy;
      dpdz_avec_voiset = dpdz;

      denum_avec_voiset = denum;
      denom_avec_voiset = denom;

    }
  }
  else {
#endif /*_CS_HAVE_MPI */

    pvar_avec_voiset = pvar;

    dpdx_avec_voiset = dpdx;
    dpdy_avec_voiset = dpdy;
    dpdz_avec_voiset = dpdz;

    denum_avec_voiset = denum;
    denom_avec_voiset = denom;
#if defined(_CS_HAVE_MPI)
  }
#endif /*_CS_HAVE_MPI */


  /* Identification du pointeur sur les coordonnées allongées
     (de longueur ncelet + voisinage étendu séparé éventuel)*/

  if (   maillage->ind_type_voiset == CS_MAILLAGE_TYPE_VOISET_SEPARE
      && maillage->nbr_cel_fac_par_voiset > 0) {
    xyzcen_avec_voiset = maillage->coord_cel_avec_voiset;
  }
  else
    xyzcen_avec_voiset = xyzcen;


  /*  Pointeurs sur le voisinage étendu */

  ipcvse = maillage->pos_connect_voiset_cel_cel_dom;
  ielvse = maillage->val_connect_voiset_cel_cel_dom;

  /*  Appel de la limitation du gradient */

  CS_PROCF(limgrd, LIMGRD)
    (ndim, ncelet, ncel, nfac,
     imrgra, imligp, idimte, itenso, iwarnp, nfecra, climgp,
     ifacel, ipcvse, ielvse,
     xyzcen_avec_voiset, pvar_avec_voiset,
     dpdx_avec_voiset, dpdy_avec_voiset, dpdz_avec_voiset,
     faclim, denom_avec_voiset, denum_avec_voiset);

  /*  Transfert du gradient limité si on a utilisé un tableau différent
      puis libération des tableaux locaux
      (voisinage étendu séparé, non vide et avec imligp == 1) */

#if defined(_CS_HAVE_MPI)
  if (   (maillage->ind_type_voiset == CS_MAILLAGE_TYPE_VOISET_SEPARE)
      && (maillage->nbr_cel_fac_par_voiset > 0)
      && (*imligp == 1)) {

    for (iel = 0 ; iel < (*ncelet) ; iel++)
      dpdx[iel] = dpdx_avec_voiset[iel];
    for (iel = 0 ; iel < (*ncelet) ; iel++)
      dpdy[iel] = dpdy_avec_voiset[iel];
    for (iel = 0 ; iel < (*ncelet) ; iel++)
      dpdz[iel] = dpdz_avec_voiset[iel];

    BFT_FREE(dpdx_avec_voiset_tmp);
    BFT_FREE(dpdy_avec_voiset_tmp);
    BFT_FREE(dpdz_avec_voiset_tmp);

    BFT_FREE(denum_avec_voiset_tmp);
    BFT_FREE(denom_avec_voiset_tmp);
  }
#endif /*_CS_HAVE_MPI */

}

/*----------------------------------------------------------------------------
 * Appel de la routine FORTRAN CFILTR pour le calcul des filtres avec voisinages
 * etendus pour le modele dynamique.
 *
 * Interface Fortran :
 *
 * SUBROUTINE CFILTR
 * *****************
 *
 *    & ( NCELET , NCEL   , NFAC   , NFABOR ,
 *    &   IFACEL , IFABOR ,
 *    &   VOLUME ,
 *    &   PVAR   , VARFIL ,
 *    &   W1     , W2     )
 *
 *----------------------------------------------------------------------------*/

void CS_PROCF (cfiltr, CFILTR)
(
 const cs_int_t   *const ncelet,      /* --> nombre de cellules étendu        */
 const cs_int_t   *const ncel,        /* --> nombre de cellules               */
 const cs_int_t   *const nfac,        /* --> nombre de faces internes         */
 const cs_int_t   *const nfabor,      /* --> nombre de faces de bord          */
 const cs_int_t          ifacel[],    /* --> liste des faces internes         */
 const cs_int_t          ifabor[],    /* --> liste des faces de bord          */
 const cs_real_t         volume[],    /* --> volumes des cellules             */
 const cs_real_t         pvar[],      /* --> variable a filter                */
       cs_real_t         varfil[],    /* --> resultat de la variable filtree  */
       cs_real_t         w1[],        /*  -  tab. de travail local            */
       cs_real_t         w2[]         /*  -  tab. de travail local            */
)
{
#if defined(_CS_HAVE_MPI)
        cs_int_t         iel;
        cs_int_t         nceltot;
#endif /*_CS_HAVE_MPI */
        cs_int_t        *ipcvse;
        cs_int_t        *ielvse;
  const cs_real_t       *pvar_avec_voiset;
  const cs_real_t       *volume_avec_voiset;
#if defined(_CS_HAVE_MPI)
        cs_real_t       *pvar_voiset;
        cs_real_t       *volume_voiset;
        cs_real_t       *volume_avec_voiset_tmp;
#endif /*_CS_HAVE_MPI */

  cs_maillage_t         *maillage = cs_glob_maillage;
#if defined(_CS_HAVE_MPI)
  cs_maillage_tmp_t     *maillage_tmp = cs_glob_maillage_tmp;
#endif /*_CS_HAVE_MPI */


  /*
     Remarque sur les tests faisant intervenir _CS_HAVE_MPI

     Dans cs_maillage.h, la structure cs_maillage_tmp_t n'est definie
       definie que si on _CS_HAVE_MPI (si le parallélisme est possible)
     Ici, on utilise cs_maillage_tmp_t s'il y a un voisinage separe. C'est
       correct car il ne peut y avoir un voisinage separé que en parallele
       etant donne le choix fait dans l'enveloppe (un voisinage separe
       en non parallele n'apporterait rien car il ne permettrait pas de
       reduire le volume des echanges).
     On ne risque donc pas de faire appel à cs_maillage_tmp_t si _CS_HAVE_MPI
       n'est pas defini, mais il faut cependant proteger cs_maillage_tmp_t
       par des tests  sur _CS_HAVE_MPI pour permettre la compilation
       lorsque _CS_HAVE_MPI n'est pas defini
     Par coherence, on fait de meme avec les variables pvar_voiset
        nceltot, iel
  */


  /* La complexite du traitement dans cette fonction est liée à la prise en
     compte d'un voisinage étendu séparé du halo classique */


  /* On va transmettre au gradient par moindre carrés une variable
     pvar_avec_voiset qui sera :
       maillage_tmp->var_cel_avec_voiset lorsque le voisinage étendu existe
                                         et qu'il est traité séparément
       pvar                              sinon (cas classiques)

     La variable  maillage_tmp->var_cel_avec_voiset devra contenir
       pvar dans les ncelet premieres cases et les valeurs actualisees
       du voisinage étendu séparé dans le reste, noté pvar_voiset

     Il faut de même transmettre la variable allongée correcte pour les
       centres de gravité des cellules et,
       lorsque la pression hydrostatique est prise en compte, transmettre
       également des variables allongées pour la fext.

     L'utilisation des pointeurs fext?_avec_voiset_tmp au lieu de
       fext?_avec_voiset directement permet d'éviter un warning et de conserver
       le caractere const des fext? en argument.
  */





  /*
     Echanges pour mise à jour du voisinage étendu séparé
     ----------------------------------------------------
  */


#if defined(_CS_HAVE_MPI)
  if (maillage->ind_type_voiset == CS_MAILLAGE_TYPE_VOISET_SEPARE) {

    /* Identification des pointeurs sur le voisinage étendu séparé */

    if (maillage->nbr_cel_fac_par_voiset > 0) {
      pvar_voiset = maillage_tmp->var_cel_avec_voiset + (*ncelet);
    }
    else {
      pvar_voiset = NULL;
    }

    if (maillage->nbr_cel_fac_par_voiset > 0) {
      nceltot = (*ncelet)+maillage->nbr_cel_fac_par_voiset;
      BFT_MALLOC(volume_avec_voiset_tmp,nceltot,cs_real_t);
      volume_voiset = volume_avec_voiset_tmp + (*ncelet);
    }
    else {
      volume_avec_voiset_tmp = NULL;
      volume_voiset = NULL;
    }

    /* Echanges parallèles et périodiques
       Attention, un processeur peut n'avoir rien à recevoir
       (maillage->nbr_cel_fac_par_voiset = 0), mais avoir cependant
       des éléments à envoyer : on doit donc entrer dans parcve
       pour tous les processeurs lorsque l'on est en parallèle   */

    if (maillage->nbr_dom > 1)
      cs_parallel_parcve(pvar, pvar_voiset);
      cs_parallel_parcve(volume, volume_voiset);

    if (maillage->nbr_per > 0)
      cs_perio_percve   (pvar, pvar_voiset);
      cs_perio_percve   (volume, volume_voiset);
  }

#endif /*_CS_HAVE_MPI */

#if defined(_CS_HAVE_MPI)
  if (   maillage->ind_type_voiset == CS_MAILLAGE_TYPE_VOISET_SEPARE
      && maillage->nbr_cel_fac_par_voiset > 0) {
    for (iel = 0 ; iel < (*ncelet) ; iel++)
      maillage_tmp->var_cel_avec_voiset[iel] = pvar[iel];
    pvar_avec_voiset = maillage_tmp->var_cel_avec_voiset;
    for (iel = 0 ; iel < (*ncelet) ; iel++)
      volume_avec_voiset_tmp[iel] = volume[iel];
    volume_avec_voiset = volume_avec_voiset_tmp;
  }

  else{
#endif /*_CS_HAVE_MPI */

    pvar_avec_voiset = pvar;
    volume_avec_voiset = volume;
#if defined(_CS_HAVE_MPI)
  }
#endif /*_CS_HAVE_MPI */

  /*  Pointeurs sur le voisinage étendu */

  ipcvse = maillage->pos_connect_voiset_cel_cel_dom;
  ielvse = maillage->val_connect_voiset_cel_cel_dom;


  /* Identification du pointeur sur la variable allongée
     (de longueur ncelet + voisinage étendu séparé éventuel)
     et copie des valeurs 1-ncelet si nécessaire */

  CS_PROCF(filtre, FILTRE)
    (ncelet, ncel  , nfac  , nfabor,
     ifacel, ifabor, ipcvse, ielvse,
     volume_avec_voiset,
     pvar_avec_voiset, varfil,
     w1  , w2 );


  /* Libération si on a réservé quelque chose */

#if defined(_CS_HAVE_MPI)
  if (   (maillage->ind_type_voiset == CS_MAILLAGE_TYPE_VOISET_SEPARE)
      && (maillage->nbr_cel_fac_par_voiset > 0)) {
    BFT_FREE(volume_avec_voiset_tmp);
  }
#endif /*_CS_HAVE_MPI */
}

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

#ifdef __cplusplus
}
#endif /* __cplusplus */

