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