MAGiC  V5.0
Mailleurs Automatiques de Géometries intégrés à la Cao
mailleur3d_intersection.cpp
Aller à la documentation de ce fichier.
1 //####//------------------------------------------------------------
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 
23 
24 #include "gestionversion.h"
25 
26 
27 
28 #include <math.h>
29 #include "mailleur3d.h"
30 #include "ot_mathematique.h"
31 #include "ot_doubleprecision.h"
32 
33 
34 
35 
36 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  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  if (commun1==0)
112  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  if (commun2==0)
126  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  if (commun3==0)
140  if (normal1*n1n7<eps1)
141  if (normal2*n1n7<eps2)
142  if (normal3*n2n7<eps3)
143  if (normal4*n1n7<eps4) return true;
144  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  // 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  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 eps=18.*eps*eps*eps*EPS_BASE_RELATIVE;
238 if (OPERATEUR::egal(det,0.0,eps)==true)
239  {
240  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  //eps2=18.*eps2*eps2*eps2*1e-6;
245  eps2=18.*eps2*eps2*eps2*EPS_BASE_RELATIVE;
246  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  }
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  //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  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 /*int MAILLEUR3D::noeud_est_triangle(MG_NOEUD* noeud1,MG_NOEUD* noeud2,MG_NOEUD* noeud3,double x,double y,double z)
291 {
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 double eps=0.25*(2.*n1n2.diff()+n1n3.diff()+n1p.diff());
301 eps=16.*eps*eps*eps*eps*EPS_BASE_RELATIVE;
302 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  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  if (ps2>(-eps))
310  {
311  OT_VECTEUR_3D pv3=n1n3&n3p;
312  double ps3=(-1)*(pv3*n);
313  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  }
317  arete 66
318 
319  }
320 return false;
321 }
322 
323 */
324 
325 
326 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 
338 
339 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  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  else return 1;//std::cout << "intersection" << std::endl;
379 
380  }
381  else return 0;//std::cout << " segment parrallele" << std::endl;
382 
383  }
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 
428 #define PSCA(a,b) (a[0]*b[0]+a[1]*b[1]+a[2]*b[2])
429 #define EGAL(x,y,eps) (double)fabs((double)(x-y))<eps
430 #define DETER(a,b,c,d) (a*d-b*c)
431 
432 
433 int MAILLEUR3D::inter_segment_segment2(MG_NOEUD* noeud1,MG_NOEUD* noeud2,MG_NOEUD* noeud3,MG_NOEUD* noeud4)
434 {
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 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 /* recherche du nombre d'equation -> inter franche ou para ou confondu */
457 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 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  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  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  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  }
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  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  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  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  {
528  double det=1/det2;
529  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  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  }
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  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  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  double tmp=1./ab[ne1];
573  double sol1=am[ne1]*tmp;
574  double sol2=an[ne1]*tmp;
575  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  }
579 return 0;
580 }
581 
582 
583 
584 
585 
586 
587 
588 int MAILLEUR3D::examine_solution(double sol1,double sol2,double epssol1,double epssol2,int type)
589 {
590 epssol1=0.0001;
591 epssol2=0.0001;
592 
593 if ((epssol1>1e-0) || (epssol2>1e-0))
594  std::cout << epssol1 << " " << epssol2 << std::endl;
595 if (type==1)
596  {
597  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 
603  }
604 if (type==2)
605  {
606  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  }
611 return false;
612 /*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 }
632 #undef EGAL
633 #undef PSCA
634 #undef DETER
635 
OT_MATRICE_3D::get_determinant
double get_determinant()
Definition: ot_mathematique.cpp:735
gestionversion.h
MG_NOEUD::v
double v
Definition: mg_noeud.h:107
MAILLEUR3D::inter_tetra_triangle
virtual int inter_tetra_triangle(MG_NOEUD *noeud1, MG_NOEUD *noeud2, MG_NOEUD *noeud3, MG_NOEUD *noeud4, MG_NOEUD *noeud5, MG_NOEUD *noeud6, MG_NOEUD *noeud7)
Definition: mailleur3d_intersection.cpp:36
MG_NOEUD::get_z
virtual double get_z(void)
Definition: mg_noeud.cpp:87
MAILLEUR3D::inter_segment_segment
virtual int inter_segment_segment(MG_NOEUD *noeud1, MG_NOEUD *noeud2, MG_NOEUD *noeud3, MG_NOEUD *noeud4)
Definition: mailleur3d_intersection.cpp:326
DETER
#define DETER(a, b, c, d)
Definition: mailleur3d_intersection.cpp:430
OT_VECTEUR_3DD::get_longueur
virtual double2 get_longueur(void)
Definition: ot_mathematique.cpp:1173
OT_VECTEUR_3D::get_x
virtual double get_x(void) const
Definition: ot_mathematique.cpp:417
MAILLEUR3D::inter_triangle_triangle
virtual int inter_triangle_triangle(MG_NOEUD *noeud1, MG_NOEUD *noeud2, MG_NOEUD *noeud3, MG_NOEUD *noeud4, MG_NOEUD *noeud5, MG_NOEUD *noeud6)
Definition: mailleur3d_intersection.cpp:154
double2
Definition: ot_doubleprecision.h:29
ot_doubleprecision.h
MG_NOEUD::u
double u
Definition: mg_noeud.h:106
OT_VECTEUR_3DD::get_z
virtual double2 get_z(void) const
Definition: ot_mathematique.cpp:1199
OT_VECTEUR_3DD::get_x
virtual double2 get_x(void) const
Definition: ot_mathematique.cpp:1183
OT_VECTEUR_3DD
Definition: ot_mathematique.h:350
EPS_BASE_RELATIVE
#define EPS_BASE_RELATIVE
Definition: mailleur3d.h:48
OPERATEUR::egal
static int egal(double a, double b, double eps)
Definition: ot_mathematique.cpp:1629
MG_NOEUD
Definition: mg_noeud.h:41
OT_MATRICE_3D
Definition: ot_mathematique.h:160
MG_NOEUD::get_coord
virtual double * get_coord(void)
Definition: mg_noeud.cpp:92
MAILLEUR3D::inter_segment_triangle
virtual int inter_segment_triangle(MG_NOEUD *noeud1, MG_NOEUD *noeud2, MG_NOEUD *noeud3, MG_NOEUD *noeud4, MG_NOEUD *noeud5)
Definition: mailleur3d_intersection.cpp:223
ot_mathematique.h
OT_VECTEUR_3D::get_y
virtual double get_y(void) const
Definition: ot_mathematique.cpp:423
MG_NOEUD::get_x
virtual double get_x(void)
Definition: mg_noeud.cpp:77
OT_VECTEUR_3D
Definition: ot_mathematique.h:94
MAILLEUR3D::inter_segment_segment2
virtual int inter_segment_segment2(MG_NOEUD *noeud1, MG_NOEUD *noeud2, MG_NOEUD *noeud3, MG_NOEUD *noeud4)
Definition: mailleur3d_intersection.cpp:433
MAILLEUR3D::examine_solution
virtual int examine_solution(double sol1, double sol2, double epssol1, double epssol2, int type)
Definition: mailleur3d_intersection.cpp:588
BOITE_3D
Definition: ot_boite_3d.h:27
mailleur3d.h
OT_VECTEUR_3DD::get_y
virtual double2 get_y(void) const
Definition: ot_mathematique.cpp:1191
OT_VECTEUR_3D::get_z
virtual double get_z(void) const
Definition: ot_mathematique.cpp:429
OT_VECTEUR_3D::diff
virtual double diff(void)
Definition: ot_mathematique.cpp:656
MAILLEUR3D::inter_segment_segment1
virtual int inter_segment_segment1(MG_NOEUD *noeud1, MG_NOEUD *noeud2, MG_NOEUD *noeud3, MG_NOEUD *noeud4)
Definition: mailleur3d_intersection.cpp:339
EGAL
#define EGAL(x, y, eps)
Definition: mailleur3d_intersection.cpp:429
MG_NOEUD::get_y
virtual double get_y(void)
Definition: mg_noeud.cpp:82