/*============================================================================
*
*                    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 utilitaires pour le module diphasique Lagrangien
 *============================================================================*/

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

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

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

#include "cs_base.h"
#include "cs_msg.h"

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

#include "cs_lagrang.h"

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


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

/*----------------------------------------------------------------------------
 * Vérification que la norme IEEE 754 de stockage des flottants est respectée
 * par l'architecture, à l'aide de tests sur la précision du calcul.
 *----------------------------------------------------------------------------*/

cs_int_t cs_loc_lagrang__ieee
(
 void
);

/*----------------------------------------------------------------------------
 * The argument is an upperbound on the coordinates
 *----------------------------------------------------------------------------*/

void Orient3D_set_maxvalue
(
 cs_real_t d
);

/*----------------------------------------------------------------------------
 * Round a coordinate on the grid
 *----------------------------------------------------------------------------*/

void Orient3D_normalize
(
 float *value
);

/*----------------------------------------------------------------------------
 * Split a in two numbers. a=a1+a0, a1 get the most significant digit
 * format specify the range of the input, and how to split
 * if format is _split_i_j it means that a is of degree i
 * (in terms of original coordinates) the splitting bit is MAX0^i*2^j
 *----------------------------------------------------------------------------*/

/* inline is not an ansi feature and not recognized everywhere
 * BUT it looks as if -O optimization level lead to inlining anyway on Linux
 *                                                                   machines */

#if defined(__GNUC__)
__inline__ void Orient3D_split
#else
void Orient3D_split
#endif
(
cs_real_t  a      ,
cs_real_t *a1     ,
cs_real_t *a0     ,
cs_real_t  format
);

/*----------------------------------------------------------------------------
 * Compute the orientation of the four points
 * assumed to be rounded on the grid.
 * The last argument specify if the perturbation technique is used
 * in case of coplanarity
 *----------------------------------------------------------------------------*/


cs_int_t Orientation3D
(
  float ax   , float ay   , float az   ,
  float bx   , float by   , float bz   ,
  float cx   , float cy   , float cz   ,
  float dx   , float dy   , float dz   ,
  cs_int_t PERTURBED
);


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

/*----------------------------------------------------------------------------
 * Vérification que la norme IEEE 754 est respectée par l'architecture.
 * Dans le cas contraire la trajectographie des particules du module
 * Lagrangien est suceptible de ne pas fonctionner correctement.
 *----------------------------------------------------------------------------*/

void CS_PROCF (csieee,CSIEEE)
(
 void
)
{

  cs_int_t error  = 0 ;

  if (cs_loc_lagrang__ieee()) {

    cs_msg_flush() ;

    cs_msg_err (__FILE__, __LINE__, CS_MSG_TYPE_ERR, error,
    _("L'arithmétique IEEE 754 n'est pas respectée par l'architecture courante")) ;

  }

  return ;
}

/*----------------------------------------------------------------------------
 * Verification de la position relative de deux points, on cherche a savoir
 * s'ils sont confondus ou non.
 *----------------------------------------------------------------------------*/

void CS_PROCF (coloca,COLOCA)
(
  cs_real_t * pvalmax ,
  cs_real_t * px  , cs_real_t * py  , cs_real_t * pz  ,
  cs_real_t * qx  , cs_real_t * qy  , cs_real_t * qz  ,
  cs_int_t  * sign
)
{

  float xx0  ;  float yy0  ;  float zz0  ;
  float xx1  ;  float yy1  ;  float zz1  ;

  xx0 = (float) *px ;
  yy0 = (float) *py ;
  zz0 = (float) *pz ;

  xx1 = (float) *qx ;
  yy1 = (float) *qy ;
  zz1 = (float) *qz ;

/* on recherche la puissance de 2 la plus proche et superieure a valmax */

  Orient3D_set_maxvalue(*pvalmax) ;

  Orient3D_normalize(&xx0) ;
  Orient3D_normalize(&yy0) ;
  Orient3D_normalize(&zz0) ;

  Orient3D_normalize(&xx1) ;
  Orient3D_normalize(&yy1) ;
  Orient3D_normalize(&zz1) ;

/* on verifie si les 2 points passes en argument sont confondus */

  if ( (xx0==xx1) && (yy0==yy1) && (zz0==zz1) ) {
    *sign = 1 ;
  }
  else {
    *sign = 0 ;
  }


}

/*----------------------------------------------------------------------------
 * Recherche d'orientation des triedres pour le reperage des particules
 * par rapport aux faces
 *----------------------------------------------------------------------------*/

void CS_PROCF (coturn,COTURN)
(
  cs_real_t * pvalmax ,
  cs_real_t * px    , cs_real_t * py    , cs_real_t * pz    ,
  cs_real_t * qx    , cs_real_t * qy    , cs_real_t * qz    ,
  cs_real_t * cdgx  , cs_real_t * cdgy  , cs_real_t * cdgz  ,
  cs_real_t * crdx  , cs_real_t * crdy  , cs_real_t * crdz  ,
  cs_int_t  * sign  , cs_int_t  * pturb
)
{

  float xx0  ;  float yy0  ;  float zz0  ;
  float xx1  ;  float yy1  ;  float zz1  ;
  float xx2  ;  float yy2  ;  float zz2  ;
  float xx3  ;  float yy3  ;  float zz3  ;

/* cast en serie */

  xx0 = (float) * px ;
  yy0 = (float) * py ;
  zz0 = (float) * pz ;

  xx1 = (float) * qx ;
  yy1 = (float) * qy ;
  zz1 = (float) * qz ;

  xx2 = (float) * cdgx ;
  yy2 = (float) * cdgy ;
  zz2 = (float) * cdgz ;

  xx3 = (float) * crdx ;
  yy3 = (float) * crdy ;
  zz3 = (float) * crdz ;


/* on recherche la puissance de 2 la plus proche et superieur a valmax */

  Orient3D_set_maxvalue(*pvalmax) ;

  Orient3D_normalize(&xx0) ;
  Orient3D_normalize(&yy0) ;
  Orient3D_normalize(&zz0) ;

  Orient3D_normalize(&xx1) ;
  Orient3D_normalize(&yy1) ;
  Orient3D_normalize(&zz1) ;

  Orient3D_normalize(&xx2) ;
  Orient3D_normalize(&yy2) ;
  Orient3D_normalize(&zz2) ;

  Orient3D_normalize(&xx3) ;
  Orient3D_normalize(&yy3) ;
  Orient3D_normalize(&zz3) ;

  *sign = Orientation3D (xx0 , yy0 , zz0 ,
                         xx1 , yy1 , zz1 ,
                         xx2 , yy2 , zz2 ,
                         xx3 , yy3 , zz3 ,
                         *pturb          ) ;

}

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

/*----------------------------------------------------------------------------
 * Vérification que la norme IEEE 754 de stockage des flottants est respectée
 * par l'architecture, à l'aide de tests sur la précision du calcul.
 *----------------------------------------------------------------------------*/

cs_int_t cs_loc_lagrang__ieee
(
 void
)
{

  cs_real_t d = 9007199254740992.0 ;  /* 2^53 */
  cs_real_t dd ;

  dd = d + 1 ;

  if (dd != d) return EXIT_FAILURE ;  /* compute with too much bits */

  dd = d - 1 ;

  if (dd == d) return EXIT_FAILURE ;  /* compute with not enough bits */

  return EXIT_SUCCESS ;

}

/*----------------------------------------------------------------------------
 * The argument is an upperbound on the coordinates
 *----------------------------------------------------------------------------*/

static cs_real_t MAX0;                /* precision on coordinates                  */
static cs_real_t MAX24;               /* bounding box size                         */
static cs_real_t SNAP_TO_GRID;        /* for rounding of coordinates               */
static cs_real_t Orient3D_split_2_25; /* split to bit MAX0^2 * 2^25                */
static cs_real_t C1_ORIENTATION_3D;   /* static orientation filter                 */
static cs_real_t C2_ORIENTATION_3D;   /* semi-static orientation filter            */
static cs_real_t C3_ORIENTATION_3D;   /* half degree 3 unit for semi-static filter */
static cs_real_t C4_PERTURB_OR_3D;    /* static filter for perturbation            */


/*                                  DATA RANGE      */
/*  change the data conditionning                   */
/*  parameter is an upper bound on data values      */

void Orient3D_set_maxvalue
(
 cs_real_t d
)
{
  float D ;

  MAX0   = 1.0 ;                                     /* 2^0 */
  MAX24  = 16777216.0;                               /* 2^24 */
  SNAP_TO_GRID = 6755399441055744.0;                 /* 3 * 2^51 */
  Orient3D_split_2_25 = 453347182355485940514816.0 ; /* 3 * 2^77 */
  C1_ORIENTATION_3D = 75497472.0 ;                   /* 9 * 2^23 */
  C2_ORIENTATION_3D = 3.0/9007199254740992.0 ;       /* 3 * 2^(-53) */
  C3_ORIENTATION_3D = 0.5 ;                          /* 1/2 */
  C4_PERTURB_OR_3D  = 5066549580791808.0 ;           /* 9 * 2^49 */

  D = d * 4503599627370497.0 ;         /* 2^52 + 1 */
  D = ( d + D ) -D ;                   /* get the power of two closer to d */
  if ( D < ((float)d) ) D *= 2 ;       /* grid max size 2^b */

  MAX0  = D / 16777216.0 ;             /* grid step 2^k = 2^b / 2^24 */
  MAX24                 = D ;
  SNAP_TO_GRID         *= MAX0 ;
  Orient3D_split_2_25  *= MAX0*MAX0 ;
  C1_ORIENTATION_3D    *= MAX0*MAX0*MAX0 ;
  C3_ORIENTATION_3D    *= MAX0*MAX0*MAX0 ;
  C4_PERTURB_OR_3D     *= MAX0*MAX0*MAX0*MAX0 ;
}

/*----------------------------------------------------------------------------
 * Round a coordinate on the grid
 *----------------------------------------------------------------------------*/

void Orient3D_normalize
(
 float *value
)
{

/* round a value with fixed precision */

 if ( (*value>MAX24) || (*value<-MAX24))
    {
    fprintf(stderr,"overflow |%g| > %g\nVerify the bounding box for your data\n",*value,MAX24) ;
    *value=0.0 ;
    exit(0) ;
    }
  else
    *value = ( *value + SNAP_TO_GRID ) - SNAP_TO_GRID ;
}

/*----------------------------------------------------------------------------
 * Split a in two numbers. a=a1+a0, a1 get the most significant digit
 * format specify the range of the input, and how to split
 * if format is _split_i_j it means that a is of degree i
 * (in terms of original coordinates) the splitting bit is MAX0^i * 2^j
 *----------------------------------------------------------------------------*/

#if defined(__GNUC__)
__inline__ void Orient3D_split
#else
void Orient3D_split
#endif
(
cs_real_t  a      ,
cs_real_t *a1     ,
cs_real_t *a0     ,
cs_real_t  format
)
{
 *a1 = a+format ;
 *a1-= format   ;
 *a0 = a-*a1    ;
}

/*----------------------------------------------------------------------------
 * Compute the orientation of the four points
 * assumed to be rounded on the grid.
 * The last argument specify if the perturbation technique is used
 * in case of coplanarity
 *----------------------------------------------------------------------------*/

/*                     PREDICATES */


static const cs_int_t POSITIVE   =  1 ;
static const cs_int_t NEGATIVE   = -1 ;
static const cs_int_t COPLANAR   =  0 ;


/*                     ORIENTATION PREDICATES */
/* return the orientation of abcd */

cs_int_t Orientation3D
(
  float ax , float ay , float az ,
  float bx , float by , float bz ,
  float cx , float cy , float cz ,
  float dx , float dy , float dz ,
  cs_int_t PERTURBED
)
{
  cs_real_t error ;
  cs_real_t M10 , M11 , M20 , M21 , M30 , M31 ;
  cs_real_t RA , RB , RC , RD ;
  cs_real_t RB0 , RB1 , RC0 , RC1 , RD0 , RD1 ;

  /* points are assumed to be distincts */

  cs_real_t A = (cs_real_t)bx - (cs_real_t)ax ;    /*    x y z   */
  cs_real_t B = (cs_real_t)by - (cs_real_t)ay ;
  cs_real_t C = (cs_real_t)bz - (cs_real_t)az ;    /*  | A B C | */
  cs_real_t D = (cs_real_t)cx - (cs_real_t)ax ;    /*  | D E F | */
  cs_real_t E = (cs_real_t)cy - (cs_real_t)ay ;    /*  | G H I | */
  cs_real_t F = (cs_real_t)cz - (cs_real_t)az ;
  cs_real_t G = (cs_real_t)dx - (cs_real_t)ax ;
  cs_real_t H = (cs_real_t)dy - (cs_real_t)ay ;
  cs_real_t I = (cs_real_t)dz - (cs_real_t)az ;

  /* minors computation */

  cs_real_t M1 = D*H - E*G ;
  cs_real_t M2 = A*H - B*G ;
  cs_real_t M3 = A*E - B*D ;

  /* determinant computation */

  cs_real_t det = (C*M1 - F*M2) + I*M3 ;

  if (det >   C1_ORIENTATION_3D ) return POSITIVE ;
  if (det <  -C1_ORIENTATION_3D ) return NEGATIVE ;

  /* det is small */
  /* Orientation 3D, continuing */

  error = (fabs(C*M1) + fabs(F*M2)) + fabs(I*M3) ;
  error *= C2_ORIENTATION_3D ;

  if (error < C3_ORIENTATION_3D) error = 0.0 ;
  if (det >   error) return POSITIVE ;
  if (det <  -error) return NEGATIVE ;

  /* det is very small, exact computation must be done */
  /* Orientation 3D, continuing */

  if (error != 0.0) /* otherwise, exact value 0 already certified */
    {

    Orient3D_split( M1, &M11, &M10, Orient3D_split_2_25 ) ;
    Orient3D_split( M2, &M21, &M20, Orient3D_split_2_25 ) ;
    Orient3D_split( M3, &M31, &M30, Orient3D_split_2_25 ) ;

    det  =  C*M11 - F*M21 + I*M31  ;
    det += (C*M10 - F*M20 + I*M30) ;     /* less significant */

    if (det>0) return POSITIVE ;
    if (det<0) return NEGATIVE ;
    }

  /* points are coplanar */

  if (! PERTURBED)  return COPLANAR ;

  /* start perturbation scheme */
  /* (x,y,x)  |-->(x+e^3(x^2+y^2+z^2), y+e^2(x^2+y^2+z^2), z +e(x^2+y^2+z^2))*/

  RA = ax*ax + ay*ay + az*az      ;
  RB = bx*bx + by*by + bz*bz - RA ;
  RC = cx*cx + cy*cy + cz*cz - RA ;
  RD = dx*dx + dy*dy + dz*dz - RA ;

  /* epsilon term */

  det = (RB*M1 - RC*M2) + RD*M3 ; /* reuse minors */     /* A B RB */
  if (det >  C4_PERTURB_OR_3D) return POSITIVE ;         /* D E RC */
  if (det < -C4_PERTURB_OR_3D) return NEGATIVE ;         /* G H RD */

  /* det is small, doing exact computation */

   /* modif EDF deb */
  Orient3D_split( M1, &M11, &M10, Orient3D_split_2_25 ) ;
  Orient3D_split( M2, &M21, &M20, Orient3D_split_2_25 ) ;
  Orient3D_split( M3, &M31, &M30, Orient3D_split_2_25 ) ;
   /* modif EDF fin */

  Orient3D_split ( RB, &RB1, &RB0, Orient3D_split_2_25 ) ;
  Orient3D_split ( RC, &RC1, &RC0, Orient3D_split_2_25 ) ;
  Orient3D_split ( RD, &RD1, &RD0, Orient3D_split_2_25 ) ;

  det  =  RB1*M11 - RC1*M21 + RD1*M31 ;
  det += (RB1*M10 + RB0*M11) - (RC1*M20 + RC0*M21) + (RD1*M30 + RD0*M31) ;
  det += (RB0*M10 - RC0*M20 + RD0*M30) ;      /* less significant */
  if (det>0) return POSITIVE ;
  if (det<0) return NEGATIVE ;

  /* points coplanar in special direction go to e^2 term */

  M1 = D*I - F*G ;                           /* A RB C */
  M2 = A*I - C*G ;                           /* D RC F */
  M3 = A*F - C*D ;                           /* G RD I */
  det = (-RB*M1 + RC*M2) - RD*M3 ;
  if (det >  C4_PERTURB_OR_3D) return POSITIVE ;
  if (det < -C4_PERTURB_OR_3D) return NEGATIVE ;

  /* det is small, doing exact computation */

  Orient3D_split ( M1, &M11, &M10, Orient3D_split_2_25 ) ;
  Orient3D_split ( M2, &M21, &M20, Orient3D_split_2_25 ) ;
  Orient3D_split ( M3, &M31, &M30, Orient3D_split_2_25 ) ;

  det  = -RB1*M11 + RC1*M21 - RD1*M31 ;
  det += -(RB1*M10 + RB0*M11) + (RC1*M20 + RC0*M21) - (RD1*M30 + RD0*M31) ;
  det += (-RB0*M10 + RC0*M20 - RD0*M30) ;      /* less significant */
  if (det>0) return POSITIVE ;
  if (det<0) return NEGATIVE ;

  /* points coplanar in special direction go to e^3 term */

  M1 = E*I - F*H ;                           /* RB B C */
  M2 = B*I - C*H ;                           /* RC E F */
  M3 = B*F - C*E ;                           /* RD H I */
  det = (RB*M1 - RC*M2) + RD*M3 ;
  if (det >  C4_PERTURB_OR_3D) return POSITIVE ;
  if (det < -C4_PERTURB_OR_3D) return NEGATIVE ;

  /* det is small, doing exact computation */

  Orient3D_split ( M1, &M11, &M10, Orient3D_split_2_25 ) ;
  Orient3D_split ( M2, &M21, &M20, Orient3D_split_2_25 ) ;
  Orient3D_split ( M3, &M31, &M30, Orient3D_split_2_25 ) ;
  det  =  RB1*M11 - RC1*M21 + RD1*M31 ;
  det += (RB1*M10 + RB0*M11) - (RC1*M20 + RC0*M21) + (RD1*M30 + RD0*M31) ;
  det += (RB0*M10 - RC0*M20 + RD0*M30) ;      /* less significant */
  if (det>0) return POSITIVE ;
  if (det<0) return NEGATIVE ;

  /* points are colinear */

  return COPLANAR ;

}

#ifdef __cplusplus
}
#endif /* __cplusplus */
