ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/REPOS_ERICCA/magic/lib/mailleur/src/mailleur3d_optimisation.cpp
Revision: 35
Committed: Thu Aug 16 21:24:14 2007 UTC (17 years, 9 months ago) by francois
Original Path: magic/lib/mailleur/mailleur/src/mailleur3d_optimisation.cpp
File size: 21017 byte(s)
Log Message:
Mailleur 3D + divers bugs + origine des elements de maillages

File Contents

# User Rev Content
1 francois 35 //------------------------------------------------------------
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     #include "mailleur3d.h"
27     #include "ot_mathematique.h"
28     #include "m3d_tetra.h"
29     #include "m3d_noeud.h"
30     //#include "mg_gestionnaire.h"
31     #include <math.h>
32    
33     //---------------------------------------------------------------------------
34    
35     //#pragma package(smart_init)
36    
37     double get_volume(MG_TETRA* tet)
38     {
39     double *xyz1=tet->get_noeud1()->get_coord();
40     double *xyz2=tet->get_noeud2()->get_coord();
41     double *xyz3=tet->get_noeud3()->get_coord();
42     double *xyz4=tet->get_noeud4()->get_coord();
43     OT_VECTEUR_3D ab(xyz1,xyz2);
44     OT_VECTEUR_3D ac(xyz1,xyz3);
45     OT_VECTEUR_3D ad(xyz1,xyz4);
46     double vol=(ab&ac)*ad;
47     vol=vol/6.;
48     return vol;
49     }
50    
51     double get_volume(double *xyz1,double *xyz2,double *xyz3,double *xyz4)
52     {
53     OT_VECTEUR_3D ab(xyz1,xyz2);
54     OT_VECTEUR_3D ac(xyz1,xyz3);
55     OT_VECTEUR_3D ad(xyz1,xyz4);
56     double vol=(ab&ac)*ad;
57     vol=vol/6.;
58     return vol;
59     }
60    
61     void MAILLEUR3D::optimise(MG_VOLUME* mgvol)
62     {
63     //preparation
64     LISTE_MG_TETRA::iterator it;
65     for (MG_TETRA* tet=mg_maillage->get_premier_tetra(it);tet;tet=mg_maillage->get_suivant_tetra(it))
66     {
67     M3D_TETRA* mtet=(M3D_TETRA*)tet;
68     double qual=OPERATEUR::qualite_tetra(mtet->get_noeud1()->get_coord(),mtet->get_noeud2()->get_coord(),mtet->get_noeud3()->get_coord(),mtet->get_noeud4()->get_coord());
69     mtet->change_qualite(qual);
70     ajouter_ordre_tetra(mtet,0);
71     }
72     o3d_data();
73     o3d_data2();
74     //on commence
75    
76     for (int phase=0;phase<2;phase++)
77     {
78     char mess[255];
79     sprintf(mess," Phase %d",phase+1);
80     affiche(mess);
81     while (lst_tetra[phase].size()>0)
82     {
83     ORDRE_TETRA::iterator i=lst_tetra[phase].begin();
84     M3D_TETRA* tet=(*i).second;
85     supprimer_ordre_tetra(tet);
86     double crit_avant=tet->get_qualite();
87     if (crit_avant<0.5)
88     {
89     double crit[4]={-1.,-1.,-1.,-1.};
90     M3D_NOEUD* no[4];
91     no[0]=(M3D_NOEUD*)tet->get_noeud1();
92     no[1]=(M3D_NOEUD*)tet->get_noeud2();
93     no[2]=(M3D_NOEUD*)tet->get_noeud3();
94     no[3]=(M3D_NOEUD*)tet->get_noeud4();
95     double x[4],y[4],z[4];
96     bouge_point(mgvol,no[0],crit[0],x[0],y[0],z[0]);
97     bouge_point(mgvol,no[1],crit[1],x[1],y[1],z[1]);
98     bouge_point(mgvol,no[2],crit[2],x[2],y[2],z[2]);
99     bouge_point(mgvol,no[3],crit[3],x[3],y[3],z[3]);
100     std::vector<double> vecteur_crit(crit,crit+4);
101     std::vector<double>::iterator it1 = max_element(vecteur_crit.begin(), vecteur_crit.end());
102     double crit_opt=*it1;
103     if (crit_opt>crit_avant)
104     {
105     int num=it1-vecteur_crit.begin();
106     no[num]->change_x(x[num]);
107     no[num]->change_y(y[num]);
108     no[num]->change_z(z[num]);
109     int nb=no[num]->get_lien_tetra()->get_nb();
110     for (int j=0;j<nb;j++)
111     {
112     M3D_TETRA* mtet=(M3D_TETRA*)no[num]->get_lien_tetra()->get(j);
113     double qual=OPERATEUR::qualite_tetra(mtet->get_noeud1()->get_coord(),mtet->get_noeud2()->get_coord(),mtet->get_noeud3()->get_coord(),mtet->get_noeud4()->get_coord());
114     mtet->change_qualite(qual);
115     }
116     crit_avant=crit_opt;
117     }
118     }
119     int coquille_non_optimise=1;
120     if (crit_avant<0.5)
121     {
122     COQUILLE coque[6];
123     double crit[6]={-1,-1,-1,-1,-1,-1};
124     remaille_coquille(tet->get_noeud1(),tet->get_noeud2(),crit[0],coque[0]);
125     remaille_coquille(tet->get_noeud1(),tet->get_noeud3(),crit[1],coque[1]);
126     remaille_coquille(tet->get_noeud1(),tet->get_noeud4(),crit[2],coque[2]);
127     remaille_coquille(tet->get_noeud2(),tet->get_noeud3(),crit[3],coque[3]);
128     remaille_coquille(tet->get_noeud2(),tet->get_noeud3(),crit[4],coque[4]);
129     remaille_coquille(tet->get_noeud3(),tet->get_noeud4(),crit[5],coque[5]);
130     std::vector<double> vecteur_crit(crit,crit+6);
131     std::vector<double>::iterator it1 = max_element(vecteur_crit.begin(), vecteur_crit.end());
132     double crit_opt=*it1;
133     if (crit_opt>crit_avant)
134     {
135     coquille_non_optimise=0;
136     int num=it1-vecteur_crit.begin();
137     for (int i=0;i<coque[num].taille;i++)
138     {
139     M3D_TETRA *tettmp=(M3D_TETRA*)coque[num].tet[i];
140     if (tet!=tettmp) supprimer_ordre_tetra(tettmp);
141     mg_maillage->supprimer_mg_tetraid(coque[num].tet[i]->get_id());
142     }
143     for (int i=0;i<2*coque[num].taille-4;i++)
144     {
145     MG_NOEUD* noeud1=coque[num].new_tetra[4*i];
146     MG_NOEUD* noeud2=coque[num].new_tetra[4*i+1];
147     MG_NOEUD* noeud3=coque[num].new_tetra[4*i+2];
148     MG_NOEUD* noeud4=coque[num].new_tetra[4*i+3];
149     MG_TRIANGLE* triangle1=mg_maillage->get_mg_triangle(noeud1->get_id(),noeud3->get_id(),noeud2->get_id());
150     MG_TRIANGLE* triangle2=mg_maillage->get_mg_triangle(noeud1->get_id(),noeud2->get_id(),noeud4->get_id());
151     MG_TRIANGLE* triangle3=mg_maillage->get_mg_triangle(noeud2->get_id(),noeud3->get_id(),noeud4->get_id());
152     MG_TRIANGLE* triangle4=mg_maillage->get_mg_triangle(noeud1->get_id(),noeud4->get_id(),noeud3->get_id());
153     if (triangle1==NULL) triangle1=insere_triangle(mgvol,noeud1,noeud3,noeud2);
154     if (triangle2==NULL) triangle2=insere_triangle(mgvol,noeud1,noeud2,noeud4);
155     if (triangle3==NULL) triangle3=insere_triangle(mgvol,noeud2,noeud3,noeud4);
156     if (triangle4==NULL) triangle4=insere_triangle(mgvol,noeud1,noeud4,noeud3);
157     M3D_TETRA* mtet=new M3D_TETRA(mgvol,noeud1,noeud2,noeud3,noeud4,triangle1,triangle2,triangle3,triangle4);
158     mg_maillage->ajouter_mg_tetra(mtet);
159     double qual=OPERATEUR::qualite_tetra(mtet->get_noeud1()->get_coord(),mtet->get_noeud2()->get_coord(),mtet->get_noeud3()->get_coord(),mtet->get_noeud4()->get_coord());
160     mtet->change_qualite(qual);
161     if (phase==0)
162     ajouter_ordre_tetra(mtet,1);
163     }
164    
165    
166     }
167    
168     }
169     if ((phase==0) && (coquille_non_optimise==1))
170     ajouter_ordre_tetra(tet,1);
171    
172     }
173    
174     }
175    
176    
177    
178     }
179    
180    
181    
182    
183    
184     int MAILLEUR3D::bouge_point(MG_VOLUME* mgvol,M3D_NOEUD* noeud,double& crit,double& x,double& y, double& z)
185     {
186     if (noeud->get_lien_topologie()!=mgvol) return 0;
187    
188     double xopt=0.;
189     double yopt=0.;
190     double zopt=0.;
191     double qual_dep=1.;
192     int nb_tet=noeud->get_lien_tetra()->get_nb();
193     double tab_coord[3000];
194     double gamma=0;
195     for (int i=0;i<nb_tet;i++)
196     {
197     M3D_TETRA* tet=(M3D_TETRA*)noeud->get_lien_tetra()->get(i);
198     qual_dep=std::min(qual_dep,tet->get_qualite());
199     MG_NOEUD* no1=tet->get_noeud1();
200     MG_NOEUD* no2=tet->get_noeud2();
201     MG_NOEUD* no3=tet->get_noeud3();
202     MG_NOEUD* no4=tet->get_noeud4();
203     MG_NOEUD *noeud1,*noeud2,*noeud3;
204     if (no1==noeud)
205     {
206     noeud1=no2;
207     noeud2=no4;
208     noeud3=no3;
209     }
210     if (no2==noeud)
211     {
212     noeud1=no1;
213     noeud2=no3;
214     noeud3=no4;
215     }
216     if (no3==noeud)
217     {
218     noeud1=no1;
219     noeud2=no4;
220     noeud3=no2;
221     }
222     if (no4==noeud)
223     {
224     noeud1=no1;
225     noeud2=no2;
226     noeud3=no3;
227     }
228     double xyzi[3];
229     double *xyz1=noeud1->get_coord();
230     double *xyz2=noeud2->get_coord();
231     double *xyz3=noeud3->get_coord();
232     xyzi[0]=(xyz1[0]+xyz2[0]+xyz3[0])/3.;
233     xyzi[1]=(xyz1[1]+xyz2[1]+xyz3[1])/3.;
234     xyzi[2]=(xyz1[2]+xyz2[2]+xyz3[2])/3.;
235     OT_VECTEUR_3D v12(xyz1,xyz2);
236     OT_VECTEUR_3D v13(xyz1,xyz3);
237     OT_VECTEUR_3D v23(xyz2,xyz3);
238     OT_VECTEUR_3D normal=v12&v13;
239     double perimetre=v12.get_longueur()+v13.get_longueur()+v23.get_longueur();
240     double hauteur = (perimetre/3.) * 0.8 ;
241     normal.norme();
242     tab_coord[3*i]=xyzi[0]+normal.get_x()*hauteur;
243     tab_coord[3*i+1]=xyzi[1]+normal.get_y()*hauteur;
244     tab_coord[3*i+2]=xyzi[2]+normal.get_z()*hauteur;
245     if (tet->get_qualite()> 1e-10) gamma=gamma+1./tet->get_qualite()/tet->get_qualite();
246     }
247     gamma=1./gamma;
248     for (int i=0;i<nb_tet;i++)
249     {
250     M3D_TETRA* tet=(M3D_TETRA*)noeud->get_lien_tetra()->get(i);
251     if (tet->get_qualite()> 1e-10)
252     {
253     double alpha=gamma/tet->get_qualite()/tet->get_qualite();
254     xopt=xopt+alpha*tab_coord[3*i];
255     yopt=yopt+alpha*tab_coord[3*i+1];
256     zopt=zopt+alpha*tab_coord[3*i+2];
257     }
258    
259     }
260     double delta[3]={xopt-noeud->get_x(),yopt-noeud->get_y(),zopt-noeud->get_z()};
261     double bmin=0.,bmax=1.;
262     double vieuxx=noeud->get_x();
263     double vieuxy=noeud->get_y();
264     double vieuxz=noeud->get_z();
265     double qualini=0.;
266     for (int iteration=0;iteration<5;iteration++)
267     {
268     double alpha=(bmin+bmax)*0.5;
269     noeud->change_x(vieuxx+alpha*delta[0]);
270     noeud->change_y(vieuxy+alpha*delta[1]);
271     noeud->change_z(vieuxz+alpha*delta[2]);
272     double qualcour=1.;
273     for (int i=0;i<nb_tet;i++)
274     {
275     M3D_TETRA* mtet=(M3D_TETRA*)noeud->get_lien_tetra()->get(i);
276     double qual=OPERATEUR::qualite_tetra(mtet->get_noeud1()->get_coord(),mtet->get_noeud2()->get_coord(),mtet->get_noeud3()->get_coord(),mtet->get_noeud4()->get_coord());
277     qualcour=std::min(qualcour,qual);
278     }
279     double alpha_eps=(bmin+bmax)*0.5-(bmax-bmin)/50.;
280     noeud->change_x(vieuxx+alpha_eps*delta[0]);
281     noeud->change_y(vieuxy+alpha_eps*delta[1]);
282     noeud->change_z(vieuxz+alpha_eps*delta[2]);
283     double qualcour_eps=1.;
284     for (int i=0;i<nb_tet;i++)
285     {
286     M3D_TETRA* mtet=(M3D_TETRA*)noeud->get_lien_tetra()->get(i);
287     double qual=OPERATEUR::qualite_tetra(mtet->get_noeud1()->get_coord(),mtet->get_noeud2()->get_coord(),mtet->get_noeud3()->get_coord(),mtet->get_noeud4()->get_coord());
288     qualcour_eps=std::min(qualcour_eps,qual);
289     }
290     if (qualcour>qualcour_eps) bmin =alpha; else bmax=alpha;
291     qualini=std::max(qualini,qualcour);
292     }
293     noeud->change_x(vieuxx);
294     noeud->change_y(vieuxy);
295     noeud->change_z(vieuxz);
296     if (qualini>qual_dep)
297     {
298     x=vieuxx+(bmin+bmax)*0.5*delta[0];
299     y=vieuxy+(bmin+bmax)*0.5*delta[1];
300     z=vieuxz+(bmin+bmax)*0.5*delta[2];
301     crit=qualini;
302     return 1;
303     }
304     return 0;
305     }
306    
307    
308    
309     void MAILLEUR3D::remaille_coquille(MG_NOEUD* noeud1,MG_NOEUD* noeud2, double& crit, COQUILLE& coque)
310     {
311     //recherche du segment a supprimer
312     MG_NOEUD* noeud;
313     if (noeud1->get_id()>noeud2->get_id()) noeud=noeud2; else noeud=noeud1;
314     int nb=noeud->get_lien_petit_segment()->get_nb();
315     for (int i=0;i<nb;i++)
316     {
317     MG_SEGMENT* seg=noeud->get_lien_petit_segment()->get(i);
318     if ( ((seg->get_noeud1()==noeud1) && (seg->get_noeud2()==noeud2)) || ((seg->get_noeud1()==noeud2) && (seg->get_noeud2()==noeud1)))
319     if (seg->get_lien_topologie()->get_dimension()!=3)
320     {
321     coque.taille=0;
322     crit=0.;
323     return;
324     }
325     }
326     //recherche des tetra qui s enroule autour du segment
327     int nb1=noeud1->get_lien_tetra()->get_nb();
328     int nb2=noeud2->get_lien_tetra()->get_nb();
329     MG_TETRA* coq[100];
330     int nb_coq=0;
331     for (int i=0;i<nb1;i++)
332     for (int j=0;j<nb2;j++)
333     if (noeud1->get_lien_tetra()->get(i)==noeud2->get_lien_tetra()->get(j))
334     {
335     coq[nb_coq]=noeud1->get_lien_tetra()->get(i);
336     nb_coq++;
337     }
338     if ((nb_coq<4) || (nb_coq>10))
339     {
340     coque.taille=0;
341     crit=0.;
342     return;
343     }
344     //recherche des noeuds qui compose la coquille
345     MG_NOEUD* tet_noeud[50];
346     int nb_tet_noeud=0;
347     coque.taille=nb_coq;
348     for (int i=0;i<nb_coq;i++)
349     {
350     coque.tet[i]=coq[i];
351     if ((coq[i]->get_noeud1()!=noeud1) && (coq[i]->get_noeud1()!=noeud2)) tet_noeud[nb_tet_noeud++]=coq[i]->get_noeud1();
352     if ((coq[i]->get_noeud2()!=noeud1) && (coq[i]->get_noeud2()!=noeud2)) tet_noeud[nb_tet_noeud++]=coq[i]->get_noeud2();
353     if ((coq[i]->get_noeud3()!=noeud1) && (coq[i]->get_noeud3()!=noeud2)) tet_noeud[nb_tet_noeud++]=coq[i]->get_noeud3();
354     if ((coq[i]->get_noeud4()!=noeud1) && (coq[i]->get_noeud4()!=noeud2)) tet_noeud[nb_tet_noeud++]=coq[i]->get_noeud4();
355     coque.volume=coque.volume+fabs(get_volume(coq[i]));
356     }
357     //ordonnancement des noeuds qui compose la coquille en polygone
358     MG_NOEUD* polygone[20];
359     int nb_poly=2;
360     polygone[0]=tet_noeud[0]; tet_noeud[0]=NULL;
361     polygone[1]=tet_noeud[1];tet_noeud[1]=NULL;
362     while (nb_poly!=nb_coq+1)
363     {
364     for (int j=0;j<nb_coq;j++)
365     {
366     if (tet_noeud[2*j]==polygone[nb_poly-1])
367     {
368     polygone[nb_poly++]=tet_noeud[2*j+1];
369     tet_noeud[2*j]=NULL;
370     tet_noeud[2*j+1]=NULL;
371     }
372     if (tet_noeud[2*j+1]==polygone[nb_poly-1])
373     {
374     polygone[nb_poly++]=tet_noeud[2*j];
375     tet_noeud[2*j]=NULL;
376     tet_noeud[2*j+1]=NULL;
377     }
378    
379     }
380     }
381     //etude du positionnement du polygone
382     double vol1=get_volume(polygone[0]->get_coord(),polygone[1]->get_coord(),polygone[nb_poly-1]->get_coord(),noeud1->get_coord());
383     double vol2=get_volume(polygone[0]->get_coord(),polygone[nb_poly-1]->get_coord(),polygone[1]->get_coord(),noeud2->get_coord());
384     if (vol1*vol2<0.)
385     {
386     crit=0.;
387     coque.taille=0;
388     return;
389     }
390     MG_NOEUD *noeuda=noeud1;
391     MG_NOEUD *noeudb=noeud2;
392     if (vol1<0.)
393     {
394     noeuda=noeud2;
395     noeudb=noeud1;
396     }
397     //Examen des solutions
398     int nb_solution,nb_triangle;
399     if (nb_coq==4) {nb_solution=2;nb_triangle=4;}
400     if (nb_coq==5) {nb_solution=5;nb_triangle=10;}
401     if (nb_coq==6) {nb_solution=14;nb_triangle=20;}
402     if (nb_coq==7) {nb_solution=42;nb_triangle=35;}
403     if (nb_coq==8) {nb_solution=132;nb_triangle=56;}
404     if (nb_coq==9) {nb_solution=429;nb_triangle=84;}
405     if (nb_coq==10) {nb_solution=1430;nb_triangle=120;}
406     double crit_triangle[120];
407     for (int i=0;i<nb_triangle;i++)
408     {
409     double crit1=OPERATEUR::qualite_tetra(polygone[tab_face[nb_coq-4][i][0]]->get_coord(),polygone[tab_face[nb_coq-4][i][1]]->get_coord(),polygone[tab_face[nb_coq-4][i][2]]->get_coord(),noeuda->get_coord());
410     double crit2=OPERATEUR::qualite_tetra(polygone[tab_face[nb_coq-4][i][0]]->get_coord(),polygone[tab_face[nb_coq-4][i][2]]->get_coord(),polygone[tab_face[nb_coq-4][i][1]]->get_coord(),noeudb->get_coord());
411     crit_triangle[i]=std::min(crit1,crit2);
412     }
413     /* examen de chaque solution */
414     double crit_opt=0.;
415     int numero_solution= -1;
416     double crit_solution[1430];
417     for (int i=0;i<nb_solution;i++)
418     {
419     double volume=0.;
420     for (int j=0;j<nb_coq-2;j++)
421     {
422     MG_NOEUD* noa=polygone[tab_face[nb_coq-4][tab_solution[nb_coq-4][i][j]][0]];
423     MG_NOEUD* nob=polygone[tab_face[nb_coq-4][tab_solution[nb_coq-4][i][j]][1]];
424     MG_NOEUD* noc=polygone[tab_face[nb_coq-4][tab_solution[nb_coq-4][i][j]][2]];
425     MG_NOEUD* nod=noeuda;
426     volume=volume+fabs(get_volume(noa->get_coord(),nob->get_coord(),noc->get_coord(),nod->get_coord()));
427     nod=noeudb;
428     volume=volume+fabs(get_volume(noa->get_coord(),nob->get_coord(),noc->get_coord(),nod->get_coord()));
429     }
430     double eps=0.0018*pow(volume,0.666666666);
431     if (OPERATEUR::egal(volume,coque.volume,eps))
432     {
433     crit_solution[i]=1.;
434     for (int j=0;j<nb_coq-2;j++)
435     crit_solution[i]=std::min(crit_solution[i],crit_triangle[tab_solution[nb_coq-4][i][j]]);
436     }
437     else crit_solution[i]=0.;
438     if (crit_opt<crit_solution[i])
439     {
440     crit_opt=crit_solution[i];
441     numero_solution=i;
442     }
443    
444     }
445     if (numero_solution==(-1))
446     {
447     crit=0.;
448     coque.taille=0;
449     return;
450     }
451     crit=crit_opt;
452     for (int j=0;j<nb_coq-2;j++)
453     {
454     MG_NOEUD* noa=polygone[tab_face[nb_coq-4][tab_solution[nb_coq-4][numero_solution][j]][0]];
455     MG_NOEUD* nob=polygone[tab_face[nb_coq-4][tab_solution[nb_coq-4][numero_solution][j]][1]];
456     MG_NOEUD* noc=polygone[tab_face[nb_coq-4][tab_solution[nb_coq-4][numero_solution][j]][2]];
457     coque.new_tetra[8*j]=noa;
458     coque.new_tetra[8*j+1]=nob;
459     coque.new_tetra[8*j+2]=noc;
460     coque.new_tetra[8*j+3]=noeuda;
461     coque.new_tetra[8*j+4]=noa;
462     coque.new_tetra[8*j+5]=noc;
463     coque.new_tetra[8*j+6]=nob;
464     coque.new_tetra[8*j+7]=noeudb;
465     }
466     }
467    
468    
469     void MAILLEUR3D::ajouter_ordre_tetra(M3D_TETRA* tet,int num)
470     {
471     std::pair<double,M3D_TETRA*> tmp(pow(tet->get_qualite(),1-2*num),tet);
472     ORDRE_TETRA::iterator p=lst_tetra[num].insert(tmp);
473     lst_tetraid[num][tet->get_id()]=p;
474     }
475    
476     void MAILLEUR3D::supprimer_ordre_tetra(M3D_TETRA* tet)
477     {
478     int num=0;
479     ORDRE_TETRA::iterator it=lst_tetraid[num][tet->get_id()];
480     if (it._M_node==NULL)
481     {
482     num=1;
483     it=lst_tetraid[num][tet->get_id()];
484     }
485     if (it._M_node==NULL) return;
486     lst_tetra[num].erase(it);
487     lst_tetraid[num].erase(lst_tetraid[num].find(tet->get_id()));
488     }