1 // Copyright (C) 2007-2008 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.
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 #ifndef __INTERPOLATIONUTILS_HXX__
20 #define __INTERPOLATIONUTILS_HXX__
22 #include "INTERPKERNELDefines.hxx"
23 #include "InterpKernelException.hxx"
25 #include "NormalizedUnstructuredMesh.hxx"
36 namespace INTERP_KERNEL
38 template<class ConnType, NumberingPolicy numPol>
39 class OTT//OffsetToolTrait
43 template<class ConnType>
44 class OTT<ConnType,ALL_C_MODE>
47 static ConnType indFC(ConnType i) { return i; }
48 static ConnType ind2C(ConnType i) { return i; }
49 static ConnType conn2C(ConnType i) { return i; }
50 static ConnType coo2C(ConnType i) { return i; }
53 template<class ConnType>
54 class OTT<ConnType,ALL_FORTRAN_MODE>
57 static ConnType indFC(ConnType i) { return i+1; }
58 static ConnType ind2C(ConnType i) { return i-1; }
59 static ConnType conn2C(ConnType i) { return i-1; }
60 static ConnType coo2C(ConnType i) { return i-1; }
63 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
64 /* calcul la surface d'un triangle */
65 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
67 inline double Surf_Tri(const double* P_1,const double* P_2,const double* P_3)
69 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]);
70 double Surface = 0.5*fabs(A);
74 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
75 /* fonction qui calcul le déterminant */
76 /* de deux vecteur(cf doc CGAL). */
77 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
79 //fonction qui calcul le déterminant des vecteurs: P3P1 et P3P2
82 inline double mon_determinant(const double* P_1,
86 double mon_det=(P_1[0]-P_3[0])*(P_2[1]-P_3[1])-(P_2[0]-P_3[0])*(P_1[1]-P_3[1]);
90 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
91 //calcul la norme du vecteur P1P2
93 inline double norme_vecteur(const double* P_1,const double* P_2)
95 double X=P_1[0]-P_2[0];
96 double Y=P_1[1]-P_2[1];
97 double norme=sqrt(X*X+Y*Y);
101 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
102 /* calcul le cos et le sin de l'angle P1P2,P1P3 */
103 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
105 inline std::vector<double> calcul_cos_et_sin(const double* P_1,
110 std::vector<double> Vect;
111 double P1_P2=norme_vecteur(P_1,P_2);
112 double P2_P3=norme_vecteur(P_2,P_3);
113 double P3_P1=norme_vecteur(P_3,P_1);
115 double N=P1_P2*P1_P2+P3_P1*P3_P1-P2_P3*P2_P3;
116 double D=2.0*P1_P2*P3_P1;
118 if (COS>1.0) COS=1.0;
119 if (COS<-1.0) COS=-1.0;
121 double V=mon_determinant(P_2,P_3,P_1);
122 double D_1=P1_P2*P3_P1;
124 if (SIN>1.0) SIN=1.0;
125 if (SIN<-1.0) SIN=-1.0;
133 * This method builds a quadrangle built with the first point of 'triIn' the barycenter of two edges starting or ending with
134 * the first point of 'triIn' and the barycenter of 'triIn'.
136 * @param triIn is a 6 doubles array in full interlace mode, that represents a triangle.
137 * @param quadOut is a 8 doubles array filled after the following call.
139 template<int SPACEDIM>
140 inline void fillDualCellOfTri(const double *triIn, double *quadOut)
143 std::copy(triIn,triIn+SPACEDIM,quadOut);
144 double tmp[SPACEDIM];
145 std::transform(triIn,triIn+SPACEDIM,triIn+SPACEDIM,tmp,std::plus<double>());
147 std::transform(tmp,tmp+SPACEDIM,quadOut+SPACEDIM,std::bind2nd(std::multiplies<double>(),0.5));
148 std::transform(tmp,tmp+SPACEDIM,triIn+2*SPACEDIM,tmp,std::plus<double>());
150 std::transform(tmp,tmp+SPACEDIM,quadOut+2*SPACEDIM,std::bind2nd(std::multiplies<double>(),1/3.));
152 std::transform(triIn,triIn+SPACEDIM,triIn+2*SPACEDIM,tmp,std::plus<double>());
153 std::transform(tmp,tmp+SPACEDIM,quadOut+3*SPACEDIM,std::bind2nd(std::multiplies<double>(),0.5));
156 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
157 /* calcul les coordonnées du barycentre d'un polygone */
158 /* le vecteur en entrée est constitué des coordonnées */
159 /* des sommets du polygone */
160 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
162 inline std::vector<double> bary_poly(const std::vector<double>& V)
164 std::vector<double> Bary;
165 long taille=V.size();
169 for(long i=0;i<taille/2;i++)
176 Bary.push_back(A);//taille vecteur=2*nb de points.
184 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
185 /* Calculate barycentric coordinates of a 2D point p */
186 /* with respect to the triangle verices. */
187 /* triaCoords are in full interlace */
188 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
190 inline void barycentric_coords(const double* triaCoords, const double* p, double* bc)
194 T11 = triaCoords[0]-triaCoords[4], T12 = triaCoords[2]-triaCoords[4],
195 T21 = triaCoords[1]-triaCoords[5], T22 = triaCoords[3]-triaCoords[5];
196 // matrix determinant
197 double Tdet = T11*T22 - T12*T21;
198 if ( fabs( Tdet ) < std::numeric_limits<double>::min() ) {
199 bc[0]=1; bc[1]=0; bc[2]=0;
203 double t11 = T22, t12 = -T12, t21 = -T21, t22 = T11;
205 double r11 = p[0]-triaCoords[4], r12 = p[1]-triaCoords[5];
206 // barycentric coordinates: mutiply matrix by vector
207 bc[0] = (t11 * r11 + t12 * r12)/Tdet;
208 bc[1] = (t21 * r11 + t22 * r12)/Tdet;
209 bc[2] = 1. - bc[0] - bc[1];
212 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
213 /* Calculate barycentric coordinates of a point p */
214 /* with respect to triangle or tetra verices. */
215 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
217 inline void barycentric_coords(const std::vector<const double*>& n, const double* p, double* bc)
220 if ( n.size() == 3 ) // TRIA3
224 T11 = n[0][_X]-n[2][_X], T12 = n[1][_X]-n[2][_X],
225 T21 = n[0][_Y]-n[2][_Y], T22 = n[1][_Y]-n[2][_Y];
226 // matrix determinant
227 double Tdet = T11*T22 - T12*T21;
228 if ( std::fabs( Tdet ) < std::numeric_limits<double>::min() ) {
229 bc[0]=1; bc[1]=bc[2]=0; // no solution
233 double t11 = T22, t12 = -T12, t21 = -T21, t22 = T11;
235 double r11 = p[_X]-n[2][_X], r12 = p[_Y]-n[2][_Y];
236 // barycentric coordinates: mutiply matrix by vector
237 bc[0] = (t11 * r11 + t12 * r12)/Tdet;
238 bc[1] = (t21 * r11 + t22 * r12)/Tdet;
239 bc[2] = 1. - bc[0] - bc[1];
243 bc[3]=0; // for no solution
245 // Find bc by solving system of 3 equations using Gaussian elimination algorithm
246 // bc1*( x1 - x4 ) + bc2*( x2 - x4 ) + bc3*( x3 - x4 ) = px - x4
247 // bc1*( y1 - y4 ) + bc2*( y2 - y4 ) + bc3*( y3 - y4 ) = px - y4
248 // bc1*( z1 - z4 ) + bc2*( z2 - z4 ) + bc3*( z3 - z4 ) = px - z4
249 const int nbCol=4, nbRow=3;
251 double T[nbRow][nbCol]=
252 {{ n[0][_X]-n[3][_X], n[1][_X]-n[3][_X], n[2][_X]-n[3][_X], p[_X]-n[3][_X] },
253 { n[0][_Y]-n[3][_Y], n[1][_Y]-n[3][_Y], n[2][_Y]-n[3][_Y], p[_Y]-n[3][_Y] },
254 { n[0][_Z]-n[3][_Z], n[1][_Z]-n[3][_Z], n[2][_Z]-n[3][_Z], p[_Z]-n[3][_Z] }};
256 // make upper triangular matrix (forward elimination)
258 int iR[nbRow] = { 0, 1, 2 };
260 for ( int i = 0; i < 2; ++i ) // nullify 2 rows
262 // swap rows to have max value of i-th column in i-th row
263 double max = std::fabs( T[ iR[i] ][i] );
264 for ( int r = i+1; r < nbRow; ++r ) {
265 double t = std::fabs( T[ iR[r] ][i] );
268 std::swap( iR[r], iR[i] );
271 if ( max < std::numeric_limits<double>::min() ) {
272 bc[0]=1; bc[1]=bc[2]=bc[3]=0;
273 return; // no solution
275 // make 0 below T[i][i] (actually we do not modify i-th column)
276 double* tUpRow = T[ iR[i] ];
277 for ( int r = i+1; r < nbRow; ++r ) {
278 double* tRow = T[ iR[r] ];
279 double coef = tRow[ i ] / tUpRow[ i ];
280 for ( int c = i+1; c < nbCol; ++c )
281 tRow[ c ] -= tUpRow[ c ] * coef;
284 double* tRow = T[ iR[2] ];
285 if ( std::fabs( tRow[ 2 ] ) < std::numeric_limits<double>::min() ) {
286 bc[0]=1; bc[1]=bc[2]=bc[3]=0;
287 return; // no solution
289 tRow[ 3 ] /= tRow[ 2 ];
291 // calculate solution (back substitution)
296 bc[ 1 ] = (tRow[ 3 ] - bc[2]*tRow[ 2 ]) / tRow[ 1 ];
299 bc[ 0 ] = (tRow[ 3 ] - bc[2]*tRow[ 2 ] - bc[1]*tRow[ 1 ]) / tRow[ 0 ];
301 bc[ 3 ] = 1. - bc[0] - bc[1] - bc[2];
305 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
306 /* calcul la surface d'un polygone. */
307 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
309 inline double Surf_Poly(const std::vector<double>& Poly)
313 for(unsigned long i=0; i<(Poly.size())/2-2; i++)
315 double Surf=Surf_Tri( &Poly[0],&Poly[2*(i+1)],&Poly[2*(i+2)] );
316 Surface=Surface + Surf ;
321 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
322 /* fonction qui teste si un point est dans une maille */
324 /* P_1, P_2, P_3 sommet des mailles */
325 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
327 inline bool point_dans_triangle(const double* P_0,const double* P_1,
328 const double* P_2,const double* P_3,
333 double det_1=mon_determinant(P_1,P_3,P_0);
334 double det_2=mon_determinant(P_3,P_2,P_0);
335 double det_3=mon_determinant(P_2,P_1,P_0);
336 if( (det_1>=-eps && det_2>=-eps && det_3>=-eps) || (det_1<=eps && det_2<=eps && det_3<=eps) )
344 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
345 /*fonction pour vérifier qu'un point n'a pas déja été considérer dans */
346 /* le vecteur et le rajouter au vecteur sinon. */
347 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
349 inline void verif_point_dans_vect(const double* P, std::vector<double>& V, double absolute_precision )
351 long taille=V.size();
352 bool isPresent=false;
353 for(long i=0;i<taille/2;i++)
355 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)
367 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
368 /* fonction qui rajoute les sommet du triangle P dans le vecteur V */
369 /* si ceux-ci sont compris dans le triangle S et ne sont pas déjà dans */
371 /*sommets de P: P_1, P_2, P_3 */
372 /*sommets de S: P_4, P_5, P_6 */
373 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
375 inline void rajou_sommet_triangl(const double* P_1,const double* P_2,const double* P_3,
376 const double* P_4,const double* P_5,const double* P_6,
377 std::vector<double>& V, double dim_caracteristic, double precision)
380 double absolute_precision = precision*dim_caracteristic;
381 bool A_1=INTERP_KERNEL::point_dans_triangle(P_1,P_4,P_5,P_6,absolute_precision);
383 verif_point_dans_vect(P_1,V,absolute_precision);
384 bool A_2=INTERP_KERNEL::point_dans_triangle(P_2,P_4,P_5,P_6,absolute_precision);
386 verif_point_dans_vect(P_2,V,absolute_precision);
387 bool A_3=INTERP_KERNEL::point_dans_triangle(P_3,P_4,P_5,P_6,absolute_precision);
389 verif_point_dans_vect(P_3,V,absolute_precision);
393 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
394 /* calcul de l'intersection de deux segments: segments P1P2 avec P3P4 */
395 /* . Si l'intersection est non nulle et si celle-ci n'est */
396 /* n'est pas déjà contenue dans Vect on la rajoute à Vect */
397 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
399 inline void inters_de_segment(const double * P_1,const double * P_2,
400 const double * P_3,const double * P_4,
401 std::vector<double>& Vect,
402 double dim_caracteristic, double precision)
404 // calcul du déterminant de P_1P_2 et P_3P_4.
405 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]);
407 double absolute_precision = dim_caracteristic*precision;
408 if(fabs(det)>absolute_precision)
410 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;
412 if (k_1 >= -absolute_precision && k_1 <= 1+absolute_precision)
413 //if( k_1 >= -precision && k_1 <= 1+precision)
415 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;
417 if (k_2 >= -absolute_precision && k_2 <= 1+absolute_precision)
418 //if( k_2 >= -precision && k_2 <= 1+precision)
421 P_0[0]=P_1[0]+k_1*(P_2[0]-P_1[0]);
422 P_0[1]=P_1[1]+k_1*(P_2[1]-P_1[1]);
423 verif_point_dans_vect(P_0,Vect,absolute_precision);
431 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
432 /* calcul l'intersection de deux triangles */
433 /* P_1, P_2, P_3: sommets du premier triangle */
434 /* P_4, P_5, P_6: sommets du deuxième triangle */
435 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
437 inline void intersec_de_triangle(const double* P_1,const double* P_2, const double* P_3,
438 const double* P_4,const double* P_5,const double* P_6,
439 std::vector<double>& Vect, double dim_caracteristic, double precision)
441 inters_de_segment(P_1,P_2,P_4,P_5,Vect, dim_caracteristic, precision);
442 inters_de_segment(P_1,P_2,P_5,P_6,Vect, dim_caracteristic, precision);
443 inters_de_segment(P_1,P_2,P_6,P_4,Vect, dim_caracteristic, precision);
444 inters_de_segment(P_2,P_3,P_4,P_5,Vect, dim_caracteristic, precision);
445 inters_de_segment(P_2,P_3,P_5,P_6,Vect, dim_caracteristic, precision);
446 inters_de_segment(P_2,P_3,P_6,P_4,Vect, dim_caracteristic, precision);
447 inters_de_segment(P_3,P_1,P_4,P_5,Vect, dim_caracteristic, precision);
448 inters_de_segment(P_3,P_1,P_5,P_6,Vect, dim_caracteristic, precision);
449 inters_de_segment(P_3,P_1,P_6,P_4,Vect, dim_caracteristic, precision);
450 rajou_sommet_triangl(P_1,P_2,P_3,P_4,P_5,P_6,Vect, dim_caracteristic, precision);
451 rajou_sommet_triangl(P_4,P_5,P_6,P_1,P_2,P_3,Vect, dim_caracteristic, precision);
454 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
455 /* fonction pour vérifier qu'un n°de maille n'a pas déja été considérer */
456 /* dans le vecteur et le rajouter au vecteur sinon. */
457 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
459 inline void verif_maill_dans_vect(int Num, std::vector<int>& V)
461 long taille=V.size();
463 for(long i=0;i<taille;i++)
475 /*! Function that compares two angles from the values of the pairs (sin,cos)*/
476 /*! Angles are considered in [0, 2Pi] bt are not computed explicitely */
480 bool operator()(std::pair<double,double>theta1, std::pair<double,double> theta2)
482 double norm1 = sqrt(theta1.first*theta1.first +theta1.second*theta1.second);
483 double norm2 = sqrt(theta2.first*theta2.first +theta2.second*theta2.second);
485 double epsilon = 1.e-12;
487 if( norm1 < epsilon || norm2 < epsilon )
488 std::cout << "Warning InterpolationUtils.hxx: AngleLess : Vector with zero norm, cannot define the angle !!!! " << std::endl;
490 return theta1.second*(norm2 + theta2.first) < theta2.second*(norm1 + theta1.first);
496 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
497 /* fonction pour reconstituer un polygone convexe à partir */
498 /* d'un nuage de point. */
499 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
501 inline std::vector<double> reconstruct_polygon(const std::vector<double>& V)
512 double *COS=new double[taille/2];
513 double *SIN=new double[taille/2];
514 //double *angle=new double[taille/2];
515 std::vector<double> Bary=bary_poly(V);
519 for(int i=0; i<taille/2-1;i++)
521 std::vector<double> Trigo=calcul_cos_et_sin(&Bary[0],&V[0],&V[2*(i+1)]);
525 // {angle[i+1]=atan2(SIN[i+1],COS[i+1]);}
527 // {angle[i+1]=-atan2(SIN[i+1],COS[i+1]);}
530 //ensuite on ordonne les angles.
531 std::vector<double> Pt_ordonne;
532 Pt_ordonne.reserve(taille);
533 // std::multimap<double,int> Ordre;
534 std::multimap<std::pair<double,double>,int, AngleLess> CosSin;
535 for(int i=0;i<taille/2;i++)
537 // Ordre.insert(std::make_pair(angle[i],i));
538 CosSin.insert(std::make_pair(std::make_pair(SIN[i],COS[i]),i));
540 // std::multimap <double,int>::iterator mi;
541 std::multimap<std::pair<double,double>,int, AngleLess>::iterator micossin;
542 // for(mi=Ordre.begin();mi!=Ordre.end();mi++)
544 // int j=(*mi).second;
545 // Pt_ordonne.push_back(V[2*j]);
546 // Pt_ordonne.push_back(V[2*j+1]);
548 for(micossin=CosSin.begin();micossin!=CosSin.end();micossin++)
550 int j=(*micossin).second;
551 Pt_ordonne.push_back(V[2*j]);
552 Pt_ordonne.push_back(V[2*j+1]);
561 template<int DIM, NumberingPolicy numPol, class MyMeshType>
562 inline void getElemBB(double* bb, const double *coordsOfMesh, int iP, int nb_nodes)
564 bb[0]=std::numeric_limits<double>::max();
565 bb[1]=-std::numeric_limits<double>::max();
566 bb[2]=std::numeric_limits<double>::max();
567 bb[3]=-std::numeric_limits<double>::max();
568 bb[4]=std::numeric_limits<double>::max();
569 bb[5]=-std::numeric_limits<double>::max();
571 for (int i=0; i<nb_nodes; i++)
573 double x = coordsOfMesh[3*(iP+i)];
574 double y = coordsOfMesh[3*(iP+i)+1];
575 double z = coordsOfMesh[3*(iP+i)+2];
576 bb[0]=(x<bb[0])?x:bb[0];
577 bb[1]=(x>bb[1])?x:bb[1];
578 bb[2]=(y<bb[2])?y:bb[2];
579 bb[3]=(y>bb[3])?y:bb[3];
580 bb[4]=(z<bb[4])?z:bb[4];
581 bb[5]=(z>bb[5])?z:bb[5];
585 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
586 /* Computes the dot product of a and b */
587 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
589 inline double dotprod( double * a, double * b)
592 for(int idim = 0; idim < dim ; idim++) result += a[idim]*b[idim];
595 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
596 /* Computes the norm of vector v */
597 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
599 inline double norm( double * v)
602 for(int idim =0; idim<dim; idim++) result+=v[idim]*v[idim];
605 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
606 /* Computes the square norm of vector a-b */
607 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
609 inline double distance2( const double * a, const double * b)
612 for(int idim =0; idim<dim; idim++) result+=(a[idim]-b[idim])*(a[idim]-b[idim]);
615 template<class T, int dim>
616 inline double distance2( T * a, int inda, T * b, int indb)
619 for(int idim =0; idim<dim; idim++) result += ((*a)[inda+idim] - (*b)[indb+idim])* ((*a)[inda+idim] - (*b)[indb+idim]);
622 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
623 /* Computes the determinant of a and b */
624 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
625 inline double determinant ( double * a, double * b)
627 return a[0]*b[1]-a[1]*b[0];
629 inline double determinant ( double * a, double * b, double * c)
631 return a[0]*determinant(b+1,c+1)-b[0]*determinant(a+1,c+1)+c[0]*determinant(a+1,b+1);
634 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
635 /* Computes the cross product of AB and AC */
636 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
638 template<int dim> inline void crossprod( const double * A, const double * B, const double * C, double * V);
641 void crossprod<2>( const double * A, const double * B, const double * C, double * V)
645 for(int idim =0; idim<2; idim++) AB[idim] = B[idim]-A[idim];//B-A
646 for(int idim =0; idim<2; idim++) AC[idim] = C[idim]-A[idim];//C-A;
648 V[0]=determinant(AB,AC);
652 void crossprod<3>( const double * A, const double * B, const double * C, double * V)
656 for(int idim =0; idim<3; idim++) AB[idim] = B[idim]-A[idim];//B-A
657 for(int idim =0; idim<3; idim++) AC[idim] = C[idim]-A[idim];//C-A;
659 V[0]=AB[1]*AC[2]-AB[2]*AC[1];
660 V[1]=-AB[0]*AC[2]+AB[2]*AC[0];
661 V[2]=AB[0]*AC[1]-AB[1]*AC[0];
664 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
665 /* Checks wether point A is inside the quadrangle BCDE */
666 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
668 template<int dim> inline double check_inside(const double* A,const double* B,const double* C,const double* D,
669 const double* E,double* ABC, double* ADE)
671 crossprod<dim>(A,B,C,ABC);
672 crossprod<dim>(A,D,E,ADE);
673 return dotprod<dim>(ABC,ADE);
677 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
678 /* Computes the geometric angle (in [0,Pi]) between two non zero vectors AB and AC */
679 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
680 template<int dim> inline double angle(const double * A, const double * B, const double * C, double * n)
686 for(int idim =0; idim<dim; idim++) AB[idim] = B[idim]-A[idim];//B-A;
687 for(int idim =0; idim<dim; idim++) AC[idim] = C[idim]-A[idim];//C-A;
689 double normAB= norm<dim>(AB);
690 for(int idim =0; idim<dim; idim++) AB[idim]/=normAB;
692 double normAC= norm<dim>(AC);
693 double AB_dot_AC=dotprod<dim>(AB,AC);
694 for(int idim =0; idim<dim; idim++) orthAB[idim] = AC[idim]-AB_dot_AC*AB[idim];
696 double denom= normAC+AB_dot_AC;
697 double numer=norm<dim>(orthAB);
699 return 2*atan2(numer,denom);
702 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
703 /* Tells whether the frame constituted of vectors AB, AC and n is direct */
704 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
705 template<int dim> inline double direct_frame(const double * A, const double * B, const double * C, double * n);
707 double direct_frame<2>(const double * A, const double * B, const double * C, double * n)
711 for(int idim =0; idim<2; idim++) AB[idim] = B[idim]-A[idim];//B-A;
712 for(int idim =0; idim<2; idim++) AC[idim] = C[idim]-A[idim];//C-A;
714 return determinant(AB,AC)*n[0];
717 double direct_frame<3>(const double * A, const double * B, const double * C, double * n)
721 for(int idim =0; idim<3; idim++) AB[idim] = B[idim]-A[idim];//B-A;
722 for(int idim =0; idim<3; idim++) AC[idim] = C[idim]-A[idim];//C-A;
724 return determinant(AB,AC,n)>0;
727 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
728 /* calcul l'intersection de deux polygones COPLANAIRES */
729 /* en dimension DIM (2 ou 3). Si DIM=3 l'algorithme ne considère*/
730 /* que les deux premières coordonnées de chaque point */
731 /*_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _*/
732 template<int DIM> inline void intersec_de_polygone(const double * Coords_A, const double * Coords_B,
733 int nb_NodesA, int nb_NodesB,
734 std::vector<double>& inter,
735 double dim_caracteristic, double precision)
737 for(int i_A = 1; i_A<nb_NodesA-1; i_A++)
739 for(int i_B = 1; i_B<nb_NodesB-1; i_B++)
741 INTERP_KERNEL::intersec_de_triangle(&Coords_A[0],&Coords_A[DIM*i_A],&Coords_A[DIM*(i_A+1)],
742 &Coords_B[0],&Coords_B[DIM*i_B],&Coords_B[DIM*(i_B+1)],
743 inter, dim_caracteristic, precision);
746 int nb_inter=((int)inter.size())/DIM;
747 if(nb_inter >3) inter=INTERP_KERNEL::reconstruct_polygon(inter);
751 *_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _
752 * fonctions qui calcule l'aire d'un polygone en dimension 2 ou 3
753 *_ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ */
754 template<int DIM> inline double polygon_area(std::vector<double>& inter)
759 for(int i = 1; i<(int)inter.size()/DIM-1; i++)
761 INTERP_KERNEL::crossprod<DIM>(&inter[0],&inter[DIM*i],&inter[DIM*(i+1)],area);
762 result +=0.5*norm<DIM>(area);
767 template<int DIM> inline double polygon_area(std::deque<double>& inter)
772 for(int i = 1; i<(int)inter.size()/DIM-1; i++)
774 INTERP_KERNEL::crossprod<DIM>(&inter[0],&inter[DIM*i],&inter[DIM*(i+1)],area);
775 result +=0.5*norm<DIM>(area);
780 /*! Computes the triple product (XA^XB).XC (in 3D)*/
781 inline double triple_product(const double* A, const double*B, const double*C, const double*X)
797 (XA[1]*XB[2]-XA[2]*XB[1])*XC[0]+
798 (XA[2]*XB[0]-XA[0]*XB[2])*XC[1]+
799 (XA[0]*XB[1]-XA[1]*XB[0])*XC[2];
802 /*! Subroutine of checkEqualPolygins that tests if two list of nodes (not necessarily distincts) describe the same polygon, assuming they share a comon point.*/
803 /*! Indexes istart1 and istart2 designate two points P1 in L1 and P2 in L2 that have identical coordinates. Generally called with istart1=0.*/
804 /*! Integer sign ( 1 or -1) indicate the direction used in going all over L2. */
805 template<class T, int dim>
806 bool checkEqualPolygonsOneDirection(T * L1, T * L2, int size1, int size2, int istart1, int istart2, double epsilon, int sign)
810 int i1next = ( i1 + 1 ) % size1;
811 int i2next = ( i2 + sign +size2) % size2;
815 while( i1next!=istart1 && distance2<T,dim>(L1,i1*dim, L1,i1next*dim) < epsilon ) i1next = ( i1next + 1 ) % size1;
816 while( i2next!=istart2 && distance2<T,dim>(L2,i2*dim, L2,i2next*dim) < epsilon ) i2next = ( i2next + sign +size2 ) % size2;
818 if(i1next == istart1)
820 if(i2next == istart2)
825 if(i2next == istart2)
829 if(distance2<T,dim>(L1,i1next*dim, L2,i2next*dim) > epsilon )
835 i1next = ( i1 + 1 ) % size1;
836 i2next = ( i2 + sign + size2 ) % size2;
842 /*! Tests if two list of nodes (not necessarily distincts) describe the same polygon.*/
843 /*! Existence of multiple points in the list is considered.*/
844 template<class T, int dim>
845 bool checkEqualPolygons(T * L1, T * L2, double epsilon)
847 if(L1==NULL || L2==NULL)
849 std::cout << "Warning InterpolationUtils.hxx:checkEqualPolygonsPointer: Null pointer " << std::endl;
850 throw(Exception("big error: not closed polygon..."));
853 int size1 = (*L1).size()/dim;
854 int size2 = (*L2).size()/dim;
858 while( istart2 < size2 && distance2<T,dim>(L1,istart1*dim, L2,istart2*dim) > epsilon ) istart2++;
862 return (size1 == 0) && (size2 == 0);
865 return checkEqualPolygonsOneDirection<T,dim>( L1, L2, size1, size2, istart1, istart2, epsilon, 1)
866 || checkEqualPolygonsOneDirection<T,dim>( L1, L2, size1, size2, istart1, istart2, epsilon, -1);