ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/mailleur3d_outil.cpp
Revision: 253
Committed: Tue Jul 13 19:40:46 2010 UTC (14 years, 10 months ago) by francois
File size: 57434 byte(s)
Log Message:
changement de hiearchie et utilisation de ccmake + mise a jour

File Contents

# User Rev Content
1 5 //------------------------------------------------------------
2     //------------------------------------------------------------
3     // MAGiC
4 francois 129 // Jean Christophe Cuilli�re et Vincent FRANCOIS
5     // D�partement de G�nie M�canique - UQTR
6 5 //------------------------------------------------------------
7 francois 129 // 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 5 // des auteurs (contact : francois@uqtr.ca)
12     //------------------------------------------------------------
13     //------------------------------------------------------------
14     //
15     // mailleur3d_outil.cpp
16     //
17     //------------------------------------------------------------
18     //------------------------------------------------------------
19     // COPYRIGHT 2000
20 francois 129 // Version du 02/03/2006 � 11H23
21 5 //------------------------------------------------------------
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 francois 35 #include "mailleur3d_front.h"
34 5 #include "ot_mathematique.h"
35 francois 102 #include "ot_boite_3D.h"
36 5 #include "m3d_noeud.h"
37     #include "m3d_tetra.h"
38     #include "m3d_triangle.h"
39     #include "const.h"
40    
41     //---------------------------------------------------------------------------
42    
43     //#pragma package(smart_init)
44    
45    
46    
47    
48    
49    
50    
51     double MAILLEUR3D::angle_front(MG_FRONT_3D* ft1,MG_FRONT_3D* ft2)
52     {
53     MG_NOEUD* noeud1=ft1->get_noeud1();
54     MG_NOEUD* noeud2=ft1->get_noeud2();
55     MG_NOEUD* noeud3=ft1->get_noeud3();
56     MG_NOEUD* noeud4=ft2->get_noeud1();
57     MG_NOEUD* noeud5=ft2->get_noeud2();
58     MG_NOEUD* noeud6=ft2->get_noeud3();
59    
60 francois 54 MG_NOEUD* noeuda=NULL;
61     MG_NOEUD* noeudb=NULL;
62     MG_NOEUD* noeudc=NULL;
63     MG_NOEUD* noeudd=NULL;
64 5
65    
66     if (!(noeud1==noeud4))
67     if (!(noeud1==noeud5))
68     if (!(noeud1==noeud6))
69     {
70     noeuda=noeud1;
71     noeudb=noeud2;
72     noeudc=noeud3;
73     }
74     if (!(noeud2==noeud4))
75     if (!(noeud2==noeud5))
76     if (!(noeud2==noeud6))
77     {
78     noeuda=noeud2;
79     noeudb=noeud1;
80     noeudc=noeud3;
81     }
82     if (!(noeud3==noeud4))
83     if (!(noeud3==noeud5))
84     if (!(noeud3==noeud6))
85     {
86     noeuda=noeud3;
87     noeudb=noeud1;
88     noeudc=noeud2;
89     }
90 francois 54 if (noeuda==NULL) return 2*M_PI;
91 5 if (!(noeud4==noeudb))
92     if (!(noeud4==noeudc)) noeudd=noeud4;
93     if (!(noeud5==noeudb))
94     if (!(noeud5==noeudc)) noeudd=noeud5;
95     if (!(noeud6==noeudb))
96     if (!(noeud6==noeudc)) noeudd=noeud6;
97    
98     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()));
99     OT_VECTEUR_3D bc(noeudc->get_x()-noeudb->get_x(),noeudc->get_y()-noeudb->get_y(),noeudc->get_z()-noeudb->get_z());
100     bc.norme();
101     double xyz[3];
102     xyz[0]=noeuda->get_x();
103     xyz[1]=noeuda->get_y();
104     xyz[2]=noeuda->get_z();
105     OT_VECTEUR_3D n1(xyz,milieu.get_xyz());
106     n1.norme();
107     OT_VECTEUR_3D n=bc&n1;
108     n1=n&bc;
109     n1.norme();
110     xyz[0]=noeudd->get_x();
111     xyz[1]=noeudd->get_y();
112     xyz[2]=noeudd->get_z();
113     OT_VECTEUR_3D n2(xyz,milieu.get_xyz());
114     n2.norme();
115     n=bc&n2;
116     n2=n&bc;
117     n2.norme();
118    
119     double ps=n1*n2;
120     if (ps>1.) ps=1.;
121     if (ps<-1.) ps=-1.;
122     double angle=acos(ps);
123     OT_VECTEUR_3D n1n2(noeud2->get_x()-noeud1->get_x(),noeud2->get_y()-noeud1->get_y(),noeud2->get_z()-noeud1->get_z());
124     OT_VECTEUR_3D n1n3(noeud3->get_x()-noeud1->get_x(),noeud3->get_y()-noeud1->get_y(),noeud3->get_z()-noeud1->get_z());
125     OT_VECTEUR_3D nbnd(noeudd->get_x()-noeudb->get_x(),noeudd->get_y()-noeudb->get_y(),noeudd->get_z()-noeudb->get_z());
126     n=n1n2&n1n3;
127     n.norme();
128     nbnd.norme();
129    
130     ps=n*nbnd;
131     if (ps<-0.0001) angle=2*M_PI-angle;
132     return angle;
133    
134    
135     }
136    
137    
138    
139     double MAILLEUR3D::angle_front(MG_FRONT_3D* ft,MG_SEGMENT* seg)
140     {
141     MG_NOEUD* noeud1;
142     MG_NOEUD* noeud2;
143     MG_NOEUD* noeud3;
144     MG_NOEUD* noeud4;
145    
146    
147     if (seg->get_noeud1()==ft->get_noeud1())
148     {
149     noeud1=ft->get_noeud1();
150     noeud2=ft->get_noeud2();
151     noeud3=ft->get_noeud3();
152     noeud4=seg->get_noeud2();
153     }
154     if (seg->get_noeud1()==ft->get_noeud2())
155     {
156     noeud1=ft->get_noeud2();
157     noeud2=ft->get_noeud3();
158     noeud3=ft->get_noeud1();
159     noeud4=seg->get_noeud2();
160     }
161     if (seg->get_noeud1()==ft->get_noeud3())
162     {
163     noeud1=ft->get_noeud3();
164     noeud2=ft->get_noeud1();
165     noeud3=ft->get_noeud2();
166     noeud4=seg->get_noeud2();
167     }
168     if (seg->get_noeud2()==ft->get_noeud1())
169     {
170     noeud1=ft->get_noeud1();
171     noeud2=ft->get_noeud2();
172     noeud3=ft->get_noeud3();
173     noeud4=seg->get_noeud1();
174     }
175     if (seg->get_noeud2()==ft->get_noeud2())
176     {
177     noeud1=ft->get_noeud2();
178     noeud2=ft->get_noeud3();
179     noeud3=ft->get_noeud1();
180     noeud4=seg->get_noeud1();
181     }
182     if (seg->get_noeud2()==ft->get_noeud3())
183     {
184     noeud1=ft->get_noeud3();
185     noeud2=ft->get_noeud1();
186     noeud3=ft->get_noeud2();
187     noeud4=seg->get_noeud1();
188     }
189    
190     OT_VECTEUR_3D n1n2(noeud2->get_x()-noeud1->get_x(),noeud2->get_y()-noeud1->get_y(),noeud2->get_z()-noeud1->get_z());
191     OT_VECTEUR_3D n1n3(noeud3->get_x()-noeud1->get_x(),noeud3->get_y()-noeud1->get_y(),noeud3->get_z()-noeud1->get_z());
192     OT_VECTEUR_3D n1n4(noeud4->get_x()-noeud1->get_x(),noeud4->get_y()-noeud1->get_y(),noeud4->get_z()-noeud1->get_z());
193     OT_VECTEUR_3D n1n=n1n2&n1n3;
194     OT_MATRICE_3D mat(n1n2,n1n3,n1n);
195     OT_MATRICE_3D matx(n1n4,n1n3,n1n);
196     OT_MATRICE_3D maty(n1n2,n1n4,n1n);
197     OT_MATRICE_3D matz(n1n2,n1n3,n1n4);
198     double det=mat.get_determinant();
199     double xsi=matx.get_determinant()/det;
200     double eta=maty.get_determinant()/det;
201     double psi=matz.get_determinant()/det;
202    
203     double proj[3];
204     proj[0]=(1.-xsi-eta)*noeud1->get_x()+xsi*noeud2->get_x()+eta*noeud3->get_x();
205     proj[1]=(1.-xsi-eta)*noeud1->get_y()+xsi*noeud2->get_y()+eta*noeud3->get_y();
206     proj[2]=(1.-xsi-eta)*noeud1->get_z()+xsi*noeud2->get_z()+eta*noeud3->get_z();
207     OT_VECTEUR_3D n1nproj(noeud1->get_x()-proj[0],noeud1->get_y()-proj[1],noeud1->get_z()-proj[2]);
208     //n1n4.norme();
209     //n1nproj.norme();
210     double ps=n1n4*n1nproj;
211     if (!(OPERATEUR::egal(ps,0.0,0.0001))) ps=ps/n1n4.get_longueur()/n1nproj.get_longueur(); else ps=0.;
212     if (ps>1.) ps=1.;
213     if (ps<-1.) ps=-1.;
214     double angle=acos(ps);
215     if ((xsi<0.) || (eta<0.) ) angle=angle+M_PI/2.;
216     if (psi<0.) angle=2*M_PI-angle;
217     return angle;
218    
219     }
220    
221    
222    
223    
224    
225     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)
226     {
227     clock_t cpu1=clock(); //debug
228     liste_intersection.vide();
229     MG_NOEUD* noeud1=ft->get_noeud1();
230     MG_NOEUD* noeud2=ft->get_noeud2();
231     MG_NOEUD* noeud3=ft->get_noeud3();
232    
233     // existence du tetra
234     MG_NOEUD* noeud;
235     MG_TETRA* tet=mg_maillage->get_mg_tetra(noeud1->get_id(),noeud2->get_id(),noeud3->get_id(),noeud4->get_id());
236     if (tet!=NULL)
237     {
238     clock_t cpu2=clock(); //debug
239     cpu[11]=cpu[11]+(cpu2-cpu1)/CLK_TCK;
240     nbfois[11]++;
241    
242     return NULL;
243     }
244 francois 129 // validit� de la solution
245 5 TPL_MAP_ENTITE<MG_FRONT_3D*> liste;
246     double x=0.25*(noeud1->get_x()+noeud2->get_x()+noeud3->get_x()+noeud4->get_x());
247     double y=0.25*(noeud1->get_y()+noeud2->get_y()+noeud3->get_y()+noeud4->get_y());
248     double z=0.25*(noeud1->get_z()+noeud2->get_z()+noeud3->get_z()+noeud4->get_z());
249     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()));
250     octree_de_front->rechercher(x,y,z,longueur,liste);
251     int nb_liste=liste.get_nb();
252     for (int i=0;i<nb_liste;i++)
253     {
254     MG_FRONT_3D* fttmp=liste.get(i);
255     clock_t cpu11=clock(); //debug
256    
257     int inter=inter_tetra_triangle(noeud1,noeud2,noeud3,noeud4,fttmp->get_noeud1(),fttmp->get_noeud2(),fttmp->get_noeud3());
258     clock_t cpu22=clock(); //debug
259     cpu[12]=cpu[12]+(cpu22-cpu11)/CLK_TCK;
260     nbfois[12]++;
261    
262     if (inter==true)
263     liste_intersection.ajouter(fttmp->get_triangle());
264     }
265     int nb_inter=liste_intersection.get_nb();
266     if (nb_inter>0)
267     {
268     clock_t cpu2=clock(); //debug
269     cpu[11]=cpu[11]+(cpu2-cpu1)/CLK_TCK;
270     nbfois[11]++;
271     return NULL;
272     }
273     // evaluation de la situation locale on ne cree pas deux tetra voisin de volume presque nul
274     MG_TRIANGLE* triangle1=mg_maillage->get_mg_triangle(noeud1->get_id(),noeud3->get_id(),noeud2->get_id());
275     MG_TRIANGLE* triangle2=mg_maillage->get_mg_triangle(noeud1->get_id(),noeud2->get_id(),noeud4->get_id());
276     MG_TRIANGLE* triangle3=mg_maillage->get_mg_triangle(noeud2->get_id(),noeud3->get_id(),noeud4->get_id());
277     MG_TRIANGLE* triangle4=mg_maillage->get_mg_triangle(noeud1->get_id(),noeud4->get_id(),noeud3->get_id());
278     double vol1=-1.;
279     double vol2=-1.;
280     double vol3=-1.;
281     double vol4=-1.;
282     if (type==GENERATION_NOEUD)
283     {
284     if (triangle1!=NULL)
285     if (triangle1->get_type_entite()==IDM3D_TRIANGLE)
286     {
287     M3D_TRIANGLE* tri=(M3D_TRIANGLE*)triangle1;
288     MG_TETRA* tet=tri->get_voisin(0);
289     if (tet->get_type_entite()==IDM3D_TETRA)
290     {
291     M3D_TETRA* mtet=(M3D_TETRA*)tet;
292     vol1=mtet->get_volume();
293     }
294     }
295     if (triangle2!=NULL)
296     if (triangle2->get_type_entite()==IDM3D_TRIANGLE)
297     {
298     M3D_TRIANGLE* tri=(M3D_TRIANGLE*)triangle2;
299     MG_TETRA* tet=tri->get_voisin(0);
300     if (tet->get_type_entite()==IDM3D_TETRA)
301     {
302     M3D_TETRA* mtet=(M3D_TETRA*)tet;
303     vol2=mtet->get_volume();
304     }
305     }
306     if (triangle3!=NULL)
307     if (triangle3->get_type_entite()==IDM3D_TRIANGLE)
308     {
309     M3D_TRIANGLE* tri=(M3D_TRIANGLE*)triangle3;
310     MG_TETRA* tet=tri->get_voisin(0);
311     if (tet->get_type_entite()==IDM3D_TETRA)
312     {
313     M3D_TETRA* mtet=(M3D_TETRA*)tet;
314     vol3=mtet->get_volume();
315     }
316     }
317     if (triangle4!=NULL)
318     if (triangle4->get_type_entite()==IDM3D_TRIANGLE)
319     {
320     M3D_TRIANGLE* tri=(M3D_TRIANGLE*)triangle4;
321     MG_TETRA* tet=tri->get_voisin(0);
322     if (tet->get_type_entite()==IDM3D_TETRA)
323     {
324     M3D_TETRA* mtet=(M3D_TETRA*)tet;
325     vol4=mtet->get_volume();
326     }
327     }
328     double vol=calcule_volume(noeud1,noeud2,noeud3,noeud4);
329     // calcul du vol du tetra ideal
330     OT_VECTEUR_3D ab(noeud2->get_x()-noeud1->get_x(),noeud2->get_y()-noeud1->get_y(),noeud2->get_z()-noeud1->get_z());
331     OT_VECTEUR_3D ac(noeud3->get_x()-noeud1->get_x(),noeud3->get_y()-noeud1->get_y(),noeud3->get_z()-noeud1->get_z());
332     OT_VECTEUR_3D da(noeud1->get_x()-noeud4->get_x(),noeud1->get_y()-noeud4->get_y(),noeud1->get_z()-noeud4->get_z());
333     OT_VECTEUR_3D vec=ab&ac;
334     vec.norme();
335     double lg_tet_ideal=fabs(vec*da*ft->get_ifail()*sqrt(3./2.));
336     double vol_ideal=sqrt(2.)/12.*lg_tet_ideal*lg_tet_ideal*lg_tet_ideal;
337     if (vol<0.1*vol_ideal)
338     {
339     if ((vol1>0.) && (vol1<0.1*vol_ideal)) {clock_t cpu2=clock(); //debug
340     cpu[11]=cpu[11]+(cpu2-cpu1)/CLK_TCK;
341     nbfois[11]++;
342     return NULL;}
343     if ((vol2>0.) && (vol2<0.1*vol_ideal)) {clock_t cpu2=clock(); //debug
344     cpu[11]=cpu[11]+(cpu2-cpu1)/CLK_TCK;
345     nbfois[11]++;
346     return NULL;}
347     if ((vol3>0.) && (vol3<0.1*vol_ideal)) {clock_t cpu2=clock(); //debug
348     cpu[11]=cpu[11]+(cpu2-cpu1)/CLK_TCK;
349     nbfois[11]++;
350     return NULL;}
351     if ((vol4>0.) && (vol4<0.1*vol_ideal)) {clock_t cpu2=clock(); //debug
352     cpu[11]=cpu[11]+(cpu2-cpu1)/CLK_TCK;
353     nbfois[11]++;
354     return NULL;}
355     }
356     }
357    
358    
359    
360    
361     // insertion du tetra
362 francois 54 if (triangle1==NULL) triangle1=insere_triangle(mgvol,noeud1,noeud3,noeud2,MAILLEUR_AUTO);
363     if (triangle2==NULL) triangle2=insere_triangle(mgvol,noeud1,noeud2,noeud4,MAILLEUR_AUTO);
364     if (triangle3==NULL) triangle3=insere_triangle(mgvol,noeud2,noeud3,noeud4,MAILLEUR_AUTO);
365     if (triangle4==NULL) triangle4=insere_triangle(mgvol,noeud1,noeud4,noeud3,MAILLEUR_AUTO);
366     MG_TETRA* tetra=cree_tetra(mgvol,noeud1,noeud2,noeud3,noeud4,triangle1,triangle2,triangle3,triangle4,MAILLEUR_AUTO);
367 5 mg_maillage->ajouter_mg_tetra(tetra);
368     if (triangle1->get_type_entite()==IDM3D_TRIANGLE)
369     {
370     M3D_TRIANGLE* tri=(M3D_TRIANGLE*)triangle1;
371     tri->ajouter_voisin(tetra);
372     }
373     if (triangle2->get_type_entite()==IDM3D_TRIANGLE)
374     {
375     M3D_TRIANGLE* tri=(M3D_TRIANGLE*)triangle2;
376     tri->ajouter_voisin(tetra);
377     }
378     if (triangle3->get_type_entite()==IDM3D_TRIANGLE)
379     {
380     M3D_TRIANGLE* tri=(M3D_TRIANGLE*)triangle3;
381     tri->ajouter_voisin(tetra);
382     }
383     if (triangle4->get_type_entite()==IDM3D_TRIANGLE)
384     {
385     M3D_TRIANGLE* tri=(M3D_TRIANGLE*)triangle4;
386     tri->ajouter_voisin(tetra);
387     }
388     clock_t cpu2=clock(); //debug
389     cpu[11]=cpu[11]+(cpu2-cpu1)/CLK_TCK;
390     nbfois[11]++;
391    
392     return tetra;
393    
394     }
395    
396    
397    
398 francois 54 class MG_TRIANGLE* MAILLEUR3D::insere_triangle(class MG_VOLUME* mgvol,MG_NOEUD* noeud1,MG_NOEUD* noeud2,MG_NOEUD* noeud3,int origine)
399 5 {
400     clock_t cpu1=clock(); //debug
401     // existence du triangle
402     MG_TRIANGLE* triangle=mg_maillage->get_mg_triangle(noeud1->get_id(),noeud3->get_id(),noeud2->get_id());
403     if (triangle!=NULL) return NULL;
404     // insertion du triangle
405     MG_SEGMENT* segment1=mg_maillage->get_mg_segment(noeud1->get_id(),noeud2->get_id());
406     MG_SEGMENT* segment2=mg_maillage->get_mg_segment(noeud1->get_id(),noeud3->get_id());
407     MG_SEGMENT* segment3=mg_maillage->get_mg_segment(noeud2->get_id(),noeud3->get_id());
408 francois 54 if (segment1==NULL) segment1=cree_segment(mgvol,noeud1,noeud2,origine);
409     if (segment2==NULL) segment2=cree_segment(mgvol,noeud1,noeud3,origine);
410     if (segment3==NULL) segment3=cree_segment(mgvol,noeud2,noeud3,origine);
411     MG_TRIANGLE* tri=cree_triangle(mgvol,noeud1,noeud2,noeud3,segment1,segment2,segment3,origine);
412 5 clock_t cpu2=clock(); //debug
413     cpu[6]=cpu[6]+(cpu2-cpu1)/CLK_TCK;
414     nbfois[6]++;
415     return tri;
416     }
417    
418    
419    
420 francois 54 class MG_SEGMENT* MAILLEUR3D::insere_segment(class MG_VOLUME* mgvol,MG_NOEUD* noeud1,MG_NOEUD* noeud2,int origine)
421 5 {
422     clock_t cpu1=clock(); //debug
423     //existence du segment
424     MG_SEGMENT* segment=mg_maillage->get_mg_segment(noeud1->get_id(),noeud2->get_id());
425     if (segment!=NULL) return NULL;
426     // insertion du segment
427 francois 54 MG_SEGMENT* seg=cree_segment(mgvol,noeud1,noeud2,origine);
428 5 clock_t cpu2=clock(); //debug
429     cpu[7]=cpu[7]+(cpu2-cpu1)/CLK_TCK;
430     nbfois[7]++;
431     return seg;
432     }
433    
434    
435    
436    
437 francois 54 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,int origine)
438 5 {
439 francois 54 M3D_TETRA* tetra=new M3D_TETRA(mgvol,noeud1,noeud2,noeud3,noeud4,triangle1,triangle2,triangle3,triangle4,origine);
440 5 mg_maillage->ajouter_mg_tetra(tetra);
441     return tetra;
442     }
443    
444    
445 francois 54 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,int origine)
446 5 {
447 francois 54 M3D_TRIANGLE* tri=new M3D_TRIANGLE(mgvol,noeud1,noeud2,noeud3,segment1,segment2,segment3,origine);
448 5 mg_maillage->ajouter_mg_triangle(tri);
449     return tri;
450     }
451    
452    
453    
454 francois 54 MG_SEGMENT* MAILLEUR3D::cree_segment(class MG_VOLUME* mgvol,MG_NOEUD* noeud1,MG_NOEUD* noeud2,int origine)
455 5 {
456 francois 54 MG_SEGMENT* seg=mg_maillage->ajouter_mg_segment(mgvol,noeud1,noeud2,origine);
457 5 return seg;
458     }
459    
460    
461    
462     double MAILLEUR3D::calcule_volume(MG_NOEUD* noeud1,MG_NOEUD* noeud2,MG_NOEUD* noeud3,MG_NOEUD* noeud4)
463     {
464     OT_VECTEUR_3D ab(noeud2->get_x()-noeud1->get_x(),noeud2->get_y()-noeud1->get_y(),noeud2->get_z()-noeud1->get_z());
465     OT_VECTEUR_3D ac(noeud3->get_x()-noeud1->get_x(),noeud3->get_y()-noeud1->get_y(),noeud3->get_z()-noeud1->get_z());
466     OT_VECTEUR_3D da(noeud1->get_x()-noeud4->get_x(),noeud1->get_y()-noeud4->get_y(),noeud1->get_z()-noeud4->get_z());
467    
468     OT_VECTEUR_3D pvec=ab&ac;
469     double volume=pvec*da;
470     volume=fabs(volume)/3.;
471     return volume;
472     }
473    
474    
475    
476     double MAILLEUR3D::calcule_longueur_caracteristique(MG_NOEUD* noeud1,MG_NOEUD* noeud2,MG_NOEUD* noeud3,MG_NOEUD* noeud4)
477     {
478    
479     OT_VECTEUR_3D ab(noeud2->get_x()-noeud1->get_x(),noeud2->get_y()-noeud1->get_y(),noeud2->get_z()-noeud1->get_z());
480     double lg=ab.get_longueur();
481     int nb=1;
482     if (noeud3!=NULL)
483     {
484     OT_VECTEUR_3D ac(noeud3->get_x()-noeud1->get_x(),noeud3->get_y()-noeud1->get_y(),noeud3->get_z()-noeud1->get_z());
485     OT_VECTEUR_3D bc(noeud3->get_x()-noeud2->get_x(),noeud3->get_y()-noeud2->get_y(),noeud3->get_z()-noeud2->get_z());
486     lg=lg+ac.get_longueur()+bc.get_longueur();
487     nb=3;
488     }
489     if (noeud4!=NULL)
490     {
491     OT_VECTEUR_3D da(noeud1->get_x()-noeud4->get_x(),noeud1->get_y()-noeud4->get_y(),noeud1->get_z()-noeud4->get_z());
492     OT_VECTEUR_3D bd(noeud4->get_x()-noeud2->get_x(),noeud4->get_y()-noeud2->get_y(),noeud4->get_z()-noeud2->get_z());
493     OT_VECTEUR_3D cd(noeud4->get_x()-noeud3->get_x(),noeud4->get_y()-noeud3->get_y(),noeud4->get_z()-noeud3->get_z());
494     lg=lg+da.get_longueur()+bd.get_longueur()+cd.get_longueur();
495     nb=6;
496     }
497     lg=lg*1./nb;
498     return lg;
499    
500     }
501    
502    
503    
504    
505 francois 35 int MAILLEUR3D::genere_noeud(MG_VOLUME* mgvol,MG_FRONT_3D* ft,TPL_MAP_ENTITE<MG_NOEUD*> &liste_noeud,std::vector<CAS_FRONT*> &liste_cas)
506 5 {
507     clock_t cpu1=clock(); //debug
508     MG_NOEUD* noeud1=ft->get_noeud1();
509     MG_NOEUD* noeud2=ft->get_noeud2();
510     MG_NOEUD* noeud3=ft->get_noeud3();
511     OT_VECTEUR_3D n1n2(noeud2->get_x()-noeud1->get_x(),noeud2->get_y()-noeud1->get_y(),noeud2->get_z()-noeud1->get_z());
512     OT_VECTEUR_3D n1n3(noeud3->get_x()-noeud1->get_x(),noeud3->get_y()-noeud1->get_y(),noeud3->get_z()-noeud1->get_z());
513     OT_VECTEUR_3D normal=n1n2&n1n3;
514     normal.norme();
515     double xyz1[3]={noeud1->get_x(),noeud1->get_y(),noeud1->get_z()};
516     double xyz2[3]={noeud2->get_x(),noeud2->get_y(),noeud2->get_z()};
517     double xyz3[3]={noeud3->get_x(),noeud3->get_y(),noeud3->get_z()};
518 francois 35 double hauteur=0.81649658*ft->get_triangle()->get_longueur();
519     double xyzg[3]={0.333333333333*(xyz1[0]+xyz2[0]+xyz3[0]),0.333333333333*(xyz1[1]+xyz2[1]+xyz3[1]),0.333333333333*(xyz1[2]+xyz2[2]+xyz3[2])};
520 francois 61 double xyz4[3];
521     do
522     {
523     xyz4[0]=xyzg[0]+hauteur*normal.get_x();
524     xyz4[1]=xyzg[1]+hauteur*normal.get_y();
525     xyz4[2]=xyzg[2]+hauteur*normal.get_z();
526     hauteur=hauteur*0.9;
527     }
528     while (metrique->valide_parametre(xyz4)==0);
529 francois 35 double l1=calcul_distance_metrique(xyz1,xyz2);
530     double l2=calcul_distance_metrique(xyz1,xyz3);
531     double l3=calcul_distance_metrique(xyz2,xyz3);
532     double l=0.333333333333333*(l1+l2+l3);
533     double longueur_desiree=0.81649658*(priorite_metrique+(1-priorite_metrique)*l);
534     longueur_desiree=longueur_desiree/ft->get_ifail();
535     ajuste_distance_metrique(xyzg,xyz4,longueur_desiree);
536     OT_VECTEUR_3D hauteurajustee(xyzg[0]-xyz4[0],xyzg[1]-xyz4[1],xyzg[2]-xyz4[2]);
537     hauteur=hauteurajustee.get_longueur();
538 5
539     // recherche des voisins
540     TPL_MAP_ENTITE<MG_FRONT_3D*> liste_front;
541     TPL_MAP_ENTITE<MG_FRONT_3D*> liste_front_proche;
542 francois 35 octree_de_front->rechercher(xyz4[0],xyz4[1],xyz4[2],hauteur*1.837117307,liste_front);
543 5 int nb_liste=liste_front.get_nb();
544 francois 35 double disref=(priorite_metrique+(1-priorite_metrique)*l)/ft->get_ifail();
545 5 MG_NOEUD* noeudref=NULL;
546     for (int i=0;i<nb_liste;i++)
547     {
548     MG_FRONT_3D* ft2=liste_front.get(i);
549     for (int j=0;j<3;j++)
550     {
551     MG_NOEUD* noeud_tmp;
552     if (j==0) noeud_tmp=ft2->get_noeud1();
553     if (j==1) noeud_tmp=ft2->get_noeud2();
554     if (j==2) noeud_tmp=ft2->get_noeud3();
555 francois 35 double *xyz_tmp=noeud_tmp->get_coord();
556     double dis=calcul_distance_metrique(xyz_tmp,xyz4,1);
557 5 if (dis<disref)
558     {
559     disref=dis;
560     noeudref=noeud_tmp;
561     }
562     }
563     if (ft2!=ft)
564     {
565 francois 35 double dismin=evaluer_distance_noeud_triangle(xyz4[0],xyz4[1],xyz4[2],ft2->get_triangle());
566     if (dismin<hauteur) liste_front_proche.ajouter(ft2);
567 5 }
568     }
569     // il y a un noeud proche
570     if (noeudref!=NULL)
571 francois 35 {
572     double dis2=calcul_distance_metrique(noeudref->get_coord(),xyzg,2);
573     if (dis2<longueur_desiree*1.5)
574 5 if (noeudref!=noeud1)
575     if (noeudref!=noeud2)
576     if (noeudref!=noeud3)
577     {
578     liste_noeud.ajouter(noeudref);
579     clock_t cpu2=clock(); //debug
580     cpu[8]=cpu[8]+(cpu2-cpu1)/CLK_TCK;
581     nbfois[8]++;
582     return EXISTE;
583     }
584 francois 35 }
585 5 // il n y a pas de noeud proche mais il y a une face proche
586     int nb_proche=liste_front_proche.get_nb();
587 francois 35 int nb_cas=liste_cas.size();
588 5 if (nb_proche>0)
589     {
590     for (int i=0;i<nb_proche;i++)
591     {
592     MG_FRONT_3D* fttmp=liste_front_proche.get(i);
593 francois 35 int ok=0;
594     for (int j=0;j<nb_cas;j++)
595     {
596     MG_NOEUD* casnoeud=liste_cas[j]->mgnoeud;
597     if ((fttmp->get_noeud1()==casnoeud) || (fttmp->get_noeud2()==casnoeud) || (fttmp->get_noeud3()==casnoeud) )
598     {
599     ok=1;
600     break;
601     }
602     }
603     if (ok==1) continue;
604 5 MG_NOEUD* noeudtmp=fttmp->get_noeud1();
605     if (noeudtmp!=noeud1)
606     if (noeudtmp!=noeud2)
607     if (noeudtmp!=noeud3) liste_noeud.ajouter(noeudtmp);
608     noeudtmp=fttmp->get_noeud2();
609     if (noeudtmp!=noeud1)
610     if (noeudtmp!=noeud2)
611     if (noeudtmp!=noeud3) liste_noeud.ajouter(noeudtmp);
612     noeudtmp=fttmp->get_noeud3();
613     if (noeudtmp!=noeud1)
614     if (noeudtmp!=noeud2)
615     if (noeudtmp!=noeud3) liste_noeud.ajouter(noeudtmp);
616    
617     }
618     clock_t cpu2=clock(); //debug
619     cpu[8]=cpu[8]+(cpu2-cpu1)/CLK_TCK;
620     nbfois[8]++;
621 francois 35 if (liste_noeud.get_nb()>0) return MULTIEXISTE;
622 5 }
623     // il faut creer un noeud c est parfait c est qu on veut.
624 francois 54 M3D_NOEUD* nouveau_noeud=new M3D_NOEUD(mgvol,xyz4[0],xyz4[1],xyz4[2],MAILLEUR_AUTO);
625 5 mg_maillage->ajouter_mg_noeud(nouveau_noeud);
626     liste_noeud.ajouter(nouveau_noeud);
627     clock_t cpu2=clock(); //debug
628     cpu[8]=cpu[8]+(cpu2-cpu1)/CLK_TCK;
629     nbfois[8]++;
630     return CREATION;
631     }
632    
633    
634    
635     double MAILLEUR3D::evaluer_distance_noeud_triangle(double x,double y,double z,MG_TRIANGLE* tri)
636     {
637     double eps=0.0001;
638     MG_NOEUD* noeud1=tri->get_noeud1();
639     MG_NOEUD* noeud2=tri->get_noeud2();
640     MG_NOEUD* noeud3=tri->get_noeud3();
641     OT_VECTEUR_3D n1n2(noeud2->get_x()-noeud1->get_x(),noeud2->get_y()-noeud1->get_y(),noeud2->get_z()-noeud1->get_z());
642     OT_VECTEUR_3D n1n3(noeud3->get_x()-noeud1->get_x(),noeud3->get_y()-noeud1->get_y(),noeud3->get_z()-noeud1->get_z());
643     OT_VECTEUR_3D n1n4(x-noeud1->get_x(),y-noeud1->get_y(),z-noeud1->get_z());
644     OT_VECTEUR_3D n1n=n1n2&n1n3;
645     OT_MATRICE_3D matbase(n1n2,n1n3,n1n);
646     OT_MATRICE_3D matxsi(n1n4,n1n3,n1n);
647     OT_MATRICE_3D mateta(n1n2,n1n4,n1n);
648     double xsi=matxsi.get_determinant()/matbase.get_determinant();
649     double eta=mateta.get_determinant()/matbase.get_determinant();
650     double dis;
651     if ((eta>-eps) && (xsi>-eps) && (eta+xsi<1.0+eps))
652     {
653     double xproj=(1-xsi-eta)*noeud1->get_x()+xsi*noeud2->get_x()+eta*noeud3->get_x();
654     double yproj=(1-xsi-eta)*noeud1->get_y()+xsi*noeud2->get_y()+eta*noeud3->get_y();
655     double zproj=(1-xsi-eta)*noeud1->get_z()+xsi*noeud2->get_z()+eta*noeud3->get_z();
656     OT_VECTEUR_3D vec(x-xproj,y-yproj,z-zproj);
657     dis=vec.get_longueur();
658     }
659     else
660     {
661     OT_VECTEUR_3D nn1(x-noeud1->get_x(),y-noeud1->get_y(),z-noeud1->get_z());
662     OT_VECTEUR_3D nn2(x-noeud2->get_x(),y-noeud2->get_y(),z-noeud2->get_z());
663     OT_VECTEUR_3D nn3(x-noeud3->get_x(),y-noeud3->get_y(),z-noeud3->get_z());
664     dis=0.3333333333333*(nn1.get_longueur()+nn2.get_longueur()+nn3.get_longueur());
665     }
666     return dis;
667     }
668    
669    
670     void MAILLEUR3D::detruit_tetra(MG_TRIANGLE* tri)
671     {
672     clock_t cpu1=clock(); //debug
673    
674     M3D_TRIANGLE* mtri;
675     //if (tri->get_type_entite()==IDMG_TRIANGLE) return;
676     //else mtri=(M3D_TRIANGLE*)tri;
677     mtri=(M3D_TRIANGLE*)tri;
678    
679     M3D_TETRA* tet=(M3D_TETRA*)mtri->get_voisin(0);
680 francois 54 M3D_TRIANGLE* mtriangle1=(M3D_TRIANGLE*)tet->get_triangle1();
681     M3D_TRIANGLE* mtriangle2=(M3D_TRIANGLE*)tet->get_triangle2();
682     M3D_TRIANGLE* mtriangle3=(M3D_TRIANGLE*)tet->get_triangle3();
683     M3D_TRIANGLE* mtriangle4=(M3D_TRIANGLE*)tet->get_triangle4();
684    
685    
686 5 tet->etat_detruit();
687     TPL_MAP_ENTITE<MG_FRONT_3D*> liste_de_front;
688     TPL_MAP_ENTITE<MG_FRONT_3D*> liste_de_front2;
689     MG_FRONT_3D* fttet1=NULL;
690     MG_FRONT_3D* fttet2=NULL;
691     MG_FRONT_3D* fttet3=NULL;
692     MG_FRONT_3D* fttet4=NULL;
693    
694 francois 54 if (mtriangle1->get_etat_front()==FRONT_ACTIF)
695 5 {
696 francois 54 MG_FRONT_3D* fttmp=mtriangle1->get_mgfront();
697 5 fttet1=fttmp;
698     liste_de_front.ajouter(fttmp->get_front_voisin(0));
699     liste_de_front.ajouter(fttmp->get_front_voisin(1));
700     liste_de_front.ajouter(fttmp->get_front_voisin(2));
701     //supprimer_front_en_reculant(fttmp);
702 francois 54 mtriangle1->supprimer_voisin(tet);
703 5 }
704     else
705     {
706 francois 54 MG_FRONT_3D* ftn1=ajouter_front_courant(PRIORITAIRE,tet->get_noeud1(),tet->get_noeud2(),tet->get_noeud3(),mtriangle1);
707 5 ftn1->incremente_ifail();
708 francois 54 mtriangle1->supprimer_voisin(tet);
709 5 liste_de_front2.ajouter(ftn1);
710     }
711 francois 54 if (mtriangle2->get_etat_front()==FRONT_ACTIF)
712 5 {
713 francois 54 MG_FRONT_3D* fttmp=mtriangle2->get_mgfront();
714 5 fttet2=fttmp;
715     liste_de_front.ajouter(fttmp->get_front_voisin(0));
716     liste_de_front.ajouter(fttmp->get_front_voisin(1));
717     liste_de_front.ajouter(fttmp->get_front_voisin(2));
718     //supprimer_front_en_reculant(fttmp);
719 francois 54 mtriangle2->supprimer_voisin(tet);
720 5 }
721     else
722     {
723 francois 54 MG_FRONT_3D* ftn1=ajouter_front_courant(PRIORITAIRE,tet->get_noeud1(),tet->get_noeud4(),tet->get_noeud2(),mtriangle2);
724 5 ftn1->incremente_ifail();
725 francois 54 mtriangle2->supprimer_voisin(tet);
726 5 liste_de_front2.ajouter(ftn1);
727     }
728 francois 54 if (mtriangle3->get_etat_front()==FRONT_ACTIF)
729 5 {
730 francois 54 MG_FRONT_3D* fttmp=mtriangle3->get_mgfront();
731 5 fttet3=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 francois 54 mtriangle3->supprimer_voisin(tet);
737 5 }
738     else
739     {
740 francois 54 MG_FRONT_3D* ftn1=ajouter_front_courant(PRIORITAIRE,tet->get_noeud2(),tet->get_noeud4(),tet->get_noeud3(),mtriangle3);
741 5 ftn1->incremente_ifail();
742 francois 54 mtriangle3->supprimer_voisin(tet);
743 5 liste_de_front2.ajouter(ftn1);
744     }
745 francois 54 if (mtriangle4->get_etat_front()==FRONT_ACTIF)
746 5 {
747 francois 54 MG_FRONT_3D* fttmp=mtriangle4->get_mgfront();
748 5 fttet4=fttmp;
749     liste_de_front.ajouter(fttmp->get_front_voisin(0));
750     liste_de_front.ajouter(fttmp->get_front_voisin(1));
751     liste_de_front.ajouter(fttmp->get_front_voisin(2));
752     //supprimer_front_en_reculant(fttmp);
753 francois 54 mtriangle4->supprimer_voisin(tet);
754 5 }
755     else
756     {
757 francois 54 MG_FRONT_3D* ftn1=ajouter_front_courant(PRIORITAIRE,tet->get_noeud1(),tet->get_noeud3(),tet->get_noeud4(),mtriangle4);
758 5 ftn1->incremente_ifail();
759 francois 54 mtriangle4->supprimer_voisin(tet);
760 5 liste_de_front2.ajouter(ftn1);
761     }
762    
763     if (fttet1!=NULL)
764     {
765     liste_de_front.supprimer(fttet1);
766     supprimer_front_en_reculant(fttet1);
767     }
768     if (fttet2!=NULL)
769     {
770     liste_de_front.supprimer(fttet2);
771     supprimer_front_en_reculant(fttet2);
772     }
773     if (fttet3!=NULL)
774     {
775     liste_de_front.supprimer(fttet3);
776     supprimer_front_en_reculant(fttet3);
777     }
778     if (fttet4!=NULL)
779     {
780     liste_de_front.supprimer(fttet4);
781     supprimer_front_en_reculant(fttet4);
782     }
783     int nb_liste1=liste_de_front.get_nb();
784     for (int i=0;i<nb_liste1;i++)
785     if (liste_de_front.get(i)!=NULL)
786     mise_a_jour_voisin_front(liste_de_front.get(i));
787     int nb_liste2=liste_de_front2.get_nb();
788     for (int i=0;i<nb_liste2;i++)
789     {
790     mise_a_jour_voisin_front(liste_de_front2.get(i));
791     mise_a_jour_voisin_front(liste_de_front2.get(i)->get_front_voisin(0));
792     mise_a_jour_voisin_front(liste_de_front2.get(i)->get_front_voisin(1));
793     mise_a_jour_voisin_front(liste_de_front2.get(i)->get_front_voisin(2));
794     }
795    
796     for (int k=0;k<4;k++)
797     {
798     MG_NOEUD* noeud;
799     if (k==0) noeud=tet->get_noeud1();
800     if (k==1) noeud=tet->get_noeud2();
801     if (k==2) noeud=tet->get_noeud3();
802     if (k==3) noeud=tet->get_noeud4();
803     int nb_tet=noeud->get_lien_tetra()->get_nb();
804    
805     if (noeud->get_type_entite()==IDM3D_NOEUD)
806     {
807     int test=false;
808     M3D_NOEUD* mnoeud=(M3D_NOEUD*)noeud;
809     for (int i=0;i<nb_tet;i++)
810     {
811     MG_TETRA* tet=noeud->get_lien_tetra()->get(i);
812     if (tet->get_type_entite()==IDM3D_TETRA)
813     {
814     M3D_TETRA* mtet=(M3D_TETRA*)tet;
815     if (mtet->get_etat()==ACTIF) test=true;
816     }
817     else test=true;
818     }
819     if (test==false) mnoeud->etat_detruit();
820     }
821     }
822     clock_t cpu2=clock();
823     cpu[9]=cpu[9]+(cpu2-cpu1)/CLK_TCK;
824     nbfois[9]++;
825    
826     }
827    
828 francois 35 double MAILLEUR3D::calcul_distance_metrique(MG_SEGMENT* seg,int pas)
829     {
830     double *xyz1=seg->get_noeud1()->get_coord();
831     double *xyz2=seg->get_noeud2()->get_coord();
832     return calcul_distance_metrique(xyz1,xyz2,pas);
833     }
834 5
835 francois 35 double MAILLEUR3D::calcul_distance_metrique(double *xyz1,double* xyz2,int pas)
836     {
837     double longueur=0.;
838     double dxyz[3];
839     dxyz[0]=xyz2[0]-xyz1[0];
840     dxyz[1]=xyz2[1]-xyz1[1];
841     dxyz[2]=xyz2[2]-xyz1[2];
842     for (int i=0;i<pas;i++)
843     {
844     double ti=1.0*i/pas;
845     double tii=1.0*(i+1)/pas;
846     double tgauss1=0.7886751345*tii+0.2113248654*ti;
847     double xyz[3];
848     xyz[0]=xyz1[0]+tgauss1*(xyz2[0]-xyz1[0]);
849     xyz[1]=xyz1[1]+tgauss1*(xyz2[1]-xyz1[1]);
850     xyz[2]=xyz1[2]+tgauss1*(xyz2[2]-xyz1[2]);
851     double tenseur[9];
852     metrique->evaluer(xyz,tenseur);
853     double val= sqrt(dxyz[0]*dxyz[0]*tenseur[0]+dxyz[0]*dxyz[1]*tenseur[3]+dxyz[0]*dxyz[2]*tenseur[6]+dxyz[1]*dxyz[0]*tenseur[1]+dxyz[1]*dxyz[1]*tenseur[4]+dxyz[1]*dxyz[2]*tenseur[7]+dxyz[2]*dxyz[0]*tenseur[2]+dxyz[2]*dxyz[1]*tenseur[5]+dxyz[2]*dxyz[2]*tenseur[8]);
854     longueur=longueur+0.5*(tii-ti)*val;
855     double tgauss2=0.7886751345*ti+0.2113248654*tii;
856     xyz[0]=xyz1[0]+tgauss2*(xyz2[0]-xyz1[0]);
857     xyz[1]=xyz1[1]+tgauss2*(xyz2[1]-xyz1[1]);
858     xyz[2]=xyz1[2]+tgauss2*(xyz2[2]-xyz1[2]);
859     metrique->evaluer(xyz,tenseur);
860     val= sqrt(dxyz[0]*dxyz[0]*tenseur[0]+dxyz[0]*dxyz[1]*tenseur[3]+dxyz[0]*dxyz[2]*tenseur[6]+dxyz[1]*dxyz[0]*tenseur[1]+dxyz[1]*dxyz[1]*tenseur[4]+dxyz[1]*dxyz[2]*tenseur[7]+dxyz[2]*dxyz[0]*tenseur[2]+dxyz[2]*dxyz[1]*tenseur[5]+dxyz[2]*dxyz[2]*tenseur[8]);
861     longueur=longueur+0.5*(tii-ti)*val;
862     }
863     return longueur;
864     }
865 5
866 francois 35 void MAILLEUR3D::ajuste_distance_metrique(double *xyz1,double *xyz2,double longueur_desiree,int pas)
867     {
868     double longueur=0,ancienne_longueur;
869     double ti,tii;
870     double dxyz[3];
871     dxyz[0]=xyz2[0]-xyz1[0];
872     dxyz[1]=xyz2[1]-xyz1[1];
873     dxyz[2]=xyz2[2]-xyz1[2];
874     int i=0;
875     do
876     {
877     ancienne_longueur=longueur;
878     ti=1.0*i/pas;
879     tii=1.0*(i+1)/pas;
880     double tgauss1=0.7886751345*tii+0.2113248654*ti;
881     double xyz[3];
882     xyz[0]=xyz1[0]+tgauss1*(xyz2[0]-xyz1[0]);
883     xyz[1]=xyz1[1]+tgauss1*(xyz2[1]-xyz1[1]);
884     xyz[2]=xyz1[2]+tgauss1*(xyz2[2]-xyz1[2]);
885 francois 61 if (i+1>pas)
886     {
887     int val=metrique->valide_parametre(xyz);
888     if (!val) return;
889     }
890 francois 35 double tenseur[9];
891     metrique->evaluer(xyz,tenseur);
892     double val= sqrt(dxyz[0]*dxyz[0]*tenseur[0]+dxyz[0]*dxyz[1]*tenseur[3]+dxyz[0]*dxyz[2]*tenseur[6]+dxyz[1]*dxyz[0]*tenseur[1]+dxyz[1]*dxyz[1]*tenseur[4]+dxyz[1]*dxyz[2]*tenseur[7]+dxyz[2]*dxyz[0]*tenseur[2]+dxyz[2]*dxyz[1]*tenseur[5]+dxyz[2]*dxyz[2]*tenseur[8]);
893     longueur=longueur+0.5*(tii-ti)*val;
894     double tgauss2=0.7886751345*ti+0.2113248654*tii;
895     xyz[0]=xyz1[0]+tgauss2*(xyz2[0]-xyz1[0]);
896     xyz[1]=xyz1[1]+tgauss2*(xyz2[1]-xyz1[1]);
897     xyz[2]=xyz1[2]+tgauss2*(xyz2[2]-xyz1[2]);
898 francois 61 if (i+1>pas)
899     {
900     int val=metrique->valide_parametre(xyz);
901     if (!val) return;
902     }
903 francois 35 metrique->evaluer(xyz,tenseur);
904     val= sqrt(dxyz[0]*dxyz[0]*tenseur[0]+dxyz[0]*dxyz[1]*tenseur[3]+dxyz[0]*dxyz[2]*tenseur[6]+dxyz[1]*dxyz[0]*tenseur[1]+dxyz[1]*dxyz[1]*tenseur[4]+dxyz[1]*dxyz[2]*tenseur[7]+dxyz[2]*dxyz[0]*tenseur[2]+dxyz[2]*dxyz[1]*tenseur[5]+dxyz[2]*dxyz[2]*tenseur[8]);
905     longueur=longueur+0.5*(tii-ti)*val;
906     i++;
907     }
908     while (longueur<longueur_desiree);
909     double t=ti+(tii-ti)*(longueur_desiree-ancienne_longueur)/(longueur-ancienne_longueur);
910     xyz2[0]=xyz1[0]+t*(xyz2[0]-xyz1[0]);
911     xyz2[1]=xyz1[1]+t*(xyz2[1]-xyz1[1]);
912     xyz2[2]=xyz1[2]+t*(xyz2[2]-xyz1[2]);
913     }
914    
915 5 void MAILLEUR3D::detruit_element_inactif(void)
916     {
917     clock_t cpu1=clock();
918 francois 146 TPL_MAP_ENTITE<MG_TETRA*> lsttet;
919     TPL_MAP_ENTITE<MG_NOEUD*> lstnoeud;
920     LISTE_MG_TETRA::iterator it;
921    
922     for (MG_TETRA* tet=mg_maillage->get_premier_tetra(it);tet!=NULL;tet=mg_maillage->get_suivant_tetra(it))
923 5 {
924     M3D_TETRA* mtet=(M3D_TETRA*)tet;
925     if (mtet->get_etat()==DETRUIT)
926 francois 146 lsttet.ajouter(tet);
927     }
928     LISTE_MG_NOEUD::iterator it2;
929     for (MG_NOEUD* noeud=mg_maillage->get_premier_noeud(it2);noeud!=NULL;noeud=mg_maillage->get_suivant_noeud(it2))
930 5 {
931     if (noeud->get_type_entite()==IDM3D_NOEUD)
932     {
933     M3D_NOEUD* mnoeud=(M3D_NOEUD*)noeud;
934     if (mnoeud->get_etat()==DETRUIT)
935 francois 146 lstnoeud.ajouter(noeud);
936     }
937 5 }
938 francois 146 TPL_MAP_ENTITE<MG_TETRA*>::ITERATEUR it3;
939     for (MG_TETRA* tet=lsttet.get_premier(it3);tet!=NULL;tet=lsttet.get_suivant(it3))
940     mg_maillage->supprimer_mg_tetraid(tet->get_id());
941     TPL_MAP_ENTITE<MG_NOEUD*>::ITERATEUR it4;
942     for (MG_NOEUD* noeud=lstnoeud.get_premier(it4);noeud!=NULL;noeud=lstnoeud.get_suivant(it4))
943     mg_maillage->supprimer_mg_noeudid(noeud->get_id());
944 5
945     clock_t cpu2=clock();
946     cpu[10]=cpu[10]+(cpu2-cpu1)/CLK_TCK;
947     nbfois[10]++;
948     }
949    
950 francois 35
951    
952     double MAILLEUR3D::calcul_volume_tetra_metrique(MG_TETRA* tet)
953     {
954     double *xyz1=tet->get_noeud1()->get_coord();
955     double *xyz2=tet->get_noeud2()->get_coord();
956     double *xyz3=tet->get_noeud3()->get_coord();
957     double *xyz4=tet->get_noeud4()->get_coord();
958     /*double xyzm[3];
959     double xyzg[3];
960     xyzm[0]=(xyz1[0]+xyz2[0])/2.;
961     xyzm[1]=(xyz1[1]+xyz2[1])/2.;
962     xyzm[2]=(xyz1[2]+xyz2[2])/2.;
963     xyzg[0]=(xyz1[0]+xyz2[0]+xyz3[0])/3.;
964     xyzg[1]=(xyz1[1]+xyz2[1]+xyz3[1])/3.;
965     xyzg[2]=(xyz1[2]+xyz2[2]+xyz3[2])/3.;
966     double base=calcul_distance_metrique(xyz1,xyz2);
967     double hauteur=calcul_distance_metrique(xyzm,xyz3);
968     double hauteur2=calcul_distance_metrique(xyzg,xyz4);
969     double vol=base*hauteur*hauteur2*1.414213562;*/
970     double a=0.1381966012;
971     double b=0.5854101965;
972     double poids_gauss[4],xsi_gauss[4],eta_gauss[4],zeta_gauss[4];
973     poids_gauss[0]=0.041666666666666;xsi_gauss[0]=a;eta_gauss[0]=a;zeta_gauss[0]=a;
974     poids_gauss[1]=0.041666666666666;xsi_gauss[1]=a;eta_gauss[1]=a;zeta_gauss[1]=b;
975     poids_gauss[2]=0.041666666666666;xsi_gauss[2]=a;eta_gauss[2]=b;zeta_gauss[2]=a;
976     poids_gauss[3]=0.041666666666666;xsi_gauss[3]=b;eta_gauss[3]=a;zeta_gauss[3]=a;
977     OT_VECTEUR_3D vec1(xyz1,xyz2);
978     OT_VECTEUR_3D vec2(xyz1,xyz3);
979     OT_VECTEUR_3D vec3(xyz1,xyz4);
980     double facteur=(vec1&vec2)*vec3;
981     double vol=0.;
982     for (int j=0;j<4;j++)
983     {
984     double xyz[3];
985     xyz[0]=(1-xsi_gauss[j]-eta_gauss[j]-zeta_gauss[j])*xyz1[0]+xsi_gauss[j]*xyz2[0]+eta_gauss[j]*xyz3[0]+zeta_gauss[j]*xyz4[0];
986     xyz[1]=(1-xsi_gauss[j]-eta_gauss[j]-zeta_gauss[j])*xyz1[1]+xsi_gauss[j]*xyz2[1]+eta_gauss[j]*xyz3[1]+zeta_gauss[j]*xyz4[1];
987     xyz[2]=(1-xsi_gauss[j]-eta_gauss[j]-zeta_gauss[j])*xyz1[2]+xsi_gauss[j]*xyz2[2]+eta_gauss[j]*xyz3[2]+zeta_gauss[j]*xyz4[2];
988     double tenseur[9];
989     metrique->evaluer(xyz,tenseur);
990     double dens=1/sqrt(tenseur[0]);
991     dens=dens*dens*dens;
992     vol=vol+poids_gauss[j]*facteur*8.485281374/dens;
993     }
994    
995     return vol;
996     }
997    
998    
999     void MAILLEUR3D::analyse_maillage_obtenu(double &vol)
1000     {
1001     int nb=mg_maillage->get_nb_mg_tetra();
1002     vol=0.;
1003     for (int i=0;i<nb;i++)
1004     {
1005     MG_TETRA* tet=mg_maillage->get_mg_tetra(i);
1006     vol=vol+calcul_volume_tetra_metrique(tet);
1007     }
1008     }
1009    
1010 francois 54
1011     void MAILLEUR3D::insere_contrainte_tetra(MG_VOLUME* mgvol,TPL_LISTE_ENTITE<MG_TETRA*> *lsttet)
1012     {
1013 francois 155 map<long,long,less<long> > lstcorrespondancenoeud;
1014 francois 54 int nbtet=lsttet->get_nb();
1015     for (int i=0;i<nbtet;i++)
1016     {
1017     M3D_TETRA* tet=(M3D_TETRA*)lsttet->get(i);
1018     double *xyz1=tet->get_noeud1()->get_coord();
1019     double *xyz2=tet->get_noeud2()->get_coord();
1020     double *xyz3=tet->get_noeud3()->get_coord();
1021     double *xyz4=tet->get_noeud4()->get_coord();
1022     double xyz[3];
1023     xyz[0]=(xyz1[0]+xyz2[0]+xyz3[0]+xyz4[0])*0.25;
1024     xyz[1]=(xyz1[1]+xyz2[1]+xyz3[1]+xyz4[1])*0.25;
1025     xyz[2]=(xyz1[2]+xyz2[2]+xyz3[2]+xyz4[2])*0.25;
1026     TPL_MAP_ENTITE<MG_FRONT_3D*> liste_trouvee;
1027     OT_VECTEUR_3D vec1(xyz,xyz1);
1028     OT_VECTEUR_3D vec2(xyz,xyz2);
1029     OT_VECTEUR_3D vec3(xyz,xyz3);
1030     OT_VECTEUR_3D vec4(xyz,xyz4);
1031     double rayonrecherche=vec1.get_longueur();
1032     if (vec2.get_longueur()>rayonrecherche) rayonrecherche=vec2.get_longueur();
1033     if (vec3.get_longueur()>rayonrecherche) rayonrecherche=vec3.get_longueur();
1034     if (vec4.get_longueur()>rayonrecherche) rayonrecherche=vec4.get_longueur();
1035     rayonrecherche=rayonrecherche*1.1;
1036     octree_de_front->rechercher(xyz[0],xyz[1],xyz[2],rayonrecherche,liste_trouvee);
1037     MG_NOEUD* noeudref1=NULL;
1038 francois 155 if (lstcorrespondancenoeud.find(tet->get_noeud1()->get_nouveau_numero())!=lstcorrespondancenoeud.end())
1039     noeudref1=mg_maillage->get_mg_noeudid(lstcorrespondancenoeud[tet->get_noeud1()->get_nouveau_numero()]);
1040     MG_NOEUD* noeudref2=NULL;
1041     if (lstcorrespondancenoeud.find(tet->get_noeud2()->get_nouveau_numero())!=lstcorrespondancenoeud.end())
1042     noeudref2=mg_maillage->get_mg_noeudid(lstcorrespondancenoeud[tet->get_noeud2()->get_nouveau_numero()]);
1043     MG_NOEUD* noeudref3=NULL;
1044     if (lstcorrespondancenoeud.find(tet->get_noeud3()->get_nouveau_numero())!=lstcorrespondancenoeud.end())
1045     noeudref3=mg_maillage->get_mg_noeudid(lstcorrespondancenoeud[tet->get_noeud3()->get_nouveau_numero()]);
1046     MG_NOEUD* noeudref4=NULL;
1047     if (lstcorrespondancenoeud.find(tet->get_noeud4()->get_nouveau_numero())!=lstcorrespondancenoeud.end())
1048     noeudref4=mg_maillage->get_mg_noeudid(lstcorrespondancenoeud[tet->get_noeud4()->get_nouveau_numero()]);
1049     int nbtrouve=liste_trouvee.get_nb();
1050 francois 54 for (int j=0;j<nbtrouve;j++)
1051     {
1052 francois 155 if ((noeudref1!=NULL) && (noeudref2!=NULL) && (noeudref3!=NULL) && (noeudref4!=NULL) ) break;
1053     MG_FRONT_3D* ft=liste_trouvee.get(j);
1054 francois 54 TPL_MAP_ENTITE<MG_NOEUD*> lstnoeud;
1055     lstnoeud.ajouter(ft->get_noeud1());
1056     lstnoeud.ajouter(ft->get_noeud2());
1057     lstnoeud.ajouter(ft->get_noeud3());
1058 francois 155 if (noeudref1==NULL)
1059     {
1060     int nb=lstnoeud.get_nb();
1061     for (int k=0;k<nb;k++)
1062     {
1063     double *xyztmp=lstnoeud.get(k)->get_coord();
1064     OT_VECTEUR_3D vec(xyztmp,xyz1);
1065     double longueur=vec.get_longueur();
1066     if (longueur<1e-6*rayonrecherche)
1067     {
1068     noeudref1=lstnoeud.get(k);
1069     lstnoeud.supprimer(lstnoeud.get(k));
1070     break;
1071     }
1072     }
1073     }
1074     if (noeudref2==NULL)
1075     {
1076     int nb=lstnoeud.get_nb();
1077     for (int k=0;k<nb;k++)
1078     {
1079     double *xyztmp=lstnoeud.get(k)->get_coord();
1080     OT_VECTEUR_3D vec(xyztmp,xyz2);
1081     double longueur=vec.get_longueur();
1082     if (longueur<1e-6*rayonrecherche)
1083     {
1084     noeudref2=lstnoeud.get(k);
1085     lstnoeud.supprimer(lstnoeud.get(k));
1086     break;
1087     }
1088     }
1089     }
1090     if (noeudref3==NULL)
1091     {
1092     int nb=lstnoeud.get_nb();
1093     for (int k=0;k<nb;k++)
1094     {
1095     double *xyztmp=lstnoeud.get(k)->get_coord();
1096     OT_VECTEUR_3D vec(xyztmp,xyz3);
1097     double longueur=vec.get_longueur();
1098     if (longueur<1e-6*rayonrecherche)
1099     {
1100     noeudref3=lstnoeud.get(k);
1101     lstnoeud.supprimer(lstnoeud.get(k));
1102     break;
1103     }
1104     }
1105     }
1106     if (noeudref4==NULL)
1107     {
1108     int nb=lstnoeud.get_nb();
1109     for (int k=0;k<nb;k++)
1110     {
1111     double *xyztmp=lstnoeud.get(k)->get_coord();
1112     OT_VECTEUR_3D vec(xyztmp,xyz4);
1113     double longueur=vec.get_longueur();
1114     if (longueur<1e-6*rayonrecherche)
1115     {
1116     noeudref4=lstnoeud.get(k);
1117     lstnoeud.supprimer(lstnoeud.get(k));
1118     break;
1119     }
1120     }
1121     }
1122 francois 54 }
1123 francois 155 if (noeudref1==NULL)
1124     {
1125     noeudref1=mg_maillage->ajouter_mg_noeud(mgvol,xyz1[0],xyz1[1],xyz1[2],IMPOSE);
1126 francois 156 if (tet->get_noeud1()->get_nouveau_numero()!=(-1))
1127 francois 155 {
1128     noeudref1->change_nouveau_numero(tet->get_noeud1()->get_nouveau_numero());
1129     pair<long,long> tmp(noeudref1->get_nouveau_numero(),noeudref1->get_id());
1130     lstcorrespondancenoeud.insert(tmp);
1131     }
1132     }
1133     if (noeudref2==NULL)
1134     {
1135     noeudref2=mg_maillage->ajouter_mg_noeud(mgvol,xyz2[0],xyz2[1],xyz2[2],IMPOSE);
1136 francois 156 if (tet->get_noeud2()->get_nouveau_numero()!=(-1))
1137 francois 155 {
1138     noeudref2->change_nouveau_numero(tet->get_noeud2()->get_nouveau_numero());
1139     pair<long,long> tmp(noeudref2->get_nouveau_numero(),noeudref2->get_id());
1140     lstcorrespondancenoeud.insert(tmp);
1141     }
1142     }
1143     if (noeudref3==NULL)
1144     {
1145     noeudref3=mg_maillage->ajouter_mg_noeud(mgvol,xyz3[0],xyz3[1],xyz3[2],IMPOSE);
1146 francois 156 if (tet->get_noeud3()->get_nouveau_numero()!=(-1))
1147 francois 155 {
1148     noeudref3->change_nouveau_numero(tet->get_noeud3()->get_nouveau_numero());
1149     pair<long,long> tmp(noeudref3->get_nouveau_numero(),noeudref3->get_id());
1150     lstcorrespondancenoeud.insert(tmp);
1151     }
1152     }
1153     if (noeudref4==NULL)
1154     {
1155     noeudref4=mg_maillage->ajouter_mg_noeud(mgvol,xyz4[0],xyz4[1],xyz4[2],IMPOSE);
1156 francois 156 if (tet->get_noeud4()->get_nouveau_numero()!=(-1))
1157 francois 155 {
1158     noeudref4->change_nouveau_numero(tet->get_noeud4()->get_nouveau_numero());
1159     pair<long,long> tmp(noeudref4->get_nouveau_numero(),noeudref4->get_id());
1160     lstcorrespondancenoeud.insert(tmp);
1161     }
1162     }
1163 francois 54 M3D_TRIANGLE* triangle1=(M3D_TRIANGLE*)mg_maillage->get_mg_triangle(noeudref1->get_id(),noeudref3->get_id(),noeudref2->get_id());
1164     M3D_TRIANGLE* triangle2=(M3D_TRIANGLE*)mg_maillage->get_mg_triangle(noeudref1->get_id(),noeudref2->get_id(),noeudref4->get_id());
1165     M3D_TRIANGLE* triangle3=(M3D_TRIANGLE*)mg_maillage->get_mg_triangle(noeudref2->get_id(),noeudref3->get_id(),noeudref4->get_id());
1166     M3D_TRIANGLE* triangle4=(M3D_TRIANGLE*)mg_maillage->get_mg_triangle(noeudref1->get_id(),noeudref4->get_id(),noeudref3->get_id());
1167     if (triangle1==NULL) triangle1=(M3D_TRIANGLE*)insere_triangle(mgvol,noeudref1,noeudref3,noeudref2,IMPOSE);
1168     if (triangle2==NULL) triangle2=(M3D_TRIANGLE*)insere_triangle(mgvol,noeudref1,noeudref2,noeudref4,IMPOSE);
1169     if (triangle3==NULL) triangle3=(M3D_TRIANGLE*)insere_triangle(mgvol,noeudref2,noeudref3,noeudref4,IMPOSE);
1170     if (triangle4==NULL) triangle4=(M3D_TRIANGLE*)insere_triangle(mgvol,noeudref1,noeudref4,noeudref3,IMPOSE);
1171     MG_TETRA* tetra=cree_tetra(mgvol,noeudref1,noeudref2,noeudref3,noeudref4,triangle1,triangle2,triangle3,triangle4,IMPOSE);
1172 francois 194 tetra->change_nouveau_numero(tet->get_nouveau_numero());
1173     mg_maillage->ajouter_mg_tetra(tetra);
1174 francois 54 triangle1->ajouter_voisin(tetra);
1175     triangle2->ajouter_voisin(tetra);
1176     triangle3->ajouter_voisin(tetra);
1177     triangle4->ajouter_voisin(tetra);
1178     MG_FRONT_3D *ftref1=NULL;
1179     MG_FRONT_3D *ftref2=NULL;
1180     MG_FRONT_3D *ftref3=NULL;
1181     MG_FRONT_3D *ftref4=NULL;
1182     if (triangle1->get_etat_front()==FRONT_ACTIF) ftref1=triangle1->get_mgfront();
1183     if (triangle2->get_etat_front()==FRONT_ACTIF) ftref2=triangle2->get_mgfront();
1184     if (triangle3->get_etat_front()==FRONT_ACTIF) ftref3=triangle3->get_mgfront();
1185     if (triangle4->get_etat_front()==FRONT_ACTIF) ftref4=triangle4->get_mgfront();
1186     TPL_MAP_ENTITE<MG_FRONT_3D*> lstfront;
1187     MG_FRONT_3D *ft1=NULL;
1188     MG_FRONT_3D *ft2=NULL;
1189     MG_FRONT_3D *ft3=NULL;
1190     MG_FRONT_3D *ft4=NULL;
1191     if (ftref1==NULL) ft1=mise_a_jour_front(tetra,0,NONFORCE);
1192     else
1193     {
1194     lstfront.ajouter(ftref1->get_front_voisin(0));
1195     lstfront.ajouter(ftref1->get_front_voisin(1));
1196     lstfront.ajouter(ftref1->get_front_voisin(2));
1197 francois 129 supprimer_front_en_avancant_sans_delete(ftref1);
1198 francois 54 }
1199     if (ftref2==NULL) ft2=mise_a_jour_front(tetra,1,NONFORCE);
1200     else
1201     {
1202     lstfront.ajouter(ftref2->get_front_voisin(0));
1203     lstfront.ajouter(ftref2->get_front_voisin(1));
1204     lstfront.ajouter(ftref2->get_front_voisin(2));
1205 francois 129 supprimer_front_en_avancant_sans_delete(ftref2);
1206 francois 54 }
1207     if (ftref3==NULL) ft3=mise_a_jour_front(tetra,2,NONFORCE);
1208     else
1209     {
1210     lstfront.ajouter(ftref3->get_front_voisin(0));
1211     lstfront.ajouter(ftref3->get_front_voisin(1));
1212     lstfront.ajouter(ftref3->get_front_voisin(2));
1213 francois 129 supprimer_front_en_avancant_sans_delete(ftref3);
1214     }
1215 francois 54 if (ftref4==NULL) ft4=mise_a_jour_front(tetra,3,NONFORCE);
1216     else
1217     {
1218     lstfront.ajouter(ftref4->get_front_voisin(0));
1219     lstfront.ajouter(ftref4->get_front_voisin(1));
1220     lstfront.ajouter(ftref4->get_front_voisin(2));
1221 francois 129 supprimer_front_en_avancant_sans_delete(ftref4);
1222     }
1223     if (ftref1!=NULL) {lstfront.supprimer(ftref1);delete ftref1;}
1224     if (ftref2!=NULL) {lstfront.supprimer(ftref2);delete ftref2;}
1225     if (ftref3!=NULL) {lstfront.supprimer(ftref3);delete ftref3;}
1226     if (ftref4!=NULL) {lstfront.supprimer(ftref4);delete ftref4;}
1227 francois 54 if (ft1!=NULL)
1228     {
1229     mise_a_jour_voisin_front(ft1);
1230     lstfront.ajouter(ft1->get_front_voisin(0));
1231     lstfront.ajouter(ft1->get_front_voisin(1));
1232     lstfront.ajouter(ft1->get_front_voisin(2));
1233     }
1234     if (ft2!=NULL)
1235     {
1236     mise_a_jour_voisin_front(ft2);
1237     lstfront.ajouter(ft2->get_front_voisin(0));
1238     lstfront.ajouter(ft2->get_front_voisin(1));
1239     lstfront.ajouter(ft2->get_front_voisin(2));
1240     }
1241     if (ft3!=NULL)
1242     {
1243     mise_a_jour_voisin_front(ft3);
1244     lstfront.ajouter(ft3->get_front_voisin(0));
1245     lstfront.ajouter(ft3->get_front_voisin(1));
1246     lstfront.ajouter(ft3->get_front_voisin(2));
1247     }
1248     if (ft4!=NULL)
1249     {
1250     mise_a_jour_voisin_front(ft4);
1251     lstfront.ajouter(ft4->get_front_voisin(0));
1252     lstfront.ajouter(ft4->get_front_voisin(1));
1253     lstfront.ajouter(ft4->get_front_voisin(2));
1254     }
1255    
1256     for (int l=0;l<lstfront.get_nb();l++) mise_a_jour_voisin_front(lstfront.get(l));
1257    
1258     }
1259     }
1260    
1261     void MAILLEUR3D::insere_contrainte_triangle(MG_VOLUME* mgvol,TPL_LISTE_ENTITE<MG_TRIANGLE*> *lsttri)
1262     {
1263     int nbtri=lsttri->get_nb();
1264     for (int i=0;i<nbtri;i++)
1265     {
1266     M3D_TRIANGLE* tri=(M3D_TRIANGLE*)lsttri->get(i);
1267     double *xyz1=tri->get_noeud1()->get_coord();
1268     double *xyz2=tri->get_noeud2()->get_coord();
1269     double *xyz3=tri->get_noeud3()->get_coord();
1270     double xyz[3];
1271     xyz[0]=(xyz1[0]+xyz2[0]+xyz3[0])/3.;
1272     xyz[1]=(xyz1[1]+xyz2[1]+xyz3[1])/3.;
1273     xyz[2]=(xyz1[2]+xyz2[2]+xyz3[2])/3.;
1274     TPL_MAP_ENTITE<MG_FRONT_3D*> liste_trouvee;
1275     OT_VECTEUR_3D vec1(xyz,xyz1);
1276     OT_VECTEUR_3D vec2(xyz,xyz2);
1277     OT_VECTEUR_3D vec3(xyz,xyz3);
1278     double rayonrecherche=vec1.get_longueur();
1279     if (vec2.get_longueur()>rayonrecherche) rayonrecherche=vec2.get_longueur();
1280     if (vec3.get_longueur()>rayonrecherche) rayonrecherche=vec3.get_longueur();
1281     rayonrecherche=rayonrecherche*1.1;
1282     octree_de_front->rechercher(xyz[0],xyz[2],xyz[2],rayonrecherche,liste_trouvee);
1283     MG_NOEUD* noeudref1=NULL;
1284     MG_NOEUD* noeudref2=NULL;
1285     MG_NOEUD* noeudref3=NULL;
1286     int nbtrouve=liste_trouvee.get_nb();
1287     for (int j=0;j<nbtrouve;j++)
1288     {
1289     MG_FRONT_3D* ft=liste_trouvee.get(j);
1290     TPL_MAP_ENTITE<MG_NOEUD*> lstnoeud;
1291     lstnoeud.ajouter(ft->get_noeud1());
1292     lstnoeud.ajouter(ft->get_noeud2());
1293     lstnoeud.ajouter(ft->get_noeud3());
1294     int nb=lstnoeud.get_nb();
1295     for (int k=0;k<nb;k++)
1296     {
1297     double *xyztmp=lstnoeud.get(k)->get_coord();
1298     OT_VECTEUR_3D vec(xyztmp,xyz1);
1299     double longueur=vec.get_longueur();
1300     if (longueur<1e-10*rayonrecherche)
1301     {
1302     noeudref1=lstnoeud.get(k);
1303     lstnoeud.supprimer(lstnoeud.get(k));
1304     break;
1305     }
1306     }
1307     nb=lstnoeud.get_nb();
1308     for (int k=0;k<nb;k++)
1309     {
1310     double *xyztmp=lstnoeud.get(k)->get_coord();
1311     OT_VECTEUR_3D vec(xyztmp,xyz2);
1312     double longueur=vec.get_longueur();
1313     if (longueur<1e-10*rayonrecherche)
1314     {
1315     noeudref2=lstnoeud.get(k);
1316     lstnoeud.supprimer(lstnoeud.get(k));
1317     break;
1318     }
1319     }
1320     nb=lstnoeud.get_nb();
1321     for (int k=0;k<nb;k++)
1322     {
1323     double *xyztmp=lstnoeud.get(k)->get_coord();
1324     OT_VECTEUR_3D vec(xyztmp,xyz3);
1325     double longueur=vec.get_longueur();
1326     if (longueur<1e-10*rayonrecherche)
1327     {
1328     noeudref3=lstnoeud.get(k);
1329     lstnoeud.supprimer(lstnoeud.get(k));
1330     break;
1331     }
1332     }
1333    
1334     }
1335 francois 155 if (noeudref1==NULL)
1336     {
1337     noeudref1=mg_maillage->ajouter_mg_noeud(mgvol,xyz1[0],xyz1[1],xyz1[2],IMPOSE);
1338     noeudref1->change_nouveau_numero(tri->get_noeud1()->get_nouveau_numero());
1339     }
1340     if (noeudref2==NULL)
1341     {
1342     noeudref2=mg_maillage->ajouter_mg_noeud(mgvol,xyz2[0],xyz2[1],xyz2[2],IMPOSE);
1343     noeudref2->change_nouveau_numero(tri->get_noeud2()->get_nouveau_numero());
1344     }
1345     if (noeudref3==NULL)
1346     {
1347     noeudref3=mg_maillage->ajouter_mg_noeud(mgvol,xyz3[0],xyz3[1],xyz3[2],IMPOSE);
1348     noeudref3->change_nouveau_numero(tri->get_noeud3()->get_nouveau_numero());
1349     }
1350    
1351 francois 54 MG_SEGMENT* segment1=mg_maillage->get_mg_segment(noeudref1->get_id(),noeudref2->get_id());
1352     MG_SEGMENT* segment2=mg_maillage->get_mg_segment(noeudref1->get_id(),noeudref3->get_id());
1353     MG_SEGMENT* segment3=mg_maillage->get_mg_segment(noeudref2->get_id(),noeudref3->get_id());
1354     if (segment1==NULL) segment1=cree_segment(mgvol,noeudref1,noeudref2,IMPOSE);
1355     if (segment2==NULL) segment2=cree_segment(mgvol,noeudref1,noeudref3,IMPOSE);
1356     if (segment3==NULL) segment3=cree_segment(mgvol,noeudref2,noeudref3,IMPOSE);
1357     MG_TRIANGLE* nvtri=cree_triangle(mgvol,noeudref1,noeudref2,noeudref3,segment1,segment2,segment3,IMPOSE);
1358     MG_TRIANGLE* nvtri2=cree_triangle(mgvol,noeudref1,noeudref3,noeudref2,segment1,segment2,segment3,IMPOSE);
1359     //MG_TRIANGLE* nvtri=insere_triangle(mgvol,noeudref1,noeudref2,noeudref3,IMPOSE);
1360     //MG_TRIANGLE* nvtri2=insere_triangle(mgvol,noeudref1,noeudref3,noeudref2,IMPOSE);
1361     MG_FRONT_3D *ft=ajouter_front_courant(NONFORCE,noeudref1,noeudref2,noeudref3,nvtri);
1362     MG_FRONT_3D *ft2=ajouter_front_courant(NONFORCE,noeudref1,noeudref3,noeudref2,nvtri2);
1363     mise_a_jour_voisin_front(ft);
1364     mise_a_jour_voisin_front(ft2);
1365     mise_a_jour_voisin_front(ft->get_front_voisin(0));
1366     mise_a_jour_voisin_front(ft->get_front_voisin(1));
1367     mise_a_jour_voisin_front(ft->get_front_voisin(2));
1368     mise_a_jour_voisin_front(ft2->get_front_voisin(0));
1369     mise_a_jour_voisin_front(ft2->get_front_voisin(1));
1370     mise_a_jour_voisin_front(ft2->get_front_voisin(2));
1371    
1372     }
1373    
1374    
1375    
1376    
1377     }
1378