ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur_auto/src/mailleur3d.cpp
Revision: 78
Committed: Thu Apr 10 00:05:47 2008 UTC (17 years, 1 month ago) by francois
Original Path: magic/lib/mailleur/mailleur/src/mailleur3d.cpp
File size: 22569 byte(s)
Log Message:
Ajout du mailleur octal pour XFEM et ajout de fusion de face pour obtenir les cofaces avec vectorisation ameliore

File Contents

# User Rev Content
1 5 //------------------------------------------------------------
2     //------------------------------------------------------------
3     // MAGiC
4     // Jean Christophe Cuilli�e et Vincent FRANCOIS
5     // D�artement de G�ie M�anique - UQTR
6     //------------------------------------------------------------
7     // Le projet MAGIC est un projet de recherche du d�artement
8     // de g�ie m�anique de l'Universit�du Qu�ec �
9     // Trois Rivi�es
10     // Les librairies ne peuvent �re utilis�s sans l'accord
11     // des auteurs (contact : francois@uqtr.ca)
12     //------------------------------------------------------------
13     //------------------------------------------------------------
14     //
15     // mailleur3d.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     #include <time.h>
30     #include "tpl_set.h"
31     #include "mailleur3d.h"
32     #include "m3d_triangle.h"
33 francois 35 #include "mg_gestionnaire.h"
34 5 //---------------------------------------------------------------------------
35    
36     //#pragma package(smart_init)
37    
38    
39    
40 francois 78 MAILLEUR3D::MAILLEUR3D(MG_MAILLAGE* mgmai,MG_GEOMETRIE *mggeo,MG_VOLUME* mgvol,FCT_TAILLE* fct_taille):MAILLEUR(),mg_maillage(mgmai),mg_geometrie(mggeo),mg_volume(mgvol),metrique(fct_taille),activelog(0),niveau_optimisation(2),type_mailleur(FRONTAL)
41 5 {
42    
43     for (int i=0;i<100;i++)
44     {
45     cpu[i]=0.;
46     nbfois[i]=0;
47     }
48     }
49    
50    
51 francois 78 MAILLEUR3D::MAILLEUR3D(FCT_GENERATEUR_3D<4> *carte,MG_MAILLAGE* mgmai):metrique(carte),type_mailleur(OCTAL),mg_maillage(mgmai),coef(2.)
52     {
53     }
54 5
55 francois 78
56 5 MAILLEUR3D::~MAILLEUR3D()
57     {
58     }
59    
60    
61    
62    
63     void MAILLEUR3D::maille(void)
64     {
65 francois 78 if (type_mailleur==FRONTAL)
66 5 {
67 francois 78 if (mg_volume!=NULL) maille(mg_volume);
68     else
69 5 {
70 francois 78 int nb_vol=mg_geometrie->get_nb_mg_volume();
71     for (int i=0;i<nb_vol;i++)
72     {
73     MG_VOLUME* mgvol=mg_geometrie->get_mg_volume(i);
74     maille(mgvol);
75     }
76 5 }
77     }
78 francois 78 if (type_mailleur==OCTAL)
79     {
80     FCT_GENERATEUR_3D<4> *carte=(FCT_GENERATEUR_3D<4> *)metrique;
81     maille(carte);
82     }
83 5
84     }
85    
86    
87    
88 francois 54 void MAILLEUR3D::maille(MG_VOLUME* mgvol,TPL_LISTE_ENTITE<MG_TRIANGLE*> *lsttri,TPL_LISTE_ENTITE<MG_TETRA*> *lsttet)
89 5 {
90     initialise_frontiere(mgvol);
91     cree_octree();
92     int creation_metrique=0;
93     if (metrique==NULL)
94     {
95     metrique=new FCT_TAILLE_VOLUME(distance_maximale,octree_de_frontiere);
96     creation_metrique=1;
97     }
98 francois 35 else distance_maximale=metrique->get_valeur_maximale(0);
99 francois 54 if (affichageactif==1) (*affiche)(" Construction du front");
100 francois 35 initialise_front(mgvol);
101 francois 54 if (lsttri!=NULL) insere_contrainte_triangle(mgvol,lsttri);
102     if (lsttet!=NULL) insere_contrainte_tetra(mgvol,lsttet);
103     if (affichageactif==1) (*affiche)(" Progression du front");
104 francois 61 //mg_maillage->get_gestionnaire()->enregistrer("c:\\void.magic");
105 5 progresse_front(mgvol);
106 francois 54 if (affichageactif==1) (*affiche)(" Construction du maillage final");
107 5 detruit_element_inactif();
108 francois 54 if (affichageactif==1) (*affiche)(" Optimisation du maillage final");
109 francois 61 int nbaoptimiser;
110     int nbaoptimiserapres=mg_maillage->get_nb_mg_tetra();
111     do
112     {
113     nbaoptimiser=nbaoptimiserapres;
114     optimise(mgvol,nbaoptimiserapres);
115     }
116     while (nbaoptimiserapres!=nbaoptimiser);
117 5 if (creation_metrique==1) delete metrique;
118    
119     delete octree_de_triangle;
120     delete octree_de_frontiere;
121     delete octree_de_front;
122    
123     triangle_frontiere.vide();
124     noeud_frontiere.vide();
125     noeud_arete_frontiere.vide();
126     }
127    
128     void MAILLEUR3D::initialise_frontiere(MG_VOLUME* mgvol)
129     {
130     clock_t cpu1=clock(); //debug
131     int nbcoq=mgvol->get_nb_mg_coquille();
132     for (int i=0;i<nbcoq;i++)
133     {
134     MG_COQUILLE* mgcoq=mgvol->get_mg_coquille(i);
135     int nbcoface=mgcoq->get_nb_mg_coface();
136     for (int j=0;j<nbcoface;j++)
137     {
138     MG_COFACE* coface=mgcoq->get_mg_coface(j);
139     MG_FACE* mgface=coface->get_face();
140     TPL_SET<MG_ELEMENT_MAILLAGE*> * lien_maillage = mgface->get_lien_maillage();
141     TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR it;
142     MG_ELEMENT_MAILLAGE* element;
143     for (element = lien_maillage->get_premier(it); element; element = lien_maillage->get_suivant(it) )
144     {
145     MG_TRIANGLE* mgtri=(MG_TRIANGLE*)element;
146     MG_TRIANGLE* tritemp=mg_maillage->get_mg_triangleid(mgtri->get_id());
147     if (tritemp==NULL) continue;
148     triangle_frontiere.ajouter(mgtri);
149     noeud_frontiere.ajouter(mgtri->get_noeud1());
150     noeud_frontiere.ajouter(mgtri->get_noeud2());
151     noeud_frontiere.ajouter(mgtri->get_noeud3());
152     }
153     int nbbou=mgface->get_nb_mg_boucle();
154     for (int l=0;l<nbbou;l++)
155     {
156     MG_BOUCLE* mgbou=mgface->get_mg_boucle(l);
157     int nbcoarete=mgbou->get_nb_mg_coarete();
158     for (int m=0;m<nbcoarete;m++)
159     {
160     MG_COARETE* mgcoarete=mgbou->get_mg_coarete(m);
161     MG_ARETE* mgarete=mgcoarete->get_arete();
162     TPL_SET<MG_ELEMENT_MAILLAGE*> * lien_maillage = mgarete->get_lien_maillage();
163     TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR it;
164     MG_ELEMENT_MAILLAGE* element;
165     for (element = lien_maillage->get_premier(it); element; element = lien_maillage->get_suivant(it) )
166     {
167     MG_SEGMENT* mgseg=(MG_SEGMENT*)element;
168     MG_SEGMENT* segtemp=mg_maillage->get_mg_segmentid(mgseg->get_id());
169     if (segtemp==NULL) continue;
170     noeud_arete_frontiere.ajouter(mgseg->get_noeud1());
171     noeud_arete_frontiere.ajouter(mgseg->get_noeud2());
172     }
173     }
174     }
175     }
176     }
177    
178     clock_t cpu2=clock(); //debug
179     cpu[0]=cpu[0]+(cpu2-cpu1)/CLK_TCK; //debug
180     nbfois[0]++; //debug
181     }
182    
183     void MAILLEUR3D::cree_octree(void)
184     {
185     double xmin,ymin,zmin;
186     double xmax,ymax,zmax;
187    
188     clock_t cpu1=clock(); // debug
189     int nb_noeud=noeud_arete_frontiere.get_nb();
190     for (int i=0;i<nb_noeud;i++)
191     {
192     MG_NOEUD* mgnoeud=noeud_arete_frontiere.get(i);
193     if (i==0)
194     {
195     xmin=mgnoeud->get_x();
196     xmax=mgnoeud->get_x();
197     ymin=mgnoeud->get_y();
198     ymax=mgnoeud->get_y();
199     zmin=mgnoeud->get_z();
200     zmax=mgnoeud->get_z();
201     }
202     else
203     {
204     xmin=std::min(xmin,mgnoeud->get_x());
205     ymin=std::min(ymin,mgnoeud->get_y());
206     zmin=std::min(zmin,mgnoeud->get_z());
207     xmax=std::max(xmax,mgnoeud->get_x());
208     ymax=std::max(ymax,mgnoeud->get_y());
209     zmax=std::max(zmax,mgnoeud->get_z());
210     }
211     double xdiff=xmax-xmin;
212     double ydiff=ymax-ymin;
213     double zdiff=zmax-zmin;
214     xmin=xmin-0.125*xdiff;
215     xmax=xmax+0.125*xdiff;
216     ymin=ymin-0.125*ydiff;
217     ymax=ymax+0.125*ydiff;
218     zmin=zmin-0.125*zdiff;
219     zmax=zmax+0.125*zdiff;
220    
221     }
222     octree_de_triangle=new TPL_OCTREE<MG_TRIANGLE*,MG_NOEUD*>;
223     octree_de_frontiere=new TPL_OCTREE<MG_TRIANGLE*,MG_NOEUD*>;
224     octree_de_front=new TPL_OCTREE<MG_FRONT_3D*,MG_NOEUD*>;
225     octree_de_triangle->initialiser(&noeud_arete_frontiere,1,xmin,ymin,zmin,xmax,ymax,zmax);
226     octree_de_frontiere->initialiser(octree_de_triangle);
227     octree_de_front->initialiser(octree_de_triangle);
228    
229     int nb_triangle=triangle_frontiere.get_nb();
230     for (int j=0;j<nb_triangle;j++)
231     {
232     octree_de_frontiere->inserer(triangle_frontiere.get(j));
233     octree_de_triangle->inserer(triangle_frontiere.get(j));
234     }
235     clock_t cpu2=clock(); //debug
236     cpu[1]=cpu[1]+(cpu2-cpu1)/CLK_TCK; //debug
237     nbfois[1]++; //debug
238     }
239    
240     void MAILLEUR3D::initialise_front(MG_VOLUME* mgvol)
241     {
242     clock_t cpu1=clock(); //debug
243     // constitution du front
244     int nbcoq=mgvol->get_nb_mg_coquille();
245     for (int i=0;i<nbcoq;i++)
246     {
247     MG_COQUILLE* mgcoq=mgvol->get_mg_coquille(i);
248     int nbcoface=mgcoq->get_nb_mg_coface();
249     for (int j=0;j<nbcoface;j++)
250     {
251     MG_COFACE* coface=mgcoq->get_mg_coface(j);
252     MG_FACE* mgface=coface->get_face();
253    
254     TPL_SET<MG_ELEMENT_MAILLAGE*> * lien_maillage = mgface->get_lien_maillage();
255     TPL_SET<MG_ELEMENT_MAILLAGE*>::ITERATEUR it;
256     MG_ELEMENT_MAILLAGE* element;
257     for (element = lien_maillage->get_premier(it); element; element = lien_maillage->get_suivant(it) )
258     {
259     MG_TRIANGLE* mgtri=(MG_TRIANGLE*)element;
260     MG_TRIANGLE* tritemp=mg_maillage->get_mg_triangleid(mgtri->get_id());
261     if (tritemp==NULL) continue;
262     if (coface->get_orientation()==1) ajouter_front_courant(NONFORCE,mgtri->get_noeud1(),mgtri->get_noeud2(),mgtri->get_noeud3(),mgtri);
263     else ajouter_front_courant(NONFORCE,mgtri->get_noeud1(),mgtri->get_noeud3(),mgtri->get_noeud2(),mgtri);
264     }
265     }
266     }
267    
268     // recherche des voisins
269     for (int l=0;l<11;l++)
270     {
271     int nb_front=get_nb_front(front_courant[l]);
272     for (int i=0;i<nb_front;i++)
273     {
274     MG_FRONT_3D* ft=get_front(front_courant[l],i);
275     MG_NOEUD* mgnoeud1=ft->get_noeud1();
276     MG_NOEUD* mgnoeud2=ft->get_noeud2();
277     MG_NOEUD* mgnoeud3=ft->get_noeud3();
278     ft->changer_front_voisin(NULL,0);
279     ft->changer_front_voisin(NULL,1);
280     ft->changer_front_voisin(NULL,2);
281     for (int j=0;j<mgnoeud1->get_lien_triangle()->get_nb();j++)
282     for (int k=0;k<mgnoeud2->get_lien_triangle()->get_nb();k++)
283     {
284     M3D_TRIANGLE* tri1=(M3D_TRIANGLE*)mgnoeud1->get_lien_triangle()->get(j);
285     M3D_TRIANGLE* tri2=(M3D_TRIANGLE*)mgnoeud2->get_lien_triangle()->get(k);
286     if (tri1==tri2)
287     if (tri1->get_etat_front()==FRONT_ACTIF)
288     if (tri1->get_mgfront()!=ft)
289     {
290     MG_FRONT_3D* ft2=tri1->get_mgfront();
291     if (ft->get_front_voisin(0)==NULL)
292     {
293     ft->changer_front_voisin(ft2,0);
294     ft->changer_angle_voisin(angle_front(ft,ft2),0);
295     }
296     else
297     {
298     double angle=angle_front(ft,ft2);
299     if (angle<ft->get_angle_voisin(0))
300     {
301     ft->changer_front_voisin(ft2,0);
302     ft->changer_angle_voisin(angle,0);
303     }
304     }
305     }
306     }
307     for (int j=0;j<mgnoeud1->get_lien_triangle()->get_nb();j++)
308     for (int k=0;k<mgnoeud3->get_lien_triangle()->get_nb();k++)
309     {
310     M3D_TRIANGLE* tri1=(M3D_TRIANGLE*)mgnoeud1->get_lien_triangle()->get(j);
311     M3D_TRIANGLE* tri2=(M3D_TRIANGLE*)mgnoeud3->get_lien_triangle()->get(k);
312     if (tri1==tri2)
313     if (tri1->get_etat_front()==FRONT_ACTIF)
314     if (tri1->get_mgfront()!=ft)
315     {
316     MG_FRONT_3D* ft2=tri1->get_mgfront();
317     if (ft->get_front_voisin(1)==NULL)
318     {
319     ft->changer_front_voisin(ft2,1);
320     ft->changer_angle_voisin(angle_front(ft,ft2),1);
321     }
322     else
323     {
324     double angle=angle_front(ft,ft2);
325     if (angle<ft->get_angle_voisin(1))
326     {
327     ft->changer_front_voisin(ft2,1);
328     ft->changer_angle_voisin(angle,1);
329     }
330     }
331     }
332     }
333     for (int j=0;j<mgnoeud2->get_lien_triangle()->get_nb();j++)
334     for (int k=0;k<mgnoeud3->get_lien_triangle()->get_nb();k++)
335     {
336     M3D_TRIANGLE* tri1=(M3D_TRIANGLE*)mgnoeud2->get_lien_triangle()->get(j);
337     M3D_TRIANGLE* tri2=(M3D_TRIANGLE*)mgnoeud3->get_lien_triangle()->get(k);
338     if (tri1==tri2)
339     if (tri1->get_etat_front()==FRONT_ACTIF)
340     if (tri1->get_mgfront()!=ft)
341     {
342     MG_FRONT_3D* ft2=tri1->get_mgfront();
343     if (ft->get_front_voisin(2)==NULL)
344     {
345     ft->changer_front_voisin(ft2,2);
346     ft->changer_angle_voisin(angle_front(ft,ft2),2);
347     }
348     else
349     {
350     double angle=angle_front(ft,ft2);
351     if (angle<ft->get_angle_voisin(2))
352     {
353     ft->changer_front_voisin(ft2,2);
354     ft->changer_angle_voisin(angle,2);
355     }
356     }
357     }
358     }
359     if ( (ft->get_front_voisin(0)==NULL) || (ft->get_front_voisin(1)==NULL)|| (ft->get_front_voisin(2)==NULL))
360     refresh();
361     }
362    
363     }
364     clock_t cpu2=clock(); //debug
365     cpu[2]=cpu[2]+(cpu2-cpu1)/CLK_TCK; //debug
366     nbfois[2]++; //debug
367     }
368    
369    
370    
371     // FONCTIONS GERANT L ENTITE FRONT (ajout suppression et tout le bordel)
372     MG_FRONT_3D* MAILLEUR3D::ajouter_front_courant(int numero_front,MG_NOEUD* noeud1,MG_NOEUD* noeud2,MG_NOEUD* noeud3,MG_TRIANGLE* triangle)
373     {
374     int numfront;
375     if (numero_front==NONFORCE)
376     {
377     numfront=(int)(triangle->get_longueur()/distance_maximale*10.);
378     if (numfront>10) numfront=10;
379     }
380     else numfront=numero_front;
381     MG_FRONT_3D* mgfront;
382     if (numfront>11) mgfront=ajouter_front(front_attente,noeud1,noeud2,noeud3,triangle);
383     else mgfront=ajouter_front(front_courant[numfront],noeud1,noeud2,noeud3,triangle);
384     mgfront->changer_numero_front(numfront);
385     if (triangle->get_type_entite()==IDM3D_TRIANGLE)
386     {
387     M3D_TRIANGLE* mtri=(M3D_TRIANGLE*)triangle;
388     mtri->change_etat_front(FRONT_ACTIF);
389     mtri->change_mgfront(mgfront);
390     }
391     return mgfront;
392     }
393    
394     MG_FRONT_3D* MAILLEUR3D::ajouter_front(FRONT& front,MG_NOEUD* noeud1,MG_NOEUD* noeud2,MG_NOEUD* noeud3,MG_TRIANGLE* triangle)
395     {
396     MG_FRONT_3D* ft=new MG_FRONT_3D(noeud1,noeud2,noeud3,triangle);
397     ajouter_front(front,ft);
398     octree_de_front->inserer(ft);
399     return ft;
400     }
401    
402     void MAILLEUR3D::ajouter_front(FRONT& front,MG_FRONT_3D *ft)
403     {
404     std::pair<const double,MG_FRONT_3D*> tmp(ft->get_triangle()->get_longueur(),ft);
405     front.insert(tmp);
406     }
407    
408     MG_FRONT_3D* MAILLEUR3D::get_front(FRONT& front,unsigned int num)
409     {
410     FRONT::iterator i=front.begin();
411     for (unsigned long j=0;j<num;j++) i++;
412     return ((*i).second);
413     }
414    
415    
416     unsigned int MAILLEUR3D::get_nb_front(FRONT& front)
417     {
418     return front.size();
419     }
420    
421     unsigned int MAILLEUR3D::get_nb_front_courant(void)
422     {
423    
424     return front_courant[0].size()+front_courant[1].size()+front_courant[2].size()+front_courant[3].size()+front_courant[4].size()+front_courant[5].size()+front_courant[6].size()+front_courant[7].size()+front_courant[8].size()+front_courant[9].size()+front_courant[10].size();
425     }
426    
427     void MAILLEUR3D::supprimer_front_en_avancant(MG_FRONT_3D* ft)
428     {
429     octree_de_front->supprimer(ft);
430     int numfront=ft->get_numero_front();
431     FRONT* front_original;
432     if (numfront==ATTENTE)
433     front_original=&front_attente;
434     else
435     front_original=front_courant+numfront;
436     FRONT::iterator j=front_original->lower_bound(ft->get_triangle()->get_longueur());
437     while (ft!=(*j).second) j++;
438     front_original->erase(j);
439     MG_TRIANGLE* mgtri=ft->get_triangle();
440     if (mgtri->get_type_entite()==IDM3D_TRIANGLE)
441     {
442     M3D_TRIANGLE* m3dtri=(M3D_TRIANGLE*)mgtri;
443     m3dtri->change_etat_front(FRONT_PASSE);
444     m3dtri->change_mgfront(NULL);
445     }
446     delete ft;
447     }
448    
449     void MAILLEUR3D::supprimer_front_en_reculant(MG_FRONT_3D* ft)
450     {
451     octree_de_front->supprimer(ft);
452     int numfront=ft->get_numero_front();
453     FRONT* front_original;
454     if (numfront==ATTENTE)
455     front_original=&front_attente;
456     else
457     front_original=front_courant+numfront;
458     FRONT::iterator j=front_original->lower_bound(ft->get_triangle()->get_longueur());
459     while (ft!=(*j).second) j++;
460     front_original->erase(j);
461     MG_TRIANGLE* mgtri=ft->get_triangle();
462     if (mgtri->get_type_entite()==IDM3D_TRIANGLE)
463     {
464     M3D_TRIANGLE* m3dtri=(M3D_TRIANGLE*)mgtri;
465     m3dtri->change_etat_front(FRONT_NONACTIF);
466     m3dtri->change_mgfront(NULL);
467     }
468     delete ft;
469     }
470    
471     void MAILLEUR3D::echange_de_front(FRONT& front_original,FRONT& front_destination,MG_FRONT_3D* ft)
472     {
473     FRONT::iterator j=front_original.lower_bound(ft->get_triangle()->get_longueur());
474     while (ft!=(*j).second) j++;
475     front_original.erase(j);
476     ajouter_front(front_destination,ft);
477     }
478    
479     void MAILLEUR3D::echange_de_front(MG_FRONT_3D* ft,int num)
480     {
481     int numfront;
482     if (num==NONFORCE)
483     {
484     numfront=(int)(ft->get_triangle()->get_longueur()/distance_maximale*10.);
485     if (numfront>10) numfront=10;
486     }
487     else if (num==ATTENTE) numfront=ATTENTE;
488     else numfront=num;
489     if (numfront>10) numfront=ATTENTE;
490     FRONT* front_original;
491     int num_front_original=ft->get_numero_front();
492     if (num_front_original==ATTENTE)
493     front_original=&front_attente;
494     else
495     front_original=front_courant+num_front_original;
496     FRONT* front_destination;
497     if (numfront==ATTENTE)
498     front_destination=&front_attente;
499     else
500     front_destination=front_courant+numfront;
501    
502    
503    
504     FRONT::iterator j=front_original->lower_bound(ft->get_triangle()->get_longueur());
505     while (ft!=(*j).second) j++;
506     front_original->erase(j);
507     ft->changer_numero_front(numfront);
508     ajouter_front(*front_destination,ft);
509    
510     }
511    
512    
513     void MAILLEUR3D::swap_front(FRONT& front_original,FRONT& front_destination)
514     {
515     front_original.swap(front_destination);
516     }
517    
518    
519    
520     int MAILLEUR3D::premier_front_courant(void)
521     {
522     int num=-1;
523     if (front_courant[0].size()>0) num=0;
524     else if (front_courant[1].size()>0) num=1;
525     else if (front_courant[2].size()>0) num=2;
526     else if (front_courant[3].size()>0) num=3;
527     else if (front_courant[4].size()>0) num=4;
528     else if (front_courant[5].size()>0) num=5;
529     else if (front_courant[6].size()>0) num=6;
530     else if (front_courant[7].size()>0) num=7;
531     else if (front_courant[8].size()>0) num=8;
532     else if (front_courant[9].size()>0) num=9;
533     else if (front_courant[10].size()>0) num=10;
534     return num;
535     }
536    
537    
538    
539     void MAILLEUR3D::active_log(char* nomfile)
540     {
541     activelog=1;
542     nomlog=nomfile;
543     }
544    
545     void MAILLEUR3D::desactive_log(void)
546     {
547     activelog=0;
548     }
549    
550 francois 35
551 francois 54