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 explicitely */
 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   /* Computes the dot product of a and b */
 951   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
 952   template<int dim>
 953   inline double dotprod( const double * a, const double * b)
 954   {
 955     double result=0;
 956     for(int idim = 0; idim < dim ; idim++) result += a[idim]*b[idim];
 957     return result;
 958   }
 959   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
 960   /* Computes the norm of vector v */
 961   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
 962   template<int dim>
 963   inline double norm(const double * v)
 964   {
 965     double result =0;
 966     for(int idim =0; idim<dim; idim++) result+=v[idim]*v[idim];
 967     return sqrt(result);
 968   }
 969   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
 970   /* Computes the square norm of vector a-b */
 971   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
 972   template<int dim>
 973   inline double distance2( const double * a, const double * b)
 974   {
 975     double result =0;
 976     for(int idim =0; idim<dim; idim++) result+=(a[idim]-b[idim])*(a[idim]-b[idim]);
 977     return result;
 978   }
 979   template<class T, int dim>
 980   inline double distance2(  T * a, int inda, T * b, int indb)
 981   {
 982     double result =0;
 983     for(int idim =0; idim<dim; idim++) result += ((*a)[inda+idim] - (*b)[indb+idim])* ((*a)[inda+idim] - (*b)[indb+idim]);
 984     return result;
 985   }
 986   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
 987   /* Computes the determinant of a and b */
 988   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
 989   inline double determinant ( double *  a, double * b)
 990   {
 991     return a[0]*b[1]-a[1]*b[0];
 992   }
 993   inline double determinant ( double *  a, double * b, double * c)
 994   {
 995     return a[0]*determinant(b+1,c+1)-b[0]*determinant(a+1,c+1)+c[0]*determinant(a+1,b+1);
 996   }
 997
 998   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
 999   /* Computes the cross product of AB and AC */
1000   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
1001
1002   template<int dim> inline void crossprod( const double * A, const double * B, const double * C, double * V);
1003
1004   template<> inline
1005   void crossprod<2>( const double * A, const double * B, const double * C, double * V)
1006   {
1007     double AB[2];
1008     double AC[2];
1009     for(int idim =0; idim<2; idim++) AB[idim] = B[idim]-A[idim];//B-A
1010     for(int idim =0; idim<2; idim++) AC[idim] = C[idim]-A[idim];//C-A;
1011
1012     V[0]=determinant(AB,AC);
1013     V[1]=0;
1014   }
1015   template<> inline
1016   void crossprod<3>( const double * A, const double * B, const double * C, double * V)
1017   {
1018     double AB[3];
1019     double AC[3];
1020     for(int idim =0; idim<3; idim++) AB[idim] = B[idim]-A[idim];//B-A
1021     for(int idim =0; idim<3; idim++) AC[idim] = C[idim]-A[idim];//C-A;
1022
1023     V[0]=AB[1]*AC[2]-AB[2]*AC[1];
1024     V[1]=-AB[0]*AC[2]+AB[2]*AC[0];
1025     V[2]=AB[0]*AC[1]-AB[1]*AC[0];
1026   }
1027   template<> inline
1028   void crossprod<1>( const double * /*A*/, const double * /*B*/, const double * /*C*/, double * /*V*/)
1029   {
1030     // just to be able to compile
1031   }
1032
1033   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
1034   /* Checks wether point A is inside the quadrangle BCDE */
1035   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
1036
1037   template<int dim> inline double check_inside(const double* A,const double* B,const double* C,const double* D,
1038                                                const double* E,double* ABC, double* ADE)
1039   {
1040     crossprod<dim>(A,B,C,ABC);
1041     crossprod<dim>(A,D,E,ADE);
1042     return dotprod<dim>(ABC,ADE);
1043   }
1044
1045
1046   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
1047   /* Computes the geometric angle (in [0,Pi]) between two non zero vectors AB and AC */
1048   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
1049   template<int dim> inline double angle(const double * A, const double * B, const double * C, double * n)
1050   {
1051     double AB[dim];
1052     double AC[dim];
1053     double orthAB[dim];
1054
1055     for(int idim =0; idim<dim; idim++) AB[idim] = B[idim]-A[idim];//B-A;
1056     for(int idim =0; idim<dim; idim++) AC[idim] = C[idim]-A[idim];//C-A;
1057
1058     double normAB= norm<dim>(AB);
1059     for(int idim =0; idim<dim; idim++) AB[idim]/=normAB;
1060
1061     double normAC= norm<dim>(AC);
1062     double AB_dot_AC=dotprod<dim>(AB,AC);
1063     for(int idim =0; idim<dim; idim++) orthAB[idim] = AC[idim]-AB_dot_AC*AB[idim];
1064
1065     double denom= normAC+AB_dot_AC;
1066     double numer=norm<dim>(orthAB);
1067
1068     return 2*atan2(numer,denom);
1069   }
1070
1071   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
1072   /* Tells whether the frame constituted of vectors AB, AC and n is direct */
1073   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
1074   template<int dim> inline double direct_frame(const double * A, const double * B, const double * C, double * n);
1075   template<> inline
1076   double direct_frame<2>(const double * A, const double * B, const double * C, double * n)
1077   {
1078     double AB[2];
1079     double AC[2];
1080     for(int idim =0; idim<2; idim++) AB[idim] = B[idim]-A[idim];//B-A;
1081     for(int idim =0; idim<2; idim++) AC[idim] = C[idim]-A[idim];//C-A;
1082
1083     return  determinant(AB,AC)*n[0];
1084   }
1085   template<> inline
1086   double direct_frame<3>(const double * A, const double * B, const double * C, double * n)
1087   {
1088     double AB[3];
1089     double AC[3];
1090     for(int idim =0; idim<3; idim++) AB[idim] = B[idim]-A[idim];//B-A;
1091     for(int idim =0; idim<3; idim++) AC[idim] = C[idim]-A[idim];//C-A;
1092
1093     return determinant(AB,AC,n)>0;
1094   }
1095
1096   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
1097   /*      calcul l'intersection de deux polygones COPLANAIRES */
1098   /* en dimension DIM (2 ou 3). Si DIM=3 l'algorithme ne considere*/
1099   /* que les deux premieres coordonnees de chaque point */
1100   /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
1101   template<int DIM> inline void intersec_de_polygone(const double * Coords_A, const double * Coords_B,
1102                                                      int nb_NodesA, int nb_NodesB,
1103                                                      std::vector<double>& inter,
1104                                                      double dim_caracteristic, double precision)
1105   {
1106     for(int i_A = 1; i_A<nb_NodesA-1; i_A++)
1107       {
1108         for(int i_B = 1; i_B<nb_NodesB-1; i_B++)
1109           {
1110             INTERP_KERNEL::intersec_de_triangle(&Coords_A[0],&Coords_A[DIM*i_A],&Coords_A[DIM*(i_A+1)],
1111                                                 &Coords_B[0],&Coords_B[DIM*i_B],&Coords_B[DIM*(i_B+1)],
1112                                                 inter, dim_caracteristic, precision);
1113           }
1114       }
1115     int nb_inter=((int)inter.size())/DIM;
1116     if(nb_inter >3) inter=INTERP_KERNEL::reconstruct_polygon(inter);
1117   }
1118
1119   /*_ _ _ _ _ _ _ _ _
1120    *_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
1121    *  fonctions qui calcule l'aire d'un polygone en dimension 2 ou 3
1122    *_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
1123   template<int DIM> inline double polygon_area(std::vector<double>& inter)
1124   {
1125     double result=0.;
1126     double area[DIM];
1127
1128     for(int i = 1; i<(int)inter.size()/DIM-1; i++)
1129       {
1130         INTERP_KERNEL::crossprod<DIM>(&inter[0],&inter[DIM*i],&inter[DIM*(i+1)],area);
1131         result +=0.5*norm<DIM>(area);
1132       }
1133     return result;
1134   }
1135
1136   template<int DIM> inline double polygon_area(std::deque<double>& inter)
1137   {
1138     double result=0.;
1139     double area[DIM];
1140
1141     for(int i = 1; i<(int)inter.size()/DIM-1; i++)
1142       {
1143         INTERP_KERNEL::crossprod<DIM>(&inter[0],&inter[DIM*i],&inter[DIM*(i+1)],area);
1144         result +=0.5*norm<DIM>(area);
1145       }
1146     return result;
1147   }
1148
1149   /*! Computes the triple product (XA^XB).XC (in 3D)*/
1150   inline double triple_product(const double* A, const double*B, const double*C, const double*X)
1151   {
1152     double XA[3];
1153     XA[0]=A[0]-X[0];
1154     XA[1]=A[1]-X[1];
1155     XA[2]=A[2]-X[2];
1156     double XB[3];
1157     XB[0]=B[0]-X[0];
1158     XB[1]=B[1]-X[1];
1159     XB[2]=B[2]-X[2];
1160     double XC[3];
1161     XC[0]=C[0]-X[0];
1162     XC[1]=C[1]-X[1];
1163     XC[2]=C[2]-X[2];
1164
1165     return
1166       (XA[1]*XB[2]-XA[2]*XB[1])*XC[0]+
1167       (XA[2]*XB[0]-XA[0]*XB[2])*XC[1]+
1168       (XA[0]*XB[1]-XA[1]*XB[0])*XC[2];
1169   }
1170
1171   /*! Subroutine of checkEqualPolygins that tests if two list of nodes (not necessarily distincts) describe the same polygon, assuming they share a comon point.*/
1172   /*! Indexes istart1 and istart2 designate two points P1 in L1 and P2 in L2 that have identical coordinates. Generally called with istart1=0.*/
1173   /*! Integer sign ( 1 or -1) indicate the direction used in going all over L2. */
1174   template<class T, int dim>
1175   bool checkEqualPolygonsOneDirection(T * L1, T * L2, int size1, int size2, int istart1, int istart2, double epsilon, int sign)
1176   {
1177     int i1 = istart1;
1178     int i2 = istart2;
1179     int i1next = ( i1 + 1 ) % size1;
1180     int i2next = ( i2 + sign +size2) % size2;
1181
1182     while(true)
1183       {
1184         while( i1next!=istart1 && distance2<T,dim>(L1,i1*dim, L1,i1next*dim) < epsilon ) i1next = (  i1next + 1 ) % size1;
1185         while( i2next!=istart2 && distance2<T,dim>(L2,i2*dim, L2,i2next*dim) < epsilon ) i2next = (  i2next + sign +size2 ) % size2;
1186
1187         if(i1next == istart1)
1188           {
1189             if(i2next == istart2)
1190               return true;
1191             else return false;
1192           }
1193         else
1194           if(i2next == istart2)
1195             return false;
1196           else
1197             {
1198               if(distance2<T,dim>(L1,i1next*dim, L2,i2next*dim) > epsilon )
1199                 return false;
1200               else
1201                 {
1202                   i1 = i1next;
1203                   i2 = i2next;
1204                   i1next = ( i1 + 1 ) % size1;
1205                   i2next = ( i2 + sign + size2 ) % size2;
1206                 }
1207             }
1208       }
1209   }
1210
1211   /*! Tests if two list of nodes (not necessarily distincts) describe the same polygon.*/
1212   /*! Existence of multiple points in the list is considered.*/
1213   template<class T, int dim>
1214   bool checkEqualPolygons(T * L1, T * L2, double epsilon)
1215   {
1216     if(L1==NULL || L2==NULL)
1217       {
1218         std::cout << "Warning InterpolationUtils.hxx:checkEqualPolygonsPointer: Null pointer " << std::endl;
1219         throw(Exception("big error: not closed polygon..."));
1220       }
1221
1222     int size1 = (*L1).size()/dim;
1223     int size2 = (*L2).size()/dim;
1224     int istart1 = 0;
1225     int istart2 = 0;
1226
1227     while( istart2 < size2  && distance2<T,dim>(L1,istart1*dim, L2,istart2*dim) > epsilon ) istart2++;
1228
1229     if(istart2 == size2)
1230       {
1231         return (size1 == 0) && (size2 == 0);
1232       }
1233     else
1234       return   checkEqualPolygonsOneDirection<T,dim>( L1, L2, size1, size2, istart1, istart2, epsilon,  1)
1235         || checkEqualPolygonsOneDirection<T,dim>( L1, L2, size1, size2, istart1, istart2, epsilon, -1);
1236
1237   }
1238 }
1239
1240
1241 #endif