1 // Copyright (C) 2007-2016 CEA/DEN, EDF R&D
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.
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.
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
17 // See http://www.salome-platform.org/ or email : webmaster.salome@opencascade.com
19 // Author : Anthony Geay (CEA/DEN)
21 #ifndef __INTERPOLATIONUTILS_HXX__
22 #define __INTERPOLATIONUTILS_HXX__
24 #include "INTERPKERNELDefines.hxx"
25 #include "InterpKernelException.hxx"
26 #include "VolSurfUser.hxx"
28 #include "NormalizedUnstructuredMesh.hxx"
40 namespace INTERP_KERNEL
42 template<class ConnType, NumberingPolicy numPol>
43 class OTT//OffsetToolTrait
47 template<class ConnType>
48 class OTT<ConnType,ALL_C_MODE>
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; }
57 template<class ConnType>
58 class OTT<ConnType,ALL_FORTRAN_MODE>
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; }
67 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
68 /* calcul la surface d'un triangle */
69 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
71 inline double Surf_Tri(const double* P_1,const double* P_2,const double* P_3)
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);
78 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
79 /* fonction qui calcul le determinant */
80 /* de deux vecteur(cf doc CGAL). */
81 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
83 //fonction qui calcul le determinant des vecteurs: P3P1 et P3P2
86 inline double mon_determinant(const double* P_1,
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]);
94 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
95 //calcul la norme du vecteur P1P2
97 inline double norme_vecteur(const double* P_1,const double* P_2)
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);
105 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
106 /* calcul le cos et le sin de l'angle P1P2,P1P3 */
107 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
109 inline std::vector<double> calcul_cos_et_sin(const double* P_1,
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);
119 double N=P1_P2*P1_P2+P3_P1*P3_P1-P2_P3*P2_P3;
120 double D=2.0*P1_P2*P3_P1;
122 if (COS>1.0) COS=1.0;
123 if (COS<-1.0) COS=-1.0;
125 double V=mon_determinant(P_2,P_3,P_1);
126 double D_1=P1_P2*P3_P1;
128 if (SIN>1.0) SIN=1.0;
129 if (SIN<-1.0) SIN=-1.0;
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'.
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.
143 template<int SPACEDIM>
144 inline void fillDualCellOfTri(const double *triIn, double *quadOut)
147 std::copy(triIn,triIn+SPACEDIM,quadOut);
148 double tmp[SPACEDIM];
149 std::transform(triIn,triIn+SPACEDIM,triIn+SPACEDIM,tmp,std::plus<double>());
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>());
154 std::transform(tmp,tmp+SPACEDIM,quadOut+2*SPACEDIM,std::bind2nd(std::multiplies<double>(),1/3.));
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));
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'.
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.
167 template<int SPACEDIM>
168 inline void fillDualCellOfPolyg(const double *polygIn, int nPtsPolygonIn, double *polygOut)
171 std::copy(polygIn,polygIn+SPACEDIM,polygOut);
172 std::transform(polygIn,polygIn+SPACEDIM,polygIn+SPACEDIM,polygOut+SPACEDIM,std::plus<double>());
174 std::transform(polygOut+SPACEDIM,polygOut+2*SPACEDIM,polygOut+SPACEDIM,std::bind2nd(std::multiplies<double>(),0.5));
175 double tmp[SPACEDIM];
177 for(int i=0;i<nPtsPolygonIn-2;i++)
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.));
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 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
192 inline std::vector<double> bary_poly(const std::vector<double>& V)
194 std::vector<double> Bary;
195 long taille=V.size();
199 for(long i=0;i<taille/2;i++)
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.
214 * Given 6 coeffs of a Tria6 returns the corresponding value of a given pos
216 inline double computeTria6RefBase(const double *coeffs, const double *pos)
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];
222 * Given xsi,eta in refCoo (length==2) return 6 coeffs in weightedPos.
224 inline void computeWeightedCoeffsInTria6FromRefBase(const double *refCoo, double *weightedPos)
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]);
235 * Given 10 coeffs of a Tetra10 returns the corresponding value of a given pos
237 inline double computeTetra10RefBase(const double *coeffs, const double *pos)
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];
245 * Given xsi,eta,z in refCoo (length==3) return 10 coeffs in weightedPos.
247 inline void computeWeightedCoeffsInTetra10FromRefBase(const double *refCoo, double *weightedPos)
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
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
270 template<unsigned nbRow>
271 bool solveSystemOfEquations(double M[nbRow][nbRow+1], double* sol)
273 const int nbCol=nbRow+1;
275 // make upper triangular matrix (forward elimination)
277 int iR[nbRow];// = { 0, 1, 2 };
278 for ( int i = 0; i < (int) nbRow; ++i )
280 for ( int i = 0; i < (int)(nbRow-1); ++i ) // nullify nbRow-1 rows
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 )
286 double m = std::fabs( M[ iR[r] ][i] );
290 std::swap( iR[r], iR[i] );
293 if ( max < std::numeric_limits<double>::min() )
295 //sol[0]=1; sol[1]=sol[2]=sol[3]=0;
296 return false; // no solution
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 )
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;
308 double* mRow = M[ iR[nbRow-1] ];
309 if ( std::fabs( mRow[ nbRow-1 ] ) < std::numeric_limits<double>::min() )
311 //sol[0]=1; sol[1]=sol[2]=sol[3]=0;
312 return false; // no solution
314 mRow[ nbRow ] /= mRow[ nbRow-1 ];
316 // calculate solution (back substitution)
318 sol[ nbRow-1 ] = mRow[ nbRow ];
320 for ( int i = nbRow-2; i+1; --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 ];
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
339 template<unsigned SZ, unsigned NB_OF_RES>
340 bool solveSystemOfEquations2(const double *matrix, double *solutions, double eps)
347 double B[SZ*(SZ+NB_OF_RES)];
348 std::copy(matrix,matrix+SZ*(SZ+NB_OF_RES),B);
360 if(fabs(B[nr*k+n])>eps)
361 {/* Rows permutation */
363 std::swap(B[nr*k+m],B[nr*n+m]);
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));
376 B[j*nr+mb]-=B[k*nr+mb]*g;
380 for(j=0;j<NB_OF_RES;j++)
382 solutions[j*SZ+k]=B[nr*k+SZ+j];
387 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
388 /* Calculate barycentric coordinates of a 2D point p */
389 /* with respect to the triangle verices. */
390 /* triaCoords are in full interlace */
391 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
393 template<int SPACEDIM>
394 inline void barycentric_coords(const double* triaCoords, const double* p, double* bc)
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;
407 double t11 = T22, t12 = -T12, t21 = -T21, t22 = T11;
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];
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.
423 inline void barycentric_coords(const std::vector<const double*>& n, const double *p, double *bc)
425 enum { _XX=0, _YY, _ZZ };
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);
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()) )
445 bc[0]=1; bc[1]=bc[2]=0; // no solution
449 double t11 = T22, t12 = -T12, t21 = -T21, t22 = T11;
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];
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
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] }};
470 if ( !solveSystemOfEquations<3>( T, bc ) )
471 bc[0]=1., bc[1] = bc[2] = bc[3] = 0;
473 bc[ 3 ] = 1. - bc[0] - bc[1] - bc[2];
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};
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];
494 solveSystemOfEquations2<6,2>(matrix2,res,std::numeric_limits<double>::min());
496 refCoo[0]=computeTria6RefBase(res,p);
497 refCoo[1]=computeTria6RefBase(res+6,p);
498 computeWeightedCoeffsInTria6FromRefBase(refCoo,bc);
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++)
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];
526 solveSystemOfEquations2<10,3>(matrix2,res,std::numeric_limits<double>::min());
528 refCoo[0]=computeTetra10RefBase(res,p);
529 refCoo[1]=computeTetra10RefBase(res+10,p);
530 refCoo[2]=computeTetra10RefBase(res+20,p);
531 computeWeightedCoeffsInTetra10FromRefBase(refCoo,bc);
535 throw INTERP_KERNEL::Exception("INTERP_KERNEL::barycentric_coords : unrecognized simplex !");
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.
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.
557 inline void quad_mapped_coords(const std::vector<const double*>& n, const double *p, double *bc)
559 double prec = 1.0e-14;
560 enum { _XX=0, _YY, _ZZ };
563 throw INTERP_KERNEL::Exception("INTERP_KERNEL::quad_mapped_coords : unrecognized geometric type! Only QUAD4 supported.");
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]};
571 // degenerated case: a rectangle:
572 if (fabs(N[0]) < prec && fabs(N[1]) < prec)
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;
581 double b,c ,a = A[1]*N[0]-A[0]*N[1];
583 if (fabs(a) > 1.0e-14)
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];
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];
596 double delta = b*b - 4.0*a*c;
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!");
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!");
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!");
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]
630 * - for HEXA8, point F (5) is taken to be the origin (see med file ref connec):
644 inline void cuboid_mapped_coords(const std::vector<const double*>& n, const double *p, double *bc)
646 double prec = 1.0e-14;
649 throw INTERP_KERNEL::Exception("INTERP_KERNEL::cuboid_mapped_coords: unrecognized geometric type! Only HEXA8 supported.");
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]);
655 dy1 = OrthoDistanceFromPtToPlaneInSpaceDim3(p, n[5],n[6],n[2]);
656 dy2 = OrthoDistanceFromPtToPlaneInSpaceDim3(p, n[4],n[0],n[3]);
658 dz1 = OrthoDistanceFromPtToPlaneInSpaceDim3(p, n[5],n[4],n[7]);
659 dz2 = OrthoDistanceFromPtToPlaneInSpaceDim3(p, n[1],n[2],n[3]);
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");
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);
669 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
670 /* calcul la surface d'un polygone. */
671 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
673 inline double Surf_Poly(const std::vector<double>& Poly)
677 for(unsigned long i=0; i<(Poly.size())/2-2; i++)
679 double Surf=Surf_Tri( &Poly[0],&Poly[2*(i+1)],&Poly[2*(i+2)] );
680 Surface=Surface + Surf ;
685 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
686 /* fonction qui teste si un point est dans une maille */
688 /* P_1, P_2, P_3 sommet des mailles */
689 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
691 inline bool point_dans_triangle(const double* P_0,const double* P_1,
692 const double* P_2,const double* P_3,
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) )
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 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
713 inline void verif_point_dans_vect(const double* P, std::vector<double>& V, double absolute_precision )
715 long taille=V.size();
716 bool isPresent=false;
717 for(long i=0;i<taille/2;i++)
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)
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 */
735 /*sommets de P: P_1, P_2, P_3 */
736 /*sommets de S: P_4, P_5, P_6 */
737 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
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)
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);
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);
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);
753 verif_point_dans_vect(P_3,V,absolute_precision);
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 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
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)
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]);
771 double absolute_precision = dim_caracteristic*precision;
772 if(fabs(det)>absolute_precision)
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;
776 if (k_1 >= -absolute_precision && k_1 <= 1+absolute_precision)
777 //if( k_1 >= -precision && k_1 <= 1+precision)
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;
781 if (k_2 >= -absolute_precision && k_2 <= 1+absolute_precision)
782 //if( k_2 >= -precision && k_2 <= 1+precision)
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);
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 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
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)
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);
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 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
823 inline void verif_maill_dans_vect(int Num, std::vector<int>& V)
825 long taille=V.size();
827 for(long i=0;i<taille;i++)
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 */
844 bool operator()(std::pair<double,double>theta1, std::pair<double,double> theta2) const
846 double norm1 = sqrt(theta1.first*theta1.first +theta1.second*theta1.second);
847 double norm2 = sqrt(theta2.first*theta2.first +theta2.second*theta2.second);
849 double epsilon = 1.e-12;
851 if( norm1 < epsilon || norm2 < epsilon )
852 std::cout << "Warning InterpolationUtils.hxx: AngleLess : Vector with zero norm, cannot define the angle !!!! " << std::endl;
854 return theta1.second*(norm2 + theta2.first) < theta2.second*(norm1 + theta1.first);
860 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
861 /* fonction pour reconstituer un polygone convexe a partir */
862 /* d'un nuage de point. */
863 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
865 inline std::vector<double> reconstruct_polygon(const std::vector<double>& V)
868 int taille((int)V.size());
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);
883 for(int i=0; i<taille/2-1;i++)
885 std::vector<double> Trigo=calcul_cos_et_sin(&Bary[0],&V[0],&V[2*(i+1)]);
889 // {angle[i+1]=atan2(SIN[i+1],COS[i+1]);}
891 // {angle[i+1]=-atan2(SIN[i+1],COS[i+1]);}
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++)
901 // Ordre.insert(std::make_pair(angle[i],i));
902 CosSin.insert(std::make_pair(std::make_pair(SIN[i],COS[i]),i));
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++)
908 // int j=(*mi).second;
909 // Pt_ordonne.push_back(V[2*j]);
910 // Pt_ordonne.push_back(V[2*j+1]);
912 for(micossin=CosSin.begin();micossin!=CosSin.end();micossin++)
914 int j=(*micossin).second;
915 Pt_ordonne.push_back(V[2*j]);
916 Pt_ordonne.push_back(V[2*j+1]);
925 template<int DIM, NumberingPolicy numPol, class MyMeshType>
926 inline void getElemBB(double* bb, const double *coordsOfMesh, int iP, int nb_nodes)
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();
935 for (int i=0; i<nb_nodes; i++)
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];
950 * Find a vector orthogonal to the input vector
952 inline void orthogonalVect3(const double inpVect[3], double outVect[3])
954 std::vector<bool> sw(3,false);
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)));
959 std::size_t posMax(std::distance(inpVect2,std::max_element(inpVect2,inpVect2+3)));
961 { posMax=(posMin+1)%3; }
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];
967 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
968 /* Computes the dot product of a and b */
969 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
971 inline double dotprod( const double * a, const double * b)
974 for(int idim = 0; idim < dim ; idim++) result += a[idim]*b[idim];
977 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
978 /* Computes the norm of vector v */
979 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
981 inline double norm(const double * v)
984 for(int idim =0; idim<dim; idim++) result+=v[idim]*v[idim];
987 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
988 /* Computes the square norm of vector a-b */
989 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
991 inline double distance2( const double * a, const double * b)
994 for(int idim =0; idim<dim; idim++) result+=(a[idim]-b[idim])*(a[idim]-b[idim]);
997 template<class T, int dim>
998 inline double distance2( T * a, int inda, T * b, int indb)
1001 for(int idim =0; idim<dim; idim++) result += ((*a)[inda+idim] - (*b)[indb+idim])* ((*a)[inda+idim] - (*b)[indb+idim]);
1004 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
1005 /* Computes the determinant of a and b */
1006 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
1007 inline double determinant ( double * a, double * b)
1009 return a[0]*b[1]-a[1]*b[0];
1011 inline double determinant ( double * a, double * b, double * c)
1013 return a[0]*determinant(b+1,c+1)-b[0]*determinant(a+1,c+1)+c[0]*determinant(a+1,b+1);
1016 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
1017 /* Computes the cross product of AB and AC */
1018 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
1020 template<int dim> inline void crossprod( const double * A, const double * B, const double * C, double * V);
1023 void crossprod<2>( const double * A, const double * B, const double * C, double * V)
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;
1030 V[0]=determinant(AB,AC);
1034 void crossprod<3>( const double * A, const double * B, const double * C, double * V)
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;
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];
1046 void crossprod<1>( const double * /*A*/, const double * /*B*/, const double * /*C*/, double * /*V*/)
1048 // just to be able to compile
1051 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
1052 /* Checks wether point A is inside the quadrangle BCDE */
1053 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
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)
1058 crossprod<dim>(A,B,C,ABC);
1059 crossprod<dim>(A,D,E,ADE);
1060 return dotprod<dim>(ABC,ADE);
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)
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;
1076 double normAB= norm<dim>(AB);
1077 for(int idim =0; idim<dim; idim++) AB[idim]/=normAB;
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];
1083 double denom= normAC+AB_dot_AC;
1084 double numer=norm<dim>(orthAB);
1086 return 2*atan2(numer,denom);
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);
1094 double direct_frame<2>(const double * A, const double * B, const double * C, double * n)
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;
1101 return determinant(AB,AC)*n[0];
1104 double direct_frame<3>(const double * A, const double * B, const double * C, double * n)
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;
1111 return determinant(AB,AC,n)>0;
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)
1124 for(int i_A = 1; i_A<nb_NodesA-1; i_A++)
1126 for(int i_B = 1; i_B<nb_NodesB-1; i_B++)
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);
1133 int nb_inter=((int)inter.size())/DIM;
1134 if(nb_inter >3) inter=INTERP_KERNEL::reconstruct_polygon(inter);
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)
1146 for(int i = 1; i<(int)inter.size()/DIM-1; i++)
1148 INTERP_KERNEL::crossprod<DIM>(&inter[0],&inter[DIM*i],&inter[DIM*(i+1)],area);
1149 result +=0.5*norm<DIM>(area);
1154 template<int DIM> inline double polygon_area(std::deque<double>& inter)
1159 for(int i = 1; i<(int)inter.size()/DIM-1; i++)
1161 INTERP_KERNEL::crossprod<DIM>(&inter[0],&inter[DIM*i],&inter[DIM*(i+1)],area);
1162 result +=0.5*norm<DIM>(area);
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)
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];
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)
1197 int i1next = ( i1 + 1 ) % size1;
1198 int i2next = ( i2 + sign +size2) % size2;
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;
1205 if(i1next == istart1)
1207 if(i2next == istart2)
1212 if(i2next == istart2)
1216 if(distance2<T,dim>(L1,i1next*dim, L2,i2next*dim) > epsilon )
1222 i1next = ( i1 + 1 ) % size1;
1223 i2next = ( i2 + sign + size2 ) % size2;
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)
1234 if(L1==NULL || L2==NULL)
1236 std::cout << "Warning InterpolationUtils.hxx:checkEqualPolygonsPointer: Null pointer " << std::endl;
1237 throw(Exception("big error: not closed polygon..."));
1240 int size1 = (*L1).size()/dim;
1241 int size2 = (*L2).size()/dim;
1245 while( istart2 < size2 && distance2<T,dim>(L1,istart1*dim, L2,istart2*dim) > epsilon ) istart2++;
1247 if(istart2 == size2)
1249 return (size1 == 0) && (size2 == 0);
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);