ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/mailleur3d_intersection.cpp
Revision: 1158
Committed: Thu Jun 13 22:18:49 2024 UTC (10 months, 4 weeks ago) by francois
Original Path: magic/lib/mailleur_auto/src/mailleur3d_intersection.cpp
File size: 28476 byte(s)
Log Message:
compatibilité Ubuntu 22.04
Suppression des refeences à Windows
Ajout d'une banière

File Contents

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