/*============================================================================
 *
 *                    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 a la périodicité pour les particules
 *============================================================================*/

/* includes système */

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

/* Includes BFT */

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

/* Includes librairie */

#include "cs_base.h"
#include "cs_msg.h"
#include "cs_maillage.h"
#include "cs_perio.h"
#include "cs_lagr_perio.h"

#ifdef __cplusplus
extern "C" {
#if 0
} /* Fausse accolade pour ramener l'auto-indentation d'Emacs à la colonne 0 */
#endif
#endif /* __cplusplus */

/*============================================================================
 * Fonctions publiques : API FORTRAN
 *============================================================================*/

/*----------------------------------------------------------------------------
 * Construction des tableaux permettant d'associer à une cellule halo donnée :
 *  - la cellule réelle correspondante,
 *  - un numéro et sens de périodicité.
 *
 * Interface Fortran :
 *
 * SUBROUTINE PERLOC
 * *****************
 *
 * INTEGER NCELET              :  -> : nombre d'elements halo compris
 * INTEGER NCEL                :  -> : nombre d'elements actifs
 * INTEGER ICELCR(NCELET-NCEL) : <-  : numero de cellule correspondante
 * INTEGER IPERCR(NCELET-NCEL) : <-  : numero de periodicite associe,
 *                                     * 2 + (0 : direct, 1 : indirect)
 *----------------------------------------------------------------------------*/

void CS_PROCF (perloc, PERLOC)(const cs_int_t  *ncelet,
                               const cs_int_t  *ncel,
                                     cs_int_t  *icelcr,
                                     cs_int_t  *ipercr)
{
  cs_int_t ind_loc;
  cs_int_t dbloc;

  cs_int_t iper;
  cs_int_t ielt;

  cs_maillage_t  *maillage = cs_glob_maillage;


  for (ind_loc = 0; ind_loc < maillage->nbr_dom_fac_par; ind_loc++) {

    dbloc = (2*maillage->nbr_per + 1)*ind_loc;

    /* On traite les cellules halos periodiques 'simples' du domaine */

    if (maillage->nbr_dom == 1 ||
        maillage->num_dom_fac_par[ind_loc] == maillage->num_dom) {

      for (iper = 0; iper < maillage->nbr_per; iper++) {

        /* Sens direct */

        for (ielt = maillage->pos_per_fac_par[dbloc + 2*iper + 1] - 1;
             ielt < maillage->pos_per_fac_par[dbloc + 2*iper + 2] - 1;
             ielt++) {

          icelcr[ielt] = maillage->num_cel_fac_par[ielt];
          ipercr[ielt] = ((iper + 1) * 2);
        }

        /* Sens inverse */

        for (ielt = maillage->pos_per_fac_par[dbloc + 2*iper + 2] - 1;
             ielt < maillage->pos_per_fac_par[dbloc + 2*iper + 3] - 1;
             ielt++) {

          icelcr[ielt] = maillage->num_cel_fac_par[ielt];
          ipercr[ielt] = ((iper + 1) * 2) + 1;
        }

      } /* Fin boucle periodicites */

    } /* Fin domaine courant */

  }


  /* Uniquement dans le cas d'un voisinage étendu traité séparément
       et non vide sur le processeur courant */

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

    for (ind_loc = 0; ind_loc < maillage->nbr_dom_fac_par_voiset; ind_loc++) {

      dbloc = (2*maillage->nbr_per + 1)*ind_loc;

      /* On traite les cellules halos periodiques 'simples' du domaine*/

      if (maillage->nbr_dom == 1 ||
          maillage->num_dom_fac_par_voiset[ind_loc] == maillage->num_dom) {

        for (iper = 0; iper < maillage->nbr_per; iper++) {

          /* Sens direct */

          for (ielt = maillage->pos_per_fac_par_voiset[dbloc + 2*iper + 1] - 1;
               ielt < maillage->pos_per_fac_par_voiset[dbloc + 2*iper + 2] - 1;
               ielt++) {

            icelcr[ielt] = maillage->num_cel_fac_par_voiset[ielt];
            ipercr[ielt] = ((iper + 1) * 2);
          }

          /* Sens inverse */

          for (ielt = maillage->pos_per_fac_par_voiset[dbloc + 2*iper + 2] - 1;
               ielt < maillage->pos_per_fac_par_voiset[dbloc + 2*iper + 3] - 1;
               ielt++) {

            icelcr[ielt] = maillage->num_cel_fac_par_voiset[ielt];
            ipercr[ielt] = ((iper + 1) * 2) + 1;
          }

        } /* Fin boucle periodicites */

      } /* Fin domaine courant */

    } /* Fin boucle domaines */

  } /* Fin voisinage étendu séparé */

}


/*----------------------------------------------------------------------------
 * Retourne la position d'une particule après passage par une periodicite
 *
 * Interface Fortran :
 *
 * SUBROUTINE LAGPER
 * *****************
 *
 * INTEGER          NCELET        :  -> : nombre d'elements halo compris
 * INTEGER          NCEL          :  -> : nombre d'elements actifs
 * INTEGER          ISENS         :  -> : sens de la periodicite
 * INTEGER          IPER          :  -> : numero de la periodicite
 * INTEGER          IELT          :  -> : numero de la cellule de halo
 * DOUBLE PRECISION POINTA        :  -> : position avant transformation
 * DOUBLE PRECISION POINTB        : <-  : position après transformation
 *----------------------------------------------------------------------------*/

 
void CS_PROCF (lagper, LAGPER)(const cs_int_t   *ncelet,
                               const cs_int_t   *ncel,
                               const cs_int_t   *isens,
                               const cs_int_t   *iper, 
                               const cs_int_t   *ielt,
                               const cs_real_t   pointa[],
                                     cs_real_t   pointb[])
{
  cs_int_t i, j;
  cs_int_t ind_loc;
  cs_int_t dbloc;
  cs_int_t ind_elt_per;

  cs_real_t vect[CS_DIM_3];

  cs_maillage_t *maillage = cs_glob_maillage;
  cs_param_perio_t param_per;

  /*-----------------------------------*/
  /* Calcul des cdg des cellules halos */
  /*-----------------------------------*/

  param_per = maillage->liste_param_per[*iper];

  if ( *isens == 0 )  {

    /* Sens direct */

    ind_elt_per = maillage->num_cel_fac_par[*ielt] - 1;

    for (i = 0; i < CS_DIM_3; i++)
      vect[i] = 0.;

    for (i = 0; i < CS_DIM_3; i++)
      for (j = 0; j < CS_DIM_3; j++)
        vect[i] = vect[i] + param_per.matrice[i][j] *
                  (pointa[j] - param_per.point_inv[j]
                  + param_per.translation[j]);

     for (i = 0; i < CS_DIM_3; i++)
       pointb[i] =  vect[i] + param_per.point_inv[i];
  
  } else {

    /* Sens inverse */

    ind_elt_per = maillage->num_cel_fac_par[*ielt] - 1;

    for (i = 0; i < CS_DIM_3; i++)
      vect[i] = 0. ;

    for (i = 0; i < CS_DIM_3; i++)
      for (j = 0; j < CS_DIM_3; j++)
        vect[i] = vect[i] + param_per.matrice[j][i] *
                  (pointa[j] - param_per.point_inv[j]);

    for (i = 0; i < CS_DIM_3; i++)
      pointb[i]
      = vect[i] + param_per.point_inv[i] - param_per.translation[i];
  }
    
  /* Fin domaine courant */

}


/*----------------------------------------------------------------------------
 * Retourne la vitesse d'une particule après passage par une periodicite
 *
 * Interface Fortran :
 *
 * SUBROUTINE LAGVEC
 * *****************
 *
 * INTEGER          ISENS         :  -> : sens de la periodicite
 * INTEGER          IPER          :  -> : numero de la periodicite
 * DOUBLE PRECISION VECTI         :  -> : vecteur avant transformation
 * DOUBLE PRECISION VECTF         : <-  : vecteur après transformation
 *----------------------------------------------------------------------------*/

 
void CS_PROCF (lagvec, LAGVEC)(const cs_int_t   *isens,
                               const cs_int_t   *iper,  
                               const cs_real_t  vecti[],
                                     cs_real_t  vectf[])
{
  cs_maillage_t    *maillage = cs_glob_maillage;
  cs_param_perio_t  param_per;

  /*-----------------------------------*/
  /* Calcul des cdg des cellules halos */
  /*-----------------------------------*/

  param_per = maillage->liste_param_per[*iper];

  if ( *isens == 0 )  {

    /* Sens direct */

    vectf[0] = param_per.matrice[0][0] * vecti[0] +
               param_per.matrice[0][1] * vecti[1] +
               param_per.matrice[0][2] * vecti[2];
    vectf[1] = param_per.matrice[1][0] * vecti[0] +
               param_per.matrice[1][1] * vecti[1] +
               param_per.matrice[1][2] * vecti[2];
    vectf[2] = param_per.matrice[2][0] * vecti[0] +
               param_per.matrice[2][1] * vecti[1] +
               param_per.matrice[2][2] * vecti[2];

  } else {

    /* Sens inverse */

    vectf[0] = param_per.matrice[0][0] * vecti[0] +
               param_per.matrice[1][0] * vecti[1] +
               param_per.matrice[2][0] * vecti[2];
    vectf[1] = param_per.matrice[0][1] * vecti[0] +
               param_per.matrice[1][1] * vecti[1] +
               param_per.matrice[2][1] * vecti[2];
    vectf[2] = param_per.matrice[0][2] * vecti[0] +
               param_per.matrice[1][2] * vecti[1] +
               param_per.matrice[2][2] * vecti[2];

  } 
    
  /* Fin domaine courant */

}

#ifdef __cplusplus
}
#endif /* __cplusplus */
   
