src/INTERP_KERNEL/InterpolationUtils.hxx

   1 // Copyright (C) 2007-2013  CEA/DEN, EDF R&D
   2 //
   3 // This library is free software; you can redistribute it and/or
   4 // modify it under the terms of the GNU Lesser General Public
   5 // License as published by the Free Software Foundation; either
   6 // version 2.1 of the License.
   7 //
   8 // This library is distributed in the hope that it will be useful,
   9 // but WITHOUT ANY WARRANTY; without even the implied warranty of
  10 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  11 // Lesser General Public License for more details.
  12 //
  13 // You should have received a copy of the GNU Lesser General Public
  14 // License along with this library; if not, write to the Free Software
  15 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA  02111-1307 USA
  16 //
  17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
  18 //
  19 // Author : Anthony Geay (CEA/DEN)
  20
  21 #ifndef __INTERPOLATIONUTILS_HXX__
  22 #define __INTERPOLATIONUTILS_HXX__
  23
  24 #include "INTERPKERNELDefines.hxx"
  25 #include "InterpKernelException.hxx"
  26
  27 #include "NormalizedUnstructuredMesh.hxx"
  28
  29 #include <deque>
  30 #include <map>
  31 #include <cmath>
  32 #include <string>
  33 #include <vector>
  34 #include <algorithm>
  35 #include <iostream>
  36 #include <limits>
  37
  38 namespace INTERP_KERNEL
  39 {
  40   template<class ConnType, NumberingPolicy numPol>
  41   class OTT//OffsetToolTrait
  42   {
  43   };
  44
  45   template<class ConnType>
  46   class OTT<ConnType,ALL_C_MODE>
  47   {
  48   public:
  49     static ConnType indFC(ConnType i) { return i; }
  50     static ConnType ind2C(ConnType i) { return i; }
  51     static ConnType conn2C(ConnType i) { return i; }
  52     static ConnType coo2C(ConnType i) { return i; }
  53   };
  54
  55   template<class ConnType>
  56   class OTT<ConnType,ALL_FORTRAN_MODE>
  57   {
  58   public:
  59     static ConnType indFC(ConnType i) { return i+1; }
  60     static ConnType ind2C(ConnType i) { return i-1; }
  61     static ConnType conn2C(ConnType i) { return i-1; }
  62     static ConnType coo2C(ConnType i) { return i-1; }
  63   };
  64
  65   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
  66   /*   calcul la surface d'un triangle                  */
  67   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
  68
  69   inline double Surf_Tri(const double* P_1,const double* P_2,const double* P_3)
  70   {
  71     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]);
  72     double Surface = 0.5*fabs(A);
  73     return Surface;
  74   }
  75
  76   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
  77   /*     fonction qui calcul le determinant            */
  78   /*      de deux vecteur(cf doc CGAL).                */
  79   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
  80
  81   //fonction qui calcul le determinant des vecteurs: P3P1 et P3P2
  82   //(cf doc CGAL).
  83
  84   inline double mon_determinant(const double* P_1,
  85                                 const double*  P_2,
  86                                 const double* P_3)
  87   {
  88     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]);
  89     return mon_det;
  90   }
  91
  92   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
  93   //calcul la norme du vecteur P1P2
  94
  95   inline double norme_vecteur(const double* P_1,const double* P_2)
  96   {
  97     double X=P_1[0]-P_2[0];
  98     double Y=P_1[1]-P_2[1];
  99     double norme=sqrt(X*X+Y*Y);
 100     return norme;
 101   }
 102
 103   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
 104   /*         calcul le cos et le sin de l'angle P1P2,P1P3     */
 105   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
 106
 107   inline std::vector<double> calcul_cos_et_sin(const double* P_1,
 108                                                const double* P_2,
 109                                                const double* P_3)
 110   {
 111
 112     std::vector<double> Vect;
 113     double P1_P2=norme_vecteur(P_1,P_2);
 114     double P2_P3=norme_vecteur(P_2,P_3);
 115     double P3_P1=norme_vecteur(P_3,P_1);
 116
 117     double N=P1_P2*P1_P2+P3_P1*P3_P1-P2_P3*P2_P3;
 118     double D=2.0*P1_P2*P3_P1;
 119     double COS=N/D;
 120     if (COS>1.0) COS=1.0;
 121     if (COS<-1.0) COS=-1.0;
 122     Vect.push_back(COS);
 123     double V=mon_determinant(P_2,P_3,P_1);
 124     double D_1=P1_P2*P3_P1;
 125     double SIN=V/D_1;
 126     if (SIN>1.0) SIN=1.0;
 127     if (SIN<-1.0) SIN=-1.0;
 128     Vect.push_back(SIN);
 129
 130     return Vect;
 131
 132   }
 133
 134   /*!
 135    * This method builds a quadrangle built with the first point of 'triIn' the barycenter of two edges starting or ending with
 136    * the first point of 'triIn' and the barycenter of 'triIn'.
 137    *
 138    * @param triIn is a 6 doubles array in full interlace mode, that represents a triangle.
 139    * @param quadOut is a 8 doubles array filled after the following call.
 140    */
 141   template<int SPACEDIM>
 142   inline void fillDualCellOfTri(const double *triIn, double *quadOut)
 143   {
 144     //1st point
 145     std::copy(triIn,triIn+SPACEDIM,quadOut);
 146     double tmp[SPACEDIM];
 147     std::transform(triIn,triIn+SPACEDIM,triIn+SPACEDIM,tmp,std::plus<double>());
 148     //2nd point
 149     std::transform(tmp,tmp+SPACEDIM,quadOut+SPACEDIM,std::bind2nd(std::multiplies<double>(),0.5));
 150     std::transform(tmp,tmp+SPACEDIM,triIn+2*SPACEDIM,tmp,std::plus<double>());
 151     //3rd point
 152     std::transform(tmp,tmp+SPACEDIM,quadOut+2*SPACEDIM,std::bind2nd(std::multiplies<double>(),1/3.));
 153     //4th point
 154     std::transform(triIn,triIn+SPACEDIM,triIn+2*SPACEDIM,tmp,std::plus<double>());
 155     std::transform(tmp,tmp+SPACEDIM,quadOut+3*SPACEDIM,std::bind2nd(std::multiplies<double>(),0.5));
 156   }
 157
 158   /*!
 159    * This method builds a potentially non-convex polygon cell built with the first point of 'triIn' the barycenter of two edges starting or ending with
 160    * the first point of 'triIn' and the barycenter of 'triIn'.
 161    *
 162    * @param triIn is a 6 doubles array in full interlace mode, that represents a triangle.
 163    * @param quadOut is a 8 doubles array filled after the following call.
 164    */
 165   template<int SPACEDIM>
 166   inline void fillDualCellOfPolyg(const double *polygIn, int nPtsPolygonIn, double *polygOut)
 167   {
 168     //1st point
 169     std::copy(polygIn,polygIn+SPACEDIM,polygOut);
 170     std::transform(polygIn,polygIn+SPACEDIM,polygIn+SPACEDIM,polygOut+SPACEDIM,std::plus<double>());
 171     //2nd point
 172     std::transform(polygOut+SPACEDIM,polygOut+2*SPACEDIM,polygOut+SPACEDIM,std::bind2nd(std::multiplies<double>(),0.5));
 173     double tmp[SPACEDIM];
 174     //
 175     for(int i=0;i<nPtsPolygonIn-2;i++)
 176       {
 177         std::transform(polygIn,polygIn+SPACEDIM,polygIn+(i+2)*SPACEDIM,tmp,std::plus<double>());
 178         std::transform(tmp,tmp+SPACEDIM,polygOut+(2*i+3)*SPACEDIM,std::bind2nd(std::multiplies<double>(),0.5));
 179         std::transform(polygIn+(i+1)*SPACEDIM,polygIn+(i+2)*SPACEDIM,tmp,tmp,std::plus<double>());
 180         std::transform(tmp,tmp+SPACEDIM,polygOut+(2*i+2)*SPACEDIM,std::bind2nd(std::multiplies<double>(),1./3.));
 181       }
 182   }
 183
 184   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
 185   /*     calcul les coordonnees du barycentre d'un polygone   */
 186   /*     le vecteur en entree est constitue des coordonnees   */
 187   /*     des sommets du polygone                              */
 188   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
 189
 190   inline std::vector<double> bary_poly(const std::vector<double>& V)
 191   {
 192     std::vector<double> Bary;
 193     long taille=V.size();
 194     double x=0;
 195     double y=0;
 196
 197     for(long i=0;i<taille/2;i++)
 198       {
 199         x=x+V[2*i];
 200         y=y+V[2*i+1];
 201       }
 202     double A=2*x/((double)taille);
 203     double B=2*y/((double)taille);
 204     Bary.push_back(A);//taille vecteur=2*nb de points.
 205     Bary.push_back(B);
 206
 207
 208     return Bary;
 209   }
 210
 211   /*!
 212    * Given 6 coeffs of a Tria6 returns the corresponding value of a given pos
 213    */
 214   inline double computeTria6RefBase(const double *coeffs, const double *pos)
 215   {
 216     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];
 217   }
 218
 219   /*!
 220    * Given xsi,eta in refCoo (length==2) return 6 coeffs in weightedPos.
 221    */
 222   inline void computeWeightedCoeffsInTria6FromRefBase(const double *refCoo, double *weightedPos)
 223   {
 224     weightedPos[0]=(1.-refCoo[0]-refCoo[1])*(1.-2*refCoo[0]-2.*refCoo[1]);
 225     weightedPos[1]=refCoo[0]*(2.*refCoo[0]-1.);
 226     weightedPos[2]=refCoo[1]*(2.*refCoo[1]-1.);
 227     weightedPos[3]=4.*refCoo[0]*(1.-refCoo[0]-refCoo[1]);
 228     weightedPos[4]=4.*refCoo[0]*refCoo[1];
 229     weightedPos[5]=4.*refCoo[1]*(1.-refCoo[0]-refCoo[1]);
 230   }
 231
 232   /*!
 233    * Given 10 coeffs of a Tetra10 returns the corresponding value of a given pos
 234    */
 235   inline double computeTetra10RefBase(const double *coeffs, const double *pos)
 236   {
 237     return coeffs[0]+coeffs[1]*pos[0]+coeffs[2]*pos[1]+coeffs[3]*pos[2]+
 238       coeffs[4]*pos[0]*pos[0]+coeffs[5]*pos[0]*pos[1]+coeffs[6]*pos[0]*pos[2]+
 239       coeffs[7]*pos[1]*pos[1]+coeffs[8]*pos[1]*pos[2]+coeffs[9]*pos[2]*pos[2];
 240   }
 241
 242   /*!
 243    * Given xsi,eta,z in refCoo (length==3) return 10 coeffs in weightedPos.
 244    */
 245   inline void computeWeightedCoeffsInTetra10FromRefBase(const double *refCoo, double *weightedPos)
 246   {
 247     //http://www.cadfamily.com/download/CAE/ABAQUS/The%20Finite%20Element%20Method%20-%20A%20practical%20course%20abaqus.pdf page 217
 248     //L1=1-refCoo[0]-refCoo[1]-refCoo[2]
 249     //L2=refCoo[0] L3=refCoo[1] L4=refCoo[2]
 250     weightedPos[0]=(-2.*(refCoo[0]+refCoo[1]+refCoo[2])+1)*(1-refCoo[0]-refCoo[1]-refCoo[2]);//(2*L1-1)*L1
 251     weightedPos[1]=(2.*refCoo[0]-1.)*refCoo[0];//(2*L2-1)*L2
 252     weightedPos[2]=(2.*refCoo[1]-1.)*refCoo[1];//(2*L3-1)*L3
 253     weightedPos[3]=(2.*refCoo[2]-1.)*refCoo[2];//(2*L4-1)*L4
 254     weightedPos[4]=4.*(1-refCoo[0]-refCoo[1]-refCoo[2])*refCoo[0];//4*L1*L2
 255     weightedPos[5]=4.*refCoo[0]*refCoo[1];//4*L2*L3
 256     weightedPos[6]=4.*(1-refCoo[0]-refCoo[1]-refCoo[2])*refCoo[1];//4*L1*L3
 257     weightedPos[7]=4.*(1-refCoo[0]-refCoo[1]-refCoo[2])*refCoo[2];//4*L1*L4
 258     weightedPos[8]=4.*refCoo[0]*refCoo[2];//4*L2*L4
 259     weightedPos[9]=4.*refCoo[1]*refCoo[2];//4*L3*L4
 260   }
 261
 262   /*!
 263    * \brief Solve system equation in matrix form using Gaussian elimination algorithm
 264    *  \param M - N x N+1 matrix
 265    *  \param sol - vector of N solutions
 266    *  \retval bool - true if succeeded
 267    */
 268   template<unsigned nbRow>
 269   bool solveSystemOfEquations(double M[nbRow][nbRow+1], double* sol)
 270   {
 271     const int nbCol=nbRow+1;
 272
 273     // make upper triangular matrix (forward elimination)
 274
 275     int iR[nbRow];// = { 0, 1, 2 };
 276     for ( int i = 0; i < (int) nbRow; ++i )
 277       iR[i] = i;
 278     for ( int i = 0; i < (int)(nbRow-1); ++i ) // nullify nbRow-1 rows
 279       {
 280         // swap rows to have max value of i-th column in i-th row
 281         double max = std::fabs( M[ iR[i] ][i] );
 282         for ( int r = i+1; r < (int)nbRow; ++r )
 283           {
 284             double m = std::fabs( M[ iR[r] ][i] );
 285             if ( m > max )
 286               {
 287                 max = m;
 288                 std::swap( iR[r], iR[i] );
 289               }
 290           }
 291         if ( max < std::numeric_limits<double>::min() )
 292           {
 293             //sol[0]=1; sol[1]=sol[2]=sol[3]=0;
 294             return false; // no solution
 295           }
 296         // make 0 below M[i][i] (actually we do not modify i-th column)
 297         double* tUpRow = M[ iR[i] ];
 298         for ( int r = i+1; r < (int)nbRow; ++r )
 299           {
 300             double* mRow = M[ iR[r] ];
 301             double coef = mRow[ i ] / tUpRow[ i ];
 302             for ( int c = i+1; c < nbCol; ++c )
 303               mRow[ c ] -= tUpRow[ c ] * coef;
 304           }
 305       }
 306     double* mRow = M[ iR[nbRow-1] ];
 307     if ( std::fabs( mRow[ nbRow-1 ] ) < std::numeric_limits<double>::min() )
 308       {
 309         //sol[0]=1; sol[1]=sol[2]=sol[3]=0;
 310         return false; // no solution
 311       }
 312     mRow[ nbRow ] /= mRow[ nbRow-1 ];
 313
 314     // calculate solution (back substitution)
 315
 316     sol[ nbRow-1 ] = mRow[ nbRow ];
 317
 318     for ( int i = nbRow-2; i+1; --i )
 319       {
 320         mRow = M[ iR[i] ];
 321         sol[ i ] = mRow[ nbRow ];
 322         for ( int j = nbRow-1; j > i; --j )
 323           sol[ i ] -= sol[j]*mRow[ j ];
 324         sol[ i ] /= mRow[ i ];
 325       }
 326
 327     return true;
 328   }
 329
 330
 331   /*!
 332    * \brief Solve system equation in matrix form using Gaussian elimination algorithm
 333    *  \param M - N x N+NB_OF_VARS matrix
 334    *  \param sol - vector of N solutions
 335    *  \retval bool - true if succeeded
 336    */
 337   template<unsigned SZ, unsigned NB_OF_RES>
 338   bool solveSystemOfEquations2(const double *matrix, double *solutions, double eps)
 339   {
 340     unsigned k,j;
 341     int nr,n,m,np;
 342     double s,g;
 343     int mb;
 344     //
 345     double B[SZ*(SZ+NB_OF_RES)];
 346     std::copy(matrix,matrix+SZ*(SZ+NB_OF_RES),B);
 347     //
 348     nr=SZ+NB_OF_RES;
 349     for(k=0;k<SZ;k++)
 350       {
 351         np=nr*k+k;
 352         if(fabs(B[np])<eps)
 353           {
 354             n=k;
 355             do
 356               {
 357                 n++;
 358                 if(fabs(B[nr*k+n])>eps)
 359                   {/* Rows permutation */
 360                     for(m=0;m<nr;m++)
 361                       std::swap(B[nr*k+m],B[nr*n+m]);
 362                   }
 363               }
 364             while (n<(int)SZ);
 365           }
 366         s=B[np];//s is the Pivot
 367         std::transform(B+k*nr,B+(k+1)*nr,B+k*nr,std::bind2nd(std::divides<double>(),s));
 368         for(j=0;j<SZ;j++)
 369           {
 370             if(j!=k)
 371               {
 372                 g=B[j*nr+k];
 373                 for(mb=k;mb<nr;mb++)
 374                   B[j*nr+mb]-=B[k*nr+mb]*g;
 375               }
 376           }
 377       }
 378     for(j=0;j<NB_OF_RES;j++)
 379       for(k=0;k<SZ;k++)
 380         solutions[j*SZ+k]=B[nr*k+SZ+j];
 381     //
 382     return true;
 383   }
 384
 385   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
 386   /*     Calculate barycentric coordinates of a 2D point p */
 387   /*     with respect to the triangle verices.             */
 388   /*     triaCoords are in full interlace                  */
 389   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
 390
 391   template<int SPACEDIM>
 392   inline void barycentric_coords(const double* triaCoords, const double* p, double* bc)
 393   {
 394     // matrix 2x2
 395     double
 396       T11 = triaCoords[0]-triaCoords[2*SPACEDIM], T12 = triaCoords[SPACEDIM]-triaCoords[2*SPACEDIM],
 397       T21 = triaCoords[1]-triaCoords[2*SPACEDIM+1], T22 = triaCoords[SPACEDIM+1]-triaCoords[2*SPACEDIM+1];
 398     // matrix determinant
 399     double Tdet = T11*T22 - T12*T21;
 400     if ( fabs( Tdet ) < std::numeric_limits<double>::min() ) {
 401       bc[0]=1; bc[1]=0; bc[2]=0;
 402       return;
 403     }
 404     // matrix inverse
 405     double t11 = T22, t12 = -T12, t21 = -T21, t22 = T11;
 406     // vector
 407     double r11 = p[0]-triaCoords[2*SPACEDIM], r12 = p[1]-triaCoords[2*SPACEDIM+1];
 408     // barycentric coordinates: mutiply matrix by vector
 409     bc[0] = (t11 * r11 + t12 * r12)/Tdet;
 410     bc[1] = (t21 * r11 + t22 * r12)/Tdet;
 411     bc[2] = 1. - bc[0] - bc[1];
 412   }
 413
 414   /*!
 415    * Calculate barycentric coordinates of a point p with respect to triangle or tetra verices.
 416    * This method makes 2 assumptions :
 417    *    - this is a simplex
 418    *    - spacedim == meshdim. For TRI3 and TRI6 spaceDim is expected to be equal to 2 and for TETRA4 spaceDim is expected to be equal to 3.
 419    *      If not the case (3D surf for example) a previous projection should be done before.
 420    */
 421   inline void barycentric_coords(const std::vector<const double*>& n, const double *p, double *bc)
 422   {
 423     enum { _X, _Y, _Z };
 424     switch(n.size())
 425       {
 426       case 2:
 427         {// SEG 2
 428           double delta=n[0][0]-n[1][0];
 429           bc[0]=fabs((*p-n[1][0])/delta);
 430           bc[1]=fabs((*p-n[0][0])/delta);
 431           break;
 432         }
 433       case 3:
 434         { // TRIA3
 435           // matrix 2x2
 436           double
 437             T11 = n[0][_X]-n[2][_X], T12 = n[1][_X]-n[2][_X],
 438             T21 = n[0][_Y]-n[2][_Y], T22 = n[1][_Y]-n[2][_Y];
 439           // matrix determinant
 440           double Tdet = T11*T22 - T12*T21;
 441           if ( (std::fabs( Tdet) ) < (std::numeric_limits<double>::min()) )
 442             {
 443               bc[0]=1; bc[1]=bc[2]=0; // no solution
 444               return;
 445             }
 446           // matrix inverse
 447           double t11 = T22, t12 = -T12, t21 = -T21, t22 = T11;
 448           // vector
 449           double r11 = p[_X]-n[2][_X], r12 = p[_Y]-n[2][_Y];
 450           // barycentric coordinates: mutiply matrix by vector
 451           bc[0] = (t11 * r11 + t12 * r12)/Tdet;
 452           bc[1] = (t21 * r11 + t22 * r12)/Tdet;
 453           bc[2] = 1. - bc[0] - bc[1];
 454           break;
 455         }
 456       case 4:
 457         { // TETRA4
 458           // Find bc by solving system of 3 equations using Gaussian elimination algorithm
 459           // bc1*( x1 - x4 ) + bc2*( x2 - x4 ) + bc3*( x3 - x4 ) = px - x4
 460           // bc1*( y1 - y4 ) + bc2*( y2 - y4 ) + bc3*( y3 - y4 ) = px - y4
 461           // bc1*( z1 - z4 ) + bc2*( z2 - z4 ) + bc3*( z3 - z4 ) = px - z4
 462
 463           double T[3][4]=
 464             {{ n[0][_X]-n[3][_X], n[1][_X]-n[3][_X], n[2][_X]-n[3][_X], p[_X]-n[3][_X] },
 465              { n[0][_Y]-n[3][_Y], n[1][_Y]-n[3][_Y], n[2][_Y]-n[3][_Y], p[_Y]-n[3][_Y] },
 466              { n[0][_Z]-n[3][_Z], n[1][_Z]-n[3][_Z], n[2][_Z]-n[3][_Z], p[_Z]-n[3][_Z] }};
 467
 468           if ( !solveSystemOfEquations<3>( T, bc ) )
 469             bc[0]=1., bc[1] = bc[2] = bc[3] = 0;
 470           else
 471             bc[ 3 ] = 1. - bc[0] - bc[1] - bc[2];
 472           break;
 473         }
 474       case 6:
 475         {
 476           // TRIA6
 477           double matrix2[48]={1., 0., 0., 0., 0., 0., 0., 0.,
 478                               1., 0., 0., 0., 0., 0., 1., 0.,
 479                               1., 0., 0., 0., 0., 0., 0., 1.,
 480                               1., 0., 0., 0., 0., 0., 0.5, 0.,
 481                               1., 0., 0., 0., 0., 0., 0.5, 0.5,
 482                               1., 0., 0., 0., 0., 0., 0.,0.5};
 483           for(int i=0;i<6;i++)
 484             {
 485               matrix2[8*i+1]=n[i][0];
 486               matrix2[8*i+2]=n[i][1];
 487               matrix2[8*i+3]=n[i][0]*n[i][0];
 488               matrix2[8*i+4]=n[i][0]*n[i][1];
 489               matrix2[8*i+5]=n[i][1]*n[i][1];
 490             }
 491           double res[12];
 492           solveSystemOfEquations2<6,2>(matrix2,res,std::numeric_limits<double>::min());
 493           double refCoo[2];
 494           refCoo[0]=computeTria6RefBase(res,p);
 495           refCoo[1]=computeTria6RefBase(res+6,p);
 496           computeWeightedCoeffsInTria6FromRefBase(refCoo,bc);
 497           break;
 498         }
 499       case 10:
 500         {//TETRA10
 501           double matrix2[130]={1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
 502                                1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.,
 503                                1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.,
 504                                1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.,
 505                                1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.5, 0., 0.,
 506                                1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.5, 0.5, 0.,
 507                                1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,0.5, 0.,
 508                                1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.5,
 509                                1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.5, 0., 0.5,
 510                                1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.5, 0.5};
 511           for(int i=0;i<10;i++)
 512             {
 513               matrix2[13*i+1]=n[i][0];
 514               matrix2[13*i+2]=n[i][1];
 515               matrix2[13*i+3]=n[i][2];
 516               matrix2[13*i+4]=n[i][0]*n[i][0];
 517               matrix2[13*i+5]=n[i][0]*n[i][1];
 518               matrix2[13*i+6]=n[i][0]*n[i][2];
 519               matrix2[13*i+7]=n[i][1]*n[i][1];
 520               matrix2[13*i+8]=n[i][1]*n[i][2];
 521               matrix2[13*i+9]=n[i][2]*n[i][2];
 522             }
 523           double res[30];
 524           solveSystemOfEquations2<10,3>(matrix2,res,std::numeric_limits<double>::min());
 525           double refCoo[3];
 526           refCoo[0]=computeTetra10RefBase(res,p);
 527           refCoo[1]=computeTetra10RefBase(res+10,p);
 528           refCoo[2]=computeTetra10RefBase(res+20,p);
 529           computeWeightedCoeffsInTetra10FromRefBase(refCoo,bc);
 530           break;
 531         }
 532       default:
 533         throw INTERP_KERNEL::Exception("INTERP_KERNEL::barycentric_coords : unrecognized simplex !");
 534       }
 535   }
 536
 537   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
 538   /*         calcul la surface d'un polygone.                 */
 539   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
 540
 541   inline double  Surf_Poly(const std::vector<double>& Poly)
 542   {
 543
 544     double Surface=0;
 545     for(unsigned long i=0; i<(Poly.size())/2-2; i++)
 546       {
 547         double Surf=Surf_Tri( &Poly[0],&Poly[2*(i+1)],&Poly[2*(i+2)] );
 548         Surface=Surface + Surf ;
 549       }
 550     return Surface ;
 551   }
 552
 553   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
 554   /*   fonction qui teste si un point est dans une maille     */
 555   /*   point: P_0                                             */
 556   /*   P_1, P_2, P_3 sommet des mailles                       */
 557   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
 558
 559   inline bool point_dans_triangle(const double* P_0,const double* P_1,
 560                                   const double* P_2,const double* P_3,
 561                                   double eps)
 562   {
 563
 564     bool A=false;
 565     double det_1=mon_determinant(P_1,P_3,P_0);
 566     double det_2=mon_determinant(P_3,P_2,P_0);
 567     double det_3=mon_determinant(P_2,P_1,P_0);
 568     if( (det_1>=-eps && det_2>=-eps && det_3>=-eps) || (det_1<=eps && det_2<=eps && det_3<=eps) )
 569       {
 570         A=true;
 571       }
 572
 573     return A;
 574   }
 575
 576   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
 577   /*fonction pour verifier qu'un point n'a pas deja ete considerer dans   */
 578   /*      le vecteur et le rajouter au vecteur sinon.                     */
 579   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
 580
 581   inline void verif_point_dans_vect(const double* P, std::vector<double>& V, double absolute_precision )
 582   {
 583     long taille=V.size();
 584     bool isPresent=false;
 585     for(long i=0;i<taille/2;i++)
 586       {
 587         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)
 588           isPresent=true;
 589
 590       }
 591     if(!isPresent)
 592       {
 593
 594         V.push_back(P[0]);
 595         V.push_back(P[1]);
 596       }
 597   }
 598
 599   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
 600   /* fonction qui rajoute les sommet du triangle P dans le vecteur V        */
 601   /* si ceux-ci sont compris dans le triangle S et ne sont pas deja dans    */
 602   /* V.                                                                     */
 603   /*sommets de P: P_1, P_2, P_3                                             */
 604   /*sommets de S: P_4, P_5, P_6                                             */
 605   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
 606
 607   inline void rajou_sommet_triangl(const double* P_1,const double* P_2,const double* P_3,
 608                                    const double* P_4,const double* P_5,const double* P_6,
 609                                    std::vector<double>& V, double dim_caracteristic, double precision)
 610   {
 611
 612     double absolute_precision = precision*dim_caracteristic;
 613     bool A_1=INTERP_KERNEL::point_dans_triangle(P_1,P_4,P_5,P_6,absolute_precision);
 614     if(A_1)
 615       verif_point_dans_vect(P_1,V,absolute_precision);
 616     bool A_2=INTERP_KERNEL::point_dans_triangle(P_2,P_4,P_5,P_6,absolute_precision);
 617     if(A_2)
 618       verif_point_dans_vect(P_2,V,absolute_precision);
 619     bool A_3=INTERP_KERNEL::point_dans_triangle(P_3,P_4,P_5,P_6,absolute_precision);
 620     if(A_3)
 621       verif_point_dans_vect(P_3,V,absolute_precision);
 622   }
 623
 624
 625   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _  _ _ _ _ _ _ _ _*/
 626   /*  calcul de l'intersection de deux segments: segments P1P2 avec P3P4      */
 627   /*  . Si l'intersection est non nulle et si celle-ci n'est                  */
 628   /*  n'est pas deja contenue dans Vect on la rajoute a Vect                  */
 629   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _  _ _ _ _ _ _ _ _*/
 630
 631   inline void inters_de_segment(const double * P_1,const double * P_2,
 632                                 const double * P_3,const double * P_4,
 633                                 std::vector<double>& Vect,
 634                                 double dim_caracteristic, double precision)
 635   {
 636     // calcul du determinant de P_1P_2 et P_3P_4.
 637     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]);
 638
 639     double absolute_precision = dim_caracteristic*precision;
 640     if(fabs(det)>absolute_precision)
 641       {
 642         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;
 643
 644         if (k_1 >= -absolute_precision && k_1 <= 1+absolute_precision)
 645           //if( k_1 >= -precision && k_1 <= 1+precision)
 646           {
 647             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;
 648
 649             if (k_2 >= -absolute_precision && k_2 <= 1+absolute_precision)
 650               //if( k_2 >= -precision && k_2 <= 1+precision)
 651               {
 652                 double P_0[2];
 653                 P_0[0]=P_1[0]+k_1*(P_2[0]-P_1[0]);
 654                 P_0[1]=P_1[1]+k_1*(P_2[1]-P_1[1]);
 655                 verif_point_dans_vect(P_0,Vect,absolute_precision);
 656               }
 657           }
 658       }
 659   }
 660
 661
 662
 663   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
 664   /*      calcul l'intersection de deux triangles            */
 665   /* P_1, P_2, P_3: sommets du premier triangle              */
 666   /* P_4, P_5, P_6: sommets du deuxi�me triangle             */
 667   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
 668
 669   inline void intersec_de_triangle(const double* P_1,const double* P_2, const double* P_3,
 670                                    const double* P_4,const double* P_5,const double* P_6,
 671                                    std::vector<double>& Vect, double dim_caracteristic, double precision)
 672   {
 673     inters_de_segment(P_1,P_2,P_4,P_5,Vect, dim_caracteristic, precision);
 674     inters_de_segment(P_1,P_2,P_5,P_6,Vect, dim_caracteristic, precision);
 675     inters_de_segment(P_1,P_2,P_6,P_4,Vect, dim_caracteristic, precision);
 676     inters_de_segment(P_2,P_3,P_4,P_5,Vect, dim_caracteristic, precision);
 677     inters_de_segment(P_2,P_3,P_5,P_6,Vect, dim_caracteristic, precision);
 678     inters_de_segment(P_2,P_3,P_6,P_4,Vect, dim_caracteristic, precision);
 679     inters_de_segment(P_3,P_1,P_4,P_5,Vect, dim_caracteristic, precision);
 680     inters_de_segment(P_3,P_1,P_5,P_6,Vect, dim_caracteristic, precision);
 681     inters_de_segment(P_3,P_1,P_6,P_4,Vect, dim_caracteristic, precision);
 682     rajou_sommet_triangl(P_1,P_2,P_3,P_4,P_5,P_6,Vect, dim_caracteristic, precision);
 683     rajou_sommet_triangl(P_4,P_5,P_6,P_1,P_2,P_3,Vect, dim_caracteristic, precision);
 684   }
 685
 686   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
 687   /* fonction pour verifier qu'un node maille n'a pas deja ete considerer  */
 688   /*  dans le vecteur et le rajouter au vecteur sinon.                     */
 689   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
 690
 691   inline void verif_maill_dans_vect(int Num, std::vector<int>& V)
 692   {
 693     long taille=V.size();
 694     int A=0;
 695     for(long i=0;i<taille;i++)
 696       {
 697         if(Num==V[i])
 698           {
 699             A=1;
 700             break;
 701           }
 702       }
 703     if(A==0)
 704       {V.push_back(Num); }
 705   }
 706
 707   /*! Function that compares two angles from the values of the pairs (sin,cos)*/
 708   /*! Angles are considered in [0, 2Pi] bt are not computed explicitely */
 709   class AngleLess
 710   {
 711   public:
 712     bool operator()(std::pair<double,double>theta1, std::pair<double,double> theta2)
 713     {
 714       double norm1 = sqrt(theta1.first*theta1.first +theta1.second*theta1.second);
 715       double norm2 = sqrt(theta2.first*theta2.first +theta2.second*theta2.second);
 716
 717       double epsilon = 1.e-12;
 718
 719       if( norm1 < epsilon || norm2 < epsilon  )
 720         std::cout << "Warning InterpolationUtils.hxx: AngleLess : Vector with zero norm, cannot define the angle !!!! " << std::endl;
 721
 722       return theta1.second*(norm2 + theta2.first) < theta2.second*(norm1 + theta1.first);
 723
 724     }
 725   };
 726
 727
 728   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
 729   /* fonction pour reconstituer un polygone convexe a partir  */
 730   /*              d'un nuage de point.                        */
 731   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
 732
 733   inline std::vector<double> reconstruct_polygon(const std::vector<double>& V)
 734   {
 735
 736     std::size_t taille=V.size();
 737
 738     //VB : why 6 ?
 739
 740     if(taille<=6)
 741       {return V;}
 742     else
 743       {
 744         double *COS=new double[taille/2];
 745         double *SIN=new double[taille/2];
 746         //double *angle=new double[taille/2];
 747         std::vector<double> Bary=bary_poly(V);
 748         COS[0]=1.0;
 749         SIN[0]=0.0;
 750         //angle[0]=0.0;
 751         for(std::size_t i=0; i<taille/2-1;i++)
 752           {
 753             std::vector<double> Trigo=calcul_cos_et_sin(&Bary[0],&V[0],&V[2*(i+1)]);
 754             COS[i+1]=Trigo[0];
 755             SIN[i+1]=Trigo[1];
 756             //if(SIN[i+1]>=0)
 757             //    {angle[i+1]=atan2(SIN[i+1],COS[i+1]);}
 758             //             else
 759             //               {angle[i+1]=-atan2(SIN[i+1],COS[i+1]);}
 760           }
 761
 762         //ensuite on ordonne les angles.
 763         std::vector<double> Pt_ordonne;
 764         Pt_ordonne.reserve(taille);
 765         //        std::multimap<double,int> Ordre;
 766         std::multimap<std::pair<double,double>,int, AngleLess> CosSin;
 767         for(std::size_t i=0;i<taille/2;i++)
 768           {
 769             //  Ordre.insert(std::make_pair(angle[i],i));
 770             CosSin.insert(std::make_pair(std::make_pair(SIN[i],COS[i]),i));
 771           }
 772         //        std::multimap <double,int>::iterator mi;
 773         std::multimap<std::pair<double,double>,int, AngleLess>::iterator   micossin;
 774         //         for(mi=Ordre.begin();mi!=Ordre.end();mi++)
 775         //           {
 776         //             int j=(*mi).second;
 777         //             Pt_ordonne.push_back(V[2*j]);
 778         //             Pt_ordonne.push_back(V[2*j+1]);
 779         //           }
 780         for(micossin=CosSin.begin();micossin!=CosSin.end();micossin++)
 781           {
 782             int j=(*micossin).second;
 783             Pt_ordonne.push_back(V[2*j]);
 784             Pt_ordonne.push_back(V[2*j+1]);
 785           }
 786         delete [] COS;
 787         delete [] SIN;
 788         //        delete [] angle;
 789         return Pt_ordonne;
 790       }
 791   }
 792
 793   template<int DIM, NumberingPolicy numPol, class MyMeshType>
 794   inline void getElemBB(double* bb, const double *coordsOfMesh, int iP, int nb_nodes)
 795   {
 796     bb[0]=std::numeric_limits<double>::max();
 797     bb[1]=-std::numeric_limits<double>::max();
 798     bb[2]=std::numeric_limits<double>::max();
 799     bb[3]=-std::numeric_limits<double>::max();
 800     bb[4]=std::numeric_limits<double>::max();
 801     bb[5]=-std::numeric_limits<double>::max();
 802
 803     for (int i=0; i<nb_nodes; i++)
 804       {
 805         double x = coordsOfMesh[3*(iP+i)];
 806         double y = coordsOfMesh[3*(iP+i)+1];
 807         double z = coordsOfMesh[3*(iP+i)+2];
 808         bb[0]=(x<bb[0])?x:bb[0];
 809         bb[1]=(x>bb[1])?x:bb[1];
 810         bb[2]=(y<bb[2])?y:bb[2];
 811         bb[3]=(y>bb[3])?y:bb[3];
 812         bb[4]=(z<bb[4])?z:bb[4];
 813         bb[5]=(z>bb[5])?z:bb[5];
 814       }
 815   }
 816
 817   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
 818   /* Computes the dot product of a and b */
 819   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
 820   template<int dim>
 821   inline double dotprod( const double * a, const double * b)
 822   {
 823     double result=0;
 824     for(int idim = 0; idim < dim ; idim++) result += a[idim]*b[idim];
 825     return result;
 826   }
 827   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
 828   /* Computes the norm of vector v */
 829   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
 830   template<int dim>
 831   inline double norm(const double * v)
 832   {
 833     double result =0;
 834     for(int idim =0; idim<dim; idim++) result+=v[idim]*v[idim];
 835     return sqrt(result);
 836   }
 837   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
 838   /* Computes the square norm of vector a-b */
 839   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
 840   template<int dim>
 841   inline double distance2( const double * a, const double * b)
 842   {
 843     double result =0;
 844     for(int idim =0; idim<dim; idim++) result+=(a[idim]-b[idim])*(a[idim]-b[idim]);
 845     return result;
 846   }
 847   template<class T, int dim>
 848   inline double distance2(  T * a, int inda, T * b, int indb)
 849   {
 850     double result =0;
 851     for(int idim =0; idim<dim; idim++) result += ((*a)[inda+idim] - (*b)[indb+idim])* ((*a)[inda+idim] - (*b)[indb+idim]);
 852     return result;
 853   }
 854   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
 855   /* Computes the determinant of a and b */
 856   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
 857   inline double determinant ( double *  a, double * b)
 858   {
 859     return a[0]*b[1]-a[1]*b[0];
 860   }
 861   inline double determinant ( double *  a, double * b, double * c)
 862   {
 863     return a[0]*determinant(b+1,c+1)-b[0]*determinant(a+1,c+1)+c[0]*determinant(a+1,b+1);
 864   }
 865
 866   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
 867   /* Computes the cross product of AB and AC */
 868   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
 869
 870   template<int dim> inline void crossprod( const double * A, const double * B, const double * C, double * V);
 871
 872   template<> inline
 873   void crossprod<2>( const double * A, const double * B, const double * C, double * V)
 874   {
 875     double AB[2];
 876     double AC[2];
 877     for(int idim =0; idim<2; idim++) AB[idim] = B[idim]-A[idim];//B-A
 878     for(int idim =0; idim<2; idim++) AC[idim] = C[idim]-A[idim];//C-A;
 879
 880     V[0]=determinant(AB,AC);
 881     V[1]=0;
 882   }
 883   template<> inline
 884   void crossprod<3>( const double * A, const double * B, const double * C, double * V)
 885   {
 886     double AB[3];
 887     double AC[3];
 888     for(int idim =0; idim<3; idim++) AB[idim] = B[idim]-A[idim];//B-A
 889     for(int idim =0; idim<3; idim++) AC[idim] = C[idim]-A[idim];//C-A;
 890
 891     V[0]=AB[1]*AC[2]-AB[2]*AC[1];
 892     V[1]=-AB[0]*AC[2]+AB[2]*AC[0];
 893     V[2]=AB[0]*AC[1]-AB[1]*AC[0];
 894   }
 895   template<> inline
 896   void crossprod<1>( const double * /*A*/, const double * /*B*/, const double * /*C*/, double * /*V*/)
 897   {
 898     // just to be able to compile
 899   }
 900
 901   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
 902   /* Checks wether point A is inside the quadrangle BCDE */
 903   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
 904
 905   template<int dim> inline double check_inside(const double* A,const double* B,const double* C,const double* D,
 906                                                const double* E,double* ABC, double* ADE)
 907   {
 908     crossprod<dim>(A,B,C,ABC);
 909     crossprod<dim>(A,D,E,ADE);
 910     return dotprod<dim>(ABC,ADE);
 911   }
 912
 913
 914   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
 915   /* Computes the geometric angle (in [0,Pi]) between two non zero vectors AB and AC */
 916   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
 917   template<int dim> inline double angle(const double * A, const double * B, const double * C, double * n)
 918   {
 919     double AB[dim];
 920     double AC[dim];
 921     double orthAB[dim];
 922
 923     for(int idim =0; idim<dim; idim++) AB[idim] = B[idim]-A[idim];//B-A;
 924     for(int idim =0; idim<dim; idim++) AC[idim] = C[idim]-A[idim];//C-A;
 925
 926     double normAB= norm<dim>(AB);
 927     for(int idim =0; idim<dim; idim++) AB[idim]/=normAB;
 928
 929     double normAC= norm<dim>(AC);
 930     double AB_dot_AC=dotprod<dim>(AB,AC);
 931     for(int idim =0; idim<dim; idim++) orthAB[idim] = AC[idim]-AB_dot_AC*AB[idim];
 932
 933     double denom= normAC+AB_dot_AC;
 934     double numer=norm<dim>(orthAB);
 935
 936     return 2*atan2(numer,denom);
 937   }
 938
 939   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
 940   /* Tells whether the frame constituted of vectors AB, AC and n is direct */
 941   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
 942   template<int dim> inline double direct_frame(const double * A, const double * B, const double * C, double * n);
 943   template<> inline
 944   double direct_frame<2>(const double * A, const double * B, const double * C, double * n)
 945   {
 946     double AB[2];
 947     double AC[2];
 948     for(int idim =0; idim<2; idim++) AB[idim] = B[idim]-A[idim];//B-A;
 949     for(int idim =0; idim<2; idim++) AC[idim] = C[idim]-A[idim];//C-A;
 950
 951     return  determinant(AB,AC)*n[0];
 952   }
 953   template<> inline
 954   double direct_frame<3>(const double * A, const double * B, const double * C, double * n)
 955   {
 956     double AB[3];
 957     double AC[3];
 958     for(int idim =0; idim<3; idim++) AB[idim] = B[idim]-A[idim];//B-A;
 959     for(int idim =0; idim<3; idim++) AC[idim] = C[idim]-A[idim];//C-A;
 960
 961     return determinant(AB,AC,n)>0;
 962   }
 963
 964   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
 965   /*      calcul l'intersection de deux polygones COPLANAIRES */
 966   /* en dimension DIM (2 ou 3). Si DIM=3 l'algorithme ne considere*/
 967   /* que les deux premieres coordonnees de chaque point */
 968   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
 969   template<int DIM> inline void intersec_de_polygone(const double * Coords_A, const double * Coords_B,
 970                                                      int nb_NodesA, int nb_NodesB,
 971                                                      std::vector<double>& inter,
 972                                                      double dim_caracteristic, double precision)
 973   {
 974     for(int i_A = 1; i_A<nb_NodesA-1; i_A++)
 975       {
 976         for(int i_B = 1; i_B<nb_NodesB-1; i_B++)
 977           {
 978             INTERP_KERNEL::intersec_de_triangle(&Coords_A[0],&Coords_A[DIM*i_A],&Coords_A[DIM*(i_A+1)],
 979                                                 &Coords_B[0],&Coords_B[DIM*i_B],&Coords_B[DIM*(i_B+1)],
 980                                                 inter, dim_caracteristic, precision);
 981           }
 982       }
 983     int nb_inter=((int)inter.size())/DIM;
 984     if(nb_inter >3) inter=INTERP_KERNEL::reconstruct_polygon(inter);
 985   }
 986
 987   /*_ _ _ _ _ _ _ _ _
 988    *_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
 989    *  fonctions qui calcule l'aire d'un polygone en dimension 2 ou 3
 990    *_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
 991   template<int DIM> inline double polygon_area(std::vector<double>& inter)
 992   {
 993     double result=0.;
 994     double area[DIM];
 995
 996     for(int i = 1; i<(int)inter.size()/DIM-1; i++)
 997       {
 998         INTERP_KERNEL::crossprod<DIM>(&inter[0],&inter[DIM*i],&inter[DIM*(i+1)],area);
 999         result +=0.5*norm<DIM>(area);
1000       }
1001     return result;
1002   }
1003
1004   template<int DIM> inline double polygon_area(std::deque<double>& inter)
1005   {
1006     double result=0.;
1007     double area[DIM];
1008
1009     for(int i = 1; i<(int)inter.size()/DIM-1; i++)
1010       {
1011         INTERP_KERNEL::crossprod<DIM>(&inter[0],&inter[DIM*i],&inter[DIM*(i+1)],area);
1012         result +=0.5*norm<DIM>(area);
1013       }
1014     return result;
1015   }
1016
1017   /*! Computes the triple product (XA^XB).XC (in 3D)*/
1018   inline double triple_product(const double* A, const double*B, const double*C, const double*X)
1019   {
1020     double XA[3];
1021     XA[0]=A[0]-X[0];
1022     XA[1]=A[1]-X[1];
1023     XA[2]=A[2]-X[2];
1024     double XB[3];
1025     XB[0]=B[0]-X[0];
1026     XB[1]=B[1]-X[1];
1027     XB[2]=B[2]-X[2];
1028     double XC[3];
1029     XC[0]=C[0]-X[0];
1030     XC[1]=C[1]-X[1];
1031     XC[2]=C[2]-X[2];
1032
1033     return
1034       (XA[1]*XB[2]-XA[2]*XB[1])*XC[0]+
1035       (XA[2]*XB[0]-XA[0]*XB[2])*XC[1]+
1036       (XA[0]*XB[1]-XA[1]*XB[0])*XC[2];
1037   }
1038
1039   /*! Subroutine of checkEqualPolygins that tests if two list of nodes (not necessarily distincts) describe the same polygon, assuming they share a comon point.*/
1040   /*! Indexes istart1 and istart2 designate two points P1 in L1 and P2 in L2 that have identical coordinates. Generally called with istart1=0.*/
1041   /*! Integer sign ( 1 or -1) indicate the direction used in going all over L2. */
1042   template<class T, int dim>
1043   bool checkEqualPolygonsOneDirection(T * L1, T * L2, int size1, int size2, int istart1, int istart2, double epsilon, int sign)
1044   {
1045     int i1 = istart1;
1046     int i2 = istart2;
1047     int i1next = ( i1 + 1 ) % size1;
1048     int i2next = ( i2 + sign +size2) % size2;
1049
1050     while(true)
1051       {
1052         while( i1next!=istart1 && distance2<T,dim>(L1,i1*dim, L1,i1next*dim) < epsilon ) i1next = (  i1next + 1 ) % size1;
1053         while( i2next!=istart2 && distance2<T,dim>(L2,i2*dim, L2,i2next*dim) < epsilon ) i2next = (  i2next + sign +size2 ) % size2;
1054
1055         if(i1next == istart1)
1056           {
1057             if(i2next == istart2)
1058               return true;
1059             else return false;
1060           }
1061         else
1062           if(i2next == istart2)
1063             return false;
1064           else
1065             {
1066               if(distance2<T,dim>(L1,i1next*dim, L2,i2next*dim) > epsilon )
1067                 return false;
1068               else
1069                 {
1070                   i1 = i1next;
1071                   i2 = i2next;
1072                   i1next = ( i1 + 1 ) % size1;
1073                   i2next = ( i2 + sign + size2 ) % size2;
1074                 }
1075             }
1076       }
1077   }
1078
1079   /*! Tests if two list of nodes (not necessarily distincts) describe the same polygon.*/
1080   /*! Existence of multiple points in the list is considered.*/
1081   template<class T, int dim>
1082   bool checkEqualPolygons(T * L1, T * L2, double epsilon)
1083   {
1084     if(L1==NULL || L2==NULL)
1085       {
1086         std::cout << "Warning InterpolationUtils.hxx:checkEqualPolygonsPointer: Null pointer " << std::endl;
1087         throw(Exception("big error: not closed polygon..."));
1088       }
1089
1090     int size1 = (*L1).size()/dim;
1091     int size2 = (*L2).size()/dim;
1092     int istart1 = 0;
1093     int istart2 = 0;
1094
1095     while( istart2 < size2  && distance2<T,dim>(L1,istart1*dim, L2,istart2*dim) > epsilon ) istart2++;
1096
1097     if(istart2 == size2)
1098       {
1099         return (size1 == 0) && (size2 == 0);
1100       }
1101     else
1102       return   checkEqualPolygonsOneDirection<T,dim>( L1, L2, size1, size2, istart1, istart2, epsilon,  1)
1103         || checkEqualPolygonsOneDirection<T,dim>( L1, L2, size1, size2, istart1, istart2, epsilon, -1);
1104
1105   }
1106 }
1107
1108
1109 #endif