ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/mailleur3d_intersection.cpp
Revision: 553
Committed: Fri Sep 26 20:09:41 2014 UTC (10 years, 7 months ago) by francois
Original Path: magic/lib/mailleur_auto/src/mailleur3d_intersection.cpp
File size: 28893 byte(s)
Log Message:
nouvelle intersection dans le mailleur3d avec meilleur controle de l'erreur numerique

File Contents

# User Rev Content
1 5 //------------------------------------------------------------
2     //------------------------------------------------------------
3     // MAGiC
4     // Jean Christophe Cuillière et Vincent FRANCOIS
5     // Département de Génie Mécanique - UQTR
6     //------------------------------------------------------------
7     // Le projet MAGIC est un projet de recherche du département
8     // de génie mécanique de l'Université du Québec à
9     // Trois Rivières
10     // Les librairies ne peuvent être utilisées sans l'accord
11     // des auteurs (contact : francois@uqtr.ca)
12     //------------------------------------------------------------
13     //------------------------------------------------------------
14     //
15     // mailleur3d_intersection.cpp
16     //
17     //------------------------------------------------------------
18     //------------------------------------------------------------
19     // COPYRIGHT 2000
20     // Version du 02/03/2006 à 11H23
21     //------------------------------------------------------------
22     //------------------------------------------------------------
23    
24    
25     #include "gestionversion.h"
26    
27    
28     //#pragma hdrstop
29    
30     #include <math.h>
31     #include "mailleur3d.h"
32     #include "ot_mathematique.h"
33 francois 553 #include "ot_doubleprecision.h"
34 5
35    
36    
37 francois 553
38 5 int MAILLEUR3D::inter_tetra_triangle(MG_NOEUD* noeud1,MG_NOEUD* noeud2,MG_NOEUD* noeud3,MG_NOEUD* noeud4,MG_NOEUD* noeud5,MG_NOEUD* noeud6,MG_NOEUD* noeud7)
39     {
40     // filtre sur numero
41     int commun1;
42     int commun2;
43     int commun3;
44     if ( (noeud5==noeud1) || (noeud5==noeud2) || (noeud5==noeud3) || (noeud5==noeud4) ) commun1=1; else commun1=0;
45     if ( (noeud6==noeud1) || (noeud6==noeud2) || (noeud6==noeud3) || (noeud6==noeud4) ) commun2=1; else commun2=0;
46     if ( (noeud7==noeud1) || (noeud7==noeud2) || (noeud7==noeud3) || (noeud7==noeud4) ) commun3=1; else commun3=0;
47     int nb_commun=commun1+commun2+commun3;
48     if (nb_commun==3) return(false);
49     // filtre sur position
50     int inter;
51     if (nb_commun>0) inter=true;
52     else
53     {
54     double bxmin=std::min(noeud1->get_x(),noeud2->get_x());
55     bxmin=std::min(bxmin,noeud3->get_x());
56     bxmin=std::min(bxmin,noeud4->get_x());
57     double bxmax=std::max(noeud1->get_x(),noeud2->get_x());
58     bxmax=std::max(bxmax,noeud3->get_x());
59     bxmax=std::max(bxmax,noeud4->get_x());
60     double bymin=std::min(noeud1->get_y(),noeud2->get_y());
61     bymin=std::min(bymin,noeud3->get_y());
62     bymin=std::min(bymin,noeud4->get_y());
63     double bymax=std::max(noeud1->get_y(),noeud2->get_y());
64     bymax=std::max(bymax,noeud3->get_y());
65     bymax=std::max(bymax,noeud4->get_y());
66     double bzmin=std::min(noeud1->get_z(),noeud2->get_z());
67     bzmin=std::min(bzmin,noeud3->get_z());
68     bzmin=std::min(bzmin,noeud4->get_z());
69     double bzmax=std::max(noeud1->get_z(),noeud2->get_z());
70     bzmax=std::max(bzmax,noeud3->get_z());
71     bzmax=std::max(bzmax,noeud4->get_z());
72     BOITE_3D boite_tetra(bxmin,bymin,bzmin,bxmax,bymax,bzmax);
73     double b2xmin=std::min(noeud5->get_x(),noeud6->get_x());
74     b2xmin=std::min(b2xmin,noeud7->get_x());
75     double b2xmax=std::max(noeud5->get_x(),noeud6->get_x());
76     b2xmax=std::max(b2xmax,noeud7->get_x());
77     double b2ymin=std::min(noeud5->get_y(),noeud6->get_y());
78     b2ymin=std::min(b2ymin,noeud7->get_y());
79     double b2ymax=std::max(noeud5->get_y(),noeud6->get_y());
80     b2ymax=std::max(b2ymax,noeud7->get_y());
81     double b2zmin=std::min(noeud5->get_z(),noeud6->get_z());
82     b2zmin=std::min(b2zmin,noeud7->get_z());
83     double b2zmax=std::max(noeud5->get_z(),noeud6->get_z());
84     b2zmax=std::max(b2zmax,noeud7->get_z());
85     BOITE_3D boite_triangle(b2xmin,b2ymin,b2zmin,b2xmax,b2ymax,b2zmax);
86     inter=boite_tetra*boite_triangle;
87     }
88     if (inter==true)
89     {
90     inter=false;
91     OT_VECTEUR_3D n1n2(noeud2->get_x()-noeud1->get_x(),noeud2->get_y()-noeud1->get_y(),noeud2->get_z()-noeud1->get_z());
92     OT_VECTEUR_3D n1n3(noeud3->get_x()-noeud1->get_x(),noeud3->get_y()-noeud1->get_y(),noeud3->get_z()-noeud1->get_z());
93     OT_VECTEUR_3D n1n4(noeud4->get_x()-noeud1->get_x(),noeud4->get_y()-noeud1->get_y(),noeud4->get_z()-noeud1->get_z());
94     OT_VECTEUR_3D n2n3(noeud3->get_x()-noeud2->get_x(),noeud3->get_y()-noeud2->get_y(),noeud3->get_z()-noeud2->get_z());
95     OT_VECTEUR_3D n2n4(noeud4->get_x()-noeud2->get_x(),noeud4->get_y()-noeud2->get_y(),noeud4->get_z()-noeud2->get_z());
96     OT_VECTEUR_3D n1n5(noeud5->get_x()-noeud1->get_x(),noeud5->get_y()-noeud1->get_y(),noeud5->get_z()-noeud1->get_z());
97     OT_VECTEUR_3D n1n6(noeud6->get_x()-noeud1->get_x(),noeud6->get_y()-noeud1->get_y(),noeud6->get_z()-noeud1->get_z());
98     OT_VECTEUR_3D n1n7(noeud7->get_x()-noeud1->get_x(),noeud7->get_y()-noeud1->get_y(),noeud7->get_z()-noeud1->get_z());
99     OT_VECTEUR_3D n2n5(noeud5->get_x()-noeud2->get_x(),noeud5->get_y()-noeud2->get_y(),noeud5->get_z()-noeud2->get_z());
100     OT_VECTEUR_3D n2n6(noeud6->get_x()-noeud2->get_x(),noeud6->get_y()-noeud2->get_y(),noeud6->get_z()-noeud2->get_z());
101     OT_VECTEUR_3D n2n7(noeud7->get_x()-noeud2->get_x(),noeud7->get_y()-noeud2->get_y(),noeud7->get_z()-noeud2->get_z());
102     OT_VECTEUR_3D normal1=n1n3&n1n2;
103     OT_VECTEUR_3D normal2=n1n2&n1n4;
104     OT_VECTEUR_3D normal3=n2n3&n2n4;
105     OT_VECTEUR_3D normal4=n1n4&n1n3;
106 francois 553 double eps1=0.33333333333*(n1n3.diff()+n1n2.diff()+n1n5.diff());
107     double eps2=0.33333333333*(n1n4.diff()+n1n2.diff()+n1n5.diff());
108     double eps3=0.33333333333*(n2n4.diff()+n2n3.diff()+n2n5.diff());
109     double eps4=0.33333333333*(n1n4.diff()+n1n3.diff()+n1n5.diff());
110     //eps=18.*eps*eps*eps*1e-6;
111     eps1=18.*eps1*eps1*eps1*EPS_BASE_RELATIVE;
112     eps2=18.*eps2*eps2*eps2*EPS_BASE_RELATIVE;
113     eps3=18.*eps3*eps3*eps3*EPS_BASE_RELATIVE;
114     eps4=18.*eps4*eps4*eps4*EPS_BASE_RELATIVE;
115 5 if (commun1==0)
116 francois 553 if (normal1*n1n5<eps1)
117     if (normal2*n1n5<eps2)
118     if (normal3*n2n5<eps3)
119     if (normal4*n1n5<eps4) return true;
120     //eps=18.*eps*eps*eps*1e-6;
121     eps1=0.33333333333*(n1n3.diff()+n1n2.diff()+n1n6.diff());
122     eps2=0.33333333333*(n1n4.diff()+n1n2.diff()+n1n6.diff());
123     eps3=0.33333333333*(n2n4.diff()+n2n3.diff()+n2n6.diff());
124     eps4=0.33333333333*(n1n4.diff()+n1n3.diff()+n1n6.diff());
125     eps1=18.*eps1*eps1*eps1*EPS_BASE_RELATIVE;
126     eps2=18.*eps2*eps2*eps2*EPS_BASE_RELATIVE;
127     eps3=18.*eps3*eps3*eps3*EPS_BASE_RELATIVE;
128     eps4=18.*eps4*eps4*eps4*EPS_BASE_RELATIVE;
129 5 if (commun2==0)
130 francois 553 if (normal1*n1n6<eps1)
131     if (normal2*n1n6<eps2)
132     if (normal3*n2n6<eps3)
133     if (normal4*n1n6<eps4) return true;
134     //eps=18.*eps*eps*eps*1e-6;
135     eps1=0.33333333333*(n1n3.diff()+n1n2.diff()+n1n7.diff());
136     eps2=0.33333333333*(n1n4.diff()+n1n2.diff()+n1n7.diff());
137     eps3=0.33333333333*(n2n4.diff()+n2n3.diff()+n2n7.diff());
138     eps4=0.33333333333*(n1n4.diff()+n1n3.diff()+n1n7.diff());
139     eps1=18.*eps1*eps1*eps1*EPS_BASE_RELATIVE;
140     eps2=18.*eps2*eps2*eps2*EPS_BASE_RELATIVE;
141     eps3=18.*eps3*eps3*eps3*EPS_BASE_RELATIVE;
142     eps4=18.*eps4*eps4*eps4*EPS_BASE_RELATIVE;
143 5 if (commun3==0)
144 francois 553 if (normal1*n1n7<eps1)
145     if (normal2*n1n7<eps2)
146     if (normal3*n2n7<eps3)
147     if (normal4*n1n7<eps4) return true;
148 5 if (inter_triangle_triangle(noeud1,noeud3,noeud2,noeud5,noeud6,noeud7)==true) return true;
149     if (inter_triangle_triangle(noeud1,noeud2,noeud4,noeud5,noeud6,noeud7)==true) return true;
150     if (inter_triangle_triangle(noeud2,noeud3,noeud4,noeud5,noeud6,noeud7)==true) return true;
151     if (inter_triangle_triangle(noeud1,noeud4,noeud3,noeud5,noeud6,noeud7)==true) return true;
152     }
153    
154     return inter;
155     }
156    
157    
158     int MAILLEUR3D::inter_triangle_triangle(MG_NOEUD* noeud1,MG_NOEUD* noeud2,MG_NOEUD* noeud3,MG_NOEUD* noeud4,MG_NOEUD* noeud5,MG_NOEUD* noeud6)
159     {
160     // filtre sur numero
161     int nb_commun=0;
162     if ( (noeud4==noeud1) || (noeud4==noeud2) || (noeud4==noeud3) ) nb_commun++;
163     if ( (noeud5==noeud1) || (noeud5==noeud2) || (noeud5==noeud3) ) nb_commun++;
164     if ( (noeud6==noeud1) || (noeud6==noeud2) || (noeud6==noeud3) ) nb_commun++;
165     if (nb_commun==3) return(false);
166     // filtre sur position
167     double bxmin=std::min(noeud1->get_x(),noeud2->get_x());
168     bxmin=std::min(bxmin,noeud3->get_x());
169     double bxmax=std::max(noeud1->get_x(),noeud2->get_x());
170     bxmax=std::max(bxmax,noeud3->get_x());
171     double bymin=std::min(noeud1->get_y(),noeud2->get_y());
172     bymin=std::min(bymin,noeud3->get_y());
173     double bymax=std::max(noeud1->get_y(),noeud2->get_y());
174     bymax=std::max(bymax,noeud3->get_y());
175     double bzmin=std::min(noeud1->get_z(),noeud2->get_z());
176     bzmin=std::min(bzmin,noeud3->get_z());
177     double bzmax=std::max(noeud1->get_z(),noeud2->get_z());
178     bzmax=std::max(bzmax,noeud3->get_z());
179     BOITE_3D boite_triangle1(bxmin,bymin,bzmin,bxmax,bymax,bzmax);
180     double b2xmin=std::min(noeud5->get_x(),noeud6->get_x());
181     b2xmin=std::min(b2xmin,noeud4->get_x());
182     double b2xmax=std::max(noeud5->get_x(),noeud6->get_x());
183     b2xmax=std::max(b2xmax,noeud4->get_x());
184     double b2ymin=std::min(noeud5->get_y(),noeud6->get_y());
185     b2ymin=std::min(b2ymin,noeud4->get_y());
186     double b2ymax=std::max(noeud5->get_y(),noeud6->get_y());
187     b2ymax=std::max(b2ymax,noeud4->get_y());
188     double b2zmin=std::min(noeud5->get_z(),noeud6->get_z());
189     b2zmin=std::min(b2zmin,noeud4->get_z());
190     double b2zmax=std::max(noeud5->get_z(),noeud6->get_z());
191     b2zmax=std::max(b2zmax,noeud4->get_z());
192     BOITE_3D boite_triangle2(b2xmin,b2ymin,b2zmin,b2xmax,b2ymax,b2zmax);
193     int inter=boite_triangle1*boite_triangle2;
194     if (inter==true)
195     {
196     inter=false;
197     OT_VECTEUR_3D n1n2(noeud2->get_x()-noeud1->get_x(),noeud2->get_y()-noeud1->get_y(),noeud2->get_z()-noeud1->get_z());
198     OT_VECTEUR_3D n1n3(noeud3->get_x()-noeud1->get_x(),noeud3->get_y()-noeud1->get_y(),noeud3->get_z()-noeud1->get_z());
199     OT_VECTEUR_3D n1n4(noeud4->get_x()-noeud1->get_x(),noeud4->get_y()-noeud1->get_y(),noeud4->get_z()-noeud1->get_z());
200     OT_VECTEUR_3D n1n5(noeud5->get_x()-noeud1->get_x(),noeud5->get_y()-noeud1->get_y(),noeud5->get_z()-noeud1->get_z());
201     OT_VECTEUR_3D n1n6(noeud6->get_x()-noeud1->get_x(),noeud6->get_y()-noeud1->get_y(),noeud6->get_z()-noeud1->get_z());
202     OT_VECTEUR_3D normal=n1n2&n1n3;
203     double ps1=normal*n1n4;
204     double ps2=normal*n1n5;
205     double ps3=normal*n1n6;
206     double eps1=0.333333333*(n1n2.diff()+n1n3.diff()+n1n4.diff());
207     double eps2=0.333333333*(n1n2.diff()+n1n3.diff()+n1n5.diff());
208     double eps3=0.333333333*(n1n2.diff()+n1n3.diff()+n1n6.diff());
209 francois 553 // eps1=18.*eps1*eps1*eps1*1e-6;
210     // eps2=18.*eps2*eps2*eps2*1e-6;
211     // eps3=18.*eps3*eps3*eps3*1e-6;
212     eps1=18.*eps1*eps1*eps1*EPS_BASE_RELATIVE;
213     eps2=18.*eps2*eps2*eps2*EPS_BASE_RELATIVE;
214     eps3=18.*eps3*eps3*eps3*EPS_BASE_RELATIVE;
215 5 if ((ps1<-eps1) && (ps2<-eps2) && (ps3<-eps3) ) return false;
216     if ((ps1>eps1) && (ps2>eps2) && (ps3>eps3) ) return false;
217     if (inter_segment_triangle(noeud1,noeud2,noeud3,noeud4,noeud5)==true) return true;
218     if (inter_segment_triangle(noeud1,noeud2,noeud3,noeud5,noeud6)==true) return true;
219     if (inter_segment_triangle(noeud1,noeud2,noeud3,noeud4,noeud6)==true) return true;
220     if (inter_segment_triangle(noeud4,noeud5,noeud6,noeud1,noeud2)==true) return true;
221     if (inter_segment_triangle(noeud4,noeud5,noeud6,noeud1,noeud3)==true) return true;
222     if (inter_segment_triangle(noeud4,noeud5,noeud6,noeud2,noeud3)==true) return true;
223     }
224     return inter;
225    
226    
227     }
228    
229     int MAILLEUR3D::inter_segment_triangle(MG_NOEUD* noeud1,MG_NOEUD* noeud2,MG_NOEUD* noeud3,MG_NOEUD* noeud4,MG_NOEUD* noeud5)
230     {
231     // filtre sur numero
232     int nb_commun=0;
233     if ( (noeud4==noeud1) || (noeud4==noeud2) || (noeud4==noeud3) ) nb_commun++;
234     if ( (noeud5==noeud1) || (noeud5==noeud2) || (noeud5==noeud3) ) nb_commun++;
235     if (nb_commun==2) return(false);
236    
237     int inter=false;
238     OT_VECTEUR_3D g1(noeud2->get_x()-noeud1->get_x(),noeud2->get_y()-noeud1->get_y(),noeud2->get_z()-noeud1->get_z());
239     OT_VECTEUR_3D g2(noeud3->get_x()-noeud1->get_x(),noeud3->get_y()-noeud1->get_y(),noeud3->get_z()-noeud1->get_z());
240     OT_VECTEUR_3D g3(noeud5->get_x()-noeud4->get_x(),noeud5->get_y()-noeud4->get_y(),noeud5->get_z()-noeud4->get_z());
241     OT_MATRICE_3D systeme(g1,g2,g3);
242     double det=systeme.get_determinant();
243     double eps=0.333333333333*(g1.diff()+g2.diff()+g3.diff());
244 francois 553 //eps=eps*eps*0.0018;
245     eps=18.*eps*eps*eps*EPS_BASE_RELATIVE;
246 francois 35 if (OPERATEUR::egal(det,0.0,eps)==true)
247 5 {
248 francois 35 OT_VECTEUR_3D g3b(noeud5->get_x()-noeud1->get_x(),noeud5->get_y()-noeud1->get_y(),noeud5->get_z()-noeud1->get_z());
249     OT_MATRICE_3D systeme2(g1,g2,g3b);
250     double det2=systeme2.get_determinant();
251     double eps2=0.333333333333*(g1.diff()+g2.diff()+g3b.diff());
252 francois 553 //eps2=18.*eps2*eps2*eps2*1e-6;
253     eps2=18.*eps2*eps2*eps2*EPS_BASE_RELATIVE;
254 francois 35 if (OPERATEUR::egal(det2,0.0,eps)==true) // cas 2D
255     {
256     if (inter_segment_segment(noeud1,noeud2,noeud4,noeud5)==true) return true;
257     if (inter_segment_segment(noeud1,noeud3,noeud4,noeud5)==true) return true;
258     if (inter_segment_segment(noeud2,noeud3,noeud4,noeud5)==true) return true;
259     }
260     else return false;
261 5 }
262     else
263     {
264     if ( (noeud4==noeud1) || (noeud4==noeud2) || (noeud4==noeud3) ) return false;
265     if ( (noeud5==noeud1) || (noeud5==noeud2) || (noeud5==noeud3) ) return false;
266     double x,y,z;
267     x=1.0/det*(g2.get_y()*g3.get_z()-g2.get_z()*g3.get_y());
268     y=1.0/det*(g3.get_x()*g2.get_z()-g2.get_x()*g3.get_z());
269     z=1.0/det*(g2.get_x()*g3.get_y()-g2.get_y()*g3.get_x());
270     OT_VECTEUR_3D g1b(x,y,z);
271     x=1.0/det*(g3.get_y()*g1.get_z()-g1.get_y()*g3.get_z());
272     y=1.0/det*(g1.get_x()*g3.get_z()-g3.get_x()*g1.get_z());
273     z=1.0/det*(g3.get_x()*g1.get_y()-g1.get_x()*g3.get_y());
274     OT_VECTEUR_3D g2b(x,y,z);
275     x=1.0/det*(g1.get_y()*g2.get_z()-g1.get_z()*g2.get_y());
276     y=1.0/det*(g2.get_x()*g1.get_z()-g1.get_x()*g2.get_z());
277     z=1.0/det*(g1.get_x()*g2.get_y()-g1.get_y()*g2.get_x());
278     OT_VECTEUR_3D g3b(x,y,z);
279     OT_VECTEUR_3D n1n4(noeud4->get_x()-noeud1->get_x(),noeud4->get_y()-noeud1->get_y(),noeud4->get_z()-noeud1->get_z());
280     double alpha1=n1n4*g1b;
281     double alpha2=n1n4*g2b;
282     double alpha3=-(n1n4*g3b);
283     double alpha4=1-alpha1-alpha2;
284 francois 553 //double eps=0.000001;
285     double eps=0.25*(g1.diff()+g2.diff()+g3.diff()+n1n4.diff());
286     eps=EPS_BASE_RELATIVE*(eps+1./eps);
287 5 if ((alpha1>-eps) && (alpha1<1.+eps))
288     if ((alpha2>-eps) && (alpha2<1.+eps))
289     if ((alpha3>-eps) && (alpha3<1.+eps))
290     if ((alpha4>-eps) && (alpha4<1.+eps)) return true;
291     return false;
292     }
293     return inter;
294     }
295    
296    
297    
298 francois 553 /*int MAILLEUR3D::noeud_est_triangle(MG_NOEUD* noeud1,MG_NOEUD* noeud2,MG_NOEUD* noeud3,double x,double y,double z)
299 5 {
300     OT_VECTEUR_3D n1n2(noeud2->get_x()-noeud1->get_x(),noeud2->get_y()-noeud1->get_y(),noeud2->get_z()-noeud1->get_z());
301     OT_VECTEUR_3D n1n3(noeud3->get_x()-noeud1->get_x(),noeud3->get_y()-noeud1->get_y(),noeud3->get_z()-noeud1->get_z());
302     OT_VECTEUR_3D n=n1n2&n1n3;
303     OT_VECTEUR_3D n1p(x-noeud1->get_x(),y-noeud1->get_y(),z-noeud1->get_z());
304     OT_VECTEUR_3D n2p(x-noeud2->get_x(),y-noeud2->get_y(),z-noeud2->get_z());
305     OT_VECTEUR_3D n3p(x-noeud3->get_x(),y-noeud3->get_y(),z-noeud3->get_z());
306     OT_VECTEUR_3D pv1=n1n2&n1p;
307     double ps1=pv1*n;
308 francois 553 //double eps=n1n2*n1n2;
309     //eps=eps*eps*4*0.0001;
310     double eps=0.25*(2.*n1n2.diff()+n1n3.diff()+n1p.diff());
311     eps=16.*eps*eps*eps*eps*EPS_BASE_RELATIVE;
312 5 if (ps1>(-eps))
313     {
314     OT_VECTEUR_3D n2n3(noeud3->get_x()-noeud2->get_x(),noeud3->get_y()-noeud2->get_y(),noeud3->get_z()-noeud2->get_z());
315     OT_VECTEUR_3D pv2=n2n3&n2p;
316 francois 553 double ps2=pv2*n;
317     double eps=0.25*(n2n3.diff()+n2p.diff()+n1n2.diff()+n1n3.diff());
318     eps=16.*eps*eps*eps*eps*EPS_BASE_RELATIVE;
319 5 if (ps2>(-eps))
320     {
321     OT_VECTEUR_3D pv3=n1n3&n3p;
322     double ps3=(-1)*(pv3*n);
323 francois 553 double eps=0.25*(n1n3.diff()+n3p.diff()+n1n2.diff()+n1n3.diff());
324     eps=16.*eps*eps*eps*eps*EPS_BASE_RELATIVE;
325     if (ps3>(-eps)) return true;
326 5 }
327 francois 553 arete 66
328 5
329     }
330     return false;
331     }
332    
333 francois 553 */
334 5
335    
336 francois 553 int MAILLEUR3D::inter_segment_segment(MG_NOEUD *A, MG_NOEUD *B,MG_NOEUD *C,MG_NOEUD *D)
337     {
338     /*static int passe=0;
339     passe++;
340     if (passe==366677)
341     std::cout << " pas pareil a venir" << std::endl;
342     */
343     int res1=inter_segment_segment1(A,B,C,D);
344     //int res2=inter_segment_segment2(A,B,C,D);
345     //if (res1!=res2)
346     //std::cout << " pas pareil" << std::endl;
347     return res1;
348     }
349 5
350    
351 francois 553 int MAILLEUR3D::inter_segment_segment1(MG_NOEUD *A, MG_NOEUD *B,MG_NOEUD *C,MG_NOEUD *D)
352     {
353     double2 a0(A->get_coord()[0]);
354     double2 a1(A->get_coord()[1]);
355     double2 a2(A->get_coord()[2]);
356     OT_VECTEUR_3DD OA(a0,a1,a2);
357     double2 b0(B->get_coord()[0]);
358     double2 b1(B->get_coord()[1]);
359     double2 b2(B->get_coord()[2]);
360     OT_VECTEUR_3DD OB(b0,b1,b2);
361     double2 c0(C->get_coord()[0]);
362     double2 c1(C->get_coord()[1]);
363     double2 c2(C->get_coord()[2]);
364     OT_VECTEUR_3DD OC(c0,c1,c2);
365     double2 d0(D->get_coord()[0]);
366     double2 d1(D->get_coord()[1]);
367     double2 d2(D->get_coord()[2]);
368     OT_VECTEUR_3DD OD(d0,d1,d2);
369     OT_VECTEUR_3DD u(OA,OB);
370     OT_VECTEUR_3DD v(OC,OD);
371     OT_VECTEUR_3DD AC(OA,OC);
372    
373     OT_VECTEUR_3DD uv=u&v;
374     double2 ZERO(0.);
375     double2 UN(1.);
376     double2 longueur=uv.get_longueur();
377     if (uv.get_longueur()==ZERO)
378     {
379     //std::cout << " segment colineaire" << std::endl;
380     OT_VECTEUR_3DD test=AC&u;
381     if (test.get_longueur()==ZERO)
382     {
383     //std::cout << " support segment confondu" << std::endl;
384     OT_VECTEUR_3DD AC(OA,OC);
385     OT_VECTEUR_3DD AD(OA,OD);
386     OT_VECTEUR_3DD BC(OB,OC);
387     OT_VECTEUR_3DD BD(OB,OD);
388     if ((AC*AD>=ZERO) && (AC*u<=ZERO) ) return 0;//std::cout << " pas intersection" << std::endl;
389     else if ((BC*BD>=ZERO) && (BC*u>=ZERO)) return 0;//std::cout << " pas intersection" << std::endl;
390     else return 1;//std::cout << "intersection" << std::endl;
391    
392     }
393     else return 0;//std::cout << " segment parrallele" << std::endl;
394 5
395 francois 553 }
396     else
397     {
398     OT_VECTEUR_3DD P2=v&uv;
399     double2 a2=P2.get_x();
400     double2 b2=P2.get_y();
401     double2 c2=P2.get_z();
402     double2 d2=ZERO-a2*C->get_coord()[0]-b2*C->get_coord()[1]-c2*C->get_coord()[2];
403     double2 t1=ZERO-a2*A->get_coord()[0]-b2*A->get_coord()[1]-c2*A->get_coord()[2]-d2;
404     t1=t1/(a2*u.get_x()+b2*u.get_y()+c2*u.get_z());
405     OT_VECTEUR_3DD H=OA+t1*u;
406    
407     OT_VECTEUR_3DD P1=u&uv;
408     double2 a1=P1.get_x();
409     double2 b1=P1.get_y();
410     double2 c1=P1.get_z();
411     double2 d1=ZERO-a1*A->get_coord()[0]-b1*A->get_coord()[1]-c1*A->get_coord()[2];
412     double2 t2=ZERO-a1*C->get_coord()[0]-b1*C->get_coord()[1]-c1*C->get_coord()[2]-d1;
413     t2=t2/(a1*v.get_x()+b1*v.get_y()+c1*v.get_z());
414     OT_VECTEUR_3DD H2=OC+t2*v;
415    
416     OT_VECTEUR_3DD PPC2(H,H2);
417    
418     if (PPC2.get_longueur()==ZERO)
419     {
420     //std::cout << " segment coplanaire" << std::endl;
421     //std::cout << " point d'intersection" << H << std::endl;
422     //std::cout << " parametre t1 " << t1 << std::endl;
423     //std::cout << " parametre t2 " << t2 << std::endl;
424     if (t1>ZERO)
425     if (t1<UN)
426     if (t2>ZERO)
427     if (t2<UN) return 1;//std::cout << " intersection " << std::endl;
428     else return 0;//std::cout << " pas intersection " << std::endl;
429    
430    
431     }
432     else
433     return 0;//std::cout << " segment 3D" << std::endl;
434     }
435     return 0;
436     }
437    
438    
439 5 // PARTIE DE CODE RECUPERER DE MON PROJET DE DEA
440     // TRES ILLISIBLE MAIS FONCTIONNEL
441     // POIL AU MARRON
442    
443     #define PSCA(a,b) (a[0]*b[0]+a[1]*b[1]+a[2]*b[2])
444 francois 507 #define EGAL(x,y,eps) (double)fabs((double)(x-y))<eps
445 5 #define DETER(a,b,c,d) (a*d-b*c)
446    
447 francois 553
448     int MAILLEUR3D::inter_segment_segment2(MG_NOEUD* noeud1,MG_NOEUD* noeud2,MG_NOEUD* noeud3,MG_NOEUD* noeud4)
449 5 {
450    
451     double ab[3];
452     double nm[3];
453     double am[3];
454     ab[0]=noeud2->get_x()-noeud1->get_x();
455     ab[1]=noeud2->get_y()-noeud1->get_y();
456     ab[2]=noeud2->get_z()-noeud1->get_z();
457     nm[0]=noeud3->get_x()-noeud4->get_x();
458     nm[1]=noeud3->get_y()-noeud4->get_y();
459     nm[2]=noeud3->get_z()-noeud4->get_z();
460     am[0]=noeud3->get_x()-noeud1->get_x();
461     am[1]=noeud3->get_y()-noeud1->get_y();
462     am[2]=noeud3->get_z()-noeud1->get_z();
463     int equation[4];
464     equation[0]=1; /* etat de l'equation 0 */
465     equation[1]=1;
466     equation[2]=1;
467     equation[3]=3; /* cette variable comporte le bilan du nombre d'equation */
468 francois 553 double epsab=0.3333333333333*(fabs(ab[0])+fabs(ab[1])+fabs(ab[2]))*EPS_BASE_RELATIVE;
469     double epsnm=0.3333333333333*(fabs(nm[0])+fabs(nm[1])+fabs(nm[2]))*EPS_BASE_RELATIVE;
470     double epsam=0.3333333333333*(fabs(am[0])+fabs(am[1])+fabs(am[2]))*EPS_BASE_RELATIVE;
471     //eps=eps*0.0001;
472     //eps2=eps2*0.0001;
473 5 /* recherche du nombre d'equation -> inter franche ou para ou confondu */
474 francois 553 if ( (EGAL(ab[0],0,epsab)) && (EGAL(nm[0],0,epsnm)) )
475     if (EGAL(am[0],0,epsam)) equation[0]=0; else return false;
476     if ( (EGAL(ab[1],0,epsab)) && (EGAL(nm[1],0,epsnm)) )
477     if (EGAL(am[1],0,epsam)) equation[1]=0; else return false;
478     if ( (EGAL(ab[2],0,epsab)) && (EGAL(nm[2],0,epsnm)) )
479     if (EGAL(am[2],0,epsam)) equation[2]=0; else return false;
480 5 equation[3]=equation[0]+equation[1]+equation[2];
481     if (equation[3]==3)
482     {
483     double det=DETER(ab[0],nm[0],ab[1],nm[1]);
484 francois 553 double eps2=0.25*(fabs(ab[0])+fabs(nm[0])+fabs(ab[1])+fabs(nm[1]));
485     eps2=4*eps2*eps2*EPS_BASE_RELATIVE;
486 5 if (fabs(det)>eps2)
487     {
488     det=1/det;
489     double sol1=det*DETER(am[0],nm[0],am[1],nm[1]);
490     double sol2=det*DETER(ab[0],am[0],ab[1],am[1]);
491 francois 553 double epssol1=0.25*(fabs(am[0])+fabs(nm[0])+fabs(am[1])+fabs(nm[1]));
492     epssol1=4*epssol1*epssol1*det*EPS_BASE_RELATIVE+2*epssol1*epssol1*det*det*eps2;
493     double epssol2=0.25*(fabs(ab[0])+fabs(am[0])+fabs(ab[1])+fabs(am[1]));
494     epssol2=4*epssol2*epssol2*det*EPS_BASE_RELATIVE+2*epssol2*epssol2*det*det*eps2;
495     double eps3=0.5*(ab[2]+nm[2]);
496     eps3=(sol1+sol2)*eps3*EPS_BASE_RELATIVE+eps3*epssol1+eps3*epssol2;
497     if (fabs((double)(sol1*ab[2]-sol2*nm[2]-am[2]))>eps3) return false;
498     return(examine_solution(sol1,sol2,epssol1,epssol2,1));
499 5 }
500     else
501     {
502     equation[0]=0;
503     equation[3]=2;
504     /* on verifie la compatibilite des deux equations dont le det est nul*/
505 francois 553 double tmp,eps;
506     if (!(EGAL(ab[0],0.,(0.333333333*(fabs(ab[0])+fabs(ab[1])+fabs(ab[2]))*EPS_BASE_RELATIVE))))
507     {
508     tmp=ab[1]*am[0]/ab[0];
509     eps=0.5*(fabs(ab[1])+fabs(am[0]));
510     eps=2*eps*eps/ab[0]+eps*eps/ab[0];
511     eps=eps*EPS_BASE_RELATIVE;
512     }
513     else
514     {
515     tmp=nm[1]*am[0]/nm[0];
516     eps=0.5*(fabs(nm[1])+fabs(am[0]));
517     eps=2*eps*eps/nm[0]+eps*eps/nm[0];
518     eps=eps*EPS_BASE_RELATIVE;
519    
520     }
521 5 if (!(EGAL(tmp,am[1],eps))) return false;
522     }
523     }
524     if (equation[3]==2)
525     {
526     /* on repere les equations qui existent */
527     int ne1;
528     int ne2;
529     if (equation[0]!=0)
530     {
531     ne1=0;
532     if (equation[1]!=0) ne2=1; else ne2=2;
533     }
534     else
535     {
536     ne1=1;
537     ne2=2;
538     }
539    
540 francois 553 double det2=DETER(ab[ne1],nm[ne1],ab[ne2],nm[ne2]);
541     double eps2=0.25*(fabs(ab[ne1])+fabs(nm[ne1])+fabs(ab[ne2])+fabs(nm[ne2]));
542     eps2=4*eps2*eps2*EPS_BASE_RELATIVE;
543     if (fabs(det2)>eps2)
544 5 {
545 francois 553 double det=1/det2;
546 5 double sol1=det*DETER(am[ne1],nm[ne1],am[ne2],nm[ne2]);
547     double sol2=det*DETER(ab[ne1],am[ne1],ab[ne2],am[ne2]);
548 francois 553 double epssol1=0.25*(fabs(am[ne1])+fabs(nm[ne1])+fabs(am[ne2])+fabs(nm[ne2]));
549     epssol1=4*epssol1*epssol1*det*EPS_BASE_RELATIVE+2*epssol1*epssol1*det*det*eps2;
550     double epssol2=0.25*(fabs(ab[ne1])+fabs(am[ne1])+fabs(ab[ne2])+fabs(am[ne2]));
551     epssol2=4*epssol2*epssol2*det*EPS_BASE_RELATIVE+2*epssol2*epssol2*det*det*eps2;
552     return(examine_solution(sol1,sol2,epssol1,epssol2,1));
553 5 }
554     else
555     {
556     equation[ne1]=0;
557     equation[3]=1;
558     /* on verifie la compatibilite des deux equations dont le det est nul */
559 francois 553 double tmp,eps;
560     if (!(EGAL(ab[ne1],0.,(0.3333333333333*(fabs(ab[0])+fabs(ab[1])+fabs(ab[2]))*EPS_BASE_RELATIVE))))
561     {
562     tmp=ab[ne2]*am[ne1]/ab[ne1];
563     eps=0.5*(fabs(ab[ne2])+fabs(am[ne1]));
564     eps=2*eps*eps/ab[ne1]+eps*eps/ab[ne1];
565     eps=eps*EPS_BASE_RELATIVE;
566     }
567     else
568     {
569     tmp=nm[ne2]*am[ne1]/nm[ne1];
570     eps=0.5*(fabs(nm[ne2])+fabs(am[ne1]));
571     eps=2*eps*eps/nm[ne1]+eps*eps/nm[ne1];
572     eps=eps*EPS_BASE_RELATIVE;
573    
574     }
575 5 if (!(EGAL(tmp,am[ne2],eps))) return false;
576     }
577    
578     }
579     if (equation[3]==1)
580     {
581     /* on repere l' equation qui existe */
582     int ne1;
583     if (equation[0]!=0) ne1=0; else
584     if (equation[1]!=0) ne1=1; else ne1=2;
585     double an[3];
586     an[0]=noeud4->get_x()-noeud1->get_x();
587     an[1]=noeud4->get_y()-noeud1->get_y();
588     an[2]=noeud4->get_z()-noeud1->get_z();
589 francois 553 double tmp=1./ab[ne1];
590 5 double sol1=am[ne1]*tmp;
591     double sol2=an[ne1]*tmp;
592 francois 553 double epssol1=2*am[ne1]*tmp*EPS_BASE_RELATIVE;
593     double epssol2=2*an[ne1]*tmp*EPS_BASE_RELATIVE;
594     return(examine_solution(sol1,sol2,epssol1,epssol2,2));
595 5 }
596     return 0;
597     }
598    
599    
600    
601 francois 553
602    
603    
604    
605     int MAILLEUR3D::examine_solution(double sol1,double sol2,double epssol1,double epssol2,int type)
606 5 {
607 francois 553 epssol1=0.0001;
608     epssol2=0.0001;
609 5
610 francois 553 if ((epssol1>1e-0) || (epssol2>1e-0))
611     std::cout << epssol1 << " " << epssol2 << std::endl;
612 5 if (type==1)
613     {
614 francois 553 if ( (sol1>epssol1) && ((sol1)<(1-epssol1)) && (sol2>epssol2) && ((sol2)<(1-epssol2)) ) return true;
615     if ( ( (EGAL(sol1,0,epssol1)) || (EGAL(sol1,1,epssol1))) && ( (sol2>epssol2) && ((sol2)<(1-epssol2)) ) ) return true;
616     if ( ( (EGAL(sol2,0,epssol2)) || (EGAL(sol2,1,epssol2))) && ( (sol1>epssol1) && ((sol1)<(1-epssol1)) ) ) return true;
617     // if ( (sol1>epssol1) && ((sol1)<(1-epssol1)) && (sol2>-epssol2) && ((sol2)<(1.-epssol2)) ) return true;
618     // if ( (sol2>epssol2) && ((sol2)<(1-epssol2)) && (sol1>-epssol1) && ((sol1)<(1.-epssol1)) ) return true;
619 5
620     }
621     if (type==2)
622     {
623 francois 553 if ( (sol1>epssol1) && ((sol1)<(1-epssol1)) ) return true;
624     if ( (sol2>epssol2) && ((sol2)<(1-epssol2)) ) return true;
625     if ( ((sol1)>(1+epssol1)) && ((-sol2)>epssol2) ) return true;
626     if ( ((sol2)>(1+epssol2)) && ((-sol1)>epssol1) ) return true;
627 5 }
628     return false;
629 francois 553 /*double epsilon=0.0001;
630    
631     if (type==1)
632     {
633     if ( (sol1>epsilon) && ((sol1)<(1-epsilon)) && (sol2>epsilon) && ((sol2)<(1-epsilon)) ) return 1;
634     if ( ( (EGAL(sol1,0,epsilon)) || (EGAL(sol1,1,epsilon))) && ( (sol2>epsilon) && ((sol2)<(1-epsilon)) ) ) return 1;
635     if ( ( (EGAL(sol2,0,epsilon)) || (EGAL(sol2,1,epsilon))) && ( (sol1>epsilon) && ((sol1)<(1-epsilon)) ) ) return 1;
636     if ( (sol1>epsilon) && ((sol1)<(1-epsilon)) && (sol2>(-0.1-epsilon)) && ((sol2)<(1.1-epsilon)) ) return 1;
637     if ( (sol2>epsilon) && ((sol2)<(1-epsilon)) && (sol1>(-0.1-epsilon)) && ((sol1)<(1.1-epsilon)) ) return 1;
638    
639     }
640     if (type==2)
641     {
642     if ( (sol1>epsilon) && ((sol1)<(1-epsilon)) ) return 1;
643     if ( (sol2>epsilon) && ((sol2)<(1-epsilon)) ) return 1;
644     if ( ((sol1)>(1+epsilon)) && ((-sol2)>epsilon) ) return 1;
645     if ( ((sol2)>(1+epsilon)) && ((-sol1)>epsilon) ) return 1;
646     }
647     return 0;*/
648 5 }
649     #undef EGAL
650     #undef PSCA
651     #undef DETER
652     ///FIN DU CODE DE DEA
653     //---------------------------------------------------------------------------
654    
655     //#pragma package(smart_init)