1 // Copyright (C) 2007-2020 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 return sqrt(X*X+Y*Y);
104 inline void mid_of_seg2(const double P1[2], const double P2[2], double MID[2])
106 MID[0]=(P1[0]+P2[0])/2.;
107 MID[1]=(P1[1]+P1[1])/2.;
110 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
111 /* calcul le cos et le sin de l'angle P1P2,P1P3 */
112 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
114 inline std::vector<double> calcul_cos_et_sin(const double* P_1,
119 std::vector<double> Vect;
120 double P1_P2=norme_vecteur(P_1,P_2);
121 double P2_P3=norme_vecteur(P_2,P_3);
122 double P3_P1=norme_vecteur(P_3,P_1);
124 double N=P1_P2*P1_P2+P3_P1*P3_P1-P2_P3*P2_P3;
125 double D=2.0*P1_P2*P3_P1;
127 if (COS>1.0) COS=1.0;
128 if (COS<-1.0) COS=-1.0;
130 double V=mon_determinant(P_2,P_3,P_1);
131 double D_1=P1_P2*P3_P1;
133 if (SIN>1.0) SIN=1.0;
134 if (SIN<-1.0) SIN=-1.0;
142 * This method builds a quadrangle built with the first point of 'triIn' the barycenter of two edges starting or ending with
143 * the first point of 'triIn' and the barycenter of 'triIn'.
145 * @param triIn is a 6 doubles array in full interlace mode, that represents a triangle.
146 * @param quadOut is a 8 doubles array filled after the following call.
148 template<int SPACEDIM>
149 inline void fillDualCellOfTri(const double *triIn, double *quadOut)
152 std::copy(triIn,triIn+SPACEDIM,quadOut);
153 double tmp[SPACEDIM];
154 std::transform(triIn,triIn+SPACEDIM,triIn+SPACEDIM,tmp,std::plus<double>());
156 std::transform(tmp,tmp+SPACEDIM,quadOut+SPACEDIM,std::bind2nd(std::multiplies<double>(),0.5));
157 std::transform(tmp,tmp+SPACEDIM,triIn+2*SPACEDIM,tmp,std::plus<double>());
159 std::transform(tmp,tmp+SPACEDIM,quadOut+2*SPACEDIM,std::bind2nd(std::multiplies<double>(),1/3.));
161 std::transform(triIn,triIn+SPACEDIM,triIn+2*SPACEDIM,tmp,std::plus<double>());
162 std::transform(tmp,tmp+SPACEDIM,quadOut+3*SPACEDIM,std::bind2nd(std::multiplies<double>(),0.5));
166 * 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
167 * the first point of 'triIn' and the barycenter of 'triIn'.
169 * @param triIn is a 6 doubles array in full interlace mode, that represents a triangle.
170 * @param quadOut is a 8 doubles array filled after the following call.
172 template<int SPACEDIM>
173 inline void fillDualCellOfPolyg(const double *polygIn, mcIdType nPtsPolygonIn, double *polygOut)
176 std::copy(polygIn,polygIn+SPACEDIM,polygOut);
177 std::transform(polygIn,polygIn+SPACEDIM,polygIn+SPACEDIM,polygOut+SPACEDIM,std::plus<double>());
179 std::transform(polygOut+SPACEDIM,polygOut+2*SPACEDIM,polygOut+SPACEDIM,std::bind2nd(std::multiplies<double>(),0.5));
180 double tmp[SPACEDIM];
182 for(mcIdType i=0;i<nPtsPolygonIn-2;i++)
184 std::transform(polygIn,polygIn+SPACEDIM,polygIn+(i+2)*SPACEDIM,tmp,std::plus<double>());
185 std::transform(tmp,tmp+SPACEDIM,polygOut+(2*i+3)*SPACEDIM,std::bind2nd(std::multiplies<double>(),0.5));
186 std::transform(polygIn+(i+1)*SPACEDIM,polygIn+(i+2)*SPACEDIM,tmp,tmp,std::plus<double>());
187 std::transform(tmp,tmp+SPACEDIM,polygOut+(2*i+2)*SPACEDIM,std::bind2nd(std::multiplies<double>(),1./3.));
191 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
192 /* calcul les coordonnees du barycentre d'un polygone */
193 /* le vecteur en entree est constitue des coordonnees */
194 /* des sommets du polygone */
195 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
197 inline std::vector<double> bary_poly(const std::vector<double>& V)
199 std::vector<double> Bary;
200 std::size_t taille=V.size();
204 for(std::size_t i=0;i<taille/2;i++)
209 double A=2*x/(static_cast<double>(taille));
210 double B=2*y/(static_cast<double>(taille));
211 Bary.push_back(A);//taille vecteur=2*nb de points.
219 * Given 6 coeffs of a Tria6 returns the corresponding value of a given pos
221 inline double computeTria6RefBase(const double *coeffs, const double *pos)
223 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];
227 * Given xsi,eta in refCoo (length==2) return 6 coeffs in weightedPos.
229 inline void computeWeightedCoeffsInTria6FromRefBase(const double *refCoo, double *weightedPos)
231 weightedPos[0]=(1.-refCoo[0]-refCoo[1])*(1.-2*refCoo[0]-2.*refCoo[1]);
232 weightedPos[1]=refCoo[0]*(2.*refCoo[0]-1.);
233 weightedPos[2]=refCoo[1]*(2.*refCoo[1]-1.);
234 weightedPos[3]=4.*refCoo[0]*(1.-refCoo[0]-refCoo[1]);
235 weightedPos[4]=4.*refCoo[0]*refCoo[1];
236 weightedPos[5]=4.*refCoo[1]*(1.-refCoo[0]-refCoo[1]);
240 * Given 10 coeffs of a Tetra10 returns the corresponding value of a given pos
242 inline double computeTetra10RefBase(const double *coeffs, const double *pos)
244 return coeffs[0]+coeffs[1]*pos[0]+coeffs[2]*pos[1]+coeffs[3]*pos[2]+
245 coeffs[4]*pos[0]*pos[0]+coeffs[5]*pos[0]*pos[1]+coeffs[6]*pos[0]*pos[2]+
246 coeffs[7]*pos[1]*pos[1]+coeffs[8]*pos[1]*pos[2]+coeffs[9]*pos[2]*pos[2];
250 * Given xsi,eta,z in refCoo (length==3) return 10 coeffs in weightedPos.
252 inline void computeWeightedCoeffsInTetra10FromRefBase(const double *refCoo, double *weightedPos)
254 //http://www.cadfamily.com/download/CAE/ABAQUS/The%20Finite%20Element%20Method%20-%20A%20practical%20course%20abaqus.pdf page 217
255 //L1=1-refCoo[0]-refCoo[1]-refCoo[2]
256 //L2=refCoo[0] L3=refCoo[1] L4=refCoo[2]
257 weightedPos[0]=(-2.*(refCoo[0]+refCoo[1]+refCoo[2])+1)*(1-refCoo[0]-refCoo[1]-refCoo[2]);//(2*L1-1)*L1
258 weightedPos[1]=(2.*refCoo[0]-1.)*refCoo[0];//(2*L2-1)*L2
259 weightedPos[2]=(2.*refCoo[1]-1.)*refCoo[1];//(2*L3-1)*L3
260 weightedPos[3]=(2.*refCoo[2]-1.)*refCoo[2];//(2*L4-1)*L4
261 weightedPos[4]=4.*(1-refCoo[0]-refCoo[1]-refCoo[2])*refCoo[0];//4*L1*L2
262 weightedPos[5]=4.*refCoo[0]*refCoo[1];//4*L2*L3
263 weightedPos[6]=4.*(1-refCoo[0]-refCoo[1]-refCoo[2])*refCoo[1];//4*L1*L3
264 weightedPos[7]=4.*(1-refCoo[0]-refCoo[1]-refCoo[2])*refCoo[2];//4*L1*L4
265 weightedPos[8]=4.*refCoo[0]*refCoo[2];//4*L2*L4
266 weightedPos[9]=4.*refCoo[1]*refCoo[2];//4*L3*L4
270 * \brief Solve system equation in matrix form using Gaussian elimination algorithm
271 * \param M - N x N+1 matrix
272 * \param sol - vector of N solutions
273 * \retval bool - true if succeeded
275 template<unsigned nbRow>
276 bool solveSystemOfEquations(double M[nbRow][nbRow+1], double* sol)
278 const int nbCol=nbRow+1;
280 // make upper triangular matrix (forward elimination)
282 int iR[nbRow];// = { 0, 1, 2 };
283 for ( int i = 0; i < (int) nbRow; ++i )
285 for ( int i = 0; i < (int)(nbRow-1); ++i ) // nullify nbRow-1 rows
287 // swap rows to have max value of i-th column in i-th row
288 double max = std::fabs( M[ iR[i] ][i] );
289 for ( int r = i+1; r < (int)nbRow; ++r )
291 double m = std::fabs( M[ iR[r] ][i] );
295 std::swap( iR[r], iR[i] );
298 if ( max < std::numeric_limits<double>::min() )
300 //sol[0]=1; sol[1]=sol[2]=sol[3]=0;
301 return false; // no solution
303 // make 0 below M[i][i] (actually we do not modify i-th column)
304 double* tUpRow = M[ iR[i] ];
305 for ( int r = i+1; r < (int)nbRow; ++r )
307 double* mRow = M[ iR[r] ];
308 double coef = mRow[ i ] / tUpRow[ i ];
309 for ( int c = i+1; c < nbCol; ++c )
310 mRow[ c ] -= tUpRow[ c ] * coef;
313 double* mRow = M[ iR[nbRow-1] ];
314 if ( std::fabs( mRow[ nbRow-1 ] ) < std::numeric_limits<double>::min() )
316 //sol[0]=1; sol[1]=sol[2]=sol[3]=0;
317 return false; // no solution
319 mRow[ nbRow ] /= mRow[ nbRow-1 ];
321 // calculate solution (back substitution)
323 sol[ nbRow-1 ] = mRow[ nbRow ];
325 for ( int i = nbRow-2; i+1; --i )
328 sol[ i ] = mRow[ nbRow ];
329 for ( int j = nbRow-1; j > i; --j )
330 sol[ i ] -= sol[j]*mRow[ j ];
331 sol[ i ] /= mRow[ i ];
339 * \brief Solve system equation in matrix form using Gaussian elimination algorithm
340 * \param M - N x N+NB_OF_VARS matrix
341 * \param sol - vector of N solutions
342 * \retval bool - true if succeeded
344 template<unsigned SZ, unsigned NB_OF_RES>
345 bool solveSystemOfEquations2(const double *matrix, double *solutions, double eps)
352 double B[SZ*(SZ+NB_OF_RES)];
353 std::copy(matrix,matrix+SZ*(SZ+NB_OF_RES),B);
365 if(fabs(B[nr*k+n])>eps)
366 {/* Rows permutation */
368 std::swap(B[nr*k+m],B[nr*n+m]);
373 s=B[np];//s is the Pivot
374 std::transform(B+k*nr,B+(k+1)*nr,B+k*nr,std::bind2nd(std::divides<double>(),s));
381 B[j*nr+mb]-=B[k*nr+mb]*g;
385 for(j=0;j<NB_OF_RES;j++)
387 solutions[j*SZ+k]=B[nr*k+SZ+j];
392 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
393 /* Calculate barycentric coordinates of a 2D point p */
394 /* with respect to the triangle verices. */
395 /* triaCoords are in full interlace */
396 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
398 template<int SPACEDIM>
399 inline void barycentric_coords(const double* triaCoords, const double* p, double* bc)
403 T11 = triaCoords[0]-triaCoords[2*SPACEDIM], T12 = triaCoords[SPACEDIM]-triaCoords[2*SPACEDIM],
404 T21 = triaCoords[1]-triaCoords[2*SPACEDIM+1], T22 = triaCoords[SPACEDIM+1]-triaCoords[2*SPACEDIM+1];
405 // matrix determinant
406 double Tdet = T11*T22 - T12*T21;
407 if ( fabs( Tdet ) < std::numeric_limits<double>::min() ) {
408 bc[0]=1; bc[1]=0; bc[2]=0;
412 double t11 = T22, t12 = -T12, t21 = -T21, t22 = T11;
414 double r11 = p[0]-triaCoords[2*SPACEDIM], r12 = p[1]-triaCoords[2*SPACEDIM+1];
415 // barycentric coordinates: multiply matrix by vector
416 bc[0] = (t11 * r11 + t12 * r12)/Tdet;
417 bc[1] = (t21 * r11 + t22 * r12)/Tdet;
418 bc[2] = 1. - bc[0] - bc[1];
421 inline void barycentric_coords_seg2(const std::vector<const double*>& n, const double *p, double *bc)
423 double delta=n[0][0]-n[1][0];
424 bc[0]=fabs((*p-n[1][0])/delta);
425 bc[1]=fabs((*p-n[0][0])/delta);
428 inline void barycentric_coords_tri3(const std::vector<const double*>& n, const double *p, double *bc)
430 enum { _XX=0, _YY, _ZZ };
433 T11 = n[0][_XX]-n[2][_XX], T12 = n[1][_XX]-n[2][_XX],
434 T21 = n[0][_YY]-n[2][_YY], T22 = n[1][_YY]-n[2][_YY];
435 // matrix determinant
436 double Tdet = T11*T22 - T12*T21;
437 if ( (std::fabs( Tdet) ) < (std::numeric_limits<double>::min()) )
439 bc[0]=1; bc[1]=bc[2]=0; // no solution
443 double t11 = T22, t12 = -T12, t21 = -T21, t22 = T11;
445 double r11 = p[_XX]-n[2][_XX], r12 = p[_YY]-n[2][_YY];
446 // barycentric coordinates: multiply matrix by vector
447 bc[0] = (t11 * r11 + t12 * r12)/Tdet;
448 bc[1] = (t21 * r11 + t22 * r12)/Tdet;
449 bc[2] = 1. - bc[0] - bc[1];
452 inline void barycentric_coords_quad4(const std::vector<const double*>& n, const double *p, double *bc)
454 enum { _XX=0, _YY, _ZZ };
455 // Find bc by solving system of 3 equations using Gaussian elimination algorithm
456 // bc1*( x1 - x4 ) + bc2*( x2 - x4 ) + bc3*( x3 - x4 ) = px - x4
457 // bc1*( y1 - y4 ) + bc2*( y2 - y4 ) + bc3*( y3 - y4 ) = px - y4
458 // bc1*( z1 - z4 ) + bc2*( z2 - z4 ) + bc3*( z3 - z4 ) = px - z4
461 {{ n[0][_XX]-n[3][_XX], n[1][_XX]-n[3][_XX], n[2][_XX]-n[3][_XX], p[_XX]-n[3][_XX] },
462 { n[0][_YY]-n[3][_YY], n[1][_YY]-n[3][_YY], n[2][_YY]-n[3][_YY], p[_YY]-n[3][_YY] },
463 { n[0][_ZZ]-n[3][_ZZ], n[1][_ZZ]-n[3][_ZZ], n[2][_ZZ]-n[3][_ZZ], p[_ZZ]-n[3][_ZZ] }};
465 if ( !solveSystemOfEquations<3>( T, bc ) )
466 bc[0]=1., bc[1] = bc[2] = bc[3] = 0;
468 bc[ 3 ] = 1. - bc[0] - bc[1] - bc[2];
471 inline void barycentric_coords_tri6(const std::vector<const double*>& n, const double *p, double *bc)
473 double matrix2[48]={1., 0., 0., 0., 0., 0., 0., 0.,
474 1., 0., 0., 0., 0., 0., 1., 0.,
475 1., 0., 0., 0., 0., 0., 0., 1.,
476 1., 0., 0., 0., 0., 0., 0.5, 0.,
477 1., 0., 0., 0., 0., 0., 0.5, 0.5,
478 1., 0., 0., 0., 0., 0., 0.,0.5};
481 matrix2[8*i+1]=n[i][0];
482 matrix2[8*i+2]=n[i][1];
483 matrix2[8*i+3]=n[i][0]*n[i][0];
484 matrix2[8*i+4]=n[i][0]*n[i][1];
485 matrix2[8*i+5]=n[i][1]*n[i][1];
488 solveSystemOfEquations2<6,2>(matrix2,res,std::numeric_limits<double>::min());
490 refCoo[0]=computeTria6RefBase(res,p);
491 refCoo[1]=computeTria6RefBase(res+6,p);
492 computeWeightedCoeffsInTria6FromRefBase(refCoo,bc);
495 inline void barycentric_coords_tetra10(const std::vector<const double*>& n, const double *p, double *bc)
497 double matrix2[130]={1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
498 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.,
499 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1., 0.,
500 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 1.,
501 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.5, 0., 0.,
502 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.5, 0.5, 0.,
503 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,0.5, 0.,
504 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.5,
505 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.5, 0., 0.5,
506 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.5, 0.5};
507 for(int i=0;i<10;i++)
509 matrix2[13*i+1]=n[i][0];
510 matrix2[13*i+2]=n[i][1];
511 matrix2[13*i+3]=n[i][2];
512 matrix2[13*i+4]=n[i][0]*n[i][0];
513 matrix2[13*i+5]=n[i][0]*n[i][1];
514 matrix2[13*i+6]=n[i][0]*n[i][2];
515 matrix2[13*i+7]=n[i][1]*n[i][1];
516 matrix2[13*i+8]=n[i][1]*n[i][2];
517 matrix2[13*i+9]=n[i][2]*n[i][2];
520 solveSystemOfEquations2<10,3>(matrix2,res,std::numeric_limits<double>::min());
522 refCoo[0]=computeTetra10RefBase(res,p);
523 refCoo[1]=computeTetra10RefBase(res+10,p);
524 refCoo[2]=computeTetra10RefBase(res+20,p);
525 computeWeightedCoeffsInTetra10FromRefBase(refCoo,bc);
529 * Calculate barycentric coordinates of a point p with respect to triangle or tetra vertices.
530 * This method makes 2 assumptions :
531 * - this is a simplex
532 * - 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.
533 * If not the case (3D surf for example) a previous projection should be done before.
535 inline void barycentric_coords(const std::vector<const double*>& n, const double *p, double *bc)
541 barycentric_coords_seg2(n,p,bc);
546 barycentric_coords_tri3(n,p,bc);
551 barycentric_coords_quad4(n,p,bc);
557 barycentric_coords_tri6(n,p,bc);
562 barycentric_coords_tetra10(n,p,bc);
566 throw INTERP_KERNEL::Exception("INTERP_KERNEL::barycentric_coords : unrecognized simplex !");
570 inline void barycentric_coords(INTERP_KERNEL::NormalizedCellType ct, const std::vector<const double*>& n, const double *p, double *bc)
576 barycentric_coords_seg2(n,p,bc);
581 barycentric_coords_tri3(n,p,bc);
586 barycentric_coords_quad4(n,p,bc);
592 barycentric_coords_tri6(n,p,bc);
597 barycentric_coords_tetra10(n,p,bc);
601 throw INTERP_KERNEL::Exception("INTERP_KERNEL::barycentric_coords : unrecognized simplex !");
606 * Compute barycentric coords of point \a point relative to segment S [segStart,segStop] in 3D space.
607 * If point \a point is further from S than eps false is returned.
608 * If point \a point projection on S is outside S false is also returned.
609 * If point \a point is closer from S than eps and its projection inside S true is returned and \a bc0 and \a bc1 output parameter set.
611 inline bool IsPointOn3DSeg(const double segStart[3], const double segStop[3], const double point[3], double eps, double& bc0, double& bc1)
613 double AB[3]={segStop[0]-segStart[0],segStop[1]-segStart[1],segStop[2]-segStart[2]},AP[3]={point[0]-segStart[0],point[1]-segStart[1],point[2]-segStart[2]};
614 double l_AB(sqrt(AB[0]*AB[0]+AB[1]*AB[1]+AB[2]*AB[2]));
615 double AP_dot_AB((AP[0]*AB[0]+AP[1]*AB[1]+AP[2]*AB[2])/(l_AB*l_AB));
616 double projOfPOnAB[3]={segStart[0]+AP_dot_AB*AB[0],segStart[1]+AP_dot_AB*AB[1],segStart[2]+AP_dot_AB*AB[2]};
617 double V_dist_P_AB[3]={point[0]-projOfPOnAB[0],point[1]-projOfPOnAB[1],point[2]-projOfPOnAB[2]};
618 double dist_P_AB(sqrt(V_dist_P_AB[0]*V_dist_P_AB[0]+V_dist_P_AB[1]*V_dist_P_AB[1]+V_dist_P_AB[2]*V_dist_P_AB[2]));
620 return false;//to far from segment [segStart,segStop]
621 if(AP_dot_AB<-eps || AP_dot_AB>1.+eps)
623 AP_dot_AB=std::max(AP_dot_AB,0.); AP_dot_AB=std::min(AP_dot_AB,1.);
624 bc0=1.-AP_dot_AB; bc1=AP_dot_AB;
629 * Calculate pseudo barycentric coordinates of a point p with respect to the quadrangle vertices.
630 * This method makes the assumption that:
631 * - spacedim == meshdim (2 here).
632 * - the point is within the quad
633 * Quadratic elements are not supported yet.
635 * A quadrangle can be described as 3 vectors, one point being taken as the origin.
636 * Denoting A, B, C the three other points, any point P within the quad is written as
637 * P = xA+ yC + xy(B-A-C)
638 * This method solve those 2 equations (one per component) for x and y.
646 inline void quad_mapped_coords(const std::vector<const double*>& n, const double *p, double *bc)
648 double prec = 1.0e-14;
649 enum { _XX=0, _YY, _ZZ };
652 throw INTERP_KERNEL::Exception("INTERP_KERNEL::quad_mapped_coords : unrecognized geometric type! Only QUAD4 supported.");
654 double A[2] = {n[1][_XX] - n[0][_XX], n[1][_YY] - n[0][_YY]};
655 double B[2] = {n[2][_XX] - n[0][_XX], n[2][_YY] - n[0][_YY]};
656 double C[2] = {n[3][_XX] - n[0][_XX], n[3][_YY] - n[0][_YY]};
657 double N[2] = {B[_XX] - A[_XX] - C[_XX], B[_YY] - A[_YY] - C[_YY]};
658 double P[2] = {p[_XX] - n[0][_XX], p[_YY] - n[0][_YY]};
660 // degenerated case: a rectangle:
661 if (fabs(N[0]) < prec && fabs(N[1]) < prec)
663 double det = C[0]*A[1] -C[1]*A[0];
664 if (fabs(det) < prec)
665 throw INTERP_KERNEL::Exception("MappedBarycentric intersection type: quad_mapped_coords() has a degenerated 2x2 system!");
666 bc[0] = (P[0]*A[1]-P[1]*A[0])/det;
667 bc[1] = (P[1]*C[0]-P[0]*C[1])/det;
670 double b,c ,a = A[1]*N[0]-A[0]*N[1];
672 if (fabs(a) > 1.0e-14)
674 b = A[1]*C[0]+N[1]*P[0]-N[0]*P[1]-A[0]*C[1];
675 c = P[0]*C[1] - P[1]*C[0];
680 a = -C[1]*N[0]+C[0]*N[1];
681 b = A[1]*C[0]-N[1]*P[0]+N[0]*P[1]-A[0]*C[1];
682 c = -P[0]*A[1] + P[1]*A[0];
685 double delta = b*b - 4.0*a*c;
687 throw INTERP_KERNEL::Exception("MappedBarycentric intersection type: quad_mapped_coords(): imaginary solutions!");
688 bc[1] = 0.5*(-b+sqrt(delta))/a;
689 if (bc[1] < -prec || bc[1] > (1.0+prec))
690 bc[1] = 0.5*(-b-sqrt(delta))/a;
691 if (bc[1] < -prec || bc[1] > (1.0+prec))
692 throw INTERP_KERNEL::Exception("MappedBarycentric intersection type: quad_mapped_coords(): point doesn't seem to be in quad4!");
695 double denom = C[0]+bc[1]*N[0];
696 if (fabs(denom) < prec)
697 throw INTERP_KERNEL::Exception("MappedBarycentric intersection type: quad_mapped_coords(): point doesn't seem to be in quad4!");
698 bc[0] = (P[0]-bc[1]*A[0])/denom;
699 if (bc[0] < -prec || bc[0] > (1.0+prec))
700 throw INTERP_KERNEL::Exception("MappedBarycentric intersection type: quad_mapped_coords(): point doesn't seem to be in quad4!");
705 double denom = A[1]+bc[0]*N[1];
706 if (fabs(denom) < prec)
707 throw INTERP_KERNEL::Exception("MappedBarycentric intersection type: cuboid_mapped_coord(): point doesn't seem to be in quad4!");
708 bc[1] = (P[1]-bc[0]*C[1])/denom;
709 if (bc[1] < -prec || bc[1] > (1.0+prec))
710 throw INTERP_KERNEL::Exception("MappedBarycentric intersection type: cuboid_mapped_coord(): point doesn't seem to be in quad4!");
715 * Doing as in quad_mapped_coords() would lead to a 4th order equation ... So go simpler here:
716 * orthogonal distance to each pair of parallel faces is computed. The ratio gives a number in [0,1]
719 * - for HEXA8, point F (5) is taken to be the origin (see med file ref connec):
733 inline void cuboid_mapped_coords(const std::vector<const double*>& n, const double *p, double *bc)
735 double prec = 1.0e-14;
738 throw INTERP_KERNEL::Exception("INTERP_KERNEL::cuboid_mapped_coords: unrecognized geometric type! Only HEXA8 supported.");
740 double dx1, dx2, dy1, dy2, dz1, dz2;
741 dx1 = OrthoDistanceFromPtToPlaneInSpaceDim3(p, n[4],n[5],n[1]);
742 dx2 = OrthoDistanceFromPtToPlaneInSpaceDim3(p, n[7],n[3],n[2]);
744 dy1 = OrthoDistanceFromPtToPlaneInSpaceDim3(p, n[5],n[6],n[2]);
745 dy2 = OrthoDistanceFromPtToPlaneInSpaceDim3(p, n[4],n[0],n[3]);
747 dz1 = OrthoDistanceFromPtToPlaneInSpaceDim3(p, n[5],n[4],n[7]);
748 dz2 = OrthoDistanceFromPtToPlaneInSpaceDim3(p, n[1],n[2],n[3]);
750 if (dx1 < -prec || dx2 < -prec || dy1 < -prec || dy2 < -prec || dz1 < -prec || dz2 < -prec)
751 throw INTERP_KERNEL::Exception("INTERP_KERNEL::cuboid_mapped_coords: point outside HEXA8");
753 bc[0] = dx1+dx2 < prec ? 0.5 : dx1/(dx1+dx2);
754 bc[1] = dy1+dy2 < prec ? 0.5 : dy1/(dy1+dy2);
755 bc[2] = dz1+dz2 < prec ? 0.5 : dz1/(dz1+dz2);
758 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
759 /* calcul la surface d'un polygone. */
760 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
762 inline double Surf_Poly(const std::vector<double>& Poly)
766 for(unsigned long i=0; i<(Poly.size())/2-2; i++)
768 double Surf=Surf_Tri( &Poly[0],&Poly[2*(i+1)],&Poly[2*(i+2)] );
769 Surface=Surface + Surf ;
774 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
775 /* fonction qui teste si un point est dans une maille */
777 /* P_1, P_2, P_3 sommet des mailles */
778 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
780 inline bool point_dans_triangle(const double* P_0,const double* P_1,
781 const double* P_2,const double* P_3,
786 double det_1=mon_determinant(P_1,P_3,P_0);
787 double det_2=mon_determinant(P_3,P_2,P_0);
788 double det_3=mon_determinant(P_2,P_1,P_0);
789 if( (det_1>=-eps && det_2>=-eps && det_3>=-eps) || (det_1<=eps && det_2<=eps && det_3<=eps) )
797 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
798 /*fonction pour verifier qu'un point n'a pas deja ete considerer dans */
799 /* le vecteur et le rajouter au vecteur sinon. */
800 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
802 inline void verif_point_dans_vect(const double* P, std::vector<double>& V, double absolute_precision )
804 std::size_t taille=V.size();
805 bool isPresent=false;
806 for(std::size_t i=0;i<taille/2;i++)
808 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)
820 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
821 /* fonction qui rajoute les sommet du triangle P dans le vecteur V */
822 /* si ceux-ci sont compris dans le triangle S et ne sont pas deja dans */
824 /*sommets de P: P_1, P_2, P_3 */
825 /*sommets de S: P_4, P_5, P_6 */
826 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
828 inline void rajou_sommet_triangl(const double* P_1,const double* P_2,const double* P_3,
829 const double* P_4,const double* P_5,const double* P_6,
830 std::vector<double>& V, double dim_caracteristic, double precision)
833 double absolute_precision = precision*dim_caracteristic;
834 bool A_1=INTERP_KERNEL::point_dans_triangle(P_1,P_4,P_5,P_6,absolute_precision);
836 verif_point_dans_vect(P_1,V,absolute_precision);
837 bool A_2=INTERP_KERNEL::point_dans_triangle(P_2,P_4,P_5,P_6,absolute_precision);
839 verif_point_dans_vect(P_2,V,absolute_precision);
840 bool A_3=INTERP_KERNEL::point_dans_triangle(P_3,P_4,P_5,P_6,absolute_precision);
842 verif_point_dans_vect(P_3,V,absolute_precision);
846 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
847 /* calcul de l'intersection de deux segments: segments P1P2 avec P3P4 */
848 /* . Si l'intersection est non nulle et si celle-ci n'est */
849 /* n'est pas deja contenue dans Vect on la rajoute a Vect */
850 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
852 inline void inters_de_segment(const double * P_1,const double * P_2,
853 const double * P_3,const double * P_4,
854 std::vector<double>& Vect,
855 double dim_caracteristic, double precision)
857 // calcul du determinant de P_1P_2 et P_3P_4.
858 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]);
860 double absolute_precision = dim_caracteristic*precision;
861 if(fabs(det)>absolute_precision)
863 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;
865 if (k_1 >= -absolute_precision && k_1 <= 1+absolute_precision)
866 //if( k_1 >= -precision && k_1 <= 1+precision)
868 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;
870 if (k_2 >= -absolute_precision && k_2 <= 1+absolute_precision)
871 //if( k_2 >= -precision && k_2 <= 1+precision)
874 P_0[0]=P_1[0]+k_1*(P_2[0]-P_1[0]);
875 P_0[1]=P_1[1]+k_1*(P_2[1]-P_1[1]);
876 verif_point_dans_vect(P_0,Vect,absolute_precision);
884 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
885 /* calcul l'intersection de deux triangles */
886 /* P_1, P_2, P_3: sommets du premier triangle */
887 /* P_4, P_5, P_6: sommets du deuxi�me triangle */
888 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
890 inline void intersec_de_triangle(const double* P_1,const double* P_2, const double* P_3,
891 const double* P_4,const double* P_5,const double* P_6,
892 std::vector<double>& Vect, double dim_caracteristic, double precision)
894 inters_de_segment(P_1,P_2,P_4,P_5,Vect, dim_caracteristic, precision);
895 inters_de_segment(P_1,P_2,P_5,P_6,Vect, dim_caracteristic, precision);
896 inters_de_segment(P_1,P_2,P_6,P_4,Vect, dim_caracteristic, precision);
897 inters_de_segment(P_2,P_3,P_4,P_5,Vect, dim_caracteristic, precision);
898 inters_de_segment(P_2,P_3,P_5,P_6,Vect, dim_caracteristic, precision);
899 inters_de_segment(P_2,P_3,P_6,P_4,Vect, dim_caracteristic, precision);
900 inters_de_segment(P_3,P_1,P_4,P_5,Vect, dim_caracteristic, precision);
901 inters_de_segment(P_3,P_1,P_5,P_6,Vect, dim_caracteristic, precision);
902 inters_de_segment(P_3,P_1,P_6,P_4,Vect, dim_caracteristic, precision);
903 rajou_sommet_triangl(P_1,P_2,P_3,P_4,P_5,P_6,Vect, dim_caracteristic, precision);
904 rajou_sommet_triangl(P_4,P_5,P_6,P_1,P_2,P_3,Vect, dim_caracteristic, precision);
907 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
908 /* fonction pour verifier qu'un node maille n'a pas deja ete considerer */
909 /* dans le vecteur et le rajouter au vecteur sinon. */
910 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
912 inline void verif_maill_dans_vect(int Num, std::vector<int>& V)
914 std::size_t taille=V.size();
916 for(std::size_t i=0;i<taille;i++)
928 /*! Function that compares two angles from the values of the pairs (sin,cos)*/
929 /*! Angles are considered in [0, 2Pi] bt are not computed explicitly */
933 bool operator()(std::pair<double,double>theta1, std::pair<double,double> theta2) const
935 double norm1 = sqrt(theta1.first*theta1.first +theta1.second*theta1.second);
936 double norm2 = sqrt(theta2.first*theta2.first +theta2.second*theta2.second);
938 double epsilon = 1.e-12;
940 if( norm1 < epsilon || norm2 < epsilon )
941 std::cout << "Warning InterpolationUtils.hxx: AngleLess : Vector with zero norm, cannot define the angle !!!! " << std::endl;
943 return theta1.second*(norm2 + theta2.first) < theta2.second*(norm1 + theta1.first);
949 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
950 /* fonction pour reconstituer un polygone convexe a partir */
951 /* d'un nuage de point. */
952 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
954 inline std::vector<double> reconstruct_polygon(const std::vector<double>& V)
957 int taille((int)V.size());
965 double *COS=new double[taille/2];
966 double *SIN=new double[taille/2];
967 //double *angle=new double[taille/2];
968 std::vector<double> Bary=bary_poly(V);
972 for(int i=0; i<taille/2-1;i++)
974 std::vector<double> Trigo=calcul_cos_et_sin(&Bary[0],&V[0],&V[2*(i+1)]);
978 // {angle[i+1]=atan2(SIN[i+1],COS[i+1]);}
980 // {angle[i+1]=-atan2(SIN[i+1],COS[i+1]);}
983 //ensuite on ordonne les angles.
984 std::vector<double> Pt_ordonne;
985 Pt_ordonne.reserve(taille);
986 // std::multimap<double,int> Ordre;
987 std::multimap<std::pair<double,double>,int, AngleLess> CosSin;
988 for(int i=0;i<taille/2;i++)
990 // Ordre.insert(std::make_pair(angle[i],i));
991 CosSin.insert(std::make_pair(std::make_pair(SIN[i],COS[i]),i));
993 // std::multimap <double,int>::iterator mi;
994 std::multimap<std::pair<double,double>,int, AngleLess>::iterator micossin;
995 // for(mi=Ordre.begin();mi!=Ordre.end();mi++)
997 // int j=(*mi).second;
998 // Pt_ordonne.push_back(V[2*j]);
999 // Pt_ordonne.push_back(V[2*j+1]);
1001 for(micossin=CosSin.begin();micossin!=CosSin.end();micossin++)
1003 int j=(*micossin).second;
1004 Pt_ordonne.push_back(V[2*j]);
1005 Pt_ordonne.push_back(V[2*j+1]);
1014 template<int DIM, NumberingPolicy numPol, class MyMeshType>
1015 inline void getElemBB(double* bb, const double *coordsOfMesh, mcIdType iP, int nb_nodes)
1017 bb[0]=std::numeric_limits<double>::max();
1018 bb[1]=-std::numeric_limits<double>::max();
1019 bb[2]=std::numeric_limits<double>::max();
1020 bb[3]=-std::numeric_limits<double>::max();
1021 bb[4]=std::numeric_limits<double>::max();
1022 bb[5]=-std::numeric_limits<double>::max();
1024 for (int i=0; i<nb_nodes; i++)
1026 double x = coordsOfMesh[3*(iP+i)];
1027 double y = coordsOfMesh[3*(iP+i)+1];
1028 double z = coordsOfMesh[3*(iP+i)+2];
1029 bb[0]=(x<bb[0])?x:bb[0];
1030 bb[1]=(x>bb[1])?x:bb[1];
1031 bb[2]=(y<bb[2])?y:bb[2];
1032 bb[3]=(y>bb[3])?y:bb[3];
1033 bb[4]=(z<bb[4])?z:bb[4];
1034 bb[5]=(z>bb[5])?z:bb[5];
1039 * Find a vector orthogonal to the input vector
1041 inline void orthogonalVect3(const double inpVect[3], double outVect[3])
1043 std::vector<bool> sw(3,false);
1045 std::transform(inpVect,inpVect+3,inpVect2,std::ptr_fun<double,double>(fabs));
1046 std::size_t posMin(std::distance(inpVect2,std::min_element(inpVect2,inpVect2+3)));
1048 std::size_t posMax(std::distance(inpVect2,std::max_element(inpVect2,inpVect2+3)));
1050 { posMax=(posMin+1)%3; }
1052 std::size_t posMid(std::distance(sw.begin(),std::find(sw.begin(),sw.end(),false)));
1053 outVect[posMin]=0.; outVect[posMid]=1.; outVect[posMax]=-inpVect[posMid]/inpVect[posMax];
1056 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
1057 /* Computes the dot product of a and b */
1058 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
1060 inline double dotprod( const double * a, const double * b)
1063 for(int idim = 0; idim < dim ; idim++) result += a[idim]*b[idim];
1066 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
1067 /* Computes the norm of vector v */
1068 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
1070 inline double norm(const double * v)
1073 for(int idim =0; idim<dim; idim++) result+=v[idim]*v[idim];
1074 return sqrt(result);
1076 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
1077 /* Computes the square norm of vector a-b */
1078 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
1080 inline double distance2( const double * a, const double * b)
1083 for(int idim =0; idim<dim; idim++) result+=(a[idim]-b[idim])*(a[idim]-b[idim]);
1086 template<class T, int dim>
1087 inline double distance2( T * a, int inda, T * b, int indb)
1090 for(int idim =0; idim<dim; idim++) result += ((*a)[inda+idim] - (*b)[indb+idim])* ((*a)[inda+idim] - (*b)[indb+idim]);
1093 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
1094 /* Computes the determinant of a and b */
1095 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
1096 inline double determinant ( double * a, double * b)
1098 return a[0]*b[1]-a[1]*b[0];
1100 inline double determinant ( double * a, double * b, double * c)
1102 return a[0]*determinant(b+1,c+1)-b[0]*determinant(a+1,c+1)+c[0]*determinant(a+1,b+1);
1105 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
1106 /* Computes the cross product of AB and AC */
1107 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
1109 template<int dim> inline void crossprod( const double * A, const double * B, const double * C, double * V);
1112 void crossprod<2>( const double * A, const double * B, const double * C, double * V)
1116 for(int idim =0; idim<2; idim++) AB[idim] = B[idim]-A[idim];//B-A
1117 for(int idim =0; idim<2; idim++) AC[idim] = C[idim]-A[idim];//C-A;
1119 V[0]=determinant(AB,AC);
1123 void crossprod<3>( const double * A, const double * B, const double * C, double * V)
1127 for(int idim =0; idim<3; idim++) AB[idim] = B[idim]-A[idim];//B-A
1128 for(int idim =0; idim<3; idim++) AC[idim] = C[idim]-A[idim];//C-A;
1130 V[0]=AB[1]*AC[2]-AB[2]*AC[1];
1131 V[1]=-AB[0]*AC[2]+AB[2]*AC[0];
1132 V[2]=AB[0]*AC[1]-AB[1]*AC[0];
1135 void crossprod<1>( const double * /*A*/, const double * /*B*/, const double * /*C*/, double * /*V*/)
1137 // just to be able to compile
1140 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
1141 /* Checks whether point A is inside the quadrangle BCDE */
1142 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
1144 template<int dim> inline double check_inside(const double* A,const double* B,const double* C,const double* D,
1145 const double* E,double* ABC, double* ADE)
1147 crossprod<dim>(A,B,C,ABC);
1148 crossprod<dim>(A,D,E,ADE);
1149 return dotprod<dim>(ABC,ADE);
1153 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
1154 /* Computes the geometric angle (in [0,Pi]) between two non zero vectors AB and AC */
1155 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
1156 template<int dim> inline double angle(const double * A, const double * B, const double * C, double * n)
1162 for(int idim =0; idim<dim; idim++) AB[idim] = B[idim]-A[idim];//B-A;
1163 for(int idim =0; idim<dim; idim++) AC[idim] = C[idim]-A[idim];//C-A;
1165 double normAB= norm<dim>(AB);
1166 for(int idim =0; idim<dim; idim++) AB[idim]/=normAB;
1168 double normAC= norm<dim>(AC);
1169 double AB_dot_AC=dotprod<dim>(AB,AC);
1170 for(int idim =0; idim<dim; idim++) orthAB[idim] = AC[idim]-AB_dot_AC*AB[idim];
1172 double denom= normAC+AB_dot_AC;
1173 double numer=norm<dim>(orthAB);
1175 return 2*atan2(numer,denom);
1178 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
1179 /* Tells whether the frame constituted of vectors AB, AC and n is direct */
1180 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
1181 template<int dim> inline double direct_frame(const double * A, const double * B, const double * C, double * n);
1183 double direct_frame<2>(const double * A, const double * B, const double * C, double * n)
1187 for(int idim =0; idim<2; idim++) AB[idim] = B[idim]-A[idim];//B-A;
1188 for(int idim =0; idim<2; idim++) AC[idim] = C[idim]-A[idim];//C-A;
1190 return determinant(AB,AC)*n[0];
1193 double direct_frame<3>(const double * A, const double * B, const double * C, double * n)
1197 for(int idim =0; idim<3; idim++) AB[idim] = B[idim]-A[idim];//B-A;
1198 for(int idim =0; idim<3; idim++) AC[idim] = C[idim]-A[idim];//C-A;
1200 return determinant(AB,AC,n)>0;
1203 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
1204 /* calcul l'intersection de deux polygones COPLANAIRES */
1205 /* en dimension DIM (2 ou 3). Si DIM=3 l'algorithme ne considere*/
1206 /* que les deux premieres coordonnees de chaque point */
1207 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
1208 template<int DIM> inline void intersec_de_polygone(const double * Coords_A, const double * Coords_B,
1209 int nb_NodesA, int nb_NodesB,
1210 std::vector<double>& inter,
1211 double dim_caracteristic, double precision)
1213 for(int i_A = 1; i_A<nb_NodesA-1; i_A++)
1215 for(int i_B = 1; i_B<nb_NodesB-1; i_B++)
1217 INTERP_KERNEL::intersec_de_triangle(&Coords_A[0],&Coords_A[DIM*i_A],&Coords_A[DIM*(i_A+1)],
1218 &Coords_B[0],&Coords_B[DIM*i_B],&Coords_B[DIM*(i_B+1)],
1219 inter, dim_caracteristic, precision);
1222 int nb_inter=((int)inter.size())/DIM;
1223 if(nb_inter >3) inter=INTERP_KERNEL::reconstruct_polygon(inter);
1227 *_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
1228 * fonctions qui calcule l'aire d'un polygone en dimension 2 ou 3
1229 *_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
1230 template<int DIM> inline double polygon_area(std::vector<double>& inter)
1235 for(int i = 1; i<(int)inter.size()/DIM-1; i++)
1237 INTERP_KERNEL::crossprod<DIM>(&inter[0],&inter[DIM*i],&inter[DIM*(i+1)],area);
1238 result +=0.5*norm<DIM>(area);
1243 template<int DIM> inline double polygon_area(std::deque<double>& inter)
1248 for(int i = 1; i<(int)inter.size()/DIM-1; i++)
1250 INTERP_KERNEL::crossprod<DIM>(&inter[0],&inter[DIM*i],&inter[DIM*(i+1)],area);
1251 result +=0.5*norm<DIM>(area);
1256 /*! Computes the triple product (XA^XB).XC (in 3D)*/
1257 inline double triple_product(const double* A, const double*B, const double*C, const double*X)
1273 (XA[1]*XB[2]-XA[2]*XB[1])*XC[0]+
1274 (XA[2]*XB[0]-XA[0]*XB[2])*XC[1]+
1275 (XA[0]*XB[1]-XA[1]*XB[0])*XC[2];
1278 /*! Subroutine of checkEqualPolygins that tests if two list of nodes (not necessarily distincts) describe the same polygon, assuming they share a common point.*/
1279 /*! Indexes istart1 and istart2 designate two points P1 in L1 and P2 in L2 that have identical coordinates. Generally called with istart1=0.*/
1280 /*! Integer sign ( 1 or -1) indicate the direction used in going all over L2. */
1281 template<class T, int dim>
1282 bool checkEqualPolygonsOneDirection(T * L1, T * L2, int size1, int size2, int istart1, int istart2, double epsilon, int sign)
1286 int i1next = ( i1 + 1 ) % size1;
1287 int i2next = ( i2 + sign +size2) % size2;
1291 while( i1next!=istart1 && distance2<T,dim>(L1,i1*dim, L1,i1next*dim) < epsilon ) i1next = ( i1next + 1 ) % size1;
1292 while( i2next!=istart2 && distance2<T,dim>(L2,i2*dim, L2,i2next*dim) < epsilon ) i2next = ( i2next + sign +size2 ) % size2;
1294 if(i1next == istart1)
1296 if(i2next == istart2)
1301 if(i2next == istart2)
1305 if(distance2<T,dim>(L1,i1next*dim, L2,i2next*dim) > epsilon )
1311 i1next = ( i1 + 1 ) % size1;
1312 i2next = ( i2 + sign + size2 ) % size2;
1318 /*! Tests if two list of nodes (not necessarily distincts) describe the same polygon.*/
1319 /*! Existence of multiple points in the list is considered.*/
1320 template<class T, int dim>
1321 bool checkEqualPolygons(T * L1, T * L2, double epsilon)
1323 if(L1==NULL || L2==NULL)
1325 std::cout << "Warning InterpolationUtils.hxx:checkEqualPolygonsPointer: Null pointer " << std::endl;
1326 throw(Exception("big error: not closed polygon..."));
1329 int size1 = (*L1).size()/dim;
1330 int size2 = (*L2).size()/dim;
1334 while( istart2 < size2 && distance2<T,dim>(L1,istart1*dim, L2,istart2*dim) > epsilon ) istart2++;
1336 if(istart2 == size2)
1338 return (size1 == 0) && (size2 == 0);
1341 return checkEqualPolygonsOneDirection<T,dim>( L1, L2, size1, size2, istart1, istart2, epsilon, 1)
1342 || checkEqualPolygonsOneDirection<T,dim>( L1, L2, size1, size2, istart1, istart2, epsilon, -1);