src/INTERP_KERNEL/InterpolationUtils.hxx

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