ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/mailleur3d_outil.cpp
Revision: 419
Committed: Tue Sep 10 20:10:09 2013 UTC (11 years, 8 months ago) by francois
File size: 51167 byte(s)
Log Message:
Suppression d'opération superflue dans le mailleur 3D frontal

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