ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/mailleur3d_outil.cpp
Revision: 5
Committed: Tue Jun 12 20:26:34 2007 UTC (17 years, 11 months ago)
Original Path: magic/lib/mailleur/mailleur/src/mailleur3d_outil.cpp
File size: 31934 byte(s)
Log Message:

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_outil.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 <time.h>
32     #include "mailleur3d.h"
33     #include "ot_mathematique.h"
34     #include "m3d_noeud.h"
35     #include "m3d_tetra.h"
36     #include "m3d_triangle.h"
37     #include "const.h"
38    
39     //---------------------------------------------------------------------------
40    
41     //#pragma package(smart_init)
42    
43    
44    
45    
46    
47    
48    
49     double MAILLEUR3D::angle_front(MG_FRONT_3D* ft1,MG_FRONT_3D* ft2)
50     {
51     MG_NOEUD* noeud1=ft1->get_noeud1();
52     MG_NOEUD* noeud2=ft1->get_noeud2();
53     MG_NOEUD* noeud3=ft1->get_noeud3();
54     MG_NOEUD* noeud4=ft2->get_noeud1();
55     MG_NOEUD* noeud5=ft2->get_noeud2();
56     MG_NOEUD* noeud6=ft2->get_noeud3();
57    
58     MG_NOEUD* noeuda;
59     MG_NOEUD* noeudb;
60     MG_NOEUD* noeudc;
61     MG_NOEUD* noeudd;
62    
63    
64     if (!(noeud1==noeud4))
65     if (!(noeud1==noeud5))
66     if (!(noeud1==noeud6))
67     {
68     noeuda=noeud1;
69     noeudb=noeud2;
70     noeudc=noeud3;
71     }
72     if (!(noeud2==noeud4))
73     if (!(noeud2==noeud5))
74     if (!(noeud2==noeud6))
75     {
76     noeuda=noeud2;
77     noeudb=noeud1;
78     noeudc=noeud3;
79     }
80     if (!(noeud3==noeud4))
81     if (!(noeud3==noeud5))
82     if (!(noeud3==noeud6))
83     {
84     noeuda=noeud3;
85     noeudb=noeud1;
86     noeudc=noeud2;
87     }
88    
89     if (!(noeud4==noeudb))
90     if (!(noeud4==noeudc)) noeudd=noeud4;
91     if (!(noeud5==noeudb))
92     if (!(noeud5==noeudc)) noeudd=noeud5;
93     if (!(noeud6==noeudb))
94     if (!(noeud6==noeudc)) noeudd=noeud6;
95    
96     OT_VECTEUR_3D milieu(0.5*(noeudb->get_x()+noeudc->get_x()),0.5*(noeudb->get_y()+noeudc->get_y()),0.5*(noeudb->get_z()+noeudc->get_z()));
97     OT_VECTEUR_3D bc(noeudc->get_x()-noeudb->get_x(),noeudc->get_y()-noeudb->get_y(),noeudc->get_z()-noeudb->get_z());
98     bc.norme();
99     double xyz[3];
100     xyz[0]=noeuda->get_x();
101     xyz[1]=noeuda->get_y();
102     xyz[2]=noeuda->get_z();
103     OT_VECTEUR_3D n1(xyz,milieu.get_xyz());
104     n1.norme();
105     OT_VECTEUR_3D n=bc&n1;
106     n1=n&bc;
107     n1.norme();
108     xyz[0]=noeudd->get_x();
109     xyz[1]=noeudd->get_y();
110     xyz[2]=noeudd->get_z();
111     OT_VECTEUR_3D n2(xyz,milieu.get_xyz());
112     n2.norme();
113     n=bc&n2;
114     n2=n&bc;
115     n2.norme();
116    
117     double ps=n1*n2;
118     if (ps>1.) ps=1.;
119     if (ps<-1.) ps=-1.;
120     double angle=acos(ps);
121     OT_VECTEUR_3D n1n2(noeud2->get_x()-noeud1->get_x(),noeud2->get_y()-noeud1->get_y(),noeud2->get_z()-noeud1->get_z());
122     OT_VECTEUR_3D n1n3(noeud3->get_x()-noeud1->get_x(),noeud3->get_y()-noeud1->get_y(),noeud3->get_z()-noeud1->get_z());
123     OT_VECTEUR_3D nbnd(noeudd->get_x()-noeudb->get_x(),noeudd->get_y()-noeudb->get_y(),noeudd->get_z()-noeudb->get_z());
124     n=n1n2&n1n3;
125     n.norme();
126     nbnd.norme();
127    
128     ps=n*nbnd;
129     if (ps<-0.0001) angle=2*M_PI-angle;
130     return angle;
131    
132    
133     }
134    
135    
136    
137     double MAILLEUR3D::angle_front(MG_FRONT_3D* ft,MG_SEGMENT* seg)
138     {
139     MG_NOEUD* noeud1;
140     MG_NOEUD* noeud2;
141     MG_NOEUD* noeud3;
142     MG_NOEUD* noeud4;
143    
144    
145     if (seg->get_noeud1()==ft->get_noeud1())
146     {
147     noeud1=ft->get_noeud1();
148     noeud2=ft->get_noeud2();
149     noeud3=ft->get_noeud3();
150     noeud4=seg->get_noeud2();
151     }
152     if (seg->get_noeud1()==ft->get_noeud2())
153     {
154     noeud1=ft->get_noeud2();
155     noeud2=ft->get_noeud3();
156     noeud3=ft->get_noeud1();
157     noeud4=seg->get_noeud2();
158     }
159     if (seg->get_noeud1()==ft->get_noeud3())
160     {
161     noeud1=ft->get_noeud3();
162     noeud2=ft->get_noeud1();
163     noeud3=ft->get_noeud2();
164     noeud4=seg->get_noeud2();
165     }
166     if (seg->get_noeud2()==ft->get_noeud1())
167     {
168     noeud1=ft->get_noeud1();
169     noeud2=ft->get_noeud2();
170     noeud3=ft->get_noeud3();
171     noeud4=seg->get_noeud1();
172     }
173     if (seg->get_noeud2()==ft->get_noeud2())
174     {
175     noeud1=ft->get_noeud2();
176     noeud2=ft->get_noeud3();
177     noeud3=ft->get_noeud1();
178     noeud4=seg->get_noeud1();
179     }
180     if (seg->get_noeud2()==ft->get_noeud3())
181     {
182     noeud1=ft->get_noeud3();
183     noeud2=ft->get_noeud1();
184     noeud3=ft->get_noeud2();
185     noeud4=seg->get_noeud1();
186     }
187    
188     OT_VECTEUR_3D n1n2(noeud2->get_x()-noeud1->get_x(),noeud2->get_y()-noeud1->get_y(),noeud2->get_z()-noeud1->get_z());
189     OT_VECTEUR_3D n1n3(noeud3->get_x()-noeud1->get_x(),noeud3->get_y()-noeud1->get_y(),noeud3->get_z()-noeud1->get_z());
190     OT_VECTEUR_3D n1n4(noeud4->get_x()-noeud1->get_x(),noeud4->get_y()-noeud1->get_y(),noeud4->get_z()-noeud1->get_z());
191     OT_VECTEUR_3D n1n=n1n2&n1n3;
192     OT_MATRICE_3D mat(n1n2,n1n3,n1n);
193     OT_MATRICE_3D matx(n1n4,n1n3,n1n);
194     OT_MATRICE_3D maty(n1n2,n1n4,n1n);
195     OT_MATRICE_3D matz(n1n2,n1n3,n1n4);
196     double det=mat.get_determinant();
197     double xsi=matx.get_determinant()/det;
198     double eta=maty.get_determinant()/det;
199     double psi=matz.get_determinant()/det;
200    
201     double proj[3];
202     proj[0]=(1.-xsi-eta)*noeud1->get_x()+xsi*noeud2->get_x()+eta*noeud3->get_x();
203     proj[1]=(1.-xsi-eta)*noeud1->get_y()+xsi*noeud2->get_y()+eta*noeud3->get_y();
204     proj[2]=(1.-xsi-eta)*noeud1->get_z()+xsi*noeud2->get_z()+eta*noeud3->get_z();
205     OT_VECTEUR_3D n1nproj(noeud1->get_x()-proj[0],noeud1->get_y()-proj[1],noeud1->get_z()-proj[2]);
206     //n1n4.norme();
207     //n1nproj.norme();
208     double ps=n1n4*n1nproj;
209     if (!(OPERATEUR::egal(ps,0.0,0.0001))) ps=ps/n1n4.get_longueur()/n1nproj.get_longueur(); else ps=0.;
210     if (ps>1.) ps=1.;
211     if (ps<-1.) ps=-1.;
212     double angle=acos(ps);
213     if ((xsi<0.) || (eta<0.) ) angle=angle+M_PI/2.;
214     if (psi<0.) angle=2*M_PI-angle;
215     return angle;
216    
217     }
218    
219    
220    
221    
222    
223     class MG_TETRA* MAILLEUR3D::insere_tetra(class MG_VOLUME* mgvol,class MG_FRONT_3D* ft,MG_NOEUD* noeud4,int type,TPL_MAP_ENTITE<class MG_TRIANGLE*>& liste_intersection)
224     {
225     clock_t cpu1=clock(); //debug
226     liste_intersection.vide();
227     MG_NOEUD* noeud1=ft->get_noeud1();
228     MG_NOEUD* noeud2=ft->get_noeud2();
229     MG_NOEUD* noeud3=ft->get_noeud3();
230    
231     // existence du tetra
232     MG_NOEUD* noeud;
233     MG_TETRA* tet=mg_maillage->get_mg_tetra(noeud1->get_id(),noeud2->get_id(),noeud3->get_id(),noeud4->get_id());
234     if (tet!=NULL)
235     {
236     clock_t cpu2=clock(); //debug
237     cpu[11]=cpu[11]+(cpu2-cpu1)/CLK_TCK;
238     nbfois[11]++;
239    
240     return NULL;
241     }
242     // validité de la solution
243     TPL_MAP_ENTITE<MG_FRONT_3D*> liste;
244     double x=0.25*(noeud1->get_x()+noeud2->get_x()+noeud3->get_x()+noeud4->get_x());
245     double y=0.25*(noeud1->get_y()+noeud2->get_y()+noeud3->get_y()+noeud4->get_y());
246     double z=0.25*(noeud1->get_z()+noeud2->get_z()+noeud3->get_z()+noeud4->get_z());
247     double longueur=1.2*sqrt((x-noeud1->get_x())*(x-noeud1->get_x())+(y-noeud1->get_y())*(y-noeud1->get_y())+(z-noeud1->get_z())*(z-noeud1->get_z()));
248     octree_de_front->rechercher(x,y,z,longueur,liste);
249     int nb_liste=liste.get_nb();
250     for (int i=0;i<nb_liste;i++)
251     {
252     MG_FRONT_3D* fttmp=liste.get(i);
253     clock_t cpu11=clock(); //debug
254    
255     int inter=inter_tetra_triangle(noeud1,noeud2,noeud3,noeud4,fttmp->get_noeud1(),fttmp->get_noeud2(),fttmp->get_noeud3());
256     clock_t cpu22=clock(); //debug
257     cpu[12]=cpu[12]+(cpu22-cpu11)/CLK_TCK;
258     nbfois[12]++;
259    
260     if (inter==true)
261     liste_intersection.ajouter(fttmp->get_triangle());
262     }
263     int nb_inter=liste_intersection.get_nb();
264     if (nb_inter>0)
265     {
266     clock_t cpu2=clock(); //debug
267     cpu[11]=cpu[11]+(cpu2-cpu1)/CLK_TCK;
268     nbfois[11]++;
269     return NULL;
270     }
271     // evaluation de la situation locale on ne cree pas deux tetra voisin de volume presque nul
272     MG_TRIANGLE* triangle1=mg_maillage->get_mg_triangle(noeud1->get_id(),noeud3->get_id(),noeud2->get_id());
273     MG_TRIANGLE* triangle2=mg_maillage->get_mg_triangle(noeud1->get_id(),noeud2->get_id(),noeud4->get_id());
274     MG_TRIANGLE* triangle3=mg_maillage->get_mg_triangle(noeud2->get_id(),noeud3->get_id(),noeud4->get_id());
275     MG_TRIANGLE* triangle4=mg_maillage->get_mg_triangle(noeud1->get_id(),noeud4->get_id(),noeud3->get_id());
276     double vol1=-1.;
277     double vol2=-1.;
278     double vol3=-1.;
279     double vol4=-1.;
280     if (type==GENERATION_NOEUD)
281     {
282     if (triangle1!=NULL)
283     if (triangle1->get_type_entite()==IDM3D_TRIANGLE)
284     {
285     M3D_TRIANGLE* tri=(M3D_TRIANGLE*)triangle1;
286     MG_TETRA* tet=tri->get_voisin(0);
287     if (tet->get_type_entite()==IDM3D_TETRA)
288     {
289     M3D_TETRA* mtet=(M3D_TETRA*)tet;
290     vol1=mtet->get_volume();
291     }
292     }
293     if (triangle2!=NULL)
294     if (triangle2->get_type_entite()==IDM3D_TRIANGLE)
295     {
296     M3D_TRIANGLE* tri=(M3D_TRIANGLE*)triangle2;
297     MG_TETRA* tet=tri->get_voisin(0);
298     if (tet->get_type_entite()==IDM3D_TETRA)
299     {
300     M3D_TETRA* mtet=(M3D_TETRA*)tet;
301     vol2=mtet->get_volume();
302     }
303     }
304     if (triangle3!=NULL)
305     if (triangle3->get_type_entite()==IDM3D_TRIANGLE)
306     {
307     M3D_TRIANGLE* tri=(M3D_TRIANGLE*)triangle3;
308     MG_TETRA* tet=tri->get_voisin(0);
309     if (tet->get_type_entite()==IDM3D_TETRA)
310     {
311     M3D_TETRA* mtet=(M3D_TETRA*)tet;
312     vol3=mtet->get_volume();
313     }
314     }
315     if (triangle4!=NULL)
316     if (triangle4->get_type_entite()==IDM3D_TRIANGLE)
317     {
318     M3D_TRIANGLE* tri=(M3D_TRIANGLE*)triangle4;
319     MG_TETRA* tet=tri->get_voisin(0);
320     if (tet->get_type_entite()==IDM3D_TETRA)
321     {
322     M3D_TETRA* mtet=(M3D_TETRA*)tet;
323     vol4=mtet->get_volume();
324     }
325     }
326     double vol=calcule_volume(noeud1,noeud2,noeud3,noeud4);
327     // calcul du vol du tetra ideal
328     OT_VECTEUR_3D ab(noeud2->get_x()-noeud1->get_x(),noeud2->get_y()-noeud1->get_y(),noeud2->get_z()-noeud1->get_z());
329     OT_VECTEUR_3D ac(noeud3->get_x()-noeud1->get_x(),noeud3->get_y()-noeud1->get_y(),noeud3->get_z()-noeud1->get_z());
330     OT_VECTEUR_3D da(noeud1->get_x()-noeud4->get_x(),noeud1->get_y()-noeud4->get_y(),noeud1->get_z()-noeud4->get_z());
331     OT_VECTEUR_3D vec=ab&ac;
332     vec.norme();
333     double lg_tet_ideal=fabs(vec*da*ft->get_ifail()*sqrt(3./2.));
334     double vol_ideal=sqrt(2.)/12.*lg_tet_ideal*lg_tet_ideal*lg_tet_ideal;
335     if (vol<0.1*vol_ideal)
336     {
337     if ((vol1>0.) && (vol1<0.1*vol_ideal)) {clock_t cpu2=clock(); //debug
338     cpu[11]=cpu[11]+(cpu2-cpu1)/CLK_TCK;
339     nbfois[11]++;
340     return NULL;}
341     if ((vol2>0.) && (vol2<0.1*vol_ideal)) {clock_t cpu2=clock(); //debug
342     cpu[11]=cpu[11]+(cpu2-cpu1)/CLK_TCK;
343     nbfois[11]++;
344     return NULL;}
345     if ((vol3>0.) && (vol3<0.1*vol_ideal)) {clock_t cpu2=clock(); //debug
346     cpu[11]=cpu[11]+(cpu2-cpu1)/CLK_TCK;
347     nbfois[11]++;
348     return NULL;}
349     if ((vol4>0.) && (vol4<0.1*vol_ideal)) {clock_t cpu2=clock(); //debug
350     cpu[11]=cpu[11]+(cpu2-cpu1)/CLK_TCK;
351     nbfois[11]++;
352     return NULL;}
353     }
354     }
355    
356    
357    
358    
359     // insertion du tetra
360     if (triangle1==NULL) triangle1=insere_triangle(mgvol,noeud1,noeud3,noeud2);
361     if (triangle2==NULL) triangle2=insere_triangle(mgvol,noeud1,noeud2,noeud4);
362     if (triangle3==NULL) triangle3=insere_triangle(mgvol,noeud2,noeud3,noeud4);
363     if (triangle4==NULL) triangle4=insere_triangle(mgvol,noeud1,noeud4,noeud3);
364     MG_TETRA* tetra=cree_tetra(mgvol,noeud1,noeud2,noeud3,noeud4,triangle1,triangle2,triangle3,triangle4);
365     mg_maillage->ajouter_mg_tetra(tetra);
366     if (triangle1->get_type_entite()==IDM3D_TRIANGLE)
367     {
368     M3D_TRIANGLE* tri=(M3D_TRIANGLE*)triangle1;
369     tri->ajouter_voisin(tetra);
370     }
371     if (triangle2->get_type_entite()==IDM3D_TRIANGLE)
372     {
373     M3D_TRIANGLE* tri=(M3D_TRIANGLE*)triangle2;
374     tri->ajouter_voisin(tetra);
375     }
376     if (triangle3->get_type_entite()==IDM3D_TRIANGLE)
377     {
378     M3D_TRIANGLE* tri=(M3D_TRIANGLE*)triangle3;
379     tri->ajouter_voisin(tetra);
380     }
381     if (triangle4->get_type_entite()==IDM3D_TRIANGLE)
382     {
383     M3D_TRIANGLE* tri=(M3D_TRIANGLE*)triangle4;
384     tri->ajouter_voisin(tetra);
385     }
386     clock_t cpu2=clock(); //debug
387     cpu[11]=cpu[11]+(cpu2-cpu1)/CLK_TCK;
388     nbfois[11]++;
389    
390     return tetra;
391    
392     }
393    
394    
395    
396     class MG_TRIANGLE* MAILLEUR3D::insere_triangle(class MG_VOLUME* mgvol,MG_NOEUD* noeud1,MG_NOEUD* noeud2,MG_NOEUD* noeud3)
397     {
398     clock_t cpu1=clock(); //debug
399     // existence du triangle
400     MG_TRIANGLE* triangle=mg_maillage->get_mg_triangle(noeud1->get_id(),noeud3->get_id(),noeud2->get_id());
401     if (triangle!=NULL) return NULL;
402     // insertion du triangle
403     MG_SEGMENT* segment1=mg_maillage->get_mg_segment(noeud1->get_id(),noeud2->get_id());
404     MG_SEGMENT* segment2=mg_maillage->get_mg_segment(noeud1->get_id(),noeud3->get_id());
405     MG_SEGMENT* segment3=mg_maillage->get_mg_segment(noeud2->get_id(),noeud3->get_id());
406     if (segment1==NULL) segment1=cree_segment(mgvol,noeud1,noeud2);
407     if (segment2==NULL) segment2=cree_segment(mgvol,noeud1,noeud3);
408     if (segment3==NULL) segment3=cree_segment(mgvol,noeud2,noeud3);
409     MG_TRIANGLE* tri=cree_triangle(mgvol,noeud1,noeud2,noeud3,segment1,segment2,segment3);
410     clock_t cpu2=clock(); //debug
411     cpu[6]=cpu[6]+(cpu2-cpu1)/CLK_TCK;
412     nbfois[6]++;
413     return tri;
414     }
415    
416    
417    
418     class MG_SEGMENT* MAILLEUR3D::insere_segment(class MG_VOLUME* mgvol,MG_NOEUD* noeud1,MG_NOEUD* noeud2)
419     {
420     clock_t cpu1=clock(); //debug
421     //existence du segment
422     MG_SEGMENT* segment=mg_maillage->get_mg_segment(noeud1->get_id(),noeud2->get_id());
423     if (segment!=NULL) return NULL;
424     // insertion du segment
425     MG_SEGMENT* seg=cree_segment(mgvol,noeud1,noeud2);
426     clock_t cpu2=clock(); //debug
427     cpu[7]=cpu[7]+(cpu2-cpu1)/CLK_TCK;
428     nbfois[7]++;
429     return seg;
430     }
431    
432    
433    
434    
435     MG_TETRA* MAILLEUR3D::cree_tetra(class MG_VOLUME* mgvol,MG_NOEUD* noeud1,MG_NOEUD* noeud2,MG_NOEUD* noeud3,MG_NOEUD* noeud4,MG_TRIANGLE* triangle1,MG_TRIANGLE* triangle2,MG_TRIANGLE* triangle3,MG_TRIANGLE* triangle4)
436     {
437     M3D_TETRA* tetra=new M3D_TETRA(mgvol,noeud1,noeud2,noeud3,noeud4,triangle1,triangle2,triangle3,triangle4);
438     mg_maillage->ajouter_mg_tetra(tetra);
439     return tetra;
440     }
441    
442    
443     MG_TRIANGLE* MAILLEUR3D::cree_triangle(class MG_VOLUME* mgvol,MG_NOEUD* noeud1,MG_NOEUD* noeud2,MG_NOEUD* noeud3,MG_SEGMENT* segment1,MG_SEGMENT* segment2,MG_SEGMENT* segment3)
444     {
445     M3D_TRIANGLE* tri=new M3D_TRIANGLE(mgvol,noeud1,noeud2,noeud3,segment1,segment2,segment3);
446     mg_maillage->ajouter_mg_triangle(tri);
447     return tri;
448     }
449    
450    
451    
452     MG_SEGMENT* MAILLEUR3D::cree_segment(class MG_VOLUME* mgvol,MG_NOEUD* noeud1,MG_NOEUD* noeud2)
453     {
454     MG_SEGMENT* seg=mg_maillage->ajouter_mg_segment(mgvol,noeud1,noeud2);
455     return seg;
456     }
457    
458    
459    
460     double MAILLEUR3D::calcule_volume(MG_NOEUD* noeud1,MG_NOEUD* noeud2,MG_NOEUD* noeud3,MG_NOEUD* noeud4)
461     {
462     OT_VECTEUR_3D ab(noeud2->get_x()-noeud1->get_x(),noeud2->get_y()-noeud1->get_y(),noeud2->get_z()-noeud1->get_z());
463     OT_VECTEUR_3D ac(noeud3->get_x()-noeud1->get_x(),noeud3->get_y()-noeud1->get_y(),noeud3->get_z()-noeud1->get_z());
464     OT_VECTEUR_3D da(noeud1->get_x()-noeud4->get_x(),noeud1->get_y()-noeud4->get_y(),noeud1->get_z()-noeud4->get_z());
465    
466     OT_VECTEUR_3D pvec=ab&ac;
467     double volume=pvec*da;
468     volume=fabs(volume)/3.;
469     return volume;
470     }
471    
472    
473    
474     double MAILLEUR3D::calcule_longueur_caracteristique(MG_NOEUD* noeud1,MG_NOEUD* noeud2,MG_NOEUD* noeud3,MG_NOEUD* noeud4)
475     {
476    
477     OT_VECTEUR_3D ab(noeud2->get_x()-noeud1->get_x(),noeud2->get_y()-noeud1->get_y(),noeud2->get_z()-noeud1->get_z());
478     double lg=ab.get_longueur();
479     int nb=1;
480     if (noeud3!=NULL)
481     {
482     OT_VECTEUR_3D ac(noeud3->get_x()-noeud1->get_x(),noeud3->get_y()-noeud1->get_y(),noeud3->get_z()-noeud1->get_z());
483     OT_VECTEUR_3D bc(noeud3->get_x()-noeud2->get_x(),noeud3->get_y()-noeud2->get_y(),noeud3->get_z()-noeud2->get_z());
484     lg=lg+ac.get_longueur()+bc.get_longueur();
485     nb=3;
486     }
487     if (noeud4!=NULL)
488     {
489     OT_VECTEUR_3D da(noeud1->get_x()-noeud4->get_x(),noeud1->get_y()-noeud4->get_y(),noeud1->get_z()-noeud4->get_z());
490     OT_VECTEUR_3D bd(noeud4->get_x()-noeud2->get_x(),noeud4->get_y()-noeud2->get_y(),noeud4->get_z()-noeud2->get_z());
491     OT_VECTEUR_3D cd(noeud4->get_x()-noeud3->get_x(),noeud4->get_y()-noeud3->get_y(),noeud4->get_z()-noeud3->get_z());
492     lg=lg+da.get_longueur()+bd.get_longueur()+cd.get_longueur();
493     nb=6;
494     }
495     lg=lg*1./nb;
496     return lg;
497    
498     }
499    
500    
501    
502    
503     int MAILLEUR3D::genere_noeud(MG_VOLUME* mgvol,MG_FRONT_3D* ft,TPL_MAP_ENTITE<MG_NOEUD*> &liste_noeud)
504     {
505     clock_t cpu1=clock(); //debug
506     MG_NOEUD* noeud1=ft->get_noeud1();
507     MG_NOEUD* noeud2=ft->get_noeud2();
508     MG_NOEUD* noeud3=ft->get_noeud3();
509     OT_VECTEUR_3D n1n2(noeud2->get_x()-noeud1->get_x(),noeud2->get_y()-noeud1->get_y(),noeud2->get_z()-noeud1->get_z());
510     OT_VECTEUR_3D n1n3(noeud3->get_x()-noeud1->get_x(),noeud3->get_y()-noeud1->get_y(),noeud3->get_z()-noeud1->get_z());
511     OT_VECTEUR_3D normal=n1n2&n1n3;
512     normal.norme();
513     OT_VECTEUR_3D nmilieu(noeud1->get_x()+noeud2->get_x(),noeud1->get_y()+noeud2->get_y(),noeud1->get_z()+noeud2->get_z());
514     nmilieu=0.5*nmilieu;
515     OT_VECTEUR_3D mediane(noeud3->get_x()-nmilieu.get_x(),noeud3->get_y()-nmilieu.get_y(),noeud3->get_z()-nmilieu.get_z());
516     double d1a=ft->get_triangle()->get_longueur();
517     double dens1,dens2,dens3;
518     double xyz1[3]={noeud1->get_x(),noeud1->get_y(),noeud1->get_z()};
519     double xyz2[3]={noeud2->get_x(),noeud2->get_y(),noeud2->get_z()};
520     double xyz3[3]={noeud3->get_x(),noeud3->get_y(),noeud3->get_z()};
521     metrique->evaluer(xyz1,&dens1);
522     metrique->evaluer(xyz2,&dens2);
523     metrique->evaluer(xyz3,&dens3);
524     double d1b=0.3333333333*(dens1+dens2+dens3);
525     double d1=0.35*d1a+0.65*d1b;
526     double d2=d1*d1-4./9.*mediane.get_longueur()*mediane.get_longueur();
527     if (d2<0.) d2=d1b;
528     d2=sqrt(d2);
529     d2=d2/ft->get_ifail();
530     double x=0.333333333*(xyz1[0]+xyz2[0]+xyz3[0])+d2*normal.get_x();
531     double y=0.333333333*(xyz1[1]+xyz2[1]+xyz3[1])+d2*normal.get_y();
532     double z=0.333333333*(xyz1[2]+xyz2[2]+xyz3[2])+d2*normal.get_z();
533    
534     // recherche des voisins
535     TPL_MAP_ENTITE<MG_FRONT_3D*> liste_front;
536     TPL_MAP_ENTITE<MG_FRONT_3D*> liste_front_proche;
537     octree_de_front->rechercher(x,y,z,d2*1.2,liste_front);
538     int nb_liste=liste_front.get_nb();
539     double disref=100.*d2;
540     MG_NOEUD* noeudref=NULL;
541     for (int i=0;i<nb_liste;i++)
542     {
543     MG_FRONT_3D* ft2=liste_front.get(i);
544     for (int j=0;j<3;j++)
545     {
546     MG_NOEUD* noeud_tmp;
547     if (j==0) noeud_tmp=ft2->get_noeud1();
548     if (j==1) noeud_tmp=ft2->get_noeud2();
549     if (j==2) noeud_tmp=ft2->get_noeud3();
550     OT_VECTEUR_3D vec(noeud_tmp->get_x()-x,noeud_tmp->get_y()-y,noeud_tmp->get_z()-z);
551     double dis=vec.get_longueur();
552     if (dis<disref)
553     {
554     disref=dis;
555     noeudref=noeud_tmp;
556     }
557     }
558     if (ft2!=ft)
559     {
560     double dismin=evaluer_distance_noeud_triangle(x,y,z,ft2->get_triangle());
561     if (dismin<d2) liste_front_proche.ajouter(ft2);
562     }
563     }
564     // il y a un noeud proche
565     if (noeudref!=NULL)
566     if (disref<d2)
567     if (noeudref!=noeud1)
568     if (noeudref!=noeud2)
569     if (noeudref!=noeud3)
570     {
571     liste_noeud.ajouter(noeudref);
572     clock_t cpu2=clock(); //debug
573     cpu[8]=cpu[8]+(cpu2-cpu1)/CLK_TCK;
574     nbfois[8]++;
575     return EXISTE;
576     }
577     // il n y a pas de noeud proche mais il y a une face proche
578     int nb_proche=liste_front_proche.get_nb();
579     if (nb_proche>0)
580     {
581     for (int i=0;i<nb_proche;i++)
582     {
583     MG_FRONT_3D* fttmp=liste_front_proche.get(i);
584     MG_NOEUD* noeudtmp=fttmp->get_noeud1();
585     if (noeudtmp!=noeud1)
586     if (noeudtmp!=noeud2)
587     if (noeudtmp!=noeud3) liste_noeud.ajouter(noeudtmp);
588     noeudtmp=fttmp->get_noeud2();
589     if (noeudtmp!=noeud1)
590     if (noeudtmp!=noeud2)
591     if (noeudtmp!=noeud3) liste_noeud.ajouter(noeudtmp);
592     noeudtmp=fttmp->get_noeud3();
593     if (noeudtmp!=noeud1)
594     if (noeudtmp!=noeud2)
595     if (noeudtmp!=noeud3) liste_noeud.ajouter(noeudtmp);
596    
597     }
598     clock_t cpu2=clock(); //debug
599     cpu[8]=cpu[8]+(cpu2-cpu1)/CLK_TCK;
600     nbfois[8]++;
601     return MULTIEXISTE;
602     }
603     // il faut creer un noeud c est parfait c est qu on veut.
604     M3D_NOEUD* nouveau_noeud=new M3D_NOEUD(mgvol,x,y,z);
605     mg_maillage->ajouter_mg_noeud(nouveau_noeud);
606     liste_noeud.ajouter(nouveau_noeud);
607     clock_t cpu2=clock(); //debug
608     cpu[8]=cpu[8]+(cpu2-cpu1)/CLK_TCK;
609     nbfois[8]++;
610     return CREATION;
611     }
612    
613    
614    
615     double MAILLEUR3D::evaluer_distance_noeud_triangle(double x,double y,double z,MG_TRIANGLE* tri)
616     {
617     double eps=0.0001;
618     MG_NOEUD* noeud1=tri->get_noeud1();
619     MG_NOEUD* noeud2=tri->get_noeud2();
620     MG_NOEUD* noeud3=tri->get_noeud3();
621     OT_VECTEUR_3D n1n2(noeud2->get_x()-noeud1->get_x(),noeud2->get_y()-noeud1->get_y(),noeud2->get_z()-noeud1->get_z());
622     OT_VECTEUR_3D n1n3(noeud3->get_x()-noeud1->get_x(),noeud3->get_y()-noeud1->get_y(),noeud3->get_z()-noeud1->get_z());
623     OT_VECTEUR_3D n1n4(x-noeud1->get_x(),y-noeud1->get_y(),z-noeud1->get_z());
624     OT_VECTEUR_3D n1n=n1n2&n1n3;
625     OT_MATRICE_3D matbase(n1n2,n1n3,n1n);
626     OT_MATRICE_3D matxsi(n1n4,n1n3,n1n);
627     OT_MATRICE_3D mateta(n1n2,n1n4,n1n);
628     double xsi=matxsi.get_determinant()/matbase.get_determinant();
629     double eta=mateta.get_determinant()/matbase.get_determinant();
630     double dis;
631     if ((eta>-eps) && (xsi>-eps) && (eta+xsi<1.0+eps))
632     {
633     double xproj=(1-xsi-eta)*noeud1->get_x()+xsi*noeud2->get_x()+eta*noeud3->get_x();
634     double yproj=(1-xsi-eta)*noeud1->get_y()+xsi*noeud2->get_y()+eta*noeud3->get_y();
635     double zproj=(1-xsi-eta)*noeud1->get_z()+xsi*noeud2->get_z()+eta*noeud3->get_z();
636     OT_VECTEUR_3D vec(x-xproj,y-yproj,z-zproj);
637     dis=vec.get_longueur();
638     }
639     else
640     {
641     OT_VECTEUR_3D nn1(x-noeud1->get_x(),y-noeud1->get_y(),z-noeud1->get_z());
642     OT_VECTEUR_3D nn2(x-noeud2->get_x(),y-noeud2->get_y(),z-noeud2->get_z());
643     OT_VECTEUR_3D nn3(x-noeud3->get_x(),y-noeud3->get_y(),z-noeud3->get_z());
644     dis=0.3333333333333*(nn1.get_longueur()+nn2.get_longueur()+nn3.get_longueur());
645     }
646     return dis;
647     }
648    
649    
650     void MAILLEUR3D::detruit_tetra(MG_TRIANGLE* tri)
651     {
652     clock_t cpu1=clock(); //debug
653    
654     M3D_TRIANGLE* mtri;
655     //if (tri->get_type_entite()==IDMG_TRIANGLE) return;
656     //else mtri=(M3D_TRIANGLE*)tri;
657     mtri=(M3D_TRIANGLE*)tri;
658    
659     M3D_TETRA* tet=(M3D_TETRA*)mtri->get_voisin(0);
660     tet->etat_detruit();
661    
662     TPL_MAP_ENTITE<MG_FRONT_3D*> liste_de_front;
663     TPL_MAP_ENTITE<MG_FRONT_3D*> liste_de_front2;
664     MG_FRONT_3D* fttet1=NULL;
665     MG_FRONT_3D* fttet2=NULL;
666     MG_FRONT_3D* fttet3=NULL;
667     MG_FRONT_3D* fttet4=NULL;
668    
669     MG_TRIANGLE* triangle1=tet->get_triangle1();
670     M3D_TRIANGLE* mtriangle=(M3D_TRIANGLE*)triangle1;
671     if (mtriangle->get_etat_front()==FRONT_ACTIF)
672     {
673     MG_FRONT_3D* fttmp=mtriangle->get_mgfront();
674     fttet1=fttmp;
675     liste_de_front.ajouter(fttmp->get_front_voisin(0));
676     liste_de_front.ajouter(fttmp->get_front_voisin(1));
677     liste_de_front.ajouter(fttmp->get_front_voisin(2));
678     //supprimer_front_en_reculant(fttmp);
679     mtriangle->supprimer_voisin(tet);
680     }
681     else
682     {
683     MG_FRONT_3D* ftn1=ajouter_front_courant(PRIORITAIRE,tet->get_noeud1(),tet->get_noeud2(),tet->get_noeud3(),triangle1);
684     ftn1->incremente_ifail();
685     mtriangle->supprimer_voisin(tet);
686     liste_de_front2.ajouter(ftn1);
687     }
688     MG_TRIANGLE* triangle2=tet->get_triangle2();
689     mtriangle=(M3D_TRIANGLE*)triangle2;
690     if (mtriangle->get_etat_front()==FRONT_ACTIF)
691     {
692     MG_FRONT_3D* fttmp=mtriangle->get_mgfront();
693     fttet2=fttmp;
694     liste_de_front.ajouter(fttmp->get_front_voisin(0));
695     liste_de_front.ajouter(fttmp->get_front_voisin(1));
696     liste_de_front.ajouter(fttmp->get_front_voisin(2));
697     //supprimer_front_en_reculant(fttmp);
698     mtriangle->supprimer_voisin(tet);
699     }
700     else
701     {
702     MG_FRONT_3D* ftn1=ajouter_front_courant(PRIORITAIRE,tet->get_noeud1(),tet->get_noeud4(),tet->get_noeud2(),triangle2);
703     ftn1->incremente_ifail();
704     mtriangle->supprimer_voisin(tet);
705     liste_de_front2.ajouter(ftn1);
706     }
707     MG_TRIANGLE* triangle3=tet->get_triangle3();
708     mtriangle=(M3D_TRIANGLE*)triangle3;
709     if (mtriangle->get_etat_front()==FRONT_ACTIF)
710     {
711     MG_FRONT_3D* fttmp=mtriangle->get_mgfront();
712     fttet3=fttmp;
713     liste_de_front.ajouter(fttmp->get_front_voisin(0));
714     liste_de_front.ajouter(fttmp->get_front_voisin(1));
715     liste_de_front.ajouter(fttmp->get_front_voisin(2));
716     //supprimer_front_en_reculant(fttmp);
717     mtriangle->supprimer_voisin(tet);
718     }
719     else
720     {
721     MG_FRONT_3D* ftn1=ajouter_front_courant(PRIORITAIRE,tet->get_noeud2(),tet->get_noeud4(),tet->get_noeud3(),triangle3);
722     ftn1->incremente_ifail();
723     mtriangle->supprimer_voisin(tet);
724     liste_de_front2.ajouter(ftn1);
725     }
726     MG_TRIANGLE* triangle4=tet->get_triangle4();
727     mtriangle=(M3D_TRIANGLE*)triangle4;
728     if (mtriangle->get_etat_front()==FRONT_ACTIF)
729     {
730     MG_FRONT_3D* fttmp=mtriangle->get_mgfront();
731     fttet4=fttmp;
732     liste_de_front.ajouter(fttmp->get_front_voisin(0));
733     liste_de_front.ajouter(fttmp->get_front_voisin(1));
734     liste_de_front.ajouter(fttmp->get_front_voisin(2));
735     //supprimer_front_en_reculant(fttmp);
736     mtriangle->supprimer_voisin(tet);
737     }
738     else
739     {
740     MG_FRONT_3D* ftn1=ajouter_front_courant(PRIORITAIRE,tet->get_noeud1(),tet->get_noeud3(),tet->get_noeud4(),triangle4);
741     ftn1->incremente_ifail();
742     mtriangle->supprimer_voisin(tet);
743     liste_de_front2.ajouter(ftn1);
744     }
745    
746     if (fttet1!=NULL)
747     {
748     liste_de_front.supprimer(fttet1);
749     supprimer_front_en_reculant(fttet1);
750     }
751     if (fttet2!=NULL)
752     {
753     liste_de_front.supprimer(fttet2);
754     supprimer_front_en_reculant(fttet2);
755     }
756     if (fttet3!=NULL)
757     {
758     liste_de_front.supprimer(fttet3);
759     supprimer_front_en_reculant(fttet3);
760     }
761     if (fttet4!=NULL)
762     {
763     liste_de_front.supprimer(fttet4);
764     supprimer_front_en_reculant(fttet4);
765     }
766     int nb_liste1=liste_de_front.get_nb();
767     for (int i=0;i<nb_liste1;i++)
768     if (liste_de_front.get(i)!=NULL)
769     mise_a_jour_voisin_front(liste_de_front.get(i));
770     int nb_liste2=liste_de_front2.get_nb();
771     for (int i=0;i<nb_liste2;i++)
772     {
773     mise_a_jour_voisin_front(liste_de_front2.get(i));
774     mise_a_jour_voisin_front(liste_de_front2.get(i)->get_front_voisin(0));
775     mise_a_jour_voisin_front(liste_de_front2.get(i)->get_front_voisin(1));
776     mise_a_jour_voisin_front(liste_de_front2.get(i)->get_front_voisin(2));
777     }
778    
779     for (int k=0;k<4;k++)
780     {
781     MG_NOEUD* noeud;
782     if (k==0) noeud=tet->get_noeud1();
783     if (k==1) noeud=tet->get_noeud2();
784     if (k==2) noeud=tet->get_noeud3();
785     if (k==3) noeud=tet->get_noeud4();
786     int nb_tet=noeud->get_lien_tetra()->get_nb();
787    
788     if (noeud->get_type_entite()==IDM3D_NOEUD)
789     {
790     int test=false;
791     M3D_NOEUD* mnoeud=(M3D_NOEUD*)noeud;
792     for (int i=0;i<nb_tet;i++)
793     {
794     MG_TETRA* tet=noeud->get_lien_tetra()->get(i);
795     if (tet->get_type_entite()==IDM3D_TETRA)
796     {
797     M3D_TETRA* mtet=(M3D_TETRA*)tet;
798     if (mtet->get_etat()==ACTIF) test=true;
799     }
800     else test=true;
801     }
802     if (test==false) mnoeud->etat_detruit();
803     }
804     }
805     clock_t cpu2=clock();
806     cpu[9]=cpu[9]+(cpu2-cpu1)/CLK_TCK;
807     nbfois[9]++;
808    
809     }
810    
811    
812    
813     void MAILLEUR3D::detruit_element_inactif(void)
814     {
815     clock_t cpu1=clock();
816     int nb_tetra=mg_maillage->get_nb_mg_tetra();
817     int num=0;
818     for (int i=0;i<nb_tetra;i++)
819     {
820     MG_TETRA* tet=mg_maillage->get_mg_tetra(num);
821     M3D_TETRA* mtet=(M3D_TETRA*)tet;
822     if (mtet->get_etat()==DETRUIT)
823     mg_maillage->supprimer_mg_tetra(num);
824     else num++;
825     }
826     int nb_noeud=mg_maillage->get_nb_mg_noeud();
827     num=0;
828     for (int i=0;i<nb_noeud;i++)
829     {
830     MG_NOEUD* noeud=mg_maillage->get_mg_noeud(num);
831     if (noeud->get_type_entite()==IDM3D_NOEUD)
832     {
833     M3D_NOEUD* mnoeud=(M3D_NOEUD*)noeud;
834     if (mnoeud->get_etat()==DETRUIT)
835     mg_maillage->supprimer_mg_noeud(num);
836     else num++;
837     }
838     else num++;
839     }
840    
841     clock_t cpu2=clock();
842     cpu[10]=cpu[10]+(cpu2-cpu1)/CLK_TCK;
843     nbfois[10]++;
844     }
845