/*============================================================================
*
* 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 */
syntax highlighted by Code2HTML, v. 0.9.1