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

/*============================================================================
 * Fonctions associées à la méthode ALE
 *============================================================================*/

/* includes système */

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


/* Includes BFT */

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

/* Includes FVM */

#include <fvm_interface.h>
#include <fvm_nodal.h>
#include <fvm_nodal_extract.h>
#include <fvm_writer.h>

/* Includes librairie */

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


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


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

static fvm_interface_set_t  *_interface_ale = NULL;


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



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

/*----------------------------------------------------------------------------
 * Calcul des centres de gravité des faces et cellules, et des volumes
 *
 * Interface Fortran :
 *
 * SUBROUTINE ALGRMA
 * *****************
 *
 *----------------------------------------------------------------------------*/

void CS_PROCF (algrma, ALGRMA)
(
)
{
  /* Calcul des centres de gravité des faces et cellules, et des volumes */

  cs_maillage_grd_calc(cs_glob_maillage ,
                       cs_glob_maillage_grd);
}



/*----------------------------------------------------------------------------
 * Projection aux noeuds du maillage du deplacement calcule au centre des
 *  cellules
 *
 * Interface Fortran :
 *
 * SUBROUTINE ALDEPL
 * *****************
 *
 * INTEGER          IFACEL(2,NFAC)  : --> : Connectivite faces internes/cellules
 * INTEGER          IFABOR(NFABOR)  : --> : Connectivite faces de bord/cellules
 * INTEGER          IPNFAC(NFAC+1)  : --> : Position du premier noeud de chaque face interne
 * INTEGER          NODFAC(LNDFAC)  : --> : Connectivite faces internes/noeuds
 * INTEGER          IPNFBR(NFABOR+1): --> : Position du premier noeud de chaque face de bord
 * INTEGER          NODFBR(LNDFBR)  : --> : Connectivite faces de bord/noeuds
 * DOUBLE PRECISION UMA(NCELET)    : --> : Vitesse de maillage selon X
 * DOUBLE PRECISION VMA(NCELET)    : --> : Vitesse de maillage selon X
 * DOUBLE PRECISION WMA(NCELET)    : --> : Vitesse de maillage selon X
 * DOUBLE PRECISION COEFAU(NCELET) : --> : Condition aux limites A pour UMA
 * DOUBLE PRECISION COEFAV(NCELET) : --> : Condition aux limites A pour VMA
 * DOUBLE PRECISION COEFAW(NCELET) : --> : Condition aux limites A pour WMA
 * DOUBLE PRECISION COEFBU(NCELET) : --> : Condition aux limites B pour UMA
 * DOUBLE PRECISION COEFBV(NCELET) : --> : Condition aux limites B pour VMA
 * DOUBLE PRECISION COEFBW(NCELET) : --> : Condition aux limites B pour WMA
 * DOUBLE PRECISION DT(NCELET)     : --> : Pas de temps
 * DOUBLE PRECISION DEPROJ(NNOD,3)): <-- : Deplacement projete aux noeuds

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


void CS_PROCF (aldepl, ALDEPL)
(
 const cs_int_t ifacel[],
 const cs_int_t ifabor[],
 const cs_int_t ipnfac[],
 const cs_int_t nodfac[],
 const cs_int_t ipnfbr[],
 const cs_int_t nodfbr[],
 cs_real_t *uma,
 cs_real_t *vma,
 cs_real_t *wma,
 cs_real_t *coefau,
 cs_real_t *coefav,
 cs_real_t *coefaw,
 cs_real_t *coefbu,
 cs_real_t *coefbv,
 cs_real_t *coefbw,
 cs_real_t *dt,
 cs_real_t *deproj
)
{
  cs_int_t  ifac;
  cs_int_t  idim;
  cs_int_t  inod, inodfa, inodfb;
  cs_int_t  nbr_fbr, nbr_fac, nbr_cel;
  cs_int_t  ndim;
  cs_int_t  nbr_som;
  cs_int_t  iel, iel1,iel2;
  cs_real_t *xnbfacnod = NULL;

  nbr_som = cs_glob_maillage->nbr_som;
  nbr_cel = cs_glob_maillage->nbr_cel;
  nbr_fbr = cs_glob_maillage->nbr_fbr;
  nbr_fac = cs_glob_maillage->nbr_fac;
  ndim    = cs_glob_maillage->dim;

  if (cs_glob_maillage->num_som != NULL && _interface_ale == NULL) {
    _interface_ale
      = fvm_interface_set_create(nbr_som, NULL, cs_glob_maillage->num_som);
  }

  BFT_MALLOC(xnbfacnod, nbr_som, cs_real_t);

  for (inod = 0 ; inod < nbr_som ; inod++)
  {
    xnbfacnod[inod] = 0.;
    for (idim = 0 ; idim < ndim ; idim++)
      deproj[nbr_som*idim + inod] = 0.;
  }

  /* Traitement des faces internes d'abord */

  for (ifac = 0 ; ifac < nbr_fac ; ifac++) {
    iel1 = ifacel[ifac*2]-1;
    iel2 = ifacel[ifac*2+1]-1;
    if(iel1 <= nbr_cel) /* test pour ne prendre en compte les faces qu'une fois*/
    for (inodfa = ipnfac[ifac]; inodfa < ipnfac[ifac+1]; inodfa++) {
      inod = nodfac[inodfa-1] -1; /* numero de noeuds pour les faces */
      deproj[inod] = deproj[inod] + 0.5 *(dt[iel1]*uma[iel1]+dt[iel2]*uma[iel2]);
      deproj[inod+nbr_som]
        = deproj[inod+nbr_som] + 0.5 *(dt[iel1]*vma[iel1]+dt[iel2]*vma[iel2]);
      deproj[inod+2*nbr_som]
        = deproj[inod+2*nbr_som] + 0.5 *(dt[iel1]*wma[iel1]+dt[iel2]*wma[iel2]);
      xnbfacnod[inod]=xnbfacnod[inod]+1.;
      }
  }

  /* Traitement des faces de bord.
     On remet a zero pour les noeuds des faces de bord, pour n'y garder que des
     contributions des faces de bord */

  for (ifac = 0 ; ifac < nbr_fbr ; ifac++) {
    for (inodfb = ipnfbr[ifac]; inodfb < ipnfbr[ifac+1]; inodfb++) {
      inod = nodfbr[inodfb-1]-1;
      xnbfacnod[inod] = 0.;
      for (idim = 0 ; idim < ndim ; idim++) deproj[inod+idim*nbr_som]=0.;
    }
  }

  for (ifac = 0 ; ifac < nbr_fbr ; ifac++) {
    iel = ifabor[ifac] - 1;
    for (inodfb = ipnfbr[ifac]; inodfb < ipnfbr[ifac+1]; inodfb++) {
      inod = nodfbr[inodfb-1]-1;
      deproj[inod] =   deproj[inod]
                     + dt[iel]*(coefau[ifac] + coefbu[ifac]*uma[iel]);
      deproj[inod+nbr_som] =   deproj[inod+nbr_som]
                             + dt[iel]*(coefav[ifac] + coefbv[ifac]*vma[iel]);
      deproj[inod+2*nbr_som] =   deproj[inod+2*nbr_som]
                               + dt[iel]*(coefaw[ifac] + coefbw[ifac]*wma[iel]);
      xnbfacnod[inod]=xnbfacnod[inod]+1.;
    }
  }

  if (cs_glob_maillage->num_som != NULL) {
    cs_parallel_interface_sr(_interface_ale, nbr_som, 3, deproj);
    cs_parallel_interface_sr(_interface_ale, nbr_som, 1, xnbfacnod);
  }


  for (inod = 0 ; inod < nbr_som ; inod++) {
      for (idim = 0 ; idim < ndim ; idim++)
          deproj[inod + idim*nbr_som] = deproj[inod + idim*nbr_som]/xnbfacnod[inod];
  }

  BFT_FREE(xnbfacnod);
}
