ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/mailleur3d_outil.cpp
Revision: 283
Committed: Tue Sep 13 21:11:20 2011 UTC (13 years, 8 months ago) by francois
File size: 55234 byte(s)
Log Message:
structure de l'écriture

File Contents

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