src/INTERP_KERNEL/InterpolationUtils.hxx

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