ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/mailleur3d_outil.cpp
Revision: 54
Committed: Fri Sep 28 19:43:11 2007 UTC (17 years, 7 months ago) by francois
Original Path: magic/lib/mailleur/mailleur/src/mailleur3d_outil.cpp
File size: 53893 byte(s)
Log Message:

File Contents

# User Rev Content
1 5 //------------------------------------------------------------
2     //------------------------------------------------------------
3     // MAGiC
4     // Jean Christophe Cuillière et Vincent FRANCOIS
5     // Département de Génie Mécanique - UQTR
6     //------------------------------------------------------------
7     // Le projet MAGIC est un projet de recherche du département
8     // de génie mécanique de l'Université du Québec à
9     // Trois Rivières
10     // Les librairies ne peuvent être utilisées sans l'accord
11     // des auteurs (contact : francois@uqtr.ca)
12     //------------------------------------------------------------
13     //------------------------------------------------------------
14     //
15     // mailleur3d_outil.cpp
16     //
17     //------------------------------------------------------------
18     //------------------------------------------------------------
19     // COPYRIGHT 2000
20     // Version du 02/03/2006 à 11H23
21     //------------------------------------------------------------
22     //------------------------------------------------------------
23    
24    
25     #include "gestionversion.h"
26    
27    
28     //#pragma hdrstop
29    
30     #include <math.h>
31     #include <time.h>
32     #include "mailleur3d.h"
33 francois 35 #include "mailleur3d_front.h"
34 5 #include "ot_mathematique.h"
35 francois 54 #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     // validité de la solution
245     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     double xyz4[3]={xyzg[0]+hauteur*normal.get_x(),xyzg[1]+hauteur*normal.get_y(),xyzg[2]+hauteur*normal.get_z()};
521     double l1=calcul_distance_metrique(xyz1,xyz2);
522     double l2=calcul_distance_metrique(xyz1,xyz3);
523     double l3=calcul_distance_metrique(xyz2,xyz3);
524     double l=0.333333333333333*(l1+l2+l3);
525     double longueur_desiree=0.81649658*(priorite_metrique+(1-priorite_metrique)*l);
526     longueur_desiree=longueur_desiree/ft->get_ifail();
527     ajuste_distance_metrique(xyzg,xyz4,longueur_desiree);
528     OT_VECTEUR_3D hauteurajustee(xyzg[0]-xyz4[0],xyzg[1]-xyz4[1],xyzg[2]-xyz4[2]);
529     hauteur=hauteurajustee.get_longueur();
530 5
531     // recherche des voisins
532     TPL_MAP_ENTITE<MG_FRONT_3D*> liste_front;
533     TPL_MAP_ENTITE<MG_FRONT_3D*> liste_front_proche;
534 francois 35 octree_de_front->rechercher(xyz4[0],xyz4[1],xyz4[2],hauteur*1.837117307,liste_front);
535 5 int nb_liste=liste_front.get_nb();
536 francois 35 double disref=(priorite_metrique+(1-priorite_metrique)*l)/ft->get_ifail();
537 5 MG_NOEUD* noeudref=NULL;
538     for (int i=0;i<nb_liste;i++)
539     {
540     MG_FRONT_3D* ft2=liste_front.get(i);
541     for (int j=0;j<3;j++)
542     {
543     MG_NOEUD* noeud_tmp;
544     if (j==0) noeud_tmp=ft2->get_noeud1();
545     if (j==1) noeud_tmp=ft2->get_noeud2();
546     if (j==2) noeud_tmp=ft2->get_noeud3();
547 francois 35 double *xyz_tmp=noeud_tmp->get_coord();
548     double dis=calcul_distance_metrique(xyz_tmp,xyz4,1);
549 5 if (dis<disref)
550     {
551     disref=dis;
552     noeudref=noeud_tmp;
553     }
554     }
555     if (ft2!=ft)
556     {
557 francois 35 double dismin=evaluer_distance_noeud_triangle(xyz4[0],xyz4[1],xyz4[2],ft2->get_triangle());
558     if (dismin<hauteur) liste_front_proche.ajouter(ft2);
559 5 }
560     }
561     // il y a un noeud proche
562     if (noeudref!=NULL)
563 francois 35 {
564     double dis2=calcul_distance_metrique(noeudref->get_coord(),xyzg,2);
565     if (dis2<longueur_desiree*1.5)
566 5 if (noeudref!=noeud1)
567     if (noeudref!=noeud2)
568     if (noeudref!=noeud3)
569     {
570     liste_noeud.ajouter(noeudref);
571     clock_t cpu2=clock(); //debug
572     cpu[8]=cpu[8]+(cpu2-cpu1)/CLK_TCK;
573     nbfois[8]++;
574     return EXISTE;
575     }
576 francois 35 }
577 5 // il n y a pas de noeud proche mais il y a une face proche
578     int nb_proche=liste_front_proche.get_nb();
579 francois 35 int nb_cas=liste_cas.size();
580 5 if (nb_proche>0)
581     {
582     for (int i=0;i<nb_proche;i++)
583     {
584     MG_FRONT_3D* fttmp=liste_front_proche.get(i);
585 francois 35 int ok=0;
586     for (int j=0;j<nb_cas;j++)
587     {
588     MG_NOEUD* casnoeud=liste_cas[j]->mgnoeud;
589     if ((fttmp->get_noeud1()==casnoeud) || (fttmp->get_noeud2()==casnoeud) || (fttmp->get_noeud3()==casnoeud) )
590     {
591     ok=1;
592     break;
593     }
594     }
595     if (ok==1) continue;
596 5 MG_NOEUD* noeudtmp=fttmp->get_noeud1();
597     if (noeudtmp!=noeud1)
598     if (noeudtmp!=noeud2)
599     if (noeudtmp!=noeud3) liste_noeud.ajouter(noeudtmp);
600     noeudtmp=fttmp->get_noeud2();
601     if (noeudtmp!=noeud1)
602     if (noeudtmp!=noeud2)
603     if (noeudtmp!=noeud3) liste_noeud.ajouter(noeudtmp);
604     noeudtmp=fttmp->get_noeud3();
605     if (noeudtmp!=noeud1)
606     if (noeudtmp!=noeud2)
607     if (noeudtmp!=noeud3) liste_noeud.ajouter(noeudtmp);
608    
609     }
610     clock_t cpu2=clock(); //debug
611     cpu[8]=cpu[8]+(cpu2-cpu1)/CLK_TCK;
612     nbfois[8]++;
613 francois 35 if (liste_noeud.get_nb()>0) return MULTIEXISTE;
614 5 }
615     // il faut creer un noeud c est parfait c est qu on veut.
616 francois 54 M3D_NOEUD* nouveau_noeud=new M3D_NOEUD(mgvol,xyz4[0],xyz4[1],xyz4[2],MAILLEUR_AUTO);
617 5 mg_maillage->ajouter_mg_noeud(nouveau_noeud);
618     liste_noeud.ajouter(nouveau_noeud);
619     clock_t cpu2=clock(); //debug
620     cpu[8]=cpu[8]+(cpu2-cpu1)/CLK_TCK;
621     nbfois[8]++;
622     return CREATION;
623     }
624    
625    
626    
627     double MAILLEUR3D::evaluer_distance_noeud_triangle(double x,double y,double z,MG_TRIANGLE* tri)
628     {
629     double eps=0.0001;
630     MG_NOEUD* noeud1=tri->get_noeud1();
631     MG_NOEUD* noeud2=tri->get_noeud2();
632     MG_NOEUD* noeud3=tri->get_noeud3();
633     OT_VECTEUR_3D n1n2(noeud2->get_x()-noeud1->get_x(),noeud2->get_y()-noeud1->get_y(),noeud2->get_z()-noeud1->get_z());
634     OT_VECTEUR_3D n1n3(noeud3->get_x()-noeud1->get_x(),noeud3->get_y()-noeud1->get_y(),noeud3->get_z()-noeud1->get_z());
635     OT_VECTEUR_3D n1n4(x-noeud1->get_x(),y-noeud1->get_y(),z-noeud1->get_z());
636     OT_VECTEUR_3D n1n=n1n2&n1n3;
637     OT_MATRICE_3D matbase(n1n2,n1n3,n1n);
638     OT_MATRICE_3D matxsi(n1n4,n1n3,n1n);
639     OT_MATRICE_3D mateta(n1n2,n1n4,n1n);
640     double xsi=matxsi.get_determinant()/matbase.get_determinant();
641     double eta=mateta.get_determinant()/matbase.get_determinant();
642     double dis;
643     if ((eta>-eps) && (xsi>-eps) && (eta+xsi<1.0+eps))
644     {
645     double xproj=(1-xsi-eta)*noeud1->get_x()+xsi*noeud2->get_x()+eta*noeud3->get_x();
646     double yproj=(1-xsi-eta)*noeud1->get_y()+xsi*noeud2->get_y()+eta*noeud3->get_y();
647     double zproj=(1-xsi-eta)*noeud1->get_z()+xsi*noeud2->get_z()+eta*noeud3->get_z();
648     OT_VECTEUR_3D vec(x-xproj,y-yproj,z-zproj);
649     dis=vec.get_longueur();
650     }
651     else
652     {
653     OT_VECTEUR_3D nn1(x-noeud1->get_x(),y-noeud1->get_y(),z-noeud1->get_z());
654     OT_VECTEUR_3D nn2(x-noeud2->get_x(),y-noeud2->get_y(),z-noeud2->get_z());
655     OT_VECTEUR_3D nn3(x-noeud3->get_x(),y-noeud3->get_y(),z-noeud3->get_z());
656     dis=0.3333333333333*(nn1.get_longueur()+nn2.get_longueur()+nn3.get_longueur());
657     }
658     return dis;
659     }
660    
661    
662     void MAILLEUR3D::detruit_tetra(MG_TRIANGLE* tri)
663     {
664     clock_t cpu1=clock(); //debug
665    
666     M3D_TRIANGLE* mtri;
667     //if (tri->get_type_entite()==IDMG_TRIANGLE) return;
668     //else mtri=(M3D_TRIANGLE*)tri;
669     mtri=(M3D_TRIANGLE*)tri;
670    
671     M3D_TETRA* tet=(M3D_TETRA*)mtri->get_voisin(0);
672 francois 54 M3D_TRIANGLE* mtriangle1=(M3D_TRIANGLE*)tet->get_triangle1();
673     M3D_TRIANGLE* mtriangle2=(M3D_TRIANGLE*)tet->get_triangle2();
674     M3D_TRIANGLE* mtriangle3=(M3D_TRIANGLE*)tet->get_triangle3();
675     M3D_TRIANGLE* mtriangle4=(M3D_TRIANGLE*)tet->get_triangle4();
676    
677    
678 5 tet->etat_detruit();
679     TPL_MAP_ENTITE<MG_FRONT_3D*> liste_de_front;
680     TPL_MAP_ENTITE<MG_FRONT_3D*> liste_de_front2;
681     MG_FRONT_3D* fttet1=NULL;
682     MG_FRONT_3D* fttet2=NULL;
683     MG_FRONT_3D* fttet3=NULL;
684     MG_FRONT_3D* fttet4=NULL;
685    
686 francois 54 if (mtriangle1->get_etat_front()==FRONT_ACTIF)
687 5 {
688 francois 54 MG_FRONT_3D* fttmp=mtriangle1->get_mgfront();
689 5 fttet1=fttmp;
690     liste_de_front.ajouter(fttmp->get_front_voisin(0));
691     liste_de_front.ajouter(fttmp->get_front_voisin(1));
692     liste_de_front.ajouter(fttmp->get_front_voisin(2));
693     //supprimer_front_en_reculant(fttmp);
694 francois 54 mtriangle1->supprimer_voisin(tet);
695 5 }
696     else
697     {
698 francois 54 MG_FRONT_3D* ftn1=ajouter_front_courant(PRIORITAIRE,tet->get_noeud1(),tet->get_noeud2(),tet->get_noeud3(),mtriangle1);
699 5 ftn1->incremente_ifail();
700 francois 54 mtriangle1->supprimer_voisin(tet);
701 5 liste_de_front2.ajouter(ftn1);
702     }
703 francois 54 if (mtriangle2->get_etat_front()==FRONT_ACTIF)
704 5 {
705 francois 54 MG_FRONT_3D* fttmp=mtriangle2->get_mgfront();
706 5 fttet2=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 francois 54 mtriangle2->supprimer_voisin(tet);
712 5 }
713     else
714     {
715 francois 54 MG_FRONT_3D* ftn1=ajouter_front_courant(PRIORITAIRE,tet->get_noeud1(),tet->get_noeud4(),tet->get_noeud2(),mtriangle2);
716 5 ftn1->incremente_ifail();
717 francois 54 mtriangle2->supprimer_voisin(tet);
718 5 liste_de_front2.ajouter(ftn1);
719     }
720 francois 54 if (mtriangle3->get_etat_front()==FRONT_ACTIF)
721 5 {
722 francois 54 MG_FRONT_3D* fttmp=mtriangle3->get_mgfront();
723 5 fttet3=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 francois 54 mtriangle3->supprimer_voisin(tet);
729 5 }
730     else
731     {
732 francois 54 MG_FRONT_3D* ftn1=ajouter_front_courant(PRIORITAIRE,tet->get_noeud2(),tet->get_noeud4(),tet->get_noeud3(),mtriangle3);
733 5 ftn1->incremente_ifail();
734 francois 54 mtriangle3->supprimer_voisin(tet);
735 5 liste_de_front2.ajouter(ftn1);
736     }
737 francois 54 if (mtriangle4->get_etat_front()==FRONT_ACTIF)
738 5 {
739 francois 54 MG_FRONT_3D* fttmp=mtriangle4->get_mgfront();
740 5 fttet4=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 francois 54 mtriangle4->supprimer_voisin(tet);
746 5 }
747     else
748     {
749 francois 54 MG_FRONT_3D* ftn1=ajouter_front_courant(PRIORITAIRE,tet->get_noeud1(),tet->get_noeud3(),tet->get_noeud4(),mtriangle4);
750 5 ftn1->incremente_ifail();
751 francois 54 mtriangle4->supprimer_voisin(tet);
752 5 liste_de_front2.ajouter(ftn1);
753     }
754    
755     if (fttet1!=NULL)
756     {
757     liste_de_front.supprimer(fttet1);
758     supprimer_front_en_reculant(fttet1);
759     }
760     if (fttet2!=NULL)
761     {
762     liste_de_front.supprimer(fttet2);
763     supprimer_front_en_reculant(fttet2);
764     }
765     if (fttet3!=NULL)
766     {
767     liste_de_front.supprimer(fttet3);
768     supprimer_front_en_reculant(fttet3);
769     }
770     if (fttet4!=NULL)
771     {
772     liste_de_front.supprimer(fttet4);
773     supprimer_front_en_reculant(fttet4);
774     }
775     int nb_liste1=liste_de_front.get_nb();
776     for (int i=0;i<nb_liste1;i++)
777     if (liste_de_front.get(i)!=NULL)
778     mise_a_jour_voisin_front(liste_de_front.get(i));
779     int nb_liste2=liste_de_front2.get_nb();
780     for (int i=0;i<nb_liste2;i++)
781     {
782     mise_a_jour_voisin_front(liste_de_front2.get(i));
783     mise_a_jour_voisin_front(liste_de_front2.get(i)->get_front_voisin(0));
784     mise_a_jour_voisin_front(liste_de_front2.get(i)->get_front_voisin(1));
785     mise_a_jour_voisin_front(liste_de_front2.get(i)->get_front_voisin(2));
786     }
787    
788     for (int k=0;k<4;k++)
789     {
790     MG_NOEUD* noeud;
791     if (k==0) noeud=tet->get_noeud1();
792     if (k==1) noeud=tet->get_noeud2();
793     if (k==2) noeud=tet->get_noeud3();
794     if (k==3) noeud=tet->get_noeud4();
795     int nb_tet=noeud->get_lien_tetra()->get_nb();
796    
797     if (noeud->get_type_entite()==IDM3D_NOEUD)
798     {
799     int test=false;
800     M3D_NOEUD* mnoeud=(M3D_NOEUD*)noeud;
801     for (int i=0;i<nb_tet;i++)
802     {
803     MG_TETRA* tet=noeud->get_lien_tetra()->get(i);
804     if (tet->get_type_entite()==IDM3D_TETRA)
805     {
806     M3D_TETRA* mtet=(M3D_TETRA*)tet;
807     if (mtet->get_etat()==ACTIF) test=true;
808     }
809     else test=true;
810     }
811     if (test==false) mnoeud->etat_detruit();
812     }
813     }
814     clock_t cpu2=clock();
815     cpu[9]=cpu[9]+(cpu2-cpu1)/CLK_TCK;
816     nbfois[9]++;
817    
818     }
819    
820 francois 35 double MAILLEUR3D::calcul_distance_metrique(MG_SEGMENT* seg,int pas)
821     {
822     double *xyz1=seg->get_noeud1()->get_coord();
823     double *xyz2=seg->get_noeud2()->get_coord();
824     return calcul_distance_metrique(xyz1,xyz2,pas);
825     }
826 5
827 francois 35 double MAILLEUR3D::calcul_distance_metrique(double *xyz1,double* xyz2,int pas)
828     {
829     double longueur=0.;
830     double dxyz[3];
831     dxyz[0]=xyz2[0]-xyz1[0];
832     dxyz[1]=xyz2[1]-xyz1[1];
833     dxyz[2]=xyz2[2]-xyz1[2];
834     for (int i=0;i<pas;i++)
835     {
836     double ti=1.0*i/pas;
837     double tii=1.0*(i+1)/pas;
838     double tgauss1=0.7886751345*tii+0.2113248654*ti;
839     double xyz[3];
840     xyz[0]=xyz1[0]+tgauss1*(xyz2[0]-xyz1[0]);
841     xyz[1]=xyz1[1]+tgauss1*(xyz2[1]-xyz1[1]);
842     xyz[2]=xyz1[2]+tgauss1*(xyz2[2]-xyz1[2]);
843     double tenseur[9];
844     metrique->evaluer(xyz,tenseur);
845     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]);
846     longueur=longueur+0.5*(tii-ti)*val;
847     double tgauss2=0.7886751345*ti+0.2113248654*tii;
848     xyz[0]=xyz1[0]+tgauss2*(xyz2[0]-xyz1[0]);
849     xyz[1]=xyz1[1]+tgauss2*(xyz2[1]-xyz1[1]);
850     xyz[2]=xyz1[2]+tgauss2*(xyz2[2]-xyz1[2]);
851     metrique->evaluer(xyz,tenseur);
852     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]);
853     longueur=longueur+0.5*(tii-ti)*val;
854     }
855     return longueur;
856     }
857 5
858 francois 35 void MAILLEUR3D::ajuste_distance_metrique(double *xyz1,double *xyz2,double longueur_desiree,int pas)
859     {
860     double longueur=0,ancienne_longueur;
861     double ti,tii;
862     double dxyz[3];
863     dxyz[0]=xyz2[0]-xyz1[0];
864     dxyz[1]=xyz2[1]-xyz1[1];
865     dxyz[2]=xyz2[2]-xyz1[2];
866     int i=0;
867     do
868     {
869     ancienne_longueur=longueur;
870     ti=1.0*i/pas;
871     tii=1.0*(i+1)/pas;
872     double tgauss1=0.7886751345*tii+0.2113248654*ti;
873     double xyz[3];
874     xyz[0]=xyz1[0]+tgauss1*(xyz2[0]-xyz1[0]);
875     xyz[1]=xyz1[1]+tgauss1*(xyz2[1]-xyz1[1]);
876     xyz[2]=xyz1[2]+tgauss1*(xyz2[2]-xyz1[2]);
877     double tenseur[9];
878     metrique->evaluer(xyz,tenseur);
879     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]);
880     longueur=longueur+0.5*(tii-ti)*val;
881     double tgauss2=0.7886751345*ti+0.2113248654*tii;
882     xyz[0]=xyz1[0]+tgauss2*(xyz2[0]-xyz1[0]);
883     xyz[1]=xyz1[1]+tgauss2*(xyz2[1]-xyz1[1]);
884     xyz[2]=xyz1[2]+tgauss2*(xyz2[2]-xyz1[2]);
885     metrique->evaluer(xyz,tenseur);
886     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]);
887     longueur=longueur+0.5*(tii-ti)*val;
888     i++;
889     }
890     while (longueur<longueur_desiree);
891     double t=ti+(tii-ti)*(longueur_desiree-ancienne_longueur)/(longueur-ancienne_longueur);
892     xyz2[0]=xyz1[0]+t*(xyz2[0]-xyz1[0]);
893     xyz2[1]=xyz1[1]+t*(xyz2[1]-xyz1[1]);
894     xyz2[2]=xyz1[2]+t*(xyz2[2]-xyz1[2]);
895     }
896    
897 5 void MAILLEUR3D::detruit_element_inactif(void)
898     {
899     clock_t cpu1=clock();
900     int nb_tetra=mg_maillage->get_nb_mg_tetra();
901     int num=0;
902     for (int i=0;i<nb_tetra;i++)
903     {
904     MG_TETRA* tet=mg_maillage->get_mg_tetra(num);
905     M3D_TETRA* mtet=(M3D_TETRA*)tet;
906     if (mtet->get_etat()==DETRUIT)
907     mg_maillage->supprimer_mg_tetra(num);
908     else num++;
909     }
910     int nb_noeud=mg_maillage->get_nb_mg_noeud();
911     num=0;
912     for (int i=0;i<nb_noeud;i++)
913     {
914     MG_NOEUD* noeud=mg_maillage->get_mg_noeud(num);
915     if (noeud->get_type_entite()==IDM3D_NOEUD)
916     {
917     M3D_NOEUD* mnoeud=(M3D_NOEUD*)noeud;
918     if (mnoeud->get_etat()==DETRUIT)
919     mg_maillage->supprimer_mg_noeud(num);
920     else num++;
921     }
922     else num++;
923     }
924    
925     clock_t cpu2=clock();
926     cpu[10]=cpu[10]+(cpu2-cpu1)/CLK_TCK;
927     nbfois[10]++;
928     }
929    
930 francois 35
931    
932     double MAILLEUR3D::calcul_volume_tetra_metrique(MG_TETRA* tet)
933     {
934     double *xyz1=tet->get_noeud1()->get_coord();
935     double *xyz2=tet->get_noeud2()->get_coord();
936     double *xyz3=tet->get_noeud3()->get_coord();
937     double *xyz4=tet->get_noeud4()->get_coord();
938     /*double xyzm[3];
939     double xyzg[3];
940     xyzm[0]=(xyz1[0]+xyz2[0])/2.;
941     xyzm[1]=(xyz1[1]+xyz2[1])/2.;
942     xyzm[2]=(xyz1[2]+xyz2[2])/2.;
943     xyzg[0]=(xyz1[0]+xyz2[0]+xyz3[0])/3.;
944     xyzg[1]=(xyz1[1]+xyz2[1]+xyz3[1])/3.;
945     xyzg[2]=(xyz1[2]+xyz2[2]+xyz3[2])/3.;
946     double base=calcul_distance_metrique(xyz1,xyz2);
947     double hauteur=calcul_distance_metrique(xyzm,xyz3);
948     double hauteur2=calcul_distance_metrique(xyzg,xyz4);
949     double vol=base*hauteur*hauteur2*1.414213562;*/
950     double a=0.1381966012;
951     double b=0.5854101965;
952     double poids_gauss[4],xsi_gauss[4],eta_gauss[4],zeta_gauss[4];
953     poids_gauss[0]=0.041666666666666;xsi_gauss[0]=a;eta_gauss[0]=a;zeta_gauss[0]=a;
954     poids_gauss[1]=0.041666666666666;xsi_gauss[1]=a;eta_gauss[1]=a;zeta_gauss[1]=b;
955     poids_gauss[2]=0.041666666666666;xsi_gauss[2]=a;eta_gauss[2]=b;zeta_gauss[2]=a;
956     poids_gauss[3]=0.041666666666666;xsi_gauss[3]=b;eta_gauss[3]=a;zeta_gauss[3]=a;
957     OT_VECTEUR_3D vec1(xyz1,xyz2);
958     OT_VECTEUR_3D vec2(xyz1,xyz3);
959     OT_VECTEUR_3D vec3(xyz1,xyz4);
960     double facteur=(vec1&vec2)*vec3;
961     double vol=0.;
962     for (int j=0;j<4;j++)
963     {
964     double xyz[3];
965     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];
966     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];
967     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];
968     double tenseur[9];
969     metrique->evaluer(xyz,tenseur);
970     double dens=1/sqrt(tenseur[0]);
971     dens=dens*dens*dens;
972     vol=vol+poids_gauss[j]*facteur*8.485281374/dens;
973     }
974    
975     return vol;
976     }
977    
978    
979     void MAILLEUR3D::analyse_maillage_obtenu(double &vol)
980     {
981     int nb=mg_maillage->get_nb_mg_tetra();
982     vol=0.;
983     for (int i=0;i<nb;i++)
984     {
985     MG_TETRA* tet=mg_maillage->get_mg_tetra(i);
986     vol=vol+calcul_volume_tetra_metrique(tet);
987     }
988     }
989    
990 francois 54
991     void MAILLEUR3D::insere_contrainte_tetra(MG_VOLUME* mgvol,TPL_LISTE_ENTITE<MG_TETRA*> *lsttet)
992     {
993     int nbtet=lsttet->get_nb();
994     for (int i=0;i<nbtet;i++)
995     {
996     M3D_TETRA* tet=(M3D_TETRA*)lsttet->get(i);
997     double *xyz1=tet->get_noeud1()->get_coord();
998     double *xyz2=tet->get_noeud2()->get_coord();
999     double *xyz3=tet->get_noeud3()->get_coord();
1000     double *xyz4=tet->get_noeud4()->get_coord();
1001     double xyz[3];
1002     xyz[0]=(xyz1[0]+xyz2[0]+xyz3[0]+xyz4[0])*0.25;
1003     xyz[1]=(xyz1[1]+xyz2[1]+xyz3[1]+xyz4[1])*0.25;
1004     xyz[2]=(xyz1[2]+xyz2[2]+xyz3[2]+xyz4[2])*0.25;
1005     TPL_MAP_ENTITE<MG_FRONT_3D*> liste_trouvee;
1006     OT_VECTEUR_3D vec1(xyz,xyz1);
1007     OT_VECTEUR_3D vec2(xyz,xyz2);
1008     OT_VECTEUR_3D vec3(xyz,xyz3);
1009     OT_VECTEUR_3D vec4(xyz,xyz4);
1010     double rayonrecherche=vec1.get_longueur();
1011     if (vec2.get_longueur()>rayonrecherche) rayonrecherche=vec2.get_longueur();
1012     if (vec3.get_longueur()>rayonrecherche) rayonrecherche=vec3.get_longueur();
1013     if (vec4.get_longueur()>rayonrecherche) rayonrecherche=vec4.get_longueur();
1014     rayonrecherche=rayonrecherche*1.1;
1015     octree_de_front->rechercher(xyz[0],xyz[1],xyz[2],rayonrecherche,liste_trouvee);
1016     MG_NOEUD* noeudref1=NULL;
1017     MG_NOEUD* noeudref2=NULL;
1018     MG_NOEUD* noeudref3=NULL;
1019     MG_NOEUD* noeudref4=NULL;
1020     int nbtrouve=liste_trouvee.get_nb();
1021     for (int j=0;j<nbtrouve;j++)
1022     {
1023     MG_FRONT_3D* ft=liste_trouvee.get(j);
1024     TPL_MAP_ENTITE<MG_NOEUD*> lstnoeud;
1025     lstnoeud.ajouter(ft->get_noeud1());
1026     lstnoeud.ajouter(ft->get_noeud2());
1027     lstnoeud.ajouter(ft->get_noeud3());
1028     int nb=lstnoeud.get_nb();
1029     for (int k=0;k<nb;k++)
1030     {
1031     double *xyztmp=lstnoeud.get(k)->get_coord();
1032     OT_VECTEUR_3D vec(xyztmp,xyz1);
1033     double longueur=vec.get_longueur();
1034     if (longueur<1e-10*rayonrecherche)
1035     {
1036     noeudref1=lstnoeud.get(k);
1037     lstnoeud.supprimer(lstnoeud.get(k));
1038     break;
1039     }
1040     }
1041     nb=lstnoeud.get_nb();
1042     for (int k=0;k<nb;k++)
1043     {
1044     double *xyztmp=lstnoeud.get(k)->get_coord();
1045     OT_VECTEUR_3D vec(xyztmp,xyz2);
1046     double longueur=vec.get_longueur();
1047     if (longueur<1e-10*rayonrecherche)
1048     {
1049     noeudref2=lstnoeud.get(k);
1050     lstnoeud.supprimer(lstnoeud.get(k));
1051     break;
1052     }
1053     }
1054     nb=lstnoeud.get_nb();
1055     for (int k=0;k<nb;k++)
1056     {
1057     double *xyztmp=lstnoeud.get(k)->get_coord();
1058     OT_VECTEUR_3D vec(xyztmp,xyz3);
1059     double longueur=vec.get_longueur();
1060     if (longueur<1e-10*rayonrecherche)
1061     {
1062     noeudref3=lstnoeud.get(k);
1063     lstnoeud.supprimer(lstnoeud.get(k));
1064     break;
1065     }
1066     }
1067     nb=lstnoeud.get_nb();
1068     for (int k=0;k<nb;k++)
1069     {
1070     double *xyztmp=lstnoeud.get(k)->get_coord();
1071     OT_VECTEUR_3D vec(xyztmp,xyz4);
1072     double longueur=vec.get_longueur();
1073     if (longueur<1e-10*rayonrecherche)
1074     {
1075     noeudref4=lstnoeud.get(k);
1076     lstnoeud.supprimer(lstnoeud.get(k));
1077     break;
1078     }
1079     }
1080     }
1081     if (noeudref1==NULL) noeudref1=mg_maillage->ajouter_mg_noeud(mgvol,xyz1[0],xyz1[1],xyz1[2],IMPOSE);
1082     if (noeudref2==NULL) noeudref2=mg_maillage->ajouter_mg_noeud(mgvol,xyz2[0],xyz2[1],xyz2[2],IMPOSE);
1083     if (noeudref3==NULL) noeudref3=mg_maillage->ajouter_mg_noeud(mgvol,xyz3[0],xyz3[1],xyz3[2],IMPOSE);
1084     if (noeudref4==NULL) noeudref4=mg_maillage->ajouter_mg_noeud(mgvol,xyz4[0],xyz4[1],xyz4[2],IMPOSE);
1085     M3D_TRIANGLE* triangle1=(M3D_TRIANGLE*)mg_maillage->get_mg_triangle(noeudref1->get_id(),noeudref3->get_id(),noeudref2->get_id());
1086     M3D_TRIANGLE* triangle2=(M3D_TRIANGLE*)mg_maillage->get_mg_triangle(noeudref1->get_id(),noeudref2->get_id(),noeudref4->get_id());
1087     M3D_TRIANGLE* triangle3=(M3D_TRIANGLE*)mg_maillage->get_mg_triangle(noeudref2->get_id(),noeudref3->get_id(),noeudref4->get_id());
1088     M3D_TRIANGLE* triangle4=(M3D_TRIANGLE*)mg_maillage->get_mg_triangle(noeudref1->get_id(),noeudref4->get_id(),noeudref3->get_id());
1089     if (triangle1==NULL) triangle1=(M3D_TRIANGLE*)insere_triangle(mgvol,noeudref1,noeudref3,noeudref2,IMPOSE);
1090     if (triangle2==NULL) triangle2=(M3D_TRIANGLE*)insere_triangle(mgvol,noeudref1,noeudref2,noeudref4,IMPOSE);
1091     if (triangle3==NULL) triangle3=(M3D_TRIANGLE*)insere_triangle(mgvol,noeudref2,noeudref3,noeudref4,IMPOSE);
1092     if (triangle4==NULL) triangle4=(M3D_TRIANGLE*)insere_triangle(mgvol,noeudref1,noeudref4,noeudref3,IMPOSE);
1093     MG_TETRA* tetra=cree_tetra(mgvol,noeudref1,noeudref2,noeudref3,noeudref4,triangle1,triangle2,triangle3,triangle4,IMPOSE);
1094     mg_maillage->ajouter_mg_tetra(tetra);
1095     triangle1->ajouter_voisin(tetra);
1096     triangle2->ajouter_voisin(tetra);
1097     triangle3->ajouter_voisin(tetra);
1098     triangle4->ajouter_voisin(tetra);
1099     MG_FRONT_3D *ftref1=NULL;
1100     MG_FRONT_3D *ftref2=NULL;
1101     MG_FRONT_3D *ftref3=NULL;
1102     MG_FRONT_3D *ftref4=NULL;
1103     if (triangle1->get_etat_front()==FRONT_ACTIF) ftref1=triangle1->get_mgfront();
1104     if (triangle2->get_etat_front()==FRONT_ACTIF) ftref2=triangle2->get_mgfront();
1105     if (triangle3->get_etat_front()==FRONT_ACTIF) ftref3=triangle3->get_mgfront();
1106     if (triangle4->get_etat_front()==FRONT_ACTIF) ftref4=triangle4->get_mgfront();
1107     TPL_MAP_ENTITE<MG_FRONT_3D*> lstfront;
1108     MG_FRONT_3D *ft1=NULL;
1109     MG_FRONT_3D *ft2=NULL;
1110     MG_FRONT_3D *ft3=NULL;
1111     MG_FRONT_3D *ft4=NULL;
1112     if (ftref1==NULL) ft1=mise_a_jour_front(tetra,0,NONFORCE);
1113     else
1114     {
1115     lstfront.ajouter(ftref1->get_front_voisin(0));
1116     lstfront.ajouter(ftref1->get_front_voisin(1));
1117     lstfront.ajouter(ftref1->get_front_voisin(2));
1118     supprimer_front_en_avancant(ftref1);
1119     }
1120     if (ftref2==NULL) ft2=mise_a_jour_front(tetra,1,NONFORCE);
1121     else
1122     {
1123     lstfront.ajouter(ftref2->get_front_voisin(0));
1124     lstfront.ajouter(ftref2->get_front_voisin(1));
1125     lstfront.ajouter(ftref2->get_front_voisin(2));
1126     supprimer_front_en_avancant(ftref2);
1127     }
1128     if (ftref3==NULL) ft3=mise_a_jour_front(tetra,2,NONFORCE);
1129     else
1130     {
1131     lstfront.ajouter(ftref3->get_front_voisin(0));
1132     lstfront.ajouter(ftref3->get_front_voisin(1));
1133     lstfront.ajouter(ftref3->get_front_voisin(2));
1134     supprimer_front_en_avancant(ftref3);
1135     }
1136     if (ftref4==NULL) ft4=mise_a_jour_front(tetra,3,NONFORCE);
1137     else
1138     {
1139     lstfront.ajouter(ftref4->get_front_voisin(0));
1140     lstfront.ajouter(ftref4->get_front_voisin(1));
1141     lstfront.ajouter(ftref4->get_front_voisin(2));
1142     supprimer_front_en_avancant(ftref4);
1143     }
1144     if (ftref1!=NULL) lstfront.supprimer(ftref1);
1145     if (ftref2!=NULL) lstfront.supprimer(ftref2);
1146     if (ftref3!=NULL) lstfront.supprimer(ftref3);
1147     if (ftref4!=NULL) lstfront.supprimer(ftref4);
1148     if (ft1!=NULL)
1149     {
1150     mise_a_jour_voisin_front(ft1);
1151     lstfront.ajouter(ft1->get_front_voisin(0));
1152     lstfront.ajouter(ft1->get_front_voisin(1));
1153     lstfront.ajouter(ft1->get_front_voisin(2));
1154     }
1155     if (ft2!=NULL)
1156     {
1157     mise_a_jour_voisin_front(ft2);
1158     lstfront.ajouter(ft2->get_front_voisin(0));
1159     lstfront.ajouter(ft2->get_front_voisin(1));
1160     lstfront.ajouter(ft2->get_front_voisin(2));
1161     }
1162     if (ft3!=NULL)
1163     {
1164     mise_a_jour_voisin_front(ft3);
1165     lstfront.ajouter(ft3->get_front_voisin(0));
1166     lstfront.ajouter(ft3->get_front_voisin(1));
1167     lstfront.ajouter(ft3->get_front_voisin(2));
1168     }
1169     if (ft4!=NULL)
1170     {
1171     mise_a_jour_voisin_front(ft4);
1172     lstfront.ajouter(ft4->get_front_voisin(0));
1173     lstfront.ajouter(ft4->get_front_voisin(1));
1174     lstfront.ajouter(ft4->get_front_voisin(2));
1175     }
1176    
1177     for (int l=0;l<lstfront.get_nb();l++) mise_a_jour_voisin_front(lstfront.get(l));
1178    
1179     }
1180     }
1181    
1182     void MAILLEUR3D::insere_contrainte_triangle(MG_VOLUME* mgvol,TPL_LISTE_ENTITE<MG_TRIANGLE*> *lsttri)
1183     {
1184     int nbtri=lsttri->get_nb();
1185     for (int i=0;i<nbtri;i++)
1186     {
1187     M3D_TRIANGLE* tri=(M3D_TRIANGLE*)lsttri->get(i);
1188     double *xyz1=tri->get_noeud1()->get_coord();
1189     double *xyz2=tri->get_noeud2()->get_coord();
1190     double *xyz3=tri->get_noeud3()->get_coord();
1191     double xyz[3];
1192     xyz[0]=(xyz1[0]+xyz2[0]+xyz3[0])/3.;
1193     xyz[1]=(xyz1[1]+xyz2[1]+xyz3[1])/3.;
1194     xyz[2]=(xyz1[2]+xyz2[2]+xyz3[2])/3.;
1195     TPL_MAP_ENTITE<MG_FRONT_3D*> liste_trouvee;
1196     OT_VECTEUR_3D vec1(xyz,xyz1);
1197     OT_VECTEUR_3D vec2(xyz,xyz2);
1198     OT_VECTEUR_3D vec3(xyz,xyz3);
1199     double rayonrecherche=vec1.get_longueur();
1200     if (vec2.get_longueur()>rayonrecherche) rayonrecherche=vec2.get_longueur();
1201     if (vec3.get_longueur()>rayonrecherche) rayonrecherche=vec3.get_longueur();
1202     rayonrecherche=rayonrecherche*1.1;
1203     octree_de_front->rechercher(xyz[0],xyz[2],xyz[2],rayonrecherche,liste_trouvee);
1204     MG_NOEUD* noeudref1=NULL;
1205     MG_NOEUD* noeudref2=NULL;
1206     MG_NOEUD* noeudref3=NULL;
1207     int nbtrouve=liste_trouvee.get_nb();
1208     for (int j=0;j<nbtrouve;j++)
1209     {
1210     MG_FRONT_3D* ft=liste_trouvee.get(j);
1211     TPL_MAP_ENTITE<MG_NOEUD*> lstnoeud;
1212     lstnoeud.ajouter(ft->get_noeud1());
1213     lstnoeud.ajouter(ft->get_noeud2());
1214     lstnoeud.ajouter(ft->get_noeud3());
1215     int nb=lstnoeud.get_nb();
1216     for (int k=0;k<nb;k++)
1217     {
1218     double *xyztmp=lstnoeud.get(k)->get_coord();
1219     OT_VECTEUR_3D vec(xyztmp,xyz1);
1220     double longueur=vec.get_longueur();
1221     if (longueur<1e-10*rayonrecherche)
1222     {
1223     noeudref1=lstnoeud.get(k);
1224     lstnoeud.supprimer(lstnoeud.get(k));
1225     break;
1226     }
1227     }
1228     nb=lstnoeud.get_nb();
1229     for (int k=0;k<nb;k++)
1230     {
1231     double *xyztmp=lstnoeud.get(k)->get_coord();
1232     OT_VECTEUR_3D vec(xyztmp,xyz2);
1233     double longueur=vec.get_longueur();
1234     if (longueur<1e-10*rayonrecherche)
1235     {
1236     noeudref2=lstnoeud.get(k);
1237     lstnoeud.supprimer(lstnoeud.get(k));
1238     break;
1239     }
1240     }
1241     nb=lstnoeud.get_nb();
1242     for (int k=0;k<nb;k++)
1243     {
1244     double *xyztmp=lstnoeud.get(k)->get_coord();
1245     OT_VECTEUR_3D vec(xyztmp,xyz3);
1246     double longueur=vec.get_longueur();
1247     if (longueur<1e-10*rayonrecherche)
1248     {
1249     noeudref3=lstnoeud.get(k);
1250     lstnoeud.supprimer(lstnoeud.get(k));
1251     break;
1252     }
1253     }
1254    
1255     }
1256     if (noeudref1==NULL) noeudref1=mg_maillage->ajouter_mg_noeud(mgvol,xyz1[0],xyz1[1],xyz1[2],IMPOSE);
1257     if (noeudref2==NULL) noeudref2=mg_maillage->ajouter_mg_noeud(mgvol,xyz2[0],xyz2[1],xyz2[2],IMPOSE);
1258     if (noeudref3==NULL) noeudref3=mg_maillage->ajouter_mg_noeud(mgvol,xyz3[0],xyz3[1],xyz3[2],IMPOSE);
1259     MG_SEGMENT* segment1=mg_maillage->get_mg_segment(noeudref1->get_id(),noeudref2->get_id());
1260     MG_SEGMENT* segment2=mg_maillage->get_mg_segment(noeudref1->get_id(),noeudref3->get_id());
1261     MG_SEGMENT* segment3=mg_maillage->get_mg_segment(noeudref2->get_id(),noeudref3->get_id());
1262     if (segment1==NULL) segment1=cree_segment(mgvol,noeudref1,noeudref2,IMPOSE);
1263     if (segment2==NULL) segment2=cree_segment(mgvol,noeudref1,noeudref3,IMPOSE);
1264     if (segment3==NULL) segment3=cree_segment(mgvol,noeudref2,noeudref3,IMPOSE);
1265     MG_TRIANGLE* nvtri=cree_triangle(mgvol,noeudref1,noeudref2,noeudref3,segment1,segment2,segment3,IMPOSE);
1266     MG_TRIANGLE* nvtri2=cree_triangle(mgvol,noeudref1,noeudref3,noeudref2,segment1,segment2,segment3,IMPOSE);
1267     //MG_TRIANGLE* nvtri=insere_triangle(mgvol,noeudref1,noeudref2,noeudref3,IMPOSE);
1268     //MG_TRIANGLE* nvtri2=insere_triangle(mgvol,noeudref1,noeudref3,noeudref2,IMPOSE);
1269     MG_FRONT_3D *ft=ajouter_front_courant(NONFORCE,noeudref1,noeudref2,noeudref3,nvtri);
1270     MG_FRONT_3D *ft2=ajouter_front_courant(NONFORCE,noeudref1,noeudref3,noeudref2,nvtri2);
1271     mise_a_jour_voisin_front(ft);
1272     mise_a_jour_voisin_front(ft2);
1273     mise_a_jour_voisin_front(ft->get_front_voisin(0));
1274     mise_a_jour_voisin_front(ft->get_front_voisin(1));
1275     mise_a_jour_voisin_front(ft->get_front_voisin(2));
1276     mise_a_jour_voisin_front(ft2->get_front_voisin(0));
1277     mise_a_jour_voisin_front(ft2->get_front_voisin(1));
1278     mise_a_jour_voisin_front(ft2->get_front_voisin(2));
1279    
1280     }
1281    
1282    
1283    
1284    
1285     }
1286