Back to index

salome-med  6.5.0
InterpolationUtils.hxx
Go to the documentation of this file.
00001 // Copyright (C) 2007-2012  CEA/DEN, EDF R&D
00002 //
00003 // This library is free software; you can redistribute it and/or
00004 // modify it under the terms of the GNU Lesser General Public
00005 // License as published by the Free Software Foundation; either
00006 // version 2.1 of the License.
00007 //
00008 // This library is distributed in the hope that it will be useful,
00009 // but WITHOUT ANY WARRANTY; without even the implied warranty of
00010 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00011 // Lesser General Public License for more details.
00012 //
00013 // You should have received a copy of the GNU Lesser General Public
00014 // License along with this library; if not, write to the Free Software
00015 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
00016 //
00017 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
00018 //
00019 
00020 #ifndef __INTERPOLATIONUTILS_HXX__
00021 #define __INTERPOLATIONUTILS_HXX__
00022 
00023 #include "INTERPKERNELDefines.hxx"
00024 #include "InterpKernelException.hxx"
00025 
00026 #include "NormalizedUnstructuredMesh.hxx"
00027 
00028 #include <deque>
00029 #include <map>
00030 #include <cmath>
00031 #include <string>
00032 #include <vector>
00033 #include <algorithm>
00034 #include <iostream>
00035 #include <limits>
00036 
00037 namespace INTERP_KERNEL
00038 {
00039   template<class ConnType, NumberingPolicy numPol>
00040   class OTT//OffsetToolTrait
00041   {
00042   };
00043   
00044   template<class ConnType>
00045   class OTT<ConnType,ALL_C_MODE>
00046   {
00047   public:
00048     static ConnType indFC(ConnType i) { return i; }
00049     static ConnType ind2C(ConnType i) { return i; }
00050     static ConnType conn2C(ConnType i) { return i; }
00051     static ConnType coo2C(ConnType i) { return i; }
00052   };
00053   
00054   template<class ConnType>
00055   class OTT<ConnType,ALL_FORTRAN_MODE>
00056   {
00057   public:
00058     static ConnType indFC(ConnType i) { return i+1; }
00059     static ConnType ind2C(ConnType i) { return i-1; }
00060     static ConnType conn2C(ConnType i) { return i-1; }
00061     static ConnType coo2C(ConnType i) { return i-1; }
00062   };
00063 
00064   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
00065   /*   calcul la surface d'un triangle                  */
00066   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
00067 
00068   inline double Surf_Tri(const double* P_1,const double* P_2,const double* P_3)
00069   {
00070     double A=(P_3[1]-P_1[1])*(P_2[0]-P_1[0])-(P_2[1]-P_1[1])*(P_3[0]-P_1[0]);
00071     double Surface = 0.5*fabs(A);
00072     return Surface;
00073   }
00074 
00075   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
00076   /*     fonction qui calcul le determinant            */
00077   /*      de deux vecteur(cf doc CGAL).                */
00078   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
00079 
00080   //fonction qui calcul le determinant des vecteurs: P3P1 et P3P2
00081   //(cf doc CGAL).
00082 
00083   inline double mon_determinant(const double* P_1,
00084                                 const double*  P_2,
00085                                 const double* P_3)
00086   {
00087     double mon_det=(P_1[0]-P_3[0])*(P_2[1]-P_3[1])-(P_2[0]-P_3[0])*(P_1[1]-P_3[1]);
00088     return mon_det;
00089   }
00090 
00091   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
00092   //calcul la norme du vecteur P1P2
00093 
00094   inline double norme_vecteur(const double* P_1,const double* P_2)
00095   {
00096     double X=P_1[0]-P_2[0];
00097     double Y=P_1[1]-P_2[1];
00098     double norme=sqrt(X*X+Y*Y);
00099     return norme;
00100   }
00101 
00102   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
00103   /*         calcul le cos et le sin de l'angle P1P2,P1P3     */
00104   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
00105 
00106   inline std::vector<double> calcul_cos_et_sin(const double* P_1,
00107                                                const double* P_2,
00108                                                const double* P_3)
00109   {
00110        
00111     std::vector<double> Vect;
00112     double P1_P2=norme_vecteur(P_1,P_2);
00113     double P2_P3=norme_vecteur(P_2,P_3);
00114     double P3_P1=norme_vecteur(P_3,P_1);
00115        
00116     double N=P1_P2*P1_P2+P3_P1*P3_P1-P2_P3*P2_P3;
00117     double D=2.0*P1_P2*P3_P1;
00118     double COS=N/D;
00119     if (COS>1.0) COS=1.0;
00120     if (COS<-1.0) COS=-1.0;
00121     Vect.push_back(COS);
00122     double V=mon_determinant(P_2,P_3,P_1);
00123     double D_1=P1_P2*P3_P1;
00124     double SIN=V/D_1;
00125     if (SIN>1.0) SIN=1.0;
00126     if (SIN<-1.0) SIN=-1.0;
00127     Vect.push_back(SIN);
00128        
00129     return Vect;
00130        
00131   }
00132 
00140   template<int SPACEDIM>
00141   inline void fillDualCellOfTri(const double *triIn, double *quadOut)
00142   {
00143     //1st point
00144     std::copy(triIn,triIn+SPACEDIM,quadOut);
00145     double tmp[SPACEDIM];
00146     std::transform(triIn,triIn+SPACEDIM,triIn+SPACEDIM,tmp,std::plus<double>());
00147     //2nd point
00148     std::transform(tmp,tmp+SPACEDIM,quadOut+SPACEDIM,std::bind2nd(std::multiplies<double>(),0.5));
00149     std::transform(tmp,tmp+SPACEDIM,triIn+2*SPACEDIM,tmp,std::plus<double>());
00150     //3rd point
00151     std::transform(tmp,tmp+SPACEDIM,quadOut+2*SPACEDIM,std::bind2nd(std::multiplies<double>(),1/3.));
00152     //4th point
00153     std::transform(triIn,triIn+SPACEDIM,triIn+2*SPACEDIM,tmp,std::plus<double>());
00154     std::transform(tmp,tmp+SPACEDIM,quadOut+3*SPACEDIM,std::bind2nd(std::multiplies<double>(),0.5));
00155   }
00156 
00164   template<int SPACEDIM>
00165   inline void fillDualCellOfPolyg(const double *polygIn, int nPtsPolygonIn, double *polygOut)
00166   {
00167     //1st point
00168     std::copy(polygIn,polygIn+SPACEDIM,polygOut);
00169     std::transform(polygIn,polygIn+SPACEDIM,polygIn+SPACEDIM,polygOut+SPACEDIM,std::plus<double>());
00170     //2nd point
00171     std::transform(polygOut+SPACEDIM,polygOut+2*SPACEDIM,polygOut+SPACEDIM,std::bind2nd(std::multiplies<double>(),0.5));
00172     double tmp[SPACEDIM];
00173     //
00174     for(int i=0;i<nPtsPolygonIn-2;i++)
00175       {
00176         std::transform(polygIn,polygIn+SPACEDIM,polygIn+(i+2)*SPACEDIM,tmp,std::plus<double>());
00177         std::transform(tmp,tmp+SPACEDIM,polygOut+(2*i+3)*SPACEDIM,std::bind2nd(std::multiplies<double>(),0.5));
00178         std::transform(polygIn+(i+1)*SPACEDIM,polygIn+(i+2)*SPACEDIM,tmp,tmp,std::plus<double>());
00179         std::transform(tmp,tmp+SPACEDIM,polygOut+(2*i+2)*SPACEDIM,std::bind2nd(std::multiplies<double>(),1./3.));
00180       }
00181   }
00182 
00183   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
00184   /*     calcul les coordonnees du barycentre d'un polygone   */ 
00185   /*     le vecteur en entree est constitue des coordonnees   */
00186   /*     des sommets du polygone                              */                             
00187   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
00188 
00189   inline std::vector<double> bary_poly(const std::vector<double>& V)
00190   {
00191     std::vector<double> Bary;
00192     long taille=V.size();
00193     double x=0;
00194     double y=0;
00195 
00196     for(long i=0;i<taille/2;i++)
00197       {
00198         x=x+V[2*i];
00199         y=y+V[2*i+1];
00200       }
00201     double A=2*x/((double)taille);
00202     double B=2*y/((double)taille);
00203     Bary.push_back(A);//taille vecteur=2*nb de points.
00204     Bary.push_back(B);
00205 
00206 
00207     return Bary;
00208   }
00209   
00213   inline double computeTria6RefBase(const double *coeffs, const double *pos)
00214   {
00215     return coeffs[0]+coeffs[1]*pos[0]+coeffs[2]*pos[1]+coeffs[3]*pos[0]*pos[0]+coeffs[4]*pos[0]*pos[1]+coeffs[5]*pos[1]*pos[1];
00216   }
00217   
00221   inline void computeWeightedCoeffsInTria6FromRefBase(const double *refCoo, double *weightedPos)
00222   {
00223     weightedPos[0]=(1.-refCoo[0]-refCoo[1])*(1.-2*refCoo[0]-2.*refCoo[1]);
00224     weightedPos[1]=refCoo[0]*(2.*refCoo[0]-1.);
00225     weightedPos[2]=refCoo[1]*(2.*refCoo[1]-1.);
00226     weightedPos[3]=4.*refCoo[0]*(1.-refCoo[0]-refCoo[1]);
00227     weightedPos[4]=4.*refCoo[0]*refCoo[1];
00228     weightedPos[5]=4.*refCoo[1]*(1.-refCoo[0]-refCoo[1]);
00229   }
00230 
00234   inline double computeTetra10RefBase(const double *coeffs, const double *pos)
00235   {
00236     return coeffs[0]+coeffs[1]*pos[0]+coeffs[2]*pos[1]+coeffs[3]*pos[2]+
00237       coeffs[4]*pos[0]*pos[0]+coeffs[5]*pos[0]*pos[1]+coeffs[6]*pos[0]*pos[2]+
00238       coeffs[7]*pos[1]*pos[1]+coeffs[8]*pos[1]*pos[2]+coeffs[9]*pos[2]*pos[2];
00239   }
00240 
00244   inline void computeWeightedCoeffsInTetra10FromRefBase(const double *refCoo, double *weightedPos)
00245   {
00246     //http://www.cadfamily.com/download/CAE/ABAQUS/The%20Finite%20Element%20Method%20-%20A%20practical%20course%20abaqus.pdf page 217
00247     //L1=1-refCoo[0]-refCoo[1]-refCoo[2]
00248     //L2=refCoo[0] L3=refCoo[1] L4=refCoo[2]
00249     weightedPos[0]=(-2.*(refCoo[0]+refCoo[1]+refCoo[2])+1)*(1-refCoo[0]-refCoo[1]-refCoo[2]);//(2*L1-1)*L1
00250     weightedPos[1]=(2.*refCoo[0]-1.)*refCoo[0];//(2*L2-1)*L2
00251     weightedPos[2]=(2.*refCoo[1]-1.)*refCoo[1];//(2*L3-1)*L3
00252     weightedPos[3]=(2.*refCoo[2]-1.)*refCoo[2];//(2*L4-1)*L4
00253     weightedPos[4]=4.*(1-refCoo[0]-refCoo[1]-refCoo[2])*refCoo[0];//4*L1*L2
00254     weightedPos[5]=4.*refCoo[0]*refCoo[1];//4*L2*L3
00255     weightedPos[6]=4.*(1-refCoo[0]-refCoo[1]-refCoo[2])*refCoo[1];//4*L1*L3
00256     weightedPos[7]=4.*(1-refCoo[0]-refCoo[1]-refCoo[2])*refCoo[2];//4*L1*L4
00257     weightedPos[8]=4.*refCoo[0]*refCoo[2];//4*L2*L4
00258     weightedPos[9]=4.*refCoo[1]*refCoo[2];//4*L3*L4
00259   }
00260 
00267   template<unsigned nbRow>
00268   bool solveSystemOfEquations(double M[nbRow][nbRow+1], double* sol)
00269   {
00270     const int nbCol=nbRow+1;
00271 
00272     // make upper triangular matrix (forward elimination)
00273 
00274     int iR[nbRow];// = { 0, 1, 2 };
00275     for ( int i = 0; i < (int) nbRow; ++i )
00276       iR[i] = i;
00277     for ( int i = 0; i < (int)(nbRow-1); ++i ) // nullify nbRow-1 rows
00278       {
00279         // swap rows to have max value of i-th column in i-th row
00280         double max = std::fabs( M[ iR[i] ][i] );
00281         for ( int r = i+1; r < (int)nbRow; ++r )
00282           {
00283             double m = std::fabs( M[ iR[r] ][i] );
00284             if ( m > max )
00285               {
00286                 max = m;
00287                 std::swap( iR[r], iR[i] );
00288               }
00289           }
00290         if ( max < std::numeric_limits<double>::min() )
00291           {
00292             //sol[0]=1; sol[1]=sol[2]=sol[3]=0;
00293             return false; // no solution
00294           }
00295         // make 0 below M[i][i] (actually we do not modify i-th column)
00296         double* tUpRow = M[ iR[i] ];
00297         for ( int r = i+1; r < (int)nbRow; ++r )
00298           {
00299             double* mRow = M[ iR[r] ];
00300             double coef = mRow[ i ] / tUpRow[ i ];
00301             for ( int c = i+1; c < nbCol; ++c )
00302               mRow[ c ] -= tUpRow[ c ] * coef;
00303           }
00304       }
00305     double* mRow = M[ iR[nbRow-1] ];
00306     if ( std::fabs( mRow[ nbRow-1 ] ) < std::numeric_limits<double>::min() )
00307       {
00308         //sol[0]=1; sol[1]=sol[2]=sol[3]=0;
00309         return false; // no solution
00310       }
00311     mRow[ nbRow ] /= mRow[ nbRow-1 ];
00312 
00313     // calculate solution (back substitution)
00314 
00315     sol[ nbRow-1 ] = mRow[ nbRow ];
00316 
00317     for ( int i = nbRow-2; i+1; --i )
00318       {
00319         mRow = M[ iR[i] ];
00320         sol[ i ] = mRow[ nbRow ];
00321         for ( int j = nbRow-1; j > i; --j )
00322           sol[ i ] -= sol[j]*mRow[ j ];
00323         sol[ i ] /= mRow[ i ];
00324       }
00325 
00326     return true;
00327   }
00328 
00329   
00336   template<unsigned SZ, unsigned NB_OF_RES>
00337   bool solveSystemOfEquations2(const double *matrix, double *solutions, double eps)
00338   {
00339     unsigned k,j;
00340     int nr,n,m,np;
00341     double s,g;
00342     int mb;
00343     //
00344     double B[SZ*(SZ+NB_OF_RES)];
00345     std::copy(matrix,matrix+SZ*(SZ+NB_OF_RES),B);
00346     //
00347     nr=SZ+NB_OF_RES;
00348     for(k=0;k<SZ;k++)
00349       {
00350         np=nr*k+k;
00351         if(fabs(B[np])<eps)
00352           {
00353             n=k;
00354             do
00355               {
00356                 n++;
00357                 if(fabs(B[nr*k+n])>eps)
00358                   {/* Rows permutation */
00359                     for(m=0;m<nr;m++)
00360                       std::swap(B[nr*k+m],B[nr*n+m]);
00361                   }
00362               }
00363             while (n<(int)SZ);
00364           }
00365         s=B[np];//s is the Pivot
00366         std::transform(B+k*nr,B+(k+1)*nr,B+k*nr,std::bind2nd(std::divides<double>(),s));
00367         for(j=0;j<SZ;j++)
00368           {
00369             if(j!=k)
00370               {
00371                 g=B[j*nr+k];
00372                 for(mb=k;mb<nr;mb++)
00373                   B[j*nr+mb]-=B[k*nr+mb]*g;
00374               }
00375           }
00376       }
00377     for(j=0;j<NB_OF_RES;j++)
00378       for(k=0;k<SZ;k++)
00379         solutions[j*SZ+k]=B[nr*k+SZ+j];
00380     //
00381     return true;
00382   }
00383 
00384   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
00385   /*     Calculate barycentric coordinates of a 2D point p */ 
00386   /*     with respect to the triangle verices.             */
00387   /*     triaCoords are in full interlace                  */
00388   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
00389 
00390   template<int SPACEDIM>
00391   inline void barycentric_coords(const double* triaCoords, const double* p, double* bc)
00392   {
00393     // matrix 2x2
00394     double
00395       T11 = triaCoords[0]-triaCoords[2*SPACEDIM], T12 = triaCoords[SPACEDIM]-triaCoords[2*SPACEDIM],
00396       T21 = triaCoords[1]-triaCoords[2*SPACEDIM+1], T22 = triaCoords[SPACEDIM+1]-triaCoords[2*SPACEDIM+1];
00397     // matrix determinant
00398     double Tdet = T11*T22 - T12*T21;
00399     if ( fabs( Tdet ) < std::numeric_limits<double>::min() ) {
00400       bc[0]=1; bc[1]=0; bc[2]=0;
00401       return;
00402     }
00403     // matrix inverse
00404     double t11 = T22, t12 = -T12, t21 = -T21, t22 = T11;
00405     // vector
00406     double r11 = p[0]-triaCoords[2*SPACEDIM], r12 = p[1]-triaCoords[2*SPACEDIM+1];
00407     // barycentric coordinates: mutiply matrix by vector
00408     bc[0] = (t11 * r11 + t12 * r12)/Tdet;
00409     bc[1] = (t21 * r11 + t22 * r12)/Tdet;
00410     bc[2] = 1. - bc[0] - bc[1];
00411   }
00412 
00420   inline void barycentric_coords(const std::vector<const double*>& n, const double *p, double *bc)
00421   {
00422     enum { _X, _Y, _Z };
00423     switch(n.size())
00424       {
00425       case 2:
00426         {// SEG 2
00427           double delta=n[0][0]-n[1][0];
00428           bc[0]=fabs((*p-n[1][0])/delta);
00429           bc[1]=fabs((*p-n[0][0])/delta);
00430           break;
00431         }
00432       case 3:
00433         { // TRIA3
00434           // matrix 2x2
00435           double
00436             T11 = n[0][_X]-n[2][_X], T12 = n[1][_X]-n[2][_X],
00437             T21 = n[0][_Y]-n[2][_Y], T22 = n[1][_Y]-n[2][_Y];
00438           // matrix determinant
00439           double Tdet = T11*T22 - T12*T21;
00440           if ( std::fabs( Tdet ) < std::numeric_limits<double>::min() )
00441             {
00442               bc[0]=1; bc[1]=bc[2]=0; // no solution
00443               return;
00444             }
00445           // matrix inverse
00446           double t11 = T22, t12 = -T12, t21 = -T21, t22 = T11;
00447           // vector
00448           double r11 = p[_X]-n[2][_X], r12 = p[_Y]-n[2][_Y];
00449           // barycentric coordinates: mutiply matrix by vector
00450           bc[0] = (t11 * r11 + t12 * r12)/Tdet;
00451           bc[1] = (t21 * r11 + t22 * r12)/Tdet;
00452           bc[2] = 1. - bc[0] - bc[1];
00453           break;
00454         }
00455       case 4:
00456         { // TETRA4
00457           // Find bc by solving system of 3 equations using Gaussian elimination algorithm
00458           // bc1*( x1 - x4 ) + bc2*( x2 - x4 ) + bc3*( x3 - x4 ) = px - x4
00459           // bc1*( y1 - y4 ) + bc2*( y2 - y4 ) + bc3*( y3 - y4 ) = px - y4
00460           // bc1*( z1 - z4 ) + bc2*( z2 - z4 ) + bc3*( z3 - z4 ) = px - z4
00461           
00462           double T[3][4]=
00463             {{ n[0][_X]-n[3][_X], n[1][_X]-n[3][_X], n[2][_X]-n[3][_X], p[_X]-n[3][_X] },
00464              { n[0][_Y]-n[3][_Y], n[1][_Y]-n[3][_Y], n[2][_Y]-n[3][_Y], p[_Y]-n[3][_Y] },
00465              { n[0][_Z]-n[3][_Z], n[1][_Z]-n[3][_Z], n[2][_Z]-n[3][_Z], p[_Z]-n[3][_Z] }};
00466           
00467           if ( !solveSystemOfEquations<3>( T, bc ) )
00468             bc[0]=1., bc[1] = bc[2] = bc[3] = 0;
00469           else
00470             bc[ 3 ] = 1. - bc[0] - bc[1] - bc[2];
00471           break;
00472         }
00473       case 6:
00474         {
00475           // TRIA6
00476           double matrix2[48]={1., 0., 0., 0., 0., 0., 0., 0.,
00477                               1., 0., 0., 0., 0., 0., 1., 0., 
00478                               1., 0., 0., 0., 0., 0., 0., 1.,
00479                               1., 0., 0., 0., 0., 0., 0.5, 0.,
00480                               1., 0., 0., 0., 0., 0., 0.5, 0.5,
00481                               1., 0., 0., 0., 0., 0., 0.,0.5};
00482           for(int i=0;i<6;i++)
00483             {
00484               matrix2[8*i+1]=n[i][0];
00485               matrix2[8*i+2]=n[i][1];
00486               matrix2[8*i+3]=n[i][0]*n[i][0];
00487               matrix2[8*i+4]=n[i][0]*n[i][1];
00488               matrix2[8*i+5]=n[i][1]*n[i][1];
00489             }
00490           double res[12];
00491           solveSystemOfEquations2<6,2>(matrix2,res,std::numeric_limits<double>::min());
00492           double refCoo[2];
00493           refCoo[0]=computeTria6RefBase(res,p);
00494           refCoo[1]=computeTria6RefBase(res+6,p);
00495           computeWeightedCoeffsInTria6FromRefBase(refCoo,bc);
00496           break;
00497         }
00498       case 10:
00499         {//TETRA10
00500           double matrix2[130]={1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
00501                                1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.,
00502                                1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.,
00503                                1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.,
00504                                1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.5, 0., 0.,
00505                                1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.5, 0.5, 0.,
00506                                1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,0.5, 0.,
00507                                1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.5,
00508                                1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.5, 0., 0.5,
00509                                1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.5, 0.5};
00510           for(int i=0;i<10;i++)
00511             {
00512               matrix2[13*i+1]=n[i][0];
00513               matrix2[13*i+2]=n[i][1];
00514               matrix2[13*i+3]=n[i][2];
00515               matrix2[13*i+4]=n[i][0]*n[i][0];
00516               matrix2[13*i+5]=n[i][0]*n[i][1];
00517               matrix2[13*i+6]=n[i][0]*n[i][2];
00518               matrix2[13*i+7]=n[i][1]*n[i][1];
00519               matrix2[13*i+8]=n[i][1]*n[i][2];
00520               matrix2[13*i+9]=n[i][2]*n[i][2];
00521             }
00522           double res[30];
00523           solveSystemOfEquations2<10,3>(matrix2,res,std::numeric_limits<double>::min());
00524           double refCoo[3];
00525           refCoo[0]=computeTetra10RefBase(res,p);
00526           refCoo[1]=computeTetra10RefBase(res+10,p);
00527           refCoo[2]=computeTetra10RefBase(res+20,p);
00528           computeWeightedCoeffsInTetra10FromRefBase(refCoo,bc);
00529           break;
00530         }
00531       default:
00532         throw INTERP_KERNEL::Exception("INTERP_KERNEL::barycentric_coords : unrecognized simplex !");
00533       }
00534   }
00535 
00536   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
00537   /*         calcul la surface d'un polygone.                 */
00538   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
00539 
00540   inline double  Surf_Poly(const std::vector<double>& Poly)
00541   { 
00542 
00543     double Surface=0;
00544     for(unsigned long i=0; i<(Poly.size())/2-2; i++)
00545       {
00546         double Surf=Surf_Tri( &Poly[0],&Poly[2*(i+1)],&Poly[2*(i+2)] ); 
00547         Surface=Surface + Surf ;
00548       }
00549     return Surface ;
00550   }
00551 
00552   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
00553   /*   fonction qui teste si un point est dans une maille     */
00554   /*   point: P_0                                             */
00555   /*   P_1, P_2, P_3 sommet des mailles                       */   
00556   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
00557 
00558   inline bool point_dans_triangle(const double* P_0,const double* P_1,
00559                                   const double* P_2,const double* P_3,
00560                                   double eps)
00561   {
00562 
00563     bool A=false;
00564     double det_1=mon_determinant(P_1,P_3,P_0);
00565     double det_2=mon_determinant(P_3,P_2,P_0);
00566     double det_3=mon_determinant(P_2,P_1,P_0);
00567     if( (det_1>=-eps && det_2>=-eps && det_3>=-eps) || (det_1<=eps && det_2<=eps && det_3<=eps) )
00568       {
00569         A=true;
00570       }
00571 
00572     return A;
00573   }
00574 
00575   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
00576   /*fonction pour verifier qu'un point n'a pas deja ete considerer dans   */ 
00577   /*      le vecteur et le rajouter au vecteur sinon.                     */
00578   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
00579 
00580   inline void verif_point_dans_vect(const double* P, std::vector<double>& V, double absolute_precision )
00581   {
00582     long taille=V.size();
00583     bool isPresent=false;
00584     for(long i=0;i<taille/2;i++) 
00585       {
00586         if (sqrt(((P[0]-V[2*i])*(P[0]-V[2*i])+(P[1]-V[2*i+1])*(P[1]-V[2*i+1])))<absolute_precision)
00587           isPresent=true;
00588       
00589       }
00590     if(!isPresent)
00591       {
00592       
00593         V.push_back(P[0]);
00594         V.push_back(P[1]);    
00595       }
00596   }
00597 
00598   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
00599   /* fonction qui rajoute les sommet du triangle P dans le vecteur V        */ 
00600   /* si ceux-ci sont compris dans le triangle S et ne sont pas deja dans    */
00601   /* V.                                                                     */
00602   /*sommets de P: P_1, P_2, P_3                                             */
00603   /*sommets de S: P_4, P_5, P_6                                             */  
00604   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */ 
00605 
00606   inline void rajou_sommet_triangl(const double* P_1,const double* P_2,const double* P_3,
00607                                    const double* P_4,const double* P_5,const double* P_6,
00608                                    std::vector<double>& V, double dim_caracteristic, double precision)
00609   {
00610 
00611     double absolute_precision = precision*dim_caracteristic;
00612     bool A_1=INTERP_KERNEL::point_dans_triangle(P_1,P_4,P_5,P_6,absolute_precision);
00613     if(A_1)
00614       verif_point_dans_vect(P_1,V,absolute_precision);
00615     bool A_2=INTERP_KERNEL::point_dans_triangle(P_2,P_4,P_5,P_6,absolute_precision);
00616     if(A_2)
00617       verif_point_dans_vect(P_2,V,absolute_precision);
00618     bool A_3=INTERP_KERNEL::point_dans_triangle(P_3,P_4,P_5,P_6,absolute_precision);
00619     if(A_3)
00620       verif_point_dans_vect(P_3,V,absolute_precision);
00621   }
00622 
00623 
00624   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _  _ _ _ _ _ _ _ _*/
00625   /*  calcul de l'intersection de deux segments: segments P1P2 avec P3P4      */
00626   /*  . Si l'intersection est non nulle et si celle-ci n'est                  */
00627   /*  n'est pas deja contenue dans Vect on la rajoute a Vect                  */
00628   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _  _ _ _ _ _ _ _ _*/ 
00629 
00630   inline void inters_de_segment(const double * P_1,const double * P_2,
00631                                 const double * P_3,const double * P_4,
00632                                 std::vector<double>& Vect, 
00633                                 double dim_caracteristic, double precision)
00634   {
00635     // calcul du determinant de P_1P_2 et P_3P_4.
00636     double det=(P_2[0]-P_1[0])*(P_4[1]-P_3[1])-(P_4[0]-P_3[0])*(P_2[1]-P_1[1]);
00637 
00638     double absolute_precision = dim_caracteristic*precision;
00639     if(fabs(det)>absolute_precision)
00640       {
00641         double k_1=-((P_3[1]-P_4[1])*(P_3[0]-P_1[0])+(P_4[0]-P_3[0])*(P_3[1]-P_1[1]))/det;
00642 
00643         if (k_1 >= -absolute_precision && k_1 <= 1+absolute_precision)
00644           //if( k_1 >= -precision && k_1 <= 1+precision)
00645           {
00646             double k_2= ((P_1[1]-P_2[1])*(P_1[0]-P_3[0])+(P_2[0]-P_1[0])*(P_1[1]-P_3[1]))/det;
00647 
00648             if (k_2 >= -absolute_precision && k_2 <= 1+absolute_precision)
00649               //if( k_2 >= -precision && k_2 <= 1+precision)
00650               {
00651                 double P_0[2];
00652                 P_0[0]=P_1[0]+k_1*(P_2[0]-P_1[0]);
00653                 P_0[1]=P_1[1]+k_1*(P_2[1]-P_1[1]);
00654                 verif_point_dans_vect(P_0,Vect,absolute_precision);
00655               }
00656           }
00657       }
00658   }
00659 
00660 
00661 
00662   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/  
00663   /*      calcul l'intersection de deux triangles            */
00664   /* P_1, P_2, P_3: sommets du premier triangle              */
00665   /* P_4, P_5, P_6: sommets du deuxi�me triangle             */
00666   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/ 
00667 
00668   inline void intersec_de_triangle(const double* P_1,const double* P_2, const double* P_3,
00669                                    const double* P_4,const double* P_5,const double* P_6, 
00670                                    std::vector<double>& Vect, double dim_caracteristic, double precision)
00671   {
00672     inters_de_segment(P_1,P_2,P_4,P_5,Vect, dim_caracteristic, precision);
00673     inters_de_segment(P_1,P_2,P_5,P_6,Vect, dim_caracteristic, precision);
00674     inters_de_segment(P_1,P_2,P_6,P_4,Vect, dim_caracteristic, precision);
00675     inters_de_segment(P_2,P_3,P_4,P_5,Vect, dim_caracteristic, precision);
00676     inters_de_segment(P_2,P_3,P_5,P_6,Vect, dim_caracteristic, precision);
00677     inters_de_segment(P_2,P_3,P_6,P_4,Vect, dim_caracteristic, precision);
00678     inters_de_segment(P_3,P_1,P_4,P_5,Vect, dim_caracteristic, precision);
00679     inters_de_segment(P_3,P_1,P_5,P_6,Vect, dim_caracteristic, precision);
00680     inters_de_segment(P_3,P_1,P_6,P_4,Vect, dim_caracteristic, precision);
00681     rajou_sommet_triangl(P_1,P_2,P_3,P_4,P_5,P_6,Vect, dim_caracteristic, precision);
00682     rajou_sommet_triangl(P_4,P_5,P_6,P_1,P_2,P_3,Vect, dim_caracteristic, precision);
00683   }
00684 
00685   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
00686   /* fonction pour verifier qu'un node maille n'a pas deja ete considerer  */
00687   /*  dans le vecteur et le rajouter au vecteur sinon.                     */
00688   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
00689 
00690   inline void verif_maill_dans_vect(int Num, std::vector<int>& V)
00691   {
00692     long taille=V.size();
00693     int A=0;
00694     for(long i=0;i<taille;i++)
00695       {
00696         if(Num==V[i])
00697           {
00698             A=1;
00699             break;
00700           } 
00701       }
00702     if(A==0)
00703       {V.push_back(Num); }
00704   }
00705 
00708   class AngleLess
00709   {
00710   public:
00711     bool operator()(std::pair<double,double>theta1, std::pair<double,double> theta2) 
00712     {
00713       double norm1 = sqrt(theta1.first*theta1.first +theta1.second*theta1.second);
00714       double norm2 = sqrt(theta2.first*theta2.first +theta2.second*theta2.second);
00715       
00716       double epsilon = 1.e-12;
00717       
00718       if( norm1 < epsilon || norm2 < epsilon  ) 
00719         std::cout << "Warning InterpolationUtils.hxx: AngleLess : Vector with zero norm, cannot define the angle !!!! " << std::endl;
00720       
00721       return theta1.second*(norm2 + theta2.first) < theta2.second*(norm1 + theta1.first);
00722     
00723     }
00724   };
00725 
00726 
00727   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */  
00728   /* fonction pour reconstituer un polygone convexe a partir  */
00729   /*              d'un nuage de point.                        */
00730   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */  
00731 
00732   inline std::vector<double> reconstruct_polygon(const std::vector<double>& V)
00733   {
00734 
00735     std::size_t taille=V.size();
00736 
00737     //VB : why 6 ?
00738 
00739     if(taille<=6)
00740       {return V;}
00741     else
00742       {
00743         double *COS=new double[taille/2];
00744         double *SIN=new double[taille/2];
00745         //double *angle=new double[taille/2];
00746         std::vector<double> Bary=bary_poly(V);
00747         COS[0]=1.0;
00748         SIN[0]=0.0;
00749         //angle[0]=0.0;
00750         for(std::size_t i=0; i<taille/2-1;i++)
00751           {
00752             std::vector<double> Trigo=calcul_cos_et_sin(&Bary[0],&V[0],&V[2*(i+1)]);
00753             COS[i+1]=Trigo[0];
00754             SIN[i+1]=Trigo[1];
00755             //if(SIN[i+1]>=0)
00756             //    {angle[i+1]=atan2(SIN[i+1],COS[i+1]);}
00757             //             else
00758             //               {angle[i+1]=-atan2(SIN[i+1],COS[i+1]);}
00759           }
00760                      
00761         //ensuite on ordonne les angles.
00762         std::vector<double> Pt_ordonne;
00763         Pt_ordonne.reserve(taille);
00764         //        std::multimap<double,int> Ordre;
00765         std::multimap<std::pair<double,double>,int, AngleLess> CosSin;
00766         for(std::size_t i=0;i<taille/2;i++)       
00767           {
00768             //  Ordre.insert(std::make_pair(angle[i],i));
00769             CosSin.insert(std::make_pair(std::make_pair(SIN[i],COS[i]),i));
00770           }
00771         //        std::multimap <double,int>::iterator mi;
00772         std::multimap<std::pair<double,double>,int, AngleLess>::iterator   micossin;
00773         //         for(mi=Ordre.begin();mi!=Ordre.end();mi++)
00774         //           {
00775         //             int j=(*mi).second;
00776         //             Pt_ordonne.push_back(V[2*j]);
00777         //             Pt_ordonne.push_back(V[2*j+1]);
00778         //           }
00779         for(micossin=CosSin.begin();micossin!=CosSin.end();micossin++)
00780           {
00781             int j=(*micossin).second;
00782             Pt_ordonne.push_back(V[2*j]);
00783             Pt_ordonne.push_back(V[2*j+1]);
00784           }
00785         delete [] COS;
00786         delete [] SIN;
00787         //        delete [] angle;
00788         return Pt_ordonne;
00789       }
00790   }
00791 
00792   template<int DIM, NumberingPolicy numPol, class MyMeshType>
00793   inline void getElemBB(double* bb, const double *coordsOfMesh, int iP, int nb_nodes)
00794   {
00795     bb[0]=std::numeric_limits<double>::max();
00796     bb[1]=-std::numeric_limits<double>::max();
00797     bb[2]=std::numeric_limits<double>::max();
00798     bb[3]=-std::numeric_limits<double>::max();
00799     bb[4]=std::numeric_limits<double>::max();
00800     bb[5]=-std::numeric_limits<double>::max();
00801     
00802     for (int i=0; i<nb_nodes; i++)
00803       {
00804         double x = coordsOfMesh[3*(iP+i)];
00805         double y = coordsOfMesh[3*(iP+i)+1];
00806         double z = coordsOfMesh[3*(iP+i)+2];
00807         bb[0]=(x<bb[0])?x:bb[0];
00808         bb[1]=(x>bb[1])?x:bb[1];
00809         bb[2]=(y<bb[2])?y:bb[2];
00810         bb[3]=(y>bb[3])?y:bb[3];
00811         bb[4]=(z<bb[4])?z:bb[4];
00812         bb[5]=(z>bb[5])?z:bb[5];
00813       }              
00814   }
00815 
00816   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
00817   /* Computes the dot product of a and b */
00818   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
00819   template<int dim> 
00820   inline double dotprod( const double * a, const double * b)
00821   {
00822     double result=0;
00823     for(int idim = 0; idim < dim ; idim++) result += a[idim]*b[idim];
00824     return result;
00825   }
00826   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
00827   /* Computes the norm of vector v */
00828   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/  
00829   template<int dim> 
00830   inline double norm(const double * v)
00831   {   
00832     double result =0;
00833     for(int idim =0; idim<dim; idim++) result+=v[idim]*v[idim];
00834     return sqrt(result);
00835   }
00836   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
00837   /* Computes the square norm of vector a-b */
00838   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/  
00839   template<int dim> 
00840   inline double distance2( const double * a, const double * b)
00841   {   
00842     double result =0;
00843     for(int idim =0; idim<dim; idim++) result+=(a[idim]-b[idim])*(a[idim]-b[idim]);
00844     return result;
00845   }
00846   template<class T, int dim> 
00847   inline double distance2(  T * a, int inda, T * b, int indb)
00848   {   
00849     double result =0;
00850     for(int idim =0; idim<dim; idim++) result += ((*a)[inda+idim] - (*b)[indb+idim])* ((*a)[inda+idim] - (*b)[indb+idim]);
00851     return result;
00852   }
00853   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
00854   /* Computes the determinant of a and b */
00855   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
00856   inline double determinant ( double *  a, double * b)
00857   {
00858     return a[0]*b[1]-a[1]*b[0];
00859   }
00860   inline double determinant ( double *  a, double * b, double * c)
00861   {
00862     return a[0]*determinant(b+1,c+1)-b[0]*determinant(a+1,c+1)+c[0]*determinant(a+1,b+1);
00863   }
00864   
00865   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
00866   /* Computes the cross product of AB and AC */
00867   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/  
00868 
00869   template<int dim> inline void crossprod( const double * A, const double * B, const double * C, double * V);
00870   
00871   template<> inline
00872   void crossprod<2>( const double * A, const double * B, const double * C, double * V)
00873   {   
00874     double AB[2];
00875     double AC[2];
00876     for(int idim =0; idim<2; idim++) AB[idim] = B[idim]-A[idim];//B-A
00877     for(int idim =0; idim<2; idim++) AC[idim] = C[idim]-A[idim];//C-A;
00878 
00879     V[0]=determinant(AB,AC);
00880     V[1]=0;
00881   }
00882   template<> inline
00883   void crossprod<3>( const double * A, const double * B, const double * C, double * V)
00884   {   
00885     double AB[3];
00886     double AC[3];
00887     for(int idim =0; idim<3; idim++) AB[idim] = B[idim]-A[idim];//B-A
00888     for(int idim =0; idim<3; idim++) AC[idim] = C[idim]-A[idim];//C-A;
00889 
00890     V[0]=AB[1]*AC[2]-AB[2]*AC[1];
00891     V[1]=-AB[0]*AC[2]+AB[2]*AC[0];
00892     V[2]=AB[0]*AC[1]-AB[1]*AC[0];    
00893   }
00894   template<> inline
00895   void crossprod<1>( const double * /*A*/, const double * /*B*/, const double * /*C*/, double * /*V*/)
00896   {
00897     // just to be able to compile
00898   }
00899   
00900   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
00901   /* Checks wether point A is inside the quadrangle BCDE */
00902   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/  
00903 
00904   template<int dim> inline double check_inside(const double* A,const double* B,const double* C,const double* D,
00905                                                const double* E,double* ABC, double* ADE)
00906   {
00907     crossprod<dim>(A,B,C,ABC);
00908     crossprod<dim>(A,D,E,ADE);
00909     return dotprod<dim>(ABC,ADE);
00910   }   
00911 
00912   
00913   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
00914   /* Computes the geometric angle (in [0,Pi]) between two non zero vectors AB and AC */
00915   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/  
00916   template<int dim> inline double angle(const double * A, const double * B, const double * C, double * n)
00917   {   
00918     double AB[dim];
00919     double AC[dim];
00920     double orthAB[dim];
00921 
00922     for(int idim =0; idim<dim; idim++) AB[idim] = B[idim]-A[idim];//B-A;
00923     for(int idim =0; idim<dim; idim++) AC[idim] = C[idim]-A[idim];//C-A;
00924 
00925     double normAB= norm<dim>(AB);
00926     for(int idim =0; idim<dim; idim++) AB[idim]/=normAB;
00927 
00928     double normAC= norm<dim>(AC);
00929     double AB_dot_AC=dotprod<dim>(AB,AC);
00930     for(int idim =0; idim<dim; idim++) orthAB[idim] = AC[idim]-AB_dot_AC*AB[idim];
00931 
00932     double denom= normAC+AB_dot_AC;
00933     double numer=norm<dim>(orthAB);
00934     
00935     return 2*atan2(numer,denom);
00936   }    
00937   
00938   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
00939   /* Tells whether the frame constituted of vectors AB, AC and n is direct */
00940   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/  
00941   template<int dim> inline double direct_frame(const double * A, const double * B, const double * C, double * n);
00942   template<> inline
00943   double direct_frame<2>(const double * A, const double * B, const double * C, double * n)
00944   {
00945     double AB[2];
00946     double AC[2];
00947     for(int idim =0; idim<2; idim++) AB[idim] = B[idim]-A[idim];//B-A;
00948     for(int idim =0; idim<2; idim++) AC[idim] = C[idim]-A[idim];//C-A;
00949     
00950     return  determinant(AB,AC)*n[0];
00951   }
00952   template<> inline
00953   double direct_frame<3>(const double * A, const double * B, const double * C, double * n)
00954   {
00955     double AB[3];
00956     double AC[3];
00957     for(int idim =0; idim<3; idim++) AB[idim] = B[idim]-A[idim];//B-A;
00958     for(int idim =0; idim<3; idim++) AC[idim] = C[idim]-A[idim];//C-A;
00959     
00960     return determinant(AB,AC,n)>0;
00961   }
00962 
00963   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/  
00964   /*      calcul l'intersection de deux polygones COPLANAIRES */
00965   /* en dimension DIM (2 ou 3). Si DIM=3 l'algorithme ne considere*/
00966   /* que les deux premieres coordonnees de chaque point */
00967   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/ 
00968   template<int DIM> inline void intersec_de_polygone(const double * Coords_A, const double * Coords_B, 
00969                                                      int nb_NodesA, int nb_NodesB,
00970                                                      std::vector<double>& inter,
00971                                                      double dim_caracteristic, double precision)
00972   {
00973     for(int i_A = 1; i_A<nb_NodesA-1; i_A++)
00974       {
00975         for(int i_B = 1; i_B<nb_NodesB-1; i_B++)
00976           {
00977             INTERP_KERNEL::intersec_de_triangle(&Coords_A[0],&Coords_A[DIM*i_A],&Coords_A[DIM*(i_A+1)],
00978                                                 &Coords_B[0],&Coords_B[DIM*i_B],&Coords_B[DIM*(i_B+1)],
00979                                                 inter, dim_caracteristic, precision);
00980           }
00981       }
00982     int nb_inter=((int)inter.size())/DIM;
00983     if(nb_inter >3) inter=INTERP_KERNEL::reconstruct_polygon(inter);
00984   }
00985 
00986   /*_ _ _ _ _ _ _ _ _
00987    *_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ 
00988    *  fonctions qui calcule l'aire d'un polygone en dimension 2 ou 3    
00989    *_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
00990   template<int DIM> inline double polygon_area(std::vector<double>& inter)
00991   {
00992     double result=0.;
00993     double area[DIM];
00994                   
00995     for(int i = 1; i<(int)inter.size()/DIM-1; i++)
00996       {
00997         INTERP_KERNEL::crossprod<DIM>(&inter[0],&inter[DIM*i],&inter[DIM*(i+1)],area);
00998         result +=0.5*norm<DIM>(area);
00999       }
01000     return result;
01001   }
01002          
01003   template<int DIM> inline double polygon_area(std::deque<double>& inter)
01004   {
01005     double result=0.;
01006     double area[DIM];
01007                   
01008     for(int i = 1; i<(int)inter.size()/DIM-1; i++)
01009       {
01010         INTERP_KERNEL::crossprod<DIM>(&inter[0],&inter[DIM*i],&inter[DIM*(i+1)],area);
01011         result +=0.5*norm<DIM>(area);
01012       }
01013     return result;
01014   }
01015   
01017   inline double triple_product(const double* A, const double*B, const double*C, const double*X)
01018   {
01019     double XA[3];
01020     XA[0]=A[0]-X[0];
01021     XA[1]=A[1]-X[1];
01022     XA[2]=A[2]-X[2];
01023     double XB[3];
01024     XB[0]=B[0]-X[0];
01025     XB[1]=B[1]-X[1];
01026     XB[2]=B[2]-X[2];
01027     double XC[3];
01028     XC[0]=C[0]-X[0];
01029     XC[1]=C[1]-X[1];
01030     XC[2]=C[2]-X[2];
01031     
01032     return 
01033       (XA[1]*XB[2]-XA[2]*XB[1])*XC[0]+
01034       (XA[2]*XB[0]-XA[0]*XB[2])*XC[1]+
01035       (XA[0]*XB[1]-XA[1]*XB[0])*XC[2];
01036   }
01037   
01041   template<class T, int dim> 
01042   bool checkEqualPolygonsOneDirection(T * L1, T * L2, int size1, int size2, int istart1, int istart2, double epsilon, int sign)
01043   {
01044     int i1 = istart1;
01045     int i2 = istart2;
01046     int i1next = ( i1 + 1 ) % size1;
01047     int i2next = ( i2 + sign +size2) % size2;
01048     
01049     while(true)
01050       {
01051         while( i1next!=istart1 && distance2<T,dim>(L1,i1*dim, L1,i1next*dim) < epsilon ) i1next = (  i1next + 1 ) % size1;  
01052         while( i2next!=istart2 && distance2<T,dim>(L2,i2*dim, L2,i2next*dim) < epsilon ) i2next = (  i2next + sign +size2 ) % size2;  
01053         
01054         if(i1next == istart1)
01055           {
01056             if(i2next == istart2)
01057               return true;
01058             else return false;
01059           }
01060         else
01061           if(i2next == istart2)
01062             return false;
01063           else 
01064             {
01065               if(distance2<T,dim>(L1,i1next*dim, L2,i2next*dim) > epsilon )
01066                 return false;
01067               else
01068                 {
01069                   i1 = i1next;
01070                   i2 = i2next;
01071                   i1next = ( i1 + 1 ) % size1;
01072                   i2next = ( i2 + sign + size2 ) % size2;
01073                 }
01074             }
01075       }
01076   }
01077 
01080   template<class T, int dim> 
01081   bool checkEqualPolygons(T * L1, T * L2, double epsilon)
01082   {
01083     if(L1==NULL || L2==NULL) 
01084       {
01085         std::cout << "Warning InterpolationUtils.hxx:checkEqualPolygonsPointer: Null pointer " << std::endl;
01086         throw(Exception("big error: not closed polygon..."));
01087       }
01088     
01089     int size1 = (*L1).size()/dim;
01090     int size2 = (*L2).size()/dim;
01091     int istart1 = 0;
01092     int istart2 = 0;
01093     
01094     while( istart2 < size2  && distance2<T,dim>(L1,istart1*dim, L2,istart2*dim) > epsilon ) istart2++;
01095   
01096     if(istart2 == size2)
01097       {  
01098         return (size1 == 0) && (size2 == 0);
01099       }
01100     else 
01101       return   checkEqualPolygonsOneDirection<T,dim>( L1, L2, size1, size2, istart1, istart2, epsilon,  1)
01102         || checkEqualPolygonsOneDirection<T,dim>( L1, L2, size1, size2, istart1, istart2, epsilon, -1);
01103 
01104   }
01105 }
01106 
01107 
01108 #endif